diff --git a/NSolver/CMakeLists.txt b/NSolver/CMakeLists.txt index f9160f84b85e753f59f0caacb47928ccb7984c2a..82d7d7aaba894ec492859bc9d3175ad9c22589e5 100644 --- a/NSolver/CMakeLists.txt +++ b/NSolver/CMakeLists.txt @@ -12,6 +12,7 @@ SET(my_headers ${CMAKE_CURRENT_SOURCE_DIR}/FiniteElementProblem.hh ${CMAKE_CURRENT_SOURCE_DIR}/GurobiHelper.hh ${CMAKE_CURRENT_SOURCE_DIR}/GUROBISolver.hh + ${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolver.hh ${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolverLean.hh ${CMAKE_CURRENT_SOURCE_DIR}/LeastSquaresProblem.hh ${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraint.hh @@ -58,6 +59,7 @@ SET(my_sources ${CMAKE_CURRENT_SOURCE_DIR}/FiniteElementProblem.cc ${CMAKE_CURRENT_SOURCE_DIR}/GurobiHelper.cc ${CMAKE_CURRENT_SOURCE_DIR}/GUROBISolver.cc + ${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolver.cc ${CMAKE_CURRENT_SOURCE_DIR}/IPOPTSolverLean.cc ${CMAKE_CURRENT_SOURCE_DIR}/LeastSquaresProblem.cc ${CMAKE_CURRENT_SOURCE_DIR}/LinearConstraint.cc diff --git a/NSolver/IPOPTSolver.cc b/NSolver/IPOPTSolver.cc index f1457ba9a160287eebd7f8c8e479ff60be6316a9..a905f61120b54d9bd51411016b0b53293f8f6535 100644 --- a/NSolver/IPOPTSolver.cc +++ b/NSolver/IPOPTSolver.cc @@ -428,771 +428,6 @@ solve(NProblemGmmInterface* _problem, std::vector& _const } -//== IMPLEMENTATION PROBLEM INSTANCE========================================================== - - -void -NProblemIPOPT:: -split_constraints(const std::vector& _constraints) -{ - // split user-provided constraints into general-constraints and bound-constraints - constraints_ .clear(); constraints_.reserve(_constraints.size()); - bound_constraints_.clear(); bound_constraints_.reserve(_constraints.size()); - - for(unsigned int i=0; i<_constraints.size(); ++i) - { - BoundConstraint* bnd_ptr = dynamic_cast(_constraints[i]); - - if(bnd_ptr) - bound_constraints_.push_back(bnd_ptr); - else - constraints_.push_back(_constraints[i]); - } -} - - -//----------------------------------------------------------------------------- - - -void -NProblemIPOPT:: -analyze_special_properties(const NProblemInterface* _problem, const std::vector& _constraints) -{ - hessian_constant_ = true; - jac_c_constant_ = true; - jac_d_constant_ = true; - - if(!_problem->constant_hessian()) - hessian_constant_ = false; - - for(unsigned int i=0; i<_constraints.size(); ++i) - { - if(!_constraints[i]->constant_hessian()) - hessian_constant_ = false; - - if(!_constraints[i]->constant_gradient()) - { - if(_constraints[i]->constraint_type() == NConstraintInterface::NC_EQUAL) - jac_c_constant_ = false; - else - jac_d_constant_ = false; - } - - // nothing else to check? - if(!hessian_constant_ && !jac_c_constant_ && !jac_d_constant_) - break; - } - - //hessian of Lagrangian is only constant, if all hessians of the constraints are zero (due to lambda multipliers) - if(!jac_c_constant_ || !jac_d_constant_) - hessian_constant_ = false; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, - Index& nnz_h_lag, IndexStyleEnum& index_style) -{ - // number of variables - n = problem_->n_unknowns(); - - // number of constraints - m = constraints_.size(); - - // get non-zeros of hessian of lagrangian and jacobi of constraints - nnz_jac_g = 0; - nnz_h_lag = 0; - - // get nonzero structure - std::vector x(n); - problem_->initial_x(P(x)); - - - // nonzeros in the jacobian of C_ and the hessian of the lagrangian - SMatrixNP HP; - SVectorNC g; - SMatrixNC H; - if(!hessian_approximation_) - { - problem_->eval_hessian(P(x), HP); - - // get nonzero structure of hessian of problem - for(int i=0; i= it.col()) - ++nnz_h_lag; - } - - // get nonzero structure of constraints - for( int i=0; ieval_gradient(P(x),g); - - nnz_jac_g += g.nonZeros(); - - if(!hessian_approximation_) - { - // count lower triangular elements - constraints_[i]->eval_hessian (P(x),H); - - SMatrixNC::iterator m_it = H.begin(); - for(; m_it != H.end(); ++m_it) - if( m_it.row() >= m_it.col()) - ++nnz_h_lag; - } - } - - // We use the standard fortran index style for row/col entries - index_style = C_STYLE; - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u, - Index m, Number* g_l, Number* g_u) -{ - // check dimensions - if( n != (Index)problem_->n_unknowns()) - std::cerr << "Warning: IPOPT #unknowns != n " << n << problem_->n_unknowns() << std::endl; - if( m != (Index)constraints_.size()) - std::cerr << "Warning: IPOPT #constraints != m " << m << constraints_.size() << std::endl; - - - // first clear all variable bounds - for( int i=0; iidx()) < n) - { - switch(bound_constraints_[i]->constraint_type()) - { - case NConstraintInterface::NC_LESS_EQUAL: - { - x_u[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound(); - }break; - - case NConstraintInterface::NC_GREATER_EQUAL: - { - x_l[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound(); - }break; - - case NConstraintInterface::NC_EQUAL: - { - x_l[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound(); - x_u[bound_constraints_[i]->idx()] = bound_constraints_[i]->bound(); - }break; - } - } - else - std::cerr << "Warning: invalid bound constraint in IPOPTSolver!!!" << std::endl; - } - - // set bounds for constraints - for( int i=0; iconstraint_type()) - { - case NConstraintInterface::NC_EQUAL : g_u[i] = 0.0 ; g_l[i] = 0.0 ; break; - case NConstraintInterface::NC_LESS_EQUAL : g_u[i] = 0.0 ; g_l[i] = -1.0e19; break; - case NConstraintInterface::NC_GREATER_EQUAL : g_u[i] = 1.0e19; g_l[i] = 0.0 ; break; - default : g_u[i] = 1.0e19; g_l[i] = -1.0e19; break; - } - } - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::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) -{ - // get initial value of problem instance - problem_->initial_x(x); - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_value) -{ - // return the value of the objective function - obj_value = problem_->eval_f(x); - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f) -{ - problem_->eval_gradient(x, grad_f); - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g) -{ - // evaluate all constraint functions - for( int i=0; ieval_constraint(x); - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::eval_jac_g(Index n, const Number* x, bool new_x, - Index m, Index nele_jac, Index* iRow, Index *jCol, - Number* values) -{ - if (values == NULL) - { - // get x for evaluation (arbitrary position should be ok) - std::vector x_rnd(problem_->n_unknowns(), 0.0); - - int gi = 0; - SVectorNC g; - for( int i=0; ieval_gradient(&(x_rnd[0]), g); - SVectorNC::InnerIterator v_it(g); - for( ; v_it; ++v_it) - { - iRow[gi] = i; - jCol[gi] = v_it.index(); - ++gi; - } - } - } - else - { - // return the values of the jacobian of the constraints - - // return the structure of the jacobian of the constraints - // global index - int gi = 0; - SVectorNC g; - - for( int i=0; ieval_gradient(x, g); - - SVectorNC::InnerIterator v_it(g); - - for( ; v_it; ++v_it) - { - values[gi] = v_it.value(); - ++gi; - } - } - - if( gi != nele_jac) - std::cerr << "Warning: number of non-zeros in Jacobian of C is incorrect: " - << gi << " vs " << nele_jac << std::endl; - } - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemIPOPT::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) -{ - if (values == NULL) - { - // return structure - - // get x for evaluation (arbitrary position should be ok) - std::vector x_rnd(problem_->n_unknowns(), 0.0); - - // global index - int gi = 0; - // get hessian of problem - SMatrixNP HP; - problem_->eval_hessian(&(x_rnd[0]), HP); - - for(int i=0; i= it.col()) - { - // it.value(); - iRow[gi] = it.row(); - jCol[gi] = it.col(); - ++gi; - } - } - - // Hessians of Constraints - for(unsigned int j=0; jeval_hessian(&(x_rnd[0]), H); - - SMatrixNC::iterator m_it = H.begin(); - SMatrixNC::iterator m_end = H.end(); - - for(; m_it != m_end; ++m_it) - { - // store lower triangular part only - if( m_it.row() >= m_it.col()) - { - iRow[gi] = m_it.row(); - jCol[gi] = m_it.col(); - ++gi; - } - } - } - - // error check - if( gi != nele_hess) - std::cerr << "Warning: number of non-zeros in Hessian of Lagrangian is incorrect while indexing: " - << gi << " vs " << nele_hess << std::endl; - } - else - { - // return values. - - // global index - int gi = 0; - // get hessian of problem - SMatrixNP HP; - problem_->eval_hessian(x, HP); - - for(int i=0; i= it.col()) - { - values[gi] = obj_factor*it.value(); - ++gi; - } - } - - // Hessians of Constraints - for(unsigned int j=0; jeval_hessian(x, H); - - SMatrixNC::iterator m_it = H.begin(); - SMatrixNC::iterator m_end = H.end(); - - for(; m_it != m_end; ++m_it) - { - // store lower triangular part only - if( m_it.row() >= m_it.col()) - { - values[gi] = lambda[j]*(*m_it); - ++gi; - } - } - } - - // error check - if( gi != nele_hess) - std::cerr << "Warning: number of non-zeros in Hessian of Lagrangian is incorrect2: " - << gi << " vs " << nele_hess << std::endl; - } - return true; -} - - -//----------------------------------------------------------------------------- - - -void NProblemIPOPT::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_value, - const IpoptData* ip_data, - IpoptCalculatedQuantities* ip_cq) -{ - // problem knows what to do - problem_->store_result(x); - - if(store_solution_) - { - x_.resize(n); - for( Index i=0; in_unknowns(); - - // number of constraints - m = constraints_.size(); - - // get nonzero structure - std::vector x(n); - problem_->initial_x(&(x[0])); - // ToDo: perturb x - - // nonzeros in the jacobian of C_ and the hessian of the lagrangian - SMatrixNP HP; - SVectorNC g; - SMatrixNC H; - problem_->eval_hessian(&(x[0]), HP); - nnz_jac_g = 0; - nnz_h_lag = 0; - - // clear old data - jac_g_iRow_.clear(); - jac_g_jCol_.clear(); - h_lag_iRow_.clear(); - h_lag_jCol_.clear(); - - // get non-zero structure of initial hessian - // iterate over rows - for( int i=0; i= (int)v_it.index()) - { - h_lag_iRow_.push_back(i); - h_lag_jCol_.push_back(v_it.index()); - ++nnz_h_lag; - } - } - } - - - // get nonzero structure of constraints - for( int i=0; ieval_gradient(&(x[0]),g); - constraints_[i]->eval_hessian (&(x[0]),H); - - // iterate over sparse vector - SVectorNC::InnerIterator v_it(g); - for(; v_it; ++v_it) - { - jac_g_iRow_.push_back(i); - jac_g_jCol_.push_back(v_it.index()); - ++nnz_jac_g; - } - - // iterate over superSparseMatrix - SMatrixNC::iterator m_it = H.begin(); - SMatrixNC::iterator m_end = H.end(); - for(; m_it != m_end; ++m_it) - if( m_it.row() >= m_it.col()) - { - h_lag_iRow_.push_back(m_it.row()); - h_lag_jCol_.push_back(m_it.col()); - ++nnz_h_lag; - } - } - - // store for error checking... - nnz_jac_g_ = nnz_jac_g; - nnz_h_lag_ = nnz_h_lag; - - // We use the standard fortran index style for row/col entries - index_style = C_STYLE; - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemGmmIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u, - Index m, Number* g_l, Number* g_u) -{ - // first clear all variable bounds - for( int i=0; iconstraint_type()) - { - case NConstraintInterface::NC_EQUAL : g_u[i] = 0.0 ; g_l[i] = 0.0 ; break; - case NConstraintInterface::NC_LESS_EQUAL : g_u[i] = 0.0 ; g_l[i] = -1.0e19; break; - case NConstraintInterface::NC_GREATER_EQUAL : g_u[i] = 1.0e19; g_l[i] = 0.0 ; break; - default : g_u[i] = 1.0e19; g_l[i] = -1.0e19; break; - } - } - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemGmmIPOPT::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) -{ - // get initial value of problem instance - problem_->initial_x(x); - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemGmmIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_value) -{ - // return the value of the objective function - obj_value = problem_->eval_f(x); - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemGmmIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f) -{ - problem_->eval_gradient(x, grad_f); - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemGmmIPOPT::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g) -{ - // evaluate all constraint functions - for( int i=0; ieval_constraint(x); - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemGmmIPOPT::eval_jac_g(Index n, const Number* x, bool new_x, - Index m, Index nele_jac, Index* iRow, Index *jCol, - Number* values) -{ - if (values == NULL) - { - // return the (cached) structure of the jacobian of the constraints - gmm::copy(jac_g_iRow_, VectorPTi(iRow, jac_g_iRow_.size())); - gmm::copy(jac_g_jCol_, VectorPTi(jCol, jac_g_jCol_.size())); - } - else - { - // return the values of the jacobian of the constraints - - // return the structure of the jacobian of the constraints - // global index - int gi = 0; - SVectorNC g; - - for( int i=0; ieval_gradient(x, g); - - SVectorNC::InnerIterator v_it(g); - - for( ; v_it; ++v_it) - { - if(gi < nele_jac) - values[gi] = v_it.value(); - ++gi; - } - } - - if( gi != nele_jac) - std::cerr << "Warning: number of non-zeros in Jacobian of C is incorrect: " - << gi << " vs " << nele_jac << std::endl; - } - - return true; -} - - -//----------------------------------------------------------------------------- - - -bool NProblemGmmIPOPT::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) -{ - if (values == NULL) - { - // return the (cached) structure of the hessian - gmm::copy(h_lag_iRow_, VectorPTi(iRow, h_lag_iRow_.size())); - gmm::copy(h_lag_jCol_, VectorPTi(jCol, h_lag_jCol_.size())); - } - else - { - // return values. - - // global index - int gi = 0; - - // get hessian of problem - problem_->eval_hessian(x, HP_); - - for( int i=0; i= (int)v_it.index()) - { - if( gi < nele_hess) - values[gi] = obj_factor*(*v_it); - ++gi; - } - } - } - - // Hessians of Constraints - for(unsigned int j=0; jeval_hessian(x, H); - - SMatrixNC::iterator m_it = H.begin(); - SMatrixNC::iterator m_end = H.end(); - - for(; m_it != m_end; ++m_it) - { - // store lower triangular part only - if( m_it.row() >= m_it.col()) - { - if( gi < nele_hess) - values[gi] = lambda[j]*(*m_it); - ++gi; - } - } - } - - // error check - if( gi != nele_hess) - std::cerr << "Warning: number of non-zeros in Hessian of Lagrangian is incorrect: " - << gi << " vs " << nele_hess << std::endl; - } - return true; -} - - -//----------------------------------------------------------------------------- - - -void NProblemGmmIPOPT::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_value, - const IpoptData* ip_data, - IpoptCalculatedQuantities* ip_cq) -{ - // problem knows what to do - problem_->store_result(x); -} - - - //============================================================================= } // namespace COMISO //============================================================================= diff --git a/NSolver/NProblemIPOPT.hh b/NSolver/NProblemIPOPT.hh index e6bd44eee63959c74e75064534c18488a07cb48d..9f8e4a6e50220fc9e2a058c012b56089c0c50d3e 100644 --- a/NSolver/NProblemIPOPT.hh +++ b/NSolver/NProblemIPOPT.hh @@ -21,7 +21,6 @@ #include "IPOPTSolverLean.hh" #include "NProblemGmmInterface.hh" #include "NProblemInterface.hh" -#include "NProblemIPOPT.hh" #include "NConstraintInterface.hh" #include "BoundConstraint.hh" #include "CoMISo/Utils/CoMISoError.hh" @@ -321,6 +320,6 @@ private: //============================================================================= #endif // COMISO_IPOPTLEAN_AVAILABLE //============================================================================= -#endif // ACG_IPOPTSOLVER_HH defined +#endif // COMISO_NPROBLEMIPOPT_HH //=============================================================================