MATLAB: HOW I plotted the x1 and x2 and u in optimal control poblem

optimal control - structure - plot

% IN THE NAME OF GOD 
% HAMED GHAFFARI 
% page 198 - 202 
% Book : Kirk
% exp 5.1-1
clc
clear 
% state equation 
syms x1 x2 u p1 p2 ;
DX1 = x2 ;
DX2 = -x2 + u ;
% cost function inside the integral 
syms j 
j = 1/2 * u^2 ;
% Hamiltonian 
syms p1 p2 H
H = j + p1*DX1 +p2*DX2 ;
% costate equations 
DP1 = -diff(H,x1) 
DP2 = -diff(H,x2)
% solve for u control 
du = diff(H,u)
sol_u = solve(du,u)
% substitute u to state eq DX2
DX2 = subs(DX2 , u , sol_u)
% Cocatenate String horizontal (الحاق کردن عبارات به صورت رشته ای)
% convert symbolic adjective to string for using 'dsolve'
eq1 = strcat('DX1 = ' , char(DX1))
eq2 = strcat('DX2 = ' , char(DX2))
eq3 = strcat('DP1 = ' , char(DP1))
eq4 = strcat('DP2 = ' , char(DP2))
% solve Hamiltonian 
sol_h = dsolve(eq1,eq2,eq3,eq4)
sol_Dx1 = sol_h.X2
sol_Dx2 = sol_h.X1
sol_DP1 = sol_h.P1
sol_DP2 = sol_h.P2
% case Boundry condition x1(0) = x2(0) = 0 ; x1(2) = 5 , x2(2) = 2
cond1 = 'X1(0)= 0' ;
cond2 = 'X2(0)= 0' ;
cond3 = 'X1(2)= 5' ;
cond4 = 'X2(2)= 2' ;
sol_Hnew = dsolve(eq1,eq2,eq3,eq4,cond1,cond2,cond3,cond4)
sol_eq1 = sol_Hnew.X1
sol_eq2 = sol_Hnew.X2
sol_eq3 = sol_Hnew.P1
sol_eq4 = sol_Hnew.P2

Best Answer

% IN THE NAME OF GOD 
% HAMED GHAFFARI 
% page 198 - 202 
% Book : Kirk
% exp 5.1-1
%clc
%clear 
% state equation 
syms x1(t) x2(t) u(t) p1(t) p2(t)
DX1 = x2 ;
DX2 = -x2 + u ;
% cost function inside the integral 
syms j 
j = 1/2 * u^2 ;
% Hamiltonian 
H = j + p1*DX1 +p2*DX2 ;
% costate equations 
DP1 = -functionalDerivative(H,x1) 
DP2 = -functionalDerivative(H,x2)
% solve for u control 
du = functionalDerivative(H,u)
syms ut
sol_u(t) = solve( subs(du, u, ut), ut)
% substitute u to state eq DX2
DX2 = subs(DX2 , u , sol_u)
% Cocatenate String horizontal (الحاق کردن عبارات به صورت رشته ای)
% convert symbolic adjective to string for using 'dsolve'
dx1 = diff(x1);
dx2 = diff(x2);
dp1 = diff(p1);
dp2 = diff(p2);
eq1 = dx1 == DX1;
eq2 = dx2 == DX2;
eq3 = dp1 == DP1;
eq4 = dp2 == DP2;
% solve Hamiltonian 
sol_h = dsolve(eq1,eq2,eq3,eq4)
sol_Dx1 = simplify(sol_h.x1)
sol_Dx2 = simplify(sol_h.x2)
sol_DP1 = simplify(sol_h.p1)
sol_DP2 = simplify(sol_h.p2)
% case Boundry condition x1(0) = x2(0) = 0 ; x1(2) = 5 , x2(2) = 2
cond1 = dx1(0) == 0 ;
cond2 = dx2(0) == 0 ;
cond3 = x1(2) == 5;
cond4 = x2(2) == 2;
sol_Hnew = dsolve(eq1,eq2,eq3,eq4,cond1,cond2,cond3,cond4)
sol_eq1 = simplify(sol_Hnew.x1)
sol_eq2 = simplify(sol_Hnew.x2)
sol_eq3 = simplify(sol_Hnew.p1)
sol_eq4 = simplify(sol_Hnew.p2)
Related Question