function [ t, x, u, l, cost ] = opt07 ( a, b, c, d, e, q, S0, E0, I0, R0, ...
  P0, A, T, N )

%*****************************************************************************80
%
%% code07() performs optimal control for case #07.
%
%  Discussion:
%
%    Epidemic Model
%
%    minimize
%      integral ( 0 <= t <= T ) A I(t) + u(t)^2 dt
%    subject to
%      S'(t) = b P(t) - d S(t) - c S(t) I(t) - u(t) S(t)
%      E'(t) = c S(t) I(t) - ( e + d ) E(t)
%      I'(t) = e E(t) - ( q + a + d ) * I(t)
%      R'(t) = q I(t) - d R(t) + u(t) S(t)
%      P'(t) = ( b - d ) P(t) - a I(t)
%      S(0) = S0
%      E(0) = E0
%      I(0) = I0
%      R(0) = R0
%      P(0) = P0 = S0 + E0 + I0 + R0
%    with Hamiltonian:
%      H = A I(t) + u(t)^2 
%         + lambdaS(t) ( b P(t) - d S(t) - c S(t) I(t) - u(t) S(t) )
%         + lambdaE(t) ( c S(t) I(t) - ( e + d ) E(t) )
%         + lambdaI(t) ( e E(t) - ( q + a + d ) * I(t) )
%         + lambdaR(t) ( q I(t) - d R(t) + u(t) S(t) )
%         + lambdaP(t) ( ( b - d ) P(t) - a I(t) )
%    with adjoint equations:
%      lambdaS(t)' = lambdaS(t) * ( d + c I(t) + u(t) ) - lambdaE(t) c I(t) 
%        - lambdaR(t) * u(t)
%      lambdaS(T) = 0.0
%      lambdaE(t)' = lambdaE(t) ( e + d ) - lambdaI(t) e
%      lambdaE(T) = 0.0
%      lambdaI(t)' = - A + lambdaS(t) c S(t) - lambdaE(t) c S(t) + lambdaI(t) ( q + a + d )
%        - lambdaR(t) q + lambdaP(t) a
%      lambdaI(T) = 0.0
%      lambdaR(t)' = lambdaR(t) d
%      lambdaR(T) = 0.0
%      lambdaP(t)' = - lambdaS(t) b - lambdaP(t) ( b - d ) )
%      lambdaP(T) = 0.0
%    and optimality condition:
%      2 u(t) - lambdaS(t) S(t) + lambdaR(t) S(t) = 0
%    and
%      0 <= u(t) <= 0.9
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 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 a: the death rate due to disease.
%
%    real b: the natural birth rate.
%
%    real c: the disease transmission rate.
%
%    real d: the natural death rate.
%
%    real e: the rate at which exposed individuals become infectious.
%
%    real q: the rate at which infectious individuals recover.
%
%    real S0, E0, I0, R0, P0: the initial values of the state variables.
%    It should be the case that P0 = S0 + E0 + I0 + R0;
%
%    real A: a weight for 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,5): the state estimates S, E, I, R, P.
%
%    real u(N+1): the optimal control estimate.
%
%    real l(N+1,5): the adjoint estimates.
%
%    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,5);
  x(1,1:5) = [ S0, E0, I0, R0, P0 ];
  x0 = [ S0, E0, I0, R0, P0 ];
  u=zeros(N+1,1);
  l=zeros(N+1,5);
  l1 = zeros(1,5);
  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, a, b, c, d, e, q );
%
%  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:5) = x(N+1:-1:1,1:5);
    ub = u(end:-1:1);
    lb(1:N+1,1:5) = l(N+1:-1:1,1:5);
%
%  Define the right hand side of the adjoint equation.
%
    dldt = @(tval,lval) h ( tval, lval, tb, xb, ub, lb, A, a, b, c, d, e, q );
%
%  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:5);
%
%  Update the control.
%
    u = 0.5 * x(:,1) .* ( l(:,1) - l(:,4) );
    u = max ( u, 0.0 );
    u = min ( u, 0.9 );
    u = 0.5 * ( u + oldu );
%
%  Evaluate the cost.
%
    f = A * x(:,3) + 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, x(:,4), oldx(:,4) );
    ratio = convergence ( ratio, x(:,5), oldx(:,5) );
    ratio = convergence ( ratio, l(:,1), oldl(:,1) );
    ratio = convergence ( ratio, l(:,2), oldl(:,2) );
    ratio = convergence ( ratio, l(:,3), oldl(:,3) );
    ratio = convergence ( ratio, l(:,4), oldl(:,4) );
    ratio = convergence ( ratio, l(:,5), oldl(:,5) );
    it = it + 1;
    if ( 0.0 <= ratio )
      fprintf ( 1, 'opt07(): Convergence after %d iterations\n', it );
      break;
    end

  end

  return
end
function gval = g ( tval, xval, t, x, u, a, b, c, d, e, q )

%*****************************************************************************80
%
%% g() evaluates the right hand side of the state equation.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 October 2022
%
%  Author:
%
%    John Burkardt
%
  sval = xval(1);
  eval = xval(2);
  ival = xval(3);
  rval = xval(4);
  pval = xval(5);
  uval = interp1 ( t, u, tval );
  gval = zeros ( 5, 1 );

  gval(1) = b * pval - d * sval - c * sval * ival - uval * sval;
  gval(2) = c * sval * ival - ( e + d ) * eval;
  gval(3) = e * eval - ( q + a + d ) * ival;
  gval(4) = q * ival - d * rval + uval * sval;
  gval(5) = ( b - d ) * pval - a * ival;

  return
end
function hval = h ( tval, lval, tb, xb, ub, lb, A, a, b, c, d, e, q )

%*****************************************************************************80
%
%% h() evaluates the right hand side of the adjoint equation.
%
%  Modified:
%
%    06 October 2022
%
%  Author:
%
%    John Burkardt
%
  lsval = lval(1);
  leval = lval(2);
  lival = lval(3);
  lrval = lval(4);
  lpval = lval(5);
  sval = interp1 ( tb, xb(:,1), tval );
  eval = interp1 ( tb, xb(:,2), tval );
  ival = interp1 ( tb, xb(:,3), tval );
  rval = interp1 ( tb, xb(:,4), tval );
  pval = interp1 ( tb, xb(:,5), tval );
  uval = interp1 ( tb, ub, tval );
  hval = zeros ( 5, 1 );
 
  hval(1) = lsval * ( d + c * ival + uval ) - leval * c * ival - lrval * uval;
  hval(2) = leval * ( e + d ) - lival * e;
  hval(3) = - A + lsval * c * sval - leval * c * sval ...
    + lival * ( q + a + d ) - lrval * q + lpval * a;
  hval(4) = lrval * d;
  hval(5) = - lsval * b - lpval * ( b - d );

  return
end

