Emmanuel Jeannot
email: emmanuel.jeannot@inria.fr

Guillaume Mercier
email: guillaume.mercier@bordeaux-inp.fr

François Tessier
email: ftessier@anl.gov

Topology and affinity aware hierarchical and distributed load-balancing in Charm++

come   L'archive ouverte pluridisciplinaire

Introduction

Scientific applications needs for computing power are continuously increasing. In order to match these needs, larger and larger parallel computers are designed and built. Optimizing the use of the underlying physical resources is an essential objective in order to take advantage of the current complex multicore architectures of the nodes. One other issue stems from the fact that parallel applications have to deal with problems featuring irregular structures that can easily induce a load imbalance at the CPU level. Such an imbalance can in return yield counterproductive effects on the application performance. For instance, a process that has finished its work will become idle if it has to wait (synchronize) for another more loaded process still computing. In such a case, this idleness leads to a considerable waste of CPU cycles. This issue is even more prominent as the 1 number of processes grows. Therefore, balancing the load among the processes is an efficient and relevant optimization technique to decrease an application execution time.

However, load balancing can be detrimental to data locality since it involves processes/tasks movements and migrations across the computing platform. In case of internode movements, it implies network communications whilst in case of intranode movements, Non-Uniform Memory Access (NUMA) effects are likely to degrade performance. Since data locality is a key-factor for large-scale applications performance, it should also be taken into account by load balancers.

Improving or managing efficiently data locality is not an easy task, especially in the context of large clusters of multicore nodes. Indeed such computers feature intricate and hierarchical memory systems with several cache levels and banks physically scattered across the whole machine. This trend is going to gain momentum with the forthcoming advent of new memory technologies (e.g. NVRAM) because main memory will now be made of several levels of different memory types. Data locality can be managed in various fashions, one being to create a match between the physical topology (that is, the CPUs and memory hierarchy layout) and a virtual topology that describes the application bahaviour. For instance, an application communication pattern (in case of processes) or the pattern of memory accesses (in case of threads) are pertinent instanciations of virtual topologies.

We developed an algorithm, called TreeMatch [START_REF] Jeannot | Process Placement in Multicore Clusters: Algorithmic Issues and Practical Techniques[END_REF] that computes such a matching between the physical entities delivering the computing resources and the virtual entities (e.g processes or threads or tasks) handling the data and performing the computation. TreeMatch uses a qualitative approach as it relies on the structure of the underlying physical topology. It has been successfully integrated within MPI libraries and runtimes in order to perform MPI process placement [START_REF] Jeannot | Near-Optimal Placement of MPI processes on Hierarchical NUMA Architectures[END_REF] and rank reordering [START_REF] Mercier | Improving MPI Applications Performance on Multicore Clusters with Rank Reordering[END_REF]. Therefore, an intuitive idea is to expand load balancing algorithms and mechanisms with TreeMatch in order to make them data locality-aware.

For instance, the Charm++ parallel programming environment [START_REF] Kale | CHARM++: A Portable Concurrent Object Oriented System Based on C++[END_REF] natively features load balancing mechanisms. In Charm++, the load is divided among chares that are able to migrate across the platform. Hence, migration is used to balance the load, as opposed to work-stealing (for instance). Another interesting and useful feature of Charm++ is its ability to provide the user with various pieces of information during an application execution. For instance, information dealing with the workload evolution or the communication pattern can be exploited.

The contribution of this paper is a scalable parallel and distributed load balancer fully topologyaware (both at the nodes and at the network levels) that takes into account the computation as well as the communication cost of the application.

This paper is organized as follows: a state of the art is given in Section 2. The topology mapping problem is exposed in Section 3. The proposed load balancer is detailed in Section 4, while Section 5 describes the experiments carried out and the results obtained. Finally, Section 6 concludes this paper.

State of the Art

Charm++ [START_REF] Kale | CHARM++: A Portable Concurrent Object Oriented System Based on C++[END_REF][START_REF] Kale | Charm++: Parallel Programming with Message-Driven Objects[END_REF][START_REF] Laxmikant | Charm++ and AMPI: Adaptive Runtime Strategies via Migratable Objects[END_REF] is a message-passing based programming environment that uses an objectoriented approach and relies on the C++ language. However, whilst MPI considers processes in its model (usually with a rather coarse granularity), Charm++ rather uses finer-grain objects and dispatches the computation onto small migratable tasks called chares. These chares, in addition to their assigned data and their ability to exchange messages asynchronously, are also characterized by their CPU load, their I/O communication volume and some other useful parameters. To manage these processing entities, an adaptive runtime support is essential [START_REF] Laxmikant | Programming Models at Exascale: Adaptive Runtime Systems, Incomplete Simple Languages, and Interoperability[END_REF]. The Charm++ runtime system is able to handle the various physical resources, to enforce fault tolerance or manage load balancing [START_REF] Kale | Migratable objects + active messages + adaptive runtime = productivity + performance a submission to 2012 HPC class II challenge[END_REF][START_REF] Kale | Charm++ for productivity and performance: A submission to the 2011 HPC class II challenge[END_REF]. Regarding the load balancing aspects, Charm++ makes it possible to design, plug in and test load-balancers transparently without modifying the application code. Also the design of the Charm++ runtime system allows to perform load balancing at different frequencies during the execution of the application [START_REF] Brunner | Handling application-induced load imbalance using parallel objects[END_REF][START_REF] Zheng | Achieving high performance on extremely large parallel machines: performance prediction and load balancing[END_REF].

Several generic load balancing mechanisms natively exist in Charm++ [START_REF] Rodrigues | A Comparative Analysis of Load Balancing Algorithms Applied to a Weather Forecast Model[END_REF], but are tailored and give the best results for a specific optimization area. Common load balancing schemes that consider the CPU load on each processing unit have been extended in some works to also take into account the topology of the underlying architecture [START_REF] Bhatele | Topology Aware Task Mapping[END_REF], [START_REF] Bhatelé | Benefits of Topology Aware Mapping for Mesh Interconnects[END_REF]. For instance, Bhatelé and al. [START_REF] Bhatelé | Benefits of Topology Aware Mapping for Mesh Interconnects[END_REF] show the gains obtained by applying a chare placement depending on the structure of a 3D torus in the context of an adaptive mesh refinement applications. Hierarchical load balancing (i.e. multilevel topologies) has also been investigated [START_REF] Zheng | Hierarchical Load Balancing for Charm++ Applications on Large Supercomputers[END_REF] [START_REF] Zheng | Periodic Hierarchical Load Balancing for Large Supercomputers[END_REF] while other methods are designed to be distributed in order to achieve a satisfactory level of scalability on large-scale platforms [START_REF] Menon | A distributed dynamic load balancer for iterative applications[END_REF]. NucoLB1 [START_REF] Laércio L Pilla | A Hierarchical Approach for Load Balancing on Parallel Multi-core Systems[END_REF] and HwTopoLB [START_REF] Laercio L Pilla | Asymptotically optimal load balancing for hierarchical multi-core systems[END_REF] apply a load balancing scheme based on a quantitative assessment of the topology links (latency and bandwidth values are mandatory). As a consequence, this type of strategy requires to gauge the target architecture communication performance with the appropriate tools before running any application. Our TreeMatch-based solution is more flexible and dynamic since we fully rely on a qualitative approach for our representation of the hardware topology. TreeMatch does not require to assess the performance of the system on which the application is running. We believe that this is a strong advantage, as gathering such information is error-prone, might be incomplete and is subject to inaccuracy.

We already worked on a topology-aware load balacing algorithm in Charm++ using Tree-Match [START_REF] Jeannot | Communication and Topology-aware Load Balancing in Charm++ with TreeMatch[END_REF] called TMLB TreeBased. However, this load-balancer suffered from critical limitations (handled topology, scalability, load balancing techniques, etc.). The completely redesigned algorithm described in this paper is no longer limited to tree-shaped networks which is undeniably a more realistic approach. Most of the algorithm is distributed and hierarchical, improving significantly the capacity of scaling while the distribution of the work is parallelized in a master-worker fashion. Besides, we modifed the part of the algorithm in charge of the CPU load balancing into a tunable refinement approach.

Context and Problem Definition

We consider a Charm++ application with both load balancing and migration capabilities enabled. The load balancing phase is executed periodically by the Charm++ runtime system. The application is composed of multiple processes that are bound on their dedicated CPU cores. Each process executes a given set of Charm++ chares. The load balancer possesses the following pieces of information. [START_REF] Aochi | Finite difference simulations of seismic wave propagation for the 2007 mw 6.6 niigata-ken chuetsu-oki earthquake: Validity of models and reliable input ground motion in the near-field[END_REF] The load of the different chares (keep in mind that the number of chares is greater than the number of processes). This load actually represents the amount of CPU used by each chare since the last load-balancing phase. [START_REF] Bhatele | Topology Aware Task Mapping[END_REF] The affinity between chares. This affinity is represented by a communication matrix which values are the amount of communication between a pair of chares since the last load-balancing phase. (3) The system topology. Here, we deal with two kinds of topologies: first, the node internal topology, where the memory hierarchy is represented by a tree and second, the network topology which forms an arbitrary graph where each vertex represents an allocated compute node for the application.

Both the load and the affinity of the chares are computed by the Charm++ runtime system and are available in the load-balancing code through the Charm++ API as a set of standard arrays or matrices. The node topology is retrieved thanks to the hwloc tool [START_REF] Broquedis | Hwloc: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF], following the assumption that all nodes feature a similar topology across the supercomputer (which is the case in most of today's top-end HPC systems). The network topology is obtained either through the batch scheduler or a standard system service that can be queried by the application to collect the graph of allocated resources (i.e. compute nodes). For instance, on the BlueWaters system, such information is obtained thanks to the xtdb2proc command, provided as part of the Cray environment of the machine.

Based on this information, our goal is to compute a new chares allocation for the processes so that the application overall execution time is decreased2 . Experiments carried out in related works [START_REF] Laércio L Pilla | A Hierarchical Approach for Load Balancing on Parallel Multi-core Systems[END_REF][START_REF] Laercio L Pilla | Asymptotically optimal load balancing for hierarchical multi-core systems[END_REF][START_REF] Jeannot | Communication and Topology-aware Load Balancing in Charm++ with TreeMatch[END_REF] have shown that taking solely into account the load to perform the load balancing leaves room for improvement. Indeed, due to both the network topology and the memory hierarchy, the communication costs between processes or chares have also to be considered. As a consequence, in order to tackle the issue of load balancing, this work considers not only the load but also both the affinity between the various computing entities of the application and the underlying target hardware topology.

4 Scalable Topology Aware Load Balancing at System Scale

Overview

In order to solve our problem we need to design a scalable load-balancing algorithm that takes into account the chares load, the chares affinity and the topology of the target machine. The chares distribution is computed at each load-balancing step. This computation is two-fold: globally at the nodes level and then inside each node, chares are assigned to cores in parallel. For the nodes step, we balance the CPU load at this level of hierarchy first and then we move groups of chares on nodes to reduce the communication costs. For the cores step, we first compute an affinity-aware chares placement then we apply a refinement algorithm to level the CPU load.

On the Need of a Global Load-balancing Phase

We first discuss the need for a global chare allocation and load-balancing step. One could devise a fully distributed algorithm with no global load-balancing phase (that is, chares would stay on their nodes and would only move from one core to an other). But the following theorem shows that the probability can be extremely high that at least one pair of two nodes has an imbalance factor greater than a given ratio.

Theorem 1 Let m be the number of nodes of a run. Let n be the number of chares per node. Assuming that chares have a duration time uniformly distributed in [a, b] then the probability q that there is two nodes which imbalance ratio is greater than α is

q = 1 -(1 -p) m(m+1) 2 where p = 2 -2Φ α √ 6n 2
b+a b-a and Φ is the cumulative distributed function (CDF) of the standard normal distribution N (0, 1).

Proof.

Let us consider that we have n chares per node. Let us suppose that the duration time of the chares is uniformly distributed between [a, b]. Hence, by definition of the uniform distribution, the mean duration time of the chares is: µ = b+a 2 and the variance of the duration time of the chares is:

σ 2 = (b-a) 2
12 . On a given node, the sum of the duration times of the chares has a mean of µ s = nµ and a variance of σ 2 s = nσ 2 . If n is large enough3 , we can apply the central limit theorem to compute the distribution X i of the sum of the duration times of the chares on node i. Assuming that the duration time of the chares is independent and identically distributed, then, X i does not depend on i and, by the central limit theorem,

X i ∼ N (µ s , σ 2 s ) = N n(b + a) 2 , n(b -a) 2 12
where N is the gaussian (normal) distribution and X ∼ D means that the random variable X follows distribution D. Let us first answer the following question: what is the probability p that the ratio of duration times of all the chares of two nodes (i and j) differs by more than a given α ≥ 0? This can be written:

p = P X i -X j nµ ≥ α = P(|Z| ≥ α) where Z = X i -X j nµ = X i -X j n(b+a) 2
From the rule of normal distribution difference, we have:

X j -X i ∼ N (0, 2 n(b-a) 2 12
). Then, thanks to the scaling property of the gaussian distribution 4 ,

Z ∼ N 0, 2n(b-a) 2 12 4 n 2 (b+a) 2 = N 0, 2(b-a) 2 3n(b+a) 2 . Moreover, if we put: Y = Z 2 √ 6n b-a b+a ∼ N (0, 2(b-a) 2 3n(b+a) 2 6n(b+a) 2 4(b-a) 2 ) = N (0, 1)
, we can simplify the solution even further:

p = P(|Z| ≥ α) = P |Y | ≥ α √ 6n 2 b + a b -a = 1 -P -α √ 6n 2 b + a b -a < Y < α √ 6n 2 b + a b -a )
Therefore, as Y is the standard normal distribution we have:

p = 1 -Φ α √ 6n 2 b + a b -a -Φ -α √ 6n 2 b + a b -a
Where Φ is the cumulative distributed function (CDF) of N (0, 1). As N (0, 1) is symmetric around

0 we have Φ(-x) = 1 -Φ(x) therefore, p = 2 -2Φ α √ 6n 2
b+a b-a and α ≥ 0. Interestingly, since Φ is increasing, we see that p is decreasing (as expected) with n and α. Note also that even if p depends on a and b that if b = r × a (b is a multiple of a then, b+a b-a = r+1 r-1 is constant.

However, p is simply the probability that two nodes have a load ratio greater than α. If we call q the probability at least one pair of nodes (among the m used in the run) have a load ratio greater than α we have:

q = 1 -(1 -p) m(m+1) 2
Indeed we have m(m+1) 2 pairs of nodes and (1 -p)

m(m+1) 2
is the probability that all pairs have a ratio lower than α.

Table 1 outlines the results for different realistic scenarios. We see that in most of the cases, there is a high probability that at least two nodes are highly imbalanced. This advocates for the use of a global balancing step where we allow to move any chare to any node. Note that the first line of the Table 1 corresponds to an experimental case of this paper.

Computing the Chares Placement at Nodes Level

As shown above, global chare placement is often mandatory. Such a step is based on LibTopoMap which is able to compute a mapping on any arbitrary topology. To do so, we first need to compute the node topology (based on the job allocation) and then call LibTopoMap to compute the new allocation.

Computing the Allocated Nodes Topology

In order to balance the workload by taking into account both the network topology and the memory hierarchy and layout, we need to compute the topology formed by the nodes allocated by the batch scheduler. Then we apply to these nodes a local numbering scheme (starting from zero). In the case Figure 1: Cray Gemeni topology management of a fat-tree network, the nodes numbering is performed in topological order and the network spans the part of the tree used by the allocated nodes. If not, building the topology of the allocated nodes requires to extract the global coordinates of the nodes as well as their interconnections and then to renumber the nodes in a topological order. For instance, the BlueWaters Cray Gemini interconnect topology is a 3D torus. However, some nodes in the torus are only routers and compute nodes are connected to the routers as there are two compute nodes per router (c.f. Fig. 1a). We proceed as follows: (1) Extraction with the Cray tools of the smallest 3D grid of routers including all the compute nodes allocated by the resource manager. Each router features a set of (x, y, z) coordinates in the 3D torus; (2) Renumbering of the routers in topological order (x, then y, then z); (3) Creation of an adjacency list of every routers of the grid; (4) Creation of the link between each compute node and its router; (5) Exclusion of the nodes not part of the resource allocation; (6) Output a file in a LibTopoMap-readable format.

As opposed to TreeMatch, LibTopoMap relies on a quantitative approach. To fulfill this requirement, we defined the weights of the topology edges as follows : the x-axis and z-axis weights are twice as large as the y-axis'. This choice is motivated by a particular feature of the Cray Gemini interconnect routers as illustrated by Figure 1a. All axes possess two interconnection links except for the y-axis where one link is devoted to the connection of compute nodes.

Figure 1b depicts how this step works. We can see a 3D grid which includes the allocated routers (blue vertices). The compute nodes are not shown but are indeed connected to these routers temporarily renumbered following the (x, y, z) order. Then, the corresponding adjacency list is given to LibTopoMap to compute the relevant placement.

Global chares allocation

As explained previously, the global chares allocation (i.e. assigning balanced groups of chares on nodes) is a two-step process.

From theorem 1 it is likely that at least two nodes are highly imbalanced. Therefore, we first smooth the load so that the imbalance factor α between two arbitrary nodes is lower than 5%. To do so, we perform a very simple set of refinements of the load by swapping the highest loaded chares from overloaded nodes for the lowest loaded one from underloaded nodes. After this step, all the nodes roughly possess an equivalent load. Of course, this step is not performed if the load is already evenly distributed among the different nodes.

The second part consists in moving groups of chares by taking into account their affinity and the communication cost of the network. We use the LibTopoMap [START_REF] Hoefler | Generic Topology Mapping Strategies for Large-Scale Parallel Architectures[END_REF] library during this step because of its following pros: it works in a distributed fashion (as an MPI-based library) and it is able to find an efficient placement of MPI processes on a cluster of compute nodes regardless of the network topology (that is, not only tree-like topologies).

This library makes use of the virtual topologies management features of the MPI standard to determine this process placement. In the load-balancing algorithm presented in this paper, we use LibTopoMap in order to perform a first placement of the groups of chares on the nodes. Since Charm++ also uses MPI in its internal communication layer, we can access to the native MPI functions. As LibTopoMap only needs one MPI process per node to compute its network-aware placement, we select a rank on each node and create a dedicated communicator with MPI Comm split. LibTopoMap considers the node communication matrix. This matrix describes the affinity between groups of chares that are already assigned to a given node. LibTopoMap uses Metis [START_REF] Karypis | METIS -Unstructured Graph Partitioning and Sparse Matrix Ordering System[END_REF] to assign groups of chares to the nodes according to the topology computed in Section 4.3.1 and the communication matrix.

Fine Placement of the Chares on the Cores

In the previous step, groups of chares are migrated from one node to an other according to their load, their affinity and the topology. In this step we assign chares to cores on each node, according to the same criteria but in a distributed fashion.

In this step, we need each core to be assigned an equal load and that the communication between chares is taken into account as well as the node topology. For taking into account the topology, we use the hwloc tool which provides a simple way to apprehend the memory and cache hierarchy as well as the core numbering. Remind that the affinity between chares and their load is provided through the Charm++ API. step is needed. Indeed, each node can perform its own step independently from the others and in parallel. To distribute the load among all the nodes, we use a master/worker scheme. The master algorithm is described in Algorithm 1 and features two phases: the first one is to decompose the global communication matrix into local communication matrices to be used by each node. To do so, we extract the element of the m global chares global communication matrix that corresponds to the chares of node i. This is done thanks to the knowledge of the current mapping of the chares M. As a master node possesses in general several cores, we use these cores to perform this phase in parallel. To do so, we create a parallel OpenMP section that extracts the communication pattern of the group of chares of each compute node (lines 2 to 4). The second phase is to distribute the chares mapping of each node to the corresponding node. Therefore, each node computes the placement of its own set of chares. To be more specific, the call to compute lb on line 5 is distributed thanks to Charm++ mechanisms and each node computes its chares placement in parallel. The call to compute lb is done asynchronously. When each local result of node i is returned (line 7), it is aggregated to compute the final global placement of the chares (line 8). This global scheme is depicted in Fig 2.

Worker Algorithm

Each node receives the communication matrix corresponding to its own set of chares and computes the placement of these chares on its own cores. The worker algorithm is depicted in Algorithm 2 and is executed on each node, including the master one. Each worker knows the load and the affinity of the chares. The first step of the algorithm consists in getting the topology tree T of the node. Each leaf of T thus corresponds to a core. The problem is that generally we have much more chares than cores. We therefore need TreeMatch to perform the assignment with oversubscribing (line 3). To this end, we compute the ratio r between the number of chares (n) and the number of cores (c). We extend the tree T by adding a new level to the tree with r leaves to each core of T , leading to a tree T with c = r * c leaves. TreeMatch then assigns each chare to a leaf of T . In order to balance the number of chares to cores, we constraint the TreeMatch algorithm so that chares are evenly distributed.

Then the solution p (line 3) is such that chares assigned to leaves (of T ) corresponding to the same core (of T ) are actually assigned to this core. In this case we have at most r chares per core of the node. Then solution p is refined (line 4). The goal of this step is to take into account the load imbalance between different cores. This refinement step works as follows: chares are topologically sorted according to the core they have been assigned to. Then cores are considered one after the other: if a considered core is the least loaded one, then it steals other chares starting from the closest leaf until it is no longer the least loaded one. In this way, we keep the chares that communicate a lot close to each other and the above refinement allows for a good repartition of the load among the cores. Every node returns its result back to the master so it can build a global view of the chares mapping and commit the chare migration using internal Charm++ features. No actual migration occurs before this point.

Experimental Validation

In this section, we present the results obtained with our load balancing algorithm. We chose three applications to evaluate our solution. The first one, called ChaNGa, is a cosmological application designed to perform collisionless N-body simulations. The second application, Ondes3D, is a simulator of three-dimensional seismical waves propagation. The third one named CommBench is a benchmark simulating irregular communications.

ChaNGa

The first selected application is a cosmological simulator. Its main goal is to perform collisionless N-body simulations while using a Barnes-Hut tree to calculate gravity [START_REF] Jetley | Massively parallel cosmological simulations with ChaNGa[END_REF][START_REF] Jetley | Scaling hierarchical n-body simulations on gpu clusters[END_REF][START_REF] Menon | Adaptive techniques for clustered n-body cosmological simulations[END_REF]. This application has been written in Charm++ in order to exploit the Charm++ features and particularly its ability to migrate chares to balance the CPU load. As a use case, we worked on Lambb [START_REF] Menon | Adaptive techniques for clustered N-body cosmological simulations[END_REF] which is an 80 million particles representation of a 70 Mpc (Megaparcec) volume. According to the ChaNGa documentation, this simulation is used to calculate the mass function of dark matter halos in a dark energy dominated universe down to the scale of dwarf galaxies. The experiments on this real application were carried out both on the PlaFRIM cluster and the Blue Waters supercomputer. For each experiment, we compared our load balancing algorithm to an execution without load balancing (NoLB) and to other load balancers. RefineLB and GreedyLB are standard load balancers available within Charm++. RefineLB migrates chares from overloaded cores to underloaded cores to reach an average load. This strategy has the advantage to limit the number of chares migrations. GreedyLB re-assign all the chares by mapping the highest loaded chare to the least loaded core. Even if these load balancers were not designed for this particular use case, we choose to compare our solution with Orb3dLB notopo (based on recursive bisection to find a balanced stat) and MultistepLB notopo (considers a range of different timesteps to take its decision). Indeed, these strategies were developed for specific use cases of ChaNGa and help us to show that TreeMatchLB is a generic, yet relevant approach. It has to be noted that the load balancing time is included in each measurement of these experiments.

PlaFRIM

The PlaFRIM cluster is composed of nodes featuring two Quad-core-Intel Xeon Nehalem X5550 (2.66 GHz) processors. Each node offers 24 GB of 1.33GHz DDR3 RAM and the 8 MB of L3 cache are shared between the four cores of a CPU. Nodes are connected through an Infiniband fat-tree network. We ran experiments on 64, 128 and 256 cores (respectively 4, 8 and 16 nodes) of the PlaFRIM platform.

The results we obtained on this platform are depicted in Figure 3. This figure shows the average execution time of one step of the Lambb case for each number of cores used. The first remark we can make is that TreeMatchLB yields a performance level that is at least on par with the other strategies. On 64 cores, TreeMatchLB is the fastest approach while for the 128 and 256 cores cases, it competes equally with Orb3dLB notopo and MultistepLB notopo. As far as the load balancing time is concerned, our solution is able to compute a chares placement in 600ms for 2048 chares distributed among 256 cores while the other solutions perform this in less than 50ms. In spite of this duration, the benefits in term of performance counterbalance this additional cost. Finally, it has to be noted that even if the results provided by GreedyLB are really subpar, we observed this behavior on this platform only but for several use-cases and applications.

Blue Waters

Blue Waters is a supercomputer from the University of Illinois at Urbana-Champaign providing a peak performance of 13.34 Petaflops. For these experiments, we used XE6 nodes made of two AMD 6276 Interlagos processors (16 cores each). Each node provides 64GB of main memory. The interconnection network is a Cray Gemini 3D torus as described in 4.3.1. Considering that the Lambb case is able to scale up to 1024 cores, we carried out experiments on 512 and 1024 cores (resp. 16 and 32 nodes) on Blue Waters.

Figure 4a shows the results obtained on 512 cores. For these experiments we carried out several versions of the Lambb case while varying the load imbalance which makes sense scientifically for ChaNGa. On the x-axis, this imbalance is measured by the percentage of active particles (an active particle generates CPU load). As for the experiments on PlaFRIM, we are able at worst to equalize the performances obtained with the best solution. Except compared to GreedyLB, the more imbalance the better performance for our load balancing algorithm. Another remark is that the gains obtained by RefineLB decrease as the the imbalance increases.

Figure 4b presents the execution time of the Lambb case on 1024 cores. In this use case, TreeMatchLB outperforms GreedyLB by 17%. Besides, it also has to be noted that the solution is computed considering eight chares per core, that is to say 8192 chares. This result shows the scalability of our algorithm and the possibility to adapt it to significant cases for coarse-grain paradigms.

Ondes3D

We also did run Ondes3D, a simulator of three-dimen-sional seismical wave propagation. This application has been ported on top of AMPI [START_REF] Huang | Adaptive MPI[END_REF] in order to exploit the Charm++ features [START_REF] Keller Tesser | Improving the performance of seismic wave simulations with dynamic load balancing[END_REF]. As a use case, we worked on a simulation based on the Mw6.6 2007 Niigata Chuetsu-Oki earthquake (Japan) [START_REF] Aochi | Finite difference simulations of seismic wave propagation for the 2007 mw 6.6 niigata-ken chuetsu-oki earthquake: Validity of models and reliable input ground motion in the near-field[END_REF]. The experiments on this real application were carried out on the PlaFRIM cluster.

In Fig 5, we compare TreeMatchLB with the standard Charm++ load balancers GreedyLB and RefineLB using 128 cores. The results show that these load balancers perform very poorly compared to our solution. Though in some cases, an execution without load-balancing step is comparable to a one using TreeMatchLB, we see that our approach is never outperformed and that its gain increases with the number of chares per core. Moreover, compared to previous results we see that the gain has increased with the size of the target platform.

CommBench

The third test application, called CommBench, emulates a large amount of communication between chares. This benchmark is designed as follows: each chare sends a large chunk of data to two chares chosen randomly and the data size is heterogeneous. These communications are particularly irregular. This benchmark was tried out on 8192, 16384 and 32768 cores (resp. (256, 512 and 1024 XE6 nodes) on the Blue Waters platform. Figure 6 shows the results obtained with this benchmark by using our load balancing algorithm TreeMatchLB. The plotted time measures correspond to the mean of ten executions. In order to show the algorithm scalability, we executed the application up to 32768 cores and 32 chares per core. It has to be noted that the standard native Charm++ load balancers (e.g. GreedyLB or RefineLB) did not work at this scale. This is why they are missing on the results. We compared the execution with TreeMatchLB to the one without load balancing. We can notice that the execution time is significantly improved while applying the TreeMatchLB algorithm. On the largest case, we see an improvement of 16.6%. The second remark is that our algorithm is able to scale for an important number of cores and a large amount of computing objects. On 32768 cores, we are able to load balance more than one million chares. This is due to the hierarchy and the distribution of TreeMatchLB. Adapted to a larger-grain paradigm, our algorithm should be able to deal with a substential number of tasks on very large-scale platforms.

Statistical comparison of the different strategies

Based on the above results, it appears that TreeMatchLB is in general the best strategy but that sometimes an other strategy provides similar performance. However, the competitive strategy varies depending on different parameters (Network topology, level of imbalance, problem size, etc.). An interesting question is therefore: is TreeMatchLB statistically better than the other load balancers? To have a statistical significance, of these experiments we measure, in addition of 5 lambb cases the performance of other use-cases: 4 different ChaNGa dwf 5 simulations and 8 Ondes3D runs (these later runs are excluded when we compare Multistep LB or Orb3dLB). All these runs have been chosen such that the amount of communication exchanged during execution 5 dwf demonstrates how disk galaxies can form in a cosmological context. Results are displayed in Figure 7. The figure is read as follows. On the diagonal (in shaded boxes), we have the different load-balancing strategies that we compare. On the upper right, each box is the mean of the difference of the runtime with the strategy on the row and the strategy on the column. For instance, we see that, runs using TreeMatchLB are, on average, 343.51s faster than the runs using RefineLB. On the lower left panel is the 95% confidence interval of the corresponding mean for the strategy on the row and the strategy on the column using the paired Student's ttest (see [START_REF] Jain | The art of computer systems performance analysis[END_REF] chap. 13 p. 209). For instance, for TreeMatchLB vs. RefineLB, the interval is [102s, 585s]. The interpretation of the interval is the following. If the interval is positive there is a 95% probability that strategy on the row outperforms the strategy on the line. If the interval is negative then the strategy on the column outperforms the strategy on the row with probability 95%. Otherwise (the interval contains negative and positive values), we cannot conclude on which strategy is the best with 95% confidence.

Here, we see that TreeMatchLB outperforms all the compared strategies. This backup the intuition that even if in some cases other strategies may perform comparably with TreeMatchLB, on the average TreeMatchLB is the best choice for the type of use-case we are dealing with (i.e. highly communicating simulations). If we compare the other strategies (NoLB, RefineLB, GreedyLB, etc.) between them we see there that RefineLB is better than GreedyLB, Multi- StepLB and Orb3dLB. We also see no statistical difference with NoLB and the 4 other ones. However if we reduce the confidence interval to 80%, the load balancers are almost totally ordered: TreeMatchLB RefineLB MultistepLB Orb3dLB GreedyLB ≈ NoLB.

Sensitivity to initial placement

The initial placement of the processing entities (here the chares) on cores actually impacts the application execution. An example is given on Figure 8 where we ran experiments on the eight cores of one PlaFRIM node. The application is balanced in terms of load: only communications impact the execution time. In this Figure , No LB Round-robin corresponds to an execution without a load-balancing algorithm where the groups of chares are affected on cores by applying the physical identity: group 0 of n chares on the physical core number 0, group 1 of n chares on the physical core number 1, and so on. On the contrary, No LB Packed assigns groups of chares by following the logical identity: group 0 of n chares on logical core number 0, group 1 of n chares on logical core number 1, and so on. On the PlaFRIM platform where the logical numbering differs from the physical one, this initial placement is decisive 6 . Therefore, the user has to know the characteristics of the underlying architecture to assign groups of chares optimally: we have a difference of up to 30% between the two cases. Thanks to TreeMatchLB, we are able to be oblivious to this initial placement. We computed an average of TreeMatchLB runtime starting with the two initial placements described earlier. This result is depicted on Figure 8 with the orange line. This graph shows that regardless of the initial placement, TreeMatchLB is able to move the chares to converge to an optimal placement.

Conclusion

Being able to dynamically balance the load of applications is necessary in the context of highperformance computation due to the hierarchical nature of modern architectures and the cost of communication and the irregular nature of many applications. In this paper, we propose a distributed and hierarchical load-balancer. Our approach is a two-level one. At the network level, we rely on LibTopoMap heuristics for dealing with arbitrary network topologies. For the fine chare-to-core mapping, we use TreeMatch. Nonetheless, we showed that in general, a global load-balancing step is mandatory and we provide such step with a refinement algorithm. At the end, our algorithm optimizes both the communication and the computation costs.

We carried out many experiments on different architectures (PlaFRIM with an Infiniband fattree network and Blue Waters with a Gemeni 3D torus) and various applications and benchmarks. Results show that our topology-aware approach is in general the best strategy and, in the worst case, on par with the best solution. Statistically, our solution is the best on average. Moreover, thanks to its decentralized approach it is able to scale up to one million of chares while standard load balancers do not work for the largest setting we tested. An other interesting feature of our approach is that it relieves the user from the burden of taking into account the core numbering when launching its application.

In future work, we could improve the network awareness. First, some links are not considered. Concerning a 3d torus, two ways are possible in the same axis to link two nodes while we only consider one. Secondly, it could be interesting to evaluate the routing rules for each routers. This would require to manage the way the routing is performed. Another prospect could be to work on load-balancing time and more particularly to let TreeMatchLB chose the two levels of hierarchy to work on. For now, our solution first work at the network level then on the cores level. Should the number of cores per node be very substantial, an interesting feature could be to automatically chose to compute a placement of chares on groups of cores sharing the same cache. This improvement should increase the scalability of this hierarchical approach.

  Cray Gemini router from a 3D torus. One link on the y-axis is reserved to connect the compute nodes.

  Extraction of the network topology of a sub-part of the Gemini 3D torus. Each vertex is a router to which computing nodes are connected.

4. 4 . 1 6 for every node i do in parallel 7 when result[i] has arrived; 8 resultFigure 2 :

 41782 Figure 2: Master-Worker scheme of the Distributed Load Balancer

Algorithm 2 :

 2 Worker Algorithm. Implementation of the compute lb function Input: m chares The communication matrix between chares inside the considered node Input: l chares Array of the load of the Chares inside the current node Input: n Number of chares to be assigned on the considered node Input: c Number of cores of the considered node 1 T ← hwloc get topology(); 2 r ← n/c // Over subscribing ratio 3 p ←tm oversubscribe chares to cores (m chare, T , r); 4 return refine solution(p,l chares);

Figure 3 :

 3 Figure 3: ChaNGa application walltime using different load-balancers for the lambb use-case when varying the number of cores on PlaFRIM.

Figure 4 :

 4 Figure 4: ChaNGa Blue Waters experiments

Figure 5 :

 5 Figure 5: Ondes3D walltime on PlaFRIM when varying the number of chares per cores. 2, 4 and 8 chares/core: 16M cells. 16 chares/core: 64M cells

Figure 6 :

 6 Figure 6: Average execution time of CommBench according to the number of cores with 32 chares per core and 1 MB per Carried out on XE6 nodes of the Blue Waters platform.

Figure 7 :

 7 Figure 7: Statistical comparison of TreeMatchLB with standard load balancers. Upper right panel: mean of the difference of the run duration. Lower left panel 95% confidence interval of the mean using the paired Student's t-test.

2 Figure 8 :

 28 Figure8: Evaluation of the initial placement at launch time for a balanced application whose optimal placement is known.

Table 1 :

 1 p: probability that a given pair of nodes has a load imbalance ratio greater than α when having m nodes, n chares per nodes and the load of the chares is uniformly distributed in [a, b]. q: probability that at least one pair of nodes has an imbalance ratio greater than α.

	n	m	α	a	b	p	q
	64 128 0.05 50	100 14,2% 100%
	32 16	0.1	75	125	0.6% 94.8%
	64 8	0.1	60	140	1.4% 40.4%
	64 256 0.05 100 130 0.02% 99,7%

Please note that NucoLB is designed for shared-memory machines.

Our proposed algorithm and solution have been implemented in Charm++. However, this approach is generic and could find its place in any other runtime system featuring migration, provided the ability to gather the pieces of information described above.

in[START_REF] Jain | The art of computer systems performance analysis[END_REF] Chap. 13, the central limit theorem is applied if n > 30 which is the case for all the runs we are doing

If X ∼ N (µ, σ 2 ) then X-µ σ ∼ N (0, 1)

Indeed, the cores of the Intel Xeon Nehalem included in a PlaFRIM node present this physical numbering scheme: 0,2,4,6,1,3,5,7

Acknowledgment

This research is partially supported by the NCSA-Inria-ANL-BSC-JSC-Riken Joint-Laboratory on Extreme Scale Computing (JLESC) in which we had allocated hours to the Blue Waters machine. We would like to thank Torsten Hoefler for providing us with the LibTopoMap code. The Pla-FRIM experimental testbed, is being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr). Part of this work has also been funded by the ANR MOEBUS project (ANR-13-INFR-0001). We would like to thanks Fabrice Dupros et Rafael Keller Tesser for providing us with the Ondes3D code.


Zoi Christoforou
email: zoi.christoforou@enpc.fr

Etienne Corbillé
email: etienne.corbille@eleves.enpc.fr

Nadir Farhi
email: nadir.farhi@ifsttar.fr

Professor, Fabien Leurent
email: fabien.leurent@enpc.fr

MANAGING PLANNED DISRUPTIONS OF MASS TRANSIT SYSTEMS

Unplanned disruptions of rail transit networks have been extensively studied. Planned disruptions for works are essentially different mainly due to their longer duration that allows passengers to build alternative route choice strategies. The literature on this topic remains scarce. We propose a novel methodology that enables operators to evaluate different disruption management schemes and to obtain explicit estimations of travel times, passenger comfort flows, and levels of service. Statistical tools are used for the evaluation of the different strategies. We illustrate the methodology through a large-scale application to a real line disruption in Paris, France. The disruption took place in July 2015 and was due to network maintenance works. It affected the major suburban railway line RER-A which counts over 1,000,000 trips on a typical working day. Results indicate that the disruption would have significantly increased the generalized cost (GC) of passengers if no action was taken. The operator's disruption management scheme included bus bridging and frequency increase on alternative routes. The evaluation of this plan shows that it restored the average GC over the whole network. Passengers initially using the disrupted line experienced increased GC by using alternative longer routes. Passengers initially using those alternative routes experienced lower GC thanks to the frequency increase. Finally, capacity problems are observed on the buses assuring the bridge of the disrupted link.

INTRODUCTION

On one hand, urban population is rapidly growing and travel demand is increasing in city centers worldwide. On the other hand, urban space is scarce and new road infrastructure becomes costly. As a result, urban agglomerations suffer from heavy congestion, long delays, and serious travel time variability. Commuters turn to public transit systems that also begin to suffer from capacity phenomena. Heavy reliance on mass transit systems and operations at capacity levels make transit systems very vulnerable to disturbances. Simple service disruptions lead to undesirable outcomes such as delay propagation and passenger inconvenience. Reliability and robustness of transit systems are major issues that receive the attention of both researchers and practitioners.

The link between transportation and city planning is known and well documented. Regarding public transport, planning includes a range of problems varying from strategic to operational level: network design, line planning, and timetabling [START_REF] Jin | Disruption Response Planning for an Urban Mass Rapid Transit Network[END_REF]. Traffic assignment models have long be used for planning purposes. These models provide an estimation of individual route choices based on a set of network and demand parameters such as origin-destination (OD) pairs, mode availability and capacity. We should note, though, that when used for planning purposes, assignment models are mainly static and consider only normal operating conditions. Regular maintenance, works and other disruptions are not integrated. In many cases, infrastructure capacity is not taken into consideration either.

However, disruptions occur regularly. Unplanned (or unforeseen) disruptions may be caused by natural disasters, terrorist attacks, accidents, medical emergencies, adverse weather, infrastructure or rolling stock failures [START_REF] Zhu | Disruptions to Transportation Networks: A Review[END_REF]. Planned (or foreseen) disruptions occur mainly due to strikes, special mega-events, maintenance and upgrade works. Disruptions differ significantly in respect to their magnitude and scale, their duration and impact. The impact concerns both demand and supply. In fact, behavioral changes (i.e. daily activities, departure time, modal and route choice) are observed during disruptions that may even alter mode choices in the long-term [START_REF] Zhu | Disruptions to Transportation Networks: A Review[END_REF]. On the supply side, operators have to assure the static and dynamic resilience of the network by developing and implementing disruption management schemes and programs (DRPs). Dynamic resilience refers to the capacity and time to recover after a shock and is mainly relevant to unplanned disruptions. Static resilience describes the ability of the system to operate under degraded conditions and is most relevant to planned events.

A key element to sustainable operations is robustness that assures low sensitivity to disturbances. Indexes for the network's robustness evaluation for the case of random failures ('mean indexes') and intentional attacks ('worst-case indexes') can be found in the literature [START_REF] De-Los-Santos | Evaluating passenger robustness in a rail transit network[END_REF]. Resilience should also be included in economic models. In (4), the authors add at the cost function a term representing operator and passenger costs due to unplanned delays. The authors value the GC of nonscheduled delay more highly than either in-vehicle travel time or scheduled waiting time. Following the earlier indication in [START_REF] Small | The scheduling of consumer activities: work trips[END_REF], they take a coefficient of 2.4 times the in-vehicle travel time.

Nevertheless, most agencies apply DRPs that do not account for user costs and do not assure system optimization. Unplanned short-term events are handled real-time with operators trying to minimize the disruption impact based on their experience and ad-hoc decision support tools. On the contrary, DRPs can be prepared and tested before a planned disruption occurs. DRPs are also relevant in the case of unplanned disruptions of long duration. As an example, the I-35W bridge over the Mississippi River collapsed in 2007 and interrupted the usual route of 140,000 daily vehicles for 66 days. The daily cost to users was estimated between $71,000 and $220,000 [START_REF] Zhu | Disruptions to Transportation Networks: A Review[END_REF]. Despite the numerous studies on unplanned disruption management, the literature on planned disruptions of mass transit systems remains scarce. Transit assignment tools have rarely been used for this purpose.

This paper proposes a novel methodology for the evaluation of DRPs using a traffic assignment model that accounts for capacity constraints and passenger comfort. This methodology can be implemented at all stages of the planning procedure to compare the robustness of different alternatives and the system's sensitivity to disturbances. We , first, adapt the OD matrix to the disruption characteristics. Then, we perform assignment under capacity constraints and estimate not only delay (due to frequency reduction and/or link removals)but also passenger comfort states (invehicle congestion, on-platform waiting times). Finally, we statistically compare results and quantitatively evaluate different DRPs. The proposed methodology is tested on a real transit network, the Ile-de-France network, which accommodates 8.3 million trips on a typical day [START_REF]Public transport organizingauthority of Paris region (STIF), Les nouvellesmobilités des Franciliens[END_REF]. A major disruption on the suburban railway was planned for summer 2015 due to maintenance works (7). We apply the methodology to evaluate the DRP proposed by the operator.

The remainder of this paper is organized as follows. Section 2 reviews the literature, section 3 explains the proposed methodology, section 4 presents the application to the Ile-de-France network, and section 5 discusses the results and draws conclusions.

BACKGROUND

Many researchers work on road traffic operations during short and long disruptions [START_REF] Fan | Modeling urban hazmat transportation with road closure consideration[END_REF][START_REF] Pillac | An event-driven optimization framework for dynamic vehicle routing[END_REF]. When a road link fails, it is usually quite simple to find alternative paths and reroute traffic. Severe delays may occur but the connection is assured. The objective is, then, to optimize the assignment procedure by minimizing delays. In the case of planned works, the problem becomes even simpler as alternative schemes can be tested and evaluated off-line through simulations as in [START_REF] Kattan | Travel behavior changes and responses to advanced traveler information in prolonged and large-scale network disruptions: A case study of west LRT line construction in the city of Calgary[END_REF].

In the case of rail transit networks, additional difficulties exist due to limited interconnections and overtaking possibilities. Unplanned disturbances may vary from minor incidents and slow-moving delays to major incidents that close entire track sections [START_REF] Schmocker | Metro Service Delay Recovery: Comparison of Strategies and Constraints across Systems[END_REF]. In real-time contexts, operators face two main challenges: (i) the incident duration is not known in advance and(ii) the status of the network changes quickly [START_REF] Cadarso | Recovery of disruptions in Rapid Transit Networks with Origin-Destination Demand[END_REF]. In [START_REF] Schranil | Forecasting the Duration of Rail Operation Disturbances[END_REF], the authors propose a methodology to forecast the duration of disturbances and, thus, reduce their impact on the scheduled timetable. This framework can be used by operators as a decision aid tool. In [START_REF] Corman | Disruption Handling in Large Railway Networks[END_REF], the authors propose an advanced decision support system to help dispatchers make real-time changes and to limit delay propagation. The network is divided into dispatching areas of limited size and various disruption scenarios are tested off-line. Performance is evaluated in terms of travel time, waiting times, train delay.

The optimal strategy also depends on the nature of the unplanned disruption. If a train is disabled, operators re-route commuters to the next train or to another mode and balance out the frequency of the remaining train service [START_REF] Pender | Planning for the unplanned: an international Review of Current Approaches to Service Disruption Management of Railways[END_REF]. In the case of complete line blockages, the re-routing also concerns trains while the frequencies of substitute lines increase. Buses or even taxis are used to assure the disrupted links. In [START_REF] Kepaptsoglou | Bus bridging Problem in Metro Operations: Mathematical Framework and Decision Support System[END_REF], an overview of ad-hoc bus bridging services is provided and an efficient decision support system is proposed. If reserve buses are strategically located over the network, the response to disturbances can be further ameliorated (i.e. shortened). The optimal location can be obtained as a function of disruption type likelihood, passenger volumes, travel time to bus bridging termini, and so on [START_REF] Pender | Improving response to Unplanned Passenger Rail Disruption: a New Method to Locate Rail-Replacement Bus Reserves[END_REF]. Another proactive recovery strategy is to optimally locate track crossovers as in [START_REF] Pender | Proactive recovery from rail Disruptions Through Provision of Track Crossovers and Bus Bridging[END_REF] or to predefine measures in case of certain disruptions that operators can apply directly and achieve time savings [START_REF] Chu | Modeling Capacity Consumption Considering Disruption Program Characteristics and the Transition Phase to Steady Operations During Disruptions[END_REF]. These techniques provide solutions for short-term disruptions or low-demand networks or require space that is not always available in dense urban areas.

The unplanned disruption management problem can be divided in three sub-problems: (i) Timetable adjustment. It is commonly achieved by control strategies such as holding and short-turning [START_REF] Krasemann | Design of an effective algorithm for fast response to the re-scheduling of railway traffic during disturbances[END_REF][START_REF] Lorek | Disruption Management with Rescheduling of Trips and Vehicle Circulations[END_REF][START_REF] Sato | A MIP-based timetable rescheduling formulation and algorithm minimizing further inconvenience to passengers[END_REF]. The objective is either schedule-adherence for long headways or headwayadherence for high frequencies. Many researchers have proposed optimization-based algorithms for rescheduling including various elements in the objective function. In [START_REF] Krasemann | Design of an effective algorithm for fast response to the re-scheduling of railway traffic during disturbances[END_REF], the algorithm performs a depth-first search using an evaluation function to prioritize when conflicts arise and then branches according to predefined criteria achieving significant time savings. In [START_REF] Lorek | Disruption Management with Rescheduling of Trips and Vehicle Circulations[END_REF], the objective function maximizes the number of recovered trips and guarantees a conflict-free timetable during the disruption. Delay is not constrained though. In [START_REF] Sato | A MIP-based timetable rescheduling formulation and algorithm minimizing further inconvenience to passengers[END_REF], train delay is not constrained either but the objective function minimizes inconvenience to passengers that is defined as the sum of traveling time on-board, on-platform waiting time and number of transfers. This formulation may result to prolonged dwell times in order for trains to wait for passenger arrivals. Train capacity constraints are considered for both interstation segments and station tracks.

(ii) Evacuation of the 'blocked' passengers. Interactions between supply and demand are complex and the optimal strategy may be affected by demand conditions [START_REF] D'acierno | Evaluation of Travel Demand Impacts in the Case of Rail System Failure[END_REF]. An international survey across 71 transit agencies suggests that only 20% of the agencies have parallel transit systems that can evacuate passengers [START_REF] Pender | Disruption Recovery in Passenger Railways: International Survey[END_REF]. Bus bridging seems to be the most common response but rarely do agencies keep a strategic vehicle reserve for this purpose. Thus, some researchers have expanded rescheduling algorithms to include passenger evacuation. In [START_REF] Kumakawa | A Novel Train Rescheduling Algorithm for Correcting Disrupted Train Operations in a Dense Urban Environment[END_REF], the authors propose a rescheduling algorithm that has two interesting features: it accounts for the passengers' waiting while queuing on platforms and the time needed to implement a plan is quantified in the delay estimation. However, the user perspective and the in-vehicle comfort state are not taken into account. In [START_REF] Cadarso | Recovery of disruptions in Rapid Transit Networks with Origin-Destination Demand[END_REF], the passenger decision making is integrated and capacity is considered. Passengers make decisions about itineraries according to a logit model and then try to board a train serving that itinerary. This relaxation allows passengers to change their initial decisions and adapt to the changing environment. However, the information leading to mind-changes remains limited. On the contrary, traveler perception is shaped by their individual prior-knowledge, accumulated experience as well as information availability [START_REF] Cats | Multi-agent Transit Operations and Assignment Model[END_REF].

(iii) Rolling stock and crew rescheduling. Many optimization algorithms under availability constraints can be found in the literature (see for example [START_REF] Cadarso | Recovery of disruptions in Rapid Transit Networks with Origin-Destination Demand[END_REF][START_REF] Corman | Disruption Handling in Large Railway Networks[END_REF] but are beyond the scope of this paper.

Planned disruptions differ from unplanned disruptions in many ways. A first 'chaotic' phase is observed immediately after an unplanned disruption. During this phase, operators and passengers have limited or no visibility of the cause, magnitude, and duration of the failure [START_REF] Chu | Modeling Capacity Consumption Considering Disruption Program Characteristics and the Transition Phase to Steady Operations During Disruptions[END_REF]. This phase does not exist in the case of planned works especially when communication campaigns are in place. When unexpected incidents occur, an important segment of the demand is obliged to make en route rerouting and real-time decisions with limited information at-hand; whereas, during planned disturbances, the route is decided in advance. Information availability changes for the operators as well. Furthermore, the time horizon is entirely different and, so, system equilibrium may be possible for planned works but not for short disruptions. Also, planned disturbances allow operators to predict and assure the availability of necessary resources; both human and technical. Concluding, best strategies and DRPs differ significantly from one case to the other.

Despite the significant research on unplanned disruptions, very few researchers have worked on the optimization and evaluation of DRPs for short planned disruptions and even fewer on longer planned disruptions. In [START_REF] D'acierno | Evaluation of Travel Demand Impacts in the Case of Rail System Failure[END_REF], the authors analyze effects on travel demand of different levels of degraded rail services during short disruptions. Their contribution is important as the authors consider both supply and demand. In particular, they introduce capacity constraints of rail vehicles and a FIFO approach for simulating the sequence of users boarding trains. They argue that mingling is not relevant in the case of major incidents on the network due to high crowding. In [START_REF] Sanchez-Martinez | Event-Driven Holding Control for High-Frequency Transit[END_REF], the authors propose a holding control model based on rolling-horizon optimization capturing dynamics of running times and demand that can be applied to high-frequency transit disruptions. They consider both foreseen and unforeseen events under two control scenarios: naïve and informed control. They test these scenarios in a toy instance of a bus route having 20 stops per direction and find important performance improvements achieved by the informed control model. This model uses preemptive holding and other assumptions that can handle short disruptions but are unsuitable in the case of long disruptions of several weeks. In [START_REF] Forsgren | Maintaining tracks and traffic flow at the same time[END_REF] we find a study for longer disruptions. The authors consider planned works in a mixed integer programming model that optimizes a production plan with regard to both trains and preventive maintenance. Thus, they estimate best traffic flow schemes given a fixed set of planned maintenance activities. However, passengers' point of view is not considered.

In conclusion, there is a complete absence of studies dedicated to disruption response planning and to the application of optimization techniques despite their potential utility to operators [START_REF] Jin | Disruption Response Planning for an Urban Mass Rapid Transit Network[END_REF].Furthermore, studies on unplanned disruptions are macroscopic and do not always account for the infrastructure capacity and the resulting passenger comfort state. If capacity is not considered, the simulated failure effects tend to be calmed more rapidly than in real cases [START_REF] D'acierno | Evaluation of Travel Demand Impacts in the Case of Rail System Failure[END_REF]. Our contribution with respect to the existing literature may be summarized as follows:

-We study planned rail transit disruptions; -We perform re-assignment considering both the supply and demand structure; -We account for network capacity and passenger comfort; -We propose a statistical evaluation of different strategies; -We illustrate our findings for the case of a real disruption of a large-scale network.

METHODOLOGY

The proposed methodology consists of the following steps: 1) Supply 'deformation': adjusting the network to the planned disruption by removing links, adapting frequencies, adding new links for bus bridging, and so on; 2) Demand 'deformation': estimation of public transport modal share during disruption period; 3) Scenario building: alternative DRPs to be evaluated; 4) Assignment for each scenario using the model CapTA; 5) Statistical comparison of results.

In [START_REF] Voith | Consumer response to Service Interruption on the Washington Metro[END_REF] the authors present the case of a disruption at the Washington metro and study its impact upon demand. They conclude that a significant shift occurs towards other modes. Reassignment is thus necessary in order to estimate delays more accurately. We give below an overview of CapTA (Capacitated Transit Assignment), the transit assignment model used in this study for the reassignment. The reader can refer to [START_REF] Leurent | Poulhès traffic assignment model for passenger transit on a capacitated network: Bi-layer framework, line sub-models and large-scale application[END_REF] for a detailed description of CapTA, together with an application to the greater Paris area on a typical working day of 2010, including a special focus on the RER-A line: at the morning peak, the westbound passenger flow in central Paris amounts to more than 50,000 passengers per hour i.e. the line is equivalent to 25 motorway lanes at capacity load. It is thus even more important to model the full set of alternative routes and their own sensitivity to congestion.

The CapTA model has a three-layered architecture: (i) The upper layer represents passenger route choice on a network of pedestrian links and transit legsa leg being a pair of stations for access and egress along a given transit line. Every passenger is assigned to its least-cost hyperpath, in the framework of static, frequency-based, optimal strategy assignment.

(ii) The lower layer deals with vehicle runs along the service routes. Each such route is addressed on an individual basis, with an average vehicle that has given seat capacity, total in-vehicle capacity, and exchange capacity. The flows of boarding and alighting passengers are assigned to the vehicle along the stations that it serves, yielding in-vehicle passenger loads and comfort (with sitting priority according to chronology of passenger access to the vehicle). Moreover, the station dwell time depends on alighting and boarding flows in relation to exchange capacityhence the run times by trip leg.

(iii) The middle layer addresses specific line issues. First comes the coordination of the service routes that make up the line, especially so for the RER-A line that has two eastern branches and three western branches grafted in a central trunk: this yields combined frequency and residual capacity by line leg. Then, the residual capacity at a given station is faced to the flow of willing-toboard people, assumedly mingled, in a special bottleneck model that yields average passenger wait time and vehicle boarding probability by leg, thereby determining the passenger load of each vehicle serving the station. Next, still by station, the occupancy of the vehicle track due to the dwell times associated to the routes that serve the station, together with the clearance times of the routes that do not stop there and the safety margins between the vehicles, is evaluated for the assignment period on the basis of the target frequency assigned to the routes. Whenever the demanded occupancy exceeds the period duration, a reduction ratio is applied to the target frequency, yielding a reduced frequency which is propagated downstream. Lastly, frequency reduction entails a related delay to subsequent trains: this is evaluated on the basis of a bottleneck formula and propagated upstream in a specific way.

To conclude, the CapTA model captures a wide range of traffic phenomena. These are integrated in a consistent way, owing to the leg representation. The state vector involves (i) passenger flows by link and destination, (ii) service frequency by route and node along it and (iii) vectors of strategy assignment by destination and current node. An equilibrium state must exist in theory. In practice, about one hundred iterations of a Method of Successive Averages enable one to get a traffic state with satisfactory convergence to equilibrium.

In the following section, we apply step-by-step the proposed methodology to a real disruption that took place in summer 2015 (7).

APPLICATION TO THE RER-A PLANNED WORKS Supply 'deformation'

Greater Paris counts 41 million trips on a typical working day and over 12 million inhabitants that spend on average 1.3 hours for transportation every day. The public transport network is extensive and absorbs 20% of the total trips made. It consists of 14 metro lines, 4 suburban railway lines, 1,500 bus lines, 8 tramway lines, and 8 regional train lines. The most congested line is the suburban railway line RER-A that counts over 1,000,000 passengers on a typical working day. It crosses the region in an east-west direction and serves the major business center of La Défense.

RER-A infrastructure is ageing and serious maintenance works must be undertaken. A program of works was decided as shown in Figure 1. More information can be retrieved on the operator's website (7).In 2015, the RER-A closed between July 25 th and August 23 rd . At the same time, another less important suburban line, the RER-C, partially closed between the stations Javel and Gare d'Austerlitz. The operator used replacement buses for the two closures and increased frequencies on substitute modes (i.e. mainly metro line 1) at winter frequency levels and beyond. This was not the case for other lines that operated at reduced summer frequencies as usual. The annual capacity reduction when passing from winter to summer schedules is generally estimated at 30%. All of these supply changes have been integrated in the network building of CapTA.

FIGURE 1Planned works on the RER-A.

Demand 'deformation'

The timing of works during French summer school vacation imposes the adjustment of the OD 'winter' matrix downwards. We estimated summer demand using the recent household travel survey (EGT 2010) that is publicly available and a series of hypotheses: 'studies' were removed from travel motives, 'work' was reduced by 35%, and tourist levels were upgraded. In total, we found a 40% reduction when passing from winter to summer demand. Figure 2 illustrates the origins of all trips to 'La Défense' station after the adjustment of demand. Figure 2 reveals that most incoming flows start from the inner city center and will be, thus, affected by the disruption.

FIGURE 2 Passenger flows to 'La Défense' per origin during disruption

Scenario Building

The following scenarios were simulated once without capacity constraints (sub-scenarios-a) and once with capacity constraints (sub-scenarios-b):

-S0: Winter demand, winter supply, without disruption; -S1: Summer demand, summer supply, without disruption; -S2: Summer demand, summer supply, with 2015 disruption; -S3: Summer demand, with 2015 disruption with winter supply and bus bridging. S0 represents the business-as-usual operating plan for wintertime and S1 represents the business-as-usual operating plan for summertime when both demand and supply are decreased. The supply is decreased mainly by frequency reduction. S2 refers to the 2015 summer disruption if no action is taken. S3 represents the official operational plan that has been prepared by the operator (i.e. the RATP). This plan includes the maintenance of the winter frequency for Metro Line 1 and the operation of a bus bridging line linking Auber to La Défense.

Assignment and Statistical Analysis

Simulation results provide passenger volume (v)and volume-to-capacity (v/c) ratio per link as well as the GG per OD-pair. In CapTA, the GC is calculated as the sum of travel time (including transfer), waiting time, and comfort state. Waiting time estimates are a function of frequencies and vehicle capacities. Two in-vehicle comfort states are considered: seated vs. standing. In total, 3,048,317 ODpairs were analyzed with 62,728 network arcs. Figure3 shows the frequency distributions of the GC and the flows for S0a. Table1 presents some descriptive statistics for the datasets obtained. We first observe that the sample is not normally distributed but very skewed. Let us now compare between the a-and b-variants of each scenario. Table 1 reveals that if link capacity is not considered (a-variants) the v/c ratio largely exceeds the 100% which is unrealistic. The latter is mainly the case when the disruption takes place (under S2 and S3) where we obtain a v/c ratio of 721%. Consequently, ignoring capacity when performing transit assignment leads to erroneous results, especially when we have line closures. This is also illustrated in Figure4 where we plot the GC of the a-variants over the GC of the b-variants. As a result, in the remainder of this paper we only consider the b-variants of the four scenarios. Focusing to the analysis of the capacitated assignments, we obtain Figure5 that gives the frequency distributions of the GC. Under normal operating conditions (S0,S1), the GC is higher during wintertime versus summertime but the two distributions have a similar form. This can be intuitively explained by the lower level of service and the congestion phenomena due to increased demand. If the line closure happens and no action is taken (S2), the GC takes its maximal value. It is not only higher than the summer GC (S1) but even higher than the winter GC (S0). It also shows a significantly greater variance as many OD pairs exist that have a GC of around 10. We can reasonably assume that these OD-pairs correspond to the affected routes. This variance disappears under S3 when the DRP is implemented and the average GC becomes even lower than the one of normal operating conditions. The mean percentage GC error for S2b versus S3b is of 10.23%, for S2b versus S1b of 9.52%, and of S1b versus S3b of only 0.49%. The difference between S2 and S3 can be easily visualized in Figure6 that was produced with the TransCAD software. In Table 2, pairs of scenarios are compared in term of GC using two statistical indicators: MPE (mean relative prediction error) and RMSE (root mean squared relative prediction error).

FIGURE 5Frequency distributions of GC for

Focus on selected itineraries

In this section, we focus on the OD-pairs that are expected to be mostly affected by the disruption and we explore the efficiency of the planned management scheme (S3). Table 2 gives the GC for selected origins and destinations: La Défense (Paris major business district), Nanterre Université (important attraction pole), Les Halles (the major transit hub of Paris), and the 8 th arrondissement (where the disruption is located). Comparing S2 to S3, we observe that, indeed, the planned scheme seems to alleviate the problem for all itineraries. Nevertheless, it becomes obvious that the passengers using the links that are directly affected by the disruption will experience increased GC compared to the normal operating conditions under summer (S1) or even wintertime (S0). They will opt for alternative longer routes as a second best choice given that public transport demand has a low elasticity and that trave l time differences (between RER and metro) are not very high. However, the capacity of the links serving alternate routes will not be exceeded as previously. The fact that the average GC over the network is equal to the GC under normal operating conditions (see previous section) implies that the GC of some other passengers will decrease. These passengers are the ones who were initially using the alternative routes whose frequency will increase according to the management disruption scheme. Their travel time decreases (increased frequency) without significant loss of level of service as capacity is not exceeded.

CONCLUSIONS

Planned disruptions of railway transit systems occur frequently and are costly to both users and operators. Despite the numerous studies on unplanned railway disruptions, few are the studies investigating planned railway disruptions. In this paper, we propose a novel methodology that enables operators to evaluate different disruption management schemes and to obtain explicit estimations of travel times, passenger comfort flows, and levels of service. We apply this methodology to a real large-scale disruption. Results indicate that the operator's plan will indeed restore the network's performance on average. However, when looking closer, we find that some passengers will experience higher GC and some others lower.

The added value of this paper is manifold. We provide quantitative evidence of the necessity of performing re-assignment when elaborating a DRP. Otherwise, important network effects may be neglected. Also, we prove that, if capacity constraints are ignored, results are entirely unrealistic. Furthermore, we show that all macroscopic analyses at the network should be combined with microscopic analyses at the level of the affected routes. Moreover, we show the inefficiency of bus bridging strategies in the case of large-scale and long disruptions when passengers can plan in advance their routes. Besides, bus capacities remain marginal compared to the metro and cannot accommodate demand. Bridging buses are interesting options for passengers that are not familiar to the network (such as tourists) or in cases where the landscape does not allow for any other alternative (such as for a major water crossing).

Finally, our research suffers from some shortcomings. Modal choice was not considered even if shifts towards private cars could occur. The role of information campaigns was not discussed and its impact upon the assignment was not investigated. At this point of time, we ignore the real-time guidance that the operator will choose to provide (itinerary calculators and so on). On the modeling side, we cannot be sure that the assignment is exact as the system's equilibrium will not necessarily be attained during the first days of the disruption. We intend to explore these issues and, also, compare our simulation results to what actually happened during the disruption using smart card validation data and other measurements.

FIGURE 3

 3 FIGURE 3 Frequency distribution for (a) passenger flows per arc and (b) GC per OD.

FIGURE 4

 4 FIGURE 4 Comparison of the GC of the a-over the b-variants for S0, S1, S2, S3.

  FIGURE 5Frequency distributions of GC for S0b, S1b, S2b, S3b.

TABLE 1Descriptive Statistics for Simulation Results per Scenario.

 1Descriptive

	Scenario	V/C per link	GC per OD
		average	SD	average	SD
	S0a	153,7%	8,4%	2,30	1,15
	S0b	108,1%	8,0%	2,24	1,43
	S1a	645,2%	33,7%	2,11	1,04
	S1b	132,5%	8,1%	2,21	1,43
	S2a	721,2%	39,3%	2,08	1,02
	S2b	134,9%	9,1%	2,48	1,75
	S3a	115,0%	6,7%	2,08	1,02
	S3b	94,0%	6,3%	2,17	1,41

TABLE 2 Generalized cost (GC) per scenario for selected itineraries.

 2

		Destination:	Origin:	Destination:	Destination:	Destination:	Origin:
		La Défense	La Défense	8 th	Nanterre Univ	Les Halles	Les
				arrondissement			Halles
	S0b	1.24	0.80	1.17	1.67	1.10	0.69
	S1b	1.20	0.78	1.14	1.63	1.08	0.67
	S2b	1.44	0.90	1.24	1.85	1.17	0.76
	S3b	1.29	0.83	1.13	1.69	1.07	0.71

ACKNOWLEDGMENTS

The authors would like to thank the Public Transport Organizing Authority of the Paris Region (STIF) for partially funding this research. We also thank Shaoqing Wang for the preparation of maps.