Vitaliy Feoktistov
email: vitaliy.feoktistov@ema.fr

Stefan Janaqi
email: stefan.janaqi@ema.fr

NEW ENERGETIC SELECTION PRINCIPLE IN DIFFERENTIAL EVOLUTION

Keywords: differential evolution, evolutionary algorithms, heuristics, optimization, selection

The Differential Evolution (DE) algorithm goes back to the class of Evolutionary Algorithms and inherits its philosophy and concept. Possessing only three control parameters (size of population, differentiation and recombination constants) DE has promising characteristics of robustness and convergence. In this paper we introduce a new principle of Energetic Selection. It consists in both decreasing the population size and the computation efforts according to an energetic barrier function which depends on the number of generation. The value of this function acts as an energetic filter, through which can pass only individuals with lower fitness. Furthermore, this approach allows us to initialize a population of a sufficient (large) size. This method leads us to an improvement of algorithm convergence.

INTRODUCTION

Evolutionary Algorithms increasingly become the primary method of choice for optimization problems that are too complex to be solved by deterministic techniques. They are universal, robust, easy to use and inherently parallel. The huge number of applications and continuous interest prove it during several decades [START_REF] Heitkötter | Hitch Hiker's Guide to Evolutionary Computation : A List of Frequently Asked Questions (FAQ)[END_REF][START_REF] Beasley | Possible applications of evolutionary computation[END_REF]. In comparison with the deterministic methods Evolutionary Algorithm require superficial knowledge about the problem being solved. Generally, it is only necessary that the algorithm is able to evaluate the objective function for a given set of input parameters. Nevertheless, in most cases such heuristics take less time to find the optimum than, for example, gradient methods. One of the latest breakthroughs in the evolutionary computation is the Differential Evolution algorithm.

DIFFERENTIAL EVOLUTION

Differential Evolution is a recently invented global optimization technique [START_REF] Storn | Differential evolution -a simple and efficient adaptive scheme for global optimization over continuous spaces[END_REF]. It can be classified as an iterative stochastic method. Enlarging the Evolutionary Algorithms' group, DE turns out to be one of the best population-based optimizers [START_REF] Storn | Minimizing the real functions of the ICEC'96 contest by differential evolution[END_REF]. In the following lines we give a brief description of DE algorithm.

An optimization problem is represented by a set of variables. Let these variables form a D-dimensional vector in continuous space X = (x 1 , . . . , x D ) ∈ IR D . Let there be some criterion of optimization f (X) : IR D → IR, usually named either fitness or cost function. Then the goal of optimization is to find the values of the variables that minimize the criterion, i.e. to find

X * : f (X * ) = min X f (X) (1)
Often, the variables satisfy boundary constraints

L ≤ X ≤ H : L, H ∈ IR D (2)
As all Evolutionary Algorithms, DE deals with a population of solutions. The population IP of a generation g has N P vectors, so-called individuals of population. Each such individual represents a potential optimal solution.

IP g = X g i , i = 1, . . . , N P (3)
In turn, the individual contains D variables, so called genes.

X g i = x g i,j , j = 1, . . . , D (4)
The population is initialized by randomly generating individuals within the boundary constraints,

IP 0 = x 0 i,j = rand i,j • (h j -l j ) + l j (5)
where rand function generates values uniformly in the interval [0, 1].

Then, for each generation the individuals of a population are updated by means of a reproduction scheme. Thereto for each individual ind a set of other individuals π is randomly extracted from a population. To produce a new one the operations of Differentiation and Recombination are applied one after another. Next, the Selection is used to choose the best. Now briefly consider these operations.

Here, we show a typical model of the Differentiation, others can be found in [START_REF] Storn | Differential evolution -a simple and efficient adaptive scheme for global optimization over continuous spaces[END_REF]. For that, three different individuals π = {ξ 1 , ξ 2 , ξ 3 } are randomly extracted from a population. So, the result, a trial individual, is

τ = ξ 3 + F • (ξ 2 -ξ 1 ) , (6)
where F > 0 is the constant of differentiation.

After, the trial individual τ is recombined with updated one ind. The Recombination represents a typical case of a genes' exchange. The trial one inherits genes with some probability. Thus,

ω j = τ j if rand j < Cr ind j otherwise (7)
where j = 1, . . . , D and Cr ∈ [0, 1) is the constant of recombination.

The Selection is realized by comparing the cost function values of updated and trial individuals. If the trial individual better minimizes the cost function, then it replaces the updated one.

ind = ω if f (ω) ≤ f (ind) ind otherwise (8)
Notice that there are only three control parameters in this algorithm. These are N P -population size, F and Cr -constants of differentiation and recombination accordingly. As for the terminal conditions, one can either fix the number of generations g max or a desirable precision of a solution V T R (value to reach).

The pattern of DE algorithm is presented in the following way (see algorithm 1).

ENERGETIC APPROACH

We introduce a new energetic approach which can be applied to population-based optimization algorithms including DE. This approach may be associated with the processes taking place in physics.

Let there exist a population IP consisting of N P individuals. Let us define the potential of individual as its cost function value ϕ = f (ind). Such potential shows the remoteness from the optimal solution ϕ * = f (ind * ), i.e. some energetic distance (potential) that should be overcome to reach the optimum.

Algorithm 1 Differential Evolution

Require: F, Cr, N P -control parameters initialize IP 0 ← {ind 1 , . . . , ind N P } evaluate f (IP 0 ) while (terminal condition) do for all ind ∈ IP g do

IP g → π = {ξ 1 , ξ 2 , . . . , ξ n } τ ← Dif f erentiate(π, F ) ω ← Recombine(τ, Cr) ind ← Select(ω, ind) end for g ← g + 1 end while
Then, the population can be characterized by superior and inferior potentials ϕ max = max f (ind i ) and ϕ min = min f (ind i ). As the population evaluates the individuals take more optimal energetic positions, the closest possible to the optimum level. So if t → ∞ then ϕ max (t) → ϕ min (t) → ϕ * , where t is an elementary step of evolution. Approaching the optimum, apart from stagnation cases [START_REF] Lampinen | On stagnation of the differential evolution algorithm[END_REF], can be as well expressed by ϕ max → ϕ min or (ϕ maxϕ min ) → 0. By introducing the potential difference of population △ϕ(t) = ϕ max (t)-ϕ min (t) the condition of optimality is represented as

△ϕ(t) → 0 (9)
In other words, the optimum is achieved when the potential difference is closed to 0 or to some desired precision ε. The value △ϕ(t) is proportional to the algorithmic efforts, which are necessary to find the optimal solution. Thus, the action A done by the algorithm in order to pass from one state t 1 to another t 2 is

A(t 1 , t 2 ) = t2 t1 △ϕ(t)dt (10)
We introduce then the potential energy of population E p that describes total computational expenses.

E p = ∞ 0 △ϕ(t)dt (11)
Notice that the equation ( 11) graphically represents the area S p between two functions ϕ max (t) and ϕ min (t).

Let us remind that our purpose is to increase the algorithm convergence. Logically, the convergence is proportional to computational efforts. It is obvious the less is potential energy E p the less computational efforts are needed. Thus, by decreasing the potential energy E p ≡ S p we augment the convergence rate of the algorithm. Hence, the convergence increasing aim transforms into a problem of potential energy minimization (or S p minimization).

E * p = min △ϕ(t) E p (△ϕ(t)) (12)

NEW ENERGETIC SELECTION PRINCIPLE

Idea

We apply the above introduced Energetic Approach to the DE algorithm. As an elementary evolution step t we choose a generation g.

In order to increase the convergence rate we minimize the potential energy of population E p (fig. 1). For that a supplementary procedure is introduced at the end of each generation g. The main idea is to replace the superior potential ϕ max (g) by so called energetic barrier function β(g). Such function artificially plays down the potential difference of generation △ϕ(g).

β(g) ≤ ϕ max (g), ∀g ∈ [1, g max ] (13)
From an algorithmic point of view this function β(g) serves as an energetic filter for the individuals passing into the next generation. Thus, only the individuals with potentials less than the current energetic barrier value can participate in the next evolutionary cycle (fig. 2).

Practically, it leads to the decrease of the population size N P by rejecting individuals which violate the condition (14).

f (ind) ≤ β(g) (14)

Energetic Barriers

Here, we show some examples of the energetic barrier function. At the beginning we outline the variables which this function should depend on. Firstly, this is the generation variable g, which provides a passage from one evolutionary cycle to the other. Secondly, it should be the superior potential ϕ max (g) that presents the upper bound of the barrier function. And thirdly, it should be the inferior potential ϕ min (g) ensuring the low bound of the barrier function (fig. 3). Linear energetic barriers. The simplest example is the use of a proportional function. It is easy to obtain by multiplying either ϕ min (g) or ϕ max (g) on a constant K.

In the first case, the value ϕ min (g) is always stored in the program as the current best value of the cost function. So, the energetic barrier looks like

β 1 (g) = K • ϕ min (g), K > 1 (15)
The constant K is selected to satisfy the energetic barrier condition (13).

In the second case, firstly a little procedure is necessary to find superior potential (maximal cost function value of the population) ϕ max (g). Here, the energetic barrier is

β 2 (g) = K • ϕ max (g), K < 1 (16)
K should not be too small in order to provide a smooth decrease of the population size N P .

An advanced example would be a superposition of the potentials.

β 3 (g) = K • ϕ min (g) + (1 -K) • ϕ max (g) (17)
So, with K < 1 the energetic barrier function is always found between the potential functions. Now, by adjusting K it is easier to get the smoothed reduction of the population without condition violation ( 13).

Examples of the energetic barrier functions are shown on the figure (fig. 4).

Nonlinear energetic barriers.

As we can see the main difficulty of using the linear barriers appears when we try to define correctly the barrier function in order to provide a desired dynamics of the population reduction. Taking into consideration that ϕ max → ϕ min when the algorithm converges, the ideal choice for the barrier function is a function which begins at a certain value between ϕ min (0) and ϕ max (0) and converges to ϕ max (g max ).

Thereto, we propose an exponential function K(g)

K(g) = K l + (K h -K l ) • e (-T gmax •g) (18)
This function smoothly converges from K h to K l . The constant T , so called temperature, controls the convergence rate. The functional dependence on the temperature constant K(T ) is represented on the figure (fig. 5). By substituting the constant K in the equations (15-17) for the exponential function (18) we can supply the energetic barrier function with improved turning (fig. 6).

Advantages

Firstly, such principle of energetic selection permits to initialize the population of a sufficiently large size. This fact leads to better (careful) exploration of a search space on the initial generations as well as it augments the probability of finding the global optimum.

Secondly, the introduction of the energetic barrier function decreases the potential energy of the population and thereby increases the algorithm convergence.

Thirdly, a double selection principle is applied. The first one is a usual DE selection for each individual of a population. Here, there is no reduction of the population size. And the second one is a selection of the best individuals which pass in the next generation, according to the energetic barrier function. It leads to the reduction of the population size.

Remark. Notice that a considerable reduction of the population size occurs at the beginning of the evolutionary process. For more efficient exploitation of this fact a population should be initialized with greatly larger size N P 0 than usually. Then, when the population shrinks to a certain size N P f , it is necessary to stop the energetic selection procedure. This forced stopping is explained by possible stagnation and not enough efficient search in a small size population [START_REF] Lampinen | On stagnation of the differential evolution algorithm[END_REF][START_REF] Price | New Ideas in Optimization, Part 2 : Differential Evolution[END_REF].

COMPARISON OF RESULTS

In order to test our approach we chose three test functions (19) from a standard test suite for Evolutionary Algorithms [START_REF] Whitley | Evaluating evolutionary algorithms[END_REF]. The first two functions, Sphere f 1 and Rosenbrock's function f 2 , are classical De Jong testbads [START_REF] Jong | An analysis of the behavior of a class of genetic adaptive systems[END_REF]. Sphere is a "dream" of every optimization algorithm. It is smooth, unimodal and symmetric function. The performance on the Sphere function is a measure of the general efficiency of the algorithm. Whereas the Rosenbrock's function is a nightmare. It has a very narrow ridge. The tip of the ridge is very sharp and it runs around a parabola. The third function, Rotated Ellipsoid f 3 , is a true quadratic non separable optimization problem.

f 1 (X) = 3 i=1 x 2 i f 2 (X) = 100(x 2 1 -x 2 ) 2 + (1 -x 1 ) 2 f 3 (X) = 20 i=1   i j=1 x j   2 (19)
We fixed the differentiation F and recombination Cr constants the same for all functions. F = 0.5. Recombination Cr = 0 (there is no recombination) in order to make the DE algorithm rotationally invariant [START_REF] Salomon | Re-evaluating genetic algorithm performance under coordinate rotation of benchmark functions : A survey of some theoretical and practical aspects of genetic algorithms[END_REF][START_REF] Price | New Ideas in Optimization, Part 2 : Differential Evolution[END_REF]. The terminal condition of algorithm is a desirable precision of optimal solution V T R (value to reach). It is fixed for all tests as V T R = 1.0e -6. We count the number of function evaluations N F E needed to reach the V T R. The initial data are shown in the table 1. For DE with energetic selection principle the initial population size was chosen three times larger than in the classical DE scheme : N P 0 = 3 • N P . The forced stopping was applied if the current population became smaller than N P . Hence N P f ≤ N P . As an energetic barrier function the linear barrier β 3 (g) was selected (17). So, K is an adjusting parameter for barrier turning. D is the dimension of test functions. The numbers of function evaluations (N F E's) were compared. It is considered that N F E cl value is equal to 100% therefore the relative convergence amelioration in percentage wise can be defined as

δ = 1 - N F E es N F E cl (20)
Thus, δ may be interpreted as the algorithm improvement.

CONCLUSION

The variation of the population size of populationbased search procedures presents a rather promising trend. In this article we have examined its decrease. The proposed energetic approach explains a theoretical aspect of such population reduction. The efficiency of the new energetic selection principle based on this energetic approach is illustrated by the example of the DE algorithm. The given innovation provides more careful exploration of a search space and leads to the convergence improvement. Thus, the probability of the global optimum finding is increased. Further works are carried on the methods of increasing the population size.

  FIG. 1 -Energetic approach.

  FIG. 2 -Energetic filter.

  FIG. 3 -Energetic barrier function.

  FIG.4-Linear energetic barriers.

FIG. 5 -

 5 FIG. 5 -Exponential function K(g, T ).

  FIG. 6 -Nonlinear energetic barrier.

  The average results of 10 runs for both the classical DE scheme and DE with energetic selection principle are summarized in the table 2. Comparison the classical DE scheme (cl) and DE with energetic selection principle (es).

	f i	N F E cl N F E es	δ, %
	1	1088.7	912.4 16,19
	2	1072.9	915.3 14,69
	3 106459.8 94955.6 10,81
	TAB. 2 -


UMRAE, Cerema Jean-Rémy Gloaguen
email: jean-remy.gloaguen@cerema.fr

David Ecotière

Benoit Gauvreau

Estimation of the sound emergence of wind turbines by semisupervised learning technique

The sound emergence is the main regulatory estimator for wind turbine noise in France. This criterion aims to limit their noise impact on local residents and is highly dependent on the variation of residual noise over time. Therefore, initially defined curtailment plans can sometimes become inadequate, in which case they cannot easily be updated without leading to significant production losses. Machine learning techniques allow today to consider the continuous estimation by measurements of the sound contribution of wind turbine noise in the ambient noise and thus its noise emergence, without needing to stop the wind farm. This operation makes it possible not only to regularly adapt these reduction plans, thus optimizing electricity production, but also limiting the possible noise annoyance for local residents. For this purpose, semi-supervised Non-negative Matrix Factorization method is considered, enhanced by a temporal regularity constraint. This approach combines a wind turbine dictionary designed on a learning basis and a free dictionary that allows the adaptation of the method to the variability of residual noise. Tests conducted on simulated measurements reveal satisfactory performances with mean estimation errors lower than 2 dBA for wind noise emergences lower than 5 dBA. Finally, the presence of these two types of dictionaries makes it possible to estimate the wind noise emergence according to one or the other depending on the predominance of the estimated wind turbine noise.

Introduction

In some European countries, the sound impact of wind farms is assessed by their sound emergences. This indicator is estimated by the difference between the A-statistical sound pressure level (SPL) of the ambient noise at the receiving point 𝐿 𝐴,50,𝑎𝑚𝑏. (i.e. when the wind farm is operating) and the residual noise 𝐿 𝐴,50,𝑟𝑒𝑠. (i.e. when the wind turbine is shutdown): 𝐸 = 𝐿 𝐴,50,𝑎𝑚𝑏. -𝐿 𝐴,50,𝑟𝑒𝑠.

(1)

The statistical SPL presents the advantage to be more suited for the wind turbine (WT) noise, which remains quite constant during time, while equivalent SPL is more sensitive to extreme sound pressure levels that might result from emergences of residual noises. These sound emergences are then limited to regulation thresholds to not exceed. To respect these criteria, the most common solution is the application of a curtailment plan on WTs to limit their operation and thus their noise emission [START_REF] Rogers | Optimal strategies for wind turbine environmental curtailment[END_REF]. However, this plan affects strongly the production capacity of the machines. Furthermore, this solution has the main drawback to depend on the variability of the residual noise, which is a time varying component that changes over the meteorological conditions, the ground characteristics and the surrounding infrastructures. Consequently, this may result in two situations where the curtailment plan might be too strong (i.e. the electrical production is unnecessarily limited) or too low (i.e. the sound emergences exceed the regulation thresholds).

Currently, to update this plan, on/off cycle measurements are made. These measurements consist in assessing alternatively the SPLs of the ambient noise and the residual noise by stopping the wind farm periodically. If needed, a new curtailment plan can then be formulated.

However, these measurements require the shutdown of the machines, which also significantly affects the electrical production. Furthermore, they are only carried out during a couple of weeks and are not sufficiently representative over a long period. Consequently, this updated plan does not guarantee to be still valid few months later. To repeat this process regularly will then increase again the loss of electrical production and is not a suitable solution.

Consequently, there is an interest to be able to estimate the sound emergence of wind farms continuously without stopping them. Thus, it would be possible to control regularly the noise annoyance for neighboring inhabitants, while optimizing the electrical production. Instead of seeking to improve the sound source emission of WT and the models for outdoor sound propagation [START_REF] Cotté | Extended source models for wind turbine noise propagation[END_REF], this emergence could be estimated by measurements made in situ. Recently, [START_REF] Gloaguen | Automatic estimation of the sound emergence of wind turbines using Non-negative Matrix Factorization : A preliminary study[END_REF] propose a first tool to estimate the sound emergence of WTs. This tool is based on the Non-negative Matrix Factorization (NMF) used as a source separation method to first extract the WT component from simulated ambient scenes and then to estimate the sound emergences. This preliminary study has been extended in [START_REF] Gloaguen | Automatic estimation of the sound emergence of wind turbine noise with Non-negative Matrix Factorization[END_REF] on a larger number of situations where the WT component has been propagated at different distances in different propagation conditions. Supervised learning associated with a temporal constraint reached satisfactory results with a mean error inferior to 2 dBA on many cases. In this study, only WT noise was considered in NMF. It might be then helpful to better consider the residual noise in NMF to bring more flexibility but also robustness to this method. Thus, this paper extends these works and considers a semi-supervised learning, which allows NMF to consider labelled data (i.e. WT noise) and unknown data that can be adapted to the diversity of the residual noise. Section 2 introduces Semi-supervised NMF, Section 3 presents the different corpora used in this experiment and Section 4 summarizes the main results and details the behavior of the proposed method.

Semi-supervised Non-negative Matrix Factorization

In audio field, Non-negative Matrix Factorization (NMF) is usually introduced as a linear approximation of an audio spectrogram 𝑉 𝐹×𝑁 by the product of two matrices, 𝑊 and 𝐻, such as 𝑉 ≈ 𝑉 ̃= 𝑊𝐻. 𝑊 𝐹×𝐾 , called dictionary, includes audio spectra and 𝐻 𝐾×𝑁 , called the activation matrix, corresponds to the temporal evolution of each of the spectra such as 𝑉 ̂= 𝑊𝐻 [START_REF] Dikmen | Sound event detection using non-negative dictionaries learned from annotated overlapping events[END_REF][START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. As only amplitude (or power) spectra are considered in 𝑉 and 𝑊, only additive combinations are possible. The dimensions of the matrices are chosen such as 𝐾 < min (𝐹, 𝑁) to avoid overcomplete representation. In supervised learning, 𝑊 is learned on labelled samples and 𝐻 is the unknown to estimate. In the present case, 𝑊 includes WT noise spectra.

For semi-supervised learning [START_REF] Kitamura | Music Signal Separation Based on Supervised Nonnegative Matrix Factorization with Orthogonality and Maximum-Divergence Penalties[END_REF], a second dictionary and activation matrix are considered, respectively 𝑌 𝐹×𝐽 and 𝑍 𝐽×𝑁 . The dimension 𝐽 corresponds to the number of considered spectra, which is often set such as 𝐽 < 𝐾 in order to allow NMF to stay focused on the data learned in 𝑊.

These two matrices are learned on each spectrogram to represent, the best as possible, the residual noise. This second part makes it possible to increase the degree of freedom of NMF and to include residual noise directly learn on the spectrogram 𝑉 ̃. In supervised learning, it would have been necessary to learn a residual dictionary on labelled data, which is more complex to do, as it is a noise composed of multiple different kind of sound sources not easy to define precisely. Furthermore, instead of a linear sum between [𝑊𝐻] and [𝑌𝑍], the energetic sum between the two components is considered here:

𝑉 ≈ 𝑉 ̃= √[𝑊𝐻] 2 + [𝑌𝑍] 2 .
(2) (𝑥 -𝑦) 2 . The added regularization term, weighted by the coefficient 𝛼 𝑠𝑚 , corresponds to a temporal regularity constraint applied only on 𝐻 [START_REF] Févotte | Single-Channel Audio Source Separation with NMF : Divergences, Constraints and Algorithms[END_REF]. This constraint forces the shape of the activations in 𝐻 to adapt a smoother behaviour and thus to better corresponds to the time signature of WT noise, which is quite constant in practice. As the NMF formulation has been modified on Eq. 2, the update algorithms of 𝐻, 𝑌 and 𝑍 are changed.

The new update algorithms are estimated as in [START_REF] Lee | Algorithms for Non-negative Matrix Factorization[END_REF] with a descent gradient. Equation 4 presents the update algorithms for the Euclidean distance for this new formulation:

𝑌 → 𝑌. [[𝑌𝑍]. 𝑉] 𝑍 𝑇 [[𝑌𝑍]. 𝑉 ̃] 𝑍 𝑇 (4.1) 𝑍 → 𝑍. 𝑌 𝑇 [[𝑌𝑍]. 𝑉] 𝑌 𝑇 [[𝑌𝑍]. 𝑉 ̃] (4.2) 𝐻 → 𝐻. 𝑊 𝑇 [[𝑊𝐻]. 𝑉] + 2𝛼[𝐻 +1 + 𝐻 -1 ] 𝑊 𝑇 [[𝑊𝐻]. 𝑉 ̃] + 2𝛼 [𝐻 + 𝐻 12 ] (4.3) with 𝐻 -1 = [𝟎 𝐻 1:𝑁-1 ], 𝐻 +1 = [𝐻 2:𝑁 𝟎] and 𝐻 12 = [𝟎 𝐻 2:𝑁-1 𝟎]
where 𝐻 𝑛:𝑚 refers to the selection of columns 𝑛 to 𝑚 of 𝐻 and 𝟎 a column of 𝐾 zeros. More details on NMF and some applications can be found in [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the β-divergence[END_REF] , [START_REF] Heittola | Sound event detection in multisource environments using source separation[END_REF][START_REF] Ludeña-Choez | Bird sound spectrogram decomposition through Non-Negative Matrix Factorization for the acoustic classification of bird species[END_REF].

In summary, among different parameters, each taking several possible values, one tries to find an optimal NMF formulation: the number of basis in 𝑊 (𝐾 ∈ {5, 10, 20}), the number of basis in 𝑌 (𝐽 ∈ {2, 5, 10}) and the weight of the temporal regulation (𝛼 𝑠𝑚 ∈ {0, 0.5, 1, 2}). Finally, NMF is performed for 1000 iterations where 𝐻, 𝑌 and 𝑍 are initialized randomly. Sound emergences can be estimated in two different ways with semi-supervised NMF (SS NMF). First, one estimates from [𝑊𝐻] the WT noise component in 𝑉, deduces the residual component, its estimated statistical SPL 𝐿 ̃𝐴,50,𝑟𝑒𝑠. and then the estimated emergence 𝐸 ̃𝑊𝐻 , 𝐸 ̃𝑊𝐻 = 𝐿 𝐴,50,𝑎𝑚𝑏. -𝐿 ̃𝐴,50,𝑟𝑒𝑠,𝑊𝐻 .

(5)

In a second way, one considers [𝑌𝑍] as the exclusively residual noise and directly use it to estimate the sound emergence, called 𝐸 ̃𝑌𝑍 . Figure 1 summarizes this process. In Section 4.1.1, the influence of this choice on the setting on the estimation of the sound emergence is exposed.

Corpora of environmental sounds

To perform the calculation, two different corpora are used. The first one is a corpus of ambient scenes, built in order to simulate the method based on on/off cycle measurements (see Section 3.1). However, instead of considering long periods, one considers short periods of 10 minutes of ambient noise (i.e. on cycles when the WTs operate) and only residual noise (i.e. off cycles when the WTs are shutdown). This second period simulates the current process used in situ that estimates the emergence 𝐸 𝑜𝑛/𝑜𝑓𝑓 . This indicator is a baseline method, which will help to compare the performances of NMF. A second corpus composed of exclusively of WT noise spectra is dedicated to the learning of the dictionary 𝑊 and is presented in section 3.2. All the following data are expressed in Pa for each third octave band between 20 Hz and 10 kHz.

Corpus of simulated ambient scenes

A corpus of 30 simulated ambient scenes is built for this experiment. This simulation process consists in the sum of a WT noise component with a residual noise component. This process makes it possible to know the 'exact' sound emergence 𝐸 𝑒𝑥. (see Eq. 1) and can then be compared to the one estimated by NMF. By using in situ measurements, this 'exact' value would have been unknown.

Wind turbine noise samples

The WT noise components are collected during on/off cycle measurements at 300 m from the front line of a wind farm at 1.5 m height [START_REF] Kayser | A new experimental database for wind turbine noise propagation in an outdoor inhomogeneous medium[END_REF] (see Figure 2). To consider the most accurate samples, the residual noise spectra, estimated with the 30 minutes preceding and following the off measurements, filter the on measurements. From these cleaned samples, 30 samples of 10 minutes are extracted (see Figure 3). To generate ambient scenes accurately, an attenuation filter must filter the WT noise in order to simulate its propagation at a certain distance through an inhomogeneous media. 3 distances between the WTs and a receiving point are considered (𝑑 ∈ {500, 1000,1500} m) through one inhomogeneous media, which is considered here as a downward propagation condition. Propagation conditions are characterized although the vertical profile of the effective sound celerity, defined as: 𝑐 𝑒𝑓𝑓 (𝑧) = √𝛾𝑅𝑇(𝑧) + 𝑈(𝑧) cos(𝜃) (6) with 𝛾 the heat capacity ratio of dry air at constant pressure and volume (𝛾 = 1.4), 𝑅, the specific gas constant for dry air (𝑅 = 287 J.kg -1 .K -1 ), 𝜃 the angle between the wind direction and the propagation direction and 𝑇(𝑧) and 𝑈(𝑧) the mean vertical profiles of air and temperature and wind speed respectively. As turbulence is not considered, these profiles are only expressed according to their mean part:

𝑇(𝑧) = 𝑇 0 + 𝑎 𝑡 log ( 𝑧 -𝑑 𝑧 0 ) , (7.1) 𝑈(𝑧) = 𝑎 𝑢 log ( 𝑧 -𝑑 𝑧 0 ) (7.2)
where 𝑇 0 (K) is the ground surface temperature, 𝑑 = 0.66ℎ 𝑣 is the displacement height accounting for the influence of the vegetation height ℎ 𝑣 (m), 𝑧 0 = 0.13ℎ 𝑣 is the roughness height of flux profiles and 𝑎 𝑡 and 𝑎 𝑢 are the coefficients that determine the shape of the temperature and wind vertical profiles, respectively. The relative humidity is set to 70 %. The settings of the downward refraction condition used in the study are 𝑇 0 = 15°C, 𝑎 𝑢 = 100 m.s -1 , 𝑎 𝑡 = 0.15 K.m -1 , 𝜃 = 45. Sound propagation is modelled through the use of Parabolic Equation (PE) [START_REF] Kayser | Environmental parameters sensitivity analysis for the modeling of wind turbine noise in downwind conditions[END_REF] and resolved by the split-step Padé approach with an average ground floor impedance (airflow resistivity 𝜎 = 500 kN.s.m -4 ). To improve the estimation of this filter, the sound emission of WT noise is defined by the extended source model proposed by Cotté [START_REF] Cotté | Extended source models for wind turbine noise propagation[END_REF]. PE is performed for all the third octave bands from 20 Hz to 3150 Hz. This process makes it possible to obtain an accurate estimation of the SPL of WT to ensure the validity of the results of NMF.

From this resolution, the attenuation filters at the distances 𝑑 are estimated with a reference microphone at 300 m (corresponding to the position where the WT samples in the ambient scenes have been collected) (see Figure 4). Due to computation limit, the attenuations for the third octave bands from 4 kHz to 10 kHz are modelled only with the geometrical divergence and the atmospheric absorption. This limitation is not prohibitive as the WT noise contribution in these bands are very low.

Residual noise samples

The residual noise samples come from two different locations (called WF1 and WF2) and have been collected during on/off cycle measurements, but this time during the off periods. For each location, 30 samples are collected. The residual noises from WF1 account for weak residual SPL as the wind farm is located far from any anthropogenic noise sources. In the opposite, WF2 is located close to a highly frequently road.

The last 5 minutes (between the 10 th minutes and the 15 th minutes) can be seen as the transition time before the WT are shut down and are then discarded.

Building of the ambient scenes

From the WT and residual noise samples, the 30 simulated ambient scenes are built. These scenes length 10 minutes with a 1 second step. The WT noise is first propagated to 3 different distances (d ∈ {500, 1000, 1500} m) for one propagation condition (see Section 3.1.1). The WT and the residual noise samples are then summed following the Signal-to-Residual-Ratio (𝑆𝑅𝑅) defined as: 𝑆𝑅𝑅 = 𝐿 𝐴50,𝑊𝑇 -𝐿 𝐴50,𝑟𝑒𝑠.,𝑜𝑛 ,

where 𝐿 𝐴,50,𝑊𝑇 is the A-statistical SPL of the WT component and 𝐿 𝐴,50,𝑟𝑒𝑠.,𝑜𝑛 the A-statistical SPL of the residual component used for the ambient scene. The 𝑆𝑅𝑅 is defined for {-9; -6; -3; 0; 3; 6; 9} dBA. When 𝑆𝑅𝑅 < 0 dBA, the residual noise is predominant over the WT noise and inversely when 𝑆𝑅𝑅 > 0 dBA. These 𝑆𝑅𝑅 values correspond to 'exact' sound emergences 𝐸 𝑒𝑥. of {0.5; 1.0; 1.8; 3; 4.8; 7.0; 9.5} dBA respectively. Finally, the simulated ambient scenes are the spectrograms 𝑉 of NMF with dimensions F = 28 (i.e. number of third octave bands between 20 Hz and 10 kHz) and N = 600. The adjustments made to set the residual noise at a SPL corresponding to the 𝑆𝑅𝑅 are also made on the last 10 minutes of the residual noise to assure its continuity. The Figure 5 displays for 2 𝑆𝑅𝑅 values the different parts of an ambient scene and Figure 7 summarizes the steps to build the ambient scenes. Finally, a pre-processing step is performed on the ambient scenes before considering it in NMF. This pre-process consists in the attenuation of the most energetic frames. One considers here that these frames are due to high emergences of residual noises. Thus, to help the approximation made by NMF and limit the influence of the residual noise, it might be beneficial to remove them.

To do so, the temporal frames where the equivalent SPLs are superior to the 𝐿 𝐴,10 are discarded and replaced by interpolated values to ensure coherent continuity in the spectrogram. Figure 6 displays an example. To know its impact, NMF is also performed without this operation.

6-a 6-b

Computation and discussions

The objective of the experiment is to find an optimal setting of NMF (𝛽, 𝐾, 𝐽, 𝛼 𝑠𝑚 , number of iteration, estimation by [𝑊𝐻] or [𝑌𝑍], pre-processing step) on 30 simulated ambient scenes with the WT noise propagated at different distances (𝑑 ∈ {500, 1000, 1500} m), for one downward propagation condition and for two locations of residual noises (WF1 and WF2). For each ambient scene, one gets the 'exact' sound emergence 𝐸 𝑒𝑥. , thanks to the simulation process, the estimated emergence by the simulated on/off cycle scenes 𝐸 𝑜𝑛/𝑜𝑓𝑓 and the estimated one by NMF 𝐸 ̃𝑊𝐻/𝑌𝑍. . To compare these indicators, the Mean Absolute Error 𝑀𝐴𝐸 is considered. It is expressed such as:

𝑀𝐴𝐸 𝑊𝐻/𝑌𝑍 = ∑ |𝐸 ̃𝑊𝐻/𝑌𝑍,𝑚 -𝐸 𝑒𝑥.,𝑚 | 𝑀 𝑚=1 𝑀 ( 10
)
This metric is also computed with 𝐸 𝑜𝑛/𝑜𝑓𝑓 instead of 𝐸 ̃. 𝑀 is the number of considered scenes. It can be computed for the total number of scenes (𝑀= 180, see Section 4.1.1) and for each distance and residual noise (see Section 4.1.2).

Global results on all the corpus of simulated ambient scenes

The 𝑀𝐴𝐸 error is computed for each association of setting (384 associations in all) on all the corpus (𝑀 = 180 scenes). Instead of summarized all of these, only the ones with the lowest errors are displayed in Table 1 according to the dimension 𝐽, the use of the pre-processing step (no or 𝐿 𝐴,10 ) and the choice of the component in 𝑉 ̃ used to estimate the sound emergence ([𝑊𝐻] or [𝑌𝑍] ). Then, the association of settings reaching the lowest error is the optimal SS NMF formulation the most efficient that is considered. The lowest error is reached when emergences are estimated from [𝑊𝐻], with the pre-processing step and for J ∈ {2; 5}. More generally, the lowest 𝑀𝐴𝐸 errors are obtained when the emergence is estimated from [𝑊𝐻] instead of [𝑌𝑍]. This choice happens to be the most influential as the NMF settings are for the most part similar (i.e. the dimensions 𝐾 and 𝐽, the number of iteration, the smoothness weight 𝛼 𝑠𝑚 stay the same).

The addition of the pre-processing step improves the results, even if its impact is more reduced when emergences are estimated from [𝑊𝐻]. This result can be understandable as it limits on the spectrogram 𝑉 the presence of the residual noise, which can naturally help the approximation of the WT component by NMF.

Furthermore, the increase of basis in 𝑌 only improves the results when [𝑌𝑍] is considered. The approximate WT component is then constant for 𝐽 = 2 or 𝐽 = 5. This behavior reveals that the estimation of the [𝑊𝐻] component is quite robust and the increase of the dimension 𝐽 does not result in more confusion. For the next, we consider 𝐽 = 5 as it is also where the estimations of the sound emergence with [𝑌𝑍] are the best. Finally, the proposed method shows lower errors with reduced standard deviations than the baseline method on this corpus, which reveals that the estimation by SS NMF is more accurate.

The baseline method, based on on/off measurement process, makes the hypothesis that residual noise is constant over time, which might not be the case. This is particularly the case for the samples that come from the WF2 wind farm. Applied directly to the ambient scenes, this versatility does not affect SS NMF, in addition to not having to shut down the wind farm.

Behaviors of the optimal SS NMF detailed

From this main results, the behavior of the optimal SS NMF (𝐾 = 5, 𝐽 = 5, 𝛼 𝑠𝑚 = 0.5, 50 iterations, with the pre-processing step and with the emergence deduced from [𝑊𝐻] component) is now according to the distance and the location of the residual noises in Figure 9. First, one notices the strong dependence of the results to the 𝑆𝑅𝑅 values. The evolution of the errors according to the residual noise location or the distance follows the same pattern. The 𝑀𝐴𝐸 𝑊𝐻 tends to decrease when 𝑆𝑅𝑅 increases and is inferior to 0 dBA. Then, it increases strongly with the 𝑆𝑅𝑅 as the WT become more and more predominant over the residual noise.

To apprehend these results, the 1 second equivalent SPL are displayed in Figure 10 for 𝑆𝑅𝑅 = -6 dBA and for 𝑆𝑅𝑅 = 6 dBA in Figure 11.

In these representative examples, when 𝑆𝑅𝑅 is negative, the estimated WT noise is overestimated, which results in the underestimation of the estimated residual noise and so the overestimation of the sound emergence. Conversely, when 𝑆𝑅𝑅 = 6 dBA, the estimated WT noise is more equivalent to its 'exact' value. Consequently, the strong errors observed in Figure 9 are not due to a bad estimation of the WT component, but to a higher sensitivity of the emergence indicator in these cases when 𝑆𝑆𝑅 is negative. For instance, an underestimation estimation of the WT SPL of 1 dBA for 𝑆𝑅𝑅 = 9 dBA results of an error on the residual noise SPL of 4 dBA. This behavior happens specifically when the sound emergence is high (above 5 dBA) and can still be detected. If the 𝑆𝑅𝑅 is the main influence parameter, the location and the distance have also a significant impact (but reduced). According to their evolution on Figure 9, it might be some interaction phenomena between these parameters. At 500 m for 𝑆𝑅𝑅< 0 dBA, for WF1 residual samples, the errors are higher than at 1500 m. They decrease with WF2 residual noise samples. When, 𝑆𝑅𝑅 becomes positive, these patterns reverse. At 𝑑 = 500 m, residual samples from WF1 are less noisy, which results in more similar spectral shapes to the WT noise and then generates confusion. With the distance, the attenuation filters altered the spectra, which results in a better fit of the estimated and the 'exact' WT noise spectra. This observation can be seen in Figure 12 where the mean spectra of each component are displayed. Finally, when 𝑆𝑅𝑅 increases, the error becomes higher with the distance. This might be due to the difference in the propagation (± 2.0) dBA with 𝐽 = 2 and 𝑀𝐴𝐸 𝑌𝑍 = 2.5 (± 1.8) dBA with 𝐽 = 5). However, increasing this dimension also increases the risk that the cost function (Eq. 3) might be minimized mainly with this component at the expanse of [𝑊𝐻], which is moreover subject to a temporal constraint. Consequently, the dimension 𝐽 must stay low in order to not include in [𝑌𝑍] some WT noise component. Furthermore, the reduced number of spectra in 𝑌 might be limited in some cases to model a component that includes many kinds of sound sources and might be insufficient to model it properly. From [𝑊𝐻], the WT noise is more constant over time, if the dimension 𝐾 is also low, it is still sufficient to approximate this component. Consequently, deduced directly from the ambient signal, the sound diversity of the residual component is better considered even if the estimated WT component presents some mistakes. This way limits the error and then improves the estimation.

Conclusions

Semi-Supervised Non-negative Matrix Factorization (SS NMF) has been considered to estimate the sound emergence of wind turbine noise. To estimate the optimal settings, simulated ambient scenes have been built. This process makes it possible to compare the estimated emergences with their 'exact' values, which is not possible with in situ measurements without stopping the machines. The semi-supervised learning consists in adding in the dictionary a free part that can be adapted on each ambient scene. In addition, a regulation term is added through a temporal continuity constraint applied only on 𝐻 (i.e. the temporal activations of the wind turbine spectra). This constraint makes it possible to have more continuous shapes of the activations and then more similar to the wind turbine noise behavior.

The common formulation of SS NMF has been adapted to correspond to the form of the input data through an energetic sum instead of a linear sum. This change required generating new update algorithms. After computation, the settings of the optimal SS NMF have been found (𝐾 = 5, 𝐽 = 5, 𝛼 𝑠𝑚 = 0.5, 50 iterations, with the pre-processing step and with the emergence deduced from [𝑊𝐻] component), which result in a global 𝑀𝐴𝐸 error of less than 2 dBA. The simulated ambient scenes enables the study of this result according to the 𝑆𝑅𝑅, the distance and the location of the residual noise. The 𝑆𝑅𝑅 reveals to be the more influential setting as it impacts directly the sound emergence. For negative 𝑆𝑅𝑅 , SS NMF performs appropriately, but the estimation errors increase significantly as the wind turbine noise becomes more predominant. This behavior is more due to a high sensitivity of the emergence indicator than to a poor estimation of wind turbine noise. The different locations of residual noise in the scene have also an impact on the results

Figure 1 :

 1 Figure 1: block diagram of the estimation of the sound emergence by NMF.

Figure 2 :

 2 Figure 2: schema of the different measurement points used for the dictionary corpus (150 m) and the corpus of simulated ambient scenes (300 m). The distance 𝑑 corresponds to the distance where the ambient scene are located (500 m, 1000 m, 1500 m).

Figure 3 :

 3 Figure 3: example of the extraction of wind turbine spectra. The emergence on cycle is considered and filtered by the mean residual spectra estimated from the measurements carried out between the dashed lines.

Figure 4 :

 4 Figure 4: attenuation filters used to propagate the wind turbine noise at distance d with a reference receiver located at 300 m from the front line of the wind farm and at 1.5 m height above the ground.

  The residual noise samples related to this place are then much more dynamic. Each sample lasts 25 minutes. The first 10 minutes are used to build the ambient noise sample with the WT noise sample. The last 10 minutes are used as only residual noise to simulate the shutdown of the wind farm and to estimate the baseline emergence 𝐸 𝑜𝑛/𝑜𝑓𝑓 , 𝐸 𝑜𝑛/𝑜𝑓𝑓 = 𝐿 𝐴,50,𝑎𝑚𝑏. -𝐿 𝐴,50,𝑟𝑒𝑠.,𝑜𝑓𝑓

Figure 5 :

 5 Figure 5: details on the building of the simulated ambient scene for SRR = -6 dBA (5-a) and SRR = 6 dBA (5-b) (scene 6, WF1, 𝑑 = 500 m).

Figure 6 :

 6 Figure 6: influence of the pre-processing step on the spectrogram and on the A-SPL. On the initial spectrogram (6-a), the residual noise provokes high emergences. On the pre-processed spectrogram (6-b), these emergences are reduced, for instance in the intervals {327; 334} s and {386; 396} s

Figure 9 :

 9 Figure 9: errors of the optimal SS NMF details according to the distance 𝑑, the location of the residual noises and the 𝑆𝑅𝑅.

  Figure 10: a 120 seconds extract of the 1 second equivalent SPL for SRR = -6 dBA with residual noise from WF1 (scene 4, 𝑑 = 500 m) (10-Erreur ! Source du renvoi introuvable.) and WF2 (scene 4, 𝑑 = 500 m) (10-Erreur ! Source du renvoi introuvable.).

Figure 14 :

 14 Figure 14: extract of 120 seconds of the exact residual components and obtained with SS NMF from [𝑊𝐻] and from [𝑌𝑍](scene 10, WF2, 𝑑 = 500 m).

Table 1 :

 1 lower errors of SS NMF among all the tested association of settings according to the dimension 𝐽, the choice of the component in 𝑉 ̃ estimation and the pre-processing step. In bold letters, the best results.

	Method	𝑲	𝑱	𝜶 𝒔𝒎	Iteration	Pre-processing step	Component in 𝑽 ̃	𝑴𝑨𝑬 (dBA)
	Baseline	-	-	-	-	-	-	2.7 (± 2.8)
		5	2	0.5	50	no	[WH]	2.1 (± 1.8)
		5	2	0.5	50	𝐿 𝐴,10	[WH]	1.9 (± 1.7)
	SS NMF	5 5	5 5	0.5 0.5	50 50	no 𝐿 𝐴,10	[WH] [WH]	2.1 (± 1.8) 1.9 (± 1.7)
		5	2	1	5	no	[YZ]	6.2 (± 5.8)
		5	2	0.5	1000	𝐿 𝐴,10	[YZ]	3.1 (± 2.0)

. Figure 7: block diagram of the ambient scene corpus building.

Corpus for the dictionary 𝑾

To build a dictionary of WT spectra, measurements made during on/off cycle measurements on the same wind farm than in Section 3.1 are collected. To ensure to have spectra the less polluted by surrounding residual noise samples, only measurements made following the IEC standard protocol are considered as they are closed to WTs (see Figure 2). In addition, only the most emerging measurements are considered. They are then filtered by the estimated residual noise spectra measured during the off cycle as in Section 3.1.1, to ensure the most accurate WT spectra. In all, 14200 third octave band spectra (ranking from 20 to 10 kHz) are extracted. This large number of spectra cannot be directly in the multiplicative update algorithms (Eqs. 4) as it will highly increase the computation cost. Furthermore, these spectra present redundant information. Consequently, a 𝐾-mean algorithm reduces these 14200 spectra to 𝐾 ∈ {5, 10, 20}. As the ambient scenes simulate the WT noise contribution at 3 distances from the front line of a wind farm, an attenuation 'transfer function' filters is also applied on the initial dictionary at the considered distances. Instead of using PE, one chooses a simpler approach to facilitate the implementation of the proposed method in an industrial process. The Weyl's Van der Pool equation is then considered and solved by the sound rays model [START_REF] Salomons | Computational Atmospheric Acoustics[END_REF]. The sound emission models the WT noise by 3 moving monopoles (hub height = 80 m, speed rotation = 20 rpm). A downward meteorological condition is chosen with a linear vertical profile of the effective celerity: 𝑐 𝑒𝑓𝑓 (𝑧) = 𝑐 0 + 𝑏 𝑙𝑖𝑛 𝑧 with 𝑐 0 the sound celerity (m/s) with an air temperature of 10°C, 𝑏 𝑙𝑖𝑛 = 0.1 m -1 and 𝑧 the height (m). The ground impedance is defined by the Miki's model [START_REF] Miki | Acoustical properties of porous materials-Modifications of Delany-Bazley models[END_REF] with an air flow resistivity representative of countryside (𝜎 = 500 kN.s.m -4 ). The corresponding attenuation filters are displayed on Figure 8. They are generated for the 3 distances with a reference microphone corresponding to the one on the ground located at 150 m from the front line of the wind farm (see Figure 2). Finally, each filtered spectra of the dictionary is normalized such as ‖𝑊 𝑘 ‖ = 1. filter applied on the dictionary 𝑊 (see Section 3.2) emphasized with the distance and the high value of the 𝑆𝑅𝑅 where the sound emergence is more sensitive.

Comparison of the estimated residual noises

In this last part, the residual noise deduced from the estimated WT noise and the residual noise directly computed with [𝑌𝑍] are compared. To do so, the best setting with [𝑊𝐻] and [𝑌𝑍] (see Table 1) are selected and their errors according to 𝑆𝑅𝑅 are displayed in Figure 13.

Figure 13: 𝑀𝐴𝐸 errors according the 𝑆𝑅𝑅 of the optimal error of SS NMF and the one that reached the lowest 𝑀𝐴𝐸 error with the residual noise estimated with [𝑌𝑍] (see Table 1).

For both settings, the errors follow a similar respond according to the 𝑆𝑅𝑅, meaning that the estimation of the residual component is similar too. To illustrate this observation, Figure 14 summarizes 120 seconds of the residual components of one scene.

The behavior of the two residual components are similar in this example. The one deduced from [𝑊𝐻] does not present the effects of the pre-processing step (between 46 s and 66 s, see Figure 14) as it is deduced from the ambient noise and from [𝑊𝐻]. Conversely, [𝑌𝑍] is learned on processed spectrograms. Finally, both approaches reveal to be very similar and to be mostly influenced by the approximation of the WT component by [𝑊𝐻] , which provokes the overestimation of the residual component (between 74 s and 83 s for instance, see Figure 14).

By using [𝑌𝑍], the main limit is the reduced number of spectra to approximate this sound source.

The results in Table 1 reveal that the increase of the dimension 𝐽 is beneficial (𝑀𝐴𝐸 𝑌𝑍 = 3.1 according to its similarity of its audio spectra with the wind turbine noise: the more similar they are, the higher the risk of confusion.

Finally, this proposed method, despite its limitations, still generates estimation errors lower than the ones based on on/off cycle measurements, which is a great improvement for wind farm operators. Now, these results and observations make it possible to consider this method using real in situ measurements and allow a better apprehension of the estimations according to the dynamic of the ambient noise or the distance between the wind farm and the measurement point.