Detection and Quantification of the Influence of Time-Variation in ...
Detection and Quantification of the Influence of Time-Variation in Closed Loop Frequency Response Function Measurements R. Pintelon, E. Louarroudi, and J. Lataire, Vrije Universiteit Brussel, dept. ELEC, Pleinlaan 2, 1050 Brussel, Belgium Tel. (+32) 2.629.29.44, Fax. (+32) 2.629.28.50, E-mail: [email protected]
Abstract—Recently a method has been developed for detecting and quantifying the time-variation in frequency response function (FRF) measurements using arbitrary excitations [1]. The following basic assumptions have been made: (i) the input is known exactly (generalized output error stochastic framework), and (ii) the time-variant system operates in open loop. The latter excludes any interaction between the time-variant system and the generator/actuator. In this paper we extend the results of [1] to noisy input-output observations (errors-in-variables stochastic framework) of time-variant systems operating in closed loop. Index Terms—frequency response function, best linear timeinvariant approximation, slowly time-varying, feedback, errorsin-variables, continuous-time, nonparametric modeling
I. I NTRODUCTION In all kinds of engineering applications frequency response function (FRF) measurements are a basic tool to obtain quickly insight into the dynamic behavior of complex systems [2], [3], [4]. Examples of complex systems that are studied via FRF measurements are (i) corrosion of metals where the time-variation is induced by the changing surface condition [5], (ii) airplane dynamics during take-off and landing which depend on the changing flight speed and height [6], and (iii) myocardium tissue impedance measurements that evolve during the cardiac cycle [7]. In practise, these real-life systems satisfy only approximately the linearity and time-invariance assumptions inherent to the concept of an FRF. While the influence of nonlinear distortions on FRF measurements using random excitations has been studied in detail in the literature [4], [8], [9], [10], much less is known about the impact of time-variations. Recently it has been shown that time-variations cause errors on FRF measurements that are correlated over the frequency [11]. Since the leakage errors have exactly the same property [4], [12], it is difficult to distinguish and quantify these two sources of errors [11]. Assuming that the input is known exactly and that the linear time-variant system operates in open loop, this problem has been solved in [1] by modeling nonparametrically the timevariation using multiple-input, single-output FRF estimates. The contribution of this paper is to extend the results of [1] to noisy input-output observations of time-variant systems operating in closed loop (see Figure 1).
Figure 1. Linear time-variant (LTV) system operating in feedback. The actuator and feedback dynamics are linear and time-invariant (LTI). r (t), w (t), u (t), and y (t) are, respectively, the reference signal, the output of the actuator, and the input and output of the LTV system.
Feedback in a measurement setup is due to either an explicit control action (see Fig. 1) or an interaction between an nonideal generator/actuator and the device under test. The latter is illustrated by the problem of measuring a time-varying impedance Z (s, t) [13] from noisy current i (t) and voltage v (t) observations (see Fig. 2). The relationship between the current through and the voltage across the impedance Ze (s) in Figure 2 is given by I(s) = (E(s)
V (s)) /Ze (s)
(1)
where I(s) = L{i(t)}, V (s) = L{v(t)}, and E(s) = L{e(t)}, with L{} the Laplace transform. Formula (1), with E(s) = R(s) (e (t) = r (t)), I(s) = U (s) (i (t) = u (t)), and V (s) = Y (s) (v (t) = y (t)), is exactly the feedback law in Figure 1, where the actuator and feedback dynamics are both equal to 1/Ze (s). Note that the feedback in the setup of Figure 2 is solely due to the non-zero output impedance of the voltage source. Hence, feedback is present in any measurement setup with a non-ideal generator/actuator. For example, a signal generator operating above 10 MHz has typically a 50 ⌦ output impedance, and a force actuator (shaker) for mechanical impedance measurements has a frequency dependent output impedance. It emphasizes the importance of handling timevariant systems operating in closed loop. The outline of the paper is as follows. First, the definition of the best linear time-invariant (BLTI) approximation of a linear time-variant system operating in feedback is given, and its properties are discussed (Section II). Next, it is shown in Section III how the BLTI approximation can be obtained from noisy input-output measurements using multiple-input, multiple-output nonparametric FRF estimators. Further, the
This is the author’s version of the article published as R. Pintelon, E. Louarroudi, and J. Lataire. Detection and quantification of the influence of time-variation in closed loop frequency response function measurements. IEEE Transactions on Instrumentation & Measurement, 62(4):853–863, 2013.
G (s, t) =
1 X
Gp (s) fp (t)
p=0
Figure 2. Measurement of a time-varying impedance Z (s, t) using a voltage source e (t) with output impedance Ze (s). v (t) and i (t) are the voltage across and the current through Z (s, t) respectively.
• •
t 2 (0, T )
(3)
with fp (t) a complete set of basis functions over the interval (0, T ). The excitation u (t) is either periodic noise or a random phase multisine with period length T . The LTI systems Gp (s) in (3) operate in steady state.
Under Assumption 1 it is shown in [11] that the direct frequency response function (FRF) estimate obtained via crossand autopower spectra equals the BLTI approximation (2) GBLTI (j!k ) = Figure 3. Best linear time-invariant (BLTI) approximation GBLTI (s) (bottom block diagram – frequency domain) of a linear time-variant (LTV) system G (s, t) operating in open loop (top block diagram – time domain). Gact (s) is the transfer function of the linear time-invariant actuator. r (t), u (t), and y (t) are, respectively, the reference signal, and the input and output of the LTV system. YTV (k) is the DFT spectrum of the difference between the actual output of the LTV system and the output of the BLTI approximation. It is uncorrelated with - but not independent of - the DFT spectrum U (k) of the input u (t).
case of time-variant actuator and/or feedback dynamics is handled in Section IV. The theory is illustrated on a simulation (Section V) and a measurement (Section VI) example. Finally, some conclusions are drawn in Section VII. II. D EFINITION OF THE B EST L INEAR T IME -I NVARIANT A PPROXIMATION First, the difficulty of estimating the best linear timeinvariant (BLTI) approximation of a linear time-variant system operating in closed loop using the open loop definition introduced in [1], [11] is revealed. Next, a generalized definition that can handle the closed loop case is given, and its properties are discussed. A. Open Loop Case If the system operates in open loop (see Fig. 3, top block diagram), then the BLTI approximation of the linear timevariant system over the time interval [0, T ] is defined as [1] 1 GBLTI (s) = T
ˆ
T
G (s, t) dt
(2)
0
where G (s, t) = L {g (t, t ⌧ )} is the time-varying transfer function, with L {} the Laplace transform w.r.t. ⌧ , and g (t, ⌧ ) the time-varying impulse response (i.e. response of the system to a Dirac impulse at time ⌧ ). Assumption 1. (Open loop case) •
The time-varying transfer function G (s, t) can be expanded in series as
E Y (k) U (k) n o 2 E |U (k)|
(4)
where the expected value E {} is taken w.r.t. the random realization of u (t), and where X (k), X = Y or U , is the discrete Fourier transform (DFT) of x (nTs ), x = y or u, N 1 1 X X (k) = p x (nTs ) e N n=0
j2⇡kn/N
(5)
with Ts the sample period and T = N Ts . Hence, the actual output Y (k) of the LTV system can be written as (see also Fig. 3, bottom block diagram) Y (k) = GBLTI (j!k ) U (k) + YTV (k)
(6)
where the residual YTV (k) is uncorrelated with – but not independent of – the input U (k). It shows that under Assumption 1, GBLTI (s) (2) is the best (in least squares sense) LTI approximation of the LTV system. If definition (2) is applied to the closed loop case (see Fig. 4), then we obtain biased estimates of the BLTI approximation. Indeed, combining (4) and (6) gives E Y (k) U (k) E YTV (k) U (k) n o = GBLTI (j!k ) + n o 2 2 E |U (k)| E |U (k)|
where E YTV (k) U (k) 6= 0 due to the feedback loop. This is the key difficulty of estimating the BLTI approximation of an LTV system operating in closed loop. In the next section we present a solution to this problem. B. Closed Loop Case Consider a linear time-variant (LTV) system operating in feedback (see Fig. 4, top block diagram). Inspired by the indirect method for measuring the FRF of an LTI system operating in closed loop [19], [20], we propose the following generalized definition for the BLTI approximation GBLTI (s) =
1 T 1 T
´T 0
´T 0
Gry (s, t) dt Gru (s, t) dt
(7)
where the residual YTV (k) is uncorrelated with – but not independent of – the reference R (k) (proof: see Appendix A). It shows that under Assumption 2, the best linear time invariant approximation GBLTI (s) (7) of the time-variant system operating in closed loop, is the ratio of the best (in least squares sense) LTI approximations of two time-variant open loop systems: that from reference to output to that from reference to input. Properties of the BLTI Approximation • •
Figure 4. Best linear time-invariant (BLTI) approximation GBLTI (s) (bottom block diagram – frequency domain) of a linear time-variant (LTV) system G (s, t) operating in feedback (top block diagram – time domain). Gact (s) and Gfb (s) represent the linear time-invariant actuator and feedback dynamics, respectively. r (t), w (t), u (t), y (t), and q (t) are, respectively, the reference, the output of the actuator, the input and output of the LTV system, and the output of the feedback branch. YTV (k) is the DFT spectrum of the difference between the actual output of the LTV system and the output of the BLTI approximation. It is uncorrelated with - but not independent of the DFT spectrum R (k) of the reference signal r (t).
where Gru (s, t) and Gry (s, t) are the time-varying transfer functions from the reference r (t) to, respectively, the input u (t) and the output y (t). The 2 ⇥ 1 time-varying transfer function [Gry (s, t) Gru (s, t)]T , with Grz (s, t) = z = [y u]T , can be expanded in series as Grz (s, t) =
1 X
Grz,p (s) fp (t)
p=0
• •
t 2 (0, T )
(8)
with fp (t) a complete set of basis functions over the interval (0, T ). The reference signal r (t) is either periodic noise or a random phase multisine with period length T . The 2 ⇥ 1 LTI systems Grz,p (s) in (8) operate in steady state.
Under Assumption 2 it is shown in Appendix A that the indirect frequency response function (FRF) estimate obtained via crosspower spectra equals the BLTI approximation (7) GBLTI (j!k ) =
f0 (t) = 1 and
1 T
ˆ
T
fp (t) dt = 0 for p > 0. (11)
0
With this choice it follows that GBLTI (s) =
Gry,0 (s) Gru,0 (s)
(12)
(proof: see Appendix B).
Assumption 2. (Closed loop case) •
•
For LTV systems operating in open loop, the generalized definition (7) reduces to (2) (proof: see Appendix B). The BLTI approximation (7) depends on the feedback dynamics Gfb (s), but is independent of the actuator characteristics Gact (s) (proof: see Appendix B). For example, the BLTI approximation of the time-varying impedance in Figure 2 depends on the value of the resistor R. Without any loss in generality, the basis functions fp (t) in (8) can always be chosen such that
E Y (k) R (k) E U (k) R (k)
Applying Mason’s rule (see [21]) to the bottom block diagram of Figure 4, it can easily be verified that, under Assumption 2, the input-output DFT spectra can be written as ˜TV (k) U (k) = UR (k) + U
(13)
Y (k) = GBLTI (j!k ) UR (k) + Y˜TV (k)
(14)
with UR (k) that part of the input that is correlated with the reference R (k), UR (k) =
(10)
Gact (j!k ) R (k) 1 + GBLTI (j!k ) Gfb (j!k )
˜TV (k), Y˜TV (k) the time-varying input-output contribuand U tions that are uncorrelated with R (k)
(9)
where the expected value E {} is taken w.r.t. the random realization of r (t). Moreover, the actual output Y (k) of the LTV system can be written as (see also Fig. 4, bottom block diagram) Y (k) = GBLTI (j!k ) U (k) + YTV (k)
Interpretation of the Output Residual
˜TV (k) U
=
Y˜TV (k)
=
Gfb (j!k ) YTV (k) (15) 1 + GBLTI (j!k ) Gfb (j!k ) 1 YTV (k) (16) 1 + GBLTI (j!k ) Gfb (j!k )
Since the input U (k) (13) is also subject to time-variation ˜TV (k), the output time-variation Y˜TV (k) in (14) cannot U be used to quantify the time-variation of the closed loop
z (t) = [ y (t) u (t) ]T . Proceeding in this way the singleinput, single-output closed loop problem is transformed into a single-input, two-output open loop problem for which the results of [1] apply. First, we show that there exists a multipleinput, two-output linear time-invariant equivalent model of the single-input r (t), two-output z (t) slowly time-varying system. Next, the relationship between the DFT spectra of the reference and the input-output signals is established. Finally, an algorithm for estimating the BLTI approximation of an LTV system operating in feedback is presented. A. Linear Time-Invariant Equivalence of a Linear Slowly Time-Varying System A linear time-variant system operating in closed loop is called slowly-time varying if the time-varying transfer function Grz (s, t) of the open loop system from reference r to inputoutput z = [y u]T is slowly time-varying: Definition 3. A linear time-variant system operating in closed loop is slowly time-varying over the interval [0, T ] if the time-varying transfer function Grz (s, t) (8) can be written as Grz (s, t) =
Nb X
Grz,p (s) fp (t)
p=0
(18)
t 2 (0, T )
with fp (t) a complete set of basis functions over the interval (0, T ) satisfying (11).
Figure 5. Linear time-invariant (LTI) equivalent open loop model (bottom) of a linear slowly time-varying system operating in feedback (top). Middle block diagram: single-input r (t), two-output z (t) = [ y (t) u (t) ]T linear time-variant open loop equivalent representation of a single-input u (t), single-output y (t) linear slowly time-varying system operating in feedback (top block diagram). Bottom block diagram: multiple-input r (t), r1 (t), . . . , rNb (t), two-output z (t) = [ y (t) u (t) ]T linear time-invariant equivalent model (bottom block diagram) of a single-input r (t), two-output z (t) = [ y (t) u (t) ]T linear time-variant open loop system (middle block diagram).
LTV system. However, multiplying (15) with GBLTI (j!k ) and subtracting the result from (16) gives exactly YTV (k) Y˜TV (k)
˜TV (k) = YTV (k) GBLTI (j!k ) U
(17)
which is consistent with Eq. (43) of Appendix A. It shows that the output residual YTV (k) of the BLTI approximation is the time-varying part of the output Y˜TV (k) minus the input time˜TV (k) filtered through the BLTI approximation variation U GBLTI (j!k ). III. E STIMATING THE BLTI A PPROXIMATION OF S LOWLY T IME -VARYING S YSTEMS O PERATING IN F EEDBACK The basic idea consists in modeling the linear time-variant (LTV) system from the reference r (t) to the input and output
Figure 5 shows a block diagram of the slowly timevarying system (18). If the time-variant system operates in open loop (Gfb (s) = 0 in Fig. 5), then Grz (s, t) = Gact (s) [ G (s, t) 1 ]T (see Appendix B, Eqs. (44) and (45)), and Definition 3 reduces to that for open loop systems [1]. Using Theorem 2 of [1], it follows that the output z (t) = [ y (t) u (t) ]T of the single-input r (t) open loop slowly time-varying system (18) can be calculated via an equivalent multiple-input rp (t) = r (t) fp (t), p = 0, 1, . . . , Nb , two-output z (t) = [ y (t) u (t) ]T open loop linear timeinvariant model (see also Fig. 5)
z (t)
=
L
1
(N b X p=0
=
L
1
(N b X p=0
Grz,p (s) L {r (t)} fp (t) )
Hrz,p (s) L {rp (t)}
)
(19) (20)
with L 1 {} and L{}, respectively, the inverse Laplace transform and the Laplace transform operators. Choosing Legendre polynomials Pp (2t/T 1), t 2 [0, T ], satisfying (11) as basis functions fp (t) (see [14]), the relationship between Grz,0 (s), the first term in the sum (19), and the transfer functions Hrz,p (s), p = 0, 1, . . . , Nb , of the equivalent LTI model (20), is given by
Grz,0 (s) = Hrz,0 (s) + j
Nb 2
k
4 X T 2 p=0
2 T
j
(2) 2p Hrz,2p
Nb 1 2
X
k
(1)
Hrz,2p+1 (s) + . . .
p=0
(s) + O T
3
(21)
2
with 2p = 1.5+2.5 (p 1)+(p 1) , bxc the largest integer (m) smaller than or equal to x, Hrz,p (s) the mth derivative of Hrz,p (s) w.r.t. s, and where O(T m ) with m > 0 means that limT !1 T m O(T m ) < 1 (proof: apply Theorem 3 of [1] to Eqs. (20) and (19)). The bias term O(T 3 ) depends (m) on the higher order derivatives Hrz,p (s), with m, p > 3 and m 6 p. Finally, combining Eqs. (12) and (21) establishes the relationship between the BLTI approximation (7) of the linear slowly time-varying system operating in closed loop and the transfer functions Hrz,p (s), p = 0, 1, . . . , Nb , of the open loop LTI equivalent model. It is the basis of the nonparametric estimation procedure developed in Section III-C. B. Relationship Between the DFT Spectra The relationship between the sampled known reference signal and the sampled noiseless input-output signals is established for band-limited reference signals . Assumption 4. The reference signal r (t) is band-limited: the power spectral density of the reference Srr (j!) is zero for f > fB , where fB < fs /2, with fs = 1/Ts the sampling frequency. Under Assumption 4, the reference R (k) and input-output Z (k) = [ Y (k) U (k) ]T DFT spectra (5) of a slowly timevarying system operating in feedback (see Definition 3) are related by Z (k) =
Nb X
Hrz,p (j!k ) Rp (k) + TZ (j!k )
(22)
p=0
with Rp (k) = DFT{r (t) fp (t)}, and where TZ (j!k ) models the sum of the spectral leakage and the residual alias errors (proof: apply Theorem 5 of [1] to Eq. (20)). C. Estimation Procedure In practise the nonparametric estimate of the best linear time-invariant (BLTI) approximation is calculated from known reference r (nTs ) and noisy input-output u (nTs ), y (nTs ) samples, n = 0, 1, . . . , N 1 (see Fig. 6). It is assumed that the input-output measurement errors nu (t), ny (t) in Figure 6 can be modelled as stationary (mutually correlated) zero mean, filtered (band-limited) white noise [4], [15]. The reference signal r (t) satisfies Assumption 4, and the input-output signals u (t) and y (t) are lowpass filtered before sampling.
Figure 6. Known reference r (t), noisy input-output u (t), y (t) observations of a linear time-variant system operating in feedback. nu (t) and ny (t) are the (mutually correlated) input-output measurement errors, and u0 (t), y0 (t) are the unknown true input-output signals.
The procedure for estimating the BLTI approximation GBLTI (j!k ) and the output residual YTV (k) consists of the following four steps: 1) Nonparametric estimation of Hrz,p (j!k ): From Eq. (22), it follows that the Nb +1 FRFs Hrz,p (j!k ) can be estimated using nonparametric multiple-input, multipleoutput (MIMO) methods. The local polynomial method ([1], [16], [17]) is used here because of its excellent leakage (TZ (j!k ) in Eq. (22)) suppression property, and because it estimates the FRFs at the full frequency resolution 1/T of the experiment duration T . 2) Estimation of Nb : The number Nb of parallel branches modeling the time-variation in (22) is determined as follows. Starting with Nb = 1, the number Nb is ˆ rz,N (j!k ) is increased till the estimated 2 ⇥ 1 FRF H b zero within its uncertainty (the LPM method calculates ˆ rx,N (j!k ) | ⇠ also the covariance of the estimates): |H b ˆ std(Hrx,Nb (j!k )) for x = y, u. The number of significant parallel branches is then Nb 1. 3) Nonparametric estimation of Grz,0 (j!k ): First, a nonparametric estimate of the 2 ⇥ 1 FRF Grz,0 (j!k ) is obtained via (21) ˆ rz,0 (j!k ) = H ˆ rz,0 (j!k ) + 2 G T
j
Nb 1 2
X
k
ˆ (1) H rz,2p+1 (j!k )
p=0
(23) where the first order derivatives are replaced by central differences ˆ rz,m (j!k+1 ) H (1) ˆ rz,m H (j!k ) = j!k+1
ˆ rz,m (j!k H j!k 1
1)
(24) Next, using (12), an estimate of the BLTI approximation is calculated as ˆ ˆ BLTI (j!k ) = Gry,0 (j!k ) G ˆ ru,0 (j!k ) G
(25)
Since the expected value of (23) equals the true value within an O T 2 bias error if at least a second order polynomial approximation is used in the LPM
ˆ BLTI (j!) (25) is method (proof: see [1]), the bias of G asymptotically (T ! 1) given by n o ˆ BLTI (j!k ) E G GBLTI (j!k ) = O T 2 (26)
(proof: see Appendix C). 4) Nonparametric estimation of YTV (k): First, the output ˜TV (k) ]T (see residual Z˜TV (k) = [ Y˜TV (k) U Appendix A, Eq. (40)) of the BLTI approximation from reference r (t) to input-output z (t) = [ y (t) u (t) ]T is estimated as
Figure 7. Linear time-invariant equivalent open loop model of a linear timeinvariant system with time-variant actuator and/or feedback dynamics.
Nb
X ˜TV (k) = ˆ rz,p (j!k ) Rp (k) Zˆ H ⇣
p=1
ˆ rz,0 (j!k ) G
···
⌘
ˆ rz,0 (j!k ) R (k) H
A. Time-Variant System (27)
where the second term in the difference has an O(T 1 ) contribution (see [1], Eq. (18)). Next, using (17), we obtain an estimate for the output residual YTV (k) of the LTV system operating in closed loop ˜TV (k) YˆTV (k) = Yˆ
ˆ ˆ BLTI (j!k ) U ˜TV (k) G
(28)
The bias of (28) is an O(T ) because the bias of the ˜TV (k) and G ˆ BLTI (j!k ) are both an O(T 2 ). estimates Zˆ 2
Notes •
Using (21), the O(T 2 ) bias (26) of the BLTI approximation estimate (25) can be reduced to an O(T 3 ). First, an improved estimate of Grz,0 (j!k ) is calculated as ˆ ˆ rz,0 (j!k ) = G ˆ rz,0 (j!k ) + · · · G j
Nb 2
k
4 X T 2 p=1
ˆ (2) 2p Hrz,2p
(j!k )
(29)
where the second order derivatives are replaced by second order central differences (see [1], [18]). Next, combining (12) and (29) gives an improved BLTI approximation estimate ˆ ˆ ˆ ˆ BLTI (j!k ) = Gry,0 (j!k ) G ˆ ˆ ru,0 (j!k ) G •
(30)
ˆ ˆ BLTI (j!k ) (25) by G ˆ BLTI (j!k ) (30) in Eq. Replacing G 2 (28) reduces the O(T ) bias of YˆTV (k) to an O(T 3 ). IV. T IME -VARIANT ACTUATOR AND / OR F EEDBACK DYNAMICS
In this section the results are extended to time-variant actuator and/or feedback dynamics. Two cases are distinghuised: the system is time-variant or time-invariant.
Eqs. (9) and (10) remain valid for time-variant actuator and/or feedback dynamics because the proof in Appendix A does not use the time-invariance of the actuator and the feedback loop. It shows that Figure 4 remains valid for timevariant actuator Gact (s, t) and feedback Gfb (s, t) transfer functions. However, the BLTI approximation (7) depends now on the time-variant actuator dynamics; even in open loop. B. Time-Invariant System If the system is time-invariant, G (s, t) = G (s), then one would expect that the BLTI approximation (7) equals G (s). This is not true because Gry (s, t) 6= Gru (s, t) G (s) [22]. However, in the sequel of this section we will show that the difference between GBLTI (s) and G (s) is an O(T 1 ). Define the linear time-invariant (LTI) approximation of the linear time-variant (LTV) system G (s, t) as GLTI (s) =
Hry,0 (s) Hru,0 (s)
(31)
where Hrz,0 (s) is the top branch in the LTI equivalent model of the LTV system (see Fig. 5). The output contribution of the time-varying branches in the LTI equivalent model is given by Z˘TV (k) =
Nb X
Hrz,p (j!k ) Rp (k)
(32)
p=1
Using (32), one finds the output residual YˇTV (k) of the LTI approximation (31) YˇTV (k) = Y˘TV (k)
˘TV (k) GLTI (j!k ) U
(33)
Applying the results of Figure 5 to the time-variant dynamics from reference r (t) to the input u (t), taking into account that the system is time-invariant G (s, t) = G (s), gives the linear time-invariant (LTI) equivalent open loop model shown in Figure 7. Shifting the LTI block G (s) in Figure 7 backwards in the parallel branches, transforms the block diagram in an LTI equivalent open loop model from reference to output. It shows that the dynamics Hry,p (s) are related to Hru,p (s) as Hry,p (s) = Hru,p (s) G (s)
(34)
Reference DFT spectrum
Poles (x) and zeroes (o) Amplitude (dB)
0.5 0.0 −0.5 −0.6
−0.4 −0.2 Real part
Combining (31)–(34), it can easily be verified that =
G (s)
(35)
=
0
(36)
if the system is time-invariant. Property (36) allows us to verify whether the system is time-invariant or not. If (36) is true within the noise standard deviation, then the system is timeinvariant and its frequency response function is calculated via Eq. (31). Otherwise, the system is time-variant and its BLTI approximation is estimated as explained in Section III-C. Since Grz,0 (s) = Hrz,0 (s) + O(T 1 ) (see Eq. (21)), it can be concluded from (12) and (31) that GBLTI (s) = GLTI (s) + O(T 1 ) (proof: follow the same lines of Appendix C). Hence, in the time-invariant case, the BLTI approximation equals the FRF within an O(T 1 ) bias error. V. S IMULATION E XAMPLE In this section we simulate a linear time-variant (LTV) system operating within a linear time-invariant (LTI) feedback loop and driven by an LTI actuator. Following the procedure of Section III-C, the best linear time-invariant (BLTI) approximation is estimated from noisy input-output observations of the LTV system (see the setup of Figure 6). The simulated LTV system is a third order differential equation with time-varying coefficients 3 X p=0
(p)
ap (t) y0 (t) =
3 X
(p)
bp (t) u0 (t)
0 10 20 0
0
Figure 8. Poles (⇥) and zeroes ( ) of the frozen transfer function (38) – simulation example. The real zero at 4 is fixed and not shown.
GLTI (s) YˇTV (k)
10
(37)
p=0
where x(p) (t) = dp x (t) /dtp , x = u0 , y0 . The poles and zeroes of the corresponding frozen transfer function P3 p p=0 bp (t) s Gf (s, t) = P3 (38) p p=0 ap (t) s
(transfer function obtained from (37) by fixing the time in ap (t) and bp (t)) are shown in Fig. 8 for t 2 [0, T ], with T = 3840 s. The system operates in closed loop (see Fig. 6) with Gfb (s) = 1, and the actuator transfer function Gact (s) is a sixth order Chebyshev filter with a passband ripple of 6 dB and a cut-off frequency of 0.3125 Hz. The reference signal r (t) is a zero mean band-limited random excitation with standard deviation equal to one satisfying Assumption 4 with fB = 0.2930 Hz (see Fig. 9, top plot). Using numerical integration (ODE45 solver in MATLAB®),
0.1 0.2 Frequency (Hz)
0.3
Input DFT spectrum
Output DFT spectrum
0
Amplitude (dB)
−1.0
Amplitude (dB)
Imaginary part
1.0
20 40 60 0
0.1 0.2 Frequency (Hz)
0.3
0 20 40 60 0
0.1 0.2 Frequency (Hz)
0.3
Figure 9. DFT spectra (black) of the known reference (top), noisy input (bottom left), and noisy output (bottom right) signals – simulation example. ˆ ˜TV (k) and Top light gray lines: estimated input-output time-variation U ˆ Y˜TV (k) (27); middle dark gray lines: estimated input-output leakage errors TˆU (j!k ) and TˆY (j!k ) (22); bottom medium gray lines: estimated inputoutput noise variances; and black dashed lines: true input-output noise variances.
the response y0 (t) is calculated for t 2 [0, T ] at the rate fs = 2 Hz, giving N = 7680 samples. The noiseless inputoutput samples u0 (t) and y0 (nTs ) are disturbed by white noise nu (nTs ) and ny (nTs ) (see Fig. 6), and the black dashed lines in the bottom plots of Fig. 9 show var (DFT(nx (nTs ))), x = u, y, as a function of the frequency. The multiple-input, two-output (MITO) model (22), with fp (t) Legendre polynomials satisfying (11), is estimated from the N known reference r (nTs ), noisy input u (nTs ), and noisy output y (nTs ) samples using the procedure of Section III-C (local polynomial method [16] with a sixth order local polynomial approximation of the FRFs and the leakage term, and six degrees of freedom for the variance estimates). Figure 10 shows the estimated dynamics of the time-varying branches ˆ rz,p (j!k ), p > 1, for Nb = 2. It can be seen that the H estimates coincide – within their uncertainty – with the true values. The latter are obtained by estimating (22) from the noiseless samples. Increasing Nb further makes no sense since ˆ rx,3 (j!k ) | ⇠ std(H ˆ rx,3 (j!k )) for x = u, y. |H Next, a nonparametric estimate of the best linear timeinvariant (BLTI) approximation of (37) is obtained via Eqs. (25) and (30). The results are shown in Figure 11, left plot. It can be seen that the bias of the estimated BLTI approximation ˆ BLTI (j!k ) is smaller than its uncertainty. G Further, the time-varying parts of the input-output DFT spectra Z˜TV (k) and the output residual YTV (k) of the BLTI approximation are estimated via, respectively, (27) and (28). The results are shown in Figures 9 and 11. From Figure 9 it can be seen that the time-variation (top light gray lines) is well above the estimated noise level (bottom medium gray lines), and is of the same order of magnitude as the leakage errors (middle dark gray lines). Note that the input time-variation-to-signal ratio is much larger than that at the
H
(jω )
ru,1
H
k
−80 0.1 0.2 Frequency (Hz)
−40
−80 −100 0
0.3
Hru,2(jωk)
0.1 0.2 Frequency (Hz)
0.3
Hry,2(jωk)
−40 −60 −80 0.1 0.2 Frequency (Hz)
−40
0.1 0.2 Frequency (Hz)
0.3
−40 −60 −80 −100 0
0.3
0
−80 0
−20 Amplitude (dB)
−20 Amplitude (dB)
GSA(jωk) 40
−60 Amplitude (dB)
Amplitude (dB)
Amplitude (dB)
−60
−100 0
k
−20
−40
−100 0
(jω )
ry,1
−20
0.1 0.2 Frequency (Hz)
0.3
Figure 12. Spectral analysis estimate GSA (j!k ) (9) of the best linear time-invariant approximation – simulation example. Black: estimate ˆ SA (j!k ); dark gray: variance estimate var(G ˆ SA (j!k )); light gray: differG ˆ SA (j!k ) G ˆ BLTI (j!k ) |; and black dashed line: var(G ˆ BLTI (j!k )). ence |G
Figure 10. Frequency response functions Hrz,p (j!k ), p = 1, 2, of the MITO equivalent model (22) – simulation example. Black: input (left) and ˆ ry,p (j!k ); light gray: difference ˆ ru,p (j!k ) and H output (right) estimates H between the estimated and the true value of Hrz,p (j!k ); and dark gray: ˆ rz,p (j!k ). estimated variance of H
G
(j
BLTI
k
)
Output residual Amplitude (dB)
Amplitude (dB)
50 0 50 100 0
0.1 0.2 Frequency (Hz)
0.3
0 20 40 60 0
0.1 0.2 Frequency (Hz)
0.3
Figure 11. Estimated best linear time invariant (BLTI) approximation (left) and its output residual (right) – simulation example. Left plot: the estimate ˆ BLTI (j!k ) (black), its variance var(G ˆ BLTI (j!k )) (dark gray), and its bias G ˆ ˆ ˆ |GBLTI (j!k ) GBLTI (j!k ) | (light gray). Right plot: comparison between ˜TV (k) (27) of the time-varying part Y˜TV (k) of the output the estimate Yˆ DFT spectrum (black), and the estimated output residual YˆTV (k) (28) of the BLTI approximation (light gray).
output below 0.15 Hz. From the right plot of Figure 11 it can be seen that YTV (k) is significantly larger than the time-varying part Y˜TV (k) of the output DFT spectrum in the band [0 Hz, 0.15 Hz]. Note that in this frequency band ˆ BLTI (j!k )| |G 1 (see Fig. 11, left plot). It nicely illustrates that a feedback loop suppresses the time-variation at those frequencies where the open loop gain is significantly larger than one. Finally, Figure 12 compares the proposed approach with the indirect spectral analysis estimate (12) obtained via (22) with Nb = 0 (see [11]). It can be seen that the standard deviation of ˆ SA (j!k ) (dark the classical indirect spectral analysis estimate G gray line) is 30 dB above that of the estimated BLTI approximation (black dashed line). It follows that the input-output ˜TV (k) and Y˜TV (k) is the dominant error term time-variation U ˆ SA (j!k ). in G
Figure 13. Time-variant electronic circuit with input u (t), output y (t), and scheduling parameter p (t). Schematic (a): high gain second order bandpass filter with time-varying resonance frequency and damping ratio. It consists of the following components: three high gain operational amplifiers (TL071), three time-varying resistors R (p(t)), four resistors (R1 = R2 = 10 k⌦, R3 = 5.31 k⌦, and R4 = 100.8 k⌦), and two capacitors (C1 = C2 = 10 nF). Schematic (b): practical realization of the time-varying resistor using an operational amplifier, a 330 ⌦ resistor, and an electro-optical component (VTL5C1) consisting of a light dependent resistor (LDR) and a light emitting diode (LED).
Reference DFT spectrum Amplitude (dB)
0
50
100 0
10 20 30 Frequency (kHz)
Input DFT spectrum
Output DFT spectrum
r (t) =
F X
Ak sin (2⇡kf0 t +
k)
(39)
k=1
with uniformly distributed phases k , and equal amplitudes in the frequency band [228.9 Hz, 39.98 kHz] chosen such that the rms value is 50 mV (f0 = fs /P , fs = 625 kHz, P = 16384, A1 = A2 = · · · = A5 = 0, and A6 = A7 = · · · = AF ). A periodic ramp p (t) is applied to the electro-optical components in Figure 13 (period 1/f0 and voltage increasing linearly from 0.8412 V to 0.9317 V). N = 15P/16 = 15360 inputoutput samples of the steady state response are acquired using a band-limited measurement setup (u (t) and y (t) are lowpass filtered before sampling). Figure 15 shows the corresponding reference and input-output DFT spectra (5). During the acquisition time the voltage p (t) increased linearly from 0.8412 V to 0.9259 V. B. Results Starting from N known reference, noisy input-output samples r (nTs ), u (nTs ) and y (nTs ), the multiple-input, twooutput (MITO) model (22), with fp (t) Legendre polynomials satisfying (11), is estimated following the procedure of Section III-C (local polynomial method [16] with a fourth order local polynomial approximation of the FRFs and the leakage
Amplitude (dB)
100 0
40
Hru,1(jωk)
10 20 30 Frequency (kHz)
40
Hry,1(jωk) −20 Amplitude (dB)
−20 −40 −60 −80 0
10 20 30 Frequency (kHz)
−40 −60 −80 0
40
Hru,2(jωk)
40
−20
−40 −60 −80 0
10 20 30 Frequency (kHz) Hry,2(jωk)
−20 Amplitude (dB)
The device under test is a time-variant electronic circuit (see Fig. 13) that behaves as a high gain second order bandpass filter with time-varying resonance frequency and damping ratio. It operates in feedback as shown in Figure 14, and the time-variation is controlled via the voltage p (t) (scheduling parameter) of the electro-optical components (see Fig. 13, bottom schematic). The reference signal r (t) is a random phase multisine consisting of the sum of F = 1043 harmonically related sinewaves
10 20 30 Frequency (kHz)
50
Figure 15. DFT spectra (black) of the known reference (top), noisy input (bottom left), and noisy output (bottom right) signals – time-variant electronic ˆ ˜TV (k) and circuit. Light gray lines: estimated input-output time-variation U ˆ Y˜TV (k) (27); top dark gray lines: estimated input-output leakage errors TˆU (j!k ) and TˆY (j!k ) (22); and bottom medium gray lines: estimated input-output noise variances.
Amplitude (dB)
A. Measurement Setup
0
50
100 0
Amplitude (dB)
VI. M EASUREMENT E XAMPLE
Amplitude (dB)
0
Figure 14. Time-variant electronic circuit (TV circuit) of Figure 13 operating in feedback. r (t), u (t), and y (t) are, respectively, the reference signal stored in the arbitrary waveform generator, the input of the time-variant circuit, and the output of the time-variant circuit. A 50 ⌦ resistor is added at the output to match the 50 ⌦ input impedance of the acquisition channel. Rg is the 50 ⌦ output impedance of the arbitrary waveform generator.
40
10 20 30 Frequency (kHz)
40
−40 −60 −80 0
10 20 30 Frequency (kHz)
40
Figure 16. Frequency response functions Hrz,p (j!k ), p = 1, 2, of the MITO equivalent model (22) – time-variant electronic circuit. Black: input ˆ ru,p (j!k ) and H ˆ ry,p (j!k ); and light (left) and output (right) estimates H ˆ rz,p (j!k ). gray: estimated variance of H
term, and six degrees of freedom for the variance estimates). Figure 16 shows the estimated dynamics Hru,p (j!) and Hry,p (j!) of the time-varying branches for Nb = 2. It can be seen that the time-variation is significant over the whole frequency band. Increasing Nb further makes no sense ˆ rx,3 (j!k ) | ⇠ std(H ˆ rx,3 (j!k )), because for all values of k, |H x = u, y. Figure 15 shows the corresponding estimates of the inputoutput leakage errors TU (j!k ) and TY (j!k ), and the inputoutput noise variances. It can be seen that the output signal-tonoise ratio (SNR) is about 60 dB, while the input SNR varies between 40 dB and 50 dB. Note that the leakage errors lie well above the noise floor. Figure 17 shows the nonparametric estimate of the best lin-
G
(j
GSA(jωk)
Output residual
40
20
10 20 30 Frequency (kHz)
40
Amplitude (dB)
40 20 0 20 40 60 80 0
)
k
Amplitude (dB)
Amplitude (dB)
BLTI
40 60 80 100 0
10 20 30 Frequency (kHz)
40
Figure 17. Estimated best linear time invariant (BLTI) approximation (left) and its output residual (right) – time-variant electronic circuit. Left plot: the ˆ BLTI (j!k ) (black), its variance var(G ˆ BLTI (j!k )) (dark gray), and estimate G ˆ ˆ ˆ its bias |GBLTI (j!k ) GBLTI (j!k ) | (light gray). Right plot: comparison ˜TV (k) (27) of the time-varying part Y˜TV (k) of the between the estimate Yˆ output DFT spectrum (black), and the estimated output residual YˆTV (k) (28) of the BLTI approximation (light gray).
ear time-invariant (BLTI) approximation obtained via Eqs. (25) ˆ ˆ BLTI (j!k ) G ˆ BLTI (j!k ) is and (30). Note that the difference G ˆ an estimate of the bias (26) of GBLTI (j!k ). From Fig. 17 it follows that this bias is smaller than the standard deviˆ BLTI (j!k ) and, hence, can be neglected. Using ation of G ˆ GBLTI (j!k ), the time-varying parts of the input-output DFT spectra Z˜TV (k) are calculated via (27). From Figure 15 it can be seen that the input and output time-variations are well above the noise level and are of the same order of magnitude as, respectively, the input and output leakage errors. Note ˜TV (k) (light also that the observed output time-variation Yˆ gray line) is about 40 dB below the signal level. From this number one could wrongly conclude that the electronic circuit is slowly time-varying. Indeed, the true output time-variation YTV (k), estimated via (28), is much larger than Y˜TV (k) at those frequencies where |GBLTI (j!k) | 1 (see Fig. 17). From Figures 15 and 17 it can be seen that YTV (k) is about 20 dB below Y (k), indicating a large (10%) time-variation of the electronic circuit over the experiment duration. Figure 18 compares the indirect spectral analysis estimate (9), obtained by applying the procedure of Section III-C with Nb = 0 (see [11]), with the estimated BLTI approximation (procedure of Section III-C with Nb = 2). The difference between the variances of both estimates is maximal in the passband of the system. This is exactly the frequency band where the time-variation is the largest. It nicely illustrates the benefit of the proposed procedure: it eliminates the random influence of the time-variation on the FRF measurement. VII. C ONCLUSIONS The main difficulty of estimating the best linear timeinvariant (BLTI) approximation of a linear time-variant (LTV) system operating in closed loop is that both the output and the input signals are disturbed by the time-variation. Indeed, the input time-variation causes a bias in the estimation procedure designed for open loop systems [1]. To solve this problem the concept of the best linear time-invariant (BLTI) approximation of a linear time-variant (LTV) system has been generalized. This generalized definition can also handle time-variant actuator and feedback dynamics. For open loop
20 0 −20 −40 0
10 20 30 Frequency (kHz)
40
Figure 18. Spectral analysis estimate GSA (j!k ) (9) of the best linear time-invariant approximation – time-variant electronic circuit. Black: estimate ˆ SA (j!k ); dark gray: variance estimate var(G ˆ SA (j!k )); light gray: differG ˆ SA (j!k ) G ˆ BLTI (j!k ) |; and black dashed line: var(G ˆ BLTI (j!k )). ence |G
LTV systems driven by linear time-invariant actuators, the generalized definition simplifies to the classical definition [1]. The difference yTV (t) between the actual output of the LTV system and the output of the BLTI approximation quantifies the amount of time-variation of the system. It allows users to decide whether the time-invariant framework (the BLTI approximation) is accurate enough for their particular application. It has been shown in this paper that yTV (t) is equal to the output time-variation minus the input time-variation filtered by the BLTI approximation. This linear correction remains exact for time-variant actuators and feedback dynamics. A nonparametric procedure for estimating the BLTI approximation and the output residual yTV (t) from known reference and noisy input-output samples has been proposed. It is based on a multiple-input, two-output linear time-invariant equivalent model of the open loop linear time-variant system from reference to input-output. In those frequency bands where the time-variation is the dominant disturbance, the variance of the estimated BLTI approximation is much smaller than that of the classical indirect spectral analysis method. If the linear system operating in closed loop is timeinvariant, and the actuator and/or feedback dynamics are timevariant, then the BLTI approximation equals the frequency response function (FRF) within an O(T 1 ) bias error, where T is the experiment duration. However, using the multipleinput, two-output linear time-invariant equivalent model, the time-invariant behaviour can be detected and an unbiased estimate of the FRF can be obtained. ACKNOWLEDGEMENT This work is sponsored by the Research Council of the Vrije Universiteit Brussel, the Research Foundation Flanders (FWO-Vlaanderen), the Flemish Government (Methusalem Fund METH1), and the Belgian Federal Government (Interuniversity Attraction Poles programme IAP VII DYSCO). A PPENDIX A. Proof of Eqs. (9) and (10) The key idea consists in modeling the closed loop system in Fig. 4 from the reference signal r (t) to the input-output signals z (t) = [y (t) u (t)]T simultaneously. As such, the single-input, single-output closed loop problem is transformed into a single-input, two-output open loop problem, for which
Laplace transform of u (t), and U (s) = Gact (s) R (s), it follows that the time-varying transfer function from reference to output is given by Gry (s, t) = G (s, t) Gact (s) . Figure 19. Equivalent one-input, two-output open loop representation of a linear time-variant system operating in closed loop (see Fig. 4).
the results of [1], [11] apply. Using the results of Section II-A we get Z (k) = Grz,BLTI (j!k ) R (k) + Z˜TV (k)
(40)
with Grz,BLTI = [Gry,BLTI Gru,BLTI ]T the 2 ⇥ 1 BLTI approximation of the single-input, two-output LTV, and where ˜TV (k)]T is uncorrelated with – but Z˜TV (k) = [ Y˜TV (k) U not independent of – the reference R (k). Using (40) we find E Y (k) R (k) E U (k) R (k)
=
Gry,BLTI (j!k ) Gru,BLTI (j!k )
(41)
Since Grx,BLTI (j!k ), x = y, u, is – by definition – the BLTI approximation of the open loop LTV system with time-varying transfer function Grx (s, t), we have that ˆ 1 T Grx,BLTI (j!k ) = Grx (j!k , t) dt (42) T 0 (see Eq. (2)). Combining (7), (41), and (42) proves (9). Calculating the difference YTV (k) between the actual output of the closed loop LTV system and the output of the BLTI approximation (7) gives YTV (k)
=
Y (k) GBLTI (j!k ) U (k) ˜TV (k) Y˜TV (k) GBLTI (j!k ) U
=
(43)
where the second equality uses Eq. (40). It shows that ˜TV (k) are YTV (k) is uncorrelated with R (k) (Y˜TV (k) and U both uncorrelated with R (k)), which concludes the proof of Eq. (10). B. Proof of the Properties of the BLTI Approximation (7) In this appendix we first prove that the generalized definition (7) reduces to (2) for open loop systems. Next, we show that the BLTI approximation (7) is independent of the actuator characteristics. Finally, the relationship (12) between the BLTI approximation (7) and the series expansion (8) is proven. If the LTV system operates in open loop (see Fig. 3, top block diagram), then the time-varying transfer function from reference to input is equal to the time-invariant actuator transfer function Gru (s, t) = Gact (s) .
(44)
Using y (t) = L {G (s, t) U (s)} (see [23], [24]), with L 1 {} the inverse Laplace transform operator and U (s) the 1
(45)
Combining Eqs. (7), (44), and (45) gives (2). Figure 19 gives an equivalent one-input, two-output open loop representation of the closed loop LTV system in Figure 4 (top block diagram). Using x (t) = L 1 {Grx (s, t) X (s)} (see [23], [24]) with x = y, u, X = Y, U , and X (s) = Gact (s) R (s) it follows that Gry (s, t)
=
Gwy (s, t) Gact (s)
(46)
Gru (s, t)
=
Gwu (s, t) Gact (s)
(47)
where w (t) is the output of the actuator (see Figures 4 and 19). Combining Eqs. (7), (46), and (47) shows that GBLTI (s) is independent of the actuator characteristics, viz. ´ 1 T T 0 Gwy (s, t) dt GBLTI (s) = ´ T . (48) 1 T 0 Gwu (s, t) dt
In general, the BLTI approximation (48) depends on the feedback transfer function Gfb (s) because, in general, the timevarying transfer function Guq (s, t) from the input u (t) of the LTV system to the output q (t) of the feedback branch (see Figure 4) is different from G (s, t) Gfb (s) (see [23], [24]). Integrating (8) over [0, T ] using (11) gives ˆ 1 T Grz (s, t) dt = Grz,0 (s) (49) T 0 Combining Eqs. (7) and (49) proves (12). C. Proof of the Bias (26) of the BLTI Approximation Estimate ˆ rz,0 (j!k ) (23) can be written The nonparametric estimate G as (to simplify the notations we drop the arguments) ˆ rz,0 = Grz,0 + nG,z G with Grz,0 the true value, and where E{nG,z } = O(T Using (50), Eq. (25) becomes ˆ BLTI = Gry,0 + nG,y = GBLTI 1 + nG,y /Gry,0 G Gru,0 + nG,u 1 + nG,u /Gru,0
(50) 2
).
(51)
For T sufficiently large, one can always choose the order of the polynomial approximation and the width of the local frequency band in the LPM method [16], [17] such that nG,y nG,u ,
Abstract—Recently a method has been developed for detecting and quantifying the time-variation in frequency response function (FRF) measurements using arbitrary excitations [1]. The following basic assumptions have been made: (i) the input is known exactly (generalized output error stochastic framework), and (ii) the time-variant system operates in open loop. The latter excludes any interaction between the time-variant system and the generator/actuator. In this paper we extend the results of [1] to noisy input-output observations (errors-in-variables stochastic framework) of time-variant systems operating in closed loop. Index Terms—frequency response function, best linear timeinvariant approximation, slowly time-varying, feedback, errorsin-variables, continuous-time, nonparametric modeling
I. I NTRODUCTION In all kinds of engineering applications frequency response function (FRF) measurements are a basic tool to obtain quickly insight into the dynamic behavior of complex systems [2], [3], [4]. Examples of complex systems that are studied via FRF measurements are (i) corrosion of metals where the time-variation is induced by the changing surface condition [5], (ii) airplane dynamics during take-off and landing which depend on the changing flight speed and height [6], and (iii) myocardium tissue impedance measurements that evolve during the cardiac cycle [7]. In practise, these real-life systems satisfy only approximately the linearity and time-invariance assumptions inherent to the concept of an FRF. While the influence of nonlinear distortions on FRF measurements using random excitations has been studied in detail in the literature [4], [8], [9], [10], much less is known about the impact of time-variations. Recently it has been shown that time-variations cause errors on FRF measurements that are correlated over the frequency [11]. Since the leakage errors have exactly the same property [4], [12], it is difficult to distinguish and quantify these two sources of errors [11]. Assuming that the input is known exactly and that the linear time-variant system operates in open loop, this problem has been solved in [1] by modeling nonparametrically the timevariation using multiple-input, single-output FRF estimates. The contribution of this paper is to extend the results of [1] to noisy input-output observations of time-variant systems operating in closed loop (see Figure 1).
Figure 1. Linear time-variant (LTV) system operating in feedback. The actuator and feedback dynamics are linear and time-invariant (LTI). r (t), w (t), u (t), and y (t) are, respectively, the reference signal, the output of the actuator, and the input and output of the LTV system.
Feedback in a measurement setup is due to either an explicit control action (see Fig. 1) or an interaction between an nonideal generator/actuator and the device under test. The latter is illustrated by the problem of measuring a time-varying impedance Z (s, t) [13] from noisy current i (t) and voltage v (t) observations (see Fig. 2). The relationship between the current through and the voltage across the impedance Ze (s) in Figure 2 is given by I(s) = (E(s)
V (s)) /Ze (s)
(1)
where I(s) = L{i(t)}, V (s) = L{v(t)}, and E(s) = L{e(t)}, with L{} the Laplace transform. Formula (1), with E(s) = R(s) (e (t) = r (t)), I(s) = U (s) (i (t) = u (t)), and V (s) = Y (s) (v (t) = y (t)), is exactly the feedback law in Figure 1, where the actuator and feedback dynamics are both equal to 1/Ze (s). Note that the feedback in the setup of Figure 2 is solely due to the non-zero output impedance of the voltage source. Hence, feedback is present in any measurement setup with a non-ideal generator/actuator. For example, a signal generator operating above 10 MHz has typically a 50 ⌦ output impedance, and a force actuator (shaker) for mechanical impedance measurements has a frequency dependent output impedance. It emphasizes the importance of handling timevariant systems operating in closed loop. The outline of the paper is as follows. First, the definition of the best linear time-invariant (BLTI) approximation of a linear time-variant system operating in feedback is given, and its properties are discussed (Section II). Next, it is shown in Section III how the BLTI approximation can be obtained from noisy input-output measurements using multiple-input, multiple-output nonparametric FRF estimators. Further, the
This is the author’s version of the article published as R. Pintelon, E. Louarroudi, and J. Lataire. Detection and quantification of the influence of time-variation in closed loop frequency response function measurements. IEEE Transactions on Instrumentation & Measurement, 62(4):853–863, 2013.
G (s, t) =
1 X
Gp (s) fp (t)
p=0
Figure 2. Measurement of a time-varying impedance Z (s, t) using a voltage source e (t) with output impedance Ze (s). v (t) and i (t) are the voltage across and the current through Z (s, t) respectively.
• •
t 2 (0, T )
(3)
with fp (t) a complete set of basis functions over the interval (0, T ). The excitation u (t) is either periodic noise or a random phase multisine with period length T . The LTI systems Gp (s) in (3) operate in steady state.
Under Assumption 1 it is shown in [11] that the direct frequency response function (FRF) estimate obtained via crossand autopower spectra equals the BLTI approximation (2) GBLTI (j!k ) = Figure 3. Best linear time-invariant (BLTI) approximation GBLTI (s) (bottom block diagram – frequency domain) of a linear time-variant (LTV) system G (s, t) operating in open loop (top block diagram – time domain). Gact (s) is the transfer function of the linear time-invariant actuator. r (t), u (t), and y (t) are, respectively, the reference signal, and the input and output of the LTV system. YTV (k) is the DFT spectrum of the difference between the actual output of the LTV system and the output of the BLTI approximation. It is uncorrelated with - but not independent of - the DFT spectrum U (k) of the input u (t).
case of time-variant actuator and/or feedback dynamics is handled in Section IV. The theory is illustrated on a simulation (Section V) and a measurement (Section VI) example. Finally, some conclusions are drawn in Section VII. II. D EFINITION OF THE B EST L INEAR T IME -I NVARIANT A PPROXIMATION First, the difficulty of estimating the best linear timeinvariant (BLTI) approximation of a linear time-variant system operating in closed loop using the open loop definition introduced in [1], [11] is revealed. Next, a generalized definition that can handle the closed loop case is given, and its properties are discussed. A. Open Loop Case If the system operates in open loop (see Fig. 3, top block diagram), then the BLTI approximation of the linear timevariant system over the time interval [0, T ] is defined as [1] 1 GBLTI (s) = T
ˆ
T
G (s, t) dt
(2)
0
where G (s, t) = L {g (t, t ⌧ )} is the time-varying transfer function, with L {} the Laplace transform w.r.t. ⌧ , and g (t, ⌧ ) the time-varying impulse response (i.e. response of the system to a Dirac impulse at time ⌧ ). Assumption 1. (Open loop case) •
The time-varying transfer function G (s, t) can be expanded in series as
E Y (k) U (k) n o 2 E |U (k)|
(4)
where the expected value E {} is taken w.r.t. the random realization of u (t), and where X (k), X = Y or U , is the discrete Fourier transform (DFT) of x (nTs ), x = y or u, N 1 1 X X (k) = p x (nTs ) e N n=0
j2⇡kn/N
(5)
with Ts the sample period and T = N Ts . Hence, the actual output Y (k) of the LTV system can be written as (see also Fig. 3, bottom block diagram) Y (k) = GBLTI (j!k ) U (k) + YTV (k)
(6)
where the residual YTV (k) is uncorrelated with – but not independent of – the input U (k). It shows that under Assumption 1, GBLTI (s) (2) is the best (in least squares sense) LTI approximation of the LTV system. If definition (2) is applied to the closed loop case (see Fig. 4), then we obtain biased estimates of the BLTI approximation. Indeed, combining (4) and (6) gives E Y (k) U (k) E YTV (k) U (k) n o = GBLTI (j!k ) + n o 2 2 E |U (k)| E |U (k)|
where E YTV (k) U (k) 6= 0 due to the feedback loop. This is the key difficulty of estimating the BLTI approximation of an LTV system operating in closed loop. In the next section we present a solution to this problem. B. Closed Loop Case Consider a linear time-variant (LTV) system operating in feedback (see Fig. 4, top block diagram). Inspired by the indirect method for measuring the FRF of an LTI system operating in closed loop [19], [20], we propose the following generalized definition for the BLTI approximation GBLTI (s) =
1 T 1 T
´T 0
´T 0
Gry (s, t) dt Gru (s, t) dt
(7)
where the residual YTV (k) is uncorrelated with – but not independent of – the reference R (k) (proof: see Appendix A). It shows that under Assumption 2, the best linear time invariant approximation GBLTI (s) (7) of the time-variant system operating in closed loop, is the ratio of the best (in least squares sense) LTI approximations of two time-variant open loop systems: that from reference to output to that from reference to input. Properties of the BLTI Approximation • •
Figure 4. Best linear time-invariant (BLTI) approximation GBLTI (s) (bottom block diagram – frequency domain) of a linear time-variant (LTV) system G (s, t) operating in feedback (top block diagram – time domain). Gact (s) and Gfb (s) represent the linear time-invariant actuator and feedback dynamics, respectively. r (t), w (t), u (t), y (t), and q (t) are, respectively, the reference, the output of the actuator, the input and output of the LTV system, and the output of the feedback branch. YTV (k) is the DFT spectrum of the difference between the actual output of the LTV system and the output of the BLTI approximation. It is uncorrelated with - but not independent of the DFT spectrum R (k) of the reference signal r (t).
where Gru (s, t) and Gry (s, t) are the time-varying transfer functions from the reference r (t) to, respectively, the input u (t) and the output y (t). The 2 ⇥ 1 time-varying transfer function [Gry (s, t) Gru (s, t)]T , with Grz (s, t) = z = [y u]T , can be expanded in series as Grz (s, t) =
1 X
Grz,p (s) fp (t)
p=0
• •
t 2 (0, T )
(8)
with fp (t) a complete set of basis functions over the interval (0, T ). The reference signal r (t) is either periodic noise or a random phase multisine with period length T . The 2 ⇥ 1 LTI systems Grz,p (s) in (8) operate in steady state.
Under Assumption 2 it is shown in Appendix A that the indirect frequency response function (FRF) estimate obtained via crosspower spectra equals the BLTI approximation (7) GBLTI (j!k ) =
f0 (t) = 1 and
1 T
ˆ
T
fp (t) dt = 0 for p > 0. (11)
0
With this choice it follows that GBLTI (s) =
Gry,0 (s) Gru,0 (s)
(12)
(proof: see Appendix B).
Assumption 2. (Closed loop case) •
•
For LTV systems operating in open loop, the generalized definition (7) reduces to (2) (proof: see Appendix B). The BLTI approximation (7) depends on the feedback dynamics Gfb (s), but is independent of the actuator characteristics Gact (s) (proof: see Appendix B). For example, the BLTI approximation of the time-varying impedance in Figure 2 depends on the value of the resistor R. Without any loss in generality, the basis functions fp (t) in (8) can always be chosen such that
E Y (k) R (k) E U (k) R (k)
Applying Mason’s rule (see [21]) to the bottom block diagram of Figure 4, it can easily be verified that, under Assumption 2, the input-output DFT spectra can be written as ˜TV (k) U (k) = UR (k) + U
(13)
Y (k) = GBLTI (j!k ) UR (k) + Y˜TV (k)
(14)
with UR (k) that part of the input that is correlated with the reference R (k), UR (k) =
(10)
Gact (j!k ) R (k) 1 + GBLTI (j!k ) Gfb (j!k )
˜TV (k), Y˜TV (k) the time-varying input-output contribuand U tions that are uncorrelated with R (k)
(9)
where the expected value E {} is taken w.r.t. the random realization of r (t). Moreover, the actual output Y (k) of the LTV system can be written as (see also Fig. 4, bottom block diagram) Y (k) = GBLTI (j!k ) U (k) + YTV (k)
Interpretation of the Output Residual
˜TV (k) U
=
Y˜TV (k)
=
Gfb (j!k ) YTV (k) (15) 1 + GBLTI (j!k ) Gfb (j!k ) 1 YTV (k) (16) 1 + GBLTI (j!k ) Gfb (j!k )
Since the input U (k) (13) is also subject to time-variation ˜TV (k), the output time-variation Y˜TV (k) in (14) cannot U be used to quantify the time-variation of the closed loop
z (t) = [ y (t) u (t) ]T . Proceeding in this way the singleinput, single-output closed loop problem is transformed into a single-input, two-output open loop problem for which the results of [1] apply. First, we show that there exists a multipleinput, two-output linear time-invariant equivalent model of the single-input r (t), two-output z (t) slowly time-varying system. Next, the relationship between the DFT spectra of the reference and the input-output signals is established. Finally, an algorithm for estimating the BLTI approximation of an LTV system operating in feedback is presented. A. Linear Time-Invariant Equivalence of a Linear Slowly Time-Varying System A linear time-variant system operating in closed loop is called slowly-time varying if the time-varying transfer function Grz (s, t) of the open loop system from reference r to inputoutput z = [y u]T is slowly time-varying: Definition 3. A linear time-variant system operating in closed loop is slowly time-varying over the interval [0, T ] if the time-varying transfer function Grz (s, t) (8) can be written as Grz (s, t) =
Nb X
Grz,p (s) fp (t)
p=0
(18)
t 2 (0, T )
with fp (t) a complete set of basis functions over the interval (0, T ) satisfying (11).
Figure 5. Linear time-invariant (LTI) equivalent open loop model (bottom) of a linear slowly time-varying system operating in feedback (top). Middle block diagram: single-input r (t), two-output z (t) = [ y (t) u (t) ]T linear time-variant open loop equivalent representation of a single-input u (t), single-output y (t) linear slowly time-varying system operating in feedback (top block diagram). Bottom block diagram: multiple-input r (t), r1 (t), . . . , rNb (t), two-output z (t) = [ y (t) u (t) ]T linear time-invariant equivalent model (bottom block diagram) of a single-input r (t), two-output z (t) = [ y (t) u (t) ]T linear time-variant open loop system (middle block diagram).
LTV system. However, multiplying (15) with GBLTI (j!k ) and subtracting the result from (16) gives exactly YTV (k) Y˜TV (k)
˜TV (k) = YTV (k) GBLTI (j!k ) U
(17)
which is consistent with Eq. (43) of Appendix A. It shows that the output residual YTV (k) of the BLTI approximation is the time-varying part of the output Y˜TV (k) minus the input time˜TV (k) filtered through the BLTI approximation variation U GBLTI (j!k ). III. E STIMATING THE BLTI A PPROXIMATION OF S LOWLY T IME -VARYING S YSTEMS O PERATING IN F EEDBACK The basic idea consists in modeling the linear time-variant (LTV) system from the reference r (t) to the input and output
Figure 5 shows a block diagram of the slowly timevarying system (18). If the time-variant system operates in open loop (Gfb (s) = 0 in Fig. 5), then Grz (s, t) = Gact (s) [ G (s, t) 1 ]T (see Appendix B, Eqs. (44) and (45)), and Definition 3 reduces to that for open loop systems [1]. Using Theorem 2 of [1], it follows that the output z (t) = [ y (t) u (t) ]T of the single-input r (t) open loop slowly time-varying system (18) can be calculated via an equivalent multiple-input rp (t) = r (t) fp (t), p = 0, 1, . . . , Nb , two-output z (t) = [ y (t) u (t) ]T open loop linear timeinvariant model (see also Fig. 5)
z (t)
=
L
1
(N b X p=0
=
L
1
(N b X p=0
Grz,p (s) L {r (t)} fp (t) )
Hrz,p (s) L {rp (t)}
)
(19) (20)
with L 1 {} and L{}, respectively, the inverse Laplace transform and the Laplace transform operators. Choosing Legendre polynomials Pp (2t/T 1), t 2 [0, T ], satisfying (11) as basis functions fp (t) (see [14]), the relationship between Grz,0 (s), the first term in the sum (19), and the transfer functions Hrz,p (s), p = 0, 1, . . . , Nb , of the equivalent LTI model (20), is given by
Grz,0 (s) = Hrz,0 (s) + j
Nb 2
k
4 X T 2 p=0
2 T
j
(2) 2p Hrz,2p
Nb 1 2
X
k
(1)
Hrz,2p+1 (s) + . . .
p=0
(s) + O T
3
(21)
2
with 2p = 1.5+2.5 (p 1)+(p 1) , bxc the largest integer (m) smaller than or equal to x, Hrz,p (s) the mth derivative of Hrz,p (s) w.r.t. s, and where O(T m ) with m > 0 means that limT !1 T m O(T m ) < 1 (proof: apply Theorem 3 of [1] to Eqs. (20) and (19)). The bias term O(T 3 ) depends (m) on the higher order derivatives Hrz,p (s), with m, p > 3 and m 6 p. Finally, combining Eqs. (12) and (21) establishes the relationship between the BLTI approximation (7) of the linear slowly time-varying system operating in closed loop and the transfer functions Hrz,p (s), p = 0, 1, . . . , Nb , of the open loop LTI equivalent model. It is the basis of the nonparametric estimation procedure developed in Section III-C. B. Relationship Between the DFT Spectra The relationship between the sampled known reference signal and the sampled noiseless input-output signals is established for band-limited reference signals . Assumption 4. The reference signal r (t) is band-limited: the power spectral density of the reference Srr (j!) is zero for f > fB , where fB < fs /2, with fs = 1/Ts the sampling frequency. Under Assumption 4, the reference R (k) and input-output Z (k) = [ Y (k) U (k) ]T DFT spectra (5) of a slowly timevarying system operating in feedback (see Definition 3) are related by Z (k) =
Nb X
Hrz,p (j!k ) Rp (k) + TZ (j!k )
(22)
p=0
with Rp (k) = DFT{r (t) fp (t)}, and where TZ (j!k ) models the sum of the spectral leakage and the residual alias errors (proof: apply Theorem 5 of [1] to Eq. (20)). C. Estimation Procedure In practise the nonparametric estimate of the best linear time-invariant (BLTI) approximation is calculated from known reference r (nTs ) and noisy input-output u (nTs ), y (nTs ) samples, n = 0, 1, . . . , N 1 (see Fig. 6). It is assumed that the input-output measurement errors nu (t), ny (t) in Figure 6 can be modelled as stationary (mutually correlated) zero mean, filtered (band-limited) white noise [4], [15]. The reference signal r (t) satisfies Assumption 4, and the input-output signals u (t) and y (t) are lowpass filtered before sampling.
Figure 6. Known reference r (t), noisy input-output u (t), y (t) observations of a linear time-variant system operating in feedback. nu (t) and ny (t) are the (mutually correlated) input-output measurement errors, and u0 (t), y0 (t) are the unknown true input-output signals.
The procedure for estimating the BLTI approximation GBLTI (j!k ) and the output residual YTV (k) consists of the following four steps: 1) Nonparametric estimation of Hrz,p (j!k ): From Eq. (22), it follows that the Nb +1 FRFs Hrz,p (j!k ) can be estimated using nonparametric multiple-input, multipleoutput (MIMO) methods. The local polynomial method ([1], [16], [17]) is used here because of its excellent leakage (TZ (j!k ) in Eq. (22)) suppression property, and because it estimates the FRFs at the full frequency resolution 1/T of the experiment duration T . 2) Estimation of Nb : The number Nb of parallel branches modeling the time-variation in (22) is determined as follows. Starting with Nb = 1, the number Nb is ˆ rz,N (j!k ) is increased till the estimated 2 ⇥ 1 FRF H b zero within its uncertainty (the LPM method calculates ˆ rx,N (j!k ) | ⇠ also the covariance of the estimates): |H b ˆ std(Hrx,Nb (j!k )) for x = y, u. The number of significant parallel branches is then Nb 1. 3) Nonparametric estimation of Grz,0 (j!k ): First, a nonparametric estimate of the 2 ⇥ 1 FRF Grz,0 (j!k ) is obtained via (21) ˆ rz,0 (j!k ) = H ˆ rz,0 (j!k ) + 2 G T
j
Nb 1 2
X
k
ˆ (1) H rz,2p+1 (j!k )
p=0
(23) where the first order derivatives are replaced by central differences ˆ rz,m (j!k+1 ) H (1) ˆ rz,m H (j!k ) = j!k+1
ˆ rz,m (j!k H j!k 1
1)
(24) Next, using (12), an estimate of the BLTI approximation is calculated as ˆ ˆ BLTI (j!k ) = Gry,0 (j!k ) G ˆ ru,0 (j!k ) G
(25)
Since the expected value of (23) equals the true value within an O T 2 bias error if at least a second order polynomial approximation is used in the LPM
ˆ BLTI (j!) (25) is method (proof: see [1]), the bias of G asymptotically (T ! 1) given by n o ˆ BLTI (j!k ) E G GBLTI (j!k ) = O T 2 (26)
(proof: see Appendix C). 4) Nonparametric estimation of YTV (k): First, the output ˜TV (k) ]T (see residual Z˜TV (k) = [ Y˜TV (k) U Appendix A, Eq. (40)) of the BLTI approximation from reference r (t) to input-output z (t) = [ y (t) u (t) ]T is estimated as
Figure 7. Linear time-invariant equivalent open loop model of a linear timeinvariant system with time-variant actuator and/or feedback dynamics.
Nb
X ˜TV (k) = ˆ rz,p (j!k ) Rp (k) Zˆ H ⇣
p=1
ˆ rz,0 (j!k ) G
···
⌘
ˆ rz,0 (j!k ) R (k) H
A. Time-Variant System (27)
where the second term in the difference has an O(T 1 ) contribution (see [1], Eq. (18)). Next, using (17), we obtain an estimate for the output residual YTV (k) of the LTV system operating in closed loop ˜TV (k) YˆTV (k) = Yˆ
ˆ ˆ BLTI (j!k ) U ˜TV (k) G
(28)
The bias of (28) is an O(T ) because the bias of the ˜TV (k) and G ˆ BLTI (j!k ) are both an O(T 2 ). estimates Zˆ 2
Notes •
Using (21), the O(T 2 ) bias (26) of the BLTI approximation estimate (25) can be reduced to an O(T 3 ). First, an improved estimate of Grz,0 (j!k ) is calculated as ˆ ˆ rz,0 (j!k ) = G ˆ rz,0 (j!k ) + · · · G j
Nb 2
k
4 X T 2 p=1
ˆ (2) 2p Hrz,2p
(j!k )
(29)
where the second order derivatives are replaced by second order central differences (see [1], [18]). Next, combining (12) and (29) gives an improved BLTI approximation estimate ˆ ˆ ˆ ˆ BLTI (j!k ) = Gry,0 (j!k ) G ˆ ˆ ru,0 (j!k ) G •
(30)
ˆ ˆ BLTI (j!k ) (25) by G ˆ BLTI (j!k ) (30) in Eq. Replacing G 2 (28) reduces the O(T ) bias of YˆTV (k) to an O(T 3 ). IV. T IME -VARIANT ACTUATOR AND / OR F EEDBACK DYNAMICS
In this section the results are extended to time-variant actuator and/or feedback dynamics. Two cases are distinghuised: the system is time-variant or time-invariant.
Eqs. (9) and (10) remain valid for time-variant actuator and/or feedback dynamics because the proof in Appendix A does not use the time-invariance of the actuator and the feedback loop. It shows that Figure 4 remains valid for timevariant actuator Gact (s, t) and feedback Gfb (s, t) transfer functions. However, the BLTI approximation (7) depends now on the time-variant actuator dynamics; even in open loop. B. Time-Invariant System If the system is time-invariant, G (s, t) = G (s), then one would expect that the BLTI approximation (7) equals G (s). This is not true because Gry (s, t) 6= Gru (s, t) G (s) [22]. However, in the sequel of this section we will show that the difference between GBLTI (s) and G (s) is an O(T 1 ). Define the linear time-invariant (LTI) approximation of the linear time-variant (LTV) system G (s, t) as GLTI (s) =
Hry,0 (s) Hru,0 (s)
(31)
where Hrz,0 (s) is the top branch in the LTI equivalent model of the LTV system (see Fig. 5). The output contribution of the time-varying branches in the LTI equivalent model is given by Z˘TV (k) =
Nb X
Hrz,p (j!k ) Rp (k)
(32)
p=1
Using (32), one finds the output residual YˇTV (k) of the LTI approximation (31) YˇTV (k) = Y˘TV (k)
˘TV (k) GLTI (j!k ) U
(33)
Applying the results of Figure 5 to the time-variant dynamics from reference r (t) to the input u (t), taking into account that the system is time-invariant G (s, t) = G (s), gives the linear time-invariant (LTI) equivalent open loop model shown in Figure 7. Shifting the LTI block G (s) in Figure 7 backwards in the parallel branches, transforms the block diagram in an LTI equivalent open loop model from reference to output. It shows that the dynamics Hry,p (s) are related to Hru,p (s) as Hry,p (s) = Hru,p (s) G (s)
(34)
Reference DFT spectrum
Poles (x) and zeroes (o) Amplitude (dB)
0.5 0.0 −0.5 −0.6
−0.4 −0.2 Real part
Combining (31)–(34), it can easily be verified that =
G (s)
(35)
=
0
(36)
if the system is time-invariant. Property (36) allows us to verify whether the system is time-invariant or not. If (36) is true within the noise standard deviation, then the system is timeinvariant and its frequency response function is calculated via Eq. (31). Otherwise, the system is time-variant and its BLTI approximation is estimated as explained in Section III-C. Since Grz,0 (s) = Hrz,0 (s) + O(T 1 ) (see Eq. (21)), it can be concluded from (12) and (31) that GBLTI (s) = GLTI (s) + O(T 1 ) (proof: follow the same lines of Appendix C). Hence, in the time-invariant case, the BLTI approximation equals the FRF within an O(T 1 ) bias error. V. S IMULATION E XAMPLE In this section we simulate a linear time-variant (LTV) system operating within a linear time-invariant (LTI) feedback loop and driven by an LTI actuator. Following the procedure of Section III-C, the best linear time-invariant (BLTI) approximation is estimated from noisy input-output observations of the LTV system (see the setup of Figure 6). The simulated LTV system is a third order differential equation with time-varying coefficients 3 X p=0
(p)
ap (t) y0 (t) =
3 X
(p)
bp (t) u0 (t)
0 10 20 0
0
Figure 8. Poles (⇥) and zeroes ( ) of the frozen transfer function (38) – simulation example. The real zero at 4 is fixed and not shown.
GLTI (s) YˇTV (k)
10
(37)
p=0
where x(p) (t) = dp x (t) /dtp , x = u0 , y0 . The poles and zeroes of the corresponding frozen transfer function P3 p p=0 bp (t) s Gf (s, t) = P3 (38) p p=0 ap (t) s
(transfer function obtained from (37) by fixing the time in ap (t) and bp (t)) are shown in Fig. 8 for t 2 [0, T ], with T = 3840 s. The system operates in closed loop (see Fig. 6) with Gfb (s) = 1, and the actuator transfer function Gact (s) is a sixth order Chebyshev filter with a passband ripple of 6 dB and a cut-off frequency of 0.3125 Hz. The reference signal r (t) is a zero mean band-limited random excitation with standard deviation equal to one satisfying Assumption 4 with fB = 0.2930 Hz (see Fig. 9, top plot). Using numerical integration (ODE45 solver in MATLAB®),
0.1 0.2 Frequency (Hz)
0.3
Input DFT spectrum
Output DFT spectrum
0
Amplitude (dB)
−1.0
Amplitude (dB)
Imaginary part
1.0
20 40 60 0
0.1 0.2 Frequency (Hz)
0.3
0 20 40 60 0
0.1 0.2 Frequency (Hz)
0.3
Figure 9. DFT spectra (black) of the known reference (top), noisy input (bottom left), and noisy output (bottom right) signals – simulation example. ˆ ˜TV (k) and Top light gray lines: estimated input-output time-variation U ˆ Y˜TV (k) (27); middle dark gray lines: estimated input-output leakage errors TˆU (j!k ) and TˆY (j!k ) (22); bottom medium gray lines: estimated inputoutput noise variances; and black dashed lines: true input-output noise variances.
the response y0 (t) is calculated for t 2 [0, T ] at the rate fs = 2 Hz, giving N = 7680 samples. The noiseless inputoutput samples u0 (t) and y0 (nTs ) are disturbed by white noise nu (nTs ) and ny (nTs ) (see Fig. 6), and the black dashed lines in the bottom plots of Fig. 9 show var (DFT(nx (nTs ))), x = u, y, as a function of the frequency. The multiple-input, two-output (MITO) model (22), with fp (t) Legendre polynomials satisfying (11), is estimated from the N known reference r (nTs ), noisy input u (nTs ), and noisy output y (nTs ) samples using the procedure of Section III-C (local polynomial method [16] with a sixth order local polynomial approximation of the FRFs and the leakage term, and six degrees of freedom for the variance estimates). Figure 10 shows the estimated dynamics of the time-varying branches ˆ rz,p (j!k ), p > 1, for Nb = 2. It can be seen that the H estimates coincide – within their uncertainty – with the true values. The latter are obtained by estimating (22) from the noiseless samples. Increasing Nb further makes no sense since ˆ rx,3 (j!k ) | ⇠ std(H ˆ rx,3 (j!k )) for x = u, y. |H Next, a nonparametric estimate of the best linear timeinvariant (BLTI) approximation of (37) is obtained via Eqs. (25) and (30). The results are shown in Figure 11, left plot. It can be seen that the bias of the estimated BLTI approximation ˆ BLTI (j!k ) is smaller than its uncertainty. G Further, the time-varying parts of the input-output DFT spectra Z˜TV (k) and the output residual YTV (k) of the BLTI approximation are estimated via, respectively, (27) and (28). The results are shown in Figures 9 and 11. From Figure 9 it can be seen that the time-variation (top light gray lines) is well above the estimated noise level (bottom medium gray lines), and is of the same order of magnitude as the leakage errors (middle dark gray lines). Note that the input time-variation-to-signal ratio is much larger than that at the
H
(jω )
ru,1
H
k
−80 0.1 0.2 Frequency (Hz)
−40
−80 −100 0
0.3
Hru,2(jωk)
0.1 0.2 Frequency (Hz)
0.3
Hry,2(jωk)
−40 −60 −80 0.1 0.2 Frequency (Hz)
−40
0.1 0.2 Frequency (Hz)
0.3
−40 −60 −80 −100 0
0.3
0
−80 0
−20 Amplitude (dB)
−20 Amplitude (dB)
GSA(jωk) 40
−60 Amplitude (dB)
Amplitude (dB)
Amplitude (dB)
−60
−100 0
k
−20
−40
−100 0
(jω )
ry,1
−20
0.1 0.2 Frequency (Hz)
0.3
Figure 12. Spectral analysis estimate GSA (j!k ) (9) of the best linear time-invariant approximation – simulation example. Black: estimate ˆ SA (j!k ); dark gray: variance estimate var(G ˆ SA (j!k )); light gray: differG ˆ SA (j!k ) G ˆ BLTI (j!k ) |; and black dashed line: var(G ˆ BLTI (j!k )). ence |G
Figure 10. Frequency response functions Hrz,p (j!k ), p = 1, 2, of the MITO equivalent model (22) – simulation example. Black: input (left) and ˆ ry,p (j!k ); light gray: difference ˆ ru,p (j!k ) and H output (right) estimates H between the estimated and the true value of Hrz,p (j!k ); and dark gray: ˆ rz,p (j!k ). estimated variance of H
G
(j
BLTI
k
)
Output residual Amplitude (dB)
Amplitude (dB)
50 0 50 100 0
0.1 0.2 Frequency (Hz)
0.3
0 20 40 60 0
0.1 0.2 Frequency (Hz)
0.3
Figure 11. Estimated best linear time invariant (BLTI) approximation (left) and its output residual (right) – simulation example. Left plot: the estimate ˆ BLTI (j!k ) (black), its variance var(G ˆ BLTI (j!k )) (dark gray), and its bias G ˆ ˆ ˆ |GBLTI (j!k ) GBLTI (j!k ) | (light gray). Right plot: comparison between ˜TV (k) (27) of the time-varying part Y˜TV (k) of the output the estimate Yˆ DFT spectrum (black), and the estimated output residual YˆTV (k) (28) of the BLTI approximation (light gray).
output below 0.15 Hz. From the right plot of Figure 11 it can be seen that YTV (k) is significantly larger than the time-varying part Y˜TV (k) of the output DFT spectrum in the band [0 Hz, 0.15 Hz]. Note that in this frequency band ˆ BLTI (j!k )| |G 1 (see Fig. 11, left plot). It nicely illustrates that a feedback loop suppresses the time-variation at those frequencies where the open loop gain is significantly larger than one. Finally, Figure 12 compares the proposed approach with the indirect spectral analysis estimate (12) obtained via (22) with Nb = 0 (see [11]). It can be seen that the standard deviation of ˆ SA (j!k ) (dark the classical indirect spectral analysis estimate G gray line) is 30 dB above that of the estimated BLTI approximation (black dashed line). It follows that the input-output ˜TV (k) and Y˜TV (k) is the dominant error term time-variation U ˆ SA (j!k ). in G
Figure 13. Time-variant electronic circuit with input u (t), output y (t), and scheduling parameter p (t). Schematic (a): high gain second order bandpass filter with time-varying resonance frequency and damping ratio. It consists of the following components: three high gain operational amplifiers (TL071), three time-varying resistors R (p(t)), four resistors (R1 = R2 = 10 k⌦, R3 = 5.31 k⌦, and R4 = 100.8 k⌦), and two capacitors (C1 = C2 = 10 nF). Schematic (b): practical realization of the time-varying resistor using an operational amplifier, a 330 ⌦ resistor, and an electro-optical component (VTL5C1) consisting of a light dependent resistor (LDR) and a light emitting diode (LED).
Reference DFT spectrum Amplitude (dB)
0
50
100 0
10 20 30 Frequency (kHz)
Input DFT spectrum
Output DFT spectrum
r (t) =
F X
Ak sin (2⇡kf0 t +
k)
(39)
k=1
with uniformly distributed phases k , and equal amplitudes in the frequency band [228.9 Hz, 39.98 kHz] chosen such that the rms value is 50 mV (f0 = fs /P , fs = 625 kHz, P = 16384, A1 = A2 = · · · = A5 = 0, and A6 = A7 = · · · = AF ). A periodic ramp p (t) is applied to the electro-optical components in Figure 13 (period 1/f0 and voltage increasing linearly from 0.8412 V to 0.9317 V). N = 15P/16 = 15360 inputoutput samples of the steady state response are acquired using a band-limited measurement setup (u (t) and y (t) are lowpass filtered before sampling). Figure 15 shows the corresponding reference and input-output DFT spectra (5). During the acquisition time the voltage p (t) increased linearly from 0.8412 V to 0.9259 V. B. Results Starting from N known reference, noisy input-output samples r (nTs ), u (nTs ) and y (nTs ), the multiple-input, twooutput (MITO) model (22), with fp (t) Legendre polynomials satisfying (11), is estimated following the procedure of Section III-C (local polynomial method [16] with a fourth order local polynomial approximation of the FRFs and the leakage
Amplitude (dB)
100 0
40
Hru,1(jωk)
10 20 30 Frequency (kHz)
40
Hry,1(jωk) −20 Amplitude (dB)
−20 −40 −60 −80 0
10 20 30 Frequency (kHz)
−40 −60 −80 0
40
Hru,2(jωk)
40
−20
−40 −60 −80 0
10 20 30 Frequency (kHz) Hry,2(jωk)
−20 Amplitude (dB)
The device under test is a time-variant electronic circuit (see Fig. 13) that behaves as a high gain second order bandpass filter with time-varying resonance frequency and damping ratio. It operates in feedback as shown in Figure 14, and the time-variation is controlled via the voltage p (t) (scheduling parameter) of the electro-optical components (see Fig. 13, bottom schematic). The reference signal r (t) is a random phase multisine consisting of the sum of F = 1043 harmonically related sinewaves
10 20 30 Frequency (kHz)
50
Figure 15. DFT spectra (black) of the known reference (top), noisy input (bottom left), and noisy output (bottom right) signals – time-variant electronic ˆ ˜TV (k) and circuit. Light gray lines: estimated input-output time-variation U ˆ Y˜TV (k) (27); top dark gray lines: estimated input-output leakage errors TˆU (j!k ) and TˆY (j!k ) (22); and bottom medium gray lines: estimated input-output noise variances.
Amplitude (dB)
A. Measurement Setup
0
50
100 0
Amplitude (dB)
VI. M EASUREMENT E XAMPLE
Amplitude (dB)
0
Figure 14. Time-variant electronic circuit (TV circuit) of Figure 13 operating in feedback. r (t), u (t), and y (t) are, respectively, the reference signal stored in the arbitrary waveform generator, the input of the time-variant circuit, and the output of the time-variant circuit. A 50 ⌦ resistor is added at the output to match the 50 ⌦ input impedance of the acquisition channel. Rg is the 50 ⌦ output impedance of the arbitrary waveform generator.
40
10 20 30 Frequency (kHz)
40
−40 −60 −80 0
10 20 30 Frequency (kHz)
40
Figure 16. Frequency response functions Hrz,p (j!k ), p = 1, 2, of the MITO equivalent model (22) – time-variant electronic circuit. Black: input ˆ ru,p (j!k ) and H ˆ ry,p (j!k ); and light (left) and output (right) estimates H ˆ rz,p (j!k ). gray: estimated variance of H
term, and six degrees of freedom for the variance estimates). Figure 16 shows the estimated dynamics Hru,p (j!) and Hry,p (j!) of the time-varying branches for Nb = 2. It can be seen that the time-variation is significant over the whole frequency band. Increasing Nb further makes no sense ˆ rx,3 (j!k ) | ⇠ std(H ˆ rx,3 (j!k )), because for all values of k, |H x = u, y. Figure 15 shows the corresponding estimates of the inputoutput leakage errors TU (j!k ) and TY (j!k ), and the inputoutput noise variances. It can be seen that the output signal-tonoise ratio (SNR) is about 60 dB, while the input SNR varies between 40 dB and 50 dB. Note that the leakage errors lie well above the noise floor. Figure 17 shows the nonparametric estimate of the best lin-
G
(j
GSA(jωk)
Output residual
40
20
10 20 30 Frequency (kHz)
40
Amplitude (dB)
40 20 0 20 40 60 80 0
)
k
Amplitude (dB)
Amplitude (dB)
BLTI
40 60 80 100 0
10 20 30 Frequency (kHz)
40
Figure 17. Estimated best linear time invariant (BLTI) approximation (left) and its output residual (right) – time-variant electronic circuit. Left plot: the ˆ BLTI (j!k ) (black), its variance var(G ˆ BLTI (j!k )) (dark gray), and estimate G ˆ ˆ ˆ its bias |GBLTI (j!k ) GBLTI (j!k ) | (light gray). Right plot: comparison ˜TV (k) (27) of the time-varying part Y˜TV (k) of the between the estimate Yˆ output DFT spectrum (black), and the estimated output residual YˆTV (k) (28) of the BLTI approximation (light gray).
ear time-invariant (BLTI) approximation obtained via Eqs. (25) ˆ ˆ BLTI (j!k ) G ˆ BLTI (j!k ) is and (30). Note that the difference G ˆ an estimate of the bias (26) of GBLTI (j!k ). From Fig. 17 it follows that this bias is smaller than the standard deviˆ BLTI (j!k ) and, hence, can be neglected. Using ation of G ˆ GBLTI (j!k ), the time-varying parts of the input-output DFT spectra Z˜TV (k) are calculated via (27). From Figure 15 it can be seen that the input and output time-variations are well above the noise level and are of the same order of magnitude as, respectively, the input and output leakage errors. Note ˜TV (k) (light also that the observed output time-variation Yˆ gray line) is about 40 dB below the signal level. From this number one could wrongly conclude that the electronic circuit is slowly time-varying. Indeed, the true output time-variation YTV (k), estimated via (28), is much larger than Y˜TV (k) at those frequencies where |GBLTI (j!k) | 1 (see Fig. 17). From Figures 15 and 17 it can be seen that YTV (k) is about 20 dB below Y (k), indicating a large (10%) time-variation of the electronic circuit over the experiment duration. Figure 18 compares the indirect spectral analysis estimate (9), obtained by applying the procedure of Section III-C with Nb = 0 (see [11]), with the estimated BLTI approximation (procedure of Section III-C with Nb = 2). The difference between the variances of both estimates is maximal in the passband of the system. This is exactly the frequency band where the time-variation is the largest. It nicely illustrates the benefit of the proposed procedure: it eliminates the random influence of the time-variation on the FRF measurement. VII. C ONCLUSIONS The main difficulty of estimating the best linear timeinvariant (BLTI) approximation of a linear time-variant (LTV) system operating in closed loop is that both the output and the input signals are disturbed by the time-variation. Indeed, the input time-variation causes a bias in the estimation procedure designed for open loop systems [1]. To solve this problem the concept of the best linear time-invariant (BLTI) approximation of a linear time-variant (LTV) system has been generalized. This generalized definition can also handle time-variant actuator and feedback dynamics. For open loop
20 0 −20 −40 0
10 20 30 Frequency (kHz)
40
Figure 18. Spectral analysis estimate GSA (j!k ) (9) of the best linear time-invariant approximation – time-variant electronic circuit. Black: estimate ˆ SA (j!k ); dark gray: variance estimate var(G ˆ SA (j!k )); light gray: differG ˆ SA (j!k ) G ˆ BLTI (j!k ) |; and black dashed line: var(G ˆ BLTI (j!k )). ence |G
LTV systems driven by linear time-invariant actuators, the generalized definition simplifies to the classical definition [1]. The difference yTV (t) between the actual output of the LTV system and the output of the BLTI approximation quantifies the amount of time-variation of the system. It allows users to decide whether the time-invariant framework (the BLTI approximation) is accurate enough for their particular application. It has been shown in this paper that yTV (t) is equal to the output time-variation minus the input time-variation filtered by the BLTI approximation. This linear correction remains exact for time-variant actuators and feedback dynamics. A nonparametric procedure for estimating the BLTI approximation and the output residual yTV (t) from known reference and noisy input-output samples has been proposed. It is based on a multiple-input, two-output linear time-invariant equivalent model of the open loop linear time-variant system from reference to input-output. In those frequency bands where the time-variation is the dominant disturbance, the variance of the estimated BLTI approximation is much smaller than that of the classical indirect spectral analysis method. If the linear system operating in closed loop is timeinvariant, and the actuator and/or feedback dynamics are timevariant, then the BLTI approximation equals the frequency response function (FRF) within an O(T 1 ) bias error, where T is the experiment duration. However, using the multipleinput, two-output linear time-invariant equivalent model, the time-invariant behaviour can be detected and an unbiased estimate of the FRF can be obtained. ACKNOWLEDGEMENT This work is sponsored by the Research Council of the Vrije Universiteit Brussel, the Research Foundation Flanders (FWO-Vlaanderen), the Flemish Government (Methusalem Fund METH1), and the Belgian Federal Government (Interuniversity Attraction Poles programme IAP VII DYSCO). A PPENDIX A. Proof of Eqs. (9) and (10) The key idea consists in modeling the closed loop system in Fig. 4 from the reference signal r (t) to the input-output signals z (t) = [y (t) u (t)]T simultaneously. As such, the single-input, single-output closed loop problem is transformed into a single-input, two-output open loop problem, for which
Laplace transform of u (t), and U (s) = Gact (s) R (s), it follows that the time-varying transfer function from reference to output is given by Gry (s, t) = G (s, t) Gact (s) . Figure 19. Equivalent one-input, two-output open loop representation of a linear time-variant system operating in closed loop (see Fig. 4).
the results of [1], [11] apply. Using the results of Section II-A we get Z (k) = Grz,BLTI (j!k ) R (k) + Z˜TV (k)
(40)
with Grz,BLTI = [Gry,BLTI Gru,BLTI ]T the 2 ⇥ 1 BLTI approximation of the single-input, two-output LTV, and where ˜TV (k)]T is uncorrelated with – but Z˜TV (k) = [ Y˜TV (k) U not independent of – the reference R (k). Using (40) we find E Y (k) R (k) E U (k) R (k)
=
Gry,BLTI (j!k ) Gru,BLTI (j!k )
(41)
Since Grx,BLTI (j!k ), x = y, u, is – by definition – the BLTI approximation of the open loop LTV system with time-varying transfer function Grx (s, t), we have that ˆ 1 T Grx,BLTI (j!k ) = Grx (j!k , t) dt (42) T 0 (see Eq. (2)). Combining (7), (41), and (42) proves (9). Calculating the difference YTV (k) between the actual output of the closed loop LTV system and the output of the BLTI approximation (7) gives YTV (k)
=
Y (k) GBLTI (j!k ) U (k) ˜TV (k) Y˜TV (k) GBLTI (j!k ) U
=
(43)
where the second equality uses Eq. (40). It shows that ˜TV (k) are YTV (k) is uncorrelated with R (k) (Y˜TV (k) and U both uncorrelated with R (k)), which concludes the proof of Eq. (10). B. Proof of the Properties of the BLTI Approximation (7) In this appendix we first prove that the generalized definition (7) reduces to (2) for open loop systems. Next, we show that the BLTI approximation (7) is independent of the actuator characteristics. Finally, the relationship (12) between the BLTI approximation (7) and the series expansion (8) is proven. If the LTV system operates in open loop (see Fig. 3, top block diagram), then the time-varying transfer function from reference to input is equal to the time-invariant actuator transfer function Gru (s, t) = Gact (s) .
(44)
Using y (t) = L {G (s, t) U (s)} (see [23], [24]), with L 1 {} the inverse Laplace transform operator and U (s) the 1
(45)
Combining Eqs. (7), (44), and (45) gives (2). Figure 19 gives an equivalent one-input, two-output open loop representation of the closed loop LTV system in Figure 4 (top block diagram). Using x (t) = L 1 {Grx (s, t) X (s)} (see [23], [24]) with x = y, u, X = Y, U , and X (s) = Gact (s) R (s) it follows that Gry (s, t)
=
Gwy (s, t) Gact (s)
(46)
Gru (s, t)
=
Gwu (s, t) Gact (s)
(47)
where w (t) is the output of the actuator (see Figures 4 and 19). Combining Eqs. (7), (46), and (47) shows that GBLTI (s) is independent of the actuator characteristics, viz. ´ 1 T T 0 Gwy (s, t) dt GBLTI (s) = ´ T . (48) 1 T 0 Gwu (s, t) dt
In general, the BLTI approximation (48) depends on the feedback transfer function Gfb (s) because, in general, the timevarying transfer function Guq (s, t) from the input u (t) of the LTV system to the output q (t) of the feedback branch (see Figure 4) is different from G (s, t) Gfb (s) (see [23], [24]). Integrating (8) over [0, T ] using (11) gives ˆ 1 T Grz (s, t) dt = Grz,0 (s) (49) T 0 Combining Eqs. (7) and (49) proves (12). C. Proof of the Bias (26) of the BLTI Approximation Estimate ˆ rz,0 (j!k ) (23) can be written The nonparametric estimate G as (to simplify the notations we drop the arguments) ˆ rz,0 = Grz,0 + nG,z G with Grz,0 the true value, and where E{nG,z } = O(T Using (50), Eq. (25) becomes ˆ BLTI = Gry,0 + nG,y = GBLTI 1 + nG,y /Gry,0 G Gru,0 + nG,u 1 + nG,u /Gru,0
(50) 2
).
(51)
For T sufficiently large, one can always choose the order of the polynomial approximation and the width of the local frequency band in the LPM method [16], [17] such that nG,y nG,u ,
