Optimización univariada: Método de fase de acotamiento en R

El método de fase de acotamiento es usado para "encerrar" el mínimo de una función. Este método garantiza encerrar el mínimo de una función unimodal. El algoritmo comienza con una apuesta inicial y partir de ahí encuentra una dirección de búsqueda al efectuar dos evaluaciones más en la vecindad del punto inicial. Posteriormente, una estrategia de búsqueda exponencial es adoptada to alcanzar el óptimo. Asumiendo que buscamos el mínimo, a continuación se presenta el algoritmo:

Algoritmo 

Paso 1: Elige un punto inicial x⁽⁰⁾ y un incremento ∆. Asigna k=0.
Paso 2: Si f(x⁽⁰⁾ − |∆|) ≥ f(x⁽⁰⁾ ) ≥ f(x⁽⁰⁾ + |∆|), luego ∆ es positivo;
si no, si f(x⁽⁰⁾ − |∆|) ≤ f(x⁽⁰⁾ ) ≤ f(x⁽⁰⁾ + |∆|), luego ∆ es negativo;

si no, ir al paso 1.

Paso 3: Asigna x^(k+1) = x^(k) + 2^k ∆.
Paso 4: Si f(x^(k+1) ) < f(x^(k)), asigna k = k + 1 y ve al paso 3;
si no, el mínimo está en el intervalo (x^(k−1) , x^(k+1)) y termina.

Si el valor de ∆ es grande, la precisión del intervalo al punto mínimo es baja pero es encontrado rápidamente. Y por el contrario, si el valor de ∆ es pequeño, la precisión del intervalo es mejor, pero más evaluaciones de la función pueden ser necesarias para encerrar al mínimo. Se presenta un ejemplo para encontrar el intervalo que encierra al mínimo de la función:

f(x) = x² + 54/x

 Ahora presentamos la implementación del algoritmo en R.
fase_acotamiento <- function(x0 = 0.6, delta = 0.5, evalua = function(x){x^2 + 54/x}, trace = FALSE){
  solucion <- NULL
  x <- x0
 
  #  Paso 1:
  k = 0

  # Paso 2:
  fx_minus_delta = evalua(x0 - abs(delta))
  fx0 = evalua(x0)
  fx_plus_delta = evalua(x0 + abs(delta))
  if (fx_minus_delta  >= fx0  & fx0 >= fx_plus_delta){
    delta = abs(delta)
  }else if(fx_minus_delta  <= fx0  & fx0 <= fx_plus_delta){
    delta = abs(delta) * -1
  } else{
    stop("Elige un nuevo punto inicial")
  }

  while(is.null(solucion)){
    # Paso 3:
    assign( paste0("x", k+1), get(paste0("x",k)) + 2^k*delta)

    # Paso 4:
    if (evalua( get(paste0("x",k+1)) ) < evalua( get(paste0("x",k))) ){
      k = k + 1
      if (trace) cat("Iteracion:", k, "\tx: ", get(paste0("x",k)) , "\tf(x)", evalua(get(paste0("x",k))) ,"\n")
    }else{
      solucion <- c( get(paste0("x", k-1)),   get(paste0("x", k+1))  )
    }
  }
  if (trace)  cat("\nEl mínimo se encuentra en el intervalo:\n")
  sort(solucion)
}

En el siguiente ejemplo, aplicamos nuestro algoritmo sobre la función presentada, con un valor de incremento ∆ relativamente pequeño, así como un punto de inicio x⁽⁰⁾ en 0.6.

> # Ejemplo
> fase_acotamiento(x0=0.6, delta = 0.000001, trace = TRUE)
Iteracion: 1     x:  0.600001     f(x) 90.35985
Iteracion: 2     x:  0.600003     f(x) 90.35955
Iteracion: 3     x:  0.600007     f(x) 90.35896
Iteracion: 4     x:  0.600015     f(x) 90.35777
Iteracion: 5     x:  0.600031     f(x) 90.35539
Iteracion: 6     x:  0.600063     f(x) 90.35063
Iteracion: 7     x:  0.600127     f(x) 90.34111
Iteracion: 8     x:  0.600255     f(x) 90.32207
Iteracion: 9     x:  0.600511     f(x) 90.28403
Iteracion: 10     x:  0.601023     f(x) 90.20804
Iteracion: 11     x:  0.602047     f(x) 90.05645
Iteracion: 12     x:  0.604095     f(x) 89.75484
Iteracion: 13     x:  0.608191     f(x) 89.15779
Iteracion: 14     x:  0.616383     f(x) 87.9878
Iteracion: 15     x:  0.632767     f(x) 85.73986
Iteracion: 16     x:  0.665535     f(x) 81.58067
Iteracion: 17     x:  0.731071     f(x) 74.3987
Iteracion: 18     x:  0.862143     f(x) 63.37791
Iteracion: 19     x:  1.124287     f(x) 49.29446
Iteracion: 20     x:  1.648575     f(x) 35.47336
Iteracion: 21     x:  2.697151     f(x) 27.29575

El mínimo se encuentra en el intervalo:
[1] 1.648575 4.794303

Optimización univariada: Método de búsqueda exhaustiva en R

Se implementa el método de búsqueda exhaustiva siguiendo el algoritmo presentado en "Optimization for Engineering Design" de Kalyanmoy Deb. Es un método de optimización univariada, en donde se hace el supuesto de que la función a optimizar cumple los requerimientos de concavidad/convexidad, unimodal y continua. Se trata del método más simple entre el resto de los métodos de optimización univariada, por lo que es un buen punto de partida.

En el método de búsqueda exhaustiva, el óptimo de la función es encerrado al calcular los valores de la función de puntos igualmente espaciados (vea figura). Usualmente, la búsqueda comienza desde un límite inferior de la variable y tres valores de la función consecutivos son comparados al mismo tiempo basado en el supuesto de unimodalidad de la función. Basado en el resultado de la comparación, la búsqueda se termina o continúa reemplazando uno de los tres puntos por un  nuevo punto. La búsqueda continua hasta que el mínimo es encerrado.

Fig. Método de búsqueda exhaustiva usa tres puntos igualmente espaciados

Implementación del algoritmo

exhaustiva <- function(a, b, n=10, evalua = function(x){x^2 + 54/x}){
  ## Parámetros
  # a: Límite inferior del intervalo.
  # b: Límite superior del intervalo.
  # n: Número de intervalos (default 10).
  # evalua: Función objetivo (default x² + 54/x).

  # Almacenar solución
  solucion <- c()
  # Contador
  iter <- 1

  # Paso 1
  x1 = a
  delta_x <- (b-a)/n
  x2 = x1 + delta_x
  x3 = x2 + delta_x

  # Paso 3a 
  while(x3 <= b  & length(solucion) == 0){
    fx1 = evalua(x1)
    fx2 = evalua(x2)
    fx3 = evalua(x3)
    # Paso 2
    if( fx1  >= fx2  & fx2 <= fx3 ){
      cat("\nFinalizado en ",iter,"Iteraciones\n")
      solucion <- c(x1,x3)
      return(list(Solucion=solucion))
    }else{
      cat("Iteración:",iter,"\t  x1:",x1,"\t  x2:",x2,"\t  x3:",x3,"\n")
      iter <- iter + 1

      # Actualiza valores
      x1 = x2
      x2 = x3
      x3 = x2 + delta_x
    } # Fin paso 2
  } # Fin paso 3a
  # Paso 3b
  cat("\nFinalizado en ",iter,"Iteraciones: Sin solución\n")
  return(list(Rango =  c(a,b)))
}

# Ejemplo
exhaustiva(0,5,500,  function(x){x^2 + 54/x})

> exhaustiva(0,5,25,  function(x){x^2 + 54/x})
Iteración: 1       x1: 0.0       x2: 0.2       x3: 0.4
Iteración: 2       x1: 0.2       x2: 0.4       x3: 0.6
Iteración: 3       x1: 0.4       x2: 0.6       x3: 0.8
Iteración: 4       x1: 0.6       x2: 0.8       x3: 1.0
Iteración: 5       x1: 0.8       x2: 1.0       x3: 1.2
Iteración: 6       x1: 1.0       x2: 1.2       x3: 1.4
Iteración: 7       x1: 1.2       x2: 1.4       x3: 1.6
Iteración: 8       x1: 1.4       x2: 1.6       x3: 1.8
Iteración: 9       x1: 1.6       x2: 1.8       x3: 2.0
Iteración: 10       x1: 1.8       x2: 2.0       x3: 2.2
Iteración: 11       x1: 2.0       x2: 2.2       x3: 2.4
Iteración: 12       x1: 2.2       x2: 2.4       x3: 2.6
Iteración: 13       x1: 2.4       x2: 2.6       x3: 2.8
Iteración: 14       x1: 2.6       x2: 2.8       x3: 3.0

Finalizado en  15 Iteraciones
$Solucion
[1] 2.8 3.2

>

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

Algoritmo de Dijkstra en R

Con el algoritmo de Dijkstra es posible determinar la distancia más corta (de menos esfuerzo o costo) entre un nodo inicial  y cualquier otro nodo en un grafo. La idea del algoritmo es calcular continuamente la distancia más corta desde un punto inicial y excluyendo las distancias más largas cuando se efectúa una actualización. Como ejemplo, en el siguiente grid, se asume que los nodos en morado son sumamente costosos, por lo que el algoritmo de Dijkstra deberá hallar la ruta que minimice el costo total entre cualquier tupla de nodos, como por ejemplo, pasar del nodo "A" al nodo "W".

¿Cómo instalar Ocaml para Windows con MSYS2?

En primer lugar instala MSYS2, para lo cual puedes seguir las instrucciones de instalación que hay en la página del programa http://www.msys2.org/, las intrucciones son muy claras.

Posteriormente confirmaremos si el paquete existe con el commando

$ pacman -Ss ocaml

Lo que imprimirá en pantalla:

Como puedes obervar, en el caso de OCaml, existen dos repositorios mingw32 y mingw64, los sufijos determinan la versión de 32 o 64bits respectivamente. De igual forma el prefijo mingw-w64-i686- es para la versión de 32-bits y; mingw-w64-x86_64- es para la versión de 64-bits.

Nota además que en este ejemplo, la versión 4.04.0-1 ya se encuentra instalada.

Una vez elegida la versión (32/64-bits), procedemos a ejecutar el siguiente comando:

pacman -S mingw-w64-i686_64-ocaml

para la versión de 64-bits o:

pacman -S mingw-w64-x86_64-ocaml

para la versión de 32-bits. La traza final de la instalación es la siguiente:

Y ahora verificamos que OCaml efectivamente esta instalado.