Estrategias Evolutivas¶
La librería Pyristic incluye una clase llamada EvolutionStrategy, inspirada en la metaheurística de Estrategias evolutivas (EE), para resolver problemas de minimización. Para trabajar con esta clase se requiere hacer lo siguiente:
Definir:
- La función objetivo $f$.
- La lista de restricciones.
- Lista de límites inferiores y superiores.
- Configuración de operadores de la metaheurística (opcional).
Crear una clase que herede de
EvolutionStrategy.Sobreescribir las siguientes funciones de la clase
EvolutionStrategy:- initialize_step_weights (opcional)
- initialize_population (opcional)
- fixer (opcional)
- mutation_operator (opcional)
- crossover_operator (opcional)
- adaptive_crossover (opcional)
- adaptive_mutation (opcional)
- survivor_selection (opcional)
A continuación se mostrará cómo resolver dos problemas de optimización continua usando la clase EvolutionStrategy. El primer paso es importar la librería y los módulos que se van a utilizar.
Librerías externas¶
from pprint import pprint
import math
import numpy as np
import copy
from IPython.display import Image
from IPython.core.display import HTML
Componentes de pyristic¶
La estructura que está organizada la librería es:
- Las metaheurísticas están ubicadas en
heuristic. - Las funciones de prueba están ubicadas en
utils.test_function. - Las clases auxiliares para mantener la información de los operadores que serán empleados para alguna de las metaheurísticas basadas en los paradigmas del cómputo evolutivo están ubicadas en
utils.helpers. - Las metaheurísticas basadas en los paradigmas del cómputo evolutivo dependen de un conjunto de operadores (selección, mutación y cruza). Estos operadores están ubicados en
utils.operators.
Para demostrar el uso de nuestra metaheurística basada en algoritmos geneticos tenemos que importar la clase llamada Genetic que se encuentra en heuristic.GeneticAlgorithm_search.
from pyristic.heuristic.EvolutionStrategy_search import EvolutionStrategy
from pyristic.utils.operators import selection, mutation, crossover
from pyristic.utils.test_function import beale_, ackley_
from pyristic.utils.helpers import EvolutionStrategyConfig,get_stats
Función de Beale¶
\begin{equation} \label{eq:BF} \begin{array}{rll} \text{minimizar:} & f(x_1, x_2) = (1.5 - x_1 + x_1x_2)^2 + (2.25 - x_1 + x_1x_2^2)^2 + (2.625 - x_1 + x_1x_2^3)^2 & \\ \text{Tal que: } & -4.5 \leq x_1,x_2 \leq 4.5 & \end{array} \end{equation}El mínimo global se encuentra en $x^* = (3, 0.5)$ y $f(x^*) = 0$.
Image(filename="include/beale.png", width=500, height=300)
La librería Pyristic tiene implementados algunos problemas de prueba en utils.helpers.test_function, entre ellos la función de Beale. Los problemas de prueba están definidos como diccionarios con las siguientes llaves:
function.Función objetivo.constraints.Lista con las funciones que validan las diferentes restricciones del problema. Todas las funciones deben regresar un valor booleano.bounds.Límites para cada una de las variables del problema. En el caso de que todas las variables del problema compartan el mismo intervalo de búsqueda, se utiliza una lista con dos valores numéricos. El primero indica el límite inferior y el segundo el límite superior. En caso contrario, se tendrá una lista con dos listas internas. La primera lista almacena los límites inferiores y la segunda los límites superiores.decision_variables.Número de variables de decisión.
beale_
{'function': CPUDispatcher(<function beale_function at 0x7f018684ac20>),
'constraints': [<function pyristic.utils.test_function.constraint1_beale(X: numpy.ndarray) -> bool>],
'bounds': [-4.5, 4.5],
'decision_variables': 2}
Función de aptitud¶
Los algoritmos evolutivos, se define una función llamada aptitud que describe a cada individuo que tan bueno es en el problema a resolver. Los valores más grandes de aptitud corresponden a las mejores soluciones. Para el problema de Beale, la función de aptitud es: $F_a(i) = \frac{1}{f(x^{(i)}_1, x^{(i)}_2) + 1}$
def aptitudeBeale(X: np.array) -> np.array:
return 1/ ( beale_['function'](X)+1)
Declaración de EvolutionStrategy¶
La metaheurística de EE implementada en la librería de Pyristic se puede utilizar de las siguientes maneras:
- Crear una clase que herede de la clase
EvolutionStrategyy sobreescribir alguno de los métodos. - Declarar un objeto de tipo EvolutionStrategyConfig, ubicado en utils.helpers, y utilizarlo para inicializar un objeto de tipo
EvolutionStrategy. - Realizar una combinación de las dos anteriores.
Es importante resaltar que se puede hacer uso de la metaheurística sin modificar los operadores que tiene por defecto.
Ejecución de la metaheurística¶
Como mencionamos antes, una forma de utilizar la metaheurística es hacer una instancia de la clase EvolutionStrategy dejando su configuración por defecto.
Los argumentos que se deben indicar al momento de inicializar son:
- función objetivo
- restricciones del problema
- límites inferior y superior (por cada variable de decisión)
- número de variables que tiene el problema.
Beale = EvolutionStrategy(
function= aptitudeBeale,\
decision_variables=beale_['decision_variables'],\
constraints= beale_['constraints'],\
bounds= beale_['bounds'])
Recordemos que beale_ es un diccionario con la información requerida por el constructor de la clase EvolutionStrategy.
($\mu + \lambda$) - EE¶
Por defecto, la metaheurística trabaja de la siguiente forma:
- Para la selección de sobrevivientes utiliza un esquema $(\mu + \lambda)$.
- El operador de cruza para las variables de decisión es recombinación discreta (operador llamado
discrete). - El operador de cruza para los tamaños de paso es recombinación intermedia (operador llamado
intermediate). - Cada individuo tiene un tamaño de paso por variable. (operador llamado
mult_sigma_adaptive_mutator).
Una vez creado el objeto, se manda a llamar a su método optimize el cual recibe los siguientes parámetros:
- generations = $300$
- population_size = $80$
- offspring_size = $160$
Beale.optimize(300,80,160,verbose=True)
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 300/300 [00:02<00:00, 121.70it/s]
print(Beale)
Evolution Strategy search: F_a(X) = 1.0 X = [2.99999999 0.5 ] Constraints: x1: -4.5 <= 3.00 <= 4.5 x2: -4.5 <= 0.50 <= 4.5
Para realizar un estudio estadístico del comportamiento de la metaheurística, la librería Pyristic cuenta con una función llamada get_stats. Esta función se encuentra en utils.helpers y recibe como parámetros:
- Objeto que realiza la ejecución de la metaheurística.
- Número de veces que se quiere ejecutar la metaheurística.
- Tupla con los argumentos que recibe la función
optimize. - Argumentos adicionales a la búsqueda (opcional).
La función get_stats considera la solución devuelta por la metaheurística en cada ejecución y retorna un diccionario con la mejor y peor solución encontrada y la media y desviación estándar del valor de la función objetivo.
args = (200, 80, 160, 0.0001, False)
statistics = get_stats(Beale, 21, args, transformer=beale_['function'])
pprint(statistics)
{'Best solution': {'f': 5.232655460094757e-18, 'x': array([3. , 0.5])},
'Mean': 7.500236923765416e-13,
'Median': 6.295559358095275e-17,
'Standard deviation': 3.1450540830911966e-12,
'Worst solution': {'f': 1.47876718588774e-11,
'x': array([3.00000693, 0.50000227])},
'averageTime': 0.19353040059407553,
'objectiveFunction': [6.974537963159639e-17,
1.5829469345976208e-17,
5.232655460094757e-18,
3.55142232876361e-17,
3.659867353206695e-14,
6.295559358095275e-17,
9.368025066127844e-17,
2.3280831331374633e-17,
9.248440100693195e-13,
5.228335842212396e-17,
9.484958935479025e-18,
4.58002399942093e-17,
7.223764853904115e-17,
2.4479100052993973e-17,
1.411336805961557e-17,
9.124874382780466e-17,
1.47876718588774e-11,
4.83541251795498e-16,
3.808890680308518e-17,
1.58427128576792e-16,
8.705432027964274e-17]}
(1+1) - EE¶
A continuación vamos a implementar la EE más simple conocida como $(1+1)-$EE para resolver la función de Beale. $(1+1)-$EE trabaja de la siguiente forma:
- Utiliza un único valor de $\sigma$ para mutar todas las variables de decisión.
- El valor de $\sigma$ se autoadapta utilizando la regla del éxito del 1/5. \begin{equation} \sigma = \begin{cases} \sigma / c & \text{Si } p_s > 1/5, \\ \sigma \cdot c & \text{Si } p_s < 1/5, \\ \sigma & \text{Si } p_s = 1/5 \\ \end{cases} \end{equation}
- Se utiliza un solo individuo, éste es mutado para generar un hijo.
- El individuo más apto es el que pasará a la siguiente generación.
Esta versión no considera un operador de cruza. Por este motivo, vamos a sobreescribir crossover_operator y adaptive_crossover para que retornen el individuo original sin cambio alguno. Las funciones que debemos sobreescribir son:
- crossover_operator
- adaptive_crossover
- adaptive_mutation
- mutation_operator
- initializing_step_weights
class ESBasic(EvolutionStrategy):
def __init__(self, function ,\
decision_variables:int,\
constraints: list,\
bounds: np.ndarray):
super().__init__(function, decision_variables, constraints, bounds)
self.successful = 0
def crossover_operator(self, parent_ind1,\
parent_ind2,\
**kargs):
return self.logger['parent_population_x']
def adaptive_crossover(self, parent_ind1,\
parent_ind2,\
**kargs):
return self.logger['parent_population_sigma']
def adaptive_mutation(self, **kargs):
if self.logger['current_iter'] % kargs['k'] != 0:
return self.logger['offspring_population_sigma']
ps = self.successful / kargs['k']
self.successful = 0
if( ps > 1/5):
return self.logger['offspring_population_sigma']/kargs['c']
elif( ps < 1/5):
return self.logger['offspring_population_sigma']*kargs['c']
else:
return self.logger['offspring_population_sigma']
def mutation_operator(self, **kargs):
X_mutated = copy.deepcopy(self.logger['offspring_population_x'])
X_mutated += self.logger['offspring_population_sigma'] * np.random.normal(0,1, size=X_mutated.shape)
if(self.f(self.logger['offspring_population_x'][0]) > self.f(X_mutated[0])):
self.successful+=1
return X_mutated
def initializing_step_weights(self, eps_sigma:int, **kwargs) -> np.ndarray:
return np.random.uniform(0,1, size=(self.logger['parent_population_size'], 1))
Beale_basic = ESBasic(
function= aptitudeBeale,\
decision_variables=beale_['decision_variables'],\
constraints= beale_['constraints'],\
bounds= beale_['bounds'])
Como se puede observar, estamos haciendo uso del argumento llamado kargs, el cual es un diccionario que contiene información adicional. En este caso kargs contiene dos constantes: la primera es el parámetro $k$ que indica cada cuantas generaciones se va a actualizar el valor de $\sigma$ y la segunda es el valor del parámetro $c$ que se utiliza para actualizar el valor de $\sigma$. El diccionario tiene que ser incluido siempre con ** antes del nombre de la variable que lo almacena.
additional_arguments = {'k':20,'c':0.83}
Beale_basic.optimize(generations=200,\
population_size=1,\
offspring_size=1,\
**additional_arguments)
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 200/200 [00:00<00:00, 1739.90it/s]
print(Beale_basic)
Evolution Strategy search: F_a(X) = 0.9969506868038376 X = [2.98335859 0.50729749] Constraints: x1: -4.5 <= 2.98 <= 4.5 x2: -4.5 <= 0.51 <= 4.5
args = (500, 1, 1, 0.001,False)
statistics = get_stats(Beale_basic, 30, args, additional_arguments, aptitudeBeale)
pprint(statistics)
{'Best solution': {'f': 0.09378523089864498,
'x': array([ 0.03287226, -3.80633412])},
'Mean': 0.7862535058145801,
'Median': 0.9967553646899534,
'Standard deviation': 0.3029875418531697,
'Worst solution': {'f': 0.9999381625668178,
'x': array([2.9834752 , 0.49499717])},
'averageTime': 0.043139831225077314,
'objectiveFunction': [0.9990169200052803,
0.09378523089864498,
0.5672019868347583,
0.9981812538791579,
0.9979996162555546,
0.9977846479786782,
0.9937017740501066,
0.5163703703126666,
0.5472588628918187,
0.49558365225972467,
0.9959051556100147,
0.9972898297178926,
0.9987279969858004,
0.5645729629595143,
0.09596677560117627,
0.996932108411761,
0.9993865654887663,
0.9991981127606973,
0.9983601706985583,
0.9970646756951872,
0.5675010831608579,
0.9982544529717703,
0.9999381625668178,
0.9965786209681458,
0.5395095891789323,
0.9892498155559739,
0.9985622337077901,
0.09551406458177113,
0.5544542372737761,
0.9977542451758096]}
Función de Ackley¶
\begin{equation} \min f(\vec{x}) = -20\exp \left( -0.2 \sqrt{\frac{1}{n} \sum_{i=1}^n x_i^2} \right) - exp \left( \frac{1}{n} \sum_{i=1}^n \cos (2\pi x_i) \right) + 20 + e \end{equation}El mínimo global está en $x^* = 0 $, $f(\vec{x}) = 0$ y su dominio es $|x_{i}| < 30$.
Image(filename="include/ackley.jpg", width=500, height=300)
Al igual que la función de Beale, la librería Pyristic tiene implementada la función de Ackley en utils.helpers.test_function.
ackley_
{'function': CPUDispatcher(<function ackley_function at 0x7f018561c950>),
'constraints': [<function pyristic.utils.test_function.constraint1_ackley(X: numpy.ndarray) -> bool>],
'bounds': [-30.0, 30.0],
'decision_variables': 10}
Función de aptitud¶
Los algoritmos evolutivos, se define una función llamada aptitud que describe a cada individuo que tan bueno es en el problema a resolver. Los valores más grandes de aptitud corresponden a las mejores soluciones. Para el problema de Ackley, la función de aptitud es: $F_a(\vec{X}) = -f(\vec{X})$
def aptitudeAckley(x):
return -1* ackley_['function'](x)
Este problema no se encuentra restringido a un número de variables de decisión. Para modificar el número de variables de decisión hacemos:
ackley_['decision_variables'] = 5
($\mu+\lambda$)-EE¶
Emplearemos la metaheurística de EE, con la configuración por defecto, para resolver la función de Ackley con 10 variables de decisión.
Ackley = EvolutionStrategy(
function= aptitudeAckley,\
decision_variables=ackley_['decision_variables'],\
constraints= ackley_['constraints'],\
bounds= ackley_['bounds'])
'''
Número de generaciones: 250
Tamaño de población de padres: 100
Tamaño de población de hijos: 200
'''
Ackley.optimize(250,100,200)
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 250/250 [00:00<00:00, 289.85it/s]
print(Ackley)
Evolution Strategy search: F_a(X) = -0.28845250921960686 X = [-0.05837 -0.10053953 -0.00708623 -0.03774568 -0.03316804 -0.03318999 0.02397432 -0.00125536 -0.00243173 -0.0546666 ] Constraints: x1: -30 <= -0.06 <= 30 x2: -30 <= -0.10 <= 30 x3: -30 <= -0.01 <= 30 x4: -30 <= -0.04 <= 30 x5: -30 <= -0.03 <= 30 x6: -30 <= -0.03 <= 30 x7: -30 <= 0.02 <= 30 x8: -30 <= -0.00 <= 30 x9: -30 <= -0.00 <= 30 x10: -30 <= -0.05 <= 30
args = (250, 100, 200, 0.001 ,False)
statistics = get_stats(Ackley, 30, args, transformer=ackley_['function'])
pprint(statistics)
{'Best solution': {'f': 0.010249799629580725,
'x': array([-0.00086015, -0.00161265, -0.00215718, 0.00065416, -0.00228198,
-0.00559011, 0.00025284, -0.00344554, 0.00059995, -0.0020872 ])},
'Mean': 0.1644989262654607,
'Median': 0.07446471752602557,
'Standard deviation': 0.28754471244926094,
'Worst solution': {'f': 1.628041385916028,
'x': array([-0.26780161, -0.21837777, -0.1178678 , -0.09132986, -0.26427773,
0.10544494, -0.13091937, 0.10397274, -0.12819718, 0.03343717])},
'averageTime': 0.8125235954920451,
'objectiveFunction': [0.034206526420007766,
1.628041385916028,
0.12596106222010706,
0.22109914950204912,
0.046531228197921326,
0.04410535059719267,
0.09213135864640565,
0.2619662918114547,
0.015350036094688324,
0.0419903892174287,
0.031469512589055615,
0.11911705593879462,
0.27359047848434104,
0.2391134627754572,
0.06432248377842997,
0.21975810167185772,
0.07006921608206129,
0.0369945594033827,
0.031840050386872765,
0.2748363430488152,
0.010249799629580725,
0.03167359762541766,
0.09315407594888869,
0.0777059003499656,
0.07122353470208553,
0.33507869438961224,
0.020629669655939598,
0.036820161759220316,
0.1742926822490216,
0.21164562887173988]}
($\mu,\lambda$)-EE¶
Ahora vamos a utilizar una configuración del tipo EvolutionStrategyConfig para definir los métodos que serán empleados en nuestro objeto del tipo EvolutionStrategy.
configuration_ackley = (EvolutionStrategyConfig()
.survivor_selection(selection.replacement_selector())
.adaptive_mutation(
mutation.single_sigma_adaptive_mutator(
ackley_['decision_variables'])
)
)
En este caso estamos modificando el esquema de selección de sobrevivientes y emplearemos un único tamaño de paso para cada individuo, donde, el argumento que tiene single_sigma_adaptive_mutator es el número de variables de decisión del problema.
print(configuration_ackley)
-------------------------------- Configuration -------------------------------- Survivor selection: Replacement population Adaptive mutation: Single Sigma --------------------------------
solver_ackley_custom = EvolutionStrategy(
function= aptitudeAckley,\
decision_variables=ackley_['decision_variables'],\
constraints= ackley_['constraints'],\
bounds= ackley_['bounds'],config=configuration_ackley)
solver_ackley_custom.optimize(250,100,200)
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 250/250 [00:01<00:00, 179.18it/s]
print(solver_ackley_custom)
Evolution Strategy search: F_a(X) = -21.55196699129465 X = [27.23771294 27.23771294 27.23771294 27.23771294 27.23771294 27.23771294 27.23771294 27.23771294 27.23771294 27.23771294] Constraints: x1: -30 <= 27.24 <= 30 x2: -30 <= 27.24 <= 30 x3: -30 <= 27.24 <= 30 x4: -30 <= 27.24 <= 30 x5: -30 <= 27.24 <= 30 x6: -30 <= 27.24 <= 30 x7: -30 <= 27.24 <= 30 x8: -30 <= 27.24 <= 30 x9: -30 <= 27.24 <= 30 x10: -30 <= 27.24 <= 30
args = (250, 100, 200,0.001 ,False)
statistics = get_stats(solver_ackley_custom, 30, args, transformer=ackley_['function'])
pprint(statistics)
{'Best solution': {'f': 3.914908634326306,
'x': array([ 0.22236554, -0.78155697, -0.90025871, 0.27553164, 1.10221801,
-0.91454935, 0.07856534, 0.78949362, -0.93589365, 0.36879387])},
'Mean': 17.869116590215647,
'Median': 19.950424956466673,
'Standard deviation': 5.142625044782836,
'Worst solution': {'f': 22.13188156750359,
'x': array([-29.6417392, -29.6417392, -29.6417392, -29.6417392, -29.6417392,
-29.6417392, -29.6417392, -29.6417392, -29.6417392, -29.6417392])},
'averageTime': 1.0821852445602418,
'objectiveFunction': [10.214464573662518,
21.68444010317522,
3.914908634326306,
19.950424956466673,
19.950424956466673,
21.773356498999437,
22.13188156750359,
19.950424956466673,
20.745762410131494,
7.161152173042002,
19.950424956466673,
19.647852175141484,
19.950424956466673,
19.950424956466673,
20.08725543332284,
19.950424956466673,
20.22782779904124,
19.950424956466673,
19.950424956466673,
19.950424956466673,
19.950424956466673,
5.0998160615647485,
19.950424956466673,
19.950424956466673,
20.368000242171377,
20.506531700594646,
19.950424956466673,
8.02794632911078,
14.300142756300238,
20.876209857848135]}
(1+1)-EE¶
A continación usaremos la EE más simple para intentar resolver la función de Ackley con 10 variables de decisión.
Ackley_basic = ESBasic(
function= aptitudeAckley,\
decision_variables=ackley_['decision_variables'],\
constraints= ackley_['constraints'],\
bounds= ackley_['bounds'])
Ackley_basic.optimize(10000,1,1,**additional_arguments )
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10000/10000 [00:01<00:00, 7829.57it/s]
print(Ackley_basic)
Evolution Strategy search: F_a(X) = -19.713742965447526 X = [-30. -6.06618645 28.00763708 -4.01995382 9.98348065 10.95156463 -12.98341022 30. -0.06114509 -19.18521633] Constraints: x1: -30 <= -30.00 <= 30 x2: -30 <= -6.07 <= 30 x3: -30 <= 28.01 <= 30 x4: -30 <= -4.02 <= 30 x5: -30 <= 9.98 <= 30 x6: -30 <= 10.95 <= 30 x7: -30 <= -12.98 <= 30 x8: -30 <= 30.00 <= 30 x9: -30 <= -0.06 <= 30 x10: -30 <= -19.19 <= 30
args = (10000, 1, 1, 0.001,False)
statistics = get_stats(Ackley_basic, 30, args, additional_arguments, transformer=ackley_['function'])
pprint(statistics)
{'Best solution': {'f': 10.310536059391417,
'x': array([ 0.03343847, -0.78729196, -0.13713295, -0.07919603, -9.08132812,
2.9751608 , 1.23007554, -0.91011728, 0.84557144, -2.29842616])},
'Mean': 19.29082310550591,
'Median': 19.663565578248548,
'Standard deviation': 1.7281086637591332,
'Worst solution': {'f': 20.26451261448206,
'x': array([ 22.10012059, -14.99663106, 20.77574185, 22.88045243,
-19.97996349, 11.1043873 , 27.08426041, -27.00899638,
16.85222655, 27.07655944])},
'averageTime': 0.9502643664677938,
'objectiveFunction': [19.766594434620053,
19.69738945159391,
20.087764430839382,
19.532432501793185,
19.2855996760406,
19.476403133473898,
19.81576868565518,
19.34851356423621,
19.35177435923925,
19.76326438683161,
19.64038637966428,
20.21808540377747,
20.02292698192497,
18.808488740334017,
20.05904845457353,
19.123067934794904,
19.44640999478022,
19.44180567590634,
19.93611334585588,
19.70559439438051,
19.06727228643281,
20.125716751767744,
17.976236933257045,
19.576789328814215,
19.551319019341047,
19.701331874439493,
19.936801590103148,
10.310536059391417,
19.686744776832814,
20.26451261448206]}