Ludovic Jacquin
email: ludovic.jacquin@inria.fr

Vincent Roca
email: vincent.roca@inria.fr

Mohamed Ali Kaafar
email: mohamed-ali.kaafar@inria.fr

Fabrice Schuler
email: fabrice.schuler@inria.fr

Jean-Louis Roch
email: jean-louis.roch@imag.fr

IBTrack: An ICMP Black holes Tracker

ICMP is a fundamental part of the Internet as it handles the control and error messages. ICMP's treatment by the network and in particular by different routers it may cross is therefore a key aspect driving troubleshooting and diagnosis processes. In this paper we present IBTrack, a tool that aims at characterizing how the network actually treats different ICMP messages from an user point of view. Specifically, we detail a classification algorithm to categorize router behaviors and we introduce its associated refining method which exploits multiple probing protocols. We illustrate the average Internet router behavior and path composition through results gathered from Planet-Lab nodes using a large CAIDA's snapshot of routed /24. We further show that our refining method improves the routers behavior characterization up to 10% for more than 1% of the total number of observed routers.

I. INTRODUCTION

The Internet Control Message Protocol (ICMP) is one of the main protocols used on Internet and more generally on IPv4 and IPv6 networks. It is responsible for transmitting control and error messages over the network (such as routing control, packet treatment error, etc). Although it operates at the second level of the OSI stacks, it is encapsulated on an IP packet.

An example of the ICMP importance is its key implication in the Path Maximun Transmission Unit Discovery (PMTUD) mechanism [START_REF]Rfc 1191: Path mtu discovery[END_REF], [START_REF] Luckie | Inferring and debugging path mtu discovery failures[END_REF], [START_REF] Luckie | Measuring path mtu discovery behaviour[END_REF], which has been developed to avoid IP fragmentation. For instance in IPv4, it sets the don't fragment bit (which is not needed in IPv6 since the protocol does not support router fragmentation). If a router cannot transmit the packet because of its size, it must send back to the source an ICMP "Too Big" (type 3 code 4 on IPv4 and type 2 code 0 on IPv6) packet. Iteratively, the source will lower the size of the packet to match the lowest MTU on a path. The importance of a well-chosen PMTU has been discussed in several previous works, and one key aspect to consider from routers perspective is the number of packets per second to handle [START_REF] Egi | Understanding the packet processing capabilities of multi-core servers[END_REF], [START_REF] Dobrescu | Routebricks: exploiting parallelism to scale software routers[END_REF]. In essence, using the highest possible PMTU value results in a significant bandwidth improvement [START_REF]Myricom 10-gigabit ethernet performance measurements[END_REF] since the packets treatment overhead remains the same regardless of the packet's size. ICMP routing issues can therefore lead to serious connection problems because of their tight links with the PMTU discovery protocol for example.

This motivates our study that aims at analyzing issues related to ICMP packet processing by routers, and in particular so-called "ICMP black holes". The tool presented in this paper, IBTrack, provides users with a thorough analysis of the routers behavior that lie along the path from a source host controlled by the user and any given destination (that is not assumed to be under control). In particular, IBTrack characterizes the routers forging of ICMP error packets and their transport the way back to the source host. It is also important to note that IBTrack should be a "lightweight" tool which only relies on measures performed by the user at the sending host, without requiring any additional resource (i.e. there is no external monitoring and/or vantage point nor any collaboration from the destinations). We further require that IBTrack does not involve any long-term measurement, so it could be used in a timely basis, when a user needs to establish a diagnosis and locate any possible problem origin (similarly to the ping tool).

Let us first compare IBTrack to some related works. In [START_REF] Katz-Bassett | Reverse traceroute[END_REF], authors describe Reverse traceroute a tool that can be used to perform measurement (using ICMP) at the level of Internet routers back to a specified host. The tool is capable of examining the router-level topology and does not look into the behaviors of routers to forge different types of ICMP packets when needed. Authors of Reverse Traceroute also assume users do have access to multiple vantage points distributed across the Internet. Our approach is similar to the concept explored in [START_REF] Katz-Bassett | Studying black holes in the internet with hubble[END_REF], where authors propose a system called Hubble which detects IP-level black holes at a global level of IP connectivity. Besides the fact that Hubble needs periodic probing campaigns of the Internet, which IBTrack tries to avoid, our proposed tool aims to further characterize the reasons to which any ICMP -connectivity related issues are due and provides a fine-grained analysis.

The rest of this paper is organized as follows. Section II presents the model used to describe routers and the assumptions made in our work. section III describes the IBTrack algorithms. Section IV focuses on the methodology and tools used for our measurements. Next, section V presents the results from both coarse-grained Internet and fine-grained ISP levels. Finally, section VI concludes the paper.

II. ASSUMPTIONS AND TAXONOMY

In this section, we define the terminology used to define the source and destination path components. We also take a closer look at the restrictions implied by IBTrack's main goals and our practical assumptions.

A. Path definitions

Let us consider a path from a source S to a destination D, and a router R along this path (there are typically several routers and R is one of them). We assume that the packet flow generated by S and destined to D uses the probing protocol PP. More specifically in this work we consider either ICMP/IP, UDP/IP, or TCP/IP, where IP denotes either IPv4 or IPv6. Figure 1 illustrates the following definitions:

• the initial forward path is the path between S and R;

• the final forward path is the path between R and D;

• and the ICMP return path is the path taken by ICMP packets destined to S and either generated by R or by a router on the final forward path and that go through R;

It is important to note that the initial forward paths and the return ICMP paths are not necessarily identical: routers often use different forwarding strategies depending of routing policies. The forward initial or final paths also potentially depend on the nature of PP, since routers often use protocoldependant forwarding rules. Finally, there are potentially as many ICMP return paths as the number of routers on the forwarding paths.

B. Assumptions

The algorithms we introduce only take into account the information gathered by the source. The main goal of our approach is to characterize routers behavior along a path in order to help users when troubleshooting connectivity issues that are potentially caused by ICMP packets. Therefore we deliberately chose not to rely on external data from vantage point or from the destination during the analysis.

We assume that the routing only depends on the packet's protocol and destination. This is a strong assumption that is however reasonnable in current Internet.

Furthermore, the property used to characterize a router is considered as global, independently from the router's IP interface that has been used during the measurement inferring this property. In other words, we assume that all IP interfaces of that router behave identically.

Finally, since we assume we only control the source S, it is worth noticing that the only ICMP return packet type we can trigger from any router on the network is the Time Exceeded type, which is in theory forged when a packet reaches a router with the IP field "Time To Live" equal to zero). During our measurements, will examine ICMP return paths only by using ICMP TTL packets.

C. Router properties taxonomy

Let us now consider a certain probing protocol PP used by the source. For this protocol PP, each router along the given path is characterized by the following three key properties:

1) Property P1: "R forwards all packets of type PP towards D". Said differently, at R, each incoming packet on the initial forward path is forwarded on the final forward path. This property, of course, does not imply that these packets reach destination D.

2) Property P2: "R is cooperative for packets of type PP". In case R should send an ICMP packet back to the source, either because of an error (e.g. a packet that exceeds the forwarding link MTU) or because the packet asks for a reply (e.g. in case of an ICMP Echo Request), then the ICMP packet is correctly initialized and sent by R on the ICMP return path.

3) Property P3: "ICMP packets are not filtered by R" and by any router on the ICMP return path from R. When R emits an ICMP packet on the return path, this ICMP packet is routed all the way back to S. This property implies that this ICMP packet arrive to S, i.e. they are not filtered by any router on the ICMP return path from R, which is a strong property. In the particular case where the return path contains one or several routers that are also part of the forward path, this property on R implies the same property on these routers.

The router behavior can now be expressed as three logical expressions, one for each property. For instance, for a packet type P: P1.(P2.!P3 + !P2) means that router R forwards these packets and is either cooperative (i.e. sends ICMP packets to S) but these ICMP packets are filtered on the ICMP return path, or is non cooperative (i.e. does not generates ICMP packets back to S).

III. THE ICMP BLACK HOLES TRACKING (IBTRACK) TOOL

In this section we describe our IBTrack tracking algorithm, whose goal is to refine as much as possible the routers behavior, considering only the initial traffic and the backward ICMP path. In a first step, we start by considering a Probing Protocol PP and try to determine the router's behavior. Then, in a second step, we explain how IBTrack can sometimes refine the analysis by crossing the results achieved independently with several probing protocols.

A. IBTrack base algorithm

The base algorithm is described in figure 2. At the beginning there is no knowledge at all for router R, which corresponds to the state at the top of the graph. In that case each property is "undecided" (i.e. it can be true or false), which is expressed by: P + !P.

Then the algorithm is composed of two main steps:

• test the forwarding path through R;

• test the ICMP return path from R; We now detail these two steps:

P1 + !P1 P2 + !P2 P3 + !P3 S sends PP packet to D Did D receive PP packets from S? P1 P2 + !P2 P3 + !P3 S sends PP packet with TTL=0 on R Did S receive ICMP packets from R? P1 + !P1 P2 + !P2 P3 + !P3
S sends PP packet with TTL>1 on R Did a router "after" R receive PP packets from S?

P1 + !P1 P2 + !P2 P3 + !P3 S sends PP packet with TTL=0 on R Did S receive ICMP packets from R? P1 P2 P3 Case A P1 P2.!P3 + !P2 P3.!P2 + !P3 Case B P1 + !P1 P2 P3 Case C P1 + !P1 P2 + !P2 P3 + !P3 Case D Yes No Yes No Yes No Yes No
Fig. 2: Base algorithm for router's behavior characterization 1) Forwarding path through R: We first want to test the P1 property (forwarding property towards D). It can be verified in two ways: either the destination replied, which implies that it has received an incoming packet, or we can assert that a router on the final forward path received an incoming packet. This second method is typical of situations where source S received ICMP packets from a router on the final forward path (i.e. after router R).

2) ICMP return path from R: We now want to check the P2 and P3 properties. For the P2 property we must verify that the router forges the ICMP error packet and emits it on the ICMP return path. This could be tested by monitoring the ICMP return path. Nevertheless, as we do not control any router of the ICMP return path, we only rely on the reception by the source of the ICMP error packet. For the P3 property, since it concerns the ICMP return path itself, it can only be evaluated with the reception of ICMP error packets by the source.

3) Classification details: We now present the result of this classification, which consists in four cases, depending on the properties that could be verified or not. a) Case A: P1.P2.P3 This is the "ideal" router which forwards the initial packets and forges ICMP error packets that all arrive to the source. b) Case B: P1.(P2.!P3 + !P2).(P3.!P2 + P3) P1 is verified, but the interesting part of this case comes from the P2 property. In this case, the router forwarded the initial packets and also dealt with them (i.e. filtering and packet processing). However the source did not receive any ICMP error packet whereas such a packet should have been forged since the source S sent a packet whose TTL ("Hop Count" in case of IPv6) has been decremented to 0 by R). It leads to two possibilities: This base algorithm is applied to every routers in a path (defined by its source and destination) for protocol PP, by increasing the TTL value It results in a collection of cases (A, B, C or D) for each router along the forward path.

• P2 ⇒ !P3: the ICMP error packet is forged but is blocked along the ICMP return path. P1 P2.!P3 + !P2 P3.!P2 + !P3 PP 1 -Case B P1 P2 P3 PP 2 -Case A P1 P2.!P3 + !P2 P3.!P2 + !P3 PP 2 -Case B P1 + !P1 P2 P3 PP 2 -Case C P1 + !P1 P2 + !P2 P3 + !P3 PP 2 -Case D P1 !P2 P3 PP 1 -Case B' P1 P2.!P3 + !P2 P3.!P2 + !P3 PP 1 -Case B P1 !P2 P3 PP 1 -Case B' P1 P2.!P3 + !P2 P3.!P2 + !P3 PP 1 -Case B

B. IBTrack refining algorithm

Let us now explain how we can refine the previous algorithm that depends on the probing protocol PP that has been used for probing the network. Let us first notice that the P3 property concerns exclusively the ICMP return path which is the same independently from the probing protocol PP 1 . Thus, for a given router R, we can use the deductions made on the P3 property using a probing protocol PP1 to refine the deductions made using another probing protocol PP2. This is the key idea of our refining algorithm. However, in order to be able to apply this refinement, we must first identify the same router for both probing protocols, whereas we do not necessarily have its identity (it can be a wildcats "*" in a traceroute trace). This aspect is addressed in section IV-C, and for the moment we assume that this assertion as true.

The combination of the results of the algorithm for two distinct probing protocols is detailed in figures 3 and 4, for router R. Since we can only use deductions for the P3 property, cases A or C cannot be refined (property P3 is already known in that case). Therefore we focus on cases B and D. 1) Initial case B for probing protocol PP1 (Fig. 3): This is the most interesting case. When property P3 is verified with the second probing protocol, PP2, from the second conditional we can deduce that P3 ⇒ !P2. Thus, if probing with protocol PP2 ends in cases A or C, we can conclude that the P3 property is verified for protocol PP1, which implies that P2 is not. This creates a subcase for B that we call B':

P1 + !P1 P2 + !P2 P3 + !P3 PP 1 -Case D P1 P2 P3 PP 2 -Case A P1 P2.!P3 + !P2 P3.!P2 + !P3 PP 2 -Case B P1 + !P1 P2 P3 PP 2 -Case C P1 + !P1 P2 + !P2 P3 + !P3 PP 2 -Case D P1 + !P1 P2 + !P2 P3 PP 1 -Case D' P1 + !P1 P2 + !P2 P3 + !P3 PP 1 -Case D P1 + !P1 P2 + !P2 P3 PP 1 -Case D' P1 + !P1 P2 + !P2 P3 + !P3 PP 1 -
P1+!P2+P3.
2) Initial case D for probing protocol PP1 (Fig. 4): This case can only be refined on the P3 property itself since there is no material implication or link with another property. This only happens when we end in cases A or C with the second probing protocol.

3) Refining using multiple probing protocols: It is of course possible to use more than one probing protocol to refine the results obtained with PP1. Nevertheless, for this to be true, we must assert that the paths are the same for every combination of the probing protocol and of the refining protocol in order to consider that P3 property is common. In the traceroute outputs of figure 5a, the 16 th router is a B case for both protocols since: (i) it forwards the initial protocol (the 17 th router replied which implies that it has received at least one initial packet); but (ii) no ICMP return packet is received by the source. There is no refining possible in this case and: P1.(P2.!P3 + !P2).(P3.!P2 + P3)

C. Examples

In the traceroutes outputs of figure 5b), the 7 th router is a B case for UDP probing whereas it is an A case for ICMP probing (the ICMP return packet is received by the source). As the ICMP return path is the same (source, protocol and destination are identical), we deduce that for UDP packet the router does not forge the ICMP return packets. The categorization for UDP probing is therefore refined to the B' case: P1+!P2+P3

The formal model and algorithms being described, the next section will focus on the methodology and on the tools used for Internet measurements.

IV. MEASUREMENT METHODOLOGY

Before describing the tools and methodology used during the measurements, we detail how we chose the IP addresses to probe.

A. Selection of destination's IP

To probe a meaningful set of Internet, we used a CAIDA [START_REF] Hyun | The caida ipv4 routed /24 topology dataset[END_REF] snapshot of the routed /24 IPv4. This snapshot is composed of a part of the routed /24 (approximately 10,000) and gives for each a random IP in the /24. It is important to notice that there is not necessarily a running host behind each IP in the snapshot.

B. Tools

To perform probing, we used scamper [START_REF] Luckie | Scamper: a scalable and extensible packet prober for active measurement of the internet[END_REF] that implements, among others, the paris-traceroute [START_REF] Augustin | Avoiding traceroute anomalies with paris traceroute[END_REF] functionnality over the Planet-Lab [START_REF]Planet-lab[END_REF] nodes.

1) Paris-traceroute: This is an improved version of the classical traceroute tool. In particular this tool reveals more routers and links over the path explored and removes some false links inferred by the usual traceroute. This is therefore highly benefic to our needs.

2) Planet-Lab: This is a world wide network aiming to help academic and industrial researchers focusing on new network services. In our test, it provides multiple servers across the world that we use to probe the same destination set of IP addresses. Therefore it helps us improving the coverage of the possible paths to the selected set of IP addresses.

3) Scamper: This is a program that implements most of the classical Internet measurement tools (like ping or paris-traceroute) that can be launched to run in parallel. In particular it is developed to run on Planet-Lab nodes.

One of the features of scamper we use is that after five unidentified routers (i.e. not responding, thus appearing as "*" in the traceroute output) it stops. Figure 6c shows a traceroute terminating like this. This situation can result from different causes, like for instance a single router filtering packets, or multiple routers not replying. Nevertheless, we will deliberately take into account only the first "*" during our tests.

C. Identifying common path

One of the main challenges before applying the refining algorithm which uses multiple probing protocols is to be able to match paths that seem to be different. In principle the presence of an unidentified router on a path (i.e. a "*" in the traceroute output) prevents path matching for the final forwarding path. But this situation must be improved. Indeed, if the router has to be identified to allow the application of our refining algorithm, then we would only perform refining between A and C cases, which is of little interesting because these cases already share the same state for the property P3.

1) Handling unknow routers:

The first and straightforward assumption we make consists in ignoring unknown routers from a path when matching two paths with the same source and destination but with different probing protocols.

Formally, there is no proof for this assumption, but we are convinced that given the maximum number of consecutive unknown routers we allowed (5 before traceroute halts) for a protocol, if two paths share the same number of routers and if the identified routers are the same, then the two paths are most likely the same.

For instance, figure 6a exhibits two situations considered as belonging to the same path by the path matching algorithm.

2) Handling tunnels: We also wanted to add to the path matching algorithm the capability to identify tunnels. More particularly, we want to isolate path divergences due to perflow routing (for example through MPLS tunnel [START_REF] Donnet | Revealing mpls tunnels obscured from traceroute[END_REF]). Therefore, we improved our algorithm by allowing some consecutive identified routers to be different on the same path. This feature adds itself to the previous one; thus "consecutive identified routers" do not involve potential unknown routers. Nevertheless, the number of routers and their position on the Figure 6b illustrates two paths being considered as identical by the path matching algorithm when allowing two different consecutive identified routers.

The following section presents large scale measurements gathered from Planet-Lab nodes and a fine-grained analysis of the algorithm results.

V. MEASUREMENT RESULTS

We first detail the results obtained from 30 Planet-Lab nodes to a set of CAIDA's routed /24 composed of more than 10,000 IP addresses. Then we detail the results obtained from a typical ISP user probing the same set of IP addresses, using ICMP, UDP, but also TCP.

A. Coarse-grain Internet results

1) Paths characterization: figure 7 represents the normalized number of router per case. On average, a traceroute is composed of 16.78 routers distributed as followed:

• A case -14.5 routers • B case -0.45 router • C case -1 router • D case -0.83 router Figure 8 shows the distribution of cases depending on the distance between the router and the initial source. Note that when handling a "premature end" of a traceroute (i.e. when finishing with five "*"), the D cases are only counted once (it corresponds to the closest router from the source) as explained in section IV-B3.

This figure shows that the closer a router is to the source, S, the mostly it will be of case A. On the opposite, the further a router is, the highest the probability to fall in C and D cases, which in turn comes at the expense of a lower probability of the router would be in case A. Perhaps not surprisingly, this shows that the longer the path is the higher the chances the routers would behave according to the D unexplicative case. 2) Routers characterization: figure 9 summarizes the classification of routers encountered during the probing of every IP of our destination set, from each Planet-Lab nodes, using the ICMP and UDP protocols. We remind that for the D case, results are divided by 5 due to the "premature end" of a traceroute (section IV-B3). In opposition to the two previous figures, we add that the router seen are counted only once as we want to survey the Internet's average router behavior (previously we were focusing on the path instead).

Figure 10 details the results of the refining algorithm applied to the Planet-Lab outputs. We first configured the path comparison algorithm to perform an exact path matching and then changed it to allow a certain number of different routers along the path. With a strict path matching configuration, we refine 0.16% (ICMP) and 0.31% (UDP) of the B case. This ratio grows up to 0.87% (ICMP) and 1.12% (UDP) when allowing up to five different consecutive identified routers. This measure (particularly the small difference when allowing five or ten different routers) can be associated to the fact that more than 90% of MPLS tunnels (which are present in at least 30% of Internet's traceroute) are at most five hops long [START_REF] Donnet | Revealing mpls tunnels obscured from traceroute[END_REF].

B. Fine-grain results

Let us now consider the results obtained by a single user using a single host, connected through a given ISP. The set of IP addresses is the same, and we use the ICMP, UDP, but also TCP probing protocols.

Figure 11 presents the results of the classification algorithm applied on the three probing protocolsICMP, UDP and TCP. Although the general distribution is globally the same, there are more A cases and fewer B cases than with Planet-Lab results (figure 9).

Figure 12 details the refining algorithm used for TCP, which refines TCP results either using ICMP solely, or UDP, or both probing protocols. It illustrates well the fact that, when using multiple probing protocols to refine a single other protocol-based properties output, the results are lower than what could be expected from the sum of the refining using single protocols outputs, as explained in section III-B3. Indeed, we observe that the improvement of the refining algorithm while using two protocols to refine a single one is small. For instance, when allowing five different consecutive identified routers, the improvement is: 5.2% when comparing ICMP to ICMP+UDP, and 14.5% between UDP and ICMP+UDP. We interpret this at a router level by bearing that there is an explicit policy on these routers that forbids it to forge ICMP Time Exceeded packets when the initial protocol is TCP. Nevertheless, when using both ICMP and UDP to refine the B case of TCP, we see that we are able to completely specify the behavior of 10% of the routers encountered. Globally, we manage to refine the classification of 1% of the routers seen.

VI. CONCLUSION

This work has investigated the problem of ICMP black holes tracking. Our contributions are fourfold: we first introduced a taxonomy for router properties, regarding ICMP. We have shown that three fundamental properties are sufficient to classify Internet routers with respect to the ICMP protocol processing.

Secondly we have introduced a methodology to classify the routers of a path, using a given probing protocol, into four cases. These cases represent the possible logical relationships between the three fundamental properties. Then we introduced a refinement procedure, using two or more probing protocols, that enables to better classify the routers. This approach does not require any control from within the Internet or at destination. It can therefore be immediately ported to any operating system and used by any user.

Third, we designed the IBTrack tool that implements the above algorithm.

Fourth, we performed several large scale experiments to analyze our tool and methodology in real conditions. In these experiments, we also report several insights on the Internet routers ICMP processing.

We therefore believe that this work is a significant step forward. In particular it will help many network administrators to better understand traceroute outputs and debug complex ICMP related problems within routers.

This work deliberately makes very minimalist assumptions on what users might control while processing and diagnosing network failures messages (only the source is controlled). As such, the proposed approach can be easily extended by relaxing some constraints (e.g. the user also controls the destination) or by revisiting some of our assumptions on network functionalities (e.g. the presence of vantage points).

Fig. 1 :

 1 Fig. 1: Paths definitions.

Fig. 3 :

 3 Fig. 3: Refining case B with another probing protocol

Fig. 5 :

 5 Fig. 5: Two examples of multiple probing protocol results

Figure 5

 5 Figure 5 illustrates the refining algorithm through two examples extracted from our experiments.In the traceroute outputs of figure5a, the 16 th router is a B case for both protocols since: (i) it forwards the initial protocol (the 17 th router replied which implies that it has received at least one initial packet); but (ii) no ICMP return

Fig. 7 :

 7 Fig. 7: Normalized (per traceroute) number of routers per case (with standard deviation)

Fig. 8 :Fig. 9 :

 89 Fig. 8: Case percentage with regards to the distance to the source S

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Mean number of router refined from B to B' case depending of the path matching algorithm configuration (with standard deviation)

Fig. 12 :

 12 Fig. 12: Refining of TCP probing protocol B case

  Case D

	Fig. 4: Refining case D with another probing protocol
	ICMP	UDP
	15 213.248.90.250 15 213.248.90.250
	16 *	16 *
	17 205.171.11.98	17 205.171.11.98
	(a) two B cases
	ICMP	UDP
	6 195.83.166.161 6 195.83.166.161
	7 193.51.181.222 7 *
	8 193.51.189.202 8 193.51.189.202

(b) A case and B case

This is a direct consequence of our assumption that routing only depends on the protocol along the return path, here ICMP, and the destination (here S).


Hai Canh Vu

Phuc Do Van b Anne Barros

Christophe Bérenguer

Maintenance activities planning and grouping for complex structure systems

Keywords: Grouping maintenance, Economic dependence, Multi-component system, Optimization

This paper presents a dynamic grouping maintenance strategy for complex systems whose structure may lead to both positive and negative economic dependence which imply that combining maintenance activities is cheaper (or more expensive respectively) than performing maintenance on components separately. Binary Particle Swarm Optimization (BPSO) algorithm is used to find optimal grouping planning which is NP-hard combinatorial problem. The proposed grouping maintenance strategy based on the rolling horizon approach can help to update the maintenance planning by taking into account short-term information which could be changed with time. A numerical example of a 10 components system is finally introduced to illustrate the use and the advantages of the proposed approach in the maintenance optimization framework.

INTRODUCTION

For multi-component systems, economic dependence which can be either positive or negative have been considered in maintenance optimization. Positive economic dependence implies that costs can be saved when several components are jointly maintained instead of separately. Negative economic dependence between components occurs when maintaining components simultaneously is more expensive than maintaining components individually. A number of maintenance optimization methods and maintenance strategies have been proposed and developed in the literature, see for example [START_REF] Nicolai | Optimal maintenance of multi-component systems: a review[END_REF][START_REF] Cho | A survey of maintenance models for multi-unit systems[END_REF][START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF][START_REF] Bouvard | Condition-based dynamic maintenance operations planning & grouping. application to commercial heavy vehicles[END_REF]. In such papers, they focus mainly on series structures which lead directly to positive economic dependence [START_REF] Dekker | Integrating optimisation, priority setting, planning and combining of maintenance activities[END_REF][START_REF] Dekker | Joint replacement in an operational planning phase[END_REF][START_REF] Laggoune | Preventive maintenance scheduling for a multicomponent system with non-negligible replacement time[END_REF][START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF]. Recently, in [START_REF] Pham | Optimal (τ , t) opportunistic maintenance of a k-out-of-n: G system with imperfect pm and partial failure[END_REF][START_REF] Sheu | Optimization problems in k-out-of-n systems with minimal repair[END_REF], both positive and negative economic dependence have been studied by introducing an opportunistic maintenance which attempts to balance the negative aspects and the positive aspects of grouping maintenance. The proposed maintenance models are however only applicable on limited class of systems with particular structure. From a practical point of view, with growing up quickly of industry, systems become more and more complex. Hence, the proposed models may no longer be used. The grouping maintenance problem remains widely open.

Moreover, in the grouping maintenance framework, static models with an infinite planning horizon are usually used in case of stable situations [START_REF] Dekker | Integrating optimisation, priority setting, planning and combining of maintenance activities[END_REF][START_REF] Dekker | Joint replacement in an operational planning phase[END_REF][START_REF] Laggoune | Preventive maintenance scheduling for a multicomponent system with non-negligible replacement time[END_REF]. Recently, dynamic models have been introduced in order to change the planning rules according to short-term information (e.g. failures and varying deterioration of components) by using a rolling horizon approach [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF]. The primary objective of this paper is thus to develop the rolling horizon approach for grouping maintenance of complex structure systems with both positive and negative dependence. For such systems, to find optimal grouping planning, dynamic programming (see [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF][START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF]) is no longer usable since the combinatorial problem can be formulated as a set partitioning problem, which however can be NP-hard due to the complex structure impact. To this end, the BPSO [START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Franken | Investigating binary pso parameter influence on the knights cover problem[END_REF][START_REF] Kennedy | A discrete binary version of the particle swarm algorithm[END_REF], recognized as a general search strategy which is simple for implementation and very efficient for global search, is used. This paper is organized as follows. Section 2 is devoted to the description of general assumptions and grouping maintenance problem. The impact of system structure on grouping maintenance is also analyzed. Section 3 focuses on the development of the rolling horizon approach in the context of complex structure systems. A simple numerical example is introduced in Section 4. Finally, the last section presents the conclusions drawn from this work.

PROBLEM FORMULATION

General assumptions

Consider a complex structure coherent system consisting of N independent components in which component state is either operational or failure one. For the rate of occurrence of failures, denoted r i (t), we take a polynomial function with scale parameter λ i > 0, and shape parameter β i > 1:

r i (t) = β i λ i ( t λ i ) β i -1 (1)
Both corrective and preventive maintenance are considered for each component and we assume also that preventive and corrective maintenance durations can be neglected. According to the system structure, we consider two kinds of components: non-critical component (the system can be still functioning when the component stops); critical component (a shutdown of the component for whatever raison leads to a shutdown of the system). By definition, maintenance of a critical component may lead to an additional cost that could be relying on, for example, production quality loss or/and restart system costs... see [START_REF] Wmj Geraerds | The cost of downtime for maintenance: preliminary considerations[END_REF]. Two kinds of shutdown costs are here considered:(i) planned shutdown cost, denoted C p sys , that is incurred when executing preventive maintenance on a critical component;(ii) unplanned shutdown cost, denoted C u sys , that is incurred when executing corrective maintenance on a critical component. Since the failure is random and preventive maintenance is planned, an unplanned shutdown is therefore more costly than a planned shutdown, C c sys ≥ C p sys .

For component i (i = 1, ..., N ), we assume that after a preventive maintenance (PM) action the component is considered "as good as new". The preventive maintenance cost can be divided into three parts: a fixed setup cost, denoted S, represents the common cost for all preventive maintenance activities. The setup cost depends only on the system. For example, the setup cost can be composed the cost of crew traveling and the preparation cost (erecting a scaffolding or opening a machine); a specific cost, denoted c p i , depending on the component; an additional cost C p sys is incurred if component i is critical. That mean that this additional cost depends on the component position in the system structure and the system characteristics. Let π i is an indicator function (π i = 1 if component i is a critical one, and otherwise π i = 0 if component i is a not-critical one). As consequence, if a PM of component i is carried out, we have to pay the following cost:

C p i = S + c p i + π i • C p sys (2)
To define the preventive maintenance cycle for each component, a long-term average cost based approach [START_REF] Dekker | Maintenance of light-standards-a case-study[END_REF] is used and presented in Section 3. Between two preventive maintenance intervals, if component fails, a minimal repair is then immediately performed to restore the component to an operational state "as bad as old". Let c r i denote the specific repair cost of component i. As in the PM case, when a corrective action is carried out for component i, it requires a corrective maintenance cost, denoted C c i .

C c i = S + c r i + π i • C u sys (3)

Grouping maintenance and the impact of the system structure

In the spirit of grouping maintenance strategies see [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF][START_REF] Nicolai | Optimal maintenance of multi-component systems: a review[END_REF], grouping maintenance can reduce the maintenance cost. The setup cost can be shared when several maintenance activities are performed simultaneously since simultaneous maintenance executions on several (group) components require usually only one set-up cost [START_REF] Dekker | Maintenance of light-standards-a case-study[END_REF][START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF]. Note well, however that when a group of components are preventively performed together, the maintenance cost could be indirectly penalized with the reduction of components useful life if the maintenance dates are advanced and with the increasing of components failure probability which could imply a system immobilization if the maintenance dates are too late.

Moreover, the maintenance cost of a group may depend on the position of the group's components.

We consider two kinds of groups: non-critical group (the system can be still functioning when all components of the group stop); critical group (the maintenance of the group leads to system stop).

As consequence, maintenance on a critical group leads to a planned shutdown cost. Three following cases are considered: (i) if the critical group contains at least two critical components, the planned shutdown cost can then be shared; (ii) if the critical group contains one critical component, the planned shutdown cost remains then unchanged; (iii) if the critical group does not contain any critical component, the planned shutdown cost can then be penalized.

In order to find the optimal groups which could balance to minimize the system maintenance costs on the scheduling horizon, the rolling horizon approach introduced recently in [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF][START_REF] Dekker | Maintenance of light-standards-a case-study[END_REF], will be developed by taking into account the system structure, see Section 3. An optimization algorithm, BPSO [START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Franken | Investigating binary pso parameter influence on the knights cover problem[END_REF], will be also used.

EXTENSION OF ROLLING HORIZON APPROACH

The developed rolling horizon approach consists of 5 phases: decomposition, tentative planning, economic profit formulations, grouping optimization using BPSO and rolling horizon.

Phase 1: Decomposition

For each individual component, we formulate an infinite horizon maintenance model to obtain the optimal individual preventive maintenance cycle that minimizes the component maintenance costs.

Let M i (x) denote the expected deterioration cost incurred in x time units since the latest operation on component i. According to the minimal repair policy, M i (x) can be expressed as the following:

M i (x) = C r i • ∫ x 0 r i (y)dy. (4)
From equations ( 1) and ( 4), we obtain:

M i (x) = C r i • ( x λ i ) β i If component i is preventively maintained at x, the expected cost within interval [0, x] is: E i (x) = C p i + M i (x) = C p i + C r i • ( x λ i ) β i (5)
By applying the renewal theory when executing the preventive maintenance of component i every x time units amount, the long-term average cost of component i can be determined as the following:

ϕ i (x) = E i (x) x = C p i + C r i • ( x λ i ) β i x (6)
Since ϕ i (x) is a convex function, the optimal interval length for the preventive maintenances on component i, denoted x * i , can be obtained by setting the derivative of ϕ i (x) to zero:

x * i = arg min x ϕ i (x) = λ i β i √ C p i C r i (β i -1) . ( 7
)
And the minimal long-term average maintenance cost of component i:

ϕ * i = ϕ i (x * i ) = C p i β i x * i (β i -1) (8)
The minimal long-term average maintenance cost of the system assuming that all components are optimally and individually maintained:

ϕ * sys = N ∑ i=1 ϕ * i . ( 9
)
The optimal interval length x * i which represents a nominal PM frequency of component i will be used to define tentative execution times in the tentative planning phase.

Phase 2: Tentative planning

The idea of this phase is to establish a tentative maintenance planning based on the individual preventive maintenance cycles.

Let t begin denote the current date and i j be jth maintenance of component i since t begin and t e i is the operational time elapsed from the last preventive maintenance of component i before t begin . Without loss generality we can set t begin = 0. Based on the individual maintenance rules of phase 1, the current state of component i, the tentative execution time of operation (or activity) i j , denoted t i j , is determined as follows, see [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF]:

t i 1 = t begin -t e i + x * i , t i j = t * i j-1 + x * i if j > 1 (10)
In order to evaluate the performance of grouping maintenance strategy, a finite planning horizon is usually defined according to the current date t begin and the ending date t end which guarantees that all components are preventively maintained at least one time in the horizon interval. It is shown in [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF], t end can be determined as follows:

t end = max i t i 1 .
After this phase, the tentative execution times of all individual maintenance activities in the scheduling horizon are identified. In the next phases we try to perform simultaneously several maintenance activities by changing their maintenance execution times to minimize total maintenance cost.

Phase 3: Economic profit formulations

The idea of this phase is to formulate the economic profits when preventive maintenance activities are simultaneously carried out. To this end, we will firstly find all cost-effective groups which can help to identify secondly an optimal grouping maintenance planning.

Cost-effective group Assume now a group of several different maintenance operations i j (i, j = 1, 2, ...), denoted G k , are simultaneously performed at time t. The corresponding economic profit of this group can be divided into three parts as follows:

Since the execution of a group of m maintenance operations requires only one set-up cost, the group G k yields a cost reduction:

U G k = (m -1) • S ( 11
)
The penalty costs due to the changing of components execution maintenance date:

△H 1 G k (t) = ∑ i j ∈G k h i (t -t i j ) = ∑ i j ∈G k h i (△t i j ), (12)
h i (△t i j ) is the penalty cost due to the movement execution date of the maintenance operation i j . Since operation i j is actually executed at time

t = t i j + △t i j (△t i j > -x * i ) instead of t i j . h i (△t i j ) = E i (x * i +△t i j )- ( E i (x * i )+△t i j .ϕ * i ) = C r i • ( x * i + △t i j λ i ) β i -C r i • ( x * i λ i ) β i -△t i j • C p i β i x * i (β i -1) (13)
The optimal execution time of the group G k , denoted t G k , can be found when the △H 1 G k (.) searches its minimal value △H 1 * G k . That is:

△H 1 * G k = △H 1 G k (t G k ) = min t ∑ i j ∈G k h i (t G k -t i j ) (14)
If G k is a critical one, the grouping can conduct to an additional cost that is expressed as follows:

△H 2 G k = C G k -C notG k , (15)
The total planned shutdown costs when all components of the group are separately maintained:

C notG k = C p sys • ∑ i j ∈G k π i ; (16) Let π G k is an indicator function (π G k = 1 if G k is a critique group, and π G k = 0 if G k is a not-critique group).
The planned shutdown cost when all components of the group are performed together:

C G k = C p sys • π G k (17)
Thanks to equations ( 15), ( 16) and (17), △H 2 G k can be written as:

△H 2 G k = C p sys • (π G k - ∑ i j ∈G k π i ). ( 18
)
Note well that if the group G k is not critical, all components of the group are then not critical, △H 2 G k is hence equal to zero. If G k is a critical group, △H 2 G k can be positive or negative depending on the number of critical components of the group (see again Section 2).

From Equations ( 11), [START_REF] Wmj Geraerds | The cost of downtime for maintenance: preliminary considerations[END_REF], and (18), the economic profit of group G k , denoted EP (G k ), can be written as the following:

EP (G k ) = U G k -△H 1 * G k -△H 2 G k (19) (i) If EP (G k
) < 0 the grouping maintenance is more expensive than performing maintenance on components separately, i.e. the grouping implies a negative economic dependence. As consequence, the components of the group should not be simultaneously maintained;

(ii) If EP (G k ) ≥ 0, the maintenance cost can be saved when components are jointly maintained, i.e. grouping leads to a positive economic dependence. And the group is called cost-effective.

Grouping structure Based on all cost-effective groups, a grouping structure (or partition of all maintenance operations in the scheduling interval), denoted SGM , can be identified. As consequence, the total economic profit is calculated as follows:

EP S(SGM ) = ∑ G k ∈SGM EP (G k ) = ∑ G k ∈SGM (U G k -△H 1 * G k -△H 2 G k ). ( 20
)
An optimal grouping structure when the total economic profit searches its maximal value:

SGM optimal = arg max SGM EP S(SGM ) (21)
Note well that the maximum number of groups to be considered is here 2 n -1 for n maintenance activities. The finding of optimal grouping structure becomes difficult since the combinatorial problem can be formulated as a set partitioning problem, which however can be NP-hard. The dynamic programming algorithm proposed in [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF][START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF] can not therefore be used. To solve this problem, BPSO algorithm will be used and described in next paragraph.

Phase 4: Grouping optimization using BPSO

Binary particle swarm optimization was proposed by Kennedy and Eberhart in 1997, see [START_REF] Kennedy | A discrete binary version of the particle swarm algorithm[END_REF], to solve optimization problems in binary search spaces (0 or 1). This method is loosely modeled on the flocking behavior of birds. A solution to a specific problem is presented by a position of a particle. A swarm of fixed number of particles is generated and each particle is initialized with a random position in the search space. Each particle flies through the search space with a velocity. At each iteration, a new velocity value of each particle is calculated based on its current velocity, the distance from its personal best (the best position of particle so far), and the distance from the global best (the best position of all particles). Once the velocity of each particle is updated, the particles are then moved to the new positions. The general structure of BPSO is presented in Fig 1 .   Coding A coding phase determines how the problem is structured in the algorithm. In our problem, each grouping structure (a solution) is presented by a binary matrix of (n -1) rows and n columns. Each row represents a group in which entry (k, j) = 1 if maintenance activity j is performed in group k and entry (k, j) = 0 for otherwise. For example considering a feasible solution of 5 preventive

NUMERICAL EXAMPLES

The purpose of this section is to show how the proposed grouping maintenance strategy can be used in preventive maintenance optimization of complex systems. The impact of the system structure on the grouping maintenance planning will be analyzed. A short comparison with previous grouping approach [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF] will be also discussed.

Consider a 10 components system whose structure is shown in Figure 2. When a component fails, it is immediately maintained according to a minimal-repair policy. Corrective maintenance restores the component involved into a state "as bad as old". We assume failure rate of a component i (i = 1, . . . , 10) is described by a Weibull distribution with scale parameter λ i > 0, and shape parameter β i > 1. Table 1 reports the random data for 10 components. For set-up cost and shutdown system costs, we take S = 10 and C p sys = 40, C u sys = 45 respectively.

Maintenance planning with negligible system structure

Consider now maintenance strategies in which the system structure is not considered. We assume that the system is stopped when performing maintenance on any component, i.e. π i = 1, ∀i = (1, . . . , 10).

Individual maintenance planning

We consider an optimal maintenance planning in which all components' preventive maintenance are separately performed. To define the optimal preventive maintenance cycle of components, the long-term maintenance cost model was used. C p i , C c i , x * i , ϕ * i and t 1 i are calculated by substitution of the data input and π i in Equations ( 2),( 3),( 7),( 6),(10) respectively. All results are shown in table 2. The average maintenance cost of the system when all components are individually maintained is finally:

ϕ SG sys = N ∑ i=1 ϕ * i = 9.2662 (23)
Grouping maintenance planning To establish an optimal grouping maintenance, the proposed grouping maintenance strategy was used. Based on the components' individual maintenance execution date, the scheduling horizon is [0, 212.22] (t begin = 0, t end = max t 1 i = 212.22). After 500 iterations, BPSO optimization algorithm provided an optimal grouping maintenance planning. The results are reported in Table 3. The total cost saving is EP S = 393.9172 which is equivalent to the average saving ∆ϕ sys = EP S t end = 393.9172 212.22 = 1.8562. Hence, the average maintenance cost of the system with proposed grouping maintenance planning is:

ϕ SG sys = ϕ SG sys -∆ϕ sys = 7.41 (24)
The result shows that grouping maintenance is cheaper than performing maintenance on components separately. Note well that under the assumptions above, the system can be considered as a series structure for which the rolling horizon approach [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF] can be used. By applying this approach, we obtain the same results shown in Table 3.

Maintenance planning with taking into account the system structure

Individual maintenance planning According to the system structure, indicator functions π i (i = 1, ..., 10) can be easily determined. In the same manner, C i p , C r i , x * i , ϕ * i , t 1 i are then calculated and showed in Table 4.

The average maintenance cost of the system when all components are individually maintained is finally:

ϕ SG sys = N ∑ i=1 ϕ * i = 6.3897 (25)
The result shows that the consideration of the system structure in maintenance can save the maintenance cost. If the system structure is taken into account, the individual maintenance planning is even cheaper than an optimal grouping maintenance planning in which the system structure is ignored. Grouping maintenance planning In order to establish a grouping maintenance planning, we used again the proposed grouping approach. According to the components' individual maintenance execution date, the planning horizon is [0, 300]. Table 5 reports the results obtained by using BPSO optimization algorithm with 500 iterations. From ( 23), (24), ( 25) and (26), we obtain an important ranking for the average maintenance cost of the system: ϕ SG sys > ϕ SG sys > ϕ SG sys > ϕ SG sys . The results assert that grouping can save the maintenance cost and that the system structure can take an important rule and should be taken into account in maintenance planning optimization.

CONCLUSIONS

In this work, the rolling horizon approach introduced recently in [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF] is developed for grouping maintenance planning of complex systems whose structure can take both positive and negative impact on maintenance grouping. The numerical results show that taking into account the system structure into grouping maintenance procedure may reduce significantly maintenance cost. In order to find an optimal grouping planning, Binary Particle Swarm Optimization algorithm is proposed to use. This optimization algorithm can help to solve combinatorial grouping problem which can be NP-hard.

Our future research work will focus on the development of the proposed method for systems which have the stochastically and structurally dependent between components and also with the taking into account the duration of maintenance activities.

Figure 2 :

 2 Figure 2: Structure system

Table 1 :

 1 Data of 10 components

	Component i	1	2	3	4	5	6	7	8	9	10
	λ i	259	270	280	249	297	260	285	285	250	260
	β i c p i c r i	1.90 115 42	2.00 125 32	2.00 125 22	1.95 105 35	1.95 145 40	2.00 125 35	2.00 165 20	2.00 145 20	1.90 125 40	2.00 135 35
	t e i	184.37 214.07 295.11 153.61 364.70 150.33 472.54 459.55 155.71 226.71

Table 2 :

 2 Values

				of C p i , C c i , x * i , ϕ * i and t 1 i with negligible system structure
	Component i	1	2	3	4	5	6	7	8	9	10
	π i C p i C r i	1 165 97	1 175 87	1 175 77	1 155 90	1 195 95	1 175 90	1 215 75	1 195 75	1 175 95	1 185 90
	x * i	362.09 382.93 422.12 337.83 440.90 362.55 482.54 459.55 364.47 372.77
	ϕ * i	0.9620 0.9140 0.8292 0.9418 0.9078 0.9654 0.8911 0.8487 1.0136 0.9926
	t 1 i	177.75 168.87 127.01 184.22 76.20 212.22 10.00	0.00	208.76 146.05

Table 3 :

 3 Optimal grouping planning with negligible system structure Group Group components Group execution date (t G )

	Group cost saving

Table 4 :

 4 Values of individual optimization with taking into account the system structure

	Component i	1	2	3	4	5	6	7	8	9	10
	π i C p i C r i	0 125 52	0 135 42	0 135 32	0 115 45	0 155 50	0 135 45	1 215 75	1 195 75	0 135 50	0 145 45
	x * i	434.37 484.07 575.11 413.61 544.70 450.33 482.54 459.55 445.71 466.71
	ϕ * i	0.6075 0.5578 0.4695 0.5707 0.5841 0.5996 0.8911 0.8487 0.6394 0.6214
	t 1 i	250	270	280	260	180	300	10	0	290	240

Table 5 :

 5 Optimal grouping planning with taking into account the system structure The total maintenance savings cost is EP S = 18.4307 + 39.3498 + 49.9538 = 107.7343 and the maintenance cost saving par time unit is:

	Group Group components Group execution date (t G ) Group cost saving
	1	1,5,10	226.60	18.4307
	2	2,3,4,6,9	280.31	39.3498
	3	7,8	5.00		49.9538
		ϕ SG sys = ϕ SG sys -	107.7343 300	= 6.0306.	(26)

Initialize algorithm constants

This solution can be presented by the matrix (22a).

Create randomly a swarm of particles and its velocities A swarm of s particles is constructed for the BPSO. A small number of particles in the swarm may result in local convergence, on the contrary, a large number of particles will increase computational efforts and may make slow convergence. Therefore, the s is usually chosen between 60 and 100. The initial position of the lth particle is identified by a binary matrix, denoted X 0 l , whose elements are randomly assigned by 1 or 0.

To adjust the position of particles, the initial velocities of all particles are also randomly generated (V 0 1 , . . . , V 0 s ). Velocity values are restricted to the limit values (V max , V min ) to prevent the particle from moving too rapidly from one region in search space to another. We usually set V max = 4, V min = -4. The initial velocity of the particle l, denoted V 0 l , is represented as matrix (22b).

Evaluate the swarm

The goal of this step is to determine the personal best and global best in the swarm. To do that, we must firstly evaluate the position of each particle l at iteration pth according to its fitness function which is defined as follows:

EP S(X p l ) otherwise The personal best, denoted ibest p l , is the best position associated with the best fitness value of the particle l obtained until iteration p. ibest p l is defined as follows:

, where q = 1 ÷ p

The global best at iteration p, denoted gbest p , is the best position among all particles in the swarm, which is achieved so far and can be expressed as follows:

, where l = 1 ÷ s Update the velocity of each particle The objective of this step is to calculate the new value of velocity of each particle to determine the next position of the particle. Based on the current velocity value, the personal best and the global best, the velocity of each particle can be updated as follows:

uniform random numbers between (0, 1); r 1 , r 2 are social and cognitive parameters which are taken as r 1 = r 2 = 2 consistent with the literature, see [START_REF] Franken | Investigating binary pso parameter influence on the knights cover problem[END_REF].

In BPSO, the velocities of the particles are defined in terms of probabilities that a element of matrix X p l will change to one. Using this definition a velocity must be restricted within the range [0, 1]. So whenever a velocity value is computed, the following piece-wise function, whose range is closed interval [V min , V max ], is used to restrict them to minimum and maximum value.

After applying the piece-wise linear function, the following sigmoid function is used to scale the velocities between 0 and 1, which is then used for converting them to the binary values.

Update position for all particles Based on new velocities obtained and let U is uniform random numbers between (0, 1), we update the new position for each particle by using the equation below:

Stopping criterion Stopping criterion is introduced to stop the algorithm process. Herein, limited iterations number (p max ) is used as a stopping criteria.

Phase 5: Rolling horizon

The maintenance manager can change the planning when some new information appears or when a planning for a new period and go back to step 4. Due to the previous phase, we have an optimal grouping structure within the finite planning horizon [t begin , t end ]. However, with time some new information (like maintenance resources constraints relying on for example management and new technology, opportunities, . . . ) by which this optimal grouping structure can be impacted may become available. To update the maintenance planning, a new optimal grouping structure within a new period must be identified. To this end, we simply go back to phase 2.


Phuc Do Van

Hai Canh Vu

Anne Barros

Christophe Berenguer

Grouping maintenance strategy with availability constraint under limited repairmen

Keywords: Grouping maintenance, availability, genetic algorithm, multi-unit system

The paper deals with a development of the rolling horizon approach introduced recently for grouping maintenance strategies of multi-component systems by integrating two efficient optimization algorithms. The proposed approach can help to construct an optimal maintenance planning with a given availability constraint under limited repairmen. Thanks to the rolling horizon spirit, the maintenance planning can be updated to take into account shortterm information which could be changed with time. A numerical example of a multi-component system is finally introduced to illustrate the use and the advantages of the proposed approach in the maintenance optimization framework.

INTRODUCTION

In the maintenance optimization framework, grouping strategy has been recently introduced for multi-component system with positive economic dependence which applies that combining maintenance activities is cheaper than performing maintenance on components separately, see for example [START_REF] Cho | A survey of maintenance models fo rmulti-unit systems[END_REF], [START_REF] Galante | An exact algorithm for preventive maintenance planning of series-parallel systems[END_REF], [START_REF] Radouane | Opportunistic policy for optimal preventive maintenance of a multicomponent system in continuous operating units[END_REF], [START_REF] Bouvard | Condition-based dynamic maintenance operations planning and grouping. application to comercial heavy vehicules[END_REF].

A major challenge of the maintenance optimization then consists in joining the stochastic processes regarding to the components (time-dependent probabilities of failure) with the combinatorial problems regarding to the grouping of maintenance activities. While a long term or infinite planning horizon can be assumed to solve this problem in case of stable situations, dynamic models have been introduced in order to change the planning rules according to short-term information (e.g. failures and varying deterioration of components), using a rolling (finite) horizon approach [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF], [START_REF] Nicolai | A review of multi-component maintenance models[END_REF]. This approach is however limited to a class of problems with significant assumptions: maintenance durations are neglected and only one preventive maintenance for each component is preventively maintained only is considered in horizon interval. To solve these problems, an extension of the rolling horizon approach is recently introduced in Do [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF]. Nevertheless in such papers, maintenance constraints have not yet been considered. From a practical point of view, it is often impossible to perform all the desirable maintenance actions due to the limitations on maintenance resources, such as maintenance budget, duration of maintenance time, and limited repairmen. The grouping maintenance problem remains widely open.

The aim of this paper is to propose an approach based on the rolling horizon spirit for grouping maintenance planning with availability constraint under limited repairmen. For example, according to a specific mission or a production campaign, the system may be required to operate with limited breaks durations, i.e. the maintenance duration time of the system have to be limited. With these constraints, to find optimal planning, dynamic programming presented in [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF], Do [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF] is no longer usable since the combinatorial problem can be formulated as a set partitioning problem, which however can be NP-hard. To this end, the Genetic Algorithm (GA), recognized as a general search strategy which is often useful for solving combinatorial problems, see [START_REF] Holland | Adaptation in natural and artificial systems[END_REF], [START_REF] Rahman | Performance analysis of genetic algorithm for solving the multiple-choice multidimensional knapsack problem[END_REF], is used. Moreover, the MULTIFIT algorithm introduced recently for bin-packing problems by [START_REF] Jr | An application of bin-packing to multiprocessor scheduling[END_REF] is also considered to determine the minimum maintenance duration time of each group under limited repairmen. This paper is organized as follows. Section 2 is devoted to the description of general assumptions and maintenance policies. Second 3 focuses on the development of the rolling horizon approach under maintenance constraints. Two optimization algorithms (GA and MULTIFIT) are presented in Section 4. To illustrate the proposed grouping maintenance strategy, a simple numerical example is introduced in Section 5. In addition some numerical results are discussed here. Finally, the last section presents the conclusions drawn from this work.

PROBLEM FORMULATION

General assumptions

Consider a series system consisting of n independent components. With this kind of systems, a preventive/corrective maintenance of one or more components leads to the system being unavailable. We assume here that the failure rate of component i (i = 1, ..., n) is described by a Weibull distribution with scale parameter λ i > 0, and shape parameter β i > 1, then

r i (t) = β i λ i ( t λ i ) βi-1 . (1)
We assume next that the deterioration of a component remains unchanged during the maintenance of other ones.

If a component fails, it is then immediately repaired. A corrective action restores the component involved into a state as good as before (minimal repair policy). Within this corrective policy, corrective maintenance durations are usually very small with respect to maintenance planning horizon and can be neglected. Therefore, when a corrective action is carried out for component i, it requires a corrective maintenance cost denoted C c i . We assume that preventive maintenance actions can be planned in advance to prevent failures or/and to reduce maintenance costs. Herein, each component is preventively maintained after a fixed interval to be optimized. After a preventive maintenance action, the maintained component is considered as good as new. The preventive maintenance cost of the component i can be divided in three parts:

• a specific component cost c p i ; • a setup-cost, denoted S, represents the preparation cost (or logistic cost) and can be shared when several maintenance activities are performed in the same maintenance occasion since execution of a group of maintenance activities requires usually only one setup; • an unavailability cost: since the system is unavailable during maintenance, therefore, if a preventive maintenance of component i, so-called activity i, leads the system to be unavailable during d i units of time, an additional cost or unavailability cost c d i = d i • C d relying on the production lost is incurred (C d is the unavailability cost rate of the system). This additional cost can be also shared when several maintenance activities are simultaneously carried.

As consequence, if a preventive maintenance of component i, namely activity i, is separately performed, we have to pay the following preventive cost:

C p i = S + c p i + c d i .
(2)

Based on the preventive/corrective maintenance costs and the reliability parameters of components, it is possible to define for each individual component i a nominal preventive maintenance periodicity x * i that minimizes the induced long-term average cost [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF], assuming that the maintenance activities are separately performed.

Maintenance constraints

If we consider an mission interval of time P H, the average availability of the system can be then calculated by:

A = P H -D P H , ( 3
)
where D is the total maintenance durations within this period or the total preventive maintenance durations since corrective maintenance durations are zero.

Assume now that the the average availability of the system is required to be higher a given level A 0 for this interval of time, i.e. A ≥ A 0 . From Equation (3), we have

D ≤ (1 -A 0 ) * P H or D ≤ D 0 , (4)
where D = (1 -A 0 ) * P H is the maximal total duration allocated for all preventive activities in the considered scheduling interval P H. To execute preventive activities, we suppose that that only m repairmen are available and each repairman can take only one preventive activity at a time.

Grouping maintenance

We are here interested in a grouping maintenance strategy according to two following reasons:

(1) it is shown in a number of papers, see for example [START_REF] Cho | A survey of maintenance models fo rmulti-unit systems[END_REF], [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF], that grouping can save maintenance costs since the set-up cost and the unavailability cost can be shared when several maintenance activities are executed together. Note well however that when several maintenance activities are performed together, maintenance cost could be indirectly penalized • with the reduction of components useful life if the maintenance dates are advanced; • with the increasing of components failure probability which could imply a system immobilization if the maintenance dates are too late. (2) according to the availability constraint, we have to execute simultaneously preventive activities to reduce maintenance durations even this grouping action may lead to a higher maintenance cost in some cases. For example, when the set-up cost and unavailability cost rate are equal to zero, grouping is then more costly but without grouping the required availability level may not be reached.

In order to find the optimal groups which could balance to minimize the system maintenance costs under given maintenance constraints (total maintenance duration allowed under limited repairmen) on the scheduling horizon, the rolling horizon approach introduced recently in Wildeman et al. [1997], Do [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF] will be developed and described in Section 3.

GROUPING MAINTENANCE STRATEGY

The developed grouping maintenance approach is divided into 4 phases:

• individual optimization;

• tentative planning;

• grouping optimization;

• update and decision.

Phase 1: Individual optimization

The objective is to find optimal individual preventive maintenance cycle based on an infinite-horizon maintenance model in which we assume an average use of components and in which the interactions between components are neglected.

Let M i (x) denote the expected deterioration cost for component i, i.e. the expected costs incurred in x time units since the latest execution of activity i. According to the minimal repair policy, M i (x) can be expressed as the following:

M i (x) = C c i • x 0 r i (y)dy, (5)
where r i (.) is the rate of occurrence of failure of component i. From equations ( 1) and ( 5), we obtain:

M i (x) = C c i • ( x λ i ) βi . (6) If component i is preventively maintained at x, the ex- pected cost within interval [0, x + d i ] is: Γ i (x) = C p i + M i (x) = C p i + C c i • ( x λ i ) βi . (7)
By using the renewal theory, the long-term average cost of component i can be determined as the following:

φ i (x) = Γ i (x) x + d i = C p i + C c i • ( x λi ) βi x + d i . ( 8
)
The optimal interval length for the preventive maintenances of the component i, denoted x * i , can be obtained when it leads to the minimum value of φ i (x).

x * i = arg min

x φ i (x) =⇒ dφ i (x) dx x=x * i = 0 ⇐⇒ C c i • (β i -1) • (x * i ) βi + C c i • β i • d i • (x * i ) βi-1 -C p i • λ βi i (x * i + d i ) 2 = 0 (9)
Let x * i be the optimal interval length for the preventive maintenances of component i when the maintenance duration is negligible (d i = 0), x * i can be calculated by (see [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF]):

x * i = λ i β i C p i C c i (β i -1) (10) Let g(x) = C c i • (β i -1) • x βi + C c i • β i • d i • x βi-1 - C p i • λ βi i . It is easy to show that g ′ (x) > 0, ∀β i > 1, ∀x > 0 ; g(0) = -C p i • λ βi i < 0 and g( x * i ) = C c i • β i • d i • ( x * i ) βi-1 > 0.
Then the equation ( 9) has only one solution (x * i ) in the interval (0, x * i ). Moreover, it is not difficult to show that φ ′′ i (x * i ) > 0. Therefore, we can deduce that x * i exists and it is unique,

x * i ∈ (0, x * i ).
The minimal long-term average maintenance cost of component i can be calculated as follows:

φ * i = φ i (x * i ) (11)
The optimal interval length x * i , which can be numerically calculated using Equation ( 9), represents a nominal preventive maintenance frequency of component i (with i = 1, ..., n) and can be used to define tentative execution times which help to identify optimal maintenance groups.

Phase 2: Tentative planning

The aim is to establish all tentative maintenance dates, assuming that maintenance activities are separately performed, within the interval of time P H.

Based on the nominal preventive maintenance frequencies, the first tentative maintenance execution time of component i (i = 1, ...n) from the current date denoted t 0 can be calculated by

t i 1 = t 0 -t e i + x * i + d Σ i , (12)
with t e i is the operational time elapsed from the last preventive maintenance of component i before t 0 , and d Σ i is the cumulative maintenance durations before the execution of component i. d Σ i is added since the system is stopped during maintenance.

It is shown in [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF], Do [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF] that to ensure that all components are take into account in the maintenance decision, the planning horizon should be chosen so that within this horizon interval each component is preventively performed at least one time. Therefore, if P H is shorter than (t j 1 + d j -t 0 ) with t j 1 = max i=1:n t i 1 , it is then set equal to (t j 1 + d j -t 0 ).

Since useful lifetime of components may be different, several components may be preventively performed more than one time. Let i j be the jth occurrence of maintenance activity i (or component i) in the scheduling horizon, the tentative execution time of operation i j , denoted t i j (j ≥ 2), depends on the executed time t * i j-1 of operation i j-1 , the cumulative maintenance durations d Σ i j-1 from t * i j-1 , and the nominal periodicity

x * i . t i j =t * i j-1 + d Σ i j-1 + x * i . ( 13
)

Phase 3: Grouping optimization

In the grouping spirit, each maintenance group is defined by three parameters:

• maintenance operations;

• group executed date;

• group maintenance duration.

The main idea is to find a grouping structure (or partition of all maintenance operations within the scheduling horizon) which minimizes the expected maintenance cost and scopes with an availability level constraint by using the Genetic Algorithm (GA). The MULTIFIT algorithm will be also used to determine the minimum maintenance duration of each group under limited repairmen, see Section 4.

Mathematical formulation

Penalty cost calculation If the jth (j = 1, ...) occurrence of maintenance activity i is actually executed at time t * i j = t i j + ∆t i j instead of the tentative execution time t i j , then the resulting additional cost is expressed by the penalty cost h i (∆t i j ) which is written as, see [START_REF] Cho | A survey of maintenance models fo rmulti-unit systems[END_REF], [START_REF] Wildeman | A dynamic policy for grouping maintenance activities[END_REF].

h i (∆t i j ) = Γ i (x * i + ∆t i j ) -Γ i (x * i ) + ∆t i j • φ * i ,
with ∆t i j > -x * i . By using Equation ( 7), we obtain:

h i (∆t i j ) = M i (x * i + ∆t i j ) -M i (x * i ) -∆t i j • φ * i ,
Another kind of penalty cost function can be found in [START_REF] Bouvard | Condition-based dynamic maintenance operations planning and grouping. application to comercial heavy vehicules[END_REF].

From equations ( 6) and ( 9), we obtain finally:

h i (∆t i j ) = C c i • x * i + ∆t i j λ i βi - x * i λ i βi -∆t i j C p i β i x * i (β i -1) . ( 14
)
Maintenance cost savings calculation Assume now that several different maintenance operations i j (i, j = 1, 2, ...) are performed in a group G k with k ∈ N . Remember here that operations i j and i l (j = l) are identical operations since they are respectively the jth and the lth occurrence of maintenance activity i. Let H G k (t) be the group penalty cost function at time t. The optimal execution time of the group t G k can be found when the H G k (.) searches its minimal value H * G k . That is:

H * G k = H G k (t G k ) = min t i j ∈G k h i (t -t i j ) ( 15
)
As mentioned above, the execution of a group of maintenance operations requires only one set-up cost, the group G k yields a cost reduction:

V 1 (G k ) = (card(G k ) -1) • S. (16)
An additional cost saving relying on the reduction of maintenance duration when maintenance activities are simultaneously executed by m repairmen.

V 2 (G k , m) = ( i j ∈G k d i -d G k (m))C d , (17)
where d G k (m) is the total duration of the group G k . Note well that d G k (m) depends on the duration of all group members and the number of repairmen m. d G k (m) will be optimally minimized by using the MULTIFIT algorithm, see Section 4.1.

The group savings is written as:

Q G k = V 1 (G k ) + V 2 (G k , m) -H * G k .
(18) Based on all different groups, a grouping structure, namely GS, can be identified. The total maintenance duration can be is determined by:

D GS = G k ∈GS d G k (m), ( 19
)
and the total cost savings can be calculated as follows:

Q GS = G k ∈GS Q G k . ( 20
)
According to availability contraint presented in Section 4, An optimal grouping structure can be searched by the following:

GS op = GS, max{Q GS }|D GS ≤ D 0 . ( 21
)
Note well that the maximum number of groups to be considered is here 2 n -1 for n maintenance activities. The finding of optimal grouping structure becomes difficult since the combinatorial problem can be formulated as a set partitioning problem, which however can be NP-hard.

To solve this problem, the GA algorithm will be used and described in Section 4.

Phase 4: Update and decision

Due to the previous phase, we have an optimal grouping planning within the finite planning horizon P H. However, with time new information by which the current planning can be impacted may occur, such as changing of maintenance constraints (increasing/decreasing of repairmen) or maintenance opportunities (see Do [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF]), ... which could lead the current planning to be no longer optimal or even unusable. A new planning is needed. To this end, we simply go back to phase 2.

OPTIMIZATION ALGORITHMS

MULTIFIT algorithm

Consider now a group consisting of k maintenance activities. To perform the maintenance of the group, only m (m ≤ k) repairmen are available. The objective is to find the minimum execution duration of this group. To this end, the MULTIFIT algorithm is used. A short description of this algorithm is the following:

Step 0 Without loss generality, k maintenance activities are renumbered by increasing order according to their decreasing duration, i.e.

d 1 ≥ d 2 ≥ . . . ≥ d k . Next, let dg up = max(d 1 , 2 k i=1 d i /m) and dg low = max(d 1 , k i=1 d i /m).
And finally, set w = 0.

Step 1 Set dg = 1 2 (dg up + dg low ) and w = w + 1.

Step 2 Assign successively maintenance activities to repairmen with time limited dg.

Step 3 If all activities can be assigned to m repairmen, then set dg up = dg and go to step 4, otherwise set dg low = dg and go to step 4.

Step 4 If w = w max then stop, otherwise go to step 1.

It is shown that the worst-case performance ratio is:

α = dg dg optimal = 13 11 + 2 -wmax
The iterations number can be fixed in advance, as recommended in [START_REF] Jr | An application of bin-packing to multiprocessor scheduling[END_REF], w max = 7.

Genetic algorithm

The Genetic Algorithm has been recognized as a general search strategy which is often useful for solving combi-natorial problems [START_REF] Holland | Adaptation in natural and artificial systems[END_REF]. Recently, GA has been widely studied, experimented and applied in many fields in engineering worlds. A GA starts by creating an initial population of solutions. Each solution is next evaluated using an objective function. During each iteration step, genetic operations (crossover and mutation operations) are applied in order to search potential better solutions. The GA structure is in presented in Fig. 1. Coding A coding step determines how the problem is structured in the algorithm. In our problem, each solution is represented by an matrix, namely X, with ngmax rows and n columns. Where, ngmax is the maximum number of groups in a feasible solution. ngmax depends on the requested availability level A 0 , or the allowed total maintenance duration D 0 . If N (N ≥ n) preventive maintenance activities with durations d 1 ≥ d 2 ≥ ... ≥ d N are considered in the scheduling interval, the maximum number of groups can be then determined as the following:

Calculation of fitness function

ngmax = {l + 1|D 0 -d 1 ≥ l i=1 d N -i+1 }. (22)
If maintenance activity j is in group i, X ij = 1 otherwise X ij = 0. Here, X corresponds to a grouping structure.

Example: considering 7 preventive activities with

d 1 = 5, d 2 = 1, d 3 = 1, d 4 = 1, d 5 = 1, d 6 = 1, d 7 = 1.
If the total maintenance duration allowed is 8, by applying ( 22), the maximum number of groups is then 4. A feasible solution containing 4 groups

G 1 = {1, 2, 7}, G 2 = {3},G 3 = 4, G 4 = {5,
6} is represented by matrix X as follows:

X =    1 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0   
Generation of initial population GA makes randomly an initial population of solutions. If there are too few solutions in the population, GA have a few possibilities to perform crossover and only a small part of search space is explored. On the other hand, if there are too many solutions, GA slows down. The population size is usually chosen between 60 and 100 solutions. To generate an initial solution, the number of groups in the solution is randomly chosen between 1 and ngmax. Next, all maintenance activities are randomly put into the chosen groups.

Calculation of fitness function Each solution is evaluated by its fitness function which is defined as the following:

f (X) = Q X if D X ≤ D 0 -∞ otherwise (23)

Selection to crossover

The goal is to select pairs of parent solutions to Crossover phase. It is divided into 2 steps:

• Solution selection: the selection process is based on fitness. Array of parent solutions are sorted according to their fitness values in ascending order. The population is categorized into s groups. Next, a parent solution is randomly chosen from these s groups with s groups probabilities p 1 ≤ p 2 ≤ ... ≤ p s ( s i=1 p i = 1) respectively. Thus, the fitter a solution is, the more chance it have to be chosen for a parent in the next generation. This combinatorial grouping was found most effective;

• Selection to Crossover: each pair of selected solutions is now randomly selected for Crossover step with Crossover probability that is introduced to leave some part of population survive to next generation. Crossover probability is usually chosen between [70% ÷ 90%]. If a pair of selected solutions is not selected for Crossover phase, it is moved directly to Mutation step.

Crossover

The objective is to combine selected solutions to generate next better generation solutions for Mutation step by preserving their characteristics. Single point crossover is used in this algorithm which randomly chooses a locus and exchanges the subsequences before and after that locus between two selected solutions to create two children. For example, in Fig. 2, given two parent solutions we choose randomly maintenance activity 4 as the crossover point, i.e. solution 1= {Ia,Ib}, solution 2={IIa,IIb}. As consequence, two new solutions are {Ia,IIb} and {IIa,Ib}.

Mutation Mutation is made to prevent falling GA into local extreme, but it should not occur very often, because then GA will in fact change to random search. So, we usually choose Mutation probability in a interval [0.2% ÷ 1.5%].

For each selected solution, a random maintenance activity in a group is next moved to another group in order to generate a new solution. An illustration of the mutation process is shown in Fig. 3.

Final condition Final condition is introduced to stop the algorithm process. Herein, limited generations number is Finding optimal solution the optimal solution with highest fitness value is searched among the last generation solutions.

NUMERICAL EXAMPLE

The purpose of this section is to show how the proposed grouping strategy can be used in preventive maintenance optimization with availability constraint under limited repairmen through a simple of 20 series components system in which each component is preventively replaced. When a component fails, it is immediately maintained according to a minimal-repair policy. A failure repair restores the component involved into a state as good as before. We assume failure rate of a component i (i = 1, ..., 20) is described by a Weibull distribution with scale parameter λ i > 0, and shape parameter β i > 1. Table 1 reports the random data for five components. For set-up cost and unavailability cost rate, we take S = 10 and C d = 5 respectively.

By assuming that all components are individually maintained, nominal maintenance periodicity x * i , minimum average maintenance cost φ * i and the next preventive replacement date t i 1 (with i = 1, ..., 20) are calculated by substitution of the input data in Equations ( 9) and ( 12). The results are given in Table 2.

Maintenance planning without availability constraint

Consider firstly a grouping maintenance planning without availability constraint using the rolling horizon approach. Only one repairman is considered for all maintenance activities. To define the scheduling horizon, the current date is set equal to 0 (t begin = 0) and the ending horizon date corresponds to t 20 1 + d 20 = 606 (t end = 606), i.e. P H = 606. Herein, 20 individual preventive activities are therefore considered in this interval of time. To find the optimal grouping planning, the dynamic programming presented in Do [START_REF] Van | Dynamic grouping maintenance strategy with time limited opportunities[END_REF] was used and the results are reported in Table 3. The total maintenance cost savings is Q * Σ = 148.64. And the total maintenance durations is D = 14 + 26 + 18 + 15 = 73, the average availability of the system is hence A = (P H -D)/P H = 0.8795.

  Fig. 1. Structure of GA

FigFig. 3 .

 3 Fig. 2. Illustration of Crossover step

Table 1 .

 1 Data of 20 components system.

	Unit	λ i	β i	c p i	C c i	d i	t e i
	1	288	1.9334	168	42	1	638.98
	2	281	1.9797	159	47	2	508.91
	3	276	1.9965	155	47	6	486.87
	4	264	1.9930	111	36	2	405.53
	5	288	1.9244	142	38	3	538.48
	6	270	1.9209	108	59	4	252.64
	7	265	1.9417	127	46	3	313.84
	8	278	1.9271	180	31	5	580.86
	9	294	1.9992	180	50	4	378.94
	10	265	1.9098	159	58	3	259.66
	11	265	1.9223	106	57	2	161.57
	12	286	1.9731	123	51	5	239.35
	13	268	1.9246	123	38	5	275.31
	14	291	1.9777	117	41	1	206.62
	15	252	1.9480	103	57	5	53.79
	16	291	1.9196	146	44	2	248.10
	17	256	1.9585	158	53	5	129.09
	18	277	1.9358	134	58	5	56.08
	19	267	1.9859	171	58	4	15.83
	20	266	1.9509	187	31	6	224.36

Table 2 .

 2 Values of x * i , φ * i , and t i 1 .

	Unit	x * i	φ * i	t i 1	Unit	x * i	φ * i	t i 1
	1	638.98	0.5932	0	11	417.57	0.6289	289
	2	557.91	0.6483	50	12	514.35	0.6229	310
	3	563.87	0.6929	80	13	585.31	0.5619	350
	4	506.53	0.5191	110	14	531.62	0.5023	370
	5	647.48	0.5369	120	15	407.79	0.6954	400
	6	438.64	0.6562	200	16	607.10	0.5708	410
	7	505.84	0.6196	210	17	506.09	0.7792	430
	8	789.86	0.5658	230	18	498.08	0.7019	500
	9	602.94	0.6969	250	19	502.83	0.8052	550
	10	509.66	0.7578	280	20	757.36	0.6149	600

Table 3 .

 3 Grouping without constraint.

	Group units	Optimal date t G	duration d G
	{1,...,5}	67.16	14
	{6,...,12}	242.96	26
	{13,...,17}	389.25	18
	{18,...,20}	538.56	15

ACKNOWLEDGMENTS

This work has been supported and partly financed by the European project MoDe -Maintenance on Demand (http://fp7-mode.eu) funded from the European Community's 7th Framework Program (Project SCP8-GA-2009-233890 -FP7 Sustainable Surface Transport).

We are interested now to construct a grouping maintenance planning with a given availability level under limited repairmen. To find optimal groups, the dynamic programming approach is no longer applicable and the GA algorithm presented in Section 4 is used. GA simulation was generated with 1500 generations from 60 initial solutions. The MULTIFIT algorithm is also used to find the minimum duration of each maintenance group with limited repairmen. Different values of the availability level constraint A 0 and the number of repairmen m are considered for the scheduling interval [0,606]. The results presented in Table 4 show that:

• to establish a maintenance planning under a given availability level constraint, the number of necessary repairmen is bounded. The lower bound corresponds to the minimum number of repairmen with that a maintenance planning can be constructed. The upper bound corresponds to the minimum number of repairmen that leads to a maintenance planning with highest cost savings. For example, when A 0 = 0.9719 (the system availability is required to be not lower than 0.9719), the lower and upper bound are 5 and 7 respectively. If the number of repairmen is lower than 5 (m < 5), no maintenance planning can be then established. Maintenance planning and cost savings remain unchanged when m ≥ 7. Note well that the upper bound of repairmen is lower than the maximum number of maintenance activities in a group, see for example, when A 0 = 0.9901) we need only 13 repairmen for executing a group of 20 maintenance activities with a minimum execution duration; • when the number of repairmen is in the bounded interval, an increasing of repairmen leads to a maintenance planning with an increasing of maintenance cost savings.

The results show also that when only single repairman is considered for preventive actions and the availability level constraint is low enough, the GA algorithm and the dynamic programming lead to the same grouping maintenance planning, see for example the case with A 0 = 0.8795 and m = 1.

CONCLUSIONS

In this work the rolling horizon approach introduced recently is developed by integrating both GA and MUL-TIFIT algorithm in order to establish an optimal maintenance planning of multi-component systems under an availability level constraint with limited repairmen. The proposed approach can also help (i) to determine the minimum number of repairmen to ensure that a establishable maintenance planning, which scopes with a given availability level constraint, can be constructed; (ii) to determine the minimum number of repairmen that leads to an optimal maintenance planning satisfying the availability constraint with highest cost savings. Moreover, thanks to the rolling horizon spirit, the grouping maintenance planning can be easily updated when short-term information is available, e.g an increasing/decreasing of repairmen or opportunities. Our future research work will focus on the development of the proposed approach for systems with inter-component dependencies. Furthermore, condition based maintenance will be also developed in this grouping maintenance framework.