用于轻量化结构设计的可微分力密度方法

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Differentiable force density method for the design of lightweight structures

　　拉斐尔·帕斯特拉纳 | 丹尼斯·奥克塔伊 | 凯-乌韦·布莱廷格 | 瑞安·P·亚当斯 | 西格丽德·阿德里安森斯
普林斯顿大学建筑学院，美国普林斯顿

**摘要**
轻质屋顶、桥梁和塔楼凭借其机械效率高的形状，能够以细长的截面跨越长距离。绳索几何结构（即外部载荷主要通过内部轴向力传递的空间配置）可以使这些结构具备这种效率。然而，由于非线性力学行为和高维搜索空间，将绳索几何结构转化为满足建筑、制造和结构要求的可行设计仍然具有挑战性。本文提出了可微力密度方法（DFDM），这是一种用于自动设计轻质结构高效且可行形状的方法，这些结构被建模为销接杆系统。DFDM利用基于梯度的优化来计算受设计约束的绳索几何形状，同时考虑了大位移效应和形状依赖的载荷（如自重和压力）。通过扩展自动微分与解析伴随函数，DFDM解决高维逆问题的速度比无梯度算法快四个数量级，比通用自动微分快一个数量级。其完全可微性使得载荷查找和支持位置查找成为可能。这些丰富的参数化允许联合优化载荷和支持位置与力密度，从而发现比单独使用力密度优化更优的解决方案。这项工作将结构模拟与可微编程结合起来，为加快开发降低建筑环境中材料消耗的米级系统开辟了途径。

**1. 引言**
像膜结构屋顶、悬索桥和网格壳这样的轻质结构在建筑尺度上是机械效率高的系统。它们的高刚度重量比使它们能够在数十或数百米的跨度上承受大载荷，且截面细长，因此所需材料体积小于依赖框架或梁柱结构的系统[1][2]。它们的机械效率常与自然界中的系统相媲美，例如鸡蛋壳的厚度仅为其半径的2%左右[3][4]。几何形状是这些结构卓越性能的核心。轻质结构主要是形状活跃的[5]，它们通过平面内的轴向力抵抗外部载荷，同时消散弯曲和扭矩，从而减少结构的应变能[6][7]。在连续系统（如壳体）中，这种力学响应被称为膜作用；而在离散情况下（如杆件和缆绳的组合，或更一般的杆件），则称为桁架作用。

**1.1. 寻找高效形状**
高效形状是设计轻质结构的关键，但无法从一开始就确定。这种形状取决于结构的行为、施加的载荷和支持条件，因此必须通过计算来找到[8]。
一种常见的计算高效形状的方法是使用形式寻找方法，该方法将所需的轴向主导行为施加在销接杆系统上，这是一种离散的结构表示。杆系统可以直接模拟缆网和网格壳，并通过将平面内的膜力浓缩成类似杆-索模型的力结果网络来近似膜结构和壳体[9][10]。在一般情况下，无论是支撑还是无支撑的销接杆系统中，强制采用桁架作用并规定力的方向（拉伸或压缩）都会引入几何非线性。因此，杆系统必须经历大的位移才能达到平衡状态。由此产生的配置是绳索几何形状：每个杆都承受一个指定方向的轴向力，且系统在给定边界条件下与施加的载荷处于平衡状态，无论其拓扑结构如何。具有绳索几何形状的结构可以降低其应变能并实现机械效率。
存在多种形式寻找方法，其历史可追溯到Varignon[11]、Rankine[12]和Maxwell[13]关于图解静力学的开创性工作，以及Pucher的膜理论公式和Airy应力函数等基本连续介质方法[14][15][16][17]。一种已建立的分类法[18]将销接杆系统的数值形式寻找方法分为三类：基于刚度的[19][20]、动态的[21][22][23]和几何方法[24][25][26][27]。几何方法的优势在于它们避免了显式的材料模型，为设计提供了更紧凑的参数化，同时不影响模拟精度。在这种分类法中，力密度方法（DFDM）是所有其他几何方法的基础[18]。
力密度方法（FDM）[28][29]是一种用于计算销接杆系统绳索几何形状的形式寻找方法。自20世纪70年代初慕尼黑奥林匹克体育场建设以来[30]，FDM已被扩展用于处理压缩结构（如混凝土壳和砌体拱顶[31][32][33][34]）以及拉伸-压缩结构（如张力网格和悬索桥[35][36][37][38]）。尽管FDM最初是为模拟拉伸缆网而开发的，但其后续扩展处理压缩力状态的过程直观地遵循了胡克关于悬挂链和刚性拱的对偶性形式[39][40][41]。对FDM的进一步研究揭示了力密度与Piola-Kirchhoff力之间的等价性[42]，这为将该方法扩展到连续结构（如膜和壳体）提供了重要基础[25][35][43]。FDM在形式寻找领域之所以占据重要地位，是因为其在模拟方面的数值可行性。它使用力密度（即内部力与构件长度的比率）而非力值作为模型输入，重新表述了原始的非线性平衡问题[28][29]。这种重新参数化使平衡状态的计算线性化，从而可以直接访问丰富的绳索几何形状空间。

**1.2. 结构设计是一个逆问题**
像FDM这样的形式寻找方法支持轻质结构设计，但仅限于一定程度。FDM擅长通过解决正向问题（将参数映射到几何形状θ→x，见图1a）来计算绳索几何形状，并通过简单改变θ来探索众多高效形状（见图1b）。然而，结构设计本质上是一个逆问题[44][45][46]，其中绳索几何形状必须满足设计规范规定的其他约束。除了效率之外，结构还应使用特定材料和制造工艺建造[47]，符合施工预算限制[7]，并响应建筑和场地要求[38]。因此，设计过程中需要寻找满足这些约束的参数θ，从而产生满足机械效率要求的几何形状x^。然而，这种关系与FDM和模拟的整体正向性质相反（见图1a）。

**下载：** 下载高分辨率图像（558KB）
**下载：** 下载全尺寸图像
**图1.** 像FDM这样的结构模拟方法根据给定参数θ∈θ计算绳索几何形状x∈X，将Rp中的模拟输入连接到Rn中的机械行为。FDM建立了一个正向映射，为θ的不同变化生成多种绳索几何形状，如图1b所示。所有这些形状都是绳索几何形状，包括xC。然而，轻质结构设计的方向与模拟相反，需要寻找产生满足设计规范的绳索几何形状x的参数θ。在这个例子中，该规范是匹配一个编码建筑视觉意图的目标形状x^，但x^不在绳索几何形状的空间X中。因此，设计构成了一个逆问题，即从Rn中的行为到Rp中的模拟输入。

**1.3. 通过优化寻找高效且可行的形状**
优化为自动搜索FDM诱导空间中的特定绳索几何形状提供了一种原理性的方法。通过将设计规范转换为标量值损失函数，优化算法可以迭代调整模拟参数θ以最小化损失值，并收敛到满足设计规范的绳索几何形状[46][52]。因此，优化建立了从绳索几何形状到参数x→θ的逆映射，但只是近似的，其中x≈x^。损失值越低，发现的几何形状与规范匹配得越好（见图1）。
无梯度方法（如遗传算法和进化策略）在结构优化中经常被使用，因为它们具有灵活性[53][54][55][56][57][58][59][60]。它们只需要评估损失函数以进行适应性评估，因此能够灵活应对一般的设计规范和损失函数公式[61]。然而，无梯度方法在模拟成本和问题维度方面的扩展性较差[62]，因此最适合小规模问题。基于梯度的优化在快速、高维搜索方面具有计算优势。这种方法利用损失函数的梯度（即损失函数对θ的导数方向）来指导搜索，从而在更少的计算步骤中找到最优解[52][63]。许多先例证明了基于梯度的优化在解决高层建筑[64][65][66]、张力网格[35][67]、悬索桥[68]和壳体[44][69]等复杂设计任务方面的强大能力。然而，由于其在适应设计变更方面的刚性，这种方法在一般结构设计中的应用仍然不足。主要障碍在于推导损失函数梯度的解析表达式。这个专门的手动过程劳动强度高，每当设计规范更新时都需要重复。高昂的人工成本阻碍了其采用，促使人们寻求自动导数计算方法，使基于梯度的设计优化不仅更具适应性，而且更适合更广泛的研究人员群体。

**1.4. 自动计算精确导数**
自动微分（AD）是一系列算法技术，可以在不进行手动推导的情况下计算计算机程序中表达的函数的导数[70][71]。它是可微编程的基石[72]。AD背后的主要思想是将函数表示为组装成计算图的基本操作的组合[73][74]。导数计算通过遍历这个图并递归应用微分的链式法则来完成。当从输入到函数的输出进行遍历时，称为正向模式AD；当反向进行时，称为反向模式AD或反向传播[75][76]。
与其他计算微分技术（如数值微分或符号微分）相比，AD在数值上更便宜且更具通用性。数值微分只能近似梯度，并且每个需要微分的参数至少需要调用一次函数，使得在大规模问题中每一步优化成本都很高[72][77]。符号微分虽然提供了具体的数学梯度表达式，但仅限于可微分解析函数，并且存在表达式膨胀的问题[34][78]。相比之下，AD能够以机器精度提供精确的梯度，计算工作量仅为被微分函数的执行成本的一小部分[70][72][79]。由于其可扩展性、精确性和灵活性，自动微分（AD）成为了现代深度学习的支柱，在深度学习中，反向传播（backpropagation）提供了训练各种神经网络所需的梯度，以解决复杂的视觉识别或自然语言处理任务[78]、[80]、[81]。除了神经网络训练之外，可微分模拟（differentiable simulation）是另一个领域，在这个领域中，AD作为一种工具，用于解决受物理定律支配的复杂优化问题。1.5. 用于快速结构设计的可微分模拟可微分模拟将AD的计算优势从数据驱动模型扩展到了基于物理的模型，如有限差分法（FDM）。可微分模拟方法在AD框架中实现基于物理的模型，从而可以自动获得模型输出对其输入的导数，以实现基于梯度的设计优化，适用于选定的损失函数和问题参数化[82]。需要注意的是，虽然这些模型在数学上是可微分的，但它们的数值实现必须在计算机程序中精心设计，以确保与AD兼容。可微分模拟方法正在推动工程和物理科学各个领域的设计优化，例如机器人技术[83]、[84]、[85]、[86]、流体和分子动力学[85]、[87]、[88]以及超材料[89]、[90]。结构设计也不例外。可微分结构模拟方法的先例主要集中在平面域的拓扑优化[91]、[92]、[93]以及受应力和屈曲约束的3D固体[94]上。相关工作还提出了用于减少柔顺性的可微分结构分析框架[95]、[96]和减少嵌入式碳含量的方法[97]。然而，这些方法仅限于小位移、线性弹性模型，这些模型不支持大位移，并且需要拓扑三角剖分（例如支撑结构）以及明确的截面和材料参数来模拟杆系[98]。此外，线性弹性模型无法预先施加目标拉伸或压缩力状态，而这正是设计轻质结构（如缆网（拉伸状态）或砌体壳（压缩状态）所必需的。对于杆系的逆向形态寻找，存在推力网络分析[99]、[100]和组合平衡建模框架[77]的可微分变体。然而，它们专注于特定结构，如壳和层状桁架，并且需要对杆系的拓扑进行仔细的预处理才能进行模拟，这最终限制了其通用性。为了提高泛化能力，先前的工作已经展示了将FDM与AD相结合的方法[34]、[101]、[102]。但这些努力主要集中在软件实现上，而不是FDM和AD背后的数学原理。通过利用关于FDM分析导数的丰富文献[29]、[42]，并研究如何将其融入AD的反向传播过程中，可以实现更紧密的耦合，从而有利于快速的结构优化。此外，这些FDM的先例主要集中在优化输入力密度上。这种限制缩小了设计空间，相比之下，当可能时联合优化载荷和支持条件可以探索更广泛的设计空间，从而在严格的设计规范下发现最佳结构形状。无论是专注于结构分析还是形态寻找，一个更普遍的问题是：当前的可微分结构模拟方法只考虑恒定的外部载荷[34]、[95]、[96]、[97]、[101]、[102]。虽然这种简化便于计算，因为它允许在一步模拟中预测结构行为，但它排除了考虑那些基于结构在模拟过程中不断变化的几何形状的载荷，而对此有分析方法[42]、[103]、[104]。这些依赖于形状的载荷包括自重、风或水的作用力（即流固耦合），或者用于加固气动系统（例如垫子或气球）所需的内部压力。开发能够考虑大位移效应和形状依赖载荷的可微分形态寻找方法将填补这一空白，并更好地支持现实世界中的轻质结构设计。1.6. 贡献在本文中，我们提出了可微分力密度方法（DFDM），这是一种利用自动和分析微分来计算满足建筑、制造和结构约束的缆状几何形状的可微分形态寻找方法。因此，DFDM简化了作为受形状依赖载荷影响的铰接杆系模型的形状设计。本文为此目的做出了三项具体贡献，扩展了经典FDM[28]、[29]的能力，并在FDM和AD交叉领域的前期研究[34]、[101]、[102]的基础上进行了发展：• 用于结构设计优化的可微分机械模型。我们在大位移范围内对铰接杆系的机械行为进行了建模。与现有的可微分结构模拟方法不同，我们的可微分模型通过基于更新参考策略的迭代平衡求解器支持形状依赖载荷，如自重和压力（第2.1-2.4节）。重要的是，整个优化过程中都保持了力平衡，而不仅仅是在收敛时，从而增强了结构设计的鲁棒性（第2.5节）。• 利用分析伴随矩阵和自动微分加速梯度计算。AD简化了导数计算，但由于重复操作而产生了计算开销。我们在DFDM中识别出这些瓶颈，并通过隐式微分导出了分析向量-雅可比矩阵乘积来解决它们。通过将这些专门的伴随矩阵注入通用AD中，我们利用所提出的机械模型的稀疏性和固定点结构来更快地计算梯度，同时保持了AD的灵活性（第2.6-2.8节）。• 为受限任务中的最佳形状发现扩展参数化。我们引入了由DFDM的全微分性启用的新的载荷寻找和支持寻找参数化。我们证明了它们在具有多个约束的复杂设计任务中的优势，例如通过外部后张力修复砌体壳和设计双曲面网格壳和平面板。这些任务，加上形状匹配和长度控制问题，突显了DFDM为所有模拟参数提供精确梯度以最小化通用损失函数的能力，从而发现高效且可行的结构形式（第3节）。我们通过复现过去研究[31]、[105]、[106]、[107]、[108]的分析和数值实验来验证DFDM的有效性和正确性。这些贡献共同建立了一个统一且可扩展的框架，用于基于梯度的轻质结构形状优化，解决了当前结构模拟方法的局限性，并更好地调整了形态寻找，特别是在概念阶段，快速迭代至关重要，设计决策对结构的性能和成本有最大影响[51]、[109]。1.7. 符号说明本文使用小写斜体a表示标量和函数，小写粗体字母a∈Rn表示列向量，大写粗体字母A∈Rn×m表示矩阵。集合和集合由大写斜体A表示。导数用?表示，无论是总导数还是偏导数。当一个变量依赖于一个输入时，导数应从总体的意义上解释。vec(A)运算符表示矩阵A的列向量化。我们通过将矩阵小写来表示这种重塑操作的结果，即a=vec(A)∈Rmn。相反，我们通过按顺序填充矩阵列将向量a重塑为矩阵A，并用大写向量符号表示输出矩阵，即A=matrix(a)∈Rn×m。当向量带有小写斜体下标时，如ai，它表示A的第i行。我们使用罗马数字下标来表示变量类别，如amax或aA。下标(??·??)0用于表示变量的初始值。2. 可微分力密度方法在本节中，我们介绍了DFDM和本工作的主要方法论贡献。首先，我们介绍了形态寻找问题（第2.1节），然后从能量角度解释了力密度方法的机械模型，以证明其与结构力学的一致性（第2.2节）。然后，我们介绍了将FDM扩展到一般铰接杆系的方法（第2.3节和2.4节），并将其作为可微分模型来计算符合约束的缆状几何形状（第2.5节和2.6节）。随后开发了线性系统求解器的自定义伴随矩阵和更新参考策略求解器的公式（第2.7节和2.8节）。附录提供了本工作的数值验证和简洁的DFDM实现。2.1. 形态寻找问题考虑一个由四个杆组成的铰接系统，其位置为x1∈R3，支撑位置为sj∈R3；内部力为fj∈R（图2）。这些杆在x1处连接到一个无旋转的关节，并且在s1到s4的位置固定在地面上。在x1处有一个向下的载荷p1∈R3。对于给定的fj，形态寻找问题的解决方案是找到x1的位置，以消除残余力r1∈R3：(1)r1(x1)=∑jfjlj(x1?sj)?p1=0如果r1的范数为零，则结构处于平衡状态。注意，杆的长度lj=∥x1?sj∥∈R是未知x1的非线性函数，因为它涉及到计算杆端点之间坐标差的范数。可以通过规定内部力与杆长度之间的比率作为输入来更易地解决方程(1)：(2)r1(x1)=∑jqj(x1?sj)?p1=0其中比率qj∈R是第j根杆的力密度，定义为：(3)qj=fjlj我们可以使用牛顿法解决方程(2)。相应的牛顿步骤Δx为：(4)?r1?x1Δx=?r1由于r1对x1的线性，我们可以在一步中找到方程(2)的根：(5)x1★=x1?K?1r1其中x1表示图2a中初始配置中自由顶点的位置，K∈R3×3是雅可比矩阵?r1/?x1，即系统的刚度矩阵：(6)K=[∑jqj000∑jqj000∑jqj]如果展开方程(5)，我们得到了一个解析表达式来解决问题：(7)x1★=∑jqjsj+p1∑jqj这个表达式提供了一个关键见解：平衡解x1★，如图2b所示，与x1无关。这个特性避免了由于x1的初始猜测不佳而在计算残差根时可能出现的数值发散，只要K是可逆的（例如，其对角线元素都不为零，∑j?qj?≠?0）。使用力密度作为输入来解决形态寻找问题，并将方程(2)扩展到一般铰接杆系，是FDM的基础，我们将在第2.3节中详细介绍。下载：下载高分辨率图像（207KB）下载：下载全尺寸图像图2. 解决形态寻找问题的目标是计算铰接杆系的自由3D坐标，使残余力消失。所得到的几何形状达到平衡状态。对于图中受到载荷向量p1影响的4杆系统，初始配置不是平衡的，因为位置x1处的残余力向量r1的范数不为零。图形静力学[110]、[111]提供了一种视觉方法来评估系统的非平衡状态，因为如果没有r1，其相应的力图将是开放的。只有将x1平移到x1★后，杆系才能达到平衡，使其成为缆状结构，如图下方的封闭力图所示。2.2. 机械模型在介绍通用FDM之前，我们从能量的角度对其进行解释，补充了第2.1节中介绍的传统力平衡方法。这种观点对于解释基于结构力学的几何非线性分析的FDM机械模型是必要的，在这种分析中，缆网或一般情况下的铰接杆系受到预应力作用和大位移的影响。我们的阐述整合了之前研究[42]、[43]、[112]中的思想。我们考虑图2a中所示的其中一个杆，并对其施加轴向预应力f，同时忽略弹性变形。在其初始配置中，杆的内能πint是第二Piola-Kirchhoff力fPK2和Green-Lagrange应变εGL的乘积，积分 over 杆的初始长度l0：(8)πint=∫0l0fPK2εGLdl=fPK22l0(l2?l02)与一阶应变不同，几何非线性应变εGL考虑了大位移，并确保了与杆的变形长度l的机械兼容性：(9)εGL=l2?l022l02根据最小势能原理，可以通过最小化系统的总能量来找到图2所示杆系统的平衡状态：(10)Π=Πint+Πext→min其中Πint是各个杆的能量之和，Πext是由施加的载荷p1产生的外部能量。在离散化杆之后，通过将能量对虚拟位移δx1的梯度（即第一变分）设为零来获得方程(10)的静止点：(11)?x1Π=∑j?x1πint+?x1Πext=∑j(fPK2ll0)j?lj?x1δx1?p1δx1=0位置x1编码了自由顶点的坐标，(?·?)j表示与第j根杆相关的量。长度lj关于x1的导数是变形后，杆端位置标准化差异与杆长度的比值：
(12) ?lj/?x1 = (x1 ? sj) / lj
其中sj表示第j根杆的固定端位置。将方程(12)代入方程(11)并假设能量的第一变分对所有兼容的虚拟位移δx1都成立，我们得到平衡条件：
(13) ?x1Π = ∑j(fPK2l0)j(x1 ? sj) ? p1 = 0
注意，在变形配置中的长度lj在代数上相互抵消，从而简化了上述表达式。方程(13)等同于方程(1)中的非线性形式寻找问题，并且在力向量共面时（例如，当结构中的三个顶点共面时）会出现奇异性[42]。如果我们通过定义第二Piola-Kirchhoff力与未变形长度之间的恒定比率来重新参数化方程，即
(14) q = fPK2l0
我们可以线性化能量泛函的变分，克服原始公式中出现的奇异性，并再次得到FDM所基于的力密度q。经过这种重新参数化后，方程(13)的解会得到一个悬链线几何形状，如图2b所示。
FDM的结构力学处理解释了为什么即使需要较大的位移，销接杆系统也能达到平衡状态，同时保持预应力力和非线性应变的机械一致性。从能量最小化的角度来看，势能的梯度?x1Π（方程(11)）对应于图2中所示的杆系统中的残余力r1（方程(2)），这是结构平衡的核心度量。此外，该能量的Hessian矩阵?x12Π给出了结构的刚度K，它等同于残余力r1的雅可比矩阵。这两点观察揭示了FDM与其他形式寻找方法以及更广泛的结构力学领域之间的密切关系，正如之前的研究[42]、[112]所指出的。

2.3. 力密度方法
在理解了FDM所基于的力学模型的基础上，我们现在将其扩展到三维结构。为此，我们简要概述了如何使用FDM为给定参数θ∈Rnp（定义在图G上）计算悬链线几何形状x∈R3?。图3总结了输入和输出。本工作中FDM的表述采用了向量和矩阵运算，即线性代数，遵循Schek[29]的原始处理方法。

下载：下载高分辨率图像（282KB）
下载：下载全尺寸图像

图3. 力密度方法（FDM）计算由θ参数化的销接杆系统的悬链线几何形状x，其拓扑结构由图G编码。参数包括力密度q、施加的载荷p和支持位置s。阴影区域A和B进一步详细说明了与G的边和顶点相关的参数。在A中，我们展示了FDM的输入：边3的力密度q3、支撑顶点2的位置s2以及顶点1处的载荷p1。在B中，特写显示了FDM的输出：平衡位置x1、边3的内力f3、边4的长度l4以及支撑顶点2的反作用力t2。

2.3.1. 数据表示
我们将模型化为销接杆系统的结构的拓扑结构表示为一个具有nv个顶点和ne条边的图G。顶点编码无旋转的关节，边则捕捉这些关节处连接的杆。每个关节只有一个顶点vi，每根杆只有一条边ej。在模拟中，G的拓扑结构嵌入到一个连通性矩阵C?∈{?1,0,1}ne×nv中，其中(cij)项定义如下：
(15) cij = {1 if i是边j的起始顶点；0 otherwise }
虽然C?假设G的边有方向，但FDM生成的几何形状对此是不变的。因此，只要任意两个顶点之间只有一条边相连，边的方向可以任意选择。

对于三维结构，连通性矩阵扩展为一个块对角矩阵C∈{?1,0,1}3?×3?v，块的数量等于笛卡尔维度的数量[113]：
(16) C = I??C?
其中I?∈R3×3是三维情况下的单位矩阵，?表示克罗内克积[114]。

2.3.2. 边界条件
为了施加边界条件，我们选择图G中的一部分顶点ns使其在平移方向上固定。这些固定顶点代表结构与地面锚定的位置，也是反作用力传递到结构基础的位置。我们假设固定顶点在所有三个维度上都被约束。其余的nu个顶点是自由顶点。悬链线几何形状就是在这些顶点上定义的。施加边界条件的影响体现在连通性矩阵C上，C现在沿着其列分为Cu和Cs两个子矩阵，分别对应自由顶点和固定顶点的自由度：
(17) C = [CuCs]
通常，C的每一列对应于一个顶点在一个笛卡尔维度上的自由度（如方程(16)所示）。虽然我们处理的是在所有三个维度上都受限的固定顶点，但C的分割可以更细致地进行，以选择性地在一个、两个或三个维度上约束一个顶点，从而提高我们方法的通用性。这种分割包括将每个连通性块C?分为自由和固定子矩阵，即C?=[C?uC?s]。然后将所有得到的子矩阵组装到C中：
(18) C = [[C?uC?s]x??? [C?uC?s]y??? [C?uC?s]z]
其中下标x、y、z表示每个笛卡尔维度的单独连通性块。然后根据方程(17)选择自由度和固定度的列来进行悬链线几何形状的计算，以确保与方程(20)和(22)的兼容性。需要注意的是，在特定维度上施加平移约束时，区分自由度和固定度比使用自由顶点和固定顶点的约定更为准确。

2.3.3. 仿真参数
FDM由定义在G的顶点和边上的三组属性参数化：
• 力密度q∈R3?：一个力密度向量，用于建立每根杆的轴力f与长度l之间的比率。q的符号决定了杆系统平衡时轴力的方向。负的q表示压缩轴力，正的q表示拉伸轴力。
• 施加的载荷p∈R3??：作用在自由顶点上的静态外力。由于结构被表示为杆系统，载荷只能施加在顶点上，而不能施加在边上。
• 支撑位置s∈R3?s：固定顶点的三维坐标；边界条件。
为了优化，我们可以将输入参数扁平化并连接成参数向量θ?=[q?p?s?]∈R3?p。输入FDM的参数总数np为np=ne+3ν+3ns。

2.3.4. 悬链线几何形状的计算
悬链线几何形状x定义为x?=[x???x???]∈R3?，完全通过轴力承受外载荷。这一属性源于在系统的变形配置中找到平衡状态。如2.2节所述，当能量Π达到最小值时达到平衡，这等同于自由顶点上的残余力r?=[r???r??]∈R3??的范数为零：
(19) r(x,θ) = K(θ)x ? b(θ) = 0
这里，K∈R3?×3?是系统的刚度矩阵[18]、[113]：
(20) K(θ) = Cu?QCu
其中Q∈R3?×3?是一个块对角矩阵，其元素为力密度q，Q=I??diag(q)。有效载荷向量b∈R3??从载荷p中减去连接到支撑点的边力的贡献：
(21) b(θ) = p ? d(θ)
(22) d(θ) = Cu?QCs
由于刚度和载荷与几何形状无关，因此给定θ时x的计算归结为解线性系统：
(23) x(θ) = K(θ)?1b(θ)
得到的x是由平衡状态下自由顶点的三维坐标定义的悬链线几何形状。方程(23)可以通过共轭梯度方法[115]迭代求解，或者直接通过矩阵分解方法（如乔列斯基分解）求解。如果杆系统的力密度既有正值也有负值，K的主对角线上可能会出现零值，导致奇异性，从而在桁架和张力结构的形式寻找过程中引入数值挑战。在这种奇异系统中，可以通过添加刚体特征向量的外积的缩放版本来正则化K[112]、[116]来处理方程(23)的求解。然而，这些正则化技术的应用超出了本工作的范围。

2.3.5. 完成平衡状态
使用FDM获得悬链线几何形状x后，可以计算出完成系统平衡状态的其余特征。它们包括：
• 边长l∈R3?：全局位移和局部变形后的杆长度。
• 边力f∈R3?：基于q的符号，杆的内部轴向力的大小（拉伸或压缩）。在结构力学中，这些力对应于伸长杆中的柯西力，考虑了泊松效应。
• 支撑反作用力t∈R3?s：支撑点处的力向量，即结构将对其基础施加的力。
上述属性通过以下公式计算：
(24) l(x,θ) = (U⊙U)1
(25) f(x,θ) = q⊙l
(26) t(x,θ) = ps ? Cs?Qu
其中u=Cux+Cs∈R3?，U=matrix(u)∈R3?，包含边坐标差；1∈R3是单位向量；ps∈R3?是支撑点处的载荷（如果有的话）；⊙表示哈达玛积[114]。

2.4. 在形状依赖载荷下的平衡
FDM使用方程(23)为恒定载荷生成悬链线几何形状。然而，载荷的不变性往往无法捕捉到载荷本身依赖于结构形状时的机械行为。这是因为在形式寻找过程中，杆系统必须在施加的载荷下变形以达到平衡，而这种变形又依赖于载荷。一个典型的形状依赖载荷的例子是自重，通常被认为是轻量化结构概念设计中的主要驱动因素[2]。虽然可以使用自重的初始估计来计算悬链线几何形状，但这些载荷很少会根据杆系统的变形更新以反映形状和质量分布的变化。图5通过一个简单示例说明了不更新载荷的后果：在恒定载荷下，悬链线拱的高度大约是其更新重量下真实平衡几何形状高度的一半。因此，不更新载荷可能会导致生成的几何形状不准确。

下载：下载高分辨率图像（322KB）
下载：下载全尺寸图像

图4. 一个FDM步骤的计算图。该图详细展示了在给定θ和G的情况下，计算悬链线几何形状x所涉及的线性代数层面的主要数值操作或原语νi，包括前向仿真和导数反向传播的过程。为了完整性，我们还指出了在多个FDM步骤求解形状依赖载荷下的平衡时载荷p依赖于先前几何形状x(t?1)的情况。反向模式AD和解析伴随矩阵提供了x关于θ的导数。反向传播路径由虚线箭头表示。

下载：下载高分辨率图像（384KB）
下载：下载全尺寸图像

图5. 平面拱的升起高度zmax随力密度q的变化。结构的两个极端顶点被固定，其余顶点可以自由移动。施加在自由顶点上的初始垂直载荷为?1/2。图5a表明q与zmax之间的关系是非线性的，并且在q=0附近表现出渐近行为。图5b和c展示了q从?1变化到?2时，zmax的变化幅度有几个单位。受形状独立载荷影响的结构比受形状依赖载荷影响的结构更矮。

为了克服这些挑战，我们扩展了FDM以支持在形状依赖作用下的几何形状计算。我们通过迭代生成悬链线几何形状来实现这一点，与其他形式寻找方法[27]、[117]保持一致。载荷p对几何形状x的依赖性在方程(19)中的b部分体现出来：(27)r(x,θ)=K(θ)x?b(x,θ)=0，其中b(x,θ)=p(x,θ)?d(θ)。我们不直接求解方程(27)的根，而是通过应用更新后的参考策略（URS）[6]将任务分解为一系列线性平衡问题。URS是一种通过迭代更新杆系几何形状来求解非线性形状优化问题的通用框架，直到达到平衡约束的固定点。这种策略的优势不仅在于其快速收敛性，还在于每次迭代都能确保系统的平衡（见图8）。然而，与原始的URS不同，我们保持力密度不变，只改变载荷。

下载：下载高分辨率图片（239KB）
下载：下载全尺寸图片

图6. 生成形状依赖载荷的汇流区计算。在每次形状优化迭代t中，归属于顶点i的面k的面积ak是两个相邻三角形(xi, mj, ck)的面积之和，其中xi是顶点位置，mj是k和i共有的两条边j的中点，ck是面的质心。在这个例子中，与顶点i相连的四个面用灰色阴影表示，第一个面的a1值等于两个粉色三角形面积之和。白色三角形完成了与顶点i相邻的其他三个面的汇流区。该图还指出了面法向量nk，这些法向量是在根据载荷方向逆时针排列面顶点后得到的。

下载：下载高分辨率图片（381KB）
下载：下载全尺寸图片

图7. 在形状依赖载荷下的FDM计算图。通过从初始猜测x0开始迭代更新几何形状和载荷，直到在迭代τ时收敛到固定点，同时θ和G保持不变。注意初始几何形状和收敛后的几何形状之间的变化。对于设计，我们制定了一个包含目标g的损失函数L，该函数衡量x与代表设计规范的目标x^之间的差异，然后通过基于梯度的优化来最小化损失值L(θ)。梯度值?θL(θ)是用反向模式AD计算的，隐式微分仅通过图中的?x/?θ反向传播到最后一步FDM。

下载：下载高分辨率图片（753KB）
下载：下载全尺寸图片

图8. 在形状依赖载荷下的平衡收敛。我们使用URS最小化两个连续版本x之间的成对差异，反复求解方程(32)，直到找到一个固定点。这里，过程从t=0时的平面网格开始，我们可视化了几次迭代以说明进展。在每次迭代中，标注显示了顶点i处的载荷pi以及从前一次迭代的位置xi计算出的残差ri。当只有载荷大小依赖于形状时，在τ=6次迭代时达到平衡；当大小和方向都随几何形状变化时，在τ=12次迭代时达到平衡。

在每次迭代t+1中，我们根据最后的几何形状x(t)重新计算施加的载荷p(t+1)，并使用方程(21)更新b。p和x之间的关系原则上是任意的，甚至可能不在顶点上直接定义。我们关注两种类型的形状依赖性载荷：仅大小依赖和大小及方向都依赖。后一种类型也称为跟随载荷，这对于模拟可充气和空气支撑结构是必要的。我们在第3.2节和第3.4节中展示了这两种载荷的设计示例。

为了处理跟随载荷，我们利用基于面的面积载荷到基于顶点的载荷的离散投影[108]、[118]，使用依赖于杆系几何形状的汇流区。面是一个多边形表面块，其边界是一个具有三条或更多边的拓扑循环。面的边界顶点的位置在R3中构成一个封闭的多边形[119]。我们使用墙跟随者（wall-follower）图遍历算法[120]在G的一个实例上计算不重叠的面。

在迭代t+1时，顶点i的形状依赖载荷pi(t+1)是其邻居在迭代t时几何形状的函数：(28)pi(t+1)=∑k∈F(i)ak(t)ψk(t)，其中ψk(t)=ψknk(t)。这里，F(i)对应于与顶点i相邻的面；ak∈R是面k的汇流区；ψk∈R是面k上的载荷大小，nk∈R3是其单位法向量。我们用鞋带公式得到nk：(29)nk(t)=n?k(t)∥n?k(t)∥，其中n?k(t)=1/2∑i=1nkxi(t)×xi+1(t)，并且满足闭合条件xnk+1=x1。在上述方程中，nk是每个面的顶点数，xi是相对于载荷施加方向逆时针排列的面顶点的3D坐标。当载荷不是跟随载荷，并且更感兴趣的是根据顶点的汇流区面积来调整载荷大小时，我们用R3中的任意单位向量替换nk。在自重的情况下，我们将这个向量设置为垂直向下方向。我们通过计算连接顶点i到面质心ck和k的两个边界边在中点mj的三角形对面积之和来估计面k的汇流区ak。图6示例了这种计算。在迭代t时，质心ck是通过平均面顶点的当前位置xi得到的：(30)ck(t)=1/nk∑i=1nkxi(t)。同样，边j的中点mj（其端点为xi和xi+1）是它们位置的平均值：(31)mj(t)=xi(t)+xi+1(t)/2

一旦重新生成了载荷，我们就应用FDM通过方程(23)找到下一个形状依赖的绳索几何形状x(t+1)：(32)x(t+1)(x(t),θ)=K(θ)?1b(x(t),θ)。我们重复这个过程直到收敛。因此，一个在形状依赖载荷下平衡的绳索几何形状是由一系列从x(t)开始并在τ次迭代后收敛到固定点的几何形状组成的：(33)x=x(τ)(x(τ?1)(?(x(t+1)(x(t),θ),θ),θ)。当所有自由顶点在连续两次迭代t和t+1之间的3D坐标xi的平均成对差异Δx(34)Δx=∑i=1nu∥xi(t+1)?xi(t)∥小于或等于收敛阈值η时，URS收敛，即Δx?≤?η。图7展示了在计算图中递归计算绳索几何形状的示意图，图8说明了电缆网上对大小和大小及方向依赖性的URS收敛情况。

2.5. 使用基于梯度的优化解决逆向形状优化问题
本工作的目标是自动化生成满足建筑、制造和结构约束的轻量结构的高效形状。找到这样的可行且高效的形状需要解决一个逆向形状优化问题，因为绳索几何形状的空间X无法直接访问。必须通过执行正向结构模拟（由FDM提供）后，通过参数θ的变化间接探索这个空间。由于缺乏可解释性和这个空间的高维度，以及在销接杆系中的机械行为非线性（见图5），导航由θ跨越的设计空间以解决此类逆问题具有挑战性。

为了应对这些挑战，我们使用基于梯度的优化来解决逆向形状优化问题。这种方法涉及在迭代优化算法的内循环中反复执行FDM以最小化损失函数。我们选择这种方法是因为基于梯度的方法在导航高维、非线性性能空间以优化物理系统以达到目标行为方面已被证明是成功的[90]、[121]、[122]、[123]。基于梯度的优化方法也适用于欠约束或过约束问题，这在轻量结构设计中很常见，其中目标条件的数量与可优化的自由度不同[113]、[124]、[125]、[126]。

对于给定的设计任务，我们根据形状优化参数θ和绳索几何形状x，按照惩罚方法制定一个标量值损失函数L。损失函数衡量x与设计规范的匹配程度：(35)L(θ,x)=∑i=1ngαigi(x(θ))，其中gi是一个可微的、通常是非线性的目标函数，用于计算应该在x上最小化的感兴趣的属性。一个设计规范可能包含ng个不同的目标函数，它们相对于任务的相对重要性由非负权重因子αi控制。L的一般公式有助于组合不同的gi来定义自定义规范，正如我们在第3节中所展示的。然而，具体目标可能会影响优化的收敛速度，特别是对于相互冲突的目标，其中一个目标的改进会恶化另一个目标的性能，在这种情况下，最优解将在目标之间产生权衡。由于x(θ)是通过FDM依赖于θ的平衡状态变量，损失值L实际上只是θ的函数。因此，为了表示简便，我们在这里将L(θ)写作L(θ,x(θ))。

为了解决逆向形状优化问题，我们寻找一个最优的θ?，使其最小化L(θ)，从而产生一个足够接近设计规范的绳索几何形状x?。我们通过一个受限优化问题得到θ?：(36)minθL(θ)s.t.hi(θ)=0?i=1,…,nhmi(θ)≤0?i=1,…,nmθlb≤θi≤θub?i=1,…,np，其中r(θ)=0。这个问题可以受到nh个等式约束h、nm个不等式约束m以及参数的下限θlb和上限θub的限制；所有这些都预先限定了设计空间的可探索区域。我们明确列出平衡约束r(θ)=0，以确保得到的几何形状是绳索形状。只要K是可逆的并且URS收敛，FDM（第2.1节）就会强制执行这一条件。请注意，虽然θ的所有条目可以同时优化，但通常将优化限制在子集上以施加已知偏差，例如固定的支撑位置或代表自重或其他外部作用的预定外部载荷。

我们可以使用梯度下降法求解方程36，沿着负梯度值?θL(θ)∈Rnp的方向迭代更新θ，直到找到局部最优参数θ?：(37)θ(s+1)=θ(s)?μ(s)?θL(θ(s))。优化过程从初始猜测开始，直到收敛，如图9所示。高阶基于梯度的优化算法，如LBFGSB [127]、SLSQP [128]和Levenberg-Marquardt [129]，可以通过计算合适的步长μ或基于L结构的增强参数更新规则来加速优化。从数值上讲，当损失值满足L(θ)≤?或梯度范数∥?θL(θ)∥≤κ时，优化会收敛。

下载：下载高分辨率图片（566KB）
下载：下载全尺寸图片

图9. 梯度下降是一种在轻量结构设计中数值上高效解决逆向形状优化问题的方法。这种优化算法通过迭代更新参数θ和损失函数的梯度?θL(θ)来寻找尽可能接近目标x^的绳索几何形状x。损失值L(θ)量化了x满足设计规范的程度。在找到最小化L(θ)的最优参数θ?和最优几何形状x?后，搜索结束，通常有x★≈x^。

我们解决逆向形状优化问题的策略与其他工作的主要区别在于我们将力学和设计分离，将问题框架化为双层优化。内循环通过FDM求解平衡，产生一个最小化能量泛函Π的绳索几何形状x，外循环调整θ以使x满足L中的目标。原始的FDM公式[28]及其后续工作[113]、[125]、[130]将平衡约束和规范的目标合并到同一个损失函数中。然后，它们直接优化x，如果任务规范不可行且优化算法无法收敛，则既无法得到平衡的几何形状也无法满足目标。相比之下，我们通过分离力学和设计来建立层次结构，优化θ而不是x。这种分离允许FDM专注于在每次迭代中强制平衡以生成有效的x，而优化算法则在绳索几何形状的空间内寻找满足设计规范的解决方案。这两种方法在结构优化中都是众所周知的，分别是同时分析和设计（SAND）以及嵌套分析和设计（NAND）[131]。我们的工作遵循NAND方法。**具有自动微分的导数计算**

为了找到方程（35）的最小值，需要计算梯度?θL，它由损失函数对FDM参数θ的偏导数组成：
$$
\nabla\theta L = \frac{\partial L}{\partial x} \frac{\partial x}{\partial \theta} + \frac{\partial L}{\partial \theta}
$$
特别是，雅可比矩阵$\frac{\partial x}{\partial \theta} \in \mathbb{R}^{3nu \times np}$，在结构优化中表示状态敏感性，它连接了θ每个元素的偏导数：
$$
\frac{\partial x}{\partial \theta} = \left[ \frac{\partial x}{\partial q}, \frac{\partial x}{\partial p}, \frac{\partial x}{\partial s} \right]
$$
之前的工作可能提供了获取上述方程中导数的方法，但仅是初步的。FDM是一种已知的形态寻找方法，已有超过五十年的历史[28], [29]，并且从那时起已经为各种设计任务发布了许多分析导数表达式。例如，规定缆网的目标边长和力[29]，设计张力结构[67]，匹配桥梁设计的反作用力[132]；近似目标形状[37], [133]；减少二维桁架的柔顺性[134]；以及最小化仅受压缩的壳体的能量[107], [135]。虽然这些信息很有用，但将这些公式重新用于与其最初制定目的不同的任务需要谨慎的数学处理。此外，这些公式主要表示对力密度的导数（例如，Schek的原始工作[29]推导出了雅可比矩阵$\frac{\partial x}{\partial q}$，即方程（39）中的第一个项），使得优化剩余的FDM参数成为一个开放的研究课题。

我们没有为每种可能的设计任务生成分析导数表达式，而是在AD框架[136]中实现FDM来自动化这一过程。我们的完整软件实现可以在参考文献[137]中找到。在这里，我们关注实现背后的数学原理，这些原理可以推广到不同的AD框架[138], [139]，并在附录中提供了一个紧凑的实现。特别是，我们使用反向模式AD来自动获得最小化方程（35）所需的梯度值$\nabla\theta L(\theta)$，适用于任意的损失函数，以及所有的FDM输入：力密度$q$、施加的载荷$p$和支持位置$s$。同时优化力密度和施加的载荷或支持位置，可以找到载荷和支持位置，从而扩展索形设计空间，达到更低的损失值，从而获得性能更好的几何形状（第3节）。同时，AD使我们不必关心损失函数L，而不必为L的每个变化重新推导$\frac{\partial L}{\partial x}$和$\frac{\partial L}{\partial \theta}$的分析表达式。

反向模式AD的工作原理是首先将损失函数评估过程中执行的中间操作或原语存储在一个有向无环计算图中。每个原语$\nu_i$对应计算图中的一个顶点，根据数据流与其他顶点相连。然后AD从损失值反向传播导数，回到输入，每一步都应用为每个原语定义的向量-雅可比积（VJP）[70], [136]。这种反向传播，称为反向传播，可以在单次图遍历中计算出完整的梯度$\nabla\theta L(\theta)$，其时间复杂度与一次正向传播相同，与用于标量损失最小化的正向模式AD相比，后者需要为每个参数进行单独的正向传播[72]。

对于可以用单个输入和输出表示为简单原语序列的函数，反向传播算法可以表示为链式法则：
$$
\nabla\theta L(\theta) = \nu_{n-1} \cdot \left( \nu_n \cdot \frac{\partial \nu_n}{\partial \nu_{n-1} \cdot \frac{\partial \nu_{n-2}} \cdot \ldots \cdot \frac{\partial \nu_1}{\partial \nu_0} \right)
$$
其中，上游余切$\nu_i = \frac{\partial \nu_n}{\partial \nu_i}$是链的输出$\nu_n$直到原语$\nu_i$输出的累积导数。当余切是一个向量$v_i$且$\nu$的输入和输出也是向量时，反向传播序列中的每一步都是一个VJP：
$$
v_{i-1}^T = v_i^T \cdot \frac{\partial \nu_i}{\partial \nu_{i-1}
$$
在每一步中，余切与原语的雅可比矩阵相乘，将原语在其输入处线性化。这个VJP是这种线性化的伴随映射，将输出空间的敏感性映射回输入空间的敏感性[72]：
$$
\left[ \frac{\partial \nu_i}{\partial \nu_{i-1} \right]^T v_i = v_{i-1}
$$
反向传播依次应用伴随映射，将导数从链的输出拉回到其输入（方程（40））。当输出原语$\nu_n$是损失函数L时，$\nu_{n-1} = 1$，反向传播递归得到梯度$\nabla\theta L(\theta)$。将反向传播算法从计算链扩展到一般图（如图4和7中的FDM图）是直接的。大致来说，它涉及处理原语的分支和合并，并通过扇入和扇出规则对VJP进行求和，重复方程（40），但现在是对图中从输出到输入的所有路径。有关深入解释，请参阅[72], [140]。

反向模式AD的一个数值高效特性是它避免了为了计算梯度$\nabla\theta L(\theta)$而实际计算原语的完整雅可比矩阵。因为L是标量值，所以只需要执行一系列VJP就可以解决方程（40），将微分视为一系列矩阵-向量乘积，而不是更昂贵的矩阵-矩阵乘法。然而，反向传播的速度受到原语抽象层次的限制：从基本算术（加法、乘法和除法），到线性代数内核（在AD框架中很常见），再到特定领域的模拟操作符[141]。更高的抽象层次可以通过减少需要反向传播的VJP数量来降低计算开销[142]。在接下来的章节中，我们将开发特定于力学的伴随映射，以加速反向传播，从而加速DFDM的优化。

**2.7. 线性求解器伴随映射**

我们方法中的一个瓶颈是求解方程（23）和（32）中的线性系统。我们希望避免通过求解系统所涉及的所有操作进行反向传播，因为我们只需要在优化时解$x$处的导数，而不需要其任何中间步骤的导数。此外，尽管对于小问题我们可以实际计算完整的密集刚度矩阵K，但对于较大问题，我们希望利用其固有的稀疏性来提高计算速度，特别是在优化过程中多次求解线性系统时。

因此，我们应用隐式微分[52], [143]来推导线性系统求解器的伴随映射，以便与反向模式AD兼容。这个伴随映射隐藏了求解器的内部机制，使得反向传播更加高效，并允许我们在不改变微分逻辑的情况下灵活选择求解器。虽然之前的工作已经为结构优化开发了线性求解器伴随映射[95], [97]，但它们假设$x$仅依赖于线性系统的左侧K，而忽略了解依赖于右侧边界条件。在我们的情况下，这个假设相当于仅对力密度$q$进行微分。这里提出的伴随映射比之前的更全面，考虑了$x$对K和有效载荷$b$的依赖性。这种能力使我们能够优化$q$，以及施加的载荷$p$和支持位置$s$，所有这些都被$\theta$包含在内。

给定一个损失函数L，AD为我们提供了一个余切$v^T = \frac{\partial L}{\partial x}$及其上游导数，我们的目标是实现线性求解的VJP，以获得方程（38）中的第一项梯度：
$$
\nabla\theta L = v^T \cdot \frac{\partial x}{\partial \theta} + \frac{\partial L}{\partial \theta}
$$
其中$v^T \cdot \frac{\partial x}{\partial \theta}$是我们对伴随映射感兴趣的项。为了推导它，我们首先重述平衡约束（方程（19）），注意到$x$是通过这个约束在解$\theta$处的隐函数：
$$
r(x, \theta) = K(\theta)x(\theta) - b(\theta) = 0
$$
接下来，我们对平衡约束求全导数以找到$\frac{\partial x}{\partial \theta}$的表达式：
$$
\frac{\partial r}{\partial \theta} + \frac{\partial r}{\partial x} \frac{\partial x}{\partial \theta} = 0 \Rightarrow \frac{\partial x}{\partial \theta} = -\left[ \frac{\partial r}{\partial x} \right]^{-1} \frac{\partial r}{\partial \theta}
$$
$$
\frac{\partial x}{\partial \theta} = -K^{-1} \frac{\partial r}{\partial \theta}
$$
计算完整的雅可比矩阵需要求解每个$\theta$参数的线性系统，其中K在左侧。我们只需要知道VJP就可以用反向传播来评估$\nabla\theta L$；不需要完整的雅可比矩阵，因此我们利用伴随映射来简化计算：
$$
v^T \cdot \frac{\partial x}{\partial \theta} = -v^T K^{-1} \frac{\partial r}{\partial \theta} = -\lambda^T \frac{\partial r}{\partial \theta}
$$
其中我们定义$\lambda \in \mathbb{R}^{3nu}$为伴随系统的解：
$$
K^T \lambda = v
$$
由于K是对称的，我们可以跳过转置，使用与正向传播中相同的线性求解器来求解方程（47），只是右侧不同。一旦我们在使用方程（23）进行模拟时找到$x$，并在反向传播中找到$\lambda$，我们就应用AD来自动计算VJP $\lambda^T \cdot \frac{\partial r}{\partial \theta}$，从而快速通过FDM进行反向传播。

**2.8. 通过形状依赖载荷的高效反向传播**

在形状依赖载荷下计算索形几何形状$x$包括重复求解方程（32），直到收敛到一个固定点（第2.4节）。虽然标准AD可以通过所有迭代来反向传播以评估梯度$\nabla\theta L$，但展开这些迭代会导致计算成本增加和数值不稳定[144]，例如梯度消失[145]，这可能会阻碍优化。这些问题可以通过再次利用隐函数定理[52], [143]来缓解。

我们制定了第二个带有隐式微分的分析伴随映射，以便通过形状依赖载荷进行反向传播。这种方法的好处是仅使用来自URS最后一次迭代$\tau$的信息来评估$\nabla\theta L$，而无需知道找到它所涉及的任何中间迭代（图7）。设$x$是方程（32）的一个固定点。我们重写这个方程以暴露$x$对$\theta$的依赖性：
$$
x(\theta) = K(\theta) - 1b(x, \theta) = y(x, \theta)
$$
其中$y$是固定点处平衡解的简写。我们对两边求导以获得固定点处的雅可比矩阵$\frac{\partial x}{\partial \theta}$：
$$
\frac{\partial x}{\partial \theta} = \frac{\partial y}{\partial \theta} + \frac{\partial y}{\partial x} \frac{\partial x}{\partial \theta} \Rightarrow \frac{\partial x}{\partial \theta} = \left[ I - \frac{\partial y}{\partial x} \right]^{-1} \frac{\partial y}{\partial \theta}
$$
其中$I \in \mathbb{R}^{3nu \times 3nu}$是单位矩阵。我们寻找$\frac{\partial x}{\partial \theta}$和余切$v^T = \frac{\partial L}{\partial x}$的乘积的表达式用于反向传播。这个VJP由下式给出：
$$
v^T \cdot \frac{\partial x}{\partial \theta} = v^T \left[ I - \frac{\partial y}{\partial x} \right]^{-1} \frac{\partial y}{\partial \theta} = \lambda^T \frac{\partial y}{\partial \theta}
$$
其中$\lambda \in \mathbb{R}^{3nu}$是伴随线性系统的解：
$$
[I - \frac{\partial y}{\partial x}]^T \lambda = v
$$
我们利用方程（48）的结构，并专门化伴随映射以获得更明确的表达式。雅可比矩阵$\frac{\partial y}{\partial x}$简化为：
$$
\frac{\partial y}{\partial x} = K^{-1} \frac{\partial p}{\partial x}
$$
因为K和$d$与$x$无关（第2.3节），并且雅可比矩阵$\frac{\partial p}{\partial x}$编码了载荷对几何形状扰动的敏感性。然而，在方程（52）中找到$\frac{\partial y}{\partial x}$需要为$x$的每个条目单独进行线性求解。我们可以再次使用伴随映射来避免这一点，因为我们不需要实际计算这个雅可比矩阵，而只需要知道如何与之相乘。如果展开方程（51）的左侧，线性系统变为：
$$
\lambda - \frac{\partial p}{\partial x}^T \alpha = v
$$
其中$\alpha \in \mathbb{R}^{3nu}$是以下方程的解：
$$
K^T \alpha = \lambda
$$
从这一点开始，我们分两步计算URS在固定点的目标伴随映射$v^T \cdot \frac{\partial x}{\partial \theta}$。首先，我们通过迭代求解器（如共轭梯度[146]）求解方程（51）和（53）来找到$\lambda$。这种方法允许我们通过重复应用VJP $\frac{\partial p}{\partial x}^T \alpha$来求解线性系统，而不是在内存中实际计算载荷的雅可比矩阵，后者会更加昂贵。找到$\lambda$后，我们再次利用AD来计算方程（50）中的$\lambda^T \cdot \frac{\partial y}{\partial \theta}$，从而获得我们所需的URS伴随映射。

**3. 结果**

在解释了DFDM如何与基于梯度的优化结合工作之后，我们现在在四个设计任务中展示了它的效果。每个任务都在不同的设计规范下突出并基准测试了DFDM的不同特性。所有呈现的数值实验都在具有64 GB RAM的M2硅芯片的一个CPU线程上运行。除非另有说明，我们使用直接稀疏求解器来解决涉及K的所有线性系统的前向和后向传递。硬件和求解器的选择基于我们在附录中呈现的计算扩展研究。我们设置$\eta = 10^{-6}$，并使用LBFGSB [147]，一种基于梯度的拟牛顿算法，以$\epsilon = 10^{-6}$和$\kappa = 10^{-6}$的收敛阈值来解决优化问题，除非我们在任务描述中指定了不同的配置。对于所有任务，我们显示由收敛最快的优化方法生成的解决方案。

**3.1. 匹配目标形状**

计算近似任意形状的索形几何形状有助于设计轻质结构，如网格壳和薄壳，同时保持视觉意图[27], [37], [77], [99]。这个任务也可以被构建为一个校正问题，其目标是找到任意形状必须经历的变形，以减少其总应变能量，从而提高其机械效率[148]。解决形状匹配问题也适用于评估传统建筑的机械完整性。其中一个例子是分析未加固的砖石拱顶，它们主要在轴向压缩应力下抵抗外部载荷[31], [33], [49], [149]。我们匹配目标形状的方法涉及选择一个目标表面并定义一个图G来近似它。然后，我们最小化每个顶点xi在平衡状态下的3D坐标与其目标位置x^i之间的平均距离：(55)g1(θ)=1nu∑i=1nu∥xi(θ)?x^i∥2 在平衡状态下的位置与方程(55)中的目标位置完全匹配时，该值为零。然而，对于任意形状来说，这种情况不太可能发生，因为这样的匹配意味着目标从一开始就是悬链线形状。因此，这项任务的目标是生成一种悬链线几何结构，使得自由顶点在最小二乘意义上最接近目标形状。我们近似图10中显示的三个双曲面形状，以评估DFDM的准确性和速度。所选形状在之前的研究[31]、[106]中出现过。它们具有同向和反向曲率，并且在其周长上得到支撑。优化问题的表述如下：(56)minqL(θ)=g1(θ)s.t.qlb≤qj≤qub?j=1,…,ne 我们特别寻找能够产生压缩轴向力的力密度q，并通过为其设定负界限来实现这一点。我们使用qlb=?50和qub=?10?3作为每个力密度的下限和上限。我们在目标形状上施加一个单位面积载荷ψ=1 kN/m2作为自重，并使用顶点在目标上的汇流面积（方程(28)）将其分配到各个顶点载荷上。为了模拟之前研究的实验设置，我们将模拟限制在一个FDM步骤内，并在优化过程中保持顶点载荷不变，使用阈值?=κ=10?5。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像图10. 形状匹配任务。我们用仅限压缩的替代方案近似三个任意双曲面形状。第一行显示了优化过程之前的悬链线几何结构，这些结构是用每个边的恒定力密度q生成的，并且在一个FDM步骤中完成。第二行展示了优化后的几何结构。为了完成任务，最优力密度需要有一定的分布，这可以通过彩色边来看出。如标注和表1中详细说明的，平衡状态下各个顶点位置xi与目标位置x^i之间的距离很小，优化后的解与期望的形状非常吻合。表1总结了结果。我们报告了所有边的初始力密度q0、开始时的损失L值、最优值、弹性应变能量Γ以及悬链线几何结构与目标形状之间的Hausdorff距离dH。Hausdorff距离通过测量两组点之间的最大成对距离来衡量它们的相似性。我们将其标准化，并以目标边界框对角线长度的百分比表示。表1还报告了离散化细节、优化迭代次数和总收敛时间。表1. 形状匹配实验的统计数据。Shapenenu?q0L0↓L↓dH?↓Γ^ax↑Γax?↑Γ/Γ^↓迭代次数 ↓时间 ↓时间/迭代次数 ↓提升 ↑空白单元空白单元空白单元空白单元空白单元空白单元[%][%][%][%]空白单元[s]AD [ms]ADA [ms]空白单元Bumps4622046.021.10.032.243.287.440.92770.32.00.92.2?×?Tristar7993464.52.10.021.023.196.524.11750.35.11.53.4?×?Cone12005293.01.00.010.942.698.435.8830.211.32.25.1?×?以图10a中显示的几何结构为起点，我们找到了将损失值至少减少两个数量级的最优力密度。图10b显示了最优的悬链线几何结构。在每个示例中，最优力密度遵循非平凡的分布以满足设计规范，这突显了自动寻找这些密度的重要性。我们的方法能够在不到半秒的时间内为所有三个目标生成紧密的仅限压缩的匹配。x与目标x^之间的最大偏差，由dH表示，对于Bumps仅为2.2%，对于Tristar和Cone则在1%以内。我们认为第一个形状相对于其他两个形状的拟合较为宽松，是因为前者具有尖锐的折痕，而后两个形状没有这种曲率不连续性。准确性对于找到高效且可行的形状很重要，但速度对于快速设计迭代也同样重要。为了量化这一点，我们将基于梯度的方法的 wall time 与两种代表性的无梯度方法进行了基准测试：Powell算法[150]和模拟退火[151]。我们使用这两种算法的标准实现和默认参数[152]，并将总优化迭代次数限制在1000次以内。图11显示了所有优化方法和三个目标形状的收敛历史。Powell算法的收敛损失值几乎与DFDM相当，但耗时超过了一千倍。模拟退火是所有算法中最慢的，仅在最小问题（Bumps）中在迭代限制内收敛。通过利用导数信息来指导搜索，基于梯度的DFDM方法比无梯度方法提供了显著的速度提升。下载：下载高分辨率图像（294KB）下载：下载全尺寸图像图11. 三个目标形状匹配任务的收敛历史。标注部分放大了优化过程前第一秒内的损失值变化。基于梯度的（GB）优化使用LBFGSB比无梯度（GF）方法（Powell方法和双退火）快三个数量级。在基于梯度的优化中，利用AD和解析伴随（ADA）的DFDM方法始终比通用AD更快。我们进一步通过将我们开发的线性求解器伴随（第2.7节）与公共AD框架[136]提供的实现进行比较，来提炼利用该伴随的计算优势。如表1所报告的，我们在每次迭代的优化速度和总wall time方面都优于通用AD，快两到五倍（图11）。特别是，随着优化问题的规模增大，基于梯度的优化变得更快，这对于在更细的离散化下建模结构非常有益。最后，我们根据结果设计的弹性应变能量Γ[148]来确定机械效率的改进。为了测量这种能量，我们对目标形状和优化后的几何结构进行了线性弹性结构分析，应用了与形状寻找过程中相同的载荷和边界条件。我们为每个结构的边缘分配了由名义钢材S235制成的管道截面。管道的直径是每个设计平均边长的十分之一，壁厚是直径的十分之一。优化后，总能量Γ中的轴向力份额Γax对于Bumps超过了87%，对于Cone和Tristar超过了96%，是优化前目标份额Γ^ax的两倍，对于最后一个形状则是四倍。这种行为证实了DFDM生成的几何结构在轴向力方面的主导地位，即使在不同的线性弹性力学模型下也是如此。此外，所有这些几何结构都将目标形状的总应变能量至少减少了一半。Tristar的减少幅度最大，为75%，而Bumps最小，约为60%。所有设计中的能量减少都强调了几何形状在结构性能中的重要性，以及微妙的几何变化可以带来巨大的效率提升。我们使用DFDM自动且几乎实时地生成了这些修改。3.2. 设计受压力载荷作用的结构设计轻质结构与流体相互作用涉及预测流体的压力如何影响结构的几何形状和内部力。由于压力载荷依赖于结构的形状，反之亦然[103]，因此模拟这种行为很复杂。有限元模拟可以准确地表征这些系统[104]、[153]、[154]，但需要二阶方法和精确的材料模型，增加了计算成本和问题的维度。将这些模拟集成到优化方案中会增加数值工作量[155]、[156]。像DFDM这样的形状寻找方法为设计过程的概念阶段提供了一种更快、与材料无关的替代方案，用于模拟这种类型的流体-结构相互作用[108]、[118]。这种能力使得DFDM非常适合探索可行的初步设计，特别是对于那些没有封闭形式解来表示的形状[157]。本工作的第二个任务转向设计一种空气支撑结构。空气支撑结构是一种单层气动系统，其刚性来自膜的內部加压。在长跨度上，这种结构通过缆网进行加固[157]、[158]。这些张力结构由于安装速度快且能够用仅几毫米厚的膜封闭大面积区域，因此在临时建筑和工业大厅及体育设施的屋顶中得到应用。在这里，我们考虑设计一个模型为半环面的气动结构的加固缆网（图12）。该缆网由24根径向缆绳和12根切向缆绳组成。每根径向缆绳被离散为12条边，而切向缆绳被离散为13条边，从而形成一个图G，具有ne=600条边和nv=312个顶点。缆网的顶点分别固定在半径为r1=3米和r2=9米的环面的内边界和外边界上，共有ns=48个固定顶点。其余顶点可以自由移动。我们施加一个内部压力载荷ψ=0.25 kN/m2，垂直于结构表面，忽略自重和外部载荷[158]。我们将结构表面的内部压力载荷转换为等效的跟随载荷，通过汇流面积施加到缆网的顶点上[108]、[118]。在正向模拟过程中，我们根据第2.4节中介绍的方法迭代更新缆网几何结构和载荷。我们设置URS的η=10?4。在这个任务中，我们在反向传播过程中仅使用第2.8节中的URS伴随，以将其计算性能与第2.7节中开发的线性求解器伴随区分开来。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像图12. 气动张力结构的设计。这种结构是一个半环面，其内边界和外边界的半径分别为r1和r2。我们计算了这个充气系统的加固缆网的几何形状，使其能够在仅受轴向拉力作用下抵抗0.25 kN/m2的内部压力。此外，其径向缆绳被限制在一个平面Φ上，并达到目标长度l^。目标长度在两个边界值β1和β2之间插值。DFDM为三组不同的β值对计算的几何形状都满足施加的约束，在图中显示。径向边根据其长度l进行了着色。设计任务包括计算满足三个属性的气动形状：径向缆绳达到特定长度l^j，径向缆绳位于平面Φi上，整个结构在轴向拉力状态下抵抗内部压力的作用。为了找到满足这些规格的悬链线几何结构，我们同时最小化三个目标：(57)minqL(θ)=∑i=13αigi(θ)s.t.qlb≤qj≤qub?j=1,…,ne 优化问题用力密度q来表述，共有600个优化变量，受限于箱形约束qlb=10?6和qub=20，以确保缆网的轴向力为拉力。第一个目标衡量nrad=288条径向边中的每条边的长度lj与其对应的目标长度l^j之间的差异：(58)g1(θ)=1nrad∑j=1nrad∥lj(θ)?l^j∥2 我们通过标量因子βj来规定每个l^j。这个因子沿着环面外边界逆时针方向在参考值β1和β2之间线性插值（图12a）。第一根和第十三根径向缆绳分别代表目标长度的极限情况。例如，β1=1/2和β2=3/2分别导致第一根和第十三根径向缆绳的边具有目标长度l^j=1/2和l^j=3/2；第七根缆绳的目标长度为l^j=1。如果β1=β2，则所有径向缆绳的目标长度相同。第二个目标量化每个自由顶点的位置xi与约束平面Φi=(oi,ni)之间的平方距离。该平面由一个沿环面边界等距分布的原点oi和一个切向于环面的单位法向量ni定义：(59)g2(θ)=1nu∑i=1nu∥(xi(θ)?oi)?ni∥2 在上述方程中，xi与Φi之间的距离是从oi到xi的向量在平面法线ni上的垂直投影的范数。第三个目标是一个L2正则化项，用于促进力密度均匀分布的解：(60)g3(θ)=∥q∥2 我们通过改变β1和β2的大小来探索三种设计选项。在第一种选项中，我们设置β1=β2=3/4；在第二种选项中，我们设置β1=β2=3/2；在第三种选项中，我们结合了前两种选项的系数β1=3/4和β2=3/2。在所有三种变体中，我们分别将q0=2和q0=3/2作为径向和切向边缘的力密度的初始值。我们从平坦状态开始充气气动结构，此时自由顶点位置xi与固定顶点位于同一平面上。权重分别为α1=5、α2=1和α3=6×10^-5。解决方程(57)中给出的优化问题后，得到了表2中列出的结果以及图12中显示的绳索几何形状。优化后的几何结构中的径向边缘满足长度和平面度目标，误差在10^-4的数量级，除了第三个实验，其中β1和β2不同，平面度误差高一个数量级，但仍保持在小数点后三位以内。长度目标影响了加压结构的几何形状和内部力。前两个长度规格产生了径向对称的几何形状，而第三个实验产生的解虽然不对称，但仍与β1和β2之间的目标长度对齐。图12b-12d展示了期望形状与实际形状之间的良好一致性。

表2. 气动设计任务的统计数据。

β1 β2 L0↓ L↓ g1↓ g2↓ azmax(1) zmax(1) fmax(1) fmax(1) 迭代次数 ↓ 时间 ↓ 提升 ↑
空单元 空单元 [×10^-2] [×10^-2] [×10^-4] [×10^-4] [m^2] [m] [m] [kN] [kN]
空单元 AD [s] ADA [s]
3/4 3/4 13 4.7 3.5 2.2 0.2 33 3.4 3.0 3.0 1.6 1.6 49 20.2 2.1 9.6 ×3/2 3/2 46 3.4 2.9 5.1 1.9 56 6.8 6.6 6.6 2.0 2.0
47 21.6 3.6 6.0 ×3/4 3/2 25 2.5 3.4 2.3 18.8 47 5.4 3.0 6.3 1.6 1.4 178 7 1.6 13.9 5.1

这些实验的一个逻辑结果是，更长的长度目标会导致具有更大表面积的扩展结构。一个有些反直觉的结果是，在第三个实验中，第十三根缆绳的高度比第二个实验中的短12.5%，尽管两个实验中这根缆绳的目标长度是相同的（见图12，底部行）。我们将高度差异归因于第三个实验中较低的表面积，这降低了总施加的压力，以及径向缆绳在高度上的轮廓变化。为了进一步说明这一点，图12底部行显示的所有径向缆绳都向半环的中心凸起，除了当β1=3/4和β2=3/2时的第一根径向缆绳。由于设计规范只规定了目标缆绳长度，而没有明确限制它们在高度上的几何形状，因此在满足任务要求的同时，径向缆绳的轮廓可以有所不同。这一结果还表明，可能存在多个最小化上述损失函数的绳索解决方案，寻找更具体的设计需要规定关于缆绳高度或曲率的额外目标。径向缆绳的内部力在绝对值上对任务规范的不敏感性低于结构的几何形状。例如，图12底部行显示的缆绳的最大内部力范围在三个实验中最多只变化了0.6 kN，介于fmax=1.4 kN和fmax=2.0 kN之间。轴向力的减少与表2和图13中的数据一致，证明了DFDM与第2.8节中开发的基于分析的伴随方法（ADA）在处理形状依赖载荷时的速度优势。通过避免在找到与施加的内部压力平衡的解的过程中进行所有迭代的反向传播，我们的方法在最昂贵的优化运行中收敛速度提高了5.1倍，将wall时间从71.6秒缩短到13.9秒；在最便宜的情况下，收敛速度提高了9.6倍，将wall时间从20.2秒缩短到2.1秒。特别是，图13表明，尽管两种微分方法最终都找到了解，但ADA需要的wall时间要长得多，其收敛曲线在时间轴上比ADA的结果拉伸得多。这种性能提升不仅加快了每次设计的优化速度，还扩大了在形状依赖载荷下的设计探索范围，因为可以在标准ADA找到单个解决方案所需的相同时间内生成五到十个不同的选项。

图13. 气动设计任务的收敛图。所提出的微分方法通过添加用于反向传播的分析导数（ADA）来增强AD，其收敛速度始终快于现成的AD，因为它缩短了反向传递中的迭代次数。

3.3. 载荷查找
在本节中，我们展示了载荷查找的应用，即同时计算最优外部载荷和力密度。虽然这种方法具有通用性，但在这里我们展示了它修复裂缝砖石圆顶的能力。这个半球形圆顶在[159]中有所描述，半径为r=5.0米，厚度为0.5米。圆顶在其支撑点处发生水平位移δi∈R^2，导致结构的内侧和外侧出现裂缝（见图14a）。

图14. 使用载荷查找修复半球形砖石圆顶。由于支撑点的水平沉降，圆顶出现裂缝，将圆顶从中部分裂（见图14a）。我们生成了两个压缩推力网络和一组修复载荷，一个没有平滑处理，另一个对各个载荷因子λ进行了平滑处理，分别如图14b和c所示。载荷查找有助于设计不同的后张缆网几何形状，这些缆网都施加相同的修复载荷（见图14d和e）。

砖石壳体的机械行为受全局不稳定性的控制，而不是局部应力失效[160]，[161]。因此，可以通过极限分析[32]，[33]来评估这些结构，这是一种离散表示的压缩应力场，采用绳索几何形状。通过极限分析，可以确定砖石壳体的坍塌载荷下限[162]，前提是存在一个与施加的载荷平衡的压缩推力网络，并且满足Heyman安全定理的无张力、无滑动和无压碎条件[163]。靠近壳体中轴的推力网络可以减少内部压缩力的偏心并降低不稳定性。在参考文献中，圆顶推力网络的离散化是一个具有20条子午线和16个环带的径向网格，包含nv=321个顶点，其中ns=20个是固定的，ne=640条边。我们使用相同的配置作为表示推力网络的销接杆系统的输入图G，用于载荷查找。该结构的代表性密度为1 kN/m^3，其重量为ψ=0.5 kN/m^2。我们通过在半球形圆顶的内侧表面上计算它们的汇流面积后，将自重汇总到网络的顶点上来转换自重（使用方程（28））。在形状查找过程中，载荷保持不变（即在模拟过程中没有形状依赖性）。总自重为77.4 kN。

受到加固拱门方法[164]，[165]的启发，我们提出了一种非侵入性和可逆的砖石壳体修复系统，通过附加外部后张缆网来修复圆顶。该网会对圆顶施加力，形成新的压缩应力状态，阻止支撑点的沉降并防止裂缝扩展[166]。然而，我们没有同时优化缆网和推力网络，而是将寻找作用在推力网络上的最优载荷的任务与获得施加这些载荷的具体缆网的任务分开。这种分离使得实现修复解决方案更加灵活，因为不同的缆网几何形状可以对圆顶施加相同的载荷向量。因此，我们寻找一组最优的外部修复载荷，使推力网络位于砖石圆顶的中心。

为了解决这个任务，我们首先用λi∈R参数化每个外部载荷，λi是一个标量，用于控制其在圆顶中间表面上单位法线ni∈R^3上的大小：
(61) pi = wi + λini
在顶点i处施加的总载荷pi是自重向量wi∈R^3加上缩放后的外部载荷λini。

我们通过逐点匹配定义在半球面上的目标高度z^来将推力网络引向圆顶的内侧表面：
(62) g1(θ) = 1/nu ∑i=1/nu ∥zi(θ) ? z^i∥2
我们使用[33]，[159]中描述的目标场。在所有可能的保持推力网络在砖石壳体厚度范围内的外部载荷分布中，我们关注一种载荷大小均匀分布的分布，以减少预应力和缆绳直径的变化，从而便于构建后张缆网系统。因此，我们制定了第二个目标来控制载荷因子λ的方差：
(63) g2(θ) = ∑i=1/nu (λi ? λavg)^2
λavg = 1/nu ∑i=1/nu λi

为了使推力网络居中并控制载荷分布的均匀性，我们最小化这两个目标：
(64a) min q, λL(θ) = ∑i=1/2 αigi(θ)
(64b) s.t. txi(θ) = 0 ?i=1,…,ns
(64c) tyi(θ) = 0 ?i=1,…,ns
(64d) xi(θ) ? x^i = 0 ?i=1,…,nu
(64e) yi(θ) ? y^i = 0 ?i=1,…,nu
(64f) z^i ? zi(θ) ≤ 0 ?i=1,…,nu
(64g) zi(θ) ? z^i ≤ 0 ?i=1,…,nu
(64h) qlb ≤ qj ≤ qub ?j=1,…,ne
(64i) λlb ≤ λi ≤ λub ?i=1,…,nv

这些约束确保了DFDM计算的绳索几何形状的可行性。根据标准的推力网络分析[26]，[33]，绳索几何形状被限制在一个固定的水平投影范围内，以保持汇总的自重载荷与重力方向一致，并与它们在圆顶内侧表面上的分布一致（方程（64d）和（64e））。同样，几何形状被限制在砖石壳体的厚度范围内，以满足安全定理[160]，[163]。我们通过将自由顶点的z坐标限制在它们在半球形圆顶内侧和外侧表面上的垂直投影zlb和zub之间来实现这一要求（方程（64f）和（64g）。我们使用[33]中提供的圆顶内侧、中间和外侧表面的分析描述来指定这些几何约束。最后，我们将支撑点处的反作用力的水平分量txy设置为零，以抵消由自重和修复载荷共同产生的推力，这可能会在圆顶支撑点处引起额外的损坏位移（方程（64b）和（64c）。

优化问题的变量是力密度向量q和载荷因子向量λ。我们只优化网络的一个径向切片的变量，以保持对称性并加快收敛速度。这个选择包括一个子午线上的边缘和顶点，每个圆环一个边缘，总共48个自由度：32个力密度和16个载荷因子。其他边缘和顶点被约束为具有与优化后的径向切片中的对应部分相同的力密度和载荷因子。我们对力密度施加盒形约束，以确保只有压缩力的推力网络（方程（64h）和（64i））。力密度约束为qlb=?50和qub=?10^-3。此外，我们规定最小载荷因子λlb=0.25，并限制最大载荷因子λub=2.5，以控制修复干预的后张量。

我们使用方程（64a）定义的优化问题进行求解，分别考虑了载荷均匀分布的目标（α2=0）和不均匀分布的目标（α2=1/2），以捕捉平滑项对解空间的影响。为了进行基准测试，我们解决了相同的优化问题，但只考虑力密度而不考虑修复载荷，并将这种方法作为基线。由于推力网络受到恒定外部载荷的影响，我们仅使用第2.7节中开发的线性求解器伴随方法来加速梯度计算。然后，使用支持等式和不等式约束的基于梯度的顺序二次规划算法（[167）来解决方程64a。优化从[159]中报告的裂缝圆顶的力密度开始，所有载荷的λ0=1。形状目标的权重为α1=1。用于载荷查找的优化后的推力网络显示在图14b和c中，损失函数和两个目标的收敛历史显示在图15中。表3报告了关键统计数据。所有三种方法的优化在不到一百次迭代和五秒的计算时间内就达到了收敛，使用的是通用的AD算法。采用线性求解器伴随（ADA）的AD算法将计算时间几乎减半，降至三秒以内。
下载：下载高分辨率图片（264KB）
下载：下载全尺寸图片
图15. 使用负载查找方法对砌体进行改造的任务的收敛历史，包括有无负载因子平滑处理的情况。基线方法（FD）仅优化力密度。我们最小化了损失值（左侧）、几何拟合度（中间）和负载方差（右侧）。在最后一个图中，由于损失函数L中的方差项α2的权重为零，因此不会受到惩罚，其方差值会波动。我们展示它只是为了与LFS进行比较。

表3. 基线方法（FD）、负载查找方法（LF）以及带有平滑处理的负载查找方法（LFS）的改造任务统计数据。

| 方法 | L0 | L | g1 | g2 | λavg | λmax | dM | ||txy||/w | ?favg | 迭代次数 | 时间 | 提升率 |
|--------------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|---------|
| AD | | | | | | | | | | | |
| ADA | | | | | | | | | | |
| FD | 1.24 | 4.74 | | | | | | | | |
| ×100 | | | | | | | | | |
| ADA | 1.24 | 4.74 | | | | | | | | |
| ×10?3 | | | | | | | | | |
| ×10?3 | | | | | | | | | |
| ×10?1 | | | | | | | | | |
| kN | | | | | | | | | |
| % | | | | | | | | |
| % | | | | | | | | |
| kN | | | | | | | | |

这三种方法都将推力网络拉得更靠近圆顶的内表面，但改进程度取决于优化问题的具体变量。在负载查找方法中，损失值的减少更为显著，损失值降低了六个数量级。相比之下，基线方法仅将损失值降低了两个数量级（见图15）。在形状目标g1上的表现对于减少损失值和使推力网络与圆顶内表面对齐至关重要。基线解中，与内表面的最大距离dM为壳层厚度的39.1%。而负载查找方法将这一距离缩短到了3.8%，提高了十倍，显著改善了内表面的拟合度。然而，通过负载查找方法获得的负载分布不均匀，λmax=2.5的峰值与λmin=0.25在圆顶周长上交替出现（见图14b）。例如，最大的负载出现在靠近圆顶支撑点的顶点处，而最小的负载则出现在圆顶的顶部和三分之一处。这种不规则的模式最终可能导致由于预应力和缆绳直径的变化而使得缆绳网的建造成本增加。

在损失函数中加入平滑项g2可以缓解这一问题。结果负载分布使得负载因子的方差从0.49降低到0.22，平均负载从0.6降低到0.3，从而实现了更加合理的后张拉设计。图14c证实了负载分布的均匀性。负载从圆顶顶部的λmin=0.25逐渐减小到支撑点的λmax=1.4。值得注意的是，虽然平滑处理会降低匹配精度，并导致推力网络与中间表面的贴合度降低，但这种方法仍使dM相比基线方法提高了三倍（39.1% vs 11.5%），同时调整了解决方案以适应制造要求。平滑项实际上是一个可以调节的参数，根据需要调整负载的均匀性。后张拉复杂性与推力网络贴合度之间的权衡取决于具体项目条件，DFDM再次提供了一个探索这些问题的框架。

在更新后的砌体圆顶的压缩应力状态下，两种改造负载集都将网络边缘的平均内力favg提高了五倍。但由于这两种解决方案都只在支撑点处产生垂直推力（满足我们施加的水平反力约束（方程（64b）和（64c）），因此它们仅在支撑点处产生推力。如果没有这些负载，圆顶的支撑点仍然会承受至少五分之一的自重作为水平推力（见表3中的||txy||/w），从而保持圆顶水平位移的趋势并加剧裂缝扩展。

图14d和e展示了两种不同的缆绳网设计，用于改造圆顶，以突出负载查找方法的灵活性。其中一个从上方修复圆顶，另一个则从下方、即内曲面下方进行修复。选择哪种缆绳网取决于现场条件；例如，圆顶的外曲面是否可进入进行施工。使用负载查找方法进行改造设计对于修复砌体壳体非常有用。它能够发现将推力网络居中于砌体结构中的绳索几何形状，并同时考虑后张拉系统的制造要求。由于额外的自由度，负载查找方法能够生成更好的设计，从而实现更低的损失值。虽然这项任务主要集中在加固拱门方法[164]、[165]上，但负载查找方法也可以支持基于推力线和推力网络分析的其他改造策略，例如使用复合材料[168]、[169]来加固砌体结构。

3.4. 支撑点寻找
在最后一个设计任务中，我们展示了在具有多个目标的任务中同时优化力密度和支撑点位置的方法。为此，我们将一个任意四主导布局的刚性网格壳的几何形状平面化。网格壳是一种能够通过模仿连续膜的最小化弯曲行为来跨越大距离的梁网络[170]。带有钢构件的玻璃网格壳结合了透明性和材料效率，为博物馆[22]、购物中心[172]、[173]、[174]、交通枢纽[175]、[176]、[177]等公共建筑创建长跨度且自然照明的屋顶，通过重新利用现有结构来减少新建需求并促进城市密集化。

我们以大英博物馆庭院上的标志性网格壳为例来展示我们的方法。该网格壳的形状是通过分析生成的，以在内外边界之间平滑过渡[171]，并采用三角化布局进行离散化，从而获得平面表面，但没有绳索几何形状的限制[178]、[179]。由于其支撑条件允许结构在外围水平滑动，现有网格壳通过弯曲和桁架作用共同承受其重量。因此，任务的目标是提出一种替代设计，使其在无弯曲状态下支撑自身重量，同时减少支撑点的水平推力，并与平板覆层兼容。任务设置见图16a。

下载：下载高分辨率图片（2MB）
下载：下载全尺寸图片
图16. 使用支撑点寻找方法对绳索网格壳进行平面化处理。第一列显示了结构的绳索几何形状、支撑点上的反力以及玻璃面板的平面性。第二列展示了初始几何形状，这是对实际建造的网格壳形状的最佳拟合。另外两列展示了优化后的设计，这些设计最小化了外边界南北两侧的推力，并最大化了面板的平面性。支撑点寻找方法得到的仅受压缩的几何形状将推力减少了90%以上，且比最佳拟合和基线选项拥有更多的平板（平均偏度仅为?avg=0.2%），从而在降低制造成本和简化施工的同时保持了机械效率。

大英博物馆的矩形庭院南北方向长97米，东西方向宽73米[171]。网格壳在两个高度相同的平面边界上得到支撑。外边界位于庭院周边的建筑物上，而内圆边界半径为22米，位于庭院中心以北3米处，从而在南北轴上引入了不对称性。与现有屋顶不同，我们允许在优化过程中内边界的几何形状发生变化。我们的假设设置还假设东西侧的建筑物能够抵抗水平推力，而南北侧的建筑物必须尽可能无推力。我们采用的四主导网格壳布局是通过中轴网格划分生成的[180]。我们在布局中引入了八个高 valency 顶点，以促进负载向结构角落的传递，正如现有网格壳所做的那样[178]、[179]。最终得到的拓扑结构包含ne=599条边、nf=282个面和nv=317个顶点。30%的顶点得到支撑。每个面对应一个玻璃面板。

我们将这个设计任务定义为四个目标的组合。第一个目标量化了需要在nout=25个支撑点上抑制的各个反力向量ti的大小：
(65) g1(θ) = 1/nout ∑i=1^nout ∥ti(θ)∥2

第二个目标通过测量代表该面的边向量与面法线的正交程度来编码面板的平面性约束。对于给定的面k，我们计算单位面法线nk与每个归一化边向量ej之间的点积的绝对值，并对面的所有边E(k)求和。最小化这个和可以促使边位于面内：
(66) g2(θ) = 1/nf ∑k=1^nf ∑j=1 |E(k)| ∥nk(θ)?ej(θ)∥
其中|E(k)|表示边集E(k)的大小。我们添加了第三个目标作为正则化项，以降低顶点的公平性能量，使解决方案倾向于具有均匀内角的平滑几何形状：
(67) g3(θ) = 1/nv ∑i=1 nv χi(θ)

顶点i的公平性能量χi考虑了不同的 valency，如我们在四主导布局中所使用的那样，遵循Tang等人的处理方法[106]。用V(i)表示顶点i的一环邻居的有序循环，|V(i)|表示其大小，能量是其位置xi与其邻居位置之间的距离的函数：
(68) χi(θ) = {∥xi(θ)?xi+1(θ)+xi+3(θ)2∥2 + ∥xi(θ)?xi+2(θ)+xi+4(θ)2∥2
如果|V(i)|=4，则 ?χi(θ) = ∑j=1|V(i)|xi+j(θ)∥2；否则

第四个也是最后一个目标限制网格壳几何形状的面积，使其与现有屋顶的面积相匹配：
(69) g4(θ) = ∥a^?∑k=1^nf ak(θ)∥2
通过测量实际建造的目标面积a^=6275平方米与新几何形状的面面积ak之和来确保适当的比较，从而隔离几何变化对性能的影响，而不是面积差异。

在定义了目标后，我们设置了优化问题，将力密度和内边界上的支撑点位置作为优化变量：
(70) minq,δL(θ) = ∑i=1^4 αigi(θ) s.t. qlb≤qj≤qub ?j=1,…,ne δlb≤δi≤δub ?i=1,…,nin

我们通过将单位平移向量ni乘以位移因子δi来参数化nin=31个内支撑点的位置：
(71) si(δi) = s0i+δini
平移向量从每个支撑点的圆形边界质心指向，并位于边界平面内。然后将缩放后的平移向量加到初始支撑点位置s0i上，如图16a所示。

优化变量共有630个，所有变量都受到盒约束。我们设置qlb=?100和qub=?10^?3，以确保仅受压缩的内部力状态（因此形状向上凸起，而不是像仅受拉伸状态那样向下凸起），并设置δlb=?0.3r和δub=0.3r，以保持内边界与其初始形状相对接近。通过求解方程27，可以保证网格壳几何形状的无弯曲特性，正如我们在第3节中介绍的其他任务一样。我们应用面积为ψ=1 kN/m2的自重作为结构的自重，这在现有屋顶的重量范围内[178]，并使用URS（第2.4节）求解方程(27)。我们注意到，我们的工作与其他同时优化弯曲和平面性的方法[106]、[181]、[182]、[183]、[184]不同，因为它们产生的形状可能是平面的，但不一定是无弯曲的。

我们比较了支撑点寻找方法与仅优化力密度的基线方法的性能。这三种方法都应用了DFDM和基于梯度的优化，但使用了不同的损失函数结构和优化变量。在最佳拟合方法中，我们计算了一个简单的绳索几何形状，该形状遵循我们在第3.1节中描述的程序来近似实际建造的网格壳形状。得到的力密度被用作支撑点寻找和力密度基线的初始猜测值，两者都求解方程(70)，目标权重设置为α1=10^?3、α2=12、α3=1、α4=10^?5，迭代次数恰好为一千次。我们固定迭代次数以提炼在相同计算预算下将支撑点位置作为额外优化变量的好处。

图16展示了优化后的网格壳、它们的承载行为和平面性。表4报告了几项统计数据，图17显示了支撑点寻找方法和基线方法的优化收敛历史。将线性求解器和URS伴随矩阵（ADA-LU）集成到反向传播中，与标准AD相比，可以将计算成本降低近四倍。然而，在支持寻找和基线方法中，大部分节省的成本来自于URS伴随矩阵。如表4所示，线性求解器伴随矩阵（ADA-L）最多只能将计算时间减少16%。由于URS伴随矩阵在前向和后向传播过程中多次调用线性求解器，因此它主导了计算成本。优化后的几何形状与目标形状之间的面积差异在2.5%的相对误差范围内，而分析得到的真实形状与三个优化几何形状之间的Hausdorff距离dH小于真实形状边界框对角线长度的4%。这些指标表明，优化后的几何形状在全局和局部上都高度逼真地再现了实际建造的网格壳结构。最佳拟合的几何形状具有最高的轮廓zmax，其次是支持寻找和基线变体，分别浅0.4米和0.8米。

表4. 网格壳设计任务的统计信息，包括最佳拟合（BF）解决方案、基线（FD）和支持寻找（SF）。

图17. 带有支持寻找功能的网格壳任务的收敛历史。我们在固定的1000次迭代计算预算下优化损失值（左）、反应力大小（中）和面板平面度（右）。我们将力密度优化（FD）与支持寻找（SF）进行了对比，在后者中，我们同时修改了力密度和支持位置。支持寻找方法得到的解决方案的反应力与基线相同，并且在平面度和最终损失值方面表现更好。支持寻找后，网格壳的内边界从圆形变为更复杂的椭圆形，并在庭院的南北两侧被压平（图16c）。这种变形是各向异性的，最大程度仅为内边界半径的五分之一：北侧和南侧的支持点向外移动，最大位移为δ=0.2r（位于顶点6和12处），而东侧和西侧的支持点向内移动，最大位移为δ=?0.1r（位于顶点1和17处）。这些变化响应了任务要求，在需要最小化推力的方向上减小跨度，在允许有推力的方向上增大跨度，以满足目标面积约束。

尽管内边界形状不同，但通过支持寻找和力密度基线优化的两种几何结构都将北侧和南侧支持上的总反应力（表示为结构总重量||t||/w的百分比）从最佳拟合方案的46.1%降低到了6.9%。这种减少是通过重新引导作用在网格壳角之间的杆件上的力传递来实现的，激活了内部拱形作用。最大杆件力fmax增加了四倍，平均杆件力favg增加了一倍（表4），以及角部更大的推力（图16，第二行），这些都证明了这种新的承载行为。拱形机制伴随着形状变化：两种优化后的几何结构的纵向轮廓在接近支持点之前向上凸起（图16b和c）。

尽管支持寻找和力密度基线在计算时间上仅相差一秒，但它们在面板平面度方面的输出差异显著。如果我们测量面板k的偏度?k（对面顶点之间的距离与面边平均长度的比率），并将其限制在1%的严格范围内，那么通过支持寻找得到的几何结构的面板最为平坦，达到96.8%，比最佳拟合和基线分别提高了36.9%和4.6%。同样，通过支持寻找得到的几何结构的面板平均平坦度是最佳拟合结构的五倍（?avg=0.2%），但最大变形仅限于六块面板。基线方法在面板平面度方面也有显著改进，但平均而言每个面板的偏度是支持寻找方法的两倍，在最坏情况下，如果选择后者进行制造，生产成本会更高。图17揭示了这种差异。单独的力密度优化在优化过程早期就达到了平面度目标g2的稳定状态，陷入了损失函数的次优区域。而将支持位置作为额外的优化变量，则能够持续降低损失值，最终超过基线水平。

总之，支持寻找方法产生的几何结构没有弯曲，且面板更加平坦，同时满足了推力、均匀性和面积约束。它提供了一种网格壳设计，减少了制造工作量，并且与现有场地条件兼容。通过同时校准力密度和支持位置，DFDM找到了在结构性能和制造便利性之间取得平衡的几何形状。

4. 结论

轻质结构能够实现长跨度系统，支持人类活动的许多方面，从跨越河流的桥梁到为公共空间提供遮蔽的屋顶。然而，设计这些结构的形状以高效承载负载并满足其他技术要求一直是一个长期挑战。为此，我们开发了DFDM，这是一种几何非线性的、完全可微分的方法，用于解决作为销接杆系统建模且在形状依赖载荷作用下的逆向形态寻找问题。通过利用AD和特定于力学的伴随矩阵，DFDM成为了一种通用方法，用于计算符合建筑、制造和结构要求的缆索结构形状，加速了轻质结构在大位移范围内的高效且可行的设计。结合基于梯度的优化，DFDM的设计速度比无梯度算法快3到4个数量级。开发的伴随矩阵结合了分析和自动导数，进一步提高了DFDM的速度，将优化计算时间相对于标准AD减少了2到10倍。它们展示了利用DFDM力学模型的计算图来开发特定领域AD实现的价值，从而加快了设计速度。DFDM的端到端可微分性支持了多种结构的设计，包括气动拉伸系统、压缩砌体壳和未受应力作用的网格壳，每种结构都有独特的任务要求、损失函数和高维参数化。通过负载寻找和支持寻找，我们发现了创新的缆索解决方案，用于修复损坏的半球形穹顶，以及设计一种与平面玻璃覆盖层兼容的无弯曲格子屋顶，同时减少了其对支撑建筑物的推力。

本文的工作扩展了经典FDM在可微分编程方面的能力。更广泛地说，DFDM促进了可微分形态寻找方法的发展，这些方法结合了稳健的力学模型和现代深度学习工具，加速了旨在满足人类需求的最小化材料结构的设计。

4.1. 局限性

本工作的局限性在于，所提出的设计任务仅考虑了单一载荷情况，而不是多种载荷场景。这种简化在概念阶段的形态寻找中很常见，尽管为初始设计选择最具代表性的载荷情况仍然是该领域的一个基本问题[185]。在这里，我们专注于最确定的一种载荷情况，即自重，除了气动设计任务，其中我们仅考虑了充气压力。尽管如此，DFDM仍然为后续分析和针对其他载荷情况的细化提供了坚实的基础。第二个挑战是，基于梯度的优化通常会在非凸损失函数中收敛到局部最优解，而不是全局最优解。这个最优解还取决于初始参数值的选择。正如我们的所有实验所证明的，如果参数初始化得当，局部最优解就足以准确解决设计任务。另一个缺点是损失函数中目标权重的校准是临时的，需要一些实验。这个问题可以通过利用DFDM的可微分性来解决，同时优化权重和损失函数值[186]，[187]，或者使用新的基于梯度的优化算法来无冲突地调整模拟参数，以达到损失函数中所有项共享的局部最优解[188]。

4.2. 未来工作

展望未来，DFDM可以从离散结构扩展到连续结构，以便直接设计膜结构和壳结构，正如其他研究对标准FDM所做的那样[25]，[35]，[42]，[43]。DFDM的结构力学基础为此扩展奠定了基础。虽然我们专注于解决具有多个目标和约束的问题，但我们将所有这些目标整合到了一个损失函数中，因此将开发的形状优化框架扩展到适当的多目标优化设置中，以提供帕累托前沿的解决方案将是下一步。通过使用二阶方法提取损失函数的曲率信息，可以进一步提高我们的方法的收敛速度，尽管计算Hessian矩阵每一步可能比计算梯度更耗费计算资源。确定这两种方法之间的权衡对于确定哪种方法更适合特定任务非常重要。最后一个有前景的方向是将像DFDM这样的可微分结构模拟方法集成到神经网络中，以进一步加速设计优化。在深度学习时代，可微分模拟有望将物理知识直接注入人工智能系统[93]，[96]，[189]，从而创建更快的设计替代模型。

CRediT作者贡献声明：

Rafael Pastrana：撰写——原始草稿、可视化、验证、软件、方法论、调查、概念化。
Deniz Oktay：撰写——原始草稿、软件、方法论、调查。
Kai-Uwe Bletzinger：撰写——审阅与编辑、验证、方法论。
Ryan P. Adams：撰写——审阅与编辑、方法论、资金获取。
Sigrid Adriaenssens：撰写——审阅与编辑、监督、资金获取。