Abstract
Thousands of victims and millions of affected people are hurt by natural disasters every year. Therefore, it is essential to prepare proper response programs that consider early activities of disaster management. In this paper, a multiobjective model for distribution centers which are located and allocated periodically to the damaged areas in order to distribute relief commodities is offered. The main objectives of this model are minimizing the total costs and maximizing the least rate of the satisfaction in the sense of being fair while distributing the items. The model simultaneously determines the location of relief distribution centers and the allocation of affected areas to relief distribution centers. Furthermore, an efficient solution approach based on genetic algorithm has been developed in order to solve the proposed mathematical model. The results of genetic algorithm are compared with the results provided by simulated annealing algorithm and LINGO software. The computational results show that the proposed genetic algorithm provides relatively good solutions in a reasonable time.
1. Introduction
It is inevitable to have an integrated scientific system for crisis’ logistics management with clearly defined functions and duties. Optimizing logistics and relief chain can affect issues such as relief logistics management [1]. Logistics can create greater coordination for delivering relief commodities while accelerating the delivery speed and reducing the response time [2, 3]. One of the logistics subproblems with strategic roles in a critical situation is to locate the depots and their supply savings [2]. In the depot locating, factors such as cost, response time, and justice in distribution have to be considered in preparedness phase [4]. Although researches in facilities location problems have been widely done both theoretically and practically, the problems have not been used in humanitarian logistics [2].
The application of operation research in disaster management programs has been one of the main issues in recent decades. Supply chain and logistics management are recently being used as analytical tools and techniques to provide efficient and effective relief to people who need help in devastated areas with optimized functions and activities. Thus, emergency logistics management has been emerged as a worldwide noticeable theme as disasters may occur anytime around the world with enormous consequences. The aim of response to disasters in humanitarian relief chains is to provide a quick relief to the affected areas in order to minimize the death and pain of people. Proper design and operation of a relief chain are an essential element to achieve an effective and efficient response, but, only in recent years, humanitarian organizations have paid attention to its importance [5].
In urban areas, earthquakes can cause serious damage to inhabitance and force people to leave their homes. Therefore, at the preparation planning phase of disaster management, municipal authorities should consider emergency shelters, proper equipment, and supplies for the affected people in order to reduce the number of casualties and people’s sufferings and bring relief to survivors [5].
The number and location of the facilities and distribution centers decrease the response time and the operation costs in relief chain. Therefore, Location-Allocation is an important issue in this process. Some of the applications of these locations in the relief logistics problem are listed as follows:(i)the location of critical management centers in the city or province and the allocation of emergency commodities such as water, food packages, and tents from these centers to damaged areas,(ii)Red Crescent depots’ locating and these depots allocation to damaged cities and provinces in the situation of facing natural crisis or even war,(iii)locating ambulance centers inside and between the cities and provinces.
Actually in crisis management, we deal with locating process already in preparedness and response phase. Locating process and facility capacity determination are the two main topics in the disaster management, and a few researches have been focusing on this subject [6].
Barbarosoǧlu and Arda [7] modeled the transportation planning in earthquake response by a scenario-based two stage stochastic programming over a multicommodity and multimodal network flow with uncertainty in demand, supply, and route capacities. The inclusion of uncertainties is an important advance in the analysis, but the focus is still postdisaster response. They did not consider the location of facilities or the inventory decisions in their model. Akkihal [8] has proposed a model for crisis management-center location in order to manage relief items. Tzeng et al. [9] have presented a certain multicriteria model in order to distribute relief commodities to damaged areas considering costs, response time, and damaged people satisfaction and solved it by fuzzy multiobjective programming. Jia et al. [10] presented an incapacitated facility location model to locate emergency service facilities in the event of a large-scale emergency.
Balcik and Beamon [2] considered the facility location problem for humanitarian relief chains in order to respond to quick onset disasters. Their proposed model aimed to locate and determine the situation and the number of distribution centers in relief network and the number of stored commodities in order to meet the demand of affected people.
Beraldi and Bruni [11] also offered a probable model to have a relief facility to optimize locating in an uncertain environment. Mete and Zabinsky [12] introduced a stochastic optimization model for disaster preparedness and response under demand/cost uncertainty in order to assist deciding on the location and allocation of medical supplies which are used during emergencies. They also offered a mixed-integer programming transportation model that is potentially useful in routing decisions during the response phase. Rawls and Turnquist [13] developed a two-stage, stochastic, mixed-integer program that determined the locations and quantities of various types of emergency commodities; their model also contained the transportation network availability following a disaster under demand/cost uncertainty.
Moreover, Jabal-Ameli et al. [14] presented a location and distribution model and solved it by fuzzy methods. Duran et al. [15] studied inventory and location model and offered a mathematical model which investigates the effect of location on the average response time. Zeng and Wu [16] considered location and routing in relief situation in two steps. It means that they first solved the location problem and directly solved the routing problem in a way to minimize total costs and then they utilized a two-step heuristic algorithm to solve it.
Bozorgi-Amiri et al. [17] considered an uncertain model for relief logistics. They suggested that some prepositioned relief items might be destroyed after disasters based on this idea; they developed a robust multiobjective approach to their uncertain location problem and solved it for a real case study in Iran. Bozorgi-Amiri et al. [18] investigated uncertainty in many parameters of a relief operation like demand, supply, and operational costs associated with it. Location of relief centers and allocation of affected areas to these centers can be determined in situation described in their model. Yazdian and Shahanaghi [19] presented a multiobjective possibility programming approach for locating distribution centers and allocating customers’ demands in supply chains.
In 2011, Barzinpour and Esmaeili [5] developed a multiobjective MILP model for a two-echelon relief chain in order to maximize coverage of urban populations and minimize logistics costs. Tancrez et al. [20] considered inventory-location problem for three-level supply chain in which distribution decisions of distribution centers, allocation, and number of the transported commodities were made altogether.
In 2012, Schmid [21] offered a probability model for relocation problem and dynamic distribution and solved it by dynamic programming. Chrétienne et al. [22] suggested location and dispatching problem which included minimizing setup and availability costs and then solved it through branch and cut algorithm. Toro-Diaz et al. [23] offered a mathematical model for location and distribution problem and solved it through genetic algorithm. In 2012, Guillermo Cabrera et al. [24] presented a hybrid approach using an artificial bee algorithm with mixed-integer programming applied to a large-scale capacitated facility location problem. Murali et al. [25] considered a facility location problem to determine the points in a city where medicine should be handed out to the population. They consider locating capacitated facilities in order to maximize the coverage, considering a distance-dependent coverage function.
In 2013, Xi et al. [26] developed a modified -median problem model that accounts for rescue time limitations. A variable neighborhood search- (VNS-) based algorithm is developed for the model considered. Gharegozloo Hamedani et al. [27] studied a multiobjective location problem in a three-level supply chain network under uncertain environment, considering inventory decisions. The proposed model presents a robust optimization model, which specifies locations of distribution centers to be opened, inventory control parameters, and allocation of supply chain components, concurrently.
Abounacer et al. [28] proposed a three-objective location-transportation problem for disaster response. The location problem aims at determining the number, the position, and the mission of required humanitarian aid distribution centers (HADC) within the disaster region. Three conflicting objectives are considered. The first objective minimizes the total transportation duration of needed products from the distribution centers to the demand points. The second objective minimizes the number of agents (first-aiders) required to open and operate the selected distribution centers. The third objective minimizes the noncovered demand for all demand points within the affected area. They proposed an epsilon-constraint method for this problem. Esmaeili and Barzinpour [29] based on a real world case study for a municipal district in Tehran developed a multiobjective mathematical model for the location-distribution problem.
Most of the researches in location and distribution of relief logistics are dependent on preparedness phase on the basis of predicting programming conditions. Now, the question is whether the same criterion in preparedness phase can be operation criterion or, based on the new conditions and available information, new programs should be suggested. Usually after any casualty, in the first 72 hours, based on the available commodities in the distribution centers, services are offered. Therefore, it would be possible to program for damaged areas during that time based on the received information. In this paper, the goal is to present a model to locate distribution centers periodically and to consider good distribution manner toward damaged areas.
This paper focuses on the logistics aspect of the response phase and more precisely on two important related problems: location and allocation. The location problem aims at designing a network for distributing humanitarian aids (e.g., water, food, medical commodities, and survival equipment). It mainly consists of determining the position and the mission of required humanitarian aid distribution centers (HADC) within the disaster region. The allocation problem deals with the distribution of humanitarian aid from HADCs to demand points. When both problems are solved simultaneously, we discuss a location-allocation problem. Based on this particular context, this paper has provided a biobjective location-allocation problem. The first objective is to minimize the total costs. The second objective is to maximize the covered demand for all demand points.
We propose a genetic algorithm and a simulated annealing algorithm to biobjective location–allocation problem addressed.
2. Problem Definition and Model Assumptions
A two-level system including distributers and damaged areas is considered. There are regions as candidate regions to establish distribution centers in this system. And also the number of commodities in distribution centers and demands of damaged areas in different periods are available. Our goal is to determine active distribution centers and their distribution manner in different periods. Location and allocation have to be done so that total costs become minimized and commodities have to be distributed in a fair manner.
2.1. Model Assumptions
(i)The inventories of distribution centers in different periods are available in the form of triangular fuzzy numbers.(ii)The damaged area requests in different periods are available in the form of triangular fuzzy numbers.(iii)The maximum number of distribution centers that can be active in each period is determined.(iv)Some types of relief commodities are available.(v)It is possible to transfer inventories from one period to the next one.(vi)Unsatisfied request does not transfer from one period to the next one.
2.2. Model Indexes
: set of distribution centers in period , : set of damaged areas in period , : set of relief commodities, : indexes related to distribution centers, : indexes related to damaged area, : indexes related to relief commodities.
2.3. Model Parameters
: setup cost of distribution center in period that is a triangular fuzzy number, : the maximum number of distribution centers that can be active in period , : transporting cost of commodity unit between distribution center and damaged area in period that is a triangular fuzzy number, : transporting cost of commodity unit between distribution center and distribution center in period that is a triangular fuzzy number, : the number of commodities in distribution center in period that is a triangular fuzzy number, : the number demand commodities in damaged area in period that is a triangular fuzzy number, : penalty costs for the shortage of good , : a big number.
2.4. Decision Variables
: number of commodities that can be delivered in period from distribution center to damaged area , : binary variable, it is 1 if commodity transfers from distribution center to distribution center in period ; otherwise, it is 0, : binary variable, it is 1 if distribution center becomes active in period ; otherwise, it is 0, : amount of transferred supply of commodity in period to period in distribution center .
2.5. The Proposed Mathematical Model
Consider There are two objective functions in the above model.
Objective 1. It aims at minimizing total transportation costs between distribution centers and damaged areas and also minimizing setup costs and penalty costs for the shortage of goods.
Objective 2. It aims at maximizing damaged area satisfaction through maximizing satisfied requests of damaged areas.
Constraint (3) indicates the number of distribution centers that can be active. Constraint (4) shows that it will be possible to send goods from one distribution center to a damaged area if that center becomes active in that period. Constraint (5) shows the way in which it will be possible to calculate transported inventories from one period to the next period. Constraint (6) indicates that the number of received commodities to one damaged area should not be more than its demand. Constraint (7) shows that it will be possible to send goods from one distribution center to another distribution center if that destination center becomes active in that period. Constraints (8)–(11) express the nature of decision variables used in the model.
2.6. Linearization of the Proposed Model
The second objective function is nonlinear term. To make it linear, we act as follows:
In constraint (5), is nonlinear. So, for linearization, it can be replaced with and three constraints are added to models as follows:
2.7. Multipurpose Linear Programming Approach
Considering the following multiobjective model, using Zimmerman method, this problem can be changed into a single objective linear programming model using max-min operator [14, 30]:
Membership functions for objective functions are defined as follows:
In the above equations, , are positive ideal solutions and , are negative ideal solutions. Now, by max-min operator and satisfied degree, the multiobjective linear programming problem has been written as a single-objective programming problem
2.8. Solving Proposed Model
The first objective function of the proposed model is a function with triangular fuzzy coefficient. Therefore, we can change it to 3 objective functions as follows:
There are three fuzzy constraints in the presented model. Therefore, we can change it to 3 certain constraints as follows [10]:
Now, a certain linear multiobjective model can be achieved as follows:
To solve this linear multiobjective model through Zimmermann fuzzy sets (15), membership functions are defined as follows for the objective functions of the problem [30]:
To calculate and , a procedure the same as that of and would be followed respectively. Now, we use the following equations to obtain an ideal positive value (pis) and a negative value (nis):
Now, we change the multiobjective problem (19) to the single-objective problem (22) via (16). The model is defined as follows:
3. Genetic Algorithm
In the last decade, we have seen a growing interest in biologically motivated approaches such as evolutionary strategies and genetic algorithms (GAs) being applied to many complex optimization problems. The processes occurring in the natural systems have inspired the development of these algorithms. Evolutionary strategies and genetic algorithms are constructed based on the observation of evolutionary processes in biological systems. Evolutionary processes such as adaptation, selection, reproduction, mutation, and competition are closely studied and translated into the form of computer simulations. Although these algorithms are a crude simplification of the natural processes, they have been successfully applied to many complex problems that were once intractable [31, 32].
GA is a stochastic search technique that explores the problem domain by maintaining a population of individuals, which represents a set of potential solutions in the search space. GA attempts to combine the good features found in each individual using a structured yet randomized information exchange in order to construct individuals who are better suited to their environment than the individuals that were created through the evolution of better individuals; it is anticipated that the desired solution will be found. The following steps are defined to design this algorithm.
3.1. Solution Representation
To show the chromosome, we used three matrices; the first matrix amount is zero or one and the rest of them are positive integer numbers.
The first matrix is called location matrix that is related to distribution centers. Its lines are related to a time period and its columns are related to distribution centers. Its dimension is in which lines are periods and columns are distribution centers. This matrix’s elements are zero and one in a way that zero and one correspond to inactiveness and activeness of the distribution center, respectively,
The second matrix is called Commodity_Allocated matrix. To each commodity, there is a combined matrix including submatrix with dimension. Each element in these submatrices indicates the number of commodities which are sent to corresponding damaged areas from the distribution center
The third matrix is called inventory matrix that indicates the amount of transported inventory from one period to the next one in each distribution center
3.2. Generating Initial Population
After determining a technique to assign a chromosome to each solution, one can create an initial population. In this process, the location matrix is generated randomly.
To initialize the Commodity_Allocated matrices, we use a heuristic method so that first, for each commodity type, the total damaged area requests and total supply of active distribution centers have been evaluated, and, then, considering justice in distribution, commodities are distributed among damaged areas equally. Now, to allocate active centers to damaged areas, an area has been chosen randomly and it would be allocated to an active distribution center obtained by the roulette wheel, which is generated based on the allocation costs. If this distribution center cannot satisfy all predicted demands, the next center obtained by the roulette wheel will meet the remaining demands of this damaged area. Then, one area is chosen randomly from the other damaged areas. All the above-mentioned activities are executed until allocation manner of centers to damaged areas is determined for all the areas.
3.3. Sampling Mechanisms
This mechanism is related to selection chromosomes’ method. Tournament selection has been used in this research. Based on this method, first, the value of each chromosome’s fitness functions in a population has been calculated, and a tournament size is chosen (here, tournament size is 3). Then, chromosome from one population has been chosen randomly equal with tournament size and, at the end, the best chromosome based on suitability will be transferred to mating pool. This has been repeated up to the number of population. It is possible to create several copies from a chromosome with higher fitness function. This shows that better chromosomes have greater chance to be chosen.
3.4. Genetic Operators
To generate a new generation of the present chromosome, the following genetic operators could be used.
3.4.1. Crossover Operator
In this crossover operator, a row of location matrix of one chromosome is exchanged with the same row of location matrix of another one and also corresponding elements are exchanged in Commodity_Allocated matrices and inventory matrices. In other words, we choose a random number from ( is the number of time periods), and we exchange corresponding row of location matrix of parent 1 with the same row of location matrix of parent 2 and also corresponding elements are exchanged in Commodity_Allocated matrices and inventory matrices in two chromosomes. Figure 1 shows the crossover operator.
(a)
(b)
(c)
(d)
3.4.2. Mutation Operators
Mutation operator operates on location matrix in the developed algorithm. In the first mutation operator, elements corresponding to two distribution centers exchange in a row of location matrix. As distribution center locating conditions change, Commodity_Allocated matrices are valued by the stated heuristic method in Section 3.2.
Figure 2 shows the mutation operator.
(a)
(b)
In the second mutation operator, one element corresponding to the distribution center in location matrix has been considered. It will be zero if it is one and it will be one if it is zero (provided that the number of active centers does not violate the number of authorized ones).
Now, as the location conditions of distribution center change, Commodity_Allocated matrices are updated by explained heuristic method in Section 3.2.
3.5. Elitism Selection
When genetic operators (mutation and crossover) are used, it would be possible to lose the best chromosomes. Elitism is a method to keep a copy from the best chromosomes in the new generations. The above mechanism makes genetic algorithm keep some of the best solutions in each generation. Previous experiences have proved that this mechanism optimizes genetic algorithm operation and shortens coverage time.
3.6. Algorithm Stop Criterion
Since metaheuristic algorithms do not have any criteria for the global optima, a maximum number of iterations have been considered in order to stop the algorithm. Figure 3 shows the genetic algorithm flowchart.
4. Simulated Annealing Algorithm
Simulated annealing is a probabilistic metaheuristic algorithm which is a local search method. The simulated annealing begins its search from a random initial solution. The iteration loop that characterizes the main procedure randomly generates in each iteration only one neighbor of the current solution . The variation for the value of the objective function is tested for each neighbor generation. To test this variation, is computed. If the value of is less than 0 (zero), then the new solution is automatically accepted to replace . Otherwise, accepting the new solution will depend on the probability established by the Metropolis criteria, which is given by , where is a temperature parameter, a key variable for the method.
Figure 4 shows the simulated annealing flowchart. In Figure 4, n is the number of internal circles, is the number of external circles, and is the objective function.
In the following, we discuss the proposed SA heuristic in detail, including the solution representation, the generation of the initial solution, and various types of neighborhood.
Solution representation: solution representation is similar to Section 3.1.
Generation of the initial solution: this is similar to Section 3.2.
Neighborhood generation mechanism: we produce neighborhood by mutation operators that is exhibited in Section 3.4.2.
5. Numerical Results
In this section, the developed algorithms have been used on the sample problems that have been generated randomly and the presented genetic algorithm efficacy has been studied. The characteristics of instances are shown in Table 3.
To solve the model used of LINGO 11 on a computer with the following characteristics: intel(R) core(Tm) i3 [email protected] Ghz, 4 G RAM. Considering the intense effect of parameter configuration on the performance of algorithms, the Taguchi method is used for the configuration of parameters. Before calibration of the employed algorithms, we run some preliminary tests to find appropriate parameter levels. To obtain more accurate, as well as better sustained, results for the offered algorithm, the following four parameters were configured: , , , and for simulated annealing and , , , and for the genetic algorithm. These parameters and their levels are given in Tables 1 and 2. The square matrix with 4 parameters in 3 levels used in the Taguchi method is L9, which is given in Table 3. In Taguchi method, the variation of the output results is measured by means of signal-to-noise () ratio. Here, the larger value of ratio leads to the smaller variation of the response variable. Factor levels that maximize the appropriate ratio are optimal [33]. The value of ratio is computed by where is the number of repetitions under the same experimental conditions and denotes the experimental result at each repetition.
Characteristics of sample problems that have been used for Taguchi method are shown in Table 5. Calculating all test results using the Taguchi method, the mean rate of for small-scale problem is presented in Figure 5 and the mean of the objective function for small-scale problem is shown in Figure 6.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Optimum factors for genetic algorithm and simulated annealing are given in Table 6. The final configurations of parameters for all algorithms are summarized in Tables 7 and 8.
Table 9 shows that the proposed genetic algorithm and simulated annealing algorithm operate faster than LINGO with little errors and the standard deviation is not considerable in ten iterations indicating that the solutions do not have deviation. Now, we generate a number of sample problems with small and large size that are exhibited in Table 4 and then operate the proposed genetic algorithm and simulated annealing algorithm on these problems. Table 9 shows more details on computational results for all test problems. The CPU times corresponding to the exact solution, as seen in Table 9, show that computational time grows exponentially by size of the problem. Unfortunately, for the large-scale calculations, LINGO cannot even find the feasible solution in the limited run time.
In Table 9, we present the following: average value, standard deviation, best solution, computational time, PRAS%, and PRBS% for each problem. As observed in Table 9, there are only 2 cases out of 25 for which the SA algorithm performs better than the GA. Also, in 2 cases out of 25, the SA algorithm and GA behave identically. Moreover, in the other cases than those mentioned, GA performs better than the SA. Finally, in an average sense, the genetic algorithm gives less standard deviation than the simulated annealing. These observations show that the GA gives more stable results than the SA. Solution average in genetic algorithm meets more suitable conditions especially in the case of large size problems. As can be seen, the percent reduction best solution (%PRBS) and the percent reduction average solution (%PRAS) are both positive values, so the GA performs better than the SA.
The percent reduction best solution (%PRBS) and the percent reduction average solution (%PRAS) are calculated using
For the proposed algorithms, the average results of the ten simulations are also presented. The efficiency of the two algorithms is measured by the quality of the produced solutions. The quality is given in terms of the relative deviation from the optimum solution that the best solution’s error = (LINGO solution − the best solution)/LINGO solution, where LINGO solution denotes the solution of the LINGO software and the best solution is considered to be the solution of the metaheuristic algorithms. As it can be seen from Table 9, the algorithms lead to very good and stable results for most cases. Also, the genetic algorithm performs better than the simulated annealing, in each and every instance.
The comparison of LINGO with the proposed genetic algorithm shows that the GA can find approximately an optimal solution in a shorter time compared to the LINGO as presented in Table 9. In Table 9, it is also shown that the average gap between the optimal and the GA solutions is 1.421%. The maximum gap is related to the test problem .
As observed in Table 9, the standard deviation has increased as the problem’s dimension increases. To obtain better solution in large dimensions, population size and the number of generations should both increase.
A small size example is considered to test the reliability of the model. Figure 7 shows the map of south-khorasan province of Iran that has been used at this example. Three cities (Birjand, Qaen, and Sarayan) are selected to set up relief commodities distribution centers for different periods. Furthermore, nine cities including Birjand, Qaen, Darmiyan, Nehbandan, Ferdows, Zirkoh, Sarayan, Sarbisheh, and Boshruyeh are considered as damaged areas. The fixed cost to open distribution centers for all the areas is $25000, $26000, and $27500. The shortage penalty for each commodity is $1.5 for commodity 1 and $2.5 for commodity 2. These values are confirmed with experts’ estimations. In this example, two types of relief commodities have been considered. In addition, there are two active centers for each period. The model has been solved for a situation that the number of planning periods is 3. The demand of damaged areas is available as triangular fuzzy numbers in Table 11. The amount of inventory in distribution centers is presented in Table 10. The proposed model tries to minimize the cost of commodities relocation and maximize the demand coverage ratio of damaged area for fair distribution. This is done by a multiobjective possibility linear programming model.
As shown in Tables 12 and 13, considering the inventories of distribution centers and damaged area demands in periods 1, 2, and 3, Birjand and Qaen are considered as active distribution centers. The amount of inventory in active distribution centers is presented in Table 14. Tables 12 and 13 show how to allocate goods from distribution centers to the affected areas. As demonstrated, it has been tried as much as possible to use one distribution center as a closest center to damaged area. However, if the inventory of a distribution center is not sufficient, several distribution centers will be used to fulfill all demands. Reviewing the tables shows that almost all damaged areas receive commodities equally to have fair distribution of commodities. For instance, almost 52% of the 1st commodity demands and 77% of the 2nd commodity demands are met in period 1. Therefore, accuracy and reliability of the model are both acceptable.
6. Conclusions
It is clear that logistic activities are very important in response phase of the disaster management. Also, precise programming can improve the efficacy of the system. According to the crucial role of location-allocation in reducing cost and time, a multiobjective possibility programming model has been developed for relief logistics. This model follows two main objectives. The first one is minimizing good distribution cost and the second one is justice in good distribution process, while the efforts are to maximize the minimum number of satisfied requests. Considering estimated information from problem conditions, some of the model parameters such as the number of damaged area requests and distribution center supply in different periods and transportation cost have been seeded as triangular fuzzy numbers. As observed earlier, the quality of the genetic algorithm is relatively higher, especially in large-scale problems.
Finally, for future research, the relief chain can be reconfigured with three echelons: central disaster facilities of the city as main suppliers, local emergency facilities as distributors, and urban areas as customer. Also, location and routing can be considered simultaneously.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.