How to solve in MATLAB 2018b ???
72 views (last 30 days)
Show older comments
thiti prasertjitsun
on 9 Jan 2019
Commented: NADA RIFAI
on 16 Sep 2020
Hello, I was trying to solve a system equation.
clear all
% State equations
syms x1 x2 p1 p2 u;
Dx1 = x2;
Dx2 = -x2 + u;
% Cost function inside the integral
syms g;
g = 0.5*u^2;
% Hamiltonian
syms p1 p2 H;
H = g + p1*Dx1 + p2*Dx2;
% Costate equations
Dp1 = -diff(H,x1);
Dp2 = -diff(H,x2);
% solve for control u
du = diff(H,u);
sol_u = solve(du,u);
% Substitute u to state equations
Dx2 = subs(Dx2,u,sol_u);
% convert symbolic objects to strings for using 'dsolve'
eq1 = strcat('Dx1=',char(Dx1));
eq2 = strcat('Dx2=',char(Dx2));
eq3 = strcat('Dp1=',char(Dp1));
eq4 = strcat('Dp2=',char(Dp2));
sol_h = dsolve(eq1,eq2,eq3,eq4);
% use boundary conditions to determine the coefficients
% case a: (a) x1(0)=x2(0)=0; x1(2) = 5; x2(2) = 2;
conA1 = 'x1(0) = 0';
conA2 = 'x2(0) = 0';
conA3 = 'x1(2) = 5';
conA4 = 'x2(2) = 2';
sol_a = dsolve(eq1,eq2,eq3,eq4,conA1,conA2,conA3,conA4);
% plot both solutions
figure(1);
ezplot(sol_a.x1,[0 2]); hold on;
ezplot(sol_a.x2,[0 2]);
ezplot(-sol_a.p2,[0 2]); % plot the control: u=-p2
axis([0 2 -1.6 7]);
text(0.6,0.5,'x_1(t)');
text(0.4,2.5,'x_2(t)');
text(1.6,0.5,'u(t)');
xlabel('time');
ylabel('states');
title('Solutions comparison (case a)');
% case b: (a) x1(0)=x2(0)=0; p1(2) = x1(2) - 5; p2(2) = x2(2) -2;
eq1b = char(subs(sol_h.x1,'t',0));
eq2b = char(subs(sol_h.x2,'t',0));
eq3b = strcat(char(subs(sol_h.p1,'t',2)),...
'=',char(subs(sol_h.x1,'t',2)),'-5');
eq4b = strcat(char(subs(sol_h.p2,'t',2)),...
'=',char(subs(sol_h.x2,'t',2)),'-2');
sol_b = solve(eq1b,eq2b,eq3b,eq4b);
C2 = double(sol_b.C2);
C3 = double(sol_b.C3);
C4 = double(sol_b.C4);
C5 = double(sol_b.C5);
sol_b2 = struct('x1',{subs(sol_h.x1)},'x2',{subs(sol_h.x2)}, ...
'p1',{subs(sol_h.p1)},'p2',{subs(sol_h.p2)});
figure(2);
ezplot(sol_b2.x1,[0 2]); hold on;
ezplot(sol_b2.x2,[0 2]);
ezplot(-sol_b2.p2,[0 2]); % plot the control: u=-p2
axis([0 2 -0.5 3]);
text(0.9,0.5,'x_1(t)');
text(0.4,1,'x_2(t)');
text(0.2,2.5,'u(t)');
xlabel('time');
ylabel('states');
title('Solutions comparison (case b)');
On Matlab 2015a, I can get the final results.
But on Matlab 2018b, only error returns sol_b = solve(eq1b,eq2b,eq3b,eq4b);
So are there some updates or changes between 2015a and 2018b?
And how can I solve algebraic equations correctly in Matlab 2018b?
Thanks!
2 Comments
Siddhartha Ganguly
on 10 Jun 2020
Edited: Siddhartha Ganguly
on 10 Jun 2020
This problem is for fixed final state and final time, do you have any idea how to change this if i have final time t_f and state x_f both free?
NADA RIFAI
on 16 Sep 2020
Hello,
I have the same type of problem but with constraints, how can I include them in the solution?
Can you please help me solve this problem.
thank you
Accepted Answer
Stephan
on 9 Jan 2019
Edited: Stephan
on 9 Jan 2019
Hi,
the following runs for me in 2018b:
clear all
% State equations
syms x1 x2 p1 p2 u;
Dx1 = x2;
Dx2 = -x2 + u;
% Cost function inside the integral
syms g;
g = 0.5*u^2;
% Hamiltonian
syms p1 p2 H;
H = g + p1*Dx1 + p2*Dx2;
% Costate equations
Dp1 = -diff(H,x1);
Dp2 = -diff(H,x2);
% solve for control u
du = diff(H,u);
sol_u = solve(du,u);
% Substitute u to state equations
Dx2 = subs(Dx2,u,sol_u);
% convert symbolic objects to strings for using 'dsolve'
eq1 = strcat('Dx1=',char(Dx1));
eq2 = strcat('Dx2=',char(Dx2));
eq3 = strcat('Dp1=',char(Dp1));
eq4 = strcat('Dp2=',char(Dp2));
sol_h = dsolve(eq1,eq2,eq3,eq4);
% use boundary conditions to determine the coefficients
% case a: (a) x1(0)=x2(0)=0; x1(2) = 5; x2(2) = 2;
conA1 = 'x1(0) = 0';
conA2 = 'x2(0) = 0';
conA3 = 'x1(2) = 5';
conA4 = 'x2(2) = 2';
sol_a = dsolve(eq1,eq2,eq3,eq4,conA1,conA2,conA3,conA4);
% plot both solutions
figure(1);
ezplot(sol_a.x1,[0 2]); hold on;
ezplot(sol_a.x2,[0 2]);
ezplot(-sol_a.p2,[0 2]); % plot the control: u=-p2
axis([0 2 -1.6 7]);
text(0.6,0.5,'x_1(t)');
text(0.4,2.5,'x_2(t)');
text(1.6,0.5,'u(t)');
xlabel('time');
ylabel('states');
title('Solutions comparison (case a)');
% case b: (a) x1(0)=x2(0)=0; p1(2) = x1(2) - 5; p2(2) = x2(2) -2;
eq1b = subs(sol_h.x1,'t',0);
eq2b = subs(sol_h.x2,'t',0);
eq3b = subs(sol_h.p1,'t',2) == subs(sol_h.x1,'t',2)-5;
eq4b = subs(sol_h.p2,'t',2) == subs(sol_h.x2,'t',2)-2;
sol_b = solve(eq1b,eq2b,eq3b,eq4b);
C1 = double(sol_b.C1);
C2 = double(sol_b.C2);
C3 = double(sol_b.C3);
C4 = double(sol_b.C4);
sol_b2 = struct('x1',{subs(sol_h.x1)},'x2',{subs(sol_h.x2)}, ...
'p1',{subs(sol_h.p1)},'p2',{subs(sol_h.p2)});
figure(2);
ezplot(sol_b2.x1,[0 2]); hold on;
ezplot(sol_b2.x2,[0 2]);
ezplot(-sol_b2.p2,[0 2]); % plot the control: u=-p2
axis([0 2 -0.5 3]);
text(0.9,0.5,'x_1(t)');
text(0.4,1,'x_2(t)');
text(0.2,2.5,'u(t)');
xlabel('time');
ylabel('states');
title('Solutions comparison (case b)');
Best regards
Stephan
More Answers (1)
madhan ravi
on 9 Jan 2019
Edited: madhan ravi
on 9 Jan 2019
Remove the ' ' single quote in the equation and change your second equal to sign as == (2018b doesn't support string for equations lookup https://www.mathworks.com/help/symbolic/solve.html clearly states the proper usage).
Also see https://www.mathworks.com/help/symbolic/solve-a-system-of-algebraic-equations.html for more examples.
0 Comments
See Also
Categories
Find more on Equation Solving in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!