code for bang bang control

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

Answer 1

[Ameer Hamza]

Use multiplication operator after P4

P4*(V^2-U*(1-U))
 %^

[Zofit Allouche]

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]

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]

You may apply chain rule here

Therefore, change the line to

DP4 = -1*diff(H,U)*(1/diff(V,U));

[Zofit Allouche]

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]

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]

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]

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]

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

[Santanu Bhattacharya]

Did you solve this problem?

I also have some similar problem, can we discuss about it?