电压稳定极限诱导分岔点快速追踪方法

钟 浩，姚 丹

（三峡大学 梯级水电站运行与控制湖北省重点实验室，湖北 宜昌 443002）

0 引言

电力系统在当前运行状态下的电压稳定程度[1-2]是运行管理人员最为关心的问题之一。负荷裕度指标[3]在体现电压稳定程度上具有直观性强和线性度好等优点，因此得到了广泛的应用。负荷裕度指标是指在功率注入空间中，从当前运行点出发，按给定方向增长负荷直至电压崩溃时，当前运行点与电压分岔点之间的距离。因此，要获得负荷裕度指标关键是要计算得到系统的电压稳定分岔点[4-6]。

若只考虑电力系统的静态特性，电压静态分岔类型[7]主要有鞍结分岔 SNB（Saddle Node Bifurcation）和极限诱导分岔LIB（Limit-Induced Bifurcation）。目前大多数静态电压分岔点的计算都是针对SNB点。但在某些情况下，系统可能在SNB点之前就发生了LIB，导致电压崩溃。只计算系统的SNB点可能会使指标过于乐观，给系统电压稳定带来隐患。因此，研究LIB的电压稳定裕度计算方法有重要意义[8-9]。

随着负荷参数或其他参数的逐渐增大，系统的无功源发生无功越限，可能会伴随LIB发生[10-12]。现阶段求取诱导分岔点的方法主要有连续潮流方法、直接法和优化方法。连续潮流方法[13]可追踪到LIB点，但是步长难以控制且追踪时间较长；直接法[14]的应用不如连续潮流方法普遍，因为直接法求解的非线性方程组维数高、占用内存空间大且不易利用稀疏技术；优化方法[15]易于加入各种约束，但随着约束条件的增多，模型的收敛性难以保证，计算时间越来越长。

针对LIB点追踪的计算精度和计算速度问题，本文提出了一种LIB点的快速追踪方法。首先给出了无功越限导致系统电压失稳的判据，在远离分岔点时利用局部曲线拟合技术，自适应确定负荷增长步长，快速逼近分岔点，并将无功越限的节点由PV节点转换成PQ节点；在靠近分岔点处采用二阶灵敏度方法对增长步长内PV节点到PQ节点转换的PV节点集的触发无功上限排序，校正增长步长，再结合步长折半搜索方法和分岔点判据易追踪到系统的LIB点。

1 负荷裕度求解模型

负荷的有功和无功以及发电机PV节点的有功沿某一方向连续增长，则表征电力系统运行情况的代数方程可表示为：

其中，x为状态变量；ρ为发电机功率分配因子；λ为负荷参数。

可将式（1）详细表示为：

其中，PiG0、QiG0分别为节点i的初始发电机有功、无功功率；PiL0、QiL0分别为初始负荷有功、无功功率；Pi、Qi分别为所计算出来的节点i有功、无功功率；ρiG为发电机i的分配因子；nL为参与增长的负荷个数。

负荷参数λ的变动范围为：

λ=1对应基态负荷，λ=λcir对应最大负荷。本文假设负荷按原功率因数增长，同时也考虑PV节点的无功极限问题。若某一PV节点的无功超出其限值，则将其设置成PQ节点；若在减小步长后其计算无功在极限值之内，则将其重新置成PV节点。在这种PV节点与PQ节点转换的过程中，可能会伴随LIB发生。

2 LIB点的识别

发电机无功越限时，发电机的运行状态将发生变化，节点类型由PV转换成PQ，此时会有2种情况出现，如图1和图2所示。图1中发电机无功越限减小了电压稳定裕度但系统未发生电压失稳，运行点位于节点转换后λ-U曲线上半支，此时dUj/dλ＜0（j为电压弱节点号）；图2中发电机无功越限后的LIB导致了电压崩溃，运行点位于节点转换后λ-U曲线下半支，此时dUj/dλ＞0。本文重点分析图2情况。

图1 极限诱导转换过程Fig.1 Limit-induced transition

图2 极限诱导分岔Fig.2 Limit-induced bifurcation

在图2中PV节点集中有转换节点a和b，预测点i+1已经在可行域之外，使得潮流不收敛，因此需通过式（15）对PV节点集中的触发无功上限排序，将运行点拉回转换点a。在转换节点a，发电机a达到无功极限，即QGa=QGliam，Ua=Ua0。随着负荷参数增大，其电压也将增大，此时运行点位于λ-U曲线（QGa=QlGima）下半支，诱发了系统的突然电压崩溃，分岔类型为LIB，满足QGa=QGlima、dUj/dλ＞0 判据。 若分岔类型为SNB，则可在逼近点斜率绝对值大于某一值情况下基于局部曲线拟合技术求取SNB点[4]。

将式（1）写成潮流修正方程的形式，表示为：

式（4）对负荷参数求导：

在稳定运行点处雅可比矩阵J为已知，因此很容易得到dU/dλ。可快速判断运行点是在λ-U曲线的上半支还是下半支。

3 LIB点追踪方法

3.1 局部曲线拟合确定增长步长

节点电压与负荷参数呈非线性，为了确保增长步长预测的相对准确性，有必要计算节点电压对负荷参数的2阶导数，通过局部曲线拟合获得增长步长。继续对式（5）求导，得到：

获得当前运行点的弱节点电压对负荷参数二阶灵敏度信息之后，可用局部曲线拟合技术求取当前运行点到预测点的负荷增长步长。局部曲线拟合式为：

其中，α、β为常系数；Δλ和ΔUc分别为当前运行点到预测点处的负荷增长步长和电压弱节点的电压幅值变化量。只要计算出α、β和ΔUc即可得到负荷增长步长Δλ。从式（7）可知，越靠近分岔点确定的增长步长越小，具有自适应性。

式（7）对负荷增长步长Δλ求1阶导：

继续对式（8）求2阶导：

由于运行点 i的 dUc/dλ、d2Uc/dλ2已知，联立式（8）、（9）可以解出参数 α、β。

然后式（7）对节点电压变化量求导，可得：

3.2 PV节点集触发无功上限排序

LIB往往由无功越限所引起，因此在增长步长内须计算出发电机输出无功，以确定发电机无功是否越限。如有越限，将其纳入PV节点集，并对节点集触发无功上限排序，识别LIB点。

发电机PV节点发出的无功功率可表示为：

由式（11）可得到在运行点i处发电机PV节点的无功出力QPV对负荷参数λ的各阶导数。其中1阶导数为：

由式（12）再次对负荷参数λ求导，可得到QPV对负荷系数λ的2阶导数：

根据式（14）很容易确定增长步长内的PV节点集：

LIB点最终是通过某台发电机无功极限约束来确定的，为了追踪LIB点，需要对PV节点集触发无功上限顺序进行排序，然后取最先触发的节点对应的增长步长Δλ，其求解式如式（15）所示：

其中，ΩPV为当前运行点与预测点之间PV节点到PQ节点转换的PV节点集；nG为PV节点到PQ节点转换的节点个数；QGlimj为发电机j无功出力上限；对应中相应值；由于在当前运行点已经求出发电机PV节点的无功出力QGj，所以根据发电机无功出力对负荷参数导数很容易确定节点集中触发无功上限的顺序。由于Q与λ的非线性，采用二阶灵敏度法预测的增长步长可能并不准确，后面将通过负荷增长步长折半搜索方法和设定阈值来重新识别。

系统的无功约束转换点数目众多，在远离分岔点时，无功约束转换点不会使系统发生电压崩溃，因此只需将无功越限的PV节点转换成PQ节点即可。当离分岔点较近（潮流不收敛情况），才需对预测步长内PV节点集中节点进行识别。因此，避免了对所有转换点识别所带来的巨大计算量。

4 算法步骤

下面给出系统电压稳定LIB点追踪基本步骤。

a.初始化。初始负荷系数λ0=1.0 p.u.，设置斜率常数 r、斜率比较值 k、最大限值 Δλmax、阈值 ε，运行点斜率 S=0（S=dUj/dλ），标志位 F=0。

b.潮流计算。如果潮流收敛，则令F=0，计算运行点斜率S，跳到步骤c；如果潮流不收敛则跳到步骤e。

c.判断斜率S正负：如果S为正，则为LIB点，结束程序；如果S为负且绝对值小于某一设定值，则继续向下执行；如果S为负且绝对值大于某一设定值，则根据曲线拟合技术计算SNB点值，程序结束。

e.F=0，则取PV节点集中最先无功越界发电机的无功/电压约束转换点作为预测点，并根据式（15）校正负荷增长步长Δλi，其他转换节点置回PV节点，令 F=1，回到步骤 b；F=1 且 Δλi＞ε，则将上步所取的最先转换点还原为PV节点且负荷增长步长折半，即 Δλi=Δλi/2，回到步骤 b；F=1 且 Δλi＜ε，则该点为LIB点，程序结束。

靠近分岔点处，PV节点集中节点的转换使系统结构发生突变，可能步长过大会使潮流不收敛，此时用折半搜索方法确保数据具有一定连续性，使本算法具有很好的稳定性。另外，在整个过程中，求取λ-U曲线上运行点时具有一定连续性，对系统计算数据随时存储和随时读取，使得潮流数据具有一定的连续性。因此，算法特点决定了本方案是可行的。

5 算例

采用IEEE 118节点系统作为仿真算例验证本文所提方法的快速性和有效性。表1给出了斜率常数r为1.15和1.20时追踪LIB点的计算过程，表中n为潮流计算步数，粗体数字即为LIB点。

表1 不同参数逼近过程比较Table 1 Comparison of approaching process between different parameters

当r=1.15、n=8时，预测步长使得潮流不收敛，斜率标记为“—”；在n=9时，通过式（15）取PV节点集中最先转换的节点对应的增长步长，潮流不收敛；在n=10时，通过折半搜索方法将转换点置回PV节点重新识别；当n=11时，将转换点置成PQ节点，从而S=0.3852为正，说明运行点已位于λ-U曲线下半支，转换点触发了LIB，整个追踪过程如图3所示。当r=1.20、n=6时，同样运行点超出可行域，潮流不收敛，此时根据PV节点集，通过式（15）确定n=7时的负荷增长步长，将系统运行点拉回可行域；当n=9时，因S=0.0292为正，说明运行点已位于λ-U曲线下半支，转换点触发了LIB，整个追踪过程见图4。

表2为不同方法的计算值比较，其中NP和NCPF分别为本文所提方法和连续潮流法的潮流计算次数，λP和λCPF分别为本文所提方法和连续潮流法的最大负荷值。本文所提方法计算IEEE118节点系统，在选取不同斜率常数r时与连续潮流法所得最大负荷值[16]的误差都在2%以内，但本文方法所用潮流计算次数NP分别为11次和8次，而连续潮流法计算次数为18次。从算例可见，本文所提方法在计算精度相当的情况下效率要高出连续潮流法1倍左右。

从表2还可以看出，斜率常数r的取值对潮流计算次数有一定的影响。常数r的取值大体思路是：r取值须大于1；当前运行点斜率较小且离分岔点较远（可根据算法中增长步长后潮流计算是否收敛来确定当前运行点离分岔点远近）的时候选择较大的r值；当前运行点斜率较大且离分岔点较近的时候选择较小的r值；特别是接近分岔点处，r值越小则计算越精确；在后续研究中需考虑选择最佳斜率常数r的方法。

图3 r为1.15时逼近LIB点示意图Fig.3 Schematic diagram of LIB point approaching，when r=1.15

图4 r为1.20时逼近LIB点示意图Fig.4 Schematic diagram of LIB point approaching，when r=1.20

表2 不同方法计算值比较Table 2 Comparison of calculative results between different methods

6 结论

本文通过二阶灵敏度确定负荷增长步长和增长步长内PV节点到PQ节点转换的PV节点集。由于采用步长限制策略且步长预测具有自适应能力，因此二阶灵敏度具有较高的准确性。二阶灵敏度只需在潮流计算基础上增加非常少的计算量即可获得。此外，本文只需在靠近分岔点处（潮流不收敛情况）才对预测步长内PV节点集中节点进行识别，避免了对所有转换点识别所带来的巨大计算量，因此具有非常高的计算效率。本文方法的潮流计算次数明显少于连续潮流方法，且计算精度较高，适用于大规模电力系统的在线求解。