Ipopt
trunk
|
00001 // Copyright (C) 2005, 2006 International Business Machines and others. 00002 // All Rights Reserved. 00003 // This code is published under the Eclipse Public License. 00004 // 00005 // $Id$ 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_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 // Set overall bounds on the variables 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 // Boundary condition for t=0 00290 for (Index it=0; it<=Nt_; it++) { 00291 Index iy = y_index(0,it); 00292 x_u[iy] = x_l[iy] = p.a(t_grid(it)); 00293 } 00294 */ 00295 // Boundary condition for t=0 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 // all discretized PDE constraints have right hand side zero 00302 for (Index i=0; i<Nt_*(Nx_-1) + Nt_; i++) { 00303 g_l[i] = 0.; 00304 g_u[i] = 0.; 00305 } 00306 00307 // but we put b on the right hand side for the x=L boundary condition 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 // Set starting point for y 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 // for u 00333 for (Index i=1; i<=Nt_; i++) { 00334 x[u_index(i)] = (ub_u_+lb_u_)/2.; 00335 } 00336 00337 /* 00338 // DELETEME 00339 for (Index i=0; i<n; i++) { 00340 x[i] += 0.01*i; 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 // Deviation of y from target 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 // smoothing for control 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 // third term 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 // First set all y entries to zero 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 // y entries from y target 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 // y entries from thrid term 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 // u entries from smoothing and third term 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 // y values from objective 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 // u from objective 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 // constraint 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 // y values from objective 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 // u from objective 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 // constrainT 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