几何数值积分方法在暂态稳定性计算中的应用

时间：2024-07-28

汪芳宗，郭梦芳，宋墩文，杨学涛，张磊

(1.三峡大学电气与新能源学院，湖北宜昌443002；2.中国电力科学研究院，北京100040)

0 引言

传统的数值积分方法主要包括线性多步法(linear multistep formulae, LMF)[1-2]以及Runge-Kutta(RK)系列方法[2-3]。电力系统暂态仿真计算中应用最为广泛的隐式梯形积分，是RK系列中的一种单级、2阶方法[4]；而另一种应用较多的方法是2阶的向后差分方法(backward differentiation formulae, BDF)，此方法属于LMF类方法[2]。

关于传统的数值积分方法，研究人员很早就建立了有关数值方法的稳定性以及其计算精度分析的理论体系，并给出了明确、具体的结论[3]。从工程实际应用的角度看，研究人员在选择传统的数值积分方法时主要是关注数值积分方法的数值稳定性，其次就是计算精度[5-7]。事实上，隐式梯形积分方法之所以在电力系统暂态仿真计算中得到广泛的应用，主要依据是该方法为具有最小局部误差截断系数的A-稳定的数值计算方法，而2阶的BDF方法则是LMF系列方法中具有A-稳定性以及L-稳定性的最高阶方法。

与传统的数值积分方法不同，几何积分方法主要是关注动力系统的物理/几何属性。如果数值积分方法能精确地保留动力系统的一个或多个物理/几何属性，则称为几何数值积分方法。几何属性的具体示例大致包括动力系统的对称性以及反对称性、辛结构、相空间中的体积(phase-space volume)等，而物理属性主要涉及动力系统的能量守恒性，例如首次积分、Hamilton能量函数、动量守恒等，因此，保几何属性的数值积分方法通常称为保结构计算方法(structure-preserving methods)[8,9]，而保物理属性的数值积分方法通常称为保能量守恒方法，简称为保能量方法(energy-preserving methods)[10]。

几何积分方法的兴起与发展得益于我国学者冯康先生的卓越贡献。早在20世纪80年代中期，冯康先生率先提出并系统建立了Hamilton动力系统的辛几何算法，也就是保Hamilton动力系统所具有的辛结构的保辛方法[11-12]。由于Hamilton动力系统具有辛结构与Hamilton能量守恒2大特征，因此，在冯康先生提出辛几何方法以后，研究人员自然就会想到并开始研究保Hamilton能量的数值积分方法[13]。换言之，冯康先生不仅创始了保结构计算方法，对保能量积分方法的研究也有重要的启发和推动作用。迄今为止，辛几何方法已在诸多领域得到了广泛应用。除辛几何方法以外，冯康先生还率先提出了无源动力系统的保积方法(volume-preserving methods)[8,14]，亦称为无散度或零散度动力系统。由于无源动力系统的真解也就是相流在相空间中是保体积的，因此，保积方法也是一类重要的保结构数值方法[15-16]。

关于几何积分方法的应用，一种思路是用几何积分方法来替换传统的数值积分方法。需要说明的是，尽管几何积分方法是基于动力系统的物理/几何属性而发展起来的，但这并不意味着几何积分方法不能用于一些不具备物理/几何属性的一般动力系统。事实上，在辛几何方法发明以前，RK系列方法中的Gauss类方法很早就在科学及工程计算中得到了应用，而Gauss方法就是一类辛几何方法[17]。另外一种思路就是：对具有特定的物理/几何属性微分动力系统，数值积分方法应该也必须保持相应的特征。换言之，对具有特定的物理/几何属性的微分动力系统，最好或必须采用相应的几何积分方法。这一思路是冯康先生提出的一个重要的学术思想[8]。

基于冯康先生的学术思想，本文尝试将几何积分方法应用于电力系统暂态稳定性计算。在暂态稳定性计算中，单机—无穷大系统既是一个经典的Hamilton系统，也是一个可分的无源动力系统。基于此思想，本文首先将Störmer-Verlet方法(也称为Leapfrog方法)、显辛RKN(Runge-Kutta-Nyström)方法[8]、简化的Takahashi-Imada方法(以下简称STI方法)[18]、平均向量场方法(averaged vector field method, AVF)[10]等应用于单机—无穷大系统的暂态稳定性计算中，并将上述方法与传统的数值积分方法进行对比分析。对单机—无穷大系统，Leapfrog方法以及显辛RKN方法既是保辛的，也是保积的；STI方法则只是保积的方法；而AVF方法则是保Hamilton能量或首次积分的数值方法。

在上述基础上，本文尝试将几何积分方法应用于多机系统的暂态稳定性计算。在采用经典模型的前提下，多机系统暂态稳定性计算的数学模型是一个可分的无源动力系统，但不是一个Hamilton系统，也就不具备辛结构及Hamilton能量守恒性。本文将Leapfrog方法、显辛RKN方法以及STI方法应用于多机系统的暂态稳定性计算，并将其与传统的隐式梯形积分方法进行对比。针对可分无源动力系统，这3种数值方法均是保积算法。

1 几何积分方法

1.1 可分Hamilton系统与显辛方法

设有2n个自变量p=[p1,p2,…,pn]T，q=[q1,q2,…,qn]T其中pi=pi(t)，qi=qi(t)。给定一个连续可微函数H=H(p,q)，称其为Hamilton函数或能量，则由该函数所确定的Hamilton系统为

(1)

或写成它的简洁形式(亦称正则方程)：

(2)

(3)

其中，In为n维的单位阵，J通常称为标准辛矩阵。若H=H(p,q)=U(p)+V(q)，则相应的Hamilton系统称为可分Hamilton系统。在可分Hamilton系统中，有一类特殊的Hamilton系统，具体可表述为

(4)

相应的Hamilton函数为

(5)

(6)

在不知道具体的Hamilton函数H(p,q)的情况下，式(6)是判断动力系统式(1)是否为Hamilton系统的准则。

针对Hamilton系统，研究人员已构造出了多种多样的辛积分格式。对一般的Hamilton系统，不存在显式的辛RK方法[8]，但对可分Hamilton系统，研究人员构造并发现了几类显式辛格式。由于隐式方法涉及的计算量较大，因此，本文主要考虑显式辛方法。以下介绍的显式方法均是辛方法。

(1)Leapfrog方法，具体可表述为

(7)

对形如式(4)的可分Hamilton系统，上述方法是2阶、对称的辛方法，其共轭或伴随方法也是2阶的辛方法。

(2)显辛RKN方法，其相应的Butcher表述为

3+3603-362+31203+3603605-33243+3121+3243-23121123+2312

(8)

对形如式(4)的可分Hamilton系统，上述方法是3级4阶的辛方法，其共轭或伴随方法也是4阶的辛方法。此类方法是由秦孟兆先生率先提出的[19]。

1.2 可分无源动力系统与保积方法

令n个自变量x≜[x1,x2,…,xn]T，f(x)≜[f1(x),f2(x),…,fn(x)]T。若有

(9)

无源动力系统具有一个重要的几何属性：其相流在相空间中是保体积的，即保持Ω=dx1∧dx2∧…∧dxn不变，因此，对无源动力系统，相应的数值积分方法也应该是保积的。对一个数值积分方法xn→xn+1，若它所生成的映射满足

(10)

则该数值积分方法就是保积方法[14-15]。式中det (■)表示矩阵的行列式。

关于无源动力系统的保积方法，冯康等在文献[14]中已证明：对一般的无源动力系统，传统的数值方法例如LMF方法以及RK方法均不可能是保积的，对一般的无源动力系统，构造保积算法是比较困难的一件事，但是，对一类特殊的无源动力系统也就是可分无源动力系统，研究人员发现并构造了众多的保积方法，而且均是显式方法。若一个动力系统可表述为

(11)

则该系统是可分无源动力系统，并且可分无源动力系统满足下列方程：

(12)

关于可分无源动力系统的保积方法，本文主要介绍以下几种保积方法：

①对形如式(4)的动力系统，Leapfrog方法式(7)是2阶的保积方法[18]。

②对形如式(4)的动力系统，RKN方法式(8)是4阶的保积方法[20]。

③对形如式(4)的动力系统，STI方法是4阶的保积方法[18]，其具体计算格式可表述为

(13)

式中α=1/12。这里需要说明的是：对可分Hamilton系统(也是可分无源动力系统)，STI方法只是一种保积方法，但一般情况下它不是一种辛方法[18]。

1.3 平均向量场方法

针对Hamilton系统的能量守恒性，研究人员提出了保Hamilton能量的平均向量场方法[10]。以Hamilton正则式(2)为例，平均向量场方法可表述为

(14)

平均向量场方法是2阶的隐式方法，理论上它可以精确保持Hamilton系统的能量守恒性。相较于显式积分方法，平均向量场方法的计算量较大，并且，该方法只是保能量的数值方法，不是辛方法。

2 单机—无穷大系统暂态稳定性计算

2.1 数学模型及其几何/物理特性

图1 单机—无穷大系统Fig.1 Single machine-infinity system

如图1所示的单机—无穷大系统，在发电机采用经典模型的情况下，相应的暂态稳定性计算的数学模型可表述为

(15)

式中各变量的物理意义或定义可参见有关电力系统分析的教科书，其中：M=11.28 s；Pm=1.0 p.u.；Ev/x∑=1.384 p.u.。

上述单机—无穷大系统结构虽然简单，但具备了一些重要的特征：此系统满足式(6)，既是一个经典的可分Hamilton系统，也是一个可分无源动力系统。其Hamilton能量函数表示为

(16)

式中δs是系统的稳定平衡点。将δs代入式(16)，可知Hamilton能量函数是守恒的[21]。即

(17)

可知，上述Hamilton能量函数式(16)也是单机—无穷大系统的首次积分。

2.2 数值测试结果

由于单机—无穷大系统具备上述几何/物理属性，因此，本文首先将Leapfrog方法式(7)、显辛RKN方法式(8)、STI方法式(13)、AVF方法式(14)分别应用于该系统的暂态稳定性计算，并将这4种方法的仿真结果与传统的隐式梯形积分方法的仿真结果进行了对比分析。

图2是3种不同的2阶数值方法计算出的相空间轨迹图。数值试验中，步长设定为h=0.02 s，积分初始值设定为δ0=1.111 2 rad，ω0=0.012 7 p.u.，此时系统处于临界稳定状态。理论上，若单机—无穷大系统是暂态稳定的，则其相空间轨迹是一个封闭的曲线[21]。

为测试数值方法的稳定性，对上述单机—无穷大系统进行长时间数值积分计算。得到的相空间轨迹如图2所示。如图2(a)所示，是隐式梯形积分法的相空间轨迹曲线，在经过157.4 s的长时间积分计算后，系统相空间轨迹发生破损(相空间轨迹不再是封闭的)，这是由于长时间仿真计算导致了错误的结果。如图2(b)所示，是利用AVF方法所计算出的相空间轨迹曲线，在经过14 172.8 s的长时间积分计算后，系统相空间轨迹发生破损，系统从稳定状态变为不稳定状态。如图2(c)所示，是利用Leapfrog方法所计算出的相空间轨迹曲线，在经过14 200 s的长时间积分计算后，系统仍然能够保持稳定状态。因此，从长时间计算数值稳定性的角度可知，尽管Leapfrog方法是显式方法，但该方法在保结构稳定性方面效果最好，AVF在保结构稳定性方面方法效果次之，隐式梯形积分方法在保结构稳定性方面效果最差。

(a) 隐式梯形积分方法

(b) AVF方法

图3 (2阶方法)Hamilton函数差值曲线 Fig.3 (2nd order method) Hamilton function difference curve

图3是3种不同方法计算出的Hamilton函数差值(用ΔH(t)表示)曲线，其中，ΔH(t)≜H(δ(t),ω(t))-H0，H0≜(δ0,ω0)。Hamilton函数差值曲线反映数值积分方法的保能量特性。如图3所示，根据能量差值的波动幅度可知，在保能量守恒方面，AVF方法的效果最好，Leapfrog方法的效果次之，隐式梯形积分方法的效果最差。

图4是4阶的RKN方法、STI方法以及传统的4阶显式RK方法所计算出的相空间轨迹图。其中，步长设定为h=0.05 s。测试结果显示，经过3 000 s的长时间积分计算后，提出的3种方法仍能保持相空间轨迹的封闭性。如图4(c)所示，传统显式RK方法的相空间轨迹在长时间积分后产生了“漂移”现象，可知其在保结构稳定性方面比保积的RKN方法以及STI方法效果要差。

图5是3种不同方法计算出的Hamilton函数差值曲线。图5(b)是RKN方法以及传统显式RK方法的Hamilton函数差值曲线。由图5可知，根据能量差值的波动幅度可知，在保能量守恒方面，RKN方法效果最好，传统RK方法效果次之，STI方法的能量差值变化幅度最大、其效果最差。

(a) RKN方法

(b) STI方法

(a) 3种方法对比结果

与传统的概念有所不同，通过对单机—无穷大系统的数值仿真分析，发现并验证了可分Hamilton系统或可分无源动力系统，显式保结构计算方法具备良好数值稳定性。

3 基于保积方法的多机系统暂态稳定性计算

基于显式保结构计算方法的优势，本节将显式保积方法应用于多机系统暂态稳定性数值计算，并将其与常用的隐式梯形积分方法进行对比。

3.1 经典模型及其几何属性

在负荷采用恒阻抗模型、发电机采用经典模型的前提下，基于收缩导纳矩阵的电力系统暂态稳定性计算的经典模型为

(18)

式中：m表示发电机的数量；Pei为电磁功率，其表达式可以写为

(19)

与单机—无穷大系统不同，上述多机系统的经典模型不是一个Hamilton系统，经验证可知，式(18)并不满足式(6)。多机系统的首次积分满足暂态稳定性分析中常用的暂态能量函数[21]，但不是一个守恒量。运用辛几何方法求解多机系统，虽然不具备“保结构”和“保能量”的特征，但上述经典模型仍具有重要的几何属性，即该系统是一个经典的可分无源动力系统。

3.2 数值测试结果

暂态稳定性计算的经典模型是一个可分的无源动力系统，本节将Leapfrog方法、RKN方法以及STI方法分别应用于多机系统暂态稳定性计算，并将计算出的结果与传统的积分方法进行对比分析。为此，利用隐式梯形方法采用小步长h=0.001 s进行数值仿真计算的结果作为基准值，对比所使用算法进行数值仿真计算的结果与基准值的差值Δδ(t)。

选用含50台发电机的IEEE 145节点系统[22]作为算例。暂态稳定性计算中，设定7号母线处发生三相短路，经0.1 s后切除故障。这种情况下，系统是暂态稳定的，其中，20号发电机受扰动影响最大，因此，将该台发电机的功角计算结果用于分析、比较不同的数值计算方法。其中隐式梯形方法设置小步长h=0.001 s。

图6是分别采用步长为h=0.02 s以及h=0.03 s时Leapfrog方法的差值曲线，图7是采用步长为h=0.05 s时RKN方法、STI方法以及传统4阶显式RK方法的差值曲线，图8是采用步长为h=0.06 s时3种方法的差值曲线。

图6 2阶Leapfrog方法差值曲线Fig.6 Difference curve of the second-order Leapfrog method

图7 4阶方法差值曲线(h=0.05 s)Fig.7 Difference curve of 4th order method (h=0.05 s)

图8 4阶方法差值曲线(h=0.06 s) Fig.8 Difference curve of 4th order method (h=0.06 s)

从图6可以看出，2阶的Leapfrog方法在采用较大步长(h=0.03 s)的情况下，仍然可以获得足够精确的、稳定的数值结果。从图7、8可以看出，在采用大步长的情况下，保积的RKN方法在计算精度上效果最好。

表1是关于计算效率的测试结果。根据表1中的计算结果可知，在采用同样步长的情况下，Leapfrog方法比隐式梯形积分方法约快3.3倍；在采用步长为隐式梯形积分方法的2倍情况下，RKN方法与STI方法大致相当，但两者均比传统显式RK方法略快，比隐式梯形积分方法约快5.8倍。

表1 计算效率测试结果Tab.1 Calculation efficiency test results