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: Andreas Waechter IBM 2005-10-18 00008 // Olaf Schenk (Univ. of Basel) 2007-08-01 00009 // modified MittelmannBndryCntrlDiri.hpp for 3-dim problem 00010 00011 #ifndef __MITTELMANNBNDRYCNTRLDIRI3D_HPP__ 00012 #define __MITTELMANNBNDRYCNTRLDIRI3D_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 : public RegisteredTNLP 00054 { 00055 public: 00057 MittelmannBndryCntrlDiriBase3D(); 00058 00060 virtual ~MittelmannBndryCntrlDiriBase3D(); 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(const MittelmannBndryCntrlDiriBase3D&); 00152 MittelmannBndryCntrlDiriBase3D& operator=(const MittelmannBndryCntrlDiriBase3D&); 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 //return 0.5*t*t; 00215 if (t > B_) { 00216 return B_*B_/2. + C_*(t - B_); 00217 } 00218 else if (t < -B_) { 00219 return B_*B_/2. + C_*(-t - B_); 00220 } 00221 else { 00222 const Number t2 = t*t; 00223 const Number t4 = t2*t2; 00224 const Number t6 = t4*t2; 00225 return PenA_*t2 + PenB_*t4 + PenC_*t6; 00226 } 00227 } 00229 inline Number PenObj_1(Number t) const 00230 { 00231 //return t; 00232 if (t > B_) { 00233 return C_; 00234 } 00235 else if (t < -B_) { 00236 return -C_; 00237 } 00238 else { 00239 const Number t2 = t*t; 00240 const Number t3 = t*t2; 00241 const Number t5 = t3*t2; 00242 return 2.*PenA_*t + 4.*PenB_*t3 + 6.*PenC_*t5; 00243 } 00244 } 00246 inline Number PenObj_2(Number t) const 00247 { 00248 //return 1.; 00249 if (t > B_) { 00250 return 0.; 00251 } 00252 else if (t < -B_) { 00253 return 0.; 00254 } 00255 else { 00256 const Number t2 = t*t; 00257 const Number t4 = t2*t2; 00258 return 2.*PenA_ + 12.*PenB_*t2 + 30.*PenC_*t4; 00259 } 00260 } 00262 00265 Number B_; 00266 Number C_; 00267 Number PenA_; 00268 Number PenB_; 00269 Number PenC_; 00271 }; 00272 00274 class MittelmannBndryCntrlDiri3D : public MittelmannBndryCntrlDiriBase3D 00275 { 00276 public: 00277 MittelmannBndryCntrlDiri3D() 00278 {} 00279 00280 virtual ~MittelmannBndryCntrlDiri3D() 00281 {} 00282 00283 virtual bool InitializeProblem(Index N) 00284 { 00285 if (N<1) { 00286 printf("N has to be at least 1."); 00287 return false; 00288 } 00289 Number alpha = 0.01; 00290 Number lb_y = -1e20; 00291 Number ub_y = 3.5; 00292 Number lb_u = 0.; 00293 Number ub_u = 10.; 00294 Number d_const = -20.; 00295 Number B = .5; 00296 Number C = 0.01; 00297 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, d_const, B, C); 00298 return true; 00299 } 00300 protected: 00302 virtual Number y_d_cont(Number x1, Number x2, Number x3) const 00303 { 00304 return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.)); 00305 } 00306 private: 00309 MittelmannBndryCntrlDiri3D(const MittelmannBndryCntrlDiri3D&); 00310 MittelmannBndryCntrlDiri3D& operator=(const MittelmannBndryCntrlDiri3D&); 00312 00313 }; 00314 00315 00316 #endif