- Main Page
- Namespaces
- Classes
- Files

- File List
- File Members

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