Posted on by

Traveling Sales Man Problem

# route inversion
TSP_inv<-function(route,len) {
  len_route=length(route)
  invP=sample(1:len_route,1)
  route[invP:min(len_route,(invP+len-1))]
    =rev(route[invP:min(len_route,(invP+len-1))])
  return(route)
}
> TSP_inv(1:10,4)
 [1]  1  2  3  4  5  9  8  7  6 10

# route switching
TSP_switch<-function(route){
  len_route=length(route)
  switchP=sample(1:len_route,2)
  route[switchP]=route[rev(switchP)]
  return(route)
}
> TSP_switch(1:10)
 [1]  1  2  3  4  5  6  7 10  9  8

# route translation
TSP_trans<-function(route){
  len_route=length(route)
  transP=sample(1:len_route,2);
  tempR=route[-transP[1]]
  if (transP[2]==1){
    tempR=c(route[transP[1]],tempR)
  } else if (transP[2]==len_route) {
    tempR=c(tempR,route[transP[1]])
  } else {
    tempR=c(tempR[1:(transP[2]-1)],route[transP[1]],tempR[(transP[2]):  
     (len_route-1)])
  }
  return(tempR)
}
> TSP_trans(1:10)
 [1]  1  3  4  5  6  7  8  9  2 10

# initialize cities' locations
FUN_initTSP<-function(n_city=10) {
  return(matrix(runif(n_city*2,1,100),nrow=n_city,ncol=2))
}
> print(loc<-FUN_initTSP(10))
          [,1]     [,2]
 [1,] 36.78996 40.33464
 [2,] 91.67296 97.33156
 [3,] 87.15730 58.82665
 [4,] 56.19289 44.84425
 [5,] 46.06971 43.22932
 [6,] 50.91508 30.14634
 [7,] 84.18521 56.17124
 [8,] 77.68784 43.77393
 [9,] 42.48220 11.74252
[10,] 12.99097 21.11037

# calc. total distance
FUN_costTSP<-function(route,cities) {
  route=c(route,route[1]);n_cities=nrow(cities);totDist=0;
  for (i_cities in 1:n_cities) {
    totDist=totDist+dist(cities[route[i_cities:(i_cities+1)],])
  }
  return(totDist)
}
> FUN_costTSP(1:10,loc)
         1
2 341.4417

# TSP demo
TSP_demo<-function(n_city=20, maxItr=1000) {
  loc=FUN_initTSP(n_city);route=1:n_city
  ## param. for simulated annealing - change values if necessary 
  C=1;eta=0.99;TEMP=50; 
  ##
  optDist=FUN_costTSP(route,loc)
  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)}
    new.dist=FUN_costTSP(new.route,loc)
    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(loc,loc[1,]),type='o',pch=20,cex=2.5, col='red',
    xlab='location @X',ylab='location @Y',main='Initial route')
  plot(loc[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)
}

TSP_plot

# 古いのバージョン
# initializing cities
FUN_initTSP<-function(n_city=10) {
  return(matrix(runif(n_city*2,1,100),nrow=n_city,ncol=2))
}
# calculating total (Euclidean) distance
FUN_costTSP<-function(route,cities) {
  route=c(route,route[1]);n_cities=nrow(cities);totDist=0;
  for (i_cities in 1:n_cities) {
    totDist=totDist+dist(cities[route[i_cities:(i_cities+1)],])
  }
  return(totDist)
}

# route inversion
TSP_inv<-function(route,len) {
  len_route=length(route)
  invP=sample(1:len_route,1)
  route[invP:min(len_route,(invP+len-1))]=rev(route[invP:min(len_route,(invP+len-1))])
  return(route)
}

# route switching
TSP_switch<-function(route){
  len_route=length(route)
  switchP=sample(1:len_route,2)
  route[switchP]=route[rev(switchP)]
  return(route)
}

# route translation / insertion
TSP_trans<-function(route){
  len_route=length(route)
  transP=sample(1:len_route,2);
  tempR=route[-transP[1]]
  if (transP[2]==1){
    tempR=c(route[transP[1]],tempR)
  } else if (transP[2]==len_route) {
    tempR=c(tempR,route[transP[1]])
  } else {
    tempR=c(tempR[1:(transP[2]-1)],route[transP[1]],tempR[(transP[2]):(len_route-1)])
  }
  return(tempR)
}

## main TSP script ###
n_city=20;loc=FUN_initTSP(n_city);route=1:n_city
maxItr=5000;C=1;eta=0.99;TEMP=20;
optDist=FUN_costTSP(route,loc)
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)}
  new.dist=FUN_costTSP(new.route,loc)
  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(loc,loc[1,]),type='b',pch=20,cex=2,
  xlab='location @X',ylab='location @Y',main='Initial route')
plot(loc[optHist[1000,c(1:n_city,1)],],type='b',pch=20,
  cex=2,xlab='location @X',ylab='location @Y',main='Optimized route')

## DEMO function ###
TSP_demo<-function(n_city=20, maxItr=1000) {
  loc=FUN_initTSP(n_city);route=1:n_city
  ## param. for simulated annealing - change values if necessary 
  C=1;eta=0.99;TEMP=20; 
  ##
  optDist=FUN_costTSP(route,loc)
  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)}
    new.dist=FUN_costTSP(new.route,loc)
    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(loc,loc[1,]),type='o',pch=20,cex=2.5, col='red',
    xlab='location @X',ylab='location @Y',main='Initial route')
  plot(loc[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)
}

Leave a Reply