Solving a problem is done through several steps:
The following example defines a cost function F and two constraints G0 and G1.
The problem that will be solved in this tutorial is the 71th problem of Hock-Schittkowski: $$min_{x \in \mathbb{R}^4} x_0 x_3 (x_0 + x_1 + x_2) + x_2$$ with the following constraints: $$x_0 x_1 x_2 x_3 \geq 25$$ $$x_0^2 + x_1^2 + x_2^2 + x_3^2 = 40$$ $$1 \leq x_0, x_1, x_2, x_3 \leq 5$$
The library contains the following hierarchy of functions:
roboptim::Functionroboptim::DifferentiableFunctionroboptim::TwiceDifferentiableFunctionroboptim::QuadraticFunctionroboptim::LinearFunctionThese types correspond to dense vectors and matrices relying on Eigen. RobOptim also support sparse matrices and you can even extend the framework to support your own types.
When defining a new function, you have to derive your new function from one of these classes. Depending on the class you derive from, you will have to implement one or several methods:
impl_compute that returns the function's result has
to be defined for all functions.impl_gradient which returns the function's gradient
is to be defined for DifferentiableFunction and its subclasses.impl_hessian for TwiceDifferentiableFunction functions and its subclasses.It is usually recommended to derive from the deepest possible class of the hierarchy (deriving from TwiceDifferentiableFunction is better than DifferentiableFunction).
Keep in mind that the type of the function represents the amount of information the solver will get, not the real nature of a function (it is possible to avoid defining a Hessian by deriving from DifferentiableFunction, even if your function can be derived twice).
In the following sample, a TwiceDifferentiableFunction will be defined.
#include <boost/shared_ptr.hpp>
#include <roboptim/core/twice-differentiable-function.hh>
#include <roboptim/core/io.hh>
#include <roboptim/core/solver.hh>
#include <roboptim/core/solver-factory.hh>
using namespace roboptim;
struct F : public TwiceDifferentiableFunction
{
F () : TwiceDifferentiableFunction (4, 1, "x₀ * x₃ * (x₀ + x₁ + x₂) + x₂")
{
}
void
impl_compute (result_ref result, const_argument_ref x) const
{
result[0] = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];
}
void
impl_gradient (gradient_ref grad, const_argument_ref x, size_type) const
{
grad << x[0] * x[3] + x[3] * (x[0] + x[1] + x[2]),
x[0] * x[3],
x[0] * x[3] + 1,
x[0] * (x[0] + x[1] + x[2]);
}
void
impl_hessian (hessian_ref h, const_argument_ref x, size_type) const
{
h << 2 * x[3], x[3], x[3], 2 * x[0] + x[1] + x[2],
x[3], 0., 0., x[0],
x[3], 0., 0., x[1],
2 * x[0] + x[1] + x[2], x[0], x[0], 0.;
}
};
A constraint is no different from a cost function and can be defined in the same way as a cost function.
The following sample defines two constraints which are twice-differentiable functions.
struct G0 : public TwiceDifferentiableFunction
{
G0 () : TwiceDifferentiableFunction (4, 1, "x₀ * x₁ * x₂ * x₃")
{
}
void
impl_compute (result_ref result, const_argument_ref x) const
{
result[0] = x[0] * x[1] * x[2] * x[3];
}
void
impl_gradient (gradient_ref grad, const_argument_ref x, size_type) const
{
grad << x[1] * x[2] * x[3],
x[0] * x[2] * x[3],
x[0] * x[1] * x[3],
x[0] * x[1] * x[2];
}
void
impl_hessian (hessian_ref h, const_argument_ref x, size_type) const
{
h << 0., x[2] * x[3], x[1] * x[3], x[1] * x[2],
x[2] * x[3], 0., x[0] * x[3], x[0] * x[2],
x[1] * x[3], x[0] * x[3], 0., x[0] * x[1],
x[1] * x[2], x[0] * x[2], x[0] * x[1], 0.;
}
};
struct G1 : public TwiceDifferentiableFunction
{
G1 () : TwiceDifferentiableFunction (4, 1, "x₀² + x₁² + x₂² + x₃²")
{
}
void
impl_compute (result_ref result, const_argument_ref x) const
{
result[0] = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3];
}
void
impl_gradient (gradient_ref grad, const_argument_ref x, size_type) const
{
grad = 2 * x;
}
void
impl_hessian (hessian_ref h, const_argument_ref x, size_type) const
{
h << 2., 0., 0., 0.,
0., 2., 0., 0.,
0., 0., 2., 0.,
0., 0., 0., 2.;
}
};
The last part of this tutorial covers how to build a problem and solve it. The steps are:
int run_test ()
{
typedef Solver<EigenMatrixDense> solver_t;
// Create cost function.
boost::shared_ptr<F> f (new F ());
// Create problem.
solver_t::problem_t pb (f);
// Set bounds for all optimization parameters.
// 1. < x_i < 5. (x_i in [1.;5.])
for (Function::size_type i = 0; i < pb.function ().inputSize (); ++i)
pb.argumentBounds ()[i] = Function::makeInterval (1., 5.);
// Set the starting point.
Function::vector_t start (pb.function ().inputSize ());
start[0] = 1., start[1] = 5., start[2] = 5., start[3] = 1.;
// Create constraints.
boost::shared_ptr<G0> g0 (new G0 ());
boost::shared_ptr<G1> g1 (new G1 ());
F::intervals_t bounds;
solver_t::problem_t::scaling_t scaling;
// Add constraints
bounds.push_back(Function::makeLowerInterval (25.));
scaling.push_back (1.);
pb.addConstraint
(boost::static_pointer_cast<TwiceDifferentiableFunction> (g0),
bounds, scaling);
bounds.clear ();
scaling.clear ();
bounds.push_back(Function::makeInterval (40., 40.));
scaling.push_back (1.);
pb.addConstraint
(boost::static_pointer_cast<TwiceDifferentiableFunction> (g1),
bounds, scaling);
// Initialize solver.
// Here we are relying on a dummy solver.
// You may change this string to load the solver you wish to use:
// - Ipopt: "ipopt", "ipopt-sparse", "ipopt-td"
// - Eigen: "eigen-levenberg-marquardt"
// etc.
// The plugin is built for a given solver type, so choose it adequately.
SolverFactory<solver_t> factory ("dummy-td", pb);
solver_t& solver = factory ();
// Compute the minimum and retrieve the result.
solver_t::result_t res = solver.minimum ();
// Display solver information.
std::cout << solver << std::endl;
// Check if the minimization has succeeded.
// Process the result
switch (res.which ())
{
case solver_t::SOLVER_VALUE:
{
// Get the result.
Result& result = boost::get<Result> (res);
// Display the result.
std::cout << "A solution has been found: " << std::endl
<< result << std::endl;
return 0;
}
case solver_t::SOLVER_ERROR:
{
std::cout << "A solution should have been found. Failing..."
<< std::endl
<< boost::get<SolverError> (res).what ()
<< std::endl;
return 2;
}
case solver_t::SOLVER_NO_SOLUTION:
{
std::cout << "The problem has not been solved yet."
<< std::endl;
return 2;
}
}
// Should never happen.
assert (0);
return 42;
}
This is the last piece of code needed to instantiate and resolve an optimization problem with this package. Running this piece of code with an appropriate solver set, you will get:
Solver:
Problem:
x₀ x₃ (x₀ + x₁ + x₂) + x₂ (twice differentiable function)
Argument's bounds: (1, 5), (1, 5), (1, 5), (1, 5)
Argument's scaling: 1, 1, 1, 1
Number of constraints: 2
Constraint 0
x₀ x₁ x₂ x₃ (twice differentiable function)
Bounds: (25, inf)
Scaling: 1
Initial value: [1](25)
Constraint 1
x₀² + x₁² + x₂² + x₃² (twice differentiable function)
Bounds: (40, 40)
Scaling: 1
Initial value: [1](52) (constraint not satisfied)
Starting point: [4](1,5,5,1)
Starting value: [1](16)
Infinity value (for all functions): inf
A solution has been found:
Result:
Size (input, output): 4, 1
X: [4](1,4.743,3.82115,1.37941)
Value: [1](17.014)
Constraints values: [2](25,40)
Lambda: [2](-0.552294,0.161469)
To see more usage examples, consider looking at the test directory of the project which contains the project test suite.
A tutorial with exercises is also available on GitHub.