# Formation Mission Design for Commercial Aircraft using Switched Optimal
Control Techniques

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

In this paper, the formation mission design problem for commercial aircraft is
studied. Given the departure times and the departure and arrival locations of
several commercial flights, the relevant weather forecast, and the expected
fuel savings during formation flight, the problem consists in establishing how
to organize them in formation or solo flights and in finding the trajectories
that minimize the overall direct operating cost of the flights. Since each
aircraft can fly solo or in different positions inside 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 and
logical constraints in disjunctive form establish the switching logic among
the discrete states of the system. The formation mission design problem has
been formulated as an optimal control problem of a switched dynamical system
and solved using an embedding approach, which allows switching decision among
discrete states to be modeled without relying on binary variables. The
resulting problem is a classical optimal control problem which has been solved
using a knotting pseudospectral method. Several numerical experiments have
been conducted to demonstrate the effectiveness of this approach. The obtained
results show that formation flight has great potential to reduce fuel
consumption and emissions.

###### Index Terms:

Formation Flight, Formation Mission Design, Commercial Aircraft, Switched
Systems, Optimal Control, Embedding Approach, Logical Constraints, Knotting
Pseudospectral Method.

##  I Introduction

The aviation sector represents a great challenge for the environment, being a
top-ten global GreenHouse Gases (GHG) emitter with a 4.9$\%$ contribution to
the Earth's warming [1]. According to a recent special report of the
Intergovernmental Panel on Climate Change (IPCC) [2], a leading climate
science body, aviation is a crucial sector to the threat of climate change and
global warming. Indeed, although air transport is only responsible for 14% of
the dioxide carbon emissions of the transport sector, much less than other
transport modes like light-duty vehicles (30%) and heavy-duty vehicles (36%),
it seems to be harder to decarbonized. The reason is the air traffic predicted
growth rate, which is projected to be greater than other transport modes.
Indeed, the air traffic predicted growth rate is another concerning aspect.
The Compound Annual Growth Rate (CAGR) forecast over the next two decades,
assessed by the International Air Transport Association (IATA), estimates that
the number of flights is expected to continue increasing by the 3.5$\%$
annually, which means doubling the current number of passengers by 2040 [3].
This continuous growth in air traffic demand imposes the search for new
solutions to mitigate its potential negative effects on our planet.
Unfortunately, electrification of large aircraft is unlikely to happen anytime
soon due to problems related to power to weight ratio and power density levels
of the batteries [4].

Hence, the predicted effects of the aviation sector to climate change are
greater than desired [2], [5], and a great effort is required to mitigate
them. In this sense, formation flight has the potential to make a significant
contribution. Even more, formation flight offers a great promise not only in
mitigating the environmental impact of aviation, but also, in increasing the
Air Traffic Management (ATM) capacity. According to [6], “the formation
operates as a single aircraft with regard to navigation and position
reporting”, which means that two or three aircraft flying in formation would
be treated as one single aircraft for air traffic control purposes. This would
help to reduce their workload.

The considered formation mission design problem can be stated as follows:
given the departure times and the departure and arrival locations of several
commercial flights, the relevant weather forecast, and the expected fuel
savings during formation flight, the problem consists in establishing how to
organize them in formation or solo flights and in finding the trajectories
that minimize the overall Direct Operating Costs (DOC) of the flights. The
rendezvous and splitting locations and times must also be determined.  DOC are
all costs incurred directly in the operation of the aircraft. These expenses
are usually divided into three groups: flight operating costs, such as flight
crew salaries, fuel and oil, and airport and en-route charges, costs related
to maintenance and overhaul, and aircraft depreciation [7].

Only formations of up to three aircraft have been considered in this paper.
This choice is motivated by the fact that the difficulties in synchronizing
more than three flights make formations of more than three aircraft
operationally impractical. Besides, increasing the number of aircraft in the
formation reduces asymptotically the benefits obtained from the formation [8].

In this paper, the formation mission design problem has been formulated as an
optimal control problem of a switched dynamical system. These 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, each aircraft has different flight modes, namely solo flight and
flight in different positions inside 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 will be represented by
different dynamical equations which may include or not formation flight
benefits in terms of fuel 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.

###  I-A Previous Approaches

The notion of formation flight arose from the observation of nature, in which
many species of bird fly in formation during their migrations [9]. Later, it
was demonstrated that formation flight led to notable efficiency improvements
for the birds and that substantial fuel savings were also possible for
aircraft flying in formation. There are three types of formations: the echelon
or in-line formation, the V formation and the inverted V formation. These
first studies mainly focused on the aerodynamic of formation flight to assess
the potential benefits in terms of induced drag reductions and, consequently,
in potential fuel burn reductions [8, 10, 11]. Experimental studies on
formation flight have also been carried out. One of the most exhaustive works
has been conducted within the Surfing Aircraft Vortices for Energy ($AVE)
project [12]. The results obtained using two C-17s aircraft, at longitudinal
distances of between 18 and 70 wingspans, have shown fuel savings of up to 11%
[13, 14] for the trailing aircraft. The results of another experimental study
[15], obtained using two F/A-18s flying closer, at a longitudinal distance of
up to 6.6 wingspans, have shown fuel savings of just over 18% for the trailing
aircraft.

Besides these studies on aerodynamics of formation flight focused on
estimating the potential induced drag and fuel burn reductions, other research
activities have been conducted to study different practical aspects of
formation flight, such as the formation tracking control problem under
actuator and sensor faults [16] and the collision avoidance in formation
flight [17].

Regarding the optimal positioning, the lateral and the stream-wise distance
gain significance, in particular, the longitudinal spacing between aircraft,
usually expressed in terms of wingspans, is used to classify formations. In
the literature, the formations in which the longitudinal distance between
leader and trailing aircraft is less than 10 wingspans are called close
formations whereas those in which this distance is between 10 and 40 wingspans
are called extended formations [18]. Close formations provide the greatest
benefits in terms of induced drag reduction while, in extended formations, the
potential induced drag reduction is lower, especially in separations of more
than 20 wingspans. However, extended formations are intrinsically safer than
close formations, a crucial aspect in commercial flights.

Due to its potential in reducing fuel consumption and GHG emissions, formation
flight has recently attracted increasing interest and the studies on
aerodynamics of formation flight have given way to studies on formation flight
planning [19, 20, 21, 22, 23, 24], most of which focus on extended formations.
In these works, the formation flight planning problem has been solved as a bi-
level optimization problem, in which the higher level problem is the optimal
partner allocation problem and the lower level problem is the mission design
problem. While the former problem consists in establishing how to group the
elements of a large-size set of flights into smaller subsets which contain
potential partners of beneficial formation missions in terms of overall DOC,
the latter problem consists in determining the optimal trajectories of the
aircraft that compose candidate formation missions.

Although the mission design and the partner allocation problems are
intrinsically very different, due to the existing interconnection between
them, they have been addressed in one go in several works such as [19] and
[20], in which the Breguet range equation and the Fermat point extension
problem have been employed, respectively. Other examples are [21, 22], in
which geometric methods have been used. In these approaches, the inherent
combinatorial complexity of the partner allocation problem, which quickly
grows with the number of flights, forced to introduce simplifications, such as
using approximated dynamic models of the aircraft and neglecting the
meteorological forecast, which are essential factors to achieve a satisfactory
predictability of the trajectories and realistic estimations of the fuel
savings. The natural framework to take into consideration all these elements
is optimal control. In this work only the mission design problem will be
addressed using optimal control techniques and accurate dynamic models of the
aircraft. The partner allocation problem is left for future research.

Even though it is possible, choosing the type of formation and the relative
position of each aircraft in the formation has not been included in the
formulation of the formation mission design problem. Thus, without loss of
generality, the echelon formation has been assumed and the position of each
aircraft in the formation has been established in advance.

In [23] the mission design problem has been studied for two-aircraft
formations using a multi-phase optimal control approach. In this case, the
number of phases and the structure of the formation missions is unique and
known in advance. There are five phases. The first two correspond to the
flights from the departure locations to the rendezvous point, after which the
formation flight, which corresponds to the third phase, starts. The other two
flight phases correspond to the flights from the splitting point to the
arrival locations. In this case, a solution is found using multi-phase optimal
control techniques and compared with the solution of a solo flight mission to
establish which gives rise to minimum DOC.  The same framework used in [23] is
extended in [24] to design three-aircraft formation missions. However, only in
the latter a wind field model, obtained from a weather database using a
polynomial regression, has been taken into account. In contrast to two-
aircraft formation missions, the structure of three-aircraft formation
missions is not unique and, consequently, many different cases must be
considered, including the case in which one or all the aircraft flight solo.
The obtained solutions must be then compared to establish which corresponds to
the minimum DOC. This drawback is due to the fact that the multi-phase optimal
control approach is only able to tackle optimal control problems with a fixed
switching structure, that is, with transitions among phases established in
advance.

This approach is feasible if the number of flights considered is low and the
type of formation and the relative position of each aircraft in the formation
are fixed. Otherwise, the number of combinations quickly grows with the number
of flights which makes the previous approach hardly scalable.

A more general framework is needed, which is able to design formation missions
solving a single optimal control problem in which neither the number of phases
nor the switching structure are known in advance. Moreover, this approach
should be easily scalable to design formation missions involving more than
three flights and should allow a more accurate aerodynamic model to be used.

###  I-B Contributions of the paper

In this paper, an innovative approach to the formation mission design problem
is proposed which is formulated as an optimal problem problem of a switched
dynamical system with logical constraints in disjunctive form and solved using
an embedding approach [25, 26].

The optimal control formulation of the problem allows accurate dynamic models
of aircraft and meteorological forecast to be included in the problem
formulation in order to improve the predictability of the trajectories and the
estimation of the fuel savings. This represents a first improvement with
respect to previous approaches [19], [20], [21], [22], and [23]. Indeed, in
the latter, which is based on optimal control, the meteorological forecast has
not been taken into account.

Modeling the formation mission as a switched dynamical systems allows the
multiphase approach adopted in [23] and [24] to be avoided. This represents
another improvement of the method presented in this paper with respect to
these previous approaches, since the embedding method employed to solve the
resulting optimal control problem of a switched dynamical system is able to
deal with both the switching dynamics of the system and the logical
constraints in disjunctive form in an efficient way. Additionally, the
switching logic among discrete states can be modeled without using binary
variables, the number of switches does not have to be established in advance,
and the optimal values of the switching times between discrete states are
obtained without introducing them as unknowns of the optimal control problem.
Therefore, the resulting problem is a classical optimal control problem, which
has been solved using a pseudospectral knotting method [27, 28].

The advantages in terms of computational time to solve the formation mission
design problem given by the method presented in this paper with respect the
multiphase approach have been quantified comparing the computational times
obtained solving the same instances of the problem with both techniques.

The embedding approach has been used in [29] to tackle logical constraints in
disjunctive form introduced to enforce separation among aircraft in a multiple
aircraft trajectory planning problem. However, in [29] there are no switches
in the dynamic model of the system, whereas, as already mentioned, in the
formation mission design problem studied in this paper, the embedding approach
is used to handle both, the switched dynamical system and the logical
constraints in disjunctive form.

Several numerical simulations have been carried out to show the effectiveness
of the proposed technique. In the first numerical experiment, a two-aircraft
transoceanic mission design problem has been solved, in which the departure
time of one flight has been left free. The benefits of formation flight with
respect to solo flights have also been quantified. In a second numerical
experiment, an analysis of the solution obtained in the previous experiment
has been carried out to determine how delays in the departure times of each
flight affect the formation flight. In the third numerical experiment, a
three-aircraft transoceanic mission design problem has been solved. An
analysis of the obtained solution has also been performed to determine how the
fuel savings scheme affects the formation flight.

The method for formation mission design presented in the paper is intended to
be incorporated in a ground-based system whose schematic diagram is given in
Fig. 1, which could be used by flight dispatchers as a support system to
design formation missions for airlines or airline alliances. Flight
dispatchers are employed by airlines that are responsible for creating and
modifying the flight plan. They also monitor the aircraft during flight. The
computational time of the proposed algorithm makes it suitable for the
strategical and tactical planning phases of the flight, as well as for the
operational phase, that is, during flight.

Fig. 1: Schematic diagram of the formation mission design system.

###  I-C Organization of the paper

The paper is organized as follows. In Section II the model of the aircraft is
introduced and the modeling of the wind forecast is described. In the same
section, the peculiarities of the formation mission design problem are
discussed. In Section III, the switched optimal control problem used to solve
the mission design problem is stated. In Section IV, the embedding approach
employed to solve the resulting switched optimal control problem is described.
In Section V, the technique employed to model logical constraints in
disjunctive form is presented and then particularized to model the logical
constraints arising in the mission design problem. In Section VI, the
pseudospectral knotting method used to solve the resulting classical optimal
control problem is outlined. The results of the numerical experiments are
reported and analyzed in Section VII. Finally, in Section VIII, some
conclusions are drawn.

##  II Model of the system

###  II-A Equations of motion

The mission design problem will be studied only in the cruise phase of the
flight under the assumption that, if the flights have the same origin or the
same destination, formation flight would be forbidden for security reasons
during departure, climb, descent, and approach phases of flight. A simplified
two-degree-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 has been considered and all the aircraft forces are
supposed to be in the plane of symmetry of the aircraft. Wind effects have
also been considered. All aspects associated to the rotational dynamics are
neglected. The set of kinematics and dynamics Differential-Algebraic Equations
(DAE) [30] that describe each aircraft motion is

$\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-A) has been used in
this paper.

The lift force is $L=qSC_{L}$, 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. A parabolic drag polar $C_{D}={C_{D}}_{0}+KC_{L}^{2}$ has been
assumed, where ${C_{D}}_{0}$ is the zero-lift drag component and $K$ is the
induced drag coefficient. Other parameters and variables include the Earth
radius $R_{E}$, the gravitational acceleration $g$, and the cruise altitude
$h$. ${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 to be added to the above set of DAE 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 [31], has been used to determine it. In
particular for jet engines the following equation applies

$\eta(t)=C_{f_{1}}\left(1+\frac{V(t)}{C_{f_{2}}}\right),$

where $C_{f_{1}}$ and $C_{f_{2}}$ are empirical thrust specific fuel
consumption coefficients.

Thus, for aircraft $p$, the state vector will be
$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 will be $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.

###  II-B Flight envelope

Flight envelope constraints model aircraft performance limitations. These
constraints involve, among others, flight altitude, load factor and airspeed,
and can be stated in the form $\eta_{l}(t)\leq\eta(t)\leq\eta_{u}(t)$, where
$\eta(t)$ is the state or control variable and $\eta_{l}(t)$ and $\eta_{u}(t)$
are the minimum and maximum allowed values of the variable, respectively. BADA
3.6 has been used to impose the following constraints

$\displaystyle V_{min}(t)\leq$ $\displaystyle V_{CAS}(t)$ $\displaystyle\leq
V_{MO},$ $\displaystyle m_{min}\leq$ $\displaystyle m(t)$ $\displaystyle\leq
m_{max},$ $\displaystyle M(t)$ $\displaystyle\leq M_{MO},$ $\displaystyle
T_{min}(t)\leq$ $\displaystyle T(t)$ $\displaystyle\leq T_{max}(t),$ (3)
$\displaystyle{C_{L}}_{min}\leq$ $\displaystyle C_{L}(t)$
$\displaystyle\leq{C_{L}}_{max},$ $\displaystyle-\mu_{max}\leq$
$\displaystyle\mu(t)$ $\displaystyle\leq\mu_{max},$

where $V_{CAS}$ is the calibrated airspeed, which can be calculated as a
function of the true airspeed [31]. $V_{min}$ and $V_{MO}$ are the minimum and
maximum operating calibrated speeds, respectively. $m_{min}$ and $m_{max}$ are
the minimum and maximum aircraft masses, respectively. $M$ is the Mach number
and $M_{MO}$ is the maximum operating Mach number. $T_{min}$ and $T_{max}$ are
the minimum and maximum available engine thrusts, respectively.
${C_{L}}_{min}$ and ${C_{L}}_{max}$ are the minimum and maximum lift
coefficients, respectively. Finally, $\mu_{max}$ is the maximum bank angle set
by the air navigation regulations in civil flight. Furthermore, the BADA
aircraft model includes the equation $V_{min}(t)=C_{V_{min}}V_{s}(t)$, which
relates $C_{V_{min}}$, the minimum speed coefficient, and $V_{s}$, the stall
speed for each flight phase.

For turbofan propelled aircraft, assuming standard atmosphere conditions, the
maximum thrust, $T_{max}$, is defined by the following empirical expression

$T_{max}(t)=C_{T_{cr}}\cdot
C_{T_{C,1}}\cdot\left[1-\frac{H_{p}(t)}{C_{T_{C,2}}}+C_{T_{C,3}}\cdot
H_{p}^{2}(t)\right],$

where $C_{T_{cr}}$ is the maximum cruise thrust coefficient. $C_{T_{C,1}}$,
$C_{T_{C,2}}$, and $C_{T_{C,3}}$ are the empirical thrust coefficients and
$H_{p}$ is the geopotential pressure altitude.

###  II-C Wind model

ERA-Interim [32] is a third-generation global atmospheric reanalysis produced
by the European Centre for Medium-Range Weather Forecasts (ECMWF). Climate
reanalysis consists in systematic and coherent assimilation of the global
meteorological data constrained by the available observations obtained from
different sources during the period of reanalysis, such as radiosondes or
aircraft, and prior state estimates using an invariant consistent assimilation
scheme and an integrated meteorological model. The third generation reanalysis
data refers to the last generation, which was developed in 2010s, and
significantly enhanced the spatio-temporal resolution of the former ones.

Gridded data products, which can be found in the ERA-Interim database, include
a wide variety of atmospheric parameters relevant to the mission design
problem, such as air temperature, pressure, and wind at different altitudes
and sea-surface temperature. In this study, the components of the wind
velocity in eastwards and northwards directions have been used.

ERA-Interim is presented as a regular latitude-longitude grid, with a spatial
resolution from $0.25^{\circ}\times 0.25^{\circ}$ to $3^{\circ}\times
3^{\circ}$ for the deterministic reanalysis, at 37 atmospheric levels
corresponding to different pressure levels from 1000 to 1 hPa. The temporal
coverage is four times daily, getting six-hour data. In this study, the chosen
pressure level has been 200 hPa, corresponding to the cruise altitude. The
selected spatial resolution has been $0.5^{\circ}\times 0.5^{\circ}$. Finally,
the selected day has been April 30, 2019, at 12:00. To include atmospheric
data in the hybrid optimal control problem used to solve the mission design
problem, an analytic function that approximates the data must be determined.
For this purpose, Radial Basis Functions (RBF) [33] have been used.

Given a set of known input data, a function approximation at any point is
obtained constructing a linear space which depends on the relative position of
the evaluated point with respect to the known observations according to an
arbitrary distance measure. RBF consider interpolating functions of the form

$F(x)=c_{0}+\sum_{i=0}^{N-1}{c_{c}}\varphi(\parallel x-R_{i}\parallel),$

where $x\in\mathbb{R}^{N}$ is the input data vector and $N$ is its dimension,
$F(x)$ is the function approximation, $\varphi$ is the basis function, and
$R_{i}$ is a vector containing the centers, which are the reference points of
the basis functions. $c_{0}$ and $c_{c}$ are constant coefficients calculated
by RBF method, in particular, $c_{0}$ is the bias term. The Euclidean norm is
considered to compute the distance between the input points and the centers.
There are different options for the choice of the basis functions. In this
paper, Gaussian basis functions have been selected.

Fig. 2: Eastward wind speed on 30 April 2019, 12:00. Fig. 3: Northward wind
speed on 30 April 2019, 12:00.

Eastwards and northwards components of the wind speed taken from ERA-Interim
and then approximated by the RBF method are represented in Fig. 2 and Fig. 3,
respectively.

###  II-D Formation flight model

The DAE system consisting of Eq. (II-A) and Eq. (2) is applicable to any
aircraft in the cruise phase of the flight. In solo flights, depending on
aircraft type, the aerodynamic coefficients, the wingspan, and the reference
surface, among others, will change. In formation flights, additionally, it is
necessary to include a modification in Eq. (2) for the trailing aircraft.

As previously mentioned, in extended formations, aircraft fly with a
longitudinal separation of more than 10 wingspans but the potential fuel burn
reduction for the trailing aircraft has a significant decrease beyond 20
wingspans separation. Therefore, in this paper, the formation benefits are
considered negligible for separations of more than 20 wingspans and the
constraint introduced to model the distance which leads to fuel burn
reductions in the trailing aircraft, and in the intermediate one in case of
three-aircraft formation, during formation flight is the following

$10b\leq\mathcal{D}_{pq}(t)\leq 20b,\quad\forall
p,q\in\\{1,\dots,N_{a}\\},p<q,$ (4)

where $b$ is the wingspan of the leader aircraft and $\mathcal{D}_{pq}(t)$ is
the orthodromic distance between aircraft $p$ and $q$, at time $t$, given by
the Haversine formula.

Induced drag reduction in formation flight is due to the fact that any
aircraft moving through a fluid generates lift by imparting a swirling motion
to the air, which creates downwash and upwash fields. An aircraft flying in
the upwash field created by another one will increase its apparent angle of
attack, getting an effective rotation of the lift vector resulting in a lower
lift-dependent component of drag. In this way, a significant reduction in
induced drag is obtained [8]. In Fig. 4 a schematic representation of the
vertical component of the induced air velocity field is given.

Fig. 4: Vertical component of the induced air velocity field (not to scale).
Figure adapted from [34].

A simplification has been introduced assuming 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 manner, the fuel saving factor can be directly incorporated in 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).$ (5)

It is quite difficult to estimate the drag reduction and, hence, the fuel burn
reduction achieved during a formation flight, since there are many parameters
to take into account. The number of aircraft involved in the formation, their
type, size, and weight, the ability of the trailing aircraft to place itself
in the optimum location of the wake, and the stresses and vibrations induced,
among others, strongly influence the benefits of the formation. In this paper,
different fuel savings schemes for the trailing and the intermediate aircraft
of a two- or three-aircraft formation have been considered. The adopted
percentages are between 6% to 14% of fuel savings for the trailing and
intermediate aircraft, with respect to fuel consumption in solo fights.

###  II-E Objective functional

The objective functional is a key factor in all trajectory optimization
problems, but it is crucial in formation mission design problems, in which
opposing optimality criteria, such as, fuel consumption and total mission
time, must be considered. On the one hand, formation flight can reduce the
total fuel consumption and, on the other hand, in general, it requires
increasing the total flight time to create the formation itself. Therefore,
formation mission design problems require to leave traditional fuel
minimization models behind. In this paper, fuel consumption and flight time
are combined in the objective functional using weighting coefficients. The
weighting coefficients have been chosen based on [35], in which a detailed
study on the costs incurred by airlines due to delays has been made. No other
costs have been considered.

The objective functional $J$ for $N_{a}$ aircraft, is defined as follows

$\min\quad
J=\alpha_{t}\sum_{p=1}^{N_{a}}t_{{flight}_{p}}+\alpha_{f}\sum_{p=1}^{N_{a}}{{{m}}_{f_{p}}},$
(6)

where $\alpha_{t}$ and $\alpha_{f}$ are the time and the fuel burn 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$. It is worth mentioning that the initial
fuel loaded in an aircraft that is going to benefit from formation flight
could be considerably lower than the amount of fuel loaded in the same
aircraft but without reasonable chances to benefit from a formation flight.
Nevertheless, in this paper, a more conservative approach has been adopted and
all aircraft are loaded with the same amount of fuel as in solo flight, even
if this fact will lead to lower benefits in terms of DOC.

As said before, in the approach to the formation mission design problem
presented in this paper, the relative position of each aircraft is set in
advance. Additionally, in three-aircraft formations, the only type of
formation allowed is the in-line formation. A maximum detour time from the
solo flight time has been set in order to perform a more realistic problem
formulation according to the usual operational requirements.

##  III The Switched Optimal Control Problem

As said in the introduction, in this paper, a novel approach to the formation
mission design problem is presented, which is formulated as an Optimal Control
Problem (OCP) for a switched dynamical system.

In this section, following [26], the formulation of a Switched Optimal Control
Problem (SOCP) is introduced for two-switched dynamical systems. The dynamical
model of a two-switched dynamical system can be represented by

$\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},$ (7) $\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}.$

For the sake of clarity, in the formulation of the SOCP, the sub-index “$S$”
has been added to the symbols that denote the state and control variables.  In
this model, the continuously differentiable vector fields,
$f_{0},f_{1}:\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$,
specify 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 identifies which of the two
possible system modes, $f_{0}$ or $f_{1}$, is active. Thus, both of them,
$u_{S}(t)$ and $v_{S}(t)$, can be regarded as control variables. The initial
time $t_{I}$, the initial state $x_{S}(t_{I})$, the final time $t_{F}$, and
final state $x_{S}(t_{F})$ are assumed to be restricted to a boundary set $B$
as follows: $(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 performance index
of the SOCP has the following form

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

where the endpoint cost function $g$ is defined on a neighborhood of $B$, and
$F_{0}$ and $F_{1}$ are real-valued continuously differentiable functions that
represent the running cost of operation of the system in each mode.

The SOCP is thus 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),$
(9)

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

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 to 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 will
obtain 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
intermediate or trailing in the formation, it will obtain benefits from the
formation flight and its mass flow will be described by Eq. (5). Thus, $f_{0}$
represents Equations (II-A) and (2), whereas $f_{1}$ represents Equations
(II-A) and (5). The functions $F_{0}$ and $F_{1}$, which appear in the running
cost of Eq. (8), do not exist in the formulation of this particular problem.

##  IV The Embedding Approach

In this section, the embedding approach employed to transform the resulting
SOCP into an Embedded Optimal Control Problem (EOCP) is presented. In this
way, a smooth OCP without binary variables is obtained, reducing the
computational complexity of finding the solution.

Following again [26], Eq. (III) can be rewritten as follows

$\displaystyle\dot{x}_{E}(t)$
$\displaystyle=[1-v_{E}(t)]f_{0}(t,x_{E}(t),u_{E_{0}}(t))+v_{E}(t)f_{1}(t,x_{E}(t),u_{E_{1}}(t)).$
(10)

For the sake of clarity, in the formulation of the EOCP, the sub-index “$E$”
has been added to the symbols that denote the state and control variables.
However, the components of each vector have not changed with respect to the
SOCP.

Equations (III) and (10) formally coincide under the conditions
$v_{E}(t)=v_{S}(t)$ and $u_{E_{0}}(t)=u_{E_{1}}(t)=u_{S}(t)$. However, in Eq.
(10), $v_{E}(t)\in[0,1]$, while in Eq. (III), $v_{S}(t)\in\\{0,1\\}$, which is
the key feature of the embedding approach [25, 26]. This method is based on
solving the SOCP defined above, using only continuous control variables
$v_{E}(t)\in[0,1]$, $u_{E_{0}}(t),u_{E_{1}}(t)\in\Omega$, in which the
performance index of the SOCP from Eq.(8) is rewritten as

$\displaystyle
J_{E}\left(t,x_{E}(t),u_{E_{0}}(t),u_{E_{1}}(t),v_{E}(t)\right)=g(t_{F},x_{E}(t_{F}))+$
$\displaystyle+\int_{t_{I}}^{t_{F}}\left([1-v_{E}(t)]F_{0}(t,x_{E}(t),u_{E_{0}}(t))+v_{E}(t)F_{1}(t,x_{E}(t),u_{E_{1}}(t))\right)dt.$
(11)

The EOCP is thus formulated as

$\min_{u_{E_{0}},u_{E_{1}}\in\Omega,v_{E}\in[0,1]}J_{E}\left(t,x_{E}(t),u_{E_{0}}(t),u_{E_{1}}(t),v_{E}(t)\right),$
(12)

subject to Eq. (10) and endpoint contraints
$(t_{I},x_{E}(t_{I}),t_{F},x_{E}(t_{F}))\in B$, which is an OCP without binary
variables. Therefore, standard numerical optimal control techniques can be
applied to solve it.

It has been shown in [25] and [26] that, once a solution of the EOCP has been
obtained, either the solution is of the switched type, that is, $v_{E}$ takes
only the values $0$ and $1$, or suboptimal trajectories of the SOCP can be
constructed that can approach the value of the cost for the EOCP arbitrarily
closely, and satisfy the boundary conditions within $\epsilon$, with arbitrary
$\epsilon>0$. A thorough discussion about the relationship between the
solutions of the SOCP and the EOCP can be found in [25] and [26].

###  IV-A Specification of the Embedded Optimal Control Problem

As said in Section III.A, only the aircraft's mass flow rate equation is
subject to switches. To model this switching, equations (2) and (5) are
combined using the selection variable
$v_{E_{p}}(t),\;p\in\left\\{1,\ldots,N_{a}\right\\}$ as follows

$\displaystyle\dot{\phi}_{p}(t)$
$\displaystyle=\quad\dfrac{V_{p}(t)\cdot\cos\chi_{p}(t)+{V_{W_{N}}}(t)}{R_{E}+h_{p}},$
$\displaystyle\dot{\lambda}_{p}(t)$
$\displaystyle=\quad\dfrac{V_{p}(t)\cdot\sin\chi_{p}(t)+{V_{W_{E}}}(t)}{\cos\phi_{p}(t)\cdot\left(R_{E}+h_{p}\right)},$
$\displaystyle\dot{\chi}_{p}(t)$
$\displaystyle=\quad\dfrac{L_{p}(t)\cdot\sin\mu_{p}(t)}{V_{p}(t)\cdot
m_{p}(t)},$ $\displaystyle\dot{V}_{p}(t)$
$\displaystyle=\quad\dfrac{T_{p}(t)-D_{p}(t)}{m_{p}(t)},$
$\displaystyle\dot{m}_{p}(t)$
$\displaystyle=\quad(1-v_{E_{p}}(t))\cdot[-T_{p}(t)\cdot\eta_{p}(t)]+$
$\displaystyle\hskip
21.33955pt+\;v_{E_{p}}(t)\cdot[-(1-\mathcal{R}_{\text{fuel}_{p}})\cdot
T_{p}(t)\cdot\eta_{p}(t)],$

$\forall p\in\left\\{1,\ldots,N_{a}\right\\}$, where $v_{E_{p}}(t)=1$
corresponds to Formation Flight in an in-line formation (FF) and
$v_{E_{p}}(t)=0$ to Solo Flight (SF).

Therefore, each aircraft has different flight modes, namely SF and FF in
different positions inside a formation. Their combination is represented by
the discrete state of the switched system which models their joint dynamic
behaviour.

However, in the mission design problems, the switching control variables
$v_{E_{p}}(t),\;p\in\left\\{1,\ldots,N_{a}\right\\}$ depend on the state
variables, in particular on the geographical coordinates of the aircraft.
Indeed, not all transitions between discrete states are possible at any time,
since while the distance between aircraft is larger than 20 wingspans, FF can
not take place. Furthermore, the switching control variables depend on one
another. For instance, in a two-aircraft formation, if one aircraft is
selected as the trailing aircraft the other one is forced to be the leader and
vice versa. Similarly, in a three-aircraft formation, if one position in the
formation is assigned to one aircraft, the other ones must fly in another
position.

Thus, the switched system that represents the joint dynamic behavior of all
the aircraft involved in the mission design problem should follow some rules
that are specified using logical constraints usually expressed in disjunctive
form. The method employed to model logical constraints in disjunctive form in
the OCP used to solve the mission design problem will be described in the next
section.

##  V Logical Constraints Modeling

To transform logical constraints in disjunctive form into inequality and
equality constraints that involve only continuous auxiliary variables, the
approach proposed in [36] has been employed.

It has been shown in [37] that every Boolean expression can be transformed
into Conjunctive Normal Form (CNF). Hence, there is no loss of generality in
considering that any logical constraint can be rewritten as a CNF expression

$Q_{1}\wedge Q_{2}\wedge\ldots\wedge Q_{n},$ (14)

in conjunctive form, where the symbol $\wedge$ denotes the $``$and$"$ operator
and

$Q_{i}=P_{i}^{1}\vee P_{i}^{2}\vee\ldots\vee P_{i}^{m_{i}},\;\forall
i\in\\{1,2,\ldots,n\\},$ (15)

where $Q_{i}$ is expressed in disjunctive form, in which the symbol $\vee$
denotes the $``$or$"$ operator and the proposition $P_{i}^{j}$ is either
$X_{i}^{j}$ or $\neg X_{i}^{j}$. The term $X_{i}^{j}$ can be either True or
False and the symbol $\neg$ represents the $``$not$"$ operator.

Term $X_{i}^{j}$ represents statements such as $``$$\mathcal{D}_{12}(t)\geq
20b\;$$"$. Therefore, $P_{i}^{j}$ takes the form

$P_{i}^{j}\equiv\left\\{g_{i}^{j}\left(x(t)\right)\leq 0\right\\},$ (16)

$\forall i\in\\{1,2,\ldots,n\\}$, $\forall j\in\\{1,2,\ldots,m_{i}\\}$, where
$g_{i}^{j}:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}$ is assumed to be a
$\mathcal{C}^{1}$ function.

To incorporate logical constraints in a OCP, they must be converted into a set
of equality or inequality constraints. Additionally, if binary variables are
not used, the combinatorial complexity of integer programming is avoided.
Conjunctions from Eq. (14) can be straightforwardly included into a OCP since
they are equivalent to $Q_{i},\forall i\in\\{1,2,\ldots,n\\}$. To be able to
transform the disjunctions into a set of inequality constraints, a continuous
variable $\alpha_{i}^{j}(t)\in[0,1]$ is defined and associated with each
$P_{i}^{j}$ in Eq. (16). Thus, Eq. (15) can be rewritten as

$\displaystyle\alpha_{i}^{j}(t)\cdot g_{i}^{j}\left(x(t)\right)\leq 0,$ and
$\displaystyle 0\leq\alpha_{i}^{j}(t)\leq 1,$ (17) and
$\displaystyle\sum_{j=1}^{m_{i}}\alpha_{i}^{j}(t)=1,$

$\forall i\in\\{1,2,\ldots,n\\},\forall j\in\\{1,2,\ldots,m_{i}\\}$. The first
constraint in (17) shows that when $\alpha_{i}^{j}(t)=0$ the constraint
$g_{i}^{j}\left(x(t)\right)\leq 0$ is not forced to be fulfilled. In contrast,
if $0<\alpha_{i}^{j}(t)\leq 1$, then $\alpha_{i}^{j}(t)\cdot
g_{i}^{j}\left(x(t)\right)\leq 0$, and, consequently, the constraint
$g_{i}^{j}\left(x(t)\right)\leq 0$ is forced to be fulfilled. The last
equation in (17) guarantees that at least one of the propositions $P_{i}^{j}$
holds.

In the next section, this approach will be applied to the formation mission
design problem in which logical constraints involve the longitudinal distance
between aircraft.

###  V-A Specification of the Logical Constraints

As stated in Eq. (4), in this paper, it has been considered that the great-
circle distance between aircraft in FF is between 10 and 20 wingspans, and is
greater than 20 wingspans in case of SF. Hence, the distance between two
aircraft, $\forall p,q\in\\{1,\dots,N_{a}\\},\;p<q$, should satisfy

$\displaystyle\;\;10b\leq$ $\displaystyle\mathcal{D}_{pq}(t)\;\;\leq 20b,$
(18) or $\displaystyle\mathcal{D}_{pq}(t)\;\;\geq 20b.$

Using simple Boolean algebra, Eq.(18) can be rewritten as

$\displaystyle\;\;10b\leq\mathcal{D}_{pq}(t),$ (19) and
$\displaystyle\quad\mathcal{D}_{pq}(t)\;\;\leq
20b\quad\text{or}\quad\mathcal{D}_{pq}(t)\;\;\geq 20b.$

To transform the disjunctions into a set of equality and inequality
constraints, a new variable is defined. Thus, in two-aircraft mission design
problems, only a continuous variable related to the longitudinal distance
$\mathcal{D}_{12}$ should be defined, whereas, in three-aircraft mission
design problems three continuous variables must be introduced, which are
related to the longitudinal distances
$\mathcal{D}_{12},\mathcal{D}_{13},\mathcal{D}_{23}$. Thus, calling the new
continuous variable $\alpha_{pq}(t)$ $\in[0,1]$, $\forall
p,q\in\\{1,\dots,N_{a}\\},\;p<q$, the logical constraints in Eq.(19) can be
rewritten as

$\displaystyle\alpha_{pq}(t)\left(\mathcal{D}_{pq}(t)-20b\right)\leq 0,$ and
$\displaystyle\left(1-\alpha_{pq}(t)\right)\left(20b-\mathcal{D}_{pq}(t)\right)\leq
0,$ and $\displaystyle 0\leq\alpha_{pq}(t)\leq 1,$ (20) and
$\displaystyle\mathcal{D}_{pq}(t)\geq 10b.$

For aircraft $p$ and $q$ $\in\\{1,\dots,N_{a}\\},\;p<q$, if
$\mathcal{D}_{pq}(t)>20b$, $\alpha_{pq}(t)$ must be zero in order to satisfy
the first constraint in (20). On the contrary, when $\mathcal{D}_{pq}(t)\leq
20b$, the variable $\alpha_{pq}(t)$ must be equal to one to fulfill the second
constraint in the above set of equations. Formation flight benefits will be
achieved only when $\alpha_{pq}(t)=1$. In any case, the last constraint in
(20) ensures that the longitudinal distance between aircraft will always be
greater than $10b$, the minimum safety distance.

It can be observed that nature of the variables $v_{E_{p}}(t)$, $\forall
p\in\\{1,\dots,N_{a}\\}$, and $\alpha_{pq}(t)$, $\forall
p,q\in\\{1,\dots,N_{a}\\},\;p<q$, is similar. While the variable
$v_{E_{p}}(t)$ selects equations of the dynamic model for aircraft $p$, the
variable $\alpha_{pq}(t)$ selects constraints. Although, for the sake of
clarity, a different notation has been used in the previous sections to
describe them, the variables $\alpha_{pq}(t)$ and $v_{E_{p}}(t)$ are directly
related by the following equivalence

$\forall p,q\in\\{1,\dots,N_{a}\\},\;p<q:\quad\alpha_{pq}(t)=v_{E_{p}}(t),$
(21)

where aircraft $p$ is supposed to be the trailing aircraft. Therefore, the
variable $v_{E_{q}}(t)$, which selects equations of the dynamic model for
aircraft $q$, must be zero regardless of the value of variable
$\alpha_{pq}(t)$, as aircraft $q$ will be the leader one and, hence, it will
achieve no benefits from the formation.

It is important to notice that the equivalence in Eq. (21) between the
continuous variables obtained from the embedding approach, used to model the
switched dynamical system and transform the logical constraints in disjunctive
form into inequality and equality constraints, need not necessarily to be
fulfilled in others EOCPs.

##  VI Solution of the Switched Optimal Control Problem

Once the SOCP has been reformulated as an EOCP and the logical constraints in
disjunctive form have been transformed into inequality and equality
constraints without binary variables using the embedding approach, the
resulting optimal control problem has been transcribed into a NonLinear
Programming (NLP) problem using a PeudoSpectral Method (PSM). More
specifically, the Radau PSM has been employed [28].

The set of constraints of the resulting NLP problem includes the system
constraints that correspond to the differential constraint (IV-A) and the
discretized versions of the other constraints of the OCP, including the state
and control constraints (II-B), the algebraic constraints (20), and the
boundary conditions. The resulting NLP problem has been modeled using Pyomo
[38], an open-source Python-based software package for modeling complex
optimization problems, and solved using the Interior Point OPTimizer (IPOPT)
solver, an open-source software package suitable for large-scale nonlinear
optimization. It implements an interior-point line-search filter method and
can be used to solve general NLP problems [39].

To solve smooth OCP, PSM are more accurate than traditional collocation
methods [40]. However, as they are based on global polynomials with
predetermined location of the nodes, their use in solving SOCP can give rise
to numerical difficulties since the location of the switches may not coincide
with the preallocated nodes. Hence, to solve these difficulties, the concept
of pseudospectral knots was introduced in [27]. The knots are some special
nodes which enable the global time interval to be divided into several time
subintervals to apply PSM over each subinterval, allowing the use of PSM to
solve SOCP, such as formation mission design problems. Therefore, a
pseudospectral knotting method has been used to solve the OCP resulting from
the EOCP.

###  VI-A The Pseudospectral Method

PSM are a type of collocation methods which use global polynomials to
parameterize the state variables. The nodes are generally obtained from a
Gaussian quadrature, which is used to collocate the differential-algebraic
equations. The use of global polynomials, instead of using piecewise-
continuous polynomials as interpolant between prescribed subintervals, makes
PSM more efficient and simpler than other direct methods [41].

Due to orthogonality properties, the first step to introduce PSM is to
transform the flight time interval of each flight $p$, defined as $t_{p}\in
I_{p}=[t_{I_{p}},t_{F_{p}}],\;p\in\\{1,\dots,N_{a}\\}$, into a new set of time
variables, $\tau_{p}\in[-1,+1],\;p\in\\{1,\dots,N_{a}\\}$. The transformation
that relates the original time interval of each aircraft, $t_{p}$, is the
following

$t_{p}=\dfrac{\left(t_{F_{p}}-t_{I_{p}}\right)\tau_{p}\;+\;\left(t_{F_{p}}+t_{I_{p}}\right)}{2}.$

In trajectory optimization problems, the three most commonly used sets of
collocation points are Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR), and
Legendre-Gauss-Lobatto (LGL) points. Notice that, although in all these
methods the discretization points used in the NLP problem must lie in the
$[-1,+1]$ interval, this is not required for the collocation points. Indeed,
depending on the set of collocation points chosen, both, one or none endpoints
of the interval $[-1,+1]$ can be excluded in the approximation of the control
variables. In this paper, a PSM based on collocation at Flipped Radau points
have been used [28]. The Flipped Radau collocation points lie in the interval
$(-1,+1]$.

Additionally, the EOCP stated to solve the formation mission design problem
has some peculiarities which led to the use of the knotting method in the PSM
that will be outlined in the next section.

###  VI-B The Pseudospectral Knotting Method

The obtained EOCP is a nonsmooth problem with switches in the state and
control variables, and in the dynamic constraints. In this EOCP, the initial
and final times are supposed to be different for each flight and,
consequently, the nodes of each flight do not match. However, to enforce the
separation constraints between aircraft, it is necessary to have the same
nodes for all the flights involved, at least, during the formation flight. The
pseudospectral knotting method solves this necessity in the following way.
They generalize the spectral patching method by the exchange of information
across the patches between two standard PSM in the form of switching
conditions [27]. These switching conditions are localized at the so-called
knots of the problem. Hence, the knots represent double nodes, coinciding with
the last node of one subinterval or patch and with the first node of the
following subinterval. Different kind of knots can be defined, such as soft or
hard knots, or free or fixed knots.

Following [27], in the formation mission design problems solved in this paper,
only two interior soft knots are considered, independently of the number of
aircraft involved in the formation. Besides, initial-time and final-time
conditions are considered as part of the general framework of knots. In
particular, the knots corresponding with the initial-time and final-time are
defined as hard knots because they are intrinsic to the problem formulation.

The interior knots are defined as free knots but constrained by the condition
that the time position of the first interior knot, $t_{k_{1}}$, should be
greater than the initial times of all flights, $\text{max}(t_{I,p})$, and the
time position of the second one, $t_{k_{2}}$, should be smaller than the final
times of all flights, $\text{min}(t_{F,p})$, as schematically shown in Fig. 5.

Fig. 5: Schematic flight times representation in a three-aircraft mission
design problem.

On this basis, three different time intervals, $I^{1}$, $I^{2}$, and $I^{3}$,
have been considered for each flight $p\in\\{1,\dots,N_{a}\\}$:

  * •

$I^{1}$: before the first interior knot, $t_{p}\in[t_{I,p},\;t_{k_{1}})$,
formation is not allowed.

  * •

$I^{2}$: between the two interior knots, $t_{p}\in[t_{k_{1}},\;t_{k_{2}}]$,
formation is allowed.

  * •

$I^{3}$: after the second interior knot, $t_{p}\in(t_{k_{2}},\;t_{F,p}]$,
formation is not allowed.

In Fig. 6 a schematic representation of the three different time intervals,
$I^{1}$, $I^{2}$, and $I^{3}$, and the knots and nodes of the problem are
represented. The vertical lines correspond to the times of each node, and the
big points represent the knots.

Fig. 6: Knots and nodes. Figure adapted from [27].

It is easy to see that, in two-aircraft formation design problems, the
formation can only happen if both aircraft are flying, in other words, during
the overlapped flight time [19]. Thus, the formation is only allowed during
the second time interval, $I^{2}$, and the collocation points of all the
aircraft are forced to be the same in this interval. Notice that, as mentioned
before, the two interior knots of the problem should be greater than
$\text{max}(t_{I,p})$ and lower than $\text{min}(t_{F,p})$, but are not fixed
to be equal to those values. Hence, the time value of the interior knots
should also be determined in the solution. The same assumption has been taken
for three-aircraft formations.

In addition to the numerical difficulties associated to the nonsmoothness of
the mission design problem, the Gibbs phenomenon should also be mentioned.
This phenomenon appears when a nonsmooth function is approximated with some
smooth functions [27], and the pseudospectral knotting method helps to deal
with it.

##  VII Numerical Results

To show the effectiveness of the methodology for the formation mission design
problem described in the previous sections, the following numerical
experiments have been carried out:

  * •

Experiment A: two-aircraft transoceanic mission design with one flight's
departure time free.

  * •

Experiment B: two-aircraft transoceanic mission design with delays in the
departure times.

  * •

Experiment C: three-aircraft transoceanic mission design with different fuel
savings schemes.

All the experiments involve transoceanic eastbound flights. Wind data from the
ERA-INTERIM reanalysis database has been used. As mentioned above, in this
paper only the cruise phase has been modeled, the rest of the flight phases
are neglected. Thus, the initial and final cruise phase locations of the
cruise phase of the flights have been 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 have been considered for each
flight. The initial mass of each aircraft is assigned as well as the initial
and final velocities, which have been set at typical cruise values for the
selected aircraft models. The initial heading angle is set to the initial
heading angle of the orthodromic path between the initial and final locations
of each flight. Arrival times are left free and a constraint on the maximum
temporal deviation from the scheduled arrival time of 45 minutes has been
introduced for each flight. The time and fuel burn weighting parameters,
$\alpha_{t}$ and $\alpha_{f}$, in the objective functional (6) have been set
to 0.3 and 0.7, respectively [35].

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 initial
guesses of the solution. The initial guesses for the latitude, the longitude
and the heading angle have been generated using the orthodromic between
departure and arrival locations of each flight. Typical cruise velocity and
fuel consumption of the aircraft model selected have been used to generate the
initial guesses of the velocity and the mass of each aircraft during the
flight, respectively.

###  VII-A Experiment A: two-aircraft transoceanic mission design with one
flight's departure time free

Experiment A involves two transoceanic eastbound flights, Flight 1 and Flight
2, with given fuel savings for the trailing aircraft and boundary values of
the state variables. Flight 1 and Flight 2 are operated by Aircraft 1 and
Aircraft 2, respectively. In the SOCP stated to solve the mission design
problem, the departure time of Flight 1 is fixed, whereas the departure time
of Flight 2 is left free. In case of formation, Aircraft 1 will be the
trailing one and Aircraft 2 will be the leader. Solving this problem entails
to decide which mode of flight, i.e. formation or solo flight, is optimal, the
optimal trajectories of the aircraft and the optimal departure time of Flight
2. Additionally, in case of formation flight, the rendezvous and splitting
locations and times must also be determined. To check the optimality of the
obtained solution, a comparison between formation flight and solo flight
results has been carried out.

The two transoceanic flights considered in this experiment have the following
departure and arrival locations:

  * •

Flight 1: New York (JFK) - Madrid (MAD).

  * •

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

The fuel burn savings for the trailing aircraft is set to 10%. The departure
time of Flight 1 has been set to 10:15. The boundary conditions for the state
variables of each aircraft are listed in Table I.

TABLE I: Experiment A: boundary conditions. Symbol | Units | Flight 1 | Flight 2
---|---|---|---
$\phi_{I}$ | [deg] | 40.64 | 45.47
$\phi_{F}$ | [deg] | 40.48 | 51.47
$\lambda_{I}$ | [deg] | -73.78 | -73.74
$\lambda_{F}$ | [deg] | -3.57 | -0.45
$\chi_{I}$ | [deg] | 66.51 | 55.70
$V_{I}$ | [m/s] | 240 | 240
$V_{F}$ | [m/s] | 220 | 220
$m_{I}$ | [kg] | $220\,000$ | $215\,000$

As mentioned above, the formation configuration is selected in advance. In
case of formation flight, Aircraft 1, associated to Flight 1, is forced to be
the trailing aircraft. Hence, Aircraft 2, associated to Flight 2, will be the
leader one.

In the optimal solution, formation flight has been selected as the optimal
solution and the obtained optimal departure time for Flight 2 has been 10:50,
35 minutes after the departure of Flight 1. The obtained optimal formation
flight routes are represented in Fig. 7, together with the solo flight routes
and the wind field. In this figure, dashed lines are used to plot solo flight
routes and solid lines are used to represent formation flight routes. Green
and blue lines represent Flight 1 and Flight 2, respectively.

Flight 1 starts 35 minutes before Flight 2, and due to this difference in the
departure times, the optimal solution implies a large initial detour of Flight
1 with respect to solo flight, in order to join Flight 2 as soon as possible,
thereby, the rendezvous point is very close to the departure location of
Flight 2.

Fig. 7: Experiment A: solo and formation flight routes.

(a) Mass (b) Heading angle (c) Mach number

Fig. 8: Experiment A: state variables.

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

Fig. 9: Experiment A: control variables.

The corresponding state and control variables are represented in Fig. 8 and
Fig. 9, respectively. For a better understanding of these figures, the
portions of the plots of the variables that correspond to formation flight
have been represented on a grey background. It can be observed that the
heading angle and the Mach number of both aircraft are the same during the
formation, as opposed to the mass flow rate which is different for each
aircraft. All the state variables vary smoothly, except the heading angle and
the Mach number at the beginning and at the end of the flight. This behavior
of the state variables is typical in boundary value problems, in which a quick
maneuver is necessary to steer each aircraft from its initial state to its
optimal route and from the optimal route to the final state. On the other
hand, the control variables vary in a quite smooth way except around the
switching instants. Notice that, the Gibbs phenomenon described in Section VI
is present near the rendezvous and the splitting locations. However, this non-
smooth behavior at the switching times lasts for a very short period of time
compared with the total flight time.

TABLE II: Experiment A: solo and formation flight results for Flight 1.

| Flight | Fuel | Covered | DOC
---|---|---|---|---
|  Time [$\boldsymbol{h}$] |  Burn [$\boldsymbol{kg}$] |  Distance [$\boldsymbol{km}$] | [ $\boldsymbol{mu}$]
Formation F. | 7.68 | 46543.93 | 6716.44 | 40875.15
Solo F. | 7.47 | 48596.82 | 6205.81 | 42085.37

TABLE III: Experiment A: solo and formation flight results for Flight 2.

| Flight | Fuel | Covered | DOC
---|---|---|---|---
|  Time [$\boldsymbol{h}$] |  Burn [$\boldsymbol{kg}$] |  Distance [$\boldsymbol{km}$] | [ $\boldsymbol{mu}$]
Formation F. | 6.19 | 40130.77 | 5348.58 | 34776.74
Solo F. | 6.18 | 39683.14 | 5462.84 | 34452.60

In this experiment, the great circle distance between departure and arrival
locations of Flight 1 is 5761.08 km and the great circle distance between
departure and arrival locations of Flight 2 is 5214.72 km. For these flights,
the great circle distance between departure locations is 537.00 km and the
great circle distance between arrival locations is 1244.77 km. Thus, this
scenario and the selected fuel burn saving for the trailing aircraft can be
considered not optimistic for formation flight.

Tables II and III summarize the results obtained in formation and solo flights
for Flight 1 and Flight 2, respectively. As expected, the flight times
increase for both aircraft in formation flight with respect to solo flight,
although for Flight 2 this flight time increment is negligible. Moreover, the
fuel consumption of Aircraft 2, the leader aircraft, increases whereas the
fuel consumption of Aircraft 1, the trailing one, decreases. The DOC is
expressed in generic monetary units, mu. The total DOC of formation flight is
almost one thousand generic monetary units lower than that of the solo flight,
which amounts to more than one percent of reduction. For Flight 1, the total
flight distance covered in formation flight is more than 500 km larger than in
solo flight. Surprisingly, for Flight 2, the total flight distance covered in
formation flight is smaller than the distance covered in solo flight.
Obviously, both of them are higher than the great circle distance. The reason
behind this unexpected result for Flight 2 is the Jet Stream. Indeed, in order
to get the greatest benefits from the wind field, Aircraft 2 in solo flight
takes a route more than 200 km longer than the orthodromic route. In case of
formation, this large detour, which would involve both aircraft, becomes less
advantageous and therefore it does not appear in the optimal solution.

###  VII-B Experiment B: two-aircraft transoceanic mission design with delays
in the departure times

The flights considered in Experiment B, Flight 1 and Flight 2, have the same
departure and arrival locations, the same boundary values for the state
variables of each aircraft, and the same fuel burn savings for the trailing
aircraft as in Experiment A. As before, in case of formation flight, Aircraft
1, associated to Flight 1, is forced to be the trailing aircraft and Aircraft
2, associated to Flight 2, is the leader one. The difference with respect to
Experiment A is that several delays in the departure times of both flights
have been considered.

This experiment has been carried out to determine how delays in the departure
times of both flights affect the formation flight, in terms of routes,
rendezvous and splitting locations and times, until formation flight becomes
non-optimal and, as a consequence, solo flights are selected by the algorithm.

TABLE IV: Experiment B: delays in the departure time with respect to the baseline scenario for the different cases. | Delays [min]
---|---
Cases | Flight 1 | Flight 2
$B_{1}$ | 20 | 0
$B_{2}$ | 15 | 0
$B_{3}$ | 10 | 0
$B_{4}$ | 5 | 0
$B_{5}$ | - | -
$B_{6}$ | 0 | 5
$B_{7}$ | 0 | 10
$B_{8}$ | 0 | 15
$B_{9}$ | 0 | 20
$B_{10}$ | 0 | 25
Fig. 10: Experiment B: routes obtained in the different cases considered.

For this analysis, several numerical simulations have been conducted
introducing delays in the departure times of both flights in Experiment A,
namely, 10:15 for Flight 1 and 10:50 for Flight 2, which represents the
baseline case. In particular, five-minute increments in the delays have been
considered for each flight. The actual values of the delays are listed in
Table IV, in which each case is identified by a different symbol
$B_{1},\dots,B_{10}$, being $B_{5}$ the baseline case.

For a better understanding of the results, the optimal aircraft routes
obtained for different delays in the departure times are represented together
in Fig. 10, in which the green and blue lines represent the optimal routes of
Flight 1 and Flight 2, respectively. The optimal route obtained in the
baseline case, $B_{5}$, is represented with a thicker line. Notice that, in
this figure, the scale is not the same for latitude and longitude axes, and
the symbols on the third axis denote the different cases.

In particular, Fig. 10 shows the rendezvous and splitting locations for each
case, which are represented with small black stars and triangles,
respectively. It can be observed that there is a great dependency of the
location of the rendezvous points on the delay. For instance, in cases $B_{6}$
and $B_{7}$, the rendezvous location is very close to the departure location
of Flight 2. On the contrary, the splitting locations barely change with
delay. When the delay on one departure time increases, since it takes more
time for the delayed aircraft to reach the other one, the rendezvous point
changes. However, once the formation is created, both aircraft follow a very
similar route for all the cases regardless of the delay in the departure time,
and therefore, the optimal splitting locations do not change significantly. In
Table V, the rendezvous and splitting locations and the distance covered
flying in formation obtained in cases $B_{2},\dots,B_{9}$ are displayed.
Notice the proximity of the splitting locations in all the cases and the
closeness among the rendezvous locations in cases $B_{6}$ and $B_{7}$ and the
departure location of Flight 2.

On the contrary, the individual flight times and the duration of the formation
flight, as well as the rendezvous and splitting times significantly change due
to the initial detour made by both aircraft to create the formation. In Table
VI, the rendezvous and splitting times and the duration of the formation
flight for each case are listed. Notice that in cases $B_{1}$ and $B_{10}$
solo flight is the optimal option. Therefore, the corresponding results have
not been reported in tables and figures, where only the results obtained in
formation flight are given. Further information on the results of cases
$B_{1}$ and $B_{10}$ can be found in the results of Experiment A.

TABLE V: Experiment B: rendezvous and splitting locations and formation flight
distance.

| Rendezvous | Splitting | FF
---|---|---|---
|  (Lat, Lon) [$\boldsymbol{deg}$] |  (Lat, Lon) [$\boldsymbol{deg}$] |  Distance [$\boldsymbol{km}$]
$B_{2}$ | (51.26, -65.69) | (57.05, -17.13) | 3208.86
$B_{3}$ | (49.15, -68.85) | (57.04, -17.11) | 3536.29
$B_{4}$ | (47.10,-71.22) | (57.04, -17.06) | 3829.61
$B_{5}$ | (46.04, -72.86) | (57.05, -16.83) | 4024.45
$B_{6}$ | (45.48, -73.71) | (57.06, -16.71) | 4130.80
$B_{7}$ | (45.49, -73.72) | (57.06, -16.71) | 4130.68
$B_{8}$ | (47.59, -71.91) | (57.06, -16.69) | 3823.66
$B_{9}$ | (49.92, -69.09) | (57.06, -16.73) | 3523.89

TABLE VI: Experiment B: rendezvous and splitting times and formation flight duration. |  Time [$\boldsymbol{h}$]
---|---
| Rendezvous | Splitting | FF duration
$B_{2}$ | 1.58 | 5.21 | 3.63
$B_{3}$ | 1.22 | 5.21 | 3.99
$B_{4}$ | 0.90 | 5.22 | 4.33
$B_{5}$ | 0.70 | 5.24 | 4.54
$B_{6}$ | 0.67 | 5.32 | 4.66
$B_{7}$ | 0.75 | 5.40 | 4.66
$B_{8}$ | 1.15 | 5.45 | 4.30
$B_{9}$ | 1.61 | 5.58 | 3.96
TABLE VII: Experiment B: flight times, fuel burn and DOC variations compared to solo flights. |  $\boldsymbol{\Delta}$ Time [$\boldsymbol{\%}$] |  $\boldsymbol{\Delta}$ Fuel burn [$\boldsymbol{\%}$] |
---|---|---|---
| Flight 1 | Flight 2 | Flight 1 | Flight 2 |  $\boldsymbol{\Delta}$ DOC [$\boldsymbol{\%}$]
$B_{2}$ | +2.56 | \+ 3.96 | -3.30 | \+ 1.97 | -0.14
$B_{3}$ | +2.56 | \+ 2.60 | -3.74 | \+ 1.80 | -0.51
$B_{4}$ | +2.57 | \+ 1.26 | -4.13 | \+ 1.55 | -0.89
$B_{5}$ | +2.79 | \+ 0.05 | -4.22 | \+ 1.13 | -1.17
$B_{6}$ | +3.79 | \+ 0.14 | -4.15 | \+ 1.01 | -1.09
$B_{7}$ | +4.94 | \+ 0.13 | -3.74 | \+ 1.01 | -0.79
$B_{8}$ | +6.07 | \+ 0.08 | -3.06 | \+ 1.01 | -0.36
$B_{9}$ | +7.25 | \+ 0.07 | -2.54 | \+ 1.00 | -0.02

The main results in terms of the objective functional are reported in Table
VII. In particular, the variations in the flight time and in the fuel
consumption for each flight and the reduction in the DOC obtained in cases
$B_{2},\dots,B_{9}$, all of them compared to solo flights results, are
reported. It can be observed that, in general, the more the flight time of
Flight 1 increases, the less the flight time of Flight 2 does and vice versa.
Besides, for those cases in which Flight 2 is delayed, the flight time
increment of Flight 2 is very small. About the fuel burn, it can be seen that
the greatest benefits for the trailing aircraft are achieved in the baseline
case, which is the optimal one, as established in Experiment A. Notice the
similarity among increments in fuel burn for Flight 2 in the results of cases
$B_{6},B_{7},B_{8}$, and $B_{9}$, in which Flight 2 has a delay in the
departure time. The greatest DOC reduction is 1.17% and the smallest is 0.02%.
This tiny DOC reduction implies that flying in formation in case $B_{9}$ has
negligible benefits, which could be further reduced by any contingency during
the flight.

In general, one may conclude that the greatest flight distances covered in
formation flight correspond to the highest total benefits in the fuel burn
reduction. However, comparing the results displayed in tables V and VII, it is
easy to check that longer distances covered in formation flight and greater
duration of formation flight do not always imply more benefits in terms of
reduction of DOC. It can be seen in tables V and VI, that in the results of
cases $B_{6}$ and $B_{7}$, aircraft are flying 4.66 hours and about 4130 km in
formation, being the greatest duration of formation flights and the greatest
distance covered in formation flights, respectively. Nevertheless, as it can
be observed in Table VII, the corresponding reductions in the DOC with respect
to solo flights are 1.09% and 0.79%, respectively. The reduction of the DOC is
slightly smaller than the one obtained in the baseline case, which amounts to
1.17%.

TABLE VIII: Experiment B: extra distance covered compared to solo flight distance. |  $\boldsymbol{\Delta}$ Distance [$\boldsymbol{km}$]
---|---
| Flight 1 | Flight 2
$B_{2}$ | 451.21 | -123.32
$B_{3}$ | 455.03 | -127.23
$B_{4}$ | 467.05 | -123.25
$B_{5}$ | 510.62 | -114.26
$B_{6}$ | 543.56 | -108.73
$B_{7}$ | 544.79 | -111.76
$B_{8}$ | 516.48 | -110.18
$B_{9}$ | 507.60 | -109.73

Table VIII has been added to give information about the total detour done by
each flight. As in Experiment A, in all the cases considered in Experiment B
there is a reduction in the distance covered by Flight 2 in formation flight
compared to solo flight. Noteworthy is the great detour made by Flight 1 in
the formation flight solution: about 500 kilometers of diversion from the solo
flight route.

The average computational time to find the solution has been $4.234$ s. For
the sake of comparison, the same problem has been solved using a multiphase
method [24]. In this case, two different OCP have been solved, namely, one in
which aircraft flight solo and another one in which they are forced to fly in
formation. The average computational time has been $7.784$ s.

###  VII-C Experiment C: three-aircraft transoceanic mission design with
different fuel savings schemes

Experiment C involves three transoceanic eastbound flights, Flight 1, Flight 2
and Flight 3, with given fuel savings scheme and boundary values of the state
variables. Flight 1, Flight 2, and Flight 3 are operated by Aircraft 1,
Aircraft 2, and Aircraft 3, respectively. The three flights considered in this
experiment have the following departure and arrival locations, the first two
of them being the same as in Experiment A:

  * •

Flight 1: New York (JFK) - Madrid (MAD).

  * •

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

  * •

Flight 3: Boston (BOS) - Paris (CDG).

The departure times of Flights 1,2, and 3, are set to 10:15, 10:50, and 10:30,
respectively. The first two of them are the same as in Experiment A. The
boundary conditions for the state variables of the three aircraft are given in
Table IX. The boundary conditions for the state variables of Aircraft 1 and
Aircraft 2 are the same as in Experiment A.

TABLE IX: Experiment C: boundary conditions for the three flights. Symbol | Units | Flight 1 | Flight 2 | Flight 3
---|---|---|---|---
$\phi_{I}$ | [deg] | 40.64 | 45.47 | 42.36
$\phi_{F}$ | [deg] | 40.48 | 51.47 | 48.85
$\lambda_{I}$ | [deg] | -73.78 | -73.74 | -71.06
$\lambda_{F}$ | [deg] | -3.57 | -0.45 | 2.35
$\chi_{I}$ | [deg] | 66.51 | 55.70 | 56.46
$V_{I}$ | [m/s] | 240 | 240 | 240
$V_{F}$ | [m/s] | 220 | 220 | 220
$m_{I}$ | [kg] | $220\,000$ | $215\,000$ | $210\,000$

Table X summarizes the results obtained assuming that each flight is performed
as a solo flight. It can be observed that the flight times, the fuel
consumption, the covered distance, and the DOC of each flight are quite
different.

TABLE X: Experiment C: results for the three solo flights. | Flight 1 | Flight 2 | Flight 3
---|---|---|---
Flight Time [$\boldsymbol{h}$] | 7.47 | 6.19 | 6.90
Fuel burn [$\boldsymbol{kg}$] | 48596.82 | 39683.14 | 42877.40
Covered Distance [$\boldsymbol{km}$] | 6205.81 | 5462.84 | 5945.89
DOC [$\boldsymbol{mu}$] | 42087.53 | 34462.32 | 37466.18

An analysis is performed to determine how the fuel savings scheme affects the
formation flight, in terms of routes, rendezvous and splitting locations and
times, and flight times.

As already mentioned, the only formation configuration allowed is the in-line
formation and the relative position in the formation are fixed. The considered
benefits for the intermediate and the trailing aircraft are ranging from 6% to
14%. Several numerical simulations have been conducted introducing changes in
the fuel savings for the intermediate and trailing aircraft, in which each
case is identified by a different symbol $C_{1},\dots,C_{5}$. The fuel savings
considered in the different cases for the intermediate and trailing aircraft
have been listed in Table XI.

TABLE XI: Experiment C: fuel savings for the intermediate and the trailing aircraft, in the different cases. | Fuel savings [%]
---|---
$C_{1}$ | 6
$C_{2}$ | 8
$C_{3}$ | 10
$C_{4}$ | 12
$C_{5}$ | 14

It can be seen that in all the considered cases there are five discrete states
of the switched dynamical system that represent the joint behavior of the
aircraft. In Fig. 11 a schematic representation of these five discrete states,
from State I to State V, is given. In particular, each discrete state is
characterized by the following modes

  * •

State I: all the aircraft fly in solo mode.

  * •

State II: two aircraft fly in formation and one in solo mode.

  * •

State III: all the aircraft fly in formation.

  * •

State IV: two aircraft fly in formation and one in solo mode.

  * •

State V: all the aircraft fly in solo mode.

Fig. 11: Experiment C: discrete states representation for three-aircraft
formation.

It is easy to see that in the switches among these discrete states, there are
two rendezvous points, RV1 and RV2, where RV1 is the first one and RV2 is the
second one, and two splitting points, SP1 and SP2 , where SP1 is the first one
and SP2 is the second one. In Fig. 11, the points RV1, RV2, SP1, and SP2 are
represented with grey points. Small triangles have been used to represent the
departure and destination of each flight.

Fig. 12 and Fig. 13 represent the three-aircraft routes obtained in the cases
$C_{1},\dots,C_{5}$. In these figures, the behavior of the rendezvous and
splitting locations, respectively, is detailed. It can be observed that the
location of both rendezvous points, RV1 and RV2, and the location of the
second splitting point, SP2, do not present significant variations in the
different cases. On the contrary, the location of the first splitting point,
SP1, notably changes, differing in some cases in more than 1000 km. This
variation in the location of SP1, which is the point in which there is a
switch between State III and State IV, implies that the three-aircraft
formation time and distance largely depend on the fuel savings scheme.
Consequently, the formation time and distance in State IV also have a high
dependency on it. Results in tables XII and XIII confirm this conclusion. It
can also be observed in these tables that, in all the considered cases, the
formation times and distances corresponding to discrete State II are small
compared to other discrete states.

TABLE XII: Experiment C: formation distances for the two- and three-aircraft formation phases. | Formation distances [$\boldsymbol{km}$]
---|---
Cases | State II | State III | State IV
$C_{1}$ | 113.53 | 3008.99 | 1371.24
$C_{2}$ | 38.19 | 3323.62 | 1142.36
$C_{3}$ | 41.01 | 3743.05 | 737.69
$C_{4}$ | 44.48 | 3847.74 | 595.32
$C_{5}$ | 41.44 | 4229.41 | 257.29
TABLE XIII: Experiment C: formation times for the two- and three-aircraft formation phases. | Formation times [$\boldsymbol{h}$]
---|---
Cases | State II | State III | State IV
$C_{1}$ | 0.13 | 3.30 | 1.72
$C_{2}$ | 0.06 | 3.69 | 1.43
$C_{3}$ | 0.05 | 4.26 | 0.92
$C_{4}$ | 0.04 | 4.42 | 0.84
$C_{5}$ | 0.06 | 4.92 | 0.74
TABLE XIV: Experiment C: total flight times, fuel burn and DOC reduction compared to solo flights, in the different cases. |  $\boldsymbol{\Delta}$ Time [$\boldsymbol{\%}$] |  $\boldsymbol{\Delta}$ Fuel burn [$\boldsymbol{\%}$] |
---|---|---|---
| Flight 1 | Flight 2 | Flight 3 | Flight 1 | Flight 2 | Flight 3 |  $\boldsymbol{\Delta}$ DOC [$\boldsymbol{\%}$]
$C_{1}$ | 2.28 | 0.13 | 0.30 | -1.83 | -4.17 | 0.90 | -1.14
$C_{2}$ | 2.45 | -0.02 | 0.16 | -2.87 | -5.81 | 1.10 | -1.80
$C_{3}$ | 2.89 | -0.06 | 0.10 | -4.15 | -7.39 | 1.25 | -2.51
$C_{4}$ | 2.90 | -0.19 | 0.01 | -5.29 | -9.14 | 1.41 | -3.25
$C_{5}$ | 3.54 | -0.10 | 0.02 | -6.69 | -10.64 | 1.45 | -3.97
Fig. 12: Experiment C: routes obtained in the different cases considered.
Detail of the rendezvous locations. Fig. 13: Experiment C: routes obtained in
the different cases considered. Detail of the splitting locations.

The main results related to the objective function obtained in Experiment C
are reported in Table XIV. In particular, the time increment of each flight,
the reduction or increase of the fuel burn and the reduction in the DOC with
respect to the results obtained assuming each flight is performed as a solo
flight, are reported. It can be seen in this table that, for case $C_{5}$, the
reduction in the DOC is almost 4% comparing to the case in which formation is
not allowed.

The average computational time to find the solution has been $15.274$ s. For
the sake of comparison, the same problem has been solved using a multiphase
method [24]. In this case, 13 different OCP have been solved, since, in this
case, there are 13 possible flight phase sequencing options, as shown in [24].
The average computational time has been $90.948$ s. This shows that the method
presented in this paper is able to drastically reduce the computational time
to solve formation mission design problems.

##  VIII Conclusions

In this paper, a novel framework to solve the formation mission design problem
is presented. In the proposed framework, the mission is modeled as a switched
dynamical system, in which aircraft are assumed to have two flight modes,
namely solo and formation flight, and the discrete state of the switched
dynamical system is the result of their combination. The discrete dynamics is
modeled using logical constraints in disjunctive form. The formation mission
design problem is transformed into an optimal control problem for a switched
dynamical system with logical constraints, which is solved using the embedding
approach.

The embedding approach is a unifying technique able to efficiently tackle both
the switching dynamics and the logical constrains of the optimal control
problem, transforming it into a smooth optimal control problem, which has been
solved using a knotting pseudospectral method.

The main advantages of this approach are that the multi-phase formulation is
avoided, as well as the use of binary variables, decreasing the computational
time and effort in finding the solution. This approach is substantially
different from the previous approaches, which are based on exhaustively
analyzing every possible formation mission individually and then, comparing
the results. Additionally, this approach is easily scalable to design
formation mission problems of an arbitrary number of aircraft.

Several numerical experiments have been conducted with two and three
transoceanic flights. The results demonstrate that the proposed framework for
mission design is fast and accurate and that the formation flight has great
potential to reduce fuel consumption and emissions and, therefore, to mitigate
the environmental impact of air transport sector. An analysis of the solutions
has also been carried out in order to study how the delays in the departure
time and changes in the fuel saving scheme influence the formation mission.
The results indicate that delays and changes in the fuel savings scheme have
significant influence on the formation mission.

This suggests that a deeper understanding of the influence of uncertainties in
these and other parameters, such as the initial mass of the aircraft or
meteorological conditions, on the formation mission is needed to assess the
feasibility of formation flight in realistic scenarios usually characterized
by in the presence of uncertainties. This will be subject of future research.

##  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] D. S. Lee, D. W. Fahey, P. M. Forster, P. J. Newton, R. C. Wit, L. L. Lim, B. Owen, and R. Sausen, “Aviation and global climate change in the 21st century,” _Atmospheric Environment_ , vol. 43, no. 22-23, pp. 3520–3537, 2009.
  * [2] The Intergovernmental Panel on Climate Change, “Global Warming of 1.5C,” ICPP, IPCC, Tech. Rep., 2018.
  * [3] International Air Transport Association. (2016) Press Release No. 59. IATA forecasts passenger demand to double over 20 years. [Online]. Available: https://www.iata.org/en/pressroom/pr/2016-10-18-02/
  * [4] V. Madonna, P. Giangrande, and M. Galea, “Electrical power generation in aircraft: Review, challenges, and opportunities,” _IEEE Transactions on Transportation Electrification_ , vol. 4, no. 3, pp. 646–659, 2018.
  * [5] International Civil Aviation Organization, “ICAO Environmental Report 2016. Aviation and Climate Change,” ICAO, Tech. Rep., 2016.
  * [6] ——, “Convention on International Civil Aviation, Annex 2,” ICAO, Tech. Rep., 1944.
  * [7] M. A. Camilleri, “Aircraft operating costs and profitability,” in _Travel Marketing, Tourism Economics and the Airline Product_. Springer, 2018, ch. 12, pp. 191–204.
  * [8] A. Ning, T. C. Flanzer, and I. M. Kroo, “Aerodynamic performance of extended formation flight,” _Journal of Aircraft_ , vol. 48, no. 3, pp. 855–865, 2011\.
  * [9] P. B. S. Lissaman and C. A. Shollenberger, “Formation flight of birds,” _Science_ , vol. 168, no. 3934, pp. 1003–1005, 1970.
  * [10] J. E. Kless, M. J. Aftosmis, S. A. Ning, and M. Nemec, “Inviscid analysis of extended-formation flight,” _AIAA Journal_ , vol. 51, no. 7, pp. 1703–1715, 2013.
  * [11] J. P. Slotnick, “Computational aerodynamic analysis for the formation flight for aerodynamic benefit program,” in _Proceedings of the 52nd Aerospace Sciences Meeting, AIAA SciTech Forum_ , National Harbor, MD, USA, 2014.
  * [12] S. R. Bieniawski, R. W. Clark, S. E. Rosenzweig, and W. E. Blake, “Summary of flight testing and results for the formation flight for aerodynamic benefit program,” in _Proceedings of the 52nd Aerospace Sciences Meeting, AIAA SciTech Forum_ , National Harbor, MD, USA, 2014.
  * [13] D. J. Haalas, S. R. Bieniawski, B. Whitehead, T. Flanzer, and W. B. Blake, “Formation flight for aerodynamic benefit simulation development and validation,” in _Proceedings of the 52nd Aerospace Sciences Meeting, AIAA SciTech Forum_ , National Harbor, MD, USA, 2014.
  * [14] T. C. Flanzer, S. R. Bieniawski, and W. B. Blake, “Operational analysis for the formation flight for aerodynamic benefit program,” in _Proceedings of the 52nd Aerospace Sciences Meeting, AIAA SciTech Forum_ , National Harbor, MD, USA, 2014.
  * [15] M. J. Vachon, R. Ray, K. Walsh, and K. Ennix, “F/A-18 aircraft performance benefits measured during the autonomous formation flight project,” in _Proccedings of the AIAA Atmospheric Flight Mechanics Conference and Exhibit_ , Monterey, CA, USA, 2002.
  * [16] C. Liu, B. Jiang, R. J. Patton, and K. Zhang, “Integrated fault-tolerant control for close formation flight,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 56, no. 2, pp. 839–852, 2019.
  * [17] J. Seo, Y. Kim, S. Kim, and A. Tsourdos, “Collision avoidance strategies for unmanned aerial vehicles in formation flight,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 53, no. 6, pp. 2718–2734, 2017.
  * [18] G. J. Durango, C. Lawson, and A. Z. Shahneh, “Formation flight investigation for highly efficient future civil transport aircraft,” _The Aeronautical Journal_ , vol. 120, no. 1229, pp. 1081–1100, 2016.
  * [19] J. Xu, S. A. Ning, G. Bower, and I. Kroo, “Aircraft route optimization for formation flight,” _Journal of Aircraft_ , vol. 51, no. 2, pp. 490–501, 2014\.
  * [20] T. E. Kent and A. G. Richards, “Analytic approach to optimal routing for commercial formation flight,” _Journal of Guidance, Control, and Dynamics_ , vol. 38, no. 10, pp. 1872–1884, 2015.
  * [21] W. B. Blake and T. C. Flanzer, “Optimal routing for drag-reducing formation flight: A restricted case,” _Journal of Guidance, Control, and Dynamics_ , vol. 39, no. 1, pp. 173–176, 2016.
  * [22] G. C. Bower, T. C. Flanzer, and I. M. Kroo, “Formation geometries and route optimization for commercial formation flight,” in _Proceedings of the 27th AIAA Applied Aerodynamics Conference_ , San Antonio, TX, USA, 2009.
  * [23] S. Hartjes, M. E. van Hellenberg Hubar, and H. G. Visser, “Multiple-phase trajectory optimization for formation flight in civil aviation,” in _Advances in Aerospace Guidance, Navigation and Control_. Springer, 2018, pp. 389–405.
  * [24] S. Hartjes, H. G. Visser, and M. E. van Hellenberg Hubar, “Trajectory optimization of extended formation flights for commercial aviation,” _Aerospace_ , vol. 6, no. 9, p. 100, 2019.
  * [25] S. C. Bengea and R. A. DeCarlo, “Optimal control of switching systems,” _Automatica_ , vol. 41, pp. 11–27, 2005.
  * [26] S. Bengea, K. Uthaichana, M. Zěfran, and R. A. DeCarlo, “Optimal control of switching systems via embedding into continuous optimal control problem,” in _The Control Handbook: Advanced Methods_ , W. S. Levine, Ed. CRC Press, 2011, ch. 31, pp. 31–1 – 31–23.
  * [27] I. M. Ross and F. Fahroo, “Pseudospectral knotting methods for solving optimal control problems,” _Journal of Guidance, Control, and Dynamics_ , vol. 27, no. 3, pp. 397–405, 2004.
  * [28] D. Garg, M. A. Patterson, W. Hager, A. V. Rao, D. A. Benson, and G. T. Huntington, “An overview of three pseudospectral methods for the numerical solution of optimal control problems,” _Advances in the Astronautical Sciences_ , vol. 135, pp. 1–17, 2009.
  * [29] D. B. Seenivasan, A. Olivares, and E. Staffetti, “Multiaircraft optimal 4D trajectory planning using logical constraints,” _International Journal of Aerospace Engineering_ , vol. 2019, pp. 1–15, 2019.
  * [30] D. G. Hull, _Fundamentals of Airplane Flight Mechanics_. Springer, 2007.
  * [31] Eurocontrol, “User manual for the Base of Aircraft Data (BADA),” Eurocontrol, Tech. Rep. 130416, 2013.
  * [32] 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.
  * [33] M. D. Buhmann, _Radial Basis Functions: Theory and Implementations_. Cambridge University Press, 2003.
  * [34] S. A. Ning, “Aircraft drag reduction through extended formation flight,” Ph.D. dissertation, Stanford University, 2011.
  * [35] 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.
  * [36] S. Wei, M. Zěfran, and R. A. DeCarlo, “Optimal control of robotic systems with logical constraints: Application to UAV path planning,” in _Proceedings of the 2008 IEEE International Conference on Robotics and Automation_ , Pasadena, CA, USA, 2008.
  * [37] T. M. Cavalier, P. M. Pardalos, and A. L. Soyster, “Modeling and integer programming techniques applied to propositional calculus,” _Computers and Operations Research_ , vol. 17, no. 6, pp. 561–570, 1990.
  * [38] 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.
  * [39] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical programming_ , vol. 106, no. 1, pp. 25–57, 2006.
  * [40] J. T. Betts, _Practical Methods for Optimal Control and Estimation using Nonlinear Programming_. SIAM, 2010.
  * [41] F. Fahroo and I. M. Ross, “A spectral patching method for direct trajectory optimization,” _Journal of the Astronautical Sciences_ , vol. 48, no. 2, pp. 269–286, 2000.

|  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.
---|---