Irrigation canal system identification for control purposes

Irrigation canal system identification for control purposes Carlos Sepúlveda Toepfer, José Rodellar Benedé IOC-DT-P-2005-19 Setembre 2005 Abstract This report gives some guidance on how to obtain linear black-box models of irrigation canal reaches, using system identification techniques. First of all, some general properties of the irrigation canal reaches are deducted, based on the use of the linearized Saint-Venant equations to model the water behavior. Then different aspects of the system identification procedure like the sampling time, the model structure, the experiment design, etc., are studied, in order to avoid possible modelling problems and, in that manner, obtain a good linear model capable to be used in control systems designs. The results obtained in the time domain and in the frequency domain show that one can achieve very accurate models, if the system identification procedure is designed with care having in mind the intrinsic properties of the system. The research reveals that it is not convenient to perform a black-box irrigation canal system identification, without having a certain knowledge of the system. Contents 1 Introduction 1 1.1 About the importance of models in automatic control . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Aboutirrigationcanalmodels .................................... 2 2 Model of a reach 3 2.1 Mathematicalmodel ......................................... 3 2.2 Linearizationofthemodel ...................................... 5 2.3 Properties of the linearized model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3 System identification of the model 9 3.1 Discrete-time modelling issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.1 Samplingtime ........................................ 9 3.1.2 Discrete transfer functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2 Discrete-time model structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.1 ARXmodel.......................................... 14 3.2.2 Laguerremodel........................................ 15 3.3 Experiment design: Input signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.4 Parametricidentification ....................................... 19 3.4.1 ARXmodel.......................................... 19 3.4.2 Laguerremodel........................................ 20 4 Results 22 5 Conclusions 29 Chapter 1 Introduction Water irrigation canals are systems developed to transport water from main water storage reservoirs to several agricultural water-demanding farms, in irrigational seasons. Generally they cover very long distances; their range can go from hundreds of meters to hundred of kilometers, where farms are located normally close to them, but all the way along. To have a certain degree of command of the system, they normally have some control structures, like check gates, embedded in the water path, so as to regulate the amount of water flow in concordance with the water demands. With respect to the farmer’s water offtakes, they are generally situated a few meters away from the upstream side of the gates, and the extractions are performed by pumps, weirs or any appropriate device. As can be thought, it is not a trivial matter to manage this type of systems. They have to transport the water, minimizing the losses, and ensuring that each local farmer becomes their corresponding (exact when possible) amount of water at their corresponding frequency. Besides the inherent characteristics of the system don’t help to much in reaching the objectives. The system presents very long time lags (from minutes to hours) in the transport of the water, delaying every decision applied on the system. Moreover, there are important water dynamical effects (generated, normally, by any change in the amount of water delivered), that produce, in different degrees depending on each case, interferences in the deliveries from the whole system (coupling). Historically, these problems and the availability of water have motivated the creation, in many countries, of irrigation associations with their own irrigation statutes and rules. For these reasons many researchers have paid attention in improving the system’s operational management by the automation of the system, applying the control engineering theory tools. As commented in many papers about the matter (Clemmens et al., 1998; Malaterre et al., 1998; Sawadogo et al., 2000; G´ omez et al., 2002; Rodellar et al., 2003), the goal of the automatic control of the operation is to maintain the water depth levels in the extraction zones, as constant as possible, by moving the intermediate check gates. This goal can be explained by the following reason: either if the irrigation water is taken out of the system by pumps or weirs, a constant level assures a constant supply of the water, eliminating flow variations to the farmers, and in that way, eliminating the coupling effect produced. In that way, several water demands can be fulfilled minimizing the interference between them. 1.1 About the importance of models in automatic control The use of a good model of a process to be controlled is indispensable for almost all the currently existent control techniques (Shook et al., 1992). In addition to the possibility to test the control strategy to be implemented by means of computational simulation of the model, many control techniques use the model explicitly in the design 2 Chapter 1. Introduction stage of the controller or/and in the calculation of the control action (in this case check gates movements). In the situation where the model is used in the design stage or in the computation of the control action, unfortunately, the complexity of the characteristics of the model (nonlinearities, delays, instability, etc), is, usually, directly proportional to the complexity of the control techniques and their implementation. For example, a linear and rational model opens the possibility to apply numerous and well-known linear control theories and standard techniques, which are relatively easy to implement. A nonlinear model, in contrast, requires much more effort and time to solve the control problem. Therefore it is always desirable to have the simplest model that can reproduce the behavior of the system. However, the ”closest to exact” behavior of a system is always obtained, when it can be obtained, by very complex models. For example all the existing processes are actually nonlinear. Linearity is only a simplification of the problem. So, models for control purposes should be made as a trade-off between simplicity and accuracy of the model, in concordance with the automation goal. 1.2 About irrigation canal models Models that involve water are generally obtained making use of simplifications of the Navier-Stokes Equations, because of the complexity in dealing with them directly. For example for irrigation canals, one of the most accepted and used model in simulations, is the system given by the Saint-Venant Equations (Henderson, 1966), because of its capacity to represent the characteristics of real interest. However this system is a nonlinear partial differential equation system, which has analytical solution only in very special cases, having to utilize numerical methods to solve it properly. As model for computational simulation it is very accurate, but as model for control, it is clearly not appropriate for the reasons exposed before. That is the reason why, usually, linearizations or simplifications of the Saint-Venant equations are recurrently studied by the irrigation control research community (Schuurmans et al., 1995, 1999; Litrico and Fromion, 2004; Weyer, 2001). These models normally are physical models. That means that they are based on physical parameters of the real system. This has advantages and drawbacks. The advantages are that a physical model is reliable and has a very strong connection with the theoretical concepts. The main drawback is that it is necessary to feed the model with many parameters that have to be determined (from theoretical and/or experimental results) and an error in this determination can lead to a wrong modelling. For example, a bad determination or small variation of the Manning number (a coefficient related to the friction of the canals) can produce quite different results. The scope of this work is to study the application of black box models to the problem of modelling irrigation canals for control purposes. These models do not necessary have a physical meaning; they only focus on reproducing the behavior of a system, and are obtained from input-output experimental collected data, in a process that is called System Identification. Chapter 2 Model of a reach 2.1 Mathematical model i-1 ii+1 L i L i-1 L i+1 s i-1 s i s i+1 Figure 2.1: Irrigation canal schematic A simplified vision of a typical irrigation canal can be observed in figure 2.1. As mentioned, it receives water form a source, and lets the water flow with a small slope. The intermediate check gates, represented by vertical lines, regulate with their openings (Wi) the desired flux, so as to maintain the water depth level (Zs i) in the zones, where a water flow (QL i) is extracted for irrigational purposes. For modelling intentions, a natural way of partitioning a canal is dividing it into reaches (also called pools). A reach is a portion of a canal between two check gates. So a normal canal can have several reaches with different characteristics (length, slope, width, etc.). However all the reaches share a common structure, focusing the problem of modelling an irrigation canal, in finding a suitable model for reaches. In that manner, the problem can be solved by a sum of the same model structure with only different parameter values. In the specialized literature there are several approaches to obtain this model. In this work two facts will be taken into account: the location of the water extractions, generally near the reach’s end, and the behavior’s gradual change of the water when approaching to an obstacle like a cross-gate, resembling the characteristics of a water deposit. To put this practical knowledge in a mathematical model, it can be created an imaginary bound, which separates the reach, in an absolute manner, in a water transport area and in a water storage area (see figure 2.2). It should be noted that the location of the imaginary division is not so crucial; as it defines a particular size of the storage area, it will be correctly chosen whenever that area is remained small enough in relation to the reach’s total size. 4 Chapter 2. Model of a reach Storage Area Transport Area Gate iGate i+1 QiQi+1 QL i Zs i Qi 0 m X m X m Figure 2.2: Simplified representation of a reach Making use of the Saint-Venant equations (Henderson, 1966) to model the transport, and of the mass conservation principle to model the storage, that leads the following mathematical model: •Transport ∂A ∂t +∂Q ∂x = 0 (2.1) ∂Q ∂t +∂ ∂x Q2 A+gA∂Z ∂x =gA(S0−Sf)(2.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 (2.3) In the transport equations (2.1)-(2.2), xis the longitudinal coordinate in the flow direction, tis the time, A=A(x, t)is the wetted transversal section area of the canal, Q=Q(x, t)is the volumetric water discharge, Z=Z(x, t)is the water depth, S0is the bottom slope, Sf=Q2n2 A2(A/P )4 3 is the friction slope of the canal and gis the gravity acceleration. Besides, it should be noted that the wetted transversal section area (A) depends explicitly on the water depth (Z) and on the transversal section shape of the canal. In the storage equation (2.3), Qi(t)is the water flow that enters the area, Qi+1(t)is the water flow delivered to the next reach, QL i(t)is the water flow 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. As can be seen from figure 2.2 and equations (2.1), (2.2) and (2.3), the variable values 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 is solvable. In summary, the model consist in a system of two nonlinear PDEs and one nonlinear (sometimes linear) ODE. Because of the reasons given before, the model is little advantageous for control purposes. In the search for a more convenient model, a model linearization around an operational condition should be performed. 2.2. Linearization of the model 5 2.2 Linearization of the model Linearization is carried out replacing in the model, described by (2.1), (2.2) and (2.3), the expressions of the variables around a working point, namely, Z(x, t) = Z0+z(x, t)and Q(x, t) = Q0+q(x, t), and neglecting all the second-order terms (Litrico and Fromion, 2004). For the Saint-Venant equations (2.1) and (2.2), that yields: B0 ∂z ∂t +∂q ∂x = 0 (2.4) ∂q ∂t + 2V0 ∂q ∂x −β0q+C02−V02B0 ∂z ∂x −γ0z= 0 (2.5) with γ0=V02dB0 dx +gB0(1 + κ)S0−1 + κ−F02(κ−2)∂Z0 ∂x ,β0=−2g V0S0−∂Z0 ∂x and κ=7 3− 4S0 3B0P0 ∂P0 ∂Z . In (2.4) and (2.5), F0=V0 C0is the Froude number, P0is the wetted perimeter, C0=qgA0 B0is the water acceleration and V0=Q0 A0is the water speed; all of them evaluated at the operation condition. The system, in this manner, has boundary conditions: q(0, t) = qi(t)q(X, t) = qi(t), and z(0, t) = zi(t)z(X, t) = zi(t). One way of obtaining a solution for this system, is applying the Laplace transform and then reordering. That produces the following system of ordinary differential equations in the variable x, with a complex parameter s(the Laplace variable): d dx q(x, s) z(x, s)=A(x, s)q(x, s) z(x, s)1(2.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   . 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 x, has 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 (2004) that the problem can be solved numerically very efficiently, if it can be discretized by several ”mini uniform regimes problems”, in a way like in figure 2.3. Using this approach, the solution to the linearized Saint-Venant equation system 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)(2.7) where the transfer function matrix Γshould be calculated using 1In the following any f(s)will correspond to L {f(t)}, the laplace transform of f(t), which is a complex-valued function. 6 Chapter 2. Model of a reach Z(x) x X1X2X3 a b c d X0X4 Figure 2.3: Schematic representation of an approximation by uniform regimes Γxn=X, x0= 0= 0 Y k=n−1 eA(xk, s)hk(2.8) with hk= (xk+1 −xk). It should be noted that the term eA(xk,s)hkcorrespond to the exponentiation 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)     (2.9) where λ1,2(s) = 1 fas +b±pcs2+ds +e 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 (2.7) using (2.8) and (2.9), sis the laplace variable, Xis the xposition of the end of the transport area (beginning of the storage area) and Γis the transfer function matrix that describes exactly all the input-output relationships of the linearized Saint-Venant equations in the laplace domain. In fact, it is possible 3.1. Discrete-time modelling issues 13 where zis the z-transform complex variable. For this reason it maintains a close relationship with the continuous Laplace transfer function obtained in (2.13). In that transfer function model, although it is general for any type of reach, there are two things that will always 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 that, 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 1 z−1will always appear in the z-transform of the model. 2. The time delays between the inputs (water flows) and the output (water depth of interest) 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 z is the z-transform variable. This implies that the z-transform model (3.1) will always appear in the following manner: zs i(z) = z−d1F0 1(z)qi(z)−z−d2F0 2(z)qi+1(z)−z−d2F0 2(z)qL i(z)(3.2) where d1and d2are the time periods that takes each water flow to influence the water depth at the extraction zone, measured in discrete time instants. Point 1 presents a property that have to be taken with care in order to avoid possible modelling 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 the position will lead to a completely different dynamical model’s behavior. With a pole at z= 1.01 the model would be unstable, making the model output go to infinite values after a small time period; with a pole at z= 0.99 the model would be strictly stable, producing that, with a constant input, the model output always reaches a constant value. There are three approaches normally taken with respect to this problem when estimating a model by means of system identification: •To forget about the problem and obtain a model anyway. •To identify a model and afterwards correct the position of the pole in the estimated model. •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 the following way: 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) (3.3) The first option is generally more recommendable (especially in the presence of noise) and is equivalent to make a cumulative sum of the input data, and the second option correspond to a differentiation of the output data. However it is important to note, especially when a maximally informative input signal was especially designed to identify a model, that the first option modifies the frequency characteristic of the input signal. 14 Chapter 3. System identification of the model From the three presented approaches, it is clear that the last two are the better ones. The third one, however, has the disadvantage when modelling reaches, that because the sample time must be very small in order to reproduce the fastest dynamics adequately, larger sample periods could induce long-term modelling bias errors. So, in this case, the second option of identify and then correct would be the best when it is possible to perform. This is not always easily possible for some discrete-time model structures and, moreover, there are some types of model structures that cannot deal with processes with an integrator. For those cases the only choice is the third approach. 3.2 Discrete-time model structures Usually linear discrete-time models can be divided into three main different classes: 1. Discrete transfer function models based only on the input-output characteristics of a system by means of the z-transform. 2. Discrete state-space models that incorporate the information related to all internal dynamics of a system directly in the time domain. 3. Orthonormal basis models that make use of the special properties of some basis functions to approximate systems. Depending on the class and particular structure chosen, there is a whole list of parameter estimation methods that, with a given informative data set, can estimate the parameter values of that model in order to approximate a true system as well as possible. 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.. In this work, only two model structures for reach modelling are considered: the ARX (Autoregressive with Exogenous Input) model, a transfer function based model, and the Laguerre model, an orthonormal basis based model. Both models have a finite number of rational elements, although the process is governed by an irrational transfer function. An approximation like that can be carried out because generally an irrational term can be approximated by a linear combination of rational ones, like in a Pad` e approximation of a function. However, depending on the irrational term, a good approximation could require a great number of rational terms to achieve good results. All the model parameters are going to be estimated by a least-squares based algorithm in order to obtain a discrete-time approximation model of a reach. 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 by: qf(k) = f(k+ 1), q−1f(k) = f(k−1) Then assuming that F0 1(z)and F0 2(z)of (3.2) can be approximated by two rational quotients of polynomials in the following manner: F0 1(z) = B1(z) A(z), F0 2(z) = B2(z) A(z) 3.2. Discrete-time model structures 15 with A(z) = 1 + a1z−1+. . . +anaz−na B1(z) = b11+b12z−1+. . . +b1nb1z−nb1+1 B2(z) = b21+b22z−1+. . . +b2nb2z−nb2+1 and replacing in (3.2), that yields after reordering: A(z)zs i(z) = z−d1B1(z)qi(z)−z−d2B2(z)qi+1(z)−z−d2B2(z)qL i(z) Finally applying the inverse z-transform to move the problem backwards to the time domain, the general model structure for this case is: A(q)zs i(k) = B1(q)qi(k−d1)−B2(q)qi+1(k−d2)−B2(q)qL i(k−d2)(3.4) with A(q) = 1 + a1q−1+. . . +anaq−na B1(q) = b11+b12q−1+. . . +b1nb1q−nb1+1 B2(q) = b21+b22q−1+. . . +b2nb2q−nb2+1 As can be seen, the parameters of the ARX model (3.4) are the polynomial orders na,nb1and nb2, the polynomial coefficients, and the a priori known time delays d1and d2expressed in sampling instants. Once the parameters are determined, this model can calculate the downstream water level zs i at instant k, by a weighted sum of past values of the water level zs i and of past values of the water discharges qi,qi+1,qL i, all of them previously collected at a given sampling interval T. In some sense it is a little bit restrictive that both transfer functions, B1(z) A(z)and B2(z) A(z), have the same denominator polynomial A(z), specially in this case when by means of (2.13) it is clear that they really don’t have the same denominator. However (2.13) shows, in the same manner, that both transfer functions share some poles, so 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 manner it can be performed with a linear least squares approach. This model is capable to model stable and unstable processes so, as mentioned before, there exist 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 obtaining the roots of A(q), then replacing the root of A(q)that is near to 1by exactly a 1and afterwards forming a new A(q)polynomial keeping the other roots locations. The second option is to estimate the model parameters using the cumulative sum of each input model variable (water flows) and then multiply the estimated A(q)polynomial by (1 −q−1). 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 16 Chapter 3. System identification of the model 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 cT=c1c2··· cN lT=l1(t)l2(t)··· lN(t) LT=L1(s)L2(s)··· LN(s) 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)(3.5) where l(k)is the state vector of order N,u(k)is the system input and y(k)is the system output. Moreover, if Tis the discrete sampling time, Aand Bare 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),τ4=√2p(1−τ1) p. Making use of the Laguerre model structure for a reach, the particular model gives: 3.3. Experiment design: Input signal 17 zs i(k) = c1T[A1 l1(k−1) + B1 qi(k−1)] (3.6) −c2T[A2 l2(k−1) + B2 qi+1(k−1)] −c2T[A2 l2(k−1) + B2 qL i(k−1)] In this type of model, the model output is obtained with only the inputs and the Laguerre functions values (always known for a particular p) at time k−1. The parameters of the Laguerre-based model (3.6) are the coefficients vectors c1 and c2, and the designerchosen number of Laguerre functions Nand Laguerre pole (p) value (to calculate Ax and Bx), one for each transfer function to be approximated. It is not necessary an a priori knowledge on the process orders or on any delay, but for increasing complexity transfer functions, it is necessary a higher number of terms in order to approximate adequately the required behavior. However the dynamical responses of each input can be adjusted in a totally independent manner. This model can not approximate systems that are not strictly stable, so in this case the integrator should be passed directly to the data before proceeding with the system identification procedure. That can be realized through (3.3). Afterwards the integrator must be included in the obtained process model. 3.3 Experiment design: Input signal An input signal should excite a system and have a rich frequency content, in order to successfully identify a model that approximates a true real system. A rich frequency content means that the input signal contains sufficiently many distinct frequencies. In practise, it is suitable to decide upon an important and interesting frequency band to identify the system in question, and then select a signal with more or less flat spectrum over this band. Such an input is provided be 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 larger time horizon. 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(3.7) where Nis the number of stages of the shift register. An example of a portion of a PRBS signal is presented in figure 3.3. Assuming that u(k)is a random binary process with the 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 (3.8) The spectral density function (3.8), 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 broad-band noise. 18 Chapter 3. System identification of the model 0 50 100 150 200 −1.5 −1 −0.5 0 0.5 1 1.5 Time [samples] Magnitude Figure 3.3: PRBS Test Signal The spectrum of the pseudo-random binary signal is therefore an approximation for 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: •One should work with an integer number of periods; that means: length{u(k)}=nL, with n={1,2, . . .}. Usually one period of the signal meets the needs of the identification procedure. •The amplitude of the signal a, must assure a good signal-to-noise ratio; normally ten times greater than the noise amplitude. Additionally in order to correctly identify the steady-state gain of the process, the duration of at last one of the pulses in the PRBS must be greater than stabilization time of the process. As the maximum duration of a pulse is NTprbs,Nand Tprbs have to be chosen to adequately cover that time period. As can be seen, if the PRBS sampling period Tprbs is chosen equal to a well-chosen process sampling period Ts (see section 3.1.1), the only way to augment the maximum duration of a pulse is increasing the number of registers N. However higher values of Nlead to very large signal sequences, increasing the experiment duration and the data length; e.g. N= 10 produces a sequence length of 1023 data points, when N= 12 produces a sequence length of 4095 data points. In that case it is useful to choose the PRBS sampling period Tprbs to be a multiple of the process sampling period Ts:Tprbs =pTs. The problem of this approach is that reduces the frequency range corresponding to a constant spectral density, so usually pis chosen to be p≤4. From a system identification point of view, the data length ranges normally from 200 to 1000 data points, in order to have reliable values of the model parameters and less computational burden. On the other hand, the proposed models are supposed to work only around an operational point of the reach, so to have locally rich data, all the movements induced to the system should be maintained sufficiently small. For normal systems this can be accomplished choosing a small enough PRBS amplitude a, but in this case that the system behaves like an integrator, it is additionally convenient that the experiment duration is kept as short as possible, to avoid that the system goes to far away from the working point. 3.4. Parametric identification 19 3.4 Parametric identification System identification can be defined as the process of obtaining a model for the behavior of a plant, based on the plant input and output data. If a particular model structure is assumed, the identification problem is reduced to obtaining the parameters of the model. The usual way of obtaining the parameters of the model is optimizing a function that measures how well the model, with a particular set of parameters, fits the existing input-output data. When process variables are perturbed by noise of a stochastic nature, the identification problem is usually interpreted as a parameter estimation problem. This problem has been extensively studied in literature for the case of processes which are linear on the parameters to be estimated and perturbed with a white noise (Ljung, 1999; Camacho and Bordons, 2004). That is, processes that can be described by: zk=Θ Φk+ek(3.9) where Θis the vector of parameters to be estimated, Φkis the vector of past input and output measures, zkis the latest output measure and ekis a white noise. Once a model is written in a form like in (3.9), the parameters can be identified by using a least-squares identification algorithm. All the models proposed in this work can easily be expressed as in (3.9) as follows: 3.4.1 ARX model Assuming that the disturbances, that is, the differences between the measured output and the output calculated by the model, can be described by ei(k), a white noise zero mean sequence, model equation (3.4) can be rewritten as: A(q)zs i(k) = B1(q)qi(k−d1)−B2(q) [qi+1(k−d2) + qL i(k−d2)] + ei(k)(3.10) Then solving for zs i(k), (3.10) yields: zs i(k) = A0(q)zs i(k−1) + B1(q)qi(k−d1)−B2(q) [qi+1(k−d2) + qL i(k−d2)] + ei(k)(3.11) with A0(q) = (1 −A(q)) q=−a1−a2q−1+. . . +anaq−na−1 This can be expressed as (3.9), by making 20 Chapter 3. System identification of the model zk=zs i(k) Θ=a1a2··· ana b11b12··· b1nb1b21b22··· b2nb2 Φk=                        −zs i(k−1) −zs i(k−2) . . . −zs i(k−na) qi(k−d1) qi(k−d1−1) . . . qi(k−d1−nb1 + 1) −qi+1(k−d2)−qL i(k−d2) −qi+1(k−d2−1) −qL i(k−d2−1)) . . . −qi+1(k−d2−nb2 + 1) −qL i(k−d2−nb2 + 1)                        3.4.2 Laguerre model Assuming that the disturbances of the model can be described in the same manner as before, that is by ei(k), a white noise zero mean sequence, model equation (3.6) can be rewritten as: zs i(k) = c1T[A1 l1(k−1) + B1 qi(k−1)] (3.12) −c2T[A2 l2(k−1) + B2 (qi+1(k−1) −qL i(k−1))] + ei(k) Remembering that l(k) = A l(k−1) + Bu(k−1) (3.13) it is straightforward to show that (3.12) can be expressed as (3.9) in the following manner: zk=zs i(k) Θ=c11c12··· c1Nc21c22··· c2N Φk=               l11(k) l12(k) . . . l1N(k) −l21(k) −l22(k) . . . −l2N(k)               As can be seen, in order to perform the parameter estimation procedure, first of all it is necessary to calculate 3.4. Parametric identification 21 the values of the NLaguerre functions responses at instant kfor each model input, by means of the function (3.13). In that case the initial state vector should be l(0)T=0 0 ··· 0. Chapter 4 Results To prove the effectiveness of the system identification of a reach, different models were identified for the reach with characteristics presented in table 3.1, by means of the parameter estimation of the proposed reach-model structures. The initial water profile of the reach can be observed in figure 4.1. Figure 4.1: Initial water profile of the reach Then following the recommendations made in section 3.1.1, the sampling period was chosen to be 215 s, and an informative-rich input sequence was designed in order to excite all the relevant dynamics of the system. In this manner, two maximum length PRBS signals were designed in a form that the maximum width of the PRBSs was greater than the tendency-stabilization time of the step response of the reach: 4000 s −680 s = 3320 s. Besides having in mind ”a not too large” sequence length, the following PRBS parameters were chosen: p= 2 and N= 8. Therefore the total sequence length was: p×L= 2 ×28−1= 510 data points Chapter 5 Conclusions Considering all the results presented before, the main conclusion of this work is that the linearized behavior of a canal reach can be approximated in good manner for control purposes by linear black-box models obtained by means of system identification techniques, in spite of the ”special” characteristics of the system in question, if some appropriate designing guides related to the system are followed. However, it is very important to remark that if the system identification procedure is performed in a blind manner, the obtained models could be very inaccurate and/or unappropriate and/or unstable. This is more true when the integrator presence is forgotten or when the sampling time is badly chosen. Getting more in detail, there are some ideas crucially related to the problem in question. It has been proven that following some vastly extended modelling assumptions, the relationship between the downstream water level and the input and output water discharges has some common elements, that are independent of the reach configuration. These elements, observed under the frequency domain approach, are a particular basic structure, the irrationality of their components and an inherent integral behavior of the system. Having in mind these particular characteristics, some of the system identification techniques were revised and adapted to the reach modelling problem. The first problem to cope is the election of the sampling period. It is not a trivial matter and depends on a large degree on each particular case. It has been found that to perform a step response of the reach is useful to have an insight of the dynamical evolution of the water level, in order to, at least, keep track of the major tendency changes of the system. However this is only a point of view, and other interests, particulary those related to the control algorithm stability and complexity, can be more preponderant. The identification experiment is one of the most important parts of the procedure, because it has to induce the process to show all its main dynamical characteristics in the generated data set. The information content has to be a priority, but without making the system go away to far away from the operating condition, because then the linearity assumption is not valid. In this case, because of the integrating property of a reach, the test should be maintained as short as possible, because of the probability of leaving the operating region. A pseudo random binary sequence (PRBS) is generally a good option, as demonstrated in this work. However, depending on the particular reach configuration or depending on the structural constrains of the water management devices, another type of signals could be more appropriate. Actually normal operation data of an irrigation reach can also be used, if it reaches a minimum threshold of information quality. For all the models used in this work, the parameters can be easily obtained by linear least squares parameter estimation techniques. In fact they can easily be estimated online; that means that they can be estimated when the system is actually working, and, moreover, the algorithms can be programmed to perform adjustments in the models if the conditions of the system change, as well as in adaptive control strategies. 30 Chapter 5. Conclusions From the proposed models, the ARX models perform very well and better than the Laguerre model, despite of their inherent restrictions (linear model, same denominator polynomial, rational expressions, etc). Their structure (delays and orders) had to be chosen comparing the performance of different choices, but their simplicity makes the process very easy to accomplish. The results didn’t show a real difference between the two ARX developed variants, but perhaps it is better to identify the model and afterwards correct the integrator position, because the process is not difficult to implement and then the data set is not filtered in any way. The Laguerre model needed more parameters to cope with the ARX models, but it is a good alternative, specially because it doesn’t require the knowledge of the delays or of the appropriate system orders; only the Laguerre poles have to be chosen. It directly approximates the system, and can give totally independent responses from each water flow discharges when required. All the presented reach models can easily be used to generate a whole irrigation canal model. In that case, a multi-reach canal would be seen as a multiple input - multiple output (MIMO) system. So, the presented work is not only useful for decentralized reach controller designs: it can be still used when developing a centralized controller for an entire canal. Bibliography J. T. Bialasiewicz. Advanced system identification techniques for wind turbine structures with special emphasis on modal parameters. NASA STI/Recon Technical Report N, 96:11276–+, 1995. E. F. Camacho and C. Bordons. Model Predictive Control. Springer-Verlag, London, second edition, 2004. A. J. Clemmens, T. F. Kacerek, B. Grawitz, and W. Schuurmans. Test cases for canal control algorithms. Journal of Irrigation and Drainage Engineering, 124(1):23–30, 1998. M. G´ omez, J. Rodellar, and J. A. Mantec´ on. Predictive control method for decentralized operation of irrigation canals. Applied Mathematical Modelling, 26(11):1039–1056, 2002. F. M. Henderson. Open channel flow. MacMillan Publishing Co., Inc., New York, 1966. X. Litrico and V. Fromion. Frequency modeling of open channel flow. Journal of Hydraulic Engineering, 130(8): 806–815, 2004. L. Ljung. System identification. Theory for the user. Prentice-Hall, Inc., Upper Saddle River, New Jersey, second edition, 1999. P.-O. Malaterre, D. C. Rogers, and J. Schuurmans. Classification of canal control algorithms. Journal of Irrigation and Drainage Engineering, 124(1):3–10, 1998. R. Malti, S. B. Ekongolo, and J. Ragot. Dynamic SISO and MIMO system approximations based on optimal Laguerre models. IEEE transactions on Automatic Control, 43(9):1318–1323, 1998. J. Rodellar, C. Sep´ ulveda, D. Sbarbaro, and M. G´ omez. Constrained predictive control of irrigation canals. In Proceedings of the 2nd International Conference on Irrigation and Drainage, pages 477–486, Phoenix, Arizona, may 2003. USCID. S. Sawadogo, R. M. Faye, A. Benhammou, and K. Akouz. Decentralized adaptive predictive control of multi-reach irrigation canal. In 2000 IEEE International Conference on Systems, Man, and Cybernetics, pages 3437–3442, Nashville, Tennessee, oct 2000. J. Schuurmans, O. H. Bosgra, and R. Brouwer. Open-channel flow model approximation for controller design. Applied Mathematical Modelling, 19(9):525–530, 1995. J. Schuurmans, A. J. Clemmens, S. Dijkstra, A. Hof, and R. Brouwer. Modeling of irrigation and drainage canals for controller design. Journal of Irrigation and Drainage Engineering, 125(6):338–344, 1999. D. S. Shook, C. Mohtadi, and S. L. Shah;. A control-relevant identification strategy for GPC. IEEE Transactions on Automatic Control, 37(7):975–980, 1992. E. Weyer. System identification of an open water channel. Control Engineering Practice, 9(12):1289–1299, 2001. C. C. Zervos and G. A. Dumont. Deterministic adaptive control based on Laguerre series representation. International Journal of Control, 48(1):2333–2359, 1988.