Solve Function Doesn't Give The Solutions
26 views (last 30 days)
Show older comments
Ivan Dwi Putra
on 15 Sep 2020
Commented: Ivan Dwi Putra
on 25 Sep 2020
This is the parameter
clear; close; clc;
syms a1_head a2_head b hstar
%Parameter Massa
m1 = 8095; % massa train set 1 dalam kg
m2 = 8500; % massa train set 2 dalam kg
g = 10;
c_0_1 = 0.01176;
c_1_1 = 0.00077616;
c_2_1 = 4.48 ;
c_0_2 = 0.01176 ;
c_1_2 = 0.00077616;
c_2_2 = 4.48;
v_0 = 300;
hstar = 120;
a_1 = -1./m1.*(c_1_1 + 2.*c_2_1.*v_0);
a_2 = -1./m2.*(c_1_2 + 2.*c_2_2.*v_0);
a_1_head = 1-(a_1.*hstar);
a_2_head = 1-(a_2.*hstar);
b = 1;
This is the model
% Model data
A = sym(zeros(4,4));
A(1,2) = a_1_head;
A(3,2) = (a_2_head) - 1; A(3,4) = a_2_head;
display(A);
B = sym(zeros(4,2));
B(1,1) = -b*hstar;
B(2,1) = b;
B(3,2) = -b*hstar ;
B(4,1) = -b; B(4,2) = b;
display(B);
This is the Q and R matrix for LQR
% Q and R matrices for ARE
Q = [10 0 0 0;
0 10 0 0;
0 0 10 0;
0 0 0 10;
];
R = sym(zeros(2,2)); R(1,:) = [1 0]; R(2,:) = [0 1]; display(R);
This is the S matrix that I want to find the value
% Matrix S to find
svar = sym('s',[1 16]);
S = [svar(1:4); svar(5:8); svar(9:12); svar(13:16)];
S(2,1) = svar(2);
S(3,1) = svar(3);
S(3,2) = svar(7);
S(4,1) = svar(4);
S(4,2) = svar(8);
S(4,3) = svar(12);
display(S);
This is the Algebra Ricatti Equation (ARE) in LQR Control
% LHS of ARE: A'*S + S*A' - S*B*Rinv*B'*S
left_ARE = transpose(A)*S + S*A - S*B*inv(R)*transpose(B)*S;
display(left_ARE);
% RHS of ARE: -Q
right_ARE = -Q;
display(right_ARE);
I use the solve function to find value of S matriks
Find S
syms s1 s2 s3 s4 s6 s7 s8 s11 s12 s16
[Sol_s] = solve(left_ARE == right_ARE,s1,s2,s3,s4,s6,s7,s8,s11,s12,s16,'Real',true)
But the answer(solution) from solve like that. There is a variable x that i don't know
Sol_s =
struct with fields:
s1: [2×1 sym]
s2: [2×1 sym]
s3: [2×1 sym]
s4: [2×1 sym]
s6: [2×1 sym]
s7: [2×1 sym]
s8: [2×1 sym]
s11: [2×1 sym]
s12: [2×1 sym]
s16: [2×1 sym]
>> vpa(Sol_s.s1)
ans =
0.00000000000000000040859201997183266363736715737362*x^7 - 0.000000000001292747173772775983242070865976*x^6 + 0.0000015262175064767887069297714018443*x^5 - 0.79679406124260841586996673723117*x^4 + 155208.81351720644561857433778132*x^3 + 7112858.7612670786982068414016205*x^2 + 108634010.42148359450942426361058*x + 552988510.84561362711489334430607
0.00000000000000000020499055506656593566660625665166*x^7 - 0.00000000000064871158271130010923295137360766*x^6 + 0.00000076601975870779972515899504293451*x^5 - 0.39998727113373076657212134425168*x^4 + 77926.337795483306183297402674704*x^3 + 3571036.0433220579501871782146977*x^2 + 54537896.71647223250810434173391*x + 277607756.88722181441803051789872
>> vpa(Sol_s.s2)
ans =
- 0.00000000000000000031260583959047249814694055447793*x^7 + 0.00000000000098905577988552374034850303549741*x^6 - 0.0000011676793205823800112784869250653*x^5 + 0.60961154908164182468172656314696*x^4 - 118747.17710114932017447930936862*x^3 - 5457407.0206974194274463617331546*x^2 - 83573329.145320152948094088924873*x - 426479043.54369671880276326751334
- 0.00000000000000000032365118353566263532308345273194*x^7 + 0.0000000000010242240764674830710076551817231*x^6 - 0.0000012094370452875433539209714880716*x^5 + 0.63152332992457965547656341647457*x^4 - 123034.60895321337878520443201568*x^3 - 5655911.3774460846397440247453401*x^2 - 86636236.777808615440767863132122*x - 442230554.80769573806658012394346
>> vpa(Sol_s.s3)
ans =
0.0000000000000000000010776575438999520316538127429084*x^7 - 0.0000000000000034096078345527396375991806997797*x^6 + 0.0000000040253825881136818413774049305516*x^5 - 0.0021015349611454802346225087217245*x^4 + 409.36084570596063399343217949922*x^3 + 18938.36880217893351979479924197*x^2 + 291993.50120838165763852700250626*x + 1500485.4678457105543481903948988
- 0.00000000000000000041666832842512374615974790886976*x^7 + 0.0000000000013185854853866499151812140304409*x^6 - 0.0000015570286714848943410537775793585*x^5 + 0.8130229545475701019194852997876*x^4 - 158394.79247503034137610097150089*x^3 - 7259023.3274984798704782100693356*x^2 - 110868678.30972936220866037440766*x - 564374966.98865965783991745150745
>> vpa(Sol_s.s4)
ans =
0.0000000000000000000014228409703424497861474646268027*x^7 - 0.0000000000000045017366112576086307555233406427*x^6 + 0.000000005314750975152131845714002972785*x^5 - 0.0027746779808566367452971500674644*x^4 + 540.48408167395620003882852791153*x^3 + 24757.747980942663555794256513835*x^2 + 377958.31952099355656247974473802*x + 1923166.530892430026770775363493
0.0000000000000000000023193355511678347273166656795793*x^7 - 0.0000000000000073397520256740640487837787182755*x^6 + 0.0000000086670182256809939350882653924303*x^5 - 0.0045255973789649557650852704388279*x^4 + 881.68616552361166280706272097366*x^3 + 40393.656783854177046779159405732*x^2 + 616754.69331512724884616078670876*x + 3138681.4020773718736419352352328
>> vpa(Sol_s.s6)
ans =
- 0.000000000000000022006336529571817230209213311748*x^7 + 0.000000000069626005207446692167524264364454*x^6 - 0.000082200475821668055631814879425621*x^5 + 42.914497164630049410319474600551*x^4 - 8359387.0930143709758499635784225*x^3 - 382387662.93012000327080222222905*x^2 - 5827384455.2847427100003526481025*x - 29587579250.550523363703915448321
- 0.000000000000000025002231287155086974127450112697*x^7 + 0.000000000079121875258332337913776610719734*x^6 - 0.000093429687431695914874771958800029*x^5 + 48.785538562659557207332817154875*x^4 - 9504498.929505799236089462678773*x^3 - 435370306.00846953347659747350646*x^2 - 6644367672.3679146002218490792474*x - 33786329437.38166020805340239827
>> vpa(Sol_s.s7)
ans =
0.00000000000000000025569963231254537081649139584652*x^7 - 0.00000000000080900975738079800385855474939705*x^6 + 0.00000095511681683491535731968205794489*x^5 - 0.49863865781308244522352509310329*x^4 + 97130.498837171086515556658964247*x^3 + 4494101.0379877774836523120182835*x^2 + 69298502.116075437676103381862363*x + 356149465.9042365809352931154343
0.0000000000000000000015023043416658044818117509971356*x^7 - 0.0000000000000047541811528920007743725411473973*x^6 + 0.0000000056138919601726527683759200771913*x^5 - 0.002931367433451308848204776259538*x^4 + 571.09503266523001592434936146943*x^3 + 26159.338852936025767577334375426*x^2 + 399345.5827675731856692656648735*x + 2031937.4929783588149908825597278
>> vpa(Sol_s.s8)
ans =
- 0.000000000000000000017423118618084831929830862790716*x^7 + 0.000000000000055125128523282993141578608099266*x^6 - 0.000000065080736244752471556193477098825*x^5 + 0.033976773274815242226234938519543*x^4 - 6618.3918152997107197929705528418*x^3 - 303061.20015840970654893674015787*x^2 - 4623463.04197949737532115401195*x - 23501224.452355499378972417099829
- 0.000000000000000000018476917340887162236508871167814*x^7 + 0.000000000000058471915042480914496818710438269*x^6 - 0.000000069045541385168727596306044837525*x^5 + 0.036053035957476786546322342447713*x^4 - 7023.9263332222742397498791136158*x^3 - 321824.06004694090757346294728009*x^2 - 4912773.5723496489797387292827707*x - 24988039.623927850151456500151167
>> vpa(Sol_s.s11)
ans =
0.00000000000000000041119629783273028264519058921762*x^7 - 0.0000000000013009868665244108119791105781729*x^6 + 0.0000015359452882952913401827632335578*x^5 - 0.80187265434412110159568710275759*x^4 + 156198.08052664176034300508696723*x^3 + 7158323.8866181131831523373220071*x^2 + 109330254.65157047767529156319121*x + 556541496.93083652400446027963023
- 0.00000000000000000020771378465236870650240546526077*x^7 + 0.00000000000065732949437268293588570845859109*x^6 - 0.00000077619606862262249756751814375767*x^5 + 0.40530096411080625205677575409412*x^4 - 78961.561847579934284088187742172*x^3 - 3618623.8712567317582544282624458*x^2 - 55266816.721319493314377599392582*x - 281328390.05487908210968779383083
>> vpa(Sol_s.s12)
ans =
- 0.0000000000000000000015153372531006669942761978246883*x^7 + 0.0000000000000047943863006019005344930952094381*x^6 - 0.0000000056602533372161527651811224147036*x^5 + 0.0029550548809761743405798500431343*x^4 - 575.62005058550944282823558765466*x^3 - 26349.726554223897542645543275317*x^2 - 401979.85441516660058061899440103*x - 2043879.8478415220815018039779193
- 0.00000000000000000000081703120950203024550491468244366*x^7 + 0.0000000000000025855708727820632744112375708782*x^6 - 0.000000003053126265508341166712345315239*x^5 + 0.00159422994551364691688049417929*x^4 - 310.59113285838164688271335950423*x^3 - 14234.317930918151279201825030306*x^2 - 217409.11054755406317689512183527*x - 1106743.9090990130586510526755051
>> vpa(Sol_s.s16)
ans =
x
x
2 Comments
Accepted Answer
Walter Roberson
on 21 Sep 2020
Warning: Solutions are parameterized by the symbols: x. To include parameters and conditions in the solution, specify the
There are several possibilities.
- When you use non-polynomial equations, it is not uncommon for there to be situations under which something that is being divided by could hypothetically equal 0. solve() tries to figure out all combinations of values under which expressions it divides by can be 0, and lists out all of the combinations as separate solutions, together with the solution that would be implied when that combination is 0. It is not uncommon in such a situation for a variable to vanish under the hypothesis that some of the inputs are 0, and when that happens the number of variables in the equations can drop below the number of equations, giving you an indefinite number of solutions. MATLAB generates parametric solutions for those situations, generating a variable to take the parametric solution with regards to.
- When you use nonlinear equations, there can potentially be families of solutions instead of single solutions. This is most obvious when you are using periodic or semi-periodic functions such as trig functions or bessel functions, but it can also occur even with something like x*exp(x) == constant, which can have up to two real solutions for some constants. When MATLAB detects that multiple solutions are possible. it will generate a variable that stands in for the multiple solutions, and then in the "conditions" returned when returnconditions is true, it will describe the periodicity or constraints
- When you use inequalities, there are multiple solutions. MATLAB synthesizes a variable with that is as refined as practical, and then in the "conditions" it describes the constraints using inequalities
- Sometimes there are exact solutions based upon the coefficients, but the exact solutions would be expressed as the roots of an expression involving variables. In such a situation, MATLAB might generate the exact solutions in-line, without using an extra parameter, but in other cases the expression would be complex enough that instead MATLAB generates an extra variable and then in "conditions" it describes the polynomial whose roots would form the variable. In this case, MATLAB can essentially be "lazy", not bothering to work out something that has a perfectly good solution.
- Sometimes equations can be solved in terms of systems of polynomials, in which one variable becomes the primary variable, and the other variables can then be expressed in terms of polynomial expressions of the primary variable. This gives you a family of solutions for the variables, in which each variable other than the first has a exact expression based upon the value of the first variable, but the first variable has multiple roots. Instead of "expanding" the solutions to all of those other variables, MATLAB sometimes generates an intermediate variable and phrases the other variables in terms of that variable, and in the conditions lists the generated variable in terms of the polynomial whose roots it is. This can be much more compact, and can make the solutions for the other variables understandable when expanding the roots out would just tend to obscure the solutions.
So, MATLAB has found solutions to the equations, but those solutions either have conditions attached (e.g., avoiding division by 0) or have multiple solutions (periodic / semi-periodic), or the solutions are most naturally expressed as a parameterized family.
In such cases, it is often necessary to go through and carefully solve for the combinations of circumstances that make the listed conditions true.
In this particular case, it is the last of those that applies: the solutions it can find are families of solutions based upon all of the real roots of a degree 8 polynomial. That particular polynomial happens to have all 8 of its roots real, so the two solutions (per variable) that MATLAB finds, parameterized in terms of the roots of a polynomial of degree 8, means that there are actually a total of 16 different solutions. (The second family of solutions uses a different polynomial of degree 8, by the way -- the two are very similar but not quite the same.)
X = solve(sol_s1.conditions(1));
subs(sol_s1.s1(1), sol_s1.parameters, X)
19 Comments
Walter Roberson
on 25 Sep 2020
I was worried that due to round-off errors that the values might not match, but when I take double() of both versions, then there are exact matches (which is rare for floating point.) So, yes, both codes will produce the same double() output. The solve() version gives you information about why those particular outputs, and does not give you any warning messages; the vpasolve() version is much less bother but gives you warning messages.
Because of the warning messages that vpasolve() gives, it would be natural to wonder if the vpasolve() results are correct. The fact that double() of the results match the solutions from solve() tells you that double() of the results is as close you can get in double precision to the actual solutions.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!