Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79346536
micpsolverold.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Aug 25, 23:01
Size
10 KB
Mime Type
text/x-c++
Expires
Tue, Aug 27, 23:01 (2 d)
Engine
blob
Format
Raw Data
Handle
20137030
Attached To
rSPECMICP SpecMiCP / ReactMiCP
micpsolverold.hpp
View Options
/*-------------------------------------------------------
- Module : micpsolver
- File : micpsolver.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMIC_MICPSOLVER_MICPSOLVEROLD_HPP
#define SPECMIC_MICPSOLVER_MICPSOLVEROLD_HPP
#include <memory>
#include <Eigen/Dense>
#include "micpsolver_structs.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 MiCPSolverOLD {
public:
//! \brief Constructor
//!
//! @param prog smart pointer toward an instance of a Program to solve
MiCPSolverOLD(std::shared_ptr<Program> prog):
m_program(prog),
m_residuals(prog->total_variables()),
m_phi_residuals(Eigen::VectorXd::Zero(prog->total_variables())),
m_jacobian(prog->total_variables(),prog ->total_variables()),
m_max_taken(false),
m_consec_max(0)
{
}
//! \brief Return a const reference to the options used by the algorithm
const MiCPSolverOptions& get_options() const {return m_options;}
//! \brief Return a reference to the options used by the algorithm
MiCPSolverOptions& get_options() {return m_options;}
//! \brief Set the options
void set_options(const MiCPSolverOptions& options) {m_options = options;}
//! \brief Return the number of equations
int get_neq() const {return m_program->total_variables();}
//! \brief Return the number of equations corresponding to the free variables (size of G)
int get_neq_free() const {return m_program->nb_free_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, get_options().penalization_factor);}
//! \brief Reformulation for a free variable
double phi_free(const double& r) const {
return r;
}
//! \brief Compute the norm of the update
template <int p>
double norm_update(const Eigen::VectorXd& update,
const Eigen::VectorXd& solution) const {
return (update.array().abs()/(solution.array().max(1))
).matrix().template lpNorm<p>();
}
// Residuals and jacobian
// ----------------------
//! \brief Compute the residuals, store it in r
//!
//! @param[in] x the variables
//! @param[out] r vector to store the residuals (must be of the same size as x)
void compute_residuals(const Eigen::VectorXd& x, Eigen::VectorXd& r)
{
m_program->get_residuals(x, r);
get_perf().nb_call_residuals += 1;
}
//! \brief Compute the residuals, use internal storage
//!
//! @param[in] x the variables
void compute_residuals(const Eigen::VectorXd& x)
{
return compute_residuals(x, m_residuals);
get_perf().nb_call_jacobian += 1;
}
//! \brief Compute the jacobian
//!
//! Assumes that the residual have been computed before
//!
//! @param[in] x the variables
void compute_jacobian(Eigen::VectorXd& x)
{
m_program->get_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);
}
//! \brief Compute the factors to scale the jacobian
//!
//! @param[in] jacobian the jacobian to scale (from the reformulated problem)
//! @param[in] residual the residuals corresponding to the jacobian
//! @param[out] rscaler scaling factors of the rows
//! @param[out] cscaler scaling factors of the columns
void scaling_jacobian(const Eigen::MatrixXd& jacobian,
const Eigen::VectorXd& residual,
Eigen::VectorXd& rscaler,
Eigen::VectorXd& cscaler);
// 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 Check for convergence
//!
//! @param nb_iterations the number of iterations
//! @param update the update taken at the previous iteration
//! @param solution the current solution
//! @return a MiCPSolverReturnCode describing the state of the algorithm
MiCPSolverReturnCode check_convergence(int nb_iterations,
Eigen::VectorXd& update,
Eigen::VectorXd& solution);
//! \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[out] update the update to apply to the solution
MiCPSolverReturnCode search_direction_calculation_no_scaling(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);
//! \brief Project variables on the feasible set
void projection(Eigen::VectorXd& x);
//! \brief Return a const reference to an instance of MiCPPerformance
const MiCPPerformance& get_performance() {return m_perf;}
protected:
//! \brief Return a reference to an instance of MiCPPerformance
MiCPPerformance& get_perf() {return m_perf;}
private:
//! \brief Return the step corrected step length if it is too long
double is_step_too_long(Eigen::VectorXd& update);
std::shared_ptr<Program> m_program; //!< Smart pointer of a program
MiCPSolverOptions m_options; //!< The options
MiCPPerformance m_perf; //!< The performance
// 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_max_taken; //!< True if the 'length' of the step was equal or bigger than the maximum step length
int m_consec_max; //!< The number of consecutive step were the step length was equal or bigger than the maximum step length
bool m_gradient_step_taken; //!< True if the update was computed using the gradient
};
} // end namespace micpsolver
} // end namespace specmicp
// ###############//
// Implementation //
// ###############//
#include "micpsolverold.inl"
#endif // SPECMIC_MICPSOLVER_MICPSOLVER_HPP
Event Timeline
Log In to Comment