function [ t, x, u, l, cost ] = opt09 ( r, K, mp, mf, P0, F0, O0, cp, cf, T, N )

%*****************************************************************************80
%
%% opt09() performs optimal control for case #09.
%
%  Discussion:
%
%    Bear Populations
%
%    minimize
%      integral ( 0 <= t <= T ) O(t) + cp up(t)^2 + cf uf(t)^2 dt
%    subject to
%      P'(t) = r P(t) - r / K P(t)^2 + mf r / K ( 1 - P(t) / K ) F(t)^2 
%        - up(t) P(t)
%      F'(t) = r F(t) - r / K F(t)^2 + mp r / K ( 1 - F(t) / K ) P(t)^2
%        - uf(t) F(t)
%      O'(t) = r ( 1 - mp ) P(t)^2 / K 
%            + r ( 1 - mf ) F(t)^2 / K
%            + mf r / K^2 P(t) F(t)^2 
%            + mp r / K^2 F(t) P(t)^2
%      P(0) = P0
%      F(0) = F0
%      O(0) = O0
%    with Hamiltonian:
%      H = O(t) + cp up(t)^2 + cf uf(t)^2 
%        + lambda1(t) ( r P(t) - r / K P(t)^2 + mf r / K ( 1 - P(t) / K ) F(t)^2
%          - up(t) P(t) )
%        + lambda2(t) ( r F(t) - r / K F(t)^2 + mp r / K ( 1 - F(t) / K ) P(t)^2
%          - uf(t) F(t) )
%        + lambda3(t) ( r ( 1 - mp ) P(t)^2 / K + r ( 1 - mf ) F(t)^2 / K
%          + mf r / K^2 P(t) F(t)^2 )
%    with adjoint equation:
%      lambda1(t)' = - lambda1(t) ( r - 2 r / K P(t) - mf r / K^2 F(t)^2 - up(t) )
%                    - lambda2(t) ( 2 mp r / K ( 1 - F(t) / K ) P(t) )
%                    - lambda3(t) ( 2 * r ( 1 - mp ) P(t) / K 
%                                   + mf r / K^2 F(t)^2 
%                                   + 2 mp r / K^2 P(t) F(t) )
%      lambda1(T) = 0.0
%      lambda2(t)' = - lambda1(t) 2 mf r / K ( 1 - P(t) / K ) F(t) 
%                    - lambda2(t) ( r - 2 r / K F(t) - mp r / K^2 P(t)^2 - uf(t) ) )
%                    - lambda3(t) ( 2 r ( 1 - mf ) F(t) / K 
%                                   + mp r / K^2 P(t)^2
%                                   + 2 mf r / K^2 P(t) F(t) )
%      lambda2(T) = 0.0
%      lambda3(t)' = - 1
%      lambda3(T) = 0.0
%    and optimality conditions:
%      2 cp up(t) - lambda1(t) P(t) = 0
%      2 cf uf(t) - lambda2(t) F(t) = 0
%    and
%      0 <= up(t) <= 1
%      0 <= uf(t) <= 1
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 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
%
%  Output:
%
%    real r: the population growth rate of the bears.
%
%    real K: the carrying capacity.
%
%    real mp: the proportion of the park boundary connected to the forest.
%
%    real mf: the proportion of the forest boundary connected to the park.
%
%    real P0, F0, O0: the initial values for the state variables.
%
%    real cp, cf: weight parameters 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 P, F, O state variable values.
%
%    real u(N+1,2): the control values.
%
%    real l(N+1,3): the 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) = [ P0, F0, O0 ];
  x0 = [ P0, F0, O0 ];
  u = zeros ( N+1, 2 );
  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, mf, mp, r );
%
%  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(1:N+1,1:2) = u(N+1:-1:1,1:2);
    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, K, mf, mp, r );
%
%  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 controls.
%
    u(:,1) = 0.5 * l(:,1) .* x(:,1) / cp;
    u(:,1) = min ( u(:,1), 1.0 );
    u(:,1) = max ( u(:,1), 0.0 );
    u(:,1) = 0.5 * ( u(:,1) + oldu(:,1) );

    u(:,2) = 0.5 * l(:,2) .* x(:,2) / cf;
    u(:,2) = min ( u(:,2), 1.0 );
    u(:,2) = max ( u(:,2), 0.0 );
    u(:,2) = 0.5 * ( u(:,2) + oldu(:,2) );
%
%  Evaluate the cost.
%
    f = x(:,3) + cp * u(:,1).^2 + cf * u(:,2).^2;
    fint = trapz ( t, f );
    cost = [ cost, fint ];
%
%  Convergence test.
%
    ratio = Inf;
    ratio = convergence ( ratio, u(:,1), oldu(:,1) );
    ratio = convergence ( ratio, u(:,2), oldu(:,2) );
    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, 'opt09(): Convergence after %d iterations\n', it );
      break;
    end

  end

  return
end
function gval = g ( tval, xval, t, x, u, K, mf, mp, r )

%*****************************************************************************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
%
  pval = xval(1);
  fval = xval(2);
  oval = xval(3);
  upval = interp1 ( t, u(:,1), tval );
  ufval = interp1 ( t, u(:,2), tval );
  gval = zeros ( 3, 1 );

  gval(1) = pval * ( r - r / K * pval - upval ) ...
    + mf / K * r * ( 1.0 - pval / K ) * fval^2;
  gval(2) = fval * ( r - r / K * fval - ufval ) ...
    + mp / K * r * ( 1.0 - fval / K ) * pval^2;
  gval(3) = ( 1.0 - mp ) * r / K * pval^2 ...
          + ( 1.0 - mf ) * r / K * fval^2 ...
          + mf * r / K^2 * pval * fval^2 ...
          + mp * r / K^2 * fval * pval^2;

  return
end
function hval = h ( tval, lval, tb, xb, ub, lb, K, mf, mp, r )

%*****************************************************************************80
%
%% h() evaluates the right hand side of the adjoint equation.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 October 2022
%
%  Author:
%
%    John Burkardt
%
  lpval = lval(1);
  lfval = lval(2);
  loval = lval(3);
  pval = interp1 ( tb, xb(:,1), tval );
  fval = interp1 ( tb, xb(:,2), tval );
  oval = interp1 ( tb, xb(:,3), tval );
  upval = interp1 ( tb, ub(:,1), tval );
  ufval = interp1 ( tb, ub(:,2), tval );
  hval = zeros ( 3, 1 );
 
  hval(1) = - lpval * ( r - 2.0 * r / K * pval ...
                      - mf * r / K^2 * fval^2 - upval ) ...
            - lfval * ( 2.0 * mp * r / K * ( 1.0 - fval / K ) * pval ) ...
            - loval * ( 2.0 * r / K * ( 1.0 - mp ) * pval ...
                        + mf * r / K^2 * fval^2 ...
                        + 2.0 * mp * r / K^2 * pval * fval );

  hval(2) = - lpval * ( 2.0 * mf * r / K * ( 1.0 - pval / K ) * fval ) ...
            - lfval * ( r - 2.0 * r / K * fval ...
                      - mp * r / K^2 * pval^2 - ufval ) ...
            - loval * ( 2.0 * r / K * ( 1.0 - mf ) * fval ...
                        + mp * r / K^2 * pval^2 ...
                        + 2.0 * mf * r / K^2 * pval * fval );
  hval(3) = - 1.0;

  return
end

