00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #ifndef __MITTELMANNDISTRCNTRLDIRI_HPP__
00011 #define __MITTELMANNDISTRCNTRLDIRI_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 using namespace Ipopt;
00027
00034 class MittelmannDistCntrlDiriBase : public RegisteredTNLP
00035 {
00036 public:
00039 MittelmannDistCntrlDiriBase();
00040
00042 virtual ~MittelmannDistCntrlDiriBase();
00043
00047 virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
00048 Index& nnz_h_lag, IndexStyleEnum& index_style);
00049
00051 virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u,
00052 Index m, Number* g_l, Number* g_u);
00053
00055 virtual bool get_starting_point(Index n, bool init_x, Number* x,
00056 bool init_z, Number* z_L, Number* z_U,
00057 Index m, bool init_lambda,
00058 Number* lambda);
00059
00061 virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value);
00062
00064 virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f);
00065
00067 virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g);
00068
00073 virtual bool eval_jac_g(Index n, const Number* x, bool new_x,
00074 Index m, Index nele_jac, Index* iRow, Index *jCol,
00075 Number* values);
00076
00081 virtual bool eval_h(Index n, const Number* x, bool new_x,
00082 Number obj_factor, Index m, const Number* lambda,
00083 bool new_lambda, Index nele_hess, Index* iRow,
00084 Index* jCol, Number* values);
00085
00087
00089 virtual bool get_scaling_parameters(Number& obj_scaling,
00090 bool& use_x_scaling, Index n,
00091 Number* x_scaling,
00092 bool& use_g_scaling, Index m,
00093 Number* g_scaling);
00094
00099 virtual void finalize_solution(SolverReturn status,
00100 Index n, const Number* x, const Number* z_L, const Number* z_U,
00101 Index m, const Number* g, const Number* lambda,
00102 Number obj_value);
00104
00105 protected:
00109 void SetBaseParameters(Index N, Number alpha, Number lb_y,
00110 Number ub_y, Number lb_u, Number ub_u,
00111 Number u_init);
00112
00116 virtual Number y_d_cont(Number x1, Number x2) const =0;
00118 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const =0;
00120 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const =0;
00122 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const =0;
00124 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const =0;
00126
00127 private:
00139 MittelmannDistCntrlDiriBase(const MittelmannDistCntrlDiriBase&);
00140 MittelmannDistCntrlDiriBase& operator=(const MittelmannDistCntrlDiriBase&);
00142
00146 Index N_;
00148 Number h_;
00150 Number hh_;
00152 Number lb_y_;
00154 Number ub_y_;
00156 Number lb_u_;
00158 Number ub_u_;
00160 Number u_init_;
00163 Number alpha_;
00165 Number* y_d_;
00167
00172 inline Index y_index(Index i, Index j) const
00173 {
00174 return j + (N_+2)*i;
00175 }
00178 inline Index u_index(Index i, Index j) const
00179 {
00180 return (N_+2)*(N_+2) + (j-1) + (N_)*(i-1);
00181 }
00184 inline Index pde_index(Index i, Index j) const
00185 {
00186 return (j-1) + N_*(i-1);
00187 }
00189 inline Number x1_grid(Index i) const
00190 {
00191 return h_*(Number)i;
00192 }
00194 inline Number x2_grid(Index i) const
00195 {
00196 return h_*(Number)i;
00197 }
00199 };
00200
00202 class MittelmannDistCntrlDiri1 : public MittelmannDistCntrlDiriBase
00203 {
00204 public:
00205 MittelmannDistCntrlDiri1()
00206 {}
00207
00208 virtual ~MittelmannDistCntrlDiri1()
00209 {}
00210
00211 virtual bool InitializeProblem(Index N)
00212 {
00213 if (N<1) {
00214 printf("N has to be at least 1.");
00215 return false;
00216 }
00217 Number alpha = 0.001;
00218 Number lb_y = -1e20;
00219 Number ub_y = 0.185;
00220 Number lb_u = 1.5;
00221 Number ub_u = 4.5;
00222 Number u_init = (ub_u+lb_u)/2.;
00223
00224 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00225 return true;
00226 }
00227 protected:
00229 virtual Number y_d_cont(Number x1, Number x2) const
00230 {
00231 return 1. + 2.*(x1*(x1-1.)+x2*(x2-1.));
00232 }
00234 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
00235 {
00236 return pow(y,3) - y - u;
00237 }
00239 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
00240 {
00241 return 3.*y*y - 1.;
00242 }
00244 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
00245 {
00246 return -1.;
00247 }
00249 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
00250 {
00251 return 6.*y;
00252 }
00253 private:
00256 MittelmannDistCntrlDiri1(const MittelmannDistCntrlDiri1&);
00257 MittelmannDistCntrlDiri1& operator=(const MittelmannDistCntrlDiri1&);
00259 };
00260
00262 class MittelmannDistCntrlDiri2 : public MittelmannDistCntrlDiriBase
00263 {
00264 public:
00265 MittelmannDistCntrlDiri2()
00266 {}
00267 virtual ~MittelmannDistCntrlDiri2()
00268 {}
00269 virtual bool InitializeProblem(Index N)
00270 {
00271 if (N<1) {
00272 printf("N has to be at least 1.");
00273 return false;
00274 }
00275 Number alpha = 0.;
00276 Number lb_y = -1e20;
00277 Number ub_y = 0.185;
00278 Number lb_u = 1.5;
00279 Number ub_u = 4.5;
00280 Number u_init = (ub_u+lb_u)/2.;
00281
00282 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00283 return true;
00284 }
00285 protected:
00287 virtual Number y_d_cont(Number x1, Number x2) const
00288 {
00289 return 1. + 2.*(x1*(x1-1.)+x2*(x2-1.));
00290 }
00292 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
00293 {
00294 return pow(y,3) - y - u;
00295 }
00297 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
00298 {
00299 return 3.*y*y - 1.;
00300 }
00302 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
00303 {
00304 return -1.;
00305 }
00307 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
00308 {
00309 return 6.*y;
00310 }
00311 private:
00314 MittelmannDistCntrlDiri2(const MittelmannDistCntrlDiri2&);
00315 MittelmannDistCntrlDiri2& operator=(const MittelmannDistCntrlDiri2&);
00317 };
00318
00320 class MittelmannDistCntrlDiri3 : public MittelmannDistCntrlDiriBase
00321 {
00322 public:
00323 MittelmannDistCntrlDiri3()
00324 :
00325 pi_(4.*atan(1.))
00326 {}
00327 virtual ~MittelmannDistCntrlDiri3()
00328 {}
00329 virtual bool InitializeProblem(Index N)
00330 {
00331 if (N<1) {
00332 printf("N has to be at least 1.");
00333 return false;
00334 }
00335 Number alpha = 0.001;
00336 Number lb_y = -1e20;
00337 Number ub_y = 0.11;
00338 Number lb_u = -5;
00339 Number ub_u = 5.;
00340 Number u_init = (ub_u+lb_u)/2.;
00341
00342 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00343 return true;
00344 }
00345 protected:
00347 virtual Number y_d_cont(Number x1, Number x2) const
00348 {
00349 return sin(2.*pi_*x1)*sin(2.*pi_*x2);
00350 }
00352 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
00353 {
00354 return -exp(y) - u;
00355 }
00357 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
00358 {
00359 return -exp(y);
00360 }
00362 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
00363 {
00364 return -1.;
00365 }
00367 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
00368 {
00369 return -exp(y);
00370 }
00371 private:
00374 MittelmannDistCntrlDiri3(const MittelmannDistCntrlDiri3&);
00375 MittelmannDistCntrlDiri3& operator=(const MittelmannDistCntrlDiri3&);
00377
00378 const Number pi_;
00379 };
00380
00381 class MittelmannDistCntrlDiri3a : public MittelmannDistCntrlDiriBase
00382 {
00383 public:
00384 MittelmannDistCntrlDiri3a()
00385 :
00386 pi_(4.*atan(1.))
00387 {}
00388 virtual ~MittelmannDistCntrlDiri3a()
00389 {}
00390 virtual bool InitializeProblem(Index N)
00391 {
00392 if (N<1) {
00393 printf("N has to be at least 1.");
00394 return false;
00395 }
00396 Number alpha = 0.;
00397 Number lb_y = -1e20;
00398 Number ub_y = 0.11;
00399 Number lb_u = -5;
00400 Number ub_u = 5.;
00401 Number u_init = (ub_u+lb_u)/2.;
00402
00403 SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, u_init);
00404 return true;
00405 }
00406 protected:
00408 virtual Number y_d_cont(Number x1, Number x2) const
00409 {
00410 return sin(2.*pi_*x1)*sin(2.*pi_*x2);
00411 }
00413 virtual Number d_cont(Number x1, Number x2, Number y, Number u) const
00414 {
00415 return -exp(y) - u;
00416 }
00418 virtual Number d_cont_dy(Number x1, Number x2, Number y, Number u) const
00419 {
00420 return -exp(y);
00421 }
00423 virtual Number d_cont_du(Number x1, Number x2, Number y, Number u) const
00424 {
00425 return -1.;
00426 }
00428 virtual Number d_cont_dydy(Number x1, Number x2, Number y, Number u) const
00429 {
00430 return -exp(y);
00431 }
00432 private:
00435 MittelmannDistCntrlDiri3a(const MittelmannDistCntrlDiri3a&);
00436 MittelmannDistCntrlDiri3a& operator=(const MittelmannDistCntrlDiri3a&);
00438
00439 const Number pi_;
00440 };
00441
00442 #endif