Bilel Derbel
email: bilel.derbel@univ-lille1.fr

Arnaud Liefooghe
email: arnaud.liefooghe@univ-lille.fr

Qingfu Zhang
email: qingfu.zhang@cityu.edu.hk

Sébastien Verel
email: verel@univ-littoral.fr

Hernán Aguirre

Kiyoshi Tanaka
email: ktanaka@shinshu-u.ac.jp

A Set-oriented MOEA/D

Keywords: Theory of computation → Evolutionary algorithms, Random search heuristics, • Computing methodologies → Optimization algorithms, Discrete space search, • Applied computing multi-and many-objective optimization, decomposition, evolutionary algorithms, combinatorial optimization Bilel Derbel, Arnaud Liefooghe, Qingfu

The working principles of the well-established multi-objective evolutionary algorithm Moea/d relies on the iterative and cooperative improvement of a number of single-objective sub-problems obtained by decomposition. Besides the definition of sub-problems, selection and replacement are, like in any evolutionary algorithm, the two core elements of Moea/d. We argue that these two components are however loosely coupled with the maintained population. Thereby, we propose to re-design the working principles of Moea/d by adopting a set-oriented perspective, where a many-toone mapping between sub-problems and solutions is considered. Selection is then performed by defining a neighborhood relation among solutions in the population set, depending on the corresponding sub-problem mapping. Replacement is performed following an elitist mechanism allowing the population to have a variable, but bounded, cardinality during the search process. By conducting a comprehensive empirical analysis on a range of combinatorial multi-and many-objective NK-landscapes, we show that the proposed approach leads to significant improvements, especially when dealing with an increasing number of objectives. Our findings indicate that a set-oriented design can constitute a sound alternative for strengthening the practice of multi-and many-objective evolutionary optimization based on decomposition.

INTRODUCTION

Context. A multi-objective problem requires to optimize multiple objectives simultaneously. Since the objectives are many-often conflicting, it is unlikely to find one single solution optimizing all target objectives. One of the most challenging issues is then to find a set of solutions providing different trade-offs between the considered objectives. For this purpose, Evolutionary Multi-objective Optimization (EMO) algorithms have been proved to be extremely efficient in finding high-quality approximation sets [START_REF] Deb | A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF][START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF].

Most existing EMO approches follow a standard algorithmic scheme where a population of solutions is evolved in an iterative manner. Specific mating selection and replacement mechanisms are used to respectively (i) choose some parents from which new individuals are created by means of variation operators, and (ii) update the population by incorporating the newly-generated offspring solutions. Apart from the variation operators, it is well understood that the dynamics of an evolutionary search process in terms of convergence and diversity is tightly related to the design of the selection and replacement mechanisms. Different classes of EMO algorithms can be distinguished according to these two critically-important components. For instance, dominance-based approaches, such as Nsga-ii [START_REF] Deb | A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF] and Pls [START_REF] Paquete | Pareto Local Optimum Sets in the Biobjective Traveling Salesman Problem: An Experimental Study[END_REF], rely on a dominance relation for "comparing" different solutions at the selection and replacement steps. Indicator-based approaches, such as Ibea [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF] and Sms-Moea [START_REF] Beume | SMS-EMOA: Multiobjective selection based on dominated hypervolume[END_REF], rely on a quality indicator to guide the evolutionary search process when deciding the most suitable parents to select and the most suitable solution to maintain in the population. In this paper, we are interested in aggregation-based algorithms, e.g., [START_REF] Giagkiozis | Generalized Decomposition[END_REF][START_REF] Hughes | Multiple Single Objective Pareto Sampling[END_REF][START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF][START_REF] Murata | Cellular Genetic Local Search for Multi-objective Optimization[END_REF][START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF].

Motivations. Our main focus is on the design principles of the so-called Moea/d (Multi-objective evolutionary algorithm based on decomposition) algorithm [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF][START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF], which can be considered as a stateof-the-art framework in the class of aggregation-based approaches. In its basic variant, given a scalarizing function in the objective space, Moea/d considers the cooperative solving of different singleobjective sub-problems, each one defined using the given scalarizing function parametrized by a different weight vector. A population of solutions is then considered, where each sub-problem, equivalently weight vector or search direction, is mapped to one single solution.

The computational flow of Moea/d consists in iterating over the set of weight vectors, by performing selection and replacement in the vicinity of the current weight vector. In fact, Moea/d does not iterate explicitly over the solutions when performing selection and replacement, which might raise some concerns. For example, when the number and distribution of the weight vectors does not properly match the distribution of non-dominated solutions, Moea/d may imply a population of multiple duplicates of the same solution. This might not only impact the (small) cardinality of the output solution set and a lack of diversity in the current population, it also biases the selection of the crucially-important mating selection mechanisms. By handling a set of solutions, where duplicates are not allowed, it can be possible to offer an alternative (re-)design of Moea/d. This idea is the core motivation of our work. In other words, a notable design feature of Moea/d, to be contrasted with the previously-mentioned EMO algorithms, is that the population is not controlled in an explicit manner. We claim that this plays an important role in the behavior of Moea/d, and subsequently in the performance of the underlying evolutionary search process. Besides, we argue that investigating the impact of such a design feature shall allow us to better grasp some of the common issues when choosing Moea/d as an effective solving approach.

Contributions and Results

Overview. In this paper, we adopt a set-oriented perspective allowing us to design an enhanced EMO algorithm based on decomposition. We propose to explicitly evolve a population of solutions considered as a bounded set. We then manage to loosely couple the sub-problems solving with the iterations of the evolutionary process, while still using the underlying (fixed-size) weight vectors as a core element to perform selection and replacement. Based on the scalarized fitness values of solutions, we adopt a many-to-one mapping of the weight vectors to the population, i.e. a subset of weight vectors can be explicitly assigned to one single solution. Mating selection is performed within the neighborhood of solutions, which is defined using the distance in the objective space between the corresponding subsets of weight vectors to which the solutions are mapped. Replacement is performed following a standard elitist mechanism where the population size is not fixed, but bounded by the number of weight vectors.

To validate these design principles, we experiment the proposed algorithm using ρMNK-landscapes, considered as a standard generic combinatorial optimization problem family, with variable problem size, number of objectives, and objective correlation. Through a comprehensive empirical analysis, the proposed algorithm is shown to improve the quality of the obtained approximation sets, and to have a robust behavior in a relatively wide range of problems and configurations. In particular, a substantial improvement is reported as the number of objectives increases and as the maximum size of the population scales. More importantly, our investigations provide evidence that a set-oriented design offers an alternative perspective to the establishment of future high-quality decomposition-based EMO algorithms for both multi-and many-objective optimization.

Outline. The paper is organized as follows. Section 2 presents some necessary background and related work. Section 3 introduces the proposed set-oriented Moea/d framework. Section 4 provides the experimental analysis of our approach. The last section concludes the paper and discusses open issues.

BACKGROUND AND RELATED WORK 2.1 Multi-objective Optimization

A multi-objective optimization problem (MOP) can be defined by a set of M objective functions f = (f 1 , f 2 , . . . , f M ), and a set X of feasible solutions in the decision space. In the combinatorial case, X is a discrete set. Let Z = f (X ) ⊆ I R M be the set of feasible outcome vectors in the objective space. In a maximization context, an objective vector z ∈ Z is dominated by a vector

z ′ ∈ Z (z ≺ z ′ ) iff ∀m ∈ {1, . . . , M }, z m ⩽ z ′ m and ∃m ∈ {1, . . . , M } s.t. z m < z ′ m . A solution x ∈ X is dominated by a solution x ′ ∈ X (x ≺ x ′ ) iff f (x) ≺ f (x ′ ). A solution x ⋆ ∈ X
is Pareto optimal if there does not exist any other solution x ∈ X such that x ⋆ ≺ x. The set of all Pareto optimal solutions is the Pareto set. Its mapping in the objective space is the Pareto front.

The Moea/d Framework

Decomposition-based EMO algorithms [START_REF] Giagkiozis | Generalized Decomposition[END_REF][START_REF] Hughes | Multiple Single Objective Pareto Sampling[END_REF][START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF] seek good-performing solutions in multiple regions of the Pareto front by decomposing the original MOP into a number of scalarized singleobjective sub-problems [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. In this paper, we use the weighted Chebyshev (д te ) scalarizing function, to be minimized:

д te (x, λ) = max i ∈ {1, ..., M } λ i • z ⋆ i -f i (x)
where x ∈ X , λ = (λ 1 , . . . , λ M ) is a positive weight vector such that λ i ⩾ 0 for all i, and [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF][START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF], sub-problems are optimized cooperatively by defining a neighborhood relation between sub-problems. Given a set of µ weight vectors Λ µ = (λ 1 , . . . , λ µ ) defining µ sub-problems, Moea/d maintains a population P µ = (x 1 , . . . , x µ ) where each individual corresponds to one sub-problem. For each sub-problem i ∈ {1, . . . , µ}, a set of neighbors B(i) is defined by considering the T -closest weight vectors using the euclidian distance (denoted || • ||). The sub-problems are optimized iteratively. For each subproblem i, two sub-problems are selected at random from B(i), and the two corresponding solutions are considered as parents. An offspring x ′ is created by means of variation operators (e.g., mutation, crossover). For every j ∈ B(i), if x ′ improves j's current solution x j , then x ′ replaces it, i.e., if д te (x ′ , λ j ) < д te (x j , λ j ) then x j = x ′ . The algorithm loops over the sub-problems, i.e., weight vectors, until a stopping condition is satisfied.

z ⋆ = (z ⋆ 1 , . . . , z ⋆ M ) is a reference point such that, ∀i ∈ {1, . . . , M } , ∀x ∈ X , z ⋆ i > f i (x). In Moea/d

Discussion and Positioning

Rationale. Since mating selection and replacement in Moea/d are performed among neighbors, the neighborhood definition, B(i), is critical for an accurate diversity/convergence balance. For example, replacement can lead to the situation where a solution which improves a sub-problem that does not belong to the neighborhood, is discarded. It can also lead to a situation where several neighbors are assigned the same solution. This has the effect of decreasing the probability that different solutions are picked at the selection step, since the T different neighbors of a sub-problem, in the weight vector space, could imply much fewer than T different solutions in the search space. Such a feature can prevent the reproduction operator from producing high-quality offsprings. This might not only slow down convergence, but it might also make the search process stuck more easily. Additionally, since the Moea/d main iterative process is performed w.r.t. sub-problems, if the solutions maintained by multiple sub-problems are the same, then the corresponding genotype is likely to be considered more often for reproduction. This might not be an issue at a first sight, especially if such a scenario happened because an exceptionally-good solution w.r.t respect to a subset of sub-problems was found at some iteration. However, this might be unfair w.r.t some other sub-problems, because their solutions are then less likely to be considered for reproduction.

Related Works. Different techniques can be adopted to deal with such issues, e.g., adapting the neighborhood size, adjusting the neighborhood w.r.t. the selected sub-problem, restricting the number of replacements, see [START_REF] Trivedi | A Survey of Multiobjective Evolutionary Algorithms Based on Decomposition[END_REF]. With respect to the population management, different variants can also be found, for instance by adopting a global (generational) replacement approach [START_REF] Marquet | Shake Them All! Rethinking Selection and Replacement in MOEA/D[END_REF][START_REF] Wang | Adaptive Replacement Strategies for MOEA/D[END_REF] or a stable matching methodology [START_REF] Li | Stable Matching-Based Selection in Evolutionary Multiobjective Optimization[END_REF]. One variant, denoted Moea/d (δ,nr ) , was extensively considered in the literature [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF]. Two modifications are considered. The first one consists in an extra probability parameter δ allowing parents to be selected from the whole population instead of solely the T -neighborhood. The second one consists in an extra parameter nr bounding the number of neighbors that can be replaced by a newly-generated offspring.

Positioning. To the best of our knowledge, there are no much investigations addressing the impact of iterating over the weight vectors Λ µ in Moea/d, instead of explicitly iterating over the population P µ , considered as a solution set. One can only find few recent works aiming at understanding the behavior of the population dynamics in Moea/d [START_REF] Hugo Monzón | Closed State Model for Understanding the Dynamics of MOEAs[END_REF]. It is in fact worth-noticing that our ultimate goal is not to provide yet another variant of Moea/d, but rather to better grasp the complexity of solving a multi-objective optimization problem by evolving a population based on the concept of decomposition. In this respect, the community still lacks much understanding on the subject despite the success of the Moea/d framework. The population maintained in Moea/d cannot be really considered as a set, in particular because it could include duplicates. Consequently, in contrast to other classes of algorithms, the behavior of Moea/d cannot be easily understood from a set-oriented multi-objective search perspective [START_REF] Zitzler | On Set-Based Multiobjective Optimization[END_REF], which is known to provide a unified theory and best practices of EMO algorithms.

Viewing Moea/d from a set-oriented angle requires to revise its initial design components. In particular, this requires to re-design its main iterative process, as well as the underlying selection and replacement mechanisms, while keeping the spirit of the ingredients that made its wide success, that is: (i) the use of decomposition and (ii) the cooperation among sub-problems. The rest of this paper is devoted to showing how such a perspective can offer new alternatives to enhance the performance of Moea/d while opening a new angle for designing EMO algorithms based on decomposition.

A SET-ORIENTED DESIGN OF MOEA/D

In the following, we adopt a set-oriented perspective to the design of Moea/d. We still consider µ weight vectors defining µ sub-problems, for which we search for a set of good performing solutions. However, we consider the population as a bounded set, which is managed in a more explicit manner. The general idea is to map the weight vectors to the solutions, and to use the mapping information to perform selection and replacement following a standard evolutionary scheme. The flow of the proposed scheme, called Set-Moea/d, is illustrated in Fig. 1, and detailed in Algorithm 1.

Population Mapping. Given the current population P, we compute at every iteration a mapping ϕ : Λ µ g -→ P of the target µ subproblems into the solutions of the population. We adopt an elitist mapping function where every weight vector λ i ∈ Λ µ is mapped to exactly one solution x ∈ P having the best fitness value w.r.t. the scalarizing function д (lines 4 to 6). Since different sub-problems can have the same best solution, ϕ is a many-to-one mapping where a solution x ∈ P can be either: (i) mapped by one or more weight vectors, denoted Φ(x) in line 6, or (ii) non-mapped by any weight vector, i.e., Φ(x) = ∅. The non-mapped solutions are discarded from the population (line 8). After this step, all solutions remaining in the population are different, and duplicates cannot survive.

Neighborhood-based Selection. We consider a standard evolutionary (genetic-like) process where two parents are used for reproduction. The first parent x is selected uniformly at random from the population (line 9). The second parent x is however selected in the neighborhood BS( x) of the first parent (lines 11 and 12). Unlike the conventional Moea/d, the neighborhood is here defined explicitly w.r.t a given solution. The distance between two solutions in the population set P (denoted d in line 10) is defined according to the distance between the weight vectors to which they are respectively mapped. In more details, since ϕ is a many-to-one mapping, a solution x surviving in the population can be mapped to a subset of weight vectors Φ(x). The distance between any two solutions x and x ′ in the population is then defined as the minimum distance between any pair of weight vectors respectively in Φ(x) and Φ(x ′ ). Since the population cannot contain any duplicates, this selection mechanism ensures that the genotypes of parents are different and likely to be close-enough to produce promising individuals.

Reproduction and Population Update. An offspring solution y is produced (line 13) and injected in the population P (line 14). Nevertheless, this does not mean that it will necessarily survive. In fact, as soon as the next iteration starts, the mapping ϕ is computed again, and the population is updated by removing non-mapped solutions. Actually, the mapping enables the replacement in our design. With the newly-generated offspring y, the updated mapping may imply different scenarios. First, y can survive iff it is elected as the best solution for at least one weight vector, i.e., Φ(y) ∅. In such a scenario, if a solution x, to which a weight vector in Φ(y) was mapped in the previous iteration, is still the best for at least one other weight vector, then x survives as well, otherwise it dies. This means that if the offspring y survives, then a number of other solutions might die. In all cases, the population size cannot exceed the pre-defined number of weights µ, but it can drop to less than µ when an offspring improves simultaneously different sub-problems.

Implementation Details and Further Considerations. Notice that the mapping ϕ does not need to be recomputed from scratch at each iteration. Since only the mapping of some weights is likely to change at a given iteration w.r.t. a newly-produced offspring y, one x ← arg min

x ∈P д(x | λ i ); 6 ϕ(λ i ) ← x; Φ(x) ← Φ(x) ∪ λ i ;
// Replacement: remove non-mapped solutions simply needs to compute the aggregated fitness values of the newly created offspring, and then update the mapping and the distances required to define BS in an incremental manner. This can only imply a restricted complexity cost compared to Moea/d. More importantly, the proposed Set-Moea/d algorithm is to be viewed as one specific implementation of a set-oriented design of Moea/d. In fact, the underlying algorithmic framework can be flexibly extended with other components guided both by: (i) decomposition-specific aspects, such as extended notions of neighborhood(s) w.r.t. solutions or other mapping functions, and (ii) by more general plug-andplay EMO techniques, such as steady-state, generational, 'comma' or 'plus' evolutionary strategies. In this respect, decomposition is simply to be considered as a way of efficiently structuring the population, which can then be handled in a more explicit manner.

RESULTS AND EXPERIMENTAL ANALYSIS 4.1 Experimental Setup

ρMNK-Landscapes. we use ρMNK-landscapes considered as a family of problem-independent models for multi-objective multimodal combinatorial optimization [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF][START_REF] Kauffman | The Origins of Order[END_REF][START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. Feasible solutions are as binary strings of size N . For every k ∈ {1, . . . , M }, the k th objective (to be maximized) is defined by:

f k (x) = 1 N N i=1 c k i (x i , x i 1 , • • • , x i K ), such that c k i is the component functions.
Component-values follow a multivariate uniform distribution of dimension M, defined by a correlation matrix [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. We consider the same correlation between all pairs of objectives, given by a coefficient ρ > -1 M -1 . A positive (resp. negative) coefficient decreases (resp. increases) the degree of conflict between the objective values. In our experiments, we set ρ ∈ {-0.2, 0.0, 0.2}, M ∈ {2, 3, 4, 5}, N = 100 and K = 2. Investigating other values for K and N are left for future work.

Parameter Setting. Set-Moea/d, the conventional Moea/d [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF] as well as Moea/d (δ,nr ) [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF], are considered in our experiments. The number of weight vectors µ is set in the range {100, 200, 600} and generated using the method described in [START_REF] Saúl | On the low-discrepancy sequences and their use in MOEA/D for high-dimensional objective spaces[END_REF]. Each considered configuration is executed 20 times. The stopping condition is set to 10 6 evaluation function calls. We use an independent bit-flip mutation with a rate 1 N and one-point crossover. The neighborhood size T is set to 10% of µ. A population of µ randomly-generated individuals is used initially. The reference point z ⋆ is updated online to the best-encountered objective values. Parameters δ and nr are set respectively to 0.1 and 2 in Moea/d (δ,nr ) [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF]. We follow the performance assessment design proposed in [START_REF] Knowles | A Tutorial on the Performance Assessment of Stochastic Multiobjective Optimizers[END_REF] using the hypervolume (to be maximized) and the additive epsilon (to be minimized) indicators [START_REF] Zitzler | Performance Assessment of Multiobjective Optimizers: An Analysis and Review[END_REF]. For each instance, we compute a reference set by merging the approximation sets over all runs. The reference point is set to the worst value in each dimension of the reference set.

Overall Performance

We start our empirical study by providing in Table 1 Furthermore, we can observe a notable change in the relative performance of the two algorithms when looking more carefully at dimension M. Although there is no statistical difference for biobjective problem instances (except for ρ = 0 and µ = 600 where Set-Moea/d is better in terms of epsilon indicator), Set-Moea/d is consistently at least as good as Moea/d, and often significantly better, when the number of objectives increases M ⩾ 3 (except for 2 particular cases, w.r.t the epsilon indicator). Such an observation indicates that Set-Moea/d is a viable alternative to Moea/d for multi-and many-objective optimization. Interestingly, Set-Moea/d is always outperforming Moea/d when using the largest population size (µ = 600), which is a particularly relevant feature for manyobjective optimization. This is confirmed by the following anytime analysis.

Anytime Behavior

When studying the anytime performance of the considered algorithms, the superiority of Set-Moea/d is even more pronounced when the maximum size allowed for the population set is sufficiently large and using a reasonable budget in term of function evaluations. In fact, we depict in Fig. 2 the evolution of the average epsilon indicator value as a function of the number of function evaluation calls, while additionally including in our comparison the extended variant Moea/d (δ,nr ) , for µ = 600. One can clearly see that Set-Moea/d has a better anytime convergence profile than both Moea/d and Moea/d (δ,nr ) for multi-and many-objective problems (M ∈ {3, 4, 5}), and especially for M = 5. Actually, the results are mitigated for bi-objective problem instances (M = 2), especially starting at 10 5 evaluations, but before that Set-Moea/d is the best-performing approach. This is a desirable feature especially when dealing with an expensive optimization scenario when one is only allowed a restricted number of function evaluation calls.

These findings hold independently of the objective correlation ρ. This is a notable feature since one might typically expect that the impact of different ρ-values on the landscape, in terms of solution diversity and Pareto front size, could influence the diversity-toconvergence trade-off of the experimented algorithms. Furthermore, we can extract an interesting observation when analyzing more carefully the anytime behavior of Moea/d (δ,nr ) . At the early stages of the search process, Moea/d (δ,nr ) has the slowest convergence behavior, whereas Set-Moea/d is the most efficient. This can be attributed to the specific replacement design of Moea/d (δ,nr ) which has a tendency to preserve solution diversity. As the search progresses, Moea/d (δ,nr ) is overall able to outperform Moea/d for M ∈ {2, 3} but to a less extent for many-objective problem instances M ∈ {3, 4, 5}, while being always worse than Set-Moea/d.

In Fig. 3, we further illustrate the difference in the anytime behavior of the three algorithms w.r.t. the different values of µ (due to space restrictions, only one many-objective instance is shown). Two insightful observations can be extracted. First, Set-Moea/d appears to be relatively robust to the choice of the cardinality bound µ of the solution set. Using different µ-values leads to different indicator values. With no surprise, larger µ-values lead to better q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -0.2, 2 -0.2, 3 -0.2, 4 -0. q q q q q q q q q q q q q q q q q q 10 -1.2 performance. However, for Moea/d and Moea/d (δ,nr ) , this is only true at the late stages of the search, whereas this tends to hold at anytime for Set-Moea/d (except with very few function evaluations). At the beginning of the run, the cardinality bound of the set maintained by Set-Moea/d has almost no impact on the algorithm performance. Only at the end of the run, having a larger population size allows for more improvement, especially for many-objective problems. This difference implies that the performance of Moea/d and Moea/d (δ,nr ) is quite sensitive to the setting of the population size µ, in the sense that it is correlated to both the available budget and to µ. In contrast, in the case of Set-Moea/d, performance seems to be simply positively correlated to the available budget, which is a natural outcome, especially for many-objective problems. We attribute this feature to the fact that parameter µ plays a fundamentally different role in the design of Set-Moea/d. In fact, parameter µ is to be viewed as one way to structure the population set while imposing a bound on its cardinality. This is to be contrasted to the initial design of Moea/d, where it also influences the intrinsic composition of the population, thus impacting more critically the search behavior, e.g., by introducing duplicates. This is analyzed further in the rest of the paper. M = 2 M = 3 M = 4 M = 5 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 25 50 75

Population Dynamics

As discussed previously, unlike in Set-Moea/d, some duplicated solutions can be maintained by Moea/d. In Fig. 4, we show the evolution of the absolute (average) size of the population maintained by the different considered algorithms. In the case of Moea/d and Moea/d (δ,nr ) , we simply count each single objective vector only once (since there are no equivalent solutions in ρMNK-landscapes).

(Average) Absolute Population Size. The long-term general tendency is roughly the same for all three algorithms. The mean absolute size of the population starts at a small value within relatively few function evaluations (10 3 ), then it grows as the search process evolves. The maximum-allowed size µ can only be reached in the setting where the number of objectives and the computing budget are both relatively high. For a given budget (e.g., 10 5 ), using more weight vectors (e.g., 200) allows the set to attain an even larger size (e.g., > 150) than the maximum allowed when using fewer weight vectors (e.g., 100). Based on our anytime analysis, this does not necessarily implies that the performance, w.r.t. some budget, is better when using more weight vectors because the population is more diverse. This actually renders a more complex behavior that can only be understood by taking into account the mapping of the weight vectors into the Pareto front. In fact, fewer weight vectors, i.e., smaller µ-values, are likely to imply more intensification, which may in turn make the search stuck more easily.

Additionally, we can observe that the absolute size of the population tends to be slightly lower in Set-Moea/d at the early stage of the search process, with the exception of bi-objective instances. The dynamics of Moea/d (δ,nr ) is also different at the beginning, because the number of replacements is restricted to at most nr . However, even if the propagation of duplicates is slower, Moea/d (δ,nr ) does not properly address the issue induced by having duplicates in the population in the long run, and the absolute size of the population drops to the same level than the other algorithms after about 100 generations. As a consequence, we argue that the overall good performance of Set-Moea/d is not only due to the absolute number of solutions (which in fact was shown to be roughly the same). This should make it clear that adopting a set-oriented design allows us to have a better control on the population. For instance, the selection of duplicates can typically have a non-zero probability in Moea/d. By contrast, we are not only able to avoid such a bias, but we can also refine the selection by controlling explicitly which individuals to choose, i.e. through the neighborhood relation BS.

Population Behavior. We can refine our understanding of the population dynamics by analyzing not only its average absolute size, but its actual absolute size in one single particular execution. This is depicted in Fig. 5, corresponding to a typical run of Set-Moea/d. We can see that the absolute population size is indeed increasing overall, but in a non-smooth manner. In fact, the population size decreases periodically in an abrupt manner, and then increases again to reach an even higher value. This can be explained as follows. The algorithm operates in several periods. At each period, newly-found, hopefully diverse, solutions are produced and mapped to different weight vectors. This hence helps the algorithm to jump to an even more efficient region of the objective space, which makes the population size decrease again, an so on. The periodicity of this behavior is however decreasing with time (notice the log-scale in Fig. 5), which can indicate that the algorithm is actually converging. Interestingly, such a dynamic implies a smooth evolution of the population quality, which suggests that the population jump is periodically localized in small regions of the objective space.

CONCLUSION AND PERSPECTIVES

By adopting a a set-oriented perspective, the work conducted in this paper enables to provide an alternative design of multi-objective decomposition-based algorithms. The proposed Set-Moea/d algorithm follows a standard evolutionary search process using a many-to-one mapping between weight vectors and solutions in order to perform selection and replacement. Besides being able to improve the quality of the approximation sets, our experimental analysis provides insights into the behavior of the underlying evolutionary process, while highlighting some interesting features that need to be studied in a more systematic manner in the future. For instance, it would be nice to provide decomposition-based features to better estimate the diversity of the population set, as well as its degree of convergence. This could be helpful to control selection and replacement in an adaptive manner, or to adjust the maximum size of the population online, which we believe to be a challenging but crucially-important perspective. In our design, the replacement is relatively aggressive, since poor quality solutions die immediately. This could be handled differently in the future if diversity appears to be an issue. Moreover, we only consider combinatorial instances in our experiments. It would be interesting to confirm the benefits of such an alternative set-oriented design when tackling other instances from both discrete and continuous domains. Interestingly, the design of the proposed algorithm is flexible enough to enable the integration of other advanced EMO techniques, especially when dealing with many-objective optimization problems, for which a set-oriented perspective seems to be particularly promising.

Figure 1 :

 1 Figure 1: A schematic illustration of Set-Moea/d: many-to-one mapping ϕ, T -neighbors, mating selection, and replacement.

7 for x ∈ P do 8 if 9 x 14 P

 78914 Φ(x) = ∅ then P ← P \ {x } ; // Select one parent ← select a (random) parent in P; // Compute T -Neighboring solutions 10 for x ∈ P do d( x, x) ← min λ ∈Φ( x ), λ ′ ∈Φ(x ) ||λ, λ ′ ||; ; 11 BS( x) ← the T closest solutions in P w.r.t. d; // Select another parent 12x ← select a (random) parent in BS( x); // Reproduction, e.g. crossover and mutation 13 y ← generate an offspring from parents x and x ; // Replacement: merge population with offspring ← P ∪ {y};

Figure 2 :

 2 Figure 2: Anytime performance of Set-Moea/d, Moea/d, and Moea/d (δ,nr ) with a (maximum) population size µ = 600. Subfigures in rows and columns correspond to instances with respectively the different objective correlation ρ ∈ {-0.2, 0.0, 0.2} and the different number of objectives M ∈ {2, 3, 4, 5}. Notice the log-scales. Error bars show 95% confidence intervals.

Figure 3 :

 3 Figure 3: Anytime performance of Set-Moea/d (left), Moea/d (middle), Moea/d (δ,nr ) (right); using different values of (maximum) population size µ for an instance with ρ = 0.0 and M = 5. Error bars show 95% confidence intervals.

Figure 4 :

 4 Figure 4: Evolution of the average absolute size of the population for the different algorithms for different µ ∈ {100, 200, 600} (in rows), and for the instances with M ∈ {2, 3, 4, 5} (in columns) and ρ = 0.0. Error bars show 95% confidence intervals.

Figure 5 :

 5 Figure 5: Dynamics of the population size in a typical run of Set-Moea/d on instances with M ∈ {2, 3, 5} (M = 4 is omitted for clarity) and ρ = 0.0. The left, middle, and right subfigures correspond to a µ-value of 100, 200 and 600, respectively. The Hypervolume (in dashed font) is normalized w.r.t. the max and min values attained by the population size to better visualize any possible correlation. Only the Hypervolume evolution is shown, but not its scale, for more clarity.

Table 1 :

 1 a global overview of the relative performance of Set-Moea/d and Moea/d. Out of the 72 settings (12 ρMNK-landscapes × 3 µ values × 2 indicators), Set-Moea/d is in average better than Moea/d 62 times (i.e. more than 85% of the cases). Set-Moea/d significantly outperforms Moea/d 26 times (36%), whereas Moea/d significantly outperforms Average and standard deviation (in braces) of the indicator-values achieved by Set-Moea/d and Moea/d at 10 6 evaluations. For every (µ, ρ, M), an underlined font indicates that the algorithm in the column is better (w.r.t. considered indicator). When the difference is statistically significant according to a Wilcoxon test (p-value = 0.05), a bold font is additionally used.

			M = 2	M = 3	M = 4	M = 5
			Set-Moea/d	Moea/d	Set-Moea/d	Moea/d	Set-Moea/d	Moea/d	Set-Moea/d	Moea/d
						Hypervolume Indicator × 10 3
		ρ = 0.2	31.088 (1.14)	30.900 (1.02)	5.738 (0.22)	5.678 (0.20)	0.663 (0.05)	0.657 (0.06)	0.090 (0.01)	0.094 (0.01)
	µ = 100	ρ = 0	28.653 (0.72)	28.286 (0.75)	9.863 (0.12)	9.774 (0.18)	2.308 (0.05)	2.295 (0.06)	0.491 (0.02)	0.490 (0.03)
		ρ = -0.2	46.941 (0.95)	46.878 (1.01)	21.539 (0.52)	21.580 (0.32)	7.664 (0.22)	7.687 (0.22)	1.568 (0.04)	1.572 (0.04)
		ρ = 0.2	31.102 (1.02)	30.948 (0.73)	6.117 (0.13)	5.947 (0.19)	0.749 (0.04)	0.730 (0.05)	0.106 (0.01)	0.108 (0.01)
	µ = 200	ρ = 0	28.721 (0.68)	28.574 (0.63)	10.166 (0.13)	10.102 (0.14)	2.488 (0.05)	2.454 (0.06)	0.581 (0.03)	0.563 (0.02)
		ρ = -0.2	47.152 (0.92)	47.101 (0.82)	22.530 (0.38)	22.303 (0.43)	8.429 (0.12)	8.374 (0.13)	1.866 (0.03)	1.817 (0.03)
		ρ = 0.2	31.362 (0.95)	31.157 (0.58)	6.225 (0.11)	6.162 (0.18)	0.818 (0.04)	0.774 (0.06)	0.130 (0.01)	0.117 (0.01)
	µ = 600	ρ = 0	28.814 (0.46)	28.344 (0.79)	10.497 (0.09)	10.380 (0.12)	2.690 (0.04)	2.601 (0.07)	0.675 (0.02)	0.640 (0.02)
		ρ = -0.2	47.459 (0.71)	47.394 (0.77)	23.517 (0.52)	23.049 (0.39)	9.369 (0.15)	9.095 (0.14)	2.241 (0.04)	2.156 (0.05)
						Epsilon Indicator × 10 3
		ρ = 0.2	16.758 (5.09)	17.981 (4.32)	29.163 (5.37)	30.085 (4.14)	44.384 (5.35)	47.798 (6.86)	53.595 (7.73)	53.696 (6.50)
	µ = 100	ρ = 0	13.545 (6.36)	14.262 (5.57)	25.936 (4.08)	27.280 (2.82)	45.575 (3.78)	46.507 (5.25)	70.286 (4.25)	65.958 (5.37)
		ρ = -0.2	17.915 (5.59)	19.255 (7.44)	42.001 (5.82)	42.531 (5.51)	73.402 (3.20)	70.631 (4.04)	94.785 (4.21)	93.876 (3.85)
		ρ = 0.2	18.017 (4.63)	19.872 (4.48)	23.128 (6.85)	28.168 (6.92)	39.180 (5.07)	40.534 (5.59)	49.246 (7.08) 56.998 (10.93)
	µ = 200	ρ = 0	12.191 (4.34)	13.559 (5.70)	21.908 (2.80)	26.148 (3.75)	41.249 (3.68)	41.990 (5.66)	57.324 (4.93)	60.197 (6.31)
		ρ = -0.2	19.636 (5.33)	19.722 (6.01)	38.930 (10.50) 45.049 (11.62)	62.758 (1.66)	60.921 (4.43)	77.930 (3.22)	76.383 (3.40)
		ρ = 0.2	16.206 (8.21)	17.628 (4.39)	23.156 (5.86)	29.094 (7.35)	35.187 (5.73)	42.370 (7.40)	40.739 (8.48) 56.427 (10.79)
	µ = 600	ρ = 0	10.927 (4.19) 15.952 (6.46)	21.105 (3.33)	25.717 (4.29)	35.436 (4.96)	44.992 (4.75)	46.554 (4.59)	60.742 (8.33)
		ρ = -0.2	16.656 (5.03)	18.086 (5.81)	36.026 (13.32) 48.277 (11.54)	50.600 (4.82)	55.374 (5.47)	61.547 (4.31)	69.937 (7.06)

Set-Moea/d only 2 times (3%), and there is no statistically significant difference 44 times (60%). Actually, a similar overall superiority of Set-Moea/d is found when compared to Moea/d (δ,nr ) , but this is not reported in Table

1

due to space restrictions. Notice that Set-Moea/d is never statistically outperformed w.r.t hypervolume.

Acknowledgments. This work was supported by the French national research agency (ANR-16-CE23-0013-01) and the Research Grants Council of the Hong Kong Special Administrative Region, China (RGC Project No. A-CityU101/16).


Céline Merlet
email: celine.merlet@upmc.fr

Mathieu Salanne
email: mathieu.salanne@upmc.fr

Benjamin Rotenberg

New Coarse-Grained Models of Imidazolium Ionic Liquids for Bulk and Interfacial Molecular Simulations

Keywords: molecular dynamics, surface tension, graphite, supercapacitor, electrode, force field

We introduce new coarse-grained models for two imidazolium based ionic liquids, namely 1-butyl-3-methyl-imidazolium tetrafluoroborate [BMI][BF 4 ] and 1-ethyl-3-methylimidazolium tetrafluoroborate [EMI][BF 4 ], derived from the original force field of Roy and Maroncelli (J. Phys. Chem. B 2010, 114, 12629-12631) representing the 1-butyl-3-methylimidazolium hexafluorophosphate

[BMI][PF 6 ] ionic liquid. We evaluate static and dynamic properties between 298 K and 500 K and show that they agree with previous experimental and all-atom simulations studies. The models are used to conduct simulations of the liquid-vapor interface and accurately predict surface tensions at 400 K. Capacitive properties are also examined by doing molecular dynamics simulations of the ionic liquids in contact with graphite electrodes. The obtained structures and capacitances are consistent with all-atom simulations results reported on these systems.

I. INTRODUCTION

Molecular dynamics simulations appear as a key method in many fields to understand experimental observations and to predict properties which are difficult to obtain experimentally. Due to the broad range of potential applications envisioned for room temperature ionic liquids (RTILs), numerous force fields have been developed to model these liquids. Among them, several polarizable 1,2,3 and non polarizable 4,5,6,7 potentials are sufficiently accurate to predict bulk properties and to analyse the behavior of the liquid under confinement 8,9,10 .

One major potential application for ionic liquids is their use for supercapacitors 11,12,13 .

In these systems, energy is stored at the electrode / electrolyte interface through reversible ion adsorption. The electrodes usually consist in nanoporous carbon electrodes which give relatively high energy densities. Molecular simulations in this area will help understanding the microscopic phenomena at play and designing new materials and ionic liquids to optimize energy storage. To conduct realistic molecular dynamics simulations (MD) on these complex systems, it is necessary to include an accurate description of: i) the ionic liquid, ii) the ionic liquid / electrode interactions, iii) the electrode including the microscopic structure and polarizability of the carbon atoms. In most of the recent studies, the ionic liquid and the ionic liquid / electrode interactions are represented in an appropriate way but simple model geometries are used for the electrode structures leading to underestimated capacitances in comparison with experiments involving realistic supercapacitors 14,15 . To recover the correct capacitances and gain insight in the microscopic structure adopted by the liquid in the nanopores, one should include the complexity and polarizability of the carbon structure 16 .

Because of the high computational cost of considering the polarizability of the electrode and the large number of carbon atoms in the calculation, it is necessary to reduce the time spent on the calculation of the interaction forces. In a previous study, we showed that the recently developed coarse-grained model representing the [BMI][PF 6 ] ionic liquid 17 is performant to predict the structural and capacitive properties of the liquid-graphite interface 18 .

According to these considerations, an electrochemical cell was designed using this simple model for the electrolyte and carbon electrodes generated using quenched molecular dynamics 19 leading to the conclusion that this setup is appropriate for the study of supercapacitors systems 16 . The coarse-grained potentials that have been investigated to date perform well but all include bonded and nonbonded interaction terms 20,21,22,23,24,25 . The capabilities of the rigid-body potential of [BMI][PF 6 ] 17,18 encouraged us to look for such simple models.

In this work, we report the development of new coarse-grained models for two imidazolium based ionic liquids, 1-butyl-3-methylimidazolium tetrafluoroborate ([BMI][BF 4 ]) and 1-ethyl-3-methylimidazolium tetrafluoroborate ([EMI][BF 4 ]), and show that they are able to predict accurately a number of properties. In section II, we describe the models developed and the methodology adopted. In sections III and IV, we deal with the structural and dynamic bulk properties of the ionic liquids. In section V, we give an account of the liquid-vapor interface properties predicted. In section VI, we provide the description of the liquid-graphite interface obtained with these models. 17 . To ensure that the new and original parameters are as compatible as possible, we kept the changes to a minimum. In the representation adopted here, site-site interaction energies are given by the sum of a Lennard-Jones potential and coulombic interactions:

II. MODEL AND SIMULATION METHODS

The model parameters for [BMI][BF

u ij (r ij ) = 4ε ij σ ij r ij 12 - σ ij r ij 6 + q i q j 4πε 0 r ij , (1)
where r ij is the distance between sites and ε 0 is the vacuum permittivity. A schematic view of the correspondence between atoms and interaction sites is given in figure 1. The parameters are varied independently to match the density and diffusion coefficients of [BMI][BF 4 ] and

[EMI][BF 4 ] at 400 K. The geometry of EMI + was chosen from the centers of mass positions of the different regions of the cation in typical configurations of an all-atom model (imidazolium ring, ethyl chain, methyl chain). Charges for the EMI + cation were calculated summing the atomic charges of an all-atom non polarizable force field 5 and scaling them to recover the reduced charge of 0.78 e on the ion used in the original model. The conservation of this reduced charge for all the ions is required to allow for future works on mixtures of the ionic liquids. All the interaction parameters of the coarse-grained models are summarized in table I.

We made use of three types of simulations to calculate different properties. First, bulk simulations were employed to determine densities, radial distribution functions, diffusion given by the sum of a Lennard-Jones potential and coulombic interactions

Interaction x y z M σ i ε i q i in BMI + q i in EMI + site ( Å) ( Å) ( Å) (g.mol -1 ) ( Å) (kJ.mol -1 ) (e) (e)
u ij (r ij ) = 4ε ij [( σ ij r ij ) 12 - ( σ ij r ij ) 6 ] + q i q j 4πε 0 r ij
where r ij is the distance between sites, ε 0 is the permittivity of free space and cross-interaction parameters are calculated by Lorentz-Berthelot mixing rules. The parameters for the carbon atoms of the graphite electrodes are σ C = 3.37 Å and ε C = 0.23 kJ.mol -126 . coefficients and viscosities. Second, simulations of the liquid-vacuum interface were used to estimate the surface tension. Finally, more complex setups were introduced to test the structural and capacitive properties of the liquid-electrode interface. For all the simulations, the timestep for the integration of the rigid body equations of motion is 2 fs. For the simulations conducted in the NVT ensemble, a weak Nosé-Hoover thermostat 27 with a time constant of 10 ps is employed. Electrostatic interactions were treated using a 3D or 2D Ewald summation according to the periodic boundary conditions adopted for each simulation type.

All the bulk properties were computed by simulations on a cubic cell containing 432 (for

[BMI][BF 4 ]) or 576 (for [EMI][BF 4 ]
) ion pairs. Cubic periodic boundary conditions were employed. When the temperature is modified, the system is allowed to equilibrate for 400 ps in the NPT ensemble using time constants of 10 ps and 20 ps for the thermostat and barostat respectively 28 . At the end of this relaxation time, the equilibrated volume of the cell is used to evaluate the density of the system. The statistics for the other properties are then typically collected during 1 ns in the NVT ensemble. Longer simulations were used to evaluate the viscosity for temperatures under 373 K.

Interfacial properties were calculated from simulations on orthorhombic cells with lengths such that L x ∼ L y < L z , according to the fact that the interface is oriented perpendicularly to the z direction. For the determination of the surface tension, a slab of ionic liquid, consisting of 1280 ions pairs, is sandwiched between two vacuum regions. The entire simulation cell has the size 7.2 nm × 7.7 nm × 50 nm. These values were chosen as a result of our previous work on [BMI][PF 6 ] 18 . The simulation is started from a configuration near equilibrium, already at 400 K, and allowed to relax for another 400 ps in the NVT ensemble. The statistics are then gathered for 1 ns using the same conditions.

For liquid-graphite simulations, the system consists of 320 ion pairs sandwiched between two groups of three graphite layers, a snapshot of the simulation cell is given in supporting information. The lengths in the x and y directions are chosen to match with the lattice parameters of graphite and are equal to 3.22 nm and 3.44 nm respectively. The length in the z direction is adjusted to recover the bulk density of the ionic liquid at 400 K leading to L z = 11.26 nm for [BMI][BF 4 ] and L z = 9.68 nm for [EMI][BF 4 ]. The method used to enforce a constant potential difference between the walls, as opposed to a constant charge, following the methodology of Siepmann and Sprik 29 , is detailed in our previous work 18 . This representation of the electrodes polarizability has been successfully implemented in simulations of various systems 18,30,31,32,33 . The potential difference (∆Ψ 0 ) is applied between the positive (Ψ + ) and negative (Ψ -) electrodes in order to have:

Ψ + = -Ψ -= ∆Ψ 0 2 .
When the potential difference is modified, the system is equilibrated during 200 ps before gathering data during 1 ns. The potential difference is increased gradually from 0.0 V to 1.50 V to ensure that the system is effectively equilibrated before gathering statistics. The simulations are conducted at 400 K in the NVT ensemble applying 2D periodic boundary conditions.

III. BULK STATIC PROPERTIES

Figure 2 compares densities evaluated under a pressure of 1 bar for our models with experimental 6,34,35,36 and all-atom simulations 5 data. As expected, the density is larger for

[EMI][BF 4 ] than for [BMI][BF 4
] for all the studied temperatures. The agreement between experimental data and ours is fairly good, in particular, the variation with temperature is well reproduced by the coarse-grained models. For this property, our new models perform even better than the all-atom ones of Chaban and coworkers 5 . 6,34,35,36 and all-atom simulations 5 data. Our coarse-grained models give results in agreement with experiments and all-atom models. In particular, the variation with temperature is well reproduced.

Then, we evaluate the radial distribution functions between anions and cations in order to assess the ability of our models to predict the local structure adopted by the liquids. The radial distribution functions, plotted in figure 3, can be compared with the ones obtained by Andrade et al. 37 (with symbols) are all-atoms models results 6,37 . Some peaks characteristic of the all-atom models merge to give the resulting peaks presented here for the coarse-grained potentials.

and [EMI][BF 4

] respectively. Thus, it appears that our models slightly underestimate the anion-imidazolium ring distance. Figure 3 also shows that the anion-cation peaks are wider and with a lower maximum intensity for all-atom models compared with our models.

For [BMI][BF 4 ], there are two others striking features in going from the all-atom model to the coarse-grained model. First, for the C1-C1 curve, the two peaks at 0.44 nm and 0.77 nm in the all-atom model merge into a single flat peak centered around 0.61 nm, i.e.

the average between these two values. Second, for the anion-anion curve, in a similar way, the peaks present at 0.64 nm and 0.80 nm merge to give a single peak at 0.72 nm, but the resulting peak is broader. Comparable phenomena appear for the [EMI][BF 4 ] case.

According to the previous observations, the densities and local structure of the bulk systems seem to be correctly reproduced by the coarse-grained models developed in this work, keeping in mind that some structural features are averaged out as a result of the coarse-grained process.

IV. BULK DYNAMIC PROPERTIES

We now turn to the calculation of some dynamic bulk properties of the ionic liquids and their confrontation with results from experiments and previous models. Diffusion coefficients were calculated using the Einstein relation:

D α = lim t→∞ 1 6t |δr i (t)| 2 , (2)
where δr i (t) is the displacement of a given ion of species α in time t and the brackets denote an ensemble average. The diffusion coefficients obtained for different temperatures between 298 K and 500 K are given in figure 4. The coarse-grained models tend to underestimate the diffusion at low temperatures at which the ionic liquids are highly viscous. For EMI-BF 4 , the increasing gap between values at high temperatures could be explained by the fact that the experimental values are extrapolated from data gathered between 250 K and 380 K.

Nevertheless, the agreement is good over the range of temperatures studied.

The viscosity of the systems was determined following the Green-Kubo relation:

η = V k B T ∞ 0 σ αβ (0)σ αβ (τ ) dτ, (3)
where σ αβ is one of the components of the stress tensor and V is the simulation cell volume.

As this property is a collective one, its precise evaluation is more difficult compared to diffusion coefficients.

The correlation function of the stress tensor was taken to be the mean of its five different anisotropic components. The viscosities resulting from our simulations are summarized in table II along with experimental 34,36,38 and all-atom simulations 5 data. The models are able to predict viscosities in agreement with experimental data within the range of experimental errors except at 348 K where this property is overestimated.

At this point, we can say that our coarse-grained models are sufficiently accurate to predict static and dynamic properties of the bulk systems for the range of temperatures studied. ments 34,36,38 and all-atom simulations 5 . Results for anions and cations are given in green (full line) and red (dashed line) respectively.

V. SURFACE TENSION OF THE LIQUID-VACUUM INTERFACE

The estimation of the surface tension is a good test of the performance of a model for the prediction of interfacial properties. The surface tension is determined following a standard methodology 39,[START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Bresme | [END_REF] , i.e. using simulations of a slab of ionic liquid sandwiched between two vacuum regions. The system is simulated using 3D periodic boundary conditions and Ewald summations to compute the electrostatic interactions. We demonstrated that for [BMI][PF 6 ],

a length of 50 nm (including the liquid slab) in the direction perpendicular to the interface was sufficient to avoid interactions between the slabs 18 . As the ionic liquids studied here present similar characteristics, we assume that this length is appropriate for our present study. The surface tensions are derived from the stress tensor components 39,[START_REF] Frenkel | Understanding Molecular Simulation[END_REF] :

γ = L z 2 σ zz - σ xx + σ yy 2 . ( 4
)
We obtain a value of 33. 44.3 mN.m -1 at 400 K, is also in agreement with the experiments giving 41.9 mN.m -1 at 400 K 43 .

VI. STRUCTURAL AND CAPACITIVE PROPERTIES OF THE LIQUID-GRAPHITE SYSTEM

The coarse-grained models built here have proven their accuracy for a number of bulk and interfacial properties but, because the goal of this work is to develop models which can be used in simulations of supercapacitors, we need to check that the liquid / carbon interface behaves realistically. A simulation cell consisting of two graphite electrodes surrounding a slab of ionic liquid was used to test the structural and capacitive properties of the systems.

From our simulations results, we plot the ionic densities for the centers of mass of anions and cations when various potential differences are applied between the electrodes. Ionic density profiles corresponding to ∆Ψ 0 = Ψ + -Ψ -= 0.0 V and ∆Ψ 0 = 1.5 V are given in supporting information for the two liquids studied here. If we assume that applying a zero potential difference or a zero constant charge affects the ionic liquid in the same way, then we can compare the ionic density profiles for [EMI][BF 4 ], with ∆Ψ 0 = 0.0 V, with the published results of Shim et al. 44 . Indeed, the ionic densities generated with our model of [EMI][BF 4 ] (see figure 5) present the same features as the ones from the all-atom simulations of Shim and coworkers 44 . A layering is observed at null potential difference, induced by the molecular sizes of the ions. ) For both coarse-grained and all-atom models, a first peak is observed for cations and anions at the same distance from the interface with ionic densities equal to approximately 20 ions per nm 3 . A small peak is observed for cations at a distance around 0.25 nm from the first peak. The following peaks of the layered structure are separated by 0.40 nm, a distance for which both models also agree. The comparison between the ionic density profiles when a non zero potential difference (or a non zero charge on the electrodes) is applied cannot be made because the fixed charge in their simulations is more than twice the highest fluctuating charge in our work. Nevertheless, we also observe an enhancement of the counter-ions ionic densities near the electrodes and the beginning of the alternating of anions and cations layers.

ρ BF 4 - ρ EMI + Graphite ΔΨ 0 = 0.0 V
The ionic density profiles for [BMI][BF 4 ] are very similar to the ones of [EMI][BF 4 ] but the differences between anions and cations peaks are more pronounced for [BMI][BF 4 ], prob-ably because of the more important asymmetry between the ions geometries for this ionic liquid. When a potential difference of 1.5 V is applied, the alternating layers are more visible

for [BMI][BF 4 ] compared to [EMI][BF 4 ]
. These small distinctions lead to differences in the capacitive properties of the ionic liquids which we now underline.

The capacitance of the ionic-liquid / graphite systems can be determined using the Poisson potential and the total charge on the electrodes. The plots showing the surface charge as a function of the potential drop for the ionic liquids studied in this work are given in figure 6. The differential capacitances are extracted from these curves using:

C diff = ∂σ ∂∆Ψ , ( 5
)
where σ is the total charge on the electrode and ∆Ψ is the potential drop across the interface We note that the capacitances for the positive electrodes in the two tetrafluoroborate ionic liquids are very similar, in concordance with the fact that the ionic liquids share the same anion, but do differ a little. We notice that the potential of zero charge (PZC) is The ionic liquids / graphite interfaces that we investigate here have also been studied

(∆Ψ ± = Ψ ± -Ψ bulk ).
BMI-BF 4 C -= 5.5 µF.cm -2 C + = 3.9 µF.cm -2 EMI-BF 4 C -= 4.9 µF.cm -2 C + = 4.2 µF.cm -2 BMI-PF 6 C -= 4.8 µF.cm -2 C + = 3.9 µF.cm -2 44 give C -= 3.9 µF.cm -2 and C + = 11.0 µF.cm -2 for a constant charge of ± 6.9 µC.cm -2 . We recall here that apart from the use of all-atom models to represent the ionic liquids, the main differences between our simulations and theirs are: i) the use of constant charges instead of constant potential differences and ii) the estimation of the capacitances as integral capacitances as opposed to differential capacitances (i.e.

the capacitances are evaluated using one point charge only). As a consequence of the large charges used in the all-atom simulations and the calculation of the capacitance as an integral one, the estimated values may be out of the linear régime. Which, in addition to the fact that we employ a coarse-grained model, can be at the origin of the quantitative discrepancy between the results.

The coarse-grained models developed in this work generate the same capacitive behavior as for the [BMI][PF 6 ] ionic liquid 18 and the simulations at constant potential differences allow us to gain insight on the influence of the ions on the electrode charging. The distributions of the carbon charges for ∆Ψ 0 = 0.0 V and ∆Ψ 0 = 1.5 V are given in figure 7 for the two ionic liquids studied here. Plots obtained for [BMI][PF 6 ] are added for comparison.

For the zero potential difference, we observe that the curve is not exactly centered around zero and that, while the shapes are very similar for (values coincide for relatively high charges).

A general feature, for the three ionic liquids, is the asymmetry between the positive and negative electrodes, with in particular, broader distributions on the positive side. This statement is consistent with the different geometries of anions and cations in our systems.

The variations in the charge distributions show that the use of constant charges uniformly distributed, even if sufficient to recover the capacitances of a planar electrode / ionic liquid interface 18 , leads to an important loss of information on the systems under examination.

VII. CONCLUSION

The aim of this work was to develop coarse-grained models for two commonly studied ionic liquids, namely [BMI][BF 4 ] and [EMI][BF 4 ], in order to simulate complex systems such as supercapacitors including explicitly the structure of the carbon electrodes. We show that the parameter sets we established lead to the prediction of a number of properties in good agreement with experimental and all-atom simulations data.

First, we calculated static bulk properties, the densities and radial distribution functions of the ionic liquids, demonstrating that the models perform well for these basic quantities.

Then, we verified that the two models are also able to accurately predict dynamic bulk properties. Next, we used simulations of the ionic liquid / vacuum interfaces to estimate the surface tensions of [BMI][BF 4 ] and [EMI][BF 4 ] at 400 K which are consistent with experimental results 42,43 . Finally, we tested the behavior of the ionic liquids confined between graphite layers under various applied potential differences. The qualitative behavior of the systems is similar to the ones observed for equivalent setups 18,44,45 . These new coarsegrained models have thus proven their accuracy for a number of relevant properties which make them suitable for simulations of complex systems, such as supercapacitors, for which a high computational efficiency is needed on the ionic liquid side.

VIII. ASSOCIATED CONTENT

Supporting Information available: Densities, diffusion coefficients and viscosities calculated in this work, a snapshot of the liquid / graphite simulation cell and additional ionic density profiles. This material is available free of charge via the Internet at http://pubs.acs.org.

IX. ACKNOWLEDGMENTS

Figure 1 :

 1 Figure 1: Schematic representations of the coarse-grained models of the three ions considered: BMI + , EMI + and BF - 4 .

Figure 2 :

 2 Figure 2: Comparison of densities obtained in this work with experimental 6,34,35,36 and all-atom

4 Figure 3 :

 43 Figure 3: Radial distribution functions at 400 K calculated in this work for anion-anion, anioncation and cation-cation pairs. Black lines (without symbols) are results from this work, red lines

4 Figure 4 :

 44 Figure 4: Diffusion coefficients obtained for our new models and comparison with experi-

  8 mN.m -1 at 400 K for [BMI][BF 4 ] which compares well with the experimental value of 38.2 mN.m -1 at 393 K 42 . The surface tension for [EMI][BF 4 ],

Figure 5 :

 5 Figure 5: Ionic density profiles of the centers of mass of the anions and cations for [EMI][BF 4 ] when a zero potential difference is applied between the electrodes. The curves present the same features as the all-atom simulations of Shim et al. 44 at constant zero charge.

Figure 6 :

 6 Figure 6: Surface charge as a function of the potential drop across the interface for [BMI][BF 4 ] and [EMI][BF 4 ] investigated in this work. Results for [BMI][PF 6 ] are given for comparison 18 . The differential capacitances of the positive and negative interfaces are given by the slopes obtained in the linear régimes.

  different for the two systems, around -69 mV for [BMI][BF 4 ] and 29 mV for [EMI][BF 4 ].

Figure 7 :

 7 Figure 7: Distributions of the carbon charges on the positive and negative electrodes for three different ionic liquids and for 0.0 V and 1.5 V potential differences applied. The distributions are influenced by the nature of the ions.

  4 ] and [EMI][BF 4 ] were derived from the force field recently developed by Roy and Maroncelli for [BMI][PF 6 ]

Table II :

 II Shear viscosities calculated in this work compared to experimental data34,36,38 and allatom molecular dynamics simulations 5 .

	η (mPa.s)	348 K 373 K 400 K 500 K
	BMI-BF 4
	This work	25.4 9.7	5.1	1.7
	Experimental 36 15.1 8.6	5.4	1.9
	All-atom MD 5 13.0 8.0	5.0	-
	EMI-BF 4
	This work	-	7.5	5.3	1.8
	Experimental 34 8.8	5.8	4.0	1.7
	Experimental 38 12.4 8.0	5.4	2.1
	All-atom MD 5 10.0 6.0	4.0	-

  As for [BMI][PF 6 ] 18 , we identify two linear régimes meaning that the differential capacitances are constant over the range of potential differences studied. The capacitances are summarized in table III along with the results for [BMI][PF 6 ] 18 .

		4
	) -2 Total charge (µC.cm	-2 0 2	BMI-BF 4 EMI-BF 4 BMI-PF 6
		-4	-0.5	0 ΔΨ	+/-(V)	0.5	1

Table III :

 III Differential capacitances obtained in this work and comparison with published values for [BMI][PF 6 ] 18 . by molecular dynamics simulations with all-atom models 44,45 . For [BMI][BF 4 ], Feng et al. 45 found capacitances around 6.7 µF.cm -2 for a constant charge of ± 10 µC.cm -2 . While for [EMI][BF 4 ], Shim et al.

  [BMI][PF 6 ] and [BMI][BF 4 ], the results are different for [EMI][BF 4 ]. This suggests that the nature of the cation influences more the charge distribution than the nature of the anion at zero potential. We notice that, for ∆Ψ 0 = 1.5 V, the charges on the negative electrode are shifted toward negative values for [EMI][BF 4 ] compared to the other ionic liquids and that [BMI][BF 4 ] is intermediate between the other ionic liquids. On the positive electrode side, [BMI][BF 4 ] also appears as an

		8		ΔΨ	0 = 0.0 V
	Fraction of C atoms (%)	2 4 6					BMI-BF 4 EMI-BF 4 BMI-PF 6
		0	-0.02	-0.01		0	0.01	0.02
				Charge on atom (e)
	12		ΔΨ	0 = 1.5 V
	Fraction of C atoms (%)	3 6 9	Negative electrode			BMI-BF 4 EMI-BF 4 BMI-PF 6 Positive electrode
		0	-0.02	-0.01		0	0.01	0.02
				Charge on atom (e)

intermediate between [BMI][PF 6 ] (values coincide for relatively low charges) and [EMI][BF 4 ]

13, 14723-14734.

We acknowledge the support of the French Agence Nationale de la Recherche (ANR) under grant ANR-2010-BLAN-0933-02 ("Modeling the Ion Adsorption in Carbon Micropores"). We thank Paul Madden for fruitful discussions.