Recursive state and parameter estimation with applications in water resources

Recursive state and parameter estimation with applications in wate, r resources W. Schilling lnstitutfiir Wasserwirtschaft, Universit~ttHannover, Callinstrasse32, D-3000Hannover, West Germany

J. Martens lngenieurb•ro Dilger, lm Biittelwoog2, D6783Dahn, West Germany (Received November 1985; revised February1986)

Hydrologic models, as well as measurements of hydrologic processes, are corrupted by noise. The Kalman filter is a convenient tool to estimate the true but unknown state of a hydrologic system. It is, however, difficult to specify the necessary error covariances. A procedure is proposed to estimate the error covariances recursively in a combined state and parameter filter. Applications of the procedure yield meaningful results for two hydrologic data series of very different character. A major benefit of the proposed algorithm seems to be its robustness against instability. Key words: Kalman filter, error covariance, multivariate regression model, rainstorm forecasting; dissolved oxygen modelling Numerical modelling of hydrologic processes has become a common technique in water resources engineering. After the development of satisfactory simulation models for rainfall/run-off processes, hydraulic' routing in channels and pipes, water quality in receiving waters, etc., more recent research is focusing on the control of these processes. The task is to manipulate the process (i.e. alleviate flows, control pollution discharges, etc.) as it occurs, according to a pre-specified performance criterion. Conceptually, these real time control problems might be formulated in the framework of optimal control theory or optimization theory. The real-time environment imposes severe constraints on process models. It is necessary that: • • • • •

states have to be readily observable noisy measurements be processed parameters be adapted to the current situation small computers be used modelling results be available in real time

The latter requirement of computational efficiency, especially, often prohibits the use of physically-based models. The numerical schemes to integrate the partial differential equations of mass and energy balance are mostly too time-consuming to be solved in real time. Simplified conceptual or statistical models have to be

applied, instead. These include multivariate regression techniques where a number of independent variables such as rainfall, water quality parameters, etc. can be auto- and cross-correlated to simulate their temporal and spatial behaviours. Measurement and model errors can be taken into account by the Kalman filter method.

Modelling concept Here, a multivariate regression model in state space format is discussed: xr+, = Frr+ WT

(1)

In this process model, a first order state transition matrix F is assumed to be unknown. The state Xr is directly observable through noisy measurements Zr: Zr = Xr + Vr

(2)

Both modelling, wr, and measurement errors, Vr, are random with constant, unknown means, w and v, and covariances, Q and R. In the parameter model, the elementsfof the transition matrix F are arranged row-wise as a vector Xr*. Their time behaviour is assumed to follow a random walk pattern around their unknown means:' XT÷t* = Xr* + Wr*

(3)

0307-904X/86/06433-05/$03.00 © 1986Butterworth & Co. (Publishers) Ltd

Appl. Math. Modelling, 1986, Vol. 10, December

433

Estimation in water resources: J. Martens and W. Schilling

The current, XT, and previous, xr_l,states of the system can similarly be linked by the parameters Xr* to

T

w = 1/T Z

x r = Hr*xr* + Vr*

(4)

I-x>, r 1 "r*= [ xr-,:..Xr_,r XT-I T

(5)

The means, w*, and covariances, Q*, of the parameter model errors are assumed to be zero and constant, respectively. Since the parameter estimation procedure is to be unbiased, the means £* = E ( x r - Hr*xr*) are zero. The covariances R* describe the uncertainty of the current state of the systems. Therefore, R*= PT/r, which is defined in the following section.

State and parameter estimation Due to the simplified modelling concept, the simulated states deviate from the true states. Furthermore, since measurements ZT are corrupted by the noise VT, the true states cannot be observed. Hence, they are unknown and have to be estimated. For the linear case as formulated in equations (1) and (2), the Kalman filter yields the optimum, i.e. best possible, state estimate.-' The Kalman filter algorithm requires a priori knowledge of the error statistics, i.e. means _vand w, and covariances R and Q. Usually these quantities are assumed. For falsely specified error statistics, however, the optimality properties of the Kalman filter are lost. Since the true states are always unknown in real world systems, i.e. cannot be simulated or measured without error, the correct specification of the error statistics is difficult. To overcome this practical problem, algorithms have been proposed to estimate states and error statistics simultaneously.3-5 These algorithms have been applied by the authors to problems which will be discussed and found to yield unstable results, i.e. unbounded states and error statistics. Since the equations for the error covariances Q and R both include the error covarianees of the estimated, PT/T, and forecasted, PT/T-~, states, respectively, the error covariances Q and R are not independent. Their independence, however, is a prerequisite to apply the Kalman filter. It is believed that this dependence is the reason for the unfavourable results obtained with the procedures referred to above. A new algorithm, directly derived from the definitions of the error terms was developed 6 and successfully applied to the hydrologic modelling problems discussed shortly. The error statistics are defined by: v = E(zr-

Xr/r)

(6)

R = E ( ( v r - v ) ( v r - v) T)

(7)

w = E(xr/r-i - Fxr-l/r-l)

(8)

O = e ( ( w r - w_)(Wr-- w) r)

(9)

These quantities can be estimated by (6): T

v_ 11r~. (z,-x,/,) =

(lO)

1=1 T

R = I/T Z

(v t - v_)(vt - v) T

(12)

Q = I/T Z

(w., - w ) ( w t - w) T

(13)

1= I

where:

I

(x,,/H - Fxr-t/,-l)

1=1 T

the state observation equation(2):

(ii)

I=I

434 Appl. Math. Modelling, 1986, Vol. 10, December

In equations (14)-(27), the computational steps for simultaneous state and parameter estimation with two Kalman filters are given for one time interval. Assuming initial states and parameters, the new state of the system is forecast (equation 14). After the time span of one interval passes ( T > T + 1), the forecasts are compared to real-time measurements z r, corrected with the difference between measurement and forecast (state estimator, equations 16-18), and the error statistics are recursively updated (equations 19-22). Then the parameters in the transition matrix Fare rearranged in vector format Xr*, forecasted and corrected (parameter estimator, equations 23-27). Finally, the corrected parameters are used for the next state forecast (equation 14), etc.: XTtT_ I = FXT/T-I- W

(14)

Pr/T-I = FPT/T FT + Q

(15)

K r = Pr/T_I(PT/T + R) -1

(16)

Xr/r = Xr/T-t + K r ( z r - XT/T-I -- VT)

(17)

PT/T = ( I - KT)PzTr-i

(18)

v = ((T-

1)/T)V_T_ I + 1 / T ( z r - - x T / r )

(19)

g = ((T-

1)/T)RT_~ + 1 / T ( Z F - - X r / T - I - - V T )

x ( Z r - Xr/r-i - V_r)r E = ((T-

(20)

1)/T)WT_I + (1/T)(XT/T_ t -- Fxr_t/T_l)

(21) Q = ( ( T - 1)/T)QT-I "t'- ( 1 / T ) ( ( X T / T _

l -- F x r _ l / r

i-Wr)

X (XT/T_ I -- FXT_I/T_ I -- WT) T

(22)

Xr/r-l* = X r - l / r - l * + w*

(23)

PT/r-l* = Pr/r* + Q*

(24)

Kr* = PT/T-I*Hr*r (HT*PT/rHT*T + PT/T) -t

(25)

Pr/r* = ( I - Kr*HF*)Pr/T-j*

(26)

Xr/r* = Xr/r-l* + K r * ( X r / r - H r * x r / r - i * )

(27)

S t r e a m flow dissolved oxygen modelling The combined state and parameter estimator was applied to model dissolved oxygen in the Leine River south of Hannover in West Germany. Here, reliable oxygen forecasts could be a valuable tool to ceatrol the oxygen concentration in the river, e.g. by artificial aeration or reduced pollution discharges. The twodimensional system state comprises the oxygen concentration in the river at an upstream location, 1, and at a site, 2, located 40 km downstream. The approximate flow time between the two sites is 24h. Modelling errors are caused, for example, by pollution discharges within the reach of the river, time-varying photosynthesis, and varying flow time. As can be seen in Figure 1, the oxygen concentrations are slightly scattered around the diurnal variations caused by photosynthesis. Forecasts of oxygen concentrations 24 h ahead are generated every hour. Since the initial values of the model error covariances

Estimation in water resources: J. Martens and W. Schilling

10.8 I0 _'Z"

=o s

== o u

iI^t

---.--

Forecast Measurement

- -

Updated

concentration

9.2 8.4 7,6

6.8

"

6

i

,%-/

,'

0

"

[

)

5.2 4.4 I I0

3.6 0

I 2.0

I 50

I 40

I 50

I 60

I 70

Time

1.4 I

I 80

I 90

I I00

I I10

I 120

step (h)

b

1.2 .o

I

0.8 E

0.6 0.4

-..................../:.!.,~

c

E

.

0.2 0

":i"~'"[~"~'~"i,-'~ ' -- I

I0

0

20

--]

Iv ~ . . ~ ' ' l ~ ' " ~ r ' ~ - , : , , i

5£)-40

50

60

70

....... i"~'~~---~l

80

~/(I

I

I

.....

"f2, I

I

,.....,,

I

. ..,

I

90 I00 I10 120:'180 190 200 210 220 23,0 ?-40 Time step (h}

0.7?_

i.i

c 0.64

IE L.

o.g ,~0.56 t 0.48 -

0.8 ~,

0.4C

0.7

E

~.

u c

0.6

8

0.24

0.5

o,.

0.16.

0.4

8 0.52 I q2~ ~A

0.08

~

.

~

0.:5

,__ _ / - - - ~

I

I

I

I

I

I0

20

3,0

40

50

I

i

I

60 70 80 Time step (h}

I

I

I

I

90

I00

I10

120

Figure 1 Hourly forecasts (24 h ahead): (a) measurements and estimates of oxygen concentrations in Leine River at downstream site; (b) variations of regression parameters; (c) error statisticsq and r

Appl. Math. Modelling, 1986,Vol. 10, December

435

Estimation in water resources: J. Martens and W. Schilling 262

,I, I00 'c

90

E

80

E

70

0

60

-;

50

-

4.0

E

.s E

144

I

I|

'iit

!!

"-.'---~

!I

= = Measurement at site 5 - - o - - - o - Forecast at site 6 .-o o- Measurement at site 6

Forecast at site 5

t

30 20 I0

,'x.

o;

4

8

12

16

20

.

.

24

.

.

,

,,_

26 3,2 36 Time s t e p ( x l 5 m i n )

, 40

,/.;", 44

"91, - - , ' W 48

52

56

60

I 56

I 60

b 1.61 1.4: C e~

1.2

j-

I

.a

0.8 o= ¢3n

0.6' 0.4 0.2 0 -O.I 0

I

I

t

4

8

12

I 16

I 20

I 24

I I I 28 32 36 Time step (x 15 rain)

I 40

I 44

I 48

I 52

C

0.24

240

E

o,2a

220

~J

0.2

200

E 0.18

18o

6

g o.16

,6o

oJ

'c

% E

v0

0.14 u ¢~ 0 . 1 2

140 120

~

~"

,oo

o

o.i o.oe

80

0.06

60

o.o4

40

0.02

L

20 4

8

12

16

I 20

I 24

I

I

I

I

28 32 3,6 40 Time s t e p ( x 15 rain)

I

I

I

I

I

44

48

52

56

60

0

Figure 2 (a} one step forecasts (15 min) and measurements of rainstorm intensities at two out of 12 gauge sites; (b) variation of regression parameter; (c) error statistics q and r

436

Appl. Math. Modelling,

1986,Vo1010, December

Estimation in water resources: J. Martens and W. Schilling Table 1 Variables of state and parameter estimators: in text, italics indicate vectors or matrices; bold italics indicate estimates Variable

Definition

K

Filter {Kalman) gain for state update n× n Error covariances of forecasted states nx n Error covariances of updated states n × n Modelling error n× 1 Measurement error nx 1 Covariance matrix o f modelling errors n× n Element of O Covariance m a t r i x of measurement nx n errors Element o f R State vector n× 1 State transition matrix n× n Element o f F Observation matrix of parameter model n × n2 Filter gain for parameter update n2x n Error covariances of forecasted n2x n2 parameters Error covariances of updated parameters n 2 x n 2 State vector of parameter model (Fin n 2 x 1 vector format)

PWT-1 PT/T

w v Q q R r x F f H* K*

Pr/r- i* PT/T~ XT*

Dimension

in Q were chosen to be large (Figure 1(c)) the updated oxygen concentrations Xr/r are close to the measurements (Figure l(g)). After about 25 time steps, the error covariances in the Kalman filter have converged and the model produces meaningful forecasts. At later stages of the process, the differences between forecasted and updated concentrations are still larger than the differences between measurements and updated concentrations. Therefore, the model error covariances (equation (22)) are always greater than the measurement error covariances (equation (20) and Figure 1(c)). The error statistics converge rapidly with only minor changes at later time steps. The four parameters j~ in the state transition matrix converge more slowly. More than 200 time steps are required to reach steady-state parameters (Figure l(b)). The choice of their initial values does not influence their final values and at no time do the results become unstable.

original time invariant model derived in reference 7 constant transition matrix F, was modified to allow parameter adaption in situations of unusual storm behaviour, e.g. speed and direction of propagation. To achieve this, the denominators in equations (19)-(22) of the error statistics were replaced by a constant, e.g. T = 10. Figure 2 shows rainfall forecasts and corresponding measurements one time step ahead (15 min) for two of the 12 gauge sites and for a period of 15h. With the occurrence of the first significant change of the storm intensity, i.e. after time step 14 at site 6 (Figure 2(a)), the corresponding covariances of the modelling and measurement errors q66 and r~, respectively, rise substantially (Figure 2(c)). In periods without rainfall, i.e. T = 34-45, the parameters in the state transition matrix F hardly change. This is demonstrated by the behaviour of f56 which is one of the 144 transition parameters. It indicates how the rainfall intensity at site 6, at time T + 1, depends on the intensity at site 5 at time T (Figure 2(b)). During the dry periods, the errors w and v are relatively small. Therefore, their covariances converge towards smaller limits (Figure 2(c)). As soon as it starts to rain again the parameter filter adjusts the transition parameters. This is shown for f56 in Figure2(b). Despite the high variability of the storm intensities, the updating procedure does not become unstable at any time.

Conclusions An algorithm is proposed to estimate the unknown error statistics in the Kalman filter. The algorithm is derived from the definitions of the modelling and measurement errors, respectively. The algorithm was successfully tested in two hydrologic case studies using a combined state and parameter filter.

Acknowledgements This contribution is based on research funded by the German Research Association (DFG) under the grants Le 229/19 and Si 242/5. The computer time granted by the Regional Computing Centre of Lower Saxony (RRZN) is appreciated.

References Short term rainstorm forecasting A few municipalities recently started to control combined storm and sanitary sewer flows in real time to use available in-pipe storage, alleviate treatment plant inflows and minimize sewage diversions into receiving waters. To achieve optimum operation, the storm inflows to the sewer system have to he forecasted. Based on telemetered rain gauge data, a short-term forecasting model in the format of equation (1) is proposed. The system state xr comprises 12 rainfall depths in adjacent subcatchments during time interval T, measured with 12 rain gauges. Other than in the dissolved oxygen case, the measurements vary irregularly and considerably, as do the model parameters and error statistics. Instantaneous changes of rainfall intensities by several hundred percent are not uncommon. The

1 Mayne, D. Q. 'Optimal non-stationary estimation of the parameters of a linear system with Gaussian inputs', J. Electr. Control, 1963, 14, January, 101-112 2 Gelb, A. (ed.) 'Applied optimal estimation' MIT Press, Cambridge, MA, 1974, pp. 374 3 Sage, A. P. and Husa, G. W. 'Adaptive filteringwith unknown prior statistics', Proc.Joint Aurora. ControlConf., 1969, 760-769 4 Myers, K. A. and Tapley, B. D. 'Adaptive sequential estimation with unknown noise statistics', 1EEE Tram. Aurora. Control, 1976, August, 520-523 5 Todini, E., Szrllrsi-Nagy, A. and Wood, E. F. 'Adaptive state/ parameter estimation algorithmsfor real-timehydrologicforecasting: a case study', IIASA/WMO Workshop Recent Developments Real-Time Forecasting/Control Water Resource Sys., October 1976, Laxenburg,Austria 6 Martens, J. 'Water quality forecastingusing the Kalman filtertechnique' in Planning and control of urban storm drainage, (Balmer, P. et el. eds.), GOteborg,Sweden, 1984, 821-830 7 Schilling,W. 'Real-time estimation and forecasting of spatiallydistributed areal rainfall"War.Sci. Techn., 1984, 12 (8/9), 327-348

Appl. Math. Modelling,

1986,Vol.

10, December

437