Predictive Control for Active Model and Its Applications on Unmanned Helicopters

1. Introduction

Unmanned helicopters are increasingly popular platforms for unmanned aerial vehicles (UAVs). With the abilities such as hovering, taking off and landing vertically, unmanned helicopters extend the potential applications of UAVs. However, due to the complex mechanism and complicated aero-flow during flight, it is almost impossible to accurately model the dynamics of an unmanned helicopter in full flight envelope, and the significant model uncertainties associated with a nominal model may degrade the performance and even stability of an onboard controller.

Due to the difficulty in obtaining a high fidelity full envelope model, the multi-mode modeling technique has been proposed for rotor aircrafts, such as tilt-rotor aircraft XV-15 [1], helicopter BO-105 [2], UH-60 [3], R-50 [4] and X-Cell [5]. The mode-dependent model, which is identified and simplified according to a specific flight mode, such as hovering, cruising, taking off and landing, can be used for control design for the corresponding flight mode. However, the mode-dependent control suffers from at least two problems: one is the difficulty in accommodating the mode transition dynamics, and the other is the compensation of the ‘model drift’ due to flight dynamics change within one particular mode. Up to now, for the purpose of practical implementation, the mode transition problem can be partially dealt with by limiting the mode switching conditions [6], e.g., mode change is made through hovering mode.

Robust and adaptive control techniques [7-8], on the other hand, have been used to deal with the ‘model-shift’ within a flight mode. However, such control schemes normally need to know the boundary of internal and external uncertainties and relative noise distribution, which are difficult to identify accurately for a helicopter in full flight envelope. Although online identification technology can be used to obtain the real-time dynamics and disturbance, it is a large burden for the flight computer to reconstruct the robust controllers and reach the requested control period (>50Hz) for sampling and actuating due to the complex calculation of the robust/adaptive optimization process [9-10] and the strict weight limits of micro flight computers.

Besides the model uncertainties, another critical problem that limits the control performance of a helicopter is the time delay between the actuator command and the generation of relative aerodynamic force/torque [11], which will be called aerodynamics-delay/time-delay in the following sections. Normally, this time delay may cause reduced feedback gain of a model-based controller and result in poor robustness [12-13], i.e., sensitive to disturbances.

In recent years, the encouraging achievement in sequential estimation makes it an important direction for online modeling and model-reference control [14]. Among stochastic estimations, the most popular one is the Kalman-type filters (KFs) [15, 16, and 17]. Although widely used, the KFs suffer from sensitivity to bias and divergence in the estimates, relying on assumptions on statistic distribution such as white noise and known mean or covariance for optimal estimation. In many cases, it is more practical to assume that the noises or uncertainties are unknown but bounded (UBB). In view of this, the set-membership filter (SMF), which computes a compact feasible set in which the true state or parameter lies only under the UBB noise assumption, provides an attractive alternative [18-19].

On the control issue, model predictive control (MPC) can compensate for the aerodynamics delay and does not require a high accuracy reference nonlinear model [20]. Among these methods, linear generalized predictive control (GPC) has become one of the most popular MPC methods in industry and academia. However, the normal GPC is sensitive to process noise and model errors [21], which are unknown but bounded for helicopters when sudden ‘mode change’ happen and model-drift in full flight envelope. This makes the prediction biased, and results in the non-optimal process of controller solving.

In this paper, for realizing the coupling control of unmanned helicopters in full flight envelope, an active modeling based controller is developed based on a modified generalized predictive control and adaptive set-membership filter estimation (ASMF). The time varying model error and its boundary are estimated by the adaptive set-member filter, which is first proposed in [19]. Incremental prediction process and dimension reduction method is embedded into traditional GPC, which can decreases the computation burden and maintain prediction unbiased when ‘mode change’ happens. Based on this active estimation and the modified GPC controller, a novel optimal strategy for on-line compensation of model error is developed. Thus, aggressive flight can be achieved only based on the hovering model with time-delay terms. Using the identified hovering dynamics model as nominal model for controller, flight experiments have been conducted to test the performance of the proposed controller in full flight envelope on our UAV platform, and experimental results have demonstrated the effectiveness of the proposed method.

2. Active model based control scheme and reference model of a helicopter

Fig. 1 illustrates the active model based control scheme. The error between the reference model and the actual dynamics of the controlled plant is estimated by an on-line modeling strategy. The control, which is designed according to the reference model, should be able to compensate the estimated model error and it in real time. In the followings of this paper, we use the ASMF as the active modeling algorithm and the modified GPC as the control.

For normal missions of an unmanned helicopter, the flight modes include hovering (velocity under 5m/s), cruising (velocity above 5m/s), taking off and landing (distance to the ground is below 3m while significant ground effect exists) and the transitions among these modes. A reference model is typically obtained by linearizing the nonlinear dynamics of a helicopter at one flying mode. The model errors from linearization, external disturbance, simplification, and un-modeled dynamics can be considered as additional process noise [22]. Thus, a linearized state-space model for helicopter dynamics in full flight envelope can be formulated as

Figure 1.

The scheme of active model based control

$$\{\begin{array}{l}{\dot{X}}_t=A_0{X}_t+B_0{U}_{t-k}+B_{f}f(X_t,{\dot{X}}_t,W_t)\\ Y_t=C{X}_t\end{array}$$E1

where $$X\in R^{13}$$ is the state, including 3-axis velocity, pitch and roll angle, 3-axis angle rate, flapping angles of main rotor and stabilizer bar, and the feedback of yaw gyro. $$Y_t\in R^8$$is the output, including 3-axis velocity, pitch and roll angle and 3-axis angle rate, $$A_0$$and $$B_0$$ contain parameters that can be identified in different flight modes, and we use them to describe the parameters in hovering mode. $$U\in R^4$$is the control input vector. $$C\in R^{13\times 8}$$is the output matrix, $$k\in R$$is the time-delay for the driving system. The detail of building the nominal model and physical meaning of parameters is explanted in Appendix A.

To describe the dynamics change, in equation (1), here, we introduce $$f(X_t,{\dot{X}}_t,W_t)\in R^{13}$$ to represent the time varying model error in full flight envelope, and $$W_t\in R^{13}$$ is the process noise.

The following two sections, based on model (1) will describe the way to estimate $$f(X_t,{\dot{X}}_t,W_t)$$ and to compensate for model errors from process noise, parameters change, control delay and flight mode change in real applications.

3. ASMF based active model error estimation

As illustrated in Fig.1, adopting the active modeling process to get the model error f and system state X is the basis for elimination of the model error. Controller can only work based on nominal model and feedback of state and model error from active modeling process. In this section, the active modeling process is built based on an adaptive set-membership filter (ASMF) [19] since the UBB process noise.

First, we must obtain the reference equation for estimation. Compared with the sampling frequency (often >50Hz for flight control) of the control system, the model error $$f(X,\dot{X},W)$$ can be considered as a slow-varying vector, which means

$$f_{t+1}=f_t+h_t$$E2

where $$f_t$$is the sampling value of $$f(X,\dot{X},W)$$at sampling time t, and $$h_t$$is the assumed unknown but bounded (UBB) process noise.

Let the extended sampling state

$$X_t^a={\left(\begin{array}{cc}{X}_t^T& f_t^T\end{array}\right)}^T$$E3

Then, we can obtain the discrete equation from Eq. (1) as

$$\{\begin{array}{l}{X}_{t+1}^a=A_d^a{X}_t^a+B_d^a{U}_t+W_t^a\\ Y_t=C_d^a{X}_t^a+V_t\end{array}$$E4

where$$A_d^a=\left(\begin{array}{cc}{A}_d& B_f\\ 0_{13\times 13}& I_{13\times 13}\end{array}\right)$$, $$B_d^a=\left(\begin{array}{c}{B}_d\\ 0_{13\times 4}\end{array}\right)$$, $$C_d^a=\left(\begin{array}{cc}{C}_d& 0_{8\times 13}\end{array}\right)$$, $$W_t^a={\left(\begin{array}{cc}{W}_t^T& h_t^T\end{array}\right)}^T$$, $$B_f=I_{13\times 13}$$and $$f_t$$ is a $$13\times 1$$vector for model errors. Here, t is the sampling time, $$I_{m\times m}$$is the m×m unit matrix and $$0_{m\times n}$$ is the m×n zero matrix. $$\{A_d,B_d,C_d\}$$is the discrete expression of system$$\{A_0,B_0,C\}$$. Here, time-delay k is ignored during the estimate process, and the compensation method will be discussed in the next part on modified GPC.

The model error i.e., f in Eq. (1), comes from the linearization while neglecting the coupling dynamics and uncertainties, and also the $$A_0$$ and $$B_0$$ because they are identified with respect to a specific flight mode, here hovering mode is selected as nominal flight mode since easy identification. Therefore, both the model error and the process noise $$W^a$$ are vehicle dynamics and flight states dependent, and do the following assumption

Assumption:$$W^a$$does not necessarily have a normal distribution.

Thus, the Kalman type filter cannot be applied, and adaptive set-membership filter, which is developed for UUB process noise and can get the uncertain boundaries of the states, is considered to estimate the states and model errors here.

In this section we only present the result of ASMF and please refer to [19] for the details about ASMF. With respect to Eq. (2), we can build the adaptive set-membership filter as Eq. (3), where $$Q^a$$ and $$R^a$$ are the initial elliptical boundary of process and measurement noise respectively, $$r_m$$is the maximum eigenvalue of R, $$p_m$$is the maximum eigenvalue of$$C_d^a{P}_{t|t-1}{C}_d^a{}^T$$, Tr(•) is the trace of a matrix, $${\delta }_t$$and $${\beta }_t$$ are the adaptive parameters of the filter. We can also obtain the boundary of the ith element $${\widehat{X}}_i^a$$of extended state $${\widehat{X}}_{t|t}^a$$as$$\left({\widehat{X}}_i^a-\sqrt{P_{ii}},{\widehat{X}}_i^a+\sqrt{P_{ii}}\right)$$, where $$P_{ii}$$ is the i-th diagonal element of matrix P.

$$\{\begin{array}{l}{\rho }_t=\frac{\sqrt{r_{mt}}}{\sqrt{r_{mt}}+\sqrt{p_{mt}}}\\ W_t=\frac{C_d^a{P}_{t|t-1}{C}_d^a{}^T}{1-{\rho }_t}+\frac{R^a}{{\rho }_t}\\ K_t^e=\frac{P_{t|t-1}{C}_d^a{}^T{W}_t^{-1}}{1-{\rho }_t}\\ {\delta }_t=1-{(}^{Y_t}{W}_t^{-1}(Y_t-C_d^a{\widehat{X}}_{t|t-1}^a)\\ {\widehat{X}}_{t|t}^a={\widehat{X}}_{t|t-1}^a+K_t^e(Y_t-C_d^a{\widehat{X}}_{t|t-1}^a)\\ P_{t|t}={\delta }_t(\frac{P_{t|t-1}}{1-{\rho }_t}-\frac{P_{t|t-1}}{1-{\rho }_t}{C}_d^a{}^T{W}_t^{-1}{C}_d^a\frac{P_{t|t-1}}{1-{\rho }_t})\\ {\widehat{X}}_{t+1|t}^a=A_d^a{\widehat{X}}_{t|t}^a+B_d^a{U}_t\\ {\beta }_t=\frac{\sqrt{Tr(Q^a)}}{\sqrt{Tr(Q^a)}+\sqrt{Tr(A_d^a{P}_{t|t}{A}_d^a{}^T)}}\\ P_{t+1|t}=\frac{A_d^a{P}_{t|t}{A}_d^a{}^T}{1-{\beta }_t}+\frac{Q^a}{{\beta }_t}\end{array}$$E5

4. Modified GPC for unmanned helicopters

To eliminate the negative influence of model errors and control delay in flight, besides the active estimation algorithm like ASMF that does not require a normal distribution assumption, an effective control algorithm has to be designed according to the reference model of Eq. (1) while adopting the on-line estimation of f as compensation.

We describe the normal GPC in Section 4.1, and then, the modified scheme is proposed in Section 4.2 & 4.3 to eliminate the negative influence of model errors in real applications.

4.1. Preliminary work for generalized predictive control

Generally, for a linear system with actuator time delay like,

$$\{\begin{array}{l}{X}_{t+1}=A_d{X}_t+B_d{u}_{t-k}+W_t\\ y_t=C_d{X}_t\end{array}$$E6

where $$X_t\in R^{n\times 1}$$is the system state vector at sampling time t, $$y_t\in R^{l\times 1}$$is the output vector, $$u_t\in R^{m\times 1}$$is the control input vector, k is the actuators’ time-delay and $$W_t$$ is process noise; traditional Generalized Predictive Control (GPC) [23] can be designed as:

Step I: Make prediction

Firstly, for the case that predictive step i is less than time-delay k (i.e., the time instant that system behavior cannot be regulated through current and future control action), prediction can be denoted as following equation,

$${\widehat{X}}_{t+i|t}=A_d{\widehat{X}}_{t+i-1|t}+B_d{u}_{t+i-1-k}\triangleq {\widehat{X}}_{t+i|t}^1$$E7

where$${\widehat{X}}_{t+i}|t$$ is the prediction state at time t+i, the superscript 1 denotes that the part of predicted variable that is independent of the current and future’s control actions.

Secondly, for the case that prediction step i is larger than the time delay k,

$$\begin{array}{l}{\widehat{X}}_{t+k+i|t}=A_d{\widehat{X}}_{t+k+i-1|t}+B_d{u}_{t+i}\\ =A_d{\widehat{X}}_{t+k+i-1|t}^1+{\displaystyle \sum _{n=0}^i{A}_d^n{B}_d}{u}_{t+n}\\ =A_d{\widehat{X}}_{t+k+i|t}^1+{\displaystyle \sum _{n=0}^i{A}_d^n{B}_d}{u}_{t+n},1\le i\le p\end{array}$$E8

where p is the prediction range; similarly, $${\widehat{X}}_{t+k+i-1|t}^1$$denotes the sub-variable of $${\widehat{X}}_{t+k+i-1|t}$$that is independent of the current and future’s control actions.

Step II: Receding horizon optimization

After making prediction, the control vector can be obtained by minimize the following cost function:

$$J={(R^x{}_t-X_t^v)}^T(R^x{}_t-X_t^v)+U_t^T\gamma U_t$$E9

And the optimal control inputs can be denoted as,

$$U_t^{*}={(G_0^T{G}_0+\gamma )}^{-1}{G}_0^T(R^x{}_t-X_t^1)$$E10

where $$G_0$$is the predictive matrix, $$X_t^v$$is the predictive state vector,$$X_t^1$$ is the known vector inside$$X_t^v$$, $\lambda$is the weight of control input, and $$R_t^x$$is the reference of system states. The detailed definition of these matrixes can be referenced in [23].

Step III: Control implementation

The first element of vector $U_t^{*}$is used as the control to the real plant. After that, go back to step I at the next time instant.

However, with application to the unmanned helicopters, this kind of GPC algorithm has the following three disadvantages, which will be solved in the next two sections:

1. It cannot reject the influence of working mode changes, i.e., if

$$\begin{array}{l}{X}_t=x_t-x_0>\pi (x_0,u_0)\\ U_t=u_t-u_0>\pi (x_0,u_0)\end{array}$$E11

where $$(x_0,u_0)$$is the current operation point, which cannot be ensured on-line, $$\pi (x_0,u_0)$$is the valid range for model linearization and $$x_t$$is the absolute state at time t, $$u_t$$is the absolute control input at time t. The biased prediction, due to the changing operation point$$(x_0,u_0)$$, will bring steady errors for velocity tracking.

1. Normal GPC is sensitive to mismatch of the nominal model, which means slow change in parameters $$(A_d,B_d)$$may result in prediction error and unstable control.

2. The transient model errors of the nominal model from external disturbance, estimated by ASMF, cannot be eliminated. And this will also result in the non-minimum variance and the instability of the closed control loop.

4.2. Stationary increment predictive control

To reject the influence of working mode change and sensitivity to nominal parameters change in real application, i.e. the problem 1) and 2) in Section 4.1, we assume that the process noise$$W_t$$’s increment in Eq. (4) is a stationary random process, which means

$$W_t^0\triangleq \Delta W_t=W_t-W_{t-1}$$E12

is normal distribution. Where $$\Delta =1-q^{-1}$$ is the difference operator; $$q^{-1}$$is one-step delay factor. Thus, Eq. (4) can be rewritten as follows,

$$\Delta X_{t+1}=A_d\Delta X_t+B_d\Delta u_{t-k}+W_t^0$$E13

Consider

$$\begin{array}{l}\Delta X_t=(x_t-x_0)-(x_{t-1}-x_0)=\Delta x_t\\ \Delta U_t=(u_t-u_0)-(u_{t-1}-u_0)=\Delta u_t\end{array}$$E14

if behavior prediction is made based on Eq. (11), only the absolute state$$x_t$$and control input$$u_t$$, which can be measured or estimated directly from sensors, are used and the current operation point $$(x_0,u_0)$$ disappears in prediction. Thus, the problem of biased prediction due to changing of working point, i.e., problem 1), can be solved.

Otherwise, according to the process of traditional GPC, the set-point $R_t^x$ must be obtained for every prediction step, and this is often set as current reference states. However, for helicopter system, only measurable outputs are cared, such as position, velocity and etc; and the internal states, such as rotor’s pitch angle and yaw gyro’s feedback and so on, are coupled with the measurable states/outputs, and cannot be set independently. Others, this reference input often comes from position track planning, which changes quickly for flight and often cause a step-like signal for tracking. To avoid the step signal reference tracking, which is dangerous for unmanned helicopter system, we use a low pass filter to calculate the set-point inputs of the output in the future i-th step, i=1, …, p.

Let $$S{P}_t\in R^{l\times 1}$$ be the set-point input at time t, then we have

$$r_{t+k+i}=S{P}_t+\alpha (r_{t+k+i-1}-S{P}_t),1\le i\le p$$E15

where α is the cut-off frequency of the filter, the initial value$$r_{t+k}={\widehat{y}}_{t+k|t}$$, $$r_{t+k+i}$$is the i-th set-point input, and $${\widehat{y}}_{t+k|t}$$is the estimate of output at time t+k.

Thus, the set-point problem is solved and the output prediction can be implanted based on increment model (11) as follows:

When the prediction step i is less than time-delay k,

$${\widehat{X}}_{t+i|t}={\widehat{X}}_{t+i-1|t}+A_d\Delta {\widehat{X}}_{t+i-1|t}+B_d\Delta u_{t+i-1-k}={\widehat{X}}_{t+i|t}^1$$E16

When the prediction step is larger than time-delay k, let

$$\Delta {\widehat{X}}_{t+i|t}^1={\widehat{X}}_{t+i|t}^1-{\widehat{X}}_{t+i-1|t}^1$$E17

Then,

$$\begin{array}{l}{\widehat{X}}_{t+k+i|t}={\widehat{X}}_{t+k+i-1|t}+A_d\Delta {\widehat{X}}_{t+k+i-1|t}+B_d\Delta u_{t+i}\\ ={\widehat{X}}_{t+k+i-1|t}^1+A_d\Delta {\widehat{X}}_{t+k+i-1|t}^1\\ +{\displaystyle \sum _{m=0}^{i-1}\left\{\left\{{\displaystyle \sum _{n=0}^{i-1-m}{A}_d^n{B}_d}\right\}\Delta u_{t+m}\right\}}\\ ={\widehat{X}}_{t+k+i|t}^1+{\displaystyle \sum _{m=0}^{i-1}\left\{\left\{{\displaystyle \sum _{n=0}^{i-1-m}{A}_d^n{B}_d}\right\}\Delta u_{t+m}\right\}},1\le i\le p\end{array}$$E18

Hence, the above problem 1), which comes from working mode change, is solved because $$x_0$$ disappears in predictive equation (14).

We can obtain the following prediction matrix for the output, which is often cared in helicopter tracking problem, from Eq. (12) and (13):

$$\begin{array}{l}{\widehat{Y}}_t={\left(\begin{array}{cccc}{\widehat{y}}_{t+k+1|t}& {\widehat{y}}_{t+k+2|t}& \mathrm{...}& {\widehat{y}}_{t+k+p|t}\end{array}\right)}^T\\ ={\left(\begin{array}{cccc}{C}_d{\widehat{X}}_{t+k+1|t}^1& C_d{\widehat{X}}_{t+k+2|t}^1& \mathrm{...}& C_d{\widehat{X}}_{t+k+p|t}^1\end{array}\right)}^T\\ +G{\left(\begin{array}{cccc}\Delta u_t{}^T& \Delta u_{t+1}{}^T& \mathrm{...}& \Delta u_{t+p-1}{}^T\end{array}\right)}^T\\ =Y_t^1+G\Delta U\end{array}$$E19

where $$Y_t^1$$ is the known part of p steps’ prediction, which cannot be influenced by current control input, and matrix G has the following form:

$$G=\left(\begin{array}{cccc}{C}_d{B}_d& 0& \mathrm{...}& 0\\ C_d{B}_d+C_d{A}_d{B}_d& C_d{B}_d& \mathrm{...}& 0\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ C_d{\displaystyle \sum _{i=0}^{p-1}{A}_d^i{B}_d}& C_d{\displaystyle \sum _{i=0}^{p-2}{A}_d^i{B}_d}& \mathrm{...}& C_d{B}_d\end{array}\right)$$E20

Compared with the normal GPC, the prediction of SIPC has better characteristics that can be described by the following theorem, which solves the above problem 2) in Section IV.A.

Theorem: for nominal model (11), when the nominal model parameters$$(A_d,B_d)$$change into$$(A_{dr},B_{dr})$$.

1. $$\exists M,N>0\in R$$, let the matrix norms satisfy

   $$\Vert A_d\Vert <M,\Vert B_d\Vert <M$$E21
   $$\Vert A_{dr}\Vert <N,\Vert B_{dr}\Vert <N$$E22

2. Define$$R_{\max }\{•\}$$is the operator for the maximum of eigenvalue of matrix$$•$$.

Thus, if

$$\begin{array}{l}{R}_{\max }\left\{A_{dr}\right\}<0\\ R_{\max }\left\{A_d\right\}<0\end{array}$$E23

Then, the state prediction obtained by Eqs. (13-14) maintains unbiased, and the characteristic is also guaranteed in traditional GPC conditions, i.e. Eq. (4), where $$W_t$$is normal distribution.

Proof: See Appendix B.

In Eq. (14),$$\Delta U$$ , including p control inputs, need to be optimized, while only the first one is used for control. This will occupy a great deal of computation resource and result in very low computational efficiency, especially with respect to the fast applications.

In order to reduce the computational burden of Eq. (14), we propose here a ‘step plan’ technique,

$$\Delta u_{t+i+1}=\beta \Delta u_{t+i}$$E24

where $$\beta$$ is an $$m\times m$$ diagonal matrix presenting the length of one step, which will be a parameter to be selected. Then, we can simplify Eq. (14) by only calculating the unknown control, which has smaller dimensions.

$$\begin{array}{l}{\widehat{Y}}_t=Y_t^1+G{\left(\begin{array}{cccc}{I}_{m\times m}& \beta & \mathrm{...}& {\beta }^{p-1}\end{array}\right)}^T\Delta u_t\\ =Y_t^1+G_2\Delta u_t\end{array}$$E25

where $$I_{m\times m}$$ is an $$m\times m$$ unit matrix. Thus, the number of the unknown control input vector (from current time t to the future time t+p-1) is reduced from p to 1, and the dimension of predictive matrix is changed from $$pl\times pm$$to$$pl\times m$$. This reduction brings low computer memory consuming and simplifies the receding horizon optimization in the following calculation.

To complete the horizon optimization and obtain the control input, the cost function of the stationary increment predictive control is designed as:

$$J={(R_t-{\widehat{Y}}_t)}^{T}W(R_t-{\widehat{Y}}_t)+\Delta u_t{}^T\lambda \Delta u_t$$E26

where$R_t={\left(\begin{array}{cccc}{r}_{t+k+1}^T& r_{t+k+2}^T& \mathrm{...}& r_{t+k+p}^T\end{array}\right)}^T$, $W\in R^{lp\times lp}$is the weight matrix for tracking error, and $\lambda \in R^{m\times m}$is the weight matrix of the control increment.

In order to minimize the cost function of Eq. (19), we can calculate the control vector as follows:

$$\begin{array}{l}\Delta u_t={(}^{G_2^T}{G}_2^{T}W(R_t-Y_t^1)\\ =K_f(R_t-Y_t^1)\end{array}$$E27

where $K_f={(G_2^{T}W{G}_2+\lambda )}^{-1}{G}_2^{T}W$ can be completed offline.

Consequently, the proposed stationary increment predictive controller (SIPC) can be designed as followings.

Step I: Make increment prediction

Based on the current and history measure value, use Eqs. (13-15) to obtain the prediction for future output $${\widehat{Y}}_t$$and initial plan point

$$r_{t+k}={\widehat{y}}_{t+k|t}$$E28

Step II: Plan for the set-point input

Use Eq. (12) to plan the future set-points, and obtain

$$R_t={\left(\begin{array}{cccc}{r}_{t+k+1}^T& r_{t+k+2}^T& \mathrm{...}& r_{t+k+p}^T\end{array}\right)}^T$$E29

Step III: Receding horizon optimization

Calculate the control increment$\Delta u_t$, based on Eq. (20).

Step IV: Control implementation

Current control input$u_t=u_{t-1}+\Delta u_t$, which is used as the control to the real plant. After that, go back to step I at the next time instant.

Thus, for real implementation, only the prediction of Eq. (13-15), the intenerating of Eq. (12), and the control law (20) need to be calculated online, thus the real time computation load, and steady tracking error are both reduced greatly compared with GPC, and the real test in section V has shown its feasibility.

The model error, problem 3), will be compensated by an online optimal strategy, which will be described later.

4.3. Optimal strategy for model error compensation

In order to compensate the model error in Eq. (1), the control vector has to match the following equation, which can be directly obtained from Eq. (1):

$$B_d{U}_t+B_f{f}_t=B_d{U}_t^0$$E30

where $$U_t^0$$ is the control vector need to be calculated by the predictive controller in section 4.2, designed based on the original model (1) without the model error f.

The control input at sampling time t cannot be solved directly from Eq. (21), because:

1. Eq. (21) is difficult to be implemented because the dimension of$$U_t$$ is less than that of$$f_t$$. Thus, only the approximate solution can be obtained with respect to (21);

2. $$f_t$$is actually an uncertainty set, an static optimal problem must be considered.

Thus, we introduce the following cost function with quadratic form to solve the above problem 1).

$$\begin{array}{l}{U}_t^{*}=\underset{U_t}{\arg \min }{J}_t(U_t)\\ J_t(U_t)\triangleq {\left(B_d{U}_t+B_f{f}_t-B_d{U}_t^0\right)}^{T}H\left(B_d{U}_t+B_f{f}_t-B_d{U}_t^0\right)\end{array}$$E31

where H is a weight matrix, which can be selected.

On the other hand, $$f_t$$is obtained from the ASMF algorithm introduced in section III, thus its convergence is very important for the validity of the whole controller. Actually, the convergence of ASMF algorithm is also influenced by the control action$$U_t$$. This is because the stability of the ASMF can be represented by the filter parameter$${\delta }_t$$, while $${\delta }_t$$ in Eq. (3) can be rewritten as follows,

$$\begin{array}{l}{\delta }_t=1-{(}^{Y_t}{W}_t^{-1}(Y_t-C_d^a{\widehat{X}}_{t|t-1}^a)\\ =1-{(}^{Y_{t+1}}\\ W_t^{-1}(Y_{t+1}-C_d^a(A_d^a{\widehat{X}}_{t|t}^a+B_d^a{U}_t))]\end{array}$$E32

In [19], it has been shown the stability of the ASMF can be represented by the filter parameter$${\delta }_t$$, i.e., the ASMF is stable when$${\delta }_t>0$$.

Firstly, define

$$\begin{array}{l}{J}_t^{\delta }(U_t,Y_{t+1})\triangleq {(}^{Y_{t+1}}\\ W_t^{-1}(Y_{t+1}-C_d^a(A_d^a{\widehat{X}}_{t|t}^a+B_d^a{U}_t))]\end{array}$$E33

Thus, from Eq. (23), in order to maintain$${\delta }_{t+1}>0$$, the maximum value of $$J_t^{\delta }(U_t,Y_{t+1})$$ with respect to $${\widehat{X}}_{t|t}^a$$ should be less than or equal to 1, i.e.,

$$\begin{array}{l}{J}_t^{\delta *}(U_t,Y_{t+1})=\underset{{\widehat{X}}_{t|t}^a}{\max }{J}_t^{\delta }(U_t,Y_{t+1})\\ =\underset{{\widehat{X}}_{t|t}^a}{\max }\left\{{\left[Y_{t+1}-C_d^a\left(A_d^a{\widehat{X}}_{t|t}^a+B_d^a{U}_t\right)\right]}^T{W}_t^{-1}\left[Y_{t+1}-C_d^a\left(A_d^a{\widehat{X}}_{t|t}^a+B_d^a{U}_t\right)\right]\right\}\le 1\end{array}$$E34

In general, larger $${\delta }_t$$ often means more rapid convergence of ASMF algorithm. That is, we should select an Ut to make $$J_t^{\delta *}(U_t,Y_{t+1})$$ small as far as possible, that is,

$$J_t^{\delta *}(Y_{t+1})=\underset{U_t}{\min }{J}_t^{\delta *}(U_t,Y_{t+1})$$E35

We introduce the following cost function $$J_t(U_t)$$ with consideration of both (22) and (25) at the same time:

$$\begin{array}{l}{}^{*}{U}_t=\underset{U_t}{{\displaystyle \arg \min }}{\overline{J}}_t(U_t)\\ {\overline{J}}_t(U_t)\triangleq J_t(U_t)+\alpha J_t^{\delta *}(U_t,Y_{t+1})\end{array}$$E36

where $$\alpha =1-{\delta }_t\in R$$ are the positive definite weight matrix. To minimize$$J_t(U_t)$$, considering$$J_t(U_t)>0$$, the control can be obtained at$$\raisebox{1ex}{$\partial J_t(U_t)$}\!\left/ \!\raisebox{-1ex}{$\partial U_t$}\right.=0$$, i.e.,

$$\raisebox{1ex}{$\partial J_t(U_t)$}\!\left/ \!\raisebox{-1ex}{$\partial U_t$}\right.=2(M{U}_t+N)$$E37

where

$$M=B_d^{T}H{B}_d+\alpha B_d^a{}^T{C}_d^a{}^T{W}_t^{-1}{C}_d^a{B}_d^a$$E38

$$N=B_d^{T}H(B_f{f}_t-B_d{U}_t^0)-\alpha B_d^a{}^T{C}_d^a{}^T{W}_t^{-1}{Y}_{t+1}$$E39

Here H can be selected as$$H={\delta }_t{C}_d^T{C}_d$$. Thus, we can obtain the optimal control that minimizes $$J_t(U_t)$$ as:

$$\begin{array}{l}{U}_t(Y_{t+1})=-M^{-1}N\\ ={(}^{B_d^T}\\ \left[\alpha B_d^a{}^T{C}_d^a{}^T{W}_t^{-1}{Y}_{t+1}-B_d^{T}H(B_f{f}_t-B_d{U}_t^0)\right]\end{array}$$E40

For the unknown measurement at time t+1 in Eq. (24), we consider that the control system is stable, so,$$Y_{t+1}\in \Delta (Y_t)$$. Here, $$\Delta (Y_t)$$is the elliptical domain of$$Y_t$$. Because $$J_t^{\delta }(U_t,Y_{t+1})$$ in Eq. (25) is positive definite, its maximum value point must be on the boundary, which can be estimated by the ASMF. Thus, we first define array $$S_t^i$$to include the estimate of the i-th element’s two boundary endpoints as

$$S_t^i\triangleq \left\{{\widehat{Y}}_{t+1}^i\left|\right\{Y_t^i+{(-1)}^h(\underset{j_l=\{\pm \sqrt{p_{ll}},l=1,\mathrm{...},13\}}{Max}{\left|C_{d}Col\left\{j\right\}\right|}_i)\}\right\}$$E41

where $$Y_t^i$$ is the i-th element in the vector$$Y_t$$, $${\widehat{Y}}_{t+1}^i$$is the corresponding output$$Y_{t+1}$$’s endpoints estimation. For set$$S_t^i$$, $$i\in \{1,2,\mathrm{...},8\}$$and h is 0 or 1 for every i, $${|•|}_i$$is the operator for absolute value of the i-th element in vector$$•$$, and the function $$Col\left\{j\right\}$$ is defined as follows:

$$Col\left\{j\right\}={\left(\begin{array}{ccc}{j}_1& \mathrm{...}& j_{13}\end{array}\right)}^T$$E42

Then, we define a set $$S_t$$to describe all possible endpoint vector of the$$Y_{t+1}$$ as

$$S_t\triangleq \left\{{\widehat{Y}}_{{}_{t+1}}^{EP}|\left(\begin{array}{ccc}{S}_t^1& \mathrm{...}& S_t^{13}\end{array}\right)\right\}$$E43

where $${\widehat{Y}}_{t+1}^{EP}$$is the possible endpoint (EP) for output $$Y_{t+1}$$at next sampling time t+1.

Thus, the proposed active modeling based predictive controller can be implemented by using the following steps:

Step I: Make increment prediction

Based on the current estimated state$${\widehat{X}}_{t|t}^a$$, use the stationary increment predictive controller, as in section 4.2, to obtain the nominal control input$$U_t^0$$;

Step II: Model error estimation and elimination

Based on$$U_t^0$$, compute the optimal control input$${}^{*}{U}_t$$:

Estimate the values and boundaries of state $$X_t$$ and model error$$f_t$$, using ASMF in (3);

Calculate the corresponding $$U_t({\widehat{Y}}_{t+1}^{EP})$$ for every $${\widehat{Y}}_{t+1}^{EP}$$ in set$$S_t$$ by Eq. (29);

For every $$U_t({\widehat{Y}}_{t+1}^{EP})$$in step 1), use Eq. (24) to obtain the maximum of function$$J_t^{\delta }(U_t({\widehat{Y}}_{t+1}^{EP}),{\widehat{Y}}_{t+1}^{EP})$$, and get the $${}^{*}{\widehat{Y}}_{t+1}^{EP}$$to let

$${}^{*}{\widehat{Y}}_{t+1}^{EP}=\underset{{\widehat{Y}}_{t+1}^{EP}\in S_t}{\arg Max}\left\{J_t^{\delta }\{U_t({\widehat{Y}}_{t+1}^{EP}),{\widehat{Y}}_{t+1}^{EP}\}\right\}$$E44

The corresponding $$U_t{(}^{*}{\widehat{Y}}_{t+1}^{EP})$$ is the optimal control $${}^{*}{U}_t$$ at time t, i.e.$${}^{*}{U}_t=U_t{(}^{*}{\widehat{Y}}_{t+1}^{EP})$$.

Step III: Receding horizon strategy

Go back to step I at the next time instant.

5. Flight test

5.1. Flight test platform

All flight tests are conducted on the Servoheli-40 setup, which was developed in the State Key Laboratory, SIACAS. It is equipped with a 3-axis gyro, a 3-axis accelerometer, a compass and a GPS. The sensory data can be sampled and stored into an SD card through an onboard DSP. Tab.1 shows the physical characteristics of SERVOHELI-40 small-size helicopter. More details of this experimental platform can be found in [24].

Figure 2.

SERVOHELI-40 small-size helicopter platform

Length	2.12m
Height	0.73m
Main rotor diameter	2.15m
Stabilizer bar diameter	0.75m
Rotor speed	1450rpm
Dry weight	20kg
Engine	2-stroke, air cooled
Flight time	45 min

Table 1.

Physical characteristics of SERVOHELI-40 small-size helicopter

5.2. Experiment for the verification of model error estimate when mode-change

We use the identified hovering parameters, through frequency estimate [25], as the nominal model for hovering dynamics of the ServoHeli-40 platform. The model accuracy is verified in hovering mode (speed less than 3m/s) and cruising mode (speed more than 5m/s), the results for lateral velocity are shown in Fig.3a.

Fig. 3 further shows the model difference due to mode change, where the red lines are the results calculated by the identified model with the inputs of hovering and cruising actuations, respectively, and blue lines are the measurements of the onboard sensors. Comparison shows that the hovering model outputs match the hovering state closely, but clear differences occur while being compared to the cruising state, even though the cruising actuations are used as the model inputs. This is the model error when flight mode is changed.

Figure 3.

Model difference due to mode change: (a) hovering conditions; (b) cruising conditions

To verify the accuracy of the estimate of the model error, described in Fig.3, the following experiment is designed:

1. Actuate the longitudinal control loop to keep the speed more than 5 meter per second;

2. Get the lateral model error value and boundaries through ASMF, and add them to the hovering model we built above;

3. Compare the model output before and after compensation for model error.

This process of experiment can be described by Fig.4, and the results are shown in Fig.5. Fig.5a shows that model output (red line) cannot describe the cruising dynamics due to the model error when ‘mode-change’, similar with Fig.3b; however, after compensation, shown in Fig.5b, the model output (red line) is very close with real cruising dynamics (blue line), and the uncertain boundaries can include the changing lateral speed, which mean that the proposed estimation method can obtain the model error and range accurately by ASMF when mode-change.

Figure 4.

The experiment process for model-error estimate

Figure 5.

Model output before/after compensation: (a) before compensation; (b) after compensation

5.3. Flight experiment for the comparison of GPC SIPC and AMSIPC when sudden mode-change

In Section 5.2, the model-error occurrence and the accuracy of the proposed method for estimation are verified. So, the next is the performance of the proposed controller in real flight. In this section, the performance of the modified GPC (Generalized Predictive Control, designed in Section 4.1), SIPC (Stationary Increment Predictive Control, designed in Section 4.2) and AMSIPC (Active Modeling Based Stationary Increment Predictive Control, designed in Section 4.3), are tested in sudden mode-change, and are compared with each other on the ServoHeli-40 test-bed. To complete this mission, the following experimental process is designed:

1. Using large and step-like reference velocity, red line in Fig.6-8, input it to longitudinal loop, lateral loop and vertical loop;

2. Based on the same inputted reference velocity, using the 3 types of control method, GPC, SIPC and AMSIPC to actuate the helicopter to change flight mode quickly;

3. Record the data of position, velocity and reference speed for the 3 control loops, and obtain reference position by integrating the reference speed;

4. Compare errors of velocity and position tracking of GPC, SIPC and AMSIPC, executively, in this sudden mode-change flight.

GPC, SIPC and AMSIPC are all tested in the same flight conditions, and the comparison results are shown in Figs. 6-8. We use the identified parameters in Section 5.2 to build the nominal model, based on the model structure in Appendix A, and parameters’ selection in Appendix C for controllers

It can be seen that, when the helicopter increases its longitudinal velocity and changes flight mode from hovering to cruising, GPC (brown line) has a steady velocity error and increasing position error because of the model errors. SIPC (blue line) has a smaller velocity error because it uses increment model to reject the influence of the changing operation point and dynamics’ slow change during the flight. The prediction is unbiased and obtains better tracking performance, which is verified by Theorem. However, the increment model may enlarge the model errors due to the uncertain parameters and sensor/process noises, resulting in the oscillations in the constant velocity period (clearly seen in Fig.6&7) because the error of its prediction is only unbiased, but not minimum variance. While for AMSIPC (green line), because the model error, which makes the predictive process non-minimum variance, has

Figure 6.

Longitudinal tracking results: (a) velocity; (b) position error (<50s hovering, >50s cruising)

Figure 7.

Lateral tracking results: (a) velocity; (b) position error (25s~80s cruising, others hovering)

Figure 8.

Vertical tracking results: (a) velocity; (b) position error (<5s hovering; >5s cruising)

been online estimated by the ASMF and compensated by the strategy in section 4.3, the proposed AMSIPC successfully reduces velocity oscillations and tracking errors together.

6. Conclusion

An active model based predictive control scheme was proposed in this paper to compensate model error due to flight mode change and model uncertainties, and realize full flight envelope control without multi-mode models and mode-dependent controls.

The ASMF was adopted as an active modeling technique to online estimate the error between reference model and real dynamics. Experimental results have demonstrated that the ASMF successfully estimated the model error even though it is both helicopter dynamics and flight-state dependent.In order to overcome the aerodynamics time-delay, also with the active estimation for optimal compensation, an active modeling based stationary increment predictive controller was designed and analyzed.

The proposed control scheme was implemented on our developed ServoHeli-40 unmanned helicopter. Experimental results have demonstrated clear improvements over the normal GPC without active modeling enhancement when sudden mode-change happens.

It should be noted that, at present, we have only tested the control scheme with respect to the flight mode change from hovering to cruising, and vice versa. Further mode change conditions will be flight-tested in near future.

7. Appendix

A. Helicopter dynamics

A helicopter in flight is free to simultaneously rotate and translate in six degrees of freedom. Fig. 9 shows the helicopter variables in a body-fixed frame with origin at the vehicle’s center of gravity.

Figure 9.

Helicopter with its body-fixed reference frame

Ref.[18] developed a semi-decoupled model for small-size helicopter, i.e.,

$$\left(\begin{array}{c}\delta \dot{u}\\ \delta \dot{q}\\ \delta \dot{\theta }\\ \dot{a}\\ \dot{c}\end{array}\right)=\left(\begin{array}{ccccc}{X}_u& 0& -g& X_a& 0\\ M_u& 0& 0& M_a& 0\\ 0& 1& 0& 0& 0\\ 0& -1& 0& -1/{\tau }_f& A_c/{\tau }_f\\ 0& -1& 0& 0& -1/{\tau }_f\end{array}\right)\left(\begin{array}{c}\delta u\\ \delta q\\ \delta \theta \\ a\\ c\end{array}\right)+\left(\begin{array}{cc}{X}_{lon}& X_{lat}\\ M_{lon}& M_{lat}\\ 0& 0\\ A_{lon}& A_{lat}\\ C_{lon}& C_{lat}\end{array}\right)\left(\begin{array}{c}{\delta }_{lon}\\ {\delta }_{lat}\end{array}\right)$$E45

$$\{\begin{array}{l}{\dot{X}}_{lon}=A_{lon}\delta X_{lon}+B_{lon}\delta u_{lon}\\ y_{lon}=\left(\begin{array}{cc}{I}_{3\times 3}& 0_{3\times 2}\end{array}\right)\delta X_{lon}=C_{lon}\delta X_{lon}\end{array}$$E46

$$\left(\begin{array}{c}\delta \dot{v}\\ \delta \dot{p}\\ \delta \dot{\phi }\\ \dot{b}\\ \dot{d}\end{array}\right)=\delta {\dot{X}}_{lat}=\left(\begin{array}{ccccc}{Y}_u& 0& g& Y_a& 0\\ L_u& 0& 0& L_a& 0\\ 0& 1& 0& 0& 0\\ 0& -1& 0& -1/{\tau }_f& B_d/{\tau }_f\\ 0& -1& 0& 0& -1/{\tau }_f\end{array}\right)\left(\begin{array}{c}\delta v\\ \delta p\\ \delta \phi \\ b\\ d\end{array}\right)+\left(\begin{array}{cc}{Y}_{lon}& Y_{lat}\\ L_{lon}& L_{lat}\\ 0& 0\\ B_{lon}& B_{lat}\\ D_{lon}& D_{lat}\end{array}\right)\left(\begin{array}{c}{\delta }_{lon}\\ {\delta }_{lat}\end{array}\right)$$E47

$$\{\begin{array}{l}{\dot{X}}_{lat}=A_{lat}\delta X_{lat}+B_{lat}\delta u_{lat}\\ y_{lat}=\left(\begin{array}{cc}{I}_{3\times 3}& 0_{3\times 2}\end{array}\right)\delta X_{lat}=C_{lat}\delta X_{lat}\end{array}$$E48

$$\left(\begin{array}{c}\delta \dot{w}\\ \delta \dot{r}\\ \delta {\dot{r}}_{fb}\end{array}\right)=\delta {\dot{X}}_{yaw-heave}=\left(\begin{array}{ccc}{Z}_w& Z_r& 0\\ N_w& N_r& -N_{ped}\\ 0& K_r& -K_{rfb}\end{array}\right)\left(\begin{array}{c}\delta w\\ \delta r\\ \delta r_{fb}\end{array}\right)+\left(\begin{array}{cc}{Z}_{ped}& Z_{col}\\ N_{ped}& N_{col}\\ 0& 0\end{array}\right)\left(\begin{array}{c}{\delta }_{ped}\\ {\delta }_{col}\end{array}\right)$$E49

$$\{\begin{array}{l}{\dot{X}}_{yaw-heave}=A_{yaw-heave}\delta X_{yaw-heave}+B_{yaw-heave}\delta u_{yaw-heave}\\ y_{yaw-heave}=\left(\begin{array}{cc}{I}_{2\times 2}& 0_{2\times 1}\end{array}\right)\delta X_{yaw-heave}=C_{yaw-heave}\delta X_{yaw-heave}\end{array}$$E50

where δu, δv, δw are longitudinal, lateral and vertical velocity, δp, δq, δr are roll, pitch and yaw angle rates, δφ and δθ are the angles of roll and pitch, respectively, a and b are the first harmonic flapping angle of main rotor, c and d are the first harmonic flapping angle of stabilizer bar, $$\delta r_{fb}$$is the feedback control value of the angular rate gyro, $${\delta }_{lat}$$is the lateral control input, $${\delta }_{lon}$$is the longitudinal control input, $${\delta }_{ped}$$is the yawing control input, and $${\delta }_{col}$$ is the vertical control input. All the symbols except gravity acceleration g in$$A_{lon}$$, $$A_{lat}$$, $$A_{yaw-heave}$$, $$B_{lon}$$,$$B_{lat}$$ and $$B_{yaw-heave}$$ are unknown parameters to be identified. Thus, all of the states and control inputs in (33), (34) and (35) are physically meaningful and defined in body-axis.

B. Proof for the predictive theorem

Proof:

Assume the real dynamics is described as:

$$X_{t+1}=A_{dr}{X}_t+B_{dr}{U}_{t-k}+W_t$$E51

which is different from the reference model of Eq. (11). In Eq. (36), $$X_t$$is system state, $$A_{dr}$$is the system matrix, $$B_{dr}$$is the control matrix, $$U_t$$is control input, $$W_t$$is process noise. The one-step prediction, according to Eq. (36), can be obtained by Eq. (13-14),

$$\begin{array}{l}{\widehat{X}}_{t|t+1}=X_t+A_d\Delta X_t+B_d\Delta U_{t-k}\\ =A_{dr}{X}_{t-1}+B_{dr}{U}_{t-1-k}+W_{t-1}\\ +A_d\Delta X_t+B_d\Delta U_{t-k}\end{array}$$E52

And

$$\begin{array}{l}E\left\{X_{t+1}-{\widehat{X}}_{t+1|t}\right\}\\ =E\{A_{dr}{X}_t+B_{dr}{U}_{t-k}+W_t\\ -(A_{dr}{X}_{t-1}+B_{dr}{U}_{t-1-k}+W_{t-1}+A_d\Delta X_t+B_d\Delta U_{t-k})\}\\ =E\{(A_{dr}-A_d)\Delta X_t+(B_{dr}-B_d)\Delta U_{t-k}+\Delta W_t\}\end{array}$$E53

According to condition 1) and 2), prediction is bounded, then,

$$\Vert X_{t+1}-{\widehat{X}}_{t+1|t}\Vert <+\infty$$E54

and, when the system of Eq. (36) works around a working point in steady state, the mean value of control inputs and states should be constant, so we can obtain:

$$\begin{array}{l}E\left\{X_{t+1}-{\widehat{X}}_{t+1|t}\right\}\\ =(A_{dr}-A_d)E\left\{\Delta X_t\right\}+(B_{dr}-B_d)E\left\{\Delta U_{t-k}\right\}+E\left\{\Delta W_t\right\}\\ =(A_{dr}-A_d)•0+(B_{dr}-B_d)•0+0=0\end{array}$$E55

Eq. (39) indicates that the one step prediction of Eq. (37) is unbiased.

Assuming that prediction at time i-1 is unbiased, i.e..

$$E\{X_{t+i-1}-{\widehat{X}}_{t+i-1|t}\}=0$$E56

for the prediction at time i, there is

$$\begin{array}{l}E\{X_{t+i}-{\widehat{X}}_{t+i|t}\}\\ =E\{A_{dr}{X}_{t+i-1}+B_{dr}{U}_{t+i-1-k}+W_{t+i-1}\\ -({\widehat{X}}_{t+i-1|t}+A_d\widehat{\Delta }{X}_{t+i-1|t}+B_d\Delta U_{t+i-1-k})\}\\ =E\{A_{dr}{X}_{t+i-1}+B_{dr}{U}_{t+i-1-k}+W_{t+i-1}-X_{t+i-1}\\ +(X_{t+i-1}-{\widehat{X}}_{t+i-1|t})-W_{t+i-2}\\ -A_d\Delta {\widehat{X}}_{t+i-1|t}-B_d\Delta U_{t+i-1-k}\}\\ =E\{A_{dr}\Delta X_{t+i-1}-A_d\widehat{\Delta }{X}_{t+i-1|t}+\\ (B_{dr}-B_d)\Delta U_{t+i-1-k}+\Delta W_{t+i-1}\}\\ =(A_{dr}-A_d)E\left\{\Delta X_{t+i-1}\right\}+\\ (B_{dr}-B_d)E\left\{\Delta U_{t+i-1-k}\right\}+E\left\{\Delta W_{t+i-1}\right\}\\ =(A_{dr}-A_d)•0+(B_{dr}-B_d)•0+0=0\end{array}$$E57

Therefore, the prediction at time i is also unbiased.

C. Parameters’ selection for estimate and control in flight experiment

1. For Modeling

The identification results for hovering dynamics are listed in Tab.D-1.

Longitudinal Loop		Lateral Loop		Vertical Loop
Para.	Val.	Para.	Val.	Para.	Val.
Xu	0.2446	Yv	-0.0577	Zw	1.666
Xa	-4.962	Yb	9.812	Zr	-3.784
Xlat	-0.0686	Ylat	-1.823	Zped	2.304
Xlon	0.0896	Ylon	2.191	Zcol	-11.11
Mu	-1.258	Lv	15.84	Yaw Loop
Ma	46.06	Lb	126.6	Para.	Val.
Mlat	-0.6269	Llat	-4.875	Nw	-0.027
Mlon	3.394	Llon	28.64	Nr	-1.087
Ac	0.1628	Bd	-1.654	Nrfb	-1.845
Alat	-0.0178	Blat	0.04732	Nped	1.845
Alon	-0.2585	Blon	-9.288	Ncol	-0.972
Clat	2.238	Dlat	-0.7798	Kr	-0.040
Clon	-4.144	Dlon	-5.726	Krfb	-2.174
tf	0.5026	ts	0.5054

Table 2.

The parameters of hovering model

1. For ASMF

$$Q=\left(\begin{array}{cc}0.01I_{13\times 13}& 0_{13\times 13}\\ 0_{13\times 13}& 0.1I_{13\times 13}\end{array}\right)$$,

$$R=0.01I_{8\times 8}$$E58

where $$I_{m\times m}$$ is the m×m unit matrix and $$0_{m\times n}$$ is the m×n zero matrix.

1. For GPC

$p=10$, $\gamma =2.32I_{40\times 40}$,

$$k=10$$E59

1. For SIPC

$p=10$, $\gamma =2.32I_{4\times 4}$,

$$\alpha =0.99I_{8\times 8}$$E60

$W=I_{80\times 80}$, $k=10$,

$$\beta =0.8I_{4\times 4}$$E61

1. For AMSIPC

$p=10$, $\gamma =2.32I_{4\times 4}$,

$$\alpha =0.99I_{8\times 8}$$E62

$W=I_{80\times 80}$, $k=10$, $\beta =0.8I_{4\times 4}$,

$$H=I_{13\times 13}$$E63