source: trunk/test_more/ipopt_cppad.cpp @ 1556

Last change on this file since 1556 was 1401, checked in by bradbell, 11 years ago

trunk: Merge in changes from branches/sweep.

  • Property svn:keywords set to Id
File size: 10.8 KB
Line 
1/* $Id: ipopt_cppad.cpp 1401 2009-06-24 17:56:51Z bradbell $ */
2/* --------------------------------------------------------------------------
3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-09 Bradley M. Bell
4
5CppAD is distributed under multiple licenses. This distribution is under
6the terms of the
7                    Common Public License Version 1.0.
8
9A copy of this license is included in the COPYING file of this distribution.
10Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
11-------------------------------------------------------------------------- */
12
13# include "../ipopt_cppad/ipopt_cppad_nlp.hpp"
14
15namespace { // Begin empty namespace
16// ---------------------------------------------------------------------------
17class FG_retape : public ipopt_cppad_fg_info
18{
19public:
20        // derived class part of constructor
21        FG_retape(void)
22        { }
23        // Evaluation of the objective f(x), and constraints g(x)
24        // using an Algorithmic Differentiation (AD) class.
25        ADVector eval_r(size_t k, const ADVector&  x)
26        {       ADVector fg(2);
27
28                // f(x)
29                if( x[0] >= 1. )
30                        fg[0] = .5 * (x[0] * x[0] + x[1] * x[1]);
31                else    fg[0] = x[0] + .5 * x[1] * x[1]; 
32                // g (x)
33                fg[1] = x[0];
34
35                return fg;
36        }
37        bool retape(size_t k)
38        {       return true; }
39}; 
40
41bool ipopt_cppad_retape(void)
42{       bool ok = true;
43        size_t j;
44
45
46        // number of independent variables (domain dimension for f and g)
47        size_t n = 2; 
48        // number of constraints (range dimension for g)
49        size_t m = 1;
50        // initial value of the independent variables
51        NumberVector x_i(n);
52        x_i[0] = 2.0;
53        x_i[1] = 2.0;
54        // lower and upper limits for x
55        NumberVector x_l(n);
56        NumberVector x_u(n);
57        for(j = 0; j < n; j++)
58        {       x_l[j] = -10.; x_u[j] = +10.;
59        }
60        // lower and upper limits for g
61        NumberVector g_l(m);
62        NumberVector g_u(m);
63        g_l[0] = -1.;     g_u[0] = 1.0e19;
64
65        // object in derived class
66        FG_retape fg_retape;
67        ipopt_cppad_fg_info *fg_info = &fg_retape; 
68
69        // create the Ipopt interface
70        ipopt_cppad_solution solution;
71        Ipopt::SmartPtr<Ipopt::TNLP> cppad_nlp = new ipopt_cppad_nlp(
72                n, m, x_i, x_l, x_u, g_l, g_u, fg_info, &solution
73        );
74
75        // Create an instance of the IpoptApplication
76        using Ipopt::IpoptApplication;
77        Ipopt::SmartPtr<IpoptApplication> app = new IpoptApplication();
78
79        // turn off any printing
80        app->Options()->SetIntegerValue("print_level", 0);
81
82        // maximum number of iterations
83        app->Options()->SetIntegerValue("max_iter", 10);
84
85        // approximate accuracy in first order necessary conditions;
86        // see Mathematical Programming, Volume 106, Number 1,
87        // Pages 25-57, Equation (6)
88        app->Options()->SetNumericValue("tol", 1e-9);
89
90        // derivative testing
91        app->Options()->
92        SetStringValue("derivative_test", "second-order");
93
94        // Initialize the IpoptApplication and process the options
95        Ipopt::ApplicationReturnStatus status = app->Initialize();
96        ok    &= status == Ipopt::Solve_Succeeded;
97
98        // Run the IpoptApplication
99        status = app->OptimizeTNLP(cppad_nlp);
100        ok    &= status == Ipopt::Solve_Succeeded;
101
102        /*
103        Check some of the solution values
104        */
105        ok &= solution.status == ipopt_cppad_solution::success;
106        //
107        double check_x[]   = { -1., 0. };
108        double rel_tol     = 1e-6;  // relative tolerance
109        double abs_tol     = 1e-6;  // absolute tolerance
110        for(j = 0; j < n; j++)
111        {       ok &= CppAD::NearEqual(
112                        check_x[j],   solution.x[j],   rel_tol, abs_tol
113                );
114        }
115
116        return ok;
117}
118// ---------------------------------------------------------------------------
119/*
120This solve the same problem as
121../ipopt_cppad/ipopt_cppad_simple.cpp (repository revision
1221276) in a convoluted way in order to test the representation code.
123*/
124class FG_K_gt_1 : public ipopt_cppad_fg_info
125{
126private:
127        bool retape_;
128public:
129        // derived class part of constructor
130        FG_K_gt_1(bool retape)
131        : retape_ (retape)
132        { }
133        // Evaluation of the objective f(x), and constraints g(x)
134        // using an Algorithmic Differentiation (AD) class.
135        ADVector eval_r(size_t k, const ADVector&  u)
136        {
137
138                // Fortran style indexing
139                ADNumber x1 = u[3];
140                ADNumber x2 = u[2];
141                ADNumber x3 = u[1];
142                ADNumber x4 = u[0];
143                if( k == 0 )
144                {       ADVector r(1);
145                        // f(x)
146                        r[0] = x1 * x4 * (x1 + x2 + x3) + x3;
147                        return r;
148                }
149                ADVector r(2);
150                // g_1 (x)
151                r[0] = x1 * x2 * x3 * x4;
152                // g_2 (x)
153                r[1] = x1 * x1 + x2 * x2 + x3 * x3 + x4 * x4;
154                return r;
155        }
156        bool retape(size_t k)
157        {       return retape_; }
158        size_t number_functions(void)
159        {       return 2; }
160        size_t domain_size(size_t k)
161        {       return 4; }
162        size_t range_size(size_t k)
163        {       if( k == 0 )
164                        return 1;
165                return 2;
166        }
167        size_t number_terms(size_t k)
168        {       return 1; }
169        void index(size_t k, size_t ell, SizeVector& I, SizeVector& J)
170        {
171                if( k == 0 )
172                        I[0] = 0;
173                else
174                {       I[0] = 1;
175                        I[1] = 2; 
176                }
177                // reverse the order of the variables in u from that in x
178                for(size_t j = 0; j < 4; j++)
179                        J[j] = 3-j;
180        }
181};
182       
183bool ipopt_cppad_K_gt_1(void)
184{       bool ok = true;
185        size_t j;
186
187
188        // number of independent variables (domain dimension for f and g)
189        size_t n = 4; 
190        // number of constraints (range dimension for g)
191        size_t m = 2;
192        // initial value of the independent variables
193        NumberVector x_i(n);
194        x_i[0] = 1.0;
195        x_i[1] = 5.0;
196        x_i[2] = 5.0;
197        x_i[3] = 1.0;
198        // lower and upper limits for x
199        NumberVector x_l(n);
200        NumberVector x_u(n);
201        for(j = 0; j < n; j++)
202        {       x_l[j] = 1.0;
203                x_u[j] = 5.0;
204        }
205        // lower and upper limits for g
206        NumberVector g_l(m);
207        NumberVector g_u(m);
208        g_l[0] = 25.0;     g_u[0] = 1.0e19;
209        g_l[1] = 40.0;     g_u[1] = 40.0;
210
211        // known solution to check against
212        double check_x[]   = { 1.000000, 4.743000, 3.82115, 1.379408 };
213
214        size_t icase;
215        for(icase = 0; icase <= 1; icase++)
216        {       // Should ipopt_cppad_nlp retape the operation sequence for
217                // every new x. Can test both true and false cases because
218                // the operation sequence does not depend on x (for this case).
219                bool retape = bool(icase);
220
221                // check case where upper and lower limits are equal
222                if( icase == 1 )
223                {       x_l[2] = check_x[2];
224                        x_u[2] = check_x[2];
225                }
226
227                // object in derived class
228                FG_K_gt_1 my_fg_info(retape);
229                ipopt_cppad_fg_info *fg_info = &my_fg_info; 
230
231                // create the Ipopt interface
232                ipopt_cppad_solution solution;
233                Ipopt::SmartPtr<Ipopt::TNLP> cppad_nlp = new ipopt_cppad_nlp(
234                n, m, x_i, x_l, x_u, g_l, g_u, fg_info, &solution
235                );
236
237                // Create an instance of the IpoptApplication
238                using Ipopt::IpoptApplication;
239                Ipopt::SmartPtr<IpoptApplication> app = new IpoptApplication();
240
241                // turn off any printing
242                app->Options()->SetIntegerValue("print_level", 0);
243
244                // maximum number of iterations
245                app->Options()->SetIntegerValue("max_iter", 10);
246
247                // approximate accuracy in first order necessary conditions;
248                // see Mathematical Programming, Volume 106, Number 1,
249                // Pages 25-57, Equation (6)
250                app->Options()->SetNumericValue("tol", 1e-9);
251
252                // derivative testing
253                app->Options()->
254                SetStringValue("derivative_test", "second-order");
255
256                // Initialize the IpoptApplication and process the options
257                Ipopt::ApplicationReturnStatus status = app->Initialize();
258                ok    &= status == Ipopt::Solve_Succeeded;
259
260                // Run the IpoptApplication
261                status = app->OptimizeTNLP(cppad_nlp);
262                ok    &= status == Ipopt::Solve_Succeeded;
263
264                /*
265                Check some of the solution values
266                */
267                ok &= solution.status == ipopt_cppad_solution::success;
268                //
269                double check_z_l[] = { 1.087871, 0.,       0.,      0.       };
270                double check_z_u[] = { 0.,       0.,       0.,      0.       };
271                double rel_tol     = 1e-6;  // relative tolerance
272                double abs_tol     = 1e-6;  // absolute tolerance
273                for(j = 0; j < n; j++)
274                {       ok &= CppAD::NearEqual(
275                        check_x[j],   solution.x[j],   rel_tol, abs_tol
276                        );
277                        ok &= CppAD::NearEqual(
278                        check_z_l[j], solution.z_l[j], rel_tol, abs_tol
279                        );
280                        ok &= CppAD::NearEqual(
281                        check_z_u[j], solution.z_u[j], rel_tol, abs_tol
282                        );
283                }
284        }
285
286        return ok;
287}
288// ---------------------------------------------------------------------------
289/*
290f(x)    = x[1]; k=0, ell=0, I[0] = 0, J[0] = 1
291g_0 (x) = x[0]; k=0, ell=1, I[0] = 1, J[0] = 0
292g_1 (x) = x[1]; k=0, ell=2, I[0] = 2, J[0] = 1
293
294minimize   f(x)
295subject to -1 <= g_0(x)  <= 0
296            0 <= g_1 (x) <= 1
297
298The solution is x[1] = 0 and x[0] arbitrary.
299*/
300
301namespace
302{
303class FG_J_changes : public ipopt_cppad_fg_info
304{
305private:
306        bool retape_;
307public:
308        // constructor
309        FG_J_changes(bool retape)
310        : retape_ (retape)
311        { }
312        size_t number_functions(void)
313        {       return 1; }
314        size_t domain_size(size_t k)
315        {       size_t q;
316                switch(k)
317                {       case 0:  q = 1;  break; 
318                        default: assert(0);
319                }
320                return q;
321        }
322        size_t range_size(size_t k)
323        {       size_t p;
324                switch(k)
325                {       case 0:  p = 1;  break;
326                        default: assert(0);
327                }
328                return p;
329        }
330        size_t number_terms(size_t k)
331        {       size_t L;
332                switch(k)
333                {       case 0:  L = 3;   break;
334                        default: assert(0);
335                }
336                return L;
337        }
338        void index(size_t k, size_t ell, SizeVector&I, SizeVector& J)
339        {       assert( I.size() >= 1 );
340                assert( J.size() >= 1 );
341                I[0] = ell;
342                if( ell == 0 )
343                {       J[0] = 1;
344                        return;
345                }
346                J[0] = ell - 1;
347                return;
348        }
349        // retape function
350        bool retape(size_t k)
351        {       return retape_; }
352        ADVector eval_r(size_t k, const ADVector&  u)
353        {
354                assert( u.size() == 1 );
355                ADVector r(1);
356                r[0] = u[0] ; 
357                return r;
358        }
359};
360}
361
362bool ipopt_cppad_J_changes(void)
363{
364        bool ok = true;
365        // number of independent variables (domain dimension for f and g)
366        size_t n = 2;
367        // number of constraints (range dimension for g)
368        size_t m = 2; 
369        // initial value of the independent variables
370        NumberVector x_i(n);
371        NumberVector x_l(n);
372        NumberVector x_u(n);
373
374        size_t i = 0;
375        for(i = 0; i < n; i++)
376        {       x_i[i] = 0.;
377                x_l[i] = -1.0;
378                x_u[i] = +1.0;
379        }
380
381        // lower and upper limits for g
382        NumberVector g_l(m);
383        NumberVector g_u(m);
384        g_l[0] = -1; g_u[0] = 0.;
385        g_l[1] = 0.; g_u[1] = 1.;
386
387        // object for evaluating function
388        bool retape = false;
389        FG_J_changes my_fg_info(retape);
390        ipopt_cppad_fg_info *fg_info = &my_fg_info;
391
392        ipopt_cppad_solution solution;
393        Ipopt::SmartPtr<Ipopt::TNLP> cppad_nlp = new ipopt_cppad_nlp(
394                n, m, x_i, x_l, x_u, g_l, g_u, fg_info, &solution
395        );
396
397        // Create an instance of the IpoptApplication
398        using Ipopt::IpoptApplication;
399        Ipopt::SmartPtr<IpoptApplication> app = new IpoptApplication();
400
401        // turn off any printing
402        app->Options()->SetIntegerValue("print_level", 0);
403
404        // approximate accuracy in first order necessary conditions;
405        // see Mathematical Programming, Volume 106, Number 1,
406        // Pages 25-57, Equation (6)
407        app->Options()->SetNumericValue("tol", 1e-9);
408        app->Options()-> SetStringValue("derivative_test", "second-order");
409
410        // Initialize the IpoptApplication and process the options
411        Ipopt::ApplicationReturnStatus status = app->Initialize();
412        ok    &= status == Ipopt::Solve_Succeeded;
413
414        // Run the IpoptApplication
415        status = app->OptimizeTNLP(cppad_nlp);
416        ok    &= status == Ipopt::Solve_Succeeded;
417
418        /*
419         Check solution status
420         */
421        ok &= solution.status == ipopt_cppad_solution::success;
422        ok &= CppAD::NearEqual(solution.x[1], 0., 1e-6, 1e-6);
423
424        return ok;
425}
426// ---------------------------------------------------------------------------
427
428} // End empty namespace
429
430bool ipopt_cppad(void)
431{       bool ok = true;
432        ok &= ipopt_cppad_retape();
433        ok &= ipopt_cppad_K_gt_1();
434        ok &= ipopt_cppad_J_changes();
435        return ok;
436}
Note: See TracBrowser for help on using the repository browser.