Outputting all possible solutions to system of equations with constraints.

I'm trying to solve a system of nonlinear equations for which I have more unknowns than equations. I have conditions for a few of the variables to limit the number of solutions. And I would also like to output the solutions into an easily visually parseable table of numeric values where all values are aligned to their respective solutions.

What I have to begin with is the code below. This is giving me a parameterized solution warning, though. Assuming I haven't made a syntax mistake this should produce a solution. I have tested a solution by explicitly declaring variables within the bounds I've set to confirm this. For example, if:

V_as = 275

f_s = 140

then without conditions MATLAB would return:

C_ms = 836e-6

BL = 3.41

Q_ms = 3.8

M_md = 1.48e-3

Maybe just setting these boundaries still isn't constraining it enough? Should I try specifically declaring an array of values for f_s and V_as? Would that help? And how would I go about doing that?

Also, I would like to output all these variables into a table of numeric values aligned to each solution but I'm not sure how to do that since they're stored in a symbolic cell array. Something like:

 output_table = [double(sol.BL), double(sol.M_md), double(Q_ms), etc...]

And how to make sure that each variable in a particular row is from the same solution as the variable next to it in the same row? And how to set the columns to have headers representing each variable?

syms f_s V_as M_ms Q_es Q_ms C_ms M_md BL;

D = 4.4e-2;
S_d = pi*(D/2)^2;
SPL = 85;
Q_ts = 0.35;
R_e = 3.3;

assume(in(f_s/10,'integer') & f_s>=100 & f_s<=150);
assume(in(V_as/25,'integer') & V_as>=100 & V_as<=1000);
assume(M_md>=1e-3 & M_md<5e-3);
assume(BL>3 & BL<3.5);

assume(C_ms>0 & C_ms<1500e-6);
assume(Q_ms>0 & Q_ms<10);

e1 = f_s == 1/(2*pi*sqrt(M_ms*C_ms));
e2 = Q_es == R_e/BL^2*sqrt(M_ms/C_ms);
e3 = Q_ts == Q_ms*Q_es/(Q_ms+Q_es);
e4 = V_as*10^-6 == 1.21*343^2*S_d^2*C_ms;
e5 = M_ms == S_d^2*(M_md/S_d^2+2*8*1.21/(3*pi^2*(D/2)));
e6 = SPL == 20*log10(1.21*BL*S_d*sqrt(R_e)/(2*pi*R_e*M_ms*2e-5));

sol = solve(e1,e2,e3,e4,e5,e6,f_s,V_as,Q_es,Q_ms,M_ms,M_md,BL,C_ms);

Answer 1

[Walter Roberson]

There do appear to be algebraic solutions, but they are numerically unjustifiable for your inputs. Your data such as D = 4.4e-2 has between 1 and 4 significant figures, so you cannot justify calculations that carry more than at most 4 significant figures. However, to achieve the conditions that in(f_s/10,'integer') and in(V_as/25,'integer') then because of the Pi and the various squares and sqrt(), some of your values such as C_ms and BL must be infinitely precise, such as BL must be exactly (14235529/62500)*sqrt(330^(1/17)*sqrt(10))*330^(8/17)/(f_s_10^2*V_as_25) where f_s_10 and V_as_25 are integers

One of the algebraic solutions appears to be at f_s = 150, V_as = 250; there are other algebraic solutions as well. But they require infinite precision for the integer constraints to be met, and they require assuming that the various low-precision constants such as D = 4.4e-2 are "actually" infinitely precise such as D = 44/1000 and not D = (44/1000 +/- 5/10000) as should be implied by the significant figures.

[Matthew Robinson]

Thanks for your reply. I think I get it. By putting those constraints in there I was hoping to help the solver, but it seems I may have "constrained" it too much?

The reason I am using different levels of precisions for some of these values is because they are representing units of measure. For instance C_ms always *10^-6 because it's units of measure are um/N. And D is simply just cm. V_as and f_s don't have to be integers, but I don't care to see they're values past the decimal point, just multiples of 10s and 25s were enough precision for me for those values. I was worried not constraining them to integers would give me too many solutions. However, if I simply just take away the integer constraint MATLAB just tells me no explicit solution can be found.

What I tried doing before this was setting up a big list in Excel with an array of varying f_s and V_as values that I know typically work, rearrange all these equations manually and filtering according to the constraints listed, which worked fine for these values in this situation. But I was hoping to come up with a more robust solution that can help me in all cases. Would it be "better" to setup f_s and V_as as an array of values instead, and just iterate the solver over each possible combination of those two? Or is that basically what it's doing already? By the way, I tried doing this but I can't figure out how to get it to work, syntax-wise.

I think you can probably see what I really need help with is setting up this problem on a bit of a higher level - then with how to implement it in MATLAB. But I'm afraid that may be outside of the scope of what should be asked of you or this community. If you're willing to help I would be grateful but I understand if you can't.

[Walter Roberson]

You can solve symbolically and then re-arrange for the purposes of iteration over f_s and V_as

f_s = (500/3993)*sqrt(33)*sqrt(sqrt(33)*330^(4/17)*sqrt(sqrt(33)*330^(3/34))/(C_ms*BL))*330^(4/17)/Pi
V_as = (208422380089/6250000)*Pi^2*C_ms
Q_es = (121/100000)*sqrt(330^(5/34)*sqrt(sqrt(33)*330^(3/34))/(C_ms*BL))*330^(9/17)/BL
Q_ms = (-2964500*330^(9/17)*BL^2*C_ms*sqrt(330^(5/34)*sqrt(sqrt(33)*330^(3/34))/(C_ms*BL))-3382071*330^(7/34)*sqrt(sqrt(33)*330^(3/34)))/(-9663060*330^(7/34)*sqrt(sqrt(33)*330^(3/34))[phone]*BL^3*C_ms)
M_ms = (1331/30000000)*BL*330^(7/34)*sqrt(sqrt(33)*330^(3/34))
M_md = (1331/30000000)*BL*330^(7/34)*sqrt(sqrt(33)*330^(3/34))-161051/2343750000
BL = BL
C_ms = C_ms

Those last two indicate that BL and C_ms are "free variables".

For the purposes of iteration, given V_as you solve the V_as to get C_ms . Then using that value of C_ms, and given f_s, you can solve f_s to get BL. The other variables are all in terms of those, so you can proceed from there.

Or if you just change the list of variables you solve for, and eliminate the one of the pair that leads to some negative values, then

BL = (1294139/75)*330^(9/17)*sqrt(33)*sqrt(sqrt(33)*330^(15/34))/(V_as*f_s^2)
C_ms = (6250000/208422380089)*V_as/Pi^2
M_md = -(161051/9375000000)*(4*V_as*f_s^2-485302125)/(V_as*f_s^2)
M_ms = 208422380089/(25000000*V_as*f_s^2)
Q_es = (11/7058940000000)*Pi*V_as*f_s^3*sqrt(10890)
Q_ms = (307461*V_as^2*Pi^2*f_s^6+6341281100000*Pi*V_as*f_s^3*sqrt(10890))/(878460*Pi^2*V_as^2*f_s^6-4069338437094000000000000)

Vectorize, use ngrid to create a grid of V_as and f_s values, calculate the variables, reject any result that does not fit the constraints.

[Matthew Robinson]

I think I like your second suggestion best. But I'm not clear on two things. First, what do you mean by eliminate the one pair that leads to some negative values?

Second, I'm not sure how you obtained those equations. I can't get them. They're the contents of the returned conditions of the solver, correct? Did you remove all the constraints I had?

Also, is there a better way to visualize the conditions matrix? In the command window it's really hard to read in one line.

[Walter Roberson]

I removed all of the constraints. solve() returns two sets of solutions, which are very similar but one of the sets is provably negative for Q_* values, so we can eliminate that set, leaving the set I post second.

I did not examine the conditions, but with the above formulas you could work backwards to constrain f_s and V_as

[Matthew Robinson]

But how did you get the above formulas? I can't seem to get them. I just get M_ms = z1 and C_ms = x and so on. I could copy your formulas but I'm trying to get better at MATLAB so I would like to know how you did that.

[Walter Roberson]

I used Maple and a straight forward implementation,just making the plain syntax changes like = becomes := . The only thing that I did that is not immediately obvious was to convert the floating point to rational, but that is done automatically by the MATLAB symbolic engine by default so that is the same.

[Matthew Robinson]

Ok, thanks for the clarification. Seems there's a reason why no one has come up with a better solution for this in my industry. It's not so cut and dry.