Optimization of Speed Profile with RTA Constraints under Wind Uncertainty

Four-dimensional (4D) trajectory is considered to be one of the effective means to reduce the environmental impact of aviation while increasing capacity and safety. This paper proposes an approach to optimize cruise speed profile subject to wind uncertainty, aiming to reduce the fuel burn complying with the Required Time of Arrival (RTA) constraints. The approach is based on a probabilistic framework, and the uncertainty propagation is analyzed using a Probability Transformation Method, and the probability distributions of arrival time and fuel consumption are determined. In addition, from an airborne operation perspective, transition profiles need to be considered in the reference speed optimization problem, aiming to improve the rationale of the reference trajectory. Numerical simulations are presented, and the results demonstrate that the speed profiles optimized by this method are able to meet the RTA constraints in the presence of wind uncertainty, with an average reduction of 7.04% in fuel consumption compared with that of the flight data.


Introduction
Air traffic increasing will bring about some new challenges to the airspace, such as increasing airspace tension and environmental impact. e Trajectory-Based Operation (TBO) [1,2] is regarded as a key concept to ensure flight safety and reduce the flight cost. It requires the aircraft to fly following the given four-dimensional (4D) trajectory strictly, including the space of 3D trajectory and the Required Time of Arrival (RTA) at each waypoint. is innovative technology can greatly enhance the predictability of trajectory and the safety of flight.
It is well known that civil aviation flight is a highly planned task, which requires each aircraft to carry out a predetermined flight plan. It is an important foundation for air traffic safety, especially in 4D flying. However, there are various uncertainties in the actual flight environment, which impose great challenges in the implementation of 4D-TBO. ese sources of uncertainty include weather, navigation accuracy, and operational procedures [3]. Hence, several relevant research studies have been carried out. In [4], a 4D trajectory in descent phase, considering meteorological uncertainty and reducing fuel consumption, is proposed. In [5], the error in conflict prediction caused by various uncertainties is considered, and optimal conflict-free trajectory is generated by solving the stochastic conflict resolution problem. In [6], an approach is proposed to optimize the 4D trajectory of each aircraft so as to minimize the interaction between trajectories, and the uncertainties of aircraft position and arrival time are taken into account.
In 4D flight, the RTA function is a significant ability to manage the 4D trajectory, and a managed speed mode is proposed to achieve time-based operation. For instance, a flight deck real-time simulation is held at the EURO-CONTROL Experimental Centre [7] to investigate 4D trajectory management in the cruise and initial descent phase, and the result shows that pilots prefer to control speed to handle time constraints. us, the optimization of speed is significant during cruise. A research is also performed at the Research Laboratory in Active Controls, Avionics and Aeroservoelasticity (LARCASE), aiming to provide the optimal flight path for the given lateral flight profile and RTA constraints [8]. In particular, due to the constraints of time dimension, wind uncertainty imposes a great influence on the performance of 4D-TBO, which includes the following two aspects. (1) e wind forecast error will directly affect the estimation of aircraft ground speed and further affect the prediction of Estimated Time of Arrival (ETA). In order to fulfill the RTA constraint, airspeed should be constantly adjusted according to the wind speed. (2) Fuel consumption will be directly affected by airspeed, and fuel economy is one of the important indexes to evaluate the performance of 4D trajectory. From the above two aspects of 4D performance, it can be seen that the 4D flight is more sensitive to the weather. In the presence of additional time constrains, the 4D flight plan reserves less adjustable margin to meet the fuel economy requirement. is will make it more difficult for the crew to fly 4D in the face of wind uncertainty. en, in this limited adjustment margin, how to manage the speed to adapt to the wind uncertainty is the first core problem to be solved in the implementation in 4D.
Researches on RTA function in trajectory planning have been addressed by many authors with different technologies and methods. In [9], based on the dynamic programming method, the 4D trajectory with multiple RTA constraints is optimized, and more sophisticated time windows such as at or before or after constraints could be solved. In [10], the minimum-fuel trajectory with RTA problem arising in air Air Traffic Management is analyzed, and pseudospectral collocation method is adopted for solving optimal control problems. In [11,12], a methodology is proposed to create a 4D reference trajectory by using meta-heuristic algorithms, finding an optimal combination of Mach numbers at different waypoints to reduce fuel burn. However, in most of literatures, the flight environment is treated as an "ideal" deterministic one. e wind experienced by the aircraft is considered to be known before take-off (fully in accordance with the forecast), but when there exist deviations between the forecast and the actual dynamic weather outcome, it is difficult to ensure the arrival time and economic requirements and even influence the flight safety. erefore, it is necessary to take uncertainty into account and put forward new strategies to reduce the impact of uncertainty on trajectory planning. To the best of the authors' knowledge, few research efforts have been devoted to this issue.
Besides the research presented in the context of Air Traffic Management (ATM), avionic companies are developing the RTA function to add it to the new generations of Flight Management System (FMS) [13,14]. e fulfillment of 4D trajectory in the future will much rely on the improvement of airborne capability. e target of the flight plan is to provide a more rational expectation for FMS. In that case, the reference trajectory of a flight plan submitted to FMS should be usable from operational airborne consideration, reducing the modification of the flight plan during flight to the maximum extent. Since the aircraft cannot adjust airspeed in a stepwise manner, if the stepwise speed reference trajectory proposed in [11] is implemented for airborne use without considering the transitions when speed variations happen, the aircraft speed cannot show any stepwise discontinuity in practice. For another reason, the acceleration or deceleration will consume additional fuel, thereby changing speed frequently will increase the fuel consumption. Hence, the optimized trajectory without considering such perspective might not be the most economical trajectory. e above analyses suggest that an improvement should be made by smoothing speed guidance, and the operational capability should be considered in future trajectory planning.
In this work, we address the problem of optimizing the reference speed profile subject to wind uncertainty, while satisfying multiple RTA constraints. To this purpose, the wind forecast model is included, which is the main source of uncertainty in trajectory planning, and the uncertainty propagation regarding 4D trajectory is analyzed. Using a numerical method called Probability Transform Method (PTM), the Probability Density Functions (PDFs) of the estimated arrival time and fuel consumption are obtained.
is way, on the basis of probability, the Mach speed permitting to fulfill multiple RTA constraints with predefined probability is determined. Moreover, for airborne operation and consideration of additional fuel consumed by speed transition parts, another interest about how to optimize the transition profile for speed variation is explored as well. Taking passenger comfort into account, the jerk and acceleration are within the acceptable limits as well. Finally, several domestic flights in China are taken as simulation verification, and some conclusions are drawn.

Uncertainty Propagation Analysis
e stochastic nature of the wind is the main source of uncertainty affecting the aircraft motion. In this section, the uncertainty model of wind speed, as an uncertainty source involved in this research, is presented. en, the uncertain propagation characteristics of 4D trajectory are here detailed by analyzing how the wind affects the flight time and fuel consumption.

Wind Uncertainty Model.
In general, the atmospheric wind field can be expressed with three wind components along three perpendicular directions, where each component is a function of both position and time. In this paper, only the horizontal trajectory is concerned, and, thus, the wind vector is converted to the wind speed along the track W(P, t).
e wind uncertainty can be illustrated by two components: a deterministic term W(P, t + τ | t) representing nominal wind or mean wind and a stochastic term ξ(P, t + τ | t) accounting for the mismatch between forecast and the actual wind, which can be expressed as where (P, t + τ | t) is defined as prediction over a look-ahead interval τ at time t for position P.
From the meteorological forecast, the required predicted wind information can be analyzed to provide prior information as the deterministic term W(P, t + τ | t) in (1), and the stochastic part ξ(P, t + τ | t) can be obtained by statistical analysis of historical data.
In order to statistically analyze the wind forecast error, another important consideration for the wind field is the spatial and temporal correlation structure, which will determine the impact of wind uncertainty on the aircraft. It is assumed that the wind field is quasi-stationary and quasihomogeneous, and the wind forecast error broadcast over similar time in the same season under the similar weather condition is further assumed; in other words, the sample with spatial and temporal correlation follows the same probabilistic distribution. In literature [15], typically spatial and temporal correlation ranges are, respectively, s * � 135 nm and τ * � 120 min. us, s i and τ i are within the look-ahead time τ * and spatial range s * , respectively; the average wind speed is approximately piecewise constant, as shown in the following equation: Eventually, the probabilistic wind model in the form of grid is discretized in time and space, as shown in Figure 1. e horizontal space is divided into meteorological grids by the distance of spatial correlation, and the predicted wind of each grid conforms to the same probability distribution. Gaussian distribution is often used to model wind forecast uncertainties [16], ξ(P, t | t 0 ) ∼ N(0, σ 2 ). After superimposing the forecast value with the uncertain value, the mean value of wind speed W i is the wind forecast value W i , and finally the distribution of wind is expressed as

Uncertainty in the Cruise Trajectory.
We focus on the cruise phase in this research, for the reason that wind uncertainty has a large impact on such phase due to its large proportion for the entire route (in particular, for long-haul flight). Besides, the most significant fuel reduction can be obtained. e cruise phase is given in multiple waypoints O i according to the flight plan; then, enroute is divided into multiple segments, the 4D waypoint of the path by Air Navigation Service provider (ANSP) constraints; these constraints, referred to as RTA, request the aircraft, reaching the waypoint within predetermined time error under the allowable performance.
Only the along-track wind is considered, ignoring the vertical component. e along-track wind component V lw can be calculated from (3), and it takes a positive value for tailwind and a negative one for headwind: As shown in Figure 2, the actual ground speed V g is calculated from the airspeed V a and the average along-track wind V lw : V lw is wind speed along the track, V w is wind speed, and λ is the angle of wind direction and head direction. According to Section 2.1, the forecast wind is modeled as a random variable, and the uncertain along-track wind is expressed as e flight time t i of segment i is expressed as where L i is the length of cruise segment i. Note from (5) and (6) that the flight time t is subjected to the along-track wind, which shows random uncertainty. en, the uncertainty of flight time t can be derived. erefore, due to time constraints, the airspeed V a is no longer constant and will be used as the optimization decision variable in this study, aiming to reduce the impact of wind uncertainty on RTA performance.
Furthermore, the fuel consumption in cruise flight is studied. e fuel flow rates in cruise ff cr are estimated by submatch fuel consumption method; the details are omitted here, referring to [17]. In this method, in general, the fuel flow rate is determined by the following input variables: pressure altitude, airspeed, longitudinal acceleration, gross weight, vertical speed, flight path angle, and total temperature. In this research, it is assumed that the altitude is constant and the flight path angle and vertical speed are zero. Also, the uncertainty of temperature is omitted and assumed to be perfectly predicted. Finally, airspeed, longitudinal acceleration, and gross weight are the variables that will influence the fuel flow rates, among which the airspeed is the decision variable; then the fuel flow rates can be estimated. Next, the total fuel consumption for all segments can be expressed as the product of fuel flow and flight time:  Mathematical Problems in Engineering 3 As shown in (7), it can be seen that, due to the uncertainty of flight time, the total fuel consumption exhibits uncertainty as well.

Speed Profile Optimization Method
In this section, the optimization of the speed profile considering the transition parts is proposed. e speed profile includes two parts: stepwise Mach number selection and the transitions between two different speeds. Based on the uncertainty analysis in the previous section, the stochastic optimization model for speed selection is developed to address the uncertainty issue. A method in a probabilistic framework is then proposed to solve this problem. Afterwards, the transition profile of speed variation is optimized considering passenger comfort. Integrating the two parts, the optimal speed profile is iteratively adjusted until the objective of minimum fuel consumption is achieved. In what follows, we will first discuss how to optimize the speed profile in the face of the random uncertainty of the wind environment.

Optimization Model.
A single RTA may lead to trajectory prediction error, especially when there is wind uncertainty, and the error will further increase. To improve the predictability of trajectories, multiple RTA can be specified.
Given a certain aircraft type, flight altitude, and flight environment in the traditional 3D trajectory, there exists a most economical cruise speed, and the economic requirement of minimum fuel consumption can be achieved by maintaining this speed. However, in addition to economic requirements, compliance with RTAs is another essential requirement for the 4D trajectory optimization task. Obviously, flying at a constant speed that satisfies the RTA constraint may deviate from the economic cruise speed, and there is a contradiction between the RTA constraint and the fuel economy. In addition, due to the dynamic uncertainty of weather, it is difficult to guarantee the constant airspeed required by the RTA. e varying airspeed approach is adopted in this paper. For the convenience of airborne operation and separation management for ATM, virtual waypoints P k (j � 1, 2, . . ., n) placed are set along the route with a fixed interval l (defined as a virtual segment), as shown in Figure 3. Airspeed V a K is assumed to be constant in each virtual segment, which can not only ensure the RTA constraints but also further improve the fuel economy if the airspeed combination is optimized. To ensure sufficient speed adjustment margin, the length of the shortest segment should be longer than twice of the virtual waypoint interval; that is, L min ≥ 2l, and, in each virtual segment k, the average along-track wind V lw k is calculated. en the cruise route is divided into K subsegments by virtual waypoints and the RTA waypoints together. e arrival time t ETA j of RTA waypoint j is as follows: where k is the index of virtual segment and RTA waypoint, j is the index of RTA waypoint, D j is the distance of RTA waypoint j from TOC point, a j is the number of virtual segments contained in D j , and a j � floor(D j /l) + 1.
For solving the stochastic optimization problem to minimize the fuel consumption, the following stochastic cost function J expressed in (9) is minimized: where fuel k is the fuel consumed by virtual segment k, ff k is the fuel flow rate, t k is the flight time, and K is the total number of virtual segments. Because the RTA constraint involves flight safety, it is represented as one of the constraints in the optimization process. Because of the uncertainty of arrival time, the time window is one of the important definitions in 4D trajectory concept, which provides a certain degree of error tolerance. Time error t Δ , defined as the difference between ETA and RTA, should be limited within the predefined window T: Since time error t Δ depends additively on wind disturbance ξ, chance constraint is introduced to describe RTA constraints. e constraint may not be feasible, allowing the decision to violate the constraint to some extent before a random variable is observed. e arrival time of each waypoint has to be satisfied with a given probability α j of arriving in a time window: where t RTA j and t ETA j are time of RTA and ETA at waypoint j, respectively, and T j is the time window. Mach number is adopted during the cruise, defined as Mach � cV a , where c is the speed of sound at the given flight pressure altitude, and its value should be kept within aircraft performance constraints: 3.1.2. Solution of the Optimization Problem. We will discuss how to solve the proposed model. ere are two uncertain variables contained in the speed selection optimization problem: the arrival time and fuel consumption, both of which indicate 4D trajectory performance. e optimization is based on a probabilistic method, which is described in this section. To reduce the search space of the problem, a genetic algorithm is applied to solve the combinational optimization problem in dynamic uncertain search space.
Probabilistic transformation method (PTM) [18] is used to estimate the probability distribution of the ETA. Given the random input variable y with PDF p Y (y), z is the response variable z � f (y), and its inverse transformation exists in the domain of y, that is, f −1 (·) � h(·), as shown in Figure 4; PDF of system output z is calculated by the product of the input PDF and Jacobian matrix of the inverse function: For the problem studied in this paper, the relationship between flight time and wind speed in subsegment i is as follows: e inverse transformation is as follows: and, thus, the PDF of the flight time of subsegment can be derived by PTM: For independent random variables X and Y, with the function relationship Z � X + Y, the PDF of Z can be obtained by the convolution formula of probability: Assume that one segment is divided into q subsegments, and the PDF of the flight time of the segment p t (t) can be computed from a set of convolutions: e aircraft can meet the RTA constraint by arriving within the time window with a probability of no less than α(0 < α < 1), as shown in Figure 5. In general, due to the existence of aircraft performance constraints and the possibility of adverse weather (such as headwind or greater uncertainty), usually the probability α < 1, the probability is expressed as where t early and t late are the flight times corresponding to the allowable earliest and latest flight time, respectively. Given a certain airspeed, the fuel flow is a determinate variable.
erefore, the PDF of fuel consumption in the subsegment is expressed as the product of fuel flow and the PDF of flight time: From the probability convolution (19), the PDF p m fuel (m fuel ) of total fuel consumption is as follows: p m fuel m fuel � p m fuel 1 m fuel 1 * p m fuel 2 m fuel 2 * · · · * p m fuel K m fuel K .

(22)
Eventually, given the wind PDF for each grid, the total fuel can be calculated, as shown in Figure 6. e expectation and variance of total consumed fuel can be written as e airspeed V a k (k � 1, 2, . . ., n), which is the variable to optimize, will determine the fuel flow rates and distribution of uncertain variable, and then genetic algorithm (GA) is used to search optimal speed combination. e details are described as follows: the horizontal trajectory of the cruise phase is divided into n + 1 virtual segments, and a gene on each chromosome is defined as the Mach number corresponding to the virtual segment. erefore, the population in this paper is a random Mach number combination, encoded by floating-point numbers, and a selection on roulette is carried out. For floating-point coding, an intermediate recombination is adopted; that is, the gene segments of any two individuals are exchanged. Sequence numbers of   Mathematical Problems in Engineering the exchange genes are generated randomly, and the new individual genes generated are located in the value range of the two parent genes. e probability-based method for searching the optimal airspeed combination satisfies the probability constraints of arrival time; we adopt the following procedures: Step 1: generate the chromosome of u � (Mach 1 , Step 2: compute Pr |t RTA − t ETA | ≤ T by (8), (14)-(20); if constraint (20) is satisfied, then go to Step 3; if not, return to Step 1 Step 3: use submatch method to calculate ff k Step 4: calculate p m fuel k (m fuel k ) using (21); then, PDF of total fuel consumption p m fuel (m fuel ) is calculated by (22) Step 5: the expectation of consumed fuel E[m fuel ] can be computed using (23) Step 6: search the optimal U opt � (V a 1 , V a 2 , . . . , V a K ) corresponding to the minimum E[m fuel ] by GA

Mathematical Formulation.
Since sporadic variation from the reference speed profile happens, the acceleration or deceleration should not be discarded for both operations might consume more fuel. erefore, an inspiration is motivated to focus on how to manage these speed transition phases so as to minimize fuel consumption caused by speed variation and to improve the maneuverability of the reference speed profile on operational level.
According to Newton's second law, the resultant force on the aircraft is the product of gross weight and acceleration. e acceleration, as the derivative of velocity, should exist. Nevertheless, it is not differentiable at the step change point for the reference profile. Taking acceleration for instance, the analysis of the deceleration process is similar to acceleration, for it is an acceleration process connecting two constant-speed phases; thus the acceleration will exhibit the following two or three phases: (i) Increasing from zero and decreasing to zero (ii) Increasing from zero, keeping constant, and decreasing to zero For circumstance (ii), suppose the jerk in each phase is kept constant: j 1 , 0, and j 2 , respectively. en, the acceleration can be deduced as follows: en, the acceleration of each phase is given by Furthermore, the velocity can be derived by integrating the acceleration equation: Combining (25)-(27), the detailed velocity can be expressed as . .
p m fuel (m fuel ) Figure 6: Computation of the PDF of the total consumed fuel for the trajectory with m segments.
From (26)-(28), it can be seen that the acceleration and velocity will be different with different values of j 1 , j 2 , t 1 , t 2 , t 3 , and t 4 , leading to diverse curves of acceleration and speed. For circumstance (i), the discrepancy in (ii) is only in default in the second phase of constant acceleration when t 2 � t 3 , so the expressions of acceleration and velocity are omitted at present in this paper. e fuel burn calculation method described in the previous section indicates that the fuel flow rate is closely related to the acceleration and airspeed. us, it is essential to find the appropriate values for these variables. Ultimately, the optimization problems are presented to search the optimal values of these six variables. e optimization objective is to minimize the fuel consumption of two segments: Considering the searching range, RTA requirements, aircraft performance, and passenger comfort, the inequality constraint sets are as follows: where t f is the RTA at the final point of the second segment: where t c is the spending time of flying the first segment at the initial speed, a max is the maximum acceleration of the aircraft, and t max is the maximum acceleration time. j min and j max are the minimum and maximum jerk of aircraft, respectively. e constraints in (30)-(32) are regional constraints, which will reduce the searching region, restricting the aircraft to stay within the flight envelope. e constraint in (33) stipulates that the acceleration of the aircraft is limited within predefined values considering passenger comfort. e constraint in (34) gives the maximum acceleration time. Also, the constraint in (35) guarantees passengers' comfort of the acceleration process, and the jerk is the derivative of the acceleration. Jerkiness is often used in engineering, especially in vehicle design and materials. When a vehicle accelerates, it will make passengers feel uncomfortable, which is not only caused by the acceleration but also related to the degree of jerkiness, serving as an index to judge passengers' discomfort. Here, considering the regulation of civil aviation industry, the jerk is ranging from 0.3 to 0.5 m/ s 3 in this paper. Equality constraints are as follows: e constraint in (37) shows that the final transition profile speed equals the target speed. e constraint in (38) indicates that the end acceleration of the transition profile equals zero. e constraint in (39) is the RTA constraint which guarantees that the Estimated Time of Arrival (ETA) equals RTA. Equality constraints are usually handled by converting them into inequality constraints as follows: where δ is a small positive value, J and R are the numbers of inequality and equality constraints, respectively, where J � 6 and R � 3 in this paper, and h r ( x → ) are r th equality constraints.

Constraint Handling Optimization Algorithm.
e constrained optimization problem in this paper is a nonlinear programming (NLP) problem. is problem has four variables, five inequality constraints, and three quality constraints, which indicate strong nonlinear characteristics. Multiple complicated constraints will lead the feasible region to occupy a tiny space in the whole search space, imposing more difficulty on searching. e method of solving the constraint optimization problem can be classified into two groups: (a) traditional methods such as Lagrange multiplier, Newton method, and quasi-Newton method and (b) heuristic algorithms, among which genetic algorithm (GA) is chosen in this paper as the effective way in solving such problems. Compared to the traditional algorithms, one big advantage of GA is that it does not rely on the initial value; moreover, the computation burden will not increase rapidly with an increasing number of variables.
In order to solve constrained optimization problems, penalty function methods have been commonly used. However, the toughest work of those methods is to find a penalty parameter, and an inappropriate one could lead to an unsatisfactory result. Here, the genetic algorithm incorporated with an efficient constraint handling method proposed by Deb is used in this paper to solve this constrained optimization problem, which does not require any penalty parameter. A detailed algorithmic description of this algorithm can be found in [19], which is only briefly described in this paper.
(1) Coding for individuals: the real-coded of floatingpoint coding instead of binary coding is adopted for all the variables, and the creation of an individual should respect the constraints in (28)-(33). A singlepoint crossover used in binary GAs may generate a large portion of infeasible solutions, since such Mathematical Problems in Engineering problem only has a tiny feasible region in the entire search space. Hence, floating-point representation of variables will be more efficient for our coding operation of GA.
(2) Evaluation: the fitness function F( x → ) shown in (41) is utilized to evaluate the optimization in terms of both objective function value and constraint violation degree: is objective function, f max is the objective function value of the worst feasible solution in the population, 〈g j ( x → )〉 is constraint violation, and m is the number of constraints. No penalty parameter is needed in this method, which is the major advantage of this method compared to other penalty functions. e feasibility of the solution is checked firstly. If the solution is infeasible, the objective function will never be computed.
(3) Crossover: simulated binary crossover (SBX) is used as the crossover operator. e spread of children solutions around parent solutions is controlled using a distribution index η c . (4) Selection: pairwise comparison in tournament selection is used as a selection method, leading the solution to come closer to the optimizer direction.

Numerical Results
In this section, an experiment of a typical domestic flight in China has been carried out. e aircraft type of the flight is Boeing 737-800, and the cruise phase is selected for case verification with the air distance of 1063.1 nm. e cruising altitude is 25,600 ft, the initial speed is 0.745 M, and the initial TOC weight is 151,60 lb. e actual horizontal trajectory drawn from the QAR data during the cruise is shown in Figure 7(a), and nine RTA waypoints (TOC is the initial point which is not included in the RTA waypoint) taken from the flight plan are marked in Figure 7(b), dividing the route into nine segments. In this experiment, every specific time corresponding to each waypoint given in the waypoint as RTA is looked up from the QAR data, as listed in Table 1. e requirements for this simulation are as follows: (1) Defined by BADA, the tighter time window Controlled Time of Arrival (CTA) [7] with a refined time window tolerance (±30 s) is adopted in this paper. (2) Frequent speed adjustment is not conducive to maintaining the separation between aircraft. In addition, computation load must also be considered. erefore, the length of the virtual segment is 30 nm, the Mach speed remains unchanged in each virtual segment, and the difference between two adjacent Mach is within the value of 0.05.

Uncertain Wind Scenario.
In this simulation study, the weather forecast data have been retrieved from the International Grand Global Ensemble (TIGGE) database of European Centre for Medium-Range Weather Forecasts (ECWMF), and the spatial resolution is 0.75°× 0.75°. e longitude and latitude ranges of the selected data in this research are 88°∼112°and 34°∼44°, respectively. e selected pressure altitude includes 200 mb, 250 mb, and 300 mb. e UTC time of this flight is from 05:30 to 08:30, and the forecast is released every 6h; hence, the latest forecast time which is 06 UTC forecast data released at 00:00 UTC is retrieved. en, the probabilistic forecast wind scenario is generated according to the following instructions: (1) e forecast data is spatial interpolated by the bilinear interpolation method. (2) According to the spatial correlations, the probabilistic forecast grids with the size of 135 nm × 135 nm are constructed; the historical prediction errors are analyzed for modeling wind uncertainties. e rootmean-square (RMS) error of the forecast between the historical forecast data and the observed data for nearly a month (i.e., 00 UTC data published at 06:00) in the statistical grid is calculated. (3) e forecast data for the day is used as the mean value of the probability distribution for each probabilistic forecast grid, and, ultimately, the Gaussian model is chosen as the distribution function. Finally, the probabilistic forecast grids are obtained.

Speed Profile Optimization
Results. First, the experiment settings for speed transition profile are as follows: a max � 1 m/s 2 , t max � 60 s, j min � 0.08 m/s 3 , and j max � 0.5 m/ s 3 . Following simple procedure for calculating the population size, N t � 10a, where a is the number of variables, suppose that j 1 � −j 2 , for these two phases are the inverse of each other. us, a � 5, and the population size is 50. e floating-point coding is adopted in generating the initial population. e crossover probability value is 0.8. e genetic algorithm for optimizing the entire cruise speed is set as follows: the population size N s � 80. During initialization, the Mach number range is limited to being from 0.68 to 0.82 according to the aircraft performance range.
e mutation probability p m � 0.02, the crossover probability p c � 0.85, the termination generation is 100, and the probability of satisfying the time window α � 0.95. e results of the optimized Mach speed profile are presented in Figure 8 and Table 2, aiming to demonstrate two groups of comparisons: (1) We run a simulation where the wind uncertainty is not accounted for, and only the deterministic forecast wind is used as the wind scenario in speed optimization, and the optimized speed profile is    Figure 8. e resulting RTA performances are compared in Table 2, and ΔTA 1 and ΔTA 2 are the arrival time error of the optimized trajectory without and with considering wind uncertainty in the optimization problem, respectively. Note that the absolute values of ΔTA 1 are greater than those of ΔTA 2 ; even ΔTA 1 of waypoints DOVOP, SUNUV, and NIXUK are greater than 30 s of the time window. Obviously, operating in this way, it cannot guarantee the accuracy of arrival time when great uncertainty occurs in weather prediction. (2) e profiles with and without considering speed transition are also plotted in Figure 8. Comparing both profiles, it can be seen that the variation of the curve with transition is reduced. It is also clear that the acceleration does exist, and the speed profile provides a smoother transition when changing speed. Meanwhile, the total fuel consumptions are calculated, while the fuel consumptions of the optimized trajectory without and with transitions are 11,990 lb and 12,022 lb, respectively. As expected, more fuel is consumed by the speed profile of the transitions being considered; it can be explained by the reason that additional fuel is consumed by the transition parts.
Here, we focus on analyzing the result in which uncertainty and speed transitions are both taken into account. e arrival time at each waypoint, the probability α of arriving at a waypoint within the time window, and the difference between ETA and RTA are listed in Table 2. It can be seen that, with introducing the probabilistic information of wind forecast, the absolute values of ΔTA 2 for all the waypoints are kept within 30 s, ranging from 2 s to 20 s. In addition, all values of probability α are not less than 95%, where P63 and NIXUK are the highest and lowest, respectively, among all the waypoints. e mean, standard deviation, and coefficient of variation of ETA for all segments are also provided in Table 3. e coefficient of variation is a normalized measure of the dispersion, indicating metric of a relative measure of uncertainty. It is easy to see that the spread of the flight time is different for different segments, with the standard deviation ranging from 2.1 s for segment 6 to 21.4 s for segment 7. A general trend can be seen: the higher the flight time, the higher the standard deviation of the flight time. e distance and flight time of segment 7 are the longest among those segments, and the standard deviations of ETA and fuel of segment 7 are much higher than those of the other segments. Besides, the waypoint NIXUK corresponding to the terminal waypoint of segment 7, whose probability α is shown in Table 2, is the lowest among all the waypoints. Hence, another trend can be seen: the longer the distance, the greater the randomness (i.e., the standard deviation and coefficient of variation) and the lower the probability of arrival within the time window. Also, note that the probability α of P63 in Table 2 is high with the value of 0.994; the above analysis can further account for this. Similar conclusions for the fuel consumption can be drawn from Table 4. Hence, the bigger the uncertainty in fuel consumption is, the more contingency fuel is needed to be loaded in dispatch when larger wind uncertainty occurs.
However, note that there is a similar distance for segments 4 and 5, which are 129.0 nm and 130.1 nm, respectively, whereas a bigger standard deviation is obtained for segment 4. Notice that the wind grid does not coincide with these two segments, and wind grids 1 and 2 are the winds experiencing on segments 4 and 5, respectively, as shown in Figure 9. Consequently, wind grid 1 occupying the main part experiencing on segment 4 shows larger uncertainty; through comparing its distribution, one obtains the result that the uncertainty in the flight time is larger in the case of wind grid 1, which performs wider spread than wind grid 2. Comparing  Figures 9(a) and 9(b), we see that the distributions of wind and time show a similar shape. is result indicates the accuracy of RTA will be greatly affected by wind uncertainty. In addition, the time error of arriving at BIKNO is still large; one possible explanation could be the error accumulation effect for prior segment 7, where a dramatic speed change happens in that segment to reduce the difference between ETA and RTA.
Eventually, the fuel efficiency of the optimized trajectory is tested. A total of 13,227 lb of fuel is burned by the original       trajectory of QAR data, while 12,022 lb is consumed by the optimized speed profile in the solid line in Figure 8, resulting in the reductions of 9.12% compared to the original flight. Consequently, the optimized speed profile has a significant effect on improving flight economy. A total of eight tests are performed to further verify the RTA performance and the profile optimization. e results are shown in Table 5. It can be observed that the absolute values of arrival time error for all the flights are within a time tolerance of 30 s, and the values of probability α for all the waypoints of these flights are above 95%. Notice that the minimum value of α of flight No. 3 is as low as 0.956, and this can be explained by checking the wind prediction accuracy and as a consequence that a large difference happened during that day. Besides, in general, the longer the range of the flight, the lower mean probability α, the result of which is consistent with that of the above typical flight. e ability of the fuel savings is tested as well, and an average fuel reduction of 7.04% can be observed when compared with the original trajectory of QAR data. e fuel reduction of flight No. 7 is up to 10.44%, which is higher than the average value for its longer range, and it indicates that long range flight will benefit more from trajectory optimization.

Conclusion
In this paper, we have developed a probability method for optimizing the speed profile subject to multiple RTA constraints, allowing accounting for wind uncertainty in trajectory planning. e impact of wind forecast uncertainty on 4D trajectory optimization is examined, and the probabilistic analysis presented is provided to explain the uncertainty propagation of the arrival time and fuel consumption and further to derive probability distributions by using the Probabilistic Transformation Method. Moreover, some improvements have been made by smoother speed guidance, and the speed transition profile is considered in the optimization process. e results shown in Section 4 prove that wind uncertainty will greatly influence the performance of the 4D trajectory, and thus taking wind uncertainty into account is crucial in fulfilling the RTA constraints, and therefore it is necessary to obtain the high-quality weather forecast. e speed profile optimized by this approach can comply with an RTA constraint and improving the fuel efficiency in the actual flight environment.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.