1 // Copyright (C) 2005, 2007 International Business Machines and others. 2 // All Rights Reserved. 3 // This code is published under the Eclipse Public License. 4 // 5 // $Id$ 6 // 7 // Authors: Andreas Waechter IBM 2005-10-18 8 // Olaf Schenk (Univ. of Basel) 2007-08-01 9 // modified MittelmannBndryCntrlDiri.hpp for 3-dim problem 10 11 #ifndef __MITTELMANNBNDRYCNTRLDIRI3D_HPP__ 12 #define __MITTELMANNBNDRYCNTRLDIRI3D_HPP__ 13 14 #include "RegisteredTNLP.hpp" 15 16 #ifdef HAVE_CONFIG_H 17 #include "config.h" 18 #else 19 #include "configall_system.h" 20 #endif 21 22 #ifdef HAVE_CMATH 23 # include <cmath> 24 #else 25 # ifdef HAVE_MATH_H 26 # include <math.h> 27 # else 28 # error "don't have header file for math" 29 # endif 30 #endif 31 32 #ifdef HAVE_CSTDIO 33 # include <cstdio> 34 #else 35 # ifdef HAVE_STDIO_H 36 # include <stdio.h> 37 # else 38 # error "don't have header file for stdio" 39 # endif 40 #endif 41 42 using namespace Ipopt; 43 44 /** Base class for boundary control problems with Dirichlet boundary 45 * conditions, as formulated by Hans Mittelmann as Examples 1-4 in 46 * "Optimization Techniques for Solving Elliptic Control Problems 47 * with Control and State Constraints. Part 2: Boundary Control" 48 * 49 * Here, the control variables are identical to the values of y on 50 * the boundary, and therefore we don't need any explicit 51 * optimization variables for u. 52 */ 53 class MittelmannBndryCntrlDiriBase3D : public RegisteredTNLP 54 { 55 public: 56 /** Constructor. */ 57 MittelmannBndryCntrlDiriBase3D(); 58 59 /** Default destructor */ 60 virtual ~MittelmannBndryCntrlDiriBase3D(); 61 62 /**@name Overloaded from TNLP */ 63 //@{ 64 /** Method to return some info about the nlp */ 65 virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, 66 Index& nnz_h_lag, IndexStyleEnum& index_style); 67 68 /** Method to return the bounds for my problem */ 69 virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u, 70 Index m, Number* g_l, Number* g_u); 71 72 /** Method to return the starting point for the algorithm */ 73 virtual bool get_starting_point(Index n, bool init_x, Number* x, 74 bool init_z, Number* z_L, Number* z_U, 75 Index m, bool init_lambda, 76 Number* lambda); 77 78 /** Method to return the objective value */ 79 virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value); 80 81 /** Method to return the gradient of the objective */ 82 virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f); 83 84 /** Method to return the constraint residuals */ 85 virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g); 86 87 /** Method to return: 88 * 1) The structure of the jacobian (if "values" is NULL) 89 * 2) The values of the jacobian (if "values" is not NULL) 90 */ 91 virtual bool eval_jac_g(Index n, const Number* x, bool new_x, 92 Index m, Index nele_jac, Index* iRow, Index *jCol, 93 Number* values); 94 95 /** Method to return: 96 * 1) The structure of the hessian of the lagrangian (if "values" is NULL) 97 * 2) The values of the hessian of the lagrangian (if "values" is not NULL) 98 */ 99 virtual bool eval_h(Index n, const Number* x, bool new_x, 100 Number obj_factor, Index m, const Number* lambda, 101 bool new_lambda, Index nele_hess, Index* iRow, 102 Index* jCol, Number* values); 103 104 //@} 105 106 /** Method for returning scaling parameters */ 107 virtual bool get_scaling_parameters(Number& obj_scaling, 108 bool& use_x_scaling, Index n, 109 Number* x_scaling, 110 bool& use_g_scaling, Index m, 111 Number* g_scaling); 112 113 /** @name Solution Methods */ 114 //@{ 115 /** This method is called after the optimization, and could write an 116 * output file with the optimal profiles */ 117 virtual void finalize_solution(SolverReturn status, 118 Index n, const Number* x, const Number* z_L, const Number* z_U, 119 Index m, const Number* g, const Number* lambda, 120 Number obj_valu, 121 const IpoptData* ip_data, 122 IpoptCalculatedQuantities* ip_cq); 123 //@} 124 125 protected: 126 /** Method for setting the internal parameters that define the 127 * problem. It must be called by the child class in its 128 * implementation of InitializeParameters. */ 129 void SetBaseParameters(Index N, Number alpha, Number lb_y, 130 Number ub_y, Number lb_u, Number ub_u, 131 Number d_const, Number B, Number C); 132 133 /**@name Functions that defines a particular instance. */ 134 //@{ 135 /** Target profile function for y */ 136 virtual Number y_d_cont(Number x1, Number x2, Number x3) const =0; 137 //@} 138 139 private: 140 /**@name Methods to block default compiler methods. 141 * The compiler automatically generates the following three methods. 142 * Since the default compiler implementation is generally not what 143 * you want (for all but the most simple classes), we usually 144 * put the declarations of these methods in the private section 145 * and never implement them. This prevents the compiler from 146 * implementing an incorrect "default" behavior without us 147 * knowing. (See Scott Meyers book, "Effective C++") 148 * 149 */ 150 //@{ 151 MittelmannBndryCntrlDiriBase3D(const MittelmannBndryCntrlDiriBase3D&); 152 MittelmannBndryCntrlDiriBase3D& operator=(const MittelmannBndryCntrlDiriBase3D&); 153 //@} 154 155 /**@name Problem specification */ 156 //@{ 157 /** Number of mesh points in one dimension (excluding boundary) */ 158 Index N_; 159 /** Step size */ 160 Number h_; 161 /** h_ squared */ 162 Number hh_; 163 /** h_ to the third power */ 164 Number hhh_; 165 /** overall lower bound on y */ 166 Number lb_y_; 167 /** overall upper bound on y */ 168 Number ub_y_; 169 /** overall lower bound on u */ 170 Number lb_u_; 171 /** overall upper bound on u */ 172 Number ub_u_; 173 /** Constant value of d appearing in elliptical equation */ 174 Number d_const_; 175 /** Weighting parameter for the control target deviation functional 176 * in the objective */ 177 Number alpha_; 178 /** Array for the target profile for y */ 179 Number* y_d_; 180 //@} 181 182 /**@name Auxilliary methods */ 183 //@{ 184 /** Translation of mesh point indices to NLP variable indices for 185 * y(x_ijk) */ y_index(Index i,Index j,Index k) const186 inline Index y_index(Index i, Index j, Index k) const 187 { 188 return k + (N_+2)*j + (N_+2)*(N_+2)*i; 189 } 190 /** Translation of interior mesh point indices to the corresponding 191 * PDE constraint number */ pde_index(Index i,Index j,Index k) const192 inline Index pde_index(Index i, Index j, Index k) const 193 { 194 return (k-1) + N_*(j-1) + N_*N_*(i-1); 195 } 196 /** Compute the grid coordinate for given index in x1 direction */ x1_grid(Index i) const197 inline Number x1_grid(Index i) const 198 { 199 return h_*(Number)i; 200 } 201 /** Compute the grid coordinate for given index in x2 direction */ x2_grid(Index i) const202 inline Number x2_grid(Index i) const 203 { 204 return h_*(Number)i; 205 } 206 /** Compute the grid coordinate for given index in x3 direction */ x3_grid(Index i) const207 inline Number x3_grid(Index i) const 208 { 209 return h_*(Number)i; 210 } 211 /** value of penalty function term */ PenObj(Number t) const212 inline Number PenObj(Number t) const 213 { 214 //return 0.5*t*t; 215 if (t > B_) { 216 return B_*B_/2. + C_*(t - B_); 217 } 218 else if (t < -B_) { 219 return B_*B_/2. + C_*(-t - B_); 220 } 221 else { 222 const Number t2 = t*t; 223 const Number t4 = t2*t2; 224 const Number t6 = t4*t2; 225 return PenA_*t2 + PenB_*t4 + PenC_*t6; 226 } 227 } 228 /** first derivative of penalty function term */ PenObj_1(Number t) const229 inline Number PenObj_1(Number t) const 230 { 231 //return t; 232 if (t > B_) { 233 return C_; 234 } 235 else if (t < -B_) { 236 return -C_; 237 } 238 else { 239 const Number t2 = t*t; 240 const Number t3 = t*t2; 241 const Number t5 = t3*t2; 242 return 2.*PenA_*t + 4.*PenB_*t3 + 6.*PenC_*t5; 243 } 244 } 245 /** second derivative of penalty function term */ PenObj_2(Number t) const246 inline Number PenObj_2(Number t) const 247 { 248 //return 1.; 249 if (t > B_) { 250 return 0.; 251 } 252 else if (t < -B_) { 253 return 0.; 254 } 255 else { 256 const Number t2 = t*t; 257 const Number t4 = t2*t2; 258 return 2.*PenA_ + 12.*PenB_*t2 + 30.*PenC_*t4; 259 } 260 } 261 //@} 262 263 /** @name Data for penalty function term */ 264 //@{ 265 Number B_; 266 Number C_; 267 Number PenA_; 268 Number PenB_; 269 Number PenC_; 270 //@} 271 }; 272 273 /** Class implementating Example 1 */ 274 class MittelmannBndryCntrlDiri3D : public MittelmannBndryCntrlDiriBase3D 275 { 276 public: MittelmannBndryCntrlDiri3D()277 MittelmannBndryCntrlDiri3D() 278 {} 279 ~MittelmannBndryCntrlDiri3D()280 virtual ~MittelmannBndryCntrlDiri3D() 281 {} 282 InitializeProblem(Index N)283 virtual bool InitializeProblem(Index N) 284 { 285 if (N<1) { 286 printf("N has to be at least 1."); 287 return false; 288 } 289 Number alpha = 0.01; 290 Number lb_y = -1e20; 291 Number ub_y = 3.5; 292 Number lb_u = 0.; 293 Number ub_u = 10.; 294 Number d_const = -20.; 295 Number B = .5; 296 Number C = 0.01; 297 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, d_const, B, C); 298 return true; 299 } 300 protected: 301 /** Target profile function for y */ y_d_cont(Number x1,Number x2,Number x3) const302 virtual Number y_d_cont(Number x1, Number x2, Number x3) const 303 { 304 return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.)); 305 } 306 private: 307 /**@name hide implicitly defined contructors copy operators */ 308 //@{ 309 MittelmannBndryCntrlDiri3D(const MittelmannBndryCntrlDiri3D&); 310 MittelmannBndryCntrlDiri3D& operator=(const MittelmannBndryCntrlDiri3D&); 311 //@} 312 313 }; 314 315 316 #endif 317