CHAPTER 5. NUMERICAL OPTIMAL CONTROL METHODS

In this chapter, a general overview to numerical methods for optimal control is given. They can be classified in indirect and direct numerical methods based on whether they are based on optimality conditions or not. Direct methods can be classified into shooting and collocation methods. The latter, can be classified into local and global collocation methods. The main advantages and disadvantages of each class of methods is discussed.

5.1. Introduction: direct and indirect numerical optimal control methods

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.

5.2. Direct methods

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.

5.3. Pseudospectral methods

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

$$t=\frac{(t_F-t_I)\tau +(t_F+t_I)}{2}.\mathrm{ (5.1)}$$

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

$$dt=\frac{t_F-t_I}{2}d\tau .\mathrm{ (5.2)}$$

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.

5.3.1. Pseudospectral method using f-LGR points

In this section, following [Garg et al., 2011, Garg et al., 2009], the PSM using the f-LGR collocation points will be presented. Consider K flipped Radau collocation points lying on the interval τ ϵ (−1,1] and N discretization points lying on the interval τ ϵ [−1,1], which include the mentioned K collocation points together with the initial endpoint −1. In this case, K = N - 1. For a better understanding of their position, a schematic representation of the position of the f-LGR collocation points and the corresponding discretization points, is given in Fig. 5.5.

Thus, considering the time transformation defined in Eq. (5.1), each component j-th of the state variable evaluated at the k-th collocation point, xj(τ_k), can be approximated using a basis of N Lagrange polynomials L_i,i = 1,…, N, by the following expression:

$$x_j({\tau }_k)\approx {\displaystyle \sum _{i=1}^N{x}_j({\tau }_i)}{L}_i({\tau }_k),fork=1,\dots ,K,\mathrm{ (5.3)}$$

where the nodes τ_i,i = 1 ,…,N are the discretization points, τ_k,k = 1 ,…,K are the collocation points, and Li(τ_k) is a Lagrange polynomials basis particularized for each collocation point k, given by the Lagrange polynomial interpolation:

$$L_i(\tau )={\displaystyle \prod _{\begin{array}{l}j=1\\ j\ne 1\end{array}}^N\frac{\tau -{\tau }_j}{{\tau }_i-{\tau }_j}.}\mathrm{ (5.4)}$$

Differentiating Eq. (5.3), and defining a new matrix $D_{ki}={\dot{L}}_i({\tau }_k)$, yields

$${\dot{x}}_j({\tau }_k)\approx {\displaystyle \sum _{i=1}^N{x}_j({\tau }_i)}{\dot{L}}_i({\tau }_k)={\displaystyle \sum _{i=1}^N{D}_{ki}{x}_j({\tau }_i)},fork=1,\dots ,K.\mathrm{ (5.5)}$$

PSM are characterized by the differentiation matrix D ϵ ℝ^K×N, which is used to approximate the derivative of the state at the collocation points. The differentiation matrix for the f-LGR collocation points is defined as follows:

$$D_{ki}=\{\begin{array}{l}\frac{\dot{g}({\tau }_k)}{({\tau }_k-{\tau }_i)\dot{g}({\tau }_i)},\text{for }k\ne i,\\ \frac{\ddot{g}({\tau }_i)}{2\dot{g}({\tau }_i)},\text{for }k\ne i,\end{array}\mathrm{ (5.6)}$$

where the function g(τ_i) = (1 + τ_i)[P_K(τ_i) -P_K−1(τ_i)], based on the roots of the combination of Legendre polynomials for the f-LGR points, shown in Fig. 5.3.

Notice that, with the method explained here, increasing the degree of the Lagrange polynomials allows the state variables to be approximated at every discretization point, including both endpoints.

However, the control variables are only discretized and approximated at the f-LGR collocation points using N - 1 Lagrange polynomials. Therefore, the j-th component of the control can be approximated as follows:

$$u_j({\tau }_k)\approx {\displaystyle \sum _{i=2}^N{u}_j({\tau }_i)}{L}_i({\tau }_k),fork=1,\dots ,K.\mathrm{ (5.7)}$$

Using Eq. (5.2), the first time derivative of the j-th of the state variable, with respect to t ϵ I = [t_I,t_F], ${\dot{x}}_j(t)$, can be expressed as follows:

$${\dot{x}}_j(t)=\frac{2}{t_F-t_I}{\displaystyle \sum _{i=1}^N{D}_{ki}{x}_j({\tau }_i)},fork=1,\dots ,K.\mathrm{ (5.8)}$$

Thus, the differential equation (4.1) of the OCP, ẋ(t) = f (t, x(t), u(t)), will be incorporated in the NLP problem as:

$${\displaystyle \sum _{i=1}^N{D}_{ki}{x}_j({\tau }_i)}-\frac{t_F-t_I}{2}f({\tau }_k,x({\tau }_k),u({\tau }_k))=0.\mathrm{ (5.9)}$$

As mentioned before, global polynomials are used in PSM. Discontinuities in the time derivatives of the state or control variables can not be accurately approximated using this type of polynomials. The EOCP described in Section 4.2.2, which arises in solving the formation mission design problem, is a nonsmooth OCP with discontinuities in the time derivatives of the state and control variables caused by the switches in the dynamic constraints of the problem. Additionally, in the EOCP, the initial and final times of each flight are, in general, different and, consequently, the collocation nodes of each flight do not match. However, to model the logical constraints from Eq. (4.16) that describe the switching logic among discrete states of the switched dynamical system that represent the formation flight, it is necessary to have the same set nodes for all the flights involved, at least during the flight phases in which aircraft fly in formation. These reasons make PSM not very suitable in dealing with the non-smooth EOCP. Therefore, an extension of the PSM, called pseudospectral knotting method has been employed in this thesis, which is outlined in the next section.

5.3.2. Pseudospectral knotting method

Pseudospectral knotting methods generalize the spectral patching method by the exchange of information across the different patches among standard PSM in the form of switching conditions [Ross and Fahroo, 2004]. These constraints are localized at the so-called knots of the problem. 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. So, 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.

In relation 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 [Ross and Fahroo, 2004], and the pseudospectral knotting method allows this phenomenon to be prevented.

5.4. Solution of the Embedded Optimal Control Problem

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:

•I¹: before the first interior knot: t_p ϵ [ t_I,_p, t_k1), ∀ p ϵ {1 ,…,N_a}, formation is not allowed.

•I²: between the two interior knots: t_p ϵ [t_k1, t_k2], ∀p ϵ {1,…, N_a}, formation is allowed.

•I³: after the second interior knot: t_p ϵ (t_k2, t_F,p,], ∀p ϵ {1,…, N_a}, formation is not allowed.

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

5.5. Conclusions

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.

Bibliography

[Abeyratne, 2020] Abeyratne, R. (2020). The outcome of the 40th ICAO assembly: A new look at ICAO? Air and Space Law, 45(1).

[Abrantes et al., 2021] Abrantes, I., Ferreira, A. F., Silva, A., and Costa, M. (2021). Sustainable aviation fuels and imminent technologies-CO₂ emissions evolution towards 2050. Journal of Cleaner Production, 313:127937.

[Airbus, a] Airbus. https://www.airbus.com/en/newsroom/press-releases/2021-11-airbus-and-its-partners-demonstrate-how-sharing-the-skies-can-save. Accessed: 2021-11-18.

[Airbus, b] Airbus. https://www.airbus.com/en/innovation/disruptive-concepts/biomimicry/fellofly. Accessed: 2022-01-10.

[Alwan and Liu, 2018] Alwan, M. S. and Liu, X. (2018). Theory of Hybrid Systems: Deterministic and Stochastic. Springer.

[Ansari et al., 2020] Ansari, M. Y., Ahmad, A., Khan, S. S., Bhushan, G., et al. (2020). Spatiotemporal clustering: A review. Artificial Intelligence Review, 53(4):2381–2423.

[Assembly 40th, 2019] Assembly 40th, I. https://www.icao.int/Meetings/a40/Documents/WP/wp_317_en.pdf. Accessed: 2021-11-21.

[Aung and Tan, 2010] Aung, H. H. and Tan, K.-L. (2010). Discovery of evolving convoys. In Gertz, M. and Ludäscher, B., editors, International Conference on Scientific and Statistical Database Management, pages 196–213. Springer.

[Bengea et al., 2011] Bengea, S., Uthaichana, K., Zěfran, M., and DeCarlo, R. A. (2011). Optimal control of switching systems via embedding into continuous optimal control problem. In Levine, W. S., editor, The Control Handbook: Advanced Methods, chapter 31, pages 31–1 – 31–23. CRC Press.

[Bengea and DeCarlo, 2005] Bengea, S. C. and DeCarlo, R. A. (2005). Optimal control of switching systems. Automatica, 41:11–27.

[Betts, 2010] Betts, J. T. (2010). Practical Methods for Optimal Control and Estimation using Nonlinear Programming. SIAM.

[Bieniawski et al., 2014] Bieniawski, S. R., Clark, R. W., Rosenzweig, S. E., and Blake, W. E. (2014). Summary of flight testing and results for the formation flight for aerodynamic benefit program. In Proceedings of the 52nd Aerospace Sciences Meeting, National Harbor, MD, USA.

[Blake and Flanzer, 2016] Blake, W. B. and Flanzer, T. C. (2016). Optimal routing for dragreducing formation flight: A restricted case. Journal of Guidance, Control, and Dynamics, 39(1):173–176.

[Blatman and Sudret, 2010] Blatman, G. and Sudret, B. (2010). An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis. Probabilistic Engineering Mechanics, 25(2):183–197.

[Bower et al., 2009] Bower, G. C., Flanzer, T. C., and Kroo, I. M. (2009). Formation geometries and route optimization for commercial formation flight. In Proceedings of the 27th AIAA Applied Aerodynamics Conference, San Antonio, TX, USA.

[Buhmann, 2003] Buhmann, M. D. (2003). Radial Basis Functions: Theory and Implementations. Cambridge University Press.

[Camilleri, 2018] Camilleri, M. A. (2018). Aircraft operating costs and profitability. In Travel Marketing, Tourism Economics and the Airline Product, chapter 12, pages 191–204. Springer.

[Caprace et al., 2019] Caprace, D.-G., Winckelmans, G., Chatelain, P., and Eldredge, J. D. (2019). Wake vortex detection and tracking for aircraft formation flight. In Proceedings of the AIAA Aviation 2019 Forum, Dallas, TX, USA.

[Cavalier et al., 1990] Cavalier, T. M., Pardalos, P. M., and Soyster, A. L. (1990). Modeling and integer programming techniques applied to propositional calculus. Computers and Operations Research, 17(6):561–570.

[Cook et al., 2015] Cook, A., Blom, H. A., Lillo, F., Mantegna, R. N., Micciche, S., Rivas, D., Vázquez, R., and Zanin, M. (2015). Applying complexity science to air traffic management. Journal of Air Transport Management, 42:149–158.

[Dee et al., 2011] Dee, D. P., Uppala, S. M., et al. (2011). The ERA-Interim reanalysis: Configuration and performance of the data assimilation system. Quarterly Journal of the Royal Meteorological Society, 137(656):553–597.

[Durango et al., 2016] Durango, G. J., Lawson, C., and Shahneh, A. Z. (2016). Formation flight investigation for highly efficient future civil transport aircraft. The Aeronautical Journal, 120(1229):1081–1100.

[Eurocontrol, 2013] Eurocontrol (2013). User manual for the Base of Aircraft Data (BADA). Technical Report 130416, Eurocontrol.

[Fahroo and Ross, 2000] Fahroo, F. and Ross, I. M. (2000). A spectral patching method for direct trajectory optimization. Journal of the Astronautical Sciences, 48(2):269–286.

[Fahroo and Ross, 2008] Fahroo, F. and Ross, I. M. (2008). Advances in pseudospectral methods for optimal control. In Proceedings of the AIAA Guidance, Navigation and Control Conference and Exhibit, Honolulu, HI, USA.

[Falck et al., 2021] Falck, R., Gray, J. S., Ponnapalli, K., and Wright, T. (2021). Dymos: A Python package for optimal control of multidisciplinary systems. Journal of Open Source Software, 6(59):2809.

[Flanzer et al., 2014] Flanzer, T. C., Bieniawski, S. R., and Blake, W. B. (2014). Operational analysis for the formation flight for aerodynamic benefit program. In Proceedings of the 52nd Aerospace Sciences Meeting, National Harbor, MD, USA.

[Flanzer et al., 2020] Flanzer, T. C., Bieniawski, S. R., and Brown, J. A. (2020). Advances in cooperative trajectories for commercial applications. In Proceedings of the AIAA SciTech 2020 Forum, Orlando, FL, USA.

[Garg, 2011] Garg, D. (2011). Advances in Global Pseudospectral Methods for Optimal Control. PhD thesis, University of Florida, Gainesville, FL, USA.

[Garg et al., 2011] Garg, D., Patterson, M. A., Francolin, C., Darby, C. L., Huntington, G. T., Hager, W. W., and Rao, A. V. (2011). Direct trajectory optimization and costate estimation of finite-horizon and infinite-horizon optimal control problems using a Radau pseudospectral method. Computational Optimization and Applications, 49(2):335–358.

[Garg et al., 2009] Garg, D., Patterson, M. A., Hager, W., Rao, A. V., Benson, D. A., and Huntington, G. T. (2009). An overview of three pseudospectral methods for the numerical solution of optimal control problems. Advances in the Astronautical Sciences, 135:1–17.

[Haalas et al., 2014] Haalas, D. J., Bieniawski, S. R., Whitehead, B., Flanzer, T., and Blake, W. B. (2014). Formation flight for aerodynamic benefit simulation development and validation. In Proceedings of the 52nd Aerospace Sciences Meeting, National Harbor, MD, USA.

[Hallock and Holzapfel, 2018] Hallock, J. N. and Holzäpfel, F. (2018). A review of recent wake vortex research for increasing airport capacity. Progress in Aerospace Sciences, 98:27–36.

[Hart et al., 2012] Hart, W. E., Laird, C. D., Watson, J.-P., Woodruff, D. L., Hackebeil, G. A., Nicholson, B. L., and Siirola, J. D. (2012). Pyomo-Optimization Modeling in Python. Springer.

[Hartjes et al., 2018] Hartjes, S., van Hellenberg Hubar, M. E., and Visser, H. G. (2018). Multiple-phase trajectory optimization for formation flight in civil aviation. In Dolega, B., Glebocki, R., Kordos, D., and Zugaj, M., editors, Advances in Aerospace Guidance, Navigation and Control, pages 389–405. Springer.

[Hartjes et al., 2019] Hartjes, S., Visser, H. G., and van Hellenberg Hubar, M. E. (2019). Trajectory optimization of extended formation flights for commercial aviation. Aerospace, 6(9):100.

[Héder, 2017] Héder, M. (2017). From NASA to EU: the evolution of the TRL scale in Public Sector Innovation. The Innovation Journal, 22(2):1–23.

[Hu et al., 2005] Hu, J., Prandini, M., and Sastry, S. (2005). Aircraft conflict prediction in the presence of a spatially correlated wind field. IEEE Transactions on Intelligent Transportation Systems, 6(3):326–340.

[Hull, 2007] Hull, D. G. (2007). Fundamentals of Airplane Flight Mechanics. Springer.

[Huntington, 2007] Huntington, G. T. (2007). Advancement and analysis of a Gauss pseudospectral transcription for optimal control problems. PhD thesis, Massachusetts Institute of Technology, Cambridge, MA, USA.

[Huntington and Rao, 2008] Huntington, G. T. and Rao, A. V. (2008). Comparison of global and local collocation methods for optimal control. Journal of Guidance, Control, and Dynamics, 31(2):432–436.

[IATA, 2013] IATA (2013). Technology roadmap, 4th Edition. Technical report, International Air Transport Association.

[IATA, 2016] IATA (2016). Press Release No. 59. IATA forecasts passenger demand to double over 20 years. https://www.iata.org/en/pressroom/pr/2016-10-18-02/. Accessed: 2021-12-19.

[ICAO, 1944] ICAO (1944). Convention on international civil aviation. Annex 2. Technical report, International Civil Aviation Organization.

[ICAO, 2016] ICAO (2016). Environmental report 2016. Aviation and climate change. Technical report, International Civil Aviation Organization.

[ICAO, 2017] ICAO (2017). Airline operating costs and productivity. Technical report, International Civil Aviation Organization.

[IPCC, 2018] IPCC (2018). Global warming of 1.5C. Technical report, The Intergovernmental Panel on Climate Change.

[Jeung et al., 2010] Jeung, H., Yiu, M. L., Zhou, X., Jensen, C. S., and Shen, H. T. (2010). Discovery of convoys in trajectory databases. Proceedings of the VLDB Endowment, 1(1):1068–1080.

[Kent and Richards, 2014] Kent, T. and Richards, A. (2014). Accounting for the effect of ground delay on commercial formation flight. In Proceedings of the 2014 UKACC International Conference on Control, Loughborough, United Kingdom.

[Kent and Richards, 2015] Kent, T. E. and Richards, A. G. (2015). Analytic approach to optimal routing for commercial formation flight. Journal of Guidance, Control, and Dynamics, 38(10):1872–1884.

[Kent and Richards, 2021] Kent, T. E. and Richards, A. G. (2021). Potential of formation flight for commercial aviation: Three case studies. Journal of Aircraft, 58(2):1–14.

[Kim and Braatz, 2012] Kim, K. K. K. and Braatz, R. D. (2012). Probabilistic analysis and control of uncertain dynamic systems: Generalized polynomial chaos expansion approaches. In Proceedings of the 2012 American Control Conference, Montreal, QC, Canada.

[Klavins, 2004] Klavins, E. (2004). Communication complexity of multi-robot systems. In Boissonnat, J., Burdick, J., Goldberg, K., and Hutchinson, S., editors, Algorithmic Foundations of Robotics V, pages 275–291. Springer.

[Kless et al., 2013] Kless, J. E., Aftosmis, M. J., Ning, S. A., and Nemec, M. (2013). Inviscid analysis of extended-formation flight. AIAA Journal, 51(7):1703–1715.

[Klower et al., 2021] Klower, M., Allen, M., Lee, D., Proud, S., Gallagher, L., and Skowron, A. (2021). Quantifying aviation’s contribution to global warming. Environmental Research Letters, 16(10):104027.

[Lai et al., 2022] Lai, Y. Y., Christley, E., Kulanovic, A., Teng, C.-C., Bjorklund, A., Nordensvard, J., Karakaya, E., and Urban, F. (2022). Analysing the opportunities and challenges for mitigating the climate impact of aviation: A narrative review. Renewable and Sustainable Energy Reviews, 156:111972.

[Lee et al., 2009] Lee, D. S., Fahey, D. W., Forster, P. M., Newton, P. J., Wit, R. C., Lim, L. L., Owen, B., and Sausen, R. (2009). Aviation and global climate change in the 21st century. Atmospheric Environment, 43(22–23):3520–3537.

[Li et al., 2014] Li, X., Nair, P. B., Zhang, Z., Gao, L., and Gao, C. (2014). Aircraft robust trajectory optimization using nonintrusive polynomial chaos. Journal of Aircraft, 51(5):1592–1603.

[Lissaman and Shollenberger, 1970] Lissaman, P. B. S. and Shollenberger, C. A. (1970). Formation flight of birds. Science, 168(3934):1003–1005.

[Liu et al., 2019] Liu, C., Jiang, B., Patton, R. J., and Zhang, K. (2019). Integrated fault-tolerant control for close formation flight. IEEE Transactions on Aerospace and Electronic Systems, 56(2):839–852.

[Marks and Gollnick, 2016] Marks, T. and Gollnick, V. (2016). Influence of aircraft type and order on fuel savings gained by two aircraft formations. In Proceedings of the 30th Congress of the International Council of the Aeronautical Science, Deajeon, South Korea.

[Matsuno et al., 2015] Matsuno, Y., Tsuchiya, T., Wei, J., Hwang, I., and Matayoshi, N. (2015). Stochastic optimal control for aircraft conflict resolution under wind uncertainty. Aerospace Science and Technology, 43:77–88.

[Michel et al., 2020] Michel, D.-T., Dolfi-Bouteyre, A., Goular, D., Augere, B., Planchat, C., Fleury, D., Lombard, L., Valla, M., and Besson, C. (2020). Onboard wake vortex localization with a coherent 1.5 μm Doppler LiDAR for aircraft in formation flight configuration. Optics Express, 28(10):14374–14385.

[Murphy, 2012] Murphy, K. P. (2012). Machine learning: A probabilistic perspective. The MIT Press.

[Ning et al., 2011] Ning, A., Flanzer, T. C., and Kroo, I. M. (2011). Aerodynamic performance of extended formation flight. Journal of Aircraft, 48(3):855–865.

[Ning, 2011] Ning, S. A. (2011). Aircraft Drag Reduction Through Extended Formation Flight. PhD thesis, Stanford University, Standord, CA, USA.

[Patterson and Rao, 2014] Patterson, M. A. and Rao, A. V. (2014). GPOPS-II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming. ACM Transactions on Mathematical Software, 41(1):1–37.

[Prandini and Hu, 2009] Prandini, M. and Hu, J. (2009). Application of reachability analysis for stochastic hybrid systems to aircraft conflict prediction. IEEE Transactions on Automatic Control, 54(4):913–917.

[Rao et al., 2009] Rao, A. V., Benson, D., Huntington, G. T., Francolin, C., Darby, C. L., and Patterson, M. (2009). User’s Manual for GPOPS Version 2.1: A MATLAB Package for Dynamic Optimization Using the Gauss Pseudospectral Method.

[Rivas and Vazquez, 2016] Rivas, D. and Vazquez, R. (2016). Uncertainty. In Cook, A. and Rivas, D., editors, Complexity Science in Air Traffic Management. Routledge.

[Ross, 2020] Ross, I. (2020). Enhancements to the DIDO optimal control toolbox. arXiv preprint arXiv:2004.13112.

[Ross and Fahroo, 2004] Ross, I. M. and Fahroo, F. (2004). Pseudospectral knotting methods for solving optimal control problems. Journal of Guidance, Control, and Dynamics, 27(3):397–405.

[Saltelli et al., 2008] Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., and Tarantola, S. (2008). Global Sensitivity Analysis. The Primer. John Wiley & Sons.

[Seo et al., 2017] Seo, J., Kim, Y., Kim, S., and Tsourdos, A. (2017). Collision avoidance strategies for unmanned aerial vehicles in formation flight. IEEE Transactions on Aerospace and Electronic Systems, 53(6):2718–2734.

[Shapiro, 2003] Shapiro, A. (2003). Monte Carlo Sampling Methods. In Ruszczynski, A. and Shapiro, A., editors, Stochastic Programming, volume 10 of Handbooks in Operations Research and Management Science, pages 353–425. Elsevier.

[Shone et al., 2021] Shone, R., Glazebrook, K., and Zografos, K. G. (2021). Applications of stochastic modeling in air traffic management: Methods, challenges and opportunities for solving air traffic problems under uncertainty. European Journal of Operational Research, 292:1–26.

[Slotnick, 2014] Slotnick, J. P. (2014). Computational aerodynamic analysis for the formation flight for aerodynamic benefit program. In Proceedings of the 52nd Aerospace Sciences Meeting, National Harbor, MD, USA.

[Sudret, 2008] Sudret, B. (2008). Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering and System Safety, 93:964–979.

[Tierno et al., 2012] Tierno, M. Á. G., Cortés, M. P., and Márquez, C. P. (2012). Mecánica del vuelo. Ibergaceta.

[Torre et al., 2019] Torre, E., Marelli, S., Embrechts, P., and Sudret, B. (2019). A general framework for data-driven uncertainty quantification under complex input dependencies using vine copulas. Probabilistic Engineering Mechanics, 55:1–16.

[Trucchia et al., 2019] Trucchia, A., Egorova, V., Pagnini, G., and Rochoux, M. (2019). On the merits of sparse surrogates for global sensitivity analysis of multi-scale nonlinear problems: Application to turbulence and fire-spotting model in wildland fire simulators. Communications in Nonlinear Science and Numerical Simulation, 73:120–145.

[Tu et al., 2008] Tu, Y., Ball, M. O., and Jank, W. S. (2008). Estimating flight departure delay distributions—a statistical approach with long-term trend and short-term pattern. Journal of the American Statistical Association, 103(481):112–125.

[Vachon et al., 2002] Vachon, M. J., Ray, R., Walsh, K., and Ennix, K. (2002). F/A-18 aircraft performance benefits measured during the autonomous formation flight project. In Proceedings of the AIAA Atmospheric Flight Mechanics Conference and Exhibit, Monterey, CA, USA.

[Vechtel et al., 2018] Vechtel, D., Fischenberg, D., and Schwithal, J. (2018). Flight dynamics simulation of formation flight for energy saving using LES-generated wake flow fields. CEAS Aeronautical Journal, 9(4):735–746.

[Vieira et al., 2009] Vieira, M. R., Bakalov, P., and Tsotras, V. J. (2009). On-line discovery of flock patterns in spatio-temporal data. In Proceedings of the 17th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems, Seattle, WA, USA.

[Voskuijl, 2017] Voskuijl, M. (2017). Cruise range in formation flight. Journal of Aircraft, 54(6):2184–2191.

[Wachter and Biegler, 2006] Wachter, A. and Biegler, L. T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical programming, 106(1):25–57.

[Wei et al., 2008] Wei, S., Zěfran, M., and DeCarlo, R. A. (2008). 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.

[Wieselsberger, 1914] Wieselsberger, C. (1914). Beitrag zur Erklärung des Winkelfluges einiger Zugvögel. Zeitschrift für Flugtechnik und Motorluftschiffahrt., 5:225–229.

[Xiu, 2010] Xiu, D. (2010). Numerical Methods for Stochastic Computations. A Spectral Method Approach. Princeton University Press.

[Xu et al., 2014] Xu, J., Ning, S. A., Bower, G., and Kroo, I. (2014). Aircraft route optimization for formation flight. Journal of Aircraft, 51(2):490–501.

[Zhu, 2019] Zhu, Y. (2019). Uncertain Optimal Control. Springer.