Computation of Viability Kernels on Grid Computers for Air- craft Control in Windshear

Article history: Received: 29 November, 2017 Accepted: 24 January, 2018 Online: 28 February, 2018


Introduction
Atmospheric conditions such as windshear continue to be considered as a source of potentially severe consequences. They are dangerous for aircraft during landing or take-off, because the wind gusts can occur at relatively low altitudes. Nevertheless, windshear is also dangerous during the cruise flight phase because it can lead to violation of the prescribed flight level.
In view of threats related to wind disturbances, there is permanent interest in designing robust aircraft guidance and control schemes (possibly for use with autopilots). The related question consists in finding safety domains, i.e. sets of initial states from which the control problem can be solved in the case of worst wind disturbance whose components lie in a known range.
There exist a large number of works devoted to the problem of aircraft control in the presence of severe windshears. In particular, papers [1][2][3][4][5][6] address the problem of aircraft control during take-off in windshear conditions. In works [1] and [2], the wind velocity field is assumed to be known. It is shown that open loop controls obtained as solutions of appropriate optimization problems provide satisfactory results for rather severe wind disturbances. Nevertheless, it is clear that the spatial distribution of wind velocity cannot be measured with appropriate accuracy, and therefore feedback principles of control design are more realistic. Different types of feedback controls are proposed in papers [3][4][5][6]. In [3], the design of a feedback robust control is based on the construction of an appropriate Lyapunov function. Robust control theory is used in [4] to develop feedback controls stabilizing the relative path inclination and (in [5] and [6]) for the design of feedback controls stabilizing the climb rate.
An approach based on differential game theory (see e.g. [7]) is presented in paper [8] in connection with the problem of landing. A high dimensional nonlinear system of dynamic equations is linearized and reduced to a two-dimensional differential game using a transformation of variables. The resulting differential game is numerically solved, and optimal feedback controls are constructed and tested in the nonlinear model against a downburst.
Another method, also based on differential game theory, is used in papers [9] and [10] to find feedback controls that are effective against downbursts. This approach assumes the computation of the value function, which is a viscosity solution (see. e.g. [11] and [12]) of an appropriate Hamilton-Jacobi equation. The numerical implementation utilizes dynamic programming techniques described in [13] and [14]. The case of known wind velocity field as well as the case of unknown wind disturbances are considered.
Recent investigations [15] and [16] utilize a method (see [17]) for the fast computation of rough approximations of solvability sets in linear conflict control problems. Using techniques of sequential linearization, this method is applied to nonlinear aircraft dynamics to design an appropriate control for take-off in the presence of downbursts.
The current paper addresses the problem of retaining trajectories in an appropriate flight domain (AFD) corresponding to the cruise flight phase (cf. [18] and [19]). Viability theory (see e.g. [20]) provides numerical methods (see e.g. [21], [22], and the Appendix) for finding the viability kernel, i.e. the largest set of initial states lying in the AFD from which viable trajectories emanate. More precisely, it includes all initial states for which there exists a feedback control that generates trajectories remaining in the viability kernel for all possible admissible wind gusts. In the case where the initial state does not belong to the viability kernel, there exists a method of designing a wind disturbance such that all trajectories violate the AFD for all possible controls.
As for wind conditions, it is assumed that only bounds on the wind velocity components are imposed. The dynamics of the aircraft will be considered as a differential game (cf. [9] and [10]) where the first player chooses control inputs, whereas the second player forms the worst wind disturbance. It is assumed that the first player is able to measure the current state vector, whereas the second player can measure both the current state vector and the current control ("future" values are not available) of the first player. Therefore, the second player may use the so-called feedback counter strategies (see [7]).
The current paper has common features with the open-access publication [23] concerning the model description and solution method. It should be noted that the publication [23] is mainly focused on theoretical fundamentals of the differential game approach. Regarding computational results, paper [23] formulates a problem of constructing the viability kernel in seven dimensions and performs several steps of the algorithm to show its feasibility. In contrast, the current paper addresses aspects of implementation on large scale grid computers and completely solves the above mentioned seven-dimensional problem including simulation of optimal trajectories. The paper is structured as follows: Section 2 outlines a point-mass aircraft model describing the vertical motion of a generic modern regional jet transport aircraft. The model is closely re-lated to the one described in paper [23]. The difference consists in a more clear method of deriving the dynamics equations.
In Section 3, state constraints related to the cruise flight phase are formulated, and the corresponding computed viability kernels are demonstrated through their three-dimensional sections. Additionally, trajectories yielded by an optimal feedback strategy, working against an optimal control of the disturbance, are shown. It should be noted that an optimal control of the disturbance can be constructed either as counter or pure feedback strategy because the so called saddle point condition holds for the differential game under consideration.
Section 4 outlines some aspects of parallel implementation of the computational method on a supercomputer system. The parallelization principles and data flow inside and between compute nodes are sketched. The novelty of our approach and comparison with existing software tools are discussed.
Section 6 (Appendix) briefly outlines the concept of differential games and viability kernels. Grid schemes for the computation of them are sketched. The details can be found in [23].

Model equations
In this section, a point-mass model representing the vertical motion of a generic modern regional jet transport aircraft is considered. Table 1 introduces euclidean coordinate systems (COS) that are necessary to compute the forces exerted on the aircraft. The origins of COSs are located either at the aircraft gravity center (CG) or at a fixed reference point (O) on the earth surface.
In positive direction of the symmetry axis of the turbine (aft looking forward) CG Body Fixed B In direction of the nose and in the symmetry plane of the aircraft CG Here, V K and V A are kinematic and aerodynamic aircraft velocities, respectively. The angles defining the relationship between the coordinate systems are the following (see also Figure 1): α K the kinematic angle of attack, α A the aerodynamic angle of attack, γ K the kinematic path inclination angle, σ the thrust inclination angle.
Parallel to the ground surface The position propagation is given in the local coordinate system (N ), whereas the translation dynamics are derived in the kinematic coordinate system (K). The model equations read as follows: Here, X T and Z T denote the components of the total force F T represented in the kinematic coordinate system (K), and the notation m stands for the aircraft mass. As usually, F T comprises aerodynamic, propulsion, and gravitation forces: Aerodynamic forces. They are defined as follows: where C D and C L are the drag and lift coefficients, respectively, ρ = ρ(h) is the air density (depends on the altitude), V A the aerodynamic velocity, and S the wing reference area. The lift and drag coefficients C D (α A , M) and C L (α A , M) are taken in the form: where the coefficients c D i and c L i , i = 1, ...9 are determined from least square fitting to experimental data.
The absolute value, V A , of the aerodynamic velocity can be derived using its relation to the kinematic velocity V K in the Local frame (N ) and the wind velocities W x and W z in the x N and z N directions, respectively. Therefore, and finally, using the matrix M N K , this implies the formula The aerodynamic angle of attack α A is computed as follows: where u A and w A are x B and z B -components of the aerodynamic velocity, respectively. The Mach number M is defined as follows: where c is the speed of sound, κ the adiabatic index for air, R the gas constant for ideal gases, and T (h) the temperature of air at the altitude h. See [23] for more details.
Propulsion Forces. Thrust forces are modeled considering a two-engine setup. Thus, where δ T ∈ [0, 1] is the thrust setting, and the functions f V and f γ are approximated similar to formulas (5) and (6), with δ T instead of α A . See [23] for more details.
Gravitation Force. For a simple gravitational model with constant acceleration g, the corresponding force is computed as: www.astesj.com The following two equations are added to the dynamics (1)-(4) to exclude jumps in the controls: Furthermore, the following equations produce smoothing of wind disturbances: (9) where the time constant, k w , is chosen as k w = 1 s −1 .

Problem setting and simulation results
Problem. The model consists of equations (2), (3), (4), (8), and (9). Thus, the state vector has seven variables: The following constraints on the controls and disturbances are prescribed: and the following state constraints are imposed: where h 0 = 10000 m being the cruise flight altitude. Additionally, the state constraints |W x | ≤ 5 m/s and |W z | ≤ 5 m/s hold automatically because of equations (9) and the constraints (10).
In  (19) were done. The computation has been done on the SuperMUC system at the Leibniz Supercomputing Centre of the Bavarian Academy of Sciences and Humanities. The computation was distributed over 100 compute nodes with 16 cores per node, which is regarded as "middle task" on the SuperMUC system. The runtime was about 40 h. Figures 2-4 show different three-dimensional sections of the seven-dimensional viability kernel. Figures 5-7 respectively demonstrate the same threedimensional sections and the corresponding projections of two optimal trajectories emanating from points lying in the viability kernel. The start point of trajectory 0 lies near to the boundary of the viability kernel, whereas trajectory 1 starts from a point lying deep inside of the viability kernel. The trajectories are computed when the control uses its optimal feedback strategy, and the disturbance utilizes its optimal counter feedback strategy (see the Appendix). The time step size used in the simulation of the trajectories was equal to 0.01. It is seen that the trajectories go to their attraction cycles and remain there. The simulation time interval is [0, 15 min]. Figure 8 shows a three-dimensional section of a smaller viability kernel corresponding to the shrinkage of the state constraints (11) by the factor 0.7. The corresponding three-dimensional projection of two optimal trajectories are shown. It is seen that the disturbance can keep trajectory 0 outside the reduced viability kernel because the start point lies outside of it.  Figure 2 and projections of two trajectories generated by an optimal feedback strategy of the control and an optimal counter feedback strategy of the disturbance. The start points are marked with bullets. Figure 6: The same section as in Figure 3 and projections of the same two trajectories as in Figure 5. Figure 7: The same section as in Figure 4 and projections of the same two trajectories as in Figure 5. 4 Implementation Aspects

Grid Computing, Parallelization and Scaling
We use a self developed software code to implement the algorithms of the grid scheme outlined in Appendix 6.3. The code is parallelized using a mixed MPI/OpenMP technique. The first two dimensions of the grid are decomposed, whereas the other dimensions remain unmodified. A compute node cartesian topology (see Figure 9) is then created such that each compute node corresponds to a grid cylinder. Additionally, each grid cylinder is supplied with ghost nodes that allow us to compute divided differences (18) used e.g. in the algorithm (19). In Figure 9, the www.astesj.com lines connecting the compute nodes show the data flow supported by MPI. The parallelization inside of each compute node (i.e. inside a grid cylinder) is supported by OpenMP. The problem considered in this paper deals with a seven-dimensional grid of 200×120×30×24×14×12×12 cells. Each of the first two dimensions was divided into 10 parts. Thus, there were 100 grid cylinders, each of size 20 × 12 × 30 × 24 × 14 × 12 × 12, plus the necessary ghost grid nodes. Therefore, the required memory per compute node was equal to about 8 GB.
Our observations show a good scaling behavior (see Table 2). The results were obtained for the problem of computing the viability kernel in five dimensions on a 200 × 120 × 30 × 24 × 14 grid. The relative speedup normalized to 32 cores and absolute timing of 30000 steps are shown.

Novelty of the method used
This paper utilizes a new method for computing viability kernels, which is based on the results of paper [21]. It is proven in [21] that the viability kernel is the Hausdorff limit of the sets {x : is the value function (see [7]) of a state constrained differential game . The numerical implementation of this method requires a theoretical basis and a stable numerical procedure for the treatment of transient Hamilton-Jacobi equations related to differential games with state constraints. Such a theoretical basis is given in paper [13] where the conventional conditions for viscosity solutions are modified to account for state constraints. This theoretical background is numerically implemented in [13] and [14], which results in monotone stable grid scheme (17) and (19). This enables us to perform a large number of time steps, usually several thousands.
Another new feature is that the corresponding software code is deeply parallelized using hybrid MPI/OpenMP techniques and adapted to run on a supercomputer system. Moreover, diverse modified variants of the code, based on sparse representations of grid functions, are tested.

Comparison with existing software
The well known software for solving Hamilton-Jacobi equations is the Level Set Method Toolbox (LSMT) described in [24]. This tool is really appropriate for solving rather general problems. However, it is not indicated in the manual whether the LSMT can compute viability kernels for differential games with state constraints. Moreover, according to the manual, the LSMT is not parallelized, whereas our software runs on a multi-node system.

Conclusion
This investigation shows that methods of differential games theory and viability theory can by applied to nonlinear aircraft models to investigate potential control abilities in the presence of wind disturbances. The new feature of our approach is the consideration of viability kernels for differential games. Feedback strategies of the players can be found from limiting grid "value functions" defining viability kernels (cf. Appendix). It is important that the amount of stored data is relatively small, which permits to implement the computed feedback strategies on a flight simulator. Further research will be focused on the treatment of models with more state variables, which will allow us to consider more realistic problems. Moreover, accounting for additional sources of uncertainty such as sensor errors or modeling uncertainties is planned.
Acknowledgment This work was supported by the DFG grants TU427/2-1 and HO4190/8-1. Computer resources for this project have been provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre under grant: pr74lu.

Appendix
This section briefly summarizes a method presented in [23] for computing viability kernels. Such a summary should help to provide a self-contained presentation.

Differential game
Consider a conflict control system with the autonomous dynamicṡ www.astesj.com Here; x stands for the state vector; u and v denote control inputs of the first and second players, respectively; and the compact sets P and Q describe constraints imposed on the control and disturbance variables, respectively. Further, it is assumed that all functions of x have global properties. For example, the right-hand side f is supposed to be bounded, continuous, and Lipschitzian in x on R n × P × Q.
In the following, it is assumed that the Isaacs saddle point condition (13) is true for all s ∈ R n and x ∈ R n . Note that this condition holds for the problem under consideration because controls and disturbances are additively separated in the model equations.

Viability kernel
For any v ∈ Q, consider the differential inclusioṅ Let G ⊂ R n be a compact set such that G = int G, this set will play the role of the state constraint. Let Definition 1 (u-stability property [7]) A set W ⊂ N is called u−stable on (−∞, T ] if for any position (t * , x * ) ∈ W , for any time instant t * ∈ [t * , T ], for any fixed v ∈ Q, there exists a solution x(·) of the differential inclusion (14) with the initial condition x(t * ) = x * such that (t * , x(t * )) ∈ W .
The next proposition is taken as a basis of the definition of viability kernels.
Proposition 1 [see [21] for the proof] Let W be a maxi- is nonempty, and W (t) → K in the Hausdorff metric as t → −∞. The set K is called the viability kernel of G and denoted by V iab(G).

Grid method for computing viability kernels
For the implementation of numerical method, it is convenient to represent viability kernels as level sets of an appropriate function. Let G λ be a family of state constraint sets defined by the relation where a continuous function g is chosen in such a way that, e.g., G 0 being the desired state constraint. It is necessary to construct a function V representing the viability kernels as follows: Such a function can be computed as a limiting solution, as t → −∞, of an appropriate Hamilton-Jacobi equation arising from conflict control problems with state constraints (see [13]). A grid approximation of V is computed as described below (cf. [21], [22], and [23]).
Let δ > 0 be a time step length, and the tuple h := (h 1 , ..., h n ) defines space step sizes. Set |h| := max{h 1 , ..., h n } and introduce the following upwind operator defined on grid functions related to the discretization h: Let {δ } be a sequence of positive reals such that δ → 0 and ∞ =0 δ = ∞. Consider the following grid scheme: where g h is the restriction of g to the grid defined by h.
It can be proven that V h monotonically converges point-wise to a grid function V h , and this function define approximations of the viability kernels according to formula (16), see [23] for more details. Remark 1 In (17), the operation min u∈P max v∈Q can be changed for max v∈Q min u∈P to obtain almost the same result. The difference tends to zero with |h|. The proof follows from the fact that the original operator (17) and the modified one satisfy the same consistency condition (see [13] www.astesj.com

Control design
This section outlines one of possible methods of control design. Consider the grid scheme (19) assuming that is large enough so that the required approximation is reached, i.e. V h +1 − V h L ∞ ≤ . The optimal control u and the worst disturbance v(u) at the current state x of the game can be found as solutions of the following program: Here, L h is an interpolation operator (see e.g. [14]) defined on the corresponding grid functions, and τ being a parameter which should be several times larger than the time step size of the simulation procedure to provide some stabilization. Note that the function V h can be transferred to a sparse grid (see e.g. [25] and [26]), which may essentially reduce the storage space. The disadvantage of such a technique is a certain loss of accuracy and a slower performance.
Remark 2 As it is described above, the second player uses the so-called feedback counter strategies, i.e. functions of x and u, where u being the current control action of the first player. Since the saddle point condition (13) holds, the theory of differential games says that the second player can achieve the same result using pure feedback strategies, i.e. functions of x. For example, a near optimal strategy of the second player can be obtained as a solution of the problem see also Remark 1. Thus, the both players achieve optimal results using pure feedback strategies.
Remark 3 If u and v appear linearly in the right-hand side of (12), then min and max operations in (17) and (20) can be computed only over extreme points of the sets P and Q respectively. This can be proven using the same arguments as in Remark 1.
Remark 4 If x 0 being a start point of the game, then x 0 lies in the approximate viability kernel defined as: and all trajectories approximately remain in this set. However, if the second player works non-optimally for a while, then, most likely, V h (x(t)) < V h (x 0 ) for somet > 0, and therefore, the state vector x(t) lies now in the smaller viability kernel {x ∈ R n : V h (x) ≤ V h (x(t))}.
Thus, faults of the second player improve the result of the first one.