Explica cómo maximizar la función:

f <- function (x) sin(10/(.03-x)) + 1e4*(-x^2+0.03*x)

Es una parábola convexa hacia arriba más un seno, por lo que el máximo estará cerca de \(-b/2a = -0.03/(-2) = -0.015\), o sea, cerca del cero:

plot (f, -5, 5, type="l")
plot (f, -1, 1, type="l")
plot (f, -.1, .1, type="l")
Warning message:
In sin(10/(0.03 - x)) : Se han producido NaNs
plot (f, -.01, .01, type="l")
plot (f, 0, .02, type="l")
plot (f, 0, .02, type="l", n=1001)
plot (f, 0.014, .016, type="l", n=1001)
plot (f, 0.014, .016, type="l", n=1001, ylim=c(3,3.5))
plot (f, 0.014, .016, type="l", n=1001, ylim=c(3.2,3.5))
plot (f, 0.014, .016, type="l", n=1001, ylim=c(3.23,3.26))
plot (f, 0.014, .016, type="l", n=1001, ylim=c(3.24,3.255))
plot (f, 0.014, .016, type="l", n=10001, ylim=c(3.24,3.255))
plot (f, 0.0147, .0155, type="l", n=10001, ylim=c(3.24,3.255))
plot (f, 0.015, .0151, type="l", n=10001, ylim=c(3.24,3.255))
## cerca de 0.01502

Vamos a buscar cerca:

> optimize (f, c(0.015, 0.0151), maximum=TRUE)
$maximum
[1] 0.0150382

$objective
[1] 2.960766

> ## mal, demasiado lejos

> optimize (f, c(0.01501, 0.01503), maximum=TRUE)
$maximum
[1] 0.01501764

$objective
[1] 3.240644

> ## esa parece buena, pero
> abline(v=0.01501764,col=2) # no vale

> optim (0.01502, f, control=list(fnscale=-1)) $ par
[1] 0.01502071
Warning message:
In optim(0.01502, f, control = list(fnscale = -1)) :
  one-dimensional optimization by Nelder-Mead is unreliable:
use "Brent" or optimize() directly
> optim (0.01502, f, control=list(fnscale=-1), method="Brent") $ par
Error in optim(0.01502, f, control = list(fnscale = -1), method = "Brent") : 
  'lower' and 'upper' must be finite values
> optim (0.01502, f, control=list(fnscale=-1), method="Brent", lower=0.01501, upper=0.01503) $ par
[1] 0.01502071
> ## esta sí vale
> abline (v=0.01502071, col=3)

Otra posibilidad es usar recocido desde un principio, por ejemplo mediante optim,

> optim (0.015, f, method="SANN", control=list(fnscale=-1, maxit=1e3)) [1:2]
$par
[1] 0.0154359

$value
[1] 3.231594

> optim (0.015, f, method="SANN", control=list(fnscale=-1, maxit=1e4)) [1:2]
$par
[1] 0.014878

$value
[1] 3.249706

> optim (0.015, f, method="SANN", control=list(fnscale=-1, maxit=1e5)) [1:2]
$par
[1] 0.01502001

$value
[1] 3.249512

> optim (0.015, f, method="SANN", control=list(fnscale=-1, maxit=1e6)) [1:2]
$par
[1] 0.01502044

$value
[1] 3.249922

> optim (0.015, f, method="SANN", control=list(fnscale=-1, maxit=1e6)) [1:2]
$par
[1] 0.015021

$value
[1] 3.249912

> optim (0, f, method="SANN", control=list(fnscale=-1, maxit=1e6)) [1:2]
$par
[1] 0.01502062

$value
[1] 3.249988

Lo mismo, pero aprovechando un ordenador con dos hilos de concurrencia, o sea, de procesamiento en paralelo.

Por ejemplo, usando la última versión:

unavez <- function (x) # "x" no se usa, pero hace falta para mclapply
  optim (0.015, f, method="SANN", control=list(fnscale=-1, maxit=1e6)) [1:2]
library(parallel)
dosveces <- mclapply (1:2, unavez)
## como SANN es aleatorio, hay que combinar soluciones,
## por ejemplo, promediando:
mean (c (dosveces[[1]]$value, dosveces[[2]]$value))
## o buscando directamente la mejor de las dos: 
i <- which.max(c(dosveces[[1]]$value, dosveces[[2]]$value))
dosveces[[i]]$par 

El siguiente programa resuelve el problema de las reinas para un tablero n×n. Adáptalo para resolver el problema de las n torres (como el de las reinas, pero con torres; las reinas comen en vertical, horizontal y diagonal; las torres, sólo en vertical y horizontal).

Simplemente, eliminamos las restricciones asociadas a la captura diagonal:

## variables: Xij i=1..n j=1..n
objetivo <- rep (1, n*n)
restric.fil <- kronecker (diag(n), t(rep(1,n)))
restric.col <- kronecker (t(rep(1,n)), diag(n))
## restric.di1 <- matrix (0, n+n-1, n*n)
## for (i in 2:(2*n))
## restric.di1 [i-1, apply(expand.grid(1:n, 1:n), 1, sum) == i] <- 1
## restric.di2 <- matrix (0, n+n-1, n*n)
## for (i in (-n+1):(n-1))
## restric.di2 [i+n, apply(expand.grid(1:n, 1:n), 1, diff) == i] <- 1
restricciones <- rbind (restric.fil, restric.col # , restric.di1, restric.di2
                        )
lado.derecho <- rep (1, nrow(restricciones))
sol <- lpSolve::lp ("max", objetivo, restricciones, "<", lado.derecho, all.int =TRUE)
print (sol)
print.table (matrix(sol$solution,n), zero.print="-")

La solución óptima es trivial, colocando las torres en diagonal.

Explica cómo hallar la ruta más corta que recorra varias ciudades. Por ejemplo, resuélvelo para:

> ciu = c("Madrid","Barcelona","Gibraltar","Lisbon","Lyons","Paris")
> as.matrix(eurodist)[ciu,ciu] # distancias kilométricas
>           Madrid Barcelona Gibraltar Lisbon Lyons Paris
Madrid         0       636       698    668  1281  1273
Barcelona    636         0      1172   1305   645  1033
Gibraltar    698      1172         0    676  1817  1971
Lisbon       668      1305       676      0  1178  1799
Lyons       1281       645      1817   1178     0   471
Paris       1273      1033      1971   1799   471     0

Los siguientes renglones aparecen en la ayuda de optim:

fw <- function (x) 10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x+80
res <- optim(50, fw, method = "SANN",
             control = list(maxit = 20000, temp = 20, parscale = 20))

Determina la relevancia de los parámetros maxit, temp y parscale para la obtención de ese res.

maxit

es el número de iteraciones del recocido; cuanto mayor sea, mayor será la probabilidad de que la solución se encuentre cerca del óptimo.

temp

indica el valor inicial del parámetro temperatura; cuanto más alto, son necesarias más iteraciones para llegar al estado de equilibrio, es decir, que converja la solución.

parscale

al igual que fnscale es un coeficiente asociado a la función objetivo (de forma que, si es negativo, permite maximizar), parscale es un coeficiente asociado al par (parámetro), es decir, a la solución (vector de valores sobre el que se aplica el vector objetivo).

fw <- function (x) 10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x+80
res <- optim(50, fw, method = "SANN",
             control = list(maxit = 20000, temp = 20, parscale = 20))
res[c("par","value")]
plot (fw, -16, -15, ylim=c(67.45,67.5), n=1001) # el mínimo cerca de -15'8...
## ...pero a veces (casi 20%) el recocido cae en el mínimo local de -15'66:
table (round (replicate (10000,
                         optim(50, fw, method = "SANN",
                               control = list(maxit = 20000,
                                              temp = 20,
                                              parscale = 20))$par), 1))
##  -16 -15.8 -15.7 
##   11  8058  1931 

resf <- function (maxit=20000, temp=20, parscale=20)
  optim(50, fw, method = "SANN",
    control = list(maxit = maxit, temp = temp, parscale = parscale))$par
table (round (replicate (10000, resf(maxit=50000)), 1))
## -15.8 -15.7 
##  9826   174 => casi se asegura convergencia correcta con 50k iteraciones
table (round (replicate (1000, resf(maxit=5e4,temp=100)), 1)) # 80%
table (round (replicate (1000, resf(maxit=5e4,temp=50)), 1)) # 88%
table (round (replicate (1000, resf(maxit=5e4,temp=30)), 1)) # 94%
table (round (replicate (10000, resf(maxit=5e4,temp=10)), 1)) # 99% mejor que temp=20
table (round (replicate (10000, resf(maxit=5e4,temp=5)), 1)) # 98%
table (round (replicate (1000, resf(maxit=5e4,temp=1)), 1)) # 85%
## para 50k iteraciones, bien con 5<temp<20, mejor en torno a 10
table (round (replicate (1000, resf(maxit=5e4,temp=10,parscale=10)), 1)) # 77%
table (round (replicate (1000, resf(maxit=5e4,temp=10,parscale=30)), 1)) # 99%
table (round (replicate (1000, resf(maxit=5e4,temp=10,parscale=50)), 1)) # 95%
table (round (replicate (1000, resf(maxit=5e4,temp=10,parscale=40)), 1)) # 99%
## bien con 20<parscale<40 ; es un parámetro muy relevante en este problema
## sorprendentemente, porque la explicación dada en ?optim sugiere que sirve
## sólo para problemas multidimensionales ; debe de afectar a la distancia
## de la nueva solución a la anterior