Optimización con algoritmo genético y Nelder-Mead
Joaquín Amat Rodrigo
Febrero, 2019
Más sobre ciencia de datos en cienciadedatos.net
Versión PDF: Github
Introducción
Los algoritmos genéticos son métodos de optimización heurística que, entre otras aplicaciones, pueden emplearse para encontrar el valor o valores que consiguen maximizar o minimizar una función. Su funcionamiento está inspirado en la teoría evolutiva de selección natural propuesta por Darwin y Alfred Russel: los individuos de una población se reproducen generando nuevos descendientes, cuyas características, son combinación de las características de los progenitores (más ciertas mutaciones). De todos ellos, únicamente los mejores individuos sobreviven y pueden reproducirse de nuevo, transmitiendo así sus características a las siguientes generaciones.
Los algoritmos genéticos son solo una de las muchas estrategias de optimización que existen, y no tiene por qué ser la más adecuada en todos los escenarios. Por ejemplo, si el problema en cuestión puede optimizarse de forma analítica, suele ser más adecuado resolverlo de esta forma.
La implementación de algoritmo genético que se muestra en este documento pretende ser lo más explicativa posible aunque para ello no sea la más eficiente.
El código de las funciones desarrolladas a lo largo del documento puede descargarse en el siguiente link.
Algoritmo
Aunque existen variaciones, algunas de las cuales se describen a lo largo de este documento, en términos generales, la estructura de un algoritmo genético para optimizar (maximizar o minimizar) una función con una o múltiples variables sigue los siguientes pasos:
Crear una población inicial aleatoria de \(P\) individuos. En este caso, cada individuo representa una combinación de valores de las variables.
Calcular la fortaleza (fitness) de cada individuo de la población. El fitness está relacionado con el valor de la función para cada individuo. Si se quiere maximizar, cuanto mayor sea el valor de la función para el individuo, mayor su fitness. En el caso de minimización, ocurre lo contrario.
Crear una nueva población vacía y repetir los siguientes pasos hasta que se hayan creado \(P\) nuevos individuos.
3.1 Seleccionar dos individuos de la población existente, donde la probabilidad de selección es proporcional al fitness de los individuos.
3.2 Cruzar los dos individuos seleccionados para generar un nuevo descendiente (crossover).
3.3 Aplicar un proceso de mutación aleatorio sobre el nuevo individuo.
3.4 Añadir el nuevo individuo a la nueva población.
Reemplazar la antigua población por la nueva.
Si no se cumple un criterio de parada, volver al paso 2.
En los siguientes apartados se describe cada una de las etapas del proceso para, finalmente, combinarlas todas en una única función.
Población inicial
En el contexto de algoritmos genéticos, el término individuo hace referencia a cada una de las posibles soluciones del problema que se quiere resolver. En el caso de maximización o minimización de una función, cada individuo representa una posible combinación de valores de las variables. Para representar dichas combinaciones, se pueden emplear vectores, cuya longitud es igual al número total de variables, y cada posición toma un valor numérico. Por ejemplo, supóngase que la función objetivo \(J(x,y,z)\) depende de las variables \(x, y, z\). El individuo \(3, 9.5, -0.5\), equivale a la combinación de valores \(x = 3, y = 9.5, z = -0.5\).
El primer paso del algoritmo genético consiste en crear una población inicial aleatoria de individuos. La siguiente función crea una matriz en la que, cada fila, está formada por una combinación de valores numéricos aleatorios. Además, el valor para cada variable puede estar acotado dentro de un rango. Esta acotación resulta útil para agilizar el proceso de optimización, pero requiere disponer de información que permita acotar el intervalo de valores dentro del cual se encuentra la solución óptima.
crear_poblacion <- function(n_poblacion, n_variables, limite_inf = NULL,
limite_sup = NULL, verbose = TRUE) {
# Esta función crea una matriz en la que, cada fila, está formada por una
# combinación de valores numéricos aleatorios. El rango de posibles valores
# para cada variable puede estar acotado.
#
# ARGUMENTOS
# ============================================================================
# n_poblacion: número total de individuos de la población.
# n_variables: longitud de los individuos.
# limite_inf: vector con el límite inferior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quiere acotar.
# limite_sup: vector con el límite superior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quieren acotar.
# verbose: mostrar información del proceso por pantalla.
#
# RETORNO
# ============================================================================
# Una matriz de tamaño n_poblacion x n_variables que representa una población.
# COMPROBACIONES
# ----------------------------------------------------------------------------
if (!is.null(limite_inf) & (length(limite_inf) != n_variables)) {
stop(paste(
"limite_inf debe tener un valor por cada variable.",
"Si para alguna variable no se quiere límite, emplear NA.",
"Ejemplo: lim_sup = c(10, NA, 10)"
))
} else if (!is.null(limite_sup) & length(limite_sup) != n_variables) {
stop(paste(
"limite_sup debe tener un valor por cada variable.",
"Si para alguna variable no se quiere límite, emplear NA.",
"Ejemplo: lim_sup = c(10, NA, 10)"
))
} else if (is.null(limite_sup) | is.null(limite_inf)) {
warning(paste(
"Es altamente recomendable indicar los límites dentro de los",
"cuales debe buscarse la solución de cada variable.",
"Por defecto se emplea [-10^3, 10^3]."
))
} else if (any(any(is.na(limite_sup)), any(is.na(limite_inf)))) {
warning(paste(
"Los límites empleados por defecto cuando no se han definido son:",
" [-10^3, 10^3]."
))
cat("\n")
}
# Si no se especifica limite_inf, el valor mínimo que pueden tomar las variables
# es -10^3.
if (is.null(limite_inf)) {
limite_inf <- rep(x = -10^3, times = n_variables)
}
# Si no se especifica limite_sup, el valor máximo que pueden tomar las variables
# es 10^3.
if (is.null(limite_sup)) {
limite_sup <- rep(x = 10^3, times = n_variables)
}
# Si los límites no son nulos, se reemplazan aquellas posiciones NA por el valor
# por defecto -10^3 y 10^3
if (!is.null(limite_inf)) {
limite_inf[is.na(limite_inf)] <- -10^3
}
if (!is.null(limite_sup)) {
limite_sup[is.na(limite_sup)] <- 10^3
}
# CREAR POBLACIÓN
# ----------------------------------------------------------------------------
# Matriz donde almacenar los individuos generados.
poblacion <- matrix(data = NA, nrow = n_poblacion, ncol = n_variables)
# Bucle para crear cada individuo.
for (i in 1:n_poblacion) {
# Se crea un vector de NA que representa el individuo.
individuo <- rep(NA, times = n_variables)
for (j in 1:n_variables) {
# Para cada posición, se genera un valor aleatorio dentro del rango permitido
# para cada variable.
individuo[j] <- runif(n = 1, min = limite_inf[j], max = limite_sup[j])
}
# Se añade el nuevo individuo a la población.
poblacion[i, ] <- individuo
}
# INFORMACIÓN ALMACENADA EN LOS ATRIBUTOS
# ----------------------------------------------------------------------------
attr(poblacion, 'fecha_creacion') <- Sys.time()
attr(poblacion, 'numero_individuos') <- n_poblacion
attr(poblacion, "class") <- c("matrix", "poblacion")
if (verbose) {
cat("Población inicial creada", "\n")
cat("------------------------", "\n")
cat("Fecha creación:", as.character(Sys.time()), "\n")
cat("Número de individuos =", n_poblacion, "\n")
cat("Límites inferiores de cada variable =", paste(limite_inf, collapse = ", "), "\n")
cat("Límites superiores de cada variable =", paste(limite_sup, collapse = ", "), "\n")
cat("\n")
}
return(poblacion)
}Ejemplo
Se crea una población de 10 individuos de longitud 2, con los valores de la primera variable acotados entre [-100, +100] y la segunda con únicamente el límite inferior [-20, NA].
poblacion <- crear_poblacion(
n_poblacion = 10,
n_variables = 2,
limite_inf = c(-100, -20),
limite_sup = c(+100, NA),
verbose = TRUE
)##
## Población inicial creada
## ------------------------
## Fecha creación: 2020-11-12 12:21:44
## Número de individuos = 10
## Límites inferiores de cada variable = -100, -20
## Límites superiores de cada variable = 100, 1000## [,1] [,2]
## [1,] -62.28138 335.615078
## [2,] -98.78874 240.589471
## [3,] -16.62883 5.461185
## [4,] 77.30608 -18.495540
## [5,] 43.67225 709.264950
## [6,] 82.62979 549.810212
## [7,] -62.29848 427.149130
## [8,] 27.80573 933.015135
## [9,] 54.55204 522.964812
## [10,] -60.75474 990.104857
## attr(,"fecha_creacion")
## [1] "2020-11-12 12:21:44 UTC"
## attr(,"numero_individuos")
## [1] 10
## attr(,"class")
## [1] "matrix" "poblacion"Fitness de un individuo
Cada individuo de la población debe ser evaluado para cuantificar cómo de bueno es como solución al problema, a esta cuantificación se le llama (fitness). Dependiendo de si se trata de un problema de maximización o minimización, la relación del fitness con la función objetivo \(f\) puede ser:
Maximización: el individuo tiene mayor fitness cuanto mayor es el valor de la función objetivo \(f(individuo)\).
Minimización: el individuo tiene mayor fitness cuanto menor es el valor de la función objetivo \(f(individuo)\), o lo que es lo mismo, cuanto mayor es el valor de la función objetivo, menor el fitness. Tal y como se describe más adelante, el algoritmo genético selecciona los individuos de mayor fitness, por lo que, para problemas de minimización, el fitness puede calcularse como \(-f(individuo)\) o también \(\frac{1}{1+f(individuo)}\).
calcular_fitness_individuo <- function(individuo, funcion_objetivo, optimizacion,
verbose = TRUE, ...) {
# Esta función devuelve el fitness de cada individuo de una población.
#
# ARGUMENTOS
# ============================================================================
# individuo: vector con los valores de cada variable. El orden de los
# valores debe coincidir con el de los argumentos de la
# función.
# funcion_objetivo: nombre de la función que se desea optimizar. Debe de haber
# sido definida previamente.
# optimizacion: "maximizar" o "minimizar". Dependiendo de esto, la relación
# del fitness es directamente o indirectamente proporcional
# al valor de la función.
# verbose: mostrar información del proceso por pantalla.
#
# RETORNO
# ============================================================================
# fitness del individuo.
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
if (length(individuo) != length(names(formals(funcion_objetivo)))) {
stop(paste("Los individuos deben tener tantos valores como argumentos tiene",
"la función objetivo."))
}
# CÁLCULO FITNESS
# ----------------------------------------------------------------------------
if (optimizacion == "maximizar") {
fitness <- do.call(funcion_objetivo, args = as.list(individuo))
} else if (optimizacion == "minimizar") {
fitness <- -(do.call(funcion_objetivo, args = as.list(individuo)))
} else {
stop("El argumento optimización debe ser maximizar o minimizar.")
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("El individuo ha sido evaluado", "\n")
cat("-----------------------------", "\n")
cat("Optimización =", optimizacion, "\n")
cat("Individuo =", paste(individuo, collapse = " "), "\n")
cat("Fitness =", fitness, "\n")
cat("\n")
}
return(fitness)
}Ejemplo
Se calcula el fitness del individuo \((x_1= 10, x_2 = 10)\) para los casos de maximización y minimización de la función \(f(x_1,x_2) = x_1 + x_2\).
# Función objetivo a optimizar.
funcion <- function(x1, x2) {
return(x1 + x2)
}
fitness <- calcular_fitness_individuo(
individuo = c(10, 10),
funcion_objetivo = funcion,
optimizacion = "maximizar",
verbose = TRUE
)## El individuo ha sido evaluado
## -----------------------------
## Optimización = maximizar
## Individuo = 10 10
## Fitness = 20fitness <- calcular_fitness_individuo(
individuo = c(10, 10),
funcion_objetivo = funcion,
optimizacion = "minimizar",
verbose = TRUE
)## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = 10 10
## Fitness = -20Fitness de todos los individuos de una población
Esta función recibe como argumentos una población de individuos, una función objetivo y el tipo de optimización, y devuelve el fitness de todos los individuos.
calcular_fitness_poblacion <- function(poblacion, funcion_objetivo, optimizacion,
verbose = TRUE, ...) {
# Esta función devuelve el fitness de cada individuo de una población.
#
# ARGUMENTOS
# ============================================================================
# poblacion: matriz que representa la población de individuos.
# funcion_objetivo: nombre de la función que se desea optimizar. Debe de haber
# sido definida previamente.
# optimizacion: "maximizar" o "minimizar". Dependiendo de esto, la relación
# del fitness es directamente o indirectamente proporcional
# al valor de la función.
# verbose: mostrar información del proceso por pantalla.
#
# RETORNO
# ============================================================================
# Vector con el fitness de todos los individuos de la población. El orden de
# los valores se corresponde con el orden de las filas de la matriz población.
# CÁLCULO DEL FITNESS DE CADA INDIVIDUO DE LA POBLACIÓN
# ----------------------------------------------------------------------------
# Vector donde almacenar el fitness de cada individuo.
fitness_poblacion <- rep(NA, times = nrow(poblacion))
for (i in 1:nrow(poblacion)) {
individuo <- poblacion[i, ]
fitness_individuo <- calcular_fitness_individuo(
individuo = individuo,
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
verbose = verbose
)
fitness_poblacion[i] <- fitness_individuo
}
# MEJOR INDIVIDUO DE LA POBLACIÓN
# ----------------------------------------------------------------------------
# Se identifica el mejor individuo de toda la población, el de mayor
# fitness.
indice_mejor_individuo <- which.max(fitness_poblacion)
# Se identifica el valor de la función objetivo para el mejor individuo.
if (optimizacion == "maximizar") {
valor_funcion <- fitness_poblacion[indice_mejor_individuo]
} else {
valor_funcion <- -1*fitness_poblacion[indice_mejor_individuo]
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("------------------", "\n")
cat("Población evaluada", "\n")
cat("------------------", "\n")
cat("Optimización =", optimizacion, "\n")
cat("Mejor fitness encontrado =", fitness_poblacion[indice_mejor_individuo], "\n")
cat("Mejor solución encontrada =",
paste(poblacion[indice_mejor_individuo,], collapse = " "), "\n")
cat("Valor función objetivo =", valor_funcion, "\n")
cat("\n")
}
return(fitness_poblacion)
}Ejemplo
Se calcula el fitness de todos los individuos de una población formada por 5 individuos.
# Función objetivo a optimizar.
funcion <- function(x1, x2) {
return(x1 + x2)
}
# Población simulada.
poblacion <- crear_poblacion(
n_poblacion = 5,
n_variables = 2,
limite_inf = c(-10, -10),
limite_sup = c(+10, +10),
verbose = TRUE
)## Población inicial creada
## ------------------------
## Fecha creación: 2020-11-12 12:21:44
## Número de individuos = 5
## Límites inferiores de cada variable = -10, -10
## Límites superiores de cada variable = 10, 10# Cálculo del fitness de todos los individuos.
fitness_poblacion <- calcular_fitness_poblacion(
poblacion = poblacion,
funcion_objetivo = funcion,
optimizacion = "minimizar",
verbose = TRUE
)## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = -7.06204414367676 -0.714093213900924
## Fitness = 7.776137
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = -0.186551273800433 6.2824933975935
## Fitness = -6.095942
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = 0.552894314751029 0.834178593941033
## Fitness = -1.387073
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = 7.59759249631315 -1.4382198266685
## Fitness = -6.159373
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = -6.80632072035223 0.723222000524402
## Fitness = 6.083099
##
## ------------------
## Población evaluada
## ------------------
## Optimización = minimizar
## Mejor fitness encontrado = 7.776137
## Mejor solución encontrada = -7.06204414367676 -0.714093213900924
## Valor función objetivo = -7.776137## [1] 7.776137 -6.095942 -1.387073 -6.159373 6.083099El vector devuelto contiene el fitness de cada uno de los individuos en el mismo orden que se encuentran en la matriz de la población.
Seleccionar individuos
La forma en que se seleccionan los individuos que participan en cada cruce difiere en las distintas implementaciones de los algoritmos genéticos. Por lo general, todas ellas tienden a favorecer la selección de aquellos individuos con mayor fitness. Algunas de las estrategias más comunes son:
Método de ruleta: la probabilidad de que un individuo sea seleccionado es proporcional a su fitness relativo, es decir, a su fitness dividido por la suma del fitness de todos los individuos de la población. Si el fitness de un individuo es el doble que el de otro, también lo será la probabilidad de que sea seleccionado. Este método presenta problemas si el fitness de unos pocos individuos es muy superior (varios órdenes de magnitud) al resto, ya que estos serán seleccionados de forma repetida y casi todos los individuos de la siguiente generación serán “hijos” de los mismos “padres” (poca variación).
Método rank: la probabilidad de selección de un individuo es inversamente proporcional a la posición que ocupa tras ordenar todos los individuos de mayor a menor fitness. Este método es menos agresivo que el método ruleta cuando la diferencia entre los mayores fitness es varios órdenes de magnitud superior al resto.
Selección competitiva (tournament): se seleccionan aleatoriamente dos parejas de individuos de la población (todos con la misma probabilidad). De cada pareja se selecciona el que tenga mayor fitness. Finalmente, se comparan los dos finalistas y se selecciona el de mayor fitness. Este método tiende a generar una distribución de la probabilidad de selección más equilibrada que las dos anteriores.
Selección truncada (truncated selection): se realizan selecciones aleatorias de individuos, habiendo descartado primero los n individuos con menor fitness de la población.
seleccionar_individuo <- function(vector_fitness, metodo_seleccion = "tournament",
verbose = FALSE) {
# Esta función recibe como argumento un vector con el fitness de cada individuo
# y selecciona una de las posiciones, donde la probabilidad de selección es
# proporcional al fitness.
# ARGUMENTOS
# ============================================================================
# vector_fitness: un vector con el fitness de cada individuo.
# metodo_seleccion: método para establecer la probabilidad de selección. Puede
# ser: "ruleta", "rank", o "tournament".
# verbose: mostrar información del proceso por pantalla.
#
# RETORNO
# ============================================================================
# El índice que ocupa el individuo seleccionado.
# COMPROBACIONES INICIALES
# ---------------------------------------------------------------------------
if (!metodo_seleccion %in% c("ruleta", "rank", "tournament")) {
stop("El método de selección debe de ser ruleta, rank o tournament.")
}
# SELECCIÓN DE INDIVIDUOS
# ----------------------------------------------------------------------------
# Se calcula la probabilidad de selección de cada individuo en función
# de su fitness.
if (metodo_seleccion == "ruleta") {
probabilidad_seleccion <- (vector_fitness) / sum(vector_fitness)
ind_seleccionado <- sample(
x = 1:length(vector_fitness),
size = 1,
prob = probabilidad_seleccion
)
} else if (metodo_seleccion == "rank") {
probabilidad_seleccion <- 1 / rank(-vector_fitness)
ind_seleccionado <- sample(
x = 1:length(vector_fitness),
size = 1,
prob = probabilidad_seleccion
)
} else if (metodo_seleccion == "tournament") {
# Se seleccionan aleatoriamente dos parejas de individuos.
ind_candidatos_a <- sample(x = 1:length(vector_fitness), size = 2)
ind_candidatos_b <- sample(x = 1:length(vector_fitness), size = 2)
# De cada pareja se selecciona el de mayor fitness.
ind_ganador_a <- ifelse(
vector_fitness[ind_candidatos_a[1]] > vector_fitness[ind_candidatos_a[2]],
ind_candidatos_a[1],
ind_candidatos_a[2]
)
ind_ganador_b <- ifelse(
vector_fitness[ind_candidatos_b[1]] > vector_fitness[ind_candidatos_b[2]],
ind_candidatos_b[1],
ind_candidatos_b[2]
)
# Se comparan los dos ganadores de cada pareja.
ind_seleccionado <- ifelse(
vector_fitness[ind_ganador_a] > vector_fitness[ind_ganador_b],
ind_ganador_a,
ind_ganador_b
)
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("----------------------", "\n")
cat("Individuo seleccionado", "\n")
cat("----------------------", "\n")
cat("Método selección =", metodo_seleccion, "\n")
cat("Índice seleccionado =", ind_seleccionado, "\n")
}
return(ind_seleccionado)
}Ejemplo
Se compara el patrón de selección del mejor individuo entre los métodos ruleta, rank y tournament. En primer lugar, se muestra un caso en el que la diferencia entre el mayor y el menor de los fitness no es muy acusada, y un segundo caso en el que sí lo es.
library(tidyverse)
fitness_poblacion <- c(20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5,
4, 3, 2, 1)
selecciones_ruleta <- rep(NA, times = 500)
for (i in 1:500) {
selecciones_ruleta[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "ruleta"
)
}
selecciones_ruleta <- data.frame(seleccion = selecciones_ruleta) %>%
mutate(metodo_seleccion = "ruleta")
selecciones_rank <- rep(NA, times = 500)
for (i in 1:500) {
selecciones_rank[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "rank"
)
}
selecciones_rank <- data.frame(seleccion = selecciones_rank) %>%
mutate(metodo_seleccion = "rank")
selecciones_tournament <- rep(NA, times = 500)
for (i in 1:500) {
selecciones_tournament[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "tournament"
)
}
selecciones_tournament <- data.frame(seleccion = selecciones_tournament) %>%
mutate(metodo_seleccion = "tournament")
bind_rows(selecciones_ruleta, selecciones_rank, selecciones_tournament) %>%
ggplot(aes(x= seleccion)) +
geom_bar(fill = "#3182bd") +
facet_grid(. ~ metodo_seleccion) +
theme_bw()
Cuando no existe una gran diferencia entre el individuo de mayor fitness y el resto, con el método rank, el individuo con mayor fitness se selecciona con mucha más frecuencia que el resto. Con los otros dos métodos, la probabilidad de selección decae de forma más gradual.
fitness_poblacion <- c(200, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5,
4, 3, 2, 1)
selecciones_ruleta <- rep(NA, times = 500)
for (i in 1:500) {
selecciones_ruleta[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "ruleta"
)
}
selecciones_ruleta <- data.frame(seleccion = selecciones_ruleta) %>%
mutate(metodo_seleccion = "ruleta")
selecciones_rank <- rep(NA, times = 500)
for (i in 1:500) {
selecciones_rank[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "rank"
)
}
selecciones_rank <- data.frame(seleccion = selecciones_rank) %>%
mutate(metodo_seleccion = "rank")
selecciones_tournament <- rep(NA, times = 500)
for (i in 1:500) {
selecciones_tournament[i] <- seleccionar_individuo(
vector_fitness = fitness_poblacion,
metodo_seleccion = "tournament"
)
}
selecciones_tournament <- data.frame(seleccion = selecciones_tournament) %>%
mutate(metodo_seleccion = "tournament")
bind_rows(selecciones_ruleta, selecciones_rank, selecciones_tournament) %>%
ggplot(aes(x = seleccion)) +
geom_bar(fill = "#3182bd") +
facet_grid(. ~ metodo_seleccion) +
theme_bw()
Cuando existe una gran diferencia entre el individuo de mayor fitness y el resto (uno o varios órdenes de magnitud), con el método ruleta, el individuo con mayor fitness se selecciona con mucha más frecuencia que el resto. A diferencia del caso anterior, en esta situación, la probabilidad de selección decae de forma más gradual con los métodos rank y tournament.
Teniendo en cuenta los comportamientos de selección de cada método, el método tournament parece ser la opción más equilibrada.
Cruzar dos individuos (crossover, recombinación)
El objetivo de esta etapa es generar, a partir de individuos ya existentes (parentales), nuevos individuos (descendencia) que combinen las características de los anteriores. Este es otro de los puntos del algoritmo en los que se puede seguir varias estrategias. Tres de las más empleadas son:
Cruzamiento a partir de uno solo punto: se selecciona aleatoriamente una posición que actúa como punto de corte. Cada individuo parental se divide en dos partes y se intercambian las mitades. Como resultado de este proceso, por cada cruce, se generan dos nuevos individuos.
Cruzamiento a partir múltiples puntos: se seleccionan aleatoriamente varias posiciones que actúan como puntos de corte. Cada individuo parental se divide por los puntos de corte y se intercambian las partes. Como resultado de este proceso, por cada cruce, se generan dos nuevos individuos.
Cruzamiento uniforme: el valor que toma cada posición del nuevo individuo se obtiene de uno de los dos parentales. Por lo general, la probabilidad de que el valor proceda de cada parental es la misma, aunque podría, por ejemplo, estar condicionada al fitness de cada uno. A diferencia de las anteriores estrategias, con esta, de cada cruce se genera un único descendiente.
cruzar_individuos <- function(parental_1,
parental_2,
metodo_cruce = "uniforme",
verbose = TRUE) {
# Esta función devuelve un individuo resultado de cruzar dos individuos
# parentales con el método de cruzamiento uniforme o punto simple.
#
# ARGUMENTOS
# ============================================================================
# parental_1: vector que representa a un individuo.
# parental_2: vector que representa a un individuo.
# metodo_cruce: estrategia de cruzamiento (uniforme", "punto_simple")
# RETORNO
# ============================================================================
# Un vector que representa a un nuevo individuo.
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
if (length(parental_1) != length(parental_2)) {
stop(paste0(
"La longitud de los dos vectores que representan a los ",
"individuos debe ser la misma."
))
}
if (!(metodo_cruce %in% c("uniforme", "punto_simple"))) {
stop("El método de cruzamiento debe ser: uniforme o punto_simple.")
}
# CRUCE
# ----------------------------------------------------------------------------
# Se crea el vector que representa el nuevo individuo
descendencia <- rep(NA, times = length(parental_1))
if (metodo_cruce == "uniforme") {
# Se seleccionan aleatoriamente las posiciones que se heredan del parental_1.
herencia_parent_1 <- sample(
x = c(TRUE, FALSE),
size = length(parental_1),
replace = TRUE
)
# El resto de posiciones se heredan del parental_2.
herencia_parent_2 <- !herencia_parent_1
descendencia[herencia_parent_1] <- parental_1[herencia_parent_1]
descendencia[herencia_parent_2] <- parental_2[herencia_parent_2]
} else {
punto_cruce <- sample(
x = 2:length(parental_1),
size = 1
)
descendencia <- c(
parental_1[1:(punto_cruce - 1)],
parental_2[punto_cruce:length(parental_1)]
)
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("---------------", "\n")
cat("Cruce realizado", "\n")
cat("---------------", "\n")
cat("Método =", metodo_cruce, "\n")
cat("Parental 1 = ", "\n")
cat(parental_1, "\n")
cat("Parental 2 = ", "\n")
cat(parental_2, "\n")
cat("Descendencia = ", "\n")
cat(descendencia, "\n")
cat("\n")
}
return(descendencia)
}Ejemplo
Se obtiene un nuevo individuo a partir del cruce de los individuos c(T, T, T, T, T) y c(F, F, F, F, F, F).
cruzar_individuos(
parental_1 = c(T, T, T, T, T),
parental_2 = c(F, F, F, F, F),
metodo = "uniforme",
verbose = TRUE
)## ---------------
## Cruce realizado
## ---------------
## Método = uniforme
## Parental 1 =
## TRUE TRUE TRUE TRUE TRUE
## Parental 2 =
## FALSE FALSE FALSE FALSE FALSE
## Descendencia =
## FALSE TRUE FALSE TRUE FALSE## [1] FALSE TRUE FALSE TRUE FALSEMutar individuo
Tras generar cada nuevo individuo de la descendencia, este se somete a un proceso de mutación en el que, cada una de sus posiciones, puede verse modificada con una probabilidad p. Este paso es importante para añadir diversidad al proceso y evitar que el algoritmo caiga en mínimos locales por que todos los individuos sean demasiado parecidos de una generación a otra.
Existen diferentes estrategias para controlar la magnitud del cambio que puede provocar una mutación.
Distribución uniforme: la mutación de la posición i se consigue sumándole al valor de i un valor extraído de una distribución uniforme, por ejemplo una entre [-1,+1].
Distribución normal: la mutación de la posición i se consigue sumándole al valor de i un valor extraído de una distribución normal, comúnmente centrada en 0 y con una determinada desviación estándar. Cuanto mayor la desviación estándar, con mayor probabilidad la mutación introducirá cambios grandes.
Aleatorio: la mutación de la posición i se consigue reemplazando el valor de i por nuevo valor aleatorio dentro del rango permitido para esa variable. Esta estrategia suele conllevar mayores variaciones que las dos anteriores.
Hay que tener en cuenta que, debido a las mutaciones, un valor que inicialmente estaba dentro del rango permitido puede salirse de él. Una forma de evitarlo es: si el valor tras la mutación excede alguno de los límites acotados, se sobrescribe con el valor del límite. Es decir, se permite que los valores se alejen como máximo hasta el límite impuesto.
mutar_individuo <- function(individuo, limite_inf, limite_sup,
prob_mut = 0.01, distribucion = "uniforme",
media_distribucion = 1, sd_distribucion = 1,
min_distribucion = -1, max_distribucion = 1,
verbose = TRUE) {
# ARGUMENTOS
# =============================================================================
# individuo: vector que representa a un individuo.
# prob_mut: probabilidad que tiene cada posición del individuo de mutar.
# distribucion: distribución de la que obtener el factor de mutación. Puede
# ser: "normal", "uniforme" o "aleatoria".
# media_distribucion: media de la distribución si se selecciona
# distribucion = "normal".
# sd_distribucion: desviación estándar de la distribución si se selecciona
# distribucion = "normal".
# min_distribucion: mínimo la distribución si se selecciona
# distribucion = "uniforme".
# max_distribucion: máximo la distribución si se selecciona
# distribucion = "uniforme".
#
# RETORNO
# ============================================================================
# Un vector que representa al individuo tras someterse a las mutaciones.
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
if (!(distribucion %in% c("normal", "uniforme", "aleatoria"))) {
stop("El argumento distribución debe ser: normal, uniforme o aleatoria.")
}
# CRUCE
# ----------------------------------------------------------------------------
# Selección de posiciones a mutar.
posiciones_mutadas <- runif(n = length(individuo), min = 0, max = 1) < prob_mut
# Se modifica el valor de aquellas posiciones que hayan sido seleccionadas para
# mutar. Si el valor de prob_mut es muy bajo, las mutaciones serán muy poco
# frecuentes y el individuo devuelto será casi siempre igual al original.
# Si se emplea distribucion = "uniforme" o distribucion = "normal":
if (distribucion == "normal" | distribucion == "uniforme") {
# Se extrae un valor aleatorio de la distribución elegida que se suma
# para modificar la/las posiciones mutadas.
if (distribucion == "normal") {
factor_mut <- rnorm(
n = sum(posiciones_mutadas),
mean = media_distribucion,
sd = sd_distribucion
)
}
if (distribucion == "uniforme") {
factor_mut <- runif(
n = sum(posiciones_mutadas),
min = min_distribucion,
max = max_distribucion
)
}
individuo[posiciones_mutadas] <- individuo[posiciones_mutadas] + factor_mut
# Se comprueba si algún valor mutado supera los límites impuestos. En tal caso
# se sobrescribe con el valor del límite correspondiente.
for (i in which(posiciones_mutadas)) {
if (individuo[i] < limite_inf[i]) {
individuo[i] <- limite_inf[i]
}
if (individuo[i] > limite_sup[i]) {
individuo[i] <- limite_sup[i]
}
}
} else if (distribucion == "aleatoria") {
for (i in which(posiciones_mutadas)) {
individuo[i] <- runif(n = 1, min = limite_inf[i], max = limite_sup[i])
}
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("-----------------", "\n")
cat("Individuo mutado", "\n")
cat("-----------------", "\n")
cat("Probabilidad =", prob_mut, "\n")
cat("Individuo = ", individuo, "\n")
cat("\n")
}
return(individuo)
}Ejemplo
Se somete a un individuo al proceso de mutación, con una probabilidad de mutación de 0.5 y con límites para el valor que puede tener cada posición.
mutar_individuo(
individuo = c(3, 3, 3, 3, 3, 3),
prob_mut = 0.5,
distribucion = "aleatoria",
media_distribucion = 1,
sd_distribucion = 1,
limite_inf = c(-5, -5, -5, -5, -5, -5),
limite_sup = c(5, 5, 5, 5, 5, 5)
)## -----------------
## Individuo mutado
## -----------------
## Probabilidad = 0.5
## Individuo = 3 3 -0.1000917 -1.04316 -2.22533 -0.9442027## [1] 3.0000000 3.0000000 -0.1000917 -1.0431598 -2.2253300 -0.9442027Algoritmo completo
En cada uno de los apartados anteriores se ha definido una de las etapas del algoritmo genético. A continuación, se combinan todas ellas dentro de una única función.
optimizar_ga <- function(
funcion_objetivo,
n_variables,
optimizacion,
limite_inf = NULL,
limite_sup = NULL,
n_poblacion = 20,
n_generaciones = 50,
elitismo = 0.1,
prob_mut = 0.01,
distribucion = "uniforme",
media_distribucion = 1,
sd_distribucion = 1,
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
verbose = 1,
...) {
# ARGUMENTOS
# =============================================================================
# funcion_objetivo: nombre de la función que se desea optimizar. Debe de haber
# sido definida previamente.
# n_variables: longitud de los individuos.
# optimizacion: "maximizar" o "minimizar". Dependiendo de esto, la relación
# del fitness es directamente o indirectamente proporcional al
# valor de la función.
# limite_inf: vector con el límite inferior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quiere acotar.
# limite_sup: vector con el límite superior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quieren acotar.
# n_poblacion: número total de individuos de la población.
# n_generaciones: número total de generaciones creadas.
# elitismo: porcentaje de mejores individuos de la población actual que
# pasan directamente a la siguiente población.
# prob_mut: probabilidad que tiene cada posición del individuo de mutar.
# distribucion: distribución de la que obtener el factor de mutación. Puede
# ser: "normal", "uniforme" o "aleatoria".
# media_distribucion: media de la distribución si se selecciona distribucion="normal".
# sd_distribucion: desviación estándar de la distribución si se selecciona
# distribucion="normal".
# min_distribucion: mínimo la distribución si se selecciona distribucion="uniforme".
# max_distribucion: máximo la distribución si se selecciona distribucion="uniforme".
# metodo_seleccion: método para establecer la probabilidad de selección. Puede
# ser: "ruleta", "rank" o "tournament".
# metodo_seleccion: método para cruzar los individuos. Puede ser: "uniforme",
# "punto_simple".
# parada_temprana: si durante las últimas "rondas_parada" generaciones la diferencia
# absoluta entre mejores individuos no es superior al valor de
# "tolerancia_parada", se detiene el algoritmo y no se crean
# nuevas generaciones.
# rondas_parada: número de generaciones consecutivas sin mejora mínima para que
# se active la parada temprana.
# tolerancia_parada: valor mínimo que debe tener la diferencia de generaciones
# consecutivas para considerar que hay cambio.
# verbose: Nivel de detalle para que se imprima por pantalla el
# resultado de cada paso del algoritmo (0, 1, 2)
# RETORNO
# =============================================================================
# La función devuelve una lista con 5 elementos:
# fitness: una lista con el fitness del mejor individuo de cada
# generación.
# mejores_individuos: una lista con la combinación de predictores del mejor
# individuo de cada generación.
# mejor_individuo: combinación de predictores del mejor individuo encontrado
# en todo el proceso.
# diferencia_abs: una lista con la diferencia absoluta entre el fitness
# del mejor individuo de generaciones consecutivas.
# df_resultados: un dataframe con todos los resultados anteriores.
start_time <- Sys.time()
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
# Si se activa la parada temprana, hay que especificar los argumentos
# rondas_parada y tolerancia_parada.
if (isTRUE(parada_temprana) &
(is.null(rondas_parada) | is.null(tolerancia_parada)) ) {
stop(paste(
"Para activar la parada temprana es necesario indicar un valor",
"de rondas_parada y de tolerancia_parada."
))
}
# ESTABLECER LOS LÍMITES DE BÚSQUEDA SI EL USUARIO NO LO HA HECHO
# ----------------------------------------------------------------------------
if (is.null(limite_sup) | is.null(limite_inf)) {
warning(paste(
"Es altamente recomendable indicar los límites dentro de los",
"cuales debe buscarse la solución de cada variable.",
"Por defecto se emplea: [-10^3, 10^3]."
))
}
if (any(
is.null(limite_sup), is.null(limite_inf), any(is.na(limite_sup)),
any(is.na(limite_inf))
)) {
warning(paste(
"Los límites empleados por defecto cuando no se han definido son:",
" [-10^3, 10^3]."
))
cat("\n")
}
# Si no se especifica limite_inf, el valor mínimo que pueden tomar las variables
# es -10^3.
if (is.null(limite_inf)) {
limite_inf <- rep(x = -10^3, times = n_variables)
}
# Si no se especifica limite_sup, el valor máximo que pueden tomar las variables
# es 10^3.
if (is.null(limite_sup)) {
limite_sup <- rep(x = 10^3, times = n_variables)
}
# Si los límites no son nulos, se reemplazan aquellas posiciones NA por el valor
# por defecto -10^3 y 10^3.
if (!is.null(limite_inf)) {
limite_inf[is.na(limite_inf)] <- -10^3
}
if (!is.null(limite_sup)) {
limite_sup[is.na(limite_sup)] <- 10^3
}
# ALMACENAMIENTO DE RESULTADOS
# ----------------------------------------------------------------------------
# Por cada generación se almacena, la población, el mejor individuo, su fitness,
# y la diferencia absoluta respecto a la última generación.
poblaciones <- vector(mode = "list", length = n_generaciones)
resultados_fitness <- vector(mode = "list", length = n_generaciones)
resultados_individuo <- vector(mode = "list", length = n_generaciones)
diferencia_abs <- vector(mode = "list", length = n_generaciones)
# ITERACIÓN DE POBLACIONES
# ----------------------------------------------------------------------------
for (i in 1:n_generaciones) {
if (verbose %in% c(1,2)) {
cat("-------------------", "\n")
cat("Generación:", paste0(i, "\\", n_generaciones), "\n")
cat("-------------------", "\n")
}
if (i == 1) {
# CREACIÓN DE LA POBLACIÓN INICIAL
# ------------------------------------------------------------------------
poblacion <- crear_poblacion(
n_poblacion = n_poblacion,
n_variables = n_variables,
limite_inf = limite_inf,
limite_sup = limite_sup,
verbose = verbose %in% c(2)
)
}
# CALCULAR FITNESS DE LOS INDIVIDUOS DE LA POBLACIÓN
# --------------------------------------------------------------------------
fitness_ind_poblacion <- calcular_fitness_poblacion(
poblacion = poblacion,
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
verbose = verbose %in% c(2)
)
# SE ALMACENA LA POBLACIÓN Y SU MEJOR INDIVIDUO
# --------------------------------------------------------------------------
poblaciones[[i]] <- poblacion
fitness_mejor_individuo <- max(fitness_ind_poblacion)
mejor_individuo <- poblacion[which.max(fitness_ind_poblacion), ]
resultados_fitness[[i]] <- fitness_mejor_individuo
resultados_individuo[[i]] <- mejor_individuo
# SE CALCULA LA DIFERENCIA ABSOLUTA RESPECTO A LA GENERACIÓN ANTERIOR
# --------------------------------------------------------------------------
# La diferencia solo puede calcularse a partir de la segunda generación.
if (i > 1) {
diferencia_abs[[i]] <- abs(resultados_fitness[[i - 1]] - resultados_fitness[[i]])
}
# NUEVA POBLACIÓN
# --------------------------------------------------------------------------
nueva_poblacion <- matrix(
data = NA,
nrow = nrow(poblacion),
ncol = ncol(poblacion)
)
# ELITISMO
# --------------------------------------------------------------------------
# El elitismo indica el porcentaje de mejores individuos de la población
# actual que pasan directamente a la siguiente población. De esta forma, se
# asegura que, la siguiente generación, no sea nunca inferior.
if (elitismo > 0) {
n_elitismo <- ceiling(nrow(poblacion) * elitismo)
posicion_n_mejores <- order(fitness_ind_poblacion, decreasing = TRUE)
posicion_n_mejores <- posicion_n_mejores[1:n_elitismo]
nueva_poblacion[1:n_elitismo, ] <- poblacion[posicion_n_mejores, ]
} else {
n_elitismo <- 0
}
# CREACIÓN DE NUEVOS INDIVIDUOS POR CRUCES
# --------------------------------------------------------------------------
for (j in (n_elitismo + 1):nrow(nueva_poblacion)) {
# Seleccionar parentales
indice_parental_1 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metodo_seleccion = metodo_seleccion,
verbose = verbose %in% c(2)
)
indice_parental_2 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metodo_seleccion = metodo_seleccion,
verbose = verbose %in% c(2)
)
parental_1 <- poblacion[indice_parental_1, ]
parental_2 <- poblacion[indice_parental_2, ]
# Cruzar parentales para obtener la descendencia
descendencia <- cruzar_individuos(
parental_1 = parental_1,
parental_2 = parental_2,
metodo_cruce = metodo_cruce,
verbose = verbose %in% c(2)
)
# Mutar la descendencia
descendencia <- mutar_individuo(
individuo = descendencia,
prob_mut = prob_mut,
limite_inf = limite_inf,
limite_sup = limite_sup,
distribucion = distribucion,
media_distribucion = media_distribucion,
sd_distribucion = sd_distribucion,
min_distribucion = min_distribucion,
max_distribucion = max_distribucion,
verbose = verbose %in% c(2)
)
nueva_poblacion[j, ] <- descendencia
}
poblacion <- nueva_poblacion
# CRITERIO DE PARADA
# --------------------------------------------------------------------------
# Si durante las últimas n generaciones, la diferencia absoluta entre mejores
# individuos no es superior al valor de tolerancia_parada, se detiene el
# algoritmo y no se crean nuevas generaciones.
if (parada_temprana && (i > rondas_parada)) {
ultimos_n <- tail(unlist(diferencia_abs), n = rondas_parada)
if (all(ultimos_n < tolerancia_parada)) {
cat(
"Algoritmo detenido en la generacion", i,
"por falta cambio mínimo de", tolerancia_parada,
"durante", rondas_parada,
"generaciones consecutivas.",
"\n"
)
break()
}
}
}
# IDENTIFICACIÓN DEL MEJOR INDIVIDUO DE TODO EL PROCESO
# ----------------------------------------------------------------------------
indice_mejor_individuo_global <- which.max(unlist(resultados_fitness))
mejor_fitness_global <- resultados_fitness[[indice_mejor_individuo_global]]
mejor_individuo_global <- resultados_individuo[[indice_mejor_individuo_global]]
# Se identifica el valor de la función objetivo para el mejor individuo.
if (optimizacion == "maximizar") {
mejor_valor_global <- mejor_fitness_global
} else {
mejor_valor_global <- -1*mejor_fitness_global
}
# RESULTADOS
# ----------------------------------------------------------------------------
# Para crear el dataframe se convierten las listas a vectores del mismo tamaño.
resultados_fitness <- unlist(resultados_fitness)
diferencia_abs <- c(NA, unlist(diferencia_abs))
# Si hay parada temprana, algunas generaciones no se alcanzan: Se eliminan sus
# posiciones de las listas de resultados
resultados_individuo <- resultados_individuo[!sapply(resultados_individuo, is.null)]
poblaciones <- poblaciones[!sapply(poblaciones, is.null)]
# Para poder añadir al dataframe la secuencia variables, se concatenan.
variables <- sapply(
X = resultados_individuo,
FUN = function(x) {
paste(x, collapse = ", ")
}
)
df_resultados <- data.frame(
generacion = seq_along(resultados_fitness),
fitness = resultados_fitness,
predictores = variables,
diferencia_abs = diferencia_abs
)
resultados <- list(
mejor_individuo_global = mejor_individuo_global,
mejor_valor_global = mejor_valor_global,
mejor_fitness_por_generacion = resultados_fitness,
mejor_individuo_por_generacion = resultados_individuo,
diferencia_abs = diferencia_abs,
df_resultados = df_resultados,
poblaciones = poblaciones,
funcion_objetivo = funcion_objetivo
)
end_time <- Sys.time()
# INFORMACIÓN ALMACENADA EN LOS ATRIBUTOS
# ----------------------------------------------------------------------------
attr(resultados, "class") <- "optimizacion_ga"
attr(resultados, 'fecha_creacion') <- end_time
attr(resultados, 'duracion_optimizacion') <- paste(
difftime(end_time, start_time, "secs"),
"secs"
)
attr(resultados, 'optimizacion') <- optimizacion
attr(resultados, 'lim_inf') <- limite_inf
attr(resultados, 'lim_sup') <- limite_sup
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'generaciones') <- i
attr(resultados, 'valor_variables') <- mejor_individuo_global
attr(resultados, 'mejor_fitness') <- mejor_fitness_global
attr(resultados, 'optimo_encontrado') <- mejor_valor_global
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'elitismo') <- elitismo
attr(resultados, 'prob_mut') <- prob_mut
attr(resultados, 'metodo_seleccion') <- metodo_seleccion
attr(resultados, 'metodo_cruce') <- metodo_cruce
attr(resultados, 'parada_temprana') <- parada_temprana
attr(resultados, 'rondas_parada') <- rondas_parada
attr(resultados, 'tolerancia_parada') <- tolerancia_parada
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose %in% c(1,2)) {
cat("-----------------------", "\n")
cat("Optimización finalizada", "\n")
cat("-----------------------", "\n")
cat("Fecha finalización =", as.character(Sys.time()), "\n")
cat("Duración selección = ")
print(difftime(end_time, start_time))
cat("Número generaciones =", i, "\n")
cat("Límite inferior =", paste(limite_inf, collapse = ", "), "\n")
cat("Límite superior =", paste(limite_sup, collapse = ", "), "\n")
cat("Optimización =", optimizacion,"\n")
cat("Óptimo encontrado =", mejor_valor_global,"\n")
cat("Valor variables =", mejor_individuo_global, "\n")
cat("\n")
}
return(resultados)
}
print.optimizacion_ga <- function(obj){
# Función print para objetos optimizacion_ga
cat("----------------------------------------------", "\n")
cat("Resultados optimización por algoritmo genético", "\n")
cat("----------------------------------------------", "\n")
cat("Fecha creación =", attr(obj, 'fecha_creacion'), "\n")
cat("Duración selección = ", attr(obj, 'duracion_optimizacion'), "\n")
cat("Número generaciones =", attr(obj, 'generaciones'), "\n")
cat("Límite inferior =", attr(obj, 'lim_inf'), "\n")
cat("Límite superior =", attr(obj, 'lim_sup'), "\n")
cat("Optimización =", attr(obj, 'optimizacion'), "\n")
cat("Óptimo encontrado =", attr(obj, 'optimo_encontrado'), "\n")
cat("Valor variables =", attr(obj, 'valor_variables'), "\n")
cat("Función objetivo =", "\n")
cat("\n")
print(obj$funcion_objetivo)
}Ejemplo 1
En este ejemplo se pretende evaluar la capacidad del algoritmo genético para encontrar el mínimo de la función \(f(x_1,x_2) = x_1^2 + x_2^2\). El mínimo global de esta función puede obtenerse de forma analítica igualando las derivadas parciales a cero, lo que permite comparar el resultado obtenido. \[f(x_1=0, x_2 = 0) = 0\]
Función objetivo
Representación gráfica de la función.
x1 <- seq(-10, 10, length.out = 50)
x2 <- seq(-10, 10, length.out = 50)
datos <- expand.grid(x1 = x1, x2 = x2)
datos <- datos %>%
mutate(f_x = map2_dbl(x1, x2, .f = funcion))
ggplot(data = datos, aes(x = x1, y = x2, z = f_x)) +
geom_contour(aes(colour = stat(level)), bins = 20) +
labs(title = "f(x1,x2) = x1^2 + x2^2") +
theme_bw() +
theme(legend.position = "none")
x1 <- seq(-10, 10, length.out = 50)
x2 <- seq(-10, 10, length.out = 50)
f_x <- outer(x1, x2, FUN = funcion)
library(viridis)
colores <- viridis::magma(n = 100, alpha = 0.7)
z.facet.center <- (f_x[-1, -1] + f_x[-1, -ncol(f_x)] +
f_x[-nrow(f_x), -1] +
f_x[-nrow(f_x), -ncol(f_x)])/4
z.facet.range <- cut(z.facet.center, 100)
par(mai = c(0,0,0,0))
persp(x = x1, y = x2, z = f_x,
shade = 0.8,
phi = 30,
theta = 30,
col = colores[z.facet.range],
axes = FALSE)
Optimización
# Ejemplo ilustrativo, en un caso real, emplear como mínimo 100 individuos por
# generación.
resultados_ga <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-10, -10),
limite_sup = c(10, 10),
n_poblacion = 30,
n_generaciones = 500,
elitismo = 0.01,
prob_mut = 0.1,
distribucion = "uniforme",
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
parada_temprana = TRUE,
rondas_parada = 10,
tolerancia_parada = 10^-8,
verbose = 1
)## -------------------
## Generación: 1\500
## -------------------
## -------------------
## Generación: 2\500
## -------------------
## -------------------
## Generación: 3\500
## -------------------
## -------------------
## Generación: 4\500
## -------------------
## -------------------
## Generación: 5\500
## -------------------
## -------------------
## Generación: 6\500
## -------------------
## -------------------
## Generación: 7\500
## -------------------
## -------------------
## Generación: 8\500
## -------------------
## -------------------
## Generación: 9\500
## -------------------
## -------------------
## Generación: 10\500
## -------------------
## -------------------
## Generación: 11\500
## -------------------
## -------------------
## Generación: 12\500
## -------------------
## -------------------
## Generación: 13\500
## -------------------
## -------------------
## Generación: 14\500
## -------------------
## -------------------
## Generación: 15\500
## -------------------
## -------------------
## Generación: 16\500
## -------------------
## -------------------
## Generación: 17\500
## -------------------
## -------------------
## Generación: 18\500
## -------------------
## -------------------
## Generación: 19\500
## -------------------
## -------------------
## Generación: 20\500
## -------------------
## -------------------
## Generación: 21\500
## -------------------
## -------------------
## Generación: 22\500
## -------------------
## -------------------
## Generación: 23\500
## -------------------
## Algoritmo detenido en la generacion 23 por falta cambio mínimo de 1e-08 durante 10 generaciones consecutivas.
## -----------------------
## Optimización finalizada
## -----------------------
## Fecha finalización = 2020-11-12 12:21:46
## Duración selección = Time difference of 0.06052494 secs
## Número generaciones = 23
## Límite inferior = -10, -10
## Límite superior = 10, 10
## Optimización = minimizar
## Óptimo encontrado = 0.0006697396
## Valor variables = -0.02152595 0.01436569Resultados
El objeto devuelto por la función optimizar_ga almacena información sobre el proceso de optimización y de los resultados finales.
## ----------------------------------------------
## Resultados optimización por algoritmo genético
## ----------------------------------------------
## Fecha creación = 1605183707
## Duración selección = 0.0605249404907227 secs
## Número generaciones = 23
## Límite inferior = -10 -10
## Límite superior = 10 10
## Optimización = minimizar
## Óptimo encontrado = 0.0006697396
## Valor variables = -0.02152595 0.01436569
## Función objetivo =
##
## function(x1, x2){
## return(x1^2 + x2^2)
## }
## <bytecode: 0x55c773e1bd18>## $names
## [1] "mejor_individuo_global" "mejor_valor_global"
## [3] "mejor_fitness_por_generacion" "mejor_individuo_por_generacion"
## [5] "diferencia_abs" "df_resultados"
## [7] "poblaciones" "funcion_objetivo"
##
## $class
## [1] "optimizacion_ga"
##
## $fecha_creacion
## [1] "2020-11-12 12:21:46 UTC"
##
## $duracion_optimizacion
## [1] "0.0605249404907227 secs"
##
## $optimizacion
## [1] "minimizar"
##
## $lim_inf
## [1] -10 -10
##
## $lim_sup
## [1] 10 10
##
## $n_poblacion
## [1] 30
##
## $generaciones
## [1] 23
##
## $valor_variables
## [1] -0.02152595 0.01436569
##
## $mejor_fitness
## [1] -0.0006697396
##
## $optimo_encontrado
## [1] 0.0006697396
##
## $elitismo
## [1] 0.01
##
## $prob_mut
## [1] 0.1
##
## $metodo_seleccion
## [1] "tournament"
##
## $metodo_cruce
## [1] "uniforme"
##
## $parada_temprana
## [1] TRUE
##
## $rondas_parada
## [1] 10
##
## $tolerancia_parada
## [1] 1e-08Mejor individuo
## [1] -0.02152595 0.01436569## [1] 0.0006697396ggplot(data = datos, aes(x = x1, y = x2, z = f_x)) +
geom_contour(aes(colour = stat(level)), bins = 20) +
annotate(
geom = "point",
x = resultados_ga$mejor_individuo_global[1],
y = resultados_ga$mejor_individuo_global[2],
color = "red",
size = 3
) +
labs(title = "f(x1,x2) = x1^2 + x2^2") +
theme_bw() +
theme(legend.position = "none")
Evolución del fitness
En el siguiente gráfico se puede ver cómo evoluciona el fitness del mejor individuo a medida que avanzan las generaciones.
library(ggplot2)
ggplot(data = resultados_ga$df_resultados,
aes(x = generacion, y = fitness)) +
geom_line(aes(group = 1)) +
geom_point() +
labs(title = "Evolución del fitness a lo largo de las generaciones") +
theme_bw()
Animación de cómo avanza la búsqueda del mínimo
library(animation)
# Se extrae la posición de cada partícula en cada iteración.
extraer_individuos <- function(optimizacion_ga){
individuos <- do.call(rbind, optimizacion_ga$poblaciones)
individuos <- as.data.frame(individuos)
individuos$generacion <- rep(
x = 1:attr(optimizacion_ga, 'generaciones'),
each = attr(optimizacion_ga, 'n_poblacion')
)
return(individuos)
}
historico_inidviduos <- extraer_individuos(resultados_ga)
animation::saveGIF(
for (i in 1:max(historico_inidviduos$generacion)) {
p <- ggplot() +
geom_contour(data = datos,
aes(x = x1, y = x2, z = f_x, colour = stat(level)),
bins = 20) +
geom_point(data = historico_inidviduos %>% filter(generacion == i),
aes(x = V1, y = V2),
color = "red",
size = 1.8) +
labs(title = paste("Generación:", i)) +
theme_bw() +
theme(legend.position = "none")
# Al ser un gráfico ggplot2 es necesario hacer el print()
print(p)
},
# Nombre del gif
movie.name = "gift_ggplot_1.gif",
# Dimensiones
ani.width = 500,
ani.height = 300,
# Tiempo de duración de cada frame (segundos)
interval = 0.5
)
Ejemplo 2
En este ejemplo se pretende evaluar la capacidad del algoritmo genético para encontrar el mínimo de la función de Mishra Bird.
\[f(x_1,x_2) = sin(x_2)exp(1-cos(x_1))^2 + cos(x_1)exp(1-sin(x_2))^2 + (x_1-x_2)^2\]
Para la región acotada entre: \[-10 \leq x_1 \leq 0\] \[-6.5 \leq x_2 \leq 0\]
la función tiene múltiples mínimos locales y un único el mínimo global que se encuentra en: \[f(-3.1302468, -1.5821422) = -106.7645367\]
Función objetivo
# Función objetivo a optimizar.
funcion <- function(x1, x2){
sin(x2)*exp(1-cos(x1))^2 + cos(x1)*exp(1-sin(x2))^2 + (x1-x2)^2
}Representación gráfica de la función.
x1 <- seq(-10, 0, length.out = 100)
x2 <- seq(-6.5, 0, length.out = 100)
datos <- expand.grid(x1 = x1, x2 = x2)
datos <- datos %>%
mutate(f_x = map2_dbl(x1, x2, .f = funcion))
ggplot(data = datos, aes(x = x1, y = x2, z = f_x)) +
geom_contour(aes(colour = stat(level)), bins = 20) +
labs(title = "f(x1,x2)") +
theme_bw() +
theme(legend.position = "none",
title = element_text(size = 10))
x1 <- seq(-10, 0, length.out = 50)
x2 <- seq(-6.5, 0, length.out = 50)
f_x <- outer(x1, x2, FUN = funcion)
library(viridis)
colores <- viridis::magma(n = 100, alpha = 0.7)
z.facet.center <- (f_x[-1, -1] + f_x[-1, -ncol(f_x)] +
f_x[-nrow(f_x), -1] +
f_x[-nrow(f_x), -ncol(f_x)])/4
z.facet.range <- cut(z.facet.center, 100)
par(mai = c(0,0,0,0))
persp(x = x1, y = x2, z = f_x,
shade = 0.8,
r = 8,
phi = 10,
theta = 25,
col = colores[z.facet.range],
axes = FALSE)
Optimización
# Ejemplo ilustrativo, en un caso real, emplear como mínimo 100 individuos por
# generación.
resultados_ga <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-10, -6.5),
limite_sup = c(0, 0),
n_poblacion = 30,
n_generaciones = 500,
elitismo = 0.01,
prob_mut = 0.1,
distribucion = "uniforme",
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
parada_temprana = TRUE,
rondas_parada = 10,
tolerancia_parada = 10^-8,
verbose = 0
)## Algoritmo detenido en la generacion 33 por falta cambio mínimo de 1e-08 durante 10 generaciones consecutivas.Resultados
El objeto devuelto por la función optimizar_ga almacena la información (fitness, valor de las variables,…) del mejor individuo de cada generación.
## ----------------------------------------------
## Resultados optimización por algoritmo genético
## ----------------------------------------------
## Fecha creación = 1605183714
## Duración selección = 0.0976943969726562 secs
## Número generaciones = 33
## Límite inferior = -10 -6.5
## Límite superior = 0 0
## Optimización = minimizar
## Óptimo encontrado = -106.7445
## Valor variables = -3.103213 -1.600978
## Función objetivo =
##
## function(x1, x2){
## sin(x2)*exp(1-cos(x1))^2 + cos(x1)*exp(1-sin(x2))^2 + (x1-x2)^2
## }
## <bytecode: 0x55c7745ca440>Mejor individuo
## [1] -3.103213 -1.600978ggplot(data= datos, aes(x = x1, y = x2, z = f_x)) +
geom_contour(aes(colour = stat(level)), bins = 20) +
annotate(geom = "point",
x = resultados_ga$mejor_individuo_global[1],
y = resultados_ga$mejor_individuo_global[2],
color = "red",
size = 2.5) +
labs(title = "Posición del mínimo encontrado") +
theme_bw() +
theme(legend.position = "none",
title = element_text(size = 10))
Evolución del fitness
En el siguiente gráfico se puede ver cómo evoluciona el fitness del mejor individuo a medida que avanzan las generaciones.
library(ggplot2)
ggplot(data = resultados_ga$df_resultados,
aes(x = generacion, y = fitness)) +
geom_line(aes(group = 1)) +
geom_point() +
labs(title = "Evolución del fitness a lo largo de las generaciones") +
theme_bw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
Animación de cómo avanza la búsqueda del mínimo
library(animation)
# Se extrae la posición de cada partícula en cada iteración.
extraer_individuos <- function(optimizacion_ga){
individuos <- do.call(rbind, optimizacion_ga$poblaciones)
individuos <- as.data.frame(individuos)
individuos$generacion <- rep(
x = 1:attr(optimizacion_ga, 'generaciones'),
each = attr(optimizacion_ga, 'n_poblacion')
)
return(individuos)
}
historico_inidviduos <- extraer_individuos(resultados_ga)
animation::saveGIF(
for (i in 1:max(historico_inidviduos$generacion)) {
p <- ggplot() +
geom_contour(data = datos,
aes(x = x1, y = x2, z = f_x, colour = stat(level)),
bins = 20) +
geom_point(data = historico_inidviduos %>% filter(generacion == i),
aes(x = V1, y = V2),
color = "red",
size = 1.8) +
labs(title = paste("Generación:", i)) +
theme_bw() +
theme(legend.position = "none")
# Al ser un gráfico ggplot2 es necesario hacer el print()
print(p)
},
# Nombre del gif
movie.name = "gift_ggplot_2.gif",
# Dimensiones
ani.width = 500,
ani.height = 300,
# Tiempo de duración de cada frame (segundos)
interval = 0.5
)
Ejemplo 3
En este ejemplo se pretende evaluar la capacidad del algoritmo genético para encontrar el mínimo de la función de Ackley.
\[f(x_1,x_2) =-20exp[-0.2\sqrt{0.5(x_1^2 + x_2^2)}] - exp[0.5(\cos 2\pi x_1 + \cos 2\pi x_2)] + e + 20\]
la función tiene múltiples mínimos locales y un único el mínimo global que se encuentra en: \[f(0, 0) = 0\]
Función objetivo
# Función objetivo a optimizar.
funcion <- function(x1, x2){
-20*exp(-0.7*sqrt(0.5*(x1^2 + x2^2))) - exp(0.5*(cos(2*pi*x1) + cos(2*pi*x2))) + exp(1) + 20
}Representación gráfica de la función.
x1 <- seq(-5, 5, length.out = 100)
x2 <- seq(-5, 5, length.out = 100)
datos <- expand.grid(x1 = x1, x2 = x2)
datos <- datos %>%
mutate(f_x = map2_dbl(x1, x2, .f = funcion))
ggplot(data = datos, aes(x = x1, y = x2, z = f_x)) +
geom_contour(aes(colour = stat(level)), bins = 20) +
labs(title = "f(x1, x2)") +
theme_bw() +
theme(legend.position = "none",
title = element_text(size = 10))
x1 <- seq(-5, 5, length.out = 50)
x2 <- seq(-5, 5, length.out = 50)
f_x <- outer(x1, x2, FUN = funcion)
library(viridis)
colores <- viridis::magma(n = 100, alpha = 0.9)
z.facet.center <- (f_x[-1, -1] + f_x[-1, -ncol(f_x)] +
f_x[-nrow(f_x), -1] +
f_x[-nrow(f_x), -ncol(f_x)])/4
z.facet.range <- cut(z.facet.center, 100)
par(mai = c(0,0,0,0))
persp(x = x1, y = x2, z = f_x,
shade = 0.8,
r = 1,
phi = 25,
theta = 25,
col = colores[z.facet.range],
axes = FALSE)
Optimización
resultados_ga <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-5, -5),
limite_sup = c(5, 5),
n_poblacion = 200,
n_generaciones = 1000,
elitismo = 0.01,
prob_mut = 0.1,
distribucion = "uniforme",
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
parada_temprana = TRUE,
rondas_parada = 10,
tolerancia_parada = 10^-8,
verbose = 0
)## Algoritmo detenido en la generacion 35 por falta cambio mínimo de 1e-08 durante 10 generaciones consecutivas.Resultados
El objeto devuelto por la función optimizar_ga almacena la información (fitness, valor de las variables,…) del mejor individuo de cada generación.
## ----------------------------------------------
## Resultados optimización por algoritmo genético
## ----------------------------------------------
## Fecha creación = 1605183725
## Duración selección = 0.611183881759644 secs
## Número generaciones = 35
## Límite inferior = -5 -5
## Límite superior = 5 5
## Optimización = minimizar
## Óptimo encontrado = 0.02512501
## Valor variables = -0.002464092 -0.0005389419
## Función objetivo =
##
## function(x1, x2){
## -20*exp(-0.7*sqrt(0.5*(x1^2 + x2^2))) - exp(0.5*(cos(2*pi*x1) + cos(2*pi*x2))) + exp(1) + 20
## }
## <bytecode: 0x55c774f811c0>Mejor individuo
## [1] -0.0024640919 -0.0005389419ggplot(data = datos, aes(x = x1, y = x2, z = f_x)) +
geom_contour(aes(colour = stat(level)), bins = 20) +
annotate(geom = "point",
x = resultados_ga$mejor_individuo_global[1],
y = resultados_ga$mejor_individuo_global[2],
color = "red",
size = 2.5) +
labs(title = "Posición del mínimo encontrado") +
theme_bw() +
theme(legend.position = "none",
title = element_text(size = 10))
Evolución del fitness
En el siguiente gráfico se puede ver cómo evoluciona el fitness del mejor individuo a medida que avanzan las generaciones.
library(ggplot2)
ggplot(data = resultados_ga$df_resultados,
aes(x = generacion, y = fitness)) +
geom_line(aes(group = 1)) +
geom_point() +
labs(title = "Evolución del fitness a lo largo de las generaciones") +
theme_bw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
Paralelización
Uno de los inconvenientes de los algoritmos genéticos es su alto requerimiento computacional. Por ejemplo, si se establecen 1000 generaciones con 150 individuos por generación, la función objetivo se evalúa \(1000x150 = 150000\) veces.
Los ejemplos anteriores se ejecutan en poco tiempo por la sencillez de las funciones objetivo pero, en la práctica, suelen ser mucho más complejas. Dos de las estrategias que se pueden emplear para agilizar el proceso son:
Parada temprana: detener el algoritmo si tras n generaciones consecutivas no se ha conseguido un cambio mínimo. Esta estrategia está implementada en los ejemplos anteriores.
Paralelización:
Evaluar de forma simultánea varios individuos de la población empleando múltiples cores del ordenador.
Paralelizar la función objetivo en el caso de que pueda hacerse.
Versión paralelizada
ESTA IMPLEMENTACIÓN NO FUNCIONA EN WINDOWS
calcular_fitness_poblacion_paral <- function(poblacion, funcion_objetivo,
optimizacion, n_cores = NULL,
verbose = TRUE, ...) {
# Esta función devuelve el fitness de cada individuo de una población.
#
# ARGUMENTOS
# ============================================================================
# poblacion: matriz que representa la población de individuos.
# funcion_objetivo: nombre de la función que se desea optimizar. Debe de haber
# sido definida previamente.
# optimizacion: "maximizar" o "minimizar". Dependiendo de esto, la relación
# del fitness es directamente o indirectamente proporcional
# al valor de la función.
# verbose: mostrar información del proceso por pantalla.
# RETORNO
# ============================================================================
# vector con el fitness de todos los individuos de la población. El orden de
# los valores se corresponde con el orden de las filas de la matriz población.
# Paquetes necesarios para paralelizar.
library("doFuture")
# CÁLCULO DEL FITNESS DE CADA INDIVIDUO DE LA POBLACIÓN
# ----------------------------------------------------------------------------
registerDoFuture()
if (is.null(n_cores)) {
future::plan(
strategy = multiprocess,
workers = future::availableCores() - 1
)
} else {
future::plan(
strategy = multiprocess,
workers = n_cores
)
}
fitness_poblacion <- foreach(i = 1:nrow(poblacion), .combine = c) %dopar% {
calcular_fitness_individuo(
individuo = poblacion[i, ],
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
verbose = verbose
)
}
# MEJOR INDIVIDUO DE LA POBLACIÓN
# ----------------------------------------------------------------------------
# Se identifica el mejor individuo de toda la población, el de mayor
# fitness.
indice_mejor_individuo <- which.max(fitness_poblacion)
# Se identifica el valor de la función objetivo para el mejor individuo.
if (optimizacion == "maximizar") {
valor_funcion <- fitness_poblacion[indice_mejor_individuo]
} else {
valor_funcion <- -1*fitness_poblacion[indice_mejor_individuo]
}
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose) {
cat("------------------", "\n")
cat("Población evaluada", "\n")
cat("------------------", "\n")
cat("Optimización =", optimizacion, "\n")
cat("Mejor fitness encontrado =", fitness_poblacion[indice_mejor_individuo], "\n")
cat("Mejor solución encontrada =",
paste(poblacion[indice_mejor_individuo,], collapse = " "), "\n")
cat("Valor función objetivo =", valor_funcion, "\n")
cat("\n")
}
return(unlist(fitness_poblacion))
}Ejemplo
Se calcula el fitness de todos los individuos de una población formada por 5 individuos.
# Función objetivo a optimizar.
funcion <- function(x, y){
return(x^2 + y^2)
}
# Población simulada.
poblacion <- crear_poblacion(
n_poblacion = 5,
n_variables = 2,
limite_inf = c(-10, -10),
limite_sup = c(+10, +10),
verbose = TRUE)## Población inicial creada
## ------------------------
## Fecha creación: 2020-11-12 12:22:05
## Número de individuos = 5
## Límites inferiores de cada variable = -10, -10
## Límites superiores de cada variable = 10, 10# Cálculo del fitness de todos los individuos.
fitness_poblacion <- calcular_fitness_poblacion_paral(
poblacion = poblacion,
funcion_objetivo = funcion,
optimizacion = "minimizar",
verbose = TRUE
)## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = -2.82793533056974 1.30472700577229
## Fitness = -9.699531
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = -3.6499057803303 2.27837579790503
## Fitness = -18.51281
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = 3.48152958322316 5.6860087485984
## Fitness = -44.45174
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = -4.12048844620585 7.09563705604523
## Fitness = -67.32649
##
## El individuo ha sido evaluado
## -----------------------------
## Optimización = minimizar
## Individuo = -5.56019010487944 8.79520375747234
## Fitness = -108.2713
##
## ------------------
## Población evaluada
## ------------------
## Optimización = minimizar
## Mejor fitness encontrado = -9.699531
## Mejor solución encontrada = -2.82793533056974 1.30472700577229
## Valor función objetivo = 9.699531## [1] -9.699531 -18.512808 -44.451744 -67.326490 -108.271323Para incluir la opción de paralelizado, se repite la función del algoritmo completo, esta vez, incluyendo el argumento paralelizado con el que el usuario pueda especificar que se emplee la función calcular_fitness_poblacion o calcular_fitness_poblacion_paral.
optimizar_ga <- function(
funcion_objetivo,
n_variables,
optimizacion,
limite_inf = NULL,
limite_sup = NULL,
n_poblacion = 20,
n_generaciones = 50,
elitismo = 0.1,
prob_mut = 0.01,
distribucion = "uniforme",
media_distribucion = 1,
sd_distribucion = 1,
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
paralelizado = FALSE,
n_cores = NULL,
verbose = 1,
...) {
# ARGUMENTOS
# =============================================================================
# funcion_objetivo: nombre de la función que se desea optimizar. Debe de haber
# sido definida previamente.
# n_variables: longitud de los individuos.
# optimizacion: "maximizar" o "minimizar". Dependiendo de esto, la relación
# del fitness es directamente o indirectamente proporcional al
# valor de la función.
# limite_inf: vector con el límite inferior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quiere acotar.
# limite_sup: vector con el límite superior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quieren acotar.
# n_poblacion: número total de individuos de la población.
# n_generaciones: número total de generaciones creadas.
# elitismo: porcentaje de mejores individuos de la población actual que
# pasan directamente a la siguiente población.
# prob_mut: probabilidad que tiene cada posición del individuo de mutar.
# distribucion: distribución de la que obtener el factor de mutación. Puede
# ser: "normal", "uniforme" o "aleatoria".
# media_distribucion: media de la distribución si se selecciona distribucion="normal".
# sd_distribucion: desviación estándar de la distribución si se selecciona
# distribucion="normal".
# min_distribucion: mínimo la distribución si se selecciona distribucion="uniforme".
# max_distribucion: máximo la distribución si se selecciona distribucion="uniforme".
# metodo_seleccion: método para establecer la probabilidad de selección. Puede
# ser: "ruleta", "rank" o "tournament".
# metodo_seleccion: método para cruzar los individuos. Puede ser: "uniforme",
# "punto_simple".
# parada_temprana: si durante las últimas "rondas_parada" generaciones la diferencia
# absoluta entre mejores individuos no es superior al valor de
# "tolerancia_parada", se detiene el algoritmo y no se crean
# nuevas generaciones.
# rondas_parada: número de generaciones consecutivas sin mejora mínima para que
# se active la parada temprana.
# tolerancia_parada: valor mínimo que debe tener la diferencia de generaciones
# consecutivas para considerar que hay cambio.
# paralelizado: TRUE para paralelizar el algoritmo genético.
# n_cores: número de cores para la paralelización.
# verbose: Nivel de detalle para que se imprima por pantalla el
# resultado de cada paso del algoritmo (0, 1, 2)
# RETORNO
# =============================================================================
# La función devuelve una lista con 5 elementos:
# fitness: una lista con el fitness del mejor individuo de cada
# generación.
# mejores_individuos: una lista con la combinación de predictores del mejor
# individuo de cada generación.
# mejor_individuo: combinación de predictores del mejor individuo encontrado
# en todo el proceso.
# diferencia_abs: una lista con la diferencia absoluta entre el fitness
# del mejor individuo de generaciones consecutivas.
# df_resultados: un dataframe con todos los resultados anteriores.
start_time <- Sys.time()
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
# Si se activa la parada temprana, hay que especificar los argumentos
# rondas_parada y tolerancia_parada.
if (isTRUE(parada_temprana) &
(is.null(rondas_parada) | is.null(tolerancia_parada)) ) {
stop(paste(
"Para activar la parada temprana es necesario indicar un valor",
"de rondas_parada y de tolerancia_parada."
))
}
# ESTABLECER LOS LÍMITES DE BÚSQUEDA SI EL USUARIO NO LO HA HECHO
# ----------------------------------------------------------------------------
if (is.null(limite_sup) | is.null(limite_inf)) {
warning(paste(
"Es altamente recomendable indicar los límites dentro de los",
"cuales debe buscarse la solución de cada variable.",
"Por defecto se emplea: [-10^3, 10^3]."
))
}
if (any(
is.null(limite_sup), is.null(limite_inf), any(is.na(limite_sup)),
any(is.na(limite_inf))
)) {
warning(paste(
"Los límites empleados por defecto cuando no se han definido son:",
" [-10^3, 10^3]."
))
cat("\n")
}
# Si no se especifica limite_inf, el valor mínimo que pueden tomar las variables
# es -10^3.
if (is.null(limite_inf)) {
limite_inf <- rep(x = -10^3, times = n_variables)
}
# Si no se especifica limite_sup, el valor máximo que pueden tomar las variables
# es 10^3.
if (is.null(limite_sup)) {
limite_sup <- rep(x = 10^3, times = n_variables)
}
# Si los límites no son nulos, se reemplazan aquellas posiciones NA por el valor
# por defecto -10^3 y 10^3.
if (!is.null(limite_inf)) {
limite_inf[is.na(limite_inf)] <- -10^3
}
if (!is.null(limite_sup)) {
limite_sup[is.na(limite_sup)] <- 10^3
}
# ALMACENAMIENTO DE RESULTADOS
# ----------------------------------------------------------------------------
# Por cada generación se almacena, la población, el mejor individuo, su fitness,
# y la diferencia absoluta respecto a la última generación.
poblaciones <- vector(mode = "list", length = n_generaciones)
resultados_fitness <- vector(mode = "list", length = n_generaciones)
resultados_individuo <- vector(mode = "list", length = n_generaciones)
diferencia_abs <- vector(mode = "list", length = n_generaciones)
# ITERACIÓN DE POBLACIONES
# ----------------------------------------------------------------------------
for (i in 1:n_generaciones) {
if (verbose %in% c(1,2)) {
cat("-------------------", "\n")
cat("Generación:", paste0(i, "\\", n_generaciones), "\n")
cat("-------------------", "\n")
}
if (i == 1) {
# CREACIÓN DE LA POBLACIÓN INICIAL
# ------------------------------------------------------------------------
poblacion <- crear_poblacion(
n_poblacion = n_poblacion,
n_variables = n_variables,
limite_inf = limite_inf,
limite_sup = limite_sup,
verbose = verbose %in% c(2)
)
}
# CALCULAR FITNESS DE LOS INDIVIDUOS DE LA POBLACIÓN
# --------------------------------------------------------------------------
if (!paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion(
poblacion = poblacion,
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
verbose = verbose %in% c(2)
)
}
if (paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion_paral(
poblacion = poblacion,
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
n_cores = n_cores,
verbose = verbose %in% c(2)
)
}
# SE ALMACENA LA POBLACIÓN Y SU MEJOR INDIVIDUO
# --------------------------------------------------------------------------
poblaciones[[i]] <- poblacion
fitness_mejor_individuo <- max(fitness_ind_poblacion)
mejor_individuo <- poblacion[which.max(fitness_ind_poblacion), ]
resultados_fitness[[i]] <- fitness_mejor_individuo
resultados_individuo[[i]] <- mejor_individuo
# SE CALCULA LA DIFERENCIA ABSOLUTA RESPECTO A LA GENERACIÓN ANTERIOR
# --------------------------------------------------------------------------
# La diferencia solo puede calcularse a partir de la segunda generación.
if (i > 1) {
diferencia_abs[[i]] <- abs(resultados_fitness[[i - 1]] - resultados_fitness[[i]])
}
# NUEVA POBLACIÓN
# --------------------------------------------------------------------------
nueva_poblacion <- matrix(
data = NA,
nrow = nrow(poblacion),
ncol = ncol(poblacion)
)
# ELITISMO
# --------------------------------------------------------------------------
# El elitismo indica el porcentaje de mejores individuos de la población
# actual que pasan directamente a la siguiente población. De esta forma, se
# asegura que, la siguiente generación, no sea nunca inferior.
if (elitismo > 0) {
n_elitismo <- ceiling(nrow(poblacion) * elitismo)
posicion_n_mejores <- order(fitness_ind_poblacion, decreasing = TRUE)
posicion_n_mejores <- posicion_n_mejores[1:n_elitismo]
nueva_poblacion[1:n_elitismo, ] <- poblacion[posicion_n_mejores, ]
} else {
n_elitismo <- 0
}
# CREACIÓN DE NUEVOS INDIVIDUOS POR CRUCES
# --------------------------------------------------------------------------
for (j in (n_elitismo + 1):nrow(nueva_poblacion)) {
# Seleccionar parentales
indice_parental_1 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metodo_seleccion = metodo_seleccion,
verbose = verbose %in% c(2)
)
indice_parental_2 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metodo_seleccion = metodo_seleccion,
verbose = verbose %in% c(2)
)
parental_1 <- poblacion[indice_parental_1, ]
parental_2 <- poblacion[indice_parental_2, ]
# Cruzar parentales para obtener la descendencia
descendencia <- cruzar_individuos(
parental_1 = parental_1,
parental_2 = parental_2,
metodo_cruce = metodo_cruce,
verbose = verbose %in% c(2)
)
# Mutar la descendencia
descendencia <- mutar_individuo(
individuo = descendencia,
prob_mut = prob_mut,
limite_inf = limite_inf,
limite_sup = limite_sup,
distribucion = distribucion,
media_distribucion = media_distribucion,
sd_distribucion = sd_distribucion,
min_distribucion = min_distribucion,
max_distribucion = max_distribucion,
verbose = verbose %in% c(2)
)
nueva_poblacion[j, ] <- descendencia
}
poblacion <- nueva_poblacion
# CRITERIO DE PARADA
# --------------------------------------------------------------------------
# Si durante las últimas n generaciones, la diferencia absoluta entre mejores
# individuos no es superior al valor de tolerancia_parada, se detiene el
# algoritmo y no se crean nuevas generaciones.
if (parada_temprana && (i > rondas_parada)) {
ultimos_n <- tail(unlist(diferencia_abs), n = rondas_parada)
if (all(ultimos_n < tolerancia_parada)) {
cat(
"Algoritmo detenido en la generacion", i,
"por falta cambio mínimo de", tolerancia_parada,
"durante", rondas_parada,
"generaciones consecutivas.",
"\n"
)
break()
}
}
}
# IDENTIFICACIÓN DEL MEJOR INDIVIDUO DE TODO EL PROCESO
# ----------------------------------------------------------------------------
indice_mejor_individuo_global <- which.max(unlist(resultados_fitness))
mejor_fitness_global <- resultados_fitness[[indice_mejor_individuo_global]]
mejor_individuo_global <- resultados_individuo[[indice_mejor_individuo_global]]
# Se identifica el valor de la función objetivo para el mejor individuo.
if (optimizacion == "maximizar") {
mejor_valor_global <- mejor_fitness_global
} else {
mejor_valor_global <- -1*mejor_fitness_global
}
# RESULTADOS
# ----------------------------------------------------------------------------
# Para crear el dataframe se convierten las listas a vectores del mismo tamaño.
resultados_fitness <- unlist(resultados_fitness)
diferencia_abs <- c(NA, unlist(diferencia_abs))
# Si hay parada temprana, algunas generaciones no se alcanzan: Se eliminan sus
# posiciones de las listas de resultados
resultados_individuo <- resultados_individuo[!sapply(resultados_individuo, is.null)]
poblaciones <- poblaciones[!sapply(poblaciones, is.null)]
# Para poder añadir al dataframe la secuencia variables, se concatenan.
variables <- sapply(
X = resultados_individuo,
FUN = function(x) {
paste(x, collapse = ", ")
}
)
df_resultados <- data.frame(
generacion = seq_along(resultados_fitness),
fitness = resultados_fitness,
predictores = variables,
diferencia_abs = diferencia_abs
)
resultados <- list(
mejor_individuo_global = mejor_individuo_global,
mejor_valor_global = mejor_valor_global,
mejor_fitness_por_generacion = resultados_fitness,
mejor_individuo_por_generacion = resultados_individuo,
diferencia_abs = diferencia_abs,
df_resultados = df_resultados,
poblaciones = poblaciones,
funcion_objetivo = funcion_objetivo
)
end_time <- Sys.time()
# INFORMACIÓN ALMACENADA EN LOS ATRIBUTOS
# ----------------------------------------------------------------------------
attr(resultados, "class") <- "optimizacion_ga"
attr(resultados, 'fecha_creacion') <- end_time
attr(resultados, 'duracion_optimizacion') <- paste(
difftime(end_time, start_time, "secs"),
"secs"
)
attr(resultados, 'optimizacion') <- optimizacion
attr(resultados, 'lim_inf') <- limite_inf
attr(resultados, 'lim_sup') <- limite_sup
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'generaciones') <- i
attr(resultados, 'valor_variables') <- mejor_individuo_global
attr(resultados, 'mejor_fitness') <- mejor_fitness_global
attr(resultados, 'optimo_encontrado') <- mejor_valor_global
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'elitismo') <- elitismo
attr(resultados, 'prob_mut') <- prob_mut
attr(resultados, 'metodo_seleccion') <- metodo_seleccion
attr(resultados, 'metodo_cruce') <- metodo_cruce
attr(resultados, 'parada_temprana') <- parada_temprana
attr(resultados, 'rondas_parada') <- rondas_parada
attr(resultados, 'tolerancia_parada') <- tolerancia_parada
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose %in% c(1,2)) {
cat("-----------------------", "\n")
cat("Optimización finalizada", "\n")
cat("-----------------------", "\n")
cat("Fecha finalización =", as.character(Sys.time()), "\n")
cat("Duración selección = ")
print(difftime(end_time, start_time))
cat("Número generaciones =", i, "\n")
cat("Límite inferior =", paste(limite_inf, collapse = ", "), "\n")
cat("Límite superior =", paste(limite_sup, collapse = ", "), "\n")
cat("Optimización =", optimizacion,"\n")
cat("Óptimo encontrado =", mejor_valor_global,"\n")
cat("Valor variables =", mejor_individuo_global, "\n")
cat("\n")
}
return(resultados)
}
print.optimizacion_ga <- function(obj){
# Función print para objetos optimizacion_ga
cat("----------------------------------------------", "\n")
cat("Resultados optimización por algoritmo genético", "\n")
cat("----------------------------------------------", "\n")
cat("Fecha creación =", attr(obj, 'fecha_creacion'), "\n")
cat("Duración selección = ", attr(obj, 'duracion_optimizacion'), "\n")
cat("Número generaciones =", attr(obj, 'generaciones'), "\n")
cat("Límite inferior =", attr(obj, 'lim_inf'), "\n")
cat("Límite superior =", attr(obj, 'lim_sup'), "\n")
cat("Optimización =", attr(obj, 'optimizacion'), "\n")
cat("Óptimo encontrado =", attr(obj, 'optimo_encontrado'), "\n")
cat("Valor variables =", attr(obj, 'valor_variables'), "\n")
cat("Función objetivo =", "\n")
cat("\n")
print(obj$funcion_objetivo)
}Comparación
Se compara el tiempo necesario para ejecutar la optimización empleando paralelización y sin ella. Como la paralelización afecta al paso en el que se calcula el fitness de todos los individuos de la población, la diferencia de tiempo se notará más cuanto mayor sea el tamaño de la población, pero no se verá afectado por el número de generaciones. Para este ejemplo se emplea n_poblacion = 1000, n_generaciones = 100 y se desactiva la parada temprana para que en ambos casos se ejecuten el mismo número de iteraciones.
Función objetivo
# Función objetivo a optimizar.
funcion <- function(x, y){
sin(y)*exp(1-cos(x))^2 + cos(x)*exp(1-sin(y))^2 + (x-y)^2
}Sin paralelización
library(tictoc)
tic()
resultados <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-10, -5),
limite_sup = c(0, 0),
n_poblacion = 1000,
n_generaciones = 100,
elitismo = 0,
prob_mut = 0.1,
distribucion = "aleatoria",
media_distribucion = 1,
sd_distribucion = 1,
min_distribucion = -1,
max_distribucion = 1,
verbose = 0,
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
paralelizado = FALSE
)
toc()## 8.588 sec elapsedCon paralelización
tic()
resultados <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-10, -5),
limite_sup = c(0, 0),
n_poblacion = 1000,
n_generaciones = 100,
elitismo = 0,
prob_mut = 0.1,
distribucion = "aleatoria",
media_distribucion = 1,
sd_distribucion = 1,
min_distribucion = -1,
max_distribucion = 1,
verbose = 0,
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
paralelizado = TRUE
)
toc()## 253.405 sec elapsedPuede observarse que, para una función tan sencilla de calcular, la latencia incorporada por el proceso de paralelización no compensa.
Algoritmo genético y Nelder-Mead
Los métodos de optimización basados en algoritmo genético son buenos identificando áreas de mínimos globales ya que exploran una gran parte del dominio de la función, sin embargo, una vez que se encuentran en una región local, son poco eficientes aproximándose al mínimo en comparación a otros métodos de optimización local como el de Nelder-Mead Simplex.
En el siguiente ejemplo, se pretende minimizar la función \(f(x_1,x_2) = x_1^2 + x_2^2\). Esta función tiene un único mínimo local \((x_1=0, x_2=0)\) que conincide con el mínimo global, por lo que pueden emplearse tanto métodos de optimización global (algoritmo genético) como métodos de optimización local (Nelder-Mead Simplex). Véase la rapidez y precisión con la que encuentran cada uno la solución.
Función objetivo
# Función objetivo a optimizar. Mínimo global en (0,0).
funcion <- function(x1, x2){
return(x1^2 + x2^2)
}Optimización algoritmo genético
tic()
resultados_ga <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-10, -10),
limite_sup = c(10, 10),
n_poblacion = 150,
n_generaciones = 500,
elitismo = 0.01,
prob_mut = 0.1,
distribucion = "uniforme",
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
parada_temprana = TRUE,
rondas_parada = 10,
tolerancia_parada = 10^-8,
verbose = 0
)## Algoritmo detenido en la generacion 22 por falta cambio mínimo de 1e-08 durante 10 generaciones consecutivas.## 0.293 sec elapsedOptimización Nelder-Mead
tic()
optim(
par = c(10, 10),
fn = function(par){
do.call(funcion, args = as.list(par))
},
method = "Nelder-Mead"
)## $par
## [1] 0.0003754010 0.0005179101
##
## $value
## [1] 4.091568e-07
##
## $counts
## function gradient
## 63 NA
##
## $convergence
## [1] 0
##
## $message
## NULL## 0.005 sec elapsedEl método de optimización local Nelder-Mead converge, en mucho menos tiempo, en un valor más próximo al mínimo real.
Combinación de métodos
Dado que ninguno de los dos métodos, algoritmo genético y Nelder-Mead, necesitan las derivadas de la función objetivo, pueden combinarse para conseguir la capacidad de encontrar áreas de mínimos globales del primero y la rápida convergencia local del segundo.
La siguiente función aplica primero un algoritmo genético y después intenta mejorar el resultado final con el algoritmo Nelder-Mead.
optimizar_ga <- function(
funcion_objetivo,
n_variables,
optimizacion,
limite_inf = NULL,
limite_sup = NULL,
n_poblacion = 20,
n_generaciones = 50,
elitismo = 0.1,
prob_mut = 0.01,
distribucion = "uniforme",
media_distribucion = 1,
sd_distribucion = 1,
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
metodo_cruce = "uniforme",
parada_temprana = FALSE,
rondas_parada = NULL,
tolerancia_parada = NULL,
Nelder_Mead = TRUE,
paralelizado = FALSE,
n_cores = NULL,
verbose = 1,
...) {
# ARGUMENTOS
# =============================================================================
# funcion_objetivo: nombre de la función que se desea optimizar. Debe de haber
# sido definida previamente.
# n_variables: longitud de los individuos.
# optimizacion: "maximizar" o "minimizar". Dependiendo de esto, la relación
# del fitness es directamente o indirectamente proporcional al
# valor de la función.
# limite_inf: vector con el límite inferior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quiere acotar.
# limite_sup: vector con el límite superior de cada variable. Si solo se
# quiere imponer límites a algunas variables, emplear NA para
# las que no se quieren acotar.
# n_poblacion: número total de individuos de la población.
# n_generaciones: número total de generaciones creadas.
# elitismo: porcentaje de mejores individuos de la población actual que
# pasan directamente a la siguiente población.
# prob_mut: probabilidad que tiene cada posición del individuo de mutar.
# distribucion: distribución de la que obtener el factor de mutación. Puede
# ser: "normal", "uniforme" o "aleatoria".
# media_distribucion: media de la distribución si se selecciona distribucion="normal".
# sd_distribucion: desviación estándar de la distribución si se selecciona
# distribucion="normal".
# min_distribucion: mínimo la distribución si se selecciona distribucion="uniforme".
# max_distribucion: máximo la distribución si se selecciona distribucion="uniforme".
# metodo_seleccion: método para establecer la probabilidad de selección. Puede
# ser: "ruleta", "rank" o "tournament".
# metodo_seleccion: método para cruzar los individuos. Puede ser: "uniforme",
# "punto_simple".
# parada_temprana: si durante las últimas "rondas_parada" generaciones la diferencia
# absoluta entre mejores individuos no es superior al valor de
# "tolerancia_parada", se detiene el algoritmo y no se crean
# nuevas generaciones.
# rondas_parada: número de generaciones consecutivas sin mejora mínima para que
# se active la parada temprana.
# tolerancia_parada: valor mínimo que debe tener la diferencia de generaciones
# consecutivas para considerar que hay cambio.
# Nelder_Mead: TRUE para que el mejor individuo devuelto por el algoritmo
# genético se intente mejorar con optimización "Nelder_Mead".
# paralelizado: TRUE para paralelizar el algoritmo genético.
# n_cores: número de cores para la paralelización.
# verbose: Nivel de detalle para que se imprima por pantalla el
# resultado de cada paso del algoritmo (0, 1, 2)
# RETORNO
# =============================================================================
# La función devuelve una lista con 5 elementos:
# fitness: una lista con el fitness del mejor individuo de cada
# generación.
# mejores_individuos: una lista con la combinación de predictores del mejor
# individuo de cada generación.
# mejor_individuo: combinación de predictores del mejor individuo encontrado
# en todo el proceso.
# diferencia_abs: una lista con la diferencia absoluta entre el fitness
# del mejor individuo de generaciones consecutivas.
# df_resultados: un dataframe con todos los resultados anteriores.
start_time <- Sys.time()
# COMPROBACIONES INICIALES
# ----------------------------------------------------------------------------
# Si se activa la parada temprana, hay que especificar los argumentos
# rondas_parada y tolerancia_parada.
if (isTRUE(parada_temprana) &
(is.null(rondas_parada) | is.null(tolerancia_parada)) ) {
stop(paste(
"Para activar la parada temprana es necesario indicar un valor",
"de rondas_parada y de tolerancia_parada."
))
}
# ESTABLECER LOS LÍMITES DE BÚSQUEDA SI EL USUARIO NO LO HA HECHO
# ----------------------------------------------------------------------------
if (is.null(limite_sup) | is.null(limite_inf)) {
warning(paste(
"Es altamente recomendable indicar los límites dentro de los",
"cuales debe buscarse la solución de cada variable.",
"Por defecto se emplea: [-10^3, 10^3]."
))
}
if (any(
is.null(limite_sup), is.null(limite_inf), any(is.na(limite_sup)),
any(is.na(limite_inf))
)) {
warning(paste(
"Los límites empleados por defecto cuando no se han definido son:",
" [-10^3, 10^3]."
))
cat("\n")
}
# Si no se especifica limite_inf, el valor mínimo que pueden tomar las variables
# es -10^3.
if (is.null(limite_inf)) {
limite_inf <- rep(x = -10^3, times = n_variables)
}
# Si no se especifica limite_sup, el valor máximo que pueden tomar las variables
# es 10^3.
if (is.null(limite_sup)) {
limite_sup <- rep(x = 10^3, times = n_variables)
}
# Si los límites no son nulos, se reemplazan aquellas posiciones NA por el valor
# por defecto -10^3 y 10^3.
if (!is.null(limite_inf)) {
limite_inf[is.na(limite_inf)] <- -10^3
}
if (!is.null(limite_sup)) {
limite_sup[is.na(limite_sup)] <- 10^3
}
# ALMACENAMIENTO DE RESULTADOS
# ----------------------------------------------------------------------------
# Por cada generación se almacena, la población, el mejor individuo, su fitness,
# y la diferencia absoluta respecto a la última generación.
poblaciones <- vector(mode = "list", length = n_generaciones)
resultados_fitness <- vector(mode = "list", length = n_generaciones)
resultados_individuo <- vector(mode = "list", length = n_generaciones)
diferencia_abs <- vector(mode = "list", length = n_generaciones)
# ITERACIÓN DE POBLACIONES
# ----------------------------------------------------------------------------
for (i in 1:n_generaciones) {
if (verbose %in% c(1,2)) {
cat("-------------------", "\n")
cat("Generación:", paste0(i, "\\", n_generaciones), "\n")
cat("-------------------", "\n")
}
if (i == 1) {
# CREACIÓN DE LA POBLACIÓN INICIAL
# ------------------------------------------------------------------------
poblacion <- crear_poblacion(
n_poblacion = n_poblacion,
n_variables = n_variables,
limite_inf = limite_inf,
limite_sup = limite_sup,
verbose = verbose %in% c(2)
)
}
# CALCULAR FITNESS DE LOS INDIVIDUOS DE LA POBLACIÓN
# --------------------------------------------------------------------------
if (!paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion(
poblacion = poblacion,
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
verbose = verbose %in% c(2)
)
}
if (paralelizado) {
fitness_ind_poblacion <- calcular_fitness_poblacion_paral(
poblacion = poblacion,
funcion_objetivo = funcion_objetivo,
optimizacion = optimizacion,
n_cores = n_cores,
verbose = verbose %in% c(2)
)
}
# SE ALMACENA LA POBLACIÓN Y SU MEJOR INDIVIDUO
# --------------------------------------------------------------------------
poblaciones[[i]] <- poblacion
fitness_mejor_individuo <- max(fitness_ind_poblacion)
mejor_individuo <- poblacion[which.max(fitness_ind_poblacion), ]
resultados_fitness[[i]] <- fitness_mejor_individuo
resultados_individuo[[i]] <- mejor_individuo
# SE CALCULA LA DIFERENCIA ABSOLUTA RESPECTO A LA GENERACIÓN ANTERIOR
# --------------------------------------------------------------------------
# La diferencia solo puede calcularse a partir de la segunda generación.
if (i > 1) {
diferencia_abs[[i]] <- abs(resultados_fitness[[i - 1]] - resultados_fitness[[i]])
}
# NUEVA POBLACIÓN
# --------------------------------------------------------------------------
nueva_poblacion <- matrix(
data = NA,
nrow = nrow(poblacion),
ncol = ncol(poblacion)
)
# ELITISMO
# --------------------------------------------------------------------------
# El elitismo indica el porcentaje de mejores individuos de la población
# actual que pasan directamente a la siguiente población. De esta forma, se
# asegura que, la siguiente generación, no sea nunca inferior.
if (elitismo > 0) {
n_elitismo <- ceiling(nrow(poblacion) * elitismo)
posicion_n_mejores <- order(fitness_ind_poblacion, decreasing = TRUE)
posicion_n_mejores <- posicion_n_mejores[1:n_elitismo]
nueva_poblacion[1:n_elitismo, ] <- poblacion[posicion_n_mejores, ]
} else {
n_elitismo <- 0
}
# CREACIÓN DE NUEVOS INDIVIDUOS POR CRUCES
# --------------------------------------------------------------------------
for (j in (n_elitismo + 1):nrow(nueva_poblacion)) {
# Seleccionar parentales
indice_parental_1 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metodo_seleccion = metodo_seleccion,
verbose = verbose %in% c(2)
)
indice_parental_2 <- seleccionar_individuo(
vector_fitness = fitness_ind_poblacion,
metodo_seleccion = metodo_seleccion,
verbose = verbose %in% c(2)
)
parental_1 <- poblacion[indice_parental_1, ]
parental_2 <- poblacion[indice_parental_2, ]
# Cruzar parentales para obtener la descendencia
descendencia <- cruzar_individuos(
parental_1 = parental_1,
parental_2 = parental_2,
metodo_cruce = metodo_cruce,
verbose = verbose %in% c(2)
)
# Mutar la descendencia
descendencia <- mutar_individuo(
individuo = descendencia,
prob_mut = prob_mut,
limite_inf = limite_inf,
limite_sup = limite_sup,
distribucion = distribucion,
media_distribucion = media_distribucion,
sd_distribucion = sd_distribucion,
min_distribucion = min_distribucion,
max_distribucion = max_distribucion,
verbose = verbose %in% c(2)
)
nueva_poblacion[j, ] <- descendencia
}
poblacion <- nueva_poblacion
# CRITERIO DE PARADA
# --------------------------------------------------------------------------
# Si durante las últimas n generaciones, la diferencia absoluta entre mejores
# individuos no es superior al valor de tolerancia_parada, se detiene el
# algoritmo y no se crean nuevas generaciones.
if (parada_temprana && (i > rondas_parada)) {
ultimos_n <- tail(unlist(diferencia_abs), n = rondas_parada)
if (all(ultimos_n < tolerancia_parada)) {
cat(
"Algoritmo detenido en la generacion", i,
"por falta cambio mínimo de", tolerancia_parada,
"durante", rondas_parada,
"generaciones consecutivas.",
"\n"
)
break()
}
}
}
# IDENTIFICACIÓN DEL MEJOR INDIVIDUO DE TODO EL PROCESO
# ----------------------------------------------------------------------------
indice_mejor_individuo_global <- which.max(unlist(resultados_fitness))
mejor_fitness_global <- resultados_fitness[[indice_mejor_individuo_global]]
mejor_individuo_global <- resultados_individuo[[indice_mejor_individuo_global]]
# Si Nelder_Mead = TRUE, el mejor individuo identificado mediante el
# algoritmo genético se somete a un proceso de optimización Nelder-Mead con
# el objetivo de mejorar la convergencia final.
if (Nelder_Mead) {
optmizacion_NM <- optim(
par = mejor_individuo_global,
fn = function(par) {
if (optimizacion == "maximizar") {
do.call(-funcion, args = as.list(par))
}else{
do.call(funcion, args = as.list(par))
}
},
method = "Nelder-Mead"
)
mejor_individuo_global <- optmizacion_NM$par
mejor_fitness_global <- optmizacion_NM$value
}
# Se identifica el valor de la función objetivo para el mejor individuo.
if (optimizacion == "maximizar") {
mejor_valor_global <- mejor_fitness_global
} else {
mejor_valor_global <- -1*mejor_fitness_global
}
# RESULTADOS
# ----------------------------------------------------------------------------
# Para crear el dataframe se convierten las listas a vectores del mismo tamaño.
resultados_fitness <- unlist(resultados_fitness)
diferencia_abs <- c(NA, unlist(diferencia_abs))
# Si hay parada temprana, algunas generaciones no se alcanzan: Se eliminan sus
# posiciones de las listas de resultados
resultados_individuo <- resultados_individuo[!sapply(resultados_individuo, is.null)]
poblaciones <- poblaciones[!sapply(poblaciones, is.null)]
# Para poder añadir al dataframe la secuencia variables, se concatenan.
variables <- sapply(
X = resultados_individuo,
FUN = function(x) {
paste(x, collapse = ", ")
}
)
df_resultados <- data.frame(
generacion = seq_along(resultados_fitness),
fitness = resultados_fitness,
predictores = variables,
diferencia_abs = diferencia_abs
)
resultados <- list(
mejor_individuo_global = mejor_individuo_global,
mejor_valor_global = mejor_valor_global,
mejor_fitness_por_generacion = resultados_fitness,
mejor_individuo_por_generacion = resultados_individuo,
diferencia_abs = diferencia_abs,
df_resultados = df_resultados,
poblaciones = poblaciones,
funcion_objetivo = funcion_objetivo
)
end_time <- Sys.time()
# INFORMACIÓN ALMACENADA EN LOS ATRIBUTOS
# ----------------------------------------------------------------------------
attr(resultados, "class") <- "optimizacion_ga"
attr(resultados, 'fecha_creacion') <- end_time
attr(resultados, 'duracion_optimizacion') <- paste(
difftime(end_time, start_time, "secs"),
"secs"
)
attr(resultados, 'optimizacion') <- optimizacion
attr(resultados, 'Nelder_Mead') <- Nelder_Mead
attr(resultados, 'lim_inf') <- limite_inf
attr(resultados, 'lim_sup') <- limite_sup
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'generaciones') <- i
attr(resultados, 'valor_variables') <- mejor_individuo_global
attr(resultados, 'mejor_fitness') <- mejor_fitness_global
attr(resultados, 'optimo_encontrado') <- mejor_valor_global
attr(resultados, 'n_poblacion') <- n_poblacion
attr(resultados, 'elitismo') <- elitismo
attr(resultados, 'prob_mut') <- prob_mut
attr(resultados, 'metodo_seleccion') <- metodo_seleccion
attr(resultados, 'metodo_cruce') <- metodo_cruce
attr(resultados, 'parada_temprana') <- parada_temprana
attr(resultados, 'rondas_parada') <- rondas_parada
attr(resultados, 'tolerancia_parada') <- tolerancia_parada
# INFORMACIÓN DEL PROCESO (VERBOSE)
# ----------------------------------------------------------------------------
if (verbose %in% c(1,2)) {
cat("-----------------------", "\n")
cat("Optimización finalizada", "\n")
cat("-----------------------", "\n")
cat("Fecha finalización =", as.character(Sys.time()), "\n")
cat("Duración selección = ")
print(difftime(end_time, start_time))
cat("Número generaciones =", i, "\n")
cat("Límite inferior =", paste(limite_inf, collapse = ", "), "\n")
cat("Límite superior =", paste(limite_sup, collapse = ", "), "\n")
cat("Optimización =", optimizacion,"\n")
cat("Nelder_Mead =", Nelder_Mead,"\n")
cat("Óptimo encontrado =", mejor_valor_global,"\n")
cat("Valor variables =", mejor_individuo_global, "\n")
cat("\n")
}
return(resultados)
}
print.optimizacion_ga <- function(obj){
# Función print para objetos optimizacion_ga
cat("----------------------------------------------", "\n")
cat("Resultados optimización por algoritmo genético", "\n")
cat("----------------------------------------------", "\n")
cat("Fecha creación =", attr(obj, 'fecha_creacion'), "\n")
cat("Duración selección =", attr(obj, 'duracion_optimizacion'), "\n")
cat("Número generaciones =", attr(obj, 'generaciones'), "\n")
cat("Límite inferior =", attr(obj, 'lim_inf'), "\n")
cat("Límite superior =", attr(obj, 'lim_sup'), "\n")
cat("Optimización =", attr(obj, 'optimizacion'), "\n")
cat("Nelder_Mead =", attr(obj, 'Nelder_Mead'), "\n")
cat("Óptimo encontrado =", attr(obj, 'optimo_encontrado'), "\n")
cat("Valor variables =", attr(obj, 'valor_variables'), "\n")
cat("\n")
cat("Función objetivo :", "\n")
cat("\n")
print(obj$funcion_objetivo)
}Comparación
Función objetivo
# Función objetivo a optimizar. Mínimo global en (0,0).
funcion <- function(x1, x2){
-20*exp(-0.7*sqrt(0.5*(x1^2 + x2^2))) - exp(0.5*(cos(2*pi*x1) + cos(2*pi*x2))) + exp(1) + 20
}Sin Nelder Mead
resultados_ga <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-10, -10),
limite_sup = c(10, 10),
n_poblacion = 150,
n_generaciones = 500,
elitismo = 0.01,
prob_mut = 0.1,
distribucion = "uniforme",
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
parada_temprana = TRUE,
rondas_parada = 10,
tolerancia_parada = 10^-8,
Nelder_Mead = FALSE,
paralelizado = FALSE,
verbose = 0
)## Algoritmo detenido en la generacion 28 por falta cambio mínimo de 1e-08 durante 10 generaciones consecutivas.## [1] -0.0109890504 -0.0005063876Con Nelder Mead
resultados_ga <- optimizar_ga(
funcion_objetivo = funcion,
n_variables = 2,
optimizacion = "minimizar",
limite_inf = c(-10, -10),
limite_sup = c(10, 10),
n_poblacion = 150,
n_generaciones = 500,
elitismo = 0.01,
prob_mut = 0.1,
distribucion = "uniforme",
min_distribucion = -1,
max_distribucion = 1,
metodo_seleccion = "tournament",
parada_temprana = TRUE,
rondas_parada = 10,
tolerancia_parada = 10^-8,
Nelder_Mead = TRUE,
paralelizado = FALSE,
verbose = 0
)## Algoritmo detenido en la generacion 22 por falta cambio mínimo de 1e-08 durante 10 generaciones consecutivas.## [1] 8.097641e-11 6.811295e-11Empleando la combinación de ambos métodos de optimización se obtiene un resultado mucho más próximo al mínimo global comparado a si solo se emplea el algoritmo genético.
Información sesión
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)Bibliografía
John McCall, Genetic algorithms for modelling and optimisation, Journal of Computational and Applied Mathematics, Volume 184, Issue 1, 2005
Optimizing with GeneticAlgorithmsby Benjamin J. Lynch Feb 23, 2006
Durand, Nicolas & Alliot, Jean-Marc. (1999). A Combined Nelder-Mead Simplex and Genetic Algorithm.
Abdel-Rahman Hedar & Masao Fukushima (2003) Minimizing multimodal functions by simplex coding genetic algorithm, Optimization Methods and Software.
https://en.wikipedia.org/wiki/Genetic_algorithm
https://en.wikipedia.org/wiki/Crossover_(genetic_algorithm)
https://www.sfu.ca/~ssurjano/optimization.html
¿Cómo citar este documento?
Optimización con algoritmo genético y Nelder-Mead por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/48_optimizacion_con_algoritmo_genetico
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.