Propagation of Initial Mass Uncertainty in Aircraft Cruise Flight

Propagation of Initial Mass Uncertainty in Aircraft Cruise Flight Rafael Vazquez∗and Dami´ an Rivas† Escuela Superior de Ingenieros, Universidad de Sevilla, Seville, 41092, Spain The propagation of initial mass uncertainty in cruise flight is studied. Two cruise conditions are analyzed: one with given cruise fuel load, and the other with given cruise range. Two different distributions of initial mass are considered: uniform and of gamma type. The Generalized Polynomial Chaos method is used to study the evolution of mean and variance of the aircraft mass. To compute the mass distribution function as a function of time, two approximate methods are developed. These methods are also applied to study the distribution functions of the flight time (in the case of given fuel load), and of the fuel consumption (in the case of given range). The dynamics of mass evolution in cruise flight is defined by a nonlinear equation, which can be solved analytically; this exact solution is used to assess the accuracy of the proposed methods. Comparison of the numerical results with the exact analytical solutions shows an excellent agreement in all cases, hence verifying the methods developed in this work. Nomenclature A, B constants of the problem CD, CLdrag and lift coefficients CD0, CD2coefficients of the drag polar cspecific fuel consumption Daerodynamic drag E[.]expectation fxprobability density function of random variable x ∗Associate Professor, Department of Aerospace Engineering; rv[email protected]. †Professor, Department of Aerospace Engineering; dri[email protected]. 1 G(k; 1) gamma distribution with scale parameter equal to one ggravity acceleration hicoefficients of the GPC expansion kshape parameter of the gamma distribution Llift LnLegendre polynomials M0minimum value of m0with nonzero probability for the gamma distribution maircraft mass m0initial aircraft mass ¯m0mean of the initial mass distribution mFfuel load Swing surface area Tthrust ttime tfflight time Vaircraft speed V ar[.]variance xhorizontal distance xfrange ∆standard uniform distribution δmwidth of the uniform distribution Γ(a)Euler gamma function Γ(a, b)incomplete Euler gamma function φk−1 ngeneralized Laguerre polynomials σ[.]typical deviation I. Introduction The Air Traffic Management (ATM) system is a very complex system which contains a large number of heterogeneous components, such as airports, aircraft, navigation systems, flight management systems (FMS), traffic controllers, and weather (see Kim et al. [1]). Correspondingly, its performance is affected by numerous factors. Within the trajectory-based-operations concept of SESAR and NextGen, aircraft trajectories are key to study ATM opera- 2 tions, 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. Among those sources, weather uncertainty has perhaps the greatest impact on ATM operations, being responsible for much of the delays. Its analysis has been addressed by many authors, using different methods, for example the following. Nilim et al. [2] consider a trajectory-based air traffic management scenario to minimize delays under weather uncertainty, where the weather processes are modeled as stationary Markov chains. Pepper et al. [3] present a method, based on Bayesian decision networks, of accounting for uncertain weather information in air traffic flow management. Clarke et al. [4] develop a methodology to study airspace capacity in the presence of weather uncertainty and formulates a stochastic dynamic programming algorithm for traffic flow management. Zheng and Zhao [5] develop a statistical model of wind uncertainties and apply it to stochastic trajectory prediction in the case of straight, level aircraft flight trajectories. The framework for this work is the analysis of uncertainty propagation in aircraft trajectories, and, eventually, its effect on the ATM system. In this paper several tools are presented to analyze uncertainty propagation in a nonlinear problem, and they are applied to study the effect of initial aircraft mass uncertainty and its propagation through the cruise flight phase. The relevance of this problem resides in two facts: first, the initial mass is an important source of uncertainty in trajectory prediction, which determines mass evolution and, therefore, fuel consumption and flight cost, and, second, cruise uncertainties have a large impact on the overall flight since the cruise phase is the largest portion of the flight (at least for long-haul routes). In the applications, two cruise conditions are studied: one with given cruise fuel load, and the other with given cruise range. Several methods have been proposed to study uncertainty propagation in dynamical systems, beyond the classical Monte-Carlo methods (which can be very expensive computationally). Halder and Bhattacharya [6] classify those methods it two categories: parametric (in which one evolves the statistical moments) and non-parametric (in which the probability density function is evolved). They address the problem of uncertainty propagation in planetary entry, descent, and landing, using a non-parametric method that reduces to solving the stochastic Liouville equation. In this paper, the evolution in time of the mean and the variance of the aircraft mass is studied using the Generalized Polynomial Chaos (GPC) method (a parametric method according to Ref. [6]). The GPC representation was introduced by Wiener [7] and it is based on the fact that any second-order process (i.e., a process with finite secondorder 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 [8] and in Schoutens [9], whereas details in numerical computations are studied in Debusschere et 3 al. [10]. The method of polynomial chaos is used in the works of Prabhakar et al. [11] and Dutta and Bhattacharya [12] to study, respectively, uncertainty propagation and trajectory estimation, for hypersonic flight dynamics with uncertain initial data, and by Fisher and Bhattacharya [13] in the problem of optimal trajectory generation in the context of stochastic optimal control. Also, the distribution function of the aircraft mass is analyzed using two approximate methods developed in this paper (non-parametric methods according to Ref. [6]). One method is based on the resolution of the variational equation for the sensitivity function with respect to the initial condition, and the other is based on the computation of the probability measure of the random variable as a function of time. These two methods are also applied, first, to the analysis of the distribution function of the flight time, in the case of given fuel load, and, second, to the analysis of the distribution function of the fuel consumption, in the case of given range. In this way, the effect of the initial mass uncertainty in flight properties other than mass is studied as well. In this paper, the case of cruise at constant altitude and constant speed is considered (cruise segments defined by these two flight constraints are commonly flown by commercial aircraft, according to air traffic control procedures). In this case, the evolution of aircraft mass is defined by a nonlinear equation which can be solved analytically. Results are presented for two different distributions of initial mass (uniform and of gamma type). The analytical solutions represent benchmark solutions which are used to assess the accuracy of the proposed methods. Comparison with the exact analytical results is made, showing an excellent agreement in all cases. This paper is organized as follows. First the problem of mass evolution in cruise flight is solved. Then, in Section III, the two initial mass distributions considered are described. In Section IV mean and variance of the mass distribution are analyzed using the GPC method. In Section V the two non-parametic methods developed to study the evolution of distribution functions are presented, and applied to the mass distribution function. These two methods are used, in Section VI, to study the distribution functions of flight time and fuel consumption. Some numerical results are presented in Section VII, and some conclusions are drawn in Section VIII. Finally, the exact analytical solutions are presented in the Appendix. II. Mass evolution in cruise flight The equations of motion for symmetric flight in a vertical plane (constant heading), using a flat Earth model, for constant altitude and constant speed are (see Ref. [14]) dx dt=V, dm dt=−cT T=D, L =mg (1) 4 where xis the horizontal distance, tthe time, Vthe speed, Tthe thrust, Dthe aerodynamic drag, Lthe lift, mthe aircraft mass, gthe acceleration of gravity, and cthe specific fuel consumption, which can be taken as a function of altitude and speed, and it is therefore constant under the given cruise condition. The drag can be written as D=1 2ρV 2SCD, where ρis the density, Sthe wing surface area, and the drag coefficient CDis modeled by a parabolic polar CD=CD0+CD2C2 L, where CLis the lift coefficient given by CL=2L ρV 2S, and the coefficients CD0and CD2are constant under the given cruise condition. Using these definitions and Eqs. (1), an autonomous equation for the mass evolution is obtained: dm dt=−c1 2ρV 2SCD0+m22CD2g2 ρV 2S(2) Thus, one can write dm dt=−(A+Bm2)(3) where the constants Aand Bare defined as A=c1 2ρV 2SCD0and B=c2CD2g2 ρV 2S. Note that A, B > 0. Equation (3) is a nonlinear equation describing the evolution of mass during cruise flight, to be solved with the initial condition m(0) = m0(4) To emphasize the dependence of the mass m(t)on the initial condition, the mass is written as m(t;m0), even though often it is just denoted as mfor the sake of simplicity. The explicit solution of Eqs. (3) and (4) is m(t;m0) = rA B m0−rA Btan √ABt rA B+m0tan √ABt (5) A. Cruise with given fuel load For the case in which the cruise fuel load is given, denoting the given mass of fuel as mF< m0, the solution obtained given by Eq. (5) is valid in the time interval t∈[0, tf(m0)], where tf(m0)(the flight time) is obtained from m(tf(m0); m0) = m0−mF. From Eq. (5) one can directly compute this time as tf(m0) = 1 √AB arctan √ABmF A+Bm0(m0−mF)!(6) Note that tfis a monotonically decreasing function of m0. Thus, for a given amount of fuel, the larger m0, the smaller tf, and, as a consequence, the smaller the distance traveled by the aircraft. The initial mass m0is unbounded, and has a lower limit equal to mF(although these limits are not physically meaningful). Thus for m0∈(mF,∞)one 5 obtains, from Eq. (6), tf∈ 0,1 √AB arctan rB AmF!!. Also, since m(tf;m0) = m0−mFthe final value of the aircraft mass satisfies m(tf;m0)∈(0,∞). In the next sections, the evolution of mass and the behavior of the flight time are studied for an uncertain value of the initial mass, while the rest of the parameters (some of them embedded in the constants Aand B) have a fixed value. B. Cruise with given range For the case in which the cruise range is given, taking xas the independent variable, one has dm dx =−1 V(A+Bm2)(7) and the same initial condition Eq. (4). The explicit solution of Eqs. (7) and (4) is m(x;m0) = rA B m0−rA Btan 1 V√ABx rA B+m0tan 1 V√ABx(8) If the given cruise range is xf, then the final value of the aircraft mass m(xf;m0)is given by Eq. (8) particularized for x=xf, and the fuel consumption during the cruise is mF(m0) = m0−m(xf;m0) = m2 0+A Btan 1 V√ABxf rA B+m0tan 1 V√ABxf(9) Note that mFis a monotonically increasing function of m0: the larger m0, the larger the fuel consumption. As before, m0is unbounded, and, in order to have m(xf;m0)>0, it has a lower limit equal to rA Btan 1 V√ABxf. Thus for m0∈ rA Btan 1 V√ABxf,∞!one obtains from Eq. (9) that mF∈ rA Btan 1 V√ABxf,∞!. Also, from Eq. (8) the final value of the aircraft mass satisfies m(xf;m0)∈ 0,rA Btan 1 V√ABxf−1!. In the next sections, the behavior of the fuel consumption is studied for an uncertain value of the initial mass, while, as before, the rest of the parameters have a fixed value. In this case the flight time is known, trivially given by tf=xf V. 6 III. Initial mass distribution It is realistic to consider that the initial mass m0is not a deterministic variable which is known a priori, but rather a random variable which is not known. Then, the solution given by Eq. (5) is still valid but in a probabilistic sense, i.e., m(t;m0)is a random process. If the distribution of m0is known, it is possible to study the time evolution of the distribution of the aircraft mass m(t;m0), as well as its statistical properties (mean, variance, typical deviation). In this work, to analyze mass evolution, two probabilistic models for m0are considered: uniform and gamma distributions, which are described next. Note that a Gaussian distribution representing the initial mass uncertainty would be non-physical, since it would allow (with small but nonzero probability) negative initial mass, and, therefore, it is not considered in this paper. A. Uniform distribution First it is considered that m0is distributed as a uniform continuous variable whose probability density function is fm0(m0) = 1 2δm in the interval [ ¯m0−δm,¯m0+δm], and zero otherwise, where ¯m0is the mean and δmthe width of the uniform distribution, as shown in Fig. 1. Denoting by ∆the standardized uniform distribution taking values in the interval [−1,1], one has that m0= ¯m0+δm∆. The mean of m0is E[m0] = R∞ 0m0fm0(m0)dm0= ¯m0, where E[·]is the mathematical expectation, and the variance of m0is Var[m0] = E[m2 0]−(E[m0])2=δ2 m 3. m0 fm0 k Gamma distribution Uniform distribution ¯m0+δm ¯m0−δm¯m0 Fig. 1 Shape of the probability density functions of the initial mass. B. Gamma distribution The gamma distribution (see Ref. [15]) represents a continuous nonnegative random variable, and is denoted by G(k, θ), where k > 0is the shape parameter and θ > 0is the scale parameter. It is known that E[G(k, θ)] = kθ and 7 Var[G(k, θ)] = kθ2, and that the probability density function of G(k, θ)is f(x;k, θ) = xk−1e−x/θ θkΓ(k)for x≥0(and zero otherwise), where Γis the Euler gamma function. Using the property that for θ > 0one has G(k, θ) = θG(k, 1), the value θ= 1 is considered in this paper without loss of generality. To represent the initial mass distribution, let m0= ¯m0+δm √3k(G(k, 1) −k), where ¯m0and δmare the same values chosen for the uniform distribution. Hence, only the values m0≥M0have nonzero probability, where M0(the minimum possible value of mass for the given values of ¯m0and δm) is obtained making G(k, 1) = 0 and it is given by M0= ¯m0−δm √3√k(10) Thus, one has the following probability density function fm0(m0)=(m0−M0)k−1e −(m0−M0) √3k δm δm √3kk Γ(k) , m0≥M0(11) and zero otherwise. In this way, one has E[m0] = ¯m0and Var[m0] = δ2 m 3(independently of k), as for the previously chosen uniform distribution. Note that for k→ ∞, one has G(k, 1) −k √k→N(0,1) which implies m0→N¯m0,δ2 m 3, i.e., for large kthe gamma distribution resembles a Gaussian distribution. However, the maximum value of kis limited by the fact that M0should be greater than zero. Therefore, the value of kmust be chosen taking into account Eq. (10). In Fig. 1 the shape of the probability density function of m0is plotted for different values of k, and compared with the uniform distribution. IV. Analysis of mass mean and variance To compute the mean and variance of the mass, the Generalized Polynomial Chaos (GPC) method is used (see Ref. [7]), in which the process is represented as a Fourier-type series, with time-dependent coefficients, and orthogonal polynomials in terms of random variables are 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 Ref. [16]). 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 φiand φj, then E[φiφj] = δijE[φ2 i], where δij is the Kronecker delta. For the uniform distribution ∆, the adequate orthogonal polynomials are the Legendre polynomials Ln(∆), whereas for the gamma distribution G(k, 1) one must use the generalized Laguerre polynomials φk−1 n(G). 8 To apply the GPC method, one first write the initial mass distribution m0in terms of the orthogonal polynomials. For the uniform distribution, one can write m0= ¯m0L0(∆)+δmL1(∆), whereas for the gamma distribution it follows that m0= ¯m0φk−1 0(G)−δm √3kφk−1 1(G). In the following the uniform distribution case is considered (the gamma distribution is handled analogously). It is assumed that m(t;m0)can be written as m(t;m0) = P X i=0 hi(t)Li(∆) (12) where the coefficients hiare to be found using the mass equation (3), and Pis the order of the approximation, which is to be taken sufficiently large. Substituting Eq. (12) in Eq. (3), the following equation is obtained P X i=0 ˙ hi(t)Li(∆) = −A−B P X i=0 P X j=0 hi(t)hj(t)Li(∆)Lj(∆) (13) Now, multiplying Eq. (13) by Ll(∆) for l= 0, . . . , P, taking expectation with respect to ∆, and using the orthogonality property of the Llpolynomials, one obtains P+ 1 equations ˙ hl(t)E[L2 l(∆)] = −Aδ0l−B P X i=0 P X j=0 hi(t)hj(t)E[Li(∆)Lj(∆)Ll(∆)], l = 0, . . . , P (14) and calling Cijl =E[LiLjLl] E[L2 l](which is a number that can be exactly computed since the involved expectations are just integrals of polynomials) it follows that ˙ hl=−Aδ0l−B P X i=0 P X j=0 hihjCijl, l = 0, . . . , P (15) which is a system of P+ 1 nonlinear coupled ordinary differential equations. The same result is reached for the gamma distribution case, with the corresponding Cijl coefficients. The initial conditions of Eqs. (15) depend on the initial mass distribution. For the uniform distribution case they are h0(0) = ¯m0, h1(0) = δm, hl(0) = 0,for l= 2, . . . , P (16) whereas for the gamma distribution they are given by h0(0) = ¯m0, h1(0) = −δm √3k, hl(0) = 0,for l= 2, . . . , P (17) The advantage of the GPC method is that a small or moderate value of Pis enough to get good results, thus resulting 9 intervals of mF. Thus, it is deduced that Pr(mi F≤mF≤mi+1 F) = Pr(mi 0≤m0≤mi+1 0), hence fmF(ξi) = Pr(mi 0≤m0≤mi+1 0) mi+1 F−mi F , i = 1, . . . , n −1(47) Taking fmF(m1 F) = fmF(ξ1) fmF(mi F) = fmF(ξi−1) + fmF(ξi) 2, i = 2, . . . , n −1 fmF(mn F) = fmF(ξn−1) (48) an approximation of fmFis obtained at npoints. VII. Results Now, the methods presented in previous sections are applied to the two initial mass distributions defined in Section III. The numerical resolution of the different problems is performed using the MATLAB environment. The numerical results are compared with the exact results of the problem, so that their accuracy can be assessed; these exact results are presented in the Appendix. For the numerical application, the following values are used: Cd0= 0.015,Cd2= 0.042,ρ= 0.5ρ0,ρ0= 1.225 kg/m3,V= 200 m/s,c= 5 ·10−5s/m,S= 150 m2,g= 9.8 m/s2,mF= 25000 kg in the case of given fuel load, and xf= 2500 km in the case of given range. For the initial mass distributions, the nominal values chosen for mean and width are ¯m0= 81633 kg and δm= 5000 kg, which yields a typical deviation σ[m0] = pVar[m0] = δm √3= 2887 kg; and for the gamma distribution the nominal value k= 8.5is considered. A parametric study as function of δmand kis also presented. For the nominal values, the two initial mass distributions are shown in Fig. 3. For the uniform distribution the values of m0with nonzero probability are in the interval [m0−δm, m0+δm] = [76633,86633] in kg, and for the gamma distribution they are in [M0,∞) = [73217,∞)in kg. For the GPC method the number of terms used in the expansions is P= 3, which turns out to be enough to obtain a good representation of m. In the computation of the distribution functions, the number of discretization points considered is n= 1000, which has proven to be good enough. All the integrations have been performed using the Matlab environment. In Section VII.A the GPC method is applied to obtain the evolution of mass mean and variance. The distribution function of the mass and its evolution in time are analyzed in Section VII.B. The distribution function of the flight time in the case of given fuel load is studied in Section VII.C, and that of the fuel consumption in the case of given range in Section VII.D. 16 7 7.5 8 8.5 9 9.5 x 1 0 4 0 0.5 1 1.5 x 10−4 fm0 m0[kg] Gamma distribution Uniform distribution Fig. 3 Probability density functions of the initial mass ( ¯m0= 81633 kg,δm= 5000 kg, and k= 8.5). A. Mass mean and variance 1. Uniform distribution of the initial mass To find the mean and variance using GPC, the value P= 3 is chosen in the GPC expansion of the mass (Eq. 12), which, as already mentioned, is enough to obtain a good representation of m. The coefficients of the GPC expansion are shown in Fig. 4. Note the fast decrease of their order of magnitude (six orders of magnitude from h0to h3). 0 4000 8000 12000 5 6 7 8 9x 104 h0(t) time (s) 0 4000 8000 12000 4000 4200 4400 4600 4800 5000 h1(t) time (s) 0 4000 8000 12000 −20 −15 −10 −5 0 h2(t) time (s) 0 4000 8000 12000 0 0.02 0.04 0.06 0.08 h3(t) time (s) Fig. 4 GPC coefficients for the uniform distribution case. 17 The evolution of mean E[m(t;m0)] and typical deviation σ[m(t;m0)] = pVar[m(t;m0)] is shown in Fig. 5. Selected values of mean and typical deviation are given in Table 1. The difference between the GPC solution and the analytical solution (Eqs. (59) and (60)) of mean and variance is negligible; the absolute error is less than 10−4for the mean and less than 2·10−3for the typical deviation. Thus a low-order GPC expansion, which is very fast to compute, is enough to capture well the mean and variance evolution. While the fact that the mean mass decreases with time is to be expected (since fuel mass in consumed), it is remarkable that the standard deviation of the mass also decreases with time. Thus, the dispersion of the distribution function and, therefore, the uncertainty decreases with time. This result can be explained by noting that the larger the aircraft mass, the larger its rate of decrease (which is given at each instant by A+Bm2). Thus, if one computes the solution m(t)given by Eq. (5) for m0= ¯m0±δm, say, m+(t) = m(t; ¯m0+δm)and m−(t) = m(t; ¯m0−δm), the distance ∆m(t) = m+(t)−m−(t)decreases with time: for example, at t= 0 one has ∆m= 2δm= 10000 kg, and at t=1.2 104s, ∆m= 8321 kg. 0 2000 4000 6000 8000 10000 12000 5 6 7 8 9x 104 time (s) E[m] 0 2000 4000 6000 8000 10000 12000 2200 2400 2600 2800 3000 time (s) σ[m] Fig. 5 Evolution of mass mean and typical deviation for the uniform distribution case. 2. Gamma distribution of the initial mass As in the uniform distribution case, to find the mean and variance using GPC, choosing P= 3 in the GPC expansion of the mass (Eq. 12) is good enough. The coefficients of the GPC expansion are shown in Fig. 6. Note again the fast 18 Table 1 Values of mass mean and typical deviation at selected times for the uniform distribution case Time (s) E[m(t;m0)] (kg) σ[m(t;m0)] (kg) 2·10377485 2787 4·10373477 2696 6·10369596 2613 8·10365831 2536 10462175 2467 1.2·10458616 2402 decrease of their order of magnitude (seven orders of magnitude from h0to h3). 0 4000 8000 12000 5 6 7 8 9x 104 h0(t) time (s) 0 4000 8000 12000 −1000 −950 −900 −850 −800 h1(t) time (s) 0 4000 8000 12000 −2.5 −2 −1.5 −1 −0.5 0 h2(t) time (s) 0 4000 8000 12000 −8 −6 −4 −2 0x 10−3 h3(t) time (s) Fig. 6 GPC coefficients for the gamma distribution case (k= 8.5). The evolution of mean and typical deviation is shown in Fig. 7. Selected values of mean and typical deviation are given in Table 2. As before, the difference between the GPC solution and the analytical solution (Eqs. (63) and (64)) of mean and variance is negligible; the absolute error is less than 4·10−5for the mean and less than 2·10−3for the typical deviation. Again, both the mean and the standard deviation decrease with time. Note that the plots and values are very similar to the ones obtained with the uniform distribution. Thus, the results show that the evolution of mean and standard deviation is very weakly affected by the specific distribution function chosen for the initial mass (at least for the two cases studied). 19 0 2000 4000 6000 8000 10000 12000 5 6 7 8 9x 104 time (s) E[m] 0 2000 4000 6000 8000 10000 12000 2200 2400 2600 2800 3000 time (s) σ[m] Fig. 7 Evolution of mass mean and typical deviation for the gamma distribution case (k= 8.5). Table 2 Values of mass mean and typical deviation at selected times for the gamma distribution case (k= 8.5) Time (s) E[m(t;m0)] (kg) σ[m(t;m0)] (kg) 2·10377485 2786 4·10373477 2695 6·10369596 2610 8·10365831 2533 10462175 2462 1.2·10458616 2397 B. Distribution function of the mass 1. Uniform distribution of the initial mass The mass distribution is represented at several time instants in Fig. 8. Both approximate methods developed in Section V to approximate Eq. (23) show excellent agreement with the exact analytical results (Eq. 70) and are indistinguishable from them. The results in Fig. 8 show that as time increases (and mdecreases), the width of the distribution function decreases, while the probability density increases. Thus, uncertainty decreases with time (as it was seen in Fig. 5). Note also that the uniform shape is approximately maintained. 20 5.5 6 6.5 7 7.5 8 8.5 9 x 1 0 4 0 0.4 0.8 1.2 x 10−4 m [ kg ] fmt=12800 s t=6400 s t=0 s Fig. 8 Mass distribution at several time instants for the uniform distribution case. 2. Gamma distribution of the initial mass In this case, the mass distribution is represented at several time instants in Fig. 9. Again, both numerical methods developed in Section V to approximate Eq. (23) show excellent agreement with the exact analytical results (Eq. 72) and are indistinguishable from them. As in Fig. 8, Fig. 9 shows that uncertainty decreases as time increases. Also, the shape of the distribution function is approximately of gamma type at all times. 5 6 7 8 9 10 x 1 0 4 0 1 2 x 10−4 m [ kg ] fm t= 12800 s t= 6400 s t=0s Fig. 9 Mass distribution at several time instants for the gamma distribution case (k= 8.5). C. Distribution function of the flight time 1. Uniform distribution of the initial mass The distribution function of the flight time is represented in Fig. 10. Note that it looks approximately uniform, similarly to the initial mass distribution. As in the computation of the mass distribution function, both approximate methods developed in Section V to approximate Eq. (32) show excellent agreement with the exact analytical result (Eq. 77). 21 The values of tfwith nonzero probability are those in the interval [T1, T2] = [12625,13664] in seconds, for values of m0with nonzero probability in [m0−δm, m0+δm] = [76633,86633] in kilograms, where, as shown in the Appendix, T1=1 √AB arctan √ABmF A+B( ¯m0+δm)( ¯m0+δm−mF)!(49) T2=1 √AB arctan √ABmF A+B( ¯m0−δm)( ¯m0−δm−mF)!(50) 1.2 1.3 1.4 x 1 0 4 0 0.4 0.8 1.2 1.6 x 10−3 ftf tf(s) a b Fig. 10 Distribution functions of the flight time: a) uniform distribution case; b) gamma distribution case (k= 8.5). The mean and the typical deviation of the flight time are obtained using the distribution function, computed numerically from E[tf] = Z∞ 0 tfftf(tf)dtf(51) (σ[tf])2=Z∞ 0 t2 fftf(tf)dtf−(E[tf])2(52) The results are given in Table 3. Table 3 Computed values of E[tf]and σ[tf]for the uniform distribution case Exact Method 1 Method 2 E[tf](s) 13144 13144 13146 σ[tf](s) 299.8 296.9 290.1 22 Now the effect of δmon the results is analyzed. Values of σ[tf]for different values of δm(obtained using the exact solution) are given in Fig. 11, where it is seen that there is a proportionality between the two parameters. The values of E[tf]are not significantly affected by changing δm. 2000 3000 4000 5000 6000 7000 8000 100 200 300 400 500 σ[tf](s) δm(kg) Fig. 11 Typical deviation of the flight time vs. δmin the uniform distribution case. 2. Gamma distribution of the initial mass The distribution function in this case is represented in Fig. 10. Note that this distribution function is somewhat different from a gamma distribution, because the values of tfwith nonzero probability are those in the finite interval (0, T] = (0,14019] in seconds, for values of m0with nonzero probability in [M0,∞) = [73217,∞)in kilograms, where, as shown in the Appendix, T=1 √AB arctan √ABmF A+BM0(M0−mF)!(53) Moreover, since tfdecreases when m0increases, the bell of the distribution is sort of inverted (with respect to the bell of the initial mass distribution). As in the computation of the mass distribution function, both numerical methods developed in Section V to approximate Eq. (32) show excellent agreement with the exact analytical result (Eq. 79). Again, the mean and the typical deviation are computed numerically from Eqs. (51) and (52) using the distribution function. The results are given in Table 4. Table 4 Computed values of E[tf]and σ[tf]for the gamma distribution case (k= 8.5) Exact Method 1 Method 2 E[tf](s) 13144 13144 13143 σ[tf](s) 299.3 301.6 303.9 23 Note that these results of mean and standard deviation are very close to the ones obtained before for the uniform distribution (especially for the exact distribution functions), showing again that the initial mass distribution chosen affects the results very weakly. Now the effect of kon the results is analyzed. Values of σ[tf]for different values of kand δm(obtained using the exact solution) are given in Fig. 12, where it is seen that there is no significant effect from changing kand, as in the uniform distribution case, there is a proportionality between the values of σ[tf]and δm. The values of E[tf]are not significantly affected by changing kor δm. 0 10 20 30 40 50 60 0 100 200 300 400 500 σ[tf](s) k δm= 3500 kg δm= 2000 kg δm= 5000 kg δm= 6500 kg δm= 8000 kg Fig. 12 Typical deviation of the flight time vs. k, for different values of δm, in the gamma distribution case. D. Distribution function of the fuel consumption 1. Uniform distribution of the initial mass The distribution function of the fuel consumption is represented in Fig. 13. Note that it looks approximately uniform, similarly to the initial mass distribution, although smaller values of mFshow a slightly higher probability. As before, both approximate methods developed in Section V to approximate Eq. (40) show excellent agreement with the exact analytical results (Eq. 83). 24 2.3 2.4 2.5 2.6 2.7 x 1 0 4 0 0.2 0.4 0.6 0.8 x 10−3 fmF mF[kg] b a Fig. 13 Distribution functions of the fuel consumption: a) uniform distribution case; b) gamma distribution case (k= 8.5). The values of mFwith nonzero probability are those in the interval [M1, M2] = [23043,24775] in kilograms, for values of m0with nonzero probability in [m0−δm, m0+δm] = [76633,86633] in kilograms, where, as shown in the Appendix, M1=( ¯m0−δm)2+A Btan 1 V√ABxf rA B+ ( ¯m0−δm) tan 1 V√ABxf(54) M2=( ¯m0+δm)2+A Btan 1 V√ABxf rA B+ ( ¯m0+δm) tan 1 V√ABxf(55) As for the flight time, the mean and the typical deviation of the fuel consumption are obtained using the distribution function, computed numerically from E[mF] = Z∞ 0 mFfmF(mF)dmF(56) (σ[mF])2=Z∞ 0 m2 FfmF(mF)dmF−(E[mF])2(57) The results are given in Table 5. Table 5 Computed values of E[mF]and σ[mF]for the uniform distribution case Exact Method 1 Method 2 E[mF](kg) 23892 23892 23892 σ[mF](kg) 499.96 499.81 500.04 25 xs−1e−x1+(s−1) 1 x+ (s−1)(s−2) 1 x2+. . ., valid for x→ ∞ (see Ref. [18]). It follows that E[m(t;m0)] ≈c2M0−c2 2c1 c1M0+c2 +kc2 2c2 1+ 1 c1M0+c2       δm √3k c1M0+c2   −(k+ 1)c1    δm √3k c1M0+c2    2 +(k+ 1)(k+ 2)c2 1    δm √3k c1M0+c2    3 −(k+ 1)(k+ 2)(k+ 3)c3 1    δm √3k c1M0+c2    4 +(k+ 1)(k+ 2)(k+ 3)(k+ 4)c4 1    δm √3k c1M0+c2    5 #(65) V ar[m(t;m0)] ≈δ2 m 3 (c2 1+ 1)2c4 2 (c1M0+c2)4"1−4(k+ 1)     c1 δm √3k c1M0+c2    + 2(k+ 1)(5k+ 9)     c1 δm √3k c1M0+c2    2 −4(1 + k)(2 + k)(12 + 5k)    c1 δm √3k c1M0+c2    3 +(1 + k)(2 + k)(300 + 7k(29 + 5k))     c1 δm √3k c1M0+c2    4 #(66) B. Distribution function of the mass To compute the distribution function of the mass, note that ∂m ∂m0 can be exactly computed from Eq. (5) as ∂m ∂m0 (t;m0) = A B 1 + tan2√ABt rA B+m0tan √ABt!2=c2 2 1 + c2 1(t) (c2+m0c1(t))2(67) Also from Eq. (5), m0can be written in terms of m(t)as follows m0=c2 m(t) + c2c1(t) c2−c1(t)m(t)(68) 32 Thus ∂m ∂m0 in terms of mis written as ∂m ∂m0 (t;m) = (c2−c1(t)m)2 c2 2(1 + c2 1(t)) (69) 1. Uniform distribution of the initial mass Since fm0=1 2δm , using Eqs. (23) and (69), the exact distribution function of the mass as a function of time is fm(m, t) = c2 2(1 + c2 1(t)) 2δm(c2−c1(t)m)2(70) if m∈c2( ¯m0−δm)−c2 2c1(t) c2+c1(t)( ¯m0−δm),c2( ¯m0+δm)−c2 2c1(t) c2+c1(t)( ¯m0+δm), and zero otherwise. The limit points in the interval have been found from Eq. (5) evaluated at the limit points in the initial mass distribution ( ¯m0−δmand ¯m0+δm). 2. Gamma distribution of the initial mass In this case the distribution function fm0given by Eq. (11) has to be written in terms of m(t)using Eq. (68), as follows fm0(m, t) = c2m+c2 2c1(t) c2−c1(t)m−M0k−1e − √3k δm  c2m+c2 2c1(t) c2−c1(t)m−M0  δm √3kk Γ(k) (71) Then, using Eqs. (23), (69), and (71), the exact distribution function of the mass as a function of time is fm(m, t) = c2m+c2 2c1(t) c2−c1(t)m−M0k−1e − √3k δm  c2m+c2 2c1(t) c2−c1(t)m−M0  δm √3kk Γ(k) c2 2(1 + c2 1(t)) (c2−c1(t)m)2(72) if m≥c2 M0−c2c1(t) c2+M0c1(t), and zero otherwise. The lower limit is found evaluating Eq. (5) at m0=M0. 33 C. Distribution function of the flight time To compute Eq. (32), using Eq. (6), the values of ϕand its inverse can be explicitly obtained as tf=ϕ(m0) = 1 √AB arctan √ABmF A+Bm0(m0−mF)!(73) m0=ϕ−1(tf) = mF 2+v u u u u t m2 F 4−A B+rA BmF tan √ABtf=mF 2+ Φ(tf)(74) where Φ(tf) = v u u u u t m2 F 4−A B+rA BmF tan √ABtfis defined to simplify the expressions. Also ϕ0(m0)is given by ϕ0(m0) = −BmF(2m0−mF) (A+Bm0(m0−mF))2+ABm2 F (75) Hence |ϕ0(ϕ−1(tf))|= 2 sin2√ABtf AmF Φ(tf)(76) These results are now used to derive an explicit expression for ftf, for the two initial mass distributions under consideration. 1. Uniform distribution of the initial mass From Eq. (32), using fm0=1 2δm and Eq. (76), the resulting expression for the exact distribution function of the flight time is ftf(tf) = AmF 4δmsin2√ABtfΦ(tf) (77) if tf∈[T1, T2], and zero otherwise, where the endpoints of this interval are found evaluating Eq. (73) at the endpoints of the uniform distribution of m0(namely, ¯m0−δmand ¯m0+δm), and are given by T1=1 √AB arctan √ABmF A+B( ¯m0+δm)( ¯m0+δm−mF)!and T2=1 √AB arctan √ABmF A+B( ¯m0−δm)( ¯m0−δm−mF)!. 34 2. Gamma distribution of the initial mass To find ftfnow, the distribution function fm0for the gamma case (Eq. 11) has to be written in terms of tfusing Eq. (74), as follows fm0(tf) = mF 2+ Φ(tf)−M0k−1e − √3k δmmF 2+Φ(tf)−M0 δm √3kk Γ(k) (78) Then, from Eq. (32), using Eqs. (78) and (76), the resulting expression for the exact distribution function of the flight time is ftf(tf) = mF 2+ Φ(tf)−M0k−1AmF 2 sin2√ABtf e − √3k δmmF 2+Φ(tf)−M0 δm √3kk Γ(k)Φ(tf) (79) for tf∈ 0,1 √AB arctan √ABmF A+BM0(M0−mF)!#, and zero otherwise. The upper limit value is found evaluating Eq. (73) at m0=M0. D. Distribution function of the fuel consumption To compute Eq. (40), the inverse of ψ(m0)is necessary. For that, one has to solve for m0in Eq. (9), finding m0=ψ−1(mF) = mF 2+sm2 F 4+mFc2 c3−c2 2=mF 2+ Ψ(mF)(80) where Ψ(mF) = sm2 F 4+mFc2 c3−c2 2has been defined. Also, taking the derivative with respect to m0in Eq. (9), the value of ψ0(m0)is found as ψ0(m0) = c3 2c2m0+m2 0−c2 2c3 (c2+m0c3)2(81) For Eq. (40), it is necessary to explicitly compute |ψ0(ψ−1(mF))|using Eqs. (80) and (81), finding |ψ0(ψ−1(mF))|=c3mF 2+ Ψ(mF)2c3+ 2 mF 2+ Ψ(mF)c2−c3c2 2 hc2+c3mF 2+ Ψ(mF)i2(82) This result is now used to derive an explicit expression for fmF, for the two initial mass distributions under consideration. 35 1. Uniform distribution of the initial mass From Eq. (40), using fm0=1 2δm and Eq. (76), the resulting expression for the exact distribution function of the flight time is fmF(mF) = 1 2δmc3hc2+c3mF 2+ Ψ(mF)i2 mF 2+ Ψ(mF)2c3+ 2 mF 2+ Ψ(mF)c2−c3c2 2 (83) if mF∈"( ¯m0−δm)2+c2 2c3 c2+ ( ¯m0−δm)c3 ,( ¯m0+δm)2+c2 2c3 c2+ ( ¯m0+δm)c3#, and zero otherwise. The endpoints of this interval are found evaluating Eq. (9) at the endpoints of the uniform distribution of m0(¯m0−δmand ¯m0+δm). 2. Gamma distribution of the initial mass To find fmFnow, the distribution function fm0for the gamma case (Eq. 11) has to be written in terms of mFusing Eq. (80), as follows fm0(mF) = mF 2+ Ψ(mF)−M0k−1e −mF 2+Ψ(mF)−M0√3k δm δm √3kk Γ(k) (84) Then, from Eq. (40), using Eqs. (84) and (82), the resulting expression for the exact distribution function of the fuel consumption is fmF(mF) = e −mF 2+Ψ(mF)−M0√3k δm δm √3kk Γ(k)c3 mF 2+ Ψ(mF)−M0k−1hc2+c3mF 2+ Ψ(mF)i2 mF 2+ Ψ(mF)2c3+ 2 mF 2+ Ψ(mF)c2−c3c2 2 (85) for mF≥M2 0+c2 2c3 c2+M0c3 , and zero otherwise. The lower limit is found evaluating Eq. (9) at m0=M0. 36