Repeated Computation of Control Problem Gradient Using Newton Method

repeat_newton_gradient

@(@\newcommand{\B}[1]{{\bf #1}} \newcommand{\R}[1]{{\rm #1}}@)@Repeated Computation of Control Problem Gradient Using Newton Method
Syntax

repeat_newton_gradient(repeat, J, reverse, x, grad)

repeat_forward_gradient(repeat, size)

repeat_reverse_gradient(repeat, size)

repeat
is the number of times to choose a new value for the controls.

J
is the number of time points in the discrete version of the control problem.

reverse
If this is true, reverse mode is used to compute the gradient. Otherwise, forward mode is used.

x
If x has size zero, it is not used. Otherwise, it must have size 2*J . In this case it specifies the value of the control vector for the first gradient computation.

grad
This vector must have size 2*J . The input value of its elements does not matter. Upon return, it contains the value last gradient computed by repeat_kedem_gradient.

size
This is the number time points in the problem; i.e., the value of J for this test.

void repeat_newton_gradient(
     size_t                repeat    ,
     size_t                J         ,
     bool                  reverse   ,
     const VECTOR(double)& x         ,
     VECTOR(double)&       grad      )
{    assert( x.size() == 0 || x.size() == grad.size() );
     assert( size_t( grad.size() ) == 2 * J );
     double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
     typedef CppAD::AD<double> adouble;
     // ----------------------------------------------------------------------
     // control problem parameters
     double T;
     VECTOR(double) p(4), q(4);
     set_T_p_and_q(T, p, q);
     double delta_t = T / double(J - 1);
     // ----------------------------------------------------------------------
     // L_fun, F_fun
     CppAD::ADFun<double> F_fun, L_fun;
     control::rec_objective(F_fun, J,  delta_t, q);
     control::rec_constraint(L_fun, J, delta_t, p);
     size_t m = L_fun.Range();
     size_t n = L_fun.Domain() - m;
     assert( F_fun.Range() == 1 );
     assert( F_fun.Domain() == n + m );
     // -----------------------------------------------------------------------
     // aL_fun
     CppAD::ADFun<adouble> aL_fun;
     adouble adelta_t = delta_t;
     VECTOR(adouble) ap(4);
     for(size_t i = 0; i < 4; i++)
          ap[i] = p[i];
     control::rec_constraint(aL_fun, J, adelta_t, ap);
     // -----------------------------------------------------------------------
     // control_solve
     double criteria = double(n + m) * eps99;
     control::implicit_solver control_solve(L_fun, aL_fun, criteria);
     // -----------------------------------------------------------------------
     // R_newton
     size_t num_step  = 1;
     bool   full_step = false; // does not matter when num_step = 1
     implicit_newton<control::implicit_solver> R_newton(
          full_step, num_step, aL_fun, F_fun, control_solve
     );
     // x0, w
     VECTOR(double) x0(n), w(1), x1(n), dR(1);
     w[0] = 1.0;
     for(size_t j = 0; j < n; j++)
          x1[j] = 0.0;
     if( x.size() == 0 )
          CppAD::uniform_01(n, x0);
     else
          x0 = x;
     for(size_t count = 0; count < repeat; count++)
     {
          //
          // zero order forward
          R_newton.Forward(0, x0);
          //
          if( reverse )
               grad = R_newton.Reverse(1, w);
          else
          {    // use first order forward to calculate gradient
               for(size_t j = 0; j < n; j++)
               {    x1[j] = 1.0;
                    dR    = R_newton.Forward(1, x1);
                    grad[j] = dR[0];
                    x1[j] = 0.0;
               }
          }
          //
          // pick a random value for next x0
          CppAD::uniform_01(n, x0);
     }
}
void time_forward_gradient(size_t size, size_t repeat)
{    size_t J       = size;
     bool   reverse = false;
     VECTOR(double) grad(2 * J), x(0);
     repeat_newton_gradient(repeat, J, reverse, x, grad);
}
void time_reverse_gradient(size_t size, size_t repeat)
{    size_t J       = size;
     bool   reverse = true;
     VECTOR(double) grad(2 * J), x(0);
     repeat_newton_gradient(repeat, J, reverse, x, grad);
}

Input File: src/time.cpp