aylin

20 Reputation

3 Badges

10 years, 64 days

MaplePrimes Activity


These are replies submitted by aylin

@ 

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

@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, 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 eficiency of drug therapy in blocking
new infection, and the control u2(t) represents the eficiency 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 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 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;

L := x[1]-(1/2)*(A[1]*u[1]^2)-(1/2)*(A[2]*u[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]);

H := x[1]-(1/2)*(A[1]*u[1]^2)-(1/2)*(A[2]*u[2]^2)+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]);

du1 := -A[1]*u[1]+psi[1](t)*beta*x[1]*x[3]-psi[2](t)*beta*x[1]*x[3]

(3)

du2:=diff(H,u[2]);

du2 := -A[2]*u[2]-psi[3](t)*k*x[2]

(4)

ddu1:=-A[1]*u[1]+psi[1](t)*beta*x[1]*x[3]-psi[2](t)*beta*x[1]*x[3];

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];

ddu2 := -A[2]*u[2]-psi[3](t)*k*x[2]

(6)

sol_u1 := solve(ddu1, u[1]);

sol_u1 := beta*x[1]*x[3]*(psi[1](t)-psi[2](t))/A[1]

(7)

sol_u2 := solve(ddu2, u[2]);

sol_u2 := -psi[3](t)*k*x[2]/A[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 );

Dx2 := x[1]-(1/2)*(beta^2*x[1]^2*x[3]^2*(psi[1](t)-psi[2](t))^2/A[1])-(1/2)*(psi[3](t)^2*k^2*x[2]^2/A[2])+psi[1](t)*(lambda-mu*x[1]-(1-beta*x[1]*x[3]*(psi[1](t)-psi[2](t))/A[1])*beta*x[1]*x[3]+delta*x[2])+psi[2](t)*((1-beta*x[1]*x[3]*(psi[1](t)-psi[2](t))/A[1])*beta*x[1]*x[3]-sigma*x[2])+psi[3](t)*((1+psi[3](t)*k*x[2]/A[2])*k*x[2]-gamma*x[3])

(9)

ode1:=diff(psi[1](t),t)=-diff(H,x[1]);

ode1 := diff(psi[1](t), t) = -1-psi[1](t)*(-mu-(1-u[1])*beta*x[3])-psi[2](t)*(1-u[1])*beta*x[3]

(10)

ode2:=diff(psi[2](t),t)=-diff(H,x[2]);

ode2 := diff(psi[2](t), t) = -psi[1](t)*delta+psi[2](t)*sigma-psi[3](t)*(1-u[2])*k

(11)

ode3:=diff(psi[3](t),t)=-diff(H,x[3]);

ode3 := diff(psi[3](t), t) = psi[1](t)*(1-u[1])*beta*x[1]-psi[2](t)*(1-u[1])*beta*x[1]+psi[3](t)*gamma

(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);

diff(x[1](t), t) = 500000-0.3e-2*x[1](t)-(1/2500000000)*(1-(1/25000000000)*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t)))*x[1](t)*x[3](t), diff(x[2](t), t) = (1/2500000000)*(1-(1/25000000000)*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t)))*x[1](t)*x[3](t)-0.43e-1*x[2](t), diff(x[3](t), t) = 6.24*(1+3.120000000*psi[3](t)*x[2](t))*x[2](t)-.65*x[3](t), diff(psi[1](t), t) = -1+0.3e-2*psi[1](t)+(1/2500000000)*x[3](t)*(1-(1/25000000000)*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t)))*(psi[1](t)-psi[2](t)), diff(psi[2](t), t) = 0.43e-1*psi[2](t)-6.24*psi[3](t)*(1+3.120000000*psi[3](t)*x[2](t)), diff(psi[3](t), t) = (1/2500000000)*x[1](t)*(psi[1](t)-psi[2](t))*(1-(1/25000000000)*x[1](t)*x[3](t)*(psi[1](t)-psi[2](t)))+.65*psi[3](t)

(13)

sol:=dsolve([ode1,ics],numeric, method = bvp[midrich]);

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

 

 


Download mapleprime.mws

 

@Carl Love 

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

λi(t_end)=0 

and t_end=70.

 

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".

Please help me to solve this equations on Maple.

@Carl Love 

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

 

 

@Preben Alsholm 

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

@Preben Alsholm 

Dear Preben,

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

@Carl Love 

Sorry I forgot

A[1]=0.2;

A[2]=0.2

@Preben Alsholm 

Sorry I forgot

A[1]=0.2;

A[2]=0.2

@Markiyan Hirnyk 

 

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.

@Markiyan Hirnyk 

 

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