Franck Gilbert
email: franck.gilbert@univ-tlse3.fr

Georges Stora

Philippe Cuny

Functional response of an adapted subtidal macrobenthic community to an oil spill: macrobenthic structure and bioturbation activity over time throughout an 18-month field experiment

Keywords: Subtidal shallow sediments, Field experiment, Acute oil contamination, Macrobenthic community, Bioturbation, Sediment reworking

An experimental oil spill was carried out in order to assess in situ responses of a macrobenthic community of shallow subtidal sediments historically exposed to petroleum contamination. Both structural and functional (bioturbation activity) parameters of the community, subjected or not to a pulse acute contamination (25,000 ppm), were studied for 18 months. No difference in the community structure was detected between contaminated and control sediments, from 6 to 18 months of experimentation. Vertical distributions of organisms, however, were affected by the presence of oil contamination leading to a deeper burial of some polychaete species. In the same time, changes in sediment-reworking activity and more especially a deeper particle burying in sediments subjected to acute oil contamination were shown. These results highlight the need to complete the analysis of community structure by assessing functional aspects, such as bioturbation activity, a process integrating various aspects of benthic behaviour (e.g. feeding, locomotion, burrow building) in order to estimate real (structural and functional) and longterm effects of oil contamination on benthic communities.

Introduction

Changes in the structure and composition of macrobenthic communities driven by pollution may have marked effects on biogeochemical cycles and other important benthic ecosystem processes and functions [START_REF] Hale | Characterizations of how species mediate ecosystem properties require more comprehensive functional effect descriptors[END_REF]). As such, crude and refined petroleum hydrocarbon pollution of coastal sediments has generally been shown to induce pronounced community structural changes. These changes, which may appear at hydrocarbon concentrations above 50 ppm [START_REF] Kingston | Impact of offshore oil production installations on the benthos of the North Sea[END_REF], are generally associated to an impoverishment of the diversity concomitant with the selection of opportunistic and tolerant species (e.g. [START_REF] Gesteira | Amphipods are good bioindicators of the impact of oil spills on soft-bottom macrobenthic communities[END_REF][START_REF] Nikitik | Patterns in benthic populations in the Milford Haven waterway following the 'Sea Empress' oil spill with special reference to amphipods[END_REF][START_REF] Gesteira | Impact of the Aegean Sea oil spill on the Rîa of Ares and Betanzos fine sand community (Northwest Spain)[END_REF]. Several studies also reported a decrease in species abundances with increasing or elevated hydrocarbon contamination levels such as those observed after an oil spill (e.g. [START_REF] Carman | Does historical exposure to hydrocarbon contamination alter the response of benthic communities to diesel contamination?[END_REF][START_REF] Je | Changes in benthic communities along a presumed pollution gradient in Vancouver Harbour[END_REF][START_REF] Venturini | Petroleum contamination impact on macrobenthic communities under the influence of an oil refinery: Integrating chemical and biological multivariate data[END_REF][START_REF] Seo | The macrofaunal communities in the shallow subtidal areas for the first 3 years after the Hebei Spirit oil spill[END_REF][START_REF] Zabbey | Community responses of intertidal soft-bottom macrozoobenthos to oil pollution in a tropical mangrove ecosystem, Niger Delta, Nigeria[END_REF][START_REF] Ferrando | ): experimental evidence of low resistance capacities of benthic systems without history of pollution[END_REF]. Apart from the concentration, many other factors can, however, modulate the response of the community like the composition and toxicity of crude oils, the characteristics of the sedimentary matrix (e.g. grain size range, organic matter content; [START_REF] Je | Changes in benthic communities along a presumed pollution gradient in Vancouver Harbour[END_REF][START_REF] Venturini | Polycyclic aromatic hydrocarbons and changes in the trophic structure of polychaete assemblages in sediments of Todos os Santos Bay, Northeastern, Brazil[END_REF] and the bioavailability of hydrophobic contaminants [START_REF] Toro | Technical basis for establishing sediment quality criteria for nonionic organic chemicals using equilibrium partitioning[END_REF]. Initial benthic community structure and variable tolerances within and among species can also strongly affect the resistance capacities of the benthic system [START_REF] Bickman | Somatic and heritable effects of environmental genotoxins and the emergence of evolutionary toxicology[END_REF][START_REF] Depledge | Genetic ecotoxicology: an overview[END_REF][START_REF] Hawkins | Population responses to contaminant exposure in marine animals: infuences of genetic diversity measured as allozyme polymorphism[END_REF][START_REF] Carman | Does historical exposure to hydrocarbon contamination alter the response of benthic communities to diesel contamination?[END_REF]. This variability may be in part explained by differences in the life Responsible editor: Robert Duran Electronic supplementary material The online version of this article (doi:10.1007/s11356-014-3906-4) contains supplementary material, which is available to authorized users.

style and feeding strategies between deposit feeders and carnivores [START_REF] Venturini | Polycyclic aromatic hydrocarbons and changes in the trophic structure of polychaete assemblages in sediments of Todos os Santos Bay, Northeastern, Brazil[END_REF]. Historical or chronic exposures to hydrocarbon contamination of benthos lead to an adapted community with higher proportion of more tolerant species and/or increased tolerance among individual species [START_REF] Carman | Does historical exposure to hydrocarbon contamination alter the response of benthic communities to diesel contamination?[END_REF]. Instead, communities from pristine ecosystems seem to be particularly sensitive to oil spills [START_REF] Ferrando | ): experimental evidence of low resistance capacities of benthic systems without history of pollution[END_REF]. Whatever the previous history of contamination, the community shifts and the removal of organisms with specific biological traits induced by the contamination may have important consequences on the functioning of the whole ecosystem. For instance, the experimental removal from marine soft sediments of large suspension feeders or deposit feeders was shown to influence the flux of nitrogen and oxygen, surficial sediment characteristics, and also community composition [START_REF] Thrush | Functional role of large organisms in intertidal communities: community effects and ecosystem function[END_REF]. [START_REF] Kristensen | Influence of benthic macrofauna community shifts on ecosystem functioning in shallow estuaries[END_REF] also demonstrated that shifts in benthic fauna composition and species functional traits from a marine fjord could affect biogeochemical cycling resulting in modulation of primary productivity in the overlying water column with a feedback effect on the benthic system. Indeed, macrobenthic organisms have a pivotal role in the functioning of benthic ecosystems notably through their bioturbation activities (sediment reworking and ventilation, i.e. transport and mixing of particles and solutes, respectively; [START_REF] Kristensen | What is bioturbation? The need for a precise definition for fauna in aquatic sciences[END_REF]. Due to their reworking activity, dwelling macroorganisms introduce temporal and spatial heterogeneity in the sediment system resulting in the modification of the expression of the main microbial activities, which in turn influences the degradation of organic matter (organic contaminants) (for a review, see [START_REF] Cuny | Biodegradation of hydrocarbons in bioturbated marine sediments[END_REF]. They notably ensure the oxygenation of the superficial sediments by the construction and ventilation of burrows and tubes [START_REF] Aller | Bioturbation and remineralization of sedimentary organic matter: effects of redox oscillation[END_REF][START_REF] Pischedda | Imaging oxygen distribution in marine sediments. The importance of bioturbation and sediment heterogeneity[END_REF][START_REF] Pischedda | Spatial oxygen heterogeneity in a Nereis diversicolor irrigated burrow[END_REF].

To date, most of the studies dealing with the effects of hydrocarbon contamination on macrobenthic communities were focused on the changes in species densities and faunal composition. In this study, time changes of macrobenthic community as well as of reworking activity (functional descriptor of macroorganism activity) induced by a field experimental oil spill in subtidal chronically contaminated sediments was assessed during 18 months.

Materials and methods

Experimental site

The experimental was carried out in the Gulf of Fos (Mediterranean Sea), a highly urbanized and industrialized site, notably surrounded by chemical and petroleum plants, characterized by heavy ships and tanker traffic (Fig. 1). In the Gulf of Fos, sediments are recognized having a long-time, from moderate to high level, contamination history by pyrogenic and petrogenic hydrocarbons [START_REF] Mille | Hydrocarbons in coastal sediments from the Mediterranean Sea (Gulf of Fos area, France)[END_REF]). More especially, in the experimental subtidal site (northern Carteau Cove; 43°23′ 40″ N, 4°51′ 35″ E) organic rich muddy sand sediments contain ~1240 mg kg -1 sed. and ~70 mg kg -1 of extractable organic matter (EOM) and total hydrocarbon content (THC), respectively [START_REF] Mille | Hydrocarbons in coastal sediments from the Mediterranean Sea (Gulf of Fos area, France)[END_REF].

Experimental and analytical procedures

At the beginning of experiment (December 1994), 18 PVC cores (diameter 10 cm, height 25 cm) were separately embedded into sediments by divers, at 5-m water depth. Then, a frozen "cake" (diameter 10 cm; 1 cm thick) was deposited at the surface of each portion of sediment ecosystem delimited by a core. For half of the cores (OIL cores), the deposited cakes were a mixing of sediment and Brut Arabian Light topped at 250 °C petroleum (25.2 g oil kg -1 wet sediment). The rest of cores received a clean sediment cake and served as uncontaminated controls (UNC-Cores). In order to quantify the sediment reworking, particulate inert tracers (luminophores; 63-125 μm; 8 g) were added to the different cakes.

Then, after 6, 12, and 18 months, three cores from each sediment type (oil-contaminated, control) were sampled by divers, transported to the laboratory, and sliced in 2-cm-thick segments. Then, each sediment slice was divided in eight different equal parts.

For macrofauna species identification, four parts of each slice were randomly selected and pooled while sieved through a 500-μm mesh. Retained organisms were then stored into a plastic bottle with fixating solution (4 % buffered formaldehyde with Rose Bengal) and further sorted and counted under optic microscopes. For luminophore quantification, two other parts of each slice were randomly selected, pooled, freeze-dried, and homogeneously mixed. Luminophores' counting in each layer was then realized under UV light using an epifluorescence microscope. The quantification of sediment reworking by the gallery-diffusor model [START_REF] François | Functional approach to sediment reworking by gallery-forming macrobenthic organisms: modeling and application with the polychaete Nereis diversicolor[END_REF] was based on each vertical luminophore profile measured in each core [START_REF] Duport | Benthic macrofauna and sediment reworking quantification in contrasted environment in the Thau Lagoon[END_REF]). This model allows describing both the biodiffusion-like transport (D b coefficient) due to the continuous displacement of the tracers and the nonlocal bioadvective displacement of the tracers (r coefficient). The best fit between observed and modelled tracer distribution with depth (i.e. producing the best Db and r coefficient couple) was estimated by the least squares method.

For hydrocarbon quantification, the two remaining parts of each slice were pooled. The extraction of petroleum hydrocarbons from the sediment was achieved by alkaline digestion as described by [START_REF] Mille | Effects of oxygen on hydrocarbon degradation studies in vitro in surficial sediments[END_REF]. The fractionation and capillary gas chromatography analyses of the saturated fraction of hydrocarbons (SF), used a as proxy to calculate oil budget and vertical repartition with time, were those used by [START_REF] Gilbert | Hydrocarbon influence on denitrification in bioturbated Mediterranean coastal sediments[END_REF].

Similarity between macrobenthos samples was estimated using the Bray-Curtis index, a measure of the Bray-Curtiss dissimilarity [START_REF] Bray | An ordination of upland forest communities of southern Wisconsin[END_REF].

Differences between biodiffusion coefficients (Db) and bioadvection coefficients (r) were each analyzed by two-way ANOVA with time and conditions (contaminated or control) as factors, after assessing homogeneity of variance and normally distributed residuals using Barlett's and Shapiro-Wilk's tests, respectively.

Results

Contamination state of sediments

As indicated earlier, in the frame of our work, the saturated fraction of hydrocarbons (SF) was used as a proxy to calculate oil budget and vertical repartition with time. This was based on the assumption that all components of the Arabian Light crude oil (i.e. saturates, aromatics, resins, asphaltenes and waxes) were submitted to the same transport processes within sediments. To be consistent with this simple model, and despite previous studies carried out in the same area have shown that biodegradation processes were active in the sediments [START_REF] Grossi | Burial, exportation and degradation of acyclic petroleum hydrocarbons following a simulated oil spill in bioturbated Mediterranean coastal sediments[END_REF][START_REF] Miralles | Effects of spilled oil on bacterial communities of mediterranean coastal anoxic sediments chronically subjected to oil hydrocarbon contamination[END_REF], we also considered the biodegradation as negligible compared to abiotic losses. Therefore, by now, we will only refer to oil in the text.

In order to allow long-time lasting contamination at the experimental site (Plante- [START_REF] Plante-Cuny | Experimental field study of the effects of crude oil, drill cuttings and natural biodeposits on microphyto-and macrozoobenthic communities in a Mediterranean area[END_REF]), but without confining sediments [START_REF] Gilbert | In situ bioturbation and hydrocarbon fate in an experimental contaminated Mediterranean coastal ecosystem[END_REF], an initial 25.2 g oil kg -1 wet sediment oil contamination was applied. Indeed, results showed an important loss (i.e. resuspension and release into the water column) of added oil in the experimental oiled sediments of 81, 91 and 96 %, after 6, 12 and 18 months. This indicated, however, that highly polluted conditions (i.e. from 25.2 to 1 g oil kg -1 ) were present throughout the experiment within oil-contaminated cores (Table 1). Remaining oil was found down to 10 cm deep but with a vertical of remaining hydrocarbons different with time (Fig. 2). After 6 months, oil presented a regular decrease with depth down to 6-cm depth (99 % of remaining oil). Then, after 12 and 18 months, oil repartition highlighted a subsurface peak, layer 2-4 cm containing 63 % of remaining oil.

Macrobenthic community structure

Within the experimental sediments, the polychaetes (56.4 to 89.9 % of the total organisms) and the crustaceans (4.7 to 35.6 %) represented the two dominant macrobenthic groups. Molluscs and echinoderms were also present (see Online Resource 1). Throughout the experiment, global macrofauna density and specific richness in experimental sediments varied from 1336 to 2413 ind. m -2 and from 10 to 38, respectively. None of these parameters of the community structure did show any significant difference with both time and type of sediments (Figs. 3 and4).

Cluster analysis based on Bray-Curtis similarity distances indicated communities with affinities close or higher than 50 %, however, suggesting that the communities grouped with sampling time (Fig. 5).

Vertical repartitions of the organisms indicated that they were mainly (more than 75 %) located in the first 8 cm of sediment but with some individuals found down to 22-cm depth (Fig. 3). On a general way, in both contaminated or control sediment, the presence of organisms decreased regularly with depth. Contaminated sediments, however, showed a reduced presence of organisms in the surface layer, but a subsurface density peak centred on the 4-6-cm depth layer, from the 6th to the 18th month of experimentation.

Macrobenthic community sediment reworking

Luminophores deposited at the sediments' surface were tracked within the sedimentary columns, with time. Whatever the absence or presence of oil contamination or the time of experiment, the tracers were never found deeper than 16 cm deep and the majority of them (more than 95 %) in the first 8 cm (Fig. 6).

In contaminated sediments, a subsurface repartition peak (layer 2-4 cm) was present throughout the experiment. On the other hand, after 6 months, the distribution of tracers initially deposited at the sediment surface decreased exponentially with depth in control sediments. After 12 months, a subsurface peak was detected in the 2-4 layer but was finally less marked at the end of experiment (18 months). Quantification of sediment reworking showed mean Db and r coefficients, respectively, ranging from 2.1 to 10.1 cm 2 year -1 and 1.1 to 3.7 year -1 , respectively (Table 2). Biodiffusive coefficient Db, but not the bioadvective coefficient r, was found to significantly decrease with time (P>0.001; F 22.451; df 15). No significant difference in both coefficients was found between control (UNC) and contaminated (OIL) sediments. The mean bioadvective coefficient, however, was higher in presence of oil contamination (2.98 year -1 ) than for the control sediments (1.98 year -1 ).

Discussion

Our results indicated very similar macrobenthic communities, from both qualitative and quantitative points of view, in controls and oiled contaminated sediments 6, 12 and 18 months following the initial experimental contamination. Both [START_REF] Gandra | Experimental study of oil impact and its consequences upon benthic macrofauna in the Estuary of Patos Lagoon (RS-Brazil)[END_REF] and [START_REF] Egres | Effects of an experimental in situ diesel oil spill on the benthic community of unvegetated tidal flats in a subtropical estuary (Paranaguá Bay, Brazil)[END_REF] 2006), followed by recoveries of population levels as soon as after 8 or 2 days, respectively. Such very fast recovery was suggested resulting from larval recruitment and surrounding migration of adults and juveniles [START_REF] Carman | Does historical exposure to hydrocarbon contamination alter the response of benthic communities to diesel contamination?[END_REF][START_REF] Filho | Recolonization of infauna on a tidal flat: an experimental analysis of modes of dispersal[END_REF]. The contaminations applied were, however, far lower than our 25,000-ppm oil pulse. They resulted in pre-impact oil concentration reestablishment after only 73 days [START_REF] Egres | Effects of an experimental in situ diesel oil spill on the benthic community of unvegetated tidal flats in a subtropical estuary (Paranaguá Bay, Brazil)[END_REF]) while our experimental sediments were still highly contaminated 18 months after initial oil deposition. This 'pulse' disturbance [START_REF] Glasby | Sampling to differentiate between pulse and press perturbations[END_REF] allowing rapid recovery may have also occurred in our experimental conditions. Our 6-month sampling interval, however, was not appropriate to allow such detection. Furthermore, the investigated Carteau Cove community had a quite long history of contamination exposure and could have been potentially able to better endure contamination stress than communities with lower or almost absent previous contamination [START_REF] Carman | Does historical exposure to hydrocarbon contamination alter the response of benthic communities to diesel contamination?[END_REF]. Community structures were similar between controls and contaminated sediments, but organisms were retrieved deeper in presence of oil. This macrofaunal vertical migration activity in presence of deposited material in sediments has been described as an efficient way for organisms to temporarily avoid unsuitable conditions [START_REF] Bolam | Burial survival of benthic macrofauna following deposition of simulated dredged material[END_REF]. Within a community, this capacity is, however, not shared by all species (e.g. [START_REF] Miller | Detrimental effects of sedimentation on marine benthos: what can be learned from natural processes and rates?[END_REF][START_REF] Hinchey | Responses of estuarine benthic invertebrates to sediment burial: the importance of mobility and adaptation[END_REF]. As a matter of fact, the presence of high amounts of oil in the surface leads to a reduction of the oxygen diffusion to the deeper layers of sediments that could be critical for some species. In the frame of our work, the polychaetes Nereis caudata, Notomastus latericeus, Paradoneis lyra and Tharyx marioni were the species that successfully vertically migrated in presence of oil.

Changes in organisms' repartition into the sedimentary column were also accompanied by changes in global community mixing pattern: the control uncontaminated sediments presented a rather biodiffusive-like mixing pattern while the mixing of particulate tracers in oiled sediments tended to be non-local and more especially of convoying type (i.e. with production of a subsurface peak). It is possible that the deeper burial of organisms may have affected sediment reworking modes and rates (e.g. by a deeper extension of tracers spreading), as demonstrated in case of predation avoidance [START_REF] Maire | Indirect effects of non-lethal predation on bivalve activity and sediment reworking[END_REF]. Switches in feeding modes have been previously demonstrated for some infaunal species as a function of biological or environmental factors [START_REF] Buchanan | A comparative study of some features of the biology of Amphiura filiformis and Amphiura chiajei (Ophiuroidea) considered in relation to their distribution[END_REF][START_REF] Lopez | Ecology of deposit-feeding animals in marine sediments[END_REF][START_REF] Lindsay | Tissue loss induces switching of feeding mode in spionid polychaetes[END_REF][START_REF] Fodrie | Intraspecific density regulates positioning and feeding mode selection of the sand dollar Dendraster excentricus[END_REF][START_REF] Riisgård | Particle capture mechanisms in marine suspension-feeding invertebrates[END_REF]. Similarly, the general change in global community mixing in presence of surface oil contamination may have been linked to switches of individual modes of mixing by the totality or part of species. Actually, among the four species that were found buried deeper in oiled sediments, Notomastus latericeus, Paradoneis lyra and Tharyx marioni have the ability of both surface or subsurface feeding modes [START_REF] Fauchald | The diet of worms: a study of polychaete feeding guilds[END_REF][START_REF] Sánchez-Moyano | Distribution and trophic structure of annelid assemblages in a Caulerpa prolifera bed from southern Spain[END_REF]. However, without any available supporting data on individual species reworking, caution must be applied. Indeed, as a function of the stressor applied, species which have the capacity to switch in feeding mode may not ex abrupto, or even never perform this ability. Moreover, a switch in feeding mode may not affect the sediment-reworking mode noticeably depending on factors such as mixing intensity or measurement time. For instance, Hediste diversicolor, described as a multiple feeding mode polychaete (e.g. [START_REF] Riisgård | Filter-feeding in the polychaete Nereis diversicolor: a review[END_REF][START_REF] Scaps | A review of the biology, ecology and potential use of the common ragworm Hediste diversicolor (O.F. Müller) (Annelida: Polychaeta)[END_REF][START_REF] Riisgård | Particle capture mechanisms in marine suspension-feeding invertebrates[END_REF], did not show any change in reworking mode when experimentally exposed to a surface oil contamination (Brut Arabian Light 250; from initial 4.5 g kg -1 dry sed. to 1 g kg -1 dry sed. after 45 days; [START_REF] Gilbert | The in vitro influence of the burrowing polychaete Nereis diversicolor on the fate of petroleum hydrocarbons in marine sediments[END_REF]). On the other hand, according to biodiffusion and bioadvection mixing coefficients recently calculated from the [START_REF] Gilbert | The in vitro influence of the burrowing polychaete Nereis diversicolor on the fate of petroleum hydrocarbons in marine sediments[END_REF] data (mean control Db/r 0.8/7.6; mean contaminated Db/r 0.1/ 4.1), sediment-reworking intensity by this experimental H. diversicolor population was significantly lowered in presence of oil. In the frame of our experimental in situ work, no significant reduction of sediment mixing by the macrobenthic community was found under contaminated conditions. Nevertheless, results clearly suggested a deeper burying of tracers in presence of oil as indicated by higher bioadvective coefficients compared with uncontaminated controls.

Conclusion

Despite high oil concentration still present, no difference in the macrobenthic community was detected between experimentally contaminated and natural control sediments historically exposed to petroleum contamination. Contrary to the structure of community, organisms' vertical distribution was affected by the presence of oil contamination leading to a deeper burial (burial survival) of some polychaete species. This was accompanied by changes in sediment-reworking activity and more specifically by a deeper particle burying in sediments subjected to acute oil contamination. These results highlighted the need to complete the analysis of community structure by assessing functional aspects such as bioturbation activity, a process integrating various aspects of benthic behaviour (e.g. feeding, locomotion, burrow building), in order to estimate real (structural and functional) and long-term effects of oil contamination on benthic communities. Even if long-term toxicity effects could not be discarded, shallow subtidal macrobenthic communities adapted to chronic hydrocarbon inputs seemed to be particularly resistant to elevated hydrocarbon concentrations. In case of an oil spill, they would ensure the functioning of the benthic system by adapting their behaviour and maintaining a similar level of bioturbation activity. Thus, the oxygenation of superficial sediments and the stimulation of aerobic oil-degrading bacteria would notably favour the natural attenuation of the sediments [START_REF] Cuny | Influence of bioturbation by the polychaete Nereis diversicolor on the structure of bacterial communities in oil contaminated coastal sediments[END_REF]).

Fig. 1

 1 Fig.1Experimental site. Schematic map of the Gulf of Fos (Mediterranean Sea) indicating the experimental site (black diamond) and nearby oil and gas refining and chemicals main facilities (grey spots)

Fig. 4

 4 Fig. 2 Sediment contamination. Vertical repartition of remaining oil within contaminated (OIL) sediments, after 6, 12 and 18 months. Values are mean±SD for triplicates

Table 1

 1 Sediment contamination

	Period of experiment	Estimated oil content (g oil kg -1 sed)
	Start to 6 months	25.2 to 4.8
	6 to 12 months	4.8 to 2.3
	12 to 18 months	2.3 to 1.0
	Oil content, within contaminated (OIL) sediments, after 6, 12 and
	18 months. Oil content was estimated from the saturated hydrocarbon
	fraction budget for an initial oil content of 25.2 g kg -1 sed

Acknowledgments Thanks are due to the divers Claude Jalong and Jacques Millet, and the crew of the N/O Antedon who provided vital support in the field. We are also grateful to Sean Fitzpatrick and anonymous referees for constructive comments on a previous draft of the manuscript. This work was funded by the programme "GDR HYCAR n°1123: Cycles biogéochimiques des hydrocarbures naturels et anthropiques en milieu matin" and supported by the Centre National de la Recherche Scientifique (CNRS), Universities and the Société Elf Aquitaine.


Etienne Fournier

Member, IEEE Antoine Picot

Member, IEEE Jérémi Régnier

Mathias Tientcheuyamdeu

Jean-Marie Andréjak

Member, IEEE Pascal Maussion

M Picot

Current-Based Detection of Mechanical Unbalance in an Induction Machine Using Spectral Kurtosis with Reference

Keywords: Spectral Kurtosis, Condition Monitoring, Statistics, Healthy Reference, Normalized Fault Indicator, Mechanical Unbalance, Induction Motor, Instantaneous Frequency, Fault Diagnosis

This article explores the design, on-line, of an electrical machine's healthy reference by means of statistical tools. The definition of a healthy reference enables the computation of normalized fault indicators whose value is independent of the system's characteristics. This is a great advantage when diagnosing a broad range of systems with different power, coupling, inertia, load, etc. In this paper, an original method called spectral kurtosis with reference is presented in order to design a system's healthy reference. Its principle is first explained on a synthetic signal. This approach is then evaluated for mechanical unbalance detection in an induction machine using the stator currents instantaneous frequency. The normalized behaviour of the proposed indicator is then confirmed for different operating conditions and its robustness with respect to load variations is demonstrated. Finally, the advantages of using a statistical indicator based on a healthy reference compared to a raw fault signature are discussed.

D IAGNOSIS and condition monitoring of electrical ma- chines has become a major topic of concern among industrial and academic research. Unexpected failures in electromechanical systems may indeed lead to important losses of production, safety issues and additional costs. For these reasons, developing efficient and robust diagnosis systems is a priority for electrical machines' manufacturers.

Numerous monitoring methods are based on vibratory signals' analysis and enable an efficient detection of mechanical faults [START_REF] Prieto | Bearing Fault Detection by a Novel Condition-Monitoring Scheme Based on Statistical-Time Features and Neural Networks[END_REF]- [START_REF] Immovilli | Bearing Fault Model for Induction Motor With Externally Induced Vibration[END_REF]. However, these methods require accelerometers which are expensive devices. They are only profitable for high power systems where vibration sensors represent a negligible part of the total cost. Therefore, research has focused on diagnosing machine failures through the analysis of current signals due to their availability for control purposes. These methods are mainly based on the Motor Current Signature Analysis (MCSA). A review of MCSA techniques can be found in [START_REF] Bellini | Advances in Diagnostic Techniques for Induction Machines[END_REF]. These methods generally focus on failures generating fault signatures at specific frequencies. These harmonics are predicted by MCSA and can be tracked with signal processing tools in order to obtain specific fault indicators [START_REF] Ebrahimi | Advanced Eccentricity Fault Recognition in Permanent Magnet Synchronous Motors Using Stator Current Signature Analysis[END_REF]- [12]. In the case of mechanical unbalance detection, the current spectral signatures at frequencies f f ± f r proved to be efficient fault indicators [13]- [14]. MCSA may also be used on other electrical quantities extracted from stator currents such as active and reactive power [15], current instantaneous frequency [16]- [17] or demodulated current signals [18].

However, recent researches pay little attention to unavoidable constraints faced in diagnosing a manufacturer's broad range of machines. Most works focus on the study of a single machine operating at a particular functioning point but MCSA fault indicators may vary a lot for healthy and faulty conditions depending on the system's characteristics and operating point. For example, Bellini et al. showed in [19] the influence of the load inertia level on fault signatures related to broken bar detection. It has also been shown in [13] that the load level could significantly affect current signatures linked to mechanical unbalance. Thus, classification systems based on neural networks (NN) or fuzzy logic (FL) have been developed in order to automatically determine a machine's health state while taking into account its characteristics and operating conditions. Relevant applications of NNs based on stator currents are presented in [20]- [23]. Bearing fault detection is also studied by Raj et al. in [24] with a FL-based surveillance system. All these classification methods require faulty and healthy learning sets at different operating points for training. During the functioning time of the system, each new data set is compared to the different learning sets in order to determine the machine's health condition. Fault severity may as well be evaluated if faulty learning sets have been realized for different levels of the considered degradation. Unfortunately, it seems impossible for a machine's manufacturer to obtain such an amount of information, especially faulty learning sets, for every machine model. Moreover, these information may depend on the machine environment such as the load inertia or the coupling for example.

In this context, it seems more profitable to create a machine's healthy reference during its initial working period by statistical means and with a protocol which takes into account the influence of the operating point. After this reference period, each new recording will be tested to find out whether the system has drifted from its referenced behaviour or not. The aim is therefore to detect faulty workings and the estimation of the fault severity is not targeted here since no faulty learning sets are used. The present article explains how new implementations of the spectral kurtosis (SK) make this kind of monitoring possible. SK is a signal processing tool complementing the Power Spectrum Density (PSD) representation. Its rules and properties are explained in [25]. Recent works [26]- [28] have shown that SK may be an efficient tool for monitoring faults generating impulses in vibration or current signals such as bearing or gear faults. Moreover, SK can also be successfully used to detect generalizedroughness bearing fault, as shown in [29]. However, SK is unsuitable for detecting long term changes within a signal in its classic implementation. Indeed, a slowly increasing spectral component is considered as stationary on each local part of the signal by the SK. A different approach is presented in a previous work [30] with new protocols of computation of the SK. These SK calculations no longer use only a single signal's recording but several recordings in order to detect a developing mechanical unbalance in an induction machine (IM) via stator phase currents.

In this context, this article aims to enhance the design of normalized fault indicators independent of a machine's type, power and operating point by using a statistical reference. This healthy reference is provided thanks to the spectral kurtosis with reference (SK R ) presented in [30] and detailed in this paper. Moreover, new experimental tests have been carried out at several load conditions and current instantaneous frequency (IF) is now used to monitor mechanical unbalance development. Applications targeted in this paper are indus-trial systems such as fans, pumps or compressors driven by induction motors. These latter are fed by open-loop inverters which enable the selection of one or several supply frequencies dedicated to the monitoring phases. In Section II, definition and properties of the classic implementation SK are presented to remind its principle of operation. In section III, the SK R computation protocol is designed in order to overcome the limitations of the SK. SK R performance are then tested in Section IV for detection of a developing mechanical unbalance via current IF. These experiments are performed at several load conditions in order to verify the normalized character of SK R indicators. SK R robustness when faced to load variations is also tested. Finally, advantages of SK R for the diagnosis of a broad range of machine are discussed and perspectives are detailed.

II. MATHEMATICAL BACKGROUND AND CLASSICAL USE OF SPECTRAL KURTOSIS A. Kurtosis

Let x be a random variable with mean µ and standard deviation σ. In high order statistics, the excess kurtosis of the random variable x is defined in (1) as

K(x) = κ 4 κ 2 2 (1)
where κ i is the i th cumulant of x and is defined as

κ i = m i - i-1 k=0 k -1 i -1 κ k m i-k (2)
with m i is the i th moment of x. For clarity reasons, the excess kurtosis will be simply referred as kurtosis from now on.

In statistics, kurtosis is the second shape factor (after the skewness) and reflects the sharpness of a distribution. The sharper the distribution of x, the higher its kurtosis. Reciprocally, the flatter the distribution, the lower its kurtosis. However, kurtosis may also be used to detect outlier elements within a distribution. For instance, let us consider a sample of a random variable x whose probability density function (PDF) is illustrated in Fig. 1 (top). The kurtosis of this sample, where there is no outliers and whose elements are well gathered around the sample's mean, has a low value (examples: 0 for a Gaussian sample, -1,2 for a uniform sample). Now, let us consider the same sample in which were added some extreme (and rare) elements around 10. Its PDF is plotted in Fig. 1 (bottom). This time the sample's kurtosis has a high value caused by the outliers' presence. This property is the basic principle used in SK computation for detecting non-stationary components within a signal.

B. Classic Spectral Kurtosis

Let x[k] ∈ R be a digital signal and X(f ) its Discrete Fourier Transform. Its spectral component at a specific frequency f i is a complex random variable, noted X(f i ). Thus, the spectral kurtosis of x(k) is noted SK x (f ) and is defined as the kurtosis of the random variable X(f i ) at each frequency Therefore, a minimum number of spectra is required to properly estimate a signal's SK. In practice, this number is reached by using the STFT whose principle is to split the signal x into N segments (with or without overlap) and to compute a FFT on each one. When a spectral component |X(f )| has a stationary behavior, its distribution is similar to the one presented in Fig. 1 (top). Thus, SK x (f ) has a low value. On the contrary, when this spectral component is non-stationary, its distribution is closer to the one presented in Fig. 1 (bottom). More explanations about the principle of operation of SK and its computation protocol are detailed in [30].

f i [26] SK x (f i ) = K(X(f i )) = κ 4 (X(f i )) κ 2 2 (X(f i )) (3)
This algorithm is therefore suitable to detect fault signatures with a non-stationary behaviour during a recording period, such as bearing faults. However, an innovative SK implementation is presented in the following section in order to design normalized fault indicators for other types of degradation.

III. DEFINITION OF A MACHINE'S HEALTHY REFERENCE VIA A NEW IMPLEMENTATION OF SPECTRAL KURTOSIS

In the case of the monitoring of an electrical machine, phase currents (or other useful measures) are recorded for a certain amount of time T recording at regular intervals, once a day for example. Each recording contains L samples. As shown in section II, SK would be an efficient indicator to detect fault signatures with a non-stationary behaviour over time T recording , such as some bearing faults signatures. However, many faults occurring in electro-mechanical systems are more likely to produce signatures with a stationary behaviour. Indeed, a mechanical unbalance slowly increasing in a ventilation system will induce fault signatures growing over large time scale but locally constant on each recording. A shaft misalignment provoked by a maintenance phase would do the same but more suddenly. Therefore, this section focuses on the design of normalized fault indicators with the aim of detecting two successive increases of the level of a simulated mechanical unbalance.

In order to simulate a mechanical unbalance in a rotating system, let us consider a synthetic signal IF (t) representing the current instantaneous frequency. To obtain this variable, I α (t) and I β (t) are first computed by applying Concordia matrix on phase currents I 1 (t), I 2 (t) and I 3 (t). Instantaneous frequency (IF) and instantaneous amplitude (IA) of stator currents are then defined as

IF (t) = 1 2π dφ(t) dt (4)
IA(t) = |I c (t)| (5)
with

I c (t) = I α (t) + iI β (t) and φ(t) = arg(I c (t)
). The authors of [14] indeed proved that faults producing load torque oscillation also cause phase modulations (PM) of stator currents. Thus, IF is a relevant quantity for detecting these kinds of faults because it reveals PM harmonics of currents.

The impact of a mechanical unbalance on the current IF can reasonably be modelled as

IF (t) = f f + β • cos(2πf r t + φ r ) + n(t) (6)
This signal is composed of: 1) a constant component of amplitude f f → the current fundamental frequency; 2) a fault component with an amplitude equal to β at the rotation frequency f r → relative to the unbalance level; 3) a noise component n(t). In this example, the signal IF (t) is composed of 400 5srecordings and the fault severity β evolves over the recordings as shown in Fig. 2. This profile is representative of two

A. Spectral Kurtosis

First, the classical implementation of the SK is tested on the synthetic signal IF . It is computed on each IF recording by using STFT windows of 8000pts with 75 % overlap. Specific explanations on the choice of the window length and the overlap are given in [25] and [27]. The SK of the simulated current IF at frequency f r is presented in Fig. 3. It can be observed that the indicator SK(f r ) decreases from about 0 to -0.8 after the first increase of the modulation excessive because its length also conditions the duration of the reference phase. A good compromise between a trustworthy reference and a short referencing phase is N ref = 50. Now let us consider the choice of the parameter N mov . As it is visible in Fig. 5, the reaction of indicator SK R is strongly related to the ratio R =

N ref
Nmov . Indeed, SK R reacts stronger for high R values but it is also noisier. In any case, the detection sensitivity is not affected by parameter N mov as long as the ratio R is much higher than 1. N mov = 1 is a good choice since it ensures the highest possible value of R for a given value of N ref . In this way, the indicator SK R tests the membership of the latest value of the considered fault signature to its reference achieved during healthy working. The couple {N ref = 50, N mov = 1} is not the unique possibility but the indicator SK R thus created is robust and easy to interpret. This choice will therefore be retained in the rest of the study.

Finally, an improvement of the SK R can be easily achieved. It may indeed be useful to subtract the kurtosis value of the reference set to the SK R value in order to correct the non-Gaussianity of the considered variable. In this way, SK R has always a quasi-zero value during healthy period, even if the random variable is not Gaussian.

IV. EXPERIMENTAL RESULTS

A. Test Bench Description

An experimental test bench has been set up in order to illustrate the efficiency of SK R for the detection of a mechanical unbalance occurring in an electro-mechanical system. The experimental bench is composed of:

1) a 5.5 kW Leroy Somer IM with two poles pairs, 28 rotor bars and a wye connection; 2) a direct current motor (DCM) used as a generator in order to produce the desired load torque level; 3) an iron disk placed on the shaft between both machines.

Several weights may be associated with it in order to create different levels of mechanical unbalance. An overall view of the test bench is presented in Fig. 6. The IM is fed by an open-loop inverter and all tests presented in this section are realized with a supply frequency f f = 40Hz. A data acquisition system is used to measure the IM phase the system's rotation axis generates load torque oscillations Γ osc (t) at the frequency rotation f r which are defined by

Γ osc (t) = m • R • g • cos(2πf r t) (7)
with g the gravitational acceleration. For example, this kind of degradation usually happens in industrial fans due to blade clogging by grease and dust accumulation or the break of a fan blade. Once detected, the mechanical unbalance may easily be fixed in order to prevent a premature wear of the mechanical components, especially bearings.

B. Tests at no load condition

250 5s-recordings of the IM currents have been measured on the machine. They are composed of:

1) 100 recordings with no mechanical unbalance and corresponding to the healthy state of the machine; 2) 50 recordings with the low level of mechanical unbalance B 1 created by a small weight m 1 = 77.5g fixed to the iron disk at a distance R = 7.5cm from the rotation axis. This level of mechanical unbalance produces load torque oscillations at f r with a theoretical amplitude of 0.15% of the IM rated torque; 3) 50 recordings with the medium level of mechanical unbalance B 2 created by a small weight m 2 = 133g fixed to the iron disk at a distance R = 7.5cm from the rotation axis. This level of mechanical unbalance produces load torque oscillations at f r with a theoretical amplitude of 0.27% of the IM rated torque; 4) 50 recordings with the high level of mechanical unbalance B 3 created by a small weight m 3 = 274g fixed to the iron disk at a distance R = 7.5cm from the rotation axis. This level of mechanical unbalance produces load torque oscillations at f r with a theoretical amplitude of 0.55% of the IM rated torque. The spectral component of the current IF at frequency f is noted ÎF (f ). It is shown in [14] that a mechanical unbalance produces a phase modulation (PM) of the current fundamental (of frequency f f ) at the rotation frequency f r . So, mechanical unbalance signature appears in the current IF spectrum at frequency f r . In order to tolerate small variations of the system operating point while monitoring spectral components, we redefine the variable IF (f ) as

IF (f ) = max{ ÎF (f -∆f : f + ∆f )} (8)
with ∆f chosen to be 0. variable for the detection of load torque oscillations due to a mechanical unbalance. Although the harmonic IF (f r ) is noisy, its value remains stable during the healthy working period. Then, its evolution well reflects the different stages of mechanical unbalance. Even the lowest unbalance level visibly affects the mean level of IF (f r ) while it produces extremely weak load oscillations. The harmonic IA(f r ) has also been illustrated in Fig. 8. It can be noticed that IA also reflects the increase of the mechanical unbalance level (especially the two highest levels), but to a much lower extent than IF. Indeed, IA is more sensitive to degradation generating amplitude modulation (AM) of phase currents, such as eccentricity faults. Thus, a joint study of IF and IA may enable the detection and differentiation of many mechanical faults occurring in a system.

However, although it may seem easy to choose a fault threshold for the fault indicator IF (f r ), the healthy level of this harmonic depends on many factors. Electrical machines with different rated power, poles number, inertia, etc... may have different healthy levels of fault signatures. In order to illustrate this issue, let us monitor the harmonic IF (f r ) during two experiments, both in healthy conditions and at the same load level but with two different total inertia. The difference of inertia in the two experiments is only due to the addition of the iron disk on the shaft. Of course, no unbalance weights were fastened on the iron disk during this test. Results of both indicators are shown in Fig. 9. It is well perceptible in Fig. 9 that both healthy indicators have different mean values µ 1 = 0.0068 and µ 2 = 0.0094 ≈ 1.4µ 1 , while the same machine was monitored at the same load level. When considering diagnosis of a broad range of machines and applications, choosing an a priori robust fault threshold is therefore tricky. It seems unavoidable to study the behavior of fault signatures during the initial working period of the system, when it can be considered as healthy. In this way, the system's reference created is thus specific to each monitored system, which is a key point of the proposed approach.

It has been shown in section III that SK R can be used both to create a system's reference and to compare each new recording to it. This second step enables the user to evaluate if the system has drifted from its healthy reference or not. In order to test SK R efficiency for the detection of a mechanical unbalance, this algorithm has been computed on the fault signature IF (f r ). This computation has been achieved according to the protocol presented in Fig. 4 with N ref = 50 and N mov = 1. Results are shown in Fig. 10 (green solid line). First, it is visible in Fig. 10 that SK R is an efficient indicator for detecting the different levels of mechanical unbalance. This statistical treatment holds the information contained in the fault signature IF (f r ). However it normalizes the healthy values around 0 which makes the detection of abnormal behaviours easier. After the reference period, the SK R (f r ) value remains close to 0 during the healthy operation period. Its value then increases with the severity of the mechanical unbalance. The application of the SK R on the fault signature IF (f r ) thus enables the creation of a normalized fault indicator. Indeed, SK R -based indicators will always have a value about equal to zero during healthy condition whereas classical energy indicators may have many different levels, depending on the system characteristics (type of motor, operating point, etc.). This point is major when considering the diagnosis of a wide range of electro-mechanical systems. Finally, the influence of the ratio R =

N ref
Nmov on the SK R is experimentally verified. SK R sensitivity is indeed directly related to the R level. The choice N ref = 50 and N mov = 1 guarantees a good detection sensitivity and an easy interpretation. It will be conserved in the rest of the study.

C. Tests at different load conditions

The same experiments have been carried out for multiple load conditions in order to confirm the previous results. As mentioned in section IV-B, the application of SK R on fault signatures enables the creation of normalized fault indicators. But also, a single fault threshold may be associated to these normalized indicators, whatever the system or the fault considered. Indeed, let SK R (f r ) be the SK R applied on current IF at frequency f r . The fault threshold t 1% is defined by

P(SK R (f r ) > t 1% ) = 1% (9)
Therefore, the normalized fault indicator SK R (f r ) exceeds t 1% only 1% of the time during healthy working conditions. If the exceeding rate of SK R (f r ) above t 1% is much more higher than 1%, the system has then drifted from its healthy reference. A maintenance phase is thus necessary. For N ref = 50 and N mov = 1, the fault threshold t 1% is about equal to 0.2.

For every experiment at the different load levels, SK R has been computed on the fault component IF (f r ). The percentage by which SK R exceeds the fault threshold t 1% is then calculated for each health state (healthy, lowest unbalance (B 1 ), middle unbalance (B 2 ) and highest unbalance (B 3 )). The results of this test campaign are presented in Table I. First, it can noticed that, as expected, the exceeding rate of SK R is always about equal to 1 % during healthy conditions. It is a key point of this reference-based monitoring system because no false alarms are generated, for all load conditions. Then, it can also be observed that the indicator sensitivity seems to decrease with the load level. Indeed, the lowest unbalance level B 1 , which is extremely low, is only clearly detectable at no load condition. Two reasons may explain the loss of SK R (f r ) sensitivity with the load level. First, energy indicators are more noisy for high load levels. Then, their faulty levels also decrease with the increase of the load level. However, both unbalance levels B 2 and B 3 seems well detectable for every load conditions (except the middle unbalance level at 70% load).

D. Toward an automatic referencing and monitoring protocol in the torque-speed plane

The previous tests proved the normalization and sensitivity properties of the SK R at different load conditions. However, a reference was necessary per operating point in order to ensure the proper functioning of the indicator. It is impossible, in practice, to associate a reference per every possible working point in the torque-speed plane because it would require a infinite number of references.

The aim here is therefore to associate a reference to a relatively large area in the torque-speed plane, and not to a singularity. To achieve this, areas must be designed such that the considered fault signature has a stable behaviour on each of them. Experimental results showed that the healthy level of the fault signature IF (f r ) varies moderately with the load level. At a specific supply frequency f f , the load axis may therefore be split into M relatively large zones and a reference is created per load area.

Experimental results also showed that the healthy level of the considered fault signature could have large variations along the speed axis. To ensure a good robustness and sensitivity of the SK R indicator, it should be split into a great number of speed zones. However, most industrial systems such as fans, pumps, or compressors, which are targeted in this study, are fed by open-loop inverters. This way, the supply frequency f f may be changed during the monitoring phases in order to always realize the diagnosis at one or a few supply frequencies.

The tests presented in section IV-C can be reorganized by selecting an adequate division of the load axis and a specific supply frequency chosen to perform the diagnosis phases. The monitoring parameters have been selected as follows:

1) the supply frequency selected to achieve the monitoring phases is equal to f f = 40Hz; 2) the load axis is split into M = 4 areas wide by 25% of the rated load.

In this way, a unique reference of the fault signature IF (f r ) is designed for each quarter of the load axis. The reference of the first quarter is built by using N ref recordings acquired at both 0% and 25% of the rated load. In the same way, the reference of the second quarter is designed with N ref values of the fault signature IF (f r ) extracted from recordings at both 35% and 45% of the rated load. Similarly, the reference of the third quarter of the load axis is achieved by using the recordings obtained at 60% and 70% of the rated load. At each new recording, the SK R applied on the fault signature IF (f r ) is computed by using the reference corresponding to the load level of the system. alarm rate is not affected at all by using a single reference per load area. The robustness of the proposed approach is therefore maintained. Moreover, the SK R sensitivity is hardly affected by the use of a reference covering a wide zone of the load axis. Indeed, the exceeding rate in a load area is often close to the mean of the exceeding rates of each functioning point included to this area. For example, the exceeding rate under the lowest unbalance level B 1 in the first load quarter is equal to 20. The exceeding rate at 0% of the rated load and the one at 25%, presented in Table I are respectively equal to 44% and 7% giving a mean exceeding rate equal to 25.5%. Both values are close and the little loss of sensitivity observed is due to the slight non-homogeneity of the variable IF (f r ) over the considered load area. Finally, this referencing process could be repeated for different supply frequencies in order to minimize the impact on the application. To conclude this experimental section, SK R seems to be an efficient monitoring protocol for the diagnosis of a broad range of induction machines operating with different characteristics (type of load, inertia, coupling, etc.). It should be recalled that the tested unbalance levels remains rather low, from B 1 = 0.15% to B 3 = 0.55% of the IM nominal torque. The test campaign presented in Table I shows that a single fault threshold with a statistical meaning is sufficient to ensure a robust detection of these different levels of mechanical unbalance without false alarms. The creation of several references (one for each load condition) enables a good detection of low levels of mechanical unbalances. Because references are created on-line, SK R indicator is therefore independent of the type of machine, of load characteristics and others parameters which can affect classical MCSA indicators as the ones presented in Fig. 8. Moreover, tests realized in section IV-D showed that this method may easily be adapted to ensure a robust diagnosis in the torque-speed plane.

V. CONCLUSION In this paper, an efficient diagnosis protocol based on the computation of a system's healthy reference has been presented. Normalized indicators well adapted for the condition monitoring of a broad range of machines were designed thanks to this statistical reference. In section II, the kurtosis has been defined and its ability to detect outliers within a reference distribution has been demonstrated. In section III, an innovative indicator, the SK R , based on the kurtosis has been designed by including a reference set. SK R ability to generate a system's healthy reference and to detect any drift from it has been tested on synthetic signals. Finally, an experimental test bench has been set-up in order to verify SK R efficiency, when applied on current IF, for the detection of low levels of mechanical unbalance. Tests performed at several load conditions have proven the SK R detection capacity with a single fault threshold for all load conditions. The creation of the system's healthy reference thus enabled a robust detection of weak mechanical unbalances without any false alarms for different operating conditions. Moreover, its robustness to load variations has been demonstrated. Future work will focus on the definition a generic monitoring protocol in the torque-speed plane which specifies the references' area of validity and optimizes the indicator's sensitivity while conserving its robustness and minimizing the impact of the monitoring protocol on the application. Moreover, the proposed approach showed good performance for mechanical unbalance detection in an IM and it will be interesting to confirm these results for other kinds of faults (eccentricity, bearing faults, etc.) and for different systems.

Fig. 1 .

 1 Fig.1. Probability density function and kurtosis value of a normal distribution without outliers (top) and with some outlier elements around 10 (bottom). The presence of outliers within a distribution is responsible for the increase of its kurtosis value.

Fig. 2 .

 2 Fig. 2. Evolution of the modulation index (or fault severity) β over the 400 recordings. This profile models two successive increases of the simulated fault severity.

Fig. 5 .

 5 Fig. 5. SK R applied on the simulated fault harmonic IF (fr) for N ref = 50 and Nmov = 1 (dark blue line with circles), N ref = 25 and Nmov = 1 (light blue line with triangles) and N ref = 50 and Nmov = 4 (magenta line with stars). The ratio R = N ref Nmov greatly impacts the indicator's behaviour.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Experimental test bench composed of a 5.5 kW induction machine (right), a direct current motor to produce the load torque (left) and a iron disk placed on the shaft to create the different levels of mechanical unbalance.

Fig. 8 .

 8 Fig. 8. Illustration of |IF (fr)| (green solid line) and |IA(fr)| (red dotted line) evolutions for the different health states of the induction machine (healthy, B 1 , B 2 and B 3 ).

2 μ 1 +40%HealthyFig. 9 .

 219 Fig. 9. Fault signature |IF (fr)| during healthy conditions with different configuration of the monitored system : with the inertia disk (light blue solid line) and without the inertia disk (dark blue dotted line).

Fig. 10 .

 10 Fig. 10. SK R applied on current IF at fault frequency fr in the healthy state and for the different unbalance levels B 1 , B 2 and B 3 . Different values of the ratio R = N ref Nmov are used and simulation results are confirmed.

  IF (t) Stator current instantaneous frequency K Kurtosis N mov Length of the tested set N ref Length of the reference set s Slip of the induction machine SK Spectral kurtosis SK R Spectral kurtosis with reference X(f ) Fourier Transform of variable X(t) at frequency f

	I. INTRODUCTION

  Etienne Fournier got his MSc in electrical engineering from Supélec Gifsur-Yvette, France in 2012. He is now a PhD Student at the Laboratory of Plasma and Energy Conversion (LAPLACE) in Toulouse, France. His Scholarship is funded under an industrial agreement with Moteurs Leroy Somer in Angoulême, France. His work focuses on preventive diagnostic of electrical machines under variable-speed control. Antoine Picot graduated from the Telecom Department of the Institut National Polytechnique (INP) Grenoble, France in 2006. He received the MSc degree in signal, image, speech processing and telecommunications in 2006 and his PhD in automatic control and signal processing in 2009 from the INP Grenoble. He is actually an associate professor at the INP Toulouse. He is also a Researcher with the Laboratory of Plasma and Energy Conversion (LAPLACE), Toulouse. His research interests are in monitoring and diagnosis of complex systems with signal processing and artificial intelligence techniques. Jérémi Régnier was born in 1975. He received the PhD degree in electrical engineering from the Institut National Polytechnique (INP) Toulouse, France in 2003. Since 2004, he has been an Assistant Professor with the Electrical Engineering and Control Systems Department, INP Toulouse. He is also a Researcher with the Laboratory of Plasma and Energy Conversion (LAPLACE), Toulouse. His research interests include modelling and simulation of faulty electrical systems as well as the development of monitoring techniques using signal-based and model-based methods. Mathias Tientcheu Yamdeu got his MSc in Electrical Engineering and Automatic control from Institut National Polytechnique (INP) Toulouse, France in 1990. He is currently a senior engineer at Moteurs Leroy Somer in Angoulême, France where he is in charge of the control and diagnostic of variable speed drives and power electronics converters. Jean-Marie Andréjak got his MSc in Electrical Engineering from EN-SEEIHT Toulouse, France in 1978. After having worked in JEUMONT-SCHNEIDER, in electronics E&D department for 10 years, he joined Moteurs Leroy Somer in Angoulême, France, in 1989. He is the Power Electronics E&D manager. Pascal Maussion got his MSc and PhD in Electrical Engineering in 1985 and 1990 from Institut National Polytechnique (INP) Toulouse, France. He is currently full Professor with the University of Toulouse and with the Laboratory of Plasma and Energy Conversion (LAPLACE), Toulouse. His research activities deal with control and diagnosis of electrical systems (power converters, drives, lighting) and with the design of experiments for optimisation in control and diagnosis. He is currently Head of Control and Diagnosis group in LAPLACE. He teaches control and diagnosis in a school of engineers.

TABLE I PERCENTAGE

 I BY WHICH THE INDICATOR SK R APPLIED ON IF (fr) EXCEEDS THE FAULT THRESHOLD t 1% FOR EACH HEALTH STATE AND

			LOAD LEVEL.
				Load level (% of rated load)
		Health state	0% 25% 35% 45% 60% 70%
		Healthy	0	0	1	0	1	0
	SK R applied on	B 1	44	7	0	0	1	2
	IF (fr)	B 2	100	100	68	42	72	9
		B 3	100	100	100	100	100	100

  The percentage by which SK R exceeds the fault threshold t 1% is then calculated for each health state (healthy, lowest unbalance (B 1 ), middle unbalance (B 2 ) and highest unbalance (B 3 )) and for the different load areas. Results are presented in Table II. It can be noticed that the false

TABLE II PERCENTAGE

 II BY WHICH THE INDICATOR SK R APPLIED ON IF (fr) EXCEEDS THE FAULT THRESHOLD t 1% FOR EACH HEALTH STATE AND

			LOAD AREA.
			Load area (% of rated load)
		Health state	[0%→25%] [25%→50%] [50%→75%]
		Healthy	0	1	1
	SK R applied on	B 1	20	1	3
	IF (fr)	B 2	99	47	41
		B 3	100	100	100


Rafael Blanquero
email: rblanquero@us.es

Emilio Carrizosa
email: ecarrizosa@us.es

Amaya Nogales Gómez

Frank Plastria
email: frank.plastria@vub.ac.be

Amaya Nogales-Gómez

Single-Facility Huff Location Problems on Networks

Keywords: Huff location models, location on networks, DC optimization, Interval Analysis

  L'archive ouverte pluridisciplinaire

Introduction

In this paper we address two competitive location models [START_REF] Blanquero | Continuous location problems and big triangle small triangle: constructing better bounds[END_REF][START_REF] Drezner | Finding the optimal solution to the Huff competitive location model[END_REF][START_REF] Fernández | Solving a Huff-like competitive location and design model for profit maximization in the plane[END_REF][START_REF] Huff | Defining and estimating a trading area[END_REF][START_REF] Huff | A programmed solution for approximating an optimum retail location[END_REF][START_REF] Plastria | Static competitive facility location: An overview of optimisation approaches[END_REF] on a network [START_REF] Dearing | Optimal locations for a class of nonlinear single-facility location problems on a network[END_REF][START_REF] Labbé | Location on Networks[END_REF]. Let N = (V, E) be a network, with node set V and edge set E. The length of the edge e ∈ E is given, and it is denoted by l e . The distance between two nodes a i , a j is calculated as the shortest path [START_REF] Labbé | Location on Networks[END_REF] from a i to a j . For each e ∈ E, with end nodes a i , a j , we identify each x ∈ [0, l e ] with the point in the edge e at distance x from a i and distance l e -x from a j . This way, we obtain that, for any vertex a k ∈ V , the distance d(x, a k ) from x to a k is, as a function of x, a concave piecewise linear function, given by: d(x, a k ) = min{f ia k (x), f ja k (x)}

(1)

f ia k (x) = x + d(a i , a k ) f ja k (x) = (l e -x) + d(a j , a k )
We refer the reader to [START_REF] Labbé | Location on Networks[END_REF] for a comprehensive introduction to location models on networks.

The remainder of this paper is structured as follows. Two different location models on networks are considered. The first model (Section 2) is the classic Huff location model, as addressed in [START_REF] Berman | Big segment small segment global optimization algorithm on networks[END_REF], in which customers perceive the facility attractiveness in terms of their distance to the facility, while the second model (Section 3) is new, called the Huff origindestination (OD) trip model. In the OD trip model the facility attraction is a function of the length of the shortest path from the origin to the destination through the facility. In Section 4 a branch and bound algorithm is designed to solve both problems and two different procedures to obtain bounds are presented. Computational results are given in Section 5, comparing the two bounding strategies implemented. Some concluding remarks are presented in Section 6.

Huff location model

In this model, the finite set V of vertices of the network represents users, asking for a certain service. Each user a ∈ V has demand ω a ≥ 0. Such demand is being patronized by different existing facilities, located at points x 1 , . . . , x r on the network, so that the demand captured by facility at x i from user a is inversely proportional to a positive nondecreasing function of the distance d(a, x i ) from the user at a to the facility at x i . In other words, the demand captured by the facility at x i from the user at a is given by

ω a 1/ϕ a (d(a, x i )) r j=1 1/ϕ a (d(a, x j ))
.

(

) 2
The usual choice for each ϕ a has the form ϕ a (d) = d λa . When λ a = 2 for all a ∈ V, we have the so-called gravitational model. A new firm is entering the market, by locating one facility at some point on the network. This perturbs market share, since the new facility at x will capture a demand from a ∈ V equal to ω a 1/ϕ a (d(a, x)) 1/ϕ a (d(a, x)) + r j=1 1/ϕ a (d(a, x j ))

.

Here ϕ a is assumed to be non-negative, non-decreasing and twice continuously differentiable in R + . Expression (3) must be carefully considered when ϕ a (0) = 0. Indeed, in such case, if some x j exists with x j = a ∈ V , then the full demand of a will be captured by x j . Hence, such a will not be taken into account, and thus we assume in what follows without loss of generality that x j / ∈ V, j = 1, . . . , r. The goal of the entering firm is the maximization of its market share. This is written as the following optimization problem:

max x∈[0,le], e∈E a∈V ω a 1/ϕ a (d(a, x)) 1/ϕ a (d(a, x)) + r j=1 1/ϕ a (d(a, x j )) . (4)
Defining for each a ∈ V the positive constant β a ,

β a = r j=1 1 ϕ a (d(a, x j )) , (5)
it follows that Problem (4) can be rewritten as

max x∈[0,le], e∈E F (x) (6)
with F defined as

F (x) = a∈V ω a 1 1 + β a ϕ a (d(a, x)) . (7)
This problem has been addressed in [START_REF] Berman | Big segment small segment global optimization algorithm on networks[END_REF], in which a branch and bound algorithm with Interval Analysis bounds is proposed.

Huff OD trip model

In this model, we have the set {{a, b}, a ∈ V, b ∈ V }, of origin-destination pairs. The trip from origin a to destination b will imply a demand ω ab > 0. Consumers in their way from origin a to destination b will stop at one facility; they choose among the r existing facilities, x 1 , . . . , x r , and the new facility at x as in the model described in Section 2: the demand from a user in his way from a to b captured by each facility is inversely proportional to a positive nondecreasing function of the length of the path from the origin to the destination via the facility.

In other words, the demand captured by x from users in their trip from origin a to destination b is given by

ω ab 1/ϕ ab (d(a, x) + d(x, b)) 1/ϕ ab (d(a, x) + d(x, b)) + r j=1 1/ϕ ab (d(a, x j ) + d(x j , b)) . ( 8
)
As in the other model, the goal of the entering firm is the maximization of its market share. This is written as the following optimization problem:

max x∈[0,le], e∈E a,b∈V ω ab 1/ϕ ab (d(a, x) + d(x, b)) 1/ϕ ab (d(a, x) + d(x, b)) + r j=1 1/ϕ ab (d(a, x j ) + d(x j , b)) . (9)
Defining for each a, b ∈ V the positive constant β ab ,

β ab = r j=1 1 ϕ ab (d(a, x j ) + d(x j , b)) ,
it follows that Problem (9) can be rewritten as

max x∈[0,le], e∈E F (x)
with F defined as

F (x) = a,b∈V ω ab 1 1 + β ab ϕ ab (d(a, x) + d(x, b)) . (10)
We see that both the location and OD trip models yield an optimization problem of the same form, namely, [START_REF] Blanquero | On covering methods for d.c. optimization[END_REF], with a rather similar function F , as given by ( 7) and [START_REF] Drezner | The big triangle small triangle method for the solution of non-convex facility location problems[END_REF] respectively. Both functions F can be written in the form

F (x) = δ∈∆ ω δ 1 1 + β δ ϕ δ (d δ (x)) . ( 11
)
where ∆ = V and, for δ ∈ ∆,

d δ (x) = d(a, x)
for the location model, and

∆ = {{a, b}, a, b ∈ V } and, for δ ∈ ∆, d δ (x) = d(a, x) + d(x, b)
for the OD trip model.

Solving the models 4.1 Multimodality

The optimization problems described in Section 2 and 3 are, in general, multimodal, and standard optimization methods get stuck at local optima. This can be seen in simple examples, even when the network is a segment and is illustrated in the following examples for the Huff location problem introduced in Section 2.

First, data for one problem with two users and r = 1 facility on a segment were randomly generated. The objective function F of such instance is plotted in Figure 1 (left). One can see that the problem is bimodal. 100 runs of a local search procedure starting with a random point were performed.

In Figure 1 (right) one can see the histogram of the objective values provided by the optimizer: below a 50% of the runs yielded the global optimum, whereas the remaining runs stopped at the local not globally optimal solution.

Another instance, with 500 user and 100 facilities was also generated. In Figure 2 the problem is shown to be multimodal, and just below a 14% of the runs solved with the optimizer yielded the global optimum.

In the right part of Figures 1-2 the x axis is normalized, so that 1 corresponds to the best found objective value, z * , and a bar at x ∈ [0, 1] indicates that an objective value x • z * was found.

Algorithm

We will use a branch and bound algorithm [START_REF] Drezner | The big triangle small triangle method for the solution of non-convex facility location problems[END_REF][START_REF] Hansen | The minisum and minimax location problems revisited[END_REF][START_REF] Plastria | GBSSS: The generalized big square small square method for planar single-facility location[END_REF] to solve the problems under consideration. We first outline the algorithm that finds an optimal solution within a relative accuracy of ε, and later give the details on the bounding process, which will exploit monotonicity properties or the fact that the objective function is DC on each edge, [START_REF] Blanquero | On covering methods for d.c. optimization[END_REF][START_REF] Horst | Dc programming: Overview[END_REF][START_REF] Tuy | DC Optimization: Theory, Methods and Algorithms, Handbook of Global Optimization[END_REF][START_REF] Tuy | A d.c. optimization method for single facility location problems[END_REF]. We remind the reader that a function h is DC if it can be written as h = h + -h -, where h + , h -are both convex; the expression h + -h -is called a DC decomposition of h. In the description of the algorithm, we have ∆ = V for the Huff location model and ∆ = {{a, b}, a, b ∈ V } for the Huff OD trip model.

• Phase 1: Initialization -Fix the required accuracy ε > 0.

-Set LB = 0.

-Compute the all-pairs distance matrix.

-Calculate β δ , ∀δ ∈ ∆.

-Set the list H of remaining segments as empty.

• Phase 2: Prepare the list of segments -Consider the edge e as segment with its nodes as the segment vertices.

-The value of the objective function is evaluated at the segment midpoint. If the value is greater than LB, then LB is updated to such value and the midpoint stored as incumbent.

-Calculate an upper bound for the segment e, U B(e).

-

In case U B(e) ≥ LB • (1 + ε) insert e into H.
• Phase 3: Branch and bound process.

Repeat as long as no stop is reached:

-Select from H the highest upper bound segment, U B max , with LB as ε-optimal value and the incumbent as an ε-optimal solution. If U B max ≤ LB • (1 + ε), stop the algorithm with LB as the solution.

-The highest upper bound segment, U B max , is selected for a split at its midpoint into two smaller segments.

-The value of the objective function at the midpoint of the two small segments is calculated. If any of these values is greater than LB, then LB is updated and all segments from H whose upper bound is lower than LB are discarded.

-An upper bound for each small segment is calculated.

-All segments whose upper bound is greater than LB • (1 + ε) are added to H.

In the algorithm above, the segment midpoint yielding the best upper bound is given as ε-optimal solution. Observe that the full list of segments in H contain all ε-optimal solutions, and can thus be used in a two-phase process, as suggested in the GBSSS algorithm, [START_REF] Plastria | GBSSS: The generalized big square small square method for planar single-facility location[END_REF].

Let us detail the algorithm steps previously outlined. To calculate the all-pairs distance matrix we use the Floyd algorithm [START_REF] Bertsekas | Network Optimization: Continuous and Discrete Models[END_REF], which uses the edge length matrix to build recursively the distance matrix.

The computation of the bounds requires more detail. The algorithm needs the calculation of an upper bound, U B(s), for each segment s. We present two procedures, one based on Interval Analysis, and the other in properties of DC functions.

Interval analysis bound

When the distance d δ (x) decreases, the market share as given by the objective function [START_REF] Fernández | Solving a Huff-like competitive location and design model for profit maximization in the plane[END_REF] increases. Hence we obtain an upper bound on the market share F (x) for any location x on a segment s = [x 0 , x 1 ] ⊂ e ∈ E by replacing d δ (x) by the lowest possible of these distances on the segment. Defining therefore such lowest distance for the location model as in [START_REF] Berman | Big segment small segment global optimization algorithm on networks[END_REF]

d δ (s) = min{d(a, x 0 ), d(a, x 1 )} for δ = a,
and for the OD trip model as

d δ (s) = min{d(a, x 0 ) + d(b, x 0 ), d(a, x 1 ) + d(b, x 1 )} for δ = {a, b}, we obtain the Interval Analysis bound U B IA (s) = δ∈∆ ω δ 1 1 + β δ ϕ δ (d δ (s)) (12)

DC bound

An upper bound obtained making use of the fact that the objective function is DC on each edge exploits the following properties:

Proposition 4.1 Let I ⊂ R be an interval. Let d : I → R be a concave function on I, and let g : R → R be DC, with a DC decomposition given by g(x) = g + (x) -g -(x), with both g + and g -non-increasing functions. Then, the function f :

I → R defined as f (x) = g(d(x)) is DC on I and a DC decomposition is given by f (x) = f + (x) -f -(x), where f + (x) = g + (d(x)) and f -(x) = g -(d(x)).

Proof

The proof follows directly from the fact that the compositions g + (d(x)) and g -(d(x)) are also convex functions.

Remark 4.1 Proposition 4.1 makes use of a function g which can be written as the difference of two convex functions, g + and g -. Since such convex functions are also non-increasing, it turns out that g belongs to a subclass of DC functions, namely DCM functions, as introduced in [START_REF] Blanquero | Continuous location problems and big triangle small triangle: constructing better bounds[END_REF], which are those functions expressed as the difference of two convex monotonic functions, as g + , g -are. See [START_REF] Blanquero | Continuous location problems and big triangle small triangle: constructing better bounds[END_REF] for further properties.

We are now in position to give a bound for F on an edge exploiting the fact that F is DC. Let d δ (x) be the concave function given in [START_REF]Data instances for arc routing problems[END_REF]. Assuming ϕ δ (d) = d λ δ , (11) can be rewritten as

F (x) = δ∈∆ ω δ 1 1 + β δ d λ δ (x) (13)
Let us define a simpler function:

g(t) = 1 1 + βt λ (14)
The following DC decomposition is known for function g, [START_REF] Bello | On minimax-regret Huff location models[END_REF]:

c = λ -1 (λ + 1)β 1/λ g + (t) = g(c) + g (c)(t -c) if t ≤ c g(t) if t > c g -(t) = g(c) + g (c)(t -c) -g(t) if t ≤ c 0 if t > c That means that g(t) = g + (t) -g -(t).
Applying Proposition 4.1 we have that a DC decomposition for F as defined in ( 13) is:

F + δ (x) = ω δ g + (d δ (x)) F - δ (x) = ω δ g -(d δ (x)) F (x) = δ∈∆ F + δ (x) - δ∈∆ F - δ (x) = δ∈∆ (F + δ (x) -F - δ (x)) (15)
To construct an upper bound, UB DC , first we obtain a convex minorant of F - δ (x) as in [START_REF] Blanquero | Continuous location problems and big triangle small triangle: constructing better bounds[END_REF]:

F - δ (x) ≥ F - δ (x 0 ) + ξ δ (x -x 0 ) F (x) ≤ δ∈∆ (F + δ (x) -F - δ (x 0 ) -ξ δ (x -x 0 )) for ξ δ ∈ ∂F - δ (x 0 ), (16)
where ∂F - δ (x 0 ) denotes the set of subgradients of F - δ at x 0 [START_REF] Tuy | DC Optimization: Theory, Methods and Algorithms, Handbook of Global Optimization[END_REF].

Define for each δ ∈ ∆ the function from [START_REF] Labbé | Location on Networks[END_REF]:

U (x) = δ∈∆ (F + δ (x) -F - δ (x 0 ) -ξ δ (x -x 0 )) for ξ δ ∈ ∂F - δ (x 0 ) (17)
An upper bound for a segment s is obtained:

U B DC (s) = max{U (v 1 ), U (v 2 )} v 1 , v 2 being vertices of s.

Computational results

The algorithm described in Section 4 was programmed in an Intel Fortran Compiler XE 12.0 and executed with an Intel Core i7 computer with 8.00 Gb of RAM memory at 2.8 Ghz. The solutions were found within an accuracy of 10 -10 . The Huff location model was first tested on the 55-node and 134-edge Swain's (1971) network [START_REF] Serra | Surviving in a competitive spatial market: The threshold capture model[END_REF][START_REF] Marianov | Location-allocation of multiple-server service centers with constrained queues or waiting times[END_REF], see appendix, with both bounding strategies.

Several instances of the problem were generated using different values for the number r of existing facilities, ranging from low-saturated markets (r = 10% of the number of edges of the network, |E|) to high-saturated markets (r ≤ 90%|E|). For each value of r, 10 different problems were solved. Each problem is obtained by randomly and independently generating the demands (each vertex of the network is assumed to have a demand uniformly distributed in the interval (0, 1)) and the location of the existing facilities. To generate the r facility locations, r edges are randomly chosen with replacement; on each edge, the facility location is generated following a uniform distribution. The results are shown in Table 1.

The percentage r of existing facilities is shown in the first column. Then, it is reported the minimum, maximum, mean and standard deviation (std) for the number of iterations, i.e., the number of executions of Phase 3, the branch and bound (B&B) list size and the CPU time. The branch and bound list size is the maximum size of the data structure used for storage reached during the algorithm execution.

In all cases, the best solution is found in less than 0.02 seconds for DC bounds, and 0.75 seconds for Interval Analysis bounds. It is remarkable that DC bounds lead to a very stable procedure, as can be seen for all values of r in memory requirements as well as computational time. On the other hand, when using Interval Analysis bounds one can see very extreme cases. There is always a huge difference between the minimum and maximum for the number of iterations, branch and bound list size and time. This means that in the ten runs that are solved for each value of r, Interval Analysis bounds are quite erratic for problems of the same difficulty. Therefore, we have an algorithm that when using DC bounds spends the same resources (computing time) than when using Interval Analysis bounds, but is much more stable and reliable.

Then the algorithm was tested for both models on a battery of test instances of larger dimensions (up to 1000 nodes and 3083 edges) to analyze the dependence of running time and memory requirements with respect to the size of the network and the number of facilities for both bounding strategies.

To attain this end, the models were tested on 43 test networks obtained from [START_REF]Data instances for arc routing problems[END_REF][START_REF] Corberán | A branch & cut algorithm for the windy general routing problem and special cases[END_REF][START_REF] Reinelt | Tsplib -a traveling salesman problem library[END_REF]. For r = 10%|E| and r = 90%|E| of existing facilities, 10 instances of the problem were solved. Results for the comparison between the two methods for obtaining bounds are shown in the following tables. Results for the Huff location model are shown in Tables 2 and3 comparing both bounding procedures applied in Phase 2 and Phase 3 of the algorithm. Results for the Huff OD trip model are shown in Tables 4 and5. Results from Table 1 show that DC bounds seem to be sharper than Interval Analysis bounds but more computationally expensive; therefore, it seems desirable to delete segments in the initialization (Phase 2 of the algorithm) using easy bounds (Interval Analysis), and then use DC bounds for the rest of the segments (Phase 3). Hence, let us note that for the OD trip model, Interval Analysis bounds are used always in Phase 2 of the algorithm, and then, in Phase 3, the two different bounding procedures are compared. This strategy was also tested for the Huff location model with no improvement and therefore not used. The name of the network, number of nodes and edges are reported in the first columns of Tables 2345. The remaining columns have the same meaning than those in Table 1. In Table 2, when dealing with r = 10%|E| facilities in the Huff location model, both bounding methods are comparable, but when dealing with a saturated market, Table 3, DC bounds clearly outperform Interval Analysis bounds. This is due to the simplicity of the Interval Analysis bounds and because DC bounds are sharper. One can see that, when using DC bounds, the branch and bound size and the number of iterations is much smaller than with the other bounds, which means that DC bounds are harder to calculate, but sharper.

DC bounds

For the Huff OD trip model, as the number of computations of the bounds increases, both methods become comparable, as there is an equilibrium between the computational cost of bounds and how sharp they are.

One can see that, although the OD trip model is, in terms of computing time, much harder to solve than the location model, the number of iterations is relatively smaller, while the size of the branch and bound tree is comparable for most networks. Note that the relevant increase of the computing time for the OD trip model with respect to the Huff location model is due to the complexity of the objective function from the model, which involves evaluating many more terms. In both cases all problems can be solved in reasonable time.

Concluding remarks

In this paper we have addressed two Huff models on a network, by considering that the users choose the facility according to the distance from their location (Huff location model, [START_REF] Berman | Big segment small segment global optimization algorithm on networks[END_REF]) or according to the length of the shortest path from the origin to the destination visiting the facility (Huff OD trip model), which is a new model. Since the problems are shown to be multimodal, in order to obtain a global optimum, a branch and bound algorithm based on two different bounding procedures, Interval Analysis and DC optimization, is proposed.

The computational experience reported shows that large networks can be successfully handled with both bounding procedures where the DC procedure seems to be more stable in both time and memory requirements.

Extensions of these problems to the multifacility case deserve further study.

APPENDIX

The Swain data set, [START_REF] Marianov | Location-allocation of multiple-server service centers with constrained queues or waiting times[END_REF][START_REF] Serra | Surviving in a competitive spatial market: The threshold capture model[END_REF] Initial

Figure 1 :

 1 Figure 1: users=2, facilities=1, λ a = 1 ∀a

Figure 2 :

 2 Figure 2: users=500, facilities=100, λ a = 20 ∀a

Table 1 :

 1 Results of the 55-Node and 134-Edge Swain's Network for the Huff location model

		r(%|E|)			iterations			B&B list size		time
			min	max	mean±	std min	max	mean±	std min max mean±	std
		10 111	196	144.40±	34.68	9	61	31.50±	17.25 0.00 0.02	0.00± 0.00
		20	80	193	140.10±	34.32	15	83	42.70±	22.76 0.00 0.00	0.00± 0.00
		30	51	209	117.80±	47.86	25	96	66.70±	22.42 0.00 0.00	0.00± 0.00
		40	50	256	176.40±	55.33	15	107	84.80±	30.51 0.00 0.02	0.00± 0.00
	DC	50	81	290	172.20±	67.40	74	126	107.50±	14.92 0.00 0.02	0.00± 0.00
		60 121	317	166.50±	60.54	72	121	110.10±	16.09 0.00 0.02	0.00± 0.00
		70	94	323	190.00±	64.79 113	131	123.80±	6.36 0.00 0.02	0.00± 0.00
		80 145	332	209.70±	58.43 117	133	123.70±	5.62 0.00 0.02	0.00± 0.01
		90 137	296	211.40±	55.51 111	129	121.40±	5.66 0.00 0.02	0.00± 0.00
		10 158	318	223.70±	50.81	26	120	64.20±	31.72 0.00 0.02	0.00± 0.00
		20 185 178915 18100.40± 56504.50	29 58170 5890.10± 18369.30 0.00 0.59	0.06± 0.19
		30 203	70911 17779.60± 25278.44	75 25465 6219.20±	8904.76 0.00 0.20	0.05± 0.07
		40 243 243922 24639.40± 77048.06	96 84425 8551.40± 26659.27 0.00 0.75	0.08± 0.24
	Interval	50 187	37594	9132.40± 13971.49 118 13306 3143.50±	4794.67 0.00 0.11	0.03± 0.05
		60 186	29111	7734.70± 12130.03 118	9927 2703.30±	4174.80 0.00 0.09	0.03± 0.04
		70 212	10748	3029.00±	4148.73 138	3489 1038.00±	1349.47 0.00 0.03	0.01± 0.02
		80 172	1950	866.60±	677.23 127	648	310.10±	192.99 0.00 0.02	0.00± 0.00
		90 224	22520	4243.00±	7225.49 147	7248 1442.20±	2338.18 0.00 0.09	0.02± 0.03

Table 2 :

 2 Results of test instances for the Huff location model with 10%|E| facilities

	Interval Analysis bounds

Table 3 :

 3 Results of test instances for the Huff location model with 90%|E| facilities

	13

Table 5 :

 5 Results of test instances for the Huff OD trip model with 90%|E| facilities

This work has been partially supported by IMUS (Mathematics Institute of University of Seville) and by projects MTM2012-36163, Ministerio de Economía y Competitividad, Spain, P11FQM7603, P10TIC6064 and FQM-329, Junta de Andalucía, Spain, all with EU ERF Funds.


Susana Guimaraes
email: susana.guimaraes@cg971.fr

Arnaud Lenoble
email: arnaud.lenoble@u-bordeaux.fr

Acquisition, preparation & specimen documentation

Most of Caribbeetek's specimens were obtained via cession from public structures that manage the Antillean natural environment:

Guadeloupe wildlife Office, Martinique Regional Natural Park, St Martin Nature Reserve, St Barthelemy Territorial Environment Agency.

Others were acquired from animal health centres, zoos, taxidermists, livestock farms, hunters, inter-institutional exchanges or dedicated collecting actions. All acquisitions of wild specimen are carried out with authorisations and under the agreement of the Wildlife French Office.

Contextual information includes the geographical origin and biology of the specimens. Tissue samples are sent to interested research institutions (MNHN, Inra) in order to contribute to a better knowledge of the Caribbean fauna.

Protocols specific to each faunal group have been adapted for the preparation of skeletons. The most common process is enzymatic digestion following the method of Sorensen et alii (2010).

Content

The collection includes 1243 specimens documenting over 300 species.

Birds are the most represented, followed by bats and fishes.

5% corresponds to 3D models derived from CT-scans of specimens preserved by others institutions.

This solution is used for documenting rare, endangered or extinct species.

Guadeloupe's wild fauna is the best represented, with representation rates above 60%.

This is particularly true for birds with 90 of the 118 native species present in the collection. Moreover, the Guadeloupe avifauna includes sixteen so-called "patrimonial species" that group together the endemic species of Guadeloupe with endemic species of the Eastern Caribbean. Fourteen of these species are represented in Caribbeetek.

Inter-museum specimen exchanges and 3D models have allowed the collection to be extended to representatives of other Caribbean islands.

The collection can serve as a basis for studies of vertebrate fauna carried out in all of the Lesser Antilles.

Specimens documenting intraspecific variability are used to support university studies dedicated to osteological or paleoecological examinations of Caribbean terrestrial fauna.

 Caribbeetek is a public collection, freely available to any researcher who requests access for research purpose;

 All virtual models are available on demand, free of charge, and without counterpart;

 Inter-institutional exchanges and specimen loans are possible.

Aim

A well-known limitation of Caribbean zooarchaeology and paleontology is the difficulty of making taxonomic determinations on the basis of partial or fragmentary skeletal elements from taxa whose osteology is either or not described and that are frequently absent from most comparative collections.

Caribbeetek was created in 2012 to support paleozoological studies in Guadeloupe (Bivaag project and then ECSIT project). It initially focused on terrestrial fauna (chiropterans, birds, reptiles and amphibians). The collection was recently extended to fishes to aid archeozoological studies of pre-Columbian assemblages. It should soon be enriched with turtles and marine mammals.

Location

The collection is housed in two institutions:

-the Edgar Clerc Museum stores specimens covering almost all of the wild terrestrial vertebrate biodiversity of the island (108 species in collection), thus faciliting studies of vertebrate assemblages in Guadeloupe;

-the PACEA laboratory at the University of Bordeaux hosts specimens documenting the insular vertebrate biodiversity but also the morphological and ecomorphological variability of particular taxa (bats, rodents).

Caribbeetek

The Edgar Clerc Departmental Museum (Guadeloupe Museum of Archaeology) and the PACEA laboratory (CNRS-University of Bordeaux-MCC) have collaborated to develop a collection of comparative animal anatomy focusing on Caribbean wild fauna.

Caribbeetek is named in honour of Frederik A. Ober who carried out the first major naturalist collections in the Lesser Antilles during the second half of the 19th century and described his travels in several books, including the famous "Camp in the Caribbees: the adventures of a naturalist in the Lesser Antilles".

View publication stats View publication stats

  Birds in GuadeloupeS a l t a t o r a l b i c o l l i s O r t h o r h y n c h u s c r i s t a t u s E u l a m p i s j u g u l a r i s E u l a m p i s h o l o s e r i c e u s M e l a n e r p e s h e r m i n i e r i T u r d u s l h e r m i n i e r i C i n c l o c e r t h i a r u fi c a u d a M a r g a r o p s f u s c a t u s A l l e n i a f u s c a M y i a r c h u s o b e r i C o n t o p u s l a r o s t r i s E l a e n i a m a r n i c a C o e r e b a fl a v e o l a L o x i g i l l a n o c s G e o t r y g o of the endemic Guadeloupe Woodpecker (Melanerpes herminieri) Geographical origin of the specimens (physical and virtual). Colors -green: > 100 specimens, yellow: > 50 sp., blue: > 10 sp, pink: < 10 sp. 3D model of a Antillomys rayi mandible from Guadeloupe (length of the specimen: 3.2 cm) Caribbeetek A public skeleton collection documenting the wild F.W.I. fauna

This poster was prepared as a part of the ECSIT Project: « Écosystèmes insulaires tropicaux, réponse de la faune indigène terrestre de Guadeloupe à 6 000 ans d'anthropisa on du milieu », conducted by the CNRS with financial support from the European PO-FEDER program (grant n° 2016-FED-503), the Guadeloupe Regional Council, the DAC of Guadeloupe, and the CNRS (SEEG ECSIT).

28e CONGRÉS