- Email: [email protected]
Identification of time delays using a polynomial identification method
Systems & Control Letters 5 (1985) 267-271 North-Holland
Identification identification P.J. GAWTHROP
February 1985
of time delays using a polynomial method
and M.T. NIHTILA
School O/ .hgirreerirlg urd Applied Scierrces. The Sussex, Fuher. Brighrotr, Sussex BNI 9QT. U.K.
* Utlioersi!,~
oJ
Received 2 August 1984 Revised 23 October 1984 An algorithm for the exact least-squares identilication of an approximate continuous-time time-delay system is derived and its operation verified by simulation. Kq~worc/s: Time-delay, Continuous-time linear identification.
identification,
Non-
There are at least two problems associated with estimating T. Sensitivity based methods (e.g. that of Marshall [1,2]) require the derivatives of the measured signals; and the transfer function eTsT is a transcendental function of T. In this paper, the former problem is solved by the well-known statevariable filter [4] method, and the latter problem is solved by approximating ems’ by a polynomial function of T. The resultant approximate system is amenable to the polynomial identification methods of Nihtila [5].
2. System models 1. Introduction
The parametric identification of systems (continuous or discrete) containing an unknown time delay is difficult as such systems are not linear in the time-delay parameter. A number of approximate solutions to this problem have been proposed. Marshall [1,2] has discussed a continuoustime method based on computing the first sensitivity derivative. In discrete time, a number of methods have been derived [3] based on determining the delay by examining identified coefficients for approximately zero values. However, all these methods are based on ad-hoc approximations or heuristic manipulations with no clear theoretical foundation. Given the lack of soundly based methods, there is a need for a more precise analysis of this problem. In this paper, we consider the identification of the single parameter T in the equation
y(t)=u(t-
T)
(1.1)
from (noise-free) measurements of the input and output signals u(t) and y(t) respectively.
* On leave from the Helsinki University of Technology, Otaniemi, Otakaari 5a, 02150 Espoo 15, Finland 0167-6911/85/$3.30
The analysis given in this paper depends on two assumptions. The first assumption (Al) is that the initial state of the time delay is zero: (Al)
~(t)=0,
t
Taking Laplace transforms and using assumption (Al), Y(s) = e-W(s),
(2.1) where Y(s) and U(s) are the Laplace transforms of v(f) and u(l) respectively. The second assumption (A2) is that. the transfer function emsT may be adequately modeled by (A2)
Y(s) =
’ (1 +sT/N)~
w-
The right-hand side of this expression is a wellknown approximation to a time delay [l]. In particular, for any fixed value of s, the right-hand side of the expression (A2) approaches essT as N * co. In addition, for a fixed N, the right-hand side of the expression approaches eesT as ST+ 0. Specialising to the frequency domain, the frequency response of the approximate time delay (obtained by setting s = jw) can be made close to e-jw if oT is small enough and N is large enough. Thus the algorithm based on this approximation will fail if oT is large and N small.
0 1985, Elsevier Science Publishers B.V. (North-Holland)
267
Volume 5, Number
SYSTEMS& CONTROL LETTERS
4
It is, however, worth noting that, in practice, such a multiple lag model may be a more realistic representation of a system than a pure time delay. The contribution of this paper is to give an exact algorithm for computing the single parameter T for systems described by the approximate system of assumption (A2). But first, an alternative system model is developed. Introducing a stable polynomial C(s) with real coefficients and of degree N, N
c(s)= c CiSi, i-o
the approximate system model may be rewritten as
(1 +sT/N)’ Y(s) -LU(s). = c(s) C(s)
(2.2)
This leads to the following:
Proof. cyi is the ith coefficient of (1 + sT/N)~ (considered as a polynomial in ST). Equations (2.3) and (2.4) then follow directly from (2.2). 3. The identification
problem
The system representation of the lemma is polynomial in the single parameter T and involves no derivatives. As such, it is a suitable candidate for the polynomial identification method of Nihtila [5]. The identification problem can be stated as follows: Identification problem. For every final time 1 E [0, tr] and for a given input-output pair on [0, t], find the optimal value 19of the parameter 6’ which minimises the exponentially weighted performance index J(t, ~)=+2-V(tl-f30)2
Lemma. Giuen approximations (Al) and (A2), the system (2.2) may be rewritten in implicit form as g(t, T) = e(t)
(2.3)
where e(t) is the approximation g(t, T) = ;
February1985
error and
cu,T’y,(t) -uo(t),
+ f o’eep(r-T)g( -r, 8)* d7, J
where p (> 0) is the exponential weighting coefficient, 0, is an initial estimate and 7 is a positive initial weighting coefficient. This problem is solved by the following:
(2.4)
i-0
Theorem. If the optimal estimate 6 of T exists and
with the filtered derivatives yi( t ) and u. ( t ) given by (i=O,
. . . . N),
(2.5)
is unique for every T then it satisfies the following system of differential equations: J@+l*
*dt
l
=o
’
(3.2) .
$+f3J=hi+.lj+1g where L- ’ denotes inverse Laplace transform, with
and
hi=
Remark. The quantities yi can be generated by the following state-variable filter equations:
(3.3)
i i! j-o j!( j -j)!
gjgi-j-
(3.4
and J2N+l=
(2.8)
0.
... .
N-l).
(2.9)
(3.6)
The only non-zero initial conditions are
&o)=e,,
268
. . . ,2N),
(2.7)
i)!Ni’
(i=O,
(i=2,
where
N! ai= i!(N-
(3.1)
J2 (0) = 7.
(3.7)
Volume
5, Number
4
SYSTEMS
& CONTROL
The scalar quantities J., hi and gi are
(3.9)
,,,=--. 1 ai ’ 2 aei
(3.10)
Remark. The cost function J(t, e(t)) may be evaluated at each time instant from the Taylor series
J(t, S)= ‘c” $(t, e)(e- 8)‘. i-0
(3.11)
’
J, = J is given by (3.12), J1 = 0 and 4, i = 2, . . .‘, 2N are given by (3.3). This may be used to check for the existence of local minima in addition to that at 6’= 8. Proof. Taking
the partial derivative of the cost function with respect to time gives %+/3J=$g(t,
0)‘.
Repeated partial differentiation gives
(3.12) with respect to 0
(3.13)
s+/3J==hi.
The total derivative of 4. with respect to the final time t is dJ;. aA dt=at+aTTf
a4 de
1985
The algorithm requires the user to choose five parameters: the order N of the approximation in (2.1), the filter polynomial C(s) (2.2), the exponential weighting factor p (3.1), the initial value of the cost Jand the initial estimate 8, of the time delay (3.1). As this algorithm has only been recently developed only a rudimentary guide to the choice of these parameters can be given. Further work is required to produce more complete guidelines. A large value of N gives a better approximation in (2.1), but gives more computation and can give numerical difficulties in implementing the statevariable filter equations (2.8) and (2.9). Its chaise depends on the value of o,T where w. is the highest significant frequency in the input signal u(t)* The choice of the coefficients of C(s) is not crucial in the noise-free case considered. Longer time constants give smoother filtered signals yi. One possibility, which has proved satisfactory in simulations, is to make the N roots of C real with time constants approximately equal to an a-priori guess of the time delay. The exponential weighting fi allows for discounting of old data. fl= 0 weights all data equally while a large positive value of p uses only recent data in an exponential window with time constant p-1. The initial values t9, and J allow an initial guess 0, of the delay to be incorporated. The magnitude of 7 corresponds to the weighting of this initial guess with respect to measured data. 4. Simulation of the algorithm
(3.14)
noting that aJ,/at = Ji+i gives (3.3). The optimum solution (if it exists and is unique) is given by setting J1(t, @=o.
February
Choice of algorithm parameters
(3.8)
aig gi=,ei,
LETTERS
(3.15)
This gives (3.2). Equation (3.4) follows from Leibniz’s rule and (3.5) from repeated differentiation of (2.4). Equation (3.6) follows as the cost function J contains no power of 8 greater than 2N. Differentiating the cost J at t = 0 and using (3.15) gives (3.7).
To verify the operation of the algorithm, and to illustrate its properties, it was programmed in Pascal and, together with the time-delay system of equation (l.l), was implemented on a DEC LSI 11/23 processor using double-precision floatingpoint arithmetic. Digital simulation of continuous-time systems clearly involves approximation. No attempt was made to optimise the approximation in any way; simple Euler integration was used. For the purposes of the simulation the sample interval was reduced until no appreciable change in the results occurred. 269
SYSTEMS L CONTROL LETTERS
Volume 5. Number 4
February
19X5
Table 1 Parameter
Value
N C(s) P i 43
6 (l+s)N 1.0 0.01 0.0
The simulated delay was 0.5 units. The estimator parameters were chosen as in Table 1. A sinusoidal input with period 10 units was applied to the delay element which had zero initial conditions. Figure 1 shows the input (u) and output (y) signals plotted against time. Figure 2 shows the estimated delay (4) plotted against time. Figure 3 shows the cost function J(t, 0) plotted against 8 in the range [0, l] for 1= 0, 1, 2 and 3. The optimal estimate approaches the true value in about 1.25. time units. The fact that the error does not approach zero is due to the approximation inherent in assumption (A2) together with the large exponential weighting factor p. Simulations with other parameter values gave similar results. Sinusoidal inputs with periods less than about 2 units did, however, give poor results. This is due to the approximation (A2) being poor
yJ!:
3.0
Fig. 1. The input and output signals.
1.0 t
ii olrTKy,o
4.0
5.0 Time
Fig. 3. The cost function J.
at high frequencies. This result implies that input signals with significant components at frequencies above about 1/2T would also give poor results. 5. Concluding
remarks
An exacr least-squares algorithm for an aptime-delay system has been derived, and its operation verified by simulation. Because the approximate delay model is polynomial in a single parameter (T), the algorithm is simple, compared with the general result of Nihtila [5]. The two approximations made in this paper are that the delay initial conditions are zero, and that a time delay can be adequately modelled by a multiple lag system. The effects of each of these approximations requires further analysis. The paper only considers a noise-free model. It does, however, have a differentiator-free implementation and so should have good properties in the presence of noise. Further work is required to investigate bias and variance of parameter estimates in,a stochastic environment and to examine alternative estimation schemes in these circumstances. A more complex system with more parameters (for example a time delay in series with a rational transfer function) can be treated in a similar way [5] but the resultant algorithms become complicated. The effect of non-zero system initial conditions may also be accounted for by introducing further paramters [6].
proximate
References 1.0
.
Fig. 2. The estimated time delay.
270
[l] J.E. Mt~shall, Conrrolo/ ‘fitne- De/q Sysrew (Peter Peregrinus, Stevenage, 1979).
Volume 5, Number 4
SYSTEMS Sr CONTROL LETTERS
[2] J.E. Marshall, Identification strategies for time-delay systerns, in: H. Unbehauen, Ed., Methods and Applications in Adapfiue Confrol (Springer, Berlin-New York, 1980) 141-150. [3] H. Km-z and W. Goedecke, Digital parameter-adaptive control of processes with unknown time-delay, Auromarica 17 (1981) 245-252.
February 1985
[4] P.C. Young, Parameter estimation for continuous-time models - A survey, Automuticu 17 (1981) 23-39. [5] M. Nihtila, Optimal finite-dimensional recursive identification in a polynomial output mapping class, Sysfefns Cortfrol Let!. 3 (1983) 341-348. (61 P.J. Gawthrop, Parametric identification of transient signals, IMA .I. Muth. Conrrol und Itrfortn. 1 (1984) 117-128.
2il
Identification identification P.J. GAWTHROP
February 1985
of time delays using a polynomial method
and M.T. NIHTILA
School O/ .hgirreerirlg urd Applied Scierrces. The Sussex, Fuher. Brighrotr, Sussex BNI 9QT. U.K.
* Utlioersi!,~
oJ
Received 2 August 1984 Revised 23 October 1984 An algorithm for the exact least-squares identilication of an approximate continuous-time time-delay system is derived and its operation verified by simulation. Kq~worc/s: Time-delay, Continuous-time linear identification.
identification,
Non-
There are at least two problems associated with estimating T. Sensitivity based methods (e.g. that of Marshall [1,2]) require the derivatives of the measured signals; and the transfer function eTsT is a transcendental function of T. In this paper, the former problem is solved by the well-known statevariable filter [4] method, and the latter problem is solved by approximating ems’ by a polynomial function of T. The resultant approximate system is amenable to the polynomial identification methods of Nihtila [5].
2. System models 1. Introduction
The parametric identification of systems (continuous or discrete) containing an unknown time delay is difficult as such systems are not linear in the time-delay parameter. A number of approximate solutions to this problem have been proposed. Marshall [1,2] has discussed a continuoustime method based on computing the first sensitivity derivative. In discrete time, a number of methods have been derived [3] based on determining the delay by examining identified coefficients for approximately zero values. However, all these methods are based on ad-hoc approximations or heuristic manipulations with no clear theoretical foundation. Given the lack of soundly based methods, there is a need for a more precise analysis of this problem. In this paper, we consider the identification of the single parameter T in the equation
y(t)=u(t-
T)
(1.1)
from (noise-free) measurements of the input and output signals u(t) and y(t) respectively.
* On leave from the Helsinki University of Technology, Otaniemi, Otakaari 5a, 02150 Espoo 15, Finland 0167-6911/85/$3.30
The analysis given in this paper depends on two assumptions. The first assumption (Al) is that the initial state of the time delay is zero: (Al)
~(t)=0,
t
Taking Laplace transforms and using assumption (Al), Y(s) = e-W(s),
(2.1) where Y(s) and U(s) are the Laplace transforms of v(f) and u(l) respectively. The second assumption (A2) is that. the transfer function emsT may be adequately modeled by (A2)
Y(s) =
’ (1 +sT/N)~
w-
The right-hand side of this expression is a wellknown approximation to a time delay [l]. In particular, for any fixed value of s, the right-hand side of the expression (A2) approaches essT as N * co. In addition, for a fixed N, the right-hand side of the expression approaches eesT as ST+ 0. Specialising to the frequency domain, the frequency response of the approximate time delay (obtained by setting s = jw) can be made close to e-jw if oT is small enough and N is large enough. Thus the algorithm based on this approximation will fail if oT is large and N small.
0 1985, Elsevier Science Publishers B.V. (North-Holland)
267
Volume 5, Number
SYSTEMS& CONTROL LETTERS
4
It is, however, worth noting that, in practice, such a multiple lag model may be a more realistic representation of a system than a pure time delay. The contribution of this paper is to give an exact algorithm for computing the single parameter T for systems described by the approximate system of assumption (A2). But first, an alternative system model is developed. Introducing a stable polynomial C(s) with real coefficients and of degree N, N
c(s)= c CiSi, i-o
the approximate system model may be rewritten as
(1 +sT/N)’ Y(s) -LU(s). = c(s) C(s)
(2.2)
This leads to the following:
Proof. cyi is the ith coefficient of (1 + sT/N)~ (considered as a polynomial in ST). Equations (2.3) and (2.4) then follow directly from (2.2). 3. The identification
problem
The system representation of the lemma is polynomial in the single parameter T and involves no derivatives. As such, it is a suitable candidate for the polynomial identification method of Nihtila [5]. The identification problem can be stated as follows: Identification problem. For every final time 1 E [0, tr] and for a given input-output pair on [0, t], find the optimal value 19of the parameter 6’ which minimises the exponentially weighted performance index J(t, ~)=+2-V(tl-f30)2
Lemma. Giuen approximations (Al) and (A2), the system (2.2) may be rewritten in implicit form as g(t, T) = e(t)
(2.3)
where e(t) is the approximation g(t, T) = ;
February1985
error and
cu,T’y,(t) -uo(t),
+ f o’eep(r-T)g( -r, 8)* d7, J
where p (> 0) is the exponential weighting coefficient, 0, is an initial estimate and 7 is a positive initial weighting coefficient. This problem is solved by the following:
(2.4)
i-0
Theorem. If the optimal estimate 6 of T exists and
with the filtered derivatives yi( t ) and u. ( t ) given by (i=O,
. . . . N),
(2.5)
is unique for every T then it satisfies the following system of differential equations: J@+l*
*dt
l
=o
’
(3.2) .
$+f3J=hi+.lj+1g where L- ’ denotes inverse Laplace transform, with
and
hi=
Remark. The quantities yi can be generated by the following state-variable filter equations:
(3.3)
i i! j-o j!( j -j)!
gjgi-j-
(3.4
and J2N+l=
(2.8)
0.
... .
N-l).
(2.9)
(3.6)
The only non-zero initial conditions are
&o)=e,,
268
. . . ,2N),
(2.7)
i)!Ni’
(i=O,
(i=2,
where
N! ai= i!(N-
(3.1)
J2 (0) = 7.
(3.7)
Volume
5, Number
4
SYSTEMS
& CONTROL
The scalar quantities J., hi and gi are
(3.9)
,,,=--. 1 ai ’ 2 aei
(3.10)
Remark. The cost function J(t, e(t)) may be evaluated at each time instant from the Taylor series
J(t, S)= ‘c” $(t, e)(e- 8)‘. i-0
(3.11)
’
J, = J is given by (3.12), J1 = 0 and 4, i = 2, . . .‘, 2N are given by (3.3). This may be used to check for the existence of local minima in addition to that at 6’= 8. Proof. Taking
the partial derivative of the cost function with respect to time gives %+/3J=$g(t,
0)‘.
Repeated partial differentiation gives
(3.12) with respect to 0
(3.13)
s+/3J==hi.
The total derivative of 4. with respect to the final time t is dJ;. aA dt=at+aTTf
a4 de
1985
The algorithm requires the user to choose five parameters: the order N of the approximation in (2.1), the filter polynomial C(s) (2.2), the exponential weighting factor p (3.1), the initial value of the cost Jand the initial estimate 8, of the time delay (3.1). As this algorithm has only been recently developed only a rudimentary guide to the choice of these parameters can be given. Further work is required to produce more complete guidelines. A large value of N gives a better approximation in (2.1), but gives more computation and can give numerical difficulties in implementing the statevariable filter equations (2.8) and (2.9). Its chaise depends on the value of o,T where w. is the highest significant frequency in the input signal u(t)* The choice of the coefficients of C(s) is not crucial in the noise-free case considered. Longer time constants give smoother filtered signals yi. One possibility, which has proved satisfactory in simulations, is to make the N roots of C real with time constants approximately equal to an a-priori guess of the time delay. The exponential weighting fi allows for discounting of old data. fl= 0 weights all data equally while a large positive value of p uses only recent data in an exponential window with time constant p-1. The initial values t9, and J allow an initial guess 0, of the delay to be incorporated. The magnitude of 7 corresponds to the weighting of this initial guess with respect to measured data. 4. Simulation of the algorithm
(3.14)
noting that aJ,/at = Ji+i gives (3.3). The optimum solution (if it exists and is unique) is given by setting J1(t, @=o.
February
Choice of algorithm parameters
(3.8)
aig gi=,ei,
LETTERS
(3.15)
This gives (3.2). Equation (3.4) follows from Leibniz’s rule and (3.5) from repeated differentiation of (2.4). Equation (3.6) follows as the cost function J contains no power of 8 greater than 2N. Differentiating the cost J at t = 0 and using (3.15) gives (3.7).
To verify the operation of the algorithm, and to illustrate its properties, it was programmed in Pascal and, together with the time-delay system of equation (l.l), was implemented on a DEC LSI 11/23 processor using double-precision floatingpoint arithmetic. Digital simulation of continuous-time systems clearly involves approximation. No attempt was made to optimise the approximation in any way; simple Euler integration was used. For the purposes of the simulation the sample interval was reduced until no appreciable change in the results occurred. 269
SYSTEMS L CONTROL LETTERS
Volume 5. Number 4
February
19X5
Table 1 Parameter
Value
N C(s) P i 43
6 (l+s)N 1.0 0.01 0.0
The simulated delay was 0.5 units. The estimator parameters were chosen as in Table 1. A sinusoidal input with period 10 units was applied to the delay element which had zero initial conditions. Figure 1 shows the input (u) and output (y) signals plotted against time. Figure 2 shows the estimated delay (4) plotted against time. Figure 3 shows the cost function J(t, 0) plotted against 8 in the range [0, l] for 1= 0, 1, 2 and 3. The optimal estimate approaches the true value in about 1.25. time units. The fact that the error does not approach zero is due to the approximation inherent in assumption (A2) together with the large exponential weighting factor p. Simulations with other parameter values gave similar results. Sinusoidal inputs with periods less than about 2 units did, however, give poor results. This is due to the approximation (A2) being poor
yJ!:
3.0
Fig. 1. The input and output signals.
1.0 t
ii olrTKy,o
4.0
5.0 Time
Fig. 3. The cost function J.
at high frequencies. This result implies that input signals with significant components at frequencies above about 1/2T would also give poor results. 5. Concluding
remarks
An exacr least-squares algorithm for an aptime-delay system has been derived, and its operation verified by simulation. Because the approximate delay model is polynomial in a single parameter (T), the algorithm is simple, compared with the general result of Nihtila [5]. The two approximations made in this paper are that the delay initial conditions are zero, and that a time delay can be adequately modelled by a multiple lag system. The effects of each of these approximations requires further analysis. The paper only considers a noise-free model. It does, however, have a differentiator-free implementation and so should have good properties in the presence of noise. Further work is required to investigate bias and variance of parameter estimates in,a stochastic environment and to examine alternative estimation schemes in these circumstances. A more complex system with more parameters (for example a time delay in series with a rational transfer function) can be treated in a similar way [5] but the resultant algorithms become complicated. The effect of non-zero system initial conditions may also be accounted for by introducing further paramters [6].
proximate
References 1.0
.
Fig. 2. The estimated time delay.
270
[l] J.E. Mt~shall, Conrrolo/ ‘fitne- De/q Sysrew (Peter Peregrinus, Stevenage, 1979).
Volume 5, Number 4
SYSTEMS Sr CONTROL LETTERS
[2] J.E. Marshall, Identification strategies for time-delay systerns, in: H. Unbehauen, Ed., Methods and Applications in Adapfiue Confrol (Springer, Berlin-New York, 1980) 141-150. [3] H. Km-z and W. Goedecke, Digital parameter-adaptive control of processes with unknown time-delay, Auromarica 17 (1981) 245-252.
February 1985
[4] P.C. Young, Parameter estimation for continuous-time models - A survey, Automuticu 17 (1981) 23-39. [5] M. Nihtila, Optimal finite-dimensional recursive identification in a polynomial output mapping class, Sysfefns Cortfrol Let!. 3 (1983) 341-348. (61 P.J. Gawthrop, Parametric identification of transient signals, IMA .I. Muth. Conrrol und Itrfortn. 1 (1984) 117-128.
2il