Ipopt
trunk
|
00001 // Copyright (C) 2005, 2007 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: Olaf Schenk (Univ. of Basel) 2007-08-01 00008 // modified MittelmannBndryCntrlDiri.hpp for 3-dim problem 00009 // based on MyNLP.hpp 00010 00011 #ifndef __MITTELMANNBNDRYCNTRLDIRI3D_27_HPP__ 00012 #define __MITTELMANNBNDRYCNTRLDIRI3D_27_HPP__ 00013 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 00053 class MittelmannBndryCntrlDiriBase3D_27 : public RegisteredTNLP 00054 { 00055 public: 00057 MittelmannBndryCntrlDiriBase3D_27(); 00058 00060 virtual ~MittelmannBndryCntrlDiriBase3D_27(); 00061 00065 virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, 00066 Index& nnz_h_lag, IndexStyleEnum& index_style); 00067 00069 virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u, 00070 Index m, Number* g_l, Number* g_u); 00071 00073 virtual bool get_starting_point(Index n, bool init_x, Number* x, 00074 bool init_z, Number* z_L, Number* z_U, 00075 Index m, bool init_lambda, 00076 Number* lambda); 00077 00079 virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value); 00080 00082 virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f); 00083 00085 virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g); 00086 00091 virtual bool eval_jac_g(Index n, const Number* x, bool new_x, 00092 Index m, Index nele_jac, Index* iRow, Index *jCol, 00093 Number* values); 00094 00099 virtual bool eval_h(Index n, const Number* x, bool new_x, 00100 Number obj_factor, Index m, const Number* lambda, 00101 bool new_lambda, Index nele_hess, Index* iRow, 00102 Index* jCol, Number* values); 00103 00105 00107 virtual bool get_scaling_parameters(Number& obj_scaling, 00108 bool& use_x_scaling, Index n, 00109 Number* x_scaling, 00110 bool& use_g_scaling, Index m, 00111 Number* g_scaling); 00112 00117 virtual void finalize_solution(SolverReturn status, 00118 Index n, const Number* x, const Number* z_L, const Number* z_U, 00119 Index m, const Number* g, const Number* lambda, 00120 Number obj_valu, 00121 const IpoptData* ip_data, 00122 IpoptCalculatedQuantities* ip_cq); 00124 00125 protected: 00129 void SetBaseParameters(Index N, Number alpha, Number lb_y, 00130 Number ub_y, Number lb_u, Number ub_u, 00131 Number d_const, Number B, Number C); 00132 00136 virtual Number y_d_cont(Number x1, Number x2, Number x3) const =0; 00138 00139 private: 00151 MittelmannBndryCntrlDiriBase3D_27(const MittelmannBndryCntrlDiriBase3D_27&); 00152 MittelmannBndryCntrlDiriBase3D_27& operator=(const MittelmannBndryCntrlDiriBase3D_27&); 00154 00158 Index N_; 00160 Number h_; 00162 Number hh_; 00164 Number hhh_; 00166 Number lb_y_; 00168 Number ub_y_; 00170 Number lb_u_; 00172 Number ub_u_; 00174 Number d_const_; 00177 Number alpha_; 00179 Number* y_d_; 00181 00186 inline Index y_index(Index i, Index j, Index k) const 00187 { 00188 return k + (N_+2)*j + (N_+2)*(N_+2)*i; 00189 } 00192 inline Index pde_index(Index i, Index j, Index k) const 00193 { 00194 return (k-1) + N_*(j-1) + N_*N_*(i-1); 00195 } 00197 inline Number x1_grid(Index i) const 00198 { 00199 return h_*(Number)i; 00200 } 00202 inline Number x2_grid(Index i) const 00203 { 00204 return h_*(Number)i; 00205 } 00207 inline Number x3_grid(Index i) const 00208 { 00209 return h_*(Number)i; 00210 } 00212 inline Number PenObj(Number t) const 00213 { 00214 if (B_ == 0.) { 00215 return 0.5*t*t; 00216 } 00217 else if (t > B_) { 00218 return B_*B_/2. + C_*(t - B_); 00219 } 00220 else if (t < -B_) { 00221 return B_*B_/2. + C_*(-t - B_); 00222 } 00223 else { 00224 const Number t2 = t*t; 00225 const Number t4 = t2*t2; 00226 const Number t6 = t4*t2; 00227 return PenA_*t2 + PenB_*t4 + PenC_*t6; 00228 } 00229 } 00231 inline Number PenObj_1(Number t) const 00232 { 00233 if (B_ == 0.) { 00234 return t; 00235 } 00236 else if (t > B_) { 00237 return C_; 00238 } 00239 else if (t < -B_) { 00240 return -C_; 00241 } 00242 else { 00243 const Number t2 = t*t; 00244 const Number t3 = t*t2; 00245 const Number t5 = t3*t2; 00246 return 2.*PenA_*t + 4.*PenB_*t3 + 6.*PenC_*t5; 00247 } 00248 } 00250 inline Number PenObj_2(Number t) const 00251 { 00252 if (B_ == 0.) { 00253 return 1.; 00254 } 00255 else if (t > B_) { 00256 return 0.; 00257 } 00258 else if (t < -B_) { 00259 return 0.; 00260 } 00261 else { 00262 const Number t2 = t*t; 00263 const Number t4 = t2*t2; 00264 return 2.*PenA_ + 12.*PenB_*t2 + 30.*PenC_*t4; 00265 } 00266 } 00268 00271 Number B_; 00272 Number C_; 00273 Number PenA_; 00274 Number PenB_; 00275 Number PenC_; 00277 }; 00278 00280 class MittelmannBndryCntrlDiri3D_27 : public MittelmannBndryCntrlDiriBase3D_27 00281 { 00282 public: 00283 MittelmannBndryCntrlDiri3D_27() 00284 {} 00285 00286 virtual ~MittelmannBndryCntrlDiri3D_27() 00287 {} 00288 00289 virtual bool InitializeProblem(Index N) 00290 { 00291 if (N<1) { 00292 printf("N has to be at least 1."); 00293 return false; 00294 } 00295 Number alpha = 1e-2; 00296 Number lb_y = -1e20; 00297 Number ub_y = 3.5; 00298 Number lb_u = 0.; 00299 Number ub_u = 10.; 00300 Number d_const = -20.; 00301 Number B = 0.; // convex case (quadratic penalty) 00302 Number C = 0.; 00303 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, d_const, B, C); 00304 return true; 00305 } 00306 protected: 00308 virtual Number y_d_cont(Number x1, Number x2, Number x3) const 00309 { 00310 return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.)); 00311 } 00312 private: 00315 MittelmannBndryCntrlDiri3D_27(const MittelmannBndryCntrlDiri3D_27&); 00316 MittelmannBndryCntrlDiri3D_27& operator=(const MittelmannBndryCntrlDiri3D_27&); 00318 00319 }; 00320 00323 class MittelmannBndryCntrlDiri3D_27BT : public MittelmannBndryCntrlDiriBase3D_27 00324 { 00325 public: 00326 MittelmannBndryCntrlDiri3D_27BT() 00327 {} 00328 00329 virtual ~MittelmannBndryCntrlDiri3D_27BT() 00330 {} 00331 00332 virtual bool InitializeProblem(Index N) 00333 { 00334 if (N<1) { 00335 printf("N has to be at least 1."); 00336 return false; 00337 } 00338 Number alpha = 1e-2; 00339 Number lb_y = -1e20; 00340 Number ub_y = 3.5; 00341 Number lb_u = 0.; 00342 Number ub_u = 10.; 00343 Number d_const = -20.; 00344 Number B = .25; // nonconves case with beaton-tukey-type penalty function 00345 Number C = 0.01; 00346 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, d_const, B, C); 00347 return true; 00348 } 00349 protected: 00351 virtual Number y_d_cont(Number x1, Number x2, Number x3) const 00352 { 00353 return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.)); 00354 } 00355 private: 00358 MittelmannBndryCntrlDiri3D_27BT(const MittelmannBndryCntrlDiri3D_27BT&); 00359 MittelmannBndryCntrlDiri3D_27BT& operator=(const MittelmannBndryCntrlDiri3D_27BT&); 00361 00362 }; 00363 00364 #endif