- Email: [email protected]
On-line estimation of the time delay via orthogonal collocation
SIMULATION PRACTICE Simulation
Practice
and Theory
= THEORY
5 (1997) lOlL120
On-line estimation of the time delay via orthogonal collocation M. NihtilZ a,*, T. Damak b. J.P. Babary ’ a University of Kuopio, Faculty of Natural and Environmental Sciences, Department of Computer Science and Applied Mathematics, P.O. Box 1627, FIN-7021 1 Kuopio, Finland b Ecole Nationale d’lngtnieurs de Sfax. Department de Genie Micanique, B.P.W. 3038 Sfax, Tunisie ’ Centre National de la Recherche Scienti$que, Laboratoire d’Analyse et d’Architecture des Syst?mes, 7 avenue du Colonel Roche, F-31077 Toulouse Cedex, France Received
12 February
1995, revised 11 November
1995
Abstract The objective of this work is to develop a new real-time time-delay estimator for the input delay of the otherwise finite-dimensional single-input single-output linear differential system. The method applies also to tracking of possible time-variations in the delay. The key idea in the construction of the estimator is to transform the delay part of the system into a transport system described by a linear partial differential equation (PDE). The PDE is linearly parametrized by the delay argument. The PDE is approximated by a finite-dimensional ordinary differential equation (ODE) system. The approximation is obtained by applying orthogonal collocation in the spatial descretization of the PDE. A fixed-interval exponentially weighted cost functional is then applied in the ODE system in formulating the delay estimation problem in mathematical terms. The data forgetting issue allows tracking of changes in the delay. Minimization of the cost functional is carried out approximately to obtain a finite-dimensional recursive delay estimator for increasing values of the final time of the integral in the cost. Simulations for a fixed as well as for a time-varying delay are given to demonstrate estimation and tracking properties of the identifier. Key words: simulation
Delay
estimation;
Least-squares
identification;
Collocation
method;
Model
1. Introduction The destabilizing effect of time delays is a known fact in closed-loop control the delays are ignored in the control law calculations, see e.g. [lo]. Different * Corresponding
author.
Email: [email protected].
0928-4869/97/$17.00 Copyright simpraOO82XX PZZ 0928-4869(95)00052-6
0 1997 - Elsevier Science B.V. All rights reserved
when types
102
M. Nihtilii et al. J Simulation Practice and Theory 5 (1997) 101-120
of predictive control algorithms have been developed to overcome these problems by taking into account the system delays in the controller design. Stable systems can be controlled by applying Smith’s predictor [5,21]. Unstable systems need so-called analytic predictors or finite spectrum assignment controllers, which have been studied among others by Manitius and Olbrot [13], Furukawa and Shimemura [6], and Ichikawa [ 111. Predictive methods result in closed loop responses whose performance is comparable with those of the controlled systems without delays. Most predictive methods directly utilize numerical values of the delays. Consequently, it is necessary to know, estimate or track on-line values of the delay to achieve sufficiently good closed-loop performance. The transcendental nature of the transfer functions of linear input delay systems causes problems in the application of ordinary control techniques. However, there are several methods to approximate the transcendental transfer function of the delay by rational ones as it has been described e.g. by Agarwal and Canudas [ 11, Gawthorp and Nihtila [7], and Paynter [ 181. Then these approximations are parametrized by some functions of the delay argument. In using rational approximations of the delay it is essential from the identification viewpoint that the approximations do not include too many free parameters. PadC approximations, tank series models and allpass filters are among the best known rational approximations. These can be obtained by applying frequency-response techniques. For discrete-time models the inclusion of a (variable) delay increases generally the number of free parameters. In some special cases the delay approximation can be based on single parameter models also in a discrete formulation as presented by Bjorklund et al. [ 2,3]. The goal of this paper is to construct an estimator for the input delay of a linear system based on given measurements of the input and the output. Furthermore, the identification problem should be recursively solvable to be applicable in real-time purposes. Due to the nonlinear character of the effect of the input delay, it can be expected that only finite-dimensional recursive estimators, which are approximate with respect to some feasible criteria, can be constructed. The identification of time-delays is not an easy task. Some different approaches can, however, be found. Among the first studies describing on-line identification of a variable delay is the one presented in [20]. Pearson and Wuu [19] applied a decomposition technique in the off-line environment to identify a constant delay separately from the other system parameters. Several recursive variants, applying delay approximation or direct gradient-based approximate recursion, are studied by Gawthorp et al. [S]. Modern state-space identification packages are also applicable in delay identification as described recently in [ 151. In this paper an alternative approach for recursive input delay identification is presented. It is based on a specific rational approximation of the delay. The input delay is considered as a transport delay. Then it is exactly described by a linear firstorder partial differential equation (PDE) including the delay as a single parameter. This PDE model is approximated by an ordinary differential equation (ODE) system. The approximation procedure is based on the method of orthogonal collocation, see e.g. [23]. Specific applications of this approximation method have been studied
M. Nihtilii et al. 1 Simulation
103
Practice and Theory 5 (1997) 101-120
among others by Molander [ 141, Dochain et al. [4], and Nihtila et al. [ 171. In the method proposed here the resulting ODE system is parametrized in a polynomial way by the delay. An approximate recursive least-squares method in continuoustime form is then applied in the delay estimation. The efficiency of any numerical method, in our case recursive identification of the time-delay, can be evaluated either via the theoretical analysis of the algorithms at hand or via simulations for some specific parameter values. Due to the nonlinear character of the problem we have chosen the simulation way. The simulations actually show the feasibility of the technique developed. The organization of this paper is as follows. The second section treats the system models and introduces the procedure to approximate the delay system by a rational transfer function parametrized by the single delay value. In the third section the identification problem is stated in mathematical terms as an optimization problem. Next the approximate on-line identifier is constructed. In the fourth section a secondorder underdamped system with input delay is studied via simulations to test the algorithm and its feasibility. Constant as well as time-varying delay cases are examined. Some concluding remarks complete the paper. In the Appendix the orthogonal collocation method is revisited with some details.
2. System models The following is considered: 4P)Y(t)
differential
= &+4t
A and B are polynomials given by
-
system with a constant
(or time-varying)
input
(1)
h).
of known
A(p)=a,+a,p+a,p2+
delay h
degrees
in the derivative
operator
p. They are
... +a$“-‘+p”,
B(p) = b, + b2p + b3p2 + ... + bmpm-l,
(2) m d n,
(3)
in which IZ and m - 1 are the degrees of A and B, respectively. In transfer function formalism the system is depicted in Fig. 1. A state variable model corresponding to
u
>
,-b
X
>
B(P) A(P)
Fig. 1. Original
input delay system
Y
)
M. Nihtiliiet al. / SimulationPracticeand Theory 5 (1997) 101-120
104
the system (1) is in canonical observable form given by i(t) = A,x(t) + B,u(t - h),
(4)
y(t) = Cx(t),
(5)
in which the system matrices of appropriate 0
A, =
0
...
0
1 0
...
0
0
f.. ..
0 :
-a2 -a3
1
-a,
1
OO...
c=
dimensions are given by
-a,
’
B, =
[O ... 0 11.
(6 7)
(8)
The first step in system identification considered is to estimate the delay in the case of a known rational part, i.e. G,(p) = B( p)/A(p), of the transfer function. In order to be able to apply input-output techniques without inverting G,(p) the order of the transfer functions is reversed in Fig. 1. If the two parts of the system are timeinvariant then the order can be reversed without affecting the input-output behaviour of the system. This is given in Fig. 2. 2.1. Delay approximation In the framework of flow processes, studied thoroughly e.g. by Zenger [24], the delay part of the system in Fig. 2 is modelled as a distributed parameter system described by the following partial differential equation with a relevant boundary condition and the given output equation (cf. [9]) Wt(Z,t) = -( l/h)w,(z, t),
0 < z < 1,
(9)
NO, t) = v(t),
(10)
Y(L)= w( 1, 0.
(11)
The subindices t and z denote partial derivatives with respect to the time and space
bp+j-iqL> Fig. 2. Reversed order model.
M. Nihtilii et al. J Simulation Practice and Theory 5 (1997) 101-120
105
variables, respectively. The input v of the delay system is obtained from the ordinary linear system 4pMt)
= NPMO.
(12)
It must be emphasized that the system (9)-( 11) is an equivalent representation of the delay part of Fig. 2. Furthermore, w(z, t) corresponds to v(t - zh), and ~(0, t) (and w,(t) in (13)) to v(t). The orthogonal collocation method, cf. [23], which will be applied in the PDE system (9)-( 11) is a spatial discretization method. When applied, it converts the PDE system into an ordinary differential equation system which is an approximate input-output description of the delay part of Fig. 2. The details are given in the Appendix. By applying the series N+l w(z,
t) =
1
BdzMt),
(13)
i=O
with appropriately chosen base functions pi in the PDE system the following finitedimensional approximate state variable system is obtained for the vector 4, where
k(t) = F,q(t) + 4v(t),
(15)
Yd(C h) =
(16)
Hdt).
A new notation yd(t; h) is introduced to emphasize that it is an approximation depending on the delay. The following expressions for the elements of the matrices are obtained (F,)ij=
-p(i(Zi), i,j= 1, . . . . N + 1,
(O1)i = -p;(zi), H=[O
i = 1, . . . . N + 1,
... 0 11,
(17) (18) (19)
in which the prime ’stands for derivation with respect to the space variable z. The points zo, zl, . . . . zN+i are appropriately chosen successive points on the spatial interval [0, 11. N is the number of internal points in the discretization. Based on the specific position of the delay in (15), i.e. the derivative operator p is everywhere multiplied by the delay h, the following unique polynomial operator representation WP)YdC
h) = WPMt)
(20)
is obtained for the delay approximation from the state-variable representation (15)-( 16). The coefficients of the polynomials F(s) and D(S) corresponding to the matrices F, and D, depend only on the tuning parameters of the collocation method a and j3 and on the number of the collocation points N (cf. the Appendix). The
M. Nihtilli et al. / Simulation Practice and Theory 5 (1997) 101-120
106
polynomials do not depend on the specific value of the delay. Consequently, the system (20) is a rational single parameter approximation of the delay system. A practical point is also that the single parameter is the delay itself.
3. Least-squares identification Minimization of the error between the measured and model generated output, i.e. of y(t) - yd(t; h), with respect to the delay h is formulated as an exponentially weighted least-squares identification problem. Then, without stochastic concepts, measure-theoretic problems arising in stochastic time-continuous formalism can be avoided. In the following the formulation of the problem and an approximate recursive solution will be given. 3.1. Problem statement The identification of the delay is stated in mathematical terms as a set of fixedinterval optimization problems. Consequently, for each fixed final time t 2 0 solve the following task. Fixed interval optimization problem. For the given input v and the output y on the interval [0, t] find the value of the delay, denoted by h(t), which minimizes the cost
functional t V(t,h)=~e-“‘J(h-h,)‘+~
ea(s-r)[y(s) - y,,(s; h)]’ ds, s0
(21)
subject to the constraint Us) = F14(s) + 444, Y&
4 = Rl(d.
4(O) = 40,
(22) (23)
for 0 < s < t. The given scalars a, ho, and J denote the forgetting exponent, the initial delay approximation and the (positive) weight of the initial cost, respectively. If the problem statement were given by using filtered signals, then a parameterpolynomial identification scheme could have been used [ 161. This is due to the fact that the operator polynomials F(hp) and D(hp) are polynomials also in the parameter h. Then the filtered model will be in a polynomial form with respect to the delay amenable for the parameter-polynomial identification. The polynomial identifier is optimal in the sense that a finite-dimensional recursive algorithm already gives the optimal, not approximate, least-squares fit in the approximate filtered model. Due to relative complexity of the parameter-polynomial method, however, here an approximate algorithm is used. This approximate algorithm is not optimal in the sense mentioned above.
M. Nihtilii et al. / Simulation
Practice and Theory 5 (1997) 101-120
107
3.2. Approximate recursion for the delay The gradient-based recursion for the delay estimate i(t) is given by the two basic scalar differential equations (see derivation in Section 3.4) h(t) = M(t)_‘$(t; &l(t) = -aM(t)
&))CYW - y,(t; + $(t; L(t))‘,
fiW)ls
i(O) = h,,
M(0) = .I.
(24) (25)
The function rl/ denotes the partial derivative or gradient of the model output, i.e.
$(t; h) =
4 ydt; h).
3.3. Gradient system Calculation of the gradient is based on the following procedure. The approximate system represented in the transfer function form (27)
Wp)Yd(t; h) = N&(t) is differentiated with respect to the delay h giving the system
(28)
F,,(hp)YJt; h) + Wp)ll/(t; h) = Q,(hp)o(t). The subindex h denotes differentiation with respect to the delay. It is seen that Uhp) = M”(hp),
(29)
&(hp) = @VP),
(30)
in which the prime ’indicates differentiation
of the polynomials F and G, i.e.
df’b) F’(s)= y, D’(s) =
dW
(32)
ds.
If the number of the internal collocation (i.e. spatial discretization) points, when approximating the system (9)-( 11) by the system (15)-( 16), is N, then the gradient system (28) can be written in the form +N+
1 yd(t; h), I
(33)
in which the following definitions are applied,
QW=
W'Vvh
R(hp) = hpF’(hp) -(N
(34) + l)F(hp).
(35)
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
108
The actual realization of the gradient is obtained by using two state-variable filters
s,: 4 = FIP + !A% p(O)= 0,
(36)
t=HP, S,:
w(O)=O,
hti=F,o+Rly,,
(37)
i=Ho+(N+l)y,, with the joint output equation
(38) in which p and o are (N + 1)-vectors. The filters S1 and SZ correspond to the transfer functions of the coefficients of v(t) and yd(t; h) of (33) (without the scalar coefficient l/h, which is included in the output equation (38)). The vectors Q, and R, in (36) and (37) correspond to realizations of the Q(s) and R(s) polynomials (cf. the corresponding realization of ( 12) by the matrices (6)-( 8)) in the respective transfer functions Q(s)/F(s) and R(s)/P(s). Consequently, they are independent of the delay. It has to be noted that, in estimating the delay via the system (24)-(25), the current value of the estimate h(t) must be used in the state variable filters and in the joint output equation. 3.4. Derivation of the identijier The minimizing condition for the cost (21) is that its gradient is equal to zero for each final time t, i.e. I/h(t, h) = e-“‘J(h
- h,) -
f J
e”(S-f)$(s; h)[y(s) - yd(s; h)] ds = 0,
(39)
0
for the optimal estimate h(t). Consequently, we have for all t E [0, T] the identity qt,
(40)
h(t)) = 0.
Total differentiation of the gradient condition with respect to the final time t (considered now as a variable) gives the condition -a&At, i(t)) - $(t;
&MY(t)
-Yy,(C
&t))l
f
- Jf
+ e-“‘J+ e”‘“-‘)J/(s;L(t))’ds i s0 eU(S-‘)tj,(s; h(t))[y(s)
0
- y,(s; h(t))] ds k(t)) = 0,
(41)
Al. Nihtilii et al. J Simulation Practice and Theory 5 (1997) 101-120
109
Introduce then the following scalar functions f M(t)=
eCatJ
eacspt)$(s; A(s))~ds,
+
(42)
s0 f N(t)=
e-a'J
ea(s-t)$(s; h(t))' ds
+
f -1
c0
ea(s-t)$h(s;
h(t))[y(s)-
yd(s; 6(t))] ds.
(43)
0
Then the truly optimal estimator is given by (43) and (24) in which the gain M(t)- ’ is substituted by N(t)-’ assumed that it exists. Because the optimal inverse gain (43) is not recursively computable by a finite-dimensional ODE system but only in some special cases it is approximated by the suboptimal inverse gain (42) which satisfies the differential equation and the initial condition in (25). The approximation of N(t) by M(t) needs two simplifications. The last term in (43) is actually zero if the residual, or the output modelling error y(s) - y,(s; g(t)) is uncorrelated with $(s, G(t)) for 0 6 s < t. Furthermore, if the estimated delay value changes sufficiently slowly on the effective interval of the exponential weighting function exp[a(s - t)] then the estimated delay k(t) in (43) can be substituted by A(s).This second assumption finally results in the inverse gain equation (42). The initial conditions of the estimator are due to the initial cost in (21). In the derivation it is assumed that h is a differentiable function of its argument t. This includes that the cost has such a unique global minimum, with respect to the delay, which does not jump between possible local minima, when the time t increases from 0 to T. This issue is considered in connection with a parameter-polynomial cost functional in [ 161.
4. Simulation study 4.1. Model approximation For the delay approximation feasible values of the tuning parameters a, p and N have to be chosen. The selection of CIand p is based on a trial-and-error technique. A feasible value of N is obtained related to the accuracy of the approximation required on the spatial domain of the PDE system (9)-( 1l), which is equivalent with the delay part of the original model (cf. [22]). The specific process model studied here is given by G(p) = Go(P)eehp,
(44)
4 G,(P)
=
p2 + 2j0,p
+ w;
A step response of the approximation,
(45)
i.e. of Go and the system (15)-( 16) is visually
110
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
fitted with the step response of G by varying the tuning parameters. This was carried out in the noiseless case. The values of the system and tuning parameters are given in Table 1. In all simulations outputs of the system model, the estimator and the state-variable filters were generated at discrete intervals of length Td = 0.1 s. A fourth-order RungeKutta method with the integration step size z = 0.01 s was applied to generate the signals. The responses of the tuned approximation and the second-order model are depicted in Fig. 3 and are almost identical from the visual viewpoint. The difference y - y,(.; h) is represented in Fig. 4 but in a different amplitude scale. The coefficients of the approximation polynomials F(s) =fr +fzs + .** +f,P+
s7,
(46)
D(s) = d, + d,s + . . . + d,@, are given in Table 2. The zeros (and the discretization
(47) points) of the interpolation
polynomial L(z) (cf.
Table 1 Values of the system and tuning parameters
wo i h
0.2 rad/s 0.5 radfs 20 s
N
Fig. 3. Step responses of the input delay (-)
6 1 1
and the approximate model (- - - -).
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
111
0.015 -
0.01 -
0.005 -
o-
-0.005 -
-0.01 -
-0.015 0
20
40
60
80
time (s) Fig. 4. Difference
Table 2 Coefficients
of the polynomials
of the step responses
of Fig. 3.
F(s) and D(s)
2.1622 x lo6 1.2355 x lo6 5.544 x xlo4 3.3264 lo5 6.3 x lo3 504 28
2.1622 x lo6 - 9.2664 x 10’ - 1.782 1.98 x xlo4 lo5
d, d, d, d, d, d,
1.35 x 103 -54 1
the Appendix) are depicted in Fig. 5. Due to the tuning values CI= /? = 1 the points are symmetrically located around z = 0.5. The matrix F, and vector D, of the statevariable system are given in Table 3 in the form [F,ID,]. Visual comparison of frequency responses of the delay model with various approx-
z()=
0
z,=
Fig. 5. Location
of collocation
points
for a = p = 1.
1
112
M. Nihtilii et al. j Simulation Practice and Theory 5 (1997) 101-120
Table 3 Elements of the enlarged matrix [F, -1.0 x 1O-6 5.6130 -2.1579 1.3223 - 1.0741 1.1471 -2.5954
-9.0872 0 4.7564 - 2.2707 1.6900 - 1.7389 3.8833
1Ill]
4.2241 -5.7510 1.0 x 1O-6 4.7779 -2.7456 2.5885 - 5.6204
- 2.5885 2.7456 -4.7778 0 5.7510 -4.2241 8.5958
1.7389 - 1.6900 2.2707 - 4.7564 2.0 x 1O-6 9.0872 -15.139
- 1.1471 1.0741 - 1.3223 2.1579 -5.6130 -7.0 x 1o-6 37.875
0.4399 -0.4066 0.4867 -0.7443 1.5850 - 6.4198 -28.0
6.4198 - 1.5850 0.7443 -0.4867 0.4066 -0.4399 1.0
imations gives an insight of the accuracy of the approximations. Here, in addition to the collocated approximation, a tank series approximation of the delay system is taken for comparisons. Frequency responses of the second-order delay model (44)-(45), the collocated model (12) with the signal u generated by Go and the tank series model of the delay with Go are given in Fig. 6. The number of the tanks has been chosen equal to the number of the state variables in the collocated model. This gives the fair basis for comparing the performance of the models on the frequency band considered. The tuning values chosen for the collocation method emphasize accurate approximation of the system gain and phase on the frequency band 0
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
1.2
IGO’V 41 -r
T
1711
-
T
T
113
-
R
..,.... /
.
:
1
0.8
I. ,.j .,...., ..,
c
,
”
/
..
1.
;
:. .:: : :
;;
0.6
..i
i. :.:
.i
:
;
:
.;
‘: i
0.4
“”
:
;
:
1: . :.
.:
j.i.:_
..:. ;
.j.
:
: : : . .
‘:’
:
‘:
^
.
.I
.
^.
.
A.,
.:
.:.
,..
0.2
I
I
‘i. :
:.:.:
.:.:.:
:
:
: :
:
.:
:
;.:.
.__
:
i .:.
d PO-2
i
i
:
:: :
+_ -;
,”
i
ii
10”
AL
1
argV3.W) Meg1
-200 -
-KJO-
-loo0 ‘S.-b -.
.-._ --__ --*-_____--------__________
w W/s1 Fig. 6. Frequency responses of the models, a: input delay system C,(jw) exp(-in,h\ G,( jw)D(joh)/F(jwh). c: tank series approximation G,( jo)/( 1 + j&/7)‘.
system
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
114
timew Fig. 7. Input and output signals of the system and approximate model. u: system input. y: output of the noise corrupted system. y,: output of the approximate model with estimated delay.
serves as the variance on the estimated delay. This is due to the fact that the quantity, say, K(t) = M(t)-‘Q,, in which Qc denotes the continuous-time counterpart of the effective variance of the disturbance, stands for the variance of the estimated delay. The discrete-time and continuous-time variances of the disturbance are related by (see C121)
t&=&T.
(48)
Crude calculations based on numerical values T = 0.1, Qd = 02 = 0.01 give a bound of the standard deviation oh of the input delay a,, < dQ,/M(t)ti,
x 0.03.
(49)
This result is in a good agreement with the simulation results of Fig. 8(b). The result is obtained by using the analogy of the overall identifier equations (24)-(25) with the equations of the linear least-squares parameter estimator scaled by the variance of the measurement noise. The details are omitted. Estimation of the bias, which is clearly visible in Fig. 8(b), would need derivation of the identifier equation as a stochastic differential equation. Then, due to the nonquadratic form of the integrand in the cost functional, a bias term appears in the identifier equation (24). These calculations are a subject of the future research. The amount of the bias depends on the forgetting exponent, as it can be seen in Fig. 8(b) in the three representative constant delay cases.
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
115
-l
,Oi 25-
I--\ I I I
. I I
! ,---__!
m15 -
\a=O.O
(4
true delay
_- -.__- __-._._ _.__.___ __._._._ __. ___.___.__,____.____ -.__- _- _.-- _-.-.-.-_-.- -
u)
t
r,?;
.. .
04
/
a=O.O4 ,a?!2
time (s)
Fig. 8. (a) Time delay estimates for three different forgetting exponents in the case of a constant delay. (b) Delay estimates in the enlarged scale.
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
116
0
50
100
150
200
250
300
350
400
time 6) Fig. 9. Gain functions
M(t)-’
of the estimator
for three forgetting
exponents.
In the case of the time-varying delay upwards and downwards ramp changes in the delay were simulated for several values of the slope of the ramp. The initial values of the estimator were ho = 1 s, and M(0) = J = 100. Two representative test runs for several different values of the forgetting exponent are depicted in Figs. 10 and 11.
5. Concluding remarks Justification of the application of the collocation method from the viewpoint of the gain approximation, as compared e.g. with a tank series model with equal number of state variables, is concluded from the frequency response curves of the different models. Due to the nonlinear character of the input delay a nonlinear structure was also obtained for the rational approximation of the delay by the collocation method. A useful feature of the approximate model is that it depends only on a single parameter: the original input delay. Then a relatively simple, even if approximate, structure were found for the identifier. In principle, the parameter-polynomial identification method developed in [ 161 and applied in [7,8] could have been used also in the new collocated model approximation. A gradient-based approximate identification procedure was, however, applied to simplify the identification task. A quadratic integral of the model error was applied as a cost functional including exponential forgetting. Then a real-time identifier was obtained. A second order system including
M. Nihtilii et al. / Simulation Pructice and Theory 5 11997) 101-120
117
c_ a=O.O3
?._ 9;
T --..,: 84 I:,____
s:
A._ __b.,
_
a=O.O4
0 time (s) Fig. 10. Delay estimates variation is 0.333).
for a time-varying
delay
for five different
forgetting
exponents
(slope in the
the input delay was applied in the simulations. Here, the rational part of the transfer function was assumed known. Simulation results presented here confirm the feasibility and applicability of the methodology developed. Furthermore, approximate error analysis could have been carried out based on the analogy between our identifier and the ordinary recursive linear least-squares parameter estimator. Crude theoretical calculations on the standard deviation of the delay estimation error are in accordance with the simulation results. Some future steps in the study of identification of time delay systems deserve to be mentioned. Most often the rational part of the transfer function is also unknown. Then, identification of the coefficients of A and B polynomials complicates the task. Furthermore, the parameter-polynomial method will be applied for the pure delay and for the joint delay and A and B polynomials. This will need application of some ready-made packages to facilitate the overall identification task. Comparisons with some other delay estimation methods are also needed to see advantages and disadvantages of the methodology studied and applied here.
Acknowledgements
This work was supported by the Research Council for Technology of the Academy of Finland and by the Centre National de la Recherche Scientifique, France, in the form of exchange fellowships, which is gratefully acknowledged.
118
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
Fig. 11. Delay estimates for a time-varying delay for three different forgetting exponents (slope in the variation is 0.1).
Appendix. Orthogonal collocation method The PDE system (9)-(10) is approximated by an ODE system via a spatial discretization in the variable z. The method is called orthogonal collocation [23]. The solution w(z, t) is approximated by a finite series Nfl w(z~
t, =
C
Bi(z)wi(t)*
64.1)
i=O
In this separation-of-variables series the functions Bi are called base functions. This series is considered at the collocation points zi: *a* <~N
O=Z~
64.2)
The base functions must have the property 1 for i =j, Pitzj)
= {
0
fori#j,
i,j=O,
1) . ..) N + 1.
64.3)
Consequently, the series gives at the collocation points w(zi,t)=wi(t),
i=O,l,...,
N+l.
(A-4)
The internal collocation points Zi, i = 1, . . . . N are obtained as zeros of the Nth order Jacobi polynomial denoted by PN (‘J) , that is defined to satisfy the following set of
M. Nihtilii et al. / Simulation
orthogonality
Practice and Theory 5 (1997) 101-120
119
conditions
1 zD(
1 - z)“zkP$*B)(z)dz=O,
k=O,l,...,
N-l.
(A.5)
s0
The powers CIand fi are tuning parameters of the method. They affect the location of the internal collocation points. An interpolation polynomial L(z) is obtained by using the Jacobi polynomial and two additional zeros z = 0 and z = 1 in the following way L(z) = z(z - l)P$,B’(z).
(A.6)
Lagrange interpolation is then applied to join the values W(Zi,t) together. This means that the approximate solution is defined by the series N+l w(z,
t,
1:
1
Pi(z)w(zi,
(A.7)
t).
i=O
The base functions are defined by Bitz)
=
L(z)IC(z
-
zi)L’(zi)l,
(A.8)
in which the prime ’denotes differentiation with respect to z. It can be seen that the base functions have the desired property (A.3) which results in the identity (A.4). The series (A.7) is substituted into the PDE model and it is considered at all collocation points but the initial one. The point z. is omitted because of the boundary condition (10). The substitution means that the derivatives N+l wt(zj,
r,
=
1
(A.9)
Bi(zj)tii(t),
i=O N+l wz(zj,
r,= 1
(A.lO)
Bl(zj)wi(r),
i=O
are applied in equation (9). By defining a state vector
dr) = [wl(r) WZ(t)“’ wN+l(r)lT
q
by (A.1 1)
the state-variable equation (15)-( 16) with the definitions (17)-( 19) are obtained.
References Cl1 M.
Agarwal and C. Canudas, On-line estimation of time-delay and continuous-time process parameters, Internat. J. Control 46(l) (1987) 2955311. and control of time-delay systems, Licentiate Thesis, Rept. c21 M. Bjiirklund, Identification UPTEC8884R, Institute of Technology, Uppsala University, Sweden, 1988. of linear systems c31 M. Bjiirklund, M. Nihtila and T. Siiderstrbm, An algorithm for identification with varying time delay, in: Preprints 9th IFACIIFORS Symp. on Identification and System Parameter Estimation, Budapest, Hungary (1991) 1254-1259. Modelling and adaptive control of non-linear c41 D. Dochain, J.-P. Babary and N. Tali-Maamar, distributed parameter bioreactors via orthogonal collocation, Automatica 28(9) (1992) 873-883.
120
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
[S] G.F. Franklin, J.D. Powell and A. Emami-Naeini, Feedback Control ofDynamical Systems (AddisonWesley, Reading, MA, 1986). [6] T. Furukawa and E. Shimemura, Predictive control for systems with time delay, Internat. J. Control 37(2) (1986) 399-412. [7]
P. Gawthrop and M. Nihtill, Identification of time delays using a polynomial identification method, Systems Control Lett. S(2) (1985) 267-271. [S] P. Gawthrop, M. Nihtilti and A. Besharati-Rad, Recursive parameter estimation of continuoustime systems with unknown time delay, Control-Theory and Advanced Technology 5(3) (1989) 227-248. [9] L.A. Gould, Chemical Process Control: Theory and Applications (Addison-Wesley, Reading, MA, 1969). [ 10) T. HBgglund, A predictive PI controller for processes with long dead times, IEEE Control Systems Mag. 12( 1) (1992) 57-60. [ 111 K. Ichikawa, Frequency-domain pole assignment and exact model matching for delay systems, Internat. J. Control 41(4) (1985) 1015-1024. [12] A.H. Jazwinski, Stochastic Processes and Nonlinear Filtering Theory (Academic Press, New
York, 1970) [13] A.Z. Manitius and A.W. Olbrot, Finite spectrum assignment problem for systems with delays, IEEE Trans. Automatic Control 24(4) (1979) 541-553. [14]
[ 151
[16] [17]
[ 181 [ 191 [20]
[21] [22] [23]
[24]
M. Molander, Computer aided modelling of distributed parameter processes, Ph.D. Thesis, Tech. Rept. No. 193, School of Electrical and Computer Engineering, Chalmers University of Technology, Sweden, 1988. P.A.J. Nagy and L. Ljung, Estimating time-delays via state-space identification methods, in: Preprints 9th IFACIIFORS Symp. on Identification and System Parameter Estimation, Budapest, Hungary (1991) 1141-1144. M. NihtilB, Optimal finite dimensional recursive identification in a polynomial output mapping class, Systems Control Lett. 3(6) (1983) 341-348. M. Nihtilti, N. Tali-Maamar and J.-P. Babary, Controller design for a nonlinear distributed parameter system via a collocation method, in: Proc. IFAC Symp. on Nonlinear Control System Design, NOLCOS’92, Bordeaux, France (1992) 447-452. H.M. Paynter, The Paynter delay line, Lightning Empiricist 11 (1963) 10-11. A.E. Pearson and C.Y. Wuu, Decoupled delay estimation in the identification of differential delay systems, Automatica 20(6) (1984) 761-772. W.R. Robinson and A.C. Soudiak, A method for the identification of time delays in linear systems, IEEE Trans. Automatic Control 15( 1) (1970) 97-101. O.J.M. Smith, Feedback Control System (McGraw-Hill, New York, 1958). N. Tali-Maamar, T. Damak, J.P. Babary and M.T. NihtilB, Application of a collocation method for simulation of distributed parameter bioreactors, Math. Comput. Simulation 35(4) (1993) 303-319. J. Villadsen and M.L. Michelsen, Solution of Difirential Equation Models by Polynomial Approximation (Prentice-Hall, London, 1978). K. Zenger, Analysis and control design of a class of time-varying systems, Licentiate Thesis, Rept. 88, Control Engineering Laboratory, Helsinki University of Technology, Finland, 1992.
Practice
and Theory
= THEORY
5 (1997) lOlL120
On-line estimation of the time delay via orthogonal collocation M. NihtilZ a,*, T. Damak b. J.P. Babary ’ a University of Kuopio, Faculty of Natural and Environmental Sciences, Department of Computer Science and Applied Mathematics, P.O. Box 1627, FIN-7021 1 Kuopio, Finland b Ecole Nationale d’lngtnieurs de Sfax. Department de Genie Micanique, B.P.W. 3038 Sfax, Tunisie ’ Centre National de la Recherche Scienti$que, Laboratoire d’Analyse et d’Architecture des Syst?mes, 7 avenue du Colonel Roche, F-31077 Toulouse Cedex, France Received
12 February
1995, revised 11 November
1995
Abstract The objective of this work is to develop a new real-time time-delay estimator for the input delay of the otherwise finite-dimensional single-input single-output linear differential system. The method applies also to tracking of possible time-variations in the delay. The key idea in the construction of the estimator is to transform the delay part of the system into a transport system described by a linear partial differential equation (PDE). The PDE is linearly parametrized by the delay argument. The PDE is approximated by a finite-dimensional ordinary differential equation (ODE) system. The approximation is obtained by applying orthogonal collocation in the spatial descretization of the PDE. A fixed-interval exponentially weighted cost functional is then applied in the ODE system in formulating the delay estimation problem in mathematical terms. The data forgetting issue allows tracking of changes in the delay. Minimization of the cost functional is carried out approximately to obtain a finite-dimensional recursive delay estimator for increasing values of the final time of the integral in the cost. Simulations for a fixed as well as for a time-varying delay are given to demonstrate estimation and tracking properties of the identifier. Key words: simulation
Delay
estimation;
Least-squares
identification;
Collocation
method;
Model
1. Introduction The destabilizing effect of time delays is a known fact in closed-loop control the delays are ignored in the control law calculations, see e.g. [lo]. Different * Corresponding
author.
Email: [email protected].
0928-4869/97/$17.00 Copyright simpraOO82XX PZZ 0928-4869(95)00052-6
0 1997 - Elsevier Science B.V. All rights reserved
when types
102
M. Nihtilii et al. J Simulation Practice and Theory 5 (1997) 101-120
of predictive control algorithms have been developed to overcome these problems by taking into account the system delays in the controller design. Stable systems can be controlled by applying Smith’s predictor [5,21]. Unstable systems need so-called analytic predictors or finite spectrum assignment controllers, which have been studied among others by Manitius and Olbrot [13], Furukawa and Shimemura [6], and Ichikawa [ 111. Predictive methods result in closed loop responses whose performance is comparable with those of the controlled systems without delays. Most predictive methods directly utilize numerical values of the delays. Consequently, it is necessary to know, estimate or track on-line values of the delay to achieve sufficiently good closed-loop performance. The transcendental nature of the transfer functions of linear input delay systems causes problems in the application of ordinary control techniques. However, there are several methods to approximate the transcendental transfer function of the delay by rational ones as it has been described e.g. by Agarwal and Canudas [ 11, Gawthorp and Nihtila [7], and Paynter [ 181. Then these approximations are parametrized by some functions of the delay argument. In using rational approximations of the delay it is essential from the identification viewpoint that the approximations do not include too many free parameters. PadC approximations, tank series models and allpass filters are among the best known rational approximations. These can be obtained by applying frequency-response techniques. For discrete-time models the inclusion of a (variable) delay increases generally the number of free parameters. In some special cases the delay approximation can be based on single parameter models also in a discrete formulation as presented by Bjorklund et al. [ 2,3]. The goal of this paper is to construct an estimator for the input delay of a linear system based on given measurements of the input and the output. Furthermore, the identification problem should be recursively solvable to be applicable in real-time purposes. Due to the nonlinear character of the effect of the input delay, it can be expected that only finite-dimensional recursive estimators, which are approximate with respect to some feasible criteria, can be constructed. The identification of time-delays is not an easy task. Some different approaches can, however, be found. Among the first studies describing on-line identification of a variable delay is the one presented in [20]. Pearson and Wuu [19] applied a decomposition technique in the off-line environment to identify a constant delay separately from the other system parameters. Several recursive variants, applying delay approximation or direct gradient-based approximate recursion, are studied by Gawthorp et al. [S]. Modern state-space identification packages are also applicable in delay identification as described recently in [ 151. In this paper an alternative approach for recursive input delay identification is presented. It is based on a specific rational approximation of the delay. The input delay is considered as a transport delay. Then it is exactly described by a linear firstorder partial differential equation (PDE) including the delay as a single parameter. This PDE model is approximated by an ordinary differential equation (ODE) system. The approximation procedure is based on the method of orthogonal collocation, see e.g. [23]. Specific applications of this approximation method have been studied
M. Nihtilii et al. 1 Simulation
103
Practice and Theory 5 (1997) 101-120
among others by Molander [ 141, Dochain et al. [4], and Nihtila et al. [ 171. In the method proposed here the resulting ODE system is parametrized in a polynomial way by the delay. An approximate recursive least-squares method in continuoustime form is then applied in the delay estimation. The efficiency of any numerical method, in our case recursive identification of the time-delay, can be evaluated either via the theoretical analysis of the algorithms at hand or via simulations for some specific parameter values. Due to the nonlinear character of the problem we have chosen the simulation way. The simulations actually show the feasibility of the technique developed. The organization of this paper is as follows. The second section treats the system models and introduces the procedure to approximate the delay system by a rational transfer function parametrized by the single delay value. In the third section the identification problem is stated in mathematical terms as an optimization problem. Next the approximate on-line identifier is constructed. In the fourth section a secondorder underdamped system with input delay is studied via simulations to test the algorithm and its feasibility. Constant as well as time-varying delay cases are examined. Some concluding remarks complete the paper. In the Appendix the orthogonal collocation method is revisited with some details.
2. System models The following is considered: 4P)Y(t)
differential
= &+4t
A and B are polynomials given by
-
system with a constant
(or time-varying)
input
(1)
h).
of known
A(p)=a,+a,p+a,p2+
delay h
degrees
in the derivative
operator
p. They are
... +a$“-‘+p”,
B(p) = b, + b2p + b3p2 + ... + bmpm-l,
(2) m d n,
(3)
in which IZ and m - 1 are the degrees of A and B, respectively. In transfer function formalism the system is depicted in Fig. 1. A state variable model corresponding to
u
>
,-b
X
>
B(P) A(P)
Fig. 1. Original
input delay system
Y
)
M. Nihtiliiet al. / SimulationPracticeand Theory 5 (1997) 101-120
104
the system (1) is in canonical observable form given by i(t) = A,x(t) + B,u(t - h),
(4)
y(t) = Cx(t),
(5)
in which the system matrices of appropriate 0
A, =
0
...
0
1 0
...
0
0
f.. ..
0 :
-a2 -a3
1
-a,
1
OO...
c=
dimensions are given by
-a,
’
B, =
[O ... 0 11.
(6 7)
(8)
The first step in system identification considered is to estimate the delay in the case of a known rational part, i.e. G,(p) = B( p)/A(p), of the transfer function. In order to be able to apply input-output techniques without inverting G,(p) the order of the transfer functions is reversed in Fig. 1. If the two parts of the system are timeinvariant then the order can be reversed without affecting the input-output behaviour of the system. This is given in Fig. 2. 2.1. Delay approximation In the framework of flow processes, studied thoroughly e.g. by Zenger [24], the delay part of the system in Fig. 2 is modelled as a distributed parameter system described by the following partial differential equation with a relevant boundary condition and the given output equation (cf. [9]) Wt(Z,t) = -( l/h)w,(z, t),
0 < z < 1,
(9)
NO, t) = v(t),
(10)
Y(L)= w( 1, 0.
(11)
The subindices t and z denote partial derivatives with respect to the time and space
bp+j-iqL> Fig. 2. Reversed order model.
M. Nihtilii et al. J Simulation Practice and Theory 5 (1997) 101-120
105
variables, respectively. The input v of the delay system is obtained from the ordinary linear system 4pMt)
= NPMO.
(12)
It must be emphasized that the system (9)-( 11) is an equivalent representation of the delay part of Fig. 2. Furthermore, w(z, t) corresponds to v(t - zh), and ~(0, t) (and w,(t) in (13)) to v(t). The orthogonal collocation method, cf. [23], which will be applied in the PDE system (9)-( 11) is a spatial discretization method. When applied, it converts the PDE system into an ordinary differential equation system which is an approximate input-output description of the delay part of Fig. 2. The details are given in the Appendix. By applying the series N+l w(z,
t) =
1
BdzMt),
(13)
i=O
with appropriately chosen base functions pi in the PDE system the following finitedimensional approximate state variable system is obtained for the vector 4, where
k(t) = F,q(t) + 4v(t),
(15)
Yd(C h) =
(16)
Hdt).
A new notation yd(t; h) is introduced to emphasize that it is an approximation depending on the delay. The following expressions for the elements of the matrices are obtained (F,)ij=
-p(i(Zi), i,j= 1, . . . . N + 1,
(O1)i = -p;(zi), H=[O
i = 1, . . . . N + 1,
... 0 11,
(17) (18) (19)
in which the prime ’stands for derivation with respect to the space variable z. The points zo, zl, . . . . zN+i are appropriately chosen successive points on the spatial interval [0, 11. N is the number of internal points in the discretization. Based on the specific position of the delay in (15), i.e. the derivative operator p is everywhere multiplied by the delay h, the following unique polynomial operator representation WP)YdC
h) = WPMt)
(20)
is obtained for the delay approximation from the state-variable representation (15)-( 16). The coefficients of the polynomials F(s) and D(S) corresponding to the matrices F, and D, depend only on the tuning parameters of the collocation method a and j3 and on the number of the collocation points N (cf. the Appendix). The
M. Nihtilli et al. / Simulation Practice and Theory 5 (1997) 101-120
106
polynomials do not depend on the specific value of the delay. Consequently, the system (20) is a rational single parameter approximation of the delay system. A practical point is also that the single parameter is the delay itself.
3. Least-squares identification Minimization of the error between the measured and model generated output, i.e. of y(t) - yd(t; h), with respect to the delay h is formulated as an exponentially weighted least-squares identification problem. Then, without stochastic concepts, measure-theoretic problems arising in stochastic time-continuous formalism can be avoided. In the following the formulation of the problem and an approximate recursive solution will be given. 3.1. Problem statement The identification of the delay is stated in mathematical terms as a set of fixedinterval optimization problems. Consequently, for each fixed final time t 2 0 solve the following task. Fixed interval optimization problem. For the given input v and the output y on the interval [0, t] find the value of the delay, denoted by h(t), which minimizes the cost
functional t V(t,h)=~e-“‘J(h-h,)‘+~
ea(s-r)[y(s) - y,,(s; h)]’ ds, s0
(21)
subject to the constraint Us) = F14(s) + 444, Y&
4 = Rl(d.
4(O) = 40,
(22) (23)
for 0 < s < t. The given scalars a, ho, and J denote the forgetting exponent, the initial delay approximation and the (positive) weight of the initial cost, respectively. If the problem statement were given by using filtered signals, then a parameterpolynomial identification scheme could have been used [ 161. This is due to the fact that the operator polynomials F(hp) and D(hp) are polynomials also in the parameter h. Then the filtered model will be in a polynomial form with respect to the delay amenable for the parameter-polynomial identification. The polynomial identifier is optimal in the sense that a finite-dimensional recursive algorithm already gives the optimal, not approximate, least-squares fit in the approximate filtered model. Due to relative complexity of the parameter-polynomial method, however, here an approximate algorithm is used. This approximate algorithm is not optimal in the sense mentioned above.
M. Nihtilii et al. / Simulation
Practice and Theory 5 (1997) 101-120
107
3.2. Approximate recursion for the delay The gradient-based recursion for the delay estimate i(t) is given by the two basic scalar differential equations (see derivation in Section 3.4) h(t) = M(t)_‘$(t; &l(t) = -aM(t)
&))CYW - y,(t; + $(t; L(t))‘,
fiW)ls
i(O) = h,,
M(0) = .I.
(24) (25)
The function rl/ denotes the partial derivative or gradient of the model output, i.e.
$(t; h) =
4 ydt; h).
3.3. Gradient system Calculation of the gradient is based on the following procedure. The approximate system represented in the transfer function form (27)
Wp)Yd(t; h) = N&(t) is differentiated with respect to the delay h giving the system
(28)
F,,(hp)YJt; h) + Wp)ll/(t; h) = Q,(hp)o(t). The subindex h denotes differentiation with respect to the delay. It is seen that Uhp) = M”(hp),
(29)
&(hp) = @VP),
(30)
in which the prime ’indicates differentiation
of the polynomials F and G, i.e.
df’b) F’(s)= y, D’(s) =
dW
(32)
ds.
If the number of the internal collocation (i.e. spatial discretization) points, when approximating the system (9)-( 11) by the system (15)-( 16), is N, then the gradient system (28) can be written in the form +N+
1 yd(t; h), I
(33)
in which the following definitions are applied,
QW=
W'Vvh
R(hp) = hpF’(hp) -(N
(34) + l)F(hp).
(35)
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
108
The actual realization of the gradient is obtained by using two state-variable filters
s,: 4 = FIP + !A% p(O)= 0,
(36)
t=HP, S,:
w(O)=O,
hti=F,o+Rly,,
(37)
i=Ho+(N+l)y,, with the joint output equation
(38) in which p and o are (N + 1)-vectors. The filters S1 and SZ correspond to the transfer functions of the coefficients of v(t) and yd(t; h) of (33) (without the scalar coefficient l/h, which is included in the output equation (38)). The vectors Q, and R, in (36) and (37) correspond to realizations of the Q(s) and R(s) polynomials (cf. the corresponding realization of ( 12) by the matrices (6)-( 8)) in the respective transfer functions Q(s)/F(s) and R(s)/P(s). Consequently, they are independent of the delay. It has to be noted that, in estimating the delay via the system (24)-(25), the current value of the estimate h(t) must be used in the state variable filters and in the joint output equation. 3.4. Derivation of the identijier The minimizing condition for the cost (21) is that its gradient is equal to zero for each final time t, i.e. I/h(t, h) = e-“‘J(h
- h,) -
f J
e”(S-f)$(s; h)[y(s) - yd(s; h)] ds = 0,
(39)
0
for the optimal estimate h(t). Consequently, we have for all t E [0, T] the identity qt,
(40)
h(t)) = 0.
Total differentiation of the gradient condition with respect to the final time t (considered now as a variable) gives the condition -a&At, i(t)) - $(t;
&MY(t)
-Yy,(C
&t))l
f
- Jf
+ e-“‘J+ e”‘“-‘)J/(s;L(t))’ds i s0 eU(S-‘)tj,(s; h(t))[y(s)
0
- y,(s; h(t))] ds k(t)) = 0,
(41)
Al. Nihtilii et al. J Simulation Practice and Theory 5 (1997) 101-120
109
Introduce then the following scalar functions f M(t)=
eCatJ
eacspt)$(s; A(s))~ds,
+
(42)
s0 f N(t)=
e-a'J
ea(s-t)$(s; h(t))' ds
+
f -1
c0
ea(s-t)$h(s;
h(t))[y(s)-
yd(s; 6(t))] ds.
(43)
0
Then the truly optimal estimator is given by (43) and (24) in which the gain M(t)- ’ is substituted by N(t)-’ assumed that it exists. Because the optimal inverse gain (43) is not recursively computable by a finite-dimensional ODE system but only in some special cases it is approximated by the suboptimal inverse gain (42) which satisfies the differential equation and the initial condition in (25). The approximation of N(t) by M(t) needs two simplifications. The last term in (43) is actually zero if the residual, or the output modelling error y(s) - y,(s; g(t)) is uncorrelated with $(s, G(t)) for 0 6 s < t. Furthermore, if the estimated delay value changes sufficiently slowly on the effective interval of the exponential weighting function exp[a(s - t)] then the estimated delay k(t) in (43) can be substituted by A(s).This second assumption finally results in the inverse gain equation (42). The initial conditions of the estimator are due to the initial cost in (21). In the derivation it is assumed that h is a differentiable function of its argument t. This includes that the cost has such a unique global minimum, with respect to the delay, which does not jump between possible local minima, when the time t increases from 0 to T. This issue is considered in connection with a parameter-polynomial cost functional in [ 161.
4. Simulation study 4.1. Model approximation For the delay approximation feasible values of the tuning parameters a, p and N have to be chosen. The selection of CIand p is based on a trial-and-error technique. A feasible value of N is obtained related to the accuracy of the approximation required on the spatial domain of the PDE system (9)-( 1l), which is equivalent with the delay part of the original model (cf. [22]). The specific process model studied here is given by G(p) = Go(P)eehp,
(44)
4 G,(P)
=
p2 + 2j0,p
+ w;
A step response of the approximation,
(45)
i.e. of Go and the system (15)-( 16) is visually
110
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
fitted with the step response of G by varying the tuning parameters. This was carried out in the noiseless case. The values of the system and tuning parameters are given in Table 1. In all simulations outputs of the system model, the estimator and the state-variable filters were generated at discrete intervals of length Td = 0.1 s. A fourth-order RungeKutta method with the integration step size z = 0.01 s was applied to generate the signals. The responses of the tuned approximation and the second-order model are depicted in Fig. 3 and are almost identical from the visual viewpoint. The difference y - y,(.; h) is represented in Fig. 4 but in a different amplitude scale. The coefficients of the approximation polynomials F(s) =fr +fzs + .** +f,P+
s7,
(46)
D(s) = d, + d,s + . . . + d,@, are given in Table 2. The zeros (and the discretization
(47) points) of the interpolation
polynomial L(z) (cf.
Table 1 Values of the system and tuning parameters
wo i h
0.2 rad/s 0.5 radfs 20 s
N
Fig. 3. Step responses of the input delay (-)
6 1 1
and the approximate model (- - - -).
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
111
0.015 -
0.01 -
0.005 -
o-
-0.005 -
-0.01 -
-0.015 0
20
40
60
80
time (s) Fig. 4. Difference
Table 2 Coefficients
of the polynomials
of the step responses
of Fig. 3.
F(s) and D(s)
2.1622 x lo6 1.2355 x lo6 5.544 x xlo4 3.3264 lo5 6.3 x lo3 504 28
2.1622 x lo6 - 9.2664 x 10’ - 1.782 1.98 x xlo4 lo5
d, d, d, d, d, d,
1.35 x 103 -54 1
the Appendix) are depicted in Fig. 5. Due to the tuning values CI= /? = 1 the points are symmetrically located around z = 0.5. The matrix F, and vector D, of the statevariable system are given in Table 3 in the form [F,ID,]. Visual comparison of frequency responses of the delay model with various approx-
z()=
0
z,=
Fig. 5. Location
of collocation
points
for a = p = 1.
1
112
M. Nihtilii et al. j Simulation Practice and Theory 5 (1997) 101-120
Table 3 Elements of the enlarged matrix [F, -1.0 x 1O-6 5.6130 -2.1579 1.3223 - 1.0741 1.1471 -2.5954
-9.0872 0 4.7564 - 2.2707 1.6900 - 1.7389 3.8833
1Ill]
4.2241 -5.7510 1.0 x 1O-6 4.7779 -2.7456 2.5885 - 5.6204
- 2.5885 2.7456 -4.7778 0 5.7510 -4.2241 8.5958
1.7389 - 1.6900 2.2707 - 4.7564 2.0 x 1O-6 9.0872 -15.139
- 1.1471 1.0741 - 1.3223 2.1579 -5.6130 -7.0 x 1o-6 37.875
0.4399 -0.4066 0.4867 -0.7443 1.5850 - 6.4198 -28.0
6.4198 - 1.5850 0.7443 -0.4867 0.4066 -0.4399 1.0
imations gives an insight of the accuracy of the approximations. Here, in addition to the collocated approximation, a tank series approximation of the delay system is taken for comparisons. Frequency responses of the second-order delay model (44)-(45), the collocated model (12) with the signal u generated by Go and the tank series model of the delay with Go are given in Fig. 6. The number of the tanks has been chosen equal to the number of the state variables in the collocated model. This gives the fair basis for comparing the performance of the models on the frequency band considered. The tuning values chosen for the collocation method emphasize accurate approximation of the system gain and phase on the frequency band 0
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
1.2
IGO’V 41 -r
T
1711
-
T
T
113
-
R
..,.... /
.
:
1
0.8
I. ,.j .,...., ..,
c
,
”
/
..
1.
;
:. .:: : :
;;
0.6
..i
i. :.:
.i
:
;
:
.;
‘: i
0.4
“”
:
;
:
1: . :.
.:
j.i.:_
..:. ;
.j.
:
: : : . .
‘:’
:
‘:
^
.
.I
.
^.
.
A.,
.:
.:.
,..
0.2
I
I
‘i. :
:.:.:
.:.:.:
:
:
: :
:
.:
:
;.:.
.__
:
i .:.
d PO-2
i
i
:
:: :
+_ -;
,”
i
ii
10”
AL
1
argV3.W) Meg1
-200 -
-KJO-
-loo0 ‘S.-b -.
.-._ --__ --*-_____--------__________
w W/s1 Fig. 6. Frequency responses of the models, a: input delay system C,(jw) exp(-in,h\ G,( jw)D(joh)/F(jwh). c: tank series approximation G,( jo)/( 1 + j&/7)‘.
system
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
114
timew Fig. 7. Input and output signals of the system and approximate model. u: system input. y: output of the noise corrupted system. y,: output of the approximate model with estimated delay.
serves as the variance on the estimated delay. This is due to the fact that the quantity, say, K(t) = M(t)-‘Q,, in which Qc denotes the continuous-time counterpart of the effective variance of the disturbance, stands for the variance of the estimated delay. The discrete-time and continuous-time variances of the disturbance are related by (see C121)
t&=&T.
(48)
Crude calculations based on numerical values T = 0.1, Qd = 02 = 0.01 give a bound of the standard deviation oh of the input delay a,, < dQ,/M(t)ti,
x 0.03.
(49)
This result is in a good agreement with the simulation results of Fig. 8(b). The result is obtained by using the analogy of the overall identifier equations (24)-(25) with the equations of the linear least-squares parameter estimator scaled by the variance of the measurement noise. The details are omitted. Estimation of the bias, which is clearly visible in Fig. 8(b), would need derivation of the identifier equation as a stochastic differential equation. Then, due to the nonquadratic form of the integrand in the cost functional, a bias term appears in the identifier equation (24). These calculations are a subject of the future research. The amount of the bias depends on the forgetting exponent, as it can be seen in Fig. 8(b) in the three representative constant delay cases.
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
115
-l
,Oi 25-
I--\ I I I
. I I
! ,---__!
m15 -
\a=O.O
(4
true delay
_- -.__- __-._._ _.__.___ __._._._ __. ___.___.__,____.____ -.__- _- _.-- _-.-.-.-_-.- -
u)
t
r,?;
.. .
04
/
a=O.O4 ,a?!2
time (s)
Fig. 8. (a) Time delay estimates for three different forgetting exponents in the case of a constant delay. (b) Delay estimates in the enlarged scale.
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
116
0
50
100
150
200
250
300
350
400
time 6) Fig. 9. Gain functions
M(t)-’
of the estimator
for three forgetting
exponents.
In the case of the time-varying delay upwards and downwards ramp changes in the delay were simulated for several values of the slope of the ramp. The initial values of the estimator were ho = 1 s, and M(0) = J = 100. Two representative test runs for several different values of the forgetting exponent are depicted in Figs. 10 and 11.
5. Concluding remarks Justification of the application of the collocation method from the viewpoint of the gain approximation, as compared e.g. with a tank series model with equal number of state variables, is concluded from the frequency response curves of the different models. Due to the nonlinear character of the input delay a nonlinear structure was also obtained for the rational approximation of the delay by the collocation method. A useful feature of the approximate model is that it depends only on a single parameter: the original input delay. Then a relatively simple, even if approximate, structure were found for the identifier. In principle, the parameter-polynomial identification method developed in [ 161 and applied in [7,8] could have been used also in the new collocated model approximation. A gradient-based approximate identification procedure was, however, applied to simplify the identification task. A quadratic integral of the model error was applied as a cost functional including exponential forgetting. Then a real-time identifier was obtained. A second order system including
M. Nihtilii et al. / Simulation Pructice and Theory 5 11997) 101-120
117
c_ a=O.O3
?._ 9;
T --..,: 84 I:,____
s:
A._ __b.,
_
a=O.O4
0 time (s) Fig. 10. Delay estimates variation is 0.333).
for a time-varying
delay
for five different
forgetting
exponents
(slope in the
the input delay was applied in the simulations. Here, the rational part of the transfer function was assumed known. Simulation results presented here confirm the feasibility and applicability of the methodology developed. Furthermore, approximate error analysis could have been carried out based on the analogy between our identifier and the ordinary recursive linear least-squares parameter estimator. Crude theoretical calculations on the standard deviation of the delay estimation error are in accordance with the simulation results. Some future steps in the study of identification of time delay systems deserve to be mentioned. Most often the rational part of the transfer function is also unknown. Then, identification of the coefficients of A and B polynomials complicates the task. Furthermore, the parameter-polynomial method will be applied for the pure delay and for the joint delay and A and B polynomials. This will need application of some ready-made packages to facilitate the overall identification task. Comparisons with some other delay estimation methods are also needed to see advantages and disadvantages of the methodology studied and applied here.
Acknowledgements
This work was supported by the Research Council for Technology of the Academy of Finland and by the Centre National de la Recherche Scientifique, France, in the form of exchange fellowships, which is gratefully acknowledged.
118
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
Fig. 11. Delay estimates for a time-varying delay for three different forgetting exponents (slope in the variation is 0.1).
Appendix. Orthogonal collocation method The PDE system (9)-(10) is approximated by an ODE system via a spatial discretization in the variable z. The method is called orthogonal collocation [23]. The solution w(z, t) is approximated by a finite series Nfl w(z~
t, =
C
Bi(z)wi(t)*
64.1)
i=O
In this separation-of-variables series the functions Bi are called base functions. This series is considered at the collocation points zi: *a* <~N
O=Z~
64.2)
The base functions must have the property 1 for i =j, Pitzj)
= {
0
fori#j,
i,j=O,
1) . ..) N + 1.
64.3)
Consequently, the series gives at the collocation points w(zi,t)=wi(t),
i=O,l,...,
N+l.
(A-4)
The internal collocation points Zi, i = 1, . . . . N are obtained as zeros of the Nth order Jacobi polynomial denoted by PN (‘J) , that is defined to satisfy the following set of
M. Nihtilii et al. / Simulation
orthogonality
Practice and Theory 5 (1997) 101-120
119
conditions
1 zD(
1 - z)“zkP$*B)(z)dz=O,
k=O,l,...,
N-l.
(A.5)
s0
The powers CIand fi are tuning parameters of the method. They affect the location of the internal collocation points. An interpolation polynomial L(z) is obtained by using the Jacobi polynomial and two additional zeros z = 0 and z = 1 in the following way L(z) = z(z - l)P$,B’(z).
(A.6)
Lagrange interpolation is then applied to join the values W(Zi,t) together. This means that the approximate solution is defined by the series N+l w(z,
t,
1:
1
Pi(z)w(zi,
(A.7)
t).
i=O
The base functions are defined by Bitz)
=
L(z)IC(z
-
zi)L’(zi)l,
(A.8)
in which the prime ’denotes differentiation with respect to z. It can be seen that the base functions have the desired property (A.3) which results in the identity (A.4). The series (A.7) is substituted into the PDE model and it is considered at all collocation points but the initial one. The point z. is omitted because of the boundary condition (10). The substitution means that the derivatives N+l wt(zj,
r,
=
1
(A.9)
Bi(zj)tii(t),
i=O N+l wz(zj,
r,= 1
(A.lO)
Bl(zj)wi(r),
i=O
are applied in equation (9). By defining a state vector
dr) = [wl(r) WZ(t)“’ wN+l(r)lT
q
by (A.1 1)
the state-variable equation (15)-( 16) with the definitions (17)-( 19) are obtained.
References Cl1 M.
Agarwal and C. Canudas, On-line estimation of time-delay and continuous-time process parameters, Internat. J. Control 46(l) (1987) 2955311. and control of time-delay systems, Licentiate Thesis, Rept. c21 M. Bjiirklund, Identification UPTEC8884R, Institute of Technology, Uppsala University, Sweden, 1988. of linear systems c31 M. Bjiirklund, M. Nihtila and T. Siiderstrbm, An algorithm for identification with varying time delay, in: Preprints 9th IFACIIFORS Symp. on Identification and System Parameter Estimation, Budapest, Hungary (1991) 1254-1259. Modelling and adaptive control of non-linear c41 D. Dochain, J.-P. Babary and N. Tali-Maamar, distributed parameter bioreactors via orthogonal collocation, Automatica 28(9) (1992) 873-883.
120
M. Nihtilii et al. / Simulation Practice and Theory 5 (1997) 101-120
[S] G.F. Franklin, J.D. Powell and A. Emami-Naeini, Feedback Control ofDynamical Systems (AddisonWesley, Reading, MA, 1986). [6] T. Furukawa and E. Shimemura, Predictive control for systems with time delay, Internat. J. Control 37(2) (1986) 399-412. [7]
P. Gawthrop and M. Nihtill, Identification of time delays using a polynomial identification method, Systems Control Lett. S(2) (1985) 267-271. [S] P. Gawthrop, M. Nihtilti and A. Besharati-Rad, Recursive parameter estimation of continuoustime systems with unknown time delay, Control-Theory and Advanced Technology 5(3) (1989) 227-248. [9] L.A. Gould, Chemical Process Control: Theory and Applications (Addison-Wesley, Reading, MA, 1969). [ 10) T. HBgglund, A predictive PI controller for processes with long dead times, IEEE Control Systems Mag. 12( 1) (1992) 57-60. [ 111 K. Ichikawa, Frequency-domain pole assignment and exact model matching for delay systems, Internat. J. Control 41(4) (1985) 1015-1024. [12] A.H. Jazwinski, Stochastic Processes and Nonlinear Filtering Theory (Academic Press, New
York, 1970) [13] A.Z. Manitius and A.W. Olbrot, Finite spectrum assignment problem for systems with delays, IEEE Trans. Automatic Control 24(4) (1979) 541-553. [14]
[ 151
[16] [17]
[ 181 [ 191 [20]
[21] [22] [23]
[24]
M. Molander, Computer aided modelling of distributed parameter processes, Ph.D. Thesis, Tech. Rept. No. 193, School of Electrical and Computer Engineering, Chalmers University of Technology, Sweden, 1988. P.A.J. Nagy and L. Ljung, Estimating time-delays via state-space identification methods, in: Preprints 9th IFACIIFORS Symp. on Identification and System Parameter Estimation, Budapest, Hungary (1991) 1141-1144. M. NihtilB, Optimal finite dimensional recursive identification in a polynomial output mapping class, Systems Control Lett. 3(6) (1983) 341-348. M. Nihtilti, N. Tali-Maamar and J.-P. Babary, Controller design for a nonlinear distributed parameter system via a collocation method, in: Proc. IFAC Symp. on Nonlinear Control System Design, NOLCOS’92, Bordeaux, France (1992) 447-452. H.M. Paynter, The Paynter delay line, Lightning Empiricist 11 (1963) 10-11. A.E. Pearson and C.Y. Wuu, Decoupled delay estimation in the identification of differential delay systems, Automatica 20(6) (1984) 761-772. W.R. Robinson and A.C. Soudiak, A method for the identification of time delays in linear systems, IEEE Trans. Automatic Control 15( 1) (1970) 97-101. O.J.M. Smith, Feedback Control System (McGraw-Hill, New York, 1958). N. Tali-Maamar, T. Damak, J.P. Babary and M.T. NihtilB, Application of a collocation method for simulation of distributed parameter bioreactors, Math. Comput. Simulation 35(4) (1993) 303-319. J. Villadsen and M.L. Michelsen, Solution of Difirential Equation Models by Polynomial Approximation (Prentice-Hall, London, 1978). K. Zenger, Analysis and control design of a class of time-varying systems, Licentiate Thesis, Rept. 88, Control Engineering Laboratory, Helsinki University of Technology, Finland, 1992.