An Improved Intelligent Auction Mechanism for Emergency Material Delivery

Abstract

Emergency material delivery is vital to disaster emergency rescue. Herein, the framework of the emergency material delivery system (EMDS) with the unmanned aerial vehicle (UAV) as the vehicle is proposed, and the problem is modeled into a multi-trip time-dependent dynamic vehicle routing problem with split-delivery (MTTDDVRP-SD) in combination with the rescue reality, which provides decision support for planning disaster relief material. Due to the universality of dynamic interference in the process of material delivery, an optimization algorithm based on the traditional intelligent auction mechanism is proposed to avoid system performance degradation or even collapse. The algorithm adds pre-authorization and sequential auction mechanisms to the traditional auction mechanism, where the pre-authorization mechanism improves the capability performance of the system when there is no interference during the rescue process and the sequential auction mechanism improves the resilience performance of the system when it faces interferences. Finally, considering three types of interference comprehensively, which includes new task generations, task unexpected changes and UAV’s number decreases, the proposed algorithm is compared with DTAP (DTA based on sequential single item auctions) and CBBA-PR (consensus-based bundle algorithms-partial replanning) algorithms under different dynamic interference intensity scenarios for simulation experimental from two perspectives of the capability performance and resilience performance. The results of Friedman’s test with 99% confidence interval indicate that the proposed algorithm can effectively improve the capability performance and resilience performance of EMDS.

1. Introduction

Earthquakes, floods and other large-scale natural disasters occur frequently, causing serious damage and far-reaching impacts on modern society. How to deal with natural disasters and their consequences can benefit the maintenance of social stability, and this is an emerging topic that has garnered increasing attention from both practitioners and academics in recent years. According to relevant research reports [1], the Wenchuan earthquake in 2008 caused 69,277 deaths, 374,643 injuries, 17,923 people missing, and direct economic losses of more than 800 billion US dollars and the 2015 earthquake in Nepal caused at least 8786 deaths, 22,303 injuries, and economic losses of about 113 billion yuan. The way to reduce the impact of natural disasters on human life, the economy and the environment requires not only to be prepared for detection and early warning before disasters occur [2], but also rapid post-disaster relief supplies and post-disaster recovery work, which led to the emergence of an EMDS. Due to the destruction of ground communication facilities, road traffic, and increased environmental hazards, the selection of vehicles in EMDS presents challenges. While the use of UAVs as vehicles in EMDS can weaken the impact of ground traffic on rescue work more than manned aircraft and unmanned vehicles and avoid the loss of rescue equipment and personnel, UAVs as temporary communication relay technology are also very mature, which can replace communication base stations to restore certain communication capabilities in disaster areas [3].

The impact of the complex post-disaster environment is the urgent demand for materials by the disaster victims. Moreover, with the delay of time, the deterioration of the health of the disaster victims will further aggravate the urgency of the relief materials in the disaster area, resulting in greater casualties [4]. Song et al. [5] proposed an emergency material scheduling model with multiple disaster points and multiple rescue points, but the model is mainly aimed at the emergency material scheduling of disaster points with the same emergency indicator of material demand. Based on the weight-TOPSIS method, Guan et al. [6] determined the emergency indicator of the tasks in the disaster area, considered the penalty cost of unmet needs in the disaster area, and they established a collaborative disaster relief material scheduling model, but did not consider the change of the emergency indicator of the tasks. At the same time, due to the aftershocks after the earthquake and the crash of the UAVs, some other situations have brought great dynamics to the post-disaster rescue. On the other hand, considering the loading capacity constraints of UAVs and battery life constraints, UAVs are often required to be reused, and the material needs of a single disaster area may not be satisfied by a single UAV at one time. These are all issues that EMDS should consider when completing disaster relief, so we modeled it as a multi-trip time-dependent dynamic vehicle routing problem with split-delivery (MTTDDVRP-SD), based on the classic multi-trip vehicle routing problem (MTVRP) [7], which also takes into account the following characteristics: multi-trip per UAV, time-dependent emergency indicator, dynamics of the environment, split-delivery, and trip duration limit.

MTTDDVRP-SD takes into account various practical situations in the delivery of emergency materials, which makes it unique and difficult. From the perspective of models, MTTDDVRP-SD can be seen as a combination of MTTDVRP, DVRP and SDVRP, although these three models are more in line with the real problem, but the relevant research is not sufficient, especially the combination of the three is less.  

• In the delivery of emergency materials, it is common to allow vehicles to make multiple trips, and it is also widely used in logistics [8], public transportation scheduling [9] and other fields. Cattaruzza et al. [10] emphasized that allowing vehicles to make multiple trips can significantly improve vehicle utilization and reduce the number of vehicles required, thereby further reducing fixed costs and overall costs, and the application in life is more realistic. Moreover, they theoretically proved that the MTVRP problem was NP-hard. Dorling et al. [11] argued that existing VRPs are insufficient for planning UAV deliveries, that is because multiple trips to the depot are not permitted, resulting in too many UAV solutions, so they proposed two multi-trip VRPs for UAV deliveries to address both problems. Sun et al. [12] proposed a time-dependent multi-trip vehicle routing problem (TDMTVRP) with an improved travel speed model, arguing that speed is time-dependent. Moreover, the experimental results indicated that compared with the capacitated vehicle routing problem with the time windows model (CVRPTW), the TDMTVRP model proposed can both decrease the vehicle utilized times dramatically and shorten the vehicle travel distances slightly in dealing with the vehicle routing problem (VRP).
• Dynamic VRP (DVRP) [13] emerged as a new variant of VRP and attracted more and more attention. Li et al. [14] investigated the emergency material dispatching problem under uncertainty and built a multi-objective optimization model to minimize the total cost and maximize the affected rate of demand, but without considering the case of multiple trips of vehicles. Nan et al. [15] investigated the dynamic demand VRP with time windows for the model of mixed distribution of electric and fuel logistics vehicles, i.e., the demand of the task changes randomly or a new task is generated randomly during the execution of the task. Sabar et al. [16] considered the DVRP problem under vehicle speed uncertainty caused by traffic congestion, but this uncertainty has regularity. Ouaddi et al. [17] proposed the MTDVRP problem with a combination of MTVRP and DVRP and recognized the complexity of such problems. Their work considered the problem of vehicle overloadable transportation and distribution, with dynamics reflected in the cancellation and modification of user orders and changes in vehicle speed due to congestion, but without considering the problem of delivery splitting at a task.
• The split-delivery vehicle routing problem (SDVRP) is a variant of the Capacitated Vehicle Routing Problem (CVRP) in which a customer can be served by more than one vehicle, i.e., the requirement of a customer can be split among vehicles. Casazza et al. [18] addressed a single commodity Pickup and Delivery Vehicle Routing Problem and assumed that the requirements of both pickup and delivery nodes can be split, and each node can be visited more than once. An extensive experimental analysis showed their method to offer simultaneously more modeling flexibility and more computational effectiveness. Hartomo et al. [19] investigated the problem of distribution of relief materials from a logistics center to a network of demand points, and their assumptions included that the requirement for materials at a node would be served by several trucks or over multiple trips. However, their work does not take into account the dynamic changes in the environment and emergency indicator over time. Zhang et al. [20] considered the case of a customer requiring multiple vehicles for service, but without considering the change in capability due to the change in external information during the working process, and the case of multiple trips.  

Currently, to solve these kinds of problems, the exact algorithm [21] and memetic algorithm [17] are commonly used, but they are only suitable for solving small-scale arithmetic cases. On the other hand, most of the processing methods of DVRP are dividing time slices and decomposing into a series of static VRP problems [22,23], and then using heuristic algorithms to solve them, but under the conditions of multi-trip and delivery splittable, the solution search space is huge and the time complexity is very high, so much so that it is not suitable for emergency material delivery scenarios under environmental uncertainty. For this purpose, this paper considers an intelligent auction algorithm based on the market mechanism, such as Guo et al. [24] who proposed an auction algorithm that introduces the idea of “blockchain” to realize the task allocation problem in the “last mile” logistics distribution. Auction algorithms have been improved in the literature, Choi et al. [25] proposed a consensus-based auction algorithm (CBBA), which uses both auction and consensus mechanisms for task selection and conflict resolution, thus reducing the computational resource consumption of the algorithm and accelerating the convergence. The performance improvement of CBBA algorithm over the initial auction algorithm was experimentally verified, but the research does not address the dynamic change of task information. In a later study, Buckman et al. [26] proposed the CBBA-PR algorithm for partial replanning based on the CBBA algorithm in order to cope with dynamic changes in the environment by mixing a subset of unfinished tasks with new tasks to achieve fast reallocation. CBBA-PR ensures optimality in dynamic assignment and improves the fast response of the algorithm to dynamic situations, but the algorithm is not effective in some cases where real-time performance is required or the problem size is large. Chen et al. [27] argued for the better applicability of online algorithms for the uncertain arrival of dynamic tasks and proposed that tasks should be assigned one by one rather than in batches to ensure the efficiency of UAV execution. Farinelli et al. [28] converted multi-robot patrol into a task assignment problem, considered online decision making techniques, proposed the DTAP algorithm based on sequential single auction, and verified the performance of the algorithm on a real robot platform. However, they relied only on the task information at the current moment, without considering more global known information. Of course, centralized algorithms such as the genetic algorithm (GA) and particle swarm optimization (PSO) can also solve VRP, but they are less adaptable to the dynamic situation of changing environment in the disaster area. For example, Luo et al. [29] pointed out that when the global information of the task changes dynamically over time, the centralized algorithm tends to decompose and re-iterate the constraint solution, which is less flexible than the auction algorithm.

When EMDS is implemented to start the emergency material delivery task, it is still necessary to pay attention to two aspects, i.e., not only to ensure the capability performance to avoid secondary damage to the disaster area due to untimely rescue, but also to ensure the resilience performance to have the ability to continue the task. In fact, accurate and timely information about the disaster area can be obtained through efficient reconnaissance detection [30]. However, the uncertainty of the environment in a disaster is common and may cause the system to collapse, thus requiring the EMDS to adapt effectively to interferences, i.e., to have a good resilience performance. While resilience was first introduced to engineering systems in 1995 [31], defined as, “Resilience is the probability of a system’s recovery from failures”, in 2009, Woods [32] argued that the focus should be on the adaptability of a system facing more complex situations, thus giving the definition, “Resilience, as a form of adaptive capacity, is a system’s potential for adaptive action in the future when information varies, conditions change, or new kinds of events occur, any of which challenge the viability of previous adaptations, models, or assumptions”. As to how to quantitatively evaluate resilience, summarizing the research results over the years, it can be broadly divided into two categories: time-dependent and time-independent models.

The core of time-dependent models is that the performance of systems (e.g., power transmission systems [33] and communication systems [34]) is time-dependent and can be expressed quantitatively, and such models usually describe the resilience of the system in terms of performance change curves. In contrast, time-independent models are mostly used to describe the system capability performance only calculated after the completion of the task, such as air defense systems [35]. The integral resilience model is commonly used for quantitative description of resilience, which describes the system resilience by comparing the cumulative distributional impact of system performance over time with and without interference, and can be found in the literature [36]. Dessavre et al. [37] studied the resilience of offensive and defensive countermeasures based on the integral resilience metric, using the number of OODA cycles as a performance indicator, but without considering the multi-interference situation, so some improvements are needed in emergency material delivery.

Motivated both by the adaption of new technologies and delivering vehicle (UAV) in practical transportation industry and the theoretical gap existing in the current literature, we investigate a problem model (MTTDDVRP-SD) that is more adapted to the emergency material delivery scenario in disaster areas. On this basis, due to the existence of two layers of optimization of UAV routing and delivery quantity, which leads to a huge solution search space of the traditional evolutionary algorithm and does not adapt to the objective demand of the dynamics of the disaster area. So, we try to design a new algorithm based on the intelligent auction mechanism and propose a resilience evaluation indicator under dynamic multiple interferences. Specifically, the main contributions of this paper can be summarized as follows.  

• To meet the performance requirements of emergency material delivery, an overall framework of EMDS with UAV as the vehicle is proposed, and the emergency material delivery problem is described as a new variant of VRP (MTTDDVRP-SD);
• Based on the traditional intelligent auction algorithm, pre-authorization and sequential auction mechanisms are added, and sequential auction algorithm based on pre-authorization (SAABoPA) is proposed to solve MTTDDVRP-SD to optimize the system performance.
• Based on the integral resilience model, a multi-interference resilience evaluation indicator suitable for emergency material delivery scenarios is proposed;
• The experiments are carried out under different dynamic interference intensities through the post-earthquake emergency material delivery scenario and compared with the comparative algorithms from two perspectives of capability performance and resilience performance to verify the effectiveness of the algorithm in this paper.  

The problem MTTDDVRP-SD will be deeply analyzed and defined in the next section, while in Section 3, the proposed method of SAABoPA for solving it is presented. The experimental results are reported in Section 4. Section 5 concludes this study and discusses future research directions.

2. Problem Definition and Formulation

2.1. Problem Definition

UAVs with loading capacity were considered to deliver relief materials to the earthquake-affected areas, and information about the danger in the affected areas was obtained by satellite remote sensing technology, UAV reconnaissance technology, and ground rescue teams [38]. The relief materials are concentrated in the distribution center at the border of the disaster area, and the UAV makes task assignment according to the disaster area danger information, determines the task target, loads the materials at the distribution center, and then carries out the delivery to the disaster area target. The specific EMDS framework is shown in Figure 1, which describes the interaction between the system and the outside. The EMDS includes information processing subsystem (IPSys), command and control subsystem (CCSys), and loading and transportation subsystem (LTSys).

The subsystem functions can be described as follows:  

• IPSys: Responsible for filtering, standardizing and sharing task information. Task information provided by the disaster environment (DE) is filtered, standardized, and shared with CCSys and LTSys. The sharing of task information is discontinuous and will only be re-shared when there is a sudden change in the information.
• CCSys: Responsible for task assignment and solution sharing. Based on the task information shared by IPSys, the task assignment is executed by combining the status and capability of each LTSys component. Due to the dynamic change of task information, task assignment is not executed at once. After the allocation plan is formed, it is shared to the material distribution center (MDC).
• LTSys: Responsible for task assignment and material delivery. Based on the task information shared by IPSys, the task assignment is completed with the unified scheduling of CCSys. After the assignment plan is determined, the UAVs in LTSys load the materials stored in the MDC and execute the material delivery process.  

The relevant assumptions of this paper are as follows:

• The material requirements of the task may not be met by one UAV at a time;
• A single deployment of a UAV may execute multiple tasks;
• When the UAV is empty, it needs to be returned to the MDC and reloaded to be able to execute the next task;
• Obstacle avoidance among UAVs and against obstacles (e.g., mountains and forests) is not considered. Obstacle avoidance can be accomplished by adjusting the height of the UAV under certain conditions [39], and the increase in range due to height change can be ignored because the height difference is small relative to the range;
• In the case of emergency rescue, the UAV always maintains the maximum capacity, i.e., the UAV always maintains the maximum power during the task. The UAV endurance can be considered as fixed for a certain battery capacity;
• Ignoring the deceleration of the UAV during the delivery of the material, as well as the change in direction, the time variation from this is insignificant in comparison. Theoretically, this variation could be integrated into the velocity and loading relationship, but we do not consider it in order to simplify the problem. Moreover, the kinematic model of the UAV is not investigated in depth in this paper;
• Consider using a single base, i.e., an MDC. For a general multi-base VRP, it can be divided into sub-blocks by clustering methods, and within the sub-blocks still remains a single base VRP;
• The material in the MDC meets the material needs of all affected areas and only one type of material is considered. The current national reserve of emergency resources is sufficient, and the amount of emergency resources available can be considered adequate in the face of a catastrophic event in a local area. In fact, many types of materials are needed in disaster relief situations, including food, water, medicine, etc. However, in order to facilitate the emergency and effective implementation of relief, in many cases we do not count the requirements of various types of materials too much and make a precise delivery, and the more general practice is to synthesize various materials into a rescue kit [40].

2.2. Problem Formulation

In this section, the mathematical model proposed to solve the above problem is introduced. However, before that we need a formal description of the disaster area environment and dynamic interferences.

2.2.1. DE Formulation

Here, we first formalize spatio-temporal characteristics and dynamic properties of DE based on task information of the disaster area.

Definition 1.

Task environment G. The task environment is modeled as an undirected graph $G=(V,E)$, where V denotes the set of vertices in Euclidean space, E denotes the set of edges.

In G, V includes the task node T and the MDC node $V_B$. The edge $E\subset V\times V$ indicates that two neighboring vertices are reachable between them and are paths for the UAV to move. A task node represents a delivery task, and the set of all tasks is $T=\{T_1,T_2,\dots ,T_j,\dots ,T_{\left|T\right|}\}$.

Definition 2.

Time t. Time is represented by a set of discrete sequences, $t\in \{0,1,2,\dots \}$.

Take the time of the first UAV taking off for the task as the base time 0.

Definition 3.

Coordinates $(x,y)$. The node position is represented by the two-dimensional coordinates $(x,y)$.

The node coordinates considered in this paper include the task coordinates and MDC coordinates, and the task coordinates are used to characterize the nearby affected areas that need relief materials.

Definition 4.

Quantity required $R\left(t\right)$. At t moments, the quantity of material required by the task $T_j\in T$ is denoted by $R_j\left(t\right)$.

To simplify the problem, the quantity required is unitized, e.g., $R_j\left(t\right)=3$ means that the quantity of material required for the task $T_j$ at time t is 3 units.

Definition 5.

Emergency indicator $E\left(t\right)$. At t moments, the emergency of the task $T_j\in T$ on the material is denoted by $E_j\left(t\right)$.

The emergency indicator [6] indicates the urgency of the requirement for materials at the task, which varies with time and the delivery of materials. When the upper threshold is reached, it will cause secondary damage to the disaster area due to the lack of replenishment of materials for a long time [41]. This is similar to a time window constraint, but the difference is that the emergency indicator can describe the range of time constraints, the value of the task is varying, i.e., the sooner the task is completed the higher the value. The formula is expressed as follows.

$$E_j\left(t\right)=\left\{\begin{array}{cc}0,\hfill & R_j\left(t\right)\le 0,\hfill \\ 1,\hfill & E_j\left(t\right)\ge 1,\hfill \\ E_j\left(t_0\right)+a·(t-t_0)-\frac{E_j\left(t_0\right)}{R_j\left(t_0\right)}·M{D}_{ij}\left(t\right),\hfill & \mathrm{others}.\hfill \end{array}\right.$$

(1)

where $t_0$ denotes the moment of the last IPSys shared information. Equation (1) indicates that the emergency indicator $E_j\left(t\right)$ of a task gradually increases with time t, but the material replenishment decreases the value of $E_j\left(t\right)$. In fact, the task emergency indicator formula is a fit to the actual indicator, and IPSys shares the information every time the task information changes, i.e., the task emergency indicator can be corrected at the $t_0$ moment.

2.2.2. Dynamic Interference Formulation

The dynamic scenarios in this paper consider three possible types of interferences in which their occurrence is unpredictable, and are therefore described in terms of stochastic processes and defined as follows.

Definition 6.

Dynamic interference type 1: new task generations. This means that during the delivery of emergency materials, new tasks will emerge as the disaster continues.

Taking the stochastic process $\left.\left\{\tilde{\mathbf{X}}\left(n^1\right),n^1\in \{1,2,\dots ,\left|DI{S}^1\right|\}\right.\right\}$ denotes the occurrence of dynamic interference events of the type new task generations, where $\tilde{\mathbf{X}}\left(n^1\right)=(\tilde{t}\left(n^1\right),\tilde{x}\left(n^1\right),$$\tilde{y}\left(n^1\right),\tilde{R}\left(n^1\right),\tilde{E}\left(n^1\right))$, and each random variable $\tilde{t}\left(n^1\right),\tilde{x}\left(n^1\right),\tilde{y}\left(n^1\right),\tilde{R}\left(n^1\right),\tilde{E}\left(n^1\right)$ satisfies a certain probability distribution, respectively, $\tilde{t}\left(n^1\right)\sim DiscreteU(0,T_{end})$, $\tilde{x}\left(n^1\right)\sim U(\underline{x},\overline{x})$, $\tilde{y}\left(n^1\right)\sim U(\underline{y},\overline{y})$, $\tilde{R}\left(n^1\right)\sim DiscreteU(\underline{R},\overline{R})$, $\tilde{E}\left(n^1\right)\sim U(\underline{E},\overline{E})$. The result of the occurrence of dynamic interference events of this type causes the task T to be expanded, and the information of the new task is all generated randomly in the random vector $\tilde{\mathbf{X}}\left(n^1\right)$, which is represented by the following equation.

$$T=T\cup \left\{T_{\left|T\right|+n^1}\right\},t=\tilde{t}\left(n^1\right)$$

(2)

Definition 7.

Dynamic interference type 2: task unexpected changes. This means that during the delivery of emergency materials, information on the material requirements and emergency indicator of currently known tasks will deteriorate as the disaster continues and disaster information is reevaluated.

Taking the stochastic process $\left.\left\{\tilde{\mathbf{Y}}\left(n^2\right),n^2\in \{1,2,\dots ,\left|DI{S}^2\right|\}\right.\right\}$ denotes the occurrence of dynamic interference events of the type task unexpected changes, where $\tilde{\mathbf{Y}}\left(n^2\right)=(\tilde{t}\left(n^2\right),\tilde{j}\left(n^2\right),$$\widetilde{\Delta R}\left(n^2\right),\widetilde{\Delta E}\left(n^2\right))$, and each random variable $\tilde{t}\left(n^2\right),\tilde{j}\left(n^2\right),\widetilde{\Delta R}\left(n^2\right),$$\widetilde{\Delta E}\left(n^2\right)$ satisfies a certain probability distribution, respectively, $\tilde{t}\left(n^2\right)\sim DiscreteU(0,T_{end})$, $\tilde{j}\left(n^2\right)\sim DiscreteU(1,|T\left|\right)$, $\widetilde{\Delta R}\left(n^2\right)\sim DiscreteU(0,\overline{R}/2)$, $\widetilde{\Delta E}\left(n^2\right)\sim U(0,\overline{E}/2)$. The result of the occurrence of dynamic interference events of this type causes changes in the information about the material requirements and emergency indicator in the current task T, and the amount of changes in the information are generated randomly in the random vector $\tilde{\mathbf{Y}}\left(n^2\right)$, which is represented by the following equation.

$$R_{\tilde{j}\left(n^2\right)}\left(t\right)=R_{\tilde{j}\left(n^2\right)}\left(t\right)+\tilde{R}\left(n^2\right),t=\tilde{t}\left(n^2\right)$$

(3)

$$E_{\tilde{j}\left(n^2\right)}\left(t\right)=E_{\tilde{j}\left(n^2\right)}\left(t\right)+\tilde{E}\left(n^2\right),t=\tilde{t}\left(n^2\right)$$

(4)

Definition 8.

Dynamic interference type 3: UAV’s number decreases. This means that during the delivery of emergency materials, the number of available UAVs will be reduced due to failure and other circumstances.

Taking the stochastic process $\left.\left\{\tilde{\mathbf{Z}}\left(n^3\right),n^3\in \{1,...,\left|DI{S}^3\right|\}\right.\right\}$ denotes the occurrence of dynamic interference events of the type UAV’s number decreases, where $\tilde{\mathbf{Z}}\left(n^3\right)=(\tilde{t}\left(n^3\right),\tilde{i}\left(n^3\right))$, and each random variable $\tilde{t}\left(n^1\right)\sim DiscreteU(0,T_{end})$, $\tilde{i}\left(n^3\right)\sim DiscreteU(1,|U\left|\right)$. The result of the occurrence of dynamic interference events of this type causes the number of currently available drones is reduced, which is represented by the following equation.

$$U=U-\left\{U_{\tilde{i}\left(n^3\right)}\right\},t=\tilde{t}\left(n^3\right)$$

(5)

2.2.3. MTTDDVRP-SD

In summary, The MTTDDVRP-SD is formulated as follows.

$$f=min{\sum }_1^{N_T}{E}_j\left(T_{end}\right)$$

(6)

$$0\le L_i^k\left(t\right)<L_i^{max}$$

(7)

$$L_i^k\left(t\right)=L_i^k\left(s{T}_i^k\right)-{\sum }_{j=1}^{\left|T\right|}{x}_{ij}^k\times M{D}_{ij}\left(t\right)$$

(8)

$$L_i^k\left(s{T}_i^k\right)>0$$

(9)

$$0\le M{D}_{ij}\left(t\right)\le L_i^k\left(t\right)$$

(10)

$$a{T}_i^k-s{T}_i^k\le T_i^{duration}$$

(11)

$$v_i\left(t\right)=v_i^{unload}-{\delta }_i\times L_i^k\left(t\right)$$

(12)

$$1\le E{x}_i$$

(13)

$${\sum }_{k=1}^{E{x}_i}{x}_{ij}^k\in \{0,1\}$$

(14)

$${\sum }_{j=1}^{\left|T\right|}{x}_{ij}^k\ge 1$$

(15)

$${\sum }_{i=1}^{\left|U\right|}{\sum }_{k=1}^{E{x}_i}{x}_{ij}^k\ge 1$$

(16)

$$0\le M{D}_{ij}\left(\mathrm{t}\right)\le R_j\left(t\right)$$

(17)

$$x_{ij}^k\in \{0,1\}$$

(18)

Equation (6) aims to minimize the number of tasks whose emergency indicator reaches the threshold 1 at the end of the emergency material delivery. Different from the general objective function of the VRP and its variants problem, which considers minimizing transportation and inventory costs or maximizing transportation and inventory rewards [42,43]. In the disaster area emergency material delivery problem, this paper considers, based on the rescue work after the earthquake disaster, to ensure the safety of victims’ lives and avoid more secondary injuries. Equation (7) indicates that the loading of any UAV at any time is less than the maximum loading and is non-negative. Equation (8) represents the change in the current UAV loading. Equation (9) represents the loading capacity of the UAV $U_i$ when it is deployed from the MDC is strictly large 0. Equation (10) indicates that the quantity delivered by the UAV $U_i$ to the task $T_j$ is less than the current loading capacity of the UAV i. Equation (11) represents the flight time constraint for a single UAV trip. Equation (12) represents the relationship between the UAV speed and the loading capacity [44]. Equation (13) indicates the constraint on the number of trips of a single UAV. Equation (14) indicates the entire emergency material delivery process, the UAV $U_i$ to the task $T_j$ at most once. Equation (15) indicates that a single UAV executes at least one task in a single trip. Equation (16) indicates that the entire emergency material delivery process, task j is executed at least once. Equation (17) indicates that the quantity of UAV $U_i$ to task $T_j$ delivery is less than the current requirement of task $T_j$. Equation (18) represents the decision variables.

3. Method: SAABoPA

3.1. Overall Framework

In the EMDS framework, CCSys and LTSys use the SAABoPA algorithm for task assignment based on the task information shared by IPSys. The tasks shared by IPSys are equivalent to auction items. CCSys acts as the auctioneer and sends the task information to the UAVs in the corresponding LTSys. The UAV acts as a bidder, calculates the net benefit for each task with the task price stated by CCSys, and sends the bid price of the task back to CCSys. CCSys decides the pre-authorization and authorization of the task based on the bid price, and sends the task to the authorized UAV and the overall framework of SAABoPA is shown in Figure 2.

The SAABoPA algorithm requires the UAV to obtain at most one task authorization and one task pre-authorization in one auction. The algorithm is based on pre-authorization to ensure the system capability performance, and based on sequential auctions to effectively improve the system resilience performance and avoid the frequent occurrence of interference that causes the system performance degradation by reallocation or adjustment of the assignment results. It should be noted that since the published task is executed in the future time and has no effect on the current execution process, the UAV can participate in the auction during the execution of the task [27], but it can only obtain the pre-authorization of the task instead of the authorization.

There are two possibilities for starting an auction:  

• CCSys active start: in case of interference.
  • In EMDS, when the task information changes abruptly (due to new task $T_{new}$ generation, unexpected changes in task information and completion of one delivery), IPSys shares the current task information with CCSys. Based on the shared information and the set of previous authorization information Y, CCSys can determine the cause of the task information mutation. If a task is completed, the auction is not started; if a new task $T_{new}$ is created and the task information changes unexpectedly, the auction is started and the previous set of preauthorization information X is deleted and the task price vector $\mathbf{P}\left(t\right)$ is updated;
  • In the case that the UAV $U_i$ is lost, the UAV $U_i$ task completion time is known from the elements in the previous set of authorization information from CCSys (expected arrival time $t_i^1$). If IPSys shares the information with CCSys at this point and it matches the expected change, the task is successfully executed and the auction is not started; otherwise, it means that the UAV is malfunctioning and there is one less bidder in the next auction, so ask IPSys to share the task information at this moment, start a new auction, and delete the previous set of pre-authorization information X and reset the task price vector $\mathbf{P}\left(t\right)$.
• CCSys passive start: in case of no interference.
  After the UAV $U_i$ arrives at the MDC$V_0$, it sends an auction asking to CCSys. If $U_i$ has been pre-authorized, it is directly converted to authorization and sent to $U_i$. Meanwhile, what needs to be achieved is to add the pre-authorization quintuple relation related to $U_i$ to the authorization set Y and remove it from the pre-authorization set X; conversely, start the auction and CCSys sends the task information set $TASK$ to all UAVs that are not pre-authorized.

3.2. Mechanism Design

In MTTDDVRP-SD, after a task is executed once, if it is not completely completed, it is considered as a new task and the gain obtained from the execution is recalculated. In a dynamic environment, interference new tasks generations and task unexpected changes occur without burden to the proposed algorithm. This is because newly generated tasks or tasks with unintended changes in information are treated the same as tasks after one execution. Moreover, in the auction mechanism, UAVs lose their capabilities due to failures, which is equivalent to a reduction in the number of auctioneers. Due to the loose coupling between UAVs, no additional burden is imposed on the algorithm of this paper.

The SAABoPA mechanism consists of four phases: (1) Task Publication Phase, (2) Bidding Phase, (3) Pre-authorization Phase and (4) Authorization Phase. The CCSys performs the task publication, pre-authorization and authorization phases, and the UAV in the LTSys performs the bidding phase.

3.2.1. Task Publication Phase

At the beginning of an auction, CCSys releases to LTSys the estimated set of task information $TASK=\{T,(x,y),R^{\prime },E^{\prime },\mathbf{P}\left(\tau \right)\}$, where T denotes the set of uncompleted tasks at the current $\tau$ moment. When task $T_j$ is partially completed, $T_j$ will not be deleted from T, but the information of material requirements and emergency indicator will be updated; when task $T_j$ is fully completed, $T_j$ will be deleted, $T=T-\left\{T_j\right\}$; when a new task $T_{new}$ appears, the set of tasks T is expanded, $T=T\bigcup \left\{T_{new}\right\}$. $N_T$ denotes the number of elements in the set T. $(x,y)=\left\{(x_j,y_j)\right|T_j\in T\}$ denotes the position information of the uncompleted tasks. $R^{\prime }=\left\{R_j^{{}^{\prime }}\left(\tau \right)\right|T_j\in T\}$ denotes the estimated material requirement of the uncompleted tasks. $E^{\prime }=\left\{E_j^{{}^{\prime }}\left(\tau \right)\right|T_j\in T\}$ denotes the estimated value of the emergency indicator of the uncompleted tasks. $\mathbf{P}\left(t\right)$ is the price vector consisting of the prices of all current outstanding tasks, and the initial moment price vector $\mathbf{P}\left(0\right)=0$.

In the r-th auction, there are some tasks that have been authorized but not executed in the $r-1$-th auction. In order to accurately calculate the rewards of the UAV for the tasks in the r-th auction, we introduce the material requirement information prediction $R_j^{{}^{\prime }}\left(\tau \right)$ and the emergency indicator information prediction $E_j^{{}^{\prime }}\left(\tau \right)$. The $R_j^{{}^{\prime }}\left(\tau \right)$ and $E_j^{{}^{\prime }}\left(\tau \right)$ of the task $T_j$ are calculated based on the $R_j\left(\tau \right)$ and $E_j\left(\tau \right)$ true values and the authorized UAV $U_{i_j}$ in the $r-1$ auction delivery $M{D}_{i^{j}j}\left(\tau \right)$, expressed as follows.

$$R_j^{{}^{\prime }}\left(\tau \right)=R_j\left(\tau \right)-M{D}_{i^{j}j}\left(\tau \right)$$

(19)

$$E_j^{{}^{\prime }}\left(\tau \right)=\left\{\begin{array}{cc}{E}_j\left(\tau \right)-\frac{E_j\left(\tau \right)}{R_j\left(\tau \right)}·M{D}_{i^{j}j}\left(\tau \right),\hfill & R_j^{{}^{\prime }}\left(\tau \right)\ne 0,\hfill \\ 0,\hfill & \mathrm{others}.\hfill \end{array}\right.$$

(20)

3.2.2. Bidding Phase

Different loadings of UAVs lead to different speeds, so that UAVs may obtain different rewards for the same task. Accordingly, while completing the task assignment in EMDS, the loading and delivery of UAVs can be determined. The pseudo-code is shown in Algorithm 1.

Before bidding, the UAVs in LTSys calculate the income values based on their own computational resources based on the set of task information shared by CCSys, and form the income matrix:

$$INCOM{E}_{i,T}\left(t\right)={\left[\begin{array}{ccc}Incom{e}_{i,1}^1\left(t\right)& \cdots & Incom{e}_{i,1}^{L_i}\left(t\right)\\ ⋮& \ddots & ⋮\\ Incom{e}_{i,\left|T\right|}^1\left(t\right)& \cdots & Incom{e}_{i,\left|T\right|}^{L_i}\left(t\right)\end{array}\right]}_{\left|T\right|\times L_i}$$

(21)

where, ${INCOME}_{i,T}\left(t\right)$ denotes the income matrix of the UAV $U_i$ for the current set of outstanding tasks T at t moments. The matrix elements are calculated as:

$$Incom{e}_{ij}^{L_i}\left(t\right)=Valu{e}_{ij}^{L_i}\left(t\right)-Cos{t}_{ij}^{L_i}\left(t\right)$$

(22)

$$Valu{e}_{ij}^{L_i}\left(t\right)=\left\{\begin{array}{cc}{E}_j\left(t_0\right)+a·(\frac{D_{ij}}{v_i\left(t\right)}+t-t_0)+\frac{E_j\left(t_0\right)}{R_j\left(t_0\right)}·max(M{D}_{ij},L_i),\hfill & E_j\left(t\right)+\frac{a·D_{ij}}{v_i\left(t\right)}<1,\hfill \\ 0,\hfill & \mathrm{others}.\hfill \end{array}\right.$$

(23)

$$Cos{t}_{ij}^{L_i}\left(t\right)=\left\{\begin{array}{cc}\frac{s·(t_i^1-t_i^0+-\frac{d_{j,V_0}}{v_i\left(t_{{}_i}^1\right)})}{T_i^{duration}-T_i^{fly}},\hfill & t_i^1-t_i^0+-\frac{d_{j,V_0}}{v_i\left(t_{{}_i}^1\right)}+T_i^{fly}\le T_i^{duration},\hfill \\ \infty ,\hfill & \mathrm{others}.\hfill \end{array}\right.$$

(24)

where, $Incom{e}_{ij}^{L_i}\left(t\right)$ denotes the reward value that a UAV $U_i$ with loading $L_i$ of materials can obtain by executing a task $T_j$ at time t, which is determined by the value of the task $Valu{e}_{ij}^{L_i}\left(t\right)$ and cost $Cos{t}_{ij}^{L_i}\left(t\right)$; s is the scaling factor; $t_i^0$ denotes the expected time to start the task $T_i$; $t_i^1$ denotes the expected time to complete the task $T_j$, satisfying $D_{ij}={\int }_t^{t_i^1}{v}_i\left(\tau \right)d\tau$; $T_i^{fly}$ denotes the time that the UAV $U_i$ has flown; $M{D}_{ij}$ denotes the time that the UAV $U_i$ has performed the task $T_j$ of material delivery, satisfying $M{D}_{ij}=min\{L_i,R_j\left(t\right)\}$; $D_{ij}$ denotes the actual flight distance of the UAV $U_i$ from the current position to the task $T_j$.

Algorithm 1: Bidding phase (Part of LTSys, take $U_i$ as an example).

In Equation (23), under the constraint of $\frac{E_j\left(t\right)}{a}+\frac{D_{ij}}{v_i\left(t\right)}<1$, the first two terms of the right-hand side of the $Valu{e}_{ij}^{L_i}\left(t\right)$ calculation equation indicate the emergency indicator at the moment before the task is executed, which reflects the basic value of the UAV to execute the task; the third term of the equation reflecting the additional value of the UAV $U_i$ executing task $T_j$, i.e., the UAV loads as much materials as possible under the condition that the material is obtained before the emergency indicator of the task $T_j$ reaches 1.

The UAV $U_i$ based on the price vector $\mathbf{P}\left(t\right)$ of the uncompleted tasks, combined with the income matrix can be calculated to obtain the net income matrix ${NET}_{i,T}\left(t\right)$ as:

$$NE{T}_{i,T}\left(t\right)=INCOM{E}_{i,T}\left(t\right)-{[{\mathbf{P}}_T\left(t\right);{\mathbf{P}}_T\left(t\right);...;{\mathbf{P}}_T\left(t\right)]}_{{}_{N_T\times L_i}}$$

(25)

The maximum net income $ne{t}_{i,T}^{max}\left(t\right)$ is:

$$ne{t}_{i,T}^{max}\left(t\right)=maxNE{T}_{i,T}\left(t\right)$$

(26)

Then, the tasks $T_{j^i}$ bid by the UAV $U_i$, the material load $L_i$ and the delivery ${MD}_{i{j}^i}$ to the task $T_{j^i}$ are:

$$(m,n)=index(ne{t}_{i,T}^{max}\left(t\right))$$

(27)

$$T_{j^i}=T\left(m\right)$$

(28)

$$L_i=n$$

(29)

$$M{D}_{i{j}^i}=min\{L_i,R_{j^i}\left(t\right)\}$$

(30)

where the function $(m,n)=index\left(a\right)$ in Equation (27) means to obtain the values of the ranks of a in the matrix A, and if there are multiple values of a in the matrix A, then a randomly selected rank value.

To further determine the bidding of the UAV $U_i$ to the task $T_{j^i}$, the maximum net income value $ne{t}_{i,T-\left\{T_{j^i}\right\}}^{max}\left(t\right)$ in the net income matrix of the UAV $U_i$ other than the task $T_{j^i}$ is also needed, expressed as follows:

$$ne{t}_{i,T-\left\{T_{j^i}\right\}}^{\max }\left(t\right)=max{NET}_{i,T-\left\{T_{j^i}\right\}}\left(t\right)$$

(31)

where $T-\left\{T_{j^i}\right\}$ denotes the deletion of elements of the set $T_{j^i}$ from the set of tasks T.

Combining Equations (21)∼(31), the bidding vector ${\mathbf{P}}_{\mathbf{i}}^{{}^{\prime }}\left(t\right)$ of the UAV $U_i$ for all tasks is determined as follows.

$${{\mathbf{P}}_i}^{{}^{\prime }}\left(t\right)={(p_{{}_{i,1}}^{{}^{\prime }}\left(t\right),p_{{}_{i,2}}^{{}^{\prime }}\left(t\right),...,p_{{}_{ij}}^{{}^{\prime }}\left(t\right),...,p_{{}_{i,N_T}}^{{}^{\prime }}\left(t\right))}^T$$

(32)

$$p_{ij}^{{}^{\prime }}\left(t\right)=\left\{\begin{array}{cc}Incom{e}_{i,j^i}^{L_i}\left(t\right)-ne{t}_{i,T-\left\{T_{j^i}\right\}}^{\max }\left(t\right)+\epsilon ,\hfill & j=j^i,\hfill \\ -\infty ,\hfill & j\ne j^i.\hfill \end{array}\right.$$

(33)

where $\epsilon >0$ is a slack parameter [42] that ensures that the algorithm does not enter a dead loop when multiple UAVs compete for the same task because the bids do not improve; $P_{ij}^{{}^{\prime }}\left(t\right)$ denotes the bidding of the UAV $U_i$ for the task $T_j$ at time t.

If the elements of the net income matrix ${NET}_{i,T}\left(t\right)$ are all non-positive, which means that the UAV $U_i$ will not profit from bidding on any task at the current price, then the UAV $U_i$ abandons further bidding in this auction and sets ${\mathbf{P}}_i^{{}^{\prime }}\left(t\right)={(-\infty ,-\infty ,\dots ,-\infty )}^T$, as shown in lines 3∼6 of Algorithm 1; lines 7∼13 show that if there are positive values of the elements in the net income matrix, the UAV $U_i$ participates in the bid, and the task, loading and bid price of the bid are determined by Equations (26)∼(33).

After the UAV $U_i$ in LTSys completes the bidding phase, it sends the set of bid information $BI{D}_i=\{U_i,T_{j^i},{\mathbf{P}}_i^{{}^{\prime }},L_i,M{D}_{ij},t_i^1\}$. The elements in the set ${BID}_i$ denote the UAV order number, the bid task, the bid vector to the task, the material loading, the delivery to the task, and the expected arrival time, respectively.

3.2.3. Pre-Authorization Phase

The CCSys receives the bid information from the UAV and starts the pre-authorization phase, the pseudo-code is shown in Algorithm 2.

As in lines 1∼3 of Algorithm 2, CCSys waits to receive the bid information $BI{D}_i$ and then updates its own stored bidding matrix $P^{\prime }\left(t\right)$ based on the bidding vector in the bid information, expressed as follows.

$$P^{{}^{\prime }}\left(t\right)={\left[\begin{array}{ccc}{P}_{1,1}^{{}^{\prime }}\left(t\right)& \cdots & P_{\left|U\right|,1}^{{}^{\prime }}\left(t\right)\\ ⋮& \ddots & ⋮\\ P_{1,\left|T\right|}^{{}^{\prime }}\left(t\right)& \cdots & P_{\left|U\right|,\left|T\right|}^{{}^{\prime }}\left(t\right)\end{array}\right]}_{\left|T\right|\times \left|U\right|}$$

(34)

Pre-authorize the task $T_j$ to the UAV with the highest bid $U_{i^j}$,

$$U_{i^j}=\underset{U_i\in U}{argmax}\left\{P_{1,j}^{{}^{\prime }}\left(t\right),P_{2,j}^{{}^{\prime }}\left(t\right),...,P_{\left|U\right|,j}^{{}^{\prime }}\left(t\right),\right\}$$

(35)

We consider removing the quintet of relations $(U_i,L_i,{MD}_{ij},t_i^1,T_j)$ associated with the task $T_j$ from the set of pre-authorized information X, and adding new quintet of relations $(U_{i^j},L_{i^j},{MD}_{i^{j}j},t_{i^j}^1,T_j)$, as shown in lines 4∼8.

Where line 7 indicates that once the task $T_j$ is pre-authorized to the UAV $U_{i^j}$ at moment t, the price vector $\mathbf{P}\left(t\right)$ of all uncompleted tasks is updated, i.e., the price of the task $T_j$ is updated to the bidding price of the UAV $U_{i^j}$,

$$p_j\left(t\right)=P_{i^j,j}^{{}^{\prime }}\left(t\right)$$

(36)

As in lines 9∼13 of Algorithm 2, at the end of an auction round, CCSys sends a rebid command to the UAVs $U_i$ that have not abandoned their bids and are not pre-authorized, and sends a set of current information about all tasks $TASK$. The UAVs that receive the rebidding command start a new round of the auction process, while the other UAVs enter the waiting phase. Lines 14∼15 indicate that the iterative process terminates with the condition that all UAVs that have not abandoned their bids during the auction process are pre-authorized for the task.

Algorithm 2: Pre-authorization phase (CCSys part).

3.2.4. Authorization Phase

CCSys updates the set of authorization information Y based on the set of pre-authorization information X, adds the pre-authorization quintet relation to Y for UAVs $U_i$ that are pre-authorized but not authorized, and removes it in X, while resetting the task price vector. CCSys sends the authorization quintet $y_i$ to the UAVs authorized in this auction, sends the loading information $U_i$, $L_i$ to the MDC, and saves the set of pre-authorization information and the set of authorization information Y as the initial state for the next auction.

3.3. Algorithm Analysis

The algorithm involves internal system interaction and inevitably has communication problems. On the one hand, in the disaster area environment, ground communication base stations are damaged, and it is necessary to restore communication, and UAVs are mature and widely used as temporary communication relays [3]; on the other hand, the bidding process is performed by UAVs, which greatly reduces the state information interacted between CCSys and UAVs, and can be used in a limited bandwidth environment [27].

The computation time of one auction in the SAABoPA algorithm is much less than that of the traditional auction algorithm [25], although the former requires multiple executions of the algorithm during the completion of all tasks, resulting in a longer cumulative execution time. However, considering the case of excessive interference, the traditional auction algorithm needs to reset or partially reset the previous allocation results and reallocate them. If we consider the case of poor communication, the conventional auction algorithm may cause excessive waiting time due to waiting for communication in a single execution, which may affect the timely departure of UAVs [28], and thus affect the capability performance and resilience performance of EMDS. In the EMDS framework, a single auction execution can be performed during the UAV’s return trip or by using another gap time. As argued in the literature [45], real-time does not mean that the faster the computation time is, the more real-time it is, but that the solution to the problem is found within a limited time.

Specifically, we analyze the time complexity and space complexity of one auction round in a single auction of SAABoPA. As mentioned above, an auction is divided into the task publication phase, bidding phase, pre-authorization phase, and authorization phase, and the main complexity is concentrated in the bidding phase and pre-authorization phase. However, in the task publication phase and authorization phase, the main work is to complete the packaging, sending, and assignment of data processing operations. In the bidding phase, see (Algorithm 1), the main work is to complete the computation of the $INCOM{E}_{i,T}\left(t\right)$ matrix, and the sorting operation to find the maximum element of the matrix. Assuming that the time complexity of computing a single matrix element is I and the space complexity is J, then the average time complexity is $O\left(\right|T|\ast L_i\ast I)$, and the space complexity is $O\left(\right|T|\ast L_i\ast J)$. On the other hand, the sorting algorithms are heap sorting, bubble sorting, and shell sorting, whose time complexity is $O(log(\left|T\right|\ast L_i\left)\right)$, $O\left(\left(\right|T|\ast L_i{)}^2\right)$, and $O\left(\left(\right|T|\ast L_i{)}^{1.3}\right)$, respectively, and the space complexity is is $O\left(1\right)$. In the Pre-authorization Phase (Algorithm 2), the main work is to complete the sorting operation to find the maximum element of the bidding vector $BID$ for each task, similarly using heap sorting, bubble sorting, shell sorting, their time complexity is $O\left(\right|T|\ast log\left(L_i\right))$, $O\left(\left|T\right|\ast {\left(L_i\right)}^2\right)$ and $O\left(\left|T\right|\ast {\left(L_i\right)}^{1.3}\right)$, respectively, and the space complexity is $O\left(1\right)$.

4. Computational Experiment

In this section, the effectiveness of the proposed method is verified by experiments. The simulation experiments were performed using Visual Studio 2022 platform, C++ programming language, CPU is Intel(R) Core(TM) i9-10885H CPU @ 2.40GHz 2.40GHz, OS is Windows 10.

4.1. Experimental Settings

4.1.1. Case Scenario

The relevant data for the computational case were chosen using a combination of real and hypothetical data, since some disaster data is not published or could not be obtained through official reports [46,47]. The data of the UAVs are derived from some public sources and scaled accordingly to fit the case scenario. Meanwhile, in order to verify the performance of the algorithm, large sample experiments are required in this section, and the data for the experimental scenarios are randomly generated on the basis of known partial data.

The earthquake disaster in a city in southwest China is simulated as a scenario, and the affected area C is a rectangular area of 4000 × 4000 m, with $\left|T\right|$ tasks uniformly distributed in the area, and the quantity of material required for the tasks $R\left(0\right)$ and the emergency indicator $E\left(0\right)$ also obey uniform distribution. The outside of the rectangular area is secure, and a well-stocked MDC is deployed on the boundary of the disaster area with a total of $\left|U\right|$ UAVs involved in the task execution. The number of the three dynamic interference types considered during the task execution is determined in advance as $|DI{S}^1|$, $|DI{S}^2|$, and $|DI{S}^3|$, but the time of occurrence of the dynamic interference types is random and uncertain, and the values of the changed information obey a uniform distribution. The specific parameters are shown in

Table 1. Scenario parameters instructions.

Parameters	Value	Parameters	Value
$x,y$	$U(0,4000)$	$(x_0,y_0)$	Boundary of C
$R\left(0\right)$	$DiscreteU(6,10)$	$L^{max}$	$U(11,15)$
$E\left(0\right)$	$U(0.1,0.8)$	$V^{unload}$	$U(15,20)$
$\Delta R$	$DiscreteU(0,5)$	$T^{duration}$	$U(400,500)$
$\Delta E$	$U(0,0.4)$	$\delta$	0.5

.

In order to investigate the effects of different dynamic interference intensities on the performance of the algorithm, this paper designs experiments with seven different control variables at different scenario scales, as shown in

Table 2. Scenario scale instructions.

Name of Scenario	$\left|\mathit{T}\right|\times \left|\mathit{U}\right|\times |\mathbf{DI}{\mathit{S}}^1|\times |\mathbf{DI}{\mathit{S}}^2|\times |\mathbf{DI}{\mathit{S}}^3|$1
Scenario 1	$50\times 5\times 10\times 10\times 2$
Scenario 2	$50\times 5\times 20\times 10\times 2$
Scenario 3	$50\times 5\times 30\times 10\times 2$
Scenario 4	$50\times 5\times 20\times 20\times 2$
Scenario 5	$50\times 5\times 20\times 30\times 2$
Scenario 6	$50\times 5\times 20\times 10\times 3$
Scenario 7	$50\times 5\times 20\times 10\times 4$

¹ Taking Scenario 2 as an example, 50 × 5 × 20 × 10 × 2 means that there are 50 tasks in the scenario, 5 UAVs, and the number of three dynamic disturbances new task generations, task unexpected changes and UAV ’s number decreases are 20, 10 and 2, respectively.

.

4.1.2. Comparison Algorithms

The experiments with 1000 random samples are conducted for each scenario, respectively, considering all three interference types, and the SAABoPA algorithm under the EMDS framework is compared with the DTAP (DTA based on sequential single item auctions) [28] and CBBA-PR (consensus-based bundle algorithms-partial replanning) [26] algorithms to verify and compare the capability performance and resilience performance of EMDS under the three algorithms.

Numbered lists can be added as follows:

• DTAP: Dynamic task assignment algorithm based on sequential single-item auctions is an online algorithm that differs from the algorithm in this paper in that there is no pre-authorization phase, and the DTAP algorithm starts once the auction timing is the current target task completion.
• CBBA-PR: Dynamic task assignment based on the consistent consent bundle algorithm is an offline algorithm where one auction UAV gets one task sequence, i.e., one task bundle. When interference occurs, a certain percentage of unexecuted tasks with low revenue are selected for authorization reset and the auction is restarted.

4.2. Performance Indicators

To verify the performance of EMDS based on the SAABoPA algorithm, the performance metrics fully consider the trade-off between efficiency and resilience [48]. In this paper, two main types of performance metrics are used: (a) capability performance, i.e., the completion of the rescue mission; and (b) resilience performance, i.e., the change of the system capability performance in the case of dynamic interference.

4.2.1. Capability Performance ($CP$)

In Definition 5, it is considered that when the emergency indicator reaches 1, the task will not be completed for a long time and the actual rescue will not be timely and cause secondary injuries. In this paper, the system capability performance is described as the proportion of tasks that are not subject to secondary injuries, i.e., the proportion of tasks whose emergency indicator fails to reach 1, which, combined with Equation (6), is expressed as follows.

$$CP=\frac{N-f}{N}\times 100\%$$

(37)

where N denotes the total number of tasks that occur during the entire EMDS emergency material delivery.

4.2.2. Resilience Performance ($RP$)

The EMDS capability performance metric proposed in this paper is time-invariant, not computable at any moment, and is obtained only at the end of the whole working process. In order to describe the calculation of the resilience performance metric, this paper considers the system capability performance as time-invariant and represents the in-process capability performance by the calculated value $CP$. When an interference situation occurs, for example, $DI{S}^1$ occurs at $t_{{}_1}^s$, the system capability performance after $t_{{}_1}^s$ moment is denoted by the computed value $CP$; the computed value $CP$ in the same initial situation, under the condition that $DI{S}^1$ does not occur at $t_{{}_1}^s$ moment denotes the system capacity performance before $t_{{}_1}^s$ moments. The same is true for the case of multiple interferences. The resulting EMDS capability performance metric is converted to time-dependent and the system capability performance changes as shown in Figure 3.

Since the capability performance at any moment of the whole working process is expressed by the calculated value $CP$, this paper does not describe too much the change of the system capability performance in a short time after the interference occurs, as Figure 3 incorporates the capability degradation, maintenance and recovery phases of the system after suffering from the interference [49].

Combined with the integral elasticity model, this paper defines the resilience performance of the system expressed as follows.

$$RP\left(\tau \right|DIS)=\frac{{\int }_{t_0}^{t_{end}}CP\left(\tau \right|DIS)d\tau }{{\int }_{t_0}^{t_{end}}CP\left(\tau \right)d\tau }$$

(38)

where $RP\left(\tau \right|DIS)$ denotes the system resilience performance with interference $DIS$ occurring, $CP\left(\tau \right)$ denotes the system baseline capability performance when $\tau$ without interference, and $CP\left(\tau \right|DIS)$ denotes the system capability performance when $\tau$ with interference.

4.3. Algorithm Results of Hypothesis Testing: Friedman Test

Based on the experimental design of 1000 random instances of 7 scenarios in the previous paper, the three algorithms can obtain 7 sets of experimental results. We compare the computational results of the three algorithms pair-wisely for each group of results, from three perspectives of no interference capability performance, interference capability performance, and resilience performance, to test the performance of the SAABoPA algorithm proposed in this paper by whether there are significant differences. Current experiments on comparing the average performance of multiple algorithms on multiple problem instances are widely and well-known by methods such as Complete Block Design (CBD) ANOVA or Friedman’s test [50]. In this paper, Friedman’s test is used, but the Friedman’s test can only test whether the average performance indicator value of at least one algorithm has a significant difference from the average performance indicator value of at least one algorithm, and cannot specifically determine which two algorithms have a significant difference between them, so Multcompare is also required. The results are shown in Figure A1. As an additional note, the confidence interval of Friedman’s test in this paper is 99%.

From the results, in the seven scenarios, there are no significant differences in the capability performance of the system without interference under the SAABoPA and CBBA-PR algorithms, but there are significant differences with the DTAP algorithm. However, there are significant differences in the capability performance of the system with interferences under the three algorithms. It indicates that the SAABoPA and CBBA-PR algorithms have similar capability performance in the no interference case and are significantly better than the DTAP algorithm. However, in the case of interferences, the capability performance of the system under the CBBA-PR algorithm degrades more and is significantly worse than SAABoPA and DTAP algorithms. Analyzing Figure A1c,f,i,l,o,r,u together, it can be found that there are significant differences in the resilience performance of the system under the three algorithms. The exception is that in Figure A1f,i (scenarios 2 and 3), there is no significant difference between SAABoPA and DTAP algorithms regarding the system resilience performance, which may be because the resilience performance of them gradually approaches as the number of new tasks generation and the intensity of interference faced by the system increases. However, compared to scenarios 4 and 5, 6 and 7, the system is not yet facing the maximum intensity of interference, but the type of interference has changed, which reflects the greater adaptability of the SAABoPA algorithm to new task generations type of dynamic interferences. On the other hand, the specific performance gap also needs to be analyzed in terms of average, standard deviation, CDF, etc.

4.4. Algorithm Comparison Analysis: Scenario 2 as an Example

At the scale of Scenario 2, 1000 random sample experiments were conducted. The results of the simulation experiments are studied from the mean and variance of the metrics, as shown in

Table 3. Average of system performance with different algorithms for 1000 sample experiments.

Algorithm	$\mathit{CP}\left(\mathit{\tau }\right)$	$\mathit{CP}\left(\mathit{\tau }\right|\mathit{DIS})$	$\mathit{RP}\left(\mathit{t}\right|\mathit{DIS})$
SAABoPA	0.9491	0.7992	0.9443
DTAP	0.8867	0.7506	0.9459
CBBA-PR	0.9497	0.6663	0.8916

and

Table 4. Variance of system performance with different algorithms for 1000 sample experiments.

Algorithm	$\mathit{CP}\left(\mathit{\tau }\right)$	$\mathit{CP}\left(\mathit{\tau }\right|\mathit{DIS})$	$\mathit{RP}\left(\mathit{t}\right|\mathit{DIS})$
SAABoPA	0.0032	0.0102	0.0009
DTAP	0.0055	0.0115	0.0005
CBBA-PR	0.0031	0.0132	0.0013

. The SAABoPA has a slightly lower capability performance than the CBBA-PR and roughly 6% higher than the DTAP in the static case where no interference occurs, because the DTAP, as an online algorithm, considers only the current deterministic information when making decisions, while the SAABoPA adds a pre-authorization process that considers the current information along with the next task estimation information, and achieved better results. However, it is obvious that the CBBA-PR, as an offline algorithm, considers all the information of the whole system working process, so the derived allocation scheme is better. It should be noticed that this algorithm requires higher accuracy of information, and in the presence of errors in the model, the errors in the offline solution will gradually accumulate causing the degradation of the solution quality.The resilience performance of SAABoPA and DTAP is similar, which is more than 5% higher than that of CBBA-PR, reflecting that the former two have better adaptability to the dynamic environment and resistance to interference situations are stronger. After analysis, the reason is that when there is an interference situation, the CBBA-PR needs to reallocate some of the uncompleted tasks with the interference in a comprehensive consideration, and considering the trade-off between algorithm execution time and optimality, it often cannot handle the interference situation in time, and it causes a significant decrease in the system capability performance, which reflects the lack of system resilience performance. Analyzing the mean value of system capability performance and the mean value of resilience performance under the three algorithms, although the mean value of resilience performance is above 0.85, the resilience indicator integrally considers the time of occurrence of multiple interferences, so the change of the calculated value of the indicator is relatively insignificant, but the capability performance after corresponding to multiple interferences will have a large reduction, for example, the highest resilient performance DTAP with multiple interferences and the ability to perform without interference differ by about 14%.

Analyzing the system performance variance, it is found that the variance of both metrics is around 0.01, which reflects the stability of the capacity performance and resilience performance of EMDS under the three algorithms in the scale of Scenario 2, and the experimental results are effective.

A comparative analysis of the cumulative distribution functions (CDF) of the capability performance and resilience performance of SAABoPA, DTAP and CBBA-PR (Figure 4) shows that: (a) the CDF curves of the capability performance of the SAABoPA and CBBA-PR almost overlap in the initial case where no interference occurs, and the percentage of capability performance greater than 0.9 is 80% much higher than the 44% of the DTAP; (b) the capability performance of all three algorithms decreases in the presence of multiple disturbances. However, the SAABoPA maintains the best capability performance with 50% of capability performance greater than 0.8 higher than 31% of DTAP and 10% of CBBA-PR; (c) the CDF curves of the resilience performance of SAABoPA and DTAP are close, and the percentage of resilience performance greater than 0.9 is over 90% much higher than 41% of CBBA-PR algorithm. Moreover, the lower limit of the resilience performance of SAABoPA and DTAP still exceeds 0.8, which indicates that SAABoPA and DTAP have stronger resilience performance and better adaptability in the face of interference situations.

For further analysis, the performance metrics of each sample were compared as shown in Figure 5. To facilitate the observation of the experimental results, we adjusted the sample serial numbers and reordered the samples in the lowest to highest order of SAABoPA’s metrics. We found that although the experimental parameters are the same, the disaster area environment and UAV parameters (e.g., mission coordinates, UAV speed, etc.) are different for different samples, resulting in different performance metrics for different samples by different algorithms, which we believe is reasonable, and there is a certain applicability of different algorithms to different scenarios themselves, which is why we consider large sample experiments.

4.5. Analysis of the Results under Different Dynamic Interference Intensity

As for the seven scenarios shown in Table 2, they are divided into three groups (scenario 1, 2, 3, scenario 2, 4, 5 and scenario 2, 6, 7) according to the principle of controlled variable method in this paper to study the effects of new task generations, task unexpected changes and UAV ’s number decreases at different intensities on the performance of EMDS, as shown in Figure 6, Figure 7 and Figure 8. From the experimental results, the more new task generations type of dynamic interferences have a certain degradation on the capability performance and resilience performance of all three algorithms. The capability performance of the SAABoPA proposed in this paper is optimal under different intensities of dynamic interference. On the other hand, when the number of new tasks increases from 10 to 30 (from Scenario 1 to Scenario 3), the capability performance of SAABoPA and DTAP decreases by about 10% and the resilience performance decreases by about 4%, but the CBBA-PA decreases by about 22% and 8%, respectively, which reflects the importance of the resilience performance indicator in the scenario with dynamic interference, which can better describe the performance of EMDS. task unexpected changes and UAV’s number decreases type of dynamic interference at different intensities can still obtain the same conclusion as new task generations.

In the vertical comparison, compared with Scenario 2, the number of new tasks generations and task unexpected changes increases by 10 in Scenario 3 and Scenario 4, respectively, but the capability performance of the three algorithms decreases by approximately 7%, 7% and 13% in Scenario 3, and the resilience performance decreases by 3%, 3% and 5%, respectively; then the capability performance of the three algorithms decreases by approximately 4%, 2% and 7% in Scenario 4, and the resilience performance decreases by 1%, 1% and 4%, respectively. In Scenario 4, the capability performance of the three algorithms decreases by about 4%, 2% and 7%, and the resilience performance decreases by 1%, 1% and 4%, respectively. This reflects the fact that dynamic interference of the new tasks generations type has a greater impact on the system capability performance and resilience performance than task unexpected changes, which is reasonable. However, since UAV’s number decreases type of dynamic interference is the reduction in task execution vehicles, it is not comparable to new tasks generations and task unexpected changes.

To analyze the dynamic interference of UAV’s number decreases type individually, from scenario 2 to scenario 6 to scenario 7, the UAVs are reduced by 1 in turn, which means that the transport delivery capacity is reduced by 20% in turn (total number = 5), while the capability performance reductions of the three algorithms become larger from 8%, 7%, and 9%, respectively, with reductions of 13%, 19%, 12%; on the other hand, the reduction in resilience performance becomes 4%, 4%, 4% from 3%, 2%, 3%, respectively. This indicates that the smaller the number of available UAVs, the worse the system performance, and this gap is increasing as the number of UAVs decreases. It can also be found that the DTAP has a more drastic performance degradation compared to the SAABoPA and CBBA-PA algorithms in scenario 7 (loss of 4 UAVs), and the reason for this analysis is that the capability performance of EMDS under the DTAP is not high when no such interference occurs, and the more UAVs are lost to the DTAP brings about a further degradation of the task planning capability.

Comparing the standard deviations of the experimental results under each experimental scenario, it can be found that the resilience indicator $RP$ is smaller compared to $CP$ under no interference and $CP$ under multiple interferences, reflecting that the resilience indicator of multiple interferences proposed in this paper is more stable under large sample experiments, and it is robust and convincing to describe the system resilience performance by $RP$.

In summary, the capability performance of the SAABoPA proposed in this paper is similar to that of CBBA-PA under no interference and better than that of the DTAP; the capability performance of the SAABoPA under multiple interference is better than that of the DTAP and CBBA-PA; the resilience performance of the SAABoPA under multiple interference is similar to that of the DTAP and better than that of the CBBA-PA, and it is more adaptable to dynamic scenarios.

5. Conclusions

This paper presents a basic framework for applying UAVs as important components of EMDS to emergency material delivery scenarios, and proposes SAABoPA for capability performance and resilience performance optimization. Combined with the actual material delivery, the problem is modeled as MTTDDVRP-SD on the basis of the VRP problem, while several features are considered, including multi-trip per UAV, time-dependent emergency indicator, dynamics of the environment, split-delivery, and trip duration limit. Three types of dynamic interference scenarios are considered, including: new task generation, task unexpected changes, and UAV’s number decreases. To solve this new problem, this paper proposes the SAABoPA algorithm, which adds a pre-authorization and sequential auction mechanism. Meanwhile, in order to describe the effect of dynamic disturbances on EMDS, we improve the integral resilience model under single interference and propose a resilience performance indicator adapted to multiple interference scenarios. Finally, through comparison experiments with DTAP and CBBA-PR in scenarios with different dynamic intensities, and using Friedman’s test with 99% confidence interval, we verified the effectiveness of SAABoPA, which can improve the resilience performance of the system while maintaining a better capability performance, compensating for the shortcomings of CBBA-PR which only focuses on capability performance and DTAP algorithm which only focuses on resilience performance.

It should be noted that the scenario in this paper is material delivery under disaster area rescue, which mainly considers task survival measured by task emergency indicator as the optimization objective, and may need to trade-off more objectives if it is applied to other scenarios such as UAV logistics, combat equipment delivery, and cluster fire strike. In the future, we consider building a small hardware experiment bed to test the algorithm effectiveness and system performance, and try to solve more practical problems to extend to real applications.