Add new class QuadraticIPOptObstacleProblem which, additionally

to bound constraints, also allow to prescribe linear constraints.

The QuadraticIpOptSolver chooses this problem class if a
constraintMatrix and constraintObstacles have been set.

dune/solvers/solvers/quadraticipopt.hh

@@ -129,7 +129,7 @@ private:
 template <class MatrixType, class VectorType>
 bool QuadraticIPOptProblem<MatrixType,VectorType>::
 get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
-             Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style)
+             Ipopt::Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style)
 {
     // Number of variables
     n = x_->dim();
@@ -186,7 +186,7 @@ get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_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());
-    assert(m == 0);
+    //assert(m == 0);
 
     if (obstacles_) {
 
@@ -461,9 +461,234 @@ finalize_solution(Ipopt::SolverReturn status,
 
 }
 
+/** \brief Implementation class used to pipe quadratic problems to
+    the IPOpt interior-point solver
+*/
+template <class MatrixType, class VectorType>
+class QuadraticIPOptObstacleProblem : public QuadraticIPOptProblem<MatrixType, VectorType>
+{
+     enum {blocksize = VectorType::block_type::dimension};
+
+    typedef typename VectorType::field_type field_type;
+    typedef QuadraticIPOptProblem<MatrixType, VectorType> Base;
+
+public:
+  /** default constructor */
+    QuadraticIPOptObstacleProblem(const MatrixType* hessian, VectorType* x, const VectorType* rhs,
+                          const Dune::BitSetVector<blocksize>* dirichletNodes,
+                          std::vector<BoxConstraint<field_type, blocksize> >* boundObstacles,
+                          const MatrixType* constraintMatrix,
+                          std::vector<BoxConstraint<field_type, blocksize> >* constraintObstacles)
+        : Base(hessian,x,rhs,dirichletNodes,boundObstacles),
+          constraintMatrix_(constraintMatrix), constraintObstacles_(constraintObstacles),
+          jacobianCutoff_(1e-8)
+    {}
+
+  /** default destructor */
+    virtual ~QuadraticIPOptObstacleProblem() {};
+
+ /**@name Overloaded from QuadraticIPOptProblem */
+  //@{
+  /** 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, Ipopt::TNLP::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 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);
+  //@}
+
+
+  // /////////////////////////////////
+  //   Data
+  // /////////////////////////////////
+
+  const MatrixType* constraintMatrix_;
+
+  std::vector<BoxConstraint<field_type, blocksize> >* constraintObstacles_;
+
+private:
+  const field_type jacobianCutoff_;
+  /**@name Methods to block default compiler methods.
+  */
+  //@{
+  //  QuadraticIPOptObstacleProblem();
+  QuadraticIPOptObstacleProblem(const QuadraticIPOptObstacleProblem&);
+  QuadraticIPOptObstacleProblem& operator=(const QuadraticIPOptObstacleProblem&);
+  //@}
+};
+
+// returns the size of the problem
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::
+get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+             Ipopt::Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style)
+{
+    Base::get_nlp_info(n,m,nnz_jac_g,nnz_h_lag,index_style);
+
+    // We only need the lower left corner of the Hessian (since it is symmetric)
+    if (!constraintMatrix_)
+        DUNE_THROW(SolverError, "No constraint matrix has been supplied!");
+    assert(constraintMatrix_->N()==constraintObstacles_->size());
+
+    // only change the constraint information
+    m = constraintMatrix_->N()*blocksize;
+
+    int numJacobianNonzeros = 0;
+    for (size_t i=0; i<constraintMatrix_->N(); i++) {
+
+        const auto& row = (*constraintMatrix_)[i];
+
+        // dim for entry in the constraint matrix
+        auto cEnd = row.end();
+        for (auto cIt = row.begin(); cIt!=cEnd; cIt++)
+            for (int k=0; k<blocksize; 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 MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::
+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)
+{
+    Base::get_bounds_info(n,x_l, x_u, m, g_l ,g_u);
+
+    // 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<blocksize; k++) {
+
+                g_l[i*blocksize + k] = (*constraintObstacles_)[i].lower(k);
+                g_u[i*blocksize + k] = (*constraintObstacles_)[i].upper(k);
+        }
+    }
+
+    return true;
+}
+
+// Evaluate the constraints
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
+{
+
+    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<blocksize; 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 += blocksize;
+    }
+
+    return true;
+}
+
+
+// Evaluate gradient of the constraints
+template <class MatrixType, class VectorType>
+bool QuadraticIPOptObstacleProblem<MatrixType,VectorType>::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)
+{
+
+    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<blocksize; 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 += blocksize;
+        }
+
+    // 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<blocksize; l++)
+                    for (int k=0; k<blocksize; k++)
+                        if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+                            values[idx++] = (*cIt)[l][k];
+            }
+
+            rowCounter += blocksize;
+
+        }
+
+    }
+
+    return true;
+}
 
 /** \brief Wraps the IPOpt interior-point solver for quadratic problems with
-    bound constraints.
+ *    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]
 */
 template <class MatrixType, class VectorType>
 class QuadraticIPOptSolver
@@ -484,7 +709,8 @@ public:
     QuadraticIPOptSolver () : IterativeSolver<VectorType>(1e-8, 100, NumProc::FULL),
                               hessian_(NULL), rhs_(NULL),
                               obstacles_(NULL),
-                              linearSolverType_("ma27")
+                              linearSolverType_("ma27"),
+                              constraintMatrix_(NULL),constraintObstacles_(NULL)
     {}
 
     /** \brief Constructor for a linear problem */
@@ -495,7 +721,8 @@ public:
                        std::string linearSolverType = "ma27")
         : IterativeSolver<VectorType>(1e-8, 100, verbosity),
           hessian_(&hessian), rhs_(&rhs), obstacles_(NULL),
-          linearSolverType_(linearSolverType)
+          linearSolverType_(linearSolverType),
+          constraintMatrix_(NULL),constraintObstacles_(NULL)
     {
         this->x_ = &x;
     }
@@ -508,22 +735,41 @@ public:
         rhs_ = &rhs;
     }
 
+    void setLinearConstraints(const MatrixType& constraintMatrix,
+                    const std::vector<BoxConstraint<field_type,blocksize> >& obstacles) {
+        constraintMatrix_ = &constraintMatrix;
+        constraintObstacles_ = &obstacles;
+    }
+
     virtual void solve();
 
     ///////////////////////////////////////////////////////
     // Data
     ///////////////////////////////////////////////////////
 
+    //! The quadratic term in the quadratic energy, is assumed to be symmetric
     const MatrixType* hessian_;
 
+    //! Vector to store the solution
     VectorType* x_;
 
+    //! The linear term in the quadratic 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 ma27
     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::vector<BoxConstraint<field_type,blocksize> >* constraintObstacles_;
+
+    //! The matrix corresponding to linear constraints
+    // Can stay unset when no linear constraints exist
+    const MatrixType* constraintMatrix_;
 };
 
 template <class MatrixType, class VectorType>
@@ -531,8 +777,13 @@ void QuadraticIPOptSolver<MatrixType,VectorType>::solve()
 {
   // Create a new instance of your nlp
   //  (use a SmartPtr, not raw)
-    Ipopt::SmartPtr<Ipopt::TNLP> mynlp
-        = new QuadraticIPOptProblem<MatrixType,VectorType>(hessian_, x_, rhs_,
+    Ipopt::SmartPtr<Ipopt::TNLP> mynlp;
+    if (constraintMatrix_) {
+        mynlp= new QuadraticIPOptObstacleProblem<MatrixType,VectorType>(hessian_, x_, rhs_,
+                                                           this->ignoreNodes_, obstacles_,
+                                                            constraintMatrix_,constraintObstacles_);
+    } else
+        mynlp= new QuadraticIPOptProblem<MatrixType,VectorType>(hessian_, x_, rhs_,
                                                            this->ignoreNodes_, obstacles_);
 
   // Create a new instance of IpoptApplication
@@ -546,6 +797,13 @@ void QuadraticIPOptSolver<MatrixType,VectorType>::solve()
   app->Options()->SetStringValue("linear_solver",linearSolverType_);
   app->Options()->SetStringValue("output_file", "ipopt.out");
   app->Options()->SetStringValue("hessian_constant", "yes");
+
+  // We only support linear constraints
+  if (constraintMatrix_) {
+    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);