1 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.

2 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).

Table 1 Aircraft coordinate systems
Full size table

Here, and are kinematic and aerodynamic airfraft velocities, respectively.

The angles which define the relationship between the coordinate systems are the kinematic angle of attack \(\alpha _{K}\), the aerodynamic angle of attack \(\alpha _{A}\), the kinematic path inclination angle \(\gamma _{K}\), and the thrust inclination angle \(\sigma \) (see Fig. 1).

Fig. 1
figure 1

Aircraft coordinate systems (local, kinematic, aerodynamic, thrust, and body fixed) and associate angles. The corresponding, not shown, “z” axes are orthogonal to the “x” axes

Full size image

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, , and the position propagation is given in the local coordinate system (N):

$$\begin{aligned} m\dot{V}_K&= P_V - D - mg\sin \gamma _K, \end{aligned}$$
(1)
$$\begin{aligned} m V_K\dot{\gamma }_K&= P_{\gamma } + L - mg\cos \gamma _K, \end{aligned}$$
(2)
$$\begin{aligned} \dot{h}_N&= V_K\sin \gamma _K. \end{aligned}$$
(3)

In the above equations, \(P_V\) and \(P_{\gamma }\) 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:

$$\begin{aligned} \begin{pmatrix} D \\ L \\ \end{pmatrix} = \begin{bmatrix} \cos (\alpha _A-\alpha _K)&\quad -\sin (\alpha _A-\alpha _K) \\ \sin (\alpha _A-\alpha _K)&\quad \cos (\alpha _A-\alpha _K) \\ \end{bmatrix} \cdot \begin{pmatrix} \hat{D} \\ \hat{L} \\ \end{pmatrix}, \end{aligned}$$
(4)

where \(\hat{D}\) and \(\hat{L}\) are, respectively, the drag and lift forces in the aerodynamic reference system (A). They are given by the relations

$$\begin{aligned} \hat{D}&= \frac{1}{2}\,f_D(\alpha _A, M) \rho (h) \textit{SV}_A^2, \end{aligned}$$
(5)
$$\begin{aligned} \hat{L}&= \frac{1}{2}f_L(\alpha _A, M) \rho (h) \textit{SV}_A^2, \end{aligned}$$
(6)

where \(\alpha _{A}\) is the aerodynamic angle of attack, M the Mach number, \(\rho (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(\alpha _A, M)\) and \(f_L(\alpha _A, M)\) are taken in the form:

$$\begin{aligned} f_D(\alpha _A, M)&= c_1^D+c_2^D \alpha _A+ c_3^D M+c_4^D \alpha _A^2+c_5^D\alpha _A M\nonumber \\&\quad +c_6^D M^2+c_7^D\alpha _A^3 +c_8^D\alpha _A^2 M+c_9^D\alpha _A M^2, \end{aligned}$$
(7)
$$\begin{aligned} f_L(\alpha _A, M)&= c_1^L+c_2^L \alpha _A+c_3^L M+c_4^L \alpha _A^2+c_5^L\alpha _A M\nonumber \\&\quad +c_6^L M^2+c_7^L\alpha _A^3+c_8^L\alpha _A^2 M+c_9^L\alpha _A M^2, \end{aligned}$$
(8)

where the constants \(c_i^D\) and \(c_i^L\), \(i=1,\ldots ,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 considered in the local frame (N) and the wind velocities \(W_{x}\) and \(W_{h}\) in the \(x_N\)- and \(h_N\)-direction, respectively:

(9)

Therefore,

$$\begin{aligned} V_A=\left[ (V_K \cos \gamma _K -W_x)^2 + (V_K \sin \gamma _K -W_h)^2\right] ^{1/2}. \end{aligned}$$

The Mach number M is defined as the ratio between the absolute value of the aerodynamic velocity \(V_\mathrm{A}\) and the speed of sound c:

$$\begin{aligned} M=\frac{V_A}{c}, \quad c = \sqrt{\kappa \textit{RT}(h)}, \end{aligned}$$
(10)

with c depending on the adiabatic index for air, \(\kappa \), 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 \ldots 11\,\text {km}\), relative to sea level), the temperature of air, T(h), and air density, \(\rho (h)\), are approximated as follows:

$$\begin{aligned} T(h)&= T_\mathrm{s} \cdot \Bigl [1-\frac{n-1}{n}\frac{g}{R \cdot T_\mathrm{s}}\cdot H_\mathrm{G} \Bigr ], \end{aligned}$$
(11)
$$\begin{aligned} \rho (h)&= \rho _\mathrm{s} \cdot \Bigl [1-\frac{n-1}{n}\frac{g}{R \cdot T_\mathrm{s}}\cdot H_\mathrm{G} \Bigr ]^{\frac{1}{n-1}}, \end{aligned}$$
(12)
$$\begin{aligned} H_\mathrm{G}&= \frac{r_\mathrm{E}\cdot h}{r_\mathrm{E} + h}. \end{aligned}$$
(13)

In the formulas above, n is the polytropic exponent, g the gravitational constant, \(r_\mathrm{E}\) the earth radius, and \(T_\mathrm{s}\) and \(\rho _\mathrm{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 \(\sigma \) 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_{\gamma }\), appearing in Eqs. (1) and (2), are defined as follows in the kinematic frame (K):

$$\begin{aligned} \begin{pmatrix} P_V \\ P_{\gamma } \\ \end{pmatrix} = \begin{bmatrix} \cos \alpha _K&\quad \sin \alpha _K \\ -\sin \alpha _K&\quad \cos \alpha _K \\ \end{bmatrix} \cdot \begin{pmatrix} \hat{P}_V \\ \hat{P}_{\gamma } \\ \end{pmatrix}, \end{aligned}$$
(14)

where the components \(\hat{P}_V\) and \(\hat{P}_{\gamma }\) are computed by subtracting the drag components RD and CD from the gross thrust GT in the body-fixed frame (B):

$$\begin{aligned} \hat{P}_V&= 2 \cdot [\hbox {GT} \cdot \cos \sigma - (\hbox {RD} + \hbox {CD}) \cdot \cos \alpha _A], \end{aligned}$$
(15)
$$\begin{aligned} \hat{P}_{\gamma }&= 2 \cdot [ \hbox {GT} \cdot \sin \sigma - (\hbox {RD} + \hbox {CD}) \cdot \sin \alpha _A]. \end{aligned}$$
(16)

Here, the thrust components GT, RD, and CD are approximated by a second-order polynomial least squares fit depending on the thrust command \(\delta _T\in [0,1]\) and the Mach number M using the constants \(c_i^{\mathrm{GT}}\), \(c_i^{\mathrm{RD}}\), and \(c_i^{\mathrm{CD}}\) \(i=1,\ldots ,5\):

$$\begin{aligned} \hbox {GT}(\delta _T, M)&= c_1^{\mathrm{GT}}+c_2^{\mathrm{GT}} M + c_2^{\mathrm{GT}} \delta _T+c_3^{\mathrm{GT}} M^2 + c_4^{\mathrm{GT}} \delta _T M + c_5^{\mathrm{GT}} \delta _T^2, \end{aligned}$$
(17)
$$\begin{aligned} \hbox {RD}(\delta _T, M)&= c_1^{\mathrm{RD}}+c_2^{\mathrm{RD}} M + c_2^{\mathrm{RD}} \delta _T+c_3^{\mathrm{RD}} M^2 + c_4^{\mathrm{RD}} \delta _T M + c_5^{\mathrm{RD}} \delta _T^2, \end{aligned}$$
(18)
$$\begin{aligned} \hbox {CD}(\delta _T, M)&= c_1^{\mathrm{CD}}+c_2^{\mathrm{CD}} M + c_2^{\mathrm{CD}} \delta _T+c_3^{\mathrm{CD}} M^2 + c_4^{\mathrm{CD}} \delta _T M + c_5^{\mathrm{CD}} \delta _T^2. \end{aligned}$$
(19)

Introduce the following additional equations that smooth the controls:

$$\begin{aligned} \dot{\alpha }_K = \widetilde{\alpha }_K,\quad \dot{\delta }_T = \widetilde{\delta }_T. \end{aligned}$$
(20)

The following equations smooth the wind disturbances:

$$\begin{aligned} \dot{W_x} = -k_\mathsf{w} \left( W_x-\widetilde{W}_x\right) ,\quad \dot{W_h} = -k_\mathsf{w} \left( W_h-\widetilde{W}_h\right) ,\quad \text {with}\quad k_\mathsf{w}=1\,1/\mathrm{s}. \end{aligned}$$
(21)

3 Problem Statement and Simulation Results

Problem 1

The model consists of Eqs. (1)–(3) and (20). Thus, the state vector has five variables: \(V_K\), \(\gamma _K\), \(h_N\), \(\alpha _K\), and \(\delta _T\). The rate of the angle of attack, \(\widetilde{\alpha }_K\), and the rate of the thrust setting, \(\widetilde{\delta }_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:

$$\begin{aligned} \widetilde{\alpha }_K \in [-5,5]^{\circ }{/}\hbox {s},\quad \widetilde{\delta }_T\in [-0.3, 0.3]\,\text {1/s},\quad |W_x| \le 4\,\text{ m/s }, \quad |W_h| \le 4\,\text{ m/s }. \end{aligned}$$
(22)

Therefore, instantaneous changes in the angle of attack and thrust setting are not permitted. The following state constraints are imposed:

$$\begin{aligned} \begin{aligned}&V_K\in [100,170]\,\text{ m/s },\quad |\gamma _K|\le 20^{\circ },\quad h_N\in [-150,150]\,\hbox {m},\\&\alpha _K\in [0,16]^{\circ },\quad \delta _T\in [0.3,1]. \end{aligned} \end{aligned}$$
(23)

In the numerical construction, the box \([90,180]\times [-25,25]\times [-160,160]\times [-4,20]\times [0.2,1.2]\) of the space \((V_K, \gamma _K, h_N, \alpha _K, \delta _T)\) was divided in \(100\times 50\times 320\times 24\times 14\) grid cells, and the grid method (36) described in Sect. 4.1 is applied. The sequence of time steps, \(\{\delta _ \ell \}\), was chosen either as \(\delta _\ell =0.01/\ln (3 + \ell ^2)\) or \(\delta _\ell \equiv 0.01\), and the computations were performed until \(|\overline{\mathcal {V}}^h_{\ell +1} - \overline{\mathcal {V}}^h_\ell |\le \epsilon =10^{-6}\) for all grid nodes. About 30,000 time steps were carried out. In both cases of choosing \(\{\delta _\ell \}\), the results were identical. Moreover, the maximin viability kernel computed with the grid schema (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, \(\alpha _K=8\) and \(\delta _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:

$$\begin{aligned} W_x(t) = \left\{ \begin{array}{ll} 4,&{}\quad t\in [20,40],\\ -4,&{}\quad t\not \in [20,40], \end{array} \right. \quad W_h(t)=-W_x(t+15). \end{aligned}$$
(24)

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, \(\gamma _K=0\) and \(\delta _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, \(\gamma _K=0\) and \(\alpha _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.

Fig. 2
figure 2

Cross section (\(\alpha _K=8\) and \(\delta _T = 0.65\)) of the five-dimensional minimax viability kernel

Full size image
Fig. 3
figure 3

Cross section (\(\alpha _K=8\) and \(\delta _T = 0.65\)) of the five-dimensional minimax viability kernel. The trajectory (0–1) corresponds to an optimal feedback strategy of the control and an optimal feedback counted strategy of the disturbance. In the case of the trajectory (0–2), the impulse open-loop control (24) of the disturbance is used (Color figure online)

Full size image
Fig. 4
figure 4

Cross section (\(\gamma _K=0\) and \(\delta _T = 0.65\)) of the five-dimensional minimax viability kernel. The black and red trajectories correspond to the trajectories (0–1) and (0–2) from Fig. 3, respectively (Color figure online)

Full size image
Fig. 5
figure 5

Cross section (\(\gamma _K = 0\) and \(\alpha _K=8\) ) of the five-dimensional minimax viability kernel. The black and red trajectories correspond to the trajectories (0–1) and (0–2) from Fig. 3, respectively (Color figure online)

Full size image
Fig. 6
figure 6

Same cross section (\(\alpha _K=8\) and \(\delta _T = 0.65\)) of the five-dimensional minimax viability kernel as in 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)

Full size image

Problem 2

The model consists of Eqs. (1)–(3), (20), (21). Thus, the state vector has seven variables: \(V_K\), \(\gamma _K\), \(h_N\), \(\alpha _K\), \(\delta _T\), \(W_x\), and \(W_h\). The rate of the angle of attack, \(\widetilde{\alpha }_K\), and the rate of the thrust setting, \(\widetilde{\delta }_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 \(\widetilde{W}_x\) and \(\widetilde{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:

$$\begin{aligned} \widetilde{\alpha }_K \in [-5,5]^{\circ }{/}\hbox {s},\quad \widetilde{\delta }_T\in [-0.3, 0.3]\,\text {1/s},\quad |\widetilde{W}_x| \le 5\,\text{ m/s },\quad |\widetilde{W}_h| \le 5\,\text{ m/s }. \end{aligned}$$
(25)

The following state constraints are imposed:

$$\begin{aligned} \begin{aligned}&V_K\in [100,200]\,\text{ m/s },\quad |\gamma _K|\le 10^{\circ },\quad h_N\in [10,190]\,\hbox {m},\\&\alpha _K\in [0,16]^{\circ },\quad \delta _T\in [0.3,1]. \end{aligned} \end{aligned}$$
(26)

Notice that the state constraints \(|W_x| \le 5\,\text{ m/s }\) and \(|W_h| \le 5\,\text{ m/s }\) hold automatically because of Eq. (21) and the constraints (25).

In the numerical construction, the box \([90,210]\times [-15,15]\times [0,200]\times [-4,20]\times [0.2,1.2]\times [-6,6]\times [-6,6]\) of the space \((V_K, \gamma _K, h_N, \alpha _K, \delta _T, W_x, W_h)\) was divided in \(120\times 30\times 200\times 24\times 14\times 12\times 12\) grid cells, and 300 steps of the algorithm (36) were carried out. The same computer resources as in Problem 1 were involved. The runtime is about 1 h. The aim of this simulation is to show the feasibility of the computations in high dimensions. 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, \(\alpha _K=8\), \(\delta _T = 0.65\), \(W_x=0\), and \(W_h=0\).

Fig. 7
figure 7

Cross section (\(\alpha _K=8\), \(\delta _T = 0.65\), \(W_x=0\), and \(W_h=0\)) of a rough approximation of the minimax viability kernel. The computation was stopped after 300 steps of the algorithm (36)

Full size image

4 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 dynamics

$$\begin{aligned} \dot{x} = f(x,u,v),\quad x \in R^n,\quad u \in P\subset R^p,\quad v \in Q\subset R^q. \end{aligned}$$
(27)

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:

$$\begin{aligned} \min \limits _{u \in P}\max \limits _{v \in Q} \langle \ell , f(x,u,v) \rangle = \max \limits _{v \in Q} \min \limits _{u \in P} \langle \ell , f(x,u,v) \rangle ,\quad \ell \in R^n,\quad x \in R^n. \end{aligned}$$
(28)

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 \in Q\) the differential inclusion

$$\begin{aligned} \dot{x} \in F_v(x) = \overline{\text{ co }}\{f: f(x,u,v), u\in P\}. \end{aligned}$$
(29)

Let \(u\rightarrow v(u)\) be a Borel measurable function with values in Q. Consider the differential inclusion

$$\begin{aligned} \dot{x} \in F_{v(\cdot )}(x) = \overline{\text{ co }}\{f: f(x,u,v(u)), u\in P\}. \end{aligned}$$
(30)

Let \(G\subset R^n\) be a compact set satisfying the condition \(G=\overline{{\text {int}} G}\), T an arbitrary time instant, and \(N=(-\infty ,T]\times G\). For any subset \(W\subset N\) and any time instant \(t\le T\), denote \(W(t):=\{x\in R^n: (t,x)\in W\}\). The set G plays the role of state constraint.

The following definitions describe stability properties of subsets of N.

Definition 1

(Maximin u-Stability Property [16]) A set \(W \subset N\) is said to be u-stable in maximin sense on the interval \((-\infty ,T]\) if for any initial position \( (t_*,x_*) \in W\), for any time instant \(t^* \in [t_*, T] \), for any \(v \in Q\) there exists a solution \(x(\cdot )\) to the differential inclusion (29) with the initial state \(x(t_*) = x_*\) such that \((t^*, x(t^*)) \in W\).

Definition 2

(Minimax u-Stability Property [16]) A set \(W \subset N\) is said to be u-stable in minimax sense on the interval \((-\infty ,T]\) if for any initial position \( (t_*,x_*) \in W\), for any time instant \(t^* \in [t_*, T] \), for any Borel measurable function \(u\rightarrow v(u)\) with values in Q there exists a solution \(x(\cdot )\) to the differential inclusion (30) with the initial state \(x(t_*) = x_*\) such that \((t^*, x(t^*)) \in W\).

Properties of u-stable (both in maximin and minimax sense) sets.

  1. 1.

    The closure of an u-stable set is u-stable set.

  2. 2.

    The union of any family of u-stable sets is u-stable set.

  3. 3.

    There exists a unique (closed) maximal u-stable set if there is at least one u-stable subset of N.

  4. 4.

    If W is a maximal u-stable set, then \(W(t_1)\subset W(t_2)\) whenever \(t_1 \le t_2\le T\), and, additionally, \(W(T)=G\).

  5. 5.

    Remembering that there are two types of u-stability, the following holds for maximal u-stable sets: \(W^{\mathrm{mima}}\subset W^{\mathrm{mami}}\), where the upper indexes point out to minimax and maximin u-stability, respectively.

  6. 6.

    If the saddle point condition (28) holds, then \(W^{\mathrm{mima}} = W^{\mathrm{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=\overline{{\text {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=(-\infty ,T]\times G\). If for all \(t \le T\), then the set

$$\begin{aligned} K=\bigcap \limits _{t\le T} W(t) \end{aligned}$$

is non-empty, and \(W(t)\rightarrow K\) in the Hausdorff metric as \(t\rightarrow -\infty \). The set K is called the viability kernel of G and denoted by \(\mathrm{Viab}(G)\).

The next proposition also deals with both types of u-stability.

Proposition 2

(See [9] for the proof) The set \((-\infty ,T]\times \mathrm{Viab}(G)\) is u-stable.

The following propositions show the appropriateness of the above-given definitions.

Proposition 3

Let \(x_*\in \mathrm{Viab}^{\mathrm{mima}}(G)\), and \(\bar{t}>0\) be an arbitrary time instant. Then there exists a pure feedback strategy U(t, x) of the first player such that all trajectories generated by U and started at \(t=0\) from \(x_*\) remain in the set \(\mathrm{Viab}^{\mathrm{mima}}(G)\) for all \(t\in [0,\bar{t}]\). If \(x^*\not \in \mathrm{Viab}^{\mathrm{mima}}(G)\), then there exists a feedback counter-strategy \(V^c(t,x,u)\) and a time instant \(t_f\) such that all trajectories generated by \(V^c\) and started at \(t=0\) from \(x^*\) violate the state constraint G for \(t>t_f\).

Proposition 4

The same as Proposition 3 but with \(\mathrm{Viab}^{\mathrm{mami}}(G)\), feedback counter-strategy \(U^c(t,x,v)\), and feedback pure strategy V(t, x).

Proposition 5

The same as Proposition 3 but with \(\mathrm{Viab}(G):=\mathrm{Viab}^{\mathrm{mami}}(G)=\mathrm{Viab}^{\mathrm{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,\bar{t}]\times \hbox {Viab}^{\mathrm{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^*\not \in \hbox {Viab}^{\mathrm{mima}}(G)\), then there exists a time instant \(t^*\) such that \(x^*\not \in 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 open-loop 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).

4.1 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_\lambda \), is defined by the relation

$$\begin{aligned} G_\lambda =\{x\in R^n,\quad g(x)\le \lambda \}, \end{aligned}$$
(31)

where g is an appropriate continuous function. It is required to construct a function V representing the viability kernels as follows:

$$\begin{aligned} Viab(G_\lambda )=\left\{ x\in R^n,\quad V(x)\le \lambda \right\} . \end{aligned}$$
(32)

A grid approximation of such a function can be computed as a limiting solution, as \(t\rightarrow -\infty \), 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 \(\delta >0\) be a time step, and \(h:=(h_1,\ldots ,h_n)\) space discretization steps. Set \(|h|:=\max \{h_1,\ldots ,h_n\}\). Consider the following operators defined on grid functions:

$$\begin{aligned} \varPi ^{*}[\delta ,h;\phi ](x)= & {} \phi (x) + \delta \min \limits _{u \in P}\max \limits _{v \in Q}\sum _{i=1}^{n}\left( p_i^{ \text{ right }} f_i^+ + p_i^{ \text{ left }} f_i^-\right) , \end{aligned}$$
(33)
$$\begin{aligned} \varPi _{*}[\delta ,h;\phi ](x)= & {} \phi (x) + \delta \max \limits _{v \in Q}\min \limits _{u \in P}\sum _{i=1}^{n}\left( p_i^{ \text{ right }} f_i^+ + p_i^{ \text{ left }} f_i^-\right) , \end{aligned}$$
(34)

where \(f_i\) are the components of f(x, u, v), and

$$\begin{aligned} a^+= & {} \max {\{a,0}\},\quad a^-=\min {\{a,0}\},\nonumber \\ p_i^{ \text{ right }}= & {} \left[ \phi (x_1,\ldots ,x_i + h_i, \ldots , x_n)-\phi (x_1,\ldots ,x_i, \ldots , x_n)\right] /h_i,\nonumber \\ p_i^{ \text{ left }}= & {} \left[ \phi (x_1,\ldots ,x_i, \ldots , x_n)-\phi (x_1,\ldots ,x_i - h_i, \ldots , x_n)\right] /h_i. \end{aligned}$$
(35)

Let \(\{\delta _\ell \}\) be a sequence of positive reals such that \(\delta _\ell \rightarrow 0\) and \(\sum \nolimits _{\ell =0}^\infty \delta _\ell = \infty \). Consider the following grid schemes:

$$\begin{aligned} \overline{\mathcal {V}}_{\ell +1}^h= & {} \max {\left\{ \varPi ^*\left( \overline{\mathcal {V}}^h_\ell ;\delta _\ell \right) , g^h \right\} },\quad \overline{\mathcal {V}}^h_0 =g^h,\quad \ell = 0, 1, \ldots , \end{aligned}$$
(36)
$$\begin{aligned} \underline{\mathcal {V}}_{\ell +1}^h= & {} \max {\left\{ \varPi _*\left( \underline{\mathcal {V}}^h_\ell ;\delta _\ell \right) , g^h \right\} },\quad \underline{\mathcal {V}}^h_0 =g^h,\quad \ell = 0, 1, \dots , \end{aligned}$$
(37)

where \(g^h\) is the restriction of g to the grid.

It can be proved that \(\overline{\mathcal {V}}^h_\ell \) and \(\underline{\mathcal {V}}^h_\ell \) monotonically converge pointwise to grid functions \(\overline{V}^h\) and \(\underline{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 \(\delta _\ell /|h|\le (\sqrt{n} M)^{-1}\) should hold for all \(\ell \), where M denotes the bound of the right-hand side of (27). Moreover, the relation \(\delta _\ell /|h_\ell |\le (\sqrt{n} M)^{-1}\) should hold if the grid fineness \(|h_\ell |\) goes to zero too.

Remark 2

If the saddle point condition (28) holds, then

$$\begin{aligned} \lim _{\ell \rightarrow \infty } \overline{\mathcal {V}}^{h_\ell }_\ell = \lim _{\ell \rightarrow \infty } \underline{\mathcal {V}}^{h_\ell }_\ell = V. \end{aligned}$$

The proof follows from the fact that the both operators (33) and (34) satisfy the same consistency condition (see [7]) corresponding to the Hamiltonian

$$\begin{aligned} H(x,p):=\max _{v\in Q}\min _{u\in P}\langle p,f(x,u,v)\rangle = \min _{u\in P}\max _{v\in Q}\langle p,f(x,u,v)\rangle . \end{aligned}$$

Remark 3

Numerical simulations show that the grid functions \(\overline{\mathcal {V}}^h_\ell \) and \(\underline{\mathcal {V}}^h_\ell \) practically coincide for large \(\ell \), if the saddle point condition (28) is true.

4.2 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 \(\ell \) is sufficiently large so that the desired approximation is reached, i.e., \(|\overline{\mathcal {V}}^h_{\ell +1} - \overline{\mathcal {V}}^h_\ell |_{L_\infty }\le \epsilon \).

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:

$$\begin{aligned} u, v \rightarrow \min \limits _{u \in P}\max \limits _{v \in Q} \mathcal{L}^h\left[ \overline{\mathcal {V}}^h_\ell \right] \left( x+\delta f(x,u,v)\right) . \end{aligned}$$

Here, \(\mathcal{L}^h\) is an interpolation operator defined on the corresponding grid functions, and \(\delta \) 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. 1.

    The conventional multilinear interpolation operator (see, e.g., [8]) was used.

  2. 2.

    The grid function \(\overline{\mathcal {V}}^h_\ell \) was transferred to a sparse grid, using the conventional multilinear interpolation operator, and the operator \(\mathcal{L}^h\) was implemented as interpolation on a sparse grid (see, e.g., [30, 34]).

The first variant assumes that the whole grid function \(\overline{\mathcal {V}}^h_\ell \) 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 high-dimensional multiscale basis, which is a tensor product of one-dimensional multiscale bases. Thus, the advantage of the second variant is that the grid function \(\overline{\mathcal {V}}^h_\ell \) 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 counter-feedback strategies, respectively. If the game starts from a state \(x_0\), then the point \(x_0\) lies in the approximate viability kernel

$$\begin{aligned} \left\{ x\in R^n: \overline{\mathcal {V}}^h_\ell (x)\le \overline{\mathcal {V}}^h_\ell (x_0)\right\} , \end{aligned}$$

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) \(\overline{\mathcal {V}}^h_\ell (x(\bar{t})) < \overline{\mathcal {V}}^h_\ell (x_0)\) for some \(\bar{t}>0\), and therefore, the state vector \(x(\bar{t})\) lies now in the smaller viability kernel

$$\begin{aligned} \left\{ x\in R^n: \overline{\mathcal {V}}^h_\ell (x)\le \overline{\mathcal {V}}^h_\ell \left( x(\bar{t})\right) \right\} . \end{aligned}$$

Thus, faults of the second player improve the result of the first one.

5 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.