基于2级3阶单对角隐式Runge-Kutta法的电磁暂态计算方法
【技术领域】
[0001] 本发明涉及一种电力系统电磁暂态数值计算方法，具体是涉及一种基于2级3阶、 非线性代数稳定的单对角隐式龙格-库塔（Runge-Kutta，RK)方法的电磁暂态计算方法。
【背景技术】
[0002] 电磁暂态计算是研究电力系统各个元件中电场和磁场以及相应的电压和电流的 变化过程，其主要目的在于研究电力系统故障或操作过后可能出现的暂态过电压和过电 流。
[0003] 在电磁暂态分析计算过程中，考虑到元件的非线性、电磁耦合、长线路波过程和线 路三相不对称、线路参数的频率特性等因素，需建立系统元件的微分方程或偏微分方程，并 借助一定的数值计算方法对这些方程进行离散化，得到代数形式的差分方程，进而求解各 时间点的物理量。
[0004] 目前，经典的电磁暂态数值计算主要是采用隐式梯形积分方法。隐式梯形积分方 法具有2阶精度以及A-稳定性，但该方法不是L-稳定的。在电磁暂态数值计算过程中， 当发生电感电流或电容电压突变以及开关元件动作等情况时，使用隐式梯形法进行电磁 暂态数值计算时会产生一系列"虚假的"、持续的数值振荡。为了解决此问题，加拿大学者 J.R.Marti和我国学者林集明等将隐式梯形法与具有强阻尼特性的隐式欧拉法相结合，提 出了临界阻尼调整法（CriticalDampingAdjustment，CDA)，并将此方法运用到电磁暂态 仿真程序（ElectromagneticTransientProgram，EMTP)中。隐式欧拉法既是A-稳定的也 是L-稳定的数值方法，可以避免数值振荡问题，但它只是1阶的数值方法。在使用CDA方 法进行电磁暂态仿真时，依旧以隐式梯形法作为数值计算的主要方法，仅在突变发生时刻 利用隐式欧拉法来进行计算，以此避免数值振荡情况的发生。CDA方法可以有效避免数值振 荡问题，但前提条件是需要检测出突变现象及其发生的时刻。在电磁暂态数值计算中，突变 现象主要包括开关元件的动作，电感电流和电容电压的突变，以及非线性电感、电容元件的 运行方式由分段线性曲线转折点的一边跃变至另一边等。在实际仿真过程中，突变的种类 众多，对某些突变现象很难准确判断出突变发生的时刻。例如，当传输线始端电压或电流发 生突变时，很难准确地判定其末端电压或电流发生突变的时刻，或者是控制系统中电压源 和电流源因为限幅环节的影响发生难以检测到的突变等。因此，对一些很难检测到的突变 现象，CDA方法仍然无法避免数值振荡问题，例如，当一个阶跃突变信号从传输线始端传输 到末端时（如图1所示），若采用CDA方法进行电磁暂态计算，因不能有效检测末端的电压 突变现象（如图2 (a)所示），仍然会发生数值振荡（如图2 (b)所示）。
[0005] 为进一步解决CDA方法的问题，日本研究人员TakuNoda等将2级2阶单对角隐 式RK方法应用于电磁暂态数值计算。该方法可用Butcher表（ButcherTable)表示如下：
[0006]
[0007] 与隐式欧拉法一样，上述方法（1)也是L-稳定的数值方法。因此，该方法可以避 免数值振荡问题。与CDA方法相比较，上述方法（1)无需对突变现象进行检测或判断，这是 该方法的主要优点。然而，无论有无突变现象，方法（1)均须在每个时间步长内计算两个内 点的变量值，其每一步的积分相当于隐式欧拉方法用更小的步长连续积分2步。因此，其计 算量约为隐式梯形法的2倍，也比CDA方法的计算效率更低。
【发明内容】
[0008] 本发明所要解决的技术问题，就是提供一种基于2级3阶单对角隐式Runge-Kutta 法的电磁暂态数值计算方法，其能在不降低经典电磁暂态数值计算方法计算效率的前提 下，解决隐式梯形积分方法所存在的数值振荡问题，且计算效率比CDA方法以及基于2级2 阶单对角隐式Runge-Kutta方法的电磁暂态数值计算方法的效率更高。
[0009] 解决上述技术问题，本发明采取如下的技术方案：
[0010] -种基于2级3阶单对角隐式Runge-Kutta法的电磁暂态数值分析方法，其特征 在于：通过建立电力系统电磁暂态数值计算的时域微分方程，采用B-稳定的2级3阶单对 角隐式Runge-Kutta法进行时域数值积分计算，逐步求解出各物理量随时间的变化曲线；
[0011] 完整的步骤包括：
[0012] 1)输入原始数据，建立各元件的微分方程，形成电磁暂态数值计算的基本数学模 型（土 = /(i，JC));
[0013] 2)电磁暂态数值计算初始化
[0014]置t= 0· 0s，积分步数η= 0;
[0015] 确定数值积分步长h、电磁暂态数值计算时程Τ;
[0016] 确定各状态变量的初值，即x(t= 0) =X。；
[0017] 输入电磁暂态数值计算的故障或操作；
[0018] 3)故障或操作判断
[0019] 依据时刻t判断系统有无故障或操作；
[0020] 若有故障或操作，则修改相应的微分方程以及相应的状态变量值xn(t);
[0021] 4)数值积分
[0022] 采用B-稳定的2级3阶单对角隐式RK方法，计算出状态变量在t=t+h处的值 -^-n+1?
[0023] 5)t=t+h；n=n+1 ；
[0024] 6)数值计算是否终止判断
[0025]若t<T，则返回步骤3)，继续下一步即下一时刻的数值计算；
[0026] 若t彡Τ，则转至步骤7);
[0027] 7)数值计算结果输出。
[0028] 所述的步骤4)数值积分中，用到的2级3阶单对角隐式RK方法的Butcher表为：
[0029]
[0030] 与隐式梯形积分方法以及方法（1) 一样，方法（2)也是A-稳定的数值方法，但方 法（2)不是严格意义上的L-稳定的数值方法，而是非线性B-稳定亦即非线性代数稳定的 数值方法；
[0031] 对常微分方程初值问题：
[0033] 所述的步骤4)数值积分也即对常微分方程初值问题（3)的求解步骤具体如下：
[0034] 从tjljtn+1时亥I」，已知状态变量X⑴在t=t"时刻的值Xn，求解其在t=tn+1时 亥1J的值xn+1;令计算的时间步长为h=tn+1_tn;
[0035] 第一步：计算状态变量在第一个内点处的近似值：
[0037] 其中，充》奴可)是状态变量在内点?= 〖"+βΑ处的近似值；
[0038] 若f(t，X)是X的线性函数，则依据方程⑷可直接求解出萬，由此得/(?)的 值；
[0039] 若f(t，X)是X的非线性函数，则方程⑷的求解采用牛顿迭代法，同样可求解出 萬及相应的
[0040] 第二步：计算状态变量在第二个内点处的近似值：
[0042] 上式中，/(??)为已知量；?Χ(ξ)是状态变量在内点ξ:? 处的近似 值；
[0043] 同理，若f(t，x)是X的线性函数，则依据方程（5)直接求解出毛，由此得/(^?) 的值；若f(t，X)是X的非线性函数，则方程（5)的求解采用牛顿迭代法，同样可求解出％及 相应的/(1么）；
[0044] 第三步：计算状态变量在t=tn+1时刻的值：
[0046] 本发明的理论基础：
[0047] 众所周知，L-稳定性是线性稳定性的范畴。从理论上讲，对线性微分动力系统， L-稳定的数值方法可以避免数值振荡问题，这是隐式欧拉法以及方法（1)能够避免数值振 荡问题的主要数学机理。但对非线性微分动力系统，L-稳定的数值方法并不一定能完全避 免数值振荡问题。为此，研究人员已建立了非线性稳定性分析的概念及相关理论体系。关 于非线性稳定性，一个重要的结论就是所谓的B-稳定性以及非线性代数稳定性。研究人员 已证明：对一个非退化的Runge-Kutta方法，非线性代数稳定性等价于B-稳定性。利用非 线性代数稳定性的定义，可以验证并得出以下结论：
[0048] 2级3阶单对角隐式RK方法（即方法⑵）的Μ矩阵的表达式为：
[0050] 上式（7)中，B=diag(b)。显然，Μ矩阵的特征值为：λ丨=〇， 此，数值方法（2)的Μ矩阵是非负定的。由于有匕=132= 1/2 >0,因此，上述2级3阶单 对角隐式RK方法（2)是非线性代数稳定的，也是Β-稳定的。
[0051] 众所周知，对常微分方程初值问题（方程（3))，Β-稳定的数值方法满足单边 Lipschitz条件，即：
[0052]<f(t1;Xj)-f(t2,x2),X!-x2> ^ 〇 (8)；
[0053] 上式⑶中，〈·，·>表示内积。因此，B-稳定的数值方法具有能量耗散性。从 物理概念上讲，这里的能量耗散性也就是非线性阻尼特性。因此，在系统发生突变时，B-稳 定的数值方法不会产生数值振荡问题。换言之，在电磁暂态数值计算中，若系统发生突变现 象，2级3阶单对角隐式RK方法（2)可以避免数值振荡问题。这就是本发明的理论基础。
[0054] 下面给出2级3阶单对角隐式RK方法（2)不会产生数值振荡的几个具体实例。
[0055] 图3a是一个基本的线性R-L串联电路，图3b是施加的电流源；其中，开关K在施 加的电流i(t)降至零时（t= 0. 01秒）突然打开。图4a是利用隐式梯形法（计算步长h =0. 05ms)对该测试电路进行数值计算的结果（产生严重的数值振荡）；图4b是利用2级 3阶单对角隐式RK方法（计算步长h= 0.lms)进行数值计算的结果。显然，从图4b可以 看出：当突变发生时，2级3阶单对角隐式RK方法没有产生数值振荡。
[0056] 图5是一个R-L串联电路，其中，电感部分由一个线性电感外加一个饱和电抗器组 成；开关K在t= 0秒时突然合闸。图6a是利用隐式梯形法（计算步长h= 0. 05ms)对该 测试电路进行数值计算的结果（产生数值振荡）；图6b是利用2级3阶单对角隐式RK方 法（计算步长h= 0.lms)进行数值计算的结果。显然，从图6b可以看出：当突变发生时， 2级3阶单对角隐式RK方法没有产生数值振荡。
[0057] 由于隐式梯形法以及方法（1)均是2阶的数值方法，而2级3阶单对角隐式RK方 法则是3阶的数值方法。因此，在采用相同的计算步长的情况下，2级3阶单对角隐式RK方 法的计算精度要比隐式梯形积分方法以及方法（1)更高。很易理解，在满足相同的计算精 度的前提下，2级3阶单对角隐式RK方法可以采用比隐式梯形法更大的步长；在2级3阶 单对角隐式RK方法采用2倍于隐式梯形积分法的步长的情况下，若两者的计算精度大致相 同，则它们的计算效率大致相当。为此，以图7所示的基本测试电路（可以获得该测试电路 的精确解即解析解）为例，对2级3阶单对角隐式RK方法相对于隐式梯形法的计算效率进 行了测试和评估。图8是分别利用隐式梯形法（步长h= 0. 01ms)和2级3阶单对角隐式 RK方法（步长h= 0.02ms)对图7所示测试电路进行数值计算的误差对比曲线。显然，从 图8可以看出：