Quadrotor stabilization by Fuzzy Kalman Filter

Abstract

In this work, state vector estimation by means of the Fuzzy Kalman Filter (FKF) is used to generate a control signal that stabilizes an unmanned quadrotor aircraft. The framework for fuzzy Kalman Filter methodology has been successfully developed, and in this sense, the FKF is implemented and compared with Kalman Filter (KF) and extended Kalman Filter (EKF). It will be proved that the fuzzy version gives some advantages such as a smaller processing time and a smaller Mean Squared Error (MSE). Finally, these results are shown in graphics and tables.

Keywords

Kalman filtering quadrotor fuzzy systems stabilization control systems

1 Introduction

The are some interesting investigations about filters [1 –4]. From the mentioned methodologies, the Kalman Filter (KF) has become one of the most useful techniques for estimating the unmeasured variables even in presence of noise. An extension of KF was developed for nonlinear systems and was named Extended Kalman Filter (EKF) which is frequently called the optimal observer. It is used in a number of estimation and tracking applications as satellite navigation, robotics, missile trajectory estimation, radar control or spacecraft stabilization [5, 6]. KF and EKF have been also successfully applied in the context of evolving fuzzy systems (for updating consequent parameters over the time through the closely related to the real-beam forward-looking weighted least squares “RFWLS” algorithm) [7, 8]. Its methodology framework has been widely studied, developed and validated by many researchers over the years.

In [9, 10], the called argument Kalman filter is introduced in order to reduce the processing time of the extended Kalman filter for the learning of a neural network where the contribution was based on the substitution of the Jacobian by its argument.

There is an issue with the above mentioned filters for the reduction of the processing time, it is that their designs are rarely focused on applications, i.e., in most of the cases, the filters are employed in synthetic examples. Then, there is a necessity to design novel filters for the reduction of the processing time and to be focused on applications.

In this work, a new version of Kalman Filter called Fuzzy Kalman Filter (FKF) is presented based on Kalman Filter and Takagi-Sugeno Fuzzy models. It is an optimal fuzzy state estimator recursive algorithm for nonlinear systems. Developed with the goal to reduce the mean square error and also the processing time, the results show better performance compared to classic algorithms used for estimating the unmeasurable state of dynamical systems, yielding a reduction of the processing time.

In present work, FKF is used to estimate the state vector in an optimal way to generate a control signal to achieve the stabilization of an unmanned aerial vehicle (UAV) by means of measurement signals with white noise. The UAV refers to a quadrotor that is a four-axis multirotor aircraft whose elevation is generated by four rotors oriented vertically. Quadrotors have been highly researched and actually there exist many applications approaches, models and control systems as in [5 , 12].

Sensors used in a quadrotor such as gyroscopes and accelerometers provide signals for control altitude process and these signals are regularly added with noise, in greater proportion in gyroscopes. This represents a complication for the control algorithm because it causes the data to have precision errors or bad response and the algorithm takes more processing time in real flights, in which the sampling speed matters. This is why the Kalman Filter is a tool widely used in UAV’s control algorithms.

2 Fuzzy Kalman filter

Basically, the mathematics involved in the Kalman Filter is derived based on the least squares criterion with an optimal gain and consists of a set of equations statically optimal with respect to some quadratic function of estimation error, while the measurement noise is decreased [13].

2.1 Kalman filter

Given a linear discrete-time system with additive noise $x_{k + 1} = A_{k} x_{k} + B_{k} u_{k} + Γ_{k} w_{k}$ (1) $y_{k} = G_{k} x_{k} + v_{k}$ (2) Where u_k of dimension m is the feedback controller, x_k of dimension n represents state vector of the dynamics, and y_k of dimension q is the measured signal, A_k and G_k are constant transition matrix and B_k are the input matrix, w_k of dimension p and v_k of dimension q are the process and measurement Gaussian noises, respectively, they are assumed as completely random with covariance matrices Q and R, respectively, and Γ_k is a known constant matrix with n × p dimension. More precisely, w_k ≈ N (0, Q) , and v_k ≈ N (0, R).

The Kalman Filter algorithm has two processes stages: the state estimation and update, where it is defined by the next set of linear equations:

state estimation ${\hat{x}}_{k}^{-} = A_{k - 1} {\hat{x}}_{k} + B_{k - 1} u_{k - 1}$ (3) ${\hat{P}}_{k}^{-} = A_{k - 1} P_{k} A_{k - 1}^{T} + Q_{k - 1}$ (4) update $K_{k} = {\hat{P}}_{k}^{-} G_{k}^{T} {(G_{k} {\hat{P}}_{k}^{-} G_{k}^{T} + R_{k})}^{- 1}$ (5) ${\hat{x}}_{k} = {\hat{x}}_{k}^{-} + K_{k} (y_{k} - G_{k} {\hat{x}}_{k}^{-})$ (6) $P_{k} = (I + K_{k} G_{k}) {\hat{P}}_{k}^{-}$ (7)

Where K_k is the optimal Kalman gain that guarantees the zero convergence of the error, ${\hat{x}}_{k}^{-}$ is the estimation of the state vector from the previous estimated state starting with ${\hat{x}}_{0}$ with zero mean, ${\hat{x}}_{k}$ is the correction of the estimate through the Kalman gain and error, likewise, ${\hat{P}}_{k}^{-}$ is the estimation of the error covariance matrix and P_k is the corrected value, it is considered an initial covariance matrix P₀. The corrected values are used again for subsequent operations and the process is recursive.

2.2 Extended Kalman filter

Most of the dynamic systems are nonlinear, then the need to extend the results obtained with the Kalman filter to the nonlinear region, then the Extended Kalman Filter (EKF) is developed. So, a nonlinear discrete-time system with perturbations is considered $x_{k + 1} = f (x_{k}, u_{k}) + Γ_{k} w_{k}$ (8) $y_{k} = g (x_{k}) + v_{k}$ (9)

With w_k and v_k, as above. It is assumed that the functions f (x_k, u_k) and g (x_k) are twice continuously differentiable, so they can be approximated by a Taylor Series expansion about ${\hat{x}}_{k}^{-}$ as follows $f (x_{k}, u_{k}) \approx f ({\hat{x}}_{k}^{-}) + J_{f} ({\hat{x}}_{k}^{-}) (x_{k} - {\hat{x}}_{k}^{-}) + H . O . T .$ (10)

Where J _f is the Jacobian of f(·) and H . O . T . are considered negligible. The process is performed similarly for g(·). Using (10), the approximation error $e_{k} = x_{k} - {\hat{x}}_{k}^{-}$ and the expected value of f and g by previous measurement (see [13]), then the Extended Kalman Filter estimation and update equations results.

state estimation ${\hat{x}}_{k}^{-} = f ({\hat{x}}_{k}, u_{k})$ (11) ${\hat{P}}_{k}^{-} = J_{f} ({\hat{x}}_{k}^{-}) P_{k} J_{f}^{T} ({\hat{x}}_{k}^{-}) + Q_{k - 1}$ (12) update $K_{k} = {\hat{P}}_{k}^{-} J_{g}^{T} ({\hat{x}}_{k}^{-}) {(J_{g} ({\hat{x}}_{k}^{-}) {\hat{P}}_{k}^{-} J_{g}^{T} ({\hat{x}}_{k}^{-}) + R_{k})}^{- 1}$ (13) ${\hat{x}}_{k} = {\hat{x}}_{k}^{-} + K_{k} (y_{k} - J_{g} ({\hat{x}}_{k}^{-}) {\hat{x}}_{k}^{-})$ (14) $P_{k} = (I + K_{k} J_{g} ({\hat{x}}_{k}^{-})) {\hat{P}}_{k}^{-}$ (15)

Just as in the KF, the EKF is a recursive algorithm, but in this case, it is necessary to linearizing by calculate the Jacobian matrix for each iteration.

2.3 Fuzzy Kalman filter

A new methodology of Kalman Filter based on the TS fuzzy modeling is developed in this work and named Fuzzy Kalman Filter which use the advantages of linear models in nonlinear systems as an alternative to the classic algorithms.

Fuzzy Kalman Filter algorithm is developed based on the assumption that is possible to describe the nonlinear model as a set of local linear models and its input-output TS fuzzy discrete-time model is described by

Plant rule i

If z₁ (k) is M_i1 and ... z_p (k) is M_ip $Then \sum_{i} {\begin{matrix} x_{k + 1} = A_{i} x_{k} + B_{i} u_{k} \\ y_{k} = C_{i} x_{k} \end{matrix}$ (16)

Where, as before, x_k, u_k, and y_k are the state, input, and output vectors, respectively. Besides, M_ip are the corresponding fuzzy sets that weighs each subsystem, [z_i (k) , …, z_p (k)] are the measurable premise variables or output state vector, i is the number of subsystems and matrices A_i, B_i, C_i, are obtained by linearizing the nonlinear system around some suitable operation point [14].

The complete fuzzy model is obtained by a singleton fuzzifier, product interference and center gravity defuzzifier [14], and adding the noises. Thus, any given f (x_k, u_k) is then defined by the next fuzzy system $x_{k + 1} = \sum_{i}^{r} h_{i} (z_{k}) A_{i} x_{k} + \sum_{i}^{r} h_{i} (z_{k}) B_{i} u_{k} + Γ_{k} w_{k}$ (17) $y_{k} = \sum_{i}^{r} h_{i} (z_{k}) C_{i} x_{k} + v_{k}$ (18)

Where each subsystem corresponds to linearization around some operation points, the factor h_i is the normalized weight for each rule calculated from the membership value of z_i (k) in M_ip, which satisfies $\sum_{i}^{r} h_{k} (z_{k}) = 1, h (z_{k}) ⩾ 0$ (19)

A nonlinear dynamic system in discrete time with disturbances is described by $x_{k + 1} = f (x_{k}, u_{k}) + w_{k}$ (20) $y_{k} = g (x_{k}) + v_{k}$ (21)

Fuzzy Kalman Filter for nonlinear system described by TS Fuzzy modeling as a good approximation of nonlinear system uses the next approximation. $f (x_{k}, u_{k}) \approx \sum_{i}^{r} h_{i} (z_{k}) A_{i} x_{k} + \sum_{i}^{r} h_{i} (z_{k}) B_{i} u_{k}$ (22) $g (x_{k}) \approx \sum_{i}^{r} h_{i} (z_{k}) C_{i} x_{k}$ (23)

This allows writing the fuzzy filter system matrices as follows ${\bar{A}}_{k} = \sum_{i = 1}^{r} h_{i} (z_{k}) A_{i}$ (24) ${\bar{B}}_{k} = \sum_{i = 1}^{r} h_{i} (z_{k}) B_{i}$ (25) ${\bar{C}}_{k} = \sum_{i = 1}^{r} h_{i} (z_{k}) C_{i}$ (26)

Consequently, the TS fuzzy model (17)–(18) can be rewritten as $x_{k + 1} = {\bar{A}}_{k} x_{k} + {\bar{B}}_{k} u_{k} + w_{k}$ (27) $y_{k} = {\bar{C}}_{k} x_{k} + v_{k}$ (28)

The so-called discrete-time Fuzzy Kalman Filter is written as one original filter in two stages, prediction and updating, and it is defined by the following.

state estimation ${\hat{x}}_{k}^{-} = {\bar{A}}_{k - 1} {\hat{x}}_{k - 1} + {\bar{B}}_{k - 1} u_{k - 1}$ (29) $P_{k}^{-} = {\bar{A}}_{k - 1} P_{k - 1} {\bar{A}}_{k - 1}^{T} + Q_{k - 1}$ (30)

update $K_{k} = P_{k}^{-} {\bar{C}}_{k}^{T} {({\bar{C}}_{k} P_{k}^{-} {\bar{C}}_{k}^{T} + R_{k})}^{T}$ (31) ${\hat{x}}_{k} = {\hat{x}}_{k}^{-} + K_{k} (y_{k} - {\bar{C}}_{k} {\hat{x}}_{k}^{-})$ (32) $P_{k} = (I + K_{k} {\bar{C}}_{k}) P_{k}^{-}$ (33)

As in the original formulation of linear KF, it is assumed that the initial state vector is Gaussian with mean $\bar{x} = 0,$ and initial covariance matrix P₀.

Notice that the process is very similar to the classical KF, but the filter is now working on a nonlinear system described by a TS fuzzy model. Besides, the proposed methodology does not require linearization for each iteration unlike the EKF. Instead, the matrices are obtained at any instant by (24)–(26). This advantage will be illustrated in a later section.

3 Quadrotor model

Several mathematical models of the quadrotor have been developed in previous works because of the mechanical simplicity and due to assumption that it is a rigid body [11, 15].

In this section, a model of quadrotor previously developed is presented. Reference systems, as well as all mathematics and modeling strategies, are given in [5], where the goal was the sinusoidal reference tracking in height generated by an exosystem within an ideal noise-free situation and in [16], to compare EKF and KF with an estimator named James-Stein State Estimator.

The mathematical model described below is the one developed in [5], with six degrees of freedom and two frameworks, inertial frame of reference (ground) and a non-inertial frame reference (body) which moves with the quadrotor, it is assumed that is a rigid body in order to be described using Newton-Euler equations and proper to the aircraft with “x” type flight configuration [11, 17].

The aircraft coordinate frame axes are described relative to the inertial frame of reference (E), the general position vector ξ is defined as $ξ = {[Γ^{E} Θ^{E}]}^{T} = {[X Y Z ϕ θ ψ]}^{T}$ (34)

Where Γ^E is the linear position vector and Θ^E is the angular position vector of the vehicle with ϕ, θ, ψ, the roll, pitch, and yaw angles with respect to the inertial frame.

On the other hand, on a non-inertial frame (B), the linear velocity (V^B) of the body is defined with components (u v w), and (p q r) represents the components of angular velocity vector (ω^B), then the general velocity is given by $v = {[V^{B} ω^{B}]}^{T} = {[u υ w p q r]}^{T}$ (35)

The generalized matrix of rotation and translation coordinates that describes the dynamics of aircraft relative to ground is $J_{Θ} = [\begin{matrix} C & 0_{3 X 3} \\ 0_{3 X 3} & T_{Θ} \end{matrix}]$ (36)

Where, C is a rotation matrix and T_Θ is a translation matrix. Consequently, the kinematics of a rigid body with six degree of freedom is given by equation $\dot{ξ} = J_{Θ} v$ (37)

Based on the axioms of Euler and the second law of Newton the equations in matrix form result $[\begin{matrix} m \hat{I} & 0_{3 X 3} \\ 0_{3 X 3} & \hat{I} \end{matrix}] [\begin{matrix} {\dot{V}}^{B} \\ {\dot{ω}}^{B} \end{matrix}] + [\begin{matrix} ω^{B} \times ({mV}^{B}) \\ ω^{B} \times (\hat{I} ω^{B}) \end{matrix}] = [\begin{matrix} F^{B} \\ τ^{B} \end{matrix}]$ (38)

Where $\hat{I}$ represents the diagonal matrix of inertia. Three kind of force contributions are considered to obtain the full measured of force, the inertial matrix of the aircraft (M_B), the centripetal Coriolis matrix (C_B), the gravitational vector represented by G_B, and results $M_{B} \dot{v} + C_{B} (v) v = G_{B} (ξ) + O_{B} (v) Ω + E_{B} Ω^{2}$ (39)

Where O_B is the gyroscopic propeller matrix, Ω the full velocity of rotors and E_B is the movement matrix. This equation describes the movement of the system with respect to the non-inertial reference.

In order to get a description with respect to inertial reference (ground), using transformations matrices between both reference and introducing the forces and moments defined in (3), the dynamics of the vehicle are described by the following equations $\begin{matrix} x = (sen ψ sen ϕ + \cos ψ sen θ \cos ϕ) \frac{U_{1}}{m} \\ y = (- \cos ψ sen ϕ + sen ψ sen θ \cos ϕ) \frac{U_{1}}{m} \\ z = - g + (\cos ϕ + \cos ϕ) \frac{U_{1}}{m} \\ ϕ = \dot{θ} ψ (\frac{I_{YY} - I_{ZZ}}{I_{XX}}) - \frac{J_{TP}}{I_{XX}} \dot{θ} Ω + \frac{{lU}_{2}}{I_{XX}} \\ θ = \dot{ϕ} ψ (\frac{I_{ZZ} - I_{XX}}{I_{YY}}) + \frac{J_{TP}}{I_{YY}} \dot{ϕ} Ω + \frac{{lU}_{3}}{I_{YY}} \\ ψ = \dot{ϕ} θ (\frac{I_{XX} - I_{YY}}{I_{ZZ}}) + \frac{{lU}_{4}}{I_{ZZ}} \end{matrix}$ (40)

This group of equations includes the transformation of axes between reference frames through a matrix of rotations and defines the displacements (x, y, z) and angular (ϕ, θ, ψ) positions of the aircraft referenced to the ground system, Ω is the velocity vector [5, 13].

The full velocity vector Ω (rad/s) is represented by the algebraic sum of each of the speeds of the rotors $Ω = Ω_{1} - Ω_{2} + Ω_{3} - Ω_{4}$ (41)

The group of control inputs is described in function of velocity of rotors by $[\begin{matrix} U_{1} \\ U_{2} \\ U_{3} \\ U_{4} \end{matrix}] = [\begin{matrix} b (Ω_{1}^{2} + Ω_{2}^{2} + Ω_{3}^{2} + Ω_{4}^{2}) \\ b (Ω_{4}^{2} + Ω_{3}^{2} - Ω_{1}^{2} - Ω_{2}^{2}) \\ b (Ω_{2}^{2} + Ω_{3}^{2} - Ω_{1}^{2} - Ω_{4}^{2}) \\ d (Ω_{1}^{2} + Ω_{3}^{2} - Ω_{2}^{2} - Ω_{4}^{2}) \end{matrix}]$ (42)

The physical characteristics for the quadrotor used in this work are presented below $\begin{matrix} b = 54.2 \times 10^{- 6} N s^{2}, \\ d = 1.1 \times 10^{- 6} N m s^{2}, l = 0.24 m, \\ m = 1 kgm, \\ g = 9.81 m / s^{2}, \\ J_{tp} = 104 \times 10^{- 6} N m^{2} s^{2}, \\ I_{xx} = 8.1 \times 10^{- 3} N m^{2} s^{2}, \\ I_{yy} = 8.1 \times 10^{- 3} N m^{2} s^{2}, \\ I_{zz} = 14.2 \times 10^{- 3} N m^{2} s^{2} . \end{matrix}$

Where b is the push factor, d is the drag factor, l is the distance between the center of gravity of the quadcopter and the center of the rotors, m is the mass of quadrotor, and J_tp is the moment of inertia of the engines.

Based on the linear model, linearized by Jacobian at the equilibrium point x_e, it was verified that the system is controllable and with this result it is possible to obtain a feedback gain by Robust Pole Placement technique [5 , 18] where is used to stabilize the drone in attitude and the feedback controller results $u = ϒ (x_{k} - x_{op}) + Ω$ (43)

Where x_op is the altitude operation point vector, that was placed as altitude equal to 30 meters and with other position, rotations, velocities and accelerations values equal to zero [5]. ϒ is the state-feedback gain matrix obtained considering the desired eigenvalues $p = [0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91,$ 0.92, 0.93, 0.94, 0.95]

4 Numerical results

In order to compare the efficiency of the proposed methodology, three versions of Kalman Filter are used, Kalman Filter (KF), Extended Kalman Filter (EKF), and Fuzzy Kalman Filter (FKF), with the purpose of comparing efficiency in the altitude stabilization of the aircraft, each filter applied to obtain the state vector x_k that allows calculating the control signal with (43).

Consider the dynamic quadrotor model described in previous section to be stabilized, for this, it is possible to write it as a state vector [5, 16, 18] $\begin{matrix} {\dot{x}}_{1} = x_{2} \\ {\dot{x}}_{2} = \sin (x_{11}) \sin (x_{7}) + \frac{u_{1}}{m} [cos (x_{11}) \sin (x_{9}) cos (x_{7})] \\ {\dot{x}}_{3} = x_{4} \\ {\dot{x}}_{4} = - cos (x_{11}) \sin (x_{7}) + \frac{u_{1}}{m} [\sin (x_{11}) \sin (x_{9}) cos (x_{7})] \\ {\dot{x}}_{5} = x_{6} \\ {\dot{x}}_{6} = - g + \frac{u_{1}}{m} [cos (x_{9}) cos (x_{7})] \\ {\dot{x}}_{7} = x_{8} \\ {\dot{x}}_{8} = x_{10} x_{12} \frac{(I_{yy} - I_{zz})}{I_{xx}} - \frac{J_{tp}}{I_{xx}} Ω x_{10} + \frac{l u_{2}}{I_{xx}} \end{matrix}$ $\begin{matrix} {\dot{x}}_{9} = x_{10} \\ {\dot{x}}_{10} = x_{8} x_{12} \frac{(I_{zz} - I_{xx})}{I_{yy}} + \frac{J_{tp}}{I_{yy}} Ω x_{8} + \frac{l u_{3}}{I_{yy}} \\ {\dot{x}}_{11} = x_{12} \\ {\dot{x}}_{12} = x_{8} x_{10} \frac{(I_{xx} - I_{yy})}{I_{zz}} + \frac{u_{4}}{I_{zz}} \end{matrix}$

A discrete-time approximation is needed and is obtained by considering the first-order derivative, [5, 16], which is approximated by the following expression $\dot{x} (t) \approx \frac{x (t + T) - x (t)}{T}$ (44)

Where T is the sampling time, which is chosen sufficiently small and reordering terms $x (t + T) \approx x (t) + T \dot{x} (t)$ (45)

Then, following the Euler’s methodology, an approximation of the system in discrete-time is obtained $x_{k + 1} = x_{k} + T {\dot{x}}_{k}$ (46) with x_k = [x₁, …, x₁₂] ^T, ${\dot{x}}_{k} = {[{\dot{x}}_{1}, {\dot{x}}_{2}, \dots, {\dot{x}}_{12}]}^{T},$ and sampling time is chosen adequately small, in this case T = 0.025s.

The nonlinear system with dynamic and measurement noise is described as follows $x_{k + 1} = f (x_{k}, u_{k}) + Γ_{k} w_{k}$ (47) $y_{k} = g (x_{k}) + v_{k}$ (48)

Where f is the quadrotor function described by the equations in (40) and $f (x_{k}, u_{k}) = x_{k} + T {\dot{x}}_{k}, Γ_{k} = I_{12 \times 12}$ , where I_12 ×12 is the identity matrix of dimension 12×12 because there is one noise for each state, and r is the standard deviation of the measurement noise [16]. Besides w_k and v_k are considered as Gaussian noises. Therefore, the TS fuzzy model to approximate (47)–(48) is

Plant rule i:

If x₁ (k) is M_i1 and ... x₁₂ (k) is M_i2 $x_{k + 1} = \sum_{i = 1}^{2} (h_{i} (x_{k}) A_{i} x_{k} + h_{i} (x_{k}) B_{i} u_{k}) + Γ_{k} w_{k}$ (49) $y_{k} = \sum_{i = 1}^{2} h_{i} (x_{k}) C_{i} x_{k} + v_{k}$ (50)

With x_k, u_k, y_k, x_k, w_k, and v_k as above, where the matrices considered to define the two-rule fuzzy system are computed as follows: $A_{i} = {\frac{\partial f}{\partial x} |}_{x_{op, i}}$ , $B_{i} = {\frac{\partial f}{\partial u} |}_{x_{op, i}}$ , $C_{i} = {\frac{\partial g}{\partial x} |}_{x_{op, i}}$ , and Γ_k = I_12 ×12. Besides, for this example, the state variable x₇, which refers to the pitch angle in radians is used to define the membership function, as follows $h_{1} (x_{k}) = (\frac{x_{7 L 2} - x_{7}}{dif}), h_{2} (x_{k}) = (\frac{x_{7} - x_{7 L 1}}{dif})$

Where $x_{7 L 1} = \frac{- 10 π}{180}$ and $x_{7 L 2} = \frac{10 π}{180}$ are the linearization points and dif = x_7L2 - x_7L1.

Then C₁ = C₂ = 1_1 ×12, where 1_1 ×12 is a vector with 1 in all the components, A₁ and A₂ are upper matrices with ones in the main diagonal and zero in the other components, except the following ones

For the matrix A₁: $\begin{matrix} A_{(1, 2)} = 0.0250 & A_{(3, 4)} = 0.0250 & A_{(5, 6)} = 0.0250 \\ A_{(4, 7)} = - 0.2452 & A_{(6, 7)} = 0.0432 & A_{(7, 8)} = 0.0250 \\ A_{(2, 9)} = 0.2452 & A_{(9, 10)} = 0.0250 & A_{(2, 11)} = - 0.0432 \\ A_{(12, 11)} = 0.0250 \end{matrix}$

And for the matrix A₂: $\begin{matrix} A_{(1, 2)} = 0.0250 & A_{(3, 4)} = 0.0250 & A_{(5, 6)} = 0.0250 \\ A_{(4, 7)} = - 0.2452 & A_{(6, 7)} = - 0.0432 & A_{(7, 8)} = 0.0250 \\ A_{(2, 9)} = 0.2452 & A_{(9, 10)} = 0.0250 & A_{(2, 11)} = 0.0432 \\ A_{(12, 11)} = 0.0250 \end{matrix}$

While matrices B₁ and B₂ turn out to be $B_{1} = [\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \begin{matrix} 0.0001 \\ 0 \end{matrix} & \begin{matrix} 0.0001 \\ 0 \end{matrix} & \begin{matrix} 0.0001 \\ 0 \end{matrix} & \begin{matrix} 0.0001 \\ 0 \end{matrix} \\ 0.0006 & 0.0006 & 0.0006 & 0.0006 \\ 0 & 0 & 0 & 0 \\ - 0.0172 & - 0.0172 & 0.0172 & 0.0172 \\ 0 & 0 & 0 & 0 \\ - 0.0172 & 0.0172 & 0.0172 & - 0.0172 \\ \begin{matrix} 0 \\ 0.0008 \end{matrix} & \begin{matrix} 0 \\ - 0.0008 \end{matrix} & \begin{matrix} 0 \\ 0.0008 \end{matrix} & \begin{matrix} 0 \\ - 0.0008 \end{matrix} \end{matrix}]$ $B_{2} = [\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \begin{matrix} - 0.0001 \\ 0 \end{matrix} & \begin{matrix} - 0.0001 \\ 0 \end{matrix} & \begin{matrix} - 0.0001 \\ 0 \end{matrix} & \begin{matrix} - 0.0001 \\ 0 \end{matrix} \\ 0.0006 & 0.0006 & 0.0006 & 0.0006 \\ 0 & 0 & 0 & 0 \\ - 0.0172 & - 0.0172 & 0.0172 & 0.0172 \\ 0 & 0 & 0 & 0 \\ - 0.0172 & 0.0172 & 0.0172 & - 0.0172 \\ \begin{matrix} 0 \\ 0.0008 \end{matrix} & \begin{matrix} 0 \\ - 0.0008 \end{matrix} & \begin{matrix} 0 \\ 0.0008 \end{matrix} & \begin{matrix} 0 \\ - 0.0008 \end{matrix} \end{matrix}]$

It is assumed that dynamic noise is very small for all process, and in order to test the filter’s ability to estimate states from noisy signals, two measures of standard deviation of the measurement noise are proposed as $r = [\begin{matrix} 0.15 & 0.01 \end{matrix}]$ . The objective is to achieve stabilization at altitude, so the state that interest is the state x₅ and random initial condition is taken around this operation point. It is worth mentioning that being a random system, each Kalman Filter is applied on many times, showing an average of results.

In order to clarify, the Fig. 1 is a representative measurement signal with noise for the altitude state (x₅) that is obtained from the quadrotor through the communication signal and is with which the filter works to estimate the states.

Fig. 1

Representative Measurement noise signal.

5 The first example

First, the nonlinear nature of the system suggest the use of the Extended Kalman Filter, then this algorithm is applied by taken the standard deviation of measurement noise is equal to r = 0.15 and its results are shown in Fig. 2, and in order to provide more information, some numerical results are shown in Table 1.

Fig. 2

Numerical results for EKF.

Table 1

Numerical results for EKF, it r = 0.15

Processing time	MSE Average
t = 0.8688sec.	8.5584e-05

In this case, during the stabilization problem, only about 15 percent of the tests were successful, in the other cases, the system cannot be well approximate by the Extended Kalman Filter, because of the nature of the nonlinear model. So, the MSE average sown in Table 2 is the MSE average of available cases.

Table 2

Numerical results for KF, it r = 0.15

Processing time	MSE Average
t = 0.2233sec.	0.9587

The next step is applying the KF. The same standard deviation of measurement noise is used. In this case, during the process, the approximation was possible in all tests. The stabilization by means of the KF is not quite well as shown in Fig. 3, so the MSE average obtained is greater than applying the EKF as shown in Table 2.

Fig. 3

Numerical results for KF.

Now, with the same parameters the Fuzzy Kalman Filter is applied during the stabilization process. As in the previous case, the approximation was successful in all tests, besides FKF provides a lower MSE and processing time than applying KF and EKF as shown in Table 3. The states and estimate signals are presented in Fig. 4.

Table 3

Numerical results for FKF, it r = 0.15

Processing time	MSE Average
t = 0.2466sec.	3.3158e-05

Fig. 4

Numerical results for FKF.

6 The second example

Now, a smaller standard deviation of measurement noise, equal to r = 0.01 is considered on the same aircraft system model. First, after using the EKF for estimation state vector it results are given in Fig. 5.

Fig. 5

Numerical results for EKF.

Again, the Kalman filter works in 30 percent of the cases, although the performance with a lower standard deviation of measurement noise improved, but cannot be obtained well approximate by the Extended Kalman Filter. Once again, the MSE average shown in Table 4 is for the successful cases.

Table 4

Numerical results for EKF, it r = 0.01

Processing time	MSE Average
T = 0.8487sec.	1.7446e-06

Now, the KF is applied under the same conditions. Again, the approximation was possible in all tests.

In this case we can notice that the stabilization is still not correct, even with less measurement noise as shown in Fig. 6 and in the numerical of the Table 5, this is because KF is a linear filter and the system is nonlinear.

Fig. 6

Numerical results for KF.

Table 5

Numerical results for KF, it r = 0.01

Processing time	MSE Average
t = 0.2209sec.	0.4861

Finally, the FKF is applied to the process with all the same parameters as in EKF and KF. As we can see in Table 6, the MSE average is quite small compared with EKF that provides a good approximation. The processing time has been also reduced due to the fact that during the estimation process it is not necessary to use linearization. The stabilization graph is presented in Fig. 7, which that the stabilization is well achieved.

Table 6

Numerical results for FKF, it r = 0.01

Processing time	MSE Average
t = 0.2550sec.	1.6910e-06

Fig. 7

Numerical results for FKF.

7 Conclusion

In this work, the Fuzzy Kalman Filter (FKF) has been proposed to stabilize a vertical component quadrotor. EKF presents complications for the estimation of states in most cases due to the nature of the system. Instead, using the FKF working with nonlinear systems, a good approximation is achieved as well as compared to standard KF, besides, a quite small processing time for the FKF. It would be interesting to apply the methodology directly on a real unmanned aircraft, and to record the results giving continuity to the investigation.

Footnotes

Acknowledgments

Authors want to express their deepest thank to CONACyT and Instituto Politécnico Nacional for the scholarships granted and research projects.

References

Shao

, Zhang

, Bo

and Park

, Multi-granulation rough filters and rough fuzzy filters in Pseudo-BCI algebras, Journal of Intelligent & Fuzzy Systems 34(6) (2018), 4377–4386.

, Zhang

, Shao

and Park

, The lattice generated by hesitant fuzzy filters in pseudo-BCI algebras, Journal of Intelligent & Fuzzy Systems 35(3) (2018), 3333–3345.

Wang

, Xin

and Wang

, Integral hoops and integral filters, Journal of Intelligent & Fuzzy Systems 35(1) (2018), 1055–1062.

Pazhanimuthu

and Ramesh

, Grid integration of renewable energy sources (RES) for power quality improvement using adaptive fuzzy logic controller based series hybrid active power filter (SHAPF), Journal of Intelligent & Fuzzy Systems 35(1) (2018), 749–766.

Ferdaus

M.M.

, Pratama

, Anavatti

S.G.

, Garratt

M.A.

and Pan

, Generic evolving self-organizing neuro-fuzzy control of bio-inspired unmanned aerial vehicles, IEEE Transaction on Fuzzy Systems (2019). DOI: 10.1109/TFUZZ.2019.2917808

, Yang

, Shi

, Pan

, Chen

and Sun

, Online recorded data-based composite neural control of strict-feedback systems with application to hypersonic flight dynamics, IEEE Transaction on Neural Networks and Learning Systems 29(8) (2018), 3839–3849.

Lughofer

, Evolving Fuzzy Systems – Methodologies, Advanced Concepts and Applications, Springer, Berlin Heidelberg, 2011.

Lughofer

, Evolving Fuzzy Systems – Fundamentals, Reliability, Interpretability and Useability, Applications. Editor P. Angelov, in Handbook of Computational Intelligence, World Scientific, pp. 67–135, New York, 2016.

Rubio

J.J.

, Lughofer

, Meda-Campaña

J.A.

, Paramo

L.A.

, Novoa

J.F.

and Pacheco

, Neural network updating via argument Kalman filter for modeling of Takagi-Sugeno fuzzy models, Journal of Intelligent & Fuzzy Systems 35(2) (2018), 2585–2596.

10.

Rubio

J.J.

, Lughofer

, Angelov

, Novoa

J.F.

and Meda-Campaña

J.A.

, A novel algorithm for the modeling of complex processes, Kybernetika 54(1) (2018), 79–95.

11.

Chovancová

, Fico

, Chovanec

L’.

and Hubinský

, Mathematical Modelling and Parameter Identification of Quadrotor (a survey), Modelling of Mechanical and Mechatronic Systems MMaMS 96 (2014), 172–181.

12.

Ferdaus

, Pratama

, Anavatti

S.G.

, Garratt

M.A.

and Lughoffer

, PAC: A novel self-adaptive neuro-fuzzy controller for micro aerial vehicles, Information Sciences 512 (2020), 481–505.

13.

Grewal

M.S.

and Andrews

A.P.

, Kalman Filtering: Theory and Practice Using MATLAB, Second Edition, New York: John Wiley & Sons, Inc., 2001.

14.

Tanaka

and Wang

H. O.

, Fuzzy Control Systems Design and Analysis: A Linear Matrix Inequality, New York: John Wiley & Sons, Inc., Apr 2004.

15.

Magnussen

Ø.

and Skjønhaug

K.E.

, Modeling, Design and Experimental Study for a Quadcopter System Construction, Kristiansand: Master’s Thesis, Department of Engineering Faculty of Technology and Science, 2011.

16.

Meda Campaña

J.A.

, On the estimation and control of nonlinear systems with parametric uncertainties and noisy outputs, IEEE Access 6 (2018), 31968–31973.

17.

Reinoso

, Minchala

L.I.

, Ortiz

J.P.

, Astudillo

and Verdugo

, Trajectory tracking of a quadrotor using sliding mode control, IEEE Latin America Transactions 14(5) (2016), 2157–2166.

18.

Rubio

J.J.

, Ochoa

, Mujica-Vargas

, Garcia

, Balcazar

, Elias

, Cruz

D.R.

, Juarez

C.F.

, Aguilar

and Novoa

J.F.

, Structure regulator for the perturbations attenuation in a quadrotor, IEEE Access 7(1) (2019), 138244–138252.