Instrumentation, model identification and control of an experimental irrigation canal

UNIVERSITAT POLITÈCNICA DE CATALUNYA Doctoral Programme: AUTOMATITZACIÓ AVANÇADA I ROBÒTICA Doctoral Thesis INSTRUMENTATION, MODEL IDENTIFICATION AND CONTROL OF AN EXPERIMENTAL IRRIGATION CANAL Carlos Alberto Sepúlveda Toepfer Supervisors: PhD Manuel Gómez V. and PhD José Rodellar B. Institut d’Organització i Control de Sistemes Industrials October 2007 To my two Silvias Abstract This thesis aims to develop control algorithms for irrigation canals in an experimental framework. These water transport systems are difficult to manage and present low efficiencies in practice. As a result, an important percentage of water is lost, maintenance costs increase and water users follow a rigid irrigation schedule. All these problems can be reduced by automating the operation of irrigation canals. In order to fulfil the objectives, a laboratory canal, called Canal PAC-UPC, was equipped and instrumented in parallel with the development of this thesis. In general, the methods and solutions proposed herein were extensively tested in this canal. In a broader context, three main contributions in different irrigation canal control areas are presented. Focusing on gate-discharge measurements, many submerged-discharge calculation methods are tested and compared using Canal PAC-UPC measurement data. It has been found that most of them present errors around 10 %, but that there are notable exceptions. Specifically, using classical formulas with a constant 0.611 contraction value give very good results (MAPE<6 %), but when data is available, a very simple calibration formula recently proposed in Ferro (2001) significantly outperform the rest (MAPE<3 %). As a consequence, the latter is encouragingly proposed as the basis of any gate discharge controller. With respect to irrigation canal modeling, a detailed procedure to obtain data-driven linear irrigation canal models is successfully developed. These models do not use physical parameters of the system, but are constructed from measurement data. In this case, these models are thought to be used in irrigation canal control issues like controller tuning, internal controller model in predictive controllers, or simply as fast and simple simulation platforms. Much effort is employed in obtaining an adequate model structure from the linearized Saint-Venant equations, yielding to a mathematical procedure that verifies the existence of an integrator pole in any type of canal working under any hydraulic condition. Time-domain and frequency-domain results demonstrate the accuracy of the resulting models approximating a canal working around a particular operation condition, in either simulation or experimentation. Regarding to irrigation canal control, two research lines are exploited. First, a new water level control scheme is proposed as an alternative between decentralized and centralized control. It is called Semi-decentralized scheme and aims to resemble the centralized control performance while maintaining an almost decentralized structure. Second, different water level ii Abstract control schemes based on Proportional Integral (PI) control and Predictive Control (PC) are studied and compared. The simulation and laboratory results show that the response and performance of this new strategy against offtake discharge changes, are almost identical to the ones of the centralized control, outperforming the other tested PI-based and PC-based schemes. In addition, it is verified that PC-based schemes with good controller models can counteract offtake discharge variations with less level deviations and in almost half the time than PI-based schemes. In addition to these three main contributions, many other smaller developments, minor results and practical recommendations for irrigation canal automation are presented throughout this thesis. Acknowledgments There are lots of people that I have to thank after four years of research, study, hard work and unforgettable experiences. First, I would like to express my enormous gratitude to my two thesis supervisors, Professor Manuel Gómez and Professor José Rodellar, for offering me the great opportunity to study at UPC, a prestigious university located in a beautiful city, doing something that I really enjoy. Their wisdom and expert guidance have been indispensable in the elaboration of this thesis and without their efforts, I would not have had the support to finish this work. It was a pleasure to work under the supervision of these two true gentlemen. I would like to thank my mentor at UdeC, Professor Daniel Sbarbaro, for his help and guidance in the embryonic stage of the research. Thanks to the laboratory personnel of the Hydraulic and Hydrological Engineering Section, Juan Pomares, Jaime Ambrós, Robert McAllon and others for their help and meticulous work in the construction of the laboratory canal. Appreciation is also extended to my "occasional" laboratory assistants and friends: Silvia Arriagada, Joaquim Blesa, Rodrigo Concha, Claudiu Iurian, Gustavo Mazza and Francisco Núñez. Their help was invaluable in particular periods of this work and I think I will never be able to reward the patience they had with a so perfectionist and demanding boss. I would need more than a lifetime to thank all what my beloved wife, Silvia Arriagada, has made for me. She has worked a lot to make this possible. Without her, I would not have been capable to bring this to an end. This thesis is also hers. Additionally, I wish to thank for the generous support provided by the Spanish Ministry of Science and Education through BES-2003-2042 grant. Finally, special thanks go to my parents, Silvia Toepfer and Luis Sepúlveda; brothers, Lucho, Jorge and Rodrigo; friends and colleagues, Anaïs, Belén, Cesca, Kat, Raquel, Úrsula, Andrés, Alejandro, Antonio, Beniamino, Carles, David, Germán, Hans, Jordi, José Luis, Quim and Vicente for their continuous encouragement and support. People like this is essential to make research possible. Carlos Sepúlveda Toepfer Barcelona October 26, 2007 Contents Abstract i Acknowledgments iii Contents v List of Figures ix List of Tables xiii Nomenclature xv List of Abbreviations xvii Glossary xix 1 Introduction 1 1.1 Background.................................... 1 1.2 Objectiveofthethesis............................... 4 1.3 Methodologyused ................................ 4 1.4 Thesislayout ................................... 5 2 Literature Review 7 2.1 System characteristics and the control problem . . . . . . . . . . . . . . . . . 7 2.2 Irrigation canal models for automatic control purposes . . . . . . . . . . . . . 8 2.3 Types of control algorithms developed for irrigation canals . . . . . . . . . . . 12 2.3.1 PIDcontrol................................ 13 2.3.2 Robustcontrol .............................. 15 2.3.3 Optimalcontrol.............................. 15 2.3.4 Predictivecontrol............................. 16 2.3.5 NonlinearControl ............................ 17 2.4 Conclusion .................................... 17 xii LIST OF FIGURES 6.29 Initial condition for the performance test . . . . . . . . . . . . . . . . . . . . . 190 6.30 Control schemes where controllers do not share information: Regulation of waterlevels...................................... 193 6.31 Control schemes where controllers share information: Regulation of water levels 194 6.32 Control schemes where controllers share information and the disturbances are measured: Regulation of water levels . . . . . . . . . . . . . . . . . . . . . . . 195 6.33 Calculated gate discharges . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 6.34 Applied gate openings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 6.35 Comparison of the overall regulation performance when using different controls 201 List of Tables 4.1 Calibration results for Gate 1, Gate 3 and Gate 5 . . . . . . . . . . . . . . . . . 56 4.2 Performanceindices................................ 60 4.3 Performance of each method for Gate 1 data . . . . . . . . . . . . . . . . . . . 60 4.4 Performance of each method for Gate 3 data . . . . . . . . . . . . . . . . . . . 61 4.5 Performance of each method for Gate 5 data . . . . . . . . . . . . . . . . . . . 61 4.6 Typical calibration values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.1 Reach’s Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.2 PRBSparameters................................. 113 6.1 Determination of the ξparameters ........................ 135 6.2 PI controller tuning rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 6.3 Åström-Hägglund PI tuning rule . . . . . . . . . . . . . . . . . . . . . . . . . 149 6.4 ID model parameters for Pool 1 and Pool 2 . . . . . . . . . . . . . . . . . . . 149 6.5 Ultimate cycle parameters for Pool 3 . . . . . . . . . . . . . . . . . . . . . . . 149 6.6 PI parameters tuned using different methods . . . . . . . . . . . . . . . . . . . 150 6.7 Retuned PI and PIF controller parameters in order to ensure stability against measurementnoise ................................ 151 6.8 Tuning values for each predictive controller . . . . . . . . . . . . . . . . . . . 157 6.9 Tuning values for 3×3predictive controller . . . . . . . . . . . . . . . . . . . 169 6.10 Tuning values for each 2×2predictive controller . . . . . . . . . . . . . . . . 175 6.11 Determination of the ξparameters ........................ 187 B.1 Weir 1 calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 B.2 Weir 3 calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 B.3 Weir 4 calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222 C.1 Gate 1 calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223 C.2 Gate 3 calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224 C.3 Gate 5 calibration data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 225 Nomenclature ΦkVector of past input and output measures at instant k ΘVector of model parameters IIdentity matrix νKinematic viscosity QiReach’s i, transport area, water discharge downstream boundary condition X x position of the end of the transport area (the start of the storage area) ZiReach’s i, transport area, water level downstream boundary condition AWetted cross-sectional area BWater width bGate or weir width CWater wave celerity CcContraction coefficient CdDischarge coefficient Crw Rectangular weir discharge coefficient Cvw V-notch weir discharge coefficient dDelay in amount of instants ekModel error at instant k FFroude number F1Upstream gate discharge transfer function F2Downstream gate discharge transfer function F3Offtake discharge transfer function gGravity acceleration hHead above the weir h1Gate upstream water level xvi Nomenclature h2Water depth at vena contracta generated by the gate h3Gate downstream water depth kDiscrete time instant variable lGate opening liReach’s iupstream gate opening OWeir height or gate position PWetted perimeter QVolumetric water discharge qForward shift time operator qWater discharge around a working point q−1Backward shift time operator Q0Initial time water discharge condition QiReach’s i, transport area, water discharge upstream boundary condition Qi+1 Water discharge delivered to reach i+ 1 QL i Reach’s iofftake discharge sLaplace variable S0Bottom slope SfFriction slope TSampling period tTime uSystem input VWater velocity Vs i Water volume stored behind gate i+ 1 xLongitudinal coordinate in the flow direction ySystem output ZWater level zWater level around a working point zZ-transform complex variable Z0Initial time water level condition ZiReach’s i, transport area, water level upstream boundary condition zkOutput measure at instant k Zs i Reach’s idownstream water level List of abbreviations ARIMAX Auto-Regressive Integrated Moving Average with eXogenous Input ARIX Auto-Regressive Integrated with eXogenous input ARX Auto-Regressive with eXogenous Input ASCE American Society of Civil Engineers ASTM American Society for Testing and Materials BIBO Bounded-Input Bounded-Output Canal PAC-UPC Canal de Prueba de Algoritmos de Control - Universitat Politècnica de Catalunya CARDD Canal Automation for Rapid Demand Deliveries DAQ Data AcQuisition HEC-RAS Hydrologic Engineering Center - River Analysis System ID Integrator Delay IDZ Integrator Delay Zero ISO International Organization for Standardization LAD Least Absolute Deviations LQ Linear Quadratic LQG Linear Quadratic Gaussian LQR Linear Quadratic Regulator MAE Mean Absolute Error MAPE Mean Absolute Percentage Error ME Mean Error MIMO Multiple-Input Multiple-Output MPE Mean Percentage Error ODE Ordinary Differential Equation xviii List of Abbreviations PProportional PC Predictive Control PDE Partial Differential Equation PI Proportional Integral PID Proportional Integral Derivative PRBS Pseudo Random Binary Sequence PV Process Variable QP Quadratic Programming RMSE Root Mean Square Error SCADA Supervisory Control And Data Acquisition SIC Simulation of Irrigation Canals SISO Single-Input Single-Output SP Set Point TITO Two-Input Two-Output USBR U.S. Department of the Interior, Bureau of Reclamation WLS Weighted Least Squares Glossary B Black-box model No physical insight is available or used, but the chosen model structure belongs to families that are known to have good flexibility and have been "successful in the past", p. 10. bode diagram Logarithmic magnitude and phase plot of a transfer function, that gives information of this function evaluated in the s-plane imaginary axis. In a more practical view, it shows what happens with the amplitude and the phase of the response of a system, when it is excited with a sinusoidal input at a given frequency. C Control theory Branch of Mathematics and Engineering which deals with the design, identification and analysis of systems with a view towards controlling them, i.e., to make them perform specific tasks or make them behave in a desired way. G Grey-box model This is the case when some physical insight is available, but several parameters remain to be determined from observed data. It is useful to consider two sub-cases: - Physical Modeling: A model structure can be built on physical grounds, which has a certain number of parameters to be estimated from data. This could, e.g., be a state space model of given order and structure. - Semi-physical modeling: Physical insight is used to suggest certain nonlinear combinations of measured data signal. These new signals are then subjected to model structures of black box character, p. 10. H hydraulic jump The sudden and usually turbulent passage of water in an open channel from low stage, below critical depth, to high stage, above critical depth. During this passage, the velocity changes from supercritical to subcritical. There is considerable loss of energy during the jump, p. 46. xx GLOSSARY I irrigation canal Permanent irrigation conduit constructed to convey water from the source of supply to one or more farms, p. 8. P Parameter Estimation Estimating the values of parameters based on measured/empirical data. R radial gate Gate with a curved upstream plate and radial arms hinged to piers or other supporting structures on which the gate pivots. S sluice gate Gate that can be opened or closed by sliding it in supporting guides. state-space model Mathematical model of a physical system as a set of input, output and state variables related by first-order differential equations. To abstract from the number of inputs, outputs and states, the variables are expressed as vectors and the differential and algebraic equations are written in matrix form, p. 8. System Identification General term to describe mathematical tools and algorithms that build dynamical models from measured data, p. 86. T transfer function Mathematical representation of the relation between the input and output of a linear time-invariant system, p. 9. W weir 1. Low dam or wall built across a stream to raise the upstream water level, termed fixed-crest weir when uncontrolled. 2. Structure built across a stream or channel for the purpose of measuring flow, sometimes called a measuring weir or gauging weir. Types of weirs include broad-crested, sharp-crested, drowned, and submerged. White-box model This is the case when a model is perfectly known; it has been possible to construct it entirely from prior knowledge and physical insight. Chapter 1 Introduction 1.1 Background Irrigation is the artificial application of water to the soil usually for assisting in growing crops. In crop production, irrigation is mainly used to replace missing rainfall in periods of drought, but also to protect plants against frost. At the global scale, approximately 2 788 000 km2of agricultural land is equipped for irrigation in the world. 68 % of this area is located in Asia, 17 % in America, 9 % in Europe, 5 % in Africa and 1 % in Oceania. Most of this vast area is gridded by irrigation canals. Irrigation canals are artificial systems developed to transport water from main water reservoirs to several water-demanding agricultural farms during irrigational seasons (see figure 1.1). Generally, they cover very long distances: their length can range from hundreds of meters to hundred of kilometers. Along these canals, farms are located close to them and distributed all over the way. A typical configuration, considered as a prototype case study in this thesis, is the one depicted in figure 1.1b. A main canal transports water from a big reservoir to the farms and controls the water flow by modifying the openings of several gates. These hydraulic structures are situated in the waterway in order to regulate discharge in relation to ongoing irrigation demands. Water is supplied only in fixed locations of the canal; generally, a few meters upstream of each gate. In real situations, these water offtakes are usually performed by pumps or weirs. It is not trivial to manage this type of systems. Water must be transported, minimizing the losses and assuring that every farm receives the stipulated amount of water at its corresponding frequency. Besides, the inherent characteristics of these systems increase the complexity of the problem. These systems present very long delays in the water transport (from minutes to hours), delay that even varies depending on the provided discharge. Moreover, there are important dynamical effects produced generally by changes in the amount of supplied water, that produce 8 Chapter 2. Literature Review (1999), it is also usual to solve the control problem using only discharges, and afterwards use local controllers or discharge formulas inversion, in order to obtain the necessary actions over the actuators. In the operational point of view, generally, it is more often required a regulation effect, in front of previous known demands, than a change in the operational working point (Clemmens et al., 1998). However, the system has also unknown disturbances due to: inaccuracies in the measurement of the supplied discharge, filtration of the canals, non-authorized water extractions and changes in the demand. A frequency analysis around a given working point, as the one made in Litrico and Fromion (2004c), gives valuable information about the different types of behaviors that can appear. According to the geometrical characteristics of the canal and the hydraulic conditions of the flow that circulates through it, the system can exhibit (small slope canals) or not (high slope canals) resonant modes, can have long delays (whose values depend mainly on the canal’s length) that vary according to the circulating water discharge and can also have pure integrating dynamics (single pole in the origin). 2.2 Irrigation canal models for automatic control purposes The modeling of an irrigation canal is carried out dividing the canal into pools (section of a canal between two gates or any similar device), characterizing then the water dynamics at each reach separately and, finally, including the water regulation devices equations as boundary conditions between reaches. As detailed in Henderson (1966), the water flow through a reach can be well characterized by the Saint-Venant equations, a nonlinear hyperbolic Partial Differential Equation (PDE) system. On the other hand, the governing equations of the devices that are usually found in an irrigation canal (gates, weirs, etc.) are of nonlinear nature. In this manner, the solution of a complete irrigation canal model does not exist analitically and can only be done by means of advanced numerical methods (finite volume, finite differences, method of characteristics, etc.). Clearly these models are not adequate for their use in automatic control design and implementation. For this reason, a series of authors have proposed different and diverse simplified types of models for control. In Malaterre and Baume (1998) there is a survey about all types of models that had been used, until that date, in the canal control literature. They cover a large spectrum that includes: Saint-Venant model linearizations, infinite order linear transfer functions, finite order nonlinear models, finite order linear models (state-space models), finite order linear models (transfer functions), neural networks based models, fuzzy logic based models and petri nets based models. In all these alternatives, Single-Input Single-Output (SISO) approaches that model each 2.2. Irrigation canal models for automatic control purposes 9 reach separately and Multiple-Input Multiple-Output (MIMO) approaches that model a whole canal, have been used. In both cases, some models include the actuators dynamics and models that do not. In the last years, the literature shows the inclusion of new models and improvements to the already existing ones. First of all, we will refer to models that use some physical knowledge about the canal for its formulation. Second, we will review some black-box models along with identification techniques. In Schuurmans et al. (1999b,a) a model, proposed by the same authors in 1995, is evaluated. This one approximates each canal’s reach as a pure integrator plus a delay, the reason why it is called Integrator Delay (ID) Model. The input variables of the model are the reach’s inflow and outflow discharges and what is obtained is the water level at the end of it. In this model, the delay is obtained in an algebraic manner as a function of the physical parameters of each reach, and the storage area (integrative part) is obtained by means of the canal’s backwater curve. In order to include the actuators, linearized models of them were used. The model’s validation in the time domain, with experimental data, showed an acceptable performance when the system was operated with small movements around a working point. In the frequency domain, the model showed a good fit in the low frequencies, but a bad fit in the high frequencies. In other words, there exist some evidence that the model does not perform well in the short-term period. For this reason, the model manifested an incapability to approximate resonant modes when they exist. They put emphasis on remarking, however, that the modeling of these modes is not so important, because they are generally filtered in control applications. This model has been also used to generate state-space MIMO models and in Clemmens and Schuurmans (2004a,b); Wahlin (2004); Montazar et al. (2005) and van Overloop et al. (2005). Few years ago, improvements to the ID model have been also proposed in Litrico and Fromion (2004a,b). There, the inclusion of a transfer function to approximate better the high frequency range was proposed. The model was called Integrator Delay Zero (IDZ) Model and this additional transfer function was considered for the influences of the inflow and outflow discharges. In that work, the algebraic expressions that describe the model parameters were also modified. Another model in the same line as the preceding ones was the one presented in Rodellar et al. (1993) and in Gómez et al. (2002), where the Muskingum model was used to model the water transport and, also, an integrator was used to characterize the water level variations upstream a gate (extraction zone). Another approach was the one used in Litrico and Georges (1999a,b), where a simplification of the Saint-Venant equations was used to model a reach by means of the Hayami model. Due to the similarities that the authors observed between this model and a standard second-order plus delay one, they used the Method of Moments to obtain the parameters of this last one as a 10 Chapter 2. Literature Review function of the Hayami parameters. For the particular modeling case where rivers are used for irrigation purposes, in Litrico (2001a,b) system identification techniques were used to obtain the parameters of a Diffusive Wave model (another simplification of the Saint-Venant model) with the aid of experimental data. For canals, in Litrico and Fromion (2004c) and in Litrico et al. (2005) they developed and used a methodology for obtaining numerically the frequency response of a reach, including the gates, by means of the linearization of the Saint-Venant equations around an operation point and the knowledge of the hydraulic and geometric parameters of the canal. A different approach was used, for example, in Malaterre and Rodellar (1997) and in Malaterre (1998). In that paper, a state-space MIMO model was generated, using the linearized Preissmann method in order to solve the Saint-Venant equations directly and construct, in that manner, a state observer. A linearized version of all gates equations were also included in order to generate a model that, by knowing the gate’s openings, can calculate all water levels in the irrigation water extraction zones. This approach includes all system coupling effects and in general generates very large matrices. In the same line of thought, Reddy and Jacquot (1999) used a linearization of the SaintVenant equations using the Taylor series and a finite difference approximation to develop statespace MIMO model. A Kalman filter was also designed to estimate values for the state variables that were not measured. In Durdu (2005), they developed a state-space MIMO model using another finite difference method. The difference was that in that work they developed a state estimator based on fuzzy logic rules. Another state-space MIMO model was used in Seatzu (1999, 2000) and in Seatzu and Usai (2002). In that case the modeling was performed around a particular hydraulic regime, called uniform regime, that is the only one that has an algebraic solution by linearizing the SaintVenant equations. Additionally, the model used as inputs, gate openings, and as outputs, not the water depth levels, but rather, the storaged water volumes in each reach. Nonlinear irrigation canal models for control has also been developed. In Dulhoste et al. (2004), a model was developed by means of a nonlinear approximation with Lagrange polynomials of the Saint-Venant equations. In de Halleux et al. (2003), they went one step ahead and used the Saint-Venant model, but only for a zero-slope rectangular canal without friction. In Sanders and Katopodes (1999) the canal was modeled solving numerically the original SaintVenant equations as an adjoint problem discretized with the Leap-Frog scheme. In Soler et al. (2004) a numerical scheme solving the Saint Venant equations using the method of the characteristics has been developed to calculate desired trajectories for control gates. The modeling problem has also been faced from the Black-box model and Grey-box model 2.2. Irrigation canal models for automatic control purposes 11 identification point of view. In Akouz et al. (1998); Ruiz and Ramírez (1998); Sawadogo et al. (1998, 2000); Rivas et al. (2002) and in Rodellar et al. (2003) Auto-Regressive Integrated with eXogenous input (ARIX) and Auto-Regressive Integrated Moving Average with eXogenous Input (ARIMAX) Black-box models were used without getting too deep into the analysis or validation issues. The majority of these models used or the discharge or the gate opening at the beginning of the reach as model input and, the water level at its end, as output. In some cases the reach’s outflow discharge and the water delivered for irrigation (when it was initially known) were used as known disturbances. All of them used data obtained by computational simulation of the Saint-Venant equations. In Weyer (2001) a more deeper work was performed for model identification of the Haughton Main Channel reaches in Australia. Three gray-box models were used, based on elementary mass balances and gate equations. Because, in this case, the canal has overflow gates, the inputs to the model were water levels over the gates and its outputs, water levels in the extraction zones. Linear and nonlinear first order, second order and second order plus integrator (third order) models were proven. All of them included explicitly the delay. The obtained results, by means of model validation against real data, showed that the only model that could reproduce the effect of the waves was the nonlinear third order model. However, the first and second order ones could follow the tendencies in most of the cases. In the final conclusions they emphasized the need to study more the cases with gates in submerged regime and the use of closed loop identification, in order to estimate models using smaller variations and shorter experimentation times. The results of this work were extended in Eurén and Weyer (2007) in several aspects: •The irrigation channel was equipped with both overshot and undershot gates. •The overshot gates operated in both submerged and free flow. •There were several gates at each regulator structure and they had different positions. •The flows and pools were larger. The results presented in this paper are very encouraging, since the system identification models were able to accurately simulate the water levels more than 12 h ahead of time. System identification applied on irrigation canals was also studied in Ooi et al. (2005). The results showed that the St. Venant equations can adequately capture the dynamics of real channels, but that to estimate their parameters from real data is more accurate than using only physical knowledge. On the other hand, system identification models were as accurate as the St. Venant equations with estimated parameters and, consequently, they should be preferred over the St. Venant equations for control and prediction purposes since they are much easier to use. 12 Chapter 2. Literature Review 2.3 Types of control algorithms developed for irrigation canals Malaterre et al. (1998); Malaterre and Baume (1998) and Ruiz et al. (1998) gave a survey of the control algorithms that had been developed until 1998 for canal irrigation control. They cover a large spectrum of approaches and techniques, among which can be mentioned: monovariable heuristical methods, Proportional Integral Derivative (PID) Control, Smith Predictor scheme, Pole Placement Control, Predictive Control, Fuzzy Logic Control, Model Inversion methods, Optimization methods, Robust Control, Adaptive Control and Nonlinear Control. Due to the diversity of proposed methods and distinct performance criteria used, the American Society of Civil Engineers (ASCE) Task Committee on Canal Automation Algorithms developed in Clemmens et al. (1998) two standard cases (Test Canal 1 and Test Canal 2) to test and evaluate automatic control algorithms. These cases are based on real canals and normal operation conditions as, for example, scheduled and unscheduled water discharge offtakes and correct and incorrect knowledge of the canal physical parameters. In that work, a series of evaluation criteria are likewise given in order to standardize the evaluation of control algorithm performance. In spite of the important amount of studies that have been done about the subject (the majority of them in computational simulation), as denoted by Rogers and Goussard (1998) and by Burt and Piao (2004), until now the few real canals that are managed in an automatic form use, in their majority, at the most PID control based techniques. It has been used, aside from several heuristic techniques, in the form of Proportional Integral (PI) control in many cases, PI plus filter (PIF) in cases with resonance problems, and occasionally PID. These developments can be found in North America, Asia and Europe, but mostly in the USA. Going back to the academical knowledge developed, it is also important to mention that some of the cited methods have used only feedback strategies and others only feedforward strategies, while others have made use of combinations of both (Malaterre et al., 1998). Feedback produces a corrective control action in order to return the controlled variable to its nominal value, inclusively in presence of unknown disturbances, whereas feedforward can compensate the inherent delays of the system by anticipating the needs of the canal users. In Bautista and Clemmens (1999) they tested a classical open-loop method, called Gate Stroking. The conclusion was that an adequate irrigation canal controller should be implemented, when possible, with feedback and feedforward capabilities. That is especially true for canals that need large water volume variations, to arrive to another steady state condition, and for those characterized by a low Froude number. 2.3. Types of control algorithms developed for irrigation canals 13 2.3.1 PID control From 1998 on, the works based on PID control have focused in improving the tuning of this types of controllers. To achieve this goal, two common practices have been identified from the literature: the use of simplified mathematical models and the employment of strategies that lead to decoupling the influences, produced by the control action of a reach over all the adjacent ones. Schuurmans et al. (1999b) proposed a control where every reach was controlled by its upstream gate. In order to achieve this, a supervisory control was used. It calculated which should be the input discharge, and then a local controller moved the gate so as to obtain the required discharge. The reach’s outflow discharge and the water demands were also included as known disturbances, so as to decouple, in a better way, the interaction between reaches. The philosophy was, thus, to include the local controller in order to minimize the nonlinear effects that a gate induces on the canal operation. The tuning of these controllers was based on the ID Model and a filter was also used, so as to filter the canal inherent resonance. In Malaterre and Baume (1999) optimum PI controller tunings were calculated in conjunction with their corresponding performances for different cases. In that work, different manipulated variable choices and decoupling strategies were tested. They arrived at the conclusion that the best results are obtained with a supervisory control that calculates for each reach their optimum inflow discharge, and a series of local slave controllers that calculate the gate openings taking into account the water level variations induced by the gate movements. Better result were also obtained when the local controllers ran at a sampling time 5 times faster than the supervisory controller, but with a considerable increase in the control effort. Other work that treated the decentralized Proportional (P) and PI controller tuning was Seatzu (1999). They proposed the use of a state feedback diagonal matrix and a H2norm minimization, so as to obtain an optimal tuning. Seatzu (2000) proposed the same scheme, now seeking to place the eigenvalues and eigenvectors of the controlled MIMO system to some optimal values, obtained previously with the Linear Quadratic Regulator (LQR) method. Besides, in a later work (Seatzu and Usai, 2002), a research was made in order to know when this controller plus an observer was robust against modeling errors. In Wahlin and Clemmens (2002) three classical controllers were tested on the ASCE Test Canal 1: PI Control, PI Control with upstream and downstream decouplers as proposed by Schuurmans in 1992 (they were tested separatelly and together) and a heuristic control called Canal Automation for Rapid Demand Deliveries (CARDD). In all cases the control variables were the gate openings and feedforward control was implemented by means of a volume compensation method. The results showed that the best option was the control with both decouplers and that feedforward strategy was indispensable for all the cases. A control deterioration was also ob- 14 Chapter 2. Literature Review served when the canal parameters were not accurate and when the gate movement restrictions were included. Afterwards, in Clemmens and Schuurmans (2004a) and in Clemmens and Schuurmans (2004b), a methodology was developed and tested in which, using a modified PI structure in order to compensate the delay of each reach, they structured a state feedback matrix with some non-zero elements. Then, formulating a LQR objective function and solving the Ricatti equation, they found the parameters of this feedback. They identified that using the trick of making zero some elements of the feedback matrix, was equivalent to use different decoupling logics and that the use of the whole matrix was equal to implement a completely centralized controller. They also made a performance study for all possible controllers, going from the complete centralization to the total decentralization. The conclusion was that, the centralized controller and the PI that sends information to all upstream reaches and to the closest downstream one, were the best options in the performance point of view. Another thing worth to mention is that they observed a possible control system destabilization when there exist a minimum gate movement restriction. Other similar PI scheme working together with a centralized controller was used in Montazar et al. (2005). van Overloop et al. (2005) also performed a decentralized PI tuning solving an optimization problem, but using, instead of one model for each flow condition, a set of models. The idea behind was the obtention of a more stable controller. Litrico et al. (2005) used each reach’s frequency response to tune PI controllers. Making use of the gain margin obtained for different discharge conditions, a series of robust controllers were calculated. In this manner, they achieved robustness against operation condition changes. The test performed in a laboratory canal showed also a great correspondence between the observed and the expected performances. A more detailed robustness analysis was performed in Litrico and Fromion (2006) and in Litrico et al. (2006). They proposed a new method to tune robust distant downstream proportional integral (PI) controllers for an irrigation canal pool. This tuning rules are appropriate to obtain specific robust margins and error characteristics. Implementation issues are also addressed. Finally, in Litrico et al. (2007) a classical closed loop PI tuning method, called ATV method, was adapted for irrigation canal decentralized level controllers. The method needs to induce sustained level oscillations to characterize the stability margins of the controlled system. In this paper, these parameters were linked with many tunings rules in order to obtain a good water level regulation performance in irrigation canals. 2.3. Types of control algorithms developed for irrigation canals 15 2.3.2 Robust control In addition to the robust PI controllers previously mentioned, other robust control techniques have also been employed. Litrico and Georges (1999b) suggested two irrigation canal controllers design methodologies: the a priori computation of a robust Smith Predictor and the trial and error tuning of a robust Pole Placement controller. Both of them used a nominal linear model and multiplicative uncertainty. These schemes were compared with the performance given by a PID tuned with the Haalman method, suitable for time delay dominant systems. The research concluded that a PID without a filter was faster, but also oscillating. In this context, the robust controllers fulfilled the established performance and robust requirements without major problems. Litrico (2001b) developed another methodology for robust controllers based on internal models, in this case, for the Gimone river in France. This river, in particular, has two water irrigation offtakes. They used the Hayami model and a multiplicative model uncertainty representation. The formulation was the following: they parameterized the filter value in order to obtain it, subsequently, assuring the closed loop robust stability. 2.3.3 Optimal control In spite of the use of this technique for tuning another types of controllers, there are few recent studies about it. In Malaterre (1998) they used a MIMO Linear Quadratic (LQ)-optimal control for irrigation canals. This control was developed together with a previously mentioned state observer. It could handle unexpected and beforehand scheduled demands. Additionally, the MIMO structure of the controller exhibited big advantages to counteract canal coupling and transport delay effects. Among the disadvantages of the method, they mentioned the large dimension of vectors and matrices, that the model validity is assured only for subcritical flows and the difficulty of LQoptimal control to include gates restrictions. In Reddy and Jacquot (1999) a proportional-plus-integral controller was developed for an irrigation canal with five pools using the linear optimal control theory. Different strategies were tested and it was found that the performance of regional constant-volume control algorithms was as good as the performance of a global control algorithm, whereas the performance of regional constant-level control algorithms was marginally acceptable. More recently Durdu (2005) used a Linear Quadratic Gaussian (LQG) control strategy for irrigation canals in order to test different state observers. 16 Chapter 2. Literature Review 2.3.4 Predictive control Similarly as occurs with Optimal control, there are few recent works that address the canal control with Predictive control techniques, and the ones that do, use in general classical techniques in this area. Malaterre and Rodellar (1997) performed a multivariable predictive control of a two reaches canal using a state space model. They observed that the increase of the prediction horizon produced a change in the controller behavior, varying the control perspective from a local to a global problem. Following the research line of Rodellar et al. (1989) and Rodellar et al. (1993), Gómez et al. (2002) presented a decentralized predictive irrigation canal control. They used the Muskingum model plus a storage model in order to perform the water dynamic predictions in each reach. In order to decouple the system, the controller used an estimation of the future discharges and the hypothesis of being linearly approaching the reference, to finally reach it, at the end of the prediction horizon. Because the control law solution was given in terms of reach’s inflow discharge, they used a local controller to adjust the gate opening to the required discharge. Akouz et al. (1998) used decentralized predictive controllers to manage three reaches of the ASCE Test Canal 2, acting on each reach’s inflow discharge. They didn’t include in the control, feedforward compensation for known scheduled demands or reach’s outflow discharges. The same technique was used in Ruiz and Ramírez (1998), including the reach’s outflow discharge as a known disturbance. In Sawadogo et al. (1998), and later in Sawadogo et al. (2000), they presented a similar decentralized adaptive predictive control, but that used the reach’s head gate opening as controllable variable and the reach’s tail gate opening and the irrigation offtake discharge as known disturbances. A decentralized adaptive predictive controller was also presented in Rivas et al. (2002). Here the manipulated variable was the inflow discharge and they did not include the known disturbances. In order to achieve some kind of robustness they used dead bands and normalization in the adaptation of the model. Sometimes it is convenient to take into account actuator and process constraints when controlling a particular system. In this respect, a constrained predictive control scheme was developed in Rodellar et al. (2003) to manage irrigation canals. It was based on a linear model that used gate openings and water levels as input and output variables respectively, and one of the novelties of the method was that it takes into account explicitly in the control problem that gates should not come out of water. Constraints on the movement velocity of the gates were also considered and the results exhibited an improvement in the control performance in comparison with the predictive unconstrained case. More recently, Wahlin (2004) tested a MIMO Constrained Predictive controller using a state 2.4. Conclusion 17 space model based on Schuurmans ID model. They performed tests where the controller either knew or did not know the canal parameters and with and without the minimum gate movement restriction. While many gate operation restrictions could be included in the control law, the minimum gate movement restriction could only be applied as a dead band in the control action once calculated. The reason for this is that this type of constrains are very difficult to implement in a controller. The results showed that it was possible to control the canal in question, but with a performance not superior than a centralized PI. Nevertheless, they conjectured that the problem was attributable to the modeling errors of the ID Model. In that case, a better model would be required in order to implement the predictive control. Additionally, they observed that the minimum gate movement restriction worsened, in a high degree, the control performance. There are also some real implementations of predictive control in laboratory canals. In Begovich et al. (2004), a multivariable predictive controller with constraints was implemented in real-time to regulate the downstream levels of a four-pool irrigation canal prototype. In Silva et al. (2007), a predictive controller, based on a linearization of the Saint Venant equations, has been also implemented on an experimental water canal. 2.3.5 Nonlinear Control Because of the complexity of the original nonlinear model, there are not many research works that had faced the nonlinear control for these types of systems. In Sanders and Katopodes (1999), they used a nonlinear optimization method for controlling one canal reach adjusting its gate openings. The computation times were near the three minutes with Pentium processors. Dulhoste et al. (2004) made a controller based on the dynamic state feedback linearization. The control was tested for set-point changes, infiltration and water extraction cases. There were good results on computational simulation for different length rectangular canals. In de Halleux et al. (2003), they described and analyzed a general stability condition for water velocities and levels in open channels. With the aid of it, they proposed and applied a controller to a two-reaches no-friction horizontal computer-simulated canal. In Soler et al. (2004), the nonlinear numerical scheme has been combined with a typical predictive control performance criterion to compute gate trajectories in an open-loop operation. 2.4 Conclusion In brief, this chapter reviews the conceptual/theoretical dimension and the methodological dimension of the literature in irrigation canal control and discovers research questions or hypotheses that are worth researching in later chapters. 24 Chapter 3. The Canal PAC-UPC 3.3 Water supply The 250 m3underground laboratory reservoir supplies the water to the canal. The water follows a path that is depicted in the aerial photo of figure 3.5. P TW MV MV MV MV CLR CLR: Const. Level Reservoir MV: Manual Valve P: Pumps TW: Triangular Weirs Pipes Water Figure 3.5: Water path from the pumps to the canal head The water can be delivered to the canal using any of the three available pumps (P). These pumps have pumping capacities of 100 L/s,200 L/sand 300 L/srespectively. The water is raised to a constant level reservoir (CLR) that is connected to two parallel canals with triangular weirs (TW) at the downstream ends through two electric valves. With these valves, it is possible to regulate the discharge that is delivered, whose value can be known very accurately by measuring the head water. Two manual limnimeters are located over the weirs to perform these level measurements. At this point, the water not entering the canals is returned to the underground reservoir and the measured discharge is conveyed through two pipes up to the canal head (shown in figure 3.6). The canal has also a small constant level reservoir at its head. The objective of this element is twofold: to dissipate the flow energy and to provide the canal with enough and virtually unlimited water. The canal takes water from this reservoir through Gate1, which is normally under submerged conditions. This gate can regulate the canal inflow by adjusting the gate opening. The water that is not used is also returned to the underground reservoir. Under normal operation conditions, the discharge delivered to the canal reservoir should be between 100 L/sand 150 L/s. Naturally, this value determines the maximum discharge that can be supplied to the canal, while the minimum discharge is a value close to zero given by the 3.3. Water supply 25 GATE 1 RETURN (a) Head reservoir (b) Gate1 Figure 3.6: Canal head 26 Chapter 3. The Canal PAC-UPC minimum available gate opening. 3.4 Gates (a) Gates (b) Control boxes Figure 3.7: Sluice gates in the Canal PAC-UPC The gates of the canal are vertical sluice gates (see figure 3.7) and were designed and constructed by the laboratory personnel. They are made of methacrylate reinforced with a metal skeleton in order to provide enough rigidity and a low weight. The vertical movement of the gates is guided by metal frameworks embedded in the canal and is executed by three-phase servomotors. These servomotors are located on top of the gates and are commanded by control boxes situated next to the canal. This particular gate motorization enables only constant speed movements of about 1 cm/s. 3.5. Weirs 27 The control boxes allow either the local or the distant operation of the gates. In local operation mode, it is possible to open, close or stop a gate, while in the distant operation mode it is possible to give an external reference signal in order to approximately position the gate at a desired opening. From these control boxes, it is also possible to have a real-time measurement of all the current gate openings. 3.5 Weirs The rectangular weirs of the canal are used to extract a water discharge and to measure its value in different sections of the canal, so as to emulate the effect of offtake discharges in real irrigation canals. Figure 3.8 shows some of the canal weirs in photos. (a) Opened weir (b) Closed weir Figure 3.8: Rectangular weirs in the Canal PAC-UPC These weirs have a width of approximately 39 cm and were constructed starting from a 35 cm canal height, except for the canal end weir (Weir4) that starts from the canal base. From this minimum height, it is possible to raise the height of any weir by placing measured-height PVC pieces in specially designed metal rails. Each weir has its own set of pieces. There are pieces of 5 cm,10 cm,20 cm and 35 cm. With different combinations of these pieces, it is possible to cover a broad range of weir heights. Weir 4 plays also an important role in the overall canal operation. By raising and lowering this weir height, it is possible to change the minimum water level value throughout the whole 28 Chapter 3. The Canal PAC-UPC canal. Here, the head value acts as a reference value for the other pools, a phenomena that allows to work with different combinations of discharges and water level values, leading to different hydraulic conditions. 3.6 Sensors 3.6.1 Water level sensors The canal has currently nine level sensors located in strategic places. Specifically, these sensors are situated upstream and downstream of each gate and at every rectangular weir (see figure 3.3 for details). Their mission is twofold: to take level measurements where water is taken from the canal and to enable the calculation of gate and weir discharges by means of appropriate hydraulic relationships. These sensors are submergible pressure sensors and were acquired from the specialized manufacturer SOFREL. Their specific model is CNPI and they can measure water depths up to 2 m. One of these sensors is shown in figure 3.9. Figure 3.9: Type of water level sensor used in the Canal PAC-UPC These sensors are mainly oriented for level measurements in real canals. Consequently, this work had to deal with practical problems similar to the ones encountered in the field. The installation of the sensors inside the canal was carried out in a the following way: the bodies of the sensors were firmly attached to the internal side of the lateral walls at 3 cm from the floor and their cables were routed to the control room separately from any power cable to avoid electrical noise. Following appropriate wiring and signal filtering procedures, they are able to provide precise level measurements. The way they are connected is sketched in figure 3.10. 3.6. Sensors 29 Level Sensor I V Signal Transducer Signal Filter Voltage Source fc=10 Hz 24V Cable Shield 4-20mA WATER LEVEL MEASUREMENT SIGNAL ≈ 15 m Figure 3.10: Connection diagram of a level sensor Once installed in the canal, these sensors where calibrated against manual limnimeter measurements. The calibration procedure was as follows: the weirs were closed and the canal was filled with a certain amount of water; one hour after, level measurements were carried out in calm water and the process was repeated. As a result, several data points were collected and an experimental calibration curve was computed for each sensor using linear regression analysis. The calibration curve of one of the level sensors is presented in figure 3.11 as an example. 3000 3200 3400 3600 3800 10 20 30 40 50 60 70 80 90 100 Level sensor signal (DAQ card raw) Water level with manual limnimeter (cm) y = 0.07428*x − 198.61 data point linear Figure 3.11: Calibration curve of level sensor 11 These level sensors are based on pressure measurements. These type of measurements are 30 Chapter 3. The Canal PAC-UPC robust against turbulent water, but can be biased by effects related to the water temperature. Hence, it is recommendable to perform periodic sensor recalibrations, specially between the seasons of the year. 3.6.2 Gate position sensors Each servomotor has a position measuring gear embedded in the motor chassis. This gear rotates when the gate is opened or closed and by means of this event the gate opening is sent to the control box. What actually happens is that the gear is mechanically connected to a potentiometer, thereby giving different resistance values for different gate openings. This signal is successively transformed up to the control room computer. The flow path of this gate position measurement signal is summarized in figure 3.12. Servo Control Box Potentiometer Position Gear (Servomotor) ΩmA Resistor GATE POSITION SIGNAL V (≈ 5m) (≈ 10m) Figure 3.12: Flow path of the gate position signal If the calibration step is performed with data arising from opening and closing a gate, it is possible to see that the data points group along two different calibration curves. This peculiarity is highlighted in figure 3.13 for the calibration data of Gate 1. 2400 2500 2600 2700 2800 2900 3000 3100 0 50 100 150 200 250 300 Gate position sensor signal (DAQ card raw) Gate opening with ruler (mm) linear linear data point data point Figure 3.13: Gate 1 position calibration data One of these data subsets correspond to move the gate upwards and the other to move it 3.7. Control room 31 downwards. This is due to the inaccuracy of the position measurement method: any mechanical system based on gears present gaps between the teeth. As a result, when the rotation direction is inverted in this case, there is an interval of time during which the potentiometer do not move, generating two different electrical readings for the same position depending on the gate moving direction. This is a problem for calibration. One possible solution is to pass through a unique calibration line equidistantly from both data subsets. This would produce a permanent measurement error. In this case, this error would be approximately ±1 cm. In order to improve the measurement accuracy another option was taken: a correction algorithm was implemented in the SCADA software by using two calibration curves and a gate moving direction detection method. This approach does not entirely solve the problem, but most of the time works fine for the required precision. Nevertheless, gate opening measurements were performed always manually when calibrating gate discharge equations. 3.6.3 Flowmeter For comparison and calibration purposes, the canal is equipped with a movable flow velocity sensor. This flowmeter has been manufactured by NIVUS and consist of an OCM EM control box and a doppler-based wedge-mouse one-dimensional flow velocity sensor. This device can measure velocities up to 6 m/sand working in combination with a level sensor is able to measure discharge. This device needs a strictly parallel flow in order to give good measurement. Consequently, it should be placed at a considerable distance from any hydraulic structure. In the Canal PAC-UPC, this device has been placed downstream of Gate 1 in order to have an independent measurement of the canal inflow. The preferred way to measure gate discharges throughout this thesis will be to calibrate and use gate discharge hydraulic equations, a method that has been proven to be highly accurate and robust in our laboratory tests. Essentially, the flowmeter measurement provides a way to continuously validate these discharge calculations, but it is not used in the control system of the canal. 3.7 Control room The Canal PAC-UPC has a control room over the canal. All sensors are wired to it and from there it is also possible to operate the servomotor control boxes. The brain of the control room is a 3.06 GHz Pentiumr4 processor, 1 GiB memory desktop computer equipped with three AdvantechrData AcQuisition (DAQ) cards. Illustrative photos are presented in figure 3.14. 32 Chapter 3. The Canal PAC-UPC (a) Control room (b) Control computer Figure 3.14: Canal PAC-UPC control installations 3.8. Software 33 3.7.1 Advantech DAQ cards The models of the data acquisition cards inside the computer are PCI-1711, PCI-1720 and PCI1750. All three cards are situated in PCI-slots and perform different tasks in the canal operation. The PCI-1711 is a multifunction DAQ card that has 16 analog inputs and 2 analog outputs. In this case, these inputs are used to receive the 9 level sensor measurements, the 3 gate position signals and 1 flow velocity signal. On the other hand, the analog outputs are used to give position reference values to Gate1 and Gate3. The PCI-1720 completes the set of required analog outputs. This card dispose of 4 analog outputs, but connection incompatibilities between the control boxes and this card allow the use of only one of these outputs. Thus, this card takes over the operation of Gate5. The PCI-1750 has not been used in this work, but can manage external alarm signals and receive orders from independent control button panels. All these cards can achieve very fast sampling rates, but water measurements in irrigation canals do not need a so continuous update. Hence, they are operated at a much slower rate. In particular, level measurements and gate openings are acquired at 10 Hz and gate opening set-points are delivered every 10 s. The resolution of the measurements deserve special attention. This resolution depends on the sensor resolution and on the DAQ card resolution. The PCI-1711 has a 12 bit resolution. That means that the card is able to distinguish 4096 different values when using the complete input range. However, this is seldom the case. In this case, the sensor signal ranges are considerably smaller, yielding water level measurements and gate openings with a resolution slightly smaller than 1 mm. 3.8 Software The data acquisition cards provide the hardware to receive and send electrical signals from the computer. However, a specialized computer program is completely indispensable to manage this huge amount of information and to provide a basis where to perform discharge computations and to implement control algorithms. 3.8.1 Base software In this case, it was decided to use the software package MATLABr. In particular, real-time programs can be easily implemented by using three of its components: Simulinkr, Real-Time Workshoprand Real-Time Windows Target. Each one of them plays a particular role. Simulink is a platform for multidomain simulation and model-based design for dynamic systems. It provides an interactive graphical environment and a customizable set of block libraries. 40 Chapter 3. The Canal PAC-UPC 3.9 Conclusion This chapter has described in detail the Canal PAC-UPC. Its purpose, operation and several of their elements has been presented along with technical details and operational constrains. An important practical knowledge related to the instrumentation of canals has been acquired, touching on topics like treatment of measurement errors, signal processing techniques, calibration of sensors, etc.. All this know-how is indispensable when working with real canals. This chapter has also presented one of the product of this thesis: an own, non-commercial, SCADA software for the Canal PAC-UPC. Many of the elements, solutions and procedures used in its development are applicable to similar systems in real canals. Based on the overall chapter content, it is possible to confirm that this laboratory canal provides a good platform where to test irrigation canal automation issues. Chapter 4 Calibration of Weirs and Sluice Gates This chapter deals with some issues related to discharge measurements in canals; particulary in our laboratory canal: the Canal PAC-UPC. It describes mathematical and empirical formulations to calculate discharge, the calibration of our weirs and sluice gates, and some results obtained from our observations. 4.1 Introduction Open channel flow is defined as the flow in any channel where the water flows with a free surface. Open channel flow is not under pressure; gravity is the main force that can produce flow in open channels, and a progressive declive in water surface elevation always occurs as the flow moves downstream. Examples of open channel flow include: rivers, streams, creeks, discharges from tailings ponds, and other uncovered conduits. Closed channels, such as adits, tunnels, sewers, and ventilation shafts, can be treated as open channels when flowing partially full and not under pressure. Measuring discharge in open channels can be a difficult task. There are two types of approaches: 1. Measure directly all the involved variables (flow area and mean velocity). 2. Define a so called control section, where there is a single relationship between Q/h, so by measuring water depths we can derive the Qvalue. Example of control sections are weirs, parshall flumes, flumes and gates. This chapter focus on discharge characteristics of two types of hydraulic structures that belong to the second group: weirs and gates. 42 Chapter 4. Calibration of Weirs and Sluice Gates 4.2 Weirs Weirs are typically installed in open channels to determine discharge (flow rate). Since the geometry of the weir is known and all water flows through it, the depth of water flowing over the weir (head) is an indication of the discharge value. Two different geometries are considered below. 4.2.1 Sharp-edged v-notch (triangular) weirs The discharge is directly related to the water depth above the crotch (bottom) of the V; this distance is called head (h). The V-notch design assumes that small changes in discharge to have a large change in depth allowing more accurate head measurement than with a rectangular weir. B O h End View Water Flow θ Figure 4.1: V-notch weir scheme The 90◦V-notch weir is most accurate when measuring discharges of less than 0.03 m3/s. The maximum discharge that can be accurately measured is approximately 0.3 m3/s. The sides of the notch are inclined outwardly at 45◦from the vertical. The V-notch weir equations have become somewhat standardized. The water measurement manuals of the International Organization for Standardization (ISO), the American Society for Testing and Materials (ASTM) and the U.S. Department of the Interior, Bureau of Reclamation (USBR) all suggest using the Kindsvater-Shen equation, which is presented below: Q=8 15p2g Cvw tan θ 2(h+corrh)5/2(4.1) where Cvw is the v-notch weir discharge coefficient, his the head (m) and corrhis the head correction factor (m). Particulary for the 90◦V-notch: Cvw = 0.578 corrh= 8.847 ×10−4m 4.2. Weirs 43 In this way, we used only (4.1) to calculate the discharges given by our two v-notch weirs. 4.2.2 Sharp-edged rectangular weirs The rectangular weir is the most commonly used thin plate weir. As its name suggests, it has a rectangular opening of a certain width (b). It is more suitable for larger flows, because the width can be chosen so that it can pass the expected flow at a suitable depth. B b O h End View Water Flow Figure 4.2: Rectangular weir scheme Depending on its width the weir can be "suppressed", "partially contracted" or "fully contracted". Suppressed means there are no contractions. A suppressed weir’s notch width (b) is equal to the channel width (B); thus, there really is not a notch; the weir is flat all the way along the top. For a weir to be fully contracted, (B−b) must be greater than 4hmax, where hmax is the maximum expected head on the weir. A partially contracted weir has B−bbetween 0 and 4hmax. Weir contractions produce the water flow lines to converge through the notch. To provide a single reliable, accurate method to model all rectangular weirs (suppressed, partially contracted, and fully contracted), the ISO, the ASTM and the USBR, all recommend using the Kindsvater-Carter method for all rectangular weirs. The equation is the following: Q=2 3p2g Crw(b+corrb)(h+corrh)3/2(4.2) where Crw is the rectangular weir discharge coefficient, his the head (m), the sum b+corrbis called "effective width" (m) and the sum h+corrhis called "effective head" (m). The value of corrhis normally taken as 0.001 m. The discharge coefficient Crw is a function of b/B and h/O, and corrbis a function of b/B. In general these relationships are of linear nature and are given, in technical manuals, by graphs and equations. Particulary for our canal, we obtained these relationships empirically. We measured head values, using a limnimeter, and discharges with our v-notch weirs. 44 Chapter 4. Calibration of Weirs and Sluice Gates Putting all constant values together, and taking into account the fact that b/B is fixed for a particular weir, (4.2) yields: Q=Cgen h O(h+ 0.001)3/2(4.3) From (4.3) one can see that the discharge can be calculated by measuring h. However, the mathematical expression of function Cgen h Ois not exactly known and has to be determined for each particular weir. This task can be performed by taking head-discharge data pairs using other measuring devices and reordering (4.3) as Q (h+0.001)3/2=Cgen h O. Using this approach, Cgen h Ocan be obtained solving a linear regression problem. It should be noted that with this approach, the weir width bis included in the general coefficient Cgen and does not need to be measured. We have three rectangular weirs in operation: Weir 1, Weir 3 and Weir 4. Although not needed in this case, their widths are 0.434 m,0.44 m and 0.44 m respectively. All measurements were performed after waiting a two-hours stabilization period. The calibration data and the obtained relationships are detailed in table B.1, table B.2 and table B.3, and in figure 4.3, figure 4.4 and figure 4.5. y = 0.184215x + 0.576043 R2 = 0.956383 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 Head / Weir height, h/O (m/m) G. Coefficient, Cgen (m1.5/s) Figure 4.3: Weir 1 Calibration curve 4.2. Weirs 45 y = 0.174193x + 0.590874 R2 = 0.971928 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 Head / Weir height, h/O (m/m) G. Coefficient, Cgen (m1.5/s) Figure 4.4: Weir 3 Calibration curve y = 0.105734x + 0.653631 R2 = 0.988943 0.65 0.66 0.67 0.68 0.69 0.7 0.71 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 Head / Weir height, h/O (m/m) G. Coefficient, Cgen (m1.5/s) Figure 4.5: Weir 4 Calibration curve Now, with the information given by these graphs and (4.3), the weir discharges Qcan be accurately calculated measuring the head h. As it can be observed in figures 4.3, 4.4 and 4.5, all linear relationships have a high correlation number. That means that the relationships suggested by technical manuals are valid for our rectangular weirs too. A more detailed comparison between the curves and equations given in manuals and our empirical equations, reveals that they differ slightly. However, it is worth to remember that the curves and equations given in technical manuals must follow some strict installation guidelines to ensure their applicability. In this case, it was not always possible to follow exactly these guidelines: the main reason for the discrepancies in our opinion. 46 Chapter 4. Calibration of Weirs and Sluice Gates 4.3 Sluice gates A gate is a hydraulic structure widely used for controlling discharge and water depth in irrigation and drainage canals. However, it can also be used as a convenient discharge measuring device. Depending on the intended purpose and the particular mechanical design, there are different types of gates, many of them widely used in irrigation canal applications: vertical lift gates, radial (tainter) gates, roller gates, flap gates, overshoot gates and so forth. A sluice gate is traditionally a wooden or metal plate, vertical (vertical lift gates) or curve (radial (tainter) gates), which slides in grooves in the sides of the channel. Sluice gates are commonly used to control water levels and flow rates in rivers and canals. Raising a sluice gate allows water to flow under it. The term sluice gate refers to any gate that operates by allowing water to flow under it. When a sluice gate is fully lowered, water sometimes spills over the top, in which case the gate operates as a weir. Usually a mechanism drives the sluice gate up and down. This may be a simple, handoperated, worm drive or rack and pinion drive, or it may be electrically or hydraulically powered. This section is focused on vertical lift gates mainly because the laboratory canal has only this type of gates. However, the majority of the following developments can also be applied to other types of sluice gates. 4.3.1 Flow conditions It is very important to remark the existence of two particular working flow conditions when dealing with sluice gates: 1. The free flow condition 2. The submerged flow condition. Both conditions are sketched in figure 4.6 and figure 4.7 respectively. Under a free flow condition (see figure 4.6), a hydraulic jump occurs downstream from the sluice gate in a channel. The downstream conjugate depth of the jump h3may be calculated by taking the water depth at vena contracta h2as the upstream conjugate depth. A submerged flow occurs when the tailwater depth is greater than the downstream conjugate depth of h2. As illustrated in figure 4.7, instead of being in presence of a normal hydraulic jump, this particular condition develops a submerged hydraulic jump. The contact of this hydraulic jump with the volume of water above it induces a turbulent water recirculation. This phenomena produces typically a backwater flow in the surface layer. 4.3. Sluice gates 47 h1h3 h 2 gate hydraulic jump free flow l Figure 4.6: Sketch of free flow h1h3 h 2 gate subm. hydraulic jump submerged flow l Figure 4.7: Sketch of submerged flow To determine whether the jump will be free or submerged is another problem. In this respect, there are several formulas that have been obtained theoretically or empirically by many researchers. Moreover, there are some formulas that do not consider an absolute frontier, but a gradually changing transition. Because this topic is out of the scope of this work, the following simple condition presented in Swamee (1992) was used when needed to distinguish both working flow conditions: Free flow: h1≥0.81h3h3 l0.72 (4.4) Submerged flow: h3< h1<0.81h3h3 l0.72 (4.5) 48 Chapter 4. Calibration of Weirs and Sluice Gates 4.3.2 Flow rate formulas 4.3.2.1 Classical theoretical formulation Calculating the discharge under a sluice gate is not a trivial matter. In theory, the sluice flow rate formula can be accurately obtained if the contraction coefficient is known (Henderson, 1966). The most common flow rate expression makes use of the conservation of energy, mass and momentum in the sluice gate - hydraulic jump flow. This procedure yields the following equation: Q=Cdl b p2g h1(4.6) In (4.6), lis the gate opening, bis the gate width and h1is the upstream water level. In this context, the discharge coefficient Cdis given by two equations (one for each flow condition), functions of the contraction coefficient Cc,b,h1and, for the submerged condition, the downstream water level h3: Free flow: Cd=Cc √1 + η(4.7) Submerged flow: Cd=Cc  ξ−sξ2−1 η2−121−1 λ2  1/2 1 η−η (4.8) where η=Ccb/h1,ξ= (1/η −1)2+ 2(λ−1) and λ=h1/h3. Unfortunately, the contraction coefficient varies with: the amount of gate opening, shape of the gate lip, upstream water depth, gate type and so forth (Lin et al., 2002). Thus, it is very difficult to know its true value for all operating conditions in practice. That is why there are other approaches that combine some theoretical and some practical knowledge in order to simplify the task. Some of them are presented in the next sections. 4.3.2.2 Classical empirical formulations Free Flow A good survey that includes many of these approaches for the free flow condition can be found in Montes (1997, 1999); Webby (1999) and in Speerli and Hager (1999). One common approach is to use (4.6) and determine an empirical constant value or, moreover, a relationship for Cdor Cc. Writing this in a mathematical form leads: Q=Cdl b p2g h1 with Cd=f(h1, l)or Cd=k0. 4.3. Sluice gates 49 Submerged Flow As mentioned in Clemmens et al. (2003), few studies are available in the literature. There are two common approaches: use (4.6) as well, including the effect of the downstream water level h3in the discharge coefficient Cdor modify (4.6) in order to incorporate h3explicitly (Malaterre and Baume, 1998). That is: Q=Cdl b p2g h1 with Cd=f(h1, h3, l)or Cd=k0, or Q=Cdl b p2g(h1−h3)(4.9) with Cd=f(h1, l)or Cd=k0. It is worth to note that a particular sluice gate in a normal irrigation canal operates, the most of the time, either in the free flow condition or in the submerged flow condition. Moreover, they are usually operated in a very narrow range, which explains, in some cases, the use of fixed discharge coefficient values. However this is not always advisable, because the discharge coefficient can suffer abrupt changes around a particular working point. In this context, it is very clarifying the experimental research performed by Henry (1950). Figure 4.8 shows the variation of Cd(using (4.6)) under free and submerged flow as obtained by him. 0.6 0 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 5 6 7 8 10 12 14 16 h3 l=2 3 4 5 6 7 8 FREE FLOW SUBMERGED FLOW Discharge coefficient, Cd Upstream level / Gate opening, (m/m) h1 l Figure 4.8: Variation of Discharge coefficient From figure 4.8 one can observe that Cdcan be very sensitive to small changes in any of the involved variables. For free flow, Cdprogressively increases to a constant value of 0.611. Under submerged flow conditions, Cdis zero when h1=h3. Any increase in h1above h3results in 56 Chapter 4. Calibration of Weirs and Sluice Gates Other alternatives are Weighted Least Squares (WLS), MM-estimators, and so forth. The calibration results are summarized in table 4.1. Table 4.1: Calibration results for Gate 1, Gate 3 and Gate 5 SIC HEC-RAS Ferro Gate Id. µ0CdC0 dk0 0k0 1 Gate 1 0.3356 0.438 0.609 0.9176 0.3489 Gate 3 0.3347 0.388 0.678 0.9482 0.3202 Gate 5 0.3790 0.405 0.768 1.0097 0.3154 At this point, all methods were evaluated in terms of predicted discharge errors, using the same data sets, in order to analyze and compare the model-fitting properties of each approach. 4.3.3.3 Results These results are presented in two forms. First of all, discharge predictions and their residuals are shown graphically, in order to observe tendencies and assess performance in a more qualitative way. Next, different mean error indices are calculated and contrasted, in order to obtain numerical values that permit a quantitative evaluation of each discharge calculation method. The predicted discharges, calculated from the application of each method to the experimental data from Table C.1, Table C.2 and Table C.3, are shown for Gate 1, Gate 3 and Gate 5 in figure 4.9, figure 4.10 and figure 4.11 respectively. The discharge errors of these predictions are shown, separately, in figure 4.12, figure 4.13 and figure 4.14. 0 5 10 15 20 25 30 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Sample number Predicted Discharge, Q (m3/s) SIC Swamme Ferro Raj&Sub HEC−RAS Cc=0.611 Henry Measured Figure 4.9: Gate 1 predicted discharges with ±10 % error bars 4.3. Sluice gates 57 0 2 4 6 8 10 12 14 16 18 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Sample number Predicted Discharge, Q (m3/s) SIC Swamme Ferro Raj&Sub HEC−RAS Cc=0.611 Henry Measure Figure 4.10: Gate 3 predicted discharges with ±10 % error bars 0 2 4 6 8 10 12 14 16 18 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Sample number Predicted Discharge, Q (m3/s) SIC Swamme Ferro Raj&Sub HEC−RAS Cc=0.611 Henry Measure Figure 4.11: Gate 5 predicted discharges with ±10 % error bars 58 Chapter 4. Calibration of Weirs and Sluice Gates 0 5 10 15 20 25 30 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 Sample number Discharge error, eQ (m3/s) no error SIC Swamee Ferro Raj&Sub HEC−RAS Cc=0.611 Henry Figure 4.12: Gate 1 discharge errors 0 2 4 6 8 10 12 14 16 18 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 Sample number Discharge error, eQ (m3/s) no error SIC Swamee Ferro Raj&Sub HEC−RAS Cc=0.611 Henry Figure 4.13: Gate 3 discharge errors 4.3. Sluice gates 59 0 2 4 6 8 10 12 14 16 18 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 Sample number Discharge error, eQ (m3/s) no error SIC Swamee Ferro Raj&Sub HEC−RAS Cc=0.611 Henry Figure 4.14: Gate 5 discharge errors The first feature that shows up from observing the figures, is that, in the operational range of this study, the Swamee method gave a very poor performance, exhibiting high discharge errors in all the studied gates. In fact, the direct inspection of Henry’s nomogram (the Swamee method only gives equations for Henry’s curves) gave more accurate results in the majority of the cases, when it was possible its application (there were several points that went out of the range of the nomogram (1< h1/l ≤16)). Without taking into account the Swamee method, the rest of the methods fell, in general, within the ±10 % error margin. The Ferro method performed very well in all three gates, displaying a very good fit (this method, the SIC and the HEC-RAS method were calibrated using the same data as in this evaluation) with very small errors in most of the cases. The HEC-RAS fit, although oscillating around zero, seems to be more erratic. The SIC method has a similar behaviour. Surprisingly, the use of a constant contraction coefficient value of 0.611 in the theoretical formulas gave very good results. The discharge predictions were very accurate with low dispersion. The results given by the Rajaratnam method and the Henry method look also good. In most of the cases, they predicted discharges with small error in a reliable manner. In order to evaluate quantitatively these results, the error based performance indices of table 4.2 were used. The Mean Error (ME), the Mean Absolute Error (MAE) and the Root Mean Square Error (RMSE) are indices measured in the same units as the error. They can be influenced in some way, if the errors depend on the magnitude of the measurements. Conversely, the Mean Percentage 60 Chapter 4. Calibration of Weirs and Sluice Gates Table 4.2: Performance indices Index Formula MAPE 100 N N X i=1  Yi−ˆ Yi Yi MPE 100 N N X i=1 Yi−ˆ Yi Yi ME 1 N N X i=1 Yi−ˆ Yi MAE 1 N N X i=1 Yi−ˆ Yi RMSE v u u t1 N N X i=1 Yi−ˆ Yi2 Error (MPE) and the Mean Absolute Percentage Error (MAPE) are relative indices given in percentage. The ME and the MPE are signed measures of error which indicate whether the predictions are biased, i.e. whether they tend to be disproportionately positive or negative. Bias is normally considered a bad thing, but it is not the bottom line. On the other hand, the MAE, the MAPE and the RMSE are indices that incorporate the bias as well as the variance of the errors. The RMSE takes normally precedence over the others because has a close relationship with the standard deviation, but it is more sensitive than the other measures to the occasional large error. The results for each gate are summarized in table 4.3, table 4.4 and table 4.5. In order to simplify comparison purposes, these results are also given in bar charts. They are shown in figures 4.15 - 4.20. Table 4.3: Performance of each method for Gate 1 data Method MAPE MPE ME MAE RMSE (%) (%) (m3/s) (m3/s) (m3/s) SIC 5.11 -1.47 -0.0002 0.0017 0.0024 Swamee 25.73 25.73 0.0093 0.0093 0.0101 Ferro 1.61 -0.15 0.0001 0.0005 0.0008 Raj & Sub 5.05 3.36 0.0014 0.0018 0.0021 HEC-RAS 9.38 -1.67 0.0000 0.0034 0.0043 Cc=0.611 3.14 -0.29 0.0000 0.0010 0.0014 Henry 4.08 -0.27 0.0003 0.0014 0.0017 4.3. Sluice gates 61 Table 4.4: Performance of each method for Gate 3 data Method MAPE MPE ME MAE RMSE (%) (%) (m3/s) (m3/s) (m3/s) SIC 7.00 -4.03 -0.0006 0.0023 0.0025 Swamee 21.13 21.13 0.0058 0.0058 0.0064 Ferro 2.09 1.41 0.0008 0.0010 0.0016 Raj & Sub 1.80 1.49 0.0004 0.0005 0.0006 HEC-RAS 4.04 -2.57 -0.0004 0.0013 0.0013 Cc=0.611 2.56 -2.35 -0.0010 0.0010 0.0014 Henry 6.48 4.98 0.0020 0.0027 0.0033 Table 4.5: Performance of each method for Gate 5 data Method MAPE MPE ME MAE RMSE (%) (%) (m3/s) (m3/s) (m3/s) SIC 10.21 -5.48 -0.0008 0.0035 0.0042 Swamee 23.88 23.88 0.0067 0.0067 0.0070 Ferro 2.42 1.73 0.0010 0.0012 0.0024 Raj & Sub 9.61 9.61 0.0033 0.0033 0.0038 HEC-RAS 5.21 -2.68 -0.0004 0.0018 0.0024 Cc=0.611 5.95 5.95 0.0020 0.0020 0.0023 Henry 10.57 10.57 0.0036 0.0036 0.0038 SIC SWAM FERR R&S HRAS CC06 HENR −2 0 2 4 6 8 10 12 x 10−3 (m3/s) ME MAE RMSE Figure 4.15: ME, MAE and RMSE for Gate 1 62 Chapter 4. Calibration of Weirs and Sluice Gates SIC SWAM FERR R&S HRAS CC06 HENR −2 0 2 4 6 8 10 12 x 10−3 (m3/s) ME MAE RMSE Figure 4.16: ME, MAE and RMSE for Gate 3 SIC SWAM FERR R&S HRAS CC06 HENR −2 0 2 4 6 8 10 12 x 10−3 (m3/s) ME MAE RMSE Figure 4.17: ME, MAE and RMSE for Gate 5 4.3. Sluice gates 63 SIC SWAM FERR R&S HRAS CC06 HENR −5 0 5 10 15 20 25 30 (%) MPE MAPE Figure 4.18: MPE and MAPE for Gate 1 SIC SWAM FERR R&S HRAS CC06 HENR −5 0 5 10 15 20 25 30 (%) MPE MAPE Figure 4.19: MPE and MAPE for Gate 3 64 Chapter 4. Calibration of Weirs and Sluice Gates SIC SWAM FERR R&S HRAS CC06 HENR −10 −5 0 5 10 15 20 25 (%) MPE MAPE Figure 4.20: MPE and MAPE for Gate 5 The performance results are in good agreement with the qualitative results. The Swamee method presented a very high bias error for all three gates. In all the tests, it underestimated the gate discharges in more than 20 %. The rest of the methods can be joined into two groups under a performance point of view. The first group is integrated with methods whose results oscillate in the 2 %-10 % MAPE range, depending on the particular studied gate. This group is integrated with the SIC method, the HEC-RAS method, the Rajaratnam method and the Henry method. The second group showed a high performance consistency through gates. For example, the Cc=0.611 method had the following MAPE values: 3.14 %,2.56 % and 5.95 %. The other method in this group had the best overall performance results: the Ferro method. This fact is exemplified by its low MAPE values: 1.61 %,2.09 % and 2.42 %. The last performance test involves the evaluation of the SIC method, the HEC-RAS method and the Ferro method using typical calibration values. This values were obtained from manuals and research works, and are summarized in table 4.6. Table 4.6: Typical calibration values Method Calibration values SIC µ0= 0.4 Ferro k0 0= 1.0559 k0 1= 0.3344 HEC-RAS Cd= 0.6C0 d= 0.8 Because there was a clear pattern among all gates, these results are only presented for Gate 3 using the MPE and the MAPE for simplicity reasons. These results are shown in figure 4.21. 4.3. Sluice gates 65 SIC SWAM FERR R&S HRAS CC06 HENR −30 −25 −20 −15 −10 −5 0 5 10 15 20 25 30 (%) MPE MAPE Figure 4.21: Performance degradation without calibration - Gate 3 This test showed a high increment of the errors (essentially bias errors) when using typical calibration values (prior errors were colored in black). As a matter of fact, the performance degradation was so pronounced (MAPEs ≈20 %) that, in general, the rest of methods performed far better than them. Additionally, the performance results from the untuned Ferro method were always better than the ones obtained using the untuned, HEC-RAS and SIC, methods. 4.3.3.4 Discussion One of the most impressing issues from the results section was the poor performance exhibited in these tests by the Swamee method. However, this point can be made clear after comparing Henry’s original nomogram with Swamee’s approximation. This comparison is sketched in figure 4.22. As it can be seen in figure 4.22, Swamee’s discharge equations produce only a good fit for h3/l < 5. For h3/l ≥5, there are clear differences between both diagrams. Besides, these differences tend to increase with the quotient h3/l. Therefore, the poor results obtained in this case were only an example of these fitting errors: most of the data lied in that range. All the rest of the methods achieve a more than acceptable performance. However, there are some differences that are worth to remark. A difference should be made when comparing the results of methods that were calibrated using field information with the ones that were not. It is clear that these methods should perform much better than the others, because, in this case, they were calibrated using the same data that were, subsequently, used for model validation. However, with the particular exception of the Ferro method, the SIC method and the HEC-RAS method performed, in general, equally than 72 Chapter 5. Canal Identification for Control Purposes Storage Area Transport Area Gate iGate i+1 QiQi+1 QL i Zs i Qi 0 X X Figure 5.2: Simplified representation of a reach •Transport ∂A ∂t +∂Q ∂x = 0 (5.1) ∂Q ∂t +∂ ∂x Q2 A+gA∂Z ∂x =gA(S0−Sf)(5.2) with initial conditions: Z(x, 0) = Z0(x),Q(x, 0) = Q0(x)and boundary conditions: Z(0, t) = Zi(t),Q(0, t) = Qi(t),Z(X, t) = Zi(t),Q(X, t) = Qi(t). •Storage Qi(t)−Qi+1(t)−QL i(t) = dVs i(t) dt (5.3) In the transport equations (5.1)-(5.2), xis the longitudinal coordinate in the flow direction, t is the time, A=A(x, t)is the wetted cross-sectional area of the pool, Q=Q(x, t)is the water discharge, Z=Z(x, t)is the water depth, S0is the bottom slope, Sfis the friction slope of the canal and gis the gravity acceleration. Besides, it should be noted that the wetted cross-sectional area (A) depends explicitly on the water depth (Z) and on the shape of the cross-section of the canal. In the storage equation (5.3), Qi(t)is the water discharge that enters the zone, Qi+1(t)is the water discharge delivered to the next reach, QL i(t)is the water discharge extracted for irrigation purposes and Vs i(t)is the water volume (that is a function of the water depth Zs i(t)and the geometry of that zone) stored behind gate i+ 1. 5.2. Model of a pool 73 As can be seen from figure 5.2 and equations (5.1), (5.2) and (5.3), the values of the variables in the interface, denoted by an upper line Zi(t), Qi(t), provide the link between both areas. Normally, they are not known a priori, but making use of extra mass and energy conservation relationships among that particular point and the described zones, the problem can be solved. In summary, the model consists in a system of two nonlinear Partial Differential Equations (PDEs) and one nonlinear (sometimes linear) Ordinary Differential Equation (ODE). Because of the reasons given before, the model is little advantageous for control purposes. In the search for a more convenient model, a linearization around an operational condition will be performed. 5.2.2 Linearization of the model Linearization is carried out replacing in the model described by (5.1), (5.2) and (5.3), the expressions of the variables around a working point, namely, Z(x, t) = Z0(x) + z(x, t)and Q(x, t) = Q0+q(x, t), and neglecting all the second-order terms (Litrico and Fromion, 2004c). For the Saint-Venant equations (5.1) and (5.2), this yields: B0 ∂z ∂t +∂q ∂x = 0 (5.4) ∂q ∂t + 2V0 ∂q ∂x −β0q+C02−V02B0 ∂z ∂x −γ0z= 0 (5.5) with γ0=V02dB0 dx +gB0h(1 + κ)S0−1 + κ−F02(κ−2)∂Z0 ∂x i,β0=−2g V0S0−∂Z0 ∂x  and κ=7 3−4S0 3B0P0 ∂P0 ∂Z , in case Manning equation is used to represent the friction slope Sf. In (5.4) and (5.5), F0=V0 C0is the Froude number, P0is the wetted perimeter, C0=qgA0 B0 is the water wave celerity and V0=Q0 A0is the water velocity; all of them evaluated at the operational condition. Using the linearized form, the boundary conditions of the model are now: q(0, t) = qi(t),q(X, t) = qi(t), and z(0, t) = zi(t),z(X, t) = zi(t). One way to obtain a solution for this system is applying the Laplace transform and then reordering. This produces the following system of ordinary differential equations in the variable x, with a complex parameter s(the Laplace variable): ∂ ∂x "q(x, s) z(x, s)#=A(x, s)"q(x, s) z(x, s)#1(5.6) with A(x, s) =      0−B0(x)s −s+β0(x) B0(x)C0(x)2−V0(x)2 2V0(x)B0(x)s+γ0(x) B0(x)C0(x)2−V0(x)2      . 1In the following any f(s)will correspond to L {f(t)}, the laplace transform of f(t), which is a complex-valued function. 74 Chapter 5. Canal Identification for Control Purposes Because matrix Adepends on the variable x, there is not a closed solution to the differential equation, and therefore, it is necessary to use a numerical integration method to obtain the solution. Only the case where Ais not dependent on xhas an analytical solution. This special case is called uniform regime and is characterized by having the same water depth and the same water flow throughout a canal. It has been proven in Litrico and Fromion (2004c) that the problem can be solved numerically very efficiently, if it can be discretized by several "mini uniform regimes problems" in a way as illustrated in figure 5.3. Z x X1X2X3 a b c d X0X4 Z0(x) Figure 5.3: Schematic representation of a backwater curve approximated by uniform regimes Using this approach, the solution of (5.6) is given in the following manner: "q(X, s) z(X, s)#= Γ X, 0"q(0, s) z(0, s)#="γ11(s)γ12(s) γ21(s)γ22(s)#"q(0, s) z(0, s)#(5.7) where the transfer function matrix Γshould be calculated using: Γxn=X, x0= 0= 0 Y k=n−1 eA(xk, s)hk(5.8) with hk= (xk+1 −xk). It should be noted that the term eA(xk,s)hkcorrespond to the exponential of a matrix (evaluated at xk), which yields in this case: eA(xk, s)hk=      λ2(s)eλ1(s)hk−λ1(s)eλ2(s)hk λ2(s)−λ1(s) ceλ1(s)hk−eλ2(s)hks λ2(s)−λ1(s) λ1(s)λ2(s)eλ2(s)hk−eλ1(s)hk c(λ2(s)−λ1(s)) s λ2(s)eλ2(s)hk−λ1(s)eλ1(s)hk λ2(s)−λ1(s)      (5.9) where λ1,2(s) = 1 fas +b±pcs2+ds +e 5.2. Model of a pool 75 with a= 2B0(xk)V0(xk) b=γ0(xk) c= 4C02(xk)B02(xk) d= 4B0(xk)V0(xk)γ0(xk)−C02(xk)−V02(xk)B0(xk)β0(xk) e=γ02(xk) f= 2B0(xk)C02(xk)−V02(xk) In the solution expressed by (5.7) using (5.8) and (5.9), sis the laplace variable, Xis the x position of the end of the transport area (the start of the storage area) and Γis the transfer function matrix that describes exactly all the input-output relationships of the linearized SaintVenant equations in the laplace domain. In fact, it is possible with this formula to determine the value of the water flow and of the water depth at X, by knowing their values at the start of the canal and by knowing the transfer function matrix. Now, since usually the water depths are the outputs and the water flows the inputs of a hydraulic model, it is convenient to express this matrix relationship in the same manner. This is performed with basic algebraic matrix manipulations to (5.7), which in this particular case yields:    z(0, s) z(X, s)   =     −γ11(s) γ12(s) 1 γ12(s) γ21(s)−γ22(s)γ11(s) γ12(s) γ22(s) γ12(s)        q(0, s) q(X, s)   (5.10) With (5.10), the water depth at the beginning and at the end of the transport area depend, in a deterministic manner, on the water flow that enters and that leaves that area. This is the end of the derivation presented in Litrico and Fromion (2004c) to solve the linearized Saint-Venant equations. However, this work aims to extend this linear model to include the influence of the storage at the end of the pool and to consider explicitly the offtake discharge as an independent variable. It is thought that this addition can be carried out in the following way. To complete the model, it is necessary to linearize the storage equation given by (5.3). Following the same procedure as with the Saint-Venant equations, we obtain: qi(t)−qi+1(t)−qL i(t) = dzs i(t) dt As i (5.11) Applying the Laplace transform to (5.11) and replacing it in the expression for z(X, s)from 76 Chapter 5. Canal Identification for Control Purposes (5.10) assuming that zs i(s)≈z(X, s), the following equation is obtained after reordering: zs i(s) = γ21(s)γ12(s)−γ22(s)γ11(s) γ12(s)−s γ22(s)As i qi(s) + γ22(s) γ12(s)−s γ22(s)As i (qi+1(s) + qL i(s)) (5.12) (5.12) represents a linearized model (directly derived from the Saint-Venant equations) for a reach working around an operational point. It can be seen that the water level of interest zs i (the one where water is diverted for irrigation) can be obtained, if the water discharges that enter (qi) and that exit the reach (qi+1 and qL i) are known. However, in order to obtain the transfer functions that relate those variables, it is necessary to know a considerable amount of information, including the design parameters and the water depths of the reach, all of them at small enough longitudinal discrete positions of the reach (in order to reproduce accurately the characteristics of interest of the reach). Since the goal of this work is to design an appropriate and simple black-box (without information about physical parameters) modeling procedure for a canal reach, that can fulfill the control automation requirements, this model is not going to be used explicitly. However, it will be used to study the main properties that a simpler model should have and to decide which model structure is more adequate for the purposes established. 5.2.3 Properties of the linearized model It is extremely important to know the transfer function characteristics of a system (e.g poles, zeros, etc.) for identification and control purposes. However only the uniform regime case has a clear analytical expression that can be analyzed. This is absolutely true, but an analytical expression can always be derived, using the fact that any shape of backwater curve of any type of reach can be well approximated using different number of terms in (5.8). Using this approach and after some manipulations, it can be conjectured by empirical induction (see Appendix A) that the structure of (5.12) will always be (no matter the reach or the operational condition) of the form: zs i(s) = 1 sn1(s) d1(s)d2(s)qi(s)−n2(s) d2(s)qi+1(s)−n2(s) d2(s)qL i(s)(5.13) Expression (5.13) has been developed considering that the offtake is located at the end of the reach. However, if the offtake is located in between the pool, a similar structure can be obtained (see also Appendix A): zs i(s) = 1 sn1(s) d1(s)d2(s)d3(s)qi(s)−n3(s) d3(s)qi+1(s)−n2(s) d2(s)d3(s)qL i(s)(5.14) In (5.13) and (5.14), n1(s),n2(s)and n3(s)are irrational numerator expressions, d1(s), 5.2. Model of a pool 77 d2(s)and d3(s)are irrational denominator expressions (all of them include exponentiation and roots of spolynomials) and 1 scorrespond to an integration in the time-domain (integrator pole). From model structures (5.13) and (5.14) several conclusions can be drawn. Some of them are the following: •The existence of an integrator pole (real pole in the origin) denotes that the system is marginally stable. If any of the inputs of the system is excited with a finite impulse input, the output magnitude will be bounded. However, if the system is given a step as an input, the system’s output could increase indefinitely. So, this system is not a Bounded-Input Bounded-Output (BIBO) system. For system identification and control designs, this type of systems should be treated with special care; otherwise very bad performance behaviors could appear. •When the offtake is located at the end of the reach, its transfer function and the one of the discharge that enters the next reach, i.e. the transfer functions of qL i and qi+1 are identical. That means that for model system identification, it is enough to identify the transfer function of one of them to know the other. However, when the offtake is somewhere else, their transfer functions are different because of their numerators. •The irrational terms of all the transfer functions imply that an approximation by rational transfer functions (with a Padé approximations for example) would be more or less accurate, depending on the number of terms used to approximate the irrationality. Hence, it is expected that the rational transfer function approximation would have more terms than the irrational original one. •The denominators of all the transfer functions share common terms. That means that the dynamical responses obtained by the inputs of the models are similar. Technically speaking, the transfer functions of the model have some poles in common. The appearance of the integrator pole, or in other words, that a reach have similarities with a swimming pool or tank is not a real surprise and is, in some sense, expected. As mentioned before, this pole appears clearly in the uniform case regime and has been successfully included in several simplified models proposed by other researchers (Integrator Delay (ID) model (Schuurmans et al., 1999b), Integrator Delay Zero (IDZ) model (Litrico and Fromion, 2004a), etc.). However, there are some research works that did not include this characteristic in their models. That is why it was found important to go a step forward in the generalization of this feature for any type of reach (slope, cross section, width, length, etc.) working around any flow condition (discharge, backwater curve, etc.). 78 Chapter 5. Canal Identification for Control Purposes 5.2.4 Characteristics of some type of pools In addition to the general properties presented above, there are also some particular characteristics that are worth to review. These characteristics are going to be studied, analyzing the Bode diagram (or frequency response) of the downstream water level zs i using different discharges for the water inflow qias in Litrico and Fromion (2004c). The Bode diagram is a logarithmic magnitude and phase plot of a transfer function, that gives information of this function evaluated in the s-plane imaginary axis. In a more practical view, it shows what happens with the amplitude and the phase of the response of a system, when it is excited with a sinusoidal input at a given frequency. In this case, the Bode diagram is obtained numerically around a given operational backwater curve (by means of the linearized Saint-Venant equations), calculating the whole matrix product series (5.8) for each frequency point s=jω and replacing the results in (5.12). There is a natural question that arise at this point: why to study the Bode diagram of a pool that is going to be identified? The answer to this question is the following: there are some pool characteristics that have a direct relationship with some modeling issues. Thus, the Bode diagram can be useful to get an insight into: the required number of parameters when using a linear model structure (model order), the validity range of a particular linear model, the number of models required to cover a given operating range, etc.. For example, a pool exhibiting high resonant peaks in the Bode diagram due to the effect of water waves traveling back and forth through it, requires a model structure with more parameters (higher order) than a pool that does not exhibit this behavior. Moreover, it is known that operating points of the same pool could exhibit considerable differences in the steady state gain, amount of delay, etc.. All this information can be obtained analyzing the Bode diagram of a pool. To illustrate these points, two pool configuration are going to be studied: one short flat pool and one long slopping pool. 5.2.4.1 Short flat pools The Bode diagram of such a pool is presented in figure 5.4 for three different operational discharges, namely 14 m3/s,7 m3/sand 1.75 m3/s. In order to carry out a detailed analysis of the diagram, several areas have been marked (dotted ellipses): (1) This area shows that in the very long-term, a pool like this acts as a pure integrator (swimmingpool or tank). However, the gain of this integrator varies for different working conditions (in this case the discharge value). 5.2. Model of a pool 79 10−5 10−4 10−3 10−2 10−1 −40 −30 −20 −10 0 10 20 Frequency (rad/s) Magnitude (dB) 14 m3/s 7 m3/s 1.75 m3/s (1) (2) (3) (a) Bode Magnitude 10−5 10−4 10−3 10−2 10−1 −4500 −4000 −3500 −3000 −2500 −2000 −1500 −1000 −500 0 Frequency (rad/s) Phase (deg) 14 m3/s 7 m3/s 1.75 m3/s (4) (b) Bode Phase Figure 5.4: Bode diagram between qiand zs i for short flat pool 80 Chapter 5. Canal Identification for Control Purposes (2) The change in the slope of the diagram in this area reveals the appearance of a zero in the transfer function. In this case, this zero is related to the propagation of a shock wave through the pool when there is a change in the water inflow or in the water outflow. (3) Each oscillation in this diagram denotes the existence of a resonant mode. That means that the water depth of this type of pools have a natural tendency to oscillation. However, the magnitude and frequency of this oscillations vary depending on the discharge and on the water level. These resonant modes have a close relationship with the shock-waves traveling back and forth through the pool. As a matter of fact, they occur approximately at a frequency equal to the time that takes the shock wave to go and return from one extreme of the pool to the other. (4) A decreasing curve in the phase margin reveals the existence of a delay between the input and the output of a system. In this case, the time that takes a discharge change at the beginning of the pool to modify the downstream water level. The amount of delay also changes for different discharge values and water level conditions. 5.2.4.2 Long slopping pools A typical Bode diagram of this type of pools is presented in figure 5.5 for three different operational discharges, namely 80 m3/s,40 m3/sand 10 m3/s. The result of the analysis of the dotted areas of this Bode diagram is given below: (1) This zone shows the integral part of the behavior. Hence, this system also resembles a water tank under certain conditions. As expected, this behavior appears in the long-term response, but the integral gain changes with the discharge or the water level. (2) This area shows some changes in the slope of the diagram. First, there is an increment in the slope attributable to the presence of another pole. Finally, the curve rises because of the presence of some zeros in the system. The whole set models the presence of a shock-wave when changing the water inflow. This shock-wave travels through the whole pool until arriving to the end of the pool. The distortion of this shock-wave through the way is also determined by this zone of the diagram. (3) This type of pools does not develop resonant modes in the frequency response. That means that the shock waves that arrive to the end of the pool are not reflected back. (4) The decreasing curve in the phase margin is an evidence of the existence of delay in the system. As a matter of fact, the frequency were the phase is 360◦is the inverse of the delay in seconds. Thus, the deviation of these curves for different discharges show variations of the delay for different discharge conditions. 5.2. Model of a pool 81 10−5 10−4 10−3 10−2 10−1 −80 −60 −40 −20 0 20 Frequency (rad/s) Magnitude (dB) 80 m3/s 40 m3/s 10 m3/s (1) (2) (3) (a) Bode Magnitude 10−5 10−4 10−3 10−2 10−1 −9000 −8000 −7000 −6000 −5000 −4000 −3000 −2000 −1000 0 Frequency (rad/s) Phase (deg) 80 m3/s 40 m3/s 10 m3/s (4) (b) Bode Phase Figure 5.5: Bode diagram between qiand zs i for long slopping pool 88 Chapter 5. Canal Identification for Control Purposes Table 5.1: Reach’s Characteristics Characteristic Value Pool length 3000 m Bottom slope 0.002 m/m Bottom width 7 m Pool shape trapezoidal Side slope 1.5 m/m Manning’s n 0.014 Operational flow 10 m3/s 0 2000 4000 6000 8000 10000 12000 14000 16000 18000 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Time (s) Water Level (m) 680s 2300s 3900s (a) Downstream water level response 0 1000 2000 3000 4000 5000 6000 0 1 2 3 4 5 6 x 10−4 Time (s) Water Level Variation (m) (b) Derivative of the water level response Figure 5.9: Water level response for an inlet flow step change for the reach of table 5.1 5.3. System identification of a pool 89 10−510−410−310−210−1 −40 −35 −30 −25 −20 −15 −10 −5 0 5 10 Frequency (rad/s) Magnitude (dB) (a) Bode Magnitude 10−510−410−310−210−1 −4000 −3500 −3000 −2500 −2000 −1500 −1000 −500 0 Frequency (rad/s) Phase (degrees) (b) Bode Phase Figure 5.10: Bode diagram between qiand zs i for the reach of table 5.1 The step response is obtained simulating, in a Saint-Venant based model, the downstream water level produced by a step increment in the inflow of this pool, maintaining the outflow constant. The time history of the step response and the Bode diagram (obtained as in Section 5.2.4) for this particular configuration are presented next. First of all, figure 5.9 shows that after the inflow change, there is a delay of approximately 680 seconds until noticing a change in the water level. At that time, a step increment occurs in the water level. Then, the level rises with a variable slope (including a big sudden slope change at time 2300 s) until time 4000 seconds. At that time the water level continue infinitely raising at a constant rate. Looking in detail at the step response obtained, one can conclude that there are water level 90 Chapter 5. Canal Identification for Control Purposes variations (for example the first increment at 680 seconds) that, in order to be accurately reproduced, require the use of a very small sampling period. That would be the only way to sample the dynamical change with at least 4 points as recommended before. However, the variations do not have a clear and predictable tendency until time 4000 seconds. A linear discrete model that has to approximate this transient behavior with a so short sampling period, would have a great amount of parameters. So, a trade off should be made in the selection of the sampling period between accuracy and model complexity, having in mind all the possible control design and numerical problems discussed before. The Bode diagram of figure 5.10 reveals that the linearized model has poles from very low to very high frequencies, standing out the dominance of the pole at the origin (integrator pole) as derived in Section 5.2.3. This is due to the resonant modes of the pool. The existence of poles in an infinite range of frequencies shows that it is impossible to reproduce the exact dynamical behavior of the system by sampling it, because of the fact that the sampling frequency should lie 10 times away from the highest frequency present in the system. The only possibility is to approximate a smaller frequency band, i.e. to loose some fast changing dynamics. If this is the case, it is very important to remember that an approximation covering frequencies above the first resonant frequency, would need an exponentially increasing number of parameters. After this analysis, a good choice for this particular case is a sampling time close to 221 seconds. It was selected having in mind the time needed by the water level to achieve a fixed course divided by 15: 4000 s −680 s 15 = 221 s With this sampling time, a correctly obtained linear model should be able to: 1. Reproduce the general tendency of the response properly. 2. Cover at least the first resonant modes. 3. Approximate the response at the gain crossover frequency. 4. Have a delay made up of few sampling times (in this case 3), to facilitate the control design and avoid possible performance problems. Nevertheless, it is still probable that the model has unstable zeros making the model to have a non-minimum phase behavior. Unstable zeros can arise if the sampling time does not divide the delay period exactly. In conclusion, the choice of the sampling period should always be taken thinking in the particular system dynamics and in the intended purpose, using all the design knowledge at hand. 5.3. System identification of a pool 91 This way of proceeding is indispensable to attain a linear model with the ability to adequately reproduce the hydraulic behavior of a pool. 5.3.1.2 Discrete transfer functions It is important to remark that a discrete-time model will never be the same as a continuous one; it will only be an approximation to a model with similar characteristics, despite the technique used. In general, a discrete-time representation can be viewed as the z-transform of the discretetime impulse response of a system for a given sample time. In this case, transforming the continuous models (5.13) or (5.14) into a z-transform model would give: zs i(z) = F1(z)qi(z)−F2(z)qi+1(z)−F3(z)qL i(z)(5.16) where zis the z-transform complex variable.2 The model given by (5.16) still maintains some characteristics of the continuous one. For example, it is also applicable to any type of reach. In addition, there are two aspects that have a direct correspondence in the discrete-time representation: 1. Any Laplace domain pole will have a direct counterpart in the z-transform domain with the relation given by pz=epsT, where psis the pole in the Laplace domain, Tis the sampling period and pzis the discrete-time pole. Because of this, the pole at s= 0 (integrator pole) of 1 swill produce always a discrete-time pole at z= 1 (discrete-time integrator pole); that means that a term 1 z−1will always appear in the z-transform of the model. 2. The time delays between the inputs (water discharges) and the output (water level) of the continuous time model will satisfy the following property: Z{f(k−d)}=z−dF(z) where f(k−d)is the delayed discrete-time impulse response of the model, kis the discrete time instant variable, dis the delay expressed in amount of instants and zis the z-transform variable. This implies that the z-transform model (5.16) will always appear in the following manner: zs i(z) = z−d1F10(z)qi(z)−z−d2F20(z)qi+1(z)−z−d3F30(z)qL i(z)(5.17) 2In the following any F(z) will correspond to Z {f(kT )}, the Laplace transform of a sampled time function f(kT )(called z-transform), with k= 0,1,2, . . . the sampling instant and Tthe sampling period. Besides for notation simplicity f(kT )will normally be written as f(k)only. 92 Chapter 5. Canal Identification for Control Purposes where d1,d2and d3are the time periods, measured in discrete time instants, that takes each discharge to influence the downstream water level. Point 1 presents a property that have to be taken with care in order to avoid possible modeling problems. Dynamically and numerically it is very difficult to identify a discrete-time model with a pole exactly located at z= 1. The problem is that a small variation in its position leads to a completely different dynamical behavior. For instance, a pole at z= 1.01 would produce an unstable model, i.e. the model output can go quickly to the infinite for finite input values. Conversely, a pole at z= 0.99 would produce a strictly stable model, whose response will always be bounded. There are three approaches normally taken with respect to this problem when estimating a model by means of system identification: 1. To forget about the problem and obtain a model anyway. 2. To identify a model and afterwards correct the position of the pole in the estimated model. 3. To acknowledge the existence of the pole and apply its influence directly to the data in order to identify the other components of the model. This can be achieved in any of the two following ways: y(z) = F(z)u(z) = 1 z−1F0(z)u(z) ⇒     y(z) = F0(z)1 z−1u(z)=F0(z)u0(z) y0(z) = y(z)(z−1) = F0(z)u(z) (5.18) The first way is generally more recommendable (especially in the presence of noise) and is equivalent to make a cumulative sum of the input data. The second way correspond to a differentiation of the output data. Unfortunately, both procedures modify the frequency characteristic of the original input signal. As a consequence, a maximally informative input signal specially designed to identify a model would lose it optimal properties. From the three approaches presented above, 2 and 3 are the better ones. However, 3 has an additional drawback when modeling reaches; this procedure applied to large sample periods can induce long-term modeling bias errors. Therefore, the second option of identify and then correct would be the recommended one. This method is not always easy to apply to some discrete-time model structures. Moreover, there are some types of model structures that cannot deal with processes with integrators. For those cases the only choice is 3. 5.3. System identification of a pool 93 5.3.2 Discrete-time model structures Generally, linear discrete-time models can be divided into three main classes: Discrete transfer function models: Models based only on the input-output characteristics of a system by means of the z-transform. Discrete state-space models: Time-domain models that incorporate all the information about the internal dynamics of a system. Orthonormal basis models: Models that make use of the special approximation properties of some basis functions to mathematically represent systems. Depending on the particular class and structure chosen, there are different types of parameter estimation methods. Examples of types of parameter estimation methods are Subspace methods for estimating state-space models, Prediction Error methods and Output Error methods for estimating transfer function models, etc.. These methods can estimate the parameter values of a model from only an enough informative data set collected from the true system. In this section, only two model structures for modeling pools are going to be considered: 1. The Auto-Regressive with eXogenous Input (ARX) model, a transfer function based model. 2. The Laguerre model, an orthonormal basis based model. The chosen models can only have a finite number of rational elements. This option has been taken despite the fact that the process is governed by irrational transfer functions. An approximation like this can be carried out because an irrational term can generally be approximated by a linear combination of rational ones, like in a Padè approximation of a function. However, depending on the particular irrational term, a good approximation can require a high number of rational terms to achieve good results. In the following, all model parameters are going to be estimated using standard least-squares based algorithms. The mathematical formulation of the ARX and Laguerre discrete-time pool models are explained in detail in the next sections. 5.3.2.1 ARX model To derive the ARX model, it is first necessary to introduce the forward shift operator qand the backward shift operator q−1respectively as: qf(k) = f(k+ 1), q−1f(k) = f(k−1) 94 Chapter 5. Canal Identification for Control Purposes Then, assuming that F10(z),F20(z)and F30(z)of (5.17) can be approximated by quotients of polynomials in the following way: F10(z) = B1(z) A(z), F20(z) = B2(z) A(z), F30(z) = B3(z) A(z) with A(z) = 1 + a1z−1+. . . +anaz−na B1(z) = b11+b12z−1+. . . +b1nb1z−nb1+1 B2(z) = b21+b22z−1+. . . +b2nb2z−nb2+1 B3(z) = b31+b32z−1+. . . +b3nb3z−nb3+1 and replacing in (5.17), that yields after reordering: A(z)zs i(z) = z−d1B1(z)qi(z)−z−d2B2(z)qi+1(z)−z−d3B3(z)qL i(z) Finally, after applying the inverse z-transform to move the problem back in the time domain, the general model structure is: A(q)zs i(k) = B1(q)qi(k−d1)−B2(q)qi+1(k−d2)−B3(q)qL i(k−d3)(5.19) with A(q) = 1 + a1q−1+. . . +anaq−na B1(q) = b11+b12q−1+. . . +b1nb1q−nb1+1 B2(q) = b21+b22q−1+. . . +b2nb2q−nb2+1 B3(q) = b31+b32q−1+. . . +b3nb3q−nb3+1 As can be seen, the parameters of the ARX model (5.19) are: the polynomial orders na, nb1,nb2and nb3, the polynomial coefficients, and the a priori known time delays d1,d2and d3 expressed in sampling instants. Once the parameters are determined, this model can calculate the downstream water level zs i at instant kby a weighted sum of measurement data from past instants. Specifically, the historical data that is needed are the water levels zs i and the water discharges qi,qi+1 qL i, collected at a sampling interval T. In some sense it is a little bit restrictive to force the three transfer functions, namely B1(z) A(z), B2(z) A(z), and B3(z) A(z)to have the same denominator polynomial A(z). From (5.13) or (5.14) it is noticeable that they do not have really the same denominator. However, (5.13) or (5.14) show 5.3. System identification of a pool 95 in the same way that their transfer functions share some poles. Therefore, the assumption is not totally wrong. The reason for using the same denominator polynomial A(q)have its roots in simplifying the parameter estimation process, because in that way it can be performed using a linear least squares approach. This model can be applied to stable and unstable processes so, as mentioned before, there are two options in relation to the integrator parameter estimation problem. One is to estimate the model parameters and then fix the integrator pole location. This can be done following this procedure: 1. To perform the parameter estimation process obtaining A(q),B1(q),B2(q)and B3(q). 2. To calculate the roots of A(q). 3. To replace the root of A(q)that is close to 1by exactly a 1. 4. To form a new A(q)polynomial keeping the other roots locations. The second option is: 1. To apply the cumulative sum on each input variable (discharges). 2. To estimate the model polynomials A(q),B1(q),B2(q)and B3(q). 3. To multiply the estimated A(q)polynomial by (1 −q−1). 5.3.2.2 Laguerre model This model is based on the Laguerre functions, a complete orthonormal set of functions in L2(0,∞), the space of square Lebesgue integrable functions in the (0,∞)interval (Zervos and Dumont, 1988). These functions are described in the time domain by: li(t) = p2pept (i−1)! di−1 dti−1ti−1e−2pt where iis the order of the function (i≥1) and pis a positive parameter. The laplace transform of the Laguerre functions produces rational functions in the svariable of the following form: Li(s) = p2p(s−p)i−1 (s+p)i By using a linear combination of a truncated number of these functions, any impulse response (or its associated transfer function) that belongs to the intersection of L1(0,∞)∩L2(0,∞), can be approximated as follows: f(t) = N X i=1 cili(t) = cTlF(s) = N X i=1 ciLi(s) = cTL 96 Chapter 5. Canal Identification for Control Purposes cT=hc1c2··· cNi lT=hl1(t)l2(t)··· lN(t)i LT=hL1(s)L2(s)··· LN(s)i A discrete-time state-space version of this model can be obtained by applying a continuous network compensation method to each transfer function (Zervos and Dumont, 1988). The result of this operation yields: l(k+ 1) = A l(k) + Bu(k) y(k) = cTl(k)(5.20) where l(k)is the state vector of order N,u(k)is the system input and y(k)is the system output. In addition, if Tis the discrete sampling time, Aand Bcan be defined as: A=               τ10··· 0 −τ1τ2−τ3 Tτ1 .... . . . . .......0 (−1)N−1τN−2 2(τ1τ2+τ3) TN−1··· −τ1τ2−τ3 Tτ1               BT=τ4−τ2 Tτ4··· −τ2 TN−1τ4 with τ1=e−pT ,τ2=T+2 p(τ1−1),τ3=−Tτ1−2 p(τ1−1) and τ4=√2p(1−τ1) p. Making use of the Laguerre model structure for a pool, the particular model gives: zs i(k) = c1T[A1 l1(k−1) + B1 qi(k−1)] −c2T[A2 l2(k−1) + B2 qi+1(k−1)] −c3T[A3 l3(k−1) + B3 qL i(k−1)] (5.21) In this type of model, the model output is obtained with just the inputs and the Laguerre functions values at time k−1(always known for a particular p). The parameters of the Laguerre-based model (5.21) are the vectors of coefficients c1,c2 and c3; the designer-chosen number of Laguerre functions Nand the Laguerre pole value (p) for each transfer function (to calculate A1,A2,A3 and B1,B2,B3). It is not necessary an a priori knowledge on the process orders or on any delay and the response produced by each input can be adjusted in a totally independent manner. On the other hand, for increasing complexity transfer functions, it is necessary a higher number of terms. 5.3. System identification of a pool 97 This model can not approximate systems that are not strictly stable. Hence, there is only one way to identify the model in this case: 1. To apply the cumulative sum (5.18) on each input variable (discharges). 2. To estimate the model parameters c1,c2 and c3. 3. To augment the state-space representation in order to include the integrator. 5.3.3 Experiment design: Input signal An input signal should excite a system with a rich frequency content, in order to successfully identify a model of a true real system. A rich frequency content means that the input signal can be decomposed into many frequencies of similar magnitude. In practice, it is appropriate to decide first upon which is the frequency band where it is important to have information about the system and then select a signal that has a more or less flat spectrum over it. Such an input is provided by the Pseudo Random Binary Sequence (PRBS) signal (Bialasiewicz, 1995). The sequence is generated by a digital waveform generator, which produces a binary signal by switching randomly between two output levels (a,−a). It owes its name "pseudo-random" to the fact that it is characterized by a sequence length within which the pulse width varies randomly, while it is periodic over a large time period. The PRBSs are generated by means of shift registers with feedback (implemented in hardware or software). The period is defined by the maximum sequence length: L= 2N−1(5.22) where Nis the number of stages of the shift register. As an example, a fragment of a PRBS signal is presented in figure 5.11. Assuming that u(k)is a random binary process with a current value of aor −aand that the value of u(k)can change every Tprbs seconds, namely Tprbs is the switching period, the corresponding spectral density of the signal is: Suu =a2 π Tprbs 2sin (ωTprbs/2) ωTprbs/22 (5.23) The spectral density function (5.23), can be assumed to be approximately flat up to a frequency about 0.3fprbs rad/s. If fprbs is sufficiently high (as compared to the bandwidth of a plant to be identified), then the random binary process has a spectrum corresponding to a broad-band noise. The spectrum of the pseudo-random binary signal is therefore an approximation of a broad band noise, provided that its clock frequency and its sequence length is large enough. However, there are two thinks to satisfy in order to enjoy the good properties of this signal: 104 Chapter 5. Canal Identification for Control Purposes As can be observed in figure 5.15, the experiment was designed to achieve the following goals: •To have the most information about the system dynamics. •To have sufficient data points to feed the parameter estimation algorithms. •To keep the water level around the operational point (2.1 m) in order to avoid the violation of the linearity assumption. With the generated data set, the following models were obtained: •ARX model with integrator position correction: A(q)=1−0.7791 q−1−0.03589 q−2−0.03576 q−3−0.01423 q−4 −0.01693 q−5−0.01205 q−6−0.2234 q−7+ 0.1174 q−8 q−d1B1(q) = q−40.0229 −0.0113 q−1(5.29) q−d2B2(q) = q−1−0.02598 + 0.01497 q−1 •ARX model with a posteriori added integrator: A(q)=1−1.199 q−1+ 0.1322 q−2+ 0.01203 q−3+ 0.03275 q−4 + 0.005332 q−5+ 0.009966 q−6−0.2101 q−7+ 0.2167 q−8 q−d1B1(q) = q−40.02262 −0.01995 q−1(5.30) q−d2B2(q) = q−1−0.02617 + 0.02631 q−1−0.002822 q−2 •Laguerre model without integrator: p1=0.0115 c1T=h−0.0001 −0.0017 −0.0083 −0.0168 −0.0127 −0.0032i p2=0.0053 (5.31) c2T=h−0.0015 −0.0022 −0.0014 −0.0008 −0.0007 −0.0003i In order to make an objective judgment, all model structures had been chosen so as to have to identify 12 unknown parameters. The orders of the polynomials and the delays of models (5.29) and (5.30) were obtained comparing the results given by different combinations of them and then picking out the most appropriate ones. 5.3. System identification of a pool 105 In the case of model (5.31), the Laguerre pole pwas determined by a Newton - Raphson iterative technique (Malti et al., 1998). This algorithm performed an optimal search in order to obtain the pole value with best fit. All model parameters where obtained using least squares based parameter estimation algorithms (Ljung, 1999). With respect to the integrator pole, model (5.29) includes the correction of the integrator pole position and model (5.30) was identified discarding the integrator in a first stage and then manually including it in A(q). Model (5.31) is the Laguerre model obtained without integrator; the integrator had to be included augmenting the states of the identified model. To compare the performances of the obtained models, they were tested in two different domains: the time domain and the frequency domain. A good performance in the time domain ensures that the model can accurately model a behavior for a given input, while a good fit in the frequency domain ensures that the identified model will approximate with good results, a response induced by a more general type of input. In the time domain, the step response of a Saint-Venant modeled pool was compared against the step response of the identified linear models. This test was performed producing a step increase in the upstream water discharge, maintaining the downstream water discharge value constant. The results are given in figure 5.16. 0 10 20 30 40 50 60 70 0 0.02 0.04 0.06 0.08 0.1 0.12 Time (minutes) Water level deviation (m) "True" Laguerre model ARX model int. corr. ARX model int. add. Figure 5.16: Step response of the reach v/s step response of the identified linear models As can be observed in figure 5.16, all models performed relatively well in approximating the behavior of the Saint - Venant based simulation. However, looking in detail, there are two features that are worthy of note: •The Laguerre model was less accurate in simulating the changes due to the water waves. •After a large time period, all models tended slowly to deviate from the true response. 106 Chapter 5. Canal Identification for Control Purposes In order to study if the models can only perform well for a step-type input or if they are really representing the system in question, it was found useful to construct the frequency responses (or Bode diagrams) of the Saint - Venant modeled pool and of the identified linear models. In this case it is required the analysis of two relationships: 1. Upstream discharge (qi)→Downstream water level (zs i) 2. Downstream discharge (qi+1)→Downstream water level (zs i) 1 and 2 are presented in figure 5.17 and figure 5.18 respectively. 10−510−410−310−2 −40 −35 −30 −25 −20 −15 −10 −5 0 5 10 Frequency (rad/s) Magnitude (dB) Laguerre model ARX model int. add. ARX model int. corr. "True" (a) Bode Magnitude 10−410−310−2 −1200 −1000 −800 −600 −400 −200 0 Frequency (rad/s) Phase (degrees) Laguerre model ARX model int. add. ARX model int. corr. "True" (b) Bode Phase Figure 5.17: Bode diagram between qiand zs i for pool of Table 5.1 5.3. System identification of a pool 107 10−510−410−310−2 −40 −35 −30 −25 −20 −15 −10 −5 0 5 10 Frequency (rad/s) Magnitude (dB) Laguerre model ARX model int. add. ARX model int. corr. "True" (a) Bode Magnitude 10−510−410−310−2 −50 0 50 100 150 200 Frequency (rad/s) Phase (degrees) Laguerre model ARX model int. add ARX model int. corr "True" (b) Bode Phase Figure 5.18: Bode diagram between qi+1 and zs i for pool of Table 5.1 Figures 5.17 and 5.18 show that there are no noticeable magnitude or phase errors for frequencies below 0.0005 rad/s. Hence, all the identified models approximate the pool behavior with high accuracy in the low frequency region. In the high frequency region, the magnitude plots show that the ARX models can approximate the location and magnitude of the resonant modes of the system with a medium-to-high accuracy until the Nyquist frequency (represented with a vertical line). In the same area, the Laguerre model tends to average the magnitudes of the frequencies, suggesting that the number of Laguerre functions is insufficient. In fact, doubling the number of terms improves the approximation of the resonant modes in a high degree (see figure 5.19). On the other hand, a zoom on the high frequency region of the phase plots reveals the existence of notable phase errors. These errors appear at approximately 1/10 of the Nyquist fre- 108 Chapter 5. Canal Identification for Control Purposes 10−510−410−310−2 −40 −35 −30 −25 −20 −15 −10 −5 0 5 10 Frequency (rad/s) Magnitude (dB) "True" Laguerre model Figure 5.19: Bode diagram of Laguerre reach model with 24 parameters quency and become more important when approaching to the Nyquist frequency. Nevertheless, this is a normal behavior of some discrete-time models working near the Nyquist frequency. 5.4. Canal model synthesis 109 5.4 Canal model synthesis In the previous sections, the issues concerning to pool models and their identification have been deeply studied. Now, it is the turn to explain how to construct a complete canal model. For the purpose proposed, it is enough to describe a canal in a simplified way. A general scheme of a canal with Npools is presented in figure 5.20. q1 qL 2 qL 1 qL N-1 zs 1 zs 2 zs N-1 Pool 2 Pool N Pool 1 q2q3 zs N qN qN+1 qL N Reservoir Figure 5.20: Canal scheme with Npools The identification of each pool of the canal would produce several models with a structure like the one presented in (5.16). It is straightforward to show that the interconnection of all these discrete models would yield a discrete transfer function matrix of the form:               zs1(z) zs2(z) . . . zs N (z)               =               F11(z)F22(z) 0 ··· 0 0F12(z)F23(z).... . . . . ..........0 0··· 0F1N(z)F2N+1(z)                                  q1(z) q2(z) . . . qN(z) qN+1(z)                    +             F31(z) 0 ··· 0 0F32(z).... . . . . .......0 0··· 0F3N(z)                         qL1(z) qL2(z) . . . qL N (z)             (5.32) 110 Chapter 5. Canal Identification for Control Purposes The canal model structure presented in (5.32) is one way to express a Multiple-Input MultipleOutput (MIMO) process. Although its accurate analysis depends in a high degree on the particular canal characteristics, one can highlight the following points in general lines: •Although it is possible that the interconnection can modify some individual pool dynamics, the pool properties and characteristics presented in other sections are still extrapolative to a complete irrigation canal. In other words, the system is likely to have delays, resonance, marginal stability, etc.. •There is an additional characteristic that appears from the union of pool models in (5.32): system coupling. That means that a single system input does not only influence a single output; it produces an action over other outputs too. In this case (5.32) shows that a particular gate discharge qiacts over two water levels: the one immediately upstream (zs i−1) and the one located in the downstream end of the pool (zs i). As explained in other sections, this canal model can be easily extended to include gate openings, weirs or other hydraulic structures by replacing discharges with appropriate hydraulic formulas. Those modifications would normally increase the coupling degree of the system (less zeros in (5.32)) and change some transient behaviors. 5.5. Experimental application 111 5.5 Experimental application The general approach of this work has been to develop knowledge on a practical application basis. This section validates the use of system identification in canals by testing its application on the Canal PAC-UPC. This canal was already presented in Chapter 3 and is depicted in figure 5.21. Reservoir Gate1 Gate3 Gate5 Weir3 Weir4 Weir1 POOL1 POOL2 POOL3 87.0m 90.2m 43.5m Zs 1 Zs 2 Zs 3 QL 1 QL 2 Q1Q2Q3 Figure 5.21: Scheme of the Canal PAC-UPC Specifically, the recommendations and procedures developed in this chapter are applied in a modeling stage in order to produce an ARX MIMO model of this laboratory canal. The necessary data comes from experiments performed on the real canal and the identification procedure is explained in the following subsections. 5.5.1 Sampling time In order to determine an appropriate sampling time, it has been found helpful to obtain a step response of the system. For example, a true 20 L/sdischarge increment at the upstream end of Pool 1 produced the downstream level variation plotted in figure 5.22. Figure 5.22 marks two important time periods: •Pool delay = 28 s •Tendency stabilization = 130 s The pool delay is time needed to see a downstream water level variation after an inflow change at the beginning of the pool. The tendency stabilization refers to the subsequent time period needed to observe a clear tendency in the downstream water level. The pool sampling time selection procedure proposed in this thesis seeks a good trade off between accuracy, simplicity and control goals. One recommendation was to pick out a sam- 112 Chapter 5. Canal Identification for Control Purposes 0 30 60 90 120 150 180 210 240 0 2 4 6 8 10 Time, t (s) Down. water level variation, zs 1 (cm) 28s 130s Figure 5.22: Inflow step response - Pool 1 pling time that divides the delay 2 or 3 times. Other suggestion was to divide the stabilization time by enough points. In this case it was found a good alternative to divide the stabilization time by 13: •Sampling time = 130 s /13 =10 s This choice gives an acceptable discretization of the transient behavior and divides the delay by almost 3. Moreover, this choice gives a Nyquist frequency higher than the first resonant modes. In a broader perspective, it is convenient to have only one sampling period for all pools in a canal. This permits the conception of a canal model with a unique sampling time. The decision should be made analyzing the step responses of all pools and picking out an average value. In this case it was not necessary; the step responses of the rest of the pools are rather similar. Furthermore, the shortest pool have a delay of approximately 20 s. The sampling time is contained twice in this period. This decision is also adequate for control in this case. This value permits the use of a control time equal to the model sampling time. It was explained that this choice eliminates the need for control recalculation. This would not have been possible with a higher sampling frequency; it is not advisable to operate the servo-motorized gates of the canal at a higher rate. 5.5. Experimental application 113 5.5.2 PRBS design The guidelines presented in this chapter suggest: 1. To remain around an operational point. 2. To maintain the experiment as short as possible. 3. A maximum pulse width bigger than the stabilization time. 4. Between 200 and 1000 data points to feed the parameter estimation algorithm. In general, there are several PRBS signals that fulfill these requirements. Trial and error tests showed that a possible parameter set could be the one presented in table 5.2. Table 5.2: PRBS parameters Parameter Value PRBS time factor, p 4 Number of registers, N 6 PRBS amplitude, a (L/s) 10 These design parameters originate a PRBS signal with the following characteristics: •Data points: 252 points •Maximum pulse width: 240 s (>130 s) •Minimum pulse width: 40 s •Experiment duration: 2520 s Two PRBS signals are required to perform the identification step because the pools of the canal have only two independent manipulable inputs: qiand qi+1. (qL i is not entirely manipulable in this canal) Generally, it is found better to excite both inputs at a time in order to capture the essential dynamics of the system and to reduce the length of the experiment. Hence, two maximally shifted PRBSs were constructed so as to have two different and statistically independent signals. These signals are shown in figure 5.23. 120 Chapter 5. Canal Identification for Control Purposes properly described by the model and computes the autocorrelation function of the residuals to determine if the model correctly describes the process disturbance or, what is the same thing, to check if the residuals are mutually independent. 3. The selected ARX models were corrected in the cases where a pure integrator should appear in the model according to Section 5.2. This was carried out pursuant to page 92 as follows: (a) The roots of the Ai(q)polynomial were calculated. (b) The root qxclose to 1 was replaced with exactly an 1. (c) Ai(q)was reconstructed using this new root value and the rest of the roots. The result of the whole process gave the following ARX models: •Pool 1 ARX model: A1(q)zs1(k) = q−d11B11(q)q1(k) + q−d21B21(q)q2(k) + e1(k) A1(q) = 1 −0.6918 q−1+ 0.09173 q−2−0.07062 q−3+ 0.1638 q−4 −0.2801 q−5−0.213 q−6 q−d11B11(q) = q−31−0.2895 q−1(5.33) q−d21B21(q) = q−1(−1.154 + 0.4729 q−1−0.06765 q−2+ 0.4496 q−3 −0.4899 q−4) •Pool 2 ARX model: A2(q)zs2(k) = q−d12B12(q)q2(k) + q−d22B22(q)q3(k) + e2(k) A2(q)=1−0.04315 q−1−0.1178 q−2−0.03932 q−3+ 0.03022 q−4 −0.1962 q−5+ 0.05333 q−6−0.3033 q−7−0.3838 q−8 q−d12B12(q) = q−4(1.644) (5.34) q−d22B22(q) = q−1(−1.994 + 0.4728 q−1−0.1484 q−2+ 0.6077 q−3 −0.4835 q−4) •Pool 3 ARX model: A3(q)zs3(k) = q−d13B13(q)q3(k) + e3(k) A3(q)=1−0.6495 q−1−0.2344 q−2+ 0.142 q−3−0.3175 q−4 + 0.2112 q−5(5.35) q−d13B13(q) = q−21.422 −0.9872 q−1 5.5. Experimental application 121 5.5.6 Canal model The union of pool models (5.33), (5.34) and (5.35) with the incorporation of two intermediate weir regulated offtakes produce the following discrete multivariable ARX model (T= 10 s):     A1(q) 0 0 0A2(q) 0 0 0 A3(q)        zs1(k) zs2(k) zs3(k)    =    q−2B11(q)B21(q) 0 0q−3B12(q)B22(q) 0 0 q−1B13(q)        q1(k−1) q2(k−1) q3(k−1)     +    B21(q) 0 0B22(q) 0 0    "qL1(k−1) qL2(k−1)#+    e1(k) e2(k) e3(k)    (5.36) In some textbooks, it is also usual to rewrite a model like this in the following form:     zs1(k) zs2(k) zs3(k)    =    1−A1(q) 0 0 0 1 −A2(q) 0 0 0 1 −A3(q)        zs1(k−1) zs2(k−1) zs3(k−1)     +    q−2B11(q)B21(q) 0 0q−3B12(q)B22(q) 0 0 q−1B13(q)        q1(k−1) q2(k−1) q3(k−1)     +    B21(q) 0 0B22(q) 0 0    "qL1(k−1) qL2(k−1)#+    e1(k) e2(k) e3(k)     (5.37) This canal model can be used for several purposes. For example, it can be used to simulate downstream water level responses to beforehand known discharge curves in a fast way. Another possibility is to use it to predict future water level variations from past discharge measurements and past water level readings. Nevertheless, its simpleness makes it specially attractive for control applications. 5.5.7 Model validation A model validation stage determines the degree of validity of a model. It should be done by checking the model with data not used in its construction. In this case, the model was validated in the time domain and in the frequency domain. 122 Chapter 5. Canal Identification for Control Purposes 5.5.7.1 Time domain A totally new data set was taken from the Canal PAC-UPC in order to prove the effectiveness of the identified model. This experiment consisted in following the evolution of the canal variables after closing 3 cm Gate 1 with Weir 1 completely closed. This test was somewhat severe because it produced a33 % discharge reduction (∆Q1≈ −20 L/s). A variation like this, make the water levels go away from their initial conditions, occasionally violating the linearity assumption. However, the test was found appropriate in order to determine some model performance limits. In order to study the efficiency of the model, two types of time-domain results were calculated: simulation results and prediction results. Simulating a model means that the response of a model to a particular input is computed. On the other hand, predicting future outputs of a model from previous data over a time horizon of ksamples or kT time units, requires both past inputs and past outputs. The main difference between simulation and prediction is whether it is used measured or computed previous outputs for calculating the next output. A high simulation performance indicates that a model is a good representation of reality. However, the model purpose should also be considered when testing a model. Using a model for prediction is common in controls applications like predictive control, where it is wanted to predict output for a specific number of steps in advance. This type of control is going to be used later in this thesis, and therefore the relevance of these prediction results. Both, simulation and prediction results using the identified ARX model were compared to measured output-data. Model inputs (measured discharges) are presented in figure 5.32 while the results are presented in figure 5.33 and 5.34 for simulation and prediction respectively. 5.5. Experimental application 123 0 200 400 600 800 1000 −20 0 20 Disch. var., q1 (L/s) Gate 1 0 200 400 600 800 1000 −20 0 20 Disch. var, q2 (L/s) Gate 3 0 200 400 600 800 1000 −5 0 5 Disch. var., q3 (L/s) Gate 5 0 200 400 600 800 1000 −1 0 1 Disch. var., qL 1 (L/s) Weir 1 0 200 400 600 800 1000 −20 0 20 Time, t (s) Disch. var., qL 2 (L/s) Weir 3 Figure 5.32: Discharge sequences used to feed the ARX model 124 Chapter 5. Canal Identification for Control Purposes 0 200 400 600 800 1000 −15 −10 −5 0Pool 1 Level var., zs 1 (cm) 0 200 400 600 800 1000 −10 −5 0 5 Pool 2 Level var., zs 2 (cm) 0 200 400 600 800 1000 −10 −5 0 5 Pool 3 Time, t (s) Level var., zs 3 (cm) Measured Simulated Measured Simulated Measured Simulated Figure 5.33: Comparison between measured water level deviations and model simulation 5.5. Experimental application 125 0 200 400 600 800 1000 −15 −10 −5 0 Level var., zs 1 (cm) Pool 1 0 200 400 600 800 1000 −10 −5 0 5 Level var., zs 2 (cm) Pool 2 0 200 400 600 800 1000 −10 −5 0 5 Time, t (s) Level var., zs 3 (cm) Measured Predicted Measured Predicted Measured Predicted Figure 5.34: Comparison between measured water level deviations and a 20-step-ahead (200 s) model prediction The simulation results of figure 5.33 show that the accuracy of the identified model deteriorates when the water level deviations exceed approximately the 5 cm value, a 10 % of the initial water level. Thus, the model gives good estimations within a ±10 % operation area. This result is somewhat expected for a linearized model approximating a nonlinear process. Within this ±10 % range, delays and transients are well approximated and there is a good agreement between measured and simulated water levels. The prediction results of figure 5.34 look even better. The 20-step-ahead predictions (a 200 s forecast) are very accurate when compared to the measured water levels, inclusively beyond the ±10 % operation range. The maximum residuals are only slightly higher than 1 cm. In summary, the model has successfully approximated the linearized behavior of the canal around an operation condition and is specially suitable for medium-term predictions. A model like this is accurate enough to be used in the design of any linear controller, but is particularly attractive for predictive controllers. 126 Chapter 5. Canal Identification for Control Purposes 5.5.7.2 Frequency domain In this case, the frequency response of the model was checked against the theoretical Bode diagram obtained from an unidimensional Saint-Venant canal model around an operation point. This model was configured using the geometrical dimensions of the real canal, but its roughness was not measured. Hence, the Manning number was arbitrarily taken as 0.02. The comparison can be observed pool by pool in figures 5.35, 5.36 and 5.36 for Pool 1, Pool 2 and Pool 3 respectively. 10−3 10−2 10−1 −5 0 5 10 15 20 25 30 35 Frequency (rad/s) Magnitude (dB) Theor. ARX (a) Bode Magnitude q1→zs1 10−3 10−2 10−1 −10 0 10 20 30 Frequency (rad/s) Magnitude (dB) Theor. ARX (b) Bode Magnitude q2→zs1 10−3 10−2 10−1 −700 −600 −500 −400 −300 −200 −100 0 Frequency (rad/s) Phase (deg) Theor. ARX (c) Bode Phase q1→zs1 10−3 10−2 10−1 −50 0 50 100 150 200 250 Frequency (rad/s) Phase (deg) Theor. ARX (d) Bode Phase q2→zs1 Figure 5.35: Theoretical Bode diagram v/s Bode diagram of identified ARX model - Pool 1 5.5. Experimental application 127 10−3 10−2 10−1 −5 0 5 10 15 20 25 30 Frequency (rad/s) Magnitude (dB) ARX Theor. (a) Bode Magnitude q2→zs2 10−3 10−2 10−1 −20 −10 0 10 20 30 Frequency (rad/s) Magnitude (dB) Theor. ARX (b) Bode Magnitude q3→zs2 10−3 10−2 10−1 −800 −700 −600 −500 −400 −300 −200 −100 0 Frequency (rad/s) Phase (deg) Theor. ARX (c) Bode Phase q2→zs2 10−3 10−2 10−1 −50 0 50 100 150 200 250 Frequency (rad/s) Phase (deg) Theor. ARX (d) Bode Phase q3→zs2 Figure 5.36: Theoretical Bode diagram v/s Bode diagram of identified ARX model - Pool 2 128 Chapter 5. Canal Identification for Control Purposes 10−3 10−2 10−1 −5 0 5 10 15 Frequency (rad/s) Magnitude (dB) Theor. ARX (a) Bode Magnitude q3→zs3 10−3 10−2 10−1 −450 −400 −350 −300 −250 −200 −150 −100 −50 0 Frequency (rad/s) Phase (deg) Theor. ARX (b) Bode Magnitude q3→zs3 Figure 5.37: Theoretical Bode diagram v/s Bode diagram of identified ARX model - Pool 3 These figures show that the frequency responses of both models are in good agreement. The analysis is detailed in the following points: •The magnitude plots display very small difference in the low frequency area. •The first resonance peaks are located at identical frequencies in the magnitude plots, but their gains differ slightly, specially close to the Nyquist frequency. •The concordance of the phase plots exhibits that the delays predicted by both models are almost identical. The great similarity between the frequency responses of both models demonstrate that the ARX model does not only produce good results in the studied cases; it actually captured the behavior of the process around a working point. In fact, it is very likely that the magnitude differences correspond to the influence of dynamics that the unidimensional theoretical model has not considered and that the ARX model identified from the experimental data. 5.6. Conclusions 129 5.6 Conclusions The main conclusion of this chapter is that an irrigation canal can be well modeled for control purposes by linear black-box models obtained by means of system identification techniques. Nevertheless, it is very important to take into account the "special characteristics" of the system (delays, resonant modes, integrator dynamics, etc.) when designing the identification procedure. If the process is performed in a blind manner, identified models are likely to be inaccurate and/or unappropriate and/or unstable. Getting more in detail, the following concepts arose from this chapter: •A detailed mathematical analysis has concluded that there is a general linear model structure that is applicable to any type of canal. One of the most relevant outcomes from this analysis, is that there are some new mathematical approximation results which reinforce and support the idea already suggested by some researchers, that there is always an integrator pole within this general structure. •It is extremely important to take into account the structure, properties and characteristics of a pool before proceeding with system identification. This information is crucial to adequately select issues like: sampling time, model structure, experiment type, etc.. •A linearized pool model is only valid around a particular operation point. Among operation points, there are changes in the model gain, the time delay and in the amplitude and location of resonant modes. •Long slopping pools are likely to have big transport delays. Short flat pools have also transport delays and in addition resonant modes. Both of them act as a swimming-pool or water tank in the long term behavior. •It is recommendable to work with discharges and water levels in a pool identification procedure. Gate openings will always add additional nonlinearities and augment the identification complexity. In fact, the validity range of a pool linearization is generally shorter when using gate openings instead of gate discharges. Gates and other hydraulic structures can always be afterwards included in an already identified discharge-based model. •There is not a unique solution in the sampling time selection. However, it was found useful to obtain the step response of a pool and divide approximately the time taken to reach a clear tendency by 15. •PRBS type discharge signals are very appropriate to collect pool identification data.