A two-time-level model for mission and flight planning of an inhomogeneous fleet of unmanned aerial vehicles

We consider the mission and flight planning problem for an inhomogeneous fleet of unmanned aerial vehicles (UAVs). Therein, the mission planning problem of assigning targets to a fleet of UAVs and the flight planning problem of finding optimal flight trajectories between a given set of waypoints are combined into one model and solved simultaneously. Thus, trajectories of an inhomogeneous fleet of UAVs have to be specified such that the sum of waypoint-related scores is maximized, considering technical and environmental constraints. Several aspects of an existing basic model are expanded to achieve a more detailed solution. A two-level time grid approach is presented to smooth the computed trajectories. The three-dimensional mission area can contain convex-shaped restricted airspaces and convex subareas where wind affects the flight trajectories. Furthermore, the flight dynamics are related to the mass change, due to fuel consumption, and the operating range of every UAV is altitude-dependent. A class of benchmark instances for collision avoidance is adapted and expanded to fit our model and we prove an upper bound on its objective value. Finally, the presented features and results are tested and discussed on several test instances using GUROBI as a state-of-the-art numerical solver.


Introduction
Due to their variety and flexibility, unmanned aerial vehicles (UAVs) have many possible applications. Next to the long-studied military use [6], also many companies make efforts to incorporate them into their processes, e.g., in parcel delivery [26], but for an efficient and autonomous use, it is crucial to plan the considered task carefully since not only the technical parameters of the used UAV but also the weather and maybe other airborne UAVs must be taken into account. Furthermore, the mission area can contain airports, power plants, or mountains, restricting the airspace and the possible routes. Incorporating all these conditions into the planning process can make the resulting problem very intricate to solve.
The flight planning problem for a given number of inhomogeneous UAVs asks to calculate a flight trajectory between a set of given waypoints for any considered UAV, complying with its technical parameters and the related flight dynamics. The mission planning problem for an inhomogeneous fleet of UAVs can be seen as a version of the Team Orienteering Problem with time windows (TOPTW). Therein, a group of participants wants to maximize its score by visiting waypoints from a given set, each equipped with a score and a time window. The complexity arises from the limited number of paths that can be generated and that there are usually far more waypoints than participants. Since this problem is a variant of the NP-hard Orienteering Problem (OP) [3], it is difficult to solve instances to proven global optimality. A survey about the OP and its variants, also mentioning the TOPTW, can be found in [29].
We consider the mission and flight planning problem for an inhomogeneous fleet of fixed-winged UAVs as a variant of the TOPTW, where further constraints regarding the flight dynamics of the considered UAVs and the environment have to be taken into account. UAVs can navigate through the air almost freely with some related conditions like minimum velocity, maximum altitude, or restricted airspaces. Furthermore, safety distances in aviation are critical. Wind affects the flight of every UAV since their movement is always relative to the surrounding atmosphere. Next to environmental parameters, the mass of a flying object has a strong influence on its flight dynamics. Several characteristics, e.g., maximum acceleration, maximum reachable altitude, or fuel consumption, depend on it [25]. Finally, the operating range has to be observed if the connection of a UAV and its ground control station is not a satellite link, but a UHF/VHF connection instead.
In the literature, the mission planning problem for UAVs applies to a large number of practical tasks, e.g., military operations [17], the observation of icebergs [1], or the reconstruction of terrain from two-dimensional data [28], but often, the calculation of the related flight trajectories is simplified or neglected. An extensive survey on the evolving field of civil applications for UAVs and possible optimization approaches can be found in [20]. Another comprehensive review with a theoretical focus is given in [7], highlighting different modeling approaches for the trajectory optimization problem and providing a taxonomy of 70 papers into the discussed methods by a large set of attributes. Furthermore, several directions for future research are identified therein.
Due to the great interest in the autonomous use of UAVs and the large number of related optimization approaches, we only highlight some recent publications, where mission and flight planning is considered. Ramirez et al. [21] use a genetic algorithm to assign several tasks to a fleet of inhomogeneous UAVs considering technical parameters of the different types of UAVs, restricted airspaces, and time windows. Zhen, Xing, and Gao [32] combine an ant colony algorithm for assigning missions to a set of homogeneous UAVs with a Dubins curve to generate the trajectory of every UAV. The resulting paths are affected by technical parameters of the considered UAV, collision avoidance constraints, and restricted airspaces. Ribeiro et al. [24] use a mixed-integer linear program (MILP) to organize the observation of conveyor belts in a mining site. Their approach is rather combinatorial without exact trajectories but incorporates energy consumption and the possibility to place and use charging stations. Li et al. [18] apply an ant colony algorithm combined with the metropolis criterion to the problem of finding optimal trajectories for a fleet of homogeneous UAVs on a grid map, regarding a given safety distance between the UAVs and given obstacles. Glock and Meyer [11] present a neighborhood search algorithm to assign possible sampling locations after fire or chemical incidents to a set of UAVs and plan their trajectories to maximize the collected information in a given time horizon. Flight dynamics, i.e., maximum velocity, acceleration, and battery level of the homogeneous fleet, are translated into travel times between target locations. Coutinho, Fliege, and Battarra [8] formulate the problem of surveying a set of locations in the aftermath of a disaster by pilotless gliders as a mixed-integer nonlinear problem (MINLP). By applying several linearization techniques, they achieve an MILP computing trajectories for several aircrafts from a starting point to a set of possible landing areas within an obstacle-free airspace assuming constant weather conditions. Cheng, Adulyasak, and Rousseau [5] derive a branch-and-cut algorithm to solve an MINLP describing the problem of multitrip parcel delivery by a homogeneous fleet of UAVs. Considering constant altitude and velocity, their energy consumption is modeled using a mass-dependent, nonlinear function, and trajectories are given by distances and travel times. Thibbotuwawa et al. [27] set up an MILP to plan the supply of a set of customers with an inhomogeneous fleet of UAVs under weather uncertainty. Their proposed model incorporates several weather zones, each with related wind conditions, energy consumption, and collision avoidance. In terms of the flight dynamics, it lacks the exact trajectory, only computing the sequence of visited customers. Kai et al. [15] give an MILP to plan the trajectory of a UAV visiting a given set of waypoints taking into account detailed flight dynamics. In terms of environmental constraints, the therein presented model lacks weather conditions and restricted airspaces. Xia, Wang, and Wang [30] formulate the problem of controlling moving ships in emission control areas by UAVs, stationed at the coast, as an MILP and compare it with a Lagrangian relaxation. The UAVs are assigned to computed waypoints, where they meet a vessel. Though, besides the battery level, no flight dynamics are taken into account. Albert, Leira, and Imsland [1] present an integer linear optimization model for the observation of drifting icebergs in arctic areas to support shipping using a homogeneous fleet of UAVs deployed at a ship. Their approach uses Dubins curves to compute trajectories and can update the present path during the mission by a new run of the optimization model. But therefore, it neglects most of the flight dynamics to speed up the solution process. Chen et al. [4] consider a UAV as particle affected by different force fields and plan trajectories avoiding given obstacles. The resulting problem is solved using optimal control theory but neglects flight dynamics and weather conditions. Within the flight planning process, collision avoidance and conflict resolution is a crucial field to ensure the safe operation of aircrafts and avoid significant harm. It considers a given number of aircraft within the same airspace and asks for the best control changes for each of them to achieve collision-free trajectories. If the problem is stated in three dimensions, also altitude changes are possible to resolve conflicts. Although the process of solving potential conflicts between midair aircraft is automated [16], considering this aspect already in the operation planning phase can reduce the number of conflicts and errors later.
An extensive review about the field of conflict detection and resolution and related solution methods can be found in [23]. Due to its currency, we further mention only another forthcoming approach. In Hoch et al. [13], consider collision avoidance as a non-stop disjoint trajectory problem, where commodities are shipped through a timeexpanded graph without coming across each other. They study the general properties of this problem and examine different graph classes which can be applied as takeoff/landing phase with several aircraft or airspaces with a fixed number of airports and aircraft. Due to the underlying network, this approach lacks detailed flight trajectories.
We contribute to the state-of-the-art by considering the mission and flight planning problem for an inhomogeneous fleet of UAVs in a mixed-integer non-linear programming framework, including detailed flight dynamics of the different UAVs, the effect of wind in different subareas and altitudes, convex shaped restricted airspaces, and collision avoidance. Furthermore, we extend several of its aspects to enhance the quality of the computed solution. The novelty of our approach is the combination of assigning waypoints to the participating UAVs and calculating detailed flight trajectories at the same time. This paper extends the work presented in [10]. The newly incorporated aspects are the following: By the combination of two different time discretizations, the computed flight trajectories have an increased level of detail. We incorporate massdependent maximum velocity, maximum acceleration, maximum reachable altitude, and a mass-dependent fuel consumption into the presented model. The operating range of every UAV is defined dependent on altitude, assuming no sky propagation, and the height of the used antenna. Restricted airspaces are generalized to convex shaped areas and wind affects any UAV within specified convex subareas. Considering collision avoidance, we adapt and extend a class of benchmark problems to fit our model and derive upper bounds for their necessary amount of discrete time steps in the twoand three-dimensional case.
This paper is structured as follows. In Sect. 2, the basic model is set up by applying the two-level time discretization. Furthermore, the computation of the UAVs flight dynamics, the operating range, the restricted airspaces, and the influence of wind are expanded to allow more realistic constraints and achieve high-quality solutions. In Sect. 3, a class of benchmark instances for collision avoidance is adapted to the introduced model to study the efficiency of the collision avoidance constraints and upper bounds for the necessary time horizon are given. The effect of the derived modifications and results is tested in several scenarios in Sect. 4, including the influence of the time discretization on the computation time, comparison between the derived upper bounds and the computed optimal solution for collision avoidance problems, and a detailed discussion of particular aspects of our model. We summarize the presented results and give directions for future work in Sect. 5.

Two-time-level model
If the velocity and the acceleration can change only at some discrete time step, the step length has to decrease to smoothen the computed trajectories. Thus, more time steps are necessary to cover the same time interval, leading to more complex instances and slower solution processes. In the following, we present a possibility to achieve smooth trajectories by using two different, coupled time discretizations and apply it to the model in [10]. Afterwards, the derived two-level time discretization model is extended in several ways. A list of all used sets can be found in Table 1. The parameters and variables within the two-level time discretization model are given in Tables 2 and  3, while the additional parameters and variables of the complete extended model are described in Tables 4 and 5. Note, that for vector valued parameters and variables the first component refers to the x-direction, the second to the y-direction, and the third to the z-direction, respectively. The components of a vector valued parameter or variable are marked by their respective coordinate direction as a superscript, e.g., the position vector r contains the components (r x , r y , r z ) T .
The time horizon [0, t · T ], for a given number of time steps T and time step length t, is discretized in two different ways to achieve a better discretization without an increase in the number of time steps and discrete decision variables. Let n f denote the number of fine time steps between two adjacent time steps in {0, 1, . . . , T }. For the first discretization, the step size t is applied, but the set T contains all multiples of n f + 1 in the interval 0, (n f + 1)T , i.e., T = {0, n f + 1, . . . , (n f + 1)T }. This discretization is applied to the binary variables to model the possible decisions. Note, that |T | is constant for varying values of n f . Then, we embed T in the superset . , (n f + 1)T } with the corresponding step length t f = t n f +1 to cover the same time horizon. This fine time approximation is applied to all continuous variables to refine the computed trajectories by controlling n f . Furthermore, a function t : T f → T of the form  Figure 1 holds an illustration of this two-level time grid approach.
A detailed discussion of the following constraints can be found in [10]. Furthermore, we use the following settings assumed therein: -An inhomogeneous fleet of fixed-winged, fuel-driven UAVs is considered  -Velocities are calculated using the TC-TH-W vector triangle, where TH is the true heading vector, W is the wind vector, and by vector addition obtained is the true course vector TC=TH+W -Every UAV u ∈ U chooses the highest allowed velocity θ u, j for every velocity band j ∈ V u since the fuel consumption rate η u,i, j depends only on the considered velocity band and not on the chosen velocity within it.
Applying the time step sets T and T f to the respective group of variables and the mapping (1) to match them, the mission and flight planning problem, incorporating the two-time-level discretization, is given by where . 2 and . 3 denote the two and three-dimensional Euclidean norm, respectively. Note, that the components for maximizing altitude and minimizing velocity in (3) are not necessarily required to minimize the fuel usage, but our experience showed that it is beneficial for the performance of the simplex method to incorporate all variables with nonzero coefficients in the objective function, if they affect the optimization goal. Constraints (2.1) and (2.2) fix the respective start and end location of every UAV, while (2.3) limits their maximum operating range. The maximum and minimum altitude of the considered UAVs is restricted by (2.4). Equations (2.5) and (2.6) are a discretization of the Newton's equations of motion to compute the exact flight trajectory together with the altitude changes in (2.7). The UAVs velocity and acceleration is bounded by (2.8) and (2.9). In (2.10)-(2.12) it is ensured, that any UAV can only take off once and remain on the ground after landing. The visit of the given waypoints is managed by (2.13) and (2.14), guaranteeing that the respective UAV is inside the required operational range and every waypoint is only visited once. Constraints (2.15) and (2.16) let a UAV avoid the given restricted airspaces, while (2.17) and (2.18) manage the collision avoidance between any pair of UAVs. The correct computation of the UAVs fuel consumption is ensured by Eqs. (2.19)-(2.23). In the objective function (3), the primary goal is to maximize the score of visited waypoints within the given time horizon. At the same time, on a subordinate level the used trajectory should be most economical in terms of fuel by reducing the airborne time and the used amount of fuel, choosing low velocity, and flying in high altitude. All of these aspects reduce fuel on their own, but their importance can be weighted arbitrarily by choosing appropriate values for the constants M air , M f uel , M alt , and M vel .
By increasing the number of fine time steps n f , the calculation of the flight dynamics in model (2), consisting of the constraints (2.1)-(2.23), is done more often between two adjacent decision time steps in T . This gives the possibility to achieve more detailed flight trajectories without increasing the number of discrete decision variables and leads to a less complex model compared to an approach using a single time discretization for all variables with the same level of detail. The Euclidean norms in the presented model have to be linearized to achieve a mixed-integer linear problem. We achieve this according to the methods shown in [10].

Flight dynamics
In the given model, the left part of Eq. (2.8) forces every UAV to fly with at least minimum velocity as soon as it takes off. If one of the deployed UAVs cannot accelerate to its minimum velocity within a single time step, it would have to stay on the ground and thus would not be used for mission planning. To overcome this drawback, we define a new binary variable b u (t), indicating whether the minimum velocity has to hold.

and only if the UAV has minimum velocity.
Proof Due to constraint (2.6), the UAVs velocity is a linear function of its acceleration. Thus, the real time for reaching minimum velocity v u with maximum acceleration a u is given by v u a u . Conversion into a number of discrete time steps t u ∈ N leads to To ensure that b u (t) is set to one t u time steps after the UAV starts to move or before it attempts to land, we relate it to b u (t).
To relate all 2 t u + 1 binary variables b u (t ) in this interval to a single binary decision, 2t u is subtracted leading to the desired Eq. (2.24).
To incorporate the new variable, the minimum velocity has to hold only t u time steps after takeoff. Thus, constraint (2.8) Accordingly, the fuel consumption Eq. (2.21) is adjusted to the new variable by Otherwise, Eqs. (2.21) and (2.22) would force any UAV that cannot accelerate to its minimum velocity in one timestep to stay on the ground. Linking the fuel consumption to the new variables b u (t) by (2.21') may lead to a underestimation of the burned fuel during the t u time steps for the takeoff and the landing phase, where the UAV u ∈ U is possibly moving with a velocity lower than its minimum velocity while it holds b u (t) = 0. The fuel surplus for climbing ξ u depends only on the UAVs current climb rate and thus is not affected by this changes. By defining ξ u for every altitude band and adjusting its value for low altitudes, the underestimation of burned fuel during the takeoff could be compensated. For descending, there is no fuel used due to the assumption of descending by the force of gravity. Due to the physics of flight, a UAV slower than minimum velocity cannot ascend and therefore not take off. Thus, similar to the minimum velocity, the left part of Eq. (2.4) eliminates any UAV from the planning process if it cannot accelerate to its minimum velocity and ascend to the minimum altitude within a single time step. To avoid this, the variable b u (t) also is incorporated into Eq. (2.4), leading to We consider the minimum altitude small enough such that any UAV with minimum velocity can ascend to it within a fine time step. Another occurring phenomenon in aviation is the altitude-and velocity-dependent change of flight dynamics due to decreasing thrust in higher altitudes caused by decreasing air temperature and pressure. It affects the maximum velocity as well as the maximum climb and descent rates. For the maximum velocity, altitude dependence is incorporated into the model by defining a maximum velocity v u,i of UAV u ∈ U for every of its altitude bands i ∈ L u and replace Eq.
Besides this family of constraints, an altitude-dependent maximum velocity has to be taken into account in the calculation of the fuel consumption. By redefining parameter θ u, j to θ u,i, j , describing the velocity of UAV u ∈ U in throttle band j ∈ V u and altitude band i ∈ L u , Eq. (2.22) changes to For climb and descent rates, the velocity has to be taken into account, next to the altitude. In a similar way to the velocity, the related parameters change to v z,+ u,i, j and v z,− u,i, j , describing the maximum climb and descent rate of UAV u ∈ U for altitude band i ∈ L u and throttle band j ∈ V u , respectively. Thereby, the upper bounds for the variables v z,+ u (t) and v z,− u (t) are no longer constant. Thus, the new constraint for the upper limit of the descent rate is given by For the climb rate, the approach is a bit different since a constraint similar to (2.25) would reject any takeoff. Therefore, a new parameter v z,+,0 u is introduced for the maximum climb rate at takeoff for every UAV u ∈ U and it is adjusted to get the maximum midair climb rates v z,+ u . The constraint, limiting the climb rate for every UAV, then has the form In model (2), the maximum acceleration of the UAVs is assumed to be constant in every situation. Thus, the computed flight trajectories are inaccurate since in practice the maximally allowed acceleration at the takeoff is significantly higher than during the flight. This results either in underestimation during the takeoff phase or overestimation when in midair. We follow the approach given in the Base of Aircraft Data (BADA) [19] of the Eurocontrol Experimental Centre. Therein, a reduction factor from takeoff thrust to cruise thrust applies with a value of C T cr = 0.95. It is incorporated into the model by changing Eq. (2.9) to

Fuel-dependent flight dynamics
The flight dynamics in model (2) are independent of the mass of the considered UAVs, but this assumption only holds for electrical systems, where the weight of the battery is constant. For UAVs with liquid fuels, the mass affects the maximum velocity, maximum acceleration, maximum altitude, and fuel consumption rate [25] of an aircraft by weight reduction due to burned fuel. Since the portion of the fuel on the gross take-off weight of a UAV can go up to 43% [9], this has a great impact on the mission planning scenario. To include these fuel-dependent flight dynamics into the proposed model, several changes are necessary. In terms of the altitude, we assume an affine-linear dependency between the amount of remaining fuel and the reachable altitude. We expect the empty weight m u of UAV u to include a small fuel reserve. Let UAV u ∈ U reach altitude ϕ alt,2 u h u for its maximum amount of fuel F u and its maximum altitude h u for ϕ alt,1 The resulting linear function is displayed in Fig. 2 and gives the constraint Besides, the right part of constraint (2.4') is still necessary to limit the reachable altitude to its maximum value and set the maximum altitude to zero when the UAV is on the ground.
To include the mass of the remaining fuel into the calculation of the velocity and the acceleration, we have the following proposition.
minimizes the quadratic error of the mass-dependent maximum velocity and acceleration, respectively.
Proof We compare the given maximum velocity v u related to m u with another velocitỹ v u for a UAV with more remaining fuel. Since the generated thrust of UAV u ∈ U is limited, it can reach a certain kinetic energy E max u = 1 2 m u v 2 u for a given time. Thus, by energy preservation, we get the relation Using the equation a = v t , in the same way the expressioñ can be derived. Since we consider the maximum velocity and acceleration in (5) and (6), we restrict on the positive solution of the respective square root. The nonlinear factor for Eqs. (5) and (6) is the same. Thus, we only have to perform one least squares approximation of the form where we fix the absolute term of the linear function to 1 to ensure that the maximum velocity and acceleration is attained according to our assumption. The factor 1 m u +F u is applied to have a greater numerical value for ϕ acc u and thus avoid numerical issues. The solution of (7) gives the desired parameter ϕ acc u for the obtainable percentage of the maximum velocity and acceleration per unit of fuel for UAV u ∈ U. Replacing the factor b u ( t ) − 0.05b u ( t ) by the computed linear function ϕ acc u g u (t) F u +m u + 1 in (2.9'), it gives the desired mass-dependent maximum acceleration (2.29) for every airborne UAV u ∈ U in time step t ∈ T f . In terms of the mass-depedent maximum velocity, the above approach is used and slightly modified since the altitude dependence in (2.8") calls for the incorporation of an sufficiently large additive constant V u 1 − j∈V u s u,i, j ( t ) with V u = max i∈L u v u,i to avoid nonlinearities. This leads to the remaining Eq. (2.28).
Note, that the used approach has no underlying physical model and is just one possibility to incorporate measured and approximated data into the problem. It is similar to the method used by Eurocontrol [19], a European organization for air traffic management, except that we assume the UAV to reach maximum velocity and acceleration with only the fuel reserve incorporated in the empty weight. This assumption also yields −1 ≤ ϕ acc u ≤ 0 for F u ≥ 1 since in (7) a linear approximation of a monotonically falling function with f (0) = 1 on the interval [0, F u ] is computed for every UAV u ∈ U.
In terms of the fuel-dependent fuel consumption, a further parameter ϕ f uel u is necessary. It describes the additional fuel consumption for every remaining fuel unit at the previous time step and is incorporated as an additional factor into the constraint (2.20). To ensure that a UAV is consuming no fuel when landed, a distinction between the fuel calculation for a landed and a flying UAV is needed. Thus, the constraint (2.20) is replaced by Note, that in addition to Eqs. (2.28) and (2.29) the constraints (2.8") and (2.9') are still necessary to force the maximum velocity and acceleration to zero for grounded UAVs. They are the only remaining constraints incorporating the variables b u (t) and thus ensure the correct calculation of the fuel consumption in (2.20i) -(2.20iv). The new introduced variables b u (t) cannot be used here since they would neglect the fuel consumption for ascending at the takeoff phase before reaching the minimum velocity.
Applying these changes to the model, the flight dynamics of the UAVs are compromised since Eq. (2.19) ensures that every UAV takes off with maximum fuel, reducing the maximum altitude, velocity, and acceleration. Thus, Eq. (2.19) is dropped to make the amount of takeoff fuel variable and allow the UAVs to reduce the amount of initial fuel to achieve better flight dynamics. To ensure these fuel-saving aspects, the objective function (3) has to change to max u∈U ,w∈W,t∈T ω

Range
In model (2), the maximum operating range u is chosen constant, although it depends on the height A u of the used antenna and the altitude r z u (t) of the UAV. Assuming UHF/VHF control and no satellite link, the UAV must ensure a signal link to the ground station, represented by the antenna. Thus, for large distances, it has to fly above a certain altitude to establish its connection, due to the curvature of the earth. To achieve a more realistic setting, we deduce an altitude-dependent maximum range.

Proposition 3 Let C be a constant with C ≥ max u∈U h u . Then the linear approximation
minimizes the quadratic error to the altitude-dependent maximum operating range at time step t ∈ T f .

Proof
The maximum range of UAV u ∈ U can be displayed as in Fig. 3.
Applying the Pythagorean theorem, the maximum range is calculated by where E is the radius of the earth and A u is the height of the antenna. Due to the nonlinearity of this expression, we use the continuous least squares method to find the best approximation in terms of a linear function f (x) = alt x + init . Since the antenna height A u is constant, the first term in (9) is an additive constant and only shifts the whole function. Thus it can be neglected in the least squares approximation, leading to the problem with C ≥ max u∈U h u representing a sufficiently large altitude as upper limit of the integral. By integration and solution of the remaining minimization problem, the parameters alt and init are given by (8). Finally, the constant maximum operating range u in (2.3) is replaced by the linear function alt r z u (t) +˜ init to obtain (2.3').

Restricted airspaces
Equations (2.15) and (2.16) in model (2) allow to incorporate restricted airspaces for the UAVs, but only cubic ones which are parallel to the coordinate axes. To compare the respective coordinates of UAV u ∈ U and the restricted airspace q ∈ Q at time step t ∈ T and report if it is inside the area, two binary variables f i u,q f i u,q are necessary for every coordinate direction i ∈ {x, y, z}, one per side. Using this approach, also more complex restricted airspaces can be modeled by a union of these cubic ones, but only at the expense of many additional binary variables.
To describe restricted airspaces by arbitrary polyhedrons, we assume its N Q q describing hyperplanes to have coordinate form with normal vector c q,i (t) and righthand side c rhs q,i (t) for i = 1, . . . , N Q q and q ∈ Q. Thus, every restricted airspace q ∈ Q is represented by Similar to the parameters, the related binary variables providing a more general approach to model restricted convex airspaces with again one binary variable per side.

Wind
In the basic model (2), the wind is assumed to blow constantly everywhere within the considered area. By introducing the fine time steps, the wind can change in every fine time step t ∈ T f . To allow a more accurate influence of weather conditions, we want to include the possibility of different wind zones within the mission area. But this comes with the restriction that wind changes only in time steps t ∈ T . For simplicity, the computation of the position (2.5) is modified to wherew i u (t) is the influence of wind for UAV u ∈ U in coordinate direction i ∈ {x, y} at time step t ∈ T − f . As a first approach, altitude-dependent wind zones could be added since the wind intensity and its direction can change with increasing altitude, e.g., the jetstream has to be taken into account for mission planning in altitudes between 8 to 12 kilometers. To incorporate this into our model, the wind vector w l (t) = w x l (t), w y l (t) ∈ R 2 is now defined for every altitude band l ∈ L u . Since the variable s u,i, j (t) cannot change in every fine step, there is a discretization error between two time steps if the altitude layer is changed. To reduce its influence, we consider a convex combination of the wind between time steps t, t + n f + 1 ∈ T at every fine time step t ∈ T f . Thus, the influence of wind is calculated bỹ with μ = t mod (n f +1) n f +1 and t = t + n f + 1. But this method requires again constant wind direction within every altitude band over the whole mission area.
To overcome this drawback, a second approach is possible. In the same way as for the restricted airspaces in Sect. 2.4, a set P of wind zones is defined, with each p ∈ P described by N P p hyperplanes with normal vectors n p,i (t) and right-hand sides n rhs p,i (t). Note, that this approach allows the coverage of the mission area by a set of polyhedral non-overlapping wind zones of arbitrary size to incorporate the wind in great detail, if necessary.
Furthermore, new binary variables w u, p,i (t) ∈ {0, 1} are introduced with w u, p,i (t) = 1, if UAV u ∈ U is within wind zone p ∈ P regarding hyperplane i ∈ A P p at time step t ∈ T . In contrast to the restricted airspaces, now additional constraints and an auxiliary binary variable ω u, p (t) are necessary to map the variables w u, p,i (t) to a single binary decision and to ensure that wind is only taken into account when UAV u ∈ U is in midair. Therefore, the decision if UAV u is affected by wind zones p at time step t is modeled by Analogue to the first approach, the variable w u, p (t) can only change at time step t ∈ T . Thus, the convex combination witht is applied again and the influence of wind is calculated bỹ where again μ = t mod (n f +1) n f +1 and t = t + n f + 1. For this approach, the number of binary variables within the model increases by p∈P (|A P p | + 1) per time step and UAV. The Eqs. (2.30)-(2.34) allow the influence of multiple wind zones to a single UAV at the same time. This is not realistic since in this case completely new weather conditions arise, in contrast to just adding the wind vectors as in (2.34). The user of our model must be aware of this and has to define all wind zones in (13) as pairwise non-overlapping areas.

Collision avoidance
Compliance with safety distances in aviation is crucial due to the disastrous consequences of its failure. Therefore, the parts of a planning model computing collision-free flight trajectories should be examined excessively. In the following, we consider a two dimensional and a three dimensional setting to examine the modelled collision avoidance (2.17) and (2.18).
A class of two-dimensional benchmark instances, called Random Circle Problems (RCPs) [22], has a wide application in collision avoidance and conflict resolution. In this problem class, a set of aircrafts is randomly arranged on a circle. At the beginning, all UAVs maintain the safety distances, head towards the circles center, and have their respective destination on its opposite side. The deviation from the shortest connection is minimized for every participant. With an increasing number of aircrafts, it is a challenging problem to conflict resolution approaches due to its fast increasing violations of the safety distance if no countermeasures are taken.

Two-dimensional model
We adapt these benchmark instances to our setting, assuming every participating UAV has the same technical characteristics. Thus, there are no waypoints and the start locations of the participating UAVs are distributed randomly along a circle of radius r RC P , with their respective end location on the opposite side. The operating range is neglected since every UAV must be able to reach its end location. Thus, the whole circle is within its operating range. In the beginning, every UAV is heading towards the center of the circle and all safety distances are ensured. To incorporate this into the model, a new constraint is added, where C is the location of the center of the considered circle and the new variable γ u scales the given direction to ensure the limits of the velocity. Due to the absence of waypoints and the altitude, many constraints can be neglected and also many binary variables are unnecessary or must be redefined, i.e., s u,i, j (t) changes to s u, j (t) indicating whether UAV u ∈ U is within throttle band j ∈ V 1 at time step t ∈ T and the variables b u (t) and b u (t) are substituted by b − u (t) since all UAVs start in midair and have to stop at their end location.
We minimize the deviation of the shortest connection by the number of necessary time steps T min , such that every UAV reaches its respective end location. This number is not known a priori and increases the more UAVs are participating. Only a lower bound T min is given by the number of time steps it takes for a single UAV to fly along the diameter of the circle with maximum velocity. Setting T to a sufficiently large value and read T min from the optimal solution can be an ineffective way since the number of constraints and variables is time-dependent, leading to a possibly too large model. In [2], the same situation occurs for a scheduling problem with the objective to minimize the makespan and it is solved by an sequential solution approach increasing the number of considered time steps. We follow this approach, solving the RCP sequentially for an increasing number of time steps, starting with T = T min , until a feasible solution is found or its existence cannot be determined within a given computation time. This procedure has the advantage of working without any objective function since the first found feasible solution gives the value of T min .
To ensure that the safety distances are observed at every time step, it has to hold where V 1 = max i∈L 1 v 1,i is again the maximum of the altitude-dependent maximum velocities of the first UAV. Due to the assumption of identical technical characteristics of all participating UAVs, the first UAV can be chosen without loss of generality. Without condition (15), a pair of UAVs, complying with the safety distance between them, could swap their positions within one time step by flying through each other. The resulting resolution of the time is sufficient to obtain smooth trajectories. Thus, we renounce on fine time steps in this case.
The following result gives a heuristical solution for the two-dimensional RCP.

Proposition 4
Consider the two-dimensional RCP with center C ∈ R 2 , radius r RC P ∈ R + , safety distances ε ∈ R 2 + , time step length t ∈ R + , n ∈ N participating UAVs, and v ∈ R + their maximum velocity. Then it holds Proof Every UAV starts at its random starting location and heads towards the center C with maximum speed. To avoid conflicts near the center C, all UAVs are assigned to a smaller circle with center C and radius r ub ≤ r RC P before the first conflict would occur. Along this smaller circle, the UAVs rotate around its center on a semicircle and then fly to their respective end locations. This approach generates a feasible solution if the safety distances are observed during the rotation. Therefore, it is necessary to ensure at least the distance ε between any pair of UAVs. Let denote the set of time steps for which in the next one the safety distance is violated by some pair of UAVs. Then t = min {t | t ∈ T ε } is the last time step maintaining all safety distances and the radius r ub is given by C − r n (t ) , where the position of any UAV at time step t can be chosen since they all have equal distance to the center. Without loss of generality, we take the position of UAV n. For the described solution, every UAV must fly a distance of 2(r RC P −r ub ) + πr ub units to reach its end location. Thus, there are time steps necessary to complete this trajectory. Substitution of r ub leads to the upper bound of T min given in (16).
Furthermore, the quality of the described heuristic solution can be calculated.

Corollary 1 For the two-dimensional RCP with minimum number of time steps T min , it holds
Proof In the worst case, there are two UAVs u 1 , u 2 ∈ U starting with R 0 u 1 − R 0 u 2 = ε . Then it follows r ub = r RC P and the trajectory of every UAV is πr RC P . Since all UAVs have the same constant maximum velocity v, it holds T min = 2r RC P v and the number of time steps T min depends only on the radius r RC P , leading to T min = πr RC P v = π 2 T min . Thus, in a general setting the desired estimate (20) is obtained by combining the definition of T min with the described worst case and rearranging it.
For the case of evenly distributed UAVs, the result of proposition 4 can be computed without the knowledge of t .

Corollary 2
Consider the Circle Problem with n ∈ N evenly distributed UAVs with maximum velocity v ∈ R + , radius r RC P ∈ R + , safety distances ε ∈ R 2 + , and time step length t ∈ R + . Then it holds Proof In the evenly distributed case, every pair of UAVs has the same distance. Thus, their start locations are the corners of a regular polygon with circumradius r RC P . Since all UAVs head towards the center with the same velocity, they stay corner points of a regular polyhedron with decreasing side length. The radius r ub of the smaller circle is then the circumradius of a regular polyhedron with n corners and side length ε. It is computed by Substituting this term into (19) gives the upper bound (21).
The model in the two-dimensional case is given by the constraints ( (14), where all terms with v z,+ u (t) are neglected, due to the absence of altitude changes. Furthermore, the variables b u (t) and b u (t) have to be substituted by b − u (t) since all UAVs start in midair and should only stop at their respective end locations.

Three-dimensional model
Since in ordinary air traffic the aircrafts are not restricted to a single altitude and to fully evolve the potential of our model, we extend the concept of RCPs into three dimensions. Therefore, minimum and maximum altitudes h and h are incorporated as lower and upper bound of the altitude r z u (t) of every UAV u ∈ U at time step t ∈ T f , respectively. All UAVs are positioned analog to the two-dimensional case, but with initial altitude R 0,z ∈ h, h .
The result of proposition 4 can also be adapted for the three-dimensional RCP.

Proposition 5
Consider the three-dimensional RCP with center C ∈ R 3 , radius r RC P ∈ R + , minimum and maximum altitude h and h, safety distances ε ∈ R 3 + , and time step length t ∈ R + . Let n l = h−h ε z and the sets U 1 , . . . , U n l be a partition of U = {1, . . . , n}. Furthermore, every participating UAV has maximum velocity v ∈ R + , maximum climb and descent rate v z,+ and v z,− , and initial altitude R 0,z . Then it holds with Proof For the given altitude range h, h and vertical safety distance ε z , it is possible UAVs one above each other at the same x-and y-coordinates. Thus, every partition {U i } i∈{1,...,n l } of the set U = {1, . . . , n} decomposes the threedimensional RCP into n l two-dimensional RCPs considering only the UAVs U i , respectively. According to this, the number of time steps to perform the threedimensional trajectory of every UAV u is the sum of the number of time steps for its trajectory in the two-dimensional problem and the number of time steps necessary for the altitude changes.
To arrange all two-dimensional RPCs into the altitude range h, h , one of them is located at the lower or upper altitude limit and the vertical distance between two of them is at least ε z . Without loss of generality, we assign the problem considering the UAVs of subset U 1 to the altitude h. Then, the necessary number of time steps for the altitude changes of every UAV u ∈ U i are computed by Thus, the number of necessary time steps in the three-dimensional case is computed by with t from proposition 4. The maximum of these T i is the desired upper bound (23).
The model in the three-dimensional case is given by the constraints ( Similar to the two-dimensional case, the variables b u (t) and b u (t) are substituted by b − u (t).

Computational results
To test the derived results, we apply the extended model ( The data for the technical parameters of the UAVs are taken from [9] and can be found in Table 6 to Table 8. As length of one time step, t = 0.1h is assumed. Regarding the maximum operating range, we choose in (8) C range = 15 km, E = 6371 km, and A u = 0.005 km Table 7 Altitude and throttle band related data of UAV-1 according to [9]   for all u ∈ U . Thus, we obtain the optimal values alt := 23.30 and init := 116.61, leading to the linear approximation displayed in Fig. 4. For the altitude and throttle dependend climb rate we assume v z,+ u,i, j = v z,+,0 u for all UAVs u ∈ U, altitude bands i ∈ L u , and throttle bands j ∈ V u .
The operational range δ u,w of UAV u ∈ U to waypoint w ∈ W is set to a value within the second highest altitude band to reject the UAVs staying at maximium altitude all the time. According to this, we choose the values δ 1,w = 10 km and δ 2,w = 3 km.
The maximum acceleration a u of UAV u ∈ U is not provided by [9], but it gives a value for the military static thrust at starting or landing. In these processes, the maximum possible acceleration of the aircraft is used [25], so we use it as an upper bound for every time step. To compute the acceleration related to the given thrust, the second of Newton's axioms F = ma is chosen, rearranged to get a, and the thrust, converted to Newton, is put in. Combined with the conversion factor from m s 2 to km h 2 , The weights M i , i ∈ {air, f uel, vel, alt}, for the subordinate level of the objective function (3') must be chosen such that they do not outperform the visit of a waypoint w ∈ W with minimal score S = min w∈W {S w }. Since the components of the subordinate level in (3') have different signs, this has to hold for both the negative and the positive part. Thus, they must be chosen accordingly to ensuring that visiting a waypoint is always preferred. In this work, we choose M air = 10 2 , M f uel = 10, M vel = 10 5 , M alt = 10 3 . Furthermore, we assume the sufficiently large constant M dist = 10 3 and set S w = 1 for all w ∈ W.
All instances were generated with the modeling language AMPL (Version: 20200501) and solved with GUROBI 9.1.2 [12] on an Apple Mac Pro with an Intel Xeon W running 32 threads parallel at 3.2 GHz clock speed and 768 GB RAM, with the relative and absolute gap set to 0 and default settings otherwise. If a time limit is considered, it is set to 3600 seconds. For objective value z and upper bound z, the relative gap is in the following defined by

Time discretization
In this section, we examine the influence of the two-time-level discretization to the computation time and the solution quality of the problem. With the use of a more detailed approach, the resulting problem becomes more complex and its computation time rises. Since we apply the state-of-the-art solver GUROBI, we find the instance sizes that can be solved within a given time limit.

CPU time analysis
To examine the influence of fine time steps on the solution time, we generate 41 instances with 15 waypoints, a single UAV, 40 time steps, and two restricted airspaces and vary the number of fine time steps from 0 to 9. The results in Fig. 6 show, that the average computation time increases with a larger number of fine time steps. To solve 90% of the given instances within the considered time limit, at most five fine time steps can be used. For nine fine time steps only nine instances were solved within the time limit and also the average relative gap takes its maximum with 18% for the arithmetic mean. The total amount of CPU time for solving all 410 instances was 691,583 seconds.

Solution quality
To study the effect of a varying number of fine time steps on the quality of the solution, we perform two experiments. At first, a comparison of the computed flight trajectories of a single instance for one, three, seven, and 15 fine time steps is shown in Fig. 7. The primal feasibility tolerance, the integer feasibility tolerance, and the optimality tolerance on the reduced costs in GUROBI were set to 10 −4 to avoid numerical issues in the biggest instance. For the problem with one, three, seven, and 15 fine time steps, the solution process took 87 seconds, 537 seconds, 2, 020 seconds, and 13, 289 seconds, respectively.
One can see that the optimal trajectory gets less crooked for an increasing number of fine time steps since the UAV has more possibilities to change its velocity and acceleration. At the same time, the differences between the computed trajectories become smaller. Thus, it suggests to use not too many fine time steps to save computation time. The remaining oscillation of the trajectory, e.g., in Fig. 7b between WP-13 and WP-9, is caused by the linearization of the Euclidean norm taken from [10] since it overestimates the real norm for some points.
In the second experiment, we compare the solution quality for a given level of detail in the absence and presence of fine time steps. Therefore, an instance was solved with T = 200, t = 0.02, and n f = 0 in the case of no fine time steps and with T = 40, t = 0.1, and n f = 4 in presence of them. For the former one, 897 seconds were necessary to compute the proven global optimum, while the solution process for the latter one took 7, 032 seconds. The results are displayed in Fig. 8.
The main differences between both solutions are found in the altitude profile and the fuel consumption. Due to the fact that in the presence of fine time steps the binary decision variables can only change every n f + 1-th time step, the UAV is sometimes restricted to certain decisions unnecessarily long. Regarding the altitude this can be seen during the takeoff and the landing phase, where the UAV stays at the same altitude for some time before ascending or descending further. It also affects the trajectory since the halfspaces described by (2.15') can only be entered or left at a coarse time step. Thus, the UAV can fly around a corner of the restricted air space only when the subsequent time step is in T . Finally, the visit of a waypoint and the landing can be performed only at a coarse time step, limiting the possibilities when the UAV passes a waypoint within sensor range and when it lands. These three aspects result in longer flight duration, force the UAV to fly detours, like its trajectory between WP-1 and its hub, to choose a smaller velocity, as at the beginning of its mission, or to burn additional fuel by flying at a lower altitude to get into a waypoints sensor range, like at WP-3. The latter behavior is also encouraged by the objective function since the flight duration has a higher weight than the fuel consumption. Thus, the solution in the presence of fine time steps is an overestimate of the solution in terms of flight duration and fuel consumption.
Next to the described discretization errors, there is also another phenomenon, called "corner cutting", visible in the trajectories in Fig. 8. Since the check of the position of  the UAV regarding the restricted airspace is done only at discrete time steps, the UAV can traverse these areas within a single time step if it is again outside at the next time step.This problem occurs also in the absence of fine time steps and is a topic of our future work.

Special instances
In the following experiments, we apply the model to specially designed instances to examine some of its aspects in detail. For the effect of wind to the optimal flight trajectories, we compare the optimal solutions of an instance with and without the presence of wind. Furthermore, the altitude-dependent operating range approach the mass-dependent flight dynamics are discussed.

Influence of wind zones
We consider an instance with 10 waypoints, two UAVs, 25 time steps with 4 fine time steps each, and three wind zones. Their respective area and wind velocity are given in Table 9. The first wind zone displays the jetstream on the northern hemisphere, while the other two describe local weather phenomena with a heavy storm to emphasize its influence on the UAVs, although real UAVs would stay on the ground in this situation. In absence of wind, the total amount of CPU time to solve the instance was 181 seconds, while in presence of wind the computation time increased to 2270 seconds. The optimal flight trajectories and the obtained altitudes of the UAVs for both cases of the considered instance are given in Fig. 9.
One can see that, in the presence of wind, the assignment from waypoints to UAVs changes. UAV-1 benefits from all three wind zones. At first, it uses wind zone 1 to travel a longer way to waypoint 8, while the time step for visiting it is the same. Then, it stays below wind zone 1 and uses wind zone 2 to visit the more distanced waypoint 3 at the same time step it visits waypoint 9 in the absence of wind. Finally, wind In the presence of wind, the different wind zones are displayed in different colors zone 3 is advantageous on the way to the end location. UAV-2 also profits from the wind zones on the way to its first visited waypoint and when it flys to waypoint 2.
In the presence of wind, the objective value decreases since both UAVs are longer in midair and UAV-1 has to reduce its altitude to avoid wind zone 1. Regarding the initial fuel, UAV-1 needs 6.7% less (337.22kg instead of 361.42kg), although its trajectory is 32.3% longer (558.1km instead of 421.84km). This fuel saving is the result of the presence of advantegeous wind zones. On the other side, UAV-2 needs 15.14% more initial fuel (35.13kg instead of 30.51kg) for its 12.69% longer trajectory (435.03km instead of 386.04km). Since UAV-2 cannot benefit from tailwinds like UAV-1, it needs more fuel for the also longer trajectory.

Altitude-dependent range
As shown in Fig. 4, the altitude-dependent approach (2.3') in Sect. 2.3 results in a significant enlarged operating range compared to the constant approach (2.3) with u = 185 km.
To illustrate the new approach and its limits, we consider both UAVs, two waypoints, 200 time steps without any fine steps each, and t = 0.025h. The UAVs are stationed Fig. 10 Optimal flight trajectories for the altitude-dependent operating range approach for UAV-1 (orange) and UAV-2 (blue). Waypoints and ground controls are indicated by a red or green square, respectively. The dashed lines mark the altitude-dependent operating range of the respective UAV and the colored area its maximum operational distance to the waypoints at the location of their respective ground control, which is also their start and end location. Furthermore, we set ϕ alt,1 u = 0.05 and ϕ alt,2 u = 0.5, u ∈ U, to have a better resolution of the increasing maximum flight altitudes of the UAVs. All components of this instance are lined up along the x-axis and both waypoints can be visited at any time step. It took 3,908 seconds to solve this instance to proven global optimality. The optimal flight trajectories of the UAVs are displayed in Fig. 10.
As shown in Fig. 10, UAV-1 could reach the position of both points but it cannot visit WP-2 since for the necessary range it requires a flight altitude greater than the maximum operational distance to WP-2. So for visiting all points, UAV-2 has to be deployed, although UAV-1 is near the waypoint. Furthermore, the optimal trajectories can be divided into several phases, marked in Fig. 10: (a) First, both UAVs ascend after takeoff with their maximal ascend rate to reach their mass-dependent maximal altitude. (b) Due to the mass reduction by fuel consumption, the maximum altitude of every UAV increases over time so the UAVs keep ascending slowly. This behavior is controlled by the parameters ϕ alt,1 u and ϕ alt,2 u . (c) Near the respective waypoint, the UAVs slow down and descent to get into the operational range of the waypoint. (d) At their point of return, the UAVs achieve the maximum operational distance of the respective waypoint, visit it and start to fly back to their end location since no UAV can get into the operational distance to both points without leaving its operating range. To save fuel, they ascend again to their maximum altitude. (e) On their way back, the UAVs can reach a higher altitude compared to the beginning since they consumed most of their initial amount of fuel. (f) Every UAV descends to its end location with maximum descent rate to stay in the higher altitude bands as long as possible and benefits from its lower fuel consumption. This approach, to reach the top of descent point before going down, is the most economic descent strategy [25].

Mass-dependent flight trajectories
To illustrate the effect of mass-dependent flight dynamics to the computed trajectories, we consider an instance with 10 waypoints, two UAVs, 25 time steps with four fine time steps each, and two restricted airspaces and solve it in the presence and absence of mass-dependent flight dynamics. In terms of mass-dependent fuel consumption, we assume the additional fuel per mass factor ϕ can be set to an arbitrary value since it is no effect for the choice ϕ alt,2 u = 1. The total amount of CPU time to solve the instance neglecting the mass dependencies was 1,105 seconds, while in its presence 4,184 seconds were necessary. Figure 11 displays the optimal flight trajectories, altitude and velocity profiles, the fuel consumption, and the time windows of all waypoints for both instances.
As the results in Fig. 11 indicates, the optimal solution changes in the presence of mass-dependent flight dynamics. Both UAVs have shorter trajectories and return to their end location earlier, while only 8 of 10 waypoints were visited, one less than in the absence of mass dependencies.
While, without the influence of mass, both UAVs use their fastest throttle band and attain their respective maximum velocity, in the presence of the reduction factor ϕ acc u , they cannot reach it and must stay longer in the second-fastest throttle band. Thus, UAV-1 cannot reach waypoint 9 and UAV-2 cannot reach waypoint 3 within the given time.
The resulting shorter trajectories are beneficial to the attained altitude of UAV-1. It can now ascend higher since it visits fewer waypoints. But it can reach its maximum altitude near the end of its flight due to the mass-dependent maximum altitude. For UAV-2, one can see the influence of the abovementioned restriction on its flight dynamics. Whenever it ascends after visiting a waypoint, it can reach a higher altitude since its fuel mass reduces over time.
In terms of the initial amount of fuel, although UAV-1 has a 27.9% shorter trajectory (382.7km instead of 530.82km), it requires 6.54% more initial fuel (from 365.03kg to 388.89kg). In contrast, UAV-2 has an 9.03% lower amount of initial fuel (from 38.89kg to 35.38kg) for an also 37.21% shorter trajectory (282.46km instead of 449.84km). This contrary behavior of UAV-1 is justified in the fuel per mass factor ϕ f uel u since it applies as a percentage of the fuel mass. Thus, larger UAVs with higher fuel consumption are affected more than smaller ones and in comparison, UAV-1 can hold nearly four times the amount of fuel of UAV-2.

Collision avoidance
In this section, we take a closer look at the computation time of RCPs for an increasing number of UAVs and outline the solution quality of the optimal solutions compared to  For the two-dimensional case, we consider a circle with center C = (140, 140) and radius r RC P = 90 km and use UAVs from the type UAV-1. The safety distances in each coordinate direction are ε i = 9.26 km, i ∈ {x, y}. They are chosen accordingly to the minima for radar seperation, published by the International Civil Aviation Organization (ICAO) [14]. For t = 0.02, we get T min = 180 232·0.02 = 39, but we initialize the calculation with T = 40 since the initial maximum velocity is reduced by the present mass dependencies.
Varying the number of participating UAVs between 2 and 15, we generate 41 instances and solve them with the sequential method, see page 22, with a time limit of 3600 seconds per iteration. The results are displayed in Fig. 12.
For small problems with up to five UAVs, at least 98% of all instances were solved with a feasible solution found within the given time limit. In medium-scaled problems with six to 13 participating UAVs, this percentage decreases from 88% to 22%, while it drops below 7% for large problems considering 14 or more UAVs.
To illustrate the result of Proposition 4, we generate the upper bound for the necessary number of time steps for the 356 solved instances and compare its minimum and median value against the computed solutions for up to 14 participating UAVs. In Fig. 13, one can see that already for small instances with four UAVs, there is a significant gap between the given upper bound and the found optimal solution in both median and minimum value. Furthermore, this gap increases rapidly for larger numbers of UAVs.
In Fig. 14, the optimal trajectories of a RCP instance with 14 UAVs are displayed for two different time steps, also including the circle used to calculate the upper bound in Proposition 4. UAVs 2 and 11 are the first pair violating the given safety distances.  Since they start near to each other, the generated circle for the upper bound has radius r ub = 86.12 km, which is only about 4% smaller than the original radius r RC P .
For the three-dimensional RCP, we take the same parameters as in the twodimensional case, together with the initial altitude r z u (0) = 10 km, vertical safety distance ε z = 1 km, and altitude bounds h = 9 km and h = 10.972 km for all UAVs. Thus, they are all within their second-highest altitude band. The number of participating UAVs is varied between 2 and 20, generating 41 instances for each case. For every instance, there is a time limit of 3600 seconds per iteration. Again the sequential Fig. 15 Solution time of the three-dimensional RCP depending on the number of participating UAVs. The count of solved instances for every number of UAVs is displayed above its respective column method, described at page 22, is used to solve all instances. The results are displayed in Fig. 15.
Allowing the UAVs to avoid conflicts by choosing different altitudes, the solution process is accelerated significantly. For 676 problems the sequential method found a feasible solution within the given time limit and for instances up to 16 participating UAVs, 90% of all cases were solved. Considering at most 13 UAVs, the desired first feasible solution was found within 252 seconds for all problems. For large-scale problems with at least 17 participating UAVs, the number of solved rapidly decreases to 10% at the maximum amount of 20 UAVs.

Conclusions and future work
In this work, we considered the mission and flight planning problem for an inhomogeneous fleet of UAVs and extended an existing model from [10] to achieve solutions with a higher level of detail. A two-time-level discretization was applied to smooth the resulting trajectory without increasing the number of discrete decision variables. The maximum velocity, acceleration, reachable altitude, and the fuel consumption of the UAVs were related to their changing weight during the flight. Taking the curvature of the earth into account, the maximum operating range of the UAVs is now depending on its current altitude. Occurring restricted air spaces are now described by arbitrary polyhedrons and every UAV is affected by time-dependent wind if it is in a polyhedral wind zone. Furthermore, the model was adapted to fit the RCP for collision avoidance in two and three dimensions and a heuristic solution for this problem class in both dimensions was derived.
By applying different linearization techniques, the model got applicable to the MILP solver GUORBI and in numerical tests, several of its aspects were examined. Apply-ing the two-time-level discretization leads to smoother trajectories for an increasing number of fine time steps, but at the cost of a higher computation time. The experiments showed, that the differences between the computed trajectories get smaller for many fine time steps, thus one should find a trade-off between the level of detail and the necessary computation time. A drawback of this approach is the addition of further discretization errors since the binary variables can no longer change in every time step. Further experiments were performed to examine the influence of the model enhancements and it turned out that the fuel-dependent flight dynamics as well as the new introduced wind zones have significant impact on the optimal solution and the computation time. In a last experiment, randomly generated instances for the RCP were solved by a sequential method and the solutions were compared to these calculated by the derived heuristic, showing a large gap especially when a large number of participating UAVs is considered.
With this article, we provide an optimization approach linking the routing of an inhomogeneous fleet of UAVs with detailed trajectory calculation and environmental constraints, a main problem identified in [7]. Possible applications are surveillance or observation tasks, like tracking icebergs [1] or monitoring emission limits [31] in shipping and observing damages after a natural disaster [8]. Our future work will focus on the acceleration of the solution process to solve also larger instances within a reasonable time, the further extension of the model to cover more applications, and the reduction of discretization errors like corner cutting.
Funding Open Access funding enabled and organized by Projekt DEAL.
Data availability All models and data used to compute the presented results are available in the Zenodo repository, https://doi.org/10.5281/zenodo.6554723. The detailed solutions of the least squares problems (7) and (10) are available from the corresponding author at request.

Conflict of interest
The authors declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.