基于増维精细积分的高压输电线路电磁暂态数值计算方法

摘要

一种基于増维精细积分的高压输电线路电磁暂态数值计算方法，其特征在于包括以下步骤：将高压输电线路的输入或激励正弦信号改写为状态变量，通过状态变量增维，将描述高压输电线路电磁暂态计算的非齐次常微分方程转换为齐次形式的常微分方程，在此基础上采用精细积分法对齐次常微分方程组进行数值积分计算。本发明继承或保持了经典精细积分方法所具有的高精度这一优点，因而可以有效避免传统电磁暂态数值计算方法所存在的数值振荡问题；与已有的电磁暂态数值积分方法相比较，基于増维精细积分方法的高压输电线路电磁暂态数值计算方法具有更高的计算效率。

说明书

技术领域

本发明涉及一种电力系统电磁暂态数值计算方法，具体是涉及一种基于増维精细积 分方法的高压输电线路电磁暂态数值计算方法。

背景技术

高压输电线路是电力系统的重要元件之一。利用数值计算方法分析高压输电线路的 电磁暂态过程，可解决电力系统诸多领域的问题：可对输电线路故障进行仿真，得到短 路电流以及谐波分量等；可用于输电线路故障测距算法的研究以及输电线路空载合闸过 电压计算研究等。

概括起来，高压输电线路的电磁暂态数值计算方法大致包括频域法和时域法这2大 类方法。频域类方法通常是先利用快速傅立叶变换或拉普拉斯变换求解输电线的频率响 应，然后利用相应的反变换以得到其时域解。频域类方法概念比较简单，但涉及的计算 极为繁琐，而且此类方法的求解过程主要是基于线性叠加原理，因而很难处理非线性问 题。

时域法的基本思路是采用时间离散或空间离散方法将偏微分方程转化为常微分方 程，并在此基础上利用各类数值积分方法求解各状态变量的值。目前使用较多的时域类 方法主要包括：时域有限差分法、基于隐式梯形积分方法和隐式欧拉法相结合的临界阻 尼调整法(CriticalDampingAdjustment，CDA)等。

时域有限差分法的主要思想是把高压输电线路的电报方程进行时间、空间同时离散， 用差分方程代替一阶偏微分方程，再通过求解线性方程组以计算出各状态变量的值。时 域有限差分法算法简单，但该方法中时间步长的选取受稳定性条件的限制，当空间步长 很小时亦须选取很小的时间步长，从而导致计算效率低下。

目前，商用化的EMTP(ElectromagneticTransientProgram，EMTP)程序中主要是采用 临界阻尼调整法(CDA法)来求解高压输电线路的电磁暂态计算问题。CDA方法计算过程 简单，从理论上讲可以有效避免传统的电磁暂态计算方法中的数值振荡问题，但对一些 无法准确判断的突变现象，使用该方法进行计算仍将产生“虚假的”数值振荡现象。

除上述数值计算方法以外，还有研究人员将精细积分方法应用到高压输电线路的电 磁暂态计算中。精细积分方法是钟万勰院士针对结构动力学方程的求解问题所提出的一 类新的数值计算方法，这种方法在求解齐次常微分方程时具有明显的优势：具有极高的 计算精度；计算过程简单、计算效率高。然而，对于高压输电线路，离散后的常微分方 程组通常是非齐次的，若直接采用精细积分法进行求解则将涉及矩阵求逆问题，而这种 矩阵求逆运算不仅计算量大，有时甚至是不稳定的。

发明内容

本发明所要解决的技术问题，就是提供一种基于増维精细积分的高压输电线路电磁 暂态数值计算方法，其能根据高压输电线路的输入或激励是正弦信号这一特点，将空间 离散后所形成的非齐次常微分方程组通过适当增维转换成齐次常微分方程组，并采用精 细积分法来求解齐次常微分方程组，可以精确、高效地求解出高压输电线路的电磁暂态 响应。

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

一种基于増维精细积分的高压输电线路电磁暂态数值计算方法，其特征在于包括以 下步骤：将高压输电线路的输入或激励正弦信号改写为状态变量，通过状态变量增维， 将描述高压输电线路电磁暂态计算的非齐次常微分方程转换为齐次形式的常微分方程， 在此基础上采用精细积分法对齐次常微分方程组进行数值积分计算。

具体步骤包括：

1)将高压输电线路进行空间离散化，建立相应的电磁暂态数值计算的数学模型，即：

$\stackrel{·}{x}=Ax+w\left(t\right)---\left(1\right);$

式(1)中，x是高压输电线路电磁暂态计算中的状态变量；A是与高压输电线路参数

相关的定常系数矩阵；w(t)即是电磁暂态计算的输入信号或激励信号，也就是所谓的

非齐次项，它是一个稀疏列向量；

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

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

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

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

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

3)故障或操作判断

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

若无故障或操作，则直接转至步骤4)；

若有故障或操作，则修改系数矩阵A以及激励信号w(t)并重新形成方程(1)；

同时，若涉及状态变量的突变，则还需依据具体情况修改相应的状态变量值x_n(t)；

4)常微分方程齐次化

依据时刻t判断w(t)是否为0：

若w(t)＝0，则无需齐次化，直接转至步骤5)；

若w(t)≠0，则利用状态变量増维方法将非齐次常微分方程即式(1)转换成以下形 式的齐次常微分方程：

$\stackrel{·}{\stackrel{~}{x}}=H\tilde{x}---\left(2\right);$

上式(2)中：H是与原系数矩阵A等有关的一个新的常系数矩阵；

上式(3)中，Δx即是増维状态变量；

5)数值积分

采用经典的精细积分方法对齐次常微分方程(2)进行求解，计算出所有状态变量的值；

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

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

若t＜T，则返回步骤3)：

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

8)数值计算结果输出。

所述步骤4)常微分方程齐次化是本发明的核心部分，其具体的实施步骤如下：

第一步：依据w(t)为正弦信号这一特点，引入新的变量并对状态变量进行増维

具体过程可叙述如下：

设非齐次项w(t)＝[w_i(t)]^T中含有M个非零元素，将这些不为零的元素依据其先后顺 序记为w_j(t)，j∈(1，M)；

由于高压输电线路电磁暂态数值计算的输入信号为正弦信号，因此，用以下通用表 达式来描述w_j(t)：

上式(4)中：ω＝2π×50；α_j为常量；为相位角，是已知量。

对应上述非零元素，引入如下变量：

则有：

$\left(\begin{array}{c}{\stackrel{·}{w}}_j\left(t\right)\\ {\stackrel{·}{y}}_j\left(t\right)\end{array}\right)=\left(\begin{array}{cc}0& \omega \\ -\omega & 0\end{array}\right)\left(\begin{array}{c}{w}_j\left(t\right)\\ y_j\left(t\right)\end{array}\right),j\in (1,M)---\left(6\right);$

因此，对应w(t)中每一个非零元素，定义以下増维状态变量：

${\mathrm{\Delta x}}_j=\left(\begin{array}{c}{w}_j\left(t\right)\\ y_j\left(t\right)\end{array}\right),j\in (1,M)---\left(7\right);$

$\Delta x={\left(\begin{array}{ccc}{\mathrm{\Delta x}}_1^T& ...& {\mathrm{\Delta x}}_M^T\end{array}\right)}^T\in R^{2M\times 1}---\left(8\right);$

第二步：非齐次常微分方程增维齐次化

具体过程可叙述如下：

设原状态变量x的维数为s，

定义

$d_j=\left(\begin{array}{cc}0& \omega \\ -\omega & 0\end{array}\right)\in R^{2\times 2},j\in (1,M)---\left(9\right);$

D＝diag(d_j)，j∈(1，M)；D∈R^2M×2M(10)；

$E=\left(\begin{array}{c}{e}_1\\ .\\ .\\ .\\ e_s\end{array}\right)\in R^{s\times \left(2M\right)};e_i\in R^{1\times 2M},i\in (1,s)---\left(11\right);$

上式(11)中：若w(t)中第i个元素为零，则e_i＝0；若w(t)中第i个元素不为零，则e_i是 一个2M维的稀疏行向量，该向量的第2j-1个元素为1，其余元素全部为0；

基于上述定义，将非齐次方程即式(1)转换为形如式(2)的齐次方程，即：

显然，H是一个常系数矩阵；

第三步：计算增维状态变量的初值

第四步：置齐次化结果标志

完成上述齐次化后，置w(t)＝0。

所述步骤5)数值积分中，利用精细积分法求解各状态变量的过程如下：

第一步：计算指数矩阵T

利用精细积分法求解齐次常微分方程，需求解如下的指数矩阵：

T＝e^Hh(14)；

指数矩阵的求解过程叙述如下。

将积分步长进一步细分为2^m段，则每一小段的时间间隔为：

Δt＝h/(2^m)(15)；

依据泰勒级数展开，首先计算如下初始矩阵：

$T_h^{\left(0\right)}=H\Delta t+{\left(H\Delta t\right)}^2/2!+...+{\left(H\Delta t\right)}^K/K!---\left(16\right);$

上式(16)中：K为泰勒多项式的阶数；在上述基础上，循环求解：

$T_h^{\left(i\right)}=2T_h^{(i-1)}+T_h^{(i-1)}{T}_h^{(i-1)},i\in (1,m)---\left(17\right);$

则式(14)中的指数矩阵为：

$T\approx I+T_h^{\left(m\right)}---\left(18\right);$

第二步：利用指数矩阵T计算各状态变量在t＝t_n+1时刻的值

利用上式(18)的结果，求解各状态变量在t＝t_n+1时刻的值：

${\tilde{x}}_{n+1}=\tilde{x}(t=t_{n+1})=T{\tilde{x}}_n={\tilde{x}}_n+T_h^{\left(m\right)}{\tilde{x}}_n---\left(19\right).$

上述精细积分方法具有极高的计算精度，因而能有效避免“虚假的”数值振荡现象。 此外，指数矩阵T只与H和步长h有关：若系数矩阵H和步长h不变，则只需一次性计 算出指数矩阵T即可；换言之，只有当发生故障或操作后才需要重新计算指数矩阵T。 因此，上述精细积分方法的计算效率高于传统的数值积分方法。

本发明的重点就在于：利用高压输电线路的输入信号或激励信号为正弦信号这一特 点，通过状态变量的适当增维，将非齐次常微分方程组转化为齐次常微分方程组，并采 用精细积分法来求解齐次常微分方程组，可以精确、高效地完成高压输电线路的电磁暂 态计算。

本发明的优点和积极效果是：基于増维精细积分方法的高压输电线路电磁暂态数值 计算方法，继承或保持了经典精细积分方法所具有的高精度这一优点，因而可以有效避 免传统电磁暂态数值计算方法所存在的数值振荡问题；与已有的电磁暂态数值积分方法 相比较，基于増维精细积分方法的高压输电线路电磁暂态数值计算方法具有更高的计算 效率。

附图说明

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

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

图3：高压均匀长输电线路空载合闸的电磁暂态计算结果(初相角)；

图4：高压均匀长输电线路空载合闸的电磁暂态计算结果(初相角)；

图5：高压输电线路电磁暂态计算结果三维曲线图(初相角)；

图6：CDA方法计算高压输电线路空载合闸的电磁暂态结果(初相角)；

图7：CDA方法计算高压输电线路空载合闸的电磁暂态结果(初相角)；

图8：相同步长下两种不同方法计算结果的对比曲线(初相角)；

图9：不同步长下两种不同方法计算结果的对比曲线(初相角)；

具体实施方式

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

以如图1所示的高压工频均匀长输电线路在不同初相角情况下进行空载合闸的电磁 暂态计算为例，本发明的具体实施步骤可以概述如下：

1)将高压输电线路进行空间离散化，建立相应的电磁暂态数值计算的数学模型

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

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

$\left(\begin{array}{c}l\frac{{\mathrm{di}}_k\left(t\right)}{dt}=l{\stackrel{·}{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{·}{v}}_k\left(t\right)=i_{k-1}\left(t\right)-i_k\left(t\right)\end{array}\right),k\in (1,N)---\left(20\right);$

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

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

上式(22)中：

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

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

I＝[i₁(t)…i_N(t)]^T，U＝[v₁(t)…v_N(t)v_N+1(t)]^T(27)； 显然，方程(22)最终可以写成形如方程(1)(即)的标准形式,其中，

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

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

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

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

输入电磁暂态数值计算的故障或操作：t＝0.0s时，线路空载合闸。

3)故障或操作判断

当t＝0.0s时，线路空载合闸，为此，需要对方程(1)中的系数矩 阵A以及激励信号w(t)进行修改。修改后的系数矩阵A以及激励信号w(t)的表达式分别 为：

上式(31)中，

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

当t＞0.0s时，无任何故障或操作发生，直接进入下一步骤；

4)常微分方程齐次化

当t＝0时，w(t)≠0，因此，需利用状态变量増维方法将非齐次常微分方程(1)齐次化。

具体步骤如下：

第一步：依据w(t)为正弦信号这一特点，引入新的变量并对状态变量进行増维

式(32)中仅含有1个非零元素，其值为：

上式(34)中：为空载合闸的初相角，是已知量。

对应上式(34)，引入如下变量：

则有：

$\left(\begin{array}{c}\stackrel{·}{w}\left(t\right)\\ \stackrel{·}{y}\left(t\right)\end{array}\right)=\left(\begin{array}{cc}0& \omega \\ -\omega & 0\end{array}\right)\left(\begin{array}{c}w\left(t\right)\\ y\left(t\right)\end{array}\right)---\left(36\right);$

因此，对应w(t)中的非零元素，可以定义以下増维状态变量：

$\Delta x=\left(\begin{array}{c}w\left(t\right)\\ y\left(t\right)\end{array}\right)---\left(37\right);$

第二步：非齐次常微分方程增维齐次化

在前述过程的基础上，很容易形成如下齐次常微分方程组：

上式(38)中：

$D=\left(\begin{array}{cc}0& \omega \\ -\omega & 0\end{array}\right)---\left(39\right);$

$E=\left(\begin{array}{c}{e}_1\\ .\\ .\\ .\\ e_{2N+1}\end{array}\right)\in R^{(2N+1)\times 2}---\left(40\right);$

上式(40)中：除e_N+1＝[10]外，其余各行向量均为0。

第三步：计算增维状态变量的初值

第四步：置齐次化结果标志

完成上述齐次化后，置w(t)＝0；

当t＞0时，w(t)＝0，因此，无需齐次化，直接转至步骤5)；

5)精细积分法求解齐次常微分方程组

第一步：计算指数矩阵

取m＝15，Δt＝h/2^m；

取K＝4，计算$T_h^{\left(0\right)}=H\Delta t+{\left(H\Delta t\right)}^2/2+{\left(H\Delta t\right)}^3/6+{\left(H\Delta t\right)}^4/24;$

循环求解：$T_h^i=2T_h^{(i-1)}+T_h^{(i-1)}{T}_h^{(i-1)},i\in (1,m);$

则指数矩阵为：

第二步：求解状态变量在t＝t_n+1时刻的值

${\tilde{x}}_{n+1}=T{\tilde{x}}_n={\tilde{x}}_n+T_h^{\left(m\right)}{\tilde{x}}_n---\left(42\right);$

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

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

若t＜T，则返回步骤5)中第二步，继续计算下一时刻沿线各点处电压、电流值；

若t＞T，则转至步骤8)。

8)数值计算结果输出：

①绘制初相角时输电线路末端电压随时间变化的曲线(图3)；

②绘制初相角时输电线路末端电压随时间变化的曲线(图4)；

③绘制初相角时输电线路各点处电压随时间变化的三维曲线图(图5)。

为进一步说明本发明方法的优越性，将增维精细积分方法的计算结果与CDA方法的 计算结果进行对比。

图6为初相角时，使用CDA方法(步长h＝1μs)计算高压输电线路末端电压的 电磁暂态结果，此时未出现数值振荡现象(与图3结果一致)；图7是初相角时，使 用CDA方法(步长h＝1μs)计算高压输电线路末端电压的电磁暂态结果，此时出现轻微的 数值振荡现象。对比图7和图4可以看出：与CDA方法相比较，增维精细积分方法抑制 数值振荡问题的效果更好。

图8是在采用相同步长(步长h＝1μs)情况下，增维精细积分方法与CDA方法两者计 算结果的差值曲线从图8可以看出：两种方法的计算结果十分接近。计算过程 及结果证实：在采用相同步长的情况下，使用增维精细积分方法所需的CPU时间较少， 其计算效率高于CDA方法。

图9是在分别采用不同的计算步长的情况下，增维精细积分方法(步长h＝10μs)与 CDA方法(步长h＝1μs)两者计算结果的差值曲线从图9可以看出：即使以10 倍于CDA方法的步长进行计算，增维精细积分方法仍能得到与CDA方法十分接近的结 果。因次，增维精细积分方法可以通过增大积分步长来进一步提高计算效率。