00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #ifndef __MITTELMANNDISTRCNTRLNEUMA_HPP__
00011 #define __MITTELMANNDISTRCNTRLNEUMA_HPP__
00012
00013 #include "IpTNLP.hpp"
00014 #include "RegisteredTNLP.hpp"
00015
00016 #ifdef HAVE_CMATH
00017 # include <cmath>
00018 #else
00019 # ifdef HAVE_MATH_H
00020 # include <math.h>
00021 # else
00022 # error "don't have header file for math"
00023 # endif
00024 #endif
00025
00026 #include <cstdio>
00027
00028 using namespace Ipopt;
00029
00036 class MittelmannDistCntrlNeumABase : public RegisteredTNLP
00037 {
00038 public:
00041 MittelmannDistCntrlNeumABase();
00042
00044 virtual ~MittelmannDistCntrlNeumABase();
00045
00049 virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00050 Index& nnz_h_lag, IndexStyleEnum& index_style);
00051
00053 virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
00054 Index m, Number* g_l, Number* g_u);
00055
00057 virtual bool get_starting_point(Index n, bool init_x, Number* x,
00058 bool init_z, Number* z_L, Number* z_U,
00059 Index m, bool init_lambda,
00060 Number* lambda);
00061
00063 virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
00064
00066 virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
00067
00069 virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
00070
00075 virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
00076 Index m, Index nele_jac, Index* iRow, Index *jCol,
00077 Number* values);
00078
00083 virtual bool eval_h(Index n, const Number* x, bool new_x,
00084 Number obj_factor, Index m, const Number* lambda,
00085 bool new_lambda, Index nele_hess, Index* iRow,
00086 Index* jCol, Number* values);
00087
00089
00091 virtual bool get_scaling_parameters(Number& obj_scaling,
00092 bool& use_x_scaling, Index n,
00093 Number* x_scaling,
00094 bool& use_g_scaling, Index m,
00095 Number* g_scaling);
00096
00101 virtual void finalize_solution(SolverReturn status,
00102 Index n, const Number* x, const Number* z_L, const Number* z_U,
00103 Index m, const Number* g, const Number* lambda,
00104 Number obj_value,
00105 const IpoptData* ip_data,
00106 IpoptCalculatedQuantities* ip_cq);
00108
00109 protected:
00113 void SetBaseParameters(Index N, Number lb_y,
00114 Number ub_y, Number lb_u, Number ub_u,
00115 Number b_0j, Number b_1j, Number b_i0, Number b_i1,
00116 Number u_init);
00117
00121 virtual Number y_d_cont(Number x1, Number x2) const =0;
00123 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const =0;
00125 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const =0;
00127 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const =0;
00129 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const =0;
00132 virtual bool fint_cont_dydy_alwayszero() const =0;
00134 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const =0;
00137 virtual bool fint_cont_dudu_alwayszero() const =0;
00139 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const =0;
00142 virtual bool fint_cont_dydu_alwayszero() const =0;
00144 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const =0;
00146 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const =0;
00148 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const =0;
00150 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const =0;
00153 virtual bool d_cont_dydy_alwayszero() const =0;
00155 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const =0;
00158 virtual bool d_cont_dudu_alwayszero() const =0;
00160 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const =0;
00163 virtual bool d_cont_dydu_alwayszero() const =0;
00165
00166 private:
00178 MittelmannDistCntrlNeumABase(const MittelmannDistCntrlNeumABase&);
00179 MittelmannDistCntrlNeumABase& operator=(const MittelmannDistCntrlNeumABase&);
00181
00185 Index N_;
00187 Number h_;
00189 Number hh_;
00191 Number lb_y_;
00193 Number ub_y_;
00195 Number lb_u_;
00197 Number ub_u_;
00200 Number b_0j_;
00203 Number b_1j_;
00206 Number b_i0_;
00209 Number b_i1_;
00211 Number u_init_;
00213 Number* y_d_;
00215
00220 inline Index y_index(Index i, Index j) const
00221 {
00222 return j + (N_+2)*i;
00223 }
00226 inline Index u_index(Index i, Index j) const
00227 {
00228 return (N_+2)*(N_+2) + (j-1) + (N_)*(i-1);
00229 }
00232 inline Index pde_index(Index i, Index j) const
00233 {
00234 return (j-1) + N_*(i-1);
00235 }
00237 inline Number x1_grid(Index i) const
00238 {
00239 return h_*(Number)i;
00240 }
00242 inline Number x2_grid(Index i) const
00243 {
00244 return h_*(Number)i;
00245 }
00247 };
00248
00250 class MittelmannDistCntrlNeumA1 : public MittelmannDistCntrlNeumABase
00251 {
00252 public:
00253 MittelmannDistCntrlNeumA1()
00254 :
00255 pi_(4.*atan(1.)),
00256 alpha_(0.001)
00257 {}
00258
00259 virtual ~MittelmannDistCntrlNeumA1()
00260 {}
00261
00262 virtual bool InitializeProblem(Index N)
00263 {
00264 if (N<1) {
00265 printf("N has to be at least 1.");
00266 return false;
00267 }
00268 Number lb_y = -1e20;
00269 Number ub_y = 0.371;
00270 Number lb_u = -8.;
00271 Number ub_u = 9.;
00272 Number b_0j = 1.;
00273 Number b_1j = 1.;
00274 Number b_i0 = 1.;
00275 Number b_i1 = 1.;
00276 Number u_init = (ub_u+lb_u)/2.;
00277
00278 SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init);
00279 return true;
00280 }
00281 protected:
00283 virtual Number y_d_cont(Number x1, Number x2) const
00284 {
00285 return sin(2.*pi_*x1)*sin(2.*pi_*x2);
00286 }
00288 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const
00289 {
00290 Number diff_y = y-y_d_cont(x1,x2);
00291 return 0.5*(diff_y*diff_y + alpha_*u*u);
00292 }
00294 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const
00295 {
00296 return y-y_d_cont(x1,x2);
00297 }
00298
00300 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const
00301 {
00302 return alpha_*u;
00303 }
00305 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const
00306 {
00307 return 1.;
00308 }
00311 virtual bool fint_cont_dydy_alwayszero() const
00312 {
00313 return false;
00314 }
00316 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const
00317 {
00318 return alpha_;
00319 }
00322 virtual bool fint_cont_dudu_alwayszero() const
00323 {
00324 return false;
00325 }
00327 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const
00328 {
00329 return 0.;
00330 }
00333 virtual bool fint_cont_dydu_alwayszero() const
00334 {
00335 return true;
00336 }
00338 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
00339 {
00340 return -exp(y) - u;
00341 }
00343 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
00344 {
00345 return -exp(y);
00346 }
00348 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
00349 {
00350 return -1.;
00351 }
00353 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
00354 {
00355 return -exp(y);
00356 }
00359 virtual bool d_cont_dydy_alwayszero() const
00360 {
00361 return false;
00362 }
00364 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const
00365 {
00366 return 0.;
00367 }
00370 virtual bool d_cont_dudu_alwayszero() const
00371 {
00372 return true;
00373 }
00375 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const
00376 {
00377 return 0.;
00378 }
00381 virtual bool d_cont_dydu_alwayszero() const
00382 {
00383 return true;
00384 }
00385 private:
00388 MittelmannDistCntrlNeumA1(const MittelmannDistCntrlNeumA1&);
00389 MittelmannDistCntrlNeumA1& operator=(const MittelmannDistCntrlNeumA1&);
00391
00392 const Number pi_;
00394 const Number alpha_;
00395 };
00396
00398 class MittelmannDistCntrlNeumA2 : public MittelmannDistCntrlNeumABase
00399 {
00400 public:
00401 MittelmannDistCntrlNeumA2()
00402 :
00403 pi_(4.*atan(1.))
00404 {}
00405
00406 virtual ~MittelmannDistCntrlNeumA2()
00407 {}
00408
00409 virtual bool InitializeProblem(Index N)
00410 {
00411 if (N<1) {
00412 printf("N has to be at least 1.");
00413 return false;
00414 }
00415 Number lb_y = -1e20;
00416 Number ub_y = 0.371;
00417 Number lb_u = -8.;
00418 Number ub_u = 9.;
00419 Number b_0j = 1.;
00420 Number b_1j = 1.;
00421 Number b_i0 = 1.;
00422 Number b_i1 = 1.;
00423 Number u_init = (ub_u+lb_u)/2.;
00424
00425 SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init);
00426 return true;
00427 }
00428 protected:
00430 virtual Number y_d_cont(Number x1, Number x2) const
00431 {
00432 return sin(2.*pi_*x1)*sin(2.*pi_*x2);
00433 }
00435 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const
00436 {
00437 Number diff_y = y-y_d_cont(x1,x2);
00438 return 0.5*diff_y*diff_y;
00439 }
00441 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const
00442 {
00443 return y-y_d_cont(x1,x2);
00444 }
00445
00447 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const
00448 {
00449 return 0.;
00450 }
00452 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const
00453 {
00454 return 1.;
00455 }
00458 virtual bool fint_cont_dydy_alwayszero() const
00459 {
00460 return false;
00461 }
00463 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const
00464 {
00465 return 0.;
00466 }
00469 virtual bool fint_cont_dudu_alwayszero() const
00470 {
00471 return true;
00472 }
00474 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const
00475 {
00476 return 0.;
00477 }
00480 virtual bool fint_cont_dydu_alwayszero() const
00481 {
00482 return true;
00483 }
00485 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
00486 {
00487 return -exp(y) - u;
00488 }
00490 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
00491 {
00492 return -exp(y);
00493 }
00495 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
00496 {
00497 return -1.;
00498 }
00500 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
00501 {
00502 return -exp(y);
00503 }
00506 virtual bool d_cont_dydy_alwayszero() const
00507 {
00508 return false;
00509 }
00511 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const
00512 {
00513 return 0.;
00514 }
00517 virtual bool d_cont_dudu_alwayszero() const
00518 {
00519 return true;
00520 }
00522 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const
00523 {
00524 return 0.;
00525 }
00528 virtual bool d_cont_dydu_alwayszero() const
00529 {
00530 return true;
00531 }
00532 private:
00535 MittelmannDistCntrlNeumA2(const MittelmannDistCntrlNeumA2&);
00536 MittelmannDistCntrlNeumA2& operator=(const MittelmannDistCntrlNeumA2&);
00538
00539 const Number pi_;
00540 };
00541
00543 class MittelmannDistCntrlNeumA3 : public MittelmannDistCntrlNeumABase
00544 {
00545 public:
00546 MittelmannDistCntrlNeumA3()
00547 :
00548 pi_(4.*atan(1.)),
00549 M_(1.),
00550 K_(0.8),
00551 b_(1.)
00552 {}
00553
00554 virtual ~MittelmannDistCntrlNeumA3()
00555 {}
00556
00557 virtual bool InitializeProblem(Index N)
00558 {
00559 if (N<1) {
00560 printf("N has to be at least 1.");
00561 return false;
00562 }
00563 Number lb_y = 3.;
00564 Number ub_y = 6.09;
00565 Number lb_u = 1.4;
00566 Number ub_u = 1.6;
00567 Number b_0j = 1.;
00568 Number b_1j = 0.;
00569 Number b_i0 = 1.;
00570 Number b_i1 = 0.;
00571 Number u_init = (ub_u+lb_u)/2.;
00572
00573 SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init);
00574 return true;
00575 }
00576 protected:
00578 virtual Number y_d_cont(Number x1, Number x2) const
00579 {
00580 return 6.;
00581 }
00583 virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const
00584 {
00585 return u*(M_*u - K_*y);
00586 }
00588 virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const
00589 {
00590 return -K_*u;
00591 }
00592
00594 virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const
00595 {
00596 return 2.*M_*u - K_*y;
00597 }
00599 virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const
00600 {
00601 return 0.;
00602 }
00605 virtual bool fint_cont_dydy_alwayszero() const
00606 {
00607 return true;
00608 }
00610 virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const
00611 {
00612 return 2.*M_;
00613 }
00616 virtual bool fint_cont_dudu_alwayszero() const
00617 {
00618 return false;
00619 }
00621 virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const
00622 {
00623 return -K_;
00624 }
00627 virtual bool fint_cont_dydu_alwayszero() const
00628 {
00629 return false;
00630 }
00632 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
00633 {
00634 return y*(u + b_*y - a(x1,x2));
00635 }
00637 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
00638 {
00639 return (u + 2.*b_*y -a(x1,x2));
00640 }
00642 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
00643 {
00644 return y;
00645 }
00647 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
00648 {
00649 return 2.*b_;
00650 }
00653 virtual bool d_cont_dydy_alwayszero() const
00654 {
00655 return false;
00656 }
00658 virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const
00659 {
00660 return 0.;
00661 }
00664 virtual bool d_cont_dudu_alwayszero() const
00665 {
00666 return true;
00667 }
00669 virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const
00670 {
00671 return 1.;
00672 }
00675 virtual bool d_cont_dydu_alwayszero() const
00676 {
00677 return false;
00678 }
00679 private:
00682 MittelmannDistCntrlNeumA3(const MittelmannDistCntrlNeumA3&);
00683 MittelmannDistCntrlNeumA3& operator=(const MittelmannDistCntrlNeumA3&);
00685
00686 const Number pi_;
00687
00689 const Number M_;
00690 const Number K_;
00691 const Number b_;
00693
00694 inline Number a(Number x1, Number x2) const
00695 {
00696 return 7. + 4.*sin(2.*pi_*x1*x2);
00697 }
00698 };
00699
00700 #endif