## 25 Reputation

11 years, 330 days

## @  Thank you very much. I'm waiting...

Thank you very much. I'm waiting for yor response.

## @aylin Thanks for your kindness. I appli...

@aylin

Thanks for your kindness. I applied method in the attached paper but I doubt in my results.

please see attached files.

The_improved_GSS1_method1.mws

Optimal_control_and_infectiology.pdf

## @  Thank you very much. Here, ...

Thank you very much. Here, x[1], x[2], x[3] denote the concentration of uninfected hepatocytes, infected
hepatocytes, and free virus, respectively and they are should be positive. u1 and u2 are not parameters.The control functions u1(t) and u2(t) are bounded Lebesgue integrable functions.
The control u1(t) represents the eficiency of drug therapy in blocking
new infection, and the control u2(t) represents the eficiency of drug therapy
in inhibiting viral production.The control set de fined by

{(u1,u2)|u1 and u2 are measurable, 0≤u1≤1 and 0≤u2≤1}

## @   The problem is to maximize...

@

The problem is to maximize the objective functional

J(u1(t); u2(t)) =int(x[1]-A[1]*(u[1])^2/2-A[2]*(u[2])^2/2,t=0..t_f)

where the parameters A[1]≥  0 and A[2]≥  0 are based on the bene ts and
costs of the treatment. Our target is to maximize the objective functional
by increasing the number of the uninfected cells, decreasing
the viral load, and minimizing the cost of treatment.

## @  Thanks alot. The original contro...

Thanks alot. The original control problem is in the following

dH/dt=λ-µ H-(1-u1)β H V+δ I,

dI/dt=(1-u1)β H V-σ I,

‎dV/dt=(1-u2)k I-γ‎V,

Please see the attached maple code.

I've used Pontryagin's maximum principle.

 > restart: unprotect('gamma');
 > L:=x[1]-A[1]*(u[1])^2/2-A[2]*(u[2])^2/2;
 (1)
 > H:=L+psi[1](t)*(lambda-mu*x[1]-(1-u[1])*beta*x[1]*x[3]+delta*x[2])+psi[2](t)*((1-u[1])*beta*x[1]*x[3]-sigma*x[2]) +psi[3](t)*((1-u[2])*k*x[2]-gamma*x[3]);
 (2)
 > du1:=diff(H,u[1]);
 (3)
 > du2:=diff(H,u[2]);
 (4)
 > ddu1:=-A[1]*u[1]+psi[1](t)*beta*x[1]*x[3]-psi[2](t)*beta*x[1]*x[3];
 (5)
 > ddu2:=-A[2]*u[2]-psi[3](t)*k*x[2];
 (6)
 > sol_u1 := solve(ddu1, u[1]);
 (7)
 > sol_u2 := solve(ddu2, u[2]);
 (8)
 > Dx2:=subs(u[1]=beta*x[1]*x[3]*(psi[1](t)-psi[2](t))/A[1],u[2]=-psi[3](t)*k*x[2]/A[2], H );
 (9)
 > ode1:=diff(psi[1](t),t)=-diff(H,x[1]);
 (10)
 > ode2:=diff(psi[2](t),t)=-diff(H,x[2]);
 (11)
 > ode3:=diff(psi[3](t),t)=-diff(H,x[3]);
 (12)
 > restart: #Digits:=10: unprotect('gamma'); lambda:=5*10^5: mu:=0.003: beta:=4*10^(-10): delta:=0: alpha:=0.043: sigma:=alpha+delta: k:=6.24: gamma:=0.65: A[1]:=10: A[2]:=2: u[1]:=beta*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t))/A[1]: u[2]:=-psi[3](t)*k*x[2](t)/A[2]:
 > ics := x[1](0)=1.7*10^8, x[2](0)=0,x[3](0)=400,psi[1](20)=0,psi[2](20)=0,psi[3](20)=0:
 > ode1:=diff(x[1](t), t)=lambda-mu*x[1](t)-(1-u[1])*beta*x[1](t)*x[3](t)+delta*x[2](t), diff(x[2](t), t) =(1-u[1])*beta*x[1](t)*x[3](t)-sigma*x[2](t), diff(x[3](t), t) =(1-u[2])*k*x[2](t)-gamma*x[3](t), diff(psi[1](t), t) =-1+mu*psi[1](t)+beta*x[3](t)*(1-u[1])*(psi[1](t)-psi[2](t)), diff(psi[2](t), t) =delta*psi[1](t)+sigma*psi[2](t)-psi[3](t)*(1-u[2])*k, diff(psi[3](t), t) = beta*x[1](t)*(psi[1](t)-psi[2](t))*(1-u[1])+gamma*psi[3](t);
 (13)
 > sol:=dsolve([ode1,ics],numeric, method = bvp[midrich]);
 >

## transversality conditions...

λ1,λ2,λ3 are transversality conditions where

λi(t_end)=0

and t_end=70.

## full code...

Thanks for your response

my full code is following.

restart:

unprotect('gamma');
x[0]:=140: y[0]:=145:v[0]:=18:
A1:= 900: A2:=1000:
mu:=0.73:
gamma:=250:
beta:=0.0014:
r:=0.01:
delta:=0:
a:=0.0693:
k:=340:

#Step 1
n:= 100:
lambda__1[n]:= 0: lambda__2[n]:= 0: lambda__3[n]:= 0:
u__1[0]:= 0: u__2[0]:= 0:

#Step 2:
for i from 0 to n-1 do
x[i+1]:= x[i]+h*(r*x[i+1]*(1-(x[i+1]+y[i])/k)-(1-u__1[i])*beta*v[i]*(x[i+1]/(x[i+1]+y[i])));
y[i+1]:= y[i]+h*(((1-u__1[i])*beta*v[i]*x[i+1]/(x[i+1]+y[i+1]))-a*y[i+1]);
v[i+1]:= v[i]+h*((1-u__2[i])*gamma*y[i+1]-mu*v[i+1]);
lambda__1[n-i-1]:=lambda__1[n-i]+h*(-1-lambda__1[n-i-1]*(r*(1-(x[i+1]+y[i])/k)
-r*x[i+1]/k-(1-u__1[i])*beta*v[i+1]*y[i+1]/(x[i+1]+y[i])^2)-lambda__2[n-i]*(1-u__1[i])*beta*v[i+1]*y[i+1]/(x[i+1]+y[i])^2);
lambda__2[n-i-1]:= lambda__2[n-i]+h*(lambda__1[n-i-1]*(r*x[i+1]/k-(1-u__1[i])*beta*v[i+1]*x[i+1]/(x[i+1]+y[i])^2)-lambda__2[n-i-1]*((1-u__1[i])*beta*v[i+1]*x[i+1]/(x[i+1]+y[i])^2+a)-lambda__3[n-i]*gamma*(1-u__2[i]));
lambda__3[n-i-1]:= lambda__3[n-i]+h*((1-u__1[i])*beta*x[i+1]/(x[i+1]+y[i])*(lambda__1[n-i-1]-lambda__2[n-i-1])+mu*lambda__3[n-i-1]);
R1:= (1/A1*(x[i+1]+y[i+1]))*(lambda__1[n-i-1]-lambda__2[n-i-1])*beta*v[i+1]*x[i+1];
R2:= -(1/A2)*lambda__3[n-i-1]*gamma*y[i+1];
u__1[i+1]:= min(1, max(R1,0));
u__2[i+1]:= min(1, max(R2,0));
end do:

When I run this code in maple I am facing with "Error, recursive assignment".

## plot...

Thanks for your kindness. How do I plot  functions T,J, V?

## @Preben Alsholm  Dear Preben I writ...

Dear Preben

I write a maple code of the  following paper.

I reached  such a system  odes that I can not solve.

Optimal_Control_of_T.pdf

2.mws

## thanks...

Dear Preben,

Thanks for your kindness but these are not the results in paper that I study.

Sorry I forgot

A[1]=0.2;

A[2]=0.2

Sorry I forgot

A[1]=0.2;

A[2]=0.2

## thanks...

Thanks for the quick response.

But this doesn't look like what I'm looking for.

Variational_iteratio.pdf

Please see Eqs. (34)-(37) in the attached file.

## thanks...

Thanks for the quick response.

But this doesn't look like what I'm looking for.

Variational_iteratio.pdf

Please see Eqs. (34)-(37) in the attached file.

 Page 1 of 1
﻿