Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME Convergence Properties of the Kalman Inverse Filter K. Soltani Naveh Graduate Student School of Mechanical and Mining Engineering The University of Queensland Brisbane, Australia Email: k.soltaninaveh@uq.edu.au sc rip tN ot Co P. R. McAree ? Professor School of Mechanical and Mining Engineering The University of Queensland Brisbane, Australia Email: p.mcaree@uq.edu.au py ed ite d B. S. Petschel Post-doctoral Research Fellow School of Mechanical and Mining Engineering The University of Queensland Brisbane, Australia ABSTRACT Ac ce pt ed Ma nu The Kalman filter has a long history of use in input deconvolution where it is desired to estimate structured inputs or disturbances to a plant from noisy output measurements. However, little attention has been given to the the convergence properties of the deconvolved signal, in particular the conditions needed to estimate inputs and disturbances with zero bias. The paper draws on ideas from linear systems theory to understand the convergence properties of the Kalman filter when used for input deconvolution. The main result of the paper is to show that, in general, unbiased estimation of inputs using a Kalman filter requires both an exact model of the plant and an internal model of the input signal. We show that for unbiased estimation, an identified sub-block of the Kalman filter that we term the plant model input generator must span all possible inputs to the plant and that the robustness of the estimator with respect to errors in model parameters depends on the eigen-structure of this sub-block. We give estimates of the bias on the estimated inputs/disturbances when the model is in error. The results of this paper provide insightful guidance in the design of Kalman filters for input deconvolution. ? Corresponding DS-13-1302 author. 1 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) by ASME List2017 of Figures 3 4 5 6 7 8 Kalman Inverse Filter block diagram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The Kalman Inverse Filter, including the open-loop system governed by A1 and A2 , and the closed-loop feedback system governed by A?1 and A3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The double mass-spring-damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Response of PIG and PMIG for step disturbance in double mass-spring-damper system. . . . . . . . . . . . Response of PIG and PMIG for combined sinusoids disturbance in double mass-spring-damper system. . . The block diagram of the mass spring damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper system with imperfect model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bias and bias upper bound for imperfect plant model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ed ite d 1 2 3 4 18 19 21 23 26 27 Co py List of Tables 1 Parameters of the double mass-spring-damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2 Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3. 24 Introduction Inverse problems in systems theory involve reconstructing the inputs to a system from measured responses. The adaptation of the Kalman filter to inverse problems, sometimes termed the Kalman inverse filter (KIF), was first proposed by Bayless and Brigham [1] as a framework for removing sensor reverberation from seismic signals. The idea is to augment the states of the plant model by additional dynamics that captures known structure in the inputs to be estimated. This is equivalent to assuming the autocorrelation function of the unknown inputs or distrubances that are to be estimated. The formulation of the KIF in [1] was for continuous-time systems. A discrete-time formulation followed shortly after in [2]. Most of the applications of the KIF over the last thirty years have been for geophysical signal processing, see for example [3?5], but the method has also attracted interest for input and disturbance identification in mechanical systems. Applications include estimation of aerodynamic loads on aircraft [6, 7], friction in machines [8?10], and the estimation of machine-environment interaction forces [11?14]. The KIF is closely linked to other methods for input identification, e.g. the regularization method of Tikhonov [15]. Work by De Nicolao and Ferrari-Trecate [16, 17] has shown that under appropriate assumptions, Tikhonov regularization coincides with Bayesian estimation where the autocorrelation function of the unknown input is specified and the plant is a linear system driven by white noise. De Nicolao and Ferrari-Trecate [16] show that the inverse problem can be solved by application of the Kalman filter, revealing an equivalence between Tikhonov regularization and the KIF. As a method for input identification the KIF has a number of attractions including: (i) that it is recursive and therefore suitable for real-time implementations; (ii) the underlying Kalman filter structure accommodates model and measurement uncertainty; and (iii) the method can be adapted to non-linear systems through the extended Kalman filter. Limitations of the method are: (i) that the structure of the inputs/disturbances must be specified, and (ii) the bandwidth to which input/disturbances can be reconstructed depends on the dynamics of the plant, and the process and measurement uncertainties, as specified in the Kalman filter. The structure of the KIF is illustrated in Fig. 1, and has four components: the Plant, the Plant Input Generator (PIG), the Plant Model and the Plant Model Input Generator (PMIG). The Plant is the system whose inputs are to be estimated. The PIG represents the process that generates these to-be-estimated inputs/disturbances acting on plant. The KIF comprises the Plant Model and the PMIG. The PMIG assumes the inputs have structure. The idea is to design the PMIG so that its outputs serve as estimates of the input to the Plant. This is done through the feedback structure of the Kalman filter which provides an implicit inverse model of the Plant. This paper establishes the conditions needed for convergence of the PMIG output to the Plant input in steady state for linear plants. Specifically, the paper establishes the conditions that must be satisfied by the Plant Model and the PMIG for convergence of the plant input/disturbance estimates to their actual values. Section 2 describes the structure of the Kalman Ac ce pt ed Ma nu sc rip tN ot 1 DS-13-1302 2 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME real system Plant Input Generator (PIG) Plant plant inputs (unknown) + estimation error ed ite d estimated plant inputs system model/estimator - Plant Model Input Generator (PMIG) + py Plant Model Co Kalman gain ot Fig. 1: Kalman Inverse Filter block diagram. The Kalman Inverse Filter The KIF is described by the block diagram in Fig. 2 and governed by the equations Ac ce pt ed 2 Ma nu sc rip tN Inverse Filter and details the notation used throughout the paper. Section 3 summarizes key assumptions. The main work of the paper begins in Section 4 where we give general conditions for an asymptotically unbiased KIF by formulating the problem in terms of the solution to a constrained linear matrix equation. It is shown that if the model of the Plant and PMIG are exact, the estimates produced by the KIF converge to their true values. Section 5 presents a more detailed investigation of the requirements on the PMIG and shows that for unbiased estimation it must span all possible inputs to the plant and that the robustness of the estimator with respect to errors in parameters depends on the eigen-structure of the PIG. It is shown, moreover, that a robust PMIG can be constructed from knowledge of just the number of outputs and the minimal polynomial of the PIG dynamics. Inter alia, the analysis also yields a simple proof of the Internal Model Principle which was first proven in [18]. Section 6 gives error bounds on estimates given knowledge of model errors. The results of this section show that the estimated bias is approximately proportional to the parameter error and that relative bias is bounded. DS-13-1302 x?1 = A1 x1 + B1 u + w1 (1a) x?2 = A2 x2 + w2 (1b) y = C1 x1 + v (1c) u = C2 x2 x?? 1 = A?1 x?1 + B?1 u? + Kb e (1d) x?3 = A3 x3 + Kc e (1f) (1e) y? = C?1 x?1 (1g) u? = C3 x3 (1h) e = y ? y? (1i) z = u ? u? (1j) 3 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME where thebyvector x is the state of the Plant, x is the state of the Plant Input Generator (PIG), u is the unknown plant input 1 2 produced by the PIG and y is the measured output signal from the plant. The vector x?1 is the state of the plant model, u? is the estimate of the input which is produced by the Plant Model Input Generator (PMIG) with state vector x3 generally not the same size as x2 , and e is the innovation sequence available for feedback. The matrices A1 , A2 , B1 , C1 , C2 , A?1 , B?1 , C?1 , A3 , C3 have appropriate sizes satisfying the Eqns. (1). In particular the PMIG and the PIG have the same number of outputs and the Plant and the Plant Model have the same number of outputs. Plant PMIG Plant Model Ma nu sc rip tN ot Co py ed ite d PIG ed Fig. 2: The Kalman Inverse Filter, including the open-loop system governed by A1 and A2 , and the closed-loop feedback system governed by A?1 and A3 . Ac ce pt We note, however, that the Plant Model and PMIG may contain additional dynamics not present in the Plant and PIG under investigation or vice-versa. Hence the model state vectors x?1 , x3 may be of different sizes to the corresponding plant w (t) state vectors x1 , x2 . The measurement noise v and the process noise w(t) = w1 (t) are vector white noises with covariance 2 structure E[w(s)wT (t)] = Q?(s ? t) E[v(s)vT (t)] = R?(s ? t) E[v(s)wT (t)] = 0. The system can be thought of as regulating z, the difference between the input signal and the estimated input signal. When E[z] = 0 the estimates u? track the actual inputs to the Plant u. DS-13-1302 4 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME The by state vector of the combined system encompassing the Plant and the PIG is " # x1 x= . x2 " # x?1 x? = , x3 which shall form the model state. The state equations (1) combine to form the system y = CP x (3b) (3c) (3d) Co y? = CM x? (3a) py x? = AP x + w, x?? = AL x? + BP x + Kv, ed ite d The state vector of the model of the system, encompassing the Plant Model and the PMIG, is (3e) ot z = DP x ? DM x?, tN where " # A?1 B?1C3 AM = 0 A3 sc rip " # A1 B1C2 AP = , 0 A2 AL = AM ? KCM , h i CP = C1 0 , h i DP = 0 C2 , " # Kb K= , Kc BP = KCP h i CM = C?1 0 , h i DM = 0 C3 Ma nu (4) T ?1 K = PCM R pt ed The so-called Kalman gain K is given by T ?1 R CM P = 0. AM P + PATM + Q ? PCM Ac ce where P satisfies the Algebraic Riccati Equation [19] In order to regulate z, it is necessary to have closed loop stability, that is ?(AL ) ? C? , where ?(AL ) is the eigenvalue spectrum of the matrix AL and C? is the open left half of the complex plane. A necessary condition for existence of a stabilizing gain K is that the pair (AM ,CM ) is detectable. We assume that the pair (AM ,CM ) is observable; this guarantees existence of a stabilizing gain K such that AP and AL have no common eigenvalues. DS-13-1302 5 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) by ASME 3 2017 Assumptions The following assumptions are made at various places and referred to by number as required. They are gathered here to avoid repetition. As indicated, some of the assumptions may be replaced with less restrictive alternatives. ed ite d A1 The plant system is observable, i.e. (AP ,CP ) is observable. A2a All eigenvalues of the plant system matrix AP are unstable. A2b The plant system matrix AP and the closed loop system matrix AL have no common eigenvalues. A3 The system model is observable, i.e. (AM ,CM ) is observable. A4 All eigenvalues of the closed loop matrix AL are stable. A5a The Plant Model exactly describes the Plant, i.e. (A?1 , B?1 , C?1 ) = (A1 , B1 ,C1 ). A5b The Plant Model has the same transfer function as the Plant. A6 The Plant matrix A1 has no eigenvalues in common with the Plant Input Generator matrix A2 . py Assumption A4 is essential to the results of this paper. Assumption A3 implies the existence of stabilizing gains K such that Assumption A4 is satisfied. Convergence properties of the KIF This section gives general conditions for an asymptotically unbiased KIF, namely E[z(?)] = 0. The approach taken is to reformulate the problem in terms of the solution to a constrained linear matrix equation by decomposing the trajectories into terms arising from the state and the state estimation error. Then using Assumption A4, the conditions for E[z(?)] = 0 are derived. tN General conditions for unbiased input estimation Consider the state equations (3), then x and x? are related by sc rip 4.1 ot Co 4 x? = ? + Xx (5) Ma nu where X has appropriate dimensions and ? represents the state estimation error. Differentiating both sides with respect to time and using state equations (3) we can write ??? = x?? ? X x? = AL x? + BP x + Kv ? X(AP x + w) = AL (?? + Xx) + BP x + Kv ? X(AP x + w) ed = AL ? + (AL X + BP ? XAP )x + Kv ? Xw. XAP ? AL X = BP (6) Ac ce pt Under Assumption A4, AL is stable and hence AL ? converges to zero as t ? ?. Therefore, in order to satisfy limt?? ? = 0, the term due to x must be zero, so must hold. This relation is the well-known Sylvester matrix equation and has unique solution X under Assumption A2b that is if AP and AL have no common eigenvalues [20?22]. Assumption A2b can be satisfied in a number of ways. One sufficient condition is that Assumptions A2a and A4 hold. Using the above relations, the input estimation error z can then be expressed in terms of the state and state estimation error. Specifically z = DP x ? DM x? = (DP ? DM X)x ? DM ? . DS-13-1302 6 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME Since E[??by ]? 0 the term to due to ? converges to 0 . Therefore in order for E[z] ? 0 , the component due to x must be zero, leading to the following results: Lemma 1 (Sufficient conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if DM X = DP . (7) where X is the solution to the matrix equation (6). ed ite d Lemma 2 (Necessary conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if and only if (DP ? DM X)U = 0 (8) where X is the solution to the matrix equation (6) and U is a matrix whose columns span the unstable subspace of AP . Co py Note that when all eigenvalues of AP are stable, U = 0 and (8) holds regardless of the accuracy of the model. The KIF is asymptotically unbiased because both the input signal and the input estimate converge to zero. Conversely, if all eigenvalues of AP are unstable (Assumption A2a), U = I in (8) which reduces to to (7). Exact Plant Model, Exact PMIG Before discussing the question of robustness with respect to parameters, we remark that for an exact model the bias is zero, i.e. for tN ot 4.2 sc rip (A?1 , A3 , B?1 , C?1 ,C3 ) = (A1 , A2 , B1 ,C1 ,C2 ), A?1 = T1 A1 T1?1 , A3 = T2 A2 T2?1 , Ma nu the equations (6) and (7) are satisfied by X = I. In fact, the choice of coordinates for the model state space does not affect the bias: given invertible matrices T1 and T2 , the model described by C?1 = C1 T1?1 , B?1 = T1 B1 , C3 = T1C2 T2?1 , " # T1 0 X= . 0 T2 pt ed has zero bias, with equations (6) and (7) satisfied by Ac ce Here the converged state estimates are related to the true states by x?1 = T1 x1 x3 = T2 x2 . 4.3 Zero-bias equations To show that for unbiased input estimation, the plant model must be similar to the plant, partition the matrix X into four blocks of compatible size, with " # X11 X12 X= X21 X22 DS-13-1302 7 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) by equation ASME (6) becomes The2017 matrix X11 A1 = (A?1 ? KbC?1 )X11 + B?1C3 X21 + KbC1 (9a) X21 A1 = ?KcC?1 X11 + A3 X21 + KcC1 (9b) X11 B1C2 + X12 A2 = (A?1 ? KbC?1 )X12 + B?1C3 X22 (9c) X21 B1C2 + X22 A2 = ?KcC?1 X12 + A3 X22 , (9d) ed ite d and the zero bias equation (7) becomes C3 X21 = 0 C3 X22 = C2 . (10a) (10b) Co py Any given Plant and PIG has a very large class of models for which (9) and (10) can be solved for X. We limit our attention to models for which the Plant Model has the same dimensionality as the true Plant. In this situation the Plant and the Plant Model must necessarily obey a similarity relationship. The argument is as follows. After substitution of (10a), (9a) can be rearranged as ot A?1 X11 ? X11 A1 = Kb (C?1 X11 ?C1 ). tN In order for this equation to hold regardless of the feedback gain Kb , the left and right sides must be zero, giving A?1 X11 = X11 A1 , sc rip C?1 X11 = C1 . Similarly, (9c) can be reduced to Ma nu (X11 B1 ? B?1 )C2 = (A?1 ? KbC?1 )X12 ? X12 A2 , which in order to hold regardless of the signal generator gain C2 necessitates ed X11 B1 ? B?1 = 0. pt Since the plant (A1 ,C1 ) is observable and the model has the same number of dimensions, X11 must be square, invertible and unique (see Sec. 5), leading to the following result: Ac ce Theorem 1. Suppose the Plant (A1 , B1 ,C1 ) is observable and controllable, and that the associated zero bias Plant Model (A?1 , B?1 , C?1 ) has the same number of state space dimensions. Then the Plant and the Plant Model must be related by a unique similarity transform ?1 A?1 = X11 A1 X11 , ?1 C?1 = C1 X11 , B?1 = X11 B1 . Moreover the Plant and the Plant Model have the same transfer function C?1 (s ? A?1 )?1 B?1 = C1 (s ? A1 )?1 B1 for all s ? C. DS-13-1302 8 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME Thisby result implies a requirement for exact knowledge of the Plant transfer function: to exactly determine the plant inputs from the plant outputs requires the exact relationship between the inputs and the outputs. In general this requirement implies a need for a perfect Plant Model including exact parameter knowledge. An interesting exception is where all eigen-values of A1 are zero corresponding, for example, to free rigid-body dynamics problems. In this situation the eigen-structure is independent of the parameter values. The model used by Siegrist [13] to estimate road tyre interaction forces for off-highway mining vehicles is in this class of problems. Conditions for unbiased input estimation with a perfect plant model For algebraic simplicity, we will henceforth assume a perfect Plant Model with A?1 = A1 , B?1 = B1 , C?1 = C1 so that X11 = I. We now seek conditions for unbiased input estimation under this assumption. A sufficient condition for zero bias is the existence of a solution to the constrained Sylvester equation ed ite d 4.4 C3 X22 = C2 (11b) py A3 X22 = X22 A2 . (11a) Co The solvability of (11) becomes a necessary condition for zero bias when the signal generator has no modes in common with the plant, as demonstrated in the following result. tN ot Theorem 2. Under Assumptions A2b, A4, A5a, suppose that the PMIG matrix A3 has no eigenvalues in common with the Plant Model matrix A?1 and that Kb has been chosen so that (A?1 ? KbC?1 ) has no eigenvalues in common with A1 or A2 (this is possible because (A1 ,C1 ) is observable). The KIF has zero bias if and only if there exists a solution X22 to the simultaneous matrix equations (11). sc rip Proof. To prove sufficiency, let X22 be a solution to (11). Then the equations (9), (10) are uniquely satisfied by " # I 0 X= . 0 X22 Ma nu To prove necessity, suppose the KIF has zero bias, i.e. (10) holds. Then (9a) reduces to X11 A1 = (A?1 ? KbC?1 )X11 + KbC1 . X21 A1 = A3 X21 X12 A2 = (A?1 ? KbC?1 )X12 ce pt ed Since the plant model (A?1 , C?1 ) is exact and A1 has no eigenvalues in common with (A?1 ? KbC?1 ), the solution X11 = I is unique. The next two equations (9b) and (9c) reduce to Ac which, by the assumption of no common eigenvalues, are uniquely satisfied by X21 = 0 and X12 = 0. The other equation (9d) reduces to X22 A2 = A3 X22 as required. We now show that equations (11) imply a relationship between the observability matrices of the PMIG and the PIG. Suppose that equations (11) hold for the matrix X22 . Multiply (11a) on the right by A2 and substitute (11b) to obtain C3 A3 X22 = C2 A2 . DS-13-1302 9 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 bythe ASME Repeating process will yield a sequence of equations j j j ? 0. C3 A3 X22 = C2 A2 , Terminating the sequence at j = k ? 1,1 the system can be expressed as Ok (A3 ,C3 )X22 = Ok (A2 ,C2 ), (12) ed ite d where Ok (A,C) is the kth order observability matrix defined by ? C ? ? ? CA ? ? Ok (A,C) ? ? . ? ?. ? .. ? CAk?1 py ? Co When k = n, where A is n О n, the matrix O (A,C) is the observability matrix of control theory. The noise-free component of an input signal produced by state x2 of the signal generator satisfies ot u = C2 x2 , tN u? = C2 x?2 = C2 A2 x2 , u? = C2 x?2 = C2 A22 x2 , sc rip .. . du dx2 = C2 k?1 = C2 Ak?1 2 x2 . dt k?1 dt Ma nu Similarly the input estimate produced by state x3 satisfies u? = C3 x3 , u?? = C3 x?3 = C3 A3 x3 , pt ed u?е = C3 x?3 = C3 A23 x3 , .. . d u? dx3 = C3 k?1 = C3 Ak?1 3 x3 . dt k?1 dt ce These relations can be summarized in two systems of equations u u? u? .. . Ac ? ? ? ? ? ? ? ? ? 1 As du dt k?1 ? ? ? ? ? ? ? = O (A2 ,C2 )x2 , ? ? ? ? ? ? ? ? ? ? ? u? u?? u?е .. . d u? dt k?1 ? ? ? ? ? ? = O (A3 ,C3 )x3 . ? ? ? a consequence of the Cayley-Hamilton theorem of linear algebra, the sequence can be stopped at k = n2 n3 , where A2 is n2 О n2 and A3 is n3 О n3 ; in the next section we show that k = n3 + 1 suffices. DS-13-1302 10 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017isby ASMEand x = X x , the PIG produces an identical signal as shown by the following relation. If (12) satisfied 3 22 2 ? ? ? ? ? ? ? ? ? u u? u? .. . du dt k?1 ? ? ? ? ? ? ? ? ? ? ? = O (A3 ,C3 )X22 x2 = O (A3 ,C3 )x3 = ? ? ? ? ? ? ? u? u?? u?е .. . d u? dt k?1 ? ? ? ? ? ?. ? ? ? tN ot Co py ed ite d The interpretation we give to this result is that the signals generated by the PMIG must span all possible inputs to the plant, so that if, for example, the plant input contains a sinusoidal component of frequency ?, the PMIG must include a mode with this same frequency. Furthermore, the matrix X22 is a map from the states of the PIG to the states of the PMIG so that the inputs to the model correspond to those of the plant in a way that is consistent with the dynamics of the two systems. The above results provide a test of whether a given PMIG will produce unbiased input estimates for a given PIG. No assumptions have been made about the size of the PMIG relative to the PIG. The PMIG may have additional signal components that do not appear in the PIG, or multiple copies of components of the signal generator. The PMIG can therefore contain a library of modes whose form corresponding to all expected inputs to be applied to the plant. By way of example, if it is known that the plant can subjected to inputs/disturbances having n distinct frequencies, ?1 , . . ., ?n , applied alone or in combination, the PMIG should be designed to have eigenvalues corresponding to those frequencies. The next section looks in more detail at the relationship between the structure of the signal generator and the structure of models that produce zero-bias signal estimates. The analysis shows that the robustness of the estimator with respect to certain parameters depends largely on the eigen-structure of the PIG. This amounts to determining the general structure of pairs of matrices (A3 ,C3 ) and (A2 ,C2 ) which satisfy Equation (11) for some X22 . Conditions on the PMIG for zero bias estimates In this section, we consider the abstract properties of the relationship between the PIG and the PMIG implied by the zero bias equations (11), culminating in the result that a robust PMIG can be constructed from knowledge of just the number of outputs and the minimal polynomial of the PIG dynamics. For ease of argument, we introduce the following notation: sc rip 5 Ma nu Definition 1. For any matrices A,C, B, D of size n О n, k О n, m О m, k О m respectively, say (A,C) models (B, D), denoted (A,C) & (B, D), if (13a) AX = XB, (13b) ed CX = D pt for some matrix X of size n О m. Say (A,C) &X1 (B, D) if (11) holds for X = X1 . In particular, the PMIG will provide zero bias estimates for a given PIG if (A3 ,C3 ) & (A2 ,C2 ). ce The relation has useful ordering properties: Ac Reflexive: (A,C) &In (A,C), where In is the n О n identity matrix. This means any plant is a model of itself. Transitive: (A1 ,C1 ) &X1 (A2 ,C2 ) &X2 (A3 ,C3 ) implies (A1 ,C1 ) &X1 X2 (A3 ,C3 ). This means any model of a model of the plant is, in turn, a model of the plant. From these properties, two equivalence relations naturally arise: Similarity: Say (A,C) ? (B, D) if (A,C) &X (B, D) where the matrix X is square and invertible. This is similarity in the usual linear systems sense. Equivalence: Say (A,C) h (B, D) if (A,C) & (B, D) & (A,C). Similarity implies equivalence, but not necessarily vice-versa. Pairs of models can be related in three fundamental ways (subject, of course, to compatibility of the size of the sub blocks): DS-13-1302 11 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) by ASME R12017 (change of basis): " (P?1 AP,CP) &! (Q?1 AQ,CQ) by X"=#P?1 Q # h i A B1 I R2 (bigger model): , C D & (A,C) by X = 0 B2 0 " # ! h i h i A 0 R3 (unobservable signals): (A,C) & , C0 by X = I 0 B3 B4 ed ite d Relation R1 says that the choice of coordinate system is irrelevant, relation R2 says that the PMIG can contain additional dynamics not present in the PIG. Relation R3 says that unobservable components of the PIG can be ignored. In fact these relations characterize &, since any matrix X that relates (A,C) and (B, D) can be reduced by row and column operations to echelon form by P?1 XQ = 0I 00 for some invertible matrices P and Q2 . This is summarized in the following result: Theorem 3. (A,C) & (B, D) if and only if there are change of basis matrices P and Q that reduce (A,C) and (B, D) to upperand lower- block triangular form respectively such that (Q?1 BQ, DQ) = # ! i A12 h , C1 C2 A21 # ! h i 0 , C1 0 . B22 py (P AP,CP) = " A11 0 " A11 B21 Co ?1 Ma nu Lemma 3. If (A,C) &X (B, D), then sc rip tN ot The above result characterizes the perfect plant model as containing a subsystem that is equivalent to the observable part of the plant. In terms of structural stability, this says that the model is insensitive to perturbations that do not change the equivalence of this subsystem that models the plant. The aim of the remainder of this section is to show that any observable model (A,C) with suitable replication of the eigenstructure of B is a perfect model of (B, D), i.e. contains a subsystem equivalent to (B, D). The implication of this result is that the model is insensitive to any perturbations that do not alter the eigenstructure or compromise the observability of the model. In order to achieve this result, we re-examine the results relating the observability matrices of the systems. In the previous section, we demonstrated that (11) implies (12), and hence Ok (A,C) ? X = Ok (B, D) for all k ? 1. (14) pt ed The reverse is not necessarily true, however note that (14) implies CAi (AX ? XB) = 0 for i ? k ? 2, leading to the following result: On (A,C) ? (AX ? XB) = 0. (15) Ac ce Lemma 4. If (14) holds for k = n + 1, then (AX ? XB) satisfies Moreover if (A,C) is observable (i.e. On (A,C) has rank n), AX ? XB = 0 and hence (A,C) & (B, D). Several results immediately follow: Corollary 3.1. 2 Note that some of the sub blocks may have zero size. DS-13-1302 12 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME If (A,C) & (B, D) and (A,C) is observable, such an X is unique. X If (A,C) &X (B, D) and (A,C) is observable, any matrices E, F satisfying E ? X = F are of the form E = Q ? O (A,C), F = Q ? O (B, D) for some matrix Q. If (A,C) &X (B, D) and (B, D) is observable, then X has full column rank. If (A,C) &X (B, D) and (B, D) is observable and n = m, then X is invertible (hence (A,C) ? (B, D) and (A,C) is observable). If (A,C) h (B, D) and both are observable, then n = m and (A,C) ? (B, D). ed ite d It is well known that any two single output systems (A, cT ), (B, dT ) are equivalent if they are observable and A and B are similar (e.g. Casti [23]). A slightly more general result follows: Corollary 3.2. Suppose A ? B, and cT , dT are any two row vectors such that (A, cT ) is observable (single-output case). Then (A, cT ) &X (B, dT ) with X = (On (A, cT ))?1 On (B, dT ). py Proof. Since (A, cT ) is observable, On (A, cT ) is n О n and invertible. Since A and B have the same characteristic polynomial, if (14) holds for k = n it must hold for k = n + 1 (the final equation is the same linear combination of the previous n equations) and the result follows. ot Co The following result shows that every model of an output system contains a subsystem equivalent to that system. " # ! A1 , A2 Lemma 5. Suppose (A,C) &X (B, D) where (B, D) is observable. Then (A,C) ? , [C1 ,C2 ] for some A1 , A2 , A4 ,C1 ,C2 0, A4 such that (A1 ,C1 ) ? (B, D). tN Proof. Let (A,C) &X (B, D). Since (B, D) is observable, X has full column rank by corollary 3.1. Then sc rip " # I Q?1 , X =P 0 where P and Q are invertible matrices. Hence by Eq.(13) Ma nu " # ! i A1 , A2 ?1 h P P , C1 , C2 P?1 , 0, A4 (A,C) = where A1 = Q?1 BQ and C1 = DQ, and the result follows " # A1 0 A1 ? A2 = , 0 A2 ce pt ed Larger systems can be built up from sums of smaller models in parallel: define the direct sum of two matrices as and we have the following result: Ac Lemma 6. If (A1 ,C1 ) & (B1 , D1 ) and (A2 ,C2 ) & (B2 , D2 ), then (A1 ? A2 , [C1 ,C2 ]) & (B1 ? B2 , [D1 , D2 ]) (provided the matrices C1 and C2 have the same number of rows) and (A1 ? A2 ,C1 ?C2 ) & (B1 ? B2 , D1 ? D2 ). A generalization of this result to k-output systems forms a basis for the internal model principle: Lemma 7. Suppose A is in block-diagonal form with k identical blocks, i.e. ? ? A 0 1 k M ? . ? .. ?, A= A1 = ? ? ? 0 A1 DS-13-1302 13 (16) McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 where A by is ASME n О n , and 1 1 1 ? ? cT1 0 ? ? . .. ?, C=? ? ? T 0 ck ? dT11 ? . . D=? ? . dTk1 ? dT1k .. ? ? . ?, и и и dTkk иии .. . ed ite d where cT1 , . . . , cTk , dT11 , . . . , dTkk are row vectors. Then (A,C) & (A, D) if and only if (A1 , cTi ) & (A1 , dTij ) for all i, j, 1 ? i ? k, 1 ? j ? k. Proof. Setting Co where the sub-blocks Xi j are each of size n1 О n1 , (13) reduces to tN ot A1 Xi j = Xi j A1 cTi Xi j = dTij , py ? ? X11 и и и X1k ? . . ? . . . ... ? , X =? ? . ? Xk1 и и и Xkk and the result follows. sc rip We note in Lemma 7 that observability of all subsystems (A1 , cTi ) is equivalent to observability of (A,C). Also by Corollary 3.1, if (A, D) is observable then the systems are equivalent, so we have the following result: Ma nu Lemma 8. Suppose (A,C1 ), (A,C2 ) are systems where A is of block-diagonal form (16), and C1 ,C2 are k О n matrices. If (A,C1 ) is observable, (A,C1 ) & (A,C2 ). Moreover, if both systems are observable, they are equivalent. Proof. All that remains is to show that if (A,C1 ) is observable, there exists a vector cT1 such that (A1 , cT1 ) is observable. Since A consists of k copies of A1 , they have the same minimal polynomial and hence ed rank(On (A,C)) ? km1 , Ac ce pt where m1 is the degree of the minimal polynomial of A1 . The matrix has full rank only if m1 = n1 , i.e. A1 is cyclic (also termed nonderogatory [21]: each eigenvalue has geometric multiplicity 1). For some change of basis matrix P, such a matrix will have Jordan form ? J1 ? ?1 P A1 P = J = ? ?0 0 ? 0 0 .. ? . 0? ?, 0 Jr where the Jordan blocks ? ? ?i 1 0 ? . ? ? . Ji = ? ?0 . 0? 0 0 ?i DS-13-1302 14 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by to ASME correspond distinct eigenvalues ? . Set dT = [1, 0, . . . , 0]; it easily follows that (J , dT ) is observable, as is (J, dT ) where i i i i dT = [dT1 , . . . , dTk ]. Then set cT1 = dT P?1 , so we have (A1 , c1 ) ? (J, dT ) and hence (A1 , c1 ) is observable. With this choice of c1 , apply Lemma 7 with cTi = cT1 , 2 ? i ? k to obtain C such that (A,C) & (A,C1 ) and (A,C) & (A,C2 ). If (A,C1 ) is observable, (A,C1 ) ? (A,C) and hence (A,C1 ) & (A,C2 ). If both (A,C1 ) and (A,C2 ) are observable, they are both equivalent to (A,C), and hence (A,C1 ) ? (A,C2 ). This result will be used to show that any system (B, D) can be perfectly modelled by an observable system (A,C), where A is of the form (16). The matrix A1 is any matrix consisting of the largest cyclic component of B, i.e. a cyclic matrix satisfying the minimal polynomial of B. Lk A1 where A1 is the largest cyclic component of B. Then (A,C) & (B, D) ed ite d Lemma 9. Suppose (B, D) is observable, and A = if (A,C) is observable. Lp " # Bi Bi1 B1 ? 0 Bi2 ot for some matrices Bi1 , Bi2 , and hence we can augment the model to obtain Co py Proof. Without loss of generality, assume A1 is in Jordan form. By the cyclic decomposition theorem [20, 21], B ? i=1 Bi where the Bi are cyclic matrices in Jordan form such that B1 = A1 and each Bi satisfies the minimal polynomial of B1 . Each Bi can be augmented to a matrix similar to B1 , i.e. tN " # ! i B B01 h , D 0 & (B, D), 0 B02 sc rip such that " # p M B B01 B1 ? 0 B02 i=1 Ma nu and B01 = ed B02 = p M i=1 p M Bi1 Bi2 . i=1 ce pt Since (B, D) is observable, p ? k,3 and so the model can be augmented with k ? p copies of B1 to obtain Ac (A,C) & " # ! ! k?p h i M B B01 ? B1 , D 0 & (B, D), 0 B02 i=1 and the result follows. As an easy corollary of the previous result and relations R1 and R2, the main result follows, namely that any observable system that models k copies of the largest cyclic component of B, is a perfect model of (B, D), where k is the number of rows of C and D. This gives an alternative (and arguably simpler) proof of the Internal Model Principle [18]. 3 PROOF: Suppose p > k. Some eigenvalue ? has geometric multiplicity p, so rank B??I D ? n + k ? p < n, and hence (B, D) is not observable. Hence p ? k. DS-13-1302 15 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by4 ASME Theorem (Internal Model Principle). Suppose D has k rows, B is the largest cyclic component of B, and 1 "L # B1 A1 , 0 A2 k A? (17) for some matrices A1 and A2 . Then (A,C) models (B, D) for any matrix C such that (A,C) is observable. "L ed ite d Proof. Let X be a change of basis matrix such that X ?1 AX = A0 , where A0 is the right hand side of (17), and so # ! i B1 A1 h , C1 C2 , 0 A2 k (A,C) ? where [C1 ,C2 ] = CX. By relations R1 and R2, (A,C) & ( k B1 ,C1 ). Since (A,C) is observable, the subsystem ( k B1 ,C1 ) is observable, and by the previous lemma, models (B, D) for any D of compatible size. Hence, by the transitive property, (A,C) & (B, D). L Co py L sc rip tN ot We have proved that a robust model of an observable plant (and in particular, a robust PMIG for a given PIG) can be constructed from knowledge of only the minimal polynomial (i.e. the largest cyclic component) of the plant dynamics, and the number of outputs, by incorporating multiple copies of the largest cyclic component. This model is robust to any perturbations that do not affect the eigenstructure of the model, or the observability of the model, and moreover it may contain additional dynamics not present in the underlying plant. We now prove the converse of the internal model principle, namely that any ?universal? model of an observable plant necessarily contains a subsystem with multiple copies of the largest cyclic component. Theorem 5. Suppose (A,C) & (B, D) for all matrices D having the same number of rows as C, and let B1 be the largest cyclic component of B. Then A is equivalent to a matrix of the form (17). ed Ma nu " # B1 , B2 Proof. Let B1 be the largest cyclic component of B, and let B ? . As in the proof of Lemma 8, choose a row vector 0, B4 T dT1 such that (B Since (A,C) & (B, D) for all D, by the transitivity property and property R2, it follows ? ?1 , d1?) is ?observable. 0 ? ? ?? that (A,C) & ?B1 , ?dT1 ??. By Lemma 6, it follows that (A,C) & (?k B1 , ?k dT1 ), and by Lemma 5 it follows that (A,C) 0 contains a subsystem equivalent to (?k B1 , ?k dT1 ). Structural stability Any PMIG satisfying (11) (for some X22 ), that is subject to perturbations in some or all of its parameters will continue to produce unbiased input estimates provided (11) is still solvable. Such a PMIG is said to be structurally stable with respect to the perturbed parameters. Structurally stable parameters do not need to be known exactly for unbiased input estimates. The degree to which the solvability of (11) is stable with respect to the parameters of PIG and PMIG depends partly on the knowledge of the problem and, for multiple input systems, on the number of identical copies of the PIG within the PMIG. By way of example, consider a plant having two inputs. If the inputs are sinusoids of the same frequency but different amplitude and phase, then for unbiased input estimates the PMIG must have two undamped modes at this frequency. This PMIG does not, however, require a priori specification of the amplitude and phase of each mode. Ac 5.1 ce pt The amount of redundancy required in the model depends entirely on the knowledge of the plant; if the parameters of the plant dynamics and outputs are known exactly, a single copy of the exact model will suffice. DS-13-1302 16 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME For more complicated multiple-output signal generators, the canonical form of the system needs to be taken into consideration. As an example, the system ? ? # 0?0 " ? abc ? ?? ? ?? 0 0 0 ? , de f 000 ?? is equivalent to the system ? ? # 010 " ? y01 ? ?? ? ?? 0 0 0 ? , 100 000 ed ite d ?? 1 0 0 0 0 0 0 0 ? ? 0 " # ? 0? ? 1000 ? ?, ?, 1? 0 0 1 0 ? 0 ot 0 ?? 0 ?? ?? ?? 0 0 Co ?? py for some y depending on the parameters ?, a, b, c, d, e, f , and hence a model of this form is not robust with respect to all of these parameters. The model ? ? ? 0 x?1 ?x? ? ? 0 ? 2? ? ? ? = ? ?k1 ?k2 ?x?1 ? ? m 0 0 ?? ? ? ? 0 0 " # x1 ? ?x ? ? 0 0 ? F ? ? 2? ? ? 1 ? c2 ? ? ? + ? 1 ? ? ? ? ? F2 0 x? 1 m1 m1 ?c2 ?c3 1 x?2 0 m2 m2 ? ? x1 ?x ? ? 2? y = ? ?. ?x?1 ? x?2 1 0 0 1 Ma nu k2 ?c2 m1 m1 ?k2 ?k3 c2 m2 m2 (18) ed x?2 1 k2 m2 sc rip tN which is observable and contains two copies of the maximal cyclic component, is structurally stable with respect to the parameters ?, a, b, c, d, e, f . The examples that follow illustrate input estimation and structural stability when the Plant Model and PMIG are exact. These examples use the forced mass-spring-damper system shown in Figure 3. The system parameters are shown in Table 1. The dynamics of the system are described by Ac ce pt Further we assume that the measurement noise and process noise are both zero-mean Gaussians with covariance R and Q respectively. After substituing the parameters in Table 1 the plant is represented by A1 , B1 ,C1 matrices shown below. ? ? 0 0 1 0 ? 0 0 0 1? ? ? A1 = ? ? ??10 5 ?1.5 1.5? 5 ?10 1.5 ?3 ? ? 0 0 ?0 0? ? ? B1 = ? ? ?10 0 ? 0 10 C1 = I4О4 DS-13-1302 17 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Value k1 0.5 N/m k2 0.5 N/m k3 0.5 N/m m1 0.1 kg m2 0.1 kg c2 0.15 N.s/m c3 0.15 N.s/m py Parameter ed ite d Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME sc rip tN ot Co Table 1: Parameters of the double mass-spring-damper system. Ma nu Fig. 3: The double mass-spring-damper system. h i A2 = 0 " # 1 C2 = 0.5 ce pt ed Example 1. Step Input Consider the double mass-spring-damper system. Suppose the driving forces are unknown steps, F2 but the ratio between them is known. Suppose = 0.5, then PIG is represented using A2 and C2 matrices. F1 Ac Note that since x?2 = A2 x2 + w2 and A2 = [0], x2 (0) must be set to a non-zero value in order to generate a step input. In this example x2 (0) = 1. Assuming exact Plant Model specifying the PMIG so that it has a single eigenvalue at the zero consistent with the structure of a step input gives. A?1 = A1 B?1 = B1 C?1 = C1 A3 = A2 C3 = C2 . In PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG. Additionally we set the Q = I5О5 О 10?3 and R = 5 О I4О4 О 10?3 . Simulating the system response produces the results shown in Figure 4, with the key point to note being that the inputs converge to their correct values when the system transients have decayed. DS-13-1302 18 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME 3 4 2.5 3 4 3 2 2 1.5 2 1 1 0 0.5 1 -1 0 0 3 6 9 12 15 -0.5 0 3 6 9 12 15 time (s) time (s) 1 3 6 9 12 0 15 0 3 6 time (s) 0 0 12 15 0 3 6 Ma nu 9 time (s) 15 0 0 3 6 9 12 time (s) ot 0.5 6 12 tN 0.5 sc rip 1 3 9 time (s) 1 0 12 9 12 15 time (s) ed Fig. 4: Response of PIG and PMIG for step disturbance in double mass-spring-damper system. Ac ce pt Example 2. Compound Sinusoid-step Input Consider the double mass-spring-damper system. Suppose the driving forces in Figure 3 are a combination of two sinusoids function with known frequencies ?1 , ?2 and a step function. The magnitude and phase of these sinusoids and the magnitude of the the step function are unknown to the model, but the ratio between components in F1 and F2 are known. More specifically if the inputs are written as F1 = a1 u(t) + b1 sin(?1t + ?1 ) + c1 sin(?2t + ?2 ) F2 = a2 u(t) + b2 sin(?1t + ?1 ) + c2 sin(?2t + ?2 ) then a2 b2 c2 , , are known. a1 b1 c1 DS-13-1302 15 0.5 Co 0 9 py 0.5 -1 -2 6 1 1 0 3 1.5 3 2 0 time (s) 1.5 4 -2 ed ite d 0 19 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use 15 Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by example ASME suppose that PIG is represented by A and C matrices shown below. In this 2 2 ed ite d ? ? 0 0 0 0 0 ? ? ?0 0 ?1 0 0 ? ? ? ? A2 = ? ?0 ??1 0 0 0 ? ? ? ?0 0 0 0 ?2 ? 0 0 0 ??2 0 " # 1 0.5 0.5 1 ?1 C2 = 0.5 0.5 ?0.5 1 1 A?1 = A1 B?1 = B1 py h iT The initial state x2 (0) determines the phase and magnitude of the sinusoids in each input. We set x2 (0) = 0.5 0.25 0.25 0.5 0.5 . Additionally we choose f1 = 1Hz, f2 = 0.2Hz and set ?1 = 2? f1 , ?2 = 2? f2 . Assuming exact Plant Model and PMIG gives C?1 = C1 A3 = A2 C3 = C2 . tN ot Co In the PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG. Additionally we set the Q and R to the values of the previous example. The the system shown in Figure 2 is simulated and the responses are presented in Figure 5. The estimated inputs are seen to track the inputs after initial transients have decayed. Imperfect plant model: error bounds The results of the previous section showed that the asymptotic un-biasedness of the filter is not generally preserved under perturbations in parameters relating to the eigenvalues of the plant and signal. In practical applications, the parameters are never known exactly and the model is likely to capture only an approximation of the dynamics of interest, hence model error must be taken into consideration. This section shows that (i) the estimated bias is approximately proportional to the parameter error and (ii) the relative bias is bounded. The following result gives an estimate of the bias under the assumption that the model is ?close? to an exact model. Ma nu sc rip 6 0 , D0 ) is in error, such that Theorem 6. Suppose that the system model (A0L ,CM M A0L = AL + ?AL , D0M = DM + ?DM , pt ed 0 CM = CM + ?CM , Ac ce where (AL ,CM , DM ) is a model producing unbiased estimates of the system (AP , BP ,CP , DP ) governed by (3). Then the bias E[z0 ] of the imperfect system model satisfies E[z0 ] ? ?(?DM X + (DM + ?DM )?X)E[x] (19) where X satisfies (6) and ?X is the unique solution to ?XAP ? AL ?X = ?AL X. (20) Proof. The imperfect estimate satisfies x?? 0 = (AL + ?AL )x?0 + BP x DS-13-1302 20 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME 4 4 3 3 2 2 2 1 1 0 0 0 -1 -1 -2 -2 4 -2 0 3 6 9 12 15 -3 0 3 6 9 12 15 time (s) time (s) 2 12 -2 -2 15 0 3 6 1 1 0.5 0.5 0 0 6 9 12 time (s) -0.5 15 0 3 Ma nu 3 15 0 3 6 9 12 15 time (s) tN 1.5 12 sc rip 1.5 0 9 time (s) time (s) -0.5 15 py Co 9 ot 6 12 -1 -1 -4 3 9 0 0 0 6 1 -2 -6 3 2 1 0 0 time (s) 4 2 -6 ed ite d -3 -4 6 9 12 15 time (s) Fig. 5: Response of PIG and PMIG for combined sinusoids disturbance in double mass-spring-damper system. z0 = DP x ? (DM + ?DM )x?0 . pt ed and hence z0 satisfies Ac ce For the perfect model, apply Lemma 1 to obtain DM X = DP , where X satisfies XAP ? AL X = BP . Let ?X be the solution of (X + ?X)AP ? (AL + ?AL )(X + ?X) = BP , DS-13-1302 21 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME which simplifies to (20). Defining ? 0 = x?0 ? (X + ?X)x, we have ???0 = (AL + ?AL )??0 ? (X + ?X)w. The resulting bias E[z0 ] then satisfies E[z0 ] = (DP ? (DM + ?DM )(X + ?X))E[x] ? (DM + ?DM )E[??0 ] ed ite d = ?(?DM X + (DM + ?DM )?X)E[x] ? (DM + ?DM )E[??0 ]. Since ? 0 is stable, E[??0 ] = 0 and the result follows. Co py When the imperfect model contains multiple copies of an approximation of the maximal cyclic component of the system, the model error can be taken to be entirely due to the error in the eigenstructure, i.e. ?CM = 0, ?DM = 0, ?AL = ?AM , where ?AM is the error of the model relative to a perfect model which contains suitably replicated copies of the true maximal cyclic component of the system. The above result, although it requires the solution of the matrix equation (20), does give a measure of the error of the input estimate, modulo noise, relative to the plant state: Corollary 6.1. The bias relative to the plant state satisfies tN ot kE[u ? u?]k kE[z]k = ? k?DM k и kXk + (kDM k + k?DM k) k?Xk . kE[x]k kE[x]k sc rip Unfortunately the result does not give an error estimate relative to the input signal; for example if the input signal is sinusoidal then the relative error can be infinite at the instants when the input is zero. The next result provides an estimate of the bias when the parameter errors are large. Ma nu Lemma 10. For the plant and model governed by (3) with closed loop stability, the bias E[z] is O(t n eat ) for some a ? 0 and n ? 0. Moreover, if the plant state is bounded, with kE[x]k ? M, then E[z] is bounded with kE[z]k ? M kDP k +C1 M kDM k max(1, kBP k /х) for some positive constants C1 , х depending only on AL . ed Proof. The solution to (3) is pt x(t) = ?P (t)x(0) + 0 ?P (t ? ?)w(?)d?, (21) = exp(tAP ), x?(t) = ?L (t)x?(0) + ?L (t) (22) Rt 0 ?L (t ? ?)BP x(?)d?, (23) = exp(tAL ). (24) Ac ce ?P (t) Rt The matrix AL is stable, and thus k?L (t)k ? C1 e?хt for some positive constants C1 , х. This gives us the bound kE[x?(t)]k ? C1 e?хt kx(0)k +C1 kBP k DS-13-1302 22 Z t e?х(t??) kE[x(?)]k d?. 0 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME It easily follows that if kE[x(t)]k is O(t n eat ) as t ?? ? for some a > 0 then kE[x?(t)]k is also O(t n eat ). The bias E[z] is a linear combination of E[x?] and E[x] and thus of the same order. For the bounded case kE[x(t)]k ? M, we have kE[x?(t)]k ? C1 Me?хt +C1 M kBP k (1 ? e?хt )/х = C1 M kBP k /х + e?хt C1 M(1 ? kBP k /х) ? max(C1 M,C1 M kBP k /х), ed ite d and applying the triangle inequality to z = DP x + DM x?, py the result follows. sc rip tN ot Co The sensitivities can be difficult to calculate symbolically, so we give a numerical example here. Ma nu Fig. 6: The block diagram of the mass spring damper system. Example 3. Consider the mass spring damper system in Figure 6 with step-plus-sine and parameters w = 1, m = 1, k = 1, c = 10. Suppose that model parameters are in error, with k? = k + ?k c? = c + ?c, ed w? = w + ?w m? = m ce pt where ?w, ?k, ?c are small. For the stabilizing gains " Ac Kb = ? ? 0 ? ? Kc = ?25? , 50 # 0 100 the tracking error satisfies ku ? u?k ? (8.77 |?w| + 0.21 |?c| + 1.41 |?k|) kxk . Example 4. Step-Sinusoid Input with Imperfect Model Consider the double mass-spring-damper system shown in Figure 3. Suppose that the Plant Model and PMIG contain errors as shown in Table 2. DS-13-1302 23 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME Parameter Plant Parameter Plant Model 0.5 N/m k?1 0.55N/m k2 0.5 N/m k?2 0.55N/m k3 0.5 N/m k?3 0.55N/m m1 0.1 kg m?1 0.11kg m2 0.1 kg m?2 0.11kg c2 0.15 N.s/m c?2 0.165 N.s/m c3 0.15 N.s/m c?3 0.165 N.s/m Parameter PIG Parameter PMIG f 1Hz f? 1.1Hz py ed ite d k1 Co Table 2: Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3. B?1 6= B1 A3 6= A2 tN A?1 6= A1 ot Therefore C?1 = C1 . sc rip The initial state of the PIG, x2 (0), determines the real input. Specifically x2 (0) determines the phase of sinusoid component, h iT the magnitude of the sinusoid component, and the size of step input component. In this example x2 (0) = 0.25 0.25 1 . By substituting the model parameters into Equation 18 the following matrices are calculated for Plant Model. ? 0 0 1 0 ? ? 0 0 0 1 ? ? A?1 = ? ? ??10 5 ?1.5273 1.5273 ? 5 ?10 1.5273 ?3.0545 ? ? 0 0 ? 0 ? 0 ? ? B?1 = ? ?. ?9.0909 ? 0 0 ?9.0909 pt ed Ma nu ? ce C?1 = C1 = I4О4 Ac Furthermore PMIG is given by ? ? ? ? o 2? f? 0 0 6.9115 0 ? ? ? ? A3 = ??2? f? 0 0? = ??6.9115 0 0? 0 0 0 0 0 0 " # 0.5 0.5 1 C3 = C2 = . 0.5 ?0.5 0.5 Using the above matrices yields the results shown in Figure 7. To verify that the input estimation bias satisfies CorolDS-13-1302 24 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) lary2017 6.1, by ASME kE[uu ? u?u]k kE[xx]k is plotted at every time step. The supremum of steady state estimate bias kE[uu ? u?u]k kuu ? u?uk = lim } t?? kx? xk kE[xx]k ed ite d sup{ is compared with k?DM k и kXk + (kDM k + k?DM k) k?Xk kE[[y1 , y2 , y?1 , y?2 , x2 ]]k ot and calculate T F 1 , F? F 2 ]T E[[[F 1 , F 2 ] ? [ F? = [y?1 , y?2 , y??1 , y??2 , x 3 ] Co kE[uu ? u?u]k = kE[xx]k T F 1 , F? F 2 ]T E[[[F 1 , F 2 ] ? [ F? py . In this example we use sc rip tN T F 1 , F? F 2 ]T [[F 1 , F 2 ] ? [ F? [y?1 , y?2 , y??1 , y??2 , x 3 ] for every point obtained from the simulation. Figure 8 shows that as the response converges to steady state, revealed. The figure shows kE[uu?u?u]k kE[xx]k is kE[uu ? u?u]k } = 0.1408. kE[xx]k Ma nu sup{ Additionally by substituting the matrices into Equation 4 and solving Equation 20 the following can be obtained. Ac ce pt ed ?DM = 0 2О7 " # 00001 1 1 DM = 0 0 0 0 1 ?1 1 ? ? ?0.1288 0.1349 ?0.0198 0.0376 0.2144 0.0855 0.0348 ? 0.0844 ?0.0888 0.0236 ?0.0383 ?0.0416 ?0.2628 ?0.0223? ? ? ? ? ? 0.2623 ?0.3413 ?0.0571 0.0207 0.1697 1.2240 0.0179 ? ? ? ? ?X = ? ??0.4445 0.5172 ?0.0104 0.0645 1.6985 0.0463 0.0513 ? . ? ? ??0.0415 0.0426 ?0.0088 0.0142 0.1291 1.2262 0.0124 ? ? ? ? 0.0384 ?0.0472 ?0.0049 ?0.0002 ?1.2068 0.0869 ?0.0011? 0.0167 ?0.0175 ?0.0017 ?0.0003 ?0.0522 0.0587 0.0863 Therefore k?DM k = 0, kDM k = 1.3229, kXk = 1, k?Xk = 2.1996 k?DM k и kXk + (kDM k + k?DM k) k?Xk = kDM k k?Xk = 2.9098 DS-13-1302 25 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 bythat ASME It follows kE[uu ? u?u]k ? k?DM k и kXk + (kDM k + k?DM k) k?Xk . kE[xx]k 3 3 2 2 1 1 4 3 ed ite d 2 1 0 0 0 -1 -1 -1 -2 3 6 9 12 15 0 3 6 time (s) 9 12 time (s) 1 1 0.8 0 0.6 -1 0.4 -2 0.2 -3 0 3 6 9 12 15 0 1.5 1 0.5 3 Ma nu 0.5 0 6 9 12 15 12 15 time (s) 1 6 9 12 0.6 0.4 0.2 15 time (s) time (s) 1 3 0.8 0 0 3 6 9 time (s) sc rip 0 0 1.2 ot 1.2 2 tN 3 -3 Co 4 15 py 0 0 -0.5 -0.5 0 3 6 9 12 0 3 6 9 12 15 time (s) ed time (s) 15 Conclusion We have shown that for unbiased input or disturbance estimation using the Kalman Inverse Filter, an exact model of the plant and a suitably replicated model describing the structure of the input signals is required. Moreover, we have shown that the model of the input signal generator is structurally stable with respect to all parameters that do not change the eigen structure of the PMIG. A new proof of the internal model principle has also been given. These results provide insight to the design of Kalman Inverse Filters. In several of the recent papers using this method, e.g. [6?13] the structural form of the PMIG is a chain of integrators without clear justification for these choices. The results of this paper show that this approach amounts to approximating the input signal as a piecewise polynomial whose order corresponds to the degree of the minimal polynomial of the PMIG matrix. A single integrator corresponds to a system with DS-13-1302 Ac 7 ce pt Fig. 7: Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper system with imperfect model parameters. 26 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME Bias relative to state 10 bias relative to state upper bound 9 8 7 ed ite d 6 5 4 3 py 2 Co 1 0 4 6 8 time (s) 14 bias relative to state versus upper bound calculated by k?DM k и kXk + (kDM k + k?DM k) k?Xk. tN kE[uu?u?u]k kE[xx]k , 12 sc rip Fig. 8: 10 ot 2 Ma nu dynamics x? = 0 and thus would be consistent with an input having the structure of a constant (i.e. step) input. Similarly a chain of two integrators is consistent with the input having the structure of a ramp. For example, in [12], a double integrator is used to model the signal associated with each of the unknown forces acting on the wheels of a vehicle. Similarly a chain of three integrators yield a piecewise parabolic fit to an input. De Nicolao [16], in his work on Tikhonov regularization has similarly observed that chains of integrators can be associated with smoothing splines. ed For applications where the inputs have a higher structure that is known, e.g. harmonic inputs and/or disturbances, the insights provided by this work informs us that the the PMIG should be modelled as a resonator having the appropriate frequency. Where the signal is not a pure sinusoid, the PMIG might be extended to include harmonics of the base frequency. The efficacy of this approach has been demonstrated in [14]. If the system has multiple inputs suitable replication of the input structure is required. We note finally the possibility that where the structure of the input may be time varying, or not precisely defined, the PMIG could include appropriately rich structure to adapt to the input signal. 8 Ac ce pt From the point of view of the philosophical foundations of the design and in terms of guiding practical implementation of input-and-disturbance estimators, we have found the principles revealed in this paper extremely powerful in applications where we seek to estimate inputs and disturbances in systems using the KIF. Acknowledgements Kianoosh Soltani Naveh would like to acknowledge The University of Queensland Centennial Scholarship for the support of this work. References [1] Bayless, J. W., and Brigham, E. O., 1970. ?Application of the Kalman filter to continuous signal restoration?. Geophysics, 35(1), Feb, pp. 2?23. DS-13-1302 27 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Ac ce pt ed Ma nu sc rip tN ot Co py ed ite d Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c)[2]2017 by ASME Crump, N. D., 1974. ?A Kalman filter approach to the deconvolution of seismic signals?. Geophysics, 39(1), Feb, pp. 1?13. [3] Dimri, V., 1992. Deconvolution and Inverse Theory. Elsevier Publishing Company. [4] Leite, L. W. B., and da. Rocha, M. P. C., 2000. ?Deconvolution of non-stationary seismic process?. Rev. Bras. Geof., 18(1), pp. 75?89. [5] Rocha, M. P., and Leite, L. W., 2003. ?Treatment of geophysical data as a non-stationary process.?. Mat. Apl. Comput., 22(2), pp. 149?166. [6] Stalford, H. L., 1981. ?High-alpha aerodynamic model identification of t-2c aircraft using the ebm method?. Journal of Aircraft, 18(10), pp. 801?809. [7] Sri-Jayantha, M., and Stengel, R. F., 1988. ?Determination of nonlinear aerodynamic coefficients using the estimationbefore-modeling method?. Journal of Aircraft, September, pp. 796?804. [8] Ray, L. R., Ramasubramanian, A., and Townsend, J., 2001. ?Adaptive friction compensation using extended KalmanBucy filter friction estimation?. Control Engineering Practice, 9(2), pp. 169?179. [9] Ray, L. R., Townsend, J., and Ramasubramanian, A., 2001. ?Optimal filtering and Bayesian detection for friction-based diagnostics in machines?. ISA Transactions, 40(3), pp. 207?221. [10] Austin, K. J., and McAree, P. R., 2007. ?Transmission friction in an electric mining shovel?. Journal of Field Robotics, 24(10), pp. 863?875. [11] Ray, L. R., 1995. ?Nonlinear state and tire force estimation for advanced vehicle control?. IEEE Transactions on Control Systems Technology, 3(1), pp. 117?124. [12] Ray, L. R., 1997. ?Nonlinear tire force estimation and road friction identification: Simulation and experiments?. Automatica, 33(10), pp. 1819?1833. [13] Siegrist, P., and McAree, P. R., 2006. ?Tyre-force estimation by kalman inverse filtering: applications to off-highway mining trucks?. Vehicle system dynamics, 44(12), pp. 921?937. [14] Reid, A. W., McAree, P. R., Meehan, P., and Gurgenci, H., 2014. ?Longwall shearer cutting force estimation?. Journal of Dynamic Systems, Measurement, and Control, 136(3). [15] Tikhonov, A. N., and Arsenin, V. Y., 1977. Solutions of ill-posed problems. Washington : Winston. [16] G., D. N., and G., F.-T., 2001. ?Regularization networks: Fast weight calculation via kalman filtering?. IEEE Transactions On Neural Networks, 12(2), March, pp. 228?235. [17] Nicolao, G. D., and Ferrari-Trecate, G., 2003. ?Regularization networks for inverse problems: A state-space approach?. Automatica, 39(4), April, pp. 669?676. [18] Francis, B. A., and Wonham, W., 1976. ?The internal model principle of control theory?. Automatica, 12(5), pp. 457? 465. [19] Gelb, A., 1974. Applied optimal estimation. MIT press. [20] Gantmacher, F. R., 1959. The Theory of Matrices, Vol. 1. Chelsea Publishing Company. [21] Horn, R. A., and Johnson, C. R., 1985. Matrix Analysis. Cambridge University Press. [22] Datta, B. N., 1994. ?Linear and numerical linear algebra in control theory: Some research problems?. Linear Algebra and Its Applications, 197, pp. 755?790. [23] Casti, J. L., 1987. Linear Dynamical Systems. Academic Press. DS-13-1302 28 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use ned sinusoids disturbance in double mass-spring-damper system. . . The block diagram of the mass spring damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper system with imperfect model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bias and bias upper bound for imperfect plant model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ed ite d 1 2 3 4 18 19 21 23 26 27 Co py List of Tables 1 Parameters of the double mass-spring-damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2 Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3. 24 Introduction Inverse problems in systems theory involve reconstructing the inputs to a system from measured responses. The adaptation of the Kalman filter to inverse problems, sometimes termed the Kalman inverse filter (KIF), was first proposed by Bayless and Brigham [1] as a framework for removing sensor reverberation from seismic signals. The idea is to augment the states of the plant model by additional dynamics that captures known structure in the inputs to be estimated. This is equivalent to assuming the autocorrelation function of the unknown inputs or distrubances that are to be estimated. The formulation of the KIF in [1] was for continuous-time systems. A discrete-time formulation followed shortly after in [2]. Most of the applications of the KIF over the last thirty years have been for geophysical signal processing, see for example [3?5], but the method has also attracted interest for input and disturbance identification in mechanical systems. Applications include estimation of aerodynamic loads on aircraft [6, 7], friction in machines [8?10], and the estimation of machine-environment interaction forces [11?14]. The KIF is closely linked to other methods for input identification, e.g. the regularization method of Tikhonov [15]. Work by De Nicolao and Ferrari-Trecate [16, 17] has shown that under appropriate assumptions, Tikhonov regularization coincides with Bayesian estimation where the autocorrelation function of the unknown input is specified and the plant is a linear system driven by white noise. De Nicolao and Ferrari-Trecate [16] show that the inverse problem can be solved by application of the Kalman filter, revealing an equivalence between Tikhonov regularization and the KIF. As a method for input identification the KIF has a number of attractions including: (i) that it is recursive and therefore suitable for real-time implementations; (ii) the underlying Kalman filter structure accommodates model and measurement uncertainty; and (iii) the method can be adapted to non-linear systems through the extended Kalman filter. Limitations of the method are: (i) that the structure of the inputs/disturbances must be specified, and (ii) the bandwidth to which input/disturbances can be reconstructed depends on the dynamics of the plant, and the process and measurement uncertainties, as specified in the Kalman filter. The structure of the KIF is illustrated in Fig. 1, and has four components: the Plant, the Plant Input Generator (PIG), the Plant Model and the Plant Model Input Generator (PMIG). The Plant is the system whose inputs are to be estimated. The PIG represents the process that generates these to-be-estimated inputs/disturbances acting on plant. The KIF comprises the Plant Model and the PMIG. The PMIG assumes the inputs have structure. The idea is to design the PMIG so that its outputs serve as estimates of the input to the Plant. This is done through the feedback structure of the Kalman filter which provides an implicit inverse model of the Plant. This paper establishes the conditions needed for convergence of the PMIG output to the Plant input in steady state for linear plants. Specifically, the paper establishes the conditions that must be satisfied by the Plant Model and the PMIG for convergence of the plant input/disturbance estimates to their actual values. Section 2 describes the structure of the Kalman Ac ce pt ed Ma nu sc rip tN ot 1 DS-13-1302 2 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME real system Plant Input Generator (PIG) Plant plant inputs (unknown) + estimation error ed ite d estimated plant inputs system model/estimator - Plant Model Input Generator (PMIG) + py Plant Model Co Kalman gain ot Fig. 1: Kalman Inverse Filter block diagram. The Kalman Inverse Filter The KIF is described by the block diagram in Fig. 2 and governed by the equations Ac ce pt ed 2 Ma nu sc rip tN Inverse Filter and details the notation used throughout the paper. Section 3 summarizes key assumptions. The main work of the paper begins in Section 4 where we give general conditions for an asymptotically unbiased KIF by formulating the problem in terms of the solution to a constrained linear matrix equation. It is shown that if the model of the Plant and PMIG are exact, the estimates produced by the KIF converge to their true values. Section 5 presents a more detailed investigation of the requirements on the PMIG and shows that for unbiased estimation it must span all possible inputs to the plant and that the robustness of the estimator with respect to errors in parameters depends on the eigen-structure of the PIG. It is shown, moreover, that a robust PMIG can be constructed from knowledge of just the number of outputs and the minimal polynomial of the PIG dynamics. Inter alia, the analysis also yields a simple proof of the Internal Model Principle which was first proven in [18]. Section 6 gives error bounds on estimates given knowledge of model errors. The results of this section show that the estimated bias is approximately proportional to the parameter error and that relative bias is bounded. DS-13-1302 x?1 = A1 x1 + B1 u + w1 (1a) x?2 = A2 x2 + w2 (1b) y = C1 x1 + v (1c) u = C2 x2 x?? 1 = A?1 x?1 + B?1 u? + Kb e (1d) x?3 = A3 x3 + Kc e (1f) (1e) y? = C?1 x?1 (1g) u? = C3 x3 (1h) e = y ? y? (1i) z = u ? u? (1j) 3 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME where thebyvector x is the state of the Plant, x is the state of the Plant Input Generator (PIG), u is the unknown plant input 1 2 produced by the PIG and y is the measured output signal from the plant. The vector x?1 is the state of the plant model, u? is the estimate of the input which is produced by the Plant Model Input Generator (PMIG) with state vector x3 generally not the same size as x2 , and e is the innovation sequence available for feedback. The matrices A1 , A2 , B1 , C1 , C2 , A?1 , B?1 , C?1 , A3 , C3 have appropriate sizes satisfying the Eqns. (1). In particular the PMIG and the PIG have the same number of outputs and the Plant and the Plant Model have the same number of outputs. Plant PMIG Plant Model Ma nu sc rip tN ot Co py ed ite d PIG ed Fig. 2: The Kalman Inverse Filter, including the open-loop system governed by A1 and A2 , and the closed-loop feedback system governed by A?1 and A3 . Ac ce pt We note, however, that the Plant Model and PMIG may contain additional dynamics not present in the Plant and PIG under investigation or vice-versa. Hence the model state vectors x?1 , x3 may be of different sizes to the corresponding plant w (t) state vectors x1 , x2 . The measurement noise v and the process noise w(t) = w1 (t) are vector white noises with covariance 2 structure E[w(s)wT (t)] = Q?(s ? t) E[v(s)vT (t)] = R?(s ? t) E[v(s)wT (t)] = 0. The system can be thought of as regulating z, the difference between the input signal and the estimated input signal. When E[z] = 0 the estimates u? track the actual inputs to the Plant u. DS-13-1302 4 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME The by state vector of the combined system encompassing the Plant and the PIG is " # x1 x= . x2 " # x?1 x? = , x3 which shall form the model state. The state equations (1) combine to form the system y = CP x (3b) (3c) (3d) Co y? = CM x? (3a) py x? = AP x + w, x?? = AL x? + BP x + Kv, ed ite d The state vector of the model of the system, encompassing the Plant Model and the PMIG, is (3e) ot z = DP x ? DM x?, tN where " # A?1 B?1C3 AM = 0 A3 sc rip " # A1 B1C2 AP = , 0 A2 AL = AM ? KCM , h i CP = C1 0 , h i DP = 0 C2 , " # Kb K= , Kc BP = KCP h i CM = C?1 0 , h i DM = 0 C3 Ma nu (4) T ?1 K = PCM R pt ed The so-called Kalman gain K is given by T ?1 R CM P = 0. AM P + PATM + Q ? PCM Ac ce where P satisfies the Algebraic Riccati Equation [19] In order to regulate z, it is necessary to have closed loop stability, that is ?(AL ) ? C? , where ?(AL ) is the eigenvalue spectrum of the matrix AL and C? is the open left half of the complex plane. A necessary condition for existence of a stabilizing gain K is that the pair (AM ,CM ) is detectable. We assume that the pair (AM ,CM ) is observable; this guarantees existence of a stabilizing gain K such that AP and AL have no common eigenvalues. DS-13-1302 5 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) by ASME 3 2017 Assumptions The following assumptions are made at various places and referred to by number as required. They are gathered here to avoid repetition. As indicated, some of the assumptions may be replaced with less restrictive alternatives. ed ite d A1 The plant system is observable, i.e. (AP ,CP ) is observable. A2a All eigenvalues of the plant system matrix AP are unstable. A2b The plant system matrix AP and the closed loop system matrix AL have no common eigenvalues. A3 The system model is observable, i.e. (AM ,CM ) is observable. A4 All eigenvalues of the closed loop matrix AL are stable. A5a The Plant Model exactly describes the Plant, i.e. (A?1 , B?1 , C?1 ) = (A1 , B1 ,C1 ). A5b The Plant Model has the same transfer function as the Plant. A6 The Plant matrix A1 has no eigenvalues in common with the Plant Input Generator matrix A2 . py Assumption A4 is essential to the results of this paper. Assumption A3 implies the existence of stabilizing gains K such that Assumption A4 is satisfied. Convergence properties of the KIF This section gives general conditions for an asymptotically unbiased KIF, namely E[z(?)] = 0. The approach taken is to reformulate the problem in terms of the solution to a constrained linear matrix equation by decomposing the trajectories into terms arising from the state and the state estimation error. Then using Assumption A4, the conditions for E[z(?)] = 0 are derived. tN General conditions for unbiased input estimation Consider the state equations (3), then x and x? are related by sc rip 4.1 ot Co 4 x? = ? + Xx (5) Ma nu where X has appropriate dimensions and ? represents the state estimation error. Differentiating both sides with respect to time and using state equations (3) we can write ??? = x?? ? X x? = AL x? + BP x + Kv ? X(AP x + w) = AL (?? + Xx) + BP x + Kv ? X(AP x + w) ed = AL ? + (AL X + BP ? XAP )x + Kv ? Xw. XAP ? AL X = BP (6) Ac ce pt Under Assumption A4, AL is stable and hence AL ? converges to zero as t ? ?. Therefore, in order to satisfy limt?? ? = 0, the term due to x must be zero, so must hold. This relation is the well-known Sylvester matrix equation and has unique solution X under Assumption A2b that is if AP and AL have no common eigenvalues [20?22]. Assumption A2b can be satisfied in a number of ways. One sufficient condition is that Assumptions A2a and A4 hold. Using the above relations, the input estimation error z can then be expressed in terms of the state and state estimation error. Specifically z = DP x ? DM x? = (DP ? DM X)x ? DM ? . DS-13-1302 6 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME Since E[??by ]? 0 the term to due to ? converges to 0 . Therefore in order for E[z] ? 0 , the component due to x must be zero, leading to the following results: Lemma 1 (Sufficient conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if DM X = DP . (7) where X is the solution to the matrix equation (6). ed ite d Lemma 2 (Necessary conditions for zero bias). Under Assumptions A2b and A4, the bias is asymptotically zero if and only if (DP ? DM X)U = 0 (8) where X is the solution to the matrix equation (6) and U is a matrix whose columns span the unstable subspace of AP . Co py Note that when all eigenvalues of AP are stable, U = 0 and (8) holds regardless of the accuracy of the model. The KIF is asymptotically unbiased because both the input signal and the input estimate converge to zero. Conversely, if all eigenvalues of AP are unstable (Assumption A2a), U = I in (8) which reduces to to (7). Exact Plant Model, Exact PMIG Before discussing the question of robustness with respect to parameters, we remark that for an exact model the bias is zero, i.e. for tN ot 4.2 sc rip (A?1 , A3 , B?1 , C?1 ,C3 ) = (A1 , A2 , B1 ,C1 ,C2 ), A?1 = T1 A1 T1?1 , A3 = T2 A2 T2?1 , Ma nu the equations (6) and (7) are satisfied by X = I. In fact, the choice of coordinates for the model state space does not affect the bias: given invertible matrices T1 and T2 , the model described by C?1 = C1 T1?1 , B?1 = T1 B1 , C3 = T1C2 T2?1 , " # T1 0 X= . 0 T2 pt ed has zero bias, with equations (6) and (7) satisfied by Ac ce Here the converged state estimates are related to the true states by x?1 = T1 x1 x3 = T2 x2 . 4.3 Zero-bias equations To show that for unbiased input estimation, the plant model must be similar to the plant, partition the matrix X into four blocks of compatible size, with " # X11 X12 X= X21 X22 DS-13-1302 7 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) by equation ASME (6) becomes The2017 matrix X11 A1 = (A?1 ? KbC?1 )X11 + B?1C3 X21 + KbC1 (9a) X21 A1 = ?KcC?1 X11 + A3 X21 + KcC1 (9b) X11 B1C2 + X12 A2 = (A?1 ? KbC?1 )X12 + B?1C3 X22 (9c) X21 B1C2 + X22 A2 = ?KcC?1 X12 + A3 X22 , (9d) ed ite d and the zero bias equation (7) becomes C3 X21 = 0 C3 X22 = C2 . (10a) (10b) Co py Any given Plant and PIG has a very large class of models for which (9) and (10) can be solved for X. We limit our attention to models for which the Plant Model has the same dimensionality as the true Plant. In this situation the Plant and the Plant Model must necessarily obey a similarity relationship. The argument is as follows. After substitution of (10a), (9a) can be rearranged as ot A?1 X11 ? X11 A1 = Kb (C?1 X11 ?C1 ). tN In order for this equation to hold regardless of the feedback gain Kb , the left and right sides must be zero, giving A?1 X11 = X11 A1 , sc rip C?1 X11 = C1 . Similarly, (9c) can be reduced to Ma nu (X11 B1 ? B?1 )C2 = (A?1 ? KbC?1 )X12 ? X12 A2 , which in order to hold regardless of the signal generator gain C2 necessitates ed X11 B1 ? B?1 = 0. pt Since the plant (A1 ,C1 ) is observable and the model has the same number of dimensions, X11 must be square, invertible and unique (see Sec. 5), leading to the following result: Ac ce Theorem 1. Suppose the Plant (A1 , B1 ,C1 ) is observable and controllable, and that the associated zero bias Plant Model (A?1 , B?1 , C?1 ) has the same number of state space dimensions. Then the Plant and the Plant Model must be related by a unique similarity transform ?1 A?1 = X11 A1 X11 , ?1 C?1 = C1 X11 , B?1 = X11 B1 . Moreover the Plant and the Plant Model have the same transfer function C?1 (s ? A?1 )?1 B?1 = C1 (s ? A1 )?1 B1 for all s ? C. DS-13-1302 8 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 ASME Thisby result implies a requirement for exact knowledge of the Plant transfer function: to exactly determine the plant inputs from the plant outputs requires the exact relationship between the inputs and the outputs. In general this requirement implies a need for a perfect Plant Model including exact parameter knowledge. An interesting exception is where all eigen-values of A1 are zero corresponding, for example, to free rigid-body dynamics problems. In this situation the eigen-structure is independent of the parameter values. The model used by Siegrist [13] to estimate road tyre interaction forces for off-highway mining vehicles is in this class of problems. Conditions for unbiased input estimation with a perfect plant model For algebraic simplicity, we will henceforth assume a perfect Plant Model with A?1 = A1 , B?1 = B1 , C?1 = C1 so that X11 = I. We now seek conditions for unbiased input estimation under this assumption. A sufficient condition for zero bias is the existence of a solution to the constrained Sylvester equation ed ite d 4.4 C3 X22 = C2 (11b) py A3 X22 = X22 A2 . (11a) Co The solvability of (11) becomes a necessary condition for zero bias when the signal generator has no modes in common with the plant, as demonstrated in the following result. tN ot Theorem 2. Under Assumptions A2b, A4, A5a, suppose that the PMIG matrix A3 has no eigenvalues in common with the Plant Model matrix A?1 and that Kb has been chosen so that (A?1 ? KbC?1 ) has no eigenvalues in common with A1 or A2 (this is possible because (A1 ,C1 ) is observable). The KIF has zero bias if and only if there exists a solution X22 to the simultaneous matrix equations (11). sc rip Proof. To prove sufficiency, let X22 be a solution to (11). Then the equations (9), (10) are uniquely satisfied by " # I 0 X= . 0 X22 Ma nu To prove necessity, suppose the KIF has zero bias, i.e. (10) holds. Then (9a) reduces to X11 A1 = (A?1 ? KbC?1 )X11 + KbC1 . X21 A1 = A3 X21 X12 A2 = (A?1 ? KbC?1 )X12 ce pt ed Since the plant model (A?1 , C?1 ) is exact and A1 has no eigenvalues in common with (A?1 ? KbC?1 ), the solution X11 = I is unique. The next two equations (9b) and (9c) reduce to Ac which, by the assumption of no common eigenvalues, are uniquely satisfied by X21 = 0 and X12 = 0. The other equation (9d) reduces to X22 A2 = A3 X22 as required. We now show that equations (11) imply a relationship between the observability matrices of the PMIG and the PIG. Suppose that equations (11) hold for the matrix X22 . Multiply (11a) on the right by A2 and substitute (11b) to obtain C3 A3 X22 = C2 A2 . DS-13-1302 9 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 bythe ASME Repeating process will yield a sequence of equations j j j ? 0. C3 A3 X22 = C2 A2 , Terminating the sequence at j = k ? 1,1 the system can be expressed as Ok (A3 ,C3 )X22 = Ok (A2 ,C2 ), (12) ed ite d where Ok (A,C) is the kth order observability matrix defined by ? C ? ? ? CA ? ? Ok (A,C) ? ? . ? ?. ? .. ? CAk?1 py ? Co When k = n, where A is n О n, the matrix O (A,C) is the observability matrix of control theory. The noise-free component of an input signal produced by state x2 of the signal generator satisfies ot u = C2 x2 , tN u? = C2 x?2 = C2 A2 x2 , u? = C2 x?2 = C2 A22 x2 , sc rip .. . du dx2 = C2 k?1 = C2 Ak?1 2 x2 . dt k?1 dt Ma nu Similarly the input estimate produced by state x3 satisfies u? = C3 x3 , u?? = C3 x?3 = C3 A3 x3 , pt ed u?е = C3 x?3 = C3 A23 x3 , .. . d u? dx3 = C3 k?1 = C3 Ak?1 3 x3 . dt k?1 dt ce These relations can be summarized in two systems of equations u u? u? .. . Ac ? ? ? ? ? ? ? ? ? 1 As du dt k?1 ? ? ? ? ? ? ? = O (A2 ,C2 )x2 , ? ? ? ? ? ? ? ? ? ? ? u? u?? u?е .. . d u? dt k?1 ? ? ? ? ? ? = O (A3 ,C3 )x3 . ? ? ? a consequence of the Cayley-Hamilton theorem of linear algebra, the sequence can be stopped at k = n2 n3 , where A2 is n2 О n2 and A3 is n3 О n3 ; in the next section we show that k = n3 + 1 suffices. DS-13-1302 10 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017isby ASMEand x = X x , the PIG produces an identical signal as shown by the following relation. If (12) satisfied 3 22 2 ? ? ? ? ? ? ? ? ? u u? u? .. . du dt k?1 ? ? ? ? ? ? ? ? ? ? ? = O (A3 ,C3 )X22 x2 = O (A3 ,C3 )x3 = ? ? ? ? ? ? ? u? u?? u?е .. . d u? dt k?1 ? ? ? ? ? ?. ? ? ? tN ot Co py ed ite d The interpretation we give to this result is that the signals generated by the PMIG must span all possible inputs to the plant, so that if, for example, the plant input contains a sinusoidal component of frequency ?, the PMIG must include a mode with this same frequency. Furthermore, the matrix X22 is a map from the states of the PIG to the states of the PMIG so that the inputs to the model correspond to those of the plant in a way that is consistent with the dynamics of the two systems. The above results provide a test of whether a given PMIG will produce unbiased input estimates for a given PIG. No assumptions have been made about the size of the PMIG relative to the PIG. The PMIG may have additional signal components that do not appear in the PIG, or multiple copies of components of the signal generator. The PMIG can therefore contain a library of modes whose form corresponding to all expected inputs to be applied to the plant. By way of example, if it is known that the plant can subjected to inputs/disturbances having n distinct frequencies, ?1 , . . ., ?n , applied alone or in combination, the PMIG should be designed to have eigenvalues corresponding to those frequencies. The next section looks in more detail at the relationship between the structure of the signal generator and the structure of models that produce zero-bias signal estimates. The analysis shows that the robustness of the estimator with respect to certain parameters depends largely on the eigen-structure of the PIG. This amounts to determining the general structure of pairs of matrices (A3 ,C3 ) and (A2 ,C2 ) which satisfy Equation (11) for some X22 . Conditions on the PMIG for zero bias estimates In this section, we consider the abstract properties of the relationship between the PIG and the PMIG implied by the zero bias equations (11), culminating in the result that a robust PMIG can be constructed from knowledge of just the number of outputs and the minimal polynomial of the PIG dynamics. For ease of argument, we introduce the following notation: sc rip 5 Ma nu Definition 1. For any matrices A,C, B, D of size n О n, k О n, m О m, k О m respectively, say (A,C) models (B, D), denoted (A,C) & (B, D), if (13a) AX = XB, (13b) ed CX = D pt for some matrix X of size n О m. Say (A,C) &X1 (B, D) if (11) holds for X = X1 . In particular, the PMIG will provide zero bias estimates for a given PIG if (A3 ,C3 ) & (A2 ,C2 ). ce The relation has useful ordering properties: Ac Reflexive: (A,C) &In (A,C), where In is the n О n identity matrix. This means any plant is a model of itself. Transitive: (A1 ,C1 ) &X1 (A2 ,C2 ) &X2 (A3 ,C3 ) implies (A1 ,C1 ) &X1 X2 (A3 ,C3 ). This means any model of a model of the plant is, in turn, a model of the plant. From these properties, two equivalence relations naturally arise: Similarity: Say (A,C) ? (B, D) if (A,C) &X (B, D) where the matrix X is square and invertible. This is similarity in the usual linear systems sense. Equivalence: Say (A,C) h (B, D) if (A,C) & (B, D) & (A,C). Similarity implies equivalence, but not necessarily vice-versa. Pairs of models can be related in three fundamental ways (subject, of course, to compatibility of the size of the sub blocks): DS-13-1302 11 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) by ASME R12017 (change of basis): " (P?1 AP,CP) &! (Q?1 AQ,CQ) by X"=#P?1 Q # h i A B1 I R2 (bigger model): , C D & (A,C) by X = 0 B2 0 " # ! h i h i A 0 R3 (unobservable signals): (A,C) & , C0 by X = I 0 B3 B4 ed ite d Relation R1 says that the choice of coordinate system is irrelevant, relation R2 says that the PMIG can contain additional dynamics not present in the PIG. Relation R3 says that unobservable components of the PIG can be ignored. In fact these relations characterize &, since any matrix X that relates (A,C) and (B, D) can be reduced by row and column operations to echelon form by P?1 XQ = 0I 00 for some invertible matrices P and Q2 . This is summarized in the following result: Theorem 3. (A,C) & (B, D) if and only if there are change of basis matrices P and Q that reduce (A,C) and (B, D) to upperand lower- block triangular form respectively such that (Q?1 BQ, DQ) = # ! i A12 h , C1 C2 A21 # ! h i 0 , C1 0 . B22 py (P AP,CP) = " A11 0 " A11 B21 Co ?1 Ma nu Lemma 3. If (A,C) &X (B, D), then sc rip tN ot The above result characterizes the perfect plant model as containing a subsystem that is equivalent to the observable part of the plant. In terms of structural stability, this says that the model is insensitive to perturbations that do not change the equivalence of this subsystem that models the plant. The aim of the remainder of this section is to show that any observable model (A,C) with suitable replication of the eigenstructure of B is a perfect model of (B, D), i.e. contains a subsystem equivalent to (B, D). The implication of this result is that the model is insensitive to any perturbations that do not alter the eigenstructure or compromise the observability of the model. In order to achieve this result, we re-examine the results relating the observability matrices of the systems. In the previous section, we demonstrated that (11) implies (12), and hence Ok (A,C) ? X = Ok (B, D) for all k ? 1. (14) pt ed The reverse is not necessarily true, however note that (14) implies CAi (AX ? XB) = 0 for i ? k ? 2, leading to the following result: On (A,C) ? (AX ? XB) = 0. (15) Ac ce Lemma 4. If (14) holds for k = n + 1, then (AX ? XB) satisfies Moreover if (A,C) is observable (i.e. On (A,C) has rank n), AX ? XB = 0 and hence (A,C) & (B, D). Several results immediately follow: Corollary 3.1. 2 Note that some of the sub blocks may have zero size. DS-13-1302 12 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME If (A,C) & (B, D) and (A,C) is observable, such an X is unique. X If (A,C) &X (B, D) and (A,C) is observable, any matrices E, F satisfying E ? X = F are of the form E = Q ? O (A,C), F = Q ? O (B, D) for some matrix Q. If (A,C) &X (B, D) and (B, D) is observable, then X has full column rank. If (A,C) &X (B, D) and (B, D) is observable and n = m, then X is invertible (hence (A,C) ? (B, D) and (A,C) is observable). If (A,C) h (B, D) and both are observable, then n = m and (A,C) ? (B, D). ed ite d It is well known that any two single output systems (A, cT ), (B, dT ) are equivalent if they are observable and A and B are similar (e.g. Casti [23]). A slightly more general result follows: Corollary 3.2. Suppose A ? B, and cT , dT are any two row vectors such that (A, cT ) is observable (single-output case). Then (A, cT ) &X (B, dT ) with X = (On (A, cT ))?1 On (B, dT ). py Proof. Since (A, cT ) is observable, On (A, cT ) is n О n and invertible. Since A and B have the same characteristic polynomial, if (14) holds for k = n it must hold for k = n + 1 (the final equation is the same linear combination of the previous n equations) and the result follows. ot Co The following result shows that every model of an output system contains a subsystem equivalent to that system. " # ! A1 , A2 Lemma 5. Suppose (A,C) &X (B, D) where (B, D) is observable. Then (A,C) ? , [C1 ,C2 ] for some A1 , A2 , A4 ,C1 ,C2 0, A4 such that (A1 ,C1 ) ? (B, D). tN Proof. Let (A,C) &X (B, D). Since (B, D) is observable, X has full column rank by corollary 3.1. Then sc rip " # I Q?1 , X =P 0 where P and Q are invertible matrices. Hence by Eq.(13) Ma nu " # ! i A1 , A2 ?1 h P P , C1 , C2 P?1 , 0, A4 (A,C) = where A1 = Q?1 BQ and C1 = DQ, and the result follows " # A1 0 A1 ? A2 = , 0 A2 ce pt ed Larger systems can be built up from sums of smaller models in parallel: define the direct sum of two matrices as and we have the following result: Ac Lemma 6. If (A1 ,C1 ) & (B1 , D1 ) and (A2 ,C2 ) & (B2 , D2 ), then (A1 ? A2 , [C1 ,C2 ]) & (B1 ? B2 , [D1 , D2 ]) (provided the matrices C1 and C2 have the same number of rows) and (A1 ? A2 ,C1 ?C2 ) & (B1 ? B2 , D1 ? D2 ). A generalization of this result to k-output systems forms a basis for the internal model principle: Lemma 7. Suppose A is in block-diagonal form with k identical blocks, i.e. ? ? A 0 1 k M ? . ? .. ?, A= A1 = ? ? ? 0 A1 DS-13-1302 13 (16) McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 where A by is ASME n О n , and 1 1 1 ? ? cT1 0 ? ? . .. ?, C=? ? ? T 0 ck ? dT11 ? . . D=? ? . dTk1 ? dT1k .. ? ? . ?, и и и dTkk иии .. . ed ite d where cT1 , . . . , cTk , dT11 , . . . , dTkk are row vectors. Then (A,C) & (A, D) if and only if (A1 , cTi ) & (A1 , dTij ) for all i, j, 1 ? i ? k, 1 ? j ? k. Proof. Setting Co where the sub-blocks Xi j are each of size n1 О n1 , (13) reduces to tN ot A1 Xi j = Xi j A1 cTi Xi j = dTij , py ? ? X11 и и и X1k ? . . ? . . . ... ? , X =? ? . ? Xk1 и и и Xkk and the result follows. sc rip We note in Lemma 7 that observability of all subsystems (A1 , cTi ) is equivalent to observability of (A,C). Also by Corollary 3.1, if (A, D) is observable then the systems are equivalent, so we have the following result: Ma nu Lemma 8. Suppose (A,C1 ), (A,C2 ) are systems where A is of block-diagonal form (16), and C1 ,C2 are k О n matrices. If (A,C1 ) is observable, (A,C1 ) & (A,C2 ). Moreover, if both systems are observable, they are equivalent. Proof. All that remains is to show that if (A,C1 ) is observable, there exists a vector cT1 such that (A1 , cT1 ) is observable. Since A consists of k copies of A1 , they have the same minimal polynomial and hence ed rank(On (A,C)) ? km1 , Ac ce pt where m1 is the degree of the minimal polynomial of A1 . The matrix has full rank only if m1 = n1 , i.e. A1 is cyclic (also termed nonderogatory [21]: each eigenvalue has geometric multiplicity 1). For some change of basis matrix P, such a matrix will have Jordan form ? J1 ? ?1 P A1 P = J = ? ?0 0 ? 0 0 .. ? . 0? ?, 0 Jr where the Jordan blocks ? ? ?i 1 0 ? . ? ? . Ji = ? ?0 . 0? 0 0 ?i DS-13-1302 14 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by to ASME correspond distinct eigenvalues ? . Set dT = [1, 0, . . . , 0]; it easily follows that (J , dT ) is observable, as is (J, dT ) where i i i i dT = [dT1 , . . . , dTk ]. Then set cT1 = dT P?1 , so we have (A1 , c1 ) ? (J, dT ) and hence (A1 , c1 ) is observable. With this choice of c1 , apply Lemma 7 with cTi = cT1 , 2 ? i ? k to obtain C such that (A,C) & (A,C1 ) and (A,C) & (A,C2 ). If (A,C1 ) is observable, (A,C1 ) ? (A,C) and hence (A,C1 ) & (A,C2 ). If both (A,C1 ) and (A,C2 ) are observable, they are both equivalent to (A,C), and hence (A,C1 ) ? (A,C2 ). This result will be used to show that any system (B, D) can be perfectly modelled by an observable system (A,C), where A is of the form (16). The matrix A1 is any matrix consisting of the largest cyclic component of B, i.e. a cyclic matrix satisfying the minimal polynomial of B. Lk A1 where A1 is the largest cyclic component of B. Then (A,C) & (B, D) ed ite d Lemma 9. Suppose (B, D) is observable, and A = if (A,C) is observable. Lp " # Bi Bi1 B1 ? 0 Bi2 ot for some matrices Bi1 , Bi2 , and hence we can augment the model to obtain Co py Proof. Without loss of generality, assume A1 is in Jordan form. By the cyclic decomposition theorem [20, 21], B ? i=1 Bi where the Bi are cyclic matrices in Jordan form such that B1 = A1 and each Bi satisfies the minimal polynomial of B1 . Each Bi can be augmented to a matrix similar to B1 , i.e. tN " # ! i B B01 h , D 0 & (B, D), 0 B02 sc rip such that " # p M B B01 B1 ? 0 B02 i=1 Ma nu and B01 = ed B02 = p M i=1 p M Bi1 Bi2 . i=1 ce pt Since (B, D) is observable, p ? k,3 and so the model can be augmented with k ? p copies of B1 to obtain Ac (A,C) & " # ! ! k?p h i M B B01 ? B1 , D 0 & (B, D), 0 B02 i=1 and the result follows. As an easy corollary of the previous result and relations R1 and R2, the main result follows, namely that any observable system that models k copies of the largest cyclic component of B, is a perfect model of (B, D), where k is the number of rows of C and D. This gives an alternative (and arguably simpler) proof of the Internal Model Principle [18]. 3 PROOF: Suppose p > k. Some eigenvalue ? has geometric multiplicity p, so rank B??I D ? n + k ? p < n, and hence (B, D) is not observable. Hence p ? k. DS-13-1302 15 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by4 ASME Theorem (Internal Model Principle). Suppose D has k rows, B is the largest cyclic component of B, and 1 "L # B1 A1 , 0 A2 k A? (17) for some matrices A1 and A2 . Then (A,C) models (B, D) for any matrix C such that (A,C) is observable. "L ed ite d Proof. Let X be a change of basis matrix such that X ?1 AX = A0 , where A0 is the right hand side of (17), and so # ! i B1 A1 h , C1 C2 , 0 A2 k (A,C) ? where [C1 ,C2 ] = CX. By relations R1 and R2, (A,C) & ( k B1 ,C1 ). Since (A,C) is observable, the subsystem ( k B1 ,C1 ) is observable, and by the previous lemma, models (B, D) for any D of compatible size. Hence, by the transitive property, (A,C) & (B, D). L Co py L sc rip tN ot We have proved that a robust model of an observable plant (and in particular, a robust PMIG for a given PIG) can be constructed from knowledge of only the minimal polynomial (i.e. the largest cyclic component) of the plant dynamics, and the number of outputs, by incorporating multiple copies of the largest cyclic component. This model is robust to any perturbations that do not affect the eigenstructure of the model, or the observability of the model, and moreover it may contain additional dynamics not present in the underlying plant. We now prove the converse of the internal model principle, namely that any ?universal? model of an observable plant necessarily contains a subsystem with multiple copies of the largest cyclic component. Theorem 5. Suppose (A,C) & (B, D) for all matrices D having the same number of rows as C, and let B1 be the largest cyclic component of B. Then A is equivalent to a matrix of the form (17). ed Ma nu " # B1 , B2 Proof. Let B1 be the largest cyclic component of B, and let B ? . As in the proof of Lemma 8, choose a row vector 0, B4 T dT1 such that (B Since (A,C) & (B, D) for all D, by the transitivity property and property R2, it follows ? ?1 , d1?) is ?observable. 0 ? ? ?? that (A,C) & ?B1 , ?dT1 ??. By Lemma 6, it follows that (A,C) & (?k B1 , ?k dT1 ), and by Lemma 5 it follows that (A,C) 0 contains a subsystem equivalent to (?k B1 , ?k dT1 ). Structural stability Any PMIG satisfying (11) (for some X22 ), that is subject to perturbations in some or all of its parameters will continue to produce unbiased input estimates provided (11) is still solvable. Such a PMIG is said to be structurally stable with respect to the perturbed parameters. Structurally stable parameters do not need to be known exactly for unbiased input estimates. The degree to which the solvability of (11) is stable with respect to the parameters of PIG and PMIG depends partly on the knowledge of the problem and, for multiple input systems, on the number of identical copies of the PIG within the PMIG. By way of example, consider a plant having two inputs. If the inputs are sinusoids of the same frequency but different amplitude and phase, then for unbiased input estimates the PMIG must have two undamped modes at this frequency. This PMIG does not, however, require a priori specification of the amplitude and phase of each mode. Ac 5.1 ce pt The amount of redundancy required in the model depends entirely on the knowledge of the plant; if the parameters of the plant dynamics and outputs are known exactly, a single copy of the exact model will suffice. DS-13-1302 16 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME For more complicated multiple-output signal generators, the canonical form of the system needs to be taken into consideration. As an example, the system ? ? # 0?0 " ? abc ? ?? ? ?? 0 0 0 ? , de f 000 ?? is equivalent to the system ? ? # 010 " ? y01 ? ?? ? ?? 0 0 0 ? , 100 000 ed ite d ?? 1 0 0 0 0 0 0 0 ? ? 0 " # ? 0? ? 1000 ? ?, ?, 1? 0 0 1 0 ? 0 ot 0 ?? 0 ?? ?? ?? 0 0 Co ?? py for some y depending on the parameters ?, a, b, c, d, e, f , and hence a model of this form is not robust with respect to all of these parameters. The model ? ? ? 0 x?1 ?x? ? ? 0 ? 2? ? ? ? = ? ?k1 ?k2 ?x?1 ? ? m 0 0 ?? ? ? ? 0 0 " # x1 ? ?x ? ? 0 0 ? F ? ? 2? ? ? 1 ? c2 ? ? ? + ? 1 ? ? ? ? ? F2 0 x? 1 m1 m1 ?c2 ?c3 1 x?2 0 m2 m2 ? ? x1 ?x ? ? 2? y = ? ?. ?x?1 ? x?2 1 0 0 1 Ma nu k2 ?c2 m1 m1 ?k2 ?k3 c2 m2 m2 (18) ed x?2 1 k2 m2 sc rip tN which is observable and contains two copies of the maximal cyclic component, is structurally stable with respect to the parameters ?, a, b, c, d, e, f . The examples that follow illustrate input estimation and structural stability when the Plant Model and PMIG are exact. These examples use the forced mass-spring-damper system shown in Figure 3. The system parameters are shown in Table 1. The dynamics of the system are described by Ac ce pt Further we assume that the measurement noise and process noise are both zero-mean Gaussians with covariance R and Q respectively. After substituing the parameters in Table 1 the plant is represented by A1 , B1 ,C1 matrices shown below. ? ? 0 0 1 0 ? 0 0 0 1? ? ? A1 = ? ? ??10 5 ?1.5 1.5? 5 ?10 1.5 ?3 ? ? 0 0 ?0 0? ? ? B1 = ? ? ?10 0 ? 0 10 C1 = I4О4 DS-13-1302 17 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Value k1 0.5 N/m k2 0.5 N/m k3 0.5 N/m m1 0.1 kg m2 0.1 kg c2 0.15 N.s/m c3 0.15 N.s/m py Parameter ed ite d Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME sc rip tN ot Co Table 1: Parameters of the double mass-spring-damper system. Ma nu Fig. 3: The double mass-spring-damper system. h i A2 = 0 " # 1 C2 = 0.5 ce pt ed Example 1. Step Input Consider the double mass-spring-damper system. Suppose the driving forces are unknown steps, F2 but the ratio between them is known. Suppose = 0.5, then PIG is represented using A2 and C2 matrices. F1 Ac Note that since x?2 = A2 x2 + w2 and A2 = [0], x2 (0) must be set to a non-zero value in order to generate a step input. In this example x2 (0) = 1. Assuming exact Plant Model specifying the PMIG so that it has a single eigenvalue at the zero consistent with the structure of a step input gives. A?1 = A1 B?1 = B1 C?1 = C1 A3 = A2 C3 = C2 . In PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG. Additionally we set the Q = I5О5 О 10?3 and R = 5 О I4О4 О 10?3 . Simulating the system response produces the results shown in Figure 4, with the key point to note being that the inputs converge to their correct values when the system transients have decayed. DS-13-1302 18 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME 3 4 2.5 3 4 3 2 2 1.5 2 1 1 0 0.5 1 -1 0 0 3 6 9 12 15 -0.5 0 3 6 9 12 15 time (s) time (s) 1 3 6 9 12 0 15 0 3 6 time (s) 0 0 12 15 0 3 6 Ma nu 9 time (s) 15 0 0 3 6 9 12 time (s) ot 0.5 6 12 tN 0.5 sc rip 1 3 9 time (s) 1 0 12 9 12 15 time (s) ed Fig. 4: Response of PIG and PMIG for step disturbance in double mass-spring-damper system. Ac ce pt Example 2. Compound Sinusoid-step Input Consider the double mass-spring-damper system. Suppose the driving forces in Figure 3 are a combination of two sinusoids function with known frequencies ?1 , ?2 and a step function. The magnitude and phase of these sinusoids and the magnitude of the the step function are unknown to the model, but the ratio between components in F1 and F2 are known. More specifically if the inputs are written as F1 = a1 u(t) + b1 sin(?1t + ?1 ) + c1 sin(?2t + ?2 ) F2 = a2 u(t) + b2 sin(?1t + ?1 ) + c2 sin(?2t + ?2 ) then a2 b2 c2 , , are known. a1 b1 c1 DS-13-1302 15 0.5 Co 0 9 py 0.5 -1 -2 6 1 1 0 3 1.5 3 2 0 time (s) 1.5 4 -2 ed ite d 0 19 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use 15 Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by example ASME suppose that PIG is represented by A and C matrices shown below. In this 2 2 ed ite d ? ? 0 0 0 0 0 ? ? ?0 0 ?1 0 0 ? ? ? ? A2 = ? ?0 ??1 0 0 0 ? ? ? ?0 0 0 0 ?2 ? 0 0 0 ??2 0 " # 1 0.5 0.5 1 ?1 C2 = 0.5 0.5 ?0.5 1 1 A?1 = A1 B?1 = B1 py h iT The initial state x2 (0) determines the phase and magnitude of the sinusoids in each input. We set x2 (0) = 0.5 0.25 0.25 0.5 0.5 . Additionally we choose f1 = 1Hz, f2 = 0.2Hz and set ?1 = 2? f1 , ?2 = 2? f2 . Assuming exact Plant Model and PMIG gives C?1 = C1 A3 = A2 C3 = C2 . tN ot Co In the PMIG the initial state is set to zero i.e. x3 = 0 so that the model is unaware of the actual input generated by PIG. Additionally we set the Q and R to the values of the previous example. The the system shown in Figure 2 is simulated and the responses are presented in Figure 5. The estimated inputs are seen to track the inputs after initial transients have decayed. Imperfect plant model: error bounds The results of the previous section showed that the asymptotic un-biasedness of the filter is not generally preserved under perturbations in parameters relating to the eigenvalues of the plant and signal. In practical applications, the parameters are never known exactly and the model is likely to capture only an approximation of the dynamics of interest, hence model error must be taken into consideration. This section shows that (i) the estimated bias is approximately proportional to the parameter error and (ii) the relative bias is bounded. The following result gives an estimate of the bias under the assumption that the model is ?close? to an exact model. Ma nu sc rip 6 0 , D0 ) is in error, such that Theorem 6. Suppose that the system model (A0L ,CM M A0L = AL + ?AL , D0M = DM + ?DM , pt ed 0 CM = CM + ?CM , Ac ce where (AL ,CM , DM ) is a model producing unbiased estimates of the system (AP , BP ,CP , DP ) governed by (3). Then the bias E[z0 ] of the imperfect system model satisfies E[z0 ] ? ?(?DM X + (DM + ?DM )?X)E[x] (19) where X satisfies (6) and ?X is the unique solution to ?XAP ? AL ?X = ?AL X. (20) Proof. The imperfect estimate satisfies x?? 0 = (AL + ?AL )x?0 + BP x DS-13-1302 20 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME 4 4 3 3 2 2 2 1 1 0 0 0 -1 -1 -2 -2 4 -2 0 3 6 9 12 15 -3 0 3 6 9 12 15 time (s) time (s) 2 12 -2 -2 15 0 3 6 1 1 0.5 0.5 0 0 6 9 12 time (s) -0.5 15 0 3 Ma nu 3 15 0 3 6 9 12 15 time (s) tN 1.5 12 sc rip 1.5 0 9 time (s) time (s) -0.5 15 py Co 9 ot 6 12 -1 -1 -4 3 9 0 0 0 6 1 -2 -6 3 2 1 0 0 time (s) 4 2 -6 ed ite d -3 -4 6 9 12 15 time (s) Fig. 5: Response of PIG and PMIG for combined sinusoids disturbance in double mass-spring-damper system. z0 = DP x ? (DM + ?DM )x?0 . pt ed and hence z0 satisfies Ac ce For the perfect model, apply Lemma 1 to obtain DM X = DP , where X satisfies XAP ? AL X = BP . Let ?X be the solution of (X + ?X)AP ? (AL + ?AL )(X + ?X) = BP , DS-13-1302 21 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME which simplifies to (20). Defining ? 0 = x?0 ? (X + ?X)x, we have ???0 = (AL + ?AL )??0 ? (X + ?X)w. The resulting bias E[z0 ] then satisfies E[z0 ] = (DP ? (DM + ?DM )(X + ?X))E[x] ? (DM + ?DM )E[??0 ] ed ite d = ?(?DM X + (DM + ?DM )?X)E[x] ? (DM + ?DM )E[??0 ]. Since ? 0 is stable, E[??0 ] = 0 and the result follows. Co py When the imperfect model contains multiple copies of an approximation of the maximal cyclic component of the system, the model error can be taken to be entirely due to the error in the eigenstructure, i.e. ?CM = 0, ?DM = 0, ?AL = ?AM , where ?AM is the error of the model relative to a perfect model which contains suitably replicated copies of the true maximal cyclic component of the system. The above result, although it requires the solution of the matrix equation (20), does give a measure of the error of the input estimate, modulo noise, relative to the plant state: Corollary 6.1. The bias relative to the plant state satisfies tN ot kE[u ? u?]k kE[z]k = ? k?DM k и kXk + (kDM k + k?DM k) k?Xk . kE[x]k kE[x]k sc rip Unfortunately the result does not give an error estimate relative to the input signal; for example if the input signal is sinusoidal then the relative error can be infinite at the instants when the input is zero. The next result provides an estimate of the bias when the parameter errors are large. Ma nu Lemma 10. For the plant and model governed by (3) with closed loop stability, the bias E[z] is O(t n eat ) for some a ? 0 and n ? 0. Moreover, if the plant state is bounded, with kE[x]k ? M, then E[z] is bounded with kE[z]k ? M kDP k +C1 M kDM k max(1, kBP k /х) for some positive constants C1 , х depending only on AL . ed Proof. The solution to (3) is pt x(t) = ?P (t)x(0) + 0 ?P (t ? ?)w(?)d?, (21) = exp(tAP ), x?(t) = ?L (t)x?(0) + ?L (t) (22) Rt 0 ?L (t ? ?)BP x(?)d?, (23) = exp(tAL ). (24) Ac ce ?P (t) Rt The matrix AL is stable, and thus k?L (t)k ? C1 e?хt for some positive constants C1 , х. This gives us the bound kE[x?(t)]k ? C1 e?хt kx(0)k +C1 kBP k DS-13-1302 22 Z t e?х(t??) kE[x(?)]k d?. 0 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME It easily follows that if kE[x(t)]k is O(t n eat ) as t ?? ? for some a > 0 then kE[x?(t)]k is also O(t n eat ). The bias E[z] is a linear combination of E[x?] and E[x] and thus of the same order. For the bounded case kE[x(t)]k ? M, we have kE[x?(t)]k ? C1 Me?хt +C1 M kBP k (1 ? e?хt )/х = C1 M kBP k /х + e?хt C1 M(1 ? kBP k /х) ? max(C1 M,C1 M kBP k /х), ed ite d and applying the triangle inequality to z = DP x + DM x?, py the result follows. sc rip tN ot Co The sensitivities can be difficult to calculate symbolically, so we give a numerical example here. Ma nu Fig. 6: The block diagram of the mass spring damper system. Example 3. Consider the mass spring damper system in Figure 6 with step-plus-sine and parameters w = 1, m = 1, k = 1, c = 10. Suppose that model parameters are in error, with k? = k + ?k c? = c + ?c, ed w? = w + ?w m? = m ce pt where ?w, ?k, ?c are small. For the stabilizing gains " Ac Kb = ? ? 0 ? ? Kc = ?25? , 50 # 0 100 the tracking error satisfies ku ? u?k ? (8.77 |?w| + 0.21 |?c| + 1.41 |?k|) kxk . Example 4. Step-Sinusoid Input with Imperfect Model Consider the double mass-spring-damper system shown in Figure 3. Suppose that the Plant Model and PMIG contain errors as shown in Table 2. DS-13-1302 23 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME Parameter Plant Parameter Plant Model 0.5 N/m k?1 0.55N/m k2 0.5 N/m k?2 0.55N/m k3 0.5 N/m k?3 0.55N/m m1 0.1 kg m?1 0.11kg m2 0.1 kg m?2 0.11kg c2 0.15 N.s/m c?2 0.165 N.s/m c3 0.15 N.s/m c?3 0.165 N.s/m Parameter PIG Parameter PMIG f 1Hz f? 1.1Hz py ed ite d k1 Co Table 2: Parameters of the real system versus the system model for the double mass-spring-damper shown in Figure 3. B?1 6= B1 A3 6= A2 tN A?1 6= A1 ot Therefore C?1 = C1 . sc rip The initial state of the PIG, x2 (0), determines the real input. Specifically x2 (0) determines the phase of sinusoid component, h iT the magnitude of the sinusoid component, and the size of step input component. In this example x2 (0) = 0.25 0.25 1 . By substituting the model parameters into Equation 18 the following matrices are calculated for Plant Model. ? 0 0 1 0 ? ? 0 0 0 1 ? ? A?1 = ? ? ??10 5 ?1.5273 1.5273 ? 5 ?10 1.5273 ?3.0545 ? ? 0 0 ? 0 ? 0 ? ? B?1 = ? ?. ?9.0909 ? 0 0 ?9.0909 pt ed Ma nu ? ce C?1 = C1 = I4О4 Ac Furthermore PMIG is given by ? ? ? ? o 2? f? 0 0 6.9115 0 ? ? ? ? A3 = ??2? f? 0 0? = ??6.9115 0 0? 0 0 0 0 0 0 " # 0.5 0.5 1 C3 = C2 = . 0.5 ?0.5 0.5 Using the above matrices yields the results shown in Figure 7. To verify that the input estimation bias satisfies CorolDS-13-1302 24 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) lary2017 6.1, by ASME kE[uu ? u?u]k kE[xx]k is plotted at every time step. The supremum of steady state estimate bias kE[uu ? u?u]k kuu ? u?uk = lim } t?? kx? xk kE[xx]k ed ite d sup{ is compared with k?DM k и kXk + (kDM k + k?DM k) k?Xk kE[[y1 , y2 , y?1 , y?2 , x2 ]]k ot and calculate T F 1 , F? F 2 ]T E[[[F 1 , F 2 ] ? [ F? = [y?1 , y?2 , y??1 , y??2 , x 3 ] Co kE[uu ? u?u]k = kE[xx]k T F 1 , F? F 2 ]T E[[[F 1 , F 2 ] ? [ F? py . In this example we use sc rip tN T F 1 , F? F 2 ]T [[F 1 , F 2 ] ? [ F? [y?1 , y?2 , y??1 , y??2 , x 3 ] for every point obtained from the simulation. Figure 8 shows that as the response converges to steady state, revealed. The figure shows kE[uu?u?u]k kE[xx]k is kE[uu ? u?u]k } = 0.1408. kE[xx]k Ma nu sup{ Additionally by substituting the matrices into Equation 4 and solving Equation 20 the following can be obtained. Ac ce pt ed ?DM = 0 2О7 " # 00001 1 1 DM = 0 0 0 0 1 ?1 1 ? ? ?0.1288 0.1349 ?0.0198 0.0376 0.2144 0.0855 0.0348 ? 0.0844 ?0.0888 0.0236 ?0.0383 ?0.0416 ?0.2628 ?0.0223? ? ? ? ? ? 0.2623 ?0.3413 ?0.0571 0.0207 0.1697 1.2240 0.0179 ? ? ? ? ?X = ? ??0.4445 0.5172 ?0.0104 0.0645 1.6985 0.0463 0.0513 ? . ? ? ??0.0415 0.0426 ?0.0088 0.0142 0.1291 1.2262 0.0124 ? ? ? ? 0.0384 ?0.0472 ?0.0049 ?0.0002 ?1.2068 0.0869 ?0.0011? 0.0167 ?0.0175 ?0.0017 ?0.0003 ?0.0522 0.0587 0.0863 Therefore k?DM k = 0, kDM k = 1.3229, kXk = 1, k?Xk = 2.1996 k?DM k и kXk + (kDM k + k?DM k) k?Xk = kDM k k?Xk = 2.9098 DS-13-1302 25 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 bythat ASME It follows kE[uu ? u?u]k ? k?DM k и kXk + (kDM k + k?DM k) k?Xk . kE[xx]k 3 3 2 2 1 1 4 3 ed ite d 2 1 0 0 0 -1 -1 -1 -2 3 6 9 12 15 0 3 6 time (s) 9 12 time (s) 1 1 0.8 0 0.6 -1 0.4 -2 0.2 -3 0 3 6 9 12 15 0 1.5 1 0.5 3 Ma nu 0.5 0 6 9 12 15 12 15 time (s) 1 6 9 12 0.6 0.4 0.2 15 time (s) time (s) 1 3 0.8 0 0 3 6 9 time (s) sc rip 0 0 1.2 ot 1.2 2 tN 3 -3 Co 4 15 py 0 0 -0.5 -0.5 0 3 6 9 12 0 3 6 9 12 15 time (s) ed time (s) 15 Conclusion We have shown that for unbiased input or disturbance estimation using the Kalman Inverse Filter, an exact model of the plant and a suitably replicated model describing the structure of the input signals is required. Moreover, we have shown that the model of the input signal generator is structurally stable with respect to all parameters that do not change the eigen structure of the PMIG. A new proof of the internal model principle has also been given. These results provide insight to the design of Kalman Inverse Filters. In several of the recent papers using this method, e.g. [6?13] the structural form of the PMIG is a chain of integrators without clear justification for these choices. The results of this paper show that this approach amounts to approximating the input signal as a piecewise polynomial whose order corresponds to the degree of the minimal polynomial of the PMIG matrix. A single integrator corresponds to a system with DS-13-1302 Ac 7 ce pt Fig. 7: Response of PIG and PMIG for combined step and sinusoid disturbance in double mass-spring-damper system with imperfect model parameters. 26 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c) 2017 by ASME Bias relative to state 10 bias relative to state upper bound 9 8 7 ed ite d 6 5 4 3 py 2 Co 1 0 4 6 8 time (s) 14 bias relative to state versus upper bound calculated by k?DM k и kXk + (kDM k + k?DM k) k?Xk. tN kE[uu?u?u]k kE[xx]k , 12 sc rip Fig. 8: 10 ot 2 Ma nu dynamics x? = 0 and thus would be consistent with an input having the structure of a constant (i.e. step) input. Similarly a chain of two integrators is consistent with the input having the structure of a ramp. For example, in [12], a double integrator is used to model the signal associated with each of the unknown forces acting on the wheels of a vehicle. Similarly a chain of three integrators yield a piecewise parabolic fit to an input. De Nicolao [16], in his work on Tikhonov regularization has similarly observed that chains of integrators can be associated with smoothing splines. ed For applications where the inputs have a higher structure that is known, e.g. harmonic inputs and/or disturbances, the insights provided by this work informs us that the the PMIG should be modelled as a resonator having the appropriate frequency. Where the signal is not a pure sinusoid, the PMIG might be extended to include harmonics of the base frequency. The efficacy of this approach has been demonstrated in [14]. If the system has multiple inputs suitable replication of the input structure is required. We note finally the possibility that where the structure of the input may be time varying, or not precisely defined, the PMIG could include appropriately rich structure to adapt to the input signal. 8 Ac ce pt From the point of view of the philosophical foundations of the design and in terms of guiding practical implementation of input-and-disturbance estimators, we have found the principles revealed in this paper extremely powerful in applications where we seek to estimate inputs and disturbances in systems using the KIF. Acknowledgements Kianoosh Soltani Naveh would like to acknowledge The University of Queensland Centennial Scholarship for the support of this work. References [1] Bayless, J. W., and Brigham, E. O., 1970. ?Application of the Kalman filter to continuous signal restoration?. Geophysics, 35(1), Feb, pp. 2?23. DS-13-1302 27 McAree Downloaded From: http://dynamicsystems.asmedigitalcollection.asme.org/ on 10/27/2017 Terms of Use: http://www.asme.org/about-asme/terms-of-use Ac ce pt ed Ma nu sc rip tN ot Co py ed ite d Journal of Dynamic Systems, Measurement and Control. Received August 05, 2013; Accepted manuscript posted October 25, 2017. doi:10.1115/1.4038267 Copyright (c)[2]2017 by ASME Crump, N. D., 1974. ?A Kalman filter approach to the deconvolution of seismic signals?. Geophysics, 39(1), Feb, pp. 1?13. [3] Dimri, V., 1992. Deconvolution and Inverse Theory. Elsevier Publishing Company. [4] Leite, L. W. B., and da. Rocha, M. P. C., 2000. ?Deconvolution of non-stationary seismic process?. Rev. Bras. Geof., 18(1), pp. 75?89. [5] Rocha, M. P., and Leite, L. W., 2003. ?Treatme

1/--страниц