CN101051749A - 电力系统低频振荡的分布式分析方法 - 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值即为对应系统的特征值。

运用传统自激法进行电力系统低频振荡分析的一般步骤是：

第1步，上级调度中心选定研究的发电机，求出其振荡频率在0.2~2.5Hz之间的特征值；

第2步，上级调度中心计算该特征值对应的特征向量；

第3步，上级调度中心根据特征向量，计算该特征值与各发电机的状态变量Δω的相关因子；

第4步，上级调度中心根据相关因子，判断特征值与哪一台或几台发电机强相关，从而找出其它参与振荡的机组。

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

发明内容

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

本发明正是在传统自激法的步骤的基础上，研究出第2步、第3步和第4步的分布式计算方法和实用的分布式计算程序，在计算过程中仅需要各分区与边界分区交换边界节点状态量等少量数据，适合电力系统分布式的环境。

本文提出的分布式的低频振荡分析方法按照如下步骤进行：

第1步，分区调度中心选定研究的发电机，求出其振荡频率在0.2~2.5Hz之间的特征值；

第2步，分区调度中心计算该特征值对应的特征向量；

第3步，分区调度中心根据特征向量，计算该特征值与各发电机的状态变量Δω的相关因子；

第4步，分区调度中心根据相关因子，判断特征值与哪一台或几台发电机强相关，从而找出其它参与振荡的机组。

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

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

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

以图3所示的电力系统为例，系统包括两个分区A₁、A₂，它们通过联络线1相连接。B₁代表分区1的边界节点，记它们的集合为A₁ ^B；分区内除边界节点外的其他节点组成的网络，称为内部节点，记它们的集合为A₁ ^In。分区1的调度中心掌握有A₁ ^B和A₁ ^In的详细数据。分区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(3\right)$

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

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

步骤3计算分区1中该特征值对应的左、右特征向量

经过步骤2的计算，可以使特征值s收敛到系统的特征值，此时即可用下面的方法计算对应的右特征向量。

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

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

${\mathrm{\Δi}}_k=\left[\begin{array}{ccc}{c}_1& c_2& C_r\end{array}\right]\left[\begin{array}{c}{\mathrm{\Δ\ω}}_k\\ {\mathrm{\Δ\δ}}_k\\ {\mathrm{\Δx}}_k\end{array}\right]-Y_{D_k}\mathrm{\Δv}---\left(5\right)$

为角速度偏差的导数，

角度偏差的导数，下标k表示自激机在系统中的标号。

及Δv的系数α₁₁、α₁₂、a_1r、ω₀、a_r1、a_r2、A_rr、c₁、c₂、C_r以及Y_Dk、b₁、B_r在互联电网网络方程确定时为已知量，Δi_k为发电机机端注入电流变化量，下标k表示该发电机在系统中的标号。Δv为网络母线电压向量变化量，Δx_k为模型方程中除了Δω_k和Δδ_k之外的所有状态变量的向量变化量(取决于所用的模型，当模型确定时Δx_k为一确定的向量)。H_k为第k台发电机的惯性常数，为已知量。对上式进行拉普拉斯变换并消去Δδ_k和Δx_k，可以得到

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

步骤301计算右特征向量

右特征向量的计算流程如图4所示。

分布式自激法迭代收敛时，ΔT_m(s)≈0，右特征向量的第i个分量(i＝1，…，n)可以由式(7)求得：

${\mathrm{\Δ\ω}}_i{(s-a_{11}-\frac{a_{12}{\mathrm{\ω}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+a_{r2}))}^{-1}\·(b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r)\mathrm{\Δv}---\left(7\right)$

式中s和Δv在步骤2计算特征值的过程之中得到，其他变量定义同上文。网络参数b₁和B_r仅在所描述的发电机相连的母线的位置上有非零元素，计算Δω_k时，仅需要该母线的电压向量数据，而不需要整个全网的电压向量。对于系统中的其它的发电机的Δω，也可以用上式计算。所以，右特征向量的计算可以由各分区独立完成，即实现分布式的计算。

步骤302计算左特征向量

左特征向量的计算过程如图5所示。原系统的左特征向量的取值等于其对偶系统右特征向量。对偶系统是电力系统分析的一个基本工具。

第一步构建分区1的对偶系统

分区1的原系统可以写成以下形式：

$\mathrm{\Δ}\stackrel{\·}{x}=\mathrm{A\Δx}+\mathrm{B\Δv}$

                Δ_i＝CΔx-Y_DΔv      (8)

其中Δx为全系统的状态变量的变化量，Δi为网络注入电流向量的变化量，Δv为网络母线电压向量的变化量，A，B，C和Y_D分别为上式两个方程中Δx、Δv的系数矩阵。

那么可以相应地定义原系统的对偶系统的动态方程：

$\mathrm{\Δ}\stackrel{\·}{y}=A^T\mathrm{\Δy}+C^T\mathrm{\Δu}$

$\mathrm{\Δj}=B^T\mathrm{\Δy}-Y_D^T\mathrm{\Δu}---\left(9\right)$

其中Δy为对偶系统的状态变量，Δu为对偶系统的节点电压偏差量，Δj为对偶系统的节点注入电流的偏差量，Δy、Δu、Δj的取值范围同原系统的相应物理量。A^T、C^T、B^T和Y_D ^T在式(9)中的含义同式(8)中的A，B，C和Y_D，其取值为A，B，C和Y_D的转置(上标T为转置符号)。

同时根据原系统的电压和电流向量的约束方程

                        Δi＝Y_NΔv           (10)

式中Y_N为全系统网络导纳矩阵，形式可以表示为：

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

同时有：

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

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

相应地定义对偶系统的网络方程

$\mathrm{\Δj}=Y_N^T\mathrm{\Δu}---\left(11\right)$

其中Δy、Δu、Δj定义同上。

第二步求解对偶系统的右特征向量

左特征向量的计算过程如图5所示。左特征向量可以通过求解其对偶系统右特征向量来实现。

Step1初始化分区各参数

发电机k的线性化状态方程的代数方程如式(4)所示。

定义对偶系统中自激机的等值电流变化量(k为发电机序号)为：

${\mathrm{\ΔJ}}_{e_k}\left(s\right)=-2H{[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]}^T---\left(12\right)$

定义对偶系统中自激机的等效导纳阵(k为发电机序号)为：

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

对偶系统中，第k台发电机的注入电流的变化量定义为：

${\mathrm{\Δj}}_k={\mathrm{\ΔJ}}_{e_k}\left(s\right)-Y_{e_k}^T\left(s\right)\mathrm{\Δu}---\left(14\right)$

式(12)、(13)、(14)中其他变量的定义同上文。

对于除了发电机k以外的发电机，在对偶系统中的等效导纳阵可以表示为：

$Y_{e_j}^T\left(s\right)=-{[C_j{(\mathrm{sI}-A_j)}^{-1}{B}_j-Y_{D_j}]}^T,j\≠k---\left(15\right)$

其中j表示表示该发电机在系统中的标号。

将式(13)和式(15)代入式(11)可以得到对偶系统所有节点的等效电流：

${\mathrm{\ΔJ}}_e\left(s\right)=(Y_N^T+Y_e^T\left(s\right))\mathrm{\Δu}---\left(16\right)$

其中ΔJ_e(s)只有发电机k对应的位置具有非零元素ΔJ_ek(s)， $Y_e^T\left(s\right)=\mathrm{diag}\left\{Y_{e_j}^T\left(s\right)\right\}$ (Y_ej ^T(s)的对角矩阵)，其中j为全系统的动态设备的标号。

Step2分布式求解边界协调方程

${\mathrm{\Δu}}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}({\mathrm{\ΔJ}}_e\left(s\right){|}_i^{\mathrm{In}}-Y_i^{\mathrm{InB}}\·\mathrm{\Δ}{u}_i^B)$

${\mathrm{\ΔJ}}_e\left(s\right){|}_i^B=Y_i^{\mathrm{BIn}}{\mathrm{\Δu}}_i^{\mathrm{In}}+Y_i^{\mathrm{BB}}{\mathrm{\Δu}}_i^B$

式中Δu_i表示第i个分区的电压向量的变化值，上标In表示内部节点，上标B表示边界节点，ΔJ_e(s)|_i ^In表示第i个分区所有内部节点的等值电流变化值，其元素由第i个分区内各节点的电流变化量组合而成，ΔJ_e(s)|_i ^B表示第i个分区所有边界节点的等值电流变化值，ΔJ_e(s)|^B表示协调侧所有边界节点的等值电流变化值。推导过程如下：

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

$g\left(\mathrm{\Δv}\right)={\mathrm{\ΔJ}}_e\left(s\right)-(Y_N^T+Y_e^T\left(s\right))\mathrm{\Δu}---\left(17\right)$

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

$\left[\begin{array}{c}{g}_i^{\mathrm{In}}\left(\mathrm{\Δv}\right)\\ g_i^B\left(\mathrm{\Δv}\right)\end{array}\right]=\left[\begin{array}{c}{\mathrm{\ΔJ}}_e\left(s\right){|}_i^{\mathrm{In}}\\ {\mathrm{\ΔJ}}_e\left(s\right){|}_i^B\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{\Δu}}_i^{\mathrm{In}}\\ {\mathrm{\Δu}}_i^B\end{array}\right]=0---\left(18\right)$

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

$g^{\tilde{B}}\left(\mathrm{\Δv}\right)={\mathrm{\ΔJ}}_e\left(s\right){|}^B-Y^B{\mathrm{\Δu}}^{\tilde{B}}=0---\left(19\right)$

式中上标

的量为边界分区中的物理量，ΔJ_e(s)|^B表示边界分区的等效电流变化值，

表示边界分区的电压变化值。

式(19)是收敛时的情况，在迭代过程中满足

$g^{\tilde{B}}\left(\mathrm{\Δv}\right)={\mathrm{\ΔJ}}_e\left(s\right){|}^B-Y^B{\mathrm{\Δu}}^{\tilde{B}}={\mathrm{\δ}}^B---\left(20\right)$

其中δ^B为计算的差值。

从式(18)和(19)容易得到各子区域分解计算本区域的Δu_i ^In：

${\mathrm{\Δu}}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}({\mathrm{\ΔJ}}_e\left(s\right){|}_i^{\mathrm{In}}-Y_i^{\mathrm{InB}}\·\mathrm{\Δ}{u}_i^B)---\left(21\right)$

然后计算边界节点等效电流变化向ΔJ_e(s)|^B：

${\mathrm{\ΔJ}}_e\left(s\right){|}^B=Y_i^{\mathrm{BIn}}{\mathrm{\Δu}}_i^{\mathrm{In}}+Y_i^{\mathrm{BB}}{\mathrm{\Δu}}_i^B---\left(22\right)$

Step3判断收敛性

对式(20)计算的δ^B进行判断：

若‖δ^B‖＞ε_δ(ε_δ＝10^-6为判断阈值)，则计算不收敛，需要修正边界节点的电压变化量。修正过程采用JFNG函数(文献“A Jacobian-Free Newton-GMRES(m)Method with AdaptivePreconditioner and ItsApplication for Power Flow Calculations，IEEE TRANSACTIONS ONPOWER SYSTEMS，VOL.21，NO.3，AUGUST 2006”已经公开该函数)，输入所有边界节点电压变化量Δu^B，得到所有分区边界电压变化量Δu^B的修正量，可以表示为：Δu^B＝Δu^B+Δ(Δu^B)，其中Δ(Δu^B)＝JFNG(Δu^B)，然后重新计算step3。

若‖δ^B‖＜ε_δ，则计算收敛。程序转入下面step4。

Ste4计算对偶系统的右特征向量

通过上面的分解协调计算，各分区与边界分区交换边界节点变量数据，解得Δu。然后类似求解原系统右特征向量计算的方法，对偶系统的右特征向量的表达式为：

${\mathrm{\Δ\ψ}}_k{=\left\{{[s-a_{11}-\frac{a_{12}{\mathrm{\ω}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\ω}}_0}{s})]}^T\right\}}^{-1}\·{[c_1+\frac{c_2{\mathrm{\ω}}_0}{s}+C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\ω}}_0}{s})}^T\mathrm{\Δu}---\left(23\right)$

式中i＝1，…，t，t为分区个数。式中各量定义同上文。该特征向量即可作为原系统的左特征向量。

步骤4根据原系统左、右特征向量，计算该特征值与各发电机的状态变量Δω的相关因子

如果记Δω和Δψ分别为原系统右特征向量和左特征向量，其形式为：

$\mathrm{\Δ\ω}=\left[\begin{array}{c}{\mathrm{\Δ\ω}}_1\\ {\mathrm{\Δ\ω}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\Δ\ω}}_n\end{array}\right],$ $\mathrm{\Δ\ψ}={\left[\begin{array}{c}{\mathrm{\Δ\ψ}}_1\\ {\mathrm{\Δ\ψ}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\Δ\ψ}}_n\end{array}\right]}^T---\left(24\right)$

相关因子的定义为

                   p_k＝Δω_kΔψ_k        (25)

式(25)定义的相关因子表明了第k台发电机与所求得的特征值的相对的相关程度，|p_k|越大，则相关程度越高。找出|p_k|最大的一台或几台发电机，就可以判断它们与该振荡模式强相关，从而为分析人士提供参考，也为安装抑制低频振荡的PSS或线性最优励磁控制器提供重要的参考。

步骤5根据相关因子，判断特征值与哪一台或几台发电机强相关，从而找出其它参与振荡的机组

在步骤4计算左右特征向量和相关因子p_k后，选出|p_k|最大的一台或多台发电机的分布式计算的步骤为(流程见图6所示)：

步骤501，分区i＝1，…，n：对于本分区内的动态设备k，根据其对应的左右特征向量Δψ_k和Δω_k，计算相关因子p_k，并将结果传递给边界分区；

步骤502，边界分区：将所有的相关因子按照|p_k|排序，选出最大的一个或多个相关因子，将其编号传递给需要的分区。

得到左右特征向量之后，各分区就可以用式(25)计算各发电机的p_k，然后由边界分区选出|p_k|最大的一台或多台发电机，流程如图6所示，图中的虚线箭头表示数据流向。

在图6中，各分区将分区内发电机的p_k及对应的编号发送给边界分区的排序服务；边界分区再将排序的结果返回给各分区。根据实际需要，分区也可以只发送一部分发电机的p_k及对应的编号。分区收到边界分区返回的排序结果之后，就判断出了本次流程计算出的特征值与哪些发电机强相关，并可将结果输出给调度员或分析人员。

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

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

Table 1 Test systems and its partitions

系统名称	区域1	区域2	区域3	边界区域
					IEEE39	15	12	12	8

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

                   表2迭代收敛过程比较

Table 2Iteration comparison between serialized and distributed methods

迭代步骤n	传统算法sn	分布式算法
				外层迭代sn	边界协调方程迭代通信次数
		1	-0.0897+6.0744i			-0.0897+6.0744i	6
2	0.1182+4.6237i			-0.1182+4.6237i	6
		3	-0.1052+4.3468i			-0.1052+4.3468i	5
4	0.0983+4.3302i			0.0983+4.3302i	5
		5	-0.0977+4.3299i			0.0977+4.3299i	5
6	-0.0977+4.3299i			-0.0977+4.3299i	5

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

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

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

最后计算的结果如表3所示，系统具有特征值λ＝0.0791+7.3418i，对应了振荡频率为1.2Hz的低频振荡模式。分布式的算法与传统串行分析方法的分析结果相比，在误差容许范围内，各区域的左右特征向量及相关因子等结果均对应相等。分布式计算的相关因子的结果显示，位于区域1的母线#30和位于区域2的母线#31所连的发电机对于该特征值具有较大的相关因子，表明了该特征值对应的振荡模式与这两台发电机有关。串行计算的分析结论与之相符。以上结果表明，本文所提出的小干扰稳定特征特征向量与相关因子的分布式算法与传统串行计算相比具有一致的结果，证明了本文的分布式分析方法的正确性和有效性。

                            表3分析结果

                     Table 3 Computing Results

发电机	#31
				传统算法	分布式算法
	特征值	0.0791+7.3418i	0.0791+7.3418i
					母线	母线
	左特征向量	#39  -0.1371-0.0003i#31  0.6934+0.0000i#32  -0.2444+0.0002i#33  -0.0342+0.0008i#34  0.0183-0.0006i#35  0.1198-0.0026i#36  -0.0700+0.0016i#37  0.0414+0.0036i#38  -00040-0.0004i#30  -0.3491-00155i	区域1	#39  -0.1371-0.0003i#37  0.0414+0.0036i#38  -0.0040-0.0004i#30  -0.3491-0.0155i
区域2					#31  0.6934+0.0000i#32  -0.2444+0.0002i
			区域3	#33  -0.0342+0.0008i#34  0.0183-0.0006i#35  0.1198-0.0026i#36  -0.0700+0.0016i
右特征向量					#39  -0.0118+0.0000i#31  1.0000-0.0000i#32  -0.2861-0.0002i#33  -0.0755-0.0017i#34  0.0353+0.0011i#35  0.2779+0.0063i#36  -0.2075-0.0049i#37  0.0617-0.0073i#38  -0.0066+0.0003i#30  -0.3407+0.0154i	区域1	#39  -0.0118+0.0000i#37  0.0617-0.0073i#38  -0.0066+00003i#30  -0.3407+0.0154i
	区域2	#31  1.0000-0.0000i#32  -0.2861-0.0002i
			区域3	#33  -0.0755-0.0017i#34  0.0353+0.0011i#35  0.2779+0.0063i#36  -0.2075-0.0049i
	相关因子|pk|	#31  0.6934#30  0.1192#32  0.0699#35  0.0333#36  0.0145#33  0.0026#37  0.0026#39  0.0016#34  0.0006#38  0.0000				区域1	#30  0.1192#37  0.0026#39  0.0016#38  0.0000
区域2			#31  0.6934#32  0.0699
				区域3	#35  0.0333#36  0.0145#33  0.0026#34  0.0006

附图说明：

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

图2应用自激法进行低频振荡分析的流程图。本文的工作是分布式进行低频振荡的分析。

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

图4计算右特征向量的流程。

图5计算左特征向量的流程。

图6相关因子计算流程。

图76节点3分区系统

具体实施方式：

示例系统描述：

以一个6节点3分区系统为例说明计算的过程。系统包括2台发电机、1个负荷、6条母线、3条线路和3条变压器支路，如图7所示。系统参数(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.8    0]；

同步机信息：

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.228i 7.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计算等值电流ΔI_ek(s)

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

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

$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

步骤301计算右特征向量

-0.4489-0.0018i

1.0555-0.0961i

步骤302计算左特征向量

${\mathrm{\ΔJ}}_{e_2}\left(s\right)=$

$-0.1799$

$0.9806$

Δu＝

0.0465+0.0016i

-0.0045+0.0017i

-0.0288+0.0015i

-0.0432+0.0018i

0.1176+0.0015i

0.0366+0.0015i

0.0465+0.0016i

-0.0045+0.0017i

0.0180+0.0016i

-0.0224+0.0017i

0.0755+0.0015i

0.0170+0.0016i

Δψ＝

[-0.9428+0.0038i  1.0555+0.0961i]

步骤4根据特征向量，计算该特征值与各发电机的状态变量Δω的相关因子

p₁＝

0.4233

p₂＝

1.1234

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

Claims (1)

1、电力系统低频振荡的分布式分析方法，其特征在于，该方法是在电力系统的各分区调度中心的计算机上依次按以下步骤实现的；

步骤(1).初始化

向计算机输入分区1中所选定的研究发电机在其振荡频率在0.2HZ～2.5HZ之间的特征值s；

步骤(2).计算分区1中该特征值对应的左、右特征向量；

步骤(2.1).按照下式计算右特征向量的第i个分量Δω_i，i＝1，2，…，n；

${\mathrm{\Δ\ω}}_i={(s-{\mathrm{\α}}_{11}-\frac{{\mathrm{\α}}_{12}{\mathrm{\ω}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+a_{r2}))}^{-1}\·(b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r)\mathrm{\Δv},$ 其中：ω为转子角速度，Δω为角速度的偏差，α₁₁、α₁₂、a_1r、ω₀、a_r1、a_r2、A_rr、c₁、c₂、C_r以及Y_Dk、b₁、B_r在电力系统参数确定时为已知量，其中，ω₀为系统设定的发电机基准角速度，δ为转子角，Δδ分别是转子角偏差，k为分区1中该选定研究的发电机的序号，Δx_k为发电机完整的线性化模型方程中除了Δω_k和Δδ_k之外的所有状态变量的向量偏差量，当发电机的线性化模型确定时，Δx_k为一确定的向量，I是和A_rr维数相等的单位矩阵；

步骤(2.2).按照以下步骤计算左特征向量Δψ：

步骤(2.2.1)构造分区1的对偶系统：

设定分区1的原发电机全系统为：

$\mathrm{\Δ}\stackrel{\·}{x}=\mathrm{A\Δx}+\mathrm{B\Δv}$

                                 Δi＝CΔx-Y_DΔv

其中：Δx为所选全系统的状态变量的变化量，Δi为网络注入电流向量的变化量，Δy为网络母线电压向量的变化量，A，B，C和Y_D分别为Δx、Δv的系数矩阵，在电力网络参数确定时为已知值；

则所述发电机全系统的对偶系统的动态方程为：

$\mathrm{\Δ}\stackrel{\·}{y}=A^T\mathrm{\Δy}+C^T\mathrm{\Δu}$

$\mathrm{\Δj}=B^T\mathrm{\Δy}+Y_D^T\mathrm{\Δu}$

其中：Δy为对偶系统的状态变量，Δu为对偶系统的节点电压偏差量，Δj为对偶系统的节点注入电流的偏差量，A^T、C^T、B^T和Y_D ^T为A，B，C和Y_D的转置矩阵；

所述对偶系统的全系统网络方程为：

$\mathrm{\Δj}=Y_N^T\mathrm{\Δu}$

其中：Y_N为全系统网络导纳矩阵，形式可以表示为：

h为全系统节点个数，Y_i，j为节点i，j， $\left\{\begin{array}{c}i=\mathrm{1,2}\·\·\·h\\ j=\mathrm{1,2}\·\·\·h\end{array}\right.,$ 之间的互导纳，T为转置符号，同时满足：

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

其中：Y_i为第i个分区的导纳矩阵，B表示分区的边界节点，In表示分区i的内部节点，Y_i ^InIn表示分区i中对应内部节点间的关联导纳矩阵，Y_i ^InB表示内部节点和边界节点的之间的关联导纳矩阵；

步骤(2.2.2)按照下列步骤求解该对偶系统的右特征向量，其值等于原系统的左特征向量；

步骤(2.2.2.1)初始化分区的参数

对偶系系统中，作为自激机的发电机k的等效电流变化量ΔJ_ek(s)为：

$\mathrm{\Δ}{J}_{e_k}\left(s\right)=-2H{[b_1+a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}{B}_r]}^T$

对偶系统所有的节点的等效电流ΔJ_e(s)为：

$\mathrm{\Δ}{J}_e\left(s\right)=(Y_N^T+Y_e^T\left(s\right))\mathrm{\Δu}$

其中：Y_e ^T(s)为Y_ej ^T(s)组成的角矩阵，Y_ej ^T(s)为对偶系统中除了发电机k以外的其余发电机j的等效导纳矩阵：

$Y_{e_j}^T\left(s\right)=-{[C_j{(\mathrm{sI}-A_j)}^{-1}{B}_j-Y_{D_j}]}^T,j\≠k,$

自激机的等效导纳矩阵Y_ek ^T(s)为：

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

Δu与发电机k的注入电流变化量Δj_k有如下的关系：

$\mathrm{\Δ}{j}_k=\mathrm{\Δ}{J}_{e_k}\left(s\right)-Y_{e_k}^T\left(s\right)\mathrm{\Δu};$

步骤(2.2.2.2)按下式计算Δu_i ^In，Δu_i ^In为分区i内部节点电压向量的变化值；

$\mathrm{\Δ}{u}_i^{\mathrm{In}}={\left(Y_i^{\mathrm{InIn}}\right)}^{-1}(\mathrm{\Δ}{J}_e\left(s\right){|}_i^{\mathrm{In}}-Y_i^{\mathrm{InB}}\·\mathrm{\Δ}{u}_i^B)$

相应地边界节点的等效电流变化量ΔJ_e(s)|^B：

$\mathrm{\Δ}{J}_e\left(s\right){|}^B=Y_i^{\mathrm{BIn}}\mathrm{\Δ}{u}_i^{\mathrm{In}}+Y_i^{\mathrm{BB}}{\mathrm{\Δu}}_i^B$

步骤(2.2.2.3)判断收敛性

在计算过程中定义差值函数

$g^{\tilde{B}}\left(\mathrm{\Δv}\right)=\mathrm{\Δ}{J}_e\left(s\right){|}^B-Y^B\mathrm{\Δ}{u}^{\tilde{B}}={\mathrm{\δ}}^B,$

若‖δ^B‖＞ε_δ，ε_δ＝10^-6为判断阈值，则计算不收敛，需要修正边界节点的电压变化量，修正过程采用JFNG函数，输入所有边界节点电压变化量Δu^B，得到所有分区边界电压变化量Δu^B的修正量，表示为：Δu^B＝Δu^B+Δ(Δu^B)，其中Δ(Δu^B)＝JFNG(Δu^B)，然后重新计算步骤(2.2.1)到步骤(2.2.2.3)；

若‖δ^B‖＜ε_δ，则计算收敛，此时的Δu为系统的解电压变化量，程序转入下面；

步骤(2.2.2.4)按照下式计算对偶系统的右特征向量，其值等于原系统的左特征向量Δψ_k，

${\mathrm{\Δ\ψ}}_k={\left\{{[s-{\mathrm{\α}}_{11}-\frac{{\mathrm{\α}}_{12}{\mathrm{\ω}}_0}{s}-a_{1r}{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\ω}}_0}{s})]}^T\right\}}^{-1}\·{[c_1+\frac{c_2{\mathrm{\ω}}_0}{s}+C_r{(\mathrm{sI}-A_{\mathrm{rr}})}^{-1}(a_{r1}+\frac{a_{r2}{\mathrm{\ω}}_0}{s})]}^T\mathrm{\Δu}$ 式中i＝1，…，t，t为分区个数，该特征向量等于原系统的左特征向量；

步骤(3).根据相关因子，判断特征值与哪一台或几台发电机强相关，从而找出其它参与振荡的机组，

步骤(3.1).各分区按照下式计算该特征值与各发电机的状态变量Δω的相关因子p_k，并将结果传递给边界分区；

                               p_k＝Δω_kΔψ_k

其中： $\mathrm{\Δ\ω}=\left[\begin{array}{c}{\mathrm{\Δ\ω}}_1\\ {\mathrm{\Δ\ω}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\Δ\ω}}_n\end{array}\right],\mathrm{\Δ\ψ}={\left[\begin{array}{c}{\mathrm{\Δ\ψ}}_1\\ \mathrm{\Δ}{\mathrm{\ψ}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\Δ\ψ}}_n\end{array}\right]}^T$ 为原系统的右特征向量和左特征向量，|p_k|越大，则相关程度越高。找出|p_k|最大的一台或几台发电机，就可以判断它们与该振荡模式强相关；

步骤(3.2).各分区将所有的相关因子按照|p_k|排序，选出最大的一个或多个相关因子，将其编号传递给需要的分区。

Images