# A Stochastic Switched Optimal Control Approach to Formation Mission Design
for Commercial Aircraft

María Cerezo-Magaña, Alberto Olivares, Ernesto Staffetti Universidad Rey Juan
Carlos, Department of Telecommunication Engineering, Fuenlabrada, Madrid,
Spain (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>ernesto.staffetti@urjc.es).

###### Abstract

This paper studies the formation mission design problem for commercial
aircraft in the presence of uncertainties. Specifically, it considers
uncertainties in the departure times of the aircraft and in the fuel burn
savings for the trailing aircraft. Given several commercial flights, the
problem consists in arranging them in formation or solo flights and finding
the trajectories that minimize the expected value of the direct operating cost
of the flights. The formation mission design problem is formulated as an
optimal control problem of a stochastic switched dynamical system and solved
using nonintrusive generalized polynomial chaos based stochastic collocation.
The stochastic collocation method converts the stochastic switched optimal
control problem into an augmented deterministic switched optimal control
problem. With this approach, a small number of sample points of the random
parameters are used to jointly solve particular instances of the switched
optimal control problem. The obtained solutions are then expressed as
orthogonal polynomial expansions in terms of the random parameters using these
sample points. This technique allows statistical and global sensitivity
analysis of the stochastic solutions to be conducted at a low computational
cost. The aim of this study is to establish if, in the presence of
uncertainties, a formation mission is beneficial with respect to solo flight
in terms of the expected value of the direct operating costs. Several
numerical experiments have been conducted in which uncertainties on the
departure times and on the fuel saving during formation flight have been
considered. The obtained results demonstrate that benefits can be achieved
even in the presence of these uncertainties.

###### Index Terms:

Formation Flight, Formation Mission Design, Commercial Aircraft, Stochastic
Switched Systems, Stochastic Optimal Control, Generalized Polynomial Chaos.

##  I Introduction

This paper studies the formation mission design problem for commercial
aircraft in the presence of uncertainties, considering extended formations, in
which the longitudinal distance between aircraft is between 10 and 40
wingspans. This paper is actually a follow-up to [1], where the potential
contribution of formation flight to mitigating the environmental impact of
aviation and increasing the capacity of the Air Traffic Management (ATM)
system is discussed and the formation mission design problem is studied in the
absence of uncertainties using deterministic switched optimal control
techniques.

The solution of the formation mission design problem in the presence of
uncertainties is stochastic, i.e, its components are random processes, which
can be characterized by their mean and standard deviation functions. The
expected values and the standard deviations of the latitude and longitude of
the optimal trajectories with respect to the expected value of the timing and
the expected values and standard deviations of the timing of the trajectory as
functions of the expected value of the distance are useful for ATM purposes,
such as conflict detection and traffic synchronization. The expected values
and standard deviations of the fuel consumption and final time of each flight
of the formation mission permit the Direct Operating Costs (DOC) to be
estimated in such a way that airlines can establish to what extent a formation
mission is economically beneficial.

The relative contributions of each random variable to the variability of the
latitude and longitude of the optimal trajectories of each aircraft as
functions of the expected value of the timing and the variability of the
timing of the trajectories of each aircraft as functions of expected value of
the distance allow both the ATM staff to establish which sources of
uncertainty must be reduced to increase the predictability of the trajectories
and airlines to determine what sources of uncertainty must be reduced to
decrease the DOC.

The proposed methodology is an effective tool for solving the formation
mission design problem in the presence of uncertainties and quantifying the
effects of uncertainties in the departure times of the aircraft and in the
fuel burn savings of the trailing aircraft on the solution of the problem.

Formation flight can have a great added value in two major current air
transport-related concerns, the decarbonisation of the aviation sector and an
increase in the (ATM) capacity [2].  The key enabling factors for this concept
of operation are a formation control system to keep the aircraft flying in
formation at the optimal relative position to optimize the fuel burn savings
for the trailing aircraft and an ATM system capable of synchronizing flights
departing from different airports to ensure that the formation mission occurs
as planned. In extended formations, due to the great distance between the
leader and the follower aircraft, instabilities, such as meandering and
external factors, may affect the motion of the vortices. It is therefore
important to maintain the relative position between the follower aircraft and
the leader’s wake vortices precisely, as the fuel burn savings are very
sensitive to that relative positioning. This can be done using a formation
control system capable of continuously locating the wake vortices, maintaning
the aircraft in the optimal relative position, and, in this way, optimizing
the fuel savings during the formation flight [3].

While technical issues related to maintain an efficient and safe formation
flight have been solved in the last years, some concerns have yet to be
addressed. In particular, a major change in airworthiness standards, policies,
and procedures is required. In the 40th International Civil Aviation
Organization (ICAO) Assembly, formation flight was proposed as a strategic
objective, and the need to develop a new operational concept which includes
reduced separations between aircraft allowing formation was established [4].

Additionally, the main aircraft manufacturers are making a serious effort to
implement formation flight. Airbus, the major European aircraft manufacturer,
is undertaking an ambitious project called Fello’fly [5], which aims to
demonstrate the technical, operational, and economical viability of formation
flight for long-haul commercial flights. In November 2021, the first flight
test was carried out in transoceanic
flights111https://simpleflying.com/airbus-a350s-bird-like-flight/. In this
project, Airbus is collaborating with airlines and air navigation service
providers in order to tackle the formation flight challenges. The objective of
this project is not only to demonstrate the operational feasibility of
formation flight, but also identify safety procedures and standards for
transatlantic operations, enabling a controlled entry into service by 2025.
Boeing has conducted an extended formation flight study, in which flight tests
on Boeing’s ecoDemonstrator platform have been carried out [6].

It is well known that uncertainty in flight departure times is among the main
causes of trajectory uncertainty, which generates inefficiency in the ATM
system [7]. However, timing is a crucial factor in formation missions. Indeed,
considering the usual cruise speed of most long-haul commercial aircraft,
missing the rendezvous location by, for instance, ten minutes means spatially
missing the partner aircraft by 150 km. Such cases require catch-up maneuvers,
which result in a loss of performance compared to the planned formation
mission. Thus, uncertainties in both departure times and fuel burn savings for
the trailing aircraft must be considered in the formulation of the formation
mission design problem and their effects on the resulting aircraft
trajectories and the DOC of the formation mission must be quantified. Although
aircraft trajectories are not only affected by uncertainty in the initial
conditions but also by uncertainties in the aircraft performance models,
operational uncertainty, and weather uncertainty, mainly wind and storms, this
paper only considers uncertainties in the departure times of the aircraft and
in the fuel burn savings for the trailing aircraft. These uncertain parameters
are modeled by means of random variables, which are described by probability
distribution functions. The aircraft model is assumed to be known with
precision as is the relevant wind field.

Given several commercial flights, the problem consists in establishing how to
organize them in formation or solo flights and in finding the trajectories
that minimize the expected value of the DOC of the formation mission. This
paper only considers formations of up to three aircraft, because the
difficulties in synchronizing more than three flights make formations of more
than three aircraft operationally impractical. Moreover, increasing the number
of aircraft in the formation asymptotically reduces the benefits obtained from
the formation. Since each aircraft can fly solo or in various positions within
a formation, the mission is modeled as a stochastic switched dynamical system,
in which the flight modes of the aircraft are described by sets of stochastic
ordinary differential equations, the discrete state describes the combination
of flight modes of the individual aircraft, and logical constraints establish
the switching logic among the discrete states of the system. In this paper,
the formation mission design problem is formulated as an optimal control
problem of a stochastic switched dynamical system.

A deterministic switched dynamical system is a particular type of hybrid
system that consists of several subsystems and a switching law that specifies
the active subsystem at each time instant. Deterministic switched dynamical
systems are described by both a continuous and a discrete dynamics, in which
the transitions among discrete states are not established in advance. In
particular, in this paper, each aircraft is assumed to have different flight
modes, namely solo flight and flight in different positions within a
formation, and their combination is represented by the discrete state of the
switched dynamical system, which models their joint dynamic behavior. Each
flight mode is represented by different dynamical equations, which may include
or not include formation flight benefits in terms of fuel burn savings.
Additionally, logical constraints in disjunctive form, based on the stream-
wise distance between aircraft, model the switching logic among the discrete
states of the system. Stochastic switched dynamical systems inherit all the
features of the deterministic switched dynamical systems. Additionally, both
the discrete and the continuous dynamics are affected by uncertainties. This
general definition can encompass a wide range of stochastic phenomena. This
paper considers those stochastic hybrid systems in which random variables only
affect the continuous dynamics. More specifically, random variables represent
uncertain parameter values and uncertain initial conditions. The probability
distribution functions of these random variables are assumed to be known. In
this case, the continuous dynamics is described by a set of stochastic
differential equations but the discrete dynamics is deterministic. The
adjective stochastic as used in this paper means that the solutions of the
differential equations depend on a vector of random variables. A systematic
account of recent developments regarding deterministic and stochastic hybrid
systems is given in the monography [8].

An optimal control problem of the stochastic switched dynamical system
described above is an optimal control problem in which the continuous dynamics
of the system is represented by stochastic differential equations, the
objective functional is a stochastic functional, and the constraints, which
are defined by means of stochastic functions, must be satisfied almost surely,
i.e., with probability 1. The set of possible exceptions in which the
constraints are not satisfied may be non-empty but must have probability 0.
The adjective stochastic for these functions and the objective functional
means in this paper that they depend on a vector of random variables. This
problem is referred to as the Stochastic Switched Optimal Control Problem
(SSOCP).

Modeling, control, and optimal control of stochastic switched systems are
addressed in [8], [9], in which the random factors acting on the continuous
dynamic models of the switched system in each discrete state are represented
by some idealized processes, such as the Wiener process, and tools such as
stochastic calculus have been employed to obtain solutions.

Another approach to model uncertainties in switched dynamical systems is to
treat uncertainties as random variables or random processes and recast the
original deterministic switched dynamical system as a stochastic switched
dynamical system. This type of stochastic systems are different from those
represented by classical stochastic differential equations, where the random
inputs are idealized processes. Consider a stochastic optimal control problem
in which only random variables are present in its formulation like in the
formation mission design problem studied in this paper. One of the most
commonly used methods to solve optimal control problems of stochastic
dynamical systems of this type is the Monte Carlo sampling, in which
independent realizations of the random variables are generated based on their
probability distributions. For each realization, the optimal control problem
becomes deterministic. The need for large number of solutions for accurate
results can lead to an excessive computational cost, especially for optimal
control problems that are already computationally intensive in the
deterministic settings, as the mission design problem considered in this
paper. The Generalized Polynomial Chaos (gPC) expansion can contribute to
alleviating this drawback. A systematic and coherent presentation of numerical
strategies for uncertainty quantification and stochastic computing is given in
[10], with a focus on the methods based on the gPC approach.

In stochastic collocation methods the stochastic model equations are satisfied
at a discrete set of points, called nodes, in the corresponding random space.
Polynomial approximation theory is used to locate the nodes strategically to
increase the numerical accuracy. With this approach, a small number of sample
points of the random variables are used to solve jointly particular instances
of the SSOCP. The obtained solutions are then expressed as orthogonal
polynomial expansions in terms of the random variables using these sample
points. This is a nonintrusive methodology because the model equations are not
altered. This technique allows statistical and sensitivity analysis of the
stochastic solutions to be conducted at a low computational cost. Thus, the
gPC method converts the SSOCP into an augmented deterministic Switched Optimal
Control Problem (SOCP), in which particular instances of the SSOCP are solved
together as a single deterministic optimal control problem. In this paper, the
resulting SOCP is solved using the method described in [1]. It is important to
point out that the gPC method is applicable to solve the SSOCP when the
solution depends smoothly on the random variables. This means that the
solutions obtained for each combination of sample points of the random
variables must give rise to the same discrete solution, i.e., to the same
sequence of discrete states of the switched dynamical system that represents
the formation mission.

Aircraft trajectories are specified by a sequence of geographical coordinates
of spatial positions and the timing, i.e., a sequence of time instants at
which the corresponding point of the trajectory must be reached. Thus, any
delay of the aircraft in reaching a spatial position is considered a temporal
deviation from the trajectory as much as a spatial deviation at a given time
instant. The uncertain position of an aircraft at a given time can be
described by a region of confidence. Likewise, the uncertain timing of an
aircraft trajectory at a given spatial position can be described by an
interval of confidence. In this paper, the expected values of the geographical
coordinates together with the expected value of timing obtained in the
solution of the SSOCP are considered as the reference trajectories to be
followed by the aircraft of the formation mission.

As previously mentioned, the solution of the formation mission design problem
in the presence of uncertainties is stochastic, i.e., its components are
random processes which can be characterized by their mean and standard
deviation functions. Given the gPC expansion, these functions can be directly
computed from the coefficients of this expansion [10]. The expected value
coincides with the first coefficient, whereas the variance is the sum of the
squares of the other coefficients of the gPC expansion. The statistical
information estimated from the stochastic solutions includes the expected
values and standard deviations of the latitude and longitude of the optimal
trajectories and the other state variables with respect to time, together with
the expected values and standard deviations of the timing of the trajectory as
functions of the distance. They also include the expected values and the
standard deviations of the arrival times and, in the case of formation flight,
the expected values and the standard deviations of the latitude, longitude,
and time of the rendezvous and splitting locations. The expected values and
standard deviations of the fuel consumption of each flight of the formation
mission are also estimated. This information is combined to estimate the
expected value of the DOC.

Additionally, in the numerical experiment that involves two random variables,
a global sensitivity analysis of the stochastic solution is also conducted
[11]. In this paper, the global sensitivity analysis is based on the Sobol’
indices, which enable the determination of what proportions of the variance of
the solutions of the switched optimal control problem can be attributed to the
different random variables of the switched optimal control model. The purpose
of the sensitivity analysis is to identify the random variables that have more
influence on the variability of a component of the stochastic solution with
the aim of reducing its variability acting on the source.

###  I-A Previous approaches

In [12], a review of the literature on stochastic modeling with applications
to ATM is provided, including literature on stochastic optimal control.

In [13], the problem of aircraft conflict prediction is studied for two-
aircraft midair encounters. First, a model is presented for prediction of the
aircraft positions along a time horizon, during which each aircraft is
following a prescribed flight plan in the presence of additive wind
perturbations on its velocity. Then, a method for estimating the probability
of conflict is proposed. This method is based on a Markov chain approximation
of the stochastic processes that models the aircraft flight. In [14], the
aircraft conflict prediction problem is formulated as a reachability problem
in a stochastic hybrid system framework. Specifically, a switching diffusion
model is employed to predict the future positions of an aircraft following a
given flight plan, and the probability that the aircraft enters an unsafe
region of the airspace is estimated using a numerical algorithm for
reachability computation.

In [15], an approach to aircraft trajectory optimization in the presence of
uncertainties based on gPC is presented, where only one random variable has
been considered, which represents an uncertain aerodynamic parameter of the
dynamic model of the aircraft. This random variable has been modeled as a
uniformly distributed random variable. In [16], a stochastic optimal control
method based on gPC is developed for determining conflict-free aircraft
trajectories under wind uncertainty. The random processes that represent the
components of the wind speed are approximated as a linear combination of
deterministic functions multiplied by independent random variables using the
Karhunen-Loève expansion.

In [1], a deterministic switched optimal control method is employed to solve
the formation mission design problem. An additional analysis has been
conducted to determine how departure delays and fuel burn savings affect the
formation flight, in terms of routes, rendezvous and splitting locations and
times, and flight times. In [17], the impact of ground delays on formation
flight is studied using stochastic dynamic programming.

###  I-B Contributions of the paper

This paper proposes a methodology for the solution of the formation mission
design problem in the presence of uncertainties in some of the parameters and
boundary conditions of the problem. These uncertainties are represented by
random variables characterized by probability density functions. The formation
mission design problem in the presence of uncertainties is formulated as an
optimal control problem of a stochastic switched system. This problem is
solved using an approach based on gPC, in which the stochastic switched
optimal control problem is transformed into an equivalent deterministic
switched optimal control problems in a higher dimensional state space and
solved using a numerical method developed by the same authors [1]. The
formulation of the problem based on switched optimal control allows accurate
dynamic models of the aircraft and meteorological forecast to be included in
the problem formulation. This technique permits statistical and global
sensitivity analysis of the stochastic solution to be conducted at a low
computational cost.

###  I-C Organization of the paper

The paper is organized as follows. The model of the switched dynamical system
that represents the formation mission is introduced in Sec. II. The
deterministic and the stochastic switched optimal control problems are stated
in sections III and IV, respectively. The generalized polynomial chaos
expansion technique and the method to solve the stochastic switched optimal
control problem are described in Sec.V. The results of the numerical
experiments are reported and analyzed in Sec. VI. The outcomes of the
sensitivity analysis of the solutions are discussed and interpreted in
Sec.VII, Finally, conclusions are drawn in Sec. VIII.

##  II Model of the system

This section outlines the model of the switched dynamical system that
represents the formation mission. Further details on the aircraft equations of
motions, flight envelope, and wind model employed can be found in [1, Sect.
II].

The mission design problem is studied only in the cruise phase of the flight.
Thus, a simplified two-degrees-of-freedom point variable-mass dynamic model is
considered assuming that all aircraft are in the cruise phase. The motion is
restricted to the horizontal plane at cruise altitude over a spherical Earth
model. A symmetric flight without sideslip is considered and all the aircraft
forces are supposed to be in the plane of symmetry of the aircraft. Wind
effects are also considered. All aspects associated with the rotational
dynamics are neglected. The set of kinematic and dynamic differential-
algebraic equations (DAE) that describe the motion of each aircraft of the
formation mission are

$\displaystyle\dot{\phi}(t)$ $\displaystyle=$
$\displaystyle\dfrac{V(t)\cdot\cos\chi(t)+{V_{W_{N}}}(t)}{R_{E}+h},$
$\displaystyle\dot{\lambda}(t)$ $\displaystyle=$
$\displaystyle\dfrac{V(t)\cdot\sin\chi(t)+{V_{W_{E}}}(t)}{\cos\phi(t)\cdot\left(R_{E}+h\right)},$
(1) $\displaystyle\dot{\chi}(t)$ $\displaystyle=$
$\displaystyle\dfrac{L(t)\cdot\sin\mu(t)}{V(t)\cdot m(t)},$
$\displaystyle\dot{V}(t)$ $\displaystyle=$
$\displaystyle\dfrac{T(t)-D(t)}{m(t)},$

where the state vector has five components: the two dimensional position
variables, latitude and longitude, denoted by $\phi$ and $\lambda$,
respectively, the heading angle $\chi$, the true airspeed $V$, and the mass of
the aircraft $m$. In this set of equations, the control vector has three
components: the thrust force $T$, the lift coefficient $C_{L}$, and the bank
angle $\mu$. The normalized version of the DAE system (II) is used in this
paper. $L=qSC_{L}$ is the lift force, where $q=\frac{1}{2}\rho V^{2}$ is the
dynamic pressure, $\rho$ is the air density, and $S$ is the reference wing
surface area. The aerodynamic drag force is $D=qSC_{D}$, where $C_{D}$ is the
drag coefficient. $h$ is the cruise altitude and ${V_{W_{E}}}$ and
${V_{W_{N}}}$ are the components of the wind velocity vector in eastward and
northward directions, respectively. The mass flow rate equation is

$\dot{m}(t)=-T(t)\cdot\eta(t),$ (2)

where $\eta$ is the thrust specific fuel consumption. The Eurocontrol’ s Base
of Aircraft Data (BADA), version 3.6 [18], has been used to determine this.
Thus, for aircraft $p$, the state vector is
$x_{p}(t)=\left(\phi_{p}(t),\lambda_{p}(t),\chi_{p}(t),V_{p}(t),\right.$
$\left.m_{p}(t)\right),\forall p\in\left\\{1,\ldots,N_{a}\right\\}$, and the
control vector is $u_{p}(t)=\left(T_{p}(t),C_{L_{p}}(t),\right.$
$\left.\mu_{p}(t)\right),\forall p\in\left\\{1,\ldots,N_{a}\right\\}$, with
$N_{a}$ the number of aircraft involved in the mission design problem.

Flight envelope constraints model aircraft performance limitations. These
constraints are algebraic constraints usually expressed as simple bound
constraints that involve, for instance, flight altitude, load factor, and
airspeed. BADA 3.6 is used to model the upper and lower bounds of the flight
envelope constraints. These constraints are included in the set of path
constraints of optimal control problems.

Wind information from ERA-Interim [19], provided by the European Centre for
Medium-Range Weather Forecasts (ECMWF) is used to compute the ${V_{W_{E}}}$
and ${V_{W_{N}}}$ components of the wind speed in Eqs. (II). Since this
information is given as grid data, analytic functions that approximate the
grid data must be determined to include wind information in the stochastic
switched optimal control problem used to solve the mission design problem.
Radial basis functions [20] are employed for this purpose. Wind data from
April 30, 2019, at 12:00 are used.

The DAE system consisting of Eqs. (II) and Eq. (2) is applicable to any
aircraft in the cruise phase of solo flights. In formation flights, it is
necessary to include a modification in the mass flow rate equation (2) for the
trailing aircraft. It is assumed that the reduction in the induced drag of the
trailing and intermediate aircraft can be modeled as a percentage of the
reduction in fuel consumption, $\mathcal{R}_{\text{fuel}}$. In this way, the
fuel saving factor can be directly incorporated into Eq. (2) to obtain the
mass flow rate equation for the trailing aircraft during formation flight

$\dot{m}(t)=-[1-\mathcal{R}_{\text{fuel}}]\cdot T(t)\cdot\eta(t).$ (3)

In a two-aircraft formation, this reduction is applied to the trailing
aircraft when the distance to the leader aircraft is shorter than 20
wingspans. In a three-aircraft formation, this reduction is applied to the
trailing and intermediate aircraft when the distances between the leader and
intermediate aircraft and between the intermediate and the trailing aircraft
are shorter than 20 wingspans. Notice that the distance of 20 wingspans is
used for checking if aircraft achieve benefits from the formation flight and
not for providing flight trajectories strictly accounting for this distance.

##  III The Deterministic Switched Optimal Control Problem

Since each aircraft can fly solo or in different positions within a formation,
the mission is modeled as a switched dynamical system, in which the discrete
state describes the combination of flight modes of the individual aircraft.

In this section, following [1, Sect. III], the formulation of the SOCP is
introduced for a two-switched dynamical system. The dynamical model of a two-
switched dynamical system is

$\displaystyle\dot{x}_{S}(t)$ $\displaystyle=$ $\displaystyle
f_{v_{S}(t)}(t,x_{S}(t),u_{S}(t)),$ $\displaystyle x_{S}(t_{I})$
$\displaystyle=$ $\displaystyle x_{I}\in\mathbb{R}^{n},$ (4) $\displaystyle
x_{S}(t_{F})$ $\displaystyle=$ $\displaystyle x_{F}\in\mathbb{R}^{n},$
$\displaystyle v_{S}(t)$ $\displaystyle\in$
$\displaystyle\\{0,1\\},\;t_{I}\leq t\leq t_{F},$

where the continuously differentiable vector fields,
$f_{0},f_{1}:\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$,
describe the dynamics of the two possible modes of the system. The control
input $u_{S}(t)\in\Omega\subset\mathbb{R}^{m}$ is constrained to belong, at
each time instant, to the bounded and convex set $\Omega$. The binary variable
$v_{S}(t)$ is the mode selection variable that specifies which of the two
possible system modes, $f_{0}$ or $f_{1}$, is active. Thus, both $u_{S}(t)$
and $v_{S}(t)$ can be regarded as control variables. The initial time $t_{I}$,
the final time $t_{F}$, the initial state $x_{S}(t_{I})$, and the final state
$x_{S}(t_{F})$ are assumed to be restricted to a boundary set $B$, namely,
$(t_{I},x_{S}(t_{I}),t_{F},x_{S}(t_{F}))\in B=T_{I}\times B_{I}\times
T_{F}\times B_{F}\subset\mathbb{R}^{2n+2}$. The objective functional of the
SOCP is

$\displaystyle J_{S}\left(t,x_{S}(t),u_{S}(t),v_{S}(t)\right)=$ (5)
$\displaystyle=$ $\displaystyle
g(t_{F},x_{F})+\int_{t_{I}}^{t_{F}}F_{v_{S}(t)}(t,x_{S}(t),u_{S}(t))dt,$ (6)

where $g$ is the endpoint cost function defined on a neighborhood of $B$ and
$F_{0}$ and $F_{1}$ are real-valued continuously differentiable functions that
represent the running cost of the system in each mode. In general, other
constraints are included in the formulation of the problem, such as the path
constraints, which apply over the whole path or at intermediate points and not
only at the end points, and the logical constraints, which establish the
switching logic among the discrete states of the system.

The SOCP is stated as follows

$\min_{u_{S}\in\Omega,v_{S}\in\\{0,1\\}}J_{S}\left(t,x_{S}(t),u_{S}(t),v_{S}(t)\right),$
(7)

subject to Eq. (III), endpoint constraints
$(t_{I},x_{S}(t_{I}),t_{F},x_{S}(t_{F}))\in B$, path constraints, and logical
constraints.

In the formation mission design problem studied in this paper, not all the
dynamic equations are subject to switches. More specifically, the equations of
motion associated with the state variables ${\phi}$, ${\lambda}$, ${\chi}$,
and ${V}$ do not switch, whereas the aircraft’ s mass flow rate equation does
when the aircraft joins or leaves a formation as a trailing or an intermediate
aircraft. When an aircraft is flying solo or as a leader in the formation it
obtains no benefits from the formation flight in terms of fuel savings and Eq.
(2) describes its mass flow. On the contrary, when an aircraft is flying as
the intermediate or trailing aircraft in the formation it obtains benefits
from the formation flight and its mass flow is described by Eq. (3). Thus,
$f_{0}$ represents Equations (II) and (2), whereas $f_{1}$ represents
Equations (II) and (3). The functions $F_{0}$ and $F_{1}$, which appear in the
running cost of Eq. (6), are not present in the formulation of the problem.
The objective functional $J_{S}$ of the switched optimal problem that models
the formation mission design problem is

$J_{S}=\alpha_{t}\sum_{p=1}^{N_{a}}t_{{flight}_{p}}+\alpha_{f}\sum_{p=1}^{N_{a}}{{m}_{f_{p}}},$
(8)

where $\alpha_{t}$ and $\alpha_{f}$ are the time and the fuel consumption
weighting parameters, respectively, $t_{{flight}_{p}},\forall
p\in\left\\{1,\ldots,N_{a}\right\\}$, is the total flight time for aircraft
$p$, and ${{m}_{f_{p}}},\forall p\in\left\\{1,\ldots,N_{a}\right\\}$, is the
fuel consumption for aircraft $p$.

The switching logic among the discrete states of the system is defined by the
logical constraints. They express the fact that not all transitions between
discrete states are possible at any time. For example, in a three-aircraft
formation mission, the transition between the discrete state in which the
three aircraft are flying solo and that in which they are flying in a three-
aircraft formation is not permitted. In this case, passing through an
intermediate discrete state in which two aircraft are flying in formation and
the third is still flying solo is mandatory. Any logic expression can be
transformed to conjunctive normal form, i.e., it can be expressed as a
conjunction (sequence of ANDs) of one or more clauses, where a clause is a
disjunction (OR) of literals, which represent statements or negations of
statements. To incorporate the logical constraints in the SOCP, they must be
converted into a set of equality or inequality constraints.

Therefore, the formation mission design problem is formulated as an optimal
control problem of a switched dynamical system with logical constraints and
solved using the techniques described in [1, Sect. IV] and [1, Sect. V], in
which binary variables are not used. The result is a smooth optimal control
problem without binary variables, which reduces the computational complexity
of finding a solution. This problem is solved using a knotting pseudospectral
method as described in [1, Sect. VI]. The solution of the deterministic
optimal control problem is denoted by $z(t)$. It includes the control and
state variables and their timing, from which the rendezvous, splitting, and
final times can be read off.

##  IV The Stochastic Switched Optimal Control Problem

Let $\theta=(\theta_{1},\theta_{2},\ldots,\theta_{M})$ be the vector of random
variables that represent the random parameters of the SOCP. The presence of
random variables in the SOCP converts it into an SSOCP. In the formation
mission design problem studied in this paper, not all the elements of the SOCP
contain random variables. Specifically, the fuel consumption reduction factor
$\mathcal{R}_{\text{fuel}}$ in the mass flow rate equation (3) of the aircraft
flying in formation as a trailing or an intermediate aircraft is a random
variable. In contrast, the equations of motion associated with the state
variables ${\phi}$, ${\lambda}$, ${\chi}$, and ${V}$ do not contain random
variables. The departure times of the aircraft are also assumed to be random
variables of the SOCP. One of these departure times coincides with $t_{I}$,
the initial time of the formation mission. These random variables are assumed
to be independent and characterized by probability density functions.

The solution of the SSOCP includes both the state and control variables of the
system that represents the formation mission. Due to the presence of the
vector random variable $\theta$ in the SSOCP, the solution is stochastic,
i.e., it depends on the actual realization of the vector random variable. In
particular, it is a random function of state or time. Therefore, the solution
is actually a random process, characterized by the mean function and the
corresponding 95% confidence envelope. Specifically, the observed values of a
state variable at a given time $t$ is a random variable. In this case, the
random process is a random function of time. Likewise, the time at which a
given value of the state variable is reached is a random variable. This means
that the timing of the trajectories, which is defined as the time at which a
state is reached, is also a random process. In this case, the random process
is a random function of the state. In this paper, the timing of the
trajectories is represented as a random function of the orthodromic distance
from the departure location. For the sake of ease of exposition, the
stochastic solution of the SSOCP is denoted in the following sections as a
random function of time. All the random variables observed at specific time
instants or states are characterized by their expected values and by the
corresponding 95% confidence intervals, including the arrival times and the
fuel consumptions of the aircraft of the formation mission.

The objective functional of the SSOCP is the expected value of the objective
functional (6) of the SOCP. In the formulation of the SSOCP, the dynamic
equations, the path constraints, and the boundary conditions formally depend
on the vector random variable $\theta$ and must be satisfied almost surely.

##  V The Generalized Polynomial Chaos Expansion

In this section, following [16], the gPC expansion is introduced and the
method for determining the stochastic solution of the SSOCP and computing its
statistical information is described.

Let $z(t,\theta)$ denote the stochastic solution of the SSOCP, which includes
$x(t,\theta)$ and $u(t,\theta)$. The P-th order gPC approximation of
$z(t,\theta)$ can be expressed as:

$z_{P}(t,\theta)=\sum_{m=1}^{M}C_{m}(t)\Phi_{m}(\theta),$ (9)

where $\theta=(\theta_{1},\theta_{2},\ldots,\theta_{N})$ is a vector of
independent random variables, $C_{m}(t)$ are the coefficients of the
expansion, which can be obtained using either a nonintrusive or an intrusive
approach, and $\Phi_{m}(\theta)$ are the multivariate orthogonal polynomial
basis functions, which are calculated from the $l_{i}$-th order one-
dimensional polynomial basis function $\phi^{(l_{i})}(\theta_{i})$ of the
random variable $\theta_{i}$ by means of the tensor product rule as follows:

$\Phi_{m}(\theta)=\prod_{i=1}^{N}\phi_{i}^{(l_{i})}(\theta_{i}).$ (10)

Notice that, in Eq. (10), a unique combination of $l_{i},i=1,2,\dots,N$,
corresponds to each subscript $m$, which satisfy the condition
$\sum_{i=1}^{N}l_{i}\leq P$, where $P$ is the maximum degree of the
multivariate polynomial $\Phi_{m}$. The number of tensor product basis
functions is $M=\binom{N+P}{N}$. The orthonormal polynomials in Eq. (10)
satisfy the following orthogonality condition

$E[\phi_{i}^{(j)}(\theta_{i})\phi_{i}^{(k)}(\theta_{i})]=\int\phi_{i}^{(j)}(\theta_{i})\phi_{i}^{(k)}(\theta_{i})\rho_{i}(\theta_{i})d\theta_{i}=\delta_{jk},$

where $E$ is the expected value operator, $\rho_{i}(\theta_{i})$ is the
probability density function of the random variable $\theta_{i}$, and
$\delta_{jk}$ is the Kronecker delta function.

To improve convergence, the choice of the orthonormal polynomials should be
made on the basis of the probability density function $\rho_{i}(\theta_{i})$.
For instance, the Legendre polynomials are the best choice for uniform random
variables, whereas Hermite polynomials are the best option for Gaussian random
variables. In this paper, the coefficients $C_{m}$ of the expansion (9) are
calculated using a nonintrusive gPC based stochastic collocation method as
follows

$C_{m}=E[z(t,\theta)\Phi(\theta)]=\int
z(t,\theta)\Phi_{m}(\theta)\rho(\theta)d\theta,$ (11)

where $\rho(\theta)=\prod_{i=1}^{N}\rho_{i}(\theta_{i})$ is the joint
probability density function of the vector random variable $\theta$. A
Gaussian quadrature can be used to approximate the integral in Eq. (11). The
univariate quadrature rule approximates the polynomial
$\phi_{i}^{(l_{i})}(\theta_{i})$, $i=1,2,\ldots,N$, using a set of $q$
collocation points $\theta_{i}^{(j)}$ and calculating the corresponding set of
quadrature weights $\alpha^{(j)}$, $j=1,2,\ldots,q$. Higher precision in the
quadrature can be achieved by increasing the number of collocation points $q$.
Using the tensor product rule, the total number of collocation points is
$q^{N}$. This number, which could become large when $N$, the number of random
variables in $\theta$, grows, can be reduced using the sparse grid quadrature
based on the Smolyak rule [10].

Let $\theta^{(j)}$ be the collocation points and $\alpha^{(j)}$,
$j=1,\ldots,Q$, the corresponding weights. Then, the coefficients of the
expansion $C_{m}$ in Eq. (11) are approximated by:

$C_{m}\approx\sum_{j=1}^{Q}z(t,\theta^{(j)})\Phi_{m}(\theta^{(j)})\alpha^{(j)},$
(12)

where $z(t,\theta^{(j)})$ denotes the solution of the augmented SOCP obtained
using the $j$-th collocation point $\theta^{(j)}$ and $\Phi_{m}(\theta^{(j)})$
represents the multivariate orthogonal polynomial basis function evaluated at
the $j$-th collocation point $\theta^{(j)}$. Thus, the SSOCP is converted into
an augmented deterministic SOCP, in which particular instances of the SSOCP,
which correspond to the collocation points of the random variables, are
combined into a single optimal control problem. The resulting augmented SOCP
is solved by means of the numerical method described in Sec. III. The expected
value and variance of the stochastic solution $z_{P}(t,\theta)$ are calculated
using the coefficients of the expansion as follows:

$\displaystyle\text{E}[z_{P}(t,\theta)]$ $\displaystyle=$ $\displaystyle
C_{1}(t),$ (13) $\displaystyle\text{VAR}[z_{P}(t,\theta)]$ $\displaystyle=$
$\displaystyle\sum_{m=2}^{M}C_{m}^{2}(t).$ (14)

A complete treatment of the gPC expansion can be found in the monography [10],
whereas the procedure followed to determine the stochastic solution of the
SSOCP and to compute its statistical information is schematically represented
in Figure 1. Further details can be found in [15] and [16].

Determine the set of collocation points $\theta^{(j)}$ and compute the
associated weights $\alpha^{(j)}$, $j=1,\ldots,Q$, for the vector random
variable $\theta$. Vector random variable $\theta$ Calculate the value of the
orthonormal polynomial $\Phi_{m}(\theta^{(j)})$ at each collocation point
$\theta^{(j)}$, $j=1,\ldots,Q,m=1,\ldots,M$. Solve the augmented SOCP for
computing the deterministic solutions $z(t,\theta^{(j)})$ of the SSOCP at each
collocation point $\theta^{(j)},j=1,\ldots,Q$, using the deterministic
numerical method described in Sec. III. Aircraft model, wind forecast, flight
data Determine the approximated values of the coefficients of the expansion
$C_{m}(t),m=1,...,M$, using Eq. (12). Calculate the approximate stochastic
solution $z_{P}(t,\theta)$ of the SSOCP using Eq. (9). Compute the mean and
variance functions of the approximate stochastic solution $z_{P}(t,\theta)$
using Eqs. (13) and (14). Figure 1:  Schematic diagram of the procedure for
determining the stochastic solution of the SSOCP and computing its statistical
information.

##  VI Numerical Results

In this section, the results of the following numerical experiments are
reported to show the effectiveness of the proposed methodology to solve the
stochastic formation mission design problem:

  * •

Experiment A: Three-aircraft transoceanic mission design with uncertainty in
the fuel burn savings.

  * •

Experiment B: Two-aircraft transoceanic mission design with uncertainty in
departure times of the flights.

All the experiments involve transoceanic eastbound flights. Wind data from the
ERA-Interim reanalysis database of the ECMWF are used. As mentioned above, in
this paper only the cruise phase is modeled; the rest of the flight phases are
neglected. Thus, the initial and final locations of the cruise phase of the
flights are assumed to be the latitudes and longitudes of the departure and
arrival airports of each flight at cruise altitude. Airbus A330-200 aircraft
BADA models are considered for each flight. The initial masses of the aircraft
are assigned, as are the initial and final velocities, which are set at
typical cruise values for the selected aircraft models. The initial heading
angles are set to the initial heading angles of the orthodromic paths between
the initial and final locations of each flight. The time and fuel burn
weighting parameters, $\alpha_{t}$ and $\alpha_{f}$, in the objective
functional (8) are set to 0.3 and 0.7, respectively, based on [21], an
Eurocontrol study on the delays-related costs.

The numerical experiments have been conducted on a 3.6 GHz Intel Core i9
computer with 32 GB RAM. The computational times reported in this section
include both the time required to generate the warm-start solution and the
time to find the optimal solution of the problem. The warm-start solution is a
feasible solution of the problem generated by the NLP solver from an initial
guess of the solution.

As previously mentioned, the gPC method converts the SSOCP into an augmented
deterministic SOCP, which is solved using a knotting pseudospectral method.
This method transforms the deterministic SOCP into an NLP problem, which is
modeled using Pyomo [22], a Python-based software package designed to model
complex optimization problems, and solved using the Interior Point OPTimizer
(IPOPT) solver. In this paper, the initial guesses for the latitude, the
longitude, and the heading angle are generated using the orthodromic paths
between the departure and arrival locations of each flight. Typical cruise
velocity and fuel consumption of the aircraft model selected are used to
generate the initial guesses of the velocity and the mass of each aircraft
during the flight, respectively.

###  VI-A Experiment A: Three-aircraft transoceanic mission design with
uncertainty in the fuel burn savings.

Experiment A involves three transoceanic eastbound flights, Flight 1, Flight
2, and Flight 3, with uncertainty in the fuel burn savings for the trailing
aircraft. The three flights considered in this experiment have the following
departure and arrival locations:

  * •

Flight 1: New York (JFK) - Paris (CDG).

  * •

Flight 2: Boston (BOS) - Madrid (MAD).

  * •

Flight 3: Montreal (YUL) - London (LHR).

The departure times of Flight 1, Flight 2, and Flight 3 are set to 10:15,
10:30, and 10:50, respectively. The boundary conditions for the state
variables of the three aircraft are given in Table I. Flight 1, Flight 2, and
Flight 3 are operated by Aircraft 1, Aircraft 2, and Aircraft 3, respectively.

The formation configuration is selected in advance. In the case of a two-
aircraft formation that includes Aircraft 2, Aircraft 2 is the leader and
Aircraft 1 or Aircraft 3 is the follower. In the case of a two-aircraft
formation that does not include Aircraft 2, Aircraft 3 is the leader and,
accordingly, Aircraft 1 is the follower. The only three-aircraft formation
allowed is the in-line formation, where Aircraft 2 is the leader, Aircraft 3
is the intermediate, and Aircraft 1 is the follower.

The random variable that represents the fuel burn savings for the trailing
aircraft is denoted by $\theta$. This random variable is modeled as a Gaussian
random variable with mean $0.1$ and standard deviation $0.02$, i.e.,
$\theta\sim\mathcal{N}(0.1,0.02)$, based on several aerodynamic models and
flight tests reported in [23]. The departure times are fixed.

Solving this problem entails deciding which mode of flight, i.e., formation or
solo flight, is optimal, the expected values and standard deviations of the
latitude and longitude of the optimal trajectories of each aircraft as
functions of time, and the expected values and standard deviation of the
timing of the trajectory as functions of the orthodromic distance from the
departure locations. They include the expected values and the standard
deviations of the arrival times and, in the case of formation flight, the
expected values and the standard deviations of the latitude, longitude, and
time of the rendezvous and splitting locations.

The sequence of discrete states obtained in the solution is the following:

  * •

State 1: All the aircraft fly solo.

  * •

State 2: Aircraft 1 and Aircraft 2 fly in a two-aircraft formation and
Aircraft 3 flies solo.

  * •

State 3: All the aircraft fly in a three-aircraft formation.

  * •

State 4: Aircraft 1 and Aircraft 3 fly in a two-aircraft formation and
Aircraft 2 flies solo.

  * •

State 5: All the aircraft fly solo.

TABLE I: Experiment A: Boundary conditions of the three flights. Symbol | Units | Flight 1 | Flight 2 | Flight 3
---|---|---|---|---
$\phi_{I}$ | [deg] | 40.64 | 42.36 | 45.47
$\phi_{F}$ | [deg] | 48.85 | 40.48 | 51.47
$\lambda_{I}$ | [deg] | -73.78 | -71.06 | -73.74
$\lambda_{F}$ | [deg] | 2.35 | -3.57 | -0.12
$\chi_{I}$ | [deg] | +54.26 | 69.25 | +55.53
$V_{I}$ | [m/s] | 240 | 240 | 240
$V_{F}$ | [m/s] | 220 | 220 | 220
$m_{I}$ | [kg] | $215\,000$ | $210\,000$ | $220\,000$
Fig. 2: Experiment A: Expected values of the latitude and longitude of the
optimal trajectories of each aircraft together with the corresponding 95%
confidence envelopes represented on the relevant map together with the wind
field.

The expected values of the longitude and latitude of the optimal routes
obtained in the solution together with the corresponding 95% confidence
envelopes are represented on the relevant map together with the wind field in
Fig. 2 and as functions of time in Fig. 3.a and Fig. 3.b, respectively. The
expected values of the rest of the state variables, namely, the mass, the
heading angle, and the Mach number, together with the corresponding 95%
confidence envelopes are represented as functions of time in Fig. 3.c, Fig.
3.d, and Fig. 3.e, respectively. The expected values of the control variables
are represented for the three flights in Fig. 4. Green, blue, and red lines
correspond to Flight 1, Flight 2, and Flight 3, respectively. For a better
understanding of these figures, the portions of the plots of the variables
that correspond to formation flight are represented on a gray background. In
particular, the portions of the plots that correspond to State 2 and State 4,
which includes a two-aircraft formation, are represented on a light gray
background, whereas the portion of the plots that correspond to State 3, which
includes a three-aircraft formation, is represented on a dark gray background.

It is important to note that all the optimal solutions obtained considering
the fuel burn savings that correspond to the nodes of the gPC expansion of the
random variable $\theta\sim\mathcal{N}(0.1,0.02)$ are formation missions with
the same structure.

The expected values of the rendezvous and splitting times are reported in
Table II together with the corresponding 95% confidence interval. State 2, in
which Aircraft 1 and Aircraft 2 fly in formation and Aircraft 3 flies solo,
has the shortest duration. The control variables vary in a quite smooth way.

(a) Longitude (b) Latitude (c) Mass (d) Heading angle (e) Mach number

Fig. 3: Experiment A: Expected values of the state variables of the optimal
trajectories of each aircraft together with the corresponding 95% confidence
envelopes.

(a) Lift coefficient (b) Bank angle (c) Adimensional thrust

Fig. 4: Experiment A: Expected values of the control variables of each
aircraft.

Regarding the spatial variability of the solution, it can be observed in Fig.
2, Fig. 3.a, and Fig. 3.b, that both geographical coordinates have a similar
behavior. In the first two discrete states, the spatial variability is rather
low, and in State 3, it is almost negligible. Then, near the first splitting
location, the spatial variability starts increasing for all the flights.
Finally, it decreases when the aircraft are approaching the final locations.
Hence, the second splitting location has a significant dependency on the
random variable $\theta$ and, consequently, both the duration and the traveled
distance in State 3 have a significant dependency on it. Flight 2 has the
greatest spatial variability.

In addition to the spatial variability, the temporal variability is also
quantified in order to have complete information regarding the spatio-temporal
variability of the solution. The expected values of the timing of the solution
together with the corresponding 95% confidence envelopes are represented in
Fig. 5 as functions of the orthodromic distance from the departure location of
each flight. It can be observed in this figure that the temporal variability
is very low in all the flights. In particular, in Flight 1 and Flight 3, the
temporal variability is almost negligible during the whole flight time,
whereas, in Flight 2, it starts increasing after the first splitting location.
However, despite this increase, it is very low in Flight 2 as well, leading to
the conclusion that the temporal variability of all the flights slightly
depends on the random variable $\theta$.

(a) Flight 1 (b) Flight 2 (c) Flight 3

Fig. 5: Experiment A: Expected values of the timing of the optimal
trajectories of each aircraft together with the corresponding 95% envelopes.

It can be observed in Table II, in which the expected values and the 95%
confidence intervals of the rendezvous and splitting times are reported, that
the amplitudes of these intervals are lower for the rendezvous times than for
the splitting times, but remains fairly small for both of them.

TABLE II: Experiment A: Expected values and 95% confidence intervals of the rendezvous and splitting times. | Expected | 95% confidence |
---|---|---|---
| value | interval |
First rendezvous time [h] | 0.88 | [0.87, 0.89] |
Second rendezvous time [h] | 1.32 | [1.32, 1.32] |
First splitting time [h] | 4.19 | [4.17, 4.21] |
Second splitting time [h] | 5.95 | [5.92, 5.98] |

The expected values and the 95% confidence intervals for both the flight time,
expressed in hours, and the fuel burn, expressed in tonnes, for each flight
are listed in Table III. From the amplitudes of the 95% confidence intervals,
it is easy to see that, as already mentioned, the flight time variability is
rather low. In contrast, as expected, the uncertainty in the fuel burn savings
has a significant impact on the variability of the fuel consumption.

TABLE III: Experiment A: Expected values and 95% confidence intervals of flight times and fuel consumptions of each flight. |  Flight time [$\boldsymbol{h}$] |  Fuel burn [$\boldsymbol{t}$] |
---|---|---|---
| Expected value | 95% confidence interval | Expected value | 95% confidence interval |
Flight 1 | 7.13 | [7.12, 7.15] | 42.47 | [41.19, 43.74] |
Flight 2 | 7.47 | [7.43, 7.52] | 46.24 | [45.90, 46.58] |
Flight 3 | 6.16 | [6.15, 6.18] | 38.00 | [36.58, 39.42] |

To quantify the benefits of formation flight with respect to solo flight, the
expected values of the flight time and the fuel consumption, along with the
expected value of the corresponding DOC expressed in monetary units, mu, are
estimated assuming that the three flights are performed as solo flights. The
obtained results are reported in Table IV. For the sake of comparison, the
expected values of the DOC obtained in both formation and solo flights of the
three aircraft are reported in Table V. It is remarkable that, in spite of the
uncertainties considered in the fuel burn savings for the trailing aircraft,
the reduction in the total DOC achieved with formation flight amounts to
$2.16\%$.

TABLE IV: Experiment A: Expected values of flight times, fuel consumptions, and DOC in solo flight of the three aircraft. |  Flight time [$\boldsymbol{h}$] |  Fuel burn [$\boldsymbol{t}$] |  DOC [$\boldsymbol{mu}$]
---|---|---|---
Flight 1 | 7.06 | 45.73 | 39635.79
Flight 2 | 7.32 | 45.44 | 39713.60
Flight 3 | 6.02 | 39.52 | 34165.60
TABLE V: Experiment A: Expected values of the DOC in solo and formation flights of the three aircraft. | Flight 1 | Flight 2 | Flight 3 | Total
---|---|---|---|---
DOC solo flight [$\boldsymbol{mu}$] | 39635.79 | 39713.60 | 34165.60 | 113514.99
DOC formation flight [$\boldsymbol{mu}$] | 37429.40 | 40435.59 | 33252.80 | 111117.79
$\boldsymbol{\Delta}$DOC [$\boldsymbol{\%}$] | -5.89 | +1.79 | -2.75 | -2.16

To quantify the effects of the uncertainty in the fuel burn savings for the
trailing aircraft on the DOC, the same formation mission is designed in the
absence of uncertainty, in which the fuel burn savings for both the trailing
and intermediate aircraft are set to 10%. For the sake of comparison, the
obtained results in terms of the DOC for the deterministic and stochastic
formation missions are summarized in Table VI. It can be seen that the
increment of the total DOC due to the presence of uncertainties amounts to
1.60%.

TABLE VI: Experiment A: Values of the DOC in deterministic and stochastic formation missions of the three aircraft. | Flight 1 | Flight 2 | Flight 3 | Total
---|---|---|---|---
DOC deterministic formation mission [$\boldsymbol{mu}$] | 36178.9 | 40205.67 | 32980.96 | 109365.53
DOC stochastic formation mission [$\boldsymbol{mu}$] (expected values) | 37429.40 | 40435.59 | 33252.80 | 111117.79
$\boldsymbol{\Delta}$DOC [$\boldsymbol{\%}$] | +3.45 | +0.57 | +0.82 | +1.60

###  VI-B Experiment B: Two-aircraft transoceanic mission design with
uncertainty in the departure times of the flights.

Experiment B involves two transoceanic eastbound flights, Flight 1 and Flight
2, with uncertainty in the departure times of the flights. The flights
considered in this experiment have the same departure and arrival locations
and the same boundary values for the state variables as Flight 1 and Flight 2
in Experiment A. They have the following departure and arrival locations:

  * •

Flight 1: New York (JFK) - Paris (CDG).

  * •

Flight 2: Boston (BOS) - Madrid (MAD).

The departure times of Flight 1 and Flight 2 are set to 10:15 and 10:30,
respectively. The boundary conditions for the state variables of the two
aircraft are given in Table I. Flight 1 and Flight 2 are operated by Aircraft
1 and Aircraft 2, respectively.

As in Experiment A, the formation configuration is selected in advance. In the
case of formation flight, Aircraft 1 is the trailing aircraft and Aircraft 2
is the leader. In this experiment, the fuel burn savings for the trailing
aircraft are set to 10%.

The random variables that represent the delays with respect to the scheduled
departure times of Flight 1 and Flight 2 from JFK and BOS airports are denoted
by $\theta_{1}$ and $\theta_{2}$, respectively. These random variables have
been modeled using a Gaussian mixture distribution, the parameters of which
are estimated from real delay data obtained from the Transtats online
database222https://www.transtats.bts.gov/  of the U.S. Bureau of
Transportation Statistics. This database provides detailed information on the
U.S. transportation systems, including data on departure delays by airport and
airline. The expectation-maximization algorithm for fitting mixture-of-
Gaussian models to data is employed for this purpose [24]. Departure delay
data corresponding to January 2020 are used.

The resulting Gaussian mixture probability density functions that model the
random variables $\theta_{1}$ and $\theta_{2}$ that represent the departure
delays at JFK and BOS airports both have 4 components. The values of the
parameters of the Gaussian mixture probability density functions are given in
Table VII, where $p_{i,j},i=1,2,j=1,2,3,4$ are the weights of the Gaussian
component probability density functions of $\theta_{i}$, which are
probabilities that sum to 1, and $\mu_{i,j}$ and $\sigma_{i,j}$ are their
means and standard deviations, respectively. This result is consistent with
previous studies on estimation of flight departure delay distributions [25].
Fig. 6 shows in black the obtained Gaussian mixture distribution of the
departure delays of Flight 1 and Flight 2 from JFK and BOS airports,
respectively, together with the histogram of the data. A bin width of 2
minutes is used for the histograms. The components of the Gaussian mixture
distribution are represented by dotted lines. The mean and the standard
deviation of the random variable $\theta_{1}$ are -1.67 [min] and 7.69 [min],
respectively. The mean and the standard deviation of the random variable
$\theta_{2}$ are -0.88 [min] and 10.37 [min], respectively.

TABLE VII: Experiment B: Estimated parameters of the Gaussian mixture probability density functions that model the departure delays of the flights.  | Departure airport | Random variable | $\boldsymbol{p_{i,1},p_{i,2},p_{i,3},p_{i,4}}$ | $\boldsymbol{\mu_{i,1},\mu_{i,2},\mu_{i,3},\mu_{i,4}}$ | $\boldsymbol{\sigma_{i,1},\sigma_{i,2},\sigma_{i,3},\sigma_{i,4}}$ |
---|---|---|---|---|---|---
Flight 1 | JFK | $\theta_{1}$ | 0.39, 0.17, 0.27, 0.17 | -4.94, 11.94, -0.99, -8.91 | 2.20, 7.17, 2.93, 2.89 |
Flight 2 | BOS | $\theta_{2}$ | 0.24, 0.56, 0.06, 0.14 | -1.76, -6.61, 28.67, 10.92 | 3.39, 3.70, 5.68, 5.78 |

(a) $\theta_{2}$, departure delay of Flight 1 from JFK airport (b)
$\theta_{2}$, departure delay of Flight 2 from BOS airport

Fig. 6: Experiment B: Estimated probability density functions of the random
variables $\theta_{1}$ and $\theta_{2}$ that represent the departure delays of
each flight.

As in Experiment A, solving this problem entails deciding which mode of flight
is optimal, the expected values and standard deviations of the latitude and
longitude of the optimal trajectories of each aircraft as functions of time,
and the expected values and standard deviation of the timing of the trajectory
as functions of the distance.

The sequence of discrete states obtained in the solution is the following:

  * •

State 1: Both aircraft fly solo.

  * •

State 2: The aircraft fly in formation.

  * •

State 3: Both aircraft fly solo.

Fig. 7: Experiment B: Expected values of the longitude and latitude of the
optimal trajectories of each aircraft together with the corresponding 95%
confidence envelopes represented on the relevant map together with the wind
field.

The expected values of the longitude and latitude of the optimal routes
obtained in the solution together with the corresponding 95% confidence
envelopes are represented on the relevant map together with the wind field in
Fig. 7 and as functions of time in Fig. 8.a and Fig. 8.b, respectively. The
expected values of the rest of the state variables, namely, the mass, the
heading angle, and the Mach number, together with the corresponding 95%
confidence envelopes are represented as functions of time in Fig. 8.c, Fig.
8.d, and Fig. 8.e, respectively. The expected values of the control variables
are represented for the two flights in Fig. 9. Green and blue lines correspond
to Flight 1 and Flight 2, respectively. As in Experiment A, the portions of
the plots of the variables that correspond to State 2, in which the aircraft
fly in formation, are represented on a gray background.

It is important to note that all the optimal solutions obtained considering
the fuel burn savings that correspond to the nodes of the gPC expansion of the
two random variables are formation missions with the same structure.

The expected values of the rendezvous and splitting times are reported in
Table VIII together with the corresponding 95% confidence intervals. State 2,
in which the aircraft fly in formation, has the longest duration.

The control variables vary in a quite smooth way. Comparing the expected
values of the state variables and the corresponding 95% confidence envelopes
represented as functions of time in Fig. 8 with the corresponding ones of
Experiment A given in Fig. 3, the high variability in every state variable of
each flight is striking. In particular, the Mach number variability in the
early part of the flight is especially worthy of note. This is due to the fact
that anticipations and delays in the departure times of the flights are
compensated by increments and reductions of their velocities.

(a) Longitude (b) Latitude (c) Mass (d) Heading angle (e) Mach number

Fig. 8: Experiment B: Expected values of the state variables of the optimal
trajectories of each aircraft together with the corresponding 95% confidence
envelopes.

(a) Lift coefficient (b) Bank angle (c) Adimensional thrust

Fig. 9: Experiment B: Expected values of the control variables of each
aircraft.

Regarding the spatial variability of the solution, it can be observed in Fig.
7, Fig. 8.a, and Fig. 8.b that the longitude variability is greater in State
2, in which the aircraft fly in formation. In contrast, the latitude
variability has a minimum in State 2, which is located around the maximum
value of this geographical coordinate. Aircraft 2, the leading aircraft, has
the greatest spatial variability.

As in Experiment A, in addition to the spatial variability, the temporal
variability is also quantified in order to have complete information regarding
the spatio-temporal variability of the solution. The expected values of the
timing of the solution together with the corresponding 95% confidence
envelopes are represented in Fig. 10, as functions of the orthodromic distance
from the departure locations of each flight. It can be observed in this figure
that the temporal variability has a similar behavior for both flights,
remaining nearly constant throughout the flight. As in the case of spatial
variability, the temporal variability of the flights in Experiment B is
considerably greater than the temporal variability of the flights in
Experiment A. As expected, uncertainty in the departure times of the flights
has a significant impact on the temporal variability of the trajectories, more
than the uncertainty in the fuel burn savings for the trailing aircraft.

(a) Flight 1 (b) Flight 2

Fig. 10: Experiment B: Expected values of the timing of the optimal
trajectories of each aircraft together with the corresponding 95% envelopes.

It can be observed in Table VIII, in which the expected values and the 95%
confidence intervals of the rendezvous and splitting times are reported, that
the amplitudes of these intervals are similar for the rendezvous and splitting
times, being, in both cases, considerably greater than the amplitudes of the
95% confidence intervals of the rendezvous and splitting times obtained in
Experiment A.

TABLE VIII: Experiment B: Expected values and 95% confidence intervals of the rendezvous and splitting times.  | Expected value | 95% confidence interval
---|---|---
Rendezvous time [h] | 1.61 | [1.29, 1.94]
Splitting time [h] | 5.03 | [4.71, 5.34]

The expected values and the 95% confidence intervals for both the flight time,
expressed in hours, and the fuel burn, expressed in tonnes, for each flight
are listed in Table IX. From the amplitudes of the 95% confidence intervals,
it is easy to see that, as already mentioned, the flight time variability is
rather high. Additionally, the uncertainty in the departure times also has a
significant impact on the variability of the fuel consumption.

TABLE IX: Experiment B: Expected values and 95% confidence intervals of flight times and fuel consumptions of each flight. |  Flight time [$\boldsymbol{h}$] |  Fuel burn [$\boldsymbol{t}$] |
---|---|---|---
| Expected value | 95% confidence interval | Expected value | 95% confidence interval |
Flight 1 | 7.18 | [6.86, 7.50] | 43.58 | [42.52, 44.63] |
Flight 2 | 7.70 | [7.37, 8.02] | 46.69 | [46.27, 47.11] |

To quantify the benefits of formation flight with respect to solo flight, the
expected values of the flight time and the fuel consumption along with the
expected value of the corresponding DOC expressed in monetary units, mu, for
the two solo flights, given in the first two rows of Table IV, are compared to
the expected values of the DOC obtained in formation flight. The obtained
results are reported in Table X. It can be observed that, in the presence of
uncertainties in the departure times of the flights, the probability density
funcion of which is estimated from real departure delay data, small benefits
are expected in terms of the reduction of the total DOC, which amount to only
0.11%. However, in real-life scenarios, some mitigation measures could be
implemented to reduce the negative impact of flight delays on the benefits of
formation flight, such as adjusting the flight plan of the formation mission
before departure when one of the flights of a formation is delayed.

TABLE X: Experiment B: Expected values of the DOC in solo and formation flights of the two aircraft. | Flight 1 | Flight 2 | Total
---|---|---|---
DOC solo flight [$\boldsymbol{mu}$] | 39635.79 | 39713.60 | 79349.39
DOC formation flight [$\boldsymbol{mu}$] | 38260.39 | 40999.01 | 79259.39
$\boldsymbol{\Delta}$DOC [$\boldsymbol{\%}$] | -3.59 | +3.14% | -0.11

To quantify the effects of the uncertainty in the departure delays on the DOC,
the same formation mission is designed in the absence of uncertainty, in which
no delays are considered. For the sake of comparison, the obtained results in
terms of the DOC for the deterministic and stochastic formation missions are
summarized in Table XI. It can be seen that the increment of the total DOC due
to the presence of uncertainties amounts to 3.76%.

TABLE XI: Experiment B: Values of the DOC in deterministic and stochastic formation missions of the two aircraft. | Flight 1 | Flight 2 | Total
---|---|---|---
DOC deterministic formation flight [$\boldsymbol{mu}$] | 36178.90 | 40205.67 | 76384.57
DOC stochastic formation flight [$\boldsymbol{mu}$] (expected values) | 38260.39 | 40999.01 | 79259.39
$\boldsymbol{\Delta}$DOC [$\boldsymbol{\%}$] | +5.75 | +1.97 | +3.76

Based on the results of the two experiments, it is possible to conclude that
uncertainties have a significant impact on the formation flight benefits in
terms of the reduction of the DOC. In particular, the uncertainties in the
departure times have a far greater impact on the DOC than the uncertainties in
the fuel burn savings. These results demonstrate that formation flight is
economically beneficial in the presence of uncertainties in the fuel savings
for the trailing aircraft and in the departure times of the aircraft.

In the next section, a sensitivity analysis of the solution of this experiment
to the random variables that represent the departure times of the flights is
carried out. The purpose of this analysis is to quantify how much uncertainty
in each component of the solution is due to the different sources of
uncertainty considered in the experiment.

##  VII Sensitivity Analysis

In this section, a variance-based sensitivity analysis of the results obtained
in Experiment B is conducted. The aim of the variance-based sensitivity
analysis is to quantify what proportion of the variance of the latitude and
the longitude of flight $k$ of the solution of the formation mission design
problem, $\lambda_{k}(t,\theta)$ and $\phi_{k}(t,\theta)$, respectively, is
due to the variance of each departure delays of Flight 1 and Flight 2,
$\theta_{1}$ and $\theta_{2}$, respectively.

In this paper, this analysis relies on the computation of the so-called Sobol’
indices, under the assumption that the random variables in $\theta$ are
independent [11]. Sobol’ indices can be directly derived from the coefficients
$C_{m}(t)$ of the gPC expansion (9) as explained in [26].

(a) $S^{\text{lon}}_{1,1}$ (gray), $S^{\text{lon}}_{1,2}$ (black) (b)
$S^{\text{lat}}_{1,1}$ (gray), $S^{\text{lat}}_{1,2}$ (black) (c)
$S^{\text{lon}}_{2,1}$ (gray), $S^{\text{lon}}_{2,2}$ (black) (d)
$S^{\text{lat}}_{2,1}$ (gray), $S^{\text{lat}}_{2,2}$ (black)

Fig. 11:  Experiment B: Sobol’ indices of the geographical coordinates of the
optimal trajectories of each flight. Gray and black lines correspond to the
random variables $\theta_{1}$ and $\theta_{2}$, respectively.

The analysis of the sensitivity of the components of the solution obtained in
Experiment B to the two random variables $\theta_{1}$ and $\theta_{2}$ that
represent the departure delays of Flight 1 and Flight 2, respectively, is
carried out employing this approach. The Sobol’ indices of the latitude and
longitude of Flight $k,k=1,2$ with respect to the random variable
$\theta_{j},j=1,2$ are denoted by $S^{\text{lat}}_{k,j}$ and
$S^{\text{lon}}_{k,j}$, respectively. They represent the proportions of the
variance of the geographical coordinates of the optimal route of Flight $k$
that is due to the variance of the departure time of Flight $j$. The Sobol’
indices of the longitude and latitude of the optimal routes obtained in the
solution of Experiment B are represented as functions of time in Fig. 11, in
which Sobol’ indices associated with $\theta_{1}$ and $\theta_{2}$ are plotted
in gray and black, respectively.

It can be seen in Fig. 11.a and Fig. 11.b that during the whole flight, most
of the variance of both geographical coordinates of the optimal route of
Flight 1 is due to the variance of the random variable $\theta_{1}$. In
particular, the percentage of the variance of the longitude of Flight 1 due to
the variance of $\theta_{1}$ is around 90%, while the percentage of the
variance of the latitude of Flight 1 due to the variance of the same random
variable is between 60% and 90% during the whole flight time.

It can be seen in Fig. 11.c and Fig. 11.d that during the whole flight, most
of the variance of the longitude of the optimal route of Flight 2 is due to
the variance of the random variable $\theta_{1}$. In contrast, the relative
influence of the variance of the random variables $\theta_{1}$ and
$\theta_{2}$ on the variance of the latitude of the optimal route of Flight 2
changes along the flight. In particular, in the first part of Flight 2, up to
8000 seconds, approximately, the variance of the random variable $\theta_{2}$
has the most influence on the variance of the latitude of Flight 2. After
that, there is an intermediate part of Flight 2, between 8000 and 14000
seconds, approximately, in which the variances of $\theta_{1}$ and
$\theta_{2}$ have a similar influence on the variance of the latitude of
Flight 2. Finally, in the last part of Flight 2, the influence of the variance
of $\theta_{1}$ on the variance of the latitude of Flight 2 becomes
predominant, reaching a percentage of 90%.

The Sobol’ indices of the timing of Flight $k,k=1,2$ with respect to the
random variable $\theta_{j},j=1,2$ are denoted by
$S^{\text{timing}}_{k,j}(d_{k})$, where $d_{k}$ is the orthodromic distance of
flight $k$ from the departure location. They numerically quantify the
proportion of the variance of the timing of the optimal route of flight $k$
that is due to the variance of $\theta_{j}$. The Sobol’ indices of the timing
of the optimal routes obtained in the solution of Experiment B are represented
as functions of the orthodromic distance in Fig. 12, in which Sobol’ indices
associated with $\theta_{1}$ and $\theta_{2}$ are plotted in gray and black,
respectively. It can be seen in this figure that, during the whole flight,
most of the variance of the timing of the optimal route of Flight 1 is due to
the variance of the random variable $\theta_{1}$, with a percentage of about
80%. On the contrary, the relative influence of the variance of the random
variables $\theta_{1}$ and $\theta_{2}$ on the variance of the timing of the
optimal route of Flight 2 changes along the flight. In particular, in the
first part of Flight 2, up to the distance of 1000 km, approximately, the
variance of $\theta_{2}$ has the most influence on the variance of the timing
of the optimal route of Flight 2. After that, the influence of the variance of
$\theta_{1}$ on the variance of the timing of the optimal route of Flight 2
becomes predominant, reaching a percentage of 80%.

(a) $S^{\text{timing}}_{1,1}$ (gray), $S^{\text{timing}}_{1,2}$ (black),
Flight 1 (b) $S^{\text{timing}}_{2,1}$ (gray), $S^{\text{timing}}_{2,2}$
(black), Flight 2

Fig. 12: Experiment B: Sobol’ indices of the timing of the optimal
trajectories of each flight. Gray and black lines correspond to the random
variables $\theta_{1}$ and $\theta_{2}$, respectively.

##  VIII Conclusions

This paper has presented a methodology for solving the formation mission
design problem for commercial aircraft in the presence of uncertainties in the
fuel burn savings of the trailing aircraft and in the departure times of the
aircraft. The solution of the problem is stochastic and the optimal aircraft
trajectories are given in terms of the expected values and standard deviations
of their latitude, longitude, and timing. The expected value of the direct
operating costs of the formation mission is also estimated. Moreover, the
proposed technique permits the relative contribution of each random variable
to the variability of the components of the solution to be estimated. This
methodology has been applied to plan two- and three-aircraft formation
missions with uncertainties. Specifically, uncertainty in the fuel burn
savings is considered in Experiment A, whereas uncertainty in the departure
times of the flights is included in Experiment B. The results of the numerical
experiments indicate that uncertainties have a significant impact on the
potential benefits of a formation mission. In particular, uncertainties in the
departure times have a greater impact on the direct operating costs than the
uncertainties in the fuel burn savings. They reveal that the increment of the
direct operating costs due to the presence of uncertainties amounts to 1.60%
and 3.76% in Experiment A and Experiment B, respectively. However, even in the
presence of uncertainties, expected reductions of the direct operating costs
amounting to 2.16% and 0.11% are achieved with formation flight compared to
solo flight in Experiment A and Experiment B, respectively. The results give
interesting additional information, revealing, for instance, that one random
variable may have predominant influence on the variability of a component of
the solution. This occurs in Experiment B for the variability of the longitude
of the optimal routes of both flights, the variance of which is mostly due to
the variance of the random variable that represents the departure time of
Flight 1. They also indicate that the relative influence of the random
variables on the variability of a component of the solution may change during
the flight. This happens in Experiment B for the variances of both the
latitude and timing of Flight 2. The results of the numerical experiments
demonstrate that formation flight is economically beneficial in the presence
of realistic levels of uncertainty in the fuel savings for the trailing
aircraft and in the departure times of the aircraft and that the proposed
methodology is an effective tool for solving the formation mission design
problem for commercial aircraft in the presence of uncertainties.

##  IX Acknowledgments

This work has been partially supported by the grants number TRA2017-91203-EXP
and RTI2018-098471-B-C33 of the Spanish Government.

## References

  * [1] M. Cerezo-Magaña, A. Olivares, and E. Staffetti, “Formation mission design for commercial aircraft using switched optimal control techniques,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 57, no. 4, pp. 2540–2557, 2021.
  * [2] International Civil Aviation Organization, “Convention on International Civil Aviation, Annex 2,” ICAO, Tech. Rep., 1944.
  * [3] D.-G. Caprace, G. Winckelmans, P. Chatelain, and J. D. Eldredge, “Wake vortex detection and tracking for aircraft formation flight,” in _Proceedings of the AIAA Aviation 2019 Forum_ , Dallas, TX, USA, June 2019.
  * [4] R. Abeyratne, “The outcome of the 40th ICAO assembly: A new look at ICAO ?” _Air and Space Law_ , vol. 45, no. 1, 2020.
  * [5] Airbus, https://www.airbus.com/en/newsroom/press-releases/2021-11-airbus-and-its-partners-demonstrate-how-sharing-the-skies-can-save, accessed: 2021-11-18.
  * [6] T. C. Flanzer, S. R. Bieniawski, and J. A. Brown, “Advances in cooperative trajectories for commercial applications,” in _AIAA SciTech 2020 Forum_ , 2020, p. 1257.
  * [7] D. Rivas and R. Vazquez, “Uncertainty,” in _Complexity Science in Air Traffic Management_ , A. Cook and D. Rivas, Eds. Routledge, 2016.
  * [8] M. S. Alwan and X. Liu, _Theory of Hybrid Systems: Deterministic and Stochastic_. Springer, 2018.
  * [9] Y. Zhu, _Uncertain Optimal Control_. Springer, 2019.
  * [10] D. Xiu, _Numerical Methods for Stochastic Computations. A Spectral Method Approach_. Princeton University Press, 2010\.
  * [11] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, and S. Tarantola, _Global Sensitivity Analysis. The Primer_. John Wiley & Sons, 2008.
  * [12] R. Shone, K. Glazebrook, and K. G. Zografos, “Applications of stochastic modeling in air traffic management: Methods, challenges and opportunities for solving air traffic problems under uncertainty,” _European Journal of Operational Research_ , vol. 292, pp. 1–26, 2021.
  * [13] J. Hu, M. Prandini, and S. Sastry, “Aircraft conflict prediction in the presence of a spatially correlated wind field,” _IEEE Transactions On Intelligent Transportation Systems_ , vol. 6, no. 3, pp. 326 – 340, 2005.
  * [14] M. Prandini and J. Hu, “Application of reachability analysis for stochastic hybrid systems to aircraft conflict prediction,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 4, pp. 913 – 917, 2009.
  * [15] X. Li, P. B. Nair, Z. Zhang, L. Gao, and C. Gao, “Aircraft robust trajectory optimization using nonintrusive polynomial chaos,” _Journal of Aircraft_ , vol. 51, no. 5, pp. 1592–1603, 2014.
  * [16] Y. Matsuno, T. Tsuchiya, J. Wei, I. Hwang, and N. Matayoshi, “Stochastic optimal control for aircraft conflict resolution under wind uncertainty,” _Aerospace Science and Technology_ , vol. 43, pp. 77–88, 2015.
  * [17] T. Kent and A. Richards, “Accounting for the effect of ground delay on commercial formation flight,” in _Proceedings of the 2014 UKACC International Conference on Control_ , Loughborough, U. K., July 2014.
  * [18] Eurocontrol, “User manual for the Base of Aircraft Data (BADA),” Eurocontrol, Tech. Rep. 130416, 2013.
  * [19] D. P. Dee, S. M. Uppala _et al._ , “The ERA–Interim reanalysis: Configuration and performance of the data assimilation system,” _Quarterly Journal of the Royal Meteorological Society_ , vol. 137, no. 656, pp. 553–597, 2011.
  * [20] M. D. Buhmann, _Radial Basis Functions: Theory and Implementations_. Cambridge University Press, 2003.
  * [21] A. J. Cook, G. Tanner, and S. Anderson, “Evaluating the true cost to airlines of one minute of airborne or ground delay,” Eurocontrol, Tech. Rep., 2004.
  * [22] W. E. Hart, C. D. Laird, J.-P. Watson, D. L. Woodruff, G. A. Hackebeil, B. L. Nicholson, and J. D. Siirola, _Pyomo-Optimization Modeling in Python_. Springer, 2012.
  * [23] T. E. Kent and A. G. Richards, “Potential of formation flight for commercial aviation: Three case studies,” _Journal of Aircraft_ , vol. In press, 2021\.
  * [24] K. P. Murphy, _Machine Learning: A Probabilistic Perspective_. The MIT Press, 2012.
  * [25] Y. Tu, M. O. Ball, and W. S. Jank, “Estimating flight departure delay distributions—a statistical approach with long-term trend and short-term pattern,” _Journal of the American Statistical Association_ , vol. 103, no. 481, pp. 112–125, 2008.
  * [26] A. Trucchia, V. Egorova, G. Pagnini, and M. Rochoux, “On the merits of sparse surrogates for global sensitivity analysis of multi-scale nonlinear problems: Application to turbulence and fire-spotting model in wildland fire simulators,” _Commun. Nonlinear Sci. Numer. Simulat._ , vol. 73, pp. 120–145, 2019.

|  María Cerezo-Magaña is a Teaching Assistant of Aerospace Engineering and a
PhD student at the Universidad Rey Juan Carlos in Madrid, Spain. She received
her MSc degree in Aeronautical Engineering from the Universidad Politécnica de
Madrid. Her research is focused on deterministic and stochastic hybrid optimal
control applied to aircraft trajectory optimization.
---|---
|  Alberto Olivares is a Professor of Statistics and Vector Calculus at the
Universidad Rey Juan Carlos in Madrid, Spain. He received his MSc degree in
Mathematics and his BSc degree in Statistics from the Universidad de
Salamanca, Spain, and his PhD degree in Mathematical Engineering from the
Universidad Rey Juan Carlos. He worked with the Athens University of Economics
and Business. His research interests include statistical learning, stochastic
hybrid optimal control and model predictive control with applications to
biomedicine, robotics, aeronautics and astronautics.
---|---
|  Ernesto Staffetti is a Professor of Statistics and Control Systems at the
Universidad Rey Juan Carlos in Madrid, Spain. He received his MSc degree in
Automation Engineering from the Università degli Studi di Roma “La Sapienza,”
and his PhD degree in Advanced Automation Engineering from the Universitat
Politècnica de Catalunya. He worked with the Universitat Politècnica de
Catalunya, the Katholieke Universiteit Leuven, the Spanish Consejo Superior de
Investigaciones Científicas, and with the University of North Carolina at
Charlotte. His research interests include stochastic hybrid optimal control,
iterative learning control and model predictive control with applications to
robotics, aeronautics and astronautics.
---|---