CN100470995C - 电力系统特征值的分布式计算方法 - Google Patents

Abstract

电力系统特征值的分布式计算方法是属于电力系统分布式仿真技术领域，其特征在于，它包含带边界分区的互联电网切分方法和电力系统特征值的分布式计算方法。其中电力系统特征值的分布式计算方法包含了分布式求解发电机的等值导纳、分布式求解边界协调方程、分布式求解分区机械转矩的偏差量、分布式求解系统特征值。它可由各分区计算低频振荡模式相关的特征值，从而为研究大规模互联电力系统的低频振荡问题提供重要参考。在计算过程中仅需要各分区与边界分区交换边界节点状态量等少量数据，适用于电力系统的分布式环境，具有较好的实用性。

Description

电力系统特征值的分布式计算方法

技术领域

电力系统特征值的分布式计算方法是属于电力系统分布式仿真技术领域。

背景技术

在大规模互联电力系统中，发电机或发电机群之间常常发生转子的相对摇摆，并在缺乏阻尼时引起持续的振荡。振荡的频率范围一般在0.2~2.5Hz之间，故称为低频振荡，或机电振荡。近年来低频振荡在我国也时有发生，严重影响了电网间功率输送和安全稳定。在这种情况下，对低频振荡问题的研究备受关注。

我国大区互联电网的运行和管理具有“分级管理、分层控制、分布处理”的特点，在这种环境下，利用电力系统分布式计算技术，通过分解协调各调度中心的计算任务，来完成全网一体化的低频振荡分析工作，能更好的分析区域间振荡的模式。但是，迄今为止有关低频振荡的分布式分析方法的研究成果还比较少。一些学者在小干扰稳定分析领域以提高计算速度为目的开展了一系列并行算法的研究，然而由于并行计算和广域分布式计算的环境差别很大，并行算法并不能直接应用到电力系统分布式的计算和分析中。

本专利基于自激法的分布式特征值算法，该方法可由各分区计算低频振荡模式相关的特征值，从而为研究大规模互联电力系统的低频振荡问题提供重要参考。该方法在计算过程中仅需要各分区与边界分区交换边界节点状态量等少量数据，适用于电力系统的分布式环境，具有较好的实用性。

基于自激法的特征值计算方法

自激法是用于分析电力系统小干扰稳定的机电模式的特征值分析方法。它的基本思想来自线性系统的频域分析理论。在电力系统中，对于其中一台发电机，描述其动态特性的微分方程为：

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{.}{\mathrm{\ω}}=\frac{1}{2H}(\mathrm{\Δ}{T}_m-\mathrm{\Δ}{T}_e)\\ \mathrm{\Δ}\stackrel{.}{\mathrm{\δ}}={\mathrm{\ω}}_0\mathrm{\Δ\ω}\end{array}\right.---\left(1\right)$

其中ω为转子角速度，ω₀为系统基准角速度，δ为转子角，

为角速度偏差的导数，

角度偏差的导数(本文下文Δ都是指变量的偏差，变量上方加点都表示该变量的导数)。H为发电机的惯性常数，T_m为机械转矩，T_e为电磁转矩，ΔT_m为机械转矩的偏差，ΔT_e为电磁转矩的偏差，且ΔT_e＝K_SΔδ+K_DΔω，K_S为同步转矩系数，K_D为阻尼系数，K_S包含了整个网络对这台发电机的影响。

根据式(1)可以用拉普拉斯变换求得以ΔT_m为输入，Δω为输出的系统传递函数为

$G\left(s\right)=\frac{\mathrm{\Δ\ω}\left(s\right)}{\mathrm{\Δ}{T}_m\left(s\right)}=\frac{1}{2\mathrm{Hs}+K_D\left(s\right)+\frac{K_S\left(s\right)}{s}}---\left(2\right)$

式中的s是拉普拉斯变换的变换因子，在频域分析中当传递函数分母为零时，求出的s值即为对应系统的特征值。自激法的基本原理就是通过计算式(2)的极点来计算与该发电机机电模式相关的的特征值。该极点可以通过Newton迭代法求解，迭代形式为

$s_{k+1}=s_k-{\left[\frac{\mathrm{\Δ}{T}_m\left(s\right)}{4H_e}\right]}_{s=s_k}---\left(3\right)$

其中H_e为“等值惯性系数”，定义为

$H_e=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{H}_i{\left|\mathrm{\Δ}{\mathrm{\ω}}_i\right|}^2---\left(4\right)$

式中H_i和Δω_i分别为第i台发电机的惯性常数和转子角速度的变化量，n为系统中发电机的数目。当 $\left|\mathrm{\&Delta;}{T}_m\left(s\right)\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Delta;}{T}_m}$ (其中

位某个阈值，可取10^-6)，迭代收敛，此时的s_k就是所求的极点，即系统的特征值。求得系统特征值之后就可以进行低频振荡分析。

传统的自激法采用的是集中计算方式，这种方式对于计算能力的要求较高，计算速度相对较慢，而且集中式的自激法在迭代过程中的每一步都需要通过求解全系统的网络方程来解出机械转矩的扰动量ΔTm的值，我国电力系统的运行和管理具有“分级管理、分层控制、分布处理”的特点，各区域调度中心拥有管辖区域内电网的参数和动态数据，而上级调度中心拥有区域间联络线网络的参数和动态数据。同级调度中心之间所能共享的数据是有限的，因此整合数据的难度非常大。以上原因制约了传统的自激法的发展。在分布式的环境下，需要对传统的自激法进行改进以适应电力系统实际运行和管理的特点。本文以下章节将在带边界分区的互联电网切分方法的基础上，从自激法的基本原理出发，提出分布式的特征值计算方法以及左右特征向量的分布式计算方法，以适应电力系统分布式低频振荡分析的实际需求。

发明内容

本发明基于传统自激法进行特征值计算的一般步骤，结合我国电力系统分层分级管理的特点，提出了分布式的特征值的计算方法。该方法仅需要各分区与边界分区交换边界节点状态量等少量数据，即可分布式求解特征值。对IEEE39节点系统的实际计算表明，该方法得出了与传统方法一致的结果，通信次数也较少，适用于电力系统分布式环境的低频振荡分析，具有较好的实用性。

通常情况下，为了实现大规模互联电网分布式特征值计算，会采用包括负责协调计算的协调计算服务和多个调度中心内的分区特征值计算服务的分布式计算环境。协调计算服务和分区计算服务之间通过网络进行通信，交换必要的数据，完成一体化特征值计算。根据我国电力系统的实际情况，我们考虑如附图1的电力系统分布式环境。在该环境中，包括上下两级调度中心，上级调度中心与各区域调度中心之间通过广域的电力数据网相连接。区域调度中心掌握有区域内网络的详细信息和实时状态；上级调度中心掌握有区域间联络线以及直属厂站的的详细信息和实时状态。各区域调度中心和上级调度中心均需要部署实现了本发明所述算法的分析计算系统。本发明的特征在于，它具有以下实施步骤(具体流程如图2所示)：

步骤1 初始化(互联电网的切分)

互联电网的切分方法是实现分布式计算的基础。广域互联电网的分布式计算要求从电力系统的拓扑连接关系出发，由上级调度中心根据系统实际的运行和管理的情况对系统进行切分。本文采用文献(陈颖.电力网格中分布式计算方法研究[D].北京：清华大学，2006.04.)提出的基于功率平衡条件的带边界分区的互联电网切分方法。切分后的系统中，各区域调度中心拥有管辖区域内电网的参数和动态数据(下文也将区域调度中心称为分区)，而上级调度中心拥有区域间联络线网络的参数和动态数据(下文也将上级调度中心称为协调侧或者边界分区)。

以图3所示的电力系统为例，系统包括两个分区A₁、A₂，它们通过联络线1相连接。B₁代表分区1的边界节点，记它们的集合为

分区内除边界节点外的其他节点组成的网络，称为内部节点，记它们的集合为

分区1的调度中心掌握有

和

的详细数据。分区2与之类似。切分时，将联络线1及两端的虚拟节点

(分别对应B₁，B₂)单独看作是一个分区，称为“边界分区”。

运用上述切分方法后，第i个分区的边界节点满足下述关系

$\left\{\begin{array}{c}{u}_x^{B_i}+{\mathrm{ju}}_y^{B_i}=u_x^{{\tilde{B}}_i}+{\mathrm{ju}}_y^{{\tilde{B}}_i}\\ (i_x^{B_i}+{\mathrm{ji}}_y^{B_i})+(i_x^{{\tilde{B}}_i}+{\mathrm{ji}}_y^{{\tilde{B}}_i})=0\end{array}\right.---\left(5\right)$

其中u_x，u_y，i_x，i_y分别是同步坐标系下电压和电流在x、y坐标轴上的分量，i＝1，2。

步骤2 选定研究的发电机，求出其振荡频率在0.2~2.5Hz之间的特征值

这一特征值指的是电力系统线性化方程状态矩阵的特征值，下文统一简称特征值。自激法的关键是在每一步迭代过程中计算出机械转矩的扰动量(ΔT_m)的值。在应用“带边界分区的互联电网切分方法”划分网络后，使用分布式的自激法计算小干扰特征值。图4是用分布式方法计算电网特征值的流程图。

步骤201分布式计算机械转矩的扰动量ΔT_m

在每一步迭代过程中计算出机械转矩的扰动量ΔT_m的值是自激法的关键(公式(3))，本文按照图5所示的详细流程计算ΔT_m.

第一步 初始化分区信息

这一步的计算顺序如图6所示。假设上述流程由分区1发起。同时，根据电力系统实际运行和管理的情况，我们可以认为分区1发起自激法的计算时，选择的自激机也在本分区内。分区1选择s初值之后，通过边界分区转发给各分区。另外，因为电力系统是时变系统，所以分区1发起计算时，还需要传递时间参考坐标。各分区收到该时间坐标后，则统一选取该时刻的自己分区的网络导纳矩阵和线性化状态矩阵参与计算。初始化的流程如图6所示，图中的虚线箭头表示分区间数据的流向。

第二步 分区1需要根据当前迭代的特征值s_k计算第k台发电机的等值电流偏差量

ΔI_e(s)表示等值电流修正量，k是分区1内发电机编号，采用Prabha Kundur在“PowerSystem Stability and Control”一书提到的方法计算等值电流修正量

$\mathrm{\&Delta;}{I}_{e_k}\left(s_k\right)=c_1+\frac{c_2{\mathrm{\&omega;}}_0}{s_k}+C_r{(s_{k}I-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s_k})$

(本文公式中以黑体符号表示的变量为向量，而以普通符号表示的则为标量)

式中a₁₁、a₁₂、a_1r、ω₀、a_r1、a_r2、A_rr为下文公式(6)中状态变量 $\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]$ 的系数矩阵的元素、b₁、B_r为电压变化量Δv的系数，这些系数在网络参数确定的情况下为已知数。

第三步 各分区还需要根据当前迭代的特征值s_k分别计算本分区内各台发电机的等值导纳

第k台发电机(自激机)的等值导纳阵为：

$Y_{e_k}\left(s\right)=Y_{D_k}-C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r$

第i台发电机(i≠k，i＝1，2，…h，h为系统发电机的总台数)的等值导纳阵为：

$Y_{e_i}\left(s\right)=-C_i{(\mathrm{sI}-A_i)}^{-1}{B}_i+Y_{D_i}$

式中I是和A_rr维数相等的单位矩阵，A_i，C_i为公式(11)中Δx_i的系数矩阵，B_i，

为公式(11)中Δv的系数矩阵。

第四步 分布式求解边界协调方程

如图5按照前面的步骤分解协调计算母线电压向量、边界母线电流向量，图中用阴影部分标出了这部分迭代的过程。其中修正量的计算公式为：

对于分区i：

$\mathrm{\&Delta;}{v}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}(\mathrm{\&Delta;}{I}_e\left(s\right){|}_i^{\mathrm{In}}-Y_i^{\mathrm{InB}}\&CenterDot;\mathrm{\&Delta;}{v}_i^B)$

$\mathrm{\&Delta;}{I}_e\left(s\right){|}_i^B=Y_i^{\mathrm{BIn}}\mathrm{\&Delta;}{v}_i^{\mathrm{In}}+Y_i^{\mathrm{BB}}\mathrm{\&Delta;}{v}_i^B$

对于协调侧：

ΔI_e(s)|^B-Y^BΔv^B＝δ^B

式中Δv_i表示第i个分区的电压向量的修正值，上标In表示内部节点，上标B表示边界节点(下文统一使用这一识别原则)，

表示第i个分区所有内部节点的等值电流修正值，其元素由第i个分区内各节点的电流修正量组合而成，

表示第i个分区所有边界节点的等值电流修正值，ΔI_e(s)|^B表示协调侧所有边界节点的等值电流修正值。δ^B为计算ΔI_e(s)|^B-Y^BΔv^B的差值，该值用来判断边界电压、电流是否收敛。分区i的导纳阵(不计入边界分区线路)的导纳阵为：

m为分区内节点数目，Y_i，j为节点i，j(i，j＝1，2…m)之间的互导纳，该导纳阵的行、列序号对应节点序号，将该导纳阵按照内部节点、边界节点重新排列，可以得到：

$Y_i=\left[\begin{array}{cc}{Y}^{\mathrm{InIn}}& Y^{\mathrm{InB}}\\ Y^{\mathrm{BIn}}& Y^{\mathrm{BB}}\end{array}\right]$

和分别表示矩阵Y在分区i中对应内部节点、边界节点的矩阵以及两者之间的关联矩阵，直接取Y_i中相应元素即得。Y^B表示矩阵Y在边界分区内的矩阵。若‖δ^B‖>ε_δ，则计算不收敛，需要修正边界节点的电压、电流值。修正过程采用JFNG函数(文献“A Jacobian-Free Newton-GMRES(m)Method with Adaptive Preconditioner andItsApplication for Power Flow Calculations，IEEE TRANSACTIONS ON POWER SYSTEMS，VOL.21，NO.3，AUGUST 2006”已经公开该函数)，输入所有分区边界电压变化量Δv^B，得到所有分区边界电压变化量Δv^B的修正量，可以表示为：Δv^B＝Δv^B+Δ(Δv^B)，其中Δ(Δv^B)＝JFNG(Δv^B)，然后重新计算第四步。

若‖δ^B‖<ε_δ，则计算收敛。程序转入下面第五步。

以上公式属于本文原创，下面附上其推导过程。

发电机的完整的线性化模型方程可以记为

$\left[\begin{array}{c}\mathrm{\&Delta;}{\stackrel{.}{\mathrm{\&omega;}}}_k\\ \mathrm{\&Delta;}{\stackrel{.}{\mathrm{\&delta;}}}_k\\ \mathrm{\&Delta;}{\stackrel{.}{x}}_k\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{1r}\\ {\mathrm{\&omega;}}_0& 0& 0\\ a_{r1}& a_{r2}& A_{\mathrm{rr}}\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]+\left[\begin{array}{c}{b}_1\\ 0\\ B_r\end{array}\right]\mathrm{\&Delta;v}+\left[\begin{array}{c}\frac{1}{2H_k}\\ 0\\ 0\end{array}\right]\mathrm{\&Delta;}{T}_m---\left(6\right)$

$\mathrm{\&Delta;}{i}_k=\left[\begin{array}{ccc}{c}_1& c_2& C_r\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]-Y_{D_k}\mathrm{\&Delta;v}---\left(7\right)$

Δω_k、Δδ_k的定义同上文，Δi_k为发电机机端注入电流修正量，下标k表示该发电机在系统中的标号。Δv为网络母线电压向量修正量，Δx_k为模型方程中除了Δω_k和Δδ_k之外的所有状态变量的向量修正量(取决于所用的模型，当模型确定时Δx_k为一确定的向量)。c₁，c₂，C_r公式(7)状态变量 $\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]$ 的系数，

为电压修正量Δv的系数，这些系数在网络参数确定是为已知量(详细公式参见PRABHA KUNDUR.Power System Stability and Control.McGraw-HillCompanies，Inc.，1994.第12章公式12.8节)，对上式进行拉普拉斯变换并消去Δδ_k和Δx_k，可以得到

$\mathrm{\&Delta;}{T}_m=2H[s-a_{11}-\frac{a_{12}{\mathrm{\&omega;}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})]\mathrm{\&Delta;}{\mathrm{\&omega;}}_k-2H[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]\mathrm{\&Delta;v}---\left(8\right)$

$\mathrm{\&Delta;}{i}_k=\mathrm{\&Delta;}{I}_{e_k}\left(s\right)-Y_{e_k}\left(s\right)\mathrm{\&Delta;v}---\left(9\right)$

其中

$\mathrm{\&Delta;}{I}_{e_k}\left(s\right)=c_1+\frac{c_2{\mathrm{\&omega;}}_0}{s}+C_r{(sI-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})---\left(10\right)$

$Y_{e_k}\left(s\right)=Y_{D_k}-C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r$

式中各量的定义同上。

对于所研究的发电机，可以设Δω_k＝1.0p.u.(p.u.表示标幺值)。那么，在式(8)中，除Δv以外都是发电机k的本地变量，所以分解协调求解Δv，就成了实现分布式计算ΔT_m的关键。

对于除了发电机k以外的网络上的其他发电机，其线性化的状态方程也具有与式(8)类似的形式，式中的ΔT_m＝0。如果简记为下面的形式

$\mathrm{\&Delta;}{\stackrel{.}{x}}_i=A_i\mathrm{\&Delta;}{x}_i+B_i\mathrm{\&Delta;v}$

$\mathrm{\&Delta;}{i}_i=C_i\mathrm{\&Delta;}{x}_i-Y_{D_i}\mathrm{\&Delta;v}---\left(11\right)$

其中下标i表示该发电机在系统中的标号，i＝1，2…h，i≠k，h为系统发电机总台数，Δx_i为设备状态变量的变化量，Δi_i为设备向网络的注入电流的变化量，Δv为网络母线电压向量的变化量，A_i，C_i为Δx_i的系数矩阵，B_i，

为公式(11)中Δv的系数矩阵。那么根据拉普拉斯变换消去状态变量即可得到

$\mathrm{\&Delta;}{i}_i=[C_i{(\mathrm{sI}-A_i)}^{-1}{B}_i-Y_{D_i}]\mathrm{\&Delta;v}=-Y_{e_i}\left(s\right)\mathrm{\&Delta;v}---\left(12\right)$

全网的电压和电流向量满足下面的约束关系：

ΔI＝Y_NΔv                   (13)

其中ΔI为注入电流向量的变化量，Δv为母线电压向量的变化量，式中Y_N为全系统网络导纳矩阵，形式可以表示为：

h为全系统节点个数，Y_i，j为节点i，j(i，j＝1，2…h)之间的互导纳。

把式(9)和式(12)与式(13)合并，可得等效电流：

ΔI_e(s)＝(Y_N+Y_e(s))Δv                (14)

其中ΔI_e(s)中只有发电机k对应的位置具有非零元素Y_e(s)为

的对角阵，i＝1，2，…，h，当i＝k时，

由式(9)确定，否则由式(12)确定。

上式可以写成式(15)差值函数的形式以便于迭代求解

g(Δv)＝ΔI_e(s)-(Y_N+Y_e(s))Δv         (15)

记Y＝Y_N+Y_e(s)，按照步骤1的带边界分区的互联电网切分方法，如果系统分为边界分区和n个分区，那么对于分区i＝1，…，n，有

$\left[\begin{array}{c}{g}_i^{\mathrm{In}}\left(\mathrm{\&Delta;v}\right)\\ g_i^B\left(\mathrm{\&Delta;v}\right)\end{array}\right]=\left[\begin{array}{c}\mathrm{\&Delta;}{I}_{e_i}^{\mathrm{In}}\left(s\right)\\ \mathrm{\&Delta;}{I}_{e_i}^B\left(s\right)\end{array}\right]-\left[\begin{array}{cc}{Y}_i^{\mathrm{InIn}}& Y_i^{\mathrm{InB}}\\ Y_i^{\mathrm{BIn}}& Y_i^{\mathrm{BB}}\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{v}_i^{\mathrm{In}}\\ \mathrm{\&Delta;}{v}_i^B\end{array}\right]=0---\left(16\right)$

式中上标含In的量为分区内部节点的变量，上标含B的量为分区边界节点的变量。对于边界分区，有

$g^{\tilde{B}}\left(\mathrm{\&Delta;v}\right)=\mathrm{\&Delta;}{I}_e^{\tilde{B}}\left(s\right)-Y^B\mathrm{\&Delta;}{v}^{\tilde{B}}=0---\left(17\right)$

式中上标

的量为边界分区中的物理量，故

表示边界分区的等效电流修正值，

表示边界分区的电压修正值。

从式(16)和(17)容易得到 $\mathrm{\&Delta;}{v}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}(\mathrm{\&Delta;}{I}_{e_i}^{\mathrm{In}}\left(s\right)-Y_i^{\mathrm{InB}}\&CenterDot;\mathrm{\&Delta;}{v}_i^B)$ 以及 $\mathrm{\&Delta;}{I}_e^B\left(s\right)-Y^B\mathrm{\&Delta;}{v}^B=0$ (迭代过程中 $\mathrm{\&Delta;}{I}_e^B\left(s\right)-Y^B\mathrm{\&Delta;}{v}^B={\mathrm{\&delta;}}^B,$ δ^B为计算过程的差值)。

第五步 分区1计算ΔT_m

分区1计算ΔT_m

$\mathrm{\&Delta;}{T}_m=2H[s-a_{11}-\frac{a_{12}{\mathrm{\&omega;}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})]-2H[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]\mathrm{\&Delta;v}$

式中各变量的说明见上文。

步骤202 分区1判断ΔT_m收敛性

如果|ΔT_m|<10^-6，计算收敛，则内层迭代结束，此时求解出的s即为待求的特征值，特征值计算结束；

如果|ΔT_m|>10^-6，计算不收敛，，程序转向步骤203。

步骤203 计算全网的等值惯性常数H_e

计算过程如图7所示，各分区的 $H_{e_i}=\underset{j=1}{\overset{m_i}{\mathrm{\&Sigma;}}}{H}_j{\left|\mathrm{\&Delta;}{\mathrm{\&omega;}}_j\right|}^2,$ 全网的等值惯性常数 $H_e=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{H}_{e_i}$ 其中H_e为等值惯性系数，

为第i个分区的等值惯性系数，H_j和Δω_j分别为第j台发电机的惯性常数和转子角速度的变化量，n为分区的数目。当 $\left|\mathrm{\&Delta;}{T}_m\left(s\right)\right|<{\mathrm{\&epsiv;}}_{\mathrm{\&Delta;}{T}_m}$ (其中 ${\mathrm{\&epsiv;}}_{\mathrm{\&Delta;}{T}_m}>0$ 为设定的收敛判据，可以取 ${\mathrm{\&epsiv;}}_{\mathrm{\&Delta;}{T}_m}={10}^{-6}$ )时，迭代收敛，此时的s_k就是所求的极点，即系统的特征值。

步骤204 分区1计算特征值s

则此时的特征值可以由 $s_{k+1}=s_k-{\left[\frac{\mathrm{\&Delta;}{T}_m\left(s\right)}{4H_e}\right]}_{s=s_k}$ 算得。程序转到步骤(201)。

式中各量的定义同上。

本节利用IEEE39节点系统来测试上述的分布式算法。IEEE39节点系统包括了39条母线，46条线路，测试系统分为边界区域和3个区域，各区域包含的边界节点数目如表1所示。

表1 测试系统区域节点数目

Table 1 Test systems and its partitions

对上述系统进行分布式小干扰稳定特征值计算。选取参数为： ${\mathrm{\&epsiv;}}_{\mathrm{\&Delta;}{T}_m}=1e-6,{\mathrm{\&epsiv;}}_{\mathrm{\&delta;}}=1e-6.$ 所研究的发电机编号为39，属于区域1。初值s＝5.0j。迭代收敛的过程与传统串行计算的自激法收敛过程对比如表2所示。

表2 迭代收敛过程比较

Table 2 Iteration comparison between serialized and distributed methods

在分布式算法中，第1步外层迭代中内层的边界协调方程迭代通信了6次，随后，协调侧通过使用上一次的预处理矩阵进行计算，内层边界协调方程的迭代通信次数有所减少，最后外层迭代接近收敛时只需要5次通信。通信过程中所传递的数据量与边界区域的规模有关，在本测试系统中，边界区域节点数为8，因此通信过程中只需要传递8个节点的数据。

在电力系统分布式环境下，网络通信的延迟是分布式计算最主要的性能瓶颈。本文所提出的分布式的特征值算法，由于在迭代过程中需要各区域与边界区域多次交换数据，性能上可能不及传统的串行算法。但是，在测试算例中也可以看出，本文所提出的算法在只要求传递边界节点状态量等少量数据的条件下，仍具有较好的收敛性，通信次数也较少，一般情况下能满足实际电力系统在一个调度周期内完成计算的要求。由于算法在计算过程中不需要各区域交换区域内的详细网络数据，因而在我国电力系统分层分区的管理体制下具有实际的意义，可以作为在线的小干扰稳定分析算法的基础。分布式算法收敛得到的特征值与传统串行算法得到的特征值保持一致，证明了分布式算法的正确性。

再选择母线#31所连接的发电机作为研究对象，计算特征向量和相关因子。右特征向量在求解特征值的过程中即可计算，不需要额外的通信过程。在本算例中，左特征向量计算过程中边界协调方程通信了6次，传递了8个边界节点的状态量数据；相关因子计算过程中各区域与边界区域协调通信了2次，传递了10个相关因子结果。由于在电力系统分布式环境下，网络通信的延迟是分布式计算最主要的性能瓶颈，所以上述通信过程是影响分布式算法的计算效率的最主要因素。从上面统计的通信情况可以看出，本文所提出的算法只要求传递边界节点向量等少量数据，只需要少量的通信即可完成计算，一般情况下能满足实际电力系统在一个调度周期内完成计算的要求。由于算法在计算过程中不需要各区域交换区域内的详细网络数据，因而在我国电力系统分层分区的管理体制下具有实际的意义。

附图说明：

图1仿真的电力系统分布式环境示意图。图中表示一个上级调度中心和三个区域调度中心，它们之间通过电力网络相连。

图2应用自激法进行低频振荡分析的流程图。本文的工作是分布式求解电力系统特征值，该项工作是低频振荡分布式分析方法的核心步骤。

图3带边界分区的分区方法。图中两侧的区域是两个分区，中间的区域是协调侧。对于实际的系统而已，分区往往多于两个，而协调侧始终是一个。

图4分布式方法计算特征值的整体流程。

图5计算ΔT_m的流程。该图实际是图4中计算ΔT_m的具体实现。

图6计算特征值前，各分区和协调侧的初始化的流程。

图7计算H_e的流程。

图86节点3分区系统，用来展示具体的实施方式。

具体实施方式：

示例系统描述：

以一个6节点3分区系统为例说明计算的过程。系统包括2台发电机、1个负荷、6条母线、3条线路和3条变压器支路，如图8所示。系统参数(PSAT格式)：

节点的信息：

Bus.con＝[...

1   1.0     1   0    1    1；

2   1.0     1   0    2    1；

3   1.0     1   0    3    1；

4   1.0     1   0    1    1；

5   1.0     1   0    2    1；

6   1.0     1   0    3    1；]；

联络线信息：

Line.con＝[...

6      5      100    1.0     50     0     0      0.032      0.161

0.306      0      0      0       0     0；

5      4      100    1.0     50     0     0      0.01       0.085

0.176      0      0      0       0     0；

6      4      100    1.0     50     0     0      0.017      0.092

0.158      0      0      0       0     0；

2      5      100    1.0     50     0   1.025        0      0.0625

0      0      0      0       0      0；

3      6      100    1.0     50     0   1.025        0      0.0586

0      0      0      0       0      0；

1      4      100    1.0     50     0   1.025        0      0.0576

0      0      0      0       0      0]；

松弛节点信息：

SW.con＝[...

1     100   1.0   1.04       0      99   -99       1.1     0.9

0.8      1]；

Pv节点信息：

PV.con＝[...

2     100   1.0   1.63    1.025     99   -99       1.1     0.9

1；

3     100   1.0   0.85    1.025     99   -99       1.1     0.9

1]；

Pq节点信息：

PQ.con＝[...

6      100    1.0      0.9    0.3      1.2     0.8     0；

5      100    1.0      1      0.35     1.2     0.8     0；

4      100    1.0      1.25   0.5      1.2     0.80    ]；

同步   机信息：

Syn.   con＝[...

2      100    1.0      50       2        0.0521    0      0.8958   0.1198

0      6      0        0.8645   0.1969   0         0.535  0        12.8

0      0      0        1        1        0.002；

3      100    1.0      50       2        0.0742    0      1.3125   0.1813

0      5.89   0        1.2578   0.25     0         0.6    0        6.02

0      0      0        1        1        0.002；

]；

％ D

Syn.con(1，19)＝10.0；

Syn.con(2，19)＝10.0；

步骤1：

系统分为3个分区，分别为分区1、分区2和分区3。分区包括发电机1、母线2和母线5；分区2包括发电机2、母线3和母线6；分区3包括负荷、母线1和母线4。将各分区的边界节点：母线4、母线5和母线6作为边界分区的母线。

在下面的计算中，选定发电机2作为自激机，属于分区2。

初始迭代选取s＝10.0i.

步骤201：

第一步 初始化分区信息：

分区1导纳阵：

Y_InIn：

017.361

-17.3610

Y_InB：

0-17.361

17.3610

Y_BIn：

0-17.361

17.3610

Y_BB：

017.361

-17.3610

分区2导纳阵：

Y_InIn：

10.651+0.228i7.419-0.456i

-15.905+0.114i-10.651-0.228i

Y_InB：

0-16

160

Y_BIn：

0-16

160

Y_BB：

016

-160

分区3导纳阵：

Y_InIn：

2.35220.019

-20.422-2.352

Y_InB：

0-17.065

17.0650

Y_BIn：

0-17.065

17.0650

Y_BB：

017.065

-17.0650

边界分区：

3.307379  21.94778   -1.36519   -11.6041    -1.94219    -10.5107

-21.9478  3.307379   11.6041    -1.36519    10.51068    -1.94219

-1.36519  -11.6041   2.552792   17.33823    -1.1876     -5.97513

11.6041   -1.36519   -17.3382   2.552792    5.975135    -1.1876

-1.94219  -10.5107   -1.1876    -5.97513    3.129796    16.25382

10.51068  -1.94219   5.975135   -1.1876     -16.2538    3.129796

第二步根据当前迭代的特征值s_k计算等值电流

当s＝10.0i时， $\mathrm{\&Delta;}{I}_{e_k}\left(s\right)=\left[\begin{array}{c}-185.45i\\ -34.03i\end{array}\right]$

第三步 各分区还需要根据当前迭代的特征值s_k分别计算本分区内各台发电机的

$Y_{e_1}\left(s\right)=\left[\begin{array}{cc}11.0558-0.3401i& -9.3908+0.6808i\\ 0.2975-0.1699i& -11.0558+0.3401i\end{array}\right]$

$Y_{e_2}\left(s\right)=\left[\begin{array}{cc}2.3517& 2.9539\\ -3.3568& -2.3517\end{array}\right]$

第四步 分布式求解边界协调方程

迭代收敛后Δv＝

0.3231-0.7833i

-0.2991-8.4029i

0.3367-7.1604i

-0.2974+6.0561i

0.2933+5.1534i

-0.2834-22.5695i

0.3231-0.7833i

-0.2991-8.4029i

0.3296-2.8327i

-0.3001-2.4659i

0.3120+1.0625i

-0.2920-14.8989i

第五步计算ΔT_m

ΔT_m＝11.9908+8.2126i

步骤202 分区1判断ΔT_m收敛性

当前|ΔT_m|>10^-6，不收敛

步骤203 计算全网的惯性常数H_e

分区1： $H_{e_1}=10.0$

分区2： $H_{e_2}=11.234$

最后总的H_e＝21.234

步骤204 分区1计算特征值s

s＝-0.9959+9.3179i

以上数据是第一步迭代的结果，由于计算尚未收敛，下面需要进行多次迭代直到收敛，限于篇幅不予记录。当迭代收敛时，可以给出的最终的特征值，该特征值即可交付分析人员使用。

Claims (1)

1、电力系统分布式特征值的计算方法，其特征在于，该方法依次含有以下步骤：步骤(1)大区互联电网基于功率平衡条件的带边界分区的互联电网切分方法按下述方式把互联电网切分为：多个各称为分区的区域调度中心以及一个称为协调侧或边界分区的上级调度中心，分区管辖区域内电网的参数和动态数据，边界分区管辖各区域间联络线l所组成的网络的参数和动态数据，以联络线l连接的相邻分区分别命名为A₁和A₂，B₁、B₂分别代表分区A₁和A₂的边界节点，

分别代表各自边界节点的集合，分区内除边界节点以外的其它节点称为内部节点，相应的内部节点集合记为

把联络线l及其两端的虚拟节点

构成的区域称为边界分区，

分别对应于各分区边界节点B₁、B₂；

所述各分区的边界节点满足下述关系：

$\left\{\begin{array}{c}{u}_x^{B_i}+{\mathrm{ju}}_y^{B_i}=u_x^{{\tilde{B}}_i}+{\mathrm{ju}}_y^{{\tilde{B}}_i}\\ (i_x^{B_i}+{\mathrm{ji}}_y^{B_i})+(i_x^{{\tilde{B}}_i}+{\mathrm{ji}}_y^{{\tilde{B}}_i})=0\end{array}\right.,$

其中u_x，u_y，i_x，i_y分别是同步坐标系下电压和电流在x、y坐标轴上的分量，i＝1，2；

步骤(2)分区A₁的计算按以下步骤：选定研究的发电机，求出该发电机的振荡频率在0.2~2.5Hz之间时，电力系统线性化方程的状态矩阵的特征值s；

步骤(2.1)计算机械转矩的偏差ΔT_m，其步骤如下：

步骤(2.1.1)初始化各分区信息：

步骤(2.1.1.1)分区A₁确定t时刻下自己的网络导纳阵和线性化状态矩阵，求出特征值s的初始估计值S₀；

步骤(2.1.1.2)分区A₁把值S₀和时间坐标t通过边界分区转发给各分区；

步骤(2.1.1.3)各分区在收到该时间坐标后，统一选取该时刻t的自己分区的网络导纳阵和线性化状态矩阵参与下述计算：

步骤(2.1.2)分区A₁用当前迭代的特征值S_k计算本分区内第k台发电机的等值电流偏差量，步骤如下：

步骤(2.1.2.1)把第k台发电机的完整的线性化模型方程记为：

$\left[\begin{array}{c}\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_k\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{\mathrm{\&delta;}}}_k\\ \mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_k\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{1r}\\ {\mathrm{\&omega;}}_0& 0& 0\\ a_{r1}& a_{r2}& A_{\mathrm{rr}}\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]+\left[\begin{array}{c}{b}_1\\ 0\\ B_r\end{array}\right]\mathrm{\&Delta;v}+\left[\begin{array}{c}\frac{1}{{2H}_k}\\ 0\\ 0\\ \end{array}\right]\mathrm{\&Delta;}{T}_m$

${\mathrm{\&Delta;i}}_k=\left[\begin{array}{ccc}{c}_1& c_2& C_r\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]-Y_{D_k}\mathrm{\&Delta;v}$

其中：Δi_k为发电机机端注入电流向量偏差量，下标k表示该发电机在系统中的标号，

Δv为网络母线电压向量偏差量，

Δx_k为模型方程中除了Δω_k和Δδ_k之外的所有状态变量的向量偏差量，为设定值， $\left[\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&omega;}}_k\\ \mathrm{\&Delta;}{\mathrm{\&delta;}}_k\\ \mathrm{\&Delta;}{x}_k\end{array}\right]$ 及Δv的系数a₁₁、a₁₂、a_1r、ω₀、a_r1、a_r2、A_rr、c₁、c₂、C_r以及

b₁、B_r在互联电网网络方程确定时为已知量，

H_k为第k台发电机的惯性常数，为已知量，ω₀为发电机基准角速度，为设定值，ω_k为发电机k的转子角速度，δ_k为转子角，Δω_k、Δδ_k分别是角速度偏差和转子角偏差；

对于除了发电机k以外的网络上的其他发电机，其线性化的状态方程也为下面的形式 $\mathrm{\&Delta;}{\stackrel{\&CenterDot;}{x}}_i=A_i\mathrm{\&Delta;}{x}_i+B_i\mathrm{\&Delta;v}$

${\mathrm{\&Delta;i}}_i=C_i{\mathrm{\&Delta;x}}_i-Y_{D_i}\mathrm{\&Delta;v},$

其中下标i表示该发电机在系统中的标号，i＝1，2…h，i≠k，h为系统发电机总台数，Δx_i为设备状态变量的变化量，A_i，C_i为Δx_i的系数矩阵，B_i，

为Δv的系数矩阵；

步骤(2.1.2.2)按下式计算

${\mathrm{\&Delta;I}}_{e_k}\left(s_k\right)=c_1+\frac{c_2{\mathrm{\&omega;}}_0}{s_k}+C_r{(s_{k}I-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s_k}),$

其中，I为和A_rr维数相等的单位矩阵；

步骤(2.1.2.3)包括A₁在内的各分区根据当前迭代特征值S_k计算各分区内各台发电机的等值导纳：

作为自激机的第k台发电机的等值导纳： $Y_{e_k}\left(s\right)=Y_{D_k}-C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r,$

第i台发电机的等值导纳： $Y_{e_i}\left(s\right)=-C_i{(\mathrm{sI}-A_i)}^{-1}{B}_i+Y_{D_i},$

i≠k，i＝1，2，…h，h为系统发电机的总台数；

步骤(2.1.3)包括分区A₁在内的各分区求解边界收敛方程，其步骤如下：

步骤(2.1.3.1)边界分区设定Δv^B的值，然后通过 ${\mathrm{\&Delta;v}}^{{\tilde{B}}_i}={\mathrm{\&Delta;v}}^{B_i}$ 送达各分区；

步骤(2.1.3.2)各分区依次计算以下各量，得到δ^B：

第i分区内部，节点的电压偏差

${\mathrm{\&Delta;v}}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}({\mathrm{\&Delta;I}}_e\left(s\right){|}_i^{\mathrm{In}}-Y_i^{\mathrm{InB}}\&CenterDot;{\mathrm{\&Delta;v}}_i^B),$

表示第i个分区所有内部节点的等值电流修正值，

第i分区边界节点的等效电流偏差

${\mathrm{\&Delta;I}}_e\left(s\right){|}_i^B=Y_i^{\mathrm{BIn}}{\mathrm{\&Delta;v}}_i^{\mathrm{In}}+Y_i^{\mathrm{BB}}{\mathrm{\&Delta;v}}_i^B,$

各分区全部边界节点的等效电流偏差ΔI_e(s)|^B：

ΔI_e(s)|^B-Y^BΔv^B＝δ^B，

其中δ^B用来判断边界电压、电流是否收敛，是一个边界协调参数，ΔI_e(s)|^B为组合各分区的边界节点等效电流得到的向量；

分区i的导纳阵(不计入边界分区线路)的导纳阵为：

m为分区内节点数目，Y_i，j为节点i，j(i，j＝1，2…m)之间的互导纳，该导纳阵的行、列序号对应节点序号，将该导纳阵按照内部节点、边界节点重新排列，可以得到：

$Y_i=\left[\begin{array}{cc}{Y}^{\mathrm{InIn}}& Y^{\mathrm{InB}}\\ Y^{\mathrm{BIn}}& Y^{\mathrm{BB}}\end{array}\right]$

和

分别表示分区i中对应内部节点、边界节点的矩阵以及两者之间的关联矩阵，直接取Y_i中相应元素即得，Y^B表示边界分区的导纳阵；

步骤(2.1.3.3)判断‖δ^B‖是否小于ε_δ，ε_δ＝10^-6：

若‖δ^B‖>ε_δ，则执行步骤(2.1.3.4)，

若‖δ^B‖<ε_δ，则执行步骤(2.1.4)；

步骤(2.1.3.4)令

令Δv^B＝Δv^B+Δ(Δv^B)，其中Δ(Δv^B)＝JFNG(Δv^B)，JFNG是一个函数，再重复执行步骤(2.1.3.1)～步骤(2.1.3.3)直到‖δ^B‖<ε_δ为止；

步骤(2.1.4)分区A₁计算第k台发电机的ΔT_m

${\mathrm{\&Delta;T}}_m=2H[s-a_{11}-\frac{a_{12}{\mathrm{\&omega;}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\&omega;}}_0}{s})]-2H[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]\mathrm{\&Delta;v},$

H为发电机的惯性常数；

步骤(2.2)判断ΔT_m的收敛性：

如果|ΔT_m|<10^-6，计算收敛，此时求解出的Sk即为待求的特征值，

如果|ΔT_m|>10^-6，则不收敛，转步骤(2.3)；

步骤(2.3)按下式计算全网的等效惯性常数H_e：

$H_e=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{H}_{e_i},$

为各分区的等效惯性常数， $H_{e_i}=\underset{j=1}{\overset{m_i}{\mathrm{\&Sigma;}}}{H}_j{\left|{\mathrm{\&Delta;\&omega;}}_j\right|}^2,$

H_j和Δω_j分别为第j台发电机的惯性常数和转子角速度的变化量，n为分区的数目，m为系统总的发电机台数；

步骤(2.4)分区A₁计算迭代用的特征值：

$s_{k+1}=s_k-{\left[\frac{{\mathrm{\&Delta;T}}_m\left(s\right)}{{4H}_e}\right]}_{s=s_k}$

用该特征值重新计算步骤(2.1)～步骤(2.2)。

Images