jueves, 6 de septiembre de 2018

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.

No hay comentarios.:

Publicar un comentario