A Liakopoulos

A Foka

S Vourli

L Zerva

F Tsiapara

E Protonotariou

Z Dailiana

M Economou

E Papoutsidou

C Koutsia-Carouzou

E D Anastassiou

E Diza

E Zintzaras

I Spiliopoulou

Dr E Petinaki
email: petinaki@med.uth.gr

Aminoglycoside-resistant staphylococci in Greece: prevalence and resistance mechanisms

Keywords: aminoglycosides, S. aureus, S. epidermidis, Greece

To determine the prevalence of aminoglycoside-modifying enzymes (AMEs) genes and the rate of resistance to amikacin, gentamicin, kanamycin and tobramycin a total of 1228 Staphylococcus aureus and 1004 S. epidermidis collected during 2009 from five Greek hospitals located in different geographic areas (Αthens, Nothern and South-western Greece) were tested. Among S. aureus isolates, 48.2% carried any AMEs gene; 73.3% of the positive strains carried the aph(3')-IIIa gene, 13.7% the ant(4')-Ia and 13% the aac(6')- gene. The rate of resistance as determined by disk diffusion method and Etest was 48.2% for kanamycin, 48.2% for amikacin, 12.9% for tobramycin and 6.3% for gentamicin. MRSA were more resistant than MSSA (92.7% resistant to amikacin and kanamycin, 25% to tobramycin and 12.3% to gentamicin versus 2.6% to amikacin and kanamycin and 0.3% to tobramycin respectively). Vitek 2 P549 card failed to characterize all S. aureus that carried the aph(3')-IIIa gene. Significant differences on the rate of resistance and distribution of AME genes were observed among S. aureus strains originating in different hospitals. On the other hand, 42.8% of S. epidermidis, all methicillin-resistant, carried any AMEs genes; 70% of the positive strains carried the aac(6')-Ie-aph(2") gene in combination with both ant(4')-Ia and aph(3')-IIIa genes, 18% carried the aac(6')-Ie-aph(2") gene with ant(4')-Ia, while, 12% carried the aac(6')-Ie-aph(2") gene with aph(3')-IIIa. The rate of resistance to all aminoglycosides tested was 42.8%; 100% concordance of susceptibility results was observed between those of disk diffusion method, Etest and Vitek 2 P549 card. No significant differences on the rate of resistance or the AME gene distribution among different geographic regions were observed among S. epidermidis strains. These results suggest an alarming rate of aminoglycoside-resistant staphylococciin Greece and the necessity for a continuous surveillance.

Introduction

Aminoglycosides are potent bactericidal agents that inhibit bacterial protein synthesis by binding to the 30S ribosomal subunit. They are often used in combination with either beta-lactam or glycopeptides, especially for the treatment of complicated staphylococcal infections, as these drugs act synergically [START_REF] Lowy | Staphylococcus aureus infections N[END_REF]. The main mechanism of aminoglycoside resistance in staphylococci is drug inactivation by cellular aminoglycoside-modifying enzymes (AMEs), which are categorized into four classes depending on the modifications they induce: acetyltransferases (AACs), phosphotransferases (APHs), nucleotidyltransferases (ANTS), and adenyltransferases (AADs). Among staphylococci, the most common AME is 6'-N-acetyltransferase-2"-O-phosphotransferase [AAC(6')-APH (2")], encoded by the aac(6')-Ie-aph(2") gene inactivating gentamicin, kanamycin, tobramycin, neomycin, and amikacin. Despite the high incidence of aminoglycoside resistance among methicillinresistant staphylococci as reported before 2000, there is currently little information on the incidence and predominant types of AMEs in many countries. This lack of information is generally of great importance, particularly for Greece where, the prevalence of methicillin-resistant Staphylococcus aureus (MRSA) and methicillinresistant coagulase-negative staphylococci (MRCoNS) is high in correlation that, the last years, aminoglycosides and specifically amikacin and gentamicin are widely used in our clinical settings.

The aim of the present study is to provide information regarding the prevalence of AMEs in correlation with the rate of aminoglycoside-resistance among staphylococcal isolates collected from five different Greek hospitals during 2009, in order to define a baseline for monitoring possible future increase.

Materials and Methods

Bacterial isolates and study design

Molecular methods

Detection of mecA, aac(6')-Ie-aph(2"), ant(4')-Ia and aph(3')-IIIa genes was performed by the multiplex PCR protocol described by Choi et al [START_REF]Performance Standards for Antimicrobial Susceptibility Testing: Nineteenth Informational Supplement M100-S19[END_REF]. Clones of representative MRSA and MRSE strains with different phenotypes and genotypes were defined by multi-locus sequence typing (MLST) ( http://www.mlst.net ).

Susceptibility testing

All isolates were tested for resistance to aminoglycosides by a disk diffusion method using amikacin 30 μg, gentamicin 10 μg, tobramycin 10 μg and kanamycin 30 μg, according to the Clinical and Laboratory Standards Institute guidelines [START_REF] Choi | Multiplex PCR for the detection of genes encoding aminoglycoside modifying enzymes and methicillin resistance among Staphylococcus species[END_REF]. Determination of MICs to the above aminoglycosides was done by Etest (AB Biodisk, Solna, Sweden), as of the manufacturer's instructions. All isolates were also tested for their susceptibility to various antimicrobial agents by the use of the VITEK 2 P549 card.

Statistical analysis

The effects of distribution of AME genes and regions, as well as their interactions were examined using a log-linear model, and the significance levels of the individual comparisons (post-hoc tests) were adjusted using the Bonferroni"s correction. The analysis was performed using GLIM3.77 (Royal Statistical Society).

Results

According to PCR results, 592 S. aureus (48.2%) carried any of the AMEs genes; 434 isolates, including 420 methicillin-resistant and 14 methicillin-sensitive, (73.3% of PCR-positive) carried the aph(3')-IIIa gene, 81 isolates, including 79 methicillinresistant and 2 methicillin-sensitive, (13.7% of PCR-positive) carried the ant(4')-Ia gene, whereas, 77 isolates, all methicillin-resistant, (13% of PCR-positive) carried the aac(6')-Ie-aph (2"). No coexistence of two or more AME genes was detected in any S. aureus strain. No direct correlation between sequence types (STs) and AME gene carriage was found. The great majority of ST80 and ST7 MRSA carried the aph(3')-IIIa, while strains of ST239 (the main nosocomial clone), ST377 and ST225, carried either the ant(4')-Ia gene or the aac(6')-Ie-aph(2") one. Strains of ST1 and ST30 did not carry any AME gene.

Susceptibility results based on disk diffusion method revealed that among 1228 S. aureus 636 isolates (51.8%) were aminoglycoside-sensitive, whereas, the remaining 592 isolates (48.2%) were resistant to any of the aminoglycosides tested; 48.2% were resistant to amikacin, 48.2% to kanamycin, 12.9% to tobramycin and 6.3% to gentamicin. Among MRSA tested, 92.7% exhibited resistance to amikacin, 92.7% to kanamycin, 25% to tobramycin and 12.3% to gentamicin. Among MSSA 2.6% exhibited resistance to amikacin, 2.6% to kanamycin and 0.3% to tobramycin.

None of the aminoglycoside-resistant strains lacked any of the AMEs genes, while, a strong correlation between the resistant phenotypes and presence of AMEs genes was observed (100% concordance). Results by Etest correlated well with those obtained by the disk diffusion method. However, discrepancies on susceptibility testing results of the automated VITEK 2 were observed: whereas the system identified correctly all S. aureus strains carrying the aac(6')-Ie-aph(2") or ant(4')-Ia gene, it failed to detect stains possessing only the aph(3')-IIIa (Table 2).

An interesting finding was the significant differences observed on the distribution of AME genes among participating hospitals (Table 1). Overall, there is a significant AME gene effect (P<0.01) in terms of distribution, whereas, the regions differ significantly (P<0.01) in terms of resistance. The interaction between the AME gene distribution effect and region effect is also significant (P<0.01). More Comparing the regions for each gene, South-western and Northern Greece are different from Athens (P<0.01) as far as aph(3')-IIIa, whereas, for South-western and Northern Greece are the same (P=0.99). For ant(4')-Ia, South-western Greece differ from Athens and Northern Greece (P<0.01), whereas, Northern Greece is not different from Athens (P=0.99). For aac(6')-Ie-aph(2"), there is a significant difference only between South-western and Northern Greece (P=0.02), whereas, Northern and Southwestern Greece are not different from Athens (P=0.26 and 0.59, respectively). Among 1004 S. epidermidis 430 MRSE (42.8%) carried at least two of AMEs genes. Among the positive isolates, 70% carried the aac(6')-Ie-aph(2") gene in combination with both ant(4')-Ia and aph(3')-IIIa genes, 18% carried the aac(6')-Ie-aph(2") gene with ant(4')-Ia, while, 12% carried the aac(6')-Ie-aph(2") gene with aph(3')-IIIa. Even though the majority of aminoglycoside-resistant MRSE belonged to ST2 and ST22, the main clones in Greece, polyclonality was observed among resistant genotypes without any correlation of clones to AME gene carriage.

Susceptibility results based on disk diffusion method revealed that 430 isolates (42.8%) were simultaneously resistant to all four aminoglycosides tested. All isolates were methicillin-resistant; no MSSE was detected with resistance to any aminoglycosides. Etest and VITEK 2 system identified correctly all aminoglycosideresistant isolates. No significant differences on the rate of resistance (P=0.54) or on AME gene distribution (P=0.29) among the regions were observed.

Discussion

Although aminoglycosides retain activity against the majority of Gram-negative clinical bacterial isolates in many parts of the world, the relatively frequent occurrence of nephrotoxicity and ototoxicity during treatment make physicians reluctant to use these compounds in everyday practice [START_REF] Grammatikos | Do we still need the aminoglycosides?[END_REF]. The doubts regarding the actual clinical utility of aminoglycosides resulted to the use of other antimicrobial agents with improved safety [START_REF] Karageorgopoulos | New antibiotics: optimal use in current clinical practice[END_REF]. However, the recent emergence of infections due to Gram-negative bacteria with advanced patterns of antimicrobial resistance has prompted physicians to reconsider these "old" antibacterial agents [START_REF] Giakoupi | Greek System for the Surveillance of Antimicrobial Resistance. (2009) KPC-2-producing Klebsiella pneumoniae infections in Greek hospitals are mainly due to a hyperepidemic clone[END_REF].

Aminoglycosides are also used in orthopaedic surgery, both locally and systemically, for prophylaxis of implant-related infections and treatment of skeletal infections. For local prophylaxis, application of gentamicin in the bone cement as well as coating of metallic implants with tobramycin or gentamicin appear to be effective, resulting in improved long-term implant survival [START_REF] Lucke | Gentamicin coating of metallic implants reduces implant-related osteomyelitis in rats[END_REF]. For the treatment of skeletal infections, gentamicin-loaded bone cement beads, tobramycin impregnated beads and bone graft substitutes offer high concentrations of antibiotics at the site of infection [START_REF] Rasyid | Concepts for increasing gentamicin release from handmade bone cement beads[END_REF]. However, local usage of aminoglycosides could rapidly induce antibiotic resistance among staphylococci with subsequent reduced effect on prophylaxis [START_REF] Anguita-Alonso | High rate of aminoglycoside resistance among staphylococci causing prosthetic joint infection[END_REF]. Since CoNS and S. aureus are the leading cause of device-related infections, continued surveillance concerning the resistance of staphylococci to aminoglycosides is required.

The rate of staphylococcal resistance to aminoglycosides varies from country to country and from year to year [START_REF] Shmitz | The prevalence of aminoglycoside resistance and corresponding resistance genes in clinical isolates of staphylococci from 19 European hospitals[END_REF][START_REF] Hope | Non-susceptibility trends among staphylococci from bacteraemias in the UK and Ireland, 2001-06[END_REF][START_REF] Emaneini | Phenotypic and genotypic evaluation of aminoglycoside resistance in clinical isolates of staphylococci in Tehran, Iran[END_REF][START_REF] Yadegar | Prevalence of the genes encoding aminoglycoside-modifying enzymes and methicillin resistance among clinical isolates of Staphylococcus aureus in Tehran, Iran[END_REF]. In a European multicenter study during 1999, Schmitz et al. reported that 23%, 29% and 31% of S. aureus were resistant to gentamicin, tobramycin and kanamycin, whereas, Hope et al. reported that 9% of MRSA and 2.5% of MSSA, collected from bacteraemia cases in UK and Ireland during 2001-2006, were resistant to gentamicin [START_REF] Shmitz | The prevalence of aminoglycoside resistance and corresponding resistance genes in clinical isolates of staphylococci from 19 European hospitals[END_REF][START_REF] Hope | Non-susceptibility trends among staphylococci from bacteraemias in the UK and Ireland, 2001-06[END_REF]. On the other hand, recent studies from Middle East showed a high incidence of aminoglycoside-resistant S. aureus [START_REF] Emaneini | Phenotypic and genotypic evaluation of aminoglycoside resistance in clinical isolates of staphylococci in Tehran, Iran[END_REF][START_REF] Yadegar | Prevalence of the genes encoding aminoglycoside-modifying enzymes and methicillin resistance among clinical isolates of Staphylococcus aureus in Tehran, Iran[END_REF]. In Greece, until now there was little information about the rates of aminoglycoside-resistant S. aureus. A previous study conducted in Central Greece during 2000 revealed that 11.2% and 18% of S. aureus exhibited resistance to gentamicin and tobramycin respectively [START_REF] Petinaki | Bacterial Resistance Study Group of Thessaly. (2001) Methicillinresistant Staphylococcus aureus in the hospitals of central Greece[END_REF]. Compared these results with our current data an inversion of the resistant phenotypes is observed, with a decrease of gentamicin and tobramycin-resistance (6.3% and 12.9% respectively) followed by a significant increase of kanamycin and amikacin resistance (48.2% for both agents versus 29.2% in 2000), partially due to the predominance of ST80 MRSA clone in the country [START_REF] Petinaki | Bacterial Resistance Study Group of Thessaly. (2001) Methicillinresistant Staphylococcus aureus in the hospitals of central Greece[END_REF][START_REF] Chini | Emergence of a new clone carrying Panton-Valentine[END_REF]. In addition, in contrast to other studies, the aph(3')-IIIa was the most prevalent AMEs gene (73.3%) in our strains [START_REF] Shmitz | The prevalence of aminoglycoside resistance and corresponding resistance genes in clinical isolates of staphylococci from 19 European hospitals[END_REF][START_REF] Emaneini | Phenotypic and genotypic evaluation of aminoglycoside resistance in clinical isolates of staphylococci in Tehran, Iran[END_REF][START_REF] Yadegar | Prevalence of the genes encoding aminoglycoside-modifying enzymes and methicillin resistance among clinical isolates of Staphylococcus aureus in Tehran, Iran[END_REF]. The significant differences on the resistance rate and AME gene distribution of S. aureus among Greek regions may be due to the predominance of specific clones and to aminoglycoside usage in each clinical setting combined with a possible horizontal AME gene transfer, as reflected in the results of participating hospitals.

Concerning the rate of aminoglycoside-resistant CoNS previous studies showed that were more often aminoglycoside-resistant than S. aureus [START_REF] Shmitz | The prevalence of aminoglycoside resistance and corresponding resistance genes in clinical isolates of staphylococci from 19 European hospitals[END_REF][START_REF] Hope | Non-susceptibility trends among staphylococci from bacteraemias in the UK and Ireland, 2001-06[END_REF]. In Greece, while 31% and 35% of S. epidermidis were resistant to gentamicin and tobramycin respectively during 2000, a significant increase of resistance (42.8% to both agents) is observed in S. epidermidis isolates eight years latter [START_REF] Petinaki | Bacterial Resistance Study Group of Thessaly. (2001) Methicillinresistant Staphylococcus aureus in the hospitals of central Greece[END_REF]16]. All isolates carried at least two AMEs genes. Indeed, CoNS have been frequently been reported to carry resistance determinants, including genes encoding AMEs and as such have the capacity to function as genetic reservoir. Conjugal transfer of resistance determinants between S. aureus and S. epidermidis leads to rapid dissemination of these determinants in the hospital environment.

As previously mentioned, VITEK 2 failed to detect kanamycin and amikacin resistance among S. aureus isolates carrying the aph(3')-IIIa gene. This can be explained by the fact that kanamycin and amikacin, that are accurate indicators of aph(3')-IIIa gene presence, were not included in P549 card; unfortunately, the newer P490 card also does not include them. Misidentification of such strains may lead to a therapeutic failure in cases where kanamycin or amikacin are used. On the contrary, VITEK 2 correctly identified aminoglycoside-resistant phenotypes among S. epidermidis, since all resistant S. epidermidis carried at least the aac(6')-Ie-aph(2") gene, and thus, exhibited resistance to gentamicin and tobramycin that are included in the P549 card. Thus, we recommend that laboratories using VITEK should employ alternative susceptibility testing methods for the detection of kanamycin and amikacin resistance in S. aureus.

In conclusion, this is the first study that describes the distribution of AME genes in correlation with the rate of aminoglycoside-resistant staphylococci in Greece during one year period (2009), providing useful comparative data for future studies. Today, according to hospital pharmacy records, amikacin and gentamicin are widely used in our clinical settings, specifically in those, where the prevalence of KPC-producer gentamicin-sensitive K. pneumoniae is high [START_REF] Giakoupi | Greek System for the Surveillance of Antimicrobial Resistance. (2009) KPC-2-producing Klebsiella pneumoniae infections in Greek hospitals are mainly due to a hyperepidemic clone[END_REF]. Probably, the re-introduction of these aminoglycosides in the clinical practice may enhance the emergence and dissemination of staphylococcal clones resistant to these agents. The high rate of aminoglycoside-resistant staphylococci (48.2% S. aureus and 42.8% S. epidermidis) in Greek hospitals located at different areas of the country indicates the need of accurate detection of resistance combined with knowledge of the epidemiology of 13 such resistance determinants, in order to limit further increase and to preserve the usefulness of aminoglycosides in the treatment of complicated infections.

  4'-Oadenyltransferase I [ANT(4')-I] is encoded by ant(4')-Ia gene and inactivates kanamycin, neomycin, tobramycin and amikacin. Finally, 3'-O-phosphotransferase III [APH(3')-III] is encoded by aph(3')-IIIa and inactivates kanamycin and amikacin [2].

  During 2009, a total of 1228 Staphylococcus aureus [621 methicilin-resistant (MRSA) and 607 methicillin-sensitive (MSSA)] and 1004 Staphylococcus epidermidis [720 methicillin-resistant (MRSE) and 284 methicillin-sensitive (MSSE), representing approximately 80% of coagulase-negative staphylococci], were tested for their susceptibility to aminoglycosides. Isolates were recovered from clinically significant specimens (blood, pus, pleural fluid, etc) in five Greek tertiary care hospitals located in different areas of the country (Athens, Northern, South-western Greece). A single isolate per patient was included in the study. Identification of staphylococci was performed by the Vitek 2 Advanced Expert System (bioMerieux, Marcy l'Etoile, France). S. aureus ATCC 25923 and S. epidermidis ATCC 12228 were used as control strains for susceptibility testing. In addition, S. aureus BM 3002 [possessing ant(4')-Ia)] and Enterococcus faecalis BM 6217 [possessing aac(6')-Ie-aph(2") and aph(3')-IIIa ] served as positive controls for each resistance gene.

  specifically, AME gene distribution effect follows different pattern across regions or the region effect varies across AME gene distribution. In particular, in South-western Greece all AME gene comparisons were significant (P<0.01). In Northern Greece the distribution of aph(3')-IIIa gene was different from ant(4')-Ia gene and aac(6')-Ieaph(2") (P<0.01), whereas, the ant(4')-Ia gene was not different from aac(6')-Ieaph(2") (P=0.27). In Athens , aph(3')-IIIa was marginally different from aac(6')-Ieaph(2") (P=0.04), whereas, the remaining comparisons were not significant [for ant(4')-Ia vs aac(6')-Ie-aph(2"), P=0.30, and for aph(3')-IIIa vs ant(4')-Ia, P=0.29].

Table 1 .

 1 Distribution of AME genes among S. aureus isolates collected from three regions of Greece: South-western Greece (A), Northern Greece (B) and Athens (C).

	Region Total No		AMEs genes distribution	No of AME-
		of	aph(3')-IIIa	ant(4')-Ia	aac(6')-Ie-aph(2")	positive
		isolates	No (%)	No (%)	No (%)	isolates
						No (%)
	A	517	196 (38)	12 (2.3)	22 (4.3)	230 (44.5)
	B	577	219 (38)	57 (10)	46 (8)	322 (55.8)
	C	134	19 (14.2)	12 (8.9)	9 (6.7)	40 (29.9)
	Total 1228	434 (35.3)	81 (6.6)	77 (6.3)	592 (48.2)

Table 2 .

 2 Comparison of susceptibility results based on Disk Diffusion method (DD), Etest, Vitek 2 P549 card (VT) and sensitivity of each 1 method, towards the respective genotypes of S. aureus.

	2
	Genotypes		DD	Etest		VT
	Gene	Positive isolates	Resistant	Sensitivity	Resistant	Sensitivity	Resistant	Sensitivity
		(number)	(number)	(%)	(number)	(%)	(number)	(%)
	aph(3')-IIIa	434	434	100%	434	100%	0	-
	ant(4')-Ia	81	81	100%	81	100%	81	100%
	aac(6')-Ie-aph(2")	77	77	100%	77	100%	77	100%
	Total	592	592	100%	592	100%	158	26.7%
	3


Florian Huc

Cláudia Linhares Sales

Hervé Rivano

The Proportional Coloring Problem: Optimizing Buffers in Radio Mesh Networks 1

Keywords: edge coloring, proportional coloring, mesh networks, call scheduling

In this paper, we consider a new edge coloring problem to model call scheduling optimization issues in wireless mesh networks: the proportional coloring. It consists in finding a minimum cost edge coloring of a graph which preserves the proportion given by the weights associated to each of its edges. We show that deciding if a weighted graph admits a proportional coloring is pseudo-polynomial while determining its proportional chromatic index is NP-hard. We then give lower and upper bounds for this parameter that can be computed in pseudo-polynomial time. We finally identify a class of graphs and a class of weighted graphs for which the proportional chromatic index can be exactly determined.

Introduction

Given a weighted graph (G, w), where w is a non-negative weight function over E, several distinct coloring problems on G have been defined. In [START_REF] Carlson | Scheduling to minimize interaction cost[END_REF], one wants to color the vertices of G using a given number of colors while minimizing the sum of the weights of the edges whose extremities receive the same color. On the other hand, in [START_REF] Mcdiarmid | Discrete mathematics and radio channel assignment[END_REF], one wants to color the vertices of G so that for each edge uv, |c(v)c(u)| is at least the weight of uv, where c(u) and c(v) are the colors assigned to u and v. There is also the problem of finding a partition of the edges of G into matchings (colors), each one of weight equal to the maximum weight of its edges, so that the total weight of the partition is minimized [START_REF] Rendl | On the complexity of decomposing matrices arising in satellite communication[END_REF].

In this paper, we consider a new edge coloring problem: the proportional edge coloring. Given a weighted simple graph (G, w), where w is a function called weight function defined on the edges of G and taking value between 0 and 1 (i.e. w : E → [0, 1]), we want to find a proper coloring, assigning one or more colors to each edge, which preserves the proportion given by the weights associated to each edge. If such a coloring exists, we want to find one using a minimum number of colors. We proved that deciding if a weighted graph admits a proportional coloring is pseudo-polynomial while determining its proportional chromatic index is NP-hard. We then give a lower bound and an upper bound for this parameter, both polynomial time computable. Finally, we exhibit a class of graphs and a class of weighted graphs for which we can exactly determine the proportional chromatic index.

The remainder of this section is dedicated to introduce the terminology used in this paper as well as some classical results on the proper edge coloring problem and the fractional edge coloring problem which will be useful in what follows.

Every graph in this paper is simple, which means that it is finite, with no loops nor multiple edges, except when it is clearly stated. In this case, it will be called a multigraph. We denote a weighted graph by (G, w), where w : E → [0, 1] is the weight function. The maximum degree of G is denoted by ∆(G) (or simply ∆, if G is clear from the context). A proper edge coloring of G is an assignment of colors to the edges of G such that no adjacent edges have the same color. Observe that if G is properly edge colored, each set of edges which are assigned a same color induces a matching of G.

The classical edge coloring problem is to determine the edge chromatic index of a simple graph G, χ ′ (G), that is, the minimum integer k such that G admits a proper edge coloring using k colors.

Equivalently, the edge chromatic index can be stated as an integer linear program. Let M be the set of all the matchings M of G. A proper edge coloring of G is an assignment of a non-negative integer weight p to every matching M ∈ M, such that for every edge e ∈ E(G), M :e∈M p(M ) ≥ 1. Then, the edge chromatic index χ ′ (G) is the minimum value given by ∀M ∈M p(M ) over all proper edge colorings of G. In 1964, Vizing proved that χ ′ (G) is at most ∆ + 1 [START_REF] Vizing | On an estimate of a chromatic class of a p-graph[END_REF]. Since it is at least ∆, we can classify every graph: a graph is Class 1 if its edge chromatic index is ∆ and Class 2 otherwise. Surprisingly, deciding if a graph is Class 1 or 2 is hard [START_REF] Holyer | The NP-completeness of edge colourings[END_REF], even for cubic graphs [START_REF] Koreas | The NP-completeness of chromatic index in triangle free graphs with maximum vertex of degree 3[END_REF]. In contrast, Padberg and Rao [START_REF] Padberg | Odd minimum cutsets and b-matchings[END_REF] have shown that the fractional relaxation of the integer linear program stated before can be solved in polynomial time, despite the fact that it has an exponential number of variables. The optimal solution of the fractional relaxation is known as the fractional edge chromatic index of G, denoted by χ ′ * (G), whose formal definition follows. Definition 1.1(i) Let M be the set of all the matchings M of G. A fractional edge coloring is an assignment of a non-negative weight p to each matching M ∈ M, so that for each edge we have: M :e∈M p(M ) ≥ 1. Thus, χ ′ * (G) is the minimum ∀M ∈M p(M ) over all fractional edge colorings of G;

(ii) For any natural number k ≥ 1, let G k be the graph obtained from G by replacing each edge by k parallel edges. Then

χ ′ * (G) = min k≥1 χ ′ (G k ) k .
Observe that in an optimal solution of the relaxed linear program, the weights of the matchings may not be integers. Therefore, since for every edge e, M :e∈M p(M ) ≥ 1, it may force an edge to belong to several weighted matchings. In terms of coloring, it means that an edge could be assigned several colors. Alternatively, the fractional edge chromatic index of a graph G can also be seen as the problem of determine integers a and b, such that we can assign a colors to each edge of G of a set of b colors, respecting the requirement that there are no two adjacent edges with a same color and such that b a is minimum. Putting these observations together with Vizing's result we have that:

∆(G) ≤ χ ′ * (G) ≤ χ ′ (G) ≤ ∆(G) + 1
A characterization stating that bipartite graphs may be decomposed into perfect matchings, due to Frobenius [START_REF] Frobenius | Über zerlegbare Determinanten[END_REF], has as an immediate corollary: for every bipartite graph G, ∆(G) = χ ′ * (G) = χ ′ (G). In fact, Reed proved that this equality holds also for nearly bipartite graphs (graphs containing a vertex whose removal produces a bipartite graph) [START_REF] Reed | Edge coloring near bipartite graphs[END_REF].

The upper bound for the proportional chromatic index presented in this work is closely related to the fractional edge chromatic index. To the reader interested in fractional edge coloring theory, we recommend [START_REF] Scheinerman | Fractional Graph Theory: a rational approache to the theory of graphs[END_REF].

Proportional coloring

The proportional coloring problem is formally presented in Definition 2.1 and illustrated in Figure 1. Before giving it, we define P(S) to be the set of all subsets of a given set S.

Definition 2.1 [c-proportional coloring] Given a weighted graph (G, w), a cproportional coloring of (G, w) is a function C : E → P({1, . . . , c}) such that for all e ∈ E, we have: |C(e)| ≥ cw(e) and for all pairs of edges e = (uv), f = (vw) for some nodes u, v, w, we have C(e) ∩ C(f ) = ∅. We call the proportional chromatic index of G, χ ′ π (G, w), the minimum number c of colors for which a c-proportional coloring of G exists. If it does not exist, then χ ′ π (G, w) = ∞. By proportional coloring, we mean a c-proportional coloring for some c.

To link the proportional chromatic index to the fractional edge chromatic index, we will use a multigraph, the m-graph, constructed from a weighted graph (G, w) and a positive integer m, whose definition follows. Definition 2. 3 [mcd(w), mcd] Given a weighted graph (G, w), with w taking value in Q, we set mcd(w) as the minimum common denominator of all the values taken by w. If there are no ambiguities, we call it mcd. It can be equivalently defined as the lowest common multiple of the denominators of all the values taken by w.

Remark 2.4 Given a weighted graph (G, w) and a positive integer m. Suppose there exists a proper edge coloring of G m using a multiple of m colors, say km colors for some integer k, such that every edge receives the same number of colors. Then, this proper edge coloring can be easily transformed into a proportional coloring of (G, w) with km colors (cf Figure 1).

1 4 1 5 1 6 1 2 1 2 {2} {1} {4} {3} {1} {2} {3} {3} {3, 4} {2} {1} {1, 2} {1, 2, 3} {4, 5, 6} {1, 2} {3, 4} {5}
Fig. 1. A weighted C 5 , a proper edge coloring of (C 5 ) 4 using 1 × 4 colors, the same coloring transformed into an optimal proportional coloring of the weighted C 5 and a suboptimal proportional coloring of the weighted C 5 using 6 colors.

Figure 1 is an example of a weighted graph and its proportional coloring. The proportional coloring problem is divided into two subproblems: the first one consists in proving that, for some integer c, there exists a c-proportional coloring of (G, w). The second one is to determine the proportional chromatic index of (G, w).

We now give simple facts that enable us to say if a weighted graph has a finite or an infinite proportional chromatic index and eventually to bound it.

Fact 2.5 Let (G, w) be a weighted graph. i If there is a vertex u of G with uv∈E w(uv) > 1, then χ ′ π (G, w) = ∞. ii If for all uv ∈ E, w(uv) ≤ 1/(∆ + 1), then χ ′ π (G, w) ≤ ∆ + 1. iii Similarly, if for all uv ∈ E, w(uv) ≤ 1/χ ′ * (G mcd ), then χ ′ π (G, w) ≤ q where q is the numerator of χ ′ * (G mcd ).
Proof.

i If a vertex u is such that uv∈E(G) w(uv) > 1, any coloring satisfying the weights have to repeat at least one color at the edges incident to u. Hence, no proportional coloring exists, since in a proportional coloring any two adjacent edges have to receive different colors. Therefore,

χ ′ π (G, w) = ∞. ii Recall that G is ∆ + 1-colorable ([16]). If for all uv ∈ E(G), w(uv) ≤
1/(∆+1), then any proper edge coloring using ∆+1 colors is a proportional coloring of (G, w).

iii If for all uv ∈ E(G), w(uv) ≤ 1/χ ′ * (G mcd ) = p/q, then a proper edge coloring assigning p colors out of q to each edge of G gives a proportional coloring of (G, w). Such a proper edge coloring exists by definition of χ ′ * (G mcd ). ✷

However the conditions given in Facts 2.5 are not sufficient. This is illustrated by Figure 2. Before moving to the next section, we prove a theorem which shows that we can suppose the value of the weight function of our instance (G, w) to be rational. Indeed, any weight function for which a proportional coloring exists can be replaced by a weight function taking only rational values and that accepts the same optimal proportional coloring. However, this theorem is only an existence theorem. Given a weight function for which we do not know an optimal proportional coloring, it does not state how to compute a rational weight function accepting the same optimal proportional coloring. Theorem 2.6 Given a weighted graph (G, w), w not necessarily taking values in Q, there is a weight function w ′ taking values in Q, such that for every e ∈ E, w(e) ≤ w ′ (e) and χ ′ . It is clear that any optimal proportional coloring of (G, w) is also an optimal proportional coloring of (G, w ′ ) and vice versa.

π(G,w) = χ ′ π(G,w ′ ) Proof. If χ ′ π (G, w) = ∞,
✷

From now on, we suppose that all the weight functions will take rational values.

A practical motivation from networking

The proportional coloring problem is a resource sharing problem similar in this sense to the one proposed in [START_REF] Barbosa | Sharing Resources at Nonuniform Access Rates[END_REF][START_REF] Yeh | Resource-sharing system scheduling and circular chromatic number[END_REF]. It is motivated by the following telecom-munication problem. Wireless mesh networks (WMNs) are cost-effective solutions for ubiquitous high-speed services [START_REF] Akyildiz | Wireless mesh networks: a survey[END_REF]. They are self-organized networks with a fixed infrastructure of interconnected wireless mesh routers whose purpose is to provide Internet access to mobile network users. This infrastructure, forming a wireless backhaul network, is connected to the internet by special routers called mesh gateways, as depicted in Figure 3. Each device has a single interface allowing to send or receive packets on the only available channel. This channel is therefore shared between all the nodes [START_REF] Liu | Optimal scheduling for link assignment in trafficsensitive STDMA wireless ad-hoc networks[END_REF]. We denote by call two antennas communicating together. We suppose that a call is achievable in both directions. Therefore, the network topology defines a graph G = (V, E) whose vertices are the routers and edges are the achievable calls.

At any time, an eligible state of the network implies that a set of simultaneously achievable calls are activated. Such sets strongly depend on the radio propagation and interference models. A quite broadly used family of models, denoted as binary models, assumes that two achievable calls are either interfering or not. In these settings, an eligible set of calls is a set of pairwise non interfering calls. In this application, we assume that there is a radio mesh network connecting routers through directional antennas. Hence, radio interferences and near-far effects impose that each router is involved in at most one call at a time. In consequence, a set of simultaneously achievable calls induces a matching on the graph.

In the following, the network is assumed to be synchronous and operating a slotted time division multiplexing (TDM). During each time slot, a set of simultaneously achievable calls is activated. We focus on the steady state of the network, which is therefore periodic, with a period consisting in a sequence of T time slots.

Therefore, the classical time slot assignment problem consists, for a given set of calls to be achieved periodically, in decomposing this set into a minimum number of subsets of simultaneously achievable calls. The number of subsets used is the size of the period, which is a key factor of the quality of service provided to the users of the network [START_REF] Molle | An optimization framework for the joint routing and scheduling in wireless mesh networks[END_REF]. This problem is equivalent to the proper edge coloring problem [START_REF] Bibb Cain | A link scheduling and ad hoc networking approach using directional antennas[END_REF].

The proportional coloring problem arises when we consider Constant Bit Rates (CBR) requests. In these settings, we are given a set of communication requests, each request being a source-destination path in the network, and a bit rate. Sending this amount of data on the paths induces that each call has to be periodically activated a given proportion of the time. This is modeled by a weight function w, which is computed in the following way: for each link e, one has to sum, over all connections whose paths use edge e, their bit rates, and divide the total by the capacity of the link.

The problem is now, if possible, to find a periodical schedule of the calls satisfying the CBR requests.

An assignment of calls to the time slots satisfies the CBR requests if:

• each slot is assigned an eligible set of calls;

• each call e is assigned to a number of slots k(e) such that k(e) = w(e) T , where T is the number of time slots.

After one time slot, each node creates a number of packets equals to its request and places them in its router queue, while each destination extracts the packets that have reached its router queue.

Observe that during each period, all the packets present at the beginning of the period make exactly one hop. Consequently, a packet will take as many periods to reach its destination as there are hops in its route.

Therefore, the total number of packets stored in the router queues of the nodes is proportional to the length of the period.

However, the size of these queues induces a cost in terms of memory to install in the routers as well as QoS parameters, like end-to-end delay. Therefore, decreasing the cost of the network and improving the QoS depends on minimizing the length of the longest queue, i.e., the length of the period. To illustrate this situation, we give an example in Figure 4.

Figure 4 represents four devices positioned on a line. Nodes s 1 , s 2 and s 3 send one packet at destination of d every five units of time. The routing we want to use is the shortest path (and unique path in our example) routing. It implies that a call s 1 s 2 is made one fifth of the time, a call s 2 s 3 two fifths of the time and a call s 3 d three fifths of the time. For this set of CBR requests and this routing, we can use the following period: node s 1 and s 3 send the packet generated during previous period (one slot), then node s 2 sends the packet generated during the previous period plus the one received (two slots).

Then node s 3 sends the two packets received during the previous period (two slots). This period of length five is optimal in the sense that no period can length less. In this case, we need two buffers, one of size one at s 2 and one of size two at s 3 .

Our problem is therefore to find a proper edge coloring of G that preserves the proportions given by the weights of the edges, with a minimum number of colors, which will be equivalent the optimum length of the period.

We finish this section with the remark that the assumption "w taking value in Q", that we did before as a consequence of Theorem 2.6, makes sense because of the origin of our practical problem. In fact, any instance of the original telecommunication problem would induce rational weights over the edges since they represent a ratio of a number of time slots over the length of the period. As for an example the weights could be computed from an initial time slot assignment that we want to improve, hence yielding rational weights with, more precisely, rather small operands.

Complexity results

In this section, we prove Theorem 4.1 which states the complexity results. ii) Determining the proportional chromatic index of (G, w) is NP-hard.

Proof.

i) Let (G, w) be any weighted graph and G mcd be its mcd-graph. Observe that G mcd can be computed in pseudo-polynomial time. Therefore, we will prove this item by showing that this problem is equivalent to determining the fractional edge chromatic index of G mcd , which can also be done in pseudo-polynomial time [START_REF] Scheinerman | Fractional Graph Theory: a rational approache to the theory of graphs[END_REF]. In fact, we claim that

χ ′ * (G mcd ) ≤ mcd ⇔ χ ′ π (G, w) < ∞. Indeed, suppose that χ ′ π (G, w) = k < ∞.
Then there exists a proportional coloring using k colors, such that each edge receives at least the correct proportion of colors. If we reproduce this proportional coloring mcd times (renaming colors appropriately), we still have a proportional coloring (which may not be minimal). The number of colors assigned to uv ∈ E is at least k × mcd × w(uv). This proportional coloring gives a proper edge coloring of G mcd using k × mcd colors and each edge receives at least k colors. Hence,

χ ′ * (G mcd ) ≤ k×mcd k .
On the other hand, suppose that χ ′ * (G mcd ) = q p ≤ mcd and consider an optimal fractional edge coloring of G mcd using q × c colors for some c (such a coloring is computed in pseudo-polynomial time during the process of computing the fractional edge chromatic index [START_REF] Scheinerman | Fractional Graph Theory: a rational approache to the theory of graphs[END_REF]). This fractional edge coloring gives at least p × c colors to each edge. Furthermore, this fractional edge coloring can be extended to an edge coloring of G that uses q × c colors and gives p × c × mcd × w(e) to each edge e (c.f. Remark 2.4). This proper edge coloring is a proportional coloring, since, by assumption,

χ ′ * (G mcd ) ≤ mcd ⇔ q p ≤ mcd ⇔ mcd × w(e) × p × c ≥ q × c × w(e). In consequence, χ ′ π (G, w) = k ≤ q × c < ∞.
ii) Let G be any graph of maximum degree ∆. Now, set the weights of all the edges to 1 ∆+1 . Since G is (∆ + 1)-colorable, G is proportionally (∆ + 1)colorable. But now observe that if (G, w) admits a proper edge coloring with ∆ colors, then it admits a proportional coloring with ∆ colors. Since determining if G is ∆ or (∆ + 1)-colorable is NP-hard [START_REF] Holyer | The NP-completeness of edge colourings[END_REF], it is NP-hard to determine the proportional chromatic index of an instance (G, w). ✷

Partial characterization

Despite the difficulty of computing the proportional chromatic index of a weighted graph (G, w), we can deduce polynomially computable lower and upper bounds. Clearly, for a weighted graph (G, w), χ ′ π (G, w) ≥ ∆. We can get a better lower bound by trying to satisfy the proportions at every vertex of G. If no positive integer m satisfies all the above equations, then no proportional coloring exists. Otherwise, if there is a solution m and (G, w) admits a proportional coloring, then χ ′ π (G, w) is at least m. Proof. If there is no integer m satisfying all the Equations (1), then for every integer m, there is a vertex u such that uv∈E ⌈mw(uv)⌉ > m. In particular, this holds when m = 1 which from Fact 2.5 i), implies that no proportional coloring using m colors exists.

On the other hand, if a proportional coloring exists which uses c colors, then c needs to satisfy all the Equations (1). Therefore, m ≤ c. ✷

If (G, w) admits no proportional coloring, its proportional chromatic index is infinite. However, if it admits a proportional coloring, we can deduce an upper bound from the fractional edge coloring of the corresponding mcd-graph. Corollary 5.2 (Upper bound) Let (G, w) be a weighted graph. If a proportional coloring of (G, w) exists, then a proportional coloring of (G, w) using q colors exists, where q is the number of colors used in a fractional edge coloring giving p colors to each edge with q p = χ ′ * (G mcd(w) ) and gcd(p, q) = 1.

Proof. Corollary 5.2 is a consequence of the proof of Theorem 4.1 i). Indeed, to prove the fact that deciding if a proportional coloring of G exists is pseudopolynomial, we construct a proportional coloring of G using c×q colors from a fractional edge coloring of G mcd(w) using c × q colors and where q is the upper bound of the theorem. This construction is valid for any value of c ≥ 1. ✷

In general, given a weighted graph (G, w), the mcd of the values taken by w is not an upper bound, even if the graph has a finite proportional chromatic index. Indeed, consider (P, 1 3 ): the Petersen graph P with a weight function uniformly equal to 1 3 on all its edges, as depicted in Figure 5. Since χ ′ * (P ) = 3, we have χ ′ π (P, 1 3 ) < ∞. However χ ′ π (P, 1 3 ) > 3 = mcd( 1 3 ), since P does not have any proper edge coloring using only 3 colors.

The proportional chromatic index may differ from the upper bound given by the fractional edge chromatic index, since an edge can receive proportionally more colors than what it asks for. This is illustrated by the example of Figure 1 (C 5 , w). Indeed, in this example, we have mcd(w) = 60, and so in (C 5 ) mcd(w) , an edge with weight 1/2 is replaced by 30 parallel edges.

1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3
Fig. 5. The Petersen graph with the weight function uniformly equal to 1 3 on all edges.

Hence, in (C 5 ) mcd(w) there is a vertex whose degree is 60, which implies that χ ′ * ((C 5 ) mcd(w) ) ≥ 60 > 4 = χ ′ π (C 5 , w). This can also be illustrated by the triangle of Figure 2, or by any odd cycle with constant weight function w=1/2, as no proportional colouring of these weighted graphs using two colours exist ; two being the mcd of the values taken by w. However these examples have infinite proportional chromatic index, which is not the case with the weighted Petersen graph given above.

One can also wonder on which classes of graphs the proportional coloring problem can be solved in polynomial time. The next theorems announce two positive results.

Theorem 5.3 Let (G, w) be a weighted graph. If there is a solution to the set of Equations [START_REF] Akyildiz | Wireless mesh networks: a survey[END_REF], for m the minimum solution, we have:

χ ′ π (G, w) = m ⇔ χ ′ (G m ) = m. Proof. We start to prove χ ′ π (G, w) = m ⇒ χ ′ (G m ) = m.
Observe that a proportional coloring using m colors is a proper edge coloring of G, which gives at least ⌈m × w(e)⌉ colors to each edge e. Additionally, in G m , each edge e of G is replaced by ⌈m × w(e)⌉ edges, so the previous proper edge coloring gives a proper edge coloring of G m . By this, we proved that χ ′ (G m ) ≤ m. It remains to prove that χ ′ (G m ) ≥ m. Suppose χ ′ (G m ) < m. Then, by the definition of a proper edge coloring, we would have, for each vertex v,

v st uv∈E ⌈m×w(uv)⌉ ≤ χ ′ (G m ) ≤ m-1. This is stronger than uv∈E ⌈(m- 1) × w(uv)⌉ ≤ χ ′ (G m ) ≤ m -1, which contradicts the hypothesis of the minimality of m. Hence χ ′ (G m ) = m. We now prove χ ′ π (G, w) = m ⇐ χ ′ (G m ) = m.
If G m has a proper edge coloring using m colors, it gives a proportional coloring of (G, w) using m colors. Hence χ ′ π (G, w) ≤ m. This fact together with Theorem 5.1 implies that χ ′ π (G, w) = m. ✷

This theorem gives the proportional chromatic index of some simple classes of weighted graphs. A sufficient condition satisfied by these classes is that, given an integer m, the corresponding m-graph is of Class 1. It includes trees, grids and more generally any bipartite graph as stated by Corollary 5.5.

We now give a remark that will be useful in what follows.

Remark 5.4 Let (G, w) be a weighted bipartite graph and m a solution to the set of Equations ( 1). Then m is finite and G m is well defined. Observe that for every u ∈ G m and for edge e incident to u, in G m , e is replaced by ⌈m × w(e)⌉ edges. Because u satisfies Equation ( 1), the total number of edges incident to u is less than or equal to m. Therefore, the maximum degree of G m is at most m. In addition, observe that in the case where w is uniform and (G, w) satisfies the necessary condition of Fact 2.5, then any m ≥ ∆(G) is a solution of the set of Equations (1).

Corollary 5.5 Let (G, w) be a weighted bipartite graph. If there is a solution m to the set of Equations ( 1), then (G, w) accepts a proportional coloring using m colors.

Proof. Let (G, w) be a weighted bipartite graph and m a solution to the set of Equations ( 1), which does exist by assumption. By Remark 5.4, the maximum degree of G m is at most m. Observe that G m is also bipartite and hence its edge chromatic index is equal to m, that is, χ ′ (G m ) = m. The corollary follows as a consequence of the previous theorem. ✷

We would like to know if it is possible to extend this corollary to other classes of graphs. After trees and bipartite graphs, a natural class of graphs to investigate is the class of outerplanar graphs. A graph is outerplanar if it admits a planar representation where all the vertices are incident to the infinite face.

If G is an outerplanar other than an odd cycle and w is an uniform function on the edges of G, Theorem 5.3 also includes (G, w). Indeed, we know that a simple outerplanar graph is of class 1 if and only if it is not an odd cycle. Hence an outerplanar graph which is not an odd cycle and such that all the edges have the same multiplicity has its chromatic index equals to its maximum degree. Therefore, given a weighted outerplanar graph (G, w) satisfying the necessary condition stated in Fact 2.5, with w being a constant, we have that, for any m ≥ ∆(G), G m is an outerplanar multigraph with all the edges having the same multiplicity and with maximum degree equal to m (see Remark 5.4. In consequence, χ ′ (G m ) = m. Hence, m is a solution to the set of Equations (1) and Theorem 5.3 applies.

We observe that the requirement of w uniform over all the edges is not a sufficient condition, since in the proof of Theorem 4.1, we used a weighted graph with a constant weight function. Furthermore, if G is an outerplanar but w is not uniform on all its edges, then we can not apply Theorem 5.3. To see this, consider a simple outerplanar graph G containing a triangle, a weight function w and m, the minimum solution to the Equations ( 1). If w is such that the m-graph G m has too many multiple edges on one of its triangle, even if its maximum degree is at most m (see Remark 5.4), it may not have edge chromatic index equal to m. Figure 6 illustrates this situation, by presenting an outerplanar weighted graph G and its outerplanar multigraph G m , m = 12. In this example, we have a triangle with edge multiplicity 4, which forces the edge chromatic index to be at least 12 while its maximum degree is 10. We can also give more information on the proportional chromatic index of a weighted graph, when its weight function verifies some properties. The next theorem goes in this direction: Theorem 5.6 Let (G, w) be a weighted graph. Let e be an edge with an end vertex v such that uv∈E w(uv) = 1. Then, the denominator of w(e) divides the number of colors used in a proportional coloring of (G, w). In particular, it divides χ ′ π (G, w).

Proof. Consider a proportional coloring of (G, w) using m colors. Let e be an edge with an end vertex v such that uv∈E w(uv) = 1. We have that uv∈E ⌈m×w(uv)⌉ ≤ m = uv∈E m×w(uv). Hence ⌈m×w(e)⌉ = m×w(e). In particular, it is true for m = χ ′ π (G, w). ✷

In a weighted graph (G, w), a vertex v with uv∈E w(uv) = 1 represents a saturated node, i.e., a node through which no extra request can be routed. A weighted graph (G, w) such that every edge is incident to a saturated vertex represents a saturated network, i.e. a network in which no extra request can be routed.

The existence of a saturated node leads to a non-approximability result: Corollary 5.7 Let (G, w) be a weighted graph and v a vertex of G with uv∈E w(uv) = 1. We call q the mcd of all the w(uv) for uv ∈ E. χ ′ π (G, w) cannot be approximated within the constant q -1.

Proof. Theorem 5.6 implies that q divides the number of colors used in a proportional coloring. Hence no proportional coloring exists with χ ′ π (G, w) + i colors for 0 < i < q. ✷

Conclusion

Motivated by the applications modelled by the proportional coloring problem and its hardness, we let the following general open questions: find approximation algorithms for classes of graphs which usually occur in telecommunication, as circular-arc graphs and triangular lattices, and determine other classes of graphs (and weighted classes of graphs) for which the proportional chromatic index can be computed in polynomial time.

Definition 2 . 2

 22 [m-graph] Given a weighted graph (G, w) and a positive integer m, the m-graph G m is obtained from G by replacing each edge e = (uv) ∈ E(G) by ⌈mw(e)⌉ multiple edges uv.

Fig. 2 .

 2 Fig. 2. Example of a weighted graph satisfying ∀v, uv∈E(G) w(uv) ≤ 1, which has infinite proportional chromatic index.

  just take w ′ as the constant weight function equals to 1 on every edge. Suppose χ ′ π (G, w) = c. Given an edge e ∈ E, we set w ′ (e) = ⌈c×w(e)⌉ c

Fig. 3 .

 3 Fig. 3. A WMN topology: mobile clients access Internet through a multi-hop wireless backhaul network of routers and gateways.

Fig. 4 .

 4 Fig. 4. An example of constant requests on a network composed of 4 antennas

Theorem 4 . 1

 41 Let (G, w) be a weighted graph on n vertices. i) Determining if there exists a proportional coloring for (G, w) is pseudopolynomial in max(n, mcd).

Theorem 5 . 1 (

 51 Lower bound) Let m be the minimum positive integer sat-isfying for all u ∈ V : uv∈E ⌈mw(uv)⌉ ≤ m (1)

Fig. 6 .

 6 Fig. 6. Example of an outerplanar weighted graph (G, w) and its G 12 multigraph where χ ′ (G 12 ) = ∆(G 12 ).

This work has been partially funded by European project IST/FET AEOLUS and ICT-2008-224460 (WISEBED), and French projects ANR-JC OSERA, ARC CARMA and ANR VERSO ECOSCells.

[START_REF] Barbosa | Sharing Resources at Nonuniform Access Rates[END_REF] Contact author. Email: Florian.Huc@epfl.ch

.[START_REF] Bibb Cain | A link scheduling and ad hoc networking approach using directional antennas[END_REF] This work has been partially done while C. Linhares-Sales was visiting the Mascotte project-team with support by CNPq/Brazil.


Francesco Di Maio

Enrico Zio
email: enrico.zio@polimi.it

Tao Liu

Jiejuan Tong

PASSIVE SYSTEM ACCIDENT SCENARIO ANALYSIS BY SIMULATION

Keywords: Passive Safety System, Residual Heat Removal system (RHRs), High Temperature Gas-Cooled Reactor -Pebble Bed Modular (HTR-PM), Scenario Analysis, Monte Carlo Simulation

In this paper, a simulation framework of analysis is presented aiming at evaluating the safety performance of the Residual Heat Removal system (RHRs) of the Chinese High Temperature Gas-Cooled Reactor -Pebble Bed Modular (HTR-PM) under uncertain operation conditions, and components and equipments failures. A transparent and fast model of the passive system has been implemented in MATLAB to reproduce the three-interconnected natural circulation trains of the RHRs, for removing the residual heat of the reactor core after a reactor shut-down. The model is characterized by a one-dimensional mono-phase moving fluid, whose operation is based on thermal-hydraulic (T-H) principles. The model is coded into a Monte Carlo (MC) failure engine for sampling single and multiple components faults at random times and of random magnitudes.

Accidental transients of the system are simulated, highlighting equipment contribution to system failure.

INTRODUCTION

Generation IV reactor systems (e.g., High Temperature Gas Cooled Reactors (HTGR)) are resorting to the use of passive safety systems, in combination with active safety or operational systems. As passive systems do not need external input, especially energy, to operate [IAEA, 1991], they are expected to contribute significantly to the improvement of plant safety [START_REF] Mathews | Functional reliability analysis of safety grade decay heat removal system of Indian 500MWe PFBR[END_REF]. To achieve this, the design of passive systems must be supported by proper modeling and simulation.

In this paper, we consider the Residual Heat Removal system (RHRs) of the Chinese High Temperature Reactor-Pebble Modular (HTR-PM) [START_REF] Zhengy | Characteristics of the 250MW Pebble-Bed Modular High Temperature Gas-Cooled Reactor in Depressurized Loss of Coolant Accidents[END_REF]. The function of this system is to remove the residual heat of the reactor core after shut-down. A thermal-hydraulic (T-H) model for scenario analysis has been implemented with characteristics of transparency, for a clear interpretation of the system response, and speed of calculation, to allow for extensive scenario analysis.

Figure 1 sketches the elements of scenario analysis, starting from the collection of system information and the characteristics of the uncertainties affecting the performance of the system function, in terms of plant setting parameters (e.g., power rate), physical condition parameters (e.g., temperatures and pressures), and structural and material properties (e.g., friction and roughness of pipes). Consideration of uncertainties is fundamental for passive safety systems, because the magnitude of the natural forces which drive their operation is small, so that it could be easier for the system behavior to be influenced by the uncertainties, and its robustness poorer than for active system [START_REF] Marquès | Methodology for the reliability evaluation of a passive system and its integration into a probabilistic safety assessment[END_REF][START_REF] Pagani | The impact of uncertainties on the performance of passive systems[END_REF]. In the end, the design must guarantee that the risk of accidents leading to the dissatisfaction of some of the basic safety functions BSFs (i.e., reactivity control, residual heat removal, primary pressure control and containment release) is not increased.

Figure 1 System scenario analysis

The model of the RHRs describes a one-dimensional mono-phase moving fluid, based on thermal-hydraulic (T-H) principles [START_REF] Zhengy | Characteristics of the 250MW Pebble-Bed Modular High Temperature Gas-Cooled Reactor in Depressurized Loss of Coolant Accidents[END_REF]. A preliminary sensitivity analysis has allowed identifying the most influential parameters affecting the passive system function. These parameters have been the focus of the analysis of the system behavior in major accidents, e.g., pipe blockage, water pipe rupture and RHR trains isolation. The model is embedded in a MC scheme for sampling the uncertain operation condition and single or multiple failures, for an integrated analysis of system response.

The paper organization is as follows. In Section 2, the main system characteristics of the High Temperature Reactor-Pebble Modular (HTR-PM) are briefly introduced. The results of a variance decomposition sensitivity analysis [McKay, 1996] aimed at identifying the most influential parameters for the system function are summarized in Section 3 [Yu et al., 2010a]. Equipment malfunctions and deficiencies are identified by HAZOP [Zio, 2007] in Section 4, and representative accidental scenarios are simulated, as described in Section 5. Finally, some conclusions are drawn in Section 6.

THE RESIDUAL HEAT REMOVAL SYSTEM OF THE HTR-PM

High-temperature gas-cooled reactors have developed for nearly 50 years. Today's Chinese design of the High Temperature Gas-Cooled Reactor-Pebble bed Modular (HTR-PM) is based on the technology and experiences of the HTR-10 10MW high-temperature gas-cooled test reactor (HTR-10) designed in China in 2000.

In Figure 2, a sketch of the HTR-PM layout is shown. At a first glance, the HTR-PM design has the following key technical features [START_REF] Zhengy | Characteristics of the 250MW Pebble-Bed Modular High Temperature Gas-Cooled Reactor in Depressurized Loss of Coolant Accidents[END_REF]:

• Uranium dioxide (UO 2 ) fuel kernel coated by tri-isotropic (TRISO) ceramics such as pyrolytic carbon and silicon carbide (SiC), in order to retain fission products in the particle under a fuel cladding temperature of 1600ºC in accident cases.

• A one-zone core design is implemented, consisting of approximately 420,000 spherical fuel elements in a pebble-bed with a diameter of 3m and an average height of 11m, that are charged and discharged in a so-called "multi-pass" mode, which means that before the fuel elements reach the discharge burn-up, they go through the reactor core several times.

• Decay heat in the fuel elements is dissipated by means of heat conduction and radiation to the outside of the reactor pressure vessel, and then taken away to the ultimate heat sink, i.e., the passive RHRs, by water cooling panels on the surface of the primary concrete cell.

Therefore, no coolant flow through the reactor core is necessary for decay heat removal in case of loss of coolant flow or loss of pressure accidents. T .

The T-H code models the following steps of the process:

1. The residual heat radiates from the reactor vessel and other thermal sources to the water in the water-cooling wall;

2. Because of the difference in temperature, natural convection initiates through water, in the water-cooling wall and pipes connected with the air-cooling heat exchanger; then, heat is transferred to the water side of the heat exchanger;

3. The heat is transferred by thermal conduction from the water side to the air side of the heat exchanger, due to the difference of temperature;

4. As the air-cooling heat exchanger is located in the air-cooling tower, natural convection of air sets up and takes heat to the final heat trap-atmosphere.

Table 1 reports the 37 parameters identified as most influential for the system function, by expert judgment. To account for the associated uncertainties, these are described by probability distributions defined on the basis of previous experience and/or information obtained by skilled experts. Air temperature and pressure variability are inferred by historical data ,collected by a representative Chinese Automatic Weather Station (CAWS) in different months. For example, a typical weather trend for a northern locality of china can be shown in Figures 4 and5.

N Parameter

SENSITIVITY ANALYSIS OF THE PASSIVE RHRs MODEL

In a previous work, the variance decomposition sensitivity analysis method has been adopted to assess the impact of uncertainties on the variance of the outlet water temperature , w out T [Yu et al., 2010a]. The results, supported by an application of the Analytical Hierarchical Process [Yu et al., 2010b], reveal that the parameters most influential on the uncertainty of , T , respectively, whereas W is the heat which is absorbed by the system per unit of time: therefore, the smaller is the number of functioning trains, the lower is the overall cooling capacity of the system, the more important is the total heat to be absorbed.

The sensitivity analysis has proceeded to considering the effects of varying simultaneously all plant setting parameters (i.e., power W , the number of pipes in the air cooler a N , the number of water cooling pipes for each train w N , etc.), physical condition parameters (i.e., water pressure in the pipes w P , temperature of air in the air-cooled tower , a in T , inlet pressure of the air in the air-cooler tower , a in P , etc.), and structural and material properties (i.e., thermal resistance of pipes inside of the heat exchanger i R , thermal resistance due to the dirty of the pipes fins o R , etc.). As suspected from the results of the individual inputs, the group of physical configuration parameters is more sensitive than either the group of plant setting parameters or the group of structural and material properties.

Coherently, the following analyses are focused on the inlet air temperature in the air-cooled tower , a in T , the water pressure in the pipes w P and the plant setting parameters such as w N , a N and , f in A : all these will be sampled for simulating the system response under normal conditions and accidental scenarios.

EQUIPMENT FAILURE CHARACTERIZATION

For the analysis of interest, four deviations which may lead to the functional failure of residual heat removal are:

-Higher water temperature in the water cooling wall -Higher water temperature in the inlet of the air cooler -Lower air cooling flow in the air cooler -Lower pressure in the primary cooling pipes An engineering analysis of these situation lead to the identification of four main types of equipment failures that can run the system into severe situations.

-Reduction of the number of air cooling pipes (N a ), i.e, pipe blockage -Reduction of the number of water cooling wall pipes (N w ), i.e, water pipe rupture -The inlet shutters fails partially close, so that A f,in is smaller -One of the three trains is isolated for some reason, e.g. maintenance

SIMULATION OF ACCIDENTAL SCENARIOS

In this Section, the outcomes of the sensitivity and engineering analyses are exploited to define the simulation of a number of accidental scenarios. The calculations have been run with the system at power 1200 W =

[kW] and water pressure in the pipes 0.3 w P =

[bar], under sampled conditions of the variables initial inlet air temperature and initial air pressure in the air-cooled tower, , a in T and , a in P , respectively [Yu et al, 2010a;Yu et al., 2010b].

A Monte Carlo (MC) sampling procedure has been developed for injecting equipment faults at random times and of random magnitudes in the RHRs. The set of faults considered are:

• A random number of air pipes in the air cooler (10) for train 1/2/3 fails stuck/broken at time t 1/2/3 in [1,2880] [min], so that the total number of pipes in each air cooler is equal to 1/ 2 / 3 a N .

• The inlet shutter (12) of the air cooling tower (11) for train 1/2/3 fails stuck at time t 4/5/6 in [1,2880] [min], providing a corresponding air flow crossing area 1/ 2 / 3 , f in A in each tower.

• A random number of water cooling pipes (3) for train 1/2/3 fails stuck/broken at time t 7/8/9 in [1,2880] [min], so that the total number of pipes in each air cooler 1/ 2 / 3 w N .

• One of the three trains is unavailable at time t 10 in [1,2880] [min].

For exemplification, we comment on two of the scenarios simulated. In accidental transient 1, the set of failures in T evaluated in nominal conditions is also reported (dashed-dotted line).

At time 240 [min], a sudden unavailability of a number of air pipes in the first train causes an increase of the temperature of the coolant in all the trains, that have to endorse the load of more heat removal proportionally to their capabilities: in fact, the temperature of trains 2 and 3 increases much more than for train 1, whose heat capacity has degraded due to the water pipe rupture. At time 1140

[min], the inlet shutter of train 3 fails partially closed: , The maximum value of ,

w out
T is reached in train 2, because all its components and equipments are available along all the system life, so that less resistances are present in this train allowing for a larger heat exchange. In this accidental transient, the system is still able to react to components and equipments failures without exceeding the threshold c T . In accidental scenario 2, , w out T exceeds c T because the components and equipments states are such to limit the heat transfer capacity, so that the heat decay is larger than the heat releasable to the heat-trap atmosphere. The set of failures characteristic of this scenario are listed in

Month

Figure 2

 2 Figure 2 Sketch of HTR-PM

Figure 3

 3 Figure 3 Schematics of 1 train of the RHRs in the HTR-PM [Zhao et al., 2008]

Figure 4 Figure 5

 45 Figure 4 Seasonal variation of the mean air temperature during one day in different months

  W , the air inlet temperature of air in the air-cooled tower , a in T and the water pressure in the pipes w P . This occurs because , a in T and w P directly influence the air cooling capacity and the critical temperature

  c

  , because of its abrupt decrease in heat capacity, i.e., pipe blockage. Shortly after, at time 1200 [min], train 3 goes out of service. Subsequently, the water temperature in train 3 drops to , a in T , reaching the temperature of the cold sink of the air-water exchanger in the air cooler tower, whereas a sharp increase of , w out T is recorded in trains 1 and 2.

  Figure 6 Evolution of the outlet water temperature in RHR train 1.

  Figure 9 Evolution of the outlet water temperature in RHR train 1.

Table 1 Parameters relevant for the behavior of the passive RHRs.

 1

			Distribution	Note
	1	W	/	Residual heat power
	2	Ta,in	Seasonal/Daily	Temperature of inlet air in the air cooled tower
	3	xi1	Uniform	Resistance coefficient of elbow
	4	xi2	Uniform	Resistance coefficient of header channel
	5	xiw	Uniform	Resistance coefficient of the water tank walls
	6	xia,in	Uniform	Sum of the resistance coefficients of inlet shutter and air cooling tower and silk net
	7	xia,out	Uniform	Sum of the resistance coefficients of outlet shutter and air cooling tower and silk net
	8	xia,narrow	Uniform	Resistance coefficient of the narrowest part of the tower
	9	Pa,in	Seasonal/Daily	Pressure of the inlet air in the cooler tower
	10	dx	Uniform	Roughness of pipes
	11	Ha	Normal	Height of chimney
	12	La	Normal	Length of pipes in the exchanger
	13	Na	Normal	Total number of pipes in the air cooler
	14	Af	Normal	Air flow crossing are in the narrowest part of the tower
	15	Af,in	Normal	Inlet air flow crossing area in the tower
	16	Af,out	Normal	Outlet air flow crossing area from the tower
	17	Af,narrow	Normal	Crossing area in the narrowest part of the tower
	18	S1	Normal	Distance between centers of adjacent pipes in horizontal direction
	19	S2	Normal	Distance between centers of adjacent pipes in vertical direction
	20	S	Normal	Distance between fins in the ribbed pipe
	21	Da	Normal	Pipes inner diameter in the air cooling exchanger
	22	Do	Normal	Pipes outer diameter
	23	Douter	Normal	Rib outer diameter

  Table 2 have been randomly sampled to occur. The ,

	w out T	evolution during the
	life time of 48 [h] (2880 [min]) is plotted in Figures 6-8 for train 1, 2 and 3, respectively; for
	comparison, the , w out

Table 2 Failures in the accidental sequence 1

 2

						April
	Simulation start time			22.00
	Time of occurrence	[min]	Type of failure	Magnitude of failure
	t1	240	N	1 a	200
	t2	/	N	2 a	/
	t3	/	N	3 a	/
	t4	/	1 f in , A	/
	t5	/	A	2 f in ,	/
	t6	1140	3 f in , A	8.53
	t7	/	N	1 w	/
	t8	/	N	2 w	/
	t9	/	N	3 a	/	70	0	500	1000	1500	2000	2500	3000
	t10	1200	Failed train	Yes

Table 3

 3

	. The

Table 3 Failures in the accidental sequence 2

 3

ACKOWLEDGEMENTS

This work was supported in part by the Science & Technology Fellowship Programme in China financed by the European Commission (EuropeAid/127024/L/ACT/CN) and by the Important National Science & Technology Specific Projects in China (No. 2008ZX06902-009).


Francesco Di Maio
email: f.dimaio@mail.polimi.it

Selina S Y Ng

Kwok Leung Tsui

Enrico Zio
email: enrico.zio@ecp.frenrico.zio@supelec.frenrico.zio@polimi.it

NAÏVE BAYSIAN CLASSIFIER FOR ON-LINE REMAINING USEFUL LIFE PREDICTION OF DEGRADING BEARINGS

Keywords: Degradation, Prognostics, Condition Monitoring, Residual Useful Life, Naïve Bayesian Classifier, Thrust Ball Bearings

In this paper, the estimation of the Residual Useful Life (RUL) of degraded thrust ball bearings is made resorting to a data-driven stochastic approach that relies on an iterative Naïve Bayesian Classifier (NBC) for regression task. NBC is a simple stochastic classifier based on applying Bayes' theorem for posterior estimate updating. Indeed, the implemented iterative procedure allows for updating the RUL estimation based on new information collected by sensors located on the degrading bearing, and is suitable for an on-line monitoring of the component health status. The feasibility of the approach is shown with respect to real world vibration-based degradation data.

Introduction

Bearings degradation and the prediction of their Residual Useful Life (RUL) are of great significance for many production and manufacturing industries. Indeed, according to [START_REF] Nandi | Condition Monitoring and Fault Diagnosis of Electrical Motors -A Review‖[END_REF], almost 40-50% of all rotating machinery failures are bearing related, thus greatly affecting productivity in those sectors in which rotating machinery plays a relevant role [START_REF] Gebraeel | Residual Life Predictions in the absence of prior degradation Knowledge‖[END_REF].

Maintenance of these components is currently performed under periodic schemes of inspections, which are not capable of giving due account to their actual health status. On the other hand, experience tells us that thrust bearing failures without warning can result in even catastrophic consequences [START_REF] Kelly | Maintenance and its management‖[END_REF][START_REF] Singh | Induction machine drive condition monitoring and diagnostic research-a survey‖[END_REF]. This has motivated a movement, in recent years, towards approaches for monitoring and predicting the trend of degradation of rotating machinery components within a philosophy of Condition-Based Maintenance (CBM) [START_REF] Williams | Condition-based maintenance and machine diagnostics‖[END_REF][START_REF] Mann | Statistical-based or condition-based preventive maintenance‖[END_REF][START_REF] Marseguerra | Condition-based maintenance optimization by means of genetic algorithms and Monte Carlo simulation‖[END_REF]. CBM allows pre-scheduled maintenance based on the information on the components RUL, reduces inventory costs of spare parts, reduces the risk of catastrophic failure and eliminates unexpected outage [START_REF] Yam | Intelligent predictive decision support system for condition-based maintenance‖[END_REF].

In practice, the estimation of the RUL may be difficult to obtain because the degradation state may not be directly observable and/or the vibration measurements may be affected by large noise and disturbances, as is the case of interest in our study.

From a methodological viewpoint, approaches for predicting rotating machinery failures can be grouped into two main categories [START_REF] Heng | Rotating machinery prognostics: State of the art, challenges and opportunities‖[END_REF]: traditional reliability approaches and condition-based approaches. Traditional reliability approaches are based on the distribution of event records of a population of identical units: their historical time-to-failure data are used to estimate the population characteristics, e.g., the mean time to failure. These methods are effective for units produced in high volumes, whereas they are not effective to characterize small numbers of component or even a single component currently running on a particular machine [START_REF] Heng | Rotating machinery prognostics: State of the art, challenges and opportunities‖[END_REF].

Condition-based approaches can be divided into two categories: physicsbased models and data-based models [START_REF] Chiang | Fault detection and diagnosis in industrial systems‖[END_REF]. Physics-based approaches attempt to set up comprehensive mathematical models to describe the physical phenomena underlying the components degradation process and failure modes, and to estimate their RUL. However, these models have been found inadequate in addressing the stochastic nature of defect-propagation; correspondingly, a number of stochastic models have been developed. For example, the Bayesian updating procedure for residual-life distributions estimation from component degradation signals proposed in [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF] assumes a first model in which the degradation signal exhibits independent identically distributed (IID) errors for an exponential signal trajectory, and a second model in which the error fluctuations follow a Brownian Motion process (BM); Bayesian updating is used to improve the estimation of the stochastic parameters in the exponential model, thus eventually improving the estimate of the true signal trajectory. Both methods proposed in [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF] are parametric and model-based.

Still, uncertainty due to assumptions and simplifications in the models may pose limitations to their applicability in practical industrial applications, where the failure-specific mechanistic knowledge is often hard to gather without interrupting operation [START_REF] Heng | Rotating machinery prognostics: State of the art, challenges and opportunities‖[END_REF].

In such cases, data-driven techniques may serve the purpose as they utilize monitored operational data related to system health and routinely collected from the machine, instead of building models based on comprehensive system physics and human expertise. They can be beneficial when understanding of first principles of system operation is not straightforward or when the system is so complex that developing an accurate model is prohibitively expensive.

In this work we resort to a non-parametric stochastic approach that relies on a Naïve Bayesian classifier (NBC) [START_REF] Duda | Bayes Decision Theory‖[END_REF] to benefit both from the stochastic framework approach and the data-driven characteristic of the method; NBC consists in a simple stochastic classifier (based on applying Bayes' theorem) which iteratively allows for the on-line estimation of the RUL of a degrading thrust bearing by updating the estimates of the posterior, i.e., the RUL, based on the data collected from the sensors located on the monitored component which are compared to other data stored in a database and previously collected on degraded components of the same type. Despite its simplicity, the NBC learning scheme performs well on most classification tasks, and is often significantly more accurate than more sophisticated methods [START_REF] Frank | Technical note: Naive Bayes for regression‖[END_REF] also when applied for regression.

The paper contents are structured as follows. Section 2 presents the benchmark data of degrading thrust bearings taken from [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF]. Section 3 contains the description of the approach at the basis of the RUL estimation, with an overview of the NBC framework. In Section 4, the results of the application of the approach here proposed are presented, and an evaluation of its performance is given. Section 5 shows a comparison of the results provided by the NBC-based approach and those given by the IID-and BM-based models [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF]. Finally, some conclusions on the advantages and limitations of the approach here propounded are given in Section 6.

Case Study

The case study here illustrated has been analyzed previously with other model-based methods [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF][START_REF] Gebraeel | Sensory-Updated Residual Life Distributions for Components with Exponential Degradation Patterns[END_REF][START_REF] Gebraeel | Residual Life Predictions in the absence of prior degradation Knowledge‖[END_REF]. The database consists of degradation signals of 25 thrust ball bearings obtained from a run-to-failure experiment under accelerated testing conditions. The bearings are run under constant operating conditions and vibration frequency spectra are acquired every 2 minutes. Thus, the monitored signal consists of the average of the first seven harmonics of the amplitude of the bearing rotation frequency. Generally, the first part of the readings remains low, representing the normal operation phase where there are no cracks on the bearings. Then, after the failure onset, the vibration readings increase as defects develop and worsen, until the bearings is considered to be in a failure state and cannot be used, when a failure threshold d=0.03 V RMS is reached [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF], as shown in Figure 1. Among the 25 available degrading trajectories, 2 never reach the failure threshold, so that they are discarded. The performance of the approach will be compared with that obtained in [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF], whose trajectories are shown in Figures 2345, when a threshold equal to 2.5 mV RMS is used to cut out the stable signal when the bearing is working in nominal conditions, i.e., in our analysis we will mainly concentrate on the second phase of the degradation.

Methodology

Figure 6 shows a schematic sketch of the computational framework; for ease of illustration, a single signal   ft is considered (a generalization to the case of multidimensional trajectories is straightforward). It is continuously monitored throughout the time horizon of observation On the other hand, it is assumed that N trajectories of evolution in time of values of relevant signals (reference patterns) are available from measurements collected in dynamic degradation scenarios of the bearings under analysis. These trajectories last all the way to system failure, i.e., to the time when the signal reaches the threshold value beyond which the bearing loses its functionality. For reasons which will become clear in the following, the database constituting the reference pattern library is organized in a reference matrix

  1 NJ R    , whose generic element   , r i j is the j-th segment of length 1 j  of the values of the i-th reference pattern, 1, 2,..., iN  , 1, 2,..., . jJ 
The idea underpinning the RUL estimation is to evaluate the similarity between the test trajectory pattern   ft of the developing degradation and the N reference trajectory patterns stored in the database through a NBC [START_REF] Frank | Technical note: Naive Bayes for regression‖[END_REF] and use the RULs of these latter to estimate, at each time j T , the RUL of the former, accounting for how similar they are. In other words, the NBC calculates the probabilities of each feature given each class based on the set of reference trajectory patterns and uses these probabilities to calculate the probabilities of a test trajectory pattern belonging to each class given its feature values.

For the scope of this work, the implementation of the NBC went through three different steps: data pre-processing, training of the classifier, testing of the classifier. Data pre-processing entails the discretization of the target range, because it has been demonstrated that by so doing NBC for regression performs comparably to well known methods for time series prediction [START_REF] Kononenko | Semi-naive Bayesian classifier‖[END_REF]. The pre-processing step, or data preparation, is a key step in the non-trivial Knowledge Discovery and Data Mining process upon which the success of the entire process depends [START_REF] Fayyad | Knowledge Discovery and Data Mining: Towards a Unifying Framework‖[END_REF][START_REF] Debuse | A Methodology for Knowledge Discovery: a KDD Roadmap‖[END_REF][START_REF] Cios | Trends in data mining and knowledge discovery‖[END_REF][START_REF] Crisp-Dm | Cross industry standard process for data mining[END_REF]. The second and the third steps are necessary within a supervised learning scheme, as it is the one proposed in this work.

a. Data pre-processing

When classifiers are used for regression, it is common to discretize the continuous target variable into a number of intervals [START_REF] Torgo | Regression using classification algorithms‖[END_REF][START_REF] Yang | Discretization for naive-Bayes learning: Managing discretization bias and variance‖[END_REF]. This decreases the number of classes to handle and increases class frequencies with the same number of training trajectories, so that to avoid over-fitting, improving the accuracy and the generalization capability of the learning scheme, and reducing at the same time the computational burden of the classification algorithms [START_REF] Yang | Discretization for naive-Bayes learning: Managing discretization bias and variance‖[END_REF].

The time-point T d at which each of the N=19 reference pattern hits the predetermined failure threshold d has been discretized by regression-bydiscretization strategy [START_REF] Kononenko | Semi-naive Bayesian classifier‖[END_REF]: the time-point T d has been discretized into C equal width intervals of 40 min, each one with a reference mean life T i , i=1,2,…,C, where C <=N, as the class label representative of the interval class label i, i=1,2,…C (Figure 7). To illustrate, discretization on this particular training and testing dataset (total 23 trajectories) reduces the number of class labels from 19 to 11.

b. Training of the NBC

A Bayesian classifier is a simple probabilistic classifier with strong (naïve) independence assumptions [START_REF] Duda | Bayes Decision Theory‖[END_REF]. Although these assumptions make the model too simple to be realistic, Naïve Bayesian Classifiers generate much better than expected performances in some context [START_REF] Frank | Technical note: Naive Bayes for regression‖[END_REF][START_REF] Sahami | A Bayesian approach to filtering junk e-mail‖[END_REF]. Moreover, NBC bear the advantage that they can be trained very efficiently in a supervised learning setting. By default, NBC handle nominal classes and assign the class with maximum probability to the instance.

Using Bayes' theorem, we can write the probability that a realization of feature variables 12 , ,..., J f f f , i.e. the vibration data at inspection times T 1 , T 2 ,…, T J , belong to class i as a conditional model:

        1 1 1
,..., ,..., ,...,

J J J p i p f f i p i f f f f f 
In practice we are only interested in the numerator of that fraction, since the denominator does not depend on i and the values of the features f j are given, j=1,2,…,J , so that the denominator is effectively constant. The numerator is equivalent to the joint probability model which can be written as follows, using repeated applications of the definition of conditional probability:

          1 1 2 1 1 2 1
,..., , ... , , ,...,

J J J p C f f p i p f i p f i f p f i f f f  
The "naïve" conditional independence assumes that each feature f j is conditionally independent of every other feature f t , ti  , so that:

    , j t j p f i f p f i 
Thus, the joint model can be expressed as:

      1 1 ,..., J Jj j p i f f p i p f i   
This means that under the above independence assumptions, the conditional distribution over the class variable i can be expressed like this:

      1 1 1 ,..., J Jj j p i f f p i p f i Z   
where Z (the evidence) is a scaling factor dependent only on 12 , ,...,

J f f f , i.e., a
constant if the values of the feature variables are known. We resort to the supervised NBC classifier implemented in the Weka package [START_REF] Hall | The WEKA Data Mining Software: An Update‖[END_REF] to perform the classification of the degrading trajectories of the thrust ball bearings presented in Section 2.

In this work, given the shortage of reference pattern, we proceed at training the NBC to provide   1 ,..., J p i f f with a 10-fold cross validation scheme [19] so that to improve the robustness of the predictions provided by the classifier. Reference patterns have been broken down into k=10 disjoint subsets of approximately equal size. Then, 10 experiments have been performed and for each experiment the k-th subset is removed: the system is trained on the remaining data.

In other words, NBC is trained with 10-fold cross validation based on the N

examples in   ,1 NJ R  .
Instead of testing on the held-out subset for each of the 10 folds, we evaluate the 10 NBC models on a separate set of m test trajectory patterns. Also, in contrast to the default NBC which chooses the nominal class with maximum probability, we utilize the probabilities of each of the possible class and the fact that the classes are continuous in our context. For each one of the 10 trained NBC, the estimated RUL is the weighted sum of all the class labels (time-to-failure) based on the probability of belonging to i-th class minus the time when the prediction is made, i.e.,  

1 ˆC i k j i j i RUL T p T T   
. Finally, the estimated remaining useful life is taken equal to be the mean value of the k=10 RUL predictions provided by the 10 trained NBC, i.e.,

   

Results

Figures 8-11 are plots of the actual RUL (middle dotted line) and the predicted RUL (solid lines) for each test trajectory degrading patterns shown in Figures 2345. Note that the RUL estimation is particularly stable and accurate throughout the life of degradation patterns 1, 3 and 4 (Figures 2, 4 and5, respectively). In Figures 891011, the ±10% error bounds are also give with respect to is the time at which the m-th test pattern hits the predetermined failure threshold d: the larger is the root mean square error, the smaller is the prediction accuracy. The root mean square error values are further scaled to the number of predictions for each test pattern for comparison with each other in Table 2. It is worth pointing out that the largest error occur in predicting the RUL of test pattern 2 (Figure 3): a further analysis has pointed out that this occurs because no reference patterns are stored in the database with a life equal to

Comparison with other methods

The results published for exponential Brownian Motion error exponential process (BM) and exponential stochastic process with Independent Identically Distributed error (IID) model in [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF] are compared with that of the NBC at times equal to 0,5, 0.75 and 0.9 of the bearing failure time T d . The best prediction among BM, IID and NBC for each of the four test bearings is highlighted in bold in Table 3. Note that these do not serve as direct comparisons because the data for training the NBC in this work is a subset, rather than an identical set, of that used for estimating the prior distributions in [START_REF] Gebraeel | Residual.-life distributions from component degradation signals: a Bayesian approach‖[END_REF].

It is important to note that NBC always generates much more useful predictions at 0.5T d for all test bearings rather than the other two methods. Generally, NBC performs better at early stage of the failure onset, i.e., RUL estimates are much more accurate than those provided by BM and IID models, and when training patterns with similar life-span exists, e.g., bearings 1, 3 and 4. When data is not available for training, either IID or BM models perform much better at earlier and later stage than NBC, as it occurs for bearing 2. This means that if training data is not available and good RUL predictions are only required at later stages e.g., the last hour before reaching the failure threshold, BM and IID have to be preferred. The use of both methods in parallel (hybridization of model-based and data-based approaches, i.e., BM and IID, and NBC, respectively) could provide even better estimation of RUL of thrust ball bearings. More work can be done in this path, and some efforts can be devoted to the reduction of the dimension of data set needed for training NBC while maintaining the capability of RUL prediction.

Conclusions

In this paper, a data-driven stochastic approach that relies on an iterative Naïve Bayesian Classifier (NBC) is presented as a feasible tool for estimating the Residual Useful Life (RUL) of degraded thrust ball bearings. NBC is fed with the vibration signal collected on the component and provide the RUL of the most similar degradation trajectory among a set of reference patterns that are stored in a database. Indeed, the simple Bayesian framework upon which NBCs are founded allows for an iterative updating procedure for updating the RUL estimation based on new information collected by sensors located on the degrading bearing. The results show that the approach is suitable for an on-line monitoring of component health status in a real industrial application.

Figure 1 .

 1 Figure 1. Typical behavior of degraded bearings measurements [8]

Figure 2 .

 2 Figure 2. Degradation pattern 1

Figure 4 .Figure 5 .

 45 Figure 4. Degradation pattern 3

Figure 6 .

 6 Figure 6. schematic sketch of the computational framework

Figure 7 .

 7 Figure 7. schematic sketch of the class discretization

  of the NBC For each m-th test trajectory, m=1,2,3,4, at each inspection time T j , j=1,2,…,J , the vector of readings   ft is fed to the 10 trained NBC as input vector. Readings beyond the inspection time are set to be missing and are ignored by NBC. NBC estimates the probability of   ft belonging to each class i, i=1,2,…,C, the single RUL estimate   ˆkj RUL T is provided and the average RUL from the 10 predictions is obtained as   ˆj RUL T .

dT

  , i.e., until d T , the time at which the state hits the predetermined failure threshold d. Results indicate that the proposed method can give an high informative prediction at a sufficiently early stage of the component or structure life, so as to allow for preventive maintenance actions.

Figure 8 . RUL estimation for the degrading pattern of Figure 2 Figure 9 . RUL estimation for the degrading pattern of Figure 3 Figure 10 . RUL estimation for the degrading pattern of Figure 4 Figure 11 . RUL estimation for the degrading pattern of Figure 5

 8293104115 Figure 8. RUL estimation for the degrading pattern of Figure 2

2 dT

 2 , whereas NBC performs much better with degradation patterns whose behavior is well covered by the reference patterns stored in the training dataset.

Table 2 . Evaluation results -root mean square prediction error

 2

Table 3 . Failure times evaluated at three prediction intervals using BM, iid, RVM and NB for four degrading patterns

 3

	Bearing	Actual	Prediction		IID			BM			NBC
		life	time
	1	230	T d	5%	50%	95%	5%	50% 95%	5%	Avg. 95%
		116	0.5 T d	360	732	>10000	209	411	1797	206	230	288
		174	0.75 T d	226	356	774	187	251	811	206	234	288
		208	0.9 T d	290	290	514	213	247	731	206	234	288
	2	508	T d	5%	50%	95%	5%	50% 95%	5%	Avg. 95%
		254	0.5 T d	284	562	>10000	270	346	1062	326	414	492
		382	0.75 T d	402	774	>10000	399	485	1461	406	408	408
		458	0.9 T d	486	1166 >10000	479	589	1861	406	418	448
	3	228	T d	5%	50%	95%	5%	50% 95%	5%	Avg. 95%
		114	0.5 T d	168	502	>10000	142	254	1656	206	224	288
		172	0.75 T d	202	688	>10000	185	263	1475	206	224	288
		206	0.9 T d	226	630	>10000	221	305	1703	206	224	288
	4	220	T d	5%	50%	95%	5%	50% 95%	5%	Avg. 95%
		110	0.5 T d	238	524	>10000	144	248	1034	206	224	288
		166	0.75 T d	378	378	8166	173	215	831	206	226	288
		198	0.9 T d	338	338	2198	204	248	830	206	228	288