Machine learning-based inverse materials discovery has attracted enormous attention recently due to its flexibility in dealing with black box models. Yet, many metaheuristic algorithms are not as widely applied to materials discovery applications as machine learning methods. There are ongoing challenges in applying different optimization algorithms to discover materials with single- or multi-elemental compositions and how these algorithms differ in mining the ideal materials. We comprehensively compare 11 different optimization algorithms for the design of single- and multi-elemental crystals with targeted properties. By maximizing the bulk modulus and minimizing the Fermi energy through perturbing the parameterized elemental composition representations, we estimated the unique counts of elemental compositions, mean density scan of the objectives space, mean objectives, and frequency distributed over the materials’ representations and objectives. We found that nature-inspired algorithms contain more uncertainties in the defined elemental composition design tasks, which correspond to their dependency on multiple hyperparameters. Runge–Kutta optimization (RUN) exhibits higher mean objectives, whereas Bayesian optimization (BO) displayed low mean objectives compared with other methods. Combined with materials count and density scan, we propose that BO strives to approximate a more accurate surrogate of the design space by sampling more elemental compositions and hence have lower mean objectives, yet RUN will repeatedly sample the targeted elemental compositions with higher objective values. Our work sheds light on the automated digital design of materials with single- and multi-elemental compositions and is expected to elicit future studies on materials optimization, such as composite and alloy design based on specific desired properties.
I. INTRODUCTION
Designing materials with desired tailored properties has always been a key goal in human society, from metallurgy, biomedicine, energy storage, to nanotechnology. Traditionally, most research studies conducted in the field of designing materials with multi-elemental compositions rely on costly experimental synthesis and time-consuming trial-and-error processes. Advances in computational engineering can greatly accelerate materials design efficiency. To achieve the design goal, many optimization methods were developed for the design process. These methods can be categorized as first and zeroth orders, where the first-order optimization methods are gradient-based methods, such as gradient descent and Newton–Raphson, that rely on the problem being first-order differentiable. The first-order methods are widely applied to topology optimization coupled with finite element methods for mechanical design, including classical problems such as the design of beams1 and trusses.2 However, the primary drawback of first-order optimization is that the differential form of the problems may not exist, or be easily obtained, in many (or most) real-world application scenarios. In many cases, the governing equations for the problem do not exist as differentiable forms. Moreover, first-order optimization tends to get “trapped” in local optima close to the initial design, which heavily relies on prior physical understanding of the problem to ameliorate. Hence, metaheuristic optimization methods, such as genetic algorithms and particle swarm optimization, are more flexible and adaptive to complex problem domains. These characteristics enable metaheuristics to be applicable to multiscale materials design problems, such as molecular simulations, machine learning (ML) surrogates, and density functional theory calculations. For example, molecular dynamics (MD) simulations model the interactions of particles based on interatomic potentials,3,4 DFT simulations calculate the energy based on pseudo-potentials for the Hamiltonian,5 and ML surrogates represent the material’s structure–property map with a neural network.6 In these situations, obtaining first-order derivatives is either an extremely taxing task or may even be intractable. To overcome these challenges, recent literature highlights the utility of metaheuristics for elemental composition design, in contrast to gradient-based optimizations. For example, Ramprasad and co-workers developed genetic algorithm methods coupled with deep learning surrogate models for designing polymer dielectrics7 and used similar strategies to design polymers with high glass transition temperatures and bandgaps.8 Winter et al. combined machine learning surrogates with particle swarm optimization (PSO) to design drugs.9 Weiel et al. used PSO to select parameters in molecular dynamics simulations.10
Many machine learning methods exist for design optimization, and two of the most widely applied methods are Bayesian optimization (BO) and deep reinforcement learning (DRL). BO leverages Gaussian process regression (GPR) to fit a surrogate model of the design space and leverages Bayesian statistics to iteratively improve the accuracy of such a surrogate by sampling more points via an acquisition function. BO is frequently applied in materials design. For example, Shin et al. combined the micromechanical experiments and BO to design spiderweb nanomechanical resonators.11 Chen and co-workers developed a series of BO frameworks for mechanical design and elemental composition design using finite element simulations,12 molecular simulations,13,14 and machine learning surrogate models.15,16 Zhai and Yeo combined BO and individual-based modeling of bacteria to design antimicrobial nanosurfaces.17 In another popular method, DRL utilizes a deep neural network (DNN) as an agent to interact with the environment to receive rewards from different output properties and learn the best strategy to achieve the design goal. When applying DRL for materials design, the evaluations of the material properties are treated as the environment for the agent to learn the correct policy to achieve the design goal. For instance, Gu and co-workers developed a DRL framework for designing tough graphene oxide18 and continuum composite materials.19 Barati Farimani and co-workers developed DRL frameworks for porous graphene for water desalination.20
Can these optimization methods be uniformly applied to the same problem and how do they differ? To answer these questions, benchmark studies are urgently needed. One particular reason for the urgency is efficiency. Optimization methods are intended to improve the efficiency of the design problems. However, when companies or research groups want to employ design optimization for their problems of interest, it is of utmost importance to select optimization methods that are best suited to their specific problems. Selecting ideal optimization methods from the massive amount of the existing literature is non-trivial: one may fall into the “trial-and-error” trap, and the efficiency of solving the problem is significantly reduced. Moreover, ML-based studies for materials design are often focused on either large-scale molecular structures (i.e., polymers) or continuum structures. This is in part because it is easier to extract structure–property relationships for larger structural bodies in the context of continuum modeling, where first-order derivatives are easy to obtain. Specifically, for materials composed of multiple elements, many studies leverage molecular graphs or use natural language processing to find this structure–property map. In contrast, while there are some studies on the design of small molecules using language representation with generative models,21 a huge gap exists in applying general optimization methods for designing such multi-elemental compositions. Therefore, following our previous discussion, the specific problem to be posed herein is how different optimization methods differ in the context of materials design.
To address these challenges, we propose and conduct a benchmark study by designing a framework to implement many state-of-the-art optimization methods. We benchmark many widely employed metaheuristic optimizations, including Bayesian optimization, genetic algorithms, particle swarm optimization, and simulated annealing. We further compared novel optimization methods, such as the arithmetic optimization algorithm and Runge–Kutta optimization, for a total of 11 different optimization algorithms. We took advantage of the large number of open-source software applications for computational elemental composition design, including pymatgen,22 Materials Project,23 matminer,24 and mealpy.25 The material evaluations are enabled by surrogate models: pretrained graph neural network models to predict the elemental compositions’ properties based on input graph structure. The general schematic of our study is shown in Fig. 1. We begin with inspiration from the general premise of industrial or fundamental research problems: one needs to design elemental compositions based on specific applications. Given this application, the targeted physical properties must be identified based on the design requirements, such as high toughness or low thermal conductivity. Based on these targeted applications and properties, the design problem can then be formulated mathematically together with the corresponding optimization algorithms. The design space can be explored to extract the optimal elemental compositions from the materials database(s). Here, we will mainly focus on the optimization process [Figs. 1(c)–1(e)] to propose simplified target material properties and focus on the algorithmic implementations and benchmarking. The significance of this work lies in that there is a lack of studies that evaluate the performance of different optimization methods, in contrast with the abundance of published literature that applies specific optimization methods for tailored materials design tasks. Our study may benefit future research in the following ways: (1) Researchers can select different algorithms that correspond well with their potential applications based on our evaluation results, and (2) researchers can use or refer to our framework and benchmark other optimization algorithms of their interest. Our major contributions are to (1) propose a framework that can implement different algorithms, which generalizes the materials’ inverse design tasks, (2) evaluate different optimization methods based on defined metrics, and (3) highlight ideal algorithms for different potential applications based on our evaluations. Some results were unexpected: the commonly employed BO did not stand out in terms of mean objective evaluation, while surprisingly RUN outperformed both the stability and objective evaluations.
General overview of the materials design process. (a) Schematic for defining the intended industrial application. (b) The process of identifying the desired material properties, such as high bulk modulus and low Fermi energy. (c) The schematic for formulating a digital design optimization problem using numerical models to discover or design elemental composition with the desired properties. (d) Different optimization methods use exploration/exploitation strategies to generate a design space for exploring the best materials. (e) The optimal elemental compositions are extracted from the design space exploration and are ready for real-world applications. In this study, we focus on formulating an optimization framework that enables the application of the most current popular optimization methods to search for desired elemental compositions and compare different optimizations based on specific metrics. Finally, the “mined” elemental compositions are analyzed.
General overview of the materials design process. (a) Schematic for defining the intended industrial application. (b) The process of identifying the desired material properties, such as high bulk modulus and low Fermi energy. (c) The schematic for formulating a digital design optimization problem using numerical models to discover or design elemental composition with the desired properties. (d) Different optimization methods use exploration/exploitation strategies to generate a design space for exploring the best materials. (e) The optimal elemental compositions are extracted from the design space exploration and are ready for real-world applications. In this study, we focus on formulating an optimization framework that enables the application of the most current popular optimization methods to search for desired elemental compositions and compare different optimizations based on specific metrics. Finally, the “mined” elemental compositions are analyzed.
This paper is articulated as follows: We present an overview of the underlying methodology, describe the workflow of connecting machine-learned models of materials properties with different optimization methods, and benchmark the results obtained from 11 different optimization algorithms in Sec. III. We then discuss the algorithms that are developed to perform the automated optimizations in Sec. IV. A detailed explanation of the methodologies employed is presented in Sec. II, and further information is contained in the supplementary material.
II. METHODS
A. Inverse optimization framework
1. Problem formulation
The optimization problem can be solved using various optimization algorithms. Here, we explored 11 different state-of-the-art optimizations to conduct benchmarking for materials discovery. The methods include Bayesian optimization,27 deep reinforcement learning,28 genetic algorithm,29 particle swarm optimization,30 hybrid gray wolf optimization,31 (hybrid-improved) whale optimization algorithm,32 ant colony optimization,33 differential evolution,34 simulated annealing,35 bees optimization,36 arithmetic optimization,37 and Runge–Kutta (RUN) optimization.38 These algorithms fall into two main categories: machine learning-based optimization and metaheuristic optimization [Figs. 2(b) and 2(c)]. We use different strategies to perturb Θ to maximize , which will be elaborated in Sec. II and the supplementary material. Here, we consider single- and multi-element composition optimization scenarios. For single-element composition optimization, the design variable natom is restricted to 1, and there are only two tunable design variables ξn and η. For multi-element composition optimization, the screened elemental compositions contain natom ∈ [1, 4]. The materials libraries as indicated in Fig. 2 consist of a series of MP IDs, and the randomly generated η is used to select the corresponding MP ID.
Schematic for optimization benchmarking and different optimization algorithms employed in this study. (a) Workflow for discovering elemental compositions by combining graph neural network (GNN) surrogate models with defined design variable randomization to select optimal elemental compositions with targeted materials properties. See Sec. II and the supplementary material for details. (b) Schematic for optimization through metaheuristic algorithms, starting with a set of random particles (or population) sampled in the design space, to iteratively update and eventually cluster around globally optimal values, indicating the search for elemental compositions with desired properties. (c) Bayesian optimization methods for elemental composition discovery, seeking to build a surrogate model of the design space by employing Gaussian process regression and searching for the next material evaluations based on an acquisition function. Optimal elemental compositions can be extracted from the design space surrogate.
Schematic for optimization benchmarking and different optimization algorithms employed in this study. (a) Workflow for discovering elemental compositions by combining graph neural network (GNN) surrogate models with defined design variable randomization to select optimal elemental compositions with targeted materials properties. See Sec. II and the supplementary material for details. (b) Schematic for optimization through metaheuristic algorithms, starting with a set of random particles (or population) sampled in the design space, to iteratively update and eventually cluster around globally optimal values, indicating the search for elemental compositions with desired properties. (c) Bayesian optimization methods for elemental composition discovery, seeking to build a surrogate model of the design space by employing Gaussian process regression and searching for the next material evaluations based on an acquisition function. Optimal elemental compositions can be extracted from the design space surrogate.
Following the discussion from Eq. (1), the optimization is benchmarked and connected through the material evaluations. We utilize MEGNet for materials properties evaluation mainly for two considerations: (1) Fast inference: by using the pretrained MEGNet model (mp-2019.4.1), the properties can be obtained “on the fly” given the input graph structure, getting rid of the lengthy computational process via traditional numerical methods such as molecular dynamics (MD) or density functional theory (DFT). (2) Better screenability of the design space. Given the pretrained MEGNet model, one can infer the corresponding properties given the generated input graph , hence allowing the evaluation of a large amount of “untouched” materials. Meanwhile, due to the lack of existing data for the properties of many potential materials, one may not directly obtain the properties given the input elements just from a query. In addition, due to the lack of pseudo- and/or interatomic potentials, MD or DFT calculations may not be conducted. Here, we employ MEGNet due to its reported accuracy and stability. As reported by Chen et al.,26 the MAE for the prediction of Fermi energy and bulk modulus is around 0.028 and 0.05, respectively. MEGNet is not the only option for materials properties inference; other existing frameworks with a competitive performance include CGCNN,39 MODNet,39 H-CLMP(T),40 and crystal graph attention networks.41 Since our focus is on optimization methods, we do not emphasize the surrogate model evaluations. Notably, many recent studies focus on materials property prediction benchmarks.42,43
For this problem, the inverse optimization problem is to directly perturb the graph structure such that the materials’ properties can be tailored. However, it is difficult to represent different materials’ graph structures containing different elemental numbers in the same dimensions for optimization. Hence, we parameterize the graph structure so that the algorithms tune the design variables Θ = [natom, ξn, η] to change the elemental composition . Here, as visualized in Fig. 2(a), the optimization begins with a variable natom, as the number represents different types of elements in the elemental compositions. We may denote the different atoms in the elemental compositions as the atomic dimensions. Associated with this natom,we define a random number distributed in [0,100] as ξn to assign atomic types to the atomic dimensions, which we may define as the material basis. For example, if the atomic dimension natom = 3, we may obtain an element basis from ξn like {C–H–O}. Given this material basis, one can generate a list of potential elemental compositions from the Materials Project, for which we denote the materials library. One can then use η to select the candidate elemental compositions from the materials library. If following our previous example, the materials library may look like {CH3O, C2H6O, C3H6O, …} and our selected material may be C2H6O. For how ξn and η determine the atomic graph structure in the algorithm, please refer to Extension on Problem Formulation in the supplementary material.
2. Workflow automation
The algorithmic flow chart is visualized in Fig. 3 [Eqs. (1) and (2)]. The formulated basic optimization method is shown in the right block, while the materials generation and evaluation processes are laid out in the left blocks. In our approach, the algorithm is initialized and repeated based on different random seeds for the five attempts in Fig. 5. The hyperparameters for each algorithm are provided in the supplementary material and were selected from commonly used default values for each algorithm.
Algorithmic flow chart for the optimization framework, corresponding to Eqs. (1) and (2) and Fig. 2.
B. Optimization algorithms
1. Bayesian optimization
The core of machine learning (ML)-based optimization algorithms is to employ the explored materials data to fit a surrogate assisting the optimization process. Here, there are two main ways of adopting ML methods for design optimization: (1) building up a surrogate model of the design space for final exploration, namely the Bayesian optimization method; (2) utilizing a deep neural network (DNN) to learn the interactions of input elemental compositions and their corresponding material properties, which can be considered as the interaction of actions and their environment in which this DNN is considered as an “agent.” This method is known as deep reinforcement learning. The detailed mathematical formulation and derivations are shown in the Machine Learning Based Algorithms and Metaheuristic Optimization in the supplementary material.
2. Metaheuristic optimization
Metaheuristic optimization is a family of optimization techniques that are used to find the optimal solution by exploring a large search space. They are often inspired by natural processes or social behavior, such as animal behavior (e.g., ant, whale, and wolf), physical process (i.e., annealing of metals), biological processes (i.e., gene mutation), or fundamental mathematics.
Metaheuristic algorithms are typically applied to complex optimization problems where traditional methods are either not feasible or inefficient. They are highly flexible and can be adapted to a wide range of problem domains, in which elemental composition design is our target scenario. Metaheuristics are often considered as “black box” optimization methods, as they do not require knowledge of the problem structure. More specifically, they do not require the existence of the first-order differentiation form of the governing equations of the problem.
III. RESULTS
A. Elemental composition design benchmarks
Figure 4 shows the benchmarks of the 11 different inverse optimization algorithms (IOAs) and their objective properties explorations [column (a)], the distribution of the material evaluation numbers with the corresponding mean values [column (b)], and the distributions of scanned density with the overall mean values [column (c)], for both the single- (row 1) and multi-element (row 2) composition design cases. In the objective distributions, different IOAs are marked in different colors, and the corresponding optimization methods are abbreviated. Here, the materials count represents the total number of unique elemental compositions scanned by the different IOAs. The mean densities of objective spaces are computed to estimate the repetitive evaluations of the same elemental compositions by computing the mean value of the density scan space (see Sec. II and the supplementary material). Both metrics help us understand the “strategy” of different algorithms when they scan the chemical design space to search for the ideal elemental compositions.
Objective properties exploration, objective value evaluation, and targeted design space count for benchmarking different optimization methods for single- and multi-element composition design cases. (a1) Overall design space benchmarking for the 11 different optimization schemes used, marked in different colors for the single-element composition design (for the projections in 2D, please refer to Fig. S6 in the supplementary material). (b1) The distribution (for the five repeated experiments) and mean values of the overall material counts during the evaluation processes. (c1) The distribution (for the five repeated experiments) and mean values of the scanned material evaluation densities (see Sec. II and Density Scan in the supplementary material for the calculation of “mean density”) during the evaluation processes. (Related details can be found in the analysis of Fig. S7 in the supplementary material.) (a2) Overall design space benchmarking for the multi-element composition design (refer to Fig. S6 in the supplementary material). (b2) The distribution and mean values of the overall material counts. (c2) The distribution and mean values of the scanned material evaluation densities. The white triangular dots are the mean values for the 11 evaluations in both (b) and (c). The corresponding numbers are marked on top of the bar plot for both (b) and (c). BO, GA, PSO, GWO, HIWOA, ACO, DE, SA, BOA, AOA, and RUN denote the Bayesian optimization, genetic algorithm, particle swarm optimization, hybrid gray wolf optimization, hybrid-improved whale optimization algorithm, ant colony optimization, differential evolution, simulated annealing, bees optimization algorithm, arithmetic optimization algorithm, and Runge–Kutta optimization, respectively.
Objective properties exploration, objective value evaluation, and targeted design space count for benchmarking different optimization methods for single- and multi-element composition design cases. (a1) Overall design space benchmarking for the 11 different optimization schemes used, marked in different colors for the single-element composition design (for the projections in 2D, please refer to Fig. S6 in the supplementary material). (b1) The distribution (for the five repeated experiments) and mean values of the overall material counts during the evaluation processes. (c1) The distribution (for the five repeated experiments) and mean values of the scanned material evaluation densities (see Sec. II and Density Scan in the supplementary material for the calculation of “mean density”) during the evaluation processes. (Related details can be found in the analysis of Fig. S7 in the supplementary material.) (a2) Overall design space benchmarking for the multi-element composition design (refer to Fig. S6 in the supplementary material). (b2) The distribution and mean values of the overall material counts. (c2) The distribution and mean values of the scanned material evaluation densities. The white triangular dots are the mean values for the 11 evaluations in both (b) and (c). The corresponding numbers are marked on top of the bar plot for both (b) and (c). BO, GA, PSO, GWO, HIWOA, ACO, DE, SA, BOA, AOA, and RUN denote the Bayesian optimization, genetic algorithm, particle swarm optimization, hybrid gray wolf optimization, hybrid-improved whale optimization algorithm, ant colony optimization, differential evolution, simulated annealing, bees optimization algorithm, arithmetic optimization algorithm, and Runge–Kutta optimization, respectively.
There are five major findings from Fig. 4. (1) Comparing (a1) and (a2), the scanned objective spaces are smaller for GWO compared with other optimization methods for both the single- and multi-element composition design cases. For the single-element composition design case: (2) the evaluated elemental composition numbers are similar for most methods, whereas SA and AOA are observably smaller. Both GWO and SA display relatively high variances compared with other methods. (3) The mean scanned densities are relatively higher for the SA and AOA methods, with both displaying higher variances. Both SA and AOA possess higher variances. For the multi-element composition design cases: (4) RUN has the lowest evaluated elemental composition numbers. GWO and BWO have the highest variances. (5) For the mean scanned densities for multi-element compositions, GWO and RUN have the highest evaluated elemental composition numbers and variances.
From Fig. 5, we can draw five key observations. (1) When comparing (a1) and (a2), the mean objective distributions differ for the single-element and multi-element composition design cases. (2) For single-element design, GA, SA, and RUN display high mean objective values. SA and RUN also exhibit relatively low variances. GWO shows relatively high mean objectives but with higher variances, while BO demonstrates a stable performance with a low mean objective value and low variance. (3) DE and AOA exhibit repetitive scans in the same composition region, indicating a tendency to get “trapped” in a local space. This aligns with their relatively lower mean objectives in (a1). Moving on to multi-element design, RUN achieves the highest mean objective values with low variances. GA, GWO, and BWO show higher variances. (4) AOA and SA again exhibit repeated scanning of the same elemental compositions, suggesting a tendency for optimization methods to get trapped in local composition spaces. (5) When considering both single- and multi-element cases, different optimization algorithms demonstrate similar elemental composition search strategies, as evidenced by similar frequency distributions.
Statistical evaluations of the optimization processes. (a1) The distribution of the mean objective values (from five attempts) for the 11 evaluated optimization algorithms for single-element composition cases. Different colors mark different attempts, as indicated in the legend in the top-right corner of this figure. The white triangular dots stand for the mean values of the five attempts. The “mean objectives” are calculated as per optimization attempt (under the same random seed). (b1) The frequency distributions (i.e., repeated evaluations) using kernel fits of different objective values for the 11 different optimization algorithms (for details, see Sec. II and Additional Results in the supplementary material). Note that the original data distribution is shown in the top left corner of the inset. (a2) The distribution of the mean objective values for the 11 evaluated optimization algorithms for multi-element composition cases, where the white triangular dots are the mean values. (b2) The frequency distributions of different objective values for the 11 optimization algorithms.
Statistical evaluations of the optimization processes. (a1) The distribution of the mean objective values (from five attempts) for the 11 evaluated optimization algorithms for single-element composition cases. Different colors mark different attempts, as indicated in the legend in the top-right corner of this figure. The white triangular dots stand for the mean values of the five attempts. The “mean objectives” are calculated as per optimization attempt (under the same random seed). (b1) The frequency distributions (i.e., repeated evaluations) using kernel fits of different objective values for the 11 different optimization algorithms (for details, see Sec. II and Additional Results in the supplementary material). Note that the original data distribution is shown in the top left corner of the inset. (a2) The distribution of the mean objective values for the 11 evaluated optimization algorithms for multi-element composition cases, where the white triangular dots are the mean values. (b2) The frequency distributions of different objective values for the 11 optimization algorithms.
To further quantitatively support our observations on the uncertainties of different optimizations for repeated elemental composition design experiments in Figs. 4 and 5, Tables I and II show the variances among different experiments per different optimizations. Both Table I and Fig. 4 show that the materials count variances are very different for single- and multi-element compositions, due to the very limited single-element compositions existing in the chemistry table. For the multi-element composition design, BO has the lowest variance and GWO has the highest variance of materials count, with GWO’s variance ∼7398 times higher than BO. Other than that, GWO’s variance is ∼27 times higher than that of RUN of materials count. GA and AOA have approximately the same level of uncertainty for materials count. The relative magnitude of the variances of the materials count and the mean density scan does not follow the exact same trend. Yet for multi-element composition design, GWO still exhibits the highest variance, with 366.93% higher than that of BO and 9992.08% higher than that of BO. Generally, it can still be deduced from Table I that, for both single- and multi-element composition design, nature-inspired algorithms (i.e., GA, GWO, and BWO) display relatively higher variances, indicating higher uncertainties, which will be further discussed in Sec. IV. Table II further confirms the observation that nature-inspired algorithms present higher uncertainties. For both single- and multi-element composition design cases, GA, GWO, and BWO show evidently higher variances of the mean objectives than those of other optimizations. For the single-element composition design, GA presents 155.56% and 1093.8% relatively higher variances, GWO presents 35.24% and 531.77% relatively higher variances and BWO presents 95.20% and 811.86% relatively higher variances than those of RUN and BO, respectively. For the multi-element composition design, GA presents 2136.79% and 6380.87% relatively higher variances, GWO presents 2190.70% and 6537.04% relatively higher variances, and BWO presents 2132.71% and 6369.04% relatively higher variances than those of RUN and BO, respectively.
Variances for the five repeated experiments of different optimization algorithms for single- and multi-element composition design cases, expressed with single-digit precision. Note that MC denotes “materials count” and MDS stands for “mean density scan.” The four columns correspond to Figs. 4(b) and 4(c). For better presentation, we rescale the variances of MDS 100 times larger (indicated in the form ×10−2).
Optimization . | MC (single-element) . | MDS (single-element) . | MC (multi-element) . | MDS (multi-element) . |
---|---|---|---|---|
BO | 0.2 | 6.5 × 10−2 | 104.3 | 10.1 × 10−2 |
GA | 0.3 | 42.9 × 10−2 | 179 704.5 | 101.7 × 10−2 |
PSO | 1.7 | 28.9 × 10−2 | 12 575.7 | 7.1 × 10−2 |
GWO | 14.3 | 13.3 × 10−2 | 771 728.3 | 1019.3 × 10−2 |
BWO | 0.7 | 26.0 × 10−2 | 402 679.3 | 131.2 × 10−2 |
ACO | 0 | 9.7 × 10−2 | 427.3 | 15.8 × 10−2 |
DE | 0 | 79.6 × 10−2 | 73 394.7 | 56.8 × 10−2 |
SA | 45.2 | 3576.7 × 10−2 | 19 843.8 | 12.5 × 10−2 |
BOA | 0.2 | 8.4 × 10−2 | 29 907.2 | 5.2 × 10−2 |
AOA | 5.8 | 1156.6 × 10−2 | 150 198.5 | 44.1 × 10−2 |
RUN | 0 | 1098.7 × 10−2 | 27 895.5 | 218.3 × 10−2 |
Optimization . | MC (single-element) . | MDS (single-element) . | MC (multi-element) . | MDS (multi-element) . |
---|---|---|---|---|
BO | 0.2 | 6.5 × 10−2 | 104.3 | 10.1 × 10−2 |
GA | 0.3 | 42.9 × 10−2 | 179 704.5 | 101.7 × 10−2 |
PSO | 1.7 | 28.9 × 10−2 | 12 575.7 | 7.1 × 10−2 |
GWO | 14.3 | 13.3 × 10−2 | 771 728.3 | 1019.3 × 10−2 |
BWO | 0.7 | 26.0 × 10−2 | 402 679.3 | 131.2 × 10−2 |
ACO | 0 | 9.7 × 10−2 | 427.3 | 15.8 × 10−2 |
DE | 0 | 79.6 × 10−2 | 73 394.7 | 56.8 × 10−2 |
SA | 45.2 | 3576.7 × 10−2 | 19 843.8 | 12.5 × 10−2 |
BOA | 0.2 | 8.4 × 10−2 | 29 907.2 | 5.2 × 10−2 |
AOA | 5.8 | 1156.6 × 10−2 | 150 198.5 | 44.1 × 10−2 |
RUN | 0 | 1098.7 × 10−2 | 27 895.5 | 218.3 × 10−2 |
Mean objective values and variances for the five repeated experiments of different optimization algorithms for single- and multi-element composition design cases, expressed with single-digit precision. Note that MO denotes “mean objective” and Vr stands for “variances.” The four columns correspond to Fig. 5(a). For better presentation, we rescale the variances of MDS 100 times larger (indicated in the form ×10−2).
Optimization . | MO (single-element) . | Vr (single-element) . | MO (multi-element) . | Vr (multi-element) . |
---|---|---|---|---|
BO | 91.9 | 58.7 × 10−2 | 84.5 | 57.5 × 10−2 |
GA | 97.2 | 1097.1 × 10−2 | 91.0 | 3726.5 × 10−2 |
PSO | 93.0 | 90.7 × 10−2 | 84.0 | 8.8 × 10−2 |
GWO | 96.8 | 580.6 × 10−2 | 89.8 | 3816.3 × 10−2 |
BWO | 94.9 | 838.0 × 10−2 | 91.0 | 3719.7 × 10−2 |
ACO | 92.5 | 148.3 × 10−2 | 84.0 | 307.5 × 10−2 |
DE | 94.2 | 557.5 × 10−2 | 89.2 | 593.5 × 10−2 |
SA | 97.6 | 81.5 × 10−2 | 82.6 | 145.8 × 10−2 |
BOA | 93.9 | 404.5 × 10−2 | 87.3 | 314.1 × 10−2 |
AOA | 94.7 | 840.2 × 10−2 | 83.6 | 523.2 × 10−2 |
RUN | 96.7 | 429.3 × 10−2 | 92.8 | 166.6 × 10−2 |
Optimization . | MO (single-element) . | Vr (single-element) . | MO (multi-element) . | Vr (multi-element) . |
---|---|---|---|---|
BO | 91.9 | 58.7 × 10−2 | 84.5 | 57.5 × 10−2 |
GA | 97.2 | 1097.1 × 10−2 | 91.0 | 3726.5 × 10−2 |
PSO | 93.0 | 90.7 × 10−2 | 84.0 | 8.8 × 10−2 |
GWO | 96.8 | 580.6 × 10−2 | 89.8 | 3816.3 × 10−2 |
BWO | 94.9 | 838.0 × 10−2 | 91.0 | 3719.7 × 10−2 |
ACO | 92.5 | 148.3 × 10−2 | 84.0 | 307.5 × 10−2 |
DE | 94.2 | 557.5 × 10−2 | 89.2 | 593.5 × 10−2 |
SA | 97.6 | 81.5 × 10−2 | 82.6 | 145.8 × 10−2 |
BOA | 93.9 | 404.5 × 10−2 | 87.3 | 314.1 × 10−2 |
AOA | 94.7 | 840.2 × 10−2 | 83.6 | 523.2 × 10−2 |
RUN | 96.7 | 429.3 × 10−2 | 92.8 | 166.6 × 10−2 |
Extracted elemental composition evaluations for the optimization benchmarking. (a) The top ten elemental compositions extracted from the single-element design case. (b) The top ten elemental compositions extracted from the multi-element composition design case. (c) The top ten chemical compounds extracted from the multi-element composition design case. (d) The normalized frequency (see Additional Results in the supplementary material for the extended analysis for frequency evaluations of optimizations) for different extracted elemental compositions under the single-element composition design case. The top inset dots stand for the projection of the normalized frequency in the range [1, 2.5]. (e) The normalized frequency under the multi-element composition design case. The top inset dots stand for the projection of the normalized frequency in the range [10, 45].
Extracted elemental composition evaluations for the optimization benchmarking. (a) The top ten elemental compositions extracted from the single-element design case. (b) The top ten elemental compositions extracted from the multi-element composition design case. (c) The top ten chemical compounds extracted from the multi-element composition design case. (d) The normalized frequency (see Additional Results in the supplementary material for the extended analysis for frequency evaluations of optimizations) for different extracted elemental compositions under the single-element composition design case. The top inset dots stand for the projection of the normalized frequency in the range [1, 2.5]. (e) The normalized frequency under the multi-element composition design case. The top inset dots stand for the projection of the normalized frequency in the range [10, 45].
From Fig. 6, we can derive five main observations. (1) Comparing (a) and (b), we observe that the top evaluated elemental compositions are similar. (2) The majority of the top chemical compounds are binary alloys, consisting of elements present in the top single-element compositions, such as Cs, Ce, Ba, and others. (3) RUN demonstrates effective scanning strategies and consistently converges toward local optimal elemental compositions. (4) SA demonstrates similar search strategies as RUN, evidenced by the high normalized frequency along the materials list. (5) BO appears to struggle in identifying the “ideal” elemental compositions and tends to scan the same region repeatedly, indicating a potential trap in the local optimal area. From Figs. 4–6, we have also added qualitative comparison and suggested potential applications in Table III. Since the overall goal is to design materials digitally with low uncertainty and targeted desired properties quantified through objectives, therefore, low variance, and high mean objectives can be viewed as “advantages” of the algorithm. “Targeted materials selection” is recommended due to repeated materials sampling, and “large-scale exploration” is recommended based on more materials candidates sampled with decent objective values.
Assessment and conclusion from the comparison results for optimization methodologies benchmarked. Note that the (dis)advantages concluded from the previous results may be subjective. Note that the “low selectivity” for RUN indicates that the variance is relatively low compared with DE and AOA.
Optimization . | Advantages . | Disadvantages . | Potential use . |
---|---|---|---|
BO | Low variances | Low objective | Large-scale materials exploration |
GA | High objectives | High variances | Targeted materials selection |
PSO | Low variances | Low objective | Large-scale materials exploration |
GWO | High objective | High variances | Not recommended |
BWO | High objective | High variances | Not recommended |
ACO | Low variances | Low objective | Large-scale materials exploration |
DE | Decent objective | Relatively high | Not recommended |
variances | |||
SA | High selectivity | Problem-specific | Not recommended |
BOA | Decent objective | Low selectivity | Not recommended |
AOA | High selectivity | High variance | Targeted materials selection |
RUN | High objective | Low selectivity | Materials exploration and design |
and stable |
Optimization . | Advantages . | Disadvantages . | Potential use . |
---|---|---|---|
BO | Low variances | Low objective | Large-scale materials exploration |
GA | High objectives | High variances | Targeted materials selection |
PSO | Low variances | Low objective | Large-scale materials exploration |
GWO | High objective | High variances | Not recommended |
BWO | High objective | High variances | Not recommended |
ACO | Low variances | Low objective | Large-scale materials exploration |
DE | Decent objective | Relatively high | Not recommended |
variances | |||
SA | High selectivity | Problem-specific | Not recommended |
BOA | Decent objective | Low selectivity | Not recommended |
AOA | High selectivity | High variance | Targeted materials selection |
RUN | High objective | Low selectivity | Materials exploration and design |
and stable |
Among all the metrics we evaluated, it can be deduced that RUN achieved superior optimization results compared with the others (Figs. 4 and 5; Tables I and II). This out-performance may be due to several reasons: (1) RUN is free of hyperparameters. The update process of RUN does not depend on fine-tuned hyperparameters, making it more robust especially for highly nondeterministic situations like elemental composition design. (2) The updating strategies of RUN are very different from those of other methods. The RUN optimization searches for the next set of evaluations by mimicking the update strategies of Runge–Kutta’s discretization of differential equations (see the supplementary material), which are fundamentally different from other metaheuristic methods. (3) From (2), the enhanced solution quality (ESQ) mechanism is well suited for our formulated elemental composition discovery cases. In the authors’ original proposal,38 the ESQ mechanism is employed to avoid the local optimal solutions and increase the convergence speed. In our case, as the material evaluation numbers are fixed at 5000, the faster convergence induced by ESQ enables the robust performance of the RUN algorithm in our benchmark cases. Further comprehensive characteristics are discussed in Sec. IV.
IV. DISCUSSION
Based on the results presented earlier, our findings raise several intriguing questions, for which we proffer the following discussions. First, we observe variations in the performance of different optimization algorithms across the two elemental composition design scenarios, such as the contrasting results of material counts and the mean density of RUN and AOA [Figs. 4(b) and 4(c)]. One explanation is that the algorithm effectiveness is problem-specific. Moreover, the optimization methods may also be sensitive to the dimensionality of the problem such that changing from single-element to multi-element composition design induces the change of metrics evaluation shown in Figs. 4 and 5. In addition, in the case of single-element composition design, the limited number of available elemental compositions can impact the performance of certain algorithms, potentially influencing their search capabilities. An interesting observation pertains to the consistent success of the RUN algorithm in achieving high mean objective values in both problem scenarios, especially in the second case [Fig. 5(a)]. This achievement can be attributed to the fact that the sampling strategy employed by RUN is hyperparameter-free and adept at storing information, which proves advantageous in handling high-dimensional spaces.
Another intriguing finding is that BO exhibits lower mean objectives accompanied by low variances in both cases. This trend becomes more apparent when considering Figs. 5(b) and 6(e). A plausible explanation would be that BO is sampling a wider range of diverse elemental compositions during the optimization process. This behavior can be attributed to the uncertainty-based sampling strategy employed by BO, which facilitates the exploration of a diverse design space and aids in constructing a more accurate surrogate model. AOA tends to repeatedly sample similar elemental compositions over an extensive number of evaluations [Fig. 5(b)]. This behavior suggests that AOA may become trapped within a specific region of the design space, necessitating further investigation to understand the underlying causes and identify potential remedies. Finally, we explore the major disparity in the performance of SA between the two elemental composition discovery cases. A possible explanation for this discrepancy is that SA is sensitive to the dimension of the optimization problem, leading to its different performances tested from the designed metrics in Figs. 4 and 5.
Each optimization algorithm demonstrates unique characteristics and strategies that contribute to its specific outcomes. The low variances of the objectives in Fig. 5 could be explained by BO’s low dependence on hyperparameters. BO aims to approximate an accurate surrogate model to optimize the objective, resulting in small mean objective values. This strategy utilizes uncertainty-based sampling, which decreases the overall objective value. The low uncertainty and variance in the final results are attributed to the small number of hyperparameters and the large sampling base. In addition, BO exhibits high material counts, as observed in Fig. 4(b). GA heavily relies on hyperparameters, such as mutation, offspring generation, and randomization (for details, see the supplementary material). A strong dependence on hyperparameters leads to high uncertainty in the mean objective values for both cases [Figs. 5(a1) and 5(a2)]. Moreover, GA repeated the samples of the same elemental compositions once ideal properties are detected, resulting in high sampling density compared with BO, PSO, and other related methods [Figs. 5(b1) and (b2)]. Meanwhile, PSO exhibits a lower uncertainty and variance as it has fewer hyperparameters, contributing to its more stable performance, as depicted in Figs. 5(a1) and (a2).
GWO is similar to GA in terms of hyperparameter dependency, resulting in a high uncertainty and sampling density focused on specific objective regions (Fig. 5). Notably, GWO displays a higher uncertainty even in material counts and mean density, which differentiates it from GA (Fig. 4). Similarly, BWO exhibits a high uncertainty in mean objective values [Fig. 5(a)], which may also be attributed to its dependence on hyperparameters and very similar searching strategies of mimicking animal dynamical behavior. However, BWO displays more material counts and less sampling density scan [Figs. 4(b) and 4(c)], potentially attributed to the setting of the encircling prey parameter, which enhances the sampling stability. ACO employs continuous probability-based sampling, resulting in relatively less uncertainty in the mean objective values [Fig. 5(a)]. ACO has similar mean objective values as BO and PSO while exhibiting lower sampling density scans [Fig. 4(c1)], indicating that ACO avoids oversampling the same regime for the single-element composition design case. DE algorithmically operates as a variant of GA, yielding comparable mean objective results [Fig. 5(a)]. With fewer hyperparameters compared to GA, DE exhibits a lower uncertainty [Fig. 5(a)]. The evolutionary strategies of DE demonstrate similar material counts to GA, as shown in Fig. 4, and repeated sampling for favorable materials, as depicted in Fig. 6. The performance of SA is observed to be more problem-specific, as both the material counts and density scan vary per different design cases [Figs. 4(b), 4(c), and a5(a)]. SA exhibits a low uncertainty for the mean objective values [Fig. 5(a)] due to fewer hyperparameters and the implementation of a Markov chain for the acceptance–rejection sampling criteria (for details, see the supplementary material). BOA shares similarities with GWO and BWO in terms of high uncertainty in the mean objectives [Fig. 5(a)] and low material counts for the multi-element composition design [Fig. 4(b2)], as they use particle-based, animal behavior-inspired search criteria (for details, see the supplementary material). However, BOA exhibits a relatively lower variance [Fig. 5(a)], indicating a lower dependence on, and having fewer, hyperparameters.
AOA uses a distinct searching strategy based on mathematical symbolized operations. AOA exhibits a high uncertainty with relatively low mean objectives [Fig. 5(a)]. AOA tends to oversample the local regime, leading to low material counts [Fig. 5(b)]. This behavior can be attributed to a large number of hyperparameters and the overdependence on certain factors, such as MOA and MOP (for details, see Arithmetic Optimization Algorithm in the supplementary material), resulting in a scaling factor that underestimates previous locations. Finally, RUN achieves high objectives with a low uncertainty and avoids oversampling in specific objective regions (Fig. 5). The material count is sensitive to dimensionality, yet all dimensions yield successful search results, as judged by the high overall mean objectives [Fig. 5(a)]. RUN is characterized by being hyperparameter-free, incorporating numerous intermediate hyperparameters to increase the dimensionality of the iteration scheme, and relying on previous solutions for updates based on search history. Following our previous discussion on RUN’s ESQ search mechanism and being hyperparameter-free, we would like to examine why RUN had outperformed other algorithms by having high objective values and low uncertainties. RUN employs an update scheme that has four update parameters,44 each containing randomization coefficients45 and are connected. This randomization and “parameter-connection” could help explain the low uncertainty observed in our numerical experiments. For the high objective, RUN’s ability to avoid local optima as reported by Ahmadianfar et al.38 could be a possible explanation, in which they contend that the adaptive mechanism employed to update the parameter and the ESQ are the main mechanisms that assure a good transition from exploration to exploitation.
V. CONCLUSION
In this paper, we propose a framework that can incorporate most state-of-the-art optimization algorithms to perform multi-elemental composition discovery and design for crystals. Our framework parameterizes the input elemental graph structure to tailor the corresponding bulk modulus K and Fermi energy EFermi. The key goal is to conduct a comprehensive benchmark for many of the state-of-the-art optimization methods and provide further insights for future implementation of these algorithms for elemental composition design. We optimized both single- and multi-element composition discovery cases for maximized K and minimized EFermi by formulating a single objective optimization. We benchmarked 11 different inverse optimization algorithms by fixing the material evaluations. We found that GA, GWO, and BWO exhibit higher variances and BO and RUN display generally lower variances for materials count, density scan (Table I), and mean objectives (Table II). We further conclude that nature-inspired algorithms contain more uncertainties in elemental design cases, which can be attributed to the dependency on hyperparameters. RUN exhibits higher mean objectives, whereas BO displayed low mean objectives compared with other optimization methods. Combined with materials count and density scan, it is proposed that BO samples more elemental compositions aiming to approximate a more accurate surrogate of the design space and, hence, has lower mean objectives, yet RUN will repeatedly sample the discovered optimal elemental compositions for the successful search strategy [Figs. 5(b), 6(d), and 6(e)]. We also proffer detailed discussions on the exhibited results that correspond to each optimization’s characteristics.
SUPPLEMENTARY MATERIAL
See the supplementary material for the extension for our problem formulation for detailed explanations. We provide further details on hyperparameters and basic derivation for the 12 optimization algorithms employed (the 11 employed optimizations for repeated experiments with deep reinforcement learning). We proffer some additional results in supplementary to our proposition in the main article.
ACKNOWLEDGMENTS
J.Y. acknowledges the support from the U.S. National Science Foundation under Award Nos. CMMI-2038057, ITE-2236190, and EFMA-2223785 and the Cornell faculty startup grant. The authors also acknowledge the computational resources provided by the NSF Advanced Cyberinfrastructure Coordination Ecosystem: Services and Support (ACCESS) program under Grant No. BIO210063 and the computational resources provided by the G2 cluster from Cornell University.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
H.Z. and H.H. conceived the research. H.Z. designed and conducted the experiment(s). H.Z., H.H., and J.Y. analyzed the results. J.Y. acquired the funding and resources. All authors reviewed the manuscript.
Hanfeng Zhai: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Hongxia Hao: Conceptualization (equal); Investigation (supporting); Methodology (supporting); Supervision (equal); Writing – review & editing (equal). Jingjie Yeo: Conceptualization (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (equal); Resources (lead); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available at https://github.com/hanfengzhai/materials-discovery-benchmark/tree/master/inverse_opt.
REFERENCES
k1, k2, k3, k4 in Ref. 38.
rand1 and rand2 in Ref. 38.