Xavier Litrico
email: xavier.litrico@cemagref.fr

Vincent Fromion
email: vincent.fromion@jouy.inra.fr

MODAL DECOMPOSITION OF LINEARIZED OPEN CHANNEL FLOW

Keywords: 2000 Mathematics Subject Classification. Primary: 35B37, 93C20, 93C80, 35L65; Secondary: 35C05, 35Q35 hyperbolic systems, shallow water equations, frequency domain approach

Open channel flow is traditionally modeled as an hyperbolic system of conservation laws, which is an infinite dimensional system with complex dynamics. We consider in this paper an open channel represented by the Saint-Venant equations linearized around a non uniform steady flow regime. We use a frequency domain approach to fully characterize the open channel flow dynamics. The use of the Laplace transform enables us to derive the distributed transfer matrix, linking the boundary inputs to the state of the system. The poles of the system are then computed analytically, and each transfer function is decomposed in a series of eigenfunctions, where the influence of space and time variables can be decoupled. As a result, we can express the time-domain response of the whole canal pool to boundary inputs in terms of discharges. This study is first done in the uniform case, and finally extended to the non uniform case. The solution is studied and illustrated on two different canal pools.

1. Introduction. The Saint-Venant equations are used by hydraulic engineers to represent open channel dynamics, either for design purposes, or simulation purposes, or even control purposes [START_REF] Cunge | Practical aspects of computational river hydraulics[END_REF]. These equations have been largely studied for practical implementation in softwares using specific numerical schemes. It is well-known that open channel dynamics are influenced by several factors, among which the Froude number, the channel slope, the geometry, the friction coefficient and the boundary conditions. However, their intimate behavior still remains unclear from the literature. As an example, there is still a debate among researchers about the dynamic properties of an open channel: How to characterize the dynamics of an open channel? How does the downstream boundary condition influence the pool dynamics? How to define the time delay of an open channel? Some of these questions were addressed by the ASCE Task Committee on Canal Automation Algorithms (see [START_REF] Strelkoff | Influence of canal geometry and dynamics on controllability[END_REF] and references therein), by using a hydraulic based approach with time-domain simulations. The paper [START_REF] Strelkoff | Influence of canal geometry and dynamics on controllability[END_REF] and related studies provided a dimensionless expression of Saint-Venant equations, and time domain simulations were performed in order to evaluate the effect of downstream boundary condition on the flow routing. However, these studies did not end with general results for the characterization of open channel flow.

We propose to tackle these questions by using a frequency domain approach, which gives a better insight into dynamical systems' properties than the time domain approach. We follow in this way the approach of automatic control engineers, but we will also illustrate the analysis in the time domain, in order to show the link between the two methods. Both methods are strongly related, but this is unfortunately too rarely pointed out (see e.g. [START_REF] Ermolin | Study of open-channel dynamics as controlled process[END_REF] for frequency domain analysis of Saint-Venant equations). Interesting results were obtained by Dooge et al. [START_REF] Dooge | The linear downstream response of a generalized uniform channel[END_REF] concerning the dynamic analysis of linearized Saint-Venant equations for a semiinfinite channel around uniform flow, also using a frequency domain approach. Tsai [START_REF] Tsai | Applicability of kinematic, noninertia, and quasi-steady dynamic wave models to unsteady flow routing[END_REF] provided a detailed analysis of Saint-Venant equations and simplified models for flow routing accounting for backwater effects. However, the developed criteria uses a local characterization of the flow, and does not provide clues for developing an integrated criteria based on these results. As we will show in the sequel, the dynamics of a reach are strongly linked to the poles of the system's transfer matrix, which are themselves linked to the reach characteristics (length, geometry, etc.). Therefore, our results show that an integrated criteria is necessary to characterize the reach dynamics.

Moussa and Bocquillon [START_REF] Moussa | Criteria for the choice of flood routing methods in natural channels[END_REF] and [START_REF] Moussa | Approximation zones of the Saint-Venant equations for flood routing with overbank flow[END_REF] have used a frequency analysis of the linearized Saint-Venant equations to provide a criteria for the choice of a simplified model. They considered the relative weight of different terms in the solutions of the equations in order to justify the choice of simplified models. This approach does not characterize the flow dynamics per se, but is based on a chosen threshold to separate the different models.

Ridolfi et al. [START_REF] Ridolfi | Green's function of the linearized de Saint-Venant equations[END_REF] computed the Green's function of the linearized Saint-Venant equations, which gives interesting results in terms of characterization of the flow. However, this does not provide a characterization of a given pool.

Thirriot and Benayada [START_REF] Thirriot | Solution analytique approchée des équations de l'écoulement transitoire à surface libre en canal avec perte de charge et pente du fond[END_REF] used the Laplace transform to derive an approximate solution in the uniform case, but they could not compute the poles in an explicit way because of the boundary condition they used. With a downstream boundary condition in terms of discharge, we have shown that the poles can be computed explicitly (see [START_REF] Litrico | Frequency modeling of open channel flow[END_REF]).

The approach developed in this paper is based on the poles of the system, which enable to fully characterize its dynamic behavior. The study is focused on the linearized equations, which capture the dynamics around stationary regimes. We provide a complete characterization in the uniform flow case, where closed-form solutions are available, and extend the results to the non uniform case using a numerical approach [START_REF] Litrico | Frequency modeling of open channel flow[END_REF]. We show that the non uniform flow case is mathematically speaking very close to the uniform flow case, even if closed-form solutions are no longer available. This opens the way to a proper classification of open channel flow characteristics, which could be based on the modal decomposition proposed in this paper.

The outline of the paper is as follows: after this introduction, Section 2 presents the Saint-Venant equations of open channel flow and the linearized equations around non uniform flow. Section 3 introduces the frequency domain approach for the linearized model around a uniform steady flow, with a complete characterization of the transfer matrix in terms of poles, delays and modal factors. These results are then extended to the linearized equations around a non uniform steady flow in Section 4, using a numerical approach.

Modeling of Open Channel Flow.

2.1. Saint-Venant Equations. The Saint-Venant equations are two coupled partial derivative equations involving the discharge Q(x, t) (m 3 /s) and the water depth Y (x, t) (m). The first one is the mass conservation equation:

∂A(x, t) ∂t + ∂Q(x, t) ∂x = 0 (1)
and the second one is the momentum conservation equation:

∂Q(x, t) ∂t + ∂ ∂x Q 2 (x, t) A(x, t) + gA(x, t) ∂Y (x, t) ∂x + S f (x, t) -S b (x) = 0 (2)
where S b (x) (m/m) is the bottom slope, A(x, t) (m 2 ) is the wetted area, and g (m/s 2 ) is the gravitational acceleration. The friction slope S f (x, t) (m/m) is modeled with the classical Manning formula [START_REF] Chow | Open-channel Hydraulics[END_REF]:

S f (x, t) = Q 2 (x, t)n 2 A 2 (x, t)R 4/3 (x, t) (3)
with n the Manning coefficient (sm -1/3 ) and R(x, t) the hydraulic radius (m), defined by R(x, t) = A(x, t)/P (x, t), where P (x, t) is the wetted perimeter (m) (see Fig. 1).

Figure 1.

Section of an open channel. Definition of the wetted perimeter P , wetted area A and top width T .

To complete the equations, we need to introduce initial and boundary conditions. The initial condition is given in terms of (Q(x, 0), Y (x, 0)), for all x ∈ [0, L], with L the length of the channel.

Steady Flow Solutions.

The steady flow solutions of Saint-Venant equations are obtained by replacing ∂ ∂t by 0 in Eqs. ( 1) and (2). Then, denoting the variables corresponding to the equilibrium regime with a subscript zero (Q 0 (x), Y 0 (x), etc.), the Saint-Venant equations become:

dQ 0 (x) dx = 0 (4a) dY 0 (x) dx = S b -S f 0 (x) 1 -F 0 (x) 2 . ( 4b
)
F 0 is the Froude number F 0 (x) = V0(x) C0(x) with C 0 (x) = g A0(x) T0(x) the gravitational wave celerity (m/s), V 0 (x) = Q0(x) A0(x) the flow velocity (m/s). These two equations define an equilibrium regime given by Q 0 (x) = Q 0 and Y 0 (x) solution of the ordinary differential equation (4b), for a boundary condition in terms of downstream elevation.

When the right-hand side of equation ( 4b) is equal to zero, the water depth is constant along the channel. In this case, given Q 0 (x) = Q 0 , the equilibrium solution Y 0 (x) = Y n (also called normal depth) can be deduced by solving the following algebraic equation:

S f (Q 0 , Y n ) = S b . (5)
This specific solution corresponds to the uniform flow regime. Equation ( 5) is usually solved numerically with a fixed-point or Newton-Raphson method.

In some cases, the uniform depth can be computed analytically. For large rectangular channels, the hydraulic radius R can be approximated by the water depth Y , and the Manning equation (3) reduces to:

S f = Q 2 n 2 T 2 Y 10/3 . (6)
Combining equations ( 5) and [START_REF] Dooge | The linear downstream response of a generalized uniform channel[END_REF] give the uniform depth Y n corresponding to a discharge Q 0 in the large rectangular case:

Y n = Q 2 0 n 2 T 2 S b 3/10 . (7)
The uniform flow regime will be specifically developed as an example, since this regime leads to closed-form analytical solutions for the transfer matrix. However, we will also show that this specific flow regime is not qualitatively different from realistic non uniform flow regimes.

The article will be illustrated on two trapezoidal prismatic channels, with different characteristics (see Fig. 2).

Figure 2. Section of a trapezoidal canal

The two example canals are presented in the following.

Example 1. The channel characteristics are given in table 1, where L is the channel length (m), m the bank slope, B the bed width (m), S b the bed slope (m/m), n the Manning coefficient (m -1/3 s), Y n the normal depth (m) corresponding to the discharge Q 0 (m 3 s -1 ). Canal 1 is a short oscillating canal, and canal 2 is a long delayed canal.

Keeping a constant discharge, we compare the backwater curves obtained for different downstream boundary conditions Y 0 (L) = Y n × [0.8, 1, 1.2] (see figure 3). We have two types of backwater curves: if Y 0 (L) > Y n , then the flow is decelerating along x, this is a so-called 'M1' curve (see e.g. [START_REF] Sturm | Open-Channel Hydraulics[END_REF]), and if Y 0 (L) < Y n , then the Y n Q 0 Canal 1 3000 1.5 7 0.0001 0.02 2.12 14 Canal 2 6000 1.5 8 0.0008 0.02 2.92 80 flow is accelerating along x, this is a 'M2' curve. The M1 curves are the most widely observed in practice, since they may occur upstream of any kind of obstacle in the flow (hydraulic structure, gate, weir, bridge, etc.). The M2 type are observed when there is a sudden drop in the canal or a steepest slope downstream. 3. Linearized Equations. The Saint-Venant equations are linearized around an equilibrium steady state (Q 0 (x), Y 0 (x)). The linearized equations are obtained by putting Q(x, t) = Q 0 (x) + q(x, t), and A(x, t) = A 0 (x) + a(x, t) into equations ( 1) and (2) and expanding in series.

After collecting terms in q, a, ∂q ∂x and ∂a ∂x , the linearized equations are obtained as follows:

∂a ∂t + ∂q ∂x = 0 (8a) ∂q ∂t + (α -β) ∂q ∂x + αβ ∂a ∂x + δq -γa = 0 (8b) with α = C 0 + V 0 β = C 0 -V 0 γ = C 2 0 T 0 dT 0 dx + g (1 + κ)S b -(1 + κ -(κ -2)F 2 0 ) dY 0 dx δ = 2g V 0 S b - dY 0 dx κ = 7 3 - 4A 0 3T 0 P 0 ∂P 0 ∂Y
The solution in terms of variations of water level y(x, t) can readily be deduced by dividing by T 0 (x) the solution obtained in terms of variations of wetted area a(x, t) = T 0 (x)y(x, t).

To facilitate the mathematical analysis, we rewrite the linearized Saint-Venant equations as follows:

∂ξ ∂t + A(x) ∂ξ ∂x + B(x)ξ = 0 (9)
where ξ(x, t) = (a(x, t), q(x, t)

) T : [0, L] × [0, +∞) → Ω ∈ R 2 is the state of the system and A(x) = 0 1 α(x)β(x) α(x) -β(x) , B(x) = 0 0 -γ(x) δ(x)
.

The initial condition and the boundary conditions are given by: q(x, 0) = q 0 (x), y(x, 0) = y 0 (x) (10a)

q(0, t) = u 1 (t), q(L, t) = u 2 (t). ( 10b
)
We also assume that the water level deviations are measured at each boundary: y(0, t) and y(L, t).

Uniform Flow Case.

In the uniform flow case, the parameters α, β, γ and δ are constant and given by:

α = C 0 + V 0 β = C 0 -V 0 γ = g(1 + κ)S b δ = 2gS b V 0 , with κ = 7 3 -4A0 3T0P0 ∂P0 ∂Y .
3.1. Derivation of the Transfer Matrix. We will now derive the transfer matrix for the linearized Saint-Venant equations, taking q(0, t) and q(L, t) as boundary inputs, and the water depth deviations y(x, t) and discharge deviations q(x, t) as outputs.

We apply Laplace transform to the linear partial differential equations (9), using the classical relation df dt = s f(s)f (0) with s the Laplace variable, which yields:

∂ ξ(x, s) ∂x = A(s) ξ(x, s) + Bξ(x, 0) (11)
with A(s) = -A -1 (sI + B) and B = A -1 , i.e.:

A(s) = 1 αβ (α -β)s + γ -s -δ -αβs 0 B = 1 αβ (β -α) 1 αβ 0 .
Let us diagonalize matrix A(s):

A(s) = X (s) -1 D(s)X (s) (12)
Author

1 s λ1(s) 1 (13a) D(s) = λ 1 (s) 0 0 λ 2 (s) (13b) X (s) -1 = λ 1 (s)λ 2 (s) s(λ 1 (s) -λ 2 (s)) 1 -1 -s λ1(s) s λ2(s) , (13c)
and where λ 1 (s) and λ 2 (s) are the eigenvalues of A(s), solutions of the equation:

det(λI -A(s)) = 0.
This equation is a second order polynomial equation in λ:

αβλ 2 -[(α -β)s + γ]λ -s(s + δ) = 0, (14)
which has in general two solutions:

λ 1 (s) = (α -β)s + γ -d(s) 2αβ (15a) λ 2 (s) = (α -β)s + γ + d(s) 2αβ (15b) with d(s) = (α + β) 2 s 2 + 2[(α -β)γ + 2αβδ]s + γ 2 .
3.1.

with Φ(x, s) = X (s) -1 e D(s)x X (s) [START_REF] Strelkoff | Influence of canal geometry and dynamics on controllability[END_REF] and ξ0 (x, s) = Φ(x, s)

x 0 Φ(v, s) -1 Bξ(v, 0)dv.
The elements of matrix Φ(x, s) are given by:

φ 11 (x, s) = λ 1 (s)e λ1(s)x -λ 2 (s)e λ2(s)x λ 1 (s) -λ 2 (s) (18a) φ 12 (x, s) = λ 1 (s)λ 2 (s)(e λ1(s)x -e λ2(s)x ) s(λ 1 (s) -λ 2 (s)) (18b) φ 21 (x, s) = s(e λ2(s)x -e λ1(s)x ) λ 1 (s) -λ 2 (s) (18c) φ 22 (x, s) = λ 1 (s)e λ2(s)x -λ 2 (s)e λ1(s)x λ 1 (s) -λ 2 (s) . (18d)
Φ(x, s) is the state-transition matrix for the differential equation (11). For simplicity, we assume zero initial conditions for the following developments, therefore ξ0 (x, s) = 0.

3.1.2. Boundary Conditions. Specifying the boundary conditions leads to express the state as a function of the boundary discharges. To this end, we use Eq. ( 16) for x = 0, and get the equality:

T 0 ŷ(0, s) = - φ 22 (L, s) φ 21 (L, s) q(0, s) + 1 φ 21 (L, s) q(L, s). (19)
Once again, this is valid provided φ 21 (L, s) is not equal to zero. The values of s such that φ 21 (L, s) = 0 correspond to the poles of the Saint-Venant transfer matrix.

3.1.3. Poles. The poles are obtained as the solutions of:

φ 21 (L, s) = 0, ( 20
)
where φ 21 (x, s) is given by Eq. (18c).

If s 0 is such that λ 1 (s 0 ) = λ 2 (s 0 ), we have: 20) is equivalent to

φ 21 (L, s 0 ) = -s 0 Le λ1(s0)L therefore s 0 is not a pole. When λ 1 (s) = λ 2 (s), equation (
s(e λ2(s)L -e λ1(s)L ) = 0.
There is a pole in zero p 0 = 0, which means that the canal pool acts as an integrator and the other poles verify the equation:

λ 1 (s) -λ 2 (s) = 2jkπ L ,
with k ∈ N * . The case k = 0 leads to λ 1 (s) = λ 2 (s), therefore s = s 0 , which is not a pole. Finally, the poles are solutions of the second order equation:

d(s) = - 4α 2 β 2 k 2 π 2 L 2 , with k ∈ N * .
The poles (p ±k ) k∈N * are then given by:

p ±k = - (α -β)γ + 2αβδ (α + β) 2 ± 2α 2 β 2 (α + β) 2 Δ(k) (21) with Δ(k) = (α+β) 2 α 2 β 2 (αδ-γ)(βδ+γ) αβ(α+β) 2 -k 2 π 2 L 2 .
Let k m ∈ N * be the greatest integer such that Δ(k m ) ≥ 0. Then the poles obtained for 0 < k ≤ k m are negative real, and those obtained for k > k m are complex conjugate, with a constant real part equal to

-(α-β)γ+2αβδ (α+β) 2
, which means that they are located on a vertical line in the left half plane. Canal pools with a dominant oscillating behavior correspond to Δ(1) < 0.

Example 2 (Poles of the example canals). The poles of the canals 1 and 2 are depicted in Fig. 4. Canal 1 has an oscillating behavior, all its poles being complex conjugate. Canal 2 has two negative real poles for k m = 1, and the other ones are complex conjugate. The complex poles of canal 2 have a larger negative real part than the ones of canal 1. It should therefore exhibit a damped behavior compared to the canal 1. 19) and the state-transition matrix [START_REF] Sturm | Open-Channel Hydraulics[END_REF], we obtain the input-output transfer matrix:

ŷ(0, s) ŷ(L, s) = p 11 (s) p 12 (s) p 21 (s) p 22 (s) q(0, s) q(L, s) (22)
with p 11 (s) = -φ22(L,s) T0φ21(L,s) , p 12 (s) =

T0φ21(L,s) , p 21 (s) = φ12(L,s) T0 -φ11(L,s)φ22(L,s) 1
T0φ21(L,s)

and p 22 (s) = φ11(L,s) T0φ21(L,s) . This leads to the final expressions:

p 11 (s) = λ 2 (s)e λ1(s)L -λ 1 (s)e λ2(s)L
T 0 s e λ2(s)Le λ1(s)L (23a)

p 12 (s) = λ 1 (s) -λ 2 (s) T 0 s e λ2(s)L -e λ1(s)L (23b) p 21 (s) = (λ 2 (s) -λ 1 (s))e (λ1(s)+λ2(s))L T 0 s e λ2(s)L -e λ1(s)L (23c) p 22 (s) = λ 1 (s)e λ1(s)L -λ 2 (s)e λ2(s)L T 0 s e λ2(s)L -e λ1(s)L . ( 23d
)
This result is similar to the expressions obtained by several authors [START_REF] Constantinides | Stability of a feedback controlled distributed system by model representation[END_REF][START_REF] Corriga | A mathematical model for open-channel networks[END_REF][START_REF] Schuurmans | Open-channel flow model approximation for controller design[END_REF][START_REF] Baume | Study of irrigation canal dynamics for control purposes[END_REF]. We depict in the following example the transfer functions p 21 (s) and p 22 (s) which are useful in a context of automatic control of open channels.

Example 3 (Bode plots of the input output transfer matrix). Bode plots of transfer functions p 21 (s) and p 22 (s) for example canals 1 and 2 are depicted in Figs. 5 and6. We observe that the low frequency behavior of the system is dominated by the integrator, that there is a delay in transfer function p 21 (s). In high frequencies, the oscillating modes of canal 1 are clearly visible, while canal 2 is damped and tends towards a constant gain. After elementary manipulations, we obtain the distributed transfer matrix, relating the water depth ŷ(x, s) and the discharge q(x, s) at any point x in the canal pool to the upstream and downstream discharges:

ŷ(x, s) q(x, s) = G(x, s) q(0, s) q(L, s) (24)
Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325 with g 11 (x, s) = λ 2 (s)e λ2(s)x+λ1(s)Lλ 1 (s)e λ1(s)x+λ2(s)L T 0 s(e λ2(s)Le λ1(s)L ) (25a)

g 12 (x, s) = λ 1 (s)e λ1(s)x -λ 2 (s)e λ2(s)x T 0 s(e λ2(s)L -e λ1(s)L ) (25b) g 21 (x, s) = e λ1(s)x+λ2(s)L -e λ2(s)x+λ1(s)L e λ2(s)L -e λ1(s)L (25c) g 22 (x, s) = e λ2(s)x -e λ1(s)x e λ2(s)L -e λ1(s)L . ( 25d
)
The input-output transfer matrix P (s) can be obtained from the distributed transfer matrix G(x, s) as follows: p 11 (s) = g 11 (0, s), p 12 (s) = g 11 (L, s), p 21 (s) = g 12 (0, s) and p 22 (s) = g 12 (L, s).

Example 4 (Spatial Bode plot of the distributed transfer matrix). The spatial Bode plots of the distributed transfer function g 11 (x, s) are depicted in Figs. 7 and8 for example canals 1 and 2. The amplitude of the transfer functions is depicted in dB vs a log scale for the frequency.

In the spatial Bode plot of canal 1, the integrator is clearly visible in low frequencies. For higher frequencies, the oscillating modes can be seen with the corresponding nodes and anti-nodes of oscillations along the channel. The canal 2 is much more damped than canal 1, but there are oscillating modes at the downstream end of the canal. This means that the water level response to the upstream discharge may exhibit local oscillations close to the downstream end of the canal. The surprising thing here is that these oscillations do not appear in the input-output transfer functions (see Fig. 6). We finally get the following factorization of all elements of the transfer matrix:

Properties of the

p 11 (s) = p11 (s) p 12 (s) = p12 (s)e -τ2s p 21 (s) = p21 (s)e -τ1s p 22 (s) = p22 (s) with p12 (s) = (λ 1 (s) -λ 2 (s)) T 0 s 1 -e (λ1(s)-λ2(s))L e τ2s-λ2(s)L p21 (s) = (λ 2 (s) -λ 1 (s)) T 0 s 1 -e (λ1(s)-λ2(s))L e τ1s+λ1(s)L .
Let us note that when |s| → ∞, we have

λ 1 (s)L = -τ 1 s + O(1) λ 2 (s)L = τ 2 s + O(1).
This explains why the terms e τ1s+λ1(s)L and e τ2s-λ2(s)L are delay free, since a delay can be pointed out in high frequencies by its effect on the phase of the transfer function.

This result shows that the Saint-Venant transfer matrix includes delays only linked to the waves propagations, which correspond to the delays obtained by the characteristics.

Distributed Transfer Matrix. The delays are obtained directly as e -x α s for the distributed transfer functions g 11 (x, s) and g 21 (x, s) and e -L-x β s for the distributed transfer functions g 12 (x, s) and g 22 (x, s). The transfer functions can therefore be factorized as

g 11 (x, s) = g11 (x, s)e -x α s g 12 (x, s) = g12 (x, s)e -L-x β s g 21 (x, s) = g21 (x, s)e -x α s g 22 (x, s) = g22 (x, s)e -L-x β s
where the delay-free parts gij (x, s) are given by:

g11 (x, s) = λ 2 (s)e (λ1(s)-λ2(s))(L-x) -λ 1 (s) T 0 s(1 -e (λ1(s)-λ2(s))L ) e s α +λ1(s) x g12 (x, s) = λ 1 (s)e (λ1(s)-λ2(s))x -λ 2 (s) T 0 s(1 -e (λ1(s)-λ2(s))L ) e s β -λ2(s) (L-x) g21 (x, s) = 1 -e (λ1(s)-λ2(s))(L-x) 1 -e (λ1(s)-λ2(s))L e s α +λ1(s) x g22 (x, s) = 1 -e (λ1(s)-λ2(s))x 1 -e (λ1(s)-λ2(s))L e s β -λ2(s) (L-x) .
3.2.2. Asymptotic Estimate of the Poles. The following proposition proved in [11] provides an asymptotic estimate of the poles for high frequencies.

Proposition 1. For k ∈ N, k k h = L π(α+β) (αδ-γ)(βδ+γ) αβ
the poles of the Saint-Venant transfer matrix (22) can be approximated by:

p ±k ≈ - (r 1 + r 2 )L τ 1 + τ 2 ± 2jkπ τ 1 + τ 2 (26)
with r 1 = αδ-γ α(α+β) , r 2 = βδ+γ β(α+β) and τ 1 = L α the delay for downstream propagation, τ 2 = L β , the delay for upstream propagation. One can show that this high frequency approximation leads to a relative error smaller than 10% when k > 5.7 × k h and smaller than 5% when k > 10.7 × k h .

This proposition shows that for high frequencies the poles of Saint-Venant transfer matrix are close to the ones of the following damped wave equation:

∂ ∂x + r 1 + 1 α ∂ ∂t ∂ ∂x -r 2 - 1 β ∂ ∂t q = 0
with boundary conditions q(0, t) and q(L, t). Using Laplace transform, this equation reduces to an ODE in x, with eigenvalues equal to -r 1 -s α and r 2 + s β . The obtained transfer function has for denominator 1e -(r1+r2)L-(τ1+τ2)s , whose roots coincide with the poles approximation (26). This shows that the oscillating modes correspond to the interaction of two gravity waves, one traveling downstream at speed α = V 0 + C 0 with attenuation factor r 1 , and one traveling upstream at speed β = C 0 -V 0 with attenuation factor r 2 .

Modal Decompositions.

Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009. 4.325 Rational Series Expansion. For simplicity, we assume in the following that the poles have single multiplicity, i.e. that Δ(k) = 0, but the solution can easily be extended to the case where Δ(k) = 0. Then, the Cauchy residues theorem implies that each transfer function g ij (x, s) can be decomposed as an infinite sum:

g ij (x, s) = b ij (x) + a (0) ij (x) s + ∞ k=-∞,k =0 a (k) ij (x)s p k (s -p k ) ( 27
)
with

a (k) ij (x) = lim s→p k (s -p k )g ij (x, s) (28)
and

b ij (x) = ∂ ∂s [sg ij (x, s)]| s=0 . (29)
The coefficient a

(k) ij (x)
is the residue of transfer function g ij (x, s) at the pole p k .

Proof. The proof can be easily adapted from the one in Appendix A.

The residues can be computed analytically with Eq. ( 28), leading to, for k = 0:

a (k) 11 (x) = -sgn(k) 2k 2 π 2 e ((α-β)p k +γ) 2αβ x T 0 L 3 p k Δ(k) × cos kπx L + ((α -β)p k + γ)L 2αβkπ sin kπx L (31a) a (k) 12 (x) = sgn(k)(-1) k 2k 2 π 2 e ((α-β)p k +γ) 2αβ (x-L) T 0 L 3 p k Δ(k) × cos kπx L + ((α -β)p k + γ)L 2αβkπ sin kπx L (31b) a (k) 21 (x) = -sgn(k) 2kπe ((α-β)p k +γ) 2αβ x L 2 p k Δ(k) sin kπx L (31c) a (k) 22 (x) = sgn(k)(-1) k 2kπe ((α-β)p k +γ) 2αβ (x-L) L 2 p k Δ(k) sin kπx L , (31d)
and the b ij (x) are given by:

b 11 (x) = e ψ x L T 0 (e ψ -1) 2 α -β αβ e ψ (1 + ψ( x L -1)) -1 -ψ x L + δ γ e ψ(1-x L ) (e ψ -1) + e ψ (1 + ψ( x L -2)) -1 -ψ x L (32a) b 12 (x) = 1 T 0 (e ψ -1) 2 α -β αβ e ψ x L ψe ψ + (1 -e ψ )(1 + ψ x L ) + δ γ (1 + e ψ )ψe ψ x L + (1 -e ψ ) 1 + e ψ x L (1 + ψ x L ) (32b) b 21 (x) = e ψ -e ψ( x L -1) e ψ -1 (32c) b 22 (x) = e ψ x L -1 e ψ -1 , (32d)
with ψ = γL αβ . The modal factors a (k) ij (x) are combination of sine and cosine functions of kπx L , modified by a multiplicative exponential term which depends on the difference α-β and the coefficient γ. Coefficients γ and δ are directly linked to the slope of the open channel. The larger the slope, the larger their influence on the modal factors. Rational Decomposition of the Delay-free Part. We have shown that the transfer functions g ij (x, s) could be factorized as a delay times a delay-free part. Therefore, one can find an approximation of g ij (x, s) by applying the Cauchy residues theorem on transfer functions gij (x, s), and then adding the delay. Finally, each transfer function gij (x, s) can be decomposed as:

gij (x, s) = bij (x) + ã(0) ij (x) s + ∞ k=-∞,k =0 ã(k) ij (x)s p k (s -p k ) (33) with ã(k) ij (x) = lim s→p k (s -p k )g ij (x, s) (34)
and bij (x) = ∂ ∂s [sg ij (x, s)]| s=0 . (35)
The coefficient ã(k) ij (x) is the residue of transfer function gij (x, s) at the pole p k .

Proof. The proof can be easily adapted from the one in Appendix A.

The residues can be computed analytically with Eq. (34), leading to, for k = 0

ã(0) ij (x) = a (0) ij (x) (36)
and, for k = 0: and the bij (x) are given by: b11

ã(k) 11 (x) = e x α p k a (k) 11 (x) (37a) ã(k) 12 (x) = e L-x β p k a (k) 12 (x) (37b) ã(k) 21 (x) = e x α p k a (k) 21 (x) (37c) ã(k) 22 (x) = e L-x β p k a (k) 22 (x), (37d)
(x) = b 11 (x) + x α a (0) 11 (x) (38a) b12 (x) = b 12 (x) + (L -x) β a (0) 12 (x) (38b) b21 (x) = b 21 (x) (38c) b22 (x) = b 22 (x). ( 38d
)
Example 5 (Modal factors of the example canals). The functions a (k) 11 (x) of the distributed modal decomposition of g 11 (x, s) of the example canals 1 and 2 are depicted in Fig. 9 and10.

For canal 1, which is flat and has oscillating modes, the modes of oscillation are clearly visible. The integrator in this case is slightly modified by the slope and the friction. We clearly see the oscillating modes and the associated nodes and anti-nodes of oscillation. A node corresponds to a point where the traveling waves interfere negatively with each other. An anti-node corresponds to a positive interference: the perturbations are in phase and act additively. For canal 2, which is steep, we see that the dominant term is the integrator, which is greatly modified by the slope. The poles p ±1 are negative real, therefore not oscillating, and the ones obtained for |k| > 1 are oscillating. The corresponding modal factors are strongly modified by the slope, compared to the ones of the example canal 1.

Figure 11 depicts the value of coefficients b 11 (x) and b 12 (x) of the distributed modal decomposition of the transfer functions g 11 (x, s) and g 12 (x, s) for canals 1 and 2.

The term b ij (x) corresponds to a direct term in transfer function g ij (x, s). We see that this term varies a lot with x for both canals. It even changes sign, which is a surprising result: this means that the water level change resulting from a step 4. Non Uniform Flow Case. We now consider the general case of non uniform flow, where analytical results are no longer possible. We first show that the computation of the transfer matrix is equivalent to solving a specific ordinary differential equation parameterized by the complex Laplace variable s. The application of classical numerical integration methods, such as Runge-Kutta, faces a main difficulty. Due to structural properties of the considered differential equation, the classical compromise between precision of the approximate solution and integration step size necessarily leads to a large computational time for each value of s. The frequency domain of interest goes from the low frequency behavior of the canal (typically ω r /100 with ω r the resonant frequency ω r = 2π/(τ 1 + τ 2 ) with τ 1 = L/(V 0 + C 0 ) and τ 2 = L/(C 0 -V 0 ) respectively the downstream and upstream propagation timedelays) to the high frequency behavior (typically 10 times the resonant frequency). Since it is important to have a good representation in between these extreme frequencies (and especially to correctly reproduce the resonant modes), it is necessary to compute a large number of frequency points (typically 500 points) [START_REF] Litrico | Frequency modeling of open channel flow[END_REF]. To this end, it is therefore essential to have an efficient numerical method to compute the frequency response.

Computation of the Transfer Matrix.

In the non uniform flow case, it is no longer possible to diagonalize the system. Indeed, applying Laplace transform and after elementary manipulations, we obtain the following differential equation:

∂ ξ(x, s) ∂x = A(x, s) ξ(x, s) + B(x)ξ(x, 0) (39)
with ξ(x, t) = (T 0 (x)y(x, t), q(x, t)) T and where A(x, s) and B(x) depend on x:

A(x, s) = 1 α(x)β(x) (α(x) -β(x))s + γ(x) -s -δ(x) -α(x)β(x)s 0 B(x) = 1 α(x)β(x) (β(x) -α(x)) 1 α(x)β(x) 0 .
To simplify the exposition, we assume in the following that ξ(x, 0) = 0.

One may show that the change of variable that diagonalizes A(x, s) introduces a supplementary term due to the derivative of matrix X (x, s) with respect to x. Indeed, let us introduce ζ(x, s) = X (x, s) ξ(x, s), where X (x, s) is such that A(x, s) = X (x, s) -1 D(x, s)X (x, s). Plugging this variable into Eq. (39) leads to:

∂ ζ(x, s) ∂x = D(x, s) + ∂X (x, s) ∂x X (x, s) -1 ζ(x, s). (40)
The matrix D(x, s)

+ ∂X (x,s) ∂x X (x, s) -1 is not diagonal, because of the term ∂X (x,s) ∂x X (x, s) -1 .
Therefore we can no longer provide an analytical solution to this equation. However, since the differential equation ( 39) is linear, we know that its general solution always exists, is unique and is given by [9, p. 598

]: ξ(x, s) = Φ(x, s) ξ(0, s) = φ 11 (x, s) φ 12 (x, s) φ 21 (x, s) φ 22 (x, s) ξ(0, s) (41)
where Φ(x, s) is the state-transition matrix associated to the differential equation (39). As a matter of fact, one may think that, as in the scalar case, the state-transition matrix could be written as Φ(x, s) = e x 0 A(v,s)dv . However, this is not true in

∂ ∂x e x 0 A(v,s)dv = A(x, s) + A(x, s) 2 x 0 A(v, s)dv + x 0 A(v, s)dv A(x, s) 2 + . . . = A(x, s)e x 0 A(v,s)dv !
The only case where this holds is when A(x, s) and

x 0 A(v, s)dv commute. This is true e.g. when A(x, s) is constant with respect to x, i.e. in the uniform flow case. In the non uniform flow case, this is generally not true, this is why we need to find a numerical way to compute the state-transition matrix.

Once the state-transition matrix is obtained, the transfer matrix corresponding to the original differential equation ( 39) is then given by:

ŷ(0, s) ŷ(L, s) = p 11 (s) p 12 (s) p 21 (s) p 22 (s) q(0, s) q(L, s) (42)
with

p 11 (s) = -φ22(L,s) T0(0)φ21(L,s) , p 12 (s) = 1 T0(0)φ21(L,s) , p 21 (s) = φ12(L,s) T0(L) -φ11(L,s)φ22(L,s) T0(L)φ21(L,s)
and p 22 (s) = φ11(L,s) T0(L)φ21(L,s) , provided φ 21 (L, s) is not equal to zero. As in the uniform flow case, the values of s such that φ 21 (L, s) = 0 correspond to the poles of the Saint-Venant transfer matrix.

The problem is therefore to solve the differential equation (39) and compute the state-transition matrix Φ(x, s). This problem is not easy to solve, due to its oscillatory and unstable nature.

We proposed an efficient numerical solution to this problem in [START_REF] Litrico | Frequency modeling of open channel flow[END_REF], leading to a general method to compute the Saint-Venant transfer matrix for any kind of non uniform flow. 4.1.1. Poles Computation. As we have seen, the poles correspond to the values of s such that φ 21 (L, s) = 0. In the non uniform case, the poles are obtained numerically, by finding the zeros of equation ( 20). Since we can compute numerically the transition matrix, we know how to compute φ 21 (L, s) for any s ∈ C, and we can look for values of s such that φ 21 (L, s) = 0.

Practically, we seek the value of s that minimizes the modulus of φ 21 (L, s). This search is done starting from the values obtained analytically for the uniform regime.

We used this numerical method to find the values of s that minimize the modulus of φ 21 (L, s). This has been done for the first 7 pairs of poles for the canals in example 1.

The results are depicted in Figs. 12 and 13 for canals 1 and 2, respectively. In both figures, we see that an accelerating flow tends to dampen the oscillating poles, while a decelerating flow tends to reduce their damping. However, the first negative real pole of canal 2 reacts differently: its damping increases for decelerating flow.

This confirms our results from the asymptotic study in high frequencies, but this does not necessarily apply for low frequencies.

(s) = 1 T 0 (0) φ 12 (0, s) -φ 11 (0, s) φ 22 (L, s) φ 21 (L, s) (43a)
p 12 (s) = 1 T 0 (0) φ 11 (0, s) φ 21 (L, s) (43b)
p 21 (s) = 1 T 0 (L) φ 12 (L, s) -φ 11 (L, s) φ 22 (L, s) φ 21 (L, s) (43c)
p 22 (s) = 1 T 0 (L) φ 11 (L, s) φ 21 (L, s) . ( 43d
)
Example 6 (Bode plots for non uniform flow). We compute the Bode plots for the three situations: uniform flow, decelerating flow (M1 backwater curve) and accelerating flow (M2 backwater curve). The results are depicted in Figs. 14 and15.

The integrator gain increases for accelerating flow, while it decreases for decelerating flow. This also is consistent with the physical intuition: in a first approximation, the integrator gain is inversely proportional to the area of the pool. In the case of a trapezoidal geometry, the area increases when the water level increases. For a constant discharge, this corresponds to a decelerating flow. This explains why the integrator gain tends to decrease for a decelerating flow. The reverse occurs for accelerating flow.

For the oscillating modes, we see that their damping decreases for decelerating flow, while it increases for accelerating flow. This will be confirmed by the direct study of the poles. x ∈ [0, L]. Then, the distributed transfer matrix G(x, s) is directly obtained as follows:

ŷ(x, s) q(x, s) = G(x, s) q(0, s) q(L, s)

with

g 11 (x, s) = 1 T 0 (x) φ 12 (x, s) -φ 11 (x, s) φ 22 (L, s) φ 21 (L, s) (45a)
g 12 (x, s) = 1 T 0 (x) φ 11 (x, s) φ 21 (L, s) (45b)
g 21 (x, s) = φ 22 (x, s) -φ 21 (x, s) φ 22 (L, s) φ 21 (L, s) (45c)
g 22 (x, s) = φ 21 (x, s) φ 21 (L, s) . (45d)
These expressions generalize the results already obtained in Section 3.1.5 to the non uniform flow case.

Properties of the Transfer Matrix.

4.2.1. Delays. In the non uniform flow case, the delays of the Saint-Venant transfer functions can be computed by integrating the characteristics lines along the channel. Indeed, we can show that the delays for the transfer functions p 21 (s) and p 12 (s) in the non uniform case are given by:

τ 1 = L 0 dx V 0 (x) + C 0 (x) (46)
for p 21 (s) and

τ 2 = L 0 dx C 0 (x) -V 0 (x) (47)
Author-produced version of the article published in Networks and Heterogeneous Media (NHM), This result is consistent with our physical intuition, and confirms that an open channel in non uniform flow behaves very similarly to an open channel in uniform flow: the main difference lies in the fact that the results can no longer be expressed in a closed form, but the proposed numerical scheme enables to compute the transfer function for any realistic flow configuration.

Example 7 (Delays for non uniform flow). Let us define the dimensionless time t * = αnt L , where α n stands for the uniform value of α = V 0 + C 0 . This enables us to compare the delays τ 1 and τ 2 for various non uniform flow conditions: for the uniform flow, we have τ * 1 = 1, and τ * 2 = αn βn , where the subscript n denotes the uniform flow values.

Fig. 16 depicts the way the dimensionless delays vary with the flow in both canals. The two canals have a very different behavior: the delay τ 1 varies much more with the downstream water level in the case of canal 1 than in the case of canal 2. Indeed, the backwater curve affects only a small portion of the canal 2, while it affects the whole canal 1. The delay τ 2 is much larger than τ 1 in the case of canal 2, because β is relatively much lower than α for canal 2 than for canal 1. We start from the characteristics form:

∂χ ∂t + D(x) ∂χ ∂x + E(x)χ = 0 ( 48
)
where χ(x, t) = X(x)ξ(x, t) with ξ(x, t) = (T 0 (x)y(x, t), q(x, t)) T ,

X(x) = β(x) 1 -α(x) 1 , D(x) = α(x) 0 0 -β(x) ,
and

E(x) = [X(x)B(x) -D(x)X (x)]X(x) -1 , with B(x) = 0 0 -γ(x) δ(x)
.

The elements of matrix E(x) are given by:

e 11 = 1 α + β -γ + αδ -αβ (49a)
e 12 = 1 α + β γ + βδ + αβ (49b) e 21 = 1 α + β -γ + αδ -βα (49c)
e 22 = 1 α + β γ + βδ + βα , (49d)
where we have dropped the argument x for readability, and α and β denote the derivatives of α and β with respect to x.

Applying Laplace transform to Eq. ( 48) leads to the ODE:

∂ χ ∂x (x, s) = C(x, s) χ(x, s) (50)
with C(x, s) = -D(x) -1 (sI + E(x)).

For high frequencies, the diagonal terms in C(x, s) dominate the anti-diagonal terms. Then, the high frequency approximate solution of Eq. ( 50) is given by:

χ(x, s) ≈ exp x 0 C diag (v, s)dv χ(0, s) where C diag (x, s) = -s α(x) -r 1 (x) 0 0 s β(x) + r 2 (x) , with r 1 (x) = α(x)δ(x) -γ(x) -α(x)β (x) α(x)(α(x) + β(x)) r 2 (x) = β(x)δ(x) + γ(x) + α (x)β(x) β(x)(α(x) + β(x)) .
Therefore, using the previous developments, a high frequency approximation of the poles is obtained by solving:

exp L 0 1 α(x) + 1 β(x) s + r 1 (x) + r 2 (x) dx = 1
which leads to:

p ±k ≈ - L 0 [r 1 (x) + r 2 (x)]dx τ 1 + τ 2 ± 2jkπ τ 1 + τ 2 (51)
with

τ 1 = L 0 dx α(x) and τ 2 = L 0 dx β(x)
. We recover an expression similar to the one already obtained in the uniform case. The high frequency poles have an imaginary part which is linked to the forward and backward delays, and a real part which is linked to the forward and backward damping.

Let us note that:

r 1 + r 2 = 2αβδ + (α -β)γ αβ(α + β) + α -β α + β ,
where the argument x has been dropped for readability.

Therefore, the high frequency poles damping is the sum of a positive term plus a term whose sign depends on the derivative of α(x)β(x) = 2V 0 (x) with respect to x. This means that the system will tend to be more damped if the velocity tends to increase along the canal. On the contrary, it shows that the system will tend to be more oscillating if the velocity tends to decrease along the canal. This is consistent with our intuition: we expect more damping in an accelerating flow, while a decelerating flow will be more prone to possible oscillations.

As in the uniform case, we may show that for high frequencies, the Saint-Venant equations are close to the following damped wave equation:

∂ ∂x + r 1 (x) + 1 α(x) ∂ ∂t ∂ ∂x -r 2 (x) - 1 β(x) ∂ ∂t q = 0
with boundary conditions q(0, t) and q(L, t).

The high frequency modes are therefore the result of the interaction of two gravity waves, one traveling downstream at speed α(x) = V 0 (x)+C 0 (x) with attenuation factor r 1 (x), and one traveling upstream at speed β(x) = C 0 (x) -V 0 (x) with attenuation factor r 2 (x). The high frequency estimate of the poles is the direct extension of the result obtained in the uniform case. 4.2.3. Distributed Modal Decomposition. We extend here the results obtained in section 3.2.3 to the non uniform flow case, where the distributed transfer matrix is given by Eq. ( 44).

We know using Eq. ( 51) that the high frequency poles are complex conjugate, and close to a vertical line in the left half plane. Using similar arguments as in the uniform case, the Cauchy residues theorem implies that each transfer function g ij (x, s) can be decomposed as an infinite sum:

g ij (x, s) = b ij (x) + a (0) ij (x) s + ∞ k=-∞,k =0 a (k) ij (x)s p k (s -p k ) ( 52
)
where the coefficient a

(k) ij (x) = lim s→p k (s -p k )g ij (x, s) is the residue of transfer function g ij (x, s) at the pole p k , and b ij (x) = ∂ ∂s [sg ij (x, s)]| s=0 .
Contrarily to the uniform case, we have no analytical expression for the residues a (k) ij (x). However, the proposed numerical method also enables us to compute the terms numerically. Proposition 2. In the non uniform case, the coefficient of the distributed modal decomposition can be computed as follows:

• The coefficients a (0) ij (x) are given by:

a (0) 11 (x) = - 1 T 0 (x) φ 11 (x, 0) φ21 (L, 0) (53a) a (0) 12 (x) = 1 T 0 (x) φ 11 (x, 0) φ21 (L, 0) (53b) a (0) 21 (x) = 0 (53c) a (0) 22 (x) = 0, (53d)
For canal 2, the non uniformity has an important effect on the modal decomposition. The nodes are modified following the same pattern observed for canal 1, but the change in amplitude is much more important, especially in the case of accelerating flow.

(x, s) = bij (x) + ã(0) ij (x) s + ∞ k=-∞,k =0 ã(k) ij (x)s p k (s -p k ) , (61)
where the coefficient ã(k) ij (x) = lim s→p k (sp k )g ij (x, s) is the residue of transfer function gij (x, s) at the pole p k .

For the delay-free decomposition, we have:

ã(0) ij (x) = a (0) ij (x)
, and for k = 0:

ã(k) 11 (x) = a (k) 11 (x)e x 0 dv α(v) p k (62a) ã(k) 12 (x) = a (k) 12 (x)e L x dv β(v) p k (62b) ã(k) 21 (x) = a (k) 21 (x)e x 0 dv α(v) p k (62c) ã(k) 22 (x) = a (k) 22 (x)e L x dv β(v) p k , (62d)
and the bij (x) are given by: b11

(x) = b 11 (x) + a (0) 11 (x) x 0 dv α(v) (63a) b12 (x) = b 12 (x) + a (0) 12 (x) L x dv β(v) (63b) b21 (x) = b 21 (x) (63c) b22 (x) = b 22 (x). ( 63d

Conclusion.

We have derived the Saint-Venant transfer matrix, first for the uniform flow case, and then for the general non uniform flow case. We have analyzed the transfer matrix for boundary conditions in terms of discharges, and have computed its poles, first analytically for the uniform case, then numerically for the backwater case. We also characterized the model in terms of delays, which gives an a priori evaluation of the available bandwidth for the controlled system. These results are important for the control of open channel flow. As shown in [START_REF] Litrico | Boundary control of hyperbolic conservation laws using a frequency domain approach[END_REF], the modal decomposition can be used to show that the transfer matrix belongs to the Callier-Desoer class of transfer matrices. This fact allows us to use classical frequency domain tools such as the Nyquist criteria for closed-loop stability.

Finally, we have shown that the dynamic behavior of an open channel in non uniform flow is qualitatively very similar to its behavior in uniform flow. The only difference is that in non uniform flow, the nice analytical expressions can no longer be used, but efficient numerical methods are available to compute poles, delays and modal factors. or

f ij (x, z) = f ij (x, 0) + N n=-N,n =0 a (n) ij (x) 1 z -p n + 1 p n + z 2jπ CN f ij (x, s) s(z -s) ds
Now, since |f ij (x, s)| is bounded, the integral on the right hand side tends to zero as N tends to infinity.

Finally, going back to the original transfer functions, we obtain:

g ij (x, s) = f ij (x, 0) + a (0) ij (x) s + ∞ n=-∞,n =0 sa (n) ij (x) p n (s -p n )
which is the result we wanted to prove.

Appendix B. Proof of Proposition 2.

Proof. We use in this proof the well-known fact that the state-transition matrix Φ(x, s) obtained from Eq. (41) verifies the differential equation (39).

Therefore, the φ ij (x, s) verify the following linear differential system:

where by (41), the initial condition is given by the identity matrix Φ(0, s) = I.

Let us first note that for s = 0 this system simplifies, and the solution can be expressed explicitly as follows: The solution of Eq. ( 69) is obtained as follows:

Φ (x, s) = Φ(x, s)Φ (0, s) + Φ(x, s)

x 0 Φ(v, s) -1 ∂A(v, s) ∂s Φ(v, s)dv.

By definition, we have Φ (0, s) lim Therefore the value of Φ (x, s) at x = 0 is 0 and we finally get Eq. (57). This equation provides a way to compute Φ (x, s) for any s ∈ C.

Finally, in order to compute the coefficients a The N ij (x, 0) can be computed using Eqs. (70) and the solution Φ (x, 0) given by Eq. (57) for s = 0. In that case, the expressions simplify, and lead to Eqs. (59).

To compute φ 21 (L, 0), we differentiate Eq. ( 68) with respect to s, which gives: ∂ ∂x φ 21 (x, s) = -φ 11 (x, s),

with an initial condition given by φ 21 (0, s) = 0. Therefore, once φ 11 (x, 0) is computed, φ 21 (L, 0) is directly obtained by integration, as in Eq. (60).

Therefore, we have derived a numerical means to compute all the coefficients of the modal decomposition in the non uniform case, which ends the proof.

Figure 3 .

 3 Figure 3. Backwater curves for canal 1 and 2, and various downstream boundary conditions Y 0 (L)

Figure 4 .

 4 Figure 4. Poles of canals 1 and 2 in the complex plane

Figure 5 .Figure 6 . 2 3. 1 . 5 .

 56215 Figure 5. Bode plots of transfer functions p 21 (s) and p 22 (s) for canal 1

Figure 7 .

 7 Figure 7. Spatial Bode plots of transfer function g 11 (x, s) for canal 1

Figure 8 .

 8 Figure 8. Spatial Bode plots of transfer function g 11 (x, s) for canal 2

Figure 9 .

 9 Figure 9. Function |a (k) 11 (x)| of the distributed modal decomposition for canal 1

Figure 10 .Figure 11 .

 1011 Figure 10. Function |a (k) 11 (x)| of the distributed modal decomposition for canal 2

  Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325 general. Indeed, using the definition of e A s)dv + . . . and therefore:

-Figure 12 .

 12 Figure 12. Location of the poles of Saint-Venant transfer matrix for canal 1 for non uniform flow conditions.

Figure 13 .

 13 Figure 13. Location of the poles of Saint-Venant transfer matrix for canal 2 for non uniform flow conditions.

Figure 14 . 4 . 1 . 3 .Figure 15 .

 1441315 Figure 14. Bode plots of transfer functions p 21 (s) and p 22 (s) for canal 1 in accelerating (• • • ), uniform (-) and decelerating (-• -) flow.

2 *Figure 16 .

 216 Figure 16. Variation of the dimensionless delays τ * 1 and τ * 2 as function of Y 0 (L)/Y n

Figure 17 .

 17 Figure 17. Coefficients |a(k) 11 (x)| of the distributed modal decomposition for canal 1 in accelerating (• • • ), uniform (-) and decelerating (-• -) flow.

Figure 18 .

 18 Figure 18. Coefficients |a(k) 11 (x)| of the distributed modal decomposition for canal 2 in accelerating (• • • ), uniform (-) and decelerating (-• -) flow.

Figure 19 .

 19 Figure 19. Functions b 11 (x) and b 12 (x) of the distributed modal decomposition for canals 1 and 2 in accelerating (• • • ), uniform (-) and decelerating (-• -) flow.

  s) = (α(x)β(x))s + γ(x) α(x)β(x) φ 11 (x, s) -s + δ(x) α(x)β(x) φ 21 (x, s)(66a)∂φ 12 ∂x (x, s) = (α(x)β(x))s + γ(x) α(x)β(x) φ 12 (x, s) -s + δ(x) α(x)β(x) φ 22 (x, s)(66b)∂φ 21 ∂x (x, s) = -sφ 11 (x, s)(66c)∂φ 22 ∂x (x, s) = -sφ 12 (x, s)

  ij (x). Let us introduce the modified transfer function φ21 (x, s) such that:φ21 (x, s) = 1 s φ 21 (x, s).This transfer function is well defined in s = 0 since 0 is a root of φ 21 (x, s). Using Eq.(66c), φ21 (x, s) verifies the following differential equation:∂ ∂x φ21 (x, s) = -φ 11 (x, s),(68)with an initial condition given by φ21 (0, s) = 0. Therefore φ21 (L, 0) is obtained by integrating Eq. (67a):φ21 (L, 0) = -L 0 φ 11 (v, 0)dv.Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325 Finally Eqs. (53) provide a way to compute the integrator gains a (0) ij (x) for the general case of non uniform flow. Computation of a (k) ij (x). Let us introduce the matrix Φ (x, s) = ∂ ∂s Φ(x, s), with Φ(x, s) the state-transition matrix given by Eq. (41). Then, we have: ∂Φ (x, s) ∂x = A(x, s)Φ (x, s) + ∂A(x, s) ∂s Φ(x, s)

|δs|→0Φ

  (0, s + δs) -Φ(0, s) δs = I -I δs = 0.

N

  ij (x), we need to compute the N ij (x, p k ). For any s ∈ C, the N ij (x, s) are given by:N 11 (x, s) = 1 T 0 (x) φ 12 (x, s)φ 21 (L, s)φ 11 (x, s)φ 22 (L, s)(70a)N 12 (x, s) = 1 T 0 (x) φ 11 (x, s)(70b)N 21 (x, s) = φ 22 (x, s)φ 21 (L, s)φ 21 (x, s)φ 22 (L, s)(70c)N 22 (x, s) = φ 21 (x, s)(70d)The equations simplify since by definition of a pole, we have φ 21 (L, p k ) = 0, which yields Eqs. (56). Computation of b ij (x). The coefficients b ij (x) are defined by:b ij (x) = ∂ ∂s [sg ij (x, s)]| s=0 . ij (x, s) φ21 (L, s) s=0 ,which yields Eq. (58).

Author-produced version

  of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325

Table 1 .

 1 Parameters for the two canals

	L	m B	S b	n

  -produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325

Transfer Matrix.

3.2.1. Delays.

Author

  Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325

  Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325

	4.1.2. Input-Output Transfer Matrix. Once the state-transition matrix Φ(x, s) is
	computed, the input-output transfer matrix can easily be obtained as follows:
	ŷ(0, s) ŷ(L, s)	=	p 11 (s) p 12 (s) p 21 (s) p 22 (s)	q(0, s) q(L, s)

with p 11

  , 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325 for p 12 (s).

)

  Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325

Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4(2), 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009.4.325

Acknowledgements. We would like to thank an anonymous referee for his valuable comments and suggestions.

Author-produced version of the article published in Networks and Heterogeneous Media (NHM), 2009, 4 [START_REF] Chow | Open-channel Hydraulics[END_REF], 325 -357. The original publication is available at http://www.aimsciences.org doi:10.3934/nhm.2009. 4.325 where φ21 (L, 0) is defined by: φ21 (L, 0) = -L 0 φ 11 (v, 0)dv.

(54)

• The modal coefficients a (k) ij (x) are given by:

where the N ij (x, p k ) are given by:

and φ 21 (L, p k ) is the term of the matrix Φ (x, s) given by:

• The coefficients b ij (x) are given by:

where the N ij (x, 0) are given by:

φ21 (L, 0) is given by Eq. ( 54), and φ 21 (L, 0) is given by:

Proof. See Appendix B.

Example 8 (Coefficients a

(k) ij (x) and b ij (x) for non uniform flow). The absolute value of functions a (k) 11 (x) of the distributed modal decomposition of g 11 (x, s) of the example canals 1 and 2 are depicted in Fig. 17 and18. We see that in both canals the modes in non uniform flow are similar to the ones in uniform flow. For canal 1, the spatial functions a Proof. To show that the transfer matrix G(x, s) can be expressed as an infinite sum of simple elements, we apply the residue theorem to each element of the transfer matrix. The proof is closely related to the proof of the series decomposition of cot(z) in [START_REF] Fuchs | Functions of a complex variable and some of their applications[END_REF]. To simplify the exposition, we assume that the poles have simple multiplicity but the proof can easily be adapted to the general case.

Let {C N ; N ≥ 0} a series of nested circular contours such that there are exactly two poles p N and p -N between C N -1 and C N . N is chosen sufficiently large such that the poles p N and p -N are complex conjugate.

Let us first define the function

, with a (0) ij (x) the residue of the function g ij (x, s) in zero. This function is meromorphic and can be continuously extended in s = 0 by f ij (x, 0) = d ds [sg ij (x, s)]| s=0 . We apply the Cauchy residue theorem to the function s → fij (x,s) z-s . For all N > 1, we have:

with a s). For z = 0, equation (64) leads to:

Subtracting (65) from (64) gives:


O Roger

V Galtier

P Baeten

Dardenne
email: dardenne@cra.wallonie.be

Juan Antonio Fernández Pierna

Fabian Chauchard

Sébastien Preys

Jean Michel Roger

Oswin Galtier

Vincent Baeten

Pierre Dardenne

How to build a robust model with only a few reference values: a chemometric challenge at 'Chimiométrie 2007'

published or not. The documents may come

How to build a robust model with only a few reference values: a chemometric challenge at 'Chimiométrie 2007'

Introduction

For the fourth consecutive year and following on from the success of the chemometric contests organised during previous congresses [START_REF] Dardenne | A NIR data set is the object of a Chemometric contest at 'Chimiométrie 2004[END_REF][START_REF] Fernández Pierna | Chemometric contest at 'Chimiométrie 2005': a discrimination study[END_REF][START_REF] Fernández Pierna | Soil parameter quantification by NIRS as a Chemometric challenge at 'Chimiométrie 2006[END_REF], another dataset was proposed for the 'Chimiométrie 2007' meeting (http://www.chimiometrie.fr/) held in Lyon, France (29-30 November 2007). As for the first contest in 2004, this dataset was selected to test the ability of participants to apply regression methods to NIR data. The aim of Challenge 2007 was to perform a calibration model as robust and precise as possible when only a few reference values were available. Regression algorithms, as PLS (Partial Least Squares) are very sensible to outlying observations, which are typically expected to be present in experimental data. One of the solutions for this drawback of the classical regression techniques is the possibility of constructing robust versions of the regression algorithms [START_REF] Huber | Robust Statistics[END_REF][START_REF] Wakeling | [END_REF][START_REF] Griep | Comparison of semirobust and robust partial least squares procedures[END_REF][START_REF] Wakelinc | A robust PLS procedure[END_REF][START_REF] Serneels | Partial robust Mregression[END_REF][START_REF] Daszykowski | Robust statistics in data analysis -a review. Basic concepts[END_REF]. In the early 1990s, Dardenne et al. [START_REF] Dardenne | Reduction of wet chemistry for NIR calibrations[END_REF] studied the effects of various properties such as moisture, particle size and ambient temperature on NIR calibration models in order to reduce the wet chemistry needed for them. Artificial spectral variations were created by changing the moisture and particle size of wheat samples when predicting the protein content. These samples, measured at different temperatures on a monochromator, produced wide spectral variations that helped to develop robust models. The data used in this challenge come from that paper.

Nine participants took up the challenge with the proposed data. The results were evaluated on the basis of the best validation criteria (R 2 and RMSEP) obtained for the predicted values of the test set, and also on the quality of the approach from a methodological perspective. The three best approaches were presented during the congress and are summarised here, together with the approach put forward by the organizers of Challenge 2007.

Material and Methods

Several datasets were provided to the participants: a calibration dataset, an experimental design dataset, a standard replicate dataset and a test dataset.

Calibration Dataset

For this challenge, only 10 spectra of ground wheat acquired using a FOSS NIRSystems 4500, measured between 1300 nm and 2398 nm each 2 nm, were supplied, together with the protein content of the 10 samples measured in g/Kg Dry Matter (DM). The aim was to perform a calibration model as robust and precise as possible with the available 10 reference values for protein content.

Experimental Design Dataset

There were also the spectra from an experimental design based on other 11 samples of ground wheat (with unknown reference values). These 11 whole grain samples were separated into two homogeneous sets: one was dried to reduce the moisture content to +/-9%, and the other was moistened to reach 13-14% humidity. Each grain sample set was then divided again in two groups: the first subsample was ground finely (Cyclotec apparatus, level C) and the second one was ground more coarsely (Ika apparatus, level I). A total of 11*4 sample sets was thus available. These samples were measured by reflection NIR in duplicate at three room temperatures (18°C, 23°C and 27°C). The experimental design database therefore contained 11×4×3×2 = 264 spectra.

Standard Replicates Dataset

Additional spectral information was supplied from one set of 10 samples (sealed cells) scanned on 31 different instruments of the same type and from a second set of 10 scanned on 17 instruments. A total of 480 spectra was thus available. The reference values corresponding to the experimental design samples and to the 10 samples measured on the different instruments were unknown.

Test Dataset

Some 2,000 spectra from routine analyses were acquired from +/-10 different instruments with several levels of granulometry, humidity and temperature, and provided by the Requasud network (http://www.requasud.be/) from 1991 to 2007. The 2,000 spectra each had a reference protein value obtained by the reference method, but these values were not communicated to the participants.

The calibration dataset was kept in matrices X cal , y cal whereas the data from the experimental design in a matrix X ed and the standard replicates dataset in two matrices X instr1 (31 instruments) and X instr2 (17 instruments). The Test dataset was kept in X test .

Results

Participant 1

Two approaches were tested: exhaustive calibration on artificial data, and Orthogonal Signal Correction (OSC) [START_REF] Wold | Orthogonal signal correction of near-infrared spectra[END_REF] on these artificial data.

Approach 1: Exhaustive calibration on artificial data

This approach involved using information about perturbation and 'injecting' it into the 10 spectra of the calibration dataset (X cal ) in order to generate an artificial calibration database X artif . As this new calibration database included all possible perturbations, PLS would automatically be able to find the most robust direction for prediction.

In order to do so, initially all the possible perturbating vectors were identified using the experimental design dataset (X ed ) and the standard replicates dataset (X instr1 + X instr2 ) for each of the 24 spectra (X sample i ) of the same sample i (one sample i had 24 'j' combinations of possible perturbations). The 24 vectors of perturbations were calculated as follows: _δ δ δ δ(x i,j ) = x i,j -xi (1) where ¯xi is the mean spectra of sample I and x i,j is the j th spectra of sample i. All the perturbations were put in a new matrix D ed containing the 264 possible perturbations. The same procedure was followed in order to identify all possible perturbations from the instruments, resulting in a matrix D instr of 480 spectra. Figure 1 shows an example of the perturbation spectra that could be obtained. Finally, all the perturbations from D ed and D instr were combined (using addition), resulting in a matrix of 127,464 vectors ((264x 480) + 264 + 480)). In order to generate an artificial database, 20,000 perturbations were randomly taken and added to 2,000 repetitions of the 10 spectra from X cal , resulting in a matrix X artif . The reference values corresponding to this matrix y artif were obtained from the 2,000 repetitions of the reference value for the calibration dataset y cal , because the perturbations added to the spectra should not change the reference value. Figure 2 shows that the spectral space of the calibration database increased when comparing from X cal to X artif using PCA. This was due to the integration of perturbation into the database. The model was then calibrated using X artif and a second derivative (window 9, polynomial order 3) as pre-processing and wavelength selection (from 1950 to 2250 nm). Some non-linearity remained, but the model in cross-validation gave interesting performances. (Figure3).

Approach 2: OSC approach on artificial data

Orthogonalisation techniques are an interesting alternative for building robust models. They include the External Parameter Orthogonalisation (EPO) [START_REF] Roger | EPO-PLS external parameter orthogonalisation of PLS : Application to temperature-independent measurement of sugar content of intact fruits[END_REF] approach (see participants 2 and 3). The Orthogonal Signal Correction (OSC) approach takes account of the reference values [START_REF] Wold | Orthogonal signal correction of near-infrared spectra[END_REF]. OSC based on an experimental design was proposed to take advantages of both EPO and OSC [START_REF] Preys | Robust calibration using orthogonal projection and experimental design. Application to the correction of the light scattering effect on turbid NIR spectra[END_REF]. In this application, however, although perturbation variability was present in the experimental designs, there was no variability in the reference values. The OSC approach was therefore proposed as a pre-treatment on the previously built artificial calibration database, where reference values and perturbation variability were present.

The other pre-treatments used were a second derivative using the Savitsky and Golay algorithm ( [START_REF] Gorry | General least-squares smoothing and differentiation by the convolution (savitzky-golay) method[END_REF]) with a window of 9 and a polynomial of degree 3, followed by an SNV on the whole wavelength range. Cross-validation optimised the calibration model with only two latent variables and one OSC-factor removed, giving a parsimonious model.

Participant 2

The method chosen for calibration was PLS. The first optimisation of the model was carried out by testing different pre-processing methods, using the calibration error (SEC) and the cross-validation leave-one-out error (SECV) as selection criteria. Because the spectra were affected by a baseline and a multiplicative effect, the pre-treatment used was a second derivative using the Savitsky and Golay algorithm with a window of 9 and a polynomial of degree 3 ([14]), followed by SNV ( [START_REF] Barnes | Standard normal variate transformation and de-trending of near-infrared diffuse reflectance spectra[END_REF]), as explained in [START_REF] Davis | Back to basics: removing multiplicative effects (1)[END_REF]. Two models were retained: one was based on the raw spectra, and the other was based on the pre-processed spectra, as previously described.

Both models were then submitted to a test for each factor, independently from other factors, in order to determine its influence. To do this, the matrices of experiments were averaged according to the factors not involved, in order to retain the variability only for the matrix whose effect needed to be evaluated. Both models were then tested on this reduced matrix. For example, to study the influence of moisture, the experimental design was averaged in terms of granulometry, temperature and repetitions, which supplied 22 spectra: 11 for the dry products and 11 for the humid products. The predictions for each sample at two moisture contents were then compared.

The method used to improve the robustness of the calibration was based on EPO ( [START_REF] Roger | EPO-PLS external parameter orthogonalisation of PLS : Application to temperature-independent measurement of sugar content of intact fruits[END_REF][START_REF] Andrew | Transfer by orthogonal projection: making nearinfrared calibrations robust to between-instrument variation[END_REF]). This involves the orthogonalisation of the measurement space of the spectra, on the basis of the disturbances caused by the factors. The orthogonalisation removes the components of the sub-space spanned by the differences between the spectra repeated for the same sample. On one hand, the more the number of removed components increases, the more spectra of each sample become similar and then the calculated predictions on these spectra become more similar. On the other hand, if too many components are removed, this could alter the calibration. To choose the dimension of the space to be removed by orthogonalisation, we used the Wilks's lambda, which represents the ratio of the variance inter samples and the variance intra samples. This indicator was calculated in accordance with the number of removed factors and the number of latent variables of the model. The examination of its evolution enabled us to choose the optimal dimension to be removed. Finally, the orthogonalisation of the models was achieved in two ways:

-By individual orthogonalisation, i.e. bringing together the six projectors (4 for X ed , 1 for X instr1 and 1 for X instr2 ) calculated individually for every size of influence -By a global orthogonalisation, i.e. calculating a global projector, on all the X matrices brought together. These two methods were applied to both models (raw and pre-treated), resulting in four different predictions for the test set. The three most different predictions were retained.

As previously explained, the models were submitted to a test for each factor independently to determine its influence. The sensitivity of these models to the different size factors is illustrated in Table 1 (the left side corresponds to the raw samples and right side to the pre-treated data). As shown in the table, the factors have varying influences on the robustness of the model: -For moisture, there is clearly a positive effect of the pre-treatment: the bias is far less strong, as is the dispersion about the calibration line. However, the slope is far from 1 (+/-0.8) in both cases. -For granulometry, the pre-treatment reduces the bias and the dispersion about the calibration line. -Temperature and sample repetition seem to be slightly influential. The pretreatment improves the prediction, independently of the robustness.

The sensitivity of the models, after orthogonalisation, to the different size factors is illustrated in Table 2. It shows the clear advantage of orthogonalisation, especially for the most influential factors. The number of components that need to be removed is 3 for the factors of influence in the first experimental design, and 8 for the inter spectrometer repetitions.

Participant 3

As a first step, the goal was to show that the external parameters (e.g., moisture, grinding, temperature, and repeatability on different NIR apparatus) were significant. To demonstrate external parameter importance, PCA was performed after centering raw spectra. For each PCA, clusters were made up of each external parameter (Figure 4). different NIR apparatus. A pre-processing method was then applied, its aim being to remove from the data space the part that was most influenced by the external parameter variations. The method proposed was EPO [START_REF] Roger | EPO-PLS external parameter orthogonalisation of PLS : Application to temperature-independent measurement of sugar content of intact fruits[END_REF], which estimates the parasitic subspace by computing a PCA on a set of spectra measured on the same objects, while the external parameter varies.

Different pre-treatments were applied: SNV [START_REF] Barnes | Standard normal variate transformation and de-trending of near-infrared diffuse reflectance spectra[END_REF] (standard normal variate), MSC [START_REF] Isaksson | The Effect of Multiplicative Scatter Correction (Msc) and Linearity Improvement in Nir Spectroscopy[END_REF] (multiplicative scatter correction), first derivative spectra using Savistsky Golay algorithm, and raw spectra. The first derivative spectra pre-treatment gave the best result.

As for participant 2, when applying EPO pre-processing, there are two possibilities: by individual orthogonalisation, i.e. eliminating external parameters on an ad hoc basis and by a global orthogonalisation. Both possibilities were tested and the second method was chosen because it gave the best PCA representation in that example. Figure 5 shows the two first scores of PCA with EPO pre-processing on the first derivative spectra and on raw spectra. Comparing this with Figure 4.d clearly shows the powerful performance of EPO pre-processing. The PLS model was calibrated on the calibration matrix after EPO pretreatment of the first derivative spectra. The Jack-Knife technique [START_REF] Riu | Jack-knife technique for outlier detection and estimation of standard errors in PARAFAC models[END_REF] was used to fix the required number of factors for model construction. Cross-validation was applied in regression, so the optimal factor number was determined based on the prediction of the sample kept out of the individual model. A final model with two latent variables and an R 2 of 0.99 and RMSEC (root mean square error in calibration) of 3.30 was applied to the test dataset matrix. Figure 6 shows the spectra of the 10 available samples. The strategy proposed sought to create a huge calibration set by taking account of the experimental design and the repetitions between instruments. The main idea was to compute all the differences between all the pairs of the same sample in both the experimental design and the repetitions, and then add these differences to the 10 spectra with known Y values. The procedure can be described in two steps:

Challenge organisers

Step 1 -The experimental design (X ed ) comprised 11 samples x 2 moistures x 2 grinders x 3 temperatures x 2 replicated scans. In total, 264 spectra were available. All the differences between the pairs of samples were computed, giving a total of 552 (24 * (24-1)) differences per sample. Then, for the 11 samples, a matrix containing 6,072 differences was retained. Figure 7 shows all these differences. Step 2 -The differences between the two series of 10 samples scanned on different instruments were computed. For set 1 (X instr1 ) scanned on 31 instruments, 930 (31* (31-1)) differences were calculated for each sample. For set 2 (X instr2 ) scanned on 17 instruments, 272 (17* (17-1)) differences for each sample were obtained. For the 10 samples, a total of 12,020 differences were computed, as shown in Figure 8. All these differences were added to the 10 available spectra (X cal ), giving a total of 180,920 spectra with only 10 reference values. Due to computer and time limitations, a random selection of spectra was made and the final model was constructed using 2,548 spectra and 10 reference values. Figure 9 shows the projection of the 2,000 samples of the test dataset (X test ) on the 2,548 spectra set constructed, indicating that both sets match perfectly. Figure 10 provides a comparison in terms of RMSEC and RPD between the model using only the 10 available samples (Figure 10 a) and the model constructed using the 2,548 spectra obtained following the strategy described (Figure 10 c). In both cases, the prediction for the test dataset is also shown (Figures 10 b andd respectively).

Final results

The evaluation of the approaches was based on the best results obtained for the predicted values of the test dataset (2,000 spectra). The reference protein value obtained by the reference method for these spectra was not communicated to the participants. For the evaluation, the R 2 and RMSEP (Root Mean Square Error in Prediction) were used as validation criteria. An additional parameter, RPD (ratio of the standard deviation of the population over the standard error of prediction), was also included. The results for the different approaches are summarised in Table 3.

Conclusion

The challenge produced a wide diversity of results. However, when evaluating the 10 answers received (nine from the participants, one from the organisers) as a whole, Challenge 2007 showed that it is still possible to perform a robust and precise calibration when only a few reference values are available. The use of well-defined experimental designs, including repeated measurements under different conditions, will lead to the use of simple, easy and low-cost instruments and the construction of robust models.

During the congress the approaches summarised here were presented, together with the challenge organisers' approach. The participants found the results interesting and it was decided to include another challenge for the next congress.

The data of this challenge and previous challenges is available on the internet address of the French Chemometric Society (http://www.chimiometrie.fr/).

Figure 1 :

 1 Figure 1: Examples of perturbation spectra from a) the experimental design dataset and b) the standard replicated dataset.

Figure 2 :

 2 Figure 2 : PCA plot showing the spectral space of the increased calibration dataset (X artif ) including the original X cal .

Figure 3 :

 3 Figure 3 : Final model constructed using X artif for participant 1.

Figure 4 :

 4 Figure 4 : PCA scores calculated without pre-processing, showing importance of external parameters: (a) temperature, (b) moisture, (c) grinding (d) repeatability ondifferent NIR apparatus. A pre-processing method was then applied, its aim being to remove from the data space the part that was most influenced by the external parameter variations. The method proposed was EPO[START_REF] Roger | EPO-PLS external parameter orthogonalisation of PLS : Application to temperature-independent measurement of sugar content of intact fruits[END_REF], which estimates the parasitic subspace by computing a PCA on a set of spectra measured on the same objects, while the external parameter varies.Different pre-treatments were applied: SNV[START_REF] Barnes | Standard normal variate transformation and de-trending of near-infrared diffuse reflectance spectra[END_REF] (standard normal variate), MSC[START_REF] Isaksson | The Effect of Multiplicative Scatter Correction (Msc) and Linearity Improvement in Nir Spectroscopy[END_REF] (multiplicative scatter correction), first derivative spectra using Savistsky Golay algorithm, and raw spectra. The first derivative spectra pre-treatment gave the best result.As for participant 2, when applying EPO pre-processing, there are two possibilities: by individual orthogonalisation, i.e. eliminating external parameters on an ad hoc basis and by a global orthogonalisation. Both possibilities were tested and the second method was chosen because it gave the best PCA representation in that example. Figure5shows the two first scores of PCA with EPO pre-processing on the first derivative spectra and on raw spectra. Comparing this with Figure4.d clearly shows the powerful performance of EPO pre-processing.

Figure 5 :

 5 Figure 5 : Two first component scores of PCA with EPO pre-processing: (a) on 1st derivative spectra, (b) on raw spectra.

Figure 6 :

 6 Figure 6 : Spectra of the 10 available samples (X cal ).

Figure 7 :

 7 Figure 7 : Differences between the pairs of samples based on the experimental design dataset.Step 2 -The differences between the two series of 10 samples scanned on different instruments were computed. For set 1 (X instr1 ) scanned on 31 instruments, 930 (31* (31-1)) differences were calculated for each sample. For set 2 (X instr2 ) scanned on 17 instruments, 272 (17* (17-1)) differences for each sample were obtained. For the 10 samples, a total of 12,020 differences were computed, as shown in Figure8.

Figure 8 :

 8 Figure 8 : Differences between the two series of samples scanned based on the Standard Replicates Dataset.

Figure 9 :

 9 Figure 9 : PC1 vs. PC2 showing the projection of the 2,000 samples of the X teston the 2,548 spectra dataset constructed by the challenge organisers. The proposed model was constructed using least-squares support vector machines (LS-SVM)[START_REF] Cogdill | Least-squares support vector machines for chemometrics: an introduction and evaluation[END_REF]. The results are shown in Figure10.

Figure 10 :

 10 Figure 10 : Results obtained after the application of the SVM model for the challenge organisers. a) model using only the 10 available samples of X cal ; b) prediction of X test based on model represented in a; c) model constructed using the 2,548 spectra obtained following the strategy described and d) prediction of X test based on model represented in c.

Table 1 :

 1 Validation criteria of the models after test for each factor individually (the left side corresponds to the raw samples and right side to the pre-treated data) for participant 2 :

				Raw data					Pretreated data
		R2	Bias	SEP	Slope	Offset	R2	Bias	SEP	Slope	Offset
	Humidity	0,97	-26,90	7,10	0,82	2,85	1,00	1,93	5,66	0,80	27,60
	Granulometry	0,99	-11,70	3,00	1,02	-14,70	0,99	-2,77	2,34	0,97	1,75
	Temperature	1,00	-0,12	1,61	0,97	3,46	1,00	-0,33	1,17	1,02	-2,86
	Repetition	1,00	0,21	1,00	1,01	-1,42	1,00	-0,06	0,52	1,00	-0,68
	S1	0,83	-	5,73	-	-	0,91	-	4,12	-	-
	S2	0,95	-	2,95	-	-	0,96	-	2,55	-	-

Table 2 :

 2 Validation criteria of the models after test for each factor individually after orthogonalisation (the left side corresponds to the raw samples and right side to the pre-treated data) for participant 2.

				Raw data					Pretreated data
		R2	Bias	SEP	Slope	Offset	R2	Bias	SEP	Slope	Offset
	Humidity	0,99	0,05	2,51	1,01	-1,60	1,00	0,13	1,47	1,00	0,75
	Granulometry	0,99	-0,27	2,46	1,03	-4,85	1,00	0,00	1,40	1,01	-1,68
	Temperature	1,00	0,01	1,37	0,98	2,02	1,00	0,00	0,58	1,00	-0,44
	Repetition	1,00	-0,20	0,47	1,00	-0,73	1,00	-0,07	0,28	1,00	-0,25
	S1	0,98	-	1,77	-	-	0,98	0,00	1,69	-	-
	S2	0,99	-	0,96	-	-	1,00	-	0,66	-	-

Table 3 :

 3 Summary of the results of the different approaches in terms of R 2 , RMSEP and RPD (ratio of the standard deviation of the population over the standard error of prediction).

		R2	RMSEP	RPD
	Participant 1 (approach 1)	0,91	3,91	3,32
	Participant 1 (approach 2)	0,88	4,51	2,88
	Participant 2	0,89	4,24	3,05
	Participant 3	0,86	4,86	2,65
	Challenge organizer	0,93	3,56	3,82

Acknowledgments

We would like to thank all the participants who spent time working on the data and presenting their results. Apart from the authors of this paper, the other participants in Challenge 2007 were Marion Cuny from Eurofins, Jean-Claude Boulet from INRA, Abdelaziz Faraj from IFP, Dominique Bertrand from INRA, Frédéric Estienne from Sanofi Aventis and Ludovic Duponchel from LASIR.