CN107679287A

CN107679287A - 基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法

Info

Publication number: CN107679287A
Application number: CN201710813294.3A
Authority: CN
Inventors: 王永; 雷小林; 杨萌; 胡洁; 陶静静; 刘雪婷; 潘明帅; 饶华兴; 李强; 谭炜东; 李乾坤; 陶渊
Original assignee: China Three Gorges University CTGU
Current assignee: China Three Gorges University CTGU
Priority date: 2017-09-11
Filing date: 2017-09-11
Publication date: 2018-02-09
Anticipated expiration: 2037-09-11
Also published as: CN107679287B

Abstract

基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，通过建立电力系统电磁暂态数值计算的时域微分方程，采用A‑稳定且无限稳定的3步4阶隐式泰勒级数法进行时域数值积分计算，逐步求解各状态变量随时间的变化曲线。本发明采用的3步4阶隐式泰勒级数法是A‑稳定且无限稳定的数值方法，其对截断误差具有较快的衰减速率，可有效地抑制数值振荡，相对于隐式梯形积分法而言，其能够彻底避免数值振荡问题。此外，本发明采用的3步4阶隐式泰勒级数法的计算精度为6阶，局部截断误差为O(h⁸)，其可以通过采用较大的积分步长而提高计算效率。与CDA方法相比较，基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法不仅可以完全避免数值振荡问题，而且计算精度和效率更高、数值稳定性更强。

Description

基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法

技术领域

本发明涉及一种电力系统电磁暂态数值计算方法，具体涉及一种基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法。

背景技术

电力系统电磁暂态(Electromagnetic Transient，EMT)数值计算广泛应用于电力系统的诸多领域：过电压、绝缘配合、线路保护以及谐波分析等。其基本理论与方法已由Dommel于20世纪60年代末建立。

在电磁暂态仿真程序(Electromagnetic Transient Program，EMTP)中，具有二阶精度且A-稳定的隐式梯形积分法被广泛用于对电路元件的微分方程差分化。然而，隐式梯形积分法不是L-稳定的，在面对因网络拓扑变化而引起非状态量的突变时，其将产生非原型的数值振荡。

为了抑制数值振荡现象，国内外研究人员提出了一系列的技术途径，其大致包括两大类：1)：附加阻尼元件，在EMTP中曾通过附加阻尼元件的方法来抑制数值振荡现象。2)：算法切换，如EMTP版本3.0中采用的临界阻尼调整法(Critical Damping Adjustment，CDA)。该算法在电磁暂态计算中，正常情况依旧采用隐式梯形法进行计算，仅在系统检测到扰动时，才将数值方法切换成半步长隐式欧拉法进行计算。该方法的主要问题在于有些突变现象难以被检测到，导致CDA方法仍然无法避免数值振荡问题。

发明内容

为解决上述技术问题，本发明提供一种基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，该方法能有效地解决隐式梯形积分方法所存在的数值振荡问题，且其数值计算的稳定性和精度比CDA方法都要高。

本发明采取的技术方案为：

基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，包括数值积分步骤，所述数值积分步骤采用A-稳定且无限稳定的3步4阶隐式泰勒级数法，进行时域数值积分计算，依次逐步求解各物理随时间的变化曲线。

基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，步骤如下：

步骤一：输入原始参数，建立电路各元件的微分方程，包括常微分方程和偏微分方程，对于用偏微分方程描述的元件应首先进行空间上离散将其转化为常微分方程，从而形成电磁暂态数值计算的统一形式的数学模型：x⁽¹⁾(t)＝f(x(t),t)+g(t)；

步骤二：电磁暂态数值计算初始化，设置仿真初始时刻t＝0.0s，积分步数n＝0；设定数值积分定步长h以及电磁暂态仿真计算的总时间t_f；设定系统中各状态变量的初值，即x(t＝0)＝x₀；确定输入电磁暂态数值计算的故障或操作；

步骤三：故障或操作判断，根据检测系统在时刻t的检测结果判断系统有无故障或操作：若无故障或操作，则直接转至步骤四；

若有故障或操作，则修改系数矩阵A和激励源g^(t)中相应位置的元素，并重新形成方程(1)；步骤四：数值积分，采用A-稳定且无限稳定的3步4阶隐式泰勒级数法，计算出状态变量在t＝t_n+1＝t_n+h处的值x_n+1；

步骤五：t＝t_n+1＝t_n+h；令n＝n+1；

步骤六：数值积分过程是否终止判断，

若t<t_f，则转至步骤三，继续下一时刻的数值积分；

若t≥t_f，则转至步骤七；

步骤七：电磁暂态数值仿真结果输出。

对一阶常微分方程的初值问题：

在所述步骤四中，用到的一种3步4阶隐式泰勒级数法计算格式为：

式中：h为时间积分步长；x_n-i≈x(t_n-i),i＝0,1,2,3是状态变量在t＝t_n-i时刻的近似值；

为状态变量在t＝t_n时刻的i阶导数的近似值；

求解具体步骤如下：

显然，利用方程(3)逐步积分求取状态变量在t＝t_n时刻的值x_n，需要已知前3步状态变量的值，也就是说需要启动电磁暂态计算的附加数值积分算法。为此，采用显式4阶泰勒级数法启动上述的数值积分过程，具体的计算公式如下：

第一步：启动数值积分即采用方程(4)计算出x₁和x₂：

第二步：根据方程(3)利用前3步的值逐步求取状态变量的值：

若f(t,x(t))是x(t)的线性函数，即f(t,x(t))＝Ax(t)；其中，A是定常系数矩阵。则x(t)的1阶导数为

x⁽¹⁾(t)＝Ax(t)+g(t) (7)；

由于在电力系统的元件如高压输电线路的激励源g(t)一般是正弦信号，故可引入新的变量并对状态变量的维数进行扩增，从而达到将方程(7)转化为齐次方程的形式的目的。设齐次化后的方程(7)如下：

上式(8)中，是扩维后的状态变量；B是扩维后的定常系数矩阵。

根据方程(7)，可以得到如下的递推公式：

根据方程(3)、(9)，可推导出的一般计算公式：

上式(10)中，I是与B同阶的单位矩阵。

若f(t,x(t))是x(t)的非线性函数，则采用牛顿迭代法求解式(2)得到状态变量在t＝t_n时刻的值x_n。

本发明一种基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，技术效果如下：1：该方法是A-稳定且无限稳定的数值方法，其对截断误差具有较快的衰减速率，可有效地抑制数值振荡，相对于隐式梯形积分法而言，其能够彻底避免数值振荡问题；该方法是计算精度为6阶，局部截断误差为O(h⁸)，其可以通过采用较大的积分步长而提高计算效率。与CDA方法相比较，基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法不仅可以完全避免数值振荡问题，而且计算精度和效率更高、数值稳定性更强。

2：与传统的、基于隐式梯形积分方法的电磁暂态数值计算方法相比较，基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法可有效地避免数值振荡问题；与CDA方法相比较，基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法不仅可完全避免数值振荡问题，而且具有更高的计算精度和更强的数值稳定性。

附图说明

图1：3步4阶隐式泰勒级数法的稳定域示意图。

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

图3a：隐式梯形法用于测试电路1的末端负载电压的数值计算结果图。

图3b：3步4阶隐式泰勒级数法用于测试电路1的末端负载电压的数值计算结果图。

图4：均匀工频长输电线空载合闸示意图。

图5：均匀工频长输电线Π型级联的集总等效电路模型图。

图6a：CDA方法用于长输电线路空载合闸电磁暂态数值计算结果图，计算步长h＝1.0μs。

图6b：3步4阶隐式泰勒级数法用于长输电线路空载合闸电磁暂态数值计算结果图，计算步长h＝80μs。

具体实施方式

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

以图4所示的单相工频均匀长输电线路在不同初相角情况下空载合闸的电磁暂态仿真计算为例，本发明的具体实施步骤归纳如下：

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

众所周知，用于描述图4所示高压输电线路电磁暂态过程的数学模型是电报方程。然而，电报方程是偏微分方程，需要首先将其转化成常微分方程才能采用数值方法进行电磁暂态计算。为实施的方便，直接采用文献(崔翔.无损耗传输线物理模拟的集总电路级联数目确定方法[J].中国电机工程学报,2017,37(9):2561-2570.)中的Π型级联的集总等效电路模型，如图5所示。对高压输电线路进行数学建模。

输入输电线路的分布参数R₀、L₀和C₀，线路全长为L，线路的边界条件，如首端的激励电压源e(t)及其内电阻。

空间离散化：取线路的区段数M＝50，离散后每一段线路的电阻r、电感l以及电容c依次如下：

r＝R₀L/M,l＝L₀L/M,c＝C₀L/M (11)

上式(11)中，R₀、L₀和C₀分别表示单位长度输电线路的电阻、电感和电容的值；L表示输电线路的全长；M表示空间离散的区间段数。

由图5中的等效电路，根据基尔霍夫电压和电流定律，易建立如下的一阶线性常微分方程组：

上式(12)中，i_m表示图5中第m个支路(离散区间)段上的电流；u_m表示图5中第m个空间离散节点处对地电压。r和l表示图5中第m个支路(离散区间)段上的电阻和电感；c为图5中第m个空间离散节点处对地电容值。

将式(12)整理成如下的矩阵形式：

上式(13)中：A∈R^{(2M+1)×(2M+1)}是定常稀疏矩阵；μ(t)是(2M+1)维稀疏列向量，它是电磁暂态计算的激励源；且有

上式(14)中，S为M+1维零矩阵；y(t)∈R^(2M+1)×1为待求状态变量；μ(t)∈R^(2M+1)×1为激励源；A∈R^{(2M+1)×(2M+1)}是与输电线路相关的定常稀疏系数矩阵。

L＝diag(l₁₁,l₂₂,…,l_MM),l_ii＝l,i∈(1,M)；

R＝diag(r₁₁,r₂₂,…,r_MM),r_ii＝r,i∈(1,M)；

C＝diag(c₁₁,c₂₂,…,c_(M+1)(M+1)),c_ii＝c,i∈(1,M+1)； (15)

上式(15)中，L、R和C均为对角矩阵。

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

设置仿真初始时刻t＝0.0s，积分步数n＝0；

设定数值积分定步长h＝80μs以及电磁暂态仿真计算的总时间t_f＝0.06s；

设定系统中各状态变量的初值，即y(t＝0)＝y₀；由于是空载线路，故y₀中的元素皆为0；

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

t≤0_-时，i₀(t)＝0；

t≥0₊时，

3)：故障或操作判断：

当t≤0_-，无任何操作，电磁暂态数值计算的数学模型不变，则直接转至步骤4)；

当t≥0₊时，开关合闸，此时需要对式(13)中的相关系数矩阵进行修改，具体情况可描述如下：

上式(18)中，y(t)∈R^(2M+1)×1为待求状态变量；为输电线路故障或操作后新的激励源；是与故障或操作后输电线路相关的、新的定常稀疏系数矩阵。

上式(19)中，

上式(20)中，为输电线路故障或操作后新的激励源；

为方便下一步采用3步4阶隐式泰勒级数法数值积分，考虑到y(t)对时间求高阶导数，故采用状态变量增维的方法从而将式(18)转化为齐次线性微分方程如下：

上式(21)中，是增维后的待求状态变量所组成的列向量；B∈R^{(2M+3)×(2M+3)}为状态变量增维后新的系数矩阵，且记：

上式(22)中，α(t)和δ(t)是新增加的2个状态变量，且有

上式(23)中，△y表示新增的2维列状态变量所组成的列向量，D是其相应的系数矩阵。

上式(24)中，E是式(18)增维后其系数矩阵中一个分块矩阵。

由式(21)可知，对时间求p阶导数的一般表达式如下：

上式(25)中，是增维后的待求状态变量所组成的列向量；B为状态变量增维后新的系数矩阵。

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

步骤四：数值积分：

采用A-稳定且无限稳定的3步4阶隐式泰勒级数法，计算出状态变量在t＝t_n+1＝t_n+h处的值具体情况可描述如下

第一步：启动数值积分，先计算出和

上式(26)～(27)中，h表示时间积分步长；是新的状态变量在t＝t_i时刻的近似值；为状态变量在t＝t_n时刻的i阶导数的近似值；

第二步：根据方程(3)利用前3步的值逐步求取状态变量的值：

根据方程(3)、(25)，可推导出的一般计算公式：

上式(28)中，h表示时间积分步长；其中I是与B同维的单位矩阵；是新的状态变量在t＝t_i时刻的近似值；

步骤五：t＝t_n+1＝t_n+h；令n＝n+1；

步骤六：数值积分过程是否终止判断(t_f为总的仿真时间)

若t<t_f，则转至步骤三，继续下一时刻的数值积分；

若t≥t_f，则转至步骤七；

步骤七：电磁暂态数值仿真结果输出。

该计算实例主要是输出空载合闸时线路末端电压曲线，即u(t)＝u_M+1(t)的变化曲线，具体如图6a和图6b所示。

在文献(张芳,仇雪芳,李传栋.电力系统中长期过程动态仿真的组合积分算法[J].电力自动化设备,2017,37(2):113-120.)中，通过对传统3步4阶隐式泰勒级数法的计算格式进行改进，并采用待定系数法构造出了一种3步4阶隐式泰勒级数法计算格式，即公式(3)。采用数值的方法画出了该3步4阶隐式泰勒级数法的稳定域如图1所示。图1中阴影部分以外的区域为本发明所采用的3步4阶隐式泰勒级数法的稳定区域。显然，算法的数值稳定域包括了整个开左半复平面，因而它是A-稳定的。此外，3步4阶隐式泰勒级数法也是无限稳定的。

概括起来，本发明所采用的3步4阶隐式泰勒级数法的计算精度为6阶，局部截断误差为O(h⁸)；在数值稳定性上，该方法具有A-稳定性和无限稳定性。因此，相对于A-稳定的2阶隐式梯形积分法而言，该算法对截断误差具有较快的衰减速率，从而可有效地抑制数值振荡。以上就是本发明的理论基础。

下面给出3步4阶隐式泰勒级数法不会产生数值振荡的具体实例。

图2是一无损传输线的仿真实例，传输线的分布参数已在图中标明，激励源为阶跃电压源e_i＝60ε(t)V，电源内阻r_s＝100Ω，负载为纯电阻r_L＝100Ω，负载电压记为u_L。图3a是利用隐式梯形法(计算步长h＝5.0ns)对该无损传输线进行电磁暂态计算的结果，产生了严重的数值振荡；图3b是利用3步4阶隐式泰勒级数法(计算步长h＝50ns)进行数值计算的结果。显然，从图3b中可以看出：当突变信号从首端传到末端时，3步4阶隐式泰勒级数法没有产生数值振荡。

本发明的重点在于采用了A-稳定且无限稳定的3步4阶隐式泰勒级数法，既有效地避免了数值振荡问题，又具有较高的计算精度和较强的数值稳定性。

Claims

1.基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，其特征在于：包括数值积分步骤，所述数值积分步骤采用A-稳定且无限稳定的3步4阶隐式泰勒级数法，进行时域数值积分计算，依次逐步求解各物理随时间的变化曲线。

2.基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，其特征在于：针对电磁暂态数值计算的数学模型，归结于一阶常微分方程的初值问题：

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mi>d</mi> <mrow> <mi>d</mi> <mi>t</mi> </mrow> </mfrac> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&equiv;</mo> <mover> <mi>x</mi> <mo>&CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>f</mi> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>)</mo> <mo>+</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>=</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>x</mi> <mn>0</mn> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>:</mo> </mrow>

式中：t表示时间变量；x(t)∈R^m×1为待求的m维向量函数，称之为状态变量；f(t,x(t))是关于时间t和状态变量x(t)的一维线性或非线性函数；g(t)是非齐次项，它是一个只与时间有关的稀疏列向量，通常称作激励源，x₀是状态变量在初始时刻的值；

在所述的数值积分算法中，用到的一种3步4阶隐式泰勒级数法计算格式为：

为状态变量在t＝t_n时刻的i阶导数的近似值；

求解具体步骤如下：

利用方程(2)逐步积分求取状态变量在t＝t_n时刻的值x_n，需要已知前3步状态变量的值，也就是说需要启动电磁暂态计算的附加数值积分算法；为此，采用显式4阶泰勒级数法启动上述的数值积分过程，具体的计算公式如下：

第一步：启动数值积分即采用方程(3)计算出x₁和x₂：

第二步：根据方程(2)利用前3步的值逐步求取状态变量的值：

若f(t,x(t))是x(t)的线性函数，即f(t,x(t))＝Ax(t)；其中，A是定常系数矩阵，则x(t)的1阶导数为：

x⁽¹⁾(t)＝Ax(t)+g(t) (6)；

由于在电力系统的元件，如高压输电线路的激励源g(t)一般是正弦信号，故可引入新的变量并对状态变量的维数进行扩增，从而达到将方程(6)转化为齐次方程的形式的目的，设齐次化后的方程(6)如下：

上式(7)中，是扩维后的状态变量；B是扩维后的定常系数矩阵；

根据方程(7)，可以得到如下的递推公式：

<mrow> <msup> <mover> <mi>x</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>p</mi> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>B</mi> <msup> <mover> <mi>x</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>p</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>...</mo> <mo>=</mo> <msup> <mi>B</mi> <mi>p</mi> </msup> <mover> <mi>x</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>p</mi> <mo>&GreaterEqual;</mo> <mn>1</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

根据方程(2)、(8)，可推导出的一般计算公式：

<mrow> <msub> <mover> <mi>x</mi> <mo>~</mo> </mover> <mi>n</mi> </msub> <mo>=</mo> <msup> <mi>T</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mfrac> <mn>3888</mn> <mn>3661</mn> </mfrac> <msub> <mover> <mi>x</mi> <mo>~</mo> </mover> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <mfrac> <mn>243</mn> <mn>3661</mn> </mfrac> <msub> <mover> <mi>x</mi> <mo>~</mo> </mover> <mrow> <mi>n</mi> <mo>-</mo> <mn>2</mn> </mrow> </msub> <mo>+</mo> <mfrac> <mn>16</mn> <mn>3661</mn> </mfrac> <msub> <mover> <mi>x</mi> <mo>~</mo> </mover> <mrow> <mi>n</mi> <mo>-</mo> <mn>3</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>,</mo> <mi>n</mi> <mo>&GreaterEqual;</mo> <mn>3</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

上式(9)中，I是与B同阶的单位矩阵，若f(t,x(t))是x(t)的非线性函数，则采用牛顿迭代法求解式(2)得到状态变量在t＝t_n时刻的值x_n。

3.基于3步4阶隐式泰勒级数法的电磁暂态数值计算方法，其特征在于步骤如下：

若有故障或操作，则修改系数矩阵A和激励源g(t)中相应位置的元素，并重新形成方程(1)；

步骤四：数值积分，采用A-稳定且无限稳定的3步4阶隐式泰勒级数法，计算出状态变量在t＝t_n+1＝t_n+h处的值x_n+1；

步骤五：t＝t_n+1＝t_n+h；令n＝n+1；

步骤六：数值积分过程是否终止判断，

若t<t_f，则转至步骤三，继续下一时刻的数值积分；

若t≥t_f，则转至步骤七；

步骤七：电磁暂态数值仿真结果输出。