Ipopt  trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
MittelmannBndryCntrlDiri3D_27.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines