Using Monarch Butterfly Optimization to Solve the Emergency Vehicle Routing Problem with Relief Materials in Sudden Disasters

2.1 EmVRP

In the Introduction above, we described briefly some of the prior work focused on EmVRP. In this section, we review examples of other significant contributions to the field.

Hale and Mober [23] studied the selection of emergency logistics supply nodes, with particular attention to the number of nodes (sites) for storage of emergency supplies. Their research proposed a quantitative model. For vehicle routing within a time window (VRPTW) involving a limited number of vehicles, Lau et al. [30] constructed a model and put forward the related solutions. In addition, they implemented a substantial number of fundamental studies on the variants of the traditional VRP problem.

Based on the features of multi-resource and multipoint rescue, Dai et al. [6] constructed a mathematical model to solve the multi-resource emergency problem. Fu et al. [15] divided emergency logistics distribution issues into three categories: road conditions, vehicle conditions, and the state of materials. They proposed a static, dynamic, and functional model in accordance with each emergency issue category. To address actual challenges of VRP with time windows, Li [32] put forward a method based on affinity group match, and further proved that the immune algorithm can provide an effective solution for routing vehicle within a time frame.

Özdamar et al. [38] studied distribution decision support systems after natural disasters. A model was established by considering disaster relief materials comprehensively, regarding vehicles as both relief materials and transportation resources. In addition, Ghiani [21], Ikou [24], and Özyurt [39] proposed several different algorithms to solve the EmVRP problem.

Though many scholars have carried out several in-depth studies regarding emergency relief, few of them has used swarm intelligence algorithm to solve emergency vehicle routing problem with relief materials in sudden disasters. In this paper, one of the most representative swarm intelligence algorithms, called monarch butterfly optimization (MBO), is introduced to solve emergency vehicle routing problem with relief materials in sudden disasters. Next, the mainframe of MBO algorithm will be provided.

2.2 MBO algorithm

Since MBO [55] was proposed, many scholars have explored this approach. In this section, some of the most representative work regarding MBO is summarized and reviewed.

Yi et al. [70] combined quantum computation theory into MBO algorithm, and proposed a novel quantum-inspired MBO methodology, called QMBO, by incorporating quantum computation into the basic MBO algorithm. In QMBO, a certain number of the worst butterflies were updated by quantum operators. The path planning navigation problem for Unmanned Combat Air Vehicles (UCAVs) was modeled into an optimization problem, and then its optimal path could be obtained by the proposed QMBO algorithm. Furthermore, B-Spline curves were utilized to refine the obtained path, making it more suitable for UCAVs. The UCAV path obtained by QMBO was studied and analyzed in comparison with the basic MBO. The experimental results showed that QMBO can find a much shorter path than MBO.

Ghetas et al. [20] incorporated the harmony search (HS) algorithm into the basic MBO algorithm, and proposed a variant of MBO, called MBHS, to deal with the standard benchmark problems. In MBHS, the HS algorithm was considered as a mutation operator to improve the butterfly adjusting operator, with the aim of accelerating the convergence rate of MBO.

Feng et al. [10] presented a novel binary MBO (BMBO) method used to address the 0-1 knapsack problem (0-1 KP). In BMBO, each butterfly individual was represented as a two-tuple string. Several individual allocation techniques were used to improve BMBO’s performance. In order to keep the number of infeasible solutions to a minimum, a novel repair operator was applied. The comparative study of BMBO with other optimization techniques showed the superiority of the former in solving the 0-1 KP.

Wang et al. [62] put forward another variant of the MBO method in combination with GCMBO. In GCMBO, two modification strategies, including a self-adaptive crossover (SAC) operator and a greedy strategy, were utilized to improve its search ability.

Feng et al. [14] combined chaos theory [50] with the basic MBO algorithm, and then proposed a novel chaotic MBO (CMBO) algorithm. The proposed CMBO algorithm enhanced the search effectiveness significantly. In CMBO, in order to tune two main operators, the best chaotic map was selected from 12 maps. Meanwhile, some of the worst individuals were improved by using a Gaussian mutation operator to avoid premature convergence.

Ghanem and Jantan [19] combined ABC with elements from MBO to proposed a new hybrid metaheuristic algorithm named Hybrid ABC/MBO (HAM). The combined method used an updated butterfly adjusting operator, considered to be a mutation operator, with the aim of sharing the information with the employee bees in ABC.

Wang et al. [59] proposed a discrete version of MBO (DMBO) that was applied successfully to tackle the Chinese TSP (CTSP). They also studied and analyzed the parameter butterfly adjusting rate (BAR). The chosen BAR was used to find the best solution for the CTSP.

Feng et al. [12] proposed a type of multi-strategy MBO (MMBO) technique for the discounted 0-1 knapsack problem (DKP). In MMBO, two modifications, including neighborhood mutation and Gaussian perturbation, were utilized to retain the diversity of the population. An array of experimental results showed that the neighborhood mutation and Gaussian perturbation were quite capable of providing significant improvement in the exploration and exploitation of the MMBO approach, respectively. Accordingly, two kinds of NMBO were proposed: NCMBO and GMMBO, respectively.

Feng et al. [13] combined MBO with seven kinds of DE mutation strategies, using the intrinsic mechanism of the search process of MBO and the character of the differential mutation operator. They presented a novel DEMBO based on MBO and an improved DE mutation strategy. In this work, the migration operator was replaced by a differential mutation operator with the aim of improving its global optimization ability. The overall performance of DEMBO was fully assessed using thirty typical discounted 0-1 knapsack problem instances. The experimental results demonstrated that DEMBO could enhance the search ability while not increasing the time complexity. Meanwhile, the approximation ratio of all the 0-1 KP instances obtained by DEMBO was close to 1.0.

Wang et al. [49] proposed a new population initialization strategy in order to improve MBO’s performance. Firstly, the whole search space is equally divided into N _P (population size) parts at each dimension. Subsequently, two random distributions (T and F distribution) are used to mutate the equally divided population. Accordingly, five variants of MBOs are proposed with a new initialization strategy.

Feng et al. [11] presented OMBO, a generalized opposition-based learning (OBL) [56] MBO with Gaussian perturbation. The authors used the OBL strategy on the portion of the individuals in the late stage of evolution, and used Gaussian perturbation on the individuals with poor fitness in each evolution. OBL guaranteed the higher convergence speed of OMBO, and Gaussian perturbation avoided the possibility of falling into a local optimum. For the sake of testing and verifying the effectiveness of OMBO, three categories of 15 large-scale 0-1 KP cases from 800 to 2,000 dimensions were used. The experimental results indicated that OMBO could find high-quality solutions.

Chen et al. [4] proposed a new variant of MBO by introducing a greedy strategy to solve dynamic vehicle routing problems (DVRPs). In contrast to the basic MBO algorithm, the proposed algorithm accepted only butterfly individuals that had better fitness than before implementation of the migration and butterfly adjusting operator. Also, a later perturbation procedure was introduced to make a trade-off between global and local search.

Meng et al. [34] proposed an improved MBO (IMBO) for the sake of enhancing the optimization ability of MBO. In IMBO, the authors divided the two subpopulations in a dynamic and random fashion at each generation, instead of using the fixed strategy applied in the original MBO approach. Also, the butterfly individuals were updated in two different ways for the sake of maintaining the diversity of the population.

Faris et al. [9] modified the position updating strategy used in the basic MBO algorithm by utilizing both the previous solutions and the butterfly individuals with the best fitness at the time. For the sake of fully exploring the search behavior of the Improved MBO (IMBO), it was benchmarked by 23 functions. Furthermore, the IMBO was applied to train neural networks. The IMBO-based trainer was verified on 15 machine learning datasets from the UCI repository. The experimental results showed that the IMBO algorithm could enhance the learning ability of neural networks significantly.

Ehteram et al. [7] used the MBO algorithm to address the utilization of a multi-reservoir system for the sake of improving production of hydroelectric energy. They studied three periods of dry (1963-64), wet (1951-52), and normal (1985-86) conditions in a 4-reservoir system. The experiments indicated that MBO can generate more energy when compared with particle swarm optimization (PSO) and a genetic algorithm (GA).

Xue et al. [67] added a self-adaptive strategy to the basic ABC algorithm in accordance with the global optimal solution, so a new variant of ABC named SABC-GB was proposed. SABC-GB has shown its superiority to other meta-heuristic algorithms when dealing the complicated optimization problems.

In addition to the MBO algorithm studied in this paper, many other intelligent algorithms [65] have been proposed, such as elephant herding optimization (EHO) [33, 57], simulated annealing (SA) [27], evolutionary strategy (ES) [2], particle swarm optimization (PSO) [26, 44], moth search (MS) algorithm [47], bat algorithm (BA) [35, 68], differential evolution (DE) [43, 53], biogeography-based optimization (BBO) [42, 52], krill herd (KH) [16, 51, 54], cuckoo search (CS) [5, 69], artificial bee colony (ABC) [25, 64], genetic algorithm (GA) [22], fireworks algorithm (FWA) [46], earthworm optimization algorithm (EWA) [61], and harmony search (HS) [18, 63]. These algorithms are used widely in various engineering applications [73, 74].

In this paper, the model of emergency vehicle routing problem with relief materials in sudden disasters is provided. Also, we will further improve the performance of the basic MBO algorithm by introducing the self-adaptive strategy and crossover operator to form a novel enhanced MBO (EMBO) algorithm. The proposed EMBO algorithm is then used to tackle EmVRP problem in comparison with the basic MBO algorithm and seven other intelligent algorithms.

3.1 Problem description

If the scale of a sudden natural disaster is large, the distribution of emergency materials may require a variety of transport modes as well as the conversion of various transportation modes. In our current work, we focused on one mode of transportation only, assuming that relief supplies brought by rail, air, and water had already arrived at train stations, airports, and docks in the general disaster area. Therefore, those terminals (depots, railway stations, airports, and docks) were also seen as material storage facilities. Our goal was to determine how to assign vehicles from parking areas (such as garages) to transport relief materials and goods from the various storage sites (depots, railway stations, airports, and docks) to the primary disaster points. An overview of the whole system involved in the emergency VRP with materials in sudden disasters is shown in Fig. 1 [75].

As shown in Fig. 1, the emergency relief material distribution system is a three-layer structure. Transporting the materials is essentially a vehicle routing problem involving movement from multiple parking areas (garages) to multiple storage sites, and then to multiple disaster points. Layer 1 includes the parking areas; all the vehicles assigned to carrying relief materials start from the parking lots. The number of vehicle types may vary, but the number is known.

Figure 1

The structure of three hierarchy of vehicle scheduling for emergency relief goods.

Layer 2 is made up of the material reserve storage sites, including all wharfs, train stations, airports, and similar facilities. We assume that the relief supplies at each reserve point are varied, and that the number of certain types of materials is unchanged during the emergency period. The bottom layer, Layer 3, comprises the affected areas that are the destinations for relief distribution. After a disaster, demand for all kinds of emergency supplies will be generated at all the impact points according to the disaster assessment. This demand will determine the need to move reserves at all levels and locations, which in turn will call for a certain number of vehicles to begin distribution of aid materials.

From Layer 1 to Layer 2, there is vehicle flow only, with no material flow, i.e., the vehicles run from Layer 1 to Layer 2, but they are not loaded. Moreover, this flow is one-way, with the vehicles moving only from the parking areas to the various material reserve sites. They do not return to the parking areas during the emergency period. However, both vehicle flow and material flow occur between Layer 2 and Layer 3. Materials move in a one-way unidirectional flow from Layer 2 to Layer 3, while the vehicles form a two-way flow between these two layers. After a vehicle arrives at a disaster point, it does not immediately return to the parking areas, but instead remains in the arrival spot on standby. Once a new distribution task is ordered, the vehicle will go back to the appropriate material storage site to be loaded and start a new distribution.

Therefore, after a vehicle arrives at an affected area from a reserve area, there are two possible states for the vehicle. Either it remains in the impacted area waiting for new orders, or it may be given a new task immediately, in which case it returns right way to the reserve to begin the next round of distribution, as shown by the dotted line in Fig. 1. In our current work, we assumed that the vehicle is running along the unloaded reverse flow. In this process, the flow of vehicles is significantly different from that of common commercial logistics.

Both materials and vehicles flow between Layers 2 and 3. However, there is no connection between each point of the same layer because the demand for supplies after a sudden disasters is much greater than the normal customer demand for small batches commonly encountered in business logistics. This difference also is reflected by the characteristics of the emergency relief supplies demanded by the affected points. The primary goal of post-disaster distribution is to meet the needs of all the disaster sites as much as possible while shortening the running time of the entire distribution system. In this way, losses are reduced at each disaster point, and the cost of vehicle operations is reduced as well.Vehicles can be reused throughout the distribution system, i.e., the vehicle flow is repeated between the second and the third layers until all the needs of the disaster sites are met.

3.2 Definition and explanation of boundary conditions of model

Under the guidance of the distribution system model of the EmVRP with relief materials in sudden disasters, the boundary conditions of the vehicle distribution system are defined as follows:

1. There are multiple supply points and demand points in the network. The supply and demand of emergency materials are known at each point, and the total supply can meet the needs of the affected areas.

2. There are many parking areas, and each parking area may have a variety of vehicles. The number of parked vehicles is sufficient to meet the need for materials distribution, and each vehicle has a number.

3. There are a variety of emergency supplies to be delivered. Each is different in weight and volume, and each has different loading efficiency. It is assumed that the loading efficiency of the material at the supply site is the same as the unloading efficiency at the disaster site.

4. Each disaster site can be served by multiple vehicles.

5. Each delivery task carries only one kind of material.

6. After the completion of a distribution task, if the vehicle still has supplies, then it will continue to complete the delivery task; otherwise, it will return to the starting point at the parking area to wait for the next dispatch.

7. There is neither vehicle flow nor material flow between the points of each layer.

3.3 The model of the EmVRP with relief materials in sudden disasters

The model of the emergency VRP with relief materials in sudden disasters consists mainly of four parts: the symbol definition, objective function, constraint condition, and the model description.

(1) Symbol definition

The model of the EmVRP with relief materials in sudden disasters includes the definition of five kinds of symbols: set, parameters related to transportation, parameters related to relief materials, parameters related to distance, and decision variables.

1) Set

The set of materials/goods: G={G1, G2,…, Gp};

The set of supply sites (i.e., relief materials reserve sites/points, material storage sites/ points): S={S1,  S2,…,Sn};

The set of disaster sites (affected areas/points): D={D1,D2,…,Dm};

The set of parking areas (garages, parking lots): K={K1,K2, …,Kk};

The set of vehicles: L={L1,L2,…,Ll};

The set of edges: E={(k,i)(i,j)|k∈K,  i∈S,j∈D}.

2) Parameters related to transportation

V_l: the maximum volume of the vehicle l;

cap_l: the maximum dead weight tonnage of the vehicle l;

v_l: the speed of the vehicle l.

3) Parameters related to relief materials

w_g: the unit weight of material type g;

c_g: the unit volume of material type g;

t_g: the time needed for loading and unloading material g.

4) Parameters related to distance

d _ki: the distance from the parking area to the material supply site;

d_ij: the distance between the supply site and the demand point.

5) Decision variables

x_lijg: the quantity of material g that is conveyed from the supply site i to j by vehicle l;

(2) Objective function

(3) Constraint condition

(4) Model description

For the model of the EmVRP with relief materials in sudden disasters, the goal is to meet the needs of the affected points in the shortest total running time possible. The total running time of a vehicle includes the travel time from the parking areas to material reserve sites, travel time from material reserve areas to the affected areas, loading time at reserve sites, and unloading time at the affected areas.

For the objective function in the model of the EmVRP with relief materials in sudden disasters, Eq. 1 represents the time from the parking area to the reserve site; Eq. 2 represents the loading time at the reserve site; Eq. 3 represents the time from the material reserve site to the affected site; and Eq. 4 represents the unloading time of the material at the affected area. In our current work, we suppose the vehicles will take the same time moving from the affected point to reserve point.

For the constraint conditions, Eq. 5 means that the goods transported by each vehicle cannot exceed the maximum load capacity of the vehicle; Eq. 6 means that the goods delivered by each vehicle should not exceed the maximum volume of the vehicle; Eq. 7 means that each vehicle carries only one material at a time from the material reserve site to the affected area; Eq. 8 means that the vehicles starting from the garage can only reach a material reserve site, which is a one-way non-circular flow between the first and second layers; and Eqs. 9-10 indicate that the variables of this problem satisfy the 0-1 integer constraint. Eq. 11 represents the range of the indexes in constraints. |·| represents the number of elements contained in a set.

To solve the problem of the model constructed in this paper, it is necessary to construct a heuristic algorithm with simple operation and excellent search performance. At present, the heuristic algorithms used to solve the vehicle scheduling problem include primarily the genetic algorithm (GA), neural network method, ant colony algorithm (ACO), tabu search (TS), and simulated annealing algorithm (SA). Among them, the most popular strategies involve the application of modern intelligent algorithms to vehicle scheduling. Many scholars have proven that the performance achieved using intelligent algorithms is better than the results from other traditional optimization algorithms. Therefore, this paper employed a new intelligent algorithm, monarch butterfly optimization (MBO), to solve the above-established EmVRP model. The required solutions can be represented by a butterfly individual, while a butterfly individual can be represented by different parking areas, reserve sites, and the affected areas.

3.4 Difference between emergency VRP and regular VRP

Emergency VRP is a special vehicle routing activity caused by outbreak of emergencies. Due to its suddenness and lack of information, it is quite different from that of regular VRP (RVRP) [29].

The goal of RVRP is to minimize costs or maximize profits [75]. The distribution network is permanent, and its network structure is designed by the needs of the customer [75]. Participants are mainly all kinds of economic entities closely related to each other, such as manufacturers, distributors and transport companies [75]. The driving force of logistics is demand, and facilities planning, fleet size and vehicle routing need to be planned in the long, medium and short term, respectively [75]. The quantity and variety of materials are limited, the means of transportation are used for a long time, and the information in the external environment is full and not easy to change [75].

The goal of emergency VRP is to meet the requirements of minimizing the delay time (primary goal) and minimizing the cost (secondary goal) [75]. The distribution network of logistics facilities is temporary, the network structure is simple and the functions are simplified [75]. The main participants are not closely linked to the interests of institutions or organizations, such as government departments, non-governmental organizations, donors and institutions that are temporarily established [75]. The driving force of logistics is divided into two kinds, the reflection stage is system propelling, and the recovery stage is demand promotion [75]. Usually, resource reserve planning is made [75]. In wartime or emergency, material transportation and vehicle scheduling plan are urgent, making the best possible decisions under shorter term and limited information, and the plan scheme will often change greatly in the implementation process [75]. There are the quantity and variety of material stock, the means of transportation are temporarily collected, the information in the external environment is inadequate and easy to change [75].

4.1 MBO algorithm

4.2 EMBO algorithm

MBO has been shown to have its own advantages over other intelligent algorithms for benchmarking and other application engineering problems [55]. However, as mentioned before, sometimes MBO may be stuck at local optima on certain problems [55]. In this paper, a self-adaptive strategy and crossover operators were combined with the basic MBO approach for the sake of enhancing the searchability of MBO. This enhanced MBO (EMBO) algorithm will be described in detail later in this paper.

Initialization. Set the generation counter t = 1, and set the maximum generation t_max, NP₁, NP₂, BAR, peri, and p;

Population evaluation. Calculate the fitness according to the objective function;

Algorithm 1

Monarch Butterfly Optimization

5.1 Individual coding and initial solutions construction

In our work, we utilize an improved natural number coding strategy [75]. A butterfly individual is a string representing an emergency relief material transportation scheme. In general, a butterfly individual is also called a chromosome in EAs.

A butterfly individual is made up of two substrings. The first substring has an element that represents the vehicle number. If there are K vehicles, the first element is an integer selected from 1 to K. The second substring has 3n elements, and n represents the number of tasks completed by the vehicle. For example, vehicle 2 at parking area 1 arrives at material storage site I₃ to transport material G₁ to the affected area J₂. Then vehicle 2 goes to material storage site I₁ to transport material G₂ to the affected area J₄. This process can be represented as K₁ − I₃ − G₁ − J₂ − I₁ − G₂ − J₄, and its corresponding butterfly individual can be expressed as 2-3-1-2-1-2-4. The (3n − 1)th element represents the number of the material reserve site, the (3n)th element represents the material type number, and the (3n + 1)th element represents the number of the disaster site (affected area). The element segments of all vehicles are arranged in parallel from small to large in order to form a single butterfly individual.

Suppose there are 3 vehicles that are numbered 1, 2, 3, located in the parking lot K₂, K₁, and K₃, respectively; there are 2 materials that are numbered G1 and G2 ; there are 2 reserve sites that are numbered I₁ and I₂ ; and there are 3 disaster sites that are numbered J₁, J₂, and J₃. The following butterfly individuals can be generated:

Butterfly individual 1:

Its corresponding solution is:

Butterfly individual 2:

Its corresponding solution is:

Butterfly individual 1 represents vehicle 1 from parking lot K₂, which arrives at material storage site I₂ to transport material G₁ to the affected area J₁. Then vehicle 1 goes to material storage site I₁ to transport material G₂ to the affected area J₂. Vehicle 2 from parking lot K₁ arrives at material storage site I₁ to transport material G₂ to the affected area J₃. Vehicle 3 from parking lot K₃ arrives at material storage site I₂ to transport material G₁ to the affected area J₂, and then vehicle 3 goes to material storage site I₁ to transport material G₂ to the affected area J₂. When the above transportation tasks have been completed, all requirements of all the affected areas will have been met. Butterfly individual 2 can be explained in the same way.

Calculate the butterfly adjusting rate BAR by using Eq. 18.

Algorithm 2

SABA operator

Initialization. Set the generation counter t = 1, and set the maximum generation t_max, NP₁, NP₂, BAR₀, peri, and p.

Population evaluation. Evaluate the butterfly individuals with their objective functions.

Algorithm 3

EMBO approach

Now, we will describe how to generate the initial solutions. First, vehicle l is selected randomly from the parking areas. This selection will be considered the starting point. Second, the reserve site S_i (i = 1, 2, ..., n) is selected randomly from reserve site set S. Finally, the affected area D_j (j = 1, 2, ..., m) is selected randomly from the affected area set D. This selection process is repeated until a whole butterfly individual is constructed, so the initial solution of the EmVRP is obtained. After all the butterfly individuals have been generated, the initial butterfly population is formed.

5.2 The EMBO for emergency VRP with relief materials in sudden disasters

In this section, we discuss how to use EMBO to solve the EmVRP with relief materials in sudden disasters. First, the initialization process is implemented, including the parameters and initial population, as described previously. Next, the optimization process is implemented to update butterfly individuals in the population. Then the newly generated butterfly individuals are evaluated according to the objective functions. The implementation of the EMBO algorithm is repeated until the given requirements are met. Finally, the optimal butterfly individual is decoded in order to get the final solution scheme for the EmVRP with relief materials in sudden disasters as given in Algorithm 4.

In Algorithm 4, Line 4 encoding is essentially the population initialization process needed so that the EMBO algorithm can solve the EmVRP. Through this encoding, the initial population is generated. Line 4 describes the decoding process that transforms the butterfly individual into the actual solution scheme for the EmVRP. Encoding and decoding are opposite processes.

6.1 An EmVRP with relief materials in sudden disasters (case 1)

To test our approach, we solved a hypothetical EmVRP with relief materials in a sudden-onset disaster [75] that can be described as follows. Suppose a sudden natural disaster occurred. There are 4 disaster spots in need of emergency materials. These locations can be identified as J₁, J₂, J₃, and J₄, respectively. The relief materials have been transported to the local airport and railway station by aircraft and railway. There are a total of 3 reserve sites, denoted by I₁, I₂, and I₃, respectively, including the airport and railway station together with the local relief supply reserve area. A total of 20 vehicles are collected, and they are located randomly in 3 parking locations numbered K₁, K₂, and K₃. There are 4 kinds of materials to be delivered: tents, quilts, clothing, and food, represented as G₁, G₂, G₃, and G₄, respectively. The number of tasks to be completed by each vehicle is not more than 5. Tables 1-3 provide the related information for the above EmVRP example.

6.2 The shortest time used by nine algorithms on the EmVRP with relief materials in sudden disasters

For the sake of carrying out a fair comparison, all the approaches were compiled using MATLAB R2017a (9.2) running under the Windows 10 Enterprise operating system on a PC with an Intel(R) Core(TM) i5-4590 CPU operating at 3.30 GHz, 8.00GB of RAM, and a hard drive of 1024 GB.

In all experiments for MBO and EMBO, we used the same parameter settings: probability p = 5/12, elitism number Keep = 2, BAR0=(5−1)/2,population size N_P = 50, dimension D = 320, and maximum generation t_max = 100.

In this paper, the proposed EMBO algorithm will be compared with seven intelligent algorithms when dealing with emergency vehicle routing problem with relief materials in sudden disasters. The seven comparative algorithms

Initialization. Set the generation counter t = 1, and set the maximum generation t_max, NP₁, NP₂, BAR₀, peri, and p.

Encoding. Initialize the population according to the encoding method in Section 5.1.

Population evaluation. Calculate the shortest time used by each butterfly individual.

Algorithm 4

EMBO algorithm for emergency VRP with relief materials in sudden disasters

can be described below. ABC (artificial bee colony) [25] is an intelligent optimization algorithm based on the smart behavior of honey bee swarm. BA (bat algorithm) [17] is a new powerful and efficient meta-heuristic optimization algorithm inspired by the echolocation behavior of bats with varying pulse rates of emission and loudness. BBO (biogeography-based optimization) [42] is a new evolution algorithm developed for the global optimization inspired by the immigration and emigration of species between islands (or habitats) in search of more compatible islands. CS (cuckoo search) [69] is a meta-heuristic optimization algorithm inspired by the obligate brood parasitism of some

cuckoo species by laying their eggs in the nests of other host birds (of other species). DE (differential evolution) [43] is a simple but excellent optimization method that uses the difference between two solutions to probabilistically adapt a third solution. An ES (evolutionary strategy) [2] is an algorithm that generally distributes equal importance to mutation and recombination, and that allows two or more parents to reproduce an offspring. PSO (particle swarm optimization) [26, 71, 72] is also a swarm intelligence algorithm which is based on the swarm behavior of fish, and bird schooling in nature.

The parameters for the other seven algorithms were set as follows.

1. For ABC, the population size N _P = 50, the number of food sources FoodNumber = N _P/2, maximum search times limit = 100.

2. For BA, loudness A = 0.5, pulse rate r = 0.5, and scaling factor ε = 0.001.

3. For BBO, habitat modification probability = 1, immigration probability bounds per gene = [0, 1], step size for numerical integration of probabilities =1, maximum migration rate for each island = 1, and mutation probability = 0.005.

4. For CS, the discovery rate p_a = 0.25.

5. For DE, the weighting factor F = 0.5, and crossover constant CR = 0.5.

6. For ES, the number of offspring produced in each generation λ = 10, and the standard deviation for changing solutions is σ = 1 .

7. For PSO, the inertial constant = 0.3, the cognitive constant = 1, and the social constant for swarm interaction = 1.

For one implementation of EMBO, the convergent trend of the best and average time of the EMBO approach is shown in Fig. 2. In Fig. 2, the best and mean time after 100 generations were 109.4341 and 109.489, respectively. Furthermore, the EMBO algorithm converged to the least time after about 20 generations. These findings indicated that the EMBO algorithm could solve the emergency VRP problem well.

Figure 2

The convergent trend of the best and average time for EMBO algorithm.

In addition, the shortest time over 30 independent runs was obtained by EMBO and the eight other intelligent algorithms, and the results are displayed in Table 6. From Table 6, although EMBO had the biggest std value, it is clear that EMBO was able to complete the task with the least time among the nine intelligent algorithms for the average, best, and worst performance. For the other eight intelligent algorithms, BBO, CS, and DE had similar performances that were inferior only to the EMBO algorithm. Also, MBO performed similarly to PSO, which was better than ABC, BA, and ES.

The convergent trend of the average time obtained by the nine algorithms is shown in Fig. 3. From Fig. 3, it is clear that EMBO converged in the fastest fashion among the nine intelligent algorithms. This convergent trend was consistent with the results provided in Table 6.

Figure 3

Results obtained by nine algorithms.

6.3 Another EmVRP with relief materials in sudden disasters (case 2)

Like case 1 studied in Section 6.1, another emergency VRP problem is further used to verify our proposed EMBO algorithm. In case 2, there are 5 disaster spots in need of emergency materials. These locations can be identified as J₁, J₂, J₃, J₄, and J₅, respectively. The relief materials have been transported to the local airport and railway station by aircraft and railway. There are a total of 4 reserve sites, denoted by I₁, I₂, I₃, and I₄, respectively. A total of 25 vehicles are collected, and they are located randomly in 3 parking locations numbered K₁, K₂, and K₃. There are 4 kinds of materials to be delivered: tents, quilts, clothing, and food, represented as G₁, G₂, G₃, and G₄, respectively. The number of tasks to be completed by each vehicle is not more than 5. The parameters of emergency supplies used in case 2 is the same with case 1, as shown in Table 2. Tables 7-10 provide the other related information for case 2.

For the software, hardware environments and parameter settings used in nine intelligent algorithms, they are the same with case 1.

In addition, the shortest time over 30 independent runs was obtained by EMBO and the eight other intelligent algorithms, and the results are displayed in Table 11. From Table 11, it is clear that EMBO was able to complete the task with the least time among the nine intelligent algorithms for the average, best, and worst performance. Also, EMBO had the smallest Std value for this case. For the other eight intelligent algorithms, BBO, BA, and MBO had similar performances that were inferior only to the EMBO algorithm. Comparing with case 1, EMBO algorithm takes less time to complete the task, and the reason is that five more vehicles are added to case 2.

6.4 Parameter study

As we are aware, for all the intelligent algorithms, the parameter setting is of significant importance to their performance. In this section, the effectiveness of the initial butterfly adjusting rate BAR₀, probability p, elitism Keep, population size N_P, and maximum generation t_max are analyzed for the proposed EMBO algorithm on case 1.

6.4.1 Influence of the initial butterfly adjusting rate BAR₀ for EMBO

First, we examined the initial butterfly adjusting rate BAR₀. The other parameters were set as follows: elitism number Keep = 2, probability p = 5/12, population size N_P = 50, and maximum generation t_max = 100. The final times required by EMBO with different values of BAR₀ were recorded as shown in Table 12 and Fig. 4. From Table 12, when BAR₀ was equal to 0.6 or 0.7, EMBO could find the optimal scheduling scheme with the least time. Looking carefully at Fig. 4, we can observe that the performance of EMBO was significantly improved with the increment of BAR₀. However, EMBO had a similar performance when BAR₀ was between 0.6 and 1.0. Considering matters overall, we set the initial butterfly adjusting rate BAR0  to(5−1)/2for this present work.

Figure 4

Results obtained by EMBO with different BAR₀.

Table 12

The least time used by the proposed EMBO algorithm for the EmVRP with relief materials in sudden disasters with different BAR₀.

	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
Best	135.36	135.39	134.07	135.52	109.44	109.34	108.58	109.02	108.85	110.42	110.46
Mean	138.17	138.09	138.10	137.63	133.57	114.28	109.51	109.51	109.47	110.91	110.91
Worst	140.69	140.84	140.55	139.91	140.86	137.88	110.14	109.94	109.94	111.03	111.03
Std	1.5266	1.3504	1.3326	0.8532	9.7697	10.3059	0.3526	0.2779	0.3323	0.1675	0.1856

It should be mentioned that the initial butterfly adjusting rate BAR₀ was between 0 and 1.0. When the initial butterfly adjusting rate BAR₀ = 0, the updating Eq. 18 would degenerate into the following equation:

(26)BAR=ttmax.

In this case, the butterfly adjusting rate BAR is changing in a linear fashion. On the other hand, when the initial butterfly adjusting rate BAR₀ = 1.0, the butterfly adjusting rate BAR is equal to the initial butterfly adjusting rate BAR₀. The butterfly adjusting rate BAR will not change during the whole optimization process. Essentially, this description fits a basic MBO algorithm.

6.4.2 Influence of the probability p for EMBO

In EMBO, there is another probability p that decides the number of NP₁ and NP₂. In other words, there are NP₁ and NP₂ butterflies that implement the migration operator and butterfly adjusting operator, respectively. Here, we studied the probability p. The other parameters were set as follows: population size N_P = 50, the initial butterfly adjusting rate BAR0=(5−1)/2,elitism number Keep = 2, and maximum generation t_max = 100. The final times achieved by EMBO with different p are shown in Table 13 and Fig. 5. From Table 13, when probability p was equal to 0.7, 0.8 or 0.7, we see that EMBO could find the optimal scheduling scheme with the least time. Examining Fig. 5 carefully, we can observe that the performance of EMBO improved significantly with the increment of p. However, EMBO had a similar performance when p was between 0.5 and 1.0. Considered overall, the probability p was set to 5/12 in this paper.

Figure 5

Results obtained by EMBO with different p.

Table 13

The least time used by the proposed EMBO algorithm for the emergency VRP with relief materials in sudden disasters with different p.

	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
Best	136.05	109.57	109.43	109.16	108.89	108.92	108.62	108.81	109.33
Mean	137.97	136.86	127.41	117.89	111.30	109.73	109.58	108.58	110.08
Worst	140.56	141.69	139.68	139.36	138.08	116.14	110.47	110.54	110.67
Std	1.0755	7.4311	13.6690	12.6734	6.4444	1.2843	0.3514	0.4003	0.3664

It should be mentioned that the probability p was between 0.1 and 0.9. For MBO, when probability p = 0, all the individuals will be updated only by the butterfly adjusting operator. When probability p = 1.0, all the individuals will be updated only by the migration operator.

No special attention was paid to these two extreme cases of p = 0 and p = 1.0.

6.4.3 Influence of the elitism Keep for EMBO

As we can see, most intelligent algorithms involve an elitism strategy that has a great influence on the performance. Here, the elitism Keep was studied in the range of [0, 10]. The other parameters were set as follows: probability p = 5/12, initial butterfly adjusting rate BAR0=(5−1)/2,population size N_P = 50, and maximum generation t_max = 100. The final times attained by EMBO with different Keep values are shown in Table 14 and Fig. 6. From Table 14, when the parameter Keep was equal to 0, 3 or 8, EMBO could find the optimal scheduling scheme with the least time. Looking carefully at Fig. 6, though the trend of EMBO is less obvious, generally speaking the performance of EMBO improved with the increment of Keep. However, on average, EMBO had a similar performance when Keep was equal to 1 and 3. Therefore, we set the elitism Keep to 2 for this research.

Figure 6

Results obtained by EMBO with different Keep.

Table 14

The least time used by the proposed EMBO algorithm for the emergency VRP with relief materials in sudden disasters with different Keep.

	0	1	2	3	4	5	6	7	8	9	10
Best	108.89	109.19	109.19	109.16	109.23	109.07	109.06	108.92	109.21	109.15	109.43
Mean	114.14	112.52	115.36	112.39	112.40	126.95	120.49	121.27	123.02	118.78	121.51
Worst	137.74	137.92	139.75	137.83	138.34	138.35	138.55	138.79	137.49	137.94	139.17
Std	10.2456	8.4974	11.5944	8.4129	13.8968	12.2843	13.5096	13.5515	13.6949	13.0870	13.5949

6.4.4 Influence of the population size NP for EMBO

Subsequently, we studied the population size N_P in the range of [10, 100] with interval 10. The other parameters were set as follows: the initial butterfly adjusting rate BAR0=(5−1)/2,probability p = 5/12, elitism number Keep = 2, and maximum generation t_max = 100. The final times obtained by EMBO with different N_P were recorded as shown in Table 15 and Fig. 7. From Table 15, when population size N_P was equal to 90, 100, 20, and 10, EMBO could find the optimal scheduling scheme with the least time. Examining Fig. 7, we can observe that the performance of EMBO decreased and then improved with the increment of N_P. Furthermore, overall the proposed EMBO algorithm was capable of obtaining satisfactory solutions within a reasonable time when the population size N_P = 10. Under this condition, fewer computational resources were required, indicating that EMBO could solve the EmVRP while consuming only limited computational resources. Considered altogether, we set the population size N_P to 50 in our present work.

Figure 7

Results obtained by EMBO with different NP.

Table 15

The least time used by the proposed EMBO algorithm for the emergency VRP with relief materials in sudden disasters with different NP.

	10	20	30	40	50	60	70	80	90	100
Best	110.25	109.63	109.16	109.19	109.20	109.16	109.16	108.89	108.85	108.85
Mean	111.38	111.37	111.95	112.60	114.34	118.55	119.17	115.55	115.27	117.76
Worst	116.17	139.58	139.23	137.49	139.86	138.85	137.51	136.82	135.96	135.97
Std	1.8294	5.4166	7.2907	8.3299	10.4997	12.9187	13.0635	11.3207	10.9097	12.2106

6.4.5 Influence of the maximum generation t_max for EMBO

Last, we studied the maximum generation t_max in the range of [20, 200] with interval 20. The other parameters were set as follows: the initial butterfly adjusting rate BAR0=(5−1)/2,elitism number Keep = 2, probability p = 5/12, and population size N_P = 50. The final times used by EMBO with different t_max are shown in Table 16 and Fig. 8. From Table 16, when the maximum generation t_max was equal to 80, 60, and 60, EMBO was able to get the optimal scheduling scheme with the least time. Looking carefully at Fig. 8, generally the performance of EMBO was significantly reduced with the increment of t_max. However, when t_max was equal to 20, 40, and 60, EMBO provided better performance than at the other t_max values. This finding goes against common expectations. For most intelligent algorithms, their performance will exceed, or at least be equal to, the performance given before the increment of generations. We propose the reason for this finding is that maximum generation t_max has an influence on the value of the butterfly adjusting rate BAR, which, in turn, has a big impact on MBO’s performance. Under this condition, the situation becomes complicated. Considered as a whole, weset the maximum generation t_max to 100 in our present work.

Figure 8

Results obtained by EMBO with different t_max.

Table 16

The least time used by the proposed EMBO algorithm for the emergency VRP with relief materials in sudden disasters with different t_max.

	20	40	60	80	100	120	140	160	180	200
Best	110.09	109.65	109.39	108.59	109.06	109.20	108.89	109.12	109.12	108.85
Mean	110.68	110.32	109.85	114.43	112.41	116.01	116.56	117.99	122.81	120.91
Worst	111.03	114.64	110.30	141.23	139.51	137.50	138.50	136.32	135.83	136.52
Std	0.2450	0.8589	0.2236	10.8763	8.7851	11.6742	11.7442	12.2216	12.6813	12.3809