Nonlinearity Tests of the Saint Venant Equations
52nd IEEE Conference on Decision and Control December 10-13, 2013. Florence, Italy
Nonlinearity tests of the Saint Venant equations Mathias Foo and Erik Weyer
Abstract— The Saint Venant equations are two nonlinear partial differential equations (PDE) which are used to describe the dynamics of one-dimensional flow in open water channels. Despite being nonlinear PDEs, the Saint Venant equations seem to exhibit linear behaviour in response to sinusoidal input signals. It is therefore of interest to determine ”how nonlinear” the Saint Venant equations are. In this paper, we investigate the nonlinearity in the Saint Venant equations using several commonly used nonlinearity tests suggested in the literature. Five different open water channels are considered, and the results from the nonlinearity tests show that the Saint Venant equations are nearly linear in an operating region from at least half the nominal flow to twice the nominal flows, and many of the channels display linear behaviour in a larger operating region. This finding is useful as it further justifies the use of linear control design methodologies for open water channels.
I. I NTRODUCTION Traditionally, the dynamics of open water channels1 are modelled by two nonlinear partial differential equations (PDE) which are known as the Saint Venant equations. The Saint Venant equations are derived from mass and momentum balance equations for one-dimensional flow and assuming no contribution from lateral flows, they are given by [1], ∂A ∂Q + =0 ∂t ∂x ∂Q ∂y ∂ Q2 + gA + + gA(Sf − S0 ) = 0 ∂t ∂x A ∂x
(1)
where Q is the flow, A is the wetted cross sectional area, T is the top width, g is the gravity constant, S0 is the bottom slope and Sf is the friction slope. The friction slope is a nonlinear function given by n2 Q2 P 4/3 A−10/3 , where P is the wetted perimeter and n is the Manning friction coefficient, which encapsulates the effect of flow resistance This work was supported by The Farms Rivers and Markets Project, an initiative of Uniwater and funded by the National Water Commission, the Victorian Water Trust, The Dookie Farms 2000 Trust (Tallis Trust) and the University of Melbourne and supported by the Departments of Sustainability and Environment and Primary Industry, the Goulburn Broken Catchment Management Authority and Goulburn-Murray Water. The first author also gratefully acknowledge the financial support from National ICT Australia (NICTA). NICTA is funded by the Australian Government as represented by the Department of Broadband, Communications and the Digital Economy and the Australian Research Council through the ICT Centre of Excellence program. M. Foo is with Asia Pacific Center for Theoretical Phyiscs (APCTP), Hogil Kim Memorial Building, #501, Pohang University of Science and Technology (POSTECH), San 31, Hyoja-dong, Nam-gu, Pohang, Gyeongbuk, 790-784, Korea [email protected] E. Weyer is with the Department of Electrical and Electronic Engineering, The University of Melbourne, 3010, Victoria, Australia
[email protected] 1 Open water channels are system where the water has a free surface to the air. Examples of open water channel include rivers and irrigation channels.
978-1-4673-5716-6/13/$31.00 ©2013 IEEE
and channel roughness. Typically, the cross sectional area is assumed to be a trapezoid,√and hence, A = (b + sy)y, T = b + 2sy and P = b + 2y 1 + s2 where b, s and y are the bottom width, side slope and water level respectively. As there is no closed form solution, the Saint Venant equations have to be solved numerically and the Preissmann scheme which is a finite difference method is commonly used to discretise the Saint Venant equations in time and space (see e.g. [1] - [3]). These equations are solved together with the boundary equations which gives the flows at the upstream and downstream end. For the downstream boundary condition, a sharp crested weir is used in this paper, whereby the flow and water level are related by Q = c(y − p)3/2 , where c is a constant, p is the height of the weir, and y > p is the water level (see Fig. 1). The initial condition is usually the steady state solution of Eqn. (1). The solutions of the Saint Venant equations are the flows and the water levels in space and in time. Despite being nonlinear PDEs, the
Fig. 1: Left: Sharp crested weir. Right: Schematic side view (not to scale). Saint Venant equations seem to exhibit linear behaviour in response to a sinusoidal input. When applying a sinusoidal input flow at the upstream end to the Saint Venant equations, we observe that the output, i.e. the downstream flow, is also a sinusoid of the same frequency but with a phase shift (see Fig. 2). No harmonic or subharmonic can be visually detected in the output. Thus, a natural suggestion is that the Saint Venant equations exhibit linear behaviour within an operating region. In this paper, we investigate ”how nonlinear” the Saint Venant equations are by applying some commonly used nonlinearity tests suggested in the literature. If the Saint Venant equations are almost linear within an operating region, then it is reasonable to conduct a frequency response analysis using the Saint Venant equations to obtain a Bode plot, which can be used for control design and in general it also justifies the use of linear control methodologies. This paper is organised as follows. In Section II a brief overview of nonlinearity tests is given. In Section III the
4925
Flow
Upstream flow Downstream flow
4
If the system is linear, Rx¯y (τ ) = 0, ∀τ and this test is claimed to be robust against noise [6].
m3/s
3.5 3
C. Higher order autocorrelation test
2.5 2 0.6
0.8
1
1.2 1.4 Time (minutes)
1.6
1.8
2 4
x 10
Fig. 2: Input and output response of the Saint Venant equations given a sinusoidal input at frequency, ω0 = 0.001 rad/min.
nonlinearity tests are applied to the Saint Venant equations. A frequency response analysis of open channel system is presented in Section IV, and conclusions are given in Section V. II. OVERVIEW
Ry¯2 y¯(τ ) = E{¯ y 2 (t − τ )¯ y (t)}
(5)
is computed with y¯(t) given by Eqn. (4). If the system is linear, Ry¯2 y¯(τ ) = 0, ∀τ . This test method has become very popular as it is simple and easy to apply [10]. Despite its popularity, this test method does raise some questions since while Ry2 y (τ ) = 0, ∀τ is necessary, it is not sufficient, and the test method fails for a simple cubic relationship (e.g. y(t) = x3 (t)) [11].
ON NONLINEARITY TESTS
An extensive list of works on nonlinearity tests can be found in [4]-[8]. Some of the commonly used tests include time domain tests, nonlinear cross-correlation tests, higher order autocorrelation tests and tests based on broadband periodic signal excitation. We will discuss some of the advantages and disadvantages of some of these tests. A. Time domain test The simplest and most popular test is the time domain test [8]. The input signal, uin,1 (t) is applied to the system and its corresponding output, uout,1 (t) is obtained. Then, the input signal, uin,1 (t) is scaled by a factor α1 (i.e. uin,1 (t) is replaced by α1 uin,1 (t) = uin,2 (t)). The signal uin,2 (t) is then applied to the system and its corresponding output, uout,2 (t) is obtained. If the system is linear, the output will be scaled by the same factor, i.e. uout,2 (t)/uout,1 (t) = α2 (t) = α1 . Although this test is simple, it is not widely used in practice as at least two experiments must be carried out. The choice of α1 may be limited by process constraints. If a small value of α1 is applied, the effect of noise is more significant, jeopardising the overall test. B. Nonlinear cross-correlation test A nonlinear cross-correlation test was introduced in [9]. In this test, the input signal is a nonzero mean white Gaussian noise or a pseudorandom signal. This signal is symmetrically distributed around the mean and has zero odd order moment and nonzero even order moment [6]. Introducing the term x(t) =
The higher order autocorrelation test in [7] is another nonlinearity test. Using a similar input signal as in the nonlinear cross-correlation test but this time with zero mean, the autocorrelation given by
[¯ u2 (t) − E{¯ u2 (t)}] 2 σ(¯ u (t))
(2)
where σ(.) denotes the standard deviation and u ¯(t) = [u(t)−E{u(t)}] , the cross-correlation σ(u(t)) Rx¯y (τ ) = E{x(t − τ )¯ y (t)}
(3)
is computed, where y¯(t) =
[y(t) − E{y(t)}] σ(y(t))
(4)
D. Broadband periodic excitation A more recent nonlinearity test suggested in [8] is based on using a broadband periodic excitation input signal. The idea originated from [12], where a signal that contains only certain harmonics is used to excite the system. One example of such a signal is the odd-odd multisines signal. The frequencies that are excited by the odd-odd multisines signal are (4k + 1)f0 , k = 0, 1, . . . , F , where F + 1 is the number of frequencies and f0 is the fundamental frequency. If the frequencies, (4k + 2)f0 , (4k + 3)f0 and (4k + 4)f0 are excited, then the system is nonlinear. III. N ONLINEARITY TESTS FOR OPEN CHANNEL SYSTEMS In this section, the nonlinearity tests presented in Section II are applied to five different open water channels. The first one is the river reach from Casey’s to Gowangardie Weir (Reach CG) in the Broken River, Australia. The second one is the Haughton Main Channel, Pool 9 (HMCP9) in Australia. The third and fourth one are the irrigation test channels 1 and 5, (ITC1 and ITC5), taken from [13]. The fifth one is the Test Canal 1, Pool 3 (TC1P3) taken from [14]. The Saint Venant equations for both Reach CG and HMCP9, have been calibrated and validated against measured data (see [15] and [16] respectively). ITC1 and ITC5, have been extensively used as benchmark test channels in several works (see e.g. [3]). TC1P3 is of interest since it is a short and steep channel. The physical parameters for all the five open channels are summarised in Table I. In all these tests, we consider the relationship between upstream flow and downstream flow. A. Time domain test The input signal is a pseudorandom binary signal (PRBS) as this input signal can excite a broad range of the spectrum enabling collection of all the required information from a single measurement [8]. The levels of the PRBS are 2.50m3/s and 3.50m3 /s, while the periods are chosen to be three times the time constant obtained from a step response, i.e. 960 minutes for Reach CG, 72 minutes for HMCP9, 216 minutes
4926
TABLE I: Open channel parameters. L is the length of the channel, b is the bottom width, s is the side slope, S0 is the bottom slope, n is the Manning friction coefficient. Parameters L b s S0 n
Reach CG 26700m 9-12m 2-3 0.0008-0.0020 0.060-0.085
HMCP9 853m 6m 2 0.0002 0.02
ITC1 3000m 7m 1.5 0.0001 0.02
ITC5 6000m 8m 1.5 0.0008 0.02
TC1P3 400m 1m 1.5 0.002 0.014
for ITC1, 225 minutes for ITC5 and 12 minutes for TC1P3. The values of α1 were 0.10, 0.20, 0.25, 0.33, 0.50, 2.00, 3.00, 4.00, 5.00 and 10.00.
B. Nonlinear cross-correlation test The input signal is a PRBS signal with levels 2.50m3/s and 3.50m3/s. This signal was applied to all five open water channels. The cross-correlation is computed using Eqn. (3) and they are shown in Fig. 3. In the plot, we included the threshold value of ±0.05 from Rx¯y (τ ) = 0, where the system is considered linear if the cross-correlation stays within this threshold value. We observe that the cross-correlation stays within the threshold for all τ . We repeated the test for several mean flows2 (i.e. mean flow from 2.00ms /s to 6.00ms/s) and also retaining the same mean flow of 3.00m3/s but with different levels for the PRBS signal (i.e. 2.90m3/s and 3.10m3/s, 2.45m3/s and 3.55m3/s to 2.00m3/s and 4.00m3/s). All tests yield similar results. The nonlinear cross-correlation
In our analysis we have allowed α2 (t) to deviate with 5% from α1 . This 5% percent tolerance bound is motivated from the definition of settling time of the step response; the settling time is defined as the time from when the output response stays within ±5% of its steady state value. Table II quantifies the percentage of the time for which α2 (t) exceeds the ±5% bound.
Rxy(τ)
Rxy(τ)
Reach CG 0.05 0 −0.05
ITC5 19.27 2.92 0.12 0.00 0.00 0.00 0.00 0.00 0.00 0.01
TC1P3 26.27 14.19 11.08 6.84 0.00 0.00 4.26 7.43 7.77 14.90
Rxy(τ)
ITC1 26.78 11.40 6.06 1.62 0.00 0.00 0.00 0.98 1.53 3.76
Rxy(τ)
HMCP9 24.04 9.63 5.19 1.45 0.00 0.00 0.00 0.16 1.37 2.57
Rxy(τ)
Reach CG 0.28 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
The results shown in Table II reveal several findings. For Reach CG, which is the longest one considered, the Saint Venant equations are linear for all considered α1 except for α1 = 0.10 but the percentage of time the bound was exceeded is very small. For ITC5, which is the second longest, we can see similar trend where the Saint Venant equations are linear for most of the considered α1 . As the open water channels get shorter and steeper, (i.e. from ITC1 to HMCP9 to TC1P3), we observe that the ’linear’ operating region gets smaller. One can note that α2 (t) exceeds the ±5% bound more often for low values of α1 , i.e. under low flow condition. One reason for the larger deviation at smaller α1 is that the resonances tend to show up more clearly at low flow. Based on the results from time domain test, we can say that the Saint Venant equations are approximately linear within an operating region ranging from at least half to two times the base flow, and for most of the channels the linear operating region is even larger.
100
200
300
400
500 600 τ HMCP9
700
800
900
1000
0
100
200
300
400
500 τ ICT1
600
700
800
900
1000
0
100
200
300
400
500 τ ICT5
600
700
800
900
1000
0
100
200
300
400
500 600 τ TC1P3
700
800
900
1000
0
100
200
300
400
700
800
900
1000
0.05 0 −0.05
TABLE II: Percentage of α2 (t) exceeding the ±5% bound. α1 0.10 0.20 0.25 0.33 0.50 2.00 3.00 4.00 5.00 10.00
0
0.05 0 −0.05
0.05 0 −0.05
0.05 0 −0.05 500 τ
600
Fig. 3: Nonlinear cross-correlation test plot with mean flow 3.00m3/s. Dashed line: Threshold. test show that all open water channels behave quite linearly for different mean flows as well as different flow levels. C. Higher order autocorrelation test The same input signal as in the nonlinear cross-correlation test is used but now with zero mean. The autocorrelation was computed using Eqn. (5). Note that a zero mean input signal means that we are dealing with negative flows, which is not physically feasible. Thus, for this test, a nonzero mean input signal is applied to the Saint Venant equations and then we remove the mean at the output of the Saint Venant equations 2 The difference between the two levels of PRBS signal remains the same but with different mean flow levels. E.g. mean flow of 2.00m3 /s has the two PRBS signal level of 1.50m3 /s and 2.50m3 /s.
4927
prior to computing the autocorrelation. The autocorrelation plots for all five open water channels are shown in Fig. 4. Using the same threshold value of ±0.05 as in Section
R
2
y y
(τ)
Reach CG 0.05 0 −0.05 0
100
200
300
400
500 τ
600
700
800
900
1000
600
700
800
900
1000
R
2
y y
(τ)
HMCP9 0.05 0 −0.05 0
100
200
300
400
500 τ
R
2
y y
(τ)
ICT1 0.05 0 −0.05 0
100
200
300
400
500 τ
600
700
800
900
1000
600
700
800
900
1000
where MNH is the magnitude of nonexcited harmonics and MEH is the magnitude of excited harmonics. Nne and Ne are the number of nonexcited and excited harmonics respectively. For illustration, the spectrum of the input and output signals for Reach CG and ITC5 are shown in Fig. 5.
0.05 0 −0.05 0
100
200
300
400
500 τ
Reach CG 0.015
Magnitude
2
y y
(τ)
ICT5
R
As the characteristics of the five open water channels are different, we considered frequencies up to where the magnitude of the output spectrum drops 20dB from the magnitude of the input spectrum. For quantitative measure, we calculate the N2E Ratio, which we define as the ratio of the normalised sum of the magnitude of the nonexcited harmonics to the normalised sum of the magnitude of the excited harmonics for each of the open channels up to the frequencies where the magnitude of the output spectrum drops 20dB from the magnitude of the input spectrum. The N2E Ratio is given by, PNne 1 MNH 1 Nne N2E Ratio = 1 PNe (7) 1 MEH Ne
R
2
y y
(τ)
TC1P3 0.05 0 −0.05
0.01
0.005
0
0
100
200
300
400
500 τ
600
700
800
900
1000
1
2
3
4 5 6 Frequency (rad/min)
7
8
9 −3
x 10
ITC5 0.015
III-B, we can see that, the autocorrelation stay within the threshold for all τ . We repeat the analyses using different mean flows and different flow levels for the PRBS signal (just as we did for the nonlinear cross-correlation test) and similar results were observed. The higher order autocorrelation tests show that all open water channels behave quite linearly for different mean flows and different flow levels. D. Broadband periodic excitation The odd-odd multisines signal is given by u(t) =
F X
[AB sin(ω0 (4k + 1)t + φk )] + BB
(6)
k=1
where AB is the amplitude, BB is a constant offset to ensure positive flow, F is the number of frequencies and φk = −k(k − 1)π/F are Schroeder phases [17]. They are constructed such that the peak to peak amplitude of the sum of the sinusoid signals is minimised (see e.g. [8] and [17]). The odd-odd multisines are applied as the input to the Saint Venant equations for all five open water channels with F in the range3 from 20 to 70000, AB = 0.025, BB = 3, and ω0 = 0.0001 rad/min. 3F
takes different values due to the different characteristics of the analysed channels.
Magnitude
Fig. 4: Higher order autocorrelation test plot with mean flow 3.00m3/s. Dashed line: Threshold.
0.01
0.005
0
0.01
0.02
0.03
0.04 0.05 0.06 Frequency (rad/min)
0.07
0.08
0.09
Fig. 5: Spectrum of odd-odd multisines signal for Reach CG and ITC5. -x: Input, -o: Output. From Fig. 5 we observe that for ITC5, the output spectrums have the same excited frequencies as the input spectrum. At the other frequencies, the output spectrum is not excited as much, (at least visually). For Reach CG, the output spectrum is excited at the same frequencies as present in the input spectrum. However, some other frequencies are also present. The magnitude of those frequencies are relatively small compared to the magnitude of the frequencies present in the input signal. The N2E Ratios for the five open water channels are given in Table III. From the N2E Ratios, we can see that all the TABLE III: N2E Ratio for the five open water channels. N2E Ratio
Reach CG 0.185
HMCP9 0.016
ITC1 0.028
ITC5 0.054
TC1P3 0.013
values are small. The largest ratio is from Reach CG but it is less than 0.20. The test is repeated with larger values of AB and similar N2E Ratios as in Table III are obtained.
4928
Reach CG
RESPONSE ANALYSIS
0
Magnitude (dB)
IV. F REQUENCY
All the nonlinearity tests presented above suggest that the Saint Venant equations are linear within an operating region. This indicates that a Bode plot can be constructed using a sine sweep4 .
−10
ω
3dB
−20 −30 −40 −3 10
−2
−1
10 rad/min
10
0
Phase (°)
A. Sine sweeping method Consider a sinusoidal input, ui (t) = AS sin(ω0 t)
−2000 −3 10
(8)
Mc (N ) φˆ = tan−1 Ms (N )
HMCP9
Magnitude (dB)
0
−5 ω
3dB
−10 −3 10
−2
10 rad/min
−1
10
0
−50
−100 −3 10
−2
10 rad/min
−1
10
(b) HMCP9
Fig. 6: Bode plot relating upstream flow to downstream flow for Reach CG and HMCP9
B. Results and discussions
(11)
where N is the number of data points. |G(jω0 )| and φ = ∠G(jω0 ) can be estimated as p Mc2 (N ) + Ms2 (N ) 2 ˆ (12) |G(jω 0 )| = AS and
−1
10
(10)
and N 1 X yo (t) sin(ω0 t) Ms (N ) = N t=1
−2
10 rad/min
(a) Reach CG
Phase (°)
yo (t) = BS sin(ω0 t + φ) + D(t) + transient + nonlinearities (9) where BS = AS |G(jω0 )|, φ = ∠G(jω0 ) = Im[G(jω0 )] , D(t) is the disturbances and G(jω) is the tan−1 Re [G(jω0 )] transfer function relating the input and the output. The effect of D(t) can be reduced by using a correlation method [18]. The transient effects can be reduced by not considering the initial part of the data. Defining N 1 X yo (t) cos(ω0 t) N t=1
−1000 −1500
where AS is the amplitude and ω0 is the frequency. If the system is linear and time invariant, the output response is a sinusoidal signal with the same frequency, scaled amplitude, and a phase shift. In the real system, the output response, yo (t) to the sinusoidal input also comprises transient effects, the effect of nonlinearities and disturbances.
Mc (N ) =
−500
(13)
(see [18] and [19] for details.). Having derived the expression for the magnitude and the phase of the transfer function, we can now obtain the frequency response analysis of all the five open water channels. The sinusoidal input flow, ui (t) = 0.75 sin(ω0 t) + 3.0 m3 /s is applied to the Saint Venant equations for all five open water channels and the experiment is repeated for several frequencies, ω0 in the frequency range of interest. By obtaining the magnitude and phase at each frequencies using Eqns. (12) and (13), the Bode plot can be constructed. 4 A sine sweep is of course not possible on real open channel system, but here we only consider the Saint Venant equations so the practical difficulties in applying a sine sweep test is not an issue.
The Bode plot for all five open water channels are shown in Figs. 6 and 7. The dashed line shown in the Bode plot for the amplitude and phase represents -3dB and -180◦ respectively. From the Bode plot of Reach CG (Fig. 6(a)), we can see that the frequency range relevant for control is similar to the frequency response of a pure time delay system. The 3dB bandwidth of the system is approximately 0.0032 rad/min but we note that the phase shift is already more than -180◦ at 0.0023 rad/min indicating a dominant time delay. At higher frequency range, we observed two small peaks at around 0.04 rad/min and 0.08 rad/min, which could be due to the effect of waves. For HMCP9, the Bode plot is shown in Fig. 6(b). In the frequency range relevant for control, the frequency response is similar to a first order system. The 3dB bandwidth of the system is approximately 0.035 rad/min, and at this frequency, we note that the phase is approximately -45◦ , which according to the theory of linear system is in agreement with a first order system. For ITC1, the Bode plot is shown in Fig. 7(a). Like HMCP9, in the frequency range relevant for control, the frequency response can also be approximated by a first order system. The 3dB bandwidth of the system is approximately 0.012 rad/min and at this frequency the phase is also approximately -45◦ , indicating a first order system. From the Bode plot of ITC5 (Fig. 7(b)), the 3dB band-
4929
ITC1
V. C ONCLUSIONS
Magnitude (dB)
0 −5
We applied some of the commonly used nonlinearity tests to the Saint Venant equations to determine how ”nonlinear” the Saint Venant equations are for a variety of channels. The results from the tests suggest that the Saint Venant equations are linear within a flow regime between at least half and two times the nominal flow, and for most channels this flow regime is even larger. Based on these findings, it is reasonable to construct Bode plots using sine sweeping methods, and the results further justify the use of linear control design methods for open water channels.
ω3dB
−10 −15 −20 −3 10
−2
−1
10 rad/min
10
0
Phase (°)
−50
−100
−150 −3 10
−2
−1
10 rad/min
10
(a) ITC1
R EFERENCES
ITC5
Magnitude (dB)
0 −5
ω3dB
−10 −15 −20 −25 −3 10
−2
−1
10 rad/min
10
0
Phase (°)
−100 −200 −300 −400 −500 −3 10
−2
−1
10 rad/min
10
(b) ITC5 TC1P3
Magnitude (dB)
0 −20 ω
3dB
−40 −60 −80 −3 10
−2
10
−1
10 rad/min
0
10
1
10
0
Phase (°)
−100 −200 −300 −400 −500 −3 10
−2
10
−1
10 rad/min
0
10
1
10
(c) TC1P3
Fig. 7: Bode plot relating upstream flow to downstream flow for ITC1, ITC5 and TC1P3.
width of the system is approximately 0.013 rad/min and at this frequency, the phase is approximately -95◦ . Moreover, at magnitude of -6dB, the phase is approximately -150◦, suggesting that the frequency range relevant for control is similar to the frequency response of a second order system with delay. For TC1P3, the Bode plot is shown in Fig. 7(c). Like the case of HMCP9 and ITC1, the frequency response is similar to a first order system where the 3dB bandwidth of the system is approximately 0.20 rad/min and the phase at this frequency is approximately -45◦ . We note that at higher frequency range, we observe peaks at around 5.00rad/min, which could be due to the effect of waves.
[1] M. Chaudhry, Open-channel Flow, 2nd Edition, Prentice Hall, Englewood Cliffs, NJ; 2007. [2] V.T. Chow, D.R. Maidment and L.W. Mays, Applied Hydrology, McGraw-Hill Series in Water Resources and Environmental Engineering; 1988. [3] X. Litrico and V. Fromion, Modeling and Control of Hydrosystems, Springer; 2009. [4] H.G. Natke, J-N. Juang, and W. Gawronski, ”A brief review on the identification of nonlinear mechanical system”, Proceedings of 6th International Modal Analysis Conference, Kissimee, USA, 1988. [5] R. Haber, ”Nonlinearity test for dynamic processes”, Proceedings of 7th IFAC/IFORS Symposium on Identification and Parameter Estimation, York, UK, 1985. [6] R. Haber and L. Keviczky, Nonlinear system identification - InputOutput Modeling Approach, Volume 2: Nonlinear System Structure Identification, Kluwer Academic Publishers; 1999. [7] S.A. Billings and W.S.F. Voon, ”Structure detection and model validity tests in the identification of nonlinear systems”, IEEE Proceedings, vol. 130, no. 4, pp 193-199, 1983. [8] R. Pintelon and J. Schoukens, System identification - A Frequency Domain Approach, 2nd Edition, IEEE Press; 2012. [9] R. Haber, ”Parametric identification of nonlinear dynamic systems based on correlation functions”, Proceedings of 5th IFAC Symposium on Identification and Parameter Estimation, Darmstadt, Germany, 1979. [10] P. Varlaki, G. Terdik and V.A. Lototsky, ”Tests for linearity and bilinearity of dynamic systems”, Proceedings of 7th IFAC/IFORS Symposium on Identification and Parameter Estimation, York, UK, 1985. [11] M.J. Korenberg and I.W. Hunter, ”The identification of nonlinear biological systems: Wiener kernel approaches”, Annuals of Biomedical Engineering, vol. 18, pp 629-654, 1990. [12] D.C. Evans, D. Rees and D.L. Jones, ”Identifying linear models of systems suffering nonlinear distortions”, Proceedings of IEE Control 94, Coventry, UK, 1994. [13] J-P. Baume, J. Sau and P-O. Malaterre, ”Modelling of irrigation channel dynamics for controller design”, Proceedings of IEEE Conference on Systems, Man and Cybernetics, San Diego, USA, 1998. [14] A.J. Clemmens, T.F. Kacerek, B. Grawitz and W. Schuurmans, ”Test cases for canal control algorithm”, Journal of Irrigation and Drainage Engineering, vol. 124, no. 1, pp 23-30, 1998. [15] M. Foo, N. Bedjaoui and E. Weyer, ”Segmentation of a river using the Saint Venant equations”, Proceedings of IEEE Conference of Control Application, Yokohama, Japan, 2010. [16] S.K. Ooi, M.P.M. Krutzen and E. Weyer, ”On physical and data driven modelling of irrigation channels”, Control Engineering Practice, vol. 13, pp 461-471, 2005. [17] M.R. Schroeder, ”Synthesis of low-peak factor signals and binary sequences with low autocorrelation”, IEEE Transactions of Information Theory, vol. 16, pp 85-89, 1970. [18] L. Ljung, System Identification: Theory for the user, 2nd Edition, Prentice Hall, Englewood Cliffs, NJ; 1999. [19] T. Soderstrom and P. Stoica, System Identification, Prentice Hall, Englewood Cliffs, NJ; 1988.
4930
Nonlinearity tests of the Saint Venant equations Mathias Foo and Erik Weyer
Abstract— The Saint Venant equations are two nonlinear partial differential equations (PDE) which are used to describe the dynamics of one-dimensional flow in open water channels. Despite being nonlinear PDEs, the Saint Venant equations seem to exhibit linear behaviour in response to sinusoidal input signals. It is therefore of interest to determine ”how nonlinear” the Saint Venant equations are. In this paper, we investigate the nonlinearity in the Saint Venant equations using several commonly used nonlinearity tests suggested in the literature. Five different open water channels are considered, and the results from the nonlinearity tests show that the Saint Venant equations are nearly linear in an operating region from at least half the nominal flow to twice the nominal flows, and many of the channels display linear behaviour in a larger operating region. This finding is useful as it further justifies the use of linear control design methodologies for open water channels.
I. I NTRODUCTION Traditionally, the dynamics of open water channels1 are modelled by two nonlinear partial differential equations (PDE) which are known as the Saint Venant equations. The Saint Venant equations are derived from mass and momentum balance equations for one-dimensional flow and assuming no contribution from lateral flows, they are given by [1], ∂A ∂Q + =0 ∂t ∂x ∂Q ∂y ∂ Q2 + gA + + gA(Sf − S0 ) = 0 ∂t ∂x A ∂x
(1)
where Q is the flow, A is the wetted cross sectional area, T is the top width, g is the gravity constant, S0 is the bottom slope and Sf is the friction slope. The friction slope is a nonlinear function given by n2 Q2 P 4/3 A−10/3 , where P is the wetted perimeter and n is the Manning friction coefficient, which encapsulates the effect of flow resistance This work was supported by The Farms Rivers and Markets Project, an initiative of Uniwater and funded by the National Water Commission, the Victorian Water Trust, The Dookie Farms 2000 Trust (Tallis Trust) and the University of Melbourne and supported by the Departments of Sustainability and Environment and Primary Industry, the Goulburn Broken Catchment Management Authority and Goulburn-Murray Water. The first author also gratefully acknowledge the financial support from National ICT Australia (NICTA). NICTA is funded by the Australian Government as represented by the Department of Broadband, Communications and the Digital Economy and the Australian Research Council through the ICT Centre of Excellence program. M. Foo is with Asia Pacific Center for Theoretical Phyiscs (APCTP), Hogil Kim Memorial Building, #501, Pohang University of Science and Technology (POSTECH), San 31, Hyoja-dong, Nam-gu, Pohang, Gyeongbuk, 790-784, Korea [email protected] E. Weyer is with the Department of Electrical and Electronic Engineering, The University of Melbourne, 3010, Victoria, Australia
[email protected] 1 Open water channels are system where the water has a free surface to the air. Examples of open water channel include rivers and irrigation channels.
978-1-4673-5716-6/13/$31.00 ©2013 IEEE
and channel roughness. Typically, the cross sectional area is assumed to be a trapezoid,√and hence, A = (b + sy)y, T = b + 2sy and P = b + 2y 1 + s2 where b, s and y are the bottom width, side slope and water level respectively. As there is no closed form solution, the Saint Venant equations have to be solved numerically and the Preissmann scheme which is a finite difference method is commonly used to discretise the Saint Venant equations in time and space (see e.g. [1] - [3]). These equations are solved together with the boundary equations which gives the flows at the upstream and downstream end. For the downstream boundary condition, a sharp crested weir is used in this paper, whereby the flow and water level are related by Q = c(y − p)3/2 , where c is a constant, p is the height of the weir, and y > p is the water level (see Fig. 1). The initial condition is usually the steady state solution of Eqn. (1). The solutions of the Saint Venant equations are the flows and the water levels in space and in time. Despite being nonlinear PDEs, the
Fig. 1: Left: Sharp crested weir. Right: Schematic side view (not to scale). Saint Venant equations seem to exhibit linear behaviour in response to a sinusoidal input. When applying a sinusoidal input flow at the upstream end to the Saint Venant equations, we observe that the output, i.e. the downstream flow, is also a sinusoid of the same frequency but with a phase shift (see Fig. 2). No harmonic or subharmonic can be visually detected in the output. Thus, a natural suggestion is that the Saint Venant equations exhibit linear behaviour within an operating region. In this paper, we investigate ”how nonlinear” the Saint Venant equations are by applying some commonly used nonlinearity tests suggested in the literature. If the Saint Venant equations are almost linear within an operating region, then it is reasonable to conduct a frequency response analysis using the Saint Venant equations to obtain a Bode plot, which can be used for control design and in general it also justifies the use of linear control methodologies. This paper is organised as follows. In Section II a brief overview of nonlinearity tests is given. In Section III the
4925
Flow
Upstream flow Downstream flow
4
If the system is linear, Rx¯y (τ ) = 0, ∀τ and this test is claimed to be robust against noise [6].
m3/s
3.5 3
C. Higher order autocorrelation test
2.5 2 0.6
0.8
1
1.2 1.4 Time (minutes)
1.6
1.8
2 4
x 10
Fig. 2: Input and output response of the Saint Venant equations given a sinusoidal input at frequency, ω0 = 0.001 rad/min.
nonlinearity tests are applied to the Saint Venant equations. A frequency response analysis of open channel system is presented in Section IV, and conclusions are given in Section V. II. OVERVIEW
Ry¯2 y¯(τ ) = E{¯ y 2 (t − τ )¯ y (t)}
(5)
is computed with y¯(t) given by Eqn. (4). If the system is linear, Ry¯2 y¯(τ ) = 0, ∀τ . This test method has become very popular as it is simple and easy to apply [10]. Despite its popularity, this test method does raise some questions since while Ry2 y (τ ) = 0, ∀τ is necessary, it is not sufficient, and the test method fails for a simple cubic relationship (e.g. y(t) = x3 (t)) [11].
ON NONLINEARITY TESTS
An extensive list of works on nonlinearity tests can be found in [4]-[8]. Some of the commonly used tests include time domain tests, nonlinear cross-correlation tests, higher order autocorrelation tests and tests based on broadband periodic signal excitation. We will discuss some of the advantages and disadvantages of some of these tests. A. Time domain test The simplest and most popular test is the time domain test [8]. The input signal, uin,1 (t) is applied to the system and its corresponding output, uout,1 (t) is obtained. Then, the input signal, uin,1 (t) is scaled by a factor α1 (i.e. uin,1 (t) is replaced by α1 uin,1 (t) = uin,2 (t)). The signal uin,2 (t) is then applied to the system and its corresponding output, uout,2 (t) is obtained. If the system is linear, the output will be scaled by the same factor, i.e. uout,2 (t)/uout,1 (t) = α2 (t) = α1 . Although this test is simple, it is not widely used in practice as at least two experiments must be carried out. The choice of α1 may be limited by process constraints. If a small value of α1 is applied, the effect of noise is more significant, jeopardising the overall test. B. Nonlinear cross-correlation test A nonlinear cross-correlation test was introduced in [9]. In this test, the input signal is a nonzero mean white Gaussian noise or a pseudorandom signal. This signal is symmetrically distributed around the mean and has zero odd order moment and nonzero even order moment [6]. Introducing the term x(t) =
The higher order autocorrelation test in [7] is another nonlinearity test. Using a similar input signal as in the nonlinear cross-correlation test but this time with zero mean, the autocorrelation given by
[¯ u2 (t) − E{¯ u2 (t)}] 2 σ(¯ u (t))
(2)
where σ(.) denotes the standard deviation and u ¯(t) = [u(t)−E{u(t)}] , the cross-correlation σ(u(t)) Rx¯y (τ ) = E{x(t − τ )¯ y (t)}
(3)
is computed, where y¯(t) =
[y(t) − E{y(t)}] σ(y(t))
(4)
D. Broadband periodic excitation A more recent nonlinearity test suggested in [8] is based on using a broadband periodic excitation input signal. The idea originated from [12], where a signal that contains only certain harmonics is used to excite the system. One example of such a signal is the odd-odd multisines signal. The frequencies that are excited by the odd-odd multisines signal are (4k + 1)f0 , k = 0, 1, . . . , F , where F + 1 is the number of frequencies and f0 is the fundamental frequency. If the frequencies, (4k + 2)f0 , (4k + 3)f0 and (4k + 4)f0 are excited, then the system is nonlinear. III. N ONLINEARITY TESTS FOR OPEN CHANNEL SYSTEMS In this section, the nonlinearity tests presented in Section II are applied to five different open water channels. The first one is the river reach from Casey’s to Gowangardie Weir (Reach CG) in the Broken River, Australia. The second one is the Haughton Main Channel, Pool 9 (HMCP9) in Australia. The third and fourth one are the irrigation test channels 1 and 5, (ITC1 and ITC5), taken from [13]. The fifth one is the Test Canal 1, Pool 3 (TC1P3) taken from [14]. The Saint Venant equations for both Reach CG and HMCP9, have been calibrated and validated against measured data (see [15] and [16] respectively). ITC1 and ITC5, have been extensively used as benchmark test channels in several works (see e.g. [3]). TC1P3 is of interest since it is a short and steep channel. The physical parameters for all the five open channels are summarised in Table I. In all these tests, we consider the relationship between upstream flow and downstream flow. A. Time domain test The input signal is a pseudorandom binary signal (PRBS) as this input signal can excite a broad range of the spectrum enabling collection of all the required information from a single measurement [8]. The levels of the PRBS are 2.50m3/s and 3.50m3 /s, while the periods are chosen to be three times the time constant obtained from a step response, i.e. 960 minutes for Reach CG, 72 minutes for HMCP9, 216 minutes
4926
TABLE I: Open channel parameters. L is the length of the channel, b is the bottom width, s is the side slope, S0 is the bottom slope, n is the Manning friction coefficient. Parameters L b s S0 n
Reach CG 26700m 9-12m 2-3 0.0008-0.0020 0.060-0.085
HMCP9 853m 6m 2 0.0002 0.02
ITC1 3000m 7m 1.5 0.0001 0.02
ITC5 6000m 8m 1.5 0.0008 0.02
TC1P3 400m 1m 1.5 0.002 0.014
for ITC1, 225 minutes for ITC5 and 12 minutes for TC1P3. The values of α1 were 0.10, 0.20, 0.25, 0.33, 0.50, 2.00, 3.00, 4.00, 5.00 and 10.00.
B. Nonlinear cross-correlation test The input signal is a PRBS signal with levels 2.50m3/s and 3.50m3/s. This signal was applied to all five open water channels. The cross-correlation is computed using Eqn. (3) and they are shown in Fig. 3. In the plot, we included the threshold value of ±0.05 from Rx¯y (τ ) = 0, where the system is considered linear if the cross-correlation stays within this threshold value. We observe that the cross-correlation stays within the threshold for all τ . We repeated the test for several mean flows2 (i.e. mean flow from 2.00ms /s to 6.00ms/s) and also retaining the same mean flow of 3.00m3/s but with different levels for the PRBS signal (i.e. 2.90m3/s and 3.10m3/s, 2.45m3/s and 3.55m3/s to 2.00m3/s and 4.00m3/s). All tests yield similar results. The nonlinear cross-correlation
In our analysis we have allowed α2 (t) to deviate with 5% from α1 . This 5% percent tolerance bound is motivated from the definition of settling time of the step response; the settling time is defined as the time from when the output response stays within ±5% of its steady state value. Table II quantifies the percentage of the time for which α2 (t) exceeds the ±5% bound.
Rxy(τ)
Rxy(τ)
Reach CG 0.05 0 −0.05
ITC5 19.27 2.92 0.12 0.00 0.00 0.00 0.00 0.00 0.00 0.01
TC1P3 26.27 14.19 11.08 6.84 0.00 0.00 4.26 7.43 7.77 14.90
Rxy(τ)
ITC1 26.78 11.40 6.06 1.62 0.00 0.00 0.00 0.98 1.53 3.76
Rxy(τ)
HMCP9 24.04 9.63 5.19 1.45 0.00 0.00 0.00 0.16 1.37 2.57
Rxy(τ)
Reach CG 0.28 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
The results shown in Table II reveal several findings. For Reach CG, which is the longest one considered, the Saint Venant equations are linear for all considered α1 except for α1 = 0.10 but the percentage of time the bound was exceeded is very small. For ITC5, which is the second longest, we can see similar trend where the Saint Venant equations are linear for most of the considered α1 . As the open water channels get shorter and steeper, (i.e. from ITC1 to HMCP9 to TC1P3), we observe that the ’linear’ operating region gets smaller. One can note that α2 (t) exceeds the ±5% bound more often for low values of α1 , i.e. under low flow condition. One reason for the larger deviation at smaller α1 is that the resonances tend to show up more clearly at low flow. Based on the results from time domain test, we can say that the Saint Venant equations are approximately linear within an operating region ranging from at least half to two times the base flow, and for most of the channels the linear operating region is even larger.
100
200
300
400
500 600 τ HMCP9
700
800
900
1000
0
100
200
300
400
500 τ ICT1
600
700
800
900
1000
0
100
200
300
400
500 τ ICT5
600
700
800
900
1000
0
100
200
300
400
500 600 τ TC1P3
700
800
900
1000
0
100
200
300
400
700
800
900
1000
0.05 0 −0.05
TABLE II: Percentage of α2 (t) exceeding the ±5% bound. α1 0.10 0.20 0.25 0.33 0.50 2.00 3.00 4.00 5.00 10.00
0
0.05 0 −0.05
0.05 0 −0.05
0.05 0 −0.05 500 τ
600
Fig. 3: Nonlinear cross-correlation test plot with mean flow 3.00m3/s. Dashed line: Threshold. test show that all open water channels behave quite linearly for different mean flows as well as different flow levels. C. Higher order autocorrelation test The same input signal as in the nonlinear cross-correlation test is used but now with zero mean. The autocorrelation was computed using Eqn. (5). Note that a zero mean input signal means that we are dealing with negative flows, which is not physically feasible. Thus, for this test, a nonzero mean input signal is applied to the Saint Venant equations and then we remove the mean at the output of the Saint Venant equations 2 The difference between the two levels of PRBS signal remains the same but with different mean flow levels. E.g. mean flow of 2.00m3 /s has the two PRBS signal level of 1.50m3 /s and 2.50m3 /s.
4927
prior to computing the autocorrelation. The autocorrelation plots for all five open water channels are shown in Fig. 4. Using the same threshold value of ±0.05 as in Section
R
2
y y
(τ)
Reach CG 0.05 0 −0.05 0
100
200
300
400
500 τ
600
700
800
900
1000
600
700
800
900
1000
R
2
y y
(τ)
HMCP9 0.05 0 −0.05 0
100
200
300
400
500 τ
R
2
y y
(τ)
ICT1 0.05 0 −0.05 0
100
200
300
400
500 τ
600
700
800
900
1000
600
700
800
900
1000
where MNH is the magnitude of nonexcited harmonics and MEH is the magnitude of excited harmonics. Nne and Ne are the number of nonexcited and excited harmonics respectively. For illustration, the spectrum of the input and output signals for Reach CG and ITC5 are shown in Fig. 5.
0.05 0 −0.05 0
100
200
300
400
500 τ
Reach CG 0.015
Magnitude
2
y y
(τ)
ICT5
R
As the characteristics of the five open water channels are different, we considered frequencies up to where the magnitude of the output spectrum drops 20dB from the magnitude of the input spectrum. For quantitative measure, we calculate the N2E Ratio, which we define as the ratio of the normalised sum of the magnitude of the nonexcited harmonics to the normalised sum of the magnitude of the excited harmonics for each of the open channels up to the frequencies where the magnitude of the output spectrum drops 20dB from the magnitude of the input spectrum. The N2E Ratio is given by, PNne 1 MNH 1 Nne N2E Ratio = 1 PNe (7) 1 MEH Ne
R
2
y y
(τ)
TC1P3 0.05 0 −0.05
0.01
0.005
0
0
100
200
300
400
500 τ
600
700
800
900
1000
1
2
3
4 5 6 Frequency (rad/min)
7
8
9 −3
x 10
ITC5 0.015
III-B, we can see that, the autocorrelation stay within the threshold for all τ . We repeat the analyses using different mean flows and different flow levels for the PRBS signal (just as we did for the nonlinear cross-correlation test) and similar results were observed. The higher order autocorrelation tests show that all open water channels behave quite linearly for different mean flows and different flow levels. D. Broadband periodic excitation The odd-odd multisines signal is given by u(t) =
F X
[AB sin(ω0 (4k + 1)t + φk )] + BB
(6)
k=1
where AB is the amplitude, BB is a constant offset to ensure positive flow, F is the number of frequencies and φk = −k(k − 1)π/F are Schroeder phases [17]. They are constructed such that the peak to peak amplitude of the sum of the sinusoid signals is minimised (see e.g. [8] and [17]). The odd-odd multisines are applied as the input to the Saint Venant equations for all five open water channels with F in the range3 from 20 to 70000, AB = 0.025, BB = 3, and ω0 = 0.0001 rad/min. 3F
takes different values due to the different characteristics of the analysed channels.
Magnitude
Fig. 4: Higher order autocorrelation test plot with mean flow 3.00m3/s. Dashed line: Threshold.
0.01
0.005
0
0.01
0.02
0.03
0.04 0.05 0.06 Frequency (rad/min)
0.07
0.08
0.09
Fig. 5: Spectrum of odd-odd multisines signal for Reach CG and ITC5. -x: Input, -o: Output. From Fig. 5 we observe that for ITC5, the output spectrums have the same excited frequencies as the input spectrum. At the other frequencies, the output spectrum is not excited as much, (at least visually). For Reach CG, the output spectrum is excited at the same frequencies as present in the input spectrum. However, some other frequencies are also present. The magnitude of those frequencies are relatively small compared to the magnitude of the frequencies present in the input signal. The N2E Ratios for the five open water channels are given in Table III. From the N2E Ratios, we can see that all the TABLE III: N2E Ratio for the five open water channels. N2E Ratio
Reach CG 0.185
HMCP9 0.016
ITC1 0.028
ITC5 0.054
TC1P3 0.013
values are small. The largest ratio is from Reach CG but it is less than 0.20. The test is repeated with larger values of AB and similar N2E Ratios as in Table III are obtained.
4928
Reach CG
RESPONSE ANALYSIS
0
Magnitude (dB)
IV. F REQUENCY
All the nonlinearity tests presented above suggest that the Saint Venant equations are linear within an operating region. This indicates that a Bode plot can be constructed using a sine sweep4 .
−10
ω
3dB
−20 −30 −40 −3 10
−2
−1
10 rad/min
10
0
Phase (°)
A. Sine sweeping method Consider a sinusoidal input, ui (t) = AS sin(ω0 t)
−2000 −3 10
(8)
Mc (N ) φˆ = tan−1 Ms (N )
HMCP9
Magnitude (dB)
0
−5 ω
3dB
−10 −3 10
−2
10 rad/min
−1
10
0
−50
−100 −3 10
−2
10 rad/min
−1
10
(b) HMCP9
Fig. 6: Bode plot relating upstream flow to downstream flow for Reach CG and HMCP9
B. Results and discussions
(11)
where N is the number of data points. |G(jω0 )| and φ = ∠G(jω0 ) can be estimated as p Mc2 (N ) + Ms2 (N ) 2 ˆ (12) |G(jω 0 )| = AS and
−1
10
(10)
and N 1 X yo (t) sin(ω0 t) Ms (N ) = N t=1
−2
10 rad/min
(a) Reach CG
Phase (°)
yo (t) = BS sin(ω0 t + φ) + D(t) + transient + nonlinearities (9) where BS = AS |G(jω0 )|, φ = ∠G(jω0 ) = Im[G(jω0 )] , D(t) is the disturbances and G(jω) is the tan−1 Re [G(jω0 )] transfer function relating the input and the output. The effect of D(t) can be reduced by using a correlation method [18]. The transient effects can be reduced by not considering the initial part of the data. Defining N 1 X yo (t) cos(ω0 t) N t=1
−1000 −1500
where AS is the amplitude and ω0 is the frequency. If the system is linear and time invariant, the output response is a sinusoidal signal with the same frequency, scaled amplitude, and a phase shift. In the real system, the output response, yo (t) to the sinusoidal input also comprises transient effects, the effect of nonlinearities and disturbances.
Mc (N ) =
−500
(13)
(see [18] and [19] for details.). Having derived the expression for the magnitude and the phase of the transfer function, we can now obtain the frequency response analysis of all the five open water channels. The sinusoidal input flow, ui (t) = 0.75 sin(ω0 t) + 3.0 m3 /s is applied to the Saint Venant equations for all five open water channels and the experiment is repeated for several frequencies, ω0 in the frequency range of interest. By obtaining the magnitude and phase at each frequencies using Eqns. (12) and (13), the Bode plot can be constructed. 4 A sine sweep is of course not possible on real open channel system, but here we only consider the Saint Venant equations so the practical difficulties in applying a sine sweep test is not an issue.
The Bode plot for all five open water channels are shown in Figs. 6 and 7. The dashed line shown in the Bode plot for the amplitude and phase represents -3dB and -180◦ respectively. From the Bode plot of Reach CG (Fig. 6(a)), we can see that the frequency range relevant for control is similar to the frequency response of a pure time delay system. The 3dB bandwidth of the system is approximately 0.0032 rad/min but we note that the phase shift is already more than -180◦ at 0.0023 rad/min indicating a dominant time delay. At higher frequency range, we observed two small peaks at around 0.04 rad/min and 0.08 rad/min, which could be due to the effect of waves. For HMCP9, the Bode plot is shown in Fig. 6(b). In the frequency range relevant for control, the frequency response is similar to a first order system. The 3dB bandwidth of the system is approximately 0.035 rad/min, and at this frequency, we note that the phase is approximately -45◦ , which according to the theory of linear system is in agreement with a first order system. For ITC1, the Bode plot is shown in Fig. 7(a). Like HMCP9, in the frequency range relevant for control, the frequency response can also be approximated by a first order system. The 3dB bandwidth of the system is approximately 0.012 rad/min and at this frequency the phase is also approximately -45◦ , indicating a first order system. From the Bode plot of ITC5 (Fig. 7(b)), the 3dB band-
4929
ITC1
V. C ONCLUSIONS
Magnitude (dB)
0 −5
We applied some of the commonly used nonlinearity tests to the Saint Venant equations to determine how ”nonlinear” the Saint Venant equations are for a variety of channels. The results from the tests suggest that the Saint Venant equations are linear within a flow regime between at least half and two times the nominal flow, and for most channels this flow regime is even larger. Based on these findings, it is reasonable to construct Bode plots using sine sweeping methods, and the results further justify the use of linear control design methods for open water channels.
ω3dB
−10 −15 −20 −3 10
−2
−1
10 rad/min
10
0
Phase (°)
−50
−100
−150 −3 10
−2
−1
10 rad/min
10
(a) ITC1
R EFERENCES
ITC5
Magnitude (dB)
0 −5
ω3dB
−10 −15 −20 −25 −3 10
−2
−1
10 rad/min
10
0
Phase (°)
−100 −200 −300 −400 −500 −3 10
−2
−1
10 rad/min
10
(b) ITC5 TC1P3
Magnitude (dB)
0 −20 ω
3dB
−40 −60 −80 −3 10
−2
10
−1
10 rad/min
0
10
1
10
0
Phase (°)
−100 −200 −300 −400 −500 −3 10
−2
10
−1
10 rad/min
0
10
1
10
(c) TC1P3
Fig. 7: Bode plot relating upstream flow to downstream flow for ITC1, ITC5 and TC1P3.
width of the system is approximately 0.013 rad/min and at this frequency, the phase is approximately -95◦ . Moreover, at magnitude of -6dB, the phase is approximately -150◦, suggesting that the frequency range relevant for control is similar to the frequency response of a second order system with delay. For TC1P3, the Bode plot is shown in Fig. 7(c). Like the case of HMCP9 and ITC1, the frequency response is similar to a first order system where the 3dB bandwidth of the system is approximately 0.20 rad/min and the phase at this frequency is approximately -45◦ . We note that at higher frequency range, we observe peaks at around 5.00rad/min, which could be due to the effect of waves.
[1] M. Chaudhry, Open-channel Flow, 2nd Edition, Prentice Hall, Englewood Cliffs, NJ; 2007. [2] V.T. Chow, D.R. Maidment and L.W. Mays, Applied Hydrology, McGraw-Hill Series in Water Resources and Environmental Engineering; 1988. [3] X. Litrico and V. Fromion, Modeling and Control of Hydrosystems, Springer; 2009. [4] H.G. Natke, J-N. Juang, and W. Gawronski, ”A brief review on the identification of nonlinear mechanical system”, Proceedings of 6th International Modal Analysis Conference, Kissimee, USA, 1988. [5] R. Haber, ”Nonlinearity test for dynamic processes”, Proceedings of 7th IFAC/IFORS Symposium on Identification and Parameter Estimation, York, UK, 1985. [6] R. Haber and L. Keviczky, Nonlinear system identification - InputOutput Modeling Approach, Volume 2: Nonlinear System Structure Identification, Kluwer Academic Publishers; 1999. [7] S.A. Billings and W.S.F. Voon, ”Structure detection and model validity tests in the identification of nonlinear systems”, IEEE Proceedings, vol. 130, no. 4, pp 193-199, 1983. [8] R. Pintelon and J. Schoukens, System identification - A Frequency Domain Approach, 2nd Edition, IEEE Press; 2012. [9] R. Haber, ”Parametric identification of nonlinear dynamic systems based on correlation functions”, Proceedings of 5th IFAC Symposium on Identification and Parameter Estimation, Darmstadt, Germany, 1979. [10] P. Varlaki, G. Terdik and V.A. Lototsky, ”Tests for linearity and bilinearity of dynamic systems”, Proceedings of 7th IFAC/IFORS Symposium on Identification and Parameter Estimation, York, UK, 1985. [11] M.J. Korenberg and I.W. Hunter, ”The identification of nonlinear biological systems: Wiener kernel approaches”, Annuals of Biomedical Engineering, vol. 18, pp 629-654, 1990. [12] D.C. Evans, D. Rees and D.L. Jones, ”Identifying linear models of systems suffering nonlinear distortions”, Proceedings of IEE Control 94, Coventry, UK, 1994. [13] J-P. Baume, J. Sau and P-O. Malaterre, ”Modelling of irrigation channel dynamics for controller design”, Proceedings of IEEE Conference on Systems, Man and Cybernetics, San Diego, USA, 1998. [14] A.J. Clemmens, T.F. Kacerek, B. Grawitz and W. Schuurmans, ”Test cases for canal control algorithm”, Journal of Irrigation and Drainage Engineering, vol. 124, no. 1, pp 23-30, 1998. [15] M. Foo, N. Bedjaoui and E. Weyer, ”Segmentation of a river using the Saint Venant equations”, Proceedings of IEEE Conference of Control Application, Yokohama, Japan, 2010. [16] S.K. Ooi, M.P.M. Krutzen and E. Weyer, ”On physical and data driven modelling of irrigation channels”, Control Engineering Practice, vol. 13, pp 461-471, 2005. [17] M.R. Schroeder, ”Synthesis of low-peak factor signals and binary sequences with low autocorrelation”, IEEE Transactions of Information Theory, vol. 16, pp 85-89, 1970. [18] L. Ljung, System Identification: Theory for the user, 2nd Edition, Prentice Hall, Englewood Cliffs, NJ; 1999. [19] T. Soderstrom and P. Stoica, System Identification, Prentice Hall, Englewood Cliffs, NJ; 1988.
4930
