Algoritmo de evolución diferencial JADE

El estado del arte de los algoritmos bio-inspirados para problemas complejos de optimización se ha desarrollado con celeridad en los últimos 20 años y con gran éxito. Desde la aparición en 1997 del algoritmo de evolución diferencial con una estrategia de mutación aleatoria y sin sesgo  (DE/rand/1/bin) para la optimización en espacios continuos (R. Storn y K. Price), han surgido diversas mejoras a este algoritmo que brindan ventajas y desventajas dependiendo de la función a optimizar.

En esta entrada se implementa el algoritmo JADE (Jingqiao Zhang y Arthur C. Sanderson) sin archivo (tamaño de archivo igual a cero), el cual surgió desde 2007 pero la literatura lo reporta formalmente a partir de 2009. Este algoritmo propone una estrategia de mutación con sesgo hacia las p mejores soluciones en la población (DE/current-to-pBest) y con un archivo A opcional que ayuda a mantener cierto nivel de diversidad en la población. El efecto sobre la búsqueda, es una convergencia más rápida pero que puede estancarse en óptimos locales. Otra de las ventajas importantes de este algoritmo, es que no requiere de una calibración fina de los parámetros de mutación F y cruza CR, ya que estos parámetros se adaptan mediante retroalimentación de la búsqueda evolutiva generación tras generación.

Como benchmark, empleamos una función unimodal (esfera) y una función multimodal (rastrigin). Dejamos los detalles de la implementación como tarea para el lector, en el artículo de Zhang y Sanderson.

#-------------------------------------------------------------------
### Para adaptación de parámetros F y CR (vea algortimo JADE, 2009)
### Distribución Normal
randN = function(Q, mCR){
CR = rnorm(Q, mCR, 0.05)
if (CR < 0 ){
CR = 0
} else if( CR > 1){
CR = 1
}
CR
}

### Distribución Cauchy
randC = function(Q, mF){
F = rcauchy(Q, mF, 0.05)
if (F > 0 & F <= 1){
F
} else if(F > 1){
F = 1
F
} else if(F <= 0){
randC(Q, mF)
}
}

### Media de Lehmer 
meanL = function(sF){
sum(sF^2) / sum(sF)

}

### Obtener los mejores
getPBest = function(x, n=30) {
  which(x >= -sort(-x, partial=n)[n])

}

#-------------------------------------------------------------------
### JADE
de <- function(fobj, limits, NP=40, G=500, c = 0.1, pB){
  # PARÁMETROS DE LA FUNCIÓN
  # fobj: Función objetivo
  # NP: (Entero positivo) Población
  # G: (entero positivo) Generaciones 
  # c: (numerico) Parámetro de adaptación (se recomienda entre 0.05 y 0.2)
  # pB: (float) Parámetro de convergencia (pBest) recomendado entre 0.05 y 0.2

  # INICIA EVOLUCIÓN DIFERENCIAL
  dim  = length(limits)
  pop  = replicate(n, runif(NP)) 
  Lmin = sapply(limits, min)
  Lmax = sapply(limits, max)
  diff = abs(Lmin - Lmax)
  pop_denorm = Lmin + pop * diff # escalamiento
  aptitud = apply(pop_denorm, 1, fobj) # resultado de fobj sobre los parámetros
  
  best_idx = which.min(aptitud) # indice de la mejor solucion
  best = pop_denorm[best_idx,] # best set of scaled parameters
  mCR = 0.68; mF = 0.5 # Default muCR y muF
  pBest = NP * pB

  for(i in seq_len(G)){
    SF = c(); SCR = c() # Inicializa SF y SCR
    for(j in seq_len(NP)){
      CR = randN(1,mCR); # CR value
      F = randC(1, mF) # F value

### current-to-pbest
      idxs = seq_len(NP)[-j]
      idbs = getPBest(aptitud, pBest)   # Indices del Top p%
      R = pop[sample(idxs,2),] # selecciona 2 vectores
      xi = pop[j,] # xi
      xp = pop[sample(idbs,1),] # Selecciona al pbest
      r1 = R[1,]; r2 = R[2,]
      mutant = xi + F * (xp - xi) + F * (r1 - r2) # vector mutante

### /bin
      cross_points = runif(dim) > CR
      if(!all(cross_points)){ # si todo es FALSE, cambia 1 (random) a TRUE
        cross_points[sample(1:dim,1)] = TRUE
      }
      trial = pop[j,]; 
trial[cross_points] = mutant[cross_points] # Cross Over

### Actualización
      f = fobj(trial)
      if(f < aptitud[j]){
        aptitud[j] = f
        pop[j,] = trial # Poblacion
  SCR = c(SCR, CR); SF = c(SF, F) # Exitoso CR y F
        if(f < aptitud[best_idx]){
          best_idx = j
          best = trial
        }
      }
    } # j loop

    #Actualiza valores de mCR y mF
    if (length(SCR) > 0){
       mCR = ((1 - c)* mCR) + (c * mean(SCR))
    }else{
       mCR = 0.9
    }
    if (length(SF)> 0 ){
       mF = ((1 - c)* mF) + (c * meanL(SF))
    }else{
mF = 0.9
    }
cat(i,":",aptitud[best_idx],"\n")
  } # i loop
  return(list(
Mejor=best, # Argumento mínimo
AptMejor=aptitud[best_idx], # Aptitud mínima
pobFinal= pop_denorm, # Población final (PF)
AptPobFinal = aptitud) # Aptitud PF
  )
}

###################################################################
###################################################################
### funcion objetivo: SQUARE
square = function(x) sum(x^2)

### funcion objetivo: RASTRIGIN
rastrigin = function(x) 10*length(x)+sum(x^2-10*cos(2*pi*x))


### LIMITES
n = 10 # Cantidad de variables (D)
limits = rep(list(c(-10.0,10.0)), n)
lims_R = rep(list(c(-5.12,5.12)), n)
options(digits=10)

jades = de(square,    limits = limits,  NP = 40, G = 500, c = 0.1, pB = 0.05)
jader = de(rastrigin, limits = lims_R,  NP = 40, G = 500, c = 0.2, pB = 0.1)

Con sólo una corrida por función (lo que no constituye una estadística robusta sobre el desempeño del algoritmo), se presentan los resultados correspondientes al valor de la función de aptitud de la población final, de la cual cabe destacar que es pequeña y contribuye en poco a la diversidad que se requiere en el proceso de búsqueda (vea muJADE, un algoritmo para la optimización de hasta 1000 variables con poblaciones iniciales pequeñas).

Resultados de la función unimodal:

> jades[[4]]
 [1] 5.370292699e-09 1.087864207e-08 9.548139344e-10 2.836403831e-08
 [5] 1.494881426e-08 2.213226281e-08 1.241186240e-09 3.304930492e-08
 [9] 9.882426906e-08 7.283155956e-09 1.198585980e-09 1.305457928e-07
[13] 1.312062918e-08 9.002579032e-10 4.393782940e-09 2.661652047e-09
[17] 9.120911136e-08 5.828017147e-08 9.263151651e-09 1.768654408e-08
[21] 3.809189247e-08 5.291591193e-09 5.800521692e-08 5.163245683e-08
[25] 2.546040478e-09 2.166963844e-08 4.291270186e-08 1.550838431e-08
[29] 9.134545758e-08 2.998656884e-08 4.619834347e-08 1.248160779e-10
[33] 1.706746591e-08 4.159098415e-08 1.272015946e-07 2.467881668e-11

[37] 2.841199376e-09 5.134834451e-08 3.172126021e-08 7.270710044e-08

Resultados de rastrigin:

> jader[[4]]

 [1] 4.201688917e-05 5.609602482e-06 6.244809170e-05 1.053176746e-05

 [5] 1.495947538e-05 5.334819079e-06 2.191708933e-05 2.692694224e-06

 [9] 1.629059085e-06 7.826923863e-06 7.189288851e-06 5.646302782e-06

[13] 2.552616789e-05 2.349400091e-06 1.470918784e-05 5.497413767e-06

[17] 9.630943651e-07 3.014301072e-05 1.049271211e-05 8.983178503e-07

[21] 1.970211811e-05 4.115299404e-06 1.622325030e-06 2.539362356e-07

[25] 1.478155715e-04 4.101470950e-07 1.130157293e-05 1.634984329e-06

[29] 4.930700320e-05 2.059940211e-05 9.882458585e-08 4.562147637e-08

[33] 2.366632685e-07 4.540429686e-06 4.706199789e-06 6.752678416e-07

[37] 3.458378472e-07 1.215024881e-06 4.130576241e-05 2.420810574e-08

Sin la evidencia que nos daría la ejecución de repeticiones del experimento, y/o la comparación contra otro algoritmo de optimización, es díficil emitir un veredicto sobre el desempeño real de JADE. Por lo que, en otra entrada del blog, se mostrará un análisis estadístico detallado sobre el performance de JADE.

Plot 3D de la función rastrigin en R

La función rastrigin es una función no convexa que regularmente es empleada como test de desempeño en diversos algoritmos de optimización (mono-objetivo). Dado la gran cantidad de mínimos locales (función multimodal), algunos algoritmos de optimización pueden estancarse en el proceso de búsqueda y no convergen al mínimo global.

En esta entrada, mostramos como plotear esta función en el lenguaje R:

# Límite del valor de las variables
x1 = seq(-5.12, 5.12, length = 500)
x2 = x1

# Rastrigin
f = function(x1,x2) { 20+(x1^2-10*cos(2*3.1416*x1))+(x2^2-10*cos(2*3.1416*x2)) }

# Evaluación
z1 = outer(x1,x2,f)
z1[is.na(z1)] = 1

# Paleta de colores (200 colores en gradiente)
# Instala viridis (si aún no lo está)

colors = viridis::viridis(200)

# Tamaño de las regiones

z.facet.center = (z1[-1, -1] + z1[-1, -ncol(z1)] + z1[-nrow(z1), -1] + z1[-nrow(z1), -ncol(z1)])/4

# Rango de la región central en escala 200 (número de colores)

z.facet.range = cut(z.facet.center, 200)

# Comienza el ploteo 3D

persp(x1,x2,z1, theta=30, phi=20, expand=1.0, col=colors[z.facet.range], border=NA, scale=TRUE, ltheta=90, shade=0.1, ticktype="detailed", d=5, r=sqrt(10), nticks=5, box=TRUE, axes=TRUE, xlab="x1", ylab="x2", zlab="f(x1,x2)")

Random pic con NetLogo

NetLogo es un lenguaje de programación basado en agentes que se usa principalmente para la simulación de diversos fenómenos de distintas áreas, incluidas por su puesto las sociales. En los últimos años, y con los avances de la computación, se ha desarrollado un área nueva dentro de la economía denominada Economía Computacional Basada en Agentes (término acuñado por Leigh Tesfatsion), con la cual se pretende demostrar los supuestos de la teoría económica suprimiendo los ajustes impulsados exógenamente (como lo hacen los modelos basados en ecuaciones y que imponen un comportamiento sobre el sistema).

En esta entrada, y como una forma de motivar la incursión en este lenguaje de programación, se presenta el resultado gráfico de una simulación en la cual participan 200 agentes, los cuales se desplazan sobre un grid de manera aleatoria trazando una ruta.

Como de costumbre, se presenta el código utilizado, el cual es muy expresivo.

to setup
  clear-all
  setup-turtles
  setup-patches
  reset-ticks
end

to setup-turtles
  create-turtles 200 [ setxy random-xcor random-ycor ]
  ask turtles [pen-down]
end

to setup-patches
  ask patches [ set pcolor white ]
end

 to go
  move-turtles
  tick
end

to move-turtles
  ask turtles [ right random 30 forward 1 ]
end