Analysis of a Phase Method for Time-Delay Estimation. Short Version.

Analysis of a Phase Method for Time-Delay Estimation.

Short Version.

Svante Bj¨orklund

Control & Communication

Department of Electrical Engineering

Link¨opings universitet, SE-581 83 Link¨oping, Sweden

WWW: http://www.control.isy.liu.se

E-mail: svabj@isy.liu.se

1st September 2003

AUTOMATIC CONTROL

COMM_{UNICATION SYSTE}MS

LINKÖPING

Report no.: LiTH-ISY-R-2539

Technical reports from the Control & Communication group in Link¨oping are available

at http://www.control.isy.liu.se/publications.

Sammanfattning Abstract Nyckelord Keywords Rapporttyp Report: category Licentiatavhandling C-uppsats D-uppsats Övrig rapport Språk Language Svenska/Swedish Engelska/English ISBN

Serietitel och serienummer

Title of series, numbering

URL för elektronisk version

Titel Title Författare Author Datum Date Avdelning, Institution Division, department

Automatic Control

X

1400-3902

Department of Electrical Engineering

2539

Analysis of a Phase Method for Time-Delay Estimation. Short Version.

In this report, estimation of time delays in linear dynamical systems with additive noise is

stud-ied. The report reviews and analyses a certain estimation method: A discrete-time model of a

continuous-time system is identified. The non-minimum phase zeros form the allpass part of the

model. The dead-time is estimated by the slope at low frequencies of the phase of the allpass

part. The method has been experimentally studied with the aid of simulations in open and closed

loop. The results show that the estimation method can be used with several different model

struc-tures. The method is, however, non-robust and can totally fail in some cases. The probability of

failure depends on the model structure, the input signal type and the signal-to-noise ratio. The

estimation method can be made more robust by removing false non-minimum phase zeros,

caused by the noise, in a certain way. It is usually better not to prewhite the data applied to the

method.

time-delay, dead-time, estimation, system identification, Laguerre, linear dynamic systems,

non-minimum phase, zero, ANOVA, confidence intervals, simulations, open loop

X

2003-09-01

ANALYSIS OF A PHASE METHOD FOR TIME-DELAY ESTIMATION

Svante Bj¨orklund

Control & Communication Department of Electrical Engineering

Link¨opings universitet, SE-581 83 Link¨oping, Sweden

WWW: http://www.control.isy.liu.se E-mail: svabj@isy.liu.se

Abstract:In this paper, estimation of time delays in linear dynamical systems with additive noise is studied. The paper reviews and analyzes a certain estimation method: A discrete-time model of a continuous-time system is identified. The non-minimum phase zeros form the allpass part of the model. The dead-time is estimated by the slope at low frequencies of the phase of the allpass part. The method has been experimentally studied with the aid of simulations in open and closed loop. The results show that the estimation method can be used with several different model structures. The method is, however, non-robust and can totally fail in some cases. The probability of failure depends on the model structure, the input signal type and the signal-to-noise ratio. The estimation method can be made more robust by removing false non-minimum phase zeros, caused by the noise, in a certain way. It is usually better not to prewhite the data applied to the method.

Keywords:

time-delay, linear systems, estimation, identification, non-minimum phase

1. INTRODUCTION

In this paper, estimation of time delays (or dead-times) in linear dynamical systems with additive noise is addressed. This is a common engineering problem, e.g. in control performance monitoring of industrial processes, in design and tuning of controllers, in range estimation in radar and in direction finding in signal intelligence. It is also a necessary part in all system identification. In (Isaksson 1997) a method for dead-time es-timation is described. A discrete-time model of a continuous-time system is identified. The ze-ros of the model are translated to continuous-time. By comparing the dead-time with a Pad´e approximation, the dead-time may be estimated

from the continuous-time non-minimum phase ze-ros. An improved method is described in (Horch 2000, Isaksson et al. 2001) in which the discrete-time non-minimum phase zeros form the allpass part of the model, which directly represents the dead-time. The dead-time is estimated by study-ing the slope at low frequencies of the phase of the allpass part. In (Isaksson 1997) the discrete-time model is identified using general orthonormal bases. In (Horch 2000, Isaksson et al. 2001) a particular such basis is used, the Laguerre basis. The method in (Horch 2000, Isaksson et al. 2001) shows very good results in both open loop and closed loop for input and reference signals in the form of steps. This method will in this paper be called Laguerre DAP . The portion of the method from the discrete-time model to the dead-time

estimate is here called the DAP (Discrete-time Allpass part Phase) method.

The purpose with this work is to review the Laguerre DAP method, describe why this method appears to work so well and investigate if the method has some drawbacks. A more thorough description of the work underlying this paper is given in (Bj¨orklund 2002).

2. DEAD-TIME ESTIMATION BY THE PHASE OF THE ALLPASS PART 2.1 The DAP method

Assume the true continuous-time system is ¯G(s) =

¯

G1(s)_·e−sTd_{= ¯G1(s)· ¯Gap(s). The system ¯G1(s) is}

linear and time invariant. A discrete-time rational

linear model G(z) of ¯G(s) is estimated and

fac-torized into a minimum-phase system G1(z) and an allpass system Gap(z): G(z) = G1(z)Gap(z).

Then, Gap(eiωTs_{) is considered as an}

approxima-tion of the dead-time ¯Gap(iω) = e−iωTd _(Horch

2000, Isaksson et al. 2001). The allpass part

Gap(z) of G(z) is formed by the non-minimum phase (outside the unit circle) zeros of G(z) and poles added to Gapt(z) which are these zeros mir-rored in the unit circle.

The dead-time estimate ˆT is given by an

approx-imation of the derivative of the phase ϕ(ω) =

arg Gap(eiωTs_{) and addition of a 1 because of the}

extra time delay that is created by the sampling: ˆ

Td =₋arg Gap(e

iω1Ts₎

ω1Ts + 1; ω1 small (1)

This method is described in (Horch 2000, Isaksson

et al. 2001), where ω1 = 10−4 _{is suggested. This}

ω1 is also used in this paper. Other forms of

approximations of the derivative should also be possible.

In (Horch 2000, Isaksson et al. 2001) a Laguerre model (Wahlberg 1991) was used. Other linear model structures, for example FIR, ARX or out-put error (OE) model structures (Ljung 1999) can also be used.

2.2 Examples of using the DAP method

2.2.1. A successful estimation Now, an example

of when the DAP method works well is given. The

linear continuous-time system G2in Section 3.1.1

was simulated in Matlab by the function lsim

Real Axis Imag Axis Pole−zero map 0.5 1 1.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6

0.8 Pzmap whole model, t122b1, Trial 3

Figure 1. Pole-zero plot of the identified Laguerre model for the subsequent successful

dead-time estimation ( ˆTd = 9.7114). The setup

is as in Section 3.1.1: The system G2 was

simulated with the input signal RBS 10-30%. The signal-to-noise ratio (SNR) was 10. The Laguerre model was identified from the input-output data. No prewhitening or zero guard-ing was used (see section 3.1.1.1). The poles of the model should all be located at α = 0.8 but due to well-known numerical problems with multiple poles they are somewhat spread in the figure.

in (CSTB n.d.) with the input signal “RBS 10-30%” in Section 3.1.1. To the output signal, white Gaussian noise was added and the signal-to-noise ratio (SNR) was 10. The sampling interval was

Ts = 1. A Laguerre model was identified from

the input-output data and the DAP method was used to estimate the dead-time. The simulation resulted in the poles and zeros of the identified Laguerre model depicted in Figure 1. The dead-time estimate became 9.7714 which is a good estimate since the true dead-time is 10.

2.2.2. A failing estimation We will now study

a case when the dead-time estimation failed. The only difference to the successful trial in the setup is a different noise realization. The zeros and poles of the identified Laguerre model are shown in Fig-ure 2. We note that the zero on the real axis just inside the unit circle in the successful simulation has moved to just outside the unit circle. The dead-time estimate will now be 1321.86, which is completely wrong. If we remove the moved zero, the dead-time estimate will be 11.6586, which is an acceptable estimate.

2.3 A solution

The reason for zeros falling on the incorrect side of the unit circle is the noise. Figure 3 shows an

Real Axis Imag Axis Pole−zero map 0.5 1 1.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6

0.8 Pzmap whole model, t122a1, Trial 941, Sys 2

Figure 2. Pole-zero plot of the identified Laguerre model for a subsequent failing dead-time

esti-mation ( ˆTd = 1321.86). The simulation setup

is the same as in Figure 1. This simulation was one of only 3 out of 1024 simulations with failing dead-time estimation for SNR = 10. With a lower SNR the percentage of failing estimations is much higher, see Figure 3. example of the spread of zeros and poles due to the noise. Figure 4 displays the dead-time estimate for different locations of the zero closest to +1.

0.5 1 1.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6

0.8 Pzmap whole model, 3 std, t125a1, Trial 3

To y1

Figure 3. Pole-zero plot of a Laguerre model with estimated uncertainty regions (3 standard deviations) and with zeros and poles from 1024 simulated trials for SNR = 1. As can be seen, the risk of a zero falling on the wrong side of the unit circle is significant. The simulation setup is as in Section 3.1.1 with

system G2 in open loop, signal RBS 10-30%,

no prewhitening and no zero guarding(see section 3.1.1.1).

It appears that moving zeros located close to but outside the unit circle (back) to the inside of the unit circle is a solution to the problem with the DAP method. The motivation is that we assume that these zeros actually should be located inside the unit circle. Since we only need the allpass part in the DAP methods, we just remove some zeros outside the unit circle without putting them

0.999 0.9992 0.9994 0.9996 0.9998 1 1.0002 1.0004 1.0006 1.0008 1.001 −5 0 5 x 10−6 2983 15 24393 10 13 15 13 14 10 12692822356794449 254923192020 t123b1:Zeros closest +1, 10−30%, SNR=1.

Figure 4. Dead-time estimate for different loca-tions of the zero closest to +1 for a Laguerre model.

somewhere inside the unit circle. We call this technique for zero guarding.

3. STATISTICAL ANALYSIS

Now follows a statistical analysis of the DAP method and investigation of the choice of zero guarding by some simulations done in Matlab.

3.1 Open loop simulations

3.1.1. Open loop setup In the open loop

simula-tions the output signal y(t) was simulated as

y(t) = G(s)u(t) + v(t) (2)

where y(t), u(t) and v(t) are the output, input and noise signals, respectively. G(s) is a linear system with dead-time. The noise v(t) was white and Gaussian. The system G(s) was simulated by the function lsim in (CSTB n.d.) with the sampling

interval Ts= 1.

3.1.1.1. Varied factors. The varied factors are

given the names Method, Prewhite, Sys, ZType, ZNo, ZSize, InType and SNR. The different pos-sible choices for the same factor are called levels. The choice of level for all factors is called a fac-tor level combination. Here follows a list with all factors and their levels. The names of the factors and factor levels are used in the following plots. Method. Dead-time estimation method:

• Lag.dap. The DAP method with a Laguerre model (Wahlberg 1991) with 10 terms in the sum and the fix pole α = 0.8.

• FIRdap. The DAP method applied to a FIR model with 15 taps.

• ARXdap. The DAP method with an ARX model y(t) = (B(q)/A(q))u(t)+(1/A(q))w(t)

with na = 4 (the order of the A(q)

poly-nomial), nb = 15 (the order of the B(q)

• OEdap. The DAP method applied to an OE model y(t) = (B(q)/F (q))u(t) + w(t) with

nb= 15 (order of the B(q) polynomial), nf =

2 (order of F (q) polynomial) and nk= 0 (the

number of time delays).

Prewhite. The input and output signals were fil-tered through a filter that made the input white (pw ) or not filtered (nopw ).

It was possible to adjust the DAP methods by zero guarding. There are three factors, called ZType, ZSize and ZNo related to zero guarding:

• ZType. Remove non-minimum-phase zeros close to the point +1 or close to the unit circle (ucirc) or remove no zeros.

• ZSize. How far away from +1 or the unit circle to remove non-minimum-phase zeros. • ZNo. The maximum number of zeros to

re-move within ZSize. The zeros closest to +1 or the unit circle are removed first.

Sys. All systems were of the form Gj(s) =

e−9s_·G¯_{j(s). ¯G1had poles−0.1 & -1, no zeros and}

DC gain 1 (a slow second order system). ¯G2 had

poles _{−1 & -10, no zeros and DC gain 1 (a fast}

second order system). ¯G5 had poles _{−0.1, −0.3,}

−0.6 & −1, zeros −0.4 & −0.9 and DC gain 1

(a fourth order system with real poles). ¯G6 had

poles−0.1(1 ± i)/√2 & −(1 ± i)/√2, zeros −0.4

&−0.9 and DC gain 1 (a fourth order system with

complex poles). For all the systems the time delay was 10 after the sampling.

InType. The input signal was 500 samples long and could be of three different types:

• RBS 10-30%. (Random Binary Signal) with frequency contents between 10%-30% of the Nyquist frequency. It was generated by the function idinput in (SITB n.d.).

• RBS 0-100%. RBS with frequency contents between 0%-100% of the Nyquist frequency, i.e white noise. It was generated by the func-tion idinput in (SITB n.d.).

• Steps. Three steps of the form (in Matlab code): [zeros(50,1);ones(150,1); -ones(150,1); zeros(150,1)].

SNR. The used SNR was in the range 1 to 100.

3.1.2. Simulation results In Figure 5 the RMS

error of the DAP dead-time estimate is displayed for several model structures (Method ). The SNR

was high (SNR_{≈ 100) or low (SNR ≈ 1). The four}

systems (Sys) in Section 3.1.1.1 were employed. For each level combination of the factors Method, Prewhite, InType and Sys the RMS error was

estimated from 1024 trials. Then for each level combination of the other factors, the system with the worst, i.e. highest, RMS error (not shown here), was plotted.

5 100*10−30% 100*0−100% 100*steps 1*10−30% 1*0−100% 1*steps FIRdap.*nopw FIRdap.*pw ARXdap.*nopw ARXdap.*pw OEdap.*nopw OEdap.*pw Lagu.dap*nopw Lagu.dap*pw 0 1000 2000 6.75 30.1 992 892 365 6.86 5.71 698 85.3 387 Method*Prewhite 1.2 3.78 473 786 274 227 3.24 2.94 481 328 426 6.7 2.07 2.66 655 305 6.99 2.84 2.64 151 95.9 962 717 1.43 343 263 1.34e+03 246 169 332 1.3 861 285 1.1. MIN 2.45 SNR*InType 1.52e+03. MAX 1.55 1.53

Figure 5. RMS error for some dead-time estima-tion methods in open loop. No zero guard-ing. The axis label “SNR*InType” means that on this axis there are different (factor level) combinations of SNR and input signal type. In the same way “Method*Prewhite” means different combinations of estimation method and prewhitening. The tick mark labels tell us what factor level combinations there are in the different “rows” and “columns”. The level “10-30%” is an abbreviation of “RBS 10-30% and “0-100%” of “RBS 0-100%”. Axis labels and tick mark labels are further explained in Section 3.1.1.

As can be seen in Figure 5, some of the DAP methods exhibit very large RMS errors for some combinations of InType, Prewhite and SNR. The methods totally fail for these factor level combina-tions. Some estimates in these cases are very large, up to about 30000 (not shown here). Laguerre DAP without prewhitening is the best method for the input signal type Steps at high SNR. It is also among the best methods for low SNR. The method OEdap fails in all cases. The only system that could be estimated without failure was the

fast second order system G2 and then only with

FIRdap and ARXdap. As a summary of Figure 5 we can say that the DAP methods seem to be non-robust and they can totally fail in some cases. In order to mitigate the failures we have tried zero guarding on the Laguerre model structure with the following distances (ZSize) to be considered as “close”: 0.05, 0.10, 0.15, 0.20, 0.25, 0.25 and 0.30. We have tested different maximum number (ZNo) of zeros to remove: 1 to 4.

We have utilized the statistical method ANOVA (Analysis of Variance) and confidence intervals (see (Montgomery 1997)) to discover significant differences between the different combinations of the levels of ZType, ZSize and ZNo. See (Bj¨orklund 2002) for details.

The result showed that there were several “good” combinations with no significant differences (over-lapping confidence intervals). It was, however, clear that removing zeros close to +1 is bet-ter than removing zeros close to the unit circle (ucirc). It was also clear that removing more than one zero is often necessary. We chose one of the good combinations: ZType=+1, ZSize=0.15 and ZNo=3 for the next simulation. This is some form of statistical optimization of an estimation method.

Figure 6 shows the RMS error for some estimation methods with the zero guarding. No DAP meth-ods fail but give reasonable and low RMS errors. Even OEdap works and give good estimates. Still Laguerre DAP gives the best result for Steps at high SNR. At low SNR it is among the best. The

fast second order system G2 gave the best result

in average (over input signal type and SNR) for all methods (model structure and with/without

prewhitening). The fourth order system G4 with

complex poles gave the worst result for all meth-ods. 100*10−30% 100*0−100% 100*steps 1*10−30% 1*0−100% 1*steps FIRdap.*nopw FIRdap.*pw ARXdap.*nopw ARXdap.*pw OEdap.*nopw OEdap.*pw Lagu.dap*nopw Lagu.dap*pw 0 5 10 6.17 5.07 6.63 5.84 5.6 6.11 5.71 5.69 5.06 6.11 Method*Prewhite 1.2 3.34 5.78 5.53 6.75. MAX 6 3.24 2.94 5.7 5.08 6.52 5.34 2.07 2.59 5.58 5.48 6.31 2.69 2.63 3.41 5.63 4.89 6.6 1.43 1.6 3.91 5.41 5.03 1.67 1.69 1.3 5.86 1.98 1.1. MIN 2.45 SNR*InType 1.66 1.55 1.53

Figure 6. RMS error in open loop of different dead-time estimation methods with the zero guard-ing ZType=+1, ZNo=3 and ZSize=0.15. 1024 trials. Otherwise the setup is as in Fig-ure 5. Compare with FigFig-ure 5.

3.2 Closed loop simulations

In automatic control of systems and processes, usually feedback is used, resulting in closed loop systems. We have also studied dead-time estima-tion by DAP methods in closed loop by statistical analysis. The setup and the results are described in (Bj¨orklund 2002). The results are similar to the ones for open loop. Laguerre DAP with zero guarding but without prewhitening and using the system input as the input for the estimation is the best method for step signals.

4. DISCUSSION

The pole location of the Laguerre model probably affects the locations of the zeros and accordingly also how sensitive they are for falling on the wrong side of the unit circle due to noise.

As seen in this paper, the DAP method fails if a zero erroneously falls outside the unit circle. The error in the dead-time estimate is not proportional to the position error of the zero. A small position error which does not cause the zero to fall on the wrong side of the unit circle will result in a small estimation error (Figure 4). But when the zero falls on the incorrect side of the unit circle, the estimation error will be large. The estimation error is larger close to the point +1 and smaller further away.

If we move a zero from inside to outside of the

unit circle when ω = 0, then arg G(eiωTs_{) will}

increase by π. If ω is not zero but very small, the phase increase will be between 0 and π. When this extra phase is divided by a very small ω1 in equation (1), the dead-time estimate will get a large bias that in most cases will dominate the estimate. The maximum possible estimation error (for a single incorrect zero) occurs for the

maximum phase error in arg Gap(eiω1Ts_{), which is}

π. The estimation error will then be ˜Td = ˆTd₋

Td= π/ω1= 3.1416/10−4_{= 31416.}

An advantage of the DAP method compared to many other methods is that it can estimate time-delays that are a fraction of a sampling interval, which is necessary in many applications.

The zero guarding chosen in Section 3.1.2 for Laguerre in open loop seems to work also for closed loop. It also appears to function for the other model structures. Thus, the choice of zero guarding appears to be robust.

In Section 3.1.2 we saw that sometimes it is not enough to remove only one zero. The reason is

likely that because complex zeros come in complex conjugated pairs (for real valued systems), both zeros in the pair must be removed.

The location of the true non-minimum phase zeros probably depends on the dead-time and the sampling interval. The longer dead-time or the shorter sampling interval, the closer to the point +1 will the true zeros be. Therefore, the choice of ZSize depends on the maximum possible dead-time and on the sampling interval.

A reason to why Laguerre DAP works well for step-like inputs could be as follows. An estimated model will become better in frequency ranges where the input signal has much energy (Ljung 1999). Since a step-like input signal has its most energy at the frequencies which are used by the dead-time estimation, i.e at low frequencies, it will result in a good dead-time estimate.

In the open-loop simulations the fast second order system gave the lowest RMS error in average. This is natural because this system has a clearer start of the rise in the step response. The fourth order system with complex poles gave the highest RMS error. It seems that a more complex system makes the dead-time estimation more difficult.

For step signals in open and closed loop, we in this paper agree with (Horch 2000, Isaksson et al. 2001) that Laguerre DAP is a suitable estimation method. (Horch 2000, Isaksson et al. 2001) have, however, only tested with step signals and not discovered that the DAP method sometimes fails. Many of the factors and levels of the open loop simulations were chosen to be the same as in (Horch 2000) and (Isaksson et al. 2001), which makes comparisons easy.

5. CONCLUSIONS

We draw the following conclusions from the work presented in this paper:

• The DAP method can be used with any linear model structure.

• DAP methods are inherently non-robust and can totally fail in some cases.

• The failure probability of a DAP method depends on the used model structure, the input signal type and the SNR.

• In failing cases, DAP methods can be made more robust by zero guarding.

• An appropriate choice of zero guarding ap-pears to be robust and work for several model structures in both open loop and closed loop.

• For zero guarding of Laguerre DAP, removing zeros outside the unit circle close to +1 works better than removing zeros at other places outside but close to the unit circle. We must also allow more than one zero to be removed. • Most often the DAP methods work better

without prewhitening the data.

• Laguerre DAP without prewhitening seems to be a suitable dead-time estimation method for open loop with steps as input signal. Zero guarding should be done for safety.

• Laguerre DAP also seems to be a suitable dead-time estimation method for closed loop with steps as reference signal. The estima-tion should use the system input (not the reference signal) and no prewhitening. Zero guarding should be done for safety.

6. ACKNOWLEDGEMENTS

This work was supported by the Swedish Research Council.

7. REFERENCES

Bj¨orklund, Svante (2002). Analysis of a phase method for time-delay estimation. Techni-cal Report LiTH-ISY-R-2467. Department of Electrical Engineering, Link¨oping University. SE-581 83 Link¨oping, Sweden.

CSTB (n.d.). Matlab control system toolbox, v. 4.2.1 (R11.1). The Mathworks Inc.

Horch, A. (2000). Condition Monitoring of Con-trol Loops. Phd thesis TRITA-S3-REG-0002. Dep. of Signals, Sensors and Systems, Royal Institute of Technology. Stockholm, Sweden. Isaksson, A. J., A. Horch and G. A. Dumont

(2001). Event-triggered deadtime estimation from closed-loop data. In: Proc. of American Control Conference. Arlington, VA, USA. Isaksson, M. (1997). A comparison of some

ap-proaches to time-delay estimation.. Master’s thesis ISRN LUTFD2/TFRT–5580–SE. De-partment of Automatic Control, Lund Insti-tute of Technology. Lund, Sweden.

Ljung, Lennart (1999). System Identification: Theory for the User. 2nd ed.. Prentice-Hall. Upper Saddle River, N.J. USA.

Montgomery, D. C. (1997). Design and Analysis of Experiments. Wiley. ISBN 0-471-15746-5. SITB (n.d.). Matlab system identification toolbox

v. 5.0.1 (R12.1). The Mathworks Inc. Wahlberg, B. (1991). System identification using

Laguerre models. IEEE Transactions on Au-tomatic Control AC-36(5), 551–562.