123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714 |
- #ifndef __MITTELMANNDISTRCNTRLNEUMA_HPP__
- #define __MITTELMANNDISTRCNTRLNEUMA_HPP__
- #include "IpTNLP.hpp"
- #include "RegisteredTNLP.hpp"
- #ifdef HAVE_CONFIG_H
- #include "config.h"
- #else
- #include "configall_system.h"
- #endif
- #ifdef HAVE_CMATH
- # include <cmath>
- #else
- # ifdef HAVE_MATH_H
- # include <math.h>
- # else
- # error "don't have header file for math"
- # endif
- #endif
- #ifdef HAVE_CSTDIO
- # include <cstdio>
- #else
- # ifdef HAVE_STDIO_H
- # include <stdio.h>
- # else
- # error "don't have header file for stdio"
- # endif
- #endif
- using namespace Ipopt;
- class MittelmannDistCntrlNeumABase : public RegisteredTNLP
- {
- public:
-
- MittelmannDistCntrlNeumABase();
-
- virtual ~MittelmannDistCntrlNeumABase();
-
-
-
- virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
- Index& nnz_h_lag, IndexStyleEnum& index_style);
-
- virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
- Index m, Number* g_l, Number* g_u);
-
- virtual bool get_starting_point(Index n, bool init_x, Number* x,
- bool init_z, Number* z_L, Number* z_U,
- Index m, bool init_lambda,
- Number* lambda);
-
- virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
-
- virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
-
- virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
-
- virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
- Index m, Index nele_jac, Index* iRow, Index *jCol,
- Number* values);
-
- virtual bool eval_h(Index n, const Number* x, bool new_x,
- Number obj_factor, Index m, const Number* lambda,
- bool new_lambda, Index nele_hess, Index* iRow,
- Index* jCol, Number* values);
-
-
- virtual bool get_scaling_parameters(Number& obj_scaling,
- bool& use_x_scaling, Index n,
- Number* x_scaling,
- bool& use_g_scaling, Index m,
- Number* g_scaling);
-
-
-
- virtual void finalize_solution(SolverReturn status,
- Index n, const Number* x, const Number* z_L, const Number* z_U,
- Index m, const Number* g, const Number* lambda,
- Number obj_value,
- const IpoptData* ip_data,
- IpoptCalculatedQuantities* ip_cq);
-
- protected:
-
- void SetBaseParameters(Index N, Number lb_y,
- Number ub_y, Number lb_u, Number ub_u,
- Number b_0j, Number b_1j, Number b_i0, Number b_i1,
- Number u_init);
-
-
-
- virtual Number y_d_cont(Number x1, Number x2) const =0;
-
- virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual bool fint_cont_dydy_alwayszero() const =0;
-
- virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual bool fint_cont_dudu_alwayszero() const =0;
-
- virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual bool fint_cont_dydu_alwayszero() const =0;
-
- virtual Number d_cont(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual bool d_cont_dydy_alwayszero() const =0;
-
- virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual bool d_cont_dudu_alwayszero() const =0;
-
- virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const =0;
-
- virtual bool d_cont_dydu_alwayszero() const =0;
-
- private:
-
-
- MittelmannDistCntrlNeumABase(const MittelmannDistCntrlNeumABase&);
- MittelmannDistCntrlNeumABase& operator=(const MittelmannDistCntrlNeumABase&);
-
-
-
-
- Index N_;
-
- Number h_;
-
- Number hh_;
-
- Number lb_y_;
-
- Number ub_y_;
-
- Number lb_u_;
-
- Number ub_u_;
-
- Number b_0j_;
-
- Number b_1j_;
-
- Number b_i0_;
-
- Number b_i1_;
-
- Number u_init_;
-
- Number* y_d_;
-
-
-
-
- inline Index y_index(Index i, Index j) const
- {
- return j + (N_+2)*i;
- }
-
- inline Index u_index(Index i, Index j) const
- {
- return (N_+2)*(N_+2) + (j-1) + (N_)*(i-1);
- }
-
- inline Index pde_index(Index i, Index j) const
- {
- return (j-1) + N_*(i-1);
- }
-
- inline Number x1_grid(Index i) const
- {
- return h_*(Number)i;
- }
-
- inline Number x2_grid(Index i) const
- {
- return h_*(Number)i;
- }
-
- };
- class MittelmannDistCntrlNeumA1 : public MittelmannDistCntrlNeumABase
- {
- public:
- MittelmannDistCntrlNeumA1()
- :
- pi_(4.*atan(1.)),
- alpha_(0.001)
- {}
- virtual ~MittelmannDistCntrlNeumA1()
- {}
- virtual bool InitializeProblem(Index N)
- {
- if (N<1) {
- printf("N has to be at least 1.");
- return false;
- }
- Number lb_y = -1e20;
- Number ub_y = 0.371;
- Number lb_u = -8.;
- Number ub_u = 9.;
- Number b_0j = 1.;
- Number b_1j = 1.;
- Number b_i0 = 1.;
- Number b_i1 = 1.;
- Number u_init = (ub_u+lb_u)/2.;
- SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init);
- return true;
- }
- protected:
-
- virtual Number y_d_cont(Number x1, Number x2) const
- {
- return sin(2.*pi_*x1)*sin(2.*pi_*x2);
- }
-
- virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const
- {
- Number diff_y = y-y_d_cont(x1,x2);
- return 0.5*(diff_y*diff_y + alpha_*u*u);
- }
-
- virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const
- {
- return y-y_d_cont(x1,x2);
- }
-
- virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const
- {
- return alpha_*u;
- }
-
- virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const
- {
- return 1.;
- }
-
- virtual bool fint_cont_dydy_alwayszero() const
- {
- return false;
- }
-
- virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const
- {
- return alpha_;
- }
-
- virtual bool fint_cont_dudu_alwayszero() const
- {
- return false;
- }
-
- virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool fint_cont_dydu_alwayszero() const
- {
- return true;
- }
-
- virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
- {
- return -exp(y) - u;
- }
-
- virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
- {
- return -exp(y);
- }
-
- virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
- {
- return -1.;
- }
-
- virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
- {
- return -exp(y);
- }
-
- virtual bool d_cont_dydy_alwayszero() const
- {
- return false;
- }
-
- virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool d_cont_dudu_alwayszero() const
- {
- return true;
- }
-
- virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool d_cont_dydu_alwayszero() const
- {
- return true;
- }
- private:
-
-
- MittelmannDistCntrlNeumA1(const MittelmannDistCntrlNeumA1&);
- MittelmannDistCntrlNeumA1& operator=(const MittelmannDistCntrlNeumA1&);
-
-
- const Number pi_;
-
- const Number alpha_;
- };
- class MittelmannDistCntrlNeumA2 : public MittelmannDistCntrlNeumABase
- {
- public:
- MittelmannDistCntrlNeumA2()
- :
- pi_(4.*atan(1.))
- {}
- virtual ~MittelmannDistCntrlNeumA2()
- {}
- virtual bool InitializeProblem(Index N)
- {
- if (N<1) {
- printf("N has to be at least 1.");
- return false;
- }
- Number lb_y = -1e20;
- Number ub_y = 0.371;
- Number lb_u = -8.;
- Number ub_u = 9.;
- Number b_0j = 1.;
- Number b_1j = 1.;
- Number b_i0 = 1.;
- Number b_i1 = 1.;
- Number u_init = (ub_u+lb_u)/2.;
- SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init);
- return true;
- }
- protected:
-
- virtual Number y_d_cont(Number x1, Number x2) const
- {
- return sin(2.*pi_*x1)*sin(2.*pi_*x2);
- }
-
- virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const
- {
- Number diff_y = y-y_d_cont(x1,x2);
- return 0.5*diff_y*diff_y;
- }
-
- virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const
- {
- return y-y_d_cont(x1,x2);
- }
-
- virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const
- {
- return 1.;
- }
-
- virtual bool fint_cont_dydy_alwayszero() const
- {
- return false;
- }
-
- virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool fint_cont_dudu_alwayszero() const
- {
- return true;
- }
-
- virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool fint_cont_dydu_alwayszero() const
- {
- return true;
- }
-
- virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
- {
- return -exp(y) - u;
- }
-
- virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
- {
- return -exp(y);
- }
-
- virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
- {
- return -1.;
- }
-
- virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
- {
- return -exp(y);
- }
-
- virtual bool d_cont_dydy_alwayszero() const
- {
- return false;
- }
-
- virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool d_cont_dudu_alwayszero() const
- {
- return true;
- }
-
- virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool d_cont_dydu_alwayszero() const
- {
- return true;
- }
- private:
-
-
- MittelmannDistCntrlNeumA2(const MittelmannDistCntrlNeumA2&);
- MittelmannDistCntrlNeumA2& operator=(const MittelmannDistCntrlNeumA2&);
-
-
- const Number pi_;
- };
- class MittelmannDistCntrlNeumA3 : public MittelmannDistCntrlNeumABase
- {
- public:
- MittelmannDistCntrlNeumA3()
- :
- pi_(4.*atan(1.)),
- M_(1.),
- K_(0.8),
- b_(1.)
- {}
- virtual ~MittelmannDistCntrlNeumA3()
- {}
- virtual bool InitializeProblem(Index N)
- {
- if (N<1) {
- printf("N has to be at least 1.");
- return false;
- }
- Number lb_y = 3.;
- Number ub_y = 6.09;
- Number lb_u = 1.4;
- Number ub_u = 1.6;
- Number b_0j = 1.;
- Number b_1j = 0.;
- Number b_i0 = 1.;
- Number b_i1 = 0.;
- Number u_init = (ub_u+lb_u)/2.;
- SetBaseParameters(N, lb_y, ub_y, lb_u, ub_u, b_0j, b_1j, b_i0, b_i1, u_init);
- return true;
- }
- protected:
-
- virtual Number y_d_cont(Number x1, Number x2) const
- {
- return 6.;
- }
-
- virtual Number fint_cont(Number x1, Number x2, Number y, Number u) const
- {
- return u*(M_*u - K_*y);
- }
-
- virtual Number fint_cont_dy(Number x1, Number x2, Number y, Number u) const
- {
- return -K_*u;
- }
-
- virtual Number fint_cont_du(Number x1, Number x2, Number y, Number u) const
- {
- return 2.*M_*u - K_*y;
- }
-
- virtual Number fint_cont_dydy(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool fint_cont_dydy_alwayszero() const
- {
- return true;
- }
-
- virtual Number fint_cont_dudu(Number x1, Number x2, Number y, Number u) const
- {
- return 2.*M_;
- }
-
- virtual bool fint_cont_dudu_alwayszero() const
- {
- return false;
- }
-
- virtual Number fint_cont_dydu(Number x1, Number x2, Number y, Number u) const
- {
- return -K_;
- }
-
- virtual bool fint_cont_dydu_alwayszero() const
- {
- return false;
- }
-
- virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
- {
- return y*(u + b_*y - a(x1,x2));
- }
-
- virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
- {
- return (u + 2.*b_*y -a(x1,x2));
- }
-
- virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
- {
- return y;
- }
-
- virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
- {
- return 2.*b_;
- }
-
- virtual bool d_cont_dydy_alwayszero() const
- {
- return false;
- }
-
- virtual Number d_cont_dudu(Number x1, Number x2, Number y, Number u) const
- {
- return 0.;
- }
-
- virtual bool d_cont_dudu_alwayszero() const
- {
- return true;
- }
-
- virtual Number d_cont_dydu(Number x1, Number x2, Number y, Number u) const
- {
- return 1.;
- }
-
- virtual bool d_cont_dydu_alwayszero() const
- {
- return false;
- }
- private:
-
-
- MittelmannDistCntrlNeumA3(const MittelmannDistCntrlNeumA3&);
- MittelmannDistCntrlNeumA3& operator=(const MittelmannDistCntrlNeumA3&);
-
-
- const Number pi_;
-
-
- const Number M_;
- const Number K_;
- const Number b_;
-
-
- inline Number a(Number x1, Number x2) const
- {
- return 7. + 4.*sin(2.*pi_*x1*x2);
- }
- };
- #endif
|