# Concerning some error For Loop in my optimal control problem

1 view (last 30 days)
Mojtaba Mohareri on 16 Jan 2020
Commented: Mojtaba Mohareri on 18 Jan 2020
There's error ezplot(sol_a.x1,[0 1]) in my code below for solving a special optimal control problem
% State equations
syms x1 x2 p1 p2 u;
for i=0:0.1:1
Dx1 = -2*x2+u;
Dx2 = 2*x1;
% Cost function inside the integral
syms g;
g = 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
Dx1 = subs(Dx1, 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
conA1 = 'x1(0) = 3-i';
conA2 = 'x2(0) = 3-i';
conA3 = 'x1(1) = 0.5-(0.5)*i';
conA4 = 'x2(1) = 0.5-(0.5)*i';
sol_a = dsolve(eq1,eq2,eq3,eq4,conA1,conA2,conA3,conA4);
% plot solutions
figure(1);
ezplot(sol_a.x1,[0 1]);
axis([0 1 -4 3]);
hold on;
end
When I run the above code for every value i=0,0.1,...,1 separately, it works, but for "For Loop" it won't work.

Walter Roberson on 17 Jan 2020
Inside a quoted string 'i' the symbolic engine will always interpret 'i' as sqrt(-1) and never as the current value of i
You should assume that you are no longer permitted to pass character vectors to dsolve and rewrite everything to symbolic expressions
You are doing some invalid calculations. Dx1 only makes sense if x1 is a function, but if it is a function then you cannot take the derivative of H with respect to it. Taking the derivative of a function with respect to a function is an invalid operation in normal calculus and requires Calculus of Variations instead. Unknown functions cannot be assumed to be independent, especially in a situation like this where the functions are being defined by their relationships to each other.

Show 1 older comment
Walter Roberson on 17 Jan 2020
I do not think you can justify the way you create Dp1 and Dp2, but I implemented it anyhow.
% State equations
syms x1(t) x2(t) p1(t) p2(t) H(t) u;
syms X1 X2
dx1 = diff(x1,t);
dx2 = diff(x2,t);
dp1 = diff(p1,t);
dp2 = diff(p2,t);
Dx1 = -2*x2+u;
Dx2 = 2*x1;
% Cost function inside the integral
syms g;
g = u^2;
% Hamiltonian
H(t) = g + p1*Dx1 + p2*Dx2;
% Costate equations
Dp1 = subs( -diff(subs(H,x1,X1),X1), X1, x1); %questionable!!!
Dp2 = subs( -diff(subs(H,x1,X2),X2), X2, x2); %questionable!!!
% solve for control u
Du = diff(H,u);
sol_u = solve(Du, u);
% Substitute u to state equations
Dx1 = subs(Dx1, u, sol_u);
% convert symbolic objects to strings for using 'dsolve'
eq1 = dx1 == Dx1;
eq2 = dx2 == Dx2;
eq3 = dp1 == Dp1;
eq4 = dp2 == Dp2;
sol_h = dsolve(eq1,eq2,eq3,eq4);
%% use boundary conditions to determine the coefficients
for i=0:0.1:1
conA1 = x1(0) == 3-i;
conA2 = x2(0) == 3-i;
conA3 = x1(1) == 0.5-(0.5)*i;
conA4 = x2(1) == 0.5-(0.5)*i;
sol_a = dsolve(eq1,eq2,eq3,eq4,conA1,conA2,conA3,conA4);
% plot solutions
figure(1);
fplot(sol_a.x1,[0 1]);
axis([0 1 -4 3]);
hold on;
end
hold off
Mojtaba Mohareri on 18 Jan 2020
It works properly. I appreciate you. Thank you so much for your help.
Mojtaba Mohareri on 18 Jan 2020
I used twise "For Loop" with different conditions in my main problem. I want to use different colors for each group of solutions obtained by each loop (red color for loop 1 and blue color for loop 2). For this, I've written ez1=ezplot(sol_a.x1,[0 1]);
set(ez1,'color',[1 0 0])
at the end of loop1 and
ez2=ezplot(sol_a.x1,[0 1]);
set(ez2,'color',[0 0 1])
at the end of loop 2.
Finaly I've used legend('Lower Bounds','Upper Bounds') but my problem is that these two exprssions are red.