CN101877482A - 基于慢模式特征值对于线路参数灵敏度的弱连接识别方法 - Google Patents

Abstract

本发明公开了属于电力系统分析和控制技术领域的一种基于慢模式相关特征值对于线路参数灵敏度的弱连接识别方法，其特征在于：首先计算系统的慢模式相关特征值组对各条线路参数变化灵敏度，再选取灵敏度数值比较大的线路作为系统的弱连接。主要按照以下几个步骤进行：计算机初始化，输入包括系统网络参数、节点电压和功率信息的初始化数据；运用牛顿拉夫逊法计算系统潮流；形成系统线性化的动态方程；运用QR方法求取系统特征值，并根据大小进行排序，选取系统慢模式相关特征值组；计算慢模式相关特征值对于系统各条线路参数变化的灵敏度，并按大小进行排序；选取系统弱连接。本发明的有益效果是相比于以前首先处理节点的方法，能够直接处理线路。

Description

基于慢模式特征值对于线路参数灵敏度的弱连接识别方法

技术领域

本发明属于电力系统分析和控制技术领域，特别涉及一种基于慢模式特征值对于线路参数灵敏度的弱连接识别方法。

背景技术

解列是电力系统受到大干扰，在系统完整性无法保持的情况下将系统分成几个孤岛的控制。以在线计算、全局控制、协调动作为特点的主动解列控制是一个重要的研究方向。基于电力系统慢同调的主动解列思路认为，如果解列操作发生在同调机群之间的弱连接上对于系统的冲击最小，解列后孤岛更容易存活。基于这种思路，弱连接识别在基于电力系统慢同调的主动解列控制中占据着举足轻重的地位。对于弱连接识别，以前的做法是首先对节点进行同调分群，各个同调群之间的连接就是弱连接，但是这种做法有以下缺陷：

1)负荷节点处理会遇到困难。尽管引入了广义特征值概念，能够将负荷节点加以考虑，但是仍然面临着发电机节点和负荷节点对系统慢模式特征值的灵敏度在数值上差别很大以及阀值难以选取的困难；

2)这种基于节点灵敏度的弱连接识别方法直接对节点进行操作，并没有直接处理线路，“弱连接”的概念并不能直接地显现出来。但是，解列最后操作的对象是线路，这种通过首先处理节点再确定线路的做法并不能直观反映出所挑选线路的特性。

本发明的理论基础就是电力系统慢同调理论，电力系统中同调的定义有多种，本发明采用的是一种模式不可观性定义。

同调定义：对于系统线性

x＝[x₁ x₂ ... x_n]^T为其的状态变量组，

为系统第i个状态变量对时间的一阶导数，即

σ_r＝{λ₁ λ₂ ... λ_r}为系统的一组特征值模式，如果x中的两个状态变量x_i(t)，x_j(t)满足y_ij(t)＝x_i(t)-x_j(t)不含σ_r中的任何一个模式分量，那么就称x_i(t)，x_j(t)关于模式组σ_r是同调的。如果系统的状态变量x＝[x₁ x₂ ... x_n]^T可以分成r组，每一组内的任何两个变量之间都是关于模式组σ_r同调的，那么就称这个系统是关于模式组σ_r同调的。同时如果σ_r＝{λ₁ λ₂ ... λ_r}是系统的一组绝对值最小的特征值，那么就说系统是严格慢同调的。

这种同调定义直接应用到电力系统还有障碍，发电机取经典二阶模型同时忽略阻尼的情况下，电力系统的线性化动态方程的形式为：

Δδ＝[Δδ₁ Δδ₂ … Δδ_nG]^T为系统所有发电机功角在平衡点的偏差列向量，nG为系统发电机的台数；

为系统所有发电机角速度在平衡点的偏差列向量；

下文中，某个列向量的上标有一个黑点“·”都表示此列向量中每一个元素对时间一阶导数组成的列向量。

可以进一步忽略角速度项：

其中：A＝M^-1K，M＝diag{M₁，M₂，...，M_n}是由所有发电机惯性时间常数组成的对角矩阵，M_i＝2H_i(i＝1，2，..，n)，K是系统的网络矩阵；

下文中某个列向量上标有两个黑点“··”都表示此列向量中每一个元素对于时间的两阶导数组成的列向量。

这是一个二阶微分方程的形式，在一阶微分方程形式下的同调定义应用到二阶微分方程的形式要做相应的修改。同样经过推导可以得到如下结论：A的特征值组σ_r＝{λ₁，λ₂，..λ_r}并不是系统真正的机电振荡模式，而对应的真正的机电振荡模式是为了与模式特征值相区别，称σ_r＝{λ₁，λ₂，..λ_r}中的每个特征值为模式相关特征值。进一步，如果σ_r＝{λ₁，λ₂，..λ_r}是A的绝对值最小的一组特征值，那么称σ_r＝{λ₁，λ₂，..λ_r}为慢模式相关特征值组。如果说A关于慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}是慢同调的，意思就是系统关于模式特征值组

是慢同调的。

进一步得到在

形式下系统的同调定义：对于系统线性

X＝[x₁ x₂ ... x_n]^T为其的状态变量，

为系统第i个状态变量对时间的二阶导数，即

σ_r＝{λ₁ λ₂ ... λ_r}为系统的一组特征值，如果x中的两个状态变量x_i(t)，x_j(t)满足y_ij(t)＝x_i(t)-x_j(t)不含

中的任何一个模式分量，那么就称x_i(t)，x_j(t)关于特征值组σ_r是同调的。如果系统的状态变量x＝[x₁ x₂ ... x_n]^T可以分成r组，每一组内的任何两个变量之间都是关于模式组σ_r同调的，那么就称这个系统是关于模式组σ_r同调的。同时如果σ_r＝{λ₁ λ₂ ... λ_r}是系统的一组绝对值最小的特征值，那么就说系统是严格慢同调的。

对于A的一组慢特征值σ_r＝{λ₁，λ₂，..λ_r}，如果系统能够关于σ_r是严格慢同调的，那么A的对应于σ_r＝{λ₁，λ₂，..λ_r}的特征相量矩阵V中同调机组对应的行的各个元素都相同，这是一个很重要的结论。

由系统的严格慢同调可以推出一些有益的结论。在系统严格慢同调的条件下，假设系统的r个同调机群已知，每个同调机群选择一个参考状态变量组成X₁＝[x₁ x₂ ... x_r]^T，称X₁为参考变量相量，那么剩下的n-r个状态变量组成X₂＝[x_r+1 x_r+2 ... x_n]^T，X₂为剩余变量相量。那么系统的状态变量相量就按照如下的顺序排列： $x={\left[\begin{array}{cc}{X}_1^T& X_2^T\end{array}\right]}^T.$

下面定义一个分群矩阵L_g，L_g具有以下特征：

1)L_g是一个(n-r)×r维的矩阵，每一行代表每一个待分群的状态(X₂＝[x_r+1 x_r+2 ... x_n]^T中的某个状态)，每一列代表着以某个状态(X₁＝[x₁ x₂ ... x_r]^T)为参考状态的机群。

2)L_g的第i(i＜＝n-r)元素按照如下的方式定义：如果状态x_r+i属于以x_j(j＜＝r)为参考变量的机群，那么L_g(i，j)＝1，否则L_g(i，j)＝0。

对于进行下面两次线性变换：

$\left[\begin{array}{c}{x}_1\\ y_d\end{array}\right]=\left[\begin{array}{cc}{I}_r& 0\\ -L_g& I_{n-r}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$ 和 $\left[\begin{array}{c}{y}_a\\ y_d\end{array}\right]=\left[\begin{array}{cc}{I}_r& H\\ 0& I_{n-r}\end{array}\right]\left[\begin{array}{c}{x}_1\\ y_d\end{array}\right]$

其中： $H=M_a^{-1}{L}_g^T{M}_2,$ 而 $M_a=\left[\begin{array}{cc}{I}_r& L_g^T\end{array}\right)M\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right],$ I_r为r×r维的单位矩阵，M₂＝diag{M_r+1，M_r+2，...，M_n}。

根据系统严格慢同调条件，经过严格的推导可以得到如下结论：

1)能够实现系统的区域分解和聚合。

其中： $A_a=\left[\begin{array}{cc}{I}_r-H{L}_g& H\end{array}\right]A\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right],$ $A_d=\left[\begin{array}{cc}-L_g& I_{n-r}\end{array}\right]A\left[\begin{array}{c}-H\\ (I_{n-r}-L_{g}H)\end{array}\right].$

即实现了系统区域分解和聚合，聚合模型为

y_a是各个区域的惯性中心变量组成的相量，而区域分解模型为y_d是各个区域的内部动态。

特别地，经过进一步的推导，在区域聚合模型

中， $A_a=M_a^{-1}{K}_a,$ 其中 $K_a=\left[\begin{array}{cc}{I}_r& L_g\end{array}\right]K\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right],$ $M_a=\left[\begin{array}{cc}{I}_r& L_g^{\′}\end{array}\right]M\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right].$

2)σ_a＝λ(A_a)＝σ_r＝{λ₁，λ₂，..λ_r}。这里的σ′定义为

就是所谓的快模式特征值，有结论：σ_d＝λ(A_d)＝σ′。即系统得慢模式相关特征值σ_r与系统的区域聚合(简化)的特征值是相同的。

继续上面的推导，K作为网络矩阵具有下列的性质：

$\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{k}_{\mathrm{ij}}=0,$ k_ij＝k_ji＞0，i≠j

也就是说K的非正元素只有对角元素k_ii，并且它与本行的非对角元素之和的相反数相等。可以知道K的特征值中有0特征值，且对应的特征向量为：u＝[1 1 ，1]^T。在已知系统分群的情况下对于系统的状态变量x重新排序，以使得在排序后的状态向量Px中，同一区域的状态量在一起。假设系统分成r个区域，重新排序后的系统矩阵具有如下的形式： $K_P={\mathrm{PKP}}^{\′}=\left[\begin{array}{cccc}{K}_{11}& K_{12}& ...& K_{1r}\\ ...& ...& ...& ...\\ K_{r1}& ...& ...& K_{\mathrm{rr}}\end{array}\right].$

K_αα表示区域内部的联结，而K_αβ(β≠α)表示区域α，β之间的连接。重写K_P以使内部和外部连接区分开来：K_P＝K^I+K^E。注意：这里面的K^I，K^E不是简单的取对角块和非对角块的结果。而是K^I中减去了对于与区域之间连接的元素，这部分元素加到了K^E上。也就是说K^E在原来K^I对应的对角块位置上的元素并不全是0。

这样得到的K^I有如下特征：

$K^{1}U=K^1\left[\begin{array}{ccccc}{u}_1& 0& ...& 0& 0\\ 0& u_2& ...& ...& ...\\ ...& ...& ...& ...& ...\\ ...& ...& ...& ...& 0\\ 0& 0& ...& 0& u_r\end{array}\right]=0,$ 其中u_α＝[1 1 ... 1]^T

已知： $A_a=M_a^{-1}{K}_a,$ $K_a=\left[\begin{array}{cc}{I}_r& L_g^{\′}\end{array}\right]K\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right].$ 这里的K是没有经过重新排序的系统矩阵。而K_P＝K^I+K^E是重新排序后的系统矩阵，且排序矩阵为P，即K＝P′K_PP。

这样： $K_a=\left[\begin{array}{cc}{I}_r& L_g^{\′}\end{array}\right]K\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right]=\left[\begin{array}{cc}{I}_r& L_g^{\′}\end{array}\right]P^{\′}{K}_{P}P\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right]=\left[\begin{array}{cc}{I}_r& L_g^{\′}\end{array}\right]P^{\′}(K^I+K^E)P\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right].$

而能够推得： $P\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right]=U=\left[\begin{array}{ccccc}{u}_1& 0& ...& 0& 0\\ 0& u_2& ...& ...& ...\\ ...& ...& ...& ...& ...\\ ...& ...& ...& ...& 0\\ 0& 0& ...& 0& u_r\end{array}\right]$

又由于： $K^{1}U=K^1\left[\begin{array}{ccccc}{u}_1& 0& ...& 0& 0\\ 0& u_2& ...& ...& ...\\ ...& ...& ...& ...& ...\\ ...& ...& ...& ...& 0\\ 0& 0& ...& 0& u_r\end{array}\right]=0$

所以，就有： $K_a=\left[\begin{array}{cc}{I}_r& L_g^{\′}\end{array}\right]P^{\′}{K}^{E}P\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right]$

进而：

$A_a=M_a^{-1}{K}_a=M_a^{-1}\left[\begin{array}{cc}{I}_r& L_g^{\′}\end{array}\right]P^{\′}{K}^{E}P\left[\begin{array}{c}{I}_r\\ L_g\end{array}\right]---\left(4\right)$

又知λ(A_a)＝σ_r＝{λ₁，λ₂，..λ_r}，所以可以得到了一个很重要的结论：系统的模式特征值只与系统的同调区域之间的弱连接K^E有关，而内连接K^I对它没有任何影响。

可以把系统区域之间的弱连接定义成K^E所对应的线路，但是我们知道以上的系统模型并不是结构保留模型，而是收缩到发电机节点的网络，所以K^E并不对应着实际的物理线路，这样就无法直接通过K^E来识别弱连接了。但是，上面的结论还是能够提供很有益的启示：区域之间的弱连接线路对于系统的慢模式特征值影响比较大，而区域内部的线路对于系统的慢模式影响比较小。

发明内容

本发明提出了一种基于慢模式特征值对于线路参数灵敏度的弱连接识别方法，其特征在于，该弱连接识别方法能够直接处理线路，通过计算系统的一组慢模式对应的特征值对于各条线路参数变化灵敏度，选取灵敏度数值比较大的线路作为系统的弱连接；在计算机上按如下步骤进行：

步骤(1)计算机初始化；

1.1)输入以下参数：

电力系统的网络参数，包括：输电线路的串联电阻、串联电抗、并联电导和并联电纳；变压器的变比和阻抗；并联在输电线路上的电容器和电抗器的阻抗；发电机的动态参数包括发电机的惯性时间常数；发电机的暂态电抗；

1.2)设定以下参数：各母线节点的电压和有功功率和无功功率的实时量测值；

步骤(2)根据步骤(1)的参数运用牛顿拉夫逊方法计算系统潮流，并记录至少包括系统各个节点电压的幅值和相位当前潮流断面信息；

步骤(3)根据步骤(1)、步骤(2)的结果形成多机系统线性化的动态模型方程

其中关键是计算矩阵A，其中

表示发电机功角偏差向量的二阶导数；

形成

按以下思路进行：

为了计算和推导的方便，本发明采用了简化的电力系统模型，具体模型如下：研究电力系统同调时有一个基本的假设：发电机模型的复杂程度不改变同调机群之间基本的振荡模态；所以本发明中选取很简单的电力系统模型来设计算法，其中发电机取经典二阶模型表同时忽略阻尼；而负荷模型则采用恒阻抗负荷模型；

发电机取经典二阶模型并忽略阻尼的情况下，电力系统的机电振荡模型如下：

$\left\{\begin{array}{c}2H\frac{d^2\mathrm{\δ}}{{\mathrm{dt}}^2}=P_m-P_e\\ P_e=f(E_d,U_G,x_d^{\′},\mathrm{\δ},{\mathrm{\θ}}_{\mathrm{UG}})\end{array}\right.---\left(1\right)$

第一个等式中：δ为所有发电机功角组成的列相量， $\mathrm{\δ}={\left[\begin{array}{cccc}{\mathrm{\δ}}_1& {\mathrm{\δ}}_2& ...& {\mathrm{\δ}}_{n_G}\end{array}\right]}^T,$ n_G为系统中发电机的个数；H为所有发电机惯性时间常数组成的对角矩阵， $H=\mathrm{diag}\left\{\begin{array}{cccc}{H}_1& H_2& ...& H_{n_G}\end{array}\right\};$ P_m为所有发电机的机械输入功率组成的列向量， $P_m=\left[\begin{array}{cccc}{P}_{m1}& P_{m2}& ...& P_{{\mathrm{mn}}_G}\end{array}\right],$ 为常数相量；P_e为所有发电机的电磁输出功率组成的列向量， $P_e=\left[\begin{array}{cccc}{P}_{e1}& P_{e2}& ...& P_{{\mathrm{en}}_G}\end{array}\right];$

第二个等式中： $P_{\mathrm{ei}}=\frac{E_{\mathrm{di}}{U}_{\mathrm{Gi}}}{x_{\mathrm{di}}^{\′}}\sin ({\mathrm{\δ}}_i-{\mathrm{\θ}}_{\mathrm{UGi}}),$ E_di为发电机i的暂态电抗后电势，x′_di为发电机i的暂态电抗，U_gi为发电机i的机端电压；U_gi为发电机i的机端电压，θ_Ugi为发电机i的机端电压的相角；所以f(E_d，U_G，x′_d，δ，θ_UG)是一个矩阵表达式，其形式如下：

$f(E_d,U_G,x_d^{\′},\mathrm{\δ},{\mathrm{\δ}}_{\mathrm{UG}})={\left[\begin{array}{cccc}\frac{E_{d1}{U}_{g1}}{x_{d1}^{\′}}\sin ({\mathrm{\δ}}_1-{\mathrm{\δ}}_{\mathrm{Ug}1})& \frac{E_{d2}{U}_{g2}}{x_{d2}^{\′}}\sin ({\mathrm{\δ}}_2-{\mathrm{\δ}}_{\mathrm{Ug}2})& ...& \frac{E_{\mathrm{dnG}}{U}_{\mathrm{gnG}}}{x_{\mathrm{dnG}}^{\′}}\sin ({\mathrm{\δ}}_{\mathrm{nG}}-{\mathrm{\δ}}_{\mathrm{Ug}{n}_G})\end{array}\right]}^T$ 其中 $E_d={\left[\begin{array}{cccc}{E}_{d1}& E_{d2}& ...& E_{{\mathrm{dn}}_G}\end{array}\right]}^T,$ $x_d^{\′}={\left[\begin{array}{cccc}{x}_{d1}^{\′}& x_{d2}^{\′}& ...& x_{{\mathrm{dn}}_G}^{\′}\end{array}\right]}^T,$ $U_G={\left[\begin{array}{cccc}{U}_{g1}& U_{g2}& ...& U_{{\mathrm{gn}}_G}\end{array}\right]}^T,$ ${\mathrm{\θ}}_{\mathrm{UG}}={\left[\begin{array}{cccc}{\mathrm{\θ}}_{\mathrm{Ug}1}& {\mathrm{\θ}}_{\mathrm{Ug}2}& ...& {\mathrm{\θ}}_{U{\mathrm{gn}}_G}\end{array}\right]}^T;$

每台发电机的暂态电抗后电势及当前功角按照如下思路确定：

假设发电机节点i的功率注入为 ${\stackrel{*}{S}}_{\mathrm{gi}}=P_{\mathrm{gi}}+j{Q}_{\mathrm{gi}},$ 同时机端电压模制和相角为U_gi，δ_Ugi，那么首先求出此节点的等效注入电流 ${\stackrel{*}{I}}_{\mathrm{gi}}=\left(\frac{S_{\mathrm{gi}}}{U_{\mathrm{gi}}\∠{\mathrm{\θ}}_{\mathrm{gi}}}\right),$ 那么： $E_{\mathrm{di}}\∠{\mathrm{\δ}}_i=U_{\mathrm{gi}}\∠{\mathrm{\θ}}_{\mathrm{Ugi}}+{\stackrel{*}{I}}_{\mathrm{gi}}*\left({\mathrm{jx}}_{\mathrm{di}}\right)$ 下文中某一个变量有上标星号(*)表示此变量是一个复数矢量；

对于系统网络的处理，对于发电机和负荷分别进行简化处理；对于发电机，经典二阶模型下的发电机相当于有内阻的电压源(内阻是x′_di，而电压是E_di)，那么经过诺顿等值将这个电压源转换成一个电流源：阻抗仍为x′_di，

$I_{\mathrm{gi}}=\frac{E_{\mathrm{di}}}{x_{\mathrm{di}}^{\′}},$ 对于负荷节点k进行如下处理：假设当前潮流下，节点的电压模值为V_{Load_k}，而负荷有功和无功分别为P_k和Q_k，那么等效负荷阻抗为 $Z_{\mathrm{Load}\_k}={\left(\frac{V_{\mathrm{Load}\_k}^2}{P_k+j{Q}_k}\right)}^{*}=R_{\mathrm{Load}\_k}+j{X}_{\mathrm{Load}\_k},$ 相应得负荷等效导纳为 $Y_{\mathrm{Load}\_k}=\frac{1}{Z_{\mathrm{Load}\_k}}=G_{\mathrm{Load}\_k}+j{B}_{\mathrm{Load}\_k},$

把负荷节点的等效阻抗以及发电机的暂态电抗加以考虑之后，系统的节点导纳矩阵变为：Y＝G_Y+jB_Y；其中：Y、G_Y、B_Y都是n×n维的矩阵；令 $\stackrel{*}{U}=U_x+{\mathrm{jU}}_y$ 表示个节点电压列向量，每个元素为各个节点的电压矢量； $\stackrel{*}{I}=I_x+{\mathrm{jI}}_y$ 代表各个节点的注入电流向量，下文中，某个列向量上标是一个星号(*)都表示此列向量每个元素都是复数矢量；其中：U_x为各个节点电压实部组成的向量，U_y为各个节点电压虚部组成的向量，I_x为各个节点注入电流实部组成的向量，I_y为个节点注入电流虚部组成的向量；

经过上述简化之后，只有发电机节点具有注入电流，其他的节点没有电流注入；

令 $Y^{\′}=\left[\begin{array}{cc}{G}_Y& -B_Y\\ B_Y& G_Y\end{array}\right],$ $U^{\′}=\left[\begin{array}{c}{U}_x\\ U_y\end{array}\right],$ $I^{\′}=\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right],$ 那么有： $I^{\′}=\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right]=Y^{\′}U=\left[\begin{array}{cc}{G}_Y& -B_Y\\ B_Y& G_Y\end{array}\right]\left[\begin{array}{c}{U}_x\\ U_y\end{array}\right]$

令系统所有节点的电压模值向量为U＝[U₁ U₂ ... U_N]，相角向量为θ_U＝{θ_u1 θ_u2 ... θ_uN}^T，N为系统所有节点的个数，那么：

U_x＝[U₁cosθ_u1 U₂cosθ_u2 ... U_ncosθ_uN]^T，U_y＝[U₁sinθ_u1 U₂sinθ_u2 ... U_nsinθ_uN]^T

令： $Z^{\′}=Y^{\′-1}{\left[\begin{array}{cc}{G}_Y& -B_Y\\ B_Y& G_Y\end{array}\right]}^{-1}$ 那么有：

$\left[\begin{array}{c}{U}_x\\ U_y\end{array}\right]=Z\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right]---\left(2\right)$

I_x，I_y都只有发电机节点对应元素非零，其他元素均为零；令：

$I_{\mathrm{Gx}}={\left[\begin{array}{cccc}{I}_{G1}\sin {\mathrm{\δ}}_1& I_{G2}\sin {\mathrm{\δ}}_2& ...& I_{{\mathrm{Gn}}_G}\sin {\mathrm{\δ}}_{n_G}\end{array}\right]}^T$ (1，...，n_G表示发电机序号)

$I_{\mathrm{Gy}}={\left[\begin{array}{cccc}{I}_{G1}\cos {\mathrm{\δ}}_1& I_{G2}\cos {\mathrm{\δ}}_2& ...& I_{{\mathrm{Gn}}_G}\cos {\mathrm{\δ}}_{n_G}\end{array}\right]}^T$ (1，...，n_G表示发电机序号)

那么I_x，I_y可以用以下的形式表示：I_x＝T×I_Gx，I_y＝T×I_Gy

其中：T是一个N×n_G维的矩阵；T的每一行对应系统的一个节点(编号为k)，如果本节点没有发电机接入，即非发电机节点，那么本行元素全为0，即T(k，：)＝[0...0]；如果本节点有发电机接入，即为发电机节点，假设此发电机序号为i，那么T(k，i)＝1，本行其他元素均为0，即T(k，：)＝[0...1...0]；

对(2)式两边线性化得：

$\left[\begin{array}{cc}\cos {\mathrm{\θ}}_U& -U\sin {\mathrm{\θ}}_U\\ \sin {\mathrm{\θ}}_U& U\cos {\mathrm{\θ}}_U\end{array}\right]\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_U\end{array}\right]=Z\left[\begin{array}{c}{I}_{\mathrm{x\δ}}\\ I_{\mathrm{y\δ}}\end{array}\right]\mathrm{\Δ\δ}---\left(3\right)$

其中：

ΔU，Δθ_U，Δδ分别为U，θ_U，δ的偏列差量，N和n_G分别为系统所有节点的个数以及发电机节点的个数，表达式中的各个矩阵的形式如下：

cosθ_U＝diag{cosθ_U1 cosθ_U2 ... cosθ_UN}；sin_θU＝diag{sinθ_U1 sinθ_U2 ... sinθ_UN}；

Ucosθ_U＝diag{U₁cosθ_U1 U₂cosθ_U2 ... U_Ncosθ_UN}；

Usinθ_U＝diag{U₁sinθ_U1 U₂sinθ_U2 ... U_Nsinθ_UN}；

$I_{\mathrm{x\δ}}=\frac{\∂I_x}{\∂\mathrm{\δ}},$ $I_{\mathrm{y\δ}}=\frac{\∂I_y}{\∂\mathrm{\δ}}$ 均为N×n_G维的矩阵，可知 $I_{\mathrm{x\δ}}=\frac{\∂I_x}{\∂\mathrm{\δ}}=T\×\frac{\∂I_{\mathrm{Gx}}}{\∂\mathrm{\δ}},$ $I_{\mathrm{y\δ}}=\frac{\∂I_y}{\∂\mathrm{\δ}}=T\×\frac{\∂I_{\mathrm{Gy}}}{\∂\mathrm{\δ}},$ 而：

$\frac{\∂I_{\mathrm{Gx}}}{\∂\mathrm{\δ}}=\mathrm{diag}\left\{\begin{array}{cccc}{I}_{G1}\cos {\mathrm{\δ}}_1& I_{G2}\cos {\mathrm{\δ}}_2& ...& I_{\mathrm{GnG}}\cos {\mathrm{\δ}}_{n_G}\end{array}\right\}$

$\frac{\∂I_{\mathrm{Gy}}}{\∂\mathrm{\δ}}=\mathrm{diag}\left\{\begin{array}{cccc}{I}_{G1}\sin {\mathrm{\δ}}_1& I_{G2}\sin {\mathrm{\δ}}_2& ...& I_{\mathrm{GnG}}\sin {\mathrm{\δ}}_{n_G}\end{array}\right\}$

令： $D=\left[\begin{array}{cc}\cos {\mathrm{\θ}}_U& -U\sin {\mathrm{\θ}}_U\\ \sin {\mathrm{\θ}}_U& U\cos {\mathrm{\θ}}_U\end{array}\right],$ $E=\left[\begin{array}{c}{I}_{\mathrm{x\δ}}\\ I_{\mathrm{y\δ}}\end{array}\right],$ E是2N×n_G维的，那么：

$\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_U\end{array}\right]=D^{-1}\mathrm{ZE\Δ\δ}---\left(4\right)$

对于式(1)所示的电力系统机电振荡模型在平衡点处进行线性化，可得：

$\frac{d^2\mathrm{\Δ\δ}}{{\mathrm{dt}}^2}=B\left[\begin{array}{c}\mathrm{\Δ}{U}_G\\ \mathrm{\Δ}{\mathrm{\θ}}_{\mathrm{UG}}\end{array}\right]+\mathrm{C\Δ\δ}---\left(5\right)$

其中：ΔU_G，Δθ_UG，Δδ为U_G，θ_UG，δ的偏差量， $B={\left(2H\right)}^{-1}\left[\begin{array}{cc}-\frac{\∂g}{\∂U_G}& 0\\ 0& -\frac{\∂f}{\∂{\mathrm{\θ}}_{\mathrm{UG}}}\end{array}\right],$ $G=-{\left(2G\right)}^{-1}\frac{\∂f}{\∂\mathrm{\δ}}$ 有：ΔU_G＝T₂×ΔU，Δθ_UG＝T₂×Δθ_U，代入(5)式有：

$\frac{d^2\mathrm{\Δ\δ}}{{\mathrm{dt}}^2}=B\left[\begin{array}{c}\mathrm{\Δ}{U}_G\\ \mathrm{\Δ}{\mathrm{\θ}}_{\mathrm{uG}}\end{array}\right]+\mathrm{C\Δ\δ}=B\left[\begin{array}{cc}{T}_2& 0\\ 0& T_2\end{array}\right]\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_u\end{array}\right]+\mathrm{C\Δ\δ}={\mathrm{BT}}_2^{\′}\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_U\end{array}\right]+\mathrm{C\Δ\δ}---\left(6\right)$

其中： $T_2^{\′}=\left[\begin{array}{cc}{T}_2& 0\\ 0& T_2\end{array}\right],$ T₂＝T^T

带入(4)式就得到了多机系统线性化模型为：

$\frac{d^2\mathrm{\Δ\δ}}{{\mathrm{dt}}^2}=({\mathrm{BT}}_2^{\′}{D}^{-1}\mathrm{ZE}+C)\mathrm{\Δ\δ}=\mathrm{A\Δ\δ}---\left(7\right)$

步骤(4)求取矩阵A的特征值并选择慢模式相关特征值组；

4.1)采用QR算法求取A的全部特征值；

4.2)把A的特征值按照绝对值由小到大的顺序进行排序；

假设排序后的特征值顺序如下： $\left\{\begin{array}{cccc}{\mathrm{\λ}}_1& {\mathrm{\λ}}_2& ...& {\mathrm{\λ}}_{n_G}\end{array}\right\},$ 其中n_G为发电机的台数，也正是A的维数；值得注意的是λ₁＝0；

4.3)选取慢模式相关特征值组，关键是确定选取慢模式相关特征值的个数r；

按照如下的方法选择r：比较所有的λ_i/λ_i+1，i＝2，3，...，n_G-1，令： ${\mathrm{\λ}}_r/{\mathrm{\λ}}_{r+1}=\underset{i=\mathrm{2,3},..n_G-1}{\min }{\mathrm{\λ}}_i/{\mathrm{\λ}}_{i+1},$ 这样就确定了慢模式相关特征值的个数r；这种选择的依据是电力系统多时间尺度特性原理，λ_r/λ_r+1的值越小表示系统的多时间尺度特性越明显，系统的区域性也越明显，系统也越适合解列；经过这一步骤，就确定了慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}；

步骤(5)求取慢模式相关特征值组对于每条线路参数变化灵敏度；

线路就是指电力系统中实际的输电线路。首先对于所有节点进行编号：1，2，...，N，再对线路进行编号：1，2，...，Nl。一条线路主要有如下元素特征：起始节点编号、终了节点编号、线路电阻r、线路电抗x，线路电导g和线路电纳b，值得注意的是在选取线路起止节点时，哪一端点作为起始节点时无所谓的，所以在确定起始节点时随便指定其中一个端点作为起始节点即可。

某条线路的编号是L，其起始节点编号为I，终了节点编号为J，线路的电阻和电抗分别为r和x，电导和电纳分别为g和b，求取慢模式相关特征值组对于线路L参数变化的灵敏度，可按以下两步进行：

5.1)求取慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}中每个特征值对线路L参数变化的灵敏度，表示为{S_λ1 ^L S_λ2 ^L ... S_λr ^L}，S_λi ^L表示特征值λ_i对于线路L参数变化的灵敏度；

特征值λ_i对于线路L参数变化的灵敏度可按照以下思路进行：

λ_i为A的一个特征值，λ_i对应的左右特征向量分别为φ_i(n_G×1)和ψ_i(1×n_G)；

那么λ_i对于线路L的导纳b变化的灵敏度表达式如下：

其中： $N^{\′}=\left[\begin{array}{cc}0& M^T\\ M^T& 0\end{array}\right],$ $M^{\′}=\left[\begin{array}{cc}M& 0\\ 0& -M\end{array}\right],$ M＝[... 1 ... -1 ...]^T，其中1的位置是I，而-1的位置是J；

λ_i对于线路L的导纳g变化的灵敏度表达式如下：

其中： $N^{\′\′}=\left[\begin{array}{cc}{M}^T& 0\\ 0& M^T\end{array}\right],$ $M^{\′\′}=\left[\begin{array}{cc}M& 0\\ 0& M\end{array}\right];$

假设线路参数变化时，r和x同比例变化，那么由于 $g+\mathrm{jb}=\frac{1}{r+\mathrm{jx}}=\frac{r-\mathrm{jx}}{r^2+x^2},$ 线路的g和b是按照-r/x的比例变化的，以b的变化(高压线路b远大于g)为基准，总结出λ_i对于线路L参数变化的灵敏度变化表达式如下：

$S_{\mathrm{\λi}}^L=-S_{\mathrm{\λi}}^{\mathrm{Lb}}+\frac{r}{x}{S}_{\mathrm{\λi}}^{\mathrm{Lg}}---\left(10\right)$

5.2)确定慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}对于线路L参数变化的灵敏度，按下式确定：

$S_{\mathrm{\σr}}^L=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{S}_{\mathrm{\λi}}^L---\left(11\right)$

S_σr ^L即为慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}对于线路L参数变化的灵敏度，用这个量来衡量线路参数变化对于慢模式组的影响；

经过这一步骤就能够得到所有线路的慢模式相关特征值组对于系统所有线路的参数的灵敏度信息；

步骤(6)：根据步骤5得到的慢模式相关特征值组对于线路参数变化的灵敏度信息来确定弱连接；

6.1)把慢模式相关特征值组对于系统所有的Nl.条线路的灵敏度按照绝对值从大到小的顺序进行排序，排序后的顺序为{S_σr ¹ S_σr ² ... S_σr ^Nl}；

6.2)选取弱连接线路，基本思路就是选取灵敏度比较大的一组线路作为系统弱连接线路；至于弱连接线路的条数wl可以由用户自己按照需要设定，也可以自动生成，本发明推荐一种自动生成wl的方法；对于一个多时间尺度特性比较明显的系统，弱连接线路与非弱连接线路的灵敏度数值差异会比较大；这种方法首先对于系统所有线路的灵敏度进行从大到小的排序，然后按照上述原理来选择弱连接线路的条数wl，选取wl的公式表达形式如下：

$\mathrm{abs}(S_{\mathrm{\σr}}^{\mathrm{wl}}/S_{\mathrm{\σr}}^{\mathrm{wl}+1})=\underset{i=\mathrm{1,2},...\mathrm{Nl}}{\max }\mathrm{abs}(S_{\mathrm{\σr}}^i/S_{\mathrm{\σr}}^{i+1})---\left(12\right)$

其中，符号abs表示取一个数的绝对值。

本发明的有益效果是该弱连接识别方法能够直接处理线路，通过计算系统的一组慢模式对应的特征值对各条线路参数变化灵敏度，选取灵敏度数值比较大的线路作为系统的弱连接。

附图说明

图1算法流程图。

图2发电机诺顿等值示意图

图3网络处理效果图

图4Kundur11节点测试系统拓扑图及相应的弱连接示意图

具体实施方式

本发明提出了一种基于慢模式特征值对线路参数灵敏度的弱连接识别方法。本节利用Kundur11节点系统来测试本发明的选取系统弱连接的方法。Kundur11节点系统包括11条母线，4台发电机，12条线路(包含变压器支路)，系统拓扑图见图4。

参照附图1所示的流程步骤识别此系统的弱连接：

步骤(1)开始，计算机初始化；

1.1)输入以下参数

电力系统的网络参数，包括：输电线路的串联电阻、串联电抗、并联电导和并联电纳；变压器的变比和阻抗；并联在输电线路上的电容器和电抗器的阻抗；发电机的动态参数包括发电机的惯性时间常数；发电机的暂态电抗；

1.2)设定以下参数：各母线节点的电压和有功功率和无功功率的实时量测值；

步骤(2)根据步骤(1)的参数运用牛顿拉夫逊方法计算系统潮流，并记录至少包括系统各个节点电压的幅值和相位当前潮流断面信息；

步骤(3)本发明中采用系统简化的模型来进行计算，首先需要对网络进行处理，对于发电机和负荷分别进行简化处理；对于发电机，经典二阶模型下的发电机相当于有内阻的电压源(内阻是x′_di，而电压是E_di)，那么经过诺顿等值将这个电压源转换成一个电流源：阻抗仍为x′_di，

$I_{\mathrm{gi}}=\frac{E_{\mathrm{di}}}{x_{\mathrm{di}}^{\′}},$ 诺顿等值示意图参见图2；对于负荷节点k进行如下处理：假设当前潮流下，节点的电压模值为V_{Load_k}，而负荷有功和无功分别为P_k和Q_k，那么等效负荷阻抗为 $Z_{\mathrm{Load}\_k}={\left(\frac{V_{\mathrm{Load}\_k}^2}{P_k+j{Q}_k}\right)}^{*}=R_{\mathrm{Load}\_k}+j{X}_{\mathrm{Load}\_k},$ 相应得负荷等效导纳为 $Y_{\mathrm{Load}\_k}=\frac{1}{Z_{\mathrm{Load}\_k}}=G_{\mathrm{Load}\_k}+j{B}_{\mathrm{Load}\_k};$ 此网络处理效果的示意图参见图3；

经过以上的网络处理后，根据式(7) $\frac{d^2\mathrm{\Δ\δ}}{{\mathrm{dt}}^2}=({\mathrm{BT}}_2^{\′}{D}^{-1}\mathrm{ZE}+C)\mathrm{\Δ\δ}=\mathrm{A\Δ\δ}$ 计算矩阵A＝BT′₂D^-1ZE+C其中各个矩阵的含义参见发明内容中的详细推导过程。由此得到：

$A=\left[\begin{array}{cccc}-0.0683& 0.0588& 0.0043& 0.0051\\ 0.0660& -0.0836& 0.0077& 0.0099\\ 0.0044& 0.0052& -0.0714& 0.0618\\ 0.0080& 0.0102& 0.0698& -0.0879\end{array}\right]$

步骤(4)求取矩阵A的特征值并选择慢模式相关特征值组；

4.1)采用QR算法求取A的全部特征值及对应的左右特征向量并存储；

4.2)把A的特征值按照绝对值由小到大的顺序进行排序，所得结果如表1；

表1系统特征值

4.3)选取慢模式相关特征值组，关键是确定选取慢模式相关特征值的个数r；

按照如下的方法选择r：比较所有的λ_i/λ_i+1，i＝2，3，...，n_G-1，令： ${\mathrm{\λ}}_r/{\mathrm{\λ}}_{r+1}=\underset{i=\mathrm{2,3},..n_G-1}{\min }{\mathrm{\λ}}_i/{\mathrm{\λ}}_{i+1},$ 这样就确定了慢模式相关特征值的个数r；这种选择的依据是电力系统多时间尺度特性原理，λ_r/λ_r+1的值越小表示系统的多时间尺度特性越明显，系统的区域性也越明显，系统也越适合解列；经过这一步骤，就确定了慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}；

对于本实施例：

{λ₂/λ₃，λ₃/λ₄}＝{-0.0265/-0.1386-0.1386/-0.1461}＝{0.19120.9487}

$\underset{i=\mathrm{2,3}}{\min }\{{\mathrm{\λ}}_i/{\mathrm{\λ}}_{i+1}\}={\mathrm{\λ}}_2/{\mathrm{\λ}}_3=0.1912$

即：r＝2，所以σ_r＝{λ₁，λ₂}＝{0-0.0265}

步骤(5)求取慢模式相关特征值组对于每条线路参数变化灵敏度；

5.1)根据式(8)(9)(10)求取慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}对于各条线路参数变化的灵敏度；

5.2)根据式(11)确定慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}对于各条线路线路参数变化的灵敏度；

本实施例中，模式特征值组中各个特征值以及此特征值组整体对于各条线路参数变化灵敏度参见表2：

表2慢模式相关特征值对于线路参数变化的灵敏度

步骤(6)：根据步骤(5)得到的慢模式相关特征值组对于线路参数变化的灵敏度信息来确定弱连接；

6.1)把慢模式相关特征值组对于系统所有的Nl(本实施例中Nl＝12)条线路的灵敏度按照绝对值从大到小的顺序进行排序；排序后的结果参见表3：

表3慢模式相关特征值组对线路参数灵敏度排序后的结果

6.2)按照式(12)选取弱连接线路的条数wl：

$\mathrm{abs}(S_{\mathrm{\σr}}^{\mathrm{wl}}/S_{\mathrm{\σr}}^{\mathrm{wl}+1})=\underset{i=\mathrm{1,2},...12}{\max }\mathrm{abs}(S_{\mathrm{\σr}}^i/S_{\mathrm{\σr}}^{i+1})=S_{\mathrm{\σr}}^4/S_{\mathrm{\σr}}^5=0.5625/0.0825=6.812$

即：wl＝4，那么就得到此实施例中的弱连接线路如表4：

表4系统的弱连接线路

此实施例中系统的弱连接线路在图4中用虚线标识，参见图4。

Claims (1)

1.一种慢模式相关特征值对于线路参数灵敏度的弱连接识别方法，其特征在于，具体算法步骤如下：

步骤(1)计算机初始化；

1.1)输入以下参数：

电力系统的网络参数，包括：输电线路的串联电阻、串联电抗、并联电导和并联电纳；变压器的变比和阻抗；并联在输电线路上的电容器和电抗器的阻抗；发电机的动态参数包括发电机的惯性时间常数；发电机的暂态电抗；

1.2)设定以下参数：各母线节点的电压和有功功率和无功功率的实时量测值；

步骤(2)根据步骤(1)的参数运用牛顿拉夫逊方法计算系统潮流，并记录至少包括系统各个节点电压的幅值和相位当前潮流断面信息；

步骤(3)根据步骤(1)、步骤(2)的结果形成多机系统线性化的动态模型方程

其中关键是计算矩阵A，其中

表示发电机功角偏差向量的二阶导数；

形成

按以下思路进行：

为了计算和推导的方便，本发明采用了简化的电力系统模型，具体模型如下：研究电力系统同调时有一个基本的假设：发电机模型的复杂程度不改变同调机群之间基本的振荡模态；所以本发明中选取很简单的电力系统模型来设计算法，其中发电机取经典二阶模型表同时忽略阻尼；而负荷模型则采用恒阻抗负荷模型；

发电机取经典二阶模型并忽略阻尼的情况下，电力系统的机电振荡模型如下：

$\left\{\begin{array}{c}2H\frac{d^2\mathrm{\δ}}{{\mathrm{dt}}^2}=P_m-P_e\\ P_e=f(E_d,U_G,x_d^{\′},\mathrm{\δ},{\mathrm{\θ}}_{\mathrm{UG}})\end{array}\right.$

第一个等式中：δ为所有发电机功角组成的列相量， $\mathrm{\δ}={\left[\begin{array}{cccc}{\mathrm{\δ}}_1& {\mathrm{\δ}}_2& \·\·\·& {\mathrm{\δ}}_{n_G}\end{array}\right]}^T,$ n_G为系统中发电机的个数；H为所有发电机惯性时间常数组成的对角矩阵， $H=\mathrm{diag}\left\{\begin{array}{cccc}{H}_1& H_2& \·\·\·& H_{n_G}\end{array}\right\};$ P_m为所有发电机的机械输入功率组成的列向量， $P_m=\left[\begin{array}{cccc}{P}_{m1}& P_{m2}& \·\·\·& P_{{\mathrm{mn}}_G}\end{array}\right],$ 为常数相量；P_e为所有发电机的电磁输出功率组成的列向量， $P_e=\left[\begin{array}{cccc}{P}_{e1}& P_{e2}& \·\·\·& P_{e{n}_G}\end{array}\right];$

第二个等式中： $P_{\mathrm{ei}}=\frac{E_{\mathrm{di}}{U}_{\mathrm{Gi}}}{x_{\mathrm{di}}^{\′}}\sin ({\mathrm{\δ}}_i-{\mathrm{\θ}}_{\mathrm{UGi}}),$ E_di为发电机i的暂态电抗后电势，x′_di为发电机i的暂态电抗，U_gi为发电机i的机端电压；U_gi为发电机i的机端电压，θ_Ugi为发电机i的机端电压的相角；所以f(E_d，U_G，x′_d，δ，θ_UG)是一个矩阵表达式，其形式如下：

$f(E_d,U_G,x_d^{\′},\mathrm{\δ},{\mathrm{\δ}}_{\mathrm{UG}})={\left[\begin{array}{cccc}\frac{E_{d1}{U}_{g1}}{x_{d1}^{\′}}\sin ({\mathrm{\δ}}_1-{\mathrm{\δ}}_{\mathrm{Ug}1})& \frac{E_{d2}{U}_{g2}}{x_{d2}^{\′}}\sin ({\mathrm{\δ}}_2-{\mathrm{\δ}}_{\mathrm{Ug}2})& \·\·\·& \frac{E_{\mathrm{dnG}}{U}_{\mathrm{gnG}}}{x_{\mathrm{dnG}}^{\′}}\sin ({\mathrm{\δ}}_{\mathrm{nG}}-{\mathrm{\δ}}_{\mathrm{Ug}{n}_G})\end{array}\right]}^T$ 其由 $E_d={\left[\begin{array}{cccc}{E}_{d1}& E_{d2}& \·\·\·& E_{{\mathrm{dn}}_G}\end{array}\right]}^T,$ $x_d^{\′}={\left[\begin{array}{cccc}{x}_{d1}^{\′}& x_{d2}^{\′}& \·\·\·& x_{d{n}_G}^{\′}\end{array}\right]}^T,$ $U_G={\left[\begin{array}{cccc}{U}_{g1}& U_{g2}& \·\·\·& U_{{\mathrm{gn}}_G}\end{array}\right]}^T,$ ${\mathrm{\θ}}_{\mathrm{UG}}={\left[\begin{array}{cccc}{\mathrm{\θ}}_{\mathrm{Ug}1}& {\mathrm{\θ}}_{\mathrm{Ug}2}& \·\·\·& {\mathrm{\θ}}_{{\mathrm{Ugn}}_G}\end{array}\right]}^T;$

每台发电机的暂态电抗后电势及当前功角按照如下思路确定：

假设发电机节点i的功率注入为 ${\stackrel{*}{S}}_{\mathrm{gi}}=P_{\mathrm{gi}}+j{Q}_{\mathrm{gi}},$ 同时机端电压模制和相角为U_gi，δ_Ugi，那么首先求出此节点的等效注入电流 ${\stackrel{*}{I}}_{\mathrm{gi}}=\left(\frac{S_{\mathrm{gi}}}{U_{\mathrm{gi}}\∠{\mathrm{\θ}}_{\mathrm{gi}}}\right),$ 那么： $E_{\mathrm{di}}\∠{\mathrm{\δ}}_i=U_{\mathrm{gi}}\∠{\mathrm{\θ}}_{\mathrm{Ugi}}+{\stackrel{*}{I}}_{\mathrm{gi}}*\left({\mathrm{jx}}_{\mathrm{di}}\right),$ 下文中某一个变量有上标星号(*)表示此变量是一个复数矢量；

对于系统网络的处理，对于发电机和负荷分别进行简化处理；对于发电机，经典二阶模型下的发电机相当于有内阻的电压源(内阻是x′_di，而电压是E_di)，那么经过诺顿等值将这个电压源转换成一个电流源：阻抗仍为x′_di，

$I_{\mathrm{gi}}=\frac{E_{\mathrm{di}}}{x_{\mathrm{di}}^{\′}},$ 对于负荷节点k进行如下处理：假设当前潮流下，节点的电压模值为V_{Load_k}，而负荷有功和无功分别为P_k和Q_k，那么等效负荷阻抗为 $Z_{\mathrm{Load}\_k}={\left(\frac{V_{\mathrm{Load}\_k}^2}{P_k+{\mathrm{jQ}}_k}\right)}^{*}=R_{\mathrm{Load}\_k}+j{X}_{\mathrm{Load}\_k},$ 相应得负荷等效导纳为 $Y_{\mathrm{Load}\_k}=\frac{1}{Z_{\mathrm{Load}\_k}}=G_{\mathrm{Load}\_k}+j{B}_{\mathrm{Load}\_k};$

把负荷节点的等效阻抗以及发电机的暂态电抗加以考虑之后，系统的节点导纳矩阵变为：Y＝G_Y+jB_Y；其中：Y、G_Y、B_Y都是n×n维的矩阵；令 $\stackrel{*}{U}=U_x+j{U}_y$ 表示个节点电压列向量，每个元素为各个节点的电压矢量； $\stackrel{*}{I}=I_x+j{I}_y$ 代表各个节点的注入电流向量，下文中，某个列向量上标是一个星号(*)都表示此列向量每个元素都是复数矢量；其中：U_x为各个节点电压实部组成的向量，U_y为各个节点电压虚部组成的向量，I_x为各个节点注入电流实部组成的向量，I_y为个节点注入电流虚部组成的向量；

经过上述简化之后，只有发电机节点具有注入电流，其他的节点没有电流注入；

令 $Y^{\′}=\left[\begin{array}{cc}{G}_Y& -B_Y\\ B_Y& G_Y\end{array}\right],$ $U^{\′}=\left[\begin{array}{c}{U}_x\\ U_y\end{array}\right],$ $I^{\′}=\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right],$ 那么有： $I^{\′}=\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right]=Y^{\′}U=\left[\begin{array}{cc}{G}_Y& -B_Y\\ B_Y& G_Y\end{array}\right]\left[\begin{array}{c}{U}_x\\ U_y\end{array}\right]$

令系统所有节点的电压模值向量为U＝[U₁ U₂ ... U_N]，相角向量为θ_U＝{θ_u1 θ_u2 ... θ_uN}^T，N为系统所有节点的个数，那么：

U_x＝[U₁cosθ_u1 U₂cosθ_u2 ... U_ncosθ_uN]^T，U_y＝[U₁sinθ_u1 U₂sinθ_u2 ... U_nsinθ_uN]^T

令： $Z^{\′}=Y^{\′-1}={\left[\begin{array}{cc}{G}_Y& -B_Y\\ B_Y& G_Y\end{array}\right]}^{-1}$ 那么有：

$\left[\begin{array}{c}{U}_x\\ U_y\end{array}\right]=Z\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right]---\left(2\right)$

I_x，I_y都只有发电机节点对应元素非零，其他元素均为零；令：

$I_{\mathrm{Gx}}={\left[\begin{array}{cccc}{I}_{G1}\sin {\mathrm{\δ}}_1& I_{G2}\sin {\mathrm{\δ}}_2& \·\·\·& I_{{\mathrm{Gn}}_G}\sin {\mathrm{\δ}}_{n_G}\end{array}\right]}^T$ (1，...，n_G表示发电机序号)

$I_{\mathrm{Gy}}={\left[-,\begin{array}{cccc}{I}_{G1}\cos {\mathrm{\δ}}_1& I_{G2}\cos {\mathrm{\δ}}_2& \·\·\·& I_{{\mathrm{Gn}}_G}\cos {\mathrm{\δ}}_{n_G}\end{array}\right]}^T$ (1，...，n_G表示发电机序号)

那么I_x，I_y可以用以下的形式表示：I_x＝T×I_Gx，I_y＝T×I_Gy

其中：T是一个N×n_G维的矩阵；T的每一行对应系统的一个节点(编号为k)，如果本节点没有发电机接入，即非发电机节点，那么本行元素全为0，即T(k，：)＝[0...0]；如果本节点有发电机接入，即为发电机节点，假设此发电机序号为i，那么T(k，i)＝1，本行其他元素均为0，即T(k，：)＝[0...1...0]；

对(2)式两边线性化得：

$\left[\begin{array}{cc}\cos {\mathrm{\θ}}_U& -U\sin {\mathrm{\θ}}_U\\ \sin {\mathrm{\θ}}_U& U\cos {\mathrm{\θ}}_U\end{array}\right]\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_U\end{array}\right]=Z\left[\begin{array}{c}{I}_{\mathrm{x\δ}}\\ I_{\mathrm{y\δ}}\end{array}\right]\mathrm{\Δ\δ}---\left(3\right)$

其中：

ΔU，Δθ_U，Δδ分别为U，θ_U，δ的偏列差量，N和n_G分别为系统所有节点的个数以及发电机节点的个数，表达式中的各个矩阵的形式如下：

cosθ_U＝diag{cosθ_U1 cosθ_U2...cosθ_UN}；sinθ_U＝diag{sinθ_U1 sinθ_U2...sinθ_UN}；

Ucosθ_U＝diag{U₁cosθ_U1 U₂cosθ_U2...U_Ncosθ_UN}；

Usinθ_U＝diag{U₁sinθ_U1 U₂sinθ_U2...U_Nsinθ_UN}；

$I_{\mathrm{x\δ}}=\frac{{\∂I}_x}{\∂\mathrm{\δ}},$ $I_{\mathrm{y\δ}}=\frac{\∂I_y}{\∂\mathrm{\δ}}$ 均为N×n_G维的矩阵，可知 $I_{\mathrm{x\δ}}=\frac{\∂I_x}{\∂\mathrm{\δ}}=T\×\frac{{\∂I}_{\mathrm{Gx}}}{\∂\mathrm{\δ}},$ $I_{\mathrm{y\δ}}=\frac{\∂I_y}{\∂\mathrm{\δ}}=T\×\frac{{\∂I}_{\mathrm{Gy}}}{\∂\mathrm{\δ}},$ 而：

$\frac{{\∂I}_{\mathrm{Gx}}}{\∂\mathrm{\δ}}=\mathrm{diag}\left\{\begin{array}{cccc}{I}_{G1}\cos {\mathrm{\δ}}_1& I_{G2}\cos {\mathrm{\δ}}_2& ...& I_{\mathrm{GnG}}\cos {\mathrm{\δ}}_{n_G}\end{array}\right\}$

$\frac{{\∂I}_{\mathrm{Gy}}}{\∂\mathrm{\δ}}=\mathrm{diag}\left\{\begin{array}{cccc}{I}_{G1}\sin {\mathrm{\δ}}_1& I_{G2}\sin {\mathrm{\δ}}_2& ...& I_{\mathrm{GnG}}\sin {\mathrm{\δ}}_{n_G}\end{array}\right\}$

令： $D=\left[\begin{array}{cc}\cos {\mathrm{\θ}}_U& -U\sin {\mathrm{\θ}}_U\\ \sin {\mathrm{\θ}}_U& U\cos {\mathrm{\θ}}_U\end{array}\right],$ $E=\left[\begin{array}{c}{I}_{\mathrm{x\δ}}\\ I_{\mathrm{y\δ}}\end{array}\right],$ E是2N×n_G维的，那么：

$\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_U\end{array}\right]=D^{-1}\mathrm{ZE\Δ\δ}---\left(4\right)$

对于式(1)所示的电力系统机电振荡模型在平衡点处进行线性化，可得：

$\frac{d^2\mathrm{\Δ\δ}}{{\mathrm{dt}}^2}=B\left[\begin{array}{c}\mathrm{\Δ}{U}_G\\ \mathrm{\Δ}{\mathrm{\θ}}_{\mathrm{UG}}\end{array}\right]+\mathrm{C\Δ\δ}---\left(5\right)$

其中：ΔU_G，Δθ_UG，Δδ为U_G，θ_UG，δ的偏差量， $B={\left(2H\right)}^{-1}\left[\begin{array}{cc}-\frac{\∂g}{\∂U_G}& 0\\ 0& -\frac{\∂f}{\∂{\mathrm{\θ}}_{\mathrm{UG}}}\end{array}\right],$ $C=-{\left(2H\right)}^{-1}\frac{\∂f}{\∂\mathrm{\δ}}$ 有：ΔU_G＝T₂×ΔU，Δθ_UG＝T₂×Δθ_U，代入(5)式有：

$\frac{d^2\mathrm{\Δ\δ}}{{\mathrm{dt}}^2}=B\left[\begin{array}{c}\mathrm{\Δ}{U}_G\\ \mathrm{\Δ}{\mathrm{\θ}}_{\mathrm{uG}}\end{array}\right]+\mathrm{C\Δ\δ}=B\left[\begin{array}{cc}{T}_2& 0\\ 0& T_2\end{array}\right]\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_u\end{array}\right]+\mathrm{C\Δ\δ}={\mathrm{BT}}_2^{\′}\left[\begin{array}{c}\mathrm{\ΔU}\\ \mathrm{\Δ}{\mathrm{\θ}}_U\end{array}\right]+\mathrm{C\Δ\δ}---\left(6\right)$

其中： $T_2^{\′}=\left[\begin{array}{cc}{T}_2& 0\\ 0& T_2\end{array}\right],$ T₂＝T^T

带入(4)式就得到了多机系统线性化模型为：

$\frac{d^2\mathrm{\Δ\δ}}{{\mathrm{dt}}^2}=({\mathrm{BT}}_2^{\′}{D}^{-1}\mathrm{ZE}+C)\mathrm{\Δ\δ}=\mathrm{A\Δ\δ}---\left(7\right)$

步骤(4)求取矩阵A的特征值并选择慢模式相关特征值组；

4.1)采用QR算法求取A的全部特征值；

4.2)把A的特征值按照绝对值由小到大的顺序进行排序；

假设排序后的特征值顺序如下：{λ₁ λ₂ ... 

}，其中n_G为发电机的台数，也正是A的维数；值得注意的是λ₁＝0；

4.3)选取慢模式相关特征值组，关键是确定选取慢模式相关特征值的个数r；

按照如下的方法选择r：比较所有的λ_i/λ_i+1，i＝2，3，...，n_G-1，令： ${\mathrm{\λ}}_r/{\mathrm{\λ}}_{r+1}=\underset{i=\mathrm{2,3},..n_G-1}{\min }{\mathrm{\λ}}_i/{\mathrm{\λ}}_{i+1},$ 这样就确定了慢模式相关特征值的个数r；这种选择的依据是电力系统多时间尺度特性原理，λ_r/λ_r+1的值越小表示系统的多时间尺度特性越明显，系统的区域性也越明显，系统也越适合解列；经过这一步骤，就确定了慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}；

步骤(5)求取慢模式相关特征值组对于每条线路参数变化灵敏度；

线路就是指电力系统中实际的输电线路。首先对于所有节点进行编号：1，2，...，N，再对线路进行编号：1，2，...，Nl。一条线路主要有如下元素特征：起始节点编号、终了节点编号、线路电阻r、线路电抗x，线路电导g和线路电纳b，值得注意的是在选取线路起止节点时，哪一端点作为起始节点时无所谓的，所以在确定起始节点时随便指定其中一个端点作为起始节点即可。

某条线路的编号是L，其起始节点编号为I，终了节点编号为J，线路的电阻和电抗分别为r和x，电导和电纳分别为g和b，求取慢模式相关特征值组对于线路L参数变化的灵敏度，可按以下两步进行：

5.1)求取慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}中每个特征值对线路L参数变化的灵敏度，表示为{S_λ1 ^L S_λ2 ^L...S_λr ^L}，S_λi ^L表示特征值λ_i对于线路L参数变化的灵敏度；

特征值λ_i对于线路L参数变化的灵敏度可按照以下思路进行：

λ_i为A的一个特征值，λ_i对应的左右特征向量分别为φ_i(n_G×1)和ψ_i(1×n_G)；

那么λ_i对于线路L的导纳b变化的灵敏度表达式如下：

其中： $N^{,}=\left[\begin{array}{cc}0& M^T\\ M^T& 0\end{array}\right],$ $M^{\′}=\left[\begin{array}{cc}M& 0\\ 0& -M\end{array}\right],$ M＝[...1...-1...]^T，其中1的位置是I，而-1的位置是J；

λ_i对于线路L的导纳g变化的灵敏度表达式如下：

其中： $N^{\′\′}=\left[\begin{array}{cc}0& M^T\\ M^T& 0\end{array}\right],$ $M^{\′\′}=\left[\begin{array}{cc}M& 0\\ 0& M\end{array}\right];$

假设线路参数变化时，r和x同比例变化，那么由于 $g+\mathrm{jb}=\frac{1}{r+\mathrm{jx}}=\frac{r-\mathrm{jx}}{r^2+x^2},$ 线路的g和b是按照-r/x的比例变化的，以b的变化为基准(b＞＞g)，总结出λ_i对于线路L参数变化的灵敏度变化表达式如下：

$S_{\mathrm{\λi}}^L=-S_{\mathrm{\λi}}^{\mathrm{Lb}}+\frac{r}{x}{S}_{\mathrm{\λi}}^{\mathrm{Lg}}---\left(10\right)$

5.2)确定慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}对于线路L参数变化的灵敏度，按下式确定：

$S_{\mathrm{\σr}}^L=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{S}_{\mathrm{\λi}}^L---\left(11\right)$

S_σr ^L即为慢模式相关特征值组σ_r＝{λ₁，λ₂，..λ_r}对于线路L参数变化的灵敏度，用这个量来衡量线路参数变化对于慢模式组的影响；

经过这一步骤就能够得到所有线路的慢模式相关特征值组对于系统所有线路的参数的灵敏度信息；

步骤(6)根据步骤5得到的慢模式相关特征值组对于线路参数变化的灵敏度信息来确定弱连接；

6.1)把慢模式相关特征值组对于系统所有的Nl条线路的灵敏度按照绝对值从大到小的顺序进行排序，排序后的顺序为{S_σr ¹S_σr ²...S_σr ^Nl}；

6.2)选取弱连接线路，基本思路就是选取灵敏度比较大的一组线路作为系统弱连接线路；至于弱连接线路的条数wl可以由用户自己按照需要设定，也可以自动生成，本发明推荐一种自动生成wl的方法；对于一个多时间尺度特性比较明显的系统，弱连接线路与非弱连接线路的灵敏度数值差异会比较大；这种方法首先对于系统所有线路的灵敏度进行从大到小的排序，然后按照上述原理来选择弱连接线路的条数wl，选取wl的公式表达形式如下：

$\mathrm{abs}(S_{\mathrm{\σr}}^{\mathrm{wl}}/S_{\mathrm{\σr}}^{\mathrm{wl}+1})=\underset{i=\mathrm{1,2},..\mathrm{Nl}}{\max }\mathrm{abs}(S_{\mathrm{\σr}}^i/S_{\mathrm{\σr}}^{i+1})---\left(12\right)$

其中，符号abs表示取一个实数的绝对值。

Images