Homec4science
Diffusion SpecMiCP / ReactMiCP (gasphase) Browse

specmicp/tests/micpsolver/test_micpsolver.hpp
761f64c2833agasphase

Last Change
Raw File

test_micpsolver.hpp
View Options

/*-------------------------------------------------------
- Module : tests/micpsolver
- File : test_micpsolver.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "cxxtest/TestSuite.h"
#include "micpsolver/micpsolver.hpp"
#include "micpsolver/micpprog.hpp"
using namespace specmicp::micpsolver;
class TestLinearProgram: public MiCPProg<TestLinearProgram>
{
public:
//! Return the number of variables
int total_variables() {return 3;}
int nb_free_variables() {return 2;}
int nb_complementarity_variables() {return 1;}
//! Return the residual (R(X))
void get_residuals(const types::vector_t& x,
types::vector_t& residual) {
assert(residual.rows() == 3);
residual << 1 + x(0) + x(2), 1 + 1*x(1) + x(2), -x(0) - x(1) +x(2);
}
//! Return the jacobian (J(x))
void get_jacobian(const types::vector_t& x,
types::matrix_t& jacobian)
{
assert(jacobian.cols() == 3);
assert(jacobian.rows() == 3);
jacobian << 1, 0, 1, 0, 1, 1, -1, -1, 1;
}
};
class TestNonLinearProgram: public MiCPProg<TestNonLinearProgram>
{
public:
//! Return the number of variables
int total_variables() {return 3;}
int nb_free_variables() {return 2;}
int nb_complementarity_variables() {return 1;}
//! Return the residual (R(X))
void get_residuals(const types::vector_t& x,
types::vector_t& residual) {
assert(residual.rows() == 3);
residual << -1 + x(0)*x(0) + x(2), 1 + 1*x(1) + x(2)*x(2), -x(0) - x(1) +x(2);
}
//! Return the jacobian (J(x))
void get_jacobian(const types::vector_t& x,
types::matrix_t& jacobian)
{
assert(jacobian.cols() == 3);
assert(jacobian.rows() == 3);
jacobian << 2*x(0), 0, 1, 0, 1, 2*x(2), -1, -1, 1;
}
};
class TestKojimaProgram: public MiCPProg<TestKojimaProgram>
{
public:
//! Return the number of variables
int total_variables() {return 4;}
int nb_free_variables() {return 0;}
int nb_complementarity_variables() {return 4;}
//! Return the residual (R(X))
void get_residuals(const types::vector_t& x,
types::vector_t& residual) {
assert(residual.rows() == 4);
residual << 3*x(0)*x(0)+2*x(0)*x(1)+2*x(1)*x(1)+x(2)+3*x(3)-6,
2*x(0)*x(0)+x(1)*x(1)+x(0)+10*x(2) +2*x(3)-2,
3*x(0)+x(0)*x(1)+2*x(1)*x(1)+2*x(2)+9*x(3)-9,
x(0)*x(0)+3*x(1)*x(1)+2*x(2)+3*x(3)-3
;
}
//! Return the jacobian (J(x))
void get_jacobian(types::vector_t& x,
types::matrix_t& jacobian)
{
assert(jacobian.cols() == 4);
assert(jacobian.rows() == 4);
const int neq = total_variables();
Eigen::VectorXd res(total_variables());
Eigen::VectorXd perturbed_res(total_variables());
get_residuals(x, res);
for (int j=0; j<neq; ++j)
{
double h = 1e-8*std::abs(x(j));
//h = std::copysign(h, x(j));
if (h==0) h = 1e-8;
double tmp = x(j);
x(j) += h;
h = x(j) - tmp;
get_residuals(x, perturbed_res);
for (int i=0; i<neq; ++i)
{
jacobian(i, j) = (perturbed_res(i) - res(i))/h;
}
x(j) = tmp;
}
return;
}
};
class Test_Micp : public CxxTest::TestSuite
{
public:
void test_compiletest()
{
}
void test_linear_program_0()
{
std::shared_ptr<TestLinearProgram> ptrprog = std::make_shared<TestLinearProgram>();
MiCPSolver<TestLinearProgram> solver(ptrprog);
Eigen::VectorXd x(3);
x << -1.5, -2 , 0;
MiCPSolverReturnCode ret = solver.solve(x);
TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8);
}
void test_linear_program_10()
{
std::shared_ptr<TestLinearProgram> ptrprog = std::make_shared<TestLinearProgram>();
MiCPSolver<TestLinearProgram> solver(ptrprog);
Eigen::VectorXd x(3);
x << -1.5, -2 , 10;
MiCPSolverReturnCode ret = solver.solve(x);
TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8);
}
void test_nonlinear_program_0()
{
std::shared_ptr<TestNonLinearProgram> ptrprog = std::make_shared<TestNonLinearProgram>();
MiCPSolver<TestNonLinearProgram> solver(ptrprog);
Eigen::VectorXd x(3);
x << -1.5, -2 , 0;
MiCPSolverReturnCode ret = solver.solve(x);
std::cout << x << std::endl;
TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8);
}
void test_nonlinear_program_5()
{
std::shared_ptr<TestNonLinearProgram> ptrprog = std::make_shared<TestNonLinearProgram>();
MiCPSolver<TestNonLinearProgram> solver(ptrprog);
Eigen::VectorXd x(3);
x << -1.5, -2 , 5;
MiCPSolverReturnCode ret = solver.solve(x);
std::cout << x << std::endl;
TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8);
TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8);
}
void test_kojima_program()
{
std::shared_ptr<TestKojimaProgram> ptrprog = std::make_shared<TestKojimaProgram>();
MiCPSolver<TestKojimaProgram> solver(ptrprog);
Eigen::VectorXd x(4);
x << 0.9, 0.1 , 2.9, 0.1;
MiCPSolverReturnCode ret = solver.solve(x);
std::cout << x << std::endl;
TS_ASSERT_EQUALS(ret, MiCPSolverReturnCode::ResidualMinimized);
//TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(0)+1), 1e-8);
//TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(1)+1), 1e-8);
//TS_ASSERT_LESS_THAN_EQUALS(std::abs(x(2)), 1e-8);
}
};
c4science · Help