Algorithms for aircraft track optimization in flexible routing

. The authors consider the problem of optimization of aircraft flight tracks in air traffic management (ATM) on basis of flexible routing technologies which involve the use of satellite navigation systems (SNS). It is shown that in optimizing a trajectory it is necessary to take into account the accuracy of track holding in flight which depends on accuracy of the navigation system and external flight path disturbances, e.g. wind. For solving the problem of optimization the authors propose to use the theory of graphs. The technique of constructing a dynamic SNS accuracy field and representing it as a graph was developed. It is proposed that the SNS field could be characterized by geometric dilution of precision changing both in space and in time. Based on the theory of graphs (A-star algorithm) the technique of constructing a trajectory of optimal length under conditions of SNS accuracy variations and external flight path disturbances is proposed. The criterion of optimization based on minimizing the true track is offered. The cost function taking into account the track holding accuracy in navigating by SNS and effects of external flight disturbances is justified. The article represents the results of A-star algorithm application for optimal flight track construction under conditions of SNS accuracy field variation and presence of prohibited zones in the provide zone of airspace


Introduction
Evolution of the world civil aviation is accompanied by the development and introduction of new technologies directed to enhancement of the ATM system. From a contemporary perspective the ATM system should provide the effective use of airspace and its high capacity under conditions of increased density and intensity of air traffic when keeping or even increasing the flight safety level. This property of the ATM system is based on optimization of airspace structure and use providing for, among other things, optimization of aircraft flight tracks. Optimization directed to reducing the length of the trajectories may lead to economy of flight time and fuel consumption as well as to decrease of adverse impact of airсraft on the environment. Navigation technologies that increase the efficiency of airspace use include the already well-established Area Navigation (RNAV) and upcoming Flexible Routing and Trajectory Based Operations (TBO) [1]. These technologies are based on the use of GPS and GLONASS satellite navigation systems (SNS) as the main means of high-precision navigation sightings of an aircraft.
The Flexible Routing involves the provision of the appropriate level of crew's and controllers' situational awareness about the state of the provided airspace and ATM in it including the recommended flight trajectory. As a rule, the level of situational awareness nowadays is determined by the information about the weather conditions and prohibited zones in the airspace as well as about the accuracy of determining the coordinates and aircraft motion parameters provided by onboard navigation means (e.g. RAIM or AAIM function in users' SNS equipment; onboard integrated navigation systems meeting RNAV and RNP specifications) [2,3]. These factors are key ones for the crew to take a decision about a new, optimal for particular conditions, trajectory for the further flight providing the required level of flight safety.
In the existing approach, however, the new flight track chosen by the crew or recommended by the controller doesn't take into account probable changes of provided airspace characteristics which can have an effect on the accuracy of track holding. Thus, accuracy of positioning by SNS data in different zones of the airspace and at different points in time depends on the changing geometric dilution of precision of the system. The external flight path disturbances related to the wind and atmospheric turbulence can also change in space and time. So the trajectory chosen when entering the provided airspace can become a non-optimal one.
Hence, the solution of a track optimization problem should take into account the change of aeronautical (in terms of aircraft positioning accuracy provided by the board means), air and weather conditions during the flight. So, the construction of an optimal trajectory should be based on the forecast of track holding accuracy in the airspace with changing characteristics. This causes a problem of creating an optimality criterion allowing the forecasted accuracy of track holding during the flight to be taken into account.
In some practically important tasks the minimal distance between the initial and final points (the minimal length of the desired track) is used as a criterion of optimality for the choice of a flight trajectory. The chosen trajectory should also meet the conditions providing the required level of flight safety (avoid flying though prohibited zones, zones of dangerous moisture targets as well as near-midair collisions).
It should be noted that modern ATM technologies, e.g. those of RNAV, require desired track holding with specified accuracy which is determined by D TSE dispersion of the Total System Error (TSE) [3] where D NSE is dispersion of the Navigation System Error, D FTE is dispersion of the Flight Technical Error, D PDE is dispersion of the Path Definition Error. The optimal defined track is a desired track (DT) characterized by length L DT . Due to inaccuracy of DT holding the aircraft will fly along the true track (TT) which is characterized by length L TT . L TT = L DT +∆L TSE , where ∆L TSE is an increment of the DT length because of inaccuracy of DT holding. Let's consider two possible defined tracks ( fig.1): an optimal one, characterized by length L DT1 , and a non-optimal one, characterized by length L DT2 > L DT1 . Under particular conditions depending on the accuracy of track holding (with L TSE1 >L TSE2 ), a situation can arise when L DT2 +∆L TSE2 < L DT1 +∆L TSE1 . Then inequality L TT2 < L TT1 is correct, i.e. the length of a true track when flying along a non-optimal (in terms of length) defined track is shorter than that when flying along an optimal defined track. It is therefore important to take into account the accuracy of track holding during the further flight when choosing a new track in a flexible routing problem.
The problem of flexible routing can be solved in two steps. The first step is to construct an optimal track in the given zone of airspace. The second step involves flying along an optimal track with the fewest deviations from it.
For construction of an optimal flight track the methods and algorithms used in field of knowledge related to artificial intelligence can be applied [4]. Thus, work [5] considers the problem of constructing a flexible 4-D track of aircraft landing with use of genetic algorithms. In works [6][7] for track optimization the methods and algorithms of graph theory are used. In particular, work [7] studies the algorithm synthesized with account of weather conditions which allows the flight time and fuel consumption to be reduced.
For solving the problem of aircraft flight along an optimal track the methods of optimal control theory are widely used, different optimization criteria being applied. Indeed, in [8] the issue of track optimization with account of wind fields is researched, in [9] the problem of track multicriterion optimization is considered, in [10] the time-referenced optimization of track (4D-track) is proposed.
In considered works, however, the aircraft flight track is optimized without taking into account the accuracy of track keeping under changing aeronautical and air conditions in the provided airspace.
The modern and promising navigation technologies are based on the SNS use as the most accurate means of positioning. We will, therefore, consider the SNS receiver as a navigational device defining the value of dispersion of Navigation System Error D NSE in expression (1).
The accuracy of SNS positioning depends on location of navigation satellites relative to an aircraft and is characterized by dilution of precision which changes both in time and in space [2,11]. Correspondingly, the accuracy of SNS positioning also changes in time and space, thus leading to changes of track holding accuracy. Therefore, an optimal (in terms of length) true track of aircraft flight should be constructed with account of forecasted track holding in the SNS accuracy field changing in space and time and characterized by DOP values.
The article considers applying the graph theory А-star algorithm for solving the problem of constructing optimal tracks. The difference from the known solutions is in the use of a new optimization criterion. As such a criterion we have chosen the minimum length of the true track taking into account the forecast of track holding accuracy. For this purpose, the graph vertex weights are proposed to be formed on the basis of forecast of the SNS accuracy field characteristics in the given zone of airspace and external flight path disturbances. Besides, the article considers using the algorithm of track optimization in the airspace having prohibited zones.

Research methodology
With regard to the air transport system, the airspace or its separate elements among which the aircraft move can be represented as a route network. Such networks are modelled as graphs and a set of specifically connected vertices allows the optimal flight route to be constructed between the entry and exit points of the airspace zone [12]. Dijkstra's and Astar algorithms are the most known graph theory algorithms for searching the shortest path between two points [6].
The Dijkstra's algorithm finds the shortest track from a desired graph vertex to all the others. Each vertex is assigned by a weight characterizing the distance to it from an adjacent vertex (edge weight). The track from the initial graph vertex to the final one with minimal sum of edge weights is optimal (F(L DT ) is the cost function, L DT is the length of the route (DT)).
A-star algorithm uses the strategy of informed search combining the mathematical and heuristic approaches. The heuristic approach involves the use of knowledge related to the given problematic field (task) allowing the A-star algorithm to be applied in artificial intelligent systems thus providing significant reduction of the computational costs. In the problem of constructing an optimal DT the heuristic approach is based on the fact that in navigation the shortest distance line between two points of a route is a section of a great circle.
The A-star algorithm minimizes the cost function of type F(L DT )=G 0i,j +H i,j . G 0i,j is the cost function of reaching the desired vertex (i,j) from the initial graph vertex (the aircraft entry point in the given airspace zone) and H i,j is a heuristic estimate of the distance from the desired vertex (i, j) to the final graph vertex (the aircraft exit point). In the problem to be solved the heuristic estimate of distance is the length of the section of a great circle.
The authors suggest taking into account the random error of track holding TSE in the cost function and representing the function as where G 0i,j is the DT length from the initial graph vertex to the vertex with number (i,j); K i,j (TSE) is a coefficient taking into account DT length extension due to the track holding random error; H i,j is a heuristic estimate of the distance from the desired vertex to the final graph vertex (great circle length). The optimal track of aircraft flight should meet condition where L opt is the TT length along the optimal constructed track; e(G) is the number of graph edges used for track construction from the initial to the final graph vertex; F TTi,j is the weight of graph vertex with number (i,j). If NSE, FTE and PDE errors are mutually independent, the equation The value of K i,j (PDE) coefficient depends on accuracy of DT definition in the flight management system. Further, we will assume that the error is much smaller than the other TSE constituents and exclude it from consideration.
The value of K i,j (NSE) coefficient is defined by accuracy of the onboard navigation system. When using SNS the accuracy of aircraft positioning is determined by accuracy of the SNS onboard receiver. It is also the case if it is used for correcting the coordinates defined by the onboard dead reckoning system, e.g. inertial navigation system. The SNS receiver accuracy is defined by the error of measuring the pseudoranges to the satellites and the value of positional dilution of precision PDOP in the observation point [2].
The value of K i,j (FTE) coefficient depends on the aerodynamic performance of aircraft and its weight as well as on the air route conditions (wind, atmospheric turbulence etc).
Let's assume that DT extension due to the influence of random NSE can be approximated by the process of type where α NSE is a time constant of error correlation; σ NSE is a root-mean-square error (RMS); ω NSE is forming white Gaussian noise with zero mathematical expectation and unit intensity.
The RMS error of DT extension in navigating by SNS can be represented as σ NSE =PDOP⸱σ R where σ R is the RMS of measuring the pseudorange to the satellites.
By data of online monitoring of May 16, 2019 from the site of Russian system of differential correction and monitoring (http://www.sdcm.ru), the maximal RMS of measuring the pseudorange to the GLONASS satellites was 15.74 meters in the observation point. Given probable abnormal conditions of receiver operation and its installation on a highly dynamic object we will take σ R =50 meters.
As a typical value let's take α NSE =0.01 Hz because when smoothing the measurements the SNS receiver output error is a fairly narrow-band (slowly changing) process.
To find K i,j (NSE) values experimentally the statistical test method was used. As a result it was stated that with the chosen initial data the K i,j (NSE) values are within 0.001-0.025 while PDOP changes from 1 to 6.5. While α NSE =0.01 Hz and PDOP changes from 1 to 6.5 the K i,j (NSE) values are within 0.03-0.17. Hence, widening of error fluctuation spectrum at the SNS receiver output leads to increase of TT extension thus proving validity of the used model.
To find K i,j (NSE) coefficient values one can use the known model of controlled movement of aircraft under external flight path disturbances [13] where ∆L FTE (t) is DT extension due to flight path disturbances; ∆W(t) and α(t) are vector projections of fluctuation of ground speed and acceleration of aircraft when moving along the DT; δ and β are coefficients characterizing the spectral density of random variations of acceleration due to external disturbances, type of object and conditions of its movement; σ 2 α is acceleration fluctuation dispersion; n α (t) is forming white Gaussian noise with zero mathematical expectation and unit intensity; W 0 is reference speed.
The coefficients δ, β and parameter σ 2 α can be calculated as δ=b+v, β= bv, σ 2 α = v 2 b σ 2 u /δ where b=V/L, V is the aircraft airspeed; L=200…1000 m is a scale of atmospheric turbulence; v=0.1…0.01 sec -1 is a parameter depending on the aircraft type and flight conditions; σ 2 u =0.4…2.7 m/sec is the RMS of wind speed fluctuation. By method of statistical tests the authors stated that for typical values δ=0.34 Hz, β=0.0044 sec, σ α =0.1…0.5 m/sec the K i,j (NSE) coefficient is within 0.003-0.001. The peculiarity of constructing the shortest TT using the graph theory algorithms is in the fact that during the flight the NSE and FTE errors can change under the influence of different factors. The authors suggest taking into account this effect for forming the graphs vertex weights in accordance with expression (2).
Let's consider the technique of graph forming for the case when the weights of its vertices are determined only with account of the NSE error. For this let's construct a GLONASS accuracy field in the given zone of airspace with account of forecast of its variation during the flight.
Distribution of PDOP values in space we will consider as a feature of the GLONASS accuracy field. Determination of set of points where the PDOP values are equal or within the given limits allows the SNS accuracy field to be constructed as areas with equal PDOP values [12]. The forecast of SNS accuracy field in any moment of time and for any airspace point is possible as PDOP is a changing but determined characteristic of SNS accuracy. The PDOP value depends on the aircraft dynamic position and the working constellation of navigation satellites which can be calculated for any moment of time and in any airspace point by SNS almanac data.

Software implementation
For conducting research in the LabVIEW graphic programming environment a software complex was developed ( fig. 2). Its structure is defined by the following peculiarities: the problem of track optimization is solved in the geodetic coordinate system (latitude B, longitude L, altitude above the ellipsoid surface H); the orbital movement of navigation satellites is specified in the rectangular geocentric coordinate system OXYZ; for constructing a graph the airspace points where PDOP is defined are projected on a plane with the help of Gauss-Kruger projection. It is necessary to provide the equality of distance between the adjacent graph vertices. The software complex includes: VI -a virtual instrument providing comfortable interface user-complex; 1 -a unit of forming Г(x, y) graph vertices; 2 -a unit of geodetic coupling of graph Г(B, L) to the geodetic coordinate system; 3 -a unit of almanac transformation; 4 -a unit of calculating orbital movement of navigation satellites; 5 -a unit of calculation of PDOP values in graph vertices; 6 -a unit of track construction; 7 -a unit of calculation of initial and boundary conditions for the constructed track; 8 -a unit of choosing an optimal track. Unit 1 forms a fixed grid (graph vertices) in Gauss-Kruger projection. In the graph vertices the PDOP value is calculated. Unit 2 is intended for coupling of the graph vertices to the earth's surface with account of its curvature. Unit 3 forms vectors from the GLONASS almanac (Lam is the longitude of the ascending node, dI is correction of the orbit inclination angle; w is the perigee argument; e is the eccentricity; dT is correction to the mean value of nodical period of satellite revolution; dTT is a speed of variation of orbiting period; Tw is the time of passing the first node) each having 24 elements according to the number of navigation satellites in the orbital group.
Unit 4 forms tree vectors (x,y,z) containing the coordinates of navigation satellites in the rectangular geocentric coordinate system for the i th moment of time. Unit 5 defines the navigation satellites in view, calculates the PDOP values in the graph vertices (so it forms a weighted graph). Unit 6 finds the optimal track in the weighted graph with account of limitations on the used airspace (prohibited zones).
Unit 7 specifies the size of the airspace zone, location and size of prohibited zones in it as well as flight path parameters (speed, altitude, heading).
The virtual instrument provides download of the GLONASS almanac from the site of GLONASS Information and Analysis Center (https://www.glonass-iac.ru/), setting of the parameters of aircraft motion and boundaries of the airspace zone, limitations on the track. The output data (the optimal track, the calculated value of TT length and etc) are displayed on the front panel of the virtual instrument.

Experimental results
The track can be optimized in a static (constructed for a specific moment of time) GLONASS accuracy field. However, the results of experiments show that the PDOP structure and values for static accuracy fields vary at different moments of time. During experiments the change was from 1.57 to 2.04. So the track, optimal in static field, may not be optimal due to its variation during the en-route flight. Under these conditions only the calculation of dynamic field will allow the information about real en-route time-referenced GLONASS accuracy to be received and an optimal 4D track to be constructed.
The dynamic field can be constructed using forecast of PDOP values along the constructible flight track for calculation time moments when flying en-route. Figure 3a represents a dynamic field constructed according to calculated PDOP values at the time moments corresponding to the aircraft flight along the DT. It has some colored areas corresponding to different ranges of PDOP variations within their borders.
On the basis of the calculated dynamic GLONASS accuracy field a graph is constructed with vertices having weights dependent on PDOP value distribution in time in the given area of airspace (fig.3b). a) b) Fig.3. Representation of a GLONASS accuracy field as a graph. The results of the experiments showed that the dynamic accuracy field is less uniform that the static one, the range of possible PDOP value variations in it being also wider than in a static field.
For solving the problem of optimal track construction the A-star algorithm was used. Its working capacity and efficiency in different conditions characterizing the airspace (absence or presence of prohibited areas in it) and GLONASS accuracy field (insignificant or essential PDOP variations) were researched.
According to the modelling results, while the GLONASS field is accurate enough and its variations are insignificant (PDOP>2.0), the optimal track coincides with a great circle. Besides, the TT length depends on the PDOP value and can vary within 2-3% depending on the mean PDOP m value during the en-route flight. Choosing the time of entering the route when the flight will be characterized by the minimal mean PDOP m , the pilot can obtain the reduction of TT length. The obtained results show that if the zone of «bad» PDOP values is small, the optimal track can coincide with a great circle ( fig.4 a). If the zone is relatively big, the optimal track differs from a great circle ( fig.4 b). In the example shown the TT length is 1,416.3 km for the great-circle flight and 1,399.5 km for the flight along the optimal track, that is almost 17 km (1.2 %) shorter.
Let's consider the problem of optimal track construction with prohibited zones in the airspace. Figure 5 represents the routes constructed in the GLONASS accuracy field with minor variations of PDOP values (1.2<PDOP<1.9) ( fig.5 a, b) and with a significant variation range (1.3<PDOP<10) ( fig.6 a, b) for different configurations of prohibited zones.
The conducted research shows that the result of track optimization depends both on the structure of the GLONASS accuracy field and on the size and configuration of prohibited zones in the provided airspace. Figure 7 represents the aircraft flight tracks for two cases: a track heuristically chosen by the crew for avoiding prohibited zones (track 1) and an optimal track (track 2).
The modelling results show that the TT length for the flight along the optimal track is 1,392.9 km and along the heuristically chosen track it is 1,451.3 km, that is almost 58 km (4%) longer than for the flight along the optimal track.

Discussion of the received results and conclusion
The obtained results demonstrate the efficiency of applying the graph theory algorisms for construction of a track with optimal TT length for different characteristics of a dynamic GLONASS accuracy field and in the presence of prohibited zones in the airspace.
The suggested new approach to constructing an optimal flight track makes it possible to take into account the track holding when navigating by GLONASS as well as presence of prohibited zones in the provided airspace. The new approach can be used both in the preliminary flight planning and during the flight itself provided the corresponding software is implemented in the onboard FMS. This will allow the situational awareness to be improved when making a decision about choosing the best flight route under varying air conditions and accuracy of air navigation support.
When using the optimization criterion proposed in the paper the choice of optimal flight track provides TT length reduction. The benefit in TT length reduction depends on characteristics of the GLONASS accuracy field and is quite significant for extended flight routes. So, for 5,000-km route the TT length is reduced by 60 km, that is equivalent to reduction of flight time approximately by 5 min. Even a slight reduction of flight time in context of high intensity of air traffic gives a significant economic effect for air carriers.