Fault Tolerant Control for MIMO Nonlinear Systems via MPC Based on MIMO ARX-Laguerre Multiple Models

Abstract

In this article, we propose a fault tolerant control for multiple-input multiple-output (MIMO) nonlinear systems via model predictive control. The MIMO nonlinear systems are approximated by MIMO ARX-Laguerre multiple models. The latter is obtained by expanding a discrete-time MIMO ARX multiple model parameters on Laguerre orthonormal bases. The resulting model ensures an efficient complexity reduction with respect to the classical MIMO ARX multiple models. This parametric complexity reduction still subjects to an optimal choice of the Laguerre poles defining Laguerre bases. The parameter and structure identifications of the MIMO ARX-Laguerre multiple models are achieved by the recursive method and a metaheuristic algorithm, respectively. The proposed model is built from the system input/output observations and is used to synthesize a MIMO nonlinear fault tolerant control algorithm via MPC. So, we develop a fault detection and isolation (FDI) scheme based on the proposed model. The scheme of the fault detection is applied at every step of MPC control calculation, where we determine the actuator faults and we use it in the MPC optimization problem to determine the new control with respect to the actuator faults. The proposed strategy is tested on numerical simulation and validated on the real system.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The security of technological systems is a major concern in the last decade [1–4]. Ensuring the security and the environment of a process require the knowledge of its operating status as finely as possible at every moment. In particular, we have to be able to decide if the working system is normal or if a malfunction has occurred. In this case, it is interesting to know the nature of this dysfunction, which is the main objective of the diagnosis. So, in the context of increasing autonomy of the systems, once the dysfunction is detected and identified, this knowledge must be taken into account when calculating a closed-loop control law to counter the influence of defects. This strategy is entitled fault tolerant control (FTC) [5].

The goal of this work is to synthesize fault tolerant control for nonlinear multiple-input multiple-output (MIMO) systems via model predictive control (MPC) based on reduced complexity multiple models. Model-based predictive control is a well-established online control strategy which iteratively computes control signals by solving an optimization problem over a future time horizon under certain process constraints [6–9]. This optimization uses a prediction model of the future plant behavior. The closed-loop performance depends on the choice of an appropriate model for prediction and several tuning parameters. To model MIMO nonlinear systems, several models are used like neural network [10], MIMO nonlinear autoregressive with exogenous input model (NARX) [11], fuzzy logic model [12], and MIMO autoregressive with exogenous input (ARX) multiple models known as MIMO Takagi–Sugeno model [13]. However, these models are constrained by a high parameter number. Furthermore, the complexity of MIMO nonlinear models handicaps the synthesis of a control law by increasing the computation time of the control. To overcome this problem, several works are developed in the literature in the case of single-input single-output (SISO) linear and nonlinear models and MIMO linear model by expansion on Laguerre orthonormal bases [14–18].

In this context, we propose in this paper to reduce the parametric complexity of ARX MIMO multiple models by decomposing its parameter associated with the inputs and the outputs on independent Laguerre orthonormal bases. This decomposition can be realized since the coefficients of the ARX MIMO multiple models are absolutely summable on $\left[0,\infty \right)$ in the sense of the bounded-input bounded-output (BIBO) stability criterion of the system. The new model, entitled MIMO ARX-Laguerre multiple models, ensures the parameter number reduction with a recursive and easy representation. The proposed model is characterized by a set of poles where an optimal choice of these poles is compulsory to lessen considerably the number of MIMO ARX-Laguerre multiple model parameters. In this paper, we propose to use a recursive method to identify the Fourier coefficients and a metaheuristic algorithm to optimize the MIMO ARX-Laguerre multiple model poles.

The synthesis of a MIMO nonlinear fault tolerant control via MPC (MIMO NFTC-MPC) requires fault detection and estimation procedures. About that, we propose in this paper to use the moving horizon fault estimation (MHE) [19] based on the proposed MIMO ARX-Laguerre multiple models. The MHE is used to estimate the actuator faults from the error between the estimated outputs and the system outputs. The main contributions of this paper is basically threefold. (1) We present a new reduced complexity model for nonlinear MIMO systems by expanding the MIMO ARX multimodel on independent Laguerre bases. The resulting MIMO ARX-Laguerre multiple models ensures the parameter number reduction with a recursive and easy representation. (2) We develop a fault actuator detection and estimation based on the identified MIMO ARX-Laguerre multimodel and using the moving horizon fault estimation. (3) By combining the fault estimation procedure and the model predictive control for MIMO nonlinear system based on the MIMO ARX-Laguerre multiple models, we develop a MIMO nonlinear fault tolerant control via model predictive control.

This paper is organized as follows: in Section 2, we present the modeling of MIMO nonlinear systems where we recall the principle of MIMO multiple model approach and we give the definition of the MIMO ARX multiple models. In Section 3, we present the MIMO ARX-Laguerre multiple models obtained by the expansion of MIMO ARX multiple models on independent Laguerre bases. In Section 4, we propose the identification procedure of the MIMO ARX-Laguerre multiple models, where we develop a recursive method to identify the Fourier coefficients and we use a metaheuristic algorithm to optimize the poles. Section 5 is devoted to the development of a MHE to detect and estimate the actuator faults of the MIMO nonlinear system using the MIMO ARX-Laguerre multiple models. In Section 6, we synthesize the MIMO nonlinear fault tolerant control via MPC where we develop the j-step ahead predictor of the MIMO ARX-Laguerre multiple model outputs by taking into account the actuator faults and we present the control calculation by taking into account the constraint on the inputs and the outputs by resolving an optimization problem. Finally, Section 7 illustrates the proposed MIMO nonlinear fault tolerant control via MPC by a numerical example.

2. Modeling of MIMO Nonlinear Systems

2.1. Principle of Multiple Model Approach

A multiple model is a set of LTI (linear time invariant) and causal submodels aggregated by an interpolation mechanism to characterize the dynamic behavior of the overall nonlinear system. It is characterized by the number of submodels, their structure, and the choice of weighting functions. A multiple model structure is represented by$$\begin{array}{}\text{(1)}& y\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\gamma \left(k-1\right)\right)y_s\left(k\right),\end{array}$$where $y\left(k\right)$ is the multiple model output, L is the submodel number, ${\mu }_s\left(\gamma \left(k\right)\right)$ is the weighting function associated to the $s^{\text{th}}$ submodel, $\gamma \left(k\right)$ is the decision variable in general is selected as the input or the output of the system, and $y_s\left(k\right)$ is the output of the $s^{\text{th}}$ submodel. The weighting functions ${\mu }_s\left(\gamma \left(k\right)\right)$ allow to determine the relative contribution of each submodel according to the zone where the system operates, and they respect the convexity properties given as follow:$$\begin{array}{}\text{(2)}& \sum _{s=1}^L{\mu }_s\left(\gamma \left(k\right)\right)=1,\hspace{1em}0\le {\mu }_s\left(\gamma \left(k\right)\right)\le 1\forall s=1,\dots ,L.\end{array}$$

The weighting functions can be constructed from continuous functions derivatives such as Gaussian functions as follows:$$\begin{array}{}\text{(3)}& \left\{\begin{array}{c}{w}_s\left(\gamma \left(k\right)\right)=\exp \left(-\frac{{\left(\gamma \left(k\right)-c_s\right)}^2}{{\sigma }_s^2}\right),\\ {\mu }_s\left(\gamma \left(k\right)\right)=\frac{w_s\left(\gamma \left(k\right)\right)}{{\sum }_{j=1}^L{w}_j\left(\gamma \left(k\right)\right)},\end{array}\right.\end{array}$$where ${\sigma }_s$ and $c_s$ are, respectively, the dispersion and the center of the indexed variable $\gamma \left(k\right)$.

2.2. MIMO ARX Multiple Models

A strictly causal discrete time MIMO nonlinear system with p inputs and m outputs can be represented by a MIMO ARX multiple models where each nonlinear multiple-input single-output (MISO) system $y_i\left(k\right)$ for $\left(i=1,\dots ,m\right)$ can be represented by a MISO ARX multiple models written as$$\begin{array}{}\text{(4)}& y_i\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right)y_i^s\left(k\right),\end{array}$$where $y_i^s\left(k\right)$ is the output of the $s^{\text{th}}$ submodel for the MISO model, ${\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right)$ is the weighting function associated with the s^th submodel for the MISO model, and $\underset{¯}{\gamma }\left(k-1\right)$ is the decision variable selected as the inputs or the outputs of the system:$$\begin{array}{}\text{(5)}& \underset{¯}{\gamma }\left(k-1\right)=\left[{\gamma }_1\left(k-1\right),{\gamma }_2\left(k-1\right),\dots ,{\gamma }_h\left(k-1\right)\right],\end{array}$$where ${\gamma }_i\left(k-1\right)=u_i\left(k-1\right)$ for $\left(i=1,\dots ,p\right)$, $\left(h=p\right)$, if the decision variable is selected as the inputs and ${\gamma }_i\left(k-1\right)=y_i\left(k-1\right)$ for $\left(i=1,\dots ,m\right)$, $\left(h=m\right)$, if the decision is variable selected as the outputs. By defining the following matrices,$$\begin{array}{c}\text{(6)}& C=\left[\begin{array}{cccc}{c}_{11}& c_{12}& \cdots & c_{1L}\\ c_{21}& c_{22}& \cdots & c_{2L}\\ ⋮& ⋮& ⋮& ⋮\\ c_{h1}& c_{h2}& \cdots & c_{hL}\end{array}\right]\in {ℝ}^{\left(L\times h\right)},\sigma =\left[\begin{array}{cccc}{\sigma }_1& {\sigma }_2& \cdots & {\sigma }_h\end{array}\right],\end{array}$$where $c_{ij}$ and ${\sigma }_i$ are, respectively, the $j^{\text{th}}$ center, $\left(j=1,\dots ,L\right)$, and the dispersion of the $i^{\text{th}}$ decision variable, $\left(i=1,\dots ,h\right)$, where ${\sigma }_i$ is defined as$$\begin{array}{}\text{(7)}& {\sigma }_i=\sqrt{\frac{1}{M}\sum _{k=1}^M\left({\gamma }_i\left(k\right)-\overline{{\gamma }_i}\right).}\end{array}$$

Then, the weighting function ${\mu }_s\left(\underset{¯}{\gamma }\left(k\right)\right)$ is defined as$$\begin{array}{}\text{(8)}& \left\{\begin{array}{c}{w}_s\left(\underset{¯}{\gamma }\left(k\right)\right)=\prod _{j=1}^h\exp \left(-\frac{{\left({\gamma }_j\left(k\right)-C_{js}\right)}^2}{{\sigma }_j^2}\right),\\ {\mu }_s\left(\underset{¯}{\gamma }\left(k\right)\right)=\frac{w_s\left(\underset{¯}{\gamma }\left(k\right)\right)}{{\sum }_{j=1}^L{w}_j\left(\underset{¯}{\gamma }\left(k\right)\right)}.\end{array}\right.\end{array}$$

The $s^{\text{th}}$ MISO linear submodel can be described by its output equation $y_i^s\left(k\right)$ given by a MISO autoregressive with exogenous input (MISO ARX) model as follows:$$\begin{array}{}\text{(9)}& y_i^s\left(k\right)=\sum _{t=1}^p\sum _{j=1}^{n_a}{h}_{a_{it}}^s\left(j\right)u_t\left(k-j\right)+\sum _{i=1}^m\sum _{j=1}^{n_b}{h}_{b_{ir}}^s\left(j\right)y_i\left(k-j\right),\end{array}$$where $h_{a_{it}}^s\left(j\right)$ and $h_{b_{ir}}^s\left(j\right)$ are the model parameters of the $s^{\text{th}}$ MISO ARX submodel, $u_t\left(k\right)$ and $y_i\left(k\right)$ for $\left(t=1,\dots ,p\right)$ and $\left(i=1,\dots ,m\right)$ are, respectively, the $t^{\text{th}}$ system inputs and the $i^{\text{th}}$ system outputs, and $n_a$ and $n_b$ are the model orders associated, respectively, with the inputs and the outputs of every $s^{\text{th}}$ MISO ARX submodel. The $s^{\text{th}}$ MISO ARX submodel can be written in the matrix form as$$\begin{array}{}\text{(10)}& y_i^s\left(k\right)={\phi }^T\left(k\right){\theta }_i^s,\end{array}$$with$$\begin{array}{c}\text{(11)}& \phi \left(k\right)={\left[u_1\left(k-1\right),\dots ,u_1\left(k-n_a\right),\dots ,u_p\left(k-1\right),\dots ,u_p\left(k-n_a\right),y_1\left(k-1\right),\dots ,y_1\left(k-n_b\right),\dots ,y_m\left(k-1\right),\dots ,y_m\left(k-n_b\right)\right]}^T,{\theta }_i^s={\left[h_{a_{i1}}^s\left(1\right),\dots ,h_{a_{i1}}^s\left(n_a\right),\dots ,h_{a_{ip}}^s\left(1\right),\dots ,h_{a_{ip}}^s\left(n_a\right),h_{b_{i1}}^s\left(1\right),\dots ,h_{b_{i1}}^s\left(n_b\right),\dots ,h_{b_{ip}}^s\left(1\right),\dots ,h_{b_{ip}}^s\left(n_b\right)\right]}^T,\end{array}$$then, the MISO ARX multiple models given by relation (4) can be rewritten as$$\begin{array}{}\text{(12)}& y_i\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\gamma \left(k-1\right)\right){\phi }^T\left(k\right){\theta }_i^s.\end{array}$$

The output of the MISO ARX multiple models $y_i\left(k\right)$ can be rewritten as$$\begin{array}{}\text{(13)}& y_i\left(k\right)={\psi }^T\left(k\right){\Theta }_i,\end{array}$$with$$\begin{array}{c}\text{(14)}& \psi \left(k\right)={\left[\begin{array}{c}{\mu }_1\left(\gamma \left(k-1\right)\right)\phi \left(k\right)\\ ⋮\\ {\mu }_L\left(\gamma \left(k-1\right)\right)\phi \left(k\right)\end{array}\right]}^T,{\Theta }_i=\left[\begin{array}{c}{\theta }_i^1\\ ⋮\\ {\theta }_i^L\end{array}\right],\end{array}$$where each MISO ARX multiple models is characterized by a $N_i$ parameters number determined as follows:$$\begin{array}{}\text{(15)}& N_i=L\left(p.n_a+m.n_b\right).\end{array}$$

Then, the MIMO ARX multiple models are characterized by a $\left(m.N_i\right)$ parameter number. From relation (15), we can conclude that the complexity of the MIMO ARX multiple models increases according to the submodel orders $n_a$ and $n_b$. In order to reduce the number of parameters, we will proceed with the decomposition of these coefficients $h_{a_{it}}^s\left(j\right)$ and $h_{b_{ir}}^s\left(j\right)$ for $\left(i,r=1,\dots ,m\text{\hspace{0.17em}and\hspace{0.17em}}t=1,\dots ,p\right)$ of the MIMO ARX multiple models given by (9), for $\left(i=1,\dots ,m\right)$, on Laguerre orthonormal bases.

3. Expansion of MIMO ARX Multiple Models on Laguerre Orthonormal Bases

In this section, we use the Laguerre orthonormal bases to reduce the parametric complexity of the MIMO ARX multiple models defined by (9) [15, 18, 20]. This choice is due to the capability of Laguerre base on parametric reduction and for the classical recurrent representation. According to the stability condition of the system in the sense of bounded-input bounded-output criterion (BIBO), the coefficients $h_{a_{it}}^s\left(j\right)$ and $h_{b_{ir}}^s\left(j\right)$ are absolutely summable and they satisfy$$\begin{array}{}\text{(16)}& \left\{\begin{array}{c}\sum _{j=1}^{\infty }\left|h_{a_{it}}^s\left(j\right)\right|<\infty ,\hspace{1em}\left(t=1,\dots ,p\right),\\ \sum _{j=1}^{\infty }\left|h_{b_{ir}}^s\left(j\right)\right|<\infty ,\hspace{1em}\left(r=1,\dots ,m\right),\\ \left(s=1,\dots ,L\right),\left(i=1,\dots ,m\right).\end{array}\right.\end{array}$$

Therefore, these coefficients belong to the Lebesgue space ${\ell }^2\left[0,+\infty \right)$. Noting that the orthogonal Laguerre functions form an orthogonal base belong also to the Lebesgue space, the coefficients $h_{a_{it}}^s\left(j\right)$ and $h_{b_{ir}}^s\left(j\right)$ can be, respectively, developed on the Laguerre bases ${\Im }_{a_{it}}^s={\left\{{\ell }_{a_{it}}^s\left(j\right)\right\}}_{j=0}^{\infty }$ and ${\Im }_{b_{ir}}^s={\left\{{\ell }_{b^{ir}}^s\left(j\right)\right\}}_{j=0}^{\infty }$ as follows [20, 21]:$$\begin{array}{}\text{(17)}& \left\{\begin{array}{c}{h}_{a_{it}}^s\left(j\right)=\sum _{n=0}^{\infty }{g}_{a_{it},n}^s{\ell }_{a_{it},n}^s\left({\xi }_{a_{it}}^s,j\right),\hspace{1em}\left(t=1,\dots ,p\right),\\ \begin{array}{c}{h}_{b_{ir}}^s\left(j\right)=\sum _{n=0}^{\infty }{g}_{b_{ir},n}^s{\ell }_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,j\right),\hspace{1em}\left(r=1,\dots ,m\right),\\ \left(i=1,\dots ,m\right),\left(s=1,\dots ,L\right),\end{array}\end{array}\right.\end{array}$$where for $\left(s=1,\dots ,L\right),\left(t=1,\dots ,p\right),\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}}\left(i,r=1,\dots ,m\right),$ $g_{a_{it},n_a}^s$ and $g_{b_{ir},n_b}^s$ are the Fourier coefficients, ${\ell }_{a_{it},n}^s\left({\xi }_{a_{it}}^s,j\right)$ and i${\ell }_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,j\right)$ are the Laguerre functions, and ${\xi }_{a_{it}}^s$ and ${\xi }_{b_{ir}}^s$ are the poles defining, respectively, the orthogonal bases ${\Im }_{a_{it}}^s$ and ${\Im }_{b_{ir}}^s$. Taking into account the stability condition (16), the $s^{\text{th}}$ MISO ARX multiple models given by (9) can be written as$$\begin{array}{}\text{(18)}& y_i^s\left(k\right)=\sum _{t=1}^p\sum _{j=1}^{\infty }{h}_{a_{it}}^s\left(j\right)u_t\left(k-j\right)+\sum _{r=1}^m\sum _{j=1}^{\infty }{h}_{b_{ir}}^s\left(j\right)y_r\left(k-j\right),\end{array}$$where $h_{a_{it}}^s\left(j\right)=0$ if $j>n_a$ and $h_{b_{ir}}^s\left(j\right)=0$ if $j>n_b$.

By substituting $h_{a_{it}}^s\left(j\right)$ and $h_{b_{ir}}^s\left(j\right)$ given by (17) in the $s^{\text{th}}$ MISO ARX submodel defined by (18), the resulting submodel can be written as$$\begin{array}{}\text{(19)}& y_i^s\left(k\right)=\sum _{t=1}^p\sum _{j=1}^{\infty }\sum _{n=0}^{\infty }{g}_{a_{it},n}^s{\ell }_{a_{it},n}^s\left({\xi }_{a_{it}}^s,j\right)u_t\left(k-j\right)+\sum _{r=1}^m\sum _{j=1}^{\infty }\sum _{n=0}^{\infty }{g}_{b_{ir},n}^s{\ell }_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,j\right)y_r\left(k-j\right),\end{array}$$for $\left(s=1,\dots ,L\right)\text{\hspace{0.17em}and\hspace{0.17em}}\left(i=1,\dots ,m\right)$.

The relation (19) can be written, for $\left(s=1,\dots ,L\right)\text{and\hspace{0.17em}}\left(i=1,\dots ,m\right)$, as$$\begin{array}{}\text{(20)}& y_i^s\left(k\right)=\sum _{t=1}^p\sum _{n=0}^{\infty }{g}_{a_{it},n}^s\left(\sum _{j=1}^{\infty }{\ell }_{a_{it},n}^s\left({\xi }_{a_{it}}^s,j\right)u_t\left(k-j\right)\right)+\sum _{r=1}^m\sum _{n=0}^{\infty }{g}_{b_{ir},n}^s\left(\sum _{j=1}^{\infty }{\ell }_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,j\right)y_r\left(k-j\right)\right).\end{array}$$

By analogy to the development given by Mbarek et al. [22] for the MIMO case and from relation (A.3) and (A.4) given in Appendix A, the MISO ARX-Laguerre multiple model can be described by the following recursive representation:$$\begin{array}{}\text{(21)}& \left\{\begin{array}{c}{X}_i^s\left(k+1\right)=A_i^s{X}_i^s\left(k\right)+B_{a_i}^s\underset{¯}{u}\left(k\right)+B_{b_i}^s\underset{¯}{y}\left(k\right),\\ y_i^s\left(k\right)={\left(C_i^s\right)}^T{X}_i^s\left(k\right),\\ y_i\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right)y_i^s\left(k\right),\end{array}\right.\end{array}$$where $C_i^s$, for $\left(i=1,\dots ,m\right)$ and $\left(s=1,\dots ,L\right)$, are the parameter vectors regrouping the Fourier coefficients $g_{a_{it},n}^s$, for $\left(n=0,\dots ,N_a^s-1\right)$ and $\left(t=1,\dots ,p\right)$, and $g_{b_{ir},n}^s$, for $\left(n=0,\dots ,N_b^s-1\right)$ and $\left(r=1,\dots ,m\right)$.$$\begin{array}{}\text{(22)}& C_i^s={\left[g_{a_{i1},0}^s,\dots ,g_{a_{i1},N_a^s-1}^s,g_{a_{i2},0}^s,\dots ,g_{a_{i2},N_a^s-1}^s,\dots ,g_{a_{ip},0}^s,\dots ,g_{a_{ip},N_a^s-1}^s,g_{b_{i1},0}^s,\dots ,g_{b_{i1},N_b^s-1}^s,g_{b_{i2},0}^s,\dots ,g_{b_{i2},N_b^s-1}^s,\dots ,g_{b_{im},0}^s,\dots ,g_{b_{im},N_b^s-1}^s\right]}^T,\end{array}$$and the matrices $A_i^s$ are defined by$$\begin{array}{}\text{(23)}& A_i^s=\left[\begin{array}{cccccc}{A}_{a_{i1}}^s& \cdots & {\underset{¯}{0}}_{N_a^s\times N_a^s}& {\underset{¯}{0}}_{N_a^s\times N_b^s}& \cdots & {\underset{¯}{0}}_{N_a^s\times N_b^s}\\ ⋮& \ddots & ⋮& ⋮& ⋮& ⋮\\ {\underset{¯}{0}}_{N_i^s\times N_i^s}& \cdots & A_{a_{ip}}^s& {\underset{¯}{0}}_{N_a^s\times N_b^s}& \cdots & {\underset{¯}{0}}_{N_a^s\times N_b^s}\\ {\underset{¯}{0}}_{N_b^s\times N_a^s}& \cdots & {\underset{¯}{0}}_{N_b^s\times N_a^s}& A_{b_{i1}}^s& \cdots & {\underset{¯}{0}}_{N_b^s\times N_b^s}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ {\underset{¯}{0}}_{N_b^s\times N_a^s}& \cdots & {\underset{¯}{0}}_{N_b^s\times N_a^s}& {\underset{¯}{0}}_{N_b^s\times N_b^s}& \cdots & A_{b_{im}}^s\end{array}\right],\end{array}$$where the matrices $A_{a_{it}}^s$, $\left(t=1,\dots ,p\right)$ and $A_{b_{ir}}^s$, $\left(r=1,\dots ,m\right)$ are defined according to the poles ${\xi }_{a_{it}}^s$ and ${\xi }_{b_{ir}}^s$ as given, respectively, by relations (A.13) and (A.14).

The matrices $B_{a_i}^s$ and $B_{b_i}^s$ are defined by$$\begin{array}{c}\text{(24)}& B_{a_i}^s=\left[\begin{array}{c}{B}_{a_i}^{s^{\prime }}\\ {\underset{¯}{0}}_{\left(p.N_a\right)\times p}\end{array}\right],B_{b_i}^s=\left[\begin{array}{c}{\underset{¯}{0}}_{\left(m.N_b\right)\times m}\\ B_{b_i}^{s^{\prime }}\end{array}\right],\end{array}$$where $B_{a_i}^{s^{\prime }}$ and $B_{b_i}^{s^{\prime }}$ are defined as follows:$$\begin{array}{c}\text{(25)}& B_{a_i}^{s^{\prime }}=\left[\begin{array}{cccc}{b}_{a_{i1}}^s& {\underset{¯}{0}}_{N_a^s\times 1}& \cdots & {\underset{¯}{0}}_{N_a^s\times 1}\\ {\underset{¯}{0}}_{N_a^s\times 1}& b_{a_{i2}}^s& \dots & {\underset{¯}{0}}_{N_a^s\times 1}\\ ⋮& ⋮& \ddots & ⋮\\ {\underset{¯}{0}}_{N_a^s\times 1}& {\underset{¯}{0}}_{N_a^s\times 1}& \cdots & b_{a_{ip}}^s\end{array}\right],B_{b_i}^{s^{\prime }}=\left[\begin{array}{cccc}{b}_{b_{i1}}^s& {\underset{¯}{0}}_{N_b^s\times 1}& \cdots & {\underset{¯}{0}}_{N_b^s\times 1}\\ {\underset{¯}{0}}_{N_b^s\times 1}& b_{b_{i2}}^s& \cdots & {\underset{¯}{0}}_{N_b^s\times 1}\\ ⋮& ⋮& \ddots & ⋮\\ {\underset{¯}{0}}_{N_b^s\times 1}& {\underset{¯}{0}}_{N_b^s\times 1}& \cdots & b_{b_{im}}^s\end{array}\right],\end{array}$$where the vectors $B_{a_{it}}^s$, $\left(t=1,\dots ,p\right)$ and $B_{b_{ir}}^s$, $\left(r=1,\dots ,m\right)$ are defined according to the poles ${\xi }_{a_{it}}^s$ and ${\xi }_{b_{ir}}^s$ as given, respectively, by relations (A.15) and (A.16).

The output of the $i^{\text{th}}$ MISO ARX-Laguerre multiple models $y_i\left(k\right)$ given by (21) can be written as$$\begin{array}{}\text{(26)}& y_i\left(k\right)={\underset{¯}{C}}_i^T{X}_i\left(k\right),\end{array}$$where ${\underset{¯}{C}}_i$ and $X_i\left(k\right)$ are defined as follows:$$\begin{array}{}\text{(27)}& {\underset{¯}{C}}_i={\left[\begin{array}{cccc}{\left(C_i^1\right)}^T& {\left(C_i^2\right)}^T& \cdots & {\left(C_i^L\right)}^T\end{array}\right]}^T,\text{(28)}& X_i\left(k\right)={\left[\begin{array}{cccc}{\mu }_1\left(\underset{¯}{\gamma }\left(k-1\right)\right){\left(X_i^1\left(k\right)\right)}^T& {\mu }_2\left(\underset{¯}{\gamma }\left(k-1\right)\right){\left(X_i^2\left(k\right)\right)}^T& \cdots & {\mu }_L\left(\underset{¯}{\gamma }\left(k-1\right)\right){\left(X_i^L\left(k\right)\right)}^T\end{array}\right]}^T.\end{array}$$

The MISO ARX-Laguerre multiple models given by relation (21) is characterized by $L.\left(p.N_a+m.N_b\right)$ parameters as given by relation (27) and $\left(\left(m+p\right).L\right)$ poles defined as follows:$$\begin{array}{}\text{(29)}& {\underset{¯}{\xi }}_i={\left[{\xi }_{a_{i1}}^1,\dots ,{\xi }_{a_{ip}}^1,{\xi }_{b_{i1}}^1,\dots ,{\xi }_{b_{im}}^1,{\xi }_{a_{i1}}^2,\dots ,{\xi }_{a_{ip}}^2,{\xi }_{b_{i1}}^2,\dots ,{\xi }_{b_{im}}^2,\dots ,{\xi }_{a_{i1}}^L,\dots ,{\xi }_{a_{ip}}^L,{\xi }_{b_{i1}}^L,\dots ,{\xi }_{b_{im}}^L\right]}^T.\end{array}$$

The proposed MIMO ARX-Laguerre multiple models given, for $\left(i=1,\dots ,m\right)$, by the MISO ARX-Laguerre multiple modes (21) can be represented, as in Figure 1, by a parallel structure diagram in the case of 2 inputs/2 outputs where every MISO model is decomposed into L submodels.

[figure omitted; refer to PDF]

By defining the output vectors $\underset{¯}{y}$ and ${\underset{¯}{y}}_s$ as$$\begin{array}{c}\text{(30)}& \underset{\_}{y}\left(k\right)={\left[y_{\mathrm{1}}\left(k\right),y_{\mathrm{2}}\left(k\right),\dots ,y_m\left(k\right)\right]}^T,{\underset{\_}{y}}_s\left(k\right)={\left[y_1^s\left(k\right),y_2^s\left(k\right),\dots ,y_m^s\left(k\right)\right]}^T,\end{array}$$the MIMO ARX-Laguerre multiple models given by relation (21) can be written as$$\begin{array}{}\text{(31)}& \left\{\begin{array}{c}{\mathbf{X}}_s\left(k+1\right)={\mathbf{A}}_s{\mathbf{X}}_s\left(k\right)+{\mathbf{B}}_{a,s}\underset{¯}{U}\left(k\right)+{\mathbf{B}}_{b,s}\underset{¯}{y}\left(k\right),\\ {\underset{¯}{y}}_s\left(k\right)={\mathbf{C}}_s{\mathbf{X}}_s\left(k\right),\\ \underset{¯}{y}\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right){\underset{¯}{y}}_s\left(k\right),\end{array}\right.\end{array}$$with$$\begin{array}{c}\text{(32)}& {\mathbf{X}}_s\left(k\right)={\left[{\left(X_1^s\left(k\right)\right)}^T,{\left(X_2^s\left(k\right)\right)}^T,\dots ,{\left(X_m^s\left(k\right)\right)}^T\right]}^T,{\mathbf{C}}_s=\left[\begin{array}{cccc}{\left(C_1^s\right)}^T& 0& \cdots & 0\\ 0& {\left(C_2^s\right)}^T& \cdots & 0\\ ⋮& ⋮& \ddots & 0\\ 0& 0& \cdots & {\left(C_m^s\right)}^T\end{array}\right],\\ {\mathbf{A}}_s=\left[\begin{array}{cccc}{A}_1^s& 0& \cdots & 0\\ 0& A_2^s& \cdots & 0\\ ⋮& ⋮& \ddots & 0\\ 0& 0& \cdots & A_m^s\end{array}\right],\\ {\mathbf{B}}_{a,s}=\left[\begin{array}{c}{B}_{a_1}^s\\ B_{a_2}^s\\ ⋮\\ B_{a_m}^s\end{array}\right],\\ {\mathbf{B}}_{b,s}=\left[\begin{array}{c}{B}_{b_1}^s\\ B_{b_2}^s\\ ⋮\\ B_{b_m}^s\end{array}\right].\end{array}$$

Then, the MIMO ARX-Laguerre multiple models can be written in the following recursive representation:$$\begin{array}{}\text{(33)}& \left\{\begin{array}{c}\mathbf{X}\left(k+1\right)=\mathbf{A}\mathbf{X}\left(k\right)+{\mathbf{B}}_a\underset{¯}{u}\left(k\right)+{\mathbf{B}}_b\underset{¯}{y}\left(k\right),\\ \underset{¯}{y}\left(k\right)=\mathbf{C}\mathbf{X}\left(k\right),\end{array}\right.\end{array}$$where the state vector $\mathbf{X}\left(k\right)$ and the matrices $\mathbf{C}$, $\mathbf{A}$, ${\mathbf{B}}_a$, and ${\mathbf{B}}_b$ are defined, respectively, as$$\begin{array}{c}\text{(34)}& \mathbf{X}\left(k\right)={\left[{\left({\mathbf{X}}_1\left(k\right)\right)}^T,{\left({\mathbf{X}}_2\left(k\right)\right)}^T,\dots ,{\left({\mathbf{X}}_L\left(k\right)\right)}^T\right]}^T,\text{(35)}& \mathbf{C}=\left[\begin{array}{cccc}{\mu }_1\left(\underset{¯}{\gamma }\left(k-1\right)\right){\mathbf{C}}_1& {\mu }_2\left(\underset{¯}{\gamma }\left(k-1\right)\right){\mathbf{C}}_2& \dots & {\mu }_L\left(\underset{¯}{\gamma }\left(k-1\right)\right){\mathbf{C}}_L\end{array}\right],\text{(36)}& \mathbf{A}=\left[\begin{array}{cccc}{\mathbf{A}}_1& \underset{¯}{0}& \dots & \underset{¯}{0}\\ \underset{¯}{0}& {\mathbf{A}}_2& \dots & \underset{¯}{0}\\ ⋮& ⋮& \ddots & \underset{¯}{0}\\ \underset{¯}{0}& \underset{¯}{0}& \dots & {\mathbf{A}}_L\end{array}\right],{\mathbf{B}}_a=\left[\begin{array}{c}{\mathbf{B}}_{a,1}\\ {\mathbf{B}}_{a,2}\\ ⋮\\ {\mathbf{B}}_{a,L}\end{array}\right],\\ {\mathbf{B}}_b=\left[\begin{array}{c}{\mathbf{B}}_{b,1}\\ {\mathbf{B}}_{b,2}\\ ⋮\\ {\mathbf{B}}_{b,L}\end{array}\right].\end{array}$$

4. Identification of the MIMO ARX-Laguerre Multiple Models

A strictly causal discrete time nonlinear MIMO system with p inputs and m outputs can be described by a MIMO ARX-Laguerre multiple models where each MISO system, for $\left(i=1,\dots ,m\right)$, is represented by the recursive representation given by relation (21) and the $i^{\text{th}}$ output $y_i\left(k\right)$ is rewritten by the vectorial representation as given in relation (26). Each MISO ARX-Laguerre multiple model is characterized by the parameter vector ${\underset{¯}{C}}_i$ defined by relation (27) and the poles vector ${\underset{¯}{\xi }}_i$ given by relation (29). To ensure a significant parametric reduction, it is necessary to optimize the MIMO ARX-Laguerre multiple model poles. The identification of the parameter vectors ${\underset{¯}{C}}_i$ and the optimization of the poles ${\underset{¯}{\xi }}_i$ are presented in the next sections.

4.1. Recursive Identification of the Fourier Coefficients

This method proposed by Abdelwahed et al. [21] is based on the minimization of a regularized square error given as$$\begin{array}{}\text{(37)}& J_i=\sum _{k=1}^h{\left({\overline{y}}_i\left(k\right)-y_i\left(k\right)\right)}^2+{\alpha }_i\sum _{s=1}^L\left(\sum _{t=1}^p\sum _{n=1}^{N_a^s}{\left(g_{a_{it},n}^s\left(h-1\right)-g_{a_{it},n}^s\left(h\right)\right)}^2+\sum _{r=1}^m\sum _{n=1}^{N_b^s}{\left(g_{b_{ir},n}^s\left(h-1\right)-g_{b_{ir},n}^s\left(h\right)\right)}^2\right),\end{array}$$where

(i)

${\overline{y}}_i\left(k\right)$ and $y_i\left(h\right)$ are, respectively, the $i^{\text{th}}$ output system and the $i^{\text{th}}$ output of the MIMO ARX-Laguerre multiple models

(ii)

$g_{a_{it},n}^s\left(h-1\right)$, $g_{b_{ir},n}^s\left(h-1\right)$, and $g_{a_{it},n}^s\left(h\right)$, $g_{b_{ir},n}^s\left(h\right)$ are, respectively, the Fourier coefficients at time instants $\left(h-1\right)$ and h

(iii)

${\alpha }_i$, $\left({\alpha }_i>0\right)$, is a regularization constant which regulates the importance that we give either to the quadratic error between ${\overline{y}}_i\left(k\right)$ and $y_i\left(h\right)$ or to that between the Fourier coefficients.

From relation (26) and at time instant h, the square error $J_i$ can be written in matrix form as$$\begin{array}{}\text{(38)}& J_i\left(h\right)={\left({\overline{Y}}_i\left(h\right)-{\mathbf{\Psi }}_i\left(h\right){\underset{¯}{C}}_i\left(h\right)\right)}^T\left({\overline{Y}}_i\left(h\right)-{\mathbf{\Psi }}_i\left(h\right){\underset{¯}{C}}_i\left(h\right)\right)+{\left({\underset{¯}{C}}_i\left(h-1\right)-{\underset{¯}{C}}_i\left(h\right)\right)}^T{\alpha }_i\left({\underset{¯}{C}}_i\left(h-1\right)-{\underset{¯}{C}}_i\left(h\right)\right),\end{array}$$where ${\underset{¯}{C}}_i\left(h-1\right)$ and ${\underset{¯}{C}}_i\left(h\right)$ are the parameter vectors regrouping the Fourier coefficients associated with the $i^{\text{th}}$ output at time instant $\left(h-1\right)$ and h, respectively, and the vector ${\overline{Y}}_i\left(h\right)$ and the matrix ${\mathbf{\Psi }}_i\left(h\right)$ are defined as$$\begin{array}{c}\text{(39)}& {\overline{Y}}_i\left(h\right)=\left[\begin{array}{c}{\overline{y}}_i\left(1\right)\\ {\overline{y}}_i\left(2\right)\\ ⋮\\ {\overline{y}}_i\left(h\right)\end{array}\right],{\mathbf{\Psi }}_i\left(h\right)=\left[\begin{array}{c}{X}_i{\left(1\right)}^T\\ X_i{\left(2\right)}^T\\ ⋮\\ X_i{\left(h\right)}^T\end{array}\right]=\left[\begin{array}{c}{\mathbf{\Psi }}_i\left(h-1\right)\\ X_i{\left(h\right)}^T\end{array}\right].\end{array}$$

The optimal parameter vector ${\underset{¯}{C}}_i\left(h\right)$ of each $i^{\text{th}}$ MISO system is obtained by minimizing the quadratic regularised square error $J_i\left(h\right)$ given by relation (38). The gradient of the criterion $J_i\left(h\right)$ with respect to ${\underset{¯}{C}}_i\left(h\right)$ is written as$$\begin{array}{}\text{(40)}& \frac{\partial J_i\left(h\right)}{\partial {\underset{¯}{C}}_i\left(h\right)}=-2\left[{\mathbf{\Psi }}_i\left(h\right){\overline{Y}}_i\left(h\right)-{\mathbf{\Psi }}_i{\left(h\right)}^T{\mathbf{\Psi }}_i\left(h\right){\underset{¯}{C}}_i\left(h\right)\right]-2{\alpha }_i\left({\underset{¯}{C}}_i\left(h-1\right)-{\underset{¯}{C}}_i\left(h\right)\right).\end{array}$$

From relation (40), we can calculate the estimated parameter vector ${\underset{¯}{C}}_i\left(h\right)$ of each $i^{\text{th}}$ MISO system as$$\begin{array}{}\text{(41)}& {\underset{¯}{C}}_i\left(h\right)={\left({\mathbf{\Psi }}_{\mathbf{i}}{\left(h\right)}^T{\mathbf{\Psi }}_{\mathbf{i}}\left(h\right)+{\alpha }_{i}I\right)}^{-1}\left({\mathbf{\Psi }}_{\mathbf{i}}\left(h\right){\overline{Y}}_i\left(h\right)+{\alpha }_i{\underset{¯}{C}}_i\left(h-1\right)\right).\end{array}$$

The recursive parametric identification of the MIMO ARX-Laguerre multiple models representing a nonlinear MIMO system of p inputs and m outputs where each MISO system is decomposed into L submodel is described by the following algorithm.

4.2. Pole Optimization of the MIMO ARX-Laguerre Multiple Models Using Metaheuristic Algorithms

To guarantee a significant parametric reduction of the proposed MIMO ARX-Laguerre multiple models, it is necessary to optimize the $\left(L.\left(m+p\right)\right)$ Laguerre bases poles. Several optimization methods are proposed in the literature in order to put the Laguerre poles into their optimal values, such as the gradient algorithm, the Newton–Raphson algorithm, and the Bouzrara el al.’s method. In this paper, since we have many poles to be optimized, $\left(L.\left(m+p\right)\right)$ poles, we propose the use of genetic algorithms (GAs) as a metaheuristic optimization method. The values of the MIMO ARX-Laguerre multiple models poles are optimized by minimizing the normalized mean squared error $\left({\text{NMSE}}_i\right)$ of each MISO ARX-Laguerre multiple models given by$$\begin{array}{}\text{(42)}& {\text{NMSE}}_i=\frac{{\sum }_{k=1}^M{\left({\overline{y}}_i\left(k\right)-y_i\left(k\right)\right)}^2}{{\sum }_{k=1}^M{\left({\overline{y}}_i\left(k\right)\right)}^2},\end{array}$$where ${\overline{y}}_i\left(k\right)$ and $y_i\left(h\right)$ are the $i^{\text{th}}$ output system and the $i^{\text{th}}$ output of the MIMO ARX-Laguerre multiple models, respectively.

The pole optimization algorithm of the proposed MIMO ARX-Laguerre multiple models using the genetic algorithm is summarized in the following algorithm.

5. Moving Horizon Fault Estimation Based on MIMO ARX-Laguerre Multiple Models

In the following, we propose to use the moving horizon fault estimation (MHE) based on the MIMO ARX-Laguerre multiple models to solve the fault estimation of the MIMO nonlinear system. The MHE is used to estimate the actuator faults from the error between the estimated output $\underset{¯}{\widehat{y}}\left(k\right)$ and the system output $\overline{y}\left(k\right)$. The MIMO ARX-Laguerre multiple models taking account of actuator faults can be written from (33) as$$\begin{array}{}\text{(43)}& \left\{\begin{array}{c}\widehat{X}\left(k+1\right)=\mathbf{A}\widehat{X}\left(k\right)+{\mathbf{B}}_a\left[\underset{¯}{u}\left(k\right)+\underset{¯}{f}\left(k\right)\right]+{\mathbf{B}}_b\widehat{\underset{¯}{y}}\left(k\right),\\ \widehat{\underset{¯}{y}}\left(k\right)=\mathbf{C}\widehat{X}\left(k\right),\end{array}\right.\end{array}$$where $\widehat{X}\left(k\right)$ and $\widehat{\underset{¯}{y}}\left(k\right)$ are the estimated state vector and the estimation outputs model, respectively, when the actuator faults $\underset{¯}{f}\left(k\right)$ are defined as$$\begin{array}{}\text{(44)}& \underset{¯}{f}\left(k\right)={\left[f_1\left(k\right),f_2\left(k\right),\dots ,f_p\left(k\right)\right]}^T.\end{array}$$

The estimate actuators faults $\widehat{\underset{¯}{f}}\left(k\right)$ are obtained by an online minimization at every sample time of the quadratic cost function $J_f\left(k\right)$ defined as follows:$$\begin{array}{}\text{(45)}& J_f\left(k\right)=\sum _{j=k-N_f+1}^k{\left(\widehat{\underset{¯}{y}}\left(j\right)-\overline{y}\left(j\right)\right)}^2,\end{array}$$where $N_f$ is the estimation fault horizon, $\widehat{\underset{¯}{y}}\left(j\right)$ is the estimation outputs model when the actuator faults are considered as given by (43), and $\overline{y}\left(j\right)$ is the outputs system defined as$$\begin{array}{}\text{(46)}& \overline{y}\left(j\right)={\left[{\overline{y}}_1\left(1\right),{\overline{y}}_2\left(j\right),\dots ,{\overline{y}}_m\left(j\right)\right]}^T.\end{array}$$

The quadratic criterion $J_f\left(k\right)$ given by (45) can be rewritten as$$\begin{array}{}\text{(47)}& J_f\left(k\right)=\sum _{j=k-N_f+1}^k{\left(C\left(\mathbf{A}\widehat{X}\left(j-1\right)+{\mathbf{B}}_a\left[\underset{¯}{u}\left(j-1\right)+\widehat{\underset{¯}{f}}\left(j\right)\right]+{\mathbf{B}}_b\widehat{\underset{¯}{y}}\left(j-1\right)\right)-\overline{y}\left(j\right)\right)}^2,\end{array}$$where $\widehat{\underset{¯}{f}}$ is the fault actuators to be optimized using the MHE method:$$\begin{array}{}\text{(48)}& \widehat{\underset{¯}{f}}\left(j\right)={\left[\begin{array}{cccc}{\widehat{f}}_1\left(j\right)& {\widehat{f}}_2\left(j\right)& \cdots & {\widehat{f}}_p\left(j\right)\end{array}\right]}^T.\end{array}$$

The MHE method is formulated by minimizing the criterion $J_f\left(k\right)$ given by (47) at every sample time:$$\begin{array}{}\text{(49)}& \widehat{F}\left(k-1\right)=\underset{F\left(k-1\right)\in S_f}{\min }{J}_f\left(k\right),\end{array}$$where $\widehat{F}\left(k-1\right)$ is defined as follows:$$\begin{array}{}\text{(50)}& \widehat{F}\left(k-1\right)={\left[\begin{array}{cccc}\widehat{\underset{¯}{f}}\left(k-Nf+1\right)& \widehat{\underset{¯}{f}}\left(k-Nf+2\right)& \cdots & \widehat{\underset{¯}{f}}\left(k-1\right)\end{array}\right]}^T,\end{array}$$and $S_f$ is the admissible set of fault actuators constraints defined as$$\begin{array}{}\text{(51)}& S_f=\left\{\widehat{F}\left(k\right)/\left(F_{\min }\le \widehat{F}\left(k\right)\le F_{\max }\right)\right\},\end{array}$$where $F_{\min }$ and $F_{\max }$ are $p.N_f$-dimensional vectors defined as follows:$$\begin{array}{c}\text{(52)}& F_{\min }=\left(\begin{array}{c}{\underset{¯}{f}}_m\\ {\underset{¯}{f}}_m\\ ⋮\\ {\underset{¯}{f}}_m\end{array}\right),F_{\max }=\left(\begin{array}{c}{\underset{¯}{f}}_M\\ {\underset{¯}{f}}_M\\ ⋮\\ {\underset{¯}{f}}_M\end{array}\right).\end{array}$$where ${\underset{¯}{f}}_m$ and ${\underset{¯}{f}}_M$ are the bounds of actuator faults that can be chosen from the physical constraints defined as$$\begin{array}{c}\text{(53)}& {\underset{¯}{f}}_m=\left(\begin{array}{c}{f}_{1,m}\\ f_{2,m}\\ ⋮\\ f_{p,m}\end{array}\right),{\underset{¯}{f}}_M=\left(\begin{array}{c}{f}_{1,M}\\ f_{2,M}\\ ⋮\\ f_{p,M}\end{array}\right).\end{array}$$

6. MIMO Nonlinear Fault Tolerant Control via Model Predictive Control

In this paper, we propose a MIMO nonlinear fault tolerant control via model predictive control using MIMO ARX-Laguerre multiple models (MIMO NFTC-MPC). The proposed MIMO NFTC-MPC is applied to a nonlinear system and used to compensate online the actuator faults by including the effects of the faults in the model predictive control optimization problem. Then, the MIMO NFTC-MPC strategy requires a fault estimation unit to estimate the actuator faults and a MIMO ARX-Laguerre multiple models observer to estimate the system state. The MIMO NFTC-MPC strategy is described by a block diagram as represented in Figure 2.

[figure omitted; refer to PDF]

The MIMO NFTC-MPC is characterized by three essential steps. The first step is devoted to the calculation of the system output prediction $\underset{¯}{y}\left(k+j/k\right)$ to determine the process behavior on the prediction horizon $\left[k+1,k+N_p+1\right]$ depending on the inputs applied to the process until the instant $\left(k-1\right)T$. The second step is dedicated to the actuator fault optimization until the instant $\left(k-1\right)$, the outputs measured until the instant $kT$, and the equations describing the model, T, is the sample time. The third step is dedicated to the optimization of a criterion to determine a future control sequence over the control horizon $N_u$.

6.1. The j-Step Ahead Predictor including the Actuator Faults

The MIMO ARX-Laguerre multiple models taking into account actuator faults can be written from relation (31) as$$\begin{array}{}\text{(54)}& \left\{\begin{array}{c}{\tilde{X}}_s\left(k+1\right)={\mathbf{A}}_s{\tilde{X}}_s\left(k\right)+{\mathbf{B}}_{a,s}\left(\underset{¯}{u}\left(k\right)+\widehat{\underset{¯}{f}}\left(k\right)\right)+{\mathbf{B}}_{b,s}\widetilde{\underset{¯}{y}}\left(k\right),\\ {\widetilde{\underset{¯}{y}}}_s\left(k\right)={\mathbf{C}}_s{\tilde{X}}_s\left(k\right),\\ \widetilde{\underset{¯}{y}}\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k\right)\right){\widetilde{\underset{¯}{y}}}_s\left(k\right),\end{array}\right.\end{array}$$where $\widehat{\underset{¯}{f}}\left(k\right)$ is the vector of actuator faults estimated using the MHE method. To calculate the j-step ahead predictor on the MIMO ARX-Laguerre multiple models output taking into account actuator faults, $\tilde{y}\left(k+j/k\right)$, we use the incremental form given from (54) as$$\begin{array}{}\text{(55a)}& \delta {\tilde{X}}_s\left(k\right)={\mathbf{A}}_s\delta {\tilde{X}}_s\left(k-1\right)+{\mathbf{B}}_{a,s}\delta \underset{¯}{u}\left(k-1\right)+{\mathbf{B}}_{a,s}\delta \widehat{\underset{¯}{f}}\left(k-1\right)+{\mathbf{B}}_{b,s}\delta \underset{¯}{\tilde{y}}\left(k-1\right),\text{(55b)}& \begin{array}{c}\widetilde{\underset{¯}{y}}\left(k\right)=\widetilde{\underset{¯}{y}}\left(k-1\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right)C_i^T\delta {\tilde{X}}_s\left(k\right)+\sum _{s=1}^L\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k\right),\end{array}\end{array}$$where $\delta \underset{¯}{u}\left(k\right)$, $\delta \underset{¯}{\tilde{y}}\left(k\right)$, $\delta \widehat{\underset{¯}{f}}\left(k\right)$, $\delta {\tilde{X}}_s\left(k\right)$, and $\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k\right)\right)$ are the control increments, the output increments, the actuator faults increments, the state vector increments, and the weighting function increments, respectively, defined as$$\begin{array}{c}\text{(56)}& \delta \underset{¯}{u}\left(k\right)=\underset{¯}{u}\left(k\right)-\underset{¯}{u}\left(k-1\right);\delta \underset{¯}{\tilde{y}}\left(k\right)=\underset{¯}{\tilde{y}}\left(k\right)-\underset{¯}{\tilde{y}}\left(k-1\right);\\ \delta \widehat{\underset{¯}{f}}\left(k\right)=\widehat{\underset{¯}{f}}\left(k\right)-\widehat{\underset{¯}{f}}\left(k-1\right);\\ \delta {\tilde{X}}_s\left(k\right)={\tilde{X}}_s\left(k\right)-{\tilde{X}}_s\left(k-1\right);\\ \delta {\mu }_s\left(\underset{¯}{\gamma }\left(k\right)\right)={\mu }_s\left(\underset{¯}{\gamma }\left(k\right)\right)-{\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right).\end{array}$$

By successive substitutions and after some reorganization, we give the following proposition.

Proposition 1.

The j-step ahead predictor including the actuator faults using the MIMO ARX-Laguerre multiple models, $\widetilde{\underset{¯}{y}}\left(k+j/k\right)$, is given as$$\begin{array}{}\text{(57)}& \widetilde{\underset{¯}{y}}\left(k+j/k\right)=\widetilde{\underset{¯}{y}}\left(k\right)+\sum _{s=1}^L\sum _{l=1}^j\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k+l\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\left[K_{\left(s,j\right)}-I_{N_s}\right]\delta {\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\sum _{l=1}^j{K}_{\left(s,j-l\right)}\left({\mathbf{B}}_{a,s}\delta \underset{¯}{U}\left(k+l-1\right)+{\mathbf{B}}_{a,s}\delta \widehat{\underset{¯}{f}}\left(k+l-1\right)+{\mathbf{B}}_{b,s}\delta \underset{¯}{\tilde{y}}\left(k+l-1\right)\right),\end{array}$$where $I_{N_s}$ is the $N_s$-dimensional identity matrix, $N_s$ is the truncation order of the $s^{\text{th}}$ L submodel for each MISO submodel, and $K_{s,n}$ is defined as$$\begin{array}{c}\text{(58)}& K_{s,n}=\sum _{r=0}^n{\left({\mathbf{A}}_s\right)}^r,\hspace{1em}\text{\hspace{0.17em}for}n\ge 0,K_{s,n}=0,\hspace{1em}\text{for\hspace{0.17em}}n<0.\end{array}$$

See proof in Appendix B.

The j-step ahead predictor given by (57) is split into two components, the free and the forced components:$$\begin{array}{}\text{(59)}& \widetilde{\underset{¯}{y}}\left(k+j/k\right)={\widetilde{\underset{¯}{y}}}_{\text{fr}}\left(k+j/k\right)+{\widetilde{\underset{¯}{y}}}_{\text{fo}}\left(k+j\right),\end{array}$$where ${\widetilde{\underset{¯}{y}}}_{\text{fr}}\left(k+j/k\right)$ is the free components determined using the measured outputs until $kT$ and the controls up to $\left(k-1\right)T$:$$\begin{array}{}\text{(60)}& {\widetilde{\underset{¯}{y}}}_{\text{fr}}\left(k+j/k\right)=\widetilde{\underset{¯}{y}}\left(k\right)+\sum _{s=1}^L\sum _{l=1}^j\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k+l\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\left[K_{\left(s,j\right)}-I_{N_s}\right]\delta {\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\sum _{l=1}^j{K}_{\left(s,j-l\right)}\left({\mathbf{B}}_{a,s}\delta \widehat{\underset{¯}{f}}\left(k+l-1\right)+{\mathbf{B}}_{b,s}\delta \underset{¯}{\tilde{y}}\left(k+l-1\right)\right).\end{array}$$

The free component ${\widetilde{\underset{¯}{y}}}_{\text{fr}}\left(k+j/k\right)$ given by relation (60) can be rewritten as$$\begin{array}{}\text{(61)}& \begin{array}{c}{\widetilde{\underset{¯}{y}}}_{\text{fr}}\left(k+j/k\right)={\widetilde{\underset{¯}{y}}}_{\text{fr}1}\left(k+j/k\right)+\sum _{s=1}^L\left[{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\sum _{l=1}^j{K}_{\left(s,j-l\right)}\left({\mathbf{B}}_{a,s}\delta \widehat{\underset{¯}{f}}\left(k+l-1\right)+{\mathbf{B}}_{b,s}\delta \underset{¯}{\tilde{y}}\left(k+l-1\right)\right)\right],\end{array}\end{array}$$where ${\widetilde{\underset{¯}{y}}}_{\text{fr}1}\left(k+j/k\right)$ is defined as$$\begin{array}{}\text{(62)}& \begin{array}{c}{\widetilde{\underset{¯}{y}}}_{\text{fr}1}\left(k+j\right)=\widetilde{\underset{¯}{y}}\left(k\right)+\sum _{s=1}^L\sum _{l=1}^j\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k+l\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\left[K_{\left(s,j-l\right)}-I_{N_s}\right]\delta {\tilde{X}}_s\left(k\right)\end{array},\end{array}$$and from relation (A.25), the forced component ${\widetilde{\underset{¯}{y}}}_{\text{fo}}\left(k+j/k\right)$ resulting from the action of future control can be written as$$\begin{array}{}\text{(63)}& \begin{array}{c}{\widetilde{\underset{¯}{y}}}_{\text{fo}}\left(k+j/k\right)=\sum _{s=1}^L\left[{\mu }_s\left(\underset{¯}{\gamma }\left.k+j\right)\right){\mathbf{C}}_s\sum _{l=1}^j{K}_{\left(s,j-l\right)}{\mathbf{B}}_{a,s}\delta \underset{¯}{u}\left(k+l-1\right)\right]\end{array}.\end{array}$$

For $j=1,\dots ,N_p$, relation (59) becomes$$\begin{array}{}\text{(64)}& \tilde{Y}\left(k\right)={\tilde{Y}}_{\text{fo}}\left(k\right)+{\tilde{Y}}_{\text{fr}}\left(k\right),\end{array}$$where$$\begin{array}{}\text{(65)}& \left\{\begin{array}{c}\tilde{Y}\left(k\right)={\left[\begin{array}{ccc}\widetilde{\underset{¯}{y}}\left(k+1/k\right)& \cdots & \widetilde{\underset{¯}{y}}\left(k+N_p/k\right)\end{array}\right]}^T,\\ {\tilde{Y}}_{\text{fr}}\left(k\right)={\left[\begin{array}{ccc}{\widetilde{\underset{¯}{y}}}_{\text{fr}}\left(k+1/k\right)& \cdots & {\widetilde{\underset{¯}{y}}}_{\text{fr}}\left(k+N_p/k\right)\end{array}\right]}^T,\\ {\tilde{Y}}_{\text{fo}}\left(k\right)={\left[\begin{array}{ccc}{\widetilde{\underset{¯}{y}}}_{\text{fo}}\left(k+1/k\right)& \cdots & {\widetilde{\underset{¯}{y}}}_{\text{fo}}\left(k+N_p/k\right)\end{array}\right]}^T.\end{array}\right.\end{array}$$

From relations (71) and (73) and for $\left(j=1,\dots ,N_p\right)$, the vectors ${\tilde{Y}}_{\text{fr}}\left(k\right)$ and ${\tilde{Y}}_{\text{fo}}\left(k\right)$ can be written as$$\begin{array}{}\text{(66)}& \left\{\begin{array}{c}{\tilde{Y}}_{\text{fr}}\left(k\right)={\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right),\\ {\tilde{Y}}_{\text{fo}}\left(k\right)={\mathbf{G}}_a\delta U\left(k\right),\end{array}\right.\end{array}$$where

(i)

${\tilde{Y}}_{\text{fr}1}\left(k\right)$ is a $\left(m.N_p\right)$-dimensional vector computed from relation (72) and defined as

$$\begin{array}{}\text{(67)}& {\tilde{Y}}_{\text{fr}1}\left(k\right)={\left[\begin{array}{ccc}{\widetilde{\underset{¯}{y}}}_{\text{fr}1}\left(k+1/k\right)& \cdots & {\widetilde{\underset{¯}{y}}}_{\text{fr}1}\left(k+N_p/k\right)\end{array}\right]}^T.\end{array}$$

(ii)

$\delta \tilde{Y}\left(k\right)$ is a $\left(m.N_p\right)$-dimensional vector defined as

$$\begin{array}{}\text{(68)}& \delta \tilde{Y}\left(k\right)={\left[\begin{array}{ccc}\delta \widetilde{\underset{¯}{y}}\left(k\right)& \cdots & \delta \widetilde{\underset{¯}{y}}\left(k+N_p-1\right)\end{array}\right]}^T.\end{array}$$

(iii)

$\delta U\left(k\right)$ is a $\left(p.N_u\right)$-dimensional future control increments vector where $N_u$ is the control horizon, defined as

$$\begin{array}{}\text{(69)}& \delta U\left(k\right)={\left[\delta \underset{¯}{u}\left(k\right),\delta \underset{¯}{u}\left(k+1\right),\dots ,\delta \underset{¯}{u}\left(k+N_u-1\right)\right]}^T.\end{array}$$

(iv)

$\delta \widehat{F}\left(k\right)$ is a $\left(p.N_u\right)$-dimensional actuator faults increments vector given by

$$\begin{array}{}\text{(70)}& \delta \widehat{F}\left(k\right)={\left[\delta \widehat{\underset{¯}{f}}\left(k\right),\delta \widehat{\underset{¯}{f}}\left(k+1\right),\dots ,\delta \widehat{\underset{¯}{f}}\left(k+N_u-1\right)\right]}^T.\end{array}$$

(v)

${\mathbf{G}}_a$ and ${\mathbf{G}}_b$ are, respectively, $\left(p.N_u\right)$-dimensional matrix and $\left(m.N_p\right)$-dimensional matrix given as follows:

$$\begin{array}{}\text{(71)}& {\mathbf{G}}_a=\left[\begin{array}{cccc}{\underset{¯}{g}}_{a,10}& {\underset{¯}{0}}_{m\times p}& \cdots & {\underset{¯}{0}}_{m\times p}\\ {\underset{¯}{g}}_{a,21}& {\underset{¯}{g}}_{a,20}& \cdots & {\underset{¯}{0}}_{m\times p}\\ ⋮& ⋮& \ddots & ⋮\\ {\underset{¯}{g}}_{a,N_u\left(N_u-1\right)}& {\underset{¯}{g}}_{a,N_u{N}_u-2}& \cdots & {\underset{¯}{g}}_{a,N_{u}0}\\ ⋮& ⋮& \ddots & ⋮\\ {\underset{¯}{g}}_{a,N_p\left(N_p\right)}& {\underset{¯}{g}}_{a,N_p\left(N_p-1\right)}& \cdots & {\underset{¯}{g}}_{a,N_p\left(N_p-N_u+1\right)}\end{array}\right],\end{array}$$where ${\underset{¯}{0}}_{m\times p}$ is a $\left(m\times p\right)$-dimensional vector of zeros.$$\begin{array}{}\text{(72)}& {\mathbf{G}}_b=\left[\begin{array}{cccc}{\underset{¯}{g}}_{b,10}& {\underset{¯}{0}}_{m\times p}& \cdots & {\underset{¯}{0}}_{m\times p}\\ {\underset{¯}{g}}_{b,21}& {\underset{¯}{g}}_{b,20}& \cdots & {\underset{¯}{0}}_{m\times p}\\ ⋮& ⋮& \ddots & ⋮\\ {\underset{¯}{g}}_{b,N_p{N}_p-1}& {\underset{¯}{g}}_{b,N_p{N}_p-2}& \cdots & {\underset{¯}{g}}_{b,N_{p}0}\end{array}\right],\end{array}$$where ${\underset{¯}{0}}_{m\times m}$ is a $\left(m\times m\right)$-dimensional vector of zeros.

According to relations (71) and (73), the components ${\underset{¯}{g}}_{a,jn}$ and ${\underset{¯}{g}}_{b,jn}$ $\left(n=0,\dots ,N_p-1\right)$ and $\left(j=1,\dots ,N_p\right)$ can be written as$$\begin{array}{}\text{(73)}& {\underset{¯}{g}}_{d,jn}=\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s{K}_{s,n}{\mathbf{B}}_{d,s},\hspace{1em}d=a,b.\end{array}$$

Then, from relation (66), the prediction vector $\tilde{Y}\left(k\right)$ given by (54) can be written in matrix form as$$\begin{array}{}\text{(74)}& \tilde{Y}\left(k\right)={\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_a\delta U\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right).\end{array}$$

6.2. Control Calculation

The control calculation of a MIMO system with p inputs and m outputs is based on the minimization of the following performance quadratic criterion:$$\begin{array}{}\text{(75)}& J_2\left(\delta U\left(k\right)\right)=\sum _{i=1}^m{\rho }_i\sum _{j=1}^{N_p}{\left[{\tilde{y}}_i\left(k+j/k\right)-w_i\left(k+j\right)\right]}^2+\sum _{t=1}^p{\varrho }_t\sum _{j=0}^{N_u-1}{\left(\delta u_t\left(k+j\right)\right)}^2,\end{array}$$where $w_i\left(k+j\right)$, $\left(i=1,\dots ,m\right)$ are the $i^{\text{th}}$ reference trajectories at time instant $\left(k+j\right)$, ${\rho }_i$ are the weights on the tracking error considered constant and equal to ρ, and ${\varrho }_t$ are positive weights generally considered constant and equal to $\varrho$.

The quadratic criterion $J_2\left(\delta U\left(k\right)\right)$ given by (75) can be rewritten on the following compact matrix:$$\begin{array}{}\text{(76)}& J_2\left(\delta U\left(k\right)\right)={\left[\tilde{Y}\left(k\right)-W\left(k\right)\right]}^T{\Lambda }_y\left[\tilde{Y}\left(k\right)-W\left(k\right)\right]+{\left(\delta U\left(k\right)\right)}^T{\Lambda }_u\delta U\left(k\right),\end{array}$$where $\tilde{Y}\left(k\right)$ is the $\left(m{N}_p\right)$-dimensional vector of predicted outputs defined in relation (65) and the vector $W\left(k\right)$ is defined from the reference trajectories $w_i$ as follows:$$\begin{array}{}\text{(77)}& W\left(k\right)={\left[w_1\left(k+1\right),\dots ,w_m\left(k+1\right),\dots ,w_1\left(k+N_p\right),\dots ,w_m\left(k+N_p\right)\right]}^T.\end{array}$$

$\delta U\left(k\right)$ is the vector of future control increments given by (69), ${\Lambda }_y$ and ${\Lambda }_u$ are positive definite matrices of dimension, respectively, $\left(m.N_p\right)$ and $\left(p.N_u\right)$, defined as follows:$$\begin{array}{}\text{(78)}& \left\{\begin{array}{c}{\Lambda }_y=\text{bloc}\_\text{diag}\underset{N_p\text{\hspace{0.17em}times}}{\underbrace{\left\{\underset{¯}{\rho },\dots ,\underset{¯}{\rho }\right\}}},\\ {\Lambda }_u=\text{bloc}\_\text{diag}\underset{N_u\text{times}}{\underbrace{\left\{\varrho ,\dots ,\varrho \right\}}},\end{array}\right.\end{array}$$where$$\begin{array}{c}\text{(79)}& \underset{¯}{\rho }=\left[\begin{array}{ccc}{\rho }_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\rho }_m\end{array}\right],\underset{¯}{\varrho }=\left[\begin{array}{ccc}{\varrho }_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\varrho }_p\end{array}\right].\end{array}$$

By replacing $\tilde{Y}\left(k\right)$ by its expression given by (74), the matrix form of criterion $J_2\left(\delta \underset{¯}{U}\left(k\right)\right)$ given by (76) becomes$$\begin{array}{}\text{(80)}& J_2\left(\delta U\left(k\right)\right)={\left[{\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_a\delta U\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right)-W\left(k\right)\right]}^T{\Lambda }_y\left[{\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_a\delta U\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right)-W\left(k\right)\right]+{\left(\delta U\left(k\right)\right)}^T{\Lambda }_u\delta U\left(k\right).\end{array}$$

The model predictive control methods allow to take into account some constraints optimization problem subject to physical constraints on the inputs and the outputs due to the actuator technology, the control system security, or the quality desired for the output of the controlled process. In the case of MIMO system of m outputs and p inputs, the constraints are given as$$\begin{array}{}\text{(81)}& \delta u_{t,m}\le \delta u_t\left(k+j\right)\le \delta u_{t,M},\hspace{1em}\forall j\in \left\{0,\dots ,N_u-1\right\},\hspace{1em}t=1,\dots ,p,\text{(82)}& u_{t,m}\le u_t\left(k+j\right)\le u_{t,M},\hspace{1em}\forall j\in \left\{0,\dots ,N_u-1\right\},\hspace{1em}t=1,\dots ,p,\text{(83)}& y_{i,m}\le y_i\left(k+j\right)\le y_{i,M},\hspace{1em}\forall j\in \left\{1,\dots ,N_p\right\},\hspace{1em}i=1,\dots ,m,\end{array}$$by analogy to the matrix form of constraints given by Mbarek et al. [22], these constraints can be expressed in terms of future control increments and written in matrix form as$$\begin{array}{}\text{(84)}& \Gamma \delta U\left(k\right)\le \underset{¯}{V}\left(k\right),\end{array}$$with$$\begin{array}{c}\text{(85)}& \Gamma =\left[\begin{array}{c}{I}_{\left(p{N}_u\right)}\\ -I_{\left(p{N}_u\right)}\\ D_1\\ -D_1\\ G_a\\ -G_a\end{array}\right],\underset{¯}{V}\left(k\right)=\left[\begin{array}{c}\delta {\underset{¯}{U}}_M\\ -\delta {\underset{¯}{U}}_m\\ {\underset{¯}{U}}_M-V_1\left(k-1\right)\\ -{\underset{¯}{U}}_m+V_1\left(k-1\right)\\ {\underset{¯}{Y}}_M-{\tilde{Y}}_{\text{fr}1}\left(k\right)-{\mathbf{G}}_a\delta \widehat{F}\left(k\right)-{\mathbf{G}}_b\delta \tilde{Y}\left(k\right)\\ -{\underset{¯}{Y}}_m+{\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right)\end{array}\right],\end{array}$$where the future control increment vector is defined by relation (69), and the square matrix $D_1$ and the vector ${\underset{¯}{V}}_1\left(k-1\right)$ of dimensions $p\times N_p$, are defined as follows:$$\begin{array}{}\text{(86)}& \begin{array}{c}\begin{array}{c}{D}_1=\left[\begin{array}{cccc}{I}_p& {\underset{¯}{0}}_{\left(p\times p\right)}& \cdots & {\underset{¯}{0}}_{\left(p\times p\right)}\\ I_p& I_p& \ddots & {\underset{¯}{0}}_{\left(p\times p\right)}\\ ⋮& ⋮& \ddots & ⋮\\ I_p& I_p& \cdots & I_p\end{array}\right],\\ {\underset{¯}{V}}_1\left(k-1\right)=\left[\begin{array}{c}\underset{¯}{u}\left(k-1\right)\\ \underset{¯}{u}\left(k-1\right)\\ ⋮\\ \underset{¯}{u}\left(k-1\right)\end{array}\right],\end{array}\end{array}\end{array}$$and the vectors ${\underset{¯}{U}}_M$, ${\underset{¯}{U}}_m$, $\delta {\underset{¯}{U}}_M$, and $\delta {\underset{¯}{U}}_m$ of dimension $\left(p.N_u\right)$ are defined as follows:$$\begin{array}{c}\text{(87)}& {\underset{¯}{U}}_M={\left[{\underset{¯}{u}}_M^T,\dots ,{\underset{¯}{u}}_M^T\right]}^T,{\underset{¯}{U}}_M={\left[{\underset{¯}{u}}_m^T,\dots ,{\underset{¯}{u}}_m^T\right]}^T,\\ \delta {\underset{¯}{U}}_M={\left[\delta {\underset{¯}{u}}_M^T,\dots ,\delta {\underset{¯}{u}}_M^T\right]}^T,\\ \delta {\underset{¯}{U}}_M={\left[\delta {\underset{¯}{u}}_m^T,\dots ,\delta {\underset{¯}{u}}_m^T\right]}^T,\end{array}$$and the vectors ${\underset{¯}{y}}_M$ and ${\underset{¯}{y}}_m$ of dimension $m.N_p$ are given by$$\begin{array}{}\text{(88)}& \begin{array}{c}\begin{array}{c}{\underset{¯}{y}}_M={\left[{\underset{¯}{y}}_M^T,\dots ,{\underset{¯}{y}}_M^T\right]}^T,\\ {\underset{¯}{y}}_m={\left[{\underset{¯}{y}}_m^T,\dots ,{\underset{¯}{y}}_m^T\right]}^T.\end{array}\end{array}\end{array}$$

These constraints define a future control increments admissible set as follows:$$\begin{array}{}\text{(89)}& \Psi =\left\{\delta U\left(k\right)/\Gamma \delta U\left(k\right)\le \underset{¯}{V}\left(k\right)\right\}.\end{array}$$

Then, the optimization problem subject to the physical constraints on the inputs and the outputs is written as follows:$$\begin{array}{}\text{(90)}& \begin{array}{c}\underset{\delta \underset{¯}{U}\left(k\right)\in \Psi }{\min }\left(J_2\left(\delta U\left(k\right)\right)\right)\end{array},\end{array}$$where $\Psi$ is the future control increments admissible set given by (89) and $J_2\left(\delta U\left(k\right)\right)$ is the quadratic criterion that can be written, from (90), as$$\begin{array}{}\text{(91)}& J_2\left(\delta U\left(k\right)\right)=\delta U^T\left(k\right)\left[{\mathbf{G}}_a^T{\Lambda }_y{\mathbf{G}}_a+{\Lambda }_u\right]\delta U\left(k\right)+2{\left({\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right)-\underset{¯}{W}\left(k\right)\right)}^T{\Lambda }_y{\mathbf{G}}_a\delta U\left(k\right)+{\left({\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right)-\underset{¯}{W}\left(k\right)\right)}^T{\Lambda }_y\left({\tilde{Y}}_{\text{fr}1}\left(k\right)+{\mathbf{G}}_a\delta \widehat{F}\left(k\right)+{\mathbf{G}}_b\delta \tilde{Y}\left(k\right)-\underset{¯}{W}\left(k\right)\right).\end{array}$$

The optimization problem (90) subject to (84) is a convex quadratic programming problem (QP), for which the global solution is unique. As with all predictive control methods, where a moving horizon strategy is applied, at instant $kT$, only the first optimized control $\underset{¯}{U}\left(k\right)$ is applied to the system, with$$\begin{array}{}\text{(92)}& \underset{¯}{u}\left(k\right)=\underset{¯}{u}\left(k-1\right)+\delta \underset{¯}{u}\left(k\right).\end{array}$$

7. Numerical Example

The proposed MIMO ARX-Laguerre multiple models, the fault detection algorithm, and the fault tolerant control via MPC are validated on numerical simulation using the following nonlinear MIMO discrete-time system characterized by two inputs $p=2$ and two outputs $\left(m=2\right)$ [14, 15]:$$\begin{array}{}\text{(93)}& \left\{\begin{array}{c}{y}_{s1}\left(k\right)=0.4y_{s1}\left(k-1\right)+\frac{u_1\left(k-1\right)}{1+u_1^2\left(k-1\right)}+0.2u_1^3\left(k-1\right)+0.5u_2\left(k-1\right),\\ y_{s2}\left(k\right)=0.2y_{s2}\left(k-1\right)+\frac{u_2\left(k-1\right)}{1+u_2^2\left(k-1\right)}+0.4u_2^3\left(k-1\right)+0.2u_1\left(k-1\right).\end{array}\right.\end{array}$$

The MIMO nonlinear system given by (93) was simulated, and a set of training input-output data were obtained. These inputs and outputs systems are represented, respectively, by Figures 3 and 4.

[figure omitted; refer to PDF]

[figures omitted; refer to PDF]

To describe the MIMO nonlinear system by the proposed MIMO ARX-Laguerre multiple models, it is necessary to determine the weighting functions ${\mu }_s\left(\gamma \left(k\right)\right)$ defined by relation (3) for $\left(L=3\right)$ submodels where we use the output systems as the decision variables. The weighting functions ${\mu }_s\left(\gamma \left(k\right)\right)$ for $\left(s=\mathrm{1,2,3}\right)$ are represented in Figure 5.

[figure omitted; refer to PDF]

From the input/output couples, we use 800 couples to identify the proposed nonlinear MIMO ARX-Laguerre multiple models where we applied the Algorithm 1 to identify the Fourier coefficients and the Algorithm 2 to optimize the poles. After several simulations for different truncation orders $N_a^s$ and $N_b^s$ for $\left(s=\mathrm{1,2,3}\right)$, the optimal choice of the latter is $N_a^s=N_b^s=N=2$. Then, in our case, the nonlinear MIMO ARX-Laguerre multiple models is obtained by optimizing the $\left(m.L.\left(m+p\right)=24\right)$ poles ${\underset{¯}{\xi }}_i$, $\left(i=1,\dots ,m\right)$, given by relation (29) and identifying the $\left(m.L.N.\left(m+p\right)=48\right)$ Fourier coefficients ${\underset{¯}{C}}_i$, $\left(i=1,\dots ,m\right)$, defined by relation (27). In the following, we give the optimized values of the poles and the Fourier coefficients (Algorithm 3):$$\begin{array}{c}\text{(94)}& C_1=\left[0.7643,-0.0075,0.6586,0.0119,0.3399,-0.0028,-0.3532,-0.0014,0.8595,-0.019,0.5649,-0.0813,-0.0198,-0.0039,-0.0406,0.0197,0.7640,-0.0268,0.5633,-0.1051,0.1968,-0.0255,-0.1096,0.0314\right],C_2=\left[0.2262,0.0067,0.8839,-0.0045,-0.1075,-0.0033,0.1704,0.0016,0.2018,-0.0021,0.9055,0.0004,0.0146,-0.0002,-0.0070,-0.0000,0.2082,0.0094,0.9045,0.0102,-0.0187,-0.0093,-0.0087,0.0032\right],\\ \text{(95)}& \begin{array}{c}{\underset{¯}{\xi }}_1=\left[0.0777,0.5907,0.2871,0.6100,0.4433,0.5557,0.4709,0.5924,0.2050,0.5529,0.4413,0.6274\right],\\ {\underset{¯}{\xi }}_2=\left[0.0777,0.5907,0.2871,0.6100,0.4433,0.5557,0.4709,0.5924,0.2050,0.5529,0.4413,0.6274\right].\end{array}\end{array}$$

Algorithm 1: Recursive parametric identification.

(1)

Assume we have M inputs/outputs couples $\left(u_1,\dots ,u_p\right)$ and $\left({\overline{y}}_1,\dots ,{\overline{y}}_m\right)$

(2)

Choose the number of submodels and fix the truncating orders $N_i^s=N$

(3)

Assume that the Laguerre poles vectors ${\underset{¯}{\xi }}_i$, $\left(i=1,\dots ,m\right)$ are optimized

(4)

Calculate the matrices $A_i^s$ and $B_i^s$ for $\left(i=1,\dots ,m\right)$ and $\left(s=1,\dots ,L\right)$ given by (23) and (25), respectively

(5)

Initialize the parameter vectors ${\underset{¯}{C}}_i={\underset{¯}{0}}_{\left(N.L\left(p+m\right)\right)\times 1}$, $y_i\left(1\right)=0$ and the state vectors $X_i^s\left(1\right)={\underset{¯}{0}}_{\left(N.\left(p+m\right)\right)\times 1}$ for $\left(i=1,\dots ,m\right)$ and $\left(s=1,\dots ,L\right)$

(6)

For $h=1,\dots ,M$

(a)

Calculate $X_i\left(h\right)$ for $\left(i=1,\dots ,m\right)$ from relations (28) and (A.17)

(b)

Deduce the matrix ${\mathbf{\Psi }}_i\left(h\right)$ and form the vector ${\overline{Y}}_i\left(h\right)$ as given in relation (39)

(c)

Compute the estimated Fourier coefficients at time instant h, ${\underset{¯}{C}}_i\left(h\right)$, using (41)

(d)

$h=h+1$, go to step (a)

Algorithm 2: Pole optimization algorithm.

(1)

Assume we have M measured input/output system $\left(u_1,\dots ,u_p\right)$ and $\left({\overline{y}}_1,\dots ,{\overline{y}}_m\right)$

(2)

Choose the size of population N

(3)

Fix a stopping criterion

(4)

Fix an initial set of N pole vectors ${\underset{¯}{\xi }}_i$, $\left(i=1,\dots ,m\right)$, given by relation (29)

(5)

For each poles vector, identify the parameters vector ${\underset{¯}{C}}_i$ for $\left(i=1,\dots ,m\right)$, using Algorithm 1

(6)

Compute the outputs $y_i$ for $\left(i=1,\dots ,m\right)$ from relation (21)

(7)

Evaluate the fitness values of the initial population ${\text{NMSE}}_i$ for each MISO model, $\left(i=1,\dots ,m\right)$

(8)

do $\left\{\right.$

(a)

Select parents

(b)

Crossover

(c)

Mutation

(d)

Next population

(e)

Identify parameter vectors and compute the outputs for each poles vectors

(f)

Evaluate the fitness values of the population ${\text{NMSE}}_i$ for every MISO model, $\left.\left(i=1,\dots ,m\right)\right\}$

(9)

Until stopping criteria reached

(10)

The last population is the optimal value of the MIMO ARX-Laguerre multiple models poles

Algorithm 3: MIMO nonlinear fault tolerant control via MPC (MIMO NFTC-MPC).

(1)

The nonlinear system is modeled by the MIMO ARX-Laguerre multiple models for L fixed linear submodels for each MISO system, where the parameters vector ${\underset{¯}{C}}_i$, $\left(i=1,\dots ,m\right)$, is identified by Algorithm 1 and the pole vector ${\underset{¯}{\xi }}_i$, $\left(i=1,\dots ,m\right)$, is optimized by Algorithm 2

(2)

Choose the prediction horizon $N_p$, the control horizon $N_u$, the weights on the tracking error ρ, and the positive weights ϱ

(3)

Determine the positive definite matrices ${\Lambda }_y$ and ${\Lambda }_u$ given by relation (78)

(4)

Choose the physical constraints given by (81)–(83)

(5)

Fix the fault actuators constraints $f_m$ and $f_M$

(6)

Compute the vectors ${\underset{¯}{U}}_M$, ${\underset{¯}{U}}_m$, $\Delta {\underset{¯}{U}}_M$, and $\Delta {\underset{¯}{U}}_m$ given by (87) and the vectors ${\underset{¯}{y}}_M$ and ${\underset{¯}{y}}_m$ given by (88)

(7)

Compute these matrices $\mathbf{A}$, ${\mathbf{B}}_a$, and ${\mathbf{B}}_b$ given by (36) and the matrix $\mathbf{C}$ given by (35)

(8)

Calculate the matrices ${\mathbf{G}}_a$, ${\mathbf{G}}_b$, and $D_1$ given, respectively, by (71), (72), and (86)

(9)

Compute the matrix $\Gamma$ given by (85)

(10)

Choose the reference signals $w_i\left(k\right)$ and initialize the input $u_i\left(1\right)$, for $\left(i=1,\dots ,m\right)$

(11)

for $k=2$ until $\left(M-N_p\right)$,

(a)

Apply the inputs $u_i\left(k-1\right)$ to the system and read the output $y_i\left(k\right)$, for $\left(i=1,\dots ,m\right)$

(b)

Optimize the fault actuators $\widehat{f}\left(k-1\right)$ using the MHE method from relation (49)

(c)

Compute $\delta {\tilde{X}}_s\left(k\right)$ and ${\tilde{X}}_s\left(k\right)$ given, respectively, by relations (54) and (55a)

(d)

For $\left(j=1,\dots ,N_p\right)$, compute ${\widetilde{\underset{¯}{y}}}_{\text{fr}1}\left(k+j/k\right)$ given by relation (72)

(e)

Compute these vectors ${\tilde{Y}}_{\text{fr}1}\left(k\right)$, $\delta \tilde{Y}\left(k\right)$, and $\delta \widehat{F}\left(k\right)$ given, respectively, by relations (77), (78), and (80) and compute the vector $W\left(k\right)$ given by (77)

(f)

Compute the vector $\underset{¯}{V}\left(k-1\right)$ given by (85)

(g)

Solve the optimization problem given by (90)

(h)

Calculate the control signal $\underset{¯}{u}\left(k\right)$ from relation (92)

(12)

END

Using the optimal poles ${\underset{¯}{\xi }}_i$, $\left(i=1,\dots ,m\right)$, given by (95), we can deduce the matrices $\mathbf{A}$, ${\mathbf{B}}_a$, and ${\mathbf{B}}_b$ given by relation (36), and using the identified parameters ${\underset{¯}{C}}_i$, $\left(i=1,\dots ,m\right)$, given by (94) and the weighting function ${\mu }_s\left(\gamma \left(k\right)\right)$ for $\left(s=\mathrm{1,2,3}\right)$, we can determine the vector $\mathbf{C}$ defined by relation (32). Then, the state space representation (31) of the proposed MIMO ARX-Laguerre multiple models representing the nonlinear MIMO system given by (93) is determined. In Figures 6 and 7, we draw the validation of the identified model.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

The identified MIMO ARX-Laguerre multiple models are used to detect the fault actuator by applying the MHE technique developed in Section 5. To validate the capability of the MHE using MIMO ARX-Laguerre multiple models, we applied to the nonlinear MIMO system given by (93), a constant input $u_1=0.5$ and $u_2=0.6$ and we added to the input the fault given in Figure 8. By considering constraint to the faults given by relation (96), the optimal faults are given in Figure 9.$$\begin{array}{}\text{(96)}& \left\{\begin{array}{c}0\le f_1\le 0.25,\\ 0\le f_2\le \mathrm{0.25.}\end{array}\right.\end{array}$$

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

To validate the proposed fault tolerant control algorithm for nonlinear MIMO systems via MPC based on ARX-Laguerre multiple models Algorithm 3, we used constant reference trajectories $w_1\left(k\right)=0.4$ and $w_2\left(k\right)=0.6$ for $\left(k=1,\dots ,200\right)$ and we applied to the numerical model given by relation (93), an input fault at each sample time as given in Figure 10. In Figures 11 and 12, we represent the trajectory references and the system outputs, respectively, for the first output and second output without FTC where we observe the influence of the input faults on the output systems. In Figures 13 and 14, we draw the system outputs and the reference trajectories in the case of FTC via MPC using nonlinear MIMO ARX-Laguerre multiple models, and in Figure 15, we plot the fault optimization.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

7.1. Example 2: Communicating Two-Tank System

7.1.1. Description of the Experimental System

The proposed faults tolerant control scheme presented in this paper is tested on a MIMO nonlinear MIMO communicating two-tank system (CTTS). The picture of the MIMO CTTS is given in Figure 16, and its block diagram is represented in Figure 17. The CTTS contains two voltage-controlled motor pumps ($P_1$ and $P_2$) which deliver the flow rates $Q_1$ and $Q_2$ inside of the tank 1 and the tank 2, respectively. Each tank is equipped with a level sensor and a valve for draining. This experimental process allows the control of both water levels $h_1$ and $h_2$ in tanks 1 and 2, respectively, (outputs) by tuning water flow rates $Q_1$ and $Q_2$ provided by the motor-pumps $P_1$ and $P_2$ (inputs). The notations used in Figure 17 are given in Table 1. The experimental process shown in Figure 16 operates as follows: the motor-pump $P_1$ fills the tank 1 with a flow rate $Q_1$ and the motor-pump $P_2$ fills the tank 2 with a flow rate $Q_2$. The manual valves $V_1$ and $V_2$ ensure the draining of tanks 1 and 2, respectively, with flow rates $Q_{10}$ and $Q_{20}$. Furthermore, if the levels $h_1$ and $h_2$ are different, a draining of tank 1 or 2 with the flow rate $Q_{12}$ is performed through the interconnection valve $V_C$.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

Table 1

Notations used for the CTTS system.

$P_1$	Motor-pump 1 with variable flow rate between 0 and $7\text{l}/\text{min}$ at $12\text{V}$ DC
$P_2$	Motor-pump 2 with variable flow rate between 0 and $3\text{l}/\text{min}$ at $12\text{V}$ DC
$V_1$	Manual drain valve of tank 1 of cross section $S_1$
$V_C$	Manual interconnection valve between the two tanks of cross section $S_{12}$
$V_2$	Manual drain valve of tank 2 of cross section $S_2$
$V_E$	Electromagnetic valve of cross section $S_E$
$Q_1$	Flow rate provided by motor-pump 1
$Q_2$	Flow rate provided by motor-pump 2
$Q_{12}$	Flow rate between tanks 1 and 2
$Q_{10}$	Flow rate drain of tank 1
$Q_{20}$	Flow rate drain of tank 2
$V_{3V}$	Three-way valve
$H_1$	Level sensor in the tank 1
$H_2$	Level sensor in the tank 2

The measurement platform and experimental process control represented in Figure 18 are characterized by the interface card plugged in the computer slot generating two PWM (pulse width modulation), the DAS-1600 card for the acquisition of water levels $h_1$ and $h_2$ and the control card of the variation of the flow rates $Q_1$ and $Q_2$. This card transfers the PWM signals generated by the interface card to the motor-pumps $P_2$ and $P_2$.

[figure omitted; refer to PDF]

The nonlinear CCTS system can be described by a nonlinear MIMO model. In this paper, we will model this system by the proposed nonlinear MIMO ARX-Laguerre multiple models. In the following, the identified model is used to validate the MIMO MHE method proposed in Section 5 by the estimation of the actuator faults. Finally, we validate the proposed Algorithm 2 “MIMO nonlinear fault tolerant control via MPC” developed based on “nonlinear MIMO ARX-Laguerre multiple models.”

7.2. Identification of Nonlinear MIMO ARX-Laguerre Multiple Models

In this subsection, we validate the proposed nonlinear MIMO ARX-Laguerre multiple models to model the CTTS process. In this simulation, we fix the number of submodel to $L=3$. The process inputs $Q_1$ and $Q_2$ were chosen to be pseudo-random sequence with variable amplitude between $0\text{\%}$ and $100\text{\%}$ as represented in Figure 19.

[figures omitted; refer to PDF]

To describe the MIMO nonlinear CTTS by the proposed MIMO ARX-Laguerre multiple models, we use the output systems as decision variables $\gamma \left(k\right)$. From the input/output couples, we use 2000 couples to identify the proposed nonlinear MIMO ARX-Laguerre multiple models where we applied the Algorithm 1 to identify the Fourier coefficients and the Algorithm 2 to optimize the poles for truncation orders $N_a^s=N_b^s=N=2$. And, we use last 1000 couples to validate the identified model as represented in Figures 20 and 21.

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

7.3. Fault Tolerant Control for CTTS via MPC

To control the CTTS using the proposed fault tolerant control algorithm via MPC based on MIMO ARX-Laguerre multiple models Algorithm 3, we used the reference trajectories $w_1$ and $w_2$ and we applied to the CTTS process an input fault actuators $f_1$ and $f_2$ at each sample time as given in Figure 22. The weighting functions ${\mu }_s\left(\gamma \left(k\right)\right)$, for $\left(s=\mathrm{1,2,3}\right)$, are determined in terms of the reference trajectories. In Figures 23 and 24, we draw the trajectory references and the system outputs, respectively, for the first outputs and second outputs without FTC where we observe the influence of the input faults on the output systems. By considering the constraint to the actuator faults given by relation (97), the optimal actuator faults $f{1}_{\text{opt}}$ and $f{2}_{\text{opt}}$ are given in Figure 22. In Figures 25 and 26, we draw the system outputs and the reference trajectories in the case of FTC via MPC using nonlinear MIMO ARX-Laguerre multiple models. We compare Figures 23–26, and observe the elimination of the effect of the actuator faults on the outputs by using the proposed fault tolerant control algorithm:$$\begin{array}{}\text{(97)}& \left\{\begin{array}{c}-10\le f_1\le 0,\\ -5\le f_2\le 0.\end{array}\right.\end{array}$$

[figures omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

[figure omitted; refer to PDF]

8. Conclusion

In this paper, we have proposed a new MIMO nonlinear fault tolerant via model predictive control to control MIMO nonlinear systems by taking into account fault actuators. The MIMO nonlinear systems are described by a MIMO ARX-Laguerre multiple models. One more important advantage of the proposed MIMO ARX-Laguerre multiple models is that its structure and its reduced parameter number can facilitate the controller design problem and reduce the computational effort of the MIMO NFTC-MPC. The fault actuators estimation is assured by the MHE method with respect to these constraints on these faults. The proposed MIMO NFTC-MPC is synthesized by taking into account the estimation fault actuators. The optimization problem of the MIMO NFTC-MPC is a minimization of a quadratic criterion with respect to constraints on input signals. The control algorithm is easily solved and yields efficient results. The stability of closed-loop system is guaranteed with respect to the choice of the tuning parameters.

Appendix

Appendix

A. Formulation of Recursive Representation

Defining by $X_{it,n}^s\left({\xi }_{it}^s,k\right)$ and $X_{it,n}^s\left({\xi }_{it}^s,k\right)$, respectively, the filtered inputs and the filtered outputs by the $s^{\text{th}}$ orthonormal functions as$$\begin{array}{}\text{(A.1)}& \left\{\begin{array}{c}{x}_{a_{it},n}^s\left({\xi }_{a_{it}}^s,k\right)=\sum _{j=1}^{\infty }{\ell }_{a_{it},n}^s\left({\xi }_{a_{it}}^s,j\right)u_t\left(k-j\right)={\ell }_{a_{it},n}^s\left({\xi }_{a_{it}}^s,k\right)\ast u_t\left(k\right),\\ x_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,k\right)=\sum _{j=1}^{\infty }{\ell }_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,j\right)y_r\left(k-j\right)={\ell }_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,k\right)\ast y_r\left(k\right),\end{array}\right.\end{array}$$where $\left(\ast \right)$ denotes the convolutions product. The resulting model entitled MISO ARX-Laguerre sub-model given by (20) can be rewritten as$$\begin{array}{}\text{(A.2)}& y_i^s\left(k\right)=\sum _{t=1}^p\sum _{n=0}^{\infty }{g}_{a_{it},n}^s{x}_{a_{it},n}^s\left({\xi }_{a_{it}}^s,k\right)+\sum _{r=1}^m\sum _{n=0}^{\infty }{g}_{b_{ir},n}^s{x}_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,k\right),\end{array}$$for $\left(s=1,\dots ,L\right)\text{\hspace{0.17em}and\hspace{0.17em}}\left(i=1,\dots ,m\right)$. According to the BIBO stability condition, the infinite series of the proposed MISO ARX-Laguerre submodel given by (A.2) can be truncated, respectively, to finite orders $N_a^s$ and $N_b^s$ as$$\begin{array}{}\text{(A.3)}& y_i^s\left(k\right)=\sum _{t=1}^p\sum _{n=0}^{N_a^s-1}{g}_{a_{it},n}^s{x}_{a_{it},n}^s\left({\xi }_{a_{it}}^s,k\right)+\sum _{r=1}^m\sum _{n=0}^{N_b^s-1}{g}_{b_{ir},n}^s{x}_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,k\right).\end{array}$$

Then, we obtained a reduced complexity MIMO ARX-Laguerre multiple models where each MISO ARX-Laguerre multiple model is defined from relations (4) and (A.3) as$$\begin{array}{}\text{(A.4)}& y_i\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right)\left(\sum _{t=1}^p\sum _{n=0}^{N_a^s-1}{g}_{a_{it},n}^s{x}_{a_{it},n}^s\left({\xi }_{a_{it}}^s,k\right)+\sum _{r=1}^m\sum _{n=0}^{N_b^s-1}{g}_{b_{ir},n}^s{x}_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,k\right)\right).\end{array}$$

The Z-transform of relation (A.4) can be written as$$\begin{array}{}\text{(A.5)}& y_i\left(k\right)=\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k-1\right)\right)\left(\sum _{t=1}^p\sum _{n=0}^{N_a^s-1}{g}_{a_{it},n}^s{X}_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)+\sum _{r=1}^m\sum _{n=0}^{N_b^s-1}{g}_{b_{ir},n}^s{X}_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)\right),\end{array}$$where $X_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)$ and $X_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)$ are, respectively, the Z-transform of $X_{a_{it},n}^s\left({\xi }_{a_{it}}^s,k\right)$ and $X_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,k\right)$ written from (A.1) as$$\begin{array}{}\text{(A.6)}& \left\{\begin{array}{c}{X}_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)=L_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)\cdot U_t\left(z\right),\\ X_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)=L_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)\cdot Y_r\left(z\right),\end{array}\right.\end{array}$$where $L_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)$ and $L_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)$ are the Z-transform of the Laguerre filters given by their recursive form as follows:$$\begin{array}{}\text{(A.7)}& \left\{\begin{array}{c}{L}_{a_{it},0}^s\left({\xi }_{a_{it}}^s,z\right)=\frac{\sqrt{1-{\left({\xi }_{a_{it}}^s\right)}^2}}{z-{\xi }_{a_{it}}^s},\\ L_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)=\frac{1-{\xi }_{a_{it}}^{s}z}{z-{\xi }_{a_{it}}^s}{L}_{a_{it},n-1}^s\left({\xi }_{a_{it}}^s,z\right),\hspace{1em}n=\mathrm{1,2},\dots ,\end{array}\right.\text{(A.8)}& \left\{\begin{array}{c}{L}_{b_{ir},0}^s\left({\xi }_{b_{ir}}^s,z\right)=\frac{\sqrt{1-{\left({\xi }_{b_{ir}}^s\right)}^2}}{z-{\xi }_{b_{ir}}^s},\\ L_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)=\frac{1-{\xi }_{b_{ir}}^{s}z}{z-{\xi }_{b_{ir}}^s}{L}_{b_{ir},n-1}^s\left({\xi }_{b_{ir}}^s,z\right),\hspace{1em}n=\mathrm{1,2},\dots \end{array}\right.\end{array}$$

From the recursive form of the Laguerre filters given by relation (A.7) and (A.8), we can deduce a recursive form of $X_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)$ and $X_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)$ as follows:$$\begin{array}{}\text{(A.9)}& \left\{\begin{array}{c}{X}_{a_{it},0}^s\left({\xi }_{a_{it}}^s,z\right)=\frac{\sqrt{1-{\left({\xi }_{a_{it}}^s\right)}^2}}{z-{\xi }_{a_{it}}^s}{U}_t\left(z\right),\\ X_{a_{it},n}^s\left({\xi }_{a_{it}}^s,z\right)=\frac{1-{\xi }_{a_{it}}^{s}z}{z-{\xi }_{a_{it}}^s}{X}_{a_{it},n-1}^s\left({\xi }_{a_{it}}^s,z\right),\hspace{1em}n=\mathrm{1,2},\dots ,\end{array}\right.\text{(A.10)}& \left\{\begin{array}{c}{X}_{b_{ir},0}^s\left({\xi }_{b_{ir}}^s,z\right)=\frac{\sqrt{1-{\left({\xi }_{b_{ir}}^s\right)}^2}}{z-{\xi }_{b_{ir}}^s}{Y}_r\left(z\right),\\ X_{b_{ir},n}^s\left({\xi }_{b_{ir}}^s,z\right)=\frac{1-{\xi }_{b_{ir}}^{s}z}{z-{\xi }_{b_{ir}}^s}{X}_{b_{ir},n-1}^s\left({\xi }_{b_{ir}}^s,z\right),\hspace{1em}n=\mathrm{1,2},\dots \end{array}\right.\end{array}$$

Then, by applying the Z-transform inverse to relations (A.9) and (A.10) and by defining the following vectors,$$\begin{array}{}\text{(A.11)}& \begin{array}{c}{X}_{a_{it}}^s\left(k\right)={\left[X_{a_{it},0}^s\left(k\right),X_{a_{it},1}^s\left(k\right),\dots ,X_{a_{it},N_a^s-1}^s\left(k\right)\right]}^T,\\ X_{b_{ir}}^s\left(k\right)={\left[X_{b_{ir},0}^s\left(k\right),X_{b_{ir},1}^s\left(k\right),\dots ,X_{b_{ir},N_b^s-1}^s\left(k\right)\right]}^T,\end{array}\end{array}$$we can compute the filters outputs $X_{a_{it},n}^s\left(k\right)$ and $X_{b_{ir},n}^s\left(k\right)$ recursively as$$\begin{array}{}\text{(A.12)}& \begin{array}{c}{X}_{a_{it}}^s\left(k+1\right)=A_{a_{it}}^s{X}_{a_{it}}^s\left(k\right)+b_{a_{it}}^s{u}_t\left(k\right),\\ X_{b_{ir}}^s\left(k+1\right)=A_{b_{ir}}^s{X}_{b_{ir}}^s\left(k\right)+b_{b_{ir}}^s{y}_r\left(k\right),\end{array}\end{array}$$where the matrices $A_{a_{it}}^s$ and $A_{b_{ir}}^s$ and the vectors $B_{a_{it}}^s$ and $B_{b_{ir}}^s$ are defined according to the poles ${\xi }_{a_{it}}^s$ and ${\xi }_{b_{ir}}^s$ as$$\begin{array}{}\text{(A.13)}& A_{a_{it}}^s=\left[\begin{array}{cccccc}{\xi }_{a_{it}}^s& 0& 0& \cdots & & 0\\ 1-{\left({\xi }_{a_{it}}^s\right)}^2& {\xi }_{a_{it}}^s& 0& \cdots & & 0\\ -{\xi }_{a_{it}}^s\left(1-{\left({\xi }_{a_{it}}^s\right)}^2\right)& 1-{\left({\xi }_{a_{it}}^s\right)}^2& {\xi }_{a_{it}}^s& \cdots & & 0\\ ⋮& ⋮& \ddots & ⋮& \\ {\left(-{\xi }_{a_{it}}^s\right)}^{\left(N_a^s-2\right)}\left(1-{\left({\xi }_{a_{it}}^s\right)}^2\right)& {\left(-{\xi }_{a_{it}}^s\right)}^{\left(N_a^s-3\right)}\left(1-{\left({\xi }_{a_{it}}^s\right)}^2\right)& {\left(-{\xi }_{a_{it}}^s\right)}^{\left(N_a^s-4\right)}\left(1-{\left({\xi }_{a_{it}}^s\right)}^2\right)& \cdots & & {\xi }_{a_{it}}^s\end{array}\right],\text{(A.14)}& A_{b_{ir}}^s=\left[\begin{array}{cccccc}{\xi }_{b_{ir}}^s& 0& 0& \cdots & & 0\\ 1-{\left({\xi }_{b_{ir}}^s\right)}^2& {\xi }_{b_{ir}}^s& 0& \cdots & & 0\\ -{\xi }_{b_{ir}}^s\left(1-{\left({\xi }_{b_{ir}}^s\right)}^2\right)& 1-{\left({\xi }_{b_{ir}}^s\right)}^2& {\xi }_{b_{ir}}^s& \cdots & & 0\\ ⋮& ⋮& \ddots & ⋮& \\ {\left(-{\xi }_{b_{ir}}^s\right)}^{\left(N_b^s-2\right)}\left(1-{\left({\xi }_{b_{ir}}^s\right)}^2\right)& {\left(-{\xi }_{b_{ir}}^s\right)}^{\left(N_b^s-3\right)}\left(1-{\left({\xi }_{b_{ir}}^s\right)}^2\right)& {\left(-{\xi }_{b_{ir}}^s\right)}^{\left(N_b^s-4\right)}\left(1-{\left({\xi }_{b_{ir}}^s\right)}^2\right)& \cdots & & {\xi }_{b_{ir}}^s\end{array}\right],\text{(A.15)}& b_{a_{it}}^s=\sqrt{1-{\left({\xi }_{a_{it}}^s\right)}^2}{\left[\begin{array}{ccccc}1& -{\xi }_{a_{it}}^s& {\left(-{\xi }_{a_{it}}^s\right)}^2& \cdots & {\left(-{\xi }_{a_{it}}^s\right)}^{\left(N_a^s-1\right)}\end{array}\right]}^T,\text{(A.16)}& b_{b_{ir}}^s=\sqrt{1-{\left({\xi }_{b_{ir}}^s\right)}^2}{\left[\begin{array}{ccccc}1& -{\xi }_{b_{ir}}^s& {\left(-{\xi }_{b_{ir}}^s\right)}^2& \cdots & {\left(-{\xi }_{b_{ir}}^s\right)}^{\left(N_b^s-1\right)}\end{array}\right]}^T.\end{array}$$

Defining the state vectors $X_i^s\left(k\right)$ as$$\begin{array}{}\text{(A.17)}& X_i^s\left(k\right)={\left[{\left(X_{a_{i1}}^s\left(k\right)\right)}^T,\dots ,{\left(X_{a_{ip}}^s\left(k\right)\right)}^T,{\left(X_{b_{i1}}^s\left(k\right)\right)}^T,\dots ,{\left(X_{b_{im}}^s\left(k\right)\right)}^T\right]}^T,\end{array}$$and the input vector $\underset{¯}{U}\left(k\right)$ and the output vector $\underset{¯}{y}\left(k\right)$ as$$\begin{array}{c}\text{(A.18)}& \underset{¯}{u}\left(k\right)={\left[u_1\left(k\right),u_2\left(k\right),\dots ,u_p\left(k\right)\right]}^T,\underset{¯}{y}\left(k\right)={\left[y_1\left(k\right),y_2\left(k\right),\dots ,y_m\left(k\right)\right]}^T.\end{array}$$

B. Proof of Proposition 1

By successive substitutions taking into account $\widetilde{\underset{¯}{y}}\left(k/k\right)=\widetilde{\underset{¯}{y}}\left(k\right)$, relation (55b) leads to write $\widetilde{\underset{¯}{y}}\left(k+j/k\right)$ as$$\begin{array}{}\text{(A.19)}& \begin{array}{c}\widetilde{\underset{¯}{y}}\left(k+j/k\right)=\widetilde{\underset{¯}{y}}\left(k\right)+\sum _{l=1}^j\left[\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+l-1\right)\right){\mathbf{C}}_s\delta {\tilde{X}}_s\left(k+l\right)+\sum _{s=1}^L\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k+l\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k+l\right)\right],\end{array}\end{array}$$also, ${\tilde{X}}_s\left(k+l\right)$ can be written in terms of $s^{\text{th}}$ state vector increments $\delta {\tilde{X}}_s\left(k\right)$:$$\begin{array}{}\text{(A.20)}& {\tilde{X}}_s\left(k+l\right)={\tilde{X}}_s\left(k\right)+\sum _{m=1}^L\delta {\tilde{X}}_s\left(k+m\right).\end{array}$$

By replacing in relation (A.19) ${\tilde{X}}_s\left(k+l\right)$ by its expression given by (A.20) and after some arrangement we obtain$$\begin{array}{}\text{(A.21)}& \widetilde{\underset{¯}{y}}\left(k+j/k\right)=\widetilde{\underset{¯}{y}}\left(k\right)+\sum _{l=1}^j\sum _{s=1}^L\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k+l\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\sum _{l=1}^j\delta {\tilde{X}}_s\left(k+l\right),\end{array}$$where, from (55a) and by successive substitutions, the state vector predictor $\delta {\tilde{X}}_s\left(k+l\right)$ is$$\begin{array}{}\text{(A.22)}& \delta {\tilde{X}}_s\left(k+l\right)={\left({\mathbf{A}}_s\right)}^l\delta {\tilde{X}}_s\left(k\right)+\sum _{m=1}^L{\left({\mathbf{A}}_s\right)}^{l-m}\left[{\mathbf{B}}_{a,s}\delta \underset{¯}{u}\left(k+m-1\right)+{\mathbf{B}}_{a,s}\delta \widehat{\underset{¯}{f}}\left(k+m-1\right)+{\mathbf{B}}_{b,s}\delta \underset{¯}{\tilde{y}}\left(k+m-1\right)\right].\end{array}$$

By replacing in equation (A.21) $\delta {\tilde{X}}_s\left(k+l\right)$ by its expression given by (A.22) and after some arrangement we obtain$$\begin{array}{}\text{(A.23)}& \begin{array}{c}\widetilde{\underset{¯}{y}}\left(k+j/k\right)=\widetilde{\underset{¯}{y}}\left(k\right)+\sum _{s=1}^L\sum _{l=1}^j\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k+l\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\left(\sum _{l=1}^j{\left({\mathbf{A}}_s\right)}^l\right)\delta {\tilde{X}}_s\left(k\right)\\ +\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\sum _{l=1}^j\sum _{h=0}^{j-l}{\left({\mathbf{A}}_s\right)}^h\left[{\mathbf{B}}_{a,s}\delta \underset{¯}{u}\left(k+l-1\right)+{\mathbf{B}}_{a,s}\delta \widehat{\underset{¯}{f}}\left(k+l-1\right)+{\mathbf{B}}_{b,s}\delta \underset{¯}{\tilde{y}}\left(k+l-1\right)\right],\end{array}\end{array}$$defining$$\begin{array}{c}\text{(A.24)}& K_{s,n}=\sum _{r=0}^n{\left({\mathbf{A}}_s\right)}^r,\hspace{1em}\text{for\hspace{0.17em}}n\ge 0,K_{s,n}=0,\hspace{1em}\text{for\hspace{0.17em}}n<0,\end{array}$$equation (A.23) becomes$$\begin{array}{}\text{(A.25)}& \begin{array}{c}\widetilde{\underset{¯}{y}}\left(\frac{k+j}{k}\right)=\widetilde{\underset{¯}{y}}\left(k\right)+\sum _{s=1}^L\sum _{l=1}^j\delta {\mu }_s\left(\underset{¯}{\gamma }\left(k+l\right)\right){\mathbf{C}}_s{\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\left[K_{\left(s,j\right)}-I_{N_s}\right]\delta {\tilde{X}}_s\left(k\right)+\sum _{s=1}^L{\mu }_s\left(\underset{¯}{\gamma }\left(k+j\right)\right){\mathbf{C}}_s\sum _{l=1}^j{K}_{\left(s,j-l\right)}\left({\mathbf{B}}_{a,s}\delta \underset{¯}{u}\left(k+l-1\right)+{\mathbf{B}}_{a,s}\delta \widehat{\underset{¯}{f}}\left(k+l-1\right)+{\mathbf{B}}_{b,s}\delta \underset{¯}{\tilde{y}}\left(k+l-1\right)\right),\\ \end{array}\end{array}$$where $I_{N_s}$ is the $N_s$-dimensional identity matrix and $N_s$ is the truncation order of the $s^{\text{th}}$ of the L submodel for each MISO submodel.

References

[1] A. H. Mohamed, A. M. Nassar, "New algorithm for fault diagnosis of photovoltaic energy systems," International Journal of Computer Applications, vol. 114 no. 9, pp. 26-31, DOI: 10.5120/20008-1959, 2015.

[2] M. Davoudi, E. Kamyab, J. Sadeh, "Parameter-free fault location for transmission lines based on optimisation," IET Generation, Transmission & Distribution, vol. 9 no. 11, pp. 1061-1068, DOI: 10.1049/iet-gtd.2014.0425, 2015.

[3] N. L. Georgijevic, M. V. Jankovic, S. Srdic, Z. Radakovic, "The detection of series arc fault in photovoltaic systems based on the arc current entropy," IEEE Transactions on Power Electronics, vol. 31 no. 8, pp. 5917-5930, DOI: 10.1109/tpel.2015.2489759, 2016.

[4] A. Mellit, G. M. Tina, S. A. Kalogirou, "Fault detection and diagnosis methods for photovoltaic systems: a review," Renewable and Sustainable Energy Reviews, vol. 91,DOI: 10.1016/j.rser.2018.03.062, 2018.

[5] M. M. Morato, O. Sename, L. Dugard, "LPV-MPC fault tolerant control of automotive suspension dampers," IFAC-PapersOnLine, vol. 51 no. 26, pp. 31-36, DOI: 10.1016/j.ifacol.2018.11.172, 2018.

[6] X. Peng, S. Jia, X. Wang, "Analytic model predictive control of uncertain nonlinear systems: a fuzzy adaptive approach," Mathematical Problems in Engineering, vol. 2015,DOI: 10.1155/2015/760490, 2015.

[7] S. Gao, L. Luan, "The PVC stripping process predictive control based on the implicit algorithm," Mathematical Problems in Engineering, vol. 2014,DOI: 10.1155/2014/838404, 2014.

[8] K. Bouzrara, A. Mbarek, T. Garna, "Non-linear predictive controller for uncertain process modelled by GOBF-Volterra models," International Journal of Modelling, Identification and Control, vol. 19 no. 4, pp. 307-322, DOI: 10.1504/ijmic.2013.055649, 2013.

[9] D. Li, Y. Xi, "The synthesis of robust model predictive control with QP formulation," International Journal of Modelling, Identification and Control, vol. 13 no. 1/2,DOI: 10.1504/ijmic.2011.040484, 2011.

[10] H.-t. Zhen, X.-h. Qi, J. Li, Q.-m. Tian, "Neural network l 1 adaptive control of mimo systems with nonlinear uncertainty," The Scientific World Journal, vol. 2014,DOI: 10.1155/2014/942094, 2014.

[11] H. P. H. Anh, N. T. Nam, "Novel adaptive forward neural mimo narx model for the identification of industrial 3-dof robot arm kinematics," International Journal of Advanced Robotic Systems, vol. 9 no. 4,DOI: 10.5772/51277, 2012.

[12] Y.-X. Li, G.-H. Yang, S. Tong, "Event-triggered adaptive fuzzy tracking control of nonlinear MIMO systems," International Journal of Systems Science, vol. 49 no. 12, pp. 2618-2629, DOI: 10.1080/00207721.2018.1506523, 2018.

[13] A. Ajavon, K. M. Pekpe, J. P. Cassar, "On efficiency of identification methods for mimo non linear system using Takagi Sugeno models," IFAC Proceedings Volumes, vol. 44 no. 1, pp. 15000-15006, DOI: 10.3182/20110828-6-it-1002.03636, 2011.

[14] I. Benabdelwahed, A. Mbarek, K. Bouzrara, T. Garna, "Non-linear system modelling based on narx model expansion on Laguerre orthonormal bases," IET Signal Processing, vol. 12 no. 2, pp. 228-241, DOI: 10.1049/iet-spr.2017.0187, 2018.

[15] A. Sameh, M. Abdelkader, G. Tarek, R. José, "Optimal multimodel representation by Laguerre filters applied to a communicating two tank system," Journal of Systems Science and Complexity, vol. 31 no. 3, pp. 621-646, DOI: 10.1007/s11424-017-6047-2, 2018.

[16] G. Sarah, T. Garna, K. Bouzrara, J. Ragot, H. Messaoud, "Online identification of the bilinear model expansion on Laguerre orthonormal bases," International Journal of Control, vol. 87 no. 3, pp. 441-463, DOI: 10.1080/00207179.2013.835868, 2014.

[17] T. Garna, K. Bouzrara, J. Ragot, H. Messaoud, "Optimal expansions of discrete-time bilinear models using laguerre functions," IMA Journal of Mathematical Control and Information, vol. 31 no. 3, pp. 313-343, DOI: 10.1093/imamci/dnt013, 2013.

[18] G. Marouani, A. Mbarek, T. Garna, H. Messaoud, "Nonlinear model based predictive control using multiple models approach expanded on Laguerre bases," WSEAS Transactions on Systems and Control, vol. 10, pp. 113-126, 2015.

[19] M. Lauricella, M. Farina, R. Schneider, R. Scattolini, "A distributed fault detection and isolation algorithm based on moving horizon estimation," IFAC-PapersOnLine, vol. 50 no. 1, pp. 15259-15264, DOI: 10.1016/j.ifacol.2017.08.2406, 2017.

[20] K. Bouzrara, T. Garna, J. Ragot, H. Messaoud, "Decomposition of an arx model on laguerre orthonormal bases," ISA Transactions, vol. 51 no. 6, pp. 848-860, DOI: 10.1016/j.isatra.2012.06.005, 2012.

[21] I. B. Abdelwahed, A. Mbarek, K. Bouzrara, "Adaptive MPC based on MIMO ARX-Laguerre model," ISA Transactions, vol. 67, pp. 330-347, DOI: 10.1016/j.isatra.2016.11.017, 2017.

[22] A. Mbarek, K. Bouzrara, T. Garna, J. Ragot, H. Messaoud, "Laguerre-based modelling and predictive control of multi-input multi-output systems applied to a communicating two-tank system (CTTS)," Transactions of the Institute of Measurement and Control, vol. 39 no. 5, pp. 611-624, DOI: 10.1177/0142331215613150, 2017.