Bayesian Optimization of Equilibrium States in Elastomeric Beams

Hyperelastic lattice-based architectures have demonstrated potential for energy absorption, multi-resonant vibration isolation, and bandgap controls, due to their inherent buckling instabilities [1–6]. The bistable structures have received specific attention for applications in vibration isolation and control, because of their negative stiffness, strong energetic response, and high adaptability [7–9]. Mechanical constraints play a critical role in the accessibility of bistable behavior, as demonstrated in a previous study on prismatic, inclined beams [1]. Bistable states are advantageous for lattice network reconfiguration, as the global static and dynamic properties of the network can be tuned by the position, energy, and spatial distribution of the local equilibrium states. However, the strength and existence of the second stable configuration are highly dependent on the reference beam geometry. For instance, Shan et al. demonstrated a transition between monostable and bistable beam buckling by tuning the aspect ratio and angle of prismatic inclined beams [1]. The bifurcation between mono- and bistable designs creates discontinuities in the design space, resulting in undefined objective values and inefficiencies in the optimization process. Our study investigates these challenges in the untapped design space of non-prismatic, inclined beams, developing a penalty method to regularize between stability domains and compiling a library of modular, beam geometries for future integration into architected lattices.

The design problem of our study focuses on the optimization of a single beam element, which is a tilted, slender beam with a variable thickness cross section defined with clamped-roller boundary condition. We consider two objectives: (i) optimizing the position of the second stable configuration and (ii) optimizing the energy required for switching between the stable configurations. The first objective serves as a proxy for the local distribution of geometry and mass of the second equilibrium position, which is important for controlling the compactness and dynamic behavior of the stable buckled state. The second objective is important for understanding the overall energetics of the local beam mechanics and serves as a proxy for the stability of the second equilibrium state. In the proposed design space, the thickness of the beam along the length direction is prescribed using a Fourier series, with the coefficients serving as the design variables. The Fourier series representation guarantees a continuous shape for the beam that is directly manufacturable, without incorporating a post-process thresholding step often required in SIMP-based optimization approaches [10,11]. Fourier series have also been used to parameterize the geometry of periodic voids in soft, elastomeric lattices, with the specific goal of characterizing the role of porosity in tuning the force-displacement profile and Poisson’s ratio of architected networks [12,13]. For a given beam design, the bistability behavior is identified by the force-displacement response, and nonlinear FEA is carried out with the robust displacement-controlled Newton Raphson method [10,14–16]. Although this particular design parameter formulation is inherently smooth and continuous, the objectives will be very sensitive to the design variables and are expected to be highly non-convex. Under these conditions, gradient-based optimization is often attracted to only locally optimal designs.

To navigate this complex design space, Bayesian optimization is implemented in this study to find an optimal beam shape with respect to its bistability behavior. Bayesian optimization is a global optimization method that generally works well for highly non-convex problems and especially for those in which the response of interest contains a significant amount of noise [17–21]. Bayesian methods have previously been leveraged to classify the attenuation behavior in elastomers loaded with negative stiffness inclusions [22–24]. The primary advantage of Bayesian optimization is the potential for improved convergence over evolutionary stochastic approaches with fewer function evaluations, which is important for this study due to the required nonlinear FEA. At each iteration of the Bayesian optimization, a Gaussian process (GP) surrogate model is estimated, based on training data from the FEA model. The surrogate model provides insight into the global behavior of the design space, which can also build the designer’s intuition. Based on the distribution of the objective function, the Bayesian optimizer then specifies the design update. Challenges of implementing Bayesian optimization are the selection of the best covariance function to predict the objective function using the GP and the determination of a good balance between exploration and exploitation. Once those are determined, the hyper-parameters (i.e., kernel length scale and kernel scaling parameter) can be optimized by maximizing the log-likelihood function. A detailed analysis of the kernel and acquisition hyper-parameters for the bistable beam design space is included in this study. The use of appropriate weights and exponents in the weighted compromise programming objective function also is critical to finding good trade-off designs during a scalarized multi-objective optimization [25,26]; thus, the process used to select these values is considered in this study. The low-dimensional design optimization results are compared with a grid search in this research study, which itself provides a useful reference database for the fabrication of multistable elastomeric lattice structures with various target equilibrium states. The database can be further utilized for future research, e.g., the classification of bistable and monostable designs using a machine-learning-based approach.

A key challenge of the buckling beam design space is the occurrence of monostable solutions, which are undefined for bistable objective functions and result in discontinuities. Design space discontinuities due to bifurcations in the mechanical behavior of the system are pervasive phenomena and a computational challenge. Techniques to handle transitions in mono- and bistable performance are varied, as are their respective cost savings and ease of implementation. Previous work from the composite laminate community provides instructive context for this class of problem, as mono- to bistable transitions occur in these design spaces [27,28]. For example, Hufenbach and Gude implemented a genetic algorithm (GA) to design the number of layers and specific properties of a laminate stack to maximize the curvature in the reference configuration [27]. Had bistability been a design objective for this study, a simple approach to incorporate this criterion into the GA framework could have been to remove all mono-stable candidates during design evolution. However, this would lead to inefficiencies as ineligible designs would continue to be evaluated and consume computational resources. Follow-on studies addressed this issue for the composite laminate design space by constraining the design space to remain within the bistable domain [28], which can be formulated as an analytical expression for the laminate problem. For many problems, the domains of mono- and bistable design solutions are discontinuous, disjoint, and not amenable to the analytical description. Discovering the boundaries of these domains using machine learning classification techniques [29,30], or smart sampling strategies to identify feasible domains [31] are promising strategies to address his problem. For example, a Bayesian Network Classifier (BNC) based on a kernel density estimate was used to find the geometry group of bistable inclusions that enables high performance of a loaded elastomeric material [23]. A similar approach could be applied to classify bistable and monostable domains. Alternatively, the objective function can be formulated to target bistable designs by explicitly shaping the force-displacement response of the structure at key locations [10,16,32]. The method is effective at driving toward solutions in the bistable design domains, but may artificially constrain the design space by over-specifying the character of the force-displacement curve. In this study, we combine the ideas of a machine learning-based surrogate model with force-displacement curve manipulation by introducing a penalty term for monostable design candidates in the GP model representation. The penalty term enables monostable design candidates to be incorporated in the surrogate representation and indirectly biases the search space toward desired bistable designs.

The bistable, elastomeric beam design optimization study presented in this manuscript is structured in the following way: Sec. 2 outlines the parameterization of the design space for the beam shape and defines the constitutive model for the material and the underlying nonlinear behaviors of interest. Section 3 describes our implementation of the Bayesian optimization method and corresponding hyperparameter tuning. Section 4 defines the two objective functions of the study, specifically to optimize the position of the second stable configuration and the energy required for switching between the stable configurations. A penalty factor is also developed in Sec. 4 to incorporate monostable designs in the GP surrogate, as the position of the second stable point is undefined in monostable cases. Section 5 presents the optimization results for both of the objectives, plus a multi-objective design study where the objectives are combined using a weighted compromise objective formulation. The values of the weights are identified through calibration using a database of a smaller design domain, generated through an exhaustive search. Lastly, the qualitative trends in the beam bistability behavior are experimentally validated in Sec. 6.

As described in the Introduction, this study utilizes the inclined, elastomeric beam from Shan et al. [1] as the fundamental unit for further optimization. The bistability of the elastomeric beam within the structure exists under the specific boundary condition shown in Fig. 1(a), where q_ext is the applied external force and where u is the corresponding displacement. Here, the resulting vertical displacement and constraints against horizontal motion enable a bistable configuration. Under these specific boundary conditions, the mechanical behavior of a straight, prismatic beam depends on inclined angle (θ_o), slenderness ratio (t_o/L_o), type of elastomeric material and cross-sectional design.

The bistable behavior for a given beam design is characterized by a nonlinear force-displacement response, which can be obtained with FEA using a Newton–Raphson method or arc-length continuation scheme [10,15,16]. A two-dimensional (2D), plane-strain FEA with free bilinear quadrilateral dominated elements is performed by the commercial FEA tool ABAQUS in this study [33]. Here, quasi-static FEA is carried out while neglecting any dynamic/inertial effects. In the Newton–Raphson method, the force–displacement relationship is obtained by iteratively applying a constant-step incremental force to the structure and solving the change in the displacement. The process is repeated until the full load is applied. Nonlinear analysis is sensitive to the size of the load increment. To strike a balance between computational cost and accuracy, the total applied displacement, u_b, was divided into 200 steps of equal size (Δu = 43 μm, u_b = 8.6 mm). See Supplementary Material—Section I available on the ASME Digital Collection for more details on the implementation.

Figure 1(c) shows the force-displacement response of a bistable beam, which includes a snap-through bifurcation and two stable configurations (u₀ = 0 and $u*$⁠). The energy barriers between the two stable equilibria are annotated as E_in and E_out, and are defined as the areas under the force–displacement curve as shown in Fig. 1(c). The relative sizes of these areas serve as a proxy for the stability of each stable configuration. For this reason, we will utilize the energy ratio, $s*=Eout/Ein$⁠, as an objective function to tune the relative energies needed to switch between stable states. In the case where the beam designs are monostable, the displacement at the second stable state does not exist (⁠$u*$⁠, undefined), and correspondingly, E_out is zero. Monostable designs inherently create discontinuities in the $u*$ design space, because they are undefined. A penalty displacement approach is proposed to address this problem, and a detailed explanation is presented in Sec. 4.2.

In general, a smaller value of ν leads to a less-smooth covariance function. In this research, the different covariance functions in Eqs. (5) and (7) are studied to best predict response functions.

The Bayesian optimizer makes decisions using the distribution of the GP predicted response while balancing exploration and exploitation, which can be scalarized using the acquisition function [19,20]. Exploration corresponds to searching the design space where large uncertainties in the response exist, while exploitation consists of seeking the design space where the average response is minimized. If optimization is biased toward exploration, it can become inefficient (requiring more iterations), while if it is biased toward exploitation, finding local optima instead of global ones can occur. In short, while the GP predicts performance throughout the entire design space with quantified uncertainty based on all designs evaluated so far, the acquisition function directs the optimization process based on the GP. The acquisition function can be designed to be attracted to both regions where the GP predicts good performance and where the GP uncertainty is high.

Both accuracy and efficiency of Bayesian optimization highly depend on optimization parameter tuning. In this study, the smoothness of the Matérn kernel (ν = 1/2, 3/2, 5/2), the trade-off parameter between exploration and exploitation of the acquisition function (β = 1, 3, 5), and the penalty displacement applied to monostable designs are systematically investigated. The process for selecting the ideal ν and β values, which were determined to be ν = 1/2 and β = 3, is discussed in Supplementary Material—Section II available on the ASME Digital Collection. For multi-objective case studies concerned with balancing both the displacement and energy ratio, it is beneficial to have a complete understanding of the trade-offs involved. Therefore, a method of finding and verifying the Pareto frontier is proposed and implemented.

One of the objective functions of this study is the displacement at the second stable configuration, J_u*; however, for monostable training data, this objective is undefined. Nevertheless, it is still desirable to include the monostable training data in the BO process to increase the efficiency of the optimization and to improve the representation of the response functions. To achieve this, a large, penalty displacement is assigned to monostable designs to define a fictitious second stable configuration for incorporation in the GP surrogate. The penalty displacement approach simultaneously increases the number of data points incorporated in the GP and pushes the GP away from monostable designs, due to the high displacement values assigned. To identify appropriate values of the penalty, we performed a parameter study using two coefficients of the first Fourier series term as design variables, with t/10 ≤ r₁ ≤ t/4 and 0 ≤ p₁ ≤ 8. The bistability of the beam only exists within a certain range of slenderness ratios (t_o/L_o) and inclination angles (θ_o), where t_o/L_o and θ_o are 0.12 and 40 deg, respectively, and t_o is 1.0 mm. 100 training data points are used for GP prediction with the squared exponential covariance function to determine the penalty displacement. Different penalty displacements (5.0, 6.0, and 7.0 mm) are applied to the monostable training data, and the penalized GP with the best prediction is selected.

Figure 2(a,i) shows a contour plot of GP predicted surface for J_u* with only bistable training data (76 of 100), in terms of the magnitude change (r₁) and periodicity (p₁) of the beam design. The GP surrogate predicts a weakly bistable beam for the design with the smallest displacement, as shown in the force–displacement response of Fig. 2(a,ii). Interestingly, the GP surface indicates a basin of small displacements between the periodicities of 2 and 6, which intensifies with increasing magnitude. Figures 2(b,i)–2(d,i) show the results of the GP predictions, with different values of penalty displacements (PD = 5.0 mm, 6.0 mm, and 7.0 mm) applied to monostable training data. The force–displacement responses of the optimal designs obtained from GP predictions from Figs. 2(b,i)–2(d,i) are all shown in Figs. 2(b,ii)–2(d,ii). The obtained optimal design, with a penalty displacement of 5.0 mm, is monostable, indicating the error in the GP estimate. With both penalty displacements of 6.0 mm and 7.0 mm, the obtained optimal design is bistable (when the penalty displacement is 6.0 mm, the energy ratio at the second stable is nearly zero $(s*≅0)$⁠). The error of prediction from GP from the FE model with penalty displacement equal to 6.0 mm is found to be 6.1%.

Figures 2(a,i)–2(d,i) indicate that $u*$ is highly nonlinear, as bistable and monostable regions are not clearly divided. Moreover, while little effect of $u*$ to magnitude change is observed on GP prediction without monostable training data (Fig. 2(a,i)), GP-predicted surfaces significantly vary in terms of magnitude change, when penalty displacements are applied. The penalty approach enables the monostable data to influence the GP surrogate and improve the prediction of bistable designs. However, the penalty creates a sharp transition in the response surface, which may be difficult for the GP surrogate to capture. According to Fig. 2(c,i), the optimal region is effectively localized in the GP-predicted surface with the penalty function of 6.0 mm, which will be used in the design case studies.

Based on the optimization parameter tuning carried out in this study, the design optimization is performed with the following parameters: β = 3, ν = 1/2, and penalty displacement of 6.0 mm. At every iteration, the corresponding hyper-parameters, which are the length scale and the noise level, are all optimized by maximizing the log-likelihood function. The gradient-based method (gradient ascent) was used, and the bounds of the search space were set to be from 10⁻⁶ to 10⁶ for both parameters. It was repeated 60 times at random initial hyper-parameter values to enhance the accuracy. The baseline design and the design space for the one-term Fourier series design are identical to the ones used in Sec. 4.2. While the two coefficients of the first term are used as design variables for the one-term Fourier series design, the two coefficients for the first term and the two coefficients for the second term are used as design variables for the two-term Fourier series design. Different ranges of periodicities for the first and second terms are assigned to understand their collective effects on design objectives. The design space for the first term is identical to the one for the one-term Fourier series design, and the design space for the second term is t/20 ≤ r₂ ≤ t/10 and 10 ≤ p₂ ≤ 18 for magnitude and periodicity, respectively. Material failure is not considered in this study as the applied strains across all designs are well below the elongation strain at break of the elastomer used in the experiment (Mold Max NV29 –316%). However, for material systems where failure is a concern, a stress- or strain-based inequality constraint could be implemented into the Bayesian optimization framework to avoid failure-prone designs [35].

The first design problem is to minimize the displacement of the stable configuration (⁠$Ju*=u*$⁠). The search history for this objective is shown in Figs. 3(a)–3(c), which shows that it converges efficiently to the optimal design within 25 iterations for the one-term Fourier series design. On the other hand, it converges after 113 iterations for the two-term Fourier series design, which is much slower than the one-term design problem. While the GP-predicted objective function values are shown in Fig. 3(a), the actual objective values obtained by FEA for one-term Fourier series design are also shown in Fig. S3(a)—Section III available in the Supplemental Materials on the ASME Digital Collection. It should be noted that the actual objective function value at the optimal design is not always calculated by FEA during the Bayesian optimization; instead, the design suggested by the acquisition function is evaluated, as explained in Sec. 3.2. For the one-term Fourier series design, the obtained optimal change of magnitude is 0.25 mm (maximum allowable), and the optimal periodicity is 5.1. For the two-term design, the optimal magnitude and periodicity of the first term are 0.232 mm and 5.06, respectively, and those of the second term are 0.063 mm and 17.4, respectively. The magnitude changes were smaller than the max allowed amounts since there is a joint effect from Fourier series terms for lower and higher periodicities of the beam design to minimize u*. This is also confirmed by the fact that the optimal periodicities obtained from two Fourier series terms do not exist at the bounds of the design space.

The minimized displacement was 4.58 mm for the one-term Fourier series design, and 3.88 mm for the two-term design. Compared to that of the baseline design, the displacement of the one-term and the two-term optimal designs was reduced 42% (4.58 mm versus 7.90 mm) and 51% (3.88 mm versus 7.90 mm), respectively. The force–displacement responses of the optimal designs are shown in Fig. 3(d), and the geometries of the optimal design at two stable configurations are shown in Fig. 3(e). The summary of optimal designs and the corresponding objective function values are shown in Fig. 3(f).

The objective of the second optimization problem is to maximize the ratio of output to input energy between stable states $(Js*=s*)$⁠. As shown in Figs. 4(a)–4(c), the design search converges efficiently to the optimal configuration within 27 iterations for the one-term Fourier series design, while it converges to the optimal design within 96 iterations for the two-term design. For the one-term Fourier series design, the obtained optimal change of magnitude was 0.25 mm (maximum allowable), the optimal periodicity was 6.4, and the maximized $s*$ was 0.52. For the two-term design, the optimal magnitude and periodicity of the first term are 0.25 mm and 6.53, respectively, and those of the second term are 0.10 mm and 15.8, respectively, and the maximized $s*$ was 0.75. It can also be seen from this case study that both optimal periodicities are bounded within the design space, which indicates the collective role of lower and higher periodicities of the beam design.

Compared to the baseline beam design $(s*=0.20)$⁠, the optimized energy ratio increased 2.6× and 3.75× for the one- and two-term Fourier series designs, respectively. It was also observed that the second term for the higher periodicity of the beam design more dramatically improved the design objective, compared to the $u*$ minimization problem. The force–displacement responses of the optimal designs are shown in Fig. 4(d), and the geometries of the optimal design at the two stable configurations are shown in Fig. 4(e). The summary of optimal designs and the corresponding objective function values are shown in Fig. 4(f).

In order to understand the trade-offs required between $u*$ and $s*$ in multi-objective optimization, and especially in order to understand the effects of the two weighting scalars (w₁, w₂) in Eq. (13), the Pareto frontier was sought. The accuracy of this Pareto frontier was then validated by comparison to a grid search of the entire design space, which was taken as ground truth. Due to computational expense, the grid search is limited to the one Fourier series-term problem, since it involves a smaller design space. The resulting database of FE models and the analysis results during the grid search is also provided for use in future research (see Supplementary Material—Section IV available on the ASME Digital Collection).

Regardless of whether the Pareto frontier happens to be convex or not, all points on the Pareto frontier can be captured using Eq. (13) as the objective function when w₁ and w₂ are varied from 0 to 1 such that their sum is always 1, if a sufficiently large value of parameter η is chosen [25,26]. However, higher values of η increase the non-convexity in the objective function and therefore the response surface, and Bayesian optimization may settle on local optima under highly non-convex conditions; thus, care must be taken to choose an appropriate η value. To achieve this balance, the value of η was varied from 5 to 50. Some example results are shown in Fig. 5(d). With smaller η values, designs with midrange energy ratios (around 0.3–0.4) are not captured, indicating features in the frontier that those η values are too small to capture. With larger η values, the frontier seems to settle on many local optima, indicating that these large exponents induce excessive non-convexity that prevents proper convergence for this design space. While there are other Bayesian optimization formulations which are capable of discovering the Pareto frontier without the need for a weighted compromise programming-based objective scalarization with a weight sweep, such as those relying on hypervolume-based methods [36], we selected this formulation for consistency with the single-objective optimization procedure.

The η value of 12 used in Eq. (13) captures the most complete frontier and is therefore selected for all scalarized multi-objective case studies in Sec. 5. It is noted that this procedure involved some discretion in the selection of η value range and visual evaluation of the results, leaving open the possibility that some accuracy was sacrificed for efficiency. The results were therefore compared to a grid search which covered the entire design space, where r₁ varied from 0.10 mm to 0.25 mm (step size of 0.01 mm) and p₁ varied from 0 to 8 (step size of 0.10), which totaled 1296 FE model evaluations. Compared to a Bayesian optimization run which could require as few as 25 FE model evaluations (as discussed in Sec. 5), the grid search was quite computationally expensive but was undertaken to confirm that the Pareto frontier was accurate and in particular to provide additional confidence in the approach used to obtain it. Because the weighted compromise programming scalarization induces greater non-convexity in the response surface, it is possible that a new exploration-exploitation tradeoff could be required. However, as shown in Fig. 5(c) in which the grid search results are shown in gray, the frontier obtained through Bayesian optimization (η = 12) in the objective function and the originally selected β and ν values very accurately represents the trade-off between displacement and energy ratio. This provides the basis for all the multi-objective case studies, which will be applied to the higher-dimensional problem.

In addition to validating the Pareto frontier results, the grid search provides deeper insight into several other important items, including the effects of the penalty displacement applied to monostable beams, optimization parameter choices, and the organization of designs in the objective function space. The results of FE model evaluations carried out for all combinations of 16 by 81 discretized r₁ and p₁ values, respectively, are shown in Figs. 5(a) and 5(b) in the design space and in Figs. 5(c) and 5(d) in objective function space.

It is observed previously that incidence of monostability increases when the displacement is minimized or the energy ratio is close to zero. Figures 5(a) and 5(b) confirm this as well, as shown in the displacement-minimized Pareto point (marked as a black solid circle), for which the energy ratio is also minimized. Based on this, we can see that GP-predicted surfaces in Figs. 5(a) and 5(b) agree well with training data from the fact that both displacement and energy ratio are most minimized adjacent to the clearly defined monostable region. This characteristic is also observed in Fig. 2(ci) in the GP, which is trained on only 100 points with a penalty displacement of 6 mm. Moreover, the shape of the displacement GP in Fig. 5(a) is very similar to the one in Fig. 2(c,i), which further indicates that the optimization parameters selected were appropriate.

In Figs. 5(a) and 5(b), the grid search results shown in gray appeared to be arranged in a pattern of rays or lines of points, with a “shelf” structure near the bottom right of the plot in Fig. 5(c). Further analysis revealed that periodicity plays a particularly important role. For example, as shown by the blue- and purple-circled points, all points of periodicity 1 fell within one line, directly adjacent to another line formed by all the points of periodicity 1.1. The top of the shelf structure on the bottom right was composed of points of p₁ = 3.8, circled in brown. As indicated in Fig. 5(a), p₁ = 3.8 lies just above a bistable region, suggesting a relationship between the empty space above p₁ = 3.8 in Fig. 5(d) and the bistable region just below p₁ = 3.8 in Figs. 5(a) and 5(b). The side of the shelf structure is composed of points ranging from $p1=[3.8,5.1]$⁠, all of r₁ = 0.25, as circled in green. In Figs. 5(a) and 5(b), this range of points is directly bordered above and below by the two bistable regions, and on its side by the maximum magnitude (0.25). This suggests that improved performance (i.e., points to the left of this shelf of green-circled points) could be obtained by allowing magnitudes greater than 0.25. However, this maximum was enforced throughout the study because greater magnitudes would have led to designs with portions that were too thin to manufacture without the risk of breaking.

It was noted that all optima from one-term Fourier series beam case studies fall on the Pareto frontier, as expected. On the other hand, a limitation of penalty-term-based Bayesian optimization is also revealed. In Fig. 5(a), the Pareto point with the lowest displacement (from Fig. 5(c)) is shown at $[r1,p1]=[0.207,5.06]$⁠, while the grid search revealed that the design of $[r1,p1]=[0.22,5.1]$ had, in fact, slightly lower displacement than the lowest-displacement Pareto point, as indicated by the red stars in Figs. 5(a)–5(c). This limitation is mainly due to the fact that, during optimization, the penalty displacement likely increased overall the GP-predicted displacement values in the vicinity of the monostable designs, driving the optimizer to a more conservative design (namely, $[r1,p1]=[0.25,5.1]$⁠, as indicated by the + symbol). As a result, while the grid search locates the design directly adjacent to the monostable region, the Bayesian optimizer locates the design marginally farther away from the monostable region. While eliminating all monostable beams from the running training set during optimization would have eliminated this effect, it would also have likely increased the number of iterations before convergence as discussed in Sec. 4.2. Furthermore, as revealed by Sec. 5.1, the design, which is close to monostable, will likely be subject to weak stability.

As shown in Figs. 4 and 5, the optimal designs are significantly different, depending on the objective (⁠$u*$ or $s*$⁠). When the objective is to minimize displacement at the second stable state, the energy ratio to retain equilibrium becomes very weak, as indicated by Fig. 4(d). This can be resolved by carrying out a multi-objective optimization. Based on Sec. 5.3, the compound objective with m value of 12, which covers most of the Pareto frontier, can be used. The specific weight combinations (w₁, w₂) used in this example for one- and two-term Fourier series designs are (0.734, 0.266) and (0.900, 0.100), respectively. The goal of the multi-objective case study is to best maximize the energy ratios, while also minimizing the equilibrium states.

The search history is shown in Figs. 6(a)–6(c). The numbers of iterations taken for this problem for one- and two-term Fourier series designs are 33 and 108, respectively, which are in similar ranges to the ones taken in previous design optimization problems. For the one-term Fourier series design, the suggested magnitude change for the optimal compromised design is 0.25 mm, which is identical to the one suggested for the optimal design for $u*$⁠. The suggested periodicity is 4.74, which is only 6.3% smaller than the one suggested for the optimal design for $u*$⁠. This indicates that the sensitivity of the energy ratio with respect to the design variable is very high and also implies the non-convexity of the response surface, considering that a higher periodicity is previously suggested to maximize $s*$⁠. For the two-term Fourier series design, the optimal magnitude and periodicity of the first term are 0.25 mm and 5.06, respectively, and those of the second term are 0.05 mm and 10.2, respectively. While the periodicity of the first term remains identical to the one of the $u*$ minimized design, the periodicity of the second term changed substantially to achieve the additional energy ratio objective.

Force–displacement responses of the optimal designs are shown in Fig. 6(d), where it can be noticed that stabilities of the designs were improved compared to those in Fig. 4(d). For the one-term Fourier series design, the energy ratio indeed significantly increased from 4.61% to 14.9%, while the displacement at the second stable state increased from 4.58 mm to 5.01 mm. For the two-term design, the minimized displacement and the energy ratio were 4.52 mm and 7.01%, respectively. Compared to the one-term Fourier series design, the displacement at the second stable state of the two-term design increased much less, although the energy ratio improvement slightly decreased. This shows the contribution of the second term to the design objective and also the necessity of carefully selecting weight ratios. The geometries of the optimal designs at two stable configurations are shown in Fig. 6(e), where it can be seen that vertical displacements at buckled states are slightly higher than those of u* minimized designs. The summary of optimal designs and the corresponding objective function values are shown in Fig. 6(f). Their values depending on the target design objectives (J_u*, J_s*, and J_m*) are additionally organized in Supplementary Material—Section V available on the ASME Digital Collection for the comparison.

The experimentally tested results on the computationally predicted optimal designs in Sec. 5 are presented in this section. In the experiment, Mold Max 29NV was used as the elastomeric material. In order to ease fabrication, handling, and testing, all samples were built to 5 cm in length, which is about six times the length of the computational model. We thus intend to demonstrate the computationally predicted optimal designs for J_u*, J_s*, and J_m through relative performances of the experimental designs instead of through absolute value comparison. A mold casting method is used to fabricate the proposed beams since the purpose of this experiment is to accurately test the prototypes to demonstrate the computational study. The molds were fabricated using a fused deposition modeling 3D printer (Flashforge Creator Pro, ABS filament). To facilitate mounting of the specimen during mechanical testing, slots for glass substrate were added to the mold designs so the elastomer could be casted and bonded to the rigid substrate. The elastomer was cured at room temperature for over 6 h. See Supplementary Material—Section VI for more details on the experimental setup.

Three-dimensional (3D) printed custom adapters (Stratasys, Objet260, VeroWhite ink) were attached to a Tinius Olsen H10K-S load frame, and the beam samples were attached to the adapters via their glass substrates (Fig. S5(c) available in the Supplemental Materials. The compression test results are shown in Figs. 7(a,i)–7(c,i), where the design objectives, which are displacement at second stable state and the energy ratio, were substantially improved compared to the ones of the baseline design, as predicted by the computational study. Identical compression tests were repeated three times, and the low deviation shows the robustness of the results (Supplementary Material—Section VII. According to Fig. 7(a,i), the $u*$ minimized design obtained with one Fourier series term became bistable, while the one obtained with two Fourier series terms turned out to be monostable or very weakly bistable. This issue has been mentioned earlier in Sec. 5 and has been resolved through the scalarized multi-objective optimization. In Fig. 7(c,i), it is shown that both multi-objective optimal designs obtained with one and two Fourier series terms became sufficiently bistable by significantly improving stabilities at the second stable state. In Fig. 7(b,i), while the energy ratio maximized designs certainly improved over the baseline design in terms of the design objective, their force–displacement responses are noticeably different from the computational prediction. Unlike the computational response, a sharp peak is observed in the experimental response, which is due to the fact that the design is subject to a bifurcation phenomenon and followed an alternative path of buckling.

Geometries of the computationally predicted optimal designs and experimentally tested specimens before and after they are buckled which are shown in Figs. 7(a,ii)–7(c,ii), where it can be seen that the shapes of computationally predicted beams and experimentally tested elastomeric beams are in very good agreement. Even for the energy ratio maximized designs that are subject to buckling, the shapes in the buckled states are almost identical to those of the computationally predicted ones. On the other hand, the slight difference was also observed that displacements at the second stable states of the experimentally tested specimen were generally smaller than those for the computationally predicted designs. Overall, the good agreement between computational prediction and experimental demonstration shows the strong potential of utilizing the proposed periodic beam designs for wider applications in architected elastomers.

A design optimization method has been developed for a bistable elastomeric beam unit cell for wider applications in architected elastomers. The geometry of the beam was parametrized by Fourier series coefficients corresponding to the magnitude and periodicity, which served as design variables. Since the response of bistability behavior in terms of design variables is highly non-convex, and computationally heavy evaluation of nonlinear FEA is required during the optimization, Bayesian optimization was implemented to improve the convergence with respect to the number of high fidelity function evaluations.

Important optimization parameters of the Bayesian optimization are systematically studied for both accuracy and efficiency. When the optimizer is biased toward exploitation, it converges to the optimal design most quickly. Convergence also improves with Matérn covariance functions using a ν value of 0.5, which indicates non-smoothness of the objective function in the design space. Since the displacement objective is undefined for monostable designs, a penalty displacement was developed to effectively incorporate monostable results in the GP surrogate. Design space bifurcations, such as the transitions in mono- to bistable buckling in this study, are a pervasive challenge for optimization processes. The penalty method created a suitable proxy for bistability to mitigate the effect of the bifurcation; however, this came at the expense of reduced resolution of the phase boundary between stabilities, as discussed in Sec. 5.3. An alternative approach to the penalty function would be to identify the mono- and bistable domains and then parameterize the design space to search within the desired stability class. Bayesian classifiers have been applied to address this challenge of identifying and grouping stability domains [24] and show promise for a “partition and optimize” design strategy.

The proposed Bayesian design optimization framework efficiently finds designs that improve the objectives based on the force–displacement response, even under the design space with a single Fourier series term. It is also demonstrated in this study that the objectives can be further improved under a larger design space with an additional Fourier series term, with the trade-off for computational cost. Both the accuracy and the robustness of the computational results are demonstrated through the experiment in this study. The deviation between some of the numerically predicted and experimental results suggests that future research is required to find designs which also consider uncertainty in manufacturing as part of the optimization. The library of beam designs appended with this manuscript provides initial insights to explore these and other design questions, as well as a training dataset to leverage with other machine learning surrogate modeling and optimization approaches. In particular, the beam library may serve as a key design reference for the construction of elastomeric lattices. Exploring this question of networking and the stability of beam networks is a promising future direction of this research.

The authors acknowledge the support of the Air Force Research Laboratory. The experimental part of the work was funded by AFOSR 17RXCOR435. The manuscript is an extended version of the authors’ 2019 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference proceedings titled, “Bayesian Optimization of Equilibrium States in Elastomeric Beams.” Authors NH, SA, and KV wish to make the following acknowledgment with respect to their contributions: “This material is based on research sponsored by the Ohio Department of Higher Education and the Southwestern Council for Higher Education under Ohio House Bill 49 of the 132nd General Assembly. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Southwestern Council for Higher Education, the Ohio Department of Higher Education or the U.S. Government.”

There are no conflicts of interest.

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The data and information that support the findings of this article are freely available online.² The authors attest that all data for this study are included in the paper.

$h$ =

joint distribution of training data

$z$ =

Training inputs

k =

covariance function

N =

normal distribution

$h*$ =

joint distribution of output responses of interest

$z*$ =

test inputs

$s*$ =

energy ratio

$s^*$ =

energy ratio of baseline beam

$u*$ =

displacement at second stable state

$u^*$ =

displacement at second stable state of baseline beam

$kMat′ern$ =

Matérn of covariance function

k_SE =

squared exponential covariance function

p_i =

periodicity of ith Fourier series term

r_i =

magnitude of ith Fourier series term

t_o =

thickness of baseline beam

u_b =

applied vertical displacement

u_c =

displacement at metastable state with zero reaction force

u_o =

displacement at first stable state

E_in =

input energy

E_out =

output energy

J_m =

weighted compromise programming objective function

J_s* =

energy ratio objective function

J_u* =

displacement at second stable state objective function

L_o =

length of baseline beam

a_LCB =

lower confidence bound acquisition function

β =

trade-off parameter that controls between exploration and exploitation

η =

parameter that determines non-convexity of J_m

θ_o =

angle of inclination of baseline beam

μ =

mean of predicted response surface

ν =

parameter that determines the smoothness of Matérn covariance function

σ =

standard deviation of predicted response surface