Add IpOpt wrapper for linear problems with bound- and linear constraints

dune/solvers/solvers/linearipopt.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_LINEAR_IPOPT_SOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_LINEAR_IPOPT_SOLVER_HH
+
+/** \file
+    \brief A wrapper for the IPOpt solver for linear optimization problem
+    with box constraints and linear constraints
+*/
+#include <algorithm>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+
+#include "IpTNLP.hpp"
+#include "IpIpoptApplication.hpp"
+
+/** \brief Implementation class used to pipe linear problems to
+    the IPOpt interior-point solver
+
+    \tparam JacobianType The type of the jacobian of the linear constraints
+*/
+template <class VectorType, class JacobianType>
+class LinearIPOptProblem : public Ipopt::TNLP
+{
+    enum {blocksize = VectorType::block_type::dimension};
+    enum {constraintBlocksize = JacobianType::block_type::rows};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+    /** Setup problem if no linear constraints are present */
+    LinearIPOptProblem( VectorType* x, const VectorType* rhs,
+                          const Dune::BitSetVector<blocksize>* dirichletNodes,
+                          std::vector<BoxConstraint<field_type, blocksize> >* obstacles)
+        : x_(x), rhs_(rhs),
+          dirichletNodes_(dirichletNodes), obstacles_(obstacles),
+          constraintMatrix_(nullptr), constraintObstacles_(nullptr), jacobianCutoff_(1e-8)
+    {}
+
+    /** Setup problem full problem */
+    LinearIPOptProblem(VectorType* x, const VectorType* rhs,
+                          const Dune::BitSetVector<blocksize>* dirichletNodes,
+                          std::vector<BoxConstraint<field_type, blocksize> >* boundObstacles,
+                          std::shared_ptr<const JacobianType> constraintMatrix,
+                          std::shared_ptr<const std::vector<BoxConstraint<field_type,constraintBlocksize> > > constraintObstacles)
+        : x_(x), rhs_(rhs),
+          dirichletNodes_(dirichletNodes), obstacles_(boundObstacles),
+          constraintMatrix_(constraintMatrix), constraintObstacles_(constraintObstacles),
+          jacobianCutoff_(1e-8)
+    {}
+
+  /** default destructor */
+    virtual ~LinearIPOptProblem() {}
+
+  /**@name Overloaded from TNLP */
+  //@{
+  /** Method to return some info about the nlp */
+  virtual bool get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+                            Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style);
+
+  /** Method to return the bounds for my problem */
+  virtual bool get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+                               Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u);
+
+  /** Method to return the starting point for the algorithm */
+  virtual bool get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+                                  bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+                                  Ipopt::Index m, bool init_lambda,
+                                  Ipopt::Number* lambda);
+
+  /** Method to return the objective value */
+  virtual bool eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value);
+
+  /** Method to return the gradient of the objective */
+  virtual bool eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f);
+
+  /** Method to return the constraint residuals */
+  virtual bool eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g);
+
+  /** Method to return:
+   *   1) The structure of the jacobian (if "values" is NULL)
+   *   2) The values of the jacobian (if "values" is not NULL)
+   */
+  virtual bool eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+                          Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+                          Ipopt::Number* values);
+
+  /** Method to return:
+   *   1) The structure of the hessian of the lagrangian (if "values" is NULL)
+   *   2) The values of the hessian of the lagrangian (if "values" is not NULL)
+   */
+  virtual bool eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+                      Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+                      bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+                      Ipopt::Index* jCol, Ipopt::Number* values);
+
+  //@}
+
+  /** @name Solution Methods */
+  //@{
+  /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
+    virtual void finalize_solution(Ipopt::SolverReturn status,
+                                 Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+                                 Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+                                 Ipopt::Number obj_value,
+                                   const Ipopt::IpoptData* ip_data,
+                                   Ipopt::IpoptCalculatedQuantities* ip_cq);
+  //@}
+
+    // /////////////////////////////////
+    //   Data
+    // /////////////////////////////////
+
+    //! Vector to store the solution
+    VectorType* x_;
+
+    //! The linear term of the linear energy
+    const VectorType* rhs_;
+
+    const Dune::BitSetVector<blocksize>* dirichletNodes_;
+
+    //! Vector containing the bound constraints of the variables
+    // Can stay unset when no bound constraints exist
+    std::vector<BoxConstraint<field_type, blocksize> >* obstacles_;
+
+    //! The matrix corresponding to linear constraints
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const JacobianType> constraintMatrix_;
+
+    //! IpOpt expects constraints to be of the type g_l <= g(x) <= g_u
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const std::vector<BoxConstraint<field_type, constraintBlocksize> > > constraintObstacles_;
+
+private:
+    const field_type jacobianCutoff_;
+    /**@name Methods to block default compiler methods.*/
+    //@{
+    LinearIPOptProblem(const LinearIPOptProblem&);
+    LinearIPOptProblem& operator=(const LinearIPOptProblem&);
+    //@}
+};
+
+// returns the size of the problem
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType, JacobianType>::
+get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+             Ipopt::Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style)
+{
+    // use the C style indexing (0-based)
+    index_style = Ipopt::TNLP::C_STYLE;
+
+    // Number of variables
+    n = x_->dim();
+
+    // Number of constraints
+    // Dirichlet conditions are actually bound inequality constraints
+    // bound constraints are handled differently by IpOpt
+    m = 0;
+
+    // hence the constraint jacobian is empty
+    nnz_jac_g = 0;
+
+    // energy is linear, thus the hessian is zero
+    nnz_h_lag = 0;
+
+    // We only need the lower left corner of the Hessian (since it is symmetric)
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    assert(constraintMatrix_->N()==constraintObstacles_->size());
+
+    // only change the constraint information
+    m = constraintMatrix_->N()*constraintBlocksize;
+
+    int numJacobianNonzeros = 0;
+    for (size_t i=0; i<constraintMatrix_->N(); i++) {
+
+        const typename JacobianType::row_type& row = (*constraintMatrix_)[i];
+
+        // dim for  each entry in the constraint matrix
+        typename JacobianType::row_type::ConstIterator cIt    = row.begin();
+        typename JacobianType::row_type::ConstIterator cEndIt = row.end();
+        for (; cIt!=cEndIt; cIt++)
+            for (int k=0; k<constraintBlocksize; k++)
+                for (int l=0; l<blocksize; l++)
+                    if (std::abs((*cIt)[k][l]) > jacobianCutoff_ )
+                        numJacobianNonzeros++;
+    }
+
+    nnz_jac_g = numJacobianNonzeros;
+
+    return true;
+}
+
+
+// returns the variable bounds
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType, JacobianType>::
+get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+                Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u)
+{
+    // here, the n and m we gave IPOPT in get_nlp_info are passed back to us.
+    // If desired, we could assert to make sure they are what we think they are.
+    assert(n == (Ipopt::Index)x_->dim());
+
+    if (obstacles_) {
+
+        for (size_t i=0; i<x_->size(); i++) {
+
+            for (size_t j=0; j<blocksize; j++) {
+
+                // Simply copy the values.  If they exceed a certain limit they are
+                // considered 'off'.
+                x_l[i*blocksize+j] = (*obstacles_)[i].lower(j);
+                x_u[i*blocksize+j] = (*obstacles_)[i].upper(j);
+
+            }
+
+        }
+
+    } else {
+
+        // unset all bounds
+        for (size_t i=0; i<x_->dim(); i++) {
+            x_l[i] = -std::numeric_limits<double>::max();
+            x_u[i] =  std::numeric_limits<double>::max();
+        }
+
+    }
+
+    if (dirichletNodes_) {
+
+        for (size_t i=0; i<x_->size(); i++) {
+
+            for (size_t j=0; j<blocksize; j++) {
+
+                if ( (*dirichletNodes_)[i][j])
+                    x_l[i*blocksize+j] = x_u[i*blocksize+j] = (*this->x_)[i][j];
+
+            }
+
+        }
+
+    }
+
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    // initialize with large numbers
+    for (int i=0; i<m; i++) {
+        g_l[i] = -std::numeric_limits<field_type>::max();
+        g_u[i] =  std::numeric_limits<field_type>::max();
+    }
+
+    for (size_t i=0; i<constraintObstacles_->size(); i++) {
+        for (int k=0; k<constraintBlocksize; k++) {
+
+                g_l[i*constraintBlocksize + k] = (*constraintObstacles_)[i].lower(k);
+                g_u[i*constraintBlocksize + k] = (*constraintObstacles_)[i].upper(k);
+        }
+    }
+
+    return true;
+}
+
+// returns the initial point for the problem
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+                   bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+                   Ipopt::Index m, bool init_lambda, Ipopt::Number* lambda)
+{
+    // Here, we assume we only have starting values for x, if you code
+    // your own NLP, you can provide starting values for the dual variables
+    // if you wish
+    assert(init_x == true);
+    assert(init_z == false);
+    assert(init_lambda == false);
+
+    // initialize to the given starting point
+    for (size_t i=0; i<x_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            x[i*blocksize+j] = (*x_)[i][j];
+
+    return true;
+}
+
+// returns the value of the objective function
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value)
+{
+    // Init return value
+    obj_value = 0;
+
+    ////////////////////////////////////
+    // Compute rhs*x
+    ////////////////////////////////////
+
+    for (size_t i=0; i<rhs_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            obj_value += x[i*blocksize+j] * (*rhs_)[i][j];
+
+  return true;
+}
+
+// return the gradient of the objective function grad_{x} f(x)
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f)
+{
+    // g = g - b
+    for (size_t i=0; i<rhs_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            grad_f[i*blocksize+j] = (*rhs_)[i][j];
+
+  return true;
+}
+
+// return the value of the constraints: g(x)
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
+{
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    for (int i=0; i<m; i++)
+        g[i] = 0;
+
+    int rowCounter = 0;
+    for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+        const auto& row = (*constraintMatrix_)[rowIdx];
+        auto cIt = row.begin();
+        auto cEnd = row.end();
+
+        // constraintMatrix times x
+        for (; cIt != cEnd; cIt++)
+            for (int l=0; l<constraintBlocksize; l++)
+                for (int k=0; k<blocksize; k++)
+                    if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+                        g[rowCounter+l] += (*cIt)[l][k] * x[cIt.index()*blocksize+k];
+
+        // increment constraint counter
+        rowCounter += constraintBlocksize;
+    }
+
+    return true;
+
+}
+
+// return the structure or values of the jacobian
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+           Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+           Ipopt::Number* values)
+{
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    int idx = 0;
+    int rowCounter = 0;
+
+    // If the values are null then IpOpt wants the sparsity pattern
+    if (values==NULL) {
+
+        for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+            const auto& row = (*constraintMatrix_)[rowIdx];
+            auto cIt = row.begin();
+            auto cEnd = row.end();
+
+            for (; cIt != cEnd; ++cIt) {
+
+                for (int l=0; l<constraintBlocksize; l++)
+                    for (int k=0; k<blocksize; k++)
+                        if ( std::abs((*cIt)[l][k]) > jacobianCutoff_) {
+
+                            iRow[idx] = rowCounter + l;
+                            jCol[idx] = cIt.index()*blocksize+k;
+                            idx++;
+                        }
+            }
+
+            rowCounter += constraintBlocksize;
+        }
+
+    // If the values are not null IpOpt wants the evaluated constraints gradient
+    } else {
+
+        for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+            const auto& row = (*constraintMatrix_)[rowIdx];
+            auto cIt = row.begin();
+            auto cEnd = row.end();
+
+            for (; cIt != cEnd; ++cIt) {
+
+                for (int l=0; l<constraintBlocksize; l++)
+                    for (int k=0; k<blocksize; k++)
+                        if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+                            values[idx++] = (*cIt)[l][k];
+            }
+
+            rowCounter += constraintBlocksize;
+
+        }
+
+    }
+
+    return true;
+
+}
+
+//return the structure or values of the hessian
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+       Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+       bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+       Ipopt::Index* jCol, Ipopt::Number* values)
+{
+    return true;
+}
+
+template <class VectorType, class JacobianType>
+void LinearIPOptProblem<VectorType,JacobianType>::
+finalize_solution(Ipopt::SolverReturn status,
+                  Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+                  Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+                  Ipopt::Number obj_value,
+                  const Ipopt::IpoptData* ip_data,
+                  Ipopt::IpoptCalculatedQuantities* ip_cq)
+{
+  // here is where we would store the solution to variables, or write to a file, etc
+  // so we could use the solution.
+
+    for (size_t i=0; i<x_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            (*x_)[i][j] = x[i*blocksize+j];
+
+}
+
+/** \brief Wraps the IPOpt interior-point solver for quadratic problems with
+ *    linear constraints and bound constraints.
+ *
+ *   The problems that can be solved are of the form
+ *   \f[
+ *       \min_x \frac{1}{2}x^T A x - b^T x
+ *               x_l \leq x \leq x_u
+ *               g_l \leq G x \leq g_u
+ *   \f]
+ *
+ *  \tparam JacobianType The type of the jacobian of the linear constraints
+*/
+template <class VectorType, class JacobianType>
+class LinearIPOptSolver
+    : public IterativeSolver<VectorType>, public CanIgnore<Dune::BitSetVector<VectorType::block_type::dimension> >
+{
+
+    enum {blocksize = VectorType::block_type::dimension};
+    enum {constraintBlocksize = JacobianType::block_type::rows};
+
+    typedef typename VectorType::field_type field_type;
+
+    typedef Dune::FieldMatrix<field_type, blocksize, blocksize> MatrixBlock;
+    typedef Dune::FieldVector<field_type, blocksize>            VectorBlock;
+
+
+public:
+
+    /** \brief Default constructor */
+    LinearIPOptSolver () : IterativeSolver<VectorType>(1e-8, 100, NumProc::FULL),
+                              rhs_(NULL),
+                              obstacles_(NULL),
+                              linearSolverType_(""),
+                              constraintObstacles_(nullptr),
+                              constraintMatrix_(nullptr)
+    {}
+
+    /** \brief Constructor for a linear problem */
+    LinearIPOptSolver (VectorType& x,
+                       const VectorType& rhs,
+                       NumProc::VerbosityMode verbosity=NumProc::FULL,
+                       std::string linearSolverType = "")
+        : IterativeSolver<VectorType>(1e-8, 100, verbosity),
+          rhs_(&rhs), obstacles_(NULL),
+          linearSolverType_(linearSolverType),
+          constraintObstacles_(nullptr),
+          constraintMatrix_(nullptr)
+    {
+        this->x_ = &x;
+    }
+
+    void setProblem(VectorType& x,
+                    const VectorType& rhs) {
+        x_ = &x;
+        rhs_ = &rhs;
+    }
+
+    void setLinearConstraints(const JacobianType& constraintMatrix,
+                     const std::vector<BoxConstraint<field_type,constraintBlocksize> >& obstacles) {
+        constraintMatrix_ = Dune::stackobject_to_shared_ptr(constraintMatrix);
+        constraintObstacles_ = Dune::stackobject_to_shared_ptr(obstacles);
+    }
+
+    virtual void solve();
+
+    ///////////////////////////////////////////////////////
+    // Data
+    ///////////////////////////////////////////////////////
+
+    //! Vector to store the solution
+    VectorType* x_;
+
+    //! The linear term in the linear energy
+    const VectorType* rhs_;
+
+    //! Vector containing the bound constraints of the variables
+    // Can stay unset when no bound constraints exist
+    std::vector<BoxConstraint<field_type,blocksize> >* obstacles_;
+
+    //! The type of the linear solver to be used within IpOpt, default is ""
+    std::string linearSolverType_;
+
+    //! IpOpt expects constraints to be of the type g_l <= g(x) <= g_u
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const std::vector<BoxConstraint<field_type,constraintBlocksize> > > constraintObstacles_;
+
+    //! The matrix corresponding to linear constraints
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const JacobianType> constraintMatrix_;
+};
+
+template <class VectorType, class JacobianType>
+void LinearIPOptSolver<VectorType,JacobianType>::solve()
+{
+  // Create a new instance of your nlp
+  //  (use a SmartPtr, not raw)
+    Ipopt::SmartPtr<Ipopt::TNLP> mynlp = new LinearIPOptProblem<VectorType,JacobianType>(x_, rhs_,
+                                                           this->ignoreNodes_, obstacles_,
+                                                            constraintMatrix_,constraintObstacles_);
+
+  // Create a new instance of IpoptApplication
+  //  (use a SmartPtr, not raw)
+    Ipopt::SmartPtr<Ipopt::IpoptApplication> app = new Ipopt::IpoptApplication();
+
+  // Change some options
+  app->Options()->SetNumericValue("tol", this->tolerance_);
+  app->Options()->SetIntegerValue("max_iter", this->maxIterations_);
+  app->Options()->SetStringValue("mu_strategy", "adaptive");
+  if (linearSolverType_ != "")
+    app->Options()->SetStringValue("linear_solver",linearSolverType_);
+  app->Options()->SetStringValue("output_file", "ipopt.out");
+  app->Options()->SetStringValue("hessian_constant", "yes");
+
+  // We only support linear constraints
+  app->Options()->SetStringValue("jac_d_constant","yes");
+  app->Options()->SetStringValue("jac_c_constant","yes");
+
+  switch (this->verbosity_) {
+  case NumProc::QUIET:
+      app->Options()->SetIntegerValue("print_level", 0);
+      break;
+  case NumProc::REDUCED:
+      app->Options()->SetIntegerValue("print_level", 5);
+      break;
+  case NumProc::FULL:
+      app->Options()->SetIntegerValue("print_level", 12);
+  }
+
+  // Intialize the IpoptApplication and process the options
+  Ipopt::ApplicationReturnStatus status;
+  status = app->Initialize();
+  if (status != Ipopt::Solve_Succeeded)
+      DUNE_THROW(Dune::Exception, "IPOpt: Error during initialization!");
+
+  // Ask Ipopt to solve the problem
+  status = app->OptimizeTNLP(mynlp);
+
+  // Print helpful text for at least some cases
+  if (status == Ipopt::Invalid_Option)
+    DUNE_THROW(Dune::Exception, "IPOpt returned 'Invalid_Option' error!");
+
+  if (status == Ipopt::Solved_To_Acceptable_Level)
+      std::cout<<"WARNING: Desired tolerance could not be reached, but still accetable tolerance is reached.\n";
+  else if (status != Ipopt::Solve_Succeeded)
+      DUNE_THROW(Dune::Exception, "IPOpt: Error during optimization!");
+}
+
+#endif