Commit f026294f authored by David Bommes's avatar David Bommes

added handling of different constraints to CPLEX Solver and added...

added handling of different constraints to CPLEX Solver and added ConeConstraint and added cplex_cone_example

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@207 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 59d1d98b
......@@ -369,3 +369,6 @@ endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_quadratic_resolve_example/CMakeLists.txt" )
add_subdirectory (Examples/small_quadratic_resolve_example)
endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_cplex_soc/CMakeLists.txt" )
add_subdirectory (Examples/small_cplex_soc)
endif()
include (CoMISoExample)
if (WIN32)
acg_add_executable (small_cplex_soc WIN32 ${sources} ${headers} )
else ()
acg_add_executable (small_cplex_soc ${sources} ${headers} )
endif ()
# enable rpath linking
set_target_properties(small_cplex_soc PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
target_link_libraries (small_cplex_soc
CoMISo
${COMISO_LINK_LIBRARIES}
)
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#include <CoMISo/Config/config.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <CoMISo/NSolver/LeastSquaresProblem.hh>
#include <CoMISo/NSolver/LinearConstraint.hh>
#include <CoMISo/NSolver/NPDerivativeChecker.hh>
#include <CoMISo/NSolver/CPLEXSolver.hh>
// solve least squares problem for x=1, y=2 and x-2y+z = 1
// with hard constraints x =-3, z>=3, z^2 >= x^2+y^2
// Example main
int main(void)
{
std::cout << "---------- 1) Problem description..." << std::endl;
std::cout << "Least squares terms: x=1, y=2 and x-2y+z = 1" << std::endl;
std::cout << "Constraints : x =-3, z>=3, z^2 >= x^2+y^2" << std::endl;
std::cout << "---------- 1) Get an instance of a LeastSquaresProblem..." << std::endl;
// number of unknowns
const int n = 3;
COMISO::LeastSquaresProblem lsqp(n);
// term0
COMISO::LinearConstraint::SVectorNC coeffs0(n);
coeffs0.coeffRef(0) = 1.0;
COMISO::LinearConstraint term0(coeffs0,-1.0,COMISO::NConstraintInterface::NC_EQUAL);
lsqp.add_term(&term0);
// term1
COMISO::LinearConstraint::SVectorNC coeffs1(n);
coeffs1.coeffRef(1) = 1.0;
COMISO::LinearConstraint term1(coeffs1,-2.0,COMISO::NConstraintInterface::NC_EQUAL);
lsqp.add_term(&term1);
// term2
COMISO::LinearConstraint::SVectorNC coeffs2(n);
coeffs2.coeffRef(0) = 1.0;
coeffs2.coeffRef(1) = -2.0;
coeffs2.coeffRef(2) = 1.0;
COMISO::LinearConstraint term2(coeffs2,-1.0,COMISO::NConstraintInterface::NC_EQUAL);
lsqp.add_term(&term2);
std::cout << "---------- 2) set up constraints" << std::endl;
// set x = -3.0
COMISO::LinearConstraint lc;
lc.coeffs().coeffRef(0) = 1.0;
lc.b() = 3.0;
// set z>=3 (cone constraint requires that z>=0 !!!)
COMISO::BoundConstraint bc(2,3,3,COMISO::NConstraintInterface::NC_GREATER_EQUAL);
// set z^2 >= x^2+y^2
COMISO::ConeConstraint cc;
cc.resize(3);
cc.i() = 2;
cc.c() = 1.0;
cc.Q()(0,0) = 1.0;
cc.Q()(1,1) = 1.0;
// check if IPOPT solver available in current configuration
#if( COMISO_CPLEX_AVAILABLE)
std::cout << "---------- 3) Get CPLEX solver... " << std::endl;
COMISO::CPLEXSolver cplx;
std::cout << "---------- 4) Solve..." << std::endl;
// fill constraint vector
std::vector<COMISO::NConstraintInterface*> constraints;
constraints.push_back(&lc);
constraints.push_back(&bc);
constraints.push_back(&cc);
cplx.solve(&lsqp, constraints);
#endif
std::cout << "---------- 5) Print solution..." << std::endl;
for( int i=0; i<n; ++i)
std::cerr << "x_" << i << " = " << lsqp.x()[i] << std::endl;
return 0;
}
......@@ -66,6 +66,13 @@ public:
const bool _silent = false) // time limit in seconds
{ std::vector<PairIndexVtype> dc; return solve(_problem, _constraints, dc, _time_limit, _silent);}
// with handling of cone constrints
inline bool solve2(NProblemInterface* _problem, // problem instance
std::vector<NConstraintInterface*>& _constraints, // linear constraints
std::vector<PairIndexVtype>& _discrete_constraints, // discrete constraints
const double _time_limit = 60,
const bool _silent = false); // time limit in seconds
// void set_problem_output_path ( const std::string &_problem_output_path);
// void set_problem_env_output_path( const std::string &_problem_env_output_path);
......
......@@ -10,6 +10,10 @@
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include "CPLEXSolver.hh"
#include "LinearConstraint.hh"
#include "BoundConstraint.hh"
#include "ConeConstraint.hh"
#if COMISO_CPLEX_AVAILABLE
//=============================================================================
......@@ -27,7 +31,7 @@ namespace COMISO {
// ********** SOLVE **************** //
bool
CPLEXSolver::
solve(NProblemInterface* _problem,
solve2(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double _time_limit,
......@@ -257,6 +261,279 @@ solve(NProblemInterface* _problem,
}
//-----------------------------------------------------------------------------
bool
CPLEXSolver::
solve(NProblemInterface* _problem,
std::vector<NConstraintInterface*>& _constraints,
std::vector<PairIndexVtype>& _discrete_constraints,
const double _time_limit,
const bool _silent)
{
try
{
//----------------------------------------------
// 0. set up environment
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> get environment...\n";
IloEnv env_;
if(!_silent)
std::cerr << "cplex -> get model...\n";
IloModel model(env_);
// model.getEnv().set(GRB_DoubleParam_TimeLimit, _time_limit);
//----------------------------------------------
// 1. allocate variables and initialize limits
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> allocate variables...\n";
// determine variable types: 0->real, 1->integer, 2->bool
std::vector<char> vtypes (_problem->n_unknowns(),0);
for(unsigned int i=0; i<_discrete_constraints.size(); ++i)
if(_discrete_constraints[i].first < vtypes.size())
{
switch(_discrete_constraints[i].second)
{
case Integer: vtypes [_discrete_constraints[i].first] = 1;
break;
case Binary : vtypes[_discrete_constraints[i].first] = 2;
break;
default : break;
}
}
else
std::cerr << "ERROR: requested a discrete variable which is above the total number of variables"
<< _discrete_constraints[i].first << " vs " << vtypes.size() << std::endl;
// CPLEX variables
std::vector<IloNumVar> vars;
// first all
for( int i=0; i<_problem->n_unknowns(); ++i)
switch(vtypes[i])
{
case 0 : vars.push_back( IloNumVar(env_,-IloInfinity, IloInfinity, IloNumVar::Float) ); break;
case 1 : vars.push_back( IloNumVar(env_, -IloIntMax, IloIntMax, IloNumVar::Int) ); break;
case 2 : vars.push_back( IloNumVar(env_, 0, 1, IloNumVar::Bool) ); break;
}
// Integrate new variables
// model.update();
//----------------------------------------------
// 2. setup constraints
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> setup constraints...\n";
// get zero vector
std::vector<double> x(_problem->n_unknowns(), 0.0);
// handle constraints depending on their tyep
for(unsigned int i=0; i<_constraints.size(); ++i)
{
// is linear constraint?
LinearConstraint* lin_ptr = dynamic_cast<LinearConstraint*>(_constraints[i]);
if(lin_ptr)
{
IloExpr lin_expr(env_);
NConstraintInterface::SVectorNC gc;
_constraints[i]->eval_gradient(P(x), gc);
NConstraintInterface::SVectorNC::InnerIterator v_it(gc);
for(; v_it; ++v_it)
lin_expr += vars[v_it.index()]*v_it.value();
double b = _constraints[i]->eval_constraint(P(x));
switch(_constraints[i]->constraint_type())
{
case NConstraintInterface::NC_EQUAL : model.add(lin_expr + b == 0); break;
case NConstraintInterface::NC_LESS_EQUAL : model.add(lin_expr + b <= 0); break;
case NConstraintInterface::NC_GREATER_EQUAL : model.add(lin_expr + b >= 0); break;
}
}
else
{
BoundConstraint* bnd_ptr = dynamic_cast<BoundConstraint*>(_constraints[i]);
if(bnd_ptr)
{
switch(bnd_ptr->constraint_type())
{
case NConstraintInterface::NC_EQUAL : vars[bnd_ptr->idx()].setBounds(bnd_ptr->bound(), bnd_ptr->bound()); break;
case NConstraintInterface::NC_LESS_EQUAL : vars[bnd_ptr->idx()].setUB(bnd_ptr->bound()); break;
case NConstraintInterface::NC_GREATER_EQUAL : vars[bnd_ptr->idx()].setLB(bnd_ptr->bound()); break;
}
}
else
{
ConeConstraint* cone_ptr = dynamic_cast<ConeConstraint*>(_constraints[i]);
if(cone_ptr)
{
IloExpr soc_lhs(env_);
IloExpr soc_rhs(env_);
soc_rhs= 0.5*cone_ptr->c()*vars[cone_ptr->i()]*vars[cone_ptr->i()];
NConstraintInterface::SMatrixNC::iterator q_it = cone_ptr->Q().begin();
for(; q_it != cone_ptr->Q().end(); ++q_it)
{
soc_lhs += 0.5*(*q_it)*vars[q_it.col()]*vars[q_it.row()];
}
model.add(soc_lhs <= soc_rhs);
}
else
std::cerr << "Warning: CPLEXSolver received a constraint of unknow type!!! -> skipping it" << std::endl;
}
}
}
// model.update();
//----------------------------------------------
// 3. setup energy
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> setup energy...\n";
if(!_problem->constant_hessian())
std::cerr << "Warning: CPLEXSolver received a problem with non-constant hessian!!!" << std::endl;
// GRBQuadExpr objective;
IloExpr objective(env_);
// add quadratic part
NProblemInterface::SMatrixNP H;
_problem->eval_hessian(P(x), H);
for( int i=0; i<H.outerSize(); ++i)
for (NProblemInterface::SMatrixNP::InnerIterator it(H,i); it; ++it)
objective += 0.5*it.value()*vars[it.row()]*vars[it.col()];
// add linear part
std::vector<double> g(_problem->n_unknowns());
_problem->eval_gradient(P(x), P(g));
for(unsigned int i=0; i<g.size(); ++i)
objective += g[i]*vars[i];
// add constant part
objective += _problem->eval_f(P(x));
model.add(IloMinimize(env_,objective));
// model.set(GRB_IntAttr_ModelSense, 1);
// model.setObjective(objective);
// model.update();
//----------------------------------------------
// 4. solve problem
//----------------------------------------------
if(!_silent)
std::cerr << "cplex -> generate model...\n";
IloCplex cplex(model);
cplex.setParam(IloCplex::TiLim, _time_limit);
{ // hack
// 0 [CPX_NODESEL_DFS] Depth-first search
// 1 [CPX_NODESEL_BESTBOUND] Best-bound search
// 2 [CPX_NODESEL_BESTEST] Best-estimate search
// 3 [CPX_NODESEL_BESTEST_ALT] Alternative best-estimate search
// cplex.setParam(IloCplex::NodeSel , 0);
}
if(!_silent)
std::cerr << "cplex -> solve...\n";
// silent mode?
if(_silent)
cplex.setOut(env_.getNullStream());
IloBool solution_found = cplex.solve();
// if (solution_input_path_.empty())
// {
// if (!problem_env_output_path_.empty())
// {
// std::cout << "Writing problem's environment into file \"" << problem_env_output_path_ << "\"." << std::endl;
// model.getEnv().writeParams(problem_env_output_path_);
// }
// if (!problem_output_path_.empty())
// {
// std::cout << "Writing problem into file \"" << problem_output_path_ << "\"." << std::endl;
// GurobiHelper::outputModelToMpsGz(model, problem_output_path_);
// }
//
// model.optimize();
// }
// else
// {
// std::cout << "Reading solution from file \"" << solution_input_path_ << "\"." << std::endl;
// }
//
//----------------------------------------------
// 5. store result
//----------------------------------------------
if(solution_found != IloFalse)
{
if(!_silent)
std::cerr << "cplex -> store result...\n";
// store computed result
for(unsigned int i=0; i<vars.size(); ++i)
x[i] = cplex.getValue(vars[i]);
_problem->store_result(P(x));
}
/*
if (solution_input_path_.empty())
{
// store computed result
for(unsigned int i=0; i<vars.size(); ++i)
x[i] = vars[i].get(GRB_DoubleAttr_X);
}
*/
// else
// {
// std::cout << "Loading stored solution from \"" << solution_input_path_ << "\"." << std::endl;
// // store loaded result
// const size_t oldSize = x.size();
// x.clear();
// GurobiHelper::readSolutionVectorFromSOL(x, solution_input_path_);
// if (oldSize != x.size()) {
// std::cerr << "oldSize != x.size() <=> " << oldSize << " != " << x.size() << std::endl;
// throw std::runtime_error("Loaded solution vector doesn't have expected dimension.");
// }
// }
//
// _problem->store_result(P(x));
//
// // ObjVal is only available if the optimize was called.
// if (solution_input_path_.empty())
// std::cout << "GUROBI Objective: " << model.get(GRB_DoubleAttr_ObjVal) << std::endl;
return solution_found;
}
catch (IloException& e)
{
cerr << "Concert exception caught: " << e << endl;
return false;
}
catch (...)
{
cerr << "Unknown exception caught" << endl;
return false;
}
return false;
}
//-----------------------------------------------------------------------------
......
//=============================================================================
//
// CLASS ConerConstraint
//
//=============================================================================
#ifndef COMISO_CONECONSTRAINT_CC
#define COMISO_CONECONSTRAINT_CC
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include "NConstraintInterface.hh"
#include "ConeConstraint.hh"
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== Implementation =========================================================
/// Default constructor
ConeConstraint::ConeConstraint()
: NConstraintInterface(NConstraintInterface::NC_GREATER_EQUAL)
{
Q_.clear();
i_ = 1.0;
}
// cone constraint of the form -> 0.5*(c_ * x(i_)^2 - x^T Q_ x) >= 0
ConeConstraint::ConeConstraint(const double _c, const int _i, const SMatrixNC& _Q)
: NConstraintInterface(NConstraintInterface::NC_GREATER_EQUAL),
c_(_c), i_(_i), Q_(_Q)
{
}
/// Destructor
ConeConstraint::~ConeConstraint() {}
int ConeConstraint::n_unknowns()
{
return Q_.cols();
}
void ConeConstraint::resize(const unsigned int _n)
{
Q_.resize(_n,_n);
}
void ConeConstraint::clear()
{
Q_.clear();
}
const ConeConstraint::SMatrixNC& ConeConstraint::Q() const
{
return Q_;
}
ConeConstraint::SMatrixNC& ConeConstraint::Q()
{
return Q_;
}
const int& ConeConstraint::i() const
{
return i_;
}
int& ConeConstraint::i()
{
return i_;
}
const double& ConeConstraint::c() const
{
return c_;
}
double& ConeConstraint::c()
{
return c_;
}
double ConeConstraint::eval_constraint ( const double* _x )
{
// cone constraint of the form -> 0.5*(c_ * x(i_)^2 - x^T Q_ x) >= 0
double v = c_*_x[i_]*_x[i_];
SMatrixNC::iterator m_it = Q_.begin();
SMatrixNC::iterator m_end = Q_.end();
for(; m_it != m_end; ++m_it)
{
v -= (*m_it)*_x[m_it.row()]*_x[m_it.col()];
}
return 0.5*v;
}
void ConeConstraint::eval_gradient( const double* _x, SVectorNC& _g )
{
_g.setZero();
_g.resize(Q_.rows());
SMatrixNC::iterator m_it = Q_.begin();
SMatrixNC::iterator m_end = Q_.end();
for(; m_it != m_end; ++m_it)
{
_g.coeffRef(m_it.row()) -= (*m_it)*_x[m_it.col()];
}
_g.coeffRef(i_) += c_*_x[i_];
}
void ConeConstraint::eval_hessian ( const double* _x, SMatrixNC& _h )
{
_h = Q_;
Q_(i_,i_) += c_;
}
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // ACG_ConeConstraint_HH defined
//=============================================================================
//=============================================================================
//
// CLASS ConeConstraint
//
//=============================================================================
#ifndef COMISO_CONECONSTRAINT_HH
#define COMISO_CONECONSTRAINT_HH
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include "NConstraintInterface.hh"
#include <Eigen/StdVector>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
class COMISODLLEXPORT ConeConstraint : public NConstraintInterface
{
public:
// sparse vector type
typedef NConstraintInterface::SVectorNC SVectorNC;
typedef NConstraintInterface::SMatrixNC SMatrixNC;
/// Default constructor
ConeConstraint();
// cone constraint of the form -> 0.5*(c_ * x(i_)^2 - x^T Q_ x) >= 0
ConeConstraint(const double _c, const int _i, const SMatrixNC& _Q);
/// Destructor
virtual ~ConeConstraint();
virtual int n_unknowns();
// resize coefficient vector = #unknowns
void resize(const unsigned int _n);
// clear to zero constraint 0 =_type 0
void clear();
const double& c() const;
double& c();
const int& i() const;