Fodil Benali
email: fbenali@adwanted.comdbodenes@adwanted.com

Damien Bodénès

Cyril De Runz
email: cyril.derunz@univ-tours.fr

Nicolas Labroche
email: nicolas.labroche@univ-tours.fr

An Enhanced Adaptive Geometry Evolutionary Algorithm Using Stochastic Diversity Mechanism

Keywords: Theory of computation → Evolutionary algorithms K-Means++ initialization, Non-Euclidean-Geometry, Evolutionary Computation

Evolutionary Algorithms have been regularly used for solving multi and many objectives optimization problems. The effectiveness of such methods is determined generally by their ability to generate a well-distributed front (diversity) that is as close as possible to the optimal Pareto front (proximity). Analysis of current multiobjective evolutionary frameworks shows that they are still suboptimal and present poor versatility on different geometries and dimensionalities. For that, in this paper, we present AGE-MOEA++, a new Multi and Many Objective Evolutionary Algorithm that:

(1) incorporates the principle of Pareto Front (PF) shape fitting to enhance the convergence in different shaped high dimensional objective spaces, and (2) adapts K-means ++ fundamentals in order to best manage the diversity in non-uniform distributed PF. The empirical study shows that our proposal has better results than the state-of-the-art approaches in terms of IGD and is competitive in terms of GD.

INTRODUCTION

A multi-objective optimization (MOP) problem with M objectives to be minimized can be written as follows [START_REF] Kalyanmoy | Multi-objective optimisation using evolutionary algorithms: an introduction[END_REF]:

𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑓 ( ì 𝑥) = [𝑓 1 ( ì 𝑥), 𝑓 2 ( ì 𝑥), . . . , 𝑓 M ( ì 𝑥)] 𝑇 , (1)
𝑆𝑢𝑏 𝑗𝑒𝑐𝑡 𝑡𝑜 ì 𝑥 ∈ U ⊂ R 𝑛 (2)
where 𝑓 𝑖 (𝑥) is the 𝑖 𝑡ℎ objective function (i ∈ 1; M ), ì 𝑥 is an ndimensional decision vector, and U is the feasible region of ì 𝑥. The objective function vector 𝑓 maps a solution ì 𝑥 in an n-dimensional decision space to a point ([𝑓 1 ( ì 𝑥), 𝑓 2 ( ì 𝑥), . . . , 𝑓 M ( ì 𝑥)] 𝑇 ) in an Mdimensional objective space. A multi-objective problem with M > 3 is often referred to as many-objective problem.

In general, a MOP does not have one single optimal solution since objectives are conflicting with each other. Therefore, no single solution simultaneously optimizes all objectives. Thus, solutions are compared using the Pareto dominance principle. A solution ì 𝑥 dominates another solution ì 𝑦 (denoted by ì 𝑥 ≻ ì 𝑦), if 𝑓 𝑖 ( ì 𝑥) ≤ 𝑓 𝑖 ( ì 𝑦) ∀𝑖 ∈ 1; M and there exists at least one index 𝑗 ∈ 1; M such that 𝑓 𝑗 ( ì 𝑥) < 𝑓 𝑗 ( ì 𝑦). If a solution ì 𝑥 * is not dominated by any other solutions, ì 𝑥 * is said to be Pareto optimal solution. The set of all Pareto optimal solutions is the Pareto optimal solution set (PS). The projection of PS in the objective space is called the Pareto front (PF) or Pareto frontier.

A large number of evolutionary algorithms (EAs) have been proposed in the literature [START_REF] Bader | HypE: An algorithm for fast hypervolume-based many-objective optimization[END_REF][START_REF] Beume | SMS-EMOA: Multiobjective selection based on dominated hypervolume[END_REF][START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints[END_REF][START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF][START_REF] He | Evolutionary many-objective optimization based on dynamical decomposition[END_REF][START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF][START_REF] Sun | IGD Indicator-Based Evolutionary Algorithm for Many-Objective Optimization Problems[END_REF][START_REF] Xiang | A many-objective evolutionary algorithm with Pareto-adaptive reference points[END_REF][START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF][START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach[END_REF] to solve Multi and Many Objective optimization problems. The goal of each algorithm is to approximate the optimal Pareto front with a set of nondominated solutions that are as close as possible to PF (proximity), and that is well-distributed over the optimal PF (diversity).

Multi Objective Evolutionary Algorithms (MOEAs) are often categorized into three classes based on their fitness evaluation mechanisms [START_REF] Ishibuchi | Many-objective problems are not always difficult for Pareto dominance-based evolutionary algorithms[END_REF]. First, decomposition-based algorithms such as MOEA/D [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF], NSGA-III [START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints[END_REF] and [START_REF] He | Evolutionary many-objective optimization based on dynamical decomposition[END_REF][START_REF] Xiang | A many-objective evolutionary algorithm with Pareto-adaptive reference points[END_REF], where a MOP is decomposed into a number of single-objective problems. Each single-objective problem has the same scalarizing function and a different weight vector. A single solution is assigned to each single-objective problem. All single-objective problems are optimized in a cooperative manner towards different directions in the objective space along the weight vectors. However such methods rely on multiple parameters (e.g., setup of the reference lines) which is not adequate for a decision maker in an industrial context.

Second, Pareto dominance-based algorithms such as SPEA [START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach[END_REF] and NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF], where the Pareto dominance relation among solutions is used as the primary fitness evaluation criterion and a crowding distance to promote diversity. However, it was reported in many studies [START_REF] Deb | Towards estimating nadir objective vector using evolutionary approaches[END_REF][START_REF] Ishibuchi | Effectiveness of scalability improvement attempts on the performance of NSGA-II for many-objective problems[END_REF][START_REF] Sato | Controlling dominance area of solutions and its impact on the performance of MOEAs[END_REF] that the performance of Pareto dominance-based severely degrades with the increase in the number of objectives. This is because almost all solutions become non-dominated with each other in very early generations. As consequence, solutions are not comparable with each other, the same fitness value is assigned to all solutions by the Pareto dominance criterion. Thus, no strong selection pressure is generated to push the population towards the PF. Finally, Indicator-based algorithms, such as HV-based methods [START_REF] Bader | HypE: An algorithm for fast hypervolume-based many-objective optimization[END_REF], work well on many-objective problems. However, they have a large computation load for calculating the HV contribution of each solution.

Another problem with a number of MOEAs mentioned above is that they are built upon the implicit assumption that the PF (or its approximation) has an Euclidean geometry. For example, in NSGA-III the reference points are generated using Das and Dennis's systematic approach [START_REF] Das | Normal-boundary intersection: A new method for generating the Pareto surface in nonlinear multicriteria optimization problems[END_REF], which places points on a flat hypersurface. However, in many MOPs, the PF is convex (i.e., hyperbolic geometry) or concave (i.e., spherical geometry). As consequence, this impacts greatly the convergence and diversity of the obtained solutions.

To overcome such limitations, a new MOEA has been proposed recently, called AGE-MOEA [START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF] (Adaptive Geometry Estimation based MOEA), for evolutionary multi and many-objective optimization. The idea consists in fitting the geometry of the PF in order to adapt diversity and convergence mechanisms on the Pareto Front shape. AGE-MOA is the best in terms of convergence and diversity of solutions on the Pareto Front. However, it has a very time consuming mechanism to promote diversity (a high computational complexity O(M𝑁 2 + 𝑁 3 )). This is due to the use of deterministic-based diversity promoting mechanism that penalizes the performance of the algorithm greatly.

The problem of data selection in a non-uniform distribution in order to guarantee certain diversity has been widely discussed in several data mining tasks. In particular, in clustering methods like kmeans, where it has been demonstrated that the quality of clustering is extremely sensitive to initial centroids selection. To remedy to this problem, there has been an extremely efficient and simple algorithm for selecting data in a non-uniform distribution called k-mean++ [START_REF] Arthur | k-means++: The advantages of careful seeding[END_REF]. Augmenting k-means with this stochastic approach (k-means++) has shown that the obtained initial set of centers is provably competitive with the optimal solution. Consequently, the obtained performance is O(log k)-competitive with the optimal clustering.

For that, in this work, we propose AGE-MOEA++, a new Multi and Many Objective Evolutionary Algorithm that: (1) incorporates the principle of PF shape Fitting to enhance the convergence in different shaped high dimensional objective spaces, and (2) adapts K-means ++ fundamentals in order to best manage the diversity in non-uniform distributed PF. Our contributions are the following: (i) we formalize the problem of diversity maximization in MOP, (ii) we propose a simple, parameter-less, Multi and Many Objective Evolutionary Algorithm with a fast and simple stochastic diversity mechanism (𝑂 (M𝑁 2 ) [START_REF] Bachem | Approximate k-means++ in sublinear time[END_REF]), and (iii) we describe experiments that show to which extent the proposal is able to balance between convergence and diversity in different shaped and dimensional test problems.

The remainder of the paper is organized as follows: Section 2 describes L 𝑝 norms and their relation with curvatures and shapes of the hyper-surfaces. Section 3 details the proposed framework AGE-MOEA++ while Section 4 describes the empirical study and results. Finally, Section 5 concludes and opens future work perspectives.

PRELIMINARIES: L 𝑝 NORMS AND NON-EUCLIDEAN GEOMETRY

A norm on a real linear space X is a mapping ∥.∥ from X into R + that satisfy the following axioms: (i) ∥ ì 𝑥 ∥ ≥ 0 with equality if and only if 𝑥 = 0, (ii) ∥𝛼 ì

𝑥 ∥ = |𝛼 |.∥ ì 𝑥 ∥ (𝛼 ∈ R), (iii) ∥ ì 𝑥 + ì 𝑦 ∥ ≤ ∥ ì 𝑥 ∥ + ∥ ì 𝑦 ∥
In an M-dimensional Euclidean space R 𝑀 , the norm (the length) of a vector ì 𝜐 = (𝜐 1 , . . . , 𝜐 𝑀 ) is defined by:

∥ì 𝜐 ∥ 2 = (𝜐 2 1 + 𝜐 2 2 + • • • + 𝜐 2 𝑀 ) 1/2 (3)
In the Euclidean geometry, the distance between two points 𝐴 and 𝐵 is the length (norm) of the straight line segment connecting the two points, i.e., 𝑑 (𝐴, 𝐵) = ∥𝐴 -𝐵∥ 2 . However, it is a well known fact that the higher the dimensionality of a space, the more sparse the space is. Here, one of the concerns is that the Euclidean distance loses its ability to depict spacing in high dimensions. As a consequence, the Euclidean norm does not necessarily provide the most accurate measure of the distance between two points in a high dimensional space [START_REF] Anthony | Minkowski geometry[END_REF]. A generalization of the Euclidean norm is the L 𝑝 norm, given by:

∥ì 𝜐 ∥ 𝑝 = (𝜐 𝑝 1 + 𝜐 𝑝 2 + • • • + 𝜐 𝑝 𝑀 ) 1/𝑝
(4) In fact, the set of points equidistant from a reference point (e.g., the origin of the axes) varies depending on the used L 𝑝 norm. The set of all points with a distance ∥.∥ 𝑝 = 1 to the origin of the axes forms a unit hyper-surface. The curvature of the hypersurface strictly depends on the value of the exponent 𝑝. To better understand this aspect, let us consider a bi-dimensional space. For 𝑝 = 1, the unit curve is flat and corresponds to the straight line connecting the point (0,1) and [START_REF] Akagi | Exact algorithms for the max-min dispersion problem[END_REF]0). Therefore, all points that are equidistant to the origin of the axes lie on that straight line. For 𝑝 > 1, the unit curve is concave and when 𝑝 = 2 it corresponds to the unit circle, i.e., a circle with a radius equal to one. Instead, for 𝑝 < 1, the equidistant points lie on a hyperbolic curve. Therefore, the value of p determines the curvature of the unit hyper-surface associated with the norm.

In this work, we propose to adapt the norm used to evaluate the distance between solutions and the Pareto front to converge more efficiently towards optimal solutions.

PROPOSED ALGORITHM

With the aim of developing a tool that should be able to (i) converge to solutions optimizing a large number of objectives simultaneously, (ii) assist the Decision Maker (DM) in choosing the most suitable configurations by providing a diverse set of solutions, and (iii) be fast in order to obtain efficient solutions within a limited computational time budget.

The proposed framework, called AGE-MOEA++, inherits from the computational fast and elitist framework NSGA-II, as illustrated in Algorithm 1. It incorporates K-means++ initialization heuristic in order to serve the DM with a well distributed and diverse set of solutions and uses the AGE-MOEA Pareto Front geometry estimation technique to be adapted to different shapes of the Pareto Front. The objective is to have an efficient management for both of the diversity and convergence towards the Pareto Optimal Front. A general description of the proposed framework is presented in Section 3.1. As outlined in Algorithm 1, the framework starts with an initial set of N solutions (line 2 in Algorithm 1). This initial population of solutions could be randomly generated or empty solutions. Then an exploration mechanism is executed based on proper crossover and mutation operators, depending on the problem, to produce new offsprings (line 4). The offspring population Q is therefore combined with the current population P forming a new population P ∪ Q of size 2 × 𝑁 (line 5). After that, an evaluation process is executed to evaluate for each solution candidate the set of M objective functions (line 6). Once the evaluation process is terminated, a survival selection process is done to ensure that unfit solutions are eliminated from the population and reduce it back to N individuals (line 7). The steps 4-7 in Algorithm 1, are repeated until a stop condition (e.g., number of generations) is satisfied.

Overview

The selection process, as shown in Algorithm 2, starts by dividing the population into different levels of non-dominated fronts using the fast-non-dominated sorting (NDS) algorithm [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] (line 2 Algorithm 2). Then, the non-dominated fronts are normalized (line 3) using the normalization procedure defined in Section 3.2. Once the solutions are normalized, the first front F 1 is used to estimate the L 𝑝 norm that best fits the geometry of approximated PF (line 4) (More details are described in Section 3.3).

The new population for the next generation is created in lines 5-14 of Algorithm 2. As long as the size of the new population is not reached. Non-dominated Fronts' individuals are assigned a survival score based on the level of the non-dominated front and the estimated 𝑝 norm. Survival score computation procedure is defined in Algorithm 3.

Solutions from the best non-dominated levels are chosen frontwise (lines [START_REF] Jay | Practical nonparametric statistics[END_REF][START_REF] Das | Normal-boundary intersection: A new method for generating the Pareto surface in nonlinear multicriteria optimization problems[END_REF][START_REF] Kalyanmoy | Multi-objective optimisation using evolutionary algorithms: an introduction[END_REF][START_REF] Deb | Simulated binary crossover for continuous search space[END_REF]. The goal of this step is to keep good performance solutions (the most advanced solutions). This elitist strategy will allow to converge faster towards the Pareto-optimal Front. However, it is commonly the case that the last front could not be entirely maintained to fit the size (N ) of the next generation population. As a consequence, the remaining solutions are chosen according to the descending order of their survival scores (lines [START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints[END_REF][START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF].

We have to mention here that the survival score is computed for the first non-dominated-fronts even though they are chosen front wise. The idea behind is to use survival scores during the reproduction process. In fact, parents are selected for reproduction using the binary tournament selection mechanism: a pair of individuals is randomly selected from the population; the winner of the tournament is the solution with the best non-dominated level or the solution with the largest survival score and the smallest constraints violation at the same level of non-dominated rank.

Details about the principle components of AGE-MOEA++ are described in Sections 3.2, 3.3, and 3.4.

Normalization

Normalizing the objective space is an important step in multiobjective optimization problems. In fact, objective functions may have different scales, which leads to neglecting one or more objective functions. As a consequence, this impacts how diversity and performance measures of solutions are compared when dominance relations are not sufficient. For that, we suggest normalizing our objective functions by applying the same formula used in [START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints[END_REF] and also used in AGE-MOEA:

f 𝑖 (𝑆) = 𝑓 𝑖 (𝑆) -𝐼𝑑𝑒𝑎𝑙 𝑖 𝑎 𝑖 , ∀𝑆 ∈ F, ∀𝑖 ∈ 1; 𝑀 (5)
𝐼𝑑𝑒𝑎𝑙 𝑖 = min ∀𝑆 ∈F 1 𝑓 𝑖 (𝑆), ∀𝑖 ∈ 1; 𝑀 (6)
[𝑧 1 , . . . , 𝑧 𝑖 , . . . , 𝑧 M ] 𝑇 × [1/𝑎 1 , . . . , 1/𝑎 𝑖 , . . . , 1/𝑎 M ] 𝑇 = 1 M (7)

𝑧 𝑖 = 𝑎𝑟𝑔 max 𝑆 ∈F 1 (𝑓 𝑖 (𝑆) -𝐼𝑑𝑒𝑎𝑙 𝑖 ), ∀𝑖 ∈ 1; 𝑀 (8)
Where 𝑓 𝑖 (𝑆) denotes the objective 𝑓 𝑖 for the solution 𝑆 and 𝐼𝑑𝑒𝑎𝑙 𝑖 is 𝑖 𝑡ℎ component of the ideal point. The ideal point represents the minimum value of the 𝑖 𝑡ℎ objective across all solutions in the first front F 1 . The objectives are translated to have the ideal point equals to the origin of the axes. Thereafter, a M-dimensional linear hyperplane 𝑍 𝑚𝑎𝑥 is constructed based on the extreme points (𝑧 𝑖 , 𝑖 ∈ 1; 𝑀 ) in each objective axis of the first Front F 1 . The denominator 𝑎 𝑖 is the intercept of the M-dimensional hyperplane with the objective axis 𝑓 𝑖 , and it is obtained by solving the linear system in equation [START_REF] Coello | A study of the parallelization of a coevolutionary multi-objective evolutionary algorithm[END_REF]. In case that the system is indefinite or leading to abnormal normalization (e.g., 𝑎 𝑖 = 0), the min-max normalization is used (i.e., 𝑎 𝑖 is replaced by 𝑧 𝑖 𝑖 -𝐼𝑑𝑒𝑎𝑙 𝑖 ).

Geometry Fitting

To determine the norm L 𝑝 such that the corresponding unit hypersurface best fits the geometry of normalized objectives, we need to find the value of 𝑝 that makes all points in F 1 equally distant to the ideal point, which coincides with the origin of the axes ì 0 after the normalization.

The fitting process consists in solving the following of non-linear equations:

               (𝑓 1 (𝑆 1 ) 𝑝 + 𝑓 2 (𝑆 1 ) 𝑝 + • • • + 𝑓 𝑀 (𝑆 1 ) 𝑝 ) 1 𝑝 = 1 (𝑓 1 (𝑆 2 ) 𝑝 + 𝑓 2 (𝑆 2 ) 𝑝 + • • • + 𝑓 𝑀 (𝑆 2 ) 𝑝 ) 1 𝑝 = 1 . . . (𝑓 1 (𝑆 𝑞 ) 𝑝 + 𝑓 2 (𝑆 𝑞 ) 𝑝 + • • • + 𝑓 𝑀 (𝑆 𝑞 ) 𝑝 ) 1 𝑝 = 1 (9)
where 𝑞 is the number of points in the front F 1 . Several numerical analysis methods have been proposed to resolve such systems of nonlinear equations (e.g., Newton's iterative method [START_REF] Jr | Numerical methods for unconstrained optimization and nonlinear equations[END_REF], Levenberg-Marquardt algorithm). However, such methods are computationally expensive and not suitable for computing the value of 𝑝 in each iteration with negligible overhead.

For that reason, we approximate the value of 𝑝 using the method proposed in [START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF]. It consists in using the central point of the front F 1 for which the corresponding L 𝑝 exponential equation can be easily computed with an exact method. Moreover, the overall complexity of this method is 𝑂 (𝑀 × 𝑁 ), which makes it a fast procedure to estimate the geometry of the non-dominated front and can be easily incorporated in the cycle of evolutionary algorithms.

Survival Score

To select the best individuals that will participate in the next generation, a survival score based on proximity and diversity is calculated as the following:

𝑆𝑐𝑜𝑟𝑒 (𝑆) 𝑆 ∈F 𝑟 = 𝑃𝑟𝑜𝑥𝑖𝑚𝑖𝑡𝑦 (𝑆) × 𝐷𝑖𝑣𝑒𝑟𝑠𝑖𝑡𝑦 (𝑆), 𝑖 𝑓 𝑟 = 1 𝑃𝑟𝑜𝑥𝑖𝑚𝑖𝑡𝑦 (𝑆)
Otherwise ( 10)

𝑃𝑟𝑜𝑥𝑖𝑚𝑖𝑡𝑦 (𝑆) = 1 ∥ f (𝑆)∥ 𝑝 (11)
𝐷𝑖𝑣𝑒𝑟𝑠𝑖𝑡𝑦 (𝑆) = min 𝑆,𝑇 ∈F 1 ,𝑆≠𝑇 ∥ f (𝑆) -f (𝑇 )∥ 𝑝 (12)
In fact, for each generic solution 𝑆 ∈ F 𝑟 , the survival score is defined based on the rank of its non-dominated front as illustrated in equation [START_REF] Kalyanmoy | Multi-objective optimisation using evolutionary algorithms: an introduction[END_REF]. For solutions of the first non-dominated front F 1 , the survival score is defined by the combination of the proximity and the diversity scores. While for solutions of the other fronts, the proximity score only is used. The idea behind is to enhance the convergence to the optimal PF by increasing the chance of survival of the most advanced solutions in the rest of the fronts.

For the proximity measure defined in (Eq. ( 11)), it is evaluated as the inverse of the distance between the solution and the ideal point (see Eq. ( 6)). The diversity of the solutions 𝑆 ∈ F 𝑟 is computed as the p-norm distance to the nearest adjacent solutions in the front F 1 . Algorithm 3 details the procedure that assigns survival scores. The survival scores for the first non-dominated front F 1 are computed in lines 2-21, while the score for the other non-dominated fronts F 𝑟 (𝑟 > 1) are computed in lines 21-26. For what concerns F 1 , first, all extreme points are assigned the maximum possible survival score (+∞) with the purpose to preserve them in the next generation's population (lines 3-4). This allows to maintain the curvature of the PF during the process of convergence. Then, two sets are initialized: (i) Ω containing all solutions yet to score (line 5), and (ii) Ω keeps track of already scored solutions (line 6). After that, the proximity score for the solutions in Ω is computed according to (Eq.11) (lines 7-9), and the pairwise L 𝑝 distances between all solutions in F 1 are computed in lines 10-14. The survival score is computed within the loop in lines 16-20. In each loop, the method computes the diversity score for a random solution selected based on a stochastic diversity mechanism inspired by K-means++ initialization heuristic (line 17). More details about the method are described in Algorithm 4.

Finally, Survival scores for the solutions in (F 𝑟 , 𝑟 > 1) are calculated as proximity scores (lines 27-30). Hence, solutions closer to the hypersurface generated by L 𝑝 have larger scores.

Diversity: Problem Formalization and Stochastic Resolution

The problem of diversity promoting in MOP and MaOP, in which the objective is to have a well distributed and a good dispersion of solutions over the obtained PF, can be seen as an instance of the NP-Hard Max-Min dispersion maximization problem [START_REF] Erkut | The discrete p-dispersion problem[END_REF]. In fact, Max-Min dispersion maximization problem considers a set 𝑃 of 𝑛 elements, and a function 𝑑 that assigns a non-negative real number 𝑑 (𝑝, 𝑞) for each pair of elements 𝑝, 𝑞 ∈ 𝑃. The objective is to find a subset 𝐺 ⊆ 𝑃 such that the cost(G) = 𝑚𝑖𝑛{ 𝑑 (𝑝, 𝑞)| 𝑝, 𝑞 ∈ 𝐺 } is maximized [START_REF] Akagi | Exact algorithms for the max-min dispersion problem[END_REF]. As consequence, our sub-problem of diversity promoting can be formulated as finding the subset Ω ⊆ Ω such that:

Ω = arg max 𝐺 ⊆Ω 𝐷𝑖𝑣𝑒𝑟𝑠𝑖𝑡𝑦 (𝐺) (13)
with:

𝐷𝑖𝑣𝑒𝑟𝑠𝑖𝑡𝑦 (𝐺) = 𝑚𝑖𝑛 ||𝑆 𝑖 , 𝑆 𝑗 || 𝑝 ∀ 𝑆 𝑖 , 𝑆 𝑗 ∈ 𝐺 (14)
To deal with such problems, we suggest to use the fast and simple k-means++ initialization heuristic. The idea consists in building the subset Ω incrementally by choosing random individuals with very specific probabilities. Basically, each solution is selected with probability proportional to its contribution to the overall cost function. Noticeably, we have to mention here that in the specific context of MOEA, we need to score all the individuals of the population. As a consequence, the problem could be reformulated as finding the best selection order that maximizes the overall cost function.

For what concerns diversity score evaluation, Algorithm 4 details the mechanism used to assign diversity score for each individual (solution) of the initial population Ω. In fact, for the first individual, a uniform probability (lines 2-6 of Algorithm 4) is used to choose a point at random from Ω (line 7). After that, the procedure defines its diversity score as the minimum L 𝑝 distance with regards to the solutions in Ω (line 9). Ω is the set of scored solutions seted initially with extreme points. For the next individuals' selection, we choose a solution based on a probability proportional to its contribution to the overall diversity by calculating the ratio between the diversity gain for each solution and the sum of the gains (lines 11-13). Each time, the diversity of a solution 𝑆 is computed with regards to solutions that have already been selected (Ω) in the previous iterations rather than considering all solutions in F 1 . As consequence, this allows to well characterize the diversity of solutions.

EMPIRICAL STUDY

This section reports the conducted experiments to show the effectiveness and validity of AGE-MOEA++ to solve multi and many objective optimization problems based on different literature test problems. In order to evaluate our approach, we answer the following questions:

(1) How does AGE-MOEA++ balance between convergence and diversity in multiobjective optimization context?

To investigate this, we propose experiments where stateof-the-art methods NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF], NSGA-III [START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints[END_REF], MOEA/D [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF], AGE-MOEA [START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF] are compared to AGE-MOEA++ using literature based multiobjective test problems. (2) How does our approach scale in terms of the number of objectives? To find out, we conduct experiments on high dimensional objective spaces to assess the scalability of our approach.

Experimentation Environment Settings

Table 1: Experiments shared parameters. n denotes number of decision variables.

Parameters M = 3 M = 5 M = 10
Population Size [START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints[END_REF] 91 210 275

Number of fitness evaluation [START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF] 27300 63000 82500 PF Evaluation (IGD, GD) Size 2016 4845 5005

Number of Generations 300

Polynomial Mutation Probability [START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF] 𝑝 𝑚 = 1/𝑛 Mutation Distributed Index [START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF] 𝜂 𝑚 = 20 SBX Probability [START_REF] Deb | Simulated binary crossover for continuous search space[END_REF] 𝑝 𝑐 = 1 SBX Distributed Index [START_REF] Deb | Simulated binary crossover for continuous search space[END_REF] 𝜂 𝑐 = 30

We implemented AGE-MOEA++1 in Python using Pymoo [START_REF] Blank | pymoo: Multi-objective optimization in python[END_REF]. Pymoo provides the source code for all benchmark problems as well as the algorithms we use as baselines in our study. For all MOEAs, we used the same parameter setting reported in the related literature [START_REF] Deb | An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints[END_REF][START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF][START_REF] Panichella | An adaptive evolutionary algorithm based on noneuclidean geometry for many-objective optimization[END_REF][START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF]. Table 1 shows all parameters' values for the evaluated MOEAs. For all the other parameters we use the values suggested by their developers. To ensure a fair comparison, we use the same population size, the same number of fitness evaluations, and the same size of optimal PF for all algorithms in our study. In particular, we set the population size N=91, 210, and 275 for the number of objectives M=3, 5, and 10 respectively. The number of fitness evaluations is set to 𝑁 × 300 iterations.

To assess the effectiveness of AGE-MOEA++ in Multi and Many-Objective Problems (MaOPs) context, we considered DTLZ [START_REF] Deb | Scalable multi-objective optimization test problems[END_REF] test benchmarks, with the number of objectives M= 3, 5, and 10. It contains several test problems with different properties, such as concave (e.g., DTLZ2), convex (e.g., Convex DTLZ2), disconnected (e.g., DTLZ7), and degenerate [START_REF] Ishibuchi | Pareto Fronts of Many-Objective Degenerate Test Problems[END_REF] (e.g., DTLZ5) PFs. Therefore, such a suite is a good representation of various real-world scenarios. Fig. 1 illustrates the curvature of some DTLZ test problems used in our experiments.

For each test problem, we run each algorithm 30 times to account for their non deterministic nature. In each independent run, we collected the solutions produced by a given algorithm at the end of the search and computed:

• The inverted generational distance (IGD) [START_REF] Coello | A study of the parallelization of a coevolutionary multi-objective evolutionary algorithm[END_REF] to measure its overall quality. In fact, the IGD measures the distance from any point in optimal PF to the closest point in obtained PF. As consequence, it provides a single scalar value measuring both proximity and diversity of the obtained solutions. • The GD performance indicator [START_REF] Allen | Multiobjective evolutionary algorithms: classifications, analyses, and new innovations[END_REF] to evaluate the convergence of the algorithms because it measures the distance from obtained solutions to the optimal PF. The smaller the IGD and GD values, the better the performance of the algorithm. To evaluate the significance of the differences among the different MOEAs, we use the Wilcoxon rank-sum test [START_REF] Jay | Practical nonparametric statistics[END_REF] with significance level 𝛼 = 0.05. A significant p-value ( p-value < 𝛼) indicates that an algorithm A achieves significant performance in terms of indicator measures than another algorithm B across 30 runs for a given test problem.

Empirical Results

In this subsection, we present and discuss the results of the comparison between AGE-MOEA++ and the four MOEAs baselines. Table 2 and 3 summarize the number of benchmark problems for which AGE-MOEA++ significantly outperforms a baseline MOEA and vice versa based on the IGD and GD measures respectively. Table 4 provides the average (mean) and the standard deviation (shown between parentheses) of the IGD achieved by each MOEA across 30 independent runs. The table also reports whether a baseline (e.g., MOEA/D) statistically outperforms (denoted with ↑), or is statistically equivalent (≈) to AGE-MOEA++ according to the Wilcoxon rank-sum test. Each time, the best performance is highlighted in blue color. For what regards NSGA-II, we can notice that AGE-MOEA ++ outperforms NSGA-II in 13 out of 18 problems based on the IGD and 14 out of 18 problems based on the GD measure. For some problems, the differences between the two MOEAs are above two orders of magnitude. For example, for DTLZ2 test problem and M = 10 (Many Objective Problems), AGE-MOEA ++ obtains an average IGD value equals to 6.156e(-1) compared to an IGD value of 1.25e(+1) for NSGA-II. On the other side, NSGA-II performs significantly better than AGE-MOEA ++ in one test problem: DTLZ 4 for M = 3. For the remaining 3 test problems, there is no significant difference. This confirms the result which shows that NSGA-II is not suitable for many objective optimization problems.

From Tables 23, AGE-MOEA++ achieves significantly lower (better) values than NSGA-III in 14 out of 18 problems based on IGD indicator and 12 out of 18 based on the GD measure. In particular, NSGA-III produces globally significantly higher (worst) IGD values for the majority of test problems. The most significant difference is observed for DTLZ 4 with M = 10 where AGE-MOEA++ obtains an average IGD value equals to 4.913e(-1) while NSGA-III achieves an IGD value of 1.71e(+0). This is owing to the fact that the method used to generate reference lines is NSGA-III is not adaptable to the shape of the PF.

Finally, from the comparison between AGE-MOEA++ and MOEA/D, we can observe from Tables 2-3, the former significantly outperforms the latter in the majority of test problems (13/18 based on IGD and 11/18 based on the GD). The largest difference between the two MOEAs is observed for the concave DTLZ 2 test problem with M = 3. In this case, AGE-MOEA++ has an average IGD value equals to 8.472e(-2) while MOEA/D has an IGD value of 7.397e(-1). In particular, for M = 10, MOEA/D outperforms AGE-MOEA++ based on DTLZ1, DTLZ2, and DTLZ3 test problems. This is due to the shape of these Pareto fronts (Continuous linear and concave PF) which makes it advantageous through the use of reference lines.

CONCLUSION

In this paper, we propose a new MOEA, called AGE-MOEA++, that inherits from the computational fast and elitist framework NSGA-II, incorporates K-means++ initialization heuristic to obtain distributed and diverse set of solutions while reducing the diversity process complexity and uses the AGE-MOEA Pareto Front geometry estimation technique to be adapted to different shapes of the PF. Our solution shows better results than the state-of-the-art approaches in terms of IGD and is competitive in terms of GD. As future works, we aim to work on the integration of DM's preferences in terms of different subsets of the PF in order to build an efficient optimization/recommendation system for different combinatorial problems such as Advertising Campaigns Allocation Problem.

Table 3: Number of DTLZ benchmarks [START_REF] Deb | Scalable multi-objective optimization test problems[END_REF] in which an algorithm A (e.g., AGE-MOEA++) statistically outperforms (<) another algorithm B (e.g., MOEA/D) according to the Wilcoxon test (p-value ≤ 0.05) Based On the GD measure.

Algorithm 1 : 6 F

 16 AGE-MOEA++ Framework. Input: Number of Objectives M, Size of the Population N; Output: Final population P; ←-Evaluate(P);// with M objectives 7 P ←-Select-Survival(P, F, N);

Algorithm 2 : 2 F 3 F 5 P 9 P 14 P

 2235914 AGE-MOEA++ Selection Mecanism. Input: Population P, Target population size N; Output: Next generation population P 𝑛𝑒𝑤 ; 1 begin ←-Fast-Non-Dominated-Sort(P); ←-Normalize(F); 4 𝑝 ←-Fit-P-Norm(F 1 ); 𝑛𝑒𝑤 ←-∅; 6 Calculate-Survival-Score(F 1 , 𝑝, 1); 7 𝑟 ←-1; 8 while | P 𝑛𝑒𝑤 | + |F 𝑟 | ≤ N do 𝑛𝑒𝑤 ←-P 𝑛𝑒𝑤 ∪ F 𝑟 ; 10 𝑟 ←-𝑟 + 1; 11 Calculate-Survival-Score(F 𝑟 , 𝑝, 𝑟 ); 𝑛𝑒𝑤 ←-P 𝑛𝑒𝑤 ∪ F 𝑟 [1 : ( N -| P 𝑛𝑒𝑤 |) ]; 15 return P 𝑛𝑒𝑤 ; 16 end

Algorithm 3 : 2 if r == 1 then 3 Γ 5 Ω 6 Ω 12 D[𝑆 1 , 16 while

 3235612116 AGE-MOEA++ Survival Score Procedure. Input: Pool of Non-Dominated-Front F 𝑟 , Exponent 𝑝 of estimated 𝑛𝑜𝑟𝑚 , 𝑟 rank of the non-dominated-front; 1 begin ←-Extreme-Points(F 𝑟 ); 4 Survival-Score[Γ] ←-+∞; ←-F 𝑟 \Γ; 𝑆 2 ] ←-| |𝑓 (𝑆 1 ) -𝑓 (𝑆 2 ) | | 𝑝 ; |Ω | > 0 do 17 W, 𝑆 * , Value, Ω, Ω ←-Stochastic-Diversity(W, Ω, Ω, D); 18 Diversity[𝑆 * ] ←-Value; 19 Survival-Score[𝑆 * ] ←-Diversity[𝑆 * ] × Proximity[𝑆 * ];

Algorithm 4 :begin 2 if 4 W 7 𝑆 8 Ω ←-Ω\𝑆 * ; 9 10 Ω

 42478910 AGE-MOEA++ Stochastic Diversity Mechanism. Input: Points Sampling Probabilities 𝑊 , Set of initial solutions Ω, Set of selected solutions Ω, L 𝑝 Pairwise distance 𝐷; Output: Updated sampling probabilities weights W, Selected solution 𝑆 * , Diversity value, Updated set of initial solutions Ω, Updated set of selected solutions Ω; 1 𝑊 = ∅ then // Case Weights Not Initialized 3 foreach 𝑆 ∈ Ω do // Uniform Probability Weights * ←-Random-Solution(Ω, W[Ω]); Diversity[𝑆 * ] ←min D[𝑆 * , Ω]; ←-Ω ∪ {𝑆 * }; 11 foreach 𝑆 ∈ Ω do 12 W[𝑆] ←-min 𝑉 ∈Ω 𝐷 [𝑉 ,𝑆 ] 𝑇 ∈Ω min 𝑉 ∈Ω 𝐷 [𝑉 ,𝑇 ] ; 13 end 14 return 𝑊 , 𝑆 * , Diversity, Ω, Ω;

Figure 1 :

 1 Figure 1: Curvature of Used DTLZ test problems.

Comparaison M = 3 M = 5 M = 10 AGE

 310

Comparaison M = 3 M = 5 M = 10 AGE

 310

Table 2 :

 2 Number of DTLZ benchmarks[START_REF] Deb | Scalable multi-objective optimization test problems[END_REF] in which an algorithm A (e.g., AGE-MOEA++) statistically outperforms (<) another algorithm B (e.g., MOEA/D) according to the Wilcoxon test (p-value ≤ 0.05) Based On the IGD measure.

Table 4 :

 4 IGD values achieved by our algorithms and baselines on DTLZ benchmark with M =3, 5, and 10 objectives.

	Problem	M MOEA/D	NSGA-III	AGE-MOEA	NSGA-II	AGE-MOEA++
	DTLZ1	3	7.3784e+2 (5.05e+1) 2.484e+2(2.57e+1)2.634e+2(3.47e+1)	3.05e+2(3.35e+1)	2.15e+2(2.53e+1)
	DTLZ2	3	7.397e-1(7.15e-2)	1.148e-2(1.23e-1) 9.477e-2(5.5e-3)	1.036e-1(5.3e-3)	8.472e-2(4.41e-3)
	C-DTLZ2 3	1.755e-1(3.14e-2)	1.109e-1(1.96e-2) 7.26e-2(7.64e-3)	1.106e-1(1.04e-2)	6.26e-2(4.44e-3)
	DTLZ3	3	2.26e+3(8.04e+1)	7.027e+2(5.30e+1)≈6.51e+2(8.11e+1) 7.80e+2(9.11e+1)	6.87e+2(5.36e+1)
	DTLZ4	3	6.046e-1(4.67e-2)	3.28e-1(2.81e-1) 5.37e-1(3.21e-1)	↑1.08e-1(7.08e-3)	4.88e-1(3.03e-1)
	DTLZ5	3	5.985e-1(7.23e-2)	3.668e-2(4.39e-3) 4.202e-2(3.24e-3)	≈2.999e-2(3.884e-3) 3.439e-2(4.70e-3)
	DTLZ6	3	8.028e+1(5.93e-1)	5.931e+1(8.45e-1) 6.195e+1(7.15e-1)	5.762e+1(1.31e+0)	5.742e+1(1.19e+0)
	DTLZ7	3	7.117e-1(8.60e-2)	4.126e-1(4.30e-2) ≈ 2.652e-1(5.53e-2) 3.284e-1(6.55e-2)	3.22e-1(6.67e-2)
	DTLZ1	5	≈2.963e+2(3.10e+1) 9.746e+2(5.41e+1)9.066e+2(8.90e+1)	2.063e+3(1.07e+2)	3.162e+2(2.80e+1)
	DTLZ2	5	4.084e-1(3.25e-2)	2.647e-1(1.36e-2) 3.215e-1(1.02e-2)	1.991e+0(1.59e-1)	2.183e-1(7.17e-3)
	C-DTLZ2 5	2.341e-1(1.08-2)	1.321e-1(2.19e-2) 7.703e-1(1.72e-1)	7.186e+1(7.35e-2)	7.354e-2(3.70e-1)
	DTLZ3	5	1.164e+3(9.66e+1)	1.63e+3(9.38e+1) 1.39e+3(1.58e+2)	4.85e+3(3.92e+2)	9.31e+2(8.08e+1)
	DTLZ4	5	5.453e-1(1.45e-1)	2.861e-1(1.85e-2) 4.39e-1(2.55e-2)	5.455e+00(5.55e-1)	2.66e-1(1.62e-2)
	DTLZ1	10 ↑7.579e+1(9.87e+0) 1.191e+3(8.44e+1)2.497e+3(9.62e+1)	2.63e+3(1.30e+2)	6.95e+2(7.36e+1)
	DTLZ2	10 ↑4.22e-1(8.61e-2)	9.72e-1(9.32e-2) 9.28e+0(1.91e+0)	1.25e+1(2.20e+0)	6.156e-1(2.23e-2)
	C-DTLZ2 10 3.603e-1(2.11e-3)	3.679e-1(1.42e-1) 4.209e+1(6.52e+0)	2.19e+2(5.04e+1)	1.109e-1(3.53e-3)
	DTLZ3	10 ↑ 5.171e+2(9.5e+1)	5.353e+3(2.35e+2)1.101e+4(7.91e+2)	1.008e+4(7.14e+2)	2.44e+3(2.43e+2)
	DTLZ4	10 5.127e-1(8.07e-2)	1.71e+0(1.12e-1) 4.761e+0(8.32e-1)	7.421e+0(3.68e-1)	4.913e-1(3.53e-2)

The source codes are available at https://github.com/AdwLab/AGE-MOEA-plus-plus (a) Concave DTLZ

(b) Linear DTLZ 1. (c) Convex DTLZ 2. (d) Disconnected DTLZ 7.(e) Degenerate DTLZ 5.

ACKNOWLEDGMENTS

This work is funded by the ANRT CIFRE Program (2019/0877).


Matthias S Ober

Daniel Dermody

Mathieu Maillard

Franck Amiot

Géraldine Malet

Benjamin Burger

Caroline Woelfle-Gupta

Bruno Berge

Geŕaldine Malet

Development of Biphasic Formulations for Use in Electrowetting-Based Liquid Lenses with a High Refractive Index Difference

Keywords: liquid lenses, electrowetting effect, statistical design of experiments, simultaneous optimization ACS Combinatorial Science

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Development of Biphasic Formulations for Use in Electrowetting-

Based Liquid Lenses with a High Refractive Index Difference to -90 °C with a temperature gradient of -10 K/min, equilibrated at -90 °C for 1 min and finally heated to 45 °C with a temperature gradient of 10 K/min. Generally, the maximum of melting peak endotherms, m p , rather the onset of melting m o , was used as a melting point. If several phase transitions or melting points were observed, the one with the highest temperature was used as the data point for model fitting.

Turbidity. The turbidity of the oil phase/conductive phase system after thermal shock was determined by high-throughput nephelometry measurements. Turbidity originates from partial miscibility from liquids at high temperature, after a thermal cycling, condensation occurs in oversaturated liquids, forming light scattering dispersion of immiscible liquids. [START_REF] Kerker | The Scattering of Light and Other Electromagnetic Radiation[END_REF] Turbidity is characterized by two main parameters, the maximum turbidity related to the amount and size of liquid dispersion, and the recovery time necessary to retrieve liquids within transparency specifications.

First, 150 μL of the hydrophilic and 100 μL of the hydrophobic phase were pipetted into a 96-well optical glass plate that was sealed with a corresponding 96-well silicone lid. The plate was then heated to 85 °C (compound selection stage) or 70 °C (formulation DOE) for 20 min in a convection oven.

After temperature equilibration, the plate was removed from the oven, the lid was removed and the hot plate was immediately placed into an Ascent Nepheloskan nephelometer (Thermo Labsystems). The turbidity of the oil/conductive phase binary systems of the plate was measured at 5 min intervals with an integration time of 1000 ms/well. Proprietary formulations with a known behavior, as well as water blanks were used as internal standards. At least two, usually three independent measurements of each sample were performed and the turbidity/time curves were averaged among these measurements.

Viscosity. Viscosities were measured at variable temperatures on a Brookfield CAP 2000+ instrument with the appropriate spindle for the expected viscosity range. High-throughput measurements of viscosities were performed by using the Dow TADM Viscosity experiment, [START_REF] Deshmukh | A Novel High-Throughput Viscometer[END_REF] utilizing a Hamilton Microlab Star dispensing robot and measuring the pressure drop during aspiration. All high-throughput viscosity measurements were performed at room temperature. Solubility (Semiquantitative). Experiments were conducted combinatorially using a 96 well plate containing individual 1.5 mL vials. A small amount of any of the solid candidate compound series (∼5 mg) was combined with a series of liquids of interest. The vials were sealed and the plate was shaken overnight at room temperature. The next day, the vials were visually inspected for solubility of the solid.

High Throughput Powder Dispensing. Solid components of high throughput experiments were dispensed using an Autodose Powdernium (Symyx) solids handler. The robotic handler was used in a Many to Many configuration with on-deck weighing. Twenty-five or fifty milliliter hoppers were used as sources and the powders were dispensed into 1 mL vials.

Design of Experiments, Surface Modeling, and Data Visualization. Constrained D-optimal mixture designs were calculated with JMP, 15 which was also used to perform multivariate regression and model development with the acquired data points. JMP was further utilized to search for optimal areas in the formulation space, sorting and clustering of data, generation of bivariate plots and phase diagrams.

Numeric data, as received from turbidity or miscibility screens, was arranged into a matrix (ASCII format) and plotted utilizing GnuPlot. [START_REF] Williams | Gnuplot 4.4: An Interactive Plotting Program[END_REF] Library composition diagrams were generated with Freeslate Library Studio. 17 ■ RESULTS AND DISCUSSION Formulation Development. Biphasic formulations consist of a low refractive index conductive hydrophilic phase and a high refractive index, nonconductive hydrophobic phase. The formulation of these phases for liquid lens applications have to meet a variety of physicochemical requirements. The most important constraints are (1) low cross miscibility (<0.2%) at the temperature range of 20 °C and +85 °C, (2) identical densities of both phases, (4) turbidity after thermal shock in contact with water phase as low as possible (ideally <20 NTU), (5) viscosity as low as possible (ideally <10 cP), ( 6) melting point below -25 °C, and (7) refractive difference as high as possible. High refractive indices are obtained by introducing formulants containing structural features positively correlated with a high refractive index in existing QSPR models, [START_REF] Katritzky | Quantitative Correlation of Physical and Chemical Properties with Chemical Structure: Utility for Prediction[END_REF] such as a unsaturation, aliphatic and/or aromatic ring structures, a high Balaban index [START_REF] Balaban | Highly Discriminating Distance-Based Topological Index[END_REF] (correlated with molecular compactness and branching), a high dipole moment, [START_REF] Katritzky | Quantitative Correlation of Physical and Chemical Properties with Chemical Structure: Utility for Prediction[END_REF][START_REF] Ha | Quantitative Structure-Property Relationship (QSPR) Models for Boiling Points, Specific Gravities, and Refraction Indices of Hydrocarbons[END_REF][START_REF] Xu | Prediction of refractive indices of linear polymers by a four-descriptor QSPR model[END_REF] and high atomic number elements. Therefore, refractive index is often positively correlates with the density of the fluid. Because of its unique electronic properties, fluorine tends to lower the refractive index of a compound, [START_REF] Katritzky | Correlation and Prediction of the Refractive Indices of Polymers by QSPR[END_REF] while still increasing the material's density. The highest refractive index optical fluids (n ≥ 1.6) are often brominated or chlorinated compounds, but their lightinduced homolysis/radical formation, and sensitivity to hydrolysis usually limit their use in liquid lenses.

ν = + (2)
wherein m is the slope of the regression line. We have experimentally determined apparent refractive index/density quotients for sodium bromide, potassium acetate, sodium triflate and sodium trifluoroacetate (Figure S1). With this, the refractive index of the conductive phase can be predicted by the equation

n wn w i i i j j j D,conductivephase D, D, ∑ ∑ ν ≈ + (3)
wherein w i and w j are the weight fractions of the respective liquid or solid formulation components, n D,i are the refractive indices of the respective liquid components and ν D,j are apparent refractive index/density quotients of the respective solid components. Table 1 summarizes n D and v D values of key formulation candidates.

Density

Model. An empirical density model for the multicomponent conductive phase had been previously developed by Varioptic SA and was expanded for the components used in this study. The development and characteristics of this model is beyond the scope of this report. 24 Freezing Point Model. We used statistical design of experiments (DOE) to develop an empirical model of the phase transition behavior of the conductive phase as a function of formulation. In order to obtain a high-resolution model with a manageable number of required experiments, it was necessary to limit the number of formulation ingredients. To simplify the system, KOAc and NaOTf were excluded due to their similar formulation behavior of NaBr and NaTFA, respectively. Permitted formulation components were water, EG, TMG, NaBr and NaTFA. We designed 96 experiments to explore the freezing points of formulations containing these five candidate compounds at any weight fraction within the approximate solubility boundaries outlined in Table 2. These 96 experiments provide sufficient resolution to explore main interactions plus second and third degree interactions and allow room for a number of standards, center points and repeats (Figure S3 and Table S7).

All samples were measured by DSC, melting points were determined within a range of -90 and 25 °C (Figure S3). As expected, not all of the parameters were statistically significant in the model, and we removed parameters with the highest P-values until R adj 2 did not increase any further and was within proximity of R 2 . Data was fitted the following model equation optimized to exclude statistically insignificant coefficients: 4). Twenty-two compounds in our selection were liquids, and 11 were solids.

Turbidity in Binary Systems. One important requirement of the application is that the system needs to remain clear both upon contact of the two phases and after thermal shock.

We S10 andS11.

Additional Selection Criteria. To achieve short response times of the liquid lens, lower viscosity candidate compounds were given preference. As discussed above, the lowest refractive indices for an optimal water phase can be found at a density range of 1.25-1.40. To achieve the highest possible refractive 449 index difference, this density range is ideally matched by the 450 high refractive index oil phase. As most of the candidate 451 compounds had lower densities, both a high density and a high 452 refractive index were important selection criteria. In addition to 453 these parameters, emphasis was placed on compounds with a 454 low toxicity. To this end published MSDS toxicity data was 455 used as a selection criterion. Table S12 summarizes the assess-456 ment of all selection criteria of the evaluated compounds.

457

Solubility of Solid Compounds and Influence on 458 Refractive Index. The goal of this study was to identify solid 459 candidates that are soluble in the oil phase while simulta-460 neously raising its refractive index. Each solid (Figure 4b) 461 was combinatorially prescreened for solubility in a selection 462 of major oil phase components. In addition, we tested for 463 undesirable solubility in EG, TMG and water. Among the solid 464 candidates, only naphthalene 30, thianaphthene 31, vinyl-465 triphenylsilane, 1,3-dimethyl-1,1,3,3-tetramethyldisiloxane 9, 466 dibenzothiophene 32, and diphenyl phthalate 29 warranted 467 further investigation. The results of this prescreening study is 468 shown in Table S13. 469 We next tested the remaining candidates for their capacity 470 to increase the refractive index of the oil phase formulation. 471 To this end, each of these compounds was dissolved combina-472 torially in a series of oil-phase candidates to saturation at 473 20 °C. The refractive index offset of the saturated solutions in 474 comparison to the pure liquid compounds was then measured. 475 Table 3 summarizes the result of these experiments. Unfor-476 tunately, in solution, none of the compounds had a higher 477 refractive index than the Santovac/Santolight oils 18 and conse-478 quently lowered the refractive index in a mixture with SL-5267. For each designed formulation of the oil phase library, we calculated an ideal composition for the corresponding conductive phase by constraining the water phase model developed above to a density of d oil , a melting point of -25 °C and iteratively resolve the model for the minimum achievable refractive index. Figure 8 shows a graphical representation of the designed oil-and conductive-phase libraries. Numeric tables with the experimental formulation compositions are listed in Tables S14 andS15. Because, as shown above, the composition of an ideal low 548 refractive index water phase (consisting of water, EG, TMG, 549 NaTFA, and NaBr as the only permitted components) can 550 be unambiguously derived for a given density, the only input 551 variables needed to characterize each complete biphasic 552 system, are the weight fractions of each oil phase components. 553 This fact greatly simplifies the optimization of parameters that 554 depend on the composition of both phases, such as refractive 555 index difference, turbidity formation or phase cross miscibility.

556

Refractive Index Difference Model and Mutual Miscibility. 557 An important goal of this study was to optimize the refractive 558 index difference between oil and conductive phase and thereby 559 the optical strength of Varioptic lenses. At the same time, mis-560 cibility and extraction of components into the opposite phase 561 should be kept to a minimum. Miscibility changes induced by 562 temperature changes are problematic as demixing upon cooling 563 may cause temporary or permanent turbidity in the liquid lens. 564 To generate data to develop empirical models for these two 565 parameters, refractive indices of both phases of each library 566 sample were measured in the initial state and after equili-567 bration at 20 and 70 °C with the opposite phase.

568

For the refractive index model, the refractive index difference 569 of the equilibrated oil phases and corresponding water phases 570 at 20 °C (Table S16) was fitted against a model parametrized 571 for linear and two-way cross terms of the weight fractions of 572 potential formulation components of the oil phase. We per-573 formed multivariate regression and model simplification (i.e., 574 elimination of insignificant cross terms) until R adj 2 did not 575 increase any further. The resulting empirical model fits the data 576 tightly (R 2 = 0.9996, R adj 2 = 0.9992) and predicts the refractive 577 index difference as a function of the weight fractions of 578 all permitted oil phase formulation components (Table S18). 579 Because of the initial constraints placed on the experimental 580 design, a substantial amount of samples (∼30%) had a desir-581 able refractive index difference (Δn D ) greater than 0.290. 582 Miscibility estimates at 20 and 70 °C were calculated from 583 the change of refractive index of the pure oil phase before and 584 after equilibration (overnight) with the water phase at the 585 respective temperatures as described in the Experimental 586 Techniques. The results of these measurements are presented 587 in Table S18. We fitted these data points to the underlying 588 DOE model and removed significant terms until R adj 2 did not 589 improve any further. The final models had acceptable R-values and satisfactory statistical validity tests (Tables S19 andS20).

They predict miscibilities between the oil phase and its corre-

π π Σ ∝ ≈ - i k j j j j j y { z z z z z (6)
that is, decreases as droplets grow in size.

The overall recovery trend is now a balance between migration kinetic and droplet coalescence: it is critical that droplets coalesce to decrease the number of them, but at the same time, droplet diffusion coefficient due to Brownian motion is inversely proportional to diameter. Thus, if the droplet coalescence process is too fast, droplet migration speed becomes too slow to migrate back to their original phase, leading to a diffusing plateau because of nonmobile droplets. The optimum here is thus a compromise between nucleation, coalescence and migration speed and the complexity of the solution emphasizes the need of a high-throughput approach.

Individual data from this experiment for every library sample (average of all individual measurements taken between 240 and 300 min) are presented in Table S17. To improve the quality of the fit, turbidity data in NTU was logarithmized before regression. The final model, consisting of main terms and statistically significant second degree cross terms can be found in Table S21.

Empirical Model of Viscosity. The viscosity of each sample of the oil phase library was determined by averaging the results of three independent runs of the Dow TADM high throughput viscosity screening experiment, performed on a Hamilton Microlab Star robot (Figure S12). The accuracy of this measurement was confirmed to be within ±5% by amending the library with a series of Brookfield viscosity standards ranging between 0.9 and 81 cP. None of the library samples had a viscosity below the desired 10 cP, but the best formulations came with approximately 14 cP close to this threshold. Individual data is presented in Table S17. The final model obtained after multivariate regression is presented in Table S22.

Empirical Melting Point Model. We measured the melting temperature of each oil phase in the library by DSC. The results of these targets are presented in Table S17. Most of the melting points were below -40 °C. As before, we performed regression and model simplification until R adj 2 did not increase any further. The final model has a good fit to the experimental data (Table S23).

Optimization of the Biphasic System. To find areas in the formulation space that simultaneously meet the optimization targets, a Monte Carlo simulation was performed by calculating predicted properties of 2 × 10 5 random biphasic formulations within the original definition range (Table 4). All formulations that did not meet the desired requirements were discarded. Selection criteria were as follows: turbidity after temperature shock after 4 h cool down, T < e 1.75 ≈ 5.75 NTU; refractive index difference of oil phase and corresponding ideal water phase Δn D > 0.295; phase miscibility at room temperature, φ waterphaseinoilphase (20 °C) < 0.15%; miscibility at 70 °C, φ waterphaseinoilphase (70 °C) < 0.3%; melting point of oil phase, m p < -35 °C. As no simulated formulation had a viscosity of v < 10 cP, the requirement was relaxed to v < 16 cP. Duplicate hits were removed and the remaining data set consisted of 124 simulated formulations. The data points were sorted and grouped into seven clusters by using Ward's minimum variance method (Figure 9). [START_REF] Ward | Jr Hierarchical grouping to optimize an objective function[END_REF] Simulated ideal formulations in cluster 1 contain about 40% 1-chloronaphthalene 6 and 60% Santolight SL-5267 18. Formulations in cluster 2 include the same compounds at weight fractions of ∼30% and ∼40%, respectively. Additionally, ∼20% ethylthiobenzothiazole 21 and ∼10% of other compounds, that is, 4-bromodiphenylether 16 are included. Formulations in clusters 3 and 4 do not contain 1-chloronaphthalene 6, but instead, 1-phenylnaphthalene 8 with a weight fraction of ∼30-40%, supplemented by ∼20-30% Santolight SL-5267 18. The remaining compounds are ∼10% 2,5-dibromotoluene 33, ∼15% diphenyl sulfide 17. and thianaphthene 31. Interestingly, diphenyl sulfide 17, which performed poorly in our turbidity screen, is in most simulated formulations accompanied by 2,5-dibromotoluene 22. 2,5-Dibromotoluene 22 increases the density of the oil phase, thereby permitting a larger salt fraction and density of the water phase. The higher salt fraction increases the polarity of the conductive phase and this effect may help counteracting a predicted turbidity associated with

First, the

  refractive indices n D,pure of the pure hydrophobic and corresponding hydrophilic phases were determined at the temperature of interest as described. Next, the phases were combined and equilibrated by agitating the mixture overnight at the prescribed temperature. The refractive indices n D,eq of the top and bottom phase were remeasured at an identical temperature. The volume fraction of phase 1 in phase 2 as a function of temperature was estimated by using the equation behavior of the refractive indices (Arago-Biot approximation). 12 Miscibility of multicomponent biphasic liquid systems was estimated by treating each phase as a single component (having the refractive index of the multicomponent phase). Density. Densities were measured with an Anton Paar DMA-500 density meter at 25 °C by injecting 1 mL of sample into the measurement cell. Melting Point. Melting points were determined by DSC measurements in a single measurement cycle. A 5 mg of sample, encapsulated into a hermitic aluminum DSC pan, was equilibrated for 1 min at 25 °C. Next, the sample was cooled

Figure 1 .

 1 Figure 1. Refractive index distribution taken from the Varioptic S.A. optical fluid database. A refractive index difference Δn D ≥ 0.290 be obtained only with a limited number of compounds.

Requirements 6

 6 and 7 are counteractive because water antifreeze additives (salts, glycols) raise the refractive index of the hydrophilic phase considerably, thereby limiting the maximum possible Δn D . Matching the densities of both phases while maximizing the refractive index is not trivial. Each formulation adjustment in Part I. Empirical Model and Optimization of the 245 Conductive Phase. Formulation Component Selection. 246 Water is a key component of the conductive phase, both 247 because of its low refractive index and its conductivity after 248 dissolution of salt components. Its key disadvantage for the 249 application is its high freezing point, which needed to be 250 lowered by antifreeze components. To achieve a high refractive 251 index difference, an important goal was to find salts and addi-252 tives that only minimally increase the reactive index of the 253 water phase, while still reduce its freezing and melting point 254 below the designed application range. As freezing point depres-255 sion is, according to Blagden's law, 23 a function of the molality 256 of a solute in a solvent, a high water miscibility or solubility 257 and a low molecular weight is favored. Because of their misci-258 bility and low volatility, common antifreeze agents like 259 ethylene glycol (EG) and 1,3-propanediol (TMG) are suitable 260 agents, but they significantly increase the refractive index of the 261 water phase. 262 In our quest to identify suitable candidate compounds with a 263 lower optical impact, we first investigated low n D fluorinated or 264 perfluorinated aliphatic amines, alcohols or acids, such as triflic 265 or trifluoromethanesulfonic acid. Their high miscibility with 266 the nonconductive hydrophobic phase (resulting in turbid 267 systems) and/or corrosivity led us to quickly abandon this 268 class of compounds. On the other hand, we found that salts of 269 fluorinated acids, such as sodium triflate (NaOTf) and sodium 270 trifluoroacetate (NaTFA) showed promise because of their 271 high water solubility, compatibility with the construction 272 material of the encasing, better optical properties, and a more 273 moderate impact on the refractive index in comparison to 274 sodium bromide (NaBr) or potassium acetate (KOAc), addi-275 tives used in previous formulations. Candidate formulation 276 components of the hydrophilic phase that were ultimately 277 selected for further investigation were water, EG, TMG, NaBr, 278 KOAc, NaTFA and NaOTf. 279 Refractive Index Model. In first approximation, and 280 assuming ideal density behavior, refractive indices n D,i of indi-281 vidual compounds i in a mixture contribute to the refractive 282 index n D,mix in proportion to their respective volume fraction 283 Φ i according to the Arago-Biot model.

Figure 3

 3 Figure3shows a slice through the response hypersurface

  have screened the 21 liquid compounds 1-21 combinatorially in biphasic systems with eight model conductive phases in a full factorial library after mixing and simulated temperature shock (Figure S8). Model conductive phases were water, a 5% (w/w) aqueous NaBr solution, a 20% aqueous NaBr solution, pure ethylene glycol, pure 1,3-propanediol, a 20% aqueous NaTFA solution, a 60% aqueous NaTFA solution, and a 40% aqueous NaOTf solution (all percentages are w/w). This setup not only allowed identification of organic compounds that are robust under a variety of conditions but gave initial information about the effect of different salts on turbidity formation. The plates were heated to 85 °C as described in the Experimental Techniques. Typically, the turbidity initially increased upon cool-down (likely caused by phase demixing and microdroplet formation), followed by gradual clear up. For most systems, the maximum turbidity was reached after 30-40 min. Whereas some system (especially most of the silanes) cleared up quickly or overnight, other systems did not improve or became even more turbid (most notably, 1-bromonaphthalene). The compounds that performed best in the test were the germane Me 3 PhGe 5, and the silanes/siloxanes Me 2 Ph 2 Si 1, 1,3-diphenyl-1,1,3,3-tetramethyldisiloxane 9, and especially, SIP6827.0 11 and DMS-T15 12.Among the compounds with a particularly high refractive index were 1-chloronaphthalene 6, 1-phenylnaphthalene 8, 4-bromodiphenylether 16, and Santolight SL-5267 18. It is also evident that systems that contain salt in the water phase perform better than systems with pure water, which can be explained with a lower phase miscibility due to salting-out effects. The highly concentrated and polar 20% NaBr solution performed better than the 5% NaBr solution or pure water, consistent with published behavior.[START_REF] Maillard | Use of bromine anions in an optical electrowetting device[END_REF] NaTFA solutions performed better than pure water, and a 20% NaTFA behaved, with exceptions, similar to 20% NaBr solutions. The NaOTf solution gave higher turbidities than any other aqueous phase, possibly due to slightly emulsifying characteristics of the large organic anion. Biphasic systems with pure EG and TMG as the simulated conductive phase were among the most turbid upon cooling, with EG being slightly better than TMG in general. This observation corresponds with the highly temperaturedependent miscibility of the glycols with many of the candidate compounds. Figures 5 and S9 summarize the results of this screening study. Experimental values are shown in Table S9. Miscibility with Water and Glycols. We have determined the miscibility of the candidate compounds with water, EG, and TMG by refractive index measurements as described in the Experimental Techniques, assuming additive behavior. Again, the silanes/siloxanes and germanes performed very well in this test. In addition, 1-chloronaphthalene 6, 1-phenylnaphthalene 8, 4-bromodiphenylether 16 and Santolight SL-5267 18 gave promising results. The measurements were conducted at -20 and 70 °C.Figure 6 shows the results, raw data and miscibility calculations are presented in Tables

Figure 2 .

 2 Figure 2. Upper panel, lowest achievable refractive index over a density range of <1.05 to 1.4 for formulations with a maximal freezing point of -25 and -40 °C; lower panel, weight fractions of components of optimal low refractive index formulations with a freezing point of less than -25 °C, over the same density range. Formulations containing sodium bromide are not ideal from a reactive index perspective; TMG is only needed at densities below 1.15.

Figure 3 .

 3 Figure 3. Slice through the response hypersurface along the w EG and w NaTFA coordinates. Ideal low refractive index formulations are for a wide range of densities located on the melting point isoline that corresponds with the low temperature specification limit (e.g., -25 °C). The location of the lowest possible refractive index on the response hypersurface highlighted by a red circle.

Thianaphtene

  31 hat the best solubility in all silanes/germanes and a positive effect on their respective refractive indices. Viscosity of Santovac MCS-293/Thianaphthene Mixtures. During these studies, we have serendipitously discovered that thianaphthene reduces the viscosity of Santovac MCS-293/ Santolight SL-5267 18, without a significant reduction of its high refractive index. To investigate this effect further, we measured the viscosity of thianaphthene/MSC-293 mixtures at different weight concentrations and at various shear rates and temperatures. We have found that an addition of just 10% (w/w) of thianaphthene to Santovac MCS-293 or Santolight SL-5267 lowered the viscosity in all cases by a factor of approximately two (Figure 7). Selection of Compounds for Use in Liquid Lens Formulations. On the basis of these screening results, we have selected eight candidate compounds for potential use in improved oil phase formulations. Dimethyldiphenylsilane 1 was included, as it performed excellently in turbidity/miscibility tests but has a relatively low refractive index. Diphenyl sulfide 17 was selected because of its low viscosity despite its average performance in other screening tests. Santolight SL-5267 18 stood out because of its high reactive index. It performed decently in the turbidity test, especially in combination with high-salt water phases. 4-Bromodiphenylether 16 showed good results in all tests and has a high density, which is important to gain sufficient freedom in formulating oil phases in the desirable density range of 1.25-1.40. Ethylthiobenzothialzole 21, on the other hand, was retained because of its high refractive index and low viscosity despite its less favorable miscibility/turbidity test results. 1-Chloronaphthalene 6 showed

Figure 4 .

 4 Figure 4. Liquid compounds that were screened for use in the oil phase.

Figure 6 .

 6 Figure 6. Miscibility of oil phase liquids candidate compounds with water, EG and TMG. The miscibility was color coded on a logarithmic scale. Dark blue and black designates low miscibility of <0.2%. Violet and pink colors designate miscibility between 1% and 5%, which is less desirable. Brighter colors indicate a miscibility >5%. The miscibility of bis(indenyl)dimethylsilane was not determined at 70 °C, as the biphasic system became too opaque for refractive index measurements.

Figure 5 .

 5 Figure 5. Results of the nephelometry screen of high-refractive index candidate compounds after temperature shock. The turbidity (in NTU) was color coded on a logarithmic scale. Dark blue and black designates turbidities < 5 NTU, which is almost invisible to the naked eye. Pink colors designate turbidities between 10 and 100 NTU, which is clearly visible. Orange or yellow fields have a turbidity of approximately 1000 NTU, which is almost opaque at a depth of 1 cm. Red arrows indicate a proprietary low-turbidity oil phase formulation used as a positive plate standard. The field on the bottom right corner is a proprietary oil/conductive phase system used as negative turbidity standard.

Figure 7 .

 7 Figure 7. Viscosities of thianaphthene/Santovac MCS-293 mixtures at 5 and 20 °C at different shear rates. An addition of 10% thianaphthene reduces the viscosity of MCS-293 by a factor of 2.

a

  Numbers printed in italics indicate a solubility of >50% (w/w) of the modifier in the respective liquid. b Offset of refractive index in comparison to pure liquid after saturation with refractive index modifier. c n.d.: Not determined because of low solubility.

  sponding ideal water phase at 20 and 70 °C as a function of the weight fractions of all included oil phase formulation components. Empirical Model of Turbidity Formation of Equilibrated Biphasic Systems after Temperature Shock. Formation and clearance of turbidity upon cooling was measured by high throughput nephelometry after equilibration of both phases at 70 °C as over the course of 5 h (Figure S10). Most formulations exhibit a strong increase in turbidity as the system cools down and phase separation/droplet nucleation occurs. After the samples reached a maximum turbidity at approximately 30 min, the systems consequently clear up until a plateau is reached, which often remains above the initial turbidity. This plateau was typically reached after approximately 4 h and persisted over prolonged times (>3 days per control measurement). A photograph of the formulations after temperature treatment is shown in Figure S11. This general trend is clearly related to a migration/droplet nucleation phenomenon. Turbidity is caused by nucleation of droplets, inducing Mie scattering of light. 26 Scattering efficiency ∑ is proportional to droplets section where N d describes the number of droplet and a is the droplet radius. As droplets coalesce, the number of droplets decreases and droplets radius increases consequently. Considering as a first approximation a constant total volume of liquid inclusion v = 4N d /3πa 3 , leading to the relation a = (3v/4πN d ) 1/3 , scattering efficiency is expected to vary as

Figure 8 .

 8 Figure 8. Oil and conductive phase library designed to explore the formulation space of liquid lenses' systems with nine candidate compounds for the oil phase and three candidate compounds for the water phase.

  694 diphenyl sulfide. Cluster 5 contains five compounds at similar 695 concentrations: 1-phenylnaphthalene 8, SL-5267 18, ethyl-696 thiobenzothiazole 21, 2,5-dibromotoluene 33, and diphenyl 697 sulfide 17. Cluster 6 relies on high diphenyl sulfide 17 concen-698 trations, accompanied by SL-5267 18 and 2,5-dibromotoluene 699 33. Cluster 7 is the only solution that does not require SL-5267 700 18 to achieve the high refractive index. Instead, high con-701 centrations (∼55%) of 1-phenylnaphthalene 8 are utilized, 702 together with 2,5-dibromotoluene 33 and diphenyl sulfide 17 703 to increase the density and lower the viscosity, which goes 704 along with lower reactive indices of the corresponding ideal 705 water phase. Table 5 summarizes the weight fraction boundaries, 706 which served as guideline for future improvements of Varioptic 707 Liquid lenses. 708 On the basis of these results, prototype highly refracting 709 liquid lenses (Δn D = 0.290) were built by Varioptic SA, which 710 showed outstanding performance properties meeting critical 711 requirements such as optical transmittance after temperature 712 shock, short response time, and minimal optical drift/hysteresis. 713 Although neither specifics about the formulation ultimately used 714 in these prototypes nor detailed application testing data can be 715 disclosed herein, the resulting prototype devices had sufficient 716 optical power to focus at very close distance (<1 cm). At this 717 distance, fine features can be imaged, here demonstrated with 718 a leaf with a 130-140 μm venation pattern (Figure 10). 719 These prototype lenses spearheaded the development of a new 720 generation of commercial Varioptic Lenses with high-power zooming capabilities. 722 ■ CONCLUSION AND SUMMARY 723 Using statistically designed experiments in conjunction with 724 high throughput/automated measurements, we have developed 725 several statistical models to predict a variety of physicochem-726 ical properties of the both phases of formulations that can be 727 used in Varioptic Lenses. These parameters include refractive 728 index, melting point and density of the water phase, the viscos-729 ity and melting point of the oil phase, the cross miscibility at 20 730 and 70 °C, and the refractive index difference of the biphasic 731 system. Additionally, we have demonstrated that fluorinated salts 732 in the conductive phase, such as sodium trifluoroacetate, lower 733 the achievable refractive index while improving transmittance 734 after temperature shock. We developed a method to calculate 735 ideal low reactive index conductive phases to match the density 736 of any given high refractive oil phase. Next, nine oil phase 737 candidate compounds were identified. With these candidate

Figure 9 . 1 (Figure 10 .

 9110 Figure 9. Formulations found in the Monte Carlo simulation that meet optimization goals. Colors specify weight fractions of formulation components (blue, -0%; gray, ∼10%; gray-pink, ∼30%; bright red, ∼60%).

Table 1 .

 1 Refractive Indices of Candidate Liquids and Apparent Refractive Index/Density Quotients of Candidate Salts

	liquid candidate components	n D,i	solid candidate components	ν D,i
	certified water a	1.33273	sodium bromide c	1.4660
	ethylene glycol b	1.4310	potassium acetate c	1.4527
	1,3-propanediol b	1.4393	sodium triflate c	1.3679
			sodium trifluoroacetate c	1.3721

). Because the 361 refractive index generally increases with the content of antifreeze 362 agents and salts, ideal formulations with the lowest possible 363 refractive index can be found at the isoline of the maximum 364 allowable melting point (i.e., -25 °C) for densities between 365 approximately 1.15 and 1.31 g/mol. For densities above 1.31, 366 ideal formulations are binary mixtures of water and NaTFA at 367 the indicated ratios. The lowest achievable refractive index of a 368 conductive phase that does not freeze at -25 °C is approxi-369 mately 1.35; the intersection point of the -25 °C melting 370 point isoline intersects with the x-axis (0% EG, a solution of 371 ∼43% (w/w) NaTFA in water). 372 Part II. Candidate Compound Selection of the 373 Hydrophobic Phase. Formulation Component Preselec-374 tion. To find suitable formulation components for the 375 hydrophobic phase, we have tested 32 high-refractive index a Refractive index of distilled water available in-house against Cargille reference standard. b Remeasured in-house. c Determined by regression as described in text.

Table 2 .

 2 DOE Boundaries of the Water Phase (Mixture Design) a NaTFA , w water > 3w NaBr , and w water + w EG + w TMG + w NaBr + w NaTFA = 1

		weight fraction (w i )
	compound	low boundary	high boundary
	water	0.3	0.8
	ethylene glycol	0	0.8
	1,3-propanediol	0	0.4
	sodium bromide	0	0.333
	sodium trifluoroacetate	0	0.5

a Additional constraints: w water > w D 376 compounds that had been preselected for characteristics 377 including high refractive index, low toxicity, and low volatility. 378 Compounds featuring a high refractive index contain typically 379 either at least one heavy heteroatom or aromatic ring system in 380 their structure. Our final compound selection included 381 materials from the following groups: aromatic silanes, aromatic 382 siloxanes (except for the polydimethoxysilane DMS-T15, 383 which is completely aliphatic and has a low refractive index), 384 aromatic germanes, naphthalenes, fused aromatic rings and 385 heterocycles, phenyl ethers and halogenated aromatic compounds, or aromatic esters (Figure

Table 4 .

 4 Constraints of Mixture Design for Oil Phase Formulations

												constraint
	description											constraint equation	value
	sum of weight fractions must be 1 (inherent constraint of mixture design)	w	w 1 phenylnaph w Me Ph Si Ph S 2 2 2 + + + + w w SL5267 2,5 dibromotol w 4 BrPhOPh w thianaphthene w ethylthiob + + --+ -	+	w 1 chloronaph -	1
		w 1.5640	Me Ph Si 2 2	+	w 1.6382	Ph S 2	+	w 1.6740 SL5267	+	w 1.6086 4 BrPhOPh -
	theoretical refractive index above 1.65 to tilt design space toward high Δn D		+	w 1.6593 ethylthiob	+	w 1.6329 1 chloronaph -	+	w 1.6638 1 phenylnaph -	≥1.65
			+	w 1.602 2,5 dibromotol -	+	w 1.64 thianaphthene
		w 0.99	Me Ph Si 2 2	+	w 1.11	Ph S 2	+	w 1.20 SL5667	+	w 1.42 4 BrPhOPh -	+	w 1.23 ethylthiob
	theoretical density (assuming simple additive behavior) above 1.19, to tilt design space toward a density range permitting low n D water phases		+	w 1.19 1 chloronaph -	+	w 1.10 1 phenylnaph -	+	w 1.815 2,5 dibromotol -	≥1.19
			+	w 1.20 thianaphthene
	limit Santolight SL-5267 and 1-phenylnaphthalene weight fractions to	w SL5267							<0.6
	prevent high oil phase viscosities	w 1 phenylnaph -						<0.6
	limit thianaphthene to prevent crystallization at low temperatures	w thianaphthene						<0.4
	limit thianaphthene below the weight fraction of Santolight SL-5267	w thianaphthene	-	w SL5267	≥0

Table 3 .

 3 Change of Refractive Index of Several Liquids (Rows) after Saturation with Various Solid Refractive Index Modifiers (Columns) a

	refractive index modifier (solid)

phase 1 to change its refractive index will affect its density and

require an adjustment in phase 2. In turn the second phase's

refractive index, and thereby Δn D changes in a nonobvious

fashion. In order to address this problem, it is necessary to

DOI: 10.1021/acscombsci.8b00042 ACS Comb. Sci. XXXX, XXX, XXX-XXX