Auburn University National Geotechnical Experimentation Site in Piedmont Residuum

1.   Introduction

In this paper, we consider the planning of an unmanned aerial vehicle (UAV) fleet mission problem with highly dynamic and unpredictable environment constraints ^[1,2,3,4,5]. Typical disruptions in urban deliveries by UAVs (characterized by frequent stop-and-go movements, low consolidation, and frequent re-scheduling), may be caused by changes in the order by customers, or shifting weather conditions (e.g., a sharp drop in temperature, icing of propellers, turbulence), which affect the energy consumption of UAVs and cause their shorter range due to the depletion of batteries ^[6,7]. For that reason, the routing of a UAV fleet in a partially known and unpredictable environment should guarantee a reactive on-line determined contingency reaction.

The UAVs' routes can be determined in the course of proactive planning ^[8] and generated in offline mode or in reactive planning ^[9,10], then executed in the online mode. Routes determined in proactive planning guarantee the achievement of the planned mission's goal for the environment's parameters, which change at predetermined intervals. Scenarios corresponding to planned reactive rules do not guarantee the existence of reactive end-to-end paths employed in the routing process while adapting it to changes in the environment during mission execution ^{[2,6,10,11,12,13,14,15,16,17,18]}. In particular, the reactive routing strategies linking the "route discovery (proactive route planning)" and "route maintenance (reactive rules adopting)" concepts ^[19] are responsible for UAV robustness to the changes that appear in the urban distribution context. In this context, our study focuses on reactive planning ^[20,21] of deliveries by a UAV fleet, which are resistant to sudden changes in weather conditions and unforeseen changes in the delivery schedules.

The originality of the paper results from its merging of the proactive and reactive planning of missions of UAV fleet. The developed models allow for predictive (i.e., taking into account forecasted weather conditions changing) and reactive (i.e., enabling interruption of a drone's mission) planning of delivery missions in terms of the constraint satisfaction problem ^[6,22,23] and the easy implementation of the proposed genetic algorithm in a commercially available constraint programming environment, e.g., IBM ILOG.

The paper is structured as follows. Section 2 presents the state of art. The approach to UAV mission contingency planning is presented in Section 3. The declarative model for reactive planning of deliveries by UAV fleets according to the constraint satisfaction problem is described in Section 4. Sufficient conditions guaranteeing the existence of a non-empty space of admissible solutions and, consequently, a reduction in the time expenditure incurred on the assessment of the weather changes are outlined in Section 5. The structure and operation of the developed genetic algorithm are presented in Section 6. The results of the conducted experiments regarding implementation of an accurate and approximate approach to UAV fleet mission planning are discussed in Section 7. Final conclusions are stated in Section 8, followed by suggestions for future research.

2.   Related work

UAV mission planning issues play a key role in supply management systems in distributed networks belonging to the VRP class of NP-hard computational complexity ^[2,24]. Consequently, the sub-problems of routing and scheduling occurring in them also encounter a combinatorial explosion barrier, forcing the use of approximate solution methods. The goal is to identify a set of routes for a fleet of vehicles that optimizes some underlying objective function while subject to problem-specific constraints ^{[25,26,27,28,29]}. These efforts, according to the current taxonomy of available route planning methods, can be partitioned into two main categories: approximation and exact methods ^{[6,26,30,31,32,33,34]}.

In general planning of the UAV path comes down to finding the path connecting the starting point with the destination that guarantees the fulfillment of the UAV operational requirements within the specified UAV flight constraints. Problems of this type undertaken in the area of UAV routing, i.e., VRP class problems ^[1,5], are formulated for various assumptions covering both the topography and the presence of obstacles. The topography, related to, e.g., mountainous and/or urbanized areas, influences the changes of the flight trajectory level, while the presence of stationary and/or dynamic obstacles makes it impossible to keep the straight shape of the trajectory ^[35]. Limitations of this type, force the inclusion of many new variables multiplying the complexity of UAVs trajectory planning problems, thus implying the need to use metaheuristic algorithms. The underlying nature-inspired algorithms such as the Bath algorithm, particle swarm optimization (PSO), artificial bee colony (ABC) algorithm, ant colony algorithm (ACO), grey wolf algorithm, etc. compete in terms of their convergence accuracy and proceeding efficiency ^[36,37,38].

Many works dealing with the problem of UAV routing, in particular those related to the use of UAVs in situations related to the monitoring of reservoirs, fields and forests, as well as flood victims support, assume UAV moves over flat terrain without obstacles that exclude the linear shape of their trajectory. In situations of this type, assuming that the UAV takes off vertically with the collected goods and also lands vertically in order to unload it, the problem of UAV routing boils down to determining the appropriate wind correction angle, correcting its drift ^[24]. Assuming these assumptions, this work addresses the problem of re-planning a UAV route in situations where sudden changes in wind direction and its intensity go beyond the fluctuation ranges adopted at the stage of proactive route planning.

An approximation method attempts to find the best possible solution while providing no guarantees as to the quality of that solution, i.e., it prefers quick solutions over optimal ones. In turn, an exact method is guaranteed to find the best solution to the problem given enough time, i.e., it is oriented toward the search for the optimal solution at the expense of the time incurred to obtain it. In the approximation approach, there are methods that implement heuristic algorithms (e.g., metaheuristics-driven ones, such as variable neighborhood search (VNS), simulated annealing (SA), and tabu search (TS)); population algorithms, such as PSO, ACO, ABC and evolutionary algorithms, such as memetic algorithms (MMA), and genetic algorithm (GA)) ^{[26,30,32,34]}.

The exact approach leverages intelligent forms of enumerative search, such as dynamic programming (DP), mixed integer linear programming (MILP), branch-and-bound (BB), and constraint programming (CP) ^[6,31,33]. All of the above-mentioned approaches are supplemented by computer simulation tools, such as flight simulators ^[37,39,40].

There is a growing trend in research aimed at developing an interactive tool dedicated to online planning missions carried out by UAV teams, as reflected in numerous publications ^{[12,13,16,22]}. The aim of this research is to devise solutions that enable proactive and/or reactive planning of UAV missions in situations caused, for instance, by contingency planning in dynamic environments, i.e., decision support tools facilitating analysis and different scenarios comparison. Besides emergency situations, such tools can be helpful in various fields of application, e.g., reconnaissance and mapping, package delivery and delivery communication capabilities, healthcare and so on ^{[3,4,12,13,18]}.

The UAV mission planning problems encountered in practice are characterized by a large variety of parameters, constraints, and specific evaluation criteria. Some examples of possible differences are listed in the following exemplary sets:

●  parameters specifying UAVs (e.g., maximum loading capacity, maximum energy capacity, empty weight, payload weigh, ground speed, and aerodynamic drag coefficient), environmental conditions (e.g., weather forecasts including wind speed and direction, air density, temperature drizzle, or snowfall), and distribution network topology (e.g., number of delivery points, number of depots, and distances between delivery points and depots);

●  constraints determining UAVs trajectories (e.g., to avoid collisions, a UAV is configured with at least one flight corridor and flight path), battery durability, delivery volumes, delivery periods, and changing weather conditions;

●  evaluation criteria used to determine UAV fleet routings that guarantee the maximum benefits with the shortest total path, in particular, the maximum travel distance (limited by the capacity of the battery used), the minimum length of the sum of the flight path of each UAV, timely deliveries, resistance to changes in weather conditions, and/or service delivery deadlines.

As already mentioned, problems of this type belong to the category of NP-hard problems. This means that instances of such problems encountered in practice (due to their size as well as the number and variety of characteristics describing them) can be resolved within a reasonable time interval using approximate methods. In that context, the exact (i.e., optimal) solutions achieved in a similar time frame by exact methods refer, however, only to the small size instances out of practical significance. It is also easy to see that, in general, it is difficult to say which approach is best used for a given problem instance. The literature is rather scarce on this topic. Therefore, the research gap in this field is the motivation behind our research aimed at the implementation of selected methods representing both the approximate and exact approaches to compare them on using an arbitrarily selected example of a UAV fleet planning problem.

The comparison included CP as representative of the exact approach and GA as representative of the approximate approach. CP is a paradigm for solving combinatorial problems according to which users declaratively state the constraints on the feasible solutions for a set of decision variables. The acceptable solution sought is the set of decision variables values satisfying all the constraints ^[6,22,27]. It is worth emphasizing that CP is more general than MILP, allowing variable types beyond integer and continuous (e.g., interval and set variables), and dropping the restriction of linearity in the constraints and objective function. In turn, the GA, being the representative of evolutionary algorithms, is a search technique inspired by evolutionary biology such as inheritance, mutation, selection, and crossover used in computing to find true or approximate solutions to optimization and search problems ^[41,42,43]. Carrying out appropriate computer experiments allows to judge which approach is faster for which size of the selected instance of the UAV mission planning problem.

3.   Approach to UAVs fleet online routing

We consider a distribution network, which is modeled by the graph $G = \left(N, E\right)$ where $N = \left\{{N}_{1}, \dots, {N}_{\lambda }, \dots, {N}_{n}\right\}$ signifies the set of $n = \left|N\right|$ nodes (distinguishing ${N}_{1}$ node representing a base and $\left\{{N}_{2}, \dots, {N}_{n}\right\}$ nodes representing delivery points), and $E = \left\{\left({N}_{i}, {N}_{j}\right)|i, j\in \{1, \dots, n\}, i\ne j\right\}$ signifies the set of edges determining the possible connections between nodes. Given a fleet of UAVs $\mathcal{U} = \left\{{U}_{1}, \dots, {U}_{k}, \dots, {U}_{K}\right\}$ that delivers to the points {${N}_{2}, \dots, {N}_{n}$}, to each delivery point ${N}_{\lambda }$ an ordered quantity of goods ${z}_{\lambda }\in \mathbb{N}$ (taken from the base ${N}_{1}$) should be transported. Deliveries are made as part of the mission $S$, which consists of sub-missions ${}^{l}S$ (i.e., delivery plans that include a single course of UAVs: base-delivery points-base). $Z$ denotes a sequence consisting of variables ${z}_{\lambda }$ : $Z$ $= ({z}_{1}, \dots, {z}_{n}$). It is assumed that all required goods should be delivered in the given horizon time $H$. The amount of goods delivered during one sub-mission ${}^{l}S$ by the ${U}_{k}$ to the delivery point ${N}_{\lambda }$ is determined by the variable ${{}^{l}c}_{\lambda }^{k}\in \mathbb{N}$. ${}^{l}C$ is a sequence $:$ ${}^{l}C = ({{}^{l}c}_{1}^{1}, \dots, {{}^{l}c}_{1}^{K}, \dots, {{}^{l}c}_{n}^{1}, \dots, {{}^{l}c}_{n}^{K})$ determining the payload weight delivered by fleet $\mathcal{U}$. Variable ${Q}_{k}$ denotes the payload capacity of ${U}_{k}$ (amount of goods transported by ${U}_{k}$ cannot exceed ${Q}_{k}$). Moreover, each ${U}_{k}$ is described by the following technical parameters: battery capacity $CAP$, airspeed $va$, drag coefficient ${C}_{D}$, front surface $A$, and UAV width $b$. Time spent on take-off and landing ${U}_{k}$ on delivery point ${N}_{\lambda }$ is indicated by variable ${w}_{\lambda }\in \mathbb{N}$. Note that ${}^{l}\mathcal{U}\subseteq \mathcal{U}$ denotes a set of UAVs used during sub-mission ${}^{l}S$. The moment when the ${U}_{k}\in {}^{l}\mathcal{U}$ arrives to the delivery point ${N}_{\lambda }$ during sub-mission ${}^{l}S$ is indicated by variable ${{}^{l}y}_{\lambda }^{k}\in \mathbb{N}$ [s]. In that context, the sequence ${}^{l}Y$ consisting of moments ${{}^{l}y}_{\lambda }^{k}$, is called the schedule of the fleet ${}^{l}\mathcal{U}$ : ${}^{l}Y = ({{}^{l}y}_{1}^{1}, \dots, {{}^{l}y}_{1}^{K}, \dots, {{}^{l}y}_{n}^{1}, \dots, {{}^{l}y}_{n}^{K})$.

We assume that the variable ${t}_{\beta, \lambda }\in \mathbb{N}$ determines traveling time between nodes ${N}_{\beta }$, ${N}_{\lambda }$, where: $\left({N}_{\beta }, {N}_{\lambda }\right)\in E$, and routes of ${U}_{k}\in {}^{l}\mathcal{U}$ during sub-mission ${}^{l}S$ are represented by sequences: ${{}^{l}\pi }_{k} = \left({N}_{{k}_{1}}, \dots, {N}_{{k}_{i}}, {N}_{{k}_{i+1}}, \dots, {N}_{{k}_{\mu }}\right)$, where: ${k}_{i}\in \{1,.., n\}$, $\left({N}_{{k}_{i}}, {N}_{{k}_{i+1}}\right)\in E$. ${}^{l}\Pi$ denotes a sequence of UAVs routes executed during sub-mission ${}^{l}S$ : ${}^{l}\Pi = \left({{}^{l}\pi }_{1}, \dots, {{}^{l}\pi }_{k}, \dots, {{}^{l}\pi }_{K}\right)$ (in cases when ${U}_{k}\notin {}^{l}\mathcal{U}$ then ${{}^{l}\pi }_{k} = △$). The delivery plan of one UAV sub-mission ${}^{l}S$ is defined as a sequence: ${}^{l}S = \left({}^{l}\mathcal{U}, {}^{l}\Pi, {}^{l}Y, {}^{l}C\right)$.

It is assumed that a plan of a UAV sub-mission ${}^{l}S$ is implemented under specific weather conditions, i.e., the weather forecast is known for each sub-mission ${}^{l}S$. The forecasted weather conditions are described by the set $\mathbb{F}$ of pairs composed of direction $\theta$ and wind speed $vw$ $\left(\theta, vw\right)\in \mathcal{F}$, i.e., $\mathbb{F}$ is defined as follows: $\mathbb{F} = \left\{\left(\theta, vw\right)|\theta \in [0°, 360°), vw\in [0, \mathcal{F}\left(\theta \right)]\right\}$. Where $\mathcal{F}\left(\theta \right)$ is a function that's values determine the maximum forecasted wind speed for given direction $\theta$. The weather conditions determine the admissibility of the adopted sub-mission's plan ${}^{l}S$, i.e., they determine whether during its implementation, the batteries of UAVs will not be prematurely discharged.

A function ${\mathrm{{\rm Y}}}_{k, l}\left(\theta \right)$ determines the borderline wind speed (for a given direction $\theta$), which guarantees the successful completion of the delivery plan by the ${U}_{k}$ during sub-mission ${}^{l}S$ in the distribution network $G$ : ${\mathrm{{\rm Y}}}_{k, l}\left(\theta \right) = \mathrm{max}{\mathrm{\Gamma }}_{k, l}\left(\theta \right)$, where: ${\mathrm{\Gamma }}_{k, l}\left(\theta \right)$ – set of wind speed's values $vw$ for a given direction $\theta$, for which the batteries of ${U}_{k}$ is not discharged. The sub-mission's plan ${}^{l}S$ is assumed ^[14] to be resistant to the forecast weather conditions $\mathbb{F}$ if the boundary wind ${\mathrm{{\rm Y}}}_{k, l}\left(\theta \right)$ of all ${U}_{k}\in {}^{l}\mathcal{U}$ in any direction $\theta$ does not exceed the forecasted value $Z\left(\theta \right)$ : ${\forall }_{{U}_{k}\in {}^{l}\mathcal{U}}{{\forall }_{\theta \in [0°, 360°)}{\rm Y}}_{k, l}\left(\theta \right)\ge \mathcal{F}\left(\theta \right).$

The implementation of the designated mission $S$ may be subject to various disturbances $IS$. Among them, there are sudden changes in the weather (beyond the expected ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right)$ range), and changes in orders $Z$), order changes (increase /decrease amount of ordered deliveries ${Z}^{*}$), changes in the number of delivery points served (changing the structure of the network ${G}^{*}$). The UAV fleet, when performing the mission plan $S$, meets a disturbance $IS\left({t}^{*}\right)$- covering one of the cases: the weather ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right)$, the network ${G}^{*}$, orders ${\mathrm{Z}}^{\mathrm{*}}$, at the time ${t}^{*}.$ In such situations, it becomes necessary to answer the question: Does a re-route plan exist for mission ${S}^{*}$ that guarantees timely deliveries in a given time horizon $H$ and at acceptable battery level? To illustrate the motivation behind our approach, let us consider a distribution network from Figure 1a) covering an area of 100 km² and containing 39 delivery points (nodes ${N}_{2}, \dots, {N}_{40}$).

The goods are delivered by a fleet of UAVs, which is stationed in the base ${N}_{1}$ –technical parameters of UAVs are collected in the table from Figure 1c). The weight of individual orders is: ${\mathrm{z}}_{2}-{\mathrm{z}}_{6} = 5\;, {\mathrm{z}}_{7} = 10, {\mathrm{z}}_{18} = {\mathrm{z}}_{19} = 15, {\mathrm{z}}_{20}-{\mathrm{z}}_{26} = 5, {\mathrm{z}}_{27} = 10, {\mathrm{z}}_{28} = 15, {\mathrm{z}}_{30}-{\mathrm{z}}_{36} = 5, {\mathrm{z}}_{37} = 10, {\mathrm{z}}_{38} = {\mathrm{z}}_{39} = 15, {\mathrm{z}}_{40} = 5$. In that context, we search for the minimum fleet size that guarantees timely deliveries of the required amount of goods. It is apparent that the smallest fleet that guarantees timely delivery (within the time horizon of 2.5 h) consists of $\mathcal{U} = \left\{{U}_{1}, {U}_{2}, {U}_{3}, {U}_{4}\right\}$ UAVs, see Figure 2. The mission $S$ of Figure 2 consists of 6 sub-missions: $S = \left({}^{1}S, {}^{2}S, \dots, {}^{6}S\right)$ following conditions where the wind speed does not exceed 9 m/s.

An example illustrating the course of UAV fleet mission ${}^{1}S$ and the values of the resistance functions: ${{\rm Y}}_{\mathrm{1, 1}}\left(\theta \right)$, ${{\rm Y}}_{\mathrm{2, 1}}\left(\theta \right)$, ${{\rm Y}}_{\mathrm{3, 1}}\left(\theta \right)$, ${{\rm Y}}_{\mathrm{4, 1}}\left(\theta \right)$ is shown in Figure 2. It is evident that the UAVs' routes in the mission ${}^{1}S$ are weatherproof for the given forecasted weather (i.e., ${{\rm Y}}_{k, l}\left(\theta \right)\ge \mathcal{F}\left(\theta \right)):$

and the robustness function UAVs: ${{\rm Y}}_{\mathrm{1, 1}}\left(\theta \right)$, ${{\rm Y}}_{\mathrm{2, 1}}\left(\theta \right)$, ${{\rm Y}}_{\mathrm{3, 1}}\left(\theta \right)$, ${{\rm Y}}_{\mathrm{4, 1}}\left(\theta \right)$.

All goods should be delivered within 2.5 hours ($T = 9000\; \left[\mathrm{s}\right])$. The deliveries take place in different forecasted weather conditions (set $\mathbb{F}$), which are illustrated in Figure 1b). According to the forecast, the wind speed does not exceed $vw = 9\frac{m}{s}$.

We consider a situation in which the weather conditions of the mission carried out rapidly changed at the time ${t}^{*} = 3000$ [s], i.e., during the execution of the sub-mission ${}^{2}S$, the wind speed increased to $vw = 11m/s$ for direction $\theta = 210°-230°$.

Such a change means that this mission cannot be continued due to too much energy consumption (the mission's resistance function ${{\rm Y}}_{\mathrm{3, 2}}\left(\theta \right)$ values are below the level corresponding to speed $11\frac{m}{s}$ – see Figure 3). Figure 3 shows the location of the UAVs at time ${t}^{*} = 3000 \left[s\right]$, i.e., upon receipt of information about a change in weather, and marked the place where the battery ${U}_{3}$ will be discharged in the event of continuation of deliveries in accordance with the current plan ${}^{2}S$.

In this situation, it is necessary to correct the route of the sub-mission ${}^{2}S$ being carried out, which forces the search for an answer to the following question:

Given a UAV fleet $\mathcal{U} = \left\{{U}_{1}, {U}_{2}, {U}_{3}, {U}_{4}\right\}$ providing deliveries to the delivery points allocated in the network $G$ from Fig. 1, the UAV fleet realizes the delivery mission plan ${S}^{*}$ from fig. 2 b). At the time ${t}^{*} = 3000$ a rapid weather ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right)$ change, i.e., disturbance $IS$, occurs and results in $vw = 11m/s$; $\theta = 210°-230°$. Does a reroute plan exist for mission ${S}^{*}$ that guarantees the timely delivery of the ordered supplies in a given time horizon $H = 9000\; \left[s\right]$ and at an acceptable battery level?

4.   Approach to UAV fleet online routing

Planning drone missions is founded on setting the so-called proactive plans specifying the UAV flight routes (see Figure 4) to guarantee timely delivery to each recipient while taking into account the constraints of the assumed distribution network $G$ and forecasted weather $\mathcal{F}\left(\theta \right)$ condition changes. This means that routes and schedules determining the flights of individual UAVs (obtained in the course of proactive fleet mission design) take into account the amount and magnitude of the anticipated disturbances.

In other words, proactive planning is carried out before the mission is carried out and boils down to setting routes and schedules, guaranteeing timely deliveries in the given distribution network even in a situation of changing weather (within the set range specified by weather forecast). The implementation of such plans in a real environment may be disrupted as a result of a sudden unforeseen change in weather conditions (exceeding the range assumed at the proactive planning stage), which may lead to premature depletion of the UAV's battery. In such situations, attempts are made to replace proactive plans by implementing reactive planning principles. Reactive planning (Figure 5), implies changing a proactive plan already being implemented, e.g., due to the occurrence of a disturbance such as rapid weather and/or orders change. Reactive plans must therefore take into account the extent to which the planned deliveries were carried out and adjust the further actions of the UAVs to the prevailing weather conditions. The purpose of reactive planning is therefore an attempt to change mission plans to those that guarantee timely deliveries under new weather conditions. Due to the need to make decisions in the online mode, reactive planning methods must have a short time of determination of the solution.

The UAV mission planning process (Figure 6) therefore includes two stages. The first one involves proactive planning, which results in a UAVs mission plan that guarantees timely deliveries in the forecasted weather conditions. The second stage is responsible for reactive planning undertaken in situations where during the implementation of the proactive plan there is a disruption exceeding previous predictions.

In that context, the proposed reaction (replanning) to the occurrence of a disturbance $IS\left({t}^{\mathrm{*}}\right)$ can be reduced to dynamic re-routing and rescheduling of previously adopted routes ${}^{l}\Pi$, schedules ${}^{l}Y$, and delivered goods ${}^{l}C$ is stated in the basic proactive plan for the mission by the UAV fleet. It is a feasible adjustment of the assumed ${}^{l}\Pi$, ${}^{l}Y$, and ${}^{l}C$ values to the changes in forecasted weather ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right)$, as well as corrections introduced to the network ${G}^{*}$ or orders ${Z}^{*}$.

To formally define the concept of disturbance $IS\left({t}^{*}\right),$ we will introduce the concept of the state of mission implementation $S$. The state of the mission $S$ at the time $t$ is defined as follows: $IS\left(t\right) = \left(M\left(t\right), {\mathcal{F}}^{\mathcal{*}}\left(\theta, t\right), {}^{\mathrm{*}}G\left(t\right), {Z}^{\mathrm{*}}\left(t\right)\right)$ where: $M\left(t\right)$ is an allocation of UAVs to nodes at the time $t$ : $M\left(t\right) = \left({N}_{{a}_{1}}, \dots, {N}_{{a}_{k}}, \dots, {N}_{{a}_{K}}\right)$, where: ${a}_{k}\in \{1,.., n\}$ determines the (delivery points) node ${N}_{{a}_{k}}$ occupied by ${U}_{k}$ (or the node the ${U}_{k}$ is headed to). ${\mathcal{F}}^{\mathcal{*}}(\theta, t)$ is the weather condition forecast at the time $t$. ${}^{*}G\left(t\right)$ is the graph model of the distribution network structure at time t (number and location of delivery points). ${Z}^{*}\left(t\right)$ is the sequence of goods requested at the time $t$. The state $IS\left({t}^{*}\right)$ following condition $\left[{\mathcal{F}}^{\mathcal{*}}\left(\theta, {t}^{*}\right)\ne {\mathcal{F}}^{\mathcal{*}}\left(\theta \right)\right]\vee \left[{}^{*}G\left({t}^{*}\right)\ne G\right]\vee \left[{Z}^{*}\left({t}^{*}\right)\ne Z\right]$ is called the disturbance occurring at the time ${t}^{*}$.

Occurrence of $IS\left({t}^{*}\right)$ disturbance should be assessed in terms of its impact on the further course of the mission of S (that is, whether the value of the resistance function ${{\rm Y}}_{k, l}\left(\theta \right)$ is greater than $\mathcal{F}\left(\theta \right)$). If the implementation of the mission is at risk (${{\rm Y}}_{k, l}\left(\theta \right)\ngeqq \mathcal{F}\left(\theta \right)$), an attempt should be made to reschedule it. The following condition action (if-then) rules are used for this purpose:

1)   If the adopted mission plan $S$ is not resistant to disturbance $IS\left({t}^{*}\right)$ $({\exists }_{\mathrm{k}\in \left\{1,.., K\right\}, \mathrm{l}\in \left\{1,.., L\right\}}{\mathrm{{\rm Y}}}_{k, l}\left(\mathrm{\theta }\right)\ngeqq \mathcal{F}\left(\mathrm{\theta }\right)),$ then it should be checked whether it is possible to adapt (re-plan) it to adjust to new conditions, i.e., decide whether all UAVs in the air should continue their current missions or make the appropriate corrections.

2)   If there are UAVs (in the set $\mathcal{U}\mathrm{R}$) that cannot continue to fly due to disturbance $IS\left({t}^{*}\right),$ then they should be returned to the base and allow, if possible airborne UAVs ($\mathrm{t}\mathrm{h}\mathrm{e}\; \; \; \; \mathrm{s}\mathrm{e}\mathrm{t}\mathcal{U}\backslash \mathcal{U}\mathrm{R}$) to take over their tasks.

3)   If the tasks of the UAVs returning to the base (the set $\mathcal{U}\mathrm{R})$ cannot be taken over by UAVs still performing their missions, then it should be checked whether the reserve UAVs available at the base (the set $\mathcal{U}\mathrm{B}$) can take over their responsibilities. This means the UAVs in the air continue their existing missions, while the reserve UAVs take over the liabilities of the UAVs returned to the base.

4)   If the reserve UAVs (the set $\mathcal{U}\mathrm{B}$) are unable to take over the responsibilities of those returned to the base (the set $\mathcal{U}\mathrm{R}$), then their activity should be suspended until the disturbance is resolved.

The above rules have been used in the reactive mission planning method $S$. The idea behind this method is as follows. During the implementation of the mission $S$, there is a continuous monitoring of the state of the $IS\left(t\right)$ (for $t\in \{0\dots H\}$). If at the state $IS\left(t\right)$ the following condition holds $\left[{\mathcal{F}}^{\mathcal{*}}\left(\theta, {t}^{\mathrm{*}}\right)\ne {\mathcal{F}}^{\mathcal{*}}\left(\theta \right)\right]\vee \left[{}^{\mathrm{*}}G\left({t}^{\mathrm{*}}\right)\ne G\right]\vee \left[{Z}^{\mathrm{*}}\left({t}^{\mathrm{*}}\right)\ne Z\right]$ (i.e., there is a change in the weather forecast or in the structure of the distribution network as well as in the size of the requests, etc.) and the mission $S$ is under threat (i.e., at least one of the UAVs will not return to base due to low battery), then an attempt is made to replan it. The purpose of replanning is to select a mission ${}^{\mathrm{*}}S$ adapted to the new conditions caused by the disturbance $IS\left(t\right)$. In practice, it comes down to solving the relevant constraints optimization problem (COP) denoted as $CS\left({}_{}{}^{\circ}\mathcal{U}, S, IS\left(t\right)\right)$, where: ${}_{}{}^{\circ}\mathcal{U}$—defines the fleet designated by condition action rules 1–4. The relevant constraints distribution process follows the sequence where, firstly, an attempt is made to designate the mission ${}^{*}S$ for the fleet ${}_{}{}^{①}\mathcal{U} = \mathcal{U}$ (according to rule 1). In the event of failure, an attempt is made to designate it for the fleet ${}_{}{}^{②}\mathcal{U} = \mathcal{U}\backslash \mathcal{U}R$ (according to rule 2) and then for the fleet ${}_{}{}^{③}\mathcal{U} = (\mathcal{U}\backslash \mathcal{U}R)\cup \mathcal{U}B$ (according to rule 3). If an admissible solution ${}^{\mathrm{*}}S$ sitll does not exist, then the currently used mission plan should be modified (due to the $reduce$ function) in such a way that it removes the UAVs sub-missions which are not resistant to disturbance $IS\left(t\right)$ (according to rule 4).

It should be noted that the designation of mission ${}^{\mathrm{*}}S$ is associated with the designation of routings ${}^{l}\Pi$, schedules ${}^{l}Y,$ and delivery sequences ${}^{l}C$ throughout the remaining time horizon $\{t\dots H\}$. Due to the disturbance $IS\left({t}^{*}\right)$ occurrence, the proposed reactive planning algorithm (implemented in the IBM ILOG environment) generates the end-to-end paths that modify previously planned routes, restoring the ability to implement the designated mission delivery plan.

The mathematical formulation of the COP $CS\left({}_{}{}^{\circ}\mathcal{U}, S, IS\left(t\right)\right)$ aimed at contingency planning of UAV mission design employs the following parameters, variables, sets, and constraints.

Parameters:

${}^{l}G$ : graph of a distribution network for sub-mission ${}^{l}S$;

${}^{l}\mathcal{U}$ : the subset of UAVs ${}^{l}\mathcal{U}\subseteq \mathcal{U}$ carrying out the sub-mission ${}^{l}S$;

${z}_{\lambda }$ : demand at node ${N}_{\lambda }$, ${z}_{1} = 0$;

$K$ : the size of the fleet of UAVs;

$n{p}_{\lambda }$ : consumer priority at the point ${N}_{\lambda }$, $n{p}_{1} = 0$;

${d}_{\beta, \lambda }$ : distance between ${N}_{\beta }$, ${N}_{\lambda }$;

$A$ : the front-facing area of a UAV;

${t}_{\beta, \lambda }$ : travel time between ${N}_{\beta }$, ${N}_{\lambda }$;

${C}_{D}$ : the aerodynamic drag coefficient;

$w$ : time spent on take-off and landing of a UAV;

$ts$ : the time interval at which UAVs can take off from the base;

$Q$ : maximum loading capacity;

$ep$ : the empty weight of a UAV;

$IS\left(t\right)$ : state of UAV mission;

$D$ : an air density;

$H$ : time horizon;

${\mathrm{{\rm Y}}}_{k, l}\left(\theta \right)$ : weather resistance function;

$g$ : the gravitational acceleration;

$\mathcal{F}\left(\theta \right)$ : forecasted wind speed;

$b$ : the width of an UAV;

$CAP$ : the energy capacity of an UAV;

${va}_{\beta, \lambda }^{}$ : air speed between nodes ${N}_{\beta }$, ${N}_{\lambda }$;

${vg}_{\beta, \lambda }^{}$ : ground speed between ${N}_{\beta }$, ${N}_{\lambda }$;

${\phi }_{\beta, \lambda }$ : heading angle of vector ${va}_{\beta, \lambda }^{}$;

${\vartheta }_{\beta, \lambda }$ : the course angle of vector ${vg}_{\beta, \lambda }^{}$.

Decision variables:

$\overline{{{}^{l}x}_{\beta, \lambda }^{k}}$ : the binary variable used to indicate if ${U}_{k}$ travels between nodes ${N}_{\beta }$, ${N}_{\lambda }$, after the disturbance $IS\left({t}^{*}\right)$ occurrence (during sub-mission ${}^{l}S$);

$\overline{{{}^{l}y}_{\lambda }^{k}}$ : the time at which ${U}_{k}$ arrives at node ${N}_{\lambda }$, after the disturbance $IS\left({t}^{*}\right)$ occurrence (during sub-mission ${}^{l}S$);

$\overline{{{}^{l}c}_{\lambda }^{k}}$ : the weight of freight delivered to node ${N}_{\lambda }$ by ${U}_{k}$, after the disturbance $IS\left({t}^{*}\right)$ occurrence (during sub-mission ${}^{l}S$);

$\overline{{{}^{l}f}_{\beta, \lambda }^{k}}$ : the weight of freight carried between nodes ${N}_{\beta }$, ${N}_{\lambda }$ by ${U}_{k}$, after the disturbance $IS\left({t}^{*}\right)$ occurrence (during sub-mission ${}^{l}S$);

$\overline{{{}^{l}P}_{\beta, \lambda }^{k}}$ : the energy per unit of time consumed by ${U}_{k}$ during the flight between nodes ${N}_{\beta }$, ${N}_{\lambda }$ (after the disturbance $IS\left({t}^{*}\right)$ occurrence);

$\overline{{}^{l}ba{t}^{k}}$ : the total energy consumed by ${U}_{k}$, after the disturbance $IS\left({t}^{*}\right)$ occurrence (during sub-mission ${}^{l}S$);

$\overline{{{}^{l}s}^{k}}$ : the take-off time of ${U}_{k}$, after the disturbance $IS\left({t}^{*}\right)$ occurrence;

$\overline{{}^{l}{cp}_{\lambda }}$ : the total weight of freight delivered to node ${N}_{\lambda }$, after the disturbance $IS\left({t}^{*}\right)$ occurrence (during sub-mission ${}^{l}S$);

$\overline{{}^{l}{\pi }_{k}^{}}$ : the route of ${U}_{k}$ after the disturbance $IS\left({t}^{*}\right)$ occurrence (during sub-mission ${}^{l}S$), $\overline{{{}^{l}\pi }_{k}} = ({N}_{{k}_{1}}, \dots, {N}_{{k}_{i}}, {N}_{{k}_{i+1}}, \dots, {N}_{{k}_{\mu }})$, ${k}_{i}\in \{1,.., n\}$, $\left({N}_{{k}_{i}}, {N}_{{k}_{i+1}}\right)\in E$.

Sets:

$\overline{{}^{l}Y}$ : is a sequence of moments $\overline{{{}^{l}y}_{\lambda }^{k}}$, schedule of the fleet ${}^{l}\mathcal{U}$ after the disturbance $IS\left({t}^{*}\right)$ occurrence;

$\overline{{}^{l}C}$ : is a sequence of weights of delivered goods $\overline{{{}^{l}c}_{\lambda }^{k}}$;

$\overline{{}^{l}\Pi }$ : the set of UAV routes $\overline{{}^{l}{\pi }_{k}^{}}$;

${\overline{{}^{l}S}}_{}$ : the plan of sub-mission after the disturbance $IS\left({t}^{*}\right)$ occurrence: $\overline{{}^{l}S} = \left({}^{l}\mathcal{U}, \overline{{}^{l}\Pi }, \overline{{}^{l}Y}, \overline{{}^{l}C}\right)$;

${}^{*}S$ : the re-route plan of the mission ${}^{*}S = \left({\overline{{}^{1}S}}_{}, \dots, {\overline{{}^{l}S}}_{}, \dots, {\overline{{}^{L}S}}_{}\right).$

Constraints limiting: routes, delivery of freight, and energy consumption:

1)   Routes. Relationships between the variables describing UAV take-off times/mission start times and task order:

Delivery of freight. Relationships between variables describing already delivered and requested amount of freight:

2)   Energy consumption. To ensure waterproofness of the ${}^{l}S$ sub-mission (i.e., its robustness to weather condition changes $Z\left(\theta \right)$), the amount of energy required to complete the task carried out by an UAV must not exceed the capacity of its battery:

where, ${{}^{l}va}_{i, j}^{}(\theta, vw)$ and ${t}_{i, j}$ depend on the assumed goods delivery strategy.

If the ground speed ${vg}_{i, j}^{}$ is constant, then an air speed ${{}^{l}va}_{i, j}^{}$ is calculated from:

3)   The objective function. A mission ${\overline{{}^{l}S}}_{}$ plan to maximize customer satisfaction expressed by the following function is sought:

The adoption of such a function means that customers with the lowest level of expectations to be met are served first.

Constraints (1)–(13) describe the relationship between routes (represented by the variables: $\overline{{{}^{l}x}_{i, j}^{k}}$) and the delivery schedule (variables $\overline{{{}^{l}y}_{j}^{k}}$ and $\overline{{{}^{l}s}^{k}}$). Constraints (5)–(7) provide, that the plan before disturbance should be the same like proactive plan (represented by ${{}^{l}s}^{k}$, ${{}^{l}x}_{i, j}^{k}$, ${{}^{l}y}_{j}^{k}$). Other constrains provide, that it is not possible for multiple UAVs to take off from the base at the same time (8), a delivery points cannot by simultaneously served by several UAVs (9), deliveries are made in accordance with the adopted route (10), (11) and on time (12) and it guarantees closed loops of the routes (13).

Constraints (14)–(24), link UAV routes ($\overline{{{}^{l}x}_{i, j}^{k}}$) to the amounts of delivered goods (variables $\overline{{{}^{l}c}_{j}^{k}}$). They also ensure that the UAVs are not overloaded (16), (17), and that correct amounts are delivered (19). Constraints (20)–(24) determine the weight (${f}_{i, j}^{k}$) of the goods at each section of the taken route.

Constraints (25)–(31) describe the values of the determined weather resistance functions ${{\rm Y}}_{k, l}\left(\theta \right)$ for the fleet ${}^{l}\mathcal{U}$ and ensure that these values exceed the value of the function $\mathcal{F}\left(\theta \right)$ (forecasted wind speed). The value of the ${{\rm Y}}_{k, l}\left(\theta \right)$ function depend on the amount of energy consumed by a UAV in flight which in turn depends non-linearly on the speed value ${{}^{l}va}_{i, j}^{}(\theta, vw)$. This means that some of the constraints, e.g. (29), (30), in the adopted model, have a non-linear character, thus implying the necessity to use the capabilities of declarative environments (in particular constraints programming).

Since the re-planning of the mission delivery plan $S$ is the result of the disturbance $IS\left({t}^{\mathrm{*}}\right)$, the new set of sub-missions ${\overline{{}^{1}S}}_{}, \dots, {\overline{{}^{l}S}}_{}, \dots, {\overline{{}^{L}S}}_{}$ guaranteeing timely delivery are determined by solving the following COP (33):

where:

$\widehat{\mathcal{V}} = \left\{\overline{{}^{l}\Pi }, \overline{{}^{l}Y}, \overline{{}^{l}C}|l = 1\dots L\right\}$—the set of decision variables: $\overline{{}^{l}\Pi }$—the set of routes determining the schedule $\overline{{}^{l}Y}$; $\overline{{}^{l}Y}$—schedule of the fleet ${}_{}{}^{\circ}\mathcal{U}$ guarantees timely service of delivery points in the case of disturbance $IS\left({t}^{*}\right)$; and $\overline{{}^{l}C}$—sequence of weights of delivered goods by the fleet ${}_{}{}^{\circ}\mathcal{U}$;

$\mathcal{D}$—the finite set of decision variable domains: $\overline{{{}^{l}x}_{i, j}^{k}}\in \left\{\mathrm{0, 1}\right\}$, $\overline{{{}^{l}y}_{\lambda }^{k}}\in \mathbb{N}$, $\overline{{{}^{l}c}_{i}^{k}}\in \mathbb{N}$;

$\widehat{\mathcal{C}}$—the set of constraints that takes into account the set of routes $\overline{{}^{l}\Pi }$, schedules $\overline{{}^{l}Y}$, and the disturbance $IS\left({t}^{*}\right)$, while determining the relationships linking the operations executed by UAVs (1)–(31);

${f}_{o}$—the objective function defined by (32).

To solve COP (33), the values of the decision variables from the adopted set of domains for which the given constraints are satisfied and guarantee the maximum value of objective function must be determined.

5.   Constraint relaxation

The previous section demonstrated that solving the COP (33) problem enables the designation of an ${}^{\mathrm{*}}S$ mission resistant to $IS\left({t}^{\mathrm{*}}\right)$ disturbances caused by changing weather conditions specified by $\mathcal{F}\left(\theta \right).$ Solving this problem, however, is very time-consuming, which results from the necessity to verify the inequality (25) of the wind direction change ${{\rm Y}}_{k, l}\left(\theta \right)\ge \mathcal{F}\left(\theta \right)$ carried out for each increment $\theta \in [{0}^{\circ }, {360}^{\circ })$. To replace condition (25) with an equivalent that is less time consuming computationally, it was assumed that the continuous (smooth) ${{\rm Y}}_{k, l}\left(\theta \right)$ function will be approximated by a discrete function represented by a finite set ${\mathbb{Y}}_{k, l} = \left\{{{\rm Y}}_{k, l}\left({\theta }_{i}\right)|i = 1\dots lq; {\theta }_{i} < {\theta }_{i+1}\right\},$ where $lq$ is an arbitrarily taken number of samples. This set contains the vertices of a polygon ${\mathbb{Y}}_{k, l}^{*}\left(\theta \right)$ depicted in Figure 7.

For such assumptions, the following property holds:

Property 1: For any wind direction $\theta \in [{0}^{\circ }, {360}^{\circ })$ function of ${{\rm Y}}_{k, l}\left(\theta \right)$ mission, ${}^{*}S$ takes values no less than its discrete form ${\mathbb{Y}}_{k, l}^{*}\left(\theta \right)$ : $\forall \theta \in [{0}^{\circ }, {360}^{\circ })$, ${{\rm Y}}_{k, l}\left(\theta \right)\ge$ ${\mathbb{Y}}_{k, l}^{*}\left(\theta \right)$ .

This means that the function ${{\rm Y}}_{k, l}\left(\theta \right)$ can be replaced by the function ${\mathbb{Y}}_{k, l}^{*}\left(\theta \right)$ and, in fact, the set of vertices in ${\mathbb{Y}}_{k, l}$. In practice, it means reducing the quantity of condition (25) checks to the number of vertices of the adopted polygon. Consequently, constraint (25) can be replaced by the following one:

6.   Evolutionary algorithm

The introduction of constraints relaxation (34) of the considered problem (33) reduces the computational effort but not adequately to solve problems of the scale of those occurring in practice. This means that in online mode, problems of the type (33) involve networks not exceeding 100 nodes. To speed up the calculations, a genetic algorithm dedicated to this problem has been developed.

For the purposes of the developed algorithm, it is assumed that sub-mission $\overline{{}^{l}S}$ is represented by chromosons ${}^{l}\overline{{CH}_{}}$ describing deliveries made by the UAVs of fleet ${}^{l}\mathcal{U}$ :

In other words, ${}^{l}\overline{{CH}_{}}$ consists of the routes $\overline{{}^{l}{\pi }_{k}^{}}$ ($k = 1, \dots, K$); schedules $\overline{{}^{l}{Y}_{k}^{}}$ ($k = 1, \dots, K$); and volumes of deliveries $\overline{{}^{l}{C}_{k}^{}}$ ($k = 1, \dots, K$) describing the sub-mission $\overline{{}^{l}S}$.

The framework of the proposed genetic algorithm implementing such chromosomes is presented in Fig. 8. According to the algorithm flowchart, its input data are: $S = \left({}^{1}S, \dots, {}^{l}S, \dots, {}^{L}S\right)—$ the proactive mission plan, $IS\left({t}^{*}\right)$—the disturbance occurring at the moment ${t}^{*}$, ${}^{l}S$—the submission affected by the disturbance, ${}^{l}\mathcal{U}$—the fleet of UAVs available in reactive mode (while following rules 1–4).

The next stages of the algorithm include:

1)   Determination of the j-th population ${}^{l}{PO}_{}^{j}$ represented by the set of chromosomes describing different ways to accomplish a sub-mission ${}^{l}S:$ ${}^{l}{PO}_{}^{j} = \left\{{{}^{l}\overline{{CH}_{}}}_{i}^{j}|i = 1\dots LP\right\}$, where ${{}^{l}\overline{{CH}_{}}}_{i}^{j}$—is the i-th chromosome from j-th population following submission ${}^{l}S$ (defined according to (35)); $LP$—population size. Individuals of the initial population ($j = 0$) are selected at random from among the so-called "alive individuals"—solutions that meet the constraints (1)–(31). Individuals of population ${}^{l}{PO}_{}^{j}$ for $j > 0$ are selected from the set ${}^{l}{PO}_{}^{j-1}\cup {}^{l}{DE}_{}^{j-1}\cup {}^{l}{MU}_{}^{j-1}$ ($LP$ best individuals population).

2)   Population growth rate. For each chromosome ${{}^{l}\overline{{CH}_{}}}_{i}^{j}\in {}^{l}{PO}_{}^{j}$, the value of the objective function is determined ${f}_{o}\left({{}^{l}\overline{{CH}_{}}}_{i}^{j}\right)$.

3)   Selection. From among the individuals of the population ${}^{l}{PO}_{}^{j}$, a subset ${}^{l}S{PO}_{}^{j}\subseteq {}^{l}{PO}_{}^{j}$ is drawn (it is assumed that that the parameter $ps$ determines the probability of the event that ${{}^{l}\overline{{CH}_{}}}_{i}^{j}$ belongs to the set ${}^{l}S{PO}_{}^{j}$ : $ps = P\left({{}^{l}\overline{{CH}_{}}}_{i}^{j}\in {}^{l}S{PO}_{}^{j}\right)$).

4)   Crossover. Subsequent pairs of the set ${}^{l}S{PO}_{}^{j}$ are crossed with each other. The result of the operation of crossing two chromosomes ${{}^{l}\overline{{CH}_{}}}_{u}^{j}$ and ${{}^{l}\overline{{CH}_{}}}_{v}^{j}$ is a pair of chromosomes $\dot{{{}^{l}\overline{{CH}_{}}}_{u}^{j}}$ and $\dot{{{}^{l}\overline{{CH}_{}}}_{v}^{j}}$, in which one (randomly selected) k-th element of the route $\overline{{{}^{l}\pi }_{u, k}^{j}}$ is replaced by another selected at random k-th element of the route $\overline{{{}^{l}\pi }_{v, k}^{j}}$. The idea of crossing operation is illustrated by Figure 9. It should be noted that along with the replacement of route elements, the corresponding elements of the schedule $\overline{{{}^{l}Y}_{u, k}^{j}}$ and delivery ${{}^{l}C}_{u, k}^{j}$ sequences are replaced. The result of the crossing process is a set of descendants ${}^{l}{DE}_{}^{j}$ of the population ${}^{l}{PO}_{}^{j}$ (only alive individuals are included in the set ${}^{l}{DE}_{}^{j}).$

5)   Mutation. From among the elements of the set ${}^{l}{PO}_{}^{j}\cup {}^{l}{DE}_{}^{j}$, a set of individuals undergoing mutation is selected at random (with a probability of $pm$). Operation of the chromosome ${{}^{l}\overline{{CH}_{}}}_{u}^{j}$ mutation consists of a random change of one (randomly selected) element of the $k$ route $\overline{{{}^{l}\pi }_{u, k}^{j}}$. As a result of this operation, a set of mutated individuals ${}^{l}{MU}_{}^{j}$ of population ${}^{l}{PO}_{}^{j}$ is formed.

According to the presented algorithm, operations used for Selection, Crossover, and Mutation of individuals of the population ${}^{l}{PO}_{}^{j}$ are repeated until the stopping condition (that can be seen as a stopping rule or simply stopping criteria) is satisfied. In the considered version of the algorithm, the following types of the stopping condition are distinguished:

●  Condition 1 (C1): value of the objective function of the best individual of the population $\underset{i = 1\dots \mathit{LP}}{\mathrm{max}}\left\{{f}_{o}\left({{}^{l}\overline{{CH}_{}}}_{i}^{j}\right)\right\}$ has reached the expected value $RV$ : $\underset{i = 1\dots \mathit{LP}}{\mathrm{max}}\left\{{f}_{o}\left({{}^{l}\overline{{CH}_{}}}_{i}^{j}\right)\right\} = RV$;

●  Condition 2 (C2): value of the objective function of the best individual of the population $\underset{i = 1\dots \mathit{LP}}{\mathrm{max}}\left\{{f}_{o}\left({{}^{l}\overline{{CH}_{}}}_{i}^{j}\right)\right\}$ has not changed in $q$ generations: $\underset{i = 1\dots \mathit{LP}}{\mathrm{max}}\left\{{f}_{o}\left({{}^{l}\overline{{CH}_{}}}_{i}^{j}\right)\right\} = \underset{i = 1\dots \mathit{LP}}{\mathrm{max}}\left\{{f}_{o}\left({{}^{l}\overline{{CH}_{}}}_{i}^{j-1}\right)\right\} = \dots = \underset{i = 1\dots \mathit{LP}}{\mathrm{max}}\left\{{f}_{o}\left({{}^{l}\overline{{CH}_{}}}_{i}^{j-q}\right)\right\}$;

● Condition 3 (C3): the number of generations has reached the limit value $le$ : $j = le$.

Satisfaction of the stopping condition allows you to determine the next submissions $\overline{{}^{l}S}$ (where $\overline{{}^{l}S} = best\left({{}^{l}\overline{{CH}_{}}}_{i}^{j}\right)$) constituting the reactive UAV fleet mission plan $\overline{S} = \left({}^{1}S, \dots, {}_{}{}^{l-1}S, \overline{{}^{l}S}, \dots, \overline{{}^{L}S}\right)$.

7.   Computational experiments

We consider the network from Figure 1, in which the four UAVs $\mathcal{U} = \left\{{U}_{1}, {U}_{2}, {U}_{3}, {U}_{4}\right\}$ service delivery points ${N}_{2}–{N}_{40}$. The structure of the implementation of subsequent sub-missions ${}^{1}S, {}^{2}S, \dots, {}^{6}S$ of the adopted mission plan $S$ is presented in Figure 2.

Let us consider a situation related to the appearance of a disturbance IS(3000), specified in Section 2, where at the time ${t}^{*} = 3000$ [s], during the execution of the sub-mission ${}^{2}S$, the wind speed increased to $vw = \; 11m/s$ with the same wind intensity and direction $\theta = 210°-230°$. With such a change in weather conditions, the implementation of the adopted plan turns out to be impossible. It becomes necessary to re-plan the implemented mission, including the introduction of the new sub-missions $\overline{{}^{2}S}, \overline{{}_{}{}^{3}S}$, $\overline{{}_{}{}^{4}S}$, $\overline{{}_{}{}^{5}S}, \overline{{}^{6}S}$ to correct its course. For this purpose, the considered problem has been modeled in $COP$ (33) formalism. The Intel Core i7-M4800MQ 2.7 GHz, 32 GB RAM has been used to carry out the necessary calculations, employing:

1)   The constraint programming environment IBM ILOG,

2)   The genetic algorithm (shown on the Figure 8) implemented in Matlab environment.

Ad. 1. Implementation of considered problem in the constraint programming environment IBM ILOG has shown that the solution time for problems of size considered does not exceed 25 s. Figure 10 demonstrates the mission ${}^{*}S$ schedule being adapted to the weather conditions caused by the disturbance IS(3000). Rule 2, i.e., If there are UAVs that (within the set $\mathcal{U}R$ ) cannot continue to fly due to disturbance $IS\left({t}^{*}\right),$ then they should be returned to the base to allow, if possible, airborne UAVs (the set $\mathcal{U}\backslash \mathcal{U}R$ ) to take over their tasks, was used in the selection of mission ${}^{*}S$.

Consequently, due to disturbance IS(3000) increasing the risk of premature battery depletion, a decision to turn ${U}_{3}$ back to the base was made (see sub-mission $\overline{{}^{2}S}$). At the same time ${U}_{4}$ continued its mission unchanged. The undertaken decision forced the necessity to reschedule subsequent sub-missions $\overline{{}_{}{}^{3}S}$, $\overline{{}_{}{}^{4}S}$, $\overline{{}_{}{}^{5}S}, \overline{{}^{6}S}$, creating a new alternative mission plan ${}^{*}S$. The mission implementing such a modified plan will end at 6900 s.

Ad. 2. Solving the problem under consideration using the genetic algorithm from Figure 8 took 2 s. Figure 11 shows the schedule of obtained mission ${}^{*}S$. To determine a proactive plan, as before, rule 2 was used. According to this rule ${U}_{3}$ is to turn back to the base. The further course of the mission is carried out in 8 sub-missions $\overline{{}^{2}S}, \overline{{}_{}{}^{3}S}$, $\overline{{}_{}{}^{4}S}$, $\overline{{}_{}{}^{5}S}, \overline{{}^{6}S}, \overline{{}_{}{}^{7}S}, \overline{{}_{}{}^{8}S}, \overline{{}_{}{}^{9}S}$. It is apparent that the obtained solution is associated with a longer mission time (8500s).

This example demonstrates that the constraint programming environment allows us to obtain solutions of better quality (mission completion time is reduced by 19%) but at the expense of increased calculation time (more than 10 times longer solution determination time).

To assess the scalability of the proposed approach in terms of the possibility of its use in an online mode (i.e., to solve the problem in < 600 s) in decision support systems, a series of quantitative experiments have been carried out. Tables A1–A3 (Appendix A) contain the results of the experiments that were conducted for the three functions of forecasted weather ${\mathcal{F}}^{}\left(\theta \right) = 9, 10, 11m/s$. The experiments were carried out for the network of $n = \mathrm{40, 60}, \dots, 220$ randomly designated delivery points (on an area of 10 km × 10 km) and a fleet consisting of $K = \mathrm{2, 3}, 4$ UAVs within the technical parameters, as shown in the table from Figure 1c). For each of the considered variants of the network, a proactive mission plan has been set out to guarantee the deliveries in the time horizon $H = \mathrm{10,000}$ s. It was assumed that at the moment ${t}^{*} =$ 2000 s, there is a change in the weather forecast (disturbance IS (2000)) that lasts until the end of the considered time horizon. The change in weather involves increasing the expected wind speed by 2m/s and equals accordingly ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right)$ = 11, 12, 13 m/s. The UAVs route planning is aimed at the reactive performance of delivery missions in a given time horizon, under expected weather conditions ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right)$.

In the conducted experiments, two approaches were compared:

1)   The constraint programming environment IBM ILOG CPLEX implemented the following stopping condition:

C0. The optimal value of the objective function has been reached or the time $TC = 600s$ allotted for calculations has elapsed.

2)   The Genetic Algorithm (GA) shown in Figure 8 specified by: $LP = 1000$ (population size); $ps = 0.25$ (probability of selection); $pm = 0.01$ (probability of mutation), and implemented in the Matlab environment using three distinct stopping conditions:

C1. value of the objective function of the best individual of the population has reached the value of solution obtained from the IBM ILOG CPLEX environment (see condition C0);

C2. value of the objective function of the best individual of the population has not changed through $q = 10$ generations;

C3. the number of generations has reached the limit value $le$ : $j = 200$.

The results (i.e., the times with relaxation $TC$ determining the reactive mission plan design) are presented in Figures 12–15 (and Tables A1–A3). Figure 12 compares the time expenditure incurred in determining the reactive plan in the ILOG environment and the developed GA (implementing the stopping condition C1). In the GA variant under consideration, it is assumed that the calculation is stopped when the value of the objective function specified by the solution obtained from the ILOG environment is reached. Calculations were carried out for weather conditions $\mathcal{F}\left(\theta \right) = 9\frac{m}{s}, \; \forall \mathrm{\Theta }\in [0°, 360°)$. Disruption (for which a reactive plan is sought) consists of changing the wind speed to ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right) = 11$ $\frac{m}{s}$.

It is evident that determining the optimal solution (maximizing objective function (32)) in the ILOG environment (dashed lines in Figure 12) requires significant computational effort. The scale of problems (number of network nodes $n$) for which it is possible to determine a solution on-line ($< 600$ s) decreases with the size of the UAV fleet and amounts to: $n = 220$ for $k = 2$; $n = 140$ for $k = 3$; and $n = 100$ for $k = 4$. Determining reactive plans for the set values of the objective function (obtained from the ILOG environment) using the implemented GA requires, in turn, a much lower computation time. As such, the considered scale of the network for all solutions were received in a time below 300 s—see ribbon-like lines presented means ($mT$) and standard deviation ($mT\pm \sigma$) of time computations. Such advantages of the GA results from the assumption that the solution (mission plan) was sought at the set (optimal) value of the objective function. Subsequent experiments were carried out on the assumption that the value of the objective function is unknown (stopping conditions C2 and C3 are used).

Figure 13 shows the results of experiments in which the stopping conditions C2 (Figure 13a) and C3 (Figure 13b) were implemented in the GA. In the first case, as before, the solution designation times for the GA do not exceed 400 s for the network size $n\le 220$.

Shorter computation times, however, came at the cost of lower "quality" solutions. This is illustrated in Figure 13, which shows the changes in the value of the objective function for the solutions obtained. These values are expressed in % of optimal solutions obtained from the ILOG environment (e.g., a value of 70% means that the solution obtained from the GA has a value of the objective function equal to 0.7 of the result obtained from the ILOG environment). Mission plans determined using the GA are characterized by a lower value of the objective function than the solutions obtained in the ILOG environment, representing respectively: $\approx 80\mathrm{\%}$ (for $k = 2$); $\approx 70\mathrm{\%}$ (for $k = 3$); and $\approx 60\mathrm{\%}$ (for $k = 4$) values obtained when using the ILOG environment.

In the second case (condition C3—see Figure 13b), the mission plans determined using the GA are characterized by a higher value of the objective function: $\approx 95\mathrm{\%}$ (for $k = 2$); $\approx 85\mathrm{\%}$ (for $k = 3$); and $\approx 75\mathrm{\%}$ (for $k = 4$). However, obtaining such values requires significant computation times, which in the considered case are greater than in the ILOG environment (see Figure 13). It is also worth noting that the calculation times in both analyzed cases (C2 and C3) are characterized by a very low standard deviation value ($\sigma < 2$ for C2 and $\sigma < 25$ for C3; see Table A1), which means that despite the randomly generated population, the solution is usually obtained at the same time.

The presented experiments demonstrate that the planning of reactive missions using accurate methods (e.g., in the ILOG environment) is limited to small-scale networks (up to 100–220 nodes). This range can be increased using the proposed GA (implementing stopping condition C2); the solution is obtained on-line (<600 s). However, this result coincides with a lower value of the target function (e.g., at the level of 60%–80% of optimal solutions).

In further experiments, the influence of weather conditions on the time and quality of the solutions obtained was examined. Figures 14 and 15 show the results of experiments carried out for weather conditions $\mathcal{F}\left(\theta \right) = 10\frac{m}{s}$ and $\mathcal{F}\left(\theta \right) = 11\frac{m}{s}$, in which a disturbance in the form of changes in wind speed ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right) = 11$ $\frac{m}{s}$ oraz ${\mathcal{F}}^{\mathcal{*}}\left(\theta \right) = 12$ $\frac{m}{s}$, respectively, was considered. The results of the experiments are collected in Tables A2 and A3 (see Appendix A).

Figures 14 and 15 compare the solutions obtained in the ILOG and GA environments with the stopping condition C2. It is evident that the change in weather conditions affects the time it takes for the ILOG environment to obtain a solution:

●  for $\mathcal{F}\left(\theta \right) = 10\frac{m}{s},$ the scale of problems for which it is possible to determine a solution on-line ($< 600$ s) amounts to: $n = 200$ for $k = 2$; $n = 120$ for $k = 3$ and $n = 80$ for $k = 4$;

●  for $\mathcal{F}\left(\theta \right) = 11\frac{m}{s},$ the scale of problems for which it is possible to determine a solution on-line ($< 600$ s) amounts to: $n = 200$ for $k = 2$; $n = 100$ for $k = 3$, $n = 80$ for $k = 4$.

It is also apparent that the change in weather conditions has increased the GA calculation time to 500 s. However, the change in weather conditions does not affect the quality of the solutions obtained. The value of the target function for the missions set by the GA remains at a similar level as before:

●  for $\mathcal{F}\left(\theta \right) = 10\frac{m}{s}$, the value of the objective function is equal to: $\approx 78\%$ (for $k = 2$); $\approx 68\%$ (for $k = 3$); $\approx 60\%$ (for $k = 4$);

●  for $\mathcal{F}\left(\theta \right) = 11\frac{m}{s}$, the value of the objective function is equal to: $\approx 75\%$ (for $k = 2$); $\approx 65\%$ (for $k = 3$); $\approx 60\%$ (for $k = 4$).

The conducted experiments indicate that the greatest impact on the time of determining the solution (reactive mission) is primarily the number of nodes $n$ of the network and the size $K$ of the UAV fleet.

8.   Conclusions

We propose a reactive routing method to solve the problem of UAV fleet mission contingency planning in a dynamically changing environment. We have considered plans for UAV route missions in the event that the weather changes beyond the previously predicted situation and/or that the previously agreed order fulfillment terms change. The need to react in such situations enforces the establishment of condition-action rules that allow for the designation of appropriate possible end-to-end routes and enabling safe completion of the mission or its continuation in a modified version during an emergency scenario.

The main advantage of the proposed model is its open structure, which allows for taking into account several variables and constraints, particularly the conditions, enabling a significant acceleration of calculations related to the variants of UAV routes caused by changes in weather conditions.

The developed model was implemented in the IBM ILOG declarative programming environment and the author's Genetic Algorithm. Computational results show that the proposed approach is suitable for online applications. The best results (i.e., the shortest time to determine the solution) were obtained for the GA in which the C2 stopping condition was applied (the value of the objective function has not changed through $q$ generations). Experiments have shown that the scale of problems for which on-line, reactive UAV missions can be determined includes networks containing $n < 220$ nodes and fleets $K\le 4$ UAVs.

In the general case, however, it should be noted that the aforementioned advantage of heuristic methods is paid for by the lack of guarantee that the trajectories of the UAVs motion determined with their help are collision-free. This disadvantage does not have an exact method that implements CP techniques, so it allows planning UAVs routes carrying out their missions at the same flight ceiling.

In our future research, we want to take into account the uncertain nature of the real-world variables which are not deterministic. Thus, a fuzzy approach could be applied to the UAV mission planning problem while allowing for a more accurate estimation of the timeliness of deliveries.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix A