- Email: [email protected]
A Least Squares Method for Identification of Feedback Cascade Systems*
17th IFAC Symposium on System Identification 17th Symposium on Identification Beijing International Convention 17th IFAC IFAC Symposium on System SystemCenter Identification 17th IFAC Symposium on SystemCenter Identification Beijing Available online at www.sciencedirect.com October 19-21, 2015. Convention Beijing, China Beijing International International Convention Center Beijing International Convention Center October 19-21, 2015. Beijing, China October 19-21, 2015. Beijing, China October 19-21, 2015. Beijing, China
ScienceDirect
A A A A
IFAC-PapersOnLine 48-28 (2015) 098–103
Least Squares Method for Identification Least Squares Method for Identification Least Squares Method for Identification ⋆⋆ Least Squares Method for Identification of Feedback Cascade Systems ⋆⋆ of Feedback Cascade Systems of Cascade Systems of Feedback Feedback Cascade Systems ∗ ∗ ∗
Miguel Galrinho ∗∗ Cristian R. Rojas ∗∗ H˚ akan Hjalmarsson ∗∗ ∗ Cristian R. Rojas ∗ H˚ ∗ Miguel Galrinho akan Hjalmarsson Miguel Galrinho Cristian R. Rojas H˚ akan Hjalmarsson ∗ ∗ ∗ Miguel Galrinho Cristian R. Rojas H˚ akan Hjalmarsson ∗ Department of Automatic Control and ACCESS Linnaeus Center, ∗ ∗ Department of Automatic Control and ACCESS Linnaeus Center, Control ACCESS ∗ Department of Automatic School of Electrical Engineering, Department of Automatic Control and and ACCESS Linnaeus Linnaeus Center, Center, School of Electrical Engineering, School of Electrical Engineering, KTH – Royal Institute of Technology, SE-100 44 Stockholm, Sweden. School of Electrical Engineering, KTH – Royal Institute of Technology, SE-100 44 Stockholm, Sweden. KTH – Royal Institute of Technology, SE-100 44 Stockholm, E-mail: {galrinho, crro, hjalmars}@kth.se KTH – Royal Institute of Technology, SE-100 44 Stockholm, Sweden. Sweden. E-mail: {galrinho, crro, hjalmars}@kth.se E-mail: {galrinho, crro, hjalmars}@kth.se E-mail: {galrinho, crro, hjalmars}@kth.se Abstract: The problem of identification of systems in dynamic networks is considered. Although Abstract: The problem of identification of systems dynamic is considered. Although Abstract: The problem of of systems in dynamic networks is Although the prediction error method (PEM) can be applied toin the overallnetworks system, the non-standard model Abstract: Theerror problem of identification identification ofapplied systemsto inthe dynamic networks is considered. considered. Although the prediction method (PEM) can be overall system, the non-standard model the prediction error method (PEM) can be applied to the overall system, the non-standard model structure requires solving a non-convex optimization problem. Alternative methods have been the prediction errorsolving methoda (PEM) can be applied to the overall system, the non-standard model structure requires non-convex optimization problem. Alternative methods have been structure requires solving a non-convex optimization problem. Alternative methods have been proposed, such as instrumental variables and indirect PEM. In this paper, we first consider structure requires solving a non-convex optimization problem. Alternative methods have been proposed, such as as instrumental variables and indirect indirect PEM. In this paper, paper, we first consider proposed, such instrumental variables and In consider acyclic cascade systems, and argue that these methodsPEM. have different rangeswe of first applicability. proposed, such as instrumental variables and indirect PEM. In this this paper, we first consider acyclic cascade systems, and argue that these methods have different ranges of applicability. acyclic cascade systems, argue that have of applicability. Then, for a network withand feedback connection, we propose andifferent approachranges to deal with the fact acyclic cascade systems, and argue connection, that these these methods methods have different ranges of with applicability. Then, for a network with feedback we propose an approach to deal the fact Then, for a network with feedback connection, we propose an approach to deal with the fact that indirect PEM yields a non-convex problem in that case. A numerical simulation may Then, for a network with feedback connection, we in propose an approach to deal with the may fact that indirect PEM yields a non-convex problem that case. A numerical simulation that indirect PEM yields a non-convex problem in that case. A numerical simulation may indicate that this approach is competitive with other existing methods. that indirect PEM yields a non-convex problem in that case. A numerical simulation may indicate that this approach is competitive with other existing methods. indicate that this is with existing methods. indicate that (International this approach approach is competitive competitive with other other existing © 2015, IFAC Federation of Automatic Control) Hostingmethods. by Elsevier Ltd. All rights reserved. Keywords: System identification, Least-squares algorithm, Feedforward networks, Feedback Keywords: System identification, Keywords: System identification identification, Least-squares Least-squares algorithm, algorithm, Feedforward Feedforward networks, networks, Feedback Feedback loops, Closed-loop Keywords: System identification, Least-squares algorithm, Feedforward networks, Feedback loops, Closed-loop identification loops, Closed-loop identification loops, Closed-loop identification 1. INTRODUCTION operator. Mathematically, the system in Fig. 1 can be 1. INTRODUCTION operator. Mathematically, the system in Fig. 1 can be 1. INTRODUCTION operator. Mathematically, described by 1. INTRODUCTION operator. Mathematically, the the system system in in Fig. Fig. 1 1 can can be be described by by described Due to the rising complexity of systems encountered described by y1 (t) = G1 (q)u(t) + e1 (t) (1a) Due to complexity of encountered Dueengineering to the the rising rising complexity of systems systems encountered yy11 (t) = G (q)u(t) + e1 (t) (1a) in problems, identification of systems that (t) = =G G211 (q)G (q)u(t) + (1a) Due to the rising complexity of systems encountered + e2 (t). (1b) y21 (t) in engineering problems, identification of systems + ee11 (t) (t) (1a) 1 (q)u(t) 1 (q)u(t) in engineering problems, identification of systems that are embedded in a dynamic network has becomethat an (t) = G (q)G (q)u(t) + e (t). (1b) y 2 2 1 2 in engineering problems, identification of systems that (t) = G (q)G (q)u(t) + e (t). (1b) y 2 2 1 2 are embedded in a dynamic network has become an (t) = G (q)G (q)u(t) + e (t). (1b) y are embedded in a dynamic network has become an 2 2 1 2 increasingly relevant problem. Thus, several contributions are embedded in a problem. dynamicThus, network hascontributions become an First, notice that the transfer function G (q) can be esincreasingly increasingly relevant problem. in Thus, several contributions have recentlyrelevant been provided thisseveral area, e.g., Dankers First, notice that the transfer function G11 (q) can be esincreasingly relevant problem. Thus, several contributions (q) be esFirst, transfer function G have recently been provided in this area, e.g., Dankers withthat (1a)the from the signals u(t) y1 (t), have recently been provided in this this Van area,den e.g.,Hof Dankers (q) can can be usesFirst, notice notice that the transfer function G11and et al. recently (2014), been Everitt et al. (2014), et al. timated timated with (1a) from the signals u(t) and yyLikewise, (t), ushave provided in area, e.g., Dankers 1 timated with (1a) from the signals u(t) and using standard system identification techniques. et al. (2014), Everitt et al. (2014), Van den Hof et al. 1 (t), et al. (2014), Everitt et al. (2014), Van den Hof et al. timated with (1a) from the signals u(t) and y (t), us(2013), Dankers et al. (2013), Everitt et al. (2013), Van den 1 ing standard system identification techniques. Likewise, et al. (2014), Everitt et al. (2014), Van den Hof et al. ing standard system identification techniques. Likewise, (2013), et (2013), et Van den product =: G21 (q) can be estimated in a 2 (q)G1 (q) (2013), Dankers et al. al.H¨ (2013), Everitt et al. al. (2013), (2013), Van and den the ing standardG system identification techniques. Likewise, Hof et Dankers al. (2012), agg etEveritt al. (2011), Wahlberg the product G (q)G (q) =: G (q) can be estimated in a (2013), Dankers et al. (2013), Everitt et al. (2013), Van den 2 1 21 the product G (q)G (q) =: G (q) can be estimated in Hof et al. (2012), H¨ a gg et al. (2011), Wahlberg and similar fashion from (1b), using u(t) and y (t) as data. 2 1 21 2 Hof et et al. al.(2008), (2012), H¨ agg gg et etet al. al. (2011), Wahlberg Wahlberg and and similar the product G2 (q)G =: using G21 (q)u(t) canand be estimated in a a Sandberg Wahlberg al. (2008). 1 (q) fashion from (1b), y (t) as data. Hof (2012), H¨ a (2011), 2 similar fashion from (1b), using u(t) and y (t) as data. Sandberg (2008), Wahlberg et al. (2008). However, the input to(1b), the using transfer function indi2G 2 (q), Sandberg (2008), Wahlberg et al. (2008). similar fashion from u(t) and y (t) as data. 2 However, the input to the transfer function G (q), indiSandberg Wahlberg al. (2008). 2 the the transfer function G indiA common(2008), particular case ofetsuch networks is the identi- However, cated in Fig. as u2to not known. Therefore, 2 (q) However, the 11input input to(t), theis transfer function G22 (q), (q),G indiA common particular case of such networks is the identicated in Fig. as u (t), is not known. Therefore, G 2 2 (q) A common particular case of such networks is the identicated in Fig. 1 as u (t), is not known. Therefore, G fication of acyclic cascade structures, e.g., the system in cannot be estimated directly using a similar approach. 2 A common particular case structures, of such networks is the identicated inbeFig. 1 as u2 (t), is notusing known. Therefore, G22 (q) (q) fication of acyclic cascade e.g., the system in cannot estimated directly a similar approach. fication of acyclic cascade structures, e.g., the system in cannot be estimated using similar approach. Fig. 1. Itofcontains one external input, e.g., u(t), the andsystem two outpossible strategy todirectly obtain G from the previously 2 (q) a fication acyclic cascade structures, in A cannot be estimated directly using a similar approach. Fig. 1. It contains one external input, u(t), and two outA possibleestimates strategy to to obtain G22 (q) (q) from from the previously Fig. 1. 1.y1It It(t)contains contains onewith external input, u(t), u(t), and etwo two outputs, and y2 (t), measurement noises (t) outand A possible strategy G previously obtained of obtain G1 (q) and (q) the is to use the 21 Fig. one external input, and A possibleestimates strategy to obtain G2 (q)G from the previously puts, yyrespectively, (t) and yy22 (t), with measurement noises ee111 (t) and obtained of G (q) and G (q) is to use the 1 1 21 −1 puts, and (t), with measurement noises (t) and obtained estimates of G (q) and G (q) is to use the eputs, which, for the purpose of this paper, 1 (t) 2 (t), y 1−1 21 (q) that (q) = G (q)G (q). However, does relation G (t) and y (t), with measurement noises e (t) and 2 21 obtained estimates of G (q) and G is to use not the 1 2 1 1 1−1 (q). However, 21 eare (t), respectively, which, for the purpose of this paper, 2 (q) = G (q)G that does not relation G ee2 (t), respectively, which, for the purpose of this paper, 2 21 Gaussian, white, and uncorrelated to each other, with 1 (q) = G (q)G (q). However, that does not relation G −1 allow imposing a particular structure on G (q). Further2 (q) = G21 (q)G1 (q). However, that (t), respectively, which, for the purpose of this paper, 2 2 does not relation G are Gaussian, white, and uncorrelated to each other, with allow imposing 2 21 1 aa particular on G (q).previously Furtherare white, and uncorrelated to other, 2 variances λ1 and λ2 . A general discussion on identification allow structure on G Furtherif G1 (q) and G21 (q) arestructure estimated 2 (q). are Gaussian, Gaussian, white, and uncorrelated to each each other, with with more, allow imposing imposing a particular particular structure onin Gthe Furthervariances λ and λ . A general discussion on identification 2 (q). 1 2 more, if G (q) and G (q) are estimated in the previously variances λ and λ . A general discussion on identification 1 21 and variance analysis of this type of cascade systems is 1 and λ2 . A general discussion on identification more, if G (q) and G (q) are estimated in the previously presented way, information that could be useful for the 1 21 variances λ 1 analysis 2 more, if G (q) and G (q) are estimated in the previously and variance of this type of cascade systems is 1 21 presented way, information that could be useful for the and variance analysis of this type of cascade systems is taken in Wahlberg et al. (2008). presented way, information that could be useful for the estimation is neglected. For example, using also y (t) and variance analysis of this type of cascade systems is 2 presented way, information that could bealso useful forwhen the taken in Wahlberg et al. (2008). estimation is neglected. For example, using yy22estimates (t) when taken in Wahlberg et al. (2008). estimation is neglected. For example, using also (t) when estimating G (q) can improve the variance of the taken in Wahlberg et al. (2008). estimation G is 11neglected. For example, usingof also y2estimates (t) when estimating (q) can improve the variance the estimating improve (see EverittG al.can (2013)). 1 (q) u2 (t) estimating Get can improve the the variance variance of of the the estimates estimates 1 (q) (see Everitt et al. (2013)). u (t) 2 u(t) (see Everitt et al. (2013)). u (t) G (q) G1 (q) 2 2 (see Everitt et al. (2013)). u (t) u(t) 2 (q) G (q) G Another possibility, which solves the problem of imposu(t) G11 (q) G22 (q) u(t) Another possibility, which solves the problem of imposG G 1 (q) 2 (q) Another possibility, which solves the problem of imposing structure, is to estimate G (q) using y (t) and an 2 2 Another possibility, which solves the using problem of and imposing structure, is to estimate G (q) y (t) an e1 (t) e2 (t) 2 2 ing structure, is to estimate G (q) using y (t) and an estimate of u (t) as data. However, the presence of in2 2 2 (t) (t) ing structure, is to estimate G2 (q) using y2 (t) and an eee11 (t) eee22 (t) estimate of u (t) as data. However, the presence of in2 estimate of u (t) as data. However, the presence of in(t) (t) put noise makes an errors-in-variables (EIV) prob2 (t) this 1 2 estimate of u as data. However, the presence of in2 y1 (t) y2 (t) put noise makes this an errors-in-variables (EIV) probput noise makes an problem (S¨oderstr¨ om this (2007)). When applied to (EIV) this type of yy11 (t) yy22 (t) put noise makes this an errors-in-variables errors-in-variables (EIV) problem (S¨ ooderstr¨ oom (2007)). When applied to this type of y1 (t) (t) y2 (t) (t) lem (S¨ derstr¨ m (2007)). When applied to this type problem, standard system identification methods typically Fig. 1. Cascade system with two transfer functions. lem (S¨oderstr¨ om (2007)). When applied to thistypically type of of problem, standard system identification methods Fig. 1. Cascade system with two transfer functions. problem, standard system methods typically yield parameter that are not consistent. InFig. 1. 1. Cascade Cascade system system with with two two transfer transfer functions. functions. problem, standardestimates system identification identification methods typically Fig. yield parameter estimates that are not consistent. Inyield parameter estimates that consistent. Invariable (IV) methods (S¨onot derstr¨ om and StoThe goal of system identification is to estimate the transfer strumental yield parameter estimates that are are consistent. Instrumental variable (IV) methods (S¨ oonot derstr¨ oom and StoThe goal of system identification is to estimate the transfer strumental variable (IV) methods (S¨ derstr¨ m and Stoica (2002)) can be used to solve this problem, since The goal of system identification is to estimate the transfer (q) and G (q), where q is the forward-shift functions G strumental variable (IV) methods (S¨ o derstr¨ o m and Sto1 2 The goal of system identification is toq estimate the transfer ica (2002)) be used to solve this problem, since and G (q), where the forward-shift functions ica (2002)) can be solve problem, since 2 some choices can of instruments consistent estimates functions G G11 (q) (q) and G where qq is is the forward-shift 2 (q), ica (2002)) can be used used to toprovide solve this this problem, since (q) and G (q), where is the forward-shift functions G some choices of instruments provide consistent estimates 1 2 some choices of instruments provide consistent estimates ⋆ This work was supported by the European Research Council under (see, e.g., S¨ o derstr¨ o m and Mahata (2002) and Thil et al. some choices of instruments provide(2002) consistent estimates ⋆ (see, e.g., S¨ o derstr¨ o m and Mahata and Thil et al. work supported by the Research under ⋆ (see, e.g., ooderstr¨ oom and Mahata (2002) Thil et al. (2008)). AS¨ generalized IV approach for EIVand identification This work was was supported by contract the European European Research Council under theThis advanced grant LEARN, 267381 and byCouncil the Swedish ⋆ (see, e.g., S¨ derstr¨ m and Mahata (2002) and Thil et al. This work was supported by contract the European Research Council under (2008)). A generalized IV approach for EIV identification the advanced grant LEARN, 267381 and by the Swedish (2008)). A generalized IV approach for EIV identification the advanced grant LEARN, contract 267381 and by the Swedish in dynamic networks has been proposed in Dankers et al. Research Council under contract 621-2009-4017. (2008)). A generalized IV approach for EIV identification the advanced grant LEARN, contract 267381 and by the Swedish in dynamic networks has been proposed in Dankers et al. Research in Research Council Council under under contract contract 621-2009-4017. 621-2009-4017. in dynamic dynamic networks networks has has been been proposed proposed in in Dankers Dankers et et al. al. Research Council under contract 621-2009-4017.
Copyright © 2015, IFAC 2015 98 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © IFAC (International Federation of Automatic Control) Copyright IFAC 2015 98 Copyright ©under IFAC responsibility 2015 98 Control. Peer review© of International Federation of Automatic Copyright © IFAC 2015 98 10.1016/j.ifacol.2015.12.107
2015 IFAC SYSID October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
99
least squares (LS) problem. Defining the regression vector as ϕ⊤ (t) := [u(t − 1) u(t − 2) . . . u(t − n)] , it is possible to write y(t) = ϕ⊤ (t)θ + y˜(t), where ⊤ θ = [θ1 θ2 . . . θn ] . Further, if we define
(2014). However, using the system in Fig. 1 as an example, notice that y2 (t) is still not used when estimating G1 (q). If the prediction error method (PEM) is applied to a model structure parametrized according to (1) and to the individual structures of G1 (q) and G2 (q), the obtained estimates are asymptotically efficient (see, e.g., Ljung (1999)). However, for such a model structure, PEM requires, in general, solving a non-convex optimization problem.
⊤
y = [y(1) y(2) . . . y(N )] ,
In Wahlberg et al. (2008), indirect PEM (S¨oderstr¨om et al. (1991)) is suggested as a suitable method for identification of cascade systems. In this method, PEM is first applied to a higher-order model. In a second step, this model is reduced to the model of interest in an optimal way, in the sense that the obtained estimate has the same asymptotic properties as if PEM had been applied to the smaller model directly. If the model in the first step is easier to estimate than the model of interest, the original problem is simplified.
⊤
Φ⊤ = [ϕ(1) ϕ(2) . . . ϕ(N )] , and y˜ analogously to y, we can write y = Φ⊤ θ + y˜. An estimate of θ, which corresponds to the minimizer of (3), can be obtained by LS, computing −1 θˆ = ΦΦ⊤ Φy. (4)
We consider now the case when the true input is not known, but it can be measured, and is corrupted by measurement noise. In this case, the data is generated according to yo (t) = G(q)uo (t) u(t) = uo (t) + u ˜(t) , y(t) = yo (t) + y˜(t)
In this contribution, we restrict ourselves to the case that each transfer function is a single-input single-output (SISO) finite impulse response (FIR) model. First, in Section 2, we revisit the application of IV methods to EIV problems. In Section 3, we review indirect PEM, and extend the discussion in Wahlberg et al. (2008) regarding the application to cascade structures. In Section 4, we point out that this method can be applied even if not all the transfer function outputs in a network are measured. Then, we consider a feedback cascade structure in Section 5, for which indirect PEM does not avoid non-convexity, and propose an intermediate step using the method in Galrinho et al. (2014). A numerical simulation is presented in Section 6, followed by a discussion in Section 7.
˜(t) the input measurewhere uo (t) is the true input, u ment noise, and u(t) the measured input. This setting corresponds to an errors-in-variables (EIV) problem (S¨oderstr¨om (2007)). In this scenario, we have that y(t) = ϕ⊤ (t)θ + v(t, θ) , v(t, θ) = y˜(t) − ϕ˜⊤ (t)θ where, if ϕo is defined analogously to ϕ, but containing true input values uo , then ϕ(t) ˜ = ϕ(t) − ϕo (t). Because v(t, θ) is not white, if the parameter vector θ is estimated according to (4), the obtained estimate is not consistent.
2. ERRORS-IN-VARIABLES METHODS Consider the SISO system G(q), and assume that data is generated according to yo (t) = G(q)u(t) , (2) y(t) = yo (t) + y˜(t)
Instrumental variable (IV) methods are appropriate to deal with EIV problems. The basic idea of IV methods is to choose a vector of instruments z(t) that is uncorrelated with the error v(t, θ), while being highly correlated with ϕ(t). Then, for such an instrument vector, computing −1 θˆ = ZΦ⊤ Zy,
where yo (t) is the true system output, y(t) is the measured output, corrupted by noise y˜(t), and the input u(t) is assumed to be known. We introduce the assumption that G(q) is FIR, and parametrize it accordingly as G(q, θ) = θ1 q −1 + θ2 q −1 + · · · + θn q −n , and that y˜(t) is Gaussian white noise.
where
Z = [z(1) z(2) . . . z(N )] , yields a consistent estimate of θ under certain excitation conditions.
The prediction error method (PEM) serves as benchmark in the field, since it is well known to provide asymptotically efficient estimates if the model orders are correct (Ljung (1999)). The essential idea of PEM is to minimize a cost function of the prediction errors. In this setting, PEM consists on minimizing the cost function N 1 2 V (θ) = (y(t) − G(q, θ)u(t)) , (3) N t=1
There is no unique way to define z(t). One approach, proposed in S¨oderstr¨om and Mahata (2002), is to choose z ⊤ (t) = [u(t − 1 − du ) . . . u(t − du − nzu )] , where du ≥ n. Another possibility, proposed in Thil et al. (2008), is to also include past outputs in the instrument vector, according to z ⊤ (t) = −y(t − 1 − dy ) . . . −y(t − dy − nzy ) (5) u(t − 1 − du ) . . . u(t − du − nzu )] , where dy is at least the order of the filter (it must be a moving average (MA) filter) applied to the noise. For the considered FIR case, dy ≥ 0.
if a quadratic cost is used, and where N is the number of samples available. Then, the minimizer of (3) is an asymptotically efficient estimate of θ, if the model orders are correct. In general, PEM requires solving a non-convex optimization problem. However, for this particular model structure, the minimizer of (3) can be obtained by solving a 99
2015 IFAC SYSID 100 October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
Returning to the cascade example of Fig. 1, it is possible to observe that such methods are appropriate for the estimation of G2 (q). Using the notation introduced in this Section, uo (t) plays the role of u2 (t), the unknown input to the second transfer function, and y1 (t) corresponds to u(t), the measurable input. Furthermore, the fact that the transfer functions are part of a network creates new possibilities for choices of instruments, as discussed in Dankers et al. (2014). In this case, a natural instrument candidate when estimating G2 (q) is the external input u(t). However, note that although G1 (q) can be obtained from data {u(t), y1 (t)} using any standard system identification method, discarding the output y2 (t) in the estimation of G1 (q) increases the variance of the estimated parameters (Everitt et al. (2013)). The contribution Gunes et al. (2014) addresses a similar issue for the two-stage method. Moreover, although G2 (q) can be estimated consistently from {y1 (t), y2 (t)} using the IV approach, it was shown in Hjalmarsson et al. (2011) that the EIV setting, even for an asymptotically efficient estimator, does not reach the asymptotic properties of PEM. 3. INDIRECT PEM In the cascade case, due to the product G2 (q)G1 (q), structures that are linear in the model parameters, such as FIR or ARX, cannot be directly applied, and PEM requires solving a non-convex optimization problem. Indirect PEM (S¨ oderstr¨ om et al. (1991)) offers, in some cases, a workaround for the non-convexity of PEM. The general idea is as follows. Consider two nested model structures M1 and M2 , such that M1 ⊂ M2 . Moreover, suppose that M1 is parametrized by the parameter vector θ, while M2 is parametrized by α. In a first step, α is estimated with PEM, by minimizing N V (α) = ǫ2 (t, α), (6) t=1
where ǫ(t, α) are the prediction errors associated with the model structure M2 . In some situations, M2 can be chosen such that minimizing (6) becomes an easy problem, even if applying PEM to the model structure M1 is difficult. For examples, see S¨ oderstr¨ om et al. (1991). Then, in a second step, an estimate of θ is obtained by minimizing 1 1 V (θ) = [ˆ α − α(θ)]⊤ Pα−1 [ˆ α − α(θ)] , (7) N2 where α ˆ is the estimate of α obtained by minimizing (6), and Pα is a consistent estimate of the covariance of α ˆ, obtained from minimizing (6). Note that (7) is still non-quadratic; however, it is rather easy to handle. Asymptotically, a single step of the GaussNewton algorithm, initialized with a consistent estimate of θ, provides an efficient estimate, i.e., it has the same asymptotic properties as if PEM had been applied to the model of interest directly. In S¨ oderstr¨ om et al. (1991), more details are given on the algorithm. Indirect PEM can be applied to the cascade system in Fig. 1 as follows. First, since the transfer functions considered are FIR, G1 (q) can be estimated from {u(t), y1 (t)},
100
and G21 (q) from {u(t), y2 (t)}, by solving LS problems as in (4). Since the parameters estimated in this step were, in the exposition above, designated by α, we shall write that estimates G1 (q, α) ˆ and G21 (q, α ˆ ) have been obtained in this step. In the second step, we solve G1 (q, θ) = G1 (q, α ˆ) (8a) ˆ (8b) G2 (q, θ)G1 (q, θ) = G21 (q, α) for θ in a weighted least squares (WLS) sense, where the weighting is a consistent estimate of the inverse covariance of α, obtained from the LS estimator by 1 Φ1 Φ ⊤ 0 1 −1 λ 1 Pα = , 1 ⊤ 0 λ2 Φ2 Φ2 and Φ1 and Φ2 are the regressor matrices for the estimation of G1 (q) and G21 (q), respectively. To be able to obtain an asymptotically efficient estimate in one step of the Gauss-Newton algorithm, consistent estimates of G1 (q, θ) and G2 (q, θ) are required. An estimate of G1 (q, θ) is readily available from (8a), while a consistent estimate of G2 (q, θ) can be obtained by applying LS to (8b), replacing G1 (q, θ) by G1 (q, α ˆ ). So, if the estimates obtained in this ˆ and G2 (q, θ), ˆ are used as initializations to way, G1 (q, θ) the Gauss-Newton algorithm, the next step yields asymp¯ totically efficient estimates, θ. Thus, the asymptotic efficiency achieved by applying indirect PEM is an advantage in comparison to the EIV approach. On the other hand, for other structures, such as output-error, IV techniques can still be applied to obtain consistent estimates of G1 (q) and G2 (q), while indirect PEM suffers from non-convexity. There are, however, certain cascade structures for which the standard IV approach is not applicable, but to which indirect PEM still applies. An example is discussed in the next section. 4. ANOTHER CASCADE STRUCTURE Consider the cascade structure in Fig. 2, which can be written as y1 (t) 0 G2 (q) u1 (t) e1 (t) , = + e2 (t) y2 (t) G3 (q)G1 (q) G3 (q)G2 (q) u2 (t) (9) where we maintain the FIR model assumption. u1 (t)
G1 (q)
u2 (t)
G2 (q) e1 (t)
G3 (q) e2 (t) y2 (t)
y1 (t) Fig. 2. Cascade system with three transfer functions. An important difference when comparing to the structure in Fig. 1, besides the presence of one more transfer function block and one more input, is that not all block inputs and outputs are measured. In particular, the output of G1 (q), as well as the input of G3 (q), are unknown. Therefore, the IV approach delineated in Section 2 is not directly applicable for estimating G1 (q) and G3 (q).
2015 IFAC SYSID October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
On the other hand, indirect PEM is still applicable. By using the parametrization G2 (q) = G2 (q, α) G3 (q)G1 (q) = G31 (q, α) , G3 (q)G2 (q) = G32 (q, α) we can use (9) to estimate G2 (q, α), G31 (q, α), and G32 (q, α). Then, the second step of indirect PEM concerns solving G2 (q, θ) = G2 (q, α) ˆ (10a) ˆ (10b) G3 (q, θ)G1 (q, θ) = G31 (q, α) ˆ (10c) G3 (q, θ)G2 (q, θ) = G32 (q, α) in the discussed WLS sense. Since (10a) can be used to ˆ (10c) can obtain a consistent estimate of G2 (q, θ), G2 (q, θ), then be used to obtain a consistent estimate of G3 (q, θ), ˆ With a consistent by replacing G2 (q, θ) with G2 (q, θ). ˆ we can use (10b) to obtain a consistent estimate G3 (q, θ), estimate of G1 (q, θ). With a consistent estimate θˆ of the complete vector θ, one step of the Gauss-Newton algorithm provides an asymptotically efficient estimate of θ. 5. CASCADE WITH FEEDBACK In this section, we consider a cascade system with two blocks connected in a feedback loop, represented in Fig. 3. Mathematically, we can write it as 1 y1 (t) G1 (q) G1 (q)G2 (q) = y2 (t) G2 (q) 1 − G1 (q)G2 (q) G1 (q)G2 (q) r1 (t) e (t) + 1 . (11) r2 (t) e2 (t) e1 (t) r1 (t)
y1 (t)
G1 (q)
y2 (t)
G2 (q)
r2 (t)
101
To do so, we use the following procedure. In a first step, we estimate high-order FIR models of each closed-loop SISO transfer function in (11). In a second step, we reduce these estimates to estimates of G1 (q), G2 (q), and the product G1 (q)G2 (q), with the desired orders, by LS. Then, we reestimate by WLS, where we use the estimate obtained in step 2 to construct the weighting matrix. Finally, the indirect PEM algorithm can be applied. In detail, the first step consists on approximating (11) by CL y1 (t) (q, g12 ) r1 (t) G11 (q, g11 ) GCL e1 (t) 12 = + , CL y2 (t) r2 (t) e2 (t) GCL 21 (q, g21 ) G22 (q, g22 ) (12) where m (k) CL gij q −k , Gij (q, gij ) = k=1
for i, j = {1, 2}. If the order m is chosen large enough, we have that 1 G1 (q) G1 (q)G2 (q) ≈ G2 (q) 1 − G1 (q)G2 (q) G1 (q)G2 (q) CL (13) G11 (q, g11 ) GCL (q, g12 ) 12 ≈ . CL GCL 21 (q, g21 ) G22 (q, g22 )
Since (12) is a multi-input multi-output (MIMO) FIR model, GCL ij (q, gij ) can be estimated by LS. This is done by solving, analogously to what was done from (2) to (4), −1 gˆ Φyi , gˆi := i1 = ΦΦ⊤ gˆi2 ⊤ Φ ⊤ = Φ⊤ r 1 Φr 2 , where Φr1 and Φr2 are the regression matrices for r1 (t) and r2 (t), respectively. The estimated parameters are, then, distributed as gˆi ∼ N (gi , Pgi ) , where the covariance matrix is given by −1 Pgi = λi ΦΦ⊤ . In a second step, we use the approximation (13) to write
e2 (t) Fig. 3. Block diagram of a feedback cascade system. Notice that this case is different from the typical feedback setting, with the measured output included in the loop. Instead, it can be seen as a natural extension of the previous cascade systems, where the systems are connected physically. A well-known historical example of a system physically connected in feedback is a centrifugal governor. Application of an IV method, in the EIV setting, to estimate G1 (q) and G2 (q) is straightforward, as both outputs are measured, and these are also the inputs to the other transfer function, summed with a known external reference. Note, however, that the same limitation of Section 4 would be encountered if, for example, output y2 (t) could not be measured. Concerning indirect PEM, it is difficult to find a model structure of finite order that is easy to estimate, and that contains the true model. The purpose of this section is to propose an intermediate step, based on the method presented in Galrinho et al. (2014), which reduces the problem to a setting that allows the application of the indirect PEM algorithm as in the previous sections. 101
Zii (q, α, gii ) :=GCL ii (q, gii ) (1 − G12 (q, α12 )) − Gi (q, αi ) = 0 Zij (q, α, gij ) :=GCL ij (q, gij ) (1 − G12 (q, α12 )) − G12 (q, α12 ) = 0
(14)
(in (14) only, i �= j), where we parametrize G1 (q) and G2 (q) as (1)
(ni ) −ni
Gi (q, αi ) = αi q −1 + · · · + αi and, to avoid the product G1 (q)G2 (q), (2)
q
(n +n −1)
.
G12 (q, α12 ) = α12 q −2 + · · · + α121 2 q −(n1 +n2 −1) . Then, if gij is replaced by its estimate, gˆij , (14) can be solved by LS. Write, first, (14) in vector form, as Q(ˆ g )α = gˆ, (15) where ⊤ ⊤ ⊤, α = α⊤ 1 α2 α12 ⊤ ⊤ ⊤ ⊤ ⊤ g = g11 g12 g21 g22 , , and obtain estimates α ˆ1, α ˆ 2 , and α ˆ 12 by computing the LS estimate of α from (15), i.e., −1 g )⊤ Q(ˆ g) Q(ˆ g )⊤ gˆ. α ˆ = Q(ˆ
2015 IFAC SYSID 102 October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
In a third step, we re-solve (14), now using WLS, where the weighting is obtained as follows. Let gˆij = gij + g˜ij . Then, we can write, in transfer function form, CL GCL ˜ij ). (16) ˆij ) = GCL ij (q, gij ) + Gij (q, g ij (q, g
100 98 96 94
˜ij ). Zij (q, α, gˆij ) = (1 − G12 (q, α12 )) GCL ij (q, g
(17)
Rewriting the equations defined in (17) in vector form, we have g˜11 z11 (α, gˆ11 ) g˜12 z12 (α, gˆ12 ) g = T (α12 ) =: T (α12 )˜ z(α, gˆ) := z21 (α, gˆ21 ) g˜21 z22 (α, gˆ22 ) g˜22 Then, the optimal weighting for the WLS is the inverse of the covariance of z(α, gˆ), i.e., �−1 � W (α12 ) = T (α12 )Pg˜ T (α12 )⊤ , where Pg˜ is the covariance of g˜, given by � � Pg1 0 . Pg˜ = 0 Pg2
Since the true α12 is not available, we use, instead, the estimate obtained in step 2, �−1 � W (ˆ α12 ) = T (ˆ α12 )Pg˜ T (ˆ α12 )⊤ , and compute � �−1 α ¯ = Q(ˆ g )⊤ W (ˆ α12 )Q(ˆ g) Q(ˆ g )⊤ W (ˆ α12 )ˆ g, which is an improved estimate of α. Asymptotically, we have that these estimates have covariance �−1 � α12 )Q(ˆ g) . (18) cov (α) ¯ = Q(ˆ g )⊤ W (ˆ
With these estimates, we can construct transfer functions ¯ 1 ), G2 (q, α ¯ 2 ), and G12 (q, α ¯ 12 ), and solve G1 (q, α G1 (q, θ) = G1 (q, α ¯1 ) (19a) ¯2 ) (19b) G2 (q, θ) = G2 (q, α ¯ 12 ) (19c) G1 (q, θ)G2 (q, θ) = G12 (q, α with the algorithm previously mentioned for indirect PEM, using (18) as covariance estimate for α ¯.
We note that, like indirect PEM for the acyclic cascade systems, this method does not require all the transfer function outputs to be measured. For example, if y2 (t) were not measured, the method would be applicable without major changes. In particular, (19b) would be absent in (19). 6. NUMERICAL SIMULATION In this section, we perform a numerical simulation to evaluate the method previously described for identification of the system in Fig. 3. One hundred Monte Carlo runs are performed, for sample sizes N = {100, 300, 600, 1000, 3000, 6000}. The systems are given by � (1) (2) G1 (q) = θ1 q −1 + θ1 q −2 . (1) (2) G2 (q) = θ2 q −1 + θ2 q −2 The closed-loop, whose poles are given by 1 − G1 (q)G2 (q) = 0, (20) is, thus, a fourth order system. The poles are placed randomly inside the unit circle, under the structure imposed 102
FIT
Replacing (16) in (14) yields
92 90 88 PEM WLSiPEM IVr IVuy
86 84 102
103
N
Fig. 4. Average FIT, for 100 Monte Carlo simulations, as function of sample size, comparing PEM, two IV methods, and the proposed WLS method combined with indirect PEM (WLSiPEM). by (20). The same colored input as in S¨oderstr¨ om and Mahata (2002) and Thil et al. (2008) is used, i.e., 1 rw (t), ri (t) = −1 1 − 0.2q + 0.5q −2 i for i = 1, 2, where rwi (t) is Gaussian white noise with unitary variance, and uncorrelated for each i. The sensor noises e1 (t) and e2 (t) are also Gaussian and white, and uncorrelated to each other. The variances are changed at each run in order to achieve a signal-to-noise ratio of 5dB at the outputs. The following methods are compared: PEM, as implemented in MATLAB2014b, starting at the true parameters (note that MATLAB does not provide an estimated point to initialize the algorithm, due to the non-standard model structure); the proposed method, referred to as WLSiPEM; the IV method with the external references as instruments (IVr); the IV method with instrument vector (5) (IVuy). Concerning WLSiPEM, the order m of the high-order FIR models is optimized over a grid of ten values, linearly spaced between m = 20 and m = N/3, and the one minimizing the prediction errors is chosen. A similar procedure is used to choose the length of the instrument vector for IV methods, by considering a grid of ten values linearly spaced between 4 and 22, and selecting the one minimizing the prediction errors. All the data available is used in the estimation. The accuracy of the estimates is measured by the FIT, in percent, defined as � � ˆ 2 ||θ − θ|| FIT = 100 1 − ¯ 2 , ||θ − θ|| where θ is a vector containing the true parameters, θ¯ its mean, and θˆ the estimated parameter. The results are presented in Fig. 4, with the average FIT as function of sample size. In this simulation, the presented method performed better than both IV methods. On the other hand, other simulations showed the IV methods
2015 IFAC SYSID October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
to be more robust in the presence of extremely colored external excitations. It is also observed that using other variables in the network to construct the instrument vector improves the estimate, in comparison to using delayed measured input and output values. Finally, we verified that, in this case, the estimates obtained by any of the methods were appropriate as initial estimates for PEM, since the algorithm converged to the same parameters as when it was started at the true values.
7. DISCUSSION In this paper, we discussed two approaches for identification of network systems. One is based on the application of an IV method to an EIV setting, while the other uses indirect PEM. We point out that the applicability of each method depends on the particular problem setting and the assumptions on the network. The IV methods require that the input and output of all transfer functions are measured, while that is not a requirement for indirect PEM. Moreover, while the accuracy of IV methods depends on the choice of instruments and the correlation to the regression vector, indirect PEM guarantees asymptotic efficiency. On the other hand, standard IV methods admit rational transfer functions, while indirect PEM suffers from non-convexity in that case. If there is feedback in the network, indirect PEM does not avoid a non-convex problem. To deal with that limitation, we propose a method that uses the approach in Galrinho et al. (2014), and transforms the problem into one where the indirect PEM algorithm can be applied. Analogously to the acyclic case, measuring all outputs is not a necessary condition to obtain all the transfer functions. We perform a simulation study with two FIR transfer functions randomly generated and connected in feedback, where both outputs are measured, comparing PEM, the proposed method, and two IV methods with different instrument choices. The proposed approach performed better than the IV methods, but it did not prove as robust in the presence of extremely colored external excitations. This reinforces our discussion on the most appropriate approach depending on the settings and the network assumptions. To provide initial values for PEM, any of the techniques were successful. The noise sources were considered white, but most methods discussed allow the noise to be colored. For the IV method in Thil et al. (2008), the output noise must be given by a MA filter of white noise, although the input noise is assumed white. In the network context, the more general case, with an autoregressive moving average (ARMA) filter, is addressed in Dankers et al. (2014). Concerning indirect PEM, we still have a linear regression problem if the noise is given by an autoregressive filter driven by white noise. Finally, the proposed method is still applicable in the ARMA filter case, for which a high-order ARX model is estimated in the first step of the algorithm. Future work includes establishing conditions on the network structure for the methodology to be applicable, and extending the approach to the case where the blocks in the network are rational transfer functions. 103
103
REFERENCES Dankers, A., Van den Hof, P., and Heuberger, P. (2013). Predictor input selection for direct identification in dynamic networks. In 52nd Conference on Decision and Control, 4541–4546. Dankers, A., den Hof, P.M.V., Bombois, X., and Heuberger, P.S. (2014). Errors-in-variables identification in dynamic networks by an instrumental variable approach. In 19th IFAC World Congress, 2335–2340. Everitt, N., Rojas, C.R., and Hjalmarsson, H. (2013). A geometric approach to variance analysis of cascaded systems. In 52nd IEEE Conference on Decision and Control, 6496–6501. Everitt, N., Rojas, C.R., and Hjalmarsson, H. (2014). Variance results for parallel cascade serial systems. In 19th IFAC World Congress, 2317–2322. Galrinho, M., Rojas, C.R., and Hjalmarsson, H. (2014). A weighted least-squares method for parameter estimation in structured models. In 53rd IEEE Conference on Decision and Control, 3322–3327. Gunes, B., Dankers, A., and den Hof, P.M.V. (2014). Variance reduction for identification in dynamic networks. In 19th IFAC World Congress, 2842–2847. H¨agg, P., Wahlberg, B., and Sandberg, H. (2011). On identification of parallel cascade serial systems. In 18th IFAC World Congress, 9978–9983. Hjalmarsson, H., M˚ artensson, J., Rojas, C.R., and S¨oderstr¨om, T. (2011). On the accuracy in errorsin-variables identification compared to prediction-error identification. Automatica, 47(12), 2704–2712. Ljung, L. (1999). System Identification. Theory for the User, 2nd ed. Prentice-Hall. S¨oderstr¨om, T. (2007). Errors-in-variables methods in system identification. Automatica, 43(6), 939–958. S¨oderstr¨om, T. and Mahata, K. (2002). On instrumental variable and total least squares approaches for identification of noisy systems. International Journal of Control, 75(6), 381–389. S¨oderstr¨om, T., Stoica, P., and Friedlander, B. (1991). An Indirect Prediction Error Method for System Identification. Automatica, 27(1), 183–188. S¨oderstr¨om, T. and Stoica, P. (2002). Instrumental variable methods for system identification. Circuits, Systems and Signal Processing, 21(1), 1–9. Thil, S., Gilson, M., and Garnier, H. (2008). On instrumental variable-based methods for errors-in-variables model identification. In 17th IFAC World Congress, 426–431. Van den Hof, P.M., Dankers, A., Heuberger, P.S., and Bombois, X. (2013). Identification of dynamic models in complex networks with prediction error methods – Basic methods for consistent module estimates. Automatica, 49(10), 2994 – 3006. Van den Hof, P., Dankers, A., Heuberger, P., and Bombois, X. (2012). Identification in dynamic networks with known interconnection topology. In 51st Conference on Decision and Control, 895–900. Wahlberg, B. and Sandberg, H. (2008). Cascade structural model approximation of identified state space models. In 47th Conference on Decision and Control, 4198–4203. Wahlberg, B., Hjalmarsson, H., and M˚ artensson, J. (2008). On Identification of Cascade Systems. In Proc. 17th IFAC World Congress, 5036–5040.
ScienceDirect
A A A A
IFAC-PapersOnLine 48-28 (2015) 098–103
Least Squares Method for Identification Least Squares Method for Identification Least Squares Method for Identification ⋆⋆ Least Squares Method for Identification of Feedback Cascade Systems ⋆⋆ of Feedback Cascade Systems of Cascade Systems of Feedback Feedback Cascade Systems ∗ ∗ ∗
Miguel Galrinho ∗∗ Cristian R. Rojas ∗∗ H˚ akan Hjalmarsson ∗∗ ∗ Cristian R. Rojas ∗ H˚ ∗ Miguel Galrinho akan Hjalmarsson Miguel Galrinho Cristian R. Rojas H˚ akan Hjalmarsson ∗ ∗ ∗ Miguel Galrinho Cristian R. Rojas H˚ akan Hjalmarsson ∗ Department of Automatic Control and ACCESS Linnaeus Center, ∗ ∗ Department of Automatic Control and ACCESS Linnaeus Center, Control ACCESS ∗ Department of Automatic School of Electrical Engineering, Department of Automatic Control and and ACCESS Linnaeus Linnaeus Center, Center, School of Electrical Engineering, School of Electrical Engineering, KTH – Royal Institute of Technology, SE-100 44 Stockholm, Sweden. School of Electrical Engineering, KTH – Royal Institute of Technology, SE-100 44 Stockholm, Sweden. KTH – Royal Institute of Technology, SE-100 44 Stockholm, E-mail: {galrinho, crro, hjalmars}@kth.se KTH – Royal Institute of Technology, SE-100 44 Stockholm, Sweden. Sweden. E-mail: {galrinho, crro, hjalmars}@kth.se E-mail: {galrinho, crro, hjalmars}@kth.se E-mail: {galrinho, crro, hjalmars}@kth.se Abstract: The problem of identification of systems in dynamic networks is considered. Although Abstract: The problem of identification of systems dynamic is considered. Although Abstract: The problem of of systems in dynamic networks is Although the prediction error method (PEM) can be applied toin the overallnetworks system, the non-standard model Abstract: Theerror problem of identification identification ofapplied systemsto inthe dynamic networks is considered. considered. Although the prediction method (PEM) can be overall system, the non-standard model the prediction error method (PEM) can be applied to the overall system, the non-standard model structure requires solving a non-convex optimization problem. Alternative methods have been the prediction errorsolving methoda (PEM) can be applied to the overall system, the non-standard model structure requires non-convex optimization problem. Alternative methods have been structure requires solving a non-convex optimization problem. Alternative methods have been proposed, such as instrumental variables and indirect PEM. In this paper, we first consider structure requires solving a non-convex optimization problem. Alternative methods have been proposed, such as as instrumental variables and indirect indirect PEM. In this paper, paper, we first consider proposed, such instrumental variables and In consider acyclic cascade systems, and argue that these methodsPEM. have different rangeswe of first applicability. proposed, such as instrumental variables and indirect PEM. In this this paper, we first consider acyclic cascade systems, and argue that these methods have different ranges of applicability. acyclic cascade systems, argue that have of applicability. Then, for a network withand feedback connection, we propose andifferent approachranges to deal with the fact acyclic cascade systems, and argue connection, that these these methods methods have different ranges of with applicability. Then, for a network with feedback we propose an approach to deal the fact Then, for a network with feedback connection, we propose an approach to deal with the fact that indirect PEM yields a non-convex problem in that case. A numerical simulation may Then, for a network with feedback connection, we in propose an approach to deal with the may fact that indirect PEM yields a non-convex problem that case. A numerical simulation that indirect PEM yields a non-convex problem in that case. A numerical simulation may indicate that this approach is competitive with other existing methods. that indirect PEM yields a non-convex problem in that case. A numerical simulation may indicate that this approach is competitive with other existing methods. indicate that this is with existing methods. indicate that (International this approach approach is competitive competitive with other other existing © 2015, IFAC Federation of Automatic Control) Hostingmethods. by Elsevier Ltd. All rights reserved. Keywords: System identification, Least-squares algorithm, Feedforward networks, Feedback Keywords: System identification, Keywords: System identification identification, Least-squares Least-squares algorithm, algorithm, Feedforward Feedforward networks, networks, Feedback Feedback loops, Closed-loop Keywords: System identification, Least-squares algorithm, Feedforward networks, Feedback loops, Closed-loop identification loops, Closed-loop identification loops, Closed-loop identification 1. INTRODUCTION operator. Mathematically, the system in Fig. 1 can be 1. INTRODUCTION operator. Mathematically, the system in Fig. 1 can be 1. INTRODUCTION operator. Mathematically, described by 1. INTRODUCTION operator. Mathematically, the the system system in in Fig. Fig. 1 1 can can be be described by by described Due to the rising complexity of systems encountered described by y1 (t) = G1 (q)u(t) + e1 (t) (1a) Due to complexity of encountered Dueengineering to the the rising rising complexity of systems systems encountered yy11 (t) = G (q)u(t) + e1 (t) (1a) in problems, identification of systems that (t) = =G G211 (q)G (q)u(t) + (1a) Due to the rising complexity of systems encountered + e2 (t). (1b) y21 (t) in engineering problems, identification of systems + ee11 (t) (t) (1a) 1 (q)u(t) 1 (q)u(t) in engineering problems, identification of systems that are embedded in a dynamic network has becomethat an (t) = G (q)G (q)u(t) + e (t). (1b) y 2 2 1 2 in engineering problems, identification of systems that (t) = G (q)G (q)u(t) + e (t). (1b) y 2 2 1 2 are embedded in a dynamic network has become an (t) = G (q)G (q)u(t) + e (t). (1b) y are embedded in a dynamic network has become an 2 2 1 2 increasingly relevant problem. Thus, several contributions are embedded in a problem. dynamicThus, network hascontributions become an First, notice that the transfer function G (q) can be esincreasingly increasingly relevant problem. in Thus, several contributions have recentlyrelevant been provided thisseveral area, e.g., Dankers First, notice that the transfer function G11 (q) can be esincreasingly relevant problem. Thus, several contributions (q) be esFirst, transfer function G have recently been provided in this area, e.g., Dankers withthat (1a)the from the signals u(t) y1 (t), have recently been provided in this this Van area,den e.g.,Hof Dankers (q) can can be usesFirst, notice notice that the transfer function G11and et al. recently (2014), been Everitt et al. (2014), et al. timated timated with (1a) from the signals u(t) and yyLikewise, (t), ushave provided in area, e.g., Dankers 1 timated with (1a) from the signals u(t) and using standard system identification techniques. et al. (2014), Everitt et al. (2014), Van den Hof et al. 1 (t), et al. (2014), Everitt et al. (2014), Van den Hof et al. timated with (1a) from the signals u(t) and y (t), us(2013), Dankers et al. (2013), Everitt et al. (2013), Van den 1 ing standard system identification techniques. Likewise, et al. (2014), Everitt et al. (2014), Van den Hof et al. ing standard system identification techniques. Likewise, (2013), et (2013), et Van den product =: G21 (q) can be estimated in a 2 (q)G1 (q) (2013), Dankers et al. al.H¨ (2013), Everitt et al. al. (2013), (2013), Van and den the ing standardG system identification techniques. Likewise, Hof et Dankers al. (2012), agg etEveritt al. (2011), Wahlberg the product G (q)G (q) =: G (q) can be estimated in a (2013), Dankers et al. (2013), Everitt et al. (2013), Van den 2 1 21 the product G (q)G (q) =: G (q) can be estimated in Hof et al. (2012), H¨ a gg et al. (2011), Wahlberg and similar fashion from (1b), using u(t) and y (t) as data. 2 1 21 2 Hof et et al. al.(2008), (2012), H¨ agg gg et etet al. al. (2011), Wahlberg Wahlberg and and similar the product G2 (q)G =: using G21 (q)u(t) canand be estimated in a a Sandberg Wahlberg al. (2008). 1 (q) fashion from (1b), y (t) as data. Hof (2012), H¨ a (2011), 2 similar fashion from (1b), using u(t) and y (t) as data. Sandberg (2008), Wahlberg et al. (2008). However, the input to(1b), the using transfer function indi2G 2 (q), Sandberg (2008), Wahlberg et al. (2008). similar fashion from u(t) and y (t) as data. 2 However, the input to the transfer function G (q), indiSandberg Wahlberg al. (2008). 2 the the transfer function G indiA common(2008), particular case ofetsuch networks is the identi- However, cated in Fig. as u2to not known. Therefore, 2 (q) However, the 11input input to(t), theis transfer function G22 (q), (q),G indiA common particular case of such networks is the identicated in Fig. as u (t), is not known. Therefore, G 2 2 (q) A common particular case of such networks is the identicated in Fig. 1 as u (t), is not known. Therefore, G fication of acyclic cascade structures, e.g., the system in cannot be estimated directly using a similar approach. 2 A common particular case structures, of such networks is the identicated inbeFig. 1 as u2 (t), is notusing known. Therefore, G22 (q) (q) fication of acyclic cascade e.g., the system in cannot estimated directly a similar approach. fication of acyclic cascade structures, e.g., the system in cannot be estimated using similar approach. Fig. 1. Itofcontains one external input, e.g., u(t), the andsystem two outpossible strategy todirectly obtain G from the previously 2 (q) a fication acyclic cascade structures, in A cannot be estimated directly using a similar approach. Fig. 1. It contains one external input, u(t), and two outA possibleestimates strategy to to obtain G22 (q) (q) from from the previously Fig. 1. 1.y1It It(t)contains contains onewith external input, u(t), u(t), and etwo two outputs, and y2 (t), measurement noises (t) outand A possible strategy G previously obtained of obtain G1 (q) and (q) the is to use the 21 Fig. one external input, and A possibleestimates strategy to obtain G2 (q)G from the previously puts, yyrespectively, (t) and yy22 (t), with measurement noises ee111 (t) and obtained of G (q) and G (q) is to use the 1 1 21 −1 puts, and (t), with measurement noises (t) and obtained estimates of G (q) and G (q) is to use the eputs, which, for the purpose of this paper, 1 (t) 2 (t), y 1−1 21 (q) that (q) = G (q)G (q). However, does relation G (t) and y (t), with measurement noises e (t) and 2 21 obtained estimates of G (q) and G is to use not the 1 2 1 1 1−1 (q). However, 21 eare (t), respectively, which, for the purpose of this paper, 2 (q) = G (q)G that does not relation G ee2 (t), respectively, which, for the purpose of this paper, 2 21 Gaussian, white, and uncorrelated to each other, with 1 (q) = G (q)G (q). However, that does not relation G −1 allow imposing a particular structure on G (q). Further2 (q) = G21 (q)G1 (q). However, that (t), respectively, which, for the purpose of this paper, 2 2 does not relation G are Gaussian, white, and uncorrelated to each other, with allow imposing 2 21 1 aa particular on G (q).previously Furtherare white, and uncorrelated to other, 2 variances λ1 and λ2 . A general discussion on identification allow structure on G Furtherif G1 (q) and G21 (q) arestructure estimated 2 (q). are Gaussian, Gaussian, white, and uncorrelated to each each other, with with more, allow imposing imposing a particular particular structure onin Gthe Furthervariances λ and λ . A general discussion on identification 2 (q). 1 2 more, if G (q) and G (q) are estimated in the previously variances λ and λ . A general discussion on identification 1 21 and variance analysis of this type of cascade systems is 1 and λ2 . A general discussion on identification more, if G (q) and G (q) are estimated in the previously presented way, information that could be useful for the 1 21 variances λ 1 analysis 2 more, if G (q) and G (q) are estimated in the previously and variance of this type of cascade systems is 1 21 presented way, information that could be useful for the and variance analysis of this type of cascade systems is taken in Wahlberg et al. (2008). presented way, information that could be useful for the estimation is neglected. For example, using also y (t) and variance analysis of this type of cascade systems is 2 presented way, information that could bealso useful forwhen the taken in Wahlberg et al. (2008). estimation is neglected. For example, using yy22estimates (t) when taken in Wahlberg et al. (2008). estimation is neglected. For example, using also (t) when estimating G (q) can improve the variance of the taken in Wahlberg et al. (2008). estimation G is 11neglected. For example, usingof also y2estimates (t) when estimating (q) can improve the variance the estimating improve (see EverittG al.can (2013)). 1 (q) u2 (t) estimating Get can improve the the variance variance of of the the estimates estimates 1 (q) (see Everitt et al. (2013)). u (t) 2 u(t) (see Everitt et al. (2013)). u (t) G (q) G1 (q) 2 2 (see Everitt et al. (2013)). u (t) u(t) 2 (q) G (q) G Another possibility, which solves the problem of imposu(t) G11 (q) G22 (q) u(t) Another possibility, which solves the problem of imposG G 1 (q) 2 (q) Another possibility, which solves the problem of imposing structure, is to estimate G (q) using y (t) and an 2 2 Another possibility, which solves the using problem of and imposing structure, is to estimate G (q) y (t) an e1 (t) e2 (t) 2 2 ing structure, is to estimate G (q) using y (t) and an estimate of u (t) as data. However, the presence of in2 2 2 (t) (t) ing structure, is to estimate G2 (q) using y2 (t) and an eee11 (t) eee22 (t) estimate of u (t) as data. However, the presence of in2 estimate of u (t) as data. However, the presence of in(t) (t) put noise makes an errors-in-variables (EIV) prob2 (t) this 1 2 estimate of u as data. However, the presence of in2 y1 (t) y2 (t) put noise makes this an errors-in-variables (EIV) probput noise makes an problem (S¨oderstr¨ om this (2007)). When applied to (EIV) this type of yy11 (t) yy22 (t) put noise makes this an errors-in-variables errors-in-variables (EIV) problem (S¨ ooderstr¨ oom (2007)). When applied to this type of y1 (t) (t) y2 (t) (t) lem (S¨ derstr¨ m (2007)). When applied to this type problem, standard system identification methods typically Fig. 1. Cascade system with two transfer functions. lem (S¨oderstr¨ om (2007)). When applied to thistypically type of of problem, standard system identification methods Fig. 1. Cascade system with two transfer functions. problem, standard system methods typically yield parameter that are not consistent. InFig. 1. 1. Cascade Cascade system system with with two two transfer transfer functions. functions. problem, standardestimates system identification identification methods typically Fig. yield parameter estimates that are not consistent. Inyield parameter estimates that consistent. Invariable (IV) methods (S¨onot derstr¨ om and StoThe goal of system identification is to estimate the transfer strumental yield parameter estimates that are are consistent. Instrumental variable (IV) methods (S¨ oonot derstr¨ oom and StoThe goal of system identification is to estimate the transfer strumental variable (IV) methods (S¨ derstr¨ m and Stoica (2002)) can be used to solve this problem, since The goal of system identification is to estimate the transfer (q) and G (q), where q is the forward-shift functions G strumental variable (IV) methods (S¨ o derstr¨ o m and Sto1 2 The goal of system identification is toq estimate the transfer ica (2002)) be used to solve this problem, since and G (q), where the forward-shift functions ica (2002)) can be solve problem, since 2 some choices can of instruments consistent estimates functions G G11 (q) (q) and G where qq is is the forward-shift 2 (q), ica (2002)) can be used used to toprovide solve this this problem, since (q) and G (q), where is the forward-shift functions G some choices of instruments provide consistent estimates 1 2 some choices of instruments provide consistent estimates ⋆ This work was supported by the European Research Council under (see, e.g., S¨ o derstr¨ o m and Mahata (2002) and Thil et al. some choices of instruments provide(2002) consistent estimates ⋆ (see, e.g., S¨ o derstr¨ o m and Mahata and Thil et al. work supported by the Research under ⋆ (see, e.g., ooderstr¨ oom and Mahata (2002) Thil et al. (2008)). AS¨ generalized IV approach for EIVand identification This work was was supported by contract the European European Research Council under theThis advanced grant LEARN, 267381 and byCouncil the Swedish ⋆ (see, e.g., S¨ derstr¨ m and Mahata (2002) and Thil et al. This work was supported by contract the European Research Council under (2008)). A generalized IV approach for EIV identification the advanced grant LEARN, 267381 and by the Swedish (2008)). A generalized IV approach for EIV identification the advanced grant LEARN, contract 267381 and by the Swedish in dynamic networks has been proposed in Dankers et al. Research Council under contract 621-2009-4017. (2008)). A generalized IV approach for EIV identification the advanced grant LEARN, contract 267381 and by the Swedish in dynamic networks has been proposed in Dankers et al. Research in Research Council Council under under contract contract 621-2009-4017. 621-2009-4017. in dynamic dynamic networks networks has has been been proposed proposed in in Dankers Dankers et et al. al. Research Council under contract 621-2009-4017.
Copyright © 2015, IFAC 2015 98 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © IFAC (International Federation of Automatic Control) Copyright IFAC 2015 98 Copyright ©under IFAC responsibility 2015 98 Control. Peer review© of International Federation of Automatic Copyright © IFAC 2015 98 10.1016/j.ifacol.2015.12.107
2015 IFAC SYSID October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
99
least squares (LS) problem. Defining the regression vector as ϕ⊤ (t) := [u(t − 1) u(t − 2) . . . u(t − n)] , it is possible to write y(t) = ϕ⊤ (t)θ + y˜(t), where ⊤ θ = [θ1 θ2 . . . θn ] . Further, if we define
(2014). However, using the system in Fig. 1 as an example, notice that y2 (t) is still not used when estimating G1 (q). If the prediction error method (PEM) is applied to a model structure parametrized according to (1) and to the individual structures of G1 (q) and G2 (q), the obtained estimates are asymptotically efficient (see, e.g., Ljung (1999)). However, for such a model structure, PEM requires, in general, solving a non-convex optimization problem.
⊤
y = [y(1) y(2) . . . y(N )] ,
In Wahlberg et al. (2008), indirect PEM (S¨oderstr¨om et al. (1991)) is suggested as a suitable method for identification of cascade systems. In this method, PEM is first applied to a higher-order model. In a second step, this model is reduced to the model of interest in an optimal way, in the sense that the obtained estimate has the same asymptotic properties as if PEM had been applied to the smaller model directly. If the model in the first step is easier to estimate than the model of interest, the original problem is simplified.
⊤
Φ⊤ = [ϕ(1) ϕ(2) . . . ϕ(N )] , and y˜ analogously to y, we can write y = Φ⊤ θ + y˜. An estimate of θ, which corresponds to the minimizer of (3), can be obtained by LS, computing −1 θˆ = ΦΦ⊤ Φy. (4)
We consider now the case when the true input is not known, but it can be measured, and is corrupted by measurement noise. In this case, the data is generated according to yo (t) = G(q)uo (t) u(t) = uo (t) + u ˜(t) , y(t) = yo (t) + y˜(t)
In this contribution, we restrict ourselves to the case that each transfer function is a single-input single-output (SISO) finite impulse response (FIR) model. First, in Section 2, we revisit the application of IV methods to EIV problems. In Section 3, we review indirect PEM, and extend the discussion in Wahlberg et al. (2008) regarding the application to cascade structures. In Section 4, we point out that this method can be applied even if not all the transfer function outputs in a network are measured. Then, we consider a feedback cascade structure in Section 5, for which indirect PEM does not avoid non-convexity, and propose an intermediate step using the method in Galrinho et al. (2014). A numerical simulation is presented in Section 6, followed by a discussion in Section 7.
˜(t) the input measurewhere uo (t) is the true input, u ment noise, and u(t) the measured input. This setting corresponds to an errors-in-variables (EIV) problem (S¨oderstr¨om (2007)). In this scenario, we have that y(t) = ϕ⊤ (t)θ + v(t, θ) , v(t, θ) = y˜(t) − ϕ˜⊤ (t)θ where, if ϕo is defined analogously to ϕ, but containing true input values uo , then ϕ(t) ˜ = ϕ(t) − ϕo (t). Because v(t, θ) is not white, if the parameter vector θ is estimated according to (4), the obtained estimate is not consistent.
2. ERRORS-IN-VARIABLES METHODS Consider the SISO system G(q), and assume that data is generated according to yo (t) = G(q)u(t) , (2) y(t) = yo (t) + y˜(t)
Instrumental variable (IV) methods are appropriate to deal with EIV problems. The basic idea of IV methods is to choose a vector of instruments z(t) that is uncorrelated with the error v(t, θ), while being highly correlated with ϕ(t). Then, for such an instrument vector, computing −1 θˆ = ZΦ⊤ Zy,
where yo (t) is the true system output, y(t) is the measured output, corrupted by noise y˜(t), and the input u(t) is assumed to be known. We introduce the assumption that G(q) is FIR, and parametrize it accordingly as G(q, θ) = θ1 q −1 + θ2 q −1 + · · · + θn q −n , and that y˜(t) is Gaussian white noise.
where
Z = [z(1) z(2) . . . z(N )] , yields a consistent estimate of θ under certain excitation conditions.
The prediction error method (PEM) serves as benchmark in the field, since it is well known to provide asymptotically efficient estimates if the model orders are correct (Ljung (1999)). The essential idea of PEM is to minimize a cost function of the prediction errors. In this setting, PEM consists on minimizing the cost function N 1 2 V (θ) = (y(t) − G(q, θ)u(t)) , (3) N t=1
There is no unique way to define z(t). One approach, proposed in S¨oderstr¨om and Mahata (2002), is to choose z ⊤ (t) = [u(t − 1 − du ) . . . u(t − du − nzu )] , where du ≥ n. Another possibility, proposed in Thil et al. (2008), is to also include past outputs in the instrument vector, according to z ⊤ (t) = −y(t − 1 − dy ) . . . −y(t − dy − nzy ) (5) u(t − 1 − du ) . . . u(t − du − nzu )] , where dy is at least the order of the filter (it must be a moving average (MA) filter) applied to the noise. For the considered FIR case, dy ≥ 0.
if a quadratic cost is used, and where N is the number of samples available. Then, the minimizer of (3) is an asymptotically efficient estimate of θ, if the model orders are correct. In general, PEM requires solving a non-convex optimization problem. However, for this particular model structure, the minimizer of (3) can be obtained by solving a 99
2015 IFAC SYSID 100 October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
Returning to the cascade example of Fig. 1, it is possible to observe that such methods are appropriate for the estimation of G2 (q). Using the notation introduced in this Section, uo (t) plays the role of u2 (t), the unknown input to the second transfer function, and y1 (t) corresponds to u(t), the measurable input. Furthermore, the fact that the transfer functions are part of a network creates new possibilities for choices of instruments, as discussed in Dankers et al. (2014). In this case, a natural instrument candidate when estimating G2 (q) is the external input u(t). However, note that although G1 (q) can be obtained from data {u(t), y1 (t)} using any standard system identification method, discarding the output y2 (t) in the estimation of G1 (q) increases the variance of the estimated parameters (Everitt et al. (2013)). The contribution Gunes et al. (2014) addresses a similar issue for the two-stage method. Moreover, although G2 (q) can be estimated consistently from {y1 (t), y2 (t)} using the IV approach, it was shown in Hjalmarsson et al. (2011) that the EIV setting, even for an asymptotically efficient estimator, does not reach the asymptotic properties of PEM. 3. INDIRECT PEM In the cascade case, due to the product G2 (q)G1 (q), structures that are linear in the model parameters, such as FIR or ARX, cannot be directly applied, and PEM requires solving a non-convex optimization problem. Indirect PEM (S¨ oderstr¨ om et al. (1991)) offers, in some cases, a workaround for the non-convexity of PEM. The general idea is as follows. Consider two nested model structures M1 and M2 , such that M1 ⊂ M2 . Moreover, suppose that M1 is parametrized by the parameter vector θ, while M2 is parametrized by α. In a first step, α is estimated with PEM, by minimizing N V (α) = ǫ2 (t, α), (6) t=1
where ǫ(t, α) are the prediction errors associated with the model structure M2 . In some situations, M2 can be chosen such that minimizing (6) becomes an easy problem, even if applying PEM to the model structure M1 is difficult. For examples, see S¨ oderstr¨ om et al. (1991). Then, in a second step, an estimate of θ is obtained by minimizing 1 1 V (θ) = [ˆ α − α(θ)]⊤ Pα−1 [ˆ α − α(θ)] , (7) N2 where α ˆ is the estimate of α obtained by minimizing (6), and Pα is a consistent estimate of the covariance of α ˆ, obtained from minimizing (6). Note that (7) is still non-quadratic; however, it is rather easy to handle. Asymptotically, a single step of the GaussNewton algorithm, initialized with a consistent estimate of θ, provides an efficient estimate, i.e., it has the same asymptotic properties as if PEM had been applied to the model of interest directly. In S¨ oderstr¨ om et al. (1991), more details are given on the algorithm. Indirect PEM can be applied to the cascade system in Fig. 1 as follows. First, since the transfer functions considered are FIR, G1 (q) can be estimated from {u(t), y1 (t)},
100
and G21 (q) from {u(t), y2 (t)}, by solving LS problems as in (4). Since the parameters estimated in this step were, in the exposition above, designated by α, we shall write that estimates G1 (q, α) ˆ and G21 (q, α ˆ ) have been obtained in this step. In the second step, we solve G1 (q, θ) = G1 (q, α ˆ) (8a) ˆ (8b) G2 (q, θ)G1 (q, θ) = G21 (q, α) for θ in a weighted least squares (WLS) sense, where the weighting is a consistent estimate of the inverse covariance of α, obtained from the LS estimator by 1 Φ1 Φ ⊤ 0 1 −1 λ 1 Pα = , 1 ⊤ 0 λ2 Φ2 Φ2 and Φ1 and Φ2 are the regressor matrices for the estimation of G1 (q) and G21 (q), respectively. To be able to obtain an asymptotically efficient estimate in one step of the Gauss-Newton algorithm, consistent estimates of G1 (q, θ) and G2 (q, θ) are required. An estimate of G1 (q, θ) is readily available from (8a), while a consistent estimate of G2 (q, θ) can be obtained by applying LS to (8b), replacing G1 (q, θ) by G1 (q, α ˆ ). So, if the estimates obtained in this ˆ and G2 (q, θ), ˆ are used as initializations to way, G1 (q, θ) the Gauss-Newton algorithm, the next step yields asymp¯ totically efficient estimates, θ. Thus, the asymptotic efficiency achieved by applying indirect PEM is an advantage in comparison to the EIV approach. On the other hand, for other structures, such as output-error, IV techniques can still be applied to obtain consistent estimates of G1 (q) and G2 (q), while indirect PEM suffers from non-convexity. There are, however, certain cascade structures for which the standard IV approach is not applicable, but to which indirect PEM still applies. An example is discussed in the next section. 4. ANOTHER CASCADE STRUCTURE Consider the cascade structure in Fig. 2, which can be written as y1 (t) 0 G2 (q) u1 (t) e1 (t) , = + e2 (t) y2 (t) G3 (q)G1 (q) G3 (q)G2 (q) u2 (t) (9) where we maintain the FIR model assumption. u1 (t)
G1 (q)
u2 (t)
G2 (q) e1 (t)
G3 (q) e2 (t) y2 (t)
y1 (t) Fig. 2. Cascade system with three transfer functions. An important difference when comparing to the structure in Fig. 1, besides the presence of one more transfer function block and one more input, is that not all block inputs and outputs are measured. In particular, the output of G1 (q), as well as the input of G3 (q), are unknown. Therefore, the IV approach delineated in Section 2 is not directly applicable for estimating G1 (q) and G3 (q).
2015 IFAC SYSID October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
On the other hand, indirect PEM is still applicable. By using the parametrization G2 (q) = G2 (q, α) G3 (q)G1 (q) = G31 (q, α) , G3 (q)G2 (q) = G32 (q, α) we can use (9) to estimate G2 (q, α), G31 (q, α), and G32 (q, α). Then, the second step of indirect PEM concerns solving G2 (q, θ) = G2 (q, α) ˆ (10a) ˆ (10b) G3 (q, θ)G1 (q, θ) = G31 (q, α) ˆ (10c) G3 (q, θ)G2 (q, θ) = G32 (q, α) in the discussed WLS sense. Since (10a) can be used to ˆ (10c) can obtain a consistent estimate of G2 (q, θ), G2 (q, θ), then be used to obtain a consistent estimate of G3 (q, θ), ˆ With a consistent by replacing G2 (q, θ) with G2 (q, θ). ˆ we can use (10b) to obtain a consistent estimate G3 (q, θ), estimate of G1 (q, θ). With a consistent estimate θˆ of the complete vector θ, one step of the Gauss-Newton algorithm provides an asymptotically efficient estimate of θ. 5. CASCADE WITH FEEDBACK In this section, we consider a cascade system with two blocks connected in a feedback loop, represented in Fig. 3. Mathematically, we can write it as 1 y1 (t) G1 (q) G1 (q)G2 (q) = y2 (t) G2 (q) 1 − G1 (q)G2 (q) G1 (q)G2 (q) r1 (t) e (t) + 1 . (11) r2 (t) e2 (t) e1 (t) r1 (t)
y1 (t)
G1 (q)
y2 (t)
G2 (q)
r2 (t)
101
To do so, we use the following procedure. In a first step, we estimate high-order FIR models of each closed-loop SISO transfer function in (11). In a second step, we reduce these estimates to estimates of G1 (q), G2 (q), and the product G1 (q)G2 (q), with the desired orders, by LS. Then, we reestimate by WLS, where we use the estimate obtained in step 2 to construct the weighting matrix. Finally, the indirect PEM algorithm can be applied. In detail, the first step consists on approximating (11) by CL y1 (t) (q, g12 ) r1 (t) G11 (q, g11 ) GCL e1 (t) 12 = + , CL y2 (t) r2 (t) e2 (t) GCL 21 (q, g21 ) G22 (q, g22 ) (12) where m (k) CL gij q −k , Gij (q, gij ) = k=1
for i, j = {1, 2}. If the order m is chosen large enough, we have that 1 G1 (q) G1 (q)G2 (q) ≈ G2 (q) 1 − G1 (q)G2 (q) G1 (q)G2 (q) CL (13) G11 (q, g11 ) GCL (q, g12 ) 12 ≈ . CL GCL 21 (q, g21 ) G22 (q, g22 )
Since (12) is a multi-input multi-output (MIMO) FIR model, GCL ij (q, gij ) can be estimated by LS. This is done by solving, analogously to what was done from (2) to (4), −1 gˆ Φyi , gˆi := i1 = ΦΦ⊤ gˆi2 ⊤ Φ ⊤ = Φ⊤ r 1 Φr 2 , where Φr1 and Φr2 are the regression matrices for r1 (t) and r2 (t), respectively. The estimated parameters are, then, distributed as gˆi ∼ N (gi , Pgi ) , where the covariance matrix is given by −1 Pgi = λi ΦΦ⊤ . In a second step, we use the approximation (13) to write
e2 (t) Fig. 3. Block diagram of a feedback cascade system. Notice that this case is different from the typical feedback setting, with the measured output included in the loop. Instead, it can be seen as a natural extension of the previous cascade systems, where the systems are connected physically. A well-known historical example of a system physically connected in feedback is a centrifugal governor. Application of an IV method, in the EIV setting, to estimate G1 (q) and G2 (q) is straightforward, as both outputs are measured, and these are also the inputs to the other transfer function, summed with a known external reference. Note, however, that the same limitation of Section 4 would be encountered if, for example, output y2 (t) could not be measured. Concerning indirect PEM, it is difficult to find a model structure of finite order that is easy to estimate, and that contains the true model. The purpose of this section is to propose an intermediate step, based on the method presented in Galrinho et al. (2014), which reduces the problem to a setting that allows the application of the indirect PEM algorithm as in the previous sections. 101
Zii (q, α, gii ) :=GCL ii (q, gii ) (1 − G12 (q, α12 )) − Gi (q, αi ) = 0 Zij (q, α, gij ) :=GCL ij (q, gij ) (1 − G12 (q, α12 )) − G12 (q, α12 ) = 0
(14)
(in (14) only, i �= j), where we parametrize G1 (q) and G2 (q) as (1)
(ni ) −ni
Gi (q, αi ) = αi q −1 + · · · + αi and, to avoid the product G1 (q)G2 (q), (2)
q
(n +n −1)
.
G12 (q, α12 ) = α12 q −2 + · · · + α121 2 q −(n1 +n2 −1) . Then, if gij is replaced by its estimate, gˆij , (14) can be solved by LS. Write, first, (14) in vector form, as Q(ˆ g )α = gˆ, (15) where ⊤ ⊤ ⊤, α = α⊤ 1 α2 α12 ⊤ ⊤ ⊤ ⊤ ⊤ g = g11 g12 g21 g22 , , and obtain estimates α ˆ1, α ˆ 2 , and α ˆ 12 by computing the LS estimate of α from (15), i.e., −1 g )⊤ Q(ˆ g) Q(ˆ g )⊤ gˆ. α ˆ = Q(ˆ
2015 IFAC SYSID 102 October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
In a third step, we re-solve (14), now using WLS, where the weighting is obtained as follows. Let gˆij = gij + g˜ij . Then, we can write, in transfer function form, CL GCL ˜ij ). (16) ˆij ) = GCL ij (q, gij ) + Gij (q, g ij (q, g
100 98 96 94
˜ij ). Zij (q, α, gˆij ) = (1 − G12 (q, α12 )) GCL ij (q, g
(17)
Rewriting the equations defined in (17) in vector form, we have g˜11 z11 (α, gˆ11 ) g˜12 z12 (α, gˆ12 ) g = T (α12 ) =: T (α12 )˜ z(α, gˆ) := z21 (α, gˆ21 ) g˜21 z22 (α, gˆ22 ) g˜22 Then, the optimal weighting for the WLS is the inverse of the covariance of z(α, gˆ), i.e., �−1 � W (α12 ) = T (α12 )Pg˜ T (α12 )⊤ , where Pg˜ is the covariance of g˜, given by � � Pg1 0 . Pg˜ = 0 Pg2
Since the true α12 is not available, we use, instead, the estimate obtained in step 2, �−1 � W (ˆ α12 ) = T (ˆ α12 )Pg˜ T (ˆ α12 )⊤ , and compute � �−1 α ¯ = Q(ˆ g )⊤ W (ˆ α12 )Q(ˆ g) Q(ˆ g )⊤ W (ˆ α12 )ˆ g, which is an improved estimate of α. Asymptotically, we have that these estimates have covariance �−1 � α12 )Q(ˆ g) . (18) cov (α) ¯ = Q(ˆ g )⊤ W (ˆ
With these estimates, we can construct transfer functions ¯ 1 ), G2 (q, α ¯ 2 ), and G12 (q, α ¯ 12 ), and solve G1 (q, α G1 (q, θ) = G1 (q, α ¯1 ) (19a) ¯2 ) (19b) G2 (q, θ) = G2 (q, α ¯ 12 ) (19c) G1 (q, θ)G2 (q, θ) = G12 (q, α with the algorithm previously mentioned for indirect PEM, using (18) as covariance estimate for α ¯.
We note that, like indirect PEM for the acyclic cascade systems, this method does not require all the transfer function outputs to be measured. For example, if y2 (t) were not measured, the method would be applicable without major changes. In particular, (19b) would be absent in (19). 6. NUMERICAL SIMULATION In this section, we perform a numerical simulation to evaluate the method previously described for identification of the system in Fig. 3. One hundred Monte Carlo runs are performed, for sample sizes N = {100, 300, 600, 1000, 3000, 6000}. The systems are given by � (1) (2) G1 (q) = θ1 q −1 + θ1 q −2 . (1) (2) G2 (q) = θ2 q −1 + θ2 q −2 The closed-loop, whose poles are given by 1 − G1 (q)G2 (q) = 0, (20) is, thus, a fourth order system. The poles are placed randomly inside the unit circle, under the structure imposed 102
FIT
Replacing (16) in (14) yields
92 90 88 PEM WLSiPEM IVr IVuy
86 84 102
103
N
Fig. 4. Average FIT, for 100 Monte Carlo simulations, as function of sample size, comparing PEM, two IV methods, and the proposed WLS method combined with indirect PEM (WLSiPEM). by (20). The same colored input as in S¨oderstr¨ om and Mahata (2002) and Thil et al. (2008) is used, i.e., 1 rw (t), ri (t) = −1 1 − 0.2q + 0.5q −2 i for i = 1, 2, where rwi (t) is Gaussian white noise with unitary variance, and uncorrelated for each i. The sensor noises e1 (t) and e2 (t) are also Gaussian and white, and uncorrelated to each other. The variances are changed at each run in order to achieve a signal-to-noise ratio of 5dB at the outputs. The following methods are compared: PEM, as implemented in MATLAB2014b, starting at the true parameters (note that MATLAB does not provide an estimated point to initialize the algorithm, due to the non-standard model structure); the proposed method, referred to as WLSiPEM; the IV method with the external references as instruments (IVr); the IV method with instrument vector (5) (IVuy). Concerning WLSiPEM, the order m of the high-order FIR models is optimized over a grid of ten values, linearly spaced between m = 20 and m = N/3, and the one minimizing the prediction errors is chosen. A similar procedure is used to choose the length of the instrument vector for IV methods, by considering a grid of ten values linearly spaced between 4 and 22, and selecting the one minimizing the prediction errors. All the data available is used in the estimation. The accuracy of the estimates is measured by the FIT, in percent, defined as � � ˆ 2 ||θ − θ|| FIT = 100 1 − ¯ 2 , ||θ − θ|| where θ is a vector containing the true parameters, θ¯ its mean, and θˆ the estimated parameter. The results are presented in Fig. 4, with the average FIT as function of sample size. In this simulation, the presented method performed better than both IV methods. On the other hand, other simulations showed the IV methods
2015 IFAC SYSID October 19-21, 2015. Beijing, China
Miguel Galrinho et al. / IFAC-PapersOnLine 48-28 (2015) 098–103
to be more robust in the presence of extremely colored external excitations. It is also observed that using other variables in the network to construct the instrument vector improves the estimate, in comparison to using delayed measured input and output values. Finally, we verified that, in this case, the estimates obtained by any of the methods were appropriate as initial estimates for PEM, since the algorithm converged to the same parameters as when it was started at the true values.
7. DISCUSSION In this paper, we discussed two approaches for identification of network systems. One is based on the application of an IV method to an EIV setting, while the other uses indirect PEM. We point out that the applicability of each method depends on the particular problem setting and the assumptions on the network. The IV methods require that the input and output of all transfer functions are measured, while that is not a requirement for indirect PEM. Moreover, while the accuracy of IV methods depends on the choice of instruments and the correlation to the regression vector, indirect PEM guarantees asymptotic efficiency. On the other hand, standard IV methods admit rational transfer functions, while indirect PEM suffers from non-convexity in that case. If there is feedback in the network, indirect PEM does not avoid a non-convex problem. To deal with that limitation, we propose a method that uses the approach in Galrinho et al. (2014), and transforms the problem into one where the indirect PEM algorithm can be applied. Analogously to the acyclic case, measuring all outputs is not a necessary condition to obtain all the transfer functions. We perform a simulation study with two FIR transfer functions randomly generated and connected in feedback, where both outputs are measured, comparing PEM, the proposed method, and two IV methods with different instrument choices. The proposed approach performed better than the IV methods, but it did not prove as robust in the presence of extremely colored external excitations. This reinforces our discussion on the most appropriate approach depending on the settings and the network assumptions. To provide initial values for PEM, any of the techniques were successful. The noise sources were considered white, but most methods discussed allow the noise to be colored. For the IV method in Thil et al. (2008), the output noise must be given by a MA filter of white noise, although the input noise is assumed white. In the network context, the more general case, with an autoregressive moving average (ARMA) filter, is addressed in Dankers et al. (2014). Concerning indirect PEM, we still have a linear regression problem if the noise is given by an autoregressive filter driven by white noise. Finally, the proposed method is still applicable in the ARMA filter case, for which a high-order ARX model is estimated in the first step of the algorithm. Future work includes establishing conditions on the network structure for the methodology to be applicable, and extending the approach to the case where the blocks in the network are rational transfer functions. 103
103
REFERENCES Dankers, A., Van den Hof, P., and Heuberger, P. (2013). Predictor input selection for direct identification in dynamic networks. In 52nd Conference on Decision and Control, 4541–4546. Dankers, A., den Hof, P.M.V., Bombois, X., and Heuberger, P.S. (2014). Errors-in-variables identification in dynamic networks by an instrumental variable approach. In 19th IFAC World Congress, 2335–2340. Everitt, N., Rojas, C.R., and Hjalmarsson, H. (2013). A geometric approach to variance analysis of cascaded systems. In 52nd IEEE Conference on Decision and Control, 6496–6501. Everitt, N., Rojas, C.R., and Hjalmarsson, H. (2014). Variance results for parallel cascade serial systems. In 19th IFAC World Congress, 2317–2322. Galrinho, M., Rojas, C.R., and Hjalmarsson, H. (2014). A weighted least-squares method for parameter estimation in structured models. In 53rd IEEE Conference on Decision and Control, 3322–3327. Gunes, B., Dankers, A., and den Hof, P.M.V. (2014). Variance reduction for identification in dynamic networks. In 19th IFAC World Congress, 2842–2847. H¨agg, P., Wahlberg, B., and Sandberg, H. (2011). On identification of parallel cascade serial systems. In 18th IFAC World Congress, 9978–9983. Hjalmarsson, H., M˚ artensson, J., Rojas, C.R., and S¨oderstr¨om, T. (2011). On the accuracy in errorsin-variables identification compared to prediction-error identification. Automatica, 47(12), 2704–2712. Ljung, L. (1999). System Identification. Theory for the User, 2nd ed. Prentice-Hall. S¨oderstr¨om, T. (2007). Errors-in-variables methods in system identification. Automatica, 43(6), 939–958. S¨oderstr¨om, T. and Mahata, K. (2002). On instrumental variable and total least squares approaches for identification of noisy systems. International Journal of Control, 75(6), 381–389. S¨oderstr¨om, T., Stoica, P., and Friedlander, B. (1991). An Indirect Prediction Error Method for System Identification. Automatica, 27(1), 183–188. S¨oderstr¨om, T. and Stoica, P. (2002). Instrumental variable methods for system identification. Circuits, Systems and Signal Processing, 21(1), 1–9. Thil, S., Gilson, M., and Garnier, H. (2008). On instrumental variable-based methods for errors-in-variables model identification. In 17th IFAC World Congress, 426–431. Van den Hof, P.M., Dankers, A., Heuberger, P.S., and Bombois, X. (2013). Identification of dynamic models in complex networks with prediction error methods – Basic methods for consistent module estimates. Automatica, 49(10), 2994 – 3006. Van den Hof, P., Dankers, A., Heuberger, P., and Bombois, X. (2012). Identification in dynamic networks with known interconnection topology. In 51st Conference on Decision and Control, 895–900. Wahlberg, B. and Sandberg, H. (2008). Cascade structural model approximation of identified state space models. In 47th Conference on Decision and Control, 4198–4203. Wahlberg, B., Hjalmarsson, H., and M˚ artensson, J. (2008). On Identification of Cascade Systems. In Proc. 17th IFAC World Congress, 5036–5040.