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

>