How to solve in MATLAB 2018b ???

72 次查看（过去 30 天）
thiti prasertjitsun
评论： NADA RIFAI 2020-9-16
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 个评论
Siddhartha Ganguly
Siddhartha Ganguly 2020-6-10
编辑：Siddhartha Ganguly 2020-6-10
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
NADA RIFAI 2020-9-16
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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Stephan
Stephan 2019-1-9
编辑：Stephan 2019-1-9
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

更多回答（1 个）

madhan ravi
madhan ravi 2019-1-9
编辑：madhan ravi 2019-1-9
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.

类别

在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by