1 view (last 30 days)

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.

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

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.