function [ t, x, u, l, cost ] = opt08 ( s, m1, m2, m3, r, tmax, k, p, T0, ...
  Ti0, V0, A, T, N )

%*****************************************************************************80
%
%% opt08() performs optimal control for case #08.
%
%  Discussion:
%
%    HIV Treatment
%
%    maximize
%      integral ( 0 <= t <= T ) A T(t) - ( 1 - u(t) )^2 dt
%    subject to
%      T'(t) = s / ( 1 + V(t) ) - m1 T(t) 
%        + r T(t) ( 1 - ( T(t) - Ti(t) ) / Tmax ) - u(t) k V(t) T(t)
%      Ti'(t) = u(t) k V(t) T(t) - m2 Ti(t)
%      V'(t) = P m2 Ti(t) - m3 V(t)
%      T(0) = T0
%      Ti(0) = Ti0
%      V(0) = V0
%    with Hamiltonian:
%      H = A T(t) - ( 1 - u(t) )^2 
%        + lambda1(t) ( s / ( 1 + V(t) ) - m1 T(t) 
%          + r T(t) ( 1 - ( T(t) - Ti(t) ) / Tmax ) - u(t) k V(t) T(t) )
%        + lambda2(t) ( u(t) k V(t) T(t) - m2 Ti(t) )
%        + lambda3(t) ( P m2 Ti(t) - m3 V(t) )
%    with adjoint equations:
%      lambda1(t)' = - A - lambda1(t) [ -m1 + r ( 1 - ( T(t) - Ti(t) ) / Tmax ) 
%        - r T(t) / Tmax - u(t) k V(t) ] - lambda2(t) u(t) k V(t)
%      lambda1(tfinal) = 0.0
%      lambda2(t)' = lambda1(t) r T(t) / Tmax + m2 lambda2(t) - lambda3(t) P m2
%      lambda2(tfinal) = 0.0
%      lambda3(t)' = - lambda1(t) [ - s / ( 1 + V(t) )^2 - u(t) k T(t) ]
%        - lambda2(t) u(t) k T(t)  + m3 lambda3(t)
%      lambda3(tfinal) = 0.0
%    and optimality condition:
%      2 ( u* - 1 ) + ( lambda2* - lambda1* ) k V(t) T(t)
%    and
%      0 <= u(t) <= 1
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 October 2022
%
%  Author:
%
%    Suzanne Lenhart, John Workman
%    Modifications by John Burkardt
%
%  Reference:
%
%    Suzanne Lenhart, John Workman,
%    Optimal Control Applied to Biological Models,
%    Chapman and Hall, 2007
%    LC: QH324.8.L46,
%    ISBN: 978-1-58488-640-2
%
%  Input:
%
%    real s: a parameter in the thymus source term.
%
%    real m1, m2, m3: the natural death rates of uninfected T cells,
%    infected T cells, and free virus particles.
%
%    real r: the logistic growth rate for T cells.
%
%    real tmax: the maximum level of T cells (including infected ones).
%
%    real k: the infection rate coefficient.
%
%    real p: average number of virus particles produced by an infected T cell.
%
%    real T0, Ti0, V0: the initial values of the state variables.
%
%    real A: a weight coefficient in the cost function.
%
%    real T: the final time.
%
%    integer N: the number of time steps.
%
%  Output:
%
%    real t(N+1): the time values.
%
%    real x(N+1,3): the estimated state variables T, Ti, V.
%
%    real u(N+1): the estimated control.
%
%    real lambda(N+13): the estimated adjoint variables.
%
%    real cost(): the cost history.
%
  t0 = 0.0;
  t1 = T;
%
%  Initialize the problem variables.
%
  t = linspace ( t0, t1, N+1 );
  t = t';
  tb = t(end:-1:1);
  x = zeros ( N + 1, 3 );
  x(1,1:3) = [ T0, Ti0, V0 ];
  x0 = [ T0, Ti0, V0 ];
  u = zeros(N+1,1);
  l = zeros(N+1,3);
  l1 = [ 0.0, 0.0, 0.0 ];
  cost = [];
%
%  Forward/backward sweep.
%
  it = 0;

  while ( true )
%
%  Save previous variables for convergence test.
%
    oldu = u;
    oldx = x;
    oldl = l;
%
%  Define the right hand side of the state equation.
%
    dxdt = @(tval,xval) g ( tval, xval, t, x, u, k, m1, m2, m3, p, r, s, tmax );
%
%  Solve state equation forward in time.
%  X will be dimension (N+1,5);
%
    [ t, x ] = ode45 ( dxdt, t, x0 );
%
%  Reverse time for adjoint equation.
%
    xb(1:N+1,1:3) = x(N+1:-1:1,1:3);
    ub = u(end:-1:1);
    lb(1:N+1,1:3) = l(N+1:-1:1,1:3);
%
%  Define the right hand side of the adjoint equation.
%
    dldt = @(tval,lval) h ( tval, lval, tb, xb, ub, lb, A, k, m1, m2, m3, p, r, s, tmax );
%
%  Solve the adjoint equation "in the wrong time direction".
%
    [ tb, lb ] = ode45 ( dldt, tb, l1 );
%
%  Re-reverse adjoint variable.
%
    l = lb(N+1:-1:1,1:3);
%
%  Update the control.
%
    u = 0.5 * ( k * ( l(:,2) - l(:,1) ) .* x(:,3) .* x(:,1) + 2.0 );
    u = max ( u, 0.0 );
    u = min ( u, 1.0 );
    u =  0.5 * ( u + oldu );
%
%  Evaluate the cost.
%
    f = A * x(:,1) - ( 1.0 - u ).^2;
    fint = trapz ( t, f );
    cost = [ cost, fint ];
%
%  Convergence test.
%
    ratio = Inf;
    ratio = convergence ( ratio, u, oldu );
    ratio = convergence ( ratio, x(:,1), oldx(:,1) );
    ratio = convergence ( ratio, x(:,2), oldx(:,2) );
    ratio = convergence ( ratio, x(:,3), oldx(:,3) );
    ratio = convergence ( ratio, l(:,1), oldl(:,1) );
    ratio = convergence ( ratio, l(:,2), oldl(:,2) );
    ratio = convergence ( ratio, l(:,3), oldl(:,3) );
    it = it + 1;
    if ( 0.0 <= ratio )
      fprintf ( 1, 'opt08(): Convergence after %d iterations\n', it );
      break;
    end

  end

  return
end
function gval = g ( tval, xval, t, x, u, k, m1, m2, m3, p, r, s, tmax )

%*****************************************************************************80
%
%% g() evaluates the right hand side of the state equation.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 October 2022
%
%  Author:
%
%    John Burkardt
%
  tcval = xval(1);
  tival = xval(2);
  vval = xval(3);
  uval = interp1 ( t, u, tval );
  gval = zeros ( 3, 1 );

  gval(1) = s / ( 1.0 + vval ) - m1 * tcval ...
    + r * tcval * ( 1.0 - ( tcval - tival ) / tmax ) - k * uval * vval * tcval;
  gval(2) = k * uval * vval * tcval - m2 * tival;
  gval(3) = p * m2 * tival - m3 * vval;

  return
end
function hval = h ( tval, lval, tb, xb, ub, lb, A, k, m1, m2, m3, p, r, s, tmax )

%*****************************************************************************80
%
%% h() evaluates the right hand side of the adjoint equation.
%
%  Modified:
%
%    07 October 2022
%
%  Author:
%
%    John Burkardt
%
  ltcval = lval(1);
  ltival = lval(2);
  lvval = lval(3);
  tcval = interp1 ( tb, xb(:,1), tval );
  tival = interp1 ( tb, xb(:,2), tval );
  vval = interp1 ( tb, xb(:,3), tval );
  uval = interp1 ( tb, ub, tval );
  hval = zeros ( 3, 1 );
 
  hval(1) = - A - ltcval * ( - m1 + r * ( 1.0 - ( tcval - tival ) / tmax ) ...
    - r * tcval / tmax - k * uval * vval ) - k * ltival * uval * vval;
  hval(2) = ltcval * r * tcval / tmax + m2 * ltival - lvval * p * m2;
  hval(3) = - ltcval * ( - s / ( 1.0 + vval )^2 - k * uval * tcval ) ...
    - ltival * k * uval * tcval + m3 * lvval;

  return
end

