Dynamic multi-objective optimization for multi-period emergency logistics network

Abstract

In order to solve the problem of multi-period material supply in emergency logistics network, a dynamic multi-objective optimization mathematical model with constraints is constructed. The model takes the minimum cost and the maximum fill rate of demands as the objectives, and takes the location of distribution centers and the allocation of material as the decision variables. A dynamic self-adaptive multi-objective differential evolution algorithm is proposed to solve the mathematical model, and the feasible non-dominated solutions of the model are obtained. In the improved algorithm, on the one hand, a new environment change detect operator and a new environment change response strategy are adopted so that the traditional static optimization algorithm can be used to solve the dynamic optimization problem. On the other hand, the improved algorithm adopts adaptive mutation strategy to improve the ability of global exploration and local exploitation. Case study shows that the improved strategy greatly improves the performance of the algorithm, and can solve the dynamic multi-objective optimization problem effectively.

Keywords

Emergency logistics network multi objective optimization dynamic optimization evolutionary algorithm self-adaption

1 Introduction

In the past few decades, the emergency logistics network optimization has drawn a great deal of attention from researchers. In these works, on the one hand, emergency logistics network optimization can be divided into single-objective problem and multi-objective optimization problem based on the number of objective functions. On the other hand, the emergency logistics optimization problems can also be divided into single-period and multi-period optimization problems based on the stage of tasks. Caunhye et al. summarized the key problems as location problem, distribution problem, and transportation problem in emergency logistics network [1]. They stressed that, different from conventional logistics network optimization, the timeliness and uncertainty of emergency logistics optimization are more prominent. Li et al. studied the location problem of emergency resource allocation, and the static single-objective optimization problem was established with the objective of minimum number of open facilities. The model in their paper was solved by a heuristic algorithm [2]. Han et al. studied a large-scale emergency material supply problem. A dynamic single-objective optimization model was established, and the linear programming was used to solve the model [3]. Hu et al. established a static single-objective optimization model to solve the allocation problem in earthquake emergency shelters [4]. At the same time, Mohammadi et al. established a static multi-objective optimization model with the goal of satisfying the material demands, and used particle swarm optimization to solve the problem [5]. Zhao et al. established a three-objective optimization model with the objectives of minimum distance, minimum response time and maximum satisfaction degree. Multi-objective evolutionary algorithm was used to solve the model [6]. Wang et al. considered the problem of multi-scenes emergency logistics optimization, and established a single-target static optimization model which was solved by ant colony algorithm [7]. Dan et al. established a multi-objective dynamic optimization model for multi-stage emergency resource supply, and used ant colony algorithm to solve the problem [8]. Shen et al. considered the problem of transportation planning in case of emergency with a multi-objective and multi-stage optimization model [9]. Chang et al. established a static multi-objective optimization model for emergency logistics optimization, and used genetic algorithm to solve the problem [10].

It can be seen from the related research works that the multi-objective optimization problem is more difficult to solve than single-objective models. However, a multi-objective optimization is more suitable for emergency logistics optimization since that the emergency logistics needs to take into account the factors of costs and material satisfaction simultaneously. Because emergency material usually need to be supplied in multiple periods in the real world, the multi-period optimization is more suitable for emergency logistics optimization. It also can be see that, most of the current research works just focused on the dynamic optimization with changed demand, and few literatures have considered the changes of environmental factors such as transportation costs, supply routes, and so on. However, the changes of these parameters lead the objectives and constraints of model to be time-varying, and the solution will be more difficult to be achieved.

As stated above, the location-allocation joint optimization problem of emergency logistics in this paper considers multiple stages, multiple objectives and multiple constraints. The mathematical model takes the minimum supply costs and the maximum demand fill rate as objectives. The proposed mixed integer nonlinear model is difficult to be solved by traditional exact algorithm, thus an improved dynamic multi-objective evolutionary algorithm is also proposed.

2 Mathematic model of dynamic emergency logistics network optimization

2.1 Modeling description and assumptions

A three-echelon emergency logistics network is adopted, as shown in Fig. 1. The main functions of the nodes in logistics network are as follows:

Fig.1

Three-echelon emergency logistics network.

Supply centers, composed of central warehouses, are used to reserve emergency materials and supply them to the next level of nodes according to demands.

Distribution centers, undertake the tasks of material storage and transfer, accept the emergency materials from supply centers and then transport them to shelters. At the same time, the rest of materials in current period are stored in the distribution centers until the next period of supply begins.

In the real world, it is often not necessary to open all the distribution centers. Thus, in Fig. 1, solid line represents an opened node, and dashed line represents a closed node. If a distribution center is closed in a period, the input and output of material are 0. The material supply scheme proposed in this paper needs to determine which distribution centers are opened and which ones are closed.

Shelters, composed of the emergency shelters in affected areas, need a large number of relief materials such as food, medicine, tents, and so on. In our model, the demands of relief materials are different in each period.

The purpose of this model is to determine which distribution centers are open and determine the quantity of materials supplied between nodes.

The dynamic emergency logistics optimization model is on the base of the following assumptions: (1) The transportation costs between nodes are known and different in each period; (2) The opening cost, inventory cost and maximum capacity of distribution centers are all known and fixed; (3) The material demands and the shortage costs of each shelter are known, which are different in each period.

2.2 Parameters notation

The following notations are defined to formulate the proposed model:

Sets

I: Set of supply centers, i = 1, 2, ⋯ , I

J: Set of distribution centers, j = 1, 2, ⋯ , J

K: Set of shelters, k = 1, 2, ⋯ , K

T: Set of periods, t = 1, 2, ⋯ , T

Input variables

$C_{ijt}^{trans}$ : Transportation cost between supply center i and distribution center j in period t

$C_{jkt}^{trans}$ : Transportation cost between distribution center j and shelter k in period t

$C_{j}^{open}$ : Opening cost of distribution center j

$C_{j}^{invent}$ : Inventory cost of distribution center j

$C_{kt}^{short}$ : Shortage cost of shelter k in period t

$d_{k}^{t}$ : Material demand of shelter k in period t

U_j: Maximum capacity of distribution centers j

Decision variables

$x_{ij}^{t}$ : Quantity of material supply from supply center i to distribution center j in period t

$z_{jk}^{t}$ : Quantity of material supply from distribution center j to shelter k in period t

$y_{j}^{t}$ : Binary variable: $y_{j}^{t} = 1$ , if distribution center j is open in period t; $y_{j}^{t} = 0$ , if distribution center j is open in period t.

2.3 Model formulation

The model has two objectives, where objective 1 minimizes the total cost of supply: $min C_{t} = C_{t}^{open} + C_{t}^{trans} + C_{t}^{invent} + C_{t}^{short}$ (1) where $C_{t}^{open}$ is the opening cost, $C_{t}^{trans}$ is the transportation cost, $C_{t}^{invent}$ is the inventory cost, $C_{t}^{short}$ is the shortage cost. The formulas of thus costs are as follows: $C_{t}^{open} = \sum_{j \in J} C_{j}^{open} \cdot y_{j}^{t}, t = 1, 2, \dots, T$ (2)

$\begin{matrix} C_{t}^{trans} = \sum_{i \in I} \sum_{j \in J} C_{ijt}^{trans} \cdot x_{ij}^{t} + \sum_{j \in J} \sum_{k \in K} C_{jkt}^{trans} \cdot z_{jk}^{t}, \\ t = 1, 2, \dots, T \end{matrix}$ (3)

$\begin{matrix} C_{t}^{invent} = \sum_{j \in J} C_{j}^{invent} \cdot (\sum_{v = 1}^{t} (\sum_{i \in I} x_{ij}^{v} - \sum_{k \in K} z_{jk}^{v})), \\ t = 1, 2, \dots, T \end{matrix}$ (4) $C_{t}^{short} = \sum_{k \in K} C_{kt}^{short} \cdot (d_{k}^{t} - \sum_{j \in J} z_{jk}^{t}), t = 1, 2, \dots, T$ (5)

Objective 2 maximizes the fill rate of material demands: $max \frac{\sum_{k \in K} \sum_{j \in J} z_{jk}^{t}}{\sum_{k \in K} d_{k}^{t}}, t = 1, 2, \dots, T$ (6)

In order to transform objective (2) into a minimum value problem, the inverse of objective (2) is taken: $\min \frac{\sum_{k \in K} d_{k}^{t}}{\sum_{k \in K} \sum_{j \in J} z_{j k}^{t}}, t = 1, 2, \dots, T$ (7)

As can be seen, the above two objective functions conflict with each other. The larger quantity of material supplied, the larger the formula (1) is, and the smaller the formula (7) will be; the smaller quantity of material supplied, the smaller the formula (1) is, and the larger the formula (7) will be.

The constraints are as follows.

s.t. $\sum_{i \in I} x_{ij}^{t} ⩽ y_{j}^{t} \cdot U_{j}, t = 1, 2, \dots, T$ (8) $\sum_{k \in K} z_{jk}^{t} ⩽ y_{j}^{t} \cdot U_{j}, t = 1, 2, \dots, T$ (9) $\sum_{k \in K} z_{jk}^{t} ⩽ \sum_{i \in I} x_{ij}^{t}, t = 1$ (10)

$\begin{matrix} \sum_{k \in K} z_{jk}^{t} ⩽ \sum_{i \in I} x_{ij}^{t} + \sum_{v = 1}^{t - 1} (\sum_{i \in I} x_{ij}^{v} - \sum_{k \in K} z_{jk}^{v}), \\ t = 2, 3, \dots, T \end{matrix}$ (11)

$\begin{matrix} \sum_{i \in I} x_{ij}^{t} + \sum_{v = 1}^{t - 1} (\sum_{i \in I} x_{ij}^{v} - \sum_{k \in K} z_{jk}^{v}) ⩽ U_{j}, \\ t = 2, 3, \dots, T \end{matrix}$ (12) $x_{ij}^{t} \in N^{+} z_{jk}^{t} \in N^{+} y_{j}^{t} = {0, 1}$ (13) where constraint (8) states that there should be no materials flowing into the closed distribution centers. Constraint (9) states that there should be no materials flowing out from the closed distribution centers. Constraints (10) and (11) refer to the flow balance, while constraint (10) is the flow balance at distribution centers in the first period, and constraint (11) is the flow balance at distribution centers in every period except the first period. Constraint (12) ensures that the sum of materials inventory and input must not exceed the given capacity of distribution center in the each period, while there is inventory in distribution centers. Constraint (13) is the integrality constraint on the decision variables.

3 The proposed algorithm

The above emergency logistics network optimization model belongs to dynamic multi-objective optimization problems (DMOPs).Different from the static multi-objective optimization problems (MOPs), the variation of time-dependent parameters in the DMOPs will lead to the change of optimization environment. When the environment changes, the optimal solutions of the optimization problem will migrate, and the optimal solutions obtained in the current environment is no longer suitable for the new environment. Therefore, the dynamic optimization algorithm should be able to detect the environment changes quickly and find the optimal solutions in the new environment as soon as possible [11].

In the MOPs and DMOPs, the objectives often conflict with each other. It is usually impossible to find a solution that makes all objectives optimal simultaneously, so the concept of Pareto dominance is used to measure the merits and demerits of the solutions [12].

Differential evolution (DE) algorithm, as a meta-heuristic optimization algorithm, has been widely used in static single-objective and multi-objective optimization problems [13]. To solve the proposed dynamic emergency logistics network optimization model, an improved differential evolution algorithm named dynamic self-adaptive multi-objective differential evolutionary algorithm (DSMODEA) is proposed in this paper. On the one hand, the new algorithm can detect and respond to environmental changes in dynamic optimization. On the other hand, it improves the ability of global exploration and local exploitation of differential evolution algorithm when the environment changes. It also improves the convergence speed of the algorithm and prevents the algorithm from converging prematurely into local optimum.

Pseudo Code of the DSMODEA
1.	Initialization: Initialize the population randomly, let iter = 0, t = 1, initial environment detector equal to 0, and the fitness of elite individual is infinite.
2.	While iter < max _ iter
3.	Calculate the environment detector
4.	If environment changes (environment detector = 1)
5.	Implementing the environment change response strategy;
6.	end if
7.	Update the population according to the DE (optimization process):
8.	Step 1: Calculate all the fitness functions of individuals in the population, and put the non-dominance individuals in an archive according to the dominance relationship;
	calculate the crowding distances of the individuals in archive, and select the individuals with the smallest distances as the elite individual;
9.	Step 2: Mutation: the mutation individuals are generated from the parent individuals according to the mutation strategy and the elite individual;
10.	Step 3: Crossover: crossover the mutation individuals and the parent individuals with a given probability, then we get the crossover individuals.
11.	Step 4: Selection: mix the crossover individuals and the parent individuals, select non-dominance individuals as the offspring individuals.
12.	Update environment change counter t;
13.	iter = iter + 1;
14.	end while
15.	The individuals in archive are final optimal solutions.

3.1 Algorithm framework

The DSMODEA adopts an environment change detect operator and an environment change response strategy, which make the static differential evolutionary algorithm be able to solve the DMOPs. We also modified the evolutionary strategy of the traditional differential evolution algorithm to improve the performance of the algorithm. The pseudo code of the DSMODEA algorithm is as follows: In the rest of this paper, fitness function calculation, the environmental change detect operator, the environmental change response strategy and the differential evolution optimization process are introduced respectively.

3.2 Fitness function calculation

For the constrained optimization problem, on the one hand, the optimal solutions of the model should satisfy all constraints. On the other hand, it is necessary to use the information of some infeasible solutions. Because that the optimal solution may appear in the boundary of the feasible domain at the early time of the algorithm [14]. Therefore, a constraint handle technology should be carried out.

Constraint handle technology can be classified into three different groups according to the study in [15]: penalty functions methods, methods based on preference of feasible solutions over infeasible solutions, and methods based on multi-objective optimization techniques. It is usually hardly to make a rule to balance the objectives and constraint violation when use methods based on preference of feasible solutions over infeasible solutions. In addition, the multi-objective method transforms the constraints into objectives, which greatly increases the number of objectives and is difficult to be solved. Penalty functions methods, as the most common used method, has the characteristics of feasibility and easy to understand. The common penalty function method includes static penalty function method, dynamic penalty function method and adaptive penalty function method. In this paper, an adaptive penalty function method is adopted to transform constrained optimization problems into unconstrained optimization problems.

Define fitness functions as the sum of objective functions and penalty term: $f_{i} (\overset{⇀}{x}) = C_{i} (\overset{⇀}{x}) + M \cdot P (\overset{⇀}{x})$ (14) where $f_{i} (\overset{⇀}{x})$ is the ith fitness function, $C_{i} (\overset{⇀}{x})$ is the ith objective, M is penalty coefficient, and $P (\overset{⇀}{x})$ is constraint violation.

Define constraint violation $P (\overset{⇀}{x})$ as follows: $P (\overset{⇀}{x}) = \sum_{j = 1}^{m} c_{j} (\overset{⇀}{x})$ (15) where $c_{j} (\overset{⇀}{x}) = max (0, g_{j} (\overset{⇀}{x}))$ , while $g_{j} (\overset{⇀}{x}) ⩽ 0$ is inequality constraint. $c_{j} (\overset{⇀}{x}) = max (0, | h_{j} (\overset{⇀}{x}) - ɛ |)$ , while $h_{j} (\overset{⇀}{x}) = 0$ is equality constraint. ɛ is a very small real number.

In order to ensure the global exploration ability in the early iteration of the optimization process, M is defined as: $M = M_{0} \cdot \frac{{iter}^{t}}{{max}_{i} {ter}^{t}}$ (16) where M₀ is a constant large real number, iter^t is iterations of the algorithm in period t, $max_{i} {ter}^{t}$ is maximum iteration in period t. It can be seen that the penalty for the infeasible solution gradually strengthened with the iteration of evolutionary process.

3.3 Environment change detect operator

In the present literature, most of the changes in the environment were detected by judging the change of objective function or constraint conditions. The “detection individual” is selected to calculate the values of objectives or constraints, and if the values change, it represents the environment has changed. However, because of the particularity of some functions, the values of functions may not change even if the parameters change. Therefore, a new environment change detect operator independent of decision variables is designed in this paper: $\begin{matrix} {Detector}_{t} = sgn (| exp (C_{1}^{t} - C_{1}^{t - 1}) - 1 |) \cup \\ sgn (| exp (C_{2}^{t} - C_{2}^{t - 1}) - 1 |) \cup \\ \dots \cup sgn (| exp (C_{M}^{t} - C_{M}^{t - 1}) - 1 |) \end{matrix}$ (17) where $sgn (x) = \frac{2}{π} \int_{0}^{+ \infty} \frac{sin xt}{t} dt$ is sign function and follows: if x > 0, sgn (x) =1, and if x = 0, sgn (x) =0. It is obvious that: (1) As long as one parameter in the model changes, the value of the detector is 1; (2) Only if all the parameters in the model are not changed, the value of detector is 0. Therefore, environment change detect operator can detect the change of the environment as long as the environment changes and be independent of decision variables.

3.4 Environmental change response strategy

In dynamic evolutionary algorithms, there are three basic response strategies when optimization environment changes [16]. The first strategy is to initialize the population randomly and restart a new optimization after an environment change occurs. This kind of method decomposes the dynamic optimization problem into several static optimization problems and optimizes them one by one. In this way, the information of population in the last environment will be ignored, which is important for the convergence of the algorithm. In the second strategy, the population of last environment will not be initialized when the environment changes and be as the initial population of the new environment. However, the population of the previous environments has strong convergence, and it is easy to make the algorithm fall into local optimum when using the second strategy. The third strategy is to keep one part of individuals remain unchanged and initialize another part of individuals of the population randomly when the environment changes. Then, the mixed population will be as the initial population of the next environment. This strategy could preserve some information of the population as well as keep the diversity of populations. However, it’s difficult to determine which part of the individual remains, and which part of the individual should beinitialized.

The DSMODEA algorithm in this paper should adopt two response strategies, respectively. The second strategy is adopted in the population, that is, when the environment changes, the population of the previous environment will be continue as the initial population of new environment. Since that the current dominant relationship of the individuals in archive will be not applicable in the new environment. On the contrary, the first strategy is adopted in archive, that is, when the environment changes, the archive is initialized randomly.

3.5 Archive update

(1) Calculate the fitness functions of all individuals of the population, then seek out all the non-dominating individuals and put them into the archive.

When the next generation population is obtained, the new population is mixed with the archive. Then, seek out all the non-dominating individuals of the mixed population and put them into the archive.

(2) Calculate the crowded distances of the individuals in archive [17]. Sort all individuals of archive in ascending order according to each fitness function. Set the first and last individual’s crowded distance to be infinity. Then, the crowded distance of the ith individual is the sum of the deviations of all objective functions of the ith individual and i + 1th individual. The pseudo code is as follows:

(3) Select individual with least crowded distance in archive as elite individual.

Pseudo Code of Crowded Distance
1. Set crowded distances of all individuals to be zero: d (i) =0
2. For each objective function j = 1 : m
3. Sort the individuals based on objective m
4. Assign infinite distance boundary values for each individual:
$d (1) = + \infty, d (N) = + \infty$
5. For i = 2 : N - 1
6. $d (i) = d (i) + \frac{d (i + 1)_{m} - d (i - 1)_{m}}{f_{max}^{m} - f_{min}^{m}}$
7. End for
8. End for

3.6 Adaptive differential evolution algorithm

There are three main steps in the differential evolution: mutation, crossover, and selection.

Step 1. Mutation

This paper uses the DE/ current-to-best/1 strategy to generate mutation population: $v_{i} = x_{i} + F \times (x_{elite} - x_{i}) + F \times (x_{rand} - x_{i})$ (18) where F is mutation rate, the vector x_i = x_i1, x_i2, ⋯ , x_iD is the ith parent individual, v_i = v_i1, v_i2, ⋯ , v_iD is the ith mutation individual, x_elite is the elite individual in the archive, and x_rand is the random individual selected from the archive randomly, x_i ≠ x_elite ≠ x_rand. DE/current-to-best/1 can make full use of the information of the elite individual in the current population so as to maintain a good convergence. In the later iteration of the algorithm, a large mutation rate is not conducive to the convergence of the algorithm. So an adaptive mutation rate is designed to control the mutation rate of the algorithm: $F = F_{0} - \frac{iter - ({max}_{i} {ter}^{1} + \dots + {max}_{i} {ter}^{t - 1})}{{max}_{i} {ter}^{t}}$ (19) where F₀ ∈ [0, 1] is initial mutation rate, $max_{i} {ter}^{t}$ is the maximum iterations in the tth environment, iter^t is the current iteration in the tth environment, iter is the current iteration in the entire algorithm.

It can be seen that the adaptive strategy can not only guarantee that the mutation rate of the algorithm decreases with the iteration in every environment, but also can initialize the mutation rate when the environment changes, and ensure the diversity of the population in the new environment.

Step 2. Crossover

The crossover individuals u_i are generated by the parent individual x_i and the mutation individual v_i as the follow equation: $u_{ij} = {\begin{matrix} v_{ij}, & if ({rand}_{ij} ⩽ CR) \\ x_{ij}, & otherwise \end{matrix}$ (20) where rand_ij is a real number following uniform distribution between [0, 1], and CR ∈ [0, 1] is crossover rate. The jth variable u_ij in the crossover individual is equal to the corresponding variable v_ij in the mutation individual when rand_ij is smaller than CR, otherwise u_ij will be equal to x_ij of the parent individual.

Step 3. Selection

The parent individuals and the crossover individuals are mixed, and the fitness functions of the mixed population are calculated. Pick out the non-dominance individuals as the offspring individuals for the next iteration. If the amount of non-dominance individuals is less than the size of population, then select the rest individuals from the mixed population randomly. If the amount of non-dominance individuals exceeds the size of population, then remove the redundant non-dominance individuals randomly.

4 Case study

4.1 Case description

In an emergency relief task, there are two supply centers and four distribution centers to cover the material demands of six shelters. The emergency supply task is completed in five periods. The material demands, shortage costs and the transportation costs of each period are different. The related data is given by Tables 1 4. The simulation studies are carried out in a MATLAB 2014b platform on an ASUS laptop with 5-6300HQ 2.3GBz CPU, 4GB RAM in Windows 7.0(64-bit) environment. The parameters are as follow: population size N = 100, the maximum iteration $max_{i} ter = 500$ , period = 5, environment change count $t = ⌈ \frac{iter}{{max}_{i} ter / period} ⌉$ , t₀ = 1, initial environment change detector Detector₁ = 0, initial penalty coefficient of constraint violation M₀ = 10000, initial mutation rate F₀ = 0.9, and crossover rate CR = 0.8.

Table 1
Transportation cost of unit materials between nodes (hundred yuan)

Distribution centers 1 Distribution centers 2 Distribution centers 3 Distribution centers 4

Supply centers1 20/25/30/40/50 30/35/40/45/50 26/32/40/45/35 10/12/16/20/16

Supply centers2 30/35/45/25/30 12/15/16/20/22 35/38/40/45/55 28/36/40/45/50

Shelters1 2/2.5/3/3.2/3.8 1/1.5/1.6/1.9/2.2 1.5/1.8/2/2.6/3 2/2.2/2.6/3/3.2

Shelters2 1.5/2/2.5/3/3.2 1.6/1.8/1.8/2.3/3 0.5/0.8/1/1/1.2 1/1.2/1.5/1.6/1.8

Shelters3 0.6/1/1.2/1.6/2 3/3.2/3.5/3.8/4 1.2/1.5/2/2.2/3 1.5/1.6/2/2.5/3

Shelters4 2.6/2.8/3/3.6/4 3/3.6/4/4.8/5 2/2.6/3/3.2/3.5 0.6/1/1.2/1.5/1.8

Shelters5 2.2/2.6/3/3.2/3 1/2/2.5/3/3.5 1.2/1.6/1.8/2/2 1.6/2/2.6/3/3.5

Shelters6 2/2.3/2.6/3/3.2 3/3.2/3.5/3.6/4 1/1.2/1.2/2/5 1/1.5/1.6/2/2.2

Each grid of data are transportation costs of period 1 to period 5 from left to right respectively

	Distribution centers 1	Distribution centers 2	Distribution centers 3	Distribution centers 4
Supply centers1	20/25/30/40/50	30/35/40/45/50	26/32/40/45/35	10/12/16/20/16
Supply centers2	30/35/45/25/30	12/15/16/20/22	35/38/40/45/55	28/36/40/45/50
Shelters1	2/2.5/3/3.2/3.8	1/1.5/1.6/1.9/2.2	1.5/1.8/2/2.6/3	2/2.2/2.6/3/3.2
Shelters2	1.5/2/2.5/3/3.2	1.6/1.8/1.8/2.3/3	0.5/0.8/1/1/1.2	1/1.2/1.5/1.6/1.8
Shelters3	0.6/1/1.2/1.6/2	3/3.2/3.5/3.8/4	1.2/1.5/2/2.2/3	1.5/1.6/2/2.5/3
Shelters4	2.6/2.8/3/3.6/4	3/3.6/4/4.8/5	2/2.6/3/3.2/3.5	0.6/1/1.2/1.5/1.8
Shelters5	2.2/2.6/3/3.2/3	1/2/2.5/3/3.5	1.2/1.6/1.8/2/2	1.6/2/2.6/3/3.5
Shelters6	2/2.3/2.6/3/3.2	3/3.2/3.5/3.6/4	1/1.2/1.2/2/5	1/1.5/1.6/2/2.2
Each grid of data are transportation costs of period 1 to period 5 from left to right respectively

Table 2

Related information of distribution centers

	Distribution centers 1	Distribution centers 2	Distribution centers 3	Distribution centers 4
Opening cost (hundred yuan)	100	150	80	180
Inventory cost (hundred yuan)	2	6	8	5
Maximum quantity	25	30	20	36

Table 3

Shortage costs of unit materials at shelters (hundred yuan)

	Shelter	Shelter	Shelter	Shelter	Shelter	Shelter
	1	2	3	4	5	6
Period 1	10	12	8	15	15	11
Period 2	15	18	12	16	18	13
Period 3	20	22	15	18	23	15
Period 4	25	25	20	25	25	15
Period 5	30	28	22	28	28	18

Table 4

Materials demand of shelters in each period

	Shelter	Shelter	Shelter	Shelter	Shelter	Shelter
	1	2	3	4	5	6
Period 1	12	13	15	10	8	5
Period 2	13	14	16	11	10	5
Period 3	16	15	17	10	10	4
Period 4	16	13	15	11	9	6
Period 5	15	15	17	12	10	6

4.2 Computational results

The DSMODEA algorithm is used to solve the model, and the non-dominated solutions of every period are obtained, which is shown in Fig. 2. Figure 2 shows the distribution of non-dominated solutions in the objectives space. The horizontal axis represents the value of the first fitness function, and the longitudinal axis represents the value of the second fitness function. It can be see that, five non-dominated solutions are obtained in the first period, one non-dominated solution is obtained in the second period, three non-dominated solutions are obtained in the third period, one non-dominated solution is obtained in the fourth period, and one non-dominated solution is obtained in the fifth period.

Fig.2

Distribution of non-dominated solutions.

Every non-dominated solution of the emergency logistics network optimization model represents an optimal feasible project of the emergency material supply. For a more accurate quantitative analysis of the results, Table 7 shows the value of corresponding fitness functions and objectives of all the non-dominated solutions. The first fitness function is the sum of total supply costs (Equation 1) and the constraint violation penalty term. The second fitness function is the sum of inverse of material fill rate (Equation 7) and the constraint violation penalty term. We can see that, the first fitness functions are equal to the supply costs, and the second fitness functions are equal to the invers of fill rate. That to say, the penalty functions of constraint violation are all zero for every non-dominated solutions, that is, all non-dominated solutions are feasible. The optimal solutions obtained by DSMODEA algorithm satisfy all the constraints. Analyze the fill rate of every solution in Table 5, we can find that all the demand fill rates are all close to 1, that to say, the obtained supply programs can not only meet the demand of materials but also achieve good economic benefits without causing too much waste.

Table 5

Fitness functions and objectives of optimal programs

	Fitness function 1	Fitness function 2	Supply costs	Fill rate
Period 1	2109.7	1.2917	2109.7	0.7742
	2189.5	1.2500	2189.5	0.800
	2073.7	1.3333	2073.7	0.7500
	2081.2	1.3125	2081.2	0.7619
	2057.2	1.3542	2057.2	0.7385
Period 2	2277.1	1.0000	2277.1	1.0000
Period 3	3043.8	0.9913	3043.8	1.0088
	3069.8	0.9635	3069.8	1.0378
	3266.2	0.9497	3266.2	1.0530
Period 4	2687.7	1.0000	2687.7	1.0000
Period 5	2835.2	1.0000	2835.2	1.0000

4.3 Sensitivity analysis

4.3.1 Results with different F

To validate the effectiveness of the adaptive mutation strategy used in DSMODEA algorithm, we set the mutation rate as fixed value from F = 0.1 to F = 0.9 by 0.2 intervals. Then use these algorithms with fixed mutation rates to solve the same case, and compare their results with our adaptive strategy. Figure 3 shows the distribution of results in objective space at each period. Figure (a) corresponds to the adaptive strategy, and the rest of the figures correspond to the other mutation rates in turn. First of all, when the mutation rate is equal to 0.1 and 0.3, the values of fitness functions are very large. In these cases, the fitness functions include the objective function and the constraint violation, that is, the obtained solutions are infeasible. This is due to the fact that when the mutation rate is too small, the diversity of the population is poor, and the algorithm falls into local optimum. Secondly, it can be seen that, when F is equal to 0.5, 0.7 and 0.9, the solutions are all feasible. Meanwhile, with the increase of F, the corresponding fitness functions increase gradually. This is because the increase of mutation rate leads to the diversity of population increases, which is conducive to finding the global solution, but not conducive to the convergence of the algorithm. Thirdly, it can be found that the results obtained by using adaptive strategy are the optimal with the minimum fitness functions. This is because the mutation rate can be adjusted according to the state of the population, which can balance the ability of global exploration and local exploitation of the algorithm. To sum up, the value of mutation rate has a great influence on the performance of the algorithm. When F is too small, the algorithm is easy to fall into local optimum, and when F is too large, it is not conducive to the converge of the algorithm. The adaptive strategy in DSMODEA of this paper solves the problem very well.

Fig.3

Distribution of non-dominated solutions with different mutation rates.

4.3.2 Result with different algorithms

In order to further verify the performance of the algorithm, DNSGA-II and DMOEAD-M [18] are set up as comparison algorithms. The population size, iteration number, cross rate, and penalty coefficient of the comparison algorithms are consistent with those of the proposed algorithm. The other parameters are consistent with those in the original literature, which are given as follows. In the DNSGA-II algorithm, the multi-sample introduction ratio is set as 0.2, the mutation rate is set as 1/dim, where dim represents the dimension of the individual. In the DMOEAD-M, the size of memory pool is set to 100, the number of weight vectors in the neighborhood is set to 20, and the mutation rate is set as 1/dim.

Each algorithm is used to optimize the case 15 times, and then the non-dominant solutions are obtained according to the Pareto dominance relation. The distribution of the non-dominant solutions in the objective space is shown in Fig. 4. The black marks in the figure represent the results obtained by the proposed DSMODEA. The blue marks represent the results of DMOEAD-M algorithm, and red marks represent the results of DNSGA-II algorithm. It can be seen that, in the first/second/third/fourth period, the results of the algorithm obtained by DSMODEA completely dominate the results of the other two algorithms. The dominant relationship of the results in the fifth period is DMOEAD - M ≺ DSMODEA ≺ DNSGA–II Therefore, it can be concluded that the results obtained in this paper are better than those obtained by the other two algorithms.

Fig.4

Distribution of non-dominated solutions with different algorithms.

The average values of the fitness function obtained from the 15 times simulations are recorded in Table 6. The total cost of the whole period is the sum of the cost of each period, and the total fill rate of is the product of fill rate of each period. The smaller fitness function corresponds to the better result, and the optimal values of the fitness functions are represented by bold numbers. For the first fitness function, the optimal values in the first/second/third/fourth period, as well as the whole period, are obtained by the proposed DSMODEA algorithm. For the second fitness function, the DSMODEA algorithm obtains the optimal value only in the first and the third period, however, the fitness function value of the whole period is much smaller than those of the other two algorithms. In summary, the results of DSMODEA dominate the results of the other algorithms.

Table 6

Fitness functions with different algorithms

	DSMODEA		DMOEAD-M		DNSGA-II
	Fitness function 1	Fitness function 2	Fitness function 1	Fitness function 2	Fitness function 1	Fitness function 2
Period 1	2935.3	1.4292	3273.6	1.6667	3496.2	2.0833
Period 2	3267.6	1.2058	3636.3	1.1159	3814.6	1.3623
Period 3	3331.8	0.8660	3470.4	1.0245	4247.4	1.1883
Period 4	3868.4	1.1393	3959.3	1.1000	4363.4	1.1428
Period 5	4524.4	1.1233	4445.4	1.1066	4867.1	1.1466
Total	17927.5	1.9097	18785	2.3068	20788.7	4.4190

5 Conclusion

In this paper, the problems of real world emergency logistics network optimization are considered. The location-allocation joint optimization problem is formulated by a multi-period and multi-objective optimization model. To solve the dynamic multi-objective problem, a dynamic self-adaptive multi-objective differential evolutionary algorithm is proposed. An environment change detect operator and an environment change response strategy are used to make the classical static differential evolutionary algorithm to be able to solve the dynamic optimization. In addition, a self-adaptive strategy is used to improve the exploration and exploitation of the algorithm when optimization environment changes. The numerical examples show that, the proposed DSMODEA algorithm can solve the emergency logistics network optimization very well. The performance improvement of the modified algorithm has been proved by the case studies, too. The basic framework of emergency logistics network optimization and the meta-heuristic intelligent optimization algorithm of this paper have provided some ideas and references to the emergency materials supply in real world.

In this paper we don’t consider the problem that the parameters such as demands and cost vary over time continuously. There are also some improvements should be made in the dynamic optimization algorithm, such as prediction strategy, multi-population strategy, and so on. These improvements will be studied in our future research.

References

Caunhye

A.M.

, Nie

, Pokharel

, Optimization models in emergency logistics: A literature review, Socio-Economic Planning Sciences 46(1) (2012), 4–13.

, Zhao

, Zhu

, Wyatt

, Covering models and optimization techniques for emergency response facility location and planning: A review, Mathematical Methods of Operations Research 74(3) (2011), 281–310.

Han

, Guan

, Sh

, Optimization based method for supply location selection and routing in large-scale emergency material delivery, IEEE Transactions on Automation Science and Engineering 8(4) (2011), 683–693.

, Xu

, Li

, A modified particle swarm optimization algorithm for optimal allocation of earthquake emergency shelters, International Journal of Geographical Information Science 26(9) (2012), 1643–1666.

Mohammadi

, Ghomi

S.M.

, Jolai

, Prepositioning emergency earthquake response supplies: A new multi-objective particle swarm optimization algorithm, Applied Mathematical Modelling 40(9-10) (2016), 5183–5199.

Zhao

, Chen

, Risk-based optimization of emergency rescue facilities locations for large-scale environmental accidents to improve urban public safety, Natural Hazards 75(1) (2014), 163–189.

Wang

, Choi

T.M.

, Liu

, Yue

, A novel hybrid ant colony optimization algorithm for emergency transportation problems during post-disaster scenarios, IEEE Transactions on Systems, Man, and Cybernetics: Systems 48(4) (2018), 545–556.

Dan

, Zhu

, Li

, Sang

, Liu

, Dynamic optimization model and algorithm design for emergency materials dispatch, Mathematical Problems in Engineering (2013), 1–6.

Shen

, Nie

, Zhang

H.M.

, Dynamic network simplex method for designing emergency evacuation plans, Transportation Research Record: Journal of the Transportation Research Board 22(1) (2007), 83–93.

10.

Chang

F.S.

, Wu

J.S.

, Lee

C.N.

, Shen

H.C.

, Greedy-search-based multi-objective genetic algorithm for emergency logistics scheduling, Expert Systems with Applications 41(6) (2014), 2947–2956.

11.

, Yu

, Lu

, Yen

G.G.

, Cai

, Differential evolution mutation operators for constrained multi-objective optimization, Applied Soft Computing 67 (2018), 452–466.

12.

Tusar

, Filipic

, Visualization of pareto front approximations in evolutionary multiobjective optimization: A critical review and the prosection method, IEEE Transactions on Evolutionary Computation 19(2) (2015), 225–245.

13.

Baraldi

, Bonfanti

, Zio

, Differential evolution-based multi-objective optimization for the definition of a health indicator for fault diagnostics and prognostics, Mechanical Systems and Signal Processing 102 (2018), 382–400.

14.

Coello Coello

C.A.

, Constraint-handling techniques used with evolutionary algorithms, Proceedings of the companion publication of the 2015 on Genetic and Evolutionary Computation Conference - GECCO Companion ’15, 2015.

15.

Panda

, Pani

, A Symbiotic Organisms Search algorithm with adaptive penalty function to solve multi-objective constrained optimization problems, Applied Soft Computing 46 (2016), 344–360.

16.

, Li

, Multi-strategy ensemble particle swarm optimization for dynamic optimization, Information Sciences 178(15) (2008), 3096–3109.

17.

Deb

, Pratap

, Agarwal

, Meyarivan

, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation 6(2) (2002), 182–197.

18.

Liu

, Zeng

H.W.

, Memory enhanced dynamic multi-objective evolutionary algorithm based on decomposition, Journal of Software(7) (2013), 1571–1588.