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.
