Changeset 1566


Ignore:
Timestamp:
Oct 28, 2009 12:21:25 PM (11 years ago)
Author:
bradbell
Message:

trunk: Make ipopt_ode_simple and introduction to ipopt_ode_fast.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/ipopt_cppad/ipopt_ode_fast.cpp

    r1565 r1566  
    1010Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1111-------------------------------------------------------------------------- */
    12 
     12/*
     13$begin ipopt_ode_fast.cpp$$
     14$spell
     15        ipopt_cppad_nlp
     16        nd
     17        ny
     18        na
     19        ns
     20$$
     21
     22$section ODE Fitting Using Fast Representation: Source Code$$
     23
     24$index ipopt_cppad_nlp, ode example source$$
     25$index ode, ipopt_cppad_nlp example source$$
     26$index example, ipopt_cppad_nlp ode source$$
     27$index source, ipopt_cppad_nlp ode example$$
     28
     29$code
     30$verbatim%ipopt_cppad/ipopt_ode_fast.cpp%0%// BEGIN PROGRAM%// END PROGRAM%1%$$
     31$$
     32
     33$end
     34*/
     35
     36// BEGIN PROGRAM
    1337# include "ipopt_cppad_nlp.hpp"
    14 
    15 // include a definition for Number.
    16 typedef Ipopt::Number Number;
    1738
    1839namespace {
    1940        //------------------------------------------------------------------
     41        // This section of the code is the same in ipopt_ode_simple.cpp
    2042        // simulated data
     43        typedef Ipopt::Number Number;
    2144        Number a0 = 1.;  // simulation value for a[0]
    2245        Number a1 = 2.;  // simulation value for a[1]
     
    3356        // Simulated data for case with no noise (first point is not used)
    3457        double z[] = { 0.0,  y_one(0.5), y_one(1.0), y_one(1.5), y_one(2.0) };
    35 }
    36 // ---------------------------------------------------------------------------
    37 namespace { // Begin empty namespace
    38 
     58        // Number of time grid points for each measurement interval
     59        size_t N[] = { 0,            5,          5,          5,          5  };
     60        // S[i] = N[0] + ... + N[i]
     61        size_t S[] = { 0,            5,         10,         15,         20  };
     62        // Number of measurement values
     63        size_t Nz  = sizeof(z) / sizeof(z[0]) - 1;
     64        // Number of components in the function y(t, a)
     65        size_t Ny  = 2;
     66        // Number of components in the vectro a
     67        size_t Na  = 3;
     68
     69        // Initial Condition function, F(a) = y(t, a) at t = 0
     70        // (for this particular example)
     71        template <class Vector>
     72        Vector eval_F(const Vector &a)
     73        {       Vector F(Ny);
     74                // y_0 (t) = a[0]*exp(-a[1] * t)
     75                F[0] = a[0];
     76                // y_1 (t) =
     77                // a[0]*a[1]*(exp(-a[2] * t) - exp(-a[1] * t))/(a[1] - a[2])
     78                F[1] = 0.;
     79                return F;
     80        }
     81        // G(y, a) =  \partial_t y(t, a); i.e. the differential equation
     82        // (for this particular example)
     83        template <class Vector>
     84        Vector eval_G(const Vector &y , const Vector &a)
     85        {       Vector G(Ny);
     86                // y_0 (t) = a[0]*exp(-a[1] * t)
     87                G[0] = -a[1] * y[0]; 
     88                // y_1 (t) =
     89                // a[0]*a[1]*(exp(-a[2] * t) - exp(-a[1] * t))/(a[1] - a[2])
     90                G[1] = +a[1] * y[0] - a[2] * y[1];
     91                return G;
     92        }
     93        // H(i, y, a) = contribution to objective at i-th data point
     94        // (for this particular example)
     95        template <class Scalar, class Vector>
     96        Scalar eval_H(size_t i, const Vector &y, const Vector &a)
     97        {       // This particular H is for a case where y_1 (t) is measured
     98                Scalar diff = z[i] - y[1];
     99                return diff * diff;
     100        }
     101
     102        //------------------------------------------------------------------
     103        // This section of the code is different in ipopt_ode_simple.cpp
    39104size_t nd = sizeof(s)/sizeof(s[0]) - 1; // number of actual data values
    40105size_t ny = 2;   // dimension of y(t, a)
     
    42107size_t ns = 5;   // number of grid intervals between each data value
    43108
    44 // F(a) = y(0, a); i.e., initial condition
    45 template <class Vector>
    46 Vector eval_F(const Vector &a)
    47 {       // This particual F is a case where ny == 2 and na == 3
    48         Vector F(ny);
    49         // y_0 (t) = a[0]*exp(-a[1] * t)
    50         F[0] = a[0];
    51         // y_1 (t) = a[0]*a[1]*(exp(-a[2] * t) - exp(-a[1] * t))/(a[1] - a[2])
    52         F[1] = 0.;
    53         return F;
    54 }
    55 // G(y, a) =  y'(t, a); i.e. ODE
    56 template <class Vector>
    57 Vector eval_G(const Vector &y , const Vector &a)
    58 {       // This particular G is for a case where ny == 2 and na == 3
    59         Vector G(ny);
    60         // y_0 (t) = a[0]*exp(-a[1] * t)
    61         G[0] = -a[1] * y[0]; 
    62         // y_1 (t) = a[0]*a[1]*(exp(-a[2] * t) - exp(-a[1] * t))/(a[1] - a[2])
    63         G[1] = +a[1] * y[0] - a[2] * y[1];
    64         return G;
    65 }
    66 // H(k, y, a) = contribution to objective at k-th data point
    67 template <class Scalar, class Vector>
    68 Scalar eval_H(size_t k, const Vector &y, const Vector &a)
    69 {       // This particular H is for a case where y_1 (t) is measured
    70         Scalar diff = z[k] - y[1];
    71         return diff * diff;
    72 }
    73 
    74 // -----------------------------------------------------------------------------
    75 class FG_info : public ipopt_cppad_fg_info
    76 {
    77 private:
    78         bool retape_;
    79 public:
    80         // derived class part of constructor
    81         FG_info(bool retape)
    82         : retape_ (retape)
    83         { }
    84         // r^k for k = 0, 1, ..., nd-1 corresponds to the data value
    85         // r^k for k = nd corresponds to initial condition
    86         // r^k for k = nd+1 , ... , 2*nd is used for trapezoidal approximation
    87         size_t number_functions(void)
    88         {       return nd + 1 + nd; }
    89         ADVector eval_r(size_t k, const ADVector &u)
    90         {       size_t j;
    91                 ADVector y_M(ny), a(na);
    92                 if( k < nd )
    93                 {       // r^k for k = 0, ... , nd-1
    94                         // We use a differnent k for each data point
    95                         ADVector r(1); // return value is a scalar
    96                         size_t j;
    97                         // u is [y( s[k+1] ) , a]
    98                         for(j = 0; j < ny; j++)
    99                                 y_M[j] = u[j];
     109        // ------------------------------------------------------------------
     110        class FG_info : public ipopt_cppad_fg_info
     111        {
     112        private:
     113                bool retape_;
     114        public:
     115                // derived class part of constructor
     116                FG_info(bool retape)
     117                : retape_ (retape)
     118                { }
     119                // r^k for k = 0, 1, ..., Nz-1 used for measurements
     120                // r^k for k = Nz              use for initial condition
     121                // r^k for k = Nz+1, ..., 2*Nz used for trapezoidal approx
     122                size_t number_functions(void)
     123                {       return Nz + 1 + Nz; }
     124                ADVector eval_r(size_t k, const ADVector &u)
     125                {       size_t j;
     126                        ADVector y(Ny), a(Na);
     127                        // objective function --------------------------------
     128                        if( k < Nz )
     129                        {       // used for measurement with index k+1
     130                                ADVector r(1); // return value is a scalar
     131                                size_t j;
     132                                // u is [y( s[k+1] ) , a]
     133                                for(j = 0; j < Ny; j++)
     134                                        y[j] = u[j];
     135                                for(j = 0; j < Na; j++)
     136                                        a[j] = u[Ny + j];
     137                                r[0] = eval_H<ADNumber>(k+1, y, a);
     138                                return r;
     139                        }
     140                        // initial condition ---------------------------------
     141                        if( k == Nz )
     142                        {       ADVector r(Ny), F(Ny);
     143                                // u is [y(t), a] at t = 0
     144                                for(j = 0; j < Ny; j++)
     145                                        y[j] = u[j];
     146                                for(j = 0; j < Na; j++)
     147                                        a[j] = u[Ny + j];
     148                                F    = eval_F(a);
     149                                for(j = 0; j < Ny; j++)
     150                                        r[j]   = y[j] - F[j];
     151                                return  r;
     152                        }
     153                        // trapezoidal approximation -------------------------
     154                        ADVector ym(Ny), G(Ny), Gm(Ny), r(Ny);
     155                        // r^k for k = Nz+1, ... , 2*Nz
     156                        // interval between data samples
     157                        Number T = s[k-Nz] - s[k-Nz-1];
     158                        // integration step size
     159                        Number dt = T / Number( N[k-Nz] );
     160                        // u = [ y(t[i-1], a) , y(t[i], a), a )
     161                        for(j = 0; j < Ny; j++)
     162                        {       ym[j] = u[j];
     163                                y[j]  = u[Ny + j];
     164                        }
     165                        for(j = 0; j < Na; j++)
     166                                a[j] = u[2 * Ny + j];
     167                        Gm  = eval_G(ym, a);
     168                        G   = eval_G(y,  a);
     169                        for(j = 0; j < Ny; j++)
     170                                r[j] = y[j] - ym[j] - (G[j] + Gm[j]) * dt / 2.;
     171                        return r;
     172                }
     173                // The operations sequence for r_eval does not depend on u,
     174                // hence retape = false should work and be faster.
     175                bool retape(size_t k)
     176                {       return retape_; }
     177                // size of the vector u in eval_r
     178                size_t domain_size(size_t k)
     179                {       if( k < Nz )
     180                                return Ny + Na;   // objective function
     181                        if( k == nd )
     182                                return Ny + Na;  // initial value constraint
     183                        return 2 * Ny + Na;      // trapezodial constraints
     184                }
     185                // size of the return value from eval_r
     186                size_t range_size(size_t k)
     187                {       if( k < Nz )
     188                                return 1;
     189                        return Ny;
     190                }
     191                // number of terms that use this value of k
     192                size_t number_terms(size_t k)
     193                {       if( k <= Nz )
     194                                return 1;  // r^k used once for k <= Nz
     195                        // r^k used N[k-Nz] times for k > Nz
     196                        return N[k-Nz];
     197                }
     198                void index(size_t k, size_t ell, SizeVector& I, SizeVector& J)
     199                {       size_t i, j;
     200                        // # of components of x corresponding to values for y
     201                        size_t ny_inx = (S[Nz] + 1) * Ny;
     202                        // objective function -------------------------------
     203                        if( k < nd )
     204                        {       // index in fg corresponding to objective       
     205                                I[0] = 0;
     206                                // u = [ y(t, a) , a ]
     207                                // The first ny components of u is y(t) at
     208                                //      t = s[k+1] = t[S[k+1]]
     209                                // x indices corresponding to this value of y
     210                                for(j = 0; j < Ny; j++)
     211                                        J[j] = S[k + 1] * Ny + j;
     212                                // components of x correspondig to a
     213                                for(j = 0; j < na; j++)
     214                                        J[Ny + j] = ny_inx + j;
     215                                return;
     216                        }
     217                        // initial conditions --------------------------------
     218                        if( k == nd )
     219                        {       // index in fg for inidial condition constraint
     220                                for(j = 0; j < Ny; j++)
     221                                        I[j] = 1 + j;
     222                                // u = [ y(t, a) , a ] where t = 0
     223                                // x indices corresponding to this value of y
     224                                for(j = 0; j < Ny; j++)
     225                                        J[j] = j;
     226                                // following that, u contains the vector a
     227                                for(j = 0; j < Na; j++)
     228                                        J[Ny + j] = ny_inx + j;
     229                                return;
     230                        }
     231                        // trapoziodal approximation -------------------------
     232                        // index of first grid point in this approximation
     233                        i = S[k - Nz - 1]  + ell;
     234                        // There are Ny difference equations for each time
     235                        // point.  Add one for the objective function, and Ny
     236                        // for the initial value constraints.
     237                        for(j = 0; j < Ny; j++)
     238                                I[j] = 1 + Ny + i * Ny + j;
     239                        // u = [ y(t, a) , y(t+dt, a) , a ] at t = t[i]
     240                        for(j = 0; j < Ny; j++)
     241                        {       J[j]      = i * Ny  + j; // y^i indices
     242                                J[Ny + j] = J[j] + Ny;   // y^{i+1} indices
     243                        }
    100244                        for(j = 0; j < na; j++)
    101                                 a[j] = u[ny + j];
    102                         r[0] = eval_H<ADNumber>(k+1, y_M, a);
    103                         return r;
    104                 }
    105                 if( k == nd )
    106                 {       // r^k for k = nd corresponds to initial condition
    107                         ADVector r(ny), F(ny);
    108                         // u is [y(0), a]
    109                         for(j = 0; j < na; j++)
    110                                 a[j] = u[ny + j];
    111                         F    = eval_F(a);
    112                         for(j = 0; j < ny; j++)
    113                         {       y_M[j] = u[j];
    114                                 // y(0) - F(a)
    115                                 r[j]   = y_M[j] - F[j];
    116                         }
    117                         return  r;
    118                 }
    119                 // r^k for k = nd+1, ... , 2*nd
    120                 // corresponds to trapezoidal approximations in the
    121                 // data interval [ s[k-nd] , s[k-nd-1] ]
    122                 ADVector y_M1(ny);
    123                 // u = [y_M, y_M1, a] where y_M is y(t) at
    124                 // t = t[ (k-nd-1) * ns + ell ]
    125                 for(j = 0; j < ny; j++)
    126                 {       y_M[j]  = u[j];
    127                         y_M1[j] = u[ny + j];
    128                 }
    129                 for(j = 0; j < na; j++)
    130                         a[j] = u[2 * ny + j];
    131                 Number dt      = (s[k-nd] - s[k-nd-1]) / Number(ns);
    132                 ADVector G_M   = eval_G(y_M,  a);
    133                 ADVector G_M1  = eval_G(y_M1, a);
    134                 ADVector r(ny);
    135                 for(j = 0; j < ny; j++)
    136                         r[j] = y_M1[j] - y_M[j] - (G_M1[j] + G_M[j]) * dt/2.;
    137                 return r;
    138         }
    139         // Operation sequence does not depend on u so retape = false should
    140         // work and be faster. Allow for both for testing.
    141         bool retape(size_t k)
    142         {       return retape_; }
    143         // size of the vector u in eval_r
    144         size_t domain_size(size_t k)
    145         {       if( k < nd )
    146                         return ny + na;   // objective function
    147                 if( k == nd )
    148                         return ny + na;  // initial value constraint
    149                 return 2 * ny + na;      // trapezodial constraints
    150         }
    151         // size of the vector r in eval_r
    152         size_t range_size(size_t k)
    153         {       if( k < nd )
    154                         return 1;
    155                 return ny;
    156         }
    157         size_t number_terms(size_t k)
    158         {       if( k <= nd )
    159                         return 1;     // r^k used once for k <= nd
    160                 return ns;            // r^k used ns times for k > nd
    161         }
    162         void index(size_t k, size_t ell, SizeVector& I, SizeVector& J)
    163         {       size_t i, j;
    164                 // number of components of x corresponding to value of y
    165                 size_t ny_inx = (nd * ns + 1) * ny;
    166                 if( k < nd )
    167                 {       // r^k for k = 0 , ... , nd-1 corresponds to objective
    168                         I[0] = 0;
    169                         // The first ny components of u is y(t) at
    170                         //      t = s[k+1] = t[(k+1)*ns]
    171                         // components of x corresponding to this value of y
    172                         for(j = 0; j < ny; j++)
    173                                 J[j] = (k + 1) * ns * ny + j;
    174                         // components of x correspondig to a
    175                         for(j = 0; j < na; j++)
    176                                 J[ny + j] = ny_inx + j;
    177                         return;
    178                 }
    179                 if( k == nd )
    180                 {       // r^k corresponds to initial condition
    181                         for(i = 0; i < ny; i++)
    182                                 I[i] = 1 + i;
    183                         // u starts with the first j components of x
    184                         // (which corresponding to y(t) at t[0])
    185                         for(j = 0; j < ny; j++)
    186                                 J[j] = j;
    187                         // following that, u contains the vector a
    188                         for(j = 0; j < na; j++)
    189                                 J[ny + j] = ny_inx + j;
    190                         return;
    191                 }
    192                 // index of first grid point in ts for difference equation
    193                 size_t M = (k - nd - 1) * ns + ell;
    194                 for(j = 0; j < ny; j++)
    195                 {       J[j]          = M * ny  + j; // index of y_M in x
    196                         J[ny + j]     = J[j] + ny;   // index of y_M1
    197                 }
    198                 for(j = 0; j < na; j++)
    199                         J[2 * ny + j] = ny_inx + j;                      // a
    200                 // There are ny difference equations for each grid point.
    201                 // Add one for the objective function index, and ny for the
    202                 // initial value constraint.
    203                 for(i = 0; i < ny; i++)
    204                         I[i] = 1 + ny + M * ny + i ;
    205         }
    206 };
     245                                J[2 * Ny + j] = ny_inx + j; // a indices
     246                }
     247        };
    207248
    208249} // End empty namespace
     250
     251
    209252// ---------------------------------------------------------------------------
    210 
    211253bool ipopt_ode_fast(void)
     254// The rest of the code is the same as in ipopt_ode_simple.cpp
    212255{       bool ok = true;
    213         size_t j, I;
    214 
    215         // number of components of x corresponding to value of y
    216         size_t ny_inx = (nd * ns + 1) * ny;
     256        size_t i, j;
     257
     258        // number of components of x corresponding to values for y
     259        size_t ny_inx = (S[Nz] + 1) * Ny;
    217260        // number of constraints (range dimension of g)
    218         size_t m = ny + nd * ns * ny;
     261        size_t m      = ny_inx;
    219262        // number of components in x (domain dimension for f and g)
    220         size_t n = ny_inx + na;
     263        size_t n      = ny_inx + Na;
    221264        // the argument vector for the optimization is
    222         // y(t) at t[0] , ... , t[nd*ns] , followed by a
     265        // y(t) at t[0] , ... , t[S[Nz]] , followed by a
    223266        NumberVector x_i(n), x_l(n), x_u(n);
    224267        for(j = 0; j < ny_inx; j++)
     
    227270                x_u[j] = +1.0e19;  // no upper limit
    228271        }
    229         for(j = 0; j < na; j++)
     272        for(j = 0; j < Na; j++)
    230273        {       x_i[ny_inx + j ] = .5;       // initiali a for optimization
    231274                x_l[ny_inx + j ] =  -1.e19;  // no lower limit
     
    234277        // all of the difference equations are constrained to the value zero
    235278        NumberVector g_l(m), g_u(m);
    236         for(I = 0; I < m; I++)
    237         {       g_l[I] = 0.;
    238                 g_u[I] = 0.;
     279        for(i = 0; i < m; i++)
     280        {       g_l[i] = 0.;
     281                g_u[i] = 0.;
    239282        }
    240283        // derived class object
     
    242285        for(size_t icase = 0; icase <= 1; icase++)
    243286        {       // Retaping is slow, so only do icase = 0 for large values
    244                 // of ns.
     287                // of N[].
    245288                bool retape = icase != 0;
    246289
     
    270313
    271314                // Derivative testing is very slow for large problems
    272                 // so comment this out if you use a large value for ns.
     315                // so comment this out if you use a large value for N[].
    273316                app->Options()-> SetStringValue(
    274317                        "derivative_test", "second-order"
     
    284327
    285328                // split out return values
    286                 NumberVector a(na), y_0(ny), y_1(ny), y_2(ny);
    287                 for(j = 0; j < na; j++)
     329                NumberVector a(Na), y_0(Ny), y_1(Ny), y_2(Ny);
     330                for(j = 0; j < Na; j++)
    288331                        a[j] = solution.x[ny_inx+j];
    289                 for(j = 0; j < ny; j++)
     332                for(j = 0; j < Ny; j++)
    290333                {       y_0[j] = solution.x[j];
    291                         y_1[j] = solution.x[ny + j];
    292                         y_2[j] = solution.x[2 * ny + j];
     334                        y_1[j] = solution.x[Ny + j];
     335                        y_2[j] = solution.x[2 * Ny + j];
    293336                }
    294337
    295338                // Check some of the return values
    296                 Number rel_tol = 1e-2; // use a larger value of ns
     339                Number rel_tol = 1e-2; // use a larger value of N[]
    297340                Number abs_tol = 1e-2; // to get better accuracy here.
    298341                Number check_a[] = {a0, a1, a2}; // see the y_one function
    299                 for(j = 0; j < na; j++)
     342                for(j = 0; j < Na; j++)
    300343                {
    301344                        ok &= CppAD::NearEqual(
     
    308351                // check the initial value constraint
    309352                NumberVector F = eval_F(a);
    310                 for(j = 0; j < ny; j++)
     353                for(j = 0; j < Ny; j++)
    311354                        ok &= CppAD::NearEqual(F[j], y_0[j], rel_tol, abs_tol);
    312355
     
    314357                NumberVector G_0 = eval_G(y_0, a);
    315358                NumberVector G_1 = eval_G(y_1, a);
    316                 Number dt = (s[1] - s[0]) / Number(ns);
     359                Number dt = (s[1] - s[0]) / Number(N[1]);
    317360                Number check;
    318                 for(j = 0; j < ny; j++)
     361                for(j = 0; j < Ny; j++)
    319362                {       check = y_1[j] - y_0[j] - (G_1[j]+G_0[j])*dt/2;
    320363                        ok &= CppAD::NearEqual( check, 0., rel_tol, abs_tol);
     
    323366                // check the second trapezoidal equation
    324367                NumberVector G_2 = eval_G(y_2, a);
    325                 if( ns == 1 )
    326                         dt = (s[2] - s[1]) / Number(ns);
    327                 for(j = 0; j < ny; j++)
     368                if( N[1] == 1 )
     369                        dt = (s[2] - s[1]) / Number(N[2]);
     370                for(j = 0; j < Ny; j++)
    328371                {       check = y_2[j] - y_1[j] - (G_2[j]+G_1[j])*dt/2;
    329372                        ok &= CppAD::NearEqual( check, 0., rel_tol, abs_tol);
     
    332375                // check the objective function (specialized to this case)
    333376                check = 0.;
    334                 NumberVector y_M(ny);
    335                 for(size_t k = 0; k < nd; k++)
    336                 {       for(j = 0; j < ny; j++)
    337                         {       size_t M = (k + 1) * ns;
    338                                 y_M[j] =  solution.x[M * ny + j];
     377                NumberVector y_M(Ny);
     378                for(size_t k = 0; k < Nz; k++)
     379                {       for(j = 0; j < Ny; j++)
     380                        {       size_t M = S[k + 1];
     381                                y_M[j] =  solution.x[M * Ny + j];
    339382                        }
    340383                        check += eval_H<Number>(k + 1, y_M, a);
     
    345388        return ok;
    346389}
     390// END PROGRAM
  • trunk/ipopt_cppad/ipopt_ode_simple.cpp

    r1565 r1566  
    1010Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
    1111-------------------------------------------------------------------------- */
    12 
     12/*
     13$begin ipopt_ode_simple.cpp$$
     14$spell
     15        ipopt_cppad_nlp
     16        Nz
     17        Ny
     18        Na
     19$$
     20
     21$section ODE Fitting Using Simple Representation: Source Code$$
     22
     23$index ipopt_cppad_nlp, ode example source$$
     24$index ode, ipopt_cppad_nlp example source$$
     25$index example, ipopt_cppad_nlp ode source$$
     26$index source, ipopt_cppad_nlp ode example$$
     27
     28$code
     29$verbatim%ipopt_cppad/ipopt_ode_simple.cpp%0%// BEGIN PROGRAM%// END PROGRAM%1%$$
     30$$
     31
     32$end
     33*/
     34
     35// BEGIN PROGRAM
    1336# include "ipopt_cppad_nlp.hpp"
    1437
    1538namespace {
     39        //------------------------------------------------------------------
     40        // This section of the code is the same in ipopt_ode_fast.cpp
     41        // simulated data
    1642        typedef Ipopt::Number Number;
    17         //------------------------------------------------------------------
    18         // simulated data
    1943        Number a0 = 1.;  // simulation value for a[0]
    2044        Number a1 = 2.;  // simulation value for a[1]
     
    3155        // Simulated data for case with no noise (first point is not used)
    3256        double z[] = { 0.0,  y_one(0.5), y_one(1.0), y_one(1.5), y_one(2.0) };
    33 
    34         // F(a) = y(0, a); i.e., initial condition
     57        // Number of time grid points for each measurement interval
     58        size_t N[] = { 0,            5,          5,          5,          5  };
     59        // S[i] = N[0] + ... + N[i]
     60        size_t S[] = { 0,            5,         10,         15,         20  };
     61        // Number of measurement values
     62        size_t Nz  = sizeof(z) / sizeof(z[0]) - 1;
     63        // Number of components in the function y(t, a)
     64        size_t Ny  = 2;
     65        // Number of components in the vectro a
     66        size_t Na  = 3;
     67
     68        // Initial Condition function, F(a) = y(t, a) at t = 0
     69        // (for this particular example)
    3570        template <class Vector>
    3671        Vector eval_F(const Vector &a)
    37         {       // This particual F is a case where ny == 2 and na == 3
    38                 Vector F(2);
     72        {       Vector F(Ny);
    3973                // y_0 (t) = a[0]*exp(-a[1] * t)
    4074                F[0] = a[0];
     
    4478                return F;
    4579        }
    46         // G(y, a) =  y'(t, a); i.e. ODE
     80        // G(y, a) =  \partial_t y(t, a); i.e. the differential equation
     81        // (for this particular example)
    4782        template <class Vector>
    4883        Vector eval_G(const Vector &y , const Vector &a)
    49         {       // This particular G is for a case where ny == 2 and na == 3
    50                 Vector G(2);
     84        {       Vector G(Ny);
    5185                // y_0 (t) = a[0]*exp(-a[1] * t)
    5286                G[0] = -a[1] * y[0]; 
     
    5690                return G;
    5791        }
    58         // H(k, y, a) = contribution to objective at k-th data point
     92        // H(i, y, a) = contribution to objective at i-th data point
     93        // (for this particular example)
    5994        template <class Scalar, class Vector>
    60         Scalar eval_H(size_t k, const Vector &y, const Vector &a)
     95        Scalar eval_H(size_t i, const Vector &y, const Vector &a)
    6196        {       // This particular H is for a case where y_1 (t) is measured
    62                 Scalar diff = z[k] - y[1];
     97                Scalar diff = z[i] - y[1];
    6398                return diff * diff;
    6499        }
    65100
    66         size_t nd = sizeof(s)/sizeof(s[0]) - 1; // number of actual data values
    67         size_t ny = 2;   // dimension of y(t, a)
    68         size_t na = 3;   // dimension of a
    69         size_t ns = 5;   // number of grid intervals between each data value
    70 
     101        //------------------------------------------------------------------
     102        // This section of the code is different in ipopt_ode_fast.cpp
    71103        class FG_info : public ipopt_cppad_fg_info
    72104        {
     
    84116                        size_t i, j, k;
    85117                        // # of components of x corresponding to values for y
    86                         size_t ny_inx = (nd * ns + 1) * ny;
     118                        size_t ny_inx = (S[Nz] + 1) * Ny;
    87119                        // # of constraints (range dimension of g)
    88                         size_t m = ny + nd * ns * ny;
     120                        size_t m = ny_inx;
    89121                        // # of components in x (domain dimension for f and g)
    90                         size_t n = ny_inx + na;
     122                        size_t n = ny_inx + Na;
    91123                        assert ( x.size() == n );
    92124                        // vector for return value
    93125                        ADVector fg(m + 1);
    94126                        // vector of parameters
    95                         ADVector a(na);
    96                         for(i = 0; i < na; i++)
    97                                 a[i] = x[ny_inx + i];
     127                        ADVector a(Na);
     128                        for(j = 0; j < Na; j++)
     129                                a[j] = x[ny_inx + j];
    98130                        // vector for value of y(t)
    99                         ADVector y(ny);
    100                         // objective function
     131                        ADVector y(Ny);
     132                        // objective function -------------------------------
    101133                        fg[0] = 0.;
    102                         for(k = 0; k < nd; k++)
    103                         {       for(i = 0; i < ny; i++)
    104                                         y[i] = x[ny*(k+1)*ns + i];
     134                        for(k = 0; k < Nz; k++)
     135                        {       for(j = 0; j < Ny; j++)
     136                                        y[j] = x[Ny*S[k+1] + j];
    105137                                fg[0] += eval_H<ADNumber>(k+1, y, a);
    106138                        }
    107                         // initial condition y(0) = F(a)
     139                        // initial condition ---------------------------------
    108140                        ADVector F = eval_F(a);
    109                         for(i = 0; i < ny; i++)
    110                         {       y[i]    = x[i];
    111                                 fg[1+i] = y[i] - F[i];
     141                        for(j = 0; j < Ny; j++)
     142                        {       y[j]    = x[j];
     143                                fg[1+j] = y[j] - F[j];
    112144                        }
    113                         // compute constraints corresponding to trapezoidal
    114                         // approximation for the solution of the ODE
    115                         ADVector ym(ny), G(ny), Gm(ny);
     145                        // trapezoidal approximation --------------------------
     146                        ADVector ym(Ny), G(Ny), Gm(Ny);
    116147                        G = eval_G(y, a);
    117148                        ADNumber dy;
    118                         for(k = 0; k < nd; k++)
     149                        for(k = 0; k < Nz; k++)
    119150                        {       // interval between data points
    120                                 Number T  = (s[k+1] - s[k]);
     151                                Number T  = s[k+1] - s[k];
    121152                                // integration step size
    122                                 Number dt = T / Number(ns);
    123                                 for(j = 0; j < ns; j++)
    124                                 {       size_t index = (j + k*ns) * ny;
     153                                Number dt = T / Number( N[k+1] );
     154                                for(j = 0; j < N[k+1]; j++)
     155                                {       size_t index = (j + S[k]) * Ny;
    125156                                        // y(t) at end of last step
    126157                                        ym = y;
     
    128159                                        Gm = G;
    129160                                        // value of y(t) at end of this step
    130                                         for(i = 0; i < ny; i++)
    131                                                 y[i] = x[ny + index + i];
     161                                        for(i = 0; i < Ny; i++)
     162                                                y[i] = x[Ny + index + i];
    132163                                        // G(y, a) at end of this step
    133164                                        G = eval_G(y, a);
    134165                                        // trapezoidal approximation residual
    135                                         for(i = 0; i < ny; i++)
     166                                        for(i = 0; i < Ny; i++)
    136167                                        {       dy = (G[i] + Gm[i]) * dt / 2;
    137                                                 fg[1+ny+index+i] =
     168                                                fg[1+Ny+index+i] =
    138169                                                        y[i] - ym[i] - dy;
    139170                                        }
     
    142173                        return fg;
    143174                }
     175                // The operations sequence for r_eval does not depend on u,
     176                // hence retape = false should work and be faster.
    144177                bool retape(size_t k)
    145178                {       return retape_; }
     
    149182       
    150183
     184// ---------------------------------------------------------------------------
     185// The rest of the code is the same as in ipopt_ode_fast.cpp
    151186bool ipopt_ode_simple(void)
    152187{       bool ok = true;
    153         size_t j, I;
    154 
    155         // number of components of x corresponding to value of y
    156         size_t ny_inx = (nd * ns + 1) * ny;
     188        size_t i, j;
     189
     190        // number of components of x corresponding to values for y
     191        size_t ny_inx = (S[Nz] + 1) * Ny;
    157192        // number of constraints (range dimension of g)
    158         size_t m = ny + nd * ns * ny;
     193        size_t m      = ny_inx;
    159194        // number of components in x (domain dimension for f and g)
    160         size_t n = ny_inx + na;
     195        size_t n      = ny_inx + Na;
    161196        // the argument vector for the optimization is
    162         // y(t) at t[0] , ... , t[nd*ns] , followed by a
     197        // y(t) at t[0] , ... , t[S[Nz]] , followed by a
    163198        NumberVector x_i(n), x_l(n), x_u(n);
    164199        for(j = 0; j < ny_inx; j++)
     
    167202                x_u[j] = +1.0e19;  // no upper limit
    168203        }
    169         for(j = 0; j < na; j++)
     204        for(j = 0; j < Na; j++)
    170205        {       x_i[ny_inx + j ] = .5;       // initiali a for optimization
    171206                x_l[ny_inx + j ] =  -1.e19;  // no lower limit
     
    174209        // all of the difference equations are constrained to the value zero
    175210        NumberVector g_l(m), g_u(m);
    176         for(I = 0; I < m; I++)
    177         {       g_l[I] = 0.;
    178                 g_u[I] = 0.;
     211        for(i = 0; i < m; i++)
     212        {       g_l[i] = 0.;
     213                g_u[i] = 0.;
    179214        }
    180215        // derived class object
     
    182217        for(size_t icase = 0; icase <= 1; icase++)
    183218        {       // Retaping is slow, so only do icase = 0 for large values
    184                 // of ns.
     219                // of N[].
    185220                bool retape = icase != 0;
    186221
     
    210245
    211246                // Derivative testing is very slow for large problems
    212                 // so comment this out if you use a large value for ns.
     247                // so comment this out if you use a large value for N[].
    213248                app->Options()-> SetStringValue(
    214249                        "derivative_test", "second-order"
     
    224259
    225260                // split out return values
    226                 NumberVector a(na), y_0(ny), y_1(ny), y_2(ny);
    227                 for(j = 0; j < na; j++)
     261                NumberVector a(Na), y_0(Ny), y_1(Ny), y_2(Ny);
     262                for(j = 0; j < Na; j++)
    228263                        a[j] = solution.x[ny_inx+j];
    229                 for(j = 0; j < ny; j++)
     264                for(j = 0; j < Ny; j++)
    230265                {       y_0[j] = solution.x[j];
    231                         y_1[j] = solution.x[ny + j];
    232                         y_2[j] = solution.x[2 * ny + j];
     266                        y_1[j] = solution.x[Ny + j];
     267                        y_2[j] = solution.x[2 * Ny + j];
    233268                }
    234269
    235270                // Check some of the return values
    236                 Number rel_tol = 1e-2; // use a larger value of ns
     271                Number rel_tol = 1e-2; // use a larger value of N[]
    237272                Number abs_tol = 1e-2; // to get better accuracy here.
    238273                Number check_a[] = {a0, a1, a2}; // see the y_one function
    239                 for(j = 0; j < na; j++)
     274                for(j = 0; j < Na; j++)
    240275                {
    241276                        ok &= CppAD::NearEqual(
     
    248283                // check the initial value constraint
    249284                NumberVector F = eval_F(a);
    250                 for(j = 0; j < ny; j++)
     285                for(j = 0; j < Ny; j++)
    251286                        ok &= CppAD::NearEqual(F[j], y_0[j], rel_tol, abs_tol);
    252287
     
    254289                NumberVector G_0 = eval_G(y_0, a);
    255290                NumberVector G_1 = eval_G(y_1, a);
    256                 Number dt = (s[1] - s[0]) / Number(ns);
     291                Number dt = (s[1] - s[0]) / Number(N[1]);
    257292                Number check;
    258                 for(j = 0; j < ny; j++)
     293                for(j = 0; j < Ny; j++)
    259294                {       check = y_1[j] - y_0[j] - (G_1[j]+G_0[j])*dt/2;
    260295                        ok &= CppAD::NearEqual( check, 0., rel_tol, abs_tol);
     
    263298                // check the second trapezoidal equation
    264299                NumberVector G_2 = eval_G(y_2, a);
    265                 if( ns == 1 )
    266                         dt = (s[2] - s[1]) / Number(ns);
    267                 for(j = 0; j < ny; j++)
     300                if( N[1] == 1 )
     301                        dt = (s[2] - s[1]) / Number(N[2]);
     302                for(j = 0; j < Ny; j++)
    268303                {       check = y_2[j] - y_1[j] - (G_2[j]+G_1[j])*dt/2;
    269304                        ok &= CppAD::NearEqual( check, 0., rel_tol, abs_tol);
     
    272307                // check the objective function (specialized to this case)
    273308                check = 0.;
    274                 NumberVector y_M(ny);
    275                 for(size_t k = 0; k < nd; k++)
    276                 {       for(j = 0; j < ny; j++)
    277                         {       size_t M = (k + 1) * ns;
    278                                 y_M[j] =  solution.x[M * ny + j];
     309                NumberVector y_M(Ny);
     310                for(size_t k = 0; k < Nz; k++)
     311                {       for(j = 0; j < Ny; j++)
     312                        {       size_t M = S[k + 1];
     313                                y_M[j] =  solution.x[M * Ny + j];
    279314                        }
    280315                        check += eval_H<Number>(k + 1, y_M, a);
     
    285320        return ok;
    286321}
     322// END PROGRAM
  • trunk/ipopt_cppad/makefile.in

    r1565 r1566  
    9090CYGPATH_W = @CYGPATH_W@
    9191
    92 # $Id: makefile.am 1370 2009-05-31 05:31:50Z bradbell $
     92# $Id: makefile.am 1565 2009-10-27 15:00:55Z bradbell $
    9393# -----------------------------------------------------------------------------
    94 # CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-08 Bradley M. Bell
     94# CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
    9595#
    9696# CppAD is distributed under multiple licenses. This distribution is under
  • trunk/omh/ipopt_cppad_ode2.omh

    r1565 r1566  
    3030y_0 (0 , a) = F(a) = \left( \begin{array}{c} a_0 \\ 0 \end{array} \right)
    3131\] $$
    32 where $latex a \in \R^3 $$ is an unknown parameter vector
    33 and $latex F : \R^3 \rightarrow \R^2 $$ is defined by the equation above.
    34 Our forward problem is stated as follows: Given $latex a \in \R^3$$
    35 determine the value of $latex y ( t , a ) $$ for
    36 various values  of $latex t $$.
     32where $latex Na = 3$$ and $latex a \in \R^{Na} $$
     33is an unknown parameter vector.
     34The function and $latex F : \R^{Na} \rightarrow \R^{Ny} $$
     35is defined by the equation above (where $latex Ny = 2$$).
     36Our forward problem is stated as follows: Given $latex a \in \R^{Na}$$
     37determine the value of $latex y ( t , a ) $$, for $latex t \in R$$, that
     38solves the initial value problem above.
    3739
    3840
     
    4244        y^0 = y(0, a) = \left( \begin{array}{c} a_0 \\ 0 \end{array} \right)
    4345\] $$
    44 Given an approximation value $latex y^M $$ for $latex y ( s_M , a )$$,
    45 the a trapezoidal method approximates $latex y ( s_{M+1} , a )$$
    46 (denoted by $latex y^{M+1}$$ ) by solving the equation
    47 $latex \[
    48 y^{M+1}  =  y^M +
    49 \left[ G( y^M , a ) + G( y^{M+1} , a ) \right] * \frac{s_{M+1} - s_M }{ 2 }
    50 \] $$
    51 where $latex G : \R^2 \times \R^3 \rightarrow \R^2$$ is the
     46Given an approximation value $latex y^{i-1} $$ for $latex y ( t_{i-1} , a )$$,
     47the a trapezoidal method approximates $latex y ( t_i , a )$$
     48(denoted by $latex y^i$$ ) by solving the equation
     49$latex \[
     50y^i  =  y^{i-1} +
     51\left[ G( y^i , a ) + G( y^{i-1} , a ) \right] * \frac{t_i - t_{i-1} }{ 2 }
     52\] $$
     53where $latex G : \R^{Ny} \times \R^{Na} \rightarrow \R^{Ny}$$ is the
    5254function representing this ODE; i.e.
    5355$latex \[
     
    5961\] $$
    6062This $latex G(y, a)$$ is linear with respect to $latex y$$, hence
    61 the implicit equation defining $latex y^{M+1} $$ can be solved
     63the implicit equation defining $latex y^i $$ can be solved
    6264inverting the a set of linear equations.
    6365In the general case,
    6466where $latex G(y, a)$$ is non-linear with respect to $latex y$$,
    65 an iterative procedure is used to calculate $latex y^{M+1}$$
    66 from $latex y^M$$.
     67an iterative procedure is used to calculate $latex y^i$$
     68from $latex y^{i-1}$$.
    6769
    6870$end
     
    7678
    7779$head Problem$$
    78 We define $latex y : \R \times \R^3 \rightarrow \R^2$$ as the
     80We define $latex y : \R \times \R^{Na} \rightarrow \R^{Ny}$$ as the
    7981solution of our $cref/ode forward problem/ipopt_cppad_ode_forward/Problem/$$.
    80 Suppose we are also given measurement values $latex z_k \in \R$$
    81 for $latex k = 1, 2, 3, 4$$ that are modeled by
    82 $latex \[
    83         z_k = y_1 ( s_k , a) + e_k
    84 \] $$
    85 where
    86 $latex s_k \in \R$$,
    87 $latex e_k \sim {\bf N} (0 , \sigma^2 ) $$,
    88 and $latex \sigma \in \R_+$$.
    89 The maximum likelihood estimate for $latex a$$ given
    90 $latex ( z_1 , z_2 , z_3 , z_4 )$$ solves the following inverse problem
     82Suppose we are also given measurement values $latex z \in \R^{Nz}$$
     83and for $latex i = 1, \ldots, Nz$$, the component
     84$latex \[
     85        z_i = y_1 ( s_i , a) + e_i
     86\] $$
     87where $latex s_i \in \R$$ is the time for the $th i$$ measurement,
     88$latex e_i \sim {\bf N} (0 , \sigma^2 ) $$ is the corresponding noise,
     89and $latex \sigma \in \R_+$$ is the corresponding standard deviation.
     90The maximum likelihood estimate for $latex a$$ given $latex z$$
     91solves the following inverse problem
    9192$latex \[
    9293\begin{array}{rcl}
    9394{\rm minimize} \;
    94         & \sum_{k=1}^4 H_k [ y( s_k , a ) , a ]
    95         & \;{\rm w.r.t} \; a \in \R^3
    96 \end{array}
    97 \] $$
    98 where the function $latex H_k : \R^2 \times \R^3 \rightarrow \R$$ is
     95        & \sum_{i=1}^{Nz} H_i [ y( s_i , a ) , a ]
     96        & \;{\rm w.r.t} \; a \in \R^{Na}
     97\end{array}
     98\] $$
     99where the functions
     100$latex H_i : \R^{Ny} \times \R^{Na} \rightarrow \R$$ is
    99101defined by
    100102$latex \[
    101         H_k (y, a) = ( z_k - y_1 )^2
     103        H_i (y, a) = ( z_i - y_1 )^2
    102104\] $$
    103105
     
    108110one uses a
    109111$cref/numerical procedure/ipopt_cppad_ode_forward/Numerical Procedure/$$
    110 to solve for $latex y_1 ( s_k , a )$$ for $latex k = 1 , 2 , 3, 4$$.
     112to solve for $latex y_1 ( s_i , a )$$ for $latex i = 1 , \ldots , Nz$$.
    111113
    112114$subhead Two levels of Iteration$$
     
    139141\begin{array}{rcl}
    140142{\rm minimize}
    141         & \sum_{k=1}^4 H_k( y^{k * ns} , a )
    142         & \; {\rm w.r.t} \; y^1 \in \R^2 , \ldots , y^{4 * ns} \in \R^2 ,
    143         \; a \in \R^3
     143& \sum_{i=1}^{Nz} H_i ( y^{N(i)} , a )
     144& \; {\rm w.r.t} \; y^1 \in \R^{Ny} , \ldots , y^{S(Nz)} \in \R^{Ny} ,
     145  \; a \in \R^{Na}
    144146\\
    145147{\rm subject \; to}
    146         & y^{M+1}  =  y^M +
    147 \left[ G( y^M , a ) + G( y^{M+1} , a ) \right] * \frac{ s_{M+1} - s_M }{ 2 }
    148         & \; {\rm for} \; M = 0 , \ldots , 4 * ns - 1
     148        & y^j  =  y^{j-1} +
     149\left[ G( y^{j-1} , a ) + G( y^j , a ) \right] * \frac{ t_j - t_{j-1} }{ 2 }
     150        & \; {\rm for} \; j = 1 , \ldots , S(Nz)
    149151\\
    150152        & y^0 = F(a)
    151153\end{array}
    152154\] $$
    153 where $latex ns$$ is the number of time intervals
    154 (used by the trapezoidal approximation)
    155 between each of the measurement times.
     155where for $latex i = 1, \ldots , Nz$$,
     156$latex N(i)$$ is the number of time intervals between
     157$latex s_{i-1}$$ and $latex s_i$$ (with $latex s_0 = 0$$)
     158and $latex S(i) = N(1) + \ldots + N(i)$$.
    156159Note that, in this form, the iterations of the optimization procedure
    157160also solve the forward problem equations.
     
    193196        standard deviation of measurement noise
    194197$rnext
    195 $latex e_k = 0$$ $pre $$ $cnext
    196         simulated measurement noise, $latex k = 1 , 2 , 3, 4$$
     198$latex e_i = 0$$ $pre $$ $cnext
     199        simulated measurement noise, $latex i = 1 , \ldots , Nz$$
    197200$rnext
    198 $latex s_k = k * .5$$ $pre $$ $cnext
    199         time corresponding to the $th k$$ measurement,
    200         $latex k = 1 , 2 , 3, 4$$
     201$latex s_i = i * .5$$ $pre $$ $cnext
     202        time corresponding to the $th i$$ measurement,
     203        $latex i = 1 , \ldots , Nz$$
    201204$tend
    202205
     
    205208$latex \[
    206209\begin{array}{rcl}
    207 z_k
    208 & = &  y_1 ( s_k , \bar{a} ) + e_k
     210z_i
     211& = &  y_1 ( s_i , \bar{a} ) + e_i
    209212\\
    210213& = &
    211214\bar{a}_0 * \bar{a}_1 *
    212         \frac{\exp( - \bar{a}_2 * s_k ) - \exp( -\bar{a}_1 * s_k )}
     215        \frac{\exp( - \bar{a}_2 * s_i ) - \exp( -\bar{a}_1 * s_i )}
    213216                { \bar{a}_1 - \bar{a}_2 }
    214217\end{array}
    215218\] $$
    216 for $latex k = 1, 2, 3, 4$$.
     219for $latex k = 1, \ldots , Nz$$.
    217220
    218221$end
    219222-----------------------------------------------------------------------------
    220 $begin ipopt_cppad_ode_represent$$
     223$begin ipopt_ode_simple$$
    221224$spell
    222225        ipopt_cppad_nlp
    223226$$
    224227
    225 $section ipopt_cppad_nlp ODE Problem Representation$$
    226 $index representation, ipopt_cppad_nlp ode$$
    227 $index ipopt_cppad_nlp, ode representation$$
    228 $index ode, ipopt_cppad_nlp representation$$
     228$section ODE Fitting Using Simple Representation$$
     229
     230$index ipopt_cppad_nlp, ode simple representation$$
     231$index ode, ipopt_cppad_nlp simple representation$$
     232$index simple, ipopt_cppad_nlp ode representation$$
    229233
    230234$head Purpose$$
    231235In this section we represent the objective and constraint functions,
    232236(in the simultaneous forward and reverse optimization problem)
    233 in terms of much simpler functions that are faster to differentiate
    234 (either by hand coding or by using AD).
     237using the $cref/simple representation/ipopt_cppad_nlp/Simple Representation/$$
     238in the sense of $code ipopt_cppad_nlp$$.
    235239
    236240$head Trapezoidal Time Grid$$
    237241The discrete time grid, used for the trapezoidal approximation, is
    238 denote by $latex \{ t_M \} $$.
    239 For $latex k = 1 , 2 , 3, 4$$,
    240 and $latex \ell = 0 , \ldots , ns$$, we define
    241 $latex \[
    242 \begin{array}{rcl}
    243         \Delta_k & = & ( s_k - s_{k-1} ) / ns
     242denoted by $latex \{ t_i \} $$ which is defined by:
     243$latex t_0 = 0$$ and
     244for $latex i = 1 , \ldots , Nz$$ and for $latex j = 1 , \ldots , N(i)$$,
     245$latex \[
     246\begin{array}{rcl}
     247        \Delta t_i & = & ( s_i - s_{i-1} ) / N(i)
    244248        \\
    245         t_{ (k-1) * ns + \ell } & = &  s_{k-1} + \Delta_k * \ell
    246 \end{array}
    247 \] $$
    248 Note that for $latex M = 1 , \ldots , 4 * ns $$,
    249 $latex y^M$$ denotes our approximation for $latex y( t_M , a )$$,
    250 $latex t_0$$ is equal to 0,
    251 and $latex t_{k*ns}$$ is equal to $latex s_k$$.
     249        t_{S(i-1)+j} & = & s_{i-1} + \Delta t_i * j
     250\end{array}
     251\] $$
     252where $latex s_0 = 0$$,
     253$latex N(i)$$ is the number of time grid points between
     254$latex s_{i-1}$$ and $latex s_i$$,
     255$latex S(0) = 0$$, and
     256$latex S(i) = N(1) + \ldots + N(i)$$.
     257Note that for $latex i = 0 , \ldots , S(Nz)  $$,
     258$latex y^i$$ denotes our approximation for $latex y( t_i , a )$$
     259and $latex t_{S(i)}$$ is equal to $latex s_i$$.
    252260
    253261$head Argument Vector$$
     
    256264has the following structure
    257265$latex \[
    258         x = \left( \begin{array}{c}
    259                 y^0 \\ \vdots \\ y^{4 * ns} \\ a
    260         \end{array} \right)
     266x = ( y^0 , \cdots , y^{S(Nz)} , a )
    261267\] $$ 
    262 Note that $latex x \in \R^{2 *(4 * ns + 1) + 3}$$ and
    263 $latex \[
    264 \begin{array}{rcl}
    265         y^M & = & ( x_{2 * M} , x_{2 * M + 1} )^\T
     268Note that $latex x \in \R^{S(Nz) + Na}$$ and
     269$latex \[
     270\begin{array}{rcl}
     271        y^i & = & ( x_{Ny * i} , \ldots ,  x_{Ny * i + Ny - 1} )
    266272        \\
    267         a   & = & ( x_{8 * ns + 2} , x_{8 * ns + 3} , x_{8 * ns + 4} )^\T
    268 \end{array}
    269 \] $$
    270 
    271 $head Objective$$
     273        a   & = & ( x_{Ny *S(Nz) + Ny} , \ldots , x_{Ny * S(Nz) + Na - 1} )
     274\end{array}
     275\] $$
     276
     277$head Objective Function$$
    272278The objective function
    273279( $latex fg_0 (x)$$ in $cref/ipopt_cppad_nlp/$$ )
    274280has the following representation,
    275281$latex \[
    276         fg_0 (x)
    277         = \sum_{k=1}^4 H_k ( y^{k * ns} , a )
    278         = \sum_{k=0}^3 r^k ( u^{k,0} )
    279 \] $$
    280 where $latex r^k = H_{k+1}$$ and $latex u^{k,0} =   ( y^{k * ns} , a )$$
    281 for $latex k = 0, 1, 2, 3$$.
    282 The range index (in the vector $latex fg (x)$$ )
    283 corresponding to each term in the objective is 0.
    284 The domain components (in the vector $latex x$$)
    285 corresponding to the $th k$$ term are
    286 $latex \[
    287 (       x_{2 * k * ns} ,
    288         x_{2 * k * ns + 1} ,
    289         x_{8 * ns + 2} ,
    290         x_{8 * ns + 3} ,
    291         x_{8 * ns + 4}
    292 )
    293 = u^{k,0}
    294 = ( y_0^{k * ns} , y_1^{k * ns} , a_0, a_1, a_2 )
    295 \] $$
    296 
    297 $head Initial Condition$$
    298 We define the function
    299 $latex r^k : \R^2 \times \R^3 \rightarrow \R^2$$ for $latex k = 4$$ by
    300 $latex \[
    301         r^k ( u ) = r^k ( y , a ) = y - F(a)
    302 \] $$
    303 where $latex u  = ( y , a)$$.
    304 Using this notation,
    305 and the function $latex fg (x)$$ in $cref/ipopt_cppad_nlp/$$,
    306 the constraint becomes
    307 $latex \[
    308 \begin{array}{rcl}
    309         fg_1 (x) & = & r_0^4 ( u^{4,0} ) = r_0^4 ( y^0 , a)
    310         \\
    311         fg_2 (x) & = & r_1^4 ( u^{4,0} ) = r_1^4 ( y^0 , a)
    312 \end{array}
    313 \] $$
    314 The range indices (in the vector $latex fg (x)$$ )
    315 corresponding to this constraint are $latex ( 1, 2 )$$.
    316 The domain components (in the vector $latex x$$)
    317 corresponding to this constraint are
    318 $latex \[
    319 (       x_0 ,
    320         x_1 ,
    321         x_{8 * ns + 2} ,
    322         x_{8 * ns + 3} ,
    323         x_{8 * ns + 4}
    324 )
    325 =
    326 u^{4,0}
    327 =
    328 ( y_0^0, y_1^0 , \ldots , y_0^{4*ns}, y_1^{4*ns}, a_0 , a_1 , a_2 )
    329 \] $$
    330 
    331 $head Trapezoidal Approximation$$
    332 For $latex k = 5 , 6, 7 , 8$$,
    333 and $latex \ell = 0 , \ldots , ns$$,
    334 define $latex M = (k - 5) * ns + \ell$$.
    335 The corresponding trapezoidal approximation is represented by the constraint
    336 $latex \[
    337 0 = y^{M+1}  -  y^{M} -
    338 \left[ G( y^{M} , a ) + G( y^{M+1} , a ) \right] * \frac{\Delta_k }{ 2 }
    339 \] $$
    340 For $latex k = 5, 6, 7, 8$$, we define the function
    341 $latex r^k : \R^2 \times \R^2 \times \R^3 \rightarrow \R^2$$ by
    342 $latex \[
    343 r^k ( y , w , a ) = y - w  [ G( y , a ) + G( w , a ) ] * \frac{ \Delta_k }{ 2 }
    344 \] $$
    345 Using this notation,
    346 (and the function $latex fg (x) $$ in $cref/ipopt_cppad_nlp/$$ )
    347 the constraint becomes
    348 $latex \[
    349 \begin{array}{rcl}
    350 fg_{2 * M + 3} (x)  & = & r_0 ( u^{k,\ell} ) = r_0^k ( y^M , y^{M+1} , a )
    351 \\
    352 fg_{2 * M + 4} (x)  & = & r_1 ( u^{k,\ell} ) = r_1^k ( y^M , y^{M+1} , a )
    353 \end{array}
    354 \] $$
    355 where $latex M = (k - 5) * ns * \ell$$.
    356 The range indices (in the vector $latex fg (x)$$ )
    357 corresponding to this constraint are
    358 $latex ( 2 * M + 3 , 2 * M + 4 )$$.
    359 The domain components (in the vector $latex x$$)
    360 corresponding to this constraint are
    361 $latex \[
    362 (       x_{2 * M} ,
    363         x_{2 * M + 1} ,
    364         x_{2 * M + 2} ,
    365         x_{2 * M + 3} ,
    366         x_{8 * ns + 2} ,
    367         x_{8 * ns + 3} ,
    368         x_{8 * ns + 4}
    369 )
    370 = u^{k, \ell}
    371 = ( y_0^M, y_1^M , y_0^{M+1} , y_1^{M+1} , a_0 , a_1 , a_2 )
    372 \] $$
     282        fg_0 (x) = \sum_{i=1}^{Nz} H_i ( y^{S(i)} , a )
     283\] $$
     284
     285$head Initial Condition Constraint$$
     286For $latex i = 1 , \ldots , Ny$$,
     287we define the component functions $latex fg_i (x)$$,
     288and corresponding constraint equations, by
     289$latex \[
     290        0 = fg_i ( x ) = y_i^0 - F_i (a)
     291\] $$
     292
     293$head Trapezoidal Approximation Constraint$$
     294For $latex i = 1, \ldots , S(Nz)$$,
     295and for $latex j = 1 , \ldots , Ny$$,
     296we define the component functions $latex fg_{Ny*i + j} (x)$$,
     297and corresponding constraint equations, by
     298$latex \[
     2990 = fg_{Ny*i + j } = y_j^{i}  -  y_j^{i-1} -
     300        \left[ G_j ( y^i , a ) + G_j ( y^{i-1} , a ) \right] *
     301                \frac{t_i - t_{i-1} }{ 2 }
     302\] $$
     303
     304$children%
     305        ipopt_cppad/ipopt_ode_simple.cpp
     306%$$
     307$head Source$$
     308The file $cref ipopt_ode_simple.cpp$$
     309contains the source code for this example.
    373310
    374311$end
    375312-----------------------------------------------------------------------------
    376 $begin ipopt_ode_simple.cpp$$
     313$begin ipopt_ode_fast$$
    377314$spell
    378315        ipopt_cppad_nlp
    379         nd
    380         ny
    381         na
    382         ns
    383316$$
    384317
    385 $section ipopt_cppad_nlp ODE Example Source Code$$
    386 $index ipopt_cppad_nlp, ode example source$$
    387 $index ode, ipopt_cppad_nlp example source$$
    388 $index example, ipopt_cppad_nlp ode source$$
    389 $index source, ipopt_cppad_nlp ode example$$
    390 
    391 $head Source Code$$
    392 This example uses the a simple representation of the objective
    393 an constraints to solve the ODE fitting problem.
    394 Almost all the code below is for the general problem
    395 (where $code nd$$, $code ny$$, $code na$$, and $code ns$$ are arbitrary)
    396 but some of it for a specific case defined by the function $code y_one(t)$$
    397 and discussed in the previous sections.
    398 $code
    399 $verbatim%ipopt_cppad/ipopt_ode_simple.cpp%0%$$
    400 $$
    401 
    402 $end
    403 -----------------------------------------------------------------------------
    404 $begin ipopt_ode_fast.cpp$$
    405 $spell
    406         ipopt_cppad_nlp
    407         nd
    408         ny
    409         na
    410         ns
    411 $$
    412 
    413 $section ipopt_cppad_nlp ODE Example Source Code$$
    414 $index ipopt_cppad_nlp, ode example source$$
    415 $index ode, ipopt_cppad_nlp example source$$
    416 $index example, ipopt_cppad_nlp ode source$$
    417 $index source, ipopt_cppad_nlp ode example$$
    418 
    419 $head Source Code$$
    420 This example uses the a decomposed representation of the objective
    421 an constraints to solve the ODE fitting problem.
    422 Almost all the code below is for the general problem
    423 (where $code nd$$, $code ny$$, $code na$$, and $code ns$$ are arbitrary)
    424 but some of it for a specific case defined by the function $code y_one(t)$$
    425 and discussed in the previous sections.
    426 $code
    427 $verbatim%ipopt_cppad/ipopt_ode_fast.cpp%0%$$
    428 $$
    429 
    430 $end
     318$section ODE Fitting Using Fast Representation$$
     319
     320$index representation, ipopt_cppad_nlp ode$$
     321$index ipopt_cppad_nlp, ode representation$$
     322$index ode, ipopt_cppad_nlp representation$$
     323
     324$head Purpose$$
     325In this section we represent a more complex representation of the
     326simultaneous forward and reverse ODE fitting problem (described above).
     327The representation defines the problem using
     328simpler functions that are faster to differentiate
     329(either by hand coding or by using AD).
     330
     331$head Objective Function$$
     332We use the following representation for the
     333$cref/objective function/ipopt_ode_simple/Objective Function/$$:
     334For $latex k = 0 , \ldots , Nz - 1$$,
     335we define the function $latex r^k : \R^{Ny+Na} \rightarrow \R$$
     336by
     337$latex \[
     338\begin{array}{rcl}
     339fg_0 (x) & = & \sum_{i=1}^{Nz} H_i ( y^{S(i)} , a )
     340\\
     341fg_0 (x) & = & \sum_{k=0}^{Nz-1} r^k ( u^{k,0} )
     342\end{array}
     343\] $$
     344where for $latex k = 0 , \ldots , Nz-1$$,
     345$latex u^{k,0} \in \R^{Ny + Na}$$ is defined by
     346$latex u^{k,0} =   ( y^{S(k+1)} , a )$$
     347
     348$subhead Range Indices I(k,0)$$
     349For $latex k = 0 , \ldots , Nz - 1$$,
     350the range index in the vector $latex fg (x)$$
     351corresponding to $latex r^k ( u^{k,0} ) $$ is 0.
     352Thus, the range indices are given by
     353$latex I(k,0) = \{ 0 \}$$ for $latex k = 0 , \ldots , Nz-1$$.
     354
     355$subhead Domain Indices J(k,0)$$
     356For $latex k = 0 , \ldots , Nz - 1$$,
     357the components of the vector $latex x$$
     358corresponding to the vector $latex u^{k,0}$$ are
     359$latex \[
     360\begin{array}{rcl}
     361u^{k,0} & = & ( y^{S(k+1} , a )
     362\\
     363& = &
     364(       x_{Ny * S(k+1)} \; , \;
     365        \ldots \; , \;
     366        x_{Ny * S(k+1) + Ny - 1} \; , \;
     367        x_{Ny * S(Nz) + Ny } \; , \;
     368        \ldots \; , \;
     369        x_{Ny * S(Nz) + Ny + Na - 1}
     370)
     371\end{array}
     372\] $$
     373Thus, the domain indices are given by
     374$latex \[
     375J(k,0) = \{
     376        Ny * S(k+1) \; , \;
     377        \ldots  \; , \;
     378        Ny * S(k+1) + Ny - 1 \; , \;
     379        Ny * S(Nz) + Ny \; , \;
     380        \ldots  \; , \;
     381        Ny * S(Nz) + Ny + Na - 1
     382\}
     383\] $$
     384
     385$head Initial Condition$$
     386We use the following representation for the
     387$cref/initial condition constraint/
     388        ipopt_ode_simple/Initial Condition Constraint/$$:
     389For $latex k = Nz$$ we define the function
     390$latex r^k : \R^{Ny} \times \R^{Na + Ny}$$ by
     391$latex \[
     392\begin{array}{rcl}
     393        0 & = & fg_i ( x ) = y_i^0 - F_i (a)
     394        \\
     395        0 & = & r_{i-1}^k ( u^{k,0} ) = y_i^0 - F_i(a)
     396\end{array}
     397\] $$
     398where $latex i = 1 , \ldots , Ny$$ and
     399where $latex u^{k,0} \in \R^{Ny + Na}$$ is defined by
     400$latex u^{k,0}  = ( y^0 , a)$$.
     401
     402$subhead Range Indices I(k,0)$$
     403For $latex k = Nz$$,
     404the range index in the vector $latex fg (x)$$
     405corresponding to $latex r^k ( u^{k,0} ) $$ are
     406$latex I(k,0) = \{ 1 , \ldots , Ny \}$$.
     407
     408$subhead Domain Indices J(k,0)$$
     409For $latex k = Nz$$,
     410the components of the vector $latex x$$
     411corresponding to the vector $latex u^{k,0}$$ are
     412$latex \[
     413\begin{array}{rcl}
     414u^{k,0} & = & ( y^0 , a)
     415\\
     416& = &
     417(       x_0 \; , \;
     418        \ldots \; , \;
     419        x_{Ny-1} \; , \;
     420        x_{Ny * S(Nz) + Ny } \; , \;
     421        \ldots \; , \;
     422        x_{Ny * S(Nz) + Ny + Na - 1}
     423)
     424\end{array}
     425\] $$
     426Thus, the domain indices are given by
     427$latex \[
     428J(k,0) = \{
     429        0 \; , \;
     430        \ldots  \; , \;
     431        Ny - 1 \; , \;
     432        Ny * S(Nz) + Ny \; , \;
     433        \ldots  \; , \;
     434        Ny * S(Nz) + Ny + Na - 1
     435\}
     436\] $$
     437
     438$head Trapezoidal Approximation$$
     439We use the following representation for the
     440$cref/trapezoidal approximation constraint/
     441        ipopt_ode_simple/Trapezoidal Approximation Constraint/$$:
     442For $latex k = 1 , \ldots , Nz$$,
     443we define the function $latex r^{Nz+k} : \R^{2*Ny+Na} \rightarrow \R^{Ny}$$ by
     444$latex \[
     445r^{Nz+k} ( y , w , a )
     446=
     447y - w  [ G( y , a ) + G( w , a ) ] * \frac{ \Delta t_k }{ 2 }
     448\] $$
     449For $latex \ell = 0 , \ldots , N(k)-1$$,
     450using the notation $latex i = Ny * S(k-1) + \ell + 1$$,
     451the corresponding trapezoidal approximation is represented by
     452$latex \[
     453\begin{array}{rcl}
     4540 & = & fg_{Ny+i} (x)
     455\\
     456& = &
     457y^i  -  y^{i-1} -
     458\left[ G( y^i , a ) + G( y^{i-1} , a ) \right] * \frac{\Delta t_k }{ 2 }
     459\\
     460& = &
     461r^{Nz+k} ( u^{Nz+k , \ell} )
     462\end{array}
     463\] $$
     464where $latex u^{Nz+k,\ell} \in \R^{2*Ny + Na}$$ is defined by
     465$latex u^{Nz+k,\ell}  = ( y^{i-1} , y^i , a)$$.
     466
     467$subhead Range Indices I(k,0)$$
     468For $latex k = Nz + 1, \ldots , 2*Nz$$,
     469and $latex \ell = 0 , \ldots , N(k)-1$$,
     470the range index in the vector $latex fg (x)$$
     471corresponding to $latex r^k ( u^{k,\ell} ) $$ are
     472$latex I(k,\ell) = \{ Ny + i , \ldots , 2*Ny + i - 1  \}$$
     473where $latex i = Ny * S(k-1) + \ell + 1$$.
     474
     475$subhead Domain Indices J(k,0)$$
     476For $latex k = Nz + 1, \ldots , 2*Nz$$,
     477and $latex \ell = 0 , \ldots , N(k)-1$$,
     478define $latex i = Ny * S(k-1) + \ell + 1$$.
     479The components of the vector $latex x$$
     480corresponding to the vector $latex u^{k,\ell}$$ are
     481(and the function $latex fg (x) $$ in $cref/ipopt_cppad_nlp/$$ )
     482$latex \[
     483\begin{array}{rcl}
     484u^{k, \ell} & = & ( y^{i-1} , y^i , a )
     485\\
     486& = &
     487(       x_{Ny * (i-1)} \; , \;
     488        \ldots \; , \;
     489        x_{Ny * (i+1) - 1} \; , \;
     490        x_{Ny * S(Nz) + Ny } \; , \;
     491        \ldots \; , \;
     492        x_{Ny * S(Nz) + Ny + Na - 1}
     493)
     494\end{array}
     495\] $$
     496Thus, the domain indices are given by
     497$latex \[
     498J(k,\ell) = \{
     499        Ny * (i-1) \; , \;
     500        \ldots  \; , \;
     501        Ny * (i+1) - 1 \; , \;
     502        Ny * S(Nz) + Ny \; , \;
     503        \ldots  \; , \;
     504        Ny * S(Nz) + Ny + Na - 1
     505\}
     506\] $$
     507
     508$children%
     509        ipopt_cppad/ipopt_ode_fast.cpp
     510%$$
     511$head Source$$
     512The file $cref ipopt_ode_fast.cpp$$
     513contains the source code for this example.
     514
     515 
     516$end
    431517------------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.