# A Comparison of Recursive Least Squares Estimation and Kalman by ert634

VIEWS: 56 PAGES: 13

• pg 1
```									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

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-
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

```
To top