Page MenuHomec4science
Files F65751161

micpsolver.hpp
No OneTemporary
Actions

File Metadata

Created
Wed, Jun 5, 23:36

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 "ncp_function.hpp"
//! \file micpsolver.hpp The MiCP solver
namespace specmicp {
namespace micpsolver {
inline Eigen::VectorXd unit_vector(int n, int j) {
Eigen::VectorXd unit(n);
unit(j) = 1.0;
return unit;
}
//! Options for the MiCPSolver
struct MiCPSolverOptions
{
int max_iter; //!< Maximum number of iterations allowed
// tolerances
double fvectol; //!< Tolerance for minimization of residuals
double steptol; //!< Tolerance for minimization of error
double condition_limit; //!< Condition number limit
double penalization_factor; //!< Penalization factor for the penalized Fisher-Burmeister function
double maxstep; //!< the maximum step allowed
int maxiter_maxstep; //!< the maximum number of step of length maxstep allowed
double factor_descent_condition; //!< > 0
double power_descent_condition; //!< > 2
MiCPSolverOptions():
max_iter(100),
fvectol(1e-8),
steptol(1e-8),
condition_limit(1e6),
penalization_factor(1.0),
maxstep(1e4),
maxiter_maxstep(5),
factor_descent_condition(1e-4),
power_descent_condition(2.1)
{}
};
enum class MiCPSolverReturnCode
{
LolItsNotSupposedToHappen = -10,
MaxStepTakenTooManyTimes = -4,
FailedToSolveLinearSystem = -3,
MaxIterations = -2,
StationnaryPoint = -1,
NotConvergedYet = 0,
Success = 1,
ResidualMinimized = 2,
ErrorMinimized = 3
};
//! MiCP Solver
//!
//! References :
//! - F. Facchinei and J.-S. Pang.
//! Finite-dimensional variational inequalities and complementarity problems.
//! Springer, New York, 2003.
//! - T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow.
//! The Semismooth Algorithm for Large Scale Complementarity Problems.
//! INFORMS Journal on Computing, 13(4):294-311, 2001.
template <class Program>
class MiCPSolver {
public:
MiCPSolver(std::shared_ptr<Program> prog):
m_program(prog),
m_residuals(prog->total_variables()),
m_phi_residuals(prog->total_variables()),
m_jacobian(prog->total_variables(),prog ->total_variables())
{
}
const MiCPSolverOptions& get_options() const {return m_options;}
MiCPSolverOptions& get_options() {return m_options;}
int get_neq() const {return m_program->total_variables();}
int get_neq_free() const {return m_program->nb_free_variables();}
// Merit function
// ##############
//! First NCP-function
//double phi1(double a, double b) const {return penalized_fisher_burmeister(a, b, get_options().penalization_factor);}
double phi1(double a, double b) const {return fisher_burmeister(a, b);}
//! Reformulation for lower bounded variable
double phi_lower_bounded(const double& x, const double& r, const double& l) const {
return phi1(x-l, r);}
double phi_free(const double& r) const {
return r;
}
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
//! Compute the residuals, store it in r
void compute_residuals(const Eigen::VectorXd& x, Eigen::VectorXd& r)
{
m_program->get_residuals(x, r);
}
//! Compute the residuals, use internal storage
void compute_residuals(const Eigen::VectorXd& x)
{
return compute_residuals(x, m_residuals);
}
//! Compute the jacobian, suppose that x have been computed before
void compute_jacobian(Eigen::VectorXd& x)
{
m_program->get_jacobian(x, m_jacobian);
}
//! Reformulation of the residuals
//!
//! Reformulate the problem - suppose that r has been computed before
void reformulate_residuals(const Eigen::VectorXd& x,
const Eigen::VectorXd& r,
Eigen::VectorXd& r_phi);
void reformulate_residuals_inplace(const Eigen::VectorXd& x,
Eigen::VectorXd &r);
//! Reformulation of the jacobian
//!
//! r is the original vector of residuals (not reformulated)
void reformulate_jacobian(const Eigen::VectorXd& x,
const Eigen::VectorXd& r,
Eigen::MatrixXd& jacobian
);
void scaling_jacobian(const Eigen::MatrixXd& jacobian,
const Eigen::VectorXd& residual,
Eigen::VectorXd& rscaler,
Eigen::VectorXd& cscaler);
// Algorithm
// #########
MiCPSolverReturnCode solve(Eigen::VectorXd& x);
//! Setup the residuals
void setup_residuals(const Eigen::VectorXd& x) {
compute_residuals(x);
reformulate_residuals(x, m_residuals, m_phi_residuals);
}
//! Setup the jacobian
void setup_jacobian(Eigen::VectorXd& x) {
compute_jacobian(x);
reformulate_jacobian(x, m_residuals, m_jacobian);
m_grad_phi = m_jacobian.transpose() * m_phi_residuals;
}
//! Check for convergence
MiCPSolverReturnCode check_convergence(int nb_iterations,
Eigen::VectorXd& update,
Eigen::VectorXd& solution);
MiCPSolverReturnCode search_direction_calculation(Eigen::VectorXd& update);
//! Linesearch
//!
//! Nonmonotone linesearch (see Munson2001)
int linesearch(Eigen::VectorXd& update, Eigen::VectorXd& x);
//! Crashing
//!
//! This function improves, if possible, the initial guess
void crashing(Eigen::VectorXd& x);
//! Project variables on the feasible set
void projection(Eigen::VectorXd& x);
private:
//! Return the step corrected steplength
double is_step_too_long(Eigen::VectorXd& update);
std::shared_ptr<Program> m_program;
MiCPSolverOptions m_options;
// Residuals and jacobian
Eigen::VectorXd m_residuals;
Eigen::VectorXd m_phi_residuals;
Eigen::VectorXd m_grad_phi;
Eigen::MatrixXd m_jacobian;
std::vector<double> m_max_merits;
bool m_max_taken;
int m_consec_max;
};
} // end namespace micpsolver
} // end namespace specmicp
// ###############//
// Implementation //
// ###############//
#include "micpsolver.inl"
#endif // SPECMIC_MICPSOLVER_MICPSOLVER_HPP

Event Timeline

c4science · Help