Ipopt
trunk
|
00001 // Copyright (C) 2005, 2006 International Business Machines and others. 00002 // All Rights Reserved. 00003 // This code is published under the Eclipse Public License. 00004 // 00005 // $Id$ 00006 // 00007 // Authors: Andreas Waechter IBM 2005-10-18 00008 // based on MyNLP.hpp 00009 00010 #ifndef __MITTELMANNDISTRCNTRLNEUMB_HPP__ 00011 #define __MITTELMANNDISTRCNTRLNEUMB_HPP__ 00012 00013 #include "IpTNLP.hpp" 00014 #include "RegisteredTNLP.hpp" 00015 00016 #ifdef HAVE_CONFIG_H 00017 #include "config.h" 00018 #else 00019 #include "configall_system.h" 00020 #endif 00021 00022 #ifdef HAVE_CMATH 00023 # include <cmath> 00024 #else 00025 # ifdef HAVE_MATH_H 00026 # include <math.h> 00027 # else 00028 # error "don't have header file for math" 00029 # endif 00030 #endif 00031 00032 #ifdef HAVE_CSTDIO 00033 # include <cstdio> 00034 #else 00035 # ifdef HAVE_STDIO_H 00036 # include <stdio.h> 00037 # else 00038 # error "don't have header file for stdio" 00039 # endif 00040 #endif 00041 00042 using namespace Ipopt; 00043 00050 class MittelmannDistCntrlNeumBBase : public RegisteredTNLP 00051 { 00052 public: 00055 MittelmannDistCntrlNeumBBase(); 00056 00058 virtual ~MittelmannDistCntrlNeumBBase(); 00059 00063 virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, 00064 Index& nnz_h_lag, IndexStyleEnum& index_style); 00065 00067 virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u, 00068 Index m, Number* g_l, Number* g_u); 00069 00071 virtual bool get_starting_point(Index n, bool init_x, Number* x, 00072 bool init_z, Number* z_L, Number* z_U, 00073 Index m, bool init_lambda, 00074 Number* lambda); 00075 00077 virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value); 00078 00080 virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f); 00081 00083 virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g); 00084 00089 virtual bool eval_jac_g(Index n, const Number* x, bool new_x, 00090 Index m, Index nele_jac, Index* iRow, Index *jCol, 00091 Number* values); 00092 00097 virtual bool eval_h(Index n, const Number* x, bool new_x, 00098 Number obj_factor, Index m, const Number* lambda, 00099 bool new_lambda, Index nele_hess, Index* iRow, 00100 Index* jCol, Number* values); 00101 00103 00105 virtual bool get_scaling_parameters(Number& obj_scaling, 00106 bool& use_x_scaling, Index n, 00107 Number* x_scaling, 00108 bool& use_g_scaling, Index m, 00109 Number* g_scaling); 00110 00115 virtual void finalize_solution(SolverReturn status, 00116 Index n, const Number* x, const Number* z_L, const Number* z_U, 00117 Index m, const Number* g, const Number* lambda, 00118 Number obj_value, 00119 const IpoptData* ip_data, 00120 IpoptCalculatedQuantities* ip_cq); 00122 00123 protected: 00127 void SetBaseParameters(Index N, Number lb_y, 00128 Number ub_y, Number lb_u, Number ub_u, 00129 Number b_0j, Number b_1j, Number b_i0, Number b_i1, 00130 Number u_init); 00131 00135 virtual Number y_d_cont(Number x1, Number x2) const =0; 00137 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const =0; 00139 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const =0; 00141 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const =0; 00143 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const =0; 00146 virtual bool fint_cont_dydy_alwayszero() const =0; 00148 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const =0; 00151 virtual bool fint_cont_dudu_alwayszero() const =0; 00153 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const =0; 00156 virtual bool fint_cont_dydu_alwayszero() const =0; 00158 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const =0; 00160 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const =0; 00162 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const =0; 00164 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const =0; 00167 virtual bool d_cont_dydy_alwayszero() const =0; 00169 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const =0; 00172 virtual bool d_cont_dudu_alwayszero() const =0; 00174 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const =0; 00177 virtual bool d_cont_dydu_alwayszero() const =0; 00179 00180 private: 00192 MittelmannDistCntrlNeumBBase(const MittelmannDistCntrlNeumBBase&); 00193 MittelmannDistCntrlNeumBBase& operator=(const MittelmannDistCntrlNeumBBase&); 00195 00199 Index N_; 00201 Number h_; 00203 Number hh_; 00205 Number lb_y_; 00207 Number ub_y_; 00209 Number lb_u_; 00211 Number ub_u_; 00214 Number b_0j_; 00217 Number b_1j_; 00220 Number b_i0_; 00223 Number b_i1_; 00225 Number u_init_; 00227 Number* y_d_; 00229 00234 inline Index y_index(Index i, Index j) const 00235 { 00236 return j + (N_+2)*i; 00237 } 00240 inline Index u_index(Index i, Index j) const 00241 { 00242 return (N_+2)*(N_+2) + (j-1) + (N_)*(i-1); 00243 } 00246 inline Index pde_index(Index i, Index j) const 00247 { 00248 return (j-1) + N_*(i-1); 00249 } 00251 inline Number x1_grid(Index i) const 00252 { 00253 return h_*(Number)i; 00254 } 00256 inline Number x2_grid(Index i) const 00257 { 00258 return h_*(Number)i; 00259 } 00261 }; 00262 00264 class MittelmannDistCntrlNeumB1 : public MittelmannDistCntrlNeumBBase 00265 { 00266 public: 00267 MittelmannDistCntrlNeumB1() 00268 : 00269 pi_(4.*atan(1.)), 00270 alpha_(0.001) 00271 {} 00272 00273 virtual ~MittelmannDistCntrlNeumB1() 00274 {} 00275 00276 virtual bool InitializeProblem(Index N) 00277 { 00278 if (N<1) { 00279 printf("N has to be at least 1."); 00280 return false; 00281 } 00282 Number lb_y = -1e20; 00283 Number ub_y = 0.371; 00284 Number lb_u = -8.; 00285 Number ub_u = 9.; 00286 Number b_0j = 1.; 00287 Number b_1j = 1.; 00288 Number b_i0 = 1.; 00289 Number b_i1 = 1.; 00290 Number u_init = (ub_u + lb_u)/2.; 00291 00292 SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init); 00293 return true; 00294 } 00295 protected: 00297 virtual Number y_d_cont(Number x1, Number x2) const 00298 { 00299 return sin(2.*pi_*x1)*sin(2.*pi_*x2); 00300 } 00302 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const 00303 { 00304 Number diff_y = y-y_d_cont(x1,x2); 00305 return 0.5*(diff_y*diff_y + alpha_*u*u); 00306 } 00308 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const 00309 { 00310 return y-y_d_cont(x1,x2); 00311 } 00312 00314 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const 00315 { 00316 return alpha_*u; 00317 } 00319 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const 00320 { 00321 return 1.; 00322 } 00325 virtual bool fint_cont_dydy_alwayszero() const 00326 { 00327 return false; 00328 } 00330 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const 00331 { 00332 return alpha_; 00333 } 00336 virtual bool fint_cont_dudu_alwayszero() const 00337 { 00338 return false; 00339 } 00341 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const 00342 { 00343 return 0.; 00344 } 00347 virtual bool fint_cont_dydu_alwayszero() const 00348 { 00349 return true; 00350 } 00352 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const 00353 { 00354 return -exp(y) - u; 00355 } 00357 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const 00358 { 00359 return -exp(y); 00360 } 00362 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const 00363 { 00364 return -1.; 00365 } 00367 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const 00368 { 00369 return -exp(y); 00370 } 00373 virtual bool d_cont_dydy_alwayszero() const 00374 { 00375 return false; 00376 } 00378 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const 00379 { 00380 return 0.; 00381 } 00384 virtual bool d_cont_dudu_alwayszero() const 00385 { 00386 return true; 00387 } 00389 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const 00390 { 00391 return 0.; 00392 } 00395 virtual bool d_cont_dydu_alwayszero() const 00396 { 00397 return true; 00398 } 00399 private: 00402 MittelmannDistCntrlNeumB1(const MittelmannDistCntrlNeumB1&); 00403 MittelmannDistCntrlNeumB1& operator=(const MittelmannDistCntrlNeumB1&); 00405 00406 const Number pi_; 00408 const Number alpha_; 00409 }; 00410 00412 class MittelmannDistCntrlNeumB2 : public MittelmannDistCntrlNeumBBase 00413 { 00414 public: 00415 MittelmannDistCntrlNeumB2() 00416 : 00417 pi_(4.*atan(1.)) 00418 {} 00419 00420 virtual ~MittelmannDistCntrlNeumB2() 00421 {} 00422 00423 virtual bool InitializeProblem(Index N) 00424 { 00425 if (N<1) { 00426 printf("N has to be at least 1."); 00427 return false; 00428 } 00429 Number lb_y = -1e20; 00430 Number ub_y = 0.371; 00431 Number lb_u = -8.; 00432 Number ub_u = 9.; 00433 Number b_0j = 1.; 00434 Number b_1j = 1.; 00435 Number b_i0 = 1.; 00436 Number b_i1 = 1.; 00437 Number u_init = (ub_u + lb_u)/2.; 00438 00439 SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init); 00440 return true; 00441 } 00442 protected: 00444 virtual Number y_d_cont(Number x1, Number x2) const 00445 { 00446 return sin(2.*pi_*x1)*sin(2.*pi_*x2); 00447 } 00449 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const 00450 { 00451 Number diff_y = y-y_d_cont(x1,x2); 00452 return 0.5*diff_y*diff_y; 00453 } 00455 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const 00456 { 00457 return y-y_d_cont(x1,x2); 00458 } 00459 00461 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const 00462 { 00463 return 0.; 00464 } 00466 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const 00467 { 00468 return 1.; 00469 } 00472 virtual bool fint_cont_dydy_alwayszero() const 00473 { 00474 return false; 00475 } 00477 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const 00478 { 00479 return 0.; 00480 } 00483 virtual bool fint_cont_dudu_alwayszero() const 00484 { 00485 return true; 00486 } 00488 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const 00489 { 00490 return 0.; 00491 } 00494 virtual bool fint_cont_dydu_alwayszero() const 00495 { 00496 return true; 00497 } 00499 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const 00500 { 00501 return -exp(y) - u; 00502 } 00504 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const 00505 { 00506 return -exp(y); 00507 } 00509 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const 00510 { 00511 return -1.; 00512 } 00514 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const 00515 { 00516 return -exp(y); 00517 } 00520 virtual bool d_cont_dydy_alwayszero() const 00521 { 00522 return false; 00523 } 00525 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const 00526 { 00527 return 0.; 00528 } 00531 virtual bool d_cont_dudu_alwayszero() const 00532 { 00533 return true; 00534 } 00536 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const 00537 { 00538 return 0.; 00539 } 00542 virtual bool d_cont_dydu_alwayszero() const 00543 { 00544 return true; 00545 } 00546 private: 00549 MittelmannDistCntrlNeumB2(const MittelmannDistCntrlNeumB2&); 00550 MittelmannDistCntrlNeumB2& operator=(const MittelmannDistCntrlNeumB2&); 00552 00553 const Number pi_; 00554 }; 00555 00557 class MittelmannDistCntrlNeumB3 : public MittelmannDistCntrlNeumBBase 00558 { 00559 public: 00560 MittelmannDistCntrlNeumB3() 00561 : 00562 pi_(4.*atan(1.)), 00563 M_(1.), 00564 K_(0.8), 00565 b_(1.) 00566 {} 00567 00568 virtual ~MittelmannDistCntrlNeumB3() 00569 {} 00570 00571 virtual bool InitializeProblem(Index N) 00572 { 00573 if (N<1) { 00574 printf("N has to be at least 1."); 00575 return false; 00576 } 00577 Number lb_y = 3.;//-1e20; 00578 Number ub_y = 6.09; 00579 Number lb_u = 1.4; 00580 Number ub_u = 1.6; 00581 Number b_0j = 1.; 00582 Number b_1j = 0.; 00583 Number b_i0 = 1.; 00584 Number b_i1 = 0.; 00585 Number u_init = (ub_u + lb_u)/2.; 00586 00587 SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init); 00588 return true; 00589 } 00590 protected: 00592 virtual Number y_d_cont(Number x1, Number x2) const 00593 { 00594 return 6.; 00595 } 00597 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const 00598 { 00599 return u*(M_*u - K_*y); 00600 } 00602 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const 00603 { 00604 return -K_*u; 00605 } 00606 00608 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const 00609 { 00610 return 2.*M_*u - K_*y; 00611 } 00613 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const 00614 { 00615 return 0.; 00616 } 00619 virtual bool fint_cont_dydy_alwayszero() const 00620 { 00621 return true; 00622 } 00624 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const 00625 { 00626 return 2.*M_; 00627 } 00630 virtual bool fint_cont_dudu_alwayszero() const 00631 { 00632 return false; 00633 } 00635 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const 00636 { 00637 return -K_; 00638 } 00641 virtual bool fint_cont_dydu_alwayszero() const 00642 { 00643 return false; 00644 } 00646 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const 00647 { 00648 return y*(u + b_*y - a(x1,x2)); 00649 } 00651 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const 00652 { 00653 return (u + 2.*b_*y -a(x1,x2)); 00654 } 00656 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const 00657 { 00658 return y; 00659 } 00661 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const 00662 { 00663 return 2.*b_; 00664 } 00667 virtual bool d_cont_dydy_alwayszero() const 00668 { 00669 return false; 00670 } 00672 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const 00673 { 00674 return 0.; 00675 } 00678 virtual bool d_cont_dudu_alwayszero() const 00679 { 00680 return true; 00681 } 00683 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const 00684 { 00685 return 1.; 00686 } 00689 virtual bool d_cont_dydu_alwayszero() const 00690 { 00691 return false; 00692 } 00693 private: 00696 MittelmannDistCntrlNeumB3(const MittelmannDistCntrlNeumB3&); 00697 MittelmannDistCntrlNeumB3& operator=(const MittelmannDistCntrlNeumB3&); 00699 00700 const Number pi_; 00701 /*@name constrants appearing in problem formulation */ 00703 const Number M_; 00704 const Number K_; 00705 const Number b_; 00707 //* Auxiliary function for state equation */ 00708 inline Number a(Number x1, Number x2) const 00709 { 00710 return 7. + 4.*sin(2.*pi_*x1*x2); 00711 } 00712 }; 00713 00714 #endif