micpsolverold.hpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Jul 8, 18:03

micpsolverold.hpp
View Options

	/*-------------------------------------------------------------------------------

	Copyright (c) 2014,2015 F. Georget <fabieng@princeton.edu>, Princeton University
	All rights reserved.

	Redistribution and use in source and binary forms, with or without modification,
	are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice, this
	list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
	this list of conditions and the following disclaimer in the documentation and/or
	other materials provided with the distribution.

	3. Neither the name of the copyright holder nor the names of its contributors
	may be used to endorse or promote products derived from this software without
	specific prior written permission.

	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
	ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
	WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
	DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
	ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
	(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
	LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
	ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
	SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

	-----------------------------------------------------------------------------*/

	#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

micpsolverold.hppNo OneTemporaryActions

File Metadata

micpsolverold.hppView Options

Event Timeline

micpsolverold.hpp
No OneTemporary
Actions

micpsolverold.hpp
View Options