Changeset 1277 for trunk/test_more


Ignore:
Timestamp:
Sep 9, 2008 8:57:06 AM (12 years ago)
Author:
bradbell
Message:

trunk: Fix indexing error in ipopt_cppad_nlp.cpp.

ipopt_cppad_nlp.cpp: error was here.
gpl_license.sh: fix error message when the script ran.
ipopt_cppad.cpp: add more ipopt_cppad_nlp tests to automated suite.
whats_new_08.omh: user's view of the changes.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/test_more/ipopt_cppad.cpp

    r1276 r1277  
    1212# include "../example/ipopt_cppad_nlp.hpp"
    1313
    14 namespace {
    15 
    16         class my_FG_info : public ipopt_cppad_fg_info
    17         {
    18         public:
    19                 // derived class part of constructor
    20                 my_FG_info(void)
    21                 { }
    22                 // Evaluation of the objective f(x), and constraints g(x)
    23                 // using an Algorithmic Differentiation (AD) class.
    24                 ADVector r_eval(size_t k, const ADVector&  x)
    25                 {       ADVector fg(2);
    26 
    27                         // f(x)
    28                         if( x[0] >= 1. )
    29                                 fg[0] = .5 * (x[0] * x[0] + x[1] * x[1]);
    30                         else    fg[0] = x[0] + .5 * x[1] * x[1];
    31                         // g (x)
    32                         fg[1] = x[0];
    33 
    34                         return fg;
    35                 }
    36                 bool retape(size_t k)
    37                 {       return true; }
    38         };
    39 }
    40        
    41 bool ipopt_cppad(void)
     14namespace { // Begin empty namespace
     15
     16class FG_retape : public ipopt_cppad_fg_info
     17{
     18public:
     19        // derived class part of constructor
     20        FG_retape(void)
     21        { }
     22        // Evaluation of the objective f(x), and constraints g(x)
     23        // using an Algorithmic Differentiation (AD) class.
     24        ADVector r_eval(size_t k, const ADVector&  x)
     25        {       ADVector fg(2);
     26
     27                // f(x)
     28                if( x[0] >= 1. )
     29                        fg[0] = .5 * (x[0] * x[0] + x[1] * x[1]);
     30                else    fg[0] = x[0] + .5 * x[1] * x[1];
     31                // g (x)
     32                fg[1] = x[0];
     33
     34                return fg;
     35        }
     36        bool retape(size_t k)
     37        {       return true; }
     38};
     39
     40bool ipopt_cppad_retape(void)
    4241{       bool ok = true;
    4342        size_t j;
     
    6463
    6564        // object in derived class
    66         my_FG_info my_fg_info;
    67         ipopt_cppad_fg_info *fg_info = &my_fg_info
     65        FG_retape fg_retape;
     66        ipopt_cppad_fg_info *fg_info = &fg_retape
    6867
    6968        // create the Ipopt interface
     
    116115        return ok;
    117116}
     117
     118/*
     119This solve the same problem as ../example/ipopt_cppad.cpp (repository revision 1276) in a convoluted way in order to test the representation code.
     120*/
     121class FG_K_gt_1 : public ipopt_cppad_fg_info
     122{
     123private:
     124        bool retape_;
     125public:
     126        // derived class part of constructor
     127        FG_K_gt_1(bool retape)
     128        : retape_ (retape)
     129        { }
     130        // Evaluation of the objective f(x), and constraints g(x)
     131        // using an Algorithmic Differentiation (AD) class.
     132        ADVector r_eval(size_t k, const ADVector&  u)
     133        {
     134
     135                // Fortran style indexing
     136                ADNumber x1 = u[3];
     137                ADNumber x2 = u[2];
     138                ADNumber x3 = u[1];
     139                ADNumber x4 = u[0];
     140                if( k == 0 )
     141                {       ADVector r(1);
     142                        // f(x)
     143                        r[0] = x1 * x4 * (x1 + x2 + x3) + x3;
     144                        return r;
     145                }
     146                ADVector r(2);
     147                // g_1 (x)
     148                r[0] = x1 * x2 * x3 * x4;
     149                // g_2 (x)
     150                r[1] = x1 * x1 + x2 * x2 + x3 * x3 + x4 * x4;
     151                return r;
     152        }
     153        bool retape(size_t k)
     154        {       return retape_; }
     155        size_t number_functions(void)
     156        {       return 2; }
     157        size_t domain_size(size_t k)
     158        {       return 4; }
     159        size_t range_size(size_t k)
     160        {       if( k == 0 )
     161                        return 1;
     162                return 2;
     163        }
     164        size_t number_terms(size_t k)
     165        {       return 1; }
     166        void domain_index(size_t k, size_t ell, SizeVector& J)
     167        {       size_t j;
     168                // reversion the order of the variables in u from that in x
     169                for(j = 0; j < 4; j++)
     170                        J[j] = 3-j;
     171        }
     172        void range_index(size_t k, size_t ell, SizeVector& I)
     173        {       if( k == 0 )
     174                        I[0] = 0;
     175                else
     176                {       I[0] = 1;
     177                        I[1] = 2;
     178                }
     179        }
     180};
     181       
     182bool ipopt_cppad_K_gt_1(void)
     183{       bool ok = true;
     184        size_t j;
     185
     186
     187        // number of independent variables (domain dimension for f and g)
     188        size_t n = 4; 
     189        // number of constraints (range dimension for g)
     190        size_t m = 2;
     191        // initial value of the independent variables
     192        NumberVector x_i(n);
     193        x_i[0] = 1.0;
     194        x_i[1] = 5.0;
     195        x_i[2] = 5.0;
     196        x_i[3] = 1.0;
     197        // lower and upper limits for x
     198        NumberVector x_l(n);
     199        NumberVector x_u(n);
     200        for(j = 0; j < n; j++)
     201        {       x_l[j] = 1.0;
     202                x_u[j] = 5.0;
     203        }
     204        // lower and upper limits for g
     205        NumberVector g_l(m);
     206        NumberVector g_u(m);
     207        g_l[0] = 25.0;     g_u[0] = 1.0e19;
     208        g_l[1] = 40.0;     g_u[1] = 40.0;
     209
     210        size_t icase;
     211        for(icase = 0; icase <= 1; icase++)
     212        {       // Should ipopt_cppad_nlp retape the operation sequence for
     213                // every new x. Can test both true and false cases because
     214                // the operation sequence does not depend on x (for this case).
     215                bool retape = bool(icase);
     216
     217                // object in derived class
     218                FG_K_gt_1 my_fg_info(retape);
     219                ipopt_cppad_fg_info *fg_info = &my_fg_info; 
     220
     221                // create the Ipopt interface
     222                ipopt_cppad_solution solution;
     223                Ipopt::SmartPtr<Ipopt::TNLP> cppad_nlp = new ipopt_cppad_nlp(
     224                n, m, x_i, x_l, x_u, g_l, g_u, fg_info, &solution
     225                );
     226
     227                // Create an instance of the IpoptApplication
     228                using Ipopt::IpoptApplication;
     229                Ipopt::SmartPtr<IpoptApplication> app = new IpoptApplication();
     230
     231                // turn off any printing
     232                app->Options()->SetIntegerValue("print_level", -2);
     233
     234                // maximum number of iterations
     235                app->Options()->SetIntegerValue("max_iter", 10);
     236
     237                // approximate accuracy in first order necessary conditions;
     238                // see Mathematical Programming, Volume 106, Number 1,
     239                // Pages 25-57, Equation (6)
     240                app->Options()->SetNumericValue("tol", 1e-9);
     241
     242                // derivative testing
     243                app->Options()->
     244                SetStringValue("derivative_test", "second-order");
     245
     246                // Initialize the IpoptApplication and process the options
     247                Ipopt::ApplicationReturnStatus status = app->Initialize();
     248                ok    &= status == Ipopt::Solve_Succeeded;
     249
     250                // Run the IpoptApplication
     251                status = app->OptimizeTNLP(cppad_nlp);
     252                ok    &= status == Ipopt::Solve_Succeeded;
     253
     254                /*
     255                Check some of the solution values
     256                */
     257                ok &= solution.status == ipopt_cppad_solution::success;
     258                //
     259                double check_x[]   = { 1.000000, 4.743000, 3.82115, 1.379408 };
     260                double check_z_l[] = { 1.087871, 0.,       0.,      0.       };
     261                double check_z_u[] = { 0.,       0.,       0.,      0.       };
     262                double rel_tol     = 1e-6;  // relative tolerance
     263                double abs_tol     = 1e-6;  // absolute tolerance
     264                for(j = 0; j < n; j++)
     265                {       ok &= CppAD::NearEqual(
     266                        check_x[j],   solution.x[j],   rel_tol, abs_tol
     267                        );
     268                        ok &= CppAD::NearEqual(
     269                        check_z_l[j], solution.z_l[j], rel_tol, abs_tol
     270                        );
     271                        ok &= CppAD::NearEqual(
     272                        check_z_u[j], solution.z_u[j], rel_tol, abs_tol
     273                        );
     274                }
     275        }
     276
     277        return ok;
     278}
     279
     280
     281/*
     282f(x)    = x[1]; k=0, ell=0, I[0] = 0, J[0] = 1
     283g_0 (x) = x[0]; k=0, ell=1, I[0] = 1, J[0] = 0
     284g_1 (x) = x[1]; k=0, ell=2, I[0] = 2, J[0] = 1
     285
     286minimize   f(x)
     287subject to -1 <= g_0(x)  <= 0
     288            0 <= g_1 (x) <= 1
     289
     290The solution is x[1] = 0 and x[0] arbitrary.
     291*/
     292
     293namespace
     294{
     295class FG_J_changes : public ipopt_cppad_fg_info
     296{
     297private:
     298        bool retape_;
     299public:
     300        // constructor
     301        FG_J_changes(bool retape)
     302        : retape_ (retape)
     303        { }
     304        size_t number_functions(void)
     305        {       return 1; }
     306        size_t domain_size(size_t k)
     307        {       size_t q;
     308                switch(k)
     309                {       case 0:  q = 1;  break; 
     310                        default: assert(0);
     311                }
     312                return q;
     313        }
     314        size_t range_size(size_t k)
     315        {       size_t p;
     316                switch(k)
     317                {       case 0:  p = 1;  break;
     318                        default: assert(0);
     319                }
     320                return p;
     321        }
     322        size_t number_terms(size_t k)
     323        {       size_t L;
     324                switch(k)
     325                {       case 0:  L = 3;   break;
     326                        default: assert(0);
     327                }
     328                return L;
     329        }
     330        void domain_index(size_t k, size_t ell, SizeVector& J)
     331        {       assert( J.size() >= 1 );
     332                if( ell == 0 )
     333                {       J[0] = 1;
     334                        return;
     335                }
     336                J[0] = ell - 1;
     337                return;
     338        }
     339        void range_index(size_t k, size_t ell, SizeVector& I)
     340        {       assert( I.size() >= 1 );
     341                I[0] = ell;
     342                return;
     343        }
     344        // retape function
     345        bool retape(size_t k)
     346        {       return retape_; }
     347        ADVector r_eval(size_t k, const ADVector&  u)
     348        {
     349                assert( u.size() == 1 );
     350                ADVector r(1);
     351                r[0] = u[0] ;
     352                return r;
     353        }
     354};
     355}
     356
     357bool ipopt_cppad_J_changes(void)
     358{
     359        bool ok = true;
     360        // number of independent variables (domain dimension for f and g)
     361        size_t n = 2;
     362        // number of constraints (range dimension for g)
     363        size_t m = 2;
     364        // initial value of the independent variables
     365        NumberVector x_i(n);
     366        NumberVector x_l(n);
     367        NumberVector x_u(n);
     368
     369        size_t i = 0;
     370        for(i = 0; i < n; i++)
     371        {       x_i[i] = 0.;
     372                x_l[i] = -1.0;
     373                x_u[i] = +1.0;
     374        }
     375
     376        // lower and upper limits for g
     377        NumberVector g_l(m);
     378        NumberVector g_u(m);
     379        g_l[0] = -1; g_u[0] = 0.;
     380        g_l[1] = 0.; g_u[1] = 1.;
     381
     382        // object for evaluating function
     383        bool retape = false;
     384        FG_J_changes my_fg_info(retape);
     385        ipopt_cppad_fg_info *fg_info = &my_fg_info;
     386
     387        ipopt_cppad_solution solution;
     388        Ipopt::SmartPtr<Ipopt::TNLP> cppad_nlp = new ipopt_cppad_nlp(
     389                n, m, x_i, x_l, x_u, g_l, g_u, fg_info, &solution
     390        );
     391
     392        // Create an instance of the IpoptApplication
     393        using Ipopt::IpoptApplication;
     394        Ipopt::SmartPtr<IpoptApplication> app = new IpoptApplication();
     395
     396        // turn off any printing
     397        app->Options()->SetIntegerValue("print_level", -2);
     398
     399        // approximate accuracy in first order necessary conditions;
     400        // see Mathematical Programming, Volume 106, Number 1,
     401        // Pages 25-57, Equation (6)
     402        app->Options()->SetNumericValue("tol", 1e-9);
     403        app->Options()-> SetStringValue("derivative_test", "second-order");
     404
     405        // Initialize the IpoptApplication and process the options
     406        Ipopt::ApplicationReturnStatus status = app->Initialize();
     407        ok    &= status == Ipopt::Solve_Succeeded;
     408
     409        // Run the IpoptApplication
     410        status = app->OptimizeTNLP(cppad_nlp);
     411        ok    &= status == Ipopt::Solve_Succeeded;
     412
     413        /*
     414         Check solution status
     415         */
     416        ok &= solution.status == ipopt_cppad_solution::success;
     417        ok &= CppAD::NearEqual(solution.x[1], 0., 1e-6, 1e-6);
     418
     419        return ok;
     420}
     421
     422
     423
     424} // End empty namespace
     425
     426bool ipopt_cppad(void)
     427{       bool ok = true;
     428        ok &= ipopt_cppad_retape();
     429        ok &= ipopt_cppad_K_gt_1();
     430        ok &= ipopt_cppad_J_changes();
     431        return ok;
     432}
Note: See TracChangeset for help on using the changeset viewer.