Sage Journals: Discover world-class research

Abstract

A lifting body unmanned aerial vehicle (UAV) generates lift by its body and shows many significant advantages due to the particular shape, such as huge loading space, small wetted area, high-strength fuselage structure, and large lifting area. However, designing the control law for a lifting body UAV is quite challenging because it has strong nonlinearity and coupling, and usually lacks it rudders. In this paper, an explicit nonlinear model predictive control (ENMPC) strategy is employed to design a control law for a saucer-shaped UAV which can be adequately modeled with a rigid 6-degrees-of-freedom (DOF) representation. In the ENMPC, control signal is calculated by approximation of the tracking error in the receding horizon by its Taylor-series expansion to any specified order. It enhances the advantages of the nonlinear model predictive control and eliminates the time-consuming online optimization. The simulation results show that ENMPC is a propriety strategy for controlling lifting body UAVs and can compensate the insufficient control surface area.

1. Introduction

A lifting body is basically an aerial vehicle that chiefly or solely generates lift by its body [1]. Unlike traditional UAVs which obtain lift mainly by the wings, lifting bodies generate little or small amount of lift by wings and some of them even have no wings. The first lifting body XF5U was designed by the American engineer Charles H. Zimmerman and showed advantageous aerodynamic characteristics [2]. Then, national aeronautics and space administration (NASA) gave more superiority in loading space, low-speed controllability, structure strength, and so forth, by studying M2-F1, M2-F2, M2-F3, HL-10, HL-20, X-24A, X-24B, and many other kinds of lifting bodies [3]. But these design schemes weren't taken seriously for many reasons at that time. From a technical point of view, the dominant issue has been stability and control, which to this day continues to plague this class of vehicles [1]. With the wide application of UAVs [4], people draw more attention to the lifting bodies which can overcome some shortcomings of conventional ones and meet more aviation demands. Therefore, many practical challenges are waiting to be solved for the lifting body. One of the most important aspects is the control system.

There are many kinds of control strategies, such as fuzzy sliding mode control, adaptive control, robust control, and optimal control, that can be used for nonlinear dynamic system [5]. Nowadays, model predictive control (MPC) has been regarded as a prospective method for highly nonlinear dynamic systems such as UAVs [6]. MPC strategy uses a dynamic model to predict the future behavior of a plant and subsequently an optimal control decision can be made based on the estimated future state and the desired future state [7]. The linear model predictive control is used by modeling a discrete linear system [8]. To develop MPC for nonlinear system, many approaches have been proposed, such as the State-Dependent Riccati Equation Technique [9], using second-order model approximation [10], and solving the Hamilton-Jacobi-Bellman equation by power series expansion [11]. However, the methods are mentioned above require a plenty of online computations which is often computationally complex and time consuming and the real-time NMPC implementation is usually limited to slow processes where the sampling time is sufficient to support the computational needs [12]. Chen developed a nonlinear model predictive control (NMPC) using approximation. This method can stabilize the original nonlinear systems by properly choosing predictive times and the order for Taylor-series expansion [13]. It is assumed that control weighting is not in the performance index and that control order is the difference between the Taylor-series expansion and the relative degree of the system [14]. According to this method, Liu et al. introduced an explicit NMPC, which inherits the advantages of the NMPC but eliminates the time-consuming online optimization [15]. The control signals in this method can be calculated instantaneously to respond to the fast dynamic systems and effectively restrain the disturbances. Besides, by removing the online optimization process, the onboard hardware can be simplified and the software development is accelerated.

The controlled object we used in this paper is a saucer-shaped lifting body UAV. The body of this UAV is designed like a saucer to provide a large number of internal spaces and improve lifting area. According to the special shape of this UAV, it is easy to conclude that the lateral maneuverability is not as good as a conventional aerial vehicle. Consequently, a better control system is needed to avoid the shortcomings of this UAV including other kinds of lifting body aerial vehicles. Therefore, we design a control law using the ENMPC strategy mentioned in this paper. The results are compared with the NMPC algorithm and conventional PID control law. Besides we also induce some wind disturbances to investigate the robustness of this method which is not discussed in any references. Due to the advantages of ENMPC, the control framework has a low complexity and high reliability. Therefore, the UAV has a nice response effect. This paper expands the application range of this method and gives other lifting bodies a feasible control strategy. The wind disturbances simulation also proves that this method has good robustness.

The remainder of this paper is organized as follows. Section 2 introduces the mathematical model of the saucer-shaped lifting body UAV. Section 3 describes the algorithm of ENMPC and designs the control law for the UAV. Then, in Section 4, series of simulations are conducted and the comparison results are discussed in details. Our concluding remarks and future work are presented in the final section.

2. The Saucer-Shaped Lifting Body UAV Model

The saucer-shaped UAV used in this paper is showed in Figure 1. The UAV is composed by a saucer-shaped fuselage, two small wings, twin vertical tail with pelvic fin, and two tailerons. The fuselage provides most of lift and huge storage space for the UAV. The area of two small wings only accounts for 12.2 percent of the total area. They are mainly used to reduce the induced drag instead of providing lift. The twin vertical tail with pelvic fin can help to maintain the lateral stability. And the tailerons mounted at the tail can be used to control the longitudinal maneuvering flight by deflecting in the same direction and control the lateral maneuver by differential deflection. Compared with flying wings, this UAV has a very small aspect ratio which makes it enjoys a stronger structural strength. The X-30 and X-33 made by NASA also adopt the similar aerodynamic layout [16, 17].

Figure 1:

The saucer-shaped lifting body UAV under the wind tunnel test.

This UAV can be modeled as a rigid body, and undergoes three-dimensional motion described by three inertial position coordinates and three Euler angles. Hence, 12 state variables are required to describe motion of the UAV at a given instance in time [18]. In the equations below, the ground frame is assumed to be a satisfactory reference frame. Besides, the earth is assumed to be flat and wind is stationary. Kinematic translational and differential equations of motion for the general system are provided by

\begin{matrix} {\begin{Bmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \end{Bmatrix}} = [T_{B / I}] {\begin{Bmatrix} u \\ v \\ w \end{Bmatrix}}, \\ {\begin{Bmatrix} \dot{ψ} \\ \dot{θ} \\ \dot{φ} \end{Bmatrix}} = [\begin{bmatrix} 0 & s_{φ} / c_{θ} & c_{φ} / c_{θ} \\ 0 & c_{φ} & - s_{φ} \\ 1 & t_{θ} s_{φ} & t_{θ} c_{φ} \end{bmatrix}] {\begin{Bmatrix} p \\ q \\ r \end{Bmatrix}} . \end{matrix}

(1)

The common shorthand notation for trigonometric functions is used throughout the paper (sin (a) = s_a, cos (a) = c_a, tan (a) = t_a). The matrix $T_{B / I}$ represents the transformation matrix from the body coordinate frame to the inertial coordinate frame:

T_{B / I} = [\begin{bmatrix} c_{θ} c_{ψ} & s_{θ} c_{ψ} s_{φ} - s_{ψ} c_{φ} & s_{θ} c_{ψ} c_{φ} + s_{ψ} s_{φ} \\ c_{θ} s_{ψ} & s_{θ} s_{ψ} s_{φ} + c_{ψ} c_{φ} & s_{θ} s_{ψ} c_{φ} - c_{ψ} s_{φ} \\ - s_{θ} & c_{θ} s_{φ} & c_{θ} c_{φ} \end{bmatrix}] .

(2)

The force and moment differential equations are represented in

\begin{matrix} {\begin{Bmatrix} \dot{u} \\ \dot{v} \\ \dot{w} \end{Bmatrix}} = g {\begin{Bmatrix} - s_{φ} \\ s_{φ} c_{θ} \\ c_{φ} c_{θ} \end{Bmatrix}} - [\begin{bmatrix} 0 & - r & q \\ r & 0 & - p \\ - q & p & 0 \end{bmatrix}] {\begin{Bmatrix} u \\ v \\ w \end{Bmatrix}} + \frac{1}{m} {\begin{Bmatrix} F_{x} \\ F_{y} \\ F_{z} \end{Bmatrix}}, \\ {\begin{Bmatrix} \dot{p} \\ \dot{q} \\ \dot{r} \end{Bmatrix}} = I_{T}^{- 1} ({\begin{Bmatrix} \bar{L} \\ M \\ N \end{Bmatrix}} + [\begin{bmatrix} 0 & r & - q \\ - r & 0 & p \\ q & - p & 0 \end{bmatrix}] I_{T} {\begin{Bmatrix} p \\ q \\ r \end{Bmatrix}}), \end{matrix}

(3)

where

\begin{matrix} {\begin{Bmatrix} F_{x} \\ F_{y} \\ F_{z} \end{Bmatrix}} = {\begin{Bmatrix} δ_{T} \\ 0 \\ 0 \end{Bmatrix}} + T_{B / W} {\begin{Bmatrix} - D \\ Y \\ - L \end{Bmatrix}}, \\ T_{B / W} = [\begin{bmatrix} c_{α} c_{β} & - c_{α} s_{β} & - s_{α} \\ s_{β} & c_{β} & 0 \\ s_{α} c_{β} & - s_{α} s_{β} & c_{α} \end{bmatrix}], \\ I_{T} = [\begin{bmatrix} I_{X X} & 0 & I_{X Z} \\ 0 & I_{Y Y} & 0 \\ I_{X Z} & 0 & I_{Z Z} \end{bmatrix}], I_{T}^{- 1} = [\begin{bmatrix} I_{X X I} & 0 & I_{X Z I} \\ 0 & I_{Y Y I} & 0 \\ I_{X Z I} & 0 & I_{Z Z I} \end{bmatrix}] . \end{matrix}

(4)

The matrix $T_{B / W}$ refers to the transformation matrix from the body coordinate frame to the wind coordinate frame. The effects of the control inputs and the wing flap deflections both influence the aerodynamic forces acting on the UAV. From the equations and parameters listed above, the control inputs of this UAV only include elevator, aileron, and throttle.

The aerodynamic forces are expressed as:

{\begin{Bmatrix} D \\ Y \\ L \end{Bmatrix}} = (\frac{1}{2}) ρ S V^{2} {\begin{Bmatrix} C_{D 0} + C_{D α} α \\ C_{Y β} β \\ C_{L 0} + C_{L α} α + C_{L δ_{e}} δ e \end{Bmatrix}} .

(5)

Similarly, the aerodynamic moments are described as follows:

{\begin{Bmatrix} \bar{L} \\ M \\ N \end{Bmatrix}} = (\frac{1}{2}) ρ S V^{2} {\begin{Bmatrix} (\frac{b}{2}) (C_{l β} β + C_{l p} p + C_{l r} r + C_{l δ_{a}} δ_{a}) \\ c (C_{m 0} + C_{m α} α + C_{m q} q + C_{m δ_{e}} δ_{e}) \\ (\frac{b}{2}) (C_{n β} β + C_{n p} p + C_{n r} r + C_{n δ_{a}} δ_{a}) \end{Bmatrix}},

(6)

where the velocity and aerodynamic angles are defined as

\begin{matrix} V = {(u^{2} + v^{2} + w^{2})}^{1 / 2}, \\ α = \arctan (\frac{w}{u}), \\ β = \arcsin (\frac{v}{V}) . \end{matrix}

(7)

The UAV model uses two control surfaces elevator and aileron. The vertical tails are fixed on the UAV and cannot be used as control surface.

3. Explicit Nonlinear Model Predictive Control Law for the Saucer-Shaped UAV

3.1. Explicit Nonlinear Model Predictive Control

For simplicity, a nonlinear MIMO system can be described by

\begin{matrix} \dot{x} = f (x, u), \\ y = h (x), \end{matrix}

(8)

where x ∈ ℝℝ^m, u ∈ ℝⁿ, y ∈ ℝⁿ are state, input and output, respectively. For a nonlinear MIMO system, it is generally known that each control will appear after a specific number of derivatives of the output [19]. The relative degree ρ_i is the number of time derivatives of the ith output needed for control to appear. If continuously differentiating the output after the control input appears, time derivatives of control will appear and the number of the input derivatives r is defined as control order. The ith output of the system in the receding horizon can be approximated by its Taylor-series expansion up to order ρ_i + r:

\begin{matrix} {\hat{y}}_{i} (t + τ) \approx y_{i} (t) + τ {\dot{y}}_{i} (t) + \dots + \frac{τ^{r + ρ_{i}}}{(r + ρ_{i})!} {y_{i}}^{[r + ρ_{i}]} (t) = {\bar{T}}_{i} Y_{i} (t), \\ {\bar{T}}_{i} \equiv [\begin{bmatrix} 1 & τ & \dots & \frac{τ^{r + ρ_{i}}}{(r + ρ_{i})!} \end{bmatrix}], \\ Y_{i} (t) \equiv {[\begin{bmatrix} y_{i} (t) & {\dot{y}}_{i} (t) & \dots & {y_{i}}^{[r + ρ_{i}]} (t) \end{bmatrix}]}^{T}, \end{matrix}

(9)

where i = 1,2, …, n. For each channel in the output matrix, the control orders r are the same for each vector in the output matrix and can be decided during the control design, but the relative degrees ρ_i are different and determined by the specific model structure. The approximation of the entire outputs can be described in a matrix form as follows:

\hat{y} (t + τ) = [\begin{bmatrix} {\hat{y}}_{1} (t + τ) \\ {\hat{y}}_{2} (t + τ) \\ \dots \\ {\hat{y}}_{n} (t + τ) \end{bmatrix}] = [\begin{bmatrix} {\bar{T}}_{1} & 0 \\ {\bar{T}}_{2} \\ ⋱ \\ 0 & {\bar{T}}_{n} \end{bmatrix}] [\begin{bmatrix} Y_{1} (t) \\ Y_{2} (t) \\ ⋮ \\ Y_{n} (t) \end{bmatrix}] = \bar{T} \bar{Y} (t) .

(10)

To simplify the calculation, we change the order of the variables in (10) and separate the terms without control input from the remaining terms. Therefore, $\bar{T}$ and $\bar{Y}$ can be expressed as

\begin{matrix} \bar{T} = [\begin{bmatrix} T_{f} & T_{s} \end{bmatrix}] = [\begin{matrix} \begin{matrix} {\bar{τ}}_{1} & 0_{1 \times ρ_{2}} & \dots & 0_{1 \times ρ_{n}} \\ 0_{1 \times ρ_{1}} & {\bar{τ}}_{2} & \dots & 0_{1 \times ρ_{n}} \\ \dots & \dots & \dots & \dots \\ 0_{1 \times ρ_{1}} & 0_{1 \times ρ_{2}} & \dots & {\bar{τ}}_{n} \end{matrix} & \begin{matrix} {\tilde{τ}}_{1} & {\tilde{τ}}_{2} & \dots & {\tilde{τ}}_{r + 1} \end{matrix} \end{matrix}], \\ \bar{Y} (t) = {[\begin{matrix} \begin{matrix} {\bar{Y}}_{1} {(t)}^{T} & {\bar{Y}}_{2} {(t)}^{T} & \dots & {\bar{Y}}_{n} {(t)}^{T} \end{matrix} & \begin{matrix} {\tilde{Y}}_{1} {(t)}^{T} & {\tilde{Y}}_{2} {(t)}^{T} & \dots & {\tilde{Y}}_{r + 1} {(t)}^{T} \end{matrix} \end{matrix}]}^{T}, \end{matrix}

(11)

where

\begin{matrix} {\bar{Y}}_{i} (t) = {[\begin{bmatrix} y_{i} (t) & {\dot{y}}_{i} (t) & \dots & {y_{i}}^{[ρ_{i} - 1]} (t) \end{bmatrix}]}^{T}, i = 1,2, \dots, n, \\ {\tilde{Y}}_{i} (t) = {[\begin{bmatrix} {y_{1}}^{[ρ_{1} + i - 1]} (t) & {y_{2}}^{[ρ_{2} + i - 1]} (t) & \dots & {y_{n}}^{[ρ_{n} + i - 1]} (t) \end{bmatrix}]}^{T}, i = 1,2, \dots, r + 1, \\ {\bar{τ}}_{i} = [\begin{bmatrix} 1 & τ & \dots & \frac{τ^{ρ_{i} - 1}}{(ρ_{i} - 1)!} \end{bmatrix}], i = 1,2, \dots, n, \\ {\tilde{τ}}_{i} = diag {\begin{Bmatrix} \frac{τ^{ρ_{1} + i - 1}}{(ρ_{1} + i - 1)!} & \frac{τ^{ρ_{2} + i - 1}}{(ρ_{2} + i - 1)!} & \dots & \frac{τ^{ρ_{n} + i - 1}}{(ρ_{n} + i - 1)!} \end{Bmatrix}}, i = 1,2, \dots, r + 1, \\ T_{f} = [\begin{bmatrix} {\bar{τ}}_{1} & 0_{1 \times ρ_{2}} & \dots & 0_{1 \times ρ_{n}} \\ 0_{1 \times ρ_{1}} & {\bar{τ}}_{2} & \dots & 0_{1 \times ρ_{n}} \\ \dots & \dots & \dots & \dots \\ 0_{1 \times ρ_{1}} & 0_{1 \times ρ_{2}} & \dots & {\bar{τ}}_{n} \end{bmatrix}], \\ T_{s} = [\begin{bmatrix} {\tilde{τ}}_{1} & {\tilde{τ}}_{2} & \dots & {\tilde{τ}}_{r + 1} \end{bmatrix}] . \end{matrix}

(12)

According to the equations shown above, we now can construct matrix ${\tilde{Y}}_{i}$ , i = 1,2, …, r + 1:

{\tilde{Y}}_{1} (t) = [\begin{bmatrix} y_{1}^{[ρ_{1}]} (t) \\ y_{2}^{[ρ_{2}]} (t) \\ \dots \\ y_{n}^{[ρ_{n}]} (t) \end{bmatrix}] = [\begin{bmatrix} L_{f}^{ρ_{1}} h_{1} (x) \\ L_{f}^{ρ_{2}} h_{2} (x) \\ \dots \\ L_{f}^{ρ_{n}} h_{n} (x) \end{bmatrix}] + A (x) u,

(13)

where the nonlinear terms L_f^ρ_ih_i(x), i = 1,2, …, n are parts of ρ_ith time derivatives of the ith output without control input. By differentiating (13) with respect to time, we can get

{\tilde{Y}}_{2} (t) = [\begin{bmatrix} y_{1}^{[ρ_{1} + 1]} (t) \\ y_{2}^{[ρ_{2} + 1]} (t) \\ \dots \\ y_{n}^{[ρ_{n} + 1]} (t) \end{bmatrix}] = [\begin{bmatrix} L_{f}^{ρ_{1} + 1} h_{1} (x) \\ L_{f}^{ρ_{2} + 1} h_{2} (x) \\ \dots \\ L_{f}^{ρ_{n} + 1} h_{n} (x) \end{bmatrix}] + A (x) \dot{u} + p_{1} (x, u),

(14)

and p₁(x, u) is a nonlinear vector function of state and input. The higher derivatives of the output and ${\tilde{Y}}_{i}$ , i = 1,2, …, r + 1 can be calculated by repeating this procedure as follows:

{\tilde{Y}}_{r + 1} (t) = [\begin{bmatrix} y_{1}^{[ρ_{1} + r]} (t) \\ y_{2}^{[ρ_{2} + r]} (t) \\ \dots \\ y_{n}^{[ρ_{n} + r]} (t) \end{bmatrix}] = [\begin{bmatrix} L_{f}^{ρ_{1} + r} h_{1} (x) \\ L_{f}^{ρ_{2} + r} h_{2} (x) \\ \dots \\ L_{f}^{ρ_{n} + r} h_{n} (x) \end{bmatrix}] + A (x) u^{[r]} + p_{r} (x, u, \dot{u}, \dots, u^{[r]}) .

(15)

In the same as (10), the desired output signal w(t) can also be appropriated as follows:

w (t + τ) = [\begin{bmatrix} w_{1} (t + τ) \\ w_{2} (t + τ) \\ \dots \\ w_{n} (t + τ) \end{bmatrix}] = [\begin{bmatrix} T_{f} & T_{s} \end{bmatrix}] \bar{W} (t) .

(16)

The construction of $\bar{W} (t)$ can consult the structure of $\bar{Y} (t)$ .

Expanding the input in a similar way we have

\tilde{u} (t + τ) \approx \tilde{u} (t) + τ {\tilde{u}}^{[1]} (t) + \dots + \frac{τ^{r}}{r!} {\tilde{u}}^{[r]} (t) .

(17)

Thereby, the outputs can be decided by the control variables $\bar{u} = {\tilde{u}, {\tilde{u}}^{[1]}, \dots, {\tilde{u}}^{[r]}}$ .

To any control problems, we hope that the output y(t) follows the prescribed reference w(t) well. Therefore, the control law can be achieved by minimizing a receding horizon performance index:

J = \frac{1}{2} \int_{0}^{T} {(\hat{y} (t + τ) - w (t + τ))}^{T} Q (\hat{y} (t + τ) - w (t + τ)) d τ,

(18)

where the weighing matrix $Q = diag {\begin{Bmatrix} q_{1} & q_{2} & \dots & q_{n} \end{Bmatrix}}$ , q_i > 0, i = 1,2, …, n. Recalling the output and reference approximation equations (10) and (16), J can also be expressed as

J = \frac{1}{2} {(\bar{Y} (t) - \bar{W} (t))}^{T} [\begin{bmatrix} Γ_{1} & Γ_{2} \\ Γ_{2}^{T} & Γ_{3} \end{bmatrix}] (\bar{Y} (t) - \bar{W} (t)),

(19)

\begin{matrix} Γ_{1} = \int_{0}^{T} T_{f}^{T} Q T_{f} d τ, \\ Γ_{2} = \int_{0}^{T} T_{f}^{T} Q T_{s} d τ, \\ Γ_{3} = \int_{0}^{T} T_{s}^{T} Q T_{s} d τ . \end{matrix}

(20)

Instead of minimizing the performance index in (18) with respect to control profile u directly, we can minimize the approximated index J in (19) with respect to $\bar{u}$ ; the necessary condition for optimality is

\frac{\partial J}{\partial \bar{u}} = 0 .

(21)

The optimal control variables ${\bar{u}}^{*}$ are given by solving the nonlinear equation (21) and can be used to construct the optimal control profile defined by (17). The input u is limited in a reasonable range according to different controlled objects. As in MPC only the current control in the control profile is implemented; the explicit solution is ${\tilde{u}}^{*} = \tilde{u} (t + τ)$ , for τ = 0. The resulting controller is given by

{\tilde{u}}^{*} = - A {(x)}^{- 1} (K M_{ρ} + M_{1}),

(22)

where $K \in ℝ^{n \times (ρ_{1} + ρ_{2} + \dots + ρ_{n})}$ is the first n row of the matrix $Γ_{3}^{- 1} Γ_{2}^{T} \in ℝ^{n (r + 1) \times (ρ_{1} + ρ_{2} + \dots + ρ_{n})}$ . The matrix $M_{ρ} \in ℝ^{(ρ_{1} + ρ_{2} + \dots + ρ_{n})}$ and matrix M_i ∈ ℝⁿ are defined as

\begin{matrix} M_{ρ} = [\begin{bmatrix} {\bar{Y}}_{1} (t) \\ {\bar{Y}}_{2} (t) \\ \dots \\ {\bar{Y}}_{n} (t) \end{bmatrix}] - [\begin{bmatrix} {\bar{W}}_{1} (t) \\ {\bar{W}}_{2} (t) \\ \dots \\ {\bar{W}}_{n} (t) \end{bmatrix}], \\ M_{i} = [\begin{bmatrix} L_{f}^{ρ_{1} + i - 1} h_{1} (x) \\ L_{f}^{ρ_{2} + i - 1} h_{2} (x) \\ \dots \\ L_{f}^{ρ_{n} + i - 1} h_{n} (x) \end{bmatrix}] - {\tilde{W}}_{i} {(t)}^{T}, i = 1,2, \dots, r + 1 . \end{matrix}

(23)

In order to obtain the control signals, conventional MPC algorithm needs to solve an optimization problem at each sampling instant. The ENMPC method can effectively avoid the computationally intensive online optimization and shows better control effects [20].

3.2. Saucer-Shaped UAV Control Law

To control an aerial vehicle by this ENMPC method, we choose Euler angels as outputs. From the special shape of this UAV, we can conclude that the yaw dirigibility is not as good as conventional ones. Besides, with the absence of a rudder, the yaw motion can only be achieved by manoeuvring the roll direction. Since in this method the number of outputs needs to be equal to the number of control inputs, we determine $y = {[\begin{bmatrix} h & θ & φ \end{bmatrix}]}^{T}$ and $u = {[\begin{bmatrix} δ_{T} & δ_{e} & δ_{a} \end{bmatrix}]}^{T}$ .

After the first derivatives, we can obtain $\dot{y}$ from (1):

\dot{y} = [\begin{bmatrix} {\dot{y}}_{1} \\ {\dot{y}}_{2} \\ {\dot{y}}_{3} \end{bmatrix}] = [\begin{bmatrix} \dot{h} \\ \dot{θ} \\ \dot{φ} \end{bmatrix}] = [\begin{bmatrix} u s_{θ} - v s_{φ} c_{θ} - w c_{φ} c_{θ} \\ q c_{φ} - r s_{φ} \\ p + (r c_{φ} + q s_{φ}) t_{θ} \end{bmatrix}] .

(24)

The ENMPC method requires the output equations to be differentiated until the control inputs appear. From the model of the UAV, we perceive that the control inputs appear in $\dot{u}$ , $\dot{v}$ , $\dot{w}$ , $\dot{p}$ , $\dot{q}$ , and $\dot{r}$ . Then, all the control inputs appear after the second time derivative and

\ddot{y} = [\begin{bmatrix} {\ddot{y}}_{1} \\ {\ddot{y}}_{2} \\ {\ddot{y}}_{3} \end{bmatrix}] = [\begin{bmatrix} \ddot{h} \\ \ddot{θ} \\ \ddot{φ} \end{bmatrix}] = [\begin{bmatrix} \dot{u} s_{θ} + u c_{θ} \dot{θ} - \dot{v} s_{φ} c_{θ} - v c_{φ} c_{θ} \dot{φ} + v s_{φ} s_{θ} \dot{θ} - \dot{w} c_{φ} c_{θ} + w s_{φ} c_{θ} \dot{φ} + w c_{φ} s_{θ} \dot{θ} \\ \dot{q} c_{φ} - q s_{φ} \dot{φ} - \dot{r} s_{φ} - r c_{φ} \dot{φ} \\ \dot{p} + \frac{(r c_{φ} + q s_{φ}) \dot{θ}}{c_{θ}^{2}} + t_{θ} (\dot{r} c_{φ} - r s_{φ} \dot{φ} + \dot{q} s_{φ} + q c_{φ} \dot{φ}) \end{bmatrix}] .

(25)

Therefore, the relative degree for each input is 2, and the matrix A(x) in (13) can be described as

A (x) = [\begin{bmatrix} \frac{s_{θ}}{m} & \frac{P V^{2} C_{L δ_{e}} (c_{α} c_{φ} c_{θ} + s_{α} s_{θ})}{m} & 0 \\ 0 & \frac{P V^{2} c C_{m δ_{e}} c_{φ}}{I_{Y Y}} & \frac{- P V^{2} (b / 2) (I_{X Z} C_{l δ_{a}} + I_{X X} C_{n δ_{a}}) s_{φ}}{I_{X X} I_{Z Z} - I_{X Z}^{2}} \\ 0 & \frac{P V^{2} c C_{m δ_{e}} t_{θ} s_{φ}}{I_{Y Y}} & \frac{P V^{2} (b / 2) ((I_{X Z} c_{φ} t_{θ} + I_{Z Z}) C_{l δ_{a}} + (I_{X X} c_{φ} t_{θ} + I_{X Z}) C_{n δ_{a}})}{I_{X X} I_{Z Z} - I_{X Z}^{2}} \end{bmatrix}] .

(26)

Therefore, when the desired fight attitude of the UAV is determined, we can get the desired output $y^{*} = {[\begin{bmatrix} h^{*} & θ^{*} & φ^{*} \end{bmatrix}]}^{T}$ . According to (23),

\begin{matrix} M_{ρ} = {[\begin{bmatrix} h - h^{*} & \dot{h} - {\dot{h}}^{*} & θ - θ^{*} & \dot{θ} - {\dot{θ}}^{*} & φ - φ^{*} & \dot{φ} - {\dot{φ}}^{*} \end{bmatrix}]}^{T}, \\ M_{1} = [\begin{bmatrix} \ddot{h} \\ \ddot{θ} \\ \ddot{φ} \end{bmatrix}] - A (x) [\begin{bmatrix} T \\ δ_{e} \\ δ_{a} \end{bmatrix}] - [\begin{bmatrix} {\ddot{h}}^{*} \\ {\ddot{θ}}^{*} \\ {\ddot{φ}}^{*} \end{bmatrix}] . \end{matrix}

(27)

Thus, the resulting controller can be calculated by (22). The control structure can be shown in Figure 2.

Figure 2:

ENMPC structure for the saucer-shaped UAV.

4. Simulations Results

The main purpose of this simulation is to prove the practicability and check the performance of the ENMPC strategy on a lifting-body UAV. There are two simulations conducted in this section. The first one aims at verifying the control effect in both longitudinal and lateral directions. The other one is conducted under longitudinal wind disturbance to evaluate the robustness of the control law. The model for the simulations is the one introduced in Section 2. For comparison, a PID controller is employed. In the PID method, we also use the model given in Section 2. Then, the nonlinear equation of motion is linearized on the condition of little perturbation around the equilibrium point and divided into longitudinal and lateral equations. The control law is obtained according to the linear equation and applied to the UAV [21]. The preceding system of equations describing the UAV is numerically integrated using a fourth-order Runge-Kutta algorithm. In the following results, the control derivative is updated every 0.005 s. The thrust provided by the engine used on the UAV is no more than 500 N. The elevator and aileron are limited between −20 deg to 20 deg. In order to be more close to the practical situation, the rate of change for each control surface is limited to −20 deg/s to 20 deg/s. Tables 1, 2, and 3 list characteristics of the UAV used for the simulation, including dimensions, mass properties, and aerodynamic derivatives. Stability derivative values were taken by a series of static and dynamic derivative wind tunnel tests.

Table 1:

Dimensions and mass properties of the saucer-shaped UAV.

Parameter		Parameter
Mass (kg)	175	I_XX (kgm²)	169.92
Wing area (m²)	5.8861	I_YY (kgm²)	128.799
Fuselage length (m)	2.5	I_ZZ (kgm²)	293.258
Wingspan (m)	3.6	I_XZ (kgm²)	−0.37
Mean aerodynamic chord (m)	2.4

Table 2:

Longitudinal aerodynamic derivatives of the saucer-shaped UAV.

Lift (L)		Drag (D)		Pitching moment (M)

C _L0	0.063	C _D0	0.007	C _m0	−0.0008
C _Lα	2.51	C _Dα	0.18	C _mα	−0.09
C _{Lδ_e}	0.0061			C _mq	−3.2
				C _{mδ_e}	−0.0026

Table 3:

Lateral-directional aerodynamic derivatives of the saucer-shaped UAV.

Side force (Y)		Rolling moment ( $\bar{L}$ )		Yawing moment (N)

C _Yβ	−0.306	C _lβ	−0.08	C _nβ	0.1008
		C _lp	−0.14	C _np	−0.0491
		C _lr	0.079	C _nr	−0.093
		C _{lδ_a}	−0.00082	C _{nδ_a}	−0.00008

Figure 3 shows the 3D view of the control results by different control strategies. In the PID method, the UAV's longitudinal and lateral equations are decoupled around an equilibrium point, and the control law is given by respective linearization. In these techniques, the linearization is restricted around the equilibrium point, and mainly the stability property is considered and little insight into other aspects of dynamic behaviour is provided [22]. But the ENMPC strategy can give a strong relationship between each variable. Therefore, it can keep the height stable even when the lateral maneuver appeared, but a descent will arise when using PID control law, which is obvious in Figure 4. That means the UAV track is smoother when using ENMPC compared with PID. The outputs are approximated with a sixth-order Taylor-series expansion, and the prediction horizon is 0–6.0 s in the ENMPC strategy.

Figure 3:

3D view of flight.

Figure 4:

The saucer-shaped UAV height response.

Figures 5 –7 show simulation results of the control signals corresponding to the trajectory in Figure 3.

Figure 5:

Throttle behavior.

Figure 6:

Elevator deflection.

Figure 7:

Aileron deflection.

In PID method, throttle signal is obtained directly by longitudinal variables. However, in ENMPC strategy, it is determined by both longitudinal variables and lateral variables. Therefore, when the UAV exhibits its later maneuverability, the throttle signal is influenced by the changing of lateral variables and adjusted immediately to keep a stable flight path. In Figure 5, the changing of the throttle is more violent under the control of ENMPC strategy, which can keep other variables stable and reduces the changing of elevator signal which is obvious in Figure 6. The aileron deflection is almost the same when using the two different strategies and showed in Figure 7. Considering the special shape of this UAV, the area of the control surface is not as large as conventional airplanes. Besides, it combines the functions of an elevator and an aileron. Therefore, the ENMPC can increase the throttle signal and decrease the elevator signal which is more suitable for this kind of UAVs with small area of the control surface.

To evaluate the robustness of the control law, longitudinal wind disturbance is introduced in the simulation. Figure 8 shows the height behavior under the unexpected wind disturbances which occur separately at the twentieth and sixtieth second and with a duration of 3 seconds. The first wind disturbance has a velocity of 2 m/s and a dip angle of 20°. The second one has the same dip angle but a velocity of −2 m/s.

Figure 8:

Height behavior under wind disturbance.

In Figure 8, when the first wind disturbance appears, height behavior has a variation of 1.14 m under PID control law which is 3 times as large as the variation under ENMPC strategy. Besides, the peak value of the dotted line arises at 23.65 second; however, when ENMPC strategy is employed, the peak value appears just right at 23 second. After the disturbance disappears at 23 second, the height behavior needs about 10 seconds to return to its initial value under PID control law which is 5 seconds longer than using ENMPC strategy. When the UAV comes across the wind disturbance at 60 second, height behavior shows a similarly change. It is obvious that ENMPC strategy has a stronger robustness compared with PID method.

Figures 9 –11 separately represent the velocity, attack angle, and pitch angle behavior when wind disturbances appear.

Figure 9:

Velocity behavior under wind disturbance.

Figure 10:

Attack angle behavior under wind disturbance.

Figure 11:

Pitch angle behavior under wind disturbance.

In PID method, velocity keeps a stable value from beginning to end. However, in ENMPC method, velocity decreases when the UAV comes across the updraft and increases under the downdraught. In other words, when the wind disturbance influences the aerodynamic force of the UAV, ENMPC method can adjust the velocity to bring an opposite influence. Therefore, ENMPC method has a better control effect on the flight attitude. It is obvious in Figures 10 and 11 that attack angle and pitch angle enjoy a relatively gentle change under ENMPC method compared with PID.

5. Conclusions

An explicit nonlinear model predictive control strategy was developed and simulated for a saucer-shaped lifting body UAV. As the lifting body UAV is difficult to be controlled compared with conventional ones but enjoys a wide application nowadays, designing a proper controller for special shaped UAV is a challenging but promising work. By comparing with PID strategy, ENMPC can get better controlling effect. On one hand, the control signals are calculated from nonlinear model without decoupling or linearization. On the other hand, the on-line optimisation process is removed from the MPC framework, which can accelerate the software development and simplify the onboard hardware. Due to these advantages of ENMPC, control framework has a low complexity, high reliability, and rapid response to the fast dynamic of UAVs. Besides, the throttle signal is not calculated only from velocity but influenced by all the variables. It can be adjusted immediately which reduces the deflection of the elevator. Therefore, the ENMPC strategy is more suitable for UAVs with smaller area of the control surfaces.

Our future work will focus on applying ENMPC strategy in actual flight tests using the saucer-shaped lifting body UAV. Furthermore, we are also interested in improving the control strategy to fit lifting body UAVs better.

Footnotes

Nomenclature

References

Wood

R. M.

and Bauer

X. S.

, “Flying wings/flying fuselages,” AIAA Paper 311, 2001.

Bowers

P. M.

, “Unconventional aircraft,” Tab Books, 1984.

Horton

V. W.

Eldredge

R. C.

, and Klein

R. E.

, “Flight-determined low-speed lift and drag characteristics of the lightweight M2-F1 lifting body,” National Aeronautics and Space Administration, 1965.

Duan

Luo

Shi

, “Hybrid particle swarm optimization and genetic algorithm for multi-UAV formation reconfiguration,” IEEE Computational Intelligence Magazine, vol. 8, no. 3, pp. 16–2, 2013.

Changhao

and Duan

, “Artificial bee colony optimized controller for unmanned rotorcraft pendulum,” Aircraft Engineering and Aerospace Technology, vol. 85, no. 2, pp. 104–114, 2013.

Ollero

and Merino

, “Control and perception techniques for aerial robotics,” Annual Reviews in Control, vol. 28, no. 2, pp. 167–178, 2004.

Mayne

D. Q.

Rawlings

J. B.

Rao

C. V.

, and Scokaert

P. O. M.

, “Constrained model predictive control: stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000.

Ikonen

Najim

, and Najim

, Advanced Process Identification, Control CRC PressI Llc, 2002.

Cloutier

J. R.

D'Souza

C. N.

, and Mracek

C. P.

, “Nonlinear regulation and nonlinear H-infinity control via the state-dependent Riccati equation technique: I- Theory,” in Proceedings of the 1st International Conference on Nonlinear Problems in Aviation and Aerospace, pp. 117–130, Daytona Beach, Fla, USA, 1996.

10.

Patwardhan

S. C.

and Madhavan

K. P.

, “Nonlinear model predictive control using second-order model approximation,” Industrial and Engineering Chemistry Research, vol. 32, no. 2, pp. 334–344, 1993.

11.

Xin

and Balakrishnan

S. N.

, “A new method for suboptimal control of a class of non-linear systems,” Optimal Control Applications and Methods, vol. 26, no. 2, pp. 55–83, 2005.

12.

Grancharova

and Johansen

T. A.

, Explicit Nonlinear Model Predictive Control, vol. 429 of Lecture Notes in Control and Information Sciences, 2012.

13.

Chen

W. H.

, “Predictive control of general nonlinear systems using approximation,” IEE Proceedings Control Theory and Applications, vol. 151, no. 2, pp. 137–144, 2004.

14.

Slegers

Kyle

, and Costello

, “Nonlinear model predictive control technique for unmanned air vehicles,” Journal of Guidance, Control, and Dynamics, vol. 29, no. 5, pp. 1179–1188, 2006.

15.

Liu

Chen

W. H.

, and Andrews

, “Explicit non-linear model predictive control for autonomous helicopters,” Journal of Aerospace Engineering, vol. 226, no. 9, pp. 1171–1182, 2012.

16.

Ricketts

R. H.

Noll

T. E.

Whitlow

Jr. , and Huttsell

L. J.

, “Overview of aeroelasticity studies for the National Aero-Space Plane,” in Proceedings of the 34th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, pp. 152–162, April 1993.

17.

Murphy

K. J.

Nowak

R. J.

Thompson

R. A.

Hollis

B. R.

, and Prabhu

, “X-33 hypersonic aerodynamic characteristics,” Journal of Spacecraft and Rockets, vol. 38, no. 5, pp. 670–683, 2001.

18.

Kyle

Evans

, and Costello

, “Atmospheric wind energy extraction by a small autonomous glider,” in Proceedings of the AIAA Atmospheric Flight Mechanics Conference Conference and Exhibit, pp. 1244–1270, San Francisco, Calif, USA, August 2005.

19.

Siegers

and Costello

, “Nonlinear model predictive control of a 6 DOF air vehicle,” in Proceedings of the AIAA Atmospheric Flight Mechanics Conference, pp. 1–15, 2005.

20.

Liu

Chen

W.-H.

, and Andrews

, “Tracking control of small-scale helicopters using explicit nonlinear MPC augmented with disturbance observers,” Control Engineering Practice, vol. 20, no. 3, pp. 258–268, 2012.

21.

Seanor

Campa

, “Design and flight testing evaluation of formation control laws,” IEEE Transactions on Control Systems Technology, vol. 14, no. 6, pp. 1105–1112, 2006.

22.

Leith

D. J.

and Leithead

W. E.

, “Gain-scheduled and nonlinear systems: dynamic analysis by velocity-based linearization families,” International Journal of Control, vol. 70, no. 2, pp. 289–317, 1998.

Explicit Nonlinear Model Predictive Control for a Saucer-Shaped Unmanned Aerial Vehicle

Abstract

1. Introduction

2. The Saucer-Shaped Lifting Body UAV Model

3. Explicit Nonlinear Model Predictive Control Law for the Saucer-Shaped UAV

3.1. Explicit Nonlinear Model Predictive Control

3.2. Saucer-Shaped UAV Control Law

4. Simulations Results

5. Conclusions

Footnotes

Nomenclature

References