Uncertainty Propagation

Introduction

The Air Traffic Management (ATM) system is a complex system composed of a large number of heterogeneous components, such as airports, aircraft, navigation systems, flight management systems (FMS), traffic controllers, and weather (see Kim et al. (2009)). Therefore, its performance is affected by numerous factors. Within the trajectory-based-operations concept of SESAR and NextGen, aircraft trajectories are key to study ATM operations, which are subject to many uncertainties. Sources of uncertainty for aircraft trajectories include wind and severe weather, navigational errors, aircraft performance inaccuracies, or errors in the FMS, among others. The analysis of the impact of uncertainties in aircraft trajectories and its propagation through the flight segments is of great interest, since it might help to understand how sensitive the system is to the lack of precise data and measurement errors, and, therefore, aid in the design of a more robust ATM system, with improved safety levels.

Methods to study uncertainty propagation: classification

The methods used to study trajectory uncertainty and uncertainty propagation can be classified into two main groups:

• Monte Carlo methods (Hastings (1970)): These are computational methods that rely on repeated random sampling to compute their results. That is, one randomly selects a value for all the uncertain elements that affect the trajectory, and then computes it in a deterministic way. This is repeated as many times as necessary until the resulting trajectories are a representative sample of the uncertain trajectory itself. The main advantage of these types of methods is that they can be used with all types of uncertainty and do not require any complicated computation beyond solving the equations of flight mechanics many times. However it is a very expensive method in computational terms (suffering from “the curse of dimensionality”; if there are many sources of uncertainty one needs many sample points and therefore it is not implementable in practice) and requires random sampling of the sources of uncertainty, which is not always easy. More advanced Monte Carlo variants, such as the sequential Monte Carlo (Liu and Chen (1998)), use more advanced ideas (carefully selecting the samples of the sources) to improve the convergence and reduce the required number of sample points.
• Non-Monte Carlo methods. To avoid the computational problem of Monte Carlo methods, other techniques have been proposed to study uncertainty propagation in dynamical systems. Halder and Bhattacharya (2011) classify those methods into two categories:

• Parametric: in which one evolves in time the statistical moments, such as mean and covariance, but the probability distribution remains unknown). An example is the Generalized Polynomial Chaos(GPC) method, which is described in more depth next.
• Non-parametric: in which the probability density function itself is evolved. An example is the method presented in Halder and Bhattacharya (2011) which finds the probability density function by solving a stochastic Liouville equation. Another example is the Probability density function Transformation Method (PTM), which is described in this document.
• Another family of methods which do not fall under any of the previous two categories uses Differential Algebra (Berz et al., 1999) to study the evolution in time of a neighbouring set of initial conditions. They can be used to efficiently compute the statistics of the propagated state (Valli et al., 2013). Differential algebra has been mostly used in orbital mechanics problem with uncertain initial conditions, however they are suitable to be applied in the setting of uncertain aircraft trajectories.

Generalized Polynomial Chaos

The GPC representation was introduced by Wiener (1938) and it is based on the fact that any second-order process (i.e., a process with finite second-order moments) can be represented as a Fourier-type series, with time-dependent coefficients, and using orthogonal polynomials as GPC basis functions in terms of random variables. A general introduction to GPC can be found in Xiu and Karniadakis (2002) and in Schoutens (2000), whereas details in numerical computations are studied in Debusschere et al. (2004).

Applications of GPC to aerospace problems

The method of polynomial chaos is used in the works of Prabhakar et al. (2010) and Dutta and Bhattacharya (2010) to study, respectively, uncertainty propagation and trajectory estimation, for hypersonic flight dynamics with uncertain initial data. Vazquez and Rivas (2013) develop methods to study the impact of initial mass uncertainty on fuel consumption during the cruise flight phase. Other examples is Okamoto and Tsuchiya (2015) which use it to solve the problem of optimal trajectory generation in the context of stochastic optimal control.

Traditionally, GPC requires formulating a set of nonlinear ordinary differential equations (based on the model of the system under study) for the time-dependent coefficients. This is known as intrusive polynomial chaos, and is the methodology applied in this work. However, it is worth mentioning that recently, a different approach (known as non-intrusive polynomial chaos) has emerged. This new approach does not directly solve these equations but instead computes the GPC coefficients from a number of simulations of the system under study, which can be treated as a black box. A comparison between the intrusive and non-intrusive approaches can be found in Xiong et al. (2014) Recent applications include Jones et al. (2013) which employ nonintrusive GPC for nonlinear propagation of orbits, and Li et al. (2014) which use it for robust aircraft trajectory optimization.

Essentials of the GPC method (intrusive version)

Start from a differential equation

[math]\dot x(t)=f(x(t),\mu),\quad x(0)=x_0[/math]

where either the initial condition [math]x_0[/math] or the parameter [math]\mu[/math] are random numbers, coming from a given probability distribution. Then [math]x(t)[/math] itself becomes random.

To compute the mean and variance of [math]x(t)[/math], the Generalized Polynomial Chaos (GPC) method can be used. The basic idea is to represent [math]x(t)[/math] as a Fourier-type series, with time-dependent coefficients, and orthogonal polynomials in terms of random variables used as basis functions. The orthogonal polynomials used in GPC are chosen from the Askey scheme (a way of organizing certain orthogonal polynomials into a hierarchy, see Askey and Wilson (1985)), depending on the probability distribution of the parameter or initial condition.

The main results is that if one chooses a family of polynomials which are orthogonal, the convergence of the series is exponential. The orthogonality property implies that, when taking expectation with respect to the random variable for two polynomials of the family [math]\phi_i[/math] and [math]\phi_j[/math], then [math]\mathrm{E}[\phi_i \phi_j]=\delta_{ij} \mathrm{E}[\phi_i^2][/math], where [math]E[\cdot] [/math] is the mathematical expectation and [math]\delta_{ij}[/math] is the Kronecker's delta function. For instance, for the uniform distribution (typically denoted as [math]\Delta[/math]), the adequate orthogonal polynomials are the Legendre polynomials [math]L_n(\Delta)[/math].

To apply the GPC method, one first write [math]x(t)[/math] in terms of the orthogonal polynomials.

[math] \begin{equation} x(t) = \sum_{i=0}^P h_i(t) L_i(\Delta) \end{equation} [/math]

where [math]P[/math] is the order of the approximation, which is to be taken sufficiently large. The coefficients [math]h_i(t)[/math] are to be found using the differential equation of [math]x(t)[/math]. The advantage of the GPC method is that a small or moderate value of [math]P[/math] is enough to get good results, thus resulting in a method that is not very intensive computationally.

Once the coefficients [math]h_i(t)[/math] are found, one can compute approximate values of the fuel mass mean and variance at all times, as follows. For the mean one has

[math] \begin{eqnarray} \mathrm{E}[x(t)] &=& \sum_{i=0}^P h_i(t) \mathrm{E}[ L_i(\Delta)] \nonumber \\ &=&\sum_{i=0}^P h_i(t) \mathrm{E}[ L_i(\Delta)L_0(\Delta)] \nonumber \\ &=&h_0(t) \mathrm{E}[ L_0^2(\Delta)] \nonumber \\ &=&h_0(t) \end{eqnarray} [/math]

and for the variance

[math] \begin{eqnarray} \mathrm{Var}[x(t)] &=& \mathrm{E}[x(t)]-\mathrm{E}[x(t)]^2 \nonumber \\ &=& \sum_{i=0}^P \sum_{j=0}^P h_i(t) h_j(t) \mathrm{E}[ L_i(\Delta)L_j(\Delta)]-h_0(t)^2 \nonumber \\ &=& \sum_{i=1}^P h_i^2(t) \mathrm{E}[ L_i^2(\Delta)] \end{eqnarray} [/math]

Probability density function Transformation Method (PTM)

This method was developed in Vazquez and Rivas (2013) to study the evolution of the aircraft mass probability density function in cruise flight with uncertain initial mass.

Start, as in GPC, from a differential equation

[math]\dot x(t)=f(x(t),\mu),\quad x(0)=x_0[/math]

where either the initial condition [math]x_0[/math] or the parameter [math]\mu[/math] are random numbers, coming from a known probability distribution. Then [math]x(t)[/math] itself becomes random.

Recall that, given a random variable [math]z[/math] with probability density function [math]f_z(z)[/math], if one defines another random variable[math]y[/math] using a transformation [math]g[/math] such that [math]y=g(z)[/math], then it is known that the probability density function [math]f_y(y)[/math] of [math]y[/math] is given by (see for instance Canavos (1984) ).

[math] \begin{equation} f_y(y)=\frac{f_z(g^{-1}(y))}{\vert g'(g^{-1}(y))\vert} \end{equation} [/math]

an expression which is valid only if the function [math]g(z)[/math] is invertible on the domain of [math]z[/math].

To analyze the probability density function of [math]x(t)[/math] we consider it as a continuous (in time) transformation of the parameter [math]\mu[/math] or the initial condition [math]x_0[/math]. For instance if the initial condition is random, then we consider [math]x(t)=g_t(x_0)[/math]. Then, applying the transformation formula

[math] \begin{equation} f_{x(t)}(x(t))=\frac{f_{x_0}(g_t^{-1}(x(t)))}{\vert g_t'(g_t^{-1}(x(t)))\vert} \end{equation} [/math]

There is no problem of inverting [math]g_t[/math] since the solutions of differential equations are uniquely defined by their initial conditions (assuming the differential equation is regular enough).

The idea of the method is to numerically approximate the transformation equation. For that, take [math]n[/math] consecutive points from the domain of [math]x_0[/math], denoted as [math]x_0^i, i=1,...,n[/math]. Solving the differential equation of [math]x(t)[/math] for each [math]i[/math] with initial condition [math]x_0^i[/math], one can compute the value of [math]x^i(t)[/math].

The numerator of the transformation equation is computed for each [math]i[/math] as [math]f_{x_0}(x_0^i)[/math]. To compute the denominator of the transformation equation, the function [math]g'_t(x_0)[/math] is needed; this function is obtained in terms of

[math] \begin{equation} \phi(t;x_0)\equiv\frac{\partial x(t;x_0)}{\partial x_0} \end{equation} [/math]

which is the sensitivity function of the solution [math]x(t)[/math] with respect to the initial condition [math]x_0[/math]. The process to obtain [math]\phi(t;x_0)[/math] can be found in any ordinary differential equations textbook, see for instance Hale (1969).

Once the function [math]\phi(t;x_0)[/math] is found, one has the values of the probability density function [math] f_{x(t)}(x(t))[/math] at the [math]n[/math] points [math]x^i(t),i=1,...,n[/math], as

[math] \begin{equation} f_{x(t)}(x(t)) =\frac{f_{x_0}(x_0^i)}{\vert\phi(t;x_0^i)\vert} \end{equation} [/math]

From the distribution function, the mean and the typical deviation can be computed as

[math] \begin{eqnarray} E[x(t)]&=&\int_{0}^{\infty} x(t)f_{x(t)}(x(t)) dx(t) \\ (\sigma[x(t)])^2&=&\int_{0}^{\infty} x^2(t)f_{x(t)}(x(t)) dx(t)-\left(E[x(t)]\right)^2 \end{eqnarray} [/math]

Application: Fuel consumption uncertainty

The required fuel load of an aircraft for a given cruise range with either uncertain average wind or uncertain initial mass is studied next to exemplify the GPC and PTM methods. The fuel mass probability density function is obtained using PTM, whereas Generalized Polynomial Chaos is used to study the mean and typical deviation of the required fuel mass. The dynamics of mass evolution in cruise flight is defined by a nonlinear equation, which can be solved analytically to obtain the fuel mass; this exact solution is used to assess the accuracy of the proposed methods.

The fuel load can be found from the evolution equation of aircraft mass. In cruise flight, the equations of flight mechanics for flight in a vertical plane (constant heading) are considered, under the following hypothesis: symmetric flight, flat Earth model, constant altitude, and constant velocity. Then, the equation of mass evolution is (taking into account the equation [math]T=D[/math], [math]T[/math] being the thrust)

[math] \begin{equation} \dot m=-cD \end{equation} [/math]

where [math]m[/math] is the aircraft mass, [math]D[/math] the aerodynamic drag, and [math]c[/math] the specific fuel consumption, which is considered constant under the previously stated hypothesis. The cruise range [math]x_f[/math] is given.

The drag can be written as [math]D=\frac{1}{2} \rho V^2 S C_D[/math], where [math]\rho[/math] is the density, [math]V[/math] the velocity, [math]S[/math] the wing surface area, and the drag coefficient [math]C_D[/math] is modeled by a parabolic polar of constant coefficients [math]m[/math], where [math]C_D=C_{D_0}+k C_L^2[/math], [math]C_L[/math] is the lift coefficient given by [math]C_L=\frac{L}{\frac{1}{2} \rho V^2 S}=\frac{mg}{\frac{1}{2} \rho V^2 S}[/math], where the equation of motion [math]L=mg[/math] ( [math]g[/math] is the acceleration of gravity) has been used.

Using these definitions, an autonomous equation for mass evolution is obtained:

[math] \begin{equation} \dot m=- c\left( \frac{1}{2} \rho V^2 S C_{D_0}+ m^2\frac{k g^2}{\frac{1}{2} \rho V^2 S} \right) \end{equation} [/math]

Thus, one can write [math] \begin{equation} \dot m=-(A+Bm^2) \end{equation} [/math] where the constants [math]A[/math] and [math]B[/math] are defined as [math]A=\frac{c}{2} \rho V^2 S C_{D_0}[/math] and [math]B=\frac{2c k g^2}{ \rho V^2 S}[/math]. Note that [math]A,B\gt 0[/math]. The mass equation is a nonlinear equation describing the evolution of mass during cruise flight as a function of time.

To find the evolution of mass as a function of distance, consider the kinematic equation

[math] \begin{equation} \dot x=V + w \end{equation} [/math]

where [math]x[/math] is the horizontal distance, and [math]w[/math] is the average wind speed, considered constant. Combining the mass equation and the kinematic equation to eliminate time, one reaches

[math] \begin{equation} \frac{dm}{dx}=-\frac{A+Bm^2}{V+w} \end{equation} [/math]

The following values are used: [math]C_{d_0}=0.015[/math],[math]C_{d_2}=0.042[/math],[math]\rho=0.5 \rho_0[/math],[math]\rho_0=1.225~\mathrm{kg/m^3}[/math],[math]V=200~\mathrm{m/s}[/math],[math]c=5\cdot 10^{-5}~\mathrm{s/m}[/math],[math]S=150~\mathrm{m^2}[/math],[math]g=9.8~\mathrm{m/s^2}[/math],[math]m_F=25000~\mathrm{kg}[/math], and [math]x_f=2500~\mathrm{km}[/math].

Fuel uncertainty due to initial mass uncertainty

It is realistic to consider that the initial mass [math]m_0[/math] is not a deterministic variable which is known a priori, but rather a random variable which is not known. Then, the solution of the mass equation is still valid but in a probabilistic sense, i.e., writing the solution as [math]m(t;m_0)[/math] to emphasize its dependence on the initial mass, this solution is a random process. If the distribution of [math]m_0[/math] is known, it is possible to apply the previously reviewed methods to study the distribution of the aircraft fuel mass consumption as well as its statistical properties (mean, variance, typical deviation).

Consider in particular that [math]m_0[/math] is distributed as a uniform continuous variable

[math] \begin{equation} f_{m_0}(m_0)=\left\{ \begin{array}{lc} \frac{1}{2\delta_m} ,& m_0\in\left[\bar m_0-\delta_m, \bar m_0+\delta_m \right] \\ 0,& m_0\notin\left[\bar m_0-\delta_m, \bar m_0+\delta_m \right] \end{array} \right. \end{equation} [/math]

where [math]\bar m_0[/math] is the mean and [math]\delta_m[/math] the width of the uniform distribution.

Since the cruise range,[math]x_f[/math], is known, then the final value of the aircraft mass [math]m(x_f;m_0)[/math] is given by the mass equation particularized for [math]x=x_f[/math], and the fuel consumption during the cruise is

[math] \begin{equation} m_F(m_0)=m_0-m(x_f;m_0) \end{equation} [/math]

Assuming the values [math]\bar m_0= 81633~\mathrm{kg}[/math] and [math]\delta_m=5000~\mathrm{kg}[/math] and zero wind, the PTM method finds the following probability density function for the fuel mass consumption.

Figure 1. Distribution function of fuel mass computed with PTM.

Using GPC, the following values are found for mean and standard deviation of fuel mass consumption: [math]\mathrm{E}[m_F]=23892~\mathrm{kg}[/math], [math]\sigma[m_F]=499.96~\mathrm{kg}[/math].

Fuel uncertainty due to average wind uncertainty

In this case, [math]w[/math] is random. In particular, an uniform distribution is assumed. Hence, the probability density function is

[math] \begin{equation} f_w(w)=\left\{ \begin{array}{lc} \frac{1}{2\delta_w} ,& w\in\left[-\delta_w, \delta_w \right] \\ 0,& w\notin\left[-\delta_w, \delta_w\right] \end{array} \right. \end{equation} [/math]

To emphasize the dependence of the mass [math]m(x)[/math] on the wind, the mass is written as [math]m(x;w)[/math], even though often it is just denoted as [math]m[/math] for the sake of simplicity.

In this particular application, the final aircraft mass [math]m_f[/math] is given (fixing [math]m_f[/math], instead of [math]m_0[/math], is consistent with having a fixed landing weight). Then, the mass equation is solved backwards with the boundary condition

[math] \begin{equation} m(x_f)=m_f \end{equation} [/math]

and finally the resulting fuel mass consumed is found from

[math] \begin{equation} m_F(w)=m(0;w)-m_f \end{equation} [/math]

Assuming the values [math]\bar w= 0~\mathrm{m/s}[/math] and [math]\delta_w=50~\mathrm{m/s}[/math] and zero wind, the PTM method finds the following probability density function for the fuel mass consumption.

Using GPC, the following values are found for mean and standard deviation of fuel mass consumption: [math]\mathrm{E}[m_F]=23941.7~\mathrm{kg}[/math], [math]\sigma[m_F]=3924.9~\mathrm{kg}[/math].