Inverting the structure–property map of truss metamaterials by deep learning

Significance More than a decade of research has been devoted to leveraging the rich mechanical playground of periodically assembled truss metamaterials. The enormous design space of manufacturable unit cells, however, has made the inverse design a challenge: How does one efficiently identify a complex truss that has given target properties? We answer this question by a data-driven method, which instantly (once trained, within milliseconds) generates not one but a variety of truss unit cells, whose effective response closely matches a given (fully anisotropic) stiffness tensor. Moreover, our framework to smoothly transition between different unit cells enables the design of lightweight structures with spatially varying, locally optimized properties, for applications from wave guiding to artificial bone.


Details on the lattice design space and dataset generation
Fig. S1 shows the seven considered elementary topologies in the proposed lattice design space. Each row of matrix X corresponds to the 3D positional coordinates of a truss junction (or node) in the UC, while matrix C defines the connectivities between truss junctions (being 1-indexed, each row defines a connection between the nodes with the given two indices).
The full set of parameters describing the design space of our truss lattices is summarized in Table S1. A dataset of 3 million lattices with their corresponding homogenized stiffness response was generated by drawing with a (discrete) uniform distribution within each of the provided ranges; e.g., the relative density is randomly sampled as ρ ∼ U (0.002, 0.1). For the rotations, the angle-axis representation was used, where an angle θ ∼ U (−π, π) and a random unit vector (as the axis of rotation) on the three-dimensional unit sphere S 2 = {k ∈ R 3 : k = 1} was drawn.

Details on computational homogenization and stiffness visualization
To extract the homogenized stiffness tensor components C ijkl from a UC in our design space, we used the finite element method. We modeled each strut in the UC as a linear elastic Timoshenko beam element with circular cross-section. Any beam intersections, arising for certain combinations of elementary cells and tessellations, were resolved by introducing a new vertex at each intersection point and subsequently splitting the affected struts into adjoining sub-struts. We imposed periodic boundary conditions onto the boundary nodes of the UC (i.e., periodic displacement and anti-periodic force boundary conditions applied to all nodes on the outer UC boundary, while boundary rotations were set to be periodic). The effective stiffness tensor components C ijkl were extracted from the resulting assembled UC stiffness matrix, following (1).
To visualize the elastic surfaces, we plot the effective Young's modulus along a direction d ∈ S 2 (S 2 = {k ∈ R 3 : k = 1} denoting the unit sphere in 3D) as [1]

Machine learning framework
A. Data preparation. The continuous (geometrical) features, consisting of the relative density ρ and the six stretch tensor eigenvalues U1, U2, U3 and V1, V2, V3, were rescaled with a min-max normalization to the range [0, 1], i.e., The categorical (topological) features, consisting of the first six parameters in Table S1, were one-hot encoded to prevent an ordinal interpretation of the different topologies, while the rotations RI , RII (which form a part of the geometrical features) were left unmodified, as they directly enter the manually computed tensor rotations. The stiffness labels were consistently scaled with a min-max normalization to the range [−1, 1], i.e., where Ci stands for each of the 21 independent stiffness constants, and each Ci was treated independently to enforce equal contributions to the loss of the 9 (for F1) or 21 (for F2 and G1, G2) given stiffness parameters during training. Each stiffness parameter is then scaled to its original value via Ci ← 0.5 (Ci + 1) [max(Ci) − min(Ci)] + min(Ci), [4] before applying rotations at the intermediate stages within the forward model to correctly evaluate the rotation. Table S1. Overview of the lattice design space described by both discrete and continuous design parameters summarized in vector Θ. Rotations R I and R II can each be represented by three parameters such as, e.g., via the angle-axis representation, which was used here to generate the dataset. Note that each elementary topology is used as a single tessellation (t i = 1) or a 2 × 2 × 2 tesselation (t i = 2), which is binary-encoded. B. Details of the inverse model and tensor rotations. The inverse model has different activation functions applied to its outputs in the final layer, depending on the type of design parameter the output is corresponding to. For the one-hot encoded discrete predictions (relating to the topology of the proposed lattice), the Gumbel-softmax trick (corresponding to Eqs. (5) and (6) of the main article) is applied. For the geometrical predictions relating to ρ, U1, U2, U3 and V1, V2, V3, the sigmoid activation function sig(Θ) = 1 1 + exp (−Θ) [5] is applied to ensure that predicted values are within the same (min-max-normalized) range of values provided within our training data. To obtain the two rotation matrices RI and RII , we leverage that any rotation matrix R ∈ SO(3) can be constructed in a Gram-Schmidt-like process via [6] where N (·) represents a normalization function, and a1, a2 ∈ R 3 correspond to the so-called 6D representation (2). These six parameters are directly assembled by the output of our inverse model (i.e., no activation function is applied) and, since we sample two rotation matrices, this results in a total of 12 additional output dimensions characterizing the rotations. Most importantly, it can be shown that this rotation parameterization is a continuous representation of the space of SO(3) rotations (unlike Euler angles or quaternions), which facilitates the training of NNs (2). Although we did not conduct a rigorous comparison of the NN performance for different rotation representations, results indicated a small improvement in the accuracy of the inverse prediction using the above 6D parameterization compared to other representations. All stiffness tensor rotations were computed by using Voigt notation, i.e., considering C in Voigt notation and applying the transformation

Range
where , [9] where Rij refer to the components of the second-order rotation tensor. This formulation bypasses transformations between the Voigt stiffness matrix and the fourth-order stiffness tensor and was first described by Bond (3). Predicted geometrical design features Θ pred i were rescaled to the corresponding physical ranges (see Table S1), e.g., for the relative density via ρ pred ← ρ pred (0.1 − 0.002) + 0.002, [10] during post-processing and FEM reconstruction.

C. Protocols for NN training.
Details of the architectures and other training parameters for both NNs of the forward model (F1 and F2) and the inverse model (G1 and G2 share the same hyperparameters) are provided in Table S2. We used 1% of the generated dataset for the optimization of the presented hyperparameters. Hidden layers perform a linear transformation followed by a leaky rectified linear unit (leaky ReLU) as an activation function. The PyTorch package (4) was used throughout our implementation, which also allows for automatic differentiation and thus straightforward evaluation of the required gradients of the tensor rotations (Eq. (7)) during training via its automatic differentiation engine, autograd. Note that the learning rate was dynamically reduced during training of the inverse model, using PyTorch's ReduceLROnPlateau scheduler with a patience of 20 and reduction factor of 0.5. Network trainings were performed on a single Nvidia Quadro RTX 6000 with 24 GB GDDR6 memory. Fig. S2 show the loss of all considered NNs at every training epoch (the two NNs of the forward model, and the inverse model whose NNs were trained jointly). Train and test loss behave similarly, indicating no significant overfitting. Table S2. Input and output dimensions and hyperparameters for the optimized forward and inverse models. The feature scaling for the second NN of the forward model, F 2 , is only applied to the stretch eigenvalues V 1 , V 2 , V 3 ; the rotated stiffness given as an input is left unchanged.
Feature and label scalings are only applied to continuous design parameters (excluding rotations). 1 Learning rate decay was applied using PyTorch's ReduceLROnPlateauscheduler with a patience of 20 and reduction factor of 0.5.

D. Evaluation of accuracy.
In the main article, the C1111 reconstruction accuracy of both the forward and inverse model was shown ( Figure 2). The full reconstruction accuracy of all stiffness components (including R 2 -values) for the forward and inverse models is summarized in Figures S3 and Fig. S4, respectively. All metrics of accuracy are evaluated on a separate test dataset (previously unseen during training) comprised of 30,000 lattices and their corresponding stiffness components.  here since it is indicated in the main article). The corresponding R 2 -deviations are indicated. Note that due to the stochastic nature of our inverse model, R 2 -deviations can vary for different drawings; however, these differences are negligible for the given test set. Different colors indicate the different coupling terms of the elasticity tensor (i.e., pure extension/compression: blue, normal/normal: violet, normal/shear: magenta, pure shear: green, shear/shear: pink).

E. Variational sampling.
To showcase the generative nature of our inverse variational model, Figure S5 presents the obtained inverse predictions for four randomly selected lattices in the test set (labeled A, B, C, and D). Except for lattice B, we obtain different topologies for a different second seed. Note that we can, however, trade off some reconstruction accuracy for a larger topology variation by increasing τ in our softmax layer (see Eq. (5) of the main article), as this increasingly flattens the posterior probability distribution of the predicted topologies. While setting τ too large will drastically deteriorate the results, we noticed that a moderate raise does increase the variance without a notable loss in accuracy. The design parameters for each of the cases A, B, C, and D are given, respectively, in Tables S3, S4, S5, and S6.
(same prediction as with seed 1)

F. Effect of the dataset size.
Large datasets (such as the chosen one containing 3 million lattices and their corresponding stiffnesses) may be expensive to generate in other scenarios, which is why we assessed the performance of our framework on a smaller dataset as well. We re-trained the networks using only 10% of the lattice-stiffness pairs of the original dataset and provide the resulting loss plots in Figure S6. For the forward model, the decrease in the corresponding coefficient of determination R 2 is between 0.004 (for the more dominant terms such as C1111) and 0.031 (for the less dominant terms such as C1312), as shown in detail in Table S7. The inverse model shows a similar behavior, the decrease in R 2 ranges from around 0.01 (for C1111) to 0.054 (for C1212). As expected, the models perform slightly worse than those trained on the significantly larger dataset. However, the R 2 values indicate that the accuracy has not deteriorated drastically. We therefore conclude that the framework still performs comparably well, if only a fraction of the data is provided (which is of importance when data generation is a limiting factor). Note that, as a benefit of our physics-informed approach, data augmentation by adding random rotations of the provided lattices can easily extend the dataset for inverse training without requiring additional FEM simulations (but this was not performed in this work).  G. Performance of the inverse model compared to the best match in the training data. As the dataset we used to train the networks contains 3 million lattice-stiffness pairs, it may seem like such a large dataset should directly contain (close to) optimal designs for a given stiffness. However, the curse of dimensionality renders these 3 million samples negligible to the total number of possible combinations. Consider a very coarse sampling of the design space, e.g., taking only 10 different numerical choices for each of the 13 continuous (geometric) design parameters. Given the 262 chosen truss topologies, this gives 262 × 10 13 possible combinations. Generating such large datasets and comparing them to a given target stiffness is computationally expensive -prohibitively expensive for applications requiring repeated inverse predictions (such as, e.g., for two-scale topology optimization). To be specific, finding the closest match in our given database for a given stiffness using the Pandas library for Python (6) took on average ca. 35s (performed on an Intel Xeon Silver 4114 2.2GHz processor), whereas our inverse model predicts a design within ca. 2ms (30,000 stiffnesses are evaluated in about 60s on the same machine). This four-orders-of-magnitude difference in efficiency is a clear advantage over the purely data-driven lookup tables approach. Besides, even if a dataset of 3 million entries could be sufficient -how would one choose the 3 million samples from the (at least) 262 × 10 13 possible combinations to cover a sufficient stiffness space? Our approach, which interpolates (in the latent space) between the chosen 3 million samples, indeed provides a sufficiently rich stiffness space. Notably, we have shown in Section 3.F that the use of a reduced dataset for training (e.g., only 10% of the original dataset of 3 million lattice-stiffness pairs) still achieves good accuracy in inverse design prediction. In addition, let us demonstrate that our framework generates better designs as compared to the closest match from the training dataset. To give a specific example, we identified the closest match in our database for the hex_Z04._R145 sample from the crystallographic periodic network database (7) (presented in Figure S7a). This architecture is not contained in our design space (it is not even close to the chosen elementary topologies). We identified the best match from our training set by minimizing the squared error for all 21 elastic components, and in Figure S7 we compare it to the reported prediction of our inverse model. For this example, the inverse-model prediction is much closer to the given target stiffness (normalized mean squared error of 0.003) than the best match from the training dataset (normalized mean squared error of 0.036).

H. Computational efficiency of the inverse model.
To estimate the efficiency of our approach, we have added a detailed overview of the computational runtimes required for the training and use of our framework in Table S8. After training is done once a priori/offline, subsequent evaluations of the inverse model take a few milliseconds (2ms on average for 30,000 inverse predictions) because no subsequent solution of a complex boundary value problem or an iterative evolution simulation is required-as, e.g., in classical topological optimization and evolutionary algorithms. If one seeks only a single optimal design, then those approaches may indeed be competitive in terms of computational costs. But if one repeatedly needs to evaluate many optimal designs-e.g., when performing multiscale topology optimization (8)-then this makes a dramatic difference. Note that we can easily even reduce the training time if we consider less training data without sacrificing much accuracy (see Table S7) or consider earlier stopping for the inverse model training (since the loss decreases only marginally after around epoch 100, as shown in Figure S6). Table S8. Runtime, software, and hardware resources used for different tasks (including in parentheses the runtimes for a reduced dataset consisting of 300,000 lattice-stiffness pairs, see also Section 3.F). 2 Reported runtimes are rough estimates only. 3 Tasks were performed on a 2.2 GHz Intel Xeon E5-2650 processor with 256 GB of DDR3 memory at 2.5 GHz. 4 Tasks were performed on a single Nvidia Quadro RTX 6000 with 24 GB GDDR6 memory using CUDA 11. 5 Runtimes for predictions via the NNs are reported for one sample (averaged over a batch size of 1% of the full training data). Note that the averaged runtime partly depends on the size of the evaluated dataset due to the computational overhead.

Evaluation of the inverse model on crystallographic truss networks and bone samples
The 21 anisotropic stiffness components of the homogenized responses of the truss lattices computed by Lumpe and Stankovic (7) and of the bone samples characterized by Colabella et al. (9), including their inverse-designed counterparts, are presented in Tables S9 and S11, respectively (using Voigt notation for the stiffness components in 3D). As mentioned in the main article, we can leverage the variational nature of our framework to propose a large variety of lattice candidates and subsequently select the ones with the lowest NMSE, which can be evaluated (at negligible computational cost) by the forward model. For this study, we set τ = 100 in the Gumbel-softmax layer to enforce a large exploration of proposed lattices and selected the best predictions from 10,000 generated samples, which only required around one minute of overall runtime. The inverse model predicts a lattice unit cell based on a normalized stiffness input, i.e., C = C b /E b for a lattice with the actual homogenized stiffness response C b manufactured with an isotropic base material having Young's modulus E b . Therefore, we normalize the target stiffness by E b . As reported in the main article, we consider Ti-6Al-4V with E b = 114 GPa (10) for the application in bone implants due to its excellent biocompatibility and corrosion resistance (11), and we hence normalize the reported bone stiffness values accordingly before passing them into the inverse model. Note that, in principle, any E b can be considered, as long as the given target stiffness can be obtained by lattices within the given limits of relative density ρ.
As the stiffnesses provided by Lumpe and Stankovic (7) were already normalized by the base material's Young's modulus, no further pre-processing was required. Nevertheless, we applied a factor of 20 to the provided stiffness of hex_Z04.0_R145, as the original stiffness is not obtainable within the density range of our design space (as it would require relative densities below our chosen lower bound of ρ = 0.002, which can easily be assessed as the corresponding inverse predictions approach the bounds of our design space). Alternatively, one could also retrain the model with an even lower bound on the relative density and subsequently omit the scaling.
To give further quantitative metrics for the performance of our framework, we have applied it to the anisotropic stiffnesses of the complete catalog of crystallographic period networks reported in (7), consisting of 17,262 vastly different truss topologies (with notable shares of highly symmetric (i.e., cubic) to less symmetric (i.e., monoclinic) structures). As before, we have uniformly scaled the stiffnesses by a factor of 10, since some of the provided structures are highly compliant and would otherwise require relative densities below the lower limit of our design space. The results are shown in the histogram in Figure S8. The vast majority of predictions shows an NMSE comparable to the predictions shown in Figure 3 of the main article, highlighting that they do give an accurate representation of the performance of our inverse model on arbitrary stiffnesses. Table S9. Elastic constants (using Voigt notation) of crystallography-inspired periodic truss networks and the corresponding elastic constants of the lattices obtained by our inverse design framework, visualized in Figure 3a of the main article. The design parameters of the inverse predictions are listed in Table S10.    Table S11. Elastic constants (using Voigt notation) of trabecular bone in bovine femurs and the corresponding elastic constants of the lattices obtained by our inverse design framework, visualized in Figure 3a of the main article. The design parameters of the inverse predictions are listed in Table S12.   different topologies) provided in (7). Different percentiles are color-coded to better interpret the evaluated NMSEs.

Generation of functionally graded lattices
For the integration of different truss UCs in a functionally graded lattice with spatially varying truss topology and geometry, we must, on the one hand, transition between different UC topologies (which here were chosen to all be based on a cubic UC).
On the other hand, we must transition between differently shaped UCs arising from the application of affine transformations and rotations. In practice, we choose a total of k control points at locations {x1, . . . , x k }, at which we define a local truss UC (along with its topology and geometry), so that we seek to smoothly transition between those. All lattices are by definition periodic, and the domain Ω UC of each UC is spanned by three translation vectors forming the basis A = {a1, a2, a3} in 3D, which defines the tessellation of the UC into the truss lattice. Let us denote by K = {k1, k2, k3} the triad of vectors forming the basis of the corresponding reciprocal lattice.
Any periodic (simple Bravais) lattice is represented equally by the so-called impulse train and its Fourier representation: where δ denotes the Dirac delta function, and M > 0 is an integer. In a spatially variant lattice, the above Bravais basis vectors vary from point to point, i.e., we assume a spatially varying basis A(x) = {a1(x), a2(x), a3(x)} and consequently a spatially varying reciprocal basis K(x) = {k1(x), k2(x), k3(x)}. Note that we tacitly assume a separation of scales between the microscale (i.e., the scale of the truss UC dimensions) and the macroscale (i.e., the scale of the overall body to be filled with the spatially graded truss), which allows us to assign a unique basis A (and hence K) to each point x in the macroscale body. Of course, in reality such a separation of scales may not be realizable, yet the resulting framework serves as an effective approximation whenever the lattice bases vary sufficiently smoothly across the macroscale body.
To realize a spatially variant truss lattice, we aim to construct a (non-periodic) truss with smoothly varying topology such as to approximate A(x) locally. To this end, we relax the above description of the Fourier-transformed impulse train by introducing three projection functions φi(x), one for each of the three ki-vectors in K. Further, we note that the translational symmetry is contained within the first spatial harmonic, which is why we simplify the Fourier decomposition by truncating after its first spatial harmonic (M = 1), resulting in the approximation 1 + cos(2πµφi x) with φi(x) = arg inf Ω ||∇φi(x) − ki|| 2 dΩ. [13] φi(x) is hence a smoothly varying field that is locally tangential to ki, such that the set {φ1(x), φ2(x), φ3(x)} presents an approximate, smooth, spatially varying basis field over Ω. Now, consider a desired basis vector field A(x) (which, e.g., locally matches the translation vectors {a1, a2, a3} of the truss UCs obtained from the inverse model applied to different properties at different control points xi on a body Ω). In practice, we define a radial basis function N (x, xi) for each control point xi, centered at control point xi, as [14] ζ > 0 controls the thickness of the transition layer between different topologies. If A k is the local basis of the UC translation vectors at control point k, then we interpolate to obtain Ai N (x, xi). [15] For example, in Figure 3b in the main article, we introduced two control points on the x1-axis at x1 = (xmin, 0, 0) T and x2 = (xmax, 0, 0) T , and we chose ζ = 0.5. Given A(x), we numerically solve the least-squares problem in Eq. (13) for the projection functions φi(x), using the finite element method. Subsequently, we construct the relaxed lattice function τ (x) and choose all N local minima of τ (x) in the domain Ω as the set of the UC centroid positions, B Ω UC = {x1, . . . , xN }. Finally, we assemble the global lattice by placing a UC at each xν ∈ B Ω UC with basis A(xν ). Note that the corner nodes of the resulting adjacent UCs do not perfectly overlap. We therefore apply a smoothing step, which replaces all almost overlapping nodes by a single new node at the arithmetic mean position of the original nodes, and we adjust the beam connectivities accordingly.
While the above strategy takes care of spatially varying translation vectors, we must also transition between the different UC topologies at the k control points. To this end, we superimpose all k UC topologies and weight the diameter d s (x) (where s = 1, ..., m) of each of the m individual struts in the superimposed UC by a linear superposition of the constituent diameters by re-using radial basis functions Eq. (14), such that [16] Here, d s i corresponds to the respective strut diameter in the UC at control point xi (where d s i = di if the strut is included in the topology with diameter di and d s i = 0 if it is not). In Figure 3b in the main article, the two control points at x1 = (xmin, 0, 0) T and x2 = (xmax, 0, 0) T reduce Eq. (16) to and we again chose ζ = 0.5.