CN105404610A - 基于2级3阶单对角隐式Runge-Kutta法的电磁暂态计算方法 - Google Patents

Abstract

一种基于2级3阶单对角隐式Runge-Kutta法的电磁暂态数值分析方法，其特征在于：通过建立电力系统电磁暂态数值计算的时域微分方程，采用B-稳定的2级3阶单对角隐式Runge-Kutta法进行时域数值积分计算，逐步求解出各物理量随时间的变化曲线。本发明与经典的、基于隐式梯形积分方法的电磁暂态数值计算方法相比较，基于2级3阶单对角隐式RK方法的电磁暂态数值计算方法可以避免数值振荡问题；与CDA方法相比较，基于2级3阶单对角隐式RK方法的电磁暂态数值计算方法不仅可以完全避免数值振荡问题，而且计算效率更高；与数值方法(1)相比较，基于2级3阶单对角隐式RK方法的电磁暂态数值计算方法的计算效率明显更高。

Description

基于2级3阶单对角隐式Runge-Kutta法的电磁暂态计算方法

技术领域

本发明涉及一种电力系统电磁暂态数值计算方法，具体是涉及一种基于2级3阶、非线性代数稳定的单对角隐式龙格-库塔(Runge-Kutta，RK)方法的电磁暂态计算方法。

背景技术

电磁暂态计算是研究电力系统各个元件中电场和磁场以及相应的电压和电流的变化过程，其主要目的在于研究电力系统故障或操作过后可能出现的暂态过电压和过电流。

在电磁暂态分析计算过程中，考虑到元件的非线性、电磁耦合、长线路波过程和线路三相不对称、线路参数的频率特性等因素，需建立系统元件的微分方程或偏微分方程，并借助一定的数值计算方法对这些方程进行离散化，得到代数形式的差分方程，进而求解各时间点的物理量。

目前，经典的电磁暂态数值计算主要是采用隐式梯形积分方法。隐式梯形积分方法具有2阶精度以及A-稳定性，但该方法不是L-稳定的。在电磁暂态数值计算过程中，当发生电感电流或电容电压突变以及开关元件动作等情况时，使用隐式梯形法进行电磁暂态数值计算时会产生一系列“虚假的”、持续的数值振荡。为了解决此问题，加拿大学者J.R.Marti和我国学者林集明等将隐式梯形法与具有强阻尼特性的隐式欧拉法相结合，提出了临界阻尼调整法(CriticalDampingAdjustment，CDA)，并将此方法运用到电磁暂态仿真程序(ElectromagneticTransientProgram，EMTP)中。隐式欧拉法既是A-稳定的也是L-稳定的数值方法，可以避免数值振荡问题，但它只是1阶的数值方法。在使用CDA方法进行电磁暂态仿真时，依旧以隐式梯形法作为数值计算的主要方法，仅在突变发生时刻利用隐式欧拉法来进行计算，以此避免数值振荡情况的发生。CDA方法可以有效避免数值振荡问题，但前提条件是需要检测出突变现象及其发生的时刻。在电磁暂态数值计算中，突变现象主要包括开关元件的动作，电感电流和电容电压的突变，以及非线性电感、电容元件的运行方式由分段线性曲线转折点的一边跃变至另一边等。在实际仿真过程中，突变的种类众多，对某些突变现象很难准确判断出突变发生的时刻。例如，当传输线始端电压或电流发生突变时，很难准确地判定其末端电压或电流发生突变的时刻，或者是控制系统中电压源和电流源因为限幅环节的影响发生难以检测到的突变等。因此，对一些很难检测到的突变现象，CDA方法仍然无法避免数值振荡问题，例如，当一个阶跃突变信号从传输线始端传输到末端时(如图1所示)，若采用CDA方法进行电磁暂态计算，因不能有效检测末端的电压突变现象(如图2(a)所示)，仍然会发生数值振荡(如图2(b)所示)。

为进一步解决CDA方法的问题，日本研究人员TakuNoda等将2级2阶单对角隐式RK方法应用于电磁暂态数值计算。该方法可用Butcher表(ButcherTable)表示如下：

与隐式欧拉法一样，上述方法(1)也是L-稳定的数值方法。因此，该方法可以避免数值振荡问题。与CDA方法相比较，上述方法(1)无需对突变现象进行检测或判断，这是该方法的主要优点。然而，无论有无突变现象，方法(1)均须在每个时间步长内计算两个内点的变量值，其每一步的积分相当于隐式欧拉方法用更小的步长连续积分2步。因此，其计算量约为隐式梯形法的2倍，也比CDA方法的计算效率更低。

发明内容

本发明所要解决的技术问题，就是提供一种基于2级3阶单对角隐式Runge-Kutta法的电磁暂态数值计算方法，其能在不降低经典电磁暂态数值计算方法计算效率的前提下，解决隐式梯形积分方法所存在的数值振荡问题，且计算效率比CDA方法以及基于2级2阶单对角隐式Runge-Kutta方法的电磁暂态数值计算方法的效率更高。

解决上述技术问题，本发明采取如下的技术方案：

一种基于2级3阶单对角隐式Runge-Kutta法的电磁暂态数值分析方法，其特征在于：通过建立电力系统电磁暂态数值计算的时域微分方程，采用B-稳定的2级3阶单对角隐式Runge-Kutta法进行时域数值积分计算，逐步求解出各物理量随时间的变化曲线；

完整的步骤包括：

1)输入原始数据，建立各元件的微分方程，形成电磁暂态数值计算的基本数学模型 $(\stackrel{\·}{x}=f(t,x\left)\right);$

2)电磁暂态数值计算初始化

置t＝0.0s，积分步数n＝0；

确定数值积分步长h、电磁暂态数值计算时程T；

确定各状态变量的初值，即x(t＝0)＝x₀；

输入电磁暂态数值计算的故障或操作；

3)故障或操作判断

依据时刻t判断系统有无故障或操作；

若有故障或操作，则修改相应的微分方程以及相应的状态变量值x_n(t)；

4)数值积分

采用B-稳定的2级3阶单对角隐式RK方法，计算出状态变量在t＝t+h处的值x_n+1；

5)t＝t+h；n＝n+1；

6)数值计算是否终止判断

若t＜T，则返回步骤3)，继续下一步即下一时刻的数值计算；

若t≥T，则转至步骤7)；

7)数值计算结果输出。

所述的步骤4)数值积分中，用到的2级3阶单对角隐式RK方法的Butcher表为：

与隐式梯形积分方法以及方法(1)一样，方法(2)也是A-稳定的数值方法，但方法(2)不是严格意义上的L-稳定的数值方法，而是非线性B-稳定亦即非线性代数稳定的数值方法；

对常微分方程初值问题：

$\left\{\begin{array}{c}\frac{d}{dt}x\left(t\right)=\stackrel{\·}{x}\left(t\right)=f(t,x(t\left)\right)\\ x(t=0)=x_0\end{array}\right.---\left(3\right);$

所述的步骤4)数值积分也即对常微分方程初值问题(3)的求解步骤具体如下：

从t_n到t_n+1时刻，已知状态变量x(t)在t＝t_n时刻的值x_n，求解其在t＝t_n+1时刻的值x_n+1；令计算的时间步长为h＝t_n+1-t_n；

第一步：计算状态变量在第一个内点处的近似值：

${\tilde{x}}_1=x_n+h\mathrm{\β}f({\tilde{t}}_1,{\tilde{x}}_1)---\left(4\right);$

其中，是状态变量在内点处的近似值；

若f(t,x)是x的线性函数，则依据方程(4)可直接求解出由此得的值；

若f(t,x)是x的非线性函数，则方程(4)的求解采用牛顿迭代法，同样可求解出及相应的

第二步：计算状态变量在第二个内点处的近似值：

${\tilde{x}}_2=x_n+h(1-2\mathrm{\β})f({\tilde{t}}_1,{\tilde{x}}_1)+h\mathrm{\β}f({\tilde{t}}_2,{\tilde{x}}_2)---\left(5\right);$

上式中，为已知量；是状态变量在内点处的近似值；

同理，若f(t,x)是x的线性函数，则依据方程(5)直接求解出由此得的值；若f(t,x)是x的非线性函数，则方程(5)的求解采用牛顿迭代法，同样可求解出及相应的

第三步：计算状态变量在t＝t_n+1时刻的值：

$x_{n+1}=x_n+h\[\frac{1}{2}f({\tilde{t}}_1,{\tilde{x}}_1)+\frac{1}{2}f({\tilde{t}}_2,{\tilde{x}}_2)\]---\left(6\right).$

本发明的理论基础：

众所周知，L-稳定性是线性稳定性的范畴。从理论上讲，对线性微分动力系统，L-稳定的数值方法可以避免数值振荡问题，这是隐式欧拉法以及方法(1)能够避免数值振荡问题的主要数学机理。但对非线性微分动力系统，L-稳定的数值方法并不一定能完全避免数值振荡问题。为此，研究人员已建立了非线性稳定性分析的概念及相关理论体系。关于非线性稳定性，一个重要的结论就是所谓的B-稳定性以及非线性代数稳定性。研究人员已证明：对一个非退化的Runge-Kutta方法，非线性代数稳定性等价于B-稳定性。利用非线性代数稳定性的定义，可以验证并得出以下结论：

2级3阶单对角隐式RK方法(即方法(2))的M矩阵的表达式为：

$M\≡\[M_{ij}\]=BA+A^{T}B-{\mathrm{bb}}^T=\left[\begin{array}{cc}\mathrm{\β}-\frac{1}{4}& -\mathrm{\β}+\frac{1}{4}\\ -\mathrm{\β}+\frac{1}{4}& \mathrm{\β}-\frac{1}{4}\end{array}\right]---\left(7\right);$

上式(7)中，B＝diag(b)。显然，M矩阵的特征值为：λ₁＝0，因此，数值方法(2)的M矩阵是非负定的。由于有b₁＝b₂＝1/2＞0，因此，上述2级3阶单对角隐式RK方法(2)是非线性代数稳定的，也是B-稳定的。

众所周知，对常微分方程初值问题(方程(3))，B-稳定的数值方法满足单边Lipschitz条件，即：

<f(t₁,x₁)-f(t₂,x₂)，x₁-x₂>≤0(8)；

上式(8)中，<·，·>表示内积。因此，B-稳定的数值方法具有能量耗散性。从物理概念上讲，这里的能量耗散性也就是非线性阻尼特性。因此，在系统发生突变时，B-稳定的数值方法不会产生数值振荡问题。换言之，在电磁暂态数值计算中，若系统发生突变现象，2级3阶单对角隐式RK方法(2)可以避免数值振荡问题。这就是本发明的理论基础。

下面给出2级3阶单对角隐式RK方法(2)不会产生数值振荡的几个具体实例。

图3a是一个基本的线性R-L串联电路，图3b是施加的电流源；其中，开关K在施加的电流i(t)降至零时(t＝0.01秒)突然打开。图4a是利用隐式梯形法(计算步长h＝0.05ms)对该测试电路进行数值计算的结果(产生严重的数值振荡)；图4b是利用2级3阶单对角隐式RK方法(计算步长h＝0.1ms)进行数值计算的结果。显然，从图4b可以看出：当突变发生时，2级3阶单对角隐式RK方法没有产生数值振荡。

图5是一个R-L串联电路，其中，电感部分由一个线性电感外加一个饱和电抗器组成；开关K在t＝0秒时突然合闸。图6a是利用隐式梯形法(计算步长h＝0.05ms)对该测试电路进行数值计算的结果(产生数值振荡)；图6b是利用2级3阶单对角隐式RK方法(计算步长h＝0.1ms)进行数值计算的结果。显然，从图6b可以看出：当突变发生时，2级3阶单对角隐式RK方法没有产生数值振荡。

由于隐式梯形法以及方法(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可以看出：在采用2倍于隐式梯形法的计算步长的情况下，2级3阶单对角隐式RK方法的计算精度比隐式梯形法略高。因此，2级3阶单对角隐式RK方法的计算效率与隐式梯形积分法大致相当，但比2阶的数值方法(1)的计算效率要高。

本发明的重点，就在于采用了B-稳定的2级3阶单对角隐式RK方法(2)，既避免了数值振荡问题，又具有更好的计算效率。

本发明的技术关键点或欲保护点：

一种基于2级3阶单对角隐式Runge-Kutta方法的电磁暂态数值计算方法，其特征在于：该方法是非线性B-稳定的数值方法，在系统发生突变时可以彻底避免数值振荡问题；该方法是3阶、A-稳定的数值方法，其计算效率与传统的、基于隐式梯形积分法的电磁暂态数值计算方法的效率大致相当，但比CDA方法以及基于2级2阶单对角隐式Runge-Kutta方法的电磁暂态数值计算方法的效率更高。

本发明的优点和积极效果是：与经典的、基于隐式梯形积分方法的电磁暂态数值计算方法相比较，基于2级3阶单对角隐式RK方法的电磁暂态数值计算方法可以避免数值振荡问题；与CDA方法相比较，基于2级3阶单对角隐式RK方法的电磁暂态数值计算方法不仅可以完全避免数值振荡问题，而且计算效率更高；与数值方法(1)相比较，基于2级3阶单对角隐式RK方法的电磁暂态数值计算方法的计算效率明显更高。

附图说明

图1：测试电路1(阶跃突变信号从传输线始端传至末端)示意图；

图2a：CDA方法用于测试电路1的末端电压数值计算结果；

图2b：CDA方法用于测试电路1的电容电流数值计算结果(数值振荡)；

图3a：测试电路2(基本的R-L串联电路)示意图；

图3b：施加于测试电路2(基本的R-L串联电路)的电流源变化曲线；

图4a：隐式梯形法应用于测试电路2的计算结果(数值振荡)；

图4b：2级3阶单对角隐式RK方法应用于测试电路2的计算结果；

图5a：测试电路3(非线性R-L串联电路)示意图；

图5b：测试电路3(非线性R-L串联电路)中饱和电抗器的电流-磁通曲线；

图6a：隐式梯形法应用于测试电路3的计算结果(数值振荡)；

图6b：2级3阶单对角隐式RK方法应用于测试电路3的计算结果；

图7a：测试电路4(基本的R-C串联电路)示意图；

图7b：施加于测试电路4(基本的R-C串联电路)的电压源变化曲线；

图8：两种方法对测试电路4进行数值计算的误差结果对比曲线；

图9：均匀长输电线路空载合闸示意图；

图10：均匀长输电线路离散化等值电路示意图；

图11：对均匀长输电线路空载合闸过程进行电磁暂态数值计算的实例结果。

具体实施方式

以下结合附图对本发明的实施方式进行概要说明。

本发明所提出的基于2级3阶单对角隐式RK方法的电磁暂态数值计算方法的实施流程与传统的电磁暂态数值计算方法的流程基本相同，不同之处主要是在步骤4)数值积分中所采用的积分方法不同。

以如图9所示的均匀长输电线路空载合闸的电磁暂态计算为例，本发明的具体实施步骤可以概述如下：

1)输入原始数据，建立各元件的微分方程，形成电磁暂态数值计算的基本数学模型：

输入传输线单位长度参数R₀、L₀、C₀，线路长度L，电压源函数e(t)及其内阻R_s；

空间离散化：取N＝10，将L＝300km的输电线路均匀分为10段，表示为如图11所示的空间离散化形式；离散后每一段线路的电阻、电感、电容分别为： 由此，可以建立各元件的微分方程如下：

$\left\{\begin{array}{c}l\frac{{\mathrm{di}}_k\left(t\right)}{dt}=l{\stackrel{\&CenterDot;}{i}}_k\left(t\right)=v_k\left(t\right)-v_{k+1}\left(t\right)-{\mathrm{ri}}_k\left(t\right)\\ c\frac{{\mathrm{dv}}_k\left(t\right)}{dt}=c{\stackrel{\&CenterDot;}{v}}_k\left(t\right)=i_{k-1}\left(t\right)-i_k\left(t\right)\end{array}\right.,k\&Element;(1,N)---\left(9\right);$

$c\frac{{\mathrm{dv}}_{N+1}\left(t\right)}{dt}=c{\stackrel{\&CenterDot;}{v}}_{N+1}\left(t\right)=i_N\left(t\right)---\left(10\right);$

上述方程共含2N+1个待求状态变量。将上述方程写成矩阵形式，即是：

上式(11)中：

L＝diag(l)∈R^N×N,C＝diag(c)∈R^(N+1)×(N+1),R＝diag(r)∈R^N×N(12)；

μ(t)＝[i₀(t)0…0]^T(15)；

I＝[i₁(t)…i_N(t)]^T，U＝[v₁(t)…v_N(t)v_N+1(t)]^T(16)；

显然，方程(11)最终可以写成以下的标准形式：

$\stackrel{\&CenterDot;}{x}\left(t\right)=Ax\left(t\right)+w\left(t\right)---\left(17\right);$

上式(17)中：A是一个常系数矩阵，

方程(17)即是本实施例电磁暂态数值计算的基本数学模型。

2)电磁暂态数值计算初始化：

置t＝0.0s，积分步数n＝0；

确定数值积分步长h＝1.0μs、电磁暂态数值计算时程T＝0.05s；

确定各状态变量的初值：由于是空载线路，故有x(t＝0)＝x₀≡0；

输入电磁暂态数值计算的故障或操作：

t＜0.01ms时，i₀(t)＝0；

t≥0.01ms时，e(t)＝220sin(2π×50t)。

3)故障或操作判断：

当t＜0.01ms时，没有任何操作，电磁暂态数值计算的数学模型(方程(17))不变，直接转向步骤4)；

当t＝0.01ms时，开关合闸，此时需要对微分方程(17)中的相关系数矩阵进行修改，具体情况可描述如下：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\tilde{A}x\left(t\right)+\tilde{w}\left(t\right)---\left(21\right);$

修改电磁暂态数值计算的模型后，进入下一步骤。

4)数值积分：

采用B-稳定的2级3阶单对角隐式RK方法，计算出状态变量在t＝t+h处的值x_n+1。具体情况可概述如下：

第一步：计算状态变量在第一个内点处的近似值：

${\tilde{x}}_1=x_n+h\mathrm{\&beta;}f({\tilde{t}}_1,{\tilde{x}}_1)---\left(4\right);$

上式(4)中，当t＜0.01ms时，当t≥0.01ms时， $f({\tilde{t}}_1,{\tilde{x}}_1)=\tilde{A}{\tilde{x}}_1+\tilde{w}(t_n+\mathrm{\&beta;}h).$

显然，利用方程(4)很容易解出自然也就可以获得的值。

第二步：计算状态变量在第二个内点处的近似值：

${\tilde{x}}_2=x_n+h(1-2\mathrm{\&beta;})f({\tilde{t}}_1,{\tilde{x}}_1)+h\mathrm{\&beta;}f({\tilde{t}}_2,{\tilde{x}}_2)---\left(5\right);$

上式(5)中，当t＜0.01ms时，当t≥0.01ms时， $f({\tilde{t}}_2,{\tilde{x}}_2)=\tilde{A}{\tilde{x}}_2+\tilde{w}(t_n+(1-\mathrm{\&beta;}\left)h\right).$

同样，利用方程(26)很容易解出自然也就可以获得的值。

第三步：计算状态变量在t＝t_n+1时刻的值：

$x_{n+1}=x_n+h\&lsqb;\frac{1}{2}f({\tilde{t}}_1,{\tilde{x}}_1)+\frac{1}{2}f({\tilde{t}}_2,{\tilde{x}}_2)\&rsqb;---\left(6\right).$

5)t＝t+h；n＝n+1

6)数值计算是否终止判断：

若t＜T，则返回步骤3)，继续下一步即下一时刻的数值计算；

若t≥T，则转至步骤7)。

7)数值计算结果输出：

本实施例主要是输出空载合闸是线路末端电压曲线，即v(t)＝v_N+1(t)的变化曲线，具体如图11所示。

Claims (3)

1.一种基于2级3阶单对角隐式Runge-Kutta法的电磁暂态数值分析方法，包括有数值积分步骤，其特征在于：所述的数值积分步骤采用B-稳定的2级3阶单对角隐式Runge-Kutta法进行时域数值积分计算，逐步求解出各物理量随时间的变化曲线。

2.根据权利要求1所述的基于2级3阶单对角隐式Runge-Kutta法的电磁暂态数值分析方法，其特征在于：所述的数值积分中，用到的2级3阶单对角隐式RK方法的Butcher表为：

对常微分方程初值问题：

$\left\{\begin{array}{c}\frac{d}{dt}x\left(t\right)=\stackrel{\&CenterDot;}{x}\left(t\right)=f(t,x(t\left)\right)\\ x(t=0)=x_0\end{array}\right.---\left(3\right);$

式中：为各元件的微分方程；

求解步骤具体如下：

从t_n到t_n+1时刻，已知状态变量x(t)在t＝t_n时刻的值x_n，求解其在t＝t_n+1时刻的值x_n+1；令计算的时间步长为h＝t_n+1-t_n；

第一步：计算状态变量在第一个内点处的近似值：

${\tilde{x}}_1=x_n+h\mathrm{\&beta;}f({\tilde{t}}_1,{\tilde{x}}_1)---\left(4\right);$

其中，是状态变量在内点处的近似值；

若f(t,x)是x的线性函数，则依据方程(4)直接求解出由此得的值；

若f(t,x)是x的非线性函数，则方程(4)的求解采用牛顿迭代法，同样可求解出及相应的

第二步：计算状态变量在第二个内点处的近似值：

${\tilde{x}}_2=x_n+h(1-2\mathrm{\&beta;})f({\tilde{t}}_1,{\tilde{x}}_1)+h\mathrm{\&beta;}f({\tilde{t}}_2,{\tilde{x}}_2)---\left(5\right);$

上式中，为已知量；是状态变量在内点处的近似值；

同理，若f(t,x)是x的线性函数，则依据方程(5)直接求解出由此得的值；

若f(t,x)是x的非线性函数，则方程(5)的求解采用牛顿迭代法，同样可求解出及相应的

第三步：计算状态变量在t＝t_n+1时刻的值：

$x_{n+1}=x_n+h\&lsqb;\frac{1}{2}f({\tilde{t}}_1,{\tilde{x}}_1)+\frac{1}{2}f({\tilde{t}}_2,{\tilde{x}}_2)\&rsqb;---\left(6\right).$

3.根据权利要求1和2所述的基于2级3阶单对角隐式Runge-Kutta法的电磁暂态数值分析方法，其特征在于：完整的步骤包括：

1)输入原始数据，建立各元件的微分方程，形成电磁暂态数值计算的基本数学模型 $(\stackrel{\&CenterDot;}{x}=f(t,x\left)\right);$

2)电磁暂态数值计算初始化

置t＝0.0s，积分步数n＝0；

确定数值积分步长h、电磁暂态数值计算时程T；

确定各状态变量的初值，即x(t＝0)＝x₀；

输入电磁暂态数值计算的故障或操作；

3)故障或操作判断

依据时刻t判断系统有无故障或操作；

若有故障或操作，则修改相应的微分方程以及相应的状态变量值x_n(t)；

4)数值积分

采用B-稳定的2级3阶单对角隐式RK方法，计算出状态变量在t＝t+h处的值x_n+1；

5)t＝t+h；n＝n+1；

6)数值计算是否终止判断

若t＜T，则返回步骤3)，继续下一步即下一时刻的数值计算；

若t≥T，则转至步骤7)；

7)数值计算结果输出。