NewSite

机械超材料执行器的自动设计

摘要

机械超材料执行器通过利用单个3D打印元件中编码的架构特征，实现预定的输入-输出操作，从而无需组装不同的结构部件。尽管该领域发展迅速，但仍缺乏高效策略来优化多种功能的超材料设计。我们提出一种结合强化蒙特卡洛方法与离散元模拟的计算方法，用于自动设计机械超材料执行器。对选定机械超材料执行器的3D打印结果表明，机器生成的结构可实现高效率，超越人工设计的结构。我们还证明，通过训练深度神经网络能够从结构图像中预测效率并识别其功能区域，从而设计高效执行器。本文设计的单元执行器可组合成任意复杂度的超材料机器，满足无数工程应用需求。

引言

机械超材料是一类新型人工材料，其设计旨在获得常规材料难以具备的特殊性能与响应特性1-5。其刚度、强度重量比6-8、弹性响应9、泊松比10-14、能量捕获能力15,16及抗断裂性17,18均可调控至匹配甚至超越常规材料的水平。数字制造技术（如3D打印和自动化装配）的最新进展，使这类材料结构能够以更低且持续降低的成本快速制造，突破了许多尺度与几何形状的限制19，从而推动了对超材料的日益关注。

超材料的性能并非源于块体材料的固有特性，而是源于其人工设计的内部几何结构——通常由多个子单元（或称晶胞）按重复规则排列而成20,21。近期研究还表明，在机械超材料中引入更高程度的结构无序性而不降低效能是可行的22-25。由于晶胞可通过多种不同方式设计与排布，所得结构可呈现多种自由度，从而产生各种非常规物理特性，这些特性在工业设计、建筑装饰、纺织品及梁体等结构的增强图案中具有天然应用前景。

基于折纸与剪纸的形状变化结构是值得关注的范例，它们利用折叠与切割实现形状编程26-29，在航空航天（如卫星太阳能板展开30）等领域具有多种工程应用。基于折纸31-33和剪纸34的结构为机械超材料设计提供了灵感。在此背景下，近期有研究探索了具有可控多功能性的棱柱状构件的折叠行为，并通过应用局部驱动模式来研究其力学性能35。

超材料亦可视为真正的机械系统36，能够将输入激励转化为一系列可编程输出，从而完成机械功能。此时，组成单元以明确界定的方式协同工作，最终实现可控的定向宏观运动。超材料机械系统可用作人机交互的机械执行器，或作为机器人中的交互/响应组件。传统超材料结构与机械系统的设计策略通常基于人工操作，在特定条件下表现尚可，但无法保证在所有情况下均实现最高效率。

为克服人工设计的局限性，我们提出一种用于机械超材料执行器自动设计的优化方法，通过对结构进行迭代修改，寻找对外加输入的最佳输出响应。基于优化算法的生成式方法此前已被用于设计变构材料37-39、纤维增强执行器40，以及在周期性超材料晶格中选择最优晶胞几何形状41,42。而本方法面向有序或无序22,25超材料，通过梁单元的移除与重新插入进行迭代修改，从而选择对外加给定输入的最优输出响应。研究表明，该优化过程可通过将优化算法与离散元模拟或适当训练的深度神经网络相结合而高效实现。

结论

位移输入/输出下的自动设计。 我们考虑两种典型的执行器，第一种中期望的输入和输出是正交的， $\mathbf{t}_{\text{inp}} = -\hat{y}$ ， $\mathbf{t}_{\text{out}} = -\hat{x}$ （见图1），第二种中它们是反平行的，

$\mathbf{t}_{\text{inp}} = -\hat{y}$ ， $\mathbf{t}_{\text{out}} = \hat{y}$ 。在补充图1中，我们报告了正交功能化MMA在蒙特卡洛（MC）动力学过程中 $\eta$ 的演化（反平行情况的轨迹类似）。优化步骤中生成此类结构的一个例子可在补充影片1中可视化。我们注意到，在退火阶段之后，效率倾向于分阶段演化，穿插着噪声部分和平台期。当触发期望响应的机制最终被激活时会出现阶跃式变化，而平台期则表明该结构的效率对一条或多条键的移除或添加具有鲁棒性。因此，我们从这些 $\eta$ 的平台期对应的构型中选取了效率最高的参考样本。我们得到了最终 $\eta$ 值的大幅变化，表明可能存在需要热激发才能逃逸的局部极小值。因此，并且取决于后退火条件，整个网络的相空间探索可能需要很长时间，这会进一步增加网络的尺寸。

在最小化算法中使用简化模型需要通过更精细的模拟进行进一步验证。因此，我们将我们的结构转换为有限元（FEM）网格，以模拟材料的真实响应（见“方法”）。正交和反平行运动的模拟结果分别报告在图2和补充图3中。请注意，可实现的效率可以轻松接近并超过人类设计结构的性能（左面板报告）。在所有考虑的情况下，模型计算的效率都通过FEM模拟进行了验证，相对于离散元（DEM）估计， $\eta$ 的修正很小。从FEM或DEM模拟中获取网络内的应力传播信息是可能的，从而识别参与机制驱动的区域。值得注意的是，机器生成的结构的特点是内部应力分布广泛，表明网络的集体参与。相反，在人类设计的结构中，应力突出了用于执行运动的少数支点。

算法的尺寸缩放特性。我们从实证角度研究该算法的尺寸缩放特性。为此，我们准备了六种不同尺寸的晶格，其键数从 $N_b = 175$ 到 $N_b = 885$ 不等。按照基于位移的方案，我们固定了相对于晶格总边长1%的输入位移，以及一个与输入方向垂直的目标方向，如图1所示。我们对每个晶格运行算法，总共进行 $N_a = 10^6$ 次接受的蒙特卡洛（MC）步骤，温度从0.06降低到0，然后在零温度下进行 $2 \cdot 10^6$ 次MC步骤。我们为每个晶格尺寸生成100个不同的结构，总计600个结构。对于每个结构，我们测量其总执行时间 $T$ 和效率 $\eta$ （见图3）。

缩放指数 $\alpha$ 是通过对 $\log T$ 和 $\log N_b$ 进行线性回归拟合得到的。我们估计 $\alpha = 2.97$ 。所得构型的效率随着 $N_b$ 的增加而提高（对于最小的晶格尺寸），然后在较大尺寸时达到平台期。我们注意到，蒙特卡洛过程也随 $N_b$ 缩放，因为算法的弛豫部分会贡献剩余的 $N_b^2$ 缩放。尽管一些计算细节或使用更快的CPU设备可能会降低报告的执行时间，但当前算法可研究的晶格尺寸存在上限。不过，我们的实现仍允许我们在大约一天内找到约700个键的高效结构，这处于机械超材料执行器的研究兴趣范围内。

效率对输入位移的非线性依赖。图2和补充图3中展示的分析依赖于初始三角形网格，但该算法也可使用随机网格执行，如“方法”中所述。对于相同的输入运动，周期性晶格和随机晶格获得的效率在图4中进行比较。我们观察到，在相对较小的输入位移下，基于三角形晶格的超材料优于基于随机三角形晶格的超材料，而对于较大的输入位移，两种晶格产生的效率相似（图4b，d）。

我们通过测量效率作为输入位移的函数来探究这一观察结果的起源。如图4c所示，效率是输入位移的非线性函数，在足够大的位移下趋于恒定值。在随机晶格情况下，效率在极小位移下迅速增加，而在三角形情况下，这种增加发生在高得多的位移处。非线性源于压缩主导行为向弯曲主导行为的转变，弯曲主导行为通过局部屈曲不稳定性产生大位移，在相对较小的力下作用，作为效率的乘法因子。

力输入/输出下的自动设计。在位移输入/输出条件下寻找最有效结构的策略可扩展到力输入的情况。为了说明这一点，我们考虑一个人类设计的参考构型，即参考文献36中提供的复制钳子的超材料机器——钳子是最常用的传统手持工具，用于夹持物体。钳子的机械运动由一对杠杆组成，它们在枢轴点连接，一侧是短颚，另一侧是长颚。由Ion等人提出的超材料机器由不同的单元格组成，这些单元格在支点附近较软，而在结构的其余部分更硬，以实现更大的灵活性³⁶。按压手柄（双反平行输入点应用）会使钳口对称闭合（沿y轴穿过支点），以握住物体（双反平行输出）。

我们通过在一组力节点上施加恒定的输入力，并测量与夹持钳口节点接触的一组量规弹簧上的力，来调整我们的方法以达到目标力基效率。然后效率按照式（3）中的定义来计算。

我们在人类设计的钳子上实施离散元法（DEM）模拟，采用合适的参数使机械响应用于线性弹性范围，选择相对较高的量规弹簧刚度（ $k_{\text{ext}} = 10.0$ ）和相对较小的输入力（ $F_{\text{ext}} = 0.01$ ）。目标是实现对小输入力的响应产生大的输出力。用这些参数达到的力基效率为22%。为了可视化实空间位移，我们还使用较软的量规弹簧（即 $k_{\text{ext}} = 0.01$ ，见图5a）进行模拟。接下来，我们运行我们的方法，自动生成一对钳子，其输入参数与人类设计的钳子相同（ $k_{\text{ext}} = 10.0$ 和 $F_{\text{ext}} = 0.01$ ），但在力传播方面的效率有所提高。实际上，我们要求力基效率是人类设计的两倍以上（ $\eta_f = 50\%$ ）。鉴于系统关于穿过支点的线的对称性，我们仅模拟系统的上半部分，然后镜像下半部分的结构。起始结构是包含人类解决方案的矩形全连接梁网络。在最小化过程中，枢轴节点的位置保持不变，y坐标对所有节点也保持相同。结果显示在图5b中，我们再次设置 $k_{\text{ext}} = 0.01$ 以可视化实空间位移。

我们还研究了 $F_{\text{ext}}$ 和 $k_{\text{ext}}$ 的选择如何影响最终的力基效率。图5c报告了在固定 $F_{\text{ext}} = 0.01$ （低于我们选择的值）和 $k_{\text{ext}} = 10$ （高于我们选择的值）的情况下，输出弹簧刚度 $\eta_f$ 的变化，适用于图5a和b中报告的人类设计和机器生成的结构。对于给定的输入力，我们观察到在大输出弹簧刚度值下达到饱和，两种解决方案的趋势相似。力基效率作为输入力的函数的变化估计，对于人类和机器生成的钳子机构，报告在图5c中。对于输入力在 $10^{-4} - 10^{-2}$ 范围内的情况，人类钳子的趋势是恒定的，而机器设计的趋势略有增加，其中我们观察到当 $F_{\text{ext}}$ 超过 $10^{-2}$ 时效率下降，这是由于屈曲效应。

机器学习预测效率。我们获得的大量超材料构型为将其视为待探究的数据集提供了进一步洞察的可能。在此背景下，人们自然会产生一系列与结构变化如何关联效率变化相关的问题。在此，我们评估是否可以通过构型的静态图像推断其效率，而无需执行离散元法（DEM）或有限元法（FEM）模拟。为回答这一问题，我们训练了一个卷积神经网络（CNN）来执行图像回归，即根据其布局的图像预测构型的效率。需注意，这与CNN更常用于图像分类的场景不同。利用大量构型（ $N \approx 10^6$ ），我们实现了 $R^2 = 0.966$ 的精度，详见图6a、c及“方法”部分。

通过机器学习生成新结构。我们将机器学习框架再推进一步，探究是否可以在从头生成新构型时，用CNN模型的预测替代弹簧 - 质量效率估计步骤。换句话说，我们使用上述相同的蒙特卡洛策略，但在每一步，不是通过DEM或FEM测量所提出构型的效率，而是用CNN模型估计它。有趣的是，该过程是成功的，我们能够生成新的、高效的构型，与DEM模型相比，速度提升约100倍。然而，我们注意到，我们仍然只需要训练CNN，因此仅通过结合DEM和CNN就能获得真正的速度提升。CNN生成的构型具有与Metropolis生成的构型相似的分布，见图6b。为确保获得的构型与训练CNN所用的构型不同，我们测量每个ML生成构型与训练期间使用的所有构型之间的最小距离，发现典型距离约为65个键。

用机器学习识别功能区域。CNN模型使用空间信息（结构的图像）和信息量（效率值）进行训练。这两种信息源在训练的CNN模型中有效地耦合在一起，可进一步用于识别功能相关的区域。我们方法的要点在于将结构的轻微扰动图像输入到CNN中，详见“方法”部分。

图7展示了一个构型的示例，突出了根据CNN（图7a）添加或（图7b）移除一个键时效率会发生显著变化的区域。为了确认这些确实是结构的功能相关区域，我们使用DEM模型进行类比计算，并将结果展示在图7c中。与图7a、b的视觉对比表明，CNN模型能够识别出结构的关键功能区域。我们系统地量化了这种效应，在五十次独立运行中，发现CNN预测的效率变化与DEM模型计算的变化之间存在高度相关性（见图7i）。该模型在不属于训练数据但本质上与训练数据非常相似的结构上进行测试，因此在一定程度上泛化了。

在图7e-h展示的第二个示例中，我们更进一步，探究同一CNN模型是否能识别出由更大尺寸晶格生成的结构的功能相关区域。具体来说，CNN模型使用172个键的晶格结构进行训练，但后来被用于识别694个键晶格结构生成的功能区域。图7e-h展示了694个键晶格上的一个结构：当向左结构添加键时，CNN正确地将左下角区域识别为效率降低的区域。虽然如预期的那样，预测和实际效率变化之间的相关性明显低于第一个示例，但存在系统性的偏差，倾向于正的 $R$ 值，在某些情况下高达0.4（图7j），这对于在CNN未接触过的晶格上进行预测来说是很显著的。

值得注意的是，图7c、d、g、h中的DEM计算是验证该方法的正确依据，但它们在计算成本上也更高，尤其是当从172个键的晶格切换到694个键的晶格时，见图3。能够使用在较小、简单晶格上训练的CNN模型识别功能相关区域，为众多可能性打开了大门，从降低蒙特卡洛（MC）算法中的拒绝率，到将人工和机器的工作结合起来生成混合结构。

讨论

在本文中，我们提出了用于具有广泛可能运动和效率的超材料致动器（MMAs）自动设计算法——其效率和运动能力超过了人类设计的解决方案。在其首次实现中，该算法利用离散单元模型来获得结构力学效率的近似值，随后驱动蒙特卡洛搜索遍历可能的结构。然后对优化后的结构进行有限元法（FEM）计算，以确认可通过3D打印实现的优化结构的效率。在此，我们专注于二维致动器，但使用相同策略也可扩展到完全三维模型。

有趣的是，将我们的算法与先前优化随机弹簧网络的方法³⁷⁻³⁹进行比较。与这些研究不同，我们专注于包含角度键的系统——因为这是重现真实梁机械特性的基础。事实上，我们的离散模型结果与有限元法（FEM）计算吻合良好。角度键的缺失会导致软模式（floppy modes）和刚度损失，而这些在真实梁的晶格中是不存在的。除了这一关键点外，先前的方法³⁷⁻³⁹还在算法的几个技术方面与我们不同。参考文献58中提出的方法实现了恒定温度下的Metropolis蒙特卡洛算法，其中单个移动涉及键的交换，因此总配位数得以保留。我们的模型没有软模式，因此无需施加这些限制，这使我们能够探索更广泛的相空间。此外，当目标是找到复杂优化问题的全局最小值时，我们使用模拟退火（其性能优于恒定温度Metropolis优化）。参考文献37中提出的算法采用复杂的最小化策略，能够识别“移除后效率提升最显著”的键。然而，当效率景观复杂时，该算法原则上可能陷入局部极小值（此时移除任何一个键都无法提高效率）。

此外，我们采用深度神经网络从结构预测致动器的效率。一旦经过适当训练，神经网络可用于创建新结构，而无需执行DEM或FEM模拟。使用机器学习辅助MMAs的自动设计，在算法速度方面开启了有趣的可能性——因为它有可能设计出无法通过DEM高效模拟的更大结构。为此，我们探索了深度神经网络在“为给定结构提供机械信息”方面的潜力。一旦训练完成，网络能够从结构的图像中识别出“修改后会使效率增加或减少”的区域。从无序结构中识别重要功能区域也是最近在流网络背景下研究的目标⁴³。

总之，我们的工作代表了建立基础致动器（EA）参考库的第一步。随后可通过互连多个EA获得更复杂的致动器，在应用和灵活性方面有无限可能性。例如，这些算法可用于设计机器人和机器人中的运动部件，特别是在小尺度（此时表面积与体积比非常大，从而导致摩擦和磨损占主导地位）的情况下。其优势范围从“现成EA的可获得性”到“为非标准应用提供定制解决方案的可能性”。这将为工程师和材料科学家提供参考。

方法

基于三角形晶格的超材料

我们考虑坐标为 $\boldsymbol{\mathrm{R}}_B$ 的三角形晶格构型，可将其视为本征构型（inherent configuration） $^{44}$ 。这种构型（见图1a）是力学稳定的，由长度为 $l = r_0$ 的 $N_b$ 根梁组成，它们连接在 $N$ 个节点上。第 $t$ 个节点的坐标 $r_t = (x_t, y_t)$ ，两个节点之间的距离为 $r_{ij} = |r_i - r_j|$ 。然后选择两个相距较远的节点组 $i$ 和 $j$ ，分别代表输入和输出区域，我们定义两个归一化向量，分别标识期望的输入方向 $t_{\text{in}}$ 和输出方向 $t_{\text{out}}$ 。

基于随机晶格的超材料

随机超材料构型是从圆盘随机构型开始获得的，这些圆盘被机械稳定地压缩到模拟盒中，这是成簇算法中的标准操作。之后，每对重叠圆盘的质心通过键连接。如参考文献37所述，这种生成无序结构的方式确保了网络具有稳定且易于理解的力学性能。

效率

超材料的响应通过其效率 $\eta$ 来监测：

\eta = \frac{t_{\text{out}} \cdot (r_j - r_i)}{\sum_{k=1}^{n_{\text{in}}} t_{\text{in}} \cdot (r_k - r_0)}, \tag{1}

其中 $(r_0$ 和 $r_j)$ 分别代表输入和输出节点的原始位置，点积在输入（ $i$ ）和输出（ $j$ ）节点的数量上取平均。可以设想效率的替代定义，以增强优化结构的特定要求下的期望响应。我们可以提供两种替代方案。

基于方向的效率：搜索超材料构型集中在将输出位移向期望方向最大化。为此，我们推广方程（1）中的点积，引入权重函数 $f(y) = (2\cos(y/2)^{n - 1} - 1)$ ， $n \geq 2$ ，其中 $y$ 是期望输出方向 $t_{\text{out}}$ 与测量方向之间的夹角。得到的效率为：
$\eta_d = \frac{|r_j - r_0| f(y)}{t_{\text{in}} \cdot (r_k - r_0)}. \tag{2}$

当 $n \gg 2$ 时，这会强制输出运动沿 $t_{\text{out}}$ 方向；而当 $n = 2$ 时， $f(y) = \cos(y)$ ，因此 $\eta_d = \eta$ 。
基于力的效率：在这种实现中，要求施加在输入节点上的力能有效地向输出节点传播，以指向目标方向。当致动器预计与其他机械部件集成（形成更大的机构）时，这一点尤其可取。在这种情况下，我们对输入节点施加恒定力，并通过充当测力计的量规弹簧，测量输出节点沿 $t_{\text{out}}$ 方向的力。这对应于添加能量项 $E_{\text{imp}} = F_{\text{ext}} t_{\text{out}} \cdot (r_i - r_0)$ 到输入节点，以及 $E_{\text{out}} = \frac{1}{2} k_{\text{ext}} (t_{\text{out}} \cdot (r_j - r_0))^2$ 到输出节点（在能量最小化期间），其中 $F_{\text{ext}}$ 是恒定输入力， $k_{\text{ext}}$ 是输出量规弹簧的刚度。基于力的效率可以简单地定义为：
$\eta_f = \frac{\sum_j k_{\text{ext}} |r_j - r_0| f(y)}{\sum_i F_{\text{ext}}}, \tag{3}$

其中分子的和是对输出节点求和，分母的和是对输入节点求和。

优化

一旦选择了合适的效率函数，我们通过最小化 $MC$ 函数 $\Delta = 1 - \eta$ 来最大化效率。优化协议是每次迭代都使用 $MC$ 函数结合优化算法：从一个现有构型（ $\Delta = \Delta^0$ ）开始，通过移除或读取随机选择的键来获得一个试验构型。在修剪过程中，输入节点和输出节点以及受约束防止运动的选定节点集（即冻结，见图1e）被丢弃。然后我们沿 $t_{\text{in}}$ 方向（或在 $\eta_f$ 的情况下对其施加外力）移动输入节点，执行FIRE优化 $^{45}$ ，并通过监测输出节点的位移（或在 $\eta_f$ 的情况下监测它们上的力）来评估相应的 $\Delta^{\text{trial}}$ 。如果 $\Delta^{\text{trial}} < \Delta^0$ ，则移除键被接受，否则它以概率 $P = \exp(-(\Delta^{\text{trial}} - \Delta^0)/T)$ 被接受，其中 $T$ 是作为 $MC$ 动力学中温度的参数。

在每次最小化开始时，为了探索复杂的效率景观 $^{46}$ ，我们进行100次（接受的）退火步骤，使 $T$ 从 $T = 0.06$ （这是持续获得 $P \simeq 1$ 的阈值）线性下降到 $T = 0$ ，随后让算法在后一温度下演化。整个过程使用不同的随机数种子重复多次（见补充图1）。

值得注意的是，在一个有 $N_b$ 根梁的系统中，构型空间的数量级为 $\sim 2^{N_b}$ 。在图1所示的示例中， $N_b = 203$ ，产生了极大的构型数量： $2^{203} \approx 10^{61}$ 。显然，这种随 $N_b$ 的指数级增长严重限制了探索所需构型空间的可能性，因此需要快速方法来预测试验结构的效率。为了最大限度地提高这种探索效率，我们采用了三种不同方法在不同近似层次上协同作用，如下所述。

离散单元模型

为了快速可靠地估计给定结构的效率，我们建立了晶格的简化离散单元模型（DEM），其中总能量可表示为：

E = \sum_{i} \sum_{j>i} \phi_1(r_{ij}) + \sum_{i} \sum_{j>i} \sum_{k} \phi_2(r_{ij}, r_{ik}, \theta_{ijk}), \tag{4}

其中成对项是具有平衡长度 $r_0$ 的弹簧势能：

\phi_1(r_{ij}) = k(r_{ij} - r_0)^2, \tag{5}

而三体项在同节点连接的第一近邻梁之间引入角弹簧：

\phi_2(r_{ij}, r_{ik}, \theta_{ijk}) = \lambda [\theta_{ijk} - \theta_{ijk}^0]^2, \tag{6}

其中 $\theta_{ijk}$ 是由梁 $\vec{ji}$ 和 $\vec{ik}$ 形成的角， $\theta_{ijk}^0$ 是三角形晶格中的初始角值。 $\phi_1$ 和 $\phi_2$ 仅在每一步对第一近邻节点动态重新计算（见补充图2）。如果未另行说明，我们采用以下无量纲参数： $k = 5$ ， $\lambda = 0.1$ ， $r_0 = 1$ 。

在冻结节点（通常在结构的底部）和错位输入节点存在的情况下最小化 $E$ ，使我们能够快速预测试验结构的响应，同时在速度和可靠性之间取得良好的平衡。这将作为我们寻找最高效率结构的参考方法（见图1b）。

有限元法

通过有限元法（FEM）可以获得受外力作用的结构力学响应的真实模拟。该方法基于对象的网格表示，以便将连续边界值问题转化为代数方程组系统。该方法在弹性范围内非常准确（计算时间在几秒量级，适用于我们的参考模型），但随着系统尺寸的增加而迅速增加，因此不适合在3D打印前通过MC搜索高效结构的试用阶段。相反，FEM已被用于在通过更简单快速方法生成的选定结构的最终验证之前。

模拟结构的三维模型是通过挤压每个键生成的，对于FEM计算，我们使用COMSOL Multiphysics和COMSOL with MATLAB通过其结构力学模块 $^{47}$ 。所有研究都假设线性弹性材料，其杨氏模量和泊松比根据感兴趣的大块样品的实验估计。结果使用Euler-Bernoulli梁单元和内置稳态研究计算（准静态求解器）

机器学习。深度学习模型：我们采用Resnet50架构作为替代，在Python Keras库中对其进行修改，使其用于回归而非分类，具体做法是将顶层替换为一个单单元全连接层。

在数据准备阶段，我们生成了从一万次DEM（离散元法）模拟中获得的结构的192×128像素PNG图像，方向为期望的 $\boldsymbol{t}_{\text{out}}$ 方向，包括被接受和被拒绝的构型。CNN模型不仅需要能够识别构型在期望方向上是否存在位移，还需要能够在相反方向上移动，且其效率需理论上为负（以便在蒙特卡洛搜索期间将这些构型移走）。由于这些构型出现的频率非常低，我们通过运行另外一万次针对相反方向—— $-\boldsymbol{t}_{\text{out}}$ 方向的DEM模拟来增加其采样率。总共，我们使用Python绘图库matplotlib生成1,163,733张图像。

模型训练：我们使用标准的Adam优化器和均方误差损失函数。75%的运行结果用作训练数据，剩余的25%用作验证。我们还实现了一种简单的数据重采样策略，使效率的分布在整个训练数据中近似均匀。这缓解了这样一个事实：否则，模型对非常高效率构型的预测会不太准确，因为这些构型在我们的数据集中很少见。我们以32为批次大小，训练修改后的Resnet50 CNN共10个epoch。

计算环境和模型性能：我们使用配备Tesla K20c GPU的16核计算机进行图3、7中的计算，以及配备Tesla V100-SXM2的Tesla K20c GPU进行图6中的计算。我们使用实际效率和预测效率之间的线性回归 $R^2$ 值来衡量模型的性能。

人工设计的结构和3D打印。为了将自动设计的MMAs（机械超材料架构）与常规解决方案进行比较，人工设计的MMAs要么由我们在之前未进行的模拟中创建，要么来自文献，以避免对设计者关于导致高效率的机制产生任何偏见。

然后通过3D打印生产样品。在这种方法中，最终结构是分层制造的，每层厚度均为0.08毫米，通过连续铺设熔融塑料的薄层。喷嘴尺寸为 $n_b = 0.4$ 毫米，使用的材料是Ninjaflex，一种具有超弹性的配方热塑性聚氨酯（TPU）材料。

3D打印结构的效率是通过手动移动执行器输入节点5毫米，用卡尺测量位移，同时固定底部零件夹具来测量的。初始和变形后的输入和输出结构的图像已被捕获，输入和输出节点的位移已通过图像处理进行测量。通过对测量到的输入/输出节点位移应用公式(2)，计算了效率 $\eta_d$ ，并在图2和补充图1中报告。

功能热图的构建。我们通过将结构的扰动图像输入CNN模型来构建功能热图。功能热图显示当结构被局部修改时，总效率的预测增加或减少。对于任何给定的结构图像，我们构建两个功能热图：一个对应于添加键的影响，另一个对应于移除键的影响。然而，由于潜在训练数据的晶格和所考虑结构的晶格可能不一致，我们遵循一种与晶格无关的程序，如下：

在图像上滑动一个 $10×10$ 的扰动方块。扰动方块根据我们正在构建的热图，要么是黑色（添加）要么是白色（移除）。我们每次将方块滑动10像素，以避免重叠的扰动，但其他更精细的实现也是可能的。在200×200像素的图像中， $10×10$ 的方块可以放置在400个不同的位置，从而生成400张扰动图像。
将400张扰动图像输入CNN模型，并预测这些图像的效率。
减去原始图像的预测效率，得到400张扰动图像的效率预测变化。
构建一个 $20×20$ 的效率变化数组，恢复扰动方块的位置。
使用三次插值放大图像，得到最终的 $200×200$ 像素的功能热图。

CNN和功能热图之间的相关性。我们使用键级的相关系数作为CNN预测的效率变化与DEM模型预测的效率变化（以功能热图形式）之间相关性的简单度量。DEM模型作为基准真值。为了整合热图中与每个键对应的线条的效率变化值，我们定义相关性为Python中scipy.stats.linregress函数实现的线性回归模型的 $R$ 值。

图1 自动结构生成的示意图。a 初始三角形晶格构型 $\boldsymbol{R}_{is}$ 。b 通过DEM（离散元法）优化得到的结构及其在输入位移下的模拟响应。c 相同的结构运动由FEM（有限元法）模拟。颜色图显示冯·米塞斯应力。

图2 自动设计实现高效的反平行运动。（a）人工设计和（b）机器设计结构的有限元模拟比较，以及（c、d）相应的3D打印实物。颜色表示冯·米塞斯应力，色标与图1c相同。红色（绿色）箭头表示期望输入 $\boldsymbol{t}_{\text{inp}}$ （输出 $\boldsymbol{t}_{\text{out}}$ ）的方向。每个面板中都报告了最终的效率。关于运动的可视化请参见补充影片2。

图3 算法的缩放特性。a 六组规模逐渐增大的规则晶格的总执行时间（单位：小时）与键数的关系。大的空心圆显示中位执行时间。小的灰色点显示100次独立模拟的原始数据。红线是对数空间中的线性拟合。该面板显示，算法的缩放大致与键数的立方成正比。b 基于位移的效率随键数变化的箱线图。底层原始数据也以黑点显示。箱线图中的误差线表示最小和最大四分位数，箱子的范围为第一和第三四分位数。该面板显示，效率最初随 $N_b$ 增加而上升，然后趋于饱和，但在固定晶格尺寸下存在很大的异质性。

图4 晶格几何形状与屈曲效应。a 具有正交输入 - 输出运动的致动器示例，由三角形（顶部）和非晶态（底部）晶格得到。b 与（a）相同，但输入节点的位移更大。红色（绿色）箭头表示期望输入 $\boldsymbol{t}_{\text{inp}}$ （输出 $\boldsymbol{t}_{\text{out}}$ ）的方向。c （a）、（b）中两个相同致动器的效率随输入位移的变化关系。误差线为标准偏差。d 使用三角形或非晶态晶格生成的100个致动器的效率分布。箭头表示两种情况下的平均效率。

图5 机器设计的钳子比人工设计的钳子具有更高的基于力的效率。a 人工设计的钳子及相应的（b）通过离散元法（DEM）模拟得到的自动设计的钳子。机器设计的构型与人工构型在 $F_{\text{ext}} = 0.01$ 和 $k_{\text{ext}} = 10.0$ 时进行比较。为了能够可视化位移，我们展示了两种结构在相同输入力 $F_{\text{ext}} = 0.01$ 但更低输出弹簧刚度 $k_{\text{ext}} = 0.01$ 下的最小化结果。元素根据其相对于灰色下划线的静止构型的位移大小进行着色。输入（输出）节点在冻结枢轴点（中心标记）的右侧（左侧）用更大的圆圈标记。下方输入和输出节点的放大图突出了位移场。红色（绿色）箭头表示期望输入 $\boldsymbol{t}_{\text{inp}}$ （输出 $\boldsymbol{t}_{\text{out}}$ ）的方向。c 基于力的效率随输出弹簧刚度 $k_{\text{ext}}$ 的变化，其中 $F_{\text{ext}} = 0.01$ 为固定值，d 随期望输入力 $F_{\text{ext}}$ 的变化，其中 $k_{\text{ext}} = 10.0$ 为固定值。

图6 机器学习可用于预测高效并设计结构。a 测试集在10个训练周期内的决定系数 $R^2$ 。b 用离散元法（DEM）模型（灰色）和卷积神经网络（CNN）模型（蓝色）生成的构型的效率分布。注意，最终优化构型的效率始终用DEM测量。c DEM模型得到的效率 $\eta_{\text{DEM}}$ 和从CNN模型得到的 $\eta_{\text{CNN}}$ 的密度图，显示CNN模型能够很好地预测DEM模型的效率值。实线表示线性拟合，其决定系数为 $R^2 = 0.966$ 。虚线仅作为视觉引导。

图7 功能热图识别结构中与力学相关的区域。热图显示了在172键晶格的一个示例结构中，由CNN模型（a、b）或DEM模型（c、d）计算的、与键添加（a、c）或键移除（b、d）相关的效率变化的空间分布。蓝色阴影对应导致效率降低的区域，而红色区域导致效率提高。i CNN和DEM的 $\Delta\eta$ 值的相关系数 $R$ 的分布，详见方法部分。e-j 对694键晶格的结构进行相同分析，但仍使用用172键结构训练的CNN模型。该图表明，较小晶格的信息可以转移到较大的晶格上，并且仍然可用于识别结构的功能区域。

补充图1：效率在不同路径上的演化。正交运动情况下蒙特卡洛（MC）动力学过程中瞬时 $\eta$ 的绘图。每条曲线由不同的初始随机种子获得。阴影区域突出显示前100个MC步骤中的退火阶段（见方法）。仅报告对应于被接受的MC步骤的效率。箭头和虚线标记人工设计结构的效率。

补充图2：角弹簧重新邻接。当移除一根梁（虚线）时，连接到该梁的两个角弹簧被移除，并且在新的相邻梁之间形成一个新的弹簧。

图10

补充图3：自动设计实现高效的直交运动。针对有限元模拟的（a）人工设计和（b）机器设计结构以及（c,d）对应的3D打印实物进行比较。颜色梯度使用与图1c相同的色标报告冯·米塞斯应力。每个面板中都报告了最终的效率。有关运动可视化的内容请参见补充视频3。