## On Mac, use # quartz(antialias=FALSE) ## On Linux, use # X11(type='Xlib') ## On Windows, may work better with # windows(buffered = FALSE) ## "wild" function , global minimum at about -15.81515 fw <- function (x) 10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x+80 plot(fw, -50, 50, n=1000, main = "optim() minimising 'wild function'") optim(50, fw, method="SANN", control=list(maxit=2e5, temp=20, parscale=20)) fw2 <- function (x) { ## show proposal every 500 iterations niter <<- niter+1 if(niter %% 500 == 0) {abline(v=xold, col='white'); xold <<- x; abline(v=x);Sys.sleep(0.1)} 10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x+80 } xold <- NA; niter <- 1 plot(fw, -50, 50, n=1000, main = "optim() minimising 'wild function'") optim(50, fw2, method="SANN", control=list(maxit=2e5, temp=20, parscale=20)) ## Combinatorial optimization: Traveling salesman problem # distanes in km between 21 European cities. eurodistmat <- as.matrix(eurodist) distance <- function(sq) { # Target function sq2 <- embed(sq, 2) sum(eurodistmat[cbind(sq2[,2],sq2[,1])]) } step <- 1 genseq <- function(sq) { # Generate new candidate sequence idx <- seq(2, NROW(eurodistmat)-1) ## pick two cities at random and swap them in the sequence changepoints <- sample(idx, size=2, replace=FALSE) tmp <- sq[changepoints[1]] sq[changepoints[1]] <- sq[changepoints[2]] sq[changepoints[2]] <- tmp if(step %% 100 == 1) { ## show proposal every 100 iterations arrows(tspres[s,1], tspres[s,2], tspres[s+1,1], tspres[s+1,2], angle=10, col="white") tspres <<- loc[sq,] arrows(tspres[s,1], tspres[s,2], tspres[s+1,1], tspres[s+1,2], angle=10, col="red") Sys.sleep(0.1) } step <<- step + 1 sq } sq <- c(1:nrow(eurodistmat), 1) # Initial sequence: alphabetic # rotate for conventional orientation loc <- -cmdscale(eurodist, add=TRUE)$points x <- loc[,1]; y <- loc[,2] tspinit <- loc[sq,] s <- seq_len(nrow(eurodistmat)) dev.control('inhibit') plot(x, y, type="n", asp=1, xlab="", ylab="", main="initial solution of traveling salesman problem", axes = FALSE) text(x, y, labels(eurodist), cex=0.8) tspres <- tspinit set.seed(123) # chosen to get a good soln relatively quickly res <- optim(sq, distance, genseq, method="SANN", control = list(maxit=4e4, temp=2000, trace=TRUE, REPORT=500)) tspres <- loc[res$par,] plot(x, y, type="n", asp=1, xlab="", ylab="", main="optim() 'solving' traveling salesman problem", axes = FALSE) arrows(tspres[s,1], tspres[s,2], tspres[s+1,1], tspres[s+1,2], angle=10, col="red") text(x, y, labels(eurodist), cex=0.8)