I am unable to solve the attached optimal control problem,please any one who many help me in guideing .tnx
restart:
unprotect('gamma');
L:=b[1]*c(t)+b[2]*i(t)+w[1]*(u[1])^2/2+w[2]*(u[2])^2/2+w[3]*(u[3])^2/2;
1 2 1 2 1 2
b[1] c(t) + b[2] i(t) + - w[1] u[1] + - w[2] u[2] + - w[3] u[3]
2 2 2
H:=L+lambda[1](t)*((1-p*Psi)*tau+phi* v + delta *r-lambda*(1-u[3])*s-u[1]*varphi*s -mu*s ) +lambda[2](t)*(p*Psi*tau + u[1]*vartheta*s -gamma*lambda* (1-u[3])*v-(mu+phi)*v ) +lambda[3](t)*( (1-u[3])*rho*lambda* (s +gamma*v)+(1-q)* u[2]*eta*i -(mu +beta +chi)*c ) +lambda[4](t)* ((1-rho)*(1-u[3])*lambda*( s +gamma*v) +chi*c - u[2]*eta*i - (mu +alpha )*i) +lambda[5](t)*( beta*c + u[2]*q*eta*i -(mu +delta)*r);
1 2 1 2 1 2
b[1] c(t) + b[2] i(t) + - w[1] u[1] + - w[2] u[2] + - w[3] u[3] + lambda[1](t
2 2 2
) ((1 - p Psi) tau + phi v + delta r - lambda (1 - u[3]) s - u[1] varphi s
- mu s) + lambda[2](t) (p Psi tau + u[1] vartheta s
- gamma lambda (1 - u[3]) v - (mu + phi) v) + lambda[3](t) ((1 - u[3]) rho
lambda (s + gamma v) + (1 - q) u[2] eta i - (mu + beta + chi) c) + lambda[4](t
) ((1 - rho) (1 - u[3]) lambda (s + gamma v) + chi c - u[2] eta i
- (mu + alpha) i) + lambda[5](t) (beta c + u[2] q eta i - (mu + delta) r)
du1:=diff(H,u[1]);
w[1] u[1] - lambda[1](t) varphi s + lambda[2](t) vartheta s
du2:=diff(H,u[2]);du3:=diff(H,u[3]);
w[2] u[2] + lambda[3](t) (1 - q) eta i - lambda[4](t) eta i
+ lambda[5](t) q eta i
w[3] u[3] + lambda[1](t) lambda s + lambda[2](t) gamma lambda v
- lambda[3](t) rho lambda (s + gamma v)
- lambda[4](t) (1 - rho) lambda (s + gamma v)
ddu1 := -A[1] u[1] + psi[1](t) beta x[1] x[3] - psi[2](t) beta x[1] x[3]
ddu2 := -A[2] u[2] - psi[3](t) k x[2]
sol_u1 := solve(du1, u[1]);
s(t) (lambda[1](t) varphi - lambda[2](t) vartheta)
--------------------------------------------------
w[1]
sol_u2 := solve(du2, u[2]);sol_u3 := solve(du3, u[3]);
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q)
----------------------------------------------------------------------
w[2]
1
---- (lambda (-lambda[1](t) s - lambda[2](t) gamma v + lambda[3](t) rho s
w[3]
+ lambda[3](t) rho gamma v + lambda[4](t) s + lambda[4](t) gamma v
- lambda[4](t) rho s - lambda[4](t) rho gamma v))
Dx2:=subs(u[1]= s*(lambda[1](t)*varphi-lambda[2](t)*vartheta)/w[1] ,u[2]= eta*i*(-lambda[3](t)+lambda[3](t)*q+lambda[4](t)-lambda[5](t)*q)/w[2], u[3]=-lambda*(lambda[1](t)*s+lambda[2](t)*gamma*v-lambda[3](t)*rho*s-lambda[3](t)*rho*gamma*v-lambda[4](t)*s-lambda[4](t)*gamma*v+lambda[4](t)*rho*s+lambda[4](t)*rho*gamma*v)/w[3] ,H );
2 2
s (lambda[1](t) varphi - lambda[2](t) vartheta)
b[1] c(t) + b[2] i(t) + -------------------------------------------------
2 w[1]
2 2 2
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q)
+ ------------------------------------------------------------------------- +
2 w[2]
1 / 2
------ \lambda (lambda[1](t) s + lambda[2](t) gamma v - lambda[3](t) rho s
2 w[3]
- lambda[3](t) rho gamma v - lambda[4](t) s - lambda[4](t) gamma v
/
\ |
+ lambda[4](t) rho s + lambda[4](t) rho gamma v)^2/ + lambda[1](t) |(1
\
/ 1
- p Psi) tau + phi v + delta r - lambda |1 + ---- (lambda (lambda[1](t) s
\ w[3]
+ lambda[2](t) gamma v - lambda[3](t) rho s - lambda[3](t) rho gamma v
- lambda[4](t) s - lambda[4](t) gamma v + lambda[4](t) rho s
\
+ lambda[4](t) rho gamma v))| s
/
2 \
s (lambda[1](t) varphi - lambda[2](t) vartheta) varphi |
- ------------------------------------------------------- - mu s| +
w[1] /
/
|
lambda[2](t) |p Psi tau
\
2
s (lambda[1](t) varphi - lambda[2](t) vartheta) vartheta /
+ --------------------------------------------------------- - gamma lambda |1 +
w[1] \
1
---- (lambda (lambda[1](t) s + lambda[2](t) gamma v - lambda[3](t) rho s
w[3]
- lambda[3](t) rho gamma v - lambda[4](t) s - lambda[4](t) gamma v
\
\ |
+ lambda[4](t) rho s + lambda[4](t) rho gamma v))| v - (mu + phi) v| +
/ /
// 1
lambda[3](t) ||1 + ---- (lambda (lambda[1](t) s + lambda[2](t) gamma v
\\ w[3]
- lambda[3](t) rho s - lambda[3](t) rho gamma v - lambda[4](t) s
\
- lambda[4](t) gamma v + lambda[4](t) rho s + lambda[4](t) rho gamma v))|
/
1 / 2 2
rho lambda (s + gamma v) + ---- \(1 - q) eta i (-lambda[3](t)
w[2]
\ \
+ lambda[3](t) q + lambda[4](t) - lambda[5](t) q)/ - (mu + beta + chi) c| +
/
/
| / 1
lambda[4](t) |(1 - rho) |1 + ---- (lambda (lambda[1](t) s
\ \ w[3]
+ lambda[2](t) gamma v - lambda[3](t) rho s - lambda[3](t) rho gamma v
- lambda[4](t) s - lambda[4](t) gamma v + lambda[4](t) rho s
\
+ lambda[4](t) rho gamma v))| lambda (s + gamma v) + chi c
/
2 2
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q)
- ------------------------------------------------------------------------
w[2]
\ /
| |
- (mu + alpha) i| + lambda[5](t) |beta c
/ \
+
2 2
eta i (-lambda[3](t) + lambda[3](t) q + lambda[4](t) - lambda[5](t) q) q
--------------------------------------------------------------------------
w[2]
\
|
- (mu + delta) r|
/
ode1:=diff(lambda[1](t),t)=-diff(H,s);ode2:=diff(lambda[2](t),t)=-diff(H,v);ode3:=diff(psi[3](t),t)=-diff(H,c);ode4:=diff(lambda[4](t),t)=-diff(H,i);ode5:=diff(lambda[5](t),t)=-diff(H,r);
d
--- lambda[1](t) = -lambda[1](t) (-lambda (1 - u[3]) - u[1] varphi - mu)
dt
- lambda[2](t) u[1] vartheta - lambda[3](t) (1 - u[3]) rho lambda
- lambda[4](t) (1 - rho) (1 - u[3]) lambda
d
--- lambda[2](t) = -lambda[1](t) phi
dt
- lambda[2](t) (-gamma lambda (1 - u[3]) - mu - phi)
- lambda[3](t) (1 - u[3]) rho lambda gamma
- lambda[4](t) (1 - rho) (1 - u[3]) lambda gamma
d
--- psi[3](t) = -lambda[3](t) (-mu - beta - chi) - lambda[4](t) chi
dt
- lambda[5](t) beta
d
--- lambda[4](t) = -lambda[3](t) (1 - q) u[2] eta
dt
- lambda[4](t) (-u[2] eta - mu - alpha) - lambda[5](t) u[2] q eta
d
--- lambda[5](t) = -lambda[1](t) delta - lambda[5](t) (-mu - delta)
dt
restart:
#Digits:=10:
unprotect('gamma');
lambda:=0.51:
mu:=0.002:
beta:=0.0115:
delta:=0.003:
alpha:=0.33:
chi:=0.00274:
k:=6.24:
gamma:=0.4:
rho:=0.338:;tau=1000:;Psi:=0.1:;p:=0.6:;phi:=0.001:;eta:=0.001124:q:=0.6:varphi:=0.9:;vatheta:=0.9:
b[1]:=2:;b[2]:=3:;w[1]:=4:;w[2]:=5:;w[3]:=6:
#u[1]:=s(t)*(lambda[1](t)*varphi-lambda[2](t)*vartheta)/w[1]:
#u[2]:=eta*i*(-lambda[3](t)+lambda[3](t)*q+lambda[4](t)-lambda[5](t)*q)/w[2]:;u[3]:=lambda*(-lambda[1](t)*s-lambda[2](t)*gamma*v+lambda[3](t)*rho*s+lambda[3](t)*rho*gamma*v+lambda[4](t)*s+lambda[4](t)*gamma*v-lambda[4](t)*rho*s-lambda[4](t)*rho*gamma*v)/w[3]:
ics := s(0)=8200, v(0)=2800,c(0)=1100,i(0)=1500,r(0)=200,lambda[1](20)=0,lambda[2](20)=0,lambda[3](20)=0,lambda[4](20)=0,lambda[5](20)=0:
ode1:=diff(s(t),t)=(1-p*Psi)*tau+phi* v(t) + delta *r(t)-lambda*(1-u[3])*s(t)-u[1]*varphi*s(t) -mu*s(t),
diff(v(t), t) =p*Psi*tau + u[1]*vartheta*s(t) -gamma*lambda* (1-u[3])*v(t)-(mu+phi)*v(t) ,
diff(c(t), t) =(1-u[3])*rho*lambda* (s(t) +gamma*v(t))+(1-q)* u[2]*eta*i(t) -(mu +beta +chi)*c(t),
diff(i(t), t) =(1-rho)*(1-u[3])*lambda*( s(t) +gamma*v(t)) +chi*c(t) - u[2]*eta*i(t) - (mu +alpha )*i(t),
diff(r(t), t) = beta*c(t) + u[2]*q*eta*i(t) -(mu +delta)*r(t),
diff(lambda[1](t), t) = -lambda[1](t)*(-lambda*(1-u[3])-u[1]*varphi-mu)-lambda[2](t)*u[1]*vartheta-lambda[3](t)*(1-u[3])*rho*lambda-lambda[4](t)*(1-rho)*(1-u[3])*lambda,diff(lambda[2](t),t)=-lambda[1](t)*phi-lambda[2](t)*(-gamma*lambda*(1-u[3])-mu-phi)-lambda[3](t)*(1-u[3])*rho*lambda*gamma-lambda[4](t)*(1-rho)*(1-u[3])*lambda*gamma,diff(lambda[3](t),t)=-lambda[3](t)*(-mu-beta-chi)-lambda[4](t)*chi-lambda[5](t)*beta,diff(lambda[4](t),t)=-lambda[3](t)*(1-q)*u[2]*eta-lambda[4](t)*(-u[2]*eta-mu-alpha)-lambda[5](t)*u[2]*q*eta,diff(lambda[5](t),t)=-lambda[1](t)*delta-lambda[5](t)*(-mu-delta);
d
--- s(t) = (1 - p Psi) tau + phi v(t) + delta r(t) - lambda (1 - u[3]) s(t)
dt
d
- u[1] varphi s(t) - mu s(t), --- v(t) = p Psi tau + u[1] vartheta s(t)
dt
d
- gamma lambda (1 - u[3]) v(t) - (mu + phi) v(t), --- c(t) = (1 - u[3]) rho lambda
dt
(s(t) + gamma v(t)) + (1 - q) u[2] eta - (mu + beta + chi) c(t), 0 = (1
- rho) (1 - u[3]) lambda (s(t) + gamma v(t)) + chi c(t) - u[2] eta - mu
d d
- alpha, --- r(t) = beta c(t) + u[2] q eta - (mu + delta) r(t), ---
dt dt
lambda[1](t) = -lambda[1](t) (-lambda (1 - u[3]) - u[1] varphi - mu)
- lambda[2](t) u[1] vartheta - lambda[3](t) (1 - u[3]) rho lambda
d
- lambda[4](t) (1 - rho) (1 - u[3]) lambda, --- lambda[2](t) =
dt
-lambda[1](t) phi - lambda[2](t) (-gamma lambda (1 - u[3]) - mu - phi)
- lambda[3](t) (1 - u[3]) rho lambda gamma
d
- lambda[4](t) (1 - rho) (1 - u[3]) lambda gamma, --- lambda[3](t) =
dt
d
-lambda[3](t) (-mu - beta - chi) - lambda[4](t) chi - lambda[5](t) beta, ---
dt
lambda[4](t) = -lambda[3](t) (1 - q) u[2] eta
- lambda[4](t) (-u[2] eta - mu - alpha) - lambda[5](t) u[2] q eta,
d
--- lambda[5](t) = -lambda[1](t) delta - lambda[5](t) (-mu - delta)
dt
sol := dsolve({c(0) = 0, i(0) = 0, r(0) = .1, s(0) = 0, v(0) = 0, diff(c(t), t) = (1-u[3])*rho*lambda*(s(t)+gamma*v(t))+(1-q)*u[2]*eta*i(t)-(mu+beta+chi)*c(t), diff(i(t), t) = (1-rho)*(1-u[3])*lambda*(s(t)+gamma*v(t))+chi*c(t)-u[2]*eta*i(t)-(mu+alpha)*i(t), diff(r(t), t) = beta*c(t)+u[2]*q*eta*i(t)-(mu+delta)*r(t), diff(s(t), t) = (1-p*Psi)*tau+phi*v(t)+delta*r(t)-lambda*(1-u[3])*s(t)-u[1]*varphi*s(t)-mu*s(t), diff(v(t), t) = p*Psi*tau+u[1]*vartheta*s(t)-gamma*lambda*(1-u[3])*v(t)-(mu+phi)*v(t), diff(lambda[1](t), t) = -lambda[1](t)*(-lambda*(1-u[3])-u[1]*varphi-mu)-lambda[2](t)*u[1]*vartheta-lambda[3](t)*(1-u[3])*rho*lambda-lambda[4](t)*(1-rho)*(1-u[3])*lambda, diff(lambda[2](t), t) = -lambda[1](t)*phi-lambda[2](t)*(-gamma*lambda*(1-u[3])-mu-phi)-lambda[3](t)*(1-u[3])*rho*lambda*gamma-lambda[4](t)*(1-rho)*(1-u[3])*lambda*gamma, diff(lambda[3](t), t) = -lambda[3](t)*(-mu-beta-chi)-lambda[4](t)*chi-lambda[5](t)*beta, diff(lambda[4](t), t) = -lambda[3](t)*(1-q)*u[2]*eta-lambda[4](t)*(-u[2]*eta-mu-alpha)-lambda[5](t)*u[2]*q*eta, diff(lambda[5](t), t) = -lambda[1](t)*delta-lambda[5](t)*(-mu-delta), lambda[1](20) = 0, lambda[2](20) = 0, lambda[3](20) = 0, lambda[4](20) = 0, lambda[5](20) = 0}, type = numeric);
Error, (in dsolve/numeric/process_input) invalid specification of initial conditions, got 1 = 0
sol:=dsolve([ode1,ics],numeric, method = bvp[midrich],maxmesh=500);
Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations
dsolve[':-interactive']({});
Error, `:=` unexpected
sol:=dsolve([ode1,ics],numeric, method = bvp[midrich],maxmesh=500);
Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations
eq1:=diff(s(t), t)=(1-p*Psi)*tau+phi* v(t) + delta *r(t)-lambda*(1-u[3])*s(t)-u[1]*varphi*s(t) -mu*s(t);
eq2:diff(v(t), t) =p*Psi*tau + u[1]*vartheta*s(t) -gamma*lambda* (1-u[3])*v(t)-(mu+phi)*v(t);
eq3:=diff(c(t), t) =(1-u[3])*rho*lambda* (s(t) +gamma*v(t))+(1-q)* u[2]*eta*i(t) -(mu +beta +chi)*c(t);
eq4:=diff(i(t), t) =(1-rho)*(1-u[3])*lambda*( s(t) +gamma*v(t)) +chi*c(t) - u[2]*eta*i(t) - (mu +alpha )*i(t);
eq5:=diff(r(t), t) = beta*c(t) + u[2]*q*eta*i(t) -(mu +delta)*r(t);
d
--- s(t) = (1 - p Psi) tau + phi v(t) + delta r(t) - lambda (1 - u[3]) s(t)
dt
- u[1] varphi s(t) - mu s(t)
d
--- v(t) = p Psi tau + u[1] vartheta s(t) - gamma lambda (1 - u[3]) v(t)
dt
- (mu + phi) v(t)
d
--- c(t) = (1 - u[3]) rho lambda (s(t) + gamma v(t)) + (1 - q) u[2] eta i(t)
dt
- (mu + beta + chi) c(t)
d
--- i(t) = (1 - rho) (1 - u[3]) lambda (s(t) + gamma v(t)) + chi c(t)
dt
- u[2] eta i(t) - (mu + alpha) i(t)
d
--- r(t) = beta c(t) + u[2] q eta i(t) - (mu + delta) r(t)
dt
eq6:=diff(Q(t),t)=b[1]*c(t)+b[2]*i(t)+w[1]*(u[1])^2/2+w[2]*(u[2])^2/2+w[3]*(u[3])^2/2;
d 1 2 1 2 1 2
--- Q(t) = b[1] c(t) + b[2] i(t) + - w[1] u[1] + - w[2] u[2] + - w[3] u[3]
dt 2 2 2
ics:=s(0)=8200, v(0)=2800,c(0)=1100,i(0)=1500,r(0)=200,Q(0)=6700;
s(0) = 8200, v(0) = 2800, c(0) = 1100, i(0) = 1500, r(0) = 200, Q(0) = 6700
sol0:=dsolve({eq1,eq2,eq3,eq4,eq5,eq6,ics},type=numeric,stiff=true,'parameters'=[u[1],u[2],u[3]],abserr=1e-15,relerr=1e-12,maxfun=0,range=0..50):
Error, (in dsolve/numeric/process_input) system must be entered as a set/list of expressions/equations
with(plots):
Q0:=6700;
6700
obj:=proc(u)
global sol0,Q0;
local ob1;
try
sol0('parameters'=[u[1],u[2],u[3]]):
ob1:=subs(sol0(20.),Q(t)):
catch :
ob1:=0;
end try;
#ob1:=subs(sol0(20.),Q(t));
if ob1>Q0 then Q0:=ob1;print(Q0,u);end;
ob1;
end proc;
proc(u) ... end;
obj([1,1,1]);
0
obj([3,2.5],4);
0
u0:=Vector(3,[0.,0.,0.],datatype=float[8]);
Vector[column](%id = 85973880)
Q0:=0;
Q0 := 0
with(Optimization);
[ImportMPS, Interactive, LPSolve, LSSolve, Maximize, Minimize, NLPSolve,
QPSolve]
sol2:=NLPSolve(3,obj,initialpoint=u0,method=nonlinearsimplex,maximize,evaluationlimit=100):
sol0('parameters'=[3.18125786060723, 2.36800986932868]);
sol0(parameters = [3.18125786060723, 2.36800986932868])
for i from 1 to 3 do odeplot(sol0,[t,x[i](t)],0..20,thickness=3,axes=boxed);od;
Error, (in plots/odeplot) input is not a valid dsolve/numeric solution
