traveling salesman problem
来週までに、inversionとtranslationを行う関数を作っておいてください。

# initialisation, etc.
n.city=20
location=matrix(runif(n.city*2,min=0,max=100),nrow=n.city)
route=sample(n.city) 
calc.Dist<-function(location,route) {
  n.city=length(route)
  total.Dist=0
  route=c(route,route[1])
  for (i_city in 1:n.city){
    total.Dist=total.Dist+dist(location[c(route[i_city],route[i_city+1]),])
  }
  return(total.Dist)
}

# plotting results
plot(rbind(location[route,],location[route[1],]),type='o',pch=20,cex=2.5, col='red',
  xlab='location @X',ylab='location @Y',main='Initial route')
plot(rbind(location[best.route,],location[best.route[1],]),type='o',pch=20, 
  col='magenta',cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')

#switching 
TSP.switch<-function(route) {
  two.cities=sample(route,2,replace=F)
  route[two.cities]=route[rev(two.cities)]
  return(route)
}

# inversion
TSP.inv<-function(route,inv.length) { 
  inv.begin=sample(route,1)
  inv.end=min(inv.begin+inv.length-1,length(route))
  route[inv.begin:inv.end]=rev(route[inv.begin:inv.end])
  return(route)
}

# translation
TSP.trans<-function(route,tr.length) {
  trP1=sample(route,1)
  tr.vec=route[trP1:min(length(route),trP1+tr.length-1)]
  temp.vec=route[-(trP1:min(length(route),trP1+tr.length-1))]
  trP2=sample(1:length(temp.vec),1)
  if (trP2==length(temp.vec)) {
    new=c(temp.vec,tr.vec)
  } else {
    new=c(temp.vec[1:trP2],tr.vec,temp.vec[(trP2+1):length(temp.vec)])
  }
return(new)
}

# demo
TSP.demo<-function(n.city=20, maxItr=1000) {
  location=matrix(runif(n.city*2,min=0,max=100),nrow=n.city)
  route=sample(n.city) 
  ## param. for simulated annealing - change values if necessary 
  C=1;eta=0.99;TEMP=50; 
  ##
  optDist=calc.Dist(location,route)
  optHist=matrix(0,nrow=maxItr,ncol=(length(route)+1))
  optHist[1,]=c(route,optDist)
  for (i_loop in 2:maxItr) {
    rand.op=sample(c('inv','sw','trans'),1,prob=c(0.75,0.125,0.125))
    if (rand.op=='inv') {
      new.route=TSP.inv(route,sample(2:n.city,1))
    } else if (rand.op=='sw') {
      new.route=TSP.switch(route)
    } else { new.route=TSP.trans(route,sample(2:(round(n.city/2)),1))}
    new.dist=calc.Dist(location,new.route)
    delta=new.dist-optDist
    prob=1/(1+exp(delta/(C*TEMP)))
    if (runif(1) < prob) {
      route=new.route;optDist=new.dist;
    } 
    optHist[i_loop,]=c(route,optDist);
    TEMP=TEMP*eta
  }
  par(mfrow=c(1,2)) 
  plot(rbind(location,location[1,]),type='o',pch=20,cex=2.5, col='red',
    xlab='location @X',ylab='location @Y',main='Initial route')
  plot(location[optHist[1000,c(1:n.city,1)],],type='o',pch=20, col='magenta',
    cex=2.5,xlab='location @X',ylab='location @Y',main='Optimized route')
  return(optHist)
}

最適化問題

 
#Objective Function 
funZ<-function(x,y) {x^4-16*x^2-5*x+y^4-16*y^2-5*y}
xmin=-5;xmax=5;n_len=100;
xs=ys<-seq(xmin, xmax,length.out=n_len)
z<-outer(xs,ys,funZ)
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))

勾配法の実装例

 
# initialization
x=x.hist=runif(1,min=-5,max=5) # initial value of x
y=y.hist=runif(1,min=-5,max=5) # initial value of y
lambda=0.01                    # step size
zeta=1e-8                      # tol (stopping criterion)
t=0                            # time
g=100                          # maximum gradient
# end initialization

### note ###
# obj fun: z=x^4-16*x^2-5*x+y^4-16*y^2-5*y
# part. deriv. x: 4*x^3-32*x-5
# part. deriv. y: 4*y^3-32*y-5
#############

# main 
while (g>zeta) {
  t=t+1
  g.x=4*x^3-32*x-5
  g.y=4*y^3-32*y-5
  x=x-lambda*g.x
  y=y-lambda*g.y
  g=max(abs(g.x),abs(g.y))
  x.hist=c(x.hist,x)
  y.hist=c(y.hist,y)
}

# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

単純確率的最適化法の実装例

 
# initialization
x=runif(1,min=-5,max=5)             # initial value of x
y=runif(1,min=-5,max=5)             # initial value of y
f.xy=x^4-16*x^2-5*x+y^4-16*y^2-5*y  # initial value of obj. fun.
w=0.1                               # search width
t=0                                 # time
max.t=1000                          # max iter (stopping criterion)
stop.crit=F                         # stop criterion
x.hist=rep(0,max.t)                 # history of x
y.hist=rep(0,max.t)                 # history of x
x.hist[1]=x
y.hist[1]=y
# end initialization

# main
while (stop.crit==F) {
  t=t+1
  x.temp=x+rnorm(1,mean=0,sd=w)
  y.temp=y+rnorm(1,mean=0,sd=w)
  f.temp=funZ(x.temp,y.temp)
  if (f.temp < f.xy) {
    x=x.temp
    y=y.temp
    f.xy=f.temp
  }
  x.hist[t]=x
  y.hist[t]=y
  if (t==max.t){stop.crit=T}
}

# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

焼きなまし法の実装例

# initialization
x=runif(1,min=-5,max=5)             # initial value of x
y=runif(1,min=-5,max=5)             # initial value of y
f.xy=x^4-16*x^2-5*x+y^4-16*y^2-5*y  # initial value of obj. fun.
t=0                                 # time
max.t=1000                          # max iter (stopping criterion)
stop.crit=F                         # stop criterion
x.hist=rep(0,max.t)                 # history of x
y.hist=rep(0,max.t)                 # history of x
x.hist[1]=x
y.hist[1]=y
tau=1
c=1
g=0.999
# end initialization
 
# main
while (stop.crit==F) {
  t=t+1
  x.temp=x+rnorm(1,mean=0,sd=tau)
  y.temp=y+rnorm(1,mean=0,sd=tau)
  f.temp=funZ(x.temp,y.temp)
  deltaE=f.temp-f.xy
  p=1/(1+exp(deltaE/(c*tau)))
  if (runif(1) < p) {
    x=x.temp
    y=y.temp
    f.xy=f.temp
  }
  x.hist[t]=x
  y.hist[t]=y
  tau=tau*g
  if (t==max.t){stop.crit=T}
}
 
# plotting results
par(mfrow=c(1,3))
par(mar=c(4.0,4.1,4.1,1.1))
contour(x=xs,y=ys,z=z,nlevels=30,drawlabels=F,col="darkgreen",
 main=expression(z=x^4-16*x^2-5*x+y^4-16*y^2-5*y))
lines(x.hist,y.hist,type='o',pch=20,col='red')
plot(x.hist,type='l',xlab="time",ylab="value of x",lwd=2,
 main="how value of y changes over time")
plot(y.hist,type='l',xlab="time",ylab="value of y",lwd=2,
 main="how value of y changes over time")

３つの手法の比較

 

GD<-function(x.init,y.init,n.iter,lambda){
x=x.init;y=y.init;z.hist=rep(0,n.iter)
# main 
for (i_iter in 2:n.iter) {
  g.x=4*x^3-32*x-5
  g.y=4*y^3-32*y-5
  x=x-lambda*g.x
  y=y-lambda*g.y
  g=max(abs(g.x),abs(g.y))
  x.hist=c(x.hist,x)
  y.hist=c(y.hist,y)
  z.hist[i_iter]=funZ(x,y)
}
return(z.hist)
}

stochOpt<-function(x.init,y.init,n.iter,w) {
x=x.init;y=y.init;f.xy=funZ(x,y)
z.hist=rep(0,n.iter);z.hist[1]=f.xy
# main
for (i_iter in 2:n.iter) {
  x.temp=x+rnorm(1,mean=0,sd=w)
  y.temp=y+rnorm(1,mean=0,sd=w)
  f.temp=funZ(x.temp,y.temp)
  if (f.temp < f.xy) {
    x=x.temp;y=y.temp;f.xy=f.temp
  }
  z.hist[i_iter]=f.xy
}
return(z.hist)
}
 
simAnn<-function(x.init,y.init,n.iter,tau,c,g) {
x=x.init;y=y.init;f.xy=funZ(x,y)
z.hist=rep(0,n.iter);z.hist[1]=f.xy
# main
for (i_iter in 2:n.iter) {
  x.temp=x+rnorm(1,mean=0,sd=tau)
  y.temp=y+rnorm(1,mean=0,sd=tau)
  f.temp=funZ(x.temp,y.temp)
  deltaE=f.temp-f.xy
  p=1/(1+exp(deltaE/(c*tau)))
  if (runif(1) < p) {
    x=x.temp;y=y.temp;f.xy=f.temp
  }
  z.hist[i_iter]=f.xy
  tau=tau*g
}
return(z.hist)
}

n.rep=500
n.iter=500
xs=runif(n.rep,-3,3);ys=runif(n.rep,-3,3)
GD.hist=SA.hist=SO.hist=matrix(0,n.iter,n.rep)
for (i.rep in 1:n.rep) {
  GD.hist[,i.rep]=GD(xs[i.rep],ys[i.rep],n.iter,0.025)
  SO.hist[,i.rep]=stochOpt(xs[i.rep],ys[i.rep],n.iter,2)
  SA.hist[,i.rep]=simAnn(xs[i.rep],ys[i.rep],n.iter,3,1,0.999)
}

plot(rowMeans(GD.hist),type='l',ylim=c(-160,0),lwd=3)
lines(rowMeans(SO.hist),col='red',lwd=3)
lines(rowMeans(SA.hist),col='blue',lwd=3)