00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #ifndef __MITTELMANNBNDRYCNTRLNEUM_HPP__
00011 #define __MITTELMANNBNDRYCNTRLNEUM_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
00035 class MittelmannBndryCntrlNeumBase : public RegisteredTNLP
00036 {
00037 public:
00040 MittelmannBndryCntrlNeumBase();
00041
00043 virtual ~MittelmannBndryCntrlNeumBase();
00044
00048 virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00049 Index& nnz_h_lag, IndexStyleEnum& index_style);
00050
00052 virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
00053 Index m, Number* g_l, Number* g_u);
00054
00056 virtual bool get_starting_point(Index n, bool init_x, Number* x,
00057 bool init_z, Number* z_L, Number* z_U,
00058 Index m, bool init_lambda,
00059 Number* lambda);
00060
00062 virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
00063
00065 virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
00066
00068 virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
00069
00074 virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
00075 Index m, Index nele_jac, Index* iRow, Index *jCol,
00076 Number* values);
00077
00082 virtual bool eval_h(Index n, const Number* x, bool new_x,
00083 Number obj_factor, Index m, const Number* lambda,
00084 bool new_lambda, Index nele_hess, Index* iRow,
00085 Index* jCol, Number* values);
00086
00088
00090 virtual bool get_scaling_parameters(Number& obj_scaling,
00091 bool& use_x_scaling, Index n,
00092 Number* x_scaling,
00093 bool& use_g_scaling, Index m,
00094 Number* g_scaling);
00095
00100 virtual void finalize_solution(SolverReturn status,
00101 Index n, const Number* x, const Number* z_L, const Number* z_U,
00102 Index m, const Number* g, const Number* lambda,
00103 Number obj_value,
00104 const IpoptData* ip_data,
00105 IpoptCalculatedQuantities* ip_cq);
00107
00108 protected:
00112 void SetBaseParameters(Index N, Number alpha, Number lb_y,
00113 Number ub_y, Number lb_u, Number ub_u,
00114 Number u_init);
00115
00119 virtual Number y_d_cont(Number x1, Number x2) const =0;
00121 virtual Number d_cont(Number x1, Number x2, Number y) const =0;
00123 virtual Number d_cont_dy(Number x1, Number x2, Number y) const =0;
00125 virtual Number d_cont_dydy(Number x1, Number x2, Number y) const =0;
00128 virtual bool d_cont_dydy_alwayszero() const =0;
00130 virtual Number b_cont(Number x1, Number x2, Number y, Number u) const =0;
00132 virtual Number b_cont_dy(Number x1, Number x2, Number y, Number u) const =0;
00134 virtual Number b_cont_du(Number x1, Number x2, Number y, Number u) const =0;
00136 virtual Number b_cont_dydy(Number x1, Number x2, Number y, Number u) const =0;
00139 virtual bool b_cont_dydy_alwayszero() const =0;
00141
00142 private:
00154 MittelmannBndryCntrlNeumBase(const MittelmannBndryCntrlNeumBase&);
00155 MittelmannBndryCntrlNeumBase& operator=(const MittelmannBndryCntrlNeumBase&);
00157
00161 Index N_;
00163 Number h_;
00165 Number hh_;
00167 Number lb_y_;
00169 Number ub_y_;
00171 Number lb_u_;
00173 Number ub_u_;
00175 Number u_init_;
00178 Number alpha_;
00180 Number* y_d_;
00182
00187 inline Index y_index(Index i, Index j) const
00188 {
00189 return j + (N_+2)*i;
00190 }
00193 inline Index u0j_index(Index j) const
00194 {
00195 return (N_+2)*(N_+2) + j-1;
00196 }
00199 inline Index u1j_index(Index j) const
00200 {
00201 return (N_+2)*(N_+2) + N_ + j-1;
00202 }
00205 inline Index ui0_index(Index j) const
00206 {
00207 return (N_+2)*(N_+2) + 2*N_ + j-1;
00208 }
00211 inline Index ui1_index(Index j) const
00212 {
00213 return (N_+2)*(N_+2) + 3*N_ + j-1;
00214 }
00216 inline Number x1_grid(Index i) const
00217 {
00218 return h_*(Number)i;
00219 }
00221 inline Number x2_grid(Index j) const
00222 {
00223 return h_*(Number)j;
00224 }
00226 };
00227
00229 class MittelmannBndryCntrlNeum1 : public MittelmannBndryCntrlNeumBase
00230 {
00231 public:
00232 MittelmannBndryCntrlNeum1()
00233 {}
00234
00235 virtual ~MittelmannBndryCntrlNeum1()
00236 {}
00237
00238 virtual bool InitializeProblem(Index N)
00239 {
00240 if (N<1) {
00241 printf("N has to be at least 1.");
00242 return false;
00243 }
00244 Number alpha = 0.01;
00245 Number lb_y = -1e20;
00246 Number ub_y = 2.071;
00247 Number lb_u = 3.7;
00248 Number ub_u = 4.5;
00249 Number u_init = (ub_u+lb_u)/2.;
00250
00251 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00252 return true;
00253 }
00254 protected:
00256 virtual Number y_d_cont(Number x1, Number x2) const
00257 {
00258 return 2. - 2.*(x1*(x1-1.) + x2*(x2-1.));
00259 }
00261 virtual Number d_cont(Number x1, Number x2, Number y) const
00262 {
00263 return 0.;
00264 }
00266 virtual Number d_cont_dy(Number x1, Number x2, Number y) const
00267 {
00268 return 0.;
00269 }
00271 virtual Number d_cont_dydy(Number x1, Number x2, Number y) const
00272 {
00273 return 0.;
00274 }
00277 virtual bool d_cont_dydy_alwayszero() const
00278 {
00279 return true;
00280 }
00282 virtual Number b_cont(Number x1, Number x2, Number y, Number u) const
00283 {
00284 return u - y*y;
00285 }
00287 virtual Number b_cont_dy(Number x1, Number x2, Number y, Number u) const
00288 {
00289 return - 2.*y;
00290 }
00292 virtual Number b_cont_du(Number x1, Number x2, Number y, Number u) const
00293 {
00294 return 1.;
00295 }
00297 virtual Number b_cont_dydy(Number x1, Number x2, Number y, Number u) const
00298 {
00299 return -2.;
00300 }
00303 virtual bool b_cont_dydy_alwayszero() const
00304 {
00305 return false;
00306 }
00307 private:
00310 MittelmannBndryCntrlNeum1(const MittelmannBndryCntrlNeum1&);
00311 MittelmannBndryCntrlNeum1& operator=(const MittelmannBndryCntrlNeum1&);
00313 };
00314
00316 class MittelmannBndryCntrlNeum2 : public MittelmannBndryCntrlNeumBase
00317 {
00318 public:
00319 MittelmannBndryCntrlNeum2()
00320 {}
00321
00322 virtual ~MittelmannBndryCntrlNeum2()
00323 {}
00324
00325 virtual bool InitializeProblem(Index N)
00326 {
00327 if (N<1) {
00328 printf("N has to be at least 1.");
00329 return false;
00330 }
00331 Number alpha = 0.;
00332 Number lb_y = -1e20;
00333 Number ub_y = 2.835;
00334 Number lb_u = 6.;
00335 Number ub_u = 9.;
00336 Number u_init = (ub_u+lb_u)/2.;
00337
00338 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00339 return true;
00340 }
00341 protected:
00343 virtual Number y_d_cont(Number x1, Number x2) const
00344 {
00345 return 2. - 2.*(x1*(x1-1.) + x2*(x2-1.));
00346 }
00348 virtual Number d_cont(Number x1, Number x2, Number y) const
00349 {
00350 return 0.;
00351 }
00353 virtual Number d_cont_dy(Number x1, Number x2, Number y) const
00354 {
00355 return 0.;
00356 }
00358 virtual Number d_cont_dydy(Number x1, Number x2, Number y) const
00359 {
00360 return 0.;
00361 }
00364 virtual bool d_cont_dydy_alwayszero() const
00365 {
00366 return true;
00367 }
00369 virtual Number b_cont(Number x1, Number x2, Number y, Number u) const
00370 {
00371 return u - y*y;
00372 }
00374 virtual Number b_cont_dy(Number x1, Number x2, Number y, Number u) const
00375 {
00376 return - 2.*y;
00377 }
00379 virtual Number b_cont_du(Number x1, Number x2, Number y, Number u) const
00380 {
00381 return 1.;
00382 }
00384 virtual Number b_cont_dydy(Number x1, Number x2, Number y, Number u) const
00385 {
00386 return -2.;
00387 }
00390 virtual bool b_cont_dydy_alwayszero() const
00391 {
00392 return false;
00393 }
00394 private:
00397 MittelmannBndryCntrlNeum2(const MittelmannBndryCntrlNeum2&);
00398 MittelmannBndryCntrlNeum2& operator=(const MittelmannBndryCntrlNeum2&);
00400 };
00401
00403 class MittelmannBndryCntrlNeum3 : public MittelmannBndryCntrlNeumBase
00404 {
00405 public:
00406 MittelmannBndryCntrlNeum3()
00407 {}
00408
00409 virtual ~MittelmannBndryCntrlNeum3()
00410 {}
00411
00412 virtual bool InitializeProblem(Index N)
00413 {
00414 if (N<1) {
00415 printf("N has to be at least 1.");
00416 return false;
00417 }
00418 Number alpha = 0.01;
00419 Number lb_y = -1e20;
00420 Number ub_y = 2.7;
00421 Number lb_u = 1.8;
00422 Number ub_u = 2.5;
00423 Number u_init = (ub_u+lb_u)/2.;
00424
00425 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00426 return true;
00427 }
00428 protected:
00430 virtual Number y_d_cont(Number x1, Number x2) const
00431 {
00432 return 2. - 2.*(x1*(x1-1.) + x2*(x2-1.));
00433 }
00435 virtual Number d_cont(Number x1, Number x2, Number y) const
00436 {
00437 return y*y*y-y;
00438 }
00440 virtual Number d_cont_dy(Number x1, Number x2, Number y) const
00441 {
00442 return 3.*y*y-1.;
00443 }
00445 virtual Number d_cont_dydy(Number x1, Number x2, Number y) const
00446 {
00447 return 6.*y;
00448 }
00451 virtual bool d_cont_dydy_alwayszero() const
00452 {
00453 return false;
00454 }
00456 virtual Number b_cont(Number x1, Number x2, Number y, Number u) const
00457 {
00458 return u;
00459 }
00461 virtual Number b_cont_dy(Number x1, Number x2, Number y, Number u) const
00462 {
00463 return 0.;
00464 }
00466 virtual Number b_cont_du(Number x1, Number x2, Number y, Number u) const
00467 {
00468 return 1.;
00469 }
00471 virtual Number b_cont_dydy(Number x1, Number x2, Number y, Number u) const
00472 {
00473 return 0.;
00474 }
00477 virtual bool b_cont_dydy_alwayszero() const
00478 {
00479 return true;
00480 }
00481 private:
00484 MittelmannBndryCntrlNeum3(const MittelmannBndryCntrlNeum3&);
00485 MittelmannBndryCntrlNeum3& operator=(const MittelmannBndryCntrlNeum3&);
00487 };
00488
00490 class MittelmannBndryCntrlNeum4 : public MittelmannBndryCntrlNeumBase
00491 {
00492 public:
00493 MittelmannBndryCntrlNeum4()
00494 {}
00495
00496 virtual ~MittelmannBndryCntrlNeum4()
00497 {}
00498
00499 virtual bool InitializeProblem(Index N)
00500 {
00501 if (N<1) {
00502 printf("N has to be at least 1.");
00503 return false;
00504 }
00505 Number alpha = 0.;
00506 Number lb_y = -1e20;
00507 Number ub_y = 2.7;
00508 Number lb_u = 1.8;
00509 Number ub_u = 2.5;
00510 Number u_init = (ub_u+lb_u)/2.;
00511
00512 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00513 return true;
00514 }
00515 protected:
00517 virtual Number y_d_cont(Number x1, Number x2) const
00518 {
00519 return 2. - 2.*(x1*(x1-1.) + x2*(x2-1.));
00520 }
00522 virtual Number d_cont(Number x1, Number x2, Number y) const
00523 {
00524 return y*y*y-y;
00525 }
00527 virtual Number d_cont_dy(Number x1, Number x2, Number y) const
00528 {
00529 return 3.*y*y-1.;
00530 }
00532 virtual Number d_cont_dydy(Number x1, Number x2, Number y) const
00533 {
00534 return 6.*y;
00535 }
00538 virtual bool d_cont_dydy_alwayszero() const
00539 {
00540 return false;
00541 }
00543 virtual Number b_cont(Number x1, Number x2, Number y, Number u) const
00544 {
00545 return u;
00546 }
00548 virtual Number b_cont_dy(Number x1, Number x2, Number y, Number u) const
00549 {
00550 return 0.;
00551 }
00553 virtual Number b_cont_du(Number x1, Number x2, Number y, Number u) const
00554 {
00555 return 1.;
00556 }
00558 virtual Number b_cont_dydy(Number x1, Number x2, Number y, Number u) const
00559 {
00560 return 0.;
00561 }
00564 virtual bool b_cont_dydy_alwayszero() const
00565 {
00566 return true;
00567 }
00568 private:
00571 MittelmannBndryCntrlNeum4(const MittelmannBndryCntrlNeum4&);
00572 MittelmannBndryCntrlNeum4& operator=(const MittelmannBndryCntrlNeum4&);
00574 };
00575
00576 #endif