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);
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
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
00282
00283
00284
00285
00286
00287
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
00294 for (Index i=0; i<Nt_*(Nx_-1) + Nt_; i++) {
00295 g_l[i] = 0.;
00296 g_u[i] = 0.;
00297 }
00298
00299
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
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
00325 for (Index i=1; i<=Nt_; i++) {
00326 x[u_index(i)] = (ub_u_+lb_u_)/2.;
00327 }
00328
00329
00330
00331
00332
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
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
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
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
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
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
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
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
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
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
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
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
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
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