Control Problem Solver for Implicit Kedem or Newton Object

implicit_solver

@(@\newcommand{\B}[1]{{\bf #1}} \newcommand{\R}[1]{{\rm #1}}@)@Control Problem Solver for Implicit Kedem or Newton Object
Syntax

control::implicit_solver solve(L_fun, aL_fun, acriteria)

y = solve.function(x)

b = solve.derivative(x, y)

u = solve.linear(v)

Purpose
The object solve can be used with implicit_kedem to compute derivatives of the function @(@ Y(x) @)@ defined by @(@ L[x, Y(x)] = 0 @)@ for the control problem.

L_fun
This argument has prototype



     const CppAD::ADFun<double>& L_fun

and is the CppAD function object corresponding to the constraint function.

aL_fun
This argument has prototype



     const CppAD::ADFun< CppAD::AD<double> >& aL_fun

and is the CppAD function object corresponding to the constraint function. This function object can be empty, if solve.derivative is note used with Scalar equal to CppAD::AD<double>. An empty aL_fun object can be created with



     CppAD::ADFun< CppAD::AD<double> > aL_fun;

criteria
This is the convergence criteria for @(@ L[x, Y(x)] = 0 @)@. To be specific, a value @(@ y @)@ is accepted if the Euclidean norm squared @(@ | L(x, y) |^2 @)@ is less than criteria . This argument has prototype



     const Solve& solve

The type Solver must support the default constructor and the assignment operator. It must also support the following operations:

solve.function
The argument x and the return value y have prototypes



     const VECTOR(double)&  x

     VECTOR(double)         y

where x.size() == n and y.size() == m . The return value satisfies the relation @(@ L(x, y) = 0 @)@.

solve.derivative
The arguments have prototypes



     const VECTOR(Scalar)&  x

     const VECTOR(Scalar)&  y

The return value has prototype



     VECTOR(Scalar) b

This returns the value of @(@ L_y (x, y) @)@ for subsequent calls to solve.linear . Only the elements of @(@ L_y (x, y) @)@ that depend on @(@ (x, y) @)@ need be included in the vector @(@ b @)@. The type Scalar is either double or CppAD::AD<double>.

solve.linear
The arguments b , v and the return value u have prototypes



     const VECTOR(Scalar)& b

     const VECTOR(Scalar)& v

     VECTOR(Scalar)        u

where both vectors have size m . The return value satisfies the equation @[@ u = L_y (x, y)^{-1} v @]@ where @(@ L_y (x, y) @)@ corresponds to b . The type Scalar is either double or CppAD::AD<double>.

Example


bool test_control_implicit_solver(void)
{    bool ok      = true;
     double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
     //
     // record L_fun
     CppAD::ADFun<double>       L_fun;
     size_t                     J = 5;
     double                     delta_t = 0.1;
     VECTOR(double) p(4);
     for(size_t i = 0; i < 4; i++)
          p[i] = 0.2;
     control::rec_constraint(L_fun, J, delta_t, p);
     //
     size_t m      = L_fun.Range();
     size_t n      = L_fun.Domain() - m;
     //
     // create implicit solver
     double criteria  = double(n + m) * eps99;
     CppAD::ADFun< CppAD::AD<double> > not_used;
     control::implicit_solver solver(L_fun, not_used, criteria);
     //
     // set x
     VECTOR(double) x(n);
     for(size_t i = 0; i < n; i++)
          x[i] = 0.0;
     //
     // solve L(x, y) == 0
     VECTOR(double) y = solver.function(x);
     ok &= size_t( y.size() ) == m;
     //
     // check L(x, y) is near zero
     VECTOR(double) xy(n + m);
     join_vector(xy, x, y);
     VECTOR(double) L  = L_fun.Forward(0, xy);
     ok &= norm_squared(L) < criteria;
     //
     // solve L_y (x, y) * u = v
     VECTOR(double) v(m);
     for(size_t i = 0; i < m; i++)
          v[i] = double(i - 5.0);
     VECTOR(double) b = solver.derivative(x, y);
     VECTOR(double) u = solver.linear(b, v);
     //
     // compute L_y (x, y)
     CPPAD_SPARSE(double) L_y;
     CppAD::sparse_jac_work work;
     jac_constraint(L_y, L_fun, xy, work);
     SparseMatrix<double> L_y_eigen;
     sparse_cppad2eigen(L_y, L_y_eigen);
     //
     // check L_y (x, y) * u = v
     VECTOR(double) d = L_y_eigen * u - v;
     ok &= norm_squared(d) < criteria;
     //
     return ok;
}

Input File: src/control.hpp