- Main Page
- Classes
- Files

- File List
- File Members

# examples/ScalableProblems/MittelmannParaCntrl.hpp

Go to the documentation of this file.

    00001 // Copyright (C) 2005, 2006 International Business Machines and others.
    00002 // All Rights Reserved.
    00003 // This code is published under the Eclipse Public License.
    00004 //
    00005 // $Id: MittelmannParaCntrl.hpp 2005 2011-06-06 12:55:16Z stefan $
    00006 //
    00007 // Authors:  Andreas Waechter           IBM    2005-10-18
    00008
    00009 #ifndef __MITTELMANNPARACNTRL_HPP__
    00010 #define __MITTELMANNPARACNTRL_HPP__
    00011
    00012 #include "RegisteredTNLP.hpp"
    00013
    00014 #ifdef HAVE_CONFIG_H
    00015 #include "config.h"
    00016 #else
    00017 #include "configall_system.h"
    00018 #endif
    00019
    00020 #ifdef HAVE_CMATH
    00021 # include <cmath>
    00022 #else
    00023 # ifdef HAVE_MATH_H
    00024 #  include <math.h>
    00025 # else
    00026 #  error "don't have header file for math"
    00027 # endif
    00028 #endif
    00029
    00030 using namespace Ipopt;
    00031
    00037 template<class T>
    00038 class MittelmannParaCntrlBase : public RegisteredTNLP
    00039 {
    00040 public:
    00042   MittelmannParaCntrlBase();
    00043
    00045   virtual ~MittelmannParaCntrlBase();
    00046
    00050   virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
    00051                             Index& nnz_h_lag, IndexStyleEnum& index_style);
    00052
    00054   virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
    00055                                Index m, Number* g_l, Number* g_u);
    00056
    00058   virtual bool get_starting_point(Index n, bool init_x, Number* x,
    00059                                   bool init_z, Number* z_L, Number* z_U,
    00060                                   Index m, bool init_lambda,
    00061                                   Number* lambda);
    00062
    00064   virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
    00065
    00067   virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
    00068
    00070   virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
    00071
    00076   virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
    00077                           Index m, Index nele_jac, Index* iRow, Index *jCol,
    00078                           Number* values);
    00079
    00084   virtual bool eval_h(Index n, const Number* x, bool new_x,
    00085                       Number obj_factor, Index m, const Number* lambda,
    00086                       bool new_lambda, Index nele_hess, Index* iRow,
    00087                       Index* jCol, Number* values);
    00088
    00090
    00092   virtual bool get_scaling_parameters(Number& obj_scaling,
    00093                                       bool& use_x_scaling, Index n,
    00094                                       Number* x_scaling,
    00095                                       bool& use_g_scaling, Index m,
    00096                                       Number* g_scaling);
    00097
    00102   virtual void finalize_solution(SolverReturn status,
    00103                                  Index n, const Number* x, const Number* z_L, const Number* z_U,
    00104                                  Index m, const Number* g, const Number* lambda,
    00105                                  Number obj_value,
    00106                  const IpoptData* ip_data,
    00107                  IpoptCalculatedQuantities* ip_cq);
    00109
    00110   virtual bool InitializeProblem(Index N);
    00111
    00112 private:
    00124   MittelmannParaCntrlBase(const MittelmannParaCntrlBase<T>&);
    00125   MittelmannParaCntrlBase& operator=(const MittelmannParaCntrlBase<T>&);
    00127
    00131   Number T_;
    00133   Number l_;
    00135   Index Nt_;
    00137   Index Nx_;
    00139   Number dt_;
    00141   Number dx_;
    00143   Number lb_y_;
    00145   Number ub_y_;
    00147   Number lb_u_;
    00149   Number ub_u_;
    00152   Number alpha_;
    00154   Number beta_;
    00156   Number* y_T_;
    00158   Number* a_y_;
    00160   Number* a_u_;
    00162
    00167   inline Index y_index(Index jx, Index it) const
    00168   {
    00169     return jx + (Nx_+1)*it;
    00170   }
    00171   inline Index u_index(Index it) const
    00172   {
    00173     return (Nt_+1)*(Nx_+1) + it - 1;
    00174   }
    00176   inline Number t_grid(Index i) const
    00177   {
    00178     return dt_*(Number)i;
    00179   }
    00181   inline Number x_grid(Index j) const
    00182   {
    00183     return dx_*(Number)j;
    00184   }
    00186 };
    00187
    00188 template <class T>
    00189 MittelmannParaCntrlBase<T>::MittelmannParaCntrlBase()
    00190     :
    00191     y_T_(NULL),
    00192     a_y_(NULL),
    00193     a_u_(NULL)
    00194 {}
    00195
    00196 template <class T>
    00197 MittelmannParaCntrlBase<T>::~MittelmannParaCntrlBase()
    00198 {
    00199   delete [] y_T_;
    00200   delete [] a_y_;
    00201   delete [] a_u_;
    00202 }
    00203
    00204 template <class T>
    00205 bool MittelmannParaCntrlBase<T>::
    00206 InitializeProblem(Index N)
    00207 {
    00208   typename T::ProblemSpecs p;
    00209
    00210   if (N<1) {
    00211     printf("N has to be at least 1.");
    00212     return false;
    00213   }
    00214
    00215   T_ = p.T();
    00216   l_ = p.l();
    00217   Nt_ = N;
    00218   Nx_ = N;
    00219   dt_ = T_/Nt_;
    00220   dx_ = l_/Nx_;
    00221   lb_y_ = p.lb_y();
    00222   ub_y_ = p.ub_y();
    00223   lb_u_ = p.lb_u();
    00224   ub_u_ = p.ub_u();
    00225   alpha_ = p.alpha();
    00226   beta_ = p.beta();
    00227
    00228   y_T_ = new Number[Nx_+1];
    00229   for (Index j=0; j<=Nx_; j++) {
    00230     y_T_[j] = p.y_T(x_grid(j));
    00231   }
    00232   a_y_ = new Number[Nt_];
    00233   for (Index i=1; i<=Nt_; i++) {
    00234     a_y_[i-1] = p.a_y(t_grid(i));
    00235   }
    00236   a_u_ = new Number[Nt_];
    00237   for (Index i=1; i<=Nt_; i++) {
    00238     a_u_[i-1] = p.a_u(t_grid(i));
    00239   }
    00240
    00241   return true;
    00242 }
    00243
    00244 template <class T>
    00245 bool MittelmannParaCntrlBase<T>::
    00246 get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
    00247              Index& nnz_h_lag, IndexStyleEnum& index_style)
    00248 {
    00249   typename T::ProblemSpecs p;
    00250
    00251   n = (Nt_+1)*(Nx_+1) + Nt_;
    00252
    00253   m = Nt_*(Nx_-1) + Nt_ + Nt_;
    00254
    00255   nnz_jac_g = 6*Nt_*(Nx_-1) + 3*Nt_ + 4*Nt_;
    00256
    00257   nnz_h_lag = Nx_+1 + Nt_;
    00258   if (!p.phi_dydy_always_zero()) {
    00259     nnz_h_lag += Nt_;
    00260   }
    00261
    00262   index_style = C_STYLE;
    00263
    00264   return true;
    00265 }
    00266
    00267 template <class T>
    00268 bool MittelmannParaCntrlBase<T>::
    00269 get_bounds_info(Index n, Number* x_l, Number* x_u,
    00270                 Index m, Number* g_l, Number* g_u)
    00271 {
    00272   typename T::ProblemSpecs p;
    00273
    00274   // Set overall bounds on the variables
    00275   for (Index jx=0; jx<=Nx_; jx++) {
    00276     for (Index it=1; it<=Nt_; it++) {
    00277       Index iy = y_index(jx,it);
    00278       x_l[iy] = lb_y_;
    00279       x_u[iy] = ub_y_;
    00280     }
    00281   }
    00282   for (Index i=1; i<=Nt_; i++) {
    00283     Index iu = u_index(i);
    00284     x_l[iu] = lb_u_;
    00285     x_u[iu] = ub_u_;
    00286   }
    00287
    00288   /*
    00289   // Boundary condition for t=0
    00290   for (Index it=0; it<=Nt_; it++) {
    00291     Index iy = y_index(0,it);
    00292     x_u[iy] = x_l[iy] = p.a(t_grid(it));
    00293   }
    00294   */
    00295   // Boundary condition for t=0
    00296   for (Index jx=0; jx<=Nx_; jx++) {
    00297     Index iy = y_index(jx,0);
    00298     x_u[iy] = x_l[iy] = p.a(x_grid(jx));
    00299   }
    00300
    00301   // all discretized PDE constraints have right hand side zero
    00302   for (Index i=0; i<Nt_*(Nx_-1) + Nt_; i++) {
    00303     g_l[i] = 0.;
    00304     g_u[i] = 0.;
    00305   }
    00306
    00307   // but we put b on the right hand side for the x=L boundary condition
    00308   for (Index i=0; i<Nt_; i++) {
    00309     g_l[Nt_*(Nx_-1) + Nt_ + i]
    00310     = g_u[Nt_*(Nx_-1) + Nt_ + i]
    00311       = p.b(t_grid(i+1));
    00312   }
    00313
    00314   return true;
    00315 }
    00316
    00317 template <class T>
    00318 bool MittelmannParaCntrlBase<T>::
    00319 get_starting_point(Index n, bool init_x, Number* x,
    00320                    bool init_z, Number* z_L, Number* z_U,
    00321                    Index m, bool init_lambda,
    00322                    Number* lambda)
    00323 {
    00324   DBG_ASSERT(init_x==true && init_z==false && init_lambda==false);
    00325
    00326   // Set starting point for y
    00327   for (Index jx=0; jx<=Nx_; jx++) {
    00328     for (Index it=0; it<=Nt_; it++) {
    00329       x[y_index(jx,it)] = 0.;
    00330     }
    00331   }
    00332   // for u
    00333   for (Index i=1; i<=Nt_; i++) {
    00334     x[u_index(i)] = (ub_u_+lb_u_)/2.;
    00335   }
    00336
    00337   /*
    00338   // DELETEME
    00339   for (Index i=0; i<n; i++) {
    00340     x[i] += 0.01*i;
    00341   }
    00342   */
    00343
    00344   return true;
    00345 }
    00346
    00347 template <class T>
    00348 bool MittelmannParaCntrlBase<T>::
    00349 get_scaling_parameters(Number& obj_scaling,
    00350                        bool& use_x_scaling,
    00351                        Index n, Number* x_scaling,
    00352                        bool& use_g_scaling,
    00353                        Index m, Number* g_scaling)
    00354 {
    00355   obj_scaling = 1./Min(dx_,dt_);
    00356   use_x_scaling = false;
    00357   use_g_scaling = false;
    00358   return true;
    00359 }
    00360
    00361 template <class T>
    00362 bool MittelmannParaCntrlBase<T>::
    00363 eval_f(Index n, const Number* x,
    00364        bool new_x, Number& obj_value)
    00365 {
    00366   // Deviation of y from target
    00367   Number a = x[y_index(0,Nt_)] - y_T_[0];
    00368   Number sum = 0.5*a*a;
    00369   for (Index jx=1; jx<Nx_; jx++) {
    00370     a = x[y_index(jx,Nt_)] - y_T_[jx];
    00371     sum += a*a;
    00372   }
    00373   a = x[y_index(Nx_,Nt_)] - y_T_[Nx_];
    00374   sum += 0.5*a*a;
    00375   obj_value = 0.5*dx_*sum;
    00376
    00377   // smoothing for control
    00378   if (alpha_!=.0) {
    00379     sum = 0.5*x[u_index(Nt_)]*x[u_index(Nt_)];
    00380     for (Index it=1; it < Nt_; it++) {
    00381       sum += x[u_index(it)]*x[u_index(it)];
    00382     }
    00383     obj_value += 0.5*alpha_*dt_*sum;
    00384   }
    00385
    00386   // third term
    00387   sum = 0.;
    00388   for (Index it=1; it<Nt_; it++) {
    00389     sum += a_y_[it-1]*x[y_index(Nx_,it)] + a_u_[it-1]*x[u_index(it)];
    00390   }
    00391   sum += 0.5*(a_y_[Nt_-1]*x[y_index(Nx_,Nt_)] + a_u_[Nt_-1]*x[u_index(Nt_)]);
    00392   obj_value += dt_*sum;
    00393
    00394   return true;
    00395 }
    00396
    00397 template <class T>
    00398 bool MittelmannParaCntrlBase<T>::
    00399 eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
    00400 {
    00401   // First set all y entries to zero
    00402   for (Index jx=0; jx<=Nx_; jx++) {
    00403     for (Index it=0; it<=Nt_; it++) {
    00404       grad_f[y_index(jx,it)] = 0.;
    00405     }
    00406   }
    00407
    00408   // y entries from y target
    00409   grad_f[y_index(0,Nt_)] = 0.5*dx_*(x[y_index(0,Nt_)] - y_T_[0]);
    00410   for (Index jx=1; jx<Nx_; jx++) {
    00411     grad_f[y_index(jx,Nt_)] = dx_*(x[y_index(jx,Nt_)] - y_T_[jx]);
    00412   }
    00413   grad_f[y_index(Nx_,Nt_)] = 0.5*dx_*(x[y_index(Nx_,Nt_)] - y_T_[Nx_]);
    00414
    00415   // y entries from thrid term
    00416   for (Index it=1; it<Nt_; it++) {
    00417     grad_f[y_index(Nx_,it)] = dt_*a_y_[it-1];
    00418   }
    00419   grad_f[y_index(Nx_,Nt_)] += 0.5*dt_*a_y_[Nt_-1];
    00420
    00421   // u entries from smoothing and third term
    00422   for (Index it=1; it<Nt_; it++) {
    00423     grad_f[u_index(it)] = alpha_*dt_*x[u_index(it)] + dt_*a_u_[it-1];
    00424   }
    00425   grad_f[u_index(Nt_)] = 0.5*dt_*(alpha_*x[u_index(Nt_)] + a_u_[Nt_-1]);
    00426
    00427   return true;
    00428 }
    00429
    00430 template <class T>
    00431 bool MittelmannParaCntrlBase<T>::
    00432 eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
    00433 {
    00434   typename T::ProblemSpecs p;
    00435
    00436   Index ig=0;
    00437
    00438   Number f = 1./(2.*dx_*dx_);
    00439   for (Index jx=1; jx<Nx_; jx++) {
    00440     for (Index it=0; it<Nt_; it++) {
    00441       g[ig] = (x[y_index(jx,it)]-x[y_index(jx,it+1)])/dt_
    00442               + f*(x[y_index(jx-1,it)] - 2.*x[y_index(jx,it)]
    00443                    + x[y_index(jx+1,it)] + x[y_index(jx-1,it+1)]
    00444                    - 2.*x[y_index(jx,it+1)] + x[y_index(jx+1,it+1)]);
    00445       ig++;
    00446     }
    00447   }
    00448
    00449   for (Index it=1; it<=Nt_; it++) {
    00450     g[ig] = (x[y_index(2,it)] - 4.*x[y_index(1,it)] + 3.*x[y_index(0,it)]);
    00451     ig++;
    00452   }
    00453
    00454   f = 1./(2.*dx_);
    00455   for (Index it=1; it<=Nt_; it++) {
    00456     g[ig] =
    00457       f*(x[y_index(Nx_-2,it)] - 4.*x[y_index(Nx_-1,it)]
    00458          + 3.*x[y_index(Nx_,it)]) + beta_*x[y_index(Nx_,it)]
    00459       - x[u_index(it)] + p.phi(x[y_index(Nx_,it)]);
    00460     ig++;
    00461   }
    00462
    00463   DBG_ASSERT(ig == m);
    00464
    00465   return true;
    00466 }
    00467
    00468 template <class T>
    00469 bool MittelmannParaCntrlBase<T>::
    00470 eval_jac_g(Index n, const Number* x, bool new_x,
    00471            Index m, Index nele_jac, Index* iRow, Index *jCol,
    00472            Number* values)
    00473 {
    00474   typename T::ProblemSpecs p;
    00475
    00476   Index ijac = 0;
    00477
    00478   if (values == NULL) {
    00479     Index ig = 0;
    00480     for (Index jx=1; jx<Nx_; jx++) {
    00481       for (Index it=0; it<Nt_; it++) {
    00482         iRow[ijac] = ig;
    00483         jCol[ijac] = y_index(jx-1,it);
    00484         ijac++;
    00485         iRow[ijac] = ig;
    00486         jCol[ijac] = y_index(jx,it);
    00487         ijac++;
    00488         iRow[ijac] = ig;
    00489         jCol[ijac] = y_index(jx+1,it);
    00490         ijac++;
    00491         iRow[ijac] = ig;
    00492         jCol[ijac] = y_index(jx-1,it+1);
    00493         ijac++;
    00494         iRow[ijac] = ig;
    00495         jCol[ijac] = y_index(jx,it+1);
    00496         ijac++;
    00497         iRow[ijac] = ig;
    00498         jCol[ijac] = y_index(jx+1,it+1);
    00499         ijac++;
    00500
    00501         ig++;
    00502       }
    00503     }
    00504
    00505     for (Index it=1; it<=Nt_; it++) {
    00506       iRow[ijac] = ig;
    00507       jCol[ijac] = y_index(0,it);
    00508       ijac++;
    00509       iRow[ijac] = ig;
    00510       jCol[ijac] = y_index(1,it);
    00511       ijac++;
    00512       iRow[ijac] = ig;
    00513       jCol[ijac] = y_index(2,it);
    00514       ijac++;
    00515
    00516       ig++;
    00517     }
    00518
    00519     for (Index it=1; it<=Nt_; it++) {
    00520       iRow[ijac] = ig;
    00521       jCol[ijac] = y_index(Nx_-2,it);
    00522       ijac++;
    00523       iRow[ijac] = ig;
    00524       jCol[ijac] = y_index(Nx_-1,it);
    00525       ijac++;
    00526       iRow[ijac] = ig;
    00527       jCol[ijac] = y_index(Nx_,it);
    00528       ijac++;
    00529       iRow[ijac] = ig;
    00530       jCol[ijac] = u_index(it);
    00531       ijac++;
    00532
    00533       ig++;
    00534     }
    00535     DBG_ASSERT(ig == m);
    00536   }
    00537   else {
    00538     Number f = 1./(2.*dx_*dx_);
    00539     Number f2 = 1./dt_;
    00540     for (Index jx=1; jx<Nx_; jx++) {
    00541       for (Index it=0; it<Nt_; it++) {
    00542         values[ijac++] = f;
    00543         values[ijac++] = f2 - 2.*f;
    00544         values[ijac++] = f;
    00545         values[ijac++] = f;
    00546         values[ijac++] = -f2 - 2.*f;
    00547         values[ijac++] = f;
    00548       }
    00549     }
    00550
    00551     for (Index it=1; it<=Nt_; it++) {
    00552       values[ijac++] = 3.;
    00553       values[ijac++] = -4.;
    00554       values[ijac++] = 1.;
    00555     }
    00556
    00557     f = 1./(2.*dx_);
    00558     for (Index it=1; it<=Nt_; it++) {
    00559       values[ijac++] = f;
    00560       values[ijac++] = -4.*f;
    00561       values[ijac++] = 3.*f + beta_ + p.phi_dy(x[y_index(Nx_,it)]);
    00562       values[ijac++] = -1.;
    00563     }
    00564
    00565   }
    00566
    00567   return true;
    00568 }
    00569
    00570 template <class T>
    00571 bool MittelmannParaCntrlBase<T>::
    00572 eval_h(Index n, const Number* x, bool new_x,
    00573        Number obj_factor, Index m, const Number* lambda,
    00574        bool new_lambda, Index nele_hess, Index* iRow,
    00575        Index* jCol, Number* values)
    00576 {
    00577   typename T::ProblemSpecs p;
    00578
    00579   Index ihes = 0;
    00580
    00581   if (values == NULL) {
    00582     // y values from objective
    00583     for (Index jx=0; jx<= Nx_; jx++) {
    00584       iRow[ihes] = y_index(jx,Nt_);
    00585       jCol[ihes] = y_index(jx,Nt_);
    00586       ihes++;
    00587     }
    00588     // u from objective
    00589     for (Index it=1; it<=Nt_; it++) {
    00590       iRow[ihes] = u_index(it);
    00591       jCol[ihes] = u_index(it);
    00592       ihes++;
    00593     }
    00594
    00595     // constraint
    00596     if (!p.phi_dydy_always_zero()) {
    00597       for (Index it=1; it<=Nt_; it++) {
    00598         iRow[ihes] = y_index(Nx_,it);
    00599         jCol[ihes] = y_index(Nx_,it);
    00600         ihes++;
    00601       }
    00602     }
    00603   }
    00604   else {
    00605     // y values from objective
    00606     values[ihes++] = obj_factor*0.5*dx_;
    00607     for (Index jx=1; jx<Nx_; jx++) {
    00608       values[ihes++] = obj_factor*dx_;
    00609     }
    00610     values[ihes++] = obj_factor*0.5*dx_;
    00611     // u from objective
    00612     for (Index it=1; it<Nt_; it++) {
    00613       values[ihes++] = obj_factor*alpha_*dt_;
    00614     }
    00615     values[ihes++] = obj_factor*0.5*alpha_*dt_;
    00616
    00617     // constrainT
    00618     if (!p.phi_dydy_always_zero()) {
    00619       Index ig = (Nx_-1)*Nt_ + Nt_;
    00620       for (Index it=1; it<=Nt_; it++) {
    00621         values[ihes++] = lambda[ig++]*p.phi_dydy(x[y_index(Nx_,it)]);
    00622       }
    00623     }
    00624   }
    00625
    00626   DBG_ASSERT(ihes==nele_hess);
    00627
    00628   return true;
    00629 }
    00630
    00631 template <class T>
    00632 void MittelmannParaCntrlBase<T>::
    00633 finalize_solution(SolverReturn status,
    00634                   Index n, const Number* x, const Number* z_L,
    00635                   const Number* z_U,
    00636                   Index m, const Number* g, const Number* lambda,
    00637                   Number obj_value,
    00638           const IpoptData* ip_data,
    00639           IpoptCalculatedQuantities* ip_cq)
    00640 {}
    00641
    00642 class MittelmannParaCntrl5_1
    00643 {
    00644 public:
    00645   class ProblemSpecs
    00646   {
    00647   public:
    00648     ProblemSpecs ()
    00649         :
    00650         pi_(4.*atan(1.)),
    00651         exp13_(exp(1./3.)),
    00652         exp23_(exp(2./3.)),
    00653         exp1_(exp(1.)),
    00654         expm1_(exp(-1.)),
    00655         sqrt2_(sqrt(2.))
    00656     {}
    00657     Number T()
    00658     {
    00659       return 1.;
    00660     }
    00661     Number l()
    00662     {
    00663       return pi_/4.;
    00664     }
    00665     Number lb_y()
    00666     {
    00667       return -1e20;
    00668     }
    00669     Number ub_y()
    00670     {
    00671       return 1e20;
    00672     }
    00673     Number lb_u()
    00674     {
    00675       return 0.;
    00676     }
    00677     Number ub_u()
    00678     {
    00679       return 1.;
    00680     }
    00681     Number alpha()
    00682     {
    00683       return sqrt2_/2.*(exp23_-exp13_);
    00684     }
    00685     Number beta()
    00686     {
    00687       return 1.;
    00688     }
    00689     inline Number y_T(Number x)
    00690     {
    00691       return (exp1_ + expm1_)*cos(x);
    00692     }
    00693     inline Number a(Number x)
    00694     {
    00695       return cos(x);
    00696     }
    00697     inline Number a_y(Number t)
    00698     {
    00699       return -exp(-2.*t);
    00700     }
    00701     inline Number a_u(Number t)
    00702     {
    00703       return sqrt2_/2.*exp13_;
    00704     }
    00705     inline Number b(Number t)
    00706     {
    00707       return exp(-4.*t)/4.
    00708              - Min(1.,Max(0.,(exp(t)-exp13_)/(exp23_-exp13_)));
    00709     }
    00710     inline Number phi(Number y)
    00711     {
    00712       return y*pow(fabs(y),3);
    00713     }
    00714     inline Number phi_dy(Number y)
    00715     {
    00716       return 4.*pow(fabs(y),3);
    00717     }
    00718     inline Number phi_dydy(Number y)
    00719     {
    00720       return 12.*y*y;
    00721     }
    00722     inline bool phi_dydy_always_zero()
    00723     {
    00724       return false;
    00725     }
    00726   private:
    00727     const Number pi_;
    00728     const Number exp13_;
    00729     const Number exp23_;
    00730     const Number exp1_;
    00731     const Number expm1_;
    00732     const Number sqrt2_;
    00733   };
    00734 };
    00735
    00736 class MittelmannParaCntrl5_2_1
    00737 {
    00738 public:
    00739   class ProblemSpecs
    00740   {
    00741   public:
    00742     ProblemSpecs ()
    00743     {}
    00744     Number T()
    00745     {
    00746       return 1.58;
    00747     }
    00748     Number l()
    00749     {
    00750       return 1.;
    00751     }
    00752     Number lb_y()
    00753     {
    00754       return -1e20;
    00755     }
    00756     Number ub_y()
    00757     {
    00758       return 1e20;
    00759     }
    00760     Number lb_u()
    00761     {
    00762       return -1.;
    00763     }
    00764     Number ub_u()
    00765     {
    00766       return 1.;
    00767     }
    00768     Number alpha()
    00769     {
    00770       return 0.001;
    00771     }
    00772     Number beta()
    00773     {
    00774       return 1.;
    00775     }
    00776     inline Number y_T(Number x)
    00777     {
    00778       return .5*(1.-x*x);
    00779     }
    00780     inline Number a(Number x)
    00781     {
    00782       return 0.;
    00783     }
    00784     inline Number a_y(Number t)
    00785     {
    00786       return 0.;
    00787     }
    00788     inline Number a_u(Number t)
    00789     {
    00790       return 0.;
    00791     }
    00792     inline Number b(Number t)
    00793     {
    00794       return 0.;
    00795     }
    00796     inline Number phi(Number y)
    00797     {
    00798       return 0.;
    00799     }
    00800     inline Number phi_dy(Number y)
    00801     {
    00802       return 0.;
    00803     }
    00804     inline Number phi_dydy(Number y)
    00805     {
    00806       DBG_ASSERT(false);
    00807       return 0.;
    00808     }
    00809     inline bool phi_dydy_always_zero()
    00810     {
    00811       return true;
    00812     }
    00813   };
    00814 };
    00815
    00816 class MittelmannParaCntrl5_2_2
    00817 {
    00818 public:
    00819   class ProblemSpecs
    00820   {
    00821   public:
    00822     ProblemSpecs ()
    00823     {}
    00824     Number T()
    00825     {
    00826       return 1.58;
    00827     }
    00828     Number l()
    00829     {
    00830       return 1.;
    00831     }
    00832     Number lb_y()
    00833     {
    00834       return -1e20;
    00835     }
    00836     Number ub_y()
    00837     {
    00838       return 1e20;
    00839     }
    00840     Number lb_u()
    00841     {
    00842       return -1.;
    00843     }
    00844     Number ub_u()
    00845     {
    00846       return 1.;
    00847     }
    00848     Number alpha()
    00849     {
    00850       return 0.001;
    00851     }
    00852     Number beta()
    00853     {
    00854       return 0.;
    00855     }
    00856     inline Number y_T(Number x)
    00857     {
    00858       return .5*(1.-x*x);
    00859     }
    00860     inline Number a(Number x)
    00861     {
    00862       return 0.;
    00863     }
    00864     inline Number a_y(Number t)
    00865     {
    00866       return 0.;
    00867     }
    00868     inline Number a_u(Number t)
    00869     {
    00870       return 0.;
    00871     }
    00872     inline Number b(Number t)
    00873     {
    00874       return 0.;
    00875     }
    00876     inline Number phi(Number y)
    00877     {
    00878       return y*y;
    00879     }
    00880     inline Number phi_dy(Number y)
    00881     {
    00882       return 2.*y;
    00883     }
    00884     inline Number phi_dydy(Number y)
    00885     {
    00886       return 2.;
    00887     }
    00888     inline bool phi_dydy_always_zero()
    00889     {
    00890       return false;
    00891     }
    00892   };
    00893 };
    00894
    00895 class MittelmannParaCntrl5_2_3
    00896 {
    00897 public:
    00898   class ProblemSpecs
    00899   {
    00900   public:
    00901     ProblemSpecs ()
    00902     {}
    00903     Number T()
    00904     {
    00905       return 1.58;
    00906     }
    00907     Number l()
    00908     {
    00909       return 1.;
    00910     }
    00911     Number lb_y()
    00912     {
    00913       return 0.;
    00914     }
    00915     Number ub_y()
    00916     {
    00917       return 0.675;
    00918     }
    00919     Number lb_u()
    00920     {
    00921       return -1.;
    00922     }
    00923     Number ub_u()
    00924     {
    00925       return 1.;
    00926     }
    00927     Number alpha()
    00928     {
    00929       return 0.001;
    00930     }
    00931     Number beta()
    00932     {
    00933       return 0.;
    00934     }
    00935     inline Number y_T(Number x)
    00936     {
    00937       return .5*(1.-x*x);
    00938     }
    00939     inline Number a(Number x)
    00940     {
    00941       return 0.;
    00942     }
    00943     inline Number a_y(Number t)
    00944     {
    00945       return 0.;
    00946     }
    00947     inline Number a_u(Number t)
    00948     {
    00949       return 0.;
    00950     }
    00951     inline Number b(Number t)
    00952     {
    00953       return 0.;
    00954     }
    00955     inline Number phi(Number y)
    00956     {
    00957       return y*y;
    00958     }
    00959     inline Number phi_dy(Number y)
    00960     {
    00961       return 2.*y;
    00962     }
    00963     inline Number phi_dydy(Number y)
    00964     {
    00965       return 2.;
    00966     }
    00967     inline bool phi_dydy_always_zero()
    00968     {
    00969       return false;
    00970     }
    00971   };
    00972 };
    00973
    00974 class MittelmannParaCntrl5_try
    00975 {
    00976 public:
    00977   class ProblemSpecs
    00978   {
    00979   public:
    00980     ProblemSpecs ()
    00981         :
    00982         pi_(4.*atan(1.)),
    00983         exp13_(exp(1./3.)),
    00984         exp23_(exp(2./3.)),
    00985         exp1_(exp(1.)),
    00986         expm1_(exp(-1.)),
    00987         sqrt2_(sqrt(2.))
    00988     {}
    00989     Number T()
    00990     {
    00991       return 1.;
    00992     }
    00993     Number l()
    00994     {
    00995       return pi_/4.;
    00996     }
    00997     Number lb_y()
    00998     {
    00999       return -1e20;
    01000     }
    01001     Number ub_y()
    01002     {
    01003       return 1e20;
    01004     }
    01005     Number lb_u()
    01006     {
    01007       return 0.;
    01008     }
    01009     Number ub_u()
    01010     {
    01011       return 1.;
    01012     }
    01013     Number alpha()
    01014     {
    01015       return sqrt2_/2.*(exp23_-exp13_);
    01016     }
    01017     Number beta()
    01018     {
    01019       return 1.;
    01020     }
    01021     inline Number y_T(Number x)
    01022     {
    01023       return (exp1_ + expm1_)*cos(x);
    01024     }
    01025     inline Number a(Number x)
    01026     {
    01027       return cos(x);
    01028     }
    01029     inline Number a_y(Number t)
    01030     {
    01031       return -exp(-2.*t);
    01032     }
    01033     inline Number a_u(Number t)
    01034     {
    01035       return sqrt2_/2.*exp13_;
    01036     }
    01037     inline Number b(Number t)
    01038     {
    01039       return exp(-4.*t)/4.
    01040              - Min(1.,Max(0.,(exp(t)-exp13_)/(exp23_-exp13_)));
    01041     }
    01042     inline Number phi(Number y)
    01043     {
    01044       return -y*sin(y/10.);
    01045     }
    01046     inline Number phi_dy(Number y)
    01047     {
    01048       return -y*cos(y/10.)/10. - sin(y/10.);
    01049     }
    01050     inline Number phi_dydy(Number y)
    01051     {
    01052       return y*sin(y/10.)/100.;
    01053     }
    01054     inline bool phi_dydy_always_zero()
    01055     {
    01056       return false;
    01057     }
    01058   private:
    01059     const Number pi_;
    01060     const Number exp13_;
    01061     const Number exp23_;
    01062     const Number exp1_;
    01063     const Number expm1_;
    01064     const Number sqrt2_;
    01065   };
    01066 };
    01067
    01068 #endif

All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines