# A Multiobjective Optimization Model for Continuous Allocation of Emergency Rescue Materials

Authors: Yonghong Liu; Yucheng Li; De Huang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Emergency rescue operations play a vital role in alleviating human suffering, reducing casualties, and cutting down economic losses. One key aspect in the management of these operations is the rational allocation of emergency relief materials, where the allocation is continuous, dynamic, and concurrent. This allocation should be made not only to minimize the emergency rescue losses, but also to reduce the cost of emergency rescue work. A reasonable and effective allocation scheme for emergency relief materials can be established to adapt to the continuity, dynamics, and concurrency of material distribution. In this work, we propose a multiobjective optimization model of emergency material allocation with continuous time-varying supply and demand constraints, where the objective is to minimize the losses and the economic cost incurred by the emergency rescue operations. The constrained optimization problem is handled through sequential unconstrained minimization techniques, and the multiobjective optimization is carried out by the fast nondominated sorting genetic algorithm (NSGA-II) with an elite strategy to obtain a Pareto solution set with fairness and balance of loss and cost. The loss and cost associated with the Pareto frontier are employed to find an appropriate noninferior solution and its corresponding material allocation scheme. We verify through several simulations the model feasibility and the effectiveness of the proposed method, which can provide decision support for continuous material allocation in emergency rescue operations.


---

## Body

## 1. Introduction

Human populations are frequently impacted by natural disasters such as earthquakes, tsunamis, floods, typhoons, volcanic eruptions, and debris flow, as well as major accidents such as those of mines, traffic, and fires. Such disasters and accidents not only harm the lives and properties of people, but also seriously hinder economic growth and society development [1, 2]. The emergency rescue operations, in response to such natural disasters and accidents, should be highly active, effective, and quick. These requirements promoted extensive research and investigations on the issues of coordination among emergency rescue units, rescue route selection, and material distribution. In particular, effective allocation of emergency materials is one of the most essential and key problems in rescue operations [3, 4]. A reasonable and effective material allocation scheme can complete the emergency rescue work at the lowest cost and loss. More importantly, such a scheme can also reduce, to a certain extent, the personal and property losses caused by disasters and accidents [5, 6].Fiedrich et al. [7] introduced a dynamic optimization model to formulate the problem of minimizing the total number of deaths after an earthquake. This problem was solved using simulated annealing and Tabu search methods to determine the best resource allocation plan for earthquake sites. Balcik and Beamon [8] used mathematical models to investigate the facility location problem in humanitarian relief chains. Hence, they determined the numbers and locations of distribution centres in relief networks, as well as the quantities of needed relief supplies for disaster-affected areas. Li et al. [9] reviewed coverage models and optimization techniques for the location and planning of emergency response facilities. The authors particularly reviewed state-of-the-art models including the hypercube queuing model, the dynamic allocation model, the progressive coverage model, and the cooperative coverage model. Wex et al. [10] developed a decision-support model for rescue unit allocation and arrangement. The model seeks to minimize the sum of rescue completion times weighted by incident severity levels. A metaheuristic algorithm was used to solve this optimization problem and significantly reduce rescue times. Coordination among rescue units and independent agencies for fast material delivery was also discussed as a potential research direction. Yang et al. [11] proposed a decision-support system for public safety and emergency management in patrol service centres. This system shows a high potential for the allocation, deployment, and scheduling of resources in order to improve the efficiency of patrol-related resources and services and boost emergency response capabilities. However, this system is not suitable for real-time or dynamic scenarios. Zhou et al. [12] proposed a multiobjective optimization model for the multiperiod dynamic emergency resource scheduling problem. Schryen et al. [13] proposed a method to support decision makers in sudden disaster situations with high risk, time pressure, uncertainty, conflict, or lack of information. The experimental results showed small differences between heuristic-based and optimal solutions over short periods of time. Chen and Yu [14] used network partitioning and integer programming to find temporary locations for the transportation infrastructure of emergency medical services (EMS). Taking into consideration worst-case EMS demands that arise after the onset of a disaster, the authors expect that their analysis results can help decision makers with site selection issues and can serve as a benchmark for planning postdisaster emergency services. Duhamel et al. [15] proposed a mathematical model and a heuristic algorithm for a multiperiod location-allocation problem in postdisaster relief operations. They also related the proposed quantitative model to another model for disaster-resistant systems. Pradhananga et al. [16] proposed a three-level network model for integrated emergency preparedness, response planning, and emergency material distribution. This model seeks to minimize social costs through determining a set of potential resource supply points at the highest level. The results showed that prearrangement of some items and postdisaster procurement can reduce the shortage of emergency supplies. Alem et al. [17] introduced a new two-stage stochastic network traffic model for providing rapid humanitarian assistance to victims after a disaster. This model takes into account typically overlooked elements in rescue operations, such as budget allocation, fleet sizing, and procurement and lead times of dynamic multicycle changes. Sung and Lee [18] modelled the ambulance routing problem as a resource-constrained classification problem. Then, they determined the source and destination hospitals for evacuating patients. A column generation algorithm was used to solve the formulated fixed-priority resource allocation problem. Though the best obtained solution may require complicated operations, it may still be easy to follow and implement in a postdisaster environment. Bai [19] proposed a new two-stage optimization method for emergency material distribution with multiple vendors, disaster areas, relief resources, and vehicles. Criteria of fairness, timeliness, and economic efficiency are considered in the formulation of this multiobjective optimization model for facility location, vehicle routing, and resource allocation decisions. The goals of this model are to minimize the proportion of unsatisfied demands, minimize the response time of emergency relief, and minimize the total cost of the entire process. Deng et al. [20] investigated the problem of two-stage random-capacity-limit location-allocation in emergency logistics. To solve this problem, the authors established a random expected-value model and its deterministic counterpart. Then, they developed an improved particle swarm optimization algorithm with restart of the Gaussian cloud operator and an adaptive parameter selection strategy.Boonmee et al. [21] studied facility location problems in emergency humanitarian logistics, where they analyzed facility location types, disaster characteristics, and decision-making approaches. They also investigated data models, typical formulations of objectives and constraints, solution methods, practical application procedures, and case studies. Liu et al. [22] addressed the allocation of emergency relief materials after chemical leakage accidents in rivers. The authors proposed a material allocation framework that minimizes the response time for dynamic emergency material demands. They established a multiobjective emergency resource allocation model that seeks to minimize both the total resource allocation cost and the total loss caused by resource insufficiency. Particle swarm optimization was used to obtain the optimal solutions. Nevertheless, the emergency resource allocation process did not consider the response time, resource urgency, and disaster severity. Guo et al. [23] addressed the problem of distributing marine emergency rescue materials. The authors formulated this multiple-resource allocation problem within a multiobjective optimization framework that considers the long-range nature of maritime rescue operations and sustainable recycling of emergency resources. The fast nondominated sorting genetic algorithm (NSGA-II) was used for solving this optimization problem, and the candidate solution set approaching the Pareto frontier was exploited to find a suitable noninferior resource allocation scheme. The authors adopted a discrete resource allocation problem formulation with decreasing proportional constraints and proposed expert-system mathematical models for solving this problem. This formulation has the advantage of a decreasing ratio which can prevent unfair material allocation, but it did not account for continuous material distribution.Therefore, material allocation is generally formulated as a multiobjective and multiconstrained optimization problem, for which reasonable and effective solutions are crucially needed to ensure the success of emergency rescue work. At present, existing approaches for emergency relief material allocation focus on selecting the locations of the relief materials. An allocation scheme is typically determined based on known locations of material storehouses and material storage capacities. However, a well-designed allocation scheme for emergency relief materials needs to consider the emergency rescue time, economic costs, the loss of materials that are undeliverable on time, and the balance of the material distribution plan. Therefore, this paper comprehensively considers the continuity, dynamics, and concurrency of emergency relief material allocation. Continuity of material allocation is defined for materials that are genuinely necessary and continuously consumed in emergency rescue sites, such as water, medicines, and food. Towards continuous material allocation, we construct a multiobjective optimization model, based on a time-varying order of material transportation, for the minimization of the conflicting factors of the cost and loss of the emergency rescue operations. Based on existing methods and algorithms for emergency material allocation with time-varying order of material transportation, we choose to minimize two conflicting objectives, namely, the economic costs and the losses. Also, we adopt the NSGA-II heuristic algorithm to solve this optimization problem and obtain Pareto frontier solutions. From this frontier, a noninferior solution is selected to get an economical, reasonable, and effective emergency relief material allocation plan that achieves fairness and balance between cost and loss. Such a plan can provide decision support for decision makers on the best continuous allocation schemes for emergency rescue materials. The experimental results show that the solution of the emergency material distribution problem is feasible based on criteria of fairness and balance of material distribution. Also, the solution shows characteristics of timeliness, time-varying dynamics, and continuity.The remainder of the paper is structured as follows. The proposed multiobjective material allocation optimization model is discussed in Section2. Section 3 presents the algorithmic solution steps of that optimization problem. Section 4 provides case studies of the emergency rescue material allocation. In Section 5, the experimental results are presented and discussed. Finally, conclusions are drawn in Section 6.


## 2. Multiobjective Optimization for Emergency Material Allocation

### 2.1. Problem Description

After a natural disaster or a major accident, emergency rescue decision makers should direct rescue teams for people rescue missions and material allocation. A reasonable and effective allocation plan of relief materials should have the shortest time and the lowest cost, while minimizing the losses caused by the natural disaster or accident. Indeed, creating a multiobjective optimization model with theses optimality criteria can lead to sought a reasonable and effective emergency rescue plan.The types, total quantities, and the total number of storage locations of emergency relief materials are typically known. Letn be the number of the candidate storage units for the emergency rescue materials. This number is set based on the estimated emergency rescue time, the time when each storage unit arrives at the rescue site, and the number of units that fail to deliver materials within the expected rescue time. The rescue time T, the total number of rescue material types m, and the quantity of each material type qj(j=1,2,…,m) are estimated based on the emergency damage level. Disregarding interruptions in material allocation, let the shortest time for transporting the material storage unit Bi to the rescue location A be ti. We seek to determine n candidate storage units for emergency rescue materials and complete the allocation of the materials in time T. Based on time-varying sorting, the times ti at which the candidate storage units Bi (i = 1, 2, …, n) arrive at the rescue locations are sorted such that ti+1≥ti. The capacity for the material type qj provided by the storage unit Bi is denoted by Sij, where i=1,2,…,n;j=1,2,…,m. Let the decision variable xij indicate whether a quantity of the material type qj is assigned from the material storage unit Bi to the accident location A. The problem of the emergency relief material allocation should be solved for xij in order to achieve the shortest rescue time, the lowest cost of material allocation, and the minimization of losses caused by the emergency. The notation used in the model is presented in Table 1.Table 1
Notation and symbols.
nThe total number of storage units that can provide emergency rescue materialsmTypes of emergency rescue materialsTEstimated emergency rescue timeAEmergency rescue locationBiThei-th storage unit of emergency rescue materials, i = 1, 2, …, nqjEstimated quantity of the rescue material of typej, j=1,2,…,mtiTravel time from the storage unitBi to the rescue location A, i = 1, 2, …, nLjThe loss caused by the lack of the material of typej in unit time, j=1,2,…,mvjThe amount of the material of typej consumed per unit of time, j=1,2,…,mSijThe storage capacity for the material of typej provided by the storage unit Bi, i = 1, 2, …, n, j=1,2,…,mCijThe cost of sending one unit of the material of typej from the storage unit Bi to the rescue location A, i = 1, 2, …, n, j=1,2,…,mxijThe demand for the material of typej from the storage unit Bi to the rescue location A, i = 1, 2, …, n, j=1,2,…,mRijThe remaining amount of the material of typej from the storage unit Bi to the rescue location A, when the time is ti, i = 1, 2, …, n, j=1,2,…,m


### 2.2. Proposed Mathematical Model

Based on the above problem description, the rescue time is limited by the distance of the material rescue storage unit to the accident location. As a practical allocation scheme for emergency rescue materials at the minimum cost while also minimizing the emergency losses, we propose the following multiobjective mathematical optimization model for continuous allocation of emergency rescue materials. The objective functions are(1)minf1=∑i=1n∑j=1mCij⋅xij,(2)minf2=∑i=1n∑j=1mFvj⋅ti−ti−1−xi−1j−Ri−1j⋅Lj⋅ti−ti−1+∑j=1mLj⋅t1⋅vj,where(3)Fx=0,x≤0,x,x>0,(4)Rij=xi−1j+FRi−1j−vj⋅ti−ti−1,∀1≤i≤n,1≤j≤m,R0j=0,x0j=0,t0=0.The constraints are(5)∑i=1nSij≥qj≥T⋅vj,∀1≤i≤n,1≤j≤m,(6)0≤xij≤Sij,∀1≤i≤n,1≤j≤m,(7)∑i=1nxij≥qj,∀1≤i≤n,1≤j≤m,(8)T≥ti+1≥ti,∀1≤i≤n,(9)∀1≤i≤n,1≤j≤m,T>0,n≥1,m≥1,qj>0,Lj≥0,vj≥0,Cij≥0.The objective function (1) uses the coefficient Cij as a variable objective-dependent cost. In particular, the choice of this coefficient depends on the relative magnitude of the relief material cost and the minimum material deployment time. On the one hand, if the objective is to minimize the deployment time, the transportation cost will be higher. On the other hand, if the objective is to minimize the transportation cost, the deployment time may be longer, and hence the loss can be larger. Practically, the appropriate transportation mode should be selected according to the actual situation, in order to appropriately determine the coefficient Cij. The objective function (2) represents the losses inflected by the emergency on the people, environment, or economy. These losses are determined based on the nature and degree of the emergency. For the emergency rescue operations, the main influencing factors include the amount of available materials, the loss caused by the lack of materials, and the consumption of materials. In the model mentioned above, the loss includes not only the economic and property losses, but also the damage caused by the lack of life-saving essentials. The loss caused by the lack of materials is based on the magnitude of the material influence on people, animals, or environments. This loss is quantified as a value in the unit interval [0, 1].For the objective (2), Fvj⋅ti−ti−1−xi−1j−Ri−1j represents the amount of missing material of type j after the storage unit Bi arrives at the emergency site. The period ti−ti−1 represents the time difference between the two latest (current and previous) arrival times of the materials from the storage unit Bi. In this period, the amount of consumed material of type j is vj⋅ti−ti−1, the remaining amount of this material at the last moment is Ri−1j, and the amount of delivered material is xi−1j. From these three amounts, the amount of lost material can be obtained. The amount of the remaining material is given by (4). When the sum of the amount of the delivered material xi−1j and the amount of the surplus material Ri−1j is greater than the amount of the consumed material vj⋅ti−ti−1, then the amount of missing materials is zero. Otherwise, the amount of missing materials is the difference between the consumed materials and the sum of the delivered and remaining materials. The constraint (5) indicates that the demand for a material of type j must be more than or equal to the consumption within the rescue time T. Also, under conditions of material sufficiency, the material demand needs to be less than or equal to the sum of the materials in the storage units. The constraint (6) indicates that the amounts of various rescue materials deployed from different material storage units cannot exceed the amounts in the storage units. The constraint (7) shows that the total sum of the relief supplies should be more than the demand of the relief materials. The constraint (8) indicates that the longest time for a material storage unit to arrive at a rescue site cannot exceed the expected rescue time, where the material storage units are sorted by the material transportation times. Side constraints on the variables are given in (9).


## 2.1. Problem Description

After a natural disaster or a major accident, emergency rescue decision makers should direct rescue teams for people rescue missions and material allocation. A reasonable and effective allocation plan of relief materials should have the shortest time and the lowest cost, while minimizing the losses caused by the natural disaster or accident. Indeed, creating a multiobjective optimization model with theses optimality criteria can lead to sought a reasonable and effective emergency rescue plan.The types, total quantities, and the total number of storage locations of emergency relief materials are typically known. Letn be the number of the candidate storage units for the emergency rescue materials. This number is set based on the estimated emergency rescue time, the time when each storage unit arrives at the rescue site, and the number of units that fail to deliver materials within the expected rescue time. The rescue time T, the total number of rescue material types m, and the quantity of each material type qj(j=1,2,…,m) are estimated based on the emergency damage level. Disregarding interruptions in material allocation, let the shortest time for transporting the material storage unit Bi to the rescue location A be ti. We seek to determine n candidate storage units for emergency rescue materials and complete the allocation of the materials in time T. Based on time-varying sorting, the times ti at which the candidate storage units Bi (i = 1, 2, …, n) arrive at the rescue locations are sorted such that ti+1≥ti. The capacity for the material type qj provided by the storage unit Bi is denoted by Sij, where i=1,2,…,n;j=1,2,…,m. Let the decision variable xij indicate whether a quantity of the material type qj is assigned from the material storage unit Bi to the accident location A. The problem of the emergency relief material allocation should be solved for xij in order to achieve the shortest rescue time, the lowest cost of material allocation, and the minimization of losses caused by the emergency. The notation used in the model is presented in Table 1.Table 1
Notation and symbols.
nThe total number of storage units that can provide emergency rescue materialsmTypes of emergency rescue materialsTEstimated emergency rescue timeAEmergency rescue locationBiThei-th storage unit of emergency rescue materials, i = 1, 2, …, nqjEstimated quantity of the rescue material of typej, j=1,2,…,mtiTravel time from the storage unitBi to the rescue location A, i = 1, 2, …, nLjThe loss caused by the lack of the material of typej in unit time, j=1,2,…,mvjThe amount of the material of typej consumed per unit of time, j=1,2,…,mSijThe storage capacity for the material of typej provided by the storage unit Bi, i = 1, 2, …, n, j=1,2,…,mCijThe cost of sending one unit of the material of typej from the storage unit Bi to the rescue location A, i = 1, 2, …, n, j=1,2,…,mxijThe demand for the material of typej from the storage unit Bi to the rescue location A, i = 1, 2, …, n, j=1,2,…,mRijThe remaining amount of the material of typej from the storage unit Bi to the rescue location A, when the time is ti, i = 1, 2, …, n, j=1,2,…,m


## 2.2. Proposed Mathematical Model

Based on the above problem description, the rescue time is limited by the distance of the material rescue storage unit to the accident location. As a practical allocation scheme for emergency rescue materials at the minimum cost while also minimizing the emergency losses, we propose the following multiobjective mathematical optimization model for continuous allocation of emergency rescue materials. The objective functions are(1)minf1=∑i=1n∑j=1mCij⋅xij,(2)minf2=∑i=1n∑j=1mFvj⋅ti−ti−1−xi−1j−Ri−1j⋅Lj⋅ti−ti−1+∑j=1mLj⋅t1⋅vj,where(3)Fx=0,x≤0,x,x>0,(4)Rij=xi−1j+FRi−1j−vj⋅ti−ti−1,∀1≤i≤n,1≤j≤m,R0j=0,x0j=0,t0=0.The constraints are(5)∑i=1nSij≥qj≥T⋅vj,∀1≤i≤n,1≤j≤m,(6)0≤xij≤Sij,∀1≤i≤n,1≤j≤m,(7)∑i=1nxij≥qj,∀1≤i≤n,1≤j≤m,(8)T≥ti+1≥ti,∀1≤i≤n,(9)∀1≤i≤n,1≤j≤m,T>0,n≥1,m≥1,qj>0,Lj≥0,vj≥0,Cij≥0.The objective function (1) uses the coefficient Cij as a variable objective-dependent cost. In particular, the choice of this coefficient depends on the relative magnitude of the relief material cost and the minimum material deployment time. On the one hand, if the objective is to minimize the deployment time, the transportation cost will be higher. On the other hand, if the objective is to minimize the transportation cost, the deployment time may be longer, and hence the loss can be larger. Practically, the appropriate transportation mode should be selected according to the actual situation, in order to appropriately determine the coefficient Cij. The objective function (2) represents the losses inflected by the emergency on the people, environment, or economy. These losses are determined based on the nature and degree of the emergency. For the emergency rescue operations, the main influencing factors include the amount of available materials, the loss caused by the lack of materials, and the consumption of materials. In the model mentioned above, the loss includes not only the economic and property losses, but also the damage caused by the lack of life-saving essentials. The loss caused by the lack of materials is based on the magnitude of the material influence on people, animals, or environments. This loss is quantified as a value in the unit interval [0, 1].For the objective (2), Fvj⋅ti−ti−1−xi−1j−Ri−1j represents the amount of missing material of type j after the storage unit Bi arrives at the emergency site. The period ti−ti−1 represents the time difference between the two latest (current and previous) arrival times of the materials from the storage unit Bi. In this period, the amount of consumed material of type j is vj⋅ti−ti−1, the remaining amount of this material at the last moment is Ri−1j, and the amount of delivered material is xi−1j. From these three amounts, the amount of lost material can be obtained. The amount of the remaining material is given by (4). When the sum of the amount of the delivered material xi−1j and the amount of the surplus material Ri−1j is greater than the amount of the consumed material vj⋅ti−ti−1, then the amount of missing materials is zero. Otherwise, the amount of missing materials is the difference between the consumed materials and the sum of the delivered and remaining materials. The constraint (5) indicates that the demand for a material of type j must be more than or equal to the consumption within the rescue time T. Also, under conditions of material sufficiency, the material demand needs to be less than or equal to the sum of the materials in the storage units. The constraint (6) indicates that the amounts of various rescue materials deployed from different material storage units cannot exceed the amounts in the storage units. The constraint (7) shows that the total sum of the relief supplies should be more than the demand of the relief materials. The constraint (8) indicates that the longest time for a material storage unit to arrive at a rescue site cannot exceed the expected rescue time, where the material storage units are sorted by the material transportation times. Side constraints on the variables are given in (9).


## 3. Solution of the Multiobjective Material Allocation Problem

According to the established mathematical model (1)–(9), the two objectives (1) and (2) conflict with each other. Hence, it is difficult to achieve the minimum transportation cost and loss at the same time. For such multiobjective optimization problems, Pareto-optimal solutions can be obtained through swarm intelligence heuristic-search stochastic optimization algorithms. We analyzed and compared the improved harmony search algorithm (IHS) [24], the multiobjective evolutionary algorithm based on decomposition (MOEA/D) [25], and the fast nondominated sorting genetic algorithm (NSGA-II) with an elite strategy [26]. The NSGA-II method appears more appropriate and effective in solving the established material allocation problem. Indeed, the NSGA-II method typically shows the best performance and has been widely used in resource allocation and inventory management problems [27–29]. As well, this method has been effective in solving multiobjective optimization problems, especially for 2-3 objective functions [30, 31]. This method essentially uses fast nondominated sorting and a diversity-preserving elite strategy to sort the parent and progeny populations at different levels of the nondominated solution frontier. A density estimation method based on the crowding distance is adopted to maintain population diversity.


### 3.1. NSGA-II Algorithm

Fast nondominated sorting of NSGA-II reduces the computational complexity of this algorithm. Moreover, this low-complexity algorithm uses a congestion degree and a congestion degree comparison operator to replace the fitness sharing strategy which needs to specify a shared radius. As well, NSGA-II sets the winning criterion in the same level of comparison after sorting, so that individual solutions in the Pareto domain can be extended to the whole domain in an evenly distributed diversified manner. In addition, NSGA-II introduces an elite strategy, which enlarges the sampling space, combines the parent population with its offspring populations, and creates competition for producing a best-fit next-generation population. By layering all individuals in the population, the best individuals will survive, and the population quality will be raised rapidly [32].Based on the NSGA-II algorithm, the steps for solving the multiobjective optimization problem of emergency material allocation are given as follows:Step 1. Model Parameter Setting and Initialization. Let T be the predicted emergency rescue time, m be the number of rescue material types, n be the number of candidate material storage units, νj be the consumption of the j-th rescue material per unit time, and Lj be the loss caused by the lack of the j-th rescue material. Based on the maximum time ti required for each material repository to reach the rescue site, the material repositories Bi are ranked. Thus, for each repository Bi and material type j, the storage capacity Sij and the storage cost Cij are obtained. The material demand qj is initialized according to the material storage, material consumption, and estimated rescue completion time.Step 2. Initialization of the NSGA-II Parameters. The following NSGA-II parameters are initialized: the number of individuals NI, the number of generations NG, the crossover probability cr, the crossover distribution index Dcr, the mutation probability mr, and the mutation distribution index Dmr.Step 3. Coding Strategy and Chromosome Design. The sought design variables are all of the decision variables xij for allocating a quantity of the j-th material from the i-th rescue material storage unit. According to the parameter setting and initialization [30, 31, 33], the chromosome associated with each individual solution of the overall population is encoded as X=x11,…,x1m,x21,…,x2m,…,xij,…,xnm. Each chromosome is initialized by real numbers according to the range set by the constraints in (6).Step 4. Nonlinear Constraint Processing. The sequential unconstrained minimization technique is used to deal with the nonlinear constraints in (7), and hence obtain a multiobjective unconstrained minimization problem [34]. The unconstrained objective is obtained from the constrained objective and the nonlinear constraints as follows [34]:(10)Fx=fx,ifgjx≥0∀j=1,2,…,m,fmax+∑j=1mgjx,otherwise,wheregjx is the j-th constraint, m is the number of constraints, fx is the constrained objective function, fmax is the objective function value corresponding to the worst feasible solution in the population, <gjx> is a penalty function of gjx, and the operator <·> denotes the absolute value of the operand. The definition of the penalty function here is different from those of the traditional methods and earlier genetic algorithms. The following rules of solution preference are adopted: (1) any feasible solution is superior to any infeasible solution; (2) among two feasible solutions, the solution with the better objective function value is preferred; (3) among two infeasible solutions, the solution with less conflict of constraints is preferred.Step 5. Fitness Calculation. Before calculating the fitness of a solution, the chromosome coding in Step 3 needs to be decoded into an n×m matrix, as shown in Figure 1. The matrix rows represent the amounts of demand for each emergency supply repository, and the matrix columns indicate the amounts of each material that needs to be allocated from each repository. According to the decoded decision xij, the fitness values are obtained by calculating the objective functions (1) and (2).Step 6. Fast Nondominated Sorting. According to the order of dominance (or preference) of the Pareto solutions, the population is divided into subpopulation levels. The similarity between subpopulations is evaluated on the Pareto frontier, and the resulting subpopulation and similarity measures are used to promote the diversity of nondominant solutions. The calculated fitness values are stratified by the fast nondominated sorting method. Thus, the population density near a specific solution can be judged. The crowding degree of each individual in the nondominated frontier is calculated according to the difference between the fitness values. So, the population density near a specific solution in a population is determined, and the density can be obtained through pair-wise competition and selection.Step 7. Crossmutation Genetic Operations. According to the initial crossover and mutation rates, the crossover and mutation genetic operations are carried out on the new population with the smallest population density. New individuals are generated, subjected to a cross mutation operation, and then merged with the parent population to obtain a new population.Step 8. Termination and Output. The fitness value of Step 5 is calculated for the individual population in Step 7, and the ranking method by Step 6 is used to reorder and stratify the population, in order to generate a new population in Step 7. This process is iteratively repeated until either the preset termination condition or the maximum number of iterations is reached. Hence, the optimal Pareto frontier is output.Figure 1
Decoding process of the chromosome encoding.


### 3.2. IHS Algorithm

The IHS multiobjective optimization method [24] is used to solve nonlinear nonconvex problems. This method generates new solution vectors to improve the accuracy and convergence speed of the harmony search algorithm (HS). The HIS-based solution steps for the optimization problem of emergency rescue multitarget material allocation are as follows:Step 1. The same as that of NSGA-II.Step 2. Initialization of the IHS Parameters. The following IHS parameters are initialized: the number of individuals NI, the number of generations NG, the probability of choosing from memory pr, the minimum pitch adjustment rate pmin, and the maximum pitch adjustment rate pmax.Step 3 to Step 5. The same as those of NSGA-II.Step 6. Improvisation of a New Harmony. Based on memory consideration, pitch adjustment, and random selection, the new harmony X′=x11′,…,xnm′ is generated as follows:(1)Memory consideration:(11)xim′⟵xim′±rd⋅bw.(2)
Pitch adjustment decision:(12)xij′⟵yeswith probabilitypminnowith probabilitypmax,j=1,2,…,m.(3)
Random selection:If the pitch adjustment decision forxij′=yes, then(13)xim′⟵xim′±rd⋅bw,wherebw is an arbitrary distance bandwidth and rd is a random floating number between 0 and 1.Step 7. Solution Update. Iterate over xi, where i = 1, 2, …, n, and exclude the worst existing harmony from the solution X.Step 8. Termination and Output. If the preset termination condition or the maximum number of iterations is reached, the algorithm is terminated. Otherwise, Steps 6 to 8 are repeated.


### 3.3. MOEA/D Algorithm

In this algorithm, a multiobjective optimization problem is decomposed into several scalar optimization subproblems, which are then simultaneously optimized. The solution of each optimization subproblem uses only information from its adjacent subproblems. This significantly reduces the computational complexity of each MOEA/D generation. The MOEA/D-based solution steps for the optimization problem of emergency rescue multitarget material allocation are as follows:Step 1. The same as that of NSGA-II.Step 2. Initialization of the MOEA/D Parameters. The following MOEA/D parameters are initialized: the number of individuals NI, the number of generations NG, the size of the weight neighborhood nd, and set the external population EP=Φ.Step 3. Distance Computation. For each i = 1, 2, …, NI, set ki=i1,…,ind, where ωi1,…,ωind are the nd closest weight vectors to ωi.Step 4. Initial Population Generation. Generate randomly an initial population X1,...,XNI, set FVi=fjXi,j=1,2, and initialize z=f1,f2T using equations (1) and (2).Step 5. Reproduction. Randomly select two serial numbers b and l from ki. Use genetic operators to generate a new solution y from Xb and Xl. Based on equations (1) and (2), apply the improvement heuristic on y to produce y′.Step 6. Update of the Objective Values z and Neighbouring Solutions. For each j=1,2,ifzj<fjy′, then set zj=fjy′. If gtey′λj,z≤gtexjλj,z, j∈Bi, then set Xj=y′ and FVj=fjy′.Step 7. External Population Update. For j = 1 and 2, remove from the external population (EP) all of the vectors dominated by fjy. Also, add fjy to EP if no vectors in EP dominate fjy. For each i = 1, 2, …, NI, loop through the execution of Step 4 to Step 7.Step 8. Termination and Output. If any of the stopping criteria are satisfied, then stop and output EP. Otherwise, go to Step 3.


## 3.1. NSGA-II Algorithm

Fast nondominated sorting of NSGA-II reduces the computational complexity of this algorithm. Moreover, this low-complexity algorithm uses a congestion degree and a congestion degree comparison operator to replace the fitness sharing strategy which needs to specify a shared radius. As well, NSGA-II sets the winning criterion in the same level of comparison after sorting, so that individual solutions in the Pareto domain can be extended to the whole domain in an evenly distributed diversified manner. In addition, NSGA-II introduces an elite strategy, which enlarges the sampling space, combines the parent population with its offspring populations, and creates competition for producing a best-fit next-generation population. By layering all individuals in the population, the best individuals will survive, and the population quality will be raised rapidly [32].Based on the NSGA-II algorithm, the steps for solving the multiobjective optimization problem of emergency material allocation are given as follows:Step 1. Model Parameter Setting and Initialization. Let T be the predicted emergency rescue time, m be the number of rescue material types, n be the number of candidate material storage units, νj be the consumption of the j-th rescue material per unit time, and Lj be the loss caused by the lack of the j-th rescue material. Based on the maximum time ti required for each material repository to reach the rescue site, the material repositories Bi are ranked. Thus, for each repository Bi and material type j, the storage capacity Sij and the storage cost Cij are obtained. The material demand qj is initialized according to the material storage, material consumption, and estimated rescue completion time.Step 2. Initialization of the NSGA-II Parameters. The following NSGA-II parameters are initialized: the number of individuals NI, the number of generations NG, the crossover probability cr, the crossover distribution index Dcr, the mutation probability mr, and the mutation distribution index Dmr.Step 3. Coding Strategy and Chromosome Design. The sought design variables are all of the decision variables xij for allocating a quantity of the j-th material from the i-th rescue material storage unit. According to the parameter setting and initialization [30, 31, 33], the chromosome associated with each individual solution of the overall population is encoded as X=x11,…,x1m,x21,…,x2m,…,xij,…,xnm. Each chromosome is initialized by real numbers according to the range set by the constraints in (6).Step 4. Nonlinear Constraint Processing. The sequential unconstrained minimization technique is used to deal with the nonlinear constraints in (7), and hence obtain a multiobjective unconstrained minimization problem [34]. The unconstrained objective is obtained from the constrained objective and the nonlinear constraints as follows [34]:(10)Fx=fx,ifgjx≥0∀j=1,2,…,m,fmax+∑j=1mgjx,otherwise,wheregjx is the j-th constraint, m is the number of constraints, fx is the constrained objective function, fmax is the objective function value corresponding to the worst feasible solution in the population, <gjx> is a penalty function of gjx, and the operator <·> denotes the absolute value of the operand. The definition of the penalty function here is different from those of the traditional methods and earlier genetic algorithms. The following rules of solution preference are adopted: (1) any feasible solution is superior to any infeasible solution; (2) among two feasible solutions, the solution with the better objective function value is preferred; (3) among two infeasible solutions, the solution with less conflict of constraints is preferred.Step 5. Fitness Calculation. Before calculating the fitness of a solution, the chromosome coding in Step 3 needs to be decoded into an n×m matrix, as shown in Figure 1. The matrix rows represent the amounts of demand for each emergency supply repository, and the matrix columns indicate the amounts of each material that needs to be allocated from each repository. According to the decoded decision xij, the fitness values are obtained by calculating the objective functions (1) and (2).Step 6. Fast Nondominated Sorting. According to the order of dominance (or preference) of the Pareto solutions, the population is divided into subpopulation levels. The similarity between subpopulations is evaluated on the Pareto frontier, and the resulting subpopulation and similarity measures are used to promote the diversity of nondominant solutions. The calculated fitness values are stratified by the fast nondominated sorting method. Thus, the population density near a specific solution can be judged. The crowding degree of each individual in the nondominated frontier is calculated according to the difference between the fitness values. So, the population density near a specific solution in a population is determined, and the density can be obtained through pair-wise competition and selection.Step 7. Crossmutation Genetic Operations. According to the initial crossover and mutation rates, the crossover and mutation genetic operations are carried out on the new population with the smallest population density. New individuals are generated, subjected to a cross mutation operation, and then merged with the parent population to obtain a new population.Step 8. Termination and Output. The fitness value of Step 5 is calculated for the individual population in Step 7, and the ranking method by Step 6 is used to reorder and stratify the population, in order to generate a new population in Step 7. This process is iteratively repeated until either the preset termination condition or the maximum number of iterations is reached. Hence, the optimal Pareto frontier is output.Figure 1
Decoding process of the chromosome encoding.


## 3.2. IHS Algorithm

The IHS multiobjective optimization method [24] is used to solve nonlinear nonconvex problems. This method generates new solution vectors to improve the accuracy and convergence speed of the harmony search algorithm (HS). The HIS-based solution steps for the optimization problem of emergency rescue multitarget material allocation are as follows:Step 1. The same as that of NSGA-II.Step 2. Initialization of the IHS Parameters. The following IHS parameters are initialized: the number of individuals NI, the number of generations NG, the probability of choosing from memory pr, the minimum pitch adjustment rate pmin, and the maximum pitch adjustment rate pmax.Step 3 to Step 5. The same as those of NSGA-II.Step 6. Improvisation of a New Harmony. Based on memory consideration, pitch adjustment, and random selection, the new harmony X′=x11′,…,xnm′ is generated as follows:(1)Memory consideration:(11)xim′⟵xim′±rd⋅bw.(2)
Pitch adjustment decision:(12)xij′⟵yeswith probabilitypminnowith probabilitypmax,j=1,2,…,m.(3)
Random selection:If the pitch adjustment decision forxij′=yes, then(13)xim′⟵xim′±rd⋅bw,wherebw is an arbitrary distance bandwidth and rd is a random floating number between 0 and 1.Step 7. Solution Update. Iterate over xi, where i = 1, 2, …, n, and exclude the worst existing harmony from the solution X.Step 8. Termination and Output. If the preset termination condition or the maximum number of iterations is reached, the algorithm is terminated. Otherwise, Steps 6 to 8 are repeated.


## 3.3. MOEA/D Algorithm

In this algorithm, a multiobjective optimization problem is decomposed into several scalar optimization subproblems, which are then simultaneously optimized. The solution of each optimization subproblem uses only information from its adjacent subproblems. This significantly reduces the computational complexity of each MOEA/D generation. The MOEA/D-based solution steps for the optimization problem of emergency rescue multitarget material allocation are as follows:Step 1. The same as that of NSGA-II.Step 2. Initialization of the MOEA/D Parameters. The following MOEA/D parameters are initialized: the number of individuals NI, the number of generations NG, the size of the weight neighborhood nd, and set the external population EP=Φ.Step 3. Distance Computation. For each i = 1, 2, …, NI, set ki=i1,…,ind, where ωi1,…,ωind are the nd closest weight vectors to ωi.Step 4. Initial Population Generation. Generate randomly an initial population X1,...,XNI, set FVi=fjXi,j=1,2, and initialize z=f1,f2T using equations (1) and (2).Step 5. Reproduction. Randomly select two serial numbers b and l from ki. Use genetic operators to generate a new solution y from Xb and Xl. Based on equations (1) and (2), apply the improvement heuristic on y to produce y′.Step 6. Update of the Objective Values z and Neighbouring Solutions. For each j=1,2,ifzj<fjy′, then set zj=fjy′. If gtey′λj,z≤gtexjλj,z, j∈Bi, then set Xj=y′ and FVj=fjy′.Step 7. External Population Update. For j = 1 and 2, remove from the external population (EP) all of the vectors dominated by fjy. Also, add fjy to EP if no vectors in EP dominate fjy. For each i = 1, 2, …, NI, loop through the execution of Step 4 to Step 7.Step 8. Termination and Output. If any of the stopping criteria are satisfied, then stop and output EP. Otherwise, go to Step 3.


## 4. Experimental Setup

In this section, we verify the validity of the multiobjective mathematical optimization model for the allocation of rescue materials. We check as well the feasibility of the solutions obtained by the NSGA-II algorithm. Givenn rescue material repositories, m types of rescue materials, an estimated rescue time T, and the constraints (5)–(9), we generated randomly the material demands, the storage capacities, the material transportation costs, and the loss due to material shortage. From this simulation, an example allocation scheme for emergency rescue materials is generated. Optimal solutions of the established optimization problem are found using numerical simulations of the NSGA-II algorithm.


### 4.1. Experimental Design

According to actual situations, the number of the relief material storage unitsn that reach the rescue site within the estimated rescue time is chosen in the interval [3, 30]. Also, the range for the number of types of rescue materials m is set to [15, 120]. The minimum rescue time T is expected to be 3 hours, while the maximum rescue time should be 72 hours. Given the total number of rescue material repositories n and the number of types of the rescue materials m, the rescue time is randomly generated, and 45 test cases are obtained. As shown in Table 2, for the number of the material storage units n, the maximum and minimum values are 30 and 5, respectively. The 25% and 75% percentile values are, respectively, 10 and 28, covering the upper and lower limits of the set. Also, the 45 test cases are more evenly distributed within the set range. For the number of recue material types m, the maximum and minimum values are 115 and 43, respectively. The 25% and 75% percentile values are, respectively, 61 and 86. Though the values are mostly distributed in [61, 86], they still cover most of the preset range. For the estimated rescue time T, the maximum and minimum values are 70 and 4, respectively. The 25% and 75% percentile values are 14 and 49, respectively. For the time variables, the maximum and minimum values are close to the preset upper and lower limits. In addition, the time in the test cases is concentrated between 14 and 49 hours.Table 2
The distribution of the variables of the storage unit countn, the material type count m, and the rescue time T in 45 test cases, where the rescue time T is in hours.
TestnmT154372561213584154586655115136104377104347810612691061341010842211108450121086151310862614101156151011549162543111725433318256114192561202025842421258441222586352325865624251151025251155826284352728436328286164292861653028842931288460322886113328866134281152335281155036304383730432438306120393061604030843241308463423086224330864444301154453011570The minimum time for a material storage unit to reach a rescue site is assumed to be 0.2 hours, while the longest time is the corresponding predicted rescue timeT under each test condition. The demand qj of each type of material is in the range of [1, 2000]. For each material type j and material repository Bi, the storage capacity Sij is in the range [100, 1000], and the unit transportation cost Ci j is in [1, 50]. As well, the loss Lj caused by the lack of the j-th material per unit time is in [0, 1]. According to the data and variable ranges given in Table 2, the required variable values corresponding to each test are randomly generated under the constraints (5)–(9). The simulation data are generated as follows:(1)
For each test case, the timeti required to send rescue materials from each material storage Bi to the rescue site A is generated randomly, while satisfying the side constraint 0.2≤ti≤T and sorting ti so that ti+1≥ti, i = 1, 2, …, n.(2)
The demandqj of the j-th material type is generated randomly. The demand is related to the emergency impact and the rescue time. Especially, the demand is positively correlated with the rescue time T through qj=⌊qj′T⌋,q′∈1,2000, j=1,2,...,m.(3)
The lossLj caused by the lack of the j-th material type in a unit time is generated randomly within the range [0, 1], j=1,2,…,m.(4)
The quantityνj of the j-th material type consumed per unit time is randomly generated. Here, the material consumption is considered as a fixed value that is related to the rescue time and the material demand by vj=qj/T, j=1,2,…,m.(5)
For each material typej and material repository Bi, the unit transportation cost Cij is generated randomly within the range [1, 50], i = 1, 2, …, n, j=1,2,…,m.(6)
For each material typej and material repository Bi, the storage capacity Sij is generated randomly, such that the total amount of all materials in the material storage units meets the rescue needs, i.e., Sij=∑i=1nsij<1.2qj,sij∈100,1000, i = 1, 2, …, n, j=1,2,…,m.


### 4.2. Evaluation Indices

#### 4.2.1. Hypervolume Index

The hypervolume index is a set measure used to evaluate the convergence and diversity performance of multiobjective optimization evolutionary algorithms [35]. In our model, there are two objective functions. So, when the number of the Pareto solutions is N, the k-th solution is represented by f1k,f2k, k = 1, 2, …, N. If we set the reference point as f1N+1,f2N+1, then the hypervolume index value is defined as(14)IHN=∑k=1Nf1k+1−f1kf2N+1−f2k.The larger this value is, the better the algorithm convergence and diversity are. Inversely, the worse the algorithm convergence and diversity are, the lower the value of the hypervolume index is.


#### 4.2.2. Diversity Metric

To evaluate the diversity of the Pareto-optimal region in multiobjective optimization, the diversity metricΔ is typically used [32]. This metric measures the extensibility and nonuniformity of the distribution of the Pareto-optimal region. The diversity metric is defined as(15)Δ=df+dl+∑i=1N−1di−d¯df+dl+N−1d¯,where N is the number of solutions on the Pareto frontier; df is the Euclidean distance of the extreme solution; dl is the Euclidean distance of the obtained boundary solution; di is the Euclidean distance between the i-th point and the (i + 1)-th point in the obtained nondominated set of solutions (i = 1, 2, …, N − 1); d¯ is the average of all distances di. A perfect solution distribution will make all the Euclidean distances equal to di and satisfy df=dl=0. Therefore, the better the distribution is, the smaller the diversity metric Δ is. Inversely, a larger diversity metric will lead to a worse solution distribution.


## 4.1. Experimental Design

According to actual situations, the number of the relief material storage unitsn that reach the rescue site within the estimated rescue time is chosen in the interval [3, 30]. Also, the range for the number of types of rescue materials m is set to [15, 120]. The minimum rescue time T is expected to be 3 hours, while the maximum rescue time should be 72 hours. Given the total number of rescue material repositories n and the number of types of the rescue materials m, the rescue time is randomly generated, and 45 test cases are obtained. As shown in Table 2, for the number of the material storage units n, the maximum and minimum values are 30 and 5, respectively. The 25% and 75% percentile values are, respectively, 10 and 28, covering the upper and lower limits of the set. Also, the 45 test cases are more evenly distributed within the set range. For the number of recue material types m, the maximum and minimum values are 115 and 43, respectively. The 25% and 75% percentile values are, respectively, 61 and 86. Though the values are mostly distributed in [61, 86], they still cover most of the preset range. For the estimated rescue time T, the maximum and minimum values are 70 and 4, respectively. The 25% and 75% percentile values are 14 and 49, respectively. For the time variables, the maximum and minimum values are close to the preset upper and lower limits. In addition, the time in the test cases is concentrated between 14 and 49 hours.Table 2
The distribution of the variables of the storage unit countn, the material type count m, and the rescue time T in 45 test cases, where the rescue time T is in hours.
TestnmT154372561213584154586655115136104377104347810612691061341010842211108450121086151310862614101156151011549162543111725433318256114192561202025842421258441222586352325865624251151025251155826284352728436328286164292861653028842931288460322886113328866134281152335281155036304383730432438306120393061604030843241308463423086224330864444301154453011570The minimum time for a material storage unit to reach a rescue site is assumed to be 0.2 hours, while the longest time is the corresponding predicted rescue timeT under each test condition. The demand qj of each type of material is in the range of [1, 2000]. For each material type j and material repository Bi, the storage capacity Sij is in the range [100, 1000], and the unit transportation cost Ci j is in [1, 50]. As well, the loss Lj caused by the lack of the j-th material per unit time is in [0, 1]. According to the data and variable ranges given in Table 2, the required variable values corresponding to each test are randomly generated under the constraints (5)–(9). The simulation data are generated as follows:(1)
For each test case, the timeti required to send rescue materials from each material storage Bi to the rescue site A is generated randomly, while satisfying the side constraint 0.2≤ti≤T and sorting ti so that ti+1≥ti, i = 1, 2, …, n.(2)
The demandqj of the j-th material type is generated randomly. The demand is related to the emergency impact and the rescue time. Especially, the demand is positively correlated with the rescue time T through qj=⌊qj′T⌋,q′∈1,2000, j=1,2,...,m.(3)
The lossLj caused by the lack of the j-th material type in a unit time is generated randomly within the range [0, 1], j=1,2,…,m.(4)
The quantityνj of the j-th material type consumed per unit time is randomly generated. Here, the material consumption is considered as a fixed value that is related to the rescue time and the material demand by vj=qj/T, j=1,2,…,m.(5)
For each material typej and material repository Bi, the unit transportation cost Cij is generated randomly within the range [1, 50], i = 1, 2, …, n, j=1,2,…,m.(6)
For each material typej and material repository Bi, the storage capacity Sij is generated randomly, such that the total amount of all materials in the material storage units meets the rescue needs, i.e., Sij=∑i=1nsij<1.2qj,sij∈100,1000, i = 1, 2, …, n, j=1,2,…,m.


## 4.2. Evaluation Indices

### 4.2.1. Hypervolume Index

The hypervolume index is a set measure used to evaluate the convergence and diversity performance of multiobjective optimization evolutionary algorithms [35]. In our model, there are two objective functions. So, when the number of the Pareto solutions is N, the k-th solution is represented by f1k,f2k, k = 1, 2, …, N. If we set the reference point as f1N+1,f2N+1, then the hypervolume index value is defined as(14)IHN=∑k=1Nf1k+1−f1kf2N+1−f2k.The larger this value is, the better the algorithm convergence and diversity are. Inversely, the worse the algorithm convergence and diversity are, the lower the value of the hypervolume index is.


### 4.2.2. Diversity Metric

To evaluate the diversity of the Pareto-optimal region in multiobjective optimization, the diversity metricΔ is typically used [32]. This metric measures the extensibility and nonuniformity of the distribution of the Pareto-optimal region. The diversity metric is defined as(15)Δ=df+dl+∑i=1N−1di−d¯df+dl+N−1d¯,where N is the number of solutions on the Pareto frontier; df is the Euclidean distance of the extreme solution; dl is the Euclidean distance of the obtained boundary solution; di is the Euclidean distance between the i-th point and the (i + 1)-th point in the obtained nondominated set of solutions (i = 1, 2, …, N − 1); d¯ is the average of all distances di. A perfect solution distribution will make all the Euclidean distances equal to di and satisfy df=dl=0. Therefore, the better the distribution is, the smaller the diversity metric Δ is. Inversely, a larger diversity metric will lead to a worse solution distribution.


## 4.2.1. Hypervolume Index

The hypervolume index is a set measure used to evaluate the convergence and diversity performance of multiobjective optimization evolutionary algorithms [35]. In our model, there are two objective functions. So, when the number of the Pareto solutions is N, the k-th solution is represented by f1k,f2k, k = 1, 2, …, N. If we set the reference point as f1N+1,f2N+1, then the hypervolume index value is defined as(14)IHN=∑k=1Nf1k+1−f1kf2N+1−f2k.The larger this value is, the better the algorithm convergence and diversity are. Inversely, the worse the algorithm convergence and diversity are, the lower the value of the hypervolume index is.


## 4.2.2. Diversity Metric

To evaluate the diversity of the Pareto-optimal region in multiobjective optimization, the diversity metricΔ is typically used [32]. This metric measures the extensibility and nonuniformity of the distribution of the Pareto-optimal region. The diversity metric is defined as(15)Δ=df+dl+∑i=1N−1di−d¯df+dl+N−1d¯,where N is the number of solutions on the Pareto frontier; df is the Euclidean distance of the extreme solution; dl is the Euclidean distance of the obtained boundary solution; di is the Euclidean distance between the i-th point and the (i + 1)-th point in the obtained nondominated set of solutions (i = 1, 2, …, N − 1); d¯ is the average of all distances di. A perfect solution distribution will make all the Euclidean distances equal to di and satisfy df=dl=0. Therefore, the better the distribution is, the smaller the diversity metric Δ is. Inversely, a larger diversity metric will lead to a worse solution distribution.


## 5. Results and Discussion

In order to verify the feasibility of the proposed mathematical model of continuous material allocation as well as the effectiveness of the algorithmic solution of the multiobjective optimization problem, the improved harmony search algorithm (IHS) [36], the MOEA/D algorithm [25], and NSGA-II were used to compare and analyze the solutions for 45 test cases, where the iteration count is set to NG=n×m×4 and the population size NI=256. For the three aforementioned algorithms, the average results of the hypervolume index, the diversity metric, and the fitness values are shown in Figures 2–4, respectively. The parameters of each group of test solutions are consistent. The parameters of the IHS algorithm are set as follows: the probability of choosing from the memory pr=0.85, the minimum pitch adjustment rate pmin=0.35, and the maximum pitch adjustment rate pmax=0.99. The parameters of the MOEA/D algorithm are set as follows: the size of the weight neighborhood nd=20, the crossover parameter in the differential evolution operator Cr=1, the parameter for the differential evolution operator F=0.5, and the distribution index used by the polynomial mutation em=20. The NSGA-II parameters are set as follows: the crossover probability cr=0.95, the crossover distribution index Dcr=10, the mutation probability mr=0.01, and the mutation distribution index Dmr=50.Figure 2
Comparisons of the average hypervolume index values.Figure 3
Comparisons of the diversity metric values.Figure 4
Comparisons of fitness values for the two types of objective functions: (a) fitness for the cost objectivef1 and (b) fitness for the loss objective f2.
(a)(b)Using the IHS algorithm, we found the mean and variance of the hypervolume index to be 1.93 and 0.099, respectively. Based on the MOEA/D algorithm, the mean and the variance of the hypervolume index are 2.18 and 0.137, respectively. For the NSGA-II method, the corresponding mean and variance of the hypervolume index are 2.42 and 0.003, respectively. The evaluation of the hypervolume index shows that the best algorithmic solution of the continuous material allocation problem is obtained by NSGA-II, followed by the MOEA/D algorithm and then the IHS one (see Figure2).Using the IHS method, the mean and variance of the diversity metric values are 1.7 and 0.185, respectively. Based on the MOEA/D method, the corresponding mean and variance of the diversity metric are 1.51 and 0.112, respectively. For the NSGA-II method, we obtained the mean and the variance of the diversity metric as 1.23 and 0.109, respectively. The evaluation of the diversity metric shows the same relative performance of the three algorithms, i.e., the NSGA-II method outperforms the MOEA/D one, which is in turn better than the IHS algorithm (see Figure3).Using the IHS method, we found the fitness mean and variance of the objectivef1 to be 1.08×107 and 3.82×1014, respectively. The corresponding values for the objective f2 are 8.3×103 and 7.01×108, respectively. Using the MOEA/D method, the fitness average and variance of the objective f1 were found to be 5.7×106 and 4.62×1013, respectively. The corresponding values for f2 are 9.71×103 and 6.75×108, respectively. Using the NSGA-II method, we computed the fitness mean and variance of f1 to be 1.37×106 and 9.35×1011, respectively. The corresponding values for f2 are 2.36×103 and 7.43×106, respectively. The results in Figure 4(b) show fitness values below 15 using the NSGA-II method, while the IHS method got three larger fitness values in the cases n=10,43,47, m=43,84,115, and T=47,50,49, respectively. Also, the MOEA/D method resulted in two larger fitness values in the cases n=43,47, m=84,115, and T=50,49, respectively. In the case n=43, m=84, and T=50, the methods of NSGA-II, IHS, and MOEA/D got the largest fitness values of 12.7, 166.4, and 166.4, respectively. Due to the limitations on the number of materials and the quantity of each material, the disaster-induced loss is relatively large no matter what algorithm is used. Nevertheless, IHS and MOEA/D are more likely to fall into local optima than NSGA-II during the optimization process. Also, it is difficult to minimize the loss in the shortest time. The results show a consistent ranking of the three algorithms where the NSGA-II method gives the best objective value, followed by the MOEA/D method and then the IHS one (see Figure 4).In general, for the given parameters and among the methods considered for solving the multiobjective continuous material allocation problem, the NSGA-II method has the best performance in terms of convergence, diversity, and extensiveness, while the IHS method has the worst performance.By using the NSGA-II method, the convergence and distribution hypervolume measure of each group was calculated. For the 45 test cases, the average and variance of the hypervolume index were 2.42 and 0.003, respectively. The average and variance of the hypervolume index for each case are shown in Table3.Table 3
Hypervolume index results for 45 test cases.
TestMeanVariance12.3192.85E − 0522.4722.21E − 0732.4061.21E − 0442.4871.02E − 0752.4515.22E − 0662.4203.85E − 0672.4151.31E − 0682.2742.12E − 0392.4587.22E − 07102.4315.53E − 06112.4671.06E − 07122.3096.83E − 06132.4882.64E − 07142.4632.51E − 06152.4741.61E − 07162.4661.00E − 06172.4468.94E − 07182.4013.39E − 06192.4022.82E − 06202.3842.52E − 06212.4281.66E − 06222.3132.90E − 06232.4941.33E − 08242.4273.05E − 06252.4067.46E − 06262.4054.60E − 07272.4367.07E − 07282.3824.85E − 06292.4681.06E − 07302.3632.77E − 05312.4711.67E − 07322.2711.58E − 04332.4703.46E − 07342.4924.50E − 08352.4712.31E − 06362.4232.16E − 06372.4366.81E − 07382.3692.37E − 06392.4426.37E − 06402.4021.73E − 06412.4363.72E − 06422.4838.50E − 08432.4327.59E − 06442.4539.56E − 07452.4762.40E − 07The average, minimum, maximum variance, and minimum variance of the hypervolume index are 2.494, 2.271, 2.12E − 03, and 1.33E − 08, respectively.Using the NSGA-II method, we calculated the diversity metricΔ of each test case. For the 45 test cases, the average and variance of the diversity metric were 1.23 and 0.109, respectively. The average results for df,dl, and Δ for multiple runs are shown in Table 4.Table 4
Results ofdf, dl, and Δ for 45 test cases.
TestdfdlΔ1001.66420.02601.1873001.96140.00901.39350.00201.3516001.44270.00400.9618001.6359001.0081000.1330.73411001.37612001.2231300.0210.797140.00101.27615001.036160.00100.698170.0010.0021.3551800.0060.88819001.54320001.479210.0040.0170.85422001.105230.00201.1224001.46725001.6022600.0151.22427001.29428001.529001.72830001.59731001.614320.0010.0231.10633001.2953400.0011.327350.010.0041.2343600.0041.3383700.1551.26338001.245390.0020.0131.16640001.604410.00201.769420.00101.11943001.391440.00200.972450.0140.0231.017The maximum values ofdf and dl are 0.281 and 0.155, respectively. The average difference between df and dl is 0.013, while the maximum absolute difference is 0.155. Test cases with a difference of less than 0.0009 accounted for approximately 57% of the total test cases, while test cases with a difference of less than 0.01 accounted for approximately 80% of the total number of cases.Analysis of the distribution of the Pareto dominant set shows that this set has sufficient extensibility and uniform distribution for all test cases. The measurement magnitude has nothing to do with the number of rescue material storage unitsn, the types of rescue materials m, and the estimated time T. That demonstrates the rationality of the material allocation model and the feasibility of the solution method.Four test cases were selected to analyze the distribution of the Pareto solution set in a single run. The Pareto frontier results are shown in Figure5. The abscissa is the fitness value of the objective function f1 in (1), and the ordinate represents the fitness value of the objective function f2 in (2). In order to exhibit particular details legibly and clearly, the target values are reduced by 102–106. As shown in Figure 5, the Pareto noninferior solution is represented by a dotted curve, and the Pareto frontier is represented by a broken line.Figure 5
The distribution of the Pareto-optimal solution set for four test cases of the material allocation problem.
(a)(b)(c)(d)According to the objective functions (1) and (2), the emergency rescue material allocation is a multiobjective nonlinear optimization problem. The Pareto frontier distribution demonstrates the validity and reliability of the constraint processing method. The mathematical model described in this paper is a multiobjective optimization model with sufficient materials. Under this condition, the Pareto frontier results show that when the number of rescue storage units n and the number of rescue material types m are fixed, the longer the estimated rescue time T is and the greater the emergency rescue cost and the resulting losses are. When the number of the emergency rescue storage units n and the estimated rescue time T are ascertained, the more the number of emergency rescue material types m is and the more the emergency rescue costs and caused losses are. When the number of the emergency rescue material types m and the estimated rescue time T are constant, the more the number of the emergency rescue storage units n is and the smaller the emergency rescue cost and the resulting losses are. Therefore, in order to minimize the economic costs and rescue losses, we should establish a sufficient number of emergency rescue material repositories and ensure that there are enough material types and quantities. Any point on the Pareto frontier distribution is a noninferior solution of the multiobjective optimization problem, and there is no distinction for such solutions. There is a contradiction between the objective values f1 and f2, which stems from the fact that the smaller the emergency rescue cost is, the greater the loss is and the higher the emergency rescue cost is. Therefore, in the emergency rescue operations, decision makers should consider the actual situation, especially when the lack of materials has a little impact on the damage caused by the rescue operations. When decision makers seek to reduce the economic cost of the emergency rescue operations, the satisfactory emergency rescue solution should be selected from the lower right corner in Figure 5. When the economic cost of the emergency rescue has a little impact, then the lack of materials will cause little rescue damage, and the satisfactory solution should be obtained from the upper left corner in Figure 5. When the economic cost of the emergency rescue has the same impact as the loss of materials, a compromise could be made, and the satisfactory solution should be obtained from the middle part in Figure 5.A noninferior solution of the bottom right corner of the Pareto frontier is taken as an example to analyze the proposed scheme. In this example, there are 30 material storage warehouses that can provide materials of 86 types to the rescue site within the estimated rescue time of 22 hours (n=30,m=86,T=22). The time distribution after sorting is shown in Figure 6. The abscissa shows the number of each emergency rescue material repository after sorting the times to reach the rescue site. The ordinate is the time for each repository to arrive at the rescue site. The shortest time is 0.4 hours, and the longest time is 20 hours.Figure 6
Response times to reach the rescue site for the emergency relief material storage units.The number of material types needed is 86. The distribution of the loss caused by the lack of each material type per unit time is shown in Figure7, where the abscissa represents the material types, the left ordinate represents the demand of each material type (the black solid line), and the right ordinate represents the loss caused by the lack of each material type (the red dotted line). The least demand is that of the 80th material m80, for which 88 units of supplies are needed. The highest demand is that of the 49th material m49, for which 968 units of materials are required. The loss caused by the lack of materials hits a value of 1 for the 25th material m25 and the 56th material m56. The minimum loss of 0.003 is attained for the 77th material.Figure 7
Demand and loss for various types of emergency rescue materials.The NSGA-II method is used to solve the established multiobjective optimization problem. Under the above conditions, 234 noninferior solutions are obtained by the NSGA-II method. Different test cases may lead to equal objective function values. So, 255 alternative material allocation schemes are obtained. The material allocation scheme corresponding to the selected noninferior solution is shown in Figure8. Table 5 is a partial result corresponding to the noninferior material allocation scheme.Figure 8
The material allocation scheme corresponding to the noninferior solution atn=30,m=86,T=22.Table 5
Partial results of the material allocation scheme corresponding to the noninferior solution atn=30,m=86,T=22.
A1A2A3A4…A27A28A29A30m11309269…0184m2731771…106213m3124491811…241513425m417114222230139…………………………m8317002243…5472m84160361033…295452m85172745126…17425653m86232283…11401In Figure8, the abscissa indicates the number of 30 emergency relief material repositories sorted by time, while the ordinate indicates 86 rescue material types. The quantity of the emergency rescue material is expressed by an indicative color, which has a few variations from dark to light settings. The darker the color is, the more the allocated materials are. Inversely, the lighter the color is, the less the allocated materials are. Because the value of the objective function f2 in (2) is related to the time when each material storage unit reaches the rescue site and the material storage units are sorted by that time, the allocation of various material types corresponding to B1 in Figure 7 is the most appropriate. According to the demand size (the black solid line in Figure 6) and the loss caused by the lack of unit materials per unit time (the red dotted line in Figure 6), not all materials are dispatched in the first material storage unit B1 to ensure that the balance between the two objective values is achieved. From the material allocation scheme in Figure 7, it can be confirmed that if the demand for the 40th material type m40 is q40 = 660, the loss caused by the lack of this material per unit time is L40 = 0.06, and its allocation amounts in B1 and B2 are 21 and 258, respectively. Therefore, these results demonstrate that the emergency relief material allocation scheme meets the requirements of the established mathematical optimization model.


## 6. Conclusions

Observing the constraints of continuous supply and demand of time-varying sequences, we establish a multiobjective mathematical optimization model for the continuous allocation of emergency relief materials. We consider the dynamic characteristics of disasters and accidents and the demand continuity for emergency relief materials, in order to reduce the losses and the economic cost of the emergency relief operations. By comparing the convergence, distribution, and diversity of the solutions based on the IHS, MOEA/D, and NSGA-II methods of multiobjective optimization, we conclude that the Pareto solution set can be obtained by the NSGA-II method with fairness and balance of loss and cost. The solution set has the advantages of higher scalability, uniformity of distribution, and convergence. According to the results of the loss caused by the lack of rescue materials and the impact of the economic costs on the rescue work, the decision makers can select the appropriate noninferior material allocation scheme from the frontier of the Pareto set. Forty five sets of emergency relief material allocation examples were simulated and tested. It is found that the NSGA-II method is more effective than its competing algorithms for solving the multiobjective optimization problem of continuous rescue material allocation. For the NSGA-II method, the mean and variance of the hypervolume index are 2.42 and 0.003, respectively. The corresponding values for the diversity metric are 1.23 and 0.109, respectively. These results verify the validity of the mathematical model. Through the analysis of the noninferior solutions, the results show that the emergency rescue material allocation scheme agrees with the established mathematical optimization model. The number of storage units of the emergency rescue materials, the types of emergency rescue materials, and the estimated rescue time are independent of the solution results, an observation that demonstrates the feasibility of using the NSGA-II method. Thus, the multiobjective mathematical optimization problem of continuous material allocation for emergency relief can be solved by the NSGA-II method, and the results can provide decision support for the continuous allocation and planning of emergency relief materials.



---

*Source: 5693182-2020-07-12.xml*