A Bayesian approach for time-delay estimation of a managed river reach in imposed experimental conditions

Copyright © IFAC Time Delay Systems, Rocquencourt, France, 2003

ELSEVIER

IFAC PUBLICATIONS www.elsevier.comllocalelifac

A BAYESIAN APPROACH FOR TIME-DELAY ESTIMATION OF A MANAGED RIVER REACH IN IMPOSED EXPERIMENTAL CONDITIONS Magalie Thomassin • , Thierry Bastogne' ,Alain Richard' and Antoine Libaux ••

• Centre de Recherche en Automatique de Nancy (CRAN) CNRS UMR 7039. Universite Henri Poincare. Nancy I. BP 239, 54506 Vand(Euvre-Ies-Nancy Cedex. France, Phone: +33 (0) 3 83 684461 - Fax: +33 (0) 3 83 68 44 62 [email protected] http://www.cran.uhp-nancy.fr •• EDF ClH FCC, 73373 Le Bourget-du-Lac Cedex, France Phone: +33 (0) 479606389 [email protected]

Abstract: This article deals with the problem of estimating time-delays in the experimental modelling of river reaches managed by hydroelectric power plants. The purpose of the article is to estimate the reach time-delay from the data collected in imposed experimental conditions. The modelling of the managed river reach shows that the feedforward control performed by the operator "hides" the reach time-delay in the transfer function of the closed-loop system. So, classical time-delay estimation methods are inappropriate. The proposed solution considers the delay estimation as a problem of detecting a discontinuity in impulse response. A Bayesian approach is proposed to detect the delay and to estimate the impulse response from production data. This method is applied to one-day data sets supplied over one year in order to show its efficiency. Copyright © 2003 IFAC Keywords: Estimation, time-delay, impulse response, Bayesian approach, river reach.

the control performance and, in particular, the control stabilization (Niculescu et al., 1998). In practice, the time-delays are estimated empirically (operator knowledge) or experimentally (measurement of intumescence propagation time). But in both cases, the estimates are still characterized by a wide uncertainty. Moreover, for economic reasons and safety precautions (flood risk), the implementation of experimental protocols is problematic. However, those plants are supervised and the principal variables measurements are available. The article aim is to show that an appropriate method, exploiting data sets collected in experimental conditions imposed by production modalities, allows to estimate the time-delay.

I. INTRODUCTION This al1icle presents a solution for estimating the time-delays of river flows in hydraulic valley. The latter consists of run-of-river cascaded hydroelectric plants. The problem addressed in this paper is the estimation of the time-delay between the inflow (upstream) rate and the downstream water level of a managed river reach in a discrete-time context; the purpose being, on the one hand, to estimate its nominal value and, on the other hand, to determine its evolution over one year. Few applications concern the timedelays estimation in this type of process (Thomassin er ul., 2003). Nevertheless, this estimation influences

299

The available data are measured in a closed-loop context, the system being controlled by a human operator. From the system modelling, it is shown that the feedforward control performed by the operator "hides" the reach time-
Table I. Main process variables. F'(k} F(t} Subscript i Subscript 0 z(t} and z(k} z' (k) T,

Table 2. Main notations. q 1 x(s} Ox i ht.F,' 't.F,' 'F,'

Inflow reference

*

F;' (k) Downstream water level

Outflow reference F';(k}

x;

Upstream liarrage

2. MANAGED RIVER REACH MODELLING This section gives the discrete-time I model of a downstream water level regulation of a managed river reach (Fig. 3). A one-
" ~ ,g

<:(s)

inflow reference

-

outflow reference

300

~

Z50 200

~ ~

~ ~

0

"

20

2'

X

= ~ (e-r;SFi(s) As

e-rosFo(s))

+ 8z(s),

(I)

"'7~

12865

The reach control is performed by a human operator. He adjusts the flow rate reference of each barrage to control the water level and reject the disturbances due to the inflow rate variations. Consequently, the operator action on the inflow rate reference can be approximately described by two controllers:

1216 128

ss

1285

o

5

10

15

20

2S

lime (h)

Fig:: A one-
Cf(q-I) = Kfq-n f

A managed river reach is a process composed of two barrages (upstream and downstream), a reach (the river portion between the two barrages) and a centralized control process. Each barrage includes a dam (sluices) and a hydroelectric power plant (turbines).

1

Convolution product Sub-vector [x(i} ... x(j}]T of the vector

Ft

where A is the reach water surface (m 2 ), 'Ci denotes the unknown delay between the inflow rate and the downstream water level, while the delay between the outflow rate and the downstream water level is represented by 'Co. In this application, since the water level measurement stalion is close to the downstream barrage, 'Co is fixed to zero. Note that the delay 'Ci is a slowly time-varying parameter.

000 350

Cross--<:orrelation function between z and Autocorrelation function of Ft

The river reach can be represented by a continuoustime model in which the input variables are the flow rates F; (t) and Fa (t) and the output variable is the water level z(t). Note that the inflow and outflow rates are not measured. In the absence of affluent, the river reach dynamics around an operating point can be approximately described by a model structure whose main elements are two time-
Fig. I. Managed river reach.

~

Shift operator Laplace transform of the signal x Modelling errors and measurement uncertainties of x Estimate of x Impulse response between z and F;'

These barrages have flow control loops whose response times are negligible as compared to the sampling period Ts :::::: 133s. Consequently, these loops may be modelled by a constant gain and a zero-order-hold (ZOH) (Fig. 3).

z(k}

-...

Flow rate reference of a barrage (discrete-time) Real flow rate of a barrage (continuous-time) Variable of the upstream barrage (inflow) Variable of the downstream barrage (outflow) Downstream water level (con1.- and disc.-time) Water level reference (discrete-time) Sampling period

C~(q-I) = Kz,

(2)

where Kf and K~ are the proportional gains and n f is the time-
Th~ discrete-time variables are noted x(k) and correspond to

th~ t~mpural

and

sampling with a constant sampling period T, of the variable X(I): x(k) = x(kT,).

conlinuuus-tim~

300

F;'{k)

z'(k) +

Fig. 3. Block diagram of the managed river reach control. reach control. However, the operator usually tends to reject the disturbance effects by acting on the control variable without waiting for its effect on the water level, so it is assumed that nf = O.

l10·S

12r-----~------,

Assuming that the time-delay is an integer multiple of the sampling period: 00

Tj =

nT,.

(3)

where n denotes the delay index (n EN), the discretetime model of a managed river reach can be represented by the following discrete-time equation:

30

z(k) = He.e, (q -I )z* (k) + H,.IiJq-1 )8z(k) + He.IiF;, (q- I )8Fa(k) + Hz.liF; (q-I )8F;(k)

+ H, IiF' (q-I )8F,; (k) + H. F' (q-I )F;* (k). ~.'

n

'.'

I

discrete-time k

Fig. 4. Theoretical impulse response between z and F;* and its one-day impulse response bar.

(4)

The level reference z' (k) is equal to zero due to the disturbance rejection control of the system. Note that, in this equation, only the signals F;*(k) and z(k) are known. The parameter to estimate, n, appears in the polynomial rational transfer Hz.F; (q-I ) describing the dynamical behavior between z(k) and F;(k) of the closed-loop system:

H. ,( -I) ,.F; q

=

-rYK q-I+aq-n-I vo.o f I - (I + auKz ) q- I I

3. BAYESIAN APPROACH FOR TIME-DELAY ESTIMATION

3. f Impulse response identification: an inverse problem By multiplying equation (4) by F;* (k -l) (where l ;::: 0 is the discrete-time lag), then by taking the expectation of each term, under the hypothesis that F;* is uncorrelated with modelling errors and measurement uncertainties, equation (4) becomes a discrete-time convolution equation:

(5)

where ai = KiTs/A, Ct" = KuT,/A and the coefficients a" Ct", Kt and K: are not null. It is easy to note that the number of first zero coefficients of the transfert numerator is equal to I and is different of n because of the feedforward controller. So, this prevents us from successfully estimating n in this context with classical time-delay estimation methods.

The IR estimation is thus a discrete-time identification problem and constitutes an inverse problem (ldier, 200 I). In finite dimension, by introducing an estimation of the correlation functions, equation (6) becomes:

An example of theoretical impulse response (lR) be· tween z(k) and Fj ' (k) is presented in figure 4 with n = 7. aj = 1.2.10- 4 . a,,' K J = 5.3210- 5 and aa' K: = -0. 16 (values obtained from a first system identifica· tion and a simulator). This figure also shows another way to represent the IR. the one-day impulse response bur. The value of each IR coefficient is represented by a grey level. This representation will be interesting fur what comes next. This figure shows that the pres· ence of the time-delay n introduces a discontinuity in the IR. So, a delay estimate can be obtained from the detection of the first point of this discontinuity. In the following of this article, a Bayesian approach realizing jointly the IR (h:.F:) and the time-delay (n) estimations is presented.

M-I

f:. F: (I) =

L hJ ; (j)fF,' (1- j) + b(l).

(7)

j=O

for I = O..... N - I and with M S N. The term b(l) represents a noise term due to model, truncation and correlation estimation errors and measurement uncer· tainties on data z and F;* at time I. From Eq. (7), the following linear system can be established: r = Rh+b.

with:

301

(8)

rz.Ft (0)

(h,n) = argmax[p(h,nlr,R,8)]

b(O)

h.n

(11)

= argmin [-lnp(h,nlr,R,8)]. r=

rz.F,·(M-I)

,b=

Under the hypothesis that b is a zero-mean white Gaussian noise sequence with a covariance matrix C111 and that b is independent of h, the PDF of r = Rh + b is a Gaussian distribution centered on Rh :

b(N - I)

rz,F,· (N - I) 'F,. (0)

R=

h,n

b(M-I)

'F:(-M+ I)

I

rF: (M - I) ...

[I

p(rlh,n,R,8) = (21rC11)N/2exp -2C1lllr-RhIl

2] .

(12)

rF,.(N -I)

and h = (

'F:(N -M)

h\O)

The solution, maximizing the likelihood function obtained from the previous PDF, is named the Maximum Likelihood (ML) estimation. It produces the maximum fidelity to data and is equal to the LS solution under the previous assumptions about b.

).

h(M - I) Inverse problems are often ill-posed (in sense of Hadamard) and ill-conditioned (Idier, 2001), The use of the generalized inverse allows to recover the wellposed nature of the problem, but not to resolve the illconditioning. Indeed, if the ill-conditioning of R is high, the least-square (LS) solution, which provides maximum fidelity to data, is unacceptable: a smoother solution is a priori waited; the noise is too amplified. So, some infidelity to the data must be introduced to obtain a less noisy solution and more in concordance with the a priori infonnation about h.

The prior PDF p(h,nIR,8) is equal to p(h,nI8) since we will see that it is independent of R, and from the joint probability law, it becomes:

p(h,nI8) = p(hln,8)p(nI8),

where p(hln,8) and p(nI8) represent respectively the prior information on hand n. The PDF p(hln,8) allows to express, at n fixed, a smoothness constraint on the IR, except on both points h(n) and h(n + I) defining the discontinuity. Let the vector ho = Dnh where: 2 -I

3.2 Bayesian approach: criterion design

I

h ,11

r,

R 8) = p(rlh.n,R,8)p(h,nIR,8) , p(rlR,8)'

-I

o

The Bayesian approach makes it possible to take into account the available a priori knowledge and provides a suitable framework for the criterion design. To estimate the time-delay, the parameter n, corresponding to the discrete-time of the discontinuity first point, is introduced. The a priori information on the objects n and h is expressed in the fonn of a prior probability distribution p(h, nlR, 8) where 8 is a hyperparameter vector composed of probability distributions parameters. The posterior probability density function (PDF) p(h,nlr,R,8), detennined by the Bayes' rule, allows to combine the prior infonnation with the one contained in the data in order to obtain: p(

(13)

o

0 :::

0

:: . 2-1 -I I 0 o I -I "'. -I 2 .. ". ""

0

. ' 0 '. -I -I 2

So, the vector ho corresponds to the second derivatives of h (second order finite differences approximations), except at times nand n + I where ho(n) and ho(n + I) correspond to the first derivatives obtained by respectively forward and backward first difference approximations. The smoothness constraint is imposed by assuming that the vector ho follows a zero-mean Gaussian distribution with a covariance matrix C1lD I:

(9)

The function p(rlh,n,R,8) denotes the PDF of the data conditioned on the real solutions n, hand 8, and is thus a likelihood function. The tenn: p(rlR,8)

=

J

p(rlh,n.R.8)p(h.nlR.8),

Since ho = Dnh, it comes that h follows a zeromean Gaussian distribution with a covariance matrix

(10)

ensures the nonnalization of the posterior PDF. The manipulation of PDF being cumbersome, a popUlar choice consists in using the maximum a posteriori (MAP) estimator by choosing the solution that maximizes (or minimizes the negative log-likelihood function of) the posterior PDF:

C1lD(DrDn)-I:

302

but det((D~Dn) 1) = 1, so: p(hln,e)

= (2 Jrcr ;

)M/Z exp

ho

o For n = I ... M - I, a) computation of the matrix D n , b) computation of h at n fixed, c) computation of the criterion J(h, n).

[-~IIDnhIIZ] . cr ho

(15)

ofi=arg The PDF p(nle) is detennined from our prior information about n. A Rayleigh distribution is chosen (Fig. 5): p(nle)

0] = -exp [- n

n-

cr;

2cr;

(16)

,

where crn = argmax [p(nle)] must be fixed from prior

An application example (Basse-Isere river between Grenoble and Valence, France) is carried out to show the efficiency of the proposed method. The timedelay and IR estimations are performed for each oneday production data set over one year. The choice of the data set length (24 hours) is based on several arguments. First, this period of time contains sufficient observations (N ~ 650) in order to pennit a punctual delay estimation. Next, the delay can be considered as constant over this period. The last reason is directly due to the measurement storage fonnat.

n

008·

l::

c;:

Oi"

OG4r 002 0

,

J(h,n).

4. APPLICATION RESULTS

knowledge; in our case, it is proportional to the reach length. This distribution is one possible choice; it has the advantage of taking the delay positivity into account. Moreover, the peak being sufficiently large, the delay estimate does not depend much on the choice of crn (see section 4).


min

(nEN.I~n~M-1}

1 I I I I

-~._0

er" = 7 "

"

20

"

n

'"

Fig. 5. Rayleigh distribution. From (9), (\2), (15) and (16), the negative loglikelihood function of the posterior PDF (9) becomes: I 0 I Z -lnp(h,nlr.Re) = -zlIr-RhIl-+ - Z IIDnhl1 2crb

2crho

o

n-

+ -z -In(n) + j(r,R,e), 2crn

(\7)

where j(r, R, e) = In [p(rIR e)g(e)], where e = Z [crh Uh f) U]T and g(e) = (2JrUb2)N/2(7JrU 2 11 hD )M/2 cr11' Minimize the negative log-likelihood function of the posterior PDF with respect to hand n amounts thus to minimize the following criterion:

discrete time k (a) lRs estimated by the ML estimator.

l

J(h./l) = Ilr - RhW + allD nhl1 o /l-0

+ub"

(

2

cr; -2In(n)

)

.

(18)

3.3 Optimization procedure -=-,

50 _ _

The criterion J(h.ll) is quadratic in h and, at n fixed, has an explicit solution: h(n.a)

=10

= (RTR+aD~DntlRTr.

15

2S

discrete time k (b) lRs estimated by the MAP estimator.

but J(h.n) is not convex with respect to n. However, /l being an positive integer, a simple method, not very expensive in our case and without convergence problem, consists in computing the criterion for all values in {n E N.l :S n:S M - I}. Finally, the estimation ofh and /l may be summarized by the following algorithm.

Fig. 6. Time-day representations of impulse responses. Figure 6(a) shows the IRs estimated over one year by the ML estimator (i.e. without prior information) un-

303

~

der a new representation: the time-day representation. This representation consists in placing side by side L one-day impulse response bars, where L denotes the number of one-day data sets supplied over one year. One can easily observe that the transition varies between 7 and 9. So. by "visual inspection", the timedelay varies between 6 and 8 over the year; but one has got no time-delay estimate. Better IR and timedelay estimates, in sense of our prior knowledge, have been obtained with the Bayesian approach (MAP estimator). The contribution is summarized by the comparisons of figures 6 and 7. The regularized IRs are smoother. The discontinuity is respected in spite of the smoothing and is more marked.

•

.,

.~..,

>.

u'" ,

"i'..,

.zu.''---::---::-----;:,----,!..

--e-

(a) crn

5

n_ 7 10 -

15

2Q

discrete time k

sets collected in experimental conditions imposed by production modalities. It is firstly shown that the feedforward control performed by the operator "hides" the reach time-delay in the transfer function of the closed-loop system. Consequently, classical timedelay estimation methods are inappropriate. A new time-delay estimation method based on the detection of a discontinuity in the impulse response by a Bayesian approach is proposed. This method can be improved by taking more prior knowledge into account. In practice, this approach is particulary interesting since it is based on production data and hence does not disturb the daily management of the hydroelectric power production by any experimental protocol. An application example has shown the efficiency of the proposed method.

ML l::Stimation MNJ eslimation _

25

ti (b) crn = 10

=7

Fig. 9. Time-delay estimate distribution over one year for two values of CTn .

.2:-1 _ _~'----'_-"-:_ _-:-::-_ _~'_ _~,~_--J o

f=E #,i.WlU.l.1ilJ.ll.C,,"----_.-..-l1J

ti

x to"

IL _

----------,

..,'"

30

Fig. 7. IRs estimated by ML and MAP estimators from a one-day data set. Figure 8 allows to examine the evolution of the delay index estimates over one year. More than 30 % of the delay estimates are equal to 6 and about 60 % are contained between 5 and 7 (Fig. 9(a)). Figure 9(b) shows that the delay estimates do not vary much in respect to CTn . This result confirms the appropriate choice of the prior distribution of n. Nevertheless, a detailed observation of the time-delay and IR estimates shows false detections. Indeed, for example, most of the detected discontinuities between I and 3 are decreasing ones. However, we look for increasing discontinuities. but we do not take this information in our prior knowledge. Consequently, decreasing discontinuities are detected. These false detections can be removed by taking the discontinuity direction into account in the prior knowledge. The isolated estimates are often due to non-informative data.

ACKNOWLEDGEMENTS A special thank you to Professor D. Brie (CRAN) for pertinent remarks and ad vices that he brought to us.

REFERENCES Cuno, B. and S. Theobald (1998). The relationship between control requirements, process complexity and modelling effort in the design process of river control systems. Mathematics and Computers in Simulation 46, 611 ~ 19. Idier, 1. (200 I). Approche bayesienne pour les problemes inverses. Hermes Science Publication. Kurz, H. and W. Goedecke (1981). Digital parameter adaptative control of processes with unknown dead time. Auromatica 17,245-252. Niculescu, S., E. Verriest, 1.-M. Dion and L. Dugard (1998). Stability and robust stability of timedelay systems: a guided tour. In: Stability and control of time-delay systems (L. Dugard and E. Verriest, Eds.). Springer-Verlag, London. pp. 1-71. Thomassin, M., T. Bastogne, A. Richard and A. Libaux (2003). Time-delay estimation of a managed river reach from supervisory data. In: Proceeding of the 13th 1FAC Symposium on System Identification. Rotterdam, The Netherlands.

days Fig. 8. Time-delay estimate map.

5. CONCLUSION This article presents a solution for estimating the time-delay of a managed river reach between the inflow rate and the downstream water level from data

304