00001
00002
00003
00004
00005
00006
00007
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
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
00284
00285
00286
00287
00288
00289
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
00296 for (Index i=0; i<Nt_*(Nx_-1) + Nt_; i++) {
00297 g_l[i] = 0.;
00298 g_u[i] = 0.;
00299 }
00300
00301
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
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
00327 for (Index i=1; i<=Nt_; i++) {
00328 x[u_index(i)] = (ub_u_+lb_u_)/2.;
00329 }
00330
00331
00332
00333
00334
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
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
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
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
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
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
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
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
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
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
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
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
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
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