Return to Revisions

Cesario's answer is almost there, but I'd like to add a couple points of correction/improvement.

1. The problem formulation has constraints on the control (0<=u<=1) and positivity constraints on the states. I would suggest rather than solving for root where the derivative is zero (equ3 above), use Maximize and add the constraints there:

parms = {s -> 10, m1 -> 0.02, m2 -> 0.5, m3 -> 4.4, r -> 0.03, 
   Tm -> 1500, k -> 0.000024, M -> 300, A -> 1, T0 -> 800, 
   Ti0 -> 0.04, V0 -> 1.5};
x1 = T[t];
x2 = Ti[t];
x3 = V[t];
X = {x1, x2, x3};
Lambda = {la1[t], la2[t], la3[t]};
U = u[t];
f1 = s/(1 + x3) - m1 x1 + r x1 (1 - (x1 + x2)/Tm) - U k x1 x3;
f2 = U k x3 x1 - m2 x2;
f3 = M m2 x2 - m3 x3;
F = {f1, f2, f3};
L = A x1 - (1 - U)^2;
H = L + Lambda . F;
equ1 = Thread[D[X, t] == Grad[H, Lambda]];
equ2 = Thread[D[Lambda, t] == -Grad[H, X]];
(*equ3=D[H,U]==0;
solU=Solve[equ3,U][[1]];*)
cons = {0 <= u[t] <= 1 && T[t] > 0 && Ti[t] > 0 && V[t] > 0, k > 0};
equ3 = Last@Maximize[{H, cons}, U];
solU = FullSimplify[equ3, cons]

This (nearly) yields the piecewise function reported in the original work. Please note that the third ODE in the SE problem is missing a term found in the original problem, but the idea is the same.

2. The boundary conditions for lambda are given as lambda[0] in Cesario's answer, but these should apply to the end time tf, which is the gradient of the end cost (Mayer term) and here {0,0,0}. Programmatically:

M=0; (*Mayer cost*)
tf=100;  (* from paper *)
cinitsLambda = Thread[(Lambda /. {t -> tf}) == Grad[M, X /. t -> tf]];

Now the setup for NDSolve can be completed:

equs = Join[equ1, equ2];
cinits = Join[cinitsX, cinitsLambda];
vars = Join[X, Lambda];
DE = Join[equs, cinits] /. solU /. parms;

3. And you might hope that the following NDSolve statement would work:

solDE = NDSolve[DE, vars, {t, 0, tf}];
Plot[Evaluate[vars /. solDE], {t, 0, tf}]

But alas it does not seem to converge in an overnight run, nor does the original equations from the paper. (vide supra)

I know that BVPs are tough problems to solve, so I tried an alternative implementation in my personal fork of PSOPT pseudospectral code, as well as the trapezoidal method in Mathematica. These are direct solvers (vs. the Pontryagin indirect method) and also have problems with these equations. I have come to conclusion that there's something fishy or extreme stiffness.

I have tried this code on numerous multivariate optimal control problems and it seems to yield answers from Ross, Bryson and Ho, et al.