micpsolver.hpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, May 29, 12:58

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"
	#include <float.h>

	//! \file micpsolver.hpp The MiCP solver

	namespace specmicp {
	namespace micpsolver {

	//! \brief 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
	// misc ...
	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
	bool use_crashing; //! Use a crashing step to improve the starting guess
	double coeff_accept_newton_step; //! Accept the newton step without line search if merit(x+d)<coeff*merit(x^k)
	bool use_scaling; //! Use scaling to improve convergence
	double factor_gradient_search_direction; //! factor which multiply the gradient in the search direction system
	double projection_min_variable; //! threshold under which a mineral is considered to be 0

	//! \brief Constructor - also defines the default value used by the algorithm
	MiCPSolverOptions():
	max_iter(100),
	fvectol(1e-8),
	steptol(1e-10),
	condition_limit(1e6),
	penalization_factor(0.8),
	maxstep(1e4),
	maxiter_maxstep(5),
	factor_descent_condition(1e-4),
	power_descent_condition(2),
	use_crashing(true),
	coeff_accept_newton_step(0.95),
	use_scaling(true),
	factor_gradient_search_direction(5.0),
	projection_min_variable(DBL_EPSILON)
	{}
	};

	//! \brief Return code of the MiCP solver
	//!
	//! A positive return code means that the algorithm converged (but maybe to a stationnary points)
	//! A negative return code means that something wrong happened and was detected
	enum class MiCPSolverReturnCode
	{
	LolItsNotSupposedToHappen = -10, //!< bummer..
	MaxStepTakenTooManyTimes = -4, //!< Propably detect a divergence
	FailedToSolveLinearSystem = -3, //!< Problem in the decomposition... shouldn't be raised, use of the gradient step
	MaxIterations = -2, //!< Algorithm has reached the maximum number of iterations
	StationnaryPoint = -1, //!< Algorithm is tuck in a stationnary point of the merit function
	NotConvergedYet = 0, //!< Keep going...
	Success = 1, //!< Test success (well be careful of stationnary points)
	ResidualMinimized = 2, //!< The residual is minimized (i.e. close to zero)
	ErrorMinimized = 3 //!< The error is minimized, may indicate a stationary point
	};

	//! \brief This structure contains counter to check/query the performance of the algorithm
	//!
	//! It should be updated after each operations, not at the end
	struct MiCPPerformance
	{
	int nb_call_residuals; //! Number of calls of the residual
	int nb_call_jacobian; //! Number of calls of the jacobian
	int nb_factorization; //! Number of factorization performed (may be > number of iterations)
	int nb_gradient_step; //! Number of gradient steps (too many gradient steps indicate a bad starting guess)
	int nb_crashing_iterations; //! Number of crashing iterations
	int nb_iterations; //! Number of iterations
	MiCPSolverReturnCode return_code; //! Return code
	};

	//! \brief The MiCP Solver
	//!
	//! Solve
	//! - \f$ G(u, v) = 0\f$
	//! - \f$ 0 \leq v \perp H(u,v) \geq 0 \f$
	//!
	//! 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);
	get_perf().nb_call_residuals += 1;
	}
	//! Compute the residuals, use internal storage
	void compute_residuals(const Eigen::VectorXd& x)
	{
	return compute_residuals(x, m_residuals);
	get_perf().nb_call_jacobian += 1;
	}
	//! Compute the jacobian, suppose that x have been computed before
	void compute_jacobian(Eigen::VectorXd& x)
	{
	m_program->get_jacobian(x, m_jacobian);
	}
	//! \brief 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);
	//! \brief 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);
	MiCPSolverReturnCode search_direction_calculation_no_scaling(Eigen::VectorXd& update);
	//! \brief 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);

	const MiCPPerformance& get_performance() {return m_perf;}

	protected:
	MiCPPerformance& get_perf() {return m_perf;}

	private:

	//! Return the step corrected steplength
	double is_step_too_long(Eigen::VectorXd& update);

	std::shared_ptr<Program> m_program;

	MiCPSolverOptions m_options;
	MiCPPerformance m_perf;

	// 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;

	bool m_gradient_step_taken;

	};

	} // end namespace micpsolver
	} // end namespace specmicp

	// ###############//
	// Implementation //
	// ###############//
	#include "micpsolver.inl"

	#endif // SPECMIC_MICPSOLVER_MICPSOLVER_HPP

micpsolver.hppNo OneTemporaryActions

File Metadata

micpsolver.hppView Options

Event Timeline

micpsolver.hpp
No OneTemporary
Actions

micpsolver.hpp
View Options