- Email: [email protected]
Identification of Properties of Open Water Channels for Controller Design
16th IFAC Symposium on System Identification The International Federation of Automatic Control Brussels, Belgium. July 11-13, 2012
Identification of Properties of Open Water Channels for Controller Design P.J. van Overloop, X. Bombois
Abstract— This article describes a way to use standard System Identification techniques to determine all pool properties of open water channels for controller design. These properties are the delay time, the storage area and the frequency and magnitude of the first resonance peak. The identification procedure is tested on an actual open water channel at the Central Arizona Irrigation and Drainage District, Eloy, AZ. The identification experiment results in excellent estimations of all properties. Especially the value for the storage area improves considerably compared to previous research on this specific open water channel.
I
I. INTRODUCTION
rrigation systems in the Western part of the United States consist of multiple open water channels in series, delivering water to remote areas. To maximize the yield of these agricultural areas, it is important that the delivery of water to these sites is accurate and flexible. This can be seen as the main function of an irrigation canal. The operational management of the water system infrastructure hence, plays a vital role in the performance of an irrigation canal. Automated irrigation systems can meet the goals of being accurate and flexible, if the controllers that operate the structures are correctly tuned. A control algorithm that is simple, has a good performance over the full range of flow conditions and keeps the water levels at a predefined target level is the Proportional Integral (PI) feedback controller. Multiple successful implementations of this type of controller have been reported in the Western part of the US, Europe and Australia (Burt & Piao (2002), Malaterre & Baume (1998), Ooi & Weyer (2008)). Often, the water level measurement is located at the downstream end of the channel and this measurement signal is communicated to the upstream structure (see Figure 1). This setup is referred to as distant downstream control. It has the advantage that the offtakes, which are generally located at the downstream end of the open water channel, have an upstream water level that is (close to being) constant. The dynamic behavior of open water channels can be described with the well-known Saint-Venant equations (Chow (1959)). These are partial differential equations that consist of a mass balance and momentum equations. Some channels have structures, such as culverts, which can be described by another type of momentum equation. The momentum Peter-Jules van Overloop is with the Department of Water Management, Delft University of Technology, Delft, The Netherlands (e-mail: [email protected]). Xavier Bombois is with the Delft Center for Systems and Control, Delft University of Technology, Delft, The Netherlands (e-mail: [email protected]).
978-3-902823-06-9/12/$20.00 © 2012 IFAC
equations are non-linear. As control theory is mainly based on linear systems and ordinary differential equations, simplification of the dynamic behavior is required. In Schuurmans (1997), a division into two classes of open water channels is made. One class is defined as channels that are long and shallow. The dynamic behavior is mainly determined by the long delay times of the water that has to travel from the upstream side to the downstream end of the channel. In Schuurmans et al. (1999) and later in Litrico & Fromion (2006), research has been conducted to find stable and well-performing controllers for these channels without an extensive trial-and error tuning procedure. Their research has resulted in controller tuning rules for this class of open water channels based on the delay time τd and the storage area As of the back water part. The other class of channels can be characterized by being short and deep and is dominated by reflecting waves. After a change in the boundaries of the system, for example the opening of a gate or a pump that is turned off, a wave travels up and down the channel a number of times before the water surface settles again. Also in Litrico & Fromion (2004a), Overloop (2006) and Ooi & Weyer (2009), these so-called resonance waves are described and analyzed. The main reason for the badly damped behavior in short and deep channels is the low bed friction force along the open water channel. Figure 1 shows a representation of the properties of the two classes of open water channels. For the resonance class open water channel, only the resonance mode with the lowest frequency is shown (first peak in frequency spectrum). It is important to consider exactly this frequency, as the measurement location for distant downstream control is in counter-phase (phase shift of -180º) with the control actions. According to the Nyquist stability criterion, this control system will, in closed loop, be unstable, when the magnitude of this resonance peak, in open loop, is larger than 1. Schuurmans (1997) proposed to filter the resonance wave with a first-order low-pass filter and developed a tuning rule for the filter. This rule requires input of the frequency ωp and the magnitude of the first peak Mp of the resonance wave in the frequency spectrum. At present, stable and well-performing PI controllers for all classes of open water channels can be directly determined when the properties delay time τd, storage area As, frequency ωp and peak magnitude Mp of the lowest resonance mode can be estimated accurately at low and high flow conditions (Schuurmans (1997)). Other implementations of feedback controllers (Weyer 2003) require the first resonance frequency for selection of the cut-off frequency of the applied high-order low-pass filter.
1019
10.3182/20120711-3-BE-2027.00414
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Qin(k)
h(k) Qin(k-k d) As
Delay time τd
Qin(k)
h(k)
Storage area
Qout (k)
Qout (k)
Figure 2. Open water channel variables and parameters
Figure 1. Delay-time-dominated open water channel (left) and resonancedominated open water channel (right) including distant downstream control loop
In Overloop et al. (2010), a identification procedure for identifying the resonance properties of an open water channel in the Central Arizona Irrigation and Drainage District in Eloy, Arizona is presented. The identified model in this research gave satisfactory results for the resonance, but was inaccurate for low frequencies. The storage area, that needs to be identified from this low-frequency part, showed unlikely values twice as high as the pool’s top width multiplied by the length. In this article, the measurements of the Overloop et al. (2010)’ experiment are used to demonstrate an improved procedure that gives accurate values for both resonance frequency and the low-frequency range. It is a modification to the procedure described in Overloop et al. (2010). Therefore large parts of that article are repeated in this article. First, the general dynamic behavior of open water channels is described. Next, system identification methods for estimating the properties of an open water channel are discussed. After this, the applied system identification experiment on the actual open water channel is described. This section includes the results. The last sections are the conclusions of this article. II. OPEN WATER CHANNEL DYNAMICS The flows and water levels in an open channel can be described by the Saint-Venant equations (Chow (1959)). These nonlinear hyperbolic partial differential equations consist of a mass balance and a momentum balance. The momentum balance is a summation of the descriptions for the inertia (1), advection (2), gravitational force (3) and friction force (4): A f Q (Eq. 1) + = qlat x t (Eq. 2) g Q Q Q Q2 h t 1
+
+ =0 + g Af x Af x C 2 R f Af 2
3
4
To use the formulas in a numerical model of an open water channel, the partial differential equations are discretized in time (Δt) and space (Δx). In case these discretized formulas are simulated, the model results in time series solutions of water levels and flows at discrete locations along the channel. Also, the time series are discrete solutions in time. This numerical model is non-linear and of high order, which makes it hard to identify and use for controller design. When we consider the entire open water channel at once, on average, the water flows through an open water channel with an average velocity V (m/s) that can be calculated as the average flow Q divided by the average wetted area Af. A wave front on top of this velocity travels with a higher velocity that is characterized by the celerity c (m/s) (Chow (1959)). The celerity is commonly estimated from the average depth d (m) according to Equation 3. (Eq. 3) c g d The most compact description of an open water channel is the Integrator Delay model proposed by (Schuurmans (1997)), where the input is the upstream water flow and the output is the water level downstream in the pool. This Laplace model given in Equation 4, describes only the low-frequency part and neglects resonances. H s e d (Eq. 4) Qs As s where H(s) is the Laplace transform of the water level downstream in the pool and Q(s) is the Laplace transform of the flow of the upstream gate. For a given frequency ω the main property for controller design, storage area As, can be calculated from the magnitude M(ω) of this transfer function according to Equation 5.
1 (Eq. 5) M The equation is also valid in the frequency region of higher order models where the (low-frequency) integrator action is present. As
3
where Q represents the flow (m /s), t the time (s), x the distance (m), Af the wetted area of the flow (m2), qlat the lateral inflow per unit length (m2/s), g the gravitational acceleration (≈9.81 m/s2), h the water level (mMeanSeaLevel), C the Chézy friction coefficient (m1/2/s) and Rf the hydraulic radius (m). Rf is calculated by Af over Pf, where Pf represents the wetted perimeter (m). Figure 2 gives a schematic representation of a typical open channel with its variables and parameter.
III. IDENTIFICATION METHODS APPLIED TO OPEN WATER CHANNELS
Deriving low-order models and open water channel properties from analytical or numerical models can be risky, as there can be mismatches between actual channel and model. To identify the actual properties, field experiments can be executed. Below, various identification methods are described.
1020
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
A step test is commonly used to find the delay time τd and the storage area As as modeled in the Integrator Delay model. In this type of tests, starting from an initial steady state in the channel, a positive step change in flow is applied at the upstream structure, while the downstream flow is kept constant. The measured water level trajectory typically shows a delay time before the water level start to rise with a more or less constant slope that is determined by the size of the storage area. Figure 3 presents the result of a step test on an open water channel.
Figure 3. Step response of model of open water channel
Another field experiment to identify the delay time and the storage area is applied in the Auto Tune Variation proposed by Litrico et al. (2007). This method applies a switching pattern of upward and downwards steps. The switches are made when the water level reaches the allowed minimum or maximum of the band around the target level. The resonance properties, especially the resonance peak magnitude Mp, are harder to identify. In Figure 3, a resonance wave can be distinguished, super-positioned on the sloped trajectory. The resonance frequency can be roughly estimated according to Equation 6, assuming the celerity c is considerably higher than the flow velocity V.
p
c L
(Eq. 6)
where L (m) is the length of the open water channel. An experiment to identify the resonance peak magnitude can be executed by applying a Proportional feedback controller of which the gain is increased until instability occurs. The resonance peak is the reciprocal of the gain that gives the instability. Another experiment is the application of an input chirp signal (Miltenburg (2008)). A recent advance in estimating the resonance is described in Clemmens et al. (2012), where the ATV procedure is extended to include resonance waves. The previously described experiments aim at finding one or a subset of the properties that are required to tune a controller for resonance-dominated channels. Potentially, the engineering field of System Identification (Ljung (1999)), in which models with low order are estimated from input and output signals, can produce the full subset of required properties. Silvis et al. (1998) apply SI-techniques to identify
the Integrator Delay model (Schuurmans (1997)) from experimental data. In Eurén & Weyer (2007) the same model of an integrator and a delay in series is used. The difference between the two approaches is that the first article applies flows at the control structures as input variables, while the second article uses the positions of the undershot and overshot control gates. Rivas Perez et al. (2007) identify a second order model in series with a delay time. In Litrico & Fromion (2004b) the first step upwards is linked to the high frequency behavior of resonance-dominated water channels that consists of the first resonance and its higher harmonics. In that article, a low order model is proposed that contains a zero in the transfer function in order to model this highfrequency resonance behavior. Weyer (2001) applies a delay time in series with a third order model. The reasoning behind this is that this model captures the storage area behavior (first order) in parallel with the first resonance wave (second order oscillating mode). In Overloop et al. (2010), which was based on Miltenburg (2008), various model structures were tested on actual measurement data in order to find the structure that best captures the frequency response of the system around the first resonance mode. The best fit and most realistic resonance properties were found with a Box-Jenkins model of order five. The low-frequency part however resulted in an unlikely estimate of the storage area As of twice the surface area. IV. SYSTEM IDENTIFICATION APPLIED TO AN ACTUAL OPEN WATER CHANNEL
Experiments are conducted at the irrigation district office of the Central Arizona Irrigation and Drainage District in Arizona (www.caidd.com). The water level and flow measurements are sent to the central computer at a sample rate of once per 10 seconds (= Ts). The flow input signal is sent to the gates at a rate of once every 20 seconds. Communication is realized by means of radio-signals. As the open water channels are part of an operational delivery system, the experiments took place at night. Multiple experiments on various open water channels of the irrigation district are performed. In this article, the experiment on open water channel NB12 is presented, which is supposed to be resonance-dominated. Other tests can be found in Miltenburg (2008). The dimensions of the channel are presented in Table 1. Table 1. Dimensions of open water channel NB12 Length
Bed
Bottom
Side
Average
Base flow
(m)
slope (-)
width (m)
slope (-)
depth (m)
(m3/s)
2.0e-4
0.61
1.5
0.80
1021
786
0.2
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
These dimensions result in a situation where the entire open water channel is under backwater (Schuurmans (1997)), which indicates that the entire pool functions as a reservoir with a storage area of approximately the surface area of top width multiplied by length. Figure 4 presents the Bode diagram of a high-order linearized Saint-Venant model (Δx=10) of the resonance-dominated channel NB12 with a customary Mannings-friction coefficient of 0.02 s/m1/3 ( ≈ Chézy 43 m1/2/s ). In this Bode diagram, it can be observed that the frequency of the first resonance peak is ωest ≈ 0.0085 rad/s and that the magnitude at that frequency has an magnitude of approximately Mp,est ≈ 0.24. Note that the estimate for the frequency without using a Saint-Venant model, according to Equation 6 is 0.011 rad/s. From the lowfrequency part, the storage area can be calculated using Equation 5 as being As ≈ 2222 m2. This is close to the surface area of the top width multiplied by the length, which is 2366 m2. In the sequel, prediction error identification will be used to refine the model estimates of the first resonance frequency ωp and of its magnitude Mp, as well as the delay time τd and storage area As. The estimates of these quantities deduced from Figure 4 will be used in the design of the identification experiment. Note that the representation of the dynamic behavior as presented in Figure 4 comes from a model and can be inaccurate due to approximations and modeling assumption. Experiment setup In order to excite the system, a Random Binary Signal (RBS) is super-positioned on the upstream flow. The amplitude of this RBS is here approximately equal to 0.075 m3/s. In general, the amplitude should be of at least 25% of the base flow in order to guarantee a sufficient signal-to-noise ratio. The upstream flow perturbed in this manner is applied at the upstream gate, while keeping the downstream outflow constant. The initial condition in the channel has to be steady and the water level has to be close to target level. In order to contain enough information, the duration of the experiment is chosen to be equal to approximately 25 times the fundamental period 2.π/ωest of the estimated first resonance frequency. The duration of the experiment was therefore of approximately 5 hours and 20 minutes which corresponds to 1937 data points with the available sampling time of 10 seconds.
Figure 4. Bode diagram of high-order linearized transfer function from upstream flow to downstream water level
The to-be-identified system contains an integrator (see also low-frequency behavior in Figure 4). It is well known that identification of a system with an integrator is troublesome, as it is on the edge of being unstable. Consequently, the output data have to be differentiated using a tamed differentiator D. A tamed differentiator only differentiates the signal up to a certain frequency ωd in order to avoid any amplification of high frequency noise. In order not to affect the resonance peaks, ωd needs to be selected considerably lower than the resonance frequency. In Overloop et al. (2010) the frequency ωd was chosen equal to 0.002 rad/s, approximately five times smaller than the expected resonance frequency ωest. Recent research shows that this caused the storage area to be estimated inaccurately. A better result can be found by using ωd as 0.0001 rad/s, approximately 100 times smaller than the estimated resonance frequency. With the original sampling time Ts = 10 seconds, the frequency response of the system could be described up to 0.31 rad/s (the Nyquist frequency π/Ts). In this frequency range, the frequency response of the system is expected to present a large number of resonance peaks (see Figure 4). In order to avoid a model of a too large order and to avoid the classical problems arising when identifying dynamics in distant frequency ranges, the input-output data are decimated by a factor 10 (after the application of an anti-aliasing filter). The new sampling time is therefore equal to Ts* = 100 seconds. This choice for the decimation factor is the results of trial and error.
1022
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
^
function with sampling time Ts* = 100 seconds. The model W obtained in this way is (Ts* = 100s): ^
W z 2
0.0225 z 9 0.1365 z 8 0.06435 z 7 0.1743z 6 0.003077 z 5 0.03225 z 4 0.04087 z 3 0.001483z 2 0.001691z z 9 0.5689 z 8 0.6226 z 7 0.003737 z 6 0.02749 z 5 0.09991z 4 0.01092 z 3 0.1917 z 2 0.04559 z1 0.3161
(Eq. 9) The Bode diagram of the ninth order model is represented in Figure 6, where the integrator behavior in low frequencies and the first three resonance peaks can be observed. The model identified in Overloop et al. (2010), which is a fifth order model, is presented as well (dashed line). It can be seen that the ninth order has a higher low-frequency part that is in much better correspondence with the high-order numerical model as presented in Figure 4. Figure 5. Model residuals for the identified Box-Jenkins model
Model identification using prediction error identification Prediction error identification (Ljung (1999)) will now be used in order to deduce a model between the input signal and the (differentiated) output signal. A Box-Jenkins structure is used for this purpose. This model structure models the process transfer function as well as the noise perturbing the output. Even though a description of the noise is not required here, it was decided to model the noise as well, since it can be proven that this is more optimal statistically (see Ljung (1999)). In order to find a model order able to describe the underlying process, the order of the process model and of the noise model were increased in a stepwise manner until the identified models were deemed validated i.e. until the socalled residual tests (see Ljung (1999)) indicated that the prediction error could be both considered as uncorrelated with the input and as representing the characteristics of a white noise. At the same time the models needed to show a frequency response with three resonating peaks in the range of [0, 0.031].This procedure has led to the identification of the following process model (Ts* = 100s) of order eight for the process model, of order five for the noise model and containing two delay steps 1: ^
G z 2
Figure 6. Bode diagram of 9th (solid) and 5th (dashed) order transfer function of open water channel
Validation of the model Now, the validation data set to verify the accuracy of model ^
W is used. The input signal in this data set is simulated with ^
model W and this simulated output ysim is compared with the actual measured output signal yval in the validation data set (after de-trending). This delivers Figure 7. The fit defined as:
0.0225z 8 0.1142 z 7 0.1774 z 6 0.001342 z 5 0.004405z 4 0.02789 z 3 0.01325z 2 0.001708z1 z 8 0.4311z 7 0.1916 z 6 0.1953z 5 0.1678z 4 0.06792 z 3 0.07884 z 2 0.2706 z1 0.3161
(Eq. 7) and the following noise model:
FIT (%) 1
z 0.08885 z 0.6504 z 0.1828 z 0.3607 z 0.17 H 8 z 0.5921z 4 0.5962 z 3 0.508 z 2 0.3378 z1 0.6165 5
^
4
3
2
1
(Eq. 8)
where z-1 is the discrete backwards step of Ts* seconds. Figure 5 presents the residual tests corresponding to the ^
^
model { G , H } and shows that this model is indeed validated. ^
The model G is a model relating the applied flow input to the differentiated water level. In order to determine a model relating effectively the flow and the water level, the identified
1 N
N
y
val
k 1
1 N
(k ) ysim (k )
N
y k 1
val
2
(Eq. 10)
(k )2
is here equal to 89%. Resulting properties Finally, in this results section, the properties are derived from the identified ninth order model in comparison with the fifth order model as found in Overloop et al. (2010). Table 2 presents the determined properties.
^
model G above must be multiplied by a the inverse of D*. D* is the tamed differentiator D re-sampled to produce a transfer
1 The number of delays in the model has been determined by trial and error and by considering the dimensions of the considered open water channel as given in Table 1.
1023
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Table 2. Properties of open water channel found from identification experiment Model
5th order 9th order
Estimated delay time estimated τd (s) 200
Estimated storage area As (m2) 5217
Estimated resonance frequency ωp (rad/s) 0.0088
Estimated resonance peak magnitude Mp 0.223
200
2668
0.0088
0.249
The storage area is calculated using Equation 5 at a frequency of 1e-3 rad/s. The ninth order model is now well capable of identifying this area close to the surface area of top width multiplied by length of the open water channel (2366 m2).
Figure 7. Validation of estimated model (solid) and measurements (dashed) of the open water channel (top-graph is input signal flow around base flow at upstream side of the channel, while bottom-graph is output signal water level around average depth at downstream side).
V. CONCLUSIONS A field experiment to identify all properties of an open water channel required for controller design is successfully tested on an actual open water channel in Arizona. In this case, a ninth order model was required to capture all relevant dynamics at low frequencies and at the frequency of the first resonance peak. The main improvement compared to the procedure as presented in earlier research (Overloop et al. (2010)) was the capability of identifying the low-frequency part of the transfer function between upstream flow and downstream water level. Because of this, the storage area can now be estimated accurately.
Eurén, K., Weyer, E. (2007), ‘System identification of open water channels with undershot and overshot gates’, Control Engineering Practice, Vol 15, p813-824. Litrico, X., Fromion, V. (2004a), ‘Frequency modeling of open-channel flow’, Journal of Hydraulic Engineering, Vol130, Issue 8, p806-815 ASCE. Litrico, X., Fromion, V. (2004b), ‘Analytical approximation of openchannel flow for controller design’, Journal of Applied Mathematical Modelling, Vol 28, Issue 7, p677-695. Litrico, X., Fromion, V. (2006), ‘Tuning of robust distant downstream PI controllers for an irrigation canal pool. I: Theory’, Journal of Irrigation and Drainage Engineering, Vol 132, Issue 4, p359-368, ASCE. Litrico, X., Malaterre, P.O., Baume, J.P., Vion, P.Y., Ribot-Bruno, J. (2007), ‘Automatic tuning of PI controllers for an irrigation canal pool’, Journal of Irrigation and Drainage Engineering, Vol 133, Issue 1, p27-37, ASCE. Ljung, L. (1999), ‘System Identification, theory for the user’, New Jersey, PTR Prentice Hall. Malaterre, P.O., Baume, B.P. (1998), ‘Modeling and regulation of irrigation canals: existing applications and ongoing researches’. Proceedings of IEEE Conference on System, Man and Cybernetics, San Diego, Vol 4, p3850-3855. Miltenburg, I.J. (2008), ‘Determination of canal pool properties by experimental modeling’, Master-thesis Delft University of Technology, The Netherlands. Ooi, S.K., Weyer, E. (2008), ‘Control design for an irrigation channel from physical data’, Control Engineering Practice, Vol 16, p11321150. Ooi, S.K., Weyer, E. (2009), ‘Detection of oscillations in control loops in irrigation channels’, Control Engineering Practice, doi:10.1016/j.conengprac.2009.04.010 Overloop, P.J. van (2006), ‘Model Predictive Control on open water systems’, Ph.D.-dissertation Delft University of Technology, The Netherlands. Overloop, P.J. van, I.J. Miltenburg, X. Bombois, A.J. Clemmens, R.J. Strand, N.C. van de Giesen, R.W. Hut (2010), ‘Identification of resonance waves in open water channels’, Control Engineering Practice. Rivas Perez, R., Feliu Battle, V., Sanchez Rodriguez, L. (2007), ‘Robust system identification of an irrigation main canal’, Advances in Water Resources, Vol 30, Issue 8, p1785-1796. Schuurmans, J. (1997), ‘Control of water levels in open channels’, Ph.D.-dissertation Delft University of Technology, The Netherlands. Schuurmans, J., Clemmens, A.J., Dijkstra, S.J., Hof, A. & Brouwer, R. (1999), 'Modeling of Irrigation and Drainage Canals for Controller Design', Journal of Irrigation and Drainage Engineering, Vol 125, Issue 6, p338-6344, ASCE. Silvis, L.G., Hof, A., Hof, P.M.J. & Clemmens, A.J. (1998), ‘System identification on open channels’, Proceedings 3rd International Conference on Hydroinformatics. Weyer, E. (2001), ‘System identification of an open water channel’, Control Engineering Practice, Vol 9, Issue 12, p1289-1299. Weyer, E. (2003), ‘LQ control of an irrigation channel’, Proceedings of the 42nd IEEE Conference on Decision and Control, Maui, Hawaii USA, p750-755.
REFERENCES Burt, C.M., Piao, X. (2002), ’Advances in PLC-Based Canal Automation’, Proceedings USCID-Conference on Benchmarking Irrigation System Performance Using Water Measurement and Water Balances, San Luis Obispo, California. Chow, V.T. (1959), ‘Open-Channel hydraulics’, McGraw-Hill Book Co. Inc, New York. Clemmens, A.J., Litrico, X., Overloop, P.J. van, Strand, R.J., (2012), ‘Estimating Canal Pool Resonance with Auto Tune Variation’, Journal of Irrigation and Drainage Engineering.
1024
Identification of Properties of Open Water Channels for Controller Design P.J. van Overloop, X. Bombois
Abstract— This article describes a way to use standard System Identification techniques to determine all pool properties of open water channels for controller design. These properties are the delay time, the storage area and the frequency and magnitude of the first resonance peak. The identification procedure is tested on an actual open water channel at the Central Arizona Irrigation and Drainage District, Eloy, AZ. The identification experiment results in excellent estimations of all properties. Especially the value for the storage area improves considerably compared to previous research on this specific open water channel.
I
I. INTRODUCTION
rrigation systems in the Western part of the United States consist of multiple open water channels in series, delivering water to remote areas. To maximize the yield of these agricultural areas, it is important that the delivery of water to these sites is accurate and flexible. This can be seen as the main function of an irrigation canal. The operational management of the water system infrastructure hence, plays a vital role in the performance of an irrigation canal. Automated irrigation systems can meet the goals of being accurate and flexible, if the controllers that operate the structures are correctly tuned. A control algorithm that is simple, has a good performance over the full range of flow conditions and keeps the water levels at a predefined target level is the Proportional Integral (PI) feedback controller. Multiple successful implementations of this type of controller have been reported in the Western part of the US, Europe and Australia (Burt & Piao (2002), Malaterre & Baume (1998), Ooi & Weyer (2008)). Often, the water level measurement is located at the downstream end of the channel and this measurement signal is communicated to the upstream structure (see Figure 1). This setup is referred to as distant downstream control. It has the advantage that the offtakes, which are generally located at the downstream end of the open water channel, have an upstream water level that is (close to being) constant. The dynamic behavior of open water channels can be described with the well-known Saint-Venant equations (Chow (1959)). These are partial differential equations that consist of a mass balance and momentum equations. Some channels have structures, such as culverts, which can be described by another type of momentum equation. The momentum Peter-Jules van Overloop is with the Department of Water Management, Delft University of Technology, Delft, The Netherlands (e-mail: [email protected]). Xavier Bombois is with the Delft Center for Systems and Control, Delft University of Technology, Delft, The Netherlands (e-mail: [email protected]).
978-3-902823-06-9/12/$20.00 © 2012 IFAC
equations are non-linear. As control theory is mainly based on linear systems and ordinary differential equations, simplification of the dynamic behavior is required. In Schuurmans (1997), a division into two classes of open water channels is made. One class is defined as channels that are long and shallow. The dynamic behavior is mainly determined by the long delay times of the water that has to travel from the upstream side to the downstream end of the channel. In Schuurmans et al. (1999) and later in Litrico & Fromion (2006), research has been conducted to find stable and well-performing controllers for these channels without an extensive trial-and error tuning procedure. Their research has resulted in controller tuning rules for this class of open water channels based on the delay time τd and the storage area As of the back water part. The other class of channels can be characterized by being short and deep and is dominated by reflecting waves. After a change in the boundaries of the system, for example the opening of a gate or a pump that is turned off, a wave travels up and down the channel a number of times before the water surface settles again. Also in Litrico & Fromion (2004a), Overloop (2006) and Ooi & Weyer (2009), these so-called resonance waves are described and analyzed. The main reason for the badly damped behavior in short and deep channels is the low bed friction force along the open water channel. Figure 1 shows a representation of the properties of the two classes of open water channels. For the resonance class open water channel, only the resonance mode with the lowest frequency is shown (first peak in frequency spectrum). It is important to consider exactly this frequency, as the measurement location for distant downstream control is in counter-phase (phase shift of -180º) with the control actions. According to the Nyquist stability criterion, this control system will, in closed loop, be unstable, when the magnitude of this resonance peak, in open loop, is larger than 1. Schuurmans (1997) proposed to filter the resonance wave with a first-order low-pass filter and developed a tuning rule for the filter. This rule requires input of the frequency ωp and the magnitude of the first peak Mp of the resonance wave in the frequency spectrum. At present, stable and well-performing PI controllers for all classes of open water channels can be directly determined when the properties delay time τd, storage area As, frequency ωp and peak magnitude Mp of the lowest resonance mode can be estimated accurately at low and high flow conditions (Schuurmans (1997)). Other implementations of feedback controllers (Weyer 2003) require the first resonance frequency for selection of the cut-off frequency of the applied high-order low-pass filter.
1019
10.3182/20120711-3-BE-2027.00414
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Qin(k)
h(k) Qin(k-k d) As
Delay time τd
Qin(k)
h(k)
Storage area
Qout (k)
Qout (k)
Figure 2. Open water channel variables and parameters
Figure 1. Delay-time-dominated open water channel (left) and resonancedominated open water channel (right) including distant downstream control loop
In Overloop et al. (2010), a identification procedure for identifying the resonance properties of an open water channel in the Central Arizona Irrigation and Drainage District in Eloy, Arizona is presented. The identified model in this research gave satisfactory results for the resonance, but was inaccurate for low frequencies. The storage area, that needs to be identified from this low-frequency part, showed unlikely values twice as high as the pool’s top width multiplied by the length. In this article, the measurements of the Overloop et al. (2010)’ experiment are used to demonstrate an improved procedure that gives accurate values for both resonance frequency and the low-frequency range. It is a modification to the procedure described in Overloop et al. (2010). Therefore large parts of that article are repeated in this article. First, the general dynamic behavior of open water channels is described. Next, system identification methods for estimating the properties of an open water channel are discussed. After this, the applied system identification experiment on the actual open water channel is described. This section includes the results. The last sections are the conclusions of this article. II. OPEN WATER CHANNEL DYNAMICS The flows and water levels in an open channel can be described by the Saint-Venant equations (Chow (1959)). These nonlinear hyperbolic partial differential equations consist of a mass balance and a momentum balance. The momentum balance is a summation of the descriptions for the inertia (1), advection (2), gravitational force (3) and friction force (4): A f Q (Eq. 1) + = qlat x t (Eq. 2) g Q Q Q Q2 h t 1
+
+ =0 + g Af x Af x C 2 R f Af 2
3
4
To use the formulas in a numerical model of an open water channel, the partial differential equations are discretized in time (Δt) and space (Δx). In case these discretized formulas are simulated, the model results in time series solutions of water levels and flows at discrete locations along the channel. Also, the time series are discrete solutions in time. This numerical model is non-linear and of high order, which makes it hard to identify and use for controller design. When we consider the entire open water channel at once, on average, the water flows through an open water channel with an average velocity V (m/s) that can be calculated as the average flow Q divided by the average wetted area Af. A wave front on top of this velocity travels with a higher velocity that is characterized by the celerity c (m/s) (Chow (1959)). The celerity is commonly estimated from the average depth d (m) according to Equation 3. (Eq. 3) c g d The most compact description of an open water channel is the Integrator Delay model proposed by (Schuurmans (1997)), where the input is the upstream water flow and the output is the water level downstream in the pool. This Laplace model given in Equation 4, describes only the low-frequency part and neglects resonances. H s e d (Eq. 4) Qs As s where H(s) is the Laplace transform of the water level downstream in the pool and Q(s) is the Laplace transform of the flow of the upstream gate. For a given frequency ω the main property for controller design, storage area As, can be calculated from the magnitude M(ω) of this transfer function according to Equation 5.
1 (Eq. 5) M The equation is also valid in the frequency region of higher order models where the (low-frequency) integrator action is present. As
3
where Q represents the flow (m /s), t the time (s), x the distance (m), Af the wetted area of the flow (m2), qlat the lateral inflow per unit length (m2/s), g the gravitational acceleration (≈9.81 m/s2), h the water level (mMeanSeaLevel), C the Chézy friction coefficient (m1/2/s) and Rf the hydraulic radius (m). Rf is calculated by Af over Pf, where Pf represents the wetted perimeter (m). Figure 2 gives a schematic representation of a typical open channel with its variables and parameter.
III. IDENTIFICATION METHODS APPLIED TO OPEN WATER CHANNELS
Deriving low-order models and open water channel properties from analytical or numerical models can be risky, as there can be mismatches between actual channel and model. To identify the actual properties, field experiments can be executed. Below, various identification methods are described.
1020
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
A step test is commonly used to find the delay time τd and the storage area As as modeled in the Integrator Delay model. In this type of tests, starting from an initial steady state in the channel, a positive step change in flow is applied at the upstream structure, while the downstream flow is kept constant. The measured water level trajectory typically shows a delay time before the water level start to rise with a more or less constant slope that is determined by the size of the storage area. Figure 3 presents the result of a step test on an open water channel.
Figure 3. Step response of model of open water channel
Another field experiment to identify the delay time and the storage area is applied in the Auto Tune Variation proposed by Litrico et al. (2007). This method applies a switching pattern of upward and downwards steps. The switches are made when the water level reaches the allowed minimum or maximum of the band around the target level. The resonance properties, especially the resonance peak magnitude Mp, are harder to identify. In Figure 3, a resonance wave can be distinguished, super-positioned on the sloped trajectory. The resonance frequency can be roughly estimated according to Equation 6, assuming the celerity c is considerably higher than the flow velocity V.
p
c L
(Eq. 6)
where L (m) is the length of the open water channel. An experiment to identify the resonance peak magnitude can be executed by applying a Proportional feedback controller of which the gain is increased until instability occurs. The resonance peak is the reciprocal of the gain that gives the instability. Another experiment is the application of an input chirp signal (Miltenburg (2008)). A recent advance in estimating the resonance is described in Clemmens et al. (2012), where the ATV procedure is extended to include resonance waves. The previously described experiments aim at finding one or a subset of the properties that are required to tune a controller for resonance-dominated channels. Potentially, the engineering field of System Identification (Ljung (1999)), in which models with low order are estimated from input and output signals, can produce the full subset of required properties. Silvis et al. (1998) apply SI-techniques to identify
the Integrator Delay model (Schuurmans (1997)) from experimental data. In Eurén & Weyer (2007) the same model of an integrator and a delay in series is used. The difference between the two approaches is that the first article applies flows at the control structures as input variables, while the second article uses the positions of the undershot and overshot control gates. Rivas Perez et al. (2007) identify a second order model in series with a delay time. In Litrico & Fromion (2004b) the first step upwards is linked to the high frequency behavior of resonance-dominated water channels that consists of the first resonance and its higher harmonics. In that article, a low order model is proposed that contains a zero in the transfer function in order to model this highfrequency resonance behavior. Weyer (2001) applies a delay time in series with a third order model. The reasoning behind this is that this model captures the storage area behavior (first order) in parallel with the first resonance wave (second order oscillating mode). In Overloop et al. (2010), which was based on Miltenburg (2008), various model structures were tested on actual measurement data in order to find the structure that best captures the frequency response of the system around the first resonance mode. The best fit and most realistic resonance properties were found with a Box-Jenkins model of order five. The low-frequency part however resulted in an unlikely estimate of the storage area As of twice the surface area. IV. SYSTEM IDENTIFICATION APPLIED TO AN ACTUAL OPEN WATER CHANNEL
Experiments are conducted at the irrigation district office of the Central Arizona Irrigation and Drainage District in Arizona (www.caidd.com). The water level and flow measurements are sent to the central computer at a sample rate of once per 10 seconds (= Ts). The flow input signal is sent to the gates at a rate of once every 20 seconds. Communication is realized by means of radio-signals. As the open water channels are part of an operational delivery system, the experiments took place at night. Multiple experiments on various open water channels of the irrigation district are performed. In this article, the experiment on open water channel NB12 is presented, which is supposed to be resonance-dominated. Other tests can be found in Miltenburg (2008). The dimensions of the channel are presented in Table 1. Table 1. Dimensions of open water channel NB12 Length
Bed
Bottom
Side
Average
Base flow
(m)
slope (-)
width (m)
slope (-)
depth (m)
(m3/s)
2.0e-4
0.61
1.5
0.80
1021
786
0.2
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
These dimensions result in a situation where the entire open water channel is under backwater (Schuurmans (1997)), which indicates that the entire pool functions as a reservoir with a storage area of approximately the surface area of top width multiplied by length. Figure 4 presents the Bode diagram of a high-order linearized Saint-Venant model (Δx=10) of the resonance-dominated channel NB12 with a customary Mannings-friction coefficient of 0.02 s/m1/3 ( ≈ Chézy 43 m1/2/s ). In this Bode diagram, it can be observed that the frequency of the first resonance peak is ωest ≈ 0.0085 rad/s and that the magnitude at that frequency has an magnitude of approximately Mp,est ≈ 0.24. Note that the estimate for the frequency without using a Saint-Venant model, according to Equation 6 is 0.011 rad/s. From the lowfrequency part, the storage area can be calculated using Equation 5 as being As ≈ 2222 m2. This is close to the surface area of the top width multiplied by the length, which is 2366 m2. In the sequel, prediction error identification will be used to refine the model estimates of the first resonance frequency ωp and of its magnitude Mp, as well as the delay time τd and storage area As. The estimates of these quantities deduced from Figure 4 will be used in the design of the identification experiment. Note that the representation of the dynamic behavior as presented in Figure 4 comes from a model and can be inaccurate due to approximations and modeling assumption. Experiment setup In order to excite the system, a Random Binary Signal (RBS) is super-positioned on the upstream flow. The amplitude of this RBS is here approximately equal to 0.075 m3/s. In general, the amplitude should be of at least 25% of the base flow in order to guarantee a sufficient signal-to-noise ratio. The upstream flow perturbed in this manner is applied at the upstream gate, while keeping the downstream outflow constant. The initial condition in the channel has to be steady and the water level has to be close to target level. In order to contain enough information, the duration of the experiment is chosen to be equal to approximately 25 times the fundamental period 2.π/ωest of the estimated first resonance frequency. The duration of the experiment was therefore of approximately 5 hours and 20 minutes which corresponds to 1937 data points with the available sampling time of 10 seconds.
Figure 4. Bode diagram of high-order linearized transfer function from upstream flow to downstream water level
The to-be-identified system contains an integrator (see also low-frequency behavior in Figure 4). It is well known that identification of a system with an integrator is troublesome, as it is on the edge of being unstable. Consequently, the output data have to be differentiated using a tamed differentiator D. A tamed differentiator only differentiates the signal up to a certain frequency ωd in order to avoid any amplification of high frequency noise. In order not to affect the resonance peaks, ωd needs to be selected considerably lower than the resonance frequency. In Overloop et al. (2010) the frequency ωd was chosen equal to 0.002 rad/s, approximately five times smaller than the expected resonance frequency ωest. Recent research shows that this caused the storage area to be estimated inaccurately. A better result can be found by using ωd as 0.0001 rad/s, approximately 100 times smaller than the estimated resonance frequency. With the original sampling time Ts = 10 seconds, the frequency response of the system could be described up to 0.31 rad/s (the Nyquist frequency π/Ts). In this frequency range, the frequency response of the system is expected to present a large number of resonance peaks (see Figure 4). In order to avoid a model of a too large order and to avoid the classical problems arising when identifying dynamics in distant frequency ranges, the input-output data are decimated by a factor 10 (after the application of an anti-aliasing filter). The new sampling time is therefore equal to Ts* = 100 seconds. This choice for the decimation factor is the results of trial and error.
1022
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
^
function with sampling time Ts* = 100 seconds. The model W obtained in this way is (Ts* = 100s): ^
W z 2
0.0225 z 9 0.1365 z 8 0.06435 z 7 0.1743z 6 0.003077 z 5 0.03225 z 4 0.04087 z 3 0.001483z 2 0.001691z z 9 0.5689 z 8 0.6226 z 7 0.003737 z 6 0.02749 z 5 0.09991z 4 0.01092 z 3 0.1917 z 2 0.04559 z1 0.3161
(Eq. 9) The Bode diagram of the ninth order model is represented in Figure 6, where the integrator behavior in low frequencies and the first three resonance peaks can be observed. The model identified in Overloop et al. (2010), which is a fifth order model, is presented as well (dashed line). It can be seen that the ninth order has a higher low-frequency part that is in much better correspondence with the high-order numerical model as presented in Figure 4. Figure 5. Model residuals for the identified Box-Jenkins model
Model identification using prediction error identification Prediction error identification (Ljung (1999)) will now be used in order to deduce a model between the input signal and the (differentiated) output signal. A Box-Jenkins structure is used for this purpose. This model structure models the process transfer function as well as the noise perturbing the output. Even though a description of the noise is not required here, it was decided to model the noise as well, since it can be proven that this is more optimal statistically (see Ljung (1999)). In order to find a model order able to describe the underlying process, the order of the process model and of the noise model were increased in a stepwise manner until the identified models were deemed validated i.e. until the socalled residual tests (see Ljung (1999)) indicated that the prediction error could be both considered as uncorrelated with the input and as representing the characteristics of a white noise. At the same time the models needed to show a frequency response with three resonating peaks in the range of [0, 0.031].This procedure has led to the identification of the following process model (Ts* = 100s) of order eight for the process model, of order five for the noise model and containing two delay steps 1: ^
G z 2
Figure 6. Bode diagram of 9th (solid) and 5th (dashed) order transfer function of open water channel
Validation of the model Now, the validation data set to verify the accuracy of model ^
W is used. The input signal in this data set is simulated with ^
model W and this simulated output ysim is compared with the actual measured output signal yval in the validation data set (after de-trending). This delivers Figure 7. The fit defined as:
0.0225z 8 0.1142 z 7 0.1774 z 6 0.001342 z 5 0.004405z 4 0.02789 z 3 0.01325z 2 0.001708z1 z 8 0.4311z 7 0.1916 z 6 0.1953z 5 0.1678z 4 0.06792 z 3 0.07884 z 2 0.2706 z1 0.3161
(Eq. 7) and the following noise model:
FIT (%) 1
z 0.08885 z 0.6504 z 0.1828 z 0.3607 z 0.17 H 8 z 0.5921z 4 0.5962 z 3 0.508 z 2 0.3378 z1 0.6165 5
^
4
3
2
1
(Eq. 8)
where z-1 is the discrete backwards step of Ts* seconds. Figure 5 presents the residual tests corresponding to the ^
^
model { G , H } and shows that this model is indeed validated. ^
The model G is a model relating the applied flow input to the differentiated water level. In order to determine a model relating effectively the flow and the water level, the identified
1 N
N
y
val
k 1
1 N
(k ) ysim (k )
N
y k 1
val
2
(Eq. 10)
(k )2
is here equal to 89%. Resulting properties Finally, in this results section, the properties are derived from the identified ninth order model in comparison with the fifth order model as found in Overloop et al. (2010). Table 2 presents the determined properties.
^
model G above must be multiplied by a the inverse of D*. D* is the tamed differentiator D re-sampled to produce a transfer
1 The number of delays in the model has been determined by trial and error and by considering the dimensions of the considered open water channel as given in Table 1.
1023
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Table 2. Properties of open water channel found from identification experiment Model
5th order 9th order
Estimated delay time estimated τd (s) 200
Estimated storage area As (m2) 5217
Estimated resonance frequency ωp (rad/s) 0.0088
Estimated resonance peak magnitude Mp 0.223
200
2668
0.0088
0.249
The storage area is calculated using Equation 5 at a frequency of 1e-3 rad/s. The ninth order model is now well capable of identifying this area close to the surface area of top width multiplied by length of the open water channel (2366 m2).
Figure 7. Validation of estimated model (solid) and measurements (dashed) of the open water channel (top-graph is input signal flow around base flow at upstream side of the channel, while bottom-graph is output signal water level around average depth at downstream side).
V. CONCLUSIONS A field experiment to identify all properties of an open water channel required for controller design is successfully tested on an actual open water channel in Arizona. In this case, a ninth order model was required to capture all relevant dynamics at low frequencies and at the frequency of the first resonance peak. The main improvement compared to the procedure as presented in earlier research (Overloop et al. (2010)) was the capability of identifying the low-frequency part of the transfer function between upstream flow and downstream water level. Because of this, the storage area can now be estimated accurately.
Eurén, K., Weyer, E. (2007), ‘System identification of open water channels with undershot and overshot gates’, Control Engineering Practice, Vol 15, p813-824. Litrico, X., Fromion, V. (2004a), ‘Frequency modeling of open-channel flow’, Journal of Hydraulic Engineering, Vol130, Issue 8, p806-815 ASCE. Litrico, X., Fromion, V. (2004b), ‘Analytical approximation of openchannel flow for controller design’, Journal of Applied Mathematical Modelling, Vol 28, Issue 7, p677-695. Litrico, X., Fromion, V. (2006), ‘Tuning of robust distant downstream PI controllers for an irrigation canal pool. I: Theory’, Journal of Irrigation and Drainage Engineering, Vol 132, Issue 4, p359-368, ASCE. Litrico, X., Malaterre, P.O., Baume, J.P., Vion, P.Y., Ribot-Bruno, J. (2007), ‘Automatic tuning of PI controllers for an irrigation canal pool’, Journal of Irrigation and Drainage Engineering, Vol 133, Issue 1, p27-37, ASCE. Ljung, L. (1999), ‘System Identification, theory for the user’, New Jersey, PTR Prentice Hall. Malaterre, P.O., Baume, B.P. (1998), ‘Modeling and regulation of irrigation canals: existing applications and ongoing researches’. Proceedings of IEEE Conference on System, Man and Cybernetics, San Diego, Vol 4, p3850-3855. Miltenburg, I.J. (2008), ‘Determination of canal pool properties by experimental modeling’, Master-thesis Delft University of Technology, The Netherlands. Ooi, S.K., Weyer, E. (2008), ‘Control design for an irrigation channel from physical data’, Control Engineering Practice, Vol 16, p11321150. Ooi, S.K., Weyer, E. (2009), ‘Detection of oscillations in control loops in irrigation channels’, Control Engineering Practice, doi:10.1016/j.conengprac.2009.04.010 Overloop, P.J. van (2006), ‘Model Predictive Control on open water systems’, Ph.D.-dissertation Delft University of Technology, The Netherlands. Overloop, P.J. van, I.J. Miltenburg, X. Bombois, A.J. Clemmens, R.J. Strand, N.C. van de Giesen, R.W. Hut (2010), ‘Identification of resonance waves in open water channels’, Control Engineering Practice. Rivas Perez, R., Feliu Battle, V., Sanchez Rodriguez, L. (2007), ‘Robust system identification of an irrigation main canal’, Advances in Water Resources, Vol 30, Issue 8, p1785-1796. Schuurmans, J. (1997), ‘Control of water levels in open channels’, Ph.D.-dissertation Delft University of Technology, The Netherlands. Schuurmans, J., Clemmens, A.J., Dijkstra, S.J., Hof, A. & Brouwer, R. (1999), 'Modeling of Irrigation and Drainage Canals for Controller Design', Journal of Irrigation and Drainage Engineering, Vol 125, Issue 6, p338-6344, ASCE. Silvis, L.G., Hof, A., Hof, P.M.J. & Clemmens, A.J. (1998), ‘System identification on open channels’, Proceedings 3rd International Conference on Hydroinformatics. Weyer, E. (2001), ‘System identification of an open water channel’, Control Engineering Practice, Vol 9, Issue 12, p1289-1299. Weyer, E. (2003), ‘LQ control of an irrigation channel’, Proceedings of the 42nd IEEE Conference on Decision and Control, Maui, Hawaii USA, p750-755.
REFERENCES Burt, C.M., Piao, X. (2002), ’Advances in PLC-Based Canal Automation’, Proceedings USCID-Conference on Benchmarking Irrigation System Performance Using Water Measurement and Water Balances, San Luis Obispo, California. Chow, V.T. (1959), ‘Open-Channel hydraulics’, McGraw-Hill Book Co. Inc, New York. Clemmens, A.J., Litrico, X., Overloop, P.J. van, Strand, R.J., (2012), ‘Estimating Canal Pool Resonance with Auto Tune Variation’, Journal of Irrigation and Drainage Engineering.
1024