MittelmannBndryCntrlDiri3D.hpp

Go to the documentation of this file.

// Copyright (C) 2005, 2007 International Business Machines and others.
// All Rights Reserved.
// This code is published under the Eclipse Public License.
//
// $Id: MittelmannBndryCntrlDiri3D.hpp 2005 2011-06-06 12:55:16Z stefan $
//
// Authors:  Andreas Waechter             IBM    2005-10-18
//           Olaf Schenk   (Univ. of Basel)      2007-08-01
//              modified MittelmannBndryCntrlDiri.hpp for 3-dim problem

#ifndef __MITTELMANNBNDRYCNTRLDIRI3D_HPP__
#define __MITTELMANNBNDRYCNTRLDIRI3D_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 MittelmannBndryCntrlDiriBase3D : public RegisteredTNLP
{
public:
  MittelmannBndryCntrlDiriBase3D();

  virtual ~MittelmannBndryCntrlDiriBase3D();

  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_valu,
                                 const IpoptData* ip_data,
                                 IpoptCalculatedQuantities* ip_cq);

protected:
  void SetBaseParameters(Index N, Number alpha, Number lb_y,
                         Number ub_y, Number lb_u, Number ub_u,
                         Number d_const, Number B, Number C);

  virtual Number y_d_cont(Number x1, Number x2, Number x3) const =0;

private:
  MittelmannBndryCntrlDiriBase3D(const MittelmannBndryCntrlDiriBase3D&);
  MittelmannBndryCntrlDiriBase3D& operator=(const MittelmannBndryCntrlDiriBase3D&);

  Index N_;
  Number h_;
  Number hh_;
  Number hhh_;
  Number lb_y_;
  Number ub_y_;
  Number lb_u_;
  Number ub_u_;
  Number d_const_;
  Number alpha_;
  Number* y_d_;

  inline Index y_index(Index i, Index j, Index k) const
  {
    return k + (N_+2)*j + (N_+2)*(N_+2)*i;
  }
  inline Index pde_index(Index i, Index j, Index k) const
  {
    return (k-1) + N_*(j-1) + N_*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;
  }
  inline Number x3_grid(Index i) const
  {
    return h_*(Number)i;
  }
  inline Number PenObj(Number t) const
  {
    //return 0.5*t*t;
    if (t > B_) {
      return B_*B_/2. + C_*(t - B_);
    }
    else if (t < -B_) {
      return B_*B_/2. + C_*(-t - B_);
    }
    else {
      const Number t2 = t*t;
      const Number t4 = t2*t2;
      const Number t6 = t4*t2;
      return PenA_*t2 + PenB_*t4 + PenC_*t6;
    }
  }
  inline Number PenObj_1(Number t) const
  {
    //return t;
    if (t > B_) {
      return C_;
    }
    else if (t < -B_) {
      return -C_;
    }
    else {
      const Number t2 = t*t;
      const Number t3 = t*t2;
      const Number t5 = t3*t2;
      return 2.*PenA_*t + 4.*PenB_*t3 + 6.*PenC_*t5;
    }
  }
  inline Number PenObj_2(Number t) const
  {
    //return 1.;
    if (t > B_) {
      return 0.;
    }
    else if (t < -B_) {
      return 0.;
    }
    else {
      const Number t2 = t*t;
      const Number t4 = t2*t2;
      return 2.*PenA_ + 12.*PenB_*t2 + 30.*PenC_*t4;
    }
  }

  Number B_;
  Number C_;
  Number PenA_;
  Number PenB_;
  Number PenC_;
};

class MittelmannBndryCntrlDiri3D : public MittelmannBndryCntrlDiriBase3D
{
public:
  MittelmannBndryCntrlDiri3D()
  {}

  virtual ~MittelmannBndryCntrlDiri3D()
  {}

  virtual bool InitializeProblem(Index N)
  {
    if (N<1) {
      printf("N has to be at least 1.");
      return false;
    }
    Number alpha = 0.01;
    Number lb_y = -1e20;
    Number ub_y = 3.5;
    Number lb_u = 0.;
    Number ub_u = 10.;
    Number d_const = -20.;
    Number B = .5;
    Number C = 0.01;
    SetBaseParameters(N, alpha, lb_y, ub_y, lb_u, ub_u, d_const, B, C);
    return true;
  }
protected:
  virtual Number y_d_cont(Number x1, Number x2, Number x3)  const
  {
    return 3. + 5.*(x1*(x1-1.)*x2*(x2-1.));
  }
private:
  MittelmannBndryCntrlDiri3D(const MittelmannBndryCntrlDiri3D&);
  MittelmannBndryCntrlDiri3D& operator=(const MittelmannBndryCntrlDiri3D&);

};


#endif