Guo-Liang Chai

Mauro Boero

Zhufeng Hou

Kiyoyuki Terakura

Wendan Cheng

Indirect Four-Electron Oxygen Reduction Reaction on Carbon Materials Catalysts in Acidic Solutions

Keywords: Oxygen reduction reaction, Fuel cells, Carbon materials catalysts, Acidic Solutions, First principles simulations, Metadynamics

Developing cost-effective and high performance oxygen reduction reaction (ORR) catalysts is a fundamental issue in fuel cells and metal-air batteries. To this aim, carbon materials catalysts (CMCs) are extensively investigated because of their performance comparable to noble metal based catalysts in alkaline solution. Yet, acidic solutions are desirable for an efficient proton-exchange across Nafion membranes to yield high power density for commercial applications. However, the ORR performance of CMCs in acidic solutions is rather low because of undesirable two-electron processes and OH radical formation. By using first principles simulations, in this work we elucidate the mechanisms and identify the active sites of 2e -ORR processes for indirect 4e -ORR. We provide evidence for the fact that nitrogen doped Stone-Wales defects in graphene favor an indirect four-electron ORR upon H2O2 formation and reduction. The low ORR potential for metal free CMCs is ascribed to H2O2 formation via hydrogen abstraction and the critical point for OH radical generation on transition metal based CMCs is 0.82 V. Moreover, we provide an insight into the indirect 4e -ORR which serves as a guide for suppressing undesired 2e -ORR, avoiding OH production, and promoting direct four electron ORR on CMCs. These results disclose a new strategy for developing high efficient ORR on CMCs in acidic solutions.

INTRODUCTION

Electrochemical oxygen reduction reaction (ORR) is the common process occurring at the cathode of fuel cells and metal-air batteries. Such a reaction is kinetically rather slow and, as such, is responsible for the degree of efficiency of electrochemical devices for clean energy applications. [1][2][3][4][5][6] To date, an economically affordable cathode ORR catalysts represents still a stumbling block because of the high cost of the routinely used noble metals. 7 More recently, carbon materials catalysts (CMCs) have been pioneered as an economically appealing alternative for the expensive Pt based catalysts. [8][9] In these systems, the specific framework of the CMCs may be graphene, carbon nanotubes, fullerenes, or amorphous carbons, depending on the preparation, synthesis procedure and experimental conditions. Moreover, an appreciable enhancement of the catalytic activity of these carbon-based structures is obtained by doping with B, N, P, S and other heteroatoms. 10 To reach the stage of practical application, it is crucial to develop ORR catalysts able to operate in acidic solution, since proton transfer processes across Nafion membranes are much faster in acidic electrolytes than in alkaline ones and a fast proton transfer is crucial to generate a higher power density. Although the ORR performance of some metal free CMCs have already been shown to be comparable to the widely used Pt catalysts in alkaline solution, their performance in acidic solution is still far below the one of Pt because of the occurrence of 2e -ORR to from H2O2. 8,[10][11][12][13] Transition metal (TM) based CMCs, such as the worldwide studied FeNx systems, are favorable to ORR in acidic solution. Their drawback is the fact that these systems realize Fenton's reactions, reducing H2O2 to OH radicals and destroying the catalysts. 14 The issues briefly summarized here for ORR in acidic solutions are mainly related to an indirect four electron ORR process that can be seen as composed of distinct two electron ORR processes. 15 However, to date, no study clarified the mechanisms for this indirect four electron ORR process operated by CMCs in acidic solution. Yet, this is a serious issue that hinders the commercial application of CMCs to replace Pt catalysts in fuel cells. Hence, understanding the underlying mechanisms of such an indirect four-electron ORR in acidic solution is of primary importance for the development of efficient CMCs.

In most of the experimental studies reported so far, the discrimination between two-electron and four-electron ORR is done on the basis of the electron transfer number (ETN) deduced from the percentage of produced H2O2. However, this empirical method is ambiguous and unable to distinguish direct and indirect four-electron ORR processes. For instance, if hydrogen peroxide H2O2 is first formed via a two-electron ORR and then the H2O2 is further reduced by a subsequent two-electron ORR to water, the ETN of this indirect ORR pathway turns out to be equal to four, thus identical to a full four-electron process. A recent investigation of the ORR processes on CMCs rationalized in terms of Damjanovic and Wroblowa models claimed that both two-electron process, direct and indirect four-electron pathways, can occur in CMCs, while the corresponding mechanisms and the active sites are still a matter of debate. 15 Here we make use of first principles simulations to unravel the indirect four-electron ORR process in CMCs consisting in (i) a two-electron H2O2 formation and (ii) a two-electron H2O2 reduction in acidic solution. The catalytic models considered are nitrogen doped metal-free and Fe based CMCs. Our results indicate that H2O2 can form via a hydrogen abstraction mechanism on CMCs surfaces. We also find that at least three ORR mechanisms are possible for H2O2 reduction via the generation of a hydroxyl anion OH -, an OH radical or an H2O molecule, respectively. The N doped Stone-Wales defect can promote the formation of H2O2 via hydrogen abstraction and favor the H2O2 reduction via the mechanisms mentioned above. This facilitates a combined indirect four-electron ORR process. Our findings suggest also that the hydrogen abstraction mechanism can lower the electro-chemical potential to ~0.4 V and hence lower the efficiency of ORR. 16 Instead, the OH radical generation during the H2O2 reduction would decrease the durability of the catalyst with a critical half-wave potential of 0.82 V. Thus, the development of high performance CMCs in acidic environments requires methods and techniques suited to overcome these negative issues.

SIMULATION METHODS

Nitrogen doped non-defective graphene and an analogous system including a Stone-Wales defect were considered, on the basis of all our former works, as model templates for CMCs (see Figure 1). Simulations were performed with the Car-Parrinello molecular dynamics (CPMD) method within the density functional theory (DFT) framework. [17][START_REF] Cpmd | MPI fur Festkorperforshung[END_REF][START_REF] Kohn | [END_REF] This approach is sufficient to investigate processes occurring on a picosecond time scale for which activation barriers of a few kBT are involved. This is for instance the case of the activation of an O2 molecule by hydrogen abstraction. The simulation cell of each system is orthorhombic and contains about 200 atoms as specified elsewhere. 9 Lattice parameters for each structure considered are reported in the supporting information. As an example, the cell lattice parameters for the G-N structure are a=12.28, b=12.76, c=15.00 and α=β=γ=90°. CPMD simulations were complemented by the Blue Moon ensemble (BME) free energy sampling technique in a canonical NVT ensemble at 300 K whenever the free energy barriers involved are too large to be overcome by standard CPMD and a single reaction coordinate is sufficient to describe the activated process. The CPMD equations of motion are numerically integrated with a time step of 4.0 a.u.. 20 to ensure a good control of the constants of motion. The exchange-correlation functional was described by a gradient-corrected spin-unrestricted Hamprecht-Cohen-Tozer-Handy (GGA-HCTH) formulation. 21 The corevalence interaction was described by norm-conserving Troullier-Martins pseudopotentials (PP) for C, N, O and P, 22 while an analytical von Barth-Car one was used for H. 23 The sampling of the Brillouin zone was restricted to the Γ point. The spin state of the various systems considered here are identical to the ones treated in our former work. 9 Whenever a single reaction coordinate is insufficient, as (e.g.) for the H2O2 dissociation on CMCs, we resorted to the metadynamics approach [START_REF] Laio | Proc. Natl. Acad. Sci[END_REF][START_REF] Barducci | [END_REF] to inspect reaction pathways and to sort out free energy barriers. More precisely, we adopted the variational formulation of the metadynamics 26 in which collective variables (CVs) sbecome dynamical degrees of freedom in the Lagrangean:

2 . 2 1 1 [R,{ }] [ ( ) ] ( , ), (1) 2 2 CP i L L M k R V t s s s s                being L CP [R,{ᴪi}]
is the standard Car-Parrinello Lagrangian. 17 The second and third terms in the right hand side of Eq. ( 1) are the fictitious kinetic energy and the restraining potential, forcing the instantaneous CV s (t) to fluctuate around its corresponding mathematical expression s (R). The time-evolving term V(t,s) is the history-dependent potential where penalty functions accumulate along the metadynamics trajectory. This allows to fill progressively the free-energy minima and to force the reactions to occur on an affordable time scale. The explicit form of the potential used in this work is a superposition of small Gaussian penalty functions:

2 1 2 2 || 4 ( ) (( ).( )) ( , ) [ exp{ }exp{ }] (2) 2( ) 2( ) i i i i i i t t i s s s s s s V t s W s s            
where t is the actual simulation time, i is the index representing the metadynamics steps, and

Δs ǁ i=|s i+1 -s i |.
The first exponential gives the hill's shape in the direction perpendicular to the trajectory, whereas the second exponential tunes the shape along the trajectory. The Gaussian height Wi and the width Δs ┴ are selected by the user and determine the accuracy in the free energy resolution and in the collective variable sampling, respectively. For our purposes, these values are set to 0.016 eV and 0.1 (in units of collective variables), respectively. In each metadynamics run, a new Gaussian contribution was added to V(s,t) every 150 steps. In our metadynamics simulations, two collective variables were employed, namely the coordination number between the two O atoms in H2O2 (CV1) and the distance between one O atom in H2O2

and the approaching C site on the CMCs surfaces (CV2). The coordination number (CN) is analytically described as:

( ) 1 1 1 ( / ) 1 ( , ) [ ] (3) 1 ( / ) A B p N N ij AB coord p q i i A ij AB r d N A B N r d        
The integer exponents p and q are set to 8. For CV1, the fictitious mass of the kinetic term was tuned to Mα = 60.0 a.u., and the harmonic coupling constants were set to kα = 1.0 a.u.. For CV2, the fictitious mass for the kinetic term is set to Mα = 15.0 a.u., and the harmonic coupling constants were set to kα = 0.25 a.u.. We remark that the mathematical function (3) of the CN is positive defined. Yet, we make use of the variational formulation of the metadynamics approach.

According to this formulation, the collective variables CV(t) evolve with their own Euler-Lagrange equations of motion. 26 This motion is constrained by a bounding harmonic potential k(CV(t) -Ncoord(A,B)) 2 along the (meta)dynamical evolution and for this reason CV(t) oscillates around the (never negative) Ncoord(A,B), but occasionally, when Ncoord(A,B) = 0, can assume instantaneously slightly negative values which can appear in the reconstructed free energy plot.

Free energy variations for each elementary step were computed by auxiliary DFT calculations done with the Quantum ESPRESSO package. 27 The adsorption energies of ORR intermediates are calculated by using H2O(l) and H2(g) as references. Details for calculating limiting potentials are provided in the Supporting Information. The generalized gradient approximation of Perdew-Burke-Ernzerhof (GGA-PBE) was used for the exchange correlation functional 28 thus providing also a benchmark for different GGA (PBE and HCTH) prescriptions. Spin unrestricted calculations were done in all cases. The kinetic energy cutoffs for the wavefunction, expanded in a plane-wave basis set, and the charge density were set to 35 Ry and 350 Ry, respectively.

RESULTS AND DISCUSSION

Active sites and mechanisms for H2O2 formation on CMCs.

Chemical reactions are driven by both thermodynamics and kinetics. In the present study, we provide an insight into both thermodynamics and kinetics for the indirect four-electron ORR.

This indirect four-electron ORR consists of a two-electron ORR responsible for the formation of H2O2 and a subsequent two-electron ORR for the conversion of H2O2 into water. The ORR at the cathode is a three-phase reaction, which was thought to be activated by O2 adsorption onto the exposed surface of noble metals (Pt, Pd-Hg and Ag-Hg). 29,30 This is indeed the case in Pt nanoclusters. 31 Yet, our recent investigation indicated that the barrier for O2 activation by abstraction of a hydrogen site on metal free CMCs can be lower than that for O2 adsorption at least on surfaces. This, in turn, can result in a new mechanisms for H2O2 formation on CMCs. 16,32 As the O2 adsorption mechanism has been extensively studied, here we focus on the possible H abstraction mechanism and on the active sites at the surface of CMCs systems which are prone to form H2O2.

The CMCs systems used in our simulations are the N doped and Fe based CMCs shown in

           
where the symbol * indicates a possible active site on the CMCs surface. For the H abstraction, the first step, involving an electron transfer, is the hydrogenation of the CMCs surface. Then, the second step, also involving an electron transfer process, is the formation of an OOH - intermediate. 9,16 Finally, the OOH -intermediate is expected to react with H + to form a H2O2. The reaction free energy for the overall H2O2 formation is -1.40 eV. The free energy variation of the first step is the hydrogenation energy ΔG*H. The reaction energy of the third step is a constant having an experimental value of -0.69 eV under standard conditions. 35,36 Hence, the reaction free energy of the second step can be obtained by subtracting the free energy variations of the other two steps from the total reaction energy, specifically, -0.71-ΔG*H. Thus, the limiting potential, due to the rate determining step for the formation of H2O2 via H abstraction, can be expressed as: 16 * * { , 0.71+ } ( 7)

L H H U Min G G   
Here, UL is determined by the limiting step, corresponding to the lowest energy decrease. The limiting potential becomes then the highest potential below which all the free energies variations for each electrochemical steps are downhill, although not necessarily barrierless. This can be compared directly with the experimental half-wave potential. The free energy variation for each elementary step in the H2O2 formation and reduction processes can be obtained from the related total energies as described in the Supporting Information. Details about this mechanism have been reported by our previous study. 16 Here we focus specifically on the screening of the active sites at the surface of CMCs for H2O2 formation. More importantly, we pay special attention to the indirect four electron ORR by combining the H2O2 formation and reduction processes.

Following this roadmap, we started by computing the hydrogenation energies of the different sites on the structures shown in structures turn out to be -0.44 eV, -0.25 eV, and -0.34 eV for the most favorable sites, respectively. This implies that on these sites the H abstraction mechanism is thermodynamically possible.

To make this point clearer, in Figure 3 we report the scatter-plot data for H abstraction under standard conditions in terms of the so-called volcano plot. In any type of fuel cells and rechargeable batteries, a quantity of major interest is the limiting potential, mainly referred to charge and discharge stages. 37 As shown in the volcano plot of Figure 3, the limiting potentials for the SW-N1, SW-N2, and SW-N3 systems are 0.27 V, 0.25 V, and 0.34 V, respectively. We find that the maximum limiting potential for H abstraction mechanism is about 0.4 V, and the SW-N1, SW-N2, and SW-N3 systems are thermodynamically possible for H2O2 formation via the H abstraction mechanism. We remark that the computed limiting potential is clearly lower than the standard potential (0.70 V) for H2O2 formation due to the H abstraction process. This agrees with the experimental observations indicating that the ORR half-wave potential of metal free CMCs in acidic solution is less than 0.4 V. [38][39][40] As the O2 adsorption mechanism and H abstraction mechanism are kinetically competing processes, we proceeded to compute the H abstraction kinetic profiles for the three structures indicated above (SW-N1, SW-N2, and SW-N3), sketched in Figure 4. The H abstraction barrier for the SW-N3 structure is 0.52 eV, while the corresponding ones for the SW-N1 and SW-N2 structures are higher than 0.6 eV. We also checked the barrier for proton transfer from the solution to the catalytic surface of the SW-N3 system. These results are shown in Figure S2 of the Supporting Information. Such a barrier is lower than the corresponding one for H abstraction, as expected, and amounts to about 0.26 eV. Since the O2 adsorption and the H abstraction are competing processes, the H abstraction mechanism might occur provided that the related activation barrier is lower than 0.6 eV, as discussed above. As shown in Figure 4, the O2 adsorption barriers are larger than the H abstraction ones in all cases. Given the scenario depicted so far, only the SW-N3 system appears to be both thermodynamically and kinetically favorable for H2O2 formation via H abstraction.

Active sites and mechanisms for H2O2 reduction on CMCs.

The mechanisms for subsequent two-electron ORR responsible for the conversion of H2O2 into water has been rarely studied and is still unclear to a great extent. As a consequence, the overall indirect four-electron mechanism remains to be elucidated. For the subsequent reduction of H2O2 to a water molecule, at least three alternative mechanisms are possible.

                        
According to this reaction pathway, in the first elementary step the H2O2 hydrogen peroxide dissociates into an *OH intermediate anchored to the CMC surface and an OH -hydroxyl anion.

Then the OH -hydroxyl anion reacts with H + in solution and this step is exoenergetic, being characterized by an energy gain of -0.83 eV in standard condition. 35,36 We remark that this step does not involve any electron transfer process. The final step consists in the removal of the *OH intermediate adduct via an electron transfer process characterized by a free energy variation indicated as -ΔG*OH in the figure. The standard potential for H2O2 reduction to water is 1.77 V.

Thus, the overall free energy variation of the first step is given by -3.54+0.83+ΔG*OH. The maximum limiting potential for this specific reaction pathway is then 1.36 V, corresponding to ΔG*OH=1.36 eV.

The second reaction mechanism, OH radical mechanism, is shown in Figure 5. The elementary steps are as follows:

2 2( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2
                                    
In this pathway, the first step also generates an *OH bound to the catalytic surface. However, here the OH group released in solution is a radical instead of a hydroxyl anion, which makes this a non-electron transfer step. The second step is the removal of the *OH intermediate from the CMC surface, which is accompanied by a free energy variation -ΔG*OH. Assuming that the OH group binds to the same catalytic site, the third step would be the formation of another *OH, whereas the fourth step is again represented by the removal of the *OH intermediate. The related free energy variations for the third and fourth steps are -2.72+ΔG*OH and -ΔG*OH, respectively.

In this case, the maximum limiting potential is 0.82 V, which allows to infer that ΔG*OH=0.82 eV.

The third reaction pathway is the H2O mechanisms, involving the generation of an H2O molecule and the binding of an *O intermediate to the surface of the system. discussed in our former work. 9 In this case, the maximum limiting potential turns out to be 1.62

2 2( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 ( ) ( ) 2 (
                       
V, corresponding to ΔG*OH=1.62 eV. Actually, the value of ΔG*OH falls usually in the range 0.2~1.3 eV, and the formation free energy of ΔG*OH can barely approach 1.62 eV. 9 Qualitatively speaking, the decrease in free energy of the third step is smaller than that of the second step. S3 for the SW-N3N3' structure as an illustrative example. Although we mainly focus on metal free CMCs, the free energy profile for a FeN4 center was also studied since Fe containing CMCs represent one of the most common systems used worldwide and are believed to be better performing than metal-free CMCs for ORR in acidic solution. Yet, Fe containing CMCs usually suffer from durability issues due to the formation of radicals during the ORR process as mentioned above. 14 Here we take iron porphine chloride (FePCl) as an example as the FeN4Cl center is a typical ORR active center. 44 For the G-N structure, the H2O2 moiety undergoes dissociation into two OH groups in water solution upon overcoming an activation barrier of 0.90 eV. Such a barrier is rather high, thus making very unlikely the hydrogen peroxide reduction via this reaction pathway. For this reason the G-N structure will not be considered in the ongoing discussion. However, for the FePCl center the calculated barrier is rather modest and amounts to just 0.31 eV as shown in Figure S4, which confirms, indeed, the better performance of this system for H2O2 dissociation. As shown in Figure S1, the H2O2 hydrogen peroxide decomposes into an *OH intermediate and an OH group. The OH group can be an OH -anion or a radical OH, which corresponds obviously either to the OH -ion or the OH radical mechanism. The OH radical can be realized with a low kinetic barrier on the FePCl center in comparison with a metal free case when the limiting potential is lower than 0.82 V as indicated above, but this process would destabilize and eventually destroy the catalytic system. The barrier for H2O2 dissociation in the absence of any catalyst was probed by BME simulations both in the gas phase and in aqueous solution. Results are summarized in Figure S5 of the Supporting Information. The resulting free energy barrier for the cleavage of the HO-OH bond in aqueous solution is 1.89eV, in good agreement with former results (1.90 eV). 45 As expected, the analogous free energy barrier in gas phase is larger and amounts to 2.19

eV. These values are significantly higher than the cases in which CMCs exert their catalytic role, corroborating the importance of CMCs for H2O2 reduction.

In the case of the SW-N3N3' system, the metadynamics simulations indicated that H2O2 can bind to the catalytic surface to form an *O intermediate adduct, accompanied by the release of one H2O molecule in solution (see Figure S1). This promotes the H2O mechanism discussed above. The corresponding H2O2 dissociation barrier is 0.62 eV. Here, we summarize that the OH radical can only be generated with a low kinetic barrier on metal doped CMCs such as the FePCl center (0.31 eV). For metal free CMCs, the more favorable kinetic pathway is to form an *O and H2O, analogous to that for SW-N3N3' system (0.62 eV). This may explain the long durability of metal free CMCs, while Fe containing CMCs are easily destroyed by radicals.

As the limiting steps for the two-electron H2O2 reduction are mainly determined by the nitrogen doped Stone-Wales defect, as found in former works. 9,46 We remark that the volcano plot shown in Figure 3 for indirect ORR refers specifically to the reaction in acidic solutions. This is substantially different from the direct ORR occurring in alkaline solution previously reported. 47 3.3 Different pathways for direct and indirect ORR.

So far, we have clarified the mechanisms and identified the active sites for the formation and reduction of H2O2 on CMCs in acidic solutions. On these grounds, we can conclude that three different routes are viable reaction pathways for ORR on CMCs in acidic solutions. This is already an indication for suppressing 2e -ORR, avoiding OH production on CMCs and promoting the direct four electrons ORR. In the first one, if the H abstraction barrier is much higher than that for O2 adsorption, only a direct four-electron ORR can occur via the O2 adsorption mechanism. In this case, the half-wave potential can be as high as 0.8 V as previously discussed. 9 In the second route the H abstraction barrier is comparable to that for O2 adsorption. Thus, if the H2O2 reduction barrier is much higher than that for O2 adsorption / H abstraction, both direct four-electron via O2 adsorption and two-electron ORR processes via H abstraction can be realized. Then, the overall ETN is expected to range between 2 and 4 and the linear sweep voltammograms (LSV) curve should be different with respect to the pure direct four-electron ORR. On the other hand, if the H2O2 reduction barrier is considerably lower than that for O2 adsorption / H abstraction, both direct and indirect four-electron ORR can occur. If this is the case, the ETN should also amount to about 4. In the third viable route, the H abstraction barrier can be much lower than that for O2 adsorption. When this happens, if the H2O2 reduction barrier is high, then only a two electron process can occur. The half-wave potential is then about 0.4 V and the ETN is 2. If the H2O2 reduction barrier is also low, an indirect four-electron ORR can be realized. This would be indicated by a low half-wave potential and a value of the ETN equal to 4.

All the direct and indirect ORR pathways discussed here are summarized in Figure 7. Finally, we call the attention on the fact that the formation of the H2O2 hydrogen peroxide is determined by the H abstraction, while the H2O2 reduction is determined by the removal of *OH.

As metal free CMCs usually have high O2 adsorption barriers, the ORR better performance would correspond to the third route. However, it would suffer from a low half-wave potential whether the ETN is 2 or 4 as discussed above. Thus, for metal free systems the ORR in acidic solutions is characteristic by a low half-wave potential (~0.4 V), in full agreement with experiments. For non-noble metal doped CMCs, such as Fe containing CMCs for instance, the first and second cases should be the dominating pathways because the O2 adsorption barriers are generally low for metal doped CMCs. However, in this case, the drawback is the generation of OH radicals, discussed in the OH radical mechanism for H2O2 reduction, eventually resulting in a destabilization of the catalytic system. We can remark also that the H2O2 can be produced via both an adsorption and H abstraction mechanism. Fortunately, the critical half-wave potential for OH radical generation is rather high (0.82 V). Note also that the maximum half-wave potential obtained in experiments is also around 0.8 V (the onset potential is generally about 0.1-0.2 V higher than half-wave potential). So, a high efficiency CMCs catalyst for ORR in acidic solution should 1) show strong binding to H or hindering bind to H to suppress the H abstraction mechanism; 2) show low activation barriers for O2 adsorption to improve the direct four-electron ORR; 3) be characterized by a half-wave potential lower than 0.82 V to avoid generation of OH radicals. This provides a roadmap and potentially a new strategy for developing CMCs in acidic solutions.

Conclusions

By combining first principles molecular dynamics, enhanced with free energy sampling techniques, and static DFT simulations, we clarified the detailed mechanisms and identified the active sites responsible for the indirect four-electron ORR and related two-electron ORR on CMCs in acidic solution. The indirect four-electron ORR is a combination of a first two-electron ORR resulting in the formation of H2O2 and a second one responsible for the H2O2 reduction.

Upon careful identification of the active sites for H2O2 formation on the CMCs surfaces via H abstraction, we provide an atomic-level insight and explain the origin of the low half-wave potential (around 0.4 V) for two-electron ORR. We find that the SW-N3 structure is a favorable system for H2O2 formation via an H abstraction mechanism. We also identify three different reaction mechanisms for the H2O2 reduction on CMCs, specifically (i) a mechanism involving the formation of OH -, (ii) a second one passing across the formation of an OH radical and a (iii) third one resulting in the production of H2O. For all these mechanisms, the rate determining step for H2O2 reduction is the removal of an *OH intermediate generally anchored to the CMC surface. The SW-N3N3' structure is remarkably favorable for H2O2 reduction via the H2O

production mechanism with a limiting potential of 1.06 V. The FePCl center is favorable to the OH radical mechanism, and the critical potential for OH radical generation is 0.82 V. The OH radical is kinetically much easier to generate on Fe containing structures than that for the N doped cases. Being this a parasite process in terms of the stability of a catalytic system, this may explain why Fe containing catalysts are easily destroyed and pathologically affected by durability issues, as opposed to the stability displayed by metal free CMCs upon ORR processes.

The whole indirect four-electron processes have been summarized and directly related to the half-wave potential and ETN experimentally measurable. The results indicate that the performance of ORR in acidic solutions can be improved by suppressing undesired 2e -H2O2

Figures 1

 1 Figures 1 and S1, respectively. The G-N and G-NN AB structures indicate the presence of either

Figure 1 .

 1 The results are summarized in Figure2(a). Namely, simulations show that the hydrogenation step is endothermic for the G-N, G-NNAB, and SW-N3N3' structures; this, in turn, suggests that the H abstraction mechanism is unlikely to occur for these structures. The hydrogenation energies calculated for the SW-N1, SW-N2, and SW-N3

  This is clearly identified from the metadynamics simulation of the SW-N3N3' system as discussed below, where the H2O2 hydrogen peroxide is decomposed into an *O intermediate and one H2O molecule in the first step. The second step is the conversion of the *O intermediate into an *OH intermediate. The final third step is the removal of the *OH intermediate from the CMC surface. The corresponding free energy variations for the second and third steps are ΔG*OH-ΔG*O and -ΔG*OH, respectively. The linear relationship for *O and *OH is "ΔG*O=2ΔG*OH+0.3" as

Therefore, the limiting 2 .

 2 potential is determined by the removal of the *OH intermediate. Further details about all the three reaction pathways are provided in the Supporting Information. The corresponding limiting potentials are mainly determined by the bond strength of the *OH intermediate formed on the CMCs surface. To make this statement more quantitative, we computed the formation free energy of the *OH intermediate; the results are presented in Figure 2(b) for all the structures considered in this work. The formation energies (ΔG*OH) for the G-N, G-NN AB , SW-N1, SW-N2, SW-N3 and SW-N3N3' systems are calculated according to our proposed procedure 9 and all values are in the range 0.2 -1.3 eV, as shown in Figure 2(b).The kinetics for H2O2 reduction are extracted from the free energy landscape provided by the metadynamics simulations. Specifically, the two CVs used to inspect the reaction pathway are the coordination number between the two O atoms in H2O2 (CV1) and the distance between one of these O atoms and the approaching C site on the CMCs surface (CV2). The specific C sites selected here is based on our extensive investigations done in former works dedicated exactly to the identification of the possible active sites.9,[33][34][41][42][43] More precisely, the C sites selected here are the ones showing the highest density of states (DOS) just below the Fermi level, thus allowing for electron donation, and the strongest binding energies for the intermediates as shown in Figure On these grounds, we focus here on reaction pathways for H2O2 reduction realized by these C sites and we select our collective variables accordingly. The free energy landscapes and related barriers obtained within this approach for the G-N and SW-N3N3' systems are shown in Figure 6 (a) and (b), respectively. The configurations of the initial and final states of the H2O2 reduction process are shown in Figure S1 of the Supporting Information, while the evolution of the CVs along the metadynamics trajectory are reported in Figure

  removal of the *OH intermediate adduct bound to the CMCs surface, the volcano plot for the H2O2 reduction of Figure 3 (b) is based on the use of the quantity ΔG*OH as a descriptor. Here we set the maximum limiting potential to 1.36 V according to the second reaction pathway. This choice is based on our finding that the ΔG*OH is mainly in the range of 0.2 eV~1.3 eV. As shown in Figure 3 (b), the G-N and G-NN AB systems can approach the maximum limiting potential, but the doped perfect graphene structures show evident difficulties in promoting the dissociation of the H2O2 molecule during the first steps, as all metadynamics simulations evidenced. For the SW-N1, SW-N2, and SW-N3 structures, the H2O2 reduction limiting potentials are relatively low, namely around or lower than 0.5 V. Interestingly, the SW-N3N3' structure shows both a relatively high limiting potential (1.06 V) and an energetically affordable H2O2 dissociation barrier of 0.62 eV. The SW-N3N3' structure is also the most stable among all possible two

Figure 1 .

 1 Figure 1. Atomic structures of (a) G-N, (b) G-NN AB , (c) SW-N1, (d) SW-N2, (e) SW-N3 and (d) SW-N3N3' considered in the present work for the ORR processes. Grey and blue balls indicate C and N atoms, respectively.

Figure 2 .

 2 Figure 2. Formation energies of (a) *H and (b) *OH for G-N, G-NN AB , SW-N1, SW-N2, SW-N3 and SW-N3N3' structures. The asterisk '*' denotes a specific site on carbon materials catalysts as shown in Figure 1.

Figure 3 .

 3 Figure 3. Volcano plot for (a) H2O2 formation on SW-N1, SW-N2, and SW-N3 structures, and (b) H2O2 reduction on G-N, G-NN AB , SW-N1, SW-N2, SW-N3 and SW-N3N3' structures.

Figure 4 .

 4 Figure 4. H abstraction and O2 adsorption free energy profiles for the SW-N1, SW-N2, and SW-N3 structures. The O2 adsorption free energy profiles reported in this figure refer to the results presented in Ref. 9. Symbols refer to the sampled points, while lines are intended as a guideline to the eye.

Figure 5 .

 5 Figure 5. H2O2 reduction pathways for (a) the OH -hydroxyl anion mechanism, (b) the OH radical mechanism, and (c) the H2O mechanism.

Figure 6 .

 6 Figure 6. H2O2 activation free energy landscapes for (a) the G-N and (b) SW-N3N3' structures. CV1 is the coordination number between the two O atoms in H2O2 molecule, whereas CV2 is the distance between O atom in H2O2 molecule and the C a site in G-N and SW-N3N3' structures as explained in the text.

Figure 7 .

 7 Figure 7. Pathways for direct and indirect ORR via O2 adsorption mechanisms and hydrogen abstraction mechanisms. E(O2,ad), E(Habs) and E(H2O2) denotes for activation barriers for O2 adsorption, hydrogen abstraction and H2O2 reduction, respectively.

  These have been explored according to the results of metadynamics simulations discussed below i.e. the H2O2 molecule can be decomposed into (*OH + OH) or (*O + H2O) at the first step depend on the specific catalytic surface. The produced OH group in solution can be OH -ion or OH radical, which leads to two different mechanisms. We name them OH -ion and OH radical mechanisms, respectively. The third one is named H2O mechanism as H2O is produced in the first step. These alternative

	2 2( ) * H O aq	( ) 2(H aq 	e 	( ) ) *OH OH aq 		( ) 2H aq 	e 	(8)
	( ) *OH OH aq 		( ) 2H aq	e	2 ( ) *OH H O aq	( ) H aq	e	(9)
	2 ( ) *OH H O aq	( ) H aq	e	2 ( ) * 2H O aq	(10)

reaction mechanisms are sketched in Figure

5

. The first one, termed OH -ion mechanism hereafter, is composed of the following elementary steps:

of Science and Technology (JAIST). M.B. thanks Pôle HPC -Equipex Equip@Meso, University of Strasbourg, GENCI under allocation DARI-A2 A0020906092 and the LaBex "Nanoparticles Interacting with their Environment" ANR-11-LABX-0058_NIE.

ABBREVIATIONS

Oxygen reduction reaction (ORR); Carbon materials catalysts (CMCs).

formation via H abstraction mechanism, avoiding OH production via OH radical mechanism during H2O2 reduction, and promoting direct 4e -ORR by decreasing O2 adsorption battier. This not only disclose a new way to improve and design high efficient CMCs for ORR in acid solutions but also provides a guideline for future development of ORR catalysts for practical applications in acidic environments.


Keywords: Workflows, graphes de tâches, ordonnancement, hétérogénéité, grille de calcul

Allocation d'une série de graphes de tâches sur une grille de calcul. Résumé : Dans ce rapport, nous nous intéressons à l'ordonnancement d'applications de type workflow sur une grille de calcul. Dans notre modèle, une telle application est formée d'une grande collection de données qui doivent toutes être traitées avec le même graphe de tâche. Nous avons donc une séries de graphes de tâches à effectuer. Nous recherchons une solution efficace pour le débit et la latence, tout en nous interdisant le recours à un contrôle trop complexe. C'est pourquoi nous nous concentrons sur les solutions utilisant une seule allocation. Nous présentons un algorithme utilisant la programmation linéaire mixte qui calcule une allocation de débit optimal, et ce pour différentes politiques de routage, en fonction de la latitude que nous avons sur les communications. Puis, grâce à des simulations, nous comparons ces allocations à des heuristiques de référence. Nos résultats montrent que notre algorithme est toujours capable de trouver une allocation de bon débit et de latence limitée et qu'il surpasse les résultats des heuristiques de référence, en particulier lorsque les temps communications sont prépondérants.

Introduction

Computing Grids gather large-scale distributed and heterogeneous resources, and make them available to large communities of users [START_REF] Foster | The Grid 2: Blueprint for a New Computing Infrastructure[END_REF]. Such Grids enable large applications from various scientific fields to be deployed on large numbers of resources. These applications come from domains such as high-energy physics [START_REF] Chervenak | The data grid: Towards an architecture for the distributed management and analysis of large scientific datasets[END_REF], bioinformatics [START_REF] Thomas | Taverna: a tool for the composition and enactment of bioinformatics workflows[END_REF], medical image processing [START_REF] Germain | Gridenabling medical image analysis[END_REF], etc. Distributing an application on such a platform is an increasingly complex duty. As far as performance is concerned, we have to take into account the computing requirements of each task, the communication volume of each data transfer, as well as the platform heterogeneity: the processing resources are intrinsically heterogeneous, and run different systems and middlewares; the communication links are heterogeneous as well, due to their various bandwidths and congestion status. In this paper, we investigate the problem of mapping an application onto the computing platform. We are both interesting in optimizing the performance of the mapping (that is, process the data as fast as possible), and to keep the deployment simple, so that we do not have to deploy complex control softwares on a large number of machines.

Applications are usually described by a (directed) graph of tasks, what is usually called a workflow in the Grid literature. The nodes of this graph represent the computing tasks, while the edges between nodes stand for the dependences between these tasks, which are usually materialized by files: a task produces a file which is necessary for the processing of some other task.

In this paper we consider Grid jobs involving a large collection of input data sets that must all be processed by the same application. In other words, the Grid jobs we consider are made of the same workflow applied to a large collection of different input data sets. We can evenly consider that we have a large number of instances of the same task graph to schedule. Such a situation arises when the same computation must be performed on independent data [START_REF] Lee | Exploring protein fold space by secondary structure prediction using data distribution method on grid platform[END_REF], independent parameter sets [START_REF] Pitt-Francis | Enabling computer models of the heart for high-performance computers and the grid[END_REF], or independent models [START_REF] Daniel Haunschild | Investigating the dynamic behavior of biochemical networks using model families[END_REF]. A classical example lies in image processing: the data set is a large number of input images, and all images have to be processed the same way, for example by using several consecutive filters and encoders.

In this paper, we concentrate on how to map several instances of a same workflow onto a computing platform, that is, on how to decide on which resource a task has to be processed, and on which route a file has to be transfered, if we assume that we have some control on the routing mechanism. We will study several scenarios, with different routing policies.

We start by motivating our problem (Section 2). Then we formally define our problem (Section 4) and describe our solution (Section 5). Finally, we assess the quality of these solutions through simulations (Section 6) and conclude (Section 7).

Problem motivation

In this section we motivate the application model we work with.

Dynamic scheduling vs. static scheduling

Many scheduling strategies use a dynamic approach: task graphs, or even tasks, are processed one after the other. This is usually done by assigning priorities to waiting tasks, and then by allocating resources to the task with highest priority, as long as there are free resources. This simple strategy is the best possible in some cases: (i) when we have very little knowledge on the future workload (i.e., the tasks that will be submitted in the near future, or released by the processing of current tasks), or (ii) under a very unstable environment, where machines join and leave the system with a high churn rate. However, in the case of scheduling a series of workflows, the number of workflows is assumed to be large, and the number of resulting tasks is even larger. If a dynamic scheduler is used to schedule this large number of tasks, it would result in a large processing time by the scheduler, or even worse, in an overflow in the task buffer. This way, the scheduler would be unavailable for other users until all our tasks are scheduled, which would be unfair in a multi-user environment.

On the contrary to the typical use case of dynamic schedulers, we have much knowledge on the system when scheduling several instances of a workflow. First, we can take advantage of the regularity of the workflows: the input is made of a large collection of data sets that will result in the same task graph. Second, the computing platform is considered to be stable enough so that we can use performance measurement tools like NWS [START_REF] Wolski | The network weather service: a distributed resource performance forecasting service for metacomputing[END_REF] in order to get some information on machine speeds and link bandwidths. Taking advantage of this knowledge, we aim at using static scheduling techniques, that is to anticipate the mapping and the scheduling of the whole workload at its submission date.

Data parallelism vs. control parallelism

In the context of scheduling series of task graphs, we can take advantage of two sources of parallelism to increase performance. First, parallelism comes from the data, as we have to process a large number of instances. Second, each instance consists in a task graph which may well include some parallelism: some tasks can be processed simultaneously, or the processing of consecutive tasks of different instances can be pipelined, using some control parallelism. In such a context, several scheduling strategies may be used.

We may only make use of the data parallelism. In this case, the whole workflow corresponding to the processing of a single input data set is executed on a single resource, as if it was a large sequential task. Different workflow instances are simultaneously processed on different processors. This is potentially the solution with the best degree of parallelism, because it may well use all available resources. This imposes that all tasks of a given instance are performed on each processor, therefore that all services must be available on each participating machine. However, it is likely that some services have heterogeneous performance: many legacy codes are specialized for specific architectures and would perform very poorly if run on other machines. Some services are even likely to be unavailable on some machines. In the extreme, most specified case, it may happen that no machine can run all services; in that case the pure data-parallelism approach is infeasible. More-INRIA over, switching, on the same machine, from one service to another may well induce some latency to deploy each service, thus leading to a large overhead. At last, a single input data set may well produce a large set of data to process or require a large amount of memory. Processing the whole workflow on a single machine may lead to a large latency for this instance, and may even not be possible if the available storage capacity or memory of the machine cannot cope with the workflow requirements. For these reasons, application workflows are usually not handled using a pure data-parallelism approach.

Another approach consists in simultaneously taking advantage of both data and control parallelism. We have previously studied this approach [START_REF] Beaumont | Scheduling strategies for mixed data and task parallelism on heterogeneous clusters[END_REF][START_REF] Beaumont | Assessing the impact and limits of steady-state scheduling for mixed task and data parallelism on heterogeneous platforms[END_REF] and proved that in a large number of cases, when the application graph is not too deep, we can compute an optimal schedule, that is a schedule which maximizes the system throughput. This approach, however, asks for a lot of control as similar files produced by different data sets must follow different paths in the interconnection network.

In this paper, we focus on a simpler framework: we aim at finding a single mapping of the application workflow onto the platform, with good performance. This means that all instances of a given task must be processed on the same machine. Thus, the corresponding service has to be deployed on a single machine, and all instances are processed the same way. Thus, the control of the Grid job is much simpler, and the number of needed resources is kept low.

Steady-state operation and throughput maximization

As in our previous work for scheduling application graphs on heterogeneous platforms, this study relies on steady-state scheduling. The goal is to take advantage of the regularity of the series of workflows; as we consider that the Grid job input is made of a large number of data sets which should be processed using the same task graph, we relax the scheduling problem, and consider the steady-state operation: we assume that after some transient initialization phase, the throughput of each resource will become steady.

In scheduling, the classical objective is to minimize the running time of the job, or makespan. However, by using steady-state techniques, we relax this objective and concentrate on maximizing the system throughput. Then, the total running time is composed of one initialization phase, a steady-state phase, and a clean-up phase. As initialization and clean-up phases do not depend on the total number of instances, we end up with asymptotically optimal schedules: when the number of instances becomes large, the time needed to perform initialization and clean-up phases becomes negligible in front of the overall running time.

Related work

In this section, we report the related work link to workflow scheduling and steady-state relaxation.

Workflows. Managing and scheduling workflow on computing Grids is the subject of a wide literature. Many middlewares are develop in order to manage workflows on the Grid. One of the most comprehensive, as far as scheduling is concerned in probably MO-TEUR [START_REF] Glatard | Flexible and efficient workflow deployement of data-intensive applications on grids with MO-TEUR[END_REF]. All these tools usually include scheduling heuristics to map the tasks of the workflow onto the available resources. These heuristics are inherited from DAG scheduling, and more or less adapted to Grid environment to cope with its intrinsic heterogeneity. Among other, these techniques make use of list scheduling heuristics [START_REF] Chrétienne | Scheduling Theory and its Applications[END_REF], clustering [START_REF] Michael | Task clustering and scheduling for distributed memory parallel architectures[END_REF] and task duplication [START_REF] Ahmad | On exploiting task duplication in parallel program scheduling[END_REF]. Many meta-heuristics are derived by assembling various scheduling heuristics. Note that some workflow may include parallel tasks, requiring a (fixed or not) number of processors, and some scheduling techniques take advantage of this knowledge [START_REF] Suter | From heterogeneous task scheduling to heterogeneous mixed parallel scheduling[END_REF]. Since workflows and Grid computing are very popular nowadays, we must be careful with the vocabulary: a workflow is most of the time considered to be a single graph of tasks with dependences, but sometimes also denotes a pipeline usage of this task graph, like our "series of workflow" setting.

Steady-state.

Minimizing the makespan, i.e., the total execution time, is a NP-hard problem in most practical situations [START_REF] Garey | Computers and Intractability, a Guide to the Theory of NP-Completeness[END_REF][START_REF] Shirazi | Scheduling and load balancing in parallel and distributed systems[END_REF][START_REF] Chrétienne | Scheduling Theory and its Applications[END_REF], while it turns out that the optimal steady-state schedule can often be characterized very efficiently, with low-degree polynomial complexity.

The steady-state approach has been pioneered by Bertsimas and Gamarnik [START_REF] Bertsimas | Asymptotically optimal algorithms for job shop scheduling and packet routing[END_REF]. It has been used successfully in many situations [START_REF] Beaumont | Steady-state scheduling on heterogeneous clusters[END_REF]. In particular, steady-state scheduling has been used to schedule independent tasks on heterogeneous tree-overlay networks [START_REF] Beaumont | Bandwidth-centric allocation of independent tasks on heterogeneous platforms[END_REF][START_REF] Banino | Scheduling strategies for master-slave tasking on heterogeneous processor platforms[END_REF]. The steady-state approach has also been used by Hong et al. [START_REF] Hong | Distributed adaptive task allocation in heterogeneous computing environments to maximize throughput[END_REF] who extend the work of [START_REF] Beaumont | Bandwidth-centric allocation of independent tasks on heterogeneous platforms[END_REF] to deploy a divisible workload on a heterogeneous platform.

Notations, hypotheses, and problem complexity 4.1 Platform and application model

We denote by G P = (V P , E P ) the undirected graph representing the platform, where V P = {P 1 , . . . , P p } is the set of all processors. The edges of E P represent the communication links between these processors. The maximum bandwidth of the communication link P q → P r is denoted by bw q,r . Moreover, we suppose that processor P q has a maximum incoming bandwidth B in q and a maximum outgoing bandwidth B out q . Figure 1(a) gives an example of such a platform graph.

We use a bidirectionnal multiport model for communications: a processor can perform several communications simultaneously. In other words, a processor can simultaneously send data to multiple targets and receive data from multiple sources, as soon as the bandwidth limitation is exceeded neither on links, nor on incoming or outgoing ports.

A path from processor P q to processor P r , denoted P q ❀ P r , is a set of adjacent communication links going from P q to P r . We denote by G A = (V A , E A ) the Directed Acyclic application Graph (DAG), where V A = {T 1 , . . . , T n } is the set of tasks, and E A represents the dependences between these tasks, that is, F i,j = (T i , T j ) ∈ E A is the file produced by task T i and consumed by task T j . The dependence file F i,j has size data i,j , so that its transfer through link P q → P r takes a time datai,j bwq,r . Computation task T k needs a time w i,k to be entirely processed by processor P i . Figure 1(b) gives an example of application graph. This last notation corresponds to the so-called unrelated-machines model: a processor can be fast for a given type of task and slow for another one. Using these notations, we can model the benefits which can be drawn on some specific hardware architectures by specially optimized tasks. For example, a Cholesky factorization can be 5.5 times faster when using a GeForce 8800GTX graphic card than when using only the CPU, while a LU factorization is only 3 times faster in the same conditions [START_REF] Volkov | qr and cholesky factorizations using vector capabilities of gpus[END_REF]. Even higher speed-ups can be reached on some graph applications, as high as 70 [START_REF] Harish | Accelerating large graph algorithms on the gpu using cuda[END_REF]. Unrelated performance may also come from memory requirements. Indeed, a given task requiring a lot of memory will be completed faster when processed by a slower processor with a larger amount of memory. Grids are often composed of several clusters bought over several years, thus with very different memory capacities, even when processors are rather similar.

Remarks on platform modeling. The way we model the platform aims at taking into account all processing units and routers in the set of nodes V P . Modelling routers that should not be enrolled in computations is straightforward by using infinite computation time on these resources. Similarly, E P is the set of all physical links connecting nodes of V P . The advantage of this approach is that it allows us to correctly model potential congestions in the platform. Its main drawback is the difficulty of first acquiring a good model of the topology of the interconnection network, and then of instantiating its parameters. Although this is a tough problem, some solutions begin to arise [START_REF] Eyraud-Dubois | A first step towards automatically building network representations[END_REF].

Another way to consider the model is at an application-level. From this point of view, all nodes in V P are processors under our control -i.e, on which we can deploy our applicationand E P is the set of logical routes connecting these nodes (each logical route can then account for several physical communication links). This model has the advantage of the simplicity: it does not need information on physical links or routers which may be out of the scope of our knowledge. On the other hand, potential link contentions are not taken into account, precisely because physical links are not directly considered.

Allocations

As described in the introduction, we assume that a large set of input data sets has to be processed. These input data sets are originally available on a given source processor P source . Each of these data sets contains the data for the execution of one instance of the application workflow. For the sake of simplicity, we are looking for strategies where all tasks of a given type T i are performed on the same resource, which means that the allocation of tasks to processors is the same for all instances. We now formally define an allocation.

Definition 1 (Allocation). An allocation of the application graph to the platform graph is a function σ which associates:

• to each task T i , a processor σ(T i ) which processes all instances of T i ; • to each file F i,j , a set of communication links σ(F i,j ) which carries all instances of this file from processor σ(T i ) to processor σ(T j ).

A file F i,j may be transfered differently from σ(T i ) to σ(T j ) depending on the routing policy enforced on the platform. We distinguish three possible policies: Single path, fixed routing. The path for any transfer from P q to P r is fixed a priori. We do not have any freedom on the routing. This scenario corresponds to the classical case where we have no freedom on the routing between machines: we cannot change the routing tables of routers.

Single path, free routing. We can choose the path from P q to P r , but a single route must be used for all data originating from P q and targeting P r . This policy corresponds to protocols allowing us to choose the route for any of the data transfer, and allowing us to reserve some bandwidth on the chosen routes. Although current network protocols do not provide this feature, bandwidth reservation, and more generally resource reservation in Grid network, is the subject of a wide literature, and will probably be available in future computing platforms [START_REF] Foster | End-to-end quality of service for high-end applications[END_REF].

Multiple paths. Data from P q to P r may be split along several routes taking different paths. This corresponds to the uttermost flexible case where we can simultaneously reserve several routes and bandwidth fractions for concurrent transfers.

INRIA

Our three routing policies allow us to model a wide range of practical situations, current and future. The techniques exposed in Section 5 enable us to deal with any of this model and even with combinations of them.

In the case of single path policies, σ(F i,j ) is the set of the links constituting the path. In the case of multiple paths, σ(F i,j ) is a weighted set of paths {(w α , P a )}: for example σ(F 7,8 ) = {(0.1, P 1 → P 3 ), (0.9, P 1 → P 2 → P 3 )} means that 10% of the file F 7,8 go directly from P 1 to P 3 and 90% are transfered through P 2 .

F 2,4 F 3,5 T 1 T 6 T 4 , T 5 T 2 , T 3 P6 P5 P2 P1 P4 P3
(a) Single path, fixed routing. In Figure 2, we give, for each routing policy, an example of allocation of the application graph of Figure 1(b). In all cases, the mapping of tasks is the same: σ(T 1 ) = P 1 , σ(T 2 ) = σ(T 3 ) = P 2 , σ(T 4 ) = σ(T 5 ) = P 5 and σ(T 6 ) = P 6 . In Figure 2(a), the path for any transfer is fixed; in particular, all data from P 2 targeting P 5 must pass through P 3 . In Figure 2(b), we assume a free routing policy: we can choose to route some files via P 3 and some others via P 4 . At last, in Figure 2(c), we are allowed to use multiple paths to route a single file, which is done for file F 3,5 .

F 2,4 F 3,5 T 1 T 6 T 4 , T 5 T 2 , T 3 P6 P5 P2 P1 P4 P3 (b) Single path, free routing. F 2,4 1 2 F 3,5 1 2 F 3,5 T 1 T 6 T 4 , T 5 T 2 , T 3 P6 P5 P2 P1 P4 P3 (c) Multiple paths.
We can also look at our routing policies from an application point of view, as we did with the platform model. Then, from this perspective:

• At an application level, the standard policy would be single path with fixed routing:

the application is the passive subject of external routing decisions it does not try to interfere with.

• Under the single path with free routing policy, the application tries to improve its performance by circumventing some of the routing decisions. In practice, a communication is not made, at an application-level, directly from a participating node P q to a participating node P r , but may be routed through other participating nodes. The use of a single path for each transfer naturally derives from the need to have a simple control, and is in strict accordance with the use of a single allocation.

• Using several concurrent routes with the multiple paths policy may considerably increase the system's throughput if the application is communication intensive. However this would come at the price of more control on the participating nodes, because of the added necessity to split transfers into pieces to be separately routed.

Thus, our platform model and routing policies can be used to model a wide variety of situations, depending on how much control we are allowed, or we want, on the platform. Moreover, this is not a comprehensive list of possibilities, and we can well imagine to incorporate other policies, or a combination of them.

Throughput

We first formally define what we call the "throughput" of a schedule. Then, we will derive a tight upper bound on the throughput of any schedule.

The definition of throughput. We focus on the optimization of the steady state. Thus, we are not interested in minimizing the execution time for a given number of workflow instances, but we concentrate on maximizing the throughput of a solution, that is the average number of instances that can be processed per time-unit in steady-state. Definition 2. Assume that the number of instances to be processed is infinite, and note N (t) the number of instances totally processed by a schedule at time t. The throughput ρ of this schedule is

given by ρ = lim t→∞ N (t) t .
This definition is the most general one, as it is valid for any schedule. In this study, we are interested in very specific schedules, consisting of only one allocation. We now show how to compute an upper bound on the achievable throughput of a given allocation. We will later show that this bound is tight.

Upper bound on the achievable throughput. First, we consider the time spent by each resource on one instance of a given allocation σ. In other words, we consider the time spent by each resource for processing a single copy of our workflow under allocation σ.

• The computation time spent by a processor P q for processing a single instance is:

t comp q = i,σ(Ti)=Pq w i,q .
• The total amount of data carried by a communication link P q → P r for a single instance is d q,r = (i,j),Pq→Pr∈σ(Fi,j ) data i,j in the case of single-path policies, and it is

INRIA d q,r =
Fi,j (Wa,Pa)∈σ(Fi,j ) Pq→Pr∈Pa w a × data i,j in the case of the multiple-paths policy. This allows us to compute the time spent by each link, and each network interface, on this instance:

on link P q → P r : t q,r = d q,r /bw q,r ; -on outgoing interface of P q : t out q = r d q,r /B out q ; -on incoming interface of P q : t in q = r d r,q /B in q .

We can now compute the maximum time τ spent by any resource for the processing of one instance: τ = min min Pq {t comp q , t out q , t in q }, min Pq→Pr t q,r . This obviously gives us an upper bound on the achievable throughput: ρ ≤ ρ max = 1/τ . Indeed, as there is at least one resource which spends a time τ to process its share of a single instance, the throughput cannot be greater than 1 instance per τ units of time. We now show that this upper bound is achievable in practice, i.e., that there exists a schedule with throughput ρ max . In the following, we call "throughput of an allocation" the optimal throughput ρ max of this allocation. The upper bound is achievable. Here, we will only explain on an example how one can built a periodic schedule achieving the throughput ρ max . Indeed, we are not interested here in giving a formal definition of periodic schedules, or to formally define and prove schedules which achieve the desired throughput, as this goes far beyond the scope of this paper. The construction of such schedules, for applications modeled by DAGs, was introduced in [START_REF] Beaumont | Scheduling strategies for mixed data and task parallelism on heterogeneous clusters[END_REF], and a fully comprehensive proof can be found in [START_REF] Beaumont | Assessing the impact and limits of steady-state scheduling for mixed task and data parallelism on heterogeneous platforms[END_REF].

T 2 T 3 T 5 T 4 T 6 T 1 F 4,6 F 2,4 F 3,5 F 5,6 F 1,2 F 1,3
Figure 3 illustrates how to build a periodic schedule of period τ for the workflow described on Figure 1(b). Once the schedule has reached its steady-state, that is after 6τ in the example, during each period, each processor computes one instance of each task assigned to it. More precisely, in steady-state, during period k (k ≥ 6), that is during time-interval [kτ ; (k + 1)τ ], the following operations happens:

• P 1 computes task T 1 of instance k,

• P 1 sends F 1,2 and F 1,3 of instance k -1 to P 2 ,

• P 2 processes T 2 and T 3 of instance k -2,

• P 2 sends F 2,4 and F 3,5 of instance k -3 to P 5 (via P 3 ),

• P 4 processes tasks T 4 and T 5 of instance k -4,

• P 5 sends F 4,6 and F 5,6 of instance k -5 to P 6 , • P 6 processes task T 6 of instance k -6. One instance is thus completed after each period, achieving a throughput of 1/τ .

NP-completeness of throughput optimization

We now formally define the decision problem associated to the problem of maximizing the throughput.

Definition 3 (DAG-Single-Alloc). Given a directed acyclic application graph G A , a platform graph G P , and a bound B, is there an allocation with throughput ρ ≥ B? Theorem 1. DAG-Single-Alloc is NP-complete for all routing policies.

Proof. We first have to prove that the problem belongs to NP, that is that we can check in polynomial time that the throughput of a given allocation is greater than or equal to B. Thanks to the previous section, we know that this check can be made through the evaluation of a simple formula; DAG-Single-Alloc is thus in NP.

To prove that DAG-Single-Alloc is NP-complete, we use a reduction from the Minimum Multiprocessor Scheduling, known to be NP-complete [START_REF] Garey | Computers and Intractability, a Guide to the Theory of NP-Completeness[END_REF]. Consider an instance I 1 of Multiprocessor Scheduling, that is a set of n independent tasks T i,1≤i≤n and a set of m processors P u,1≤u≤m , where task i takes time t(i, u) to be processed on processor P u . The problem is to find a schedule with total execution time less than a given bound T . We construct a very similar instance of DAG-Single-Alloc:

• The application DAG is a simple fork, made of all tasks T i plus a task T 0 , root of the fork: for each 1 ≤ i ≤ n, there is an edge F 0,i , with data 0,i = 0.

• The platform consists of the same set of processors than I 1 , connected with a complete network where all bandwidths are equal to 1. The time needed to process task T i on processor P u is w i,u = t(i, u) for each 1 ≤ i ≤ n, and w 0,u = 0.

Note that communications need not being taken into account during performance evaluation, since all data sizes are null. Thus, this reduction applies to any routing policy. The throughput of an allocation is directly related to the total execution time of the set of tasks: an allocation has throughput ρ if and only if it completes all the tasks in time 1/ρ. Thus finding a schedule with completion time less than T is equivalent to finding an allocation with throughput greater than 1/T .

INRIA 5 Mixed linear program formulation for optimal allocations

In this section, we present a mixed linear program formulation that allows to find optimal allocation with respect to the total throughput.

Single path, fixed routing

In this section, we assume that the path to be used to transfer data from a processor P q to a processor P r is determined in advance; we have thus no freedom on its choice. We then denote by P q ❀ P r the set of edges of E P which are used by this path.

Our linear programming formulation makes use of both integer and rational variables. The resulting optimization problem, although NP-complete, is solvable by specialized softwares (see Section 6 about simulations). The integer variables can take 0 or 1 value. The only integer variables are the following: • y's variables which characterize where each task is processed: y k q = 1 if and only if task T k is processed on processor P q ; • x's variables which characterize the mapping of file transfers: x k,l q,r = 1 if and only if file F k,l is transfered using path P q ❀ P r ; note that we may well have x k,l q,q = 1 if processor P q executes both tasks T k and T l .

Obviously, these two sets of variables are related. In particular, for any allocation, x k,l q,r = y k q × y l r . This redundancy allows us to write linear constraints. Linear program (1) expresses the optimization problem for the fixed-routing policy. The objective function is to minimize the maximum time τ spent by all resources, in order to maximize the throughput 1/τ . The intuition behind the linear program is the following:

• Constraints (1a) define the domain of each variable: x, y lie in {0, 1}, while τ is rational.

• Constraint (1b) ensures that each task is processed exactly once.

• Constraint (1c) asserts that a processor can send the output file of a task only if it processes the corresponding task.

• Constraint (1d) asserts that the processor computing a task holds all necessary input data: for each predecessor task, it either received the data from that task or computed it.

• Constraint (1e) ensures that the computing time of a processor is no larger that τ .

• In Constraint (1f), we compute the amount of data carried by a given link, and the following constraints (1g,1h,1i) ensure that the time spent on each link or interface is not larger than τ , with a formulation similar to that of Section 4.3.

We denote ρ opt = 1/τ opt , where τ opt is the value of τ in any optimal solution of Linear Program [START_REF]SimGrid: a Generic Framework for Large-Scale Distributed Experimentations. Henri casanova and arnaud legrand and martin quinson[END_REF]. The following theorem states that ρ opt is the maximum achievable throughput.

Theorem 2. An optimal solution of Linear Program (1) describes an allocation with maximal throughput for the fixed routing policy.

(1)

                                                                             MINIMIZE τ UNDER THE CONSTRAINTS (1a) ∀F k,l , ∀P q ❀ P r , x k,l q,r ∈ {0, 1}, y k q ∈ {0, 1} (1b) ∀T k , Pq y k q = 1 (1c) ∀F k,l , ∀P q ❀ P r , x k,l q,r ≤ y k q (1d) ∀T l , ∀F k,l , ∀P r , y k r + Pq❀Pr x k,l q,r ≥ y l r (1e) ∀P q , T k y k q w q,k ≤ τ (1f) ∀P q → P r , d q,r = Ps❀Pt with Pq→Pr∈Ps❀Pt   F k,l x k,l s,t data k,l   (1g) ∀P q → P r , d q,r bw q,r ≤ τ (1h) ∀P q Pq→Pr∈E P d q,r B out q ≤ τ (1i) ∀P r Pq→Pr∈E P d q,r B in q ≤ τ
Proof. We first prove that for any allocation of throughput ρ, described by x's and y's variables, (x, y, τ = 1/ρ) satisfies the constraints of the linear program.

• All tasks of the workflow are processed exactly once, thus Constraint (1b) is verified.

• In any allocation, x k,l q,r = y k q × y l r , thus Constraint (1c) is verified.

• If y l r = 1, that is if P r processes T l , then P r must own all files F k,l . It can either have it because it also processed T k (in this case y k r = 1) or because it received it from some processor P q (and then x k,l q,r = 1). In both cases, Constraint (1d) is satisfied. • As the allocation has throughput ρ, it means that the occupation time of each processor, each link, and each network interface is at most 1/ρ, as explained in Section 4.3.

INRIA

This is precisely what is stated by Constraints (1e), (1g), (1h) and (1i). Thus, these constraints are satisfied.

As all allocations satisfy the constraint of the linear program, and ρ opt is the maximal value of the throughput under these constraints, then all allocations have a throughput smaller than, or equal to, ρ opt . We now prove that any solution of the linear program represents an allocation. We have to verify that all tasks are performed, all input data needed to process a task are correctly sent to the corresponding processor, and that the allocation has the expected throughput 1/τ .

• All tasks are processed exactly once due to constraint (1b).

• Thanks to Constraint (1c), a processor is allowed to send a file F k,l only if it processed T k .

• Thanks to Constraint (1d), a task T l with predecessor T k is processed on P r only if P r either processed T k , or received F k,l from some processor P q .

• Thanks to Constraints (1e), (1g), ( 1h) and (1i), we know that the maximum utilization time of any resource (processor, link or network interface) is at least equal to τ , the corresponding allocation has a throughput at most 1/τ .

Thus, any solution of the linear program describes a valid allocation. In particular, there is an allocation with throughput ρ opt .

Single path, free routing

We now move to the free routing setting. The transfer of a given file between two processors can take any path between these processors in the platform graph. We introduce a new set of variables to take this into account. For any file F k,l and link P i → P j , f k,l i,j is an integer value, with value 0 or 1: f k,l i,j = 1 if and only if the transfer of file F k,l between the processor processing T k to the one processing T l takes the link P i → P j . Using these new variables, we transform the previous linear program to take into account the free routing policy. The new program, Linear Program (2) has exactly the same equations than Linear Program (1) except for the following: 1) the new variables are introduced (Constraint (2a)); 2) the computation of the amount of data in Constraint (1f) is modified into Contraint (2f) to take into account the new definition of the routes; 3) the new set of constraints (2j) ensures that a flow of value 1 is defined by the variables f k,l from the processor executing T k to the one executing T k .

(2)

                       MINIMIZE τ UNDER THE CONSTRAINTS (2a) ∀F k,l , ∀P q ❀ P r , x k,l q,r ∈ {0, 1}, y k q ∈ {0, 1}, f k,l i,j ∈ {0, 1} (2f) ∀P q → P r , d q,r = F k,l f k,l i,j data k,l (2j) ∀P q , ∀F k,l , Pq→Pr f k,l q,r - P r ′ →Pq f k,l r ′ ,q = Pt x k,l q,t - Ps x k,l s,q
AND (1b), (1c), (1d), (1e), (1g), (1h), (1i)

In the following lemma, we clarify the role of the f variables.

Lemma 1. Given a file F k,l , the following two properties are equivalent

(i) ∀P q , Pq→Pr f k,l q,r - P r ′ →Pq f k,l r ′ ,q =      1 if P q = P prod -1 if P q = P cons 0 otherwise
(ii) the set of links P q → P r such that f k,l q,r = 1 defines a route from P prod to P cons .

Proof. The result is straightforward since Property (i) is a simple conservation law of f quantities. Note that this route may include cycles. These cycles do not change the fact that the links can be used to ship the files from P prod to P cons , but the time needed for the transportation is artificially increased. That is why in our experiments these cycles are sought and deleted to keep routes as short as possible.

The following theorem states that the linear program computes an allocation with optimal throughput: again, we denote ρ opt = 1/τ opt , where τ opt is the value of τ in any optimal solution of this linear program. Theorem 3. An optimal solution of Linear Program (2) describes an allocation with optimal throughput for the free routing policy.

Proof. Similarly to the proof of Theorem 2, we first consider an allocation, define the x, y and f variables corresponding to this allocation, and prove that they satisfy the constraints of the linear program.

• Since f k,l i,j describes if transfer F k,l uses link P i → P j , all constraints excepted (2j) are satisfied by the same justifications as in Theorem 2.

• In any allocation, file F k,l must be routed from the processor executing T k to the processor executing T l (provided that these tasks are excepted by different processors).

INRIA

Thus, f define a route between those two processors. Then, we note that Pt x k,l q,t -

Ps x k,l s,q =      1 if P q executes T k and not T l -1 if P q executes T l and not T k 0 otherwise
If T k and T k are executes on the same processor, all corresponding f and x variables are equal to 0. Thus, thanks to Lemma 1, all Constraints (2j) are verified.

We now prove that any solution of Linear Program (2) defines a valid allocation with throughput ρ under the Free routing policy.

• As in the proof of Theorem 2, we can prove that x and y variables define an allocation with throughput ρ.

• As above, we note that for a given file

F k,l Pt x k,l q,t - Ps x k,l s,q =      1 if P q executes T k and not T l -1 if P q executes T l and not T k 0 otherwise
If tasks T k and T l are not executed on the same processors, we know thanks to Lemma 1 that the f variables define a route from the processor executing T k to the one executing T l .

Thus, any solution of the linear program describes a valid allocation. In particular, there is an allocation with throughput ρ opt .

Multiple paths

Finally, we present our linear programming formulation for the most flexible case, the multiple-paths routing: any transfer may now be split into several routes in order to increase its throughput. The approach is extremely similar to the one used for the single route, free routing policy: we use the same set of f variables to define a flow from processors producing files to processors consuming them. The only difference is that we no longer restrict f to integer values: by using rational variables in [0; 1], we allow each flow to use several concurrent routes. Theorem 4 expresses the optimality of the allocation found by the linear program. Its proof is very similar to the proof of Theorem 3.

( 1c), (1d), (1e), (2f), (1g), (1h), (1i), (2j)

)      MINIMIZE τ UNDER THE CONSTRAINTS (3a) ∀F k,l , ∀P q ❀ P r , x k,l q,r ∈ {0, 1}, y k q ∈ {0, 1}, f k,l i,j ∈ [0; 1] AND (1b), ( 3
The role of the f variable is modified, as expressed by the following lemma.

Lemma 2. Given a file F k,l , the following two properties are equivalent (i) ∀P q Pq→Pr f k,l q,r -

P r ′ →Pq f k,l r ′ ,q =      1 if P q = P prod -1 if P q = P cons 0 otherwise
(ii) the set of links P q → P r such that f k,l q,r = 0 defines a set of weighted routes from P prod to P cons , with total weight 1.

Proof. Consider the subgraph of the platform graph comprising only the linksP q → P r such that f k,l q,r = 0. We construct the set of weighted routes by the following iterative process. We extract a route r from P prod to P cons from this graph (such a route exists thanks to the conservation law). We then compute the minimum weight w of the links in route r. Route r is added with weight w to the set of weighted routes, and the subgraph is pruned as followed: w is subtracted from the value of f k,l of all links included in route r, and links whose f k,l value becomes null are removed from the subgraph. We can prove that Property (i) still holds with value 1 -w instead of 1. We continue the process until there is no more link in the subgraph. Theorem 4. An optimal solution of Linear Program (3) describes an allocation with optimal throughput for the multiple paths policy.

Performance evaluation

In this section, we present the simulations performed to study the performance of our strategies. Simulations allow us to test different heuristics on the very same scenarios, and also to consider far more scenarios than real experiments would. We can even test scenarios that would be quite hard to run real experiments with. This is especially true for the flexible or multiple-path routing policies. Our simulations consist here in computing the throughput obtained by a given heuristic on a given platform, for some application graphs. We also study another metric: the latency of the heuristics, that is the time between the beginning and the end of the processing of one input data set. A large latency may lead to a bad quality of service in the case of an interactive workflow (e.g., in image processing), and to a huge amount of temporary files. This is why we intend to keep the latency low for all input data sets.

Reference heuristics

In order to assess the quality and usefulness of our strategies, we compare them against four classical task-graph scheduling heuristics. Some of these heuristics (Greedy and HEFT) are dynamic strategies: they allocate resources to tasks in the order of their arrival. As this approach is not very practical when scheduling a series of identical workflows, we transform these heuristics into static scheduling strategies: the mapping of a single allocation is INRIA computed with the corresponding strategy (Greedy or HEFT), and then this allocation is used in a pipelined way for all instances.

Greedy. This strategy maps the tasks onto the platform starting from the task with the highest w i,k value, i.e., the task which can reach the worst computation time. The processor that would process this task the fastest is then allocated to this task, and the processor is "reserved" for the time needed for the processing. The allocation proceeds until all tasks are scheduled. Communications are scheduled using either the compulsory route, or the shortest path between allocated processors in case of flexible routing.

HEFT. This heuristic builds up an allocation for a single instance using the classical Heterogeneous Earliest Finish Time [START_REF] Haluk Topcuoglu | Task scheduling algorithms for heterogeneous processors[END_REF]. Then, this allocation is used for all instances.

Pure data-parallelism. We also compare our approach to a pure data-parallelism strategy: in this case, all tasks of a given instance are processed sequentially on a given processor, as detailed in the introduction.

Multi-allocations upper bound. In addition to the previous classical heuristics, we also study the performance when mixing control-and data-parallelism. This approach uses concurrent allocations to reach the optimal throughput of the platform, as is explained in details in [START_REF] Beaumont | Assessing the impact and limits of steady-state scheduling for mixed task and data parallelism on heterogeneous platforms[END_REF]. Rather than using the complex algorithm described in that paper for task graphs with bounded dependences, we use an upper bound on the throughput based on this study, which consists of a simple linear program close to the one described in this paper, and solved over the rational numbers. This bound is tight when the task graph is an in-or out-tree, but may not take all dependences into account otherwise. This upper bound, however, has proved to be a very good comparison basis in the following, and is used as a reference to assess the quality of other heuristics. No latency can be derived from this bound on the throughput, since no real schedule is constructed.

Simulation settings

Platforms. We use several platforms representing existing computing Grids. The descriptions of the platforms were obtained through the SimGrid simulator repository [START_REF]SimGrid: a Generic Framework for Large-Scale Distributed Experimentations. Henri casanova and arnaud legrand and martin quinson[END_REF], as described in [START_REF] Frincu | Handling very large platforms with the new simgrid platform description formalism[END_REF].

• DAS-3, the Dutch Grid infrastructure [START_REF]The distributed asci supercomputer 3 (DAS-3)[END_REF],

• Egee, a large-scale European multi-disciplinary Grid infrastructure, gathering more than 68.000 CPUs [START_REF]Enabling grids for e-science (EGEE)[END_REF],

• Grid5000, a French research Grid which targets 5000, processors [START_REF] Bolze | Grid'5000: a large scale and highly reconfigurable experimental grid testbed[END_REF],

• GridPP, the UK Grid infrastructure for particle physics [START_REF] Britton | GridPP: Meeting the particlephysics computing challenge[END_REF].

Users are often limited to a small number of processors of a given platform. To simulate this behavior, a subset of the available processors is first selected and then used by all heuristics. It is composed of 10 to 20 processors. To evaluate our fixed-routing strategies, we pre-compute a shortest-path route between any pair of processors, which is used as the compulsory route.

Applications. Several workflows are used to assess the quality of our strategies, with a number of tasks between 8 and 12: 1) pipeAlign [START_REF]Pipealign[END_REF], a protein family analysis tool, 2) several random task graphs generated by the TGFF generator [START_REF] Dick | Tgff task graphs for free[END_REF]. In order to evaluate the impact of communications on the quality of the result, we artificially multiply by different factors all communications of the application graphs. There are many ways to define a communication-to-computation ratio (CCR) for a given workflow. We chose to define an average computation time t comp by dividing the sum of all tasks by the average computational power of the platform (excluding powerless nodes like routers), and an average communication time t com by dividing the sum of all files by the average bandwidth. Then the CCR is given by the ratio t com /t comp . Finally, we impose the first task and the last one to be processed on the first processor. These tasks have a size 0 and correspond to the storage of input and output data.

Results

Throughput and Latency

As said before, we compare both throughput and latency for each solution. For each of the 233 application/platform scenarios, results are normalized to the best one: the performance of a heuristic on an instance is divided by the best performance achieved by a heuristic on that instance (or by the upper bound for throughput); therefore, the closer to 1, the better. Figure 4(a) gives the throughput observed for related applications, while Figure 4(b) gives the throughput for unrelated applications. Latencies are displayed for related applications in Figure 5(a) and for unrelated applications in Figure 5(b). The two last figures Figure 6(a) and Figure 6(b) show the latency divided by τ , which is the inverse of the throughput and represents the time between the completion of two successive instances. This gives an idea on the performance for the latency: the smaller the ratio, the shorter the latency in comparison to τ . Normalizing with τ allows us to compare latencies of very different scenarios, in spite of the discrepancies in computing sizes. Note that this ratio also gives the average number of instances that are being processed at any time (started but not completed), and is thus an indication of the storage capacity needed for temporary files.

In the following, our strategies based on Mixed Linear Programing are noted MLP in the following.

Data-parallel.

When dealing with the related model and a low CCR, the data-parallel strategy offers the best throughput since all processors are working at full speed. Since sending input and output data requires communications, performance decreases with the CCR, down to 10% of the upper bound for a CCR greater than 10. When processing DAGs with unrelated computing speeds, this strategy can be the worst even with low CCRs, because other strategies take unrelated speeds into account. Moreover, as it makes use of all processors, even slow ones, the data-parallel heuristic leads to a very high latency, which increases with the CCR. In Figure 6(a) and Figure 6(b), we can see that the ratio between latency and τ (which is also the average number of unfinished instances in the system) is roughly equal to the number of used processors, which is a natural consequence of the data-parallel strategy.

Greedy and HEFT. In both related and unrelated cases, HEFT and greedy heuristics have rather good results when the CCR is very low: around 80% of the upper bound for ratios below 0.1. When the scenario becomes more communication-intensive (CCR equal to, or greater than, 1), their performance is dropping to 30% of the upper bound, and even to 10% when CCR ≥ 10. On the other hand, these strategies build schedules whose latency is often very large compared to those of own MLP strategies. Like the data-parallel strategy, relative latencies increase with the CCR, leading to latencies 100 times worse the best one when the CCR is greater than 10. The ratio between latency and τ is better than for the data-parallel strategy but roughly 3 times higher than our MLP strategies.

MLP strategies. Our three strategies based on linear programs often return similar results: the best one is MLP-multiple-paths, followed by MLP-free-routing, and MLP-fixedrouting. This is quite natural: the more flexible the routing, the better the results. For CCRs below 10, these strategies obtain throughputs between 50% and 80% of the upper bound, and their performance increases with the CCR: when this ratio is over 10, our strategies achieve more than 80% of the upper bound, and often over 90%. Moreover, the MLP strategies always give the best latency. MLP-fixed-routing almost always achieves the best latency: using other routes than the shortest-path, or concurrent routes, can increase the latency (up to a factor 2). Finally, the ratio between latency and τ is very small, often between 2 and 4, less than other strategies.

Running time of the algorithm

Our strategies relies on mixed linear programs, which can take much time to solve. Even if the use of linear programs in our heuristics can significantly slow down the computation of the schedule, it allows to reach better throughputs, and thus better running times. When dealing with very large series of workflows, this gain can be significant and justify the choice of a scheduler. Moreover, specialized tools like GPLK [24] or CPLEX [START_REF]Ilog cplex: High-performance software for mathematical programming and optimization[END_REF] can solve mixed linear programs efficiently. During our simulations, all linear programs and mixed linear programs were solved using CPLEX 11.0 on a 2.4GHz Opteron processor. As we target reasonable numbers of tasks, the processing time of the schedule is kept low: all problems were solved in less than 10 seconds.

INRIA 7 Conclusion and perspectives

In this paper, we have studied the scheduling of a series of workflows on a heterogeneous platform. We have taken advantage of the regularity due to the series to optimize the system throughput by applying steady-state techniques. We have derived single-allocation strategies, which combine good performance with simple control. Indeed, we have proven that our mixed linear program computes an optimal allocation under a number of routing scenarios. Simulations have proven that the benefit of our approach in comparison to classical reference heuristics is significant as soon as communication times are not negligible, and that our allocations lead to much smaller latencies as a side effect. Future work include simplifying our mixed linear program to cope with larger applications, and using task duplication to further improve the system throughput.

Figure 1 :

 1 Figure 1: Examples of platform and applications.

Figure 2 :

 2 Figure 2: Allocation examples for various routing policies.

P 2 → P 3 andP 3 → P 5 P 5 P 5 → P 6 P 6 P 2 P 1 → P 2 P 1

 233555662121

Figure 3 :

 3 Figure 3: Example of a periodic schedule for the allocation represented on Figure 2(a). Only the first instance is represented with task and file labels.

Figure 4 :

 4 Figure 4: Normalized throughputs of different strategies, with related and unrelated applications.

Figure 5 :

 5 Figure 5: Normalized latencies of different strategies, with related and unrelated applications.

Figure 6 :

 6 Figure 6: Ratio between latency and τ for different strategies, with related and unrelated applications.

RR n°6603