Homec4science
Diffusion SpecMiCP / ReactMiCP (chloride) Browse

specmicp/src/micpsolver/micpsolverold.hpp
64b3616eea48chloride

Last Change
Raw File

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
c4science · Help