基于s域节点导纳矩阵的谐振稳定性分析方法

2023-11-11徐政

电力自动化设备 2023年10期

徐政

（浙江大学电气工程学院，浙江杭州 310027）

0 引言

随着新型电力系统建设的深入，电力系统在源-网-荷-储4 个层面逐渐表现出电力电子化特征，使得以往较少出现的谐振不稳定问题变得越来越普遍。国内外已有大量文献报道了发生谐振不稳定现象的案例［1-3］，包括交流电网中的谐振不稳定现象［1-2］和直流电网中的谐振不稳定现象［3］。

当发生谐振不稳定时，电力系统中该谐振频率下的电流和电压会很大，有可能造成设备损坏和保护跳闸［1］；当谐振频率落在次同步频段时，谐振频率电流流过汽轮发电机的定子绕组后，会在发电机的转子上产生振荡频率与此谐振频率互补的转矩，从而有可能激发汽轮发电机的固有扭振模态［4-5］，造成以往只在汽轮发电机经串联补偿线路送出系统中才会出现的机网复合共振现象［4-5］。随着电力系统电力电子化程度的不断加深，谐振稳定性问题日益突出，并将成为继同步稳定性、电压稳定性和频率稳定性之后的第四大电力系统稳定性问题。电力系统谐振稳定性分析将与经典的电力系统三大计算，即潮流计算、暂态稳定计算和短路电流计算具有同等重要的地位。因此，迫切需要有效的理论和工具对电力系统的谐振稳定性进行分析，并在此基础上提出改善谐振稳定性的措施。

然而，关于电力系统谐振稳定性分析的工具还比较匮乏，缺乏能够应用于大系统层面的谐振稳定性分析工具。目前业界普遍使用基于Nyquist 判据的方法来分析谐振稳定性问题［6-9］。这种方法将整个系统分为待接入装置和既有系统2 个部分，分别求出待接入装置的增量阻抗频率特性和既有系统的增量阻抗频率特性，然后基于Nyquist 判据，判断装置接入后整个系统是否会发生谐振不稳定问题。这种方法对单个装置接入系统的谐振稳定性分析比较方便，但难以提供谐振模态的精确信息，不适用于大系统层面的谐振稳定性分析。这是因为对于一般性电力系统而言，谐振稳定性分析需要确定该系统存在哪种谐振模态、各谐振模态的阻尼和谐振频率、各谐振模态的节点电压振型和参与因子、各谐振模态对特定元件参数的灵敏度等。

理论上，能满足上述大系统层面谐振稳定性分析需求的方法主要有状态空间法［10］和本文将要阐述的s域节点导纳矩阵法［10-14］。采用状态空间法分析谐振稳定性主要存在如下困难：①当考虑输电线路等分布参数元件时，因描述分布参数元件特性的方程是偏微分方程，整个电力网络已不能用线性时不变网络的标准状态空间模型来描述；②进一步考虑元件参数随频率变化的特性，即使对于由集总参数元件构成的电力网络，也无法用线性时不变网络的标准状态空间模型来描述；③对于电力电子装置，不易建立其在某个工作点上的增量线性化状态空间模型。

s域节点导纳矩阵法处理分布参数元件和频变参数元件不存在任何困难，故具有更好的适用性。但应用s域节点导纳矩阵法的主要问题是如何高效求解s域节点导纳矩阵行列式的根［11-14］。以往提出的方法还不能满足大系统层面电力系统谐振稳定性分析的要求［11-15］，本文将对如何解决此问题进行详细阐述。

1 谐振稳定性的定义和物理机理

关于谐振稳定性，目前业界并没有统一的定义，笔者将其定义为电力系统中以“固有谐振频率”振荡的自由分量的衰减特性［16］。当电力系统遭受扰动后，必然进入电磁暂态振荡过程，其电压、电流响应中除了具有基波频率的强制分量外，还包含以“固有谐振频率”振荡的自由分量。如果所有以“固有谐振频率”振荡的自由分量都是衰减的，则称电力系统是谐振稳定的，否则就称电力系统是谐振不稳定的。下面以图1 所示三调谐滤波器合闸到理想电源的简单系统为例说明谐振稳定性的概念。图中：us、is分别为电源电压、流过三调谐滤波器的电流；Us为电源电压幅值；ω0为电源角频率；C1—C3、L1—L3、R1—R3分别为三调谐滤波器的电容、电感、电阻；K 为电源开关；N1—N3、Ns为节点。

图1 三调谐滤波器合闸到电源的简单系统Fig.1 Simple system with triple tuned filter closing to power source

电源开关K 右侧为1 个三调谐滤波器［17］，其由1 个RL串联电路和2个RLC并联电路串联而成。当电源开关K 合上时，图1 所示的简单系统进入电磁暂态振荡过程，其特性可以由流过三调谐滤波器的电流is来呈现。is的表达式为：

式中：is，forced、is，free分别为is的强制分量部分和自由分量部分。强制分量部分仅包含1 个分量Ismcos(ω0tφzfilter)，其角频率与电源us的角频率相同，φzfilter为三调谐滤波器在角频率ω0下的阻抗角，幅值Ism=Us/Zfilter，Zfilter为三调谐滤波器在角频率ω0下的阻抗模值。另一方面，自由分量部分包含3 个分量Ires1eσres1cos(ωres1t+φres1) —Ires3eσres3cos(ωres3t+φres3)，分别与三调谐滤波器的3 种固有谐振模态相对应，分别为σres1±jωres1、σres2±jωres2、σres3±jωres3（σres1—σres3、ωres1—ωres3分别为3 种固有谐振模态的实部、虚部），而幅值Ires1—Ires3、阻抗角φres1—φres3由三调谐滤波器的参数和初始条件决定。对应每种谐振模态，其分别由实部和虚部组成；其中实部与该分量的衰减特性相对应，虚部与该分量的振荡频率相对应。

由式（1）可知：当三调谐滤波器的3 种固有谐振模态的实部σres1—σres3为负时，式（1）中的自由分量是衰减的，此时称图1 所示的简单系统是谐振稳定的；若3 种固有谐振模态中存在实部为正的谐振模态，则意味着is的3个自由分量中存在不衰减的自由分量，此时就称图1 所示的简单系统是谐振不稳定的。

值得指出的是，以往在进行电力系统机电暂态过程分析时，假定电力网络中以“固有谐振频率”振荡的自由分量是迅速衰减的。也就是说，假定在机电暂态振荡的时间尺度内，电压、电流中以“固有谐振频率”振荡的自由分量已衰减到0［5］，即电力系统不仅是谐振稳定的，而且其谐振模态具有很强的阻尼。基于上述假设条件，在进行电力系统机电暂态过程分析时不考虑电力网络本身的电磁暂态过程，电力网络仅被看作用于功率传递的静态元件，从而采用正序基频阻抗和代数方程进行描述，网络物理量为正序基频相量［5］。

通常情况下，固有谐振模态的实部是负的，但若在某些频段的电力系统中部分元件的增量阻抗呈现为负电阻特性（如同步发电机在次同步频段会呈现出负电阻效应［4-5］，两电平换流器、双馈风机等在一定的频段会呈现出负电阻效应［18-19］），那么固有谐振模态的实部可能由负变正，导致扰动后电压、电流响应中以固有谐振频率振荡的自由分量不会衰减，系统进入谐振不稳定状态。

2 s域节点导纳矩阵法的理论基础

1969 年，N.Balabanian 等人在《电网络理论》的专著中［20］给出了这样一个定理：对于线性时不变网络，其s域回路阻抗矩阵Zloop(s)的行列式det (Zloop(s))、s域节点导纳矩阵Ynode(s)的行列式det (Ynode(s))和系统状态方程的特征多项式det(sI-A)三者之间具有相同的非零值零点。

下面基于图2 所示的RLC 串联电路，对上述定理所给出的三者之间的关系进行验证。图中：R为电阻；iL为流经电感L的电流；uC为电容器C两端电压。

图2 RLC串联电路Fig.2 RLC series circuit

图2所示RLC电路的状态空间方程为：

式中：A、b为系数矩阵；x为状态变量；u为输入变量。可以求出系统的特征方程为：

式中：I为单位矩阵。若采用s域节点导纳矩阵法，则该电路的s域节点导纳矩阵为：

从而有：

若采用s域回路阻抗矩阵法，则该电路的s域回路阻抗矩阵为：

从而有：

显然，在s≠0的条件下，式（3）、式（5）和式（7）具有相同的零点。这样，状态方程的特征根计算可以转化为s域回路阻抗矩阵行列式det(Zloop(s))或s域节点导纳矩阵行列式det(Ynode(s))的零点计算。

在基于数字计算机的电网络分析方法中，节点导纳矩阵分析法极易实现［21］。因此本文将基于s域节点导纳矩阵Ynode(s)来求解系统的特征根，进而分析系统的谐振稳定性。状态空间分析法［10］与s域节点导纳矩阵分析法［11-15］的优势比较如表1 所示。由表可知，在大系统层面进行谐振稳定性分析时，s域节点导纳矩阵法具有更强的适应性。

表1 2种大系统层面谐振稳定性分析方法比较Table 1 Comparison of two resonance stability analysis methods for large scale system

3 基于s域节点导纳矩阵的谐振模态两阶段计算方法

3.1 电力系统谐振不稳定的基本表现形式及其内在特性

根据线性系统的基本理论，1 种谐振模态与特征根中的1 对共轭复根相对应，其中该共轭复根的实部表示该谐振模态的阻尼，该共轭复根的虚部表示该谐振模态的谐振频率。电力系统中所有谐振模态的阻尼特性决定于元件的电阻（电导）特性。由于谐振不稳定是由电力系统元件在某些频段的负电阻特性引起的，但这种负电阻特性所产生的能量通常不足以引起电气量的单调发散，因此，谐振不稳定一般表现为电气量的振荡发散。当研究电力系统的谐振不稳定时，可以排除由实数特征根引起的不稳定情形，仅关注由成对共轭复根引起的不稳定情形即可。在以下的det(Ynode(s))零点分析过程中本文仅关注det(Ynode(s))的共轭复根。

一般情况下，电力系统元件的电阻特性并不会对谐振模态的谐振频率产生明显影响。即在电力系统中，谐振模态所对应特征根的实部和虚部分别与电力系统中的电阻元件和电抗元件相对应，两者之间存在一定程度的解耦关系。

3.2 基于s 域节点导纳矩阵计算谐振模态的总体思路

根据电力系统中阻尼与谐振频率分别与电力系统中的电阻和电抗相对应的解耦特性，基于s域节点导纳矩阵计算谐振模态的数值方法分为2 个阶段。①第1 阶段，在无阻尼系统中进行分析：此时电力系统所有元件中的电阻（电导）性参数置0，计算谐振模态的无阻尼谐振频率，从而确定系统在所研究频段内的谐振模态数目，再进一步确定每个谐振模态的节点电压振型和节点参与因子。②第2 阶段，在考虑所有阻尼的完整系统中进行分析：此时考虑所有元件中的电阻（电导）特性，基于测试信号法［5］精确计算谐振模态的阻尼和谐振频率，从而判断系统的稳定性，必要时可以进一步计算各谐振模态对特定元件参数的灵敏度等。

在以下的分析中，采用图1 所示的简单系统具体展示各阶段的实现过程。图1 中三调谐滤波器的调谐次数分别设置为3 次、24 次和37 次，具体参数如表2所示［17］。

表2 三调谐滤波器参数Table 2 Parameters of triple tuned filter

需注意列写s域节点导纳矩阵时，电压源按短路处理，电流源按开路处理，此时图1 所示系统中的节点Ns变为参考电压节点，图1 所示系统就变为一个3 节点系统。设节点N1—N3分别对应1 — 3 号节点，则图1所示系统的s域节点导纳矩阵为：

4 在无阻尼系统中实现第1阶段算法的过程

电力系统中所有元件中的电阻（电导）性参数置0 后，重点分析谐振模态结构、s域节点导纳矩阵结构。

4.1 无阻尼系统的谐振模态结构

在有阻尼的电力系统中，与第i种谐振模态相对应的1对共轭复根的一般形式为：

式中：si为第i种谐振模态的特征根，σi为si的实部，ωi为si的虚部；s*i为si的共轭；fi为第i种谐振模态的谐振频率。

而在无阻尼的电力系统中，代表谐振模态阻尼的系统特征根实部必然为0，这样si就变为1 个纯虚数，即：

当si取纯虚数时，Ynode(si)的计算就退化为正弦稳态分析时的常规节点导纳矩阵计算。目前正弦稳态分析时的常规节点导纳矩阵计算已有较成熟的模型和算法［21］。

4.2 无阻尼系统的s域节点导纳矩阵结构

在系统无阻尼的条件下，系统特征根为纯虚数。设si=jωi=j2πfi，则si对应的Ynode(si)变为：

式中：Ynode(jωi)为常规正弦稳态分析中的节点导纳矩阵；Gnode(jωi)为Ynode(jωi)的实部，表示节点电导矩阵；Bnode(jωi)为Ynode(jωi)的虚部，表示节点电纳矩阵。由于已设系统的阻尼为0，则节点电导矩阵Gnode(jωi)=0，从而有Ynode(si)=jBnode(jωi)，其中Bnode(jωi)为实对称矩阵。

4.3 无阻尼系统谐振频率的计算方法

根据式（11），由于si为系统的特征根，从而满足如下关系：

因而可以通过求解det(Bnode(jωi))=0来得到ωi。即可得到无阻尼谐振频率fi的计算步骤如下：

1）设谐振稳定性分析所关注的频段为[fst，fend]，fst、fend分别为起始频率和终止频率；

2）在[fst，fend] 频段内以一定的步长计算det(Bnode(j2πfk))，这里fk为离散频率点；

3）检查相邻2 个离散频率点之间det(Bnode(jω))是否存在过零点，本文中根据det(Bnode(j2πfk-1))×det(Bnode(j2πfk))的计算结果进行判断，若此值大于0，则表示没有过零点，若此值小于0，则表示存在过零点；

4）若存在过零点，则通过2 点线性插值求出过零点频率fzero。

4.4 谐振模态的节点电压振型与最佳可控可观测节点

文献［11］对谐振模态的节点电压振型与节点参与因子进行了分析，得到如下结果：

1）Bnode(jωi)的特征值λ1=0；

2）与零特征值相对应的右特征向量M1被称为节点电压振型，其标示了在谐振模态jωi下各节点电压的相对振幅大小和相对相位关系；

3）节点电压振型中振幅最大的节点既是谐振模态jωi的最佳可观测节点，又是谐振模态jωi的最佳可控制节点，后文中将这个节点称为最佳可控可观测节点。

下面讨论与零特征值相对应的右特征向量M1的数值计算方法。推荐的计算步骤如下。

1）对Bnode(jωi)进行QR 分解，得到Bnode(jωi)=QR，这里Q为正交矩阵，R为上三角矩阵。由于Bnode(jωi)具有一个零特征值，因此矩阵R为：

2）与零特征值相对应的右特征向量M1≠0，且满足如下关系：

即：

3）由于特征向量M1不是唯一的，不妨设mn=1。这样M1中的其他元素可以通过求解如下方程得到：

求解方程式（16）可得特征向量M1，找出M1中的最大元素序号，即可确定与谐振模态jωi对应的最佳绝对值可控可观测节点。

4.5 第1阶段算法的具体实施示例

以图1 所示简单系统为例展示第1 阶段算法计算过程。构造无阻尼系统的方法是将图1 系统中的所有电阻（电导）性参数置0，具体做法为：对于2 个节点之间唯一的电阻性元件将其电阻值置0；对于2 个节点之间存在其他元件的电阻性元件将其电导值置0。

假设谐振稳定性分析所关注的频段[fst，fend]=[50，2 500] Hz，在此频段内以1 Hz 为步长扫描计算det(Bnode(j2πfk))，然后检查相邻2 个离散频率点之间det(Bnode(jω))是否存在过零点，若存在过零点，则通过2点线性插值求出过零点频率fzero。

计算结果得到3个过零点如下：第1种谐振模态无阻尼谐振的过零点频率fzero1=150.0 Hz；第2 种谐振模态无阻尼谐振的过零点频率fzero2=1 200.1 Hz；第3 种谐振模态无阻尼谐振的过零点频率fzero3=1 850.0 Hz。

与第1 种谐振模态频率fzero1=150.0 Hz 对应的节点电压振型Vms1=[1.000 0 0.988 7 0.002 3]T，与第2 种谐振模态频率fzero2=1 200.1 Hz 对应的节点电压振型Vms2=[1.000 0 0.277 5 0.498 2]T，与第3 种谐振模态频率fzero3=1 850.0 Hz 对应的节点电压振型Vms3=[-1.000 0 0.717 1 ]

0.498 9T。在3组节点电压振型中，第1 个元素的绝对值均是最大的，这表明1 号节点（N1）是最佳可控可观测节点。

5 在有阻尼的完整系统中实现第2阶段算法的过程

5.1 采用测试信号法的理论依据

考虑所有元件中的电阻（电导）特性后进行第2阶段的计算过程。s域下系统的节点电压方程为：

式中：Vnode(s)和Inode(s)分别为s域节点电压向量和s域节点注入电流向量。考察Inode(s)为单个输入时的情况，设节点k为最佳可控可观测节点，在节点k上注入电流ik(s)，计算与节点k相连接的某条（或某几条）分流支路（并联支路）的电流ishunt(s)。设分流支路的一端为节点k，另一端为节点l，且分流支路的导纳为yshunt(s)。则式（17）可变为：

式中：β为n×1 维列向量，其第k个元素为1，其余元素为0。而ishunt(s)可以表示为：

式中：γ为1×n维行向量，其第k个元素为1，第l个元素为-1，其余元素为0。定义分流比ηshunt为流过分流支路的电流ishunt(s)与注入电流ik(s)之比，则根据式（19）可得：

式中：Y*node(s)为Ynode(s)的伴随矩阵。显然，ηshunt(s)是一个标量，可以理解为是一种传递函数，其分母为系统的特征方程det(Ynode(s))。因此ηshunt(s)的分母多项式一定包含系统特征方程的信息。

一般情况下，ηshunt(s)的分子和分母都是s的高次多项式。但当s处于系统特定特征根的邻域中时，ηshunt(s)的分母可以用与该特征根对应的二次多项式来表示，而ηshunt(s)的分子则同样可用1 个二次多项式来逼近。这样，在系统特定特征根的邻域中，ηshunt(s)可以表示为：

式中：a0、a1、b0、b1、b2为待定的多项式系数，均为实数。

值得指出的是，在最佳可控可观测节点k注入电流信号的条件下，可以得到多种分母为系统特征方程的传递函数，比如节点k本身的策动点阻抗Zkk(s)［11］。但就传递函数在系统特征根邻域中进行2 阶近似的可适用性而言，分流比ηshunt优于策动点阻抗Zkk(s)。

5.2 测试信号法的具体实施示例

为了确定ηshunt(s)在特定特征根邻域中的近似表达式（21），可以让s沿着虚轴在j2πfzero的邻域内变化，然后再辨识出式（21）中的多项式系数。具体的实现步骤如下。

1）对应第1 阶段已经计算出的特定无阻尼谐振频率fzero，在±3 %fzero的频率范围内构造5个临近的频率计算点，分别为f1=0.97fzero、f2=0.985fzero、f3=fzero、f4=1.015fzero、f5=1.03fzero。

2）在与fzero相对应的最佳可控可观测节点k中注入单位电流相量。然后在上述的5 个频率点中分别计算出ηshunt(j2πf1)、ηshunt(j2πf2)、ηshunt(j2πf3)、ηshunt(j2πf4)和ηshunt(j2πf5)。

3）根据s分别取j2πf1、j2πf2、j2πf3、j2πf4和j2πf5时的ηshunt(j2πf1)、ηshunt(j2πf2)、ηshunt(j2πf3)、ηshunt(j2πf4)和ηshunt(j2πf5)，辨识式（21）中的a0、a1、b0、b1、b2这5个待定实系数。

4）求解式（21）分母的二次多项式，得到与特定无阻尼谐振过零点频率fzero对应的精确谐振模态特征根σres±jωres。

值得指出的是，上述±3 %fzero的频率范围和5 个临近的频率计算点并不绝对。但改变这2 组参数，并不会对结果有太大影响。比如，将±3 %fzero的频率范围改变为±2 %fzero或 ±4 %fzero，以及将5 个临近的频率计算点变为7个或9个，对结果的影响都不大。

仍然以图1 所示简单系统为例展示第2 阶段算法的过程。计算结果如附录A 表A1 所示。对上述计算结果的准确性进行检验，方法是验证所得到的特征根是否满足det(Ynode(s))=0 的条件。将上述特征根代入式（8）中，得到：

可见，第2 阶段计算得到的谐振模态特征根，在很高的精度上满足了det(Ynode(s))=0 的条件。

5.3 谐振模态阻尼的灵敏度分析

显然，特征根的实部决定于系统中的电阻（电导）性参数，通过改变系统中的特定电阻（电导）参数，可以计算特征根的阻尼对特定电阻（电导）参数的灵敏度。仍以图1 所示简单系统为例进行分析。分别改变图1 中并联电阻R1、R2的取值，观察3 个特征根的变化情况。计算结果分别见表3、4。

表3 谐振模态阻尼关于并联电阻R1的灵敏度Table 3 Sensitivity of resonance mode damping to parallel resistance R1

由表3 可知，谐振模态2 和谐振模态3 的阻尼值对R1的改变很敏感，特别是当R1取负值时，谐振模态2、3 变得不稳定。说明谐振模态2、3 的阻尼与R1强相关，而谐振模态1 的阻尼与R1几乎不相关。由表4 可知，谐振模态1 的阻尼值对R2的改变十分敏感，特别是当R2取负值时，谐振模态1 变得不稳定。说明谐振模态1 的阻尼与R2强相关，而谐振模态2和谐振模态3的阻尼与R2仅为弱相关。

表4 谐振模态阻尼关于并联电阻R2的灵敏度Table 4 Sensitivity of resonance mode damping to parallel resistance R2

另外，对于阻尼很强的谐振模态（如表3 中的谐振模态1），其5 个频率点上的分流比模值是单调变化的，与阻尼较弱的谐振模态2、3 存在峰值有所不同。当改变R2的取值，使得谐振模态1 变为弱阻尼时，对应谐振模态1的5个频率点上的分流比模值也会呈现峰值特征。R2=1 200 Ω时，谐振模态1的阻尼比为0.049 6，谐振频率为149.816 6 Hz。对应谐振模态1的5个频率点上的分流比分别为2.219 4∠ 144.44°、2.243 5∠132.54°、2.096 7∠118.91°、1.781 9∠ 106.31°、1.408 6∠97.40°。显然，第2 个频率点上的模值是这5个值中的峰值。

综上分析，对于弱阻尼的谐振模态，其有阻尼谐振频率与无阻尼谐振频率非常相近，分流比ηshunt的模值在有阻尼谐振频率点上会出现峰值。

另外，还需指出，表3、4 中，当设置电阻为负时，基于s域节点导纳矩阵的分析方法会出现不稳定的谐振模态，但该结果在基于电磁暂态仿真方法研究系统谐振稳定性的过程中难以复现。这是因为在这种系统条件下，电磁暂态仿真方法不可能建立初始工作点。这也是s域节点导纳矩阵法相较于电磁暂态仿真方法在分析系统谐振稳定性时的优势之一。

6 大系统层面谐振稳定性分析实例

考察1 个四端直流电网的谐振稳定性问题。该四端直流电网接线如附录A 图A1 所示。该电网为双极带金属回线的直流电网，金属回线全网单点接地，接地点设在换流站C。

设图A1 所示直流电网中所有平波电抗器的电感值均为100 mH；从直流侧看进去的每个模块化多电平换流器（modular multilevel converter，MMC）的阻抗频率特性如式（25）所示；三线耦合输电线路的单位长度电阻参数如式（26）所示，单位长度电感参数如式（27）所示，单位长度电容参数如式（28）所示。

式中：f为频率。由式（25）可见，MMC的阻抗频率特性在100 Hz范围内存在负电阻频段，因此图A1所示直流电网的谐振稳定性问题是一个需关注的问题。

对图A1所示直流电网进行谐振稳定性分析，计算结果如附录A表A2所示。由表可知，在100 Hz的频率范围内，该系统存在6 个谐振模态，尽管阻尼比小于5 % 的有3个，但这6个谐振模态都是稳定的。

基于此实际案例的分析，以下2点需特别指出：

1）输电线路具有多相耦合、分布参数和参数随频率变化的特点，这些特点限制了状态空间模型在谐振稳定性分析方面的有效应用，而s域节点导纳矩阵法可以完全克服这些困难；

2）即使MMC 的阻抗频率特性存在负电阻频段，谐振不稳定也并非是必然发生的现象，实际系统是否会出现谐振不稳定问题需要基于具体的系统参数进行精确计算。