Ipopt  trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
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$
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 Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines