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 .. math:: f\left(x_1, x_2, x_3, x_4\right) = x_1x_4(x_1+x_2+x_3) + x_3, subject to the box bounds .. math:: 1 \le x_{1,2,3,4} \le 5 and to the constraints .. math:: \begin{align} x_1^2+x_2^2+x_3^2+x_4^2 - 40 &= 0, \\ 25 - x_1 x_2 x_3 x_4 &\le 0. \end{align} In pagmo's taxonomy, this optimisation problem is * *continuous* (because :math:`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: .. literalinclude:: ../../../../../tutorials/first_udp_ver0.cpp :caption: first_udp_ver0.cpp :language: c++ :linenos: Let's analyse this code. After the inclusion of the necessary header files, and a helpful .. code-block:: c++ using namespace pagmo; to reduce typing, we reach the definition of our optimisation problem via a class called ``problem_v0``. As explained in the :ref:`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 :cpp: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``. 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: .. code-block:: c++ 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: .. code-block:: c++ std::pair 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 :cpp:class:`pagmo::problem` from an instance of ``problem_v0``: .. code-block:: c++ problem p{problem_v0{}}; :cpp:class:`pagmo::problem` is pagmo's :ref:`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), :cpp:class:`pagmo::problem` provides various member functions to access the properties and capabilities of the UDP. We can, for instance, call the :cpp:func:`pagmo::problem::fitness()` member function of ``p`` to invoke the fitness function of the UDP: .. code-block:: c++ 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 :cpp:func:`pagmo::problem::get_lb()` and :cpp:func:`pagmo::problem::get_ub()` member functions: .. code-block:: c++ // 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 .. code-block:: c++ std::cout << p << '\n'; will produce a human-readable summary that may look like this: .. code-block:: none 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: .. literalinclude:: ../../../../../tutorials/first_udp_ver1.cpp :language: c++ :diff: ../../../../../tutorials/first_udp_ver0.cpp 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: .. code-block:: c++ 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 :cpp:type:`~pagmo::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 (:math:`\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, .. math:: \begin{align} x_1^2+x_2^2+x_3^2+x_4^2 - 40 &= 0, \\ 25 - x_1 x_2 x_3 x_4 &\le 0, \end{align} 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: .. code-block:: c++ 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 :cpp:class:`pagmo::problem` from our new UDP, and, if we print it to screen, we can verify how the :cpp:class:`pagmo::problem` class has correctly identified the number and type of constraints from the implementation of our UDP: .. code-block:: none 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 :cpp:func:`pagmo::problem::fitness()` function now produces a vector with three components: .. code-block:: c++ // 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';