- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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