MittelmannDistCntrlNeumB.hpp

Go to the documentation of this file.

// Copyright (C) 2005, 2006 International Business Machines and others.
// All Rights Reserved.
// This code is published under the Eclipse Public License.
//
// $Id: MittelmannDistCntrlNeumB.hpp 2005 2011-06-06 12:55:16Z stefan $
//
// Authors:  Andreas Waechter           IBM    2005-10-18
//                  based on MyNLP.hpp
 
#ifndef __MITTELMANNDISTRCNTRLNEUMB_HPP__
#define __MITTELMANNDISTRCNTRLNEUMB_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 MittelmannDistCntrlNeumBBase : public RegisteredTNLP
{
public:
  MittelmannDistCntrlNeumBBase();
 
  virtual ~MittelmannDistCntrlNeumBBase();
 
  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:
  MittelmannDistCntrlNeumBBase(const MittelmannDistCntrlNeumBBase&);
  MittelmannDistCntrlNeumBBase& operator=(const MittelmannDistCntrlNeumBBase&);
 
  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 MittelmannDistCntrlNeumB1 : public MittelmannDistCntrlNeumBBase
{
public:
  MittelmannDistCntrlNeumB1()
      :
      pi_(4.*atan(1.)),
      alpha_(0.001)
  {}
 
  virtual ~MittelmannDistCntrlNeumB1()
  {}
 
  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:
  MittelmannDistCntrlNeumB1(const MittelmannDistCntrlNeumB1&);
  MittelmannDistCntrlNeumB1& operator=(const MittelmannDistCntrlNeumB1&);
 
  const Number pi_;
  const Number alpha_;
};
 
class MittelmannDistCntrlNeumB2 : public MittelmannDistCntrlNeumBBase
{
public:
  MittelmannDistCntrlNeumB2()
      :
      pi_(4.*atan(1.))
  {}
 
  virtual ~MittelmannDistCntrlNeumB2()
  {}
 
  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:
  MittelmannDistCntrlNeumB2(const MittelmannDistCntrlNeumB2&);
  MittelmannDistCntrlNeumB2& operator=(const MittelmannDistCntrlNeumB2&);
 
  const Number pi_;
};
 
class MittelmannDistCntrlNeumB3 : public MittelmannDistCntrlNeumBBase
{
public:
  MittelmannDistCntrlNeumB3()
      :
      pi_(4.*atan(1.)),
      M_(1.),
      K_(0.8),
      b_(1.)
  {}
 
  virtual ~MittelmannDistCntrlNeumB3()
  {}
 
  virtual bool InitializeProblem(Index N)
  {
    if (N<1) {
      printf("N has to be at least 1.");
      return false;
    }
    Number lb_y = 3.;//-1e20;
    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:
  MittelmannDistCntrlNeumB3(const MittelmannDistCntrlNeumB3&);
  MittelmannDistCntrlNeumB3& operator=(const MittelmannDistCntrlNeumB3&);
 
  const Number pi_;
  /*@name constrants appearing in problem formulation */
  const Number M_;
  const Number K_;
  const Number b_;
  //* Auxiliary function for state equation */
  inline Number a(Number x1, Number x2) const
  {
    return 7. + 4.*sin(2.*pi_*x1*x2);
  }
};
 
#endif