Neural Computing Enhanced Parameter Estimation for Multi-Input and Multi-Output Total Non-Linear Dynamic Models

Abstract

In this paper, a gradient descent algorithm is proposed for the parameter estimation of multi-input and multi-output (MIMO) total non-linear dynamic models. Firstly, the MIMO total non-linear model is mapped to a non-completely connected feedforward neural network, that is, the parameters of the total non-linear model are mapped to the connection weights of the neural network. Then, based on the minimization of network error, a weight-updating algorithm, that is, an estimation algorithm of model parameters, is proposed with the convergence conditions of a non-completely connected feedforward network. In further determining the variables of the model set, a method of model structure detection is proposed for selecting a group of important items from the whole variable candidate set. In order to verify the usefulness of the parameter identification process, we provide a virtual bench test example for the numerical analysis and user-friendly instructions for potential applications.

1. Introduction

Because a total non-linear model can provide a very concise representation for complex non-linear systems and has good extrapolation characteristics, it has attracted the attention of academic research and applications. Compared with the polynomial non-linear auto-regressive moving average with exogenous input (NARMAX) model, the total non-linear model is an extension of the polynomial model, which can be defined as the ratio of two polynomial expressions [1,2,3]. The introduction of denominator polynomials makes the NARMAX model non-linear in parameters and regression terms. Therefore, compared with the polynomial model, the model identification and the controller design of the total non-linear model are much more challenging [4,5]. In view of the difficulty of parameter estimation of a total non-linear model, using simple and effective algorithm and machine learning should be considered for extracting the information from measurement data.

1.1. Literature Survey

At present, a variety of model structure detection techniques and parameter estimation algorithms are developed for non-linear models, including the orthogonal model structure detection and parameter estimation program [6], the generalized least square estimator [7,8], the prediction error estimator [9,10], the Kalman filter estimator [11,12], the genetic algorithm estimator [12,13], the artificial neural network estimator [14,15,16,17], etc. However, most of these algorithms are parameter estimators for polynomial non-linear models. Zhu and Billings have done a lot of research work on the parameter identification of a total non-linear model [7,8], and they put forward the parameter estimation method of a total non-linear model based on a back-propagation (BP) algorithm in 2003. They discussed the advantages of BP calculation in recognition of the classical model to provide the best combination of classical and neural network methods and provided a powerful tool for analyzing a large number of systems.

In [18], a back-propagation estimation formula based on neuro-computing was presented for estimating the total non-linear model parameters, where a pack of solutions were derived for the problems of parameter initialization, learning rate selection, stop criteria and model structure detection, and the convergence of a back-propagation estimator (BPE). However, Reference [18] only proposed a parameter estimation method for single-input and single-output (SISO) systems, and correspondingly the case studies. Expanding [18], this paper presents solutions for the parameter estimation of a total non-linear multi-input and multi-output (MIMO) model. Due to the complexity of a MIMO system, it is more difficult to estimate the model parameters, but they are more general in academic research and applications. For example, the parameters of a MIMO system are many more than those of a SISO system, and the parameters to be estimated each time will be multiplied, which increases the difficulty of estimation. Moreover, due to the coupling of multiple systems, the parameter values of each system also affect each other. The algorithms to estimate these parameters are not independent but interactive and complex. Because the components of different MIMO systems are different, the total connection neural network structure adopted in [18] is not suitable for estimating the parameters of MIMO systems. When the MIMO system is mapped into a neural network, the network structure is often asymmetric or non-completely connected (the neurons in the hidden layer are not connected with all the neurons in the input layer). That is to say, the network is not a common completely connected feedforward neural network, and the general BP algorithm cannot be directly applied to the estimation of the parameters. Therefore, the learning algorithm of the parameters must be properly derived. Due to the asymmetry of the network, the convergence of the network is also facing challenges. It is necessary to analyze the convergence of the network and give the specific conditions of the network convergence. A MIMO system needs to identify the parameters of a SISO system several times, and a MIMO system can have multiple inputs. In the simulation experiment, the parameters of the system should be estimated under different combinations of multiple inputs, and the performance of the network estimator should be verified. Therefore, the parameter identification of a MIMO system is much more challenging.

1.2. Motivation and Contributions

The authors of [19] presented a thorough analysis that included two kernel components, the SISO model and the orthogonal algorithms are parameter estimators for polynomial non-linear models such as predictive and back propagation computation. Since then, rational model identification has gone to diversified directions, such as more theoretical considerations of a non-linear least squares algorithm [4], a maximum likelihood estimation [3], and a biased compensation recursive least squares algorithm [2]. It has been noted that the MIMO rational model identification has seldom attracted research, probably due to the complexity in algorithm formulation and the coupling effect. However, this MIMO rational model identification should be a research agent now because of recent applications and increasing computing facilities.

The total non-linear system model, which is relatively new, is the alternative name of the NARMAX rational model, which was defined by a survey paper on the rational model identification [19]. The total non-linear model emphasizes the non-linearity in both the parameters and control inputs, and it has been taken as a challenging structure for designing non-linear dynamic control systems [1]. The rational model gives more consideration as expanded polynomials in math, structure detection, and parameter estimation in the field of system identification [2,3]. Therefore, the main contribution of the new study is to use neural computing algorithms for a MIMO model parameter estimation. The new study is a complement to those classical NAMAX approaches.

The rest of the paper is organized as follows. The total non-linear model is described in Section 2. Section 3 presents the gradient descent calculation of parameter estimation. Next, model structure detection is discussed in Section 4. A convergence analysis of an algorithm is presented in Section 5. Simulation results and discussions are demonstrated in Section 6. Finally, Section 7 includes the paper conclusions and some of the future aspects.

2. Total Non-Linear Model

In mathematics, the dynamic total non-linear model of a MIMO system with error can be defined as

$${\mathrm{y}}_i(t)={\hat{\mathrm{y}}}_i+e_i(t)=\frac{a_i(t)}{b_i(t)}+e_i(t) =\frac{a_i(u_1,u_2,\dots ,u_J,y_1,y_2,\dots ,y_I,e_1,e_2,\dots ,e_I)}{b_i(u_1,u_2,\dots ,u_J,y_1,y_2,\dots ,y_I,e_1,e_2,\dots ,e_I)}+e_i(t)\mathrm{i}=1,2,\dots ,\mathrm{I}$$

(1)

$$a_i(t)={{\displaystyle \sum }}_{k=1}^N{p}_k^n(t){\mathsf{\theta }}_k^n,b_i(t)={{\displaystyle \sum }}_{k=1}^D{p}_k^d(t){\mathsf{\theta }}_k^d\mathrm{i}=1,2,\dots ,\mathrm{I}$$

(2)

where $y(t)=\left[y_1(t),y_2(t),\dots ,y_I(t)\right]\in R^I$and $\widehat{y}(t)=\left[{\hat{\mathrm{y}}}_1(t),{\hat{\mathrm{y}}}_2(t),\dots ,{\hat{\mathrm{y}}}_I(t)\right]\in R^I$are the measured output and model output, respectively; $u(t)=\left[u_1(t),u_2(t),\dots ,u_J(t)\right]\in R^J$ is the input; $e(t)=\left[e_1(t),e_2(t),\dots ,e_I(t)\right]\in R^I$ is the model error; and $t=\left[1,2,\dots ,\mathrm{T}\right]\in {\mathit{Z}}^{+T}$ is the sampling time index. Numerator$a_i(t)\in \mathrm{R}$ and denominator $b_i(t)\in \mathrm{R}$ as represented by polynomials, regression term $p_k^n(t)$, and $p_k^d(t)$ are products of past inputs, outputs, and errors, such as $u_1(t-1)y_2(t-3)$, $u_1(t-1)e_2(t-2)$, $y_2^3(t-1)$.${\mathsf{\theta }}^n=\left[{\mathsf{\theta }}_1^n,{\mathsf{\theta }}_2^n,\dots ,{\mathsf{\theta }}_N^n\right]\in R^N$, and ${\mathsf{\theta }}^d=\left[{\mathsf{\theta }}_1^d,{\mathsf{\theta }}_2^d,\dots ,{\mathsf{\theta }}_D^d\right]\in R^D$ are the parameter sets of $a_i(t)$ and $b_i(t)$, respectively.

The task of parameter estimation is to extract the relevant parameter values from the measured input and output data for a given model structure. To form a regression expression for parameter estimation, multiplying $e_i(t)$ of both sides of Formula (1) gives

$${\mathrm{y}}_i(t)b_i(t)-a_i(t)=b_i(t)e_i(t)$$

(3)

To consider the neuro-computing approach for parameter estimation, a total non-linear model is expressed into a non-completely connected feedforward neural network, as shown in Figure 1.

We define the network with an on both sides, Formula (11) is obtainedinput layer, a hidden layer, and an output layer, where:

(i)

The input layer consists of regression terms $p_k^n(t)(k=1,\dots ,N)$ and $p_k^d(t) (k=1,\dots ,D)$; here, a neuron in the hidden layer is not connected to all the neurons in the input layer, that is, the network is a non-completely connected feedforward neural network.

(ii)

The action function of the neurons in the hidden layer is linear, and the output of the hidden layer neurons is $a_i(t)$ or $b_i(t)$.

(iii)

The action function of the output layer neurons is linear, and the output of the ith output layer neuron is $b_i(t)e_i(t)$.

(iv)

The connection weights between the input layer neurons and the hidden layer neurons are the parameters ${\mathsf{\theta }}_k^n$ and ${\mathsf{\theta }}_k^d$ of the model.

(v)

The connection weight between the hidden layer neurons and the ith output layer neurons are $-1$ and the observed output ${\mathrm{y}}_i(t)$.

Leung and Haykin proposed a rational function neural network [20] but did not define a generalized total non-linear model structure or consider the relevant errors. Therefore, their parameter estimation algorithm could not provide an unbiased estimation for noise damaged data, which was essentially a special implementation of Zhu and Billings’s [7,8] methods in the case of no noise data. The method proposed in this paper is a further study of the method in Zhu [18]. The characteristics of a total non-linear model (1) are as follows:

(i)

By setting parameter $i=1$, Zhu’s [18] model can be a special case of the model in Formula (1).

(ii)

The model is non-linear in parameters and regression terms, which was caused by denominator polynomials.

(iii)

When the denominator $b_i(t)$ of the model is close to 0, the output deviation would be large. In this paper, considering this point, division operation was avoided in the action function of the neuron when the neural network model was being built.

(iv)

The structure of the neural network corresponding to the total non-linear model is a non-completely connected feedforward neural network, or a partially connected feedforward neural network. Therefore, the convergence of the network becomes a big problem, which is the difficulty of this paper.

(v)

The model has a wide range of application prospects. In many non-linear system modeling and control applications, the total non-linear model has been gradually adopted. Some non-linear models, such as the exponential model ${\mathrm{e}}^x$, which describes the change of dynamic rate constant with temperature, cannot be directly used. The exponential model can be firstly transformed into a non-linear model (${\mathrm{e}}^x=\frac{1-\frac{x}{2}+\frac{x^2}{12}}{1+\frac{x}{2}+\frac{x^2}{12}}$), and then, system identification can be implemented [19,21,22].

3. Gradient Descent Calculation of Parameter Estimation

For the convenience of the following derivations, set the output of neuron $i$ in the output layer of the neural network as $f_i(t)$.

$$f_i(t)=b_i(t)e_i(t)$$

(4)

Define the error measure function of one iteration of network as:

$$E(t)=\frac{1}{2}{{\displaystyle \sum }}^{\text{}}{(y_i(t)-\widehat{y_i}(t))}^2=\frac{1}{2}{{\displaystyle \sum }}^{\text{}}{(e_i(t))}^2$$

(5)

The Lyapunov method is often used to analyze the stability of a neural network [23]; similarly, the network parameters are estimated by minimizing the network error based on the Lyapunov method. It should be noted that when the total non-linear model is represented in the neural network structure of Figure 1, the parameter estimation of the model can be described as the training of neural network weight by minimizing the error $E(t)$ in Formula (5).

In order to train the weights of the network, the learning algorithm based on the gradient descent is given by Formulas (6) and (7):

$$\Delta {\mathsf{\theta }}_k^n=-{\eta }_n\frac{\partial \mathrm{E}}{\partial {\mathsf{\theta }}_k^n} =-{\eta }_n{e}_i(t)\frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^n}$$

(6)

$$\Delta {\mathsf{\theta }}_k^d=-{\eta }_{\mathrm{d}}\frac{\partial \mathrm{E}}{\partial {\mathsf{\theta }}_k^d} =-{\eta }_d{e}_i(t)\frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^d}$$

(7)

where ${\eta }_n$ and ${\eta }_{\mathrm{d}}$ are learning rates.

By deriving Formula (4) from ${\mathsf{\theta }}_k^n$ on both sides, Formula (8) is obtained:

$$\begin{gathered} \frac{\partial f_i(t)}{\partial {\mathsf{\theta }}_k^n}=\frac{\partial b_i(t)}{\partial {\mathsf{\theta }}_k^n}{e}_i(t) +b_i(t)\frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^n} \\ \frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^n} =\frac{1}{b_i(t)}(\frac{\partial f_i(t)}{\partial {\mathsf{\theta }}_k^n}-\frac{\partial b_i(t)}{\partial {\mathsf{\theta }}_k^n}{e}_i(t))=\frac{1}{b_i(t)}\frac{\partial f_i(t)}{\partial {\mathsf{\theta }}_k^n}=\frac{1}{b_i(t)}\frac{\partial (y_i(t)b_i(t)-a_i(t))}{\partial {\mathsf{\theta }}_k^n}=-\frac{1}{b_i(t)}\frac{\partial a_i(t)}{\partial {\mathsf{\theta }}_k^n}=-\frac{p_k^n(t)}{b_i(t)} \end{gathered}$$

(8)

Substituting Formula (8) into Formula (6) to get Formula (9), we can then get Formula (10):

$$\Delta {\mathsf{\theta }}_k^n=-{\eta }_n{e}_i(t)\frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^n}={\eta }_n{e}_i(t)\frac{p_k^n(t)}{b_i(t)}$$

(9)

$${\mathsf{\theta }}_k^n(\mathrm{t}+1)={\mathsf{\theta }}_k^n(\mathrm{t})+\Delta {\mathsf{\theta }}_k^n={\mathsf{\theta }}_k^n(\mathrm{t})+{\eta }_n{e}_i(t)\frac{p_k^n(t)}{b_i(t)}$$

(10)

By deriving Formula (4) from ${\mathsf{\theta }}_k^d$ on both sides, Formula (11) is obtained:

$$\begin{gathered} \frac{\partial f_i(t)}{\partial {\mathsf{\theta }}_k^d}=\frac{\partial b_i(t)}{\partial {\mathsf{\theta }}_k^d}{e}_i(t) +b_i(t)\frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^d} \\ \frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^d} =\frac{1}{b_i(t)}(\frac{\partial f_i(t)}{\partial {\mathsf{\theta }}_k^d}-\frac{\partial b_i(t)}{\partial {\mathsf{\theta }}_k^d}{e}_i(t)) \\ =\frac{1}{b_i(t)}(\frac{\partial (y_i(t)b_i(t)-a_i(t))}{\partial {\mathsf{\theta }}_k^d}-\frac{\partial b_i(t)}{\partial {\mathsf{\theta }}_k^d}{e}_i(t))=\frac{1}{b_i(t)}(y_i(t)p_k^d(t)-p_k^d(t)e_i(t)) \\ =\frac{1}{b_i(t)}(y_i(t)-e_i(t))p_k^d(t)=\frac{1}{b_i(t)} \frac{a_i(t)}{b_i(t)}{p}_k^d(t)=\frac{a_i(t)}{b_i^2(t)}{p}_k^d(t) \end{gathered}$$

(11)

Substituting Formula (11) into Formula (8) to get Formula (12), we then get Formula (13):

$$\Delta {\mathsf{\theta }}_k^d=-{\eta }_d{e}_i(t)\frac{\partial e_i(t)}{\partial {\mathsf{\theta }}_k^d}=-{\eta }_d{e}_i(t)\frac{a_i(t)}{b_i^2(t)}{p}_k^d(t)$$

(12)

$${\mathsf{\theta }}_k^d(t+1)={\mathsf{\theta }}_k^d(t)+\Delta {\mathsf{\theta }}_k^d ={\mathsf{\theta }}_k^d(t)-{\eta }_d{e}_i(t)\frac{a_i(t)}{b_i^2(t)}{p}_k^d(t)$$

(13)

The gradient descent algorithm for parameter estimation of a total non-linear model is summarized in Algorithm 1.

Algorithm 1. Gradient Descent Algorithm

1: Initialization: The weights of the neural network (parameters of a total non-linear model) are set as random little numbers with uniform distribution; the average value is zero, and the variance is small. Set the maximum number of iterations T, the minimum error ε, and the maximum number of samples $P$. 2: Generate training sample set {$X$,$Y$} of the neural network according to Formula (1), where $X=\left\{X_1,X_2,\dots ,X_I\right\}$, $Y$ = {$Y_1,Y_2,\dots ,Y_I$},    $X_i\ni${$p_1^{\mathrm{n}}(t)$,$p_2^{\mathrm{n}}(t)$,…,$p_N^{\mathrm{n}}(t), p_1^d$,$p_2^d,\dots ,p_D^d$}, ${\mathrm{Y}}_i$={$y_i(t)$}. 3: Input a training sample p to the neural network. 4: Calculate the output value $a_i(\mathrm{k})$,$y_i(t)e_i(t)$ and $f_i(t)$ of the neurons in the hidden layer and the output layer according to Formulas (2), (3), and (4), respectively. 5: Adjust the weight of the neural network according to Formulas (10) and (13). 6: Calculate the error $E(t)$ according to Formula (4) and calculate the total error according to Formula (14). $$E={{\displaystyle \sum }}^{\text{}}E(t)$$ (14) 7: $p=p+1$ 8: If p > P, then t = t + 1; otherwise, run step 3. 9: If $E<$ ε or $t>T$, stop training; otherwise, run step 3.

4. Model Structure Detection

Model structure detection is to select important items from a rather large model set (usually called the whole item set) and determine the sub-model with important items [18]. Because of the powerful self-learning and associative memory function of an artificial neural network [24], it is the first-choice tool to identify the model structure. When identifying systems with unknown structures, it is important to avoid losing these important items in the final model. For the structure detection of a total non-linear model, the connection weight estimation in the neural network, that is, the parameter estimation of the total non-linear model, could be used to select the significant terms.

For the important and unimportant items in the whole model item set, the knock-out algorithm is adopted. First, remove the items that lead to the increase of network error, and then knock out the items with lighter weight according to the requirements of significance. Finally, test the error of the non-linear model composed of the remaining items. The specific algorithm is summarized in Algorithm 2.

Algorithm 2. Knock-Out Algorithm

1: Using the network structure shown in Figure 1, all the items contained in the whole items set are taken as the input of the network. 2: The algorithm in Section 3 is used to train the network, and network error $E_1$ is obtained. 3: A new network structure is obtained by randomly removing a network input. The algorithm in Section 3 is used to train the new network, and network error $E_2$ is obtained. If $E_2\le E_1$, then $E_1=E_2$. Otherwise, this operation should be invalid (the input is reserved). 4: Another input is selected, and step 3 is executed again until all the input items are executed once. 5: The N connection weights between the input layer and the hidden layer are sorted in descending order. The first n weights are selected to make the significance reach 95%. Meanwhile, Formulas (15) and (16) are met, and the network input items corresponding to the first n weights are retained. $$\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^n\left|w_i\right|}{{{\displaystyle \sum }}_{i=1}^N\left|w_i\right|}\ge 0.95$$ (15) $$\frac{{{\displaystyle \sum }}_{i=1}^{n-1}\left|w_i\right|}{{{\displaystyle \sum }}_{i=1}^N\left|w_i\right|}<0.95$$ (16)

In the above process, the neural network is not only used to estimate the parameters of the model but also to detect the structure of the model and analyze the significance of the regression term.

5. Convergence Analysis of the Algorithm

Convergence proof:

Assuming that a connection weight of the neural network shown in Figure 1 is changed, this weight can take any value. When the weight ${\mathsf{\theta }}_k^n$corresponds to the regression term parameter of the numerator of the total non-linear model, the resulting network error changes as follows (remove the lower corner marks in Formula (2) for the convenience of proof):

$$\mathrm{y}(\mathrm{t})=\frac{a(t)}{b(t)}+e(t)$$

(17)

Substitute Formula (2) into Formula (1) to get Formula (18):

$$\mathrm{y}(\mathrm{t})-\frac{{{\displaystyle \sum }}_{\mathrm{n}=1}^N{p}_k^n(t){\mathsf{\theta }}_k^n}{b(t)}=e(t)$$

(18)

When ${\mathsf{\theta }}_k^n$ is updated, (18) becomes (19):

$$\mathrm{y}(\mathrm{t})-\frac{{{\displaystyle \sum }}_{\mathrm{n}=1,n\ne j}^N{p}_k^n(t){\mathsf{\theta }}_k^n+p_k^n(t)({\mathsf{\theta }}_k^n+\Delta {\mathsf{\theta }}_k^n)}{b(t)}=\tilde{e}(t)$$

(19)

$\tilde{e}(t)$ is the new error of the neural network after the weight has been updated. Subtract Formula (18) from Formula (19) to get Formulas (20) and (21):

$$\begin{gathered} \tilde{e}(t)-e(t)=-\frac{p_k^n(t)\Delta {\mathsf{\theta }}_k^n}{b(t)} =-{\eta }_n{(\frac{p_k^n(t)}{b(t)})}^2\mathrm{e}(t) \\ \tilde{e}(t) =(1-{\eta }_n(\frac{p_k^n(t)}{b(t)}{)}^2)\mathrm{e}(t) \end{gathered}$$

(20)

$$\tilde{e}{(t)}^2=(1-{\eta }_n(\frac{p_k^n(t)}{b(t)}{)}^2){}^2\mathrm{e}{(t)}^2$$

(21)

In order to ensure $\tilde{e}{(t)}^2\le \mathrm{e}{(t)}^2$, $-1\le 1-{\eta }_n{(\frac{p_k^n(t)}{b(t)})}^2\le 1$, namely:

$$\{\begin{array}{c}{\eta }_n{(\frac{p_k^n(t)}{b(t)})}^2\le 2\\ {\eta }_n{(\frac{p_k^n(t)}{b(t)})}^2\ge 0\end{array}$$

(22)

Solving Formula (22) gives:

$$0\le {\eta }_n\le \frac{2\mathrm{b}{(t)}^2}{p_k^n{(t)}^2}$$

(23)

When the changed weight ${\mathsf{\theta }}_k^d$ corresponds to the regression parameter of the denominator of the total non-linear model, the resulting network error change is as follows:

$$\mathrm{y}(\mathrm{t})-\frac{a(t)}{{{\displaystyle \sum }}_{\mathrm{d}=1}^D{p}_k^d(t){\mathsf{\theta }}_k^d}=e(t)$$

(24)

$$\mathrm{y}(\mathrm{t})-\frac{a(t)}{{{\displaystyle \sum }}_{\mathrm{d}=1,d\ne j}^D{p}_k^d(t)+p_k^d(t)({\mathsf{\theta }}_k^d+\Delta {\mathsf{\theta }}_k^d)}=\tilde{e}(t)$$

(25)

Subtracting Formula (24) from Formula (25) gives $\tilde{b}(t)$ as the new denominator of the neural network after the weight has been updated.

$$\begin{gathered} \tilde{e}(t)-e(t)=\frac{(\tilde{b}(t)-b(t))a(t)}{\tilde{b}(t)b(t)}=p_k^d(t)\Delta {\mathsf{\theta }}_k^d\frac{a(t)}{\tilde{b}(t)b(t)}=-{\eta }_{d}e(t)\frac{\partial e(t)}{\partial {\mathsf{\theta }}_k^d}{p}_k^d(t)\frac{a(t)}{\tilde{b}(t)b(t)} \\ =-{\eta }_{d}e(t)\frac{a{(t)}^2}{\tilde{b}(t)b{(t)}^3}{p}_k^d{(t)}^2 \end{gathered}$$

(26)

$$\tilde{e}{(t)}^{2 }={(1-{\eta }_d\frac{a{(t)}^2}{\tilde{b}(t)b{(t)}^3}{p}_k^d{(t)}^2)}^2 e{(t)}^2$$

(27)

In order to satisfy $\tilde{e}{(t)}^2\le \mathrm{e}{(t)}^2$, namely, $-1\le 1-{\eta }_d\frac{a{(t)}^2}{\tilde{b}(t)b{(t)}^3}{p}_k^d{(t)}^2\le 1$, that is:

$$\{\begin{array}{c}{\eta }_d\frac{a{(t)}^2}{\tilde{b}(t)b{(t)}^3}{p}_k^d{(t)}^2\le 2\\ {\eta }_d\frac{a{(t)}^2}{\tilde{b}(t)b{(t)}^3}{p}_k^d{(t)}^2\ge 0\end{array}$$

(28)

Because the learning coefficient is too large, the training effect of the network is not effective; accordingly, we take $0\le {\eta }_n\le 1$ to get $\tilde{b}(t)b(t)>0$, and thus, it has:

$$0\le {\eta }_d\le \frac{2\tilde{b}(t)b{(t)}^3}{a{(t)}^2{p}_k^d{(t)}^2}$$

(29)

To sum up, the network is convergent when the following conditions are met:

$$1.0\le {\eta }_n\le \frac{2\mathrm{b}{(t)}^2}{p_k^n{(t)}^2}$$

$$2.0\le {\eta }_d\le \frac{2\tilde{b}(t)b{(t)}^3}{a{(t)}^2{p}_k^d{(t)}^2}$$

Under these two conditions, this algorithm provides a convergence estimate for the parameters of the total non-linear model.

6. Simulation Results and Discussions

Consider a representative example of a total non-linear model:

$$\begin{gathered} y_1(t) =\frac{0.5y_1(t-1)+0.8y_2^3(t-2)+u_1(t-1)}{1+y_1^2(t-1)+u_2^2(t-1)}+{\mathrm{r}}_1(t) \\ \hspace{1em}\hspace{1em}\hspace{1em}=\frac{{\mathsf{\theta }}_1{y}_1(t-1)+{\mathsf{\theta }}_2{y}_2^3(t-2)+{\mathsf{\theta }}_3{u}_1(t-1)}{1+{\mathsf{\theta }}_4{y}_1^2(t-1)+{\mathsf{\theta }}_5{u}_2^2(t-1)}+{\mathrm{r}}_1(t) \end{gathered}$$

(30)

$$\begin{gathered} y_2(t) =\frac{0.2y_2(t-1)-0.5y_1^2(t-2)+u_2(t-1)}{1+y_2^2(t-1)+u_2^2(t-1)}+r_2(t) \\ \hspace{1em}\hspace{1em}\hspace{1em}=\frac{{\mathsf{\theta }}_6{y}_2(t-1)-{\mathsf{\theta }}_7{y}_1^2(t-2)+{\mathsf{\theta }}_8{u}_2(t-1)}{1+{\mathsf{\theta }}_9{y}_2^2(t-1)+{\mathsf{\theta }}_{10}{u}_2^2(t-1)}+r_2(t) \end{gathered}$$

(31)

Because the disturbance of input data will cause interference to the estimation of parameters [25], in this section, the parameter estimation for different inputs was selected. Firstly, for the simulation system without noise, 2000 pairs of input/output data were used as data sets for uniform sampling in 20 cycles, and the learning rate was designed as a linear attenuation sequence (in 50 iterations, the learning rate decreases from η₀ = 0.5 to η_end = 0.02). The algorithm in this paper was used to estimate 10 parameters at the same time.

Table 1. Parameter estimation of a noiseless system.

${\mathit{u}}_1(\mathit{t})$	${\mathit{u}}_2(\mathit{t})$	${\mathsf{\theta }}_1$	${\mathsf{\theta }}_2$	${\mathsf{\theta }}_3$	${\mathsf{\theta }}_4$	${\mathsf{\theta }}_5$	${\mathsf{\theta }}_6$	${\mathsf{\theta }}_7$	${\mathsf{\theta }}_8$	${\mathsf{\theta }}_9$	${\mathsf{\theta }}_{10}$	MSE
sine	sine	0.5002	0.8025	1.0003	1.0034	1.0000	0.2006	0.5010	1.0004	1.0018	0.9991	2.351E-06
sine	square	0.5000	0.8000	1.0000	1.0000	1.0000	0.1996	0.4982	1.0182	0.9677	1.0473	0.0003
square	square	0.4973	0.8760	1.0110	1.0031	1.0153	0.2013	0.5072	1.0354	0.9744	1.0840	0.0015

, after 50 iterations, shows the estimated values and mean square deviation of parameters.

The inputs $u_1(t)$ and $u_2(t)$ of the system are either a sine wave or square wave with an amplitude of 2. Figure 2 shows the difference between the measured value $y_1(t)$ (the real output value of Formula (6.1)) and the output value ${\hat{\mathrm{y}}}_1(t)$obtained using the parameter estimator when inputs $u_1(t)$ and $u_2(t)$ are both sine waves. In the same way, Figure 3 shows the difference between measured value $y_2(t)$ and output value ${\hat{\mathrm{y}}}_2(t)$ when inputs $u_1(t)$ and $u_2(t)$ are both sine waves. Figure 4 shows the difference between measured value $y_1(t)$ and output value ${\hat{\mathrm{y}}}_1(t)$ when input $u_1(t)$ is a sine wave and $u_2(t)$ is a square wave.

Figure 5 shows the difference between measured value $y_2(t)$ and output value ${\hat{\mathrm{y}}}_2(t)$ when input $u_1(t)$ is a sine wave and $u_2(t)$ is square wave. Figure 6 shows the difference between measured value $y_1(t)$ and output value ${\hat{\mathrm{y}}}_1(t)$ when input $u_1(t)$ and $u_2(t)$ are both square waves. Figure 7 shows the difference between measured value $y_2(t)$ and output value ${\hat{\mathrm{y}}}_2(t)$ when input $u_1(t)$ and $u_2(t)$ are both square waves. It can be seen that when the inputs are both sine waves, the accuracy of parameter estimation is the highest, while when the inputs are square waves, the estimation accuracy of the parameters is relatively lower. This is because when the inputs are square waves, the output of the system has an overshoot, that is to say, the observation value of the system itself has an error. Training the network with error data will certainly lead to an error of parameter estimation; especially the estimation error of ${\mathsf{\theta }}_2$ is the highest. Here, ${\mathsf{\theta }}_2$ is the parameter of $y_2^3(t-2)$ because $y_2^3(t-2)$ is the third power of the output, which further amplifies the error. Substituting the value of $y_2^3(t-2)$ into the update of ${\mathsf{\theta }}_2$ would inevitably lead to an estimation error.

For a system with noise interference, it is more difficult to estimate its parameters because of the error of the measured value itself [26,27]. In order to further verify the performance of the estimator, a system with noise was selected for parameter estimation, that is, by adding a noise signal to the original system. The noises ${\mathrm{r}}_1(t)$ and $r_2(t)$ were random uniform sequences, where each mean value was zero, each variance was 0.1, and each signal-to-noise ratio was 10. We repeated the previous experiment for the system with noise, and the experimental results are shown in

Table 2. Parameter estimation of a noisy system.

${\mathit{u}}_1(\mathit{t})$	${\mathit{u}}_2(\mathit{t})$	${\mathsf{\theta }}_1$	${\mathsf{\theta }}_2$	${\mathsf{\theta }}_3$	${\mathsf{\theta }}_4$	${\mathsf{\theta }}_5$	${\mathsf{\theta }}_6$	${\mathsf{\theta }}_7$	${\mathsf{\theta }}_8$	${\mathsf{\theta }}_9$	${\mathsf{\theta }}_{10}$	MSE
sine	sine	0.5003	0.8041	1.0005	1.0054	1.0001	0.2008	0.5014	1.0005	1.0016	0.9987	5.342E-06
sine	square	0.5000	0.8001	1.0000	1.0001	1.0000	0.2045	0.5019	1.073	1.1364	1.0898	0.0032
square	square	0.4953	0.8765	1.0085	1.0327	1.0095	0.2969	0.7030	0.9971	1.0007	0.9953	0.0058

. The difference between the measured value and the estimated value of the system output is shown in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13.

In order to detect the structure of the model, 20 items including 10 items in Formula (6.1) were used as the whole items set of models. The newly added numerator term is of order 1, the denominator term is of order 2, and the input lag and output lag are both of order 1. Using the knock-out algorithm in Section 4, the final 10 items of the model are in good agreement with those in Formula (6.1).

From the above experimental results, the estimation accuracy of the algorithm proposed in this paper is acceptable, and the mean square deviations are all less than 0.003. This level of error is acceptable.

7. Conclusions

In this paper, the parameter estimation of a SISO rational model was extended to that of a MIMO total non-linear model. A method of parameter estimation of a MIMO non-linear rational model based on a gradient descent algorithm was proposed, and the convergence condition was proposed for the asymmetry of the network. It was proven that the estimator is properly effective by mathematical derivation and simulation. This estimation method has a strong generalization property and could be widely used in many fields, such as non-linear system modeling and control applications. Some systems that could not directly use this method, such as the exponential model describing the change of the kinetic rate constant with the temperature, could first be converted into a rational model and then use the developed estimation method. Some of the future work could be foreseen as (1) estimating the parameters of the state space model based on an artificial neural network, (2) estimating the parameters of a MIMO state space model, (3) estimating the parameters of the non-linear state space model, and (4) estimating the parameters of total non-linear spatial state models.