In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo para encontrar o maximo de uma funcao bivariada pelo metodo de
# arrefecimento simulado (simulated anneling)
#==============================================================================
# Funcao do metodo de arrefecimento simulado (x,y)
#==============================================================================
maxf <- function(hx,x0,y0,temperatura,Delta1,Delta2,maxit){
xymax <- matrix(NA,maxit,2)
x <- x0
y <- y0
for(k in 1:maxit){
xprop <- runif(1,x-Delta1,x+Delta1)
yprop <- runif(1,y-Delta2,y+Delta2)
aux <- (log(hx(xprop,yprop))-log(hx(x,y)))/temperatura
prob <- min(1,exp(aux))
if (runif(1) < prob){
xymax[k,1] <- xprop
xymax[k,2] <- yprop
x <- xprop
y <- yprop
}else{
xymax[k,1] <- x
xymax[k,2] <- y
}
if((floor(k/100)-k/100)==0){
temperatura = temperatura*0.9
}
}
return(xymax)
}
In [2]:
#==============================================================================
# Funcao bivariada
#==============================================================================
hx <- function(x,y){
saida <- ( 0.60*exp(-0.5*((x-2)^2+(y-2)^2))
+ 0.35*exp(-0.5*((x-5)^2+(y-5)^2))
+ 0.05*exp(-0.5*((x-3)^2+(y-6)^2)) )
return(saida)
}
In [3]:
#==============================================================================
# Avaliando a funcao em uma grade de valores
#==============================================================================
xseq <- seq(-2,10,length=100)
yseq <- seq(-2,10,length=100)
zseq <- matrix(NA,100,100)
for(i in 1:100)
for(j in 1:100)
zseq[i,j] <- hx(xseq[i],yseq[j])
In [4]:
#==============================================================================
# Graficos da funcao
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
mar=c(4.5,5,1,1),cex.main=2.0)
persp(xseq,yseq,zseq,col="green",xlab=expression(x),theta=30,
ylab=expression(y),zlab=expression(f(x,y)),box=F)
In [5]:
#==============================================================================
# Curvas de nivel
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
contour(xseq,yseq,zseq,col="blue",lwd=4,xlab=expression(x),
ylab=expression(y))
In [6]:
#==============================================================================
# Encontrando o maximo pelo arrefecimento simulado
#==============================================================================
set.seed(12345)
xymax.1 <- maxf(hx,1,1,100,1,1,10000)
print(head(xymax.1))
print(tail(xymax.1))
[,1] [,2]
[1,] 1,441808 1,7515464
[2,] 2,214057 1,6645083
[3,] 1,864248 1,6829570
[4,] 2,843722 0,7520278
[5,] 3,315091 -0,2456990
[6,] 3,240081 -0,4694110
[,1] [,2]
[9995,] 1,965371 1,967533
[9996,] 1,965371 1,967533
[9997,] 1,965371 1,967533
[9998,] 1,965371 1,967533
[9999,] 1,965371 1,967533
[10000,] 1,965371 1,967533
In [7]:
#==============================================================================
# Encontrando o maximo pelo arrefecimento simulado
#==============================================================================
xymax.2 <- maxf(hx,7,7,100,1,1,10000)
print(head(xymax.2))
print(tail(xymax.2))
[,1] [,2]
[1,] 6,785296 7,669410
[2,] 6,246879 6,949912
[3,] 6,655508 7,343412
[4,] 6,912796 6,832714
[5,] 7,386606 7,244876
[6,] 7,952625 6,250445
[,1] [,2]
[9995,] 1,967688 2,033776
[9996,] 1,967688 2,033776
[9997,] 1,967688 2,033776
[9998,] 1,967688 2,033776
[9999,] 1,967688 2,033776
[10000,] 1,967688 2,033776
In [8]:
#==============================================================================
# Grafico: iteracoes
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(6,8,5),
mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
plot(xymax.1[,1],type="l",col="darkgreen",xlab="Iteração",
ylab=expression(x),ylim=c(-10,30),lwd=4)
lines(c(0,10000),c(2,2),lwd=4,lty=2)
plot(xymax.1[,2],type="l",col="darkgreen",xlab="Iteração",
ylab=expression(y),ylim=c(-10,15),lwd=4)
lines(c(0,10000),c(2,2),lwd=4,lty=2)
In [9]:
#==============================================================================
# Grafico: iteracoes
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(6,10,5),
mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
plot(xymax.2[,1],type="l",col="blue",xlab="Iteração",
ylab=expression(x),ylim=c(-15,15),lwd=4)
lines(c(0,10000),c(2,2),lwd=4,lty=2)
plot(xymax.2[,2],type="l",col="blue",xlab="Iteração",
ylab=expression(y),ylim=c(-15,15),lwd=4)
lines(c(0,10000),c(2,2),lwd=4,lty=2)
In [10]:
#==============================================================================
# Grafico: iteracoes
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,5,1,1),cex.main=2.0,bty="n")
plot(xymax.1,type="l",xlim=c(-10,30),ylim=c(-10,15),lwd=2,
xlab=expression(x),ylab=expression(y),col="darkgreen")
points(2,2,pch=16,cex=2,col="red") # Alvo
points(1,1,pch=15,cex=2,col="black") # Inicio
plot(xymax.2,type="l",xlim=c(-15,15),ylim=c(-15,15),lwd=2,
xlab=expression(x),ylab=expression(y),col="blue")
points(2,2,pch=16,cex=2,col="red") # Alvo
points(7,7,pch=15,cex=2,col="black") # Inicio
In [11]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================