# Writing your first optimisation problem

## Contents

# Writing your first optimisation problem#

A fundamental requirement of any optimisation library is the ability to represent an optimisation problem to be solved. In this section, our goal is to show how this is done in pagmo.

## A simple problem#

We will be considering the minimisation of the multidimensional objective function

subject to the box bounds

and to the constraints

In pagmo’s taxonomy, this optimisation problem is

*continuous*(because \(x_{1,2,3,4}\) are real variables),*deterministic*(as neither the objectives nor the constraints depend on stochastic variables),*single-objective*(because the objective function produces a single value), and*constrained*(because there are nonlinear constraints in addition to the box bounds).

Note

In this chapter, we will focus on the implementation of a single-objective continuous problem. pagmo can also represent many other types of problems, including multi-objective and (mixed) integer problems. TODO link to the deep dive.

We will start by implementing the objective function and the box bounds first, and we will then add the constraints at a later stage.

## The bare necessities#

In pagmo, an optimisation problem must provide at the very least an objective function and the box bounds. Let’s see how this is done for our simple example problem:

```
1#include <cmath>
2#include <initializer_list>
3#include <iostream>
4#include <utility>
5
6#include <pagmo/problem.hpp>
7#include <pagmo/types.hpp>
8
9using namespace pagmo;
10
11// Our simple example problem, version 0.
12struct problem_v0 {
13 // Implementation of the objective function.
14 vector_double fitness(const vector_double &dv) const
15 {
16 return {dv[0] * dv[3] * (dv[0] + dv[1] + dv[2]) + dv[2]};
17 }
18 // Implementation of the box bounds.
19 std::pair<vector_double, vector_double> get_bounds() const
20 {
21 return {{1., 1., 1., 1.}, {5., 5., 5., 5.}};
22 }
23};
24
25int main()
26{
27 // Construct a pagmo::problem from our example problem.
28 problem p{problem_v0{}};
29
30 // Compute the value of the objective function
31 // in the point (1, 2, 3, 4).
32 std::cout << "Value of the objfun in (1, 2, 3, 4): " << p.fitness({1, 2, 3, 4})[0] << '\n';
33
34 // Fetch the lower/upper bounds for the first variable.
35 std::cout << "Lower bounds: [" << p.get_lb()[0] << "]\n";
36 std::cout << "Upper bounds: [" << p.get_ub()[0] << "]\n\n";
37
38 // Print p to screen.
39 std::cout << p << '\n';
40}
```

Let’s analyse this code.

After the inclusion of the necessary header files, and a helpful

```
using namespace pagmo;
```

to reduce typing, we reach the definition
of our optimisation
problem via a class called `problem_v0`

. As explained
in the section about type erasure,
this class does not need to derive from any base class.
It is just a “regular” class implementing two specific
member functions:

`fitness()`

, which is used to compute the value of the objective function for the input decision vector`dv`

, and`get_bounds()`

, which is used to fetch the box bounds of the problem.

As you can notice, the signatures of both `fitness()`

and `get_bounds()`

use the type `pagmo::vector_double`

. This is the type
used by pagmo to represent decision and fitness vectors, and currently
it is just an alias for `std::vector<double>`

.

The objective function is called `fitness()`

in pagmo because
it is used to return not only the value of the objective function,
but also of the constraints (thus in some sense it computes the
overall “fitness”
of the input decision vector). In this specific case, however,
our optimisation problem does not have constraints yet,
and thus the `fitness()`

implementation just returns
a vector of size 1 whose only element is the value
of the (single) objective function:

```
vector_double fitness(const vector_double &dv) const
{
return {dv[0] * dv[3] * (dv[0] + dv[1] + dv[2]) + dv[2]};
}
```

The other mandatory function, `get_bounds()`

, returns the
box bounds of the problem as a `std::pair`

of lower/upper
bounds:

```
std::pair<vector_double, vector_double> get_bounds() const
{
return {{1., 1., 1., 1.}, {5., 5., 5., 5.}};
}
```

In addition to returning the box bounds, the `get_bounds()`

function plays another important role: it also
(implicitly) establishes the dimension of the problem
via the sizes of the returned lower/upper bounds vectors
(in this specific case, 4).

## Meet pagmo::problem#

After the definition of our optimisation problem, `problem_v0`

,
we encounter the `main()`

function. In the `main()`

, the first
thing we do is to construct a `pagmo::problem`

from an
instance of `problem_v0`

:

```
problem p{problem_v0{}};
```

`pagmo::problem`

is pagmo’s
type-erased interface to optimisation
problems. It is a generic container which can store internally
an instance of any class which “acts like” an optimisation problem,
that is, any class which provides (at least) the two member functions
described earlier (`fitness()`

and `get_bounds()`

). In the pagmo
jargon, we refer to classes which “act like” optimisation problems
as *user-defined problems*, or UDPs.

In addition to storing a UDP (which, by itself, would not be
that useful), `pagmo::problem`

provides various member
functions to access the properties and capabilities of the UDP.
We can, for instance, call the `pagmo::problem::fitness()`

member function of `p`

to invoke the fitness function of the UDP:

```
std::cout << "Value of the objfun in (1, 2, 3, 4): " << p.fitness({1, 2, 3, 4})[0] << '\n';
```

We can also fetch the lower/upper box bounds of the UDP
via the `pagmo::problem::get_lb()`

and
`pagmo::problem::get_ub()`

member functions:

```
// Fetch the lower/upper bounds for the first variable.
std::cout << "Lower bounds: [" << p.get_lb()[0] << "]\n";
std::cout << "Upper bounds: [" << p.get_ub()[0] << "]\n\n";
```

Printing `p`

to screen via

```
std::cout << p << '\n';
```

will produce a human-readable summary that may look like this:

```
Problem name: 10problem_v0
Global dimension: 4
Integer dimension: 0
Fitness dimension: 1
Number of objectives: 1
Equality constraints dimension: 0
Inequality constraints dimension: 0
Lower bounds: [1, 1, 1, 1]
Upper bounds: [5, 5, 5, 5]
Has batch fitness evaluation: false
Has gradient: false
User implemented gradient sparsity: false
Has hessians: false
User implemented hessians sparsity: false
Fitness evaluations: 1
Thread safety: basic
```

Quite a mouthful! Do not worry about deciphering this output right now, as we will examine the more intricate aspects of the definition of an optimisation problem in due time.

For now, let us just point out that, from our simple UDP definition, pagmo was able to infer on its own various properties of the optimisation problem (e.g., the problem dimension, the number of objectives, the absence of constraints, etc.). pagmo is able to do this thanks to both introspection capabilities (based on template metaprogramming) and (hopefully) sensible defaults.

## Adding the constraints#

In order to implement the constraints in our UDP we have to:

add a couple of member functions which describe the type and number of constraints,

modify the fitness function to return, in addition to the objective function, also the value of the constraints for an input decision vector.

Let us see the code:

```
--- /home/runner/work/pagmo2/pagmo2/tutorials/first_udp_ver0.cpp
+++ /home/runner/work/pagmo2/pagmo2/tutorials/first_udp_ver1.cpp
@@ -8,12 +8,26 @@
using namespace pagmo;
-// Our simple example problem, version 0.
-struct problem_v0 {
+// Our simple example problem, version 1.
+struct problem_v1 {
+ // Number of equality constraints.
+ vector_double::size_type get_nec() const
+ {
+ return 1;
+ }
+ // Number of inequality constraints.
+ vector_double::size_type get_nic() const
+ {
+ return 1;
+ }
// Implementation of the objective function.
vector_double fitness(const vector_double &dv) const
{
- return {dv[0] * dv[3] * (dv[0] + dv[1] + dv[2]) + dv[2]};
+ return {
+ dv[0] * dv[3] * (dv[0] + dv[1] + dv[2]) + dv[2], // objfun
+ dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2] + dv[3] * dv[3] - 40., // equality con.
+ 25. - dv[0] * dv[1] * dv[2] * dv[3] // inequality con.
+ };
}
// Implementation of the box bounds.
std::pair<vector_double, vector_double> get_bounds() const
@@ -25,11 +39,14 @@
int main()
{
// Construct a pagmo::problem from our example problem.
- problem p{problem_v0{}};
+ problem p{problem_v1{}};
- // Compute the value of the objective function
- // in the point (1, 2, 3, 4).
- std::cout << "Value of the objfun in (1, 2, 3, 4): " << p.fitness({1, 2, 3, 4})[0] << '\n';
+ // Compute the value of the objective function, equality and
+ // inequality constraints in the point (1, 2, 3, 4).
+ const auto fv = p.fitness({1, 2, 3, 4});
+ std::cout << "Value of the objfun in (1, 2, 3, 4): " << fv[0] << '\n';
+ std::cout << "Value of the eq. constraint in (1, 2, 3, 4): " << fv[1] << '\n';
+ std::cout << "Value of the ineq. constraint in (1, 2, 3, 4): " << fv[2] << '\n';
// Fetch the lower/upper bounds for the first variable.
std::cout << "Lower bounds: [" << p.get_lb()[0] << "]\n";
```

In order to specify the type and number of constraints in our
optimisation problem, we have to implement the two member
functions `get_nec()`

, which returns the number of equality
constraints, and `get_nic()`

, which returns the number of
inequality constraints:

```
vector_double::size_type get_nec() const
{
return 1;
}
vector_double::size_type get_nic() const
{
return 1;
}
```

Note that the number of (in)equality constraints is represented
via the size type of `vector_double`

(which is an unsigned integral type, usually `std::size_t`

).

Next, we need to modify our fitness function to compute, in addition to the objective function, the (in)equality constraints for an input decision vector. pagmo adopts the following conventions:

the constraints are expressed as equations with zero on the right-hand-side,

the values returned by the fitness function are computed from the left-hand-sides of the constraint equations,

the inequality constraints are expressed via a less-than-or-equal relation (\(\leq\)),

in the fitness vector, the equality constraints follow the value(s) of the objective function and precede the inequality constraints.

In our specific example, we have 1 equality constraint and 1 inequality constraint,

and thus the fitness function will have to return a vector with 3 values, which are, in order, the objective function, the equality constraint and the inequality constraint:

```
vector_double fitness(const vector_double &dv) const
{
return {
dv[0] * dv[3] * (dv[0] + dv[1] + dv[2]) + dv[2], // objfun
dv[0] * dv[0] + dv[1] * dv[1] + dv[2] * dv[2] + dv[3] * dv[3] - 40., // equality con.
25. - dv[0] * dv[1] * dv[2] * dv[3] // inequality con.
};
}
```

Now we can create a `pagmo::problem`

from our
new UDP, and, if we print it to screen, we can verify how
the `pagmo::problem`

class has correctly
identified the number and type of constraints from the
implementation of our UDP:

```
Problem name: 10problem_v1
Global dimension: 4
Integer dimension: 0
Fitness dimension: 3
Number of objectives: 1
Equality constraints dimension: 1
Inequality constraints dimension: 1
Tolerances on constraints: [0, 0]
Lower bounds: [1, 1, 1, 1]
Upper bounds: [5, 5, 5, 5]
Has batch fitness evaluation: false
Has gradient: false
User implemented gradient sparsity: false
Has hessians: false
User implemented hessians sparsity: false
Fitness evaluations: 3
Thread safety: basic
```

We can also verify that the `pagmo::problem::fitness()`

function now produces a vector with three components:

```
// Compute the value of the objective function, equality and
// inequality constraints in the point (1, 2, 3, 4).
const auto fv = p.fitness({1, 2, 3, 4});
std::cout << "Value of the objfun in (1, 2, 3, 4): " << fv[0] << '\n';
std::cout << "Value of the eq. constraint in (1, 2, 3, 4): " << fv[1] << '\n';
std::cout << "Value of the ineq. constraint in (1, 2, 3, 4): " << fv[2] << '\n';
```