**Towards Simultaneous Actuator and Sensor Faults Estimation for a Class of Takagi-Sugeno Fuzzy**

#### **Marcin Pazera, Marcin Witczak, Norbert Kukurowski and Mariusz Buciakowski** *∗*

**Systems: A Twin-Rotor System Application**

Institute of Control and Computation Engineering, University of Zielona Góra, ul. Szafrana 2, 65-516 Zielona Góra, Poland; m.pazera@issi.uz.zgora.pl (M.P.); m.witczak@issi.uz.zgora.pl (M.W.); n.kukurowski@issi.uz.zgora.pl (N.K.)

**\*** Correspondence: m.buciakowski@issi.uz.zgora.pl

Received: 6 May 2020; Accepted: 17 June 2020; Published: 19 June 2020

**Abstract:** The paper is devoted to the problem of estimating simultaneously states, as well as actuator and sensor faults for Takagi–Sugeno systems. The proposed scheme is intended to cope with multiple sensor and actuator faults. To achieve such a goal, the original Takagi–Sugeno system is transformed into a descriptor one containing all state and fault variables within an extended state vector. Moreover, to facilitate the overall design procedure an auxiliary fault vector is introduced. In comparison to the approaches proposed in the literature, a usual restrictive assumption concerning fixed fault rate of change is removed. Finally, the robust convergence of the whole observer is guaranteed by the so-called quadratic boundedness approach which assumes that process and measurement uncertainties are unknown but bounded within an ellipsoid. The last part of the paper portrays an exemplary application concerning a nonlinear twin-rotor system.

**Keywords:** fault detection and diagnosis; faults estimation; actuator and sensor fault; observer design; Takagi-Sugeno fuzzy systems

#### **1. Introduction**

Nowadays, owing to the interest in highly efficient systems, industrial companies permanently extend the number of sensors and actuators being used. Indeed, with the advent of IoT the number of components is systematically proliferating. This is mainly due to the fact of their relative low cost and wide accessibility. Irrespective of such an appealing effect, such a growth can increase the chance that actuator and sensor faults appear simultaneously. Moreover, the probability of multiple actuator and sensor faults is also increased. Thus, fault estimation is an important actor in modern Fault Diagnosis (FD) [1–4]. Indeed, it may give a knowledge about the presence, location and size of the fault. Such information is also a necessary ingredient of the active Fault-Tolerant Control (FTC) [5–8], which can accommodate the fault using appropriate control-oriented recovery actions.

There is no doubt that the problem of fault estimation was approached from many different angles. However, most of them focus on estimating either actuator or sensor faults. For a good representative example of recent developments in the area of actuator fault estimation the reader is referred to [9,10]. The case of sensor fault estimation can be realized similarly (cf. Pazera et al. [6] and the references therein).

The number of observers for simultaneous actuator and sensor faults is proliferating as well. However, most of them suffer due to the assumption of a limited fault rate, i.e., the fault derivative (or difference in consecutive discrete samples) is assumed to be close to zero.

Among the existing fault estimation strategies, the crucial attention is focused on the ones based on the sliding-mode observers [11,12], as well as the Kalman filter [13,14] and the related applications.

Taking into account the above discussion, it is also natural that the problem of simultaneous actuator and sensor fault estimation receives growing research attention [8,15–18]. However, the researches are developed for linear systems while the number of strategies capable of handling some classes of nonlinear systems is rather limited. Indeed, a recent state-of-the-art overview clearly indicates the trends for handling multiple simultaneous actuator and sensor fault estimation for nonlinear systems. Indeed, Gu et al. [19] converted an original Lipschitz system into a Linear Parameter Varying (LPV) one by a suitable reformulation of the Lipschitz property. The resulting design approach provides optimal fault estimates according to H<sup>∞</sup> criterion within the frequency domain. Another appealing strategy for Lipschitz class of systems was introduced in Abdollahi [20]. It transforms the system into two subsystems while each of them is affected by either a sensor or an actuator fault, respectively. Finally, two separate Sliding Mode Observers (SMOs) are designed. This scheme suffers from the fact that an asymptotically convergence is guaranteed without considering any robustness to modelling uncertainities and/or external disturbances. A similar SMO was also developed in Hmidi et al. [21] and the authors remove the above drawback by ensuring robustness to bounded disturbances.

An important group of strategies convert the original nonlinear system into an equivalent Takagi-Sugeno one. A representative example of such strategies is presented in Bounemeur et al. [22]. As a result, an adaptive fuzzy estimator is developed capable of handling numerous fault scenarios. Unfortunately, it is devoted to deterministic systems neglecting uncertainty factors. Another strategy is developed in Fu et al. [23] and tackles the systems with switching non-linearities. Unfortunately, it also does not provide appropriate means for settling the robustness issue. Finally, in Shaker [24] a novel multiple integral unknown input observer is developed that is capable of decoupling disturbances. The final group of approaches are the ones for polynomial systems [25,26]. In this case, the design methodology reduces to converting the system by augmenting the state vector with fault variables. This strategy can also cope with process disturbances.

The unappealing feature of the existing approaches to simultaneous actuator and sensor fault estimation is that they do not provide sufficient information about the estimation quality. To handle this issue, two kind of approaches can be utilized: the first one uses post-fault measurements [27] and the second predicts the fault based on its historical values and past measurements [28]. It should be noted that the latter is the only scheme which can be applied in the active FTC. Unfortunately, such schemes struggle with the problem of minimising the effect of fault prediction [9,28].

The scheme proposed in this paper can handle the above mentioned issues, which constitutes the novelty of the proposed approach:


The paper is organised as follows. In Section 2, a simultaneous estimation of actuator and sensor faults problem is described along with suitable feasibility assumptions. Furthermore, a suitable convergence analysis is performed. Section 3 exhibits an illustrative example pertaining the twin-rotor system. Finally, Section 4 concludes the paper.

#### **2. Observer Design**

Let us start from the formulation of a Takagi-Sugeno (T-S) fuzzy system:

$$\begin{split} \mathbf{x}\_{k+1} &= \mathbf{A} \left( \boldsymbol{\nu}\_{k} \right) \mathbf{x}\_{k} + \mathbf{B} \left( \boldsymbol{\nu}\_{k} \right) \boldsymbol{\mu}\_{k} + \mathbf{B}\_{f} \left( \boldsymbol{\nu}\_{k} \right) \vec{f}\_{a,k} + \mathbf{W}\_{1} \boldsymbol{\nu}\_{1,k} \\ &= \sum\_{i=1}^{M} h\_{i} \left( \boldsymbol{\nu}\_{k} \right) \left[ \mathbf{A}^{i} \mathbf{x}\_{k} + \mathbf{B}^{i} \boldsymbol{\mu}\_{k} + \mathbf{B}\_{f}^{i} \vec{f}\_{a,k} \right] + \mathbf{W}\_{1} \boldsymbol{\nu}\_{1,k} \end{split} \tag{1}$$

$$y\_k = \mathcal{C}\mathbf{x}\_k + \mathcal{C}\_f f\_{s,k} + \mathcal{W}\_2 \mathbf{w}\_{2,k} \tag{2}$$

with

$$h\_i(\boldsymbol{\upsilon}\_k) \ge 0, \qquad \forall i = 1, \dots, M\_{\text{ill}} \qquad \sum\_{i=1}^{M\_{\text{UV}}} h\_i(\boldsymbol{\upsilon}\_k) = 1,\tag{3}$$

where *<sup>k</sup>* stands for the discrete time and *<sup>x</sup><sup>k</sup>* <sup>∈</sup> <sup>X</sup> <sup>⊂</sup> <sup>R</sup>*n*, *<sup>u</sup><sup>k</sup>* <sup>∈</sup> <sup>R</sup>*<sup>r</sup>* , *<sup>y</sup><sup>k</sup>* <sup>∈</sup> <sup>R</sup>*<sup>m</sup>* are state, input and output vectors, respectively. Furthermore, **¯** *<sup>f</sup> <sup>a</sup>*,*<sup>k</sup>* <sup>∈</sup> <sup>F</sup>*<sup>a</sup>* <sup>⊂</sup> <sup>R</sup>*na* and *<sup>f</sup>s*,*<sup>k</sup>* <sup>∈</sup> <sup>F</sup>*<sup>s</sup>* <sup>⊂</sup> <sup>R</sup>*ns* signify the actuator and sensor fault vectors, respectively, where *na* and *ns* stand for the number of actuator and sensor faults, respectively. Matrices *A*, *B* and *C* are the known state, input and output matrices, respectively, while *M<sup>m</sup>* stands for the number of submodels. Thus, *C<sup>f</sup>* denotes the sensor fault distribution matrix, where rank(*Cf*) = *ns* and rank(*Bf*(*vk*)) = *na* satisfy the inequality *na* <sup>+</sup> *ns* <sup>≤</sup> *<sup>m</sup>*. This means that it is impossible to estimate more faults than there are measured outputs. Finally, *W*<sup>1</sup> and *W*<sup>2</sup> are the distribution matrices of *w*1,*<sup>k</sup>* and *w*2,*k*, which are exogenous disturbance vectors for the process and measurement uncertainties, respectively. The activation functions *hi*(·) depend on the vector of premise variables *v<sup>k</sup>* = *v*1 *<sup>k</sup>*, *<sup>v</sup>*<sup>2</sup> *<sup>k</sup>*, ..., *<sup>v</sup><sup>p</sup> k T* , which is assumed to depend on measurable variables, e.g., system outputs and known inputs [29]. However, an extension towards unmeasurable premise variables is possible via direct applications of the solution proposed, e.g., in Ichalal et al. [30].

Finally, in the remaining part of the paper, the following notation is used *X* (*vk*) = ∑*<sup>M</sup> <sup>i</sup>*=<sup>1</sup> *hi* (*vk*) *<sup>X</sup><sup>i</sup>* . First, let us start with transforming the state Equation (1) into an equivalent form

$$\mathbf{x}\_{k+1} = \sum\_{i=1}^{M} h\_i \left( \mathbf{v}\_k \right) \left[ \mathbf{A}^i \mathbf{x}\_k + \mathbf{B}^i \boldsymbol{u}\_k \right] + \mathbf{B}\_f \boldsymbol{f}\_{a,k} + \mathbf{W}\_1 \mathbf{w}\_{1,k'} \tag{4}$$

with an auxiliary matrix *<sup>B</sup><sup>f</sup>* satisfying rank(*Bf*) = *na* and an auxiliary actuator fault vector *<sup>f</sup> <sup>a</sup>*,*k*. Comparing Equation (1) and Equation (4), it can be observed that

$$\mathcal{B}\_f(\upsilon\_k) f\_{a,k} = \mathcal{B}\_f f\_{a,k}.\tag{5}$$

Thus, having *<sup>f</sup> <sup>k</sup>*, the original fault vector can be determined with

$$
\vec{f}\_{a,k} = \left(\mathcal{B}\_f(\upsilon\_k)\right)^\dagger \mathcal{B}\_f f\_{a,k}.\tag{6}
$$

where † stands for the pseudo inverse operator. Thus, the objective of further deliberations is to propose a novel observer capable of providing *<sup>x</sup>k*, *<sup>f</sup> <sup>a</sup>*,*<sup>k</sup>* and *<sup>f</sup>s*,*<sup>k</sup>* estimates simultaneously. It should be also noted that the selection of *B<sup>f</sup>* is not critical, i.e., it can be set as *B<sup>f</sup>* = <sup>1</sup> *<sup>M</sup>* <sup>∑</sup>*<sup>M</sup> <sup>i</sup>*=<sup>1</sup> *<sup>B</sup><sup>i</sup> <sup>f</sup>* , which clearly preserves rank(*Bf*) = *na*.

The proposed strategy starts with transforming Equation (4) and Equation (2) into an equivalent descriptor form with the following state variable

$$\mathbf{x}\_k = \begin{bmatrix} \mathbf{x}\_k^T & f\_{a,k-1}^T & f\_{s,k}^T \end{bmatrix}^T. \tag{7}$$

Using the above state vector Equation (7), the system Equation (4) and Equation (2) can be rewritten as:

$$\begin{aligned} E\overline{\mathbf{x}}\_{k+1} &= \overline{A}\left(\overline{\boldsymbol{v}}\_{k}\right)\overline{\mathbf{x}}\_{k} + \overline{B}\left(\overline{\boldsymbol{v}}\_{k}\right)\boldsymbol{u}\_{k} + \overline{\mathbf{W}}\_{1}\overline{\boldsymbol{w}}\_{k'} \\ \boldsymbol{y}\_{k} &= \overline{\mathbf{C}}\mathbf{x}\_{k} + \overline{\mathbf{W}}\_{2}\overline{\boldsymbol{w}}\_{k'} \end{aligned} \tag{9}$$

with:

$$\begin{aligned} \boldsymbol{E} &= \begin{bmatrix} I\_n & -B\_f & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix}, \quad \boldsymbol{A} \begin{pmatrix} \boldsymbol{v}\_k \end{pmatrix} = \begin{bmatrix} \boldsymbol{A} \begin{pmatrix} \boldsymbol{v}\_k \end{pmatrix} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix}, \quad \boldsymbol{B} \begin{pmatrix} \boldsymbol{v}\_k \end{pmatrix} = \begin{bmatrix} \mathbf{B} \begin{pmatrix} \boldsymbol{v}\_k \end{pmatrix} \\ \mathbf{0} \\ \mathbf{0} \end{pmatrix}, \\\ \boldsymbol{\tilde{W}}\_1 &= \begin{bmatrix} \mathbf{W}\_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix}, \quad \boldsymbol{\tilde{W}}\_2 = \begin{bmatrix} \mathbf{0} & \mathbf{W}\_2 \end{bmatrix}, \quad \boldsymbol{\tilde{w}}\_k = \begin{bmatrix} \boldsymbol{w}\_{1,k}^T & \boldsymbol{w}\_{2,k}^T \end{bmatrix}^T, \quad \boldsymbol{\tilde{C}} = \begin{bmatrix} \mathbf{C} & \mathbf{0} & \mathbf{C}\_f \end{bmatrix}. \end{aligned}$$

This means that both faults as well as the state of the system are incorporated into a super state vector *x***¯***k*. Thus, for the purpose of the state estimation obeying Equation (8) and Equation (9), the following observer is proposed:

$$\boldsymbol{z}\_{k+1} = \mathbf{N}\left(\boldsymbol{\upsilon}\_{k}\right)\boldsymbol{z}\_{k} + \mathbf{M}\left(\boldsymbol{\upsilon}\_{k}\right)\boldsymbol{u}\_{k} + \boldsymbol{L}\left(\boldsymbol{\upsilon}\_{k}\right)\boldsymbol{y}\_{k},\tag{10}$$

$$\pounds\_k = \kern-1.72214pt}{\kern-1.72214pt} \pounds\_k + \kern-1.72214pt} \tag{11}$$

where *<sup>z</sup><sup>k</sup>* <sup>∈</sup> <sup>R</sup>*n*+*na*+*ns* is the internal state of the estimator while *<sup>x</sup>***ˆ¯***<sup>k</sup>* <sup>∈</sup> <sup>R</sup>*n*+*na*+*ns* stands for the estimate of Equation (7).

It should be noted that the proposed approach eliminates the usual assumption concerning a bounded rate of change of actuator and sensor fault, which increase the conservatives level of the approaches proposed in the literature (see, e.g., Pazera et al. [28] and the references therein.)

Let us start with assuming that there exist matrices *T*<sup>1</sup> and *T*<sup>2</sup> such that

$$T\_1 E + T\_2 \tilde{\mathcal{C}} = I,\tag{12}$$

or

$$
\begin{bmatrix} T\_1 & T\_2 \end{bmatrix} \begin{bmatrix} E \\ \mathcal{C} \end{bmatrix} = I,\tag{13}
$$

which yields Equation (10) and Equation (11) design condition

$$\text{rank}\left(\begin{bmatrix} E\\ \vec{C} \end{bmatrix}\right) = n + n\_d + n\_s. \tag{14}$$

Under the above assumption, it is possible to derive the error of state estimation with Equation (11), i.e.,

$$\mathfrak{e}\_{k} = \mathfrak{x}\_{k} - \hat{\mathfrak{x}}\_{k} = \mathfrak{x}\_{k} - \mathfrak{z}\_{k} - T\_{2}\mathfrak{y}\_{k} = \left(I - T\_{2}\tilde{\mathcal{C}}\right)\mathfrak{x}\_{k} - \mathfrak{z}\_{k} - T\_{2}\vec{\mathcal{W}}\_{2}\mathfrak{w}\_{k} \tag{15}$$

which, by using Equation (12), boils down to

$$
\mathbf{e}\_k = T\_1 \mathbf{E} \mathbf{\tilde{x}}\_k - \mathbf{z}\_k - T\_2 \mathbf{\tilde{W}}\_2 \mathbf{\tilde{w}}\_k. \tag{16}
$$

Thus, by substituting Equations (8)–(10), the dynamics of state estimation error obeys

$$\begin{aligned} \mathbf{e}\_{k+1} &= T\_1 \mathbf{E} \ddot{\mathbf{x}}\_{k+1} - \mathbf{z}\_{k+1} - T\_2 \tilde{W}\_2 \mathbf{z} \ddot{\mathbf{z}}\_{k+1} = T\_1 \tilde{A} \left( \mathbf{v}\_k \right) \ddot{\mathbf{x}}\_k + T\_1 \tilde{B} \left( \mathbf{v}\_k \right) \mathbf{u}\_k \\ &+ T\_1 \tilde{W}\_1 \ddot{\mathbf{w}}\_k - \mathcal{N} \left( \mathbf{v}\_k \right) \mathbf{z}\_k - \mathcal{M} \left( \mathbf{v}\_k \right) \mathbf{u}\_k - L \left( \mathbf{v}\_k \right) \mathbf{C} \ddot{\mathbf{x}}\_k - L \left( \mathbf{v}\_k \right) \tilde{W}\_2 \ddot{\mathbf{w}}\_k - T\_2 \tilde{W}\_2 \ddot{\mathbf{w}}\_{k+1} \\ &= \left( T\_1 \tilde{A} \left( \mathbf{v}\_k \right) - L \left( \mathbf{v}\_k \right) \tilde{C} - \left( \mathbf{v}\_k \right) T\_1 \mathbf{E} \right) \ddot{\mathbf{x}}\_k + \left( T\_1 \tilde{B} \left( \mathbf{v}\_k \right) - M \left( \mathbf{v}\_k \right) \right) \mathbf{u}\_k \\ &+ N \left( \mathbf{v}\_k \right) \mathbf{e}\_k + \left( T\_1 \tilde{W}\_1 - L \left( \mathbf{v}\_k \right) \tilde{W}\_2 + N \left( \mathbf{v}\_k \right) T\_2 \tilde{W}\_2 \right) \ddot{\mathbf{w}}\_k - T\_2 \tilde{W}\_2 \ddot{\mathbf{w}}\_{k+1}. \end{aligned} \tag{17}$$

From Equation (17), it is evident that the following relations should be satisfied:

$$\left(T\_1 \bar{A}\left(\upsilon\_k\right) - \mathcal{N}\left(\upsilon\_k\right)T\_1 E - L\left(\upsilon\_k\right)\mathcal{C} = 0,\tag{18}$$

$$T\_1 \mathcal{B} \left( \upsilon\_k \right) - M \left( \upsilon\_k \right) = 0. \tag{19}$$

Indeed, by satisfying Equation (18), the term related to *x***¯***<sup>k</sup>* vanished from Equation (17). Similarly, satisfying Equationv Equation (19) means that Equation (18) is no longer dependent on the system input *uk*. Applying Equation (12) to Equation (18) gives

$$T\_1 \bar{A}\left(\upsilon\_k\right) - N\left(\upsilon\_k\right)\left(I - T\_2 \tilde{C}\right) - L\left(\upsilon\_k\right)\tilde{C} = 0,\tag{20}$$

or

$$N\left(\boldsymbol{\upsilon}\_{k}\right) = T\_{1}\bar{A}\left(\boldsymbol{\upsilon}\_{k}\right) - \left(L\left(\boldsymbol{\upsilon}\_{k}\right) - N\left(\boldsymbol{\upsilon}\_{k}\right)T\_{2}\right)\bar{\mathsf{C}}.\tag{21}$$

Finally, by defining

$$\mathbf{K}\left(\boldsymbol{\upsilon}\_{k}\right) = \mathbf{L}\left(\boldsymbol{\upsilon}\_{k}\right) - \mathbf{N}\left(\boldsymbol{\upsilon}\_{k}\right)T\_{2} \,\tag{22}$$

equality Equation (21) boils down to

$$\mathcal{N}\left(\boldsymbol{\upsilon}\_{k}\right) = T\_{1}\boldsymbol{\bar{A}}\left(\boldsymbol{\upsilon}\_{k}\right) - \mathbb{K}\left(\boldsymbol{\upsilon}\_{k}\right)\boldsymbol{\bar{C}}\_{\prime} \tag{23}$$

which makes it possible to transform Equation (17) into

*ek*+<sup>1</sup> = *N* (*vk*) *e<sup>k</sup>* + *<sup>T</sup>*1*W***¯** <sup>1</sup> <sup>−</sup> *<sup>L</sup>* (*vk*) *<sup>W</sup>***¯** <sup>2</sup> <sup>+</sup> *<sup>N</sup>* (*vk*) *<sup>T</sup>*2*W***¯** <sup>2</sup> *<sup>w</sup>***¯** *<sup>k</sup>* <sup>−</sup> *<sup>T</sup>*2*W***¯** <sup>2</sup>*w***¯** *<sup>k</sup>*+<sup>1</sup> = *<sup>T</sup>*1*A***¯** (*vk*) <sup>−</sup> *<sup>K</sup>* (*vk*) *<sup>C</sup>***¯** *<sup>e</sup><sup>k</sup>* <sup>+</sup> *<sup>T</sup>*1*W***¯** <sup>1</sup>*w***¯** *<sup>k</sup>* <sup>−</sup> *<sup>K</sup>* (*vk*) *<sup>W</sup>***¯** <sup>2</sup>*w***¯** *<sup>k</sup>* <sup>−</sup> *<sup>N</sup>* (*vk*) *<sup>T</sup>*2*W***¯** <sup>2</sup>*w***¯** *<sup>k</sup>* <sup>+</sup> *<sup>N</sup>* (*vk*) *<sup>T</sup>*2*W***¯** <sup>2</sup>*w***¯** *<sup>k</sup>* <sup>−</sup> *<sup>T</sup>*2*W***¯** <sup>2</sup>*w***¯** *<sup>k</sup>*+<sup>1</sup> <sup>=</sup> *<sup>T</sup>*1*A***¯** (*vk*) <sup>−</sup> *<sup>K</sup>* (*vk*) *<sup>C</sup>***¯** *ek* <sup>+</sup> *<sup>T</sup>*1*W***¯** <sup>1</sup>*w***¯** *<sup>k</sup>* <sup>−</sup> *<sup>K</sup>* (*vk*) *<sup>W</sup>***¯** <sup>2</sup>*w***¯** *<sup>k</sup>* <sup>−</sup> *<sup>T</sup>*2*W***¯** <sup>2</sup>*w***¯** *<sup>k</sup>*<sup>+</sup>1. (24)

Subsequently, by defining the following super-vector

$$
\mathfrak{w}\_k = \begin{bmatrix} \mathfrak{w}\_k \\ \mathfrak{w}\_{k+1} \end{bmatrix}' \tag{25}
$$

equality Equation (24) can be rewritten into a simpler form

$$\mathfrak{e}\_{k+1} = \mathbf{X}\left(\mathfrak{v}\_k\right)\mathfrak{e}\_k + \mathbf{Z}\left(\mathfrak{v}\_k\right)\mathfrak{w}\_{k\prime} \tag{26}$$

with:

$$\mathbf{X}\left(\boldsymbol{\upsilon}\_{k}\right) = \tilde{\mathbf{A}}\left(\boldsymbol{\upsilon}\_{k}\right) - \mathbf{K}\left(\boldsymbol{\upsilon}\_{k}\right)\tilde{\mathbf{C}}\_{\prime} \qquad \mathbf{Z}\left(\boldsymbol{\upsilon}\_{k}\right) = \tilde{\mathbf{W}}\_{1} - \mathbf{K}\left(\boldsymbol{\upsilon}\_{k}\right)\tilde{\mathbf{W}}\_{2}\mathbf{V}\_{1}$$

where:

$$
\tilde{\mathbf{W}}\_1 \begin{pmatrix} \mathbf{v}\_k \end{pmatrix} = \begin{bmatrix} T\_1 \tilde{\mathbf{W}}\_1 & -T\_2 \tilde{\mathbf{W}}\_2 \end{bmatrix}, \quad \tilde{\mathbf{W}}\_2 = \begin{bmatrix} \tilde{\mathbf{W}}\_2 & \mathbf{0} \end{bmatrix}, \quad \tilde{\mathbf{A}} \begin{pmatrix} \mathbf{v}\_k \end{pmatrix} = T\_1 \tilde{\mathbf{A}} \begin{pmatrix} \mathbf{v}\_k \end{pmatrix}.
$$

For the purpose of further convergence analysis, let us start with reminding the Finsler's Lemma [31]:

**Lemma 1.** *The following expressions are equivalent:*


Let us also define the Lyapunov function

$$V\_k = \mathbf{e}\_k^T \mathbf{P} \mathbf{e}\_{k\prime} \tag{27}$$

with *<sup>P</sup>*  0. Furthermore, the estimation error Equation (26) can be rewritten in an alternative form

$$X\left(\boldsymbol{\upsilon}\_{k}\right)\mathfrak{e}\_{k} + Z\left(\boldsymbol{\upsilon}\_{k}\right)\mathfrak{w}\_{k} - \mathfrak{e}\_{k+1} = \mathbf{0},\tag{28}$$

which implies that the following extended-vector can be defined

$$
\tilde{\mathbf{x}}\_k = \begin{bmatrix} \mathbf{e}\_k^T & \mathbf{z}\mathbf{\bar{v}}\_k^T & \mathbf{e}\_{k+1}^T \end{bmatrix}^T,\tag{29}
$$

and as a consequence the following statement can be formulated

$$\mathbf{R}\left(\boldsymbol{\upsilon}\_{k}\right)\begin{bmatrix}\boldsymbol{\mathsf{e}}\_{k} \\ \boldsymbol{\mathsf{w}}\_{k} \\ \boldsymbol{\mathsf{e}}\_{k+1}\end{bmatrix} = \mathbf{0}.\tag{30}$$

Thus, based on Equation (28), Equation (29) and Equation (30), it can be shown that

$$\mathcal{R}\left(\boldsymbol{\upsilon}\_{k}\right) = \begin{bmatrix} \mathbf{X}\left(\boldsymbol{\upsilon}\_{k}\right) \ \mathbf{Z}\left(\boldsymbol{\upsilon}\_{k}\right) \ -\mathbf{I} \end{bmatrix},\tag{31}$$

where

$$\mathcal{R}\left(\boldsymbol{v}\_{k}\right)\tilde{\boldsymbol{x}}\_{k} = \mathbf{0},\tag{32}$$

in order to satisfy Equation (28). The convergence of the proposed observer is to be determined with the so-called Quadratic Boundedness (QB) approach [32]. This technique can be perceived as an extension of the usual Lyapunov approach towards the systems with external bounded disturbances. The usefulness of QB approach was proven in many papers while in Pazera et al. [28] it was proven that the standard H<sup>∞</sup> framework can be perceived as a special case of QB. To use the QB approach, it is necessary to assume that *w***˜** *<sup>k</sup>* is bounded by the following ellipsoid:

$$\mathbb{E}\_w = \{ \mathfrak{w}\_k \, : \, \mathfrak{w}\_k^T \mathbf{Q}\_w \mathfrak{w}\_k \le 1 \}, \tag{33}$$

with *<sup>Q</sup><sup>w</sup>*  0. Under such assumptions, the following definition can be recalled:

**Definition 1.** *The system Equation* (26) *is strictly quadratically bounded for all <sup>w</sup>***˜** *<sup>k</sup>* <sup>∈</sup> <sup>E</sup>*w, <sup>k</sup>* <sup>≥</sup> <sup>0</sup>*, if Vk* <sup>&</sup>gt; <sup>1</sup> <sup>=</sup><sup>⇒</sup> *Vk*<sup>+</sup><sup>1</sup> <sup>−</sup> *Vk* <sup>&</sup>lt; <sup>0</sup> *for any <sup>w</sup>***˜** *<sup>k</sup>* <sup>∈</sup> <sup>E</sup>*w.*

As shown in Alessandri et al. [32] and Pazera et al. [28], the stability condition associated with

$$\left(V\_{k+1} - \left(1 - \mathfrak{a}\right)V\_k - \mathfrak{a}\mathfrak{T}\mathfrak{v}\_k^T\mathbf{Q}\_w\mathfrak{v}\mathfrak{v}\_k < 0\right) \tag{34}$$

with 0 < *α* < 1. Based on the above considerations, the following theorem is established:

**Theorem 1.** *The observer-based system Equation* (26) *is strictly quadratically bounded for all <sup>w</sup>***˜** *<sup>k</sup>* <sup>∈</sup> <sup>E</sup>*<sup>w</sup> if there exist matrices <sup>P</sup>* <sup>0</sup>*, <sup>U</sup>, <sup>N</sup>***¯** *as well as <sup>α</sup>* <sup>∈</sup> (0, 1) *such that the following holds:*

$$
\begin{bmatrix}
a P - P & \mathbf{0} & \tilde{A}^T \begin{pmatrix} \boldsymbol{\upsilon}\_k \end{pmatrix} \boldsymbol{\mathcal{U}}^T - \tilde{\boldsymbol{\mathcal{C}}}^T \tilde{\mathbf{N}}^T \left(\boldsymbol{\upsilon}\_k\right) \\
\mathbf{0} & -a \boldsymbol{Q}\_w & \tilde{\mathbf{W}}\_1^T \boldsymbol{\mathcal{U}}^T - \tilde{\mathbf{W}}\_2^T \tilde{\mathbf{N}}^T \left(\boldsymbol{\upsilon}\_k\right) \\
\boldsymbol{\mathsf{U}}\tilde{\boldsymbol{A}}\left(\boldsymbol{\upsilon}\_k\right) - \tilde{\mathbf{N}}\left(\boldsymbol{\upsilon}\_k\right) \tilde{\mathbf{C}} & \boldsymbol{\mathsf{U}}\tilde{\mathbf{W}}\_1 - \tilde{\mathbf{N}}\left(\boldsymbol{\upsilon}\_k\right) \tilde{\mathbf{W}}\_2 & \boldsymbol{P} - \boldsymbol{\mathsf{U}} - \boldsymbol{\mathsf{U}}^T
\end{bmatrix} \prec 0,\tag{35}
$$

**Proof.** Using Equation (34) and setting

$$Q = \begin{bmatrix} aP - P & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & -aQ\_w & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & P \end{bmatrix},\tag{36}$$

$$
\tilde{\mathbf{M}} = \begin{bmatrix} \mathbf{0}^T & \mathbf{0}^T & \mathbf{U}^T \end{bmatrix}^T \,. \tag{37}
$$

along with Lemma 1 leads to

$$
\begin{bmatrix}
\alpha \mathbf{P} - \mathbf{P} & \mathbf{0} & \mathbf{X}^T \begin{pmatrix} \boldsymbol{\upsilon}\_k \end{pmatrix} \mathbf{U}^T \\
\mathbf{0} & -\alpha \mathbf{Q}\_w & \mathbf{Z}^T \begin{pmatrix} \boldsymbol{\upsilon}\_k \end{pmatrix} \mathbf{U}^T \\
\mathbf{U} \mathbf{X} \begin{pmatrix} \boldsymbol{\upsilon}\_k \end{pmatrix} & \mathbf{U} \mathbf{Z} \begin{pmatrix} \boldsymbol{\upsilon}\_k \end{pmatrix} \mathbf{U}^T \\
\mathbf{U} \mathbf{X} \begin{pmatrix} \boldsymbol{\upsilon}\_k \end{pmatrix} & \mathbf{U} \mathbf{Z} \begin{pmatrix} \boldsymbol{\upsilon}\_k \end{pmatrix} & \mathbf{P} - \mathbf{U} - \mathbf{U}^T \end{pmatrix}
\tag{38}
$$

Setting

$$\mathcal{U}X\left(\boldsymbol{\upsilon}\_{k}\right) = \mathcal{U}\left(\tilde{A}\left(\boldsymbol{\upsilon}\_{k}\right) - \mathcal{K}\left(\boldsymbol{\upsilon}\_{k}\right)\tilde{\mathcal{C}}\right) = \mathcal{U}\tilde{A}\left(\boldsymbol{\upsilon}\_{k}\right) - \tilde{\mathcal{N}}\left(\boldsymbol{\upsilon}\_{k}\right)\tilde{\mathcal{C}},\tag{39}$$

$$\mathbf{U}\mathbf{Z}\left(\boldsymbol{\upsilon}\_{k}\right) = \mathbf{U}\left(\tilde{\boldsymbol{W}}\_{1} - \boldsymbol{\mathsf{K}}\left(\boldsymbol{\upsilon}\_{k}\right)\tilde{\boldsymbol{W}}\_{2}\right) = \mathbf{U}\tilde{\boldsymbol{W}}\_{1} - \tilde{\boldsymbol{\mathsf{N}}}\left(\boldsymbol{\upsilon}\_{k}\right)\tilde{\boldsymbol{W}}\_{2} \tag{40}$$

into Equation (38) completes the proof.

Irrespective of the incontestable appeal of the proposed solution Equation (35), its numerical tractability is feasible only if it is transformed to the set of linear matrix inequalities under fixed *α*:

$$
\begin{bmatrix}
\boldsymbol{a}\boldsymbol{P} - \boldsymbol{P} & \mathbf{0} & \boldsymbol{A}^{\boldsymbol{i}^T}\boldsymbol{\mathcal{U}}^T - \boldsymbol{\mathcal{C}}^T\boldsymbol{\mathcal{N}}^T \\
\mathbf{0} & -\boldsymbol{a}\boldsymbol{Q}\_w & \boldsymbol{\mathcal{W}}\_1^T\boldsymbol{\mathcal{U}}^T - \boldsymbol{\mathcal{W}}\_2^T\boldsymbol{\mathcal{N}}^T \\
\boldsymbol{\mathsf{U}}\boldsymbol{\tilde{A}}^{\boldsymbol{i}} - \boldsymbol{\mathcal{N}}^{\boldsymbol{i}}\boldsymbol{\tilde{C}} & \boldsymbol{\mathsf{U}}\boldsymbol{\tilde{W}}\_1 - \boldsymbol{\bar{N}}^{\boldsymbol{i}}\boldsymbol{\tilde{W}}\_2 & \boldsymbol{P} - \boldsymbol{\mathcal{U}} - \boldsymbol{\mathcal{U}}^T \\
\end{bmatrix} \prec \boldsymbol{0}, \quad \boldsymbol{i} = 1, \ldots, M. \tag{41}
$$

Finally, the design procedure boils down to:

**Step 0:** Determine *T*<sup>1</sup> and *T*<sup>2</sup> by solving Equation (12). **Step 1:** Set *α*, 0 < *α* < 1 and determine *N***¯** *<sup>i</sup>* and *U* by solving Equation (41). **Step 2:** Calculate:

$$
\mathbb{K}^i = \mathcal{U}^{-1} \mathbb{N}^i,\tag{42}
$$

$$\mathbf{N}^{i} = T\_1 \tilde{\mathbf{A}}^{l} - \mathbf{K}^{i} \vec{\mathbf{C}}\_{l} \tag{43}$$

$$\mathbf{M}^i = T\_1 \mathbf{B}^i,\tag{44}$$

$$L^i = \mathbb{K}^i + \mathbf{N}^i T\_2. \tag{45}$$

While the application procedure leads to:

**Step 0:** Substitute *k* = 0 and set the initial conditions *z*0. **Step 1:** Calculate *zk*+<sup>1</sup> and *x***ˆ***<sup>k</sup>* with Equations (10) and (11). **Step 2:** Set *k* = *k* + 1 and go to *Step 1*.

As demonstrated in Pazera et al. [28], the value of *α* has a direct influence of the convergence rate of the underlying observer. On the other hand, matrix *P* shapes the ellipsoidal bound of the estimation error:

$$\mathbf{e}\_k^T \mathbf{P} \mathbf{e}\_k \le \mathbb{Z}\_k \quad \mathbb{Z}\_k = 1 + (1 - \boldsymbol{\alpha})^k (1 - \mathbf{e}\_0^T \mathbf{P} \mathbf{e}\_0),\tag{46}$$

which can be directly used to determine the uncertainty intervals of the state variables [28]:

$$
\mathfrak{Ax}\_{k,i} - \sqrt{\zeta\_k \mathbf{c}\_i^T \mathbf{P}^{-1} \mathbf{c}\_i} \le \mathfrak{x}\_{k,i} \le \mathfrak{x}\_{k,i} + \sqrt{\zeta\_k \mathbf{c}\_i^T \mathbf{P}^{-1} \mathbf{c}\_i} \tag{47}
$$

where *c<sup>i</sup>* signifies *i*th column of an identity matrix of an appropriate dimension.

#### **3. Case Study: Twin-Rotor System**

The aim of this section is to validate the performance and correctness of the novel observer. Accordingly, the proposed actuator and sensor fault estimation method has been implemented to the Twin-Rotor System (TRS) [28], illustrated in Figure 1. Note that all equations and a detailed nonlinear model of the TRS can be found in [28]. This model was further transformed into the Takagi-Sugeno form using the dedicated methodology proposed in Rotondo et al. [33]. The matrices shaping the model Equation (1) and Equation (2) are presented in the Appendix A.

**Figure 1.** A nonlinear twin-rotor MIMO system.

The state vector of the system is given by

$$\mathbf{x} = \begin{bmatrix} \omega\_{\upsilon} & \Omega\_{\upsilon} & \theta\_{\upsilon} & \omega\_{h} & \Omega\_{h} & \theta\_{h} \end{bmatrix}^{T} \tag{48}$$

hence the input vector is described by

$$\boldsymbol{\mu} = \begin{bmatrix} \boldsymbol{\mu}\_{\boldsymbol{\mathcal{U}}} & \boldsymbol{\mu}\_{\boldsymbol{\mathcal{U}}} \end{bmatrix}^{T} \text{ \,} \tag{49}$$

where *ω<sup>v</sup>* and Ω*<sup>v</sup>* stand for the rotational and angular velocities of the main rotor while *θ<sup>v</sup>* signifies the pitch angle of the beam. Furthermore, *ω<sup>h</sup>* and Ω*<sup>h</sup>* are the tail rotor rotational as well as angular velocities whilst *θ<sup>h</sup>* is the yaw angle of the beam. Finally, *uv* and *uh* stand for actuation due to the main and tail DC motors, respectively. It should be noted that both inputs *uv* and *uh* operate in range (−1, 1). Moreover, all details and parameters can be found in manufacturer's user manual [34] and in Pazera et al. [28].

The following fault scenarios have been considered:

$$f\_{a,1,k} = \begin{cases} -0.1 \cdot u\_{1,k} & 3000 \le k \le 6500 \\ 0 & \text{otherwise,} \end{cases} \tag{50}$$

$$f\_{a,2,k} = \begin{cases} 0.1 \cdot \mathfrak{u}\_{2,k} & 4500 \le k \le 8000 \\ 0 & \text{otherwise,} \end{cases} \tag{51}$$

$$f\_{s,1,k} = \begin{cases} y\_k - 2.5 & 5000 \le k \le 7000\\ 0 & \text{otherwise,} \end{cases} \tag{52}$$

$$f\_{s,2,k} = \begin{cases} y\_k + 1.6 & \text{6000} \le k \le 9500\\ 0 & \text{otherwise} \end{cases} \tag{53}$$

According to the approach proposed along with Equation (6):

$$\mathbf{B}\_{f} = \frac{1}{M\_{m}} \sum\_{i=1}^{M\_{m}} \mathbf{B}\_{f}^{i} = \begin{bmatrix} 49.4084 \cdot 10^{-6} & -66.2480 \cdot 10^{-6} \\ 14.3375 \cdot 10^{-3} & -13.2367 \cdot 10^{-3} \\ 303.3781 & 0 \\ 6.9349 \cdot 10^{-6} & 3.6497 \cdot 10^{-6} \\ 1.3860 \cdot 10^{-3} & 1.0854 \cdot 10^{-3} \\ 0 & 58.0016 \end{bmatrix} \tag{54}$$

which means that all actuators are considered as possibly faulty and their faults have to be estimated. Whilst the sensor fault distribution matrix is obtained from matrix C (see Appendix A, Equation (A25)) by extracting the rows corresponding to the first and second sensor, respectively and it is given as follows

$$\mathcal{C}\_f = \begin{bmatrix} 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{bmatrix}^T. \tag{55}$$

It can be clearly viewed that within this scenario, both actuator and sensor faults are partially at the same time. Moreover, the considered faults are constant biases, which are either positive or negative. From Equations (50)–(53), it can be clearly viewed that actuator and sensor faults appear partially at the same time. The considered actuator fault Equation (50) pertains a 10% performance decrease of the first actuator, while Equation (51) concerns a 10% performance increase of the second actuator. Moreover, the considered sensor faults Equation (52) and Equation (53) are represented by either positive or negative constant biases equal to 1.6 or −2.5, respectively. Thus, they denote a significant sensor readings inaccuracies. However, the distribution matrix of sensor fault denotes that they influence the measurements of the angular velocity of the tail rotor Ω*<sup>h</sup>* as well as the angular velocity of the main rotor Ω*v*. Firstly, let us present a comparison between nonlinear estimation approach described in Pazera et al. [28] and the proposed Takagi-Sugeno scheme, which is shown in Figure 2a,b. As it can be observed, the Takagi-Sugeno model-based approach follows the response of the original nonlinear model-based one, and hence, it does not impair the quality of state and fault estimates significantly.

Figure 3a,b illustrate the real actuator faults of the main *<sup>f</sup> <sup>a</sup>*,1 and tail *<sup>f</sup> <sup>a</sup>*,2 rotor with blue dash-dotted lines along with their estimates given with red dashed lines. Additionally, the real and estimated values are overbounded by uncertainty intervals which are indicated with black dashed lines. Moreover, the real sensor faults *<sup>f</sup>s*,1 and *<sup>f</sup>s*,2 are presented in Figure 4a,b with blue dash-dotted lines as well as their estimates depicted by red dashed lines along with the uncertainty intervals indicated with black dashed lines. From these figures it is easily seen that the actuator and sensor faults were estimated with a very good accuracy under the measurement uncertainties. Consequently, the system states have been correctly reconstructed.

**Figure 2.** Comparison between nonlinear and Takagi-Sugeno response of the system: angular velocity of the main rotor Ω*<sup>h</sup>* (**a**) and angular velocity of the tail rotor Ω*<sup>v</sup>* (**b**).

**Figure 3.** Actuator faults *<sup>f</sup> <sup>a</sup>*,1 (**a**) and *<sup>f</sup> <sup>a</sup>*,2 (**b**).

**Figure 4.** Sensor faults *<sup>f</sup>s*,1 (**a**) and *<sup>f</sup>s*,2 (**b**).

Figure 5a,b show the real rotational velocities of the main *ω<sup>v</sup>* and tail *ω<sup>h</sup>* rotor with blue dash-dotted lines. Their measured outputs are presented by light green solid lines along with their estimates given with red dashed lines as well as the uncertainty intervals indicated with black dashed lines. Moreover, Figure 6a,b present the angular velocities of the main Ω*<sup>v</sup>* and of the tail Ω*<sup>h</sup>* rotor with blue dash-dotted lines whilst their estimates are given with red dashed lines as well as the measured outputs depicted by light green solid lines. Additionally, the uncertainty intervals are indicated with black dashed lines. As can be observed in Figure 6a,b, irrespective of the intermittent fault (marked in light green) the state estimates are very close to the original states. The rotational velocities of the main and tail rotor and the pitch angle of the beam have been accurately estimated even if the actuator and sensor faults occurred simultaneously. This fact is illustrated in Figure 7a,b and Figure 8, which show the evolution of the state estimation error for the above variables. The figures clearly show that the states are properly reconstructed under the actuator glitch along with positive and negative incorrect readings of the sensor.

**Figure 5.** State variables *ω<sup>v</sup>* (**a**) and *ω<sup>h</sup>* (**b**).

**Figure 6.** State variables Ω*<sup>v</sup>* (**a**) and Ω*<sup>h</sup>* (**b**).

**Figure 7.** Evolution of the state estimation error for the rotational velocities of the main rotor (**a**) and of the tail rotor (**b**).

**Figure 8.** Evolution of the state estimation error for the pitch angle of the beam.

#### **4. Concluding Remarks**

An important aspect of the paper was to develop an observer for the Takagi-Sugeno systems which would be capable of estimating both state with both actuator and sensor. The first part of the proposed strategy was developed in such a way that no adaptive observer is required to obtain the sensor fault estimates. In other words, it is obtained directly by simply transforming an output equation. The second part of the observer concerns the adaptive actuator fault estimation and the state as well. In order to achieve robustness against process and measurement uncertainties and to ensure the stability performance a Quadratic Boundedness strategy is employed. For the purpose of obtaining the observer gain matrices, a set of LMIs is needed to be solved. The other part of the paper presents a performance validation. It has been achieved by implementing the approach to a nonlinear Twin-Rotor MIMO System. The obtained results plainly confirm that the desired properties of the observer have been accomplished. The state has been estimated properly even in the case of both the actuator and sensor fault case. While analysing an actuator fault, it can be concluded that its estimation error is at a very low level in both cases: fault and fault-free as well. Furthermore, the obtained sensor fault estimates quality strongly depends on the quality of the state estimation. The future work will be devoted to applying such an approach to a fault-tolerant control scheme for a T–S fuzzy systems.

**Author Contributions:** M.P. derived the design procedure for fault estimation. M.W. contributed to the development of the observer structure and determination of the uncertainty interval calculation rules. N.K. contributed to the experiments and described the experimental results. M.B. contributed to the development of the convergence conditions and review of the literature. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported by the National Science Centre, Poland under Grant: UMO-2017/27/B/ST7/00620. **Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*A*<sup>1</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9709 · <sup>10</sup>−<sup>3</sup> 4.5486 · <sup>10</sup>−<sup>7</sup> −5.1930 · <sup>10</sup>−<sup>7</sup> −1.7300 · <sup>10</sup>−<sup>9</sup> −1.5629 · <sup>10</sup>−<sup>14</sup> 0 9.9418 · <sup>10</sup>−<sup>1</sup> 8.5060 · <sup>10</sup>−<sup>5</sup> −1.0376 · <sup>10</sup>−<sup>4</sup> −5.1856 · <sup>10</sup>−<sup>7</sup> −6.2171 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9988 · <sup>10</sup>−<sup>1</sup> 9.9785 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −2.4758 · <sup>10</sup>−<sup>2</sup> 9.9566 · <sup>10</sup>−<sup>1</sup> 3.4967 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A1) *A*<sup>2</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9709 · <sup>10</sup>−<sup>3</sup> 4.5486 · <sup>10</sup>−<sup>7</sup> −4.6737 · <sup>10</sup>−<sup>7</sup> −1.5570 · <sup>10</sup>−<sup>9</sup> −1.4066 · <sup>10</sup>−<sup>14</sup> 0 9.9418 · <sup>10</sup>−<sup>1</sup> 8.5060 · <sup>10</sup>−<sup>5</sup> −9.3380 · <sup>10</sup>−<sup>5</sup> −4.6671 · <sup>10</sup>−<sup>7</sup> −5.5954 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9988 · <sup>10</sup>−<sup>1</sup> 9.9785 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −2.4758 · <sup>10</sup>−<sup>2</sup> 9.9566 · <sup>10</sup>−<sup>1</sup> 3.4967 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A2) *A*<sup>3</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9707 · <sup>10</sup>−<sup>3</sup> 4.5690 · <sup>10</sup>−<sup>7</sup> −5.2162 · <sup>10</sup>−<sup>7</sup> −1.7378 · <sup>10</sup>−<sup>9</sup> −1.5699 · <sup>10</sup>−<sup>14</sup> 0 9.9415 · <sup>10</sup>−<sup>1</sup> 8.5440 · <sup>10</sup>−<sup>5</sup> −1.0422 · <sup>10</sup>−<sup>4</sup> −5.2088 · <sup>10</sup>−<sup>7</sup> −6.2449 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9988 · <sup>10</sup>−<sup>1</sup> 9.9785 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −2.4758 · <sup>10</sup>−<sup>2</sup> 9.9566 · <sup>10</sup>−<sup>1</sup> 3.4967 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A3) *A*<sup>4</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9707 · <sup>10</sup>−<sup>3</sup> 4.5690 · <sup>10</sup>−<sup>7</sup> −4.6946 · <sup>10</sup>−<sup>7</sup> −1.5640 · <sup>10</sup>−<sup>9</sup> −1.4129 · <sup>10</sup>−<sup>14</sup> 0 9.9415 · <sup>10</sup>−<sup>1</sup> 8.5440 · <sup>10</sup>−<sup>5</sup> −9.3798 · <sup>10</sup>−<sup>5</sup> −4.6879 · <sup>10</sup>−<sup>7</sup> −5.6204 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9988 · <sup>10</sup>−<sup>1</sup> 9.9785 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −2.4758 · <sup>10</sup>−<sup>2</sup> 9.9566 · <sup>10</sup>−<sup>1</sup> 3.4967 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A4) *A*<sup>5</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9709 · <sup>10</sup>−<sup>3</sup> 4.5486 · <sup>10</sup>−<sup>7</sup> −5.1929 · <sup>10</sup>−<sup>7</sup> −1.7300 · <sup>10</sup>−<sup>9</sup> −1.5629 · <sup>10</sup>−<sup>14</sup> 0 9.9418 · <sup>10</sup>−<sup>1</sup> 8.5060 · <sup>10</sup>−<sup>5</sup> −1.0375 · <sup>10</sup>−<sup>4</sup> −5.1856 · <sup>10</sup>−<sup>7</sup> −6.2171 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9981 · <sup>10</sup>−<sup>1</sup> 9.9783 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −3.7387 · <sup>10</sup>−<sup>2</sup> 9.9559 · <sup>10</sup>−<sup>1</sup> 3.4966 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A5) *A*<sup>6</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9709 · <sup>10</sup>−<sup>3</sup> 4.5486 · <sup>10</sup>−<sup>7</sup> −4.6736 · <sup>10</sup>−<sup>7</sup> −1.5570 · <sup>10</sup>−<sup>9</sup> −1.4066 · <sup>10</sup>−<sup>14</sup> 0 9.9418 · <sup>10</sup>−<sup>1</sup> 8.5060 · <sup>10</sup>−<sup>5</sup> −9.3378 · <sup>10</sup>−<sup>5</sup> −4.6670 · <sup>10</sup>−<sup>7</sup> −5.5954 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9981 · <sup>10</sup>−<sup>1</sup> 9.9783 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −3.7387 · <sup>10</sup>−<sup>2</sup> 9.9559 · <sup>10</sup>−<sup>1</sup> 3.4966 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A6)


*A*<sup>14</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9709 · <sup>10</sup>−<sup>3</sup> 3.1840 · <sup>10</sup>−<sup>7</sup> −4.6737 · <sup>10</sup>−<sup>7</sup> −1.5570 · <sup>10</sup>−<sup>9</sup> −1.4066 · <sup>10</sup>−<sup>14</sup> 0 9.9418 · <sup>10</sup>−<sup>1</sup> 5.9542 · <sup>10</sup>−<sup>5</sup> −9.3381 · <sup>10</sup>−<sup>5</sup> −4.6671 · <sup>10</sup>−<sup>7</sup> −5.5954 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9989 · <sup>10</sup>−<sup>1</sup> 9.9785 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −2.1851 · <sup>10</sup>−<sup>2</sup> 9.9567 · <sup>10</sup>−<sup>1</sup> 3.4967 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A14) *A*<sup>15</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9707 · <sup>10</sup>−<sup>3</sup> 3.1983 · <sup>10</sup>−<sup>7</sup> −5.2162 · <sup>10</sup>−<sup>7</sup> −1.7378 · <sup>10</sup>−<sup>9</sup> −1.5699 · <sup>10</sup>−<sup>14</sup> 0 9.9415 · <sup>10</sup>−<sup>1</sup> 5.9808 · <sup>10</sup>−<sup>5</sup> −1.0422 · <sup>10</sup>−<sup>4</sup> −5.2088 · <sup>10</sup>−<sup>7</sup> −6.2449 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9989 · <sup>10</sup>−<sup>1</sup> 9.9785 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −2.1851 · <sup>10</sup>−<sup>2</sup> 9.9567 · <sup>10</sup>−<sup>1</sup> 3.4967 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A15) *A*<sup>16</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 9.9707 · <sup>10</sup>−<sup>3</sup> 3.1983 · <sup>10</sup>−<sup>7</sup> −4.6946 · <sup>10</sup>−<sup>7</sup> −1.5640 · <sup>10</sup>−<sup>9</sup> −1.4129 · <sup>10</sup>−<sup>14</sup> 0 9.9415 · <sup>10</sup>−<sup>1</sup> 5.9808 · <sup>10</sup>−<sup>5</sup> −9.3798 · <sup>10</sup>−<sup>5</sup> −4.6880 · <sup>10</sup>−<sup>7</sup> −5.6205 · <sup>10</sup>−<sup>12</sup> 0 0 6.6291 · <sup>10</sup>−<sup>1</sup> 00 0 0 0 0 9.9989 · <sup>10</sup>−<sup>1</sup> 9.9785 · <sup>10</sup>−<sup>3</sup> 1.7793 · <sup>10</sup>−<sup>7</sup> 00 0 −2.1851 · <sup>10</sup>−<sup>2</sup> 9.9567 · <sup>10</sup>−<sup>1</sup> 3.4967 · <sup>10</sup>−<sup>5</sup> 0 0 0 0 0 9.0333 · <sup>10</sup>−<sup>1</sup> ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A16) *B*<sup>1</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.7998 · <sup>10</sup>−<sup>5</sup> −5.4643 · <sup>10</sup>−<sup>5</sup> 1.6830 · <sup>10</sup>−<sup>2</sup> −1.0918 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>2</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.7998 · <sup>10</sup>−<sup>5</sup> −5.4643 · <sup>10</sup>−<sup>5</sup> 1.6830 · <sup>10</sup>−<sup>2</sup> −1.0918 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A17) *B*<sup>3</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.8257 · <sup>10</sup>−<sup>5</sup> −5.4887 · <sup>10</sup>−<sup>5</sup> 1.6905 · <sup>10</sup>−<sup>2</sup> −1.0967 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>4</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.8257 · <sup>10</sup>−<sup>5</sup> −5.4887 · <sup>10</sup>−<sup>5</sup> 1.6905 · <sup>10</sup>−<sup>2</sup> −1.0967 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A18) *B*<sup>5</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.7998 · <sup>10</sup>−<sup>5</sup> −1.2162 · <sup>10</sup>−<sup>4</sup> 1.6830 · <sup>10</sup>−<sup>2</sup> −2.4301 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9348 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3859 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>6</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.7998 · <sup>10</sup>−<sup>5</sup> −1.2162 · <sup>10</sup>−<sup>4</sup> 1.6830 · <sup>10</sup>−<sup>2</sup> −2.4301 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9348 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3859 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A19) *B*<sup>7</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.8257 · <sup>10</sup>−<sup>5</sup> −1.2217 · <sup>10</sup>−<sup>4</sup> 1.6905 · <sup>10</sup>−<sup>2</sup> −2.4410 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9348 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3859 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>8</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 5.8257 · <sup>10</sup>−<sup>5</sup> −1.2217 · <sup>10</sup>−<sup>4</sup> 1.6905 · <sup>10</sup>−<sup>2</sup> −2.4410 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9348 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3859 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A20) *B*<sup>9</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0598 · <sup>10</sup>−<sup>5</sup> −1.0576 · <sup>10</sup>−<sup>5</sup> 1.1781 · <sup>10</sup>−<sup>2</sup> −2.1132 · <sup>10</sup>−<sup>3</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9350 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>10</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0598 · <sup>10</sup>−<sup>5</sup> −1.0576 · <sup>10</sup>−<sup>5</sup> 1.1781 · <sup>10</sup>−<sup>2</sup> −2.1132 · <sup>10</sup>−<sup>3</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9350 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A21)

*B*<sup>11</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0780 · <sup>10</sup>−<sup>5</sup> −1.0623 · <sup>10</sup>−<sup>5</sup> 1.1834 · <sup>10</sup>−<sup>2</sup> −2.1226 · <sup>10</sup>−<sup>3</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9350 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>12</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0780 · <sup>10</sup>−<sup>5</sup> −1.0623 · <sup>10</sup>−<sup>5</sup> 1.1834 · <sup>10</sup>−<sup>2</sup> −2.1226 · <sup>10</sup>−<sup>3</sup> 3.0338 · 102 <sup>0</sup> 6.9350 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A22) *B*<sup>13</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0598 · <sup>10</sup>−<sup>5</sup> −7.7558 · <sup>10</sup>−<sup>5</sup> 1.1781 · <sup>10</sup>−<sup>2</sup> −1.5496 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>14</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0598 · <sup>10</sup>−<sup>5</sup> −7.7558 · <sup>10</sup>−<sup>5</sup> 1.1781 · <sup>10</sup>−<sup>2</sup> −1.5496 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A23) *B*<sup>15</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0780 · <sup>10</sup>−<sup>5</sup> −7.7904 · <sup>10</sup>−<sup>5</sup> 1.1834 · <sup>10</sup>−<sup>2</sup> −1.5566 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , *B*<sup>16</sup> = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 4.0780 · <sup>10</sup>−<sup>5</sup> −7.7904 · <sup>10</sup>−<sup>5</sup> 1.1834 · <sup>10</sup>−<sup>2</sup> −1.5566 · <sup>10</sup>−<sup>2</sup> 3.0338 · <sup>10</sup><sup>2</sup> <sup>0</sup> 6.9349 · <sup>10</sup>−<sup>6</sup> 3.6497 · <sup>10</sup>−<sup>6</sup> 1.3860 · <sup>10</sup>−<sup>3</sup> 1.0854 · <sup>10</sup>−<sup>3</sup> 0 58.002 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (A24)

$$\mathbf{W}\_1 = 0.01\mathbf{I}, \qquad \mathbf{W}\_2 = 0.05\mathbf{I}, \qquad \mathbf{C} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}. \tag{A25}$$

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
