/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

Generated on Thu May 15 22:17:56 2008 by  doxygen 1.4.7