0 x1 x2 Kf e x1 x2 : KfTe ^! = d dt x2 x1 Kfex1 x2 ... - Semantic Scholar
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
[27]
, “Robust control of state delayed systems with polytopic type uncertainties via parameter-dependent Lyapunov functionals,” Syst. Control Lett., vol. 50, no. 3, pp. 183–193, 2003. [28] , “Robust sliding mode control of uncertain time-delay systems : An LMI approach,” IEEE Trans. Autom. Control, vol. 48, no. 6, pp. 1986–1092, Jun. 2003. [29] C. Yuan and X. Mao, “Robust stability and controllability of stochastic differential delay equations with Markovian switching,” Automatica, vol. 40, pp. 343–354, 2004.
Noise Analysis of an Algorithm for Uncertain Frequency Identification Qing Zhang and L. J. Brown Abstract—This note presents a noise analysis for an algorithm to identify the uncertain frequency of periodic signals or disturbances. This algorithm is based on the time-varying states of an internal model principle controller which can be mapped nonlinearly to the frequency and the magnitude or energy of the periodic signal or disturbance. This note provides an analysis of the ’measurement’ of this frequency in the presence of white noise. In the case of an additive white noise, we develop some formulas to calculate the means and variances of the measured difference between the true frequency and nominal frequency for high and low signal-to-noise ratio (SNR). When an integral controller is used to eliminate this difference, we prove that this frequency estimation is unbiased. The formulae to calculate the mean and variance are also given for the output of the integral controller. The simulations verify the validity of approximations used in our noise analysis. Index Terms—Adaptive notch filter, Cramer–Rao bound, frequency estimation, internal model principle, periodic disturbance.
103
adaptive algorithm for Regalia’s method based on traditional adaptive control scheme. In [1], it was proposed that the two states of an internal model principle controller (IM) could be mapped into a ’measurement’ of the signal frequency. This measurement resembles a normalized version of Regalia’s frequency update. The error between this measurement and the IMs frequency is fed into a integrator or other estimator which updates the frequency of the IM. We can show that the two state variables x1 t and x2 t have the same relationship to the Hilbert transform as Stratonovich suggested [14]. Under appropriate assumptions transfer functions at x1 t and x2 t are narrow band. In Section II, we present the algorithm and the main results. The analysis is conducted in Section III. We first provide a noise analysis for the frequency measurement in the presence of white noise when the internal model controller’s intrinsic frequency is fixed and not equal to the true frequency of the external signal. This noise analysis will provide us a good guideline to select an appropriate estimation filter to achieve an accurate parameter estimate. Next we will give the mean and variance of the output of a simple update algorithm driven by this measurement and show that this algorithm is unbiased. This result is then compared with Cramer-Rao bound.In Section IV, simulations verify the results in Section II. In Section V, some conclusions are drawn. In particular, we conclude that the variance calculated in this note will apply in the periodic disturbance cancellation case.
()
() ()
II. ALGORITHM AND RESULTS The block diagram of an IM based frequency identification method that we use for noise analysis is shown in Fig. 1. We assume that the a !c t ' n, where a, !c and ' are uninput signal is d t known and n is a white noise with a constant spectral density function N0 . The state–space form for IM is
( ) = sin( + ) + x_ 1 x_ 2
I. INTRODUCTION This note provides an analysis of the behavior of the frequency identification algorithm of [1] in the presence of noise. Noise contamination extensively exists in a variety of fields such as control system, communication, signal processing, parameter identification, etc. When we want to estimate the parameters of periodic signals or disturbances (e.g., magnitude, frequency, and phase), it is important to understand the effects of noise. Since many of the approaches to estimate parameters are derived in the noise free case, they may not work well with the presence of noise. In particular, many algorithms fail in the presence of constant offsets. There exist several methods to identify frequency: 1) adaptive notch filter [2]–[4]; 2) time frequency representation (TFR) based method [5]; 3) phase locked loop (PLL) based method [6]; 4) eigensubspace tracking estimation [7]; 5) extended Kalman filter frequency estimation [8]; 6) internal model based method [9], [10]. Regalia [2] presented an adaptive frequency estimation algorithm most similar to our algorithm based on notch filter. He uses a gradient based algorithm to update the parameter in the notch filter. Hsu et al. [3] proposed a normalized
()
= !0 00! x y = [ 0 1] x
+ K0
x1 x2
e
f
:
1 2
(1) (2)
A discrete-time implementation of the proposed adaptive algorithm in this note was presented in [16]. For piecewise constant e, the states x1 and x2 can be calculated at time instant kT by the discrete state equation
x1 x2
;kT
;kT
sin !T = 0cossin!T!T cos !T
+ K0 T e(kT )
x1 ( 01) x2 ( 01) ; k
T
; k
T
(3)
f
where T is the sampling time. Then in steady state, an estimation of the frequency !c is
!^ = d arctan x2 dt x1
= ! + xK +exx f
2 1
1 2 2
:
(4)
The second term in (4) is the error between the true and estimated frequencies. This term is used for the noise analysis and we define
1! = !^ 0 ! = xK +exx : The transfer function, T , from d(t) to e(t) is L(s)(s + ! ) T = (s + ! ) + K sL(s) and L(s) and K are chosen so that s +! T =T 2 s + 2"!s + ! 1 2 2
f
2 1
(5)
de
2
Manuscript received June 9, 2004; revised February 23, 2005 and July 28, 2005. Recommended by Associate Editor E. Bai. This work was supported by NSERC. The authors are with the Department of Electrical and Computer Engineering, The University of Western Ontario, London, ON N6A 5B9, Canada (e-mail: [email protected]). Digital Object Identifier 10.1109/TAC.2005.861712
de
2
de
bp
2
2
(6)
f
f
2
2
2
2
(7)
by matching their corresponding numerators and denominators, as in [10], where Tbp is a symmetric bandpass filter which has the center
0018-9286/$20.00 © 2006 IEEE
104
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
frequency ! and bandwidth BW and " is a small positive number which can control the width of the notch. BW must be chosen less than !c . Similar equations apply to the discrete time case. It will be shown that the variance in e, for ideal bandpass filters, is given by
e2 = (BW 0 2"!) 2 N0 :
nc (t), and ns (t) have the same variance x2 = "!N0 and mean 0. An
alternative representation [13] for the noise is
n1 = Xc (t)c(t; 0) 0 Xs (t)s(t; 0) where
Xc (t) = nc (t)cos (! 0 !c )t
0 ns (t)sin (! 0 !c)t Xs (t) = ns (t)cos (! 0 !c )t + nc (t)sin (! 0 !c )t :
(8)
1
Further, although the measurement of ! is biased in the presence of white noise, if an integral controller is used to update the nominal frequency, ! Ke ! or discrete equivalent, the frequency estimation is unbiased and its variance is
_= 1
V ar[!] =
(
0 Ke2 "!N a2
) 2
(9)
[ ]
as long as a = "!N0 > and " V ar ! =! . Note that in [1] the local exponential stability for this algorithm was shown. Throughout we assume that we remain in the region of attraction of this result. 2
A. Noise Model To simplify our mathematical expressions we will use the notation
s(t; ) = sin(wc t + ) and c(t; ) = cos(wc t + ). We use the symbol to define new variables hk , k = f1; 2; 3; 4; 5; 6g. The transfer functions for x1 , x2 from d(t) are Tdx Tdx
()
f! = 0 Tbp 2 s2 + 2K"!s + !2 K s f = Tbp 2 s2 + 2"!s + !2 :
(10)
[ ( ) ( )] = 0
()
x = "!N0 : 2
In steady state, the responses at x1 , x2 and e
(12)
(t) are
x1 = ar c(t; ' + ) + n1 x2 = ar !c s(t; ' + ) + n2 ! e = 0 ar (!2 0 !c2 )c(t; ' + ) + ne Kf !
ne =
(20)
0 K1f ! [(!2 0 !c2) Xc (t)c(t; 0) 0 Xs (t)s(t; 0) 0 2!cX_ c (t)s(t; 0) 0 2!cX_ s (t)c(t; 0) + Xc (t)c(t; 0) 0 Xs (t)s(t; 0)]:
(22)
_ () _ ()
Since Xs t , Xc t vary slowly in relation to center frequency ! , we neglect the last two terms in (22) and obtain
ne 0
1 [(!2 0 !c2)
Kf !
Xc (t)c(t; 0) 0 Xs (t)s(t; 0)
02!cX_ c (t)s(t; 0) 0 2!cX_ s (t)c(t; 0)]:
(23)
(14)
B. Mean and Variance Estimations for
1! = f (!; t) + n(!; t)
(24)
where
A(t)(!c2 0 !2 )!c2 (t; (t)) A(t)[!2 c2 (t; (t)) + !c2 s2 (t; (t))] 2!!cc(t; (t)) X_ c (t)s(t; 0) + X_ s (t)c(t; 0) n(!; t) = A(t)[!2 c2 (t; (t)) + !c2 s2 (t; (t))] f (!; t) =
(16)
where nc t and ns t are in-phase component and quadrature component, respectively, and they are uncorrelated. If n1 is zero-mean Gaussian noise, the quadrature components nc t and ns t of the noise are statistically independent zero-mean Gaussian processes. n1 ,
1!
Substituting (13)–(15), (17), (21), and (23) into (5), we get
[13]
()
!
(13)
= (=2) + 6 D(j!c ), ar = a!jD(j!c )j, and D(s) = Kf L(s)=((s2 + !2 )+ Kf sL(s)). n1 , n2 and ne are the filtered white noises at x1 , x2 , and e(t), respectively. Since Tdx , Tdx and Tde are approximately symmetric about ! , the filtered noise at x1 is a narrow band noise [12] and can be expressed as
()
()
() 2
nc (t), ns (t) and Xc (t), Xs (t) are slowly varying envelops in relation to center frequency ! [14]. Therefore, we can omit the last two terms in n2 and get n2 !c [Xc (t)s(t; 0) + Xs (t)c(t; 0)]: (21) ! The relationship between e and x1 is e = 0(1=Kf ! )(d2 x1 =dt2 + 2 ! x1 ). Therefore
where
()
()
Since the transfer function Tde is an ideal bandpass filter with a narrow notch, we can get the variance for ne given in (8) [11].
(15)
n1 = nc (t)cos(!t) 0 ns (t)sin(!t)
(19)
()
()
[ ( ) _ ( )] = 0 [ ( ) _ ( )] = 0 = (1 )_ n2 = !c [Xc (t)s(t; 0) + Xs (t)c(t; 0)] ! _ _ 0 Xc (t) c(t; 0) + Xs (t) s(t; 0)
=2
2
()
()
(11)
Here, we assume Tbp is a ideal bandpass filter with magnitude 1 and bandwidth BW. In our design [10], we get Kf "! . Therefore,when " is small enough, Tdx and Tdx will be approximately narrow first order bandpass filters with maximal magnitude 1 and bandwidth "! . The variances of noise at x1 and x2 are approximately equal and can be expressed as [11]
()
()
(18)
In this case, Xc t and Xs t will be Gaussian if nc t and ns t are Gaussian. Obviously, by the definitions for Xc t and Xs t in (18) . Thus, Xc t and Xs t and (19), we know that E Xc t Xs t are orthogonal and slowly varying in relation to !c since nc t and ns t are slowly varying in relation to !c and j! 0 !c j < BW= . We can also get E Xc t Xc t and E Xs t Xs t . Since x2 0 =! x1 , n2 can be written as follows:
!
III. NOISE ANALYSIS
(17)
(25)
(26)
and
A(t)=
ar cos(' + )+ Xc (t)
2
+ ar sin(' + )+ Xs(t)
2
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
(t)= tan01 ar sin(' + ) + Xs (t) : ar cos(' + ) + Xc (t)
()
The joint pdf for A and is (see [12, Ch. 4])
0 21 2 [A2 + a2r 0 2Aar cos( 0 (' + )] x
A : 2 2 2 x
The pdf of the phase noise, [12]
[ ( )] () ( )
Although E f !; t can be approximated as !c 0 ! by statistical and averaging methods, in the numerator of f !; t , the magnitude A t contributed by e t can not simply be cancelled by A t in the denominator of f !; t in the case of noise. We consider them different, that is, we assume that A t and =A t are independent. Therefore, we obtain E f !; t as
and
p(A; ) = exp
105
(28)
=0 1
(29)
12
()
2 SNR = a : N0
E
(30)
[n(!; t)] can E
(33)
1 A(t)
:
1) The Expectation of 1! : When SNR = a2 =N0 ! 0, the useful signal is so weak that we can not detect !c . Therefore, we use (5), (16), (21), and (23) to estimate the mean E [1! ] directly
E [1!]=E 2cos(!t + (t)) 2 (_npc (t2)sin(2!t)+n_ s (t)cos(!t)) n c + ns = 0: (31) If SNR is high, we have E [Xs2 (t)] = E [Xc2 (t)] = E [n2s (t)] E [n2c (t)] = x2 , we get (see the Appendix for details)
E [n(!; t)] = (! 0 !c ) 2 E [Xc2 (t) + Xs2 (t)] 2 E 2 1 A(t) 2 2 2 = 0 (!c 0 !) exp 0 4arx2 I0 4arx2 0 (!c 0 !) 2 h3 :
=
(34)
[1 ] = 1 0 [1 ] = 0 [1 ] = 1 [1 ] = 0 2 [1 ] = 2( ) 1 0
E [1!] = h5 0 1 h3 2 (!c 0 !) 2 2 2 2 a r = 42 exp 0 2ar2 I0 4ar2 x x x 2r 2r a a 2 I 0 4 2 + I 1 4 2 (! c 0 ! ) x x h 6 2 (! c 0 ! ) :
()
2c(t; (t)) X_ c (t)s(t; 0) + X_ s (t)c(t; 0)
2 I0 ar2 4x
While ideally, we would like E ! ! as SNR ! , it is also true that when a , ! is undefined and it is best that E ! . ! for large SNR to E ! Thus a smooth transition of E ! for no signal is desirable. When SNR ! , h5 ! = and h3 ! , 0= !c 0 ! instead and the high SNR based estimate of E ! of 0. When SNR ! 1, h5 ! and h3 ! . Thus, scaling h3 by = in (34) will not affect the calculation for high SNR and cause the result to be consistent with the low SNR result. Thus
where I0 is the modified zeroth order Bessel function. Formula (29) is called Ricean distribution [13]. Here, we assume A t and t are uncorrelated for all input signal-to-noise ratio (SNR)
E [n(!; t)]
2
E [1!] = E [f (!; t)] + E [n(!; t)] = (h5 0 h3 ) 2 (!c 0 !):
0 21 2 (A2 + a2r ) I0 A 22ar x x
Therefore, by using statistical and averaging methods, E be approximated as
1
A
Therefore
plus narrow band Gaussian noise is
()
()
1 ()
= (!c 0 !) 2 exp 0 2ar2 x 2r 2r 2 2 a a 2 1 + 22 I0 42 + 2ar2 I1 4ar2 x x x x (! c 0 ! ) hs :
(27)
where the error function (erf ) is defined as erf (x) = (2=p) 0x exp(0z2 )dz . The pdf of the envelope A(t) of a sine wave 0
pA (A) = A2 exp x
() [ ( )]
E [f (!; t)] = (!c 0 !) 2 E [A] 2 E
1 = (t) 0 (' + ), can be obtained by
2 cos(1) p1 (1) = 1 exp 0 ar2 + ar p 2 2x 2 2x 2r 2 a 2 exp 0 22 sin (1) x cos(1) 2 1 + erf ar p 2x2
( )
1! :
2) The Variance of The variance of ! is
1
Next, we calculate the variances of
V ar[1!] = E [(1! 0 E [1!])2 ]:
=
(35)
1! . (36)
When ! 6 !c , although the average value of the (25) is a constant, there is a high frequency oscillation component ( !c ) in f !; t . This high frequency component can be obtained by
2
( )
2 2 f (!; t) 0 (!c 0 !) !c 0 ! cos2 (!c t + (t)) 0 (!c 0 !) !c (37) (!c 0 !)cos 2(!ct + (t)) : Therefore
V ar[1!] = E [f 2 (!; t)] 0 2E [f (!; t)]E [n(!; t)] + E [n2(!; t)] 0 E 2 [1!] = 32 h12h24 0 2h3 h5 0 h62 (!c 0 !)2 + 2("!)2 e2h2 (38)
(32)
where h1 , h2 , h4 are defined in the appendix.
106
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
Fig. 1. Block diagram of an internal model control system.
D. Comparison With Cramer–Rao Bound
C. Noise Analysis for Integrated Error The simplest means of updating ! is to use an integral controller. !c The integral controller will converge to the correct answer, ! and ar a, as long as the sign of h5 0 h3 is positive, that is
=
=
a2r > max x2 I0
2I0
a
4
= 2:
+ I1 This inequality means that the power a2 =2 (x21 + x22 )=2 should be 2 a 4
a 4
greater than the variance x in the states of IM in order to make the algorithm work. Since the integral operator has the function of a low pass filter, the high frequency components in the estimation of ! can be neglected. Thus, by using the result in (45) in the Appendix , the noise in the output of the integral controller is
1
X_ s ( )cos( ( )) 0 X_ c ( )sin( ( )) d !c 0 ! = Ke A( ) 01
The Cramer–Rao bound (CRB) is a lower bound on the error variance of any unbiased estimate. The CRB provides a standard against which to measure estimation algorithm performance. We are interested in comparing the behavior of the estimated frequency produced by our integral controller, which behaves like ! 0Ke ! Ke !c and Ke = . In time domain, it is the convolution of Ke e0K t and !c . In order to compare with our method, we calculate the CRB for a sinusoid signal with a time-varying SNR with !c 2 ; 1 and Kn e t 2 . The choice of this SNR results in a time varying weighting of the signal that will result in tracking behavior identical to our algorithm, at the expense of including more noise. Thus, it is possible to achieve better results than the presented CRB. If a white noise has an autocorrelation N0 , the Cramer–Rao bound in this case is [12] function R
=
()
()
where Ke is an adaptation gain. Since SNR is high, A t and t are slowly varying functions and we can approximate them as constants E =A t and E t . So
[1 ( )] [ ( )] !c 0 ! = Ke [Xs (t)cos(E [ (t)]) 0Xc (t)sin(E [ (t)])]E A1(t)
( )=
:
(40)
+ E [Xc2(t)]sin2 (E [ (t)])] 2 E 2 A1 a2r 4x2 :
(41) 2 p For high SNR, since I0 (x) ex = 2x and x2 = "!c N0 in (41), we
have the result (9). When we reduce Ke , ", N0 or increase a, we can reduce the output variance. This variance is independent of the bandwidth BW of Tde , though BW < ! for solution to controller equations to be feasible. =Ke . ThereThe convergence time constant for the algorithm is tc fore, the larger Ke , the smaller the convergence time tc . There exists a tradeoff between fast tracking speed and small noise in frequency estimation. From the previous analysis, note that we just require that the noise has constant spectral density N0 in the passband BW of the closed-loop transfer function Tde . Therefore, the noise used in this note could be a band limited white noise, with a constant component.
=1
(42)
In this section, we show some results of simulations to verify our noise analysis. In other words, we show the various approximations we have made have minimal effect on the accuracy of (34), (35), (38), and (41). First, we select a second-order Chebyshev low-pass filter and obtain a desired bandpass filter with a notch 2 2 Hbpn = 1:102510BW s Q (s )
V ar[!c 0 !] = Ke2 [E [Xs2 (t)]cos2 (E [ (t)])
2
2 !2 4 2N0 : a
)
The previous formula states that smaller will give a smaller Cramer-Rao bound. In the next section, we will show the comparison of the variance of our method with CRB.
integral controller is
2 0 2ar2 I02 x
()
(
IV. SIMULATIONS
Since Xs and Xc are Gaussian random with mean zero and variance x2 , the mean of !c 0 !o is zero. The variance of the output of the
2 = Ke exp
+
(0 )
t
(39)
_=
2
+! 2 s2 +s2"!s + !2 2
2
(43)
where
Q(s) = s4 + 1:097734BW s3 + (2!2 + 1:102510BW 2 )s2 +1:097734!2BWs + !4 : By using the designing method for the desired bandpass filter with a b1 s2 b2 s b3 =s2 , notch for Tde in [10], we get C s 2 2 Gs : BW = s a1 s a2 , L s Gs= G s C s , and an IM Kf s= s2 !2 , where a1 : BW "! , a2 : BW 2 "! !2 , 3 2 b1 , b2 "! : ! BW = : BW 2 , 4 2 b3 ! = : BW , and Kf "! . Since our bandpass filter design is based on a second-order Chebyshev low-pass filter with p and is far from a ideal bandpass filter, we approximate (see [11])
() = ( + + ) ( ) = 1 102510 ( + + ) () = ( + ) ( ) (1 + ( ) ( )) = 1 097734 + 2 = 1 097734 2 +2 = 1 = (2 + 1 097734 ) 1 102510 = 1 102510 =2 = 1 dB e2 = (BW 0 2"!) 2 N0 :
2
(44)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
Fig. 2.
Theoretical and simulated means versus the a = for a = 1, ! = 2
Fig. 3. Theoretical and simulated variance versus the a = , a = 1, ! = 2
107
2 490, !
2 490, !
= 2 2 500, BW = 300, " = 0:01.
= 2 2 500, BW = 300, " = 0:01, = 10
.
108
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
= for the output of integral controller with a = 1, ! = 2 2 100, K = 5, BW = 200, " = 0:02.
Fig. 4.
Theoretical and simulated variances versus the a
Fig. 5.
Comparison of the CRB and the variance of the proposed frequency estimation method.
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
The simulated variance was based on a window of data which has 50 s width. Figs. 2 and 3 show the simulated mean and variance versus the a2r =x2 and the corresponding theoretical mean and variance for the initial value ! = 2 2 490. We can see that the theoretical results (34), (35) and (38) tally with the simulated results for a2r =x2 2. The modified formula (35) is better than (34). Fig. 4 shows the relationship between the variance of the output of the integral controller and a2 =x2 . A discrete implementation of the algorithm given by (3) was implemented with !c = 2 2 100, Ke = 5, BW = 200, " = 0:01, signal to noise ratio a2 =No = 10 000 and sampling rate of 0.5 ms. G(s) and C(s) were converted to discrete representations via bilinear transformation. The measured variance of estimated frequency was 0.0162 rad=s2 versus the 0.0157 rad=s2 given by (9). The comparison of the CRB and the variances of our frequency estimation method with " = 0:01, 0.001 and !c = 2 2 100 is shown in Fig. 5. From Fig. 5, we know that if we reduce the value of ", the variance of the estimated frequency will approach to and then exceed the corresponding CRB for no a priori knowledge of ! . But for smaller ", the region of convergence will be reduced. In particular, our algorithm remains well behaved only for "! > V ar[!]. Therefore, there exists a tradeoff between small estimate variance and big region of convergence. We have extended our approach to the cases of harmonics and multitone signal in [15] and [1], respectively, and shown the corresponding simulated results. In [1], we also showed the performance of our algorithm to tracking a pseudo periodic signal with time-varying amplitude and frequency. In [10], [15], and [16], we compared our algorithm to the methods proposed in [5], [6], and [9]. It is found that our algorithm tracks time-varying frequency with less delay and less susceptible to noise.
109
1 E 2c(t; (t)) 2 X_ c (t)s(t; 0) + X_ s (t)c(t; 0) E A(t) = E[X_ c (t) (sin(2!c t + (t)) 0 sin( (t))) + X_ s (t) (cos(2!c t + (t)) + cos( (t)))]E
APPENDIX A. Approximation of E[n(!; t)] Since
1 TE !!c dt=1: T 0 !2 cos2 (!c t + (t)) + !c2 sin2 (!c t + (t)) E [n(!; t)] can be approximated as E[n(!; t)]
(45)
We neglect the high frequency components in the previous equation by averaging and get
E [n(!; t)] = (! 0 !c ) 2 E Xc2 (t) + Xs2 (t)
1 : 2 E 2 A(t)
(46)
B. Expectations of 1=A and 1=A2 In the general case, the v th moments of a Rice pdf can be analytically expressed as a function of the confluent hypergeometric (Kummer’s) function of the first kind 1 F1 [13]
E [Av ] =
1
0
Av pA (A)dA
2 = (2x2 )v=2 0 1 + v 1 F1 0 v ; 1; 0 ar2 (47) 2 2 2x 1 where 0(c) = 0 c01 exp(0 )d is Gamma function. From for-
mula (47), we can obtain
2 E [A] = x exp 0 ar2 2 4x 2r a a2r + a2r I a2r 2 1 + 22 I0 4 2x2 1 4x2 x x2 h4 2 E A = a2r + 2x2
V. CONCLUSION We show that although the measurement 1! in our algorithm is biased, the estimate ! of the update law is not biased and a formula for its variance is given. Even with constant offsets in d, because of the presence of the bandpass filter Tbp , signals n1 , n2 and e are guaranteed to be zero mean, ensuring our algorithm, unlike many other algorithms, is unaffected. A necessary condition to make our frequency estimation algorithm work appropriately is that the signal power must be greater than noise power at the states of the IM. When the values of the designed parameters " and Ke are reduced, the variance of the estimated frequency can be decreased. But the convergence rate of the identified algorithm declines with the reduction of " and Ke and region of rapid convergence decreases with ". The noise model presented here will allow optimal selection of gains for this algorithm. The results of this note can be applied directly to feedback problem when plant is low pass filter and disturbance is bandlimited white noise. When the disturbance is band limited, the Fourier transform of e(t) will be bandpass with notch just as we require.
1 : A(t)
(48) (49)
where I1 () is the first order modified Bessel function. When v = 01 and using the integral representation [17, eq. (13.2.1)] for 1 F1 , integral formula [18, eq. (3.364)], we have
2 E 1 = 2x2 01=2 0 1 1 F1 1 ; 1; 0 ar2 A 2 2 2x 2r 2r a a = exp 0 2 I0 2x2 4x 4x2 h1 :
(50)
In practice, we add a small positive constant real number in the denominator of 1! . For purpose of analysis, we calculate E 1=(A2 + ) . Considering very small (1008 in our application) and using the [18, eqs. (3.351), (3.352), and (8.214)], we have
1 A2 + 1 1 = pA (A)dA 0 A2 + exp 0 2a = 2x2 a2r 2 exp 2 2 E1 2 2 + Ei 2 x x x2 h2
E
1
a2r 0 C 0 ln 2 2 x
where E1 (z) = z (exp(0t)=t)dt and Ei (z) 0 lim"!+0 [ 00z" (exp(0t)=t)dt + "1 (exp(0t)=t)dt]
(51)
= are
110
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
exponential integrals with z > 0, C = 0:577 215 664 9 1 1 1 is Euler’s constant. The formula (51) is finite for > 0.
Transition-Time Optimization for Switched-Mode Dynamical Systems
REFERENCES
Magnus Egerstedt, Yorai Wardi, and Henrik Axelsson
[1] L. Brown and Q. Zhang, “Identification of periodic signals with uncertain frequency,” IEEE Trans. Signal Process., vol. 51, no. 6, pp. 1538–1545, Jun. 2003. [2] P. A. Regalia, “An improved lattice-based adaptive IIR notch filter,” IEEE Trans. Signal Process., vol. 39, no. 9, pp. 2124–2128, Sep. 1991. [3] L. Hsu, R. Ortega, and G. Damm, “A globally convergent frequency estimator,” IEEE Trans. Autom. Control, vol. 44, no. 4, pp. 698–713, Apr. 1999. [4] S. Bittanti, M. Campi, and S. Savaresi, “Unbiased estimation of a sinusoid in colored noise via adapted notch filters,” Automatica, vol. 33, no. 2, pp. 209–215, 1997. [5] H. Kowok and D. Jones, “Improved instantaneous frequency estimation using an adaptive short-time Fourier transform,” IEEE Trans. Signal Process., vol. 48, no. 10, pp. 2964–2972, Oct. 2000. [6] B. Wu and M. Bodson, “A magnitude/phase-locked loop approach to parameter estimation of periodic signals,” IEEE Trans. Signal Process., vol. 48, no. 10, pp. 2964–2972, Oct. 2000. [7] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood and Cramer-Rao bound,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 5, pp. 720–741, May 1989. [8] B. F. Scala and R. Bitmead, “Design of an extended Kalman filter frequency tracker,” IEEE Trans. Signal Process., vol. 44, no. 3, pp. 739–742, Mar. 1996. [9] T. Tsao, Y. Qian, and M. Nemani, “Repetitive control for asymptotic tracking of periodic signals with an unknown period,” J. Dyna. Syst., Measure., Control, vol. 122, pp. 364–369, Jun. 2000. [10] Q. Zhang and L. J. Brown, “Designing of adaptive bandpass filter with adjustable notch for frequency demodulation,” in Proc. 2003 Amer. Control Conf., vol. 4, Jun. 4–6, 2003, pp. 2931–2936. [11] R. G. Brown, Introduction to Random Signal Analysis and Kalman Filtering. New York: Wiley, c1983. [12] A. D. Whalen, Detection of Signals in Noise. New York: Academic, 1971. [13] N. Wax and Ludeman, Selected Papers on Noise and Stochastic Processes. New York: Dover, 1954. [14] R. L. Stratonovich, Topics in the Theory of Random Noise. New York: Gordon and Breach, pp. 1963–1967. [15] L. J. Brown and Q. Zhang, “Periodic disturbance cancellation with uncertain frequency,” Automatica, vol. 40, no. 4, pp. 631–637, Apr. 2004. [16] Z. Zhao and L. J. Brown, “Fast estimation of power system frequency using adaptive internal model control technique,” in Proc. 43rd IEEE Conf. Decision and Control, Dec. 14–17, 2004, pp. 845–850. [17] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. Washington, DC: U.S. Government, 1970. [18] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products. New York: Academic, 1980.
Abstract—This note considers the problem of determining optimal switching times at which mode transitions should occur in multimodal, hybrid systems. It derives a simple formula for the gradient of the cost functional with respect to the switching times, and uses it in a gradient-descent algorithm. Much of the analysis is carried out in the setting of optimization problems involving fixed switching-mode sequences, but a possible extension is pointed out for the case where the switching-mode sequence is a part of the variable. Numerical examples testify to the viability of the proposed approach. Index Terms—Gradient-descent algorithms, hybrid systems, optimal control, switched dynamical systems, switching-time control.
I. INTRODUCTION Switched dynamical systems are often described by differential inclusions of the form
x_ (t) 2 fg (x(t); u(t))g 2A
(1)
where x(t) 2 n , u(t) 2 k , and fg : n+k ! n g 2A is a collection of continuously differentiable functions, parameterized by belonging to some given set A. The time t is confined to a given finite-length interval [0; T ]. Such systems arise in a variety of applications, including situations where a control module has to switch its attention among a number of subsystems [12], [15], [19], or collect data sequentially from a number of sensory sources [4], [6], [11]. A supervisory controller is normally engaged for dictating the switching law, i.e., the rule for switching among the functions g in the right-hand side of (1). Recently, there has been a growing interest in optimal switching time control of such hybrid systems, where the control variable consists of a proper switching law as well as the input function u(t) (see [3], [5], [9], [10], [16]–[18], and [20]). Reference [3] establishes a framework for optimal control, and [16]–[18] present suitable variants of the maximum principle to the setting of hybrid systems. [1], [2], [9], [14] consider piecewise-linear or affine systems. The special case of autonomous systems, where the term u(t) is absent and the control variable consists solely of the switching times, is considered in [9], [11], [21], [22]. In particular, Xu and Antsaklis [21], [22] consider general nonlinear systems, and they have developed nonlinear-programming algorithms that compute the gradient and second-order derivatives of the cost functional. This note, whose preliminary version has appeared in [8], also falls in this category. [21] provides the starting point for the results presented in this note, as we initially consider a similar problem, where the sequence of switching functions as well as the number of switching times are fixed. We develop a simpler formula than the one in [21] for the gradient of the cost functional, and use it in a gradient-descent algorithm. Finally, we suggest a possible extension of the algorithm to a class of problems Manuscript received November 11, 2003; revised December 15, 2004 and June 14, 2005. Recommended by Associate Editor J. Hespanha. This work was supported in part by the National Science Foundation under Grants 0207411 and 0509064, and by a grant from the Georgia Institute of Technology Manufacturing Research Center. The authors are with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332 USA (e-mail: [email protected]; [email protected]; [email protected]). Digital Object Identifier 10.1109/TAC.2005.861711
0018-9286/$20.00 © 2006 IEEE
[27]
, “Robust control of state delayed systems with polytopic type uncertainties via parameter-dependent Lyapunov functionals,” Syst. Control Lett., vol. 50, no. 3, pp. 183–193, 2003. [28] , “Robust sliding mode control of uncertain time-delay systems : An LMI approach,” IEEE Trans. Autom. Control, vol. 48, no. 6, pp. 1986–1092, Jun. 2003. [29] C. Yuan and X. Mao, “Robust stability and controllability of stochastic differential delay equations with Markovian switching,” Automatica, vol. 40, pp. 343–354, 2004.
Noise Analysis of an Algorithm for Uncertain Frequency Identification Qing Zhang and L. J. Brown Abstract—This note presents a noise analysis for an algorithm to identify the uncertain frequency of periodic signals or disturbances. This algorithm is based on the time-varying states of an internal model principle controller which can be mapped nonlinearly to the frequency and the magnitude or energy of the periodic signal or disturbance. This note provides an analysis of the ’measurement’ of this frequency in the presence of white noise. In the case of an additive white noise, we develop some formulas to calculate the means and variances of the measured difference between the true frequency and nominal frequency for high and low signal-to-noise ratio (SNR). When an integral controller is used to eliminate this difference, we prove that this frequency estimation is unbiased. The formulae to calculate the mean and variance are also given for the output of the integral controller. The simulations verify the validity of approximations used in our noise analysis. Index Terms—Adaptive notch filter, Cramer–Rao bound, frequency estimation, internal model principle, periodic disturbance.
103
adaptive algorithm for Regalia’s method based on traditional adaptive control scheme. In [1], it was proposed that the two states of an internal model principle controller (IM) could be mapped into a ’measurement’ of the signal frequency. This measurement resembles a normalized version of Regalia’s frequency update. The error between this measurement and the IMs frequency is fed into a integrator or other estimator which updates the frequency of the IM. We can show that the two state variables x1 t and x2 t have the same relationship to the Hilbert transform as Stratonovich suggested [14]. Under appropriate assumptions transfer functions at x1 t and x2 t are narrow band. In Section II, we present the algorithm and the main results. The analysis is conducted in Section III. We first provide a noise analysis for the frequency measurement in the presence of white noise when the internal model controller’s intrinsic frequency is fixed and not equal to the true frequency of the external signal. This noise analysis will provide us a good guideline to select an appropriate estimation filter to achieve an accurate parameter estimate. Next we will give the mean and variance of the output of a simple update algorithm driven by this measurement and show that this algorithm is unbiased. This result is then compared with Cramer-Rao bound.In Section IV, simulations verify the results in Section II. In Section V, some conclusions are drawn. In particular, we conclude that the variance calculated in this note will apply in the periodic disturbance cancellation case.
()
() ()
II. ALGORITHM AND RESULTS The block diagram of an IM based frequency identification method that we use for noise analysis is shown in Fig. 1. We assume that the a !c t ' n, where a, !c and ' are uninput signal is d t known and n is a white noise with a constant spectral density function N0 . The state–space form for IM is
( ) = sin( + ) + x_ 1 x_ 2
I. INTRODUCTION This note provides an analysis of the behavior of the frequency identification algorithm of [1] in the presence of noise. Noise contamination extensively exists in a variety of fields such as control system, communication, signal processing, parameter identification, etc. When we want to estimate the parameters of periodic signals or disturbances (e.g., magnitude, frequency, and phase), it is important to understand the effects of noise. Since many of the approaches to estimate parameters are derived in the noise free case, they may not work well with the presence of noise. In particular, many algorithms fail in the presence of constant offsets. There exist several methods to identify frequency: 1) adaptive notch filter [2]–[4]; 2) time frequency representation (TFR) based method [5]; 3) phase locked loop (PLL) based method [6]; 4) eigensubspace tracking estimation [7]; 5) extended Kalman filter frequency estimation [8]; 6) internal model based method [9], [10]. Regalia [2] presented an adaptive frequency estimation algorithm most similar to our algorithm based on notch filter. He uses a gradient based algorithm to update the parameter in the notch filter. Hsu et al. [3] proposed a normalized
()
= !0 00! x y = [ 0 1] x
+ K0
x1 x2
e
f
:
1 2
(1) (2)
A discrete-time implementation of the proposed adaptive algorithm in this note was presented in [16]. For piecewise constant e, the states x1 and x2 can be calculated at time instant kT by the discrete state equation
x1 x2
;kT
;kT
sin !T = 0cossin!T!T cos !T
+ K0 T e(kT )
x1 ( 01) x2 ( 01) ; k
T
; k
T
(3)
f
where T is the sampling time. Then in steady state, an estimation of the frequency !c is
!^ = d arctan x2 dt x1
= ! + xK +exx f
2 1
1 2 2
:
(4)
The second term in (4) is the error between the true and estimated frequencies. This term is used for the noise analysis and we define
1! = !^ 0 ! = xK +exx : The transfer function, T , from d(t) to e(t) is L(s)(s + ! ) T = (s + ! ) + K sL(s) and L(s) and K are chosen so that s +! T =T 2 s + 2"!s + ! 1 2 2
f
2 1
(5)
de
2
Manuscript received June 9, 2004; revised February 23, 2005 and July 28, 2005. Recommended by Associate Editor E. Bai. This work was supported by NSERC. The authors are with the Department of Electrical and Computer Engineering, The University of Western Ontario, London, ON N6A 5B9, Canada (e-mail: [email protected]). Digital Object Identifier 10.1109/TAC.2005.861712
de
2
de
bp
2
2
(6)
f
f
2
2
2
2
(7)
by matching their corresponding numerators and denominators, as in [10], where Tbp is a symmetric bandpass filter which has the center
0018-9286/$20.00 © 2006 IEEE
104
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
frequency ! and bandwidth BW and " is a small positive number which can control the width of the notch. BW must be chosen less than !c . Similar equations apply to the discrete time case. It will be shown that the variance in e, for ideal bandpass filters, is given by
e2 = (BW 0 2"!) 2 N0 :
nc (t), and ns (t) have the same variance x2 = "!N0 and mean 0. An
alternative representation [13] for the noise is
n1 = Xc (t)c(t; 0) 0 Xs (t)s(t; 0) where
Xc (t) = nc (t)cos (! 0 !c )t
0 ns (t)sin (! 0 !c)t Xs (t) = ns (t)cos (! 0 !c )t + nc (t)sin (! 0 !c )t :
(8)
1
Further, although the measurement of ! is biased in the presence of white noise, if an integral controller is used to update the nominal frequency, ! Ke ! or discrete equivalent, the frequency estimation is unbiased and its variance is
_= 1
V ar[!] =
(
0 Ke2 "!N a2
) 2
(9)
[ ]
as long as a = "!N0 > and " V ar ! =! . Note that in [1] the local exponential stability for this algorithm was shown. Throughout we assume that we remain in the region of attraction of this result. 2
A. Noise Model To simplify our mathematical expressions we will use the notation
s(t; ) = sin(wc t + ) and c(t; ) = cos(wc t + ). We use the symbol to define new variables hk , k = f1; 2; 3; 4; 5; 6g. The transfer functions for x1 , x2 from d(t) are Tdx Tdx
()
f! = 0 Tbp 2 s2 + 2K"!s + !2 K s f = Tbp 2 s2 + 2"!s + !2 :
(10)
[ ( ) ( )] = 0
()
x = "!N0 : 2
In steady state, the responses at x1 , x2 and e
(12)
(t) are
x1 = ar c(t; ' + ) + n1 x2 = ar !c s(t; ' + ) + n2 ! e = 0 ar (!2 0 !c2 )c(t; ' + ) + ne Kf !
ne =
(20)
0 K1f ! [(!2 0 !c2) Xc (t)c(t; 0) 0 Xs (t)s(t; 0) 0 2!cX_ c (t)s(t; 0) 0 2!cX_ s (t)c(t; 0) + Xc (t)c(t; 0) 0 Xs (t)s(t; 0)]:
(22)
_ () _ ()
Since Xs t , Xc t vary slowly in relation to center frequency ! , we neglect the last two terms in (22) and obtain
ne 0
1 [(!2 0 !c2)
Kf !
Xc (t)c(t; 0) 0 Xs (t)s(t; 0)
02!cX_ c (t)s(t; 0) 0 2!cX_ s (t)c(t; 0)]:
(23)
(14)
B. Mean and Variance Estimations for
1! = f (!; t) + n(!; t)
(24)
where
A(t)(!c2 0 !2 )!c2 (t; (t)) A(t)[!2 c2 (t; (t)) + !c2 s2 (t; (t))] 2!!cc(t; (t)) X_ c (t)s(t; 0) + X_ s (t)c(t; 0) n(!; t) = A(t)[!2 c2 (t; (t)) + !c2 s2 (t; (t))] f (!; t) =
(16)
where nc t and ns t are in-phase component and quadrature component, respectively, and they are uncorrelated. If n1 is zero-mean Gaussian noise, the quadrature components nc t and ns t of the noise are statistically independent zero-mean Gaussian processes. n1 ,
1!
Substituting (13)–(15), (17), (21), and (23) into (5), we get
[13]
()
!
(13)
= (=2) + 6 D(j!c ), ar = a!jD(j!c )j, and D(s) = Kf L(s)=((s2 + !2 )+ Kf sL(s)). n1 , n2 and ne are the filtered white noises at x1 , x2 , and e(t), respectively. Since Tdx , Tdx and Tde are approximately symmetric about ! , the filtered noise at x1 is a narrow band noise [12] and can be expressed as
()
()
() 2
nc (t), ns (t) and Xc (t), Xs (t) are slowly varying envelops in relation to center frequency ! [14]. Therefore, we can omit the last two terms in n2 and get n2 !c [Xc (t)s(t; 0) + Xs (t)c(t; 0)]: (21) ! The relationship between e and x1 is e = 0(1=Kf ! )(d2 x1 =dt2 + 2 ! x1 ). Therefore
where
()
()
Since the transfer function Tde is an ideal bandpass filter with a narrow notch, we can get the variance for ne given in (8) [11].
(15)
n1 = nc (t)cos(!t) 0 ns (t)sin(!t)
(19)
()
()
[ ( ) _ ( )] = 0 [ ( ) _ ( )] = 0 = (1 )_ n2 = !c [Xc (t)s(t; 0) + Xs (t)c(t; 0)] ! _ _ 0 Xc (t) c(t; 0) + Xs (t) s(t; 0)
=2
2
()
()
(11)
Here, we assume Tbp is a ideal bandpass filter with magnitude 1 and bandwidth BW. In our design [10], we get Kf "! . Therefore,when " is small enough, Tdx and Tdx will be approximately narrow first order bandpass filters with maximal magnitude 1 and bandwidth "! . The variances of noise at x1 and x2 are approximately equal and can be expressed as [11]
()
()
(18)
In this case, Xc t and Xs t will be Gaussian if nc t and ns t are Gaussian. Obviously, by the definitions for Xc t and Xs t in (18) . Thus, Xc t and Xs t and (19), we know that E Xc t Xs t are orthogonal and slowly varying in relation to !c since nc t and ns t are slowly varying in relation to !c and j! 0 !c j < BW= . We can also get E Xc t Xc t and E Xs t Xs t . Since x2 0 =! x1 , n2 can be written as follows:
!
III. NOISE ANALYSIS
(17)
(25)
(26)
and
A(t)=
ar cos(' + )+ Xc (t)
2
+ ar sin(' + )+ Xs(t)
2
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
(t)= tan01 ar sin(' + ) + Xs (t) : ar cos(' + ) + Xc (t)
()
The joint pdf for A and is (see [12, Ch. 4])
0 21 2 [A2 + a2r 0 2Aar cos( 0 (' + )] x
A : 2 2 2 x
The pdf of the phase noise, [12]
[ ( )] () ( )
Although E f !; t can be approximated as !c 0 ! by statistical and averaging methods, in the numerator of f !; t , the magnitude A t contributed by e t can not simply be cancelled by A t in the denominator of f !; t in the case of noise. We consider them different, that is, we assume that A t and =A t are independent. Therefore, we obtain E f !; t as
and
p(A; ) = exp
105
(28)
=0 1
(29)
12
()
2 SNR = a : N0
E
(30)
[n(!; t)] can E
(33)
1 A(t)
:
1) The Expectation of 1! : When SNR = a2 =N0 ! 0, the useful signal is so weak that we can not detect !c . Therefore, we use (5), (16), (21), and (23) to estimate the mean E [1! ] directly
E [1!]=E 2cos(!t + (t)) 2 (_npc (t2)sin(2!t)+n_ s (t)cos(!t)) n c + ns = 0: (31) If SNR is high, we have E [Xs2 (t)] = E [Xc2 (t)] = E [n2s (t)] E [n2c (t)] = x2 , we get (see the Appendix for details)
E [n(!; t)] = (! 0 !c ) 2 E [Xc2 (t) + Xs2 (t)] 2 E 2 1 A(t) 2 2 2 = 0 (!c 0 !) exp 0 4arx2 I0 4arx2 0 (!c 0 !) 2 h3 :
=
(34)
[1 ] = 1 0 [1 ] = 0 [1 ] = 1 [1 ] = 0 2 [1 ] = 2( ) 1 0
E [1!] = h5 0 1 h3 2 (!c 0 !) 2 2 2 2 a r = 42 exp 0 2ar2 I0 4ar2 x x x 2r 2r a a 2 I 0 4 2 + I 1 4 2 (! c 0 ! ) x x h 6 2 (! c 0 ! ) :
()
2c(t; (t)) X_ c (t)s(t; 0) + X_ s (t)c(t; 0)
2 I0 ar2 4x
While ideally, we would like E ! ! as SNR ! , it is also true that when a , ! is undefined and it is best that E ! . ! for large SNR to E ! Thus a smooth transition of E ! for no signal is desirable. When SNR ! , h5 ! = and h3 ! , 0= !c 0 ! instead and the high SNR based estimate of E ! of 0. When SNR ! 1, h5 ! and h3 ! . Thus, scaling h3 by = in (34) will not affect the calculation for high SNR and cause the result to be consistent with the low SNR result. Thus
where I0 is the modified zeroth order Bessel function. Formula (29) is called Ricean distribution [13]. Here, we assume A t and t are uncorrelated for all input signal-to-noise ratio (SNR)
E [n(!; t)]
2
E [1!] = E [f (!; t)] + E [n(!; t)] = (h5 0 h3 ) 2 (!c 0 !):
0 21 2 (A2 + a2r ) I0 A 22ar x x
Therefore, by using statistical and averaging methods, E be approximated as
1
A
Therefore
plus narrow band Gaussian noise is
()
()
1 ()
= (!c 0 !) 2 exp 0 2ar2 x 2r 2r 2 2 a a 2 1 + 22 I0 42 + 2ar2 I1 4ar2 x x x x (! c 0 ! ) hs :
(27)
where the error function (erf ) is defined as erf (x) = (2=p) 0x exp(0z2 )dz . The pdf of the envelope A(t) of a sine wave 0
pA (A) = A2 exp x
() [ ( )]
E [f (!; t)] = (!c 0 !) 2 E [A] 2 E
1 = (t) 0 (' + ), can be obtained by
2 cos(1) p1 (1) = 1 exp 0 ar2 + ar p 2 2x 2 2x 2r 2 a 2 exp 0 22 sin (1) x cos(1) 2 1 + erf ar p 2x2
( )
1! :
2) The Variance of The variance of ! is
1
Next, we calculate the variances of
V ar[1!] = E [(1! 0 E [1!])2 ]:
=
(35)
1! . (36)
When ! 6 !c , although the average value of the (25) is a constant, there is a high frequency oscillation component ( !c ) in f !; t . This high frequency component can be obtained by
2
( )
2 2 f (!; t) 0 (!c 0 !) !c 0 ! cos2 (!c t + (t)) 0 (!c 0 !) !c (37) (!c 0 !)cos 2(!ct + (t)) : Therefore
V ar[1!] = E [f 2 (!; t)] 0 2E [f (!; t)]E [n(!; t)] + E [n2(!; t)] 0 E 2 [1!] = 32 h12h24 0 2h3 h5 0 h62 (!c 0 !)2 + 2("!)2 e2h2 (38)
(32)
where h1 , h2 , h4 are defined in the appendix.
106
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
Fig. 1. Block diagram of an internal model control system.
D. Comparison With Cramer–Rao Bound
C. Noise Analysis for Integrated Error The simplest means of updating ! is to use an integral controller. !c The integral controller will converge to the correct answer, ! and ar a, as long as the sign of h5 0 h3 is positive, that is
=
=
a2r > max x2 I0
2I0
a
4
= 2:
+ I1 This inequality means that the power a2 =2 (x21 + x22 )=2 should be 2 a 4
a 4
greater than the variance x in the states of IM in order to make the algorithm work. Since the integral operator has the function of a low pass filter, the high frequency components in the estimation of ! can be neglected. Thus, by using the result in (45) in the Appendix , the noise in the output of the integral controller is
1
X_ s ( )cos( ( )) 0 X_ c ( )sin( ( )) d !c 0 ! = Ke A( ) 01
The Cramer–Rao bound (CRB) is a lower bound on the error variance of any unbiased estimate. The CRB provides a standard against which to measure estimation algorithm performance. We are interested in comparing the behavior of the estimated frequency produced by our integral controller, which behaves like ! 0Ke ! Ke !c and Ke = . In time domain, it is the convolution of Ke e0K t and !c . In order to compare with our method, we calculate the CRB for a sinusoid signal with a time-varying SNR with !c 2 ; 1 and Kn e t 2 . The choice of this SNR results in a time varying weighting of the signal that will result in tracking behavior identical to our algorithm, at the expense of including more noise. Thus, it is possible to achieve better results than the presented CRB. If a white noise has an autocorrelation N0 , the Cramer–Rao bound in this case is [12] function R
=
()
()
where Ke is an adaptation gain. Since SNR is high, A t and t are slowly varying functions and we can approximate them as constants E =A t and E t . So
[1 ( )] [ ( )] !c 0 ! = Ke [Xs (t)cos(E [ (t)]) 0Xc (t)sin(E [ (t)])]E A1(t)
( )=
:
(40)
+ E [Xc2(t)]sin2 (E [ (t)])] 2 E 2 A1 a2r 4x2 :
(41) 2 p For high SNR, since I0 (x) ex = 2x and x2 = "!c N0 in (41), we
have the result (9). When we reduce Ke , ", N0 or increase a, we can reduce the output variance. This variance is independent of the bandwidth BW of Tde , though BW < ! for solution to controller equations to be feasible. =Ke . ThereThe convergence time constant for the algorithm is tc fore, the larger Ke , the smaller the convergence time tc . There exists a tradeoff between fast tracking speed and small noise in frequency estimation. From the previous analysis, note that we just require that the noise has constant spectral density N0 in the passband BW of the closed-loop transfer function Tde . Therefore, the noise used in this note could be a band limited white noise, with a constant component.
=1
(42)
In this section, we show some results of simulations to verify our noise analysis. In other words, we show the various approximations we have made have minimal effect on the accuracy of (34), (35), (38), and (41). First, we select a second-order Chebyshev low-pass filter and obtain a desired bandpass filter with a notch 2 2 Hbpn = 1:102510BW s Q (s )
V ar[!c 0 !] = Ke2 [E [Xs2 (t)]cos2 (E [ (t)])
2
2 !2 4 2N0 : a
)
The previous formula states that smaller will give a smaller Cramer-Rao bound. In the next section, we will show the comparison of the variance of our method with CRB.
integral controller is
2 0 2ar2 I02 x
()
(
IV. SIMULATIONS
Since Xs and Xc are Gaussian random with mean zero and variance x2 , the mean of !c 0 !o is zero. The variance of the output of the
2 = Ke exp
+
(0 )
t
(39)
_=
2
+! 2 s2 +s2"!s + !2 2
2
(43)
where
Q(s) = s4 + 1:097734BW s3 + (2!2 + 1:102510BW 2 )s2 +1:097734!2BWs + !4 : By using the designing method for the desired bandpass filter with a b1 s2 b2 s b3 =s2 , notch for Tde in [10], we get C s 2 2 Gs : BW = s a1 s a2 , L s Gs= G s C s , and an IM Kf s= s2 !2 , where a1 : BW "! , a2 : BW 2 "! !2 , 3 2 b1 , b2 "! : ! BW = : BW 2 , 4 2 b3 ! = : BW , and Kf "! . Since our bandpass filter design is based on a second-order Chebyshev low-pass filter with p and is far from a ideal bandpass filter, we approximate (see [11])
() = ( + + ) ( ) = 1 102510 ( + + ) () = ( + ) ( ) (1 + ( ) ( )) = 1 097734 + 2 = 1 097734 2 +2 = 1 = (2 + 1 097734 ) 1 102510 = 1 102510 =2 = 1 dB e2 = (BW 0 2"!) 2 N0 :
2
(44)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
Fig. 2.
Theoretical and simulated means versus the a = for a = 1, ! = 2
Fig. 3. Theoretical and simulated variance versus the a = , a = 1, ! = 2
107
2 490, !
2 490, !
= 2 2 500, BW = 300, " = 0:01.
= 2 2 500, BW = 300, " = 0:01, = 10
.
108
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
= for the output of integral controller with a = 1, ! = 2 2 100, K = 5, BW = 200, " = 0:02.
Fig. 4.
Theoretical and simulated variances versus the a
Fig. 5.
Comparison of the CRB and the variance of the proposed frequency estimation method.
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
The simulated variance was based on a window of data which has 50 s width. Figs. 2 and 3 show the simulated mean and variance versus the a2r =x2 and the corresponding theoretical mean and variance for the initial value ! = 2 2 490. We can see that the theoretical results (34), (35) and (38) tally with the simulated results for a2r =x2 2. The modified formula (35) is better than (34). Fig. 4 shows the relationship between the variance of the output of the integral controller and a2 =x2 . A discrete implementation of the algorithm given by (3) was implemented with !c = 2 2 100, Ke = 5, BW = 200, " = 0:01, signal to noise ratio a2 =No = 10 000 and sampling rate of 0.5 ms. G(s) and C(s) were converted to discrete representations via bilinear transformation. The measured variance of estimated frequency was 0.0162 rad=s2 versus the 0.0157 rad=s2 given by (9). The comparison of the CRB and the variances of our frequency estimation method with " = 0:01, 0.001 and !c = 2 2 100 is shown in Fig. 5. From Fig. 5, we know that if we reduce the value of ", the variance of the estimated frequency will approach to and then exceed the corresponding CRB for no a priori knowledge of ! . But for smaller ", the region of convergence will be reduced. In particular, our algorithm remains well behaved only for "! > V ar[!]. Therefore, there exists a tradeoff between small estimate variance and big region of convergence. We have extended our approach to the cases of harmonics and multitone signal in [15] and [1], respectively, and shown the corresponding simulated results. In [1], we also showed the performance of our algorithm to tracking a pseudo periodic signal with time-varying amplitude and frequency. In [10], [15], and [16], we compared our algorithm to the methods proposed in [5], [6], and [9]. It is found that our algorithm tracks time-varying frequency with less delay and less susceptible to noise.
109
1 E 2c(t; (t)) 2 X_ c (t)s(t; 0) + X_ s (t)c(t; 0) E A(t) = E[X_ c (t) (sin(2!c t + (t)) 0 sin( (t))) + X_ s (t) (cos(2!c t + (t)) + cos( (t)))]E
APPENDIX A. Approximation of E[n(!; t)] Since
1 TE !!c dt=1: T 0 !2 cos2 (!c t + (t)) + !c2 sin2 (!c t + (t)) E [n(!; t)] can be approximated as E[n(!; t)]
(45)
We neglect the high frequency components in the previous equation by averaging and get
E [n(!; t)] = (! 0 !c ) 2 E Xc2 (t) + Xs2 (t)
1 : 2 E 2 A(t)
(46)
B. Expectations of 1=A and 1=A2 In the general case, the v th moments of a Rice pdf can be analytically expressed as a function of the confluent hypergeometric (Kummer’s) function of the first kind 1 F1 [13]
E [Av ] =
1
0
Av pA (A)dA
2 = (2x2 )v=2 0 1 + v 1 F1 0 v ; 1; 0 ar2 (47) 2 2 2x 1 where 0(c) = 0 c01 exp(0 )d is Gamma function. From for-
mula (47), we can obtain
2 E [A] = x exp 0 ar2 2 4x 2r a a2r + a2r I a2r 2 1 + 22 I0 4 2x2 1 4x2 x x2 h4 2 E A = a2r + 2x2
V. CONCLUSION We show that although the measurement 1! in our algorithm is biased, the estimate ! of the update law is not biased and a formula for its variance is given. Even with constant offsets in d, because of the presence of the bandpass filter Tbp , signals n1 , n2 and e are guaranteed to be zero mean, ensuring our algorithm, unlike many other algorithms, is unaffected. A necessary condition to make our frequency estimation algorithm work appropriately is that the signal power must be greater than noise power at the states of the IM. When the values of the designed parameters " and Ke are reduced, the variance of the estimated frequency can be decreased. But the convergence rate of the identified algorithm declines with the reduction of " and Ke and region of rapid convergence decreases with ". The noise model presented here will allow optimal selection of gains for this algorithm. The results of this note can be applied directly to feedback problem when plant is low pass filter and disturbance is bandlimited white noise. When the disturbance is band limited, the Fourier transform of e(t) will be bandpass with notch just as we require.
1 : A(t)
(48) (49)
where I1 () is the first order modified Bessel function. When v = 01 and using the integral representation [17, eq. (13.2.1)] for 1 F1 , integral formula [18, eq. (3.364)], we have
2 E 1 = 2x2 01=2 0 1 1 F1 1 ; 1; 0 ar2 A 2 2 2x 2r 2r a a = exp 0 2 I0 2x2 4x 4x2 h1 :
(50)
In practice, we add a small positive constant real number in the denominator of 1! . For purpose of analysis, we calculate E 1=(A2 + ) . Considering very small (1008 in our application) and using the [18, eqs. (3.351), (3.352), and (8.214)], we have
1 A2 + 1 1 = pA (A)dA 0 A2 + exp 0 2a = 2x2 a2r 2 exp 2 2 E1 2 2 + Ei 2 x x x2 h2
E
1
a2r 0 C 0 ln 2 2 x
where E1 (z) = z (exp(0t)=t)dt and Ei (z) 0 lim"!+0 [ 00z" (exp(0t)=t)dt + "1 (exp(0t)=t)dt]
(51)
= are
110
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 1, JANUARY 2006
exponential integrals with z > 0, C = 0:577 215 664 9 1 1 1 is Euler’s constant. The formula (51) is finite for > 0.
Transition-Time Optimization for Switched-Mode Dynamical Systems
REFERENCES
Magnus Egerstedt, Yorai Wardi, and Henrik Axelsson
[1] L. Brown and Q. Zhang, “Identification of periodic signals with uncertain frequency,” IEEE Trans. Signal Process., vol. 51, no. 6, pp. 1538–1545, Jun. 2003. [2] P. A. Regalia, “An improved lattice-based adaptive IIR notch filter,” IEEE Trans. Signal Process., vol. 39, no. 9, pp. 2124–2128, Sep. 1991. [3] L. Hsu, R. Ortega, and G. Damm, “A globally convergent frequency estimator,” IEEE Trans. Autom. Control, vol. 44, no. 4, pp. 698–713, Apr. 1999. [4] S. Bittanti, M. Campi, and S. Savaresi, “Unbiased estimation of a sinusoid in colored noise via adapted notch filters,” Automatica, vol. 33, no. 2, pp. 209–215, 1997. [5] H. Kowok and D. Jones, “Improved instantaneous frequency estimation using an adaptive short-time Fourier transform,” IEEE Trans. Signal Process., vol. 48, no. 10, pp. 2964–2972, Oct. 2000. [6] B. Wu and M. Bodson, “A magnitude/phase-locked loop approach to parameter estimation of periodic signals,” IEEE Trans. Signal Process., vol. 48, no. 10, pp. 2964–2972, Oct. 2000. [7] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood and Cramer-Rao bound,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 5, pp. 720–741, May 1989. [8] B. F. Scala and R. Bitmead, “Design of an extended Kalman filter frequency tracker,” IEEE Trans. Signal Process., vol. 44, no. 3, pp. 739–742, Mar. 1996. [9] T. Tsao, Y. Qian, and M. Nemani, “Repetitive control for asymptotic tracking of periodic signals with an unknown period,” J. Dyna. Syst., Measure., Control, vol. 122, pp. 364–369, Jun. 2000. [10] Q. Zhang and L. J. Brown, “Designing of adaptive bandpass filter with adjustable notch for frequency demodulation,” in Proc. 2003 Amer. Control Conf., vol. 4, Jun. 4–6, 2003, pp. 2931–2936. [11] R. G. Brown, Introduction to Random Signal Analysis and Kalman Filtering. New York: Wiley, c1983. [12] A. D. Whalen, Detection of Signals in Noise. New York: Academic, 1971. [13] N. Wax and Ludeman, Selected Papers on Noise and Stochastic Processes. New York: Dover, 1954. [14] R. L. Stratonovich, Topics in the Theory of Random Noise. New York: Gordon and Breach, pp. 1963–1967. [15] L. J. Brown and Q. Zhang, “Periodic disturbance cancellation with uncertain frequency,” Automatica, vol. 40, no. 4, pp. 631–637, Apr. 2004. [16] Z. Zhao and L. J. Brown, “Fast estimation of power system frequency using adaptive internal model control technique,” in Proc. 43rd IEEE Conf. Decision and Control, Dec. 14–17, 2004, pp. 845–850. [17] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables. Washington, DC: U.S. Government, 1970. [18] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products. New York: Academic, 1980.
Abstract—This note considers the problem of determining optimal switching times at which mode transitions should occur in multimodal, hybrid systems. It derives a simple formula for the gradient of the cost functional with respect to the switching times, and uses it in a gradient-descent algorithm. Much of the analysis is carried out in the setting of optimization problems involving fixed switching-mode sequences, but a possible extension is pointed out for the case where the switching-mode sequence is a part of the variable. Numerical examples testify to the viability of the proposed approach. Index Terms—Gradient-descent algorithms, hybrid systems, optimal control, switched dynamical systems, switching-time control.
I. INTRODUCTION Switched dynamical systems are often described by differential inclusions of the form
x_ (t) 2 fg (x(t); u(t))g 2A
(1)
where x(t) 2 n , u(t) 2 k , and fg : n+k ! n g 2A is a collection of continuously differentiable functions, parameterized by belonging to some given set A. The time t is confined to a given finite-length interval [0; T ]. Such systems arise in a variety of applications, including situations where a control module has to switch its attention among a number of subsystems [12], [15], [19], or collect data sequentially from a number of sensory sources [4], [6], [11]. A supervisory controller is normally engaged for dictating the switching law, i.e., the rule for switching among the functions g in the right-hand side of (1). Recently, there has been a growing interest in optimal switching time control of such hybrid systems, where the control variable consists of a proper switching law as well as the input function u(t) (see [3], [5], [9], [10], [16]–[18], and [20]). Reference [3] establishes a framework for optimal control, and [16]–[18] present suitable variants of the maximum principle to the setting of hybrid systems. [1], [2], [9], [14] consider piecewise-linear or affine systems. The special case of autonomous systems, where the term u(t) is absent and the control variable consists solely of the switching times, is considered in [9], [11], [21], [22]. In particular, Xu and Antsaklis [21], [22] consider general nonlinear systems, and they have developed nonlinear-programming algorithms that compute the gradient and second-order derivatives of the cost functional. This note, whose preliminary version has appeared in [8], also falls in this category. [21] provides the starting point for the results presented in this note, as we initially consider a similar problem, where the sequence of switching functions as well as the number of switching times are fixed. We develop a simpler formula than the one in [21] for the gradient of the cost functional, and use it in a gradient-descent algorithm. Finally, we suggest a possible extension of the algorithm to a class of problems Manuscript received November 11, 2003; revised December 15, 2004 and June 14, 2005. Recommended by Associate Editor J. Hespanha. This work was supported in part by the National Science Foundation under Grants 0207411 and 0509064, and by a grant from the Georgia Institute of Technology Manufacturing Research Center. The authors are with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332 USA (e-mail: [email protected]; [email protected]; [email protected]). Digital Object Identifier 10.1109/TAC.2005.861711
0018-9286/$20.00 © 2006 IEEE
