# source:trunk/Ipopt/examples/ScalableProblems/MittelmannBndryCntrlDiri3D_27.hpp@1861

Last change on this file since 1861 was 1861, checked in by andreasw, 3 years ago

moved Ipopt trunk from Common Public License to successor Eclispe Public License (see e.g. http://www.ibm.com/developerworks/library/os-cplfaq.html)

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