Economic Model-Predictive Control for Aircraft Forced Landing: Framework and Two-Level Implementation

This article presents an optimization-based control framework for the autonomous forced landing of a fixed-wing unmanned aircraft (UA). A two-level model-predictive control (MPC) scheme is proposed to realize this framework, where an economic model-predictive control (EMPC) in a long piecewise constant fashion is proposed at the high level, while a short fixed-horizon linear time-varying MPC at the low level responds to fast dynamics of UA and tracks the optimal path provided by the high-level controller, alleviating computational burden compared to the high-frequency single-layer MPC scheme. Compared with a single EMPC setup with high sampling frequency, this hierarchical EMPC controller can significantly reduce the computational complexity and make it feasible to be implemented in real time. In addition, it also responds to disturbances (e.g., wind) and aircraft fast dynamics in a timely manner. The recursive feasibility and stability of the high- and low-level MPC schemes are established. The performance of the proposed EMPC forced landing function is illustrated by simulation case studies on both Aerosonde and Skywalker X8, compared favorably with competing schemes.

COEFFICIENTS C L Aerodynamic lift force coefficient. C D Aerodynamic drag force coefficient. C M Coefficient of pitch moment. C y Aerodynamic lateral force coefficient. C l Coefficient of roll moment. C n Coefficient of yaw moment.

I. INTRODUCTION
In recent years, unmanned aircraft (UA) has received considerable attention and development in both military and civilian applications due to its unique features: UA can significantly reduce the cost for training skilled pilots and guarantee human pilots' safety for long-range remote sensing, delivery, and surveillance missions in dangerous areas.
However, without pilots, many safety issues arise. Among many others, an important consideration is how a UA deals with engine failure and other critical failures and carries out an emergency landing. As fixed-wing UA systems are primarily manufactured with a single engine, they are as vulnerable to lose power (due to engine failure) as any single-engine General Aviation aircraft [1]. Therefore, with an increased presence of UA over civilian populations and infrastructure assets, a reliable and robust system for forced landing is crucial in the deployment of UA systems to reduce the risk of casualties and asset losses. This is regarded as a bottleneck in the widespread applications of UA. Furthermore, autonomous forced landing systems can also be served as an advanced assistance function for human pilots in emergency situations.
The landing site selection and the associated reachability have been studied in previous research (see, e.g., [2]- [4]). A human forced landing approach technique named "high-key low-key" technique intending to give a human pilot good visibility to study the site and time to prepare the aircraft for landing has been introduced in [5], but this technique is not particularly suitable for UA especially in an environment with obstacles and low altitude. In [3] and [6]- [9], Dubins curve techniques, which would be reshaped to trochoids if the effects of uniform wind are taken into account, were developed for forced landing, where the height loss is minimized by flying the shortest path. However, these geometric curves might not be attainable for the aircraft as some dynamic behavior and constraints were not taken into account; for example, the desired bank angle was assumed to be achieved instantaneously. The authors of [10] and [11] assumed simplified kinematic modeling of the aircraft and studied the trajectory generation and optimization regardless of the analysis on how to design control actions to drive the aircraft. When the engine loses power, such kinematic models are not adequate to cope with maneuverability constraints imposed by loss of thrust.
Our objective is to find the optimal solution such that the aircraft reaches the target region above the landing site with minimum height loss. This will either increase the reachable range during the aircraft gliding or increase the safe margin for performing a landing exercise. To improve the autonomous landing performance of a fixedwing UA, this article first proposes an optimal control framework by introducing the idea of economic modelpredictive control (EMPC). Unlike many standard applications in industrial process control, which decompose the plant's management and optimization into two levels (the real-time optimization (RTO) layer and the control layer), the EMPC looks into how to address economic optimization directly in real time (see [12]- [16] and references therein). The height loss of reaching a target region around a landing site is selected as the only criterion and consequently used as the objective function in online optimization.
For applications like UA manoeuvring, in particular forced landing, the control objective is to reach a specified target region, in which it does not necessarily contain equilibrium points, while minimizing the height loss. This is different from the objective in [17], which is to recover a fixed-wing aircraft to a steady-state altitude from a deep stall using an EMPC controller. Hence, the notion of stability is compromised by that of completion [18]: once the target region is reached (at completion), the problem ends. In addition, for most applications of completion problems, it is desirable to have a finite completion time. Variable-horizon model-predictive control (MPC) [19], [20] and shrinkinghorizon MPC [21], [22] have been shown to provide these features, and both types of MPC algorithms are able to relax the requirements on terminal ingredients 1 that are central to the stability proofs. Since the variable horizon approach introduces an additional integer decision variable in the constrained optimization problem solved at each iteration, which might become intolerable for fast dynamic systems, the variable horizon is not employed in this article. In order to further alleviate the computational intensity of nonlinear optimization, a control framework employing a two-level control structure has been proposed in [25]. The high level adopted a piecewise constant input EMPC to generate baseline control commands by exploiting the nonlinear aircraft model, while the low level was designed based on linearized models around the reference trajectories provided by the high-level controller.
By integrating the advantages of shrinking horizon and piecewise constant control, this article proposes a two-leveloptimization-based controller for a fixed-wing UA in forced landing applications. The main contributions of this article include the following: 1) For the first time, the forced landing problem is formulated in an EMPC framework, where the height loss and a target region around a forced landing site can be explicitly taken into account in the sense of optimal control. This study promotes a wider application of EMPC, particularly in aerospace sectors. Similar to the early development of MPC, most of the current applications of EMPC are in the process industry. 2) A two-level MPC scheme is proposed to implement this framework in trading off between many key factors, e.g., computational complexity, optimality, fast dynamics, and disturbance rejection. In this implementation, based on a simplified dynamic model, the high-level EMPC adopts shrinking horizon, piecewise constant control, and pure economic cost function so that it significantly reduces the workload of the online optimization; a fast linear time-varying MPC within a fixed horizon has been designed in order to reduce the effect of disturbances and uncertainties using its inherent robustness. 2 With the optimal control actions determined by the two-level controller, numerical implementations have been performed on a high-fidelity six-degree-of-freedom aircraft model. 1 Terminal ingredients include terminal constraint set, terminal cost, and local controller [23], [24]. 2 Inherent robustness refers to stability properties of a perturbed system in a closed loop with a nominal MPC design ignoring disturbances [24], [26], through repeated measurement. Meanwhile, this adoption of lowlevel MPC can also mitigate the demand of fast sampling frequency in high level while ensuring the reachability to the landing site.
3) Theoretic properties of the proposed hierarchical control scheme have been established, including recursive feasibility and closed-loop stability.
The rest of this article is organized as follows. Section II presents the mathematical model to describe the 3-D motion of a fixed-wing UA, the landing site with its associated target region, and the formal problem formulation. In Section III, the details of the proposed two-level MPC-based controllers are given, including the design of high-level shrinking horizon piecewise constant EMPC and low-level fixed horizon linear time-varying MPC. The closed-loop feasibility and stability results are analyzed in Section IV. Section V demonstrates the performance of the controller through illustrative examples. Finally, Section VI concludes this article.

A. Aircraft Modeling
This article aims to develop an optimal control framework for fixed-wing UA systems. Considering a righthanded coordinate system that coincides with the local eastnorth-up coordinates, specifically, we define the positive x-axis to be east, the positive y-axis to be north, and the positive z-axis to be up.
In the simulation, the aircraft dynamic model is described by high-dimensional nonlinear equations formulated in [27] and [28]. Nevertheless, to predict the long future aircraft trajectories within an acceptable computational complexity, a simplified aircraft dynamic model is employed to generate command guidance that drives the aircraft toward the desired region of a landing site. Specifically, the following set of equations [29]- [31] are used to describe the 3-D motion for a gliding aircraft over a flat Earth:Ẋ where state vector x = [X Y Z V γ ψ α φ] T and input command u = [ᾱφ]. The state variables X , Y , and Z denote the (x, y, z) position in a global coordinate system, V represents the airspeed of the aircraft, γ is the vertical path angle, ψ denotes the heading angle, α is the angle of attack, and φ is the roll angle. The control variablesᾱ andφ are the command angle of attack and roll angle, respectively. In particular, by adopting the techniques in [32] and [33], we employ the last two first-order differential equations in (1) to describe the dynamics of angle of attack and roll angle, relative to the piecewise constant input command u, where k α and k φ are constant values. In short notation, the nonlinear representation is denoted bẏ Normally, the aerodynamic forces and moments in both longitudinal and lateral are determined by many factors. The full expression of aerodynamic functions is presented in the book [27]. For the dynamic model (1), the lift force L and drag force D are simplified, which are expressed as where C L (α) and C D (α) are nondimensional aerodynamic coefficients. For small angles of attack, these coefficients can be modeled with acceptable accuracy using linear approximations and ρ is the air density, S is the wing area, and AR = b 2 S is the wing aspect ratio, with b denoting the wingspan.
The aircraft state variables are subject to physical and safety constraints that affect its maneuverability, which should be taken into account in the optimization problem. Mathematically, the feasible sets of system states and control inputs are denoted by x ∈ X and u ∈ U, respectively. Letting x t and u t denote the instantaneous values of the state variable x and input variable u at time instant t. Then, the (possibly coupled) pointwise-in-time state and input constraints to be considered are the following: where Z ⊆ X × U is assumed to be compact. In particular, a realistic form of the aircraft constraints is explicitly formulated as whereγ V andψV are the vertical and horizontal structural load limitations of the acceleration components, respectively. The constraint (X, Y, Z ) ≤ 0 typically represents safe distance above ground. Moreover, it is worth to notice that the constraint on the angle of attack α also relates to the accuracy of the aircraft model because the force parameters may be valid only in an interval of α.

B. Target Region
The objective of forced landing is to drive the aircraft to a specific destination. In the literature, quite often the destination landing site was simply regarded as a single point which might not be appropriate in many practical scenarios. In practical applications, it is not necessary to ask the UA to land on a point. Instead, a corridor as shown in Fig. 1 is acceptable and so considered in this article.
To save height and design an achievable flight path, this article describes a predefined landing site in the 3-D space. Let us denote the two ending locations of the runway by as shown in Fig. 1. Then, the destination can be expressed as a vector S = S end − S str . Accordingly, the horizontal angle of the landing site is calculated as Assuming that the required final straight-line landing distance is L S and the height should not be below H S , the critical point P L = [XȲZ] T , where the aircraft targets to reach and starts the straight-line autolanding process is In addition to the minimum required final height H S , which is normally 500 ft or 152.4 m for the manned aircraft [34], before landing at the predefined landing site safely, the aircraft must be heading straight toward the runway with the heading angleψ = θ S . Therefore, the targeted state for landing is [XȲZψ] T .
In summary, the control objective of the forced landing system in this article is to drive the aircraft into a specified target region in a certain range of orientation. The target region is characterized by a geometry composed of two parts: with the aforementioned site parameters. The parallelogram represents the target runway with the starting and ending locations S str and S end . The upper partial cylinder depicted with solid lines is the region of X cyld , and the lower dotted lines encircle the partial cone X cone (see Fig. 1). The first portion is for the case when the aircraft approaches the critical point P L with enough height, i.e., Z ≥Z . Within this region, the aircraft may take turns and adjust its orientation to consume extra height and energy until Z=Z. This extra height and energy consumption has been addressed in [8], which can be used after the aircraft reaches the target region, and it is out of the scope of this article. At this point, we only specify the mathematical expression of the cylinder region that is provided by where r ∈ R ≥0 and ε ∈ R ≥0 denote the threshold distance and angle, respectively, for the terminal constraints. The function β(v 1 , v 2 ) denotes the angle difference between two vectors v 1 and v 2 .
If the aircraft does not have enough energy to reach the height above the critical point P L , it can also travel directly to the area around the final straight-line path. This region forms a cone Once the aircraft enters the target region X f , it can execute the final landing step, which ensures a safe landing.

C. EMPC Formulation
If a single-engine aircraft with an engine failure has no thrust to propel itself, it is effectively a glider. Based on the glider dynamics, we consider the problem of integrating the path planning and dynamic control into a shrinking horizon EMPC optimization problem.
Let us denote the discretization sampling time of the aforementioned continuous dynamic model as T d which corresponds to the aircraft operation frequency. As a result, the corresponding discrete-time system dynamics is A UA effectively exchanges height for traveling distance. This is essential to follow a path with minimum height loss so that the aircraft has enough height to make turns in order to adjust its heading angle. Another important advantage of maintaining enough height is to preserve potential reachability to alternative landing sites in case of some unexpected emergencies. Therefore, the objective of the controller developed in this article is to maximize the terminal height Z N|t .
To complete the forced landing task, the target region X f in Section II-B is adopted as the terminal region of the EMPC algorithm. In selecting a landing site, a reachability analysis using the trochoidal turn path in [3] and [8] or path planning techniques in [35]- [37] to any candidate landing site is conducted. Based on that, an estimated time required to reach a candidate landing site can be calculated. According to the estimated time duration, the aircraft initiates its prediction horizon T = N · T d and starts the forced landing process from time instant t = 0. At each sampling time instant, a shrinking horizon EMPC at time t ∈ I [0,N−1] , formulated as follows, is solved: Note that, in this optimization problem, the number of decision variables u t:N|t = [u t|t , u t+1|t , . . . , u N|t ] are shrinking as the time instant t increases. After solving (12) at every time t, following the standard MPC algorithm, only the first input element u t|t , as a constant value, is applied to control the aircraft over the time interval T d . REMARK 1 In the standard hierarchical MPC algorithm, the traditional MPC controller is designed to track a path given by an RTO problem at the beginning of forced landing. To overcome the impact of uncertain wind along the path, many existing studies within the topic of wind estimate and compensation for small autonomous UA have been presented, e.g., disturbance observer [38], adaptive backstepping [39], vector field [40], etc. Since the primary objective of this article is to introduce the two-level EMPC framework, all of the formulations are deterministic, and combining the wind estimate techniques with the controller design is left for future research. REMARK 2 The optimization problem (12) of "pure" economic nature is solved every T d period and integrates the path planning and tracking problems. Constraint (12a) imposes the measured system state as the initial condition, and the future states are predicted according to constraint (12b). The pointwise-in-time constraint (12c) ensures the system feasibility within the prediction window. Moreover, the terminal constraint (12d) secures the aircraft to enter the desired terminal region X f at the end of the time horizon.
Although the EMPC controller generates a flyable trajectory between an original point (that includes the aircraft position, orientation and velocity) and a destination runway while minimizing the loss in altitude, such a single layer method is prohibitive for real-time applications due to the required long prediction horizon and the nonlinearity of the optimization problem. For the practical application of operating an aircraft, a discretized model with a high sampling frequency has to be chosen to approximate the continuoustime dynamic system. If the sampling time interval T d is selected to be sufficiently small, the MPC design based on the discretized model (11) can stabilize the original model (2). However, a smaller value of T d will require a faster online optimization. To address the implementation issues of optimization problem (12), we aim to propose a two-level framework in the next section.

III. TWO-LEVEL IMPLEMENTATION OF EMPC
This section presents an implementable two-level landing procedure that serves as a guideline for designing our autolanding algorithm. As shown in Fig. 2, the high-level EMPC planner adopts piecewise constant control actions to generate a feasible and optimal path for the aircraft at a lower sampling frequency. Then, a low-level time-varying MPC based on linearized models operates at a high sampling frequency to track the optimal path provided by the high-level EMPC planner.

A. High Level: Economic MPC Algorithm
We borrow the idea of piecewise constant MPC algorithm to yield a computationally efficient optimization problem [25]. A time interval T s is introduced as the integration step for the high-level EMPC. In the conventional discrete-time MPC, T s = T d , where T d is the integration step for the lower-level design. Note that T d is also the sampling time for the aircraft. However, in the piecewise constant MPC framework, the command inputs are held as constant values over several steps of T d time period. For example, as shown in Fig. 3, the holding time duration T s is calculated as T s = H · T d . In other words, for any prediction horizon N ∈ I ≥1 , the control command generated by the high-level EMPC remains constant over H low-level MPC integration steps of T d , i.e., u kH+i = u kH+ j , ∀k ∈ I [0,N −1] and ∀i, j ∈ Under the piecewise constant setup, the nonlinear discrete-time EMPC optimization problem subject to aircraft dynamics and various constraints at a high-level sampling time instant t = k 0 H can be formulated as REMARK 3 Note that the objective function (14) maximizes height only at the terminal state, which is intrinsically nonnegative and fulfills certain properties useful for stability proof argument. Specifically, since the optimal height is directly adopted in the construction of a Lyapunov function, there exists a lower bound for the actual height of the aircraft when it enters the terminal region as in Theorem 1.
The optimization problem (13)  The proposed high-level EMPC scheme generates a baseline control by: 1) providing a longer time interval for solving the optimization problem via reducing the sampling rate and 2) reducing the number of variables to be optimized. However, the EMPC strategy becomes an open-loop optimal control within the interval T s . Unfortunately, due to the mismatch between the mathematical model used in the high-level EMPC and the real aircraft system, the noises, and disturbances in the process, this kind of optimal control might not perform as designed. Within the interval T s , the MPC cannot suppress any tracking error, so the bandwidth associated with the MPC may not be adequate for stabilizing and controlling the UA that have fast dynamics.
To overcome these difficulties in the implementation of EMPC with online optimization, a low-level controller has been adopted in the proposed framework. The highlevel controller is the EMPC strategy described as before, which provides optimized state reference x r kH:(k+1)H−1 = x * kH:(k+1)H−1|kH and the corresponding baseline control u r kH:(k+1)H−1 = {u * kH|kH , . . . , u * kH|kH } after every time interval T s , whereas the low-level controller is a linear feedback controller paralleled to EMPC that can provide stability around the optimized state reference in the presence of disturbances and uncertainties. The high-level controller runs at a large sampling time T s due to the calculation time caused by solving nonlinear optimization problem (13). In contrast, the low-level controller works at a much higher sampling rate to reduce the impact of disturbances.
The low-level controller with sampling time interval T d is designed based on the linearized aircraft model (11) around the reference state x r kH:(k+1)H−1 and control input u r kH:(k+1)H−1 . The linearized system is expressed as , with δ t a parameter vector of the form δ t = [x r t , u r t ] T , and it is known at every step. Correspondingly, the state-input constraints will be shifted and depend on δ, so the time-varying constraints for the error system (15) are denoted bỹ Letting represent the closed set containing all the possible values of δ, a specific model matrix pair (A(δ), B(δ)) is denoted by λ ∈ , where the set is a multiplant description defined as Note that the linear time-varying system (15) is dependent on x r kH:(k+1)H−1 and u r kH:(k+1)H−1 , which are bounded in the high-level optimization. The parameter δ is bounded by , and each matrix pairs λ ∈ is time invariant and depends on a known parameter δ.
Assume that a full measurement or estimate of the state x t is available at the current time t ∈ I [kH:(k+1)H −1] . With a fixed prediction horizon N ∈ I ≥0 , we can formulate the following MPC problem: The objective function is to penalize the deviation of the predicted state and input sequences to the reference trajectories where (x t+N |t ) represents the terminal cost at time t + N given the current sampling time instant t. The cost function (18) is convex and the constraints (18a)-(18d) are linear. Therefore, it can be solved with efficient quadratic programming solvers. In the objective function, the matrices R and Q are positive definite to penalize the deviations from the reference input and state, respectively.
REMARK 4 According to the property of shrinking horizon EMPC, there are always solutions to optimization problem (13) in a certain environment if one solution can be found at the beginning. However, the low-level MPC may break down due to external disturbances. If this is the case, the prediction horizon will be reinitialized, and an updated shrinking horizon high-level EMPC will proceed accordingly. In worse situations when the external disturbances are quite large, the high-level EMPC may become infeasible as well, indicating that the selected target region may be no longer reachable. Hence, the landing site selection algorithm [41] will be triggered and the UA will fly to a different destination.
REMARK 5 Compared to an alternative scheme of the integration of high-level MPC and low-level proportionalintegral-derivative (PID) tracking controller, the proposed two-level MPC framework solves the high-level EMPC optimization problem by means of piecewise-constant inputs at a low frequency and solves the low-level quadratic problem with a short prediction horizon. Therefore, it significantly reduces the computational load. More importantly, this article addresses the forced landing problem where the UA loses its power and operates in an extreme situation. The satisfaction of the constraints becomes more relevant and important, compared with normal aircraft operation with power. Therefore, the low-level controller, in addition to providing robustness for the high-level EMPC controller in the presence of uncertainties and disturbances, must respect all the physical and operational constraints while maximizing the maneuver capability. This is why we adopted a lowlevel MPC controller to track the references generated by the high-level EMPC. Finally, the employment of low-level MPC (rather than a simple PID tracking controller) provides a well-structured formulation and specification of stability analysis given in Theorem 2.
Letũ * t:t+N −1|t = {ũ * t|t ,ũ * t+1|t ,ũ * t+2|t , . . . ,ũ * t+N −1|t } be the optimal solution of (18) at time t. Following the fashion of standard MPC algorithm, the sum of the first element u * t|t from the low-level MPC and the reference input u * t|t generated by the high-level EMPC is sent to the autopilot system to generate the real deflection control inputs, i.e., δ e , δ a , and δ r , as shown in Fig. 2. The proposed forced landing control scheme is then implemented to the six-degree-offreedom dynamic system described in [27] and [28]. In the next sampling time, the optimal control problem (18) is solved again using the new state measurements. When the model is linear time invariant, which corresponds to constant δ, a typical choice forX f is the maximal positive invariant set for the closed-loop system with state feedback control law u = Kx, where K is the associated linear-quadratic regulator (LQR) gain (i.e., the unconstrained infinite time optimal controller gain). However, since the model (15) is linear time varying, there are several maximal positive invariant setsX f (δ), and several different LQR feedback gain K (δ), one for each δ ∈ or equivalently λ ∈ . The maximal positive invariant setX f that is invariant for all λ ∈ is explored in [42].

IV. FEASIBILITY AND STABILITY ANALYSIS
The stability of the EMPC employing Lyapunov function was first proved in [43] under the assumption of strong duality. Such a condition was subsequently relaxed by dissipativity in [12] and [13], where Amrit et al. have shown how to transform the economic cost to a standard positive-definite function and to choose the transformed cost-to-go as a candidate Lyapunov function. In analogy, this section explores the stability of the closed-loop nominal sequence under EMPC control actions.
For the particularly defined objective function (14) in the high level, we are able to claim our first main result. THEOREM 1 Suppose the optimization problem (13) has a solution at the beginning. Then, the high-level closed-loop system in the absence of disturbances using the proposed EMPC algorithm is recursively feasible, and the target region is stable in the sense that the aircraft reaches it in finite time. Moreover, letting that the optimization problem (13) has a solution at t = 0 and the predicted terminal height is Z kH|0 , then the actual height Z f when the aircraft enters the region is no worse than Z N H|0 , viz., Z f ≥ Z N H|0 .

PROOF See Appendix A.
In the low-level tracking MPC with a fixed horizon, a tail term has to be added to the shifted input sequence to ensure recursive feasibility. The main idea is to find a local state feedback controller that stabilizes the system (15) inside an invariant terminal regionX f . Accordingly, some rigorous conditions of how to choose the terminal cost function and the terminal constraint setX f should be satisfied in order to ensure closed-loop stability and feasibility [23], [24]. In this respect, we make the following assumption: ASSUMPTION 1 (SEE[23, ASSUMPTION 2.33]) There exists a terminal setX f including the origin as an interior point, a terminal control input u f , and a terminal cost function (·), such that for all x ∈X f and δ ∈ , it holds that REMARK 6 Assumption 1 preserves the corresponding basic stability conditions of [23, Assumption 2.33] for linear time-varying systems. Following the general approach in the MPC algorithm, we have selected a quadratic terminal cost function (x) = x T Px, where P is the state weight, and the terminal control policy u f = K (δ)x. Then, Assumption 1(iii) can be written as (22) Notice that the above constraint implies that the quadratic function is a Lyapunov function, i.e., V (x) = x T Px; thus, the level sets of this Lyapunov function could be a choice of the terminal regionX f , and then, the other conditions in Assumption 1 are satisfied.
To simplify the calculation of terminal ingredients, one may consider a simple method that is to choose the following state equality as the terminal constraint (18d): x t+N |t = 0 (23) which corresponds toX f (δ) = {0}, K (δ) = 0 and (·) = 0 for all δ ∈ . This setup is applicable since the origin is an equilibrium point of the system (15) and trivially satisfies the terminal stability conditions in Assumption 1. However, the main drawback is also obvious that the constraints may cause the need of long horizon to obtain feasibility; otherwise, a control law with dead-beat behavior 3 is necessary if the horizon is short.
We are now in the position to state the main result for the low-level control setup. THEOREM 2 Supposing that Assumption 1 holds, the lowlevel closed-loop system using time-varying MPC algorithm (18) is recursively feasible and stable.
PROOF See Appendix B.
From previous analysis, we know that if there is no uncertainty and disturbance, the aircraft always follows the trajectory obtained from the high-level EMPC controller until it reaches the target region. However, in reality, disturbances (such as winds) may deviate the aircraft from the reference trajectory. For a bounded disturbance along the reference trajectory, the error dynamic is governed by the linearized system (15). Let K t denote the feedback gain, the control action isũ * t|t = K txt , and the resulting closed-loop system is expressed as Note that the system (24) can be converted to a polytopic system with the vertices determined by the known boundary values of uncertainties, and the associated robust stability can be ensured by using the parameter-dependent Lyapunov function techniques [44], [45]. Moreover, the high-level EMPC resets the error state to zero every T s time interval because it synchronizes the real aircraft state as the new initial state of the future reference trajectory. Thus, the use of low-level MPC guarantees that the divergence between the real system state and the reference trajectory caused by a bounded disturbance ω is bounded.

V. SIMULATION
Simulations were carried out to investigate the closedloop performance of the proposed two-level EMPC forced landing system. The controller was tested on a high-fidelity aircraft model of Aerosonde, the dynamics of which are detailed in [27] and [28]. First, the advantage of using the two-level EMPC controller compared with following a Dubins curve is presented. Then, the forced landing maneuver of Aerosonde has been conducted in both no-wind and wind conditions to test the robustness of the proposed algorithm.  In addition, the performance of the closed-loop system with the proposed two-level EMPC controller targeting at different predefined landing sites has also been investigated. The comparison of our approach with the traditional control method (i.e., single-level EMPC with autopilot) demonstrating the effectiveness of our method against disturbances (i.e., winds). Finally, the generality of this algorithm has been validated via another type of aircraft, i.e., Skywalker X8 fixed-wing UA.

A. Performance of the Two-Level EMPC Scheme
Consider forced landing for an Aerosonde UA with the relevant parameters in Table I Table II. The time setting for high-level shrinking horizon piecewise constant EMPC is N = 15 and H = 40. This setup implies that the high-level EMPC sampling time is T s = H · T d = 1 s and the initial high-level prediction time duration is N · T s = 15 s. To cover the same prediction time window, the high-level prediction horizon N (without considering piecewise constant control) has to be increased from 15 to 15 0.025 = 600 steps, which is extremely large. For the lowlevel linear time-varying MPC, the state and input weights are Q = I and R = 0.01I, where I is the identity matrix with proper dimensions. The low-level prediction horizon is selected as N = 20, which implies the corresponding time duration N · T d = 0.5 s. Moreover, the terminal equality constraint was adopted to avoid the specification of invariant sets and quadratic penalty functions. Our control algorithm was implemented using MAT-LAB R2020a and solved by fmincon and quadprog functions, on a computer with two-core 3.50-GHz Intel(R) Xeon(R) E5-1650 processor and 32-GB RAM. For the high-level shrinking horizon EMPC, the longest computational time could occur in the first optimization problem with the longest prediction horizon (i.e., the most number of decision variables), which is 0.9884 s. Note that the computational time is also affected by the number of active constraints. Among all the optimization problems being solved at high-level sampling time instants, the problem at t = 3 s has the largest number of active constraints. Specifically, 65 out of 385 constraints are active. In this case, the computational time is 0.8330 s. Additionally, the computational time for the optimization problem at t = 10 without active constraints is 0.0735 s. Therefore, we can conclude that the high-level EMPC optimization problem can be solved within the sampling time interval T s = 1 s. Regarding the low-level tracking MPC with a fixed prediction horizon, it is a quadratic programming problem taking around 0.0025 s to find a solution, which is significantly smaller than T d = 0.025 s. Moreover, further reduction in computation time for both levels can be easily obtained by increasing the corresponding sampling interval and utilizing more powerful machines.
Figs. 4 and 5 show the trajectories of the aircraft using the proposed two-level EMPC controller and tracking a Dubins curve. In Fig. 4, the critical point P L = [421.1210.56152.40], which is represented by the blue dot. This critical point is connected to the target region by a purple dashed line, which represents the straight-line autolanding process before the aircraft can land on the yellow runway safely. The top view of the target area X cyld is encircled by the orange dotted line. The red solid curve is obtained by maximizing the terminal height, while the khaki solid one is obtained by following the black dash-dotted Dubins curve. Fig. 4 shows that the aircraft can complete the task and reach the target region by both methods. However, as can be seen in Fig. 5, the resulting terminal aircraft height using the proposed two-level EMPC is significantly higher than that  following the Dubins curve. When the EMPC controller is used, the height loss is only about 10 m, whereas the Dubins method decreases the aircraft altitude by more than 150 m. Therefore, the EMPC method shows a great advantage in preserving aircraft height. Then, to test the performance of the proposed controller in the presence of wind uncertainties, simulations in a no-wind environment and a "1-cos" deterministic wind environment [46] have been compared. The time histories of the wind on X, Y , and Z dimensions are shown in Fig. 6, where the negative sign means that the wind was blowing along the negative direction of the axis. Simulation results are displayed in Fig. 7.
According to Fig. 7, the aircraft heights in the wind environment (represented by the red line) are similar to that in a no-wind environment (represented by the blue curve) within the initial 4 s due to the relatively small accumulative wind impact. However, after 4 s, the aircraft states are increasingly dissimilar. Particularly, the magnitude of wind in the Z-direction is always negative, which is equivalent to a dragging force to the aircraft, so the height of the aircraft in the wind environment is lower. The negative values of wind in X and Y directions hinder the aircraft from the target region, so it costs an extra second for the aircraft to  reach it in the wind condition. Moreover, it is observed that the altitude of the aircraft increased between 5 and 10 s due to the velocity loss, while the aircraft turned toward the direction of landing site. In particular, the lifted altitude is higher at the expense of decreasing more speed when there is no wind. The final speeds of these two environment scenarios are almost the same, which implies that the aircraft has to spend more energy to fly against the wind as the associated final height is lower. The time histories of the aircraft's angle of attack α and roll angle φ are shown in the bottom two subplots in Fig. 7, from which we can see that the two angles have similar patterns but different magnitudes when the aircraft was approaching to the target region because of the wind. The time histories of deflection control sequences for closed-loop implementation are shown in Fig. 8. Given the command angles from the two-level controller, an autopilot system is used to decide the deflection angles corresponding to the command. Note that the deflection angles of the aircraft are bounded by ±5 • . The elevator deflection almost keeps at a constant after 5 s, which corresponds to the variation of height and velocity with a constant rate. The aileron and rudder deflections are used to adjust the lateral angles toward the target region.

B. Different Landing Sites
Four different landing sites have been selected to test the reachability of the aircraft from a common initial state. The parameters of each landing site are shown in Table III, and the corresponding flight paths are shown in Fig. 9, showing that the aircraft is able to reach the corresponding target region for all landing sites. This generalization on the reachability to multiple landing sites is useful when the control algorithm breaks down (see more details in Section V-D) for one landing site since it can direct the aircraft to other available ones.

C. Different Tracking Controller
The low-level controller can improve the performance of the high-level EMPC in the presence of uncertainties and disturbances. To evaluate the benefits of the two-level EMPC control framework, a comparison has been conducted with single-level EMPC integrated with a simple tracking controller, e.g., PID-based autopilot. For these two cases, the same wind as shown in Fig. 6 has been implemented. The resulting aircraft's trajectories and deflection control inputs are shown in Figs. 10 and 11, respectively. As can be seen from the aircraft path view, the control method without low-level MPC cannot drive the aircraft into the target region. This is because autopilot in this case is tracking the piecewise constant command from high-level EMPC directly, rather than a command updated every T d s. Such comparison shows that adopting a low-level tracking MPC can update the command signal at a higher frequency Fig. 10. Aerosonde UA path view with the low-level MPC controller and with a PID controller. Red: two-level algorithm + autopilot system; and green: high-level EMPC + autopilot system. Fig. 11. Aerosonde UA control inputs with the low-level MPC controller and with a PID controller. Blue: two-level algorithm + autopilot system; and red: high-level EMPC + autopilot system. and provide better robustness against the uncertain wind.
Regarding the control inputs illustrated in Fig. 11, the two-level algorithm (represented by the blue line) is able to achieve better performance as the magnitudes of the deflection angles against the wind are smaller.

D. Different Aircraft Physical Parameters
Besides the Aerosonde model, the controller has also been implemented on a Skywalker X8 fixed-wing UA. Some basic parameters of this UA were collected from [47]. The X8 performed its forced landing starting from the origin to the same site as in Section V-A. Unlike Aerosonde, the X8 has a smaller size and the range of velocity is only between 10 and 20 m/s. The aircraft trajectories of these two types of UA under the zero wind condition are displayed in Fig. 12. Although the terminal height of X8 is lower, i.e., around 460 m, it can be controlled by the same two-level EMPC, demonstrating the wide applicability of the proposed control framework. Fig. 13 shows the airspeed and control surfaces of Aerosonde and X8. Note that the speed of X8 is  almost constant at 15 m/s, and the angular deflections have more input oscillations compared to those of Aerosonde. This is because the X8 took turns at the beginning and near the end of its flight by adjusting the deflection angles more frequently to enter the terminal region with the desired heading angle while maximizing the height. Moreover, due to the relatively low velocity, the X8 has to spend more time to reach the target region.

VI. CONCLUSION
In this article, a novel two-level EMPC scheme is developed for fixed-wing UA systems in order to accomplish a successful emergency landing when the aircraft loses power. The proposed scheme integrates the high-level EMPC method, which aims to minimize the height loss, and the low-level time-varying tracking MPC, so that the closed-loop feasibility and stability are guaranteed. The effectiveness of the proposed approach has been demonstrated by several case studies. In particular, the robustness and generality of the proposed controller on the aircraft flight have been tested and analyzed. The benefits of using a two-level MPC framework, rather than the traditional high-level MPC planner with low-level simple tracking controller, have also been investigated. In the future, we will combine the control algorithm with wind estimation techniques and investigate how to take advantage of favorable winds.