VIEWS: 56 PAGES: 13 POSTED ON: 5/9/2011 Public Domain
Turkish J. Eng. Env. Sci. 29 (2005) , 171 – 183. ¨ ˙ c TUBITAK A Comparison of Recursive Least Squares Estimation and Kalman Filtering for Flow in Open Channels ¨ Omer Faruk DURDU Adnan Menderes University, Faculty of Agriculture, Aydın-TURKEY e-mail: odurdu@adu.edu.tr Received 13.10.2004 Abstract An integrated approach to the design of an automatic control system for canals using a Linear Quadratic Gaussian regulator based on recursive least squares estimation was developed. The one-dimensional partial diﬀerential equations describing open channel ﬂow (Saint-Venant) equations are linearized about an average operating condition of the canal. The concept of optimal control theory is applied to drive a feedback control algorithm for constant-level control of an irrigation canal. The performance of state observers designed using the recursive least squares technique and the Kalman ﬁltering technique is compared with the results obtained using a full-state feedback controller. An example problem with a multi-pool irrigation canal is considered for evaluating the techniques used to design an observer for the system. Considering the computational complexity and accuracy of the results obtained, the recursive least squares technique is found to be adequate for irrigation canals. In addition, the recursive least squares algorithm is simpler than the Kalman technique and provides an attractive alternative to the Kalman ﬁltering. Key words: Optimal control, Least squares estimation, Kalman ﬁlters, Open-channel ﬂow. Introduction turbances or ﬂow rate changes) are not known in ad- vance, the eﬀect of the random disturbances on the system variables must be measured and used in the With increasing demand for food and competing use feedback loop to control the system. The variation in within the water sector, the pressure is on irriga- the depths of ﬂow is used in the closed loop (feedback tion professionals to manage water eﬃciently. In re- loop). During the last few decades, several control al- sponse to this, strategic decisions and interventions gorithms have been developed to derive the relation- need to be made on a continuous basis. These deci- ship between the deviations in the system variables sions should cover the full spectrum of the irrigation (ﬂow depth and ﬂow rate) and the change in gate water supply system, from diversion and distribu- opening (gate-control algorithm). Linear quadratic tion to on-farm application down to the crop root optimal control theory has been applied for deriving zone. Supply-oriented systems have not been able closed-loop control algorithms for real-time control to provide the needed ﬂexibility in terms of water of open irrigation canals (Reddy, 1999; Reddy et al., quantity and timing to achieve improved crop yields 1999; Seatzu et al., 2000; Seatzu, 2002; Durdu, 2003; and water use eﬃciency. This calls for a demand de- Durdu, 2004). However, when lumped parameter livery of water to the farmers in the command area models are used to derive control algorithms for ir- of an irrigation project. Demand delivery oﬀers the rigation canals, the number of state variables (ﬂow maximum ﬂexibility and convenience to the water depths and ﬂow rates) that must be used in the feed- user (Reddy, 1990). In a demand delivery schedule, back loop is large. Consequently, it is costly to mea- under constant level control, since the demands (dis- 171 DURDU sure ﬂow depths and ﬂow rates at several points in in which ql = lateral discharge rate, m3 /s; Cd = out- a multi-pool irrigation canal. Therefore, to mini- let discharge coeﬃcient; bl = width of outlet struc- mize the cost of implementing feedback control algo- ture, m; wl = height of gate opening of outlet struc- rithms, the number of measurements per pool must ture, m; Z= water surface elevation in the supply be kept to an absolute minimum. Since 1 or 2 ﬂow canal, m; Zl = water surface elevation in the lat- depths per pool are normally measured in practice, it eral canal, m; and Es = sill elevation of the head is preferable (and possible) to estimate values for the regulator, m. In the regulation of the main canal, state variables that are not measured. This is done decisions regarding the opening of gates in response by using an observer. Reddy (1995) demonstrated to random changes in water withdrawal rates into Kalman ﬁltering for the estimation of state variables lateral canals are required to maintain the ﬂow rate in the control system. The objective of this paper into laterals close to the desired value. This is ac- is to develop an LQG regulator using the recursive complished either by maintaining the depth of ﬂow least squares method and to demonstrate the perfor- in the immediate vicinity of the turnout structures mance of this regulator in comparison to full-state in the supply canal constant or by maintaining the feedback and a Kalman ﬁlter based LQG regulator. volume of water in the canal pools at the target value (Reddy, 1999). When the latter option is used, the Basic Equations outlets are often ﬁtted with discharge rate regulators. The water levels or the volumes of water stored in The Saint-Venant equations, presented below, are the canal pools are regulated using a series of spa- used to model ﬂow in a canal tially distributed gates (control elements). Hence, irrigation canals are modeled as distributed control systems. Therefore, in the solution of Eqs. (1) and ∂A ∂Q + − ql = 0 (1) (2), additional boundary conditions are speciﬁed at ∂t ∂x the control structures in terms of the ﬂow continuity and the gate discharge equations, which are given by (Reddy, 1990) 1 ∂Q 1 ∂ Q2 ∂x + = gS0 − gSf − g (2) A ∂t A ∂x A ∂y Qi−1,N = Qgi = Qi,1 (continuity) (6) in which Q = discharge, m3 /s; A = the cross- sectional area, m2 ; ql = lateral ﬂow, m2 /s; y = wa- ter depth, m; t = time, s; x = longitudinal direction along the channel, m; g = the acceleration due to Qgi = Cdi bi ui (2g(Zi−1,N − Zi,1 ))1/2 (gate discharge) gravity, m2 /s; So = canal bottom slope (m/m); and (7) Sf = frictional slope, m/m, and is deﬁned as in which, Qi−1,N = ﬂow rate through downstream gate (or node N ) of pool i − 1, m3 /s; Qgi = ﬂow Sf = Q|Q|/K 2 (3) rate through upstream gate of pool i, m3 /s; Qi,1 = ﬂow rate through upstream gate (or node 1) of pool in which K = hydraulic conveyance of canal = i, m3 /s; Cdi = discharge coeﬃcient of gate i; bi = AR 2/3 /n; R = hydraulic radius, m; and n = Man- width of gate i, m; ui = opening of gate i, m; Zi−1,N ning friction coeﬃcient, s/m1/3 . Lateral canals in = water surface elevation at node N of pool i − 1, the main canal are usually scattered throughout the m; Zi,1 = water surface elevation at node 1 of pool length of the supply canal. The mathematical rep- i, m; and i = pool index (i = 0 refers to the up- resentation of ﬂow through these structures is given stream constant level reservoir). The Saint-Venant as follows (Reddy, 1999): open-channel equations, Eqs. (1) and (2), are lin- earized about and average operating conditions of the canal to apply the linear control theory concepts ql = Cd bl wl (2g(Z − Zl ))1/2 for submerged ﬂow to the problem. After applying a ﬁnite-diﬀerence ap- (4) proximation and the Taylor series expansions to Eqs. (1) and (2), a set of linear, ordinary diﬀerential equa- tions is obtained for the canal with control gates and ql = Cd bl wl (2g(Z − Es ))1/2 for free ﬂow (5) turnouts (Durdu, 2004) 172 DURDU where Φ = (AL )−1 ∗ AR , Γ = (AL)−1 ∗ B, and Ψ = (AL )−1 ∗ C, δx(k) = nxn state vector, δy(k) = + A11 δQ+ + A12 δzj + A13 δQ+ + A14 δz + 1 = j j+1 rxn vector of output (measured variables), H = rxn A11 δQj + A12 δzj + A13 δQj+1 + A14 δzj+1 + C1 output matrix, and r = number of outputs. The ele- (8) ments of the matrices Φ, Γ, and Ψ depend upon the canal parameters, the sampling interval, and the as- sumed average operating condition of the canal. In + + A21 δQ+ j + A22 δzj + A23 δQ+ + A24 δzj+1 = j+1 Eq. (11), the vector of state variables is deﬁned as follows A21 δQj + A22 δzj + A23 δQj+1 + A24 δzj+1 + C2 (9) δx = (δQi,1 , δZi,2 , δQi,2 , . . . δZi,N−1 , δQi,N−1, δQi,N ) where δQ + and δz + = discharge and water-level in- j j (13) crements from time level n + 1 at node j; δQ j and δzj = discharge and water-level increments from time Full-State Feedback Regulator level n at node j; and A11 , A21 , . . . . A12 , A22 are the coeﬃcients of the continuity and momentum equa- The full-state feedback (Linear Quadratic Regulator tions, respectively, computed with known values at (LQR)) control problem is an optimization problem time level n. Similar equations are derived for chan- in which the cost function, J, to be minimized is nel segments that contain a gate structure, a weir given as follows or some other type of hydraulic structure. From the matrix form of Eqs. (8) and (9), the state of system K∞ equation at any sampling interval k can be written J= [δxT (k)Qxnxn δx(k) + δuT (k)Rmxm δu(k)] in a compact form as follows (Durdu, 2003) i=1 (14) AL δx(k + 1) = AR δx(k) + Bδu(k) + Cδq(k) (10) subject to the constraint that: where A = nxn system feedback matrix, B= nxm control distribution matrix, C = pxn disturbance −δx(k + 1) + Φδx(k) + Γδu(k) = 0k = 0, . . . , K∞ matrix, δx(k) = nxn state vector, δu(k) = mxn con- (15) trol vector, ∆δq = variation in demands (or distur- bances) at the turnouts, m2 /s, n = number of depen- where K∞ = number of sampling intervals consid- dent (state) variables in the system, m = number of ered to derive the steady state controller; Qx nxn = controls (gates) in the canal, p = number of outlets state cost weighting matrix; and Rmxm = control in the canal, and k = time increment, s. The ele- cost weighting matrix. The matrices Qx and R are ments of the matrices A, B, and C depend upon the symmetric, and to satisfy the non-negative deﬁnite initial condition. The dimensions of the control dis- condition, they are usually selected to be diagonal tribution matrix, B, depend on the number of state with all diagonal elements positive or zero. The ﬁrst variables and the number of gates in the canal. The term in Eq. (3) represents the penalty on the devia- dimensions of the disturbance matrix, C, depend on tion of the state variables from the average operating the number of disturbances acting on the canal sys- (or target) condition, where the second term repre- tem and the number of dependent state variables sents the cost of control. Equations (3) and (15) (Malaterre, 1997). Equation (10) can be written in a constitute a constrained-minimization problem that state-variable form along with the output equations can be solved using the method of Lagrange multi- as follows pliers (Reddy, 1999). This produces a set of coupled diﬀerence equations that must be solved recursively backwards in time. In the optimal steady-state case, δx(k + 1) = Φδx(k) + Γδu(k) + Ψ δq(k) (11) the solution for change in gate opening, δu(k), is of the same form as δy(k) = Hδx(k) (12) δu(k) = −Kδx(k) (16) 173 DURDU where K is given by condition of the system may be either above or below the equilibrium condition, depending upon the sign and magnitude of the disturbances. If the system de- K = [R + ΓT SΓ]−1 ΓT SΦ (17) viates signiﬁcantly from the equilibrium condition, the discharge rates into the laterals will be diﬀer- S is a solution of the discrete algebraic Riccati ent (either more or less) from the desired values. In equation (DARE) canal operations, however, the main objective is to keep these deviations to a minimum so that a nearly ΦT SΦ − ΦT SΓ[R + ΓT SΓ]−1 ΓT SΦ + Qx = S constant rate of discharge is maintained through the (18) turnouts (Reddy et al., 1992). where R = RT > 0 and Qx = Qx T = H T H ≥ 0. Gaussian Regulator with Kalman Filtering The control law deﬁned by Eq. (16) brings an ini- tially disturbed system to an equilibrium condition The Linear Quadratic Gaussian (LQG) theory pro- in the absence of any external disturbances acting on vides an integrated knowledge base for the devel- the system. In the presence of these external distur- opment of a ﬂexible controller (Figure 1). Since it bances, the system cannot be returned to the equilib- is expensive to measure all the state variables (ﬂow rium condition using Eq. (16). An integral control, rates and ﬂow depths) in a canal system, the num- in which the cumulative (or integrated) deviation of ber of measurements per pool must be kept to an a selected output variable is used in the feedback absolute minimum. Usually the ﬂow depths at the control loop, is required to return the system to the upstream and downstream ends of each pool are mea- equilibrium condition in the presence of external dis- sured. The relationship between the state variables turbances (Kwakernaak and Sivan, 1972). Integral and the measured (or output) variables is control is achieved by appending additional variables of the following form to the system dynamic equation δy(k) = Hδx(k) + η(k) (21) −δxn (k + 1) = Dδx(k) + Γδxn (k) (19) where η(k) is measurement error inputs. For a in which δx n =integral state variables; and D = in- steady-state Kalman ﬁlter, the observer gain matrix, tegral feedback matrix. This produces a new control L, is calculated as follows law to the form (Reddy et al., 1992) L = P H T [HP H T + RC]−1 (22) δu(k) = −Kδx(k) − Kn δxn (k) (20) where P is the covariance of estimation uncertainty The ﬁrst term in Eq. (20) accounts for initial disturbances, whereas the second term accounts for external disturbances. Equation (20) predicts the ΦT P Φ − ΦP H T [RC +H T P H]−1HP ΦT +Qesti = P desired gate openings as a function of the measured (23) deviations in the values of the state variables (Reddy et al. 1992). In this paper, the water surface eleva- where RC = RC T > 0 is a tolerance value for the RC tion and ﬂow rate were considered the state variables. covariance matrix, which is an identity matrix and Given initial conditions [δx(0)], δu, and δq, Eq. (3) Qesti = QT ≥ 0 is a diagonal matrix. The dis- esti can be solved for variations in ﬂow depth and ﬂow turbances δq(k) and η(k), in Eqs. (11) and (21), are rate as a function of time. If the system is really at assumed to be zero mean Gaussian white noise se- equilibrium [i.e. δx (0) = 0 at time t = 0] and there quences with symmetric positive deﬁnite covariance is no change in the lateral withdrawal rates (distur- matrices Qesti and RC, respectively. Furthermore, bances), the system would continue to be at equilib- sequences δq(k) and η(k) are assumed to be statis- rium forever and there is then no need for any control tically independent. The system dynamic equation action. Conversely, in the presence of disturbances is used to predict the state and estimation error co- (known or random), the system would deviate from variance as follows: time update equations (Tewari, the equilibrium condition (Reddy, 1990). The actual 2002) 174 DURDU Measurement noise η(k) Γq Process noise q(k) δu(k) q-1 δx(k) δy(k) + Γu H Φ Controlled system δu(k) + Γ q-1 H δy(k) ^ + - Φ δx(k) ^ -K L(k) Regulator + Figure 1. Linear Quadratic Gaussian (LQG) controller with a state estimator. If the initial conditions and the inputs (control inputs and the disturbances) are known without er- P− (k + 1) = ΦP (k)ΦT + ΨQesti ΨT (24) ror, the system dynamic equation, Eq. (2), can be used to estimate the state variables that are not mea- sured. Since parts of the disturbances are random and usually are not measured, the canal parameters δˆ − (k + 1) = Φδˆ(k) + Γδu(k) x x (25) are not known very accurately, and the estimated values of the state variables would diverge from the x in which δˆ (k) = estimated values of the state vari- actual values. This divergence can be minimized ables. As soon as measured values for the out- by utilizing the diﬀerence between measured output put variables δy(k) are available, the time-update and the estimated output (error signal), and by con- values are corrected using the measurement update stantly correcting the system model with the error equations as follows: measurement update equations signal (Reddy, 1995). Therefore, the modiﬁed state (Tewari, 2002) equations are given as L(k + 1) = P−(k + 1)H T [HP−(k + 1)H T + RC]−1 x x x δˆ (k + 1) = Φδˆ(k) + Γδu(k) + L[δy(k) + Hδˆ (k)] (26) (29) P (k + 1) = [I − L(k + 1)H]P−(k + 1) (27) Gaussian Regulator with Recursive Least Squares This estimation technique is useful in identifying δˆ (k + 1) = δˆ − (k + 1) + L(k + 1)[δy(k + 1)− x x time varying parameters and has been considerably discussed in the literature. The optimal estimation Hδˆ − (k + 1)] x criterion for recursive weighted least squares is writ- (28) ten in the following form 175 DURDU on the previous optimal estimate is obtained by (El- k Hawary, 1989) J= β(k, j)εT (j)Λ−1 (j)ε(j) (30) j=1 δˆ− (k) = Φ(k − 1)δx+ (k − 1) x (38) where ε(j) = prediction error in vector form, and β(k,j) = weighting function. The weighting function whereas, in the least squares method, the transition is assumed to satisfy the following relationship matrix is assumed to be unity and thus β(k, j) = λ(k)β(k − 1, j), 1 ≤ j ≤ k − 1, β(k, k) = 1 δˆ − (k) = δx+ (k) = δx(k) x (39) (31) In addition, the error covariance matrix is ob- where λ is the weighting or forgetting factor. These tained by relations imply that P− (k) = Φ(k − 1)P+ (k − 1)ΦT (k − 1) + ΓQesti(k − 1)ΓT k (40) β(k, j) = Π λ(i) (32) ε=j+1 For the least squares method the matrix Qesti is The recursive estimation algorithm can be writ- zero and the equivalent of Eq. (16) is given by ten in a number of diﬀerent equivalent forms. The following form of the recursive least square algorithm has many computational advantages (El-Hawary, P−(k) = P+ (k − 1) (41) 1989) 2) in the corrector stage of the Kalman ﬁlter, an updated state estimate equation is obtained δˆ = δˆ (k − 1) + K(k)r(k) x x (33) in which r(k) is the innovations sequence deﬁned by δˆ + (k) = δˆ−(k) + K(k)r(k) x x (42) in which r(k) is the innovation sequence and is given r(k) = y(k) − H(k)δˆ(k − 1) x (34) by The gain matrix K(k) is deﬁned by r(k) = y(k) − Hδx− (k) (43) K(k) = P (k − 1)H T (k)Φ−1 (k) (35) Equations (12) and (13) of the least squares method are the same as Eqs. (21) and (22). In addi- where Φ(k) is deﬁned by tion, an update of the covariance matrix is obtained Φ(k) = λ(k)Λ(k) + H(k)P (k − 1)H T (k) (36) P+ (k) = (I − K(k)H(k))P− (k) (44) The equivalent of the state error covariance ma- This deﬁnition is similar to Eq. (16) except for trix is given by the division by λ(k). The Kalman gain matrix K is given by 1 P (k) = [1 − K(k)H(k)]P (k − 1) (37) λ(k) K(k) = P−(k)H T (k)Φ(k)−1 (45) The equations for the least squares method bear many similarities to the Kalman ﬁltering equations. where The diﬀerences between Kalman ﬁltering and the least squares method are: 1) in the predictor stage of the Kalman ﬁlter, a prediction of the state based Φ(k) = R(k) + H(k)P− (k)H T (k) (46) 176 DURDU Equation (24) is similar to Eq. (14), but Eq. Ljung and Soderstrom (1983) take α to be unity. (25) diﬀers from Eq. (15), since R(k) in Eq. (25) The convergence of the ﬁlter is inﬂuenced by the is replaced by λ(k)Λ(k) in Eq. (15). In the least choice of the weighting function. Once the equations squares method, one does not need the value of the of the optimal state feedback and the recursive least state-noise covariance matrix Qesti , and the mea- squares method are obtained, and measured values surement error covariance matrix RC(k) is replaced for one or more state variables for each pool are avail- by the choice of the weighting function. The prob- able, the dynamics of the linear system can be sim- lem of deﬁning the noise statistics manifested by the ulated for any arbitrarily selected values of external Qesti and RC matrices of Kalman ﬁltering can be disturbances. In this study, a multi-pool irrigation seen from references. The approaches for deﬁning canal was considered. The algorithm predicts the Qesti and RC matrices involve a considerable com- ﬂow rate, Q(x, t), and the depth of ﬂow, y(x, t), given putational burden, which may be avoided by suitably the initial boundary conditions. The optimal state selecting the weighting function to suit the practical feedback and the least squares equations were added application (El-Hawary, 1989). The choice of weight- as subroutines to this algorithm. Given the initial ing function β(k,j) controls the way in which each ﬂow rate and the target depth at the downstream end measurement is incorporated relative to other mea- of the each pool, the algorithm computed the back- surements. The choice should clearly be such that water curve. Later on, the downstream ﬂow require- measurements that are relevant to current system ment and the withdrawal rate into the lateral were properties are included. If one chooses to assign less provided as a boundary condition. The model pre- weight to older measurements such as in the case of dicted the depths and ﬂow rates at the nodal points variable system parameters, then a popular choice is for the next time increment. The computed depths given by at the upstream and downstream ends of each pool were used with the least squares under constraints to estimate the ﬂow depths and ﬂow rates at some λ(i) = λ (47) selected intermediate nodal points. These estimated values were then used in the optimal state feedback Therefore, Eq. (11) can be written as subroutine to compute the change in the upstream gate opening in order to bring the depth at the down- k k−j stream end of the pool close to the target depth. β(k, j) = Π λ=λ (48) ε=j+1 When the estimated values of the state variables are used in the feedback loop, the controller equation, This is referred to as an Exponential Weighting Eq. (16) becomes (Reddy, 1999) into the Past (EWP), with λ being a forgetting fac- tor that shapes the estimator’s memory. Λ is chosen to be slightly less than 1. A second possible choice δu(k) = −Kδˆ(k) − Kl δxl (k) x (52) is such that the forgetting factor is time-varying and in this case (El-Hawary, 1989) Equation (52) computes the desired change in gate opening as a function of the estimated (instead of measured) deviations in the state variables. Based β(k, j) = (αk )k−j (49) upon this gate opening, the new ﬂow rate into the pool at the upstream end was calculated and used where αk is deﬁned by the ﬁrst order discrete ﬁlter as the boundary condition at the upstream end of the each pool. This process was repeated during the αk = λ0 αk−1 + (1 − λ0 )α (50) entire simulation period. Typically, α0 = 0.95, λ0 = 0.99, and 0 < α < 1. Results and Analysis The forgetting factor starts at α0 and reaches a steady state value of α (El-Hawary, 1989). For a To explore and compare the performance of the reasonably large k Kalman estimator and recursive least squares method, an LQG regulation problem for a discrete- time multi-pool irrigation canal was simulated. The λ(k) = αk−1 = λ0 [λ(k − 1) − α] + α (51) data used in the simulation study are demonstrated 177 DURDU in Table 1. These data were ﬁrst used to calculate pensive, the control algorithm ﬁrst estimated state the steady state values, which in turn were used to variables using a Kalman ﬁlter. Next, the recursive compute the initial gate openings and the elements least squares method was employed to estimate the of the Φ, Γ, and H matrices using a sampling in- values for the state variables in the algorithm. The terval of 30 s. After computing steady state values, Kalman ﬁlter for the system used the control input the control algorithm formulates an LQG controller δu(k), generated by the optimal state feedback, mea- with a Kalman ﬁlter and recursive least squares, re- sured water depths δy(k) for each pool, the distur- spectively. As a ﬁrst part of the LQG controller, bances noise δq(k), and measurement noise, η(k). In a full-state feedback controller (assuming all state the design of the Kalman ﬁlter, in lieu of actual ﬁeld variables are available) was designed to regulate the data on withdrawal rates from the turnouts, the ran- 6-pool canal system using a constant-level control dom disturbances were assumed to have some pre- approach. The system response was simulated using speciﬁed levels of variance. The variances of the dis- the controller in the feedback loop. In the deriva- turbances, Qesti , were w1 = 12 , w2 = 1.32, w3 = 0.72 , tion of the feedback gain matrix K, the control cost w4 = 1.42 , w5 = 1.32 , and w6 = 1.22. The actual time weighting matrix, R, of dimension 6, was set equal series of the demands was not used in the design of to 100, whereas the state cost weighting matrix, Qx, the ﬁlter; only the variance of the time-series was was set equal to an identity matrix of dimension 85. required in the design of the ﬁlter. Usually the sen- The matrix dimension 85 comes from the system di- sors used to measure ﬂow depths in an open-channel mension. Since the irrigation canal is divided into are reasonably accurate to a fraction of a centimeter; 49 nodes and each node has a set of 2 equations, therefore, the variance of the measurement error is the dimension of the system should be equal to 98. usually very small. A value of 0.0005 was used for However, the system has 7 gates and 6 turnouts; the variance of the measurement matrix (RC ), and therefore, the system matrix dimensions numbered this was an identity matrix. Using the given initial 85. The cost weighting matrix and the control cost values, the system response was simulated for 250 matrix must be symmetric and positive deﬁnite (i.e. time increments or 7500 s. After designing the LQG all eigenvalues of R and Qx must be positive real controller using a Kalman ﬁlter, the algorithm de- numbers). A priori, we do not know what values of signed a recursive least squares based on a deﬁned Qx and R will produce the desired eﬀect. In the weighting factor (β). The analysis was started by absence of a well-deﬁned procedure for selecting the evaluating the system stability. All the eigenvalues elements of these matrices, these values are selected of the feedback matrix were positive and had val- based upon trial and error. At ﬁrst, both Qx and ues less than 1. The system was also found to be R as identity matrices were selected. By doing so, both controllable and observable. In the derivation it was speciﬁed that all state variables and control of the control matrix elements, Γ, it was assumed inputs were equally important in the objective func- that both the upstream and downstream gates of tion, i.e. it was equally important to bring all the each reach could be manipulated to control the sys- deviations in the state variables (water surface ele- tem dynamics. The last pool’s downstream-end gate vations and ﬂow rate) and the deviations in the con- position was frozen at the original steady state value, trol inputs to zero while minimizing their overshoots. and only the upstream gates of the given canal were Note that the existence of a unique, positive deﬁnite controlled to maintain the system at the equilibrium solution to the algebraic Riccati equation (Eq. 18) condition. The eﬀect of variations in the opening is guaranteed if Qx and R are positive semi-deﬁnite of the downstream gate must be taken into account and positive deﬁnite, respectively, and the system is through real-time feedback of the actual depths im- controllable. To test whether the system was con- mediately upstream and downstream of the down- trollable, the system controllability matrix was cal- stream gate (node N ). Figure 3 demonstrates the culated and it was found that the system was control- variations in ﬂow depths for each pool and for all 3 lable. In the derivation of the feedback gain matrix, techniques. The variations in ﬂow depths for recur- R was set equal to 100, whereas Qx was set equal to sive least squares were compared to the variations in an identity matrix of dimension 85 (the dimensions ﬂow depths computed using optimal state feedback of the system). After deﬁning Qx and R matrices, as well as a steady-state Kalman ﬁlter. Since pool 1 the optimal feedback gain matrix, K, was calculated. is the ﬁrst pool of the irrigation canal, with an in- Since measurement of all the state variables was ex- crease in ﬂow rate into the lateral (turnout) or down- 178 DURDU stream demand, the depth of ﬂow at the downstream downstream pools and meeting the demand at the end of pool 1 decreased rapidly and approached a downstream end. Figure 4 demonstrates the incre- maximum deviation of -0.157 m for least squares, - mental gate openings for each design technique (opti- 0.15 m for optimal state feedback and -0.13 m for mal state feedback, Kalman ﬁlter and least squares) the Kalman ﬁlter at approximately 2000 s from the and for each gate in the canal. The deviation in the beginning of the disturbance period. By the end of gate openings for recursive least squares was com- the simulation, the system returned very close to the pared with the deviation in gate opening computed original equilibrium condition for all 3 techniques. In using optimal state feedback as well as the steady- pool 2, in the ﬁrst 1700 s of the simulation, the ﬂow state Kalman ﬁlter. At the beginning, gate 1 had depth decreased dramatically and reached a maxi- sharp peaks for all 3 design techniques. Since opti- mum deviation of -0.144 m for optimal state feed- mal state feedback has the best stability properties, back, -0.141 m for least squares and -0.118 m for the the state feedback curves will be the target loop. At Kalman ﬁlter. The variations in ﬂow depth in pool 3 gate 1, incremental gate opening for recursive least reached -0.162 m for least squares and -0.16 m for op- squares was closer to optimal state feedback (target timal state feedback and -0.127 for the Kalman ﬁlter loop function) than were those for the Kalman ﬁlter. after around 1500 s of the simulations period. Pool After 6000 s, gate 1 reached an equilibrium position 4 and pool 5 have less ﬂuctuation in comparison to for all 3 techniques. At gates 2, 3, 4, 5 and 6, in- the other pools. Pool 6 has the highest variations in cremental gate opening values for the Kalman ﬁlter ﬂow depth and the depth of ﬂow at the downstream were far away from the optimal state feedback values end decreased rapidly and approached a maximum in comparison to the least squares values. At the end deviation of -0.32 m for optimal state feedback, -0.28 of the simulation, the variations in the gate openings m for least squares and -0.25 for the Kalman ﬁlter at (for all gates) approached a constant value, indicat- around 2000 s of the simulation period. The rapid ing that a new equilibrium condition was established. decreases in the downstream depth of ﬂow in each It is obvious from the ﬁgures that the variations in pool resulted in an attendant sudden increase in the ﬂow depth and in gate openings for recursive least gate opening at the upstream end of each reach to squares are closer to the target loop function (opti- release more water into the pool. However, because mal state feedback) than are those for the Kalman of the wave travel time, the depth of ﬂow at the estimator. In other words, the recursive least squares downstream end did not start to rise until around ﬁlter has better stability properties than the Kalman 1700 s. In all the pools considered, the maximum ﬁlter in the control of irrigation canals. Moreover, in deviation in depth of ﬂow occurred at the ﬁrst and the computation of estimation, one does not need to last pools of the canal for all 3 design techniques. ﬁnd the disturbance covariance matrix, Qesti , and To meet the downstream target depth, the last pool measurement covariance matrix, RC. Therefore, a had the highest ﬂuctuations. The ﬂuctuations in the computational burden may be avoided by selecting ﬁrst pool were due to releasing more water into the the appropriate weighting function. Gate 1 Gate 2 Gate 3 Supply Pool 1 Gate 4 Pool 2 Gate 5 Reservoir Pool 3 Gate 6 Pool 4 Gate 7 Pool 5 downstream Pool 6 Q=5 m3/s 52,500 m Figure 2. Schematic of multi-pool irrigation canal. 179 DURDU Figure 3. Comparison of ﬂow depth variations for a full-state feedback regulator using Kalman ﬁltering and a regulator using recursive least squares. 180 DURDU Figure 4. Comparison of incremental gate openings for a full-state feedback regulator using Kalman ﬁltering and a regulator using recursive least squares. 181 DURDU Table 1. Data used in the simulation study. parameter pool 1 pool 2 pool 3 pool 4 pool 5 pool 6 length, m 9000 9000 9000 9000 9000 9000 width, m 5 5 5 5 5 5 bottom slope 0.0002 0.0002 0.0002 0.0002 0.0002 0.0002 side slope 1.5 1.5 1.5 1.5 1.5 1.5 initial lateral ﬂow depth, m3 /s 2.5 2.5 2.5 2.5 2.5 2.5 initial downstream depth, m 4 3.2 2.86 2.5 2.12 1.81 parameter gate 1 gate 2 gate 3 gate 4 gate 5 gate 6 width, m 5 5 5 5 5 5 discharge coeﬃcient 0.8 0.8 0.8 0.8 0.8 0.8 initial opening, m 1.13 1.37 1.16 0.97 0.85 0.63 disturbances, m3 /s 2.5 2.5 2.5 2.5 2.5 2.5 Conclusions best robustness and stability properties, it was cho- sen as a target loop function. The results obtained In this paper, recursive least squares estimation has from simulations indicate that the least squares al- been employed in the control of multi-pool irrigation gorithm provides good stability and performance in canals to estimate the state variables (ﬂow depth and the control of irrigation canals. The algorithm is sim- ﬂow rate) at intermediate nodes based on the mea- pler than Kalman ﬁltering in terms of the knowledge sured variables. The performance of the recursive of covariance matrices required and provides an at- least squares was compared to the performance of tractive alternative to Kalman ﬁltering. Overall, the the optimal state feedback and the Kalman ﬁlter in performance of the recursive least squares technique terms of variations in the depths of ﬂow and the up- for constant-level control was better than the perfor- stream gate opening. Since the full-state feedback mance of the Kalman ﬁlter. (assuming all state variables are measured) has the References Durdu, O.F., Robust Control of Irrigation Canals, Malaterre, P.O., “Multivariable Predictive Control PhD Dissertation, Colorado State University, Civil of Irrigation Canals. Design and Evaluation on a 2- Engineering Department, Hydraulic Division, Fort pool Model”, International Workshop on Regulation Collins, CO, USA, 2003. of Irrigation Canals, Morocco, 230-238, 1997. Durdu, O.F., “Estimation of State Variables for Reddy, J.M., “Local Optimal Control of Irrigation Controlled Irrigation Canals via a Singular Value Canals”, Journal of Irrigation and Drainage Engi- Based Kalman Filter”, Fresenius Environmental neering, 116, 616-631, 1990. Bulletin, 13(11), 2004. Reddy, J.M., Dia, A. and Oussou, A., “Design El-Hawary, F., “A Comparison of Recursive of Control Algorithm for Operation of Irrigation Weighted Least Squares Estimation and Kalman Canals”, Journal of Irrigation and Drainage Engi- Filtering for Source Dynamic Motion Evaluation”, neering, 118, 852-867, 1992. OCEANS ’89. Proceedings, 4, 1082-1086, Septem- ber 18-21, 1989. Reddy, J.M., “Kalman Filtering in the Control of Irrigation Canals”, Int. J. Appl. Math. Modeling, Kwakernaak, H. and Sivan, R., Linear Optimal Con- 19(4), 201–209, 1995. trol Systems. Wiley Publishers, New York, 1972. Reddy, J.M., “Simulation of Feedback Controlled Ljung, L. and Soderstrom, T., “Theory and Practice Irrigation Canals”, Proceedings USCID Workshop, of Recursive Identiﬁcation”, MIT Press, Cambridge, Modernization of Irrigation Water Delivery Sys- MA, 1983. tems, 605-617, 1999. 182 DURDU Reddy, J.M. and Jacquot R.G., “Stochastic Opti- Systems Science, 31(6), 759-770, June, 2000. mal and Suboptimal Control of Irrigation Canals”, Seatzu, C. and Usai G., “A Decentralized Volume J. Water Resources Pln. and Mgt., 125, 369-378, Variations Observer for Open-Channels”, Applied 1999. Mathematical Modelling, 26(10), 35-61, October Seatzu, C., “Design of Decentralized Constant- 2002. Volume Controllers for Open-Channels by Solving Tewari, A., Modern Control Design with Matlab and a Least Squares Problem”, International Journal of Simulink. Wiley Publishers, New York, 2002. 183