R Ravaud

G Lemarquand
email: guy.lemarquand@ieee.org

T Roussel

V Lemarquand

Ranking of the nonlinearities of electrodynamic loudspeakers December 3, 2009

Keywords: electrodynamic loudspeaker, nonlinearities, small signal parameters, nonlinear modeling

The aim of this paper is to present a way of ranking the nonlinearities of electrodynamic loudspeakers. For this purpose, we have constructed a nonlinear analytic model which takes into account the variations of the small signal parameters. The determination of these variations is based on a very precise measurement of the electrical impedance of the electrodynamic loudspeaker. First, we present the experimental method to identify the variations of these parameters, then we propose to study theoretically the importance of these nonlinearities according to the input level or the input frequency. We show that the parameter which creates the most distortions is not always the same and depends mainly on both the input level and the input frequency. Such results can be very useful for optimizing electrodynamic loudspeakers.

INTRODUCTION

The reference model describing the electrodynamic loudspeaker was designed by Thiele and Small [START_REF] Thiele | Loudspeakers in vented boxes: Part 1 and 2[END_REF]. Their model was adapted to describe an electrodynamic loudspeaker as a linear system. This model is very useful because it is very simple to use and its parameters can be conveniently presented in terms of an electric analog circuit [START_REF] Klippel | Dynamic measurement and interpretation of the nonlinear parameters of electrodynamic loudspeakers[END_REF]. However, a loudspeaker shows nonlinearities that produce distortions. These nonlinearities have three major sources [3][4]: the suspensions, the diaphragm [START_REF] Suzuki | Radiation and diffraction effects by convex and concave domes[END_REF]- [START_REF] Lemarquand | Large bandwith loudspeaker emitting coherent acoustic waves: nonlinear inter-modulation effects[END_REF] and the motor [START_REF] Dobrucki | Nontypical effects in an electrodynamic loudspeaker with a nonhomogeneous magnetic field in the air gap and nonlinear suspension[END_REF]- [START_REF] Ravaud | Time-varying non linear modeling of electrodynamic loudspeakers[END_REF].

They have been largely studied and there are many attempts to model them [START_REF] Kaizer | Modeling of the nonlinear response of an electrodynamic loudspeaker by a volterra series expansion[END_REF]- [START_REF] Vanderkooy | A model of loudspeaker driver impedance incorporating eddy currents in the pole structure[END_REF].

However, most of the papers dealing with the nonlinearities of electrodynamic loudspeakers do not take into account the variations of the mechanical damping R ms or the equivalent mass M ms . But these parameters are also nonlinear and their variations must be taken into account in order to precisely characterize the distortions created by electrodynamic loudspeakers.

It can be noted that alternative loudspeakers structure without iron have been proposed and studied [START_REF] Lemarquand | Ironless loudspeakers[END_REF]- [START_REF] Remy | Ironless and leakage free voice-coil motor made of bonded magnets[END_REF] which use ferrofluid suspensions in order to delete the classical distortions generated by electrodynamic loudspeakers [START_REF] Ravaud | Modelling an ironless loudspeaker by using three-dimensional analytical approaches[END_REF]- [START_REF] Ravaud | Magnetic pressure and shape of ferrofluid seals in cylindrical structures[END_REF]. This paper has two objectives. In the first part, we present an experimental method which allows us to determine precisely the variations of the small signal parameters according to the input current. Our experimental method is based on a very precise measurement of the electrical impedance of the electrodynamic loudspeaker. Indeed, the measurement accuracy is about 10 -4 Ohm for both the real and imaginary parts of the electrical impedance.

However, all our approach does not use the classical definition of an electrical impedance. An electrical impedance, commonly defined by its ratio U I (in the frequency domain) does not depend on input current. Indeed, this electrical impedance is only defined for linear systems. However, our experimental measurement shows that this ratio U I depends on I. This can be easily verified experimentally by measuring the experimental impedance for different input currents. Such a result is in fact known for scientifics involved in the modeling of nonlinear systems: a nonlinear system depends generally on the input level. In the case of an electrodynamic loudspeaker, the electrical impedance, that is to say, the ratio U I depends on the input level.

We also show that it exists a bijective relation between the input current I and the position of the voice coil X. This result is of great importance because it shows that the non linear parameters describing an electrodynamic loudspeaker (Bl(x), R ms (x) and L e (x)) can also be described as parameters depending on input current (in the time-domain or in the frequency-domain).

It can be noted that, strictly speaking, the force factor Bl which is commonly used in analog circuits represents the mean induction field times the length l of the voice coil inside the air gap. Therefore, when the voice coil position

x(t) is sinusoidal, in complex notations, the force factor should be seen as a mean force factor Bl.

In the second part, we discuss the behaviour of each nonlinearity accord-ing to both the input level and the input frequency. It is noted that all this part is treated theoretically. Indeed, a very good agreement between our analytical model and experimental measurements has confirmed that our model can be used for modelling the distortions created by an electrodynamic loudspeaker. We show here that the lumped parameter whose relative variation according to input current is the most important is not always the same and depends greatly on both the input level and the input frequency.

Another drawback known in electrodynamic loudspeakers is that it is a time-varying system. In fact, the electrical resistance increases with time and the mechanical compliance of the outer rims depends also on time. These dependence have been studied with a similar approach in previous papers.

In this paper, some care has been taken in the experimental measurements for omitting these temporal effects.

This paper proposes a way of ranking the nonlinearities of electrodynamic loudspeakers. To our knowledge, this way of characterizing these nonlinearties has never been used for electrodynamic loudspeakers. However, many new phenomena can be seen by using such an experimental approach.

This is why such results are very interesting for many manufacturers involved in the design of electrodynamic loudspeakers. We precise here that we use for the rest of this paper the small notations x(t) and i(t), the displacement of the voice coil and the input current in the voice coil in the time-domain and the capital letters (X, I, U ) in the frequency domain. We also precise that U is the input voltage.

2 Classical description of a loudspeaker and its limits

2.1 The small signal model using lumped parameters

A common way of characterizing an electrodynamic loudspeaker is to measure its electrical impedance Ze. Its theoretical expression is well known and is given by Eq. [START_REF] Thiele | Loudspeakers in vented boxes: Part 1 and 2[END_REF].

Z e = R e + jL e wR µ jL e w + R µ + Bl 2 R ms + jM ms w + 1 jCmsw (1)
where all the parameters are defined below. However, strictly speaking, the term -electrical impedance-should not be used to describe a nonlinear system like an electrodynamic loudspeaker.

As a consequence, the lumped parameters defined in Eq. (1) depend on the input current. The main problem is thus to know how to characterize their dependence according to the input current. According to the current state of art, cited in the introduction of this paper, the dominant nonlinearities in electrodynamic loudspeakers are Bl(x), C ms (x) and L e (x). However, it exists a relation between the voice coil position x and the input current i.

This relation is shown in next section. Consequently, it can be noted that these nonlinearities can depend on x or i: it is in fact equivalent. In this paper, we choose to work directly with the input current i because Z (N L) e (I) is easier to determine experimentally. Consequently, thanks to this very precise experimental measurement, the variations of the lumped parameters can be determined according to the input current i. It is noted here that the aim of this section is merely to show that it exists a bijective relation between the input current I and the voice coil position X.

To find this relation, we use one differential equation describing the electrodynamic loudspeaker and we use complex notations. We find:

M ms (jw) 2 + R ms (jw) + 1 C ms X = BlI (2)
So, we can write:

X = Bl M ms (jw) 2 + R ms (jw) + 1 Cms I (3)
Therefore, we see that Eq. ( 3) is the bijective relation between the input current I and the voice coil position X, in complex notation. This relation shows that we can work either with the input current I or with the voice coil position X so as to describe the nonlinear variations of the Thiele and Small parameters.

Relation between the input current I and the voice coil position X according to the nonlinear model used in this paper

Strictly speaking, the relation between the input current I and the voice coil position is nonlinear. Indeed, all the small signal parameters depend on the input current. It is emphasized here that the transfer function describing the ratio U I is considered as the ratio of the part of response with fundamental frequency w to the excitation. We call it a nonlinear electric impedance in this paper but it is useful to precise that nonlinear systems cannot be described with complex notations.

We describe in the next section a way of determining the nonlinear laws of all the small signal parameters. Therefore, we precise here for the rest of this paper that (3) can be written in the frequency-domain as follows(4):

X = Bl(I) M ms (jw)(I) 2 + R ms (I)(jw) + 1 Cms(I) I (4)
It is noted that the determination of the small signal parameters can be determined in the frequency domain. We explain in next section how to find these variations.

3 Determination of the nonlinear variations of the lumped parameters

Principle of the measurement

We have shown in the previous section that the nonlinear variations of the lumped parameters describing an electrodynamic loudspeaker depend on the input current. We explain now how to find these variations.

The electrodynamic loudspeaker used is a boomer (mark: Eminence, number: beta 15 ). We use the electrical impedance of the electrodynamic loudspeaker because its measurement is very precise. Our measurement device is a Wayne Kerr bridge which has an excellent precision (10 -4 Ω). This experimental devices is dedicated to the impedance measurement and functions as a classical impedance bridge.

In order to measure the electrical impedance of a loudspeaker, it is mounted in a normalized plane in an anechoic room [23]. Dalmont [START_REF] Dalmont | Radiation impedance of tubes with different flanges:numerical and experimental investigations[END_REF] has shown that the acoustical impedance in these conditions is the same as the one when the loudspeaker is mounted in an infinite baffle in an anechoic room. We measure the electrical impedance by varying the frequency and the coil current.

So, we build an experimental impedance layer by using a Runge Kutta algorithm to determine all the experimental measured points. In other words, this measurement algorithm chooses the best measurement point according to the gradient method. Such a measurement method allows us to detect all the subtle effects due to the nonlinearities in electrodynamic loudspeakers.

This algorithm has been improved in relation to the one used in a previous paper. In the previous paper, the algorithm took some experimental points by using constant intensities. In other words, our measurement system used the algorithm of gradient only in two dimensions for one intensity and then it was repeated for another intensity. For example, an intensity was fixed (for example 50 mA), and a two-dimensional algortihm allowed us to determine the measurement points, that is to say, our algorithm automatically settled one intensity and took some experimental points by using a method based on the gradient method. In this paper, the algorithm uses a Runge Kutta algorithm in three dimensions (Z e ,I,f ) and we can rapidly obtain an impedance layer which is very precise. Consequently, the temporal effects can be neglected. A two-dimensional representation of the electrical impedance is shown in Fig. 2: the imaginary part is a function of the real part ( it is Nyquist diagram for different coil currents). We can say that the Wayne Kerr bridge cannot supply currents greater than 0.2A. Consequently, the parameter variations are determined in this current interval.

Figure 3 confirms that the electrical impedance is a function of the input current.

Determination of the nonlinear electrical impedance

The next step is thus to determine the dependence of the small signal parameters with the coil current so as to construct the nonlinear transfer function.

Bl(I) = Bl(1 + µ Bl I + µ 2 Bl I 2 + ... + µ n Bl I n ) (5) R ms (I) = R ms (1 + µ Rms I + µ 2 Rms I 2 + ... + µ n Rms I n ) (6) C ms (I) = C ms (1 + µ Cms I + µ 2 Cms I 2 + ... + µ n Cms I n ) (7)
M ms (I) = M ms (1 + µ Mms I + µ

Mms Mms0

Figure 8: The ratio Mms M ms0 is a function of the coil current ( x : coil current 0A 0,2A)(f=100Hz);M ms0 = 0.009kg

Discussion

The most sensitive parameter to the variations of the coil current is the equivalent damping parameter R ms . The parameter which is the less sensitive to the variations of the coil current is the inductance L e of the voice coil.

However, though this inductance sensitiveness is rather weak, the relative variation of the inductance L e is rather important. In other words, even if the inductance L e is a parameter that contributes just a little to the variation of its electrical impedance, we see that its relative variation is important.

Another interesting point can be adressed. With the values found in this paper, at 200Hz, we see that the eddy current resistance R µ equals roughly wL e . We can conclude that the effects of the eddy currents become important from 200Hz for this electrodynamic loudspeaker.

5 Study of the distortions created by the electrodynamic loudspeaker

Obtaining the nonlinear differential equation

The nonlinear differential equation describing the electrodynamic loudspeaker is given by Eq. ( 12). We use the parameter k = 1

Cms to solve it.

a(i) d 3 x(t) dt 3 + b(i) d 2 x(t) dt 2 + c(i) dx(t) dt + d(i)x(t) = u(t) (12)
with

a(i) = M ms (1 + µ Mms i + µ 2 Mms i 2 ) (Le(1 + µ Le i + µ 2 Le i 2 )) Bl(1 + µ Bl i + µ 2 Bl i 2 ) (13) b(i) = (M ms (1 + µ Mms i + µ 2 Mms i 2 )Re Bl(1 + µ Bl i + µ 2 Bl i 2 ) + R ms (1 + µ Rms i + µ 2 Rms i 2 )R e Bl(1 + µ Bl i + µ 2 Bl i) (14)
c(i) = Re R ms (1 + µ Rms i + µ 2 Rms i 2 ) + (Bl(1 + µ Bl i + µ 2 Bl i 2 )) 2 Bl(1 + µ Bl i + µ 2 Bl i 2 ) + (Le(1 + µ Le i + µ 2 Le i 2 ))k(1 + µ k i + µ 2 k i 2 ) Bl(1 + µ Bl i + µ 2 Bl i 2 ) ( 15
)
d(i) = Re (k(1 + µ k i + µ 2 k i 2 )) Bl(1 + µ Bl i + µ 2 Bl i 2 ) (16)

Solving the nonlinear differential equation

The relation between the coil current i and the position of the voice coil x is nonlinear (Eq. ( 4)). Each small signal parameter is a function of the voice coil position and we obtain, for example, for one parameter :

R ms (x) = R ms (1 + μRms x + μ2 Rms x 2 ) ( 17
)
with μ2

Rms << |μ Rms | << 1 .

If we assume that the input voltage u(t) is sinusoidal, in this case, the solution of the nonlinear differential equation is a trigonometric expansion. The solution can be developed until the order 2 (µ 2 ).

x(t) = x 0 (t) + µx 1 (t) + µ 2 x 2 (t) + ... (18)
where x 0 (t) is the solution of the nonlinear differential equation of the electrodynamic loudspeaker when we neglect the terms with orders higher than zero, x 1 (t) is the solution of the nonlinear differential equation when we omit the terms with orders higher than one and smaller than one, x 2 (t) is the solution of the nonlinear differential equation of the electrodynamic loudspeaker when we neglect the terms with orders smaller than two and higher than two.

It can be noted that some methods like Volterra Series are interesting but do not show what parameters are really nonlinear. The way of solving the nonlinear differential equation is very important but we think that it is not the most important thing for characterizing electrodynamic loudspeakers.

The real problem is to know what parameters vary and how. Furthermore, the ranking of these parameters according to the input current or frequency is of great importance. Indeed, we can think that if the structure of an electrodynamic loudspeaker must be improved, we must know all the defects in the motor or the suspensions.

Let us consider now Eq. ( 12). The denominator in Eq. ( 12) contains a nonlinear term. Consequently, it is very difficult to solve this nonlinear differential equation with a nonlinear denominator. One possible solution is to approximate this denominator as follows :

1 Bl(x) = Bl 0 + Bl 1 x + Bl 2 x 2 + ... (19)
The previous relation is used as a simplification for solving numerically the nonlinear differential equation. Moreover, we can use a classical trigonometric expansion to solve Eq. [START_REF] Gander | Moving-coil loudspeaker topology as an indicator of linear excursion capability[END_REF]. In short, the solution of the nonlinear differential equation of the electrodynamic loudspeaker is

x(t) = A cos(wt) + B sin(wt) + C cos(2wt) + D sin(2wt) + ... (20)
All the terms A,B,C,D,... are found numerically but an analytical solution is possible if the force factor is approximated. Indeed, the terms A and B can be found by inserting A cos(wt) + B sin(wt) in the equation ( 12) with an excitation u(t) equals to P sin(wt) where P is an amplitude. The terms C and D can be found by taking the terms with orders higher than one into account, etc...

Theoretical results and position of the small signal parameters according to their created distortions

We present the position of the small signal parameters according to their created distortions. For this purpose, we solve the nonlinear differential equation by using the serial expansion presented in the previous section but we take into account only one variation of a parameter at a time. Figure 9 shows the created distortions by each Thiele and Small parameter. The level of the input voltage is 10V and the frequency of the excitation is 50Hz. We can see that for this electrodynamic loudspeaker, the nonlinear parameter which creates the most important second-harmonic is the inductance Le of the voice coil, and the nonlinear parameter which creates the most important third-harmonic is the equivalent damping parameter R ms . However, when the frequency of the excitation increases, the nonlinear parameters which create more distortions are not the same. Figure 10 presents the created distortions by each Thiele and Small parameter when the frequency of the excitation is 150 Hz. As shown in the previous figure, the level of the input voltage is 10V and the frequency of the excitation is 150 Hz. We see that the nonlinear parameter which creates the most important second-harmonic is the equivalent mass M ms and the nonlinear parameter which creates the more important third-harmonic is the damping parameter R ms . We sum up Table 2: Ranking of the parameters according to their created distortions( the frequency of the excitation is 50 Hz) all the results in tables 2, 3, 4 and 5.

Discussion

The previous section shows that the nonlinear parameters which create the most distortions are not those which are the most nonlinear according to the coil current. For example, figure 7 shows that the equivalent damping parameter R ms is the most nonlinear parameter ( R ms varies a lot with the coil current ) but it is not the parameter which creates the most distortions. On the contrary, the inductance L e of the voice coil varies very little with the coil Moreover, we see that the mechanical compliance C ms is not an important nonlinear parameter since it does not create important distortions. In addition, the ranking of these nonlinearities depend also on input frequency.

Such results are interesting for the design of electrodynamic loudspeakers.

BEHAVIOR OF THE CREATED DISTOR-TIONS ACCORDING TO THE VARIA-TIONS OF THE INPUT VOLTAGE

The previous section presents the ranking of the nonlinear Thiele and Small parameter according to the variation of the frequency excitation. The aim of this section is to show that the distortions created by the nonlinear parameters are more sensitive than the fundamental component according to the variation of the input voltage. To show this, we solve the nonlinear differential equation of the electrodynamic loudspeaker [START_REF] Gander | Moving-coil loudspeaker topology as an indicator of linear excursion capability[END_REF]. We solve this equation by taking four different amplitudes for the excitation ( 1V, 5V, 10V, 50V ) and we represent the spectrum in Fig. 11. The frequency of the excitation

CONCLUSION

The aim of this paper was the study of the spectrum of the electrodynamic loudspeaker. The experimental method, based on the impedance measurement of an electrodynamic loudspeaker allows us to find all the variations of the Thiele and Small parameters. We can say that this experimental method can be used to characterize many transducers which are described with their electrical impedance. Indeed, an electrical impedance can be seen as a nonlinear transfer function which varies with the input current or the input voltage [START_REF] Gille | Systemes asservis non lineaires[END_REF]. In section 2, we have presented a method to derive the coefficients of the nonlinear parameters based on Symplex algorithm. It is noted that a simplifying method has been used in a previous paper [START_REF] Ravaud | Time-varying non linear modeling of electrodynamic loudspeakers[END_REF] and has been improved in this paper.

The small signal coefficients, inserted in the differential equation of the electrodynamic loudspeaker enable us to find the generated harmonics. Many new results are discussed. The equivalent damping parameter is the parameter which is the most nonlinear if we look at its relative variation according to the input current. However, it is not the parameter which creates the most important distortions. This result is important because it gives information about the way of manufacturing an electrodynamic loudspeaker.

More generally, when we take into account the variations of the small signal parameters with the coil current, the mean difference between the experimental and the theoretical values is 0,4Ω , whereas the mean difference is 2,0Ω with Thiele and Small model with constant parameters.

Another interesting result is the weak variation of the electrical inductance with the input current. It is noted that the electrodynamic loudspeaker characterized in this paper is a good one and is less nonlinear than a bad loudspeaker. However, we see that this weak variation creates important distortions. This is why it can be very important to try to build electrodynamic loudspeakers with constant inductance.

Furthermore, we have seen that these generated harmonics become more and more important when the input voltage increases. This result is in fact consistent with all the studies dealing with the modeling of nonlinear systems.

This paper is a first step to derive and class the defects of electrodynamic loudspeakers. The experimental approach taken is certainly the more precise for characterizing the variation of a nonlinear transfer function according to the input level. We can say that such defects are very important to determine because they lower the quality of loudspeakers.

ACKNOWLEDGMENT

We thank James Blondeau for his help with the experimental manipulations.

Bl

  = electrodynamic driving parameter [T.m] R ms = mechanical damping, drag force [N.s.m -1 ] C ms = mechanical compliance of suspension( spider, outer rim)[N -1 .m] L e = inductance of voice coil [H] M ms = equivalent mass of moving voice coil, cone and air [kg] Re = electrical resistance of voice coil [Ω] Rµ = eddy current resistance [Ω] w = angular frequency [rad.s -1 ]These parameters can be represented in terms of an analog circuit in Fig.

1 .Figure 1 :Figure 2 :

 112 Figure 1: Representation of the lumped parameters in terms of an analog circuit

2. 2

 2 Relation between the input current I and the voice coil position X according to the Thiele and Small model This section presents the relation between the input curreent I and the voice coil position X with the linear approximation of the Thiele and Small model.

Figure 3 :

 3 Figure 3: Experimental three-dimensional representation of the modulus of the electrical impedance of the electrodynamic loudspeaker ( x : coil current 0 A 0,2 A)(y : frequency 0 Hz 1000 Hz)(z: impedance 0Ω 30Ω)

Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 4567 Figure 4: The ratio LeLe0 is a function of the coil current ( x : coil current 0A 0,2A) (f=100Hz); Le0 = 0.0017H

Figure 9 :

 9 Figure 9: Theoretical spectrums of Thiele and Small parameters: the frequency of the excitation is 50Hz and the input voltage is 10V

Figure 10 :

 10 Figure 10: Theoretical spectrums of Thiele and Small parameters: the frequency of the excitation is 150Hz and the input voltage is 10V

Figure 11 :

 11 Figure 11: Theoretical spectrum : influence of the input voltage on the created harmonics

Table 1 :

 1 Ranking of the parameters according to their sensitivity to the least square algorithm

	Ranking Parameter Law of variation (100 Hz) D [Ω](100 Hz)
	1	R ms	1.1(1 + 4.09I -8.36I 2 )	1.24
	2	Bl	5.5(1 + 0.33I -1.02I 2 )	1.67
	3	M ms	0.009(1 + 0.56I -0.22I 2 )	1.74
	4	C ms	0.00013(1 + 2.02I -9.3I 2 )	1.86
	5	L e	0.0017(1 -1.68I + 7.58I 2 )	1.98

2

Mms I 2 + ... + µ n Mms I n )

[START_REF] Dobrucki | Nontypical effects in an electrodynamic loudspeaker with a nonhomogeneous magnetic field in the air gap and nonlinear suspension[END_REF]

theoretical values is 0,4Ω whereas the mean difference is 2,0Ω with Thiele and Small model with constant parameters. We present in table 1 the laws of variations of the five parameters that vary according to the coil current, and we give for each parameter the sensitivity to the least square. We propose a ranking of these parameters based on the criterion D. To get this ranking, we proceed as follows : we write that one parameter ( L e , R ms , M ms , Bl or C ms ) is a function of the coil current. We input the polynomial in the theoretical impedance and we use the least square algorithm to determine at the value of the mean difference between the experimental impedance and the theoretical impedance. We obtain in this case the ranking shown in

table 1 .

 1 Moreover, the value of the eddy current resistance is 2,2Ω and the value of the electrical resistance is 3,3Ω. The table2shows the laws of variations of Thiele and Small parameters. The representations of these variations are

	given in Figs 4, 5, 6, 7 and 8.

Table 3 :

 3 Ranking of the parameters according to their created distortions( the frequency of the excitation is 50 Hz)

	Ranking ( 150 Hz) Parameter(harmonic 2) log[x]
	1	M ms	-5.3
	2	L e	-6.4
	3	R ms	-7.7
	4	C ms	-8.2
	5	Bl	-8.4

Table 4 :

 4 Ranking of the parameters according to their created distortions( the frequency of the excitation is 150 Hz)

	Ranking ( 150 Hz) Parameter(harmonic 3) log[x]
	1	R ms	-11.7
	2	L e	-12.6
	3	C ms	-14.2
	4	Bl	-15
	5	M ms	-15.1

Table 5 :

 5 Ranking of the parameters according to their created distortions( the frequency of the excitation is 150 Hz) current but at low frequency this parameter creates important distortions.

Table 6 :

 6 Ranking

	of the parameters according to their created distortions
	(the input frequency is 100Hz)
	is 100 Hz.

L e (i) = L e (1 + µ Le I + µ 2 Le I 2 + ... + µ n Le I n ) (

So the electrical impedance becomes :

We use a least square method to identify all the parameters; this method is based on Symplex algorithm. The Simplex method is a systematic procedure which selects the variable that will produce the largest change towards the minimum solution. This algorithm selects the best choice at each iteration without needing information from previous and future iterations. In our case, the principle of this algorithm is to minimize the difference D between the experimental and the theoretical Nyquist diagrams. Consequently, the two parameters which are minimized are the real part Real(I) and the imaginary part Imag(I) of the electrical impedance defined by eq. [START_REF] Kaizer | Modeling of the nonlinear response of an electrodynamic loudspeaker by a volterra series expansion[END_REF].

When the algorithm converges, we obtain the values of the parameters of the equations ( 5), ( 6), ( 7), ( 8) and ( 9).

Experimental results

When we take into account the variations of the small signal parameters with the coil current, the mean difference between the experimental and the


Vincent Auvray

Uwe Jandt

Raoul Florent

Dirk Schäfer

Improved Vessel Enhancement for Fully Automatic Coronary Modeling

3D coronary modeling extracts the centerlines and width of the coronary arteries from a rotational sequence of angiographies. This process heavily relies on a preliminary filtering of the 2D angiograms that enhances the vessels. We propose an improved vessel enhancement method specifically designed for this application. It keeps the advantages of Hessian-based extraction methods (speed, robustness, multiscale) while bypassing its more important limitations: the blurring of bifurcations, and the incomplete filling of very large vessels. The major contributions of this paper are twofold. First, the classical centered kernel used in Hessian-based methods is substituted with an elongated off-centered kernel. The new filter detects the different orientations involved at a bifurcation: it can answer properly to 'half vessels' beginning at the considered pixel (as opposed to the centered classical filter). The proposed "semi-oriented ridge" filter is also more robust to noise, and it stays multi-scale and quickly computable. Second, an original bifurcation detection and enhancement method is presented, based on the following heuristics: "bifurcations have three vessels (at least) in their immediate neighborhood". More precisely, the semi-oriented ridges answers in each tested orientation θ ∈] -π, π] are stored in a circular histogram. The proposed bifurcation energy is the height of the third peak in this histogram: it will have a significant value at bifurcations only. The performance of the complete framework is demonstrated both on the produced vessel maps and on the final modeling results.

INTRODUCTION

In clinical practice, coronary angiography is mainly limited to the direct analysis of 2D X-ray projections. This can lead to misinterpretations of the actual vessel structure, caused by vessel foreshortening and overlap. The fully automatic 3D modelling of the coronary artery tree proposed in [START_REF] Jandt | Automatic generation of 3D coronary artery centerlines using rotational X-ray angiography[END_REF] tackles this problem, based on X-ray rotational angiography. Automatic modelling heavily depends on its first step: the X-ray projections are initially filtered to enhance vessel-like structures and penalize other objects -yielding a vesselness map (Fig. 1). This paper presents an efficient vesselness map computation method specifically designed for that application. It alleviates the two major drawbacks of classical vesselness maps for 3D coronary modelling: the imperfect large vessel filling and the poor performance at bifurcations. The paper is organized as follows: Section 2 describes the considered automatic 3D modelling method [START_REF] Jandt | Automatic generation of 3D coronary artery centerlines using rotational X-ray angiography[END_REF], and details the limitations of the current vesselness map. Section 3 proposes a novel vessel enhancement method that is free of the main drawbacks of classical techniques. It partly relies on a bifurcation detection and enhancement module that is detailed in Section 4. Section 5 reports vesselness maps and modelling results, and Section 6 contains concluding remarks.

PROBLEM STATEMENT

To make this paper self-contained, we briefly review the considered automatic 3D coronary modelling technique based on rotational X-ray angiography [START_REF] Jandt | Automatic generation of 3D coronary artery centerlines using rotational X-ray angiography[END_REF]. We then focus on the expected vesselness maps properties.

3D coronary arteries modelling

A possible solution to represent the 3D coronary vascular structure is to model it (in its centerlines and vessel width -see Fig. 1.c). Such a 3D coronary modelling operation is fast, robust and automatic [START_REF] Jandt | Automatic generation of 3D coronary artery centerlines using rotational X-ray angiography[END_REF]. It analyses a subset of standard rotational angiography projections belonging to one single cardiac phase. The method is based on a region growing approach, which selects voxels that most probably belong to the vascular structures. More precisely, a measure for the probability R 3D of a voxel x 3D to belong to a vessel is calculated:

R 3D (x 3D ) = N n=1 R 2D n (x 3D ) (1)
with n denoting the current projection and N the number of total projections. R 2D n (x 3D ) refers to the likelihood that the (cone-beam) forward projection x 2D of x 3D over the detector plane (at time n) lies on a coronary artery. R 2D n is computed by previously filtering the projections: this process is the vesselness map computation -or vessel enhancement (Fig. 1.b). The probabilities R 3D (x 3D ) are used to propagate a front that aims at extracting the vessel centerlines. The centerlines representing the most significant vessels are finally selected by a geometry-based weighting criterion.

Vessel Map Properties

Vessel enhancement methods. The whole process therefore relies heavily on the quality of the vessel maps R 2D n . They must enhance vessel-like structures and strongly penalize any other objects. Vessel enhancement has been the object of an extensive work (see [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF] for a review). Three major categories of approaches can be distinguished: matched filters, Hessian operators and model-based detection methods. Which method to choose is very dependant of the final application. In our case, two specific requirements have to be considered. The selected method has to be fast since a typical rotational acquisition contains 120 to 260 projections * . Moreover, the chosen vessel enhancement method needs to be multiscale since the considered angiograms can present a variety of vessels sizes (from thin to very large). This militates for Hessian-based vessel enhancement methods.

Hessian-based vessel enhancement methods properties. Tubular objects result in high second order derivatives D 2 (., θ) of the image I (orthogonally to the dominant orientation). Hessian-based vessel enhancement methods therefore use:

G(x 2D ) = max θ D 2 (x 2D , θ)
(2) * Every one of them must be prefiltered if a modelling is needed in each cardiac phase.

as a vesselness map R(x 2D ) (with θ the orientation). Now, the second order derivative of the image I, in a given direction θ, can be analytically computed [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale detection[END_REF]:

D 2 (x 2D , θ) = -V θ T HV θ (3)
with V θ = (cos θ, sin θ) T and H the Hessian of I. Moreover, it is well known that the maximal value of (3) over every possible orientation can be analytically computed as well (as the Hessian largest eigenvalue). This makes the method particularly fast compared to other vessel enhancement techniques. Furthermore, Hessian-based vessel enhancement can be used in an efficient multi-scale framework [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale detection[END_REF]. The maximum local value (3) is simply computed at different resolutions σ † , yielding different map R σ (x 2D ). The final vessel map then reads:

R(x 2D ) = max σ σ γ R σ (x 2D ) (4)
with γ a parameter typically in the range [0.7, 2]. (in practice, we set γ = 1.5).

In fact, convolving an image with a Gaussian filter and then looking for its maximum second order derivative boils down to looking for the maximum answers of the image to elongated rotating filters (the so-called ridge filters, pictured on Fig. 2) over every possible orientation (defined modulo π). The observation of their impulse response at different resolutions makes clear that they are well-suited to favor elongated objects of different widths. Drawbacks of the Hessian-based vessel enhancement methods. However, every classical vessel enhancement method, including the Hessian-based approach, suffers from the common assumption that there locally exists only one oriented structure [START_REF] Qian | Detection of complex vascular structures using polar neighborhood intensity profile[END_REF]. Consequently, they perform poorly at bifurcations points.

This incorrect behavior at bifurcations has also adverse effects at coarse resolutions, yielding a blurring of the thin vessels that are well processed at fine resolutions. As a result, one cannot correctly enhance thin and very large vessels simultaneously. With conventional settings designed to preserve thin vessels, large vessels are not correctly detected: they appear not to be correctly 'filled'.

Both mentioned effects are illustrated on Fig. 1.d. They impact negatively the 3D modelling: attenuated bifurcations can stop front propagation, possibly even leading to shortcuts. On the other hand, imperfect large vessel fillings yield an erroneous vascular structure. † where σ is the standard deviation of the Gaussian filter applied when computing the Hessian matrix

SEMI-ORIENTED RIDGES

Building on the ridge framework, we propose to construct off-centered, asymmetrical filters able to rotate over a 2π angular amplitude (see Fig. 3). The major advantage of those new filters is that they are likely to detect the different orientations involved at a bifurcation: they can answer properly to 'half vessels' beginning at the considered pixel (as opposed to the centered classical ridges). These new filters can be easily computed: they consist in a weighted addition of shifted versions of the classical ridge, along a given direction θ (in ] -π, π]):

S σ (x 2D , θ) = K k=0 w k R σ (x 2D + kρV θ ) K k=0 w k (5)
However, if used directly this way, semi-oriented ridges would also react in the immediate continuation of elongated objects -even beyond elongated structures. To correct for that effect, the weights w i are associated with adaptive factors: In practice, R σ (x 2D , θ) is computed following Eq.3, with K = 2 and ρ = 2σ. σ GV is set to the image noise level, w ′ 0 to 1, w ′ 1 to 0.7 and w ′ 2 to 0.5. Semi-oriented ridges are used in a multiresolution framework similar to Eq.4.

w k = w ′ k exp - |I(x 2D ) -I(x 2D + kρV θ )| 2 2σ 2 GV ( 6
)
To summarize, semi-oriented ridges allow to keep the advantages of classical ridges: these filters are well-suited to enhance vessels, they can be computed quickly ‡ and they are merged in an efficient multiscale framework. Moreover, they better preserve bifurcations, as will be shown in Section 5.

BIFURCATION ENHANCEMENT

Even if semi-oriented ridges perform noticeably better than classical ridges at bifurcations, they still tend to blur them at coarse resolution. To preserve bifurcations anyway, we propose to detect and enhance them before coarse resolutions are explored. ‡ Even if less quickly than classical ridges: no analytic formula can be derived to explicitly compute the maximum of Sσ(., θ) over every orientations. However, Sσ(., θ) can be computed quickly for a given orientation, using Eq.3 and two rigid translations. The use of the steerable filter theory for the same purpose would require very high order derivatives (given the targeted elongation), and would therefore be much less robust to noise.

Bifurcation Energy

Early works on bifurcation detection [START_REF] Johanson | Representing multiple orientation in 2D with multiple channel histograms[END_REF] exploited local dominant orientations, gathered in an histogram. However, these orientations are defined modulo π and therefore hardly suited for bifurcations. Moreover, it is unclear how the considered local neighborhood is selected, though its size is critical. Finally, these methods also tend to classify curved vessels as bifurcations. Two recent articles proposed to deal with elongated structures including bifurcations without detecting them explicitly. Qian et al proposed to consider a pixel as being part of a vessel if a small circle centered on it crosses at least one vessel [START_REF] Qian | Detection of complex vascular structures using polar neighborhood intensity profile[END_REF]. The authors derived from this idea a vessel detection method able to address any local configuration (e.g. bifurcations), but it could be applied in limited situations only (no background). An elegant framework handling mono-and double-orientation detections was proposed in [START_REF] Aach | Analysis of superimposed oriented patterns[END_REF], yet it proved to lack robustness when applied in angiography. None of these three approaches is multiscale. We propose to explicitly characterize bifurcations, based on the following idea: bifurcations have three vessels (at least) in their immediate neighborhood. To study the vessel density in the neighborhood of each pixel x 2D , we resort to the semi-oriented ridges (Eq.5). Since we are interested in the vessel configuration around x 2D , we modify the weights, setting w 0 to 0, w 1 to 0.5 and w 2 to 1. The filter answer in each tested orientation θ ∈ [-π, π] is stored in a circular histogram (θ-addressed). This process is carried out at different scales: every entry in the histogram corresponds to the best filter answer for the considered orientation over the different scales σ, following Eq.4. Three kinds of circular histograms can be finally expected. If x 2D lies on a linear vessel, the histogram will present two clear peaks. If x 2D lies on a bifurcation, the histogram will present at least three peaks, corresponding to the crossed vessels in the neighborhood. Finally, if x 2D lies in the background, no clear peaks will appear; only background noise will be present in the histogram. We therefore choose the height of the third peak of the circular histogram as an indicator for the bifurcation energy B. It will have a significant value on bifurcations only. Fig. 4 demonstrates the quality of the bifurcation energy B over a typical angiogram.

Bifurcation Enhancement

We extract the bifurcation probability B 2 from the energy B through soft-thresholding. The thresholds are set to ensure that a realistic surface area of the image is labeled as bifurcations -typically 1.5% (see Fig. 4). Our objective now is to enhance the vesselness where bifurcations B 2 have been detected. To do so, we resort to a gray-value-based diffusion of the current semi-oriented ridge image S. It reads at iteration i:

T i (x 2D ) = x 2D 2 ∈V (x 2D ) f (|I(x 2D ) -I(x 2D 2 )|)S i-1 (x 2D ) x 2D 2 ∈V (x 2D ) f (|I(x 2D ) -I(x 2D 2 )|) (7) S i (x 2D ) = (1 -B 2 (x 2D ))S i-1 (x 2D ) + B 2 (x 2D ) max S i-1 (x 2D ), T i (x 2D ) (8)
with S i (.) the vesselness map at time i of the bifurcation enhancement process (S 0 (.) = S(.)), V (x 2D ) a neighborhood of x 2D and f (t) = exp(-t 2 /(2σ GV ) 2 ). Eq.7 refers to a vesselness diffusion, controlled by the image grayvalues. Vesselness is diffused only within regions having similar grayvalues. Eq.8 integrates the diffused vesselness T i into the vesselness map S i where bifurcations have been detected, only if this increases the local vesselness value .

RESULTS

The proposed method involves three steps. At fine resolutions (typically σ = 2 to 4) the vessel maps are computed using the semi-oriented filters § (Section 3). Then, bifurcations are detected and enhanced (Section 4). And finally, coarse resolutions are explored with semi-oriented ridges(σ = 5 to 7). We compare on Fig. 5 the classical and proposed vessel maps on two examples. Bifurcations are better preserved, and vessels better filled. Moreover, a better vessel to background SNR can be noticed. Fig. 6 presents 3D coronary modelling on two cases, based on the classical and the new proposed vessel maps. It shows the extracted centerline reprojected on the vessel maps. The first case shows a better modelling of secondary vessels, and the second example that one shortcut could be avoided. § And the circular histogram values are stored for the subsequent bifurcation detection.

CONCLUSION

We have presented a novel vessel enhancement method, specifically designed for 3D coronary modelling. Its major contributions are twofold. First, a novel multiscale Hessian-based vessel enhancement method was proposed. It substitutes the classical centered kernel with an elongated off-centered asymmetrical kernel, yielding much better behavior at bifurcations. Second, an original bifurcation detection and enhancement method is presented. Bifurcations are detected by extracting the third peak of the circular histogram that stores the local vessel density. This can be an interesting output per se.

The performance of the complete framework was demonstrated by visual inspection. Some clear improvements on the final modelling results were also presented. A perspective for the future is to allow for an even larger simultaneous enhancement of thin and very large vessels.

Figure 1 .

 1 Figure 1. a) original projected angiogram, b) corresponding Hessian-based vessel map, c) targeted modelling (observed from a slightly different perspective) d) details showing the imperfect large vessel filtering and the limited vesselness quality at bifurcation

Figure 2 .

 2 Figure 2. Top: Elongated filter equivalent to the convolution with a Gaussian, followed by a second derivative. Bottom: Elongated filter pictured at different resolutions.

Figure 3 .

 3 Figure 3. Top: Ridge filter (left) compared to the semi-oriented ridge filter (right). Bottom: Computation of the semioriented filter based on shifted ridge filters.

Figure 4 .

 4 Figure 4. Left: bifurcation energy superimposed to the original image, right: top percentiles of the bifurcation energy (red top 0.5%, green 0.5% to 1% and blue 1 to 1.5%).

Figure 5 .

 5 Figure 5. Left: classical ridge vessel map, right: proposed vessel map. Illustration of the better bifurcation enhancement (top), and of the better large vessel filling (bottom).

Figure 6 .

 6 Figure 6. Two modelling examples. The 3D centerlines are reprojected on the vessel maps. Left: classical Hessian-based vessel map, right: proposed vessel map. Top: example showing that thin vessels are better modelled. Bottom: example of an avoided shortcut.


Christophe Paoli
email: christophe.paoli@univ-corse.fr

Cyril Voyant
email: cyril.voyant@univ-corse.fr

Marc Muselli
email: marc.muselli@univ-corse.fr

Marie-Laure Nivet
email: marie-laure.nivet@univ-corse.fr

Solar radiation forecasting using ad-hoc time series preprocessing and neural networks 1

Keywords: Time Series, Preprocessing, Seasonality, Multi-Layer Perceptron

In this paper, we present an application of neural networks in the renewable energy domain. We have developed a methodology for the daily prediction of global solar radiation on a horizontal surface. We use an ad-hoc time series preprocessing and a Multi-Layer Perceptron (MLP) in order to predict solar radiation at daily horizon. First results are promising with nRMSE < 21% and RMSE < 998 Wh/m². Our optimized MLP presents prediction similar to or even better than conventional methods. Moreover we found that our data preprocessing approach can reduce significantly forecasting errors.

Introduction

Artificial intelligence techniques are becoming more and more popular in the renewable energy domain [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF], [START_REF] Mellit | Artificial intelligence techniques for photovoltaic applications: A review[END_REF] and particularly for the prediction of meteorological data such as solar radiation [3], [START_REF] Mubiru | Estimation of monthly average daily global solar irradiation using artificial neural networks[END_REF] [START_REF] Kaligirou | Artificial neural networks in renewable energy systems applications: a review[END_REF] [START_REF] Hocaoglu | Hourly solar forecasting using optimal coefficient 2-D linear filter and feed-forward neural networks[END_REF]. Thereby many research works have shown the ability of Artificial Neural Networks (ANNs) to predict times series of meteorological data. In this study and according to electricity suppliers, we focus on the prediction of global solar irradiation on a horizontal plane for daily horizon. In this way, we have investigated time series forecasting which is a challenge in many fields. Because it has made tremendous progress in the past twenty years in terms of theory, algorithms and applications, we have chosen to study ANNs. Moreover, if we compare to conventional algorithms based on linear models, ANNs offer an attractive alternative by providing nonlinear parametric models. Through the proposed study, we will particularly look at the Multi-Layer Perceptron (MLP) network which has been the most used of ANN architectures in the renewable energy domain [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF], [START_REF] Mellit | Artificial intelligence techniques for photovoltaic applications: A review[END_REF]. The originality of our study is to add an ad-hoc time series preprocessing step before using neural networks. Indeed, as seen in [START_REF] Zhang | Neural network forecasting for seasonal and trend time series[END_REF] a data preprocessing including deseasonalization and detrending can improve ANN forecasting performance. As global solar radiation has a deterministic part, we want to take into account this specificity.

The paper is organized as follow: the section 2 describes the physical phenomena we want to predict and introduces our ad-hoc time series preprocessing. Section 3 presents the neural network architecture we designed. Results are presented and discussed in section 4 where several conventional methods for estimation and modeling of the meteorological data are compared with our methodology. Section 5 concludes and suggests perspectives.

Data analysis and preprocessing

There are two approaches that allow quantifying solar radiation: the "physical modeling" based on physical processes occurring in the atmosphere and influencing solar radiation [START_REF] Badescu | Modelling Solar radiation at the earth surface[END_REF], and the "statistic solar climatology" mainly based on time series analysis [START_REF] Badescu | Modelling Solar radiation at the earth surface[END_REF]. As already said we have chosen to combine these two methods in a gray box approach to improve the quality of prediction. In this work, we have used the physical phenomena in an attempt to overcome the seasonality of the resource.

Meteorological data

Our study proposes to analyze the radiation time series (Wh.m-2) measured at the meteorological station of Ajaccio (METEO FRANCE, Corsica, France, 41°55'N, 8°44'E). The data representing the global solar radiation were measured on a daily basis from January 1971 to December 1989. Thus we have a X t time series to forecast for time t+1; that is at horizon 1. To achieve this, we choose to use a gray box (or semi-physical) model where time series prediction and modeling are mixed. For time series prediction only past values are used to forecast the future values at a given horizon. In the case of modeling, the different physical processes involved are taken into account in order to represent the variable however the horizon is. In the remainder of this paper, we choose the following naming convention: X t designate the time series, and X d,y is the modeling of the variables, where d is the day of the year y. In the next section, an explanation of the physical phenomenon is proposed, and then we describe our time series preprocessing.

Physical phenomenon

We can observe in Figure 1 that the global radiation consists of three types of radiation: direct radiation, diffuse radiation and ground-reflected radiation [START_REF] Liu | Daily insolation on surfaces tilted towards the equator[END_REF]. The ground-reflected radiation does not concern us because we try to predict the radiation on a horizontal surface where the ground reflected radiation does not make sense. For clear sky, ie without cloud cover, global radiation is relatively easy to model because it is primarily due to the distance from the sun sensor [9], [START_REF] Liu | Daily insolation on surfaces tilted towards the equator[END_REF], [START_REF] Hay | Calculation of the solar radiation incident on an inclined surface[END_REF], [START_REF] Perez | Modelling daylight availability and irradiance components rom direct and global irradiance[END_REF], [START_REF] Ineichen | Ground-reflected radiation and albedo[END_REF]. It is not the same, when there are clouds near the detectors. Indeed, these are mostly stochastic phenomena, which depend on the weather site. The spectral analysis of the global radiation series highlights the high periodicity of the phenomenon (almost 365 days period). As proposed in [START_REF] Zhang | Neural network forecasting for seasonal and trend time series[END_REF], it appeared wise to make the time series stationary as much as possible. Deseasonalize and detrend a series allows to eliminate seasonal and trend components without changing the other information. In the present study, we have not considered the solar irradiation like a trend process but only like a seasonal process. The choice of the methodology used depends on the nature of the seasonality. As in our case the seasonality is very pronounced and repetitive, so very deterministic and not stochastic. Moreover, it is possible to physically quantify the components of our irradiation series. Like we have seen in figure 1, the global solar irradiance on horizontal plane depends on direct and diffuse radiation. This specificity allows us to apprehend the periodicity of the phenomenon, and to reduce the non-stationarity of the series.

In the next section, we propose to stationarize radiation data to overcome the deterministic component which is easily quantifiable. Thus we devote to the prediction of cloud cover included in the global radiation of the site.

Ad-hoc time series preprocessing

In fact, by dividing the series by daily extraterrestrial radiation [START_REF] Bourbonnais | Analyse des séries temporelles[END_REF], we can quantify the annual periodicity. It is the first step of our stationarization process. In first approximation, it is possible to quantify the deterministic component of global radiation by the extraterrestrial solar radiation alone (H 0 ). Thus we apply on the original series X d,y (where "d" is the day and "y" is the year) the ratio to trend method. This leads to a new series (S d,y ), known as series of index clarity:

d y d y d t H x S S 0 , , / = = . ( 1
)
After this step a new rigid seasonality is updated, we can lift it with the use of periodic seasonal factors [START_REF] Bourbonnais | Analyse des séries temporelles[END_REF]. This treatment aims to create a new distribution without periodicity. Although this pre treatment tends to stationarize the time series, a test of Fisher shows that seasonality was not optimal. According Bourbonnais [START_REF] Bourbonnais | Analyse des séries temporelles[END_REF], after using a ratio to trend method (H o in this case) to correct rigid seasonalities, we can use a ratio to moving average. This second ratio can be applied when there is no analytical expression of the trend. In our case, we find that H o led a new seasonality which is difficult to model. That's why we made a moving average ratio to overcome the seasonality. In the case of flexible seasonality, ie random in amplitude or period, the filtering techniques by successive moving averages are recommended.

∑ - = + + = m m i y i d y d y d S m S y , , , . 1 . 2 1 /
.

(2)

In our case as 2m+1 = 365 days, we obtain m = 182. To complete the process, then we use the 365 seasonal factors (y d ). These are in fact coefficients which allow to overcome rigid seasonality by a moving average ratio described above. In order to not distort the series, we have considered that the total sum of the components of the series is the same before and after the report (final seasonal factors y d * of the equation 5). The transition coefficients (N = 18, the number of years of history) and the average coefficients of the regular 365 days are given by the equations 3 and 4. A new series seasonally adjusted that represents only the stochastic component of global radiation is given by the equation 6.

        = ∑ = N y y d d y N y 1 , 1 . (3)
In the next section, we present the RNA architecture used.

Neural network architecture and design

The search for the ideal network structure is a complex and crucial task. We have adopted a feed forward Multi-Layer Perceptron (MLP), which is the most commonly used in the renewable energy domain [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF], [START_REF] Mellit | Artificial intelligence techniques for photovoltaic applications: A review[END_REF]. In order to determine the best network configuration, we have tried to study all the parameters available in this network architecture. To perform this optimization, we have considered the practice hypothesis that parameters are orthogonal. We have optimized parameters by considering each other constants. We used the Matlab® software and its neural network toolbox to implement our network. The Matlab® training and testing data sets were set respectively to 80% and 20%. As a result of this iterative process, the selected network (see figure 2) has three neuron layers: input, hidden and output layer. There was no significant difference in the use of 1, 2 and 3 hidden layer architectures. One hidden layer was used in order to minimize the complexity of the proposed ANN model. We tried several input layer configurations: 1, 2, 4, 8 and 15 variables. Best results were obtained with 8 inputs which received as input the endogenous entries S t-1 ,.. , S t-8 normalized on {0,1}. The same process was used to fix the number of hidden layer neurons. We found that 3 neurons were sufficient. Finally, we have one neuron on the output layer t S ˆ.

Concerning the transfer functions the best results were obtained with the Gaussian (hidden layer) and linear (output layer) function. Regarding the training algorithm, many experiences have enabled us to choose the Levenberg-Marquardt optimization (second-order algorithm) with 5000 epochs and a decrease factor (mu) of 0.5. Default values have been use for all other parameters. The early stopping technique was set to the maximum validation failure (the parameter max_fail = 5).

.

St-1 St-2 St-3 St-4 St-5 St-6 St-7 St-8 Bias W 1,1 ,1 W 1, 1, 2 W1 ,1, 3 W1,2,1 b 1, 1 b1,2 Bias W 2 ,1 ,1 W2,2,1 W 2, 3, 1 b2 ,1 1 , 2 )

Results and discussions

Figure 3 summarizes the protocol that has allowed us to conduct our experiences and validate our approach. A first treatment (step 1 of the figure 3) allows to clean the series of atypical points. We have replaced them by the average over the 19 years for the hours and the days corresponding to the problems. Steps 2 and 3 have been described in the previous section and lead to a series corrected. In the step 4 we compare classical forecasting methods outlined in the next section with our optimized MLP. Finally step 5 allows to reverse the preprocessing treatment and to obtain the prediction of global irradiation.

The classical forecasting methods

In order to measure the effectiveness of our approach, we have decided to compare it with the following classical forecasting methods.

The ARIMA techniques [15], [START_REF] Poggi | Forecasting and simulating wind speed in Corsica by using an autoregressive model[END_REF] are reference estimators in the prediction of global radiation field. It is a stochastic process coupling autoregressive component (AR) to a moving average component (MA). After several experiments, we have obtained and decided to use an ARMA [START_REF] Mellit | Artificial intelligence techniques for photovoltaic applications: A review[END_REF][START_REF] Mellit | Artificial intelligence techniques for photovoltaic applications: A review[END_REF]. Bayesian inference [START_REF] Diday | Éléments d'analyse de données[END_REF], [START_REF] Celeux | Analyse discriminante sur variables qualitatives[END_REF] is another classical technique. In this method evidences or observations are used to update or to newly infer the probability that a hypothesis may be true. We have identified that the prediction was better if we had 50 classes and an order of 3. Some authors have tried to use so-called Markov process [19], specifically the Markov chains, which is a stochastic process. The mean idea of this technique is that the description of the present state fully captures all the information that could influence the future evolution of the process. Future states will be reached through a probabilistic process instead of a deterministic one. In our case, we obtained 50 for the dimension of the transition matrix (number of class) and an order of 3 for the chain (determination of the prediction lag). The k-nearest neighbors algorithm (k-NN) [START_REF] Mohammed | Simulating climate change scenarios using an improved Knearest neighbor model[END_REF] is a method for classifying objects based on closest training examples in the feature space. k-NN is a type of instance-based learning, or lazy learning where the function is only approximated locally and all computation is deferred until classification. It can also be used for regression. Unlike previous ones this tool does not use a learning base. The method consists in looking into the history of the series for the case the most resembling to the present case. In our study we choose a k equal to 10. The following section presents the results obtained.

Results

To determine whether our network was really interesting in terms of daily prediction of solar radiation, we compared its performances with the forecasting results obtained with a naive predictor (average over 18 years of the day considered), order 3 Markov chains, order 3 Bayesian inferences, an order 10 k-NN, an order 8 AR without preprocessing, an ARMA(2,2) with preprocessing. The predicted results for each combination were compared statistically using three parameters: the Root Mean Square Error (RMSE), the normalized RMSE (nRMSE), and the Mean Bias Error (MBE):

∑ + - = N i i i M C N RMSE 1 2 ) ( 1 . ( 6
) ∑ ∑ = + - = N i i N i i i M N M C N nRMSE 1 2 1 2 ) ( 1 / ) ( 1 . (7)
∑ + - = N i i i M C N MBE 1 ) ( 1 . (8)
Table 1 presents results we have obtained in the case of an annual error for daily prediction of global solar radiation. We highlight that the predictors other than ARMA and ANN give the same results, slightly higher than those obtained with a naive predictor. Even without preprocessing, ARMA and ANN are the best predictors. The preprocessing improves the quality of prediction and allows access to the 20% error. C i and M i are respectively the i th calculated and measured values and N is the total number of observation. Table 2 details in the MLP case, the annual prediction errors obtained for the years 1988 and 1989. The confidence interval is calculated after 10 simulations. Given the small size of the confidence intervals, we can say that there are very few local minimums. The monthly average error represents the error for the value of the irradiation. We obtain that the combination of the prediction of global radiation received after 1 month is different from an average of 4% of the aggregate measured. The negative MBE means that we underestimate the solar potential on averaged over the year. Since we have an atypical day of low irradiation then there is a tendency to overestimate. The determination coefficient R² is greater than 0,8. Figure 4 shows the errors of prediction and distinguishes the seasons for the years 1988 and 1989. As we can see best results; i.e. less important error, in term of forecast are obtained in summer. These results can be used for example by energy managers who need to avoid using hydraulic power plants in dry season. There should be a compromise between RMSE and nRMSE. The nRMSE are useful for comparison and optimization. But for the interpretation of energy, we must look at the RMSE. The spring season is the most difficult to predict. The absolute error is consistent. However, we find that in summer the error does not exceed 900 Wh / m², while the irradiation is important. MBE are found negative, which indicates an underestimation. The MBE is not significantly different from one season to another. Thus we will always have the same prediction error, whatever the season. Finally, Figure 5 compares the real data of solar radiation with the results obtained with our MLP with preprocessing. The error of prediction is also drawn. The increase in errors at the beginning of the cycle corresponds to the spring when the cloud disturbances are very important. We can see it is very difficult to predict the radiation cause of very noisy variable. Figure 6 shows the correlation between experimental and simulated global irradiation. We systematically overestimate the days when the irradiation is minimal (winter). The points that lie at the very top of the line y = x shows that it is very difficult to predict the days when the irradiation had to be theoretically important. We would undoubtedly have improved the results optimizing an ANN by season, but it would complicate the procedure, and tend to decrease the robustness of the procedure. The next section concludes this paper and suggests prospects.

Conclusion and perspectives

This paper has developed an ANN predictor approach to determine global irradiation at daily horizon in order to help electrical managers. We have used an ad-hoc time series preprocessing and a time series prediction designed MLP. Although the location was very specific, with the proximity to the sea and the mountain that can greatly affect nebulosity, we have obtained relevant results. Seasonal RMSE are less than 998 Wh/m-2 (nRMSE < 21%). ANN processes presents a great interest compared to classical stochastic predictor like ARIMA. Moreover we found that our data preprocessing approach reduce significantly forecasting errors.

The next step of our work will be to validate our predictor on real photovoltaic system. It was recently installed in our laboratory and we are awaiting data. In the future, it seems important to study shorter time horizons. As matter of fact, electrical managers are also interested to horizons that can range from ½ hour to several hours: from 3 hours to 24 hours. Thus others ANN architecture types have to be studied: recurrent ANNs, adaptative ANNs, etc. An ongoing study will be based on the implementation of exogenous variables on the input neurons like METARs data (pressure gradient, temperature, etc.). Determining the relevant variables could be done by the random probe method [21].

Fig. 1 .

 1 Fig. 1. Origin of the three types of radiation: direct radiation, diffuse radiation and groundreflected radiation on a detector at ground level (left of the figure). Modification of the global irradiation profile accordingly to clouds cover (right of the figure).

Fig. 2 .

 2 Fig. 2. Architecture of the optimized MLP. The learning has concerned the years 1971 to 1987 and the performance function was mean square error -MSE.

Fig. 3 .

 3 Fig. 3. Summarize of the protocol followed to obtain the predicted irradiation.

Fig. 4 .

 4 Fig. 4. Seasonal errors for the daily prediction of the years 1988 and 1989 (mean with 95% confidence interval).

Fig. 5 .

 5 Fig. 5. Real data of solar radiation and results obtained with our MLP with preprocessing and error of prediction. Dashed line is the error; solid line is the real data of solar radiation and red points are the prediction.

Fig. 6 .

 6 Fig. 6. Correlation between experimental and simulated global irradiation.

Table 1 .

 1 Annual error for all prediction methods, years 1988 and 1989, 8 simulations

	Prediction methods	nRMSE	Confidence interval
	Naïve predictor	26 %	0%
	Markov Chain (order 3)	25,11 %	0%
	Bayes (order 3)	25,16 %	0%
	k-NN (order 10)	25,20 %	0%
	AR(8) without preprocessing	21,18 %	0,2%
	ARMA(2,2) with preprocessing	20,31 %	<0,1%
	ANN[8,3,1] without preprocessing	20,97 %	0,15%
	ANN[8,3,1] with preprocessing	20,17 %	0,1%

Table 2 .

 2 Annual prediction error for the years 1988 and 1989 with our MLP

		Arithmetic Mean 95% confidence interval
	nRMSE	20,2 %	0,1 %
	RMSE	997,971	6,333
	MBE	-104,239	28,794
	R²	0,801	0,002
	Monthly average error	3,9 %	0,4 %

A part of this research is founded by the Territorial Collectivity of Corsica.


P Le Gauffre
email: pascal.le-gauffre@insa-lyon.fr

F Cherqui

Sewer rehabilitation criteria evaluated by fusion of fuzzy indicators

Keywords: performance assessment, aggregation, decision criteria, uncertainty, fuzzy, robustness

 and is part of a research program (INDIGAU, 2007(INDIGAU, -2010, www.indigau.fr) , www.indigau.fr) supported by the French National Research Agency (ANR). Three main ideas are exposed and illustrated in this paper. 1) Each rehabilitation criterion may be seen as the result of successive aggregations of complementary indicators provided by various methods or sources (visual inspection, O&M data, investigation techniques, etc.) at several spatial scales (sewer segment, catchment). 2) Combining indicators may be computed by means of a limited number of "if ... then ..." rules. The French RERAU 1 method provides 4x4 tables defining how to combine couples of indicators, each one being evaluated on a four-grade ordinal scale (G4 being the highest gravity). 3) The aggregation tables can be converted into fuzzy aggregation rules associated to the use of fuzzy indicators. Compared to crisp data and rules, the use of fuzzy indicators and fuzzy rules provides finer results and improves the robustness of the fusion process.

Introduction: the need to link the condition of assets to the quality of service [START_REF] Baptista | Target definition and assessment of performance in water services, 2 nd Leading Edge Conference on Strategic Asset Management LESAM[END_REF] have pointed out that a major part of strategic asset management (SAM) is the "capability of assessing the infrastructure performance at both asset-specific and system levels", "taking into account the utility objectives, on the one hand, and all the other components of the utility management, on the other hand." D. [START_REF] Parsons | Regulating asset management through serviceability and a common framework for investment planning[END_REF], from OFWAT, underlines that "it is essential to link asset behaviour (and failure modes) to service delivered". One basis of asset management should be the understanding of the link between the condition of the physical assets and the quality of service provided to customers. This need was addressed within the French R&D RERAU program, involving research teams, consulting companies and French utilities such as the Urban Communities of Lille, Lyon, Strasbourg, Caen, etc. The proposed methodology is detailed in a guidebook (Le Gauffre et al., 2004, in French) and is presented in (Le [START_REF] Gauffre | Performance Indicators and Multicriteria Decision Support for Sewer Asset Management[END_REF]. This contribution to sewer asset management may be summarized as follows:

-two sets of similar criteria are defined for supporting the definition of investigation programs and rehabilitation programs;

-investigation and rehabilitation criteria may be calculated from complementary performance indicators (PIs) relating to defects, dysfunctions or impacts;

-PIs may be assessed by observation (by means of visual inspection, investigation techniques, monitoring, etc.), by estimation (ageing models, etc.), or by combination of other PIs;

-PIs refer to individual assets or to subsystems;

-each PI is assessed on a four-grades scale and relationships between PIs are modelled by qualitative models (4x4 tables);

-each proposed criterion assesses a contribution of a particular dysfunction of a sewer segment to a particular impact. Each of the eight defined impacts is linked to some of the ten source dysfunctions (Figure 1).

Dysfunctions: INF: infiltration, EXF: exfiltration, HYD: decrease in hydraulic capacity, SAN: sand silting, BLO: blockage, SPD: destabilization of the ground-pipe system, COR: degradation due to corrosion, ROO: degradation due to roots intrusion, ABR: degradation from abrasion, COL: risk of collapse.

DYSFUNCTIONS INF EXF HYD SAN BLO SPD COR

ROO ABR COL POL X X X X POG X X X X NUH X X X TRA X X X DAB X X OCS X X X OCP X SLC X X X X

IMPACTS

Impacts: POL: pollution of surface waters, POG: pollution of ground and groundwater, NUH: nuisances of a hydraulic nature (flooding, etc.), TRA: disruption of surface activities (traffic, etc.), DAB: damages to the built environment, OCS: network operation surplus costs, OCP: treatment plant operating surplus costs, SLC: cost associated to shortened lifetime. -R/SLC/SPD: shortened lifetime cost (SLC) due to destabilization of the ground-pipe system (SPD).

Another example is detailed in next section.

Linking the asset condition to impacts 1) Each criterion results from the aggregation of complementary indicators

Each rehabilitation criterion may be seen as the result of successive aggregations of complementary indicators provided by various methods or sources (visual inspection, O&M data, investigation techniques, etc.) at several spatial scales (sewer segment, catchment). Figure 2 exposes a criterion defined in the RERAU method (Le [START_REF] Gauffre | Gestion patrimoniale des réseaux d'assainissement urbains. Guide méthodologique[END_REF][START_REF] Gauffre | Performance Indicators and Multicriteria Decision Support for Sewer Asset Management[END_REF]. INF2 is an indicator relating to the factors influencing water infiltration into the sewer segment: water table always or periodically above pipe, and type of backfill (watertight or permeable).

INF4 is derived from visual inspection results [START_REF] Cherqui | CCTV inspection of sewer segments: calibration of performance indicators based on experts' opinions[END_REF][START_REF] Cherqui | Condition grade assessment of sewer segments using crisp and fuzzy calibration procedures[END_REF][START_REF] Ibrahim | Gestion patrimoniale des réseaux d'assainissement urbains : évaluation et agrégation d'indicateurs de performance précis ou flous[END_REF] and is a measure of watertightness deficiency.

INF6 is the result of aggregating INF2 and INF4: it is an indirect measure of a dysfunction (infiltration).

INF3 is a direct measure of dysfunction infiltration that may be available (for instance if infiltration has been observed during the inspection) or not (if no infiltration has been observed in a non favourable context such as a dry period).

INF8 refers to the conclusion that may be obtained - INF1) is an explanation of these spillages [START_REF] Raynaud | A model-based assessment of infiltration and inflow in the scope of controlling separate sanitary overflows at pumping stations[END_REF], sewer segments that contribute to infiltration (INF8) should be listed as rehabilitation needs.

2) Fusion of indicators with "if ... then ..

." rules

Base data are expressed with many different units. For the purpose of combining these data by means of a limited number of "if ... then ..." rules, a common framework is defined and used for converting data into indicators. This framework combines two principles:

-each indicator is assigned a grade out of four possible grades, between G1 (no or low deficiency) and G4 (worst situation);

-indicators are combined by means of aggregation tables.

The RERAU guidebook (Le [START_REF] Gauffre | Gestion patrimoniale des réseaux d'assainissement urbains. Guide méthodologique[END_REF]) provides 4x4 tables defining how to combine couples of indicators (black boxes in Figure 2). [START_REF] Grabisch | Fuzzy aggregation of numerical preferences[END_REF] or [START_REF] Bloch | Fusion: general concepts and characteristics[END_REF] discuss the issue of aggregating complementary information or assessments. Possible aggregation operators may be classified into several categories: conjunctive, disjunctive, compensative, non-compensative. Aggregation tables used in the RERAU PI processing system refer to these categories:

-conjunctive: two indicators relate to two required conditions for an event or a situation (e.g. a sewer segment is said to contribute to a dysfunction if this dysfunction is observed and if the sewer segment is in a poor condition, regarding this dysfunction); the min operator is commonly used as conjunctive operator;

-disjunctive: two indicators relate to two alternative ways of getting an event (cause1 or cause2); the max operator is commonly used as disjunctive operator;

-compensative: refers to a situation where two indicators are combined for assessing a resulting effect (e.g. an impact indicator is evaluated by combining a dysfunction indicator with a vulnerability indicator), the arithmetic mean is a relevant operator for this reasoning ;

-non-compensative: a first case of non compensation is a situation where two indicators refer to two phenomena that may be reinforced by synergy; taking into account this interaction requires a specific synergy operator; a second case of non compensation is a situation where two indicators refer to two parallel sources for assessing a state or a phenomenon: reliability of assessment is reinforced when both sources agree; on the contrary if sources provide different conclusions a composite indicator will retain "all answers from all sources at the price of increased imprecision" [START_REF] Bloch | Fusion: general concepts and characteristics[END_REF].

Figure 3 provides some examples of (crisp) aggregation tables:

-PI1 and PI2 are two indicators combined with a "mean" operator;

-the result (PI3) and PI4 are combined with a "max" operator;

-finally criterion C is evaluated by combining PI5 and PI6 with a "min" operator.

Each cell of a table defines one single rule, for instance: "if PI1 = G2 and PI2 = G4, then PI3 = G3".

PI2 PI4 PI6 G1 G2 G3 G4 G1 G2 G3 G4 G1 G2 G3 G4 G1 1 1 2 2 G1 1 2 3 4 G1 1 1 1 1 PI1 G2 1 2 2 3 G2 2 2 3 4 G2 1 2 2 2 G3 2 2 3 4 G3 3 3 3 4 G3 1 2 3 3 G4 2 3 4 4 G4 4 4 4 4 G4 1 2 3 4 PI3 PI5 C Figure 3
. Three examples of aggregation tables ("1" stands for "G1", "2" for "G2", etc.)

3) Fusion of fuzzy indicators by using fuzzy rules

The RERAU framework allows combining various types of base data, but induces two issues, the first one regarding the conversion of quantitative data into a grade, the second one regarding the aggregation of two grades.

Converting a quantitative base data into a grade may be done by means of four thresholds (or cut-off values). For instance, if we denote 1 S the cut-off value between grades G1 and G2 and x a base data, then a crisp conversion means that ) ( 1

S x ≤

implies that the indicator is assigned grade G1, while ) ( 1

S x >

implies that the indicator is assigned grade G2 at least; thus ) (

1 ε - = S x and ) ( 1 ε + = S
x lead to two different assignments. This common problem may be addressed by using a fuzzy representation: within this approach a given value x may lead to assign grade G1 and grade G2 with two corresponding weights or membership values (e.g. 0.6 and 0.4).

Combining two indicators for assessing a high level indicator (or a criterion) may pose some problems: for instance, averaging grade G2 (assessment of indicator X) and grade G3 (assessment of indicator Y) may lead to two possible crisp aggregation results: G2 or G3; a pessimistic view will give G3, while an optimistic view will give G2. This question may be addressed by using a fuzzy arithmetic mean operator. Averaging G2 and G3 could lead to the following result: indicator Z (derived from X and Y) is assigned a grade which is "between G2 and G3" or "half G2 -half G3"; we will denote "0, 0.5, 0.5, 0" this aggregation result .)

0 ) 4 ( , 5 . 0 ) 3 ( , 5 . 0 ) 2 ( , 0 ) 1 ( ( = = = = G G G G μ μ μ μ
A third issue leads to prefer fuzzy operators to crisp operators. The case where two different sources provide different conclusions (in evaluating the same assumption) is a good example of a need for propagating uncertainty within the PI system: "at the combination level, the choice of an operator is crucial. Conjunctive combination operators are discontinuous in the presence of conflicts and may provide no interesting results at all. Averaging is not realistic since such operators may provide answers that are given by none of the sources, or even rejected individually by each source. If some of the sources are assumed to be reliable, disjunctive operators can be used, which will retain all answers from all sources at the price of increased imprecision" [START_REF] Bloch | Fusion: general concepts and characteristics[END_REF].

The next section is dedicated to the implementation of a fuzzy rule-based information fusion system.

Implementing a fuzzy rule-based information fusion system

Five steps and issues are presented in this section, following [START_REF] Mauris | The aggregation of complementary information via fuzzy sensors[END_REF]:

-representing an indicator with fuzzy linguistic variables; -representing each aggregation table with "if… then…" rules; -combining two premises of a rule; -fuzzy representation of the implication; -aggregating the results of several rules.

1) Representing an indicator with fuzzy linguistic variables

The four terms G1, G2, G3 and G4 may be viewed as four fuzzy subsets allowing a numeric/symbolic conversion. The fuzzy subset theory led to the development of the concept of fuzzy meanings and fuzzy description [START_REF] Mauris | The aggregation of complementary information via fuzzy sensors[END_REF]:

-fuzzy meanings are the representation of fuzzy subsets corresponding to linguistic terms: ) (x Li μ denotes the membership function associated to the linguistic term i L (Figure 4);

the fuzzy description is a simple way of describing a measurement with words: ) ( i x L μ denotes the membership value associated to the linguistic term i L , for a given numeric value (e.g. 12 = x in Figure 4). Each aggregation table can be expressed as a set of if/then rules. We denote:

G1 G2 G3 G4 1 0 G1 G2 G3 G4 x μ Li (x) μ 12 (L i ) 1 0 L i 12
-

k j i R , ,
the rule: "

k j i G Z G Y G X is then is and is if ", with { } G4 G3, G2, G1, , , ∈ k j i G G G - ) , , ( k j i Γ μ
the degree of validity (or weight) of the rule

k j i R , ,
for operator Γ , an aggregation table defined for a given indicator; Figure 5 provides an example of an aggregation table with crisp and fuzzy conclusions:

case ["X is G1" and "Y is G1"] gives a crisp conclusion:

"Z is G1" ( 1 ) 1 , 1 , 1 ( = Γ μ , and 4 , 3 , 2 0 ) , 1 , 1 ( = = Γ k for k μ ); -case ["X is G2" and "Y is G3"] gives a fuzzy conclusion: 3 & 2 5 . 0 ) , 3 , 2 ( = = Γ k for k μ and 4 & 1 0 ) , 3 , 2 ( = = Γ k for k μ . Y G1 G2 G3 G4 X G1 1

3) Combining two premises of a rule

Each aggregation table uses two indicators for defining a high level indicator. Each associated rule has two premises connected by an "and" operator ) is and is

( j i L Y L X .
We denote ) , ( . If we make the assumption that the two indicators are independent, ) , ( [START_REF] Mauris | The aggregation of complementary information via fuzzy sensors[END_REF] and can be derived from

j i L L μ is decomposable
) ( i x L μ and ) ( j y L μ by using an operator of intersection, denoted 1 O (1). ) ( ) ( ) , ( 1 j y i x j i L O L L L μ μ μ = (1)
Several operators of intersection may be used: min, product, etc. Results presented in this paper have been computed with a product operator:

) ( ) ( ) , ( j y i x j i L L L L μ μ μ × =
(2)

4) Fuzzy implication

We denote ) , , (

k j i L L L μ the membership value associated to " is " k L Z which is inferred from " is " i L X
and " is " O (e.g. min operator or product operator) [START_REF] Mauris | The aggregation of complementary information via fuzzy sensors[END_REF]:

) , , ( ) , ( ) , , ( 2 k j i O L L L L L j i k j i Γ = μ μ μ (3) Using (1) for ) , ( j i L L μ , we obtain: ( ) ) , , ( ) ( ) ( ) , , ( 2 1 k j i O L O L L L L j i k j i Γ = μ μ μ μ (4)
Results presented in this paper have been computed with a product operator:

) , , ( ) ( ) ( ) , , ( k j i L L L L L j i k j i Γ × × = μ μ μ μ (5)

5) Aggregating the results of several rules

Conclusion

" is " k L Z
, concerning a high level indicator, may be inferred from several rules

k j i R , , .
We denote ) ( k L μ the membership value associated to the linguistic term k L . This value is calculated by considering all the rules inferring this conclusion, by means of a union operator denoted 3 O :

) , , ( ) ( 4 ... 1 , 4 ... 1 3 k j i j i k L L L O L μ μ = = = (6)
Results presented in this paper have been computed with a sum operator:

∑ = = = ) , , ( ) ( 4 ... 1 , 4 ... 1 k j i j i k L L L L μ μ (7)

Benefits of using a fuzzy rule-based information fusion system

Numerical experiments have been performed in order to assess the benefits of using a fuzzy rule-based information fusion system instead of crisp indicators and criteria. We expose two demonstrative examples:

-example 1 shows that the use of fuzzy indicators and fuzzy aggregation rules provides finer results than crisp indicators;

-example 2 is used to show that fuzzy tools improve the robustness of the fusion process.

Example 1: comparing two sewer segments

This first example is using crisp and fuzzy rules that are presented in Figure 6, regarding criterion R/POL2/SPI/INF "Infiltration contributing to Pollution due to dry-weather spillages" (this criterion and relative indicators were presented in Figure 2). We compare two sewer segments (SS1 and SS2) according to criterion R/POL2/SPI/INF. Grades assigned to base indicators and grades assigned to the criterion are given in Table 1.

For SS1, combining "INF2 is G3" and "INF4 is G4" gives "INF6 is G4" (with crisp and fuzzy rules): major watertightness deficiency combined to water table above pipe leads to a high level of presumption regarding water infiltration into SS1. ) are obtained for criterion R/POL2/SPI/INF (as INF1, V-SPI and SPI have been assigned the highest grades, Table 1).

For sewer segment SS2, crisp aggregation and fuzzy aggregation give "INF8 is G3." Indicator SPI is assigned grade G2 (few events registered during 3-5 years). Receiving waters are considered as sensitive (G2) to very sensitive (G3):

50 . 0 ) 3 ( , 50 . 0 ) 2 ( = = G G μ μ
, whereas a crisp assessment leads to assign grade G3 to indicator V-SPI.

The final result for SS2 is "R/POL2 is G3" within a crisp assessment procedure; whereas a fuzzy assessment process gives 50 This example shows that a fuzzy assessment process provides finer results than those provided by crisp indicators and crisp rules. These two sewer segments are assigned the same grade (G3) with a crisp method, but can be clearly distinguished with a fuzzy assessment method: SS1 appears as a priority projects whereas the condition of SS2 appears less severe.

. 0 ) 3 ( , 50 . 0 ) 2 ( = = G G μ μ .
Example 2: improving robustness of the performance assessment Figure 7 illustrates a situation where visual inspection may lead to two different condition assessments: let us denote 1

x (respectively 2 x ) the value of a numerical score x (Le [START_REF] Gauffre | Performance Indicators and Multicriteria Decision Support for Sewer Asset Management[END_REF] obtained by interpretation of the defects reported by operator 1 (respectively operator 2). Dirksen et al. (2007) have shown that results obtained from 2 different human operators may be highly different.

A crisp assessment of indicator INF4 can be obtained by comparing 1 x or 2 x with cut-off values. Here we consider the two following scenarios:

-scenario 1:

1 x x = is close to the cut-off value 2
S between grades G2 and G3, and x x ≅ ) will provide similar conclusions.

The same applies when considering uncertainty linked to the definition of a crisp cut-off value.

Defining a crisp cut-off value is a calibration process [START_REF] Ibrahim | Sewer asset management: from visual inspection survey to dysfunction indicators, 2 nd Leading Edge Conference on Strategic Asset Management LESAM[END_REF][START_REF] Cherqui | Condition grade assessment of sewer segments using crisp and fuzzy calibration procedures[END_REF]) that is using available data (more or less reliable) together with parameters and hypotheses that are more or less precise (importance given to assignment errors such as a false negative or a false positive, etc.) For instance, ' 2 S may be an alternative value to 2 S . Figure 7 G G μ μ will slightly vary with small variations of the parameters and hypotheses that are used in the calibration process (Cherqui et al. 2009). Results in Table 2 show that a crisp assessment framework leads to switch from grade G2 to grade G3 (or the opposite) for criterion R/POL2/SPI/INF when a slight change occurs in reporting defects or in calibrating a crisp cut-off value for indicator INF4.

) ( 2 x G μ ) ( 3 x G μ x x 1 x 2 S 2 S' 2
On the contrary, a fuzzy framework allows dealing with uncertainties linked to data and to rules. This leads to smooth the assessment process and provides more robust results.

As robustness is a major criterion for disseminating models and tools towards utilities and consulting companies, this example illustrates, together with other works such as (Bagheri et al. 2006; Rajani et al. 2006; Sadiq et al. 2004), the benefits of designing performance indices and decision support tools based on a fuzzy inference approach.

Summary and conclusion

One basis of asset management should be the understanding of the link between the condition of the physical assets and the quality of service provided to customers. Within this concern, three main ideas are exposed and illustrated in this paper:

-a rehabilitation criterion may be seen as the result of successive aggregations of complementary indicators provided by various methods or sources (visual inspection, O&M data, investigation techniques, etc.) at several spatial scales (sewer segment, catchment);

-combining indicators may be computed by means of a limited number of "if .. then ..." rules; the French RERAU method provides 4x4 tables defining how to combine couples of indicators, each one being evaluated on a four-grade ordinal scale (G4 being the highest gravity);

-to deal with uncertainty linked to base data and to knowledge or decision rules, the crisp aggregation tables can be converted into fuzzy aggregation rules associated to the use of fuzzy indicators.

Fuzzy indicators allow dealing with the uncertainty linked to base data (for instance a score calculated after visual inspection, based on reported defects) or linked to PI calibration processes (for instance the uncertainty linked to the hypotheses that are used for defining a cut-off value).

Fuzzy rules allow dealing with imprecise human reasoning when combining complementary indicators into a more synthetic indicator. Additionally, fuzzy rules (rules with several conclusions) allow to compute the uncertainty arising from the use of two complementary sources of information on a given phenomenon (for instance INF3 and INF6 may reinforce themselves if they provide same results, this will lead to a low level of fuzziness for composite indicator INF8; on the contrary two different results will lead to a high level of fuzziness for INF8).

Examples provided in this paper tend to demonstrate that using a fuzzy framework, instead of crisp indicators and crisp rules, improves the assessment process by providing finer and more robust conclusions:

example 1 shows that the use of fuzzy indicators and fuzzy aggregation rules provides finer results than crisp indicators;

example 2 is used to show that fuzzy tools improve the robustness of the fusion process.

Figure 1 .

 1 Figure 1. Decision criteria defined by linking impacts and source dysfunctions.

Figure 2 .

 2 Figure 2. Indicators related to criterion R/POL2/SPI/INF "Infiltration contributing to Pollution due to dryweather spillages"

Figure 4 .

 4 Figure 4. Link between fuzzy meanings and fuzzy description

Figure 5 .

 5 Figure 5. An aggregation table with fuzzy conclusions: four values

  table Γ , by means of an intersection operator, denoted 2

Figure 6 .

 6 Figure 6. Fuzzy rules used for assessing criterion R/POL2/SPI/INF, with corresponding crisp rules represented by shaded cells, and values of indicators (in circles) for sewer segment SS1.

  Additionally, direct observation provides "INF3 is G3". INF3 and INF6 are two possible ways or sources for assessing water infiltration into pipe segment SS1. In this case these two indicators provide similar results, without conflicts. Crisp aggregation gives "INF8 is G3" whereas fuzzy aggregation gives

Figure 7 .

 7 Figure 7. Two possible values 1 x & 2 x of a score x computed after visual inspection of a sewer segment, are translated into grades by comparison to a crisp cut-off value ( 2 S or ' 2 S ) or by using

Table 1 .

 1 Base indicators and criterion R/POL2/SPI/INF evaluated for two sewer segments

	Sewer segments

RERAU: « Réhabilitation des Réseaux d'Assainissement Urbains » (rehabilitation of urban sewer networks).

IWA & AWWA, Miami, Florida, November 11-13, 2009

Acknowledgements

Works presented in this paper were initiated in the frame of the French R&D RERAU program ("Réhabilitation des Réseaux d'Assainissement Urbains") and are developed in the frame of the INDIGAU research project ("indicateurs de performance pour la gestion patrimoniale des réseaux d'assainissement urbains"; www.indigau.fr) supported by the French National Agency for Research -ANR, within the PGCU-2006 program (project # ANR-06-RGCU-003). The authors would like to thank all the participants of the RERAU and INDIGAU projects.

Aggregating results of these 3 rules provides the assessment of INF8 :