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