- Email: [email protected]
Analytical approximation of open-channel flow for controller design
Applied Mathematical Modelling 28 (2004) 677–695 www.elsevier.com/locate/apm
Analytical approximation of open-channel flow for controller design Xavier Litrico a
a,*
, Vincent Fromion
b,1
Cemagref, UR Irrigation, 361 rue-J-F Breton, B.P. 5095, 34033 Montpellier Cedex 1, France b INRA, LASB, 2 place Viala, 34060 Montpellier, France
Received 13 March 2003; received in revised form 2 October 2003; accepted 27 October 2003 Available online 21 January 2004
Abstract The integrator delay (ID) model is a popular and simple way to model a canal for control purposes. In this paper, particular attention is given to accurately model the delay and the integrator gain in any flow configuration. Based on a previous work by the authors allowing to have a very accurate frequency domain representation of Saint-Venant equations, the paper proposes a new approximate model: the integrator delay zero (IDZ) model has an integrator and a delay in low frequencies, and models the high frequencies by a constant gain and a delay. Analytical formulas are derived to compute the model parameters for a canal pool possibly in backwater conditions. The IDZ model is compared with an accurate model on two example canals for different flow conditions. The comparisons done in the frequency domain and in the time domain show the accuracy of the model. 2003 Elsevier Inc. All rights reserved. Keywords: Open-channel flow; Hydraulic model; Analytical model; Laplace transform; Saint-Venant equations
1. Introduction Canal controller design requires a good model that captures the main dynamics of the system. Open-channel dynamics are usually described using non-linear partial derivatives equations (Saint-Venant equations). In order to use classical linear control design tools (such as frequency response, Nyquist plot, etc.), a linear model of the system is necessary. Such a model can be *
Corresponding author. Tel.: +33-467-04-63-47/00; fax: +33-467-63-57-95. E-mail addresses: [email protected] (X. Litrico), [email protected] (V. Fromion). 1 Tel.: +33-499-61-22-38; fax: +33-467-52-14-27.
0307-904X/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2003.10.014
678
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
obtained in different ways from Saint-Venant equations. The paper focuses on a simple analytical model for realistic cases (including backwater flow). Analytical models can be useful in a predesign phase, in order to provide simple design rules for the engineer. It can easily be simulated to give a (good) approximation of the dynamics of the system. Since the flow in a canal pool is subject to large variations due to the demand changes and perturbations, it is important when dealing with linear models to have a set of linear models, corresponding to various flow configurations. The easiest way to obtain this is to have model parameters expressed as a function of physical variables, such as the discharge and the water depth. This is the main objective of the paper: to obtain a good analytical approximation of SaintVenant model with parameters expressed as functions of physical variables. Many analytical models have been proposed in the literature. Most of them have been obtained in the frequency domain by applying Laplace transform to linearized Saint-Venant equations [1– 4]. The obtained transcendental function can then be simplified using various methods in order to get a model expressed as a rational function of s (the Laplace variable), possibly including a time delay. This approach followed by many authors (see e.g. [5–7]) can provide acceptable models in some cases, but is restricted to pools at uniform regime. Indeed, an analytical solution of SaintVenant equations only exists in the case of constant flow and water depth, which is a situation almost never encountered in practice. More realistic models can be obtained for non-uniform regimes by linearizing Saint-Venant equations and using a finite difference scheme [8,9]. This leads to a high order model, which cannot be expressed analytically as a function of physical parameters. It is therefore important to develop simple analytical models able to accurately reproduce the dynamic behavior of the system in realistic conditions, i.e. non-uniform flow. This goal was firstly identified by Schuurmans et al. [10], who derived analytical expressions of an integrator delay model for pools in non-uniform regimes. This model is based on simple approximations concerning the canal dynamics, and captures with two parameters the main frequency behavior of the system, i.e. a delayed integrator for the transfer between upstream discharge and downstream level. This was an important contribution to the modelling of open-channel systems, and the so-called integrator delay (ID) model has become a popular way to model a canal. However, if the approximations proposed by Schuurmans et al. [10] are efficient for time-domain simulations, they may not be suited for control design. From a control point of view, it is very important to precisely know the delay of a system, since it limits the achievable closed loop bandwidth (which is directly linked to the real-time performance) [11]. It has been shown that for a canal in uniform flow, the time delay s (in terms of automatic control) is equal to X =ðV0 þ C0 Þ, where X is the length of the pool, V0 the water velocity and C0 the celerity [12]. The time delay T obtained by [10] is equal to 2X =ð1 þ j0 ÞV0 , where j0 is a geometric coefficient linked to the form of the section (j0 ¼ 7=3 in a large rectangular canal). This time delay T can be much larger than s. Such a mismatch may lead to poorly tuned controllers, that may not perform correctly when applied on the real system. Even if the model proposed by [10] is indeed an approximate model, the question of the validity of the approximation remains open, as long as there is no way to evaluate the exact model (except for the uniform case). Recent works allow to have a very accurate frequency domain representation of linearized Saint-Venant equations for a canal pool in any flow regime [12]. This result shows that two approximations done by [10] are rather crude, and can be improved:
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
679
• For a pool in uniform regime, it is assumed that the system is a pure delay and that there is no integrator. However, there is always an integrator in the linearized Saint-Venant model. • For a pool completely affected by backwater, it is assumed that the delay is zero, while the system is a pure integrator. However, the delay which can be computed by the characteristics method is never equal to zero. These limitations also apply to a pool with intermediate flow conditions (i.e. partly affected by backwater and partly in uniform regime), since the above assumptions for each part of the system are also used: the part affected by uniform flow is supposed to be a pure delay, while the downstream part affected by backwater is supposed to be a pure integrator. It is possible to derive a new approximate model from the more general one of [12]. The paper presents an analytical integrator delay zero model for a canal pool that reproduces accurately the main frequency behavior of the system. The proposed model is multivariable, i.e. it can represent the upstream influence of a downstream modification, and accurately reproduces the behavior of the system in backwater flow configurations. The model parameters (integrator gains, delays and zeros) are explicitly given as functions of the steady-state discharge Q0 and the downstream water level YX and change continuously with those parameters. The approximate model is compared with the complete model for various flow conditions on two different canals: a long steep canal and a short and flat one, taken from the benchmarks proposed by [6].
2. Complete Saint-Venant model 2.1. Linearized Saint-Venant model For a prismatic channel, linearized Saint-Venant equations around a given steady flow regime (including backwater curves) are given by the following equations [13], where the values of parameters for the steady state regime are denoted with an underscored zero: T0 ðxÞ
oy oq þ ¼ 0; ot ox
oq oq oy þ 2V0 ðxÞ b0 ðxÞq þ ðC0 ðxÞ2 V0 ðxÞ2 ÞT0 ðxÞ c0 ðxÞy ¼ 0; ot ox ox
ð1Þ ð2Þ
where y ¼ Y Y0 is the variation of water depth (m), q ¼ Q Q0 variation of flow rate (m3 /s) from the steady state regime, defined by dY0 ðxÞ Sb Sf0 ðxÞ ; ¼ dx 1 F0 ðxÞ2
ð3Þ
T0 ðxÞ is the top width for the equilibrium regime (m), A0 ðxÞ is the wetted area (m2 ), Q0 the reference discharge (m3 /s) across section A0 , V0 ðxÞ the average velocity (m/s) in section A0 , Y0 ðxÞ the water depth (m), Sf0 ðxÞ the friction slope, Sb the bed slope, g the gravitational acceleration (m/s2 ),
680
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi C0 ¼ gA0 =T0 the wave celerity (m/s) and F0 ¼ V0 =C0 the Froude number. Throughout the paper the flow is assumed to be subcritical, i.e. F0 < 1. The friction slope Sf0 is modelled with Manning–Strickler formula [13]: Sf0 ðxÞ ¼
Q20 n2 2
4=3
A0 ðxÞ R0 ðxÞ
;
ð4Þ
with n the roughness coefficient (s/m1=3 ) and R0 ðxÞ the hydraulic radius (m), defined by R0 ¼ A0 =P0 , where P0 is the wetted perimeter (m). Taking into account Eq. (3), parameters c0 and b0 are expressed as oY0 2 oT0 2 þ gT0 ð1 þ j0 ÞSb ð1 þ j0 F0 ðj0 2ÞÞ c0 ¼ V0 ; ð5Þ ox ox and
2g oY0 Sb ; b0 ¼ V0 ox
ð6Þ
with j0 ¼ 7=3 4A0 =ð3T0 P0 ÞðdP0 =dY Þ. The boundary conditions are then given by qð0; tÞ ¼ q0 ðtÞ and qðX ; tÞ ¼ qX ðtÞ. The model of the system is therefore given by two linear partial differential equations. Before showing how to obtain a frequency domain representation of these equations, the considered input and outputs are recalled. 2.2. Inputs/outputs representation of a irrigation canal 2.2.1. Representation of a pool In the following, a canal pool is considered using an inputs/outputs approach: • The considered inputs are the upstream and downstream relative discharges u1 ðtÞ ¼ qð0; tÞ and u2 ðtÞ ¼ qðX ; tÞ, with X the length of the pool. • The outputs are the upstream and downstream relative water depths y1 ðtÞ ¼ yð0; tÞ and y2 ðtÞ ¼ yðX ; tÞ. The 2 · 2 transfer matrix relating the inputs to the outputs will be denoted by P ðX ; sÞ. 2.2.2. Hydraulic structures modelling The hydraulic structures are modelled as local linear relations between the relative discharge q and the relative water depth y at a given point. These relations are obtained by linearizing the hydraulic equations around a given functioning point. Let f denote the gate equation involving the water depth Y ðX ; tÞ, the discharge QðX ; tÞ and the gate opening W ðtÞ: QðX ; tÞ ¼ f ðY ðX ; tÞ; W ðtÞÞ;
ð7Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi for example, f ðY ðX ; tÞ; W ðtÞÞ ¼ Cd W ðtÞLv 2gðY ðX ; tÞ W ðtÞ=2Þ where Cd is a discharge coefficient (close to 0.6) and Lv is the gate width.
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
681
Linearizing this relation, the gate can be modelled as a local feedback between the discharge and the water depth: qðX ; tÞ ¼ kv yðX ; tÞ þ kw wðtÞ;
ð8Þ
where kv ¼ df ðy; wÞ=dy is the ‘‘feedback’’ gain of the gate and kw ¼ df ðy; wÞ=dw the gain of the gate opening. 2.2.3. Representation of a canal Finally, a canal is represented by a series of pools interconnected with hydraulic structures. In the sequel, it is shown how to obtain a frequency domain representation of these equations. 2.3. Transfer matrix representation of Saint-Venant model 2.3.1. Transfer matrix for uniform flow It is well known that a closed-form solution of linearized Saint-Venant equations can be obtained in the frequency domain for uniform flow conditions. Applying Laplace transform to the linear partial differential equations (1) and (2), and reordering them lead to an ordinary differential equation in the variable x, with a complex parameter s (the Laplace variable): d qðx; sÞ qðx; sÞ ; ð9Þ ¼ As ðxÞ yðx; sÞ dx yðx; sÞ with
0 As ðxÞ ¼ @
0 s þ b0 ðxÞ T0 ðxÞðC0 ðxÞ2 V0 ðxÞ2 Þ
1 T0 ðxÞs 2V0 ðxÞT0 ðxÞs þ c0 ðxÞ A: T0 ðxÞðC0 ðxÞ2 V0 ðxÞ2 Þ
Boundary conditions are given by qð0; sÞ ¼ q0 ðsÞ and qðX ; sÞ ¼ qX ðsÞ. In the general case, matrix As depends on x, thus this ordinary differential equation cannot be solved analytically, except for one particular case (when As does not depend on x, i.e. for uniform flow conditions). When As is constant (i.e. for uniform flow conditions), the transfer matrix is given by (see Appendix A for details): qð0; sÞ yð0; sÞ p11 ðX ; sÞ p12 ðX ; sÞ ; ð10Þ ¼ qðX ; sÞ p21 ðX ; sÞ p22 ðX ; sÞ yðX ; sÞ with p11 ðX ; sÞ ¼
k2 ek1 X k1 ek2 X ; T0 sðek2 X ek1 X Þ
ð11Þ
p12 ðX ; sÞ ¼
k1 k2 ; T0 sðek2 X ek1 X Þ
ð12Þ
p21 ðX ; sÞ ¼
ðk2 k1 Þeðk1 þk2 ÞX ; T0 sðek2 X ek1 X Þ
ð13Þ
682
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
p22 ðX ; sÞ ¼
k1 ek1 X k2 ek2 X : T0 sðek2 X ek1 X Þ
ð14Þ
This result has a very limited domain of application, since real canals are seldom in uniform regime. The authors have proposed a way to compute the solution in more realistic cases, where backwater curves occur, i.e. when matrix As ðxÞ varies with x [12]. This method is briefly explained before exposing the approximate model. 2.3.2. Numerical solution in the general case The numerical solution uses a space discretization, and the solution is obtained in each part by evaluating the transfer matrix as if the water depth was constant: this enables to use the same formulas as in the uniform case to compute the solution. The parts are then interconnected in order to get the whole transfer matrix. Let xk be a space discretization of interval ½0; X into n subintervals: 0 ¼ x0 < x1 < < xk < < xn ¼ X : This numerical solution uses the analytical solution derived above on each subinterval [xk ; xkþ1 ], considering the matrix As constant on each subinterval, which leads to the overall transition matrix: 1 Y Cs ðX ; 0Þ ¼ Cs ðxk ; xk1 Þ; ð15Þ k¼n
where matrices Cs ðxk ; xk1 Þ are obtained with the analytical method (Eq. (A.1)), using the value of coefficients in x ¼ xk . This procedure enables to compute the transition matrix Cs ðX ; 0Þ for each desired frequency s ¼ jx with a good numerical accuracy and limited computation time [12]. This complete model however relies on a large number of computations, since the backwater curve is needed at each discretization point to obtain the transfer matrix. To obtain a tractable analytical model, simple approximations of the backwater profile and of the transfer matrix are considered in the following. To get a good numerical solution, it is important to choose the number of space steps according to the variation of the water depth. In the uniform regime, only one evaluation of matrix As is necessary. In non-uniform regimes, it is necessary to evaluate this matrix at several locations. The approximations will be done in the following by considering only two locations. 3. Analytical approximate model Two different types of approximation are considered in the following: • Approximation of the backwater profile by a piecewise linear function. • Low and high frequencies approximations of the transcendental transfer function to obtain rational models with time delay. The space approximation will firstly be exposed, then the frequency approximations will be presented, leading to the global approximate model.
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
683
3.1. Backwater profile 3.1.1. Backwater profile approximation The backwater profile can be computed accurately by solving the ordinary differential equation (3), which necessitates a space discretization in many parts. In order to get simple analytical expressions, the pool is here separated in two parts: • The upstream part, in uniform flow: the water surface is approximated by a line parallel to the bottom. • The downstream part, with backwater flow: the water surface is approximated by a line tangent to the real backwater curve at x ¼ X (downstream end of the pool). Schuurmans et al. [10] used similar assumptions, but using a horizontal line for the downstream part. Using the backwater curve definition given by Eq. (3), let SX be the deviation from bed slope of the line tangent to the backwater curve at the downstream end of the pool: SX ¼
Sb Sf0 ðX Þ : 1 F02 ðX Þ
ð16Þ
Then, the intersection between the two lines occurs at abscissa x1 : 8
YX Yn < ;0 if SX 6¼ 0; max X x1 ¼ SX : X if SX ¼ 0:
ð17Þ
Then, the approximate backwater curve satisfies the equation:
Y1 for x 2 ½0; x1 ; e Y ðxÞ ¼ Y1 þ ðx x1 ÞSb for x 2 ½x1 ; X with
Y1 ¼
Yn YX XSX
if x1 6¼ 0; if x1 ¼ 0:
ð18Þ
Here, Y1 ¼ Yn when part of the pool is in uniform flow, and when the whole pool is affected by backwater, Y1 ¼ YX XSX i.e. the water depth is a straight line of slope SX . The corresponding approximation of the backwater profile is schematized in Fig. 1. 3.1.2. Computation of the uniform depth Yn To obtain an analytical model, it is necessary to have an analytical expression for the uniform depth Yn . It is not possible to have an analytical solution in the general case, but Yn can either be approximated by an analytic expression (e.g. with the hydraulic exponent method [13]), or be given by tables for a considered canal pool. The uniform depth is computed from the solution of the algebraic equation Sf0 ¼ Sb , with Sf0 the friction slope given by Manning–Strickler formula (4).
684
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
Fig. 1. Backwater profile (––) and its approximation (- - -).
In the general case, Sf0 is not a simple function of Y0 , therefore there is no analytical solution of this equation. The hydraulic exponent method [13] enables to obtain a good approximation of the pffiffiffiffiffiffi uniform depth in prismatic channels. In this method, the conveyance K0 ¼ Q0 = Sf0 is assumed to be a power function of the depth of flow Y0 : K02 ¼ CY0N ; where C is a coefficient and N is the hydraulic exponent for uniform flow computation. A good approximation of N is obtained by computing the uniform depth for two different discharges Q1 and Q2 . N is then given by N ¼2
logðQ1 =Q2 Þ : logðYn ðQ1 Þ=Yn ðQ2 ÞÞ
ð19Þ
Given the uniform flow depth Yn ðQ1 Þ corresponding to a discharge Q1 (e.g. the maximum discharge for the considered pool), the values of the uniform flow depth Yn ðQ0 Þ for other discharges Q0 can be approximated by 2=N Q0 Yn ðQ0 Þ Yn ðQ1 Þ: ð20Þ Q1 3.1.3. Interconnection of two parts To clarify the interconnection of transition matrices (15), let us consider the case of a pool separated into two parts: • The upstream part of the pool is modelled by the transfer matrix given by Eqs. (11)–(14) denoted P ðx1 ; sÞ. • The downstream part is modelled by a transfer matrix P ðX x1 ; sÞ obtained by evaluating matrix As ðxÞ at x2 2 ½x1 ; X . Let expand As ðxÞ in Taylor series around x2 2 ½x1 ; X : As ðxÞ ¼ As ðx2 Þ þ
dAs ðx2 Þ ðx x2 Þ þ Oðx x2 Þ2 : dx
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
685
Choosing x2 ¼ x1 þX enables to minimize the effect of the first order term dAdxs ðx2 Þ ðx x2 Þ on the final 2 result. The real transfer is therefore approximated by the transfer obtained using the formulas in . uniform flow, but with matrix As ðx2 Þ, where x2 ¼ x1 þX 2 In the following, pij ðsÞ denotes the terms of transfer matrix corresponding to the upstream part, pij ðsÞ those corresponding to the downstream part and ^pij ðsÞ those corresponding to the whole pool. The equivalent transfer matrix Pb ðX ; sÞ is therefore obtained by multiplication of the corresponding transition matrices (as in (15)) and rearranging in order to get the desired inputs and outputs, which leads to ^ p11 ¼ p11 þ ^ p12 ¼ ^ p21 ¼
p12 p21 ; p11 p22
p12 p12 ; p11 p22
p21 p21 ; p11 p22
^ p22 p22 ¼
p12 p21 : p11 p22
ð21Þ ð22Þ ð23Þ ð24Þ
These formulas will be used to get the overall approximation Pb ðX ; sÞ from the approximations of P ðx1 ; sÞ and P ðX x1 ; sÞ. 3.2. Frequency domain approximations 3.2.1. About approximate models at low frequencies A classical approach followed by many authors [2,3,9,10,14] is to use the moment matching method in zero (or analog method) to obtain rational approximations of the irrational transfer function. This leads to a model that fits the low frequency response of the system. In the case of the integrator delay model of Schuurmans et al., since the phase is only linked to the delay the approximation will take the group delay of the system as the delay value (the group delay of GðjxÞ is defined by o\GðjxÞ=ox [15]). However, the group delay also incorporates phase changes due to zeros and poles. This explains why the delay T computed by [10] is greater than the real delay s. Even if this approximation is mathematically valid, a delay is qualitatively completely different from stable poles or zeros from a control point of view. In fact, a delay introduces a structural limitation on the real-time performance of a controlled dynamical system [11], which is not the case of stable poles or zeros. This is why it is important to model the delay as exactly as possible, and add zeros or poles to better fit the phase shift. As shown in [12], the exact delay can be obtained by the characteristics. In the following, low and high frequencies approximations of the transfer matrix P are derived for constant flow conditions. These approximations, together with the interconnection rules (21)– (24) enable to get the overall approximate model.
686
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
3.2.2. Low frequencies approximation For low frequencies, the behavior of the transfer matrix is dominated by the integrator and the delays. The relation between the downstream (resp. upstream) discharge and the downstream level can be considered as an integrator (resp. a delayed integrator) whose gain is the variation of volume with respect to the variation of downstream level. The situation is symmetrical for the upstream level. The integrator gains will be called ‘‘equivalent backwater areas’’ since they have the dimension of an area. They can be close to the pool water surface in specific cases (e.g. at uniform flow for small pools), but are usually different from the backwater area. Variables related to the downstream water level are denoted with a subscript d, whereas variables related to the upstream water level are denoted with a subscript u. For a pool with constant water depth, the low frequency approximation reads: p11 ðsÞ ¼
1 ; Au s
p12 ðsÞ ¼ p21 ðsÞ ¼
esu s ; Au s
ð25Þ ð26Þ
esd s ; Ad s
ð27Þ
1 Ad s
ð28Þ
p22 ðsÞ ¼
with the equivalent backwater areas, obtained by computing ðspij ðsÞÞð0Þ Ad ¼
c0 T02 ðC02 V02 Þ 2 2 x 1 e T0 ðC0 V0 Þ ; c0
c T02 ðC02 V02 Þ T0 ðC20V 2 Þx Au ¼ e 0 0 1 ; c0
ð29Þ
ð30Þ
and the delays sd ¼
x ; C0 þ V 0
ð31Þ
su ¼
x : C0 V 0
ð32Þ
Using this approximation for each part of the pool together with Eqs. (21)–(24) for the interconnection, the approximation for the transfer matrix Pb ðsÞ at low frequencies is ^ p11 ðsÞ ¼
1 ; b Aus
ð33Þ
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
^ p12 ðsÞ ¼
^ p21 ðsÞ ¼
e^su s ; b us A
e^sd s ; b ds A
^ p22 ðsÞ ¼
1 : b Ads
The equivalent areas are given by b d ¼ Ad 1 þ Ad ; A Au A u b u ¼ Au 1 þ A ; Ad
687
ð34Þ
ð35Þ
ð36Þ
ð37Þ
ð38Þ
Au , Ad are given by Eqs. (29) and (30) where the variables are evaluated at 0 and where x ¼ x1 . Au , Ad are given by the same equations (29) and (30) where the variables are evaluated at x2 and where x ¼ X x1 . The equivalent delays are given by ^sd ¼ sd þ sd ;
ð39Þ
^su ¼ su þ su :
ð40Þ
The delays for the upstream uniform part sd and su are given by Eqs. (31) and (32) where the variables are evaluated at 0 and where x ¼ x1 . The delays for the downstream backwater part sd , su are given by the same equations (31) and (32) where the variables are evaluated at x2 and where x ¼ X x1 . 3.2.3. High frequencies approximation For high frequencies, the delay and the gravity waves are predominant in the transfer matrix elements. It is not easy to obtain a simple approximation of the oscillatory behavior of the gravity waves. For simplicity, only the static behavior of the oscillation is represented. This is modelled by constant gains in high frequencies, denoted pij1 . These gains are obtained by computing the average value of the modulus of each transfer for high frequencies. Indeed, limjsj!1 jpij ðsÞj may not exist, since jpij ðsÞj may oscillate at infinity. The oscillating terms are sinusoidal, therefore the average value is computed by nullifying them. For a pool with constant water depth, one then gets the approximation: p11 ðsÞ ¼ p111 ;
ð41Þ
p12 ðsÞ ¼ p121 esu s ;
ð42Þ
p21 ðsÞ ¼ p211 esd s ;
ð43Þ
p22 ðsÞ ¼ p221 :
ð44Þ
688
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
Using this approximation for each part of the pool together with Eqs. (21)–(24) for the interconnection, the transfer matrix Pb ðsÞ at high frequencies is approximated by ^ p111 ; p11 ðsÞ ^
ð45Þ
^ p121 e^su s ; p12 ðsÞ ^
ð46Þ
^ p211 e^sd s ; p21 ðsÞ ^
ð47Þ
^ p221 : p22 ðsÞ ^
ð48Þ
The gains in high frequencies are given by: ^ p111 ¼ p111 þ
p121 p211 ; p111 þ p221
ð49Þ
^ p121 ¼
p121 p121 ; p111 þ p221
ð50Þ
^ p211 ¼
p211 p211 ; p111 þ p221
ð51Þ
^ p221 ¼ p221 þ
p121 p211 p111 þ p221
ð52Þ
with
p111 ¼
1 T0 C0 ð1 F0 Þ
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u1 þ 1F0 2 eax t 1þF0 1 þ eax
p121
c0 x 2T0 ðC 2 V 2 Þ 0 0
T0 ð2 þ ðj0 1ÞF02 ÞSb A0 F0 ð1 F02 Þ
and x ¼ x1 .
ð53Þ
x
2 e pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2 T0 C0 ð1 F0 Þ 1 þ eax vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u1 þ 1þF0 2 eax t 1F0 1 : ¼ T0 C0 ð1 þ F0 Þ 1 þ eax
In the upstream uniform part, a¼
2
2 e 2T0 ðC0 V0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ¼ T0 C0 ð1 F02 Þ 1 þ eax
p211 ¼
p221
c0 2
;
ð54Þ
ð55Þ
ð56Þ
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
689
For the downstream backwater part, the same expressions are valid, replacing x with x ¼ X x1 , a with T0 A0 dT0 2 2 a¼ ð2 þ ðj0 1ÞF0 ÞSb 2 þ ðj0 1ÞF0 þ j0 2 F04 SX 2 2 A0 F0 ð1 F0 Þ T0 dY and all the variables by their value taken at x2 . 3.3. Global approximate model The global model is simply obtained by adding the approximations in high and low frequencies: 1 þ^ p111 ; b us A
^ p11 ðsÞ ¼
ð57Þ
1 ^ þ^ p121 e^su s ; p12 ðsÞ ¼ b us A
ð58Þ
^ p21 ðsÞ ¼
1 þ^ p211 e^sd s ; b ds A
ð59Þ
^ p22 ðsÞ ¼
1 ^ p221 : b ds A
ð60Þ
The approximate model is now completely presented. Given a prismatic canal with known geometry and Manning coefficient, a discharge Q0 and a downstream limit condition YX , Eqs. (16)–(20) enable to compute the integrator delay zero model parameters: the delays ^sd and ^su given b u , given by Eqs. (37) and (38) and b d and A by Eqs. (39) and (40), the equivalent backwater areas A the gains at high frequencies ^ pij1 , given by Eqs. (49)–(52). The accuracy of this approximate model will now be evaluated.
4. Validation of the approximate model 4.1. Example canals The paper is illustrated with two trapezoidal prismatic channels, with different characteristics (see Table 1, where X is the channel length (m), m the bank slope, B the bed width (m), Sb the bed
Table 1 Parameters for the two example canals Canal 1 Canal 2
X
m
B
Sb
n
Yn
Qmax
Qmin
3000 6000
1.5 1.5
7 8
0.0001 0.0008
0.02 0.02
2.12 2.92
14 80
1.75 10
690
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
slope, n the roughness coefficient (s m1=3 ), Yn the uniform depth (m) corresponding to the maximum discharge Qmax (m3 /s)). Canal 1 is a short oscillating canal, and canal 2 is a long delayed canal. The downstream limit condition YX is chosen equal to the uniform depth Yn for Qmax . The uniform depth approximations are obtained using Eqs. (19) and (20) where Q1 ¼ Qmax and Q2 ¼ Qmax =8. This leads to N ¼ 3:51 for canal 1 and N ¼ 3:56 for canal 2. The obtained approximations of the uniform depths are very accurate for both canals (see Fig. 2). The approximate backwater curves obtained in this case for an average discharge Qmax =2 are depicted in Fig. 3 for both canals.
Fig. 2. Variation of the uniform depth with the discharge, canals 1 and 2. Comparison between the exact (––) and approximate (- - -) values.
Fig. 3. Backwater curve (––) and its approximation (- - -) for Qmax =2, canals 1 and 2.
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
691
4.2. Model evaluation 4.2.1. Comparison of frequency responses The Bode plots of canals 1 and 2 are obtained in non-uniform flow conditions for Q0 ¼ Qmax =2 using the complete numerical method [12] (see Fig. 4). One can notice the behavior of both canals, which is completely different: canal 1 is oscillating, with resonant modes, and canal 2 exhibits a long time delay. The approximate model reproduces very well the global frequency behavior of the system, especially in the case of canal 1. Canal 2 has a pole that is not taken into account by the approximate model. In order to show the accuracy of the method, it is important to see how the parameters vary with the hydraulic conditions, i.e. the discharge Q0 . 4.2.2. Parameters comparison Let us consider a constant downstream limit condition YX corresponding to the uniform flow for Qmax , and vary the discharge Q0 . The canal is therefore in uniform flow for Q0 ¼ Qmax and the backwater curve increases as the discharge diminishes. The approximate parameters are compared with their exact values when they are available. b d , ^sd and lim supx!1 jp21 ðjxÞj Figs. 5–7 respectively depict the exact and approximate values of A and lim inf x!1 jp21 ðjxÞj for both canals. It is clear from the figures that the approximation is very accurate for the delay ^sd and for b d . The gain at high frequencies ^ p211 is rather well estimated for canal 2, and is in coefficient A between min and max values for canal 1. In fact, the max value of the gain at high frequencies ^p211 increases a lot as the discharge tends towards zero, since the damping of resonant modes diminishes. This is not taken into account by the approximate model, which is sufficient for control purposes, since the resonant modes will not be controlled. These resonant modes are usually filtered in order not to destabilize the controlled system.
Fig. 4. Bode plot of transfer p21 ðsÞ (––) and its approximation e^sd s
1 bA d s
þ^ p211 (- - -), canals 1 and 2 for Q0 ¼ Qmax =2.
692
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
b d with the discharge, for a downstream limit condition YX corresponding to the uniFig. 5. Variation of coefficient A form flow for Qmax , canals 1 and 2. Comparison between the exact (––) and approximate (- - -) values.
Fig. 6. Variation of delay ^sd with the discharge, for a downstream limit condition YX corresponding to the uniform flow for Qmax , canals 1 and 2. Comparison between the exact (––) and approximate (- - -) values.
4.2.3. Comparison of time-domain simulations The approximate model is compared to a finite dimensional model obtained by rational approximation of the exact model [12]. Simulations are done on Matlab for canals 1 and 2, for Q0 ¼ Qmax =2 (see Fig. 8). The models are connected with a downstream gate (equivalent to a local feedback) with a gain kv ¼ 6:84. The simulation is very accurate for both canals: in fact, as can be seen from the Bode plot, the whole frequency response is very well fitted by the approximate model, except the resonant modes
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
693
Fig. 7. Variation of gain ^p211 with the discharge, for a downstream limit condition YX corresponding to the uniform flow for Qmax , canals 1 and 2. Comparison between min and max values (––) and the approximation (- - -).
Fig. 8. Step responses of the rational approximation of the complete model (––) and IDZ model (- - -), canals 1 and 2.
for canal 1 and a supplementary pole for canal 2. These discrepancies occur in rather high frequencies and at a small gain, this is why they are not clearly visible in the step response.
5. Conclusion The paper proposes a new analytical approximate model for an open-channel pool subject to backwater. This integrator delay zero model is simple, and very accurate. It is obtained from rigorous mathematical development on linearized Saint-Venant equations. It models the complete transfer matrix (i.e. the downstream influence is also taken into account). The inputs are the
694
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
upstream and downstream discharges, and its outputs are the upstream and downstream water levels. The model has been validated in frequency and time domain using the complete transfer matrix obtained following [12] on two different example canals. The approximate model can be used for controller design, such as distant downstream PI controllers, or even for advanced controller design (e.g. multivariable controllers).
Acknowledgements This work was partially supported by the joint research program INRA/Cemagref ASS AQUAE no 02 on the control of delayed hydraulic systems.
Appendix A. Saint-Venant transfer matrix: analytic solution for uniform flow In uniform conditions, since matrix As does no longer depend on x, the ordinary differential equation (9) can be solved analytically. To do this, let us diagonalize matrix As As ¼ Ps Ds Ps1 ; with
Ds ¼
k1 ðsÞ 0
0 ; k2 ðsÞ
Ps ¼
T0 s k1 ðsÞ
T0 s ; k2 ðsÞ
Ps1
1 ¼ T0 sðk1 ðsÞ k2 ðsÞÞ
k2 ðsÞ k1 ðsÞ
T0 s : T0 s
The eigenvalues are given by k1;2 ðsÞ ¼
2T0 V0 s þ c0
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4C02 T02 s2 þ 4T0 ðV0 c0 ðC02 V02 ÞT0 b0 Þs þ c20 ; 2T0 ðC02 V02 Þ
Once the matrix As is diagonalized, the ordinary differential equation (9) can be solved analytically: qðx; sÞ qð0; sÞ ¼ Cs ðx; 0Þ ðA:1Þ yðx; sÞ yð0; sÞ with 0
Cs ðx; 0Þ :¼ Ps eDs x Ps1
k1 ek2 x k2 ek1 x B k1 k2 B ¼B @ k1 k2 ðek1 x ek2 x Þ T0 sðk1 k2 Þ
1 T0 sðek2 x ek1 x Þ C k1 k2 C C: k1 x k2 x A k1 e k2 e k1 k2
This is the classical result already obtained by many authors [2,3,5,6,10]. Reordering the result in order to have the initial boundary conditions as inputs (qð0; sÞ and qðX ; sÞ) leads to the Saint-Venant transfer matrix:
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
yð0; sÞ yðX ; sÞ
!
0
k2 ek1 X k1 ek2 X B T sðek2 X ek1 X Þ B 0 ¼B @ ðk2 k1 Þeðk1 þk2 ÞX
1 k1 k2 T0 sðek2 X ek1 X Þ C C C k1 X k2 X A k1 e k2 e
qð0; sÞ qðX ; sÞ
695
! :
ðA:2Þ
T0 sðek2 X ek1 X Þ T0 sðek2 X ek1 X Þ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} P ðX ; sÞ ¼ ðpij ðX ; sÞÞ
References [1] M. Shand, Automatic Downstream Control Systems for Irrigation Canals, Ph.D. Thesis, University of California, Berkeley, 1971. [2] G. Corriga, D. Salembeni, S. Sanna, G. Usai, A control method for speeding up response of hydroelectric stations power canals, Appl. Math. Modell. 12 (1988) 627–633. [3] Y. Ermolin, Study of open-channel dynamics as controlled process, J. Hydraul. Eng. 118 (1) (1992) 59–71. [4] S. Hancu, P. Dan, Wave-motion stability in canals with automatic controllers, J. Hydraul. Eng. 118 (12) (1992) 1621–1638. [5] G. Corriga, A. Fanni, S. Sanna, G. Usai, A constant-volume control method for open channel operation, Int. J. Modell. Simulat. 2 (2) (1982) 108–112. [6] J. Baume, J. Sau, P.-O. Malaterre, Modeling of irrigation channel dynamics for controller design, in: Conf. on Systems, Man and Cybernetics, SMCÕ98, San Diego, CA, 1998, pp. 3856–3861. [7] C. Seatzu, Design and robustness analysis of decentralized constant volume-control for open-channels, Appl. Math. Modell. 23 (6) (1999) 479–500. [8] O. Balogun, Design of Real-Time Feedback Control for Canal Systems Using Linear Quadratic Regulator Theory, Ph.D. Thesis, Department of Mechanical Engineering, University of California at Davis, 1985, 230 p. [9] P.-O. Malaterre, Modelisation, Analyse et Commande Optimale LQR dÕun Canal dÕIrrigation, Ph.D. Thesis, ISBN 2-85362-368-8, Etude EEE no. 14, LAAS, CNRS, ENGREF, Cemagref (in French), 1994. [10] J. Schuurmans, O. Bosgra, R. Brouwer, Open-channel flow model approximation for controller design, Appl. Math. Modell. 19 (1995) 525–530. [11] K. Astr€ om, Limitations on control system performance, Eur. J. Control 6 (2000) 1–19. [12] X. Litrico, V. Fromion, Infinite dimensional modelling of open-channel hydraulic systems for control purposes, in: 41st Conf. on Decision and Control, Las Vegas, 2002, pp. 1681–1686. [13] V. Chow, Open-channel Hydraulics, McGraw-Hill Book Company, New York, 1988, 680 p. [14] X. Litrico, D. Georges, Robust continuous-time and discrete-time flow control of a dam–river system (I): Modelling, Appl. Math. Modell. 23 (11) (1999) 809–827. [15] J. Freudenberg, D. Looze, A sensitivity tradeoff for plants with time delay, IEEE Trans. Automat. Control 32 (2) (1987) 99–104.
Analytical approximation of open-channel flow for controller design Xavier Litrico a
a,*
, Vincent Fromion
b,1
Cemagref, UR Irrigation, 361 rue-J-F Breton, B.P. 5095, 34033 Montpellier Cedex 1, France b INRA, LASB, 2 place Viala, 34060 Montpellier, France
Received 13 March 2003; received in revised form 2 October 2003; accepted 27 October 2003 Available online 21 January 2004
Abstract The integrator delay (ID) model is a popular and simple way to model a canal for control purposes. In this paper, particular attention is given to accurately model the delay and the integrator gain in any flow configuration. Based on a previous work by the authors allowing to have a very accurate frequency domain representation of Saint-Venant equations, the paper proposes a new approximate model: the integrator delay zero (IDZ) model has an integrator and a delay in low frequencies, and models the high frequencies by a constant gain and a delay. Analytical formulas are derived to compute the model parameters for a canal pool possibly in backwater conditions. The IDZ model is compared with an accurate model on two example canals for different flow conditions. The comparisons done in the frequency domain and in the time domain show the accuracy of the model. 2003 Elsevier Inc. All rights reserved. Keywords: Open-channel flow; Hydraulic model; Analytical model; Laplace transform; Saint-Venant equations
1. Introduction Canal controller design requires a good model that captures the main dynamics of the system. Open-channel dynamics are usually described using non-linear partial derivatives equations (Saint-Venant equations). In order to use classical linear control design tools (such as frequency response, Nyquist plot, etc.), a linear model of the system is necessary. Such a model can be *
Corresponding author. Tel.: +33-467-04-63-47/00; fax: +33-467-63-57-95. E-mail addresses: [email protected] (X. Litrico), [email protected] (V. Fromion). 1 Tel.: +33-499-61-22-38; fax: +33-467-52-14-27.
0307-904X/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2003.10.014
678
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
obtained in different ways from Saint-Venant equations. The paper focuses on a simple analytical model for realistic cases (including backwater flow). Analytical models can be useful in a predesign phase, in order to provide simple design rules for the engineer. It can easily be simulated to give a (good) approximation of the dynamics of the system. Since the flow in a canal pool is subject to large variations due to the demand changes and perturbations, it is important when dealing with linear models to have a set of linear models, corresponding to various flow configurations. The easiest way to obtain this is to have model parameters expressed as a function of physical variables, such as the discharge and the water depth. This is the main objective of the paper: to obtain a good analytical approximation of SaintVenant model with parameters expressed as functions of physical variables. Many analytical models have been proposed in the literature. Most of them have been obtained in the frequency domain by applying Laplace transform to linearized Saint-Venant equations [1– 4]. The obtained transcendental function can then be simplified using various methods in order to get a model expressed as a rational function of s (the Laplace variable), possibly including a time delay. This approach followed by many authors (see e.g. [5–7]) can provide acceptable models in some cases, but is restricted to pools at uniform regime. Indeed, an analytical solution of SaintVenant equations only exists in the case of constant flow and water depth, which is a situation almost never encountered in practice. More realistic models can be obtained for non-uniform regimes by linearizing Saint-Venant equations and using a finite difference scheme [8,9]. This leads to a high order model, which cannot be expressed analytically as a function of physical parameters. It is therefore important to develop simple analytical models able to accurately reproduce the dynamic behavior of the system in realistic conditions, i.e. non-uniform flow. This goal was firstly identified by Schuurmans et al. [10], who derived analytical expressions of an integrator delay model for pools in non-uniform regimes. This model is based on simple approximations concerning the canal dynamics, and captures with two parameters the main frequency behavior of the system, i.e. a delayed integrator for the transfer between upstream discharge and downstream level. This was an important contribution to the modelling of open-channel systems, and the so-called integrator delay (ID) model has become a popular way to model a canal. However, if the approximations proposed by Schuurmans et al. [10] are efficient for time-domain simulations, they may not be suited for control design. From a control point of view, it is very important to precisely know the delay of a system, since it limits the achievable closed loop bandwidth (which is directly linked to the real-time performance) [11]. It has been shown that for a canal in uniform flow, the time delay s (in terms of automatic control) is equal to X =ðV0 þ C0 Þ, where X is the length of the pool, V0 the water velocity and C0 the celerity [12]. The time delay T obtained by [10] is equal to 2X =ð1 þ j0 ÞV0 , where j0 is a geometric coefficient linked to the form of the section (j0 ¼ 7=3 in a large rectangular canal). This time delay T can be much larger than s. Such a mismatch may lead to poorly tuned controllers, that may not perform correctly when applied on the real system. Even if the model proposed by [10] is indeed an approximate model, the question of the validity of the approximation remains open, as long as there is no way to evaluate the exact model (except for the uniform case). Recent works allow to have a very accurate frequency domain representation of linearized Saint-Venant equations for a canal pool in any flow regime [12]. This result shows that two approximations done by [10] are rather crude, and can be improved:
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
679
• For a pool in uniform regime, it is assumed that the system is a pure delay and that there is no integrator. However, there is always an integrator in the linearized Saint-Venant model. • For a pool completely affected by backwater, it is assumed that the delay is zero, while the system is a pure integrator. However, the delay which can be computed by the characteristics method is never equal to zero. These limitations also apply to a pool with intermediate flow conditions (i.e. partly affected by backwater and partly in uniform regime), since the above assumptions for each part of the system are also used: the part affected by uniform flow is supposed to be a pure delay, while the downstream part affected by backwater is supposed to be a pure integrator. It is possible to derive a new approximate model from the more general one of [12]. The paper presents an analytical integrator delay zero model for a canal pool that reproduces accurately the main frequency behavior of the system. The proposed model is multivariable, i.e. it can represent the upstream influence of a downstream modification, and accurately reproduces the behavior of the system in backwater flow configurations. The model parameters (integrator gains, delays and zeros) are explicitly given as functions of the steady-state discharge Q0 and the downstream water level YX and change continuously with those parameters. The approximate model is compared with the complete model for various flow conditions on two different canals: a long steep canal and a short and flat one, taken from the benchmarks proposed by [6].
2. Complete Saint-Venant model 2.1. Linearized Saint-Venant model For a prismatic channel, linearized Saint-Venant equations around a given steady flow regime (including backwater curves) are given by the following equations [13], where the values of parameters for the steady state regime are denoted with an underscored zero: T0 ðxÞ
oy oq þ ¼ 0; ot ox
oq oq oy þ 2V0 ðxÞ b0 ðxÞq þ ðC0 ðxÞ2 V0 ðxÞ2 ÞT0 ðxÞ c0 ðxÞy ¼ 0; ot ox ox
ð1Þ ð2Þ
where y ¼ Y Y0 is the variation of water depth (m), q ¼ Q Q0 variation of flow rate (m3 /s) from the steady state regime, defined by dY0 ðxÞ Sb Sf0 ðxÞ ; ¼ dx 1 F0 ðxÞ2
ð3Þ
T0 ðxÞ is the top width for the equilibrium regime (m), A0 ðxÞ is the wetted area (m2 ), Q0 the reference discharge (m3 /s) across section A0 , V0 ðxÞ the average velocity (m/s) in section A0 , Y0 ðxÞ the water depth (m), Sf0 ðxÞ the friction slope, Sb the bed slope, g the gravitational acceleration (m/s2 ),
680
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi C0 ¼ gA0 =T0 the wave celerity (m/s) and F0 ¼ V0 =C0 the Froude number. Throughout the paper the flow is assumed to be subcritical, i.e. F0 < 1. The friction slope Sf0 is modelled with Manning–Strickler formula [13]: Sf0 ðxÞ ¼
Q20 n2 2
4=3
A0 ðxÞ R0 ðxÞ
;
ð4Þ
with n the roughness coefficient (s/m1=3 ) and R0 ðxÞ the hydraulic radius (m), defined by R0 ¼ A0 =P0 , where P0 is the wetted perimeter (m). Taking into account Eq. (3), parameters c0 and b0 are expressed as oY0 2 oT0 2 þ gT0 ð1 þ j0 ÞSb ð1 þ j0 F0 ðj0 2ÞÞ c0 ¼ V0 ; ð5Þ ox ox and
2g oY0 Sb ; b0 ¼ V0 ox
ð6Þ
with j0 ¼ 7=3 4A0 =ð3T0 P0 ÞðdP0 =dY Þ. The boundary conditions are then given by qð0; tÞ ¼ q0 ðtÞ and qðX ; tÞ ¼ qX ðtÞ. The model of the system is therefore given by two linear partial differential equations. Before showing how to obtain a frequency domain representation of these equations, the considered input and outputs are recalled. 2.2. Inputs/outputs representation of a irrigation canal 2.2.1. Representation of a pool In the following, a canal pool is considered using an inputs/outputs approach: • The considered inputs are the upstream and downstream relative discharges u1 ðtÞ ¼ qð0; tÞ and u2 ðtÞ ¼ qðX ; tÞ, with X the length of the pool. • The outputs are the upstream and downstream relative water depths y1 ðtÞ ¼ yð0; tÞ and y2 ðtÞ ¼ yðX ; tÞ. The 2 · 2 transfer matrix relating the inputs to the outputs will be denoted by P ðX ; sÞ. 2.2.2. Hydraulic structures modelling The hydraulic structures are modelled as local linear relations between the relative discharge q and the relative water depth y at a given point. These relations are obtained by linearizing the hydraulic equations around a given functioning point. Let f denote the gate equation involving the water depth Y ðX ; tÞ, the discharge QðX ; tÞ and the gate opening W ðtÞ: QðX ; tÞ ¼ f ðY ðX ; tÞ; W ðtÞÞ;
ð7Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi for example, f ðY ðX ; tÞ; W ðtÞÞ ¼ Cd W ðtÞLv 2gðY ðX ; tÞ W ðtÞ=2Þ where Cd is a discharge coefficient (close to 0.6) and Lv is the gate width.
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
681
Linearizing this relation, the gate can be modelled as a local feedback between the discharge and the water depth: qðX ; tÞ ¼ kv yðX ; tÞ þ kw wðtÞ;
ð8Þ
where kv ¼ df ðy; wÞ=dy is the ‘‘feedback’’ gain of the gate and kw ¼ df ðy; wÞ=dw the gain of the gate opening. 2.2.3. Representation of a canal Finally, a canal is represented by a series of pools interconnected with hydraulic structures. In the sequel, it is shown how to obtain a frequency domain representation of these equations. 2.3. Transfer matrix representation of Saint-Venant model 2.3.1. Transfer matrix for uniform flow It is well known that a closed-form solution of linearized Saint-Venant equations can be obtained in the frequency domain for uniform flow conditions. Applying Laplace transform to the linear partial differential equations (1) and (2), and reordering them lead to an ordinary differential equation in the variable x, with a complex parameter s (the Laplace variable): d qðx; sÞ qðx; sÞ ; ð9Þ ¼ As ðxÞ yðx; sÞ dx yðx; sÞ with
0 As ðxÞ ¼ @
0 s þ b0 ðxÞ T0 ðxÞðC0 ðxÞ2 V0 ðxÞ2 Þ
1 T0 ðxÞs 2V0 ðxÞT0 ðxÞs þ c0 ðxÞ A: T0 ðxÞðC0 ðxÞ2 V0 ðxÞ2 Þ
Boundary conditions are given by qð0; sÞ ¼ q0 ðsÞ and qðX ; sÞ ¼ qX ðsÞ. In the general case, matrix As depends on x, thus this ordinary differential equation cannot be solved analytically, except for one particular case (when As does not depend on x, i.e. for uniform flow conditions). When As is constant (i.e. for uniform flow conditions), the transfer matrix is given by (see Appendix A for details): qð0; sÞ yð0; sÞ p11 ðX ; sÞ p12 ðX ; sÞ ; ð10Þ ¼ qðX ; sÞ p21 ðX ; sÞ p22 ðX ; sÞ yðX ; sÞ with p11 ðX ; sÞ ¼
k2 ek1 X k1 ek2 X ; T0 sðek2 X ek1 X Þ
ð11Þ
p12 ðX ; sÞ ¼
k1 k2 ; T0 sðek2 X ek1 X Þ
ð12Þ
p21 ðX ; sÞ ¼
ðk2 k1 Þeðk1 þk2 ÞX ; T0 sðek2 X ek1 X Þ
ð13Þ
682
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
p22 ðX ; sÞ ¼
k1 ek1 X k2 ek2 X : T0 sðek2 X ek1 X Þ
ð14Þ
This result has a very limited domain of application, since real canals are seldom in uniform regime. The authors have proposed a way to compute the solution in more realistic cases, where backwater curves occur, i.e. when matrix As ðxÞ varies with x [12]. This method is briefly explained before exposing the approximate model. 2.3.2. Numerical solution in the general case The numerical solution uses a space discretization, and the solution is obtained in each part by evaluating the transfer matrix as if the water depth was constant: this enables to use the same formulas as in the uniform case to compute the solution. The parts are then interconnected in order to get the whole transfer matrix. Let xk be a space discretization of interval ½0; X into n subintervals: 0 ¼ x0 < x1 < < xk < < xn ¼ X : This numerical solution uses the analytical solution derived above on each subinterval [xk ; xkþ1 ], considering the matrix As constant on each subinterval, which leads to the overall transition matrix: 1 Y Cs ðX ; 0Þ ¼ Cs ðxk ; xk1 Þ; ð15Þ k¼n
where matrices Cs ðxk ; xk1 Þ are obtained with the analytical method (Eq. (A.1)), using the value of coefficients in x ¼ xk . This procedure enables to compute the transition matrix Cs ðX ; 0Þ for each desired frequency s ¼ jx with a good numerical accuracy and limited computation time [12]. This complete model however relies on a large number of computations, since the backwater curve is needed at each discretization point to obtain the transfer matrix. To obtain a tractable analytical model, simple approximations of the backwater profile and of the transfer matrix are considered in the following. To get a good numerical solution, it is important to choose the number of space steps according to the variation of the water depth. In the uniform regime, only one evaluation of matrix As is necessary. In non-uniform regimes, it is necessary to evaluate this matrix at several locations. The approximations will be done in the following by considering only two locations. 3. Analytical approximate model Two different types of approximation are considered in the following: • Approximation of the backwater profile by a piecewise linear function. • Low and high frequencies approximations of the transcendental transfer function to obtain rational models with time delay. The space approximation will firstly be exposed, then the frequency approximations will be presented, leading to the global approximate model.
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
683
3.1. Backwater profile 3.1.1. Backwater profile approximation The backwater profile can be computed accurately by solving the ordinary differential equation (3), which necessitates a space discretization in many parts. In order to get simple analytical expressions, the pool is here separated in two parts: • The upstream part, in uniform flow: the water surface is approximated by a line parallel to the bottom. • The downstream part, with backwater flow: the water surface is approximated by a line tangent to the real backwater curve at x ¼ X (downstream end of the pool). Schuurmans et al. [10] used similar assumptions, but using a horizontal line for the downstream part. Using the backwater curve definition given by Eq. (3), let SX be the deviation from bed slope of the line tangent to the backwater curve at the downstream end of the pool: SX ¼
Sb Sf0 ðX Þ : 1 F02 ðX Þ
ð16Þ
Then, the intersection between the two lines occurs at abscissa x1 : 8
YX Yn < ;0 if SX 6¼ 0; max X x1 ¼ SX : X if SX ¼ 0:
ð17Þ
Then, the approximate backwater curve satisfies the equation:
Y1 for x 2 ½0; x1 ; e Y ðxÞ ¼ Y1 þ ðx x1 ÞSb for x 2 ½x1 ; X with
Y1 ¼
Yn YX XSX
if x1 6¼ 0; if x1 ¼ 0:
ð18Þ
Here, Y1 ¼ Yn when part of the pool is in uniform flow, and when the whole pool is affected by backwater, Y1 ¼ YX XSX i.e. the water depth is a straight line of slope SX . The corresponding approximation of the backwater profile is schematized in Fig. 1. 3.1.2. Computation of the uniform depth Yn To obtain an analytical model, it is necessary to have an analytical expression for the uniform depth Yn . It is not possible to have an analytical solution in the general case, but Yn can either be approximated by an analytic expression (e.g. with the hydraulic exponent method [13]), or be given by tables for a considered canal pool. The uniform depth is computed from the solution of the algebraic equation Sf0 ¼ Sb , with Sf0 the friction slope given by Manning–Strickler formula (4).
684
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
Fig. 1. Backwater profile (––) and its approximation (- - -).
In the general case, Sf0 is not a simple function of Y0 , therefore there is no analytical solution of this equation. The hydraulic exponent method [13] enables to obtain a good approximation of the pffiffiffiffiffiffi uniform depth in prismatic channels. In this method, the conveyance K0 ¼ Q0 = Sf0 is assumed to be a power function of the depth of flow Y0 : K02 ¼ CY0N ; where C is a coefficient and N is the hydraulic exponent for uniform flow computation. A good approximation of N is obtained by computing the uniform depth for two different discharges Q1 and Q2 . N is then given by N ¼2
logðQ1 =Q2 Þ : logðYn ðQ1 Þ=Yn ðQ2 ÞÞ
ð19Þ
Given the uniform flow depth Yn ðQ1 Þ corresponding to a discharge Q1 (e.g. the maximum discharge for the considered pool), the values of the uniform flow depth Yn ðQ0 Þ for other discharges Q0 can be approximated by 2=N Q0 Yn ðQ0 Þ Yn ðQ1 Þ: ð20Þ Q1 3.1.3. Interconnection of two parts To clarify the interconnection of transition matrices (15), let us consider the case of a pool separated into two parts: • The upstream part of the pool is modelled by the transfer matrix given by Eqs. (11)–(14) denoted P ðx1 ; sÞ. • The downstream part is modelled by a transfer matrix P ðX x1 ; sÞ obtained by evaluating matrix As ðxÞ at x2 2 ½x1 ; X . Let expand As ðxÞ in Taylor series around x2 2 ½x1 ; X : As ðxÞ ¼ As ðx2 Þ þ
dAs ðx2 Þ ðx x2 Þ þ Oðx x2 Þ2 : dx
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
685
Choosing x2 ¼ x1 þX enables to minimize the effect of the first order term dAdxs ðx2 Þ ðx x2 Þ on the final 2 result. The real transfer is therefore approximated by the transfer obtained using the formulas in . uniform flow, but with matrix As ðx2 Þ, where x2 ¼ x1 þX 2 In the following, pij ðsÞ denotes the terms of transfer matrix corresponding to the upstream part, pij ðsÞ those corresponding to the downstream part and ^pij ðsÞ those corresponding to the whole pool. The equivalent transfer matrix Pb ðX ; sÞ is therefore obtained by multiplication of the corresponding transition matrices (as in (15)) and rearranging in order to get the desired inputs and outputs, which leads to ^ p11 ¼ p11 þ ^ p12 ¼ ^ p21 ¼
p12 p21 ; p11 p22
p12 p12 ; p11 p22
p21 p21 ; p11 p22
^ p22 p22 ¼
p12 p21 : p11 p22
ð21Þ ð22Þ ð23Þ ð24Þ
These formulas will be used to get the overall approximation Pb ðX ; sÞ from the approximations of P ðx1 ; sÞ and P ðX x1 ; sÞ. 3.2. Frequency domain approximations 3.2.1. About approximate models at low frequencies A classical approach followed by many authors [2,3,9,10,14] is to use the moment matching method in zero (or analog method) to obtain rational approximations of the irrational transfer function. This leads to a model that fits the low frequency response of the system. In the case of the integrator delay model of Schuurmans et al., since the phase is only linked to the delay the approximation will take the group delay of the system as the delay value (the group delay of GðjxÞ is defined by o\GðjxÞ=ox [15]). However, the group delay also incorporates phase changes due to zeros and poles. This explains why the delay T computed by [10] is greater than the real delay s. Even if this approximation is mathematically valid, a delay is qualitatively completely different from stable poles or zeros from a control point of view. In fact, a delay introduces a structural limitation on the real-time performance of a controlled dynamical system [11], which is not the case of stable poles or zeros. This is why it is important to model the delay as exactly as possible, and add zeros or poles to better fit the phase shift. As shown in [12], the exact delay can be obtained by the characteristics. In the following, low and high frequencies approximations of the transfer matrix P are derived for constant flow conditions. These approximations, together with the interconnection rules (21)– (24) enable to get the overall approximate model.
686
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
3.2.2. Low frequencies approximation For low frequencies, the behavior of the transfer matrix is dominated by the integrator and the delays. The relation between the downstream (resp. upstream) discharge and the downstream level can be considered as an integrator (resp. a delayed integrator) whose gain is the variation of volume with respect to the variation of downstream level. The situation is symmetrical for the upstream level. The integrator gains will be called ‘‘equivalent backwater areas’’ since they have the dimension of an area. They can be close to the pool water surface in specific cases (e.g. at uniform flow for small pools), but are usually different from the backwater area. Variables related to the downstream water level are denoted with a subscript d, whereas variables related to the upstream water level are denoted with a subscript u. For a pool with constant water depth, the low frequency approximation reads: p11 ðsÞ ¼
1 ; Au s
p12 ðsÞ ¼ p21 ðsÞ ¼
esu s ; Au s
ð25Þ ð26Þ
esd s ; Ad s
ð27Þ
1 Ad s
ð28Þ
p22 ðsÞ ¼
with the equivalent backwater areas, obtained by computing ðspij ðsÞÞð0Þ Ad ¼
c0 T02 ðC02 V02 Þ 2 2 x 1 e T0 ðC0 V0 Þ ; c0
c T02 ðC02 V02 Þ T0 ðC20V 2 Þx Au ¼ e 0 0 1 ; c0
ð29Þ
ð30Þ
and the delays sd ¼
x ; C0 þ V 0
ð31Þ
su ¼
x : C0 V 0
ð32Þ
Using this approximation for each part of the pool together with Eqs. (21)–(24) for the interconnection, the approximation for the transfer matrix Pb ðsÞ at low frequencies is ^ p11 ðsÞ ¼
1 ; b Aus
ð33Þ
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
^ p12 ðsÞ ¼
^ p21 ðsÞ ¼
e^su s ; b us A
e^sd s ; b ds A
^ p22 ðsÞ ¼
1 : b Ads
The equivalent areas are given by b d ¼ Ad 1 þ Ad ; A Au A u b u ¼ Au 1 þ A ; Ad
687
ð34Þ
ð35Þ
ð36Þ
ð37Þ
ð38Þ
Au , Ad are given by Eqs. (29) and (30) where the variables are evaluated at 0 and where x ¼ x1 . Au , Ad are given by the same equations (29) and (30) where the variables are evaluated at x2 and where x ¼ X x1 . The equivalent delays are given by ^sd ¼ sd þ sd ;
ð39Þ
^su ¼ su þ su :
ð40Þ
The delays for the upstream uniform part sd and su are given by Eqs. (31) and (32) where the variables are evaluated at 0 and where x ¼ x1 . The delays for the downstream backwater part sd , su are given by the same equations (31) and (32) where the variables are evaluated at x2 and where x ¼ X x1 . 3.2.3. High frequencies approximation For high frequencies, the delay and the gravity waves are predominant in the transfer matrix elements. It is not easy to obtain a simple approximation of the oscillatory behavior of the gravity waves. For simplicity, only the static behavior of the oscillation is represented. This is modelled by constant gains in high frequencies, denoted pij1 . These gains are obtained by computing the average value of the modulus of each transfer for high frequencies. Indeed, limjsj!1 jpij ðsÞj may not exist, since jpij ðsÞj may oscillate at infinity. The oscillating terms are sinusoidal, therefore the average value is computed by nullifying them. For a pool with constant water depth, one then gets the approximation: p11 ðsÞ ¼ p111 ;
ð41Þ
p12 ðsÞ ¼ p121 esu s ;
ð42Þ
p21 ðsÞ ¼ p211 esd s ;
ð43Þ
p22 ðsÞ ¼ p221 :
ð44Þ
688
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
Using this approximation for each part of the pool together with Eqs. (21)–(24) for the interconnection, the transfer matrix Pb ðsÞ at high frequencies is approximated by ^ p111 ; p11 ðsÞ ^
ð45Þ
^ p121 e^su s ; p12 ðsÞ ^
ð46Þ
^ p211 e^sd s ; p21 ðsÞ ^
ð47Þ
^ p221 : p22 ðsÞ ^
ð48Þ
The gains in high frequencies are given by: ^ p111 ¼ p111 þ
p121 p211 ; p111 þ p221
ð49Þ
^ p121 ¼
p121 p121 ; p111 þ p221
ð50Þ
^ p211 ¼
p211 p211 ; p111 þ p221
ð51Þ
^ p221 ¼ p221 þ
p121 p211 p111 þ p221
ð52Þ
with
p111 ¼
1 T0 C0 ð1 F0 Þ
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u1 þ 1F0 2 eax t 1þF0 1 þ eax
p121
c0 x 2T0 ðC 2 V 2 Þ 0 0
T0 ð2 þ ðj0 1ÞF02 ÞSb A0 F0 ð1 F02 Þ
and x ¼ x1 .
ð53Þ
x
2 e pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 2 T0 C0 ð1 F0 Þ 1 þ eax vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u1 þ 1þF0 2 eax t 1F0 1 : ¼ T0 C0 ð1 þ F0 Þ 1 þ eax
In the upstream uniform part, a¼
2
2 e 2T0 ðC0 V0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ¼ T0 C0 ð1 F02 Þ 1 þ eax
p211 ¼
p221
c0 2
;
ð54Þ
ð55Þ
ð56Þ
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
689
For the downstream backwater part, the same expressions are valid, replacing x with x ¼ X x1 , a with T0 A0 dT0 2 2 a¼ ð2 þ ðj0 1ÞF0 ÞSb 2 þ ðj0 1ÞF0 þ j0 2 F04 SX 2 2 A0 F0 ð1 F0 Þ T0 dY and all the variables by their value taken at x2 . 3.3. Global approximate model The global model is simply obtained by adding the approximations in high and low frequencies: 1 þ^ p111 ; b us A
^ p11 ðsÞ ¼
ð57Þ
1 ^ þ^ p121 e^su s ; p12 ðsÞ ¼ b us A
ð58Þ
^ p21 ðsÞ ¼
1 þ^ p211 e^sd s ; b ds A
ð59Þ
^ p22 ðsÞ ¼
1 ^ p221 : b ds A
ð60Þ
The approximate model is now completely presented. Given a prismatic canal with known geometry and Manning coefficient, a discharge Q0 and a downstream limit condition YX , Eqs. (16)–(20) enable to compute the integrator delay zero model parameters: the delays ^sd and ^su given b u , given by Eqs. (37) and (38) and b d and A by Eqs. (39) and (40), the equivalent backwater areas A the gains at high frequencies ^ pij1 , given by Eqs. (49)–(52). The accuracy of this approximate model will now be evaluated.
4. Validation of the approximate model 4.1. Example canals The paper is illustrated with two trapezoidal prismatic channels, with different characteristics (see Table 1, where X is the channel length (m), m the bank slope, B the bed width (m), Sb the bed
Table 1 Parameters for the two example canals Canal 1 Canal 2
X
m
B
Sb
n
Yn
Qmax
Qmin
3000 6000
1.5 1.5
7 8
0.0001 0.0008
0.02 0.02
2.12 2.92
14 80
1.75 10
690
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
slope, n the roughness coefficient (s m1=3 ), Yn the uniform depth (m) corresponding to the maximum discharge Qmax (m3 /s)). Canal 1 is a short oscillating canal, and canal 2 is a long delayed canal. The downstream limit condition YX is chosen equal to the uniform depth Yn for Qmax . The uniform depth approximations are obtained using Eqs. (19) and (20) where Q1 ¼ Qmax and Q2 ¼ Qmax =8. This leads to N ¼ 3:51 for canal 1 and N ¼ 3:56 for canal 2. The obtained approximations of the uniform depths are very accurate for both canals (see Fig. 2). The approximate backwater curves obtained in this case for an average discharge Qmax =2 are depicted in Fig. 3 for both canals.
Fig. 2. Variation of the uniform depth with the discharge, canals 1 and 2. Comparison between the exact (––) and approximate (- - -) values.
Fig. 3. Backwater curve (––) and its approximation (- - -) for Qmax =2, canals 1 and 2.
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
691
4.2. Model evaluation 4.2.1. Comparison of frequency responses The Bode plots of canals 1 and 2 are obtained in non-uniform flow conditions for Q0 ¼ Qmax =2 using the complete numerical method [12] (see Fig. 4). One can notice the behavior of both canals, which is completely different: canal 1 is oscillating, with resonant modes, and canal 2 exhibits a long time delay. The approximate model reproduces very well the global frequency behavior of the system, especially in the case of canal 1. Canal 2 has a pole that is not taken into account by the approximate model. In order to show the accuracy of the method, it is important to see how the parameters vary with the hydraulic conditions, i.e. the discharge Q0 . 4.2.2. Parameters comparison Let us consider a constant downstream limit condition YX corresponding to the uniform flow for Qmax , and vary the discharge Q0 . The canal is therefore in uniform flow for Q0 ¼ Qmax and the backwater curve increases as the discharge diminishes. The approximate parameters are compared with their exact values when they are available. b d , ^sd and lim supx!1 jp21 ðjxÞj Figs. 5–7 respectively depict the exact and approximate values of A and lim inf x!1 jp21 ðjxÞj for both canals. It is clear from the figures that the approximation is very accurate for the delay ^sd and for b d . The gain at high frequencies ^ p211 is rather well estimated for canal 2, and is in coefficient A between min and max values for canal 1. In fact, the max value of the gain at high frequencies ^p211 increases a lot as the discharge tends towards zero, since the damping of resonant modes diminishes. This is not taken into account by the approximate model, which is sufficient for control purposes, since the resonant modes will not be controlled. These resonant modes are usually filtered in order not to destabilize the controlled system.
Fig. 4. Bode plot of transfer p21 ðsÞ (––) and its approximation e^sd s
1 bA d s
þ^ p211 (- - -), canals 1 and 2 for Q0 ¼ Qmax =2.
692
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
b d with the discharge, for a downstream limit condition YX corresponding to the uniFig. 5. Variation of coefficient A form flow for Qmax , canals 1 and 2. Comparison between the exact (––) and approximate (- - -) values.
Fig. 6. Variation of delay ^sd with the discharge, for a downstream limit condition YX corresponding to the uniform flow for Qmax , canals 1 and 2. Comparison between the exact (––) and approximate (- - -) values.
4.2.3. Comparison of time-domain simulations The approximate model is compared to a finite dimensional model obtained by rational approximation of the exact model [12]. Simulations are done on Matlab for canals 1 and 2, for Q0 ¼ Qmax =2 (see Fig. 8). The models are connected with a downstream gate (equivalent to a local feedback) with a gain kv ¼ 6:84. The simulation is very accurate for both canals: in fact, as can be seen from the Bode plot, the whole frequency response is very well fitted by the approximate model, except the resonant modes
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
693
Fig. 7. Variation of gain ^p211 with the discharge, for a downstream limit condition YX corresponding to the uniform flow for Qmax , canals 1 and 2. Comparison between min and max values (––) and the approximation (- - -).
Fig. 8. Step responses of the rational approximation of the complete model (––) and IDZ model (- - -), canals 1 and 2.
for canal 1 and a supplementary pole for canal 2. These discrepancies occur in rather high frequencies and at a small gain, this is why they are not clearly visible in the step response.
5. Conclusion The paper proposes a new analytical approximate model for an open-channel pool subject to backwater. This integrator delay zero model is simple, and very accurate. It is obtained from rigorous mathematical development on linearized Saint-Venant equations. It models the complete transfer matrix (i.e. the downstream influence is also taken into account). The inputs are the
694
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
upstream and downstream discharges, and its outputs are the upstream and downstream water levels. The model has been validated in frequency and time domain using the complete transfer matrix obtained following [12] on two different example canals. The approximate model can be used for controller design, such as distant downstream PI controllers, or even for advanced controller design (e.g. multivariable controllers).
Acknowledgements This work was partially supported by the joint research program INRA/Cemagref ASS AQUAE no 02 on the control of delayed hydraulic systems.
Appendix A. Saint-Venant transfer matrix: analytic solution for uniform flow In uniform conditions, since matrix As does no longer depend on x, the ordinary differential equation (9) can be solved analytically. To do this, let us diagonalize matrix As As ¼ Ps Ds Ps1 ; with
Ds ¼
k1 ðsÞ 0
0 ; k2 ðsÞ
Ps ¼
T0 s k1 ðsÞ
T0 s ; k2 ðsÞ
Ps1
1 ¼ T0 sðk1 ðsÞ k2 ðsÞÞ
k2 ðsÞ k1 ðsÞ
T0 s : T0 s
The eigenvalues are given by k1;2 ðsÞ ¼
2T0 V0 s þ c0
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4C02 T02 s2 þ 4T0 ðV0 c0 ðC02 V02 ÞT0 b0 Þs þ c20 ; 2T0 ðC02 V02 Þ
Once the matrix As is diagonalized, the ordinary differential equation (9) can be solved analytically: qðx; sÞ qð0; sÞ ¼ Cs ðx; 0Þ ðA:1Þ yðx; sÞ yð0; sÞ with 0
Cs ðx; 0Þ :¼ Ps eDs x Ps1
k1 ek2 x k2 ek1 x B k1 k2 B ¼B @ k1 k2 ðek1 x ek2 x Þ T0 sðk1 k2 Þ
1 T0 sðek2 x ek1 x Þ C k1 k2 C C: k1 x k2 x A k1 e k2 e k1 k2
This is the classical result already obtained by many authors [2,3,5,6,10]. Reordering the result in order to have the initial boundary conditions as inputs (qð0; sÞ and qðX ; sÞ) leads to the Saint-Venant transfer matrix:
X. Litrico, V. Fromion / Appl. Math. Modelling 28 (2004) 677–695
yð0; sÞ yðX ; sÞ
!
0
k2 ek1 X k1 ek2 X B T sðek2 X ek1 X Þ B 0 ¼B @ ðk2 k1 Þeðk1 þk2 ÞX
1 k1 k2 T0 sðek2 X ek1 X Þ C C C k1 X k2 X A k1 e k2 e
qð0; sÞ qðX ; sÞ
695
! :
ðA:2Þ
T0 sðek2 X ek1 X Þ T0 sðek2 X ek1 X Þ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} P ðX ; sÞ ¼ ðpij ðX ; sÞÞ
References [1] M. Shand, Automatic Downstream Control Systems for Irrigation Canals, Ph.D. Thesis, University of California, Berkeley, 1971. [2] G. Corriga, D. Salembeni, S. Sanna, G. Usai, A control method for speeding up response of hydroelectric stations power canals, Appl. Math. Modell. 12 (1988) 627–633. [3] Y. Ermolin, Study of open-channel dynamics as controlled process, J. Hydraul. Eng. 118 (1) (1992) 59–71. [4] S. Hancu, P. Dan, Wave-motion stability in canals with automatic controllers, J. Hydraul. Eng. 118 (12) (1992) 1621–1638. [5] G. Corriga, A. Fanni, S. Sanna, G. Usai, A constant-volume control method for open channel operation, Int. J. Modell. Simulat. 2 (2) (1982) 108–112. [6] J. Baume, J. Sau, P.-O. Malaterre, Modeling of irrigation channel dynamics for controller design, in: Conf. on Systems, Man and Cybernetics, SMCÕ98, San Diego, CA, 1998, pp. 3856–3861. [7] C. Seatzu, Design and robustness analysis of decentralized constant volume-control for open-channels, Appl. Math. Modell. 23 (6) (1999) 479–500. [8] O. Balogun, Design of Real-Time Feedback Control for Canal Systems Using Linear Quadratic Regulator Theory, Ph.D. Thesis, Department of Mechanical Engineering, University of California at Davis, 1985, 230 p. [9] P.-O. Malaterre, Modelisation, Analyse et Commande Optimale LQR dÕun Canal dÕIrrigation, Ph.D. Thesis, ISBN 2-85362-368-8, Etude EEE no. 14, LAAS, CNRS, ENGREF, Cemagref (in French), 1994. [10] J. Schuurmans, O. Bosgra, R. Brouwer, Open-channel flow model approximation for controller design, Appl. Math. Modell. 19 (1995) 525–530. [11] K. Astr€ om, Limitations on control system performance, Eur. J. Control 6 (2000) 1–19. [12] X. Litrico, V. Fromion, Infinite dimensional modelling of open-channel hydraulic systems for control purposes, in: 41st Conf. on Decision and Control, Las Vegas, 2002, pp. 1681–1686. [13] V. Chow, Open-channel Hydraulics, McGraw-Hill Book Company, New York, 1988, 680 p. [14] X. Litrico, D. Georges, Robust continuous-time and discrete-time flow control of a dam–river system (I): Modelling, Appl. Math. Modell. 23 (11) (1999) 809–827. [15] J. Freudenberg, D. Looze, A sensitivity tradeoff for plants with time delay, IEEE Trans. Automat. Control 32 (2) (1987) 99–104.