问题 此问题出自Solving optimal control problems with MATLAB.pdf

此问题matlab实现程序

function er3OC_sym

%EG3OC_sym Example 3 of optimal control tutorial.

% This example is from D.S.Naidu, "Optimal contorl systems"

% page 77-80, Example 2.14

% Symbolic toolbox is used to get the analytical solution

sol = dsolve('Dx1 = x2, Dx2 = -p2, Dp1 = 0, Dp2 = -p1',...

'x1(0) = 1, x2(0) = 2, x1(tf) = 3, p2(tf) = 0');

eq1 = subs(sol.x1) - 'x1tf';

eq2 = subs(sol.x2) - 'x2tf';

eq3 = subs(sol.p1) - 'p1tf';

eq4 = subs(sol.p2) - 'p2tf';

eq5 = sym('p1tf*x2tf - 0.5*p2tf^2');

%%

sol_2 = solve(eq1, eq2, eq3, eq4, eq5);

tf = sol_2.tf;

x1tf = sol_2.x1tf;

x2tf = sol_2.x2tf;

x1 = subs(sol.x1);

x2 = subs(sol.x2);

p1 = subs(sol.p1);

p2 = subs(sol.p2);

%%

sol_book = {@(t)((4/54).*t.^3-(2/3)*t.^2+2.*t+1),...

@(t)((4/18).*t.^2-(4/3).*t+2)};

u = @(t)((4/9).*t-(4/3));

t = double(tf);

time = linspace(0,t,20);

s_book = [sol_book{1}(time);sol_book{2}(time)];

ut = u(time);

figure(1);

ezplot(x1,[0 t]); hold on;

ezplot(x2,[0 t]);

plot(time, s_book,'*');

plot(time, ut, 'k:');

axis([0 t -1.5 3]);

text(1.3,2.5,'x_1(t)');

text(1.3,1,'x_2(t)');

text(1.3,-.5,'u(t)');

xlabel('time');

ylabel('states');

title('Analytical solution');

hold off;

现在将原问题添加内点约束条件，x1(t1)=2.4,x2(t1)=0.9.

程序如下

clear;

sol_1=dsolve('Dx1=x2,Dx2=-p2,Dp1=0,Dp2=-p1',...

'x1(0)=1,x2(0)=2,x1(t1)=2.4,x2(t1)=0.9');%前半段

sol_2=dsolve('Dx1=x2,Dx2=-p2,Dp1=0,Dp2=-p1',...

'x1(tf)=3,p2(tf)=0,x1(t1)=2.4,x2(t1)=0.9');%后半段

%前半段的系数值

%x11=sol_1.x1;

%x12=sol_1.x2;

p11=sol_1.p1;

p12=sol_1.p2;

%后半段系数值

%极值条件：H对u求导=0

eq0 = subs(sol_2.p1,'t','tf')*subs(sol_2.x2,'t','tf') - 0.5*subs(sol_2.p2,'t','tf')^2;

sol_a=solve(eq0,'tf');

tf=sol_a;

x21=subs(sol_2.x1);

x22=subs(sol_2.x2);

p21=subs(sol_2.p1);

p22=subs(sol_2.p2);

%内点约束的条件：p(t1+)=p(t1-)-n

eq1=p21-p11+'n1';

eq2=subs(p22,'t','t1')-subs(p12,'t','t1')+'n2';

%内点约束条件:H(t1+)=H(t1-)

eq3='0.9*n1'+0.5*subs(p22,'t','t1')*subs(p22,'t','t1')-0.5*subs(p12,'t','t1')*subs(p12,'t','t1');

sol_h=solve(eq1,eq2,eq3);

t1=sol_h.t1;

%t1得出含两个值

t1_1=t1(1,1);

t1_2=t1(2,1);

%由于t1有两个值，分别计算相应目标函数，取最小

tf=subs(tf);

p12=subs(p12);

p22=subs(p22);

tf_1=tf(1,1);

tf_2=tf(2,1);

p12_1=p12(1,1);

p12_2=p12(2,1);

p22_1=p22(1,1);

p22_2=p22(2,1);

J1=int(p12_1^2,0,t1_1)+int(p22_1^2,t1_1,tf_1);

J2=int(p12_2^2,0,t1_2)+int(p22_2^2,t1_2,tf_2);

%J2值最小，故t1值取第二个值t1_2

t1=t1_2;

x11=subs(sol_1.x1);

x12=subs(sol_1.x2);

x21=subs(x21);

x22=subs(x22);

u1=-p12_2;

u2=-p22_2;

t1_1=double(t1);

tf=double(tf_2);

figure(1);

ezplot(x11,[0 t1_1]); hold on;

ezplot(x21,[t1_1 tf]);hold on;

ezplot(x12,[0 t1_1]); hold on;

ezplot(x22,[t1_1 tf]);hold on;

ezplot(u1,[0 t1_1]); hold on;

ezplot(u2,[t1_1 tf]);

axis([0 3 -1.5 3]);

text(1.3,2.5,'x_1(t)');

text(1.3,1,'x_2(t)');

text(1.3,-.5,'u(t)');

xlabel('time');

ylabel('states');

title('Analytical solution');

hold off;

将等式约束改为不等式约束，x1(t1)<=2.4,x2(t1)<=1.

程序如下

function ieq

clear;

sol_1=dsolve('Dx1=x2,Dx2=-p2,Dp1=0,Dp2=-p1',...

'x1(0)=1,x2(0)=2,x1(t1)=2.4-c1^2,x2(t1)=1-c2^2');%前半段

sol_2=dsolve('Dx1=x2,Dx2=-p2,Dp1=0,Dp2=-p1',...

'x1(tf)=3,p2(tf)=0,x1(t1)=2.4-c1^2,x2(t1)=1-c2^2');%后半段

global eq1 eq2 eq3 tf p11 p12 p21 p22;

%前半段的值

p11=sol_1.p1;

p12=sol_1.p2;

%后半段

%极值条件：H对u求导=0

eq0 = subs(sol_2.p1,'t','tf')*subs(sol_2.x2,'t','tf') - 0.5*subs(sol_2.p2,'t','tf')^2;

sol_t=solve(eq0,'tf');

tf=sol_t;%解出tf代入

x21=subs(sol_2.x1);

x22=subs(sol_2.x2);

p21=subs(sol_2.p1);

p22=subs(sol_2.p2);

%内点约束的条件：p(t1+)=p(t1-)-n

eq1=p21-p11+'n1';

eq2=subs(p22,'t','t1')-subs(p12,'t','t1')+'n2';

%内点约束条件:H(t1+)=H(t1-)

eq3='n1*(1-c2^2)'+0.5*subs(p22,'t','t1')*subs(p22,'t','t1')-0.5*subs(p12,'t','t1')*subs(p12,'t','t1');

%要使目标函数极小可得条件n1*c1=0,n2*c2=0,则有四种情况

%情况1：n1=0,n2=0时，无解删除

%情况2：n2=0,c1=0时，得出tf解无穷大删除

%还剩下两种情况：n1=0,c2=0以及c1=0,c2=0如下

[Jam,t1_a,c1_a,c2_a]=mina(0,0);

[Jbm,t1_b,c1_b,c2_b]=minb(0,0);

%取两种情况下最小值及其相应系数

if(Jam>Jbm)

t1=t1_b;

c1=c1_b;

c2=c2_b;

J=Jbm;

else

t1=t1_a;

c1=c1_a;

c2=c2_a;

J=Jam;

end

%代入系数值画图

x11=subs(sol_1.x1);

x12=subs(sol_1.x2);

x21=subs(x21);

x22=subs(x22);

u1=-subs(p12);

u2=-subs(p22);

t1=double(t1);

tf=double(subs(tf));

figure(1);

h=ezplot(x11,[0 t1]); hold on;

set(h,'Color','red');

ezplot(x21,[t1 tf]);hold on;

h1=ezplot(x12,[0 t1]); hold on;

set(h1,'Color','red');

ezplot(x22,[t1 tf]);hold on;

h2=ezplot(u1,[0 t1]); hold on;

set(h2,'Color','red');

ezplot(u2,[t1 tf]);hold on;

axis([0 3 -1.5 3]);

text(1.3,2.5,'x_1(t)');

text(1.3,1,'x_2(t)');

text(1.3,-.5,'u(t)');

xlabel('time');

ylabel('states');

title('Analytical solution');

hold off;

%情况1:n1=0,c2=0

function [Jam,t1,c1,c2]=mina(n1,c2)

global eq1 eq2 eq3 tf p12 p22;

sol_a=solve(subs(eq1),subs(eq2),subs(eq3));

at1=sol_a.t1;

ac1=sol_a.c1;

for i=1:length(at1)

t1=at1(i,1);

c1=ac1(i,1);

tfa(i,1)=subs(tf);

%t1得出含12个值,由已知条件0

if (double(tfa(i,1))<0)

continue

end;

if (double(tfa(i,1))

continue

end;

if(double(at1(i,1))<=0)

continue

end;

p12a=subs(p12);

p22a=subs(p22);

Ja(i,1)=int(p12a^2,0,t1)+int(p22a^2,t1,tfa(i,1));

end

%找到Ja大于零的最小值

Jan=double(Ja);

Jam=min(Jan(Jan>0));

na=find(Jan==Jam);

na=na(1,1);

t1=at1(na,1);

c1=ac1(na,1);

%情况2：c1=0,c2=0

function [Jbm,t1,c1,c2]=minb(c1,c2)

global eq1 eq2 eq3 tf p12 p22;

sol_b=solve(subs(eq1),subs(eq2),subs(eq3));

bt1=sol_b.t1;

%t1只有两个解

for i=1:length(bt1)

t1=bt1(i,1);

tfb(i,1)=subs(tf);

if (double(tfb(i,1))<0) continue

end;

if (double(tfb(i,1))

continue

end;

if(double(bt1(i,1))<=0)

continue

end;

p12b=subs(p12);

p22b=subs(p22);

Jb(i,1)=int(p12b^2,0,t1)+int(p22b^2,t1,tfb(i,1));

end

% %找到Ja大于零的最小值

Jbn=double(Jb);

Jbm=min(Jbn(Jbn>0));

nb=find(Jbn==Jbm);

t1=bt1(nb,1); 多试几组数据发现，只有n1=0和c2=0能得到最优结果