电力系统电磁暂态仿真算法研究综述

杨明，张永明，张子骞，顾禹轩

(1. 上海电机学院 电气学院， 上海 201306； 2. 格拉茨技术大学 电力系统研究所，奥地利格拉茨)

0 引 言

目前大规模新能源的并网以及新型电力电子器件的应用，电力系统数字仿真复杂性不断提高[1-2]，电磁暂态特性也变得日趋复杂。在实时仿真系统中，使用数字计算的仿真必须无限接近物理实际系统，要求仿真时间与被仿真系统时间相对应并能在一个仿真步长内快速计算所有差分方程。实时数字仿真器(Real Time Digital Simulator，RTDS)是能将电力系统电磁暂态仿真实时进行分析的一种并行系统[3-4]。电磁暂态仿真分析根据各类模型的不同特性对电力系统中相关元器件进行精准建模，准确地获得电力系统的动态特征，在分析电力系统中不同暂态过程时，这就要求电磁暂态仿真具备高性能仿真能力，目前基于Dommel算法的电磁暂态程序EMTP(Electromagnetic Transients Program，EMTP)通常应用在电力系统电磁暂态仿真中，我国中国电力科学研究院研发的EMTPE和加拿大的PSCAD/EMTDC都是基于EMTP程序而开发。为克服未来电网电磁暂态实时仿真大规模的挑战，基于HYPERSIM的电磁暂态实时仿真技术的研究也有了新突破[5]。电力系统电磁暂态仿真程序软件一般采用定步长的数值积分算法用于描述电力系统中的基本元器件，同时对各元器件进行离散化处理，考虑到仿真计算过程中的稳定性与精度问题，梯形法成为各种数值分析方法的首选方法，但梯形法最大的问题是在数值求解过程中产生数值振荡现象。另外在电磁暂态仿真分析过程中也会因其它各种原因出现同样现象，使仿真过程不准确。数值积分代换技术是解决电磁暂态问题的普遍方法，主要特征在于其简单性、全面应用性以及计算高效性[6]。对于电力系统的刚性特性，算法会在求解稳定性上受到影响，非线性特性也会影响求解过程的计算效率。由于数值积分算法特性和性能的差异性，因此选择合适的算法是电磁暂态仿真分析的基础和核心。

近年来，一些新算法和新技术的研究使电磁暂态仿真在其稳定性、计算精度和仿真效率上得到提高，应用领域也更为广泛。文章讨论了目前应用在电磁暂态仿真中的经典积分算法和研究者提出的新型积分算法，并通过讨论算法的性能指标对比分析了所应用方法的优缺点。

1 仿真理论基础

1.1 理论基础

电力系统中各类数学模型和与之对应的数值算法的求解构成了电力系统电磁暂态仿真过程[7]，例如图1所示的电感支路。

图1 电感支路及诺顿等效电路

由图1可得其微分方程如式(1)所示：

(1)

对式(1)方程进行差分化，由此可得到差分方程[8-10]。例如对式(1)采用梯形积分法后可得：

(2)

式(2)差分方程可表示成如图1所示诺顿等效电路(伴随电路)，作用于当前时步电压项部分等效为电导，前一时步电气量并作用于当前时步电流项部分等效为电流源历史项，通过节点方程联立式(2)可得式(3)所示的电磁暂态仿真基本方程：

Gu=i

(3)

节点分析法在大规模电力系统中相较于状态方程法在计算速度上更占优势，例如传统电磁暂态仿真程序即采用了节点分析法。

1.2 问题与措施

随着电网规模的不断扩大，电力系统大量应用电力电子装置的趋势使其复杂性日趋增加，传统电磁暂态仿真性能已经无法满足需求[11]，对电磁暂态仿真分析的精准建模能力要求也不断提高。在仿真过程中，如开关动作过程中存在电感电容元件的换路过程，将会引起非状态变量的突变，若继续采用梯形积分法将会引起非原型的数值振荡[12-14]；同样，因开关动作不在整步长时刻点，定步长算法隐开关的延迟会引入非特征谐波，导致仿真失真[15]；此外控制系统与主系统之间产生的外部时延问题也将导致数值波动现象。由于控制系统中非线性环节产生的内部时延，也会导致在控制系统求解时出现数值不稳定问题[16]。

针对数值振荡产生的原因分析，研究者提出了许多解决方法：(1)求出突变后的非状态变量。但对于实 际复杂电力系统不切实际；(2)避免在网络结构突变时引入非状态变量。例如在采用后退欧拉法时，由于后退欧拉法在求解方程中不会出现t-Δt时刻的非状态变量，因此在突变过程中，不会产生不合理的等值注入电流源，从而避免了数值振荡问题，但这种解决方式精度不高[17]；(3)减小时间步长。由于数值振荡的幅值随1/Δt变化，因此较小的时间步长反而会造成数值振荡问题，并不普遍适用；(4)添加阻尼电阻。该方法中阻尼电阻Rd与步长Δt、阻尼因子α选取有关，但选取值不能保证最优，影响计算精度，并且当电阻值选取太小时，反而会加剧数值振荡[18]；(5)文献[19]通过适当的控制模型来最小化或补偿解决外部时延。另外文献[20]提出了一种基于滤波的抑制方法解决外部时延。为消除内部时延，文献[21]采用牛顿法来获得控制系统的联立求解。

围绕电力系统电磁暂态仿真需求，数值积分算法是开展电力系统非线性刚性系统的基础与核心，由于数值积分算法的不同特性，应考虑到算法的稳定性、仿真精度以及仿真效率。下一章则将着重分析所提及到的算法性能指标。

2 数值积分算法的性能指标

2.1 计算稳定性

(4)

对于线性微分动力方程，具备L稳定性的数值积分算法也能够消除数值振荡，基于对刚性系统的适应性以及仿真结果的准确性要求，选取电力系统仿真积分算法时，考虑计算稳定性是准确模拟电力系统电磁暂态仿真的前提。

2.2 计算精度

不同数值积分算法，相应的精度和误差也不尽相同。在采用数值方法求解微分方程时会存在截断误差，例如当yn=y(xn)时，用数值积分算法计算yn+1的误差：

Rn=y(xn+1)-y(xn)

(5)

式(5)中Rn即为局部截断误差。通过泰勒公式展开可得数值方法的局部截断误差为O(hP+1)，其中P为阶数。当步长h越小时，P值就越高，而局部截断误差与P值成反比，Rn将不断减小，因此计算精度就越高。显然若在同等步长h的数值积分算法中，阶数越高，计算精度则越高。

2.3 计算效率

数值积分算法的计算效率是满足电磁暂态仿真实时性要求的关键。显示积分法仿真效率较高，但稳定性差，很少应用于电磁暂态仿真中。文章主要选择隐式积分法进行计算效率的对比。可通过采用高斯消元法求解离散后的方程[23]，以求解后的浮点运算数作为计算效率的指标。

文章对数值积分算法的讨论中，对图2所示RL电路进行求解计算，以式(6)所示微分方程一般形式为基础，对运用在电力系统电磁暂态仿真的欧拉法、梯形法、均值法、CDA法、THTA法以及新型积分算法(块广义向后差分法、精细积分法、根匹配技术)进行讨论。其中x(t)为状态变量；f为时间t与x的函数。

(6)

图2 进行计算的电路

3 经典积分算法研究

3.1 欧拉法

3.1.1 前向欧拉法

前向欧拉法的积分格式如式(7)所示:

yn+1=yn+hf(tn,yn)

(7)

前向欧拉法采用固定步长，计算简单，但由于是显示积分，等效电阻为负，因此稳定性差。在选取为固定步长后，当仿真步数越多，则计算误差就越大，当步长较大时，前向欧拉法的精度并不高，因此在实际数值求解中不经常使用。前向欧拉法的稳定域如图3所示。

图3 前向欧拉法的稳定域

3.1.2 后向欧拉法

后向欧拉法采用固定步长，计算简单，计算效率高，积分格式为：

yn+1=yn+hf(tn,yn+1)

(8)

图4 后向欧拉法的稳定域

该方法在求解方程中不会出现t-Δt时刻的非状态变量，因此在突变过程中，不会产生不合理的等值注入电流源，从而不会引起数值振荡问题。后向欧拉法的等效电阻为Re=2L/Δt，具有强阻尼特性，可抑制谐振电路中的高频振荡，但对低频振荡的抑制作用不佳[24]，并且也会将原来仿真中本应该显现出来的振荡消除而不能够准确的模拟整个仿真过程。因相位畸变的存在，在仿真的过程中也不能完全使用后向欧拉法，另外因阻抗的存在也会使得离散系统产生附加损耗。

对图2中RL电路采用后向欧拉法时仿真波形如图5所示。

图5 后向欧拉法

后向欧拉法通常结合其他方法一同仿真计算发挥其稳定性优势，基于MATLAB平台的PSAT(Power System Analysis Toolbox)仿真工具即采用了此方法[25]。

3.2 梯形积分法

3.2.1 隐式梯形积分法

隐式梯形积分法的积分格式为：

(9)

局部截断误差采用泰勒级数展开方式分析如式(10)所示:

(10)

与泰勒级数相对比可知，隐式梯形积分法的局部截断误差为O(h3)，则P值为2，因此具有2阶精度。

对图2中RL电路采用隐式梯形积分法时仿真波形如图6所示。

图6 隐式梯形积分法

图7 隐式梯形积分法的稳定域

对任意仿真步长均可以保证截断误差衰减，隐式梯形积分法是A稳定方法中具有最小误差常数的一种数值方法[26]，在仿真中可利用其稳定性优势，联立同步发电机和网络的差分与代数方程计算后，可有效消除接口误差[24]。但由于梯形法不是L稳定性的，因此在仿真中若出现因网络拓扑结构变化而导致非状态变量突变时，将产生不衰减的数值振荡现象[27-28]。

梯形积分法的研究已经比较成熟，但许多研究者将其与其他方法相组合形成新方法，为目前实时仿真分析提供了更多的扩展空间。针对梯形积分法产生数值振荡问题，文献[29-30]提出了基于梯形积分法的改进算法，即阻尼梯形法，此方法将梯形法和后向欧拉法加权混合，相当于在电感上并联小电导，在电容上串联小电阻。阻尼梯形法引入了阻尼因子α，精度可介于梯形法与后向欧拉法之间，该方法的精度与阻尼因子α的大小选取有关[30]，然而阻尼因子的选取并不是能够完全确定的，都是通过估计得来，因此文献[31]在分析阻尼梯形法误差的基础上，提出了阻尼梯形法的补偿修正法，提高了仿真精度。文献[32]提出的龙-库-梯法将龙格-库塔法与梯形法相结合，相互弥补不足，使仿真更加精确可靠，对数值振荡也有良好的衰减作用。

隐式梯形积分法在求解过程中，其浮点数运算量与后向欧拉法相当，计算效率高。由于隐式梯形积分法对于电力系统这样的刚性系统具备良好的适应性，因此在实际仿真中被广泛应用。传统电磁暂态EMTP仿真程序因其精度高与稳定性好等优势而广泛应用。另外电力系统商业计算程序BPA、PSASP也釆用了隐式梯形积分法。

3.2.2 均值法

为消除非原型数值振荡，EMTP曾采用算后处理法，即在仿真计算过程中依旧采用梯形积分法，不对数值振荡做特别处理，而在输出时将振荡的平均值作为仿真结果。如式(11)所示，由梯形积分法在t-Δt时刻步长与t时刻步长得到的值取平均后输出[33]。与隐式梯形积分法对比表明具有便捷性和有效性，计算简单，对数值振荡有一定抑制效果，但稳定性较差。

(11)

对图2中RL电路采用均值法时仿真波形如图8所示。

图8 均值法

3.3 CDA法

临界阻尼调整法(Critical Damping Adjustment, CDA)是由文献[34]提出。采用变步长，计算复杂，实际上是一种积分算法相结合的方式，仍以梯形法作为方程求解计算的基础，主要计算步骤是仅在网络突变发生时，采用后向欧拉法在两个半步长中进行计算，此方法将数值振荡从一开始就得以消除，因此对数值振荡有一定抑制作用[35-36]。CDA法利用其临界阻尼特性，避免了电导矩阵的变化，为电路中的所有元件带来相同的等效电导，节省了计算工作量[37]，但算法的相互转换却导致系统编程更为繁琐，切换时刻的确定造成了计算负担，影响仿真速度，计算效率低。受组合方式的启发，文献[38]将一种新型4步线性多步法与梯形积分法结合切换计算，精度更高且能有效避免数值振荡。

采用CDA法消除数值振荡的前提条件是需检测出突变现象的发生时刻，而在实际情况下突变类型多种多样，如控制系统中电压源和电流源因限幅环节的影响，难以检测到突变现象情况[39]，此时CDA法无法准确抑制数值振荡。另外因暂态仿真中延迟出现的突变时刻，会造成CDA法无法准确检测与抑制振荡，因此文献[40]提出基于2级3阶Radau II A法的并行求解方法，无需检测突变时刻即可消除振荡，且并行求解方式也提高了计算效率，该方法也适用于线性与非线性常微分初值问题。

对图2中RL电路采用CDA法时仿真波形如图9所示。

图9 临界阻尼调整法

另外CDA法中包含后向欧拉法，在长期模拟结果中会引入相移问题，计算时会导致精度降低并产生误差[41]。在含有大量电力电子开关器件开关操作时，CDA法只有在正常计算时可保持2阶精度，在切换算法消除数值振荡过程中，算法精度会降至1阶。因此文献[42]提出了3S-DIRK变阶变步长方法，可在正常计算时具备3阶精度，其余计算可切换至具有L稳定性的算法抑制数值振荡并可保证精度不低于2阶。

实际应用中，软件Microtran实现了CDA技术[43]，EMTP的DCG-EPRI版本(即EMTP版本3.0)中也采用了临界阻尼调整法，虽然消除了数值振荡，但仿真波形中会出现一些处理痕迹，结果并不理想[44]。

3.4 THTA法

THTA法(Trapezoidal History Term Averaging)采用变步长，具备A稳定性，计算较复杂。此方法仍以梯形积分法为基础，在Δt/2的两个步长中确定历史源值的平均值。如式(12)所示，先采用梯形积分法在时间t中计算历史电流源值，再计算平均值ITHTA(t)，并将此值作为新的历史电流源值计算VLTHTA(t)，第二个半步长重复以上步骤计算，在数值稳定后，继续采用隐式梯形积分法求解。

(12)

采用THTA法对各步长中的参数分析列举，如表1所示。

表1 THTA法对各参数分析

对图2中RL电路采用THTA法时仿真波形如图10所示。

该方法在仿真过程中保持了梯形积分法精度高的优势，消除了误差，仅在Δt/2步长中出现一次抖动，之后保持稳定，对数值振荡有一定抑制作用。THTA法与CDA法不同的是，它只需上一步骤的历史源值产生一个新的值就能免受数值振荡，计算效率更高。

图10 THTA法

4 新型积分算法研究

目前大规模新能源发电接入电力系统后，呈现出电力电子化趋势[45]，大量具有宽频带响应特性的电力电子装置不断接入至电力系统中，极大改变了电力系统的动态特性[46]，使电磁暂态仿真的分析复杂性不断增加。电力系统电磁暂态仿真在经典积分算法的基础上不断发展改进，目前越来越来多的新型积分算法和技术不断提出并获得应用。

4.1 块广义向后差分法

块广义向后差分方法(Generalized Backward Differentiation Formula,GBDF)是块边界值方法中的一类，由文献[47]提出。具备A稳定性和L稳定性[48]。精度可在1阶(也称为后向欧拉法)～5阶之间变化，采用2或3阶的GBDF方法时可求解初值问题，因具有较强的阻尼特性，可以得到在时间上很平滑的解，因此在电磁暂态仿真中可消除数值振荡现象。当采用s级s阶的GBDF方法求解式(6)一阶微分方程的初值问题，计算格式可描述为式(13)～式(15)：

(13)

式中N为时间离散网格点数；h为时间积分步长；αi,i∈(0,s)的取值可参考文献[49]。

(14)

(15)

同时采用Brugnano等提出的附加方法选择策略与边界值方法联合求解待求点数N相匹配的N个方程，将不会对数值稳定性与计算精度产生影响。采用块广义向后差分法可允许较大的时间积分步长求解离散后的一阶模型得到离散点的时域解，相较于经典积分算法中的临界阻尼调整法，无需检测网络结构的突变现象，减轻了计算负担，且能有效避免数值振荡，适用于输电线路电磁暂态仿真中[50]。

电力系统规模与复杂性的增加，对电力系统电磁暂态实时仿真要求也更加严格，为达到实时仿真分析计算要求，文献[51]将块广义向后差分法与扩展的隐式梯形积分法结合应用在非线性常微分初值问题的数值计算中，解决了高维非线性初值问题面临的维数灾问题。文献[52]利用块广义向后差分法，采用矩阵分解，不会因整体雅可比矩阵或多个分块子矩阵被三角分解而降低效率，数值计算效率比经典微分求积法更高。

4.2 精细积分法

电力系统电磁暂态仿真过程需要解决线性或非线性系统，但最主要的还是先将线性时不变的系统精细的分析好[51]。文献[53]以Duhamel积分为基础，提出了矩阵指数计算的精细积分方法，在此基础上又提出解决瞬态初始问题时保持高精度、高稳定性的求解方法。对于线性常微分方程，其解基于Duhamel积分可表示为式(16)所示：

(16)

通过等时间步长η将时间划分，即可得到式(16)的递推形式：

(17)

其中，eHη=I+Ht+(Hη)2/2+(Hη)3/3!+…；H为常数矩阵，非齐次项f为时间函数。选取不同的Duhamel数值积分方法可衍生出精细积分法的其他形式，虽然本质上为显示积分，但例如高精度直接积分法利用拉格朗日高阶插值多项式，在处理非线性问题上，精度和稳定性依然有优势[54]。

精细积分法将微分变量的导数部分细分为线性部分和非线性部分，是用于计算矩阵指数的一种方法，对于非齐次方程，通过增维方法，将其扩容为齐次方程，最后采用精细积分法求解时，可避免矩阵的求逆计算，提高了计算速度，但精度不高。当在求解非线性刚柔耦合问题时，精细积分法处理刚性问题也有独特优势[55]。文献[56]通过将精细积分法进行增维的方式应用到电磁暂态计算中，一定程度上解决了数值振荡问题，但对于大规模电力系统，庞大的计算量增加了误差，计算精度降低。文献[57]提出一种基于泰勒级数的精细计算方法，整合的过程划分为更小的部分，运用泰勒级数展开计算指数矩阵，但存在严重的舍入误差，因此文献[58]提出改进的泰勒级数精细积分法，并将此运用到电磁暂态数值计算中，相较于经典积分算法，既保留了算法的优点，又能够避免数值振荡问题，计算精度也得到提高。为处理采用隐式梯形积分法产生的问题，文献[59]在精细积分算法的基础上，提出了隐式精细积分算法，比经典积分算法精度更高，更易于实现。

精细积分法有较好的稳定性，但在时间上划分的许多区间会使计算效率降低，在实际应用中影响很大，因此文献[60]采用组合稀疏矩阵技术，提出了精细积分级数解的并行算法。为解决电力系统输电线中过电压数值模拟的计算效率问题，在精细积分法的基础上，文献[61]提出高精度直接积分法，相较于CDA法小步长仿真，大步长计算效率更高，具有良好的数值稳定性。为避免状态矩阵求逆影响仿真速度，文献[62]提出基于精细积分法的多步法，具备良好的稳定性和高精度优势，也极大地提高了仿真效率。另外文献[63]将非线性部分作线性化假设，使增维的精细积分法拓展到了非线性系统研究领域，应用领域进一步拓宽。

4.3 根匹配技术

基于截断泰勒级数的积分方法在仿真阶跃响应时，会产生数值振荡问题，相较于经典积分算法，差分方程的指数形式(Root-matching Techniques，根匹配技术)是一种替代方法。

采用指数形式能够在循环计算之前，预先计算好指数项并将其存储，计算效率并比隐式梯形积分法在数值计算上更高效[64]。该方法通过确定电路结构计算得到其S域的传递函数H(s)，由传递函数可得到其差分方程的指数形式，与Dommel方法得到的差分方程类似，可将差分方程的指数形式看作一个诺顿等效形式。文章以RL电路为例，其传递函数如式(18)所示:

(18)

差分方程指数形式的诺顿表达式如式(19)所示：

(19)

传统电磁暂态仿真直接在时域内进行差分求解，而根匹配法从频域离散系统出发，通过反变换进行差分化，具有良好的稳定性，而且保持了原连续模型的动态特性，精度高于传统积分算法。然而根匹配法并非适合任何元件，针对RL并联与RC串联回路并不能有效避免非原型振荡以及根匹配法无法对独立电感、电容元件进行差分处理的问题，文献[65]提出了基于离散相似的电磁暂态仿真方法，首先将组合形式拆分为独立元件，再添加正负虚拟电阻转换为RL串联与RC并联的组合元件，经过离散相似法差分后，历史项中不含非状态变量，因此有效避免了数值振荡。

根匹配技术总结有以下优点[66]：(1)能够消除截断误差，并因此消除了数值振荡。与采用的步长无关，能够同时应用于电网络方程和控制模块；(2)能够等效成诺顿等效电路，可与数值积分代换方法完全兼容，并且矩阵求解技术保持不变；(3)能够提供高效和准确的时域仿真解，目前可应用于PSCAD和EMTDC仿真软件中。

5 结束语

文章根据电磁暂态仿真分析的高性能需求讨论了经典积分算法和新型积分算法，主要结论如下：

(1)后向欧拉法具备A稳定性和L稳定性，可避免数值振荡，但精度低，阻抗的存在也会产额外损耗。隐式梯形积分法具备A稳定性和2阶精度，计算效率高，在电力系统仿真中应用最为广泛，但不是L稳定性，因此不具备消除数值振荡的能力。临界阻尼调整法能够从源头上消除数值振荡，但需检测突变发生时刻，算法切换增加了计算负担，在计算精度和计算效率上不占优势。以梯形法为基础的THTA法仅一个时间步长即可消除数值振荡并保持稳定，比CDA法计算更简便有效;

(2)经典积分算法的研究主要集中在不同特性算法的组合运用以发挥各自优势，为电力系统电磁暂态仿真提供了更多发展空间。如受临界阻尼调整法的启发，将隐式梯形积分法与一类新的低阶、L稳定的线性多步法相结合，比较适用于含有大量电力电子装置的电力系统电磁暂态仿真中;

(3)新型积分算法中精细积分法具备高精度和良好的稳定性。为满足实时分析计算需求，在求解大规模问题中可采用并行计算技术提高精细积分法的计算效率。精细积分法本质上适用于线性定常部分的求解，利用其处理线性刚性问题的优势可构造出更高性能数值方法。如何利用精细积分法优点求解非线性问题，在高性能电力系统电磁暂态仿真中尚有研究价值。根匹配技术不仅可以形成与Dommel算法结构相似的诺顿等效电路，可较为容易的应用在现有电磁暂态程序中并有效避免了非原型数值振荡，也在精度和稳定性方面大大提高，但根匹配法为恒稳定算法，并非适用于任何元件，针对基于根匹配法的非线性以及存在耦合关系的元件处理还需进一步研究。