Design of Phononic Bandgap Metamaterials Based on Gaussian Mixture Beta Variational Autoencoder and Iterative Model Updating

Controlling the energy flow attracts great attention in the design of mechanical metamaterials [1–3]. Periodic cellular phononic bandgap metamaterials have demonstrated the capability to influence the propagation of mechanical waves within a certain frequency range [4]. Vibrations within the bandgap range are attenuated by the mechanism of wave interferences within the periodic structures, so the waves are blocked by the metamaterial. Phononic bandgap metamaterials show their advantage in a variety of applications, such as vibration isolation [5], noise reduction [6,7], and stress wave protection [8–10]. It has been demonstrated that bandgaps spanning different frequency ranges can be achieved by tailoring the topological features of the cellular structures [11–16]. A rational methodology for discovering new metamaterials with a wide bandgap is highly desired. Topology optimization methods have been proposed to pursue specific bandgap properties (e.g., maximum total bandgap width, maximum bandgap width in the desired frequency range, or achieving bandgap in a specific frequency range) [17–22]. Topology optimization shows its strength in exploitative design optimization, but it is not effective for the generative design of a highly diversified design set.

Deep learning, as a subject of artificial intelligence (AI), has been successfully applied to a variety of research fields [23–28]. Here, we focus on the deep learning-driven design methods for computational material design [29–37]. Deep generative models, such as variational autoencoders (VAEs), generative adversarial networks (GANs), and their variations, have been used to compress the high-dimensional input data into a low-dimensional latent space for design representation and design generation. A GAN model converges when the discriminator and the generator reach a Nash equilibrium [38]. The generator model learns to generate new images and the discriminator learns to decide whether the newly generated images are true or false. Unlike GAN that does not work with any explicit density estimation, VAE provides an explicit likelihood for modeling training. The quality of the VAE model can be evaluated easily based on its log-likelihood loss function, while GAN is notoriously difficult to train and requires a lot of tuning [39]. Wang et al. proposed a data-driven metamaterial design framework based on VAE to learn encoded structure and property information for the design of metamaterial structural systems that achieve desired distortion patterns [40]. Tan et al. applied a deep convolutional generative adversarial network and convolutional neural network to design microstructural materials. The optimal designs achieve desired compliance tensor matrix and satisfy the geometrical constraints [41]. GAN has also been applied to design complex architectured materials for stiffness and Poisson ratio [42]. However, only elastic or elastic-related properties are considered in the works mentioned above. The deep generative modeling methods have been applied to design phononic/photonic metamaterials as well. Li et al. proposed an autoencoder-based method to design phononic crystals [43]. Ma et al. proposed a VAE-based approach to design photonic structures [44]. The VAE-based inverse design framework proposed by Liu et al. is generally applicable to the design for various material properties [45]. GAN-based design frameworks have been developed to design photonic structures with user-defined objectives [46–48]. However, these works only consider metamaterial designs with simple topological features (e.g., single-hole structure) and a single design objective. For multiobjective design, Kim et al. proposed to integrate VAE and Gaussian process (GP) regression to find optimal microstructure designs that achieve both ultimate tensile strength and toughness [49].

The objective of this study is to find high-performance designs with a good diversity, which would incur a significant computational cost with the traditional design methods. In this work, we propose an effective design framework that enables both unsupervised generative design and exploitative optimization of cellular metamaterial structures with multiple design objectives. The unsupervised generative design can be done by exploring the latent feature space by space-filling sampling techniques (e.g., Latin hypercube sampling (LHS)) [50]. The random space-filling process may not be the most effective way with limited computational resources. To resolve this issue, we propose a new approach based on the modes of the latent features learned by the Gaussian mixture beta variational autoencoder (GM-βVAE). GM-βVAE enables unsupervised clustering through deep generative models [51]. The exploitative optimization is based on the supervised learning model of the relationship between the latent features and the design performances, such as GP regression [52] and neural network [53]. The optimal latent feature values are obtained by optimization and then mapped back to the original structure design space to reconstruct the corresponding structure (decoding). In this research, we establish and compare a single-loop supervised learning-based design approach and a new iterative model updating-based design framework. The proposed methods are expected to be generally applicable to design problems where the gradient information is not available or too expensive to obtain (e.g., structure crashworthiness, phononic bandgap).

The remainders of this paper are organized as follows. In Sec. 2, fundamental physics and the simulation model of phononic bandgap properties are introduced. A convergence study is presented to determine the optimal mesh size. In Sec. 3, we present the deep generative model-based design methodology: a single-loop approach and an iterative GM-βVAE model updating-based design framework. The deep generative model is selected based on a comparative study of the learning capabilities of the VAE, βVAE, and GM-βVAE models. In Sec. 4, we first present the GM-βVAE-based unsupervised generative design method and compare it with the traditional space filling-based unsupervised generative design method. Then we showcase the application of the single-loop approach on property-driven design optimization. Finally, the iterative model updating-based framework is demonstrated on a multiobjective design problem. The optimal designs are compared with the topologically optimized designs reported in literatures. Section 5 is the conclusion.

Generally, fine meshes of the unit cell are preferred in terms of the accuracy of the solutions to the eigenproblem. The expense of the simulation, however, increases as the meshes are refined. Therefore, a convergence test was conducted in pursuit of the mesh refinement that balances accuracy and efficiency. We randomly picked two different cellular metamaterials in our database to perform the convergence test. In the test, refinements of the meshes are quantified as the size of 4-node elements along the edges of the structures. Refinements are set as 1/50, 1/125, 1/250, and 1/375, from coarse to refined. Corresponding dispersion characteristics are estimated. We compared, in Fig. 2, for two different unit cells. Only the first five curves are shown for clarity. Assuming the 1/375 as the finest possible resolution that can be achieved with available computational resources, we evaluated the accuracy of the rest resolutions. Symbols of line-circle, line-star, line-cross, and line-triangle represent reduced natural frequencies with respect to different reduced wave vector K, for 1/50, 1/125, 1/250, and 1/375, respectively. Lines are guides to the eyes. As clearly shown, the 1/50 meshes have obvious deviations from the references 1/375 results. The 1/250 meshes, on the other hand, show good consistency to the references but the computational expense is out of our expectation in terms of the generation of a database with thousands of structures. Therefore, the choice of meshes 1/125 is taken as the compromise. Sampled with 60 discrete sampling locations on wave vector K, simulations for a single cellular metamaterial take about up to 30 min, up to 90 min, up to 180 min, and up to 300 min for 1/50, 1/125, 1/250, and 1/375, respectively. There is no surprise that simulations for sparse structure (Fig. 2(a)) are faster than those for dense structure (Fig. 2(b)).

Based on the convergence study, we choose to use 250 × 250 pixel image (correspond to the refinements of 1/125) to represent our input data for the sake of a balance between calculation time of bandgap simulation and the accuracy. Note that the size of input data is also used in the simulation of elastic properties and the generative model for the sake of consistency. We automated the simulation-postprocessing process and implemented it on our cellular metamaterial dataset. The total bandgap of a cellular metamaterial is selected as a property of interest. The statistics of all samples in the database are shown in Fig. 3(a) and Table 1.

Simulations of the cellular metamaterials’ stiffness are also implemented in ABAQUS. Using the unified periodical boundary conditions in literature [59], the elastic moduli and Poisson’s ratio of the cellular metamaterials can be calculated. In our study, we are only interested in the longitudinal moduli E₁. The statistics of all samples’ longitudinal stiffness E₁ are shown in Fig. 3(b) and Table 1.

We first present a single-loop computational design approach, which will be used as the basic building block in the proposed iterative model updating-based design framework. The single-loop approach consists of three major components (top dotted box of Fig. 4):

1. A deep generative model (e.g., VAE, βVAE, or GM-βVAE) to extract low-dimensional latent feature vectors from the input images. The cellular metamaterial images are mapped to the points in the continuous latent space. New cellular structure designs can be generated by conducting arithmetic operations in the latent space (e.g., increasing or decreasing the values of the latent features, which is equivalent to moving the points in the latent space).

2. A GP model of the latent feature–property relationship for structure optimization. The properties of training samples are obtained by high-fidelity simulations. The GP model is trained on the simulation dataset to provide an efficient emulator of the latent feature–property relationship. Metamaterial structural design is formulated as a parametric optimization problem in the latent space:

   $maxzl{GP1(zl),…,GPM(zl)}Tsubjecttomin(z)≤zl≤max(z)$

   (7)

3. z_l represents the latent feature vector, l is the dimension of the latent features, GP_i represents the GP model for the ith mechanical property of interest, i = 1, 2, …, M. z represents the latent feature vectors of the training points. As the computational cost of design evaluation by GP is negligible compared with the true simulations, a variety of optimization algorithms, such as gradient-based algorithms and nongradient evolutionary algorithms, can be applied to solve the optimal z_l values. In this work, we adopt the Genetic Algorithm to find the optimal designs.

4. A deep decoder that maps the latent feature vector to the corresponding cellular structure image in the original input space. It is used to reconstruct the optimal cellular metamaterial design based on the optimal z_l values obtained in the previous step.

However, the single-loop design framework has a significant shortcoming. The deep generative model is trained on a randomly generated cellular structure dataset. This initial dataset may not fully cover the critical design regions where the global optimum is located. Therefore, the lack-of-data uncertainties in the critical design regions limit (i) the accuracy of the latent feature-structural image mapping and (ii) the accuracy of the latent feature–property supervised learning model, and thus will deteriorate the optimal search performances.

To resolve this issue, we propose an iterative GM-βVAE model updating-based design framework (Fig. 4). Similar to the strategy of active learning [60–62], this framework improves the model predictability in the critical design regions during the search iterations.

In each iteration, the single-loop design method is implemented to generate high-performance designs based on the existing dataset. The newly discovered designs are validated by simulations and used to replace the low-performance designs in the training dataset. Then the GM-βVAE model is re-trained for the next iteration. This process is iterated to achieve progressive improvements in design performances.

The proposed strategy is a performance-driven model updating process, which adds the newly discovered high-performance designs to the training dataset. Another model updating strategy is to reduce the epistemic uncertainty in the supervised learning model (uncertainty-driven model updating). New sample points will be generated at locations with the largest epistemic uncertainties to improve the predictability of the model. In Sec. 4.3, the performance-driven and the uncertainty-driven modeling updating strategies will be compared in a multiobjective design optimization problem.

The proposed approach has two stages: “offline model training” and “online design optimization.” The computational costs are mainly induced by the offline modeling training stage, which includes the property simulations for the initial training dataset, and the simulations of newly added samples during the iterative model updating process. A sufficient number of training samples are needed to guarantee the quality of the deep generative model and the GP model. We recommend to conduct a convergence study to find the minimum number of samples.

A deep generative model can extract the low-dimensional latent features from the input image data and reconstruct images based on the latent feature values. Here we provide a review of three VAE-based deep generative models and present a comparative study that helps to select the best model for structure design.

Variational autoencoder, which was first proposed by Kingma and Welling [63], is one of the most popular methods to learn complicated data distribution such as images in an unsupervised manner. With embedded Bayesian inference in the autoencoder architecture, a VAE learns the mapping relationship between the input image x and the corresponding latent space z.

Variational autoencoder consists of an encoder (recognition model) and a decoder (generative model). The encoder of the VAE performs nonlinear dimensionality reduction and compresses the high-dimensional data into a low-dimensional latent space, which can be expressed as $Qϕ(z|x)$⁠. $Qϕ(z|x)$ is the approximate posterior that follows a normal distribution. ϕ is the vector of the recognition model parameters.

VAE/βVAE shows a great advantage in nonlinear mapping, However, the prior distribution of the latent variable is assumed to be a standard normal distribution, which may not perform well on the dataset that shows multimode characteristics.

Different from the typical VAE and βVAE, the prior over each latent variable of GM-βVAE is modeled as a mixture of Gaussian, which can capture the complex statistical distribution with multimodal characteristics. We assume that the input data are generated from different classes, which means that each mode of the latent feature corresponds to a specific class. Each class of data obeys a Gaussian distribution. The departure of each class of data is measured by Mahalanobis distance [65].

The learning capacity of a deep generative model is assessed by measuring the difference between the original image and the reconstructed image. Here, we demonstrate the advantage of GM-βVAE over VAE and βVAE in terms of reconstruction accuracy. For a fair comparison, the same deep learning network architecture and the same training dataset are used in all models.

The details of the generative network models are shown in Fig. 5. The dimensionality of the latent space is set to be four for both the VAEs and GM-βVAE model. Based on the results of test runs, we find that using a lower dimensionality will incur a large reconstruction loss while using a dimensionality higher than four will not lead to a significant improvement in the reconstruction accuracy. The differences between the two networks lie in the definitions of the loss function and the distribution functions of the latent features (Gaussian distribution versus Gaussian mixture distribution).

We created a 2D cellular metamaterial database with 400 samples collected from literature [66,67] and our previous research [68]. We divide our dataset into a training set with 360 samples and a testing set with 40 samples randomly. Model training is performed for 150 epochs with a batch size of 32 for all three models. We used Adam as the optimizer to train the model. A visual comparison is conducted on the reconstructed images. The VAE, βVAE, and GM-βVAE models all have a good reconstruction quality and preserve most of the geometrical features (Fig. 6). However, some of the reconstructed images have blurred local regions, such as the central region in the image and the structural boundaries. The blurred regions are common in the reconstructions by VAE because of their log-likelihood loss functions [69]. We tackle this issue by setting a threshold to binarize the reconstructed images. We set the binarization threshold as 0.9. Pixels higher than 0.9 are converted to solid, and pixels lower than 0.9 are converted to void. The binarization step induces minimal changes to the structural geometry. We use the binarized structures in our finite element simulations, which has been introduced in Sec. 3.2.3.

The purpose of unsupervised generative design is to discover a diverse set of cellular metamaterial designs that span a wide range of properties (e.g., bandgap width), without conducting the expensive design evaluations (e.g., phononic property simulations). This design set provides a start point for the following supervised learning-based design optimization. First, we implement a traditional approach of unsupervised generative design. With a trained βVAE/GM-βVAE model, random samples are collected by sampling the latent space using space-filling methods (e.g., Design of Experiments) [60]. The assumption of this approach is that a wide range of design performances can be achieved by filling the latent space with stochastically but evenly distributed sample points.

Here a case study is presented to illustrate the mechanism of generating diverse structural designs using a deep generative model. The basic idea is to morph the structure design by manipulating the values of the latent feature variables.

First, we choose the cellular structure whose latent features are [0, 0, 0, 0] and vary each latent variable value along the negative and positive directions. When we vary the value of a latent feature, the other three latent variables are fixed to 0. The resultant new structures are shown in Fig. 7(a). It is obvious that different latent variables control different topological patterns in the structural image. As the latent variables are continuous, we can generate an infinite number of structural patterns theoretically.

Second, we investigate the interpolation and extrapolation between two existing samples in the training dataset (highlighted by squares in Fig. 7(b)). New designs are generated by conducting linear interpolation/extrapolation based on the two existing sample points in the latent space. A continuous evolution in structural morphology can be obtained along this linear path.

We propose an unsupervised generative design approach by leveraging the multimode nature of the GM-βVAE latent feature distributions. The statistical distribution of each latent feature is a mixture of Gaussian distributions. The mean of each Gaussian distribution (mode center) is considered as the center of a cluster of samples with similar characteristics. Therefore, a “cluster exemplar” can be obtained by constructing a latent feature vector in which each feature takes the value of a mode center (Fig. 8(a)). There are four latent features and two distinctive modes are identified for each latent feature, so in total 2⁴= 16 combinations (designs) are obtained. The corresponding cellular metamaterial images of the sixteen designs are obtained by decoding, as shown in Fig. 8(b).

The effectiveness of the unsupervised generative design is evaluated by comparing the range of properties achieved by the sample set. A wider range indicates a higher diversity of the design set. We also generate 16 designs using the traditional LHS approach. The property of each design (bandgap width) is obtained by simulation. It is observed that the proposed generative design approach achieves a much wider property range (Fig. 8(c)).

The purpose of the second case study is to showcase the single-loop approach in property-driven structure design optimization. The design goal is to find cellular metamaterial designs with the maximal bandgap width. We repeat the GP modeling-optimization process three times independently and get three different optimal candidates, as shown in Fig. 9. The total bandgap width values of the three optimal candidates are verified by phononic simulations in ABAQUS. There are two important observations.

1. The single-objective optimization results in unrealistic designs for engineering applications, as shown in Fig. 9(a). The structural designs consist of several solid components connected by thin beams, which will lead to low stiffness and loss of structural integrity under loadings.

2. We also found that the performances of the training data could affect the performances of the optimization results. We regroup our existing dataset into two groups: the first group contains 320 designs with the top 80% width of bandgap width, the second group contains 320 designs with the bottom 80% width of bandgap width. The bandgap widths of the optimal designs found by the first model (trained on the high-performance data) are one order of magnitude higher than the optimal designs found by the second model (trained on the low-performance data). Therefore, we hypothesize that improving the performance of the VAE training dataset will lead to better designs.

A multiobjective design problem is formulated to maximize the total bandgap width and to maximize the structural stiffness (longitudinal stiffness E₁). The multiobjective design formulation can resolve the issue of lacking stiffness in the single-objective optimization. Two surrogate modeling strategies are compared in this study: establishing a multi-response GP model [70] that predicts both properties simultaneously, and establishing independent GP models for each response. The metamaterial samples are divided into the training dataset (364 points) and the test dataset (36 points) randomly. Multi-response GP yields a much lower accuracy (Fig. 10), because the two desired mechanical properties (width of total bandgap and longitudinal stiffness) are weakly correlated (the correlation coefficient is less than 1e-8). Therefore, we decide to train two independent GP models for each response for the following optimization studies. Nondominated Sorting Genetic Algorithm (NSGA-II) [71] is applied to find the nondominated design set.

In the first iteration, 20 nondominated designs are obtained and the corresponding design images are reconstructed by decoding. The top 10 designs with the largest bandgap width are used to replace the top 10 undesired designs with a small total bandgap width and a small elastic modulus in the training set. The GM-βVAE is re-trained with the updated dataset for the next iteration.

In total three iterations are conducted. The training dataset and the GM-βVAE model are updated once in each iteration. The nondominated sets of the three iterations are shown in Fig. 11(a). The design performances in Fig. 11(a) are predicted by the GP model. The true performances of the nondominated designs, which are obtained by finite element simulations, are shown in Fig. 11(b). Three representative designs, the two ends (#1 and #20) and the mid-point (#11) of the nondominated set (Figs. 11(a) and 11(d)), are selected and compared in Table 3. We find that:

• The nondominated design set moves closer to the utopian point (upper right corner) during the iterative design process.

• Most of the designs on the Pareto frontier (marked by the dashed line in Fig. 11(b)) are from the last iteration, which indicates the effectiveness of the iterative search process.

We compared our designs with the optimal phononic metamaterial designs reported in literature [17,21,22] (Figs. 11(c) and 11(e)). These designs are obtained by topology optimization. We simulate the stiffness and total bandgap width of these designs and compare them with the nondominated designs found by our proposed design framework. As shown in Fig. 11(c), our proposed design framework successfully generates designs that excel the literature designs in both stiffness and total bandgap width.

We also investigated another iterative model updating strategy that focuses on reducing the epistemic uncertainty in the supervised learning model (uncertainty-driven model updating). In the first two iterations, the GM-βVAE and GP models are updated by adding 15 new sample points collected from design regions with the highest epistemic uncertainties to the training set. In the last iteration, the design objective is formulated as maximizing the total bandgap width and the structural stiffness (longitudinal stiffness E₁). We compare the nondominated design set found by the uncertainty-driven model updating strategy and the proposed property-driven model updating strategy. As shown in Fig. 11(c), the proposed property-driven model updating strategy performs better than the uncertainty-driven model updating strategy.

It is to be noted that there exists a large difference between the design performances predicted by GP and the true values obtained by finite element simulation (Figs. 11(a) and 11(b), Table 3). The prediction errors are induced by the following factors.

• Small changes are made on the reconstructed design images in postprocessing, which includes eliminating the small disconnected solid parts in the design [72–74] and converting the blurry image to a binary image.

• The accuracy of the GP model is limited due to the limited amount of data.

Although the design performances predicted by GP cannot fully match the true values obtained by finite element simulation, a strong correlation is observed between them (the R squared value is 0.88), as shown in Fig. 13 in the  Appendix. It indicates that the “high-performance” design points found by the surrogate model have a high probability to possess a real high performance.

In this paper, we present a deep feature learning-based framework to design cellular metamaterials for bandgap width and stiffness. By taking advantage of the latent space learned by GM-βVAE, this design framework enables both unsupervised generative design and property-driven optimization. Our major conclusions are summarized as follows:

1. The proposed deep feature learning-based framework is effective in designing metamaterials with complex topological characteristics.

2. The proposed iterative updating-based design framework is effective in improving the design performances.

3. The GM-βVAE-based unsupervised generative design method can generate a more diversified design set than the traditional space filling-based unsupervised generative design method.

4. As the two responses are correlated weakly, it is recommended to establish independent GP models for each response (stiffness and total bandgap width), instead of an all-in-one multi-response GP, for supervised learning and design optimization.

Two caveats to designing phononic bandgap metamaterials are identified.

1. When considering only the phononic bandgap as the design objective, the optimization algorithm tends to generate structure designs with weakly connected components. Stiffness is an indispensable design objective/constraint to guarantee structural integrity under mechanical loads.

2. Postprocessing (e.g., filtering) on the reconstructed structure image may lead to significant changes in the bandgap properties, which indicates a lack of robustness in our proposed designs.

We also identified three opportunities to further improve the proposed framework.

1. The proposed design framework will be further improved and extended to the design of 3D metamaterials.

2. To resolve the issue that small changes in the optimal design’s geometry lead to a large degradation in phononic properties, robust optimization methods will be investigated.

3. The diversity of the design set can be improved in two ways. First, the VAE model can be improved by incorporating a diversity-oriented loss function that encourages the generation of novel designs. The diversity-oriented loss function was proposed based on the determinantal point processes [75], which provides a measurement of the design similarity. Second, the design diversity can also be improved in the stage of multiobjective optimization. By integrating the niching method with the Evolutionary optimization algorithms (e.g., niching Genetic Algorithm [76]), multiple optima at different locations in the design space can be obtained in the property-driven optimization.

This research is supported by Ford University Research Program. Z.W. acknowledges the support from the University of Connecticut School of Engineering General Electric (GE) Graduate Fellowship for Innovation. Y.L. gratefully acknowledges financial support from the Air Force Office of Scientific Research through the Air Force's Young Investigator Research Program (FA9550-20-1-0183; Program Manager: Dr. Ming-Jen Pan) and the National Science Foundation (CMMI-1934829). Y.L. would like to thank the support from 3M’s Non-Tenured Faculty Award. This research also benefited in part from the computational resources and staff contributions provided by the Booth Engineering Center for Advanced Technology (BECAT) at UConn. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the U.S. Department of Defense.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

We trained multiple GM-βVAE models with different numbers of modes (1, 2, 3, and 4) in the Gaussian mixture latent space for each design iteration. A convergence study is conducted to find the optimal number of modes that lead to the highest fitting accuracy. When the number of modes equals 1, the latent variable follows the normal distribution. We compared the cumulated distribution function (CDF) of each latent variable fitted by Gaussian mixture with the empirical CDF of the true latent variable values. Figure 12(a) shows the fitted CDF versus the true CDF of the four latent features when the number of modes equals 2. The fitting error is evaluated as the area enclosed by two CDF curves. Figure 12(b) shows the result of the convergence study. Increasing the number of modes will reduce the enclosed area between the true CDF and the fitted CDF (smaller error). We observe a significant improvement by increasing the number of modes increases 1 to 2, but a slight improvement from 2 to 3 or from 3 to 4. So we choose a Gaussian mixture with 2 modes in our GM-βVAE model.

For the second and third iteration of the design process, the trends are the same. So we use Gaussian mixture with two models for all three iterations.