自动微分驱动的参数并行增量谐波平衡法

《INTERNATIONAL JOURNAL OF MECHANICAL SCIENCES》：Automatic differentiation driven parameter parallel incremental harmonic balance method

　　

摘要自动微分（AD）使得增量谐波平衡（IHB）所需的雅可比矩阵无需手工微分或有限差分扰动即可构建。本研究开发了一种由 AD 驱动的参数并行 IHB 框架，称为 IHB-ADP，用于同时求解独立的非线性动力学子系统。从频域系数到残差的完整映射被封装为一个可微计算图，而响应重构、残差评估、雅可比矩阵构建和牛顿修正则作为 GPU 批处理张量操作实现。一种识别 - 退出机制独立管理子系统的收敛和失败，而一种扰动 - 重启机制避免了在孤立的不可微采样状态下评估牛顿雅可比矩阵。对范德波尔、达芬和 Mathieu 系统的测试显示，与其他 GPU 批处理实现相比，最小加速比约为 1.80，与串行基线相比为 7.92，证明了数值公式和参数级并行性的贡献。对于 Mathieu 系统，一种梯度引导的预测 - 偏差 - 校正策略沿着预定稳定性边界附近的交点构建参数集。对近零刚度隔振平台的动态测试显示计算响应与测量响应之间具有良好的一致性。最后，结合依赖坐标的傅里叶集和事件局部化的伯努利富集，分析了具有多尺度响应和单向接触的分层压电能量收集器。这些结果表明，IHB-ADP 为低至中等维度非线性动力系统的参数批处理分析提供了一个高效且灵活的框架。

引言非线性周期响应广泛存在于旋转机械 [1], [2]、航天结构 [3]、隔振器 [4]、能量收集器 [5] 和多体机械系统 [6], [7] 中。其分析通常涉及频率响应 [8]、共存解 [9]、稳定性边界 [10] 以及多个物理参数的影响 [11]。在工程应用中，相同的控制方程通常需要在数百或数千个参数组合上重复求解 [12], [13]。因此，整体效率既取决于为单个周期解采用的数值公式，也取决于用于完整参数集的执行架构。

时间积分和频域方法提供了互补的途径来进行非线性动力分析。时间积分直接重现瞬态演化，适合从指定的初始条件识别稳定吸引子 [14], [15]。然而，当只需要稳态周期响应时，长瞬态衰减和在大参数网格上的重复积分可能会变得昂贵。此外，不稳定分支和共存解通常需要打靶法 [16]、配点法 [17] 或延拓法 [18]。频域方法则直接近似周期响应，因此特别适合稳态响应、频谱成分和分支分析 [19], [20]。

在频域方法中，谐波平衡（HB）方法使用截断傅里叶级数表示周期响应，并将控制微分方程转化为非线性代数方程。直接系数平衡对低阶多项式非线性效率很高，但当涉及许多谐波、多个频率或非多项式本构关系时变得繁琐 [21]。因此，开发了配点 [22]、重构和频 - 时转换技术 [23]，以在不显式展开每个频域乘积的情况下评估非线性项 [24]。这些发展将 HB 分析扩展到了约束和多体机械系统 [7]，尽管解决由此产生的非线性代数方程仍然需要高效且可靠的雅可比矩阵。

增量谐波平衡（IHB）方法为 HB 方程提供了系统的线性化框架。Lau 等人通过在谐波投影前应用牛顿型线性化引入了增量过程 [25], [26]。该公式随后被扩展到受迫 [27]、自治和多时间尺度系统 [28], [29]。交替频域/时域方法在时域评估非线性力并将其转换回频域 [30]。尽管 IHB 和 HB-AFT 在计算顺序上有所不同，但它们的离散残差方程可以用密切相关的形式表示 [31]。

IHB 及相关 HB 方法的最近发展可分为三个主要方向。第一个涉及残差评估和代数效率。基于 FFT 的残差转换、Broyden 更新、张量收缩和 Galerkin 平均公式减少了重复数值积分和谐波组装 [32], [33], [34]。降维和 Krylov 回收技术进一步改善了大规模频域系统的求解 [35]。然而，基于 FFT 的公式对于完整且连续的傅里叶空间最有效。当不同坐标采用不同的谐波集时，当保留的频率集中在分离的频段中时，或者当傅里叶函数与其他基函数结合时，转换矩阵提供了更大的灵活性。

第二个方向涉及谐波截断和基构造。自适应谐波选择策略减少了对固定经验截断阶数的依赖 [36]，而基于重构的方法避免了复杂非线性项的显式频域展开 [37]。Zhou 等人用周期非光滑伯努利函数丰富了传统的傅里叶近似 [38]。通过表示切换事件处局部正则性的损失，这种富集可以加速收敛，而不需要长序列的高阶傅里叶谐波。

对灵活近似空间的需求得到了实验观察的支持。在频率上转换能量收集器中，低频主振荡器通过冲击、拨动或内部共振激发高频压电振荡器。因此，测量响应包含分离的低频和高频成分，主坐标和次坐标的主导谐波阶数可能存在很大差异 [39], [40]。类似的多尺度能量传递已在内部共振和多稳态收集器中被实验观察到 [41], [42]。因此，为每个坐标分配相同的连续傅里叶集可能会引入许多对实际响应贡献甚微的系数。依赖坐标且不连续的谐波集可以更经济地表示这些响应。

摩擦、间隙、碰摩和单向接触引入了不同的谱困难。实验摩擦接头结构表现出振幅相关共振、超谐波成分和失真的磁滞回线 [43], [44]，而转子 - 定子碰摩实验显示与接触状态变化相关的高阶和分数频率成分 [45]。尽管周期响应仍然具有离散谱，但非光滑切换通常导致其傅里叶系数衰减更慢，比全局光滑响应的那些系数慢。因此，纯傅里叶近似可能需要许多连续的高阶谐波。

涉及涡轮叶片接头、非光滑控制结构和颗粒 - 梁系统的实验进一步表明，正则性的损失通常局部化在有限数量的接触或切换事件附近 [46], [47], [48]。全局傅里叶基可以近似这种行为，但其收敛可能很慢，因为光滑性的局部变化必须在完整周期内表示。这些观察为混合近似空间提供了物理动机：傅里叶函数描述响应的全局光滑部分，而事件局部化的非光滑函数表示局部正则性损失。对非光滑近零刚度（QZS）隔振器的实验同样将接触约束与分段恢复特性和振动传递率变化联系起来 [48], [49]。

第三个方向涉及延拓和初始值构造。已开发了双时间尺度 IHB 公式用于耦合和准周期响应 [50], [51]，而延拓技术改善了折叠和稳定性边界的跟踪 [52]。Li 等人引入了一种数据驱动框架，用于预测强非线性 IHB 计算的可行初始值 [53]。他们随后结合了 FFT 残差评估、Broyden 雅可比近似和改进的 Levenberg-Marquardt 迭代，以提高包含非多项式非线性系统的鲁棒性 [54]。这些发展解决了重要的收敛困难，但主要是围绕单个系统或顺序连接的响应分支组织的，而不是同时管理大量独立的参数子系统。

两项最近的研究与目前的工作特别相关。Zhang 等人结合了龙格 - 库塔方法和 IHB 来跟踪非线性系统响应 [55]。首先在一个参数点获得 IHB 解，之后使用龙格 - 库塔型演化方程来跟踪响应曲线。这种策略减少了传统增量跟踪对不合适的牛顿初始值的敏感性，并提供了一种跟踪连接解分支的有效手段。其目的不同于目前的参数并行计算，在目前的计算中子系统保持相互独立，并且可能分别收敛或终止。

Li 等人开发了一种高效且鲁棒的 IHB 方法，用于移动基座倒立摆系统的分岔分析 [56]。他们的方法改善了分支跟踪和奇异或病态状态的处理，说明了 IHB 在复杂机械分岔问题方面的持续发展。目前的工作并不试图取代这种一般的响应延拓过程。在本文考虑的 Mathieu 示例中，延拓仅用于在预定稳定性边界上构建参数点；这些点的周期子系统随后使用所提出的并行框架独立组装和求解。

雅可比矩阵仍然是牛顿型 HB 和 IHB 公式的核心。解析微分可以保持模型结构和稀疏性，但对于强耦合、非多项式或分段本构关系，手动推导变得繁琐 [57], [58], [59]。有限差分避免了解析推导 [60]，但需要重复的残差评估，并引入了一个扰动步骤，其合适的大小取决于系数尺度和浮点精度 [61]。拟牛顿近似减少了完整雅可比评估的次数，但在病态状态附近可能会失去精度。随着谐波、坐标或参数点数量的增加，这些限制变得越来越重要。

自动微分（AD）为雅可比构建提供了另一条途径。它将链式法则应用于残差程序执行的基本操作，因此避免了手工闭式微分和有限差分扰动 [62]。AD 已被引入有限元残差、本构计算和非线性动力分析 [63], [64]。它也被纳入用于机械复杂非线性残差的 HB 公式 [65]。然而，AD 本身并不被声称是本研究的主要创新点。现有的 AD-HB 研究主要为串行系统或传统延拓计算构建雅可比矩阵；它们没有将端到端残差微分与参数批处理重构和子系统级收敛管理结合起来。

并行计算提供了加速非线性动力分析的额外手段。已为大规模流体、结构和优化问题开发了基于 CPU 和 GPU 的求解器 [66], [67]。也已为大型计算模型开发了并行频域求解器。例如，域分解已与谐波平衡方法结合，用于并行化具有摩擦接触的叶盘有限元模型的非线性受迫响应分析 [68]，而负载均衡和并行输入/输出已被用于提高谐波平衡流求解器的可扩展性 [69]。这些研究主要分布单个大规模模型的计算操作，而目前的工作考虑了许多相互独立的非线性子系统之间的参数级并行性。

参数级并行性引入了额外的管理问题。关于次级共振、非线性能量汇和强非线性频率响应的实验表明，不同的稳定响应可能在相同或邻近的参数条件下共存 [70], [71]。因此，相应的参数子系统可能表现出不同的收敛速度，收敛到不同的分支，或因初始值不合适而失败。有效的批处理架构不仅必须并行评估残差和雅可比矩阵，还必须保留每个子系统的状态、残差历史和终止状态。收敛或失败的点应从活动批处理中移除，而不中断剩余的计算。

尽管有这些进展，仍然存在三个相关的困难。首先，完整频域残差映射的雅可比矩阵仍然经常是手动推导、有限差分近似或为每个参数点单独构建。其次，传统的批处理计算通常对所有参数子系统施加共同的迭代历史，尽管不同的点可能在不同的迭代中收敛、发散或遇到不可微状态。第三，传统的连续傅里叶空间对于涉及依赖坐标的光谱和局部接触转换的多尺度和非光滑响应可能是低效的。

为了解决这些问题，本研究开发了一种鲁棒、公式无关、GPU 并行的增量谐波平衡框架，称为 IHB-ADP。其主要贡献总结如下：

•建立了一个端到端自动微分的残差映射。时域重构、非线性力评估、残差转换和频域雅可比构建保持在统一的计算图内，消除了手动推导的切线表达式和重复的有限差分扰动。

•独立的参数子系统被重组为批处理张量操作。子系统级收敛识别、发散退出和扰动重启允许收敛、失败和活动点被独立管理，而不中断完整的参数批处理。

•引入了一种专门的预测 - 偏差 - 校正过程来构建预定 Mathieu 稳定性边界参数集。该过程为随后的并行求解准备参数空间，并与非线性周期响应的顺序延拓不同。

•通过结合依赖坐标、不连续的傅里叶集和事件局部化的伯努利富集，将该框架扩展到了多尺度非光滑分层压电能量收集器。该示例检查了接触切换点附近的 AD，并展示了转换矩阵公式在混合基方面的灵活性。

所提出的框架使用五个示例进行了评估。范德波尔振荡器检查具有未知基频的自治响应。达芬振荡器演示了具有硬化非线性和共存解的参数并行频率响应分析。Mathieu 系统检查稳定性边界参数构建和并行周期响应求解。

使用 QZS 隔离平台的实验动态数据来评估模型 - 试验一致性，而堆叠式压电能量收集器则用于检查多尺度、坐标相关和非光滑响应。比较了 RK4、FD-IHB、HB-AD 和 IHB-ADP 的串行和 GPU 批处理实现，以区分由数值公式化引起的加速与参数级 GPU 并行提供的加速。本文其余部分安排如下。第 2 节介绍了 IHB 公式化、基于 AD 的雅可比矩阵构建、参数批处理架构和子系管理机制。第 3 节展示了五个数值和实验示例。第 4 节讨论了计算效率、AD 实现、可扩展性和局限性。最后，第 5 节总结了主要发现。
章节片段
方法论
方法论部分首先在 2.1 节回顾了 IHB 方法的基本理论，然后在 2.2 节简要介绍了 AD 技术，最后在 2.3 节详细阐述了 IHB-ADP 方法在非线性动力学问题中的技术实现。
数值示例
本节采用提出的 IHB-ADP 方法研究五个代表性非线性动力学系统。范德波尔振子是一个自治系统，其基本频率必须与傅里叶系数一起确定。杜芬振子展示了在存在硬化非线性和共存解情况下的参数并行频响分析。马丢系统用于检查沿...周期响应的并行求解。
讨论
本节讨论了所提出方法的四个特别与其数值实现相关的方面。首先，检查了识别 - 退出机制对参数批处理迭代的影响。其次，阐明了马丢示例中使用的稳定性边界跟踪策略的目的和范围。接着分析了计算效率、内存消耗和 AD 模式依赖。最后，讨论了 IHB-ADP 的当前局限性和可能的...
结论
本研究开发了一种自动微分驱动的参数并行增量谐波平衡方法（称为 IHB-ADP），用于同时求解多个非线性周期子系。从假设的频域系数到频域残差的完整映射被封装为一个可微计算图，使得无需手动推导非线性切线项或选择...即可构建所需的雅可比矩阵。
CRediT 作者贡献声明
Xiaohui Li：写作 – 初稿，调查。Yi Chen：写作 – 审查与编辑。Zeyuan Chang：概念化。Shijun Wang：数据整理。Yifan Jiang：形式分析。Yuhong Jin：方法论。Haiming Yi：方法论。Qingye Meng：调查。Lei Hou：写作 – 审查与编辑。
利益冲突声明
作者声明，他们没有已知的竞争性经济利益或个人关系，这些利益或关系可能会影响本文报告的工作。
致谢
衷心感谢国家自然科学基金（项目编号：12422213, U244120491, 12372008, 12502004），以及中国国家关键研发计划（项目编号：2023YFE0125900），中国国家重大科技专项（编号：2017-IV-0008-0045-JX），军工研究所基础研究稳定支持计划第二期项目（项目编号：03020051），中国博士后科学基金博士后 fellowship 项目（C 级）（项目编号：Xiaohui Li|Yi Chen|Zeyuan Chang|Shijun Wang|Yifan Jiang|Yuhong Jin|Haiming Yi|Qingye Meng|Lei Hou
哈尔滨工业大学航天学院，哈尔滨，150001，中国