A Type-2 Fuzzy State Observer Model for Non-Stationary Dynamic System Identification: An Incremental Learning Method with Noise Handling

Abstract

Real-world identification involves dealing with challenges such as system complexity, noise, and uncertainties. In this context, a method for incremental learning is suggested, utilizing an evolving type-2 state observer fuzzy model. The process involves structure learning through an evolving type-2 multiscaling clustering approach, eliminating the need for data normalization. The estimation of linear state observer models for each rule is achieved using observer Markov parameters computed via a Type-2 Instrumental Variable (T2-IV) algorithm. For obtaining the instruments for the T2-IV algorithm, a recursive moving-average filter is used. Benchmark and online identification tasks are conducted to demonstrate the practicality and robustness of the proposed methodology, with performance comparisons against existing methodologies.

1. Introduction

Machine learning-based identification system is a relevant approach for modeling nonlinear, uncertain, multivariable, and complex systems. This approach aims to estimate a model that represents accurately the dynamic behavior of a real plant [1]. In these terms, fuzzy identification arises as a relevant method for obtaining models which represent nonlinear systems. These techniques are shown a powerful tool in practical problems such as uncertainty, unpredictable dynamics, and noisy measurements. One reason is that fuzzy logic systems (FLS) have the ability to integrate information from different sources, such as physical laws, empirical models, or measurements [2]. For identification of problems with white noise signals, the classical fuzzy sets (type-1) present satisfactory results. However, when colored noise is considered, type-1 fuzzy sets are not able to mitigate the noise effects. In this case, type-2 fuzzy systems were proposed. The first mention of type-2 fuzzy sets as an extention of classical fuzzy sets was made by Zadeh [3], where theoretical concepts were addressed in the 1990s by Karnik [4]. In practical application, the interval type-2 fuzzy sets gained prominence in problems such as control [5] and modeling [6], as it presents less computational load. However, the development of methodologies for the experimental data analysis in order to obtain the rule-base for an interval type-2 fuzzy model in order to use the advantages of type-2 fuzzy sets described in the literature is still an open research field. In the literature, approaches such as heuristic methods [7] and incremental learning [8] have been used for this purpose. In Aissa Bencherif and Fatima Chouireb’s work [9], an incremental learning algorithm for type-2 recurrent Takagi-Sugeno neural-fuzzy network is proposed. For structure learning, the rule firing strength-based approach is used. For a new data, the type-2 firing strength is computed for each rule, where the rule with the highest firing strength is considered for creation rule mechanism. The parameter update is perfomed by gradient descent algorithm. The mobile robot trajectory tracking problem is used to show the applicability of the methodology. In [10], Morteza Montazeri-Gh and Shabnam Yazdani introduce the use of interval type-2 fuzzy logic systems for gas turbine fault diagnosis, aiming to reduce maintenance costs and downtime. Fuzzy Rule Base is estimated using Interval Type-2 Fuzzy C-Means clustering, and parameters of the IT2FLSs are optimized with a metaheuristic algorithm. The performance of the IT2FL-based FDI system is compared to other classification techniques, showing promising results in terms of online applicability, accuracy, and robustness.

In literature, linear models are commonly used in consequent part in Takagi-Sugeno models, such as vector auto-regressive and state-space models. The state-space models present an interesting feature: a compact formula that that shows the relationship between internal variables and the experimental data (output and input signals) [11]. In this context, methodologies based on fuzzy state-space models have been proposed [12, 13]. In Gil et al.’s work [14], a recurrent state-space neural-fuzzy network is introduced. For parameter adjustment of the antecedent/consequent parts, a recursive learning method based on the constrained unscented Kalman filter is employed. The applicability of this methodology is demonstrated through the online identification of a three-tank system. In Yancho et al.’s work [15], a fuzzy state-space model predictive control approach is proposed. The learning algorithm is founded on gradient descent, which is used to fit the modeling structure parameters. The trained model is subsequently applied in model-based predictive control. The applicability of this methodology is demonstrated through computational experiments.

In identification problems, efficiency in mitigating the effects of noise to adjust the consequent parameters must be ensured. In a noisy environment, the Instrumental Variable (IV) method is considered a relevant tool for system modeling [11]. When compared to other identification methods, it is noted that, in the IV method, the requirement for an accurate noise model is not essential [16]. According to literature, the accuracy of the IV method relies on the selection of an appropriate instrument, which must guarantee non-polarized estimation [11]. The fuzzy version of the IV method was introduced by Yancho. Therefore, with the aim of integrating the type-2 state-space fuzzy modeling and non-polarized consequent estimation, in this paper, an incremental learning for evolving interval type-2 state observer fuzzy model based on instrumental variables approach is proposed. For estimating the consequent-part, fuzzy observer Markov parameters are computed via type-2 fuzzy version of instrumental variable (T2-IV) identification algorithm. The observer Markov parameters are then used to compute the local system states at the current instant, which are subsequently used to compute the matrices of the local linear state observer model.

1.1 Contributions

The proposed methodology presents the following main contributions:

• Proposal of an evolving interval type-2 fuzzy state observer modeling approach with non-polarizing consequent estimation.

• Proposal for a novel composition of type-2 fuzzy rules for estimating an uncertain region. The adjustment of the uncertain region is accomplished through a proportional-integral-based adaptation rule.

• Structure learning based on the multiscale approach, where the data normalization is not required in clustering algorithm.

• Novel state-space online nonlinear identification based on interval type-2 fuzzy observer Markov parameters.

2. Overview of interval type-2 state observer fuzzy model

The proposed methodology is based on an evolving interval type-2 state observer fuzzy model (eIT2-SOFM), where its rules can be described as follows:

$$\begin{array}{ll}{\mathbf{\text{Rule}}}^i:\mathbf{\text{IF}}{z}_{1,k}\text{is}{\tilde{Z}}_1^i\mathbf{\text{AND}}\cdots \mathbf{\text{AND}}{z}_{n_z,k}\text{is}{\tilde{Z}}_{n_z}^i& \mathbf{\text{THEN}}\\ & \left\{\begin{array}{c}{\mathbf{x}}_{k+1}^i={\mathbf{A}}^i{\mathbf{x}}_k^i+{\mathbf{B}}^i{\mathbf{u}}_k+{\mathbf{K}}^i{\mathbf{e}}_k\\ {\mathbf{y}}_k^i={\mathbf{C}}^i{\mathbf{x}}_k^i+{\mathbf{D}}^i{\mathbf{u}}_k\end{array}\right.\end{array}$$E1

where i = 1, 2…, c represents the rule number, and $z_{1,k},z_{2,k},\dots ,z_{n_z,k}$ correspond to the antecedent input variables, where $n_z$ is the number of antecedent variables. Additionally, ${\mathbf{A}}^i\in {\Re }^{n\times n}$, ${\mathbf{B}}^i\in {\Re }^{n\times m}$, ${\mathbf{C}}^i\in {\Re }^{p\times n}$, ${\mathbf{D}}^i\in {\Re }^{p\times m}$, and ${\mathbf{K}}^i\in {\Re }^{n\times p}$ represent the state-space matrices of the local linear model for each rule. The local state vector for the i-th rule is denoted as ${\mathbf{x}}_k^i=\left[x_{1,k}^i,x_{2,k}^i,\dots ,x_{n,k}^i\right]\in {\Re }^n$, and the local output vector is ${\mathbf{y}}_k^i=\left[y_{1,k}^i,y_{2,k}^i,\dots ,y_{p,k}^i\right]\in {\Re }^p$. The input signal vector is represented by $\mathbf{u}k=\left[u1,k,u_{2,k},\dots ,u_{m,k}\right]\in {\Re }^m$, and ${\mathbf{e}}_k\in {\Re }^p$ is the error vector given by

$${\mathbf{e}}_k={\mathbf{y}}_k-{\widehat{\mathbf{y}}}_k$$E2

where ${\widehat{\mathbf{y}}}_k$ is the type-0 eIT2-SOFM output estimation in the eIT2-SOFM. An interval type-2 Gaussian membership function is adopted, represented as ${\tilde{\mu }}_j^i=\left[{\overline{\mu }}_j^i,{\underset{¯}{\mu }}_j^i\right]$, with uncertain dispersion denoted as $\tilde{\sigma }=\left[{\overline{\sigma }}^i,{\underset{¯}{\sigma }}^i\right]$, described by

$${\overline{\mu }}_j^i\left(z_{j,k}\right)=exp\left[-\frac{1}{2}{\left(\frac{z_j^{i,\ast }-z_j}{{\overline{\sigma }}_j^i}\right)}^2\right]$$E3

$${\underset{¯}{\mu }}_j^i\left(z_{j,k}\right)=exp\left[-\frac{1}{2}{\left(\frac{z_j^{i,\ast }-z_j}{{\underset{¯}{\sigma }}_j^i}\right)}^2\right]$$E4

where ${\overline{\mu }}_j^i$ and ${\underset{¯}{\mu }}_j^i$ are the upper and lower membership functions, respectively, ${\overline{\sigma }}_j^i$ is the upper dispersion, ${\underset{¯}{\sigma }}_j^i$ is the lower dispersion, and $z_j^{i,\ast }$ is the center of i-th cluster and j-th input axis. The proposed eIT2-SOFM adopts interval output estimation ${\tilde{\mathbf{y}}}_k=\left[{\overline{\mathbf{y}}}_k,{\underset{¯}{\mathbf{y}}}_k\right]$ as the output model. To compute the eIT2-SOFM output, first, the interval firing strength $\left[{\overline{f}}_k^i,{\underset{¯}{f}}_k^i\right]$ is computed as follows:

$${\overline{f}}_k^i=\prod _{j=1}^{n_z}{\overline{\mu }}_j^i\left(z_{j,k}\right)$$E5

$${\underset{¯}{f}}_k^i=\prod _{j=1}^{n_z}{\underset{¯}{\mu }}_j^i\left(z_{j,k}\right)$$E6

and from interval firing strength, the upper and lower normalized firing strength $\left[{\overline{\gamma }}_k^i,{\underset{¯}{\gamma }}_k^i\right]$ is computed as follows:

$${\overline{\gamma }}_k^i=\frac{{\overline{f}}_k^i}{{\sum }_{j=1}^{c_k}{\overline{f}}_k^j}$$E7

$${\underset{¯}{\gamma }}_k^i=\frac{{\underset{¯}{f}}_k^i}{{\sum }_{j=1}^{c_k}{\underset{¯}{f}}_k^j}$$E8

and $\left[{\mathbf{y}}_k^r,{\mathbf{y}}_k^l\right]$ is computed as follows:

$${\mathbf{y}}_k^r=\sum _{i=1}^{c_k}{{\overline{\gamma }}^i}_k{\mathit{y}}_k^i$$E9

$${\mathbf{y}}_k^l=\sum _{i=1}^{c_k}{\underset{¯}{{\gamma }^i}}_k{\mathit{y}}_k^i$$E10

Thus, in the function of $\left[{\mathbf{y}}_k^r,{\mathbf{y}}_k^l\right]$, the upper and lower outputs are computed using the following equations:

$${\overline{\mathbf{y}}}_k=max\left({\mathbf{y}}_k^r,{\mathbf{y}}_k^l\right)+{\dot{\mathbf{y}}}_k$$E11

$${\underset{¯}{\mathbf{y}}}_k=min\left({\mathbf{y}}_k^r,{\mathbf{y}}_k^l\right)-{\dot{\mathbf{y}}}_k$$E12

where ${\dot{\mathbf{y}}}_k$ is the adaptive output degree of uncertainty, which is adjusted based on the digital PI control algorithm, given by

$${\dot{\mathbf{y}}}_k=\left\{\begin{array}{cc}{\dot{\mathbf{y}}}_k\cdot f_f& {\left.\text{for}{\underset{¯}{\mathbf{y}}}_k\le {\mathbf{y}}_k\right|}_f\le {\overline{\mathbf{y}}}_k\\ {\dot{\mathbf{y}}}_k+g_p\mid {\mathbf{e}}_k\mid +g_i\left|{\sum }_{j=k-w}^k{\mathbf{e}}_k\right|& \text{otherwise}\\ & \end{array}\right.$$E13

such that

$${\mathbf{e}}_k={\mathbf{y}}_k-\frac{{\overline{\mathbf{y}}}_k+{\underset{¯}{\mathbf{y}}}_k}{2}$$E14

where $1\le f_f\le 0.90$ is a adjustment factor, $g_p$ and $g_i$ are the proportional gain and integral gain, respectively, w is the window size, and ${\left.{\mathbf{y}}_k\right|}_f$ is the filtered output in instant k computed in the filtering process (Section 3). The proposed incremental learning is performed by the following steps: 1) filtering process, 2) structure learning via evolving method, and 3) submodels updating via type-2 state observer fuzzy identification. In Figure 1, the block diagram of the proposed methodology is shown, where ${\left.{\mathbf{y}}_k\right|}_r$ is the corrupted output data in instant k. In the next sections, the mathematical formulation of each step is presented.

Figure 1.

Block diagram of the proposed methodology. From the dynamic system, output ${\left.{\mathbf{y}}_k\right|}_r$ and input ${\mathbf{u}}_k$ data are obtained and filtered. For the evolving process, the antecedent input vector ${\mathbf{z}}_k$ is generated from the filtered input and/or output, that is, ${\mathbf{z}}_k=\left[y_{k-1}{\left|{}_f\cdots y_{k-l}\right|}_f{u}_{k-1}{\left|{}_f\cdots j_{k-l}\right|}_f\right]$. From vector ${\mathbf{z}}_k$, the evolving mechanism (creation, type-2 adaptation and merging) performs the structure learning, which chances the number of rules in each incoming data. In the sequel, the submodel of each rule is updated by a type-2 fuzzy state observer identification algorithm.

3. Filtering process

Consider a dynamic system where its experimental data are corrupted by correlated noise. For a data-driven learning algorithm, the noise in the database is a problem, once inconsistent cluster (fuzzy partition) and polarized submodel parameter can be computed. Therefore, in the proposed methodology, a filtering process is performed in order to compute data highly correlated with system dynamic and independent of noise. In this step, a recursive moving-average filter are used, where

$$y_k{\left|{}_f={\left(1-a_f\right)}^4{y}_k\right|}_r+4a_f{y}_{k-1}{\left|{}_f-6a_f^2{y}_{k-2}\right|}_f+4a_f^3{y}_{k-3}{\left|{}_f-a_f^4{y}_{k-4}\right|}_f$$E15

where $a_f\in \left(0,1\right)$ is the filtering coefficient chosen by the user. For input data, the filtering process is similar, as follows:

$$u_k{\left|{}_f={\left(1-a_f\right)}^4{u}_k\right|}_r+4a_f{u}_{k-1}{\left|{}_f-6a_f^2{u}_{k-2}\right|}_f+4a_f^3{u}_{k-3}{\left|{}_f-a_f^4{u}_{k-4}\right|}_f$$E16

From ${\left.y_k\right|}_f$ and ${\left.u_k\right|}_f$, the vector ${\mathbf{z}}_k$ is generated, which is used in structure learning, and the regressor vector ${\mathit{\delta }}_k$, which is used for consequent estimation.

4. Structure learning via evolving type-2 fuzzy clustering method

The structure of the eIT2-SOFM is updated with each new incoming dataset, and the adopted learning method does not necessitate prior knowledge. In other words, the rule base initializes with zero rules. The structure learning relies on an evolving type-2 fuzzy clustering (eT2FC), which is employed to create a fuzzy partition in the input variable space. This clustering method projects an interval type-2 fuzzy set onto each input space axis, characterized by interval type-2 Gaussian membership functions with uncertain dispersion. The eT2FC algorithm is based on a multidimensional scaling approach, eliminating the need for data normalization and providing improved handling of non-stationary problems [17].

Initially, for instant k = 1 and the number of rules $c_k=0$, the antecedent input vector ${\mathbf{z}}_k=\left[z_{1,k}{z}_{2,k}\cdots z_{n_z,k}\right]\in {\Re }^{n_z}$ becomes the center of the first cluster, with an initial type-1 dispersion ${\sigma }_0$ defined by the user. The dispersions ${\overline{\sigma }}_j^i$ and ${\underset{¯}{\sigma }}_j^i$ are computed as a function of ${\sigma }^i$ as follows:

$${\overline{\sigma }}_j^i=\sigma +\mathit{\zeta \sigma }$$E17

$${\underset{¯}{\sigma }}_j^i=\sigma -\mathit{\zeta \sigma }$$E18

For $k>1$, the interval type-2 membership values ${\tilde{\mu }}_j^i=\left[{\overline{\mu }}_j^i,{\underset{¯}{\mu }}_j^i\right]$ are computed using Eqs. (3) and (4), and the interval firing strengths $\left[{\overline{f}}_k^i,{\underset{¯}{f}}_k^i\right]$ are calculated using Eqs. (5) and (6). The mean between ${\overline{f}}^i$ and ${\underset{¯}{f}}^i$ is computed by

$$f_k^i=\frac{{\overline{f}}_k^i+{\underset{¯}{f}}_k^i}{2}$$E19

and it is used for the cluster creation (rule creation) mechanism, which described in the sequel.

4.1 Cluster creation rule

To determine the necessity of creating a new cluster, initially, examine the cluster $\chi$ with the highest membership value, that is,

$$\chi =arg\underset{i\in \left[1,c_k\right]}{max}{f}_k^i$$E20

Therefore, the condition for rule creation is defined as follows:

$$\mathbf{\text{IF}}{f}_k^{\chi }<T_f\mathbf{\text{THEN}}{\mathbf{z}}^{c_k+1\ast }={\mathbf{z}}_k$$E21

where $T_f$ represents the firing strength threshold. When the condition for rule creation is satisfied, the vector ${\mathbf{z}}_k$ becomes the center of a new cluster (cluster $c_k+1$). The type-1 dispersion for the new cluster is determined by

$${\sigma }_j^{c_k+1}=\alpha \left|z_{j,k}-z_{j,k}^{\epsilon ,\ast }\right|$$E22

where dispersion ${\overline{\sigma }}_j^i$ and ${\underset{¯}{\sigma }}_j^i$ are computed by Eqs. (17) and (18), respectively.

4.2 Merging mechanism

Once the creation rule (21) is satisfied, the merging condition is checked. This mechanism verifies if the new membership function ${\tilde{\mu }}_j^{c_k+1}$ is redundant. First, it determines the closest membership function to ${\tilde{\mu }}_j^{c_k+1}$, i.e.,

$$\epsilon =arg\underset{i\in \left[1,c_k\right]}{max}exp\left[-\frac{1}{2}{\left(\frac{z_j^{i,\ast }-z_j}{{\sigma }_j^i}\right)}^2\right]$$E23

where $i\ne c_k+1$ and $\epsilon$ is the index of the closest membership function. Therefore, for the membership functions $c_k+1$ and $\epsilon$ along the j-th axis, the similarity degree is verified as follows:

$$S\left(z_j^{c_k+1,\ast },z_j^{\epsilon ,\ast }\right)=max\left({\mu }_j^{c_k+1}\left(z_j^{\epsilon ,\ast }\right),{\mu }_j^{\epsilon }\left(z_j^{c_k+1,\ast }\right)\right)$$E24

where ${\mu }_j^{c_k+1}\left(z_j^{\chi ,\ast }\right)$ and ${\mu }_j^{\chi }\left(z_j^{c_k+1,\ast }\right)$ are computed as follows:

$${\mu }_j^{c_k+1}\left(z_j^{\epsilon ,\ast }\right)=exp\left[-\frac{1}{2}{\left(\frac{z_j^{\epsilon ,\ast }-z_j^{c_k+1,\ast }}{{\sigma }_j^{\epsilon }}\right)}^2\right]$$E25

$${\mu }_j^{\epsilon }\left(z_j^{c_k+1,\ast }\right)=exp\left[-\frac{1}{2}{\left(\frac{z_j^{c_k+1,\ast }-z_j^{\epsilon ,\ast }}{{\sigma }_j^{c_k+1}}\right)}^2\right]$$E26

From a upper threshold $T_u$ and lower threshold $T_l$ defined by user, the following conditions are verified:

1. If $S>T_u$, the new membership function is replaced by ${\mu }_j^{c_k+1}\left(z_j^{\epsilon ,\ast }\right)$.

2. If $T_l<S<T_u$, the two membership function must be merged.

3. If $S<T_l$, The new membership function is maintained

If condition 2 is satisfied, the following equations are used for computing the new center and new dispersion

$$z^{\mathit{new}}=\frac{z_j^{c_k+1}+N_p^i{z}_j^i}{1+N_p^i}$$E27

$${\sigma }^{\mathit{new}}=\frac{{\sigma }_j^{c_k+1}+{\sigma }_j^i}{\sqrt{\pi }}$$E28

where $N_p^i$ is the number of points (z) associated with the cluster i.

4.3 Cluster adaptation mechanism

In the antecedent parameters adaptation, the approach adopted is to update the cluster center with the highest membership value when the new cluster condition is not satisfied. Thus, the updating of the center ${\mathbf{z}}^{\chi ,\ast }$ is given by

$$\mathrm{\Delta }z=\frac{z_{j,k}^{\chi ,\ast }{N}_p^{\chi }}{N_p^{\chi }+1}+\frac{z_{j,k}}{N_p^{\chi }+1}$$E29

$$z_{j,k+1}^{\chi ,\ast }=z_{j,k}^{\chi ,\ast }+\rho \left(\mathrm{\Delta }z-z_{j,k}^{\chi ,\ast }\right)$$E30

where $\rho$ is a learning rate defined by user and $\mathrm{\Delta }z$ is the center adjustment.

5. Submodels updating via IV-based type-2 fuzzy state observer identification

The state-space equations are regarded as the consequent part in the proposed eIT2-SOFM. The estimation of the matrices ${\mathbf{A}}^i,{\mathbf{B}}^i,{\mathbf{C}}^i,{\mathbf{D}}^i,$ and ${\mathbf{K}}^i$ is based on the fuzzy observer Markov parameters, where its mathematical foundations for the type-1 version are detailed in the works of [13, 18, 19]. This paper presents a type-2 state-space fuzzy modeling approach that utilizes the observer Markov parameters estimated through the IV fuzzy method. The subspace approach is employed to estimate the matrices ${\mathbf{A}}^i,{\mathbf{B}}^i,{\mathbf{C}}^i,{\mathbf{D}}^i,$ and ${\mathbf{K}}^i$ from the observer Markov parameters of each rule.

5.1 Mathematical definition of type-2 fuzzy observer markov parameters

Considering a vector auto-regressive model to estimate the local state-space model, as follows:

$${\mathbf{y}}_k^i=\sum _{j=0}^{q_p}{\dot{\mathbf{\Xi }}}_{k-j}^{i,\left(u\right)}{\mathbf{u}}_{k-j}+\sum _{j=1}^{q_p}{\dot{\mathbf{\Xi }}}_{k-j}^{i,\left(y\right)}{\mathbf{y}}_{k-j}^i$$E31

According to [20], if the local state-space model is asymptotically stable, the matrices ${\dot{\mathbf{\Xi }}}_{k-j}^{i,\left(u\right)}$ and ${\dot{\mathbf{\Xi }}}_{k-j}^{i,\left(y\right)}$ in Eq. (31) are the observer Markov parameters of i-th local state-space model, being

$${\dot{\mathbf{\Xi }}}_{k-j}^{i,\left(u\right)}=\left\{\begin{array}{cc}{\mathbf{D}}^i& \text{if}j=0\\ {\mathbf{C}}^i{\left({\mathbf{A}}^i\right)}^{j-1}{\mathbf{B}}^i& \text{if}j>0\end{array}\right.$$E32

$${\dot{\mathbf{\Xi }}}_{k-j}^{i,\left(y\right)}={\mathbf{C}}^i{\left({\mathbf{A}}^i\right)}^{j-1}{\mathbf{K}}^i$$E33

where ${\mathbf{K}}^i\in {\Re }^{n\times p}$ is the local observer matrix [13]. Thus, the matrix composed by the observer Markov parameters matrix is given by

$${\dot{\mathbf{\Xi }}}^i=\left[{\dot{\mathbf{\Xi }}}_{k-q_p}^{i,\left(u\right)},\cdots ,{\dot{\mathbf{\Xi }}}_k^{i,\left(u\right)},{\dot{\mathbf{\Xi }}}_{k-q_p}^{i,\left(y\right)}\cdots {\dot{\mathbf{\Xi }}}_{k-1}^{i,\left(y\right)}\right]$$E34

and Eq. (31) is rewritten in matrix form as follows:

$${\left({\mathbf{y}}_k^i\right)}^T={\left({\mathit{\delta }}_k^i\right)}^T{\left({\dot{\mathbf{\Xi }}}^i\right)}^T+{\mathit{\xi }}_k$$E35

with ${\left({\mathit{\delta }}_k^i\right)}^T=\left[{\stackrel{⌣}{\mathbf{u}}}_{k-q_p}^T{\mathbf{u}}_k^T{\left({\stackrel{⌣}{\mathbf{y}}}_{k-q_p}^i\right)}^T\right]$, such that

$${\stackrel{⌣}{\mathbf{u}}}_{k-q_p}=\left[\begin{array}{c}{\mathbf{u}}_{k-q_p}\\ {\mathbf{u}}_{k-q_p+1}\\ ⋮\\ {\mathbf{u}}_{k-1}\end{array}\right]{\stackrel{⌣}{\mathbf{y}}}_{k-q_p}^i=\left[\begin{array}{c}{\mathbf{y}}_{k-q_p}^i\\ {\mathbf{y}}_{k-q_p-1}^i\\ ⋮\\ {\mathbf{y}}_{k-1}^i\end{array}\right]$$E36

Assuming k samples, where $q_p$ is the past time-window and $k>q_p$, from Eq. (35), the following batch equation is derived:

$${\mathbf{Y}}_k^i\mathbf{=}{\mathbf{\Delta }}_k^i{\left({\dot{\mathbf{\Xi }}}^i\right)}^T$$E37

where

$${\mathbf{Y}}_k^i=\left[\begin{array}{c}{\left({\mathbf{y}}_{q_p+1}^i\right)}^T\\ {\left({\mathbf{y}}_{q_p+2}^i\right)}^T\\ ⋮\\ {\left({\mathbf{y}}_k^i\right)}^T\end{array}\right],{\mathbf{\Delta }}_k^i=\left[\begin{array}{c}{\left({\mathit{\delta }}_{q_p+1}^i\right)}^T\\ {\left({\mathit{\delta }}_{q_p+2}^i\right)}^T\\ ⋮\\ {\left({\mathit{\delta }}_k^i\right)}^T\end{array}\right]$$E38

Thus, considering the TSK fuzzy theory, the batch computation fuzzy model output history is given by

$${\mathbf{Y}}_k=\sum _{i=1}^c{\tilde{\mathbf{\Gamma }}}_k^i{\mathbf{\Delta }}_k^i{\left({\dot{\mathbf{\Xi }}}^i\right)}^T$$E39

where ${\tilde{\mathbf{\Gamma }}}_k^i=diag\left(\left[{\gamma }_{q_p+1}^i,{\tilde{\gamma }}_{q_p+2}^i,\cdots ,{\tilde{\gamma }}_k^i\right]\right)$, so that ${\tilde{\gamma }}_{q_p+1}^i$ is type-2 normalized firing strength computed by

$${\tilde{\gamma }}_k^i=\frac{{\overline{f}}_k^i+{\underset{¯}{f}}_k^i}{{\sum }_{j=1}^{c_k}\left({\overline{f}}_k^j+{\underset{¯}{f}}_k^j\right)}$$E40

Thus, from the batch estimation approach, the outputs from instant $q_p$ to k can be computed by:

$${\tilde{\mathbf{\Gamma }}}_k^i{\mathbf{Y}}_k\mathbf{=}{\tilde{\mathbf{\Gamma }}}_k^i{\mathbf{\Delta }}_k{\left({\dot{\mathbf{\Xi }}}^i\right)}^T$$E41

where

$${\mathbf{Y}}_k=\left[\begin{array}{c}{\left({\mathbf{y}}_{q_p+1}\right)}^T\\ {\left({\mathbf{y}}_{q+2}\right)}^T\\ ⋮\\ {\left({\mathbf{y}}_k\right)}^T\end{array}\right]$$E42

5.2 Fuzzy observer markov parameters estimation via IV approach

Assuming experimental data are corrupted by correlated noise, the vector ${\mathit{\delta }}_k$ presents noisy data, that is,

$${\mathit{\delta }}_k^r={\mathit{\delta }}_k^i+{\mathit{\upsilon }}_k$$E43

where ${\mathit{\upsilon }}_k$ is the correlated noise vector related to ${\mathit{\delta }}_k$.

According to literature [21, 22], from an instrument vector ${\left({\mathit{\delta }}_k^f\right)}^T$, which is highly correlated with output and/or input data and not correlated with the noise, the following solution for ${\dot{\mathbf{\Xi }}}^i$, from Eq. (41), is derived:

$${\left({\dot{\mathbf{\Xi }}}^i\right)}^T={\left({\mathbf{\Delta }}_k^f{\tilde{\mathbf{\Gamma }}}_k^i{\mathbf{\Delta }}_k^r\right)}^{-1}{\mathbf{\Delta }}_k^f{\tilde{\mathbf{\Gamma }}}_k^i{\mathbf{Y}}_k^f$$E44

where ${\mathbf{\Delta }}_k^f\in {\Re }^{\left(k-q_p-1\right)\times \left(q_p\left(m+p\right)+m\right)}$ is the batch instruments matrix. Extended the recursive type-1 fuzzy IV algorithm presented in [21, 22] for the type-2 case, the updation of the interval type-2 fuzzy observer Markov parameters is performed by following equations:

$${\mathbf{L}}_{k+1}^i={\tilde{\gamma }}_k^i\frac{{\mathbf{P}}_k^i{\mathit{\delta }}_k^f}{{\tilde{\gamma }}_k^i+{\left({\mathit{\delta }}_k^r\right)}^T{\mathbf{P}}_k^i{\mathit{\delta }}_k^r}$$E45

$${\dot{\mathbf{\Xi }}}_{k+1}^i\mathbf{=}{\dot{\mathbf{\Xi }}}_k^i+{\tilde{\gamma }}_k^i{\mathbf{e}}_k{\mathbf{L}}_{k+1}^i$$E46

$${\mathbf{P}}_{k+1}^i=\frac{1}{\beta }\left[\mathbf{I}-{\mathbf{L}}_{k+1}^i{\left({\mathit{\delta }}_k^r\right)}^T\right]{\mathbf{P}}_k^i$$E47

where $1\le \beta \le 0.9$ is the forgetting factor, L is the IV gain matrix, and P is the IV variance matrix.

5.3 Computation of space state matrices

According to [23, 24], the local state vector ${\mathbf{x}}_k^i$ can be computed by the following close-formula:

$${\mathbf{x}}_k^i=\mathbf{S}\left[{\mathbf{\Lambda }}_k^i{\overline{\mathbf{u}}}_{k-q_p}^f+{\mathbf{\Upsilon }}_k^i{\overline{\mathbf{y}}}_{k-q_p}^f\right]$$E48

where $\mathbf{S}\in {\Re }^{n\times \left(\mathit{mf}-n\right)}$ is a positive-defined matrix, and ${\mathbf{\Lambda }}_k^i$ and ${\mathbf{\Upsilon }}_k^i$ are formed by the observer Markov parameters, such that

$${\mathbf{\Lambda }}_k^i=\left[\begin{array}{cccccc}{\dot{\mathbf{\Xi }}}_{k-q}^{i,\left(u\right)}& {\dot{\mathbf{\Xi }}}_{k-q_p+1}^{i,\left(u\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-q_p+q_f-1}^{i,\left(u\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-1}^{i,\left(u\right)}\\ \mathbf{0}& {\dot{\mathbf{\Xi }}}_{k-q_p}^{i,\left(u\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-q_p+q_f-2}^{i,\left(u\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-2}^{i,\left(u\right)}\\ ⋮& \ddots & \ddots & ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & \mathbf{0}& {\dot{\mathbf{\Xi }}}_{k-q_p}^{i,\left(u\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-q_f}^{i,\left(u\right)}\end{array}\right]$$E49

and

$${\mathbf{\Upsilon }}_k^i=\left[\begin{array}{cccccc}{\dot{\mathbf{\Xi }}}_{k-q}^{i,\left(y\right)}& {\dot{\mathbf{\Xi }}}_{k-q_p+1}^{i,\left(y\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-q_p+q_f-1}^{i,\left(y\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-1}^{i,\left(y\right)}\\ \mathbf{0}& {\dot{\mathbf{\Xi }}}_{k-q_p}^{i,\left(y\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-q_p+q_f-2}^{i,\left(y\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-2}^{i,\left(y\right)}\\ ⋮& \ddots & \ddots & ⋮& \ddots & ⋮\\ \mathbf{0}& \cdots & \mathbf{0}& {\dot{\mathbf{\Xi }}}_{k-q_p}^{i,\left(y\right)}& \cdots & {\dot{\mathbf{\Xi }}}_{k-q_f}^{i,\left(y\right)}\end{array}\right]$$E50

Once the local state vectors ${\mathbf{x}}_k^i$ are computed, the fuzzy state vector ${\tilde{\mathbf{x}}}_k$ is obtained as follows:

$${\tilde{\mathbf{x}}}_k=\sum _{i=1}^{c_k}{\tilde{\gamma }}_k^i{\mathbf{x}}_k^i$$E51

From Eq. (51), the matrices $\left[{\mathbf{A}}_k^i,{\mathbf{B}}_k^i,{\mathbf{K}}_k^i\right]$ can be estimated using the QR solution. Thus, the state equation is formulated as follows:

$${\left({\mathbf{x}}_k^i\right)}^T=\left[\begin{array}{ccc}{\left({\mathbf{x}}_{k-1}^i\right)}^T& {\mathbf{u}}_{k-1}^T& {\mathbf{e}}_{k-1}^T\end{array}\right]\left[\begin{array}{c}{\left({\mathbf{A}}_k^i\right)}^T\\ {\left({\mathbf{B}}_k^i\right)}^T\\ {\left({\mathbf{K}}_k^i\right)}^T\end{array}\right]={\left({\mathit{\nu }}_{k-1}^{x,i}\right)}^T{\left({\mathbf{\Theta }}_k^{i,x}\right)}^T$$E52

Let the least square solution of Eq. (52), given by:

$${\left({\mathbf{\Theta }}_k^{i,x}\right)}^T={\left({\mathbf{V}}_k{\tilde{\mathbf{\Gamma }}}_k^i{\mathbf{V}}_k\right)}^{-1}{\mathbf{V}}_k{\tilde{\mathbf{\Gamma }}}_k^i{\mathbf{X}}_k$$E53

$${\left({\mathbf{\Theta }}_k^{i,x}\right)}^T=({\mathbf{P}}_k^{{\mathbf{\Theta }}_x^i}{\left|{}_{\mathit{\nu \nu }}){}^{-1}{\mathbf{P}}_k^{{\mathbf{\Theta }}_x^i}\right|}_{\nu y^f}$$E54

where

$${\mathbf{X}}_k=\left[\begin{array}{c}{\mathbf{x}}_{k-q_p}^T\\ {\mathbf{x}}_{k-q_p+1}^T\\ ⋮\\ {\mathbf{x}}_k^T\end{array}\right]$$E55

Rewritten Eq. (53), it has

$${\mathbf{P}}_k^{{\mathbf{\Theta }}_x^i}{\left|{}_{\mathit{\nu \nu }}{\left({\mathbf{\Theta }}^{i,x}\right)}^T={\mathbf{P}}_k^{{\mathbf{\Theta }}_x^i}\right|}_{\nu y^f}$$E56

where the following recursion is derived:

$${\mathbf{P}}_k^{{\mathbf{\Theta }}_x^i}{\left|{}_{\mathit{\nu \nu }}={\mathbf{P}}_{k-1}^{{\mathbf{\Theta }}_x^i}\right|}_{\mathit{\nu \nu }}+{\tilde{\gamma }}_k^i{\mathit{\nu }}_{k-1}{\left({\mathit{\nu }}_{k-1}^{x,i}\right)}^T$$E57

$${\mathbf{P}}_k^{{\mathbf{\Theta }}_x^i}{\left|{}_{\mathit{\nu y}}={\mathbf{P}}_{k-1}^{{\mathbf{\Theta }}_x^i}\right|}_{\mathit{\nu y}}+{\tilde{\gamma }}_k^i{\mathit{\nu }}_{k-1}^{x,i}{\left({\mathbf{y}}_{k-1}^f\right)}^T$$E58

Thus, by applying QR factorization to ${\mathbf{P}}^{{\mathbf{\Theta }}_x^i}k\mid \mathit{\nu \nu }$, it has

$${\left.\mathbf{R}{\left({\mathbf{\Theta }}^{i,x}\right)}^T={\mathbf{Q}}^T{\mathbf{P}}_k^{{\mathbf{\Theta }}_x^i}\right|}_{\nu y^f}$$E59

The matrix ${\mathbf{\Theta }}^{i,x}$ is computed through the backward substitution method in Eq. (59) [21]. To compute the $\left[{\mathbf{C}}_k^i,{\mathbf{D}}_k^i\right]$ matrices, the output equation can be formulated as follows:

$${\left({\mathbf{y}}_k^i\right)}^T=\left[\begin{array}{cc}{\left({\mathbf{x}}_k^i\right)}^T& {\mathbf{u}}_k^T\end{array}\right]\left[\begin{array}{c}{\left({\mathbf{C}}_k^i\right)}^T\\ {\left({\mathbf{D}}_k^i\right)}^T\end{array}\right]={\left({\mathit{\nu }}_{k-1}^{y,i}\right)}^T{\left({\mathbf{\Theta }}_k^{i,y}\right)}^T$$E60

Using the same steps for computing ${\left({\mathbf{\Theta }}_k^{i,x}\right)}^T$, it has

$${\mathbf{P}}_k^{{\mathbf{\Theta }}_y^i}{\left|{}_{\mathit{\nu \nu }}{\left({\mathbf{\Theta }}^{i,y}\right)}^T={\mathbf{P}}_k^{{\mathbf{\Theta }}_y^i}\right|}_{\nu y^f}$$E61

such that

$${\mathbf{P}}_k^{{\mathbf{\Theta }}_y^i}{\left|{}_{\mathit{\nu \nu }}={\mathbf{P}}_{k-1}^{{\mathbf{\Theta }}_y^i}\right|}_{\mathit{\nu \nu }}+{\tilde{\gamma }}_k^i{\mathit{\nu }}_{k-1}^{y,i}{\left({\mathit{\nu }}_{k-1}^{y,i}\right)}^T$$E62

$${\mathbf{P}}_k^{{\mathbf{\Theta }}_y^i}{\left|{}_{\mathit{\nu y}}={\mathbf{P}}_{k-1}^{{\mathbf{\Theta }}_y^i}\right|}_{\mathit{\nu y}}+{\tilde{\gamma }}_k^i{\mathit{\nu }}_{k-1}^{y,i}{\left({\mathbf{y}}_{k-1}^f\right)}^T$$E63

where ${\mathbf{\Theta }}_k^{i,y}$ can be computed by applying QR factorization to ${\mathbf{P}}^{{\mathbf{\Theta }}_y^i}k\mid \mathit{\nu \nu }$, such that

$${\left.\mathbf{R}{\left({\mathbf{\Theta }}^{i,y}\right)}^T={\mathbf{Q}}^T{\mathbf{P}}_k^{{\mathbf{\Theta }}_y^i}\right|}_{\nu y^f}$$E64

followed by backward substitution applied in Eq. (64).

6. Computational results

Considering the mathematical formulation of the proposed algorithm described in Sections 4 and 5, two case studies are presented: the identification of a SISO nonlinear system and online identification of a time-varying MIMO dynamic system. To validate the results, the following metrics were used:

• Non-dimensional error index (NDEI):

  $$\text{NDEI}=\frac{\sqrt{\frac{1}{N}{\sum }_{k=1}^N{\left({\tilde{e}}_{r,k}\right)}^2}}{\mathit{std}\left(\mathbf{y}\right)}$$E65

  where $\mathit{std}\left(•\right)$ is the standard deviation.

• Variance accounted for (VAF%):

  $$\mathrm{VAF}\left(\%\right)=\left[1-\frac{{\tilde{\mathbf{e}}}_r}{\mathit{var}\left(\mathbf{y}\right)}\right]\times 100$$E66

  where $\mathit{var}\left(•\right)$ is the variance.

  where ${\tilde{e}}_k$ is the confidence region error for interval type-2 estimation, which is described as follows:

$${\tilde{e}}_k=\left\{\begin{array}{cc}0& \text{if}{\underset{¯}{y}}_k<y_k<{\widehat{\overline{y}}}_k\\ y_k-{\widehat{\overline{y}}}_k& \text{if}{\overline{e}}_k<{\underset{¯}{e}}_k\\ y_k-{\underset{¯}{y}}_k& \text{otherwise}\end{array}\right.$$E67

such that

$${\overline{e}}_k=\mid y_k-{\overline{\widehat{y}}}_k\mid ,{\underset{¯}{e}}_k=\mid y_k-{\underset{¯}{\widehat{y}}}_k\mid$$E68

6.1 Nonlinear dynamic system

The identification problem under consideration is a SISO nonlinear dynamic system, commonly utilized as a benchmark in the type-2 fuzzy modeling literature. It is described by the following equation:

$$y_k=\frac{y_{k-1}{y}_{k-2}\left(y_{k-1}+0.5\right)}{1+y_{k-1}^2+y_{k-2}^2}+u_{k-1}$$E69

where the input signal is given by $u_k=sin\left(\frac{2\mathit{k\pi }}{25}\right)$. For the identification process, a dataset consisting of 1300 samples was generated. Among these samples, 1000 were allocated for the training step, while 300 were used for the validation step. The algorithm parameters were configured with the following values: $a_f=0$; $T_f=0.002$; $T_u=0.7$; $T_l=0.5$; $q_p=11$; $q_f=6$; $n=2$; $w=5$; $g_p={10}^{-2}$; $f_f=0.99$; and $g_i={10}^{-3}$. The rule structure adopted in this experiment is

$$\begin{array}{ll}{\mathbf{\text{Rule}}}^i:\mathbf{\text{IF}}{z}_{1,k}\text{is}{\tilde{Z}}_1^i\mathbf{\text{AND}}{z}_{2,k}\text{is}{\tilde{Z}}_2^i& \mathbf{\text{THEN}}\\ & \left\{\begin{array}{c}{\mathbf{x}}_{k+1}^i={\mathbf{A}}^i{\mathbf{x}}_k^i+{\mathbf{B}}^i{\mathbf{u}}_k+{\mathbf{K}}^i{\mathbf{e}}_k\\ {\mathbf{y}}_k^i={\mathbf{C}}^i{\mathbf{x}}_k^i+{\mathbf{D}}^i{\mathbf{u}}_k\end{array}\right.\end{array}$$E70

where $z_{1,k}=u_k$ and $z_{2,k}=y_k$.

For comparative analysis, the models eTS [25], xTS (cited in [26]), DENFIS [27], eTF [28], eMG [29], and RIV-NFM [22] are considered. The performances, as assessed by the NDEI metric, are presented in Table 1. Figure 2 illustrates the uncertain region estimated by eIT2-SOFM for the validation dataset.

Model	Rules	NDEI
eTS [25]	7	0.1038
xTS cited in [26]	7	0.0936
DENFIS [27]	7	0.0842
eFT [28]	7	0.0653
eMG [29]	7	0.0501
RIV-NFM [22]	6	0.0413
Proposed	3	0.0203

Table 1.

Comparative analysis of the proposed methodology with other relevant methodologies for the nonlinear dynamic system problem.

Figure 2.

Uncertain region estimation for nonlinear dynamic system identification.

6.2 Time-varying MIMO dynamic system

A time-varying nonlinear MIMO dynamic system is considered to demonstrate the adaptability of the proposed methodology for time-varying dynamic systems. The nonlinear MIMO dynamic system is described by the following equations:

$${\mathbf{v}}_k=\left\{\begin{array}{c}{v}_{1,k+1}=\frac{v_{1,k-1}}{v_{1,k}^2+1}+\frac{1}{v_{2,k}^2+5}+u_{1,k}^3\\ v_{2,k+1}=0.1v_{2,k}{v}_{1,k-1}-0.2v_{2,k}{u}_{2,k}\end{array}\right.$$E71

$${\mathbf{y}}_k={\mathbf{G}}_k{\mathbf{v}}_k$$E72

where

$${\mathbf{G}}_k=\left\{\begin{array}{c}\left[\begin{array}{cc}1& -1\\ -1& 2\end{array}\right]\text{for}k<1400\\ \left[\begin{array}{cc}3& 2\\ -1& 5\end{array}\right]\text{otherwise}\end{array}\right.$$E73

and ${\mathbf{u}}_k={\left[u_{1,k}{u}_{2,k}\right]}^T$, where $u_{1,k}$ is a multistep signal with a uniform distribution between [−2, 2], $u_{2,k}$ is a multistep signal with a uniform distribution between [−1 1], and ${\mathbf{v}}_k={\left[v_{1,k}{v}_{2,k}\right]}^T$ represents the vector of intermediate signals. The outputs signals were corrupted by correlated noises, which are given by

$${\nu }_k^{y_1}=\frac{1+0.2z^{-1}}{1+0.6z^{-1}+0.2z^{-2}}{e}_k$$E74

$${\nu }_k^{y_2}=\frac{1+0.2z^{-1}}{1+0.3z^{-1}+0.1z^{-2}}{e}_k$$E75

where $e_k$ represents white noise with a mean of zero and a variance of ${\sigma }_e^2$. The dataset comprises 2800 samples, with 500 used for initializing the eIT2-SOFM and 2400 samples used for online identification.

The algorithm parameters were set to the following values: $a_f=0.9$; $T_f=0.001$; $T_u=0.7$; $T_l=0.5$; $q_p=12$; $q_f=7$; $n=4$; $w=1$; $g_p=5\times {10}^{-3}$; $f_f=0.985$; and $g_i=3\times {10}^{-5}$. The rule structure adopted in this experiment is

$$\begin{array}{ll}{\mathbf{\text{Rule}}}^i:\mathbf{\text{IF}}{z}_{1,k}\text{is}{\tilde{Z}}_1^i\mathbf{\text{AND}}{z}_{2,k}\text{is}{\tilde{Z}}_2^i& \mathbf{\text{AND}}{z}_{3,k}\text{is}{\tilde{Z}}_3^i\mathbf{\text{AND}}{z}_{4,k}\text{is}{\tilde{Z}}_4^i\\ & \mathbf{\text{THEN}}\left\{\begin{array}{l}{\mathbf{x}}_{k+1}^i={\mathbf{A}}^i{\mathbf{x}}_k^i+{\mathbf{B}}^i{\mathbf{u}}_k+{\mathbf{K}}^i{\mathbf{e}}_k\\ {\mathbf{y}}_k^i={\mathbf{C}}^i{\mathbf{x}}_k^i+{\mathbf{D}}^i{\mathbf{u}}_k\end{array}\right.\end{array}$$E76

where $z_{1,k}=u_{1,k-1}$, $z_{2,k}=u_{2,k-1}$, $z_{3,k}=y_{1,k-1}$, and $z_{4,k}=y_{2,k-1}$.

In this case, the Monte Carlo method was employed, involving 50 experiment realizations to compute the means of the VAF% and NDEI criteria. The interval estimations of $y_1$ and $y_2$ according to the SNR variation, are shown in Table 2. The online estimations of the time-varying nonlinear MIMO dynamic system outputs, are shown in Figure 3, with an SNR of 10 dB.

	VAF% (±std^a)		NDEI
SNR	$y_1$	$y_2$	$y_1$	$y_2$
10	95.83 (±1.34%)	95.10 (±1.91%)	0.210	0.202
15	96.34 (±0.86%)	95.61 (±0.92%)	0.202	0.191
20	96.62 (±0.57%)	95.80 (±0.58%)	0.193	0.188
25	97.01 (±0.41%)	95.93 (±0.53%)	0.186	0.181
30	97.27 (±0.23%)	96.16 (±0.44%)	0.179	0.169

Table 2.

Estimation performance of the proposed methodology for different levels of noise for VAF and NDEI metrics in time-varying nolinear MIMO dynamic system modeling.

^astd: standard deviation.

Figure 3.

Upper and lower estimation of time-varying nonlinear MIMO dynamic system: (a) $y_1$ and b) $y_2$. The SNR for this experiment was 10 dB. It is noted that the estimation accuracy of the proposed methodology even in noise environments. The purple region was estimated by the eIT2-SOFM based on the experimental data.

7. Discussions

In this paper, aspects of the proposed methodology were presented. The application of eIT2-SOFM for the nonlinear identification of SISO and MIMO dynamic systems was discussed. In Section 6.1, a modeling benchmark problem was used to compare the performance of eIT2-SOFM with other methodologies presented in the literature. Upon reviewing Table 1, it becomes evident that significantly improved results are achieved by the proposed methodology. Additionally, it is noteworthy that only three rules were created by eIT2-SOFM during the identification process, making it the model with the fewest number of rules among the compared methodologies. This result highlights the methodology’s ability to track the nonlinear behavior of dynamic systems.

In Section 6.2, a case study involving a nonlinear MIMO system was presented to demonstrate the tracking capabilities of the proposed methodology in dealing with time-varying problems. The performance results, as shown in Table 2, indicate that the eIT2-SOFM achieved a performance exceeding 95% for each SNR value. This underlines the adaptability of the proposed learning algorithm, even in the presence of correlated noise within the dataset.

8. Conclusions

Considering the experimental results and the methodological aspects of the proposed modeling approach based on eIT2-SOFM, the following concluding remarks are made:

• The proposed method demonstrates robustness to outliers and noise through the incorporation of a filtering process, type-2 fuzzy sets, and the T2-IV algorithm. The filtering process precedes the structure learning step to prevent the creation of nonrelevant rules, while the T2-IV algorithm provides a nonpolarized estimation of the local state observer model parameters.

• Numerical robustness is ensured, as the QR-decomposition is applied to compute the local state observer models.

• The computational results have demonstrated that the proposed methodology is effective for modeling complex dynamic systems characterized by uncertainty, nonlinearity, and both single and multivariable aspects, even in the presence of colored noise.

Among practical projects and problems that can be solved by the algorithm, the following has been widely considered for research:

• Black-box model-based control, where the plant presents nonlinearity, uncertain behavior, and correlated noise, as satellite positioning [30], multimobile manipulator, and induction motor.

• Computational modeling of experimental data, where the data are nonlinear, uncertain, and/or corrupted by correlated noise, such as parameter estimation of vehicle dynamics, mechatronic systems, mobile robot navigation, and nonstationary processes [31].