Introduction

Tracking aircraft trajectories in adverse conditions, in particular under strong winds, can be considered a challenging task for aircraft control. This is especially critical during departure and the terminal phases of flight. First of all, these phases need to be performed in crowded airspace and ground proximity. Furthermore, they exhibit more dynamic maneuvers compared to other phases, such as cruise flight. Thus, robust tracking performance is vital for safe operation in these phases.

The investigation of methods for aircraft trajectory tracking has been the subject of numerous studies in the last decades (see e.g. [4,7,8,11,12,14,15]). In this paper, we follow a differential game based approach developed in [2,10]. It is particularly noteworthy that this approach can be applied without prior knowledge regarding the disturbance (wind). We show that the application of this algorithm in a realisitic flight simulator model is feasible and yields good tracking performance in rather severe wind conditions. The effectiveness and general applicability of this method is demonstrated for an approach trajectory which starts after the alignment turn from the holding pattern and terminates at the end of the glide-slope. In this study the reference trajectory is determined through the solution of an optimal control problem which minimizes the fuel consumption using direct optimal control methods (cf. [1,6]). This approach seems appealing as, even though, the primary goal of the trajectory tracking algorithm is to follow the reference as close as possible under adverse wind conditions the desired properties of the trajectory to be tracked are, to a certain extent, transferred to the control of the flight simulator.

This paper is structured as follows: First, the models used for the generation of optimal reference trajectories and the tracking algorithm are presented in Section 2. The following Section 3 describes the generation of a fuel optimal reference trajectory for the approach phase based on a multi-phase optimal control problem. In Section 4 the integration of the tracking algorithm with the NDI control architecture in a realistic flight simulator model is outlined. Finally, the approach is illustrated in Section 5 for tracking a reference trajectory in wind conditions.

Modeling

The dynamics implemented in the flight simulator model represent a realistic regional transport aircraft model. This model consists of a plant model with detailed aerodynamic and propulsion models as well as actuator and sensor models. For the controller we use the control architecture described in [5,9]. This architecture is based on a Nonlinear Dynamic Inversion (NDI) [13] controller with two cascaded loops, considering the attitude dynamics as middle loop and the rotation dynamics as the inner loop. Both loops feature reference dynamics of relative degree one which are modified by hedging signals and error controllers. These hedging signals represent the difference between the dynamics of the reference model and the expected reaction of the flight simulator model. It is noteworthy that in the implementation for this study we consider the Euler body angles as attitude states for the middle loop and the rotational body rates for the reference models of the innermost loop.

For the tracking algorithm and the generation of the optimal reference trajectory a reduced model is derived which contains a simplified model of the translational dynamics and the reference model of the middle (attitude) loop. The state vector of the simplified model

x = [x, y, h, V K , γ K , χ K , Φ RM , Θ RM , Ψ RM , T, M ] ,(1)

contains the position states g = [x, y, h] described in Cartesian coordinates, the translational states t = [V K , γ K , χ K ] described by the absolute kinematic velocity V K , the kinematic climb angle γ K , the kinematic course angle χ K , as well as the Euler angles Φ RM , Θ RM , and Ψ RM , the thrust state T , and the aircraft mass M . It should be noted that throughout this paper the symbol " " denotes transposition.

For the propagation of the position states we use the kinematic relations

  ẋ ẏ − ḣ  =   u K,O v K,O w K,O   = R OK (V K ) K =   cos(χ K ) cos(γ K ) sin(χ K ) cos(γ K ) − sin(γ K )   V K ,(2)

with (V K ) K = [V K , 0, 0] and the transformation matrix R OK from the kinematic (K) frame to the NED (O) frame. The propagation of the translational dynamics is expressed as

  VK χK γK   =   1 M 0 0 0 1 M V K cos(γ K ) 0 0 0 1 M V K   (F T ) K , (3)

using the total force vector denoted in the kinematic frame (F T ) K . The propagation of the attitude states is performed using first order models for the Euler angles and the dynamic model for the thrust is approximated by a first order lag (time constant τ T )

ΦRM = K Φ (V K , h)(u Φ − Φ RM ) (4) ΘRM = K Θ (V K , h)(u Θ − Θ RM ) (5) ΨRM = K Ψ (V K , h)(u Ψ − Ψ RM ) (6) Ṫ = τ −1 T (u T − T ),(7)

where the gains K Φ (V K , h), K Θ (V K , h), and K Ψ (V K , h) are scheduled over kinematic velocity V K and height h. Moreover, the control vector for the model is

defined as u = [u Φ , u Θ , u Ψ , u T ] .

It should be mentioned that from an engineering point it is advisable to schedule these gains over the aerodynamic velocity V A and height h (or similar quantities such as the dynamic and static pressure). However, note that the aerodynamic velocity depends on the disturbance inputs

v = (V W ) O = [(u W ) O , (v W ) O , (w W ) O ]

through the wind equation:

V A = R OK (V K ) K − (V W ) O . (8)

Here and below the notation • means the Euclidean norm. By taking the kinematic velocity instead of the aerodynamic velocity for scheduling the gains in the reference model this dependency is removed and the dynamics can be written in the form

d dt     g t c M     =     ġ(x) ṫ(x, v) ċ(x, u) Ṁ (x, v)     = f (x, u, v),(9)

with c = [Φ RM , Θ RM , Ψ RM , T ] , the disturbance vector v = (V W ) O , and the control vector u. Moreover, the time derivative of the aircraft's mass depends on the states and disturbances, i.e. is of the form d dt M = Ṁ (x, v), and models the fuel consumption by the engines. This dynamic equation is rather complicated as it includes tabulated data and depends on quantities such as the Mach number, thrust level, and height, which however are not directly influenced by the controls. Thus, the right-hand side is additive regarding the influence of the control and the disturbance vector

f (x, u, v) = f u (x, u) + f v (x, v),(10)

which implies that the Isaacs saddle point condition [10] is automatically fulfilled

min u∈P max v∈Q f (x, u, v) = max v∈Q min u∈P f (x, u, v). (11)

This property is desirable for the application of the differential game based tracking algorithm described in [2]. For the approach trajectory two different configurations are considered. The first configuration models the clean configuration which is used for the regular operation and is expressed as

ẋ = f c (x, u, v).(12)

In the landing configuration

ẋ = f l (x, u, v). (13)

the high lift devices and the gears (front and main gears) are deployed. This configuration is used in the terminal part of the approach. The model structure (9) remains the same in both configurations which implies that for both cases the saddle point condition (11) holds.

Generation of Optimal Reference Trajectories

In order to generate smooth reference trajectories an optimal control based approach is applied. This strategy is preferred as it does not only guarantee that the trajectory is in fact feasible regarding the aircraft dynamics but also allows for the consideration of a performance index and additional constraints. The reduced model introduced in Section 2 is employed to generate these trajectories.

The whole time interval of the problem P = t (i) , t (np) is split into n p = 4 time intervals (phases)

P (i) := t (i−1) , t (i) , i = 1, . . . , n p , (14)

with t (0) = 0 s and t (i) , i = 1, . . . , n p free. As the last phase P (np) represents the glide-slope, which is essentially pre-defined by the flight path angle and the approach velocity, this phase is not considered during the optimization. Thus, merely phases P (1) , . . . , P (np−1) are included in the optimal control problem and a final boundary condition to intercept the glide-slope is imposed.

The optimal control problem for the approach is defined as follows:

minimize u(t) ∈ U − M (t (3) ) subject to ẋ(t) − f c (x(t), u(t), 0) = 0, t ∈ 2 i=1 P (i) , ẋ(t) − f l (x(t)

, u(t), 0) = 0, t ∈ P (3) ,

φ i x(t (i) ) = 0, i = 0, . . . , n p − 1, −0.1 m ≤ h(t) − 1524 m ≤ 0.1 m, t ∈ 2 i=1 P (i) , −10 m ≤ h(t) − 1524 m ≤ 10 m, t ∈ P (3) , −10 • ≤ Θ t (0) ≤ 20 • , 0 % ≤ T t (0) ≤ 100 %. (15)

Note that the optimal control problem is formulated for the nominal case (without wind), i.e. v = 0. Moreover, the admissible control set U is defined by the following inequalities:

−25 • ≤ u Φ ≤ 25 • , (16)−10 • ≤ u Θ ≤ 20 • , (17) 0 % ≤ u T ≤ 100 %. (18)

At the horizontal start position [x, y] = 0 m the aircraft with initial mass M 0 is assumed to leave the alignment turn after the holding pattern with V K = 102.89 m/s at a course angle

χ K = 227 • φ 0 (x) = [x, y, V K − 102.89 m/s, χ K − 227 • , Φ RM , Ψ RM − 227 • , M − M 0 ] = 0,(19)

and heads over the waypoint defined by

φ 1 (x) = [x + 23839 m, y + 11197 m, χ K − 227 • ] = 0,(20)

to the next waypoint imposed through:

φ 2 (x) = [x + 25766 m, y + 73978 m, V K − 82.31 m/s, χ K − 137 • ] = 0. (21)

At this point, the aircraft switches from the clean configuration to the landing configuration. For the study considered in this paper this switch is assumed to happen instantaneously. In the last phase the aircraft in landing configuration intercepts the glide-slope with the following final boundary constraint

φ 3 (x) = [x + 22956 m, y + 12592 m, V K − 72.02 m/s, χ K − 82 • , γ K + 3 • ] = 0.(22)

The numerical solution of the optimal control problem is achieved by a direct method using the optimal control toolbox Falcon.m3 which discretizes the problem for each phase (index i) at time points: t (i) j , j = 1, . . . , N (i) = 10 3 . For the application under consideration a full discretization approach using the Backward Euler method is used which introduces equality constraints for the transcribed dynamics of the form

x (i) j+1 − x (i) j − (t j+1 − t j ) f c/l x (i) j+1 , u (i) j+1 , 0 = 0. (23)

Note that the optimal control problem (15) aims at a maximization of the final mass of the aircraft which is equivalent to a fuel minimization. The fuel minimal trajectory is depicted in Fig. 1. Observe, that the optimal trajectory is extended by the glide-slope which is assumed to be tracked with a velocity of V K = 72.02 m/s at an angle of γ K = −3 • . The reference trajectory for the position states as well as the velocity are approximated using n s piecewise polynomials. Let s(t) be the monotonically increasing distance covered

s(t) = t t (0) V K (η)dη. (24)

The optimal states x(t(s)) to be tracked are then represented by a continuously differentiable spline ζ : R → R 11 on intervals

S i,i+1 = [s i , s i+1 ] , i = 1, . . . , n s , s i = s(t i ) ζ(s) = a i,i+1 (s − s i ) 3 + b i,i+1 (s − s i ) 2 + c i,i+1 (s − s i ) + d i,i+1

, s ∈ S i,i+1 , (25)

P (1)

P (2) P ( 3)

P (4)

Fig. 1. Optimal approach trajectory (black, solid line) including the glide-slope (black, dashed line). Black dots mark the beginning and end of the phases in the optimal control problem.

which fulfills the following conditions at both endpoints of each interval:

ζ(s i ) = x(t i ) = d i,i+1 , (26)ζ(s i+1 ) = x(t i+1 ) = a i,i+1 ∆ 3 i,i+1 + b i,i+1 ∆ 2 i,i+1 + c∆ i,i+1 + d i,i+1 , (27) d ds ζ(s i ) = d ds x(t i ) = c i,i+1 ,(28)d ds ζ(s i+1 ) = d ds x(t i+1 ) = 3a i,i+1 ∆ 2 i,i+1 + 2b i,i+1 ∆ i,i+1 + c i,i+1 ,(29)

with ∆ i,i+1 = s i+1 − s i , the optimal state values x(t i ) and

d ds x(t i ) = f (x(t i ), û(t i ), 0) /V K (t i ). (30)

Integration of the Tracking Approach in the Flight Simulator Model

Within the closed-loop flight simulator model the tracking algorithm is executed at a sampling time δ u in the outermost control loop. At the beginning of each tracking step the value of the traveled distance s k at a base-point on the reference trajectory, which is defined as the closest point to the position states, needs to be determined. This point is defined by the relation

g(t) − ζ g (s) d ds ζ g (s) = 0,(31)

with

ζ g (s) = [ζ x (s), ζ y (s), ζ h (s)

] which may by solved numerically in each step, e.g. using Newton type iterations. In order to facilitate the search for the basepoint a predictor-corrector scheme is used. First, a predictor is created by propagating the base-point from the last time step

s k−1 sk = s k−1 + (t k − t k−1 ) d dt s(t k−1 ), (32) using d dt s(t) = ġ (t) d ds ζ g (s(t)) d ds ζ g (s(t)) d ds ζ g (s(t)) − g(t) − ζ g (s(t)) d 2 ds 2 ζ g (s(t)) . (33)

This predictor is corrected by at most 10 Newton-type steps using (31) starting with sk . The last iterate is then taken as the distance covered at the base-point s k in the current step.

Besides the current base-point

ζ k = ζ(s k ), the base-point ζ p = ζ(s p

) at the end of the prediction horizon with the prediction time-step δ p is required (see Alg. 1). Here, we assume that our control objective to closely follow the reference trajectory is, in fact, achievable which allows for the exact determination of the predicted base-point as s p = s(t k + δ p ).

The tracking algorithm uses a guide model

ẇ = f c/l (w, u, v),(34)

which has exactly the form of the reduced model and is initialized as

w 0 = ζ(0). (35)

This guide model represents an auxiliary model which first chooses a disturbance at the current sampling time in order to remain close to the states of the primary model x k . Then it chooses a control to aim at the reference trajectory in the end of the current sampling interval (cf. [2,10]). For the application to the flight simulator model it should be mentioned that the states x k of the primary model are taken as the measurements of the respective flight simulator states. The steps of the tracking algorithm are presented in Alg. 1.

For the application under consideration the admissible set Q of the wind

disturbances v = [(u W ) O , (v W ) O , (w W ) O ]

which are used to evaluate the maxoperator is defined by

(u W ) O ∈ {−5 m/s, 5 m/s}, (v W ) O ∈ {−5 m/s, 5 m/s}, (w W ) O ∈ {−5 m/s,w k+1 ← w k + δ u f c/l (x k , û, v * )

propagate w with û and v * 8: return w k+1 , û 9: end procedure Regarding the controls, the following procedure has shown to be effective for the application at hand: the input set P for each of the controls u i , i ∈ {Φ, Θ, Ψ, T } is restricted to values which are either the current value of the corresponding state (e.g. a good choice in case of a stationary flight condition), maximal or minimal dynamic limits ūi and u i , or values u − i and u + i leading to smaller incremental changes.

The dynamic limits are defined as

u i = max{u i,min , z i − δu i },(37) ūi= min{u i,max , z i + δu i },(38)

and are bounded from below and above by the absolute limits u i,min and u i,max . Note that these dynamic limits are relative to the current value z i of the measured states corresponding to the respective control quantity, i.e. the roll angle z Φ = Φ, pitch angle z Θ = Θ, yaw angle z Ψ = Ψ , and the thrust z T = T . The respective δ-values and the absolute (minimum and maximum) limits for all controls are defined in Table 1. Note that in Table 1 the minimum and maximum values for the roll angle (u Φ,min and u Φ,max ) differ from the limits imposed in the optimal control problem (cf. ( 16)). This is contributed to the fact that in the case with wind the roll angle is observed to saturate during simulations which does not occur in the nominal case. As such, the limits u Φ,min and u Φ,max of the roll angle are increased to ±30 • instead of ±25 • in order to provide additional control authority for compensating the wind disturbance. The incremental changes are defined based on these dynamic limits as

u − i = max{z i − (ū i − u i ) /10, u i }, (39)u + i = min{z i + (ū i − u i ) /10, ūi },(40)

i.e. at most a ±10% change regarding the dynamic limits w.r.t. the corresponding state value is allowed.

In particular, this choice of the controls appears to be less prone to the typical effect of control chattering, which is very often observed for differential Table 1. Control limits and increments

i u i,min u i,max δu i U nit Φ −30 30 15 • Θ −10 20 12 • Ψ −∞ ∞ 5 • T 0 100 6.4 %

game applications and is undesirable for most practical applications from an engineering viewpoint.

Let the optimal control returned from a step in the tracking procedure be denoted with û = [û Φ , ûΘ , ûΨ , ûT ] . The optimal thrust command ûT is directly fed to the engine model of the flight simulator. However, the optimal attitude commands ûΦ , ûΘ , and ûΨ generated by the tracking procedure are used as reference command signal r ΦΘΨ = [u Φ , u Θ , u Ψ ] for the middle loop of the NDI controller in the flight simulator model. This controller uses a modified reference model with the states ΦRM , ΘRM and ΨRM . The structure of the modified reference model in the attitude loop of this NDI controller is illustrated in Fig. For the propagation of the modified reference model states we use

   ΦRM ΘRM ΨRM    =   ν RM,Φ ν RM,Θ ν RM,Ψ   −   ν h,Φ ν h,Θ ν h,Ψ   ,(41) with  ν RM,Φ ν RM,Θ ν RM,Ψ   =   K Φ (h, V K )(u Φ − ΦRM ) K Θ (h, V K )(u Θ − ΘRM K Ψ (h, V K )(u Ψ − ΨRM )   . (42)

Here, the gains K Φ (h, V K ), K Θ (h, V K ) and K Ψ (h, V K ) are exactly the same gains used for the attitude dynamics in the reduced model (cf. ( 4)-( 6)). Moreover, the so-called hedging signals ν h,Φ , ν h,Θ , and ν h,Ψ are obtained as the difference between the pseudo commands ΦRM,e , ΘRM,e , and ΨRM,e as well as the measured time derivatives of the flight simulator states Φ, Θ, and Ψ :

  ν h,Φ ν h,Θ ν h,Ψ   =   ΦRM,e − Φ ΘRM,e − Θ ΨRM,e − Ψ   . (43)

The pseudo commands ΦRM,e , ΘRM,e and ΨRM,e in (43) are defined as the sum of the modified reference states dynamics and the error dynamics ν e,Φ , ν e,Θ and ν e,Ψ :

ν m =   ΦRM,e ΘRM,e ΨRM,e   =   ν RM,Φ ν RM,Θ ν RM,Ψ   +   ν e,Φ ν e,Θ ν e,Ψ   . (44)

The contribution from the error dynamics is defined as

  ν e,Φ ν e,Θ ν e,Ψ   =    K P e,Φ • ( ΦRM − Φ) + K I e,Φ • ( ΦRM − Φ)dt K P e,Θ • ( ΘRM − Θ) + K I e,Θ • ( ΘRM − Θ)dt K P e,Ψ • ( ΨRM − Ψ ) + K I e,Ψ • ( ΨRM − Ψ )dt    , (45)

where K P e,i i ∈ {Φ, Θ, Ψ } are the proportional and K I e,j j ∈ {Φ, Θ, Ψ } are the integral gains, respectively. The pseudo command ν m of the middle loop is then used in the inverted strap-down equation

r pqr =   1 0 − sin(Θ) 0 cos(Φ) sin(Φ) cos(Θ) 0 − sin(Φ) cos(Φ) cos(Θ)   ν m , (46)

to obtain the reference command r pqr for the rotational body rates p, q, and r. These commands are then fed to the inner loop, which has the same form illustrated in Fig. 2.

Step responses for the Euler angle commands are provided in Fig. 3. Step responses for the attitude states.

Illustrative Example

This section shows numerical results of the trajectory tracking approach using the differential game based tracking method for the realistic flight simulator model. The simulation was performed with the Euler forward method and a step size of δ s = 0.005 s. A sampling time of δ u = 0.02 s, prediction time of δ p = 2 s, and tolerance 0 = 0.01 was used to compute the optimal controls, see Alg. 1. For all states noisy measurements were used and the outputs of the flight simulator which need to be measured are estimated using an Extended Kalman Filter (EKF). The states considered for the tracking task are the position states x, y, and h as well as the kinematic velocity V K . A Dryden disturbance model [3] is used to test the implementation. The comparison of the reference and the flight simulator trajectory are depicted in Fig. 4 and the wind velocities are presented in Fig. 5. The values of the individual trajectories regarding the states to be tracked are shown in Fig. 6 and the deviations from the reference values are provided in Fig. 7. Moreover, a comparison of the optimal controls obtained from the tracking algorithm, the corresponding reference model states, as well as the flight simulator states are depicted in Fig. 8. All flight simulator states presented in these figures represent measured states using the EKF. The deviation in the height exhibits a prominent peak around t ≈ 300 s (see Fig. 7). This peak may be explained by the switch from the clean to the landing configuration which occurs approximately at this time point. Moreover, it is noteworthy that also the wind produced by the Dryden disturbance model exhibits its maximum absolute value of V W ≈ 20.19 m/s shortly before at 290.53 s. The deviations of the horizontal position states δx and δy are within ≈ 5 m and the deviation in the kinematic velocity δV K is within ≈ 3 m/s as can be seen in Fig. 7.

Conclusions and Outlook

In this paper it is shown that the differential game tracking approach developed in [2,10] can be applied in a realistic context for the robust trajectory tracking task in wind conditions. This trajectory tracking approach is integrated in the control architecture of a flight simulator model using reference models of relative degree one. For testing our implementation in realistic conditions an example problem for the approach flight phase featuring a switch from clean to landing configuration is considered. The wind disturbances imposed in this example are generated from a Dryden model and exhibit a maximum absolute value of V W ≈ 20.19 m/s. The numerical results of this numerical experiment indicate that the differential game based trajectory tracking approach can be employed for the tracking task under rather severe wind conditions. It should be noted that reference trajectory for this example problem is computed using a direct optimal control method. This approach for the generation of optimal reference trajectories is found to be very versatile for this task as it allows for the inclusion of features such as model switches (e.g. from clean to landing configuration) and the modeling of operational (path-)constraints in combination with typical performance indices (such as fuel consumption). Moreover, the properties of these optimal references implicitly transfer, to a certain extent, to the aircraft control through the robust tracking procedure.

Nevertheless, regarding future research it would be worthwhile to investigate if the tracking algorithm can be extended in such a way that constraints and performance indices can be included explicitly (e.g. in the min-operator computing the optimal control).