EM-based Global Identification of LPV ARX Models with a Noisy Scheduling Variable Xin Liu
Zeyuan Wang
Xianqiang Yang
Research Institute of Intelligent Control and System Harbin Institute of Technology Harbin, Heilongjiang 150001, China Email:
[email protected] Research Institute of Intelligent Control and System Harbin Institute of Technology Harbin, Heilongjiang 150001, China Email: wangzeyuan
[email protected] Research Institute of Intelligent Control and System Harbin Institute of Technology Harbin, Heilongjiang 150001, China Email:
[email protected] Abstract—This paper deals with the identification problem of a linear parameter varying (LPV) model with a noisy (uncertain) scheduling variable. An autoregressive exogenous (ARX) model is used to describe the system behavior and a nonlinear statespace model is used to model dynamics of the noisy scheduling variable. The uncertain scheduling variable is estimated by using a particle smoother, which is based on sequential Monte-Carlo (SMC) method, and all the unknown parameters in the LPV ARX model are identified under the framework of an expectationmaximization (EM) algorithm based on the intput-output data of the system. A numerical example and a mechanical unit are used to validate the proposed approach.
I. I NTRODUCTION A part of the industrial processes are modeled by nonlinear differential equations because of their inherent nonlinear properties and the identification of these systems is still a long term problem which needs to be studied in the next few decades [1]–[6]. Generally, most complex industrial processes are designed to perform certain task at several pre-defined operating points. The operating point variable is also referred to as scheduling variable. The transition dynamics between the operating points excite the nonlinearity of many industrial processes which poses great challenges for process modeling [7]–[9]. Many researchers have made their own contributions to identify those complex systems. A little part of the general identification methods are highlighted below and the more detailed review can be found in [10], [11]. The first principle method is often utilized to model nonlinear systems and the input-output data of the systems are recorded to identify the parameters of the models as well. Furthermore, nonlinear least squares [13] and prediction error method (PEM) [12] are proposed and regarded as the advanced approaches for solving parameter estimation problems. Unfortunately, a complete understanding of the process dynamics and high computation costs are required in this method, which makes it not realizable for complex industrial processes [14]. Data-driven modeling method is regarded as an effective way to represent the dynamics of complex processes. This This work was supported in part by the First-class General Financial Grant from the China Postdoctoral Science Foundation (Grant No. 2015M80264), and Financial Assistance under Heilongjiang Postdoctoral Fund (Grant No. LBH-Z15080). (Corresponding author: Xianqiang Yang.)
l-)))
kind of methods, such as neural networks (NN) [15] and support vector regression (SVR) [16], do not require any priori knowledge about the modeled systems. Nowadays, a class of model which has linear model structure and varying model parameters, referred to as linear parameter varying (LPV) model, is widely used in the area of nonlinear systems identification. Many researchers have paid their attention to the LPV model and identification approaches for the LPV model [17]–[21]. Generally, the parameters of the LPV model are polynomial functions of the scheduling variables and they are not constant which differs from the linear time invariant (LTI) model. The persistent excitation of the nonlinear systems throughout the whole operating range is required in LPV model identification to estimate the constant coefficients in the polynomial functions, however, this is overcost and unrealistic for majority of the industrial processes. To overcome this obstacle, Zhu [22] put forward a LPV method for nonlinear process identification in which some linear models were identified by using the data at several pre-defined operating-points and the LPV model was finally derived by interpolating the pre-identified linear models. In other words, the LPV model consists of several linear local models and weighting functions. Only the persistent excitation around several pre-defined working-points is need in this method. It is worth mentioning that the pre-chosen workingpoints and the dynamic transition between different workingpoints should be precisely known a priori in this method. Jin et al. [7] proposed a probability-based method for LPV model identification and the EM algorithm was utilized to estimate the parameters. The scheduling variable was assumed to be known precisely. This paper mainly focuses on identifying a LPV model with an uncertain scheduling variable. The scheduling variable of a system is a slow-varying physical variable which determines the behavior of the system. The scheduling variables are often assumed to be measurable or to be estimated precisely in the existing literature, which is not often realizable. In this paper, it is assumed that the scheduling variables are measured with stochastic and unknown disturbances. A state-space model is utilized to describe the time-varying scheduling variables in this paper and an ARX model is used to represent the dynamics
of the nonlinear system. The main objectives in this paper are to estimate the coefficients of the parameters’ polynomial functions in the ARX model and the unknown parameters in weighting functions. The identification problems are solved by using EM algorithm which is an iterative optimization procedure essentially [7], [14], [23]. Since the noise-free scheduling variable is unknown, it should be estimated in identification process. Because of the nonlinearity of dynamics of the scheduling variable, a particle smoother based on SMC method is utilized to estimate the noise-free scheduling variable. The rest of this paper is given as follows: Section II describes the problem for identifying the LPV model with an uncertain scheduling variable and the necessary formulations are given. Section III shows the derivation procedure of the proposed identification approach in the EM algorithm framework. Section IV gives a numerial example and a mechanical unit to verify the proposed approach. Section V draws conclusions. II. P ROBLEM D ESCRIPTION A general LPV system can be described by the following autoregressive exogenous (ARX) model A(ωk , z −1 )yk = B(ωk , z −1 )uk + Δk ,
(1)
and A(ωk , z −1 ) = 1 + a1 (ωk )z −1 + · · · + ana (ωk )z −na B(ωk , z −1 ) = b1 (ωk )z −1 + · · · + bnb (ωk )z −nb
(2)
where yk is the output data of the system at each sampling moment and uk is the input data; Δk is the system noise and is often assumed to be the Gaussian white noise with variance σ 2 ; z −1 is the backward operator; The na and nb denote the model orders which are assumed to be known; ωk is the uncertain scheduling variable of the system and its dynamics is represented by a nonlinear state-space model in this paper. The coefficients in the polynomials A(ωk , z −1 ) and B(ωk , z −1 ) are linear combinations of the meromorphic functions {Γm (ωk )}m=1,...,nα and {Λn (ωk )}n=1,...,nβ of the scheduling variable ωk with static dependence ai (ωk ) = ai,0 +
nα
ai,m Γm (ωk ), i = 1, . . . , na
m=1 nβ
bj (ωk ) = bj,0 +
bj,n Λn (ωk ), j = 1, . . . , nb
(3)
n=1
Based on Eqs. (2) and (3), the system model in Eq. (1) is rewritten as yk = ϕk (ωk )T θ + Δk
where ϕk (ωk ) = [−yk−1 − Γ1 (ωk )yk−1 · · · − Γnα (ωk )yk−1 − yk−2 − Γ1 (ωk )yk−2 · · · − Γnα (ωk )yk−na uk−1 Λ1 (ωk )uk−1 . . . Λnβ (ωk )uk−1 uk−2 Λ1 (ωk )uk−2 . . . Λnβ (ωk )uk−nb ]T θ = [a10 a11 . . . a1nα a20 a21 . . . annaα
b10 b11 . . . b1nβ b20 b21 . . . bnnbβ ]T .
(5)
The following state-space model is utilized to represent the dynamics of ωk Hk = f (Hk−1 , uk−1 , γk−1 ) zk = g(Hk , uk , νk )
(6a) (6b)
ωk = ψ(Hk ),
(6c)
where zk is the observable process variable; Hk and uk are time-varying variables which denote the state and the input of the scheduling variable state-space model; γk and νk are zero mean Gaussian noise with variance Rγ and Rν , respectively; f (·) and g(·) are nonlinear state equation and output equation; ψ(·) generates the scheduling variable data. In this paper, the input-output data {uk , yk }k=1,2...N of the LPV system and the input-output data {uk , zk }k=1,2...N of scheduling variable model are recorded. Then, the observed data set is defined as Cobs = {y1:N , u1:N , u1:N , z1:N } and the parameters need to be estimated is Θ = {θ, σ}. At each sampling time k, ωk is hidden and unknown, then the missing data set can be expressed as Cmis = {ω0:N }. The complete data set is then represented as {Cobs , Cmis }. Therefore, the main objective of this paper is to estimate the parameter Θ given Cobs and Eqs. (6a)-(6c). III. EM ALGORITHM BASED APPROACH FOR LPV SYSTEM WITH NOISY SCHEDULING VARIABLE
A. Brief review on EM algorithm Maximum likelihood estimation is a general method which is often used in the area of industrial process modeling. But it becomes complex or even unreliable to compute the likelihood function of incomplete data set directly. The EM algorithm is a method to solve this kind of problems and is useful for the problems that the data are missing or involve hidden variables [24]. The EM algorithm consists of the expectation step (Estep) and the maximization step (M-step). First, the likelihood function of the complete data is given and its expectation, referred to as Q-function, is computed with respect to Cmis given the current estimate Θ and observed data Cobs in the E-step [25], [26]. Second, the parameter Θ is updated by maximizing the Q-function. Finally, the two steps iterate until the algorithm converges [25], [26]. The EM algorithm can be performed following the steps: 1) Set initial value for unknown parameter Θ . 2) Compute Q-function based on Cobs and current parameter Θ ,
(4)
Q(Θ|Θ ) = ECmis |(Cobs ,Θ ) {log pΘ (Cobs , Cmis )}. (7)
3) Maximize the Q-function to get the new parameter estimates, Θ = arg max Q(Θ|Θ ), Θ
The value of the likelihood function increases after each iteration of the EM algorithm and the optimal parameter estimate is obtained upon convergence. B. Mathematical derivations of the proposed method Before applying the EM algorithm, the likelihood function of complete data pΘ (Cobs , Cmis ) should be computed firstly. Based on the probability chain rule, it can be written as
The first term pΘ (y1:N |u1:N , z1:N , u1:N , ω0:N ) can be further decomposed into:
k=1
(10)
k=1
Based on Eqs. (4) and (5), yk depends only on the previous input and output sequences u1:k−1 = {u1 , . . . , uk−1 } and y1:k−1 = {y1 , . . . , yk−1 }, the noisy variable ωk , and the parameter Θ. Therefore, Eq. (10) is transformed into pΘ (y1:N |u1:N , z1:N , u1:N , ω0:N ) =
N
pΘ (yk |y1:k−1 , u1:k−1 , ωk ).
(11)
The p(ωk |Cobs ) can be simplified into p(ωk |z1:N ), which is a smoothing density function for ωk given the whole observations z1:N . ωk is a continuous random variable which makes it difficult to calculate the integral term in Eq. (17) directly. In this paper, the numercial solution, Sequential Monte-Carlo (SMC) method, is used to solve the intergral problem and the distribution of ωk is approximated by a set of random weighted samples. Sequential importance sampling and resampling steps are important for the SMC method to make sure the diversity of the samples is enough [7], [14]. For example, if s is a random variable and p(s) is a density function of s, then p(s) ≈
Similarly, in Eq. (9b), the second term can be rewritten as: pΘ (ω0:N |u1:N , z1:N , u1:N ) =
pΘ (ωk |ω1:k−1 , u1:N , z1:N , u1:N ).
×
(12)
k=1
From Eqs. (6a)-(6c), it can be seen that Eq. (12) and the third term in Eq. (9b) are independent of Θ. Define C1 = C2 =
pΘ (ω0:N |u1:N , z1:N , u1:N ) pΘ (u1:N , z1:N , u1:N ),
(17)
where C0 = p(ωk |Cobs ) log C3 dωk is independent of Θ. In order to further simplify the Q-function in Eq. (17), pΘ (yk |y1:k−1 , u1:k−1 , ωk ) and p(ωk |Cobs ) should be computed. Since Δk is the Gaussian white noise, then pΘ (yk |y1:k−1 , u1:k−1 , ωk ) is a Gaussian distribution
k=1
pΘ (ω0 |u1:N , z1:N , u1:N ) N
× p(ωk |Cobs )}dωk + C0 ,
pΘ (yk |y1:k−1 , u1:k−1 , ωk ) 1 −1 T T T =√ (yk − ϕk (ωk )θ) (yk − ϕk (ωk )θ) . exp 2σ 2 2πσ (18)
× pΘ (y1:N −1 |u1:N , z1:N , u1:N , ω0:N )
pΘ (yk |y1:k−1 , u1:N , z1:N , u1:N , ω0:N ).
(16)
Compute the expectation with respect to uncertain scheduling variable ωk yields N Q(Θ|Θ ) = {log pΘ (yk |y1:k−1 , u1:k−1 , ωk )
= pΘ (yN , y1:N −1 |u1:N , z1:N , u1:N , ω0:N ) = pΘ (yN |y1:N −1 , u1:N , z1:N , u1:N , ω0:N ) N
(15)
k=1
+ log C3 .
pΘ (y1:N |u1:N , z1:N , u1:N , ω0:N )
=
pΘ (yk |y1:k−1 , u1:k−1 , ωk )C3 .
Perform the E-step in Eq. (7) yields N pΘ (yk |y1:k−1 , u1:k−1 , ωk ) Q(Θ|Θ ) = EpΘ (ω0:N |Cobs ) log
(9b)
N k=1
(9a)
= pΘ (y1:N |u1:N , z1:N , u1:N , ω0:N ) × pΘ (ω0:N |u1:N , z1:N , u1:N ) × pΘ (u1:N , z1:N , u1:N ).
pΘ (Cobs , Cmis ) =
(8)
and then set Θ = Θ. 4) Repeat 2) and 3) until convergence.
pΘ (Cobs , Cmis ) = pΘ (y1:N , u1:N , z1:N , u1:N , ω0:N )
The Eq. (9b) is finally rewritten as
(13) (14)
M 1 δ(s − sm ), M m=1
(19)
where M is the total number of samples, δ denotes the dirac function, sm is drawn from the density function p(s) and the weight of each sample is 1/M . λ(s) is assumed as a function of s, then λ(s)p(s)dt can be approximately calculated as M 1 λ(s)p(s)ds ≈ λ(sm ). (20) M m=1 Similarly, the smoother density function p(ωk |z1:N ) can be numerically approximated as
where both C1 and C2 are independent of Θ and C3 = C1 C2 .
1 δ(ωk − ωkl ), L L
p(ωk |z1:N ) ≈
l=1
(21)
where ωkl is drawn from distribution p(ωk |z1:N ). The scheduling variable dynamic is described by the nonlinear state space model and the particle smoother based on sequential Monte Carlo method is used to estimate the distribution p(ωk |z1:N ). The ditribution p(ωk |z1:N ) at each moment can be represented by the sampled particles with l their corresponding weights as {ωkl , Wk|N }. The wights are initialized to 1/L after resampling for all the particles. The detailed derivation procedure is shown in [14]. According to Eqs. (18), (21), and (20), the Q-function in Eq. (17) can be written as 1 log pΘ (yk |y1:k−1 , u1:k−1 , ωkl ) L N
Q(Θ|Θ ) =
L
k=1 l=1
+ C0 .
(22)
To perform the M-step of the EM algorithm, as shown in Eq. (8), derivatives are taken with respect to θ and σ over the Q-function. Taking the derivative of Eq. (22) with respect to θ and equating it to zero, ∂ 1 log pΘ (yk |y1:k−1 , u1:k−1 , ωkl ) = 0. ∂θ L N
L
(23)
k=1 l=1
Assume the scheduling variable dynamic is described by the nonlinear state space model: √ ωk = 0.01 ωk−1 + uk−1 + γk−1 (26a) 0.9ωk + ωk + ν k . (26b) zk = 1 + ωk The system input uk and the input of scheduling variable model uk are chosen as Gaussian signal uk = N (0, 0.01) and periodic signal uk = 0.5 sin(0.05πk) + 0.7, respectively. γk and νk are Gaussian white noise represented by N (0, Rγ ) and N (0, Rν ), respectively. In this example, Rγ = 0.01 and Rν = 0.2. A Gaussian white noise is added to the output data and the signal-to-noise ratio (SNR) is calculated. Three Monte Carlo simulations are performed under SNRs 15dB, 20dB, and 25dB and 100 different noise sequences are generated in each Monte Carlo simulation. The mean value and standard derivation of estimated parameters are calculated and results are shown in Fig. 1 in which the ∗ and red bar denote the mean value and the standard derivation of estimated parameters, respectively. For all the estimated parameters, it is clearly shown that the mean values get closer to the true values and the standard derivations decrease with the increase of SNR.
Solving Eq. (23), the formula to update the parameter θ is N L ϕk (ωkl )yk l=1 θ = N k=1 . (24) L ϕk (ωkl )ϕTk (ωkl )
B. A mass-spring-damper system
Similarily, to estimate the unknown noise variance σ, taking the derivative of Eq. (22) with respect to σ 2 and equating it to zero N L (yk − ϕTk (ωkl )θ)T (yk − ϕTk (ωkl )θ) 2 k=1 l=1 σ = . N L (1) k=1 l=1 (25)
d d d (m x(t)) = F (t) − ks (t) − cd (t) x(t) (27) dt dt dt where x(t) denotes the position of the mass, ks (t) is the stiffness of the spring, cd (t) represents the varying damping, F (t) is the external force imposed on mass m and all of them are continuous time variables. Furthermore, F (t) and x(t) are regarded as the input and output variables in this example. Since the change of the damping excites the dynamics of the system, the variable cd (t) is selected as the scheduling variable of the system. The stiffness of the spring and the weight of the mass are 0.85N/m and 0.01kg, respectively. The input is selected as the random binary signal and the model of scheduling variable dynamic is given in Eq. (26). The input-output data of the system are shown in Fig. 3 and the × denotes the noise-free system output. The proposed method is applied to identify a LPV ARX model for the mass-spring-damper system. The simulated output of the identified model is compared with the true output of the system and the result is presented in Fig. 4.
k=1
l=1
IV. S IMULATION V ERIFICATION A. A numerical simulation Consider the LPV system with an uncertain scheduling variable described by the following LPV model [29] A(ωk , z −1 )yk = B(ωk , z −1 )uk + Δk , where A(ωk , z −1 ) = 1 + a1 (ωk )z −1 + a2 (ωk )z −2 a1 (ωk ) = 1 − 0.5ωk − 0.1ωk2
a2 (ωk ) = 0.5 − 0.7ωk − 0.1ωk2 , and B(ωk , z −1 ) = b1 (ωk )z −1 + b2 (ωk )z −2 b1 (ωk ) = 0.5 − 0.4ωk + 0.01ωk2 b2 (ωk ) = 0.5 − 0.3ωk −
In this example, a mass-spring-damper system, as shown in Fig. 2, is utilized to verify the proposed approach [27]. The mass is connected to a spring and a varying damper. This system can be decribed as
For evaluating the accuracy of the identified model, the relative error between the output of the system and output of the identified model is calculated following the equation in [28] Err =
0.02ωk2 .
var(y − y ) × 100%, var(y)
(28)
0.6 0.4 0.2 0 í0.2 í0.4 í0.6 í0.8 í1
The estimated means and standard derivations
The estimated means and standard derivations
The estimated means and standard derivations
1 0.8
1 0.8 0.6 0.4 0.2 0 í0.2 í0.4 í0.6 í0.8
The estimated parameters
í1
The estimated parameters
(a) SNR=15dB
1 0.8 0.6 0.4 0.2 0 í0.2 í0.4 í0.6 í0.8 í1
(b) SNR=20dB
The estimated parameters
(c) SNR=25dB
Fig. 1: The means and standard derivations of the estimated parameters calculated during Monte Carlo simulations under different SNRs
0.6
Real Output Simulated Output
0.4
Position/m
0.2 0 í0.2 í0.4 í0.6 í0.8 0
500
1000
1500 Time schedule/s
2000
2500
3000
Fig. 4: The real output and the identified model simulated output where y is the real system output, y is the simulated model output, var(·) denotes the variance of the signal. In this example, the relative error is 9.08% which proves that the identified model captures the dynamic of the mass-springdamper system. V. C ONCLUSION
Fig. 2: The mass-spring-damper system
Position/m
1
Output data
0.5 0 í0.5 í1
Force/N
0.5
Input data
0
í0.5 0
500
1000
1500 Time/s
2000
2500
3000
Fig. 3: The input and output of the mass-spring-damper system
This paper considers the identification problem of linear parameter varying (LPV) system with an uncertain scheduling variable in the EM algorithm scheme. A LPV ARX model with an uncertain scheduling variable is used to describe the relationship between the input and output of the system. Since the nonlinear state space model is used to represent the dynamics of the hidden scheduling variable, a particle smoother based on the sequential Monte-Carlo method is utilied to estimate probability distribution of the scheduling variable. This is a key solution to calculate the Q-function in the E-step of the EM algorithm. The algorithm to estimate the unknown model parameters and noise variance iteratively is derived. A simulation example and a mass-spring-damper system are used to validate the proposed approach and the results demonstrate that the proposed method can give satisfactory performance for identifying a LPV system with a noisy scheduling variable.
R EFERENCES [1] X. Yang, Y. Lu and Z. Yan, Robust global identification of linear parameter varying systems with generalised expectation-maximisation algorithm, IET Control Theory & Applications, vol. 9, no. 7, pp. 11031110, 2015. [2] W. Tang, J. Wu, Z. Shi, Identification of reduced-order model for an aeroelastic system from flutter test data, Chinese Journal of Aeronautics, vol. 30, no. 1, pp. 337-347, 2017. [3] S. Yin, X. Li, H. Gao and O. Kaynak, Data-based techniques focused on modern industry: An overview, IEEE Transactions on Industrial Electronics, vol. 62, no. 1, pp. 657-667, 2014. [4] G. Sun, Z. Zhu, Fractional-Order Dynamics and Control of RigidFlexible Coupling Space Structures, Journal of Guidance Control and Dynamics, vol. 38, no. 7, pp. 1324-1333, 2015. [5] G. Sun, Z. Zhu, Fractional-Order Tension Control Law for Deployment of Space Tether System, Journal of Guidance Control and Dynamics, vol. 37, no. 6, pp. 157-167, 2014. [6] Z. Fei, S. Shi, C. Zhao and L. Wu, Asynchronous Control for 2-D Switched Systems with Mode-Dependent Average Dwell Time, Automatica, vol. 79, pp. 198-206, 2017. [7] X. Jin, B. Huang, D.S. Shook, Multiple model LPV approach to nonlinear process identification with EM algorithm, Journal of Process Control, vol. 21, no. 1, pp. 182-193, 2011. [8] G. Sun, S. Xu, Z. Li, Finite-time fuzzy sampled-data control for nonlinear flexible spacecraft with stochastic actuator failures, IEEE Transactions on Industrial Electronics, DOI:10.1109/TIE.2017.2652366, 2017. [9] G. Sun, S. Xu, Z. Li, Finite-time fuzzy sampled-data control for nonlinear flexible spacecraft with stochastic actuator failures, IEEE Transactions on Industrial Electronics, vol. 64, no. 5, pp. 3851-3861, 2017. [10] R.D. Nowak, Nonlinear system identification, Circuits, System, and Signal Processing, vol. 21, no.1, pp. 109-122, 2002. [11] X.Hong, R.J. Mitchell, S.Chen, C.J. Harris, K.Li, G.W. Irwin, Model selection approaches for nonlinear system identification: a review, International Journal of Systems Science, vol. 39, no. 10, pp. 925-946, 2008. [12] L. Ljung, Preciction error estimation methods, Circuits Systems Signal Processing, vol. 21, no.1, pp. 11-21, 2002. [13] D.W. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, Journal of the Society for Industrial and Applied Mathematics, vol. 11, no. 2, pp. 431-441, 1963. [14] L. Chen, A. Tulsyan, B. Huang, F. Liu, Multiple model approach to nonlinear system identification with an uncertain scheduling variable using EM algorithm, Journal of Process Control, vol. 23, no. 10, pp. 1480-1496, 2013. [15] J.J. Hopfield, Neural networks and physical systems with emergent collective computational abilities, Proceedings of the National Academy of Sciences of the United States of Ameriva, vol. 79, pp. 2254-2258, 1982. [16] A.J. Smola, B. Sch¨olkopf, A tutorial on support vector regression, Statistics and Computing, vol. 14, no. 3, pp. 199-222, 2004. [17] Z. Xu, J. Zhao, J. Qian, Y. Zhu, Nonlinear MPC using identified LPV model, Industrial & Engineering Chemistry Research, vol. 2, no. 48, pp. 3043-3051, 2009. [18] R. Murray-Smith, T.A. Johansen, Multiple Model Approaches to Modelling and Control, Taylor & Francis, 1997. [19] L.Lee, K. Poolla, Identification of linear parameter-varying systems using nonlinear programming, In Proceedings of the 35th IEEE Conference on Decision and Control, pp. 1545-1550, 1996. [20] B. Bamieh, L. Giarr, Identification of linear parameter varying models, International Journal of Robust & Nonlinear Control, vol. 12, no. 9, pp. 841-853, 2002. [21] A. Banerjee, Y. Arkun, B. Ogunnaike, R. Pearson, Estimation of nonlinear multiple systems using linear models, AIChE Journal, vol. 43, no. 5, pp. 1204-1226, 1997. [22] Y. Zhu, Z. Xu, A method of LPV model identification for control, In Proceedings of the 17th World Congress of the International Federation of Automatic Control, pp. 5018-5023, 2008. [23] Y. Zhao, B. Huang, H. Su, J. Chu, Prediction error method for identification of LPV models, Journal of Process Control, vol. 22, no. 1, pp. 180-193, 2012.
[24] R.H. Shumway, D.S. Stoffer, An approach to time series smoothing and forecasting using the EM algorithm, Journal of Time Series Analysis, vol. 3, no. 4, pp. 253-264, 1982. [25] R.B. Gopaluni, A particle filter approach to identification of nonlinear processes under missing observations, The Canadian Journal of Chemical Engineering, vol. 86, no. 6, pp. 1081-1092, 2008. [26] M. Klaas, M. Briers, N. Freitas, A. de Doucet, S. Maskell, D. Lang, Fast prticle smoothing: if I had a million particles, In Proceedings of the 23rd International Conference on Machine Learning, Pittsburgh, 2006. [27] R. Toth, V. Laurain, M. Gilson, H. Garnier, Instrumental variable scheme for closed-loop LPV model identification, Automatica, vol. 48, no. 9, pp. 2314-2320, 2012. [28] Y. Zhu, H. Telkamp, J. Wang, Q. Fu, System identification using slow and irregular output samples, Journal of Process Control, vol. 19, no. 1, pp. 58-67, 2009. [29] V. Laurain, M. Gilson, R. Toth, H. Garnier, Refined instrumental variable methods for identification of LPV box-jenkins models, Automatica, vol. 46, no. 6, pp. 959-967, 2010.
