Robust optimization of the multi-objective multi-period location-routing problem for epidemic logistics system with uncertain demand

The effective distribution of relief to an emergency logistics system plays a crucial role during the disaster response phase. Considering stochastic characteristics of relief demand, this study investigates the robust optimization of a multi-objective multi-period location-routing problem for epidemic logistics, a special emergency logistics, with uncertain scenarios. A corresponding robust multi-objective multi-period optimization model is proposed, which aims to determine the optimal location of temporary relief distribution centres and route planning simultaneously. The optimization objectives include the total travel time, the total cost, and the disutility of relief service. To solve the above optimization model, a preference-inspired co-evolutionary algorithm with Tchebycheff decomposition (PICEA-g-td) is given. The performance of the proposed PICEA-g-td is evaluated by comparing it with NSGA-II, MOEA/D and PICEA-g. The experimental results show that the proposed algorithm performs better than the other three algorithms in terms of the solution quality. Finally, some useful management insights are obtained.


I. INTRODUCTION
Outbursts of emergencies from public health events have occurred worldwide, causing a large number of deaths. Examples in recent years include the 2003 Severe Acute Respiratory Syndrome (SARS) outbreak in Canada, China, Hong Kong, Singapore, Taiwan, and Vietnam [1], the 2009 H1N1 influenza in Turkey, Canada, Florida, and China [2], the outbreak of Ebola in Uganda, Sudan, Congo and some countries in West Africa [3,4], and the Coronavirus Disease 2019 (COVID-19) Outbreak, which is currently occurring throughout the world [5]. These epidemic diseases not only have caused casualties but also had a serious impact on the national and regional economies.
Since December 2019, more than 80 thousand cases of pneumonia caused by Coronavirus Disease 2019 (COVID -19) have been reported in China. Wuhan, as a severe epidemic area in China, has treated 50 thousand people, and there have been more than 50 thousand infected people [6]. A series of measures have been taken to control the spread of the epidemic and to speed up epidemic disease recovery. From January 23rd, 2020, an effective lockdown was used in the urban area of Wuhan [7]. In other words, traffic in Wuhan was restricted to make people and vehicles unable to communicate to the outside world without permission. By May 2020, the number of infected people in China dropped below 1,000, and the death rate was under 5% (approximately 4650 deaths). In contrast, the number of COVID-19 cases globally remains at the highest levels since the beginning of the pandemic, with over 5.7 million new weekly cases, and new deaths continuing to increase for the seventh consecutive week, with over 93,000 deaths by May 2021 [8]. Compared with the current situation outside China, Wuhan epidemic prevention measures can be used as a good case study. To ensure that a large number of the epidemic infected people receive timely medical treatment, epidemic logistics, as a special type of emergency logistics, decision-makers must make optimal decisions to ensure a stable supply of medical materials during a city lockdown. Disaster temporary relief distribution centres play an important role in emergency supplies, i.e., three logistics parks were chosen as disaster temporary relief centres to provide logistics services during the lockdown of Wuhan. At the same time, sixteen module hospitals that are temporary hospitals converted from sports stadiums and exhibition centres were set up to treat the rapidly growing number of infected epidemic people. Locating disaster temporary relief centres and planning rational routes to provide relief to these hospitals have a significantly positive impact on the overall performance of epidemic logistics.
Infection materials, such as medical materials and daily supplies during the epidemic outbreak, were defined as class 6 hazardous materials by the Jefferson Lab [9]. For the transportation of hazardous materials, researchers suggest avoiding societal risk and reducing the potential exposure of the population in the logistics [10,11]. Compared with traditional logistics, epidemic logistics must respond quickly within a short time of the outbreak to effectively control an epidemic outbreak [12]. There are many challenges to the location of temporary relief distribution centres and route planning, while the external environment changes rapidly with the disaster outbreak. The location-routing problem (LRP), as the core of emergency logistics, has wide application value in socioeconomic activities [13,14]. In the traditional location-routing problem, decision makers access information about customers and make two types of decisions: the location of distribution centres and the design of the vehicle routes [15]. Considering that information about the relief needs is imprecise as the epidemic disaster progresses, it is one of the most important decisions for emergency logistics management to dynamically optimize temporary relief distribution centres and route planning [16]. Normally, in stochastic logistics, there are some uncertain parameters within the planning horizon. To reduce the error caused by parameter variations, stochastic problems can be divided into multiple periods and adopt multiple decision making to reduce decision errors [17,18]. Considering that multi-period LRPs are a much better model to address stochastic location and routing problems with uncertain parameters, a multi-period LRP model is used to solve the dynamic LRP under emergency logistics [19]. The robust predictive control approach is a useful tool to handle emergency logistics with uncertain data such as supply and demand [20]. Based on the robust optimization approach of Bertsimas and Sim, robust relief route planning is used to ensure that relief supplies are met as much as possible [21].
Different from traditional logistics, which mainly pursue the maximization of benefits, emergency logistics decision makers prefer to consider the relief process as a humanitarian issue [16]. Timely relief supplies can ensure that epidemic disaster-infected patients are treated as quickly as possible, and thus, the backing time is a key factor for emergency logistics. Humanitarian relief operations are concerned by a variety of stakeholders, such as local governments, the military and non-governmental organizations. Therefore, a major challenge of emergency logistics also focuses on providing equitable services to all aid recipients. A good example is to encourage a vehicle not to necessarily satisfy a vertex's entire demand but rather to save supply to serve another vertex when relief supplies are scarce [22]. In addition, non-profit indicators such as fair distribution (such as minimizing the absolute deviations of a fraction of unsatisfied demands between affected areas to fairly allocate recourse) reflect the humanistic concern of relief operations [23,24]. To balance the economy, society and efficiency of epidemic logistics, the problem features three objectives: the total travel time, the total cost, and the fairness of relief allocation. This study considers a multi-objective multi-period robust location-routing problem with uncertain demand (MMRLRP) to optimize epidemic logistics during public health events. In a planning horizon, the demand for relief is stochastic in the epidemic area. However, the information about the demand points is known at the beginning of each period. Several potential relief distribution centres are placed in reasonable locations, and the challenges are to make decisions on the relief distribution centre locations and on allocating demand points under several periods of the epidemic. The problem consists of locating several relief distribution centres to serve in each period and planning routes for the demand points under the condition that the following requirements are met: (1) the demand allows for a shortage of relief supplies (2) each demand point is visited once, (3) there are shortage of relief supplies at the early and mid-stage of an outbreak, and (4) each vehicle route has a carrying capacity Q.
The main contributions of this paper are shown as follows: First, a robust optimization model on a multi-objective location-routing problem for epidemic logistics system design with multi-period is proposed. Second, an improved heuristic algorithm named PICEA-g-td is developed and examined by numerical experiments, as well as a real-world case study. The proposed algorithm is embedded with a decomposition strategy, so as to improve its computational performance. Finally, the sensitivity analysis of some key parameters on the proposed model is conducted, and some useful management insights are obtained by the numerical experiments and case study on the epidemic disease outbreak of Coronavirus disease 2019  in Wuhan.
The remainder of this paper is organized as follows. Section 2 reviews the related work. Section 3 elaborates the problem description and describes the mathematical model. The details of the improved preference-inspired co-evolutionary algorithm are described in Section 4. Numerical experiments are conducted in Section 5, followed by conclusions in Section 6.

II. LITERATURE REVIEW
The most important function of emergency logistics is to ensure the supply of relief--i.e., medical mask, protection suit, and medicine--quickly and in sufficient quantities during the disaster outbreak. Therefore, epidemic logistics ensures that patients are treated and prevents the spread of the disease among affected people. For the past 20 years, emergency logistics have received more attention because of their theoretical and practical significance. There are several typical optimization models for emergency logistics: multi-objective relief distribution models and stochastic relief distribution models.
Considering the characteristics of relief supply, emergency logistics have usually been constructed as multi-objective location-routing models with timeliness indicators [25,26]. Wang et al. [27] proposed a multi-objective location-routing model for the relief distribution problem. The model considers three objectives: the total cost, the travel time, and the reliability under high uncertainty in emergency logistics. Gan and Liu [28] developed the modified non-dominated sorting genetic algorithm II to optimize the emergency logistics model in large-scale disaster relief, which minimizes the transportation cost and the total unsatisfactory time. However, the information from the model is considered to be deterministic, and emergency logistics scheduling should be regulated with uncertain and updated information. Feng et al. [29] considered the timeliness and economy of emergency logistics and proposed a location selection of emergency supplies and route planning to minimize the total transport length and cost. The model converts the two objectives mentioned above into a single objective, therefore, it is impossible to obtain Pareto solutions that can provide sufficient decision information. Vahdani et al. [19,30] suggested a multi-objective, multi-period, multi-commodity location-routing model to optimize the transportation time, the total cost, and the routing reliability. Two meta-heuristic algorithms were developed to seek the Pareto solutions and for locating relief centres, allocating demand points and planning vehicle routes. Beyond the above-mentioned multi-objective emergency logistics problems that are related to the objective function of timeliness, there are still some studies about other humanitarian factors. Bozorgi-Amiri et al. [31] considered maximizing the affected areas' satisfaction levels to be a key factor in emergency logistics. They constructed a multi-objective robust model that involved minimizing shortages in the affected areas and minimizing expected total costs. Mohamadi et al. [32] presented a multi-objective emergency logistics model that consists of three objectives: the total transportation distance, the service coverage of the facilities and the routes' availabilities. Barzinpour considered both economic and humanitarian objectives in urban disasters and proposed a multi-objective location distribution model to solve the location-allocation problem of urban relief distribution [33]. In addition to the cost, the model considers cumulative coverage of the population to be a key factor in humanitarian relief chain management. According to the papers mentioned above, multi-objective emergency logistics mainly focus on the economy, timeliness and humanitarian objectives to address the rescue process during the disaster outbreak.
Different from traditional business logistics with a stable market environment and deterministic information, decision-makers of emergency logistics must make decisions in stochastic environments [34,35]. Stochastic LRPs with random parameters are more suitable for solving real-life location problems with routing and usually divide the planning horizon into multiple periods [17]. Bozorgi-Amiri et al. [16] considered a humanitarian emergency logistics problem as a multi-objective stochastic model, and they converted the multi-objective model to a single-objective model by the ε -constraint method. Duhamel et al. [36] considered the population distribution for post-disaster situations and proposed a multi-period location-allocation model to seek optimal combinations to maximize population assistance. Moreno et al. [37] studied a multi-period location-transportation problem for delivering relief supplies under uncertain conditions in emergency logistics. Two stochastic programming models and decomposition heuristics are proposed for relief facility location and transportation decisions. Vahdani et al. [19] considered the repair of damaged roads after a disaster and suggested a multi-period location-routing model for timely relief distribution to locate the distribution centres and arrange vehicle routing. In this research, the supplies of relief are divided into several periods because the repair of roads lasts for several periods depending on the proportion of roads damaged. Yu et al. [38] proposed a multi-period reverse logistics network focused on the optimal location of facilities and vehicle routes to provide a strong response to healthcare services, to tackle medical waste during an epidemic outbreak. A real-world case study based on the outbreak of COVID-19 in Wuhan was used to illustrate the applicability of the model. In addition to the relief distribution, multi-period models are widely used in business logistics with uncertainty. Klibi et al. [39] studied stochastic multi-period location-transportation problem consisting of a two-stage stochastic programming: the strategic level is responsible for the network location and allocation decisions, and the user level is responsible for vehicle transportation decisions. Imen et al. [40] considered that the planning horizon must be partitioned into several periods to illustrate the uncertainty of the demand and cost. They proposed a two-echelon stochastic multi-period capacitated location-routing model to design a logistics network to satisfy future flexibility distribution for a company. Rabbani et al. [41] innovated a multi-period stochastic location-routing model for industrial hazardous waste management and developed the non-dominated sorting genetic algorithm-II (NSGA-II) based on a Monte Carlo simulation to tackle location, inventory and routing decisions. Rafie-Majd et al. [42] suggested a perishable product supply chain with stochastic demand and proposed a multi-period inventory-location-routing model in which the time horizon of deliveries was divided into several time periods. To solve the model, a heuristic algorithm and Lagrangian relaxation method were used to obtain the optimal solution. Peiman Ghasemi et al. [43] proposed an uncertain multi-objective multi-commodity multi-period multi-vehicle location-allocation mixed-integer mathematical programming model for earthquake rescue. The research considered two objectives: the total cost and the amount of the shortage of relief supplies under an uncertain disaster scenario. In most existing studies on emergency logistics, uncertainty is modeled using a scenario-based probability approach. However, it is difficult to accurately set scenario-based probabilities during the response phase of a disaster, which is a tough challenge for decision-makers.
According to this literature review, existing studies about emergency logistics focus on multi-objective models or stochastic multi-period models, and most studies focus on post-disaster relief distribution [19,20]. Although some studies have addressed multi-objective multi-period relief distribution models with uncertainties, the uncertainty in the demand is mainly to research the demand level, which is usually considered to be distributed within a range [16,37]. This paper studies the multi-objective multi-period location-routing model considering the uncertainty of demand points for emergency supplies at various periods of a disaster. In particular, this study studies relief distribution in the emergence, development and control of epidemic diseases.
As discussed above, the location-routing problem is one stream of epidemic logistics and is a typical NP-hard problem solved by an evolutionary strategy algorithm. The multi-objective multi-period robust location-routing problem with uncertain demand (MMRLRP), as a variant of the standard LRP, poses great challenges to the current research in multi-objective optimization. To obtain Pareto solutions for multi-objective problems (MOPs), a variety of MOEAs have been proposed for MOPs during the past few decades [36]. Some classical MOEAs, such as NSGA-II [44], MOPSO [45,46] and MOSA [47], have proven to be advantageous in obtaining Pareto solutions. When the MOPs have more objectives, it can easily exist that the current population becomes non-dominated with each other, which makes the algorithms have poor performance [48]. Many-objective MOEAs (e.g., MOEA/D [49], HypE [50], MSOPS [51]) have been proposed to solve the problem with more objectives. Purshouse et al.
[52] proposed a decision-maker preference co-evolving concept to solve multi-objective problems. In this concept, a family of decision-maker preferences is co-evolved with candidate solutions to converge towards the Pareto front. Wang et al. [53,54] developed a preference-inspired co-evolutionary algorithm (PICEA-g) based on the co-evolved concept and verified the superiority of the proposed algorithm by comparison with NSGA-II, ε-MOEA, HypE and MOEA/D on multi-objective benchmark instances.

III. MATHEMATICAL MODEL
In this study, a multi-objective multi-period robust location-routing problem(MMRLRP) is considered to deliver relief to demand points during an epidemic disaster outbreak. The decisions of epidemic logistics included the location of temporary relief distribution centres (TRDCs), relief allocation of hospitals and vehicle routing in each period of the epidemic disaster outbreak. Medical masks, protection suits, and medicine are regarded as crucial relief items required during epidemic disaster outbreaks, and the demand is stochastic in each period. Servicing to all demand points of the route, vehicles will return to the TRDCs because of the need to prevent the spread of the epidemic, which is consistent with the practical operation.

A. PROBLEM DESCRIPTION
In the MMRLRP, information on the demand points, including the positions, service time, is known at the beginning of the planning horizon. However, factual demands are revealed only after the location-routing decision has been made. The stochastic demand takes values in , and represent the demand nominal value of demand point at period and the maximum deviations from this nominal value, respectively. Relief supplies are related to the treatment of patients, and a shortage of relief supplies is allowed in this paper.
The problem aims to determine the subset of TRDCs to be used and to plan routes for the demand points considering the condition of meeting the vehicle capacity and the capacity of TRDCs. The location of the TRDCs and the allocation of the demand points are adjusted based on the stochastic demand level. Considering the stochastic demand, this paper proposes a robust optimization approach to ensure that relief supplies (as many as possible) meet the capacity of TRDCs. In addition, it is noteworthy that the multi-period problem must account for the condition in the previous period at each period. There is a fundamental difference in the decision-making about each period of the location-routing problem separately.
The following three objectives are considered for the MMRLRP: (1) minimization of the total travel time of vehicles; (2) minimization of the total costs, including the rental costs of TRDCs, the fixed vehicle costs and the transportation costs; and (3) minimization of the disutility of relief service.
Objective 1 is for the pursuit of effectiveness. An early vehicle backing time means that faster assistance is provided for the affected population. Objective 2 is for the economic value. Even in the context of epidemic logistics, the MMRLRP should consider having limited funds. Objective 3 is for the fairness of emergency logistics. In addition, non-profit indicators such as fair distribution (such as minimizing the absolute deviations of a fraction of unsatisfied demands between affected areas to fairly allocate recourse) reflect the humanistic concern of relief operations [23,24]. Four objectives of the MMRLRP are as follows:

Sets and Indices
(1) The first objective: The sooner the vehicle reaches the demand points, the better we provide treatment to the affected population. The response time is considered to be key factor in relief logistics. We focus on the total traveling time of vehicles to ensure time effectiveness.
(2) The second objective is to minimize the total cost and to balance the rental costs of TRDCs and transportation costs in the planning horizon: (2) is the rental cost of the TRDCs. Part two is used to calculate the fixed cost of the vehicles. Part three is the transportation cost. The last part is the penalty cost for shortages of relief supplies.
The third objective is the disutility of relief service. This concept applies to measure the fairness of customers' access to emergency relief [22].
Constraints (5) and (6) ensure that only the selected TRDC can allocate vehicles and provide services for demand points. Constraint (7) ensures that the amount of relief shipped out shall not exceed the supply of the TRDC. Constraint (8) ensures that vehicle k starts from TRDC and drives through demand point only when demand point is assigned to TRDC in the period . Constraints (9) and (10) ensure that each route k serviced by vehicle k starts at one TRDC. Constraint (11) ensures that each demand point can only be served by one route in period . Constraint (12) ensures that each vehicle returns to the departure TRDC in each period. Constraint (13) ensures that each route is continuous. Formula (14) calculates the number of emergency supplies out of stock at each demand point in each period. Constraint (15) ensures that the total demand assigned to each vehicle does not exceed the loading capacity of the vehicle. Constraints (16) to (21) explain the start time, delivering services, and return journey in each period, using vehicle. Constraints (22)- (27) ensure that the decision variables are binary integers and have non-negative values.

C. ROBUST COUNTERPART OF THE MMRLRP
Considering the uncertain demand of relief in the MMRLRP, a robust optimization method is used to describe the uncertainty of the model. Before deducing the robust counterpart of the MMRLRP, the robust optimization method proposed by Bertsimas and Sim [21] and Najafi et al. [55] is introduced as follows: In the constraint of the optimization problem, some elements of the coefficient are uncertain and take the value of , where b is is the nominal value of b is and b is is the maximum deviation from . Let be the total number of uncertain parameters available in constraint , and let be the uncertainty robustness budget, which takes a value in the interval . Specifically, coefficients vary in the constraint , and one coefficient varies in the range of it i . So, the right-hand side of the inequation can be rewritten as: To ensure that the inequality holds, a protective function β i is proposed as follows: Therefore, the optimization model (28) can be rewritten as Finally, the robust counterpart of the optimization problem (28) is written as follows: Since Eq. (14) involves uncertain parameters, the MMRLRP deduces its robust counterpart using the robust optimization method proposed by Bertsimas and Sim [21].
In constraint (14), the set of uncertain coefficients is represented as , which represents the number of uncertain demand points in route k. The uncertain parameter is proposed to control the robustness of the solution. Based on the constraint (14), the protective function (15) can be rewritten as According to robust optimization theory, the robust counterpart of formula (43) can be represented as where and are dual auxiliary variables of the robust counterpart model.

IV. AN IMPROVED MOEA HEURISTIC ALGORITHM
As a variant of the typical LRP, the MMRLRP is an NP-hard problem. To efficiently obtain Pareto optimal solutions, an improved MOEA heuristic algorithm is proposed, which combines a decomposition strategy with the framework of PICEA-g. The proposed algorithm outperforms other state-of-the-art MOEAs on multi-objective benchmark problems [46]. The Tchebycheff decomposition strategy is a widely used multi-objective optimization method to recombine neighbourhood solutions to improve the diversity and superiority of algorithms. The algorithm is denoted as PICEA-g-td, in which the key operators are described as follows.

A. GENETIC OPERATORS
PICEA-g-td is a general framework of evolutionary algorithms that is composed of the initial population, crossover operators, mutation operators, and selection operators based on the fitness of the solutions.

1) ENCODING SCHEME
To determine the MMRLRP, a natural number permutation encoding is used to represent the solution.

2) CROSSOVER OPERATORS
Considering the multi-period characteristic of the solution, crossover operations occur within one period of the chromosome. Crossover operations are used to operate a line vector , which is randomly selected from a chromosome, and a period of location-routing planning during the planning horizon. is as follows.

3) MUTATION OPERATORS
After the procedure of the crossover operations, mutation operations are used to operate a line vector , which is randomly selected from a chromosome and means a period of location-routing planning during the planning horizon. (1). For the sub-vectors Ch g1 τt , Ch g3 τt Ch g4 τt and , a reverse sequence mutation is used to generate offspring. We take the mutation of 1 as follows.
First, two inverse points are randomly selected within a sub-vector Ch g1 τt . 1 = 1 2 3 4 6 7 8 9 Second, the numbers between the two inverse points are inversed to obtain offspring. 1 = 1 2 3 7 6 4 8 9 (2). For the sub-vector Ch g2 τt , reverse sequence mutation, shift mutation and exchange mutation are randomly used to generate offspring. The exchange mutation operator is as follows.
First, two genes are randomly selected within the sub-vector 2 .

= 1 3 4 3 3 1 1 4 2
The shift mutation randomly selects a number within the sub-vector 2 and then shifts its content, which represents the TRDC that belongs to the vehicle in the corresponding location in the sub-vector 1 . First, one gene is randomly selected within the sub-vector 2 .

= 1 3 1 3 3 1 4 4 2
Second a random number is generated from the integers to M and shifts the selected number to obtain offspring.

1) FRAMEWORK OF PICEA-g-td
PICEA-g-td is a universal co-evolutionary algorithm in which the solutions co-evolve with a set of goal vectors. In the search process, the goal vectors are updated and guide the solutions towards the Pareto optimal front. The framework of PICEA-g-td is an μ λ elitist approach. At the beginning of the calculation, candidate solutions and goal vectors are generated at random. new candidate solutions are produced by genetic operators, and new goal vectors are generated randomly in solution space in each iteration. The last-generation solutions and goal vectors are combined with the new candidate solutions and goal vectors in a pool, with best candidate solutions and best goal vectors according to the fitness function [21,56].
Considering that the values of multiple objectives are different, PICEA-g-td does not directly use the objective value of a solution as the fitness value. The fitness assignment precept has two rules. First, a candidate solution gain a fitness value if they dominate a particular set of goal vectors, and the fitness value is shared between other solutions that also dominate these goals. Second, goal vectors gain fitness values that are inversely proportional to the number of candidate solutions that dominate the vector, and a goal vector has a low fitness when it is dominated by more candidate solutions [53,54]. The framework of PICEA-g-td is as follows: (Insert Figure 3 here)

2) A NEIGHBOURHOOD STRATEGY BASED ON THE TCHEBYCHEFF DECOMPOSITION
Researchers have applied evolutionary algorithms to handle the multi-objective optimization problems and have achieved great success. However, evolutionary algorithms have many difficulties in solving optimization problems that have more objectives [57]. The decomposition approach is demonstrated to have superiority in engineering applications [58]. To develop the treatment of PICEA-g, the Tchebycheff decomposition strategy is proposed to determine the neighbourhood of the evolutionary algorithm in this study. In the Tchebycheff decomposition approach, a set of uniformly distributed vectors λ 1 λ 2 … λ is generated and then used to convert a multi-objective optimization problem into a set of scalar optimization sub-problems, as follows: where λ = λ 1 λ 2 … λ with 1 1 m i i     , and λ is the weight vector of a sub-problem for a multi-objective problem. * = 1 * 2 * … * is an ideal reference vector, for which * . To obtain the optimal value of the solutions in each objective, the proposed algorithm, called the PICEA-g-td, constructs a neighbourhood using the Tchebycheff decomposition approach for each solution and ensures that only adjacent sub-problems can be used to optimize each other. In the PICEA-g-td, the neighbourhood strategy selects the nearest T weight vectors to build the neighbourhood for each weight vector depending on their Euclidean distance. Calculating the sub-problems of solutions, each solution selects a parent solution from the neighbourhood individuals independently. Only adjacent sub-problems can be used to optimize the solution and the excellent individuals of some dimension can be used to produce offspring.
The neighbourhood strategy based on the Tchebycheff decomposition is as follows: Step 1: Input solution set = 1 2 … and a set of uniformly distributed weight vectors λ = λ 1 λ 2 … λ ; Step 2: Calculate the objectives of all solutions and set the ideal reference point * = 1 * 2 * … * with * ; Step 3: Calculate sub-problems of solution = 1 2 … by Eq. (36); Step 4: Construct the neighbourhood = 1 2 … for each sub-problem by calculating the Euclidean distance between the weight vectors; Step 5: With of solutions as criteria, individuals in the neighbourhood of the sub-problem are selected as paired individuals for the individual = 1 2 … ; Step 6: Generate offspring by cross-operation between individual = 1 2 … and the paired individuals; Step 7: Output parent individuals and child individuals.

A. DESCRIPTION OF EXPERIMENTS AND PERFORMANCE OF PICEA-G-TD ALGORITHMS
To verify the performance of the PICEA-g-td algorithm, 24 test instances are randomly generated based on the test data sets C1, C2, R, and RC. For example, a problem called C1-3-5-50 means that five candidate locations to set up the TRDCs and 50 demand points are selected randomly for dataset C1, and the planning horizon consists of three periods. The capacity of the TRDCs is considered to be insufficient, and the rental cost for all TRDCs is 2000 Yuan. We assume that there are sufficient homogeneous vehicles to transport relief supplies and that these vehicles have a loading capacity 50. Without loss of generality, a package of relief supplies, including medicine and daily supplies, is assumed to be worth 50 Yuan and is sent to the demand points. For those problems, the nominal value of the relief demand is equal to the demand of the test data sets, and the maximum deviation from the nominal value is equal to 10% of the nominal value in each demand point. The cost of transportation is correlated with the distance between points and the transportation cost per unit distance c is equal to 1.7. Considering there could be a shortage of relief supplies, we consider that the penalty for a shortage of relief is ten times the value of the relief ( = 1 ).
The algorithm described in Section 4 is coded in MATLAB 2017, and all results are obtained using a 2.50 GHz Intel Core i5-6200U CPU with 8 GB of RAM running in Windows 10. All test instances were solved by NSGA-II, MOEA/D, PICEA-g, and PICEA-g-td. Considering that the MMRLRP problem is NP-hard, we cannot expect to obtain the set of all Pareto optimal solutions exactly. Ten runs are performed for the test instances to obtain the approximate solutions. The parameters of the four algorithms mentioned above are given in Table 1.
(Insert Table 1 here) To compare the efficiency of the algorithms, a widely used performance indicator called the set coverage (C-metric) is used to analyse Pareto approximation solutions obtained by those algorithms [59]. In a minimization multi-objective optimization problem, the C-metric can be described as follows: where and are two Pareto solutions obtained by the two algorithms. The definition of the is the percentage of the Pareto solutions in algorithm B that are dominated by at least one solution in the Pareto solutions in algorithm A. Note that might not necessarily is equal to 1 . By conducting 10 runs on 24 test instances, the other three algorithms are compared with the PICEA-g-td algorithm, and the median C-metric values are shown in Table 2.
(Insert Table 2 here) The superior results of the median C metric are in boldface. As seen from Table 2, the average values  ,  and are larger than , and , which means that more Pareto solutions obtained by the PICEA-g-td algorithm dominate Pareto solutions obtained by the other three algorithms but fewer Pareto solutions obtained by the other three algorithms dominate Pareto solutions obtained by the PICEA-g-td algorithm. Taking test instance C1-3-5-50 as an example, it can be seen that 13.83% of solutions obtained by PICEA-g are As seen from the results of the MMRLRP problem, the improved algorithm, named PICEA-g-td, is better than the MOEA/D and NSGA-II algorithms, which are two types of multi-objective algorithms that have wide application and excellent performance. Interestingly, the quality of Pareto solutions obtained by PICEA-g-td is better in terms of the total cos. Because the number of Pareto solutions obtained is small in the test instances with few demand points and planning periods, both PICEA-g-td and PICEA-g are prone to excessive convergence and obtained narrow C metric values. However, the results of the 24 test instances still show that the proposed algorithm outperforms the other algorithms.

B. TEST INSTANCES
This part addressed the epidemic logistics of Wuhan, China, during the COVID-19 outbreak (coronavirus disease 2019). To stop the spread of the disease, the Chinese government blocked Wuhan city. Providing supplies to module hospitals, which are important places for responding to the epidemic, poses a challenge to epidemic logistics when the demand for medical supplies changes over time, with the development of the epidemic. In this study, we mainly aim at the 16 module hospitals in Wuhan with 3 candidate TRDCs to supply humanitarian relief during the outbreak. The location of the candidate TRDCs and demand points are numbered from 1 to 3 and from 4 to 19, respectively, as shown in Figure 4. We consider the standard carton ( 36 × 26 × 3 cm 3 ) as a transportation unit that contains surgical masks, medicine and other relief for a patient. The relief is gathered from other cities in China, which is beyond the domain of this research.
Depending on some data released for public use by the government, relief-related information included the TRDC locations, the number of module hospital beds, and the nominal value of relief demand estimated by the number of beds in the module hospital. Although many parameters related to the MMRLRP model have been obtained, more parameters must be determined before solving TRDC locations and route planning for epidemic logistics. Because some data are not yet released for public use by the government, we insert man-made data to test the model and algorithm, which will not lead to essentially different results.
Assumptions and parameters are as follows: (1) In Wuhan city, three candidate TRDCs are Jieli logistics park, Cuiyuan cold chain logistics park, and Baowan logistics park. TRDCs existed before the COVID-19 outbreak, and the use costs can be estimated if a candidate TRDC is selected. The information on the TRDCs is shown in Table 3; (2) Considering that the emergency relief items are transported in a standard carton (unit), the loading capacity and other parameters of the vehicle are given in Table 4; (3) Taking 7 days as a planning horizon, we plan daily for relief supplies for module hospitals. Considering the development of the epidemic, we consider that the demand for module hospitals are stochastic; (4) In each period, the relief demand of each module hospital is related to the number of beds, as shown in Table 5. In the case study, we consider that the relief demand (such as masks and medicine) is uncertain. The present study considers one relief item for each patient [60]; (5) Using satellite maps, the distance between traffic vertexes can be assessed in Table 6; (6) To avoid the potential risk of virus transmission, as in practice, vehicles returns to the starting TRDC after all of the demand points are served, to prevent drivers from being quarantined. (Insert Figure 4 here) (Insert Table 3 here) (Insert Table 4 here) (Insert Table 5 here) (Insert Table 6 here) We have demonstrated the superiority of PICEA-g-td over NSGA-II, MOEA/D and PICEA-g in the previous section. Here, we take a test case to illustrate the decision design of the proposed robust optimization model when addressing uncertain demand parameters. First, we assume that the actual values of the uncertain demand are randomly assigned around the predicted nominal value and within the given interval. In other words, the actual demands are randomly generated into the interval .
The MMRLRP model considers the demand uncertainties and adopts the robust optimization model proposed by Bertsimas and Sim [21]. In the robust optimization model, both the budget of uncertainty Γ and the data variability are key parameters used to address these uncertainties. To test the influence of uncertain parameters on the robust model, a series of experiments with different uncertain parameters are designed. The budget of uncertainty Γ is set to the interval . Ten percent of the relief demand parameters in the planning window are uncertain but accurately predicted when the value of is equal to 0.1, whereas if the value of is , all of the relief demand parameters in the planning window are considered to deviate from these nominal demand values. Simultaneously, we test three demand variability (0.1, 0.2, and 0.3), which correspond to deviations of relief demand in epidemic logistics. According to the statistical results, the analysis conclusions are as follows: (1) The uncertain parameter have the great impact on the total cost. The Pareto solution obtained is better when the uncertainty coefficient α is small ( = .1 ). Increasing the budget of uncertainty can reduce the total cost when the uncertain parameter is large, which mainly benefit from reducing penalty costs in the shortage of emergency supplies; (2) increasing the budget of uncertainty has led to an increase in the total transportation time, due to the addition vehicles are used to ensure fewer shortages of relief. Because adjustment of transport scheme for each period in the MMRLRP model, the total transportation time is roughly stable; (3) Increasing the budget of uncertainty can bring a certain total cost and obtain fairer Pareto solutions. However, a big has little significance to reduce the total transportation time and improve the fairness distribution when the budget of uncertainty increases to a certain value(such as is greater than 0.3 when = . 2 ). So, it is important to accurately predict the relief demand of demand points for emergency supplies distribution. In addition, an appropriate budget of uncertainty can improve the performance of the Pareto solution set of the MMRLRP problem when the demand cannot be accurately predicted.
Taking the MMRLRP model with demand variability = .2 and budget of uncertainty = . as an example, the objectives of the Pareto solution is as follows: Fig. 8 illustrates the Pareto-optimal solution based on objective 1 (minimizing the total travel time), Fig. 9 describes the best solution of objective 2 (minimizing the total cost), and the Pareto-optimal solution depending on objective 3 is displayed in Fig. 10 (minimizing the disutility of relief service). Sixteen black circles represent the location of module hospitals, and three candidate TRDCs are represented by red, green and blue squares. In the MMRLRP model, routes are represented by closed lines with different colours. To show the results clearly, we demonstrate the decision results in the 2, 4, and 6 periods as follows.
(Insert Figure 8 here) Figure 8 shows the Pareto solution with the minimal total travel time, and the three objectives are 42.94 hours, 75318.21 Yuan, and 7.88 disutility units, respectively.
(Insert Figure 9 here) Figure 9 shows the Pareto solution with a minimal total cost, and the three objectives are 50.97 hours, 66571.45 Yuan, and 8.01 disutility units.
(Insert Figure 10 here) Figure 10 shows the Pareto solution with minimal the disutility of relief service, and the three objectives are 48.19 hours, 78238.09 Yuan, and 6.77 disutility units.
Comparing the three figures show that some routes appear in different solutions, such as the two routing decisions of Figs. 8, 9, and 10 in period 2 and the location decisions for 2 and 4 periods in three Pareto solutions. Different from the shortest path principle, some routes can be longer to select cheaper TRDC 3 in Figs. 9. However, most routing decisions are close to conventional routing decisions from considering time and economic factors. Because of the shortage of relief supplies, all of the TRDCs are used to provide relief for module hospitals at the early and mid-stage of an outbreak. We note that in the case of a shortage of relief supplies, all demand points are delivered depending on the set of fairness objective in the MMRLRP model. However, more vehicles start from TRDCs 1 and 3 which close to most module hospitals in the later period in Figs. 8 when relief supplies are sufficient.
The Pareto-optimal solutions based on objective 1 (minimizing the total travel time) use three TRDCs to provide relief to module hospitals, which is the principle that using multiple TRDCs in traditional location-routing decisions can reduce the total service time and total transport distance. In the solution given by Fig. 8, where the total travel time is accounted for 42.94 hours, and it took 8.03 hours less than the solutions given by Fig. 9. The Pareto-optimal solution based on objective 2(minimal total cost) spend 66571.45 Yuan to provide relief supplies for the module hospitals, and it is 8746.76 and 12667.45 Yuan less than the solutions given in Fig. 8 and 10 respectively. Considering the fairness of emergency logistics, the solution given by Fig. 10 tends to use more vehicles and is 14.08%, and 15.49% better than the solutions given by Figs. 8, and 9, respectively, on the objective of the disutility of relief service. By analyzing the multi-period decision of the Pareto solution, it is shown that the MMRLRP model can provide decision-makers with multiple options considering different objectives.

VI. CONCLUSIONS AND FUTURE STUDIES
In this paper, we present a multi-period multi-objective robust location-routing model with uncertain demand. During public health outbreaks, the total travel time, the total cost, and the disutility of relief service are considered as optimization objectives throughout the planning horizon. Considering the stochastic demand, this study uses the robust optimization method proposed by Bertsimas and Sim [21] to address the uncertainties. By introducing the decomposition strategy, we propose PICEA-g-td to solve the MMRLRP models. Twenty-four sets of instances with different demands and periods are generated randomly to evaluate the performance of the proposed PICEA-g-td. The comparison results showed that PICEA-g-td outperforms PICEA-g, MOEA/D and NSGA-II in most cases. To further evaluate the performance of the MMRLRP, a numerical study with the epidemic logistics of the COVID-19 epidemic in Wuhan, China, is conducted to illustrate the applicability of the proposed model. Experimental results have shown that the parameter settings of robust optimization have improved the robustness of uncertain epidemic logistics systems.
In the future, the following three research directions could be explored. First, we need to consider disasters with other characteristics, such as the reliability of routes. Second, it is necessary to provide relief supplies while evacuating the affected people, which makes the research more realistic in the application of disaster scenarios. Third, we do not propose a prediction method to forecast the uncertainty value and do not account for the uncertainty other than the relief demand. Other constraints of the relief location-distribution model, such as the capacity of different vehicles and the capacity of the TRDCs, can also be accounted for in post-disaster conditions in relief distribution.