MittelmannBndryCntrlDiri3D.hpp
Go to the documentation of this file.
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: MittelmannBndryCntrlDiri3D.hpp 2005 2011-06-06 12:55:16Z stefan $
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 
54 {
55 public:
58 
61 
65  virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
66  Index& nnz_h_lag, IndexStyleEnum& index_style);
67 
69  virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
70  Index m, Number* g_l, Number* g_u);
71 
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 
79  virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
80 
82  virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
83 
85  virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
86 
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 
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 
105 
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 
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,
124 
125 protected:
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 
136  virtual Number y_d_cont(Number x1, Number x2, Number x3) const =0;
138 
139 private:
154 
181 
186  inline Index y_index(Index i, Index j, Index k) const
187  {
188  return k + (N_+2)*j + (N_+2)*(N_+2)*i;
189  }
192  inline Index pde_index(Index i, Index j, Index k) const
193  {
194  return (k-1) + N_*(j-1) + N_*N_*(i-1);
195  }
197  inline Number x1_grid(Index i) const
198  {
199  return h_*(Number)i;
200  }
202  inline Number x2_grid(Index i) const
203  {
204  return h_*(Number)i;
205  }
207  inline Number x3_grid(Index i) const
208  {
209  return h_*(Number)i;
210  }
212  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  }
229  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  }
246  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  }
262 
271 };
272 
275 {
276 public:
278  {}
279 
281  {}
282 
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:
302  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:
312 
313 };
314 
315 
316 #endif