Compute Jacobian of Implicit Function Constraints

jac_constraint

@(@\newcommand{\B}[1]{{\bf #1}} \newcommand{\R}[1]{{\rm #1}}@)@Compute Jacobian of Implicit Function Constraints
Syntax
jac_constraint(L_y, L_fun, xy, work)

Prototype


template <class Scalar>
void jac_constraint(
     CPPAD_SPARSE(Scalar)&                             L_y         ,
     CppAD::ADFun<Scalar>&                             L_fun       ,
     VECTOR(Scalar)&                                   xy          ,
     CppAD::sparse_jac_work&                           work        )

Purpose
We are given a function @(@ L : \B{R}^n \times \B{R}^m \rightarrow \B{R}^m @)@ and define the implicit function @(@ Y : \B{R}^n \rightarrow \B{R}^m @)@ by the constraint equation @[@ L[ x , Y(x) ] = 0 @]@ This routine computes the sparse derivative of @(@ L (x , y) @)@ w.r.t @(@ y @)@.

L_y
This argument must initially be created with the empty constructor; i.e.,



     CPPAD_SPARSE(Scalar) L_y

It's return value is a sparse matrix representation of the Jacobian of the constraints with respect to the implicit variables; i.e., the partial of @(@ L(x, y) @)@ w.r.t @(@ y @)@. The return value can be used to speed up subsequent calls when the operation sequence in L_fun does not change. Note that this is an @(@ m \times m @)@ matrix; i.e., the column indices are relative to @(@ y @)@ and not for @(@ (x, y) @)@. The entries in L_y are in row-major order.

L_fun
The operation sequence for the constraint function @(@ L : \B{R}^n \times \B{R}^m \rightarrow \B{R}^m @)@ is stored in L_fun .

xy
This is the argument @(@ (x, y) @)@ at which we are evaluating the Jacobian. It size is n+m .

work
This is a work vector used to reduce the work. If the operation sequence ion L_fun changes, L_y and work should be reset to empty. Otherwise, they should remain unchanged from their previous return value.

Example

Simple


bool test_jac_constraint(void)
{    bool ok = true;
     double eps99 = 99.0 * std::numeric_limits<double>::epsilon();

     // record z = L(x, y)
     VECTOR( CppAD::AD<double> ) axy(2), az(1);
     axy[0] = 0.5;
     axy[1] = 0.5;
     CppAD::Independent( axy );
     az[0] = axy[0] * axy[0] + axy[1] * axy[1] - 1.0;
     CppAD::ADFun<double> L_fun(axy, az);

     // choose a point that satisfies the constraints
     VECTOR(double) xy(2);
     xy[0] = 0.5;
     xy[1] = std::sqrt( 1.0 - xy[0] * xy[0] );

     // Evaluate L_y (x, y)
     CPPAD_SPARSE(double) L_y;
     CppAD::sparse_jac_work work;
     jac_constraint(L_y, L_fun, xy, work);

     // check the result is L_y(x, y) = 2.0 * y
     ok &= L_y.nr()  == L_fun.Range();
     ok &= L_y.nc()  == L_fun.Range();
     ok &= L_y.nnz() == 1;
     ok &= L_y.row()[0] == 0;
     ok &= L_y.col()[0] == 0;
     ok &= CppAD::NearEqual( L_y.val()[0], 2.0 * xy[1], eps99, eps99 );
     //
     return ok;
}

Control Problem


# include "control.hpp"
bool test_control_jac_constraint(void)
{    bool ok = true;
     double eps99 = 99.0 * std::numeric_limits<double>::epsilon();
     //
     // record L_fun
     CppAD::ADFun<double>       L_fun;
     size_t                     J = 6;
     double                     delta_t = 2.0;
     VECTOR(double) p(4);
     for(size_t i = 0; i < 4; i++)
          p[i] = double(5 - i);
     control::rec_constraint(L_fun, J, delta_t, p);
     size_t m      = L_fun.Range();
     size_t n      = L_fun.Domain() - m;
     //
     // compute L_y (x, y)
     CPPAD_SPARSE(double) L_y;
     VECTOR(double) xy(n + m);
     for(size_t i = 0; i < n+m; i++)
          xy[i] = 0.5;
     CppAD::sparse_jac_work work;
     //
     // all vj are 0.5
     double vj  = 0.5;
     double r   = sqrt( (vj + 1.0)*(vj + 1.0) + vj * vj );
     double r3  = r * r * r;
     double r5  = r * r * r * r * r;
     double a   = 1.0 / r3 - 1.0;
     double ap0 = -3.0 * (vj + 1.0) / r5;
     double ap1 = -3.0 * vj / r5;
     //
     // check repeated calls to jac_constraint
     for(size_t count = 0; count < 2; count ++)
     {    jac_constraint(L_y, L_fun, xy, work);
          //
          // check Jacobian
          size_t nnz = L_y.nnz();
          const VECTOR(size_t)& row( L_y.row() );
          const VECTOR(size_t)& col( L_y.col() );
          const VECTOR(double)& val( L_y.val() );
          for(size_t k = 0; k < nnz; k++)
          {    //
               size_t i = row[k];
               size_t j = col[k];
               double v = val[k];
               //
               // check lower triangular
               ok &= j <= i;
               //
               bool diagonal_block = (i / 4) == (j / 4);
               if( diagonal_block )
               {    ok &= i == j;
                    ok &= v == 1.0;
               }
               else
               {    double check;
                    // must be the lower diagonal block
                    ok &= (i / 4) == (j / 4) + 1;
                    //
                    // derivative of this v w.r.t previous
                    size_t i_4 = i % 4;
                    size_t j_4 = j % 4;
                    if( i_4 == 0 )
                    {    // derivative of this v0 w.r.t previous v0 or v2
                         ok &= j_4 == 0 || j_4 == 2;
                         if( j_4 == i_4 )
                              ok &= v == -1.0;
                         else
                              ok &= v == - 1.0 * delta_t;
                    }
                    if( i_4 == 1 )
                    {    // derivative of this v1 w.r.t previous v1 or v3
                         ok &= j_4 == 1 || j_4 == 3;
                         if( j_4 == i_4 )
                              ok &= v == -1.0;
                         else
                              ok &= v == - 1.0 * delta_t;
                    }
                    if( i_4 == 2 )
                    {    // derivative of this v2
                         switch( j_4 )
                         {
                              case 0:
                              // w.r.t previous v0
                              check = (a + (1.0 + vj) * ap0) * delta_t;
                              ok &= CppAD::NearEqual(v, check, eps99, eps99);
                              break;

                              case 1:
                              // w.r.t previous v1
                              check = (1.0 + vj) * ap1 * delta_t;
                              ok &= CppAD::NearEqual(v, check, eps99, eps99);
                              break;

                              case 2:
                              // w.r.t previous v2
                              ok &= v == -1.0;
                              break;

                              case 3:
                              // w.r.t previous v3
                              ok &= v == - 2.0 * delta_t;
                              break;

                              default:
                              ok = false;
                         }
                    }
                    if( i_4 == 3 )
                    {    // derivative of this v3
                         switch( j_4 )
                         {
                              case 0:
                              // w.r.t previous v0
                              check = vj * ap0 * delta_t;
                              ok &= CppAD::NearEqual(v, check, eps99, eps99);
                              break;

                              case 1:
                              // w.r.t previous v1
                              check = (a + vj * ap1) * delta_t;
                              ok &= CppAD::NearEqual(v, check, eps99, eps99);
                              break;

                              case 2:
                              // w.r.t previous v2
                              ok &= v == 2.0 * delta_t;
                              break;

                              case 3:
                              // w.r.t previous v3
                              ok &= v == -1.0;
                              break;

                              default:
                              ok = false;
                         }
                    }
               }
          }
     }
     return ok;
}

Input File: src/utility.hpp