朱永强，刘光晔，曹 斌，匡一雷，康志豪，刘媛媛

（1.湖南大学电气与信息工程学院，长沙 410082；2.国网湖南省电力公司岳阳供电分公司，岳阳414000；3.国网湖南省电力公司检修公司，长沙 410000）

电力系统电压稳定鞍结分岔点的快速求取方法

朱永强1，2，刘光晔1，曹 斌3，匡一雷1，康志豪1，刘媛媛1

（1.湖南大学电气与信息工程学院，长沙 410082；2.国网湖南省电力公司岳阳供电分公司，岳阳414000；3.国网湖南省电力公司检修公司，长沙 410000）

为快速求取电力系统电压稳定鞍结分岔点，提出一种自动变步长的连续潮流法和局部曲线拟合法相结合的方法。运用广义戴维南动态等值原理，证明非线性电力网络传输功率达到极值必要条件，定义阻抗模裕度指标，判定系统静态电压稳定性最薄弱节点。由当前运行点阻抗模预测下一运行点潮流状态，得到增长步长，实现自动变步长，快速逼近而不越过鞍结分岔点。当阻抗模裕度值达到满足局部曲线拟合要求的设定值时，满足阻抗模裕度判据，对鞍结分岔点精确拟合，得到较准确的鞍结分岔点值。多个IEEE测试系统的仿真表明，该算法简单、精度高、速度快、计算次数少，且适用于规模较大的电力网络。

电力系统；电压稳定；鞍结分岔；广义戴维南等值；阻抗模裕度判据；局部曲线拟合

Abstract:In order to fast calculate the saddle node bifurcation point of voltage stability of power system，a combination method with automatic variable step is proposed based on continuation power flow method and local curve fitting meth⁃od.By using generalized Thevenin dynamic equivalence theory，the essential conditions for reaching the extreme values of transmission power in a nonlinear power network are proved.Moreover，impedance modulus margin index is defined to determine the weakest node of static voltage stability of the system.An increasing step length is obtained according to the power flow state at the next operation point predicted by the impedance modulus at the current operation point，thus an automatic variable step is realized，and the saddle node bifurcation point is approached fast but not crossed over.When the impedance modulus margin satisfies the value required by local curve fitting，i.e.，impedance modulus mar⁃gin criterion，the saddle node bifurcation point will be accurately fitted and a more accurate value will be obtained.Simulations of several IEEE test systems show that the proposed method is simple，accurate，fast and with less computa⁃tional load，thus it can be applied to large-scale power networks.

Key words:power system；voltage stability；saddle node bifurcation；generalized Thevenin equivalence；impedance modulus margin criterion；local curve fitting

随着我国电力网络的日益扩大，电压稳定问题成为电网安全运行的主要研究方向[1]。电压稳定研究的重点是反映电压稳定性的安全指标的确定。裕度指标包含信息量大，且有较好的直观性和线性，受到电力研究者的高度重视[2]。负荷裕度是指在某种负荷功率增长方式下，基态运行点到电压极限点的距离，重点是求取鞍结分岔点SNBP（saddle node bifurcation point）[3]。

求取SNBP的方法主要有直接法[4]、非线性规划法[5]以及连续潮流法 CPF（continuation power flow）[6]。直接法计算量小，且可解决雅可比矩阵奇异问题，但提供有用信息极少，难以满足运行人员要求。非线性规划法可避免潮流不收敛，但计算量过大，难以满足大系统要求。CPF法通过增加参数方程，避免临界崩溃点潮流发散，但是通过固定步长逐步计算，获取鼻形曲线，将鼻尖作为SNBP，若步长过小，计算速度太慢，在线应用困难；若步长过大，极易越过极限点，造成迭代不收敛。

为了快速精确地求取SNBP，本文提出一种自动变步长的CPF法，快速逼近却不越过SNBP，并采用局部曲线拟合法对SNBP进行精确拟合。对多个IEEE测试系统的仿真计算表明，所提方法不仅有效，而且实用。

1 基于阻抗模裕度的系统薄弱节点分析

1.1 广义戴维南动态等值方法

文献[7]将一个复杂的线性系统通过戴维南等值原理等值为图1所示的两节点简单系统，并证明电力网络输送有功功率达到极值必要条件是负荷静态阻抗的模值与负荷节点戴维南等值阻抗的模值相等。本文提出了一种命名广义戴维南的动态等值方法，证明该结论在非线性电力网络中且负荷的功率因素变化的情况下也正确，且只需一个潮流点就可求解出各种等值参数。

图1 戴维南等值两节点系统Fig.1 Two-node system with Thevenin equivalence

式中：̇LD为负荷两端电压相量；为流过负荷电流相量。将非线性电力网络的广义戴维南动态等值阻抗定义为，即

由式（4）可得，电压实部和虚部分别对电流实部和虚部求偏导均可得到RTHEV和XTHEV。

负荷节点注入的复功率为

式中，PLD、QLD分别为负荷有功和无功功率。

考虑负荷的功率因数是变化的，则限制条件为

式中：k为负荷功率因数；C为常数。当C=0时，负荷的功率因数是恒定的；当k=0时，负荷的无功功率也是恒定的。

通过构造拉格朗日函数来求取PLD的极大值，其中γ为拉格朗日乘数，即

将式（4）和式（5）代入式（7），得到关于Ix、Iy的函数，分别对Ix、Iy求偏导，并令其等于零可得

又因为电压电流还存在关系

将式（4）和式（9）代入式（8）可得到关于Ix、Iy为变量的二元齐次线性方程组。因为电流不可能恒为零，所以系数矩阵行列式必为零，由此得到负荷的等值阻抗的模值等于系统的广义戴维南动态等值阻抗的模值，即

这样证明了负荷功率因素变化的非线性电力网络输送有功功率达到极值必要条件。

1.2 阻抗模裕度的定义

通过SNBP判断电压稳定临界崩溃点，重点研究电压稳定最薄弱节点的SNBP，参考文献[2]给出了阻抗模裕度的定义为

该指标反映系统正常工作点距离电压稳定崩溃点有“多远”，且指标越小电压稳定性越差。因此由该指标可以判断电压稳定性最薄弱节点，然后去单独求取该薄弱节点的SNBP。当ΔZ＞0，节点电压稳定；当ΔZ=0，电压临界稳定，发生鞍结分岔；当ΔZ＜0，节点电压不稳定。

2 变步长的CPF法

2.1 CPF法基本思想

CPF法的基础是由状态变量和可变参数构成的潮流方程，即

式中：x为状态变量，表示电压幅值及相角；λ为控制功率变化的参数，称为负荷因子。

将式（12）更具体表述为

式中：PG,i、QG,i为发电机的输入功率；PL,i、QL,i为负荷消耗的功率；Ps,i、Qs,i为节点的注入功率；g1为PV、PQ节点的并集；g2为PQ节点的集合。若不是发电机节点，则QG,i=0。当发电机节点无功越限，才有无功平衡方程，QG,i=QGmax或QGmin。上下限分别对应最大负荷与基态负荷。

CPF法有4个主要环节：参数化、预估、步长控制、校正[8]。从基态点(x0,λ0)开始，逐步增加负荷因子，求解下一个准确解，最终求得SNBP，具体计算如下。

1）参数化

为了防止在临近鼻形曲线鼻尖时系数矩阵奇异，以状态变量的某个分量为参变量，增加一个方程，通常采用局部参数法为

式中：xk为x的第k个分量；η为xk的预测值。

2）预估步

预估步通常采用切线法，即负荷功率控制参变量沿着切线方向增长，逐步预测下一个潮流解[9]。首先确定切线方向为

式中：dx、dλ表示切线方向；±1表示崩溃点前为正，之后为负；为潮流函数对状态变量和负荷因子的导数；表示与选取状态变量对应的第k个元素为1，其余为0。再预测下一个潮流点为

式中：λk为当前状态量；σ为预估的步长。

3）控制步长

常规CPF法的定步长是有极大缺陷的，步长过大能提高计算速度，但极易越过极限点造成潮流发散，步长过小则效率过低[10]。而自动地变步长，既能大大提高运算效率，又不会越限。

4）校正步

2.2 广义戴维南动态等值预测步长

随着负荷因子增长，对广义戴维南动态等值阻抗模值与负荷静态等值阻抗模值大小及变化趋势作简要分析，以图2所示的简单两节点系统为例，线路只计电抗X，超前的相位为δ，δ为功角。

图2 两节点系统结构Fig.2 Structure of two-node system

负荷节点的电压表示为

两边同时取负号，并对I求导得

若式（19）为定义的广义戴维南动态等值阻抗，则可将式（18）改写为

式中，K为式（19）中广义戴维南动态等值阻抗的模值，K越来越大，则越来越大。

以上分析表明，随着负荷因子λ逐渐增大，负荷静态等值阻抗的模值越来越小，负荷节点的广义戴维南动态等值阻抗的模值越来越大。可见在基态运行点相差最大，越接近SNBP，差值越小。基态运行时系统稳定，阻抗模裕度应大于0，由式（11）可得，为了使预测到的负荷静态等值阻抗模能快速地向广义戴维南动态等值阻抗模靠近，即迅速向SNBP靠近，令预测到的负荷静态等值阻抗模为上个运行状态两者的中值，即

由式（10）知，负荷节点获取极限功率必要条件是静态模与动态模相等，则由式（23）知，预测到的运行状态点既向极限点靠近，而不会越过极限点。

假设采用同步负荷功率扰动方式，详见式（13），负荷的功率因数也是恒定的，所以预测的负荷视在功率的增长倍数与预测的有功功率的增长倍数相等，那么预测到的负荷因子可表示为

式中：Pi,0、Qi,0分别为节点i基态运行的有功功率与无功功率；Si,k+1为预测到的负荷视在功率。

预测节点i的视在功率平方为

预测到的负荷静态等值阻抗模与上一个运行状态的静态等值阻抗模之差为

以上分析可知，在靠近基态点时，负荷的静态等值阻抗模值与广义戴维南动态等值阻抗的模值差值比较大，越靠近SNBP，差值越小。由式（26）知，在靠近基态点时，差值较大，预测到的负荷增长因子距离较大；靠近SNBP时，两者差值较小，预测到的负荷因子距离也较小。以上方案刚好满足CPF法自动变步长的要求。

将预测到的负荷因子λk+1与当前运行点的负荷因子λk求差值，并除以预测方向分量dλ，可求出预测步长σ，代入式（16）求出预测运行点，再按照CPF法逐步进行计算。

2.3 阻抗模裕度的计算过程

注入节点的功率可表示为

系统潮流方程表示为

式中：W=[Ps,1,Qs,1,…,Ps,n-1,U2(n-1)]T；Un、U分别为平衡节点与其他节点的电压，且电压是在直角坐标表示的，U=[e1,f1,…,en-1,fn-1]T。

注入负荷节点i电流方程为

负荷节点i静态等值阻抗为

式（29）两边同对λ求导得

式（28）两边同时对λ求导可得

由式（32）可得

式（31）和式（33）通过求导链式法则可得广义戴维南动态等值阻抗为

将式（30）和式（34）代入式（11）可求出阻抗模裕度值。

3 局部曲线拟合技术

若拟合点距离SNBP较远，则经过曲线拟合得到SNBP值与准确SNBP值存在很大误差[12]。一般拟合法至少需要两个运行点，本文采用一种局部曲线拟合法，通过广义戴维南动态等值计算的阻抗模来预测CPF的负荷增长因子，实现自动变步长，从而迅速地向SNBP靠近（见第2节）。当预测的运行点潮流收敛，且满足阻抗模裕度判据ΔZ≤ε（ε为很小的正数，保证此时运行点离SNBP很近），用该运行点就可拟合出SNBP，使拟合精度大大提高。

局部曲线拟合法的数学表达式为

式中：ΔUw、Δλ分别为电压稳定性薄弱节点满足要求的拟合点到SNBP的电压大小变化量和增长的负荷因子（步长）；α、β为拟合系数。

对式（35）两边同时对Δλ求导可得

因为λ=λ′+Δλ，λ′为当前状态的负荷因子，Δλ为增长负荷因子，λ为待求的负荷因子，所以dΔUw/dΔλ=dUw/dλ，代入式（36）得

对式（37）继续对Δλ求导可得

拟合点离SNBP很近，ΔUw=0，dUw/dλ通过式（33）可得，式（32）继续对λ求导得

dJ/dλ通过J中的电压对λ求导即可求取，由式（39）就可求出d2Uw/dλ2。

将ΔUw、dUw/dλ、d2Uw/dλ2代入式（36）和式（38）可求得拟合系数α、β。式（35）对ΔUw求导得

在SNBP，dλ/dΔUw=0，所以2αΔUw+β=0，可求出ΔUw。再将α、β、ΔUw代入式（35）可求出拟合点到SNBP的负荷因子增量Δλ，从而求出SN⁃BP值。

4 算法步骤

步骤1 初始化，输入原始数据，设置常数ε，并置负荷因子λ0=1，用牛顿法求出基态点潮流数据。

步骤2 潮流若收敛，则按第2.3节步骤计算出各PQ节点的负荷的静态等值阻抗模和广义戴维南动态等值阻抗模并计算出阻抗模裕度值ΔZ，将ΔZ最小节点确定为电压稳定最薄弱节点，假设该节点为j；若潮流不收敛，则返回步骤1。

步骤3 单独针对节点j，使用局部参数法参数化，并采用切线法得出步长预测的方向dx、dλ。

步骤4 采用广义戴维南动态等值原理，按照第2.2节方法，预测出下一运行点的负荷因子λ1，并按式(λ1-λ0)/dλ预测步长σ，并计算预测的下一运行点通过垂直校正法求出准确解(x1,λ1)。若不收敛，则采用局部校正法进行校正。

步骤6 采用第3节的局部曲线拟合方法，求出该运行点到SNBP的负荷增长因子Δλcr，继而求出极限负荷因子λcr=λ+Δλcr。因为拟合点离SN⁃BP很近，到SNBP的负荷增长因子很小，因此一般不会出现无功越限的情形，该预测点便是SNBP。

5 仿真计算与分析

采用同步功率扰动方式，若节点i为PQ节点，则在直角坐标下的潮流方程为

式中：Pi,0、Qi,0分别为基态下节点i注入的有功与无功；n为潮流计算次数。若节点i为PV节点，则潮流方程表示为

式中：Us,i为节点给出的电压幅值；ei、fi分别为要求的实际节点电压实部与虚部。

对多个IEEE系统（IEEE14、IEEE30、IEEE57、IEEE118）进行仿真计算，并考虑无功越限，若无功越限，则将PV节点改成PQ节点。

在基态(λ0=1)运行时，计算4个系统各PQ节点的阻抗模裕度值，并比较大小，确定阻抗模裕度最小的节点为系统静态电压稳定性最薄弱节点。4个系统薄弱节点分别为14、30、31、41，阻抗模裕度分别为0.723 1、0.688 0、0.625 4、0.787 0，重点是求这些薄弱节点的SNBP。

先以IEEE30和IEEE57节点系统为例，详细计算结果见表1和表2。其中λ0为当前运行点的负荷因子，λ1为预测的下一个运行点的负荷因子。

由表1可知，n在1～6时是自动变步长的连续潮流计算过程，可见在基态点附近（1～1.5），预测到的相邻两个运行点的负荷因子之差较大，即步长较大。随着负荷因子逐渐向SNBP（1.5以后）靠近，差值越来越小，即步长越来越小，满足理想的步长控制要求。阻抗模裕度随着负荷因子的增大而减小，当n=6时，ΔZ＜ε（ε=0.03 p.u.），该运行点离SNBP很近，满足局部曲线拟合的要求。n=7是拟合的结果，拟合点负荷因子λ=1.546 34，从拟合点到SN⁃BP增长负荷因子很小Δλcr=0.000 25，拟合后无PV节点无功越界，计算结束。只用了7次潮流计算就求取了SNBP的值。

表1 IEEE30节点系统计算结果Tab.1 Calculation results of IEEE 30-node system

表2 IEEE57节点系统计算结果Tab.2 Calculation results of IEEE 57-node system

由表2可知，对于节点数较多的IEEE57节点系统，前6次潮流计算也是向SNBP快速逼近的过程，也满足连续潮流计算的步长控制要求。当n=6时，ΔZ＜ε（ε=0.03 p.u.），满足拟合要求，拟合后也无PV节点无功越界，只用了7次潮流计算就求取了SNBP值。可见，计算次数不会随着系统规模的扩大而明显增加。

将所提方法和文献[13]提出的改进CPF法以及固定步长（取0.01）的常规CPF法计算结果进行比较，比较结果见表3（含IEEE14、IEEE118的结果）。其中NB、λB分别为改进CPF法的计算次数与SN⁃BP的值；NG、λG分别为所提方法总共计算次数与SNBP的值；NC、λC分别为常规的定步长连续潮流法计算次数与SNBP的值。

表3 3种方法的比较Tab.3 Comparison among three methods

由表3可知，本文所提方法求取的SNBP值与改进CPF法计算值最大相对误差为0.028%，与定步长常规CPF法相比，最大相对误差为0.028%。但本文所提方法的潮流计算次数比另外两种方法少很多，并且随着系统节点数的增多，计算次数变化不明显，计算稳定，而另外两种方法却不具备这个特点。可见，本文所提方法不仅计算精度高，而且计算量小，计算速度快，适用于多节点的电力网络。

6 结语

运用广义戴维南动态等值方法，提出了一种自动变步长的CPF法，预测步长在基态点附近的概率大，而在崩溃点附近的概率小，能快速地逼近而不越过SNBP。由于阻抗模裕度反映当前运行点到SNBP距离，提出了阻抗模裕度判据，由阻抗摸裕度判据判定满足拟合要求的运行点，进行精确局部拟合，可以求取较准确的SNBP。算法简单，只用到雅可比矩阵保留的因子表，多个IEEE测试系统的仿真表明，该法计算精度高，次数少，速度快，适用于规模较大的电力网络，具有一定的在线应用价值。

[1]李娟，黄喜旺，关程宇，等（Li Juan，Huang Xiwang，Guan Chengyu，et al）.SVC和IPC联合改善异步机风电场的电压稳定性（Improvement of voltage stability of asynchro⁃nous wind farms based on SVC and IPC）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），2016，28（8）：79-84．

[2]吴政球，李日波，钟浩，等（Wu Zhengqiu，Li Ribo，Zhong Hao，et al）.电力系统静态电压稳定极限及裕度计算综述（Summary of power system's static voltage stability lim⁃itation and load margin calculation）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），2010，22（1）：126-132．

[3]牟晓明，李志民（Mou Xiaoming，Li Zhimin）.一种计算电压稳定边界的两阶段潮流方法（Two-stage power flow method for calculation voltage stability margin）[J].电力自动化设备（Electric Power Automation Equipment），2013，33（3）：72-77．

[4]杨小煜，周孝信（Yang Xiaoyu，Zhou Xiaoxin）.基于极小扩张系统方法的静态电压稳定临界点计算（Calcula⁃tion of the critical points of static voltage stability with minimally extended system method）[J].中国电机工程学报（Proceedings of the CSEE），2009，29（25）：32-36．

[5]Irisarri G D，Wang X，Tong J，et al.Maximum loadability of power systems using interior point non-linear optimiza⁃tion method[J].IEEE Trans on Power Systems，1997，12（1）：162-172．

[6]赵晋泉，张伯明（Zhao Jinquan，Zhang Boming）.改进连续潮流计算鲁棒性的策略研究（A study on the strategy for improving robustness of continuation power flow com⁃putation）[J].中国电机工程学报（Proceedings of the CSEE），2005，25（22）：7-11．

[7]汤涌.电力系统电压稳定性分析[M].北京：科学出版社，2011．

[8]张伯明，陈寿孙，严正.高等电力网络分析[M].北京：清华大学出版社，2007．

[9]Ajjarapu V，Christy C.The continuation power flow：A tool for steady state voltage stability analysis[J].IEEE Trans on Power Systems，1992，7（1）：416-423．

[10]张尧，张建设，袁世强（Zhang Yao，Zhang Jianshe，Yuan Shiqiang）.求取静态电压稳定极限的改进连续潮流法（Improved continuation power flow algorithm for obtain⁃ing the limit of static voltage stability）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），2005，17（2）：21-25.

[11]王成山，魏炜（Wang Chengshan，Wei Wei）.一种改进的步长控制连续性潮流计算方法（An improved continua⁃tion method with controlled step size）[J].电工技术学报（Transactions of China Electrotechnical Society），2004，19（2）：59-63．

[12]Zhao Jinquan，Chiang H D.Enhanced look-ahead load margin estimation for voltage security assessment[J].IEEE Trans on Power Systems，2003，6（3）：2640-2645．

[13]Zarate L A L，Castro C A.Fast computation of security margins to voltage collapse based on sensitivity analysis [J].IEE Proceedings-Generation，Transmission and Distri⁃bution，2006，153（1）：35-43.

Fast Calculation Method for Saddle Node Bifurcation Point of Voltage Stability of Power System

ZHU Yongqiang1，2，LIU Guangye1，CAO Bin3，KUANG Yilei1，KANG Zhihao1，LIU Yuanyuan1
（1.College of Electrical and Information Engineering，Hunan University，Changsha 410082，China；2.Yueyang Power Supply Company，State Grid Hunan Electric Power Company，Yueyang 414000，China；3.Maintenance Company，State Grid Hunan Electric Power Company，Changsha 410000，China）

TM712

A

1003-8930（2017）09-0086-07

10.3969/j.issn.1003-8930.2017.09.015

2015-09-23；

2017-01-11

国家自然科学基金资助项目（51577053）

朱永强（1989—），男，硕士，助理工程师，研究方向为电力系统电压稳定分析与控制。Email：1090321338@qq.com

刘光晔（1960—），男，博士，教授，博士生导师，研究方向为电力系统分析与控制、输变电技术、铁路牵引供电系统、电力系统继电保护。Email：liuguangye@21cn.com

曹 斌（1991—），男，硕士研究生，研究方向为高电压与绝缘技术。Email：791589645@qq.com