How to solve in MATLAB 2018b ???
74 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
thiti prasertjitsun
am 9 Jan. 2019
Kommentiert: NADA RIFAI
am 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 Kommentare
Siddhartha Ganguly
am 10 Jun. 2020
Bearbeitet: Siddhartha Ganguly
am 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
am 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
Akzeptierte Antwort
Stephan
am 9 Jan. 2019
Bearbeitet: Stephan
am 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
Weitere Antworten (1)
madhan ravi
am 9 Jan. 2019
Bearbeitet: madhan ravi
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Equation Solving finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!