重构谐波平衡法及其求解复杂非线性问题应用1)

代洪华 王其偲 严子朴 岳晓奎

(西北工业大学航天学院,西安 710072)

(航天飞行动力学技术国家级重点实验室,西安 710072)

引言

工程中多数动力学问题,其数学模型都是非线性的,线性系统只是在一定假设及限制条件下对非线性系统的理想化近似[1].随着控制科学与航空宇航任务日益复杂,强非线性的影响越发不容忽视.谐波平衡法(harmonic balance method,HB)在求解非线性系统周期解中应用广泛,但作为一种半解析半数值方法,随着系统自由度、方法阶数的增加,公式推导工作将变得困难[2].借助计算机数学软件可辅助推导工作,但对计算效率及性能的提升作用仍然有限.采用7 阶HB 法推导立方非线性项,Mathematica代码长达1321 行.Hall 等[3]提出了高维谐波平衡法(high-dimensional harmonic balance,HDHB),通过“频域非线性量的时域近似计算”来简化公式推导,但是,由于引入近似关系导致非物理假解问题[4-5].Dai 等[6]发现了频域和时域非线性项之间的条件等价恒等式,基于此,首创了重构谐波平衡法(reconstruction harmonic balance,RHB)实现超高阶(N＞ 100)高精计算,并给出时域配点计算的最优采样定理,从理论上消除非物理解.

HB 类方法(HB 法及其改进方法)不仅用于计算简单非线性系统(杜芬、范德波方程等),还发展到航空航天、深空探测等前沿领域.哈尔滨工业大学陈毅等[7]使用谐波平衡−交变频域/时域法(HBalternating frequency-time,HB-AFT)用于求解航空发动机中双转子−轴承−机匣系统动力学方程.中山大学陈衍茂等[8]使用增量谐波平衡法对带机外挂载的二元机翼进行动力学特性研究.RHB 法也被作为三体轨道的设计依据,计算结果符合鹊桥中继卫星实际飞行数据[6].

但是该方法仍面临着两个问题: (1)谐波平衡法的本质依赖于非线性项的傅里叶级数展开,因此,受限于多项式型非线性系统求解,对于非多项式型复杂非线性问题,谐波平衡法难以适用;(2)已有谐波平衡类方法建立在单基频的假设上进行级数展开,由于拟周期响应存在多基频的特征,因此已有方法难以直接求解高精度拟周期解.

针对非多项式型非线性系统的求解问题,目前有两类处理方法,分别为直接法和间接法.直接法包括HB-Taylor 法[9]与HB-AFT 法[10-11].HB-Taylor 法通过泰勒级数展开将非线性函数用有限阶近似多项式描述;HB-AFT 法通过对非线性项的时域值采样以离散原问题.由于直接法对原系统进行了近似,导致求解精度低,且计算性能分别受制于高阶级数描述和过采样等问题.Cochelin 等[12]提出的谐波平衡−重铸法(HB-recast)是一种间接方法,通过重铸(recast)技术,成功将复杂非线性微分动力学系统无损变换为多项式型微分代数方程,然后用HB 法加以求解[13-14].但是,受限于HB 法的高阶计算(使用重铸法会增加系统的维度,进一步增加了高阶计算的难度)与原重铸形式的二次型限制(原方法要求新系统中非线性至多为二次多项式),至今难以对复杂非线性系统周期响应进行高效高精求解.

第2 个难题是拟周期响应求解问题.拟周期响应大量出现在非线性动力学系统中[15-17],其频率响应由多个不可约基频及其线性组合描述[18].由于传统HB 类方法基于单个基频进行近似解逼近,不能简单通过基频及其整数倍频率分量对拟周期响应描述,Chua 等[19]提出了结合广义傅里叶级数的改进多频HB 法,实现拟周期响应的准确求解.然而,使用多频HB 类方法对强非线性项进行频域内高阶描述时,计算效率严重受限于高维傅里叶分析(频域分量由多重积分[20]、求和[21]计算得到).Liu 等[22]使用多频HDHB 法求解受迫范德波振子的稳态响应,避免非线性项谐波系数的直接表示以提高求解效率,但是多频HDHB 法的精度受到非物理解的破坏[23].总之,当前关于HB 类方法的研究已拓宽到对拟周期响应求解领域,然而由于多基频计算中的高阶频域描述困难,对此类复杂响应的高性能求解尚待解决.

针对上述问题,本文提出了RHB 法与重铸法、多频谐波平衡计算相结合的两种方法: (1) RHB-重铸法;(2)重构多谐波平衡法(reconstruction multiple harmonic balance,RMHB).一方面RHB-重铸法通过将一般非线性等价转化为多项式型非线性系统,再采用RHB 法以确定时域最优采样点数,实现对复杂非线性动力学系统的高阶高精求解,仿真误差达计算机精度.另一方面,RMHB 法通过筛选和补充多个基频,借鉴RHB 法的时域等价重构思想,完善了学术界在使用多频HB 类方法求解拟周期响应时消除混淆假解的理论研究.

1 谐波平衡法及重构谐波平衡法

1.1 谐波平衡法

对非线性动力学系统

HB 法用有限阶傅里叶级数来构造近似解及其导数

其中N是HB 法的阶数.x0,x1,···,x2n为未知傅里叶系数,也被称为频域变量,将频域变量组成待求解向量=[x0x1···x2n]T.假设多项式函数f(x)=xϕ,忽略由非线性项而出现的高次谐波,只需展开前N次的傅里叶分量

其中各次分量r0,r1,···,r2n为

其中n=1,2,···,N;θ=ωt.将表示非线性函数的傅里叶分量记为,是待求解向量xˆ 的非线性函数.

将式(2)～式(4)代入微分方程(1),令常数项及前n次谐波 cosnωt和 sinnωt(n=1,2,···,N) 系数配平,得到非线性代数方程组

其中,分块矩阵A为

由三角函数的求导关系得到.

非线性项的近似表示(4)中,符号运算的复杂度会随算法阶次的提高呈指数增长[6].当HB 法的阶数超过20,即便有计算机数学软件的辅助,非线性项的推导整理工作量仍难以接受.

1.2 高维谐波平衡法

为简化HB 法的符号运算量,Hall 等[3]使用时域值替代频域分量,即将N阶HB 法中的傅里叶系数作用转换矩阵EHDHB,建立与一个周期 2N+1 个等距配点上时域量间的联系,定义

将式(7)代入HB 法代数方程组(6)得到HDHB法求解微分方程(1)对应的代数方程组

HDHB 法显著提高计算效率,并被认为是HB法的一种改进,在计算流体力学领域得到应用[2,24].但该算法求解强非线性动力学问题时产生非物理解(也称假解),严重影响求解精度.如图1 所示,HDHB 法计算结果虽然是方程组(8)的数学解(使求解算法收敛),但已偏离了真实的物理响应.

图1 HB 法和HDHB 法计算杜芬方程对比(N=2)Fig.1 Comparison of the Duffing equation computed by the HB and HDHB methods (N=2)

Liu 等[4]指出,假解现象产生于对非线性项近似的过程中,存在高阶谐波对低次的混淆(污染).以立方项非线性x3为例,二阶HB 法中() 表达式

HDHB 法计算中,非线性项的频域分量(9)中不仅包含了原HB 法中的所有项,还包括附加杂项,这些杂项的混入导致求解中出现非物理解

1.3 重构谐波平衡法

Dai 等[25]证明,HDHB 法与时域配点法等价,时域配点法通过使用时域配点上函数值对微分方程进行离散,从而联立代数方程组

其中配点矩阵

混淆是造成计算出现非物理解的原因.至于高次谐波在时域配点计算中影响低次谐波的机理,则遵循“混淆规则”[25].

混淆规则:假设 α ∈[0,π] 被以h为间隔均分为离散时间点,那么配点法中可区分的最大谐波次数n∈[−L,L],其中L=π/h,L为“极限波次”.高次谐波n(|n|＞L) 被误认为是相应低次谐波na

其中na∈[−L,L],m是整数.

混淆规则指出,时域配点法可区分的最大谐波次数由配点数(离散时间点的间隔)决定.增加时域配点法中的配点数量,可区分的谐波次数越高,此时方程(10)个数多于待求解变量数,配点矩阵E为列满秩矩阵,存在Moore-Penrose 广义逆矩阵E∗

其中I2N+1为单位矩阵.在配点法方程同时作用矩阵E∗,对方程组降维,得到RHB 法代数方程组

RHB 法在保证算法效率的同时,消除非物理解,从而实现对原HB 法的最佳重构.定理1 围绕如何选择合适的配点数,给出消除混淆确定条件[6].

定理1(条件等价性):如果配点数M,方法阶数N,多项式非线性的幂次 ϕ 满足

则RHB 法与HB 法等价.

为说明RHB 法消除混淆的效果.分别使用HDHB法和RHB 法计算杜芬方程,在分岔处(频率 ω=2)进行蒙特卡洛法模拟,选取104组随机初值(各频域分量在区间[−2,2]中选取),得到如表1 所示的统计结果.统计结果表明,HDHB 法计算产生58 组解,其中55 组都是非物理解;而RHB 法只计算得到3 个具有物理含义的解[6].

表1 RHB 法与HDHB 法在分岔处的解分布Table 1 Statistical distribution of solution by the RHB and the HDHB method at bifurcation point

此外,HB-AFT 法与HDHB 法计算流程一致,但HB-AFT 法的原理是,选用配点数 2ϕN+1 来消除混淆[11],但过采样会占用计算机资源,在实际计算时会导致更大的CPU 和RAM 计算负担[2].RHB 法基于时域配点法的统一框架,根据配点数的差异将所有HB 类方法(HDHB 法、HB-AFT 法等)建立起联系.

2 重构谐波平衡法改进策略

2.1 微分方程重铸技术

针对非多项式型非线性系统的HB 法求解,Cochelin等[12]提出将原系统改写为二次型系统

其中未知向量z包含微分方程的原始变量x及一些新的变量(引入的新变量用以改写方程).c是关于未知量z的常数向量;l是关于z的线性向量值运算符;q则是二次向量值运算符.

因为对任意xϕ次非线性,RHB 法都能实现高效计算,改写后方程可以是更高次多项式,将式(13)进一步写成适配于RHB 法的微分方程重铸形式

其中n可以是关于变量z的任意 ϕ 次多项式函数.重铸格式(14)涵盖多类型非多项式型非线性问题,下面将分类加以介绍.

(1)微分项耦合型

通过将方程重铸为标准形式(14),得到

原方程转化为典型的非线性度 ϕ=3 的动力学系统.得出结论: 导数项的耦合不影响非线性度计入,在实际HB 计算中,任意阶导数只计入一个非线性度.

(2)有理分式型

以Rayleigh-Plesset 方程为例

式中,A,B,C,D均为常系数.将方程两端同除以R,并运用倒数关系,令v=,x=1/R得到[6,12]

处理有理分式型非线性项,通过额外增加方程Rx=1,将方程改写为非线性度 ϕ=4 的多项式型非线性系统.

(3)根式型

相对论谐振子方程(17)中非线性项为根式

对任意形如 αq/p的根式,令 βp=α,使原根式型非线性转化得到多项式型: αq/p=(βp)q/p=βq.

(4)初等超越函数

对初等超越函数的处理,需要导数信息来实现对原微分方程的重铸.以非线性单摆方程为例

由于三角函数是超越函数,不能通过简单的代数关系式改写转化为多项式型.首先额外引入两个变量s,c来分别表示s(t)=sinθ(t),c(t)=cosθ(t).利用三角函数的导数性质

建立补充方程,从而实现对微分方程组的重铸[13]

初等超越函数非线性的重铸,以导数关系作为方程组改写的依据(部分需要用到有理分式、根式型非线性重铸技巧).附录表格中罗列了几类常见初等超越函数改写思路与重铸形式,可供参考.

总之,本文主要针对非线性项为连续函数情形(包括微分项耦合、有理分式型、根式型以及初等超越函数4 类),重铸技术通过变量替换、利用特殊函数的导数关系(见附表1)等关键步骤将一般连续函数非线性项改写为更利于重构谐波平衡法求解的多项式形式(14).

2.2 多谐波平衡计算

“补频”(supplemental frequency)思想[17]是在原来单频假设的基础上,补充多个与响应相关的频率.假设考虑两个基频,将假设函数改写为

参数m和n满足不等式

其中p代表二维傅里叶级数的截断[19-20],类似于RHB法中阶数N.

以RHB 法理论为基础,提出的RMHB 法利用时域信息,将动力学问题(1)转化为非线性代数方程(11)求解.多基频计算中,配点矩阵E以及转换矩阵E∗形式,记cm,n=cos(mω1+nω2)t,sm,n=sin(mω1+nω2)t.在RMHB 法计算中,存在选取合适的采样周期和时域配点数用以消除混淆的结论[23].

定理2(多频计算中条件等价性):假设多项式非线性的幂次 ϕ 的非线性系统响应中包含两个基频,基频之比 ω1/ω2为有理数.则RMHB 法消除混淆需满足采样时间T=2π/GCD(ω1,ω2),配点数

其中 GCD 为两数最大公因数.

3 数值仿真结果

3.1 相对论谐振子

作为物理学中经典问题,有必要对谐振子模型进行完整而严格的相对论处理[26].考虑一个静质量为m的质点在一维谐振力F=−的作用下做相对论运动.其中k为弹性常数,为位移量.结合牛顿运动运动学方程以及动量定理,可以推导得到相对论振子方程[27-28]

图2 Mickens 变换解与真实物理解对比Fig.2 Comparison of Mickens transformation result and real physical solution

RHB-重铸法按照根式型非线性重铸原则,将相对论谐振子改写为方程(18),再使用高阶RHB 法求解.结合表2 和图3,RHB-重铸法可对高速运动的谐振子(β=0.85)直接进行高阶计算.55 阶RHB-重铸法能将计算误差控制在 10−12量级(以数值积分为参照,全文采用MATLAB 内置函数ode15 s,相对精度容差 10−13,绝对误差容限10−20),总计算耗时在1 s 内.

表2 各阶RHB-重铸法求解相对论谐振子Table 2 Solving relativistic harmonic oscillator by the RHBrecast method with different orders

图3 55 阶RHB-重铸法求解相对论谐振子的计算结果Fig.3 The computing result of the RHB-recast method (N=55) for solving relativistic harmonic oscillator

除通过方法阶数对求解精度的提高,还能通过两种方式实现对计算性能的改善.

(1) 非线性代数方程组求解算法

非线性代数方程(11)求解算法能提高HB 类方法计算性能,Zheng 等[31]结合Tikhonov 正则化求解,黄建亮等[32]引入狗腿法思想结合回溯线搜索算法求解,Thomas 等[33]使用Broyden 法以提高HB 法计算性能.本文根据算例来说明L-M (Levenberg-Marquardt)法较之传统迭代方法在求解RHB 方程组时体现出的优势.图4 展示当 β=0.85,频域量初值估计x2=1,u0=0.7,使用两种不同的方程组求解算法: 牛顿迭代 (Newton-Raphson)法与L-M 法得到的一个周期内误差曲线.不同于牛顿迭代法,L-M 法迭代格式为[34]

图4 非线性方程组求解算法对计算精度的影响(相对论谐振子)Fig.4 Influence of nonlinear equation algorithm on calculation accuracy (relativistic harmonic oscillator)

其中xk为当前迭代未知变量值,Fk为函数值,Jk为雅可比矩阵.L-M 法通过衡量每步迭代的误差是否发散,来决定撤回迭代并将参数 λk按十倍放缩.LM 法在求解非线性方程组时显示出更高的计算精度,同比牛顿迭代法,精度提高了 104以上.

(2) 合理选择代数方程

相对论谐振子(17)为保守系统,振动频率 ω 是状态变量[35].使用RHB-重构法对n自由度保守系统求解时,共需考虑n(2N+1)+1 个变量,需额外增加一个代数方程使RHB 方程组适定.本文以初值约束来探讨方程对计算性能的影响.图5 为采用55 阶RHB-重构法与L-M 法求解器,使用不同代数方程(约束条件)得到的误差曲线.分别对应: 初始位移约束x(0)=0;初始速度约束(0)=β;同时考虑初始位移与速度约束.同时考虑位移与速度约束时,RHB-重构法计算精度更高.

图5 改变代数方程对计算精度影响Fig.5 Influence of changed algebraic equation on calculation accuracy

相对论谐振子在高速运动时显示出复杂的动力学特性,需要高阶(＞20 阶)方法来进行分析.使用重铸技术,将根式非线性转化为多项式型,再使用RHB法实现快速解算,克服高阶估计的限制.如图6 所示,RHB-重铸法适用于不同初速度条件下相对论振子的计算.

3.2 非线性单摆

非线性单摆(19)是物理学中经典问题,且时域响应具有典型的周期性.但现有HB 类方法不能做到高性能求解: HB-Taylor 法需要采用高的截断阶来保证计算精度;HB-AFT 法则需进行过采样以抑制混淆误差(如图7 所示).

对非线性单摆重铸,得到方程(20).使用25 阶RHB-重铸法求解(大角度摆动,θ(0)=1.5),初值估计x1=1.423,s1=1.065,c0=1.028,辅助L-M 法求解器,计算结果如图8 所示,与数值积分参考解比较,误差控制在 10−12数量级以下.

图8 25 阶RHB-重铸法求解非线性单摆Fig.8 The computing result of the RHB-recast method (N=25) for solving nonlinear pendulum

分别采用牛顿迭代法、Tikhonov 正则法与L-M法得到图9 计算结果.牛顿迭代法中雅可比阵奇异,求解失效;Tikhonov 法与L-M 法计算所得的误差数量级相近,但是由于Tikhonov 法对正则参数的优化使求解更耗时[31].本例中,使用L-M 法计算RHB 法方程组仅需0.72 s,而Tikhonov 法耗时达1.17 s,L-M法较之Tikhonov 法效率提高了62.5%,达到了相近的求解精度.对比传统方法(牛顿迭代法)与优化方法(Tikhonov 法),L-M 法都更适合于RHB 法方程组的求解.

图9 非线性方程组求解算法对计算精度的影响(非线性单摆)Fig.9 Influence of nonlinear equation algorithm on calculation accuracy (nonlinear pendulum)

非线性单摆问题是保守系统,其周期解 θ(t) 与响应频率都依赖于初始振幅 θ(0).使用RHB 法求解重铸方程组(20)时,本文给出3 种代数方程组合方案,以助于对比分析.

方案1: 对未知量 θ,v,s和c全谐波平衡(从常数项到N次谐波项),计 4(2N+1) 多个方程.增加初始振幅约束 θ(0)=α.

方案2: 对未知量 θ 和v全谐波平衡,计2(2N+1)多个方程.s和c从1 次谐波开始,平衡系数到N次,计 2N多个方程.增加初始振幅约束 θ (0)=α 与两个对附加变量s,c的初值约束

方案3: 对未知量 θ 和v全谐波展开,计2(2N+1)多个方程.s和c从1 次谐波展开到N次,计 2N多个方程.增加初始速度约束v(0)=0 与两个对附加变量s,c的初值约束式(24).

分别采用3 种方案得到的计算结果如图10 所示.对比方案2 和方案3,方案1 的计算误差大,主要原因是没有对附加变量s和c的初值进行约束.方案3 同时利用了初始位移与速度信息,比方案2 计算精度更高.

图10 不同代数方程组合方案对计算精度影响Fig.10 Influence of different combination schemes of nonlinear algebraic equation on calculation accuracy

对比求解非线性单摆的3 种方法(HB-AFT法、RHB-Taylor 法和RHB-重铸法)计算结果.考虑同阶截断N=25,综合图11 误差曲线和表3 的计算结果,RHB-重铸法确定最优时域配点数M=76,降低采样成本.RHB-Taylor 法虽可采用高阶截断(15 次多项式)保证计算精度,但超越函数的级数表示降低计算效率,计算用时达1.80 s.RHB-重铸法与HB-AFT 法计算效率相当,但将计算精度提高了两个数量级以上.

表3 3 种方法计算非线性单摆结果对比Table 3 Comparison of corresponding results of three methods for solving nonlinear pendulum

图11 3 种方法求解非线性单摆对应误差曲线对比Fig.11 Comparison of corresponding error curves of three methods for solving nonlinear pendulum

RHB-重铸法适用于求解初等超越函数非线性问题.尤其是当单摆做大幅度振动时,传统的线性化手段无法很好地捕捉动力学响应,而RHB-重铸法则可以提供高精度解析解(图12).

图12 RHB-重铸法求解非线性单摆问题: 相平面图Fig.12 The RHB-recast method for solving nonlinear pendulum: phase plot

3.3 非线性耦合的非对称摆

用绳将质点小球悬挂于一根固定在铁架上的弹性杆末端,而弹性杆能在水平面与质点小球同步振动.此物理模型在许多其他二维摆系统中普遍存在,如球面摆、傅科摆等.特别地,当傅科摆由于不完美的悬挂或由于侧向运动引起非线性耦合而导致不对称,可能会导致额外的旋转,从而掩盖地面效应.非线性耦合的非对称单摆微分形式可以写作[36]

此方程中2 阶导数耦合入非线性项,导致使用传统方法进行直接求解变得棘手.以传统的有限差分类方法(欧拉法、龙格库塔法、MATLAB 内置ode算法等)而言,需额外计算代数方程组来辅助求解.而部分解析求解法的计算效果同样不佳,Jia 等[36]曾采用多时间尺度法推导方程的解为

令参数 κ=0.01,初始条件x(0)=0.1,y(0)=0.2,(0)=(0)=0.

以修正Chebyshev-Picard 迭代 (modified Chebyshev-Picard iteration,MCPI)法计算结果为标准解[37-38].仿真得到如图13 所示,解(26)不仅推导复杂,且与真实的物理解间有较大误差.

图13 多时间尺度法与修正Chebyshev-Picard 迭代法(参考解)求解轨迹对比Fig.13 Comparison of the method of multiple scales and the MCPI method (reference solution) for solving trajectory diagram

对多时间尺度解(26)分析得知,动力学响应包含两个基频(所有频率成分都是基频的线性组合),通过快速傅里叶变换(FFT)得到: ω1=0.985 7,ω2=0.993 5.又从微分方程重铸的角度分析非对称摆微分方程(25),因为高阶导数只算作一个非线性度,该方程是具有3 次多项式非线性的动力学系统.使用5 阶RMHB 法进行计算得到图14 所示仿真结果,计算用时8.24 s,x方向振动幅值误差 9.21×10−3,y方向振动幅值误差 3.90×10−7.

图14 5 阶RMHB 法求解非线性耦合非对称摆Fig.14 Solving nonlinear coupling asymmetric pendulum by the RMHB5

对2 阶导耦合型非线性系统,有限差分方法不能直接求解,一些解析求解方法的计算精度低.考虑两个基频的RMHB 法不仅保证计算效率,还能对拟周期运动进行准确的预测.

4 结论

围绕复杂非线性动力学系统高效高精度周期解算需求,本文基于微分方程重铸、“补频”等思想,分别提出RHB-重铸法与RMHB 法,主要的工作与结论总结如下.

(1)提出RHB-重铸法,将一般非线性问题改写为多项式型非线性方程,配合RHB 法确定消除非物理解所需的最优配点阈值,实现对非多项式型非线性动力学系统周期响应进行高阶预测.数值实验说明,对比HB-AFT 法的时域过采样和HB-Taylor 法对非线性函数进行泰勒级数截断两种处理方式,RHB-重铸法的计算精度高达 10−12,且计算时间不超过1 s,同时保证了求解的高精度与高效率.

(2)提出结合“补频”思想的RMHB 法,在原来单频假设的基础上,补充多个与物理响应相关的频率.借鉴RHB 法中时域重构等价性,给出多频谐波平衡计算的最优配点数结论,并充分利用非线性量的时域采样代替复杂的高维傅里叶分析,最终实现对拟周期响应的准确捕捉.

(3)通过解算相对论谐振子、非线性单摆以及非线性耦合的非对称摆问题验证了本文算法的有效性.针对传统有限差分类方法求解困难的多自由度耦合系统,RHB 及本文提出的两种方法都不受方程形式的限制.此外,在相同的系统参数设置下,合理地选择代数方程求解器、代数约束方案,可有效提高非线性问题的求解精度.

本文提出的RHB-重铸法与RMHB 法对复杂非线性动力学系统周期性及拟周期响应的计算效率与精度相较于现行的方法与处理方式都具有优势.在后续的研究中,将探究RHB 法在高精周期响应求解方面的工程应用.

数据可用性声明

支撑本研究的科学数据已在中国科学院科学数据银行 (science data bank) ScienceDB 平台公开发布,访问地址为: https://cstr.cn/31253.11.sciencedb.j00140.00022 或https://doi.org/10.57760/sciencedb.j00140.00022.

附录A

附表A1 常见初等超越函数的重铸Table A1 Recast form of the most common elementary transcendental functions