Quick start

Solving a problem is done through several steps:

• Define your cost function by deriving one kind of function, depending on whether or not you want to provide a Jacobian and/or a Hessian.
• Define your constraints functions in the same manner.
• Build an instance of problem matching your requirements.
• Use one of the solvers to solve your problem.

The following example defines a cost function F and two constraints G0 and G1.

### Problem definition

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$$

### Defining the cost function

The library contains the following hierarchy of functions:

• roboptim::Function
• roboptim::DifferentiableFunction
• roboptim::TwiceDifferentiableFunction
• roboptim::QuadraticFunction
• roboptim::LinearFunction

These 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
{
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.;
}
};


### Defining the constraints

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
{
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
{
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.;
}
};


### Building the problem and solving it

The last part of this tutorial covers how to build a problem and solve it. The steps are:

• Instanciate your functions (cost functions and constraints).
• Pass them to the problem.
• Optional: set a starting point.
• Instanciate a solver which solves your class of problem.
• Solve the problem by calling minimum.


int run_test ()
{
// 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;

bounds.push_back(Function::makeLowerInterval (25.));
scaling.push_back (1.);
(boost::static_pointer_cast<TwiceDifferentiableFunction> (g0),
bounds, scaling);

bounds.clear ();
scaling.clear ();

bounds.push_back(Function::makeInterval (40., 40.));
scaling.push_back (1.);
(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_VALUE_WARNINGS:
{
// Get the result.
ResultWithWarnings& result = boost::get<ResultWithWarnings> (res);

// Display the result.
std::cout << "A solution has been found: " << std::endl
<< result << std::endl;

return 0;
}

case solver_t::SOLVER_NO_SOLUTION:
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;
}
}

// 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.