Solving the humanitarian multi-trip cumulative capacitated routing problem via a grouping metaheuristic algorithm

Abstract

Every year, natural disasters such as earthquakes, floods, volcanos, etc. cause millions of victims. So a quick response to these disasters is vital to reduce their negative consequences. Vehicle routing models can make important contributions to faster response, and thus, save lives. This paper proposes a vehicle routing problem to deliver relief resources from origins to destinations in response to disasters. For this purpose, a multi-period, multi-depot, multi-trip mixed-integer linear programming model is developed. Minimizing the sum of arrival times is considered as a service-based objective function for the increase of the survival rate. For the first time, the problem is solved using a grouping metaheuristic algorithm. Then its performance is compared with two other grouping algorithms. To evaluate the solution method, the algorithms are implemented on various test problems and compared statistically. Additionally, to show the validity of the model, sensitivity analyses are performed and managerial insights are given.

Appendices

Appendix A

Notations are defined below:

Notation	Definition
Index sets
$N$	Set of all nodes indexed by $i,j\in N$
$D$	Set of demand nodes indexed by $d, {d}^{\prime}$
$S$	Set of distribution centers indexed by $s, {s}^{\prime}$
$V$	Set of vehicles indexed by $v$
$C$	Set of commodity types indexed by $c$
$T$	Set of periods indexed by $t$
$R$	Set of trips any vehicle makes indexed by $r$
Parameters
${d}_{cdt}$	Amount of demand for commodity $c$ in node $d$ in period $t$
${t}_{ijt}$	Travel time from node $i$ to node $j$ in period $t$
$VC$	Volume capacity of a vehicle
$WC$	Weight capacity of a vehicle
${C}_{c}$	Unit volume of commodity type $c$
${W}_{c}$	Unit weight of commodity type $c$
${\delta }_{ijts}$	Zero and one matrix that indicates accessibility to arc $(i,j)$ in period $t$
${t}_{max}$	Maximum working hours per vehicle per period
$M$	A large positive number

Notation	Definition	Description
Decision variables
${Z}_{vrijt}$	1 if vehicle $v$ in trip $r$ in period $t$ traverses arc $(i,j)$; otherwise, 0	This variable helps to determine the sequence of points, which each vehicle meets on each trip and gives the vehicle routing plan. Besides, it can help to determine the goods’ distribution plan
${st}_{rvt}$	Start time of trip $r$ for vehicle $v$ in period $t$	The departure time of a vehicle from its DC to make a specific trip in a specific time period is determined by this variable
${et}_{rvt}$	End time of trip $r$ for vehicle $v$ in period $t$	The return time of a vehicle to its DC after making a specific trip in a specific time period is determined by this variable
${at}_{ivrt}$	Arrival time of vehicle $v$ in trip $r$ in period $t$ in node $i$	This variable determines when a vehicle meets a node in the network on a given trip and period. With the aid of this variable, we can attain the operations’ scheduling plan

Appendix B

The steps of the F3EA algorithm are explained in the following.

2.1 The find step

This step is referred to as the selection step that imitates the process of military objects detection by a radar. To generate a new solution in each iteration, it is first assumed that the parent’s solution has a temporary role in a radar. On the other side, all other solutions in the population play the role of the enemy’s military facility and it is possible to detect it by the artificial radar. Identification of individuals is determined by the maximum radar range (${Rn}_{max}$).

Detection capability is measured based on the value of the function. The radar equation can be expressed as follows according to the function values:

$${Rn}_{max}=(f\left({\overrightarrow{F}}_{i}^{t}\right)-{f}_{min}^{t})\times \sqrt[4]{\left((\frac{{f}_{max}^{t}-f\left({\overrightarrow{F}}_{j}^{t}\right)}{{f}_{max}^{t}-{f}_{min}^{t}}){e}^{-(\frac{f\left({\overrightarrow{F}}_{j}^{t}\right)-{f}_{min}^{t}}{f\left({\overrightarrow{F}}_{i}^{t}\right)-{f}_{min}^{t}})}\right)}$$

(29)

Let’s assume that the worst and best function values in the population are shown by ${f}_{min}^{t}$ and ${f}_{max}^{t}$, respectively. The function values of the ith individual in the population (the parent solution) and the jth individual are represented by $f\left({\overrightarrow{F}}_{i}^{t}\right)$ and $f\left({\overrightarrow{F}}_{j}^{t}\right)$, respectively.

Definition 1

Let ${R}_{ij}^{t}=\left|f\left({\overrightarrow{F}}_{i}^{t}\right)-f\left({\overrightarrow{F}}_{j}^{t}\right)\right|$ be the distance between ${\overrightarrow{F}}_{i}^{t}$ and ${\overrightarrow{F}}_{j}^{t}$ in the function value space. ${\overrightarrow{F}}_{i}^{t}$ detect ${\overrightarrow{F}}_{j}^{t}$ if ${R}_{ij}^{t}\le {Rn}_{max ij}^{t}$.

Observation 1

If $f\left({\overrightarrow{F}}_{i}^{t}\right)={f}_{min}^{t}$, then ${R}_{max ij}^{t}=0, \forall i\ne j$. In a radar with the best solution, nothing is detected.

Observation 2

If $f\left({\overrightarrow{F}}_{j}^{t}\right)={f}_{min}^{t}$, then ${R}_{max ij}^{t}=f\left({\overrightarrow{F}}_{i}^{t}\right)-{f}_{min}^{t}={R}_{ij}^{t}$. All other individuals in the population can always detect the best solution in the population.

Observation 3

If $f\left({\overrightarrow{F}}_{j}^{t}\right)={f}_{max}^{t}$, then ${R}_{max ij}^{t}=0$. No individual can detect the worst solution.

Among all the detected solutions by the radar, one is selected based on any logic to generate the new solution in the Finish step. If no individual is identified, an individual of the population will be selected. Hereafter, we assume that ${\overrightarrow{F}}_{j}^{t}$ is the selected individual.

2.2 The fix step

In this step, which is a type of local search step, the imitation of the aiming process toward the target in order to monitor it is performed. As can be seen in Fig. 

Fig. 9

A schematic concept of the Fix step

9, when the military tank target has to be hit, the angle by which the rocket is launched has to be acute enough to let the rocket traverses the maximum peak. A single-variable optimization problem is used to obtain the improved solutions locally.

Let ${G}_{0}^{t}$ be defined as an individual that reveals the position of an artificial rocket launcher in the search space and ${S}^{t}$ as a search direction. The Fix step probes for the location of the deepest valley (the minimum of the scalar function) on the line passing through ${G}_{0}^{t}$, and in parallel with ${S}^{t}$, the following equations are used for one-dimensional optimization problem as follows:

$$min f({G}_{0}^{t}+\lambda {S}^{t})$$

$$s.t. {\lambda }_{min}\le \lambda \le {\lambda }_{max},$$

(30)

$${\lambda }_{min}=\underset{d=1,..,n}{\mathrm{max}}\left\{\left(\frac{\left({x}_{min,d}-{g}_{0d}^{t}\right)}{{S}_{d}^{t}}|{S}_{d}^{t}>0\right),\left(\frac{\left({x}_{max,d}-{g}_{0d}^{t}\right)}{{S}_{d}^{t}}|{S}_{d}^{t}<0\right)\right\},$$

(31)

$${\lambda }_{max}=\underset{d=1,..,n}{\mathrm{min}}\left\{\left(\frac{\left({x}_{max,d}-{g}_{0d}^{t}\right)}{{S}_{d}^{t}}|{S}_{d}^{t}>0\right),\left(\frac{\left({x}_{min,d}-{g}_{0d}^{t}\right)}{{S}_{d}^{t}}|{S}_{d}^{t}<0\right)\right\}.$$

(32)

where ${G}_{d}^{t}={G}_{d-1}^{t}+{S}_{d}^{t}{e}_{d}$ and components ${S}_{d}^{t}$, $\forall d=1,\dots ,d$ are obtained consecutively as follows:

$$min f({G}_{0}^{t}+{S}_{d}^{t}{e}_{d})$$

$$s.t. {\gamma }_{d}^{t}({x}_{min,d}-{g}_{d-1,d}^{t})\le {S}_{d}^{t}\le {\gamma }_{d}^{t}({x}_{max,d}-{g}_{d-1,d}^{t}),$$

(33)

where ${e}_{d}$ is the dth coordinate unit vector and ${\gamma }_{d}^{t}$ is every time set randomly selected as 0.05 or 1 to shrink the line search extent for local or global exploration, respectively. According to the ${\lambda }^{*}$ value obtained by Eqs. (30–32), the position of the newly generated solution is ${G}_{0}^{t}+{\lambda }^{*}{S}^{t}$.

2.3 The finish step

Let ${\overrightarrow{F}}_{i}^{t}$ and ${\overrightarrow{F}}_{j}^{t}$ represent the position of a rocket launcher and a target military facility, respectively. The main purpose of this step is to find a new solution based on ${\overrightarrow{F}}_{i}^{t}$ and benefiting from ${\overrightarrow{F}}_{j}^{t}$. It is assumed that an artificial rocket is launched from ${\overrightarrow{F}}_{i}^{t}$ toward ${\overrightarrow{F}}_{j}^{t}$ and the position of the explosion ${\overrightarrow{E}}_{ij}^{t}$ is determined via the projectile motion equations and considered as a new solution in the search space.

Let ${\left[{\overrightarrow{F}}_{i}^{t},u({\overrightarrow{F}}_{i}^{t})\right]}_{1\times (n+1)}$ and ${\left[{\overrightarrow{F}}_{j}^{t},u({\overrightarrow{F}}_{j}^{t})\right]}_{1\times (n+1)}$ be the positions of the rocket launcher and the jth artificial military object selected for devastation in the $F-u(F)$ coordinate system, respectively (see Fig. 

Fig. 10

The rocket’s move path is parabolic regardless of the air or wind drag

10).

As shown in Fig. 10, the inclination ${\beta }_{ij}^{t}$ with the horizon is calculated as follows:

$${\beta }_{ij}^{t}=\mathrm{arctan}(\frac{(u\left({\overrightarrow{F}}_{j}^{t}\right)-u({\overrightarrow{F}}_{i}^{t}))}{\Vert {\overrightarrow{F}}_{j}^{t}-{\overrightarrow{F}}_{i}^{t}\Vert })$$

(34)

Let the artificial rocket be launched from ${\overrightarrow{F}}_{i}^{t}$ toward ${\overrightarrow{F}}_{j}^{t}$ with angle ${\alpha }_{ij}^{t}$. The value of the artificial rocket’s initial velocity (${v}_{0ij}^{t}$) can be obtained for a given value of ${\alpha }_{ij}^{t}$ as bellow:

$${v}_{0ij}^{t}=\left\{\begin{array}{c}\sqrt{\frac{({AB}_{ij}^{t}\times g\times {cos}^{2}{\beta }_{ij}^{t})}{(2cos{\alpha }_{ij}^{t}\times \mathrm{sin}({\alpha }_{ij}^{t}-{\beta }_{ij}^{t}))}}, f\left({\overrightarrow{F}}_{i}^{t}\right)<f\left({\overrightarrow{F}}_{j}^{t}\right)\\ \sqrt{\frac{({AB}_{ij}^{t}\times g\times {cos}^{2}\left|{\beta }_{ij}^{t}\right|)}{\left(2cos{\alpha }_{ij}^{t}\times \mathrm{sin}\left({\alpha }_{ij}^{t}+\left|{\beta }_{ij}^{t}\right|\right)\right)}}, f\left({\overrightarrow{F}}_{i}^{t}\right)\ge f\left({\overrightarrow{F}}_{j}^{t}\right)\end{array}\right.$$

(35)

where,

$${AB}_{ij}^{t}=\Vert {\left[{\overrightarrow{F}}_{i}^{t},u({\overrightarrow{F}}_{i}^{t})\right]}_{1\times (n+1)}-{\left[{\overrightarrow{F}}_{j}^{t},u({\overrightarrow{F}}_{j}^{t})\right]}_{1\times (n+1)}\Vert ,$$

(36)

$$ \alpha _{{ij}}^{t} = \left\{ {\begin{array}{*{20}c} {rand\left( {\beta _{{ij}}^{t} ,\pi /2} \right),~~~~f\left( {\vec{F}_{i}^{t} } \right) < f\left( {\vec{F}_{j}^{t} } \right)} \\ {rand\left( {0,\pi /2} \right),~~f\left( {\vec{F}_{i}^{t} } \right) \ge f\left( {\vec{F}_{j}^{t} } \right).} \\ \end{array} } \right. $$

(37)

By Eq. (38), the time of the rocket’s flight is determined:

$${T}_{ij}^{t}=\left\{\begin{array}{c}\frac{(2{v}_{0ij}^{t} \mathrm{sin}({\alpha }_{ij}^{t}-{\beta }_{ij}^{t}))}{g cos{\beta }_{ij}^{t}}, f\left({\overrightarrow{F}}_{i}^{t}\right)<f\left({\overrightarrow{F}}_{j}^{t}\right)\\ \frac{(2{v}_{0ij}^{t} \mathrm{sin}({\alpha }_{ij}^{t}+\left|{\beta }_{ij}^{t}\right|))}{g cos\left|{\beta }_{ij}^{t}\right|}, f\left({\overrightarrow{F}}_{i}^{t}\right)\ge f\left({\overrightarrow{F}}_{j}^{t}\right).\end{array}\right.$$

(38)

After launching the rocket, it is exposed to two different conditions of an unanticipated air drag force and the wind force, which cause the deviation of explosion position from the predicted position (see Fig. 

Fig. 11

Air drag and wind force cause the rocket’s move path to be non-parabolic

11). The amount of this diversion depends on the mass of the rocket (${rm}_{ij}^{t}$).

It seems logical to assume that the wind force on the X-axis blows in an unswerving direction toward ${\overrightarrow{F}}_{i}^{t}$ to push the artificial rocket toward ${\overrightarrow{F}}_{j}^{t}$, if $f\left({\overrightarrow{F}}_{j}^{t}\right)<f\left({\overrightarrow{F}}_{i}^{t}\right)$. Otherwise, if $f\left({\overrightarrow{F}}_{i}^{t}\right)<f\left({\overrightarrow{F}}_{j}^{t}\right)$, it is logical to say that the wind force on the X-axis blows in a direction away from ${\overrightarrow{F}}_{i}^{t}$, and pushes the artificial rocket away from ${\overrightarrow{F}}_{j}^{t}$. No wind component in the z direction is assumed (${w}_{zij}^{t}=0$):

$${w}_{xij}^{t}=\left\{\begin{array}{c}-r{m}_{ij}^{t}\times {v}_{0ij}^{t}\times rand(\mathrm{0,1}), f\left({\overrightarrow{F}}_{i}^{t}\right)<f\left({\overrightarrow{F}}_{j}^{t}\right)\\ {rm}_{ij}^{t}\times {v}_{0ij}^{t}\times rand(\mathrm{0,1}), f\left({\overrightarrow{F}}_{i}^{t}\right)\ge f\left({\overrightarrow{F}}_{j}^{t}\right).\end{array}\right.$$

(39)

If ${E}_{i}^{t}$ is considered as the explosion position of the artificial rocket launched from ${\overrightarrow{F}}_{i}^{t}$ toward ${\overrightarrow{F}}_{j}^{t}$, which is a new solution to the problem, Eqs. 40–41 are used to generate ${E}_{i}^{t}$ as follows (Eq. 40 is the explosion position on the x-axis, and Eq. 41 represents the new solution generated in the Finish step):

$$ x\left( {T_{{ij}}^{t} } \right) = w_{{xij}}^{t} \times T_{{ij}}^{t} + rm_{{ij}}^{t} \left( {v_{{0ij}}^{t} ~cos\alpha _{{ij}}^{t} - w_{{xij}}^{t} } \right)\left( {1 - e^{{ - {{T_{{ij}}^{t} } \mathord{\left/ {\vphantom {{T_{{ij}}^{t} } {rm_{{ij}}^{t} }}} \right. \kern-\nulldelimiterspace} {rm_{{ij}}^{t} }}}} } \right) $$

(40)

$${E}_{i}^{t}={\overrightarrow{F}}_{i}^{t}+x\left({T}_{ij}^{t}\right)\frac{{(\overrightarrow{F}}_{j}^{t}-{\overrightarrow{F}}_{i}^{t})}{\Vert {\overrightarrow{F}}_{j}^{t}-{\overrightarrow{F}}_{i}^{t}\Vert }.$$

(41)

2.4 The exploit step

It is necessary that the outcome of the Finish step via the crossover operations be exploited in the Exploit step of the F3EA metaheuristic algorithm. The feasible solution ${E}_{i}^{t}$ differs from ${\overrightarrow{F}}_{i}^{t}$ in all dimensions. By contrast, for many functions, making changes in all directions may not be a good idea. To simulate the number of changes, a truncated geometric distribution is used.

2.5 The analyze step

There are two circumstances in the Analyze step; either when a new solution is generated during the Fix or Exploit step, which is able to provide a better function value, or when it starts to update the global best solution.

Fig. 12

An example of a new solution generation

Fig. 13

Average RPD value of the algorithms for medium-sized problems

Fig. 14

CPU time of the algorithms for medium-sized problems

Fig. 15

Average RPD value of the algorithms for large-sized problems

Fig. 16

Standard deviation of the algorithms’ RPD for large-sized problems

Fig. 17

CPU time of the algorithms for large-sized problems

Table 5 Comparison outputs of metaheuristics for medium-sized problems

Table 6 Non-parametric Friedman’s test for performance measures in medium-sized problems

Table 7 Detailed statistics of Wilcoxon’s signed-rank test for performance measures in medium-sized problems

Table 8 Comparison outputs of metaheuristics for large-sized problems

Table 9 Non-parametric Friedman’s test for performance measures in large-sized problems

Table 10 Detailed statistics of Wilcoxon’s signed-rank test for performance measures in large-sized problems