Aircraft Control During Cruise Flight in Windshear Conditions: Viability Approach

This paper addresses the analysis of aircraft control capabilities during the cruise phase (flying at the established level with practically constant configuration and speed) in the presence of windshears. The study uses a point-mass aircraft model describing flight in a vertical plane. The problem is formulated as a differential game against wind disturbances. The first player, autopilot, controls the angle of attack and the power setting, whereas the second player, wind, produces dangerous gusts. The state variables of the model are subjected to constraints expressing aircraft safety conditions. Namely, the altitude, path inclination, and velocity are constrained. Viability theory is used to find the so-called viability kernel, the maximal subset of the state constraint where the aircraft trajectories can remain arbitrary long if the first player utilizes an appropriate feedback control, and the second player generates any admissible disturbances. The computations are based on grid methods developed by the authors and implemented on a multiprocessor computer system.


Introduction
There is permanent interest in designing aircraft guidance schemes (possibly for use with autopilots) ensuring safe aircraft trajectories under conditions of windshear. In this connection, we can refer to papers [6,13,[18][19][20]22,23,29] where optimal control theory, robust control techniques, and differential game theory have been used to design appropriate controls. Both the case of known wind velocity field and the case of unknown wind disturbances have been considered.
While the papers quoted above deal with obtaining control strategies for survival enhancing, it is reasonable to address the question of the existence of such controls. In particular, it is meaningful to impose constraints on the state variables of the system to define an appropriate flight domain (AFD) and to study the question of the existence of controls keeping the system there (see [3,32]).
Such an approach is concerned with viability theory (see [1]). In recent time, there has been an essential progress in approximate computing viability kernels for control and conflict control problems. A comprehensive outline of abstract algorithms and their numerical implementations can be found in [26]. Traditionally, viability kernels are approximated by approaches proposed in [5,27,31]. Algorithms for computing approximate discriminating kernels are given in [11,12]. In particular, efficient numerical implementations are developed in the case of linear control problems (see [21]). Because of the relationship between reachable sets and viability kernels, it is easy to adopt level set algorithms based on Hamilton-Jacobi equations (see, e.g., [25]) to approximate these kernels. A realistic example of the application of such a technique is given in [28] where stabilization of a realistic vehicle is achieved using a viability approach.
The relation between differential games and viability kernels should also be mentioned. In particular, the current paper utilizes theoretical and numerical results on differential games. There are two approaches to numerical solving of differential games. The first one is related to the construction of the so-called maximal stable bridges introduced in [16]. The works [15,17] show typical techniques based on approximations of stable bridges. The second approach is concerned with viscosity solutions to Hamilton-Jacobi equations and grid-based algorithms. Here, the current state of scientific knowledge is well reflected in Refs. [2,4,14,24]. In papers [7,8], the Hamilton-Jacobi equation approach is extended to more general cost functionals.
In the current paper, a numerical method (see Sect. 4) for finding viability kernels, i.e., the largest sets of initial states lying in AFD from which viable trajectories emanate, is proposed. Moreover, feedback controls that produce trajectories remaining in the viability kernel for all possible admissible wind gusts are constructed. Thus, the dynamics of the aircraft is considered as a differential game, formalized according to [16], where the first player is associated with control inputs, whereas the second player forms the worst wind disturbance. The first player can measure the current state of the plant, whereas the second one has at his disposal both the current state and current control of the opponent ("future" values are not available). In other words, the second player uses the so-called feedback counter-strategies. The value function of this game is computed using a stable grid algorithm reported in [9,10]. Each level set of the value function converges to a viability (discriminating) kernel as the backward time tends to infinity. Since the first player is discriminated, the resulting viability kernels are called minimax ones. Taking into account that feedback counter-strategies and non-anticipative open-loop counter-controls of the second player yield the same result (cf. [33]), minimax viability kernels coincide with discriminating ones for Hamiltonians defined through the min u max v operation. The specific of minimax viability kernels is that they are defined without assuming the Isaacs (saddle point) condition. Therefore, the notation "minimax viability kernel" is used through the current paper. Similar arguments are true for maximin viability kernels related to feedback counter-strategies of the first player and pure feedback strategies of the second one.
The paper is organized as follows: Section 2 presents a point-mass model describing a generic modern regional jet transport aircraft flying in a vertical plane. The model is written in the kinematic reference system, and therefore, it does not contain the time derivatives of the wind components.
In Sect. 3, conflict control problems are formulated, state constraints corresponding to the cruise phase are imposed, and computed viability kernels are presented. Additionally, trajectories generated by an optimal feedback control, working against different disturbances, are shown. Section 4 outlines the concept of differential games and viability kernels. A grid algorithm for the computation thereof is sketched, and a method for the design of optimal feedback strategies is given.

Model Equations
This section introduces a point-mass aircraft model representing the flight of a modern generic regional jet transport aircraft in a vertical plane. To describe the influence of the relevant forces on the aircraft dynamics, the following euclidian coordinate systems (COS) with their origin either in the center of gravity of the aircraft (CG) or at a fixed reference point on the earth surface (O) are considered (see Table 1; Fig. 1).
Here, − → V K and − → V A are kinematic and aerodynamic airfraft velocities, respectively. The angles which define the relationship between the coordinate systems are the kinematic angle of attack α K , the aerodynamic angle of attack α A , the kinematic path inclination angle γ K , and the thrust inclination angle σ (see Fig. 1).
The translation dynamics are derived in the kinematic coordinate system (K ) for which the x K -axis points into the direction of the kinematic velocity of the aircraft, − → V K , and the position propagation is given in the local coordinate system (N ): In the above equations, P V and P γ denote the thrust forces, m is the aircraft mass, g the gravitational constant, and D and L represent the drag and lift forces, respectively. The last two variables are defined in the kinematic coordinate system (K ) as follows: whereD andL are, respectively, the drag and lift forces in the aerodynamic reference system (A). They are given by the relationŝ where α A is the aerodynamic angle of attack, M the Mach number, ρ(h) the air density at altitude h, S the wing area, and V A is the absolute value of the aerodynamic velocity. The lift and drag coefficients f D (α A , M) and f L (α A , M) are taken in the form: where the constants c D i and c L i , i = 1, . . . , 9 are found from least square fitting to experimental data. The absolute value of the aerodynamic velocity V A can be derived using its relation to the kinematic velocity − → V K considered in the local frame (N ) and the wind velocities W x and W h in the x N -and h N -direction, respectively: Therefore, The Mach number M is defined as the ratio between the absolute value of the aerodynamic velocity V A and the speed of sound c: with c depending on the adiabatic index for air, κ, the gas constant for ideal gases, R, and the temperature of air, T (h), at the altitude h. From the formulas of the International Standard Atmosphere ISA (DIN ISO 2533) for the troposphere layer (h = −2 . . . 11 km, relative to sea level), the temperature of air, T (h), and air density, ρ(h), are approximated as follows: In the formulas above, n is the polytropic exponent, g the gravitational constant, r E the earth radius, and T s and ρ s are the reference temperature and density of air, respectively. Note that, in all following state constraints, the altitude h N is defined relative to the altitude h = 5000 m.
For modeling thrust forces, a two engine setup with thrust inclination angle σ is considered. The following three main components influencing the thrust can be referred: the gross thrust (GT), the ram drag (RD), and the cowl drag (CD). The net thrust components P V and P γ , appearing in Eqs. (1) and (2), are defined as follows in the kinematic frame (K ): where the componentsP V andP γ are computed by subtracting the drag components RD and CD from the gross thrust GT in the body-fixed frame (B): Here, the thrust components GT, RD, and CD are approximated by a second-order polynomial least squares fit depending on the thrust command δ T ∈ [0, 1] and the Mach number M using Introduce the following additional equations that smooth the controls: The following equations smooth the wind disturbances:

Problem 1
The model consists of Eqs. (1)-(3) and (20). Thus, the state vector has five variables: V K , γ K , h N , α K , and δ T . The rate of the angle of attack, α K , and the rate of the thrust setting, δ T , are considered as controls, whereas W x and W h are regarded as disturbances; see formula (9). The following constraints are imposed on the controls and disturbances: Therefore, instantaneous changes in the angle of attack and thrust setting are not permitted.
The following state constraints are imposed:  (37) is only a few larger, which shows that the additional information about current values of the wind disturbances does not provide too much improvement of the control quality.
It should be noted that the constraints imposed on the disturbances, see (22), are near to critical. That is, the viability kernel is empty if the bound on the disturbances is larger than 5 m/s. Figure 2 shows the cross section of the minimax viability kernel if the last two state variables are fixed, α K = 8 and δ T = 0.65. Figure 3 shows the same set as in Fig. 2 and two optimal trajectories emanating from the same initial point lying in the viability kernel near to its boundary. The trajectories are computed when the first player (control) uses its optimal feedback strategy. The trajectory (0-1) corresponds to the case where the second player (disturbance) utilizes its optimal feedback counter-strategy. The trajectory (0-2) is computed when the disturbance is chosen as follows: The simulated flight time is equal to 30 min. The trajectories go to their attraction points and stay there. Figure 4 shows another cross section of the viability set. Namely, the second and fifth components of the state vector are fixed, γ K = 0 and δ T = 0.65. The projections of the same trajectories are shown: The black and red curves correspond to (0-1) and (0-2), respectively.
In the case of Fig. 5, the second and fourth components of the state vector are fixed, γ K = 0 and α K = 8. The projections of the same trajectories are shown. Figure 6 is analogous to Fig. 3. The difference is that the sparse grid approximation technique is used for the construction of an optimal feedback strategy of the control and an optimal feedback counter-strategy of the disturbance (see Sect. 4.2).
The computation was performed on the SuperMUC system at the Leibniz Supercomputing Centre of the Bavarian Academy of Sciences and Humanities. The problem was parallelized between 200 computers with 16 cores per computer and two threads per core. About 30,000 time steps were done. The runtime was about 1 h.

Problem 2
The model consists of Eqs. (1)-(3), (20), (21). Thus, the state vector has seven variables: V K , γ K , h N , α K , δ T , W x , and W h . The rate of the angle of attack, α K , and the rate of the thrust setting, δ T , are considered as controls. In the same way as in Problem 1, their instantaneous changes are not permitted. The disturbances are now associated with the artificial variables W x and W h that define the physical wind components W x and W h . Thus, the physical wind components do not exhibit instantaneous changes now.
The following constraints are imposed on the controls and disturbances:  Fig. 3. The difference is that the sparse grid approximation technique is used for the construction of an optimal feedback strategy of the control and an optimal feedback counter-strategy of the disturbance (Color figure online) The following state constraints are imposed:  It should be noted that the computer resources used are regarded as "middle task" on the SuperMUC system. Figure 7 shows the cross section of the computed seven-dimensional set if the last four state variables are fixed, α K = 8, δ T = 0.65, W x = 0, and W h = 0.

Differential Games and Viability Kernels
Let us shortly outline a method for computing viability kernels. The description will be given from the point of view of nonlinear differential games.
Consider a conflict control system with the autonomous dynamicṡ Here, x is the state vector; u and v are control parameters of the first and second players, respectively; and P and Q are compacts of the corresponding dimensions. In the following, it is assumed that all functions of x are defined on the whole R n and have global regularity properties. Thus, the right-hand side f is supposed to be globally bounded, continuous in (x, u, v), and Lipschitzian in x.
The following relation is called saddle point condition: Note that this condition does not hold for Problem 1 but holds for Problem 2 in Sect. 3. Bearing in mind the conflict control system (27), consider for any v ∈ Q the differential inclusionẋ Let u → v(u) be a Borel measurable function with values in Q. Consider the differential inclusionẋ Let G ⊂ R n be a compact set satisfying the condition G = int G, T an arbitrary time instant, and N = (−∞, T ] × G. For any subset W ⊂ N and any time instant t ≤ T , denote W (t) := {x ∈ R n : (t, x) ∈ W }. The set G plays the role of state constraint.
The following definitions describe stability properties of subsets of N . [16]) A set W ⊂ N is said to be u-stable in maximin sense on the interval (−∞, T ] if for any initial position (t * , x * ) ∈ W , for any time instant t * ∈ [t * , T ], for any v ∈ Q there exists a solution x(·) to the differential inclusion (29) with the initial state x(t * ) = x * such that (t * , x(t * )) ∈ W . [16]) A set W ⊂ N is said to be u-stable in minimax sense on the interval (−∞, T ] if for any initial position (t * , x * ) ∈ W , for any time instant t * ∈ [t * , T ], for any Borel measurable function u → v(u) with values in Q there exists a solution x(·) to the differential inclusion (30) with the initial state

Definition 2 (Minimax u-Stability Property
Properties of u-stable (both in maximin and minimax sense) sets.
1. The closure of an u-stable set is u-stable set.
2. The union of any family of u-stable sets is u-stable set.
3. There exists a unique (closed) maximal u-stable set if there is at least one u-stable subset of N . 4. If W is a maximal u-stable set, then W (t 1 ) ⊂ W (t 2 ) whenever t 1 ≤ t 2 ≤ T , and, additionally, W (T ) = G. 5. Remembering that there are two types of u-stability, the following holds for maximal ustable sets: W mima ⊂ W mami , where the upper indexes point out to minimax and maximin u-stability, respectively. 6. If the saddle point condition (28) holds, then W mima = W mami . Property 1 follows from properties of the differential inclusions (29) and (30); see, e.g., [9] for the proof. Property 2 immediately follows from the definition of u-stability. Property 3 is the consequence of the previous two properties. The first part of property 4 is proved in [9], and the last part follows from the continuous dependency of the solution sets of the inclusions (29) and (30) on the initial state and from the condition G = int G. The last two properties are true due to the following climes: the minimax u-stability property implies the maximin one, and they are equivalent if the saddle point condition holds.
The next proposition gives rise to the definition of viability kernels. The proposition deals with both types of u-stability.
Proposition 1 (See [9] for the proof) Let W be a maximal u-stable subset of N = (−∞, T ]× G. If W (t) = / for all t ≤ T , then the set is non-empty, and W (t) → K in the Hausdorff metric as t → −∞. The set K is called the viability kernel of G and denoted by Viab(G).
The next proposition also deals with both types of u-stability.
The following propositions show the appropriateness of the above-given definitions.  U c (t, x, v), and feedback pure strategy V (t, x).

Proposition 5 The same as Proposition 3 but with Viab(G) := Viab mami (G) = Viab mima (G) and feedback pure strategies U (t, x) and V (t, x).
The proof of these propositions is based on results of the monograph [16]. Consider, for example, Proposition 3. Note that the set [0,t]×Viab mima (G) is also u-stable in the minimax sense, and therefore, the first part of Proposition 3 is exactly a result of [16,Chapter 10].
If now x * / ∈ Viab mima (G), then there exists a time instant t * such that x * / ∈ W (t * ), where W is the maximal u-stable in the minimax sense subset of N . According to [16], there exists a feedback counter-strategy V c such that all trajectories generated by V c and started from x * at t = t * remain outside of a closed neighborhood of W . Since W (T ) = G (see Property 4), all trajectories violate the state constraint G at t = T . Since the system (27) is autonomous, we can set t * = 0 and t f = T − t * , which proves the second part of Proposition 3.

Remark 1
It should be mentioned that the notion of maximin and minimax viability kernels is introduced in this paper to modify the concept of discriminating kernel (see [11,12]) in the case where the saddle point condition (28) does not hold. According to the approach of [16], adopted in the current paper, a player is discriminated by prescribing him the usage of feedback strategies against feedback counter-strategies of his opponent. Therefore, each of the players can be discriminated, which results in maximin and minimax viability kernels. As it is proved in [16], the discrimination does not play any role if the saddle point condition (28) holds. In this case, each of the players cannot improve his result by using feedback counter-strategies instead of pure feedback controls, and therefore, maximin and minimax viability kernels coincide in this case.
Taking into account that the use of feedback counter-strategies and non-anticipative openloop counter-controls of the players yields the same result, maximin and minimax viability kernels can be referred as discriminating kernels for the Hamiltonians defined through the max v min u and min u max v operation, respectively. The specific is that maximin and minimax viability kernels are defined without requiring the saddle point condition (28).

Grid Method for Computing Viability Kernels
The numerical method utilizes the idea of representation of viability kernels as level sets of appropriate functions. Suppose that a family of state constraint sets, G λ , is defined by the relation where g is an appropriate continuous function. It is required to construct a function V representing the viability kernels as follows: A grid approximation of 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 [7]). The algorithm looks as follows (cf. [9,10]).
Let δ > 0 be a time step, and h := (h 1 , . . . , h n ) space discretization steps. Set |h| := max{h 1 , . . . , h n }. Consider the following operators defined on grid functions: where f i are the components of f (x, u, v), and a + = max {a, 0}, a − = min {a, 0}, Let {δ } be a sequence of positive reals such that δ → 0 and ∞ =0 δ = ∞. Consider the following grid schemes: where g h is the restriction of g to the grid. It can be proved that V h and V h monotonically converge pointwise to grid functions V h and V h , respectively. These functions define approximations of the viability kernels according to formula (32) (see [9,10] for more details). Note that the relation δ /|h| ≤ ( √ nM) −1 should hold for all , where M denotes the bound of the right-hand side of (27). Moreover, the relation δ /|h | ≤ ( √ nM) −1 should hold if the grid fineness |h | goes to zero too.

Remark 2
If the saddle point condition (28) holds, then The proof follows from the fact that the both operators (33) and (34) satisfy the same consistency condition (see [7]) corresponding to the Hamiltonian

Remark 3 Numerical simulations show that the grid functions V h and V h practically coincide
for large , if the saddle point condition (28) is true.

Control Design
This section outlines one of the possible methods (cf. [10]) of control design. Consider the case of Proposition 3 where the first player uses pure feedback strategies, whereas the second player utilizes feedback counter-controls. Consider the grid scheme (36) assuming that is sufficiently large so that the desired approximation is reached, i.e., |V Let x be a current state of the game. The control u and the disturbance v(u) are being found as solutions of the following problem: f (x, u, v)) .
Here, L h is an interpolation operator defined on the corresponding grid functions, and δ is a parameter which should be larger than the time step size of the control scheme to provide some stabilization.
The following two variants of the interpolation operator were tested: 1. The conventional multilinear interpolation operator (see, e.g., [8]) was used.
2. The grid function V h was transferred to a sparse grid, using the conventional multilinear interpolation operator, and the operator L h was implemented as interpolation on a sparse grid (see, e.g., [30,34]).
The first variant assumes that the whole grid function V h is stored on a hard disk, which can require several gigabytes of memory. On the other hand, the interpolation process runs relatively quickly in this case. The sparse grid method, introduced in [34], is based on the construction of a highdimensional multiscale basis, which is a tensor product of one-dimensional multiscale bases. Thus, the advantage of the second variant is that the grid function V h may be strongly compressed, so that the required disk space is being reduced to several tenth of megabytes. The disadvantage of this method is its slower performance.
Remark 4 Assume that the first and second players use the above feedback and counterfeedback strategies, respectively. If the game starts from a state x 0 , then the point x 0 lies in the approximate viability kernel and all trajectories remain in this set. However, if one of the players, say the second one, 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 Thus, faults of the second player improve the result of the first one.

Conclusion
The current investigation shows that methods of viability theory can by applied to realistic models of aircraft to evaluate its potential control ability in the presence of windshear. The new feature of this approach is the use of viability kernels arising from differential games. In particular, we define and compute viability kernels related to feedback and counter-feedback strategies of the players in the case where the saddle point condition does not hold. Such strategies can be computed and stored in a sparse form, which enables us to integrate them into realistic control systems.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.