Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/micpsolver/micpsolver.inl
7f13d090d6a1reactmicp/unsaturated

Last Change
Raw File

micpsolver.inl
View Options

/*-------------------------------------------------------
- Module : micpsolver
- File : micpsolver.inl
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "micpsolver.hpp" // for syntaxic coloration...
#include "estimate_cond_number.hpp"
namespace specmicp {
namespace micpsolver {
// Main algorithm
// ##############
template <class Program>
MiCPSolverReturnCode MiCPSolver<Program>::solve(Eigen::VectorXd &x)
{
int cnt = 0;
crashing(x);
MiCPSolverReturnCode retcode = MiCPSolverReturnCode::NotConvergedYet;
Eigen::VectorXd update(get_neq());
while (retcode == MiCPSolverReturnCode::NotConvergedYet)
{
std::cout << "Iteration : " << cnt << std::endl;
setup_residuals(x);
//std::cout << "Residuals \n ------- \n " << m_residuals << "\n ---- \n" << m_phi_residuals << std::endl;
retcode = check_convergence(cnt, update, x);
if (retcode != MiCPSolverReturnCode::NotConvergedYet) break;
setup_jacobian(x);
std::cout << "Solve system " << std::endl;
search_direction_calculation(update);
linesearch(update, x);
std::cout << "Solution : \n ------ \n" << x << std::endl;
++cnt;
}
return retcode;
}
template <class Program>
MiCPSolverReturnCode MiCPSolver<Program>::check_convergence(int nb_iterations,
Eigen::VectorXd& update,
Eigen::VectorXd& solution)
{
MiCPSolverReturnCode termcode = MiCPSolverReturnCode::NotConvergedYet;
if (m_phi_residuals.lpNorm<Eigen::Infinity>() < get_options().fvectol)
{
termcode = MiCPSolverReturnCode::ResidualMinimized;
}
else if (nb_iterations >0 and norm_update<Eigen::Infinity>(update, solution) < get_options().steptol)
{
termcode = MiCPSolverReturnCode::ErrorMinimized;
}
else if (nb_iterations > get_options().max_iter)
{
termcode = MiCPSolverReturnCode::MaxIterations;
}
else if (m_max_taken)
{
++m_consec_max;
if (m_consec_max == get_options().maxiter_maxstep) {
termcode = MiCPSolverReturnCode::MaxStepTakenTooManyTimes;
}
}
else
{
m_consec_max = 0;
}
return termcode;
}
template <class Program>
MiCPSolverReturnCode MiCPSolver<Program>::search_direction_calculation(Eigen::VectorXd& update)
{
m_grad_phi = m_jacobian*m_phi_residuals;
Eigen::VectorXd rscaler(get_neq());
Eigen::VectorXd cscaler(get_neq());
scaling_jacobian(m_jacobian, m_phi_residuals, rscaler, cscaler);
m_jacobian = rscaler.asDiagonal() * m_jacobian * cscaler.asDiagonal();
m_phi_residuals = - rscaler.cwiseProduct(m_phi_residuals);
double descent_cond = 0;
Eigen::ColPivHouseholderQR<Eigen::MatrixXd> solver;
solver.compute(m_jacobian);
if (solver.isInvertible())
{
double cond = estimate_condition_number(solver.matrixR().triangularView<Eigen::Upper>());
int cnt = 0;
while (cond > get_options().condition_limit and cnt < 10)
{
double delta = m_phi_residuals.norm() / 10;
m_jacobian = m_jacobian + delta*Eigen::MatrixXd::Identity(get_neq(), get_neq());
solver.compute(m_jacobian);
cond = estimate_condition_number(solver.matrixR().triangularView<Eigen::Upper>());
++cnt;
}
update = solver.solve(m_phi_residuals);
update = cscaler.cwiseProduct(update);
descent_cond = m_grad_phi.dot(update)
+ get_options().factor_descent_condition * std::pow(update.norm(), get_options().power_descent_condition);
}
std::cout << "Descent condition : " << descent_cond << std::endl;
if (not solver.isInvertible() or descent_cond > 0)
{
update = -m_grad_phi;
}
return MiCPSolverReturnCode::NotConvergedYet;
}
template <class Program>
MiCPSolverReturnCode MiCPSolver<Program>::linesearch(Eigen::VectorXd& update, Eigen::VectorXd& x)
{
double beta = 0.5;
double sigma = 1e-4;
double prev_test = HUGE_VAL;
Eigen::VectorXd xp(get_neq());
int l = 0;
// non monotone linesearch
double merit_value = 0.5*m_phi_residuals.squaredNorm();
std::cout << "Merit value : " << merit_value << std::endl;
double mmax = HUGE_VAL;
if (m_max_merits.size() > 0)
{
mmax = m_max_merits[m_max_merits.size()-1];
}
if (m_max_merits.size() < 4)
{
m_max_merits.push_back(merit_value);
}
else if (merit_value < mmax)
{
m_max_merits[3] = merit_value;
merit_value = mmax;
}
std::cout << "Merit value used : " << merit_value << std::endl;
std::sort(m_max_merits.begin(), m_max_merits.end());
while (l<20)
{
xp = x + std::pow(beta, l)*update;
Eigen::VectorXd new_res(get_neq());
m_program->get_residuals(xp, new_res);
reformulate_residuals_inplace(xp, new_res);
double test = 0.5*new_res.squaredNorm() - merit_value;
if (test > prev_test) {
xp = x + std::pow(beta, l-1)*update;
break;
}
if (test <= sigma*std::pow(beta, l)*m_grad_phi.dot(update))
{
break;
}
prev_test = test;
++l;
}
x = xp;
return MiCPSolverReturnCode::NotConvergedYet;
}
template <class Program>
void MiCPSolver<Program>::crashing(Eigen::VectorXd x)
{
std::cout << "Crashing " << std::endl;
const double beta = 0.5;
const double sigma = 1e-4;
int cnt = 0;
while (cnt < 10)
{
setup_residuals(x);
setup_jacobian(x);
m_grad_phi = m_jacobian*m_phi_residuals;
Eigen::VectorXd xp(get_neq());
int l=0;
while (l<10)
{
xp = x - std::pow(beta, l)*m_grad_phi;
for (int i=get_neq_free(); i<get_neq(); ++i)
{
if (xp(i) < 0) xp(i) = 0;
}
Eigen::VectorXd new_res(get_neq());
compute_residuals(x, new_res);
reformulate_residuals_inplace(x, new_res);
double test = 0.5*(new_res.squaredNorm() - m_phi_residuals.squaredNorm()) + sigma*m_grad_phi.dot(x - xp);
if (test <= 0) break;
++l;
}
if (l == 10) break;
x = xp;
++cnt;
}
std::cout << "Merit function after crashing " << 0.5*m_grad_phi.squaredNorm() << std::endl;
m_max_merits = std::vector<double>(4, 0.5*m_grad_phi.squaredNorm());
}
// Others
// ######
template <class Program>
void MiCPSolver<Program>::reformulate_residuals(const Eigen::VectorXd& x,
const Eigen::VectorXd& r,
Eigen::VectorXd& r_phi)
{
r_phi.block(0, 0, get_neq_free(), 1) = r.block(0, 0, get_neq_free(), 1);
for (int i = get_neq_free(); i<get_neq(); ++i)
{
r_phi(i) = phi_lower_bounded(x(i), r(i), 0);
}
}
template <class Program>
void MiCPSolver<Program>::reformulate_residuals_inplace(const Eigen::VectorXd& x,
Eigen::VectorXd& r)
{
for (int i = get_neq_free(); i<get_neq(); ++i)
{
r(i) = phi_lower_bounded(x(i), r(i), 0);
}
}
double dz(double z) {
if (z > 0)
{
return 1.0;
}
else if (z < 0)
{
return 0;
}
else // z ==0
return 0.5;
}
// page 808 Facchinei and Pang(2003) (vol 2, chap 9)
template <class Program>
void MiCPSolver<Program>::reformulate_jacobian(const Eigen::VectorXd& x,
const Eigen::VectorXd& r,
Eigen::MatrixXd& jacobian
)
{
// setp 1
std::vector<int> beta;
for (int i=get_neq_free(); i<get_neq(); ++i)
{
if (x(i) == 0 and r(i) == 0)
beta.push_back(i);
}
// step 2
Eigen::VectorXd z(get_neq());
z.setZero();
for (auto it=beta.begin(); it!=beta.end(); ++it)
{
z(*it) = 1.0;
}
// modify the jacobian*
const double lambda = get_options().penalization_factor;
for (int i=get_neq_free(); i<get_neq(); ++i)
{
const auto ei = unit_vector(get_neq(), i);
Eigen::VectorXd grad_fi = jacobian.block(i, 0, 1, get_neq()).transpose();
if (x(i) == 0 and r(i) == 0)
{
double gpdotz = grad_fi.dot(z);
double s = std::abs(z(i)) + std::abs(gpdotz);
s = s * std::sqrt(z(i)*z(i)/(s*s) + gpdotz*gpdotz/(s*s));
const double c = lambda*(z(i)/s - 1);
const double d = lambda*(gpdotz/s -1);
grad_fi = c*ei + d*grad_fi;
}
else
{
double s = std::abs(x(i)) + std::abs(r(i));
s = s * std::sqrt(x(i)*x(i)/(s*s) + r(i)*r(i)/(s*s));
const double c = lambda*(x(i)/s - 1) + (1-lambda)*(std::min(0.0, r(i)*dz(x(i))));
const double d = lambda*(r(i)/s -1) + (1-lambda)*(std::min(0.0, x(i))*dz(r(i)));
grad_fi = c*ei + d*grad_fi;
}
jacobian.block(i, 0, 1, get_neq()) = grad_fi.transpose();
}
}
// ref : Munson et al. (2001)
template <class Program>
void MiCPSolver<Program>:: scaling_jacobian(
const Eigen::MatrixXd& jacobian,
const Eigen::VectorXd& r_phi,
Eigen::VectorXd& rscaler,
Eigen::VectorXd& cscaler)
{
for (int i=0; i<get_neq(); ++i)
{
const double sumhsq = jacobian.block(i, 0, 1, get_neq()).array().square().sum();
double s = std::sqrt(r_phi(i)*r_phi(i) + sumhsq);
rscaler(i) = 1.0/std::max(s, 1e-10);
}
for (int i=0; i<get_neq(); ++i)
{
const double sumhsq = (rscaler.asDiagonal()*jacobian).block(0, i, get_neq(), 1).array().square().sum();
double s = std::sqrt(sumhsq);
cscaler(i) = 1.0/std::max(s, 1e-10);
}
}
} // end namespace micpsolver
} // end namespace specmicp
c4science · Help