Multiple Model LPV Approach to Identification of Nonlinear Dual-rate System with Random Time Delay*

17th IFAC IFAC Symposium Symposium on on System System Identification Identification 17th 17th Symposium on Identification Beijing International Convention 17th IFAC IFAC Symposium on System SystemCenter Identification Beijing International Convention Center 17th IFAC Symposium on System Identification Beijing International Center October 19-21, 2015. 2015. Convention Beijing, China China Available online at www.sciencedirect.com Beijing International Convention Center October 19-21, Beijing, Beijing International Convention Center October 19-21, 2015. Beijing, China October October 19-21, 19-21, 2015. 2015. Beijing, Beijing, China China

ScienceDirect

IFAC-PapersOnLine 48-28 (2015) 1220–1225

Multiple Model LPV Approach to Multiple LPV to Multiple Model Model LPV Approach Approach to Identification of Nonlinear Dual-rate Identification of Nonlinear Dual-rate Identification of Nonlinear Dual-rate ⋆⋆ System with Random Time Delay System with Random Time Delay System with Random Time Delay ⋆

∗ ∗∗ ∗ Lei Yongsheng Lei Chen Chen ∗∗ Biao Biao Huang Huang ∗∗ Yongsheng Ding Ding ∗∗ ∗∗ ∗ ∗∗ ∗ Lei Chen Biao Huang Yongsheng Ding Lei Chen ∗ Biao Huang ∗∗ Yongsheng Ding ∗ ∗ ∗ Engineering Research Center of Digitized Textile & Fashion Engineering Research Center of Digitized Textile & Fashion ∗ ∗ Research Center Digitized Textile ∗ Engineering Technology, Ministry of Education; Education; College of Information Information Sciences Technology, Ministry of of Sciences Engineering Research Center of ofCollege Digitized Textile & & Fashion Fashion Technology, Ministry of Education; College of Information Sciences and Technology, Donghua University, Shanghai 201620, China and Technology, Donghua University,College Shanghai 201620, P. P. R. R. China Technology, Ministry of Education; of Information Sciences and Technology, Donghua University, Shanghai 201620, P. R. China (e-mail: [email protected], [email protected]). (e-mail:Donghua [email protected], [email protected]). and Technology, University, Shanghai 201620, P. R. China ∗∗ (e-mail: [email protected], [email protected]). ∗∗ Department of Chemical and Materials Engineering, University of Department of Chemical and Materials Engineering, University of (e-mail: [email protected], [email protected]). ∗∗ ∗∗ Department of Chemical and Materials Engineering, University of ∗∗ Alberta, Edmonton, ABMaterials T6G 2G6, 2G6,Engineering, Canada (e-mail: (e-mail: Alberta, AB T6G Canada Department of Edmonton, Chemical and University of Alberta, AB [email protected]) [email protected]) Alberta, Edmonton, Edmonton, AB T6G T6G 2G6, 2G6, Canada Canada (e-mail: (e-mail: [email protected]) [email protected])

Abstract: Abstract: Abstract: This paper This paper deals deals with with the the problem problem of of nonlinear nonlinear dual-rate dual-rate system system identification identification with with random random Abstract: This paper deals with the problem of nonlinear dual-rate system identification time delay. The proposed approach adopts the multiple modeling framework, and therandom global time delay. The proposed approach adopts the multiple modeling framework, and the global This paper deals with the problem of nonlinear dual-rate system identification with with random time proposed adopts multiple modeling the LPV model is by of models weighted aa probability LPV model The is represented represented by a a combination combination of various various local modelsframework, weighted by byand probability time delay. delay. The proposed approach approach adopts the the multiplelocal modeling framework, and the global global LPV represented by various models aa process probability function. The structure of model is aa state space form, has function. Theis considered structure of the the local localof model is in inlocal state spaceweighted form, and andby the process has LPV model model isconsidered represented by a a combination combination of various local models weighted bythe probability function. The considered structure of the local model is in a state space form, and the process fast rate inputs and slow rate outputs with random time delay. The expectation maximization fast rate inputs and slow rate outputs with random time delay. The expectation maximization function. The considered structure of the local model is in a state space form, and the process has has fast and rate time The maximization algorithm is utilized utilized to formulate formulate andwith solverandom the problem problem of interest. interest. The parameters parameters of the the algorithm is to and solve the of The of fast rate rate inputs inputs and slow slow rate outputs outputs with random time delay. delay. The expectation expectation maximization algorithm is utilized to and the of parameters of local models the functions are simultaneously. smoothing local models and the weighting weighting functions are estimated estimated simultaneously. The particle particle smoothing algorithm is and utilized to formulate formulate and solve solve the problem problem of interest. interest. The parameters of the the local and are simultaneously. The particle smoothing technique is to the of functions. effectiveness of technique is adopted adopted to handle handle functions the computation computation of expectation expectation functions. The effectiveness of local models models and the the weighting weighting functions are estimated estimated simultaneously. TheThe particle smoothing technique is adopted to handle the computation of expectation functions. The effectiveness of the proposed approach is further illustrated through a simulation example. the proposed approach is further illustrated through a simulation example. technique is adopted to handle the computation of expectation functions. The effectiveness of the proposed approach illustrated through aa simulation example. the proposed approach is is further further illustrated through simulation example. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Multiple models, Expectation maximization algorithm, Keywords: Multiple models, Expectation maximization algorithm, Nonlinear Nonlinear process, process, System System Keywords: Multiple models, Expectation maximization algorithm, Nonlinear identification, Time delay identification, Time delay Keywords: Multiple models, Expectation maximization algorithm, Nonlinear process, process, System System identification, identification, Time Time delay delay 1. scheduling 1. INTRODUCTION INTRODUCTION scheduling variables. variables. Jin Jin et et al. al. (2011) (2011) proposed proposed aa nonlinear nonlinear 1. INTRODUCTION scheduling variables. Jin et al. (2011) proposed aa nonlinear system identification approach under the system identification approach under the expectationexpectation1. INTRODUCTION scheduling variables. Jin et al. (2011) proposed nonlinear approach under maximization (EM) framework, framework, wherein, the expectationlocal model model Most maximization (EM) the local Most of of the the system system identification identification approaches approaches have have focused focused system system identification identification approach wherein, under the the expectationmaximization (EM) framework, wherein, the local model Most of the system identification approaches have focused takes the form of the auto regressive exogenous (ARX) on systems whose input-output data are available at every takes the form of the auto regressive exogenous (ARX) on systems whose input-output data are available at every maximization (EM) framework, wherein, the local model Most of the system identification approaches have focused takes the form of the auto regressive exogenous (ARX) on systems whose input-output data are available at every model. Deng and Huang (2012) described an identification sampling time (S¨ o derstr¨ o m and Stoica, 1988; Ljung, 1999). model.the Deng andofHuang (2012) described an identification sampling time (S¨ oderstr¨ om and Stoica, Ljung, form the auto regressive exogenous (ARX) on systems whose input-output data are1988; available at1999). every takes model. Deng Huang (2012) identification sampling time (S¨ o o for nonlinear varying under Stoica, Ljung, 1999). For multirate these algorithms cannot be approach for and nonlinear parameter varyingan systems under For multirate systems, these algorithms cannot be directly directly model. Deng and Huangparameter (2012) described described ansystems identification sampling timesystems, (S¨ oderstr¨ derstr¨ om m and and Stoica, 1988; 1988; Ljung, 1999). approach approach for nonlinear parameter varying systems under For multirate systems, these algorithms cannot be directly the EM framework, and the state-space model is used as applied. A system with two different sampling periods the EM framework, and the state-space model is used as applied. A system with two different sampling periods For multirate systems, these algorithms cannot be directly approach for nonlinear parameter varying systems under the EM framework, and the state-space model is used as applied. A system with two different sampling periods local model. They considered the process identification between input and output is called dual-rate system. local model. They considered the process identification between input and output is called dual-rate system. the EM framework, and the state-space model is used as applied. A system with two different sampling periods the local model. They considered the process identification between input and without time delay problems. output is called dual-rate system. Over the past few decades, techniques for linear system without delay problems. Over the input past few for linear system. system the local time model. They considered the process identification between anddecades, output techniques is called dual-rate without time delay Over few techniques linear system identification been and linear dual-rate identification have been matured, matured, and for linear dual-rate time delay problems. problems. Over the the past past have few decades, decades, techniques for linear system without Chen et al. (2014b) considered Chen et al. (2014b) considered identification identification of of nonlinear nonlinear identification have been matured, and linear dual-rate system identification has also also been discussed discussed extensively system identification has been identification have been matured, and linearextensively dual-rate Chen et al. (2014b) considered identification of nonlinear systems with a single noisy scheduling variable, and systems with a single noisy scheduling variable, and the the Chen et al. (2014b) considered identification of nonlinear system identification has also been discussed extensively in literatures (Ding et ethas al.,also 2008; Ding and Chen, Chen, 2004). systems in literatures (Ding al., 2008; Ding and 2004). system identification been discussed extensively with a single noisy scheduling variable, and the measurement of the system has an unknown time delay. measurement of the system has an unknown time delay. systems with a single noisy scheduling variable, and the in literatures (Ding et al., 2008; Ding and Chen, 2004). Time delay is a long standing problem in many industrial Time delay is (Ding a long et standing problem many industrial in literatures al., 2008; Ding in and Chen, 2004). measurement of the system has an unknown time delay. The global global model model is aasystem combination ofunknown multiple time local delay. ARX The is combination of multiple local ARX measurement of the has an Time delay is a long standing problem in many industrial processes, and methods for identification of linear time processes, for identification of linear time The Time delayand is a methods long standing problem in many industrial global model is multiple local ARX models. The parameters and models. The parameters and the the of unknown time delay global model is aa combination combination ofunknown multiple time local delay ARX processes, and methods of linear delay systems have been for wellidentification developed (Bj¨ (Bj¨ orklund rklund and The delay systems been well developed and processes, and have methods for identification ofo linear time time models. The parameters and the unknown time delay are estimated simultaneously under the EM algorithm. are estimated simultaneously under the EM algorithm. models. The parameters and the unknown time delay delay systems have been well developed (Bj¨ o rklund and Ljung, 2003; Richard, Richard, 2003). Chen et et al. al. (2014a) have Ljung, 2003; 2003). Chen (2014a) delay systems have been well developed (Bj¨ orklund have and are simultaneously the algorithm. When the time time delay is is time timeunder varying, the mentioned mentioned When the delay varying, the are estimated estimated simultaneously under the EM EM algorithm. Ljung, 2003). Chen al. have considered the parameter estimation problem for considered theRichard, parameter estimation problem for aa dualdualLjung, 2003; 2003; Richard, 2003). Chen et et al. (2014a) (2014a) have When the time delay is time varying, the mentioned approaches are not not applicable. This paper is mentioned concerned approaches are applicable. This paper is concerned When the time delay is time varying, the considered the parameter estimation problem for a dualrate system with constant time delay. rate systemthe with constant estimation time delay. problem for a dual- approaches considered parameter are This paper concerned with of with with the the identification identification of aa nonlinear nonlinear dual-rate system with approaches are not not applicable. applicable. Thisdual-rate paper is issystem concerned rate system with constant time delay. rate system with constant time delay. with the identification of a nonlinear dual-rate system unknown and randomly varying time delay. Industrial processes are typically operated along certain unknown and randomly timedual-rate delay. system with Industrial processes are typically operated along certain with the identification of varying a nonlinear with unknown and randomly varying time delay. Industrial processes operated along certain fixed trajectories, and the trajectories are composed fixed trajectories, andare thetypically trajectories are often often composed and randomly varying time delay. Industrial processes are typically operated along certain unknown In this paper, a multiple model LPV approach is develthis paper, a multiple model LPV approach is develfixed trajectories, and are composed of several pre-determined operating points, points, governed by In of several pre-determined operating governed by fixed trajectories, and the the trajectories trajectories are often often composed In this paper, aa multiple LPV approach oped, the LPV represented comoped, and the global global LPV model is represented byis adevelcomIn thisand paper, multiple model is LPV approachby is a develof several pre-determined operating points, governed by of several pre-determined operating points, governed by oped, and the global LPV model is represented by a bination of various local models weighted by a probability ⋆ bination of various local models weighted by a probability oped, and the global LPV model is represented by a comcom⋆ This This work work was was supported supported by by the the Natural Natural Sciences Sciences and and EngineerEngineerbination of local models by ⋆ density function. function. The considered structure ofprobability the local local ⋆ density structure the work supported by and bination of various variousThe localconsidered models weighted weighted by aaof probability ing Research Council (NSERC) of Canada, the ReThis work was was supported by the the Natural Sciences and EngineerEngineer⋆ This ing Research Council (NSERC) of Natural Canada,Sciences the Fundamental Fundamental ReThis work was supported by the Natural Sciences and Engineerdensity function. The considered structure of the local model is aa state-space model. The EM algorithm is utilized ing Research Council (NSERC) of Canada, the Fundamental Remodel is state-space model. The EM algorithm is utilized search Funds for for the Central Central Universities, Shanghai Pujiang Program Program density function. The considered structure of the local ing Research Council (NSERC) of Canada, the Fundamental Research Funds the Universities, Shanghai Pujiang ing Research Council (NSERC) of Canada, the Fundamental Remodel is a state-space model. The EM algorithm is utilized to formulate and solve the problem of interest. The paramsearch Funds for the Central Universities, Shanghai Pujiang Program (15PJ1400100), Key Project of the National Nature Science to formulate and solve model. the problem of interest. The paramsearch Funds Funds for for the the Central Central Universities, Shanghai Pujiang Program model is a state-space The EM algorithm is utilized (15PJ1400100), Key Project of theShanghai NationalPujiang NatureProgram Science search Universities, to and the problem of The (15PJ1400100), the Key of the eters of the the local local models and the weighting weighting functions are Foundation of China China (No.Project 61134009), the National National Nature Nature Science (15PJ1400100), the Project of National Nature Science eters of models the are Foundation of (No. 61134009), the National Nature to formulate formulate and solve solve theand problem of interest. interest.functions The paramparam(15PJ1400100), the Key Key Project of the the National Nature Science eters of the local models and the weighting functions are Foundation of China (No. 61134009), the National Nature Science of China (No. 61473078), Specialized Research Fund for estimated simultaneously. The particle smoothing of China China (No. (No.61473078), 61134009),Specialized the National National NatureFund Science Foundation of Research for estimated Thetheparticle smoothing techeters of thesimultaneously. local models and weighting functionstechare Foundation of China (No. 61134009), the Nature Science Foundation of China (No. 61473078), Research Fund for Shanghai Leading Talents, Project of ofSpecialized the Shanghai Shanghai Committee of estimated simultaneously. The particle smoothing techFoundation of China (No. 61473078), Specialized Research Fund for Shanghai Leading Talents, Project the Committee of Foundation of China (No. 61473078), Specialized Research Fund for estimated simultaneously. The particle smoothing techShanghai Leading Talents, Project of the Shanghai Committee of Science and Technology (Nos. 13JC1407500). Shanghai Leading Talents, Project of Science and Technology (Nos. 13JC1407500). Shanghai Leading Talents, Project of the the Shanghai Shanghai Committee Committee of of Science and Technology (Nos. 13JC1407500). Science and Technology (Nos. 13JC1407500). Science and Technology (Nos. 13JC1407500).

Copyright © 2015, IFAC 1220 2405-8963 © IFAC (International Federation of Automatic Control) Copyright IFAC 2015 2015 1220Hosting by Elsevier Ltd. All rights reserved. Copyright © IFAC 2015 1220 Peer review under responsibility of International Federation of Automatic Copyright © IFAC 2015 1220 Copyright © IFAC 2015 1220Control. 10.1016/j.ifacol.2015.12.298

2015 IFAC SYSID October 19-21, 2015. Beijing, China

Lei Chen et al. / IFAC-PapersOnLine 48-28 (2015) 1220–1225

nique is adopted to handle the computation of expectation functions for the EM algorithm. The remainder of this paper is organized as follows: Section 2 presents the problem formulation. Section 3 solves the identification problem under the framework of the EM algorithm. Section 4 verifies the proposed method through a simulation example. Conclusions are given in Section 5. 2. PROBLEM FORMULATION Consider the following nonlinear dual-rate state-space model with random time delay: xk+1 = fθf (xk , uk , γk ),

(1a)

yt = gθg (xt−λt , vt ),

(1b)

where uk is the fast rate input, xk is the hidden state for k = 1, . . . , N ; yt is the slow rate output with random time delay λt for t = 1, . . . , H, where t = ϕk, and ϕ is a positive integer; f (·) and g(·) are the state and output mapping functions, and γk and vt are zero mean Gaussian noise represented by N (0, Rγ ) and N (0, Rv ), respectively. θf and θg are unknown parameters. T = [T1 , · · · , TM ] are M different operating points. The time delay λt is a random integer, which is bounded by λmin and λmax , and D = λmax − λmin + 1 is the range of the time delay. The random time delay λt is unknown but follows a Markov chain model, and it has the following transition probability: πd1 d1 πd1 d2  πd 2 d 1 π d 2 d 2 Π= ..  ... . π d D d 1 πd D d 2 

where

· · · πd 1 d D · · · πd 2 d D  , ..  .. . .  · · · πdD dD 

(2)

πdi dj = P r(λt = dj |λt−1 = di ).

(3)

The row elements of the transition probability matrix Π �D are summated to 1, that is j=1 πdi dj = 1; all the elements of the transition probability matrix Π are between 0 and 1, that is 0 ≤ pidi dj ≤ 1.

2

αk,i =

yt =

αt,i yt,i .

M �

i=1

,

(5)

2 k −Ti ) exp(− (w2(o ) 2 i)

Since the output is slow rate, we denote y1:N = {yobs , ymis }, where ymis = {ym1 , ym2 , . . . , yma } is the unobserved or treated as missing output data, hereinafter denoted as ym1 :ma , and yobs = {yo1 , yo2 , . . . , yob } is the observed output data, hereinafter denoted as yo1 :ob . To match the length of the unobserved output data ymis , 0 is used to extend the time delays λ1 , . . . , λH to λ1 , . . . , λN . At each sampling time k, Ik and xk are hidden, and thus the missing or hidden data set is Cmis = {x1:N , I1:N , ym1 :ma , λ1:N }. It is assumed that the model structure of Eq. (1) is known a priori, we can access the retarded data {u1:N , y1:H }, and the data w1:N , T1:M are available. Therefore, we have the observed data set Cobs = {yo1 :ob , u1:N , w1:N , T1:M }. Then the complete data set that include the missing data are {Cobs , Cmis }. The objective of this paper is using the observed data set Cobs to estimate the parameters Θ = {ΘIk , O}, where ΘIk = {θf Ik , θg Ik }. We adopt the EM algorithm to solve this problem. 3. EM ALGORITHM The EM algorithm is used for the maximum-likelihood estimation from incomplete data (Dempster et al., 1977). It consists of the expectation step (E-step) and the maximization step (M-step). In the E-step, we calculate the expectation of the log-likelihood function for the complete data set {Cobs , Cmis } with respect to the missing data Cmis based on the observed data Cobs and the current estimated parameter set Θ′ , also called the Q-function. In the M-step, through maximizing the Q-function, the parameters are re-estimated. The mathematical formulation of EM algorithm is as follows (McLachlan and Krishnan, 2008): (1) Initialization: Initialize the parameter to Θ′ . (2) E-step: Calculate the approximate Q-function as Q(Θ|Θ′ ) = ECmis |(Cobs ,Θ′ ) {log pΘ (Cobs , Cmis )}. (6) (3) M-step: Through maximizing the Q-function, the update parameter can be estimated as Θ = arg max Q(Θ|Θ′ ). (7) Θ

Then set Θ′ = Θ. (4) Iterate: Evaluate the relative change of the estimated parameters, Θ − Θ′ �. (8) Θ′ If δ is larger than a pre-determined tolerance, then repeat steps 2 and 3.

(4b)

δ=�

i=1

I1:N is introduced to represent the local model identity, and Ik can be inferred from wk and T1:M . A normalized exponential function is used to represent the weighting function:

k −Ti ) exp(− (w2(o ) 2 i)

where oi is the validity width of the Ik = ith local model which needs to be estimated from data. In Eq. (5), oi has a lower bound oi,min and a upper bound oi,max . O = {o1 , o2 , . . . , oM } is a set of validity width for the M local models respectively.

Within relatively small region of each of the operating points, a local model can be applied to approximate the process dynamics. Since the operating condition is time varying, multiple model approach is utilized, and the global nonlinear model is a weighted interpolation of local linear models, so that system state and output can be written as M � xk = αk,i xk,i , (4a) i=1 M �

1221

Under certain regularity conditions, the EM algorithm guarantees to converge to a stationary point (Wu, 1983).

1221

2015 IFAC SYSID 1222 October 19-21, 2015. Beijing, China

Lei Chen et al. / IFAC-PapersOnLine 48-28 (2015) 1220–1225

3.1 Formulation of the multiple model approach based on the EM algorithm For any generic sequence χk , let χi:j = {χi , χi+1 , . . . , χj }, with χi:j = 0 for i > j. Under the EM algorithm, the complete likelihood function pΘ (Cobs , Cmis ) needs to be estimated. pΘ (Cobs , Cmis ) can be decomposed using the probability chain rule as pΘ (Cobs , Cmis ) (9a) = pΘ (y1:N , u1:N , x1:N , w1:N , T1:M , I1:N , λ1:N ) = pΘ (y1:N |u1:N , x1:N , w1:N , T1:M , I1:N , λ1:N ) · P rΘ (x1:N |u1:N , w1:N , T1:M , I1:N , λ1:N ) · pΘ (I1:N |u1:N , w1:N , T1:M , λ1:N ) · pΘ (λ1:N |u1:N , w1:N , T1:M ) (9b) · pΘ (u1:N , w1:N , T1:M ). Since the input sequence u1:N is considered to be known and it does not play a role in the following derivations. Hereinafter, u1:N is omitted from the probability density functions for notation simplicity. Simplifying Eq. (9b), we get N  pΘ (Cobs , Cmis ) = pΘ (yk |xk−λk , λk , Ik )

(6) P rO (Ik = i|wk , T1:M ); (7) pΘ′ (Ik = i|Cobs ) ˆ k = d, Ik = i|Cobs ) can be simplified as pΘ′ (xk−λˆ k , λ ˆ k = d, Ik = i|yo :o ), which is a smoothing pΘ′ (xk−λˆ k , λ 1 b density function for the hidden state with the observed outputs. Since it is intractable to calculated Eq. (13) directly, the numerical solutions are utilized. In this work, the Sequential Monte-Carlo (SMC) method is applied (Metropolis and Ulam, 1949). The SMC methods is to approximate the distributions of interest by a set of random samples with associated weights. Following the procedure used in Gopaluni (2008) and Gopaluni (2010), ˆ k = d, Ik = the smoothing density function pΘ′ (xk−λˆ k , λ i|yo1 :ob ) can be numerically calculated as: L

 ˆ k = d, Ik = i|yo :o ) ≈ 1 pΘ′ (xk−λˆ k , λ δ(xk−d − xlk−d ), 1 b L l=1

(14)

xlk−d

represents the where δ(·) is a dirac-delta function, lth hidden state particle of the state-space model, and L is the number of particles.

k=1

As Gaussian white noise has been assumed for the stateˆ P rO (Ik |wk , T1:M ) · C1 , space model in Eq. (1), pΘ (yk |xk−λˆ k , λk = d, Ik = i), pΘ (yk |xk , Ik = i) and pΘ (xk |xk−1 , Ik = i) are Gaussian k=2 k=1 (10) probability density function. When yk is missing, we can l l where C1 = pΘ (λ1:N |u1:N , w1:N , T1:M )pΘ (u1:N , w1:N , T1:M ) obtain an estimate of yk for a given xk from the density function. P r (I = i|w , T ) and p O k k 1:M Θ′ (Ik = i|Cobs ) can is independent of Θ. be calculated using Eq. (5). By computing the E-step shown in Eq. (6), the Q-function can be calculated as Eq. (11a)(see the next page). Since Using the particle smoother and the probability chain rule, y1:N = {yo1 :ob , ym1 :ma }, and by further marginalizing the the Q-function can be determined as Eq. (15). states, the Q-function can be derived as Eq. (11b), where For the M-step, derivatives are taken with respect to Θ C2 = EpΘ′ (x1:N ,I1:N ,λ1:N ,ym1 :ma |Cobs ) log C1 . and O in the Q-function. By equating the derivatives to zero, the updated parameters are obtained. Since the time delay λk is a discrete random variable, the first term of Eq. (11b) can be simplified and written as Eq. (12)(This long equation is shown in the next page). 3.2 Summary of the proposed identification approach Since the time-varying delay λt is unknown, the possible combinations for λ1:H is DH , which creates a challenge to Multiple model LPV approach for nonlinear dual-rate sysevaluate pΘ′ (xk−λk , λk = d, Ik |Cobs ). With the knowledge tem identification with time delay using the EM algorithm of the transition probability for the time delay, we will can be summarized as: ˆ k for evaluation. Due to directly use point estimates of λ (1) Initialization: Set initial value to Θ′ . page limit, the detailed derivation is omitted here. (2) Particle smoother : According to the current statespace model, apply particle smoother to calculate the Since the local model identity Ik is a discrete random smoothed density function of the state in Eq. (14). variable, Eq. (11b) can be further expressed as Eq. (13). (3) E-step: Evaluate the approximate Q-function accordNow to compute the Q-function in Eq. (13), we need ing to Eq. (15) through particle smoothers using Θ′ . to evaluate the following probability mass and density (4) M-step: Maximize the Q-function with respect to Θ functions: to update the parameter, and then set Θ′ = Θ. (5) Iterate: Evaluate the relative change shown in Eq. ˆ k = d, Ik = i|Cobs ), when the observed (1) pΘ′ (xk−λˆ k , λ (8), and according to the pre-determined tolerance, outputs are available at time instant k, i.e. yk ∈ yobs decide whether to repeat step 2 to 4 or terminate. with yobs = yo1 :ob ; (2) pΘ′ (xk , Ik = i|Cobs ), when the observed outputs are unavailable at time instant k, i.e. yk ∈ ymis with 4. SIMULATION EXAMPLE ymis = ym1 :ma ; ˆ k = d, Ik = i), when yk ∈ yobs ; (3) pΘ (yk |xk−λˆ k , λ A first-order process with varying system parameters is (4) pΘ (yk |xk , Ik = i), when yk ∈ ymis ; described in Zhu and Xu (2008) as: (5) pΘ (xk |xk−1 , Ik = i); ·

N 

pΘ (xk |xk−1 , Ik ) · pΘ (x1 |I1 ) ·

N 

1222

2015 IFAC SYSID October 19-21, 2015. Beijing, China

Lei Chen et al. / IFAC-PapersOnLine 48-28 (2015) 1220–1225

′

Q(Θ|Θ ) = EpΘ′ (x1:N ,I1:N ,λ1:N ,ym1 :ma |Cobs ) N 

log pΘI1 (x1 |I1 ) +



N 

log pΘ (yk |xk−λk , λk , Ik ) +

k=1

ob  

+

log P rO (Ik |wk , T1:M ) + log C1

N  

k=2

+ +

log pΘ (xk |xk−1 , Ik )pΘ′ (xk , xk−1 , Ik |Cobs )dxk dxk−1 dIk

log P rO (Ik |wk , T1:M )pΘ′ (Ik |Cobs )dIk + C2 ,



log pΘ (yk |xk−λk , λk = d, Ik )pΘ′ (xk−λk , λk = d, Ik |Cobs )dxk−d dIk .

ob  M  

+

xk−d

ma  M  

k=m1 i=1

+

(12)

xk−d ,Ik

k=o1 i=1

+

(11b)

Ik

k=o1 d=λmin

+

log pΘ (yk |xk−λk , λk , Ik )pΘ′ (xk−λk , λk , Ik |Cobs )dxk−λk dλk dIk dyk

log pΘ (x1 |I1 )pΘ′ (x1 , I1 |Cobs )dx1 dI1

λ max

Q(Θ|Θ ) =

(11a)

xk ,xk−1 ,Ik

x1 ,I1 N  

′



xk−λk ,λk ,Ik ,yk



k=1 ob 

log pΘ (xk |xk−1 , Ik )+

log pΘ (yk |xk−λk , λk , Ik )pΘ′ (xk−λk , λk , Ik |Cobs )dxk−λk dλk dIk

k=o1 xk−λk ,λk ,Ik ma   k=m1

+

N 

k=2

k=1

=

1223

N  M  

k=2 i=1 M  

i=1 x1 M N  

ˆ k = d, Ik = i)pΘ′ (x ˆ , λ ˆ k = d, Ik = i|Cobs )dxk−d log pΘ (yk |xk−λˆ k , λ k−λk log pΘ (yk |xk , Ik = i)pΘ′ (xk , Ik = i|Cobs )dxk dyk

xk ,yk

log pΘ (xk |xk−1 , Ik = i)pΘ′ (xk , xk−1 , Ik = i|Cobs )dxk dxk−1 xk ,xk−1

log pΘ (x1 |I1 = i)pΘ′ (x1 , I1 = i|Cobs )dx1 log P rO (Ik = i|wk , T1:M )pΘ′ (Ik = i|Cobs ) + C2 .

(13)

k=1 i=1

Q(Θ|Θ′ ) ≈

ob  L M  1  ˆ k = d, Ik = i)P rΘ′ (Ik = i|Cobs ) log pΘ (yk |xlk−λˆ , λ k L i=1 k=o1

+

1 L

ma 

log pΘ (ykl |xlk , Ik = i)P rΘ′ (Ik = i|Cobs )

k=m1 i=1 l=1 N

+

l=1

M  L 

M

L

L

1  log pΘ (xjk |xsk−1 , Ik = i)P rΘ′ (Ik = i|Cobs ) L i=1 j=1 s=1 k=1

+

M N  

log P rO (Ik = i|wk , T1:M )P rΘ′ (Ik = i|Cobs ) + C2 .

(15)

k=1 i=1

G(s, w) =

K(w) , τ (w)s + 1

The process can be transferred as xk+1 = Axk + Buk + γk , x3 yt = C 2 t−λt + vt , xt−λt + 1

K(w) = 0.6 + w2 τ (w) = 3 + 0.5w3 , w ∈ [1, 4]. 1223

(16a) (16b)

2015 IFAC SYSID October 19-21, 2015. Beijing, China 1224

Lei Chen et al. / IFAC-PapersOnLine 48-28 (2015) 1220–1225

10

1 0.5

5

u

0 −0.5 −1

0 0

500

1000

1500

10

−5

5 0 True output Estimated output

y

−10

−5 −10 −15

0

500

1000

−15

1500

0

500

Time

1000

1500

Time

Fig. 1. The input-output data of the nonlinear process

Fig. 3. Comparison of the identified nonlinear model (selfvalidation)

4.5 1 4

0.9

3 2.5 2 1.5

0.7 0.6 0.5 0.4 0.3 0.2

1 0.5

Local model 1 Local model 2 Local model 3

0.8 The paobability of local models

Scheduling variable

3.5

0.1 0

500

1000

1500

0

Time

Fig. 2. The scheduling variable of the nonlinear process 1 τ (w) ,

1

1.5

2

2.5 3 Scheduling variable

3.5

4

Fig. 4. The probability of each local model (self-validation)

K(w) τ (w) ,

where A = 1 − B = C = 1. The operating points are T = [1, 2.25, 4]. Here we set t = 5k, that is, 80% outputs are missing. For simplicity, the time delay can change from 2 to 3 samples, and the transition probability  0.64 0.36 . of the time delays is 0.81 0.19

15

10

5

For the M-step shown in Eq. (7), differentiating Eq. (15) with respect to Θi and oi and then equating it to zero, we get updated parameters as shown in Eq. (17) and Eq. (18). The validity width cannot be derived in an explicit form when maximizing Eq. (15). The mathematical formulation of the optimization problem in searching for the optimal oi can be expressed by Eq. (19).

0

−5 True output Estimated output

−10

−15

0

100

200

300

400 500 Time

600

700

800

900

The simulated data are shown in Figures 1 and 2. The observed output has random time delay λ in each sample. Applying the proposed method to the simulated data, three local models are identified, where the convergence tolerance is 1e-5. The self-validation results are shown in Figures 3 and 4.

Fig. 5. Comparison of the identified nonlinear model (cross-validation)

Cross validation is conducted under different operating points at 1.5 and 4.5, respectively. The results are shown in figures 5 and 6.

5. CONCLUSION

The results of the example show that the predictions obtained by the identified nonlinear process models are

in close agreement with the true outputs. These results demonstrate the effectiveness of the proposed algorithm.

This article considers a multiple model LPV approach under the EM framework for nonlinear dual-rate system

1224

2015 IFAC SYSID October 19-21, 2015. Beijing, China

Lei Chen et al. / IFAC-PapersOnLine 48-28 (2015) 1220–1225

Ai =

L  N  L 

k=1 j=1 s=1

P rΘ′ (Ik = i|Cobs )(xjk xsk−1 − Bi uk−1 xsk−1 )

L  L N  

k=1 l=1 s=1

Bi =

L  L N  

1225

.

(17)

.

(18)

P rΘ′ (Ik = i|Cobs )xsk−1 xsk−1

P rΘ′ (Ik = i|Cobs )(xjk uk−1 − Ai xsk−1 uk−1 )

k=1 l=1 s=1 L L  N  

P rΘ′ (Ik = i|Cobs )uk−1 uk−1

k=1 l=1 s=1

max oi

M N  

log P rO (Ik = i|wk , T1:M )P rΘ′ (Ik = i|Cobs ). S.t. omin ≤ oi ≤ omax

15 True output Estimated output

Estimated output

10

5

0

−5

−10

−15 −15

−10

−5

0 True output

(19)

k=1 i=1

5

10

15

Fig. 6. The scatter plot comparisons of the actual outputs and estimated outputs (cross-validation) identification with unknown random time delay at each sample. The local models are identified at each process operating points, and the complete dynamics of the nonlinear system is a combination of all the local models associated with a normalized exponential function as its probability density function. The parameters of the models and the exponential functions are estimated simultaneously. The self and cross validation results of a simulation example demonstrate that the proposed method can give satisfactory results for identifying nonlinear dual-rate systems evolving along different operating points with random time delay at each sampling instant. REFERENCES S¨ oderstr¨ om, T., Stoica, P. System Identification. Englewood-Cliffs, NJ: Prentice-Hall, 1988. Ljung, L. System Identification: Theory for the User, 2nd ed. Engle-wood-Cliffs, NJ: Prentice-Hall, 1999. Ding, F., Liu, P.X., Yang, H. Parameter identification and intersample output estimation for dual-rate systems. IEEE Transactions on Systems, Man, and Cybernetics– Part A: Systems and Humans, 38(4), 966-975, 2008. Ding, F., Chen, T. Combined parameter and output estimation of dual-rate systems using an auxiliary model. Automatica, 40(10) 1739-1748, 2004.

Bj¨orklund, S., Ljung, L. A review of time-delay estimation techniques. In Proceedings of the 42nd IEEE Conference on Decision and Control. Maui, Hawaii USA. 3, 25022507, 2003. Richard, J. P. Time-delay systems: an overview of some recent advances and open problems. Automatica, 39(10), 1667-1694, 2003. Chen, L., Han, L., Huang, B., Liu, F. Parameter estimation for a dual-rate system with time delay. ISA Transactions: The Journal of Automation, 53(5), 13681376, 2014. Jin, X., Huang, B., Shook, D. S. Multiple model LPV approach to nonlinear process identification with EM algorithm. Journal of Process Control, 21(1), 182-193, 2011. Chen, L., Huang, B., Liu, F. Multi-model approach to nonlinear system identification with unknown time delay. 19th World Congress of the International Federation of Automatic Control, Cape Town, South Africa, Aug 2429, 2014. Deng, J., Huang, B. Identification of nonlinear parameter varying systems with missing output data. AIChE Journal. 58(11), 3454-3467, 2012. Dempster, A. P., Laird, N. M., Rubin, D. B. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, 39(1), 1-38, 1977. McLachlan, G. J., Krishnan, T. The EM algorithm and extensions, 2E. USA: John Wiley & Sons, Inc, 2008. Wu, J. On the convergence properties of the EM algorithm. The Annals of Statistics. 11(1), 95-103, 1983. Metropolis, N., Ulam, S. The Monte Carlo method. Journal of the American Statistical Association. 44(247), 335-341, 1949. Gopaluni, R B. Nonlinear system identification under missing observations: The case of unknown model structure. Journal of Process Control, 20(3), 314-324, 2010. Gopaluni, R B. A particle filter approach to identification of nonlinear processes under missing observations. The Canadian Journal of Chemical Engineering, 86(6), 10811092, 2008. Zhu, Y., Xu, Z. A method of LPV model identification for control. 17th IFAC World Congress, Seoul, Korea, July 6-11, 5018-5023, 2008.

1225