Multi-depot open vehicle routing problem with fuzzy time windows

Abstract

This paper presents one method and one hybrid genetic algorithm for multi-depot open vehicle routing problem with fuzzy time windows (MDOVRPFTW) without maximum time windows. For the method, the degree of customers’ willingness to accept goods (DCWAG) is firstly proposed, it’s one fuzzy vague and determines maximum time windows. Referring to methods to determine fuzzy membership function, the function between DCWAG and the starting service time is constructed. By setting an threshold for DCWAG, the starting service time that the threshold corresponds can be treated as the maximum time window, which meets the actual situation. The goal of the model is to minimize the total cost. For the algorithm, MDOVRPFTW without maximum time windows is an extension of the NP-hard problem, the hybrid genetic algorithm was designed, which is combination of genetic algorithm and Hungarian algorithm. When the hybrid genetic algorithm applied to one pharmaceutical logistics company in Beijing City, China, one optimal scheme is determined. Then the rationality and the stability of solutions by the hybrid genetic algorithm are proved. Finally, sensitivity analyses are performed to investigate the impact of someone factor on DCWAG and some suggestions are proposed.

Keywords

MDOVRPFTW maximum time window fuzzy membership function the degree of customers’ willingness to accept goods the hybrid genetic algorithm

1 Introduction

The vehicle routing problem (VRP) plays a central role in optimization of distribution networks [1]. Today, there are more diversified requirements of customers for distribution [2]. Technology development leads to improvement about distribution tools’ efficiency [3]. Merger and reorganization of logistics enterprises lead to expansion of distribution scale. The government formulates many requirements in order to promote the logistic industry developing stably [4]. To adapt to development of economic society, many distribution enterprises try to apply some frontier theories to vehicle routing problem, which strengthens connection between theory and industry and broadens researches of VRP. MDOVRPFTW as one extended VRP, will play an important role in the trend of warehousing integration and artificial intelligence.

MDOVRPFTW is combination of multi-depot open vehicle routing problem (MDOVRP) and vehicle routing problem with fuzzy time windows (VRPFTW). MDOVRP means that there are many depots in the distribution system, vehicles can travel from depots to customers and don’t need to return to any depot after completing tasks. There are two traditional solution methods for MDOVRP. One is to divide the distribution area into several independent areas and dispatch vehicles separately, the method simplifies MDOVRP and may not lead to the optimal allocation of global resources. Another is to schedule globally, which can lead to the optimal allocation of global resources, but it’s more complex. About researches of MDOVRP, Wang et al. [5] considered the cost of carbon emissions and designed a hybrid heuristic algorithm. Lim et al. [6] designed a two-stage algorithm. Salhi et al. [7] considered multi-type vehicles. Yu et al. [8] improved ant colony algorithm. Rahimi et al. [9] designed a module heuristic algorithm. Nilay et al. [10] designed a genetic clustering algorithm based on geometric shape to solve the problem, etc.

VRPFTW is special vehicle routing problem with soft time windows (VRPSTW). For VRPSTW, the starting service time can be out of customers’ ideal time windows. There are mainly two ways to constrain the starting service time. One is to set penalty cost for deviating time (the starting service time deviates from the ideal time window) [11, 12]. Another is to set customers’ satisfaction function for deviating time, which is VRPFTW. VRPFTW is mentioned in previous researches [13, 14], its essence is that, customers have subjective feelings for deviating time and the subjective feelings are vague. By introducing fuzzy theory [15], customers’ feelings towards deviating time can be quantified and defined as customers’ satisfaction. Because VRPFTW meets actual distribution situations closely, it’s widely concerned and researched. Jiafu et al. [16] studied non-linear customers’ satisfaction function. Seyed et al. [17] designed multi-objective model to solve dynamic vehicle routing problem (DVRP). Canhong et al. [18] designed an embedded hybrid neighborhood search algorithm. Wang et al. [19] improved solution algorithm of Cheng et al. [14]. Similarly, Jun et al. [20] studied that the charging operation is affected by divers’ feeling for the rate of electricity consuming. Mehdi et al. [21] constructed multi-objective model to solve multi-depot vehicle routing problem with fuzzy time windows (MDVRPFTW).

For above researches about VRPFTW, limited by the existing method to describe customer’s satisfaction, maximum time windows and ideal time windows of customers must be known in advance. However, there are some distribution situations without maximum time windows. Like replenishment problem (RP), supplementary goods aren’t used immediately and customers’ feelings towards maximum time windows are fuzzy, customers don’t set maximum time windows. How to describe customer’s satisfaction for VRPFTW without maximum time windows hasn’t been studied by anybody. This paper is focused on two aspects, one is combination of MDOVRP and VRPFTW, another is VRPFTW without maximum time windows. Contributions are as follows:

One method to describe customers’ satisfaction for MDOVRPFTW without maximum time windows was proposed.

One hybrid genetic algorithm for MDOVRPFTW without maximum time windows was designed. It is applied to one case and proved valid.

This paper is organized as follows. Section 2 presents the method and the model for MDOVRPFTW without maximum time windows. Section 3 presents the hybrid genetic algorithm. Section 4 presents one case that illustrates the model and algorithm. Some concluding remarks and possible avenues for future work are given in section 5.

2 The method and the model for MDOVRPFTW without maximum time windows

2.1 The existing method to describe customers’ satisfaction in VRPFTW

Ideal time windows and maximum time windows are known. Generally, customers’ satisfaction is expressed as follows [16 –18].

$L_{j} ({SST}_{j}) = {\begin{matrix} 0 & {SST}_{j} ⩽ E_{j} \\ ({SST}_{j} - E_{j}) / (e_{j} - E_{j}) & E_{j} < {SST}_{j} ⩽ e_{j} \\ 1 & e_{j} < {SST}_{j} ⩽ l_{j} \\ (L_{j} - {SST}_{j}) / (L_{j} - l_{j}) & l_{j} < {SST}_{j} ⩽ L_{j} \\ 0 & L_{j} < {SST}_{j} \end{matrix}$ (1)

${ET}_{j} (δ) = L_{j}^{- 1} (δ) = δ \cdot (e_{j} - E_{j}) + E_{j}, 0 ⩽ δ ⩽ 1$ (2a) ${LT}_{j} (δ) = L_{j}^{- 1} (δ) = L_{j} - δ \cdot (L_{j} - l_{j}), 0 ⩽ δ ⩽ 1$ (2b)

j represents customer, SST_j represents starting service time of j. [e_j, l_j] represents ideal time window of j. [E_j, L_j] represents maximum time window of j. L_j (SST_j) represents satisfaction of j at SST_j. δ represents minimum satisfaction of customers, it’s to prevent to lost customers because of low customers’ satisfaction. For Equations (2a) and (2b), [ET_j (δ) , LT_j (δ)] represents actual maximum time window of j that δ corresponds. The actual maximum time window is different from the maximum time window.

2.2 The method to describe customers’ satisfaction for MDOVRPFTW without maximum time windows

Maximum time windows are the boundaries between customers’ acceptance and rejection to goods. For MDOVRPFTW without maximum time windows, when the starting service time is out of the ideal time window, customers’ acceptance or rejection to goods is determined by DCWAG. Customers accept goods if DCWAG is high. DCWAG is affected by the starting service time, DCWAG is high if the starting service time is close to the ideal time window. The main steps of the method proposed are as follows.

Step 1: Refine and define the factors affecting DCWAG from MDOVRPFTW without maximum time windows and ensure the factors are independent to each other. Quantify the factors by some ways and ensure one factor is affected by the starting service time.

Step 2: DCWAG is one fuzzy concept. By setting weight coefficients for the factors, the function of DCWAG and the factors can be constructed, that is the function between DCWAG and the starting service time.

Step 3: Set an threshold for DCWAG as the boundary between customers’ acceptance and rejection to goods, the starting service time that the threshold corresponds is treated as the maximum time window. Using the existing methods to describe customers’ satisfaction for MDOVRPFTW without maximum time windows.

In this paper, three factors are refined and defined, the importance of goods towards customers, customers’ satisfaction towards starting service time, the coordination relationship between drivers and customers. Three factors are quantified by direct method [22] (one method to determine fuzzy membership function) and questionnaire survey. The rest content is shown in section 2.4 and 2.5.

2.3 Notation for MDOVRPFTW without maximum time windows

Parameters and variables in the mathematical model are shown in Table 1.

Table 1
Parameters and variables in the mathematical model

Indices:

i, j Node (depots and customers) index

k Vehicle index

Set:

D Set of depots

C Set of customers

V Set of all nodes (depots and customers), V = D ∪ C

K Set of all vehicles

Parameters:

t ₀ Earliest time of dispatching vehicles from depots

d _j Demand of customerj;∀j ∈ C

$d_{j}^{Max}$ Maximum demand of customer j in history; ∀j ∈ C

[e_j, l_j] Ideal time window of customer j;∀j ∈ C

$t_{j}^{Max}$ Maximum deviating time from ideal time window of customerjin history; ∀j ∈ C

$t_{j}^{serve}$ Time for serving customer j;∀j ∈ C

k ^sum Number of all vehicles

Q Maximum load of one vehicle

C ₁ Traveling cost for one kilometer

C ₂ Fixed cost for one vehicle

v Speed

χ Price for unit good

χ ⁰ Reference price for unit good, it’s decided according the situations

δ Threshold for customers’ satisfaction

ɛ Threshold for DCWAG

d _ij Distance from node i to node j,d_ij = d_ji;∀i, j ∈ V

t _ij Time from node i to node j.t_ij = d_ij/(w_ij · v),w_ijrepresents circuit coefficient, which is small when the traffic from node i to node jis convenient; ∀i, j ∈ V

M A large constant

Decision variables:

$x_{ij}^{k}$ It’s 1 if vehicle ktravels from node i to node j, 0 otherwise; ∀i, j ∈ V,∀k ∈ K

$y_{j}^{k}$ It’s 1 if vehicle k serves customer j, 0 otherwise; ∀j ∈ C,∀k ∈ K

$t_{j}^{k}$ Cumulative time that vehicle k travels from the depot to customer j, it’s composed of traveling time and serving time; ∀j ∈ C, ∀k ∈ K

U _jk An auxiliary (real) variable for sub-tour elimination constraint; ∀j ∈ V, ∀k ∈ K

2.4 Quantify three factors

The importance of goods towards customers is one fuzzy concept, it’s described referring to direct method [22]. For MDOVRPFTW without maximum time windows, generally goods needn’t to be used urgently and the importance of goods towards customers is mainly reflected in price and quantity. More goods customers order, higher the importance of goods towards customers is. Higher the price of good is, higher the importance of goods towards customers is. The importance of goods towards customers is expressed in Equation (3).

$A_{j} = (χ / χ^{0}) \cdot (d_{j} / d_{j}^{Max}), A_{j} \in [0, 1]$ (3)

A_j represents the importance of goods towards customer j. Generally $d_{j} \in (0, d_{j}^{Max}]$ , $d_{j} / d_{j}^{Max} \in (0, 1]$ , while $d_{j} > d_{j}^{Max}$ , $d_{j} / d_{j}^{Max} = 1$ . χ/χ⁰ ∈ (0, + ∞]. The relationship between the importance of goods towards customers and relative quantity is shown in Fig. 1.

Fig. 1

The relationship between the importance of goods towards customers and relative quantity.

When relative quantity is fixed, there’s positive correlation between the importance of goods towards customers and relative quantity. When relative price is large, the importance of goods towards customers varies fast with relative quantity.

Customers’ satisfaction towards starting service time is one fuzzy concept, it’s described referring to direct method. (Customers’ satisfaction towards starting service time is similar to customers’ satisfaction in VRPFTW, but two maximum time windows are different) It’s expressed in Equation (4).

$B_{j} = {\begin{matrix} 0 & {SST}_{j} ⩽ e_{j} - t_{j}^{Max} \\ ({SST}_{j} - e_{j} + t_{j}^{Max}) / t_{j}^{Max} & e_{j} - t_{j}^{Max} < {SST}_{j} ⩽ e_{j} \\ 1 & e_{j} < {SST}_{j} ⩽ l_{j} \\ (l_{j} + t_{j}^{Max} - {SST}_{j}) / t_{j}^{Max} & l_{j} < {SST}_{j} ⩽ l_{j} + t_{j}^{Max} \\ 0 & l_{j} + t_{j}^{Max} < {SST}_{j} \end{matrix}$ (4)

B_j represents customer j’s satisfaction towards SST_j.

The coordination relationship between drivers and customers is that, when uncertain situations affect terminal delivery of goods between drivers and customers, drivers and customers have different subjective intentions to complete delivery of goods. For example, outgoing drivers are better in completing delivery lately than incoming drivers, tolerant customers are likely to accept late delivery, etc. Questionnaire survey is introduced, the coordination relationship between drivers and customers is described by customers scoring drivers. Higher the scores are, better the coordination relationship between drivers and customers is. κ_kj represents the coordination relationship between driver (vehicle) k and customer j, κ_kj ∈ [0, 1].

2.5 The function between DCWAG and starting service time

DCWAG is positively correlated with above three factors. Three factors are independent to each other. Three factors are in [0, 1].

DCWAG is one fuzzy concept. The way to describe DCWAG is to set weight coefficients (P₁, P₂, P₃, P₁ + P₂ + P₃ = 1) for three factors (A_j, B_j, κ_kj). When starting service time of customer j is out of the ideal time window, as shown in Equation (5).

$DCWAG = P_{1} \cdot A_{j} + P_{2} \cdot B_{j} + P_{3} \cdot κ_{kj}, DCWAG \in (0, 1)$ (5)

The function between DCWAG and starting service time is expressed in Equation (6).

P_kj (SST_j) represents the degree of customer j’s wiliness to accept goods when vehicle k starts serving j at SST_j. P_kj (SST_j) ∈ [P₁ · A_j + P₃ · κ_kj, P₁ · A_j + P₂ + P₃ · κ_kj).

The mode for MDOVRPFTW without maximum time windows

$P_{kj} ({SST}_{j}) = {\begin{matrix} P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} & {SST}_{j} < e_{j} - t_{j}^{Max} \\ P_{1} \cdot A_{j} + P_{2} \cdot ({SST}_{j} - e_{j} + t_{j}^{Max}) / t_{j}^{Max} + P_{3} \cdot κ_{kj} & e_{j} - t_{j}^{Max} ⩽ {SST}_{j} < e_{j} \\ P_{1} \cdot A_{j} + P_{2} \cdot (l_{j} + t_{j}^{Max} - {SST}_{j}) / t_{j}^{Max} + P_{3} \cdot κ_{kj} & l_{j} < {SST}_{j} ⩽ l_{j} + t_{j}^{Max} \\ P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} & l_{j} + t_{j}^{Max} < {SST}_{j} \end{matrix}, j \in C, k \in K$ (6)

$E_{kj} (ɛ) = {\begin{matrix} e_{j} & P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} + P_{2} ⩽ ɛ < 1 \\ (ɛ - P_{1} \cdot A_{j} - P_{3} \cdot κ_{kj}) \cdot t_{j}^{Max} / P_{2} - t_{j}^{Max} + e_{j} & P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} < ɛ < P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} + P_{2} \\ the lower of the business time window of j & 0 < ɛ ⩽ P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} \end{matrix}, j \in C, k \in K$ (7a) $L_{kj} (ɛ) = {\begin{matrix} l_{j} & P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} + P_{2} ⩽ ɛ < 1 \\ l_{j} + t_{j}^{Max} - (ɛ - P_{1} \cdot A_{j} - P_{3} \cdot κ_{kj}) \cdot t_{j}^{Max} / P_{2} & P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} < ɛ < P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} + P_{2} \\ the upper of the business time window of j & 0 < ɛ ⩽ P_{1} \cdot A_{j} + P_{3} \cdot κ_{kj} \end{matrix}, j \in C, k \in K$ (7b)

$\begin{matrix} {ET}_{kj} (δ) = & δ \cdot (e_{j} - E_{kj} (ɛ)) + E_{kj} (ɛ), \\ 0 ⩽ δ ⩽ 1, j \in C, k \in K \end{matrix}$ (8a) $\begin{matrix} {LT}_{kj} (δ) = & L_{kj} (ɛ) - δ \cdot (L_{kj} (ɛ) - l_{j}), \\ 0 ⩽ δ ⩽ 1, j \in C, k \in K \end{matrix}$ (8b)

For Equations (7a) and (7b), [E_kj (ɛ) , L_kj (ɛ)] indicates maximum time window of customer j that ɛ corresponds when vehiclek serves customer j. For Equations (8a) and (8b), [ET_kj (δ) , LT_kj (δ)] indicates actual maximum time window of customer j that δ corresponds when vehicle k serves customer j.

$MinZ = C_{2} \cdot \sum_{i \in D} \sum_{j \in C} \sum_{k \in K} x_{ij}^{k} + C_{1} \cdot \sum_{i \in V} \sum_{j \in C} \sum_{k \in K} x_{ij}^{k} \cdot d_{ij}$ (9)

$\sum_{i \in D} \sum_{j \in D} \sum_{k \in K} x_{ij}^{k} = 0$ (10) $\sum_{i \in C} \sum_{j \in D} \sum_{k \in K} x_{ij}^{k} = 0$ (11) $U_{ik} - U_{jk} + | V | \cdot x_{ij}^{k} ⩾ | V | - 1 \forall i, j \in V \forall k \in K$ (12) $\sum_{i \in D} \sum_{j \in C} \sum_{k \in K} x_{ij}^{k} ⩽ k^{sum}$ (13) $\sum_{i \in V} \sum_{k \in K} x_{ij}^{k} = 1 \forall j \in C$ (14a) $\sum_{j \in C} \sum_{k \in K} x_{ij}^{k} ⩽ 1 \forall i \in C$ (14b)

$\sum_{i \in V} \sum_{j \in C} x_{ij}^{k} \cdot d_{j} ⩽ Q \forall k \in K$ (15) $t_{j}^{k} ⩽ t_{ij} + M \cdot (1 - x_{ij}^{k}) \forall i \in D \forall j \in C \forall k \in K$ (16a) $t_{j}^{k} ⩾ t_{ij} - M \cdot (1 - x_{ij}^{k}) \forall i \in D \forall j \in C \forall k \in K$ (16b)

$t_{j}^{k} ⩽ t_{i}^{k} + t_{ij} + s_{i}^{serve} + M \cdot (1 - x_{ij}^{k}) \forall i, j \in C \forall k \in K$ (17a) $t_{j}^{k} ⩾ t_{i}^{k} + t_{ij} + s_{i}^{serve} - M \cdot (1 - x_{ij}^{k}) \forall i, j \in C \forall k \in K$ (17b)

$t_{0} + t_{j}^{k} ⩽ L_{kj} (ɛ) + M \cdot (1 - y_{j}^{k}) \forall j \in C \forall k \in K$ (18) $\sum_{i \in V} \sum_{k \in K} x_{ij}^{k} = \sum_{k \in K} y_{j}^{k} \forall j \in C$ (19)

Algorithm 1 Decode
Input:N(one chromosome); k^sum(number of all vehicles)
Output: One decoding scheme.
1: Initialize: Add “0” at starting and ending positions of N
2: forj∈ { 1, 2, ·· · , k^sum - 1 } do
3: Find positions of all “0” inN
4: Find positions that we can insert “0” inN
5: Exclude positions by 3 from positions by 4
6: Insert “0” into someone position by 5 randomly and get newN
7: end for

Algorithm 2 Judge whether the decoding scheme meets actual maximum time windows of all customers
Input: One decoding scheme; k^sum (number of all vehicles); \|D\|(number of all depots)
Output: Judge whether the decoding scheme meets actual maximum time windows of all customers
1: Initialize: EnsureA (number of all routes of the decoding scheme); define a (route index) and a^route (route that route
index a corresponds); define k (vehicle index); define d (depot index);define E = zeros (A, k^sum)
2: fora∈ { 1, 2, ·· · , A } do
3: fork∈ { 1, 2, ·· · , k^sum } do
4: ford∈ { 1, 2, ·· · , \|D\| } do
5: if Vehicle k serves the route a^route from depot d, which meets actual maximum time windows of all customers then
6: E (a, k) = 0 and end for in step 4
7: elseE (a, k) = 1
8: end if
9: end for
10: end for
11: end for
12: Use Hungarian algorithm to judge whether the minimum sum of non-peer, non-identical data in E is 0
13: if The minimum sum is 0 then
14: The decoding scheme meets the actual maximum time window of each customer
15: else no
16: end if

Note: Step 5 in algorithm 2 is specified in algorithm 3.

The objective function is given by Equation (9), which includes two terms: one term is the cost of vehicles dispatched, another term is the cost of routes. Constraint (10) guarantees no distribution of goods between depots. Constraint (11) guarantees that depots are only treated as starting nodes. Constraint (12) is used for eliminating sub-tours in the vehicle routing. Constraint (13) guarantees that number of vehicles used can’t be larger than total number available. Constraints (14a) and (14b) guarantee that each customer must be serviced and once. Constraint (15) ensures that vehicles are not overweight. Constraints (16a) and (16b) are used for calculating cumulative time that vehicles arrive first customers. Constraints (17a) and (17b) are used for calculating cumulative time of adjacent customers i and j covered by vehicle k. Constraint (18) guarantees that the scheme can meet actual maximum time windows of all customers that ɛ corresponds. Constraint (19) indicates the relationship between variables $x_{ij}^{k}$ and $y_{j}^{k}$ .

3 The hybrid genetic algorithm

About VRP, there are genetic algorithm, ant colony algorithm, particle swarm optimization algorithm and tabu search algorithm, etc. [23 –25] In this paper, the hybrid genetic algorithm designed follows the framework of genetic algorithm (coding, decoding, selecting, crossing and mutating) and Hungarian algorithm is embedded in it. Meanwhile, some links of the hybrid genetic algorithm were optimized. The main steps are as follows.

3.1 Code

Customers, depots and vehicles are respectively represented by real number code. One chromosome is generated by arraying all customers randomly.

3.2 Decode and construct initial population

Firstly, add “0” at starting and ending positions of one chromosome. Then randomly select locations and insertk^sum - 1“0”. Finally, customers are grouped without violating total number of vehicles. While inserting “0”, that “0” is adjacent to “0” will generate duplicate decoding schemes, as shown in Fig. 2. The corresponding decoding steps are shown in algorithm 1.

Fig. 2

The situation with duplicate decoding schemes.

Algorithm 3 Judge whether actual maximum time windows of all customers can be met, when vehicle k serves the route a^route from depot d
Input: a^route(route that route index a corresponds); [ET_jk (δ) , LT_jk (δ)] (the actual maximum time window that δ corresponds when vehicle k serves customerj;t₀ (earliest time of dispatching vehicles from depots); d (depot)
Output: Judge whether actual maximum time windows of all customers can be met, when vehicle k serves the route a^route from depot d
1: Initialize: Get \|a\| (number of customers in a^route); s = 0; define t_bk (cumulative time that vehicle k traveling
from d to node b)
2: forb∈ { a^route (1) , a^route (2) , ·· · , a^route (\|a\|) } do
3: ift_bk ⩽ LT_bk (δ) - t₀then
4: s = s + 1.
5: elses = s.
6: end if
7: end for
8: forf∈ { 1, 2, ·· · , \|a\| - 1 } do
9: b = a^route (f)
10: forc∈ { a^route (f + 1) , a^route (f + 2) , ·· · , a^route (\|a\|) } do
11: ift_ck - t_bk ⩽ LT_ck (δ) - ET_bk (δ) then
12: s = s + 1.
13: elses = s.
14: end if
15: end for
16: end for
17: ifs = = \|a\| · (\|a\| + 1)/2 then
18: E (a, k) = 0, actual maximum time windows of all customers can be met
19: elseE (a, k) = 1, actual maximum time windows of all customers can’t be met
20: end if

Algorithm 4 Get one initial solution
Input: N(one chromosome); M (number of decoding schemes of N)
Output: One initial solution
1: Initialize: Define m(decoding scheme index)
2: form∈ { 1, 2, ·· · , M }do
3: DecodeN
4: Calculate total of cost of vehicles used and cost of routes for N as fitness value
5: if Decoding scheme violates Q (maximum load of vehicle) then
6: Fitness value becomes one sufficiently large positive number
7: end if
8: if Decoding scheme can’t meet actual maximum time windows of all customers then
9: Fitness value becomes one sufficiently large positive number
10: end if
11: end for
12: Choose decoding scheme with lowest fitness value from M decoding schemes as one initial solution

Algorithm 1 can meet constraints (10) — (14). To meet constraints (16), (17) and (18), the main decoding steps are shown in algorithm 2, 3 and 4.

Algorithm 4 can also meet constraint (15). Repeatedly run algorithm 3 to get the initial population.

3.3 Select

Tournament selection strategy and elitist preservation strategy are applied to the algorithm. Tournament selection strategy is to choose one with better fitness value as progeny chromosome from any two paternal chromosomes. Because tournament selection strategy possibly retains chromosomes with poor fitness value, elitist preservation strategy is to choose the best chromosome from paternal population instead of the worst chromosome from progeny population.

3.4 Cross

Order-based crossover is adopted. As shown in Fig. 3.

Fig. 3

Order-based crossover.

Step 1: Select and retain genes “345” from position 3 to 5 of primitive chromosome A.

Step 2: Eliminate genes “3”, “4” and “5” of primitive chromosome B.

Step 3: In turn, put rest genes “7”, “8”, “2”, “1”and “6” of primitive chromosome B at position 1, 2, 6, 7 and 8 of primitive chromosome A, new chromosome A is generated. In the same way, new chromosome B is generated.

3.5 Mutate

Reciprocal exchange operator, insertion operator and inversion operator are adopted. As shown in turn in Fig. 4. Reciprocal exchange operator is to randomly select two genes from one chromosome and exchange positions of two genes, as shown in Fig. 4(a). Insertion operator is to randomly select two genes from one chromosome and put one gene behind anther, as shown in Fig. 4(b). Inversion operator is to randomly select two genes from one chromosome and reverse the order of genes between the two genes, as shown in Fig. 4(c).

Fig. 4

The mutation operator.

4 Numerical example

4.1 Background of the example

The pharmaceutical logistics company is composed by 4 depots and 20 rental vehicles. The medicines are delivered to some hospitals or drugstores in Beijing and Hebei province. The medicines that hospitals or drugstores request are generally used to replenish stocks and the medicines in stocks can be used up after 2-3 days without replenishment. Hospitals or drug-stores probably reject or accept the medicines out of ideal time windows (hospitals are open for 24 hours). Relevant parameters of the example are adjusted, some parameters are shown in Table 2, Table 3 and Table 4.

Table 2
Parameters of depots and customers

Nodes q _j ST_j(h) [e_j, l_j]

Depot 1 0 0 [0:00,19:00]

Depot 2 0 0 [0:00,19:00]

Depot 3 0 0 [0:00,19:00]

Depot 4 0 0 [0:00,19:00]

Customer 1 1.2 0.36 [6:30, 10:00]

Customer 2 0.8 0.24 [7:30, 16:30]

Customer 3 1.6 0.48 [8:30, 15:30]

Customer 4 0.5 0.15 [7:30, 16:30]

Customer 5 1.2 0.36 [8:30, 13:30]

Customer 6 0.5 0.15 [7:30, 16:30]

Customer 7 1.3 0.39 [6:30, 9:30]

Customer 8 2 0.6 [8:30, 13:30]

Customer 9 1.3 0.39 [6:30, 10:00]

Customer 10 1.8 0.54 [7:30, 16:30]

Customer 11 0.7 0.21 [8:30, 13:30]

Customer 12 0.6 0.18 [7:30, 16:30]

Customer 13 0.9 0.27 [7:30, 15:30]

Customer 14 0.9 0.27 [6:30, 10:30]

Customer 15 0.4 0.12 [6:30, 9:30]

Customer 16 2.5 0.75 [7:30, 16:30]

Customer 17 0.5 0.15 [7:30, 15:30]

Customer 18 1.7 0.51 [6:30, 9:30]

Customer 19 0.3 0.09 [7:30, 16:30]

Customer 20 1.6 0.48 [6:30, 12:00]

Customer 21 2.5 0.75 [7:30, 12:00]

Customer 22 2.1 0.63 [7:30, 16:30]

Customer 23 1.4 0.42 [7:30, 16:30]

Customer 24 1.9 0.57 [6:30, 16:30]

Customer 25 1.4 0.42 [8:30, 13:30]

Customer 26 0.6 0.18 [6:30,12:00]

Customer 27 1.6 0.48 [7:30,16:30]

Customer 28 0.9 0.27 [8:30,13:30]

Customer 29 2 0.6 [7:30,16:30]

Customer 30 1.3 0.39 [6:30, 12:00]

Customer 31 1 0.3 [7:30, 15:30]

Customer 32 1.6 0.48 [6:30, 10:30]

Customer 33 1.9 0.57 [6:30, 9:30]

Customer 34 2.2 0.66 [7:30, 12:00]

Customer 35 1.4 0.42 [7:30, 16:30]

Customer 36 1 0.3 [8:30, 14:00]

Customer 37 1.1 0.33 [7:30, 15.30]

Customer 38 1.5 0.45 [7:30, 16:30]

Customer 39 1.3 0.39 [6:30, 12:00]

Customer 40 1.5 0.45 [7:30, 16:30]

Customer 41 0.8 0.24 [7:30, 16:30]

Customer 42 2.2 0.66 [7:30, 11:30]

Customer 43 2.4 0.72 [7:30, 16:30]

Customer 44 0.3 0.09 [7:30, 16:30]

Customer 45 2.5 0.75 [7:30, 11:30]

Customer 46 1.9 0.57 [6:30, 13:30]

Customer 47 2.1 0.63 [7:30, 12:00]

Customer 48 1 0.3 [7:30, 16:30]

Nodes	q _j	ST_j(h)	[e_j, l_j]
Depot 1	0	0	[0:00,19:00]
Depot 2	0	0	[0:00,19:00]
Depot 3	0	0	[0:00,19:00]
Depot 4	0	0	[0:00,19:00]
Customer 1	1.2	0.36	[6:30, 10:00]
Customer 2	0.8	0.24	[7:30, 16:30]
Customer 3	1.6	0.48	[8:30, 15:30]
Customer 4	0.5	0.15	[7:30, 16:30]
Customer 5	1.2	0.36	[8:30, 13:30]
Customer 6	0.5	0.15	[7:30, 16:30]
Customer 7	1.3	0.39	[6:30, 9:30]
Customer 8	2	0.6	[8:30, 13:30]
Customer 9	1.3	0.39	[6:30, 10:00]
Customer 10	1.8	0.54	[7:30, 16:30]
Customer 11	0.7	0.21	[8:30, 13:30]
Customer 12	0.6	0.18	[7:30, 16:30]
Customer 13	0.9	0.27	[7:30, 15:30]
Customer 14	0.9	0.27	[6:30, 10:30]
Customer 15	0.4	0.12	[6:30, 9:30]
Customer 16	2.5	0.75	[7:30, 16:30]
Customer 17	0.5	0.15	[7:30, 15:30]
Customer 18	1.7	0.51	[6:30, 9:30]
Customer 19	0.3	0.09	[7:30, 16:30]
Customer 20	1.6	0.48	[6:30, 12:00]
Customer 21	2.5	0.75	[7:30, 12:00]
Customer 22	2.1	0.63	[7:30, 16:30]
Customer 23	1.4	0.42	[7:30, 16:30]
Customer 24	1.9	0.57	[6:30, 16:30]
Customer 25	1.4	0.42	[8:30, 13:30]
Customer 26	0.6	0.18	[6:30,12:00]
Customer 27	1.6	0.48	[7:30,16:30]
Customer 28	0.9	0.27	[8:30,13:30]
Customer 29	2	0.6	[7:30,16:30]
Customer 30	1.3	0.39	[6:30, 12:00]
Customer 31	1	0.3	[7:30, 15:30]
Customer 32	1.6	0.48	[6:30, 10:30]
Customer 33	1.9	0.57	[6:30, 9:30]
Customer 34	2.2	0.66	[7:30, 12:00]
Customer 35	1.4	0.42	[7:30, 16:30]
Customer 36	1	0.3	[8:30, 14:00]
Customer 37	1.1	0.33	[7:30, 15.30]
Customer 38	1.5	0.45	[7:30, 16:30]
Customer 39	1.3	0.39	[6:30, 12:00]
Customer 40	1.5	0.45	[7:30, 16:30]
Customer 41	0.8	0.24	[7:30, 16:30]
Customer 42	2.2	0.66	[7:30, 11:30]
Customer 43	2.4	0.72	[7:30, 16:30]
Customer 44	0.3	0.09	[7:30, 16:30]
Customer 45	2.5	0.75	[7:30, 11:30]
Customer 46	1.9	0.57	[6:30, 13:30]
Customer 47	2.1	0.63	[7:30, 12:00]
Customer 48	1	0.3	[7:30, 16:30]

Table 3

Parameters of vehicles

Parameters	Q	k ^sum	C ₁	C ₂	v	w _ij
			(¥/km)	(¥)	(km/h)
Value	15	20	5	600	60	1

Table 4

Other parameters

_Parameters	Value
χ/χ⁰	1
$Q_{j}^{Max}$	1.35
$T_{j}^{Max}$	1 hour
κ _kj	One random value in [0:0.1:1],
P₁, P₂, P₃	Weights of three factors are 1/3

Notes: for $Q_{j}^{Max}$ , 1.35 is the median of 48 customers’ demands, the reason to set 1.35 is that some customers’ demands are fixed, whose relative quantities are 1.

Notes: Q represents that medicines are packed in unified boxes and one vehicle can be loaded into 15 boxes at most. k^sum,C₁ and C₂ are provided by the company. In reality, v and w_ij aren’t fixed. Aiming at adopting to the model and the algorithm (meanwhile fixing v and w_ij), we justified the real distances (between customers, between customers and depots) to ensure that the driving time (between customers, between customers and depots) meets to the actual situations. The justified distances aren’t equal to the real distances.

4.2 Solution

The population and number of iterations of the algorithm are set as 50 and 150, respectively. The crossover and mutation rates are set as 0.8 and 0.15, respectively. All the algorithms are coded in MATLAB and run on a computer with Intel Core i7-700 3.60 GHz and 16 GB RAM operating on Microsoft Windows 10 operation system. When ɛ is 0.6 and δ is 0.2, by running 10 times, the scheme with lowest cost is selected. The data are shown in Table 5, routes of the scheme with lowest cost are shown in Fig. 5.

Table 5
The scheme with lowest cost

Route index Routes

1 ①-46-21-47-42

2 ③-14-29

3 ①-32-19-28-5-8-13

4 ②-22-37-9

5 ②-10-6-7-44-4-30-43-2-31

6 ④-39-12-18-36

7 ①-35-15-38-26-24-34-3-27

8 ①-41-48-11-45-16-40-25-23-27

9 ③-20-33-1

Routes cost(¥) 8862

Vehicles cost(¥) 5400

Total cost(¥) 14262

Route index	Routes
1	①-46-21-47-42
2	③-14-29
3	①-32-19-28-5-8-13
4	②-22-37-9
5	②-10-6-7-44-4-30-43-2-31
6	④-39-12-18-36
7	①-35-15-38-26-24-34-3-27
8	①-41-48-11-45-16-40-25-23-27
9	③-20-33-1
Routes cost(¥)	8862
Vehicles cost(¥)	5400
Total cost(¥)	14262

Notes: ①, ②, ③ and ④ represent depots.

Fig. 5

Routes of the scheme with lowest cost.

The relationship between routes and vehicles is shown in Table 6. “0” represents that the vehicle serving the route can meet actual maximum time windows of all customers in the route, “1” represents the vehicle can’t serve the route.

Table 6

The relationship between routes and vehicles

k\s	1	2	3	4	5	6	7	8	9
1	0	0	0	0	0	0	0	0	0
2	0	0	0	0	0	0	0	0	0
3	0	0	1	0	0	0	0	0	0
4	0	0	0	0	0	0	0	0	0
5	0	0	1	0	0	0	0	0	0
6	0	0	0	0	0	0	0	0	0
7	0	0	1	0	0	0	0	0	0
8	0	0	0	0	0	0	0	0	0
9	0	0	0	0	0	0	0	0	0
10	0	0	0	0	0	0	0	0	0
12	0	0	0	0	0	0	0	0	0
13	0	0	1	0	0	0	0	0	0
14	0	0	0	0	0	0	0	0	0
15	0	0	1	0	0	0	0	0	0
16	0	0	0	0	0	0	0	0	0
17	0	0	1	0	0	0	0	0	0
18	0	0	0	0	0	0	0	0	0
19	0	0	0	0	0	0	0	0	0
20	0	0	0	0	0	0	0	0	0

Notes: k represents vehicle index,srepresents route index.

4.3 Prove rationality of solutions by the hybrid genetic algorithm

The solutions by the hybrid genetic algorithm are compared with the exact solutions. Firstly, 4 examples (A, B, C and D) are set. Example A is aimed at serving front 5 customers in Table 2. Example B is aimed at serving front 10 customers in Table 2. Example C is aimed at serving front 20 customers in Table 2. Example D is aimed at serving all 48 customers in Table 2. When ɛ is 0.9 and δ is 0.2, it’s to run 10 times by the hybrid genetic algorithm for each example, the scheme with lowest cost is selected as the solution by the hybrid genetic algorithm. As follows in Table 7.

Table 7
Comparison of exact solutions and solutions by the hybrid genetic algorithm

Examples A B C D

Solutions by the hybrid algorithm (¥) 1322 2015 4426 14198

Exact solutions(¥) 1322 2006 4215 13635

Deviation degree 0 0.45% 5.0% 4.1%

Examples	A	B	C	D
Solutions by the hybrid algorithm (¥)	1322	2015	4426	14198
Exact solutions(¥)	1322	2006	4215	13635
Deviation degree	0	0.45%	5.0%	4.1%

Notes: Deviation degree is the difference between the solution by the hybrid algorithm and the exact solution divided by the exact solution.

Table 6 shows that when the scale of customers is small, the solutions by the hybrid genetic algorithm are equal to the exact solutions. When the scale of customers is large, the solutions by the hybrid genetic algorithm approach the exact solutions with slight deviation.

4.4 Prove stability of solutions by the hybrid genetic algorithm

When δ is 0.2, it’s to run 10 times respectively for differentɛ. Maximum deviation degree and average deviation degree are obtained. Maximum deviation degree is the difference between maximum cost for 10 times and minimum cost for 10 times divided by minimum cost for 10 times. Average deviation degree is the average of deviation degrees for 10 times. As shown in Table 8.

Table 8
Stability of solutions by the hybrid genetic algorithm

ɛ Maximum deviation degree Average deviation degree

0.1 0.064 0.030

0.2 0.071 0.040

0.3 0.072 0.042

0.4 0.085 0.064

0.5 0.096 0.069

0.6 0.074 0.041

0.7 0.074 0.044

0.8 0.092 0.056

0.9 0.072 0.041

ɛ	Maximum deviation degree	Average deviation degree
0.1	0.064	0.030
0.2	0.071	0.040
0.3	0.072	0.042
0.4	0.085	0.064
0.5	0.096	0.069
0.6	0.074	0.041
0.7	0.074	0.044
0.8	0.092	0.056
0.9	0.072	0.041

Table 8 shows that average deviation degrees are within 6.9% and maximum deviation degrees are within 9.6%, which indicates solutions by the hybrid genetic algorithm are close to one small value and fluctuate in small range. To a certain extent, solutions by the hybrid genetic algorithm are stable.

4.5 Analysis for the coordination relationship between drivers and customers

Firstly, 5 examples (A, B, C, D and E) are set. The coordination relationship between drivers and customers in example C is one random value in [0:0.1:1]. That in example A is one value equal to the random value in example C minus 1. That in example B is one value equal to the random value in example C minus 0.5. That in example D is one value equal to the random value in example C plus 0.5, That in example B is one value equal to the random value in example C plus 1 (While values are less than 0, values are 0. While values are more than 1, values are 1.). When δ is 0.2, it’s to run 10 times for each example and the scheme with lowest cost is selected. As shown in Table 9 and Fig. 6.

Fig. 6

Influence of the coordination relationship between drivers and customers on cost.

Table 9

The cost of 5 examples

ɛ	A	B	C	D	E
0.4	14297	13815	13637	13490	13122
0.5	14510	14193	13858	13644	13453
0.6	14969	14416	14262	13822	13630

Notes: The unit is ¥.

Figure 6 reveals that improving the coordination relationship between drivers and customers can reduce distribution cost. When the coordination relationship between drivers and customers ranges from A to E, cost reduces 7.9–9.8%. The suggestion: when maximum time windows of customers are unknown and customers may refuse goods out of ideal time windows, logistics enterprises can adopt piecework wages or award drivers with high turnover, which can improve the desire that drivers complete delivery task and drivers’ attitude to serve customers. Finally logistics enterprises can improve DCWAG and reduce cost.

5 Conclusions

In this paper, one method to deal with MDOVRPFTW without maximum time windows is presented. By direct method and questionnaire method, the factors affecting DCWAG can be quantified. The function between DCWAG and the starting service time can be constructed. By setting an threshold for DCWAG, maximum time windows meeting actual situations can be determined. Meanwhile, one hybrid genetic algorithm is designed. Applying it to one example proves that the hybrid genetic algorithm is feasible. Some suggestions to improve DCWAG are proposed.

Note that factors affecting DCWAG are different at different situations, which means the factors aren’t fixed. The disadvantage of the proposed method is that the performance to optimize of the hybrid genetic algorithm is insufficient, a relatively better solution can be determined by running many times. The next research is to further improve the hybrid genetic algorithm.

6 Founding statement

This research was supported by the Emergency Management System Construction Research Special Project of the National Social Science Fund of China (20VYJ024).

References

Prins

, A simple and effective evolutionary algorithm for the vehicle routing problem, Computers and Operations Re-search 31(12) (2004), 1985–2002.

Tong

, Zhou

L.S.

, Liu

J.T.

and Zhou

X.S.

, Customized bus service design for jointly optimizing passenger-to-vehicle assignment and vehicle routing, Transportation Research Part C 85 (2017), 451–475.

Bae

and Chung

, Heuristics for Two Depot Heterogeneous Unmanned Vehicle Path Planning to Minimize Maximum Travel Cost, Sensors 19(11) (2019), 2461.

Ramakrishnan

, Ramanathan

and Ko

L.W.L.

, Adoption of RFID technologies in UK logistics: Moderating roles of size, barcode experience and government support, Expert Systems With Applications 41(1) (2014), 230–236.

Wang

, Assogba

, Fan

J.X.

, Xu

M.Z.

, Liu

and Wang.

H.Z.

, Multi-depot vehicle routing problem with shared transportation resource: Integration of time-dependent speed and piecewise penalty cost, Journal of Cleaner Production 232 (2019), 12–29.

Lim

and Wang

, Multi-Depot Vehicle Routing Problem: A One-Stage Approach, IEEE Transactions on Automation Science and Engineering 2(4) (2005), 397–402.

Salhi

, Imran

and Wassan

N.A.

, The multi-depot vehicle routing problem with heterogeneous vehicle fleet: Formulation and a variable neighborhood search implementation, Computers & Operations Research 52 (2013), 315–325.

, Y

Z.Z.

and Xie

J.X.

, A parallel improved ant colony optimization for multi-depot vehicle routing problem, The Journal of the Operational Research Society 62 (1) (2011), 183–188.

Rahimi

V.A.

, Gabriel

C.T.

and Gendreau

, Fleet-sizing for multi-depot and periodic vehicle routing problems using a modular heuristic algorithm, Computers & Operations Research 53 (2015), 9–23.

10.

Yucenur

G.N.

and Demirel

N.C.

, A new geometric shape-based genetic clustering algorithm for the multi-depot vehicle routing problem, Expert Systems with Applications 38(9) (2011), 11859–11865.

11.

Chiu

H.N.

, Lee

Y.S.

and Chang

J.H.

, Two approaches to solving the multi-depot vehicle routing problem with time windows in a time-based logistics environment, Production Planning & Control 17(5) (2006), 480–493.

12.

Keskin

, Laporte

and Catay

, Electric Vehicle Routing Problem with Time-Dependent Waiting Times at Recharging Stations, Computer & Operation Research 107 (2019), 77–94.

13.

Wang

H.F.

and Wen

Y.P.

, Time-constrained Chinese postman problems, Computers & Mathematics with Applications 44(3-4) (2002), 375–387.

14.

Cheng

R.W.

, Gen

and Tozawa

, Vehicle routing problem with fuzzy due-time using genetic algorithms, Japan Society for Fuzzy Theory and Systems 7(5) (1995), 1050–1061.

15.

J.P.

, Yan

and Li

, Vehicle routing optimization with soft time windows in a fuzzy random environment, Transportation Research Part E Logistics & Transportation Review 47(6) (2011), 1075–1091.

16.

Tang

J.F.

, Pan

Z.D.

, Fung

R.Y.K.

and Lau

, Vehicle routing problem with fuzzy time windows, Fuzzy Sets and Systems 160(5) (2009), 683–695.

17.

Ghannadpour

S.F.

, Noori

, Tavakkoli

M.R.

and Ghoseiri

, A multi-objective dynamic vehicle routing problem with fuzzy time windows: Model, solution and application, Applied Soft Computing Journal 14 (2014), 504–527.

18.

Lin

, Choy

K.L.

, Ho

G.T.S.

, Lam

H.Y.

, Pang

G.K.H.

and Chin.

K.S.

, A decision support system for optimizing dynamic courier routing operations, Expert Systems with Applications 41(15) (2014), 6917–6933.

19.

X.P.

and Zhang

, D.L., Y Disruption Management of Real-time Vehicle Routing Problem Based on Fuzzy Time Window, Journal of Systems & Management 19(6) (2010), 633–639.

20.

Yang

and Sun

, Battery swap station location-routing problem with capacitated electric vehicles, Computers & Operations Research 55 (2015), 217–232.

21.

Adelzadeh

, Asl

V.M.

and Koosha

, A mathematical model and a solving procedure for multi-depot vehicle routing problem with fuzzy time window and heterogeneous vehicle, The International Journal of Advanced Manufacturing Technology 75(5-8) (2014), 793–802.

22.

Lixin

, Fuzzy Engineering, monograph, Tsinghua University press, 2008, pp. 315–317.

23.

Liu

, Jiang

and Geng

, A hybrid genetic algorithm for the multi-depot open vehicle routing problem, OR Spectrum 36(2) (2014), 423–424.

24.

Y.F.

, Li

Z.M.

, Yan

and Feng

C.Y.

, A hybrid priority-based genetic algorithm for simultaneous pickup and delivery problems in reverse logistics with time windows and multiple decision-makers, Soft Computing 23(15) (2019), 6697–6714.

25.

Baniamerian

, Bashiri

and Tavakkoli