function [ t, x, u, l, cost ] = opt10 ( a, b, c, x10, A, D, T, N )

%*****************************************************************************80
%
%% opt10() performs optimal control for case #10.
%
%  Discussion:
%
%    Glucose Model
%
%    minimize
%      integral ( 0 <= t <= T ) A ( x1(t) - D )^2 + u(t)^2 dt
%    subject to
%      x1'(t) = - a x1(t) - b x2(t)
%      x2'(t) = - c x2(t) + u(t)
%      x1(0) = x10
%      x2(0) = 0
%    with Hamiltonian:
%      H = A ( x1(t) - D )^2 + u(t)^2 
%        + lambda1(t) ( - a x1(t) - b x2(t) )
%        + lambda2(t) ( - c x2(t) + u(t) )
%    with adjoint equation:
%      lambda1(t)' = - 2 A ( x1(t) - D ) + a lambda1(t) )
%      lambda1(T) = 0.0
%      lambda2(t)' = b lambda1(t) + c lambda2(t) )
%      lambda2(T) = 0.0
%    and optimality condition:
%      2 u* + lambda2*(t) = 0
%
%  Modified:
%
%    09 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, b, c: coefficients in the state equation.
%
%    real x10: the initial condition for x1.
%
%    real A: a weight coefficient in the cost function.
%
%    real D: the desired glucose level.
%
%    real T: the final time.
%
%    integer N: the number of time steps.
%
%  Output:
%
%    real t(N+1): the time values.
%
%    real x(N+1,2): the final state estimates.
%
%    real u(N+1): the final optimal control estimate.
%
%    real l(N+1,2): the final 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,2);
  x0 = [ x10, 0.0 ]
  x(1,1:2) = x0;
  u = zeros(N+1,1);
  l = zeros(N+1,2);
  l1 = [ 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, a, b, c );
%
%  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:2) = x(N+1:-1:1,1:2);
    ub(1:N+1) = u(N+1:-1:1);
    lb(1:N+1,1:2) = l(N+1:-1:1,1:2);
%
%  Define the right hand side of the adjoint equation.
%
    dldt = @(tval,lval) h ( tval, lval, tb, xb, ub, lb, A, a, b, c, D );
%
%  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:2);
%
%  Update the control.
%
    u = - 0.5 * l(:,2);
    u = 0.5 * ( u + oldu );
%
%  Evaluate the cost.
%
    f = A * ( x(:,1) - D ).^2 + 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, l(:,1), oldl(:,1) );
    ratio = convergence ( ratio, l(:,2), oldl(:,2) );
    it = it + 1;
    if ( 0.0 <= ratio )
      fprintf ( 1, 'opt10(): Convergence after %d iterations\n', it );
      break;
    end

  end

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

%*****************************************************************************80
%
%% g() evaluates the right hand side of the state equation.
%
%  Modified:
%
%    07 October 2022
%
%  Author:
%
%    John Burkardt
%
  x1val = xval(1);
  x2val = xval(2);
  uval = interp1 ( t, u, tval );
  gval = zeros ( 2, 1 );

  gval(1) = - a * x1val - b * x2val;
  gval(2) = - c * x2val + uval;

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

%*****************************************************************************80
%
%% h() evaluates the right hand side of the adjoint equation.
%
%  Modified:
%
%    09 October 2022
%
%  Author:
%
%    John Burkardt
%
  l1val = lval(1);
  l2val = lval(2);
  x1val = interp1 ( tb, xb(:,1), tval );
  x2val = interp1 ( tb, xb(:,2), tval );
  uval = interp1 ( tb, ub, tval );
  hval = zeros ( 2, 1 );
 
  hval(1) = - 2.0 * A * ( x1val - D ) + a * l1val;

  hval(2) = b * l1val + c * l2val;

  return
end

