Skip to content
Snippets Groups Projects
Commit f2e194e7 authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Merge QuadraticIPOptProblem and QuadraticIPOptObstacleProblem

parent edc2f6dd
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/solvers/solvers/quadraticipopt.hh
+ 163
− 245
View file @ f2e194e7
......@@ -35,12 +35,25 @@ class QuadraticIPOptProblem : public Ipopt::TNLP
typedef typename VectorType::field_type field_type;
public:
/** default constructor */
/** Setup problem if no linear constraints are present */
QuadraticIPOptProblem(const MatrixType* hessian, VectorType* x, const VectorType* rhs,
const Dune::BitSetVector<blocksize>* dirichletNodes,
std::vector<BoxConstraint<field_type, blocksize> >* obstacles)
: hessian_(hessian), x_(x), rhs_(rhs),
dirichletNodes_(dirichletNodes), obstacles_(obstacles)
dirichletNodes_(dirichletNodes), obstacles_(obstacles),
constraintMatrix_(NULL), constraintObstacles_(NULL), jacobianCutoff_(1e-8)
{}
/** Setup problem full problem */
QuadraticIPOptProblem(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)
: hessian_(hessian), x_(x), rhs_(rhs),
dirichletNodes_(dirichletNodes), obstacles_(boundObstacles),
constraintMatrix_(constraintMatrix), constraintObstacles_(constraintObstacles),
jacobianCutoff_(1e-8)
{}
/** default destructor */
......@@ -105,24 +118,36 @@ public:
// 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_;
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
const MatrixType* constraintMatrix_;
//! 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_;
private:
/**@name Methods to block default compiler methods.
*/
//@{
// QuadraticIPOptProblem();
QuadraticIPOptProblem(const QuadraticIPOptProblem&);
QuadraticIPOptProblem& operator=(const QuadraticIPOptProblem&);
//@}
const field_type jacobianCutoff_;
/**@name Methods to block default compiler methods.*/
//@{
QuadraticIPOptProblem(const QuadraticIPOptProblem&);
QuadraticIPOptProblem& operator=(const QuadraticIPOptProblem&);
//@}
};
// returns the size of the problem
......@@ -134,7 +159,8 @@ get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
// Number of variables
n = x_->dim();
// No real constraints: Dirichlet conditions are actually bound inequality constraints
// Number of constraints
// Dirichlet conditions are actually bound inequality constraints
// bound constraints are handled differently by IpOpt
m = 0;
......@@ -170,8 +196,31 @@ get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
}
// use the C style indexing (0-based)
index_style = Ipopt::TNLP::C_STYLE;
// We only need the lower left corner of the Hessian (since it is symmetric)
if (!constraintMatrix_)
return true;
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 typename MatrixType::row_type& row = (*constraintMatrix_)[i];
// dim for each entry in the constraint matrix
typename MatrixType::row_type::ConstIterator cIt = row.begin();
typename MatrixType::row_type::ConstIterator cEndIt = row.end();
for (; cIt!=cEndIt; 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;
}
......@@ -186,7 +235,6 @@ 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);
if (obstacles_) {
......@@ -228,7 +276,24 @@ get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
}
return true;
if (!constraintMatrix_)
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<blocksize; k++) {
g_l[i*blocksize + k] = (*constraintObstacles_)[i].lower(k);
g_u[i*blocksize + k] = (*constraintObstacles_)[i].upper(k);
}
}
return true;
}
// returns the initial point for the problem
......@@ -346,7 +411,32 @@ template <class MatrixType, class VectorType>
bool QuadraticIPOptProblem<MatrixType,VectorType>::
eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
{
return true;
if (!constraintMatrix_)
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<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;
}
// return the structure or values of the jacobian
......@@ -356,7 +446,61 @@ 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)
{
return true;
if (!constraintMatrix_)
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<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;
}
//return the structure or values of the hessian
......@@ -461,225 +605,6 @@ 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
* linear constraints and bound constraints.
*
......@@ -777,14 +702,9 @@ void QuadraticIPOptSolver<MatrixType,VectorType>::solve()
{
// Create a new instance of your nlp
// (use a SmartPtr, not raw)
Ipopt::SmartPtr<Ipopt::TNLP> mynlp;
if (constraintMatrix_) {
mynlp= new QuadraticIPOptObstacleProblem<MatrixType,VectorType>(hessian_, x_, rhs_,
Ipopt::SmartPtr<Ipopt::TNLP> mynlp = new QuadraticIPOptProblem<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
// (use a SmartPtr, not raw)
......@@ -799,10 +719,8 @@ void QuadraticIPOptSolver<MatrixType,VectorType>::solve()
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");
}
app->Options()->SetStringValue("jac_d_constant","yes");
app->Options()->SetStringValue("jac_c_constant","yes");
switch (this->verbosity_) {
case NumProc::QUIET:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment