Changeset 1277 for trunk/test_more
 Timestamp:
 Sep 9, 2008 8:57:06 AM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/test_more/ipopt_cppad.cpp
r1276 r1277 12 12 # include "../example/ipopt_cppad_nlp.hpp" 13 13 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) 14 namespace { // Begin empty namespace 15 16 class FG_retape : public ipopt_cppad_fg_info 17 { 18 public: 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 40 bool ipopt_cppad_retape(void) 42 41 { bool ok = true; 43 42 size_t j; … … 64 63 65 64 // 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; 68 67 69 68 // create the Ipopt interface … … 116 115 return ok; 117 116 } 117 118 /* 119 This solve the same problem as ../example/ipopt_cppad.cpp (repository revision 1276) in a convoluted way in order to test the representation code. 120 */ 121 class FG_K_gt_1 : public ipopt_cppad_fg_info 122 { 123 private: 124 bool retape_; 125 public: 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] = 3j; 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 182 bool 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 2557, Equation (6) 240 app>Options()>SetNumericValue("tol", 1e9); 241 242 // derivative testing 243 app>Options()> 244 SetStringValue("derivative_test", "secondorder"); 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 = 1e6; // relative tolerance 263 double abs_tol = 1e6; // 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 /* 282 f(x) = x[1]; k=0, ell=0, I[0] = 0, J[0] = 1 283 g_0 (x) = x[0]; k=0, ell=1, I[0] = 1, J[0] = 0 284 g_1 (x) = x[1]; k=0, ell=2, I[0] = 2, J[0] = 1 285 286 minimize f(x) 287 subject to 1 <= g_0(x) <= 0 288 0 <= g_1 (x) <= 1 289 290 The solution is x[1] = 0 and x[0] arbitrary. 291 */ 292 293 namespace 294 { 295 class FG_J_changes : public ipopt_cppad_fg_info 296 { 297 private: 298 bool retape_; 299 public: 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 357 bool 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 2557, Equation (6) 402 app>Options()>SetNumericValue("tol", 1e9); 403 app>Options()> SetStringValue("derivative_test", "secondorder"); 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., 1e6, 1e6); 418 419 return ok; 420 } 421 422 423 424 } // End empty namespace 425 426 bool 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.