Solving the Vehicle Routing Problem with Simulated Annealing

1 Introduction and Motivation

The Vehicle Routing Problem (VRP) is a classic operations research problem concerned with finding the optimal routing of a fleet of k vehicles seeking to reach n nodes. Each vehicle must begin its route at a shared source node and return to said source node at the end of its path. Optimality can be defined in a few ways but generally includes minimizing some quantity such as distance or time. For example, we might want to minimize the distance travelled across all vehicles in an effort to minimize fuel usage. Or, we may want to minimize the time until the last vehicle returns to the source thereby setting a cap on the amount of time to serve all n nodes.

Another way of thinking about this problem is as a generalization of the classic Travelling Salesman Problem (TSP). The TSP involves a wandering salesman seeking to sell his goods to customers in n far off cities. The salesman wants to minimize the amount of time to reach each of these cities and return to his starting place. Generalizing to the VRP from the TSP, instead of one salesperson, we now have k salespeople who agree that each destination will be serviced by one member of the sales team before the salespeople return home to the source. Notice that the VRP must have a fixed source unlike the TSP since the TSP’s solution generates a single cycle in the network and thus the starting point does not affect optimality. In contrast, the starting point of the VRP defines the distinct cycles each vehicle or salesperson must complete. Secondly, the VRP contains an assignment problem buried into it. To find the optimal solution, we have to find the optimal assignment of nodes to vehicles along with the optimal path for each vehicle to its assigned node.

An important shared trait of the TSP and the VRP is that neither of them can be solved simple, efficient algorithms. Just like many combinatorial optimization problems, they fall into the class of NP-Hard problems. The feasible set of solutions (or in this case routing choices) is enormous and grows exponentially with the problem size. For a given VRP with n destinations and k vehicles, there are \(k^{n-1}\) valid assignments of destinations to vehicles plus the sum of within common assignment permutations of destinations. Clearly, enumerating over the entire space is not remotely feasible for larger problems. Finally, one of the other challenges of combinatorial optimization problems is that the cost function typically contains a large number of local optima that an approach such as hill-climbing will be unable to overcome. These challenges motivate the use of simulated annealing.

Simulated annealing presents a very plausible approach to the VRP given its ability to maneuver between peaks and valleys in the cost function to avoid getting caught in a local optima. The random proposals aspect is absolutely necessary for this movement to happen. Then, as the peaks become sharper, we should expect to find the optimal routing.

2 Problem Formulation

In more formal terms, the Vehicle Routing Problem can be formulated as follows. If our objective is to minimize the total distance travelled by the vehicles, the problem becomes:

\[ \begin{aligned} \min &\sum_{i=1}^k \sum_{j=2}^{n_k} d(x_{i(j-1)},x_{ij}) \\ \text{subject to} & \sum_{i=1}^k n_k = n + 2k - 1 \\ & x_{i1} = x_{in_k} = x_1 \qquad \forall i \in 1,2,\dots,k \end{aligned} \]

where the variables are defined as follows:

• \(k\) is the number of vehicles in the fleet
• \(n\) is the number of destinations in the network
• \(n_k\) is the number of nodes vehicle \(k\) will visit including starting and ending at the source node
• \(x_{ij}\) is the \(j\)th destination in vehicle \(i\)’s route
• \(d(x_{i(j-1)},x_{ij})\) is the distance measure computing the distance between destination \(j-1\) and \(j\) on vehicle \(i\)’s route. In this report we will be using Euclidean distance.
• \(x_1\) is the source node

The functions computing the distance and cost of a given routing are available in Appendix A1.

As specified above, the problem operates over a graph with \(n\) nodes and a designated source node. For the purposes of examining the performance of simulated annealing on the VRP, we can generate sample networks by plotting \(n\) uniformly generated nodes in the unit square. We will leave the option to set a specified starting node to illustrate the impact of different source nodes on the solutions. Below are a set of four sample networks. One of them is completely random, the other three have their source nodes set at \((0,0)\), \((0.5,0)\), and \((0.5,0.5)\). The code to generate these plots can be found in Appendix A3.

The above plot highlights the earlier discussion about sensitivity to source location. For example, the (0.5,0.5) source could send a vehicle into each quadrant of the plot, however the (0,0) source plot does not have that luxury. Some vehicles sent from this source must cross the entire map to reach the far off destinations. Having formulated the problem, we can move forward with searching for the optimal routings.

3 Classical Monte Carlo

The most basic simulation based approach to solving this problem is to just uniformly generate a sample of routings and hope that the optimal solution is found within the sample. However, we know that this method relies on sheer luck of the draw as the combinatorial explosion in this problem is immense. However, just to confirm, we can try to find optimality for the random source example plotted above.

With a sample size of 10000, the best randomly generated pathing is visibly incorrect. The code to execute this method, the plot of the costs, and the plot of the routes can be found in Appendix A2.

4 Basic Metropolis Algorithm

Before investigating the simulated annealing approach, we can first attempt to use a random-walk Metropolis algorithm to find the optimal solution. The process to find the optimal solution is similar to the uniform samples in that we want to draw a sample and pick the best one. However, we will be proposing a new routing solution that neighbours the current solution and accepting based on a function of the change in costs. We will define a neighbouring solution as the current solution with two randomly chosen nodes swapping position in the ordering and one node being randomly re-assigned to a vehicle. Note that the chosen nodes cannot include the source as it is fixed for all vehicles, and that the nodes being swapped can be chosen from any vehicle. As the changes being made to the current solution are uniform, the proposals are reversible.

To model the density of a given state, we can treat the cost of the current and proposed states as the energy of that state and use the Boltzmann distribution excluding the Boltzmann constant to model the probability of being in a given state. This is the approach typically used in simulated annealing. The basic Metropolis Algorithm then reduces to a case of simulated annealing where we hold the temperature fixed at 1 (though there wouldn’t be much annealing happening)! Mathematically, a proposal is then accepted with probability \(p=\min\left(1,\frac{e^{-C(x_{n+1})}}{e^{-C(x_n)}}\right)=\min\left(1,e^{-(C(x_{n+1})-C(x_n))}\right)\) where \(C(\cdot)\) is the cost of the state. Intuitively, if we are able to find a new state with lower cost, the proposal will always be accepted. If not, there is some probability of rejecting.

Notice that this structure has produced an algorithm with the following properties. The constructed chain will be reversible and thus have a stationary distribution, aperiodic since there is a positive probability of rejecting a proposal, and irreducible since it takes at most n steps to get from one state to any other state. Thus, our Markov chain will converge to the density modelled by the energy of the states. Equipped with this, if we run the chain we should hope that for some large enough sample we will choose the optimal solution at some point. Thus to find optimality, we can try to generate a large sample and pull the minimum cost sample. The code for this algorithm can be found in Appendix A4.

Once again, we attempt to find the optimal routing for the sample network with a random source node. Just like with the simple uniform sample, we were unable to find the optimal routing despite a large sample of 100 000. Just looking at the selected routes plot below, we see that vehicle 2 is all over the map and intersecting with the routes of other vehicles. It is quite obvious that this is not the optimal routing. When we look at the trace plot and the acf plot, it doesn’t appear that the chain is performing poorly, it seems that this is a question of the algorithm being poorly suited to the optimization task. It is not finding the specific route that we are searching for.

5 Simulated Annealing Approach

Having had the first two methods fail, we turn to simulated annealing, the focus of our solution. The underlying idea behind simulated annealing is that by varying the temperature, meaning the powers of the density, and slowly transitioning from a flattened distribution to a sharp distribution allows the algorithm to explore the space easily in the flattened distribution and then be forced into the sharpest peak by the time the temperature cools at the end of the algorithm. This sharpest peak corresponds to the optimal value of our cost function. Unlike the previous two approaches, as the density sharpens, we won’t accept poor proposals and are forced into the optimal solution.

Recall previously we had that the energy of the system was represented by \(\pi=e^{-C(x_n)}\). Taking this to the power of one over the temperature \(T\), the quantity \(\pi^{\frac{1}{T}}=e^{-\frac{C(x_n)}{T}}\) is pulled towards T for large values of T and is sharpened for values of T less than 1. Introducing this new result to our accept-reject criteria produces, \(p=\min\left(1,\frac{e^{-\frac{C(x_{n+1})}{T}}}{e^{-\frac{C(x_n)}{T}}}\right)=\min\left(1,e^{-\frac{C(x_{n+1})-C(x_n)}{T}}\right)\). Once again, a proposal with a lower cost value will always be accepted. The difference here is that the likelihood of accepting a worse proposal drastically decreases as the temperature decreases. To decrease the temperature, we will test the following cooling schedules:

• Geometric: \(T_n = T_0r^n\)
• Linear: \(T_n = T_0-dn\)
• Logarithmic: \(T_n = \frac{cT_0}{\log(1+n)}\)
• Rational: \(T_n = \frac{T_0}{cn^2}\)
• Root: \(T_n = T_0-d\sqrt{n}\)

As a first test of the performance of simulated annealing relative to the first two algorithms, we can attempt to run the algorithm on the four samples of 15 destinations and 4 vehicles above. We will begin with the geometric cooling scheudle with a sample size of 20 000.

## Cost 1: 4.789146 Cost 2: 4.761482

The first run of the simulated annealing algorithm produces drastically better results than the first two solutions. In each of the four networks, each vehicle makes a distinct path towards its own group of destinations. This is in contrast to the earlier approaches where we found overlapping routes and obvious departures from optimality. However, notice in the (0,0) source network that there appears to be a slight deviation from the expected optimal routing in that the green path travels to a destination that would appear to be more aligned with the blue path. If we manually reallocate this node to the blue path, we find that the cost drops from 4.789146 to 4.761482. The details of this computation along with the code for the simulated annealing algorithm can be found in Appendix A5. While this is just a slight improvement over our current solution, it does indicate that the simulated annealing solution is not perfect.

Another interesting result is that the routings do not necessarily make use of all available vehicles. This is in an effort to minimize the impact of expensive moves from the source to the extreme regions of the network. This is particularly apparent in the (0,0) source case as we see long stretches from the source into the network. Such a result would be highly unlikely to identify without the cooling process involved in simulated annealing. As an example, if we just consider the assignment component of the problem, there are \(\binom{4}{2} 2^{n-1}=3\cdot2^n\) assignment options with only two vehicles as compared to the total of \(4^{n-1}\) total valid assignments. In the current case with \(n=15\), this represents a proportion of less than 0.0002. This is not a situation we would encounter under the previous algorithms, expecially once we grow the problem size.

While the simulated annealing approach is significantly better than the sampling approaches, it is not perfect. We’ve seen that in the difficult case of the source being placed in the extremes of the network, we did not find the optimal solution in the first run. There are a few remedies we can employ. First, increasing the sample size will offer two possible solutions. A larger sample size offers more proposal opportunities and also slows down the cooling schedule to allow more time for mixing. In Appendix A6, plots of a higher sample size are provided for the (0,0) case. It turns out that we can much further reduce the cost of this case by only using 1 vehicle - a result that required a very large sample size of 100 000 and would be extremely unlikely to find with simple random generation. Looking at the trace plot of this large sample run, the algorithm spend the first 50 000 mixing before begining to descend towards the optimal value. The extra time mixing was necessary to find the obscure scenario where three vehicles remain at the source.

6 Comparing Cooling Schedules on a Larger Test Case

Having seen the generally postive results of the simulated annealing algorithm on the vehicle routing problem, we can try to extend to a larger network size and compare the different cooling schedules. Of the five schedules described above, we will look at the geometric, linear, rational, and root methods. The convergence rate of the logarithmic cooling schedule is far too slow to be a viable method. On a problem size of \(n=30\), \(k=4\), and a random source node, we try each method with a sample size of 100 000 and a cooling schedule starting at 500 and decreasing to 0.01.

## Geometric Cost: 5.065856 Linear Cost: 11.82136

## Rational Cost: 5.404348 Root Cost: 10.36937

Looking at the proposed routes, it is clear that the linear and root cooling schedules are not able to solve this problem. The route plots exhibit obvious departures from optimality in that the routes are sending vehicles across the map to locations that could be serviced by another vehicle. This point is echoed by the optimal cost values. The optimal costs of the linear and root cooling schedules are much higher than those of the geometric and rational cases. When looking at the respective trace plots and cooling schedules in Appendix A7, we can see that the linear and root cooling schedules are still in the mixing phase and have not begun descending towards an optimal solution by the end of the sample.

In comparing the geometric and rational cooling schedules, the geometric schedule is our immediate preference given that it has a lower cost. However, the geometric solution still does not look fully optimal. Vehicle 3 forms a loop in its path which indicates that there exists an improved solution. That being said, the solution is still well chosen and certainly better than the other algorithms would have produced. Looking at the trace plots of the geometric and rational cooling runs provides insight into why the geometric cooling schedule performed better. The geometric schedule spends more time mixing before descending towards a minimum wheras the rational schedule begins its descent almost immediately. This lack of mixing pushed it towards an early convergence to a suboptimal solution. While this schedule does not find as optimal a solution, notice that it did converge using roughly half of the sample size and thus could be looked to as a more efficient proxy to the geometric schedule if it is too computationally expensive.

If we try to run the geometric schedule 10 times over and take the minimum cost, our new route avoids the loop observed earlier. The optimal cost value has also decreased. These results are in Appendix A7.

In summary, the simulated annealing algorithm constructed does scale to slightly larger problems, however we are not finding the global minimum despite the large 100 000 sample size. We are finding very close local optima especially when using the geometric cooling schedule. The geometric cooling schedule did outclass most of the other options though the rational cooling schedule did offer a valid solution as well.

7 Extending from Distance to Time

We’ve seen that simulated annealing offers promising results in finding the optimal routing to minimize the total distance travelled by a fleet of vehicles. However, often times we may be more interested in the time to complete a task as opposed to the distance. Introducing time as an objective substitutes the deterministic nature of the distance with a stochastic problem instead.

Suppose we are interested in minimizing the maximum amount of time a vehicle is out for deliveries. While our constraints on the pathing are identical to the distance problem we have been investigating so far, the objective in our problem formulation becomes, \[ \min \left( \max_{i\in(1,k)} \sum_{j=2}^{n_k} t(x_{i(j-1)},x_{ij})\right) \] where t is the time to move between destinations. Time, unlike its distance counterpart is not static. The time to travel between points could be hampered by the weather, traffic, or any number of other unexpected events. To model this uncertaintly, let’s assume that one unit of distance can on average be travelled in one unit of time and add a standard normal jitter to the average travel time with variance dependent on the distance (say standard deviation is one quarter the distance/average time). Finally, to make our model more realistic, let’s introduce dependency within a route. If a leg of a route differs from the average time to complete that segment, scale the average of the next leg by the same percentage deviation. A new cost function modelling this updated objective and cost can be found in Appendix A8.

Since the cost function is no longer constant for a given ordering, we are really searching to minimizing the expected maximum vehicle time. As we have specified the distributions of the time data, we can compute Monte Carlo estimates of the cost of a given solution at every iteration! We are restricting the evaluation of the time problem to the smaller 15 node network with a random source as the expected value computation each iteration is expensive. Looking at the results in A8, we see that the time based solution is much more complex. All 4 vehicles are allocated notes to travel to as expected but the randomness involved in the computation has slowed down convergence. Vehicle 2 seems to be intruding on the area dominated by Vehicle 1 which indicates that we are likely not at optimality.

8 Conclusions

In summary, we found that simulated annealing is able to produce high quality solutions to the vehicle routing problem and that simpler methods fail to come close to matching its performance. However, in some of the more difficult cases such as those with source nodes in the extremities of the network, the algorithm does have trouble finding the optimal routing. This was remedied by drastically increasing the sample size to slow down the cooling schedule.

The algorithm was able to extend to a larger problem but we were unable to avoid falling into local maximums near the global optimal solution. While we did consider large sample sizes to slow down the cooling, this was still not enough to guarantee convergence. We also found that in general, the geometric cooling schedule was the best of the tested options. The alternatives, such as linear, often yielded poor results that reflect their lack of convergence. The minor variances from the optimal solution found under the geometric cooling schedule should not take away from the successes of simulated annealing to the VRP. What is an explosive combinatorial optimization problem has been reduced to random proposals and temperature cooling.

Finally, we extended the VRP to a time based problem. In this situation, we are optimizing an expected value and thus introduce Classical Monte Carlo to our cost function. Given the change in objective, the results of the time and distance problems differ substantially. As a next step, we can consider updating the routing as information regarding the random variation in travel time becomes apparent.

9 Appendix

9.1 A1 Distance and Cost

#Define a function that calculates the Euclidean distance
dist <- function(a,b){
  sqrt(((a[1]-b[1])^2) + ((a[2]-b[2])^2))
}

#Define a function to compute the cost of a given allocation and ordering
#Note that x is a (n-1)x2 matrix where the first column contains a number from 
#2 to n with the name of the destination number and the second column is which 
#of the k delivery vehicles is assigned to the destination. The cost is computed 
#as the sum of the distances between destinations assigned to a specific vehicle
#in the order that they appear. 
#Note that loc is the location matrix of the destinations that characterizes our
#network. The first row corresponds to the source node which will all vehicles 
#must originate from and all vehicles must return to
cost <- function(x, loc){
  k <- max(x[,2])
  #print(x)
  d <- rep(0,k)
  cur <- matrix(data=rep(loc[1,],k), nrow=k, ncol=2, byrow=T)
  for (i in 1:nrow(x)){
    dest <- x[i,1]
    vehicle <- x[i,2]
    d[vehicle] <- d[vehicle] + dist(loc[dest,],cur[vehicle,])
    cur[vehicle,] <- loc[dest,]
  }
  #return each vehicle to the source
  for (i in 1:k){
    d[i] <- d[i] + dist(loc[1,], cur[i,])
  }
  #print(d)
  sum(d)
}

9.2 A2 IID Samples

set.seed(2)
#cmc is a function that uniformly generates possible routings and selects the sample 
#that produced the minimum cost
cmc <- function(M, dests, k, plt=FALSE){
  costlist <- rep(0,M)
  opt <- 100000
  ans <- NA
  for (i in 1:M){
    cur <- data.frame(order=sample(c(2:nrow(dests))), 
                      vehicle=sample(c(1:k), size=(nrow(dests)-1), replace=TRUE))
    costlist[i] <- cost(cur,dests)
    if (costlist[i] < opt){
      opt <- costlist[i]
      ans <- cur
    } 
  }
  if (plt){
    df <- data.frame(Index=c(1:M), Cost=costlist)
    print(ggplot(data=df, aes(Index, Cost))+geom_line()+ggtitle("IID Cost Estimates"))
  }
  ans
}


plotpath(cmc(10000, sample1, 4, plt=TRUE), sample1,"Uniform Sampling Route")

9.3 A3 Network Generation

set.seed(1)

#Define a function to generate a network with n nodes. If source is a 2 element 
#vector, it is set as the source node 
networkgen <- function(n, source=NA){
  if (any(is.na(source))){
    matrix(data=runif(2*n), nrow=n, ncol=2)
  } else {
    matrix(data=c(source,runif(2*(n-1))), nrow=n, ncol=2, byrow=TRUE)
  }
}

sample1 <- networkgen(15)
sample2 <- networkgen(15, source=c(0,0))
sample3 <- networkgen(15, source=c(0.5,0))
sample4 <- networkgen(15, source=c(0.5,0.5))

#networkplot takes a matrix of (x,y) coordinates by first converting it to a
#data frame and then producing a ggplot object that can be printed
networkplot <- function(m, title){
  df <- data.frame(xcoord=m[,1], ycoord=m[,2])
  ggplot(data=df, aes(x=xcoord,y=ycoord)) + 
    geom_point(size=3, alpha=1/3) + 
    geom_point(data=df[1,], color="steelblue", size=3) + 
    ggtitle(title)
}

#Define a function to overlay a routing onto the network
plotpath <- function(path, loc, title){
  k <- max(path[,2])
  source <- data.frame(order=rep(1,k),vehicle=seq(1,k))
  df <- source %>% 
    bind_rows(path) %>% 
    bind_rows(source) %>% 
    full_join(data.frame(order=c(1:nrow(loc)), xcoord=loc[,1], ycoord=loc[,2]), by="order")
  ggplot(data=df, aes(x=xcoord, y=ycoord, group=factor(vehicle))) + 
    geom_point(size=3, alpha=1/3) +
    geom_point(data=df[1,], color="steelblue", size=3) +
    geom_path(aes(color=factor(vehicle))) +
    ggtitle(title) +
    labs(color="Vehicle")
}

p1 <- networkplot(sample1, "Source: Random")
p2 <- networkplot(sample2, "Source: (0,0)")
p3 <- networkplot(sample3, "Source: (0.5,0)")
p4 <- networkplot(sample4, "Source: (0.5,0.5)")

grid.arrange(p1,p2,p3,p4, nrow=2, top="Sample Networks (n=15)")

9.4 A4 Random-Walk Metropolis

#metro is a function to generate a random-walk Metropolis MCMC sample of the possible routings.
metro <- function(M, B, dests, k, plt=FALSE){
  costlist <- rep(0,(M+B))
  opt <- 100000
  ans <- NA
  acceptnum <- 0 
  
  #Initialize the ordering
  cur <- data.frame(order=sample(c(2:nrow(dests))),
                    vehicle=sample(c(1:k), size=(nrow(dests)-1), replace=TRUE))
  
  #set the initial cost
  costlist[1] <- cost(cur, dests)
  
  for (i in 2:(M+B)){
    #Select a pair of nodes to swap in the ordering
    orderswap <- sample(c(1:nrow(cur)), size=2)
    proposal <- cur
    proposal[orderswap[1],] <- cur[orderswap[2],]
    proposal[orderswap[2],] <- cur[orderswap[1],]
    
    #Select a node to swap vehicle assignment
    vswap <- sample(c(1:nrow(proposal)), size=1)
    proposal[vswap,2] <- sample(c(1:k), size=1)
    
    #Compute the cost of the proposal and the uniform sample for accept reject
    propcost <- cost(proposal, dests)
    U <- runif(1)
    alpha <- exp(-(propcost - costlist[i-1]))
    if (U < alpha){
      costlist[i] <- propcost
      cur <- proposal
      acceptnum <- acceptnum + 1
      if (propcost < opt){
        ans <- proposal
      }
    } else {
      costlist[i] <- costlist[i-1]
    }
  }
  if (plt){
    df <- data.frame(Index=c(1:M), Cost=costlist[(B+1):(M+B)])
    trace <- ggplot(df, aes(x=Index, y=Cost)) + geom_line()
    acftemp <- acf(costlist[(B+1):(M+B)], plot=FALSE)
    acfdata <- with(acftemp, data.frame(lag,acf))
    acfplot <- ggplot(acfdata, aes(x=lag,y=acf)) + 
      geom_hline(aes(yintercept=0)) + 
      geom_hline(aes(yintercept=0.05), color="steelblue", linetype="dashed") + 
      geom_hline(aes(yintercept=-0.05), color="steelblue", linetype="dashed") +
      geom_segment(aes(xend=lag, yend=0))
    empty <- networkplot(dests, "Source: Random")
    routes <- plotpath(ans,dests,"Selected Routes")
    grid.arrange(trace, acfplot, 
                 empty, routes, nrow=2, top="Random Source RWM Routes")
    
  }
  ans
}

test <- metro(100000, 10000, sample1, 4, plt=TRUE)

9.5 A5 Simulated Annealing

#vrp is a function to solve the vehicle routing problem using simulated annealing.
#The parameters are:
#M, the length of the chain, 
#temps, a 2 element vector with the initial and final temperatures,
#dests, the matrix containing the coordinates of the target destinations, 
#k, the number of vehicles,
#type, the cooling schedule to follow - one of ("geometric", "linear", "logarithmic","rational","root"),
#plt, a boolean specifying whether to plot the trace plot of the costs and the cooling plot
vrp <- function(M, temps = c(500,0.01), dests, k, type="geometric", plt=FALSE){
  #Organize the relevant variables
  init <- temps[1]
  final <- temps[2]
  
  #Prepare the variables that will record cooling progress, cost changes and the number of 
  #proposals accepted
  tempschedule <- rep(0,M)
  tempschedule[1] <- init
  costlist <- rep(0,M)
  acceptnum <- 0
  
  #Initialize the ordering
  cur <- data.frame(order=sample(c(2:nrow(dests))), 
                    vehicle=sample(c(1:k), size=(nrow(dests)-1), replace=TRUE))
  
  #Set the initial cost
  costlist[1] <- cost(cur, dests)
  
  #Set the cooling factor and temperature schedule according to the type
  coolingfactor <- 0
  if (type == "geometric"){
    coolingfactor <- (init/final)^(-1/M)
    tempschedule[2:M] <- sapply(c(2:M), function(x) init*coolingfactor^x)
  } else if (type == "linear"){
    coolingfactor <- (init-final)/M
    tempschedule[2:M] <- sapply(c(2:M), function(x) init-coolingfactor*x)
  } else if (type == "logarithmic"){
    coolingfactor <- (final/init)*log(1+M)
    tempschedule[2:M] <- sapply(c(2:M), function(x) coolingfactor*init/log(1+x))
  } else if (type == "rational"){
    coolingfactor <- (init/final)/(M^2)
    tempschedule[2:M] <- sapply(c(2:M), function(x) init/(coolingfactor*(x^2)))
  } else if (type == "root"){
    coolingfactor <- (init-final)/sqrt(M)
    tempschedule[2:M] <- sapply(c(2:M), function(x) init-coolingfactor*sqrt(x))
  }
  
  for (i in 2:M){
    
    #Propose a new arrangement by swapping two destinations in the ordering and allocating one 
    #destination a different vehicle
    orderswap <- sample(c(1:nrow(cur)), size=2)
    #print(orderswap)
    proposal <- cur
    proposal[orderswap[1],] <- cur[orderswap[2],]
    proposal[orderswap[2],] <- cur[orderswap[1],]
    
    vswap <- sample(c(1:nrow(proposal)), size=1)
    proposal[vswap,2] <- sample(c(1:k), size=1)
    
    #Compute the new cost and the uniform sample for the accept/reject
    propcost <- cost(proposal, dests)
    U <- runif(1)
    alpha <- exp(-(propcost - costlist[i-1])/tempschedule[i])
    if (U < alpha){
      costlist[i] <- propcost
      cur <- proposal
      acceptnum <- acceptnum + 1
    } else {
      costlist[i] <- costlist[i-1]
    }
  }
  if (plt){
    df <- data.frame(Step=c(1:M), Cost=costlist, Temperature=tempschedule)
    p1 <- ggplot(data=df, aes(x=Step, y=Cost)) + geom_line() + ggtitle("Cost Trace Plot")
    p2 <-ggplot(data=df, aes(x=Step, y=Temperature)) + geom_line() + ggtitle("Cooling Schedule")
    grid.arrange(p1,p2, nrow=1)
  }
  #Return the optimal routing
  cur
}

s1 <- vrp(20000, c(500,0.01), sample1, 4, type="geometric", plt=TRUE)

s2 <- vrp(20000, c(500,0.01), sample2, 4, type="geometric", plt=TRUE)

s3 <- vrp(20000, c(500,0.01), sample3, 4, type="geometric", plt=TRUE)

s4 <- vrp(20000, c(500,0.01), sample4, 4, type="geometric", plt=TRUE)

r1 <- plotpath(s1, sample1, title="Source: Random")
r2 <- plotpath(s2, sample2, title="Source: (0,0)")
r3 <- plotpath(s3, sample3, title="Source: (0.5,0)")
r4 <- plotpath(s4, sample4, title="Source: (0.5,0.5)")

grid.arrange(r1,r2,r3,r4, nrow=2, top="Geometric Cooling (n=15)")

s22 <- rbind(s2[1:5,],
             c(s2[13,1],3),
             s2[6:12,],
             s2[14,])
cost(s2, sample2)

## [1] 4.789146

cost(s22, sample2)

## [1] 4.761482

9.6 A6 Higher Sample Size

s2large <- vrp(100000, c(500,0.01), sample2, 4, type="geometric", plt=TRUE)

cost(s2large, sample2)

## [1] 3.438277

plotpath(s2large, sample2, title="Source: (0,0), Sample Size = 100 000")

9.7 A7 Larger Network

set.seed(23)
largemap <- networkgen(30)

l1 <- vrp(100000, c(500,0.01), largemap, 4, type="geometric", plt=TRUE)

l2 <- vrp(100000, c(500,0.01), largemap, 4, type="linear", plt=TRUE)

l3 <- vrp(100000, c(500,0.01), largemap, 4, type="rational", plt=TRUE)

l4 <- vrp(100000, c(500,0.01), largemap, 4, type="root", plt=TRUE)

lp1 <- plotpath(l1, largemap, title="Geometric Cooling")
lp2 <- plotpath(l2, largemap, title="Linear Cooling")
lp3 <- plotpath(l3, largemap, title="Rational Cooling")
lp4 <- plotpath(l4, largemap, title="Root Cooling")

cat("Geometric Cost:", cost(l1,largemap), "\n")

## Geometric Cost: 5.065856

cat("Linear Cost:", cost(l2,largemap), "\n")

## Linear Cost: 11.82136

cat("Rational Cost:", cost(l3,largemap), "\n")

## Rational Cost: 5.404348

cat("Root Cost:", cost(l4,largemap), "\n")

## Root Cost: 10.36937

grid.arrange(lp1,lp2,lp3,lp4, nrow=2, top="Different Cooling Schedules (n=30)")

largegeo <- rep(0,10)
opt <- 100000
ans <- NA

for (i in 1:10){
  cur <-vrp(100000, c(500,0.01), largemap, 4, type="geometric", plt=FALSE)
  largegeo[i] <- cost(cur, largemap)
  if (largegeo[i] < opt) {
    ans <- cur
  }
}
largegeo

##  [1] 4.946364 5.622043 5.506219 5.556710 5.550898 5.465494 5.421263
##  [8] 4.957598 5.593201 5.210520

plotpath(ans, largemap, title="Geometric Cooling (Best of 10 Runs, n=30)")

9.8 A8 Time Based Problem

set.seed(7)
#Define a function that calculates the travel time
traveltime <- function(a,b,multiplier){
  avg <- sqrt(((a[1]-b[1])^2) + ((a[2]-b[2])^2))*multiplier
  #For practicality purposes, set a floor on the traveltime
  max(0,rnorm(1, mean=avg, sd=(avg/4)))
}

#Define a function to compute the cost of a given allocation and ordering
#Note that x is a (n-1)x2 matrix where the first column contains a number from 
#2 to n with the name of the destination number and the second column is which 
#of the k delivery vehicles is assigned to the destination. The cost is computed 
#as the sum of the distances between destinations assigned to a specific vehicle
#in the order that they appear. 
#Note that loc is the location matrix of the destinations that characterizes our
#network. The first row corresponds to the source node which will all vehicles 
#must originate from and all vehicles must return to
cost2 <- function(x, loc){
  k <- max(x[,2])
  #print(x)
  d <- rep(0,k)
  mult <- rep(1,k)
  cur <- matrix(data=rep(loc[1,],k), nrow=k, ncol=2, byrow=T)
  for (i in 1:nrow(x)){
    dest <- x[i,1]
    vehicle <- x[i,2]
    ttime <- traveltime(loc[dest,],cur[vehicle,],mult[vehicle])
    d[vehicle] <- d[vehicle] + ttime
    mult[vehicle] <- mult[vehicle] * ttime/dist(loc[dest,],cur[vehicle,])
    cur[vehicle,] <- loc[dest,]
  }
  #return each vehicle to the source
  for (i in 1:k){
    d[i] <- d[i] + traveltime(loc[1,], cur[i,], mult[i])
  }
  #Note that we are now interested in the longest travel time and not the sum
  max(d)
}

#vrp2 is a function to solve the vehicle routing problem using simulated annealing.
#This version minimizes the maximum vehicle route time. The parameters are:
#M, the length of the chain, 
#temps, a 2 element vector with the initial and final temperatures,
#dests, the matrix containing the coordinates of the target destinations, 
#k, the number of vehicles,
#N, the number of samples of the travel time to compute the expected cost
#type, the cooling schedule to follow - one of ("geometric", "linear", "logarithmic","rational","root"),
#plt, a boolean specifying whether to plot the trace plot of the costs and the cooling plot
vrp2 <- function(M, temps = c(500,0.01), dests, k, N, type="geometric", plt=FALSE){
  #Organize the relevant variables
  init <- temps[1]
  final <- temps[2]
  
  #Prepare the variables that will record cooling progress, cost changes and the number of 
  #proposals accepted
  tempschedule <- rep(0,M)
  tempschedule[1] <- init
  costlist <- rep(0,M)
  acceptnum <- 0
  
  #Initialize the ordering
  cur <- data.frame(order=sample(c(2:nrow(dests))), 
                    vehicle=sample(c(1:k), size=(nrow(dests)-1), replace=TRUE))
  
  #Set the initial cost
  costlist[1] <- mean(replicate(N, cost2(cur, dests)))
  
  #Set the cooling factor and temperature schedule according to the type
  coolingfactor <- 0
  if (type == "geometric"){
    coolingfactor <- (init/final)^(-1/M)
    tempschedule[2:M] <- sapply(c(2:M), function(x) init*coolingfactor^x)
  } else if (type == "linear"){
    coolingfactor <- (init-final)/M
    tempschedule[2:M] <- sapply(c(2:M), function(x) init-coolingfactor*x)
  } else if (type == "logarithmic"){
    coolingfactor <- (final/init)*log(1+M)
    tempschedule[2:M] <- sapply(c(2:M), function(x) coolingfactor*init/log(1+x))
  } else if (type == "rational"){
    coolingfactor <- (init/final)/(M^2)
    tempschedule[2:M] <- sapply(c(2:M), function(x) init/(coolingfactor*(x^2)))
  } else if (type == "root"){
    coolingfactor <- (init-final)/sqrt(M)
    tempschedule[2:M] <- sapply(c(2:M), function(x) init-coolingfactor*sqrt(x))
  }
  
  for (i in 2:M){
    
    #Propose a new arrangement by swapping two destinations in the ordering and allocating one 
    #destination a different vehicle
    orderswap <- sample(c(1:nrow(cur)), size=2)
    #print(orderswap)
    proposal <- cur
    proposal[orderswap[1],] <- cur[orderswap[2],]
    proposal[orderswap[2],] <- cur[orderswap[1],]
    
    vswap <- sample(c(1:nrow(proposal)), size=1)
    proposal[vswap,2] <- sample(c(1:k), size=1)
    
    #Compute the new cost and the uniform sample for the accept/reject
    propcost <- mean(replicate(N,cost2(proposal, dests)))
    U <- runif(1)
    alpha <- exp(-(propcost - costlist[i-1])/tempschedule[i])
    if (U < alpha){
      costlist[i] <- propcost
      cur <- proposal
      acceptnum <- acceptnum + 1
    } else {
      costlist[i] <- costlist[i-1]
    }
  }
  if (plt){
    df <- data.frame(Step=c(1:M), Cost=costlist, Temperature=tempschedule)
    p1 <- ggplot(data=df, aes(x=Step, y=Cost)) + geom_line() + ggtitle("Cost Trace Plot")
    p2 <-ggplot(data=df, aes(x=Step, y=Temperature)) + geom_line() + ggtitle("Cooling Schedule")
    grid.arrange(p1,p2, nrow=1)
  }
  #Return the optimal routing
  cur
}

t1 <- vrp2(20000, c(500,0.01), sample1, 4, 100, type="geometric", plt=TRUE)

tp1 <- plotpath(s1, sample1, title="Distance Solution")
tp2 <- plotpath(t1, sample1, title="Time Solution")

grid.arrange(tp1,tp2, nrow=1, top="Two Problems, Geometric Cooling (n=15)")

10 References

The following sources offered insight into the vehicle routing problem and to applying simulated annealing to combinatorial optimization problems:

• http://www.cs.cmu.edu/afs/cs.cmu.edu/project/learn-43/lib/photoz/.g/web/glossary/anneal.html
• https://developers.google.com/optimization/routing/vrp
• https://epubs.siam.org/doi/pdf/10.1137/1.9780898718515.fm