code for bang bang control

2020 6 14
1 个回答
更新时间：2022 5 30
8 次查看（30 天）

0 个投票

I have an optimal control problem with a bounded control : u is between [0,1]. in order to solve it I have use a "bang bang contorl",e.g. I have used the following trick:
V^2=(U+0)*(1-U);
and then created the following code:
syms V U S I R P1 P2 P3 P4;
beta=3/11;
gamma=1/11;
DS= -beta*U*S*I;
DI= beta*U*S*I - gamma*I;
DR= gamma*I
V=sqrt((U+0)*(1-U));
g=0.5*I^2 + 0.5*U^2;
H = g + P1*DS + P2*DI+ P3*DR+P4(V^2-U*(1-U));
However when trying to run this the system gives the following error:
"Array indices must be positive integers or logical values
error in syms/subsref
R_tilde=builtin ('subsref,L_tilde,Idx)
Error in line 9
H = g + P1*DS + P2*DI+ P3*DR+P4(V^2-U*(1-U))
I would like to mention that when removing the part of P4(V^2-U*(1-U)), it runs without Error.
Can someone explain to me how to solve the problem? please advice.
Zofit

请先登录，再回答此问题。

回答（1 个）

Ameer Hamza
Ameer Hamza 2020-6-14

0 个投票

Use multiplication operator after P4
P4*(V^2-U*(1-U))
%^

9 个评论
显示 7更早的评论 隐藏 7更早的评论

Zofit Allouche
Zofit Allouche 2020-6-15
编辑：Zofit Allouche 2020-6-15
First of all many thanks - indeed a sign of multiplication was missing...
I did enter it and then the comment disapeared (code was running without any errors) .
However - when entering the differential equations and arriveing to the differential equation for DP4 (see the following code):
syms V U S I R P1 P2 P3 P4;
beta=3/11;
gamma=1/11;
DS= -beta*U*S*I;
DI= beta*U*S*I - gamma*I;
DR= gamma*I
V=sqrt((U+0)*(1-U));
g=0.5*I^2 + 0.5*U^2;
H = g + P1*DS + P2*DI+ P3*DR+P4*(V^2-U*(1-U));
DP1 = -1*diff(H,S);
DP2= -1*diff(H,I);
DP3 = -1*diff(H,R);
DP4 = -1*diff(H,V);
DU=diff(H,U);
for DP4 = -1*diff(H,V) the following error appeared:
"Error using sym/diff (line 70)
Second argument must be a variable or a nonnegative integer specifying the number of differentiations.
Error in trialrun (line 13)
DP4 = -1*diff(H,V);"
Can you please advice how to solve this problem?
Many thanks
Zofit
Walter Roberson
Walter Roberson 2020-6-15
V=sqrt((U+0)*(1-U));
DP4 = -1*diff(H,V);
V is a symbolic expression, not a simple variable. You cannot differentiate a function with respect to an expression (or function)
Ameer Hamza
Ameer Hamza 2020-6-15
You may apply chain rule here
Therefore, change the line to
DP4 = -1*diff(H,U)*(1/diff(V,U));
Zofit Allouche
Zofit Allouche 2020-6-15
编辑：Ameer Hamza 2020-6-15
many many thanks!.
Indeed now the code is running but the matlab gives a warning that its can't find symbolic solution.
maybe you can understand why?
the code is:
syms V U S I R P1 P2 P3 P4;
beta=3/11;
gamma=1/11;
DS= -beta*U*S*I;
DI= beta*U*S*I - gamma*I;
DR= gamma*I
V=sqrt((U+0)*(1-U));
g=0.5*I^2 + 0.5*U^2;
H = g + P1*DS + P2*DI+ P3*DR+P4*(V^2-U*(1-U));
DP1 = -1*diff(H,S);
DP2= -1*diff(H,I);
DP3 = -1*diff(H,R);
DP4 = -1*diff(H,U)*(1/diff(V,U));
DU=diff(H,U);
sol_u = solve(DU,U);
DS = subs(DS,U,sol_u);
DI= subs(DI,U,sol_u);
eq1 = strcat('DS=',char(DS));
eq2 = strcat('DI=',char(DI));
eq3 = strcat('DR=',char(DR));
eq4 = strcat('DP1=',char(DP1));
eq5 = strcat('DP2=',char(DP2));
eq6 = strcat('DP3=',char(DP3));
eq7 = strcat('DP4=',char(DP4));
sol_h = dsolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7);
conA1 = 'S(0) = 1-(6/8740000)';
conA2 = 'I(0) = 6/8740000';
conA3 = 'R(0) = 0';
conA4 = 'P1(365) = 0';
conA5 = 'P2(365) = 0';
conA6 = 'P3(365) = 0';
conA7 = 'P4(365) = 0';
sol_a = dsolve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,conA1,conA2,conA3,conA4,conA5,conA6,conA7);
matlab is warning that :
"Warning: Support of character vectors and strings will be removed in a future release. Use sym objects to define differential equations instead.
> In dsolve (line 126)
In trialrun (line 25)
Warning: Unable to find symbolic solution.
> In dsolve (line 216)
In trialrun (line 25)
Warning: Support of character vectors and strings will be removed in a future release. Use sym objects to define differential equations instead.
> In dsolve (line 126)
In trialrun (line 33)
Warning: Unable to find symbolic solution.
> In dsolve (line 216)
In trialrun (line 33) "
Can you advise for solution?
Zofit
Ameer Hamza
Ameer Hamza 2020-6-15
What are you trying to do in these lines?
eq1 = strcat('DS=',char(DS));
eq2 = strcat('DI=',char(DI));
eq3 = strcat('DR=',char(DR));
eq4 = strcat('DP1=',char(DP1));
eq5 = strcat('DP2=',char(DP2));
eq6 = strcat('DP3=',char(DP3));
eq7 = strcat('DP4=',char(DP4));
Can you show your equations in mathematical form?
Zofit Allouche
Zofit Allouche 2020-6-15
编辑：Zofit Allouche 2020-6-15
What I am trying to do is the following:
% 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);
In order to solve the differential equations
Ameer Hamza
Ameer Hamza 2020-6-15
dsolve() can take symbolic equations as input.
Also, what I meant that what equation are you trying to solve? It seems that there is some confusion about how to feed those equations to dsolve(). Can you show the actual equations?
Zofit Allouche
Zofit Allouche 2020-6-15
编辑：Walter Roberson 2020-6-15
yes. Ofcourse
The equations are:
State equations are:
d/dt S= -uβSi
d/dt i=uβSi- γi
d/dt R= γi
with the following boundry conditions:
S(t=0) = 1- 6/8740000
I (t=0) = 6/8740000
R(t=0) =0
and there is a restricition on the control:
0< u<1
I have used the folloing trick for u:
ν^2=(u+0)(1-u)
I have the following functional :
J(u) = integral between 0 to 365 of ((0.5*i(t)^2 + 0.5*u(t)^2) dt'
and after entering in the state equation into it with lagrange mulitplier I have created the following hamiltonian:
= g(x(t),u(t),t)+ p(t)*(a(x(t),u(t),t)
=0.5*I(t)^2 + 0.5*u(t)^2+P1*(-uβSi- dS/dt) + P2*(uβSi- γi-dI/dt) + P3*(γi-dR/dt) + P4*(v^2-u(t)(1-u(t))
and then created the following olier lagrange equations from it:
  1. dH/dS- d/dt (dH/dS')=0
  2. dH/dI- d/dt (dH/dI')=0
  3. dH/dR- d/dt (dH/dR')=0
  4. dH/du- d/dt(dH/du')=0
  5. dH/dv- d/dt(dH/dv')=0
and the following boundry conditions:
P1(356)=0
P2(365) =0
P3(365)=0
I am afraid that as we are dealign with non linera differential equations, only noumerical method can work here (like BVP) - am I right?
Zofit
Did you solve this problem?
I also have some similar problem, can we discuss about it?

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Physics 的更多信息

产品

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by