Robust LPV models identification approach based on shifted asymmetric Laplace distribution

Abstract

This paper focuses on the robust parameters estimation algorithm of linear parameters varying (LPV) models. The classical robust identification techniques deal with the polluted training data, for example, outliers in white noise. The paper extends this robustness to both symmetric and asymmetric noise with outliers to achieve stronger robustness. Without the assumption of Gaussian white noise pollution, the paper employs asymmetric Laplace distribution to model broader noise, especially the asymmetrically distributed noise, since it is an asymmetric heavy-tailed distribution. Furthermore, the asymmetric Laplace (AL) distribution is represented as the product of Gaussian distribution and exponential distribution to decompose this complex AL distribution. Then, a shifted parameter is introduced as the regression term to connect the probabilistic models of the noise and the predict output that obeys shifted AL distribution. In this way, the posterior probability distribution of the unobserved variables could be deduced and the robust parameters estimation problem is solved in the general Expectation Maximization algorithm framework. To demonstrate the advantage of the proposed algorithm, a numerical simulation example is employed to identify the parameters of LPV models and to illustrate the convergence.

Introduction

The modern industry is developing fast, and systems of various industrial processes become larger and the internal relations become more complex. ¹ Establishing accurate processes’ models is the foundation of control, fault diagnosis and condition monitoring, and therefore higher requirements are put forward for the precision of modeling.^2–4 As a result, many existing researches focus on the complex process modeling techniques with complex structure, strong dynamic characteristics and strong nonlinearity in modern industrial process control. ⁵ Linear parameter varying (LPV) systems attract much attention since this has a simple linear structure and variable parameters, and therefore can model complex and time-varying processes. In addition, based on LPV systems, many controllers have been studied to achieve stable and robust control.^6,7

Because of its simple structure, the least square estimation method can be employed to identify the parameters of the most classical LPV models. However, this requires high quality training data, which is almost impossible in practice. As for the state space LPV models, the various parameters identification problems have been studied. For instance, ⁸ proposed the parameters estimation algorithm to identify state space LPV models and the identification approach of multiple-input single output systems with time-delay was proposed in Yang et al. ⁹ In addition, to reduce computational cost, Turk, Gillis, Pipeleers & Swevers proposed the regularization method, and used $l_{2,1}$ norm regularization to identify the LPV systems.

As for the single input and single output systems, many research focuses on the robust identification problems and addressing various challenges posed by the training data.^10–12 For example, sensors hardly collect all the information roundly and accurately. Therefore, the algorithm that can estimate the parameters when output measurements are randomly missing is proposed in Yang et al. ¹³ and Liu et al. ¹⁴ Output-error (OE) structure is employed in Yang et al. ¹⁵ as the linear structure of LPV models that has the advantages and flexibility to model various noise, and proposed a robust parameters estimation method to identify the parameters when the training data is polluted by outliers and the system has time delays. Moreover, when it comes to the LPV piecewise systems,a membership recognition method is proposed in Hamdi et al. ¹⁶ In practice, because of the different working environments, different types of sensors are usually employed and therefore the sampling frequencies may be also different. To solve this, the dual-rate parameters estimation approach is proposed in Yang and Yin ¹⁷ in which the LPV systems can be identified under the expectation maximization framework. Among these researches, Bayesian approach is of significance, and therefore the probabilistic models definitely determine the identification ability. For example, compared with Gaussian distribution, Laplace distribution and Student-T distribution are employed to model the noise since they are heavy-tailed distributions, and therefore the system identification algorithms developed based on these are robust to the outliers. ¹⁸ However, heavy-tailed distributions are not perfect since they can only model the outliers compared with Gaussian distribution. In practice, there are not only outliers but also all kinds noise in the training data, for example, asymmetric noise. To be more specific, these probability distributions are symmetric about the vertical axis ( $x=0$ ), and therefore cannot model the asymmetric noise that is common in practice.

The work based on the asymmetric distribution can model both the asymmetric and symmetric noise. Compared with Gaussian distribution, the asymmetric Laplace distribution is a heavy-tailed distribution. Thus, the algorithm studied based on this is robust to outliers. In addition, asymmetric Laplace distribution is different from Laplace distribution. Specifically, Laplace distribution is a special type of asymmetric Laplace distribution, the symmetric form. Thus, asymmetric Laplace distribution can model more general noise and therefore has stronger robustness. However, the application of asymmetric Laplace is still challenging because of the complex structure. Therefore, the expectation maximization (EM) algorithm framework is employed, and the LPV system identification algorithm can be deduced. To sum up, the contribution of the paper is summarized as follow:

1. A more robust parameters estimation method of LPV models is proposed. The algorithm identifies the parameters when the training data are polluted by both symmetric/asymmetric noise with outliers.

2. The noise is modeled by a more general probabilistic model, asymmetric Laplace distribution, which can both model the symmetric and asymmetric probabilities. In order to utilize the asymmetric Laplace distribution in parameters estimation, the asymmetric Laplace distribution with a shifted parameter is represented as the product of an exponential distribution and a normal distribution, and the posterior probability distribution of the random variable that obeys the exponential distribution is deduced.

3. The EM algorithm is utilized to deduced the parameters estimation algorithm. The random variable that obeys the exponential distribution is treated as the latent variable in the EM algorithm, and therefore the asymmetric Laplace distribution can be possible to model the noise.

The rest part of the paper is organized as follow. The LPV system identification problem is described in Section II, and the properties of asymmetric Laplace distribution are also introduced. In Section III, the parameters’ estimation algorithm of LPV models is deduced under the EM algorithm framework, and the EM algorithm is briefly introduced at the beginning. To illustrate the robustness and advantages of the proposed algorithm, the numerical example is present in Section IV. Finally, the conclusion is given in Section V.

Problem formulation

The LPV–OE model

The proposed algorithm focuses on the identification of LPV models, and the local models are output error (OE) models. The LPV–OE models are defined as equation (1).

$$\begin{array}{c}{x}_k=G(w_k,z^{-1})u_k\hfill \\ y_k=x_k+{\mathrm{\Delta }}_k\hfill \end{array}$$

(1)

where $y_k,u_k,w_k$ are the input, output and scheduling variable respectively. $z^{-1}$ is the time operator defined as equation (2).

$$x_{k-1}=z^{-1}{x}_k=G(w_k,z^{-1})u_k$$

(2)

The transfer function $G(w_k,z^{-1})$ is defined as equation (3).

$$G(w_k,z^{-1})=A^{†}(w_k,z^{-1})B(w_k,z^{-1})$$

(3)

The polynomial $A(w_k,z^{-1}),B(w_k,z^{-1})$ represent that the parameters of LPV–OE models vary with the scheduling variable $w_k$ .

$$\begin{array}{cc}A(w_k,z^{-1})\hfill & =1+\sum _{i=1}^{n_a}{a}_i\left(w_k\right)z^{-i}\hfill \\ \multicolumn{1}{c}{B(w_k,z^{-1})}& =\sum _{j=1}^{n_b}{b}_j\left(w_k\right)z^{-j}\hfill \end{array}$$

(4)

Normally, $a_i(w_k),b_j(w_k)$ are defined by a group of meromorphic functions, $\mathrm{\Xi },\mathrm{\Psi }$ , respectively.

$$\begin{array}{cc}{a}_i(w_k)\hfill & =a_0^i+\sum _{p=1}^{n_{\alpha }}{a}_P^i{\mathrm{\Xi }}_p(w_k),i=1,...,n_a\hfill \\ \multicolumn{1}{c}{b_j(w_k)}& =b_0^i+\sum _{q=1}^{n_{\beta }}{b}_q^i{\mathrm{\Psi }}_q(h_k),j=1,...,n_b\hfill \end{array}$$

(5)

where $n_{\alpha },n_{\beta }$ are the number of meromorphic functions that are the prior known information. $A^{†}(w_k,z^{-1})$ denotes the adjoint matrix of $A(w_k,z^{-1})$ .

The LPV–OE models that defined by equations (1)–(5) can be rewritten in equations (6)–(8), the autoregressive form.

$$x_k={\varphi }_k^T\theta$$

(6)

where

$${\varphi }_k={\left[\begin{array}{c}-x_{k-1}-{\mathrm{\Xi }}_1(w_k)x_{k-1}\cdots -{\mathrm{\Xi }}_{n_{\alpha }}(w_k)x_{k-1}\hfill \\ \multicolumn{1}{c}{-x_{k-2}-{\mathrm{\Xi }}_1(w_k)x_{k-2}\cdots -{\mathrm{\Xi }}_{n_{\alpha }}(w_k)x_{k-2}}\\ \multicolumn{1}{c}{\cdots \cdots \cdots \cdots }\\ \multicolumn{1}{c}{-x_{k-n_a}-{\mathrm{\Xi }}_1(w_k)x_{k-n_a}\cdots -{\mathrm{\Xi }}_{n_{\alpha }}(w_k)x_{k-n_a}}\\ \multicolumn{1}{c}{u_{k-1}-{\mathrm{\Psi }}_1(w_k)u_{k-1}\cdots -{\mathrm{\Psi }}_{n_{\beta }}(w_k)u_{k-1}}\\ \multicolumn{1}{c}{u_{k-2}-{\mathrm{\Psi }}_1(w_k)u_{k-2}\cdots -{\mathrm{\Psi }}_{n_{\beta }}(w_k)u_{k-2}}\\ \multicolumn{1}{c}{\cdots \cdots \cdots \cdots }\\ \multicolumn{1}{c}{u_{k-n_b}-{\mathrm{\Psi }}_1(w_k)u_{k-n_b}\cdots -{\mathrm{\Psi }}_{n_{\beta }}(w_k)u_{k-n_b}}\end{array}\right]}^T$$

(7)

$$\theta ={\left[\begin{array}{c}{a}_0^1{a}_1^1\cdots a_{n_{\alpha }}^1\hfill \\ \multicolumn{1}{c}{a_0^2{a}_1^2\cdots a_{n_{\alpha }}^2}\\ \multicolumn{1}{c}{\cdots \cdots \cdots \cdots }\\ \multicolumn{1}{c}{a_0^{n_a}{a}_1^{n_a}\cdots a_{n_{\alpha }}^{n_a}}\\ \multicolumn{1}{c}{b_0^1{b}_1^1\cdots b_{n_{\beta }}^1}\\ \multicolumn{1}{c}{b_0^2{b}_1^2\cdots b_{n_{\beta }}^2}\\ \multicolumn{1}{c}{\cdots \cdots \cdots \cdots }\\ \multicolumn{1}{c}{b_0^{n_b}{b}_1^{n_b}\cdots b_{n_{\beta }}^{n_b}}\end{array}\right]}^T$$

(8)

The aim of this paper is to identify the parameters of LPV–OE models, $\theta$ , based on the collected input/output data when the training data is polluted by outliers and asymmetrically distributed noise. In order to get the robustness, the noise ${\mathrm{\Delta }}_k$ in equation (1) is assumed to obey the asymmetric Laplace distribution.

Asymmetric Laplace distribution

The probability density function of a p dimensional random variable that obeys centralized asymmetric Laplace (CAL) distribution is shown in equation (9).

$$f(\mathbf{v}|\alpha ,\mathbf{\Sigma })=\frac{2\exp \left\{{\mathbf{v}}^{\prime }{\mathbf{\Sigma }}^{-1}\alpha \right\}}{{(2\pi )}^{p/2}|\mathbf{\Sigma }{|}^{1/2}}\times {\left(\frac{{\mathbf{v}}^{\prime }{\mathbf{\Sigma }}^{-1}\mathbf{v}}{2+{\mathit{\alpha }}^{\prime }{\mathbf{\Sigma }}^{-1}\alpha }\right)}^{v/2}{K}_v(u)$$

(9)

where, $v=(2-p)/2$ , $u=\sqrt{(2+{\alpha }^{\prime }{\mathbf{\Sigma }}^{-1}\alpha )({\mathbf{v}}^{\prime }{\mathbf{\Sigma }}^{-1}\mathbf{v})}$ . $\mathbf{\Sigma }$ and $\alpha \in {\mathbb{R}}^p$ denote the scale matrix and the skewness respectively. In this paper, symbol $V~\mathcal{C}\mathcal{A}{\mathcal{L}}_p(\alpha ,\mathbf{\Sigma })$ is employed to denote that a p dimensional random variable obeys CAL distribution.

However, CAL distribution cannot be utilized in clustering and system identification directly because CAL distribution forces all components connecting at the same region. Thus, by considering the shift parameter $\mu \in {\mathbb{R}}^p$ with the same dimension to CAL distribution, the random variable X that obeys shifted (non-centralized) asymmetric Laplace (SAL) distribution can be defined as $X=(V+\mu )$ . Herein, symbol $X~\mathcal{S}\mathcal{A}{\mathcal{L}}_p(\alpha ,\mathbf{\Sigma },\mu )$ is employed to denote that p dimension random variable X obeys SAL distribution, and the probability density function is equation (10).

$$\xi (\mathbf{x}|\alpha ,\mathbf{\Sigma },\mu )=\frac{2\exp \{{(\mathbf{x}-\mu )}^{\prime }{\mathbf{\Sigma }}^{-1}\alpha \}}{\begin{array}{c}\multicolumn{1}{c}{{(2\pi )}^{p/2}|\mathbf{\Sigma }{|}^{1/2}}\end{array}}\times {\left(\frac{\delta (\mathbf{x},\mu |\mathbf{\Sigma })}{2+{\alpha }^{\prime }{\mathbf{\Sigma }}^{-1}\alpha }\right)}^{v/2}{K}_v(u)$$

(10)

where $u={\left((2+{\alpha }^T{\Sigma }^{-1}\alpha )\delta (\mathbf{x},\mu |\Sigma )\right)}^{1/2}$ and $\delta (x,\mu |\Sigma )=(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )$ ; $v=(2-p)/2$ ; $\mathbf{\Sigma }$ and $\alpha$ are defined as before. The paper emphasizes the identification of the Single Input Single Output system that is one-dimensional problem.

To illustrate the asymmetry of SAL distribution, the probability density function of the SAL distribution with skewness $5$ shift $5$ and scale $100$ is shown in Figure 1.

Figure 1.

The probability density function of the SAL distribution.

However, because of the complexity of SAL distribution, it still cannot be utilized directly like Gaussian distribution. ¹⁹ has studied SAL distribution in detail that the random variable $V~\mathcal{C}\mathcal{A}{\mathcal{L}}_p(\alpha ,\mathbf{\Sigma })$ can be generated by equation (11).

$$V=z\alpha +\sqrt{z}Y$$

(11)

where Z obeys exponential distribution with rate parameters $1$ ; Y is the random variable Y obeys normal distribution with zero mean and variance $\Sigma$ , immediate $Y~\mathcal{N}(0,\mathbf{\Sigma })$ . Z and Y are independent.

Thus, the random variable $X~\mathcal{S}\mathcal{A}{\mathcal{L}}_p(\alpha ,\mathbf{\Sigma },\mu )$ can be illustrated by adding shifted parameter $\mu$ to equation (11).

$$X=\mu +z\alpha +\sqrt{z}Y$$

(12)

In other words, the random variable that obeys SAL distribution ultimately obeys normal distribution when the random variable z is observable.

$$X|Z=z~\mathcal{N}(\mu +z\alpha ,z\mathbf{\Sigma })$$

(13)

In the system identification problem described in Section, the noise obeys the CAL distribution, and therefore, the output obeys the SAL distribution with $\mu =x_k$ . In addition, when Z in SAL distribution is observable, the output radically obeys the normal distribution and therefore the identification problem can be solved. However, it is impossible to monitor the parameter of noise. According to Bayes formula, the probability density function of z can be represent by equation (14).

$$\begin{array}{l}{f}_Z(z|\mathbf{X}=\mathbf{x})=\frac{f_{\mathbf{X}}(\mathbf{x}|Z=z)h(z)}{f_{\mathbf{X}}(\mathbf{x})}\\ =\frac{z^{v-1}}{2}{\left(\frac{\delta (\mathbf{x},\mu |\Sigma )}{2+{\alpha }^{\prime }{\Sigma }^{-1}\alpha }\right)}^{-v/2}\times \frac{\exp \left\{-\frac{1}{2z}\delta (\mathbf{x},\mu |\Sigma )-\frac{z}{2}\left(2+{\alpha }^{\prime }{\Sigma }^{-1}\alpha \right)\right\}}{K_v(\sqrt{\left(2+{\alpha }^{\prime }{\Sigma }^{-1}\alpha \right)\delta (\mathbf{x},\mu |\Sigma )})}\end{array}$$

(14)

It can be seen that equation (14) is fundamentally Generalized inverse Gaussian (GIG) distribution with $a\equiv 2+{\alpha }^{\prime }{\mathbf{\Sigma }}^{-1}\alpha$ and $b\equiv \delta (\mathbf{x},\mu |\mathbf{\Sigma })$ which has been studied in detail. In this paper, the expectation of the random variable and the reciprocal of the random variable that obeys GIG distribution are useful and are shown in equations (15) and (16).

$$\mathbb{E}[X]=\sqrt{\frac{b}{a}}{R}_v(\sqrt{\mathrm{ab}})$$

(15)

$$\mathbb{E}[1/X]=\sqrt{\frac{a}{b}}{R}_v(\sqrt{\mathrm{ab}})-\frac{2v}{b}$$

(16)

where $R_v(z):=K_{v+1}(z)/K_v(z)$ .

Robust parameters estimation algorithm for LPV–OE systems

Brief introduction of the EM algorithm

The EM algorithm is an iterative method of the maximum likelihood estimation problem and its flexibility to deal with the missing part (hidden state) attracts wide attention. In the EM algorithm, Expectation step (E-Step) and Maximization step (M-Step) run alternatively until convergence. In the E-Step, the conditional expectation of the log likelihood function with respect to the latent state is calculated. In addition, in the M-Step, the conditional expectation that calculated in the E-Step is maximized. Several methods have been studied to achieve it, including the Newton-Raphson method, Derivative analysis, etc. The procedures of the EM algorithm is summarized below:

• (1) Initialization. According to the identification target, the complete data set and the parameters to be estimated are divided into observable data set $D_{\mathrm{obs}}$ , missing data set $D_{\mathrm{mis}}$ and the parameters data $\mathrm{\Theta }$ . In addition, maximum number of iteration times and expected residual are also defined before the iteration.

• (2) E-Step: The conditional expectation of the log likelihood function with respect to the latent state.

$$Q(\mathrm{\Theta }|{\mathrm{\Theta }}^l)=E_{D_{\mathrm{mis}}|D_{\mathrm{obs}}}\{\mathrm{logP}(D_{\mathrm{mis}},D_{\mathrm{obs}}|\mathrm{\Theta })\}$$

(17)

• (3) M-Step: Maximize the $Q(\mathrm{\Theta }|{\mathrm{\Theta }}^l)$ by finding suitable ${\mathrm{\Theta }}^{l+1}$

$$Q({\mathrm{\Theta }}^{l+1}|{\mathrm{\Theta }}^l)>Q({\mathrm{\Theta }}^l|{\mathrm{\Theta }}^l)$$

(18)

• (4) Let ${\mathrm{\Theta }}^l={\mathrm{\Theta }}^{l+1}$ and running steps (2) and (3) alternatively until the parameters converge.

Mathematical formulation of the robust parameters estimation approach

The data-driven system identification algorithm is proposed in order to estimate the parameters of LPV–OE systems. Compared with mechanism models, the proposed algorithm only utilizes limited prior knowledge and the input-output data to establish the process models. However, there is too much information that need to be identified, for example, the model orders and parameters, the meromorphic functions, etc. The paper emphasizes on the parameters estimation, thus, some information about the structures of LPV–OE models is assumed to be known or fixed a priori.

• Assumption 0.1. The structures of the LPV–OE systems are assumed based on the prior information. Thus, the ground truth orders of the LPV–OE models to be identified are $n_a$ and $n_b$ , and the meromorphic functions to construct the polynomials $A(\omega ,z^{-1})$ and $B(\omega ,z^{-1})$ are pre-selected.

• Remark 1. The paper focus on the parameters estimation algorithm, therefore, the structure information is assumed to be known. To find a suitable structure, especially for some complex processes, identifying the parameters and the structures iteratively is a good way. After the structures have been predefined, the parameters can be estimated by utilizing the proposed algorithm. If the identified LPV–OE model cannot meet the index, we can adjust the structures and identify the parameters repeatedly until the performing is satisfactory.

Under this assumption, the only information to be determined is the parameters of models if the structures of LPV–OE models have been identified or known. Thus, by using the collected input and output data, estimating the parameters of LPV–OE models is the emphasis. In particular, the difficulty of identification depends on the quality of the collected data. In this section, the robust parameters estimation algorithm will be studied and introduced in detail, and the robustness to outliers and asymmetric noise will be illustrated.

As a statistical learning method, the proposed algorithm will be formulated based on the EM algorithm framework. The data-driven parameters estimation algorithm identifies the parameters of LPV–OE systems by extracting the features from input $\{u_k{\}}_{k=1:L}$ output $\{y_k{\}}_{k=1:L}$ and scheduling $\{{\omega }_k{\}}_{k=1:L}$ data, all of which are collected as the observed data set $D_{\mathrm{obs}}=\{y_{1:L},{\omega }_{1:L},u_{1:L}\}$ in the EM algorithm framework. Moreover, there is a significant variable that cannot be observed, and therefore, it is the latent variable that obeys exponential distribution in the SAL distribution. In the EM algorithm framework, the expectation of the likelihood function respect to the latent variable is maximized. Thus, the unobservable data set $D_{\mathrm{mis}}=\{{ϵ}_{1:L}$ } can be formed. Finally, the parameters set $\mathrm{\Theta }$ collects all the parameters to be identified, including the LPV–OE models parameters $\theta$ and the parameters of noise $\alpha ,\Sigma$ , namely $\mathrm{\Theta }=\{\theta ,\alpha ,\Sigma \}$ .

E-Step: Firstly, the likelihood function of all data set including observed and unobserved variables $D_{\mathrm{obs}}$ and $D_{\mathrm{mis}}$ is calculated in the E-step.

$$\begin{array}{c}p(D_{\mathrm{obs}},D_{\mathrm{mis}}|\mathrm{\Theta })=p(y_{1:L},u_{1:L},{\omega }_{1:L},{ϵ}_{1:L}|\mathrm{\Theta })\hfill \\ =p(y_{1:L}|u_{1:L},{\omega }_{1:L},{ϵ}_{1:L},\mathrm{\Theta })\times p({ϵ}_{1:L}|u_{1:L},{\omega }_{1:L},\mathrm{\Theta })\times C\hfill \end{array}$$

(19)

where $C=p(u_{1:l,\_1:l|\mathrm{\Theta }})$ is a constant since both $u_{1:L}$ and ${\omega }_{1:L}$ can be observed. Thus, the log-likelihood function of equation (19) can be deduced.

$$\begin{array}{c}\mathrm{logp}(D_{\mathrm{obs}},D_{\mathrm{mis}}|\mathrm{\Theta })\hfill \\ =\mathrm{logp}(y_{1:L}|u_{1:L},{\omega }_{1:L},{ϵ}_{1:L},\mathrm{\Theta })\hfill \\ \multicolumn{1}{c}{+\mathrm{logp}({ϵ}_{1:L}|u_{1:L},{\omega }_{1:L},\mathrm{\Theta })+\mathrm{logC}}\\ =\sum _{k=1}^L\mathrm{logp}(y_k|u_{u:k-1},{\omega }_k,{ϵ}_k,\mathrm{\Theta })\hfill \\ \multicolumn{1}{c}{+\sum _{k=1}^L\mathrm{logp}({ϵ}_k|\mathrm{\Theta })+\mathrm{logC}}\end{array}$$

(20)

The Q-function in the EM algorithm, conditional expectation of the log-likelihood function with respect to the latent variables, can be represented by equation (21).

$$\begin{array}{c}Q(\mathrm{\Theta }|{\mathrm{\Theta }}^s)\hfill \\ =E_{{ϵ}_{1:L}|D_{\mathrm{obs}}}\{\sum _{k=1}^L\mathrm{logp}(y_k|u_{u:k-1},w_k,{ϵ}_k,\mathrm{\Theta })+\sum _{k=1}^L\mathrm{logp}({ϵ}_k|\mathrm{\Theta })+C_2\}\hfill \\ =\sum _{k=1}^L\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)\times \mathrm{logp}(y_k|u_{1:k-1},w_k,{ϵ}_k,\mathrm{\Theta })d{ϵ}_k\hfill \\ +\sum _{k=1}^L\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)\times \mathrm{logp}({ϵ}_k|\mathrm{\Theta })d{ϵ}_k+C_3\hfill \end{array}$$

(21)

In order to deduce the Q-function, the unknown quantities in equation (21) have to be studied, including

1. $p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)$ ;

2. $\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)\times \mathrm{logp}(y_k|u_{1:k-1},w_k,{ϵ}_k,\mathrm{\Theta })d{ϵ}_k$ ;

3. $\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)\times \mathrm{logp}({ϵ}_k|\mathrm{\Theta })d{ϵ}_k$ .

The posterior probability distribution of the latent variable ${ϵ}_{1:L}$ can be deduced by using Bayes formula.

$$\begin{array}{c}\multicolumn{1}{c}{p({ϵ}_k|y_{1:L},u_{1:L},w_{1:L},{\mathrm{\Theta }}^s)}\\ \multicolumn{1}{c}{=\frac{p(y_k|{ϵ}_k,u_{1:k-1},w_{1:k},{\mathrm{\Theta }}^s)p({ϵ}_k|{\mathrm{\Theta }}^s)}{p(y_k|y_{1:k-1},u_{1:k-1},w_{1:k},{\mathrm{\Theta }}^s)}}\end{array}$$

(22)

It has been introduced in detail in Section that the SAL distribution can be reduced to Gaussian distribution when the latent variable $ϵ$ are observable. Thus, when the latent variable is assumed to be observed, the posterior probability distribution of the output data $y_k$ is Gaussian distribution that is represented by equation (23).

$$\begin{array}{c}\multicolumn{1}{c}{p(y_k|{ϵ}_k,u_{1:k-1},w_{1:k},{\mathrm{\Theta }}^s)=\frac{1}{\sqrt{2\pi {ϵ}_k\Sigma }}\exp }\\ \multicolumn{1}{c}{(-\frac{{(y_k-G^s(w_k,z^{-1})u_k)}^2}{2{ϵ}_k\Sigma })}\end{array}$$

(23)

When it comes to the condition that the latent variable cannot be observed, the output data obeys SAL distribution of which the probability density function is equation (10), and the random variable ${\epsilon }_k$ obeys the exponential distribution with the rate of 1.

In this way, the unknown items of the posterior probability function of latent variable $\epsilon$ are obtained and can be deduced in equation (24).

$$\begin{array}{c}p({ϵ}_k|y_{k-1:1},u_{k-1:1},w_k,{\mathrm{\Theta }}^l)=\frac{{ϵ}_k^{-\frac{1}{2}}}{2}(\frac{\delta (y_k,G^s(w_k,z^{-1})u_k|\Sigma )}{2+{\alpha }^{\mathrm{T}}{\Sigma }^{-1}\alpha }{)}^{-\frac{1}{4}}\hfill \\ \times \frac{\exp \{-\frac{1}{2z_k}\delta (y_k,G^s(w_k,z^{-1})u_k|\Sigma )-\frac{z_k}{2}(2+{\alpha }^{\mathrm{T}}{\Sigma }^{-1}\alpha )\}}{K_{\frac{1}{2}}(\sqrt{(2+{\alpha }^{\mathrm{T}}{\Sigma }^{-1}\alpha )(\delta (y_k,G^s(w_k,z^{-1})u_k|\Sigma ))})}\hfill \end{array}$$

(24)

where $\delta (y_k,G^s(w_k,z^{-1})u_k|\Sigma )=(y_k-G^s(w_k,z^{-1})u_k{)}^T{\Sigma }^{-1}(y_k-G^s(w_k,z^{-1})u_k)$ . $K_v$ is the third kind of modified Bessel function with the exponent of v.

Comparing equation (24) with the generalized inverse Gaussian (GIG) distribution, it is clear that equation (24) is the GIG distribution with $a\equiv 2+{\alpha }^{\prime }{\Sigma }^T\alpha$ , $b\equiv \delta (y_k,G^s(w_k,z^{-1})u_k|\Sigma )$ . Herein, the properties of the GIG distribution are utilized to calculate the means of the variable $\mathbb{E}(x)$ and the reciprocal of the variable $\mathbb{E}(1/x)$ .

$$\mathbb{E}\left[z_k\right]=\sqrt{\frac{b}{a}}{R}_{\frac{1}{2}}(\sqrt{\mathrm{ab}})$$

(25)

$$\mathbb{E}\left[\frac{1}{z_k}\right]=\sqrt{\frac{a}{b}}{R}_t(\sqrt{\mathrm{ab}})-\frac{1}{b}$$

(26)

where $R_{1/2}(x):=K_{3/2}(x)/K_{1/2}(x)$ , $\mathbb{E}\left[z_k\right]\stackrel{\mathrm{\Delta }}{=}{\mathbb{E}}_1$ , $\mathbb{E}[1/z_k]\stackrel{\mathrm{\Delta }}{=}{\mathbb{E}}_2$ .

Since all the unknown items in equation (21) have been deduced, by substituting equations (25) and (26) into (21), we can deduce the Q-function. For the sake of presentation, the Q-function is divided into the sum of three terms. In addition, in the M-step in the EM algorithm, the Q-function will be maximized by finding suitable parameters $\mathrm{\Theta }$ , and therefore the constant term in the Q-function that does not matter is neglected. In this way, the Q-function is shown as the sum of two terms, $Q(\mathrm{\Theta }|{\mathrm{\Theta }}^S)=Q_1(\mathrm{\Theta }|{\mathrm{\Theta }}^S)+Q_2(\mathrm{\Theta }|{\mathrm{\Theta }}^S)$ .

$$\begin{array}{c}\multicolumn{1}{c}{Q_1(\mathrm{\Theta }|{\mathrm{\Theta }}^S)=\sum _{k=1}^L\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)\times \mathrm{logp}(y_k|u_{1:k-1},w_k,{ϵ}_k,\mathrm{\Theta })d{ϵ}_k}\\ \multicolumn{1}{c}{=\sum _{k=1}^L\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)}\\ \multicolumn{1}{c}{\times \left\{\log \frac{1}{\sqrt{2\pi }\sqrt{{ϵ}_k\Sigma }}\exp (-\frac{{(y_k-G^s(w_k,z^{-1})u_k-{ϵ}_k\alpha )}^2}{2{ϵ}_k\Sigma })\right\}}\\ \multicolumn{1}{c}{=\sum _{k=1}^L-\frac{1}{2}\log (2\pi )-\frac{1}{2}\log (\Sigma )-\frac{1}{2}\mathbb{E}(\log ({ϵ}_k))-}\\ \multicolumn{1}{c}{{\mathbb{E}}_2\frac{1}{2\Sigma }\left\{{(y_k-G^s(w_k,z^{-1})u_k-{ϵ}_k\alpha )}^2\right\}}\end{array}$$

(27)

$$\begin{array}{c}\multicolumn{1}{c}{Q_2(\mathrm{\Theta }|{\mathrm{\Theta }}^S)=\sum _{k=1}^L\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)\times \mathrm{logp}({ϵ}_k|\mathrm{\Theta })d{ϵ}_k}\\ \multicolumn{1}{c}{=\sum _{k=1}^L\int p({ϵ}_k|D_{\mathrm{obs}},{\mathrm{\Theta }}^s)\times \log (\exp (-{ϵ}_k))}\\ \multicolumn{1}{c}{=\sum _{k=1}^L-{\mathbb{E}}_1}\end{array}$$

(28)

M-Step: Deducing the iteration formulas of the parameters $\mathrm{\Theta }$ to maximize the log likelihood function (Q-function).

In order to deduce the skewness $\alpha$ , as there is no connection between the second item of the Q-function $Q_2(\mathrm{\Theta }|{\mathrm{\Theta }}^S)$ and the skewness $\alpha$ , taking the derivative of $\alpha$ with respect to $Q_1(\mathrm{\Theta }|{\mathrm{\Theta }}^S)$ and setting it equal to zero can get the iteration formula of the skwness $\alpha$ , which is shown in equation (29).

$${\alpha }^{s+1}=\frac{1}{L}\frac{\sum _{k=1}^L(y_k-G^s(w_k,z^{-1})u_k)}{\sum _{k=1}^L{\alpha }^s}$$

(29)

When it comes to the iteration formula of the parameters of LPV models $\theta$ , the Newton–Raphson method is utilized to search for the suitable $\theta$ that maximizes the first item of the Q-function $Q_1(\mathrm{\Theta }|{\mathrm{\Theta }}^S)$ since there is also no connection between $Q_2(\mathrm{\Theta }|{\mathrm{\Theta }}^S)$ and $\theta$ .

$${\theta }^{s+1}={\theta }^s-{(\frac{d^{2}Q(\theta )}{d{\theta }^2}+\mu I)}^{-1}\frac{\mathrm{dQ}(\theta )}{d\theta }{|}_{\theta ={\theta }^s}$$

(30)

where the first-order derivative of $Q_1(\mathrm{\Theta }|{\mathrm{\Theta }}^S)$ is deduced in equation (31), and the second-order derivative of it is deduced in equation (32).

$$\frac{\mathrm{dQ}(\theta )}{d\theta }=\frac{2}{L}\left\{\sum _{k=1}^L{\mathbb{E}}_2(y_k-{{\hat{\varphi }}_k}^T\theta -{\mathbb{E}}_1\alpha ){\hat{\varphi }}_k\right\}$$

(31)

$$\frac{d^{2}Q(\theta )}{d{\theta }^2}=\frac{2}{L}\left\{\sum _{k=1}^L{\mathbb{E}}_2{{\hat{\varphi }}_k}^T{\hat{\varphi }}_k\right\}$$

(32)

where ${\hat{\varphi }}_k$ is the estimation of ${\varphi }_k$ .

$${\hat{x}}_k=G^s(w_k,z^{-1})u_k$$

(33)

The complete LPV–OE models’ parameters estimation algorithm is summarized as in Algorithm 1.

Algorithm 1: LPV Models’ Parameters Estimation Algorithm
Input: Observed data set $D_{\mathrm{obs}}=\{y_{1:L},u_{1:L},w_{1:L}\}$
Output: Parameters to be estimated $\mathrm{\Theta }=\{\theta ,\alpha ,\Sigma \}$
1 Initializing the number of iterations $l=1$ ; Setting the target residuals E and maximum number of iterations L;
2 Initializing the parameters to be estimated ${\theta }^l,{\alpha }^l$ ;
3 forl = from 1 to $\mathrm{Res}<E$ or $l>L$ do
4  E-Step:;
5  Calculating the posterior probability density of $z_k$ based on Eq.(24);
6  Calculating the means of latent variables $1/z_i$ and $z_i$ , ${\mathrm{E}}_1$ and ${\mathrm{E}}_2$ based on Eq.(25) and Eq.(26) respectively;
7  M-Step:;
8  Calculating the skewness ${\alpha }^l$ based on Eq.(29);
9  Calculating the parameters ${\theta }^l$ of LPV models based on skewness ${\alpha }^l$ and Eq.(30,31,32);
10  Calculating the noise-free output $\widehat{\{}{x}_k{\}}_{k=1:L}$ based on parameters ${\theta }^l$ and skewness ${\alpha }^l$ ;
11  Calculating the residual $\mathrm{Res}=y-\hat{x}$ ;
12 end
13 Return $\theta$ , $\alpha$ ;

Simulation example

Considering the following LPV–OE model

$$\begin{array}{c}\multicolumn{1}{c}{C(v_k,d^{-1})z_k=E(v_k,d^{-1})u_k}\\ \multicolumn{1}{c}{y_k=z_k+{\mathrm{\Delta }}_k}\end{array}$$

(34)

where

$$\begin{array}{c}\multicolumn{1}{c}{C(v_k,d^{-1})=1+c_1\left(v_k\right)d^{-1}+c_2\left(v_k\right)d^{-2}}\\ \multicolumn{1}{c}{c_1\left(v_k\right)=1-0.5v_k-0.1v_k^2}\\ \multicolumn{1}{c}{c_2\left(v_k\right)=0.5-0.7v_k-0.1v_k^2}\end{array}$$

(35)

and

$$\begin{array}{c}\multicolumn{1}{c}{E(v_k,d^{-1})=e_1\left(v_k\right)d^{-1}+e_2\left(v_k\right)d^{-2}}\\ \multicolumn{1}{c}{e_1\left(v_k\right)=0.5-0.4v_k+0.01v_k^2}\\ \multicolumn{1}{c}{e_2\left(v_k\right)=0.5-0.3v_k-0.02v_k^2}\end{array}$$

(36)

where $z_k$ is the noise-free output. $\mathrm{\Delta }$ is the noise. As the robustness to asymmetric noise is emphasized in this paper, the noise is assumed to obey the asymmetric Laplace distribution comparing with the symmetric distribution, for example, Gaussian distribution (white noise) and Laplace distribution. $v_k$ in equation (34) is the scheduling variable in LPV models. In addition, the parameters $c_1,c_2,e_1,e_2$ vary with the scheduling variable, and the relationships between the time-varying parameters and the scheduling variable is defined in equations (35) and (36).

In the simulation example, Gaussian white noise is the input signal which is a common and effective input signal when the hardware conditions permit, and it can be utilized to fully simulate the process and obtain the data needed for identification. At the same time, in the single-input single-output system, the input signal is also the scheduling variable immediate $u_k=v_k$ .

Moreover, outliers are common in the practice. Since the SAL distribution is a special type of Laplace distribution, the proposed algorithm is also robust to the outliers. Thus, the output data is also polluted by outliers in order to verify the robustness to outliers.

The training data are shown in Figure 2 where 10% of the output data are polluted by outliers that are 10 times bigger than the noise-free output. By using the polluted data and the Algorithm1, the parameters of the LPV–OE model can be identified, and the comparison of the estimated parameters with the ground truth parameters are shown in Figure 3. The convergence procedures of the parameters are also illustrated in it. The value of the parameters at each iteration are shown in blue, and the red lines are the ground truth parameters. In the experiment, with the initialization parameters (2), all the parameters converge to the ground truth values and remain stable after 100 iterations. Thus, the proposed algorithm has local convergence, and it can identify the parameters of LPV–OE models accurately when the training data is polluted by asymmetric noise and outliers.

• Remark 2. The robustness to the outliers and asymmetric noise comes from equation (31). In terms of outliers, compared with the algorithm (equation (37)) based on Gaussian distribution, ²⁰ there is the weight coefficient ${\mathbb{E}}_2$ before the residual $y_k-{{\hat{\varphi }}_k}^T\theta$ in equation (31) that can identify the outliers in the training data. Almost by definition equation (26), when the estimation parameters closer to the ground truth parameters after several iteration, the weight coefficients decrease only if the deviation between estimation output and collected output is large, and a larger deviation means that the sampling point has a greater likelihood of becoming outliers. Thus, the sampling points with larger deviation are given small weight coefficient to weaken the influence of outliers.

Figure 2.

The input data (scheduling variable) and output data.

Figure 3.

The convergence procedures.

As for the robustness to the asymmetric noise, the item ${\mathbb{E}}_1\alpha$ in equation (31) represents the identified skewness, and therefore eliminates the effects of asymmetry. Of course, $\alpha$ equals to zeros when the data is only polluted by symmetric noise.

$$\frac{\mathrm{dQ}(\theta )}{d\theta }=\frac{2}{L}\left\{\sum _{k=1}^L(y_k-{{\hat{\varphi }}_k}^T\theta ){\hat{\varphi }}_k\right\}$$

(37)

Most of the robust LPV models identification methods focus on the robustness to outliers but various noise. Specifically, the parameters estimation algorithms are deduced based on Laplace distribution ¹⁷ and Student’s T distribution ¹⁵ since these two distributions are heavy-tailed distributions. However, they can only model the symmetric noise. If the training data is polluted by asymmetric noise, the estimated parameters are bias. Fortunately, the proposed algorithm is robust to both symmetric/asymmetric noise with outliers.

The numerical example is also employed to verify the robustness to asymmetric noise. In the meanwhile, two parameters estimation algorithms proposed in Yang and Yang ²¹ and Liu et al. ²² are employed when the training data is polluted by asymmetric noise with outliers. The first algorithm is deduced based on Gaussian distribution. Thus, the algorithm is not robust. In the second algorithm, Laplace distribution is employed to model the noise. So, the algorithm is robust to the outliers. However, these two algorithms are not robust to asymmetric noise. The parameters identified by the three algorithms as well as the ground truth parameters are shown in Table 1.

Table 1.

The comparison of estimated parameters with true parameters.

	$c_1$
	$c_{1,1}$	$c_{1,2}$	$c_{1,3}$
Ground truth	1	−0.5	−0.1
The algorithm based on SAL distribution	0.9983	−0.4901	−0.1066
The algorithm based on Gaussian distribution	1.0284	−0.5134	−0.1389
The algorithm based on Laplace distribution	0.9737	−0.4875	−0.0959
	$e_1$ $e_{2,1}$
	$e_{1,1}$	$e_{1,2}$	$e_{1,3}$
Ground truth	0.5	−0.4	0.01
The algorithm based on SAL distribution	0.4880	−0.4227	0.0245
The algorithm based on Gaussian distribution	0.3960	−0.5924	0.1245
The algorithm based on Laplace distribution	0.5496	−0.3937	−0.0849
	$c_2$
	$c_{2,1}$	$c_{2,2}$	$c_{2,3}$
Ground truth	0.5	−0.7	−0.1
The algorithm based on SAL distribution	0.4959	−0.6907	−0.0991
The algorithm based on Gaussian distribution	0.4761	−0.7678	−0.1344
The algorithm based on Laplace distribution	0.4756	−0.7131	−0.0849
	$e_2$
	$e_{2,1}$	$e_{2,2}$	$e_{2,3}$
Ground truth	0.5	−0.3	0.02
The algorithm based on SAL distribution	0.4879	−0.1854	−0.0279
The algorithm based on Gaussian distribution	0.4563	−0.3413	−0.0493
The algorithm based on Laplace distribution	0.4786	−0.1575	−0.0496

From the results, it can be seen that the performance of the proposed algorithm is better than the other two algorithms. The robustness of the algorithm based on Gaussian is the worst because Gaussian distribution a symmetric distribution and it cannot model outliers and asymmetric noise in the data. The algorithm based on Laplace distribution performs better than the counterpart based on Gaussian distribution, since Laplace distribution is a heavy-tailed distribution. Its heavy tail can model the outliers in the training data. However, it is still a symmetric distribution, so it is not robust to the asymmetric noise. The histogram of the noise in the numerical example is shown in Figure 4. The vast majority of noise are greater than zero. The asymmetry of this noise cannot be modeled by Laplace distribution but it is common in industries. Fortunately, the CAL distribution can model the noise and therefore, the proposed algorithm has stronger robustness and identifies the parameters accurately.

Figure 4.

The histogram of the noise.

Conclusion

In the work, the parameters estimation algorithm of LPV–OE models is deduced based on asymmetric Laplace distribution and the EM algorithm. The proposed algorithm is robust to both asymmetric noise and outliers, and the effectiveness of the proposed algorithm is verified by a numerical example.

The advantages of the proposed algorithm come from the asymmetric Laplace distribution. Compared with the Gaussian and the Laplace distribution, the asymmetric Laplace distribution has heavy tail and asymmetry respectively. Thus, it can model more general noise, both the asymmetric noise and symmetric noise with outliers. To make it is possible to utilized asymmetric Laplace distribution, it is represented as the product of exponential distribution and Gaussian distribution. When the random variable that obeys exponential distribution is selected as the latent variable in the EM algorithm framework. The posterior distribution of the latent variable can be deduced and be represented as the structure of GIG distribution. In this way, the parameters estimation algorithm can be deduced. Finally, by using a numerical example, it is verified that the proposed algorithm can identify the parameters of LPV–OE models when the training data is polluted by asymmetric noise and outliers.