# Identification of Wiener Model with Internal Noise Using a Cubic Spline Approximation-Bayesian Composite Quantile Regression Algorithm

Authors: Tianhong Pan; Wei Guo; Ying Song; Fujia Yin
Journal: Complexity
 (2020)

---

## Abstract

A cubic spline approximation-Bayesian composite quantile regression algorithm is proposed to estimate parameters and structure of the Wiener model with internal noise. Firstly, an ARX model with a high order is taken to represent the linear block; meanwhile, the nonlinear block (reversibility) is approximated by a cubic spline function. Then, parameters are estimated by using the Bayesian composite quantile regression algorithm. In order to reduce the computational burden, the Markov Chain Monte Carlo algorithm is introduced to calculate the expectation of parameters’ posterior distribution. To determine the structure order, the Final Output Error and the Akaike Information Criterion are used in the nonlinear block and the linear block, respectively. Finally, a numerical simulation and an industrial case verify the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

The Wiener model is a nonlinear system, which is composed of a dynamic linear block and a static nonlinear block. Due to its simple structure and strong adaptability, the Wiener model is widely used in industrial processes, such as chemical engineering processes [1], oil recovery [2], biological plants [3], and fluid control units [4]. However, how to identify the parameters and structure of the Wiener model is still a challenging issue.Most of the existing works are related to identify the Wiener model with a known structure. Xu et al. introduced a sigmoid function to improve the mutation rateF in the Adaptive Differential Evolution algorithm [5]. Furthermore, a blind identification method for the Wiener model was also investigated [6]. The input signal of the system was given by a cyclostationary signal instead of a Gaussian random signal and the internal variables were recovered only based on the output. Then, the order and parameters of the Wiener model were estimated by using the support vector machine regression.However, the structure and order of the Wiener model are usually unknown in the industrial process. It is difficult to identify the Wiener model because of its nonlinear characteristic. Additionally, the existing research studies are mostly related to the system without noise or only with the output noise [7–9]. Lamia et al. described the Wiener model using the polynomial nonlinear state space (PNLSS) model and developed an output error identification method for the nonlinear block [10]. Riccardo et al. presented a kernel-based identification to estimate parameters of the Wiener system [11]. The impulse response of the linear block and the static nonlinearity were modelled with a Gaussian process and combination of basis functions, respectively. Then, an iterative algorithm using the expectation-maximization method was developed to estimate parameters of the Wiener model. Al-Dhaifallah applied twin support vector regressions to identify the nonlinear Wiener system, including a linear dynamic block [12]. The linear block was expanded in terms of basis functions, while the nonlinear part is determined by twin support vector machine regressions.Although mentioned algorithms work well, the industrial process with the internal noise is less involved. Lindsten et al. handled the internal white noise in a state space systematic manner and used the nonparametric Gaussian process model for the static nonlinearity [13]. Jing et al. proposed a variable knot spline approximation recursive Bayesian algorithm to reduce the influence of the internal white noise [14]. Li et al. presented a neurofuzzy-based single-input-single-output (SISO) Wiener model identification method for colour noises [15]. Zhang and Mao proposed a robust recursive least squares algorithm with a dead zone weighted factor based on the inverse of the nonlinear function block, which took process noises and measurement noises into consideration [16]. Most of related articles are on the assumption that the internal process noise satisfies the Gaussian distribution or approximately satisfies the symmetrical distribution. Due to the internal noise change with the gain variation of the system, it is difficult to eliminate affections.(a)
Variables between the linear block and the nonlinear block cannot be measured, and the existing methods cannot estimate parameters accurately [9].(b)
Owing to the nonlinear characteristic of the Wiener model, the internal noise will be amplified with the gain of the nonlinear block [14].(c)
The inverse nonlinear block is commonly approximated by using the polynomial function. However, there are some oscillation phenomena in the polynomial with a high order [17]. As a result, the existing algorithms cannot achieve a satisfied estimation when the effect of external interpolation is not good.(d)
The best unbiased estimation of parameters can be obtained only if the modelling error satisfies the distribution with the mean zero and the same variance. Because the error cannot be guaranteed to satisfy the Gaussian distribution, the least squares algorithm no longer shows its robustness [18].To get an unbiased estimation of the Wiener model with the internal noise, a cubic spline approximation-Bayesian composite quantile regression (CSA-BCQR) is presented in this work. The contributions are as follows:(a)
The internal process noise with the unknown asymmetric distribution form in the Wiener model is addressed(b)
Repeated sampling by using the Markov Chain Monte Carlo (MCMC) method is implemented to achieve a faster convergence of parameters(c)
Overcome the high-order oscillation caused by using the polynomial approximation, and a structure detection framework is also consideredThe structure is organized as follows. Section2 gives the problem description. Section 3 states the principle of quantile regression and Bayesian composite quantile regression in detail. Section 4 describes the order determination of the nonlinear block and linear block of the Wiener model. Section 5 presents a numerical simulation and an industrial case to evaluate the proposed algorithm, respectively. Finally, conclusions are summarized in Section 6.


## 2. Problem Description

As shown in Figure1, the input signal uk firstly passes through the linear block of the Wiener model Bnq−1/Anq−1 and generates z1k. Under the disturbance of the internal noise εk, z1k is further developed into z2k and becomes the input of the nonlinear block fz2. The overall output yk is finally generated.Figure 1
The structure of the Wiener model with internal noise.Assume that(a)
The nonlinear functionfz2k is continuous, monotonic, and reversible(b)
εk∼N0,σv2 is considered as a Gaussian white noise and independent of the input signal uk(c)
Unknown polynomial order can be defined by the Final Output Error (FOE) and the Akaike Information Criterion (AIC)Here, an ARX model with a high order is employed to approximate the linear block of the Wiener model:(1)z2k=Bnq−1Anq−1uk+1Anq−1εk.Considering the nonlinear part is reversible, equation (1) is rewritten as follows:(2)Anq−1f−1yk=Bnq−1uk+εk.The cost function is defined as follows:(3)JARX=∑k=1NAnq−1f−1yk−Bnq−1uk2,where N is the size of samples. Because the specific form of the expression is unknown, a stable linear transfer function is approximated by the finite impulse response (FIR) model.(4)z1k=Bnq−1uk=∑j=1nbbjuk−j,where nb is the order of the linear block.For the nonlinear block of the Wiener model, a cubic spline approximation (CSA) function is applied to fit the inverse function of the nonlinear clock:(5)f−1yk=z2k=∑i=2nγ−1γiyk−yi3+γnγ+γnγ+1yk+γnγ+2y2k+γnγ+3y3k,where nγ is the order of spline function, γ2,γ3,⋯,γnγ+3 are the estimated parameters, y2,y3,⋯,ynγ−1 are the internal gathering points in the knot set y1,y2,⋯,ynγ and satisfy minyj=y1<y2<⋯ynγ−1<ynγ=maxyj, j=1,2,⋯,N, N is the size of samples, and γnγ+1=1.Becauseεk is unmeasurable, an arbitrary size gain is allocated between the linear module and the nonlinear module. Set Anq−1=1 and combine equations (2) and (5):(6)∑i=2nγ−1γiyk−yi3+γnγ+yk+γnγ+2y2k+γnγ+3y3k=Bnq−1uk+εk.Then, the Wiener model is parameterized as follows:(7)yk=φTkθ+εk,with(8)φk=−yk−y23,⋯,−yk−ynγ−13,−1,−y2k,−y3k,uk−1,⋯,uk−nbT∈ℜW×1,θ=γ2,⋯,γnγ−1,γnγ,γnγ+2,γnγ+3,b1,⋯,bnbT∈ℜW×1,W=nγ+nb+1.


## 3. Bayesian Composite Quantile Regression

### 3.1. Quantile Regression

According to the theory of the least squares method, the estimate is unbiased when the modelling errorε∼Norm0,σ2. If modelling error ε∼Normμ,σ2 and μ≠0, the estimate is unbiased with the minimum variance. When the data has a sharp/thick tail distribution and a significant heteroscedastic, the least squares algorithm no longer has the abovementioned superiority. To make up those deficiencies, a Quantile Regression (QR) has been proposed by Bang et al. [19].The nature of the QR is to describe the linear relationship between the independent variable and the dependent variable in the quantile formation. When the quantile changes from 0 to 1, the position and steering of the corresponding regression plane will be adjusted. Therefore, the dependent variable under different quantile conditions can be obtained. QR not only describes the range of variation of the dependent variable but also measures the influence of the conditional distribution shape of the regression variable. Unlike the least squares regression, QR has the best linear unbiased estimates and modelling residuals satisfy the normal distribution.If the conditional quantile ofyk is a linear function of φk, then Qτykφk=φTkθτ. The quantile estimate of θτ is shown as follows:(9)θ^τ=argminθ^∑k=1Nρτyk−φTkθ^τ,where ρτε=ετ−Iε<0,I=1,ε<00,others, ρτ is the linear loss function, and I⋅ is the indicator function.Equivalently,ρτε=τε,ε≥0τ−1ε,ε<0. Equation (9) can also be expressed as follows:(10)θ^τ=argminθ^∑k:yk≥θ^τφTkτyk−φTkθ^τ+∑k:yk<θ^τφTk1−τyk−φTkθ^τ.Whenτ takes different values, the estimated parameter is also different, as well as the expression of the regression equation.However, the objective function in QR is nondifferentiable or nonconvex, which makes the optimal solution process complicate [20]. Both equations (9) and (10) are piecewise linear functions, so the solution cannot be solved directly from the objective function. It is necessary to find another way to obtain the posterior density function of each QR’s parameter.Definition 1.
A random variablex has an Asymmetric Laplace Distribution (ALD) which is noted as ALDλ,m with m∈ℜ and λ>0 [21]. Its probability density function (pdf) is given by(11)pdfxm,λ,κ=λκ+1/κexp−λκx−m,x≥m,expλκx−m,x<m,where m is a location parameter, λ>0 is a scale parameter, and κ is an asymmetry parameter. When κ=1, the distribution simplifies to the Laplace distribution.
Assuming thatεk satisfies the ALD, the corresponding pdf is written as follows [22]:(12)pdfε=τ1−τexp−ρτε,and the likelihood function of θ is defined as follows:(13)Lyθ=τN1−τNexp−∑k=1Nρτyk−φTkθτ.
Then, minimization of the cost function equation (10) is equivalent to the maximization of the likelihood function equation (13) based on ALD [23].


### 3.2. Bayesian Composite Quantile Regression

Compared with least squares algorithm, QR has much more flexibility in assessing the effect of predictors on different locations of the response distribution [22]. It also solves the problem that the least squares regression can only describe the local influence of the dependent variable to the dependent variable. Zou and Yuan pointed out that QR can lead to arbitrarily small relative efficiency compared with the least square [24]. In addition, the median regression, a special case of QR, may not be the best choice for some abnormal errors. Therefore, the Composite Quantile Regression (CQR) can aggregate multiple quantile information together to produce a robust parameter estimation. Meanwhile, CQR can achieve an estimated efficiency gain based on a single quantile regression [25].Withτ=τ1,τ2,⋯,τl, 0<τ1<τ2<⋯<τl<1, the CQR of θτ is(14)θ^τ=argminθ∑j=1l∑k=1Nρτjyk−φTkθτ,which is also rewritten by(15)θ^τ=argmin∑j=1l∑k:yk≥θτφTkτjyk−φTkθτ+∑k:yk<θτφTk1−τjyk−θτφTk.Obviously, the estimatedθ can be calculated by the QR of a combination of τj. However, the computational burdens of CQR and QR are high. It is necessary to find a good solution.The posterior distribution ofθ, πθy is given by(16)πθy∝Lyθpθ,where pθ is the prior distribution of θ, Lyθ is the likelihood function:(17)Lyθ=τN1−τNexp−∑j=1l∑k=1Nρτjyk−φTkθτ.Lemma 1.
If the likelihood function isLyθ=τN1−τNexp−∑j=1l∑k=1Nρτjyk−φTkθτ, then the posterior distribution of θ, πθy, will have a proper distribution [25], i.e.,(18)0<∫πθydθ=∫Lyθpθdθ<∞.Ifpθ=1, −∞<θ<+∞, πθy still exists and the expectation of the posterior distribution is the Bayesian estimate of θ.In order to improve the computational efficiency, the Markov Chain Monte Carlo (MCMC) method, which is convenient for researchers to sample from complex distributions, has been developed [26]. There are two main MCMC sampling methods: Metropolis-Hastings (M-H) sampling and Gibbs sampling. Here, the M-H algorithm is used. The procedure of M-H sampling is shown in Algorithm 1.The Bayesian CQR (BCQR) algorithm used in this paper does not depend on the actual distribution of data, but on the likelihood function formed by the ALD [27]. The essence of BCQR is that the estimated parameter is regarded as a random variable, and the sampling distribution of parameter can be obtained by repeated sampling. When a stable distribution is obtained, the mean and standard deviation of the parameters at each quantile can be determined.Algorithm 1:Pseudocode of M-H sampling algorithm.
Begin
Define the maximum number of state transitionsntrans, the number of samples nsamp, the input state transition probability of  Markov Chain qtrans, input smooth distribution πx, initialize variable j and state xj=0.
while j≤ntrans+nsamp−1 do
calculateα=πx′qtransx′,x/πxqtransx,x′;
determine the probability of receptionαacceptx,x′=min1,α;
uni=unifrnd0,1;
if uni<αacceptx,x′
xj+1=x′;
else
xj+1=xj;
end if
j=j+1;
end while
end begin


## 3.1. Quantile Regression

According to the theory of the least squares method, the estimate is unbiased when the modelling errorε∼Norm0,σ2. If modelling error ε∼Normμ,σ2 and μ≠0, the estimate is unbiased with the minimum variance. When the data has a sharp/thick tail distribution and a significant heteroscedastic, the least squares algorithm no longer has the abovementioned superiority. To make up those deficiencies, a Quantile Regression (QR) has been proposed by Bang et al. [19].The nature of the QR is to describe the linear relationship between the independent variable and the dependent variable in the quantile formation. When the quantile changes from 0 to 1, the position and steering of the corresponding regression plane will be adjusted. Therefore, the dependent variable under different quantile conditions can be obtained. QR not only describes the range of variation of the dependent variable but also measures the influence of the conditional distribution shape of the regression variable. Unlike the least squares regression, QR has the best linear unbiased estimates and modelling residuals satisfy the normal distribution.If the conditional quantile ofyk is a linear function of φk, then Qτykφk=φTkθτ. The quantile estimate of θτ is shown as follows:(9)θ^τ=argminθ^∑k=1Nρτyk−φTkθ^τ,where ρτε=ετ−Iε<0,I=1,ε<00,others, ρτ is the linear loss function, and I⋅ is the indicator function.Equivalently,ρτε=τε,ε≥0τ−1ε,ε<0. Equation (9) can also be expressed as follows:(10)θ^τ=argminθ^∑k:yk≥θ^τφTkτyk−φTkθ^τ+∑k:yk<θ^τφTk1−τyk−φTkθ^τ.Whenτ takes different values, the estimated parameter is also different, as well as the expression of the regression equation.However, the objective function in QR is nondifferentiable or nonconvex, which makes the optimal solution process complicate [20]. Both equations (9) and (10) are piecewise linear functions, so the solution cannot be solved directly from the objective function. It is necessary to find another way to obtain the posterior density function of each QR’s parameter.Definition 1.
A random variablex has an Asymmetric Laplace Distribution (ALD) which is noted as ALDλ,m with m∈ℜ and λ>0 [21]. Its probability density function (pdf) is given by(11)pdfxm,λ,κ=λκ+1/κexp−λκx−m,x≥m,expλκx−m,x<m,where m is a location parameter, λ>0 is a scale parameter, and κ is an asymmetry parameter. When κ=1, the distribution simplifies to the Laplace distribution.
Assuming thatεk satisfies the ALD, the corresponding pdf is written as follows [22]:(12)pdfε=τ1−τexp−ρτε,and the likelihood function of θ is defined as follows:(13)Lyθ=τN1−τNexp−∑k=1Nρτyk−φTkθτ.
Then, minimization of the cost function equation (10) is equivalent to the maximization of the likelihood function equation (13) based on ALD [23].


## 3.2. Bayesian Composite Quantile Regression

Compared with least squares algorithm, QR has much more flexibility in assessing the effect of predictors on different locations of the response distribution [22]. It also solves the problem that the least squares regression can only describe the local influence of the dependent variable to the dependent variable. Zou and Yuan pointed out that QR can lead to arbitrarily small relative efficiency compared with the least square [24]. In addition, the median regression, a special case of QR, may not be the best choice for some abnormal errors. Therefore, the Composite Quantile Regression (CQR) can aggregate multiple quantile information together to produce a robust parameter estimation. Meanwhile, CQR can achieve an estimated efficiency gain based on a single quantile regression [25].Withτ=τ1,τ2,⋯,τl, 0<τ1<τ2<⋯<τl<1, the CQR of θτ is(14)θ^τ=argminθ∑j=1l∑k=1Nρτjyk−φTkθτ,which is also rewritten by(15)θ^τ=argmin∑j=1l∑k:yk≥θτφTkτjyk−φTkθτ+∑k:yk<θτφTk1−τjyk−θτφTk.Obviously, the estimatedθ can be calculated by the QR of a combination of τj. However, the computational burdens of CQR and QR are high. It is necessary to find a good solution.The posterior distribution ofθ, πθy is given by(16)πθy∝Lyθpθ,where pθ is the prior distribution of θ, Lyθ is the likelihood function:(17)Lyθ=τN1−τNexp−∑j=1l∑k=1Nρτjyk−φTkθτ.Lemma 1.
If the likelihood function isLyθ=τN1−τNexp−∑j=1l∑k=1Nρτjyk−φTkθτ, then the posterior distribution of θ, πθy, will have a proper distribution [25], i.e.,(18)0<∫πθydθ=∫Lyθpθdθ<∞.Ifpθ=1, −∞<θ<+∞, πθy still exists and the expectation of the posterior distribution is the Bayesian estimate of θ.In order to improve the computational efficiency, the Markov Chain Monte Carlo (MCMC) method, which is convenient for researchers to sample from complex distributions, has been developed [26]. There are two main MCMC sampling methods: Metropolis-Hastings (M-H) sampling and Gibbs sampling. Here, the M-H algorithm is used. The procedure of M-H sampling is shown in Algorithm 1.The Bayesian CQR (BCQR) algorithm used in this paper does not depend on the actual distribution of data, but on the likelihood function formed by the ALD [27]. The essence of BCQR is that the estimated parameter is regarded as a random variable, and the sampling distribution of parameter can be obtained by repeated sampling. When a stable distribution is obtained, the mean and standard deviation of the parameters at each quantile can be determined.Algorithm 1:Pseudocode of M-H sampling algorithm.
Begin
Define the maximum number of state transitionsntrans, the number of samples nsamp, the input state transition probability of  Markov Chain qtrans, input smooth distribution πx, initialize variable j and state xj=0.
while j≤ntrans+nsamp−1 do
calculateα=πx′qtransx′,x/πxqtransx,x′;
determine the probability of receptionαacceptx,x′=min1,α;
uni=unifrnd0,1;
if uni<αacceptx,x′
xj+1=x′;
else
xj+1=xj;
end if
j=j+1;
end while
end begin


## 4. Order Selection of the Wiener Model

The FOE and AIC are commonly used to determine the order of the nonlinear block and linear block, respectively.


### 4.1. Order Determination of the Nonlinear Block

As shown in equation (5), the order of the inverse nonlinear function is nγ and the corresponding FOE is as follows:(19)FOEynγ=N+WN−W×1N∑k=1Nyk−f^B^q−1uk2,where N is the size of samples and W is the number of parameters. N+W/N−W is the factor which overcomes the overfitting issue. However, equation (19) requires the estimated value f^z2k. Therefore, an equivalent criterion is introduced:(20)FOEz2nγ=N+WN−W×1N∑k=1Nf^−1yk−B^q−1uk2,where FOEynγ=FOEz2nγ.Assume(21)f^f^−1yk=yk.Since the nonlinear function is monotonic, thenminFOEz2nγ=minFOEynγ.


### 4.2. Order Determination of the Linear Block

FOE criterion needs to know the input of the system and the output of the linear block, but the intermediate variables of the system are unmeasurable. So, an indirect method is adopted to estimate the intermediate variables firstly and the order is determined lately.According to the estimated parameterθ^, the linear module is(22)z^1k=∑j=1nbbjuk−j,with(23)θnb=b1,b2,⋯,bnbT,z^nb=z^11,z^12,⋯,z^1kT,Ηnb=u0u1⋯u1+nbu1u2⋯u2+nb⋮⋮⋮uN−1uN⋯uN+nbT.Equation (22) can also be rewritten by(24)z^nb=ΗnbTθnb.The likelihood function ofθnb is Lθnb. The maximum likelihood estimation of θnb is given by(25)θ^ML=ΗnbTΗnb−1ΗnbTz^nb.Using AIC criterion,(26)AICn^b=−2logLθ^ML+2n^b.Finally, the order is determined withminAICn^b.


## 4.1. Order Determination of the Nonlinear Block

As shown in equation (5), the order of the inverse nonlinear function is nγ and the corresponding FOE is as follows:(19)FOEynγ=N+WN−W×1N∑k=1Nyk−f^B^q−1uk2,where N is the size of samples and W is the number of parameters. N+W/N−W is the factor which overcomes the overfitting issue. However, equation (19) requires the estimated value f^z2k. Therefore, an equivalent criterion is introduced:(20)FOEz2nγ=N+WN−W×1N∑k=1Nf^−1yk−B^q−1uk2,where FOEynγ=FOEz2nγ.Assume(21)f^f^−1yk=yk.Since the nonlinear function is monotonic, thenminFOEz2nγ=minFOEynγ.


## 4.2. Order Determination of the Linear Block

FOE criterion needs to know the input of the system and the output of the linear block, but the intermediate variables of the system are unmeasurable. So, an indirect method is adopted to estimate the intermediate variables firstly and the order is determined lately.According to the estimated parameterθ^, the linear module is(22)z^1k=∑j=1nbbjuk−j,with(23)θnb=b1,b2,⋯,bnbT,z^nb=z^11,z^12,⋯,z^1kT,Ηnb=u0u1⋯u1+nbu1u2⋯u2+nb⋮⋮⋮uN−1uN⋯uN+nbT.Equation (22) can also be rewritten by(24)z^nb=ΗnbTθnb.The likelihood function ofθnb is Lθnb. The maximum likelihood estimation of θnb is given by(25)θ^ML=ΗnbTΗnb−1ΗnbTz^nb.Using AIC criterion,(26)AICn^b=−2logLθ^ML+2n^b.Finally, the order is determined withminAICn^b.


## 5. Case Studies

### 5.1. Numerical Simulation

The discrete Wiener model with an internal noise is given by(27)yk=fBnq−1uk+εk,fz2k=z2k+3z23k,where Bnq−1=1−q−1+3q−2−2q−3, uk=unifrnd0,0.4, εk satisfies the ALD, and N=600.


#### 5.1.1. Estimation ofnγ

The order of the Cubic Spline Function is estimated by using the FOE criterion under different variances. The variation of the FOE of the cubic spline function with different quantiles was shown in Figure2. When nγ=6, FOE is minimum.Figure 2
Order of the nonlinear block using FOE.


#### 5.1.2. Estimation ofnb

As shown in Figure3, AIC is minimum when n^b=3, which is completely consistent with the true order of the Wiener model.Figure 3
Order of the linear block using AIC.


#### 5.1.3. Parameter Estimation Using CSA-BCQR

To estimate the error of the linear part and nonlinear part, absolute relative error (ARE) and mean square error (MSE) are defined as follows:(28)ARE=b^1,b^2,b^3−b1,b2,b3b1,b2,b3×100%,MSE=1N∑k=1Nyk−y^k2,and the comprehensive error (CE) is(29)CE=nb×AREmeanARE+nγ−1×MSEmeanMSE.The estimated results were shown in Table1 and Figure 4. To show the superiority of CAS-BCQR algorithm, CSA-Bayesian quantile regression (CSA-BQR) and CSA-recursive least squares (CSA-RLS) algorithms were taken for comparison.Table 1
Comparison among CSA-BCQR, CSA-BQR, and CSA-RLS in the numerical simulation.
Variable
CSA-BCQR
CSA-BQR
CSA-RLS
Linear block
b
1
−1.001
−1.01
−0.996
b
2
2.99
2.95
2.93
b
3
−2.0
−2.07
−2.11
ARE
2.69E − 03
2.31E − 02
3.49E − 02
Nonlinear block
MSE
3.121E − 04
5.675E − 04
9.554E − 04
Comprehensive error
CE
2.95
8.06
12.98Figure 4
Parameter identification results of the linear block.When using CSA-BCQR algorithm,b1,b2,andb3 converge at the iteration of 195, 79, and 83, respectively. However, in CSA-BQR algorithm, b1 does not converge when the iteration time is 300 and b2andb3 converge at 91 and 249, respectively. In CSA-RLS algorithm, b1 also fails to converge at 300 and b2andb3 converges at 190 and 298, respectively. Obviously, the convergence speed of CSA-BCQR algorithm is faster than those of the CSA-BQR and CSA-RLS algorithms.


### 5.2. Industrial Case

A water tank model was taken for an example (shown in Figure5). uk is the water inlet flow, z1k is the theoretical liquid level, z2k is the actual liquid level, and yk is the water outlet flow [28].Figure 5
A water tank model.The relationship between the variables of the model is as follows:(30)z1k=0.8z−1+0.4z−2+0.2z−3+0.1z−4uk,z2k=z1k+εk,yk=z22k+3z2k,where uk∼Norm0,12 is a persistent excitation signal sequence, N=500. Here, nb=4 and nγ=5 and parameters estimated by CSA-BCQR algorithm were shown in Table 2.Table 2
Comparison among CSA-BCQR, CSA-BQR, and CSA-RLS in the water tank.
Algorithm
k
b
1
b
2
b
3
b
4
ARE
True value
0.8
0.4
0.2
0.1
CSA-BCQR
1
0.9543
0.5544
0.3645
0.2459
0.3360
50
0.8953
0.4929
0.3208
0.1912
0.2186
100
0.8589
0.4746
0.2678
0.1589
0.1418
150
0.8161
0.4098
0.2087
0.1157
0.0282
200
0.8007
0.4014
0.2057
0.1009
0.0065
CSA-BQR
1
0.9735
0.5968
0.3867
0.2644
0.3922
50
0.9168
0.5115
0.3342
0.2319
0.2689
100
0.8794
0.4896
0.2794
0.2046
0.1927
150
0.8248
0.4687
0.2301
0.1687
0.1136
200
0.8091
0.4284
0.2126
0.1166
0.0395
CSA-RLS
1
0.9859
0.5888
0.3911
0.2721
0.4005
50
0.9277
0.5472
0.3558
0.2526
0.3172
100
0.9024
0.5078
0.2942
0.2311
0.2380
150
0.8585
0.4995
0.2489
0.1848
0.1642
200
0.8159
0.4232
0.2151
0.1205
0.0411As shown in Table2, the estimation errors of the CSA-BCQR, CSA-BQR, and CSA-RLS algorithms decrease gradually with k increasing. Meanwhile the estimated value was getting closer to the true value. When k=150, the ARE of CSA-BCQR algorithm was 0.0282 and the ARE of CSA-BQR and CSA-RLS were 0.1136 and 0.1642, respectively.As shown in Figure6, the MSE of CSA-BCQR algorithm has a significant reduction and converges after 100 iterations. Using CSA-BQR and CSA-RLS, the accuracy of the nonlinear block parameter identification is less than that of CSA-BCQR algorithm. It can be seen that the CSA-BCQR algorithm has higher accuracy and faster convergence.Figure 6
MSE of the nonlinear block.


## 5.1. Numerical Simulation

The discrete Wiener model with an internal noise is given by(27)yk=fBnq−1uk+εk,fz2k=z2k+3z23k,where Bnq−1=1−q−1+3q−2−2q−3, uk=unifrnd0,0.4, εk satisfies the ALD, and N=600.


### 5.1.1. Estimation ofnγ

The order of the Cubic Spline Function is estimated by using the FOE criterion under different variances. The variation of the FOE of the cubic spline function with different quantiles was shown in Figure2. When nγ=6, FOE is minimum.Figure 2
Order of the nonlinear block using FOE.


### 5.1.2. Estimation ofnb

As shown in Figure3, AIC is minimum when n^b=3, which is completely consistent with the true order of the Wiener model.Figure 3
Order of the linear block using AIC.


### 5.1.3. Parameter Estimation Using CSA-BCQR

To estimate the error of the linear part and nonlinear part, absolute relative error (ARE) and mean square error (MSE) are defined as follows:(28)ARE=b^1,b^2,b^3−b1,b2,b3b1,b2,b3×100%,MSE=1N∑k=1Nyk−y^k2,and the comprehensive error (CE) is(29)CE=nb×AREmeanARE+nγ−1×MSEmeanMSE.The estimated results were shown in Table1 and Figure 4. To show the superiority of CAS-BCQR algorithm, CSA-Bayesian quantile regression (CSA-BQR) and CSA-recursive least squares (CSA-RLS) algorithms were taken for comparison.Table 1
Comparison among CSA-BCQR, CSA-BQR, and CSA-RLS in the numerical simulation.
Variable
CSA-BCQR
CSA-BQR
CSA-RLS
Linear block
b
1
−1.001
−1.01
−0.996
b
2
2.99
2.95
2.93
b
3
−2.0
−2.07
−2.11
ARE
2.69E − 03
2.31E − 02
3.49E − 02
Nonlinear block
MSE
3.121E − 04
5.675E − 04
9.554E − 04
Comprehensive error
CE
2.95
8.06
12.98Figure 4
Parameter identification results of the linear block.When using CSA-BCQR algorithm,b1,b2,andb3 converge at the iteration of 195, 79, and 83, respectively. However, in CSA-BQR algorithm, b1 does not converge when the iteration time is 300 and b2andb3 converge at 91 and 249, respectively. In CSA-RLS algorithm, b1 also fails to converge at 300 and b2andb3 converges at 190 and 298, respectively. Obviously, the convergence speed of CSA-BCQR algorithm is faster than those of the CSA-BQR and CSA-RLS algorithms.


## 5.1.1. Estimation ofnγ

The order of the Cubic Spline Function is estimated by using the FOE criterion under different variances. The variation of the FOE of the cubic spline function with different quantiles was shown in Figure2. When nγ=6, FOE is minimum.Figure 2
Order of the nonlinear block using FOE.


## 5.1.2. Estimation ofnb

As shown in Figure3, AIC is minimum when n^b=3, which is completely consistent with the true order of the Wiener model.Figure 3
Order of the linear block using AIC.


## 5.1.3. Parameter Estimation Using CSA-BCQR

To estimate the error of the linear part and nonlinear part, absolute relative error (ARE) and mean square error (MSE) are defined as follows:(28)ARE=b^1,b^2,b^3−b1,b2,b3b1,b2,b3×100%,MSE=1N∑k=1Nyk−y^k2,and the comprehensive error (CE) is(29)CE=nb×AREmeanARE+nγ−1×MSEmeanMSE.The estimated results were shown in Table1 and Figure 4. To show the superiority of CAS-BCQR algorithm, CSA-Bayesian quantile regression (CSA-BQR) and CSA-recursive least squares (CSA-RLS) algorithms were taken for comparison.Table 1
Comparison among CSA-BCQR, CSA-BQR, and CSA-RLS in the numerical simulation.
Variable
CSA-BCQR
CSA-BQR
CSA-RLS
Linear block
b
1
−1.001
−1.01
−0.996
b
2
2.99
2.95
2.93
b
3
−2.0
−2.07
−2.11
ARE
2.69E − 03
2.31E − 02
3.49E − 02
Nonlinear block
MSE
3.121E − 04
5.675E − 04
9.554E − 04
Comprehensive error
CE
2.95
8.06
12.98Figure 4
Parameter identification results of the linear block.When using CSA-BCQR algorithm,b1,b2,andb3 converge at the iteration of 195, 79, and 83, respectively. However, in CSA-BQR algorithm, b1 does not converge when the iteration time is 300 and b2andb3 converge at 91 and 249, respectively. In CSA-RLS algorithm, b1 also fails to converge at 300 and b2andb3 converges at 190 and 298, respectively. Obviously, the convergence speed of CSA-BCQR algorithm is faster than those of the CSA-BQR and CSA-RLS algorithms.


## 5.2. Industrial Case

A water tank model was taken for an example (shown in Figure5). uk is the water inlet flow, z1k is the theoretical liquid level, z2k is the actual liquid level, and yk is the water outlet flow [28].Figure 5
A water tank model.The relationship between the variables of the model is as follows:(30)z1k=0.8z−1+0.4z−2+0.2z−3+0.1z−4uk,z2k=z1k+εk,yk=z22k+3z2k,where uk∼Norm0,12 is a persistent excitation signal sequence, N=500. Here, nb=4 and nγ=5 and parameters estimated by CSA-BCQR algorithm were shown in Table 2.Table 2
Comparison among CSA-BCQR, CSA-BQR, and CSA-RLS in the water tank.
Algorithm
k
b
1
b
2
b
3
b
4
ARE
True value
0.8
0.4
0.2
0.1
CSA-BCQR
1
0.9543
0.5544
0.3645
0.2459
0.3360
50
0.8953
0.4929
0.3208
0.1912
0.2186
100
0.8589
0.4746
0.2678
0.1589
0.1418
150
0.8161
0.4098
0.2087
0.1157
0.0282
200
0.8007
0.4014
0.2057
0.1009
0.0065
CSA-BQR
1
0.9735
0.5968
0.3867
0.2644
0.3922
50
0.9168
0.5115
0.3342
0.2319
0.2689
100
0.8794
0.4896
0.2794
0.2046
0.1927
150
0.8248
0.4687
0.2301
0.1687
0.1136
200
0.8091
0.4284
0.2126
0.1166
0.0395
CSA-RLS
1
0.9859
0.5888
0.3911
0.2721
0.4005
50
0.9277
0.5472
0.3558
0.2526
0.3172
100
0.9024
0.5078
0.2942
0.2311
0.2380
150
0.8585
0.4995
0.2489
0.1848
0.1642
200
0.8159
0.4232
0.2151
0.1205
0.0411As shown in Table2, the estimation errors of the CSA-BCQR, CSA-BQR, and CSA-RLS algorithms decrease gradually with k increasing. Meanwhile the estimated value was getting closer to the true value. When k=150, the ARE of CSA-BCQR algorithm was 0.0282 and the ARE of CSA-BQR and CSA-RLS were 0.1136 and 0.1642, respectively.As shown in Figure6, the MSE of CSA-BCQR algorithm has a significant reduction and converges after 100 iterations. Using CSA-BQR and CSA-RLS, the accuracy of the nonlinear block parameter identification is less than that of CSA-BCQR algorithm. It can be seen that the CSA-BCQR algorithm has higher accuracy and faster convergence.Figure 6
MSE of the nonlinear block.


## 6. Conclusion

A Cubic Spline Approximation-Bayesian Composite Quantile Regression algorithm is presented to estimate the structure and parameter of the Wiener system with internal noise. Using a cubic spline function to approximate the nonlinear block, overcome the high-order oscillation caused by the polynomial approximation. Then, using Bayesian Composite Quantile Regression algorithm to consider different information of quantiles can effectively improve the accuracy of parameter estimation. A numerical simulation and an industrial case show that the CSA-BCQR algorithm has faster convergence speed and higher parameter identification accuracy compared with CSA-BQR and CSA-RLS algorithms. Furthermore, the CSA-BCQR algorithm may also be applied to other block-oriented models, such as a Hammerstein system or a Hammerstein-Wiener system.



---

*Source: 9195819-2020-01-09.xml*