用于求解非线性分数阶微分-代数方程周期稳态解的数值方法

《Mathematics and Computers in Simulation》：Numerical method for obtaining periodic steady-state solutions of nonlinear fractional differential-algebraic equations

　　马辛·索瓦
西里西亚工业大学电气工程系，波兰西里西亚格利维采市阿卡德米卡2A街，邮编44-100

摘要
本文提出了一种新的数值方法，用于求解时变非线性分数阶问题中的周期性稳态解。该方法避免了传统中使用谐波平衡的思路。其独特之处在于通过选定解点处的函数值来重构稳态波形（类似时间步进法的做法）。这种方法能够更好地处理周期波形中的较大失真，同时避免在处理具有较长稳定时间的瞬态现象时出现过长的计算时间（如时间步进法所面临的问题）。所研究问题的通用形式为非线性分数阶微分-代数方程组，这一形式源于包含经典元件与非线性分数阶元件（如非线性分数阶线圈）的电路问题。本文的分析方法也可应用于其他可通过此类非线性分数阶微分-代数方程组描述的领域。文中详细阐述了分数阶导数的近似方法，以便感兴趣的研究人员能在自己的研究中加以应用。该近似方法已在存在解析解的周期波形示例上进行了验证。至于非线性分数阶微分-代数方程组的求解，所提出的方法也在含有分数阶元件和非正弦源的非线性电路示例上得到了测试，这类电路的解波形中含有显著的高次谐波成分。测试不仅关注精度，还考虑了计算时间。

引言
本研究源自分数阶微积分在电路理论中的应用，但其数学描述、数值方法本身以及计算流程也可应用于其他需要时变分析的领域。本文首先回顾了分数阶导数的卡普托定义[1]。

定义1
对于函数x(t)，其α∈[0,1]阶的卡普托分数阶导数定义为：
t0Dtαx(t)=1Γ(1?α)∫t0tx(1)(τ)(t?τ)αdτ，
其中t0为分析的初始时刻，Γ(?)为伽马函数。

由于需要求解周期性稳态解，因此还需做出某些假设。此时是在t→∞的条件下寻求解。另一种假设是令t0=?∞，即认为稳态状态已经持续了无限长时间。如[3]中所提及，有一种与魏尔导数[2]相关的有用公式。

定义2
对于α∈[0,1]，所采用的卡普托分数阶导数的变体定义为：
?∞Dtαx(t)=1Γ(1?α)∫?∞tx(1)(τ)(t?τ)αdτ。

未来还可以应用其他分数阶导数定义。通常这会大幅改变第4节中的核心内容，但不会影响整体研究思路以及后续阶段解的求解方式（第6节）。

尽管分数阶微积分已在诸多领域展现出应用价值，比如在建模[4][5]方面，或在应用数学中构建系统的一般模型[6][7]，但作者更关注其在电气工程领域的应用，尤其是电路理论领域。近年来，分数阶微积分在电气电路分析中的应用越来越普遍。它在超级电容器[8]、铁磁芯线圈[9]、燃料电池[10]，乃至高频交流激励器[11]和电气驱动装置[12]等复杂设备的建模中都显示出重要作用。

运用分数阶微积分的电路还可用于设计具有特定功能的电路，例如滤波器[13]。由于这类元件的实用性，人们多次将其作为电路理论中经典元件的扩展而引入，但这也引发了一些疑问、争议和讨论[14][15][16]。如[17]中所提到的，一种较为折中的观点认为，分数阶元件应被视为更为复杂的元件，仅在现象学层面（如建模时）有用，而非直接源自经典理论，也不具备经典电磁场理论中的直接物理解释。

分数阶线圈由以下分数阶微分方程描述：
Lαt0Dtαi(t)=u(t)，
其中Lα为一个人为设定的参数，其单位为H?sα?1。而分数阶电容器的描述则为：
Cβt0Dtβu(t)=i(t)，
这里的Cβ同样不属于经典电路理论中的参数，其单位为F?sβ?1。

在引入有助于建模的电路元件时，还可以引入非线性分数阶线圈（可用于建模铁磁芯线圈[18]），其方程为：
t0Dtαψ(t)=u(t)，
其中ψ是一个单位为Wb?sα?1的变量，它并不表示磁链，而只是一个中间变量，需结合非线性函数ψ(i)或i(ψ)一起使用——类似于经典非线性线圈中引入的磁链与电流之间的依赖关系。虽然ψ不代表磁链，但通过该模型仍可使用以下公式计算磁链Ψ：
Ψ(t)=Ψ(t0)+∫t0τu(τ)dτ。

如此定义的非线性分数阶元件依然描述了电压与电流之间的关系，只是方式更为复杂。非线性分数阶电容器（目前应用尚不十分普遍[19]）的方程为：
t0Dtβq(t)=i(t)，
这里的q并非电荷，而是一个单位为C?sβ?1的变量。与非线性分数阶电容器类似，此处也需要引入非线性函数q(u)或u(q)。

文献中提出了多种表示分数阶电路元件的符号，本文采用了最初见于[20]的符号，并稍作修改，如图1所示。修改之处在于，原本以较大字体显示在电路符号旁边的导数阶数（如α），现在改为以括号形式出现在元件名称的下标中（如L(α)ψ）。电路符号本身则与[20]中的保持一致。

本文结构如下：第2节阐述开展此项研究的具体动机；第3节详细介绍分数阶微分-代数方程的通用形式；第4节讲解新方法中分数阶导数的处理方式；第5节通过解析解对分数阶导数的单独计算进行验证；第6节说明如何将所研究的分数阶微分-代数方程转化为非线性方程组；第7节介绍一个复杂的电路示例，列出其方程，利用所提方法求出稳态解，并通过与经典时间步进法得到的解进行比较来验证结果，同时还会探讨计算时间问题；最后，第8节总结研究成果并提出未来研究方向。

本文的大部分结果都是通过GNU Octave[21]软件计算获得的，部分案例还使用了预编译的DLL动态链接库（具体将在后文提及）。本文中的计算脚本（稍作修改后）也应能在Matlab[22]中运行。

研究动机
本研究受到以下因素的启发：
• 需要能够获得参考解的方法，这些参考解在构建高级数值方法时非常有用（详见第2.1节）；
• 电气工程实践分析中观察到的稳态波形现象（包括非正弦波形）所引发的思考（这一动机实际上是本研究的初衷），相关内容见第2.2节。

在仿真过程中，还会遇到各种……

已求解分数阶微分-代数方程的通用形式
当问题的数学描述明确之后，就可以推导出需要求解的方程的通用形式。对这类方程进行一般性处理，有助于得到其通解，或者针对更复杂的问题，基于数值方法构建求解器，从而有效解决特定类型的问题。此前针对第1节中提到的分数阶电路元件（见图1）所做的分析[18][20][33]，已经促成了相应通用形式的推导。

分数阶导数的处理
为便于进一步分析，假设v(t)的每个分量都可以用周期为T的周期函数来描述（一般而言可为T/n，其中n∈N+）。同时假设每个待求未知量的周期波形，都是根据在等距节点t0、t1、t2……tN处确定的解来重构的（满足tN?t0=T）。所有待求未知量ξ都被视为满足周期性条件的有限连续函数，即ξ(t)=ξ(t?T)。

考虑到以上条件，解并非……

分数阶导数近似——测试示例
在解释如何将第4节中提出的分数阶导数近似方法应用于形式为(10)的分数阶微分-代数方程之前，我们先通过一个特定的周期函数示例来验证分数阶导数本身的近似效果。需要计算的是在t∈[0,nTT]区间内的稳态时间函数y(t)，它是给定函数χ(t)的分数阶导数：y(t)=?∞Dtαχ(t)，其中测试函数为χ(t)=tanpssin(ωt)，且要求ps<π/2以确保函数的连续性。选择这个示例是因为……

用于稳态解的非线性方程的构建
方程(10)在使用时间步进法时可以直接应用。而在本研究中，我们提出了一种方法，可以同时为所有nw个计算波形的全部n个时间节点求解（因此总共需要求解nwn个未知量）。关于如何安排这些未知量，有两种思路：第一种是将相邻节点上的元素相互关联；第二种则是按顺序排列所有的解向量，即w(t1)、w(t2)、……w(tn)。第一种思路会……

计算测试
第5节通过一个特定的稳态计算示例，展示了提高分数阶导数近似精度的趋势。在本节中，我们将对能够用形式(10)表示的问题，全面测试所提出的稳态计算方法。本研究设计的算法所得结果，将与经过充分验证的时间步进法的结果进行比较（相关内容简要介绍于第7.4节）。

结论
本文提出的新型数值方法为处理具有长稳定时间、非正弦周期源以及强非线性的非线性问题提供了替代时间步进法的解决方案，这类问题会导致解波形中高次谐波成分极为显著。需要再次强调的是，在开展本研究之前，我们已对相关文献进行了全面梳理，结果表明此前尚未有类似方法的出现。

致谢
本研究未获得公共部门、商业机构或非营利组织提供的任何专项资助。 2以确保函数的连续性。选择这个示例是因为…… 用于稳态解的非线性方程的构建 方程(10)在使用时间步进法时可以直接应用。而在本研究中，我们提出了一种方法，可以同时为所有nw个计算波形的全部n个时间节点求解（因此总共需要求解nwn个未知量）。关于如何安排这些未知量，有两种思路：第一种是将相邻节点上的元素相互关联；第二种则是按顺序排列所有的解向量，即w(t1)、w(t2)、……w(tn)。第一种思路会…… 计算测试 第5节通过一个特定的稳态计算示例，展示了提高分数阶导数近似精度的趋势。在本节中，我们将对能够用形式(10)表示的问题，全面测试所提出的稳态计算方法。本研究设计的算法所得结果，将与经过充分验证的时间步进法的结果进行比较（相关内容简要介绍于第7.4节）。 结论 本文提出的新型数值方法为处理具有长稳定时间、非正弦周期源以及强非线性的非线性问题提供了替代时间步进法的解决方案，这类问题会导致解波形中高次谐波成分极为显著。需要再次强调的是，在开展本研究之前，我们已对相关文献进行了全面梳理，结果表明此前尚未有类似方法的出现。 致谢>
用于稳态解的非线性方程的构建
方程(10)在使用时间步进法时可以直接应用。而在本研究中，我们提出了一种方法，可以同时为所有nw个计算波形的全部n个时间节点求解（因此总共需要求解nwn个未知量）。关于如何安排这些未知量，有两种思路：第一种是将相邻节点上的元素相互关联；第二种则是按顺序排列所有的解向量，即w(t1)、w(t2)、……w(tn)。第一种思路会……

计算测试
第5节通过一个特定的稳态计算示例，展示了提高分数阶导数近似精度的趋势。在本节中，我们将对能够用形式(10)表示的问题，全面测试所提出的稳态计算方法。本研究设计的算法所得结果，将与经过充分验证的时间步进法的结果进行比较（相关内容简要介绍于第7.4节）。

结论
本文提出的新型数值方法为处理具有长稳定时间、非正弦周期源以及强非线性的非线性问题提供了替代时间步进法的解决方案，这类问题会导致解波形中高次谐波成分极为显著。需要再次强调的是，在开展本研究之前，我们已对相关文献进行了全面梳理，结果表明此前尚未有类似方法的出现。

致谢
本研究未获得公共部门、商业机构或非营利组织提供的任何专项资助。>