Optimization of “vehicle-UAV” joint distribution routing for cold chain logistics considering risk of epidemic spreading and green cost

To address the epidemic, such as COVID-19, the government may implement the home quarantine policy for the infected residents. The logistics company is required to control the risk of epidemic spreading while delivering goods to residents. In this case, the logistics company often uses vehicles and unmanned aerial vehicles (UAVs) for delivery. This paper studies the distribution issue of cold chain logistics by integrating UAV logistics with epidemic risk management innovatively. At first, a "vehicle-UAV" joint distribution mode including vehicles, small UAVs and large UAVs, is proposed. The green cost for vehicles and UAVs is calculated, respectively. The formula for infection risk due to large numbers of residents gathering at distribution centers to pick up goods is then derived. Furthermore, based on the control of infection risk, an optimization model is developed to minimize the total logistics cost. A modified ant colony algorithm is designed to solve the model. The numerical results show that the maximum acceptable risk and the crowd management level of distribution centers both have significant effects on the distribution network, logistics cost and number of new infections. Our study provides a new management method and technical idea for ensuring the needs of residents during the epidemic.


Introduction
In the context of accelerating quality upgrading in the consumer market, people's quality requirements for fresh goods have gradually increased, and market demand has continued to expand, which has led to the rapid development of cold chain logistics.The normal operation of logistics system is an indispensable part of social sustainability.However, the outbreak of COVID-19 epidemic made the cold chain logistics and distribution of fresh goods face extremely severe challenges.Under the emergence of the epidemic, "national online grab vegetables" led to a surge in demand for fresh goods, and online orders needed to be guaranteed by offline cold chain logistics distribution.Traditionally, the fresh goods are carried by vehicles to the distribution center and further delivered to the residents waiting here.During this process, face-to-face contact exists widely, increasing the risk of cross-infection between deliveryman and community residents.As a result, there may be more infected residents and even more deaths.Therefore, the traditional distribution mechanism cannot well solve the problem of fresh distribution under the background of sudden epidemic [1].The sustainable development of society is under the threat of COVID-19 pandemic.
As a new tool, unmanned aerial vehicle (UAV) has achieved rapid development in the field of logistics and transportation during the past few years of COVID-19 pandemic.Particularly, during the COVID-19 pandemic, the Chinese government implemented the home quarantine policy for epidemic areas with traffic control and travel restriction.Many infected residents were put in quarantine timely and had to wait for the home delivery service offered by the logistics company.Furthermore, the logistics company was ordered to control the risk of epidemic spreading while delivering goods to residents.As a result, many companies introduced the "vehicle-UAV" joint operation service mode, which performed the "non-contact" distribution service effectively.On one hand, this mode can alleviate the shortage of distribution personnel caused by the epidemic and reduce the risk of crowd infection in the process of residents picking up their own goods.On the other hand, it can avoid the high logistics distribution cost and huge value loss in the transportation process of fresh goods.In fact, vehicles and UAVs were used together in the cold-chain logistics of many Chinese cities during the COVID-19 pandemic.The optimization of "vehicle-UAV" joint distribution scheme is valuable for reducing the cost of the logistics company and controlling the risk of epidemic spreading.
Therefore, this paper focuses on the cold chain logistics distribution of fresh goods using "vehicle-UAV" joint operation in the context of the epidemic.The main contributions are: (1) In view of the infection risk caused by the gathering of some residents at the pick-up point, this paper uses the SI modeling method of infectious disease to derive the analytical formula of risk, and takes the control of infection risk as the basis for optimal selection of distribution plan.(2) Our study integrates UAV logistics with epidemic risk management, providing a new management method and technical idea for ensuring the needs of residents during the epidemic.It promotes the application of UAV technology in the field of fresh logistics and epidemic risk infection management.(3) Based on the differences in size, load and endurance, this paper divides UAVs into small UAVs and large UAVs, and studies the "vehicle-small UAVs-large UAVs" fresh cooperative delivery operation during the epidemic period.The comparison between this paper and existing literature can be found in Table 1.The results of this paper provide new ideas for improving the efficiency of fresh distribution, timely meeting the demand for fresh commodities of residents and reducing the risk of epidemic transmission during the epidemic period.As far as we know, this is the first research work to study the optimization of "vehicle-drone" joint distribution scheme from the quantitative perspective of controlling the risk of epidemic spreading.The rest of this article is organized as follows.Section 2 introduces the literature review of LRP and "vehicle-UAV" joint operations in cold chain logistics under the background of epidemic.The third part discusses the construction of the "vehicle-UAV" integrated logistics optimization distribution model for fresh goods under the control of the spread risk of the epidemic.Section 4 introduces the improved ant colony algorithm to solve the model.A numerical example is given in Section 5. Finally, some conclusions are drawn in section 6.

Literature review
Location-routing problem (LRP) is a classical combinatorial optimization problem, which needs to determine both the location of candidate facilities and the path of vehicles [2].Since Cooper [3] first proposed LRP, it has gradually become a hot issue in the fields of supply chain, logistics and operation research optimization.Location-path planning has been applied in many fields, such as waste recycling, emergency logistics, low-carbon, closed-loop supply chain, etc.For example, Li et al. [4] developed a novel multi-objective optimization model for the LRP for biomass waste collection using tailoring evolutionary algorithm.Adarang H et al. [5] addressed a location-routing problem (LRP) under uncertainty for providing emergency medical services (EMS) during disasters, which is formulated using a robust optimization (RO) approach.Jiang et al. [6] established a low-carbon open location-routing problem model related to distribution center size and distribution path, and designed a quantum evolutionary algorithm to solve it.Govindan et al. [7] studied the location-inventory-routing problem to design a circular closed-loop supply chain network with carbon tax policy for achieving circular economy.Guo et al. [8,9] not only considered the location-inventory decisions for closedloop supply chain management in the presence of the secondary market, but also solved a multi-commodity location-inventory problem in a closed-loop supply chain with commercial product returns.At the same time, location-path planning has also been introduced into the fresh cold chain industry.In recent years, scholars have conducted a lot of research in the field of LRP of fresh commodity, mainly considering the influence of product perishability, time window limit, low-carbon environmental protection and other factors on the model.Rahmanifar et al. [10] designed a novel nonlinear multi-objective model designed to concurrently optimize warehouse facility location and vehicle routing, addressing the challenges inherent in cold chain logistics processes.Moghaddasi et al. [11] developed a mixed integer nonlinear programming routing-location model to improve the cold chain logistics network design within the Balanced Score Card pillars by maximizing the customer satisfaction and minimizing unit cost and greenhouse gas emissions in order to solve the optimization problem of the marine product distribution logistics system.Ma et al. [12] established a multi-agent optimization model of fresh agricultural distribution location-route based on conflict cooperation, considering the high distribution cost of fresh agricultural products and high product loss, and took into account customer fuzzy time window in the model.Li et al. [13] built a multi-vehicle cold chain logistics vehicle routing optimization model with carbon emission cost and time window on the basis of considering the congestion index.Govindan et al. [14] proposed a multi-objective optimization model for the network distribution of perishable food supply chain considering the time window constraints, and studied the cold chain location path from the perspective of low carbon and environmental protection, so as to minimize the impact of carbon emissions on the environment.From the perspective of low-carbon and environmental protection, Wang et al. [15] considered the carbon footprint of cold chain logistics location path optimization according to the characteristics of perishable products.Leng L et al. [16] analyzed the low-carbon LRP of cold chain logistics, established a dual-objective model with the minimization of total logistics cost and vehicle waiting time, and selected a variety of algorithms for solution comparison.Pan [17] established a cold chain logistics route optimization model from the perspective of low-carbon, and solved the model by using Wolf pack algorithm and ant colony algorithm.Mohammadi et al. [18] developed a multi-objective mixed-integer non-linear programming model for a four-level sustainable supply chain of a perishable product with price-dependent demand and deterioration rates.In addition, some scholars have carried out extended studies on the influence of factors such as facility capacity limitation [19], influence of multiple vehicle types [20], and inventory cost [21,22] on LRP.The above literature indicates that considering the influence of location problem in fresh distribution route planning can effectively reduce the overall logistics distribution cost.
In view of the particularity of the epidemic background, vehicles cannot directly complete the distribution service for isolated residents, and the collaborative distribution characteristics of vehicle-UAV should be considered in the path planning process.As for the research on vehicle-UAV cooperative distribution, Murray and Chu [23] earlier proposed the traveling salesman problem of truck and UAV cooperative distribution.Aiming at the optimal route and scheduling of UAV and delivery truck, they proposed a mixed integer linear programming model for two UAV delivery problems.This new package delivery paradigm facilitates lastmile delivery.At present, the research on vehicle-UAV collaborative delivery is mainly differentiated from the launch and landing position selection of UAVs, the number of UAVs that vehicles can carry, and the collaborative operation mode of UAVs.For example, Carlsson and Song [24] consider that the UAV takes a package from the truck, and at the same time, the truck continues its route.After delivering the package, the UAV returns to the truck to take another package, and the vehicle and the UAV can meet at any position.Chang and Lee [25] proposed a new method of nonlinear programming model, which used K-means clustering to group customers and determine the UAV launch location, so as to expand the UAV distribution area and reduce the driving route of vehicles.Boysen et al. [26] considered the situation that a truck could carry multiple UAVs, distinguished whether multiple UAVs or one UAV were placed on the truck, and whether the take-off and landing stations must be the same.Karak et al. [27] proposed a mathematical formula and an effective solution for the hybrid vehicle-drone routing problem for pick-up and delivery services, considering that vehicles only serve as the carrier of UAVs and all delivery services are completed only by UAVs.Ham [28] did not regard vehicles as UAV carriers, and UAV and trucks were distributed in parallel, and studied the comprehensive scheduling of multiple warehouses carrying truck fleets and UAV fleets to achieve excellent operation.Moshref-Javadi et al. [29] considered the load capacity and flight time limit of UAV.Although some research results have been achieved on the problem of vehicle-UAV collaborative delivery, there is a lack of research and discussion on the effective cooperation of "vehicle-UAV" and the provision of distribution services for normal and isolated residents under the epidemic environment by setting vehicle temperature control.Gao et al. [30] studied the problem of coordinated distribution of emergency resources between unmanned ground vehicles (UGVs) and unmanned aerial vehicles (UAVs) based on nested vehicle paths, including two sub-problems: how to accept the commander's operation order and how to generate nested vehicle routing planning.
In addition, a small number of literatures have carried out research on the location and routing problem of fresh food under the background of epidemic.For example, Li et al. [31] studied the retail logistics distribution vehicle routing problem considering the order release time under the emergency epidemic environment.The results show that the model and algorithm can provide effective decision support for the efficiency improvement and cost control of retail logistics distribution.Zhang and Li [32], in order to effectively solve the "contact-free" distribution problem of fresh commodities under the current epidemic background, proposed the site-routing optimization problem of fresh distribution based on the epidemic background, taking into account the impact of the epidemic, commodity perishability, temperature control characteristics, site selection and vehicle-UAV collaborative distribution route planning.Shibu and Philip [33] considered the context of COVID-19 and proposed a distribution system using UAVs, in which delivery agencies and customers work together to distribute goods, and established a new framework to find the optimal path of delivery UAVs.Based on the practical experience of COVID-19 vaccine delivery in China, Zheng et al. [34] proposed a hybrid approach of machine learning and evolutionary computing, using an evolutionary algorithm to determine vehicle routes and distribute vaccines daily from satellite/warehouse to vaccination points.
However, there are the several shortcomings in the existing literatures: (1) The impact of crowd infection risk on the "vehicle-UAV" delivery path, cost and number of infected people during the epidemic is not considered from a quantitative perspective.(2) There is lack of research about the effect of home quarantine policy on the distribution scheme.(3) The studies of "vehicle-UAV" logistics distribution during a major epidemic in the literature has several defects in the size setting of UAV, such as not considering the coordination of large and small UAVs, and the restrictions on load and endurance.(4) There is lack of research on the location and path of fresh distribution based on "vehicle-UAV" collaborative distribution under the impact of the epidemic.
In view of this, this paper focuses on the cold chain logistics distribution of fresh commodities under the home quarantine policy of an emergency epidemic and studies the risk of mass infection and logistics distribution network under the joint operation of "vehicle-UAV".The main differences between our study and the existing literature are summarized in Table 1.

Problem description
The epidemic, which spreads through the air (COVID-19), and the home quarantine policy are considered throughout this paper.If some resident is infected and shows the symptoms, he will be sent home and put in quarantine timely.
A logistics company plans to deliver the goods to residents, who live in a great number of housing estates.The company can use two distribution methods-vehicle and UAV.The location-routing logistics network of fresh green distribution in the epidemic environment considered in this paper is shown in Fig 1 .It can be seen that, there are several distribution centers and the unique logistics center.Every housing estate falls within the scope of its unique distribution center and the vehicle can provide services to the residents only through the distribution center.The vehicle carries the goods from the logistics center to a series of distribution centers and returns to the logistics center finally.When the vehicle arrives at a distribution center, the goods are unloaded and sent to the residents waiting here.If a small number of infected residents are put in quarantine at home, several small UAVs (Fig 2A) carried by the vehicle can send the goods to those residents directly.For the residents waiting at the distribution center, some of them may be infected but asymptomatic and they may spread the virus to others unintentionally.As a result, the risk of epidemic spreading rises.Moreover, if there are many infected residents in a housing estate, the company can send a team of large UAVs (Fig 2B ) to deliver goods.The large UAV flies from the logistics center to some housing estates and returns to the center finally.Because it is non-contact, the UAV delivery does not result in more infections.
This paper aims to study how the logistics company makes the optimal decision about the integrated "Vehicle-UAV" distribution network with controlling the risk of epidemic spreading.The assumptions about the studies on the spread of epidemics in this paper, which mainly stems from the classical SI epidemic model, are given as follows.
The main assumptions about the spread of epidemics A1.The population is divided into infected people and uninfected people, and the total population remains unchanged.
A2.It takes no more than one day for a logistics company to distribute goods in the city, and it takes significantly more than one day for an infected person to be cured.Therefore, this paper does not consider the situation that the infected person is cured and becomes healthy again.A3.A healthy person becomes infected with the virus by coming into close contact with an infected person while waiting in line at a distribution center to pick up goods.Eventually, the epidemic spread.
A4.The rate at which a population is infected with the virus is proportional to the base rate of infection, the average number of people each infected person comes into contact with per unit time, the probability of a single infected person coming into contact with a healthy individual, and the number of infected people.
The above assumptions will help us model the spread of epidemics in distribution logistics, calculate the number of new infections and evaluate the risk of epidemic spreading, as shown by Section 3.3.2.

Model parameters and variables
According to the needs to build the model, this paper sets the following parameters and variables, as shown in Table 2.

Objective function analysis of model. 1. Fixed costs.
The fixed costs mainly conclude the operation, maintenance and labor costs [15,35].The fixed cost of vehicles is where c V k includes the related cost of the small UAVs going with a vehicle.Similarly, the fixed cost of large UAVs is Then the total fixed cost is 2. Green cost.The green cost equals to the fuel cost plus the carbon tax, which further depend on the amount of fuel consumption.According to the literature [36][37][38][39][40], the fuel consumption of vehicles is where is the fuel consumption per unit distance of a vehicle.The fuel cost and the carbon tax for all vehicles are respectively Table 2.The explanation of parameters and variables.

Parameters and Variables
Meaning 0 The unique logistics center station

N
The set of distribution centers  Thus, the green cost for vehicles is Similarly, the green cost for large UAVs is where The total green cost of logistics system is 3. Refrigeration cost.Vehicle and large UAV are the transportation facilities of logistics system.However, most of the large UAVs being used within the city do not equip with a refrigerator because of the limited size and rapid delivery speed.Therefore, the refrigeration cost is only related to the vehicles.It is 4. Delay penalty cost.During the epidemic, the residents who are infected with the virus and have symptoms of infection are required to stay at home, while other residents can walk to the given distribution center and get the goods.Due to the seriousness of epidemic, logistics companies must deliver the goods to the specified center within the required time; otherwise, they have to pay for the delay penalty cost.When vehicle k arrives at the center n at time t k n , an asymptomatic resident has to walk from their housing estate s to center n and wait about � t n .After all these, this resident gets his goods.Thus, the delay penalty cost is

T s
The time that team w of large UAVs needs to stay at housing estate s � t � n The average waiting time of residents at center n to get the goods

R max
The maximum acceptable risk, which is the acceptable upper limit for the number of new infections caused by residents gathering at a distribution center to pick up goods https://doi.org/10.1371/journal.pone.0306127.t002 be computed as where v 1 is the average pick-up speed of residents, X k2K X s2SðnÞ ð1 À a s Þq s y ns f k n is the quantity of goods which belong to the asymptomatic residents waiting at center n.Moreover, because the large UAV has a much faster speed than vehicles in reality, we do not consider the delay case of large UAVs.

Risk analysis of the residents waiting for goods at the distribution center.
The risk of infection stems from the close contact among residents who get together at the distribution center and wait in lines to get their goods.It is because there may be infected but asymptomatic residents in the crowd.Section 3.3.2aims to quantify the risk of infection, helping us design reasonable constraint conditions later.Throughout this paper, the risk of a distribution center is denoted as its number of new infections.According to the classical SI model [41], the change of total infections in center n over time is where δ is the infection rate when a healthy resident has close contact with one infector, β n is the average number of persons that a resident closely contacts when he is waiting at center n to get the goods, I n and N n denote the number of infection and the number of residents getting together at center n, respectively.It is clear that, δ will decrease if more residents are persuaded by managers at center n to wearing masks, and β n will decline if more residents are commanded to keep distance with each other while waiting in line at center n.Thus, overall infection rate δβ n measures the management level of center n.The management level increases as δ and β n decrease.In addition, N n À I n N n in Eq (15) is the probability of a single infected person coming into contact with a healthy individual.Eq (15) means that the number of new infections, i.e., the risk of distribution center n, is That is (17) can be calculated as where coefficients a s ; âs have been explained in Table 2.
In terms of Eqs ( 14), ( 17)-( 19), we have the risk formula of each distribution center 3.3.3Model design.According to the results in Section 3.3.1-3.3.2, the optimization model for the emergency routing problem considering the risk of virus infection is proposed as follows.

Subject to
Risk n � R max ; 8n 2 N ð22Þ y ns ¼ e s ; 8s 2 SðnÞ; n 2 N ð23Þ X s2SðnÞ q s y ns � E n ; 8n 2 N ð24Þ X s2SðnÞ q s y ns ¼ r n ; 8n 2 N ð25Þ In this model, the objective function (21) aims to minimize the sum of fixed cost, green cost, refrigeration cost and delay penalty cost, as shown in Eqs (1,(11)(12)(13).Eq (22) represents a constraint that the risk of each distribution center, as shown in Eq (20), cannot exceed upper limit R max .Note that R max is the acceptable upper limit for the number of new infections caused by residents gathering at a distribution center to pick up goods.Undoubtedly, there will be no new infection if all residents are serviced by large UAVs.However, the large UAV has a much higher operation cost than the vehicle and the quantity of large UAVs is usually very limited.Therefore, the logistics company has to maintain a balance between the use of vehicles and large UAVs.The cost is that, when trucks arrive at one distribution center and unload the goods, the residents have to gather at the center to pick up their goods, resulting in more infections due to direct close contact.Eq (23) means that distribution center n will service its own housing estate s (y ns = 1) if s is not serviced by large UAVs (e s = 1).Eq (24) implies that the goods shipping to a distribution center cannot exceed its capacity.Eq (25) the demand of every housing estate must be satisfied.Eq (26) sets that each distribution center can be serviced by only one vehicle.Eq (27) guarantees that a vehicle must travel to a center if the vehicle aims to service it.Eq (28) indicates that each vehicle can pass through a distribution center only once.Eq (29) means that any vehicle cannot stop at each distribution center finally.Eqs (30)(31) are computational formulas of Q k 0 ; Q k m , the goods weight of a vehicle just arriving at a distribution center.Eq (32) represents that the goods weight of a vehicle starting from the logistics center cannot exceeds its capacity.Eqs (33,44) are the constrains of subtour-elimination for vehicles and large UAVs, respectively.Eq (34) is the formula of the waiting time the vehicle needs at a center, where v 1 is the average pick-up speed of residents, v 2 is the flying speed of small UAVs and φ 2 is the total capacity of all small UAVs for a vehicle.Eq (35) indicates that a large UAV team is assigned to service a housing estate only if the estate is not serviced by any vehicle.Eq (36) implies that a large UVA team can fly to a housing estate only if the team has been assigned to service the estate.Eq (39) sets that any large UAV team cannot stop at any housing estate finally.Eqs (40)(41) are computational formulas of D w 0 ; D w s , the goods weight of a team of large UAVs just arriving at a housing estate.Eq (43) represents that the goods weight of a large UAV team starting from the logistics center cannot exceed its capacity.Eq (45) is the formula of the waiting time the team of large UAVs needs at a housing estate, where v 3 is the flying speed and φ 3 is the total capacity of a large UAV team.
Particularly, there is an obvious overlap of service object between the vehicle and the large UAV.For example, once one housing estate has been serviced by a team of large UAVs, the related distribution center nearby will not provide the goods to this housing estate.And the goods of distribution center are delivered by vehicles.

Solution algorithm
In this section, a modified ant colony algorithm (MACO) is used to solve the new model, which belongs to the NP-Hard problem.It is well known that NP-Hard problems are often solved by intelligence algorithms, such as genetic algorithm (GA) [15], ant colony algorithm (ACO) [42,43] and simulated annealing algorithm (SA) [44].Among these algorithms, ACO is well suited and effective.This is because the foraging routes of ants are highly similar to the transport routes of vehicles and drones in VPR/LPR, helping us flexibly improve ACO according to the characteristics of distribution routes of cars or drones.
However, many versions of ant colony algorithms (i.e.[15,42,43]) cannot be applied to our model directly because there is an obvious overlap of service object between the vehicle and the large UAV in our model.For the same reason, our model cannot be solved directly by some famous commercial solution software, such as CPLEX and GUROBI.To address this problem, we set two kinds of ants-one is the vehicle ant (VEA) and another one is the team of large UAV ant (UAA).The two kinds of ants look for food in rotation, every VEA only visits distribution centers directly and services related housing estates indirectly, and every UAA only visits housing estates and provides direct service.
Our modified ant colony algorithm (MACO) integrates the foraging routes of VEA and UAA, ensuring the flexible cooperation between the two types of ants.Then the overlapping issue of service objectives between cars and drones is solved.In the same iteration, specifically, for one VEA and one UAA, VEA first forages for food from the logistics center to some distribution center.For a candidate distribution center, there are several related housing estates.VEA must choose the housing estates which he should service if VEA decides to visit the distribution center.Meanwhile, VEA must avoid servicing any housing estate that has been serviced by some UAA.Once VEA has returned to the logistics center, UAA starts from the logistics center to some housing estate.Within the same iteration, UAA must avoid any housing estate which has ever been indirectly serviced by a vehicle.The above "VEA-UAA" cycle repeats until all housing estates have been serviced.After that, the VEA-pheromone and the UAA-pheromone are updated respectively, just like the classical ant colony algorithm [43].Then this iteration ends and a new iteration starts.
According to the above discussion, we introduce the key details, subroutines and algorithm flow of MACO as follows [45]. (

1) Calculation of transition probability
In each iteration, the transition probability of an VEA Ant VEA k visiting distribution center j from another distribution center i is where where the three terms on the right-hand side of Eq (47) represent the green cost, the refrigeration cost and the delay penalty cost, which are derived in Eqs ( 8), ( 12) and ( 13), respectively.Similarly, the transition probability of an UAA Ant UAA s visiting housing estate j from another housing estate i is where UN þ i is the candidate set of the next feasible housing estate, ξ ij is the amount of pheromones, which is derived in Eq (9). (

2) Updating strategy of pheromone
At first, the pheromone τ ij of VEA Ant UAA s in each iteration is updated by where ρ2(0,1), Q>0 and TD VEA k is the length of the complete path (from the logistics center to the center itself) that the k-th VEA Ant VEA k has built.Second, the pheromone ξ ij of UAA TD UAA k is updated by where ρ2(0,1), Q>0 and TD UAA k is the length of the complete path (from the logistics center to the center itself) that the k-th UAA Ant UAA s has built.
(3) Choice of the feasible service object in the region of j for VEA Suppose one VEA Ant VEA k has reached distribution center i and starts to evaluate another distribution center j.Denote the set of housing estates belonging to the scope of service of distribution center j as S(j) = {s 1 ,s 2 ,� � �,s G }. Ant VEA k selects the housing estates in the region of j to service according to the following procedure (Note that the service object for VEA is housing estates, not distribution centers).
Step 1: Set k = 1 and the set of feasible service objects F(j) = {}.
Step 2: For s k 2S(j), if S k has not been serviced by VEA or UAA in the same iteration, go to the next step; otherwise, go to Step 5.
Step 3: If Q k i À r i À q s k > 0, go to the next step; otherwise, go to Step 5.
Step 4: Suppose that Ant VEA k visits from i to j and services the housing estates in F(j)[{s k }.If Risk j �R max according to Eq (20), set Step 5: Set k = k+1.If k>G, stop and output F(j) (the set of feasible service objects for Ant VEA k ); otherwise, go to Step 2. In the above procedure, Q k i À r i À q s k > 0 in Step 3 means that Ant VEA k has enough goods to meet the demand of housing estate s k after visiting distribution center i, Risk j �R max in Step 4 implies that the risk can be controlled if s k is added to the set of service object of Ant VEA k , i.e. s k can be set as one feasible service object of Ant VEA k .
(4) Choice of the feasible service object for UAA Suppose that one UAA Ant UAA k has reached housing estate i.Denote H = {j 1 ,j 2 ,� � �,j h } as the set of housing estates which have not been serviced by VEA or UAA in the same iteration.Ant UAA k constructs the set of feasible service objects UN þ i based on the following procedure.
Step 1: Set k = 1 and the set of feasible service objects Step 2. In the above procedure, D w i À q i À q j k > 0 in Step 2 means that Ant UAA k has enough goods to meet the demand of housing estate j k after visiting housing estate i, then housing estate j k is included in UN þ i . (

5) Stopping criterion of MACO
In this paper, MACO will be terminated if one of the following criteria is satisfied:

PLOS ONE
1. Stop if no better solution can be generated in consecutive L iterations (L is a pre-defined counter value).
2. Stop if maximum iterations MAXGEN is reached.

Algorithm flow of MACO
Step 0. Set K VEAs and K UAAs, the maximum number of iterations MAXGEN, the maximum number of consecutive iterations that can be tolerated if no better solution is produced L, iteration number gen = 0.
Step 1.If stopping criterion procedure ( 5) is satisfied, stop; otherwise, set gen = gen+1, go to the next step.
Step 4. Update the pheromone of VEA and UAA according to procedure (2).Set gen = gen +1, go to Step 1.

Numerical studies
In this section, we first test the effectiveness of the modified ant colony algorithm (MACO) proposed in Section 4. Then the distribution data of a cold chain Logistics Company is used to verify the model and obtain useful results.

Algorithm experiment
In this part, the genetic algorithm (GA), the simulated annealing algorithm (SA) and MACO are used to solve the model of Section 3.3.3.
The test data is called as n-m-k, where n represents the number of housing estates, m represents the number of distribution centers, k represents the location distribution rule of housing estates.The locations of housing estates are uniformly distributed if k = 1; randomly distributed if k = 2; half-and-half if k = 3.The unique logistics center is randomly distributed.Because each housing estate can be serviced by at most one distribution center, it is necessary to determine the locations of distribution centers.To do so, we use the square grid to divide the whole plane into several regions.Then the centre-of-gravity method is used to calculate the location of the unique distribution center corresponding to the distribution centers within each region.It is set that the capacity of each distribution center is 150 tons, the total demand of each housing estate is 1500kg, every family needs 5 kg goods and the ratio of infected residents who are sent home and put in quarantine is 0.05.In addition, maximum acceptable risk R max = 30 and overall infection rate of each distribution center δβ n = 0.1.The parameter settings of these algorithms do have obvious effects on the solution quality.Thus, we first refer to the relevant literature [15,[42][43][44][45] which focused on the intelligent optimization algorithms for solving LPR.According to the above literature, we set the maximum iterations is 300 for GA and MACO.The other parameters are set as: (1) For GA, population size Np = 100, crossover probability p c = 0.8, and mutation probability p m = 0.2; (2) For SA, initial temperature T 0 = 100, temperature cooling coefficient α = 0.95, and termination temperature T f = 0.05; (3) For MACO, number of ants N = 100, pheromone importance factor α = 3, heuristic importance factor β = 2, pheromone evaporation factor ρ = 0.2, and total pheromone deposit Q = 500.Based on the above settings, the solution quality can be guaranteed in the following experiments.
The results are shown in Table 3 and Fig 3, where t is the time spent at the end of the algorithm.

Model experiment
In this section, the distribution data of a cold chain Logistics Company is used to verify the model and obtain useful results.The company mainly provides warehousing and distribution services for refrigerated foods, such as dairy products, chilled meat, and so on.It has the only one logistics center and 11 candidate distribution centers in a certain area with a capacity of 150 tons (the locations are shown in Fig 4).Before arranging a delivery task, a total of 50 orders of demand points (housing estates) is received.The specific information such as the geographical location, latest service demand time, and demands are shown in S1 Appendix.Every family needs 5 kg goods and sends one person to take delivery of goods, if the family is serviced by a distribution center.
The effects of risk constraint R max and overall infection rates δβ n on the minimum total cost are presented in Figs 5-8.
The minimum total cost changes little when R max increases from 20 to 30, decreases rapidly with R max if R max falls into interval [40 120], then remains stable if R max is larger than 120.In fact, a higher value of risk constraint means that more infections, as well as more residents  getting together at the distribution center, are acceptable.And the large UAV has a higher unit cost than the vehicle.In order to reduce the total cost, as a result, the logistics company tends to replace the large UAVs with vehicles partly when R max increases gradually.But this result is invalid if R max is too small or large.Meanwhile, it can be seen from Fig 6 that, as R max increases, the number of total new infections keep unchanged first, then increases quickly and increases slowly at last.Therefore, Figs 5 and 6 show that there is an obvious negative correlation between the minimum total cost of logistic company and the number of total new infections if R max changes.That is to say, the decrease of R max reduces the minimum total cost at the expense of the heath of residents.
As shown in Fig 7, the minimum total cost first increase quickly and then keeps unchanged as the overall infection rate of each distribution center δβ n increases continuously.The cost reaches the highest level when δβ n slightly increases to 0.9.Meanwhile, Fig 8 shows that the number of new infections also increases with δβ n obviously.Note that a smaller value of δβ n implies a higher management level.Thus, the improvement of management can reduce the total cost of logistic company and the number of new infections significantly at the same time.For the logistic company and the residents, the improvement of management of each distribution center is a win-win strategy, while another strategy-change of risk constraint presents a situation of "zero-sum game" (Fig 6).More specifically, δβ n depends on δ (the infection rate when a healthy resident has close contact with one infector) and β n (the average number of persons that a resident closely contacts when he is waiting to get the good).So, what the mangers of each distribution center should do is to persuade or even force every resident waiting at the center to wear masks (reduce δ) and keep distance with each other (reduce β n ).For example, if everyone wears a mask, the infection rate when a healthy resident has close contact with one infector can reduce to 0.15%.
In the following, we further explore the effect of risk constraint R max on the distribution path scheme of vehicles and large UAVs.The distribution paths are compared between the high acceptable risk (R max = 140) case and the moderate acceptable risk (R max = 70), as shown by Figs 9 and 10 and Tables 4-6.
From the above results, we can see that R max has a significant effect on the joint distribution scheme of large UAVs and vehicles.As R max decreases from 140 to 70, the number of teams of large UAVs increases from 1 to 5 and the number of housing estates serviced by large UAVs increases from 3 to 19.Moreover, the vehicle distribution paths change obviously and the number of the housing estates serviced by vehicles decreases from 47 to 31.There is a positive correlation between the number of housing estates serviced by large UAVs and R max , and a negative correlation between the number of housing estates serviced by vehicles and R max .In addition, it is found that the number of vehicles and their distribution paths barely change as R max increases, the only difference is the load of vehicles and the housing estates serviced.

Managerial implications
In summary, our model in this paper mainly aims to optimize the "Vehicle-UAV" joint distribution strategy of fresh goods during the epidemic, and the objective function is minimizing the total cost of a logistics enterprise based on considering the effect of home quarantine policy and controlling the risk of epidemic spreading.Thus, the sustainable development of society in the context of a major epidemic is truly realized.However, the total cost of considering the risk of epidemic spreading is obviously higher than the total cost when the risk is ignored for the design of distribution path optimization of cold chain logistics.This implies that it is necessary to pay a certain amount of economic costs to control the risk of epidemic spreading and realize the sustainable development of society towards today's increasingly widespread threat of various infectious diseases.
From the governmental agencies' point of view, first of all, what they should do is to pay attention to the risk issue of epidemic spreading during the logistics distribution process, strictly supervise the logistics enterprises and raise their awareness of this risk.Secondly, the local government must carefully assess the spatial distribution of the epidemic and available medical resources.On this basis, the government should set standards for epidemic prevention, i.e., the acceptable upper limit for the number of new infections per day, and negotiates with logistics companies on the proportion of drones and trucks used in the delivery process.Third, because the drone-delivery service does not result in more infections, the government can carry out subsidy policies and encourage enterprises to increase the proportion of drone deliveries to a certain extent.
From the cold chain logistics companies' point of view, first of all, they should raise awareness of controlling the risk of epidemic spreading, strictly abide by the epidemic prevention policies implemented by the government and introduce the risk of cluster infection into the "Vehicle-UAV" joint distribution strategy optimization of cold chain logistics.As a result, the total cost and the number of new infections can be reduced simultaneously, and the demand of residents can also be satisfied timely.Second, the logistics enterprises must try hard to improve the management level of each distribution center because many residents gather there to get the goods.They need to effectively guide residents to wear masks and keep their distance to avoid close contact while waiting in line to pick up goods.Third, when epidemic breaks out, it may be not necessary for the logistics company to significantly change his daily vehicle distribution paths, he only needs to reduce the load of vehicles and sent more large UAVs to deliver the goods.Particularly, the paths of large UAVs must be updated.
When our model is applied in practice, there may be some operational, technical and regulatory challenges.For example, the cost of procurement, maintenance and use of drones, especially large drones, may be high; The joint delivery scheme requires the UAV to have accurate positioning function and real-time information sharing technology, and also requires the staff to be familiar with the standard operating technology of the UAV; In the process of joint distribution, the joint coordination and supervision of vehicles, small drones and large drones may have certain difficulties.

Conclusions and future research
This paper addresses the distribution routing problem for cold chain logistics during the pandemic, such as COVID-19 which spreads through the air.To control the risk of epidemic spreading, the home quarantine policy and an integrated "vehicle -UAV" transportation mode are introduced.In this case, the formula for the risk of epidemic spreading, which relates to the situation that residents get together at the distribution center to get the foods, is derived.Then the federated distribution network of vehicles and large UAVs is designed and a new logistics optimization model is developed.
Using a modified ant colony algorithm to solve the model, we study the effect of the maximum acceptable risk R max and crowd management level of each distribution center in detail.The main conclusions are: (1) There is a negative correlation between the total cost of logistics company and the new infections, i.e., the decrease of R max reduces the minimum total cost at the expense of the heath of residents.(2) the improvement of crowd management of distribution center reduces the total cost and the number of new infections significantly at the same time.(3) The maximum acceptable risk R max has a significant effect on the distribution network and number of large UAVs.But this effect is not obvious for the paths of vehicles.This implies that the logistics company may not need to change his usual delivery paths of vehicles clearly when epidemic outbreaks.Instead, the company should reduce the load of vehicles, sent more large UAVs to deliver the goods and significantly updated the distribution network of large UAVs.
In the future, the development of artificial intelligence, Internet of Things, and 6G communication technology will further promote the effective operation of the "vehicle-drone" joint distribution model.For example, under the influence of the above technologies, multiple drones can realize real-time information sharing among each other, respond to the needs of residents in a more timely manner, and the drone group can automatically arrange distribution tasks and optimize the distribution network according to the needs of residents.As a result, distribution efficiency will be significantly improved.
Although our model and numerical results provide insights to the emergency logistics during the epidemic, there are still some limitation and problems which can be further considered.First, this paper only considers the example of a single logistics center, and the case of multiple logistics centers remains to be studied.Second, some uncertain events, such as the sudden intensification of the epidemic in a community, resulting in the failure of the previously specified distribution method to meet the demand, need to be further researched.A two-stage stochastic programming model with recourse should be developed.Finally, compared with the classical ant colony algorithms, our MACO algorithm has higher theoretical time/space complexity in one iteration because two kinds of foraging ants are set.In Section 4, we only use small and medium-sized numerical examples to test the effectiveness of MACO.Therefore, a numerical experiment based on real-world large-scale urban transportation networks is necessary for MACO in the future.

0mnn0
team w of large UAVs passes from housing estate s to housing estate j, p w sj ¼ 1; otherwise p w sj ¼ 0 r n The total goods received by distribution center n Q k The goods weight of vehicle k starting from logistics center 0 Q k The goods weight of vehicle k just arriving at center m t k The time of vehicle k arriving at center n T k The time that vehicle k needs to stay at center n D w The goods weight of team w of large UAVs starting from logistics center station 0 D w s The goods weight of team w of large UAVs just arriving at housing estate s t w s The time of team k of large UAVs arriving at housing estate s (Continued )

2 . 1 .
The k-th VEA Ant VEA k leaves the logistics center 0 and visits a series of distribution centers gradually and services a series of housing estates.The operations of collecting the feasible service objects and feasible visiting center for Ant VEA k are conducted by procedure (3).The operations of determining the next visiting center and related service objects for Ant VEA k are conducted by procedure (1).When Ant VEA k returns to the logistics center 0 and builds a complete path Path VEA k , go to Step 2.2.Step 2.2.The k-th UAA Ant UAA k leaves the logistics center 0 and visits a series of housing estates gradually.The operations of collecting the feasible service objects, i.e. the feasible visiting center, are implemented by procedure (4).The operations of determining the next visiting housing estate are implemented by procedure (1).When Ant UAA k returns to the logistics center 0 and builds a complete path Path UAA k , go to Step 2.3.Step 2.3.For the "Vehicle-UAV" joint distribution scheme Path VEA k þ Path UAA k , calculate the total cost Total_cost k .

Table 1 . Main differences between our study and the existing literature.
https://doi.org/10.1371/journal.pone.0306127.t001 0 is the set of all distribution centers processing freights sjThe time of a team of large UAVs flying from housing estate s to housing estate jβ nThe average number of persons that a resident closely contacts when he is waiting at center n to get his own goods s Latest arrival time of housing estate s to get the goods θ Unit delay penalty cost e s If housing estate s is serviced by some vehicle, e s = 1; if s is serviced by large UAVs, e s = 0 y ns If housing estate s is serviced by a vehicle passing through center n, y ns = 1; if s is serviced by large UAVs, y ns = 0

Table 2 .
(Continued) is the delay time for vehicle k to service distribution center n, d sn is the distance between n and housing estate, v walk is the average walking speed of residents and LT s is the latest demand service time for housing estate s.Moreover, the time � t n can

Table 3 . Experimental results of MACO and other algorithms.
As shown in Table3 and Fig3, for the solutions of the above 9 examples, MACO always has a certain advantage over GA and SA in terms of reducing total cost of the model.According to the computation time, MACO is superior to GA and SA in more than 77% of cases.Therefore, the MACO is competitive in solving our model in Section 3.3.3.