A Monarch Butterfly Optimization for the Dynamic Vehicle Routing Problem

Abstract

The dynamic vehicle routing problem (DVRP) is a variant of the Vehicle Routing Problem (VRP) in which customers appear dynamically. The objective is to determine a set of routes that minimizes the total travel distance. In this paper, we propose a monarch butterfly optimization (MBO) algorithm to solve DVRPs, utilizing a greedy strategy. Both migration operation and the butterfly adjusting operator only accept the offspring of butterfly individuals that have better fitness than their parents. To improve performance, a later perturbation procedure is implemented, to maintain a balance between global diversification and local intensification. The computational results indicate that the proposed technique outperforms the existing approaches in the literature for average performance by at least 9.38%. In addition, 12 new best solutions were found. This shows that this proposed technique consistently produces high-quality solutions and outperforms other published heuristics for the DVRP.

1. Introduction

The dynamic vehicle routing problem (DVRP) is a hard combinatorial optimization problem which is advanced by information and communication technologies that allow information to be obtained and processed in real time, typically used in distribution logistics and transportation systems. This system facilitates quick updating of transportation system plans in unexpected or uncertain events, for example if roads between two customers are blocked off, customers can modify their orders, or when the travel time for some routes is increased due to bad weather conditions or traffic congestion, etc. On the other hand, customers in the e-market place increasingly expect quicker and more flexible fulfilment of their transportation requests. In this context, the DVRP is becoming increasingly interesting and important.

The DVRP is a variant of the VRP—the well-known NP-hard problem which is used in general cases. Because VRPs can be regarded as special cases of dynamic VRPs, dynamic VRPs are at least as hard as VRPs; therefore, research efforts mainly focus on metaheuristics and intelligent optimization algorithms. For example, Ismail [1] designed a hybrid GA and Tabu Search heuristic to solve the DVRP. Oliveira et al. [2] applied ant colony system (ACS) techniques to solve the DVRP with time window constraints and a capacitated fleet. Belfiore and Yoshizaki [3] proposed a scatter search (SS) approach to deal with real-life heterogeneous fleet vehicle routing with time windows and split deliveries in Brazil. Chen and Xu [4] proposed a dynamic column generation algorithm for the DVRP with time windows. de Armas et al. [5] presented a general variable neighborhood search algorithm to solve a real-world application of the DVRP.

In this paper, a monarch butterfly optimization (MBO) approach is proposed for the DVRP. Although MBO has been successfully applied to a great variety of hard combinatorial optimization problems ([6,7,8,9]), this paper, as far as we know, proposes the first MBO-based algorithm for the DVRP.

Biologically-inspired computation is a field focusing on the development of computational tools, modeled on the principles that exist in natural systems. The results presented in recent publications show that bio-inspired approaches are now highly competitive with other state-of-the-art heuristics.

MBO was proposed by Wang et al. [10] in 2015. Preliminary studies indicated that MBO is a very competitive metaheuristic algorithm [10], when compared with ant colony optimization (ACO) [11], biogeography-based optimization (BBO) [12] and differential evolution (DE) [13], is very easy to implement because it only needs to fine-tune migration and adjusting operators. In addition, MBO is very capable of finding the shortest paths [9]. However, MBO may fail to reach optimal performance on average and standard values in certain test cases [10].

This paper considers a version of the DVRP in which some of the customers’ locations and associated demands are unknown at the start of the working day and arrive gradually as time passes. All vehicles start from the depot to serve yesterday’s remaining customers after opening time and must return to the depot by the closing time. A similar problem was first introduced by Kilby et al. [14]. In order to overcome the MBO shortcoming, a greedy algorithm is incorporated into the migration operation and butterfly adjusting operator, and a later perturbation procedure is applied after the butterfly adjusting operator, to maintain a balance between global diversification and local intensification. To the best of our knowledge, this is the first MBO implementation for the DVRP. It is tested using data sets introduced in Kilby et al. [14] and Montemanni et al. [15], and compared to other well-known meta-heuristics. As a result, we obtain an average gap of −9.38%, provide 12 new best solutions from 21 instances, and improve the overall performance for this set of instances.

The remainder of this paper is organized as follows: In Section 2, we define the DVRP model tackled here and present a literature review of the existing papers that deal with the DVRP. Section 3 contains a description of our algorithm. In Section 4, we include the results of the experiment and the corresponding analysis to test our proposed algorithm. We summarize the major conclusions of this article in Section 5.

2. Dynamic Vehicle Routing Problem (DVRP)

2.1. Problem Description

To allow better understanding of the DVRP, Figure 1 shows a simple example of a dynamic vehicle routing problem. In the example, two uncapacitated vehicles must service both advance- and immediate-request customers. The static customers are represented by white nodes, while those that are immediate requests are depicted by black nodes. The dashed lines represent the two routes that the dispatcher has planned prior to the vehicles leaving the depot (Figure 1A). The two dash-dotted lines indicate the vehicle positions at the time the dynamic requests appear. Ideally, the dispatcher routes should be adjusted to fulfill the dynamic customers’ requests (Figure 1B). Finally, the vehicles deliver the demands of all customers and return to the depot (Figure 1C).

The study problem concerns routing a fleet of capacitated vehicles in real time to deliver goods ordered by a partially-known set of customers. There is a fleet of K homogenous vehicles, each with an identical capacity (Q), and a set of n customers {1, 2, …, n} is considered. The problem can be modeled by an undirected graph, G = (V, E), where V is a vertex set and E is an edge set (V = {0}∪S∪D, and E = {(i, j) | i, j ∈ V, i < j}). Vertex 0 denotes a depot, which has a working time from e₀ to l₀ (0 ≤ e₀ < l₀). C is a set of customers (requests) to be serviced. For each customer, (I ∈ S∪D), there is a location (x_i, y_i), a demand of goods (d_i), an appearance moment (a_i (e₀ ≤ a_i ≤ l₀)), and a service time (s_i). At the same time, each edge is defined by a non-negative distance (c_ij) and the travel time (t_ij) between all customers’ locations. For DVRPs, there are two sets of customers: S and D. S is composed of static customers, who are known from an earlier working day, and D includes all dynamic customers, who appear with time. A solution, S, is composed of a set of routes satisfying all customers’ demands once. Each route (R_k ∈ S) starts and ends at the depot, and does not exceed the vehicle capacity (Q). A route consists of a sequence of customers to visit (R_k = {0, …, j, …, 0}), where j ∈ V denotes the j-th customer. For each customer (I ∈ R_k), w_i is the vehicle’s waiting time, and the vehicle’s arriving time (b_i) should satisfy b_i = b_i−1 + s_i−1 + w_i−1 + t_i−1,i (b₀ = e₀), when we reversely calculate the latest arrival time with z_i = z_i+1 − t_i+1,i − s_i (z_j+1 = l₀).

Similar to [14], we also divide a working day (T) into time slices (n_ts), each with an equal length of time (T/n_ts, where T = l₀ − e₀), where TS = {TS₁, TS₂, …, TS_ts}. Each time slice represents a partial static VRP, where TS₁ only includes static customers in S. The equation TS\TS₁ considers the customers received at the previous time slice and those who have not been visited yet.

Hence, two decision variables, x_ijk and y_ik, are defined as follows:

x_ijk = 1 if the edge (i, j) is visited by vehicle k and 0 in all other cases.

y_ik = 1 if customer i is visited by vehicle k and 0 in all other cases.

The general mathematical model for DVRP is described as follows:

$$\min {\displaystyle \sum }_{(i,j)\in E}{\displaystyle \sum }_{k\in K}{c}_{ij}·x_{ijk}$$

(1)

s.t.

$${\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{d}_i·y_{ik}\le Q, \forall k\in K$$

(2)

$${\displaystyle \sum }_{k\in K}{y}_{ik}=1,\forall i\in V$$

(3)

$${\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{x}_{ijk}=y_{jk},\forall k\in K,\forall j\in V$$

(4)

$${\displaystyle \sum }_{j\in V\backslash \left\{0\right\}}{x}_{ijk}=y_{ik},\forall k\in K,\forall i\in V$$

(5)

$${\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{x}_{ijk}={\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{x}_{jik},\forall k\in K,\forall j\in V$$

(6)

$${\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{\displaystyle \sum }_{k\in K}{x}_{0ik}={\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{\displaystyle \sum }_{k\in K}{x}_{i0k}\le \left|K\right|$$

(7)

$${\displaystyle \sum }_{(i, j)\in E}{x}_{ijk}·t_{ij}+{\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{y}_{ik}·s_i+{\displaystyle \sum }_{i\in V\backslash \left\{0\right\}}{y}_{ik}·w_i\le l_0,\forall k\in K$$

(8)

The objective (1) is to minimize the total distance travelled during each time slice. Expression (2) ensures that during each time slice, each vehicle does not exceed its capacity. Expressions (3)–(6) ensure that each customer, apart from the depot, is visited once, and only once, by one vehicle, during each time slice. The Expression (7) ensures that all vehicles start from the depot during initialization and return to the depot at the end of the route. Meanwhile, the number of vehicles out in use cannot exceed the total quantity at the depot. Expression (8) is the travel time constraints, because every vehicle must return to the depot before the closing time.

2.2. Related Work

The Dynamic Vehicle Routing Problem (DVRP) with capacity and time duration constraints was introduced by Kilby et al. [14] and further refined by Montemanni et al. [15]. These authors proposed some benchmark instances for the DVRP and presented a study on how the degree-of-dynamism affects the final travel costs. Montemanni et al. [15], who extended Kilby et al.’s work [14], considered a DVRP as an extension to the standard VRP, by decomposing a DVRP as a sequence of static VRPs, and then solving them using an ACS algorithm. Other optimization techniques have been applied. Khouadjia et al. [16] presented particle swarm optimization. Yang et al. [17] proposed a hybrid large neighborhood search. In a recent survey, Pillac et al. [18] classified routing problems from the perspective of information quality and evolution. They introduced the notion of degree of dynamism, and presented a comprehensive review of applications and solution methods for dynamic vehicle routing problems. Bekta et al. [19] provided another survey in this area, which provided a deeper and more detailed analysis. Last but not least, Psaraftis et al. [20] shed more light into work in this area, over more than 3 decades, by developing a taxonomy of DVRP papers according to 11 criteria.

Along with advances in technology, improvements to problem-solving methods have been made in order to interact with dynamic environments, and to combine fast responses with quality solutions. Based on the quality of information, Pillac et al. [18] classified DVRPs into two categories: deterministic and stochastic. The present work falls under the deterministic category, for which approaches are based on periodic or continuous re-optimization.

Periodic re-optimization approaches typically commence at the beginning of the day, with an initial optimization that produces an initial set of routes; then the solution is re-optimized either whenever the available information changes, or at a fixed re-optimization interval. Chen and Xu [4] proposed a dynamic, column generation approach for solving a DVRP, with hard time windows, in which all requests need to be serviced. The authors used the concept of decision epochs over the planning horizon, which indicate the moments of the day when the re-optimization process is executed. Montemanni et al. [15] employed an ACS to solve the dynamic VRP by dividing the overall planning horizon into periods (time slices), as in Kilby et al. [14]. During each time slice, a static problem is solved by considering all requests known at the beginning of this time slice. A similar approach was also employed by Elhassania et al. [21], Euchi et al. [22] and Mańdziuk and Żychowski [23].

In contrast, continuous re-optimization approaches commence at the beginning of operations with an initial set of routes, and the vehicles are informed only about their next destination. It keeps high quality solutions within an adaptive memory. A decision procedure is used to update the solutions in the adaptive memory whenever the available information is updated. Updates typically occur due to either a new dynamic customer request service, or the vehicle reaches a customer. Gendreau et al. [24] were the first to employ continuous re-optimization. These authors proposed a Tabu Search heuristic to address the DVRP arising in a long-distance courier service, in which time windows could be violated at some cost. The approach keeps a pool of good solutions in its adaptive memory, which is used to generate initial solutions for a parallel Tabu Search. Whenever a new customer request arrives, it is checked against all the solutions from the adaptive memory to decide whether it should be accepted or rejected. A fast local search procedure is then applied to select the best solution. When the vehicle reaches a customer, the solution is also updated in the adaptive memory and the next destination of the vehicle is identified, based on the best solution. In another study, Bent and Van Hentenryck [25] generalized this framework and introduced the concept of the multiple plan approach to solve the DVRP. The method attempts to continuously generate different solutions based on static and known dynamic customers. A solution pool (routing plans) is used to generate an appropriate solution. When a new customer arrives, a procedure checks if they can be served or not. If they can be served, the customer request is inserted into the solution pool and incompatible solutions are discarded. The solution pool is updated during each event in order to ensure that all solutions are consistent with the current state of the system.

3. The Proposed Algorithm for DVRP

We used the periodic re-optimization strategy to solve the DVRP. This strategy has been used in many DVRP works, such as Kilby et al. [14] and Montemanni et al. [15]. The basic idea of the strategy is to divide a DVRP into a series of static VRPs for every time-slice. The initial states of the static VRPs are different from standard static VRPs—that is, vehicles start from their position in the last slice.

The modified MBO algorithm we propose includes three important parts: (1) an insert heuristic to place new requests into the current solution at the end of each time slice; (2) a greedy strategy, incorporated into the standard MBO, which accepts only the monarch butterfly individuals that have better fitness than their parents, and (3) a subsequent perturbation (2-opt*) to increase the search diversity and the chance of escaping local optima. Next, we discuss the proposed algorithm in detail.

3.1. Solution Representation

The representation used in our work is dedicated to the DVRP. We propose a simple and intuitionistic natural number encode which expresses the route of k vehicles to the n customers they serve. As previously described, 0 represents the depot, and 1, 2, …, n refers to each customer. As the total number of vehicles is k, the maximum route distribution is also k. Each route begins with the depot, and ends with the depot. In order to reflect the distribution route in coding, we increase k − 1 to represent the virtual depots, respectively expressed as n + 1, …, n + k − 1. A non-duplication of a random arrangement of the natural numbers—1, 2, …, n, n + 1, …, n + k − 1—constitutes an individual.

An example of butterfly representation and its decoding are presented in Figure 2. There are seven customers and three vehicles to complete the task of distribution; we use a non-duplication of a random arrangement of the natural numbers—1, 2, …, 9 (among them, 8 and 9 are virtual depots)—to represent the distribution routes. In Figure 1A, an individual (7 4 6 8 5 2 9 1 3) contains three distribution routes (route 1: 0-4-6-0, route 2: 0-5-2-0, and route 3: 0-1-3-0), while in Figure 2B, an individual (7 4 6 5 8 9 2 3 1) corresponds to the program of the distribution path, route 1: 0-7-4-6-5-0, route 2: 0-2-3-1-0—a total of two distribution routes.

3.2. Initial Population

In accordance with Kilby et al. [14] and many similar works, a working day is split into n_ts equal-length time slices, and in each time slice a static VRP is solved, with the aim of servicing the known customers so far. We postpone the arrival of a new request to the end of the current time slice and optimize the next time slice. So, a cut-off time (T_co), is considered, as introduced in Montemanni et al. [15], which postpones orders after T_co to the next working day. In addition, an advanced commitment time, (T_ac), is also employed to allow a driver to respond to new orders prior to the time of processing the order itself. The first partial static VRP created at the beginning of the working day (T₁) consists of all static customers which are left over from the previous working day, as defined by the T_co parameter. We apply a random approach to generate an initial feasible solution (S). Initially, S consists of k − 1 virtual depots, and then we try to update those routes by inserting customers randomly. This process is repeated until all customers have been assigned to a route. Each individual is defined as a random permutation of these known customers and existing routes. Hence, at each slice, all butterflies have the same length and contain the same numbers (customers’ and vehicles’ identifiers). All butterflies will be verified against the capacity and feasibility rules described above.

The next static problem considers all orders received during the previous time slice as well as those which have not been visited by drivers yet. In our simulation, each vehicle (k) starts from the location of the last visited customer (i), with a departure time (dt_i) corresponding to the end of the service time for this customer (dt_i = b_i + s_i), and with a remaining capacity equal to the capacity left after serving all previously visited customers.

A decision to wait or go occurs when a vehicle finishes servicing a customer. The algorithm must choose between keeping the vehicle waiting at its current location or beginning to move towards the next customer. Assume that in l-th time slice (TS_l), the vehicle (k) has finished serving the customer (i) and is ready to serve the next customer (I + 1) in the route (R_k). If it satisfies the inequalities (9), the vehicle will wait at the customer until T_l (T_l = l*T/n_ts), and the wait time (w_i) at the customer is defined in Equation (10)

$$d{t}_i<T_l\wedge b_{i+1}>T_l\wedge T_l+t_{i,i+1}<z_i$$

(9)

$$w_i=T_l-d{t}_i$$

(10)

3.3. Monarch Butterfly Optimization (MBO)

MBO is a new naturally-inspired metaheuristic algorithm, which is proposed by Wang et al. [10] in 2015. It is inspired by the behavior of the monarch butterfly during migration. In MBO, all the monarch butterfly individuals are idealized and located in two lands only: northern USA and southern Canada (Land 1) and Mexico (Land 2). Accordingly, the positions of the monarch butterflies are updated in two ways: the migration operator and the butterfly adjusting operator. In the following subsections, we discuss our modified operations for solving the DVRPs.

3.3.1. Migration Operator

The migration operator is designed to update the monarch butterfly migration between Land 1 and Land 2, on which monarch butterflies make up subpopulations 1 and 2 respectively. At first, the number of monarch butterflies in Lands 1 and 2 can be considered as NP1 = ceil(p*NP) and NP2 = NP − NP1, respectively, where NP is the total number of monarch butterflies, p is the ratio of monarch butterflies in Land 1, and ceil(x) rounds x to the nearest integer greater than or equal to x. Accordingly, migration operator can be formulated as

$$x_{i,k}^{t + 1}=\{\begin{array}{c}{x}_{r_1,k}^t |r\le p\\ x_{r_2,k}^t |r>p\end{array}$$

(11)

where $x_{i,k}^{t+1}$. represents the kth element of x_i at generation t + 1. Similarly, $x_{r1,k}^t$. indicates the kth element of x_r1 at generation t, and $x_{r2,k}^t$. represents the kth element of x_r2 at generation t. The current generation number is represented by t. Monarch butterflies (r1 and r2) are selected randomly from subpopulation 1 and subpopulation 2. The condition variable (r) is calculated as

r = rand * peri

(12)

where peri indicates the migration period and is set to 1.2 in the basic MBO method [10] and rand is a random number derived from a uniform distribution.

3.3.2. Butterfly Adjusting Operator

This operator is used to update the positions of the monarch butterflies in subpopulation 2. It can be updated as follows:

$$x_{j,k}^{t + 1}=\{\begin{array}{c}{x}_{best,k}^t |rand\le p\\ x_{r3,k}^t |rand>p\end{array}$$

(13)

where $x_{j,k}^{t+1}$ represents the kth element of x_j at generation t + 1; $x_{best,k}^t$ indicates the kth element of x_best at generation t, which represents the best location of the monarch butterflies in Lands 1 and 2; and $x_{r3,k}^t$ represents the kth element of x_r3 at generation t; the monarch butterfly r3 is selected randomly from subpopulation 2. If rand > p, there has another step. The position of the butterfly is further updated using Levy flight, if rand > BAR:

$$x_{j,k}^{t + 1}=x_{j,k}^{t + 1}+\alpha \times (dx-0.5)$$

(14)

where the variable BAR is the butterfly adjusting rate. If BAR is smaller than a random value, the kth element of x_j at generation t + 1 is updated, where α is the weighting factor, as shown in Equation (15).

$$\mathsf{\alpha }=S_{max}/t^2$$

(15)

where S_max is the maximum walk step. In Equation (14), dx is the walk step of butterfly j that can be calculated by Levy flight.

$$dx=\mathrm{Levy}(x_j^t)$$

(16)

3.4. Greedy Acceptance

In the basic MBO method, all the newly-generated butterflies are accepted, and are kept in the next generation, while a greedy strategy is used to only accept the butterfly individuals that have better fitness in our algorithm. This greedy strategy can be formulated as follows:

$$x_{i,new}^{t + 1}=\{\begin{array}{cc}\hfill x_i^{t + 1},& |f(x_i^{t + 1})<f(x_i^t)\hfill \\ \hfill x_i^t,& |otherwise\hfill \end{array}$$

(17)

where $x_{i,new}^{t+1}$ is a newly-generated butterfly individual for the next generation. The terms $f(x_i^{t+1})$ and $f(x_i^t)$ are the fitness levels of butterfly $x_i^{t+1}$ and $x_i^t$, respectively.

3.5. Later Perturbation

In order to enhance the ability of escaping from a local optimum, we employ a perturbation method to increase the search diversity. We randomly select an individual from the population, and then apply the later perturbation to improve the selected solution. If the resultant solution is feasible, it replaces the worst individual in the population.

The later perturbation uses two improving procedures, as illustrated in Figure 3. First, we use a 2-opt heuristic, which inverts sequences of one single route; second, we use a 2-opt* heuristic, which exchanges the ending segments of two different routes. A first improvement strategy is used, which means that the local search restarts immediately when an improvement is found. Each local search runs until no more improvements are possible. The 2-opt local search runs first, its output solution being the input to the 2-opt* local search. The local search terminates when no 2-opt* move exists that can further improve the solution.

3.6. Algorithm Structure

The modified MBO is shown in Algorithm 1. At first, all parameters are initialized. The variable n_ts divides the working day (T) into equal time slices. For each time slice, the following procedure is repeated: dynamic customers arriving at time slice TS_l are added to the request pool (U) (Line 4). Specially, for TS₁, U only stores all static customers. For other time slices, the location of each vehicle is updated by Equations (9) and (10) with the current best solution (S). Then, the unvisited customers from S are removed and inserted into U (Lines 5–8). The population (P) of NP (=|U|) individuals is initialized. Next, the basic MBO is repeated. First, the population is divided into two subpopulations, according to their fitness values, in which subpopulation 1 includes NP1 better individuals and subpopulation 2 includes the rest. Afterwards, an offspring of subpopulation 1 is generated, using the migration operator (Lines 14–19). Then, individuals of subpopulation 2 are updated by the butterfly adjusting operator (Lines 20–25). Note that each two operators only accept the butterfly individuals that have better fitness (Lines 18 and 24). Next, offspring of the two subpopulations are regrouped into one population. Finally, later perturbations are applied to local search for any one individual. The current best individual is used as an input in the next time slice. The above procedure is repeated until the termination condition is satisfied.

Algorithm 1 MBODVRP

1. Initializing all the parameters  2. S ← ϕ, U ←ϕ, l ←1  3. while (l <= nts) do  4.   Add the customers received at the time slice TSl to U.  5.   if (S is not null) then  6.     Update the location of each vehicle by Equations (9) and (10)  7.     Remove the unvisited customers from S and then insert into U.  8.   end if  9.   Initializing the population P of NP individuals randomly. 10.   Evaluate butterfly population 11.   while (termination criterion is not reached) do 12.     Sort all monarch butterfly individuals in non-ascending order based on fitness. 13.     Divide butterflies into two subpopulations. 14.     for (each monarch butterfly i in subpopulation 1) do 15.       for (each element k in i-th butterfly) do 16.         Generate k-th element by Equations (11) and (12). 17.       end for k 18.       Update i-th butterfly by greedy strategy as Equation (17) 19.     end for i 20.     for (each monarch butterfly j in subpopulation 2) do 21.       for (each element k in j-th butterfly) do 22.         Generate k-th element by Equations (13)–(16). 23.       end for k 24.       Update j-th butterfly by greedy strategy as Equations (17) 25.     end for j 26.     Combine the two newly-generated subpopulations into one population; 27.     Evaluate butterfly population 28.     Apply later perturbation 29.     S ← the current best individual in the population. 30.   end while 31.   l ← l + 1, U ←ϕ 32. end while 33. return S.

4. Experimental Results

In this section, some computational results are presented, in order to evaluate the performance of the improved MBO described in Section 3. We carried out extensive experiments on relevant benchmarks. We first introduce the benchmarks instances and the parameter settings for the experiments. We then show the results, along with a comparative analysis, on several public data sets derived from Montemanni et al. [15] and Kilby et al. [14].

In order to evaluate the performance of the proposed algorithm, we implemented the algorithm in Visual C# under Windows Server 2012 R2 Standard on a PC with Intel^® Xeon^® E5-2620 @ 2.40 GHz CPU and 16 GB RAM.

4.1. Benchmarks Description

Our experimental results are based on the benchmarks proposed by Kilby et al. [14] and extended by Monemanni et al. [15]. They were derived from publicly available VRP benchmark data from three separate VRP sources (datasets). The first dataset consists of 12 instances, varying from 75 to 150 customers, taken from Taillard et al. [26]. The second dataset consists of seven instances, with sizes varying from 50 to 199 customers, derived from Christofides et al. [27]. The third set consists of two instances, with sizes varying from 71 to 134 customers, derived from Fisher et al. [28]. The number of customers can be inferred from the name of each instance. The service area may consist of uniformly distributed customers, clustered customers, or a combination of both. In order to obtain dynamic problems, Kilby et al. [14] added the concept of length of the working day into to these problem, and gave each instance an appearance time which signifies when the order becomes available (during the working day) and a duration, namely the time required to perform an order once it reaches the customer. They also fixed the number of available vehicles to 50 for each problem. More details can be found in Kilby et al. [14].

In order to compare our results with existing algorithms, a number of parameters needed to be adjusted. The first was n_ts, the number of time slices in the optimization. In Montemanni et al. [15], it was found that setting the parameter to n_ts = 25 yielded the best trade-off between the objective value and computational cost. Secondly, the advanced commitment time (T_ac) was set at 0. Finally, the cutoff time (T_co), was set to 0.5T, where T is the total length of the working day. In particular, for each time slice, we limited computation time to 30 s. Then, T equals 30 × 25 = 750 s for the entire simulation.

4.2. Comparison with the Literature for DVRP

A comparison for the quality of the solution, in terms of minimizing travel distances, was done between our approach and other metaheuristics proposed previously in the literature. In order to obtain significant results, our approach was executed five times for each instance.

Table 1. Comparison of the best Monarch Butterfly Optimization (MBO) results with the hitherto best-known solutions.

	MBO	BKS	Algorithm	RE (%)	TMBO (s)
c50	570.61	551.95	AAC-2OPT	3.38	42.44
c75	897.16	962.79	GA	−6.82	83.01
c100	915.27	961.10	GA2007	−4.77	80.80
c100b	819.60	800.93	AAC-2OPT	2.33	101.77
c120	1070.18	1049.47	AAC-2OPT	1.97	122.56
c150	1118.03	1318.22	TS	−15.19	250.99
c199	1394.74	1640.40	DAPSO	−14.98	470.15
f71	271.43	279.52	DAPSO	−2.89	51.53
f134	11759.06	13015.56	AAC-2OPT	−9.77	205.78
tai75a	1748.71	1755.33	AAC-2OPT	−0.38	81.75
tai75b	1371.90	1306.47	AAC-2OPT	5.01	90.42
tai75c	1464.37	1406.27	TS	4.13	77.41
tai75d	1419.00	1334.67	AAC-2OPT	6.32	83.13
tai100a	2185.49	2194.93	AAC-2OPT	−0.43	113.82
tai100b	2057.80	2126.09	AAC-2OPT	−3.21	108.28
tai100c	1463.42	1490.58	VNS	−1.82	131.28
tai100d	1704.35	1834.60	GA2007	−7.10	119.77
tai150a	3560.28	2999.27	AAC-2OPT	18.70	336.75
tai150b	2966.19	2846.28	AAC-2OPT	4.21	333.55
tai150c	2528.04	2612.68	GA2007	−3.24	353.91
tai150d	2980.44	2950.61	GA2007	1.01	292.68
Average	2107.91	2163.70		−1.12	168.18

and

Table 2. Comparison of the average MBO results with the hitherto best known solutions.

	MBO	BKS	Algorithm	RE (%)	TMBO (s)
c50	590.73	570.89	AAC-2OPT	3.48	44.20
c75	909.56	1013.45	TS 2007	−10.25	84.50
c100	930.38	987.59	GA 2007	−5.79	82.43
c100b	839.93	841.44	AAC-2OPT	−0.18	103.90
c120	1112.58	1153.29	AAC-2OPT	−3.53	126.08
c150	1135.85	1386.93	GA 2007	−18.10	262.19
c199	1414.40	1758.51	AAC-2OPT	−19.57	492.80
f71	277.00	306.33	TS 2007	−9.57	52.48
f134	11853.29	15528.81	AAC-2OPT	−23.87	215.09
tai75a	1799.15	1782.91	AAC-2OPT	0.91	83.63
tai75b	1385.45	1452.26	AAC-2OPT	−4.60	93.61
tai75c	1483.41	1441.91	AAC-2OPT	2.88	81.09
tai75d	1435.62	1422.27	AAC-2OPT	0.94	86.42
tai100a	2212.08	2232.71	AAC-2OPT	−0.92	116.53
tai100b	2105.15	2182.61	AAC-2OPT	−3.55	111.98
tai100c	1474.32	1541.25	GA 2014	−4.04	137.05
tai100d	1722.88	1912.43	AAC-2OPT	−10.78	123.70
tai150a	3539.61	3185.727	AAC-2OPT	8.78	375.45
tai150b	3038.38	2880.57	AAC-2OPT	5.48	365.32
tai150c	2641.62	2743.55	AAC-2OPT	−3.72	377.15
tai150d	3047.30	3045.16	GA 2007	0.07	322.49
Average	2140.41	2350.98		−9.38	178.01

show the best computational results by our MBO algorithm and the average results over five runs, respectively. There are six columns in Table 1 and Table 2: the instance, our MBO algorithm, the best known solutions (BKS), the algorithm to find the BKS, the relative error (RE) (i.e., RE = (MBO − BKS)/BKS × 100%) and the computational time by MBO (note that we cannot compare with other algorithms as there is no such data available in the public domain or in the literature). There are other best and average results in Table 1 and Table 2 provided by some heuristic algorithms as follows: GA2007 and TS by Hanshar and Ombuki-Berman [29], DAPSO and VNS by [16], GA2014 by [21], HLNS by [17] and AAC-2OPT by [22].

As you can see, our approach is very competitive. Table 1 shows that MBO found 12 out of 21 new best solutions, compared to the BKS. The average distance of all instances by the MBO was 2096.55 which is lower than the BKS’s average cost of 2163.70. It made an improvement of 1.12% over the BKS. Regarding the computation time of MBO, the longest was no more than 470.15 s (much lower than the time limit of 750 s that we set), and the shortest was 42.44 s, while the average time for each instance was 168.18 s. In particular, the computation time for each time slice was only 6.73 s on average. This indicates that MBO has good computational capability.

From Table 2, our average results were superior to GRASP in almost all instances, and better than BKS in 14 instances. All these results allow us to say that our MBO is effective and shows the viability to generate very high quality solutions for the DVRP.

In the next experiment, we compared the AAC-2OPT, DAPSO, and GA2007 with our MBO. The percentage of deviation from the known optimal solution concerning problems can be seen as a chart in Figure 4. It shows that MBO is superior to other methods. The deviation of MBO is less than 5% for all instances except tail75d and tail150a.

In this experiment, we analyzed the results in detail. Figure 5 shows, for each data set, the dispersion diagrams of the solutions found by our approach. Note that the solution values are presented as the relative error from the BKS.

In most cases, our MBO gets a relative error of no greater than 10% from the BKS, and central values tend to be between 0 and 5%. It is interesting to note that the quality and stability of the solutions are not size dependent. In particular, f134, c75, c100, c150, c199, and t100d have instances with less than 5% relative error. This means those instances are better than the BKS.

5. Conclusions

This paper studied the DVRP, a method where the customers appear dynamically, and which requires a combination of the vehicles’ current states and the customers’ current information to re-plan routes. In this paper, an MBO approach, called MBO–DVRP, is proposed for the problem. The approach is based on a Monarch Butterfly Optimization algorithm and perturbed a 2-opt(*) procedure. To the best of our knowledge, this is the first time that MBO has been proposed and evaluated by the DVRP. The experiments conducted on the benchmark for the DVRP clearly demonstrate the efficiency and the competitiveness of our approach compared to the existing methods in the literature.