The aviation sector represents a major challenge to the environment, being one of the top ten emitters of greenhouse gases (GHG) in the world, and accounting for approximately 4% of human contribution to the global warming [Klower et al., 2021, Lee et al., 2009]. According to a recent report by the intergovernmental panel on climate change (IPCC) [IPCC, 2018], 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 transport sector’s carbon dioxide emissions, far less than other modes of transport like light-duty vehicles (30%) and heavy-duty vehicles (36%), it seems to be harder to decarbonize. The reason for this lies in the predicted growth rate for air traffic, which is projected to be greater than for other modes of transport.

The predicted growth rate for air traffic is another concerning aspect. Before the COVID-19 crisis, the compound annual growth rate (CAGR) forecast over the next two decades, assessed by the international air transport association (IATA), estimated that the number of flights was expected to continue increasing by 3.5% annually, which translates to twice the current number of passengers by 2040 [IATA, 2016]. It is well known that the COVID-19 pandemic has severely affected the whole aviation sector, international commercial flights in particular, and, consequently, air traffic has plummeted by more than two thirds compared with previous levels.¹ At present, vaccination rates, travel restrictions, and quarantine procedures are rendering it unclear how quickly the aviation sector will return to normal values. However, air traffic growth is expected to recover to pre-crisis levels within the next few years, while continuing to increase annually [Lai et al., 2022]. This continuous growth in air traffic demand imposes the search for new solutions to mitigate its potential negative effects on our planet.

Unfortunately, the measures taken by the aviation sector until recently are considered to be insufficient to mitigate its impact on climate change. Current policies should be strengthened and a firm commitment from all the stakeholders is essential to success. For instance, the IATA is looking for improvements in different areas in order to achieve a reduction in net aviation CO₂ emissions of 50% by 2050, relative to 2005 levels [IATA, 2013]. The different measures considered can be classified into technical improvements, including more efficient aircraft and engines, operational measures, including more efficient flight procedures and airport operations, and the search of alternative, more sustainable fuels [Abrantes et al., 2021]. Fig. 1.1 gives a schematic overview of the potential contribution of the mentioned measures to the reduction of the net CO₂ emissions of international aviation.

As such, the predicted effects of the aviation sector on climate change are greater than desired [IPCC, 2018, ICAO, 2016], and a great effort is required to mitigate them. In this regard, formation flight has the potential to make a significant contribution. Furthermore, formation flight offers great promise not only with regard to mitigating the environmental impact of aviation but also in terms of increasing air traffic management (ATM) capacity and reducing air traffic control workload, which represents another major challenge of the sector.

Formation flight is where two or more aircraft fly at a short distance from each other. The notion of formation flight arose from the observation of nature, in particular from the different species of bird that fly in formation during their migrations [Lissaman and Shollenberger, 1970]. Aircraft formation flight was first introduced in military aviation for mutual defense during missions, concentrating the fire power during attacks, and, consequently, improving combat efficiency. 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. The reduction achieved strongly depends on the longitudinal distance between aircraft, this being negligible over approximately 40 wingspans. Formations in which aircraft fly at a longitudinal distance between 10 and 40 wingspans are called extended formations [Durango et al., 2016]. They represent a trade-off between safety and efficiency. This thesis will investigate the possibility of implementing extended formations in commercial aviation. In particular, the focus is on the possibility of implementing economically and environmentally beneficial formation flights in the presence of different sources of uncertainty, namely the departure times of the flights and the fuel savings during formation flight.

The key enabling factors for this concept of operation consist of a formation control system to keep the aircraft flying in formation at the optimal relative position to optimize the fuel burn savings for the trailing aircraft and an ATM system capable of synchronizing flights departing from different airports to ensure that the formation mission occurs as planned. Within a formation flight, instabilities such as meandering and external factors may affect the motion of the vortices. It is therefore important to maintain the relative position between the follower aircraft and the leader’s wake vortices precisely, as the fuel burn savings are very sensitive to that relative positioning. This can be done using a formation control system capable of continuously locating the wake vortices, maintaining the aircraft in the optimal relative position, and, in this way, optimizing the fuel savings during the formation flight [Caprace et al., 2019]. While technical issues related to maintaining an efficient and safe formation flight have been solved in recent years, some concerns have yet to be addressed. In particular, a major change in airworthiness standards, policies, and procedures is required. At the 40^th international civil aviation organization (ICAO) Assembly, formation flight was proposed as a strategic objective, and the necessity of developing a new operational concept including reduced separations between aircraft to allow formation flights was established [Abeyratne, 2020].

Due to its potential for reducing fuel consumption and GHG emissions, formation flight has recently attracted increasing interest, and the studies on the aerodynamics of formation flight have given way to studies on formation flight planning. The majority of the commercial formation flight planning studies are based on geometric models which do not allow use of either accurate dynamic models or meteorological forecasts. The past few years have brought with them some research studies [Hartjes et al., 2018, Hartjes et al., 2019] in which the mission design problem has been studied for two- and three-aircraft formations using a multiphase optimal control approach. This methodology has clear advantages compared to previous ones since the optimal control approach allows the use of accurate models, improving the predictability of the trajectories. However, using the multiphase approach implies having a fixed switching structure in advance, that is to say, knowing the transition among phases established in advance. Since the switching structure of the optimal solution is unknown for the formation mission design problem, any combination should be addressed and, then, the obtained solutions must be compared to establish which one corresponds to the minimum operation-related aircraft costs, also known as direct operating costs (DOC) [Camilleri, 2018]. Thus, this approach is feasible if the number of flights considered is low and both the type of formation and the relative position of each aircraft in the formation are fixed. Otherwise, the number of combinations quickly increases with the number of flights, making the previous approach hardly scalable.

A more general framework is therefore needed, which is able to design formation missions solving a single optimal control problem (OCP) in which neither the number of phases nor the switching structure is known in advance. Moreover, this approach should be easily scalable to design, for instance, formation missions involving more than three flights. For this purpose, the formation mission design problem can be formulated as an OCP of a switched dynamical system. These systems are described by both continuous and 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 within a formation, and their combination is represented by the discrete state of the switched dynamical system, which models their joint dynamic behavior. Each flight mode will be represented by different dynamical equations, which may or may not include formation flight benefits in terms of fuel savings. Additionally, logical constraints in disjunctive form, based on the streamwise distance between aircraft, model the switching logic among the discrete states of the system. To solve the proposed problem, the embedding technique [Bengea et al., 2011] has been employed to obtain a new smooth formulation of the problem in which the switching dynamics are defined without binary variables together with equality and inequality constraints. The obtained problem can be solved using classical optimal control techniques.

In addition to using accurate models, the operational concept of the formation flight entails a deeper understanding of how uncertainties affect the trajectories of each aircraft involved in the formation. In particular, addressing the uncertainty quantification in formation flight design problems is essential to improving the predictability of the trajectories and achieving more realistic solutions and estimations of the costs and the formation benefits. There are different sources of uncertainties which can affect the whole ATM system. One of the main ones is uncertainty in flight departure times, which causes trajectory uncertainty and ATM-context inefficiencies [Rivas and Vazquez, 2016]. Timing is also a crucial factor in formation missions. Effectively, considering the usual cruise speed of most long-haul commercial aircraft, missing the rendezvous location by, for instance, ten minutes means spatially missing the partner aircraft by 150 km. Such cases require catch-up maneuvers, which result in a loss of performance compared to the planned formation mission. As such, uncertainties in departure times have been considered in the formation mission design problem, along with uncertainties in the fuel savings during formation flight. These uncertain parameters have been modeled as random variables described by probability density functions.

The stochastic formation flight mission design problem can be formulated as follows. Given several commercial flights together with the relevant weather forecast, the problem consists in establishing how to organize them in formation, of two- or three-aircraft, or solo flights and in finding the trajectories that minimize the expected DOC value of the formation mission. To solve it, the formation mission design problem is formulated as an OCP of a stochastic switched dynamical system and solved using nonintrusive generalized polynomial chaos (gPC) based stochastic collocation [Matsuno et al., 2015, Li et al., 2014]. The gPC method converts the stochastic switched OCP into an augmented deterministic switched OCP, which can be solved using the embedding approach together with classical numerical optimal control techniques. This technique allows statistical and global sensitivity analysis of the stochastic solutions to be conducted at a low computational cost.

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

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

It is well known that uncertainty regarding the weather forecast, in particular the wind field, is of particular interest in the context of optimization trajectories. However, the stochastic impact of the wind field in the formation mission design problem has not been considered. Being a random function of space and time, a wind field should be modeled as a stochastic process. Although generalized polynomial chaos can be employed to represent stochastic processes, specific techniques are needed for their representation, such as the Karhunen—Loeve expansion. As a consequence, the dimensionality of the problem increases, as does the complexity of dealing with it. Therefore, it is assumed that the aircraft model is accurately known, as is the corresponding weather forecast. Neither storms nor operational uncertainties have been considered. The study of the formation mission design problem in the presence of stochastic processes such as uncertain wind field, is left for future research. In this thesis, only random variables have been considered.

Recently, the aviation industry has shown renewed interest in achieving the real implementation of commercial formation flight. Indeed, during the completion of this thesis, the Airbus Fello’fly project [Airbus, b] has been carried out, in which Airbus collaborated with airlines, air navigation service providers, and civil aviation authorities to tackle the challenges of formation flight and demonstrate its operational feasibility. The first flight test in the field of transoceanic flights was conducted in November 2021.³ This test, in which a system for formation flight control proved itself in an operational environment, shows that this technology has reached the highest readiness level, according to the European Union’s scale [Héder, 2017].

The general objective of this thesis is to solve the formation mission design problem for commercial aircraft in the presence of uncertainties.

The stochastic formation mission planner should include accurate dynamic models of the aircraft involved in the formation mission, the relevant flight data such as departure and arrival times and locations, wind forecast, a model of the effects of flying in formation on the fuel consumption of the trailing aircraft, the probability density functions that characterize the uncertain parameters of the problem, and an objective functional reflecting the DOC of the formation mission.

The stochastic formation mission planner presented in this thesis, which is based on stochastic optimal control, finds the optimal solution for up to three aircraft candidates to fly in formation, establishing how to arrange them in formation or solo and determining the stochastic optimal trajectories. The stochastic solution of the formation mission planning problem includes the expected values and standard deviations of the latitude, longitude, and timing of the optimal trajectories. The solution also includes the expected values and standard deviations of the latitude, longitude, and timing of the rendezvous and splitting points as well as the expected values and standard deviations of the arrival times of the flights. Finally, by performing the sensitivity analysis of the solutions, the relative influence of the uncertain parameters of the problem on the variability of the stochastic solution components can be estimated. The schematic representation of the elements of the formation mission planner is given in Fig. 1.3.

The general objective of this thesis has been divided into several specific objectives:

The solutions obtained demonstrate that formation flight is environmentally and economically beneficial in the presence of uncertainties regarding the departure times of the aircraft currently observed at major international airports and the fuel savings achieved during formation flight.

The formation mission planner presented in this thesis is intended to be incorporated into a ground-based system 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 who are responsible for compiling and modifying the flight plan. They also monitor the aircraft during flight. The computational time of the proposed algorithm makes it suitable for the strategic and tactical planning phases of the flight, as well as for the operational phase, in other words, during flight.

The contributions of this thesis to the state of the art in formation mission planning can be summarized as follows:

The rest of the document is organized as follows. In Chapter 2, an overview of formation flight is given. In Chapter 3, the model of the system that represents the formation flight is introduced. In Chapter 4, the optimal control method employed to solve the deterministic formation mission design problem is introduced. In Chapter 5, a general overview to numerical methods for optimal control is given. In Chapter 6, the stochastic optimal control method employed to solve the stochastic formation mission design problem is introduced. In Chapter 7, the methodology introduced in Chapter 4 is employed to solve three different formation mission design problems in the absence of uncertainties. In Chapter 8, the methodology introduced in Chapter 6 is employed to solve two different formation mission design problems in the presence of uncertainties. Finally, the conclusions and a discussion of the future work are outlined in Chapter 9.

As stated in Chapter 1, the current GHG emissions caused by civil aviation have become a priority at international level. Formation flight has the potential not only to mitigate the environmental impact of aviation but also to increase the capacity of the ATM. According to [ICAO, 1944], “the formation operates as a single aircraft with regard to navigation and position reporting,” meaning that two or three aircraft flying in formation would be treated as a single aircraft for air traffic control purposes. This would help to reduce their workload. Furthermore, formation flight may even provide enhanced inter-aircraft communication, since the proximity between aircraft enables cost reductions and limits jamming sensitivity [Klavins, 2004].

Therefore, the question that naturally arises is what needs to be changed in the current paradigm to ensure effective implementation of formation flight in commercial aviation. Surprisingly, making today’s aircraft appropriate to fly in formation does not require significant changes to the current aircraft appearance. However, the possibility of implementing formation flight in commercial aviation depends upon several factors, such as defining a new concept of operation for formation flight of commercial aircraft, changing the international airworthiness regulation, updating the ATM system in such a way that formation flights and solo flights can share the same airspace, creating alliances among airlines for sharing the benefit derived from formation flight, and developing the avionics systems for formation control. These complex issues will be discussed in the following sections.

Formation flight is based on the fact that any aircraft moving through a fluid generates lift by imparting a swirling motion to the air. In particular, an aircraft flying in the cruise phase produces a vortex wake which mainly maintains its structure and intensity for several kilometers. This wake is characterized by a downwash region in the center of the wake and two similar upwash regions outside of the vortex core. A schematic representation of these upwash and downwash regions is given in Fig. 2.3. A second aircraft flying in either of these upwash regions, will increase its apparent angle of attack obtaining useful energy. The increment of the apparent angle of attack produces an effective forward rotation of the lift vector, L→, resulting in a lower lift-dependent component of drag, D_i. Thus, a significant decrease in induced drag is achieved [Ning et al., 2011]. For a better understanding, the apparent and induced angle of attack, a_e and a_i respectively, the local and induced velocities, u_i and w_i respectively, and the lift and induced drag vector for both solo and formation flight are given in Fig. 2.4, where sub-index S F has been used to indicate “Solo Flight” and sub-index F F indicates “Formation Flight”. Additionally, the components related to S F are represented in blue, whereas the components related to FF are represented in green. The total angle of attack, a, together with the free flowing fluid velocity, μ_ꝏ, have been represented in black. Notice that the total angle of attack is defined by the following sum of the apparent and the induced angle of attack,

From Eq. (2.1), it is easy to see that an increase of the apparent angle of attack is obtained due to the formation flight, which results in a decrease of the induced angle of attack. Consequently, both the induced velocity and the induced drag are reduced.

Formations can be classified based on the longitudinal distance between aircraft in close and extended formations or based on their relative positions in V, inverted-V, and echelon formations.

Formation flight has great potential to contribute significantly to reducing fuel consumption and the environmental impact of the air transport sector while increasing its capacity. These advantages can be achieved by updating aircraft avionics and air traffic management technologies and procedures. Its potentialities have been widely demonstrated from the aerodynamic point of view via numerical simulations and flight tests. Flight safety aspects have been also addressed, leading to the proposal of extended formations for commercial formation flight, which are intrinsically safer than close formations. Although formation flight is still not a reality for commercial aviation, international aviation organizations, aircraft manufacturers, airlines, air traffic service providers, and researchers are collaborating to make formation flight for commercial aircraft possible in the near future. The challenging Fello’fly Airbus project is a great indication of this.

In aircraft trajectory optimization problems developed within the ATM context, a three degrees of freedom point-mass dynamic model of the aircraft is usually employed. In this model, the aircraft is represented as a point with variable mass due to the fuel consumption, moving in a four-dimensional space characterized by the three spatial dimensions together with time.

In this section, first, the reference systems are introduced. Later, some usual assumptions made to build an aircraft model to be used in the ATM context are described. After that, some simplifications are introduced to the model. Finally, the flight envelope is presented.

The differential algebraic equations (DAE) system (3.3) is applicable to any aircraft in any phase of solo flight. However, depending on the aircraft type, the aerodynamic coefficients, the wingspan, and the reference surface, among others, the parameters of this DAE system change. In formation flight, further modifications must be introduced into the dynamic equations of the trailing aircraft.

As stated in Chapter 2, extended formations are intrinsically safer than close formations, which is why only extended formations have been considered in this thesis. 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 thesis, 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 during formation flight is the following:

where b is the wingspan of the leader aircraft and D_pq is the great-circle distance between aircraft p and q at time t.

The haversine formula has been used to determine the great-circle distance D between two points given their latitudes and longitudes:

where (ϕ_q,λ_q) and (ϕ_p,λ_p) are the latitude and longitude of aircraft q and p, respectively, at any time.

As mentioned in Chapter 2, due to the formation flight, a significant decrease in induced drag is achieved for the trailing aircraft. According to Assumption L, the expression of the parabolic drag polar is

where C_D0 is the parasitic drag, which represents the aerodynamics clean performance of the aircraft, and C_Di is the induced drag, that is, the drag component due to lift. The induced drag can be expressed by

where AR = b²/S is the aspect ratio, e is an efficiency factor, which will be equal to the unity for elliptical lift distribution and lower than the unity for any other distribution, and K = 1/(π • AR • e) is the lift-induced drag coefficient.

As mentioned in Chapter 2, due to formation flight, a significant decrease in the induced drag, KC²_L, is achieved for the trailing aircraft, hence, the parabolic drag polar for a trailing aircraft flying in a formation can be expressed as [Hartjes et al., 2018], [Hartjes et al., 2019]

where ϵ represents the rate of induced drag reduction achieved for the trailing aircraft. The value of ϵ is highly dependent on the type of aircraft involved in the formation, the streamwise and the lateral distance between aircraft, the type of formation flight, and the wind, among others. According to [Ning, 2011], in which a zero wind field has been assumed, a maximum reduction in the induced drag of 30±3% is achieved in a two-aircraft formation and a maximum reduction in the induced drag of 40 ± 6% in a three-aircraft formation (95% confidence intervals).

Therefore, during formation flight, the dynamics model of the leader aircraft is represented by the ODE system (3.3) with the standard parabolic drag polar, whereas the dynamic model of the trailing aircraft is represented by the same set of equations but replacing the standard parabolic drag polar by Eq. (3.11). The selection of the form of the drag polar is made based on the longitudinal distance between aircraft, as schematically illustrated in Table 3.1.

Considering that aircraft will not join in formation during the take-off, landing, and approach phases of flight [Hartjes et al., 2018], in the mission design problem studied in this thesis, only the cruise phase of the flights is considered. Moreover, the motion of the aircraft involved in the formation flight is restricted to the horizontal plane at cruise flight level, not allowing changes in the flight level. Thus, in addition to the assumptions of Section 3.1, the following assumption is made:

Assumption M: The flight path angle is zero, γ = 0.

Considering Assumption M, the dynamic models of the aircraft can be simplified. The resulting simplified model is a two degrees of freedom point variable-mass dynamic obtained by substituting γ = 0 in the DAE system (3.3). The motion of the aircraft restricted to the horizontal cruise plane is described by the following set of kinematics and dynamics DAE:

together with the following algebraic equation:

where the state vector has only five components, two less than DAE system (3.3): the two dimensional position variables, latitude and longitude, denoted by φ and λ, respectively; the heading angle χ; the true airspeed V; and the mass of the aircraft m.

In this set of equations, the control vector formally includes three components: the thrust force T, the lift coefficient C_L, and the bank angle μ. However, two control variables are directly related due to the requirement of equilibrium of forces given by Eq. (3.13) and, therefore, the number of control variables is two. Thus, for aircraft p, the state vector is x_p = (ϕ_p,λ_p,V_p,χ_p ,m_p), ∀ p ϵ {1 ,…,N_a}, and the control vector is u_p = (T_p,μ_p(C_Lp)), ∀p ϵ {1,…, N_a}, with N_a the number of aircraft involved in the mission design problem.

Note that, for problem conditioning and numerical stability, normalized versions of both the DAE system (3.12) and Eq. (3.13) have been used in this thesis in the numerical experiments described in chapters 7 and 8.

Concerning the formation flight modeling, the following assumption has been adopted:

Considering Assumption N, the fuel savings factor can be directly introduced in the last differential equation of the DAE system (3.12) instead of considering the reduction in the induced drag. With R_fuel being the rate of reduction in fuel burn consumption due to the formation flight, the mass flow rate equation for the trailing aircraft during the formation flight can be expressed as

A schematic description of the differences between the flight models employed in this thesis to describe the solo and formation flights is given in 3.2. The model is selected based on the longitudinal distance between aircraft.

It is quite difficult to model the drag reduction ϵ and, hence, the fuel burn reduction ℛ_fuel achieved during a formation flight, since there are many factors 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.

It is important to point out that the full three degrees of freedom point variable-mass dynamic model described in the DAE system (3.3) with the induced drag reduction modeling described by Eq. (3.11) could have been employed in the formation mission planner presented in this thesis. However, using the full model makes the statement of the problem, its numerical resolution, and the analysis of the results more complicated. For this reason, the reduced model has been used in this thesis.

Taking into account the assumptions made in this section, the obtained results are slightly less accurate. In particular, taking Assumption M and Assumption N into consideration may lead to a loss of accuracy in the speed profile and flight altitude. In [Ning, 2011], the aerodynamic performance of extended formations is analyzed in detail. The results show that aircraft within a formation fly at a slower cruise speed than in their respective solo flights, as this helps to achieve a reduction in the penalties associated with viscosity and compressibility effects induced by the formation flight. Consequently, aircraft within the formation would also need to fly at a lower altitude, at a higher lift coefficient, or a combination of both compared with solo flight. Therefore, aircraft in formation tend to cruise at a slightly lower speed than when flying solo, as shown in [Voskuijl, 2017], [Hartjes et al., 2018], and [Hartjes et al., 2019].

Other assumptions made in the formation mission design problem presented in this thesis regard the relative position of each aircraft, which is set in advance, and the type of formation allowed in three-aircraft formations, which is the in-line formation.

Finally, a maximum detour time from the solo flight time has been set in order to perform a more realistic formulation of the formation mission design problem, according to the usual operational requirements.

Consider a dynamical system described by a differential equation

in which variable t represents the time variable and the functions x(t) and u(t) represent the state and the control input of the system, respectively. Let

be a cost functional, which maps the time variable and the state and control functions into a single number.

An instance of a smooth OCP can be informally stated as follows: determine the control input that steers the system from an initial state to a final state minimizing the cost functional J. In general, the time variable and the control and state functions are subject to constraints.

This simple instance of an OCP corresponds to a single-phase OCP. In general, optimal control models include more than one phase. Each phase in an OCP represents a portion of the trajectory of a dynamical system, and may be subject to different equations of motion, different time discretization, different control parameterizations, or different path constraints.

A different set of ODEs can be employed in each phase to reflect the switching behaviour of a dynamical system. A specific time discretization fitted to the dynamics of that portion of the trajectory can be used in every phase. For instance, if one part of a trajectory has slower dynamics and another part has faster dynamics, the time interval can be split into two phases with coarser time discretization for the first phase and a finer time discretization for the second phase. In the control parameterization method, the control space is discretized by approximating the control function using a linear combination of basis functions. In this way, the OCP is reduced to a NLP with a finite number of decision variables. Different control parameterizations such as piecewise constant and piecewise linear, can be employed in the same problem by defining different phases. Different path constraints, which are bound constraints on some performance parameter, can also be imposed within each phase. Phases are also necessary to impose intermediate constraints upon some variables. This can be done by imposing them as boundary constraints at the junction between two phases. Finally, linkage constraints should be enforced in order to assemble all the phases of the problem to get a unique trajectory for the system [Falck et al., 2021].

OCP can be formulated for continuous-time dynamical systems, as is the case in this thesis, and for discrete-time dynamical systems. OCP formulations may also include binary or integer variables to model decision-making processes.

As mentioned before, in general, optimal control problems include not only differential constraints but also algebraic constraints and path constraints. A more general formulation of a continuous smooth OCP is the following:

where u(t) ϵ ℝ^m, x(t) ϵ ℝⁿ, and Eq. (4.3a) represents the performance index of the OCP, in which J is a Bolza cost functional, consisting of two terms, the endpoint cost functional M, also known as Mayer term, which is defined on a neighborhood of B, and the running cost functional, L, also known as Lagrange term, which is a continuously differentiable functional. The control input u(t) is constrained to belong, at each time instant, to the bounded and convex set Ω. Eqs. (4.3b) and (4.3c) represent the set of algebraic-differential equations, where function f determines the dynamics of the system and function g describes the right-hand side of the algebraic constraints. The problem can be also constrained by path constraints represented by Eq. (4.3d), where η_l and η_u represent the lower and upper bound for each function η. Finally, Eqs. (4.3e) and (4.3f) represents the initial and final boundary conditions, respectively. The initial time t_I, the initial state x(t_I), the final time t_F, and final state x(t_F) are assumed to be restricted to a boundary set B as follows: (t_I,x (t_I) ,t_F, x (t_F)) ϵ B = T_I × B_I × T_F × B_F ⊂ ℝ²ⁿ +².

In typical aerospace optimization problems the performance index of Eq. (4.3a) only contains the Mayer term. This is the case of this thesis, in which the objective functional in Eq. (3.18) is in Mayer form and represents a combination of fuel consumption and flight time of the aircraft.

Eqs. (4.3b) and (4.3c) correspond to the DAE systems given in Table 3.2, which describe the motion of each aircraft involved in a formation. In particular, the DAE system in the left column of the table describes the solo flight and the flight of the leader of a formation, and the DAE system in the right column of the table describes the flight of the intermediate and trailing aircraft of a formation.

The path constraints represented by Eq. (4.3d) correspond to the flight envelope given in Eq. (3.4). These constraints are usually given in terms of lower and upper bounds for several state and control variables such as flight altitude, load factor and airspeed. Notice that, Eqs.(4.3c) and (4.3d) also include logical constraints based on the stream-wise distance between aircraft, which are used in this thesis to model the switching between solo flight mode and formation flight mode for each aircraft. These constraints will be discussed later, in Section 4.2.4 of this chapter.

Eqs. (4.3e) and (4.3f) represent the initial and final conditions which are specified for each experiment in chapters 7 and 8.

For the sake of ease of exposition, in the formulation of the OCP for switched dynamical systems presented in the following section, neither algebraic equations nor path constraints have been included.

In this thesis, the formation mission design problem is formulated as an OCP for a switched dynamical system. Switched 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.

In this section, following [Bengea et al., 2011], the formulation of a simple OCP given in Section 4.1, is generalized for switched systems, introducing the switched optimal control problem (SOCP) for two-switched dynamical systems. Hereinafter, subindex "S" has been employed in the formulation of the OCP for switched dynamical systems for the sake of clarity. In particular, the dynamical model of a two-switched dynamical system can be represented by

In this model, the continuously differentiable vector fields, f₀, f₁ : ℝ × ℝⁿ × ℝ^m → Rⁿ, specify the dynamics of the two possible modes of the system. The control input u_s(t) ϵ Ω ⊂ ℝ^m is constrained to belong, at each time instant t, to the bounded and convex set Q. The binary variable υ_s (t) is the mode selection variable that identifies which of the two possible system modes, f₀ or f₁, is active. Thus, both of them, u_s (t) and υ_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 xs(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)) ϵ B = T_I × B_I × T_f × B_F ⊂ R²ⁿ⁺². The performance index of the SOCP has the following form

where the endpoint cost function M_s is defined on a neighborhood of B and formally coincide with function M, and L₀ and L₁ 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

subject to dynamical equations and endpoint constraints from the set of Eqs. (4.4).

Note that, although the statements of the OCP and the SOCP are quite similar, they are very different in their essence. The fact that the binary variable, υ_s (t) ϵ {0, 1}, appears within the SOCP formulation, requires an additional transformation of this kind of problems to enable the use of classical optimal control techniques to solve it. To do so, the embedding approach is introduced in the following section.

Two different numerical approaches have been traditionally followed to solve OCPs, indirect and direct methods.

The indirect methods follow the strategy of “ first optimize, then discretize” and are based on the Pontryagin’s Maximum Principle. By contrast, direct methods follow the strategy of “first discretize, then optimize” and transcribe an infinite-dimensional OCP into a finitedimensional optimization problem, which, in general, is a nonlinear programming (NLP) problem [Betts, 2010].

Indirect methods are based on the necessary conditions of the Pontryagin’s maximum principle. However, these conditions depend on the constraints of the OCP and change when some constraints become active, i.e., when the state or the control variables reach the boundary of their admissible sets. The main drawback of this setting is that it is not possible to know in advance when this occurs, and therefore, it is not easy to derive a general numerical method from these conditions. Moreover, accurate initial guesses are required, not only for state and control variables, as it is in many numerical methods, but also for the costate variables, which are non-intuitive and therefore, very difficult to be accurately established. For the same reasons, analytical solutions can be obtained only for the simplest OCPs.

Direct methods do not rely on the necessary conditions of the Pontryagin’s maximum principle. They are numerical methods in which the process to obtain the optimal solution of the OCP is quite easier than in indirect methods, because although, in general, large scale NLP problems are obtained with a great amount of variables and constraints, fast and reliable NLP solvers are available to solve them. In Fig. 5.1, a schematic overview of both advantages and disadvantages of each approach is given.

Since the advantages of direct methods outweigh their disadvantages compared to indirect approaches, in this thesis, a direct approach has been employed to solve OCPs.

Direct methods can be divided in two main categories, direct shooting methods and direct collocation methods, as shown in Fig. 5.2.

Direct shooting methods can be further divided into two categories: single and multiple shooting methods. Single shooting methods are control parameterization methods in which a time grid is defined and the controls are discretized on that grid as a set of NLP variables, while the states are discretized using numerical integration. These methods are the simplest technique to numerically solve OCPs. The major advantage of them is that, although a high-dimensional state vector is involved in the problem, the transcribed problem will have few degrees of freedom. Additionally, only initial guesses for controls are required. However, the main disadvantage of single shooting is that this technique is quite sensitive to initial conditions.

These drawbacks make them poorly suitable for solving most OCPs [Betts, 2010]. To reduce the sensitivity to initial conditions of the single shooting method, the multiple shooting method has been devised, in which the time interval is split into shorter subintervals, and thereafter, the problem is transcribed as in the single shooting methods, where the initial state of each subinterval is a variable of the resulting NLP problem. In the multiple shooting method, although the sensitivity-related issues are fixed, the number of variables involved in the NLP problem increases.

Direct collocation methods do not rely only on parameterization of the control, they rely on parameterization of both state and control. Within collocation methods, there are two major groups: local collocation and global collocation methods. As its name suggests, in local collocation methods local approximations of states and controls are made, and local integration techniques are applied for the differential equations that represent the dynamic system. Local collocation methods are the simplest collocation methods and are less accurate than global collocation methods [Huntington and Rao, 2008]. Some of the best-known local collocation methods are the trapezoidal and the Hermite-Simpson collocation schemes. The main difference between them is the kind of functions that are employed to approximate the system dynamics in each of them, being piecewise linear functions in the trapezoidal scheme and piecewise quadratic functions in the Hermite-Simpson collocation scheme. Unlike local collocation methods, global ones use global polynomials to approximate the states across the whole time interval. The main family of global collocation methods are the pseudospectral methods (PSM). In this thesis, PSM have been used. They will be introduced in the next section.

PSM are a type of direct global collocation methods which use global polynomials to parameterize the state variables. The nodes are generally obtained from a Gaussian quadrature, which is employed to collocate the set of 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 [Fahroo and Ross, 2000].

The first step to introduce PSM is to transform the original time variable, t ϵ I = [t_I, t_F], into a new time variable, τ ϵ [−1, +1]. The transformation that relates both time variables is the following

From Eq. (5.1), it is easy to obtain the following relationship:

In trajectory optimization problems, the four most commonly used sets of collocation points are the Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR), also known as direct LGR, Flipped Legendre-Gauss-Radau (f-LGR), and Legendre-Gauss-Lobatto (LGL) points. Notice that, although in all these methods the discretization points employed in the problem must lie in the [−1, +1] interval, this is not required for the collocation points. Indeed, although the position of the points changes from one set of collocation points to another, as described in detail below, the major difference between the methods is that, depending on the set of collocation points chosen, both, one, or none endpoints of the interval [−1, +1] will be excluded in the approximation of the control variables [Garg et al., 2011].

Collocation points are roots of orthogonal polynomials, or linear combinations of them and its derivatives [Huntington, 2007, Garg et al., 2009]. In particular, the four sets of collocations points seen above are based on Legendre polynomials. Denoting the Legendre polynomial of k-th degree as P_K, where K represents the total number of discretization points, the particular relationship between the different sets of collocation points and the Legendre polynomial is specified in Fig. 5.3.

In Fig. 5.4 a schematic overview of the positioning of the points for the different collocation schemes is given. In particular, it can be seen that both endpoints are included in the LGL points, i.e. τ ϵ [−1, 1], neither of them are included in the LG points, i.e. τ ϵ (−1, 1) and just one of the endpoints is included in the LGR points. If the initial point −1 is included, i.e. τ ϵ [−1, 1), the set of collocation points are usually called just LGR points or direct LGR points, whereas, when the terminal endpoint +1 is included, i.e. τ ϵ (−1, 1], the resulting set of points is known as the f-LGR points. Notice that, while the LG and LGL points located symmetrically with respect to the origin, the direct LGR and f-LGR points are not.

Since in the LG points, the endpoints are not included, additional quadrature constraints are necessary in order to force collocation at initial and terminal endpoints. These additional constraints make the problem resolution unnecessarily more complicated. This drawback can be overcome using the LGL points. However, the use of the LGL points has some additional drawbacks such as the slow convergence of the costate variables. Finally, due to their asymmetric positioning with respect to the origin, both direct LGR and f-LGR collocation points, are good candidates, because they permit collocation at endpoints without significant convergence issues for the costate variables [Fahroo and Ross, 2008, Garg, 2011].

In [Garg et al., 2009], it is stated that direct LGR and f-LGR collocation points are useful in infinite-horizon problems and in finite-horizon problems in which endpoints constraints are involved or the state at the terminal time is included in the objective function.Therefore, PSM based on f-LGR points have been employed in this thesis, which will be introduced in the next section.

In this section, the specification of the PSM employed to solve the formation mission design problem studied in this thesis is given. The peculiarities of the problem that make it to be solved using the knotting PSM are discussed and the solver employed to solve the resulting NLP problem is described.

The deterministic formation mission design problem presented in Chapter 4 is formulated as a SOCP, which has been transcribed into an EOCP using the embedding approach. Finally, the knotting PSM is applied to solve the EOCP.

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

The interior knots are defined as free knots but constrained by the condition that the time position of the first interior knot, t_k1, should be greater than the initial times of all flights, max(t_I,p), and the time position of the second one, t_k2, should be smaller than the final times of all flights, min(t_F,_p), as schematically shown in Fig. 5.6. On this basis, three different time intervals, I¹, I², and I³, have been considered for each flight:

In Fig. 5.7 a schematic representation of the three different time intervals, I¹, I², and I³, and the knots and nodes of the problem are presented. The vertical lines correspond to the times of each node, and the big points represent the knots.

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 [Xu et al., 2014]. Thus, the formation is only allowed during the second time interval, 1², and the collocation points of all the aircraft are forced to be the same in this interval, I², 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 max(t_I,p) and lower than min(t_F,p), but are not fixed to be equal to those values. Thus, the time value of the interior knots should be also determined in the solution. The same assumption has been taken for three-aircraft formations.

In short, the resulting EOCP together with the inequality and equality constraints given in Section 4.2.4 can be solved using traditional optimal control techniques.

There are several integrated software packages available to solve OCPs using direct collocation methods such as DIDO [Ross, 2020], GPOPS [Rao et al., 2009] and GPOPS-II [Patterson and Rao, 2014]. All of them implement PSM and are user-friendly, in the sense that they allow people with no knowledge of the underlying pseudospectral theory to use them. However, they are commercial software packages for the commercial Matlab environment. To gain a deeper insight into knotting PSM and to avoid dependance from commercial software, the knotting pseudospectral method has been implemented from scratch in this thesis. The four different collocation points schemes have been implemented and tested, obtaining that f-LGR points show the best convergence for the formation mission design problem. Therefore, the Radau knotting PSM have been selected to transcribe the EOCP into a NLP problem.

The resulting NLP problem has been modeled using Pyomo [Hart et al., 2012], an opensource Python-based software package for modeling complex optimization problems, and solved using 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 employed to solve general NLP problems [Wachter and Biegler, 2006].

In this chapter, the main numerical methods for optimal control have been described. The accuracy of global collocation methods has led to choose in this thesis a pseudospectral method to solve optimal control problems, namely the pseudospectral method based on the Flipped Legendre-Gauss-Radau collocation points, which is especially suitable due to the presence of endpoints constraints in the optimal control problem and the state at the terminal time in the objective functional. The peculiarities of the formation mission design problem studied in this thesis, has made it necessary to use a special pseudospectral method, the knotting pseudospectral method. To gain a deeper insight into these methods and to avoid dependance from commercial software, the knotting pseudospectral method has been implemented from scratch in this thesis.

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

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

Another approach to model uncertainties in switched dynamical systems is to treat uncertainties as random variables or random processes and recast the original deterministic switched dynamical system as a stochastic switched dynamical system. This type of stochastic systems are different from those represented by classical stochastic differential equations, where the random inputs are idealized processes. Consider a stochastic optimal control problem in which only random variables are present in its formulation like in the formation mission design problem studied in this thesis. One of the most commonly used methods to solve OCPs of stochastic dynamical systems of this type is the Monte Carlo sampling [Shapiro, 2003], in which independent realizations of the random variables are generated based on their probability distributions. For each realization, the OCP becomes deterministic. After solving the deterministic instances of the problem that correspond to these realizations of the random variables, an ensemble of solutions is obtained, from which statistical information can be extracted. Although Monte Carlo sampling is straightforward to apply as it only requires iterative resolutions of deterministic instances of the OCP, it requires a large number of solutions, because the solution statistics converge relatively slowly. For example, the mean value typically converges as 1/N, where N is the number of realizations. The need for large number of solutions for accurate results can lead to an excessive computational cost, especially for OCPs that are already computationally intensive in the deterministic settings, as the mission design problem considered in this thesis. The Generalized Polynomial Chaos (gPC) expansion can contribute to alleviating this drawback.

A systematic and coherent presentation of numerical strategies for uncertainty quantification and stochastic computing is given in [Xiu, 2010], with a focus on the methods based on the gPC approach. Both the intrusive stochastic Galerkin and the nonintrusive stochastic collocation approaches to gPC are illustrated in detail. The nonintrusive stochastic collocation approach based on regression is described in [Sudret, 2008] and [Blatman and Sudret, 2010].

In this thesis, a stochastic collocation method has been selected. With this approach, a small number of sample points of the random variables are used to jointly solve particular instances of the SSOCP. The obtained solutions are then expressed as orthogonal polynomial expansions in terms of the random variables using these sample points. This is a nonintrusive methodology because the model equations are not altered. Depending on the distributions of the random variables, different types of orthogonal polynomials can be chosen to achieve better numerical precision. This technique allows statistical and sensitivity analysis of the stochastic solutions to be conducted at a low computational cost. Thus, the gPC method converts the SSOCP into an augmented deterministic SOCP, in which particular instances of the SSOCP are solved together as a single deterministic OCP. It is important to point out that the gPC method is applicable to solve the SSOCP when the solution depends smoothly on the random variables. This means that the solutions obtained for each combination of sample points of the random variables must give rise to the same discrete solution, i.e., to the same sequence of discrete states of the switched dynamical system that represents the formation mission.

The presence of random variables in the SOCP converts it into an SSOCP. In the formation mission design problem studied in this thesis, not all the elements of the SSOCP contain random variables. Specifically, the fuel consumption reduction factor ℜ_fuel in the mass flow rate equation (3.14) of the aircraft flying in formation as a trailing or an intermediate aircraft is a random variable. In contrast, the equations of motion associated with the state variables ϕ, λ, χ, and V do not contain random variables. The departure times of the aircraft are also assumed to be random variables of the SSOCP. One of these departure times coincides with t_I, the initial time of the formation mission. These random variables are assumed to be independent and characterized by probability density functions.

Let θ = (θ₁,θ₂,…,θ_N) be the vector of random variables that represent the random parameters of the SOCP. For sake of clarity, hereinafter the subindex “SS” will be employed in the formulation of the OCP for stochastic switched dynamical systems. It is assumed that the control vector u_ss(t) is deterministic and, therefore, it is a function of time only, whereas the state vector x_ss(t, θ) is stochastic and, therefore, it is a function of both time and the vector of random variables, θ [Li et al., 2014].

The continuous-time SSOCP has been formulated as in [Li et al., 2014], namely all the elements of the problem depend on θ and the dynamic equations, the path constraints, and the boundary conditions must be satisfied almost surely. In the formulation of the SSOCP, the dynamical equations and endpoint constraints are:

The performance index of the SSOCP has the following form:

The SSOCP is thus stated as follows

subject to dynamical equations and endpoint constraints from the set of Eqs. (6.1). It can be seen that the difference between the formulations of the deterministic SOCP and the SSOCP is that the vector random variables θ is included in the latter.

It is important to point out that probabilistic constraints, also known as chance constraints, could be introduced in the formulation of the SSOCP [Matsuno et al., 2015]. However, in this thesis, chance constraints have not been considered in the formulation of the SSOCP employed to solve the stochastic formation mission design problem.