Strength through defects: A novel Bayesian approach for the optimization of architected materials.

We use a previously unexplored Bayesian optimization framework, “evolutionary Monte Carlo sampling,” to systematically design the arrangement of defects in an architected microlattice to maximize its strain energy density before undergoing catastrophic failure. Our algorithm searches a design space with billions of 4 × 4 × 5 3D lattices, yet it finds the global optimum with only 250 cost function evaluations. Our optimum has a normalized strain energy density 12,464 times greater than its commonly studied defect-free counterpart. Traditional optimization is inefficient for this microlattice because (i) the design space has discrete, qualitative parameter states as input variables, (ii) the cost function is computationally expensive, and (iii) the design space is large. Our proposed framework is useful for architected materials and for many optimization problems in science and elucidates how defects can enhance the mechanical performance of architected materials.

INTRODUCTION

Background

The advances in modeling, fabrication, and testing of architected materials have promulgated their utility in engineering applications, such as ultralight (–), reconfigurable (), high–energy absorption materials (), and in bioimplants (). Architected materials that mimic crystal microstructures and share their enhanced mechanical performances have captured the attention of the scientific community (, –). In particular, lattice structures consisting of different unit cells or with strategically removed beam lattice members, defined as geometric defects (), can emulate mechanisms in crystal microstructures, e.g., dislocation motion and strain hardening mechanisms (), such that they have superior strength, energy dissipation, and controllable failure. Unfortunately, a systematic, algorithmic search for the optimal designs of lattices can be prohibitively expensive (), so historically, their designs have been driven by engineering ingenuity and intuition and further improved by iterative evolution with design approaches such as addition of beam members, controlled shell buckling, chemical effects, or hybrid structures (, –, –). Although there have been successful designs based on intuition (, –, , , , , ), intuition cannot be universally applied to other problems. Furthermore, there is no guarantee that a successful design found by intuition is the true global optimal solution, with the intuition merely showing the potential for improvement within the design space. In addition, the counterintuitive strategy of removing lattice members to enhance the mechanical performance (, ) has not been explored as thoroughly as the aforementioned design approaches. This study aims to show how the removal of beam members guided by optimization can provide a remarkable mechanical enhancement.

We propose to optimize the strain energy density of a 4 × 4 × 5 microlattice structure. Each site of the microlattice can choose between four different unit cells and can have two different orientations of attachment, yielding a large combinatorial design space. In addition to qualitative, discrete input variables, the cost function for evaluating strain energy density is extremely expensive. In previously reported works, architected defects were introduced as a design strategy (–). Although the improved designs were validated by fabrication and testing, few systematic approaches have been taken to optimize them. Successful designs of tailored materials have been found with systematic searches using data mining (), particle swarm and evolutionary algorithms (–), random search algorithms (, ), and topology optimization (, –). However, in all of these studies, the input variables of the design were continuous, and these algorithms are not applicable to our qualitative design space. The optimization method in some of these studies (, –) required tens of thousands of samples to be evaluated. Recently, deep learning (DL) techniques have also been used in material design (, ), but DL generally requires millions of samples or evaluations of cost functions. In our study, where the cost of evaluating the black-box function is extremely high, large number of evaluations would be intractable, rendering these algorithms impractical (, –).

Bayesian optimization (BO) has proven to work well for optimization problems with costly black-box functions by finding global optimum with minimum number of function calls (, ). However, in previous studies, BO has been used for architected materials for problems with continuous variables (, , ). We introduce a novel BO framework, which we term “evolutionary Monte Carlo sampling” (EMCS), to optimize the prohibitively expensive qualitative input design space of our architected microlattice structures. We also introduce “stochastic Monte Carlo” (SMC) acquisition function as part of the EMCS framework, which outperforms the existing acquisition functions used in literature for qualitative variables.

Although some implementations of BO have been developed to work with qualitative input variables (–), these require training a large set of hyperparameters and modifying the underlying kernel. Our EMCS algorithm can handle mixed variable problems and find global optimum with minimum number of function calls without modifying underlying kernel for BO by using a genetic algorithm (GA) to deal with mixed variable space. On the basis of standard benchmark tests against other common optimization algorithms used in literature, we show that our EMCS framework is more efficient in terms of minimum number of function calls compared to other commonly used optimization schemes for mixed variable problems. Applying the EMCS framework to our microlattice structure problem, our optimization yields four orders of magnitude greater normalized strain energy density compared to the existing microlattice structures in literature. These notable results highlight the efficiency and usefulness of our proposed design strategy for optimization of large and expensive mixed variable design spaces.

Problem setup

The objective of our optimization problem is to design an optimal, nonmonolithic microlattice consisting of discrete unit cells to maximize strain energy density. Figure 1 illustrates the design scheme used to obtain and validate the optimum structure.

Fig. 1.

Design space of the optimization and flowchart of work.

(A) Classification of inputs. Top: Four types of unit cells: (A) unblemished (no missing members), (B) defected (four missing members shown in cyan), (C) defected (12 missing members shown in purple), and (D) defected (16 missing members shown in red). Middle: Schematic of orientations. Bottom left: Schematic of the 2D cross section of the 3D lattice in the plane perpendicular to the extrusion, showing 4 × 4 unit cells. Bottom right: Symbolic representation of the lattice, where a letter F or E in the rightmost box indicates that all the unit cells are either aligned or rotated by 45°. In both orientations, the structure is perpendicular to the load (vertical) direction. Thus, the input has 16 dimensions with four possible values and 1 dimension with two possible values. (B) Flow chart beginning with the selection of 50 initial random microlattices (and 5 intuitive solutions; see text). The two-step iteration loop, shown within the box with broken lines, consists of (i) finite element analysis (FEA) evaluation of the critical buckling load P, of a microlattice structure and (ii) sequential optimization of the design space using our proposed EMCS framework. After the two-step iteration loop is exited, the optimal design is fabricated and tested, as illustrated with SEM images of the microlattice structure before loading and after the instigation of buckling. Our objective is optimizing the impedance of fracture due to the early commencement of buckling modes, illustrated in the helium ion microscopy (HIM) image, and thus increasing the strain energy density.

In this study, we use four different unit cells, henceforth described as states. Each of these four possible states, shown in Fig. 1A, is labeled by a letter A, B, C, or D. The octet truss unit cell () without any defects is labeled with A and, hereafter, is referred to as the unblemished structure. Unit cells B, C, and D have missing members. We selected them as inputs because a lattice made up with any combination of them has sufficient rigidity and will not collapse under its own weight, as it was shown in previous work that describes these explicit design principles and their benefits with respect to the buckling and densification response (–). In addition, we tested other types of unit cells to examine their structural integrity, and it was found that these four were the best candidates based on the required dimensions realized by the multiphoton lithography (MPL) apparatus. Despite the fact that including more states is feasible, it will also exponentially increase the design space. For four states, the design space is 8.58 × 109 (including another design variable that will be described next). All of these combinations, regardless of symmetries that decrease the number of physically independent geometries, determine the required computational cost. Because of the repercussions of dimensionality, increasing the number of states or variables will have a substantial effect on the required number of function evaluations for finding the global optimum. Furthermore, lattices made of these four cells will not have a disproportional reduction in structural stiffness (–). These unit cells, or states, are positioned in a 4 × 4 × 5 array. The 16 qualitative design variables (i.e., the lattice sites of the 16 unit cells) can each have four possible states. Note that there are five layers in the “extruded” direction of the microlattice, aiming to emulate the periodic arrangement of lattice layers in a crystal lattice structure (). That is, whatever 4 × 4 array of cells is selected (e.g., the structure shown at the bottom of Fig. 1A) is exactly repeated in all five layers of the extruded direction, discussed in further detail in the subsequent section. Apart from the strategy to imitate crystal lattices, arbitrary or irregular three-dimensional (3D) arrangement of unit cells will precipitate an abrupt increase of the design space, rendering the optimization problem intractable within normal computational or experimental budget limits. For instance, the 4 × 4 × 5 array will have ∼1.46 × 1048 structures if arbitrary 3D arrangements are used. This will also render the experimental validation of internal instability events abstruse using a characterization technique with a single detector as in the scanning electron microscopy (SEM). While optical microscopy from two objective lenses and digital image correlation could capture the 3D deformation, the lower resolution and mitigation of noise would increase the complexity and cost of the experimental validation. However, there is one additional important design variable, which is the relative orientation of the unit cells. Adjacent cells can be connected either at their faces (labeled as state F in Fig. 1) or their edges (state E). This two-valued state is the 17th variable of the problem. The reason for selecting only these two orientations is that other orientations would lead to a noncuboid distortion of the lattice, and allowing that distortion in the design space would significantly increase the complexity of our search. Figure 1A shows a compact way to symbolically represent the 17 inputs of a design: the unit cells in the 4 × 4 array and their orientation.

Figure 1B shows the optimization scheme of the problem. Ideally, as a load is applied to a lattice, our optimization goal would be to maximize the strain energy density of the lattice while maintaining the lattice’s structural integrity and stiffness before fracturing. To achieve this goal, a numerical code would be required to compute the static behavior of the lattice in its elastic domain and its elastic-plastic domain (where post-contact of lattice members induces densification) at the instigation of fracture-induced collapse. For the optimization of a 3D lattice with thousands of beam members, this analysis needs to be repeatedly carried out, which would be prohibitively expensive for any search method. Furthermore, this optimization goal presents difficulties with respect to experimental fabrication and verification. The material that was used was the hybrid organic-inorganic Zr-DMAEMA [(2-dimethylaminoethyl) methacrylate] (further details about the fabrication and the material properties are provided in Materials and Methods). Using other types of photoresists would require finding the correct fabrication parameters that will enable ductile instead of brittle behavior. Otherwise, the commencement of buckling would lead to early brittle fracture before densification. However, this material and design challenge should be surpassed for applications that require high-strength materials such as ceramics made through pyrolysis (). The reason is that the reduction of the unit cell size will reduce the distance between proximal beam members, enabling densification after buckling. Despite the fact that this would lead to ultrastrong and resilient materials, the volume reduction should also be included as a design variable to maximize the densification before fracture. Since the MPL polymerization process is a near-threshold phenomenon (–), we have observed that small variations in the process can affect the plastic performance (, ). Hence, the strain energy density at fracture collapse u and other post-yield constitutive parameters of the fabricated lattices, such as the back stress modulus, yield function, tensile flow resistance, and rate of strain hardening, are likely to be too sensitive for a valid comparison between experimental and numerical results. To mitigate the problems of the high computational cost of numerically optimizing u and the inherent sensitivities in the mechanical properties of the fabricated materials, we conceived another metric to optimize the structure: the structure’s critical buckling load P. Our motivation for choosing P as the proxy cost function to be minimized is that the large deformation of lattice members instigated by buckling instability leads to the densification and stiffening of the structure by the post-contact of the deformed lattice members. The consequence of densification is that fracture is impeded despite excessive plastic deformation. Using P as the proxy performance function to be minimized is advantageous because the calculation of P does not require excessive computation since P can be determined by solving an eigenvalue problem in the elastic domain. Moreover, P depends on parameters that characterize the deformation of the structure in the easier-to-compute elastic domain. Not only is the optimization with this proxy less expensive to compute, but the laboratory measurement of P is also not overly sensitive to fabrication variations. In previous studies, we found that the elastic modulus and Poisson’s ratio are relatively less sensitive to variations of the fabrication conditions (, ). Although the critical buckling load has been used as a cost function for machine learning optimization (–), it must be noted that these investigations required thousands of training data using continuous variables for single structural members, which has substantially smaller computational cost than a 3D lattice structure. Despite the fact that the challenge of an explicit finite element analysis (FEA) simulation is alleviated, a small number of data is still a necessary requirement.

Although the effect of buckling deformation on densification is well established (), there is not a succinct, closed expression that relates the critical buckling load P to the strain energy density at fracture-induced collapse u. In part, this is due to the fact that u depends on properties in the elastoplastic contact and fracture regime, while P depends on those in the elastic regime. However, we show below that laboratory experiments demonstrate that minimizing the critical buckling load leads to an exponential increase of the structure’s strain energy density; therefore, a small P is a very good proxy for large u. Details about the methodology to obtain the critical buckling load are provided in Materials and Methods.

The flow chart in Fig. 1B shows that our optimization method is initialized with 50 randomly chosen lattice structures. The main work needed to achieve optimization is carried out in a two-step iterative loop (shown by the steps within the dashed lines). The first (and most computationally expensive) step (i) is either the FEA calculation of the P of the 50 random structures (but only for the first iteration of the loop) or the calculation of the P of the new lattice design that was selected by the BO algorithm [step (ii)] in the previous iteration. The algorithm exits the iterative loop with an “optimally designed” lattice, i.e., with the globally minimized P. As shown in Fig. 1B [step (ii)], our 17D search space has many local minima. A key strength of our framework is that it does not become “stuck” at a local minimum but moves on to find the global minimum. Once the iteration loop for optimization is complete, the optimum lattice is then fabricated and tested.

The BO framework: EMCS

Optimization techniques, in general, are used for optimizing functions with ordered continuous, input states such as real numbers or integers. In our problem, an input state is a unit cell or an orientation with no inherent numerical value i.e., a qualitative variable. Even nonconvex optimization problems in which all of the inputs are quantitative variables can often be efficiently solved with gradient and other standard methods, but for our optimization problem, the qualitative inputs rule out these methods. The defining characteristics of our optimization problem is that step (i) of the iterative loop, where the cost function (i.e., P) is computed, is much more costly to evaluate than step (ii), where the search algorithm determines the next set of input values (i.e., the 17 values in Fig. 1A) to be evaluated. In optimizations such as ours, it is highly advantageous to choose a sequential search algorithm that converges to the optimal solution with the fewest possible iterations of the two-step loop, even if it requires a slight increase in the cost of step (ii).

We propose a nascent BO framework to be used in step (ii) that we term EMCS, which can handle mixed variables and uses all of the information about the design space from previous function evaluations. Running our BO algorithm on test problems with known global minimum values (see the next section) and with complexities and dimensions similar to our microlattice optimization problem, we determined that our BO would find the global minimum of P for a 17D problem (8.58 × 109 combinations) within 250 evaluations of the cost function. BO is efficient because (see Materials and Methods) it uses all of the information from the previous evaluations of the cost function to examine the design space and determine the next point to evaluate by using Bayes’ theorem. In the BO’s exploration phase, it examines regions of the design space where the cost function is unknown and in which the uncertainties of its value are very large. In the BO’s exploitation phase, the BO is focused in regions of the design space that are most likely to contain the global optimum. An acquisition function is used to determine the next point in design space to be evaluated by balancing both exploration and exploitation and is tailored so that the evaluation of the cost function at the point will yield the most useful information about the design space. The acquisition function prevents the BO from getting “stuck” at a local extremum by exploring, whereas other optimization methods, such as simulated annealing and GAs, often use less efficient probabilistic methods (e.g., heating/cooling or mutation) to get out of a local extremum. Conventional BO methods have been developed for problems with continuous and ordered inputs. The novelty of the present study is the adaptation of the BO to the discrete, qualitative inputs of the microlattice (Fig. 1A) using GA. Acquisition functions used in literature are usually designed for continuous variables and do not perform well for qualitative variables. We introduce SMC as part of our EMCS framework, which outperforms all the existing acquisition functions in literature for qualitative variables. Our major contributions for the optimization technique can be summarized as follows: (i) A novel EMCS framework is developed that uses a constrained GA to handle mixed variable problems without modifying the underlying kernel and is more efficient than existing optimization algorithms for mixed variables. (ii) A proposed acquisition function SMC is introduced as part of EMCS framework that outperforms existing acquisition functions for categorical variables in anisotropic input parameter spaces. (iii) Benchmark tests are performed on EMCS and SMC, and the EMCS framework is shown to outperform existing optimization strategies for expensive mixed variable problems. The details of our optimization algorithm are included in Materials and Methods.

For benchmarking the efficiency of EMCS, we use a broad range of existing optimization strategies that are commonly used for expensive black-box function optimization with mixed variable capabilities: GP_skopt [a Gaussian process (GP)–based BO library]; DV_skopt is a dummy variable random search done by uniform sampling within the given bounds and is an efficient methodology for hyperparameter optimization ()); RF_skopt (random forest sequential optimization) uses decision trees; GBRT_skopt is a sequential optimization that uses gradient boosted regression trees and is a set of mixed variable optimization libraries publicly available in scikit-optimize (); erlineTPE_Hyperopt (tree-structured Parzen estimator) is a sequential method for optimizing expensive black-box functions, introduced in () and available in the Hyperopt package; and SMAC3, which is a popular BO algorithm in combination with an aggressive racing mechanism to efficiently decide which of the two configurations performs better (). Note that GP_skopt uses a hedge strategy () as an acquisition function, with expected improvement (EI), probability of improvement (PI), and upper confidence bound (UCB) acquisition functions used in the hedge portfolio. RF_skopt and GBRT_skopt use UCB as their acquisition function. EI is used as an acquisition function for TPE_Hyperopt and SMAC3. The number of points to be sampled on the surrogate surface for scikit-optimize libraries was set to 105. Our EMCS framework uses our SMC acquisition function.

There are two important items that must be addressed. The primordial fundamental problem would be an infinite array with infinite numbers of defected states. However, from the standpoint of computational budget, such a problem would require hundreds of thousands of data. In addition, the buckling load of an infinite lattice would be related to long wavelength buckling, which leads to global instability of the structure and catastrophic collapse (). From this perspective, investigating a finite array is expedient since, in this case, the buckling load is related to microscopic instability patterns that are associated with densification and impedance of catastrophic fracture (). This array can be used as the building block to create a large structure having high energy absorption and densification of its constituents. On the basis of previous results (, , ), it was verified that the size of the array will not render the boundary effects dominant, as it will also be shown by the experimental results in the next section.

Moreover, as in most design problems, geometric and loading symmetries (in our case, reflection symmetries in each of the three principal directions and the fourfold rotational symmetry about the vertical axis) reduce the number of independent input configurations in the design space. Therefore from the perspective of the optimization algorithm, our design space is larger than 8 billion, but the space of independent configurations is fewer than a billion. However, as it was addressed in the previous section, the computational cost depends on all of the design space and not the independent configurations.

Compared to other optimization techniques, 250 evaluations for this 17D problem is very small. For example, our numerical experiments on searches with our test functions (see the Supplementary Materials for a list and definitions of those functions) showed that GAs (see the Supplementary Materials) require an order of magnitude more evaluations of the cost functions to find the global minimum than our BO required. Unlike the tests of the BO described in the next section, the search used in the actual optimization of the lattice was initiated not only with the P values of 50 random lattices but also with the P values of 5 lattices that we refer to as “intuitive solutions” that were selected, fabricated, and tested because studies by others (or our own intuition) suggested that they might be near optimal. As discussed in the next section, they were not even close to optimal. These were included to augment the initial dataset. For any black-box function, it is mathematically impossible to say whether the global optimum for any problem has been achieved. However, on the basis of the benchmarking presented in the next section, we believe that the microlattice design found by our BO is the global optimum of our design space.

RESULTS

EMCS framework benchmarking results

The EMCS algorithm was validated on a suite of test functions against alternative optimization strategies commonly used for expensive mixed variable black-box functions. We present the results of the tests on our “‘encrypted amalgamated function,” an anisotropic function that we created using a combination of commonly used benchmark functions (, ). We created this test function, described in detail in Materials and Methods, for use with mixed variables. The 17D variation of the function tested here consists of 17 categorical variables with four categories each. To incorporate the noise caused from fabrication imperfections and the experimental measurement, a 0.005 noise variance (7% SD) was built into the test function on the basis of previously reported results ().

The benchmarking algorithms were run for the test function in comparison with EMCS (with SMC as the acquisition function). The function evaluation budget was set at 300 with 50 initial randomly sampled data points, except for SMAC3, which determines its own initial sample size. The evaluation budget for SMAC3 was still capped at 300. Each benchmark was run 10 times for every test function. The sequential optimization history results for all algorithms applied to the optimization of the amalgamated function with 17 categorical variables (each with four possible input values so it is similar in scope to the architected lattice optimization problem) are presented in Fig. 2A. The figure shows how the normalized cost function (i.e., the difference between current optimum and global optimum) decreases with increasing evaluation number. The best algorithm is defined to be the one that consistently reaches a cost of zero after 300 evaluations. Note that only the normalize current best optimum is shown in Fig. 2A, and the algorithm is actually exploring the 17D space searching for new optima and the value in the graph only changes when a better optimum is found by the algorithm. The EMCS framework with SMC acquisition function performs significantly better than the other algorithms for the 17D, anisotropic modified amalgamated function. Moreover, for the encrypted amalgamated function 17D, EMCS(SMC) finds the global optimum (with zero cost) 9 of 10 times. None of the other algorithms are able to find the global optimum in any of their runs. RF_skopt and GP_skopt are the next best optimization algorithms, but they are significantly worse than our EMCS framework. This shows the efficacy of our algorithm for qualitative variables against other algorithms and the strength of the EMCS framework for expensive mixed variable problems.

Fig. 2.

Convergence of the BO algorithm.

(A) Optimization progress history for algorithm comparison. Each algorithm was run 10 times for encrypted amalgamated test function for a maximum budget of 300 function evaluations including 50 initial random samples. Cost ≡ global optimum − current optimum. The optimum value of the Cost is 0. For each algorithm, the figure shows the mean Cost of the 10 runs as a function of the number of test function evaluations. The cloud around the mean value shows ±0.2 SDs. Our EMCS (SMC) algorithm outperforms the others. (B) Convergence of the BO algorithm for the FEA-computed critical buckling load P of the microlattice plotted as in (A), but for the lattice, the cost function is not shifted or normalized. Evaluations 56 to 250 for the 195 configurations are determined by the BO algorithm (see text for the other 55). The insets in (B) show SEM images of the structures with four fabricated states: the random structure with the lowest P, the monolithic state D; the intuitive solution with the lowest P; and the BO optimum or optimal lattice found with the BO algorithm (see text for details of these structures).

Optimization of architected microlattice structures

Figure 2B shows the results of using our BO algorithm to optimize the lattice problem initialized with 50 random lattices. In addition, five microlattice structures that we call the intuitive solutions were also used in the initialization. These lattices were based on the proposed design principles for the strain energy density (, ).

These five structures have defected states in their ±45∘ planes, leading to localized collapse that imitates the failure of bulk materials at the planes of maximum shear (). From those five structures, the one with the lowest P, referred to as the “best intuitive solution,” was fabricated and tested, along with eight structures that were randomly generated. They were fabricated by the MPL process [see the Supplementary Material ()]. The 17 inputs of the fabricated materials and the numerical values of their P are given in tables S1 and S2 (). These structures were tested, and the average difference between experiments and simulations with respect to P was found to be equal to 7.4%, indicating a good match between experiments and simulations. Figure 2B shows results of the optimization algorithm. After 160 evaluations of the P [or 105 ≡ (160 − 55) BO iterations], our search found a structure with a critical buckling load 42.4% less than that of the best random structure. The optimal lattice found after 250 function evaluations was defined to be the “BO optimum.” The P of the BO optimum is 85.7% less than that of the defect-free unblemished structure, as shown in Table 1.

Table 1.

Experimentally obtained values of the critical buckling P and the experimentally measured value of the strain energy density at fracture and u.

The last column is the ratio of the normalized strain energy density of each structure compared to the unblemished structure.

Structure	Pc (μN)	uf (MJ m3)	(ufi/ubi)/(ufUB/ubUB)
Unblemished	3905.8	0.071	1
Random sampling	938.6	2.85	20.12
Monolithic D	1106.0	2.51	52.8
BO optimum	547.1	14.71	12,464

However, there are four other microlattice structures that we deemed important to evaluate. We conjectured that the lattice with the lowest P should be a configuration consisting of only one of the four units cells. The logic being that if one unit cell is better than all of the others, then a lattice made up of only those superior unit cells would be the optimal lattice. Therefore, we evaluated the P values of the A, B, C, and D unit cells (see Materials and Methods, Experimental Design for the simulations and mechanical modeling). Thus, we found a hierarchy of the P values of the unit cells. In addition, the P values of 4 × 4 monolithic microlattices made up exclusively of unit cells A, B, C, or D have the same ordering as the P values of the unit cells. The 4 × 4 microlattice with only D unit cells connected at their edges (17th design variable = E) also has lower critical buckling load compared to the microlattice with D unit cells connected at their faces (F) (see Materials and Methods, Experimental Design for the simulations and mechanical modeling). Hence, we speculated that the optimum lattice might intuitively only have D cells connected at their edges. We call this the “monolithic state D.” However, Fig. 2B shows that the P of the BO optimum is 43.4% less than that of the monolithic state D, which is opposed to our initial expectation based on mechanics principles only. The comparative results are shown in Table 1.

The 17 input values of the BO optimum lattice that uniquely determine its structure along with its P are listed in the first row of table S1 (). Figure 2B shows SEM images of the fabricated best random structure, the monolithic state D, and the BO optimum. Note that 17 input states of the BO optimum have no discernible pattern that would have led to its design by intuition. The reduced P and the unexpected structure of the BO optimum support not only the efficiency of BO but also, more generally, the argument against the use of intuition in designing materials.

Properties of the BO optimum

To assess the effect of the critical buckling load on the mechanical response and strain energy density of the microlattice structures, additional specimens were fabricated (14 in total). Apart from the structures mentioned in the previous section, i.e., the three best solutions determined by the BO algorithm, the unblemished structure composed only of state A and the monolithic structure with state D were also fabricated. The arrays of the 14 different geometries are given in table S1 (). The results of the measured buckling load P, strain energy density at fracture-induced collapse u, and measured stiffness S are summarized in table S2 (). As shown in table S2, the BO optimum has substantially lower critical buckling load compared to the other geometries, and there is a close match between the experiments and simulations.

Figure 3 shows illustrative results revealing the mechanical performance of microlattice structures. The stress-strain curves shown in Figs. 3A and 4 indicate that the design determined by the BO algorithm sustained fracture at a significantly larger deformation compared to the unblemished structure, the best intuitive structure, and the monolithic structure with state D. To evaluate the delaying of fracture-induced collapse, the deformation of the unblemished and the BO optimum structures were recorded and are presented in movies S1 and S2, respectively (). For the case of the unblemished structure (movie S1), shortly after the beginning of the loading, the first layer began to undergo buckling and almost instantly collapsed, resulting in the fracture-induced collapse event shown in the SEM image of the inset of Fig. 3A and indicated by the diamond symbol. These deformation modes are consistent with buckling eigenmodes of the unit cells and the whole structure obtained by the FEA simulations and explained in previous work (, ). As the indenter continued to load the structure, the next layer that carried the load also sustained fracture-induced collapse and the same mechanism occurred in the subsequent layer. Nevertheless, it is shown in movie S2 that the unblemished structure sustained buckling almost instantly after the initiation of loading (marked by the red box in Fig. 3A), resulting in contact of the beam members. This mechanism obstructed fracture-induced collapse of the layers due to the densification of the structure. The stress-strain curves of the optimum structures show serrated profiles, indicating the conglomeration of buckling and post-contact events prompted by densification, which is the same mechanism that has been investigated in our previous work (). In the light of movie S2, the layers of the structure sustained plastic deformation and fracture commenced only at 50% of deformation, which is verified by the video recording. The deformed structure just before the instigation of fracture is shown in the inset of Fig. 3A. As it will be shown later, this was also verified by high-depth-of-focus imaging of the tested specimens. While these responses can be efficiently characterized for a limited number of experiments, such a complex 3D mechanical response beyond the plastic domain is extremely computationally expensive to be attained with FEA simulations and then compared with experimental results obtained at the microscale for the reasons interpreted earlier. The effectiveness of the choice of the proxy cost function is illustrated in Fig. 3B, where the normalized strain energy density u/u is plotted against the critical buckling load P, where u is the strain energy density at the instigation of buckling. This normalization is used to enable focusing on the effect of the critical buckling load on the mechanical response after the occurrence of buckling. The normalized strain energy density of the optimum structure is 12,464 times higher than that of the unblemished structure (Table 1). The ratio of the normalized strain energy density of the other structures with respect to the BO optimum is given in table S2 (). Even without normalization, the strain energy density of the optimum structure is 203 times higher than that of the unblemished structure, one of the most thoroughly investigated and used architected geometries (–, –, –, , , , ). The strain energy density corresponding to the three best solutions determined by the BO algorithm shows a high sensitivity to the critical buckling load. While the mechanics responsible for the observed direct correlation between u/u and P with respect to the selected states requires further investigation, Fig. 3B shows that the selected cost function yields a solution that is orders of magnitude better than that of the unblemished structure. It is noted that in previously reported results relying solely on mechanics principles, the best result was inferior by just one or two orders of magnitude (, ). Moreover, although the random structures and the rest of the defined designs also exhibit one or two orders of magnitude higher u/u than the unblemished structure, they do not approach the optimum structure because they do not have such a low P. In comparison with previous studies that did not use optimization (, –), this result illuminates the unique effectiveness of our BO algorithm. Nevertheless, it is also important to elucidate these results in the light of the perspective deformation profiles of the tested specimens.

Fig. 3.

Mechanical performance of microlattice structures.

The commencement of buckling is indicated by a square for the BO optimum (red) and a diamond for the unblemished structure (black). (A) The prematurely instigated buckling in the optimum structure found by the BO algorithm precipitated an “ävalanche" of buckling events, leading to densification of the structure and the impedance of fracture (shown in the SEM image). Because buckling occurred at a later deformation stage of the unblemished structure, fracture was not obstructed, resulting in early collapse of the first layer (shown in the SEM image) and lower strain energy density before fracture. Black scale bars, 10 μm. (B) Experimentally determined critical buckling load P versus strain energy density at the instant of fracture collapse u normalized by the strain energy density at the onset of buckling u. The results show that 1/P is a good proxy function for u; as P decreases, u increases. With only 14 nonrandomly selected points, this figure should not be used to infer to correlations between u/u and P. The strain energy density of the optimum structures is 12,464 times that of the unblemished structure. Inset: Schematic to determine the onset of buckling. The critical buckling load P at (B) is calculated upon the commencement of an instability and is consistent with the onset of deformation in the videos included in the Supplementary Materials. With the progression of deformation, fracture of the layers or conglomeration of beam members occurred, leading to the collapse of the structure (F). The blue, black-lined (red-lined) shaded area represents the strain energy density u at the onset of buckling (at fracture u). In the schematic, u and u are commensurate for illustrative purposes.

Fig. 4.

Comparison of the mechanical performance of different structures.

The stress versus strain of the BO optimum (# 1 on table S1), BO penultimate optimum (# 2), BO antepenultimate optimum (# 3), unblemished structure (# 14), best intuitive solution (# 8), monolithic state D (# 5), and representative random structure (# 11).

Figure 5 shows high-resolution images of undeformed and deformed microlattice structures obtained with helium ion microscopy (HIM). HIM was selected since it provides high resolution and high depth of focus that cannot be accomplished by other imaging techniques, enabling the effective characterization of the internal members of the structure. The unblemished structure (Fig. 5A) sustained excessive collapse of the top layer (Fig. 5B), which propagated and proliferated in the next layers and the internal beam members (Fig. 5C). However, the optimum structure (Fig. 5D) demonstrated excessive buckling rather than catastrophic collapse and almost no fracture (Fig. 5E). The deformed buckling mode shows the out-of-plane buckling mechanism of beam members that was observed in movie S1 (Fig. 5F), without the total collapse of the layers, conversely to the unblemished structure. Figure 4G shows a side view of the optimum structure before testing. The large deformation of the array induced structure densification, which, in turn, inhibited fracture, leading to partial deformation recovery of the array upon unloading (Fig. 5H). A comparison of the mechanical responses of the unblemished and optimum structures reveals the remarkably improved mechanical performance determined by the BO algorithm based on the minimum critical buckling load. It was observed that the monolithic structure with state D also sustained premature fracture due to the delay of the instigation of buckling (Fig. 6). A random orientation of defected states does not provide the same response as the optimum structure a priori because specimens with random geometries, such as that shown in movie S3, sustained catastrophic collapse at the early stage of deformation. In this structure, buckling occurs at a higher load and does not provide early densification and impedance of fracture. Last, the best intuitive structure (movie S4) exhibited localized failure (Fig. 7), in agreement with our previously reported results (, ). However, this localized mechanism did not lead to uniform densification as in the BO optimum structure but in densification of a tilted layer (Fig. 7B), and fracture could not be impeded for the whole geometry. It is also noted that the stiffness of the optimum structure is 18% larger than that of the unblemished structure (table S2). The stiffness enhancement is attributed to the design of the states, as dictated by the previous work that guided the incorporation of the architected defects in the present study (, ). The results regarding the unblemished structure, monolithic state D, the BO optimum, and the best random structure are summarized in Table 1.

Fig. 5.

HIM images of the loaded and unloaded unblemished and optimum structures.

(A) Image of the unblemished structure consisting only of units cells of type A as in Fig. 1. (B) Same as (A) but after loading, showing severe fracture and collapse of many beam members. (C) High-depth-of-focus image of the region inside the square box shown in (B) revealing several fractured beams and the internal collapse of the upper layer that subsequently instigated the accumulation of damage in the underlying layers. (D) Same as in (A) but for the unloaded optimum structure. (E) Same as in (B) but after the structure was subjected to the same maximum compressive load as the structure shown in (B). Unloading of the optimum structure showed only excessive plastic deformation without catastrophic collapse and the manifestation of the buckling mode. (F) High-depth-of-focus image of the region inside the square box shown in (E) revealing the effect of buckling that led to deformation but no fracture due to the occurrence of densification. (G) Side view of the unloaded optimum structure shown from an isometric view. (H) Side view of the unloaded optimum structure shown from an isometric view revealing that fracture was inhibited throughout the structure due to the densification precipitated by the low critical buckling load. Scale bars, (A to H) 10 μm.

Fig. 6.

HIM images of the microlattice structure with a monolithic state D.

(A) Unloaded structure, (B) loaded structure, and (C) high magnification of part of the panel shown in (B) revealing fracture at the nodes of beam members (shown by white arrows). White scale bars, (A to C) 10 μm.

Fig. 7.

SEM images of the microlattice structure with the best intuitive solution.

(A) Before loading, (B) after loading showing excessive deformation in a slopped plane that eventually led to the collapse of the beam members on that plane, and (C) at maximum deformation showing sustained global collapse in the proximal region of the sloped plane. Black scale bars, (A to C) 10 μm.

DISCUSSION

We have presented a BO framework EMCS that can be used for problems with discrete, qualitative design variables and a large input design space composed of billions of possible combinations. Using different test functions, we showed that our BO algorithm would lead to the global optimum design of our 3D geometrically defected microlattice with only 250 function calls (i.e., FEA calculations of the critical buckling load). This is an order of magnitude smaller number of function calls than that required by a GA and two orders of magnitude smaller than that required in previously reported optimizations of architected materials (, –). To apply BO to discrete and qualitative variables, a new BO framework EMCS was introduced, which restricts the variable space for mixed variable problems using GA and a new acquisition function SMC was proposed as part of the EMCS framework. This framework was shown to outperform all existing optimization techniques for expensive mixed variable problems.

Rather than trying to directly maximize the normalized critical strain energy density u/u to obtain our optimal design, we chose to minimize the critical buckling load P as a less expensive and easier to experimentally validate proxy cost function. The optimal microlattice that we obtained using this proxy has a normalized strain energy density that is four orders of magnitude greater than the unblemished microlattice structure. Our use of the critical buckling load P as a proxy during the optimization revealed an unexpected strong inverse correlation between P and u/u. While a qualitative inverse relation has been previously established (, , , , ), we have shown that near the maximum value of u/u, it increases much faster than 1/P increases. These findings illuminate the design requirement for the emergence of global microscale buckling events for uniform densification instead of the previously reported layer-by-layer or tilted layer collapse. Furthermore, our results illuminate a mathematical methodology to design optimal architected materials with exceptional mechanical properties, and they pave the way on how the BO can be critical for low budget optimization of large, multivariable problems. The BO technique presented here is suitable for use in not only other architected materials but also a wide range of problems in science and engineering where other optimization techniques either fail or are too costly.

Our study has shown that a systematic optimization can produce large increases in the mechanical performance of architected structures and that the optimum design found by an unbiased BO algorithm has no obvious or intuitive pattern or structure. We are currently applying our BO technique to achieve multi-objective optimization of a nonmonolithic structure for both stiffness and strain energy density and to explore multi-objective optimization of a nonmonolithic structure for a structure with minimal anisotropy and maximal auxeticity, which has applications ranging from stress shielding on bone implants to ultralight resilient structures in aerospace vehicles against space debris. Specifically, this algorithmic scheme could also be used for the optimization of irregular metamaterials (), structures that have a nonuniform random arrangement. Since the algorithm presented in this work can also be used for problems with mixed variables, including continuous variables such as the prebending curvature of specific beam, members for some states would lead to even more complex buckling modes and reversible nonlinear elastic behavior, as it was presented in previous work (). The computational cost for nonlinear wave mechanics for cloaking, plastic evolution, instabilities, and crack propagation impedance is prohibitively expensive, rendering the utility of thousands of data for GAs or machine learning unfeasible. Moreover, including the constitutive material behavior and unit cell size as design variables into the optimization process will lead to the expansion of this design scheme for any category of material, setting the framework for the study of brittle to ductile transitions (), scalability with smaller unit cells but larger arrays, and nonlinear response (). Therefore, using BO with either discrete qualitative variables or mixed variables would circumvent this major challenge to obtain the global optimum for large-scale problems using small-scale and cheap experimental or simulation tools. In addition, future work should also focus on the enlargement of the array size and the correlation between long wavelength buckling and microscale buckling for the transition from densification to catastrophic collapse and its effect on the required computational cost. This will also be an important stage for the design of “polycrystalline” lattice structures. All of these compelling applications would also set the framework for further investigation on strategies that can decrease the computational cost for larger or mixed problems.

The EMCS framework is further developed to deal with higher dimensional problems, which is a concern with GP-based methods. Multi-objective optimization is an active research area and, generally, evolutionary algorithms are preferred. However, the inherent framework of EMCS using GA makes multi-objective optimization and batch update at each iteration feasible tasks and a natural fit for EMCS, a possibility that we are currently exploring and applying to develop nonmonolithic structures with maximal isotropy and auxeticity. The results presented in this paper have been using SMC as our acquisition function. Studies have shown that “Hedge”-based strategies (acquisition function strategies working with a portfolio of acquisition functions at each iteration) outperform individual acquisition functions (). A hedge strategy for mixed variables thus, with SMC included as a candidate in the portfolio of acquisition function, is an attractive prospect, and we are further developing our framework to incorporate this.

On the basis of the large improvements that we achieved with systematic optimization of microlattice structures, we speculate that a similar systematic optimization would also lead to large increases in the electrical, magnetic, wetting, surface chemistry, and thermal properties of architected materials. We are already using modified versions of our BO scheme in a wide range of science and engineering problems. In particular, we are currently applying our EMCS framework to maximize the power output of vertical axis wind turbines () and for draft tubes of hydroelectric turbines. In addition, we are modifying our current BO algorithm for use in planetary sciences, specifically to determine the structures of Jupiter’s Great Red Spot and zonal jet streams (, ). We believe that many optimization problems that require discrete and qualitative input design variables, which have previously been considered intractable because of the computational expense, can now be solved with our algorithm.

MATERIALS AND METHODS

Experimental design

All of the FEA simulations were performed with the multipurpose finite element code ANSYS (Workbench 18.0). The selected material properties are typical of polymeric materials used in MPL processing and have been found to yield a good agreement between experimental and simulation results. Specifically, the following beam lattice properties were used in the FEA: 1.281 elastic modulus, 0.4999 Poisson’s ratio, and 131.99 yield strength. The structures were discretized with 10-node, tetrahedral finite elements. Each structure comprised an average of 262,008 elements with 552,042 nodes. The number of nodes and elements was determined by examining the convergence of the critical buckling load with the variation of the mesh size. The boundary conditions were chosen to closely match those in the experiments reported in our previous work (, ). The compressive load was applied at the cusps of the top-face beams that initially came into contact with the indentation tip. Because the bottom face of each structure was firmly attached to the substrate, all the degrees of freedom of the bottom nodes were fully constrained. All of the structures were designed with the ANSYS Design Modeler.

The critical buckling load was determined from an analysis of the structure’s stability using the potential energy Π, defined bywhere W is the energy function of the system, u is the displacement field, V is the volume domain of the structure, and Πext is the potential due to the external loading. The structure instability was determined by setting the variational potential energy equal to zero, i.e.

The above equation yielded a set of differential equations that dictate the buckling response of the system, depending on the beam model used, e.g., Timoshenko, Euler-Bernoulli, or Elastica. Because of the complexity of the problem, the displacement field of the beam members was expressed bywhere N is the shape function in Galerkin theory and is the displacement field at the nodes of the mesh. Using Eqs. 2 and 3 and the fact that the critical buckling load P can be extracted from the potential due to the external loading, the variational potential energy yieldedwhere Kmat is the material stiffness that depends on the energy function of the structure and Kgeom is the geometric stiffness that depends on the loading conditions of the structure and is derived from the potential due to the external loading Πext. Both matrices depend on the shape function N used to approximate the displacement field in the beam members. For this equation to hold, the global stiffness of the structure must be positive semidefinite, leading to the instigation of buckling instability. Since the term multiplying the variational must be zero, Eq. 4 leads to the following eigenvalue problem

Equation 5 gives the critical buckling load P, whereas Eq. 4 provides the corresponding buckling mode . Further details about the analytical derivation of the foregoing equation and the association of the stiffness matrices with the shape function can be found elsewhere (). The numerical computation of the critical buckling load can be implemented in ANSYS Workbench using the Lanczos algorithm (Block Lanczos). Further details about this numerical process can be found elsewhere ().

The P values of the A, B, C, and D unit cells are 338.23, 290.72, 131.68, and 32.59 μN, respectively. Moreover, a lattice structure composed only of monolithic state D has a critical buckling load equal to 964.68 and 1130.32 μN for connection at the edges (E) or at the faces (F), respectively.

To calculate the normalized strain energy density ratio between different structures, the strain energy density at the instigation of buckling u was first calculated from the following equationwhere S is the stiffness of the structure and A is the contact area due to the applied load. The experimental values of u were verified by numerical results of the strain energy density at buckling. Since all of the structures comprised the same number of unit cells and contact with the structure commenced at the cusps of the top beam members, the contact area was the same for all structures. Therefore, using Eq. 6, the normalized strain energy density ratio between the BO optimum (with index 1) and any other geometry from table S2 (with index i) was obtained from the following equation

Fabrication

The microlattice structures were fabricated with a hybrid organic-inorganic material Zr-DMAEMA [30 weight % (wt %)]. The composition of this material is 70 wt % zirconium propoxide and 10 wt % DMAEMA (Sigma-Aldrich). Further information about the material preparation can be found elsewhere (, ).

The test structures were fabricated by diffusion-assisted high-resolution direct femtosecond laser writing, which uses MPL and the aforementioned photoresist for high-resolution fabrication. The system consists of a FemtoFiber pro NIR laser with a wavelength of 780 nm, pulse width of 100 fs, and repetition rate of 80 MHz. The local photopolymerization of the photosensitive material was accomplished with a 100× microscope objective lens (Plan-Apochromat 100×/1.40 Oil M27, Zeiss).

The laser output energy for the fabrication was measured right before the objective lens opening at 10 mW, and the scanning speed used was 80 μm/s. A detailed description of the setup can be found elsewhere ().

Mechanical testing and characterization

Uniaxial compression tests were performed in situ with a nanoindentation apparatus (PI 85 SEM PicoIndenter, Hysitron) mounted inside the chamber of a SEM (FEI Quanta 3D FEG). All of the structures were aligned such that the front faces to be visible during testing to facilitate the detection of the location and instant of buckling and fracture. To capture the structure collapse and track the proliferation of damage, each structure was deformed at a rate of 500 nm/s to a maximum compressive strain of 55%. To ensure repeatability, at least five tests were performed with structures of a given design. To capture the instigation of the first buckling instability and, consequently, compute the critical buckling load, the measured force-displacement curves were frame-by-frame juxtaposed with the recordings of the deformation to identify whether excessive deformation produced a singularity in the force-displacement response. Further details can be found elsewhere (, , ).

Statistical analysis

The BO that we used to determine the optimal structure was based on the EMCS algorithm that we developed specifically for optimization problems with discrete, qualitative variables. The fact that our optimization of an architected material used a design space with nonquantitative inputs (i.e., four types of unit cells) and required an expensive cost function (FEA) presented challenges that cannot be addressed with conventional optimization techniques. Therefore, we used a version of BO, EMCS, that is adapted to work with mixed variables using a GA to restrict the design space. We use “Hamming distance” in the covariance kernel and a novel acquisition function that we call SMC for qualitative variables.

BO is a sequential optimization technique that uses Bayes’ theorem (, , ). We are given a cost function f(x) (in this case the FEA-calculated P) to be optimized such that f(x) can be evaluated at any input value x (in this case the 17 input values in Fig. 1). Assume f(x) has been evaluated at n inputs, and let those inputs and their evaluations be represented as ordered pairs as D1: = {x, f(x1:)}. We use Bayes’ theorem to express

Here, P(g) is the prior distribution, a probability that represents our existing beliefs about a function g(x) from the space of possible candidate functions that approximate f(x); P(D1:∣g) is the conditional probability or likelihood of observing D1: given a specific candidate function g; and P(g∣D1:) is the conditional probability or posterior distribution, which represents the probability that a candidate function g is the cost function f, given our observations D1: ().

To minimize the number of evaluations of f(x) required to find the optimum solution, an acquisition function () is used to determine the next input x to be evaluated. After the next input, x is determined by the acquisition function, f(x) is evaluated, the posterior distribution is updated, and the process is repeated until the global optimum is found or a computational limitis reached.

Here, a uniform GP prior is used (, ). A GP is defined as a stochastic process such that a linear combination of a finite set of the random variables is a multivariate Gaussian distribution. The GP is uniquely specified by the mean μ(x) of g(x) and the covariance function k(x, x′) and is a distribution over functions. g(x) is a function that is sampled from this GP such that

Here, k(x, x′) is the covariance between any two input variables, x and x′ (, ). Once a GP has been defined, at any x, the GP returns the mean μ(x) and variance σ(x). In our problem, the inputs x are discrete and qualitative, rather than continuous and ordered due to the fact that the inputs are states (e.g., A, B, C, and D) that cannot be numerically ordered. We define the covariance k with a modified squared exponential kernel that works with inputs that are stateswhere ∣x, x′∣ is our capped distance vector made from two qualitative input variables. Hamming distance is used to determine this vector. The matrix is the covariance hyperparameter matrix. In this problem, we have chosen it to be diagonal with

The 17 values of h and ϵ are hyperparameters. Hyperparameter values essentially control the shape of the surrogate function based on the available samples of the objective function. The capped norm is the product . There were many possible choices of norms that work with qualitative inputs, but the simple capped norm in Eq. 10 worked well with all five of our test functions.

In this study, we use leave out one cross-validation (LOO-CV) (, ), a data-driven technique of model fitting given a representative sample of data. Cross-validation (CV) () is a commonly used hyperparameter fitting method. LOO-CV is the most computationally expensive variant of CV and therefore often prohibitive in large optimization problems, but it provides the best model fit. However, we need only 250 evaluations of our cost function, so we can use LOO-CV in this problem.

BO uses acquisition function, 𝒜, for exploring the surrogate model to maximize reward by balancing exploration and exploitation. The input point where this “reward” is optimized is the next test point or black-box function call for BO. For mixed variables, optimizing 𝒜 is problematic if the surrogate surface is not smooth. GAs, however, can be constrained to work with mixed variable spaces, particularly with qualitative variables. These variables can be dealt with by using probabilistic mutation rates. The mutated genes are only allowed mutation within the prescribed categories. The same approach for GAs for mixed variables is used in most openly available evolutionary algorithm libraries. This allows us to use our GA to optimize 𝒜(μ(x), σ(x)) with x constrained to discrete qualitative space. Note that this framework, constrained GA for optimizing 𝒜, coupled with a modified distance metric can now deal with mixed variable spaces. This is the basis of our EMCS framework.

Any acquisition function can be used within this framework for BO. However, most common acquisition functions such as EI, PI, and UCB () or Hedge techniques () do not work well with ordinal and particularly categorical variables. We propose a new acquisition function, SMC, which, for minimization of an objective, is defined aswhere r is a random uniformly distributed number between [0, −2] and μ(x) and σ(x) are the mean and SD returned by the GP at x, respectively. This is equivalent to taking Monte Carlo samples from left half of the distribution. For qualitative variable problems, this acquisition function outperforms other acquistion functions.

The use of a constrained GA for optimizing 𝒜, coupled with SMC acquisition function in our EMCS framework, provides a novel strategy to handle mixed variables. Using SMC or adding it to the portfolio of hedge strategies () can significantly improve algorithm performance for mixed variable problems. Existing BO frameworks for continuous variables can thus easily be modified, to incorporate EMCS framework, to increase capability for mixed variable problems without having to compromise efficiency.

The next input x to be sampled and evaluated by the cost function (FEA) is the input that minimizes 𝒜(x). If x were continuous, then 𝒜(x) could be minimized using standard schemes, such as gradient methods. However, when x is discrete or (as in our case) discrete and qualitative, these schemes cannot be used. Therefore, we minimize 𝒜 with a GA despite our earlier criticisms of GA. The reason that we can use a GA efficiently to minimize 𝒜 is that the evaluation of 𝒜, unlike the evaluation of P with an FEA, is inexpensive, and we can easily afford to perform thousands of evaluations of the GP to find the x that minimizes 𝒜.

Once the next input x has been determined by 𝒜, the P of that lattice is computed and the GP is updated. The whole sequential procedure is then reapplied to find the next input x. This process is repeated until the black-box function is globally optimized or maximum allowed budget is reached.

Our BO algorithm was benchmarked against other algorithm, which are commonly used to optimize expensive qualitative variable problems and shown to outperform them all (Fig. 2A). We also tested our algorithm on different test functions, and our algorithm found the global optimum of each of the test functions with only 250 function calls per test, including the 50 initial random inputs. The results are demonstrated in Fig. 8.

Fig. 8.

Convergence tests of our BO algorithm with a variety of test functions.

The colored curves show the convergence of the normalized and shifted amalgamated function, the encrypted amalgamated function, the Rastringin function, the Syblinski-Tang function, and the spherical test function. [See their definitions in the Supplementary Materials ().] Each test function has 17 dimensions, and each dimension has four discrete inputs. Each of the test functions was shifted by adding a constant to it so that its global minimum in the discrete search space of input values is zero. The shifted test functions were then evaluated with 50 random configurations, and the minimum cost function value was defined as the norm. Each test function’s value was divided by its norm, so its shifted and normalized "cost" was unity after the 50th evaluation. The cost function of the optimized configuration is zero.

Note that our algorithm successfully found the global minimum of each test function with only 250 or fewer evaluations of the test function. The test functions used for benchmarking are detailed in the Supplementary Materials.