# Optimization of Speed Profile with RTA Constraints under Wind Uncertainty

Authors: Lisha Ye; Li Cao; Shuli Gong; Jiyun Lu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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.


---

## Body

## 1. Introduction

Air traffic increasing will bring about some new challenges to the airspace, such as increasing airspace tension and environmental impact. The 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. This 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. These 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 EUROCONTROL 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. Thus, 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) The 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. This will make it more difficult for the crew to fly 4D in the face of wind uncertainty. Then, 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. The 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. Therefore, 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]. The fulfillment of 4D trajectory in the future will much rely on the improvement of airborne capability. The 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. The 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. This 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.


## 2. Uncertainty Propagation Analysis

The 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. Then, the uncertain propagation characteristics of 4D trajectory are here detailed by analyzing how the wind affects the flight time and fuel consumption.


### 2.1. 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 trackW(P, t). The 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(1)WP,t=W^P,t+τt+ξP,t+τt,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 termW^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 quasi-homogeneous, 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. Thus, si 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:(2)Wisi,t+τit=W^isi,t+τit+ξisi,t+τit0.Eventually, the probabilistic wind model in the form of grid is discretized in time and space, as shown in Figure1. The 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,tt0∼N0,σ2. After superimposing the forecast value with the uncertain value, the mean value of wind speed Wi is the wind forecast value W^i, and finally the distribution of wind is expressed as Wisi,t∼NW^i,σi2.Figure 1
Probabilistic model for wind forecasting in the grid form.


### 2.2. 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. The cruise phase is given in multiple waypointsOi 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. The along-track wind componentVlw can be calculated from (3), and it takes a positive value for tailwind and a negative one for headwind:(3)Vlw=Vw⋅cosλ.As shown in Figure2, the actual ground speed Vg is calculated from the airspeed Va and the average along-track wind Vlw:(4)Vg=Va+Vlw.Vlw is wind speed along the track, Vw is wind speed, and λ is the angle of wind direction and head direction.Figure 2
Vector relationship of airspeed, ground speed, and wind.According to Section2.1, the forecast wind is modeled as a random variable, and the uncertain along-track wind is expressed as(5)Vlw=W^+ξ.The flight timeti of segment i is expressed as(6)ti=LiVgi=LiVai+W^i+ξi,where Li 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. Then, the uncertainty of flight time t can be derived. Therefore, due to time constraints, the airspeed Va 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. The fuel flow rates in cruiseffcr 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:(7)Fueltotal=∑iffcri×ti,As shown in (7), it can be seen that, due to the uncertainty of flight time, the total fuel consumption exhibits uncertainty as well.


## 2.1. 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 trackW(P, t). The 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(1)WP,t=W^P,t+τt+ξP,t+τt,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 termW^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 quasi-homogeneous, 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. Thus, si 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:(2)Wisi,t+τit=W^isi,t+τit+ξisi,t+τit0.Eventually, the probabilistic wind model in the form of grid is discretized in time and space, as shown in Figure1. The 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,tt0∼N0,σ2. After superimposing the forecast value with the uncertain value, the mean value of wind speed Wi is the wind forecast value W^i, and finally the distribution of wind is expressed as Wisi,t∼NW^i,σi2.Figure 1
Probabilistic model for wind forecasting in the grid form.


## 2.2. 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. The cruise phase is given in multiple waypointsOi 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. The along-track wind componentVlw can be calculated from (3), and it takes a positive value for tailwind and a negative one for headwind:(3)Vlw=Vw⋅cosλ.As shown in Figure2, the actual ground speed Vg is calculated from the airspeed Va and the average along-track wind Vlw:(4)Vg=Va+Vlw.Vlw is wind speed along the track, Vw is wind speed, and λ is the angle of wind direction and head direction.Figure 2
Vector relationship of airspeed, ground speed, and wind.According to Section2.1, the forecast wind is modeled as a random variable, and the uncertain along-track wind is expressed as(5)Vlw=W^+ξ.The flight timeti of segment i is expressed as(6)ti=LiVgi=LiVai+W^i+ξi,where Li 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. Then, the uncertainty of flight time t can be derived. Therefore, due to time constraints, the airspeed Va 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. The fuel flow rates in cruiseffcr 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:(7)Fueltotal=∑iffcri×ti,As shown in (7), it can be seen that, due to the uncertainty of flight time, the total fuel consumption exhibits uncertainty as well.


## 3. Speed Profile Optimization Method

In this section, the optimization of the speed profile considering the transition parts is proposed. The 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.


### 3.1. Speed Selection Subject to Uncertainty

#### 3.1.1. 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.The varying airspeed approach is adopted in this paper. For the convenience of airborne operation and separation management for ATM, virtual waypointsPk (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 VaK 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, Lmin ≥ 2l, and, in each virtual segment k, the average along-track wind Vlwk is calculated. Then the cruise route is divided into K subsegments by virtual waypoints and the RTA waypoints together. The arrival time tETAj of RTA waypoint j is as follows:(8)tETAj=∑k=1aj−1lVgk+Dj−aj−1lVgaj,where k is the index of virtual segment and RTA waypoint, j is the index of RTA waypoint, Dj is the distance of RTA waypoint j from TOC point, aj is the number of virtual segments contained in Dj, and aj = floor(Dj/l) + 1.Figure 3
Illustration of virtual waypoints and subsegments.For solving the stochastic optimization problem to minimize the fuel consumption, the following stochastic cost functionJ expressed in (9) is minimized:(9)J=Efueltotal=E∑k=1Kfuelk=E∑k=1Kffk⋅tk,where fuelk is the fuel consumed by virtual segment k, ffk is the fuel flow rate, tk 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 errortΔ, defined as the difference between ETA and RTA, should be limited within the predefined window T:(10)tΔ=tRTA−tETA≤T.Since time errortΔ depends additively on wind disturbance ξ, chance constraint is introduced to describe RTA constraints. The constraint may not be feasible, allowing the decision to violate the constraint to some extent before a random variable is observed. The arrival time of each waypoint has to be satisfied with a given probability αj of arriving in a time window:(11)PrtRTAj−tETAj≤Tj≥αj,j=1,2,…,M,where tRTAj and tETAj are time of RTA and ETA at waypoint j, respectively, and Tj is the time window.Mach number is adopted during the cruise, defined as Mach = cVa, where c is the speed of sound at the given flight pressure altitude, and its value should be kept within aircraft performance constraints:(12)Machmin≤Mach≤Machmax.


#### 3.1.2. Solution of the Optimization Problem

We will discuss how to solve the proposed model. There are two uncertain variables contained in the speed selection optimization problem: the arrival time and fuel consumption, both of which indicate 4D trajectory performance. The 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 pYy, 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:(13)pZz=pYhzJhhz.Figure 4
Probabilistic transformation method.For the problem studied in this paper, the relationship between flight time and wind speed in subsegmenti is as follows:(14)ti=ftiVwi=LiVai+Vwi.The inverse transformation is as follows:(15)htiti=f−1ti=Liti−Vai,(16)ΦVwi=∂ftiVwi∂VwiVwi=hVwiVwi=−LiVai+Vwi2Vw=hVwVw=−ti2Li,and, thus, the PDF of the flight time of subsegment can be derived by PTM:(17)ptiti=pVwihtitiΦVwi=Liti2⋅pVwiLiti−Vai.For independent random variablesX and Y, with the function relationship Z = X + Y, the PDF of Z can be obtained by the convolution formula of probability:(18)pZz=pXx∗pYy=∫−∞∞pXz−ypYydy=∫−∞∞pXxpYz−xdx.Assume that one segment is divided intoq subsegments, and the PDF of the flight time of the segment ptt can be computed from a set of convolutions:(19)ptt=pt1t1∗pt2t2∗⋯∗ptqtq.The 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(20)PrtRTA−tETA≤T=Prtearly≤t≤tlate=∫tearlytlatepttdt≥α,where tearly and tlate are the flight times corresponding to the allowable earliest and latest flight time, respectively.Figure 5
Probabilistic formulation of arriving at time window.Given a certain airspeed, the fuel flow is a determinate variable. Therefore, the PDF of fuel consumption in the subsegment is expressed as the product of fuel flow and the PDF of flight time:(21)pmfuelkmfuelk=ffk⋅ptktk.From the probability convolution (19), the PDF pmfuelmfuel of total fuel consumption is as follows:(22)pmfuelmfuel=pmfuel1mfuel1∗pmfuel2mfuel2∗⋯∗pmfuelKmfuelK.Eventually, given the wind PDF for each grid, the total fuel can be calculated, as shown in Figure6.Figure 6
Computation of the PDF of the total consumed fuel for the trajectory withm segments.The expectation and variance of total consumed fuel can be written as(23)Emfuel=∫0∞mfuelpmfuelmfueldmfuel,(24)Dmfuel=Emfuel2−Emfuel2.The airspeedVak (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. The 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. Therefore, 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 the exchange genes are generated randomly, and the new individual genes generated are located in the value range of the two parent genes.The 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 ofu=Mach1,Mach2,…,MachK randomly, and calculate U=Va1,Va2,…,VaKStep 2: computePrtRTA−tETA≤T by (8), (14)–(20); if constraint (20) is satisfied, then go to Step 3; if not, return to Step 1Step 3: use submatch method to calculateffkStep 4: calculatepmfuelkmfuelk using (21); then, PDF of total fuel consumption pmfuelmfuel is calculated by (22)Step 5: the expectation of consumed fuelEmfuel can be computed using (23)Step 6: search the optimalUopt=Va1,Va2,…,VaK corresponding to the minimum Emfuel by GA


### 3.2. Speed Transition Profile

#### 3.2.1. 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. Therefore, 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. The 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 zeroFor circumstance (ii), suppose the jerk in each phase is kept constant:j1, 0, and j2, respectively. Then, the acceleration can be deduced as follows:(25)at=a0+∫t0tfjtdt.Then, the acceleration of each phase is given by(26)at=j1t−t1,t∈t1,t2,j1t2−t1,t∈t2,t3,j1t2−t1+j2t−t3,t∈t3,t4.Furthermore, the velocity can be derived by integrating the acceleration equation:(27)Vt=V0+∫t0tfatdt.Combining (25)–(27), the detailed velocity can be expressed as(28)Vt=V1=V0+12j1t2−t12−j1t1t−t1,t∈t1,t2,V2=V1t2+j1t2−t1t−t2,t∈t2,t3,V3=V2t3+j1t2−j1t1−j2t3t−t3+12j2t2−t32,t∈t3,t4.From (26)–(28), it can be seen that the acceleration and velocity will be different with different values of j1, j2, t1, t2, t3, and t4, 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 t2 = t3, so the expressions of acceleration and velocity are omitted at present in this paper.The fuel burn calculation method described in the previous section indicates that the fuel flow rate is closely related to the acceleration and airspeed. Thus, 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.The optimization objective is to minimize the fuel consumption of two segments:(29)Obj=minJ=min∑i=t0tfffi×Δt.Considering the searching range, RTA requirements, aircraft performance, and passenger comfort, the inequality constraint sets are as follows:(30)0<t1<t2≤t3<t4<tf,(31)t1<tc,(32)t3>tc,(33)at1<amax,(34)t4−t1<tmax,(35)jmin<j<jmax,where tf is the RTA at the final point of the second segment:(36)tf=L1Vg0+L2Vgf,where tc is the spending time of flying the first segment at the initial speed, amax is the maximum acceleration of the aircraft, and tmax is the maximum acceleration time. jmin and jmax are the minimum and maximum jerk of aircraft, respectively.The constraints in (30)–(32) are regional constraints, which will reduce the searching region, restricting the aircraft to stay within the flight envelope. The constraint in (33) stipulates that the acceleration of the aircraft is limited within predefined values considering passenger comfort. The 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/s3 in this paper. Equality constraints are as follows:(37)Vtf=Vfinal,(38)at4=0,(39)tf=RTAsegment.The constraint in (37) shows that the final transition profile speed equals the target speed. The constraint in (38) indicates that the end acceleration of the transition profile equals zero.The 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:(40)gJ+rx⟶≡δ−hrx⟶≥0,r=1,…,R,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 hrx⟶ are rth equality constraints.


#### 3.2.2. Constraint Handling Optimization Algorithm

The constrained optimization problem in this paper is a nonlinear programming (NLP) problem. This 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. The 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 floating-point 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 single-point crossover used in binary GAs may generate a large portion of infeasible solutions, since such 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 functionFx⟶ shown in (41) is utilized to evaluate the optimization in terms of both objective function value and constraint violation degree:(41)Fx⟶=fx⟶,ifgjx⟶≥0,∀j=1,2,…,m,fmax+∑j=1mgjx⟶,otherwise,wherefx⟶ is objective function, fmax is the objective function value of the worst feasible solution in the population, gjx⟶ 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. The 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. The 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.


## 3.1. Speed Selection Subject to Uncertainty

### 3.1.1. 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.The varying airspeed approach is adopted in this paper. For the convenience of airborne operation and separation management for ATM, virtual waypointsPk (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 VaK 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, Lmin ≥ 2l, and, in each virtual segment k, the average along-track wind Vlwk is calculated. Then the cruise route is divided into K subsegments by virtual waypoints and the RTA waypoints together. The arrival time tETAj of RTA waypoint j is as follows:(8)tETAj=∑k=1aj−1lVgk+Dj−aj−1lVgaj,where k is the index of virtual segment and RTA waypoint, j is the index of RTA waypoint, Dj is the distance of RTA waypoint j from TOC point, aj is the number of virtual segments contained in Dj, and aj = floor(Dj/l) + 1.Figure 3
Illustration of virtual waypoints and subsegments.For solving the stochastic optimization problem to minimize the fuel consumption, the following stochastic cost functionJ expressed in (9) is minimized:(9)J=Efueltotal=E∑k=1Kfuelk=E∑k=1Kffk⋅tk,where fuelk is the fuel consumed by virtual segment k, ffk is the fuel flow rate, tk 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 errortΔ, defined as the difference between ETA and RTA, should be limited within the predefined window T:(10)tΔ=tRTA−tETA≤T.Since time errortΔ depends additively on wind disturbance ξ, chance constraint is introduced to describe RTA constraints. The constraint may not be feasible, allowing the decision to violate the constraint to some extent before a random variable is observed. The arrival time of each waypoint has to be satisfied with a given probability αj of arriving in a time window:(11)PrtRTAj−tETAj≤Tj≥αj,j=1,2,…,M,where tRTAj and tETAj are time of RTA and ETA at waypoint j, respectively, and Tj is the time window.Mach number is adopted during the cruise, defined as Mach = cVa, where c is the speed of sound at the given flight pressure altitude, and its value should be kept within aircraft performance constraints:(12)Machmin≤Mach≤Machmax.


### 3.1.2. Solution of the Optimization Problem

We will discuss how to solve the proposed model. There are two uncertain variables contained in the speed selection optimization problem: the arrival time and fuel consumption, both of which indicate 4D trajectory performance. The 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 pYy, 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:(13)pZz=pYhzJhhz.Figure 4
Probabilistic transformation method.For the problem studied in this paper, the relationship between flight time and wind speed in subsegmenti is as follows:(14)ti=ftiVwi=LiVai+Vwi.The inverse transformation is as follows:(15)htiti=f−1ti=Liti−Vai,(16)ΦVwi=∂ftiVwi∂VwiVwi=hVwiVwi=−LiVai+Vwi2Vw=hVwVw=−ti2Li,and, thus, the PDF of the flight time of subsegment can be derived by PTM:(17)ptiti=pVwihtitiΦVwi=Liti2⋅pVwiLiti−Vai.For independent random variablesX and Y, with the function relationship Z = X + Y, the PDF of Z can be obtained by the convolution formula of probability:(18)pZz=pXx∗pYy=∫−∞∞pXz−ypYydy=∫−∞∞pXxpYz−xdx.Assume that one segment is divided intoq subsegments, and the PDF of the flight time of the segment ptt can be computed from a set of convolutions:(19)ptt=pt1t1∗pt2t2∗⋯∗ptqtq.The 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(20)PrtRTA−tETA≤T=Prtearly≤t≤tlate=∫tearlytlatepttdt≥α,where tearly and tlate are the flight times corresponding to the allowable earliest and latest flight time, respectively.Figure 5
Probabilistic formulation of arriving at time window.Given a certain airspeed, the fuel flow is a determinate variable. Therefore, the PDF of fuel consumption in the subsegment is expressed as the product of fuel flow and the PDF of flight time:(21)pmfuelkmfuelk=ffk⋅ptktk.From the probability convolution (19), the PDF pmfuelmfuel of total fuel consumption is as follows:(22)pmfuelmfuel=pmfuel1mfuel1∗pmfuel2mfuel2∗⋯∗pmfuelKmfuelK.Eventually, given the wind PDF for each grid, the total fuel can be calculated, as shown in Figure6.Figure 6
Computation of the PDF of the total consumed fuel for the trajectory withm segments.The expectation and variance of total consumed fuel can be written as(23)Emfuel=∫0∞mfuelpmfuelmfueldmfuel,(24)Dmfuel=Emfuel2−Emfuel2.The airspeedVak (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. The 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. Therefore, 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 the exchange genes are generated randomly, and the new individual genes generated are located in the value range of the two parent genes.The 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 ofu=Mach1,Mach2,…,MachK randomly, and calculate U=Va1,Va2,…,VaKStep 2: computePrtRTA−tETA≤T by (8), (14)–(20); if constraint (20) is satisfied, then go to Step 3; if not, return to Step 1Step 3: use submatch method to calculateffkStep 4: calculatepmfuelkmfuelk using (21); then, PDF of total fuel consumption pmfuelmfuel is calculated by (22)Step 5: the expectation of consumed fuelEmfuel can be computed using (23)Step 6: search the optimalUopt=Va1,Va2,…,VaK corresponding to the minimum Emfuel by GA


## 3.1.1. 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.The varying airspeed approach is adopted in this paper. For the convenience of airborne operation and separation management for ATM, virtual waypointsPk (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 VaK 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, Lmin ≥ 2l, and, in each virtual segment k, the average along-track wind Vlwk is calculated. Then the cruise route is divided into K subsegments by virtual waypoints and the RTA waypoints together. The arrival time tETAj of RTA waypoint j is as follows:(8)tETAj=∑k=1aj−1lVgk+Dj−aj−1lVgaj,where k is the index of virtual segment and RTA waypoint, j is the index of RTA waypoint, Dj is the distance of RTA waypoint j from TOC point, aj is the number of virtual segments contained in Dj, and aj = floor(Dj/l) + 1.Figure 3
Illustration of virtual waypoints and subsegments.For solving the stochastic optimization problem to minimize the fuel consumption, the following stochastic cost functionJ expressed in (9) is minimized:(9)J=Efueltotal=E∑k=1Kfuelk=E∑k=1Kffk⋅tk,where fuelk is the fuel consumed by virtual segment k, ffk is the fuel flow rate, tk 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 errortΔ, defined as the difference between ETA and RTA, should be limited within the predefined window T:(10)tΔ=tRTA−tETA≤T.Since time errortΔ depends additively on wind disturbance ξ, chance constraint is introduced to describe RTA constraints. The constraint may not be feasible, allowing the decision to violate the constraint to some extent before a random variable is observed. The arrival time of each waypoint has to be satisfied with a given probability αj of arriving in a time window:(11)PrtRTAj−tETAj≤Tj≥αj,j=1,2,…,M,where tRTAj and tETAj are time of RTA and ETA at waypoint j, respectively, and Tj is the time window.Mach number is adopted during the cruise, defined as Mach = cVa, where c is the speed of sound at the given flight pressure altitude, and its value should be kept within aircraft performance constraints:(12)Machmin≤Mach≤Machmax.


## 3.1.2. Solution of the Optimization Problem

We will discuss how to solve the proposed model. There are two uncertain variables contained in the speed selection optimization problem: the arrival time and fuel consumption, both of which indicate 4D trajectory performance. The 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 pYy, 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:(13)pZz=pYhzJhhz.Figure 4
Probabilistic transformation method.For the problem studied in this paper, the relationship between flight time and wind speed in subsegmenti is as follows:(14)ti=ftiVwi=LiVai+Vwi.The inverse transformation is as follows:(15)htiti=f−1ti=Liti−Vai,(16)ΦVwi=∂ftiVwi∂VwiVwi=hVwiVwi=−LiVai+Vwi2Vw=hVwVw=−ti2Li,and, thus, the PDF of the flight time of subsegment can be derived by PTM:(17)ptiti=pVwihtitiΦVwi=Liti2⋅pVwiLiti−Vai.For independent random variablesX and Y, with the function relationship Z = X + Y, the PDF of Z can be obtained by the convolution formula of probability:(18)pZz=pXx∗pYy=∫−∞∞pXz−ypYydy=∫−∞∞pXxpYz−xdx.Assume that one segment is divided intoq subsegments, and the PDF of the flight time of the segment ptt can be computed from a set of convolutions:(19)ptt=pt1t1∗pt2t2∗⋯∗ptqtq.The 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(20)PrtRTA−tETA≤T=Prtearly≤t≤tlate=∫tearlytlatepttdt≥α,where tearly and tlate are the flight times corresponding to the allowable earliest and latest flight time, respectively.Figure 5
Probabilistic formulation of arriving at time window.Given a certain airspeed, the fuel flow is a determinate variable. Therefore, the PDF of fuel consumption in the subsegment is expressed as the product of fuel flow and the PDF of flight time:(21)pmfuelkmfuelk=ffk⋅ptktk.From the probability convolution (19), the PDF pmfuelmfuel of total fuel consumption is as follows:(22)pmfuelmfuel=pmfuel1mfuel1∗pmfuel2mfuel2∗⋯∗pmfuelKmfuelK.Eventually, given the wind PDF for each grid, the total fuel can be calculated, as shown in Figure6.Figure 6
Computation of the PDF of the total consumed fuel for the trajectory withm segments.The expectation and variance of total consumed fuel can be written as(23)Emfuel=∫0∞mfuelpmfuelmfueldmfuel,(24)Dmfuel=Emfuel2−Emfuel2.The airspeedVak (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. The 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. Therefore, 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 the exchange genes are generated randomly, and the new individual genes generated are located in the value range of the two parent genes.The 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 ofu=Mach1,Mach2,…,MachK randomly, and calculate U=Va1,Va2,…,VaKStep 2: computePrtRTA−tETA≤T by (8), (14)–(20); if constraint (20) is satisfied, then go to Step 3; if not, return to Step 1Step 3: use submatch method to calculateffkStep 4: calculatepmfuelkmfuelk using (21); then, PDF of total fuel consumption pmfuelmfuel is calculated by (22)Step 5: the expectation of consumed fuelEmfuel can be computed using (23)Step 6: search the optimalUopt=Va1,Va2,…,VaK corresponding to the minimum Emfuel by GA


## 3.2. Speed Transition Profile

### 3.2.1. 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. Therefore, 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. The 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 zeroFor circumstance (ii), suppose the jerk in each phase is kept constant:j1, 0, and j2, respectively. Then, the acceleration can be deduced as follows:(25)at=a0+∫t0tfjtdt.Then, the acceleration of each phase is given by(26)at=j1t−t1,t∈t1,t2,j1t2−t1,t∈t2,t3,j1t2−t1+j2t−t3,t∈t3,t4.Furthermore, the velocity can be derived by integrating the acceleration equation:(27)Vt=V0+∫t0tfatdt.Combining (25)–(27), the detailed velocity can be expressed as(28)Vt=V1=V0+12j1t2−t12−j1t1t−t1,t∈t1,t2,V2=V1t2+j1t2−t1t−t2,t∈t2,t3,V3=V2t3+j1t2−j1t1−j2t3t−t3+12j2t2−t32,t∈t3,t4.From (26)–(28), it can be seen that the acceleration and velocity will be different with different values of j1, j2, t1, t2, t3, and t4, 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 t2 = t3, so the expressions of acceleration and velocity are omitted at present in this paper.The fuel burn calculation method described in the previous section indicates that the fuel flow rate is closely related to the acceleration and airspeed. Thus, 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.The optimization objective is to minimize the fuel consumption of two segments:(29)Obj=minJ=min∑i=t0tfffi×Δt.Considering the searching range, RTA requirements, aircraft performance, and passenger comfort, the inequality constraint sets are as follows:(30)0<t1<t2≤t3<t4<tf,(31)t1<tc,(32)t3>tc,(33)at1<amax,(34)t4−t1<tmax,(35)jmin<j<jmax,where tf is the RTA at the final point of the second segment:(36)tf=L1Vg0+L2Vgf,where tc is the spending time of flying the first segment at the initial speed, amax is the maximum acceleration of the aircraft, and tmax is the maximum acceleration time. jmin and jmax are the minimum and maximum jerk of aircraft, respectively.The constraints in (30)–(32) are regional constraints, which will reduce the searching region, restricting the aircraft to stay within the flight envelope. The constraint in (33) stipulates that the acceleration of the aircraft is limited within predefined values considering passenger comfort. The 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/s3 in this paper. Equality constraints are as follows:(37)Vtf=Vfinal,(38)at4=0,(39)tf=RTAsegment.The constraint in (37) shows that the final transition profile speed equals the target speed. The constraint in (38) indicates that the end acceleration of the transition profile equals zero.The 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:(40)gJ+rx⟶≡δ−hrx⟶≥0,r=1,…,R,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 hrx⟶ are rth equality constraints.


### 3.2.2. Constraint Handling Optimization Algorithm

The constrained optimization problem in this paper is a nonlinear programming (NLP) problem. This 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. The 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 floating-point 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 single-point crossover used in binary GAs may generate a large portion of infeasible solutions, since such 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 functionFx⟶ shown in (41) is utilized to evaluate the optimization in terms of both objective function value and constraint violation degree:(41)Fx⟶=fx⟶,ifgjx⟶≥0,∀j=1,2,…,m,fmax+∑j=1mgjx⟶,otherwise,wherefx⟶ is objective function, fmax is the objective function value of the worst feasible solution in the population, gjx⟶ 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. The 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. The 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.


## 3.2.1. 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. Therefore, 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. The 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 zeroFor circumstance (ii), suppose the jerk in each phase is kept constant:j1, 0, and j2, respectively. Then, the acceleration can be deduced as follows:(25)at=a0+∫t0tfjtdt.Then, the acceleration of each phase is given by(26)at=j1t−t1,t∈t1,t2,j1t2−t1,t∈t2,t3,j1t2−t1+j2t−t3,t∈t3,t4.Furthermore, the velocity can be derived by integrating the acceleration equation:(27)Vt=V0+∫t0tfatdt.Combining (25)–(27), the detailed velocity can be expressed as(28)Vt=V1=V0+12j1t2−t12−j1t1t−t1,t∈t1,t2,V2=V1t2+j1t2−t1t−t2,t∈t2,t3,V3=V2t3+j1t2−j1t1−j2t3t−t3+12j2t2−t32,t∈t3,t4.From (26)–(28), it can be seen that the acceleration and velocity will be different with different values of j1, j2, t1, t2, t3, and t4, 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 t2 = t3, so the expressions of acceleration and velocity are omitted at present in this paper.The fuel burn calculation method described in the previous section indicates that the fuel flow rate is closely related to the acceleration and airspeed. Thus, 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.The optimization objective is to minimize the fuel consumption of two segments:(29)Obj=minJ=min∑i=t0tfffi×Δt.Considering the searching range, RTA requirements, aircraft performance, and passenger comfort, the inequality constraint sets are as follows:(30)0<t1<t2≤t3<t4<tf,(31)t1<tc,(32)t3>tc,(33)at1<amax,(34)t4−t1<tmax,(35)jmin<j<jmax,where tf is the RTA at the final point of the second segment:(36)tf=L1Vg0+L2Vgf,where tc is the spending time of flying the first segment at the initial speed, amax is the maximum acceleration of the aircraft, and tmax is the maximum acceleration time. jmin and jmax are the minimum and maximum jerk of aircraft, respectively.The constraints in (30)–(32) are regional constraints, which will reduce the searching region, restricting the aircraft to stay within the flight envelope. The constraint in (33) stipulates that the acceleration of the aircraft is limited within predefined values considering passenger comfort. The 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/s3 in this paper. Equality constraints are as follows:(37)Vtf=Vfinal,(38)at4=0,(39)tf=RTAsegment.The constraint in (37) shows that the final transition profile speed equals the target speed. The constraint in (38) indicates that the end acceleration of the transition profile equals zero.The 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:(40)gJ+rx⟶≡δ−hrx⟶≥0,r=1,…,R,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 hrx⟶ are rth equality constraints.


## 3.2.2. Constraint Handling Optimization Algorithm

The constrained optimization problem in this paper is a nonlinear programming (NLP) problem. This 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. The 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 floating-point 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 single-point crossover used in binary GAs may generate a large portion of infeasible solutions, since such 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 functionFx⟶ shown in (41) is utilized to evaluate the optimization in terms of both objective function value and constraint violation degree:(41)Fx⟶=fx⟶,ifgjx⟶≥0,∀j=1,2,…,m,fmax+∑j=1mgjx⟶,otherwise,wherefx⟶ is objective function, fmax is the objective function value of the worst feasible solution in the population, gjx⟶ 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. The 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. The 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.


## 4. Numerical Results

In this section, an experiment of a typical domestic flight in China has been carried out. The 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. The cruising altitude is 25,600 ft, the initial speed is 0.745 M, and the initial TOC weight is 151,60 lb. The actual horizontal trajectory drawn from the QAR data during the cruise is shown in Figure7(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.Figure 7
Horizontal trajectory of case flight. (a) Horizontal trajectory. (b) Waypoints.
(a)
(b)Table 1
Waypoints and RTA.
Waypoint
Longitude
Latitude
RTA
TOC
112.3933
34.1050
13:28
P63
110.5850
33.9367
13:41
DOVOP
108.7083
34.7567
13:56
P03
106.6933
35.9817
14:13
SUNUV
104.1133
36.5300
14:32
LAXEV
101.6850
37.5050
14:51
NOGEX
100.9183
37.8300
15:02
NIXUK
95.1450
40.2000
15:41
BIKNO
93.8483
40.9150
15:51
TOD
92.6020
41.8576
16:01The 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. Therefore, 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.(3)
The Mach number should be kept within the performance envelop of Boeing 737–800, Machmin = 0.68 and Machmax = 0.82.


### 4.1. 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°. The longitude and latitude ranges of the selected data in this research are 88°∼112° and 34°∼44°, respectively. The selected pressure altitude includes 200 mb, 250 mb, and 300 mb. The 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. Then, the probabilistic forecast wind scenario is generated according to the following instructions:(1)
The 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. The root-mean-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)
The 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.


### 4.2. Speed Profile Optimization Results

First, the experiment settings for speed transition profile are as follows:amax = 1 m/s2, tmax = 60 s, jmin = 0.08 m/s3, and jmax = 0.5 m/s3. Following simple procedure for calculating the population size, Nt = 10a, where a is the number of variables, suppose that j1 = −j2, for these two phases are the inverse of each other. Thus, a = 5, and the population size is 50. The floating-point coding is adopted in generating the initial population. The crossover probability value is 0.8.The genetic algorithm for optimizing the entire cruise speed is set as follows: the population sizeNs = 80. During initialization, the Mach number range is limited to being from 0.68 to 0.82 according to the aircraft performance range. The mutation probability pm = 0.02, the crossover probability pc = 0.85, the termination generation is 100, and the probability of satisfying the time window α = 0.95.The results of the optimized Mach speed profile are presented in Figure8 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 plotted in Figure8. The resulting RTA performances are compared in Table 2, and ΔTA1 and ΔTA2 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 ΔTA1 are greater than those of ΔTA2; even ΔTA1 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)
The profiles with and without considering speed transition are also plotted in Figure8. 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.Figure 8
Optimized speed profile.Table 2
Comparison of RTA and ETA.
Waypoint
RTA
ETA1
ΔTA1
ETA2
ΔTA2
α
TOC
13:28:00
13:28:00
0 s
13:28:00
0 s
—
P63
13:41:00
13:41:04
4 s
13:41:02
2 s
0.994
DOVOP
13:56:00
13:55:28
−32 s
13:55:57
−3 s
0.983
P03
14:13:00
14:13:20
20 s
14:13:10
10 s
0.975
SUNUV
14:32:00
14:32:37
37 s
14:31:56
−4 s
0.986
LAXEV
14:51:00
14:50:50
−10 s
14:51:03
3 s
0.994
NOGEX
15:02:00
15:02:10
10 s
15:02:09
9 s
0.995
NIXUK
15:41:00
15:40:29
−31 s
15:40:40
−20 s
0.966
BIKNO
15:51:00
15:50:41
−19 s
15:50:43
−17 s
0.971
TOD
16:01:00
16:01:12
12 s
16:01:06
6 s
0.971
Note: the negative value indicates early arrival, and positive value indicates late arrival.Here, we focus on analyzing the result in which uncertainty and speed transitions are both taken into account. The 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 ΔTA2 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.The mean, standard deviation, and coefficient of variation of ETA for all segments are also provided in Table3. The 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. The 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.Table 3
Mean, standard deviation, and coefficient of variation of ETA.
Segment
Range (nm)
E[t] (s)
S[t] (s)
S
t
/
E
t·10−3
1
80.5
787.2
4.7
5.97
2
105.3
899.1
7.4
8.23
3
123.1
1005.5
8.9
8.85
4
129.0
1134.0
8.5
7.49
5
130.1
1149.8
5.9
5.13
6
41.3
671.2
2.1
3.13
7
304.5
2322.7
21.4
9.21
8
72.9
609.3
6.4
10.50
9
78.7
1625.3
14.4
8.92Table 4
Mean, standard deviation, and coefficient of variation of fuel.
Segment
E[fuel] (lb)
S[fuel] (lb)
S
fuel
/
E
fuel·10−3
1
996.4
4.6
4.76
2
1164.1
7.2
6.18
3
1348.9
8.7
6.44
4
1418.2
8.6
6.06
5
1420.5
7.0
4.93
6
449.4
2.1
4.67
7
3358.4
25.5
7.59
8
804.6
6.1
7.58
9
906.0
6.3
6.26However, 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 Figure9. 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. This 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.Figure 9
PDFs of wind speed and flight time. (a) Wind speed PDF. (b) Flight time PDF.
(a)
(b)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 Figure8, 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. The results are shown in Table5. 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.Table 5
Fuel reduction and RTA tests.
No.
Flight
Range (nm)
Abs (ΔTAmax) (s)
Mean Abs (ΔTA) (s)
α
min
Meanα
Fuel saving (%)
1
PEK-SHA
454
18
6
0.968
0.977
6.21
2
SHA-PEK
454
12
5
0.970
0.985
6.34
3
CAN-PEK
686
20
11
0.956
0.973
6.56
4
WUH-CTU
815
18
9
0.972
0.990
7.03
5
WUH-SHA
388
7
4
0.979
0.992
5.88
6
WUH-PEK
468
11
5
0.973
0.979
5.79
7
URC-HFE
1602
19
10
0.968
0.980
10.44
8
SZX-PEK
617
15
5
0.977
0.986
8.12The 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. The 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.


## 4.1. 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°. The longitude and latitude ranges of the selected data in this research are 88°∼112° and 34°∼44°, respectively. The selected pressure altitude includes 200 mb, 250 mb, and 300 mb. The 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. Then, the probabilistic forecast wind scenario is generated according to the following instructions:(1)
The 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. The root-mean-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)
The 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.


## 4.2. Speed Profile Optimization Results

First, the experiment settings for speed transition profile are as follows:amax = 1 m/s2, tmax = 60 s, jmin = 0.08 m/s3, and jmax = 0.5 m/s3. Following simple procedure for calculating the population size, Nt = 10a, where a is the number of variables, suppose that j1 = −j2, for these two phases are the inverse of each other. Thus, a = 5, and the population size is 50. The floating-point coding is adopted in generating the initial population. The crossover probability value is 0.8.The genetic algorithm for optimizing the entire cruise speed is set as follows: the population sizeNs = 80. During initialization, the Mach number range is limited to being from 0.68 to 0.82 according to the aircraft performance range. The mutation probability pm = 0.02, the crossover probability pc = 0.85, the termination generation is 100, and the probability of satisfying the time window α = 0.95.The results of the optimized Mach speed profile are presented in Figure8 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 plotted in Figure8. The resulting RTA performances are compared in Table 2, and ΔTA1 and ΔTA2 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 ΔTA1 are greater than those of ΔTA2; even ΔTA1 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)
The profiles with and without considering speed transition are also plotted in Figure8. 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.Figure 8
Optimized speed profile.Table 2
Comparison of RTA and ETA.
Waypoint
RTA
ETA1
ΔTA1
ETA2
ΔTA2
α
TOC
13:28:00
13:28:00
0 s
13:28:00
0 s
—
P63
13:41:00
13:41:04
4 s
13:41:02
2 s
0.994
DOVOP
13:56:00
13:55:28
−32 s
13:55:57
−3 s
0.983
P03
14:13:00
14:13:20
20 s
14:13:10
10 s
0.975
SUNUV
14:32:00
14:32:37
37 s
14:31:56
−4 s
0.986
LAXEV
14:51:00
14:50:50
−10 s
14:51:03
3 s
0.994
NOGEX
15:02:00
15:02:10
10 s
15:02:09
9 s
0.995
NIXUK
15:41:00
15:40:29
−31 s
15:40:40
−20 s
0.966
BIKNO
15:51:00
15:50:41
−19 s
15:50:43
−17 s
0.971
TOD
16:01:00
16:01:12
12 s
16:01:06
6 s
0.971
Note: the negative value indicates early arrival, and positive value indicates late arrival.Here, we focus on analyzing the result in which uncertainty and speed transitions are both taken into account. The 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 ΔTA2 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.The mean, standard deviation, and coefficient of variation of ETA for all segments are also provided in Table3. The 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. The 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.Table 3
Mean, standard deviation, and coefficient of variation of ETA.
Segment
Range (nm)
E[t] (s)
S[t] (s)
S
t
/
E
t·10−3
1
80.5
787.2
4.7
5.97
2
105.3
899.1
7.4
8.23
3
123.1
1005.5
8.9
8.85
4
129.0
1134.0
8.5
7.49
5
130.1
1149.8
5.9
5.13
6
41.3
671.2
2.1
3.13
7
304.5
2322.7
21.4
9.21
8
72.9
609.3
6.4
10.50
9
78.7
1625.3
14.4
8.92Table 4
Mean, standard deviation, and coefficient of variation of fuel.
Segment
E[fuel] (lb)
S[fuel] (lb)
S
fuel
/
E
fuel·10−3
1
996.4
4.6
4.76
2
1164.1
7.2
6.18
3
1348.9
8.7
6.44
4
1418.2
8.6
6.06
5
1420.5
7.0
4.93
6
449.4
2.1
4.67
7
3358.4
25.5
7.59
8
804.6
6.1
7.58
9
906.0
6.3
6.26However, 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 Figure9. 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. This 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.Figure 9
PDFs of wind speed and flight time. (a) Wind speed PDF. (b) Flight time PDF.
(a)
(b)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 Figure8, 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. The results are shown in Table5. 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.Table 5
Fuel reduction and RTA tests.
No.
Flight
Range (nm)
Abs (ΔTAmax) (s)
Mean Abs (ΔTA) (s)
α
min
Meanα
Fuel saving (%)
1
PEK-SHA
454
18
6
0.968
0.977
6.21
2
SHA-PEK
454
12
5
0.970
0.985
6.34
3
CAN-PEK
686
20
11
0.956
0.973
6.56
4
WUH-CTU
815
18
9
0.972
0.990
7.03
5
WUH-SHA
388
7
4
0.979
0.992
5.88
6
WUH-PEK
468
11
5
0.973
0.979
5.79
7
URC-HFE
1602
19
10
0.968
0.980
10.44
8
SZX-PEK
617
15
5
0.977
0.986
8.12The 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. The 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.


## 5. 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. The 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.The results shown in Section4 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. The speed profile optimized by this approach can comply with an RTA constraint and improving the fuel efficiency in the actual flight environment.



---

*Source: 4803747-2020-05-13.xml*