function test04 ( )

%*****************************************************************************80
%
%% test04() seeks an optimal control for two cases.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 September 2022
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test04():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Seek optimal controls for two cases.\n' );

  case1 ( );
  case2 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test04():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function case1 ( )

%*****************************************************************************80
%
%% case1() seeks an optimal control for example 1.5.
%
%  Discussion:
%
%    J(t,x,u) = 1/2 int(0<=t<=1) (x(t)-y(t))^2 + alpha u(t)^2 dt
%    alpha = 1
%    y(t) = t + 1
%    x'(t) = u(t), x(0) = 1
%
  t0 = 0.0;
  t1 = 1.0;
  n = 100;
  t = linspace ( t0, t1, n + 1 );
  t = t';
  y = t + 1.0;
  x = zeros ( n + 1, 1 );
  lam = zeros ( n + 1, 1 );
  step_num = 30;
  cost = zeros ( step_num, 1 );
%
%  Using alpha=0.1 makes a huge improvement in the approximation
%  of y(t) by x(t).
%
  alpha = 1.0;

  mu = + 0.1;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Control function update rate mu = %g\n', mu );

  for step = 1 : step_num
%
%  Initialize the control function.
%  Then update it with gradient descent.
%
    if ( 1 == step )
      u = zeros ( n + 1, 1 );
    else
      u = u - mu * ( 2.0 * alpha * u + lam );
    end
%
%  Forward Euler on x.
%
    dt = ( t1 - t0 ) / n;
    for i = 1 : n + 1
      if ( i == 1 )
        x(i) = 1.0;
      else
        x(i) = x(i-1) + dt * u(i-1);
      end
    end
%
%  Evaluate cost.
%
    c = ( x - y ).^2 + alpha * u.^2;
    cost(step) = 0.5 * trapz ( t, c );
    fprintf ( '    Step %d, Cost = %g\n', step, cost(step) );
%
%  Integrate backward in time for LAMBDA.
%
    dt = ( t0 - t1 ) / n;
    for i = n + 1 : -1 : 1
      if ( i == n + 1 )
        lam(i) = 0.0;
      else
        lam(i) = lam(i+1) + dt * ( - 2.0 * ( x(i+1) - y(i+1) ) );
      end
    end
%
%  Display the current control.
%
    figure ( 1 );
    plot ( t, u, 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<-- t -->' );
    ylabel ( '<-- u(t) -->' );
    s = sprintf ( 'Case 1, u, step %d', step );
    title ( s );
    filename = sprintf ( 'test04/case1_u_step%d.png', step );

    print ( '-dpng', filename );

  end

  figure ( 2 );
  hold ( 'on' );
  plot ( t, x, 'b-', 'linewidth', 3 );
  plot ( t, y, 'r-', 'linewidth', 3 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<-- t -->' );
  ylabel ( '<-- X(t) -->' );
  legend ( 'x(t)', 'target y(t)' );
  title ( 'Case 1, Final X(t) and target Y(t)' );
  filename = 'test04/case1_x.png';
  print ( '-dpng', filename );

  figure ( 3 );
  plot ( t, lam, 'linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<-- t -->' );
  ylabel ( '<-- {\lambda}(t) -->' );
  s = sprintf ( 'Case 1, Adjoint {\\lambda}(t)' );
  title ( s );
  filename = sprintf ( 'test04/case1_lam.png' );
  print ( '-dpng', filename );

  figure ( 4 );
  plot ( 1:step_num, cost, 'linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<-- step -->' );
  ylabel ( '<-- Cost -->' );
  title ( 'Case 1, Cost history' );
  filename = sprintf ( 'test04/case1_cost.png' );
  print ( '-dpng', filename );

  return
end

function case2 ( )

%*****************************************************************************80
%
%% case2() seeks an optimal control of a variation of example 1.5.
%
%  Discussion:
%
%    J(t,x,u) = 1/2 int(0<=t<=1) (x(t)-y(t))^2 + alpha u(t)^2 dt
%    alpha = 1.0e-4
%    y(t) = 1+sin(5*t)
%    x'(t) = u(t), x(0) = 1
%
  t0 = 0.0;
  t1 = 1.0;
  n = 100;
  t = linspace ( t0, t1, n + 1 );
  t = t';
  y = 1.0 + sin ( 5.0 * t );
  x = zeros ( n + 1, 1 );
  lam = zeros ( n + 1, 1 );
  step_num = 200;
  cost = zeros ( step_num, 1 );
  alpha = 1.0e-04;

  mu = + 0.1;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Control function update rate mu = %g\n', mu );

  for step = 1 : step_num
%
%  Initialize the control function.
%  Then update it with gradient descent.
%
    if ( 1 == step )
      u = zeros ( n + 1, 1 );
    else
      u = u - mu * ( 2.0 * alpha * u + lam );
    end
%
%  Forward Euler on x.
%
    dt = ( t1 - t0 ) / n;
    for i = 1 : n + 1
      if ( i == 1 )
        x(i) = 1.0;
      else
        x(i) = x(i-1) + dt * u(i-1);
      end
    end
%
%  Evaluate cost.
%
    c = ( x - y ).^2 + alpha * u.^2;
    cost(step) = 0.5 * trapz ( t, c );
    fprintf ( '    Step %d, Cost = %g\n', step, cost(step) );
%
%  Integrate backward in time for LAMBDA.
%
    dt = ( t0 - t1 ) / n;
    for i = n + 1 : -1 : 1
      if ( i == n + 1 )
        lam(i) = 0.0;
      else
        lam(i) = lam(i+1) + dt * ( - 2.0 * ( x(i+1) - y(i+1) ) );
      end
    end
%
%  Display the current control.
%
    figure ( 5 );
    plot ( t, u, 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<-- t -->' );
    ylabel ( '<-- u(t) -->' );
    s = sprintf ( 'Case 2, u, step %d', step );
    title ( s );
    filename = sprintf ( 'test04/case2_u_step%d.png', step );

    print ( '-dpng', filename );

  end

  figure ( 6 );
  hold ( 'on' );
  plot ( t, x, 'b-', 'linewidth', 3 );
  plot ( t, y, 'r-', 'linewidth', 3 );
  hold ( 'off' );
  grid ( 'on' );
  xlabel ( '<-- t -->' );
  ylabel ( '<-- X(t) -->' );
  legend ( 'x(t)', 'target y(t)' );
  title ( 'Case 2, Final X(t) and target Y(t)' );
  filename = 'test04/case2_x.png';
  print ( '-dpng', filename );

  figure ( 7 );
  plot ( t, lam, 'linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<-- t -->' );
  ylabel ( '<-- {\lambda}(t) -->' );
  s = sprintf ( 'Case 2, Adjoint {\\lambda}(t)' );
  title ( s );
  filename = sprintf ( 'test04/case2_lam.png' );
  print ( '-dpng', filename );

  figure ( 8 );
  plot ( 1:step_num, cost, 'linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<-- step -->' );
  ylabel ( '<-- Cost -->' );
  title ( 'Case 2, Cost history' );
  filename = sprintf ( 'test04/case2_cost.png' );
  print ( '-dpng', filename );

  return
end

