In this Tutorial we will learn how to code simple optimization problems (continuous, single objective, unconstrained), so that PyGMO can then apply all of its algorithmic power to solve it. In a nutshell, we will write a class deriving from PyGMO.problem.base and reimplement some of its ‘virtual’ methods.
Let us start with defining one of the classic textbook examples of an optimization problem.
In [ ]:
from PyGMO.problem import base
class my_problem(base):
"""
De Jong (sphere) function implemented purely in Python.
USAGE: my_problem(dim=10)
* dim problem dimension
"""
def __init__(self, dim=10):
# First we call the constructor of the base class telling PyGMO
# what kind of problem to expect ('dim' dimensions, 1 objective, 0 contraints etc.)
super(my_problem, self).__init__(dim)
# We set the problem bounds (in this case equal for all components)
self.set_bounds(-5.12, 5.12)
# Reimplement the virtual method that defines the objective function.
def _objfun_impl(self, x):
# Compute the sphere function
f = sum([x[i] ** 2 for i in range(self.dimension)])
# Note that we return a tuple with one element only. In PyGMO the objective functions
# return tuples so that multi-objective optimization is also possible.
return (f, )
# Finally we also reimplement a virtual method that adds some output to the __repr__ method
def human_readable_extra(self):
return "\n\t Problem dimension: " + str(self.__dim)
Note that by default PyGMO will assume one wants to minimize the objective function. In the second part of this tutorial we will also see how it is possible to change this default behaviour.
To solve our problem we will use Artificial Bee Colony algorithm with 20 individuals.
In [ ]:
from PyGMO import algorithm, island
prob = my_problem(dim=10) # Create a 10-dimensional problem
algo = algorithm.bee_colony(gen=500) # 500 generations of bee_colony algorithm
isl = island(algo, prob, 20) # Instantiate population with 20 individuals
isl.evolve(1) # Evolve the island once
isl.join()
print(isl.population.champion.f)
And we are done! Objective value in the order of $10^{-25}$, no big deal for a sphere problem.
Let’s consider now a maximization problem. To solve such a problem, two possibilities are available to the PaGMO/PyGMO user. The first one is to code the original problem as a minimization problem by premultiplying the objective function by $-1$ (a technique wich is often used and requires no particular effort). If such a method is used, the final fitness value obtained with PyGMO has to be multiplied by $-1$ to get back to the correct value.
A second method, more elegant and most of all serving the purpose to show the use of another virtual method which can be reimplemented in python objects deriving from base, is to override the function that compares two fitness vectors. This function is used by all pagmo algorithms to compare performances of individuals. By default, this function compares the fitness $f_1$ to a fitness $f_2$ and returns true if $f_1$ dominates $f_2$ (which is single objective optimization correspond to minimization). Let us see how...
In [ ]:
class my_problem_max(base):
"""
Analytical function to maximize.
USAGE: my_problem_max()
"""
def __init__(self):
super(my_problem_max, self).__init__(2)
self.set_bounds(-10, 10)
# We provide a list of the best known solutions to the problem
self.best_x = [[1.0, -1.0], ]
# Reimplement the virtual method that defines the objective function
def _objfun_impl(self, x):
f = -(1.0 - x[0]) ** 2 - 100 * (-x[0] ** 2 - x[1]) ** 2 - 1.0
return (f, )
# Reimplement the virtual method that compares fitnesses
def _compare_fitness_impl(self, f1, f2):
return f1[0] > f2[0]
# Add some output to __repr__
def human_readable_extra(self):
return "\n\tMaximization problem"
Additionally in the constructor we provide a list of all known global minima (we will use those later for testing). The list of corresponding objective function values will be then computed and accessible through best_f of the problem's instance.
As before, we use our favorite optimization algorithm:
In [ ]:
from math import sqrt
prob = my_problem_max()
algo = algorithm.de(gen=20)
isl = island(algo, prob, 20)
isl.evolve(10)
isl.join()
print("Best individual:")
print(isl.population.champion)
print("Comparison of the best found fitness with the best known fitness:")
for best_fitness in prob.best_f:
print(best_fitness[0] - isl.population.champion.f[0])
print("L2 distance to the best decision vector:")
for best_decision in prob.best_x:
l2_norm = 0
for n in range(0, len(best_decision)):
l2_norm += (best_decision[n] - isl.population.champion.x[n]) ** 2
l2_norm = sqrt(l2_norm)
print(l2_norm)
Note here that we used the best_f and best_x methods which return the best known fitness and decision vectors. The best_f vector is automatically available as we defined best_x in the problem. With these vectors, we can have an idea of the optimizer performances. The result of this optimization should be in order of $10^{-11}$ for the comparison with the best fitness and $10^{-6}$ for the distance to the best decision vector.
As hinted before, users can also define their own multi-objective problem. In that case we need to overload the the base constructor with third argument stating the desired objective function dimension and return a tuple or a list with more than one element in the objective function implementation (both dimensions must agree).
In [ ]:
class my_mo_problem(base):
"""
A multi-objective problem.
(This is actually a Python implementation of 2-dimensional ZDT-1 problem)
USAGE: my_mo_problem()
"""
def __init__(self, dim=2):
# We call the base constructor as 'dim' dimensional problem, with 0 integer parts and 2 objectives.
super(my_mo_problem,self).__init__(dim, 0, 2)
self.set_bounds(0.0, 1.0)
# Reimplement the virtual method that defines the objective function
def _objfun_impl(self, x):
f0 = x[0]
g = 1.0 + 4.5 * x[1]
f1 = g * (1.0 - sqrt(f0 / g))
return (f0, f1, )
# Add some output to __repr__
def human_readable_extra(self):
return "\n\tMulti-Objective problem"
We instantiate our problem as before, but this time we use one of the multi-objective algorithms available in PaGMO:
In [ ]:
from PyGMO import population
prob = my_mo_problem()
algo = algorithm.sms_emoa(gen=2000) # 2000 generations of SMS-EMOA should solve it
pop = population(prob, 30)
pop = algo.evolve(pop)
Since in the Multi-Objective world the idea of a single 'champion' solution is not very well defined, we plot the Pareto front of the whole population, i.e., the two objectives $f_i^{(1)}$ and $f_i^{(2)}$ of each individual $i \in 1,\ldots,30$.
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
# Fold each objectives into vectors and print the Pareto front
F = np.array([ind.cur_f for ind in pop]).T
plt.scatter(F[0], F[1])
plt.xlabel("$f^{(1)}$")
plt.ylabel("$f^{(2)}$")
plt.show()
NOTE1: This problems of tutorial are implemented in PyGMO under the name PyGMO.problem.py_example and PyGMO.problem.py_example_max
NOTE2: When evolve is called from an island, the process is forked and transferred to another python or ipython instance. As a consequence, when writing your _obj_fun_impl you cannot use stuff like matplotlib to make interactive plots and alike. If you need, during development, to have this kind of support, use the algorithm evolve method (see the optimization of the Multi-Objective problemabove)
NOTE3: If performance is your goal, you should implement your problem in C++, and then expose it into Python.
Now let us see how to code a stochastic optimization problem, that is a problem where the objective function is stochastic.
Just as before we will write a class deriving from PyGMO.problem.base_stochastic and reimplement some of its ‘virtual’ methods. With respect to PyGMO.problem.base, this base class has a field called seed which can be used to control the pseudo random numbers and that is changed by the algorithms that are compatible with stochastic optimization problems.
In [ ]:
from PyGMO.problem import base_stochastic
class my_problem_stochastic(base_stochastic):
"""
Noisy De Jong (sphere) function implemented purely in Python.
USAGE: my_problem_stochastic(dim = 10, seed=0)
* dim problem dimension
* seed initial random seed
"""
def __init__(self, dim=10, seed=0):
#As before, we call the base constructor, but this time with random number generator seed
super(my_problem_stochastic, self).__init__(dim, seed)
self.set_bounds(-5.12, 5.12)
def _objfun_impl(self, x):
from random import random as drng
from random import seed
#We initialize the random number generator using the
#data member seed (in base_stochastic). This will be changed by suitable
#algorithms when a stochastic problem is used. The mod operation avoids overflows
seed(self.seed)
#We write the objfun using the same pseudorandonm sequence
#as long as self.seed is unchanged.
f = 0.0
for i in range(self.dimension):
noise = (2 * drng() - 1) / 10
f = f + (x[i] + noise) * (x[i] + noise)
return (f, )
def human_readable_extra(self):
return "Seed: " + str(self.seed)
Let us now try to solve our stochastic sphere problem using Particle Swarm Optimization algorithm
In [ ]:
from PyGMO import algorithm, island
prob = my_problem_stochastic(dim=10, seed=123456)
algo = algorithm.pso_gen(gen=1)
print("Seed before island initialization: {}".format(prob.seed))
isl = island(algo, prob, 20)
print("Seed after island initialization: {}".format(prob.seed))
for i in range(30):
isl.evolve(1)
print("Best F:{}, seed:{}".format(isl.population.champion.f, isl.problem.seed))
As you can see, the current seed automatically changes during the algorithm run.
In this section we will learn how to code a simple optimization algorithm. We will write the algorithm so that it manage multi-objective, mixed_int, constrained optimization as this will allow us to explain all the basic PyGMO functioning. Clearly our algorithm will not be very good, but a random search is always useful to benchmark against :)
In a nutshell ... we will write a class deriving from PyGMO.algorithm.base and reimplement some of its ‘virtual’ methods, the main one being evolve!!!
In [ ]:
from PyGMO.problem import base
class my_algorithm(base):
"""
Monte-Carlo (random sampling) algorithm implemented purely in Python.
"""
def __init__(self, iter=10):
"""
Constructs a Monte-Carlo (random sampling) algorithm
USAGE: algorithm.my_algorithm(iter=10)
NOTE: At the end of each iteration, the randomly generated
point substitutes the worst individual in the population if better
* iter: number of random samples
"""
#We start calling the base constructor
super(my_algorithm, self).__init__()
#We then define the algorithm 'private' data members
self.__iter = iter
#This is the 'juice' of the algorithm, the method where the actual optimzation is coded.
def evolve(self, pop):
#If the population is empty (i.e. no individuals) nothing happens
if len(pop) == 0:
return pop
#Here we rename some variables, in particular the problem
prob = pop.problem
#Its dimensions (total and continuous)
dim, cont_dim = prob.dimension, prob.dimension - prob.i_dimension
#And the lower/upper bounds for the chromosome
lb, ub = prob.lb, prob.ub
import random
#The algorithm now starts manipulating the population
for _ in range(self.__iter):
# We create a random vector within the bounds.
# First the continuous part:
tmp_cont = [random.uniform(lb[i], ub[i]) for i in range(cont_dim)]
# Then the integer part:
tmp_int = [float(random.randint(lb[i], ub[i])) for i in range(cont_dim, dim)]
# Assemble them into one decision vector:
tmp_x = tmp_cont + tmp_int
# Push back in the population:
pop.push_back(tmp_x)
# Remove the worst individual:
pop.erase(pop.get_worst_idx())
# return the 'evolved' population
return pop
def get_name(self):
return "Monte Carlo Algorithm (Python)"
def human_readable_extra(self):
return "n_iter=" + str(self.__n_iter)
The above code contains a lot of interesting points worth to be discussed. So, we start