Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79344968
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, 22:50
Size
10 KB
Mime Type
text/x-c++
Expires
Tue, Aug 27, 22:50 (1 d, 23 h)
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