Page MenuHomec4science
Files F60199853

micpsolver.hpp
No OneTemporary
Actions

File Metadata

Created
Sun, Apr 28, 06:13

micpsolver.hpp
View Options

/*-------------------------------------------------------
- Module : micpsolver
- File : micpsolver.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMIC_MICPSOLVER_MICPSOLVER_HPP
#define SPECMIC_MICPSOLVER_MICPSOLVER_HPP
#include <memory>
#include <Eigen/Dense>
#include "micpsolver_base.hpp"
#include "ncp_function.hpp"
//! \file micpsolver.hpp The MiCP solver
namespace specmicp {
namespace micpsolver {
//! \brief Call a NCP-function
//!
//! @tparam ncp_t the NCP function to use
//! @param a first argument
//! @param b second argument
//! @param t parameter of the NCP function
//!
template <NCPfunction ncp_t>
inline double ncp_function(double a, double b, double t);
//! \brief Reformulate the jacobian using the NCP-function ncp_t
//!
//! @tparam ncp_t the NCP-function to use
//! @param[in] neq number of equation
//! @param[in] neq_free number of free variables
//! @param[in] x the variables
//! @param[in] r the residuals
//! @param[in,out] jacobian the jacobian to reformulate
//! @param[in,out] r_reformulated the reformulated residuals
//! @param[in] t parameter of the NCP function
template <NCPfunction ncp_t>
inline void reformulate_jacobian_helper(
int neq,
int neq_free,
const Eigen::VectorXd& x,
const Eigen::VectorXd& r,
Eigen::MatrixXd& jacobian,
Eigen::VectorXd& r_reformulated,
double t
);
//! \brief reformulate the result of the Newton system
template <NCPfunction ncp_t>
inline void reformulate_result(
int neq,
int neq_free,
Eigen::VectorXd& x,
const Eigen::VectorXd& orig_r,
Eigen::VectorXd& grad_phi,
Eigen::VectorXd& update
);
//! \brief The MiCP Solver
//!
//! Solve
//! - \f$ G(u, v) = 0\f$
//! - \f$ 0 \leq v \perp H(u,v) \geq 0 \f$
//!
//! \tparam Program a subclass of MiCPProg
//!
//! References :
//! - \cite Munson2001
//! - \cite Facchinei2003
//!
template <class Program, NCPfunction ncp_t>
class MiCPSolver: public MiCPSolverBaseProgram<Program>
{
public:
using base = MiCPSolverBaseProgram<Program>;
using MiCPSolverBaseProgram<Program>::get_options;
using MiCPSolverBaseProgram<Program>::get_program;
using MiCPSolverBaseProgram<Program>::get_perf;
using MiCPSolverBaseProgram<Program>::get_neq;
using MiCPSolverBaseProgram<Program>::get_neq_free;
//! \brief Constructor
//!
//! @param prog smart pointer toward an instance of a Program to solve
MiCPSolver(std::shared_ptr<Program> prog):
MiCPSolverBaseProgram<Program>(prog),
m_residuals(prog->total_variables()),
m_phi_residuals(Eigen::VectorXd::Zero(prog->total_variables())),
m_jacobian(prog->total_variables(),prog ->total_variables())
{
}
// Merit function
// ##############
//! \brief Reformulation for lower bounded variable
double phi_lower_bounded(const double& x, const double& r, const double& l) const {
return ncp_function<ncp_t>(x-l, r, base::get_options().penalization_factor);}
//! \brief Reformulation for a free variable
double phi_free(const double& r) const {
return r;
}
// Residuals and jacobian
// ----------------------
//! \brief Compute the residuals, use internal storage
//!
//! @param[in] x the variables
void compute_residuals(const Eigen::VectorXd& x)
{
base::compute_residuals(x, m_residuals);
}
//! \brief Compute the jacobian
//!
//! Assumes that the residual have been computed before
//!
//! @param[in] x the variables
void compute_jacobian(Eigen::VectorXd& x)
{
base::compute_jacobian(x, m_jacobian);
}
//! \brief Reformulation of the residuals
//!
//! Reformulate the problem - assumes that r, the residual, has been computed before
//!
//! @param[in] x the variables
//! @param[in] r the residuals
//! @param[out] r_phi a vector of size neq, which will contain the reformulated residuals
void reformulate_residuals(const Eigen::VectorXd& x,
const Eigen::VectorXd& r,
Eigen::VectorXd& r_phi);
//! \brief Reformulation of the residuals *inplace*
//!
//! @param[in] x the variables
//! @param[in,out] r the residual, will contain the reformulated residuals
void reformulate_residuals_inplace(const Eigen::VectorXd& x,
Eigen::VectorXd &r);
//! \brief Reformulation of the jacobian
//!
//! r is the original vector of residuals (not reformulated)
void reformulate_jacobian(const Eigen::VectorXd& x
)
{
reformulate_jacobian_helper<ncp_t>(get_neq(), get_neq_free(),
x, m_residuals,
m_jacobian,
m_phi_residuals,
get_options().penalization_factor);
}
// Algorithm
// #########
//! \brief Solver the program using x as initial guess
//!
//! @param[in,out] x the initial guess, as output contains the solution (from the last iteration)
MiCPSolverReturnCode solve(Eigen::VectorXd& x);
//! \brief Setup the residuals
//!
//! @param[in] x the current solution
void setup_residuals(const Eigen::VectorXd& x) {
compute_residuals(x);
reformulate_residuals(x, m_residuals, m_phi_residuals);
}
//! \brief Setup the jacobian
//!
//! @param[in] x the current solution
void setup_jacobian(Eigen::VectorXd& x) {
compute_jacobian(x);
reformulate_jacobian(x);
m_grad_phi = m_jacobian.transpose() * m_phi_residuals;
}
//! \brief Solve the Newton system
//!
//! Assume that the Newton system has been formed
//!
//! \param[out] update the update to apply to the solution
MiCPSolverReturnCode search_direction_calculation(Eigen::VectorXd& update);
//! \brief Solve the Newton system - does not scale the jacobian
//!
//! Assume that the Newton system has been formed
//!
//! \param[in,out] x the variables
//! \param[out] update the update to apply to the solution
MiCPSolverReturnCode search_direction_calculation_no_scaling(Eigen::VectorXd &x, Eigen::VectorXd& update);
//! \brief Linesearch
//!
//! Nonmonotone linesearch
//!
//! References :
//! - \cite Dennis1983
//! - \cite Munson2001
//!
int linesearch(Eigen::VectorXd& update, Eigen::VectorXd& x);
//! \brief Crashing
//!
//! This function improves, if possible, the initial guess
//! Reference :
//! - \cite Munson2001
void crashing(Eigen::VectorXd& x);
private:
// Residuals and jacobian
Eigen::VectorXd m_residuals; //!< The residuals
Eigen::VectorXd m_phi_residuals; //!< The reformulated residuals
Eigen::VectorXd m_grad_phi; //!< The gradient of the reformulated residuals
Eigen::MatrixXd m_jacobian; //!< The jacobian
std::vector<double> m_max_merits; //!< Contains the m best value of the merit function
bool m_gradient_step_taken; //!< True if the update was computed using the gradient
};
} // end namespace micpsolver
} // end namespace specmicp
// ###############//
// Implementation //
// ###############//
#include "micpsolver.inl"
#endif // SPECMIC_MICPSOLVER_MICPSOLVER_HPP

Event Timeline

c4science · Help