function [ t, x, u, l, cost ] = opt03 ( r, A, B, C, x0, T, N )

%*****************************************************************************80
%
%% opt03() performs optimal control for case #03.
%
%  Discussion:
%
%    Bacteria
%
%    maximize
%      C x(1) - integral ( 0 <= t <= 1 ) u(t)^2 dt
%    subject to
%      x'(t) = r x(t) + A u(t) x(t) - B u(t)^2 exp(-x(t))
%      x(0) = x0
%    with Hamiltonian:
%      H = u(t)^2 + lambda(t) * (  r x(t) + A u(t) x(t) - B u(t)^2 exp(-x(t)) )
%    with adjoint equation:
%      lambda'(t) = - lambda(t) * ( r + A u(t) + B u(t)^2 exp(-x(t)) )
%      lambda(1) = C
%    and optimality condition:
%      - 2 u* + A lambda* x* + 2 B lambda* u* exp(-x*) = 0
%
%  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 r: the bacterial growth rate;
%
%    real A, B: two coefficients in the state equation.
%
%    real C: a coefficient in the cost function.
%
%    real x0: the initial value for the state variable.
%
%    real T: the final time.
%
%    integer N: the number of time steps.
%
%  Output:
%
%    real t(N+1): the time values.
%
%    real x(N+1): the state estimate.
%
%    real u(N+1): the optimal control estimate.
%
%    real l(N+1): the adjoint estimate.
%
%    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,1);
  x(1) = x0;
  u = zeros(N+1,1);  
  l = zeros(N+1,1);
  l1 = C;
  l(N+1) = C;
  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, r, A, B );
%
%  Solve state equation forward in time.
%
    [ t, x ] = ode45 ( dxdt, t, x0 );
%
%  Reverse time for adjoint equation.
%
    xb = x(end:-1:1);
    ub = u(end:-1:1);
    lb = l(end:-1:1);
%
%  Define the right hand side of the adjoint equation.
%
    dldt = @(tval,lval) h ( tval, lval, tb, xb, ub, lb, r, A, B );
%
%  Solve the adjoint equation "in the wrong time direction".
%
    [ tb, lb ] = ode45 ( dldt, tb, l1 );
%
%  Re-reverse adjoint variable.
%
    l = lb(end:-1:1);
%
%  Update the control.
%
    u = 0.5 * ( A * l .* x ) ./ ( 1.0 + B * l .* exp ( - x ) );
    u = 0.5 * ( u + oldu );
%
%  Update the cost.
%
    f = u.^2;
    fint = C * x(N+1) - trapz ( t, f );
    cost = [ cost, fint ];
%
%  Convergence test.
%
    ratio = Inf;
    ratio = convergence ( ratio, u, oldu );
    ratio = convergence ( ratio, x, oldx );
    ratio = convergence ( ratio, l, oldl );
    it = it + 1;
    if ( 0.0 <= ratio )
      fprintf ( 1, 'opt03(): Convergence after %d iterations\n', it );
      break;
    end

  end

  return
end
function gval = g ( tval, xval, t, x, u, r, A, B )

%*****************************************************************************80
%
%% g() evaluates the right hand side of the state equation.
%
%  Modified:
%
%    06 October 2022
%
%  Author:
%
%    John Burkardt
%
  uval = interp1 ( t, u, tval );
  gval =  r * xval + A * uval .* xval - B * uval.^2 .* exp ( - xval );

  return
end
function hval = h ( tval, lval, tb, xb, ub, lb, r, A, B )

%*****************************************************************************80
%
%% h() evaluates the right hand side of the adjoint equation.
%
%  Modified:
%
%    06 October 2022
%
%  Author:
%
%    John Burkardt
%
  xval = interp1 ( tb, xb, tval );
  uval = interp1 ( tb, ub, tval );
  hval = - lval * ( r + A * uval + B * uval.^2 .* exp ( - xval ) );

  return
end

