Multi-Objective Algorithm for Blood Supply via Unmanned Aerial Vehicles to the Wounded in an Emergency Situation

Unmanned aerial vehicle (UAV) has been widely used in many industries. In the medical environment, especially in some emergency situations, UAVs play an important role such as the supply of medicines and blood with speed and efficiency. In this paper, we study the problem of multi-objective blood supply by UAVs in such emergency situations. This is a complex problem that includes maintenance of the supply blood’s temperature model during transportation, the UAVs’ scheduling and routes’ planning in case of multiple sites requesting blood, and limited carrying capacity. Most importantly, we need to study the blood’s temperature change due to the external environment, the heating agent (or refrigerant) and time factor during transportation, and propose an optimal method for calculating the mixing proportion of blood and appendage in different circumstances and delivery conditions. Then, by introducing the idea of transportation appendage into the traditional Capacitated Vehicle Routing Problem (CVRP), this new problem is proposed according to the factors of distance and weight. Algorithmically, we use the combination of decomposition-based multi-objective evolutionary algorithm and local search method to perform a series of experiments on the CVRP public dataset. By comparing our technique with the traditional ones, our algorithm can obtain better optimization results and time performance.


Introduction
In natural disaster zones where there are serious casualties, it is sometimes necessary to send blood supplies to the injured on their spots instead of bringing them to the hospital for blood infusion due to inaccessibility. With the advancement of unmanned aerial vehicles (UAV) technology (like UAV surveillance [1], traffic monitoring [2], rescue mission [3] and aerial photography [4]) for use in rugged environment, it is easier to substitute UAVs for traditional transportation (such as vehicles, helicopter) in some emergency situations. However, in the event of numerous casualties, a computer program is required to effectively deploy UAVs to their designated locations. This is similar to the travelling salesman problem (TSP) in which the shortest path to be traversed by the salesman to a series of location points is to be determined. However, in real applications, it is much more complicated than the traditional route planning problem, which involves many sub problems such as multiple UAVs' scheduling, blood's temperature change during transportation and the different demand of blood in different areas.
The logistics routing problem was firstly proposed by Assad et al. in 1983[5]. The traditional Vehicle Routing Problem (VRP) focuses on how to make the lowest cost when using only one car to send goods from the warehouse to several customers. This problem has been widely studied in many fields, such as network and logistics. In reality, different conditions and situations also bring different constraints to the traditional problem, thus it is extended continuously. For instance, Capacitated VRP (CVRP) [6] limits the vehicle's carrying capacity, Timedependent VRP (TDVRP) [7] focuses on the transport time constraints, Multi-depot VRP (MDVRP) [8] introduces multiple starting points into the original problem, and Periodic VRP (PVRP) [9] requires the vehicles' route planning to be some of time periodicity. At the same time, people are concerned with not only the economic cost in the optimization process, but also the time cost, no-load ratio, etc. For example, in the green logistics [10,11,12,13], people also study the carbon emission problem where the routes in VRP are related to the vehicle emissions.
VRP is a typical NP-hard combinatorial optimization problem [14]. Specifically, when the scale of the problem is not big, the exact algorithm can be used to obtain the optimal solution, such as integer linear programming [15], dynamic programming [16], etc. But with the increment of scale, it is sometimes difficult to solve it in polynomial time. Therefore, we need to use the approximate algorithm to obtain a near-optimal solution. Currently, many researchers use the typical meta-heuristics algorithm [17] to solve the VRP, such as the tabu search [18], simulated annealing [19] in local search method, the bee colony algorithm [20], ant colony algorithm [21] and genetic algorithm [22] in population search method.
Our proposed blood supply problem is an extension to the original CVRP. In the CVRP, there are n customers in different locations on the map and each customer has a personal need d i . Also, there are several vehicles which have the same carrying capacity Q (Q > max{d i | 1 i n}). We need to arrange and delegate these vehicles to satisfy the need of the injured personnel (our customers). However, each customer can only accept service from one vehicle, each vehicle from the warehouse will finally return without overload during the trip, and the goal of this problem is to minimize vehicles' mileage. But our problem is to supply the blood to several blood needing places by UAVs which has a limited carrying capacity. However, the blood is transported in the form of blood bags and its temperature will affect its quality, so we need to use the heating agent or refrigerant to keep the temperature in an appropriate range, at the same time, the amount of blood needed or the length of path will affect the weight of the required heating agent or refrigerant. Therefore, the payload of the UAV will be reduced. Because a UAV may need to supply several places in one task, and so the longer the trip, the more heating agent or refrigerant it needs. If we only try to minimize the vehicles' mileage, it may increase the number of UAVs that is required. In real application, the distance cost is not the only cost in one trip of the UAV, so it is more meaningful to reduce the number of flights. Therefore, we may treat the distance cost and the number of flights as our objectives, which should be optimized in our model. We use the decomposition-based multi-objective evolutionary algorithm to solve this problem and compare it with the traditional multi-objective CVRP algorithm. The experimental result shows that our method can generate a better solution in an acceptance time.
The rest of this paper is organized as follows. In Section 2, we provide a problem definition of the new UAV-based Capacitated Vehicle Routing Problem (UCVRP). Decomposition-based multi-objective evolutionary algorithm with detailed explanation is discussed in Section 3 to solve the new UVCRP. Section 4 gives experimental results. In our experiments, we obtain the relationship among the weight of blood, the heat agent and the route in a certain circumstance, and we also conduct the comparison on the well-acknowledged benchmark instances. The conclusions and further work are drawn in the last section.

Problem Definition
A group of UAVs U = {u 1 ,. . .,u m } transport the blood bags to each blood needing place from the warehouse. The warehouse (starting point) and several blood needing places forms a bidirectional graph G = (V,E); V = {v 0 ,v 1 ,. . .,v n } that consists of the n + 1 nodes in the graph and E = {(i,j)|i,j 2 V,i 6 ¼ j} is the set of edges; v 0 stands for the starting point while v 1 ,. . ., v n are n blood needing places, and e ij is the flight distance between v i and v j . Each UAV has a limited carrying capacity w, the amount of blood each blood needing place needs is d i , each blood needing place can be supplied at most once, and the UAV will return to the starting point after finishing its task (Fig 1).
During live blood transportation, in order to keep the temperature of the blood within 2-10°C [23], when the environment temperature is below -5°C, some heating agent (hot water) will be placed with the blood bags during the transit. Also, when the environment temperature is higher than 15°C, the refrigerant (ice) will be needed. The hot water and the ice will be called as appendages in the rest of the paper. We hope to have a more reasonable arrangement of the appendages in order to reduce their impact on the payload of UAVs.
In the mathematical model of heat conduction [24], the exothermic process of the hot water is similar to the endothermic process of the ice. For a place v i , the former is defined as follows: The physical meaning of Eq (1) is that the hot water's releasing heat to the environment and the blood results in the reduction of its enthalpy. Eq (2) denotes the enthalpy reduction of the blood owning to the heat releasing from the blood to the environment and the hot water. The temperature of the blood will rise when the heat from the hot water is greater than the heat flow to the environment; otherwise, it will drop. When the initial temperature of the hot water, t 1 (0), the initial temperature of the blood t 2 (0) and the temperature of the environment, t 3 , are known, we can solve the differential equations. Table 1 summarizes the notations used in the heat conduction model.
In this model, c 1 , c 2 , g 11 , g 21 and g 31 are constants and relate to the physical characteristics of blood, water and the material of the container respectively. The flight time, τ, is determined by the flight distance, s i . Therefore, when t 3 , t 1 (0) and t 2 (0) are given as input values of a specific application environment, the weight of hot water, p i , can be calculated according to Eqs (1) and (2).
Usually, the typical UAV that is used in the emergency situation has a limited flight distance. Therefore, we assume that all the blood needing places are within the reach of UAV and the weight of blood (possible appendages) is not exceeding the carrying capacity. In the UAV scheduling process, m UAVs will start from the starting point at the same time and complete the task within one trip. Our goal is to minimize the total mileage and the number of UAVs.
Based on the above model, we can define a complete UAV blood supply route planning problem as follows: x ijk is a Boolean variable, x ijk = 1 indicates that the UAV u k flies from v i the v j , otherwise, y ik is a Boolean variable, y ik = 1 indicates that the UAV u k supplies the v j , otherwise, y ik = 0. S k is a set of nodes visited by the UAV u k . s:t: Objective 1: Eq (3) is used to calculate the total mileage of UAVs, and we minimize the total mileage by obtain the shortest path for each UAV. (4) is used to obtain the optimal scheduling for minimizing the number of UAVs.

Objective 2: Eq
Eq (5) is the constraint for Objective 1 and 2. Constraint 1: The UAV cannot be overloaded. Constraint 2: The number of UAVs must be not more than m. Constraint 3: Each blood needing place can only be visited and served by one UAV. Constraint 4: Each blood needing place is only accessed once. Constraint 5: Only one UAV departs from each blood needing place. Constraint 6: The number of edges passed by each UAV is equal to the number of nodes visited by itself. Constraint 7: Each UAV has to visit the starting point. Constraint 8: All of the routes of UAVs cover all nodes. Constraint 3,4,5,6,7,8: All routes should be formed as loops

Method
UCVRP is a multi-objective problem [25], which attempts to minimize the total distance and the number of required vehicles. The multi-objective optimization is different from the single objective optimization, which will generate a Pareto optimal set having a trade-off between multi objectives. Multi-objective evolutionary algorithms usually initialize a set of random solutions and iteratively generate new ones by selection, crossover, mutation, local search and etc. (as shown in Fig 2) When dealing with the multi-objective problem, we usually make use of the evolutionary algorithm in two ways. The first is based on the definition of Pareto dominance. By using the non-dominated sort [27], the population will be divided into several levels where the lower level means the better and the non-dominated individuals in the evolution will be reserved by elitist strategy. The second is the decomposition-based multi-objective evolutionary algorithm framework [28]. In this way, the original multi-objective optimization problem will be decomposed into a certain number of single objective sub-problems and each sub-problem is a single objective optimization problem. Thus, the traditional single objective optimization algorithm can be extended and applied in this framework. In this paper, we extend the MOEA/D algorithm framework [28] and use it to solve the UCVRP.

Chromosome representation
Path problem usually uses the natural number coding, that is, using a set of natural numbers H = [h 0 ,h 1 ,. . .,h n ] to represent the chromosome. There are two common coding methods [29]: the vehicle-oriented natural number coding, and the destination-oriented natural number coding. Fig 3 shows a solution to the UCVRP, where 0 is the starting point of UAVs, number 1-9 represent various blood needing places and the three routes mean that three UAVs are needed to access all the points in one round. Vehicle-oriented coding will firstly determine the number of delivery sub-routes, then place the destination points into sub-routes and determine the point access order for each subroute (i.e. determine the required number of UAVs, allocate the point required to access for each UAV and determine the order in each sub-route). For example, Table 2 is the vehicle-oriented coding result of the solution in Fig 3, and Target 1 is assigned to Route 1 as the first target to be accessed. The advantage of this method is the decomposition of the problem which reduces the complexity, but it is costly in memory consumption and not intuitive. Destination-oriented coding directly uses a one-dimensional array to represent the route. For instance, Table 3 is the destination-oriented coding result of the solution in Fig 3. Although this method can clearly obtain the optimized route and is easily to be implemented, it often generates many infeasible solutions and is not conducive to the optimization of the solution.
These two methods are commonly used when the number of vehicles or the number of routes is known. Because the total number of vehicles is one of our optimizing goals in this paper, therefore, we choose the destination-oriented coding and eliminate all 0s from it. By partitioning the route from left to right, we can get gene segment which is not exceeding the carrying capacity as a route (as shown in Table 4). This representation of the solution uses less computer memory and is more intuitive; moreover, the solution generated by this method must be a feasible one.

Crossover
In genetic algorithm, the crossover can generate new individuals by swapping the partial structure of two parents and we can use it to greatly improve the searching ability of the genetic algorithm. For the permutation encoding, there are some common crossover algorithms, such as Partially Matched Crossover (PMX), Position-based Crossover, Order Crossover, Cycle Crossover [30]. In this paper, we use the classical crossover algorithm PMX. Since the accessed points are placed in chromosome with the order from left to right, so if we try to directly swap the partial structure of two parents, there may exist some points which appear twice or disappear in children, and such chromosomes will not meet the constraints of the model. Firstly, PMX will randomly select two cutting points X and Y from the two parents, swap the segment  between X and Y from one parent with that from the other parent and record the mapping relation for this swap. For the remaining part of each parent, it's by using the relation that we map the old value to a new value. (As shown in Fig 4)

Mutation
Mutation is usually used in genetic algorithm to increase the diversity of the population. When dealing with the route encoding, in order to avoid the crossing, we use the gene segment to represent the route and do the insertion, swap and reverse to optimize it by various ways [31,32].
(As shown in Fig 5)

Local search
Although the genetic algorithm is a global optimization algorithm, its local search ability is poor since its local search ability is mainly realized by mutation, which is more suitable for searching in a large scale. Therefore, it has a poor performance in the local search, in another word, the fine-tune ability of the mutation is limited. In this paper, the local search heuristics is introduced to improve the local search ability of this algorithm.
In the Holland's schema theorem [33], the representation and the reproduction of the artificial chromosome has been qualitatively analyzed and discussed. Since the schema can also be interpreted as the same structure, it can represent the structure which has the same features in the chromosome. Based on the encoding method in this paper, the schema is a route or a subroute, and we can get the expectation of the number of individuals which contains the schema S in generation t+1 according to the HOLLAND schema theorem, where n t (S) and the n t+1 (S) represents the number of individuals, which contains the schema S in generation t and t+1 respectively; f(S) is the average fitness of individuals contains the schema, while the f avg is the average fitness of the whole population; δ(S) is the defined length of the schema; l is the length of each chromosome; the probability of the crossover and the mutation is P c and P m respectively; O(S) is the number of definite characters in the schema. The determinants in the number of individuals contains the schema in generation t+1 are described in the definition of the schema theorem, which mainly includes the ratio of the average fitness of individuals contains the schema to the average fitness of the population, probability of the crossover and the probability of the mutation. It can be seen that after the operation of selection, crossover and mutation in genetic algorithm, the appearance possibility of loworder schemata, shorter defining length and higher average fitness will increase exponentially. When we use the random connection strategy that nearest-neighbors have the priority, the generated schema and its sub-schema usually have the segment of the optimal solution with low-order schemata, shorter defining length and higher average fitness and these segments are easy to be passed to the next generation. The basic idea of the strategy is that we divide the code sequence of one chromosome into several segments from which some segments are randomly selected at a proper proportion. For each selected segment, we select the first element as the first one for the new segment after reorganization, select the nearest element to the first one from the remaining elements of the original segment as the second one for the new segment, and then select the nearest one to the second element as the third one, until all elements from the original selected segment have been processed. The fitness of the recombined gene code will be calculated and compared with the original one. If the new recombined gene code is dominated by the original one, the new one will be discarded; otherwise the new one will replace the old one. As shown in Fig 6(A), the new segment has a better fitness and will replace the old one. However in Fig 6(B), the old one is better and the new one is discarded.
Because the algorithm needs to repeatedly use the order of distances between all nodes, in order to avoid repeatedly calculating the distances and the orders, we build the adjacency matrix R. Assuming there is a N×N distance matrix D, D ij represents the distance between nodes i and j. R ij is the order of distance D ij among all the distances D ik (k2[0,N-1], k6 ¼i). For example, the R 02 = 1 represents that node 2 is the nearest to node 0 and the R 01 = 2 represents that node 1 is the second nearest to node 0. The detail process of the random nearest neighbor search is as follows: Step 1: Evenly divide the gene code into βM segments (M is the number of UAVs required in this solution while β is used to control the segment size).
Step 3: According to the adjacency matrix R, the selected gene segments are reorganized based on the random nearest neighbor connection strategy.
Step 4: Complete the recombination of β Ã M Ã θ gene segments and generate the new gene code.

Algorithm framework
In this paper, we extend the framework of decomposition-based MOEA and apply the local search with random nearest neighbor to solve the UCVRP. The framework of our proposed algorithm, MOEA/D-N-UVRP, is shown in Fig 7. Besides, we also use two other algorithms, MOEA/D [28] and NASGII [27], in order to solve the UCVRP for comparison. Here, we name them as MOEA/D-UCVRP and NAS-GII-UCVRP respectively. All three algorithms adopt the strategy of mutation and crossover; however, the random nearest neighbor search is our contribution.

Calculation of the proportion between blood and appendage
In our model, the UAV carries blood and appendage, and we can calculate the weight of required appendage by using the weight of blood required and the transport time. Assuming that the temperature of the external environment is -25°C, 2000ml blood that has an initial temperature, T 2 (0) = 5 based on 2.4 kg of hot water with an initial temperature, T 1 (0) = 50, and that the current temperature of the blood transportation container is 20°C. According to the Eqs (1) and (2), T 1 and T 2 will change as shown in Fig 8. It can be seen from Fig 8 that the temperature of the water is decreasing with time whereas the temperature of the blood is increasing during the time as the water temperature is decreasing. In the problem, we must ensure that the temperature of the blood is in an appropriate range during the transit. If the weight of the hot water is relatively heavier than the weight of the blood, the blood temperature will rise rapidly in a short period of time. In Fig 9(A), when the ratio of the weight of hot water to the weight of blood is 4.4, the blood temperature can be kept between 2-10°C (in the appropriate range) within 4.5 hours, and above 10°C between 4.5 hours and 12.6 hours. Then in Fig 9(B), we increase the initial temperature of hot water, the blood temperature will be above 10°C within just 2 hours and last for more than 20 hours. Therefore, during the transit, the change of the blood temperature depends on the initial temperature and the weight of hot water, and it's not practical to increase the initial temperature of the hot water for the goal of reducing its total weight.
In our real application, in addition to keeping the blood temperature within the appropriate range, considering that the blood may not be used immediately, the blood should have a proper temperature in a longer period of time. As shown in Fig 10, in order to use less hot water while keeping the temperature of the blood in an appropriate range after delivery, we can choose the delivery time when the temperature in the descending zone.
In the optimization of the UAVs' scheduling, the evolutionary algorithm requires a lot of iterations, and it will cost a computational time to calculate the weight of required hot water per iteration. In order to reduce the calculation time, we draw a hot water weight proportion table according to the blood weight and transit time in the real situation at a given temperature as shown in Table 5 (The horizontal axis shows weight and the vertical axis shows distance).
For example, when the flight distance to one place is 7 and the weight of blood is 11, the weight of hot water needed is 11×0.06.
Since the experiment is based on the VRP public data set, the data set has only the numerical values without unit and its distribution tends to be discrete. So when calculating the weight of required hot water, we map the weight of blood needed and the distance to the range [0,20] according to Eq 7, where the distance is linear with the time, which can be mapped directly into the time interval. We also note that Eq 7 is a linear transformation. In that equation, CoordX is the mapped value of the blood's weight, CoordY is the mapped value of the distance. Now, MaxX and MinX are maximum and minimum loadage respectively, while MaxY and MinY are maximum and minimum flight distance respectively. However, since the weight of blood and hot water may exceed the carrying capacity of the UAV during transit, there may exist some nodes which cannot be accessed in the data set.
In the experimental process, the starting point also can be considered as a destination point, since the blood requirement and the flight distance is zero, so we have MinY = MinX = 0 here. When mapping the blood weight, if the weight of blood is close to the MaxX, the minimum ratio of the weight of hot water to the weight of blood can be set as 0.02:1. In order to ensure that the weight of blood and hot water will not exceed the maximum carrying capacity of the UAV, the MaxX can be set as w / (1 + 0.02). As for the mapped value of distance MinY, we can set it directly for different data instances.
It can be seen from Table 5 that the change of the proportion is not uniform; the proportion value increases along the time axis from small to large and decreases with the blood weight.  While the MinX for blood weight has been given above, the MaxY for distance is not set yet. As shown in Fig 11, when the MaxY takes a large value, the proportion between hot water and blood is easy to be concentrated in the smaller area on the time axis; however, if the MaxY takes a small value, the proportion is easy to be concentrated in the larger area. Fig 11 shows the distribution of the proportion caused by different MaxY values. If the MaxY takes a larger value, it will result in a larger requirement of the hot water, then the total weight UAVs need to carry will increase and the more UAVs will be needed.

Parameter setting
Parameters used in the experiments are shown in Table 6. All the algorithms are written using Java programming language and operate on a computer (Core I3 CPU, 2.93GHZ, 4G memory space). All the operating computers adopt single-threaded execution.

Two-objective: total distance and number of UAVs
In traditional CVRP, the total distance and the number of vehicles cannot be the two partially conflicting objectives. If the total distance is the minimal one, we cannot reduce it by adding additional vehicles. In our problem, we should consider increasing the payload and reducing the carry of the hot water. If the number of UAVs is small, then the UAVs will need to fly a  longer time. However, that means we need more hot water. Thus, the total distance and number of UAVs are both relevant and competitive. In Table 7, the MOEA/D-N-UAV is the algorithm proposed in this paper for UCVRP, the NSGAII-CUVRP is based on the traditional NSGAII for CUVRP and the MOEA/D-CUVRP is based on the MOEA/D. We run these three algorithms on nine instances according to different proportions between hot water and blood. Each algorithm has been executed five times and we take the average of the results. It's shown that the MOEA/D-N-UAV is better than the MOEA/ D-CVRP and the NSGAII-CVRP in terms of optimization, and has an average increase of 3% compared with MOEA/D-CVRP. The detail optimal solutions of MOEA/D-N-UAV on instance E-n23-k3 and E-n101-k14 are illustrated in Table 8, Figs 12 and 13. In the actual experiments, although the total distance and the number of UAVs show some competition, they are not completely antagonistic. As shown in Table 7, for instance E-n76-k7, the MOEA/D-N-UCVRP can obtain two optimal solutions while the NSGAII-CUVRP only has one. From the results of MOEA/D-N-UCVRP, the total distance by 20 UAVs is less than the one by 19 UAVs. Also, we can adjust the MinY to generate more solutions which indicates that introducing the concept of the appendage into UCVRP makes distance and weight be associated and competitive so that generate the UCVRP.

Performance analysis
The main metrics of the multi-objective algorithm are Hypervolume (HV) and Inverted Generational distance (IGD) [34]. HV computes the size of the region that is dominated by a set of non-dominated solutions, based on a reference vector that is constructed using the worst objective values of each objective. It is defined as where the Pareto front (PF) denotes the set of non-dominated sets, O is the objective space and HV(f) is the set of objective vectors dominated by f. Distance from Representatives in the PF (IGD-metric): Let P Ã be a set of uniformly distributed points along the PF. Let A be an  11, 13, 9, 7, 0} approximation to the PF, the average distance from A to P Ã is defined as: where d(v, A) is the minimum Euclidean distance between v and the points in A.
However, in our model, the total distance and the total number of UAVs are not completely antagonistic. The optimization results of these three algorithms only have little non-dominated solutions each time, usually, the number is 1 to 2. Then the PF may have only one solution, which cannot be calculated by Eqs (8) and (9).
We also conduct an experiment about the iteration times and the time cost of these three algorithms. Fig 14 shows  iterations, the time cost of three algorithms shows a linear increase, the fluctuation of the NSGAII-CUVRP's time cost is relatively large while the other two are more stable. Also, the time cost of MOEA/D-N-UCVRP is very close to MOEA/D-UCVRP.

Conclusion
In this paper, the blood supply in emergency situation is introduced into the traditional Capacitated Vehicle Routing Problem, and a new multi-objective optimization problem is proposed according to the two factors of distance and weight. Aiming to achieve a real-life blood supply logistics, we study the blood's temperature change due to the external environment, the heating agent (refrigerant) and time factor during transit, and propose an optimal method for calculating the mixing proportion between blood and appendage in different circumstance. Algorithmically, we use the combination of decomposition-based multi-objective evolutionary algorithm and local search method, in order to conduct a series of experiments on the VRP public dataset, which is included in S1 File. By comparing our technique with the traditional ones, we proved that our algorithm can obtain better optimization results and time performance. Our model can also be used in the food warming and route planning of food delivery, the transportation of frozen seafood and etc via UAVs.