电力系统静态电压稳定临界点的快速计算方法

黎 晓， 刘崇茹， 蔡 晖， 辛蜀骏， 郑 乐

(1.新能源电力系统国家重点实验室(华北电力大学)， 北京 102206；2.国网江苏省电力有限公司经济技术研究院， 江苏 南京 210008；3.国网经济技术研究院有限公司， 北京 102209)

0 引 言

20世纪70年代以来，国内外发生过多次大的电压崩溃事故，学者们对电压稳定的评估、预警与预防控制开展了大量的研究工作[1～3]。电压稳定临界点的快速、准确计算对于电压稳定评估以及电压稳定控制具有重要意义。分岔点对应的负荷大小可表示当前运行点的负荷裕度，可用于电压稳定的预警与预防控制。临界点的灵敏度信息，可以为配置无功补偿[4]，切负荷[5]等措施提供重要信息。

常用的临界点计算方法有连续潮流法和分岔点直接计算法。连续潮流法[6]具有鲁棒性强、能同时考虑发电机的无功功率限制等约束条件[7]的优势，很容易处理极限诱导分岔。但是准确地计算鞍结型分岔点，现有的连续潮流计算方法计算量大，难以满足在线计算的需求。主要表现在以下两个方面：(1)由于PV曲线在鞍结型分岔点附近的非线性特性明显，预测步长过大会导致预测点距离PV曲线较远，使得校正环节中牛顿迭代法的迭代次数显著增加，甚至导致牛顿法不收敛[8]或扩展潮流方程无解[7]。通过步长调整策略和参数化过程能够明显改善连续潮流的收敛性，但会使得连续潮流的步长在远离鞍结型分岔点时使用较大的步长，在靠近分岔点时使用较小的步长[9，10]，这种策略需要较多的计算量来逼近负荷极限点。(2)“斜锐角型”PV曲线，可能会使得校正环节中的牛顿法不收敛[11]。这种不收敛现象在全局性参数化方法与不合适的局部参数化方法中都可能出现。出现这种情况时，计算过程中需要调整参数化方法[9]，重新计算校正步，或者采用线性搜索等求解方式[12]，会显著增加总的迭代计算次数，降低计算速度。

另外一类方法基于分岔点的判别条件，直接计算分岔点[13～17]，被称为直接法。直接法通常采用牛顿迭代求解，在合适的初值选择下能够快速准确地计算出鞍结型分岔点。但是牛顿法对初值较为敏感，不合适的初值可能导致计算发散。文献[18-20]提出的直接计算方法，不需要设置初值，提高了直接迭代方法的收敛性，但是计算节点的选择直接影响到该方法的收敛特性。这类直接法对于PV曲线呈现“抛物线”形状的薄弱节点具有较好的收敛特性，但是当初始潮流解远离负荷极限点时，薄弱节点的选择并不容易，因为PV曲线的特性在逼近临界点处才显现出来[21，22]。文献[23,24]尝试将连续潮流法和直接法的优势综合起来进行计算，通过连续潮流快速穿过鞍结型分岔点，寻找鞍结型分岔点的近似解作为直接法的迭代初值，然后用直接法计算准确的鞍结型分岔点。但是，由于前述分析可知，连续潮流方法在PV曲线上穿过鞍结型分岔点时，为了解决收敛性的问题而采用的步长调整策略或参数化方法均需要较多的计算量，限制了该方法的计算效率。

本文结合连续潮流和直接法两者的优势，提出了一种静态电压稳定临界点的快速计算方法。该方法不需要连续潮流计算穿过分岔点，避免了连续潮流为解决收敛问题而增加的计算量，同时又能够为直接法找到接近负荷极限点的合理初值，从而保证直接法的收敛性。首先，通过大步长的连续潮流快速靠近临界点，当判断临界点为极限诱导分岔时，则结束临界点计算；当判断临界点是鞍结型分岔点时，以最后一步连续潮流的解为初值，通过直接法快速准确地计算出鞍结型分岔点。

本文的创新点主要体现在：(1)提出了一种线性程度较好的鞍结型分岔点逼近指标，用以评估系统逼近鞍结型分岔点的程度，并根据这一指标确定当前连续潮流计算结果能够满足直接法初值的要求；(2)提出了连续潮流与直接法的切换策略，显著降低了连续潮流的计算时间，同时保证直接法的初始值接近负荷临界点；(3)提出了一种海森矩阵的隐式计算方法，提高了直接法的计算效率。

1 鞍结型分岔点逼近指标

1.1 指标的理论基础

鞍结型分岔点指标在第i个负荷节点处的体现记为ci，其在数值上应具备两个特征：

(1)在鞍结点分岔点处ci等于0。保证在鞍结点分岔点处ci=0能够保证在计算过程中通过ci接近于0的程度而快速判断是否已接近鞍结型分岔点。

(2)指标ci数值随着参数λ增大应呈现近似线性的变化趋势。由于某一个参数λ所对应的指标ci并不具有预测的能力，此指标的关键特征在于其数值随着参数λ增大的线性变化趋势。

定义：鞍结型分岔点逼近指标在第i个负荷节点的计算值ci：

(1)

式中：dVi表示节点i的电压幅值的微分；dλ/dVi可以通过连续潮流计算过程直接获取；Ii表示负荷节点i的负荷电流。

指标ci还可以写成式 (2)的形式

(2)

式中：w和τj为方程 (3)的解，w表示PV曲线上切向量的方向；w(i)表示向量w的第i个元素。fx表示潮流方程的雅各比矩阵，fλ表示潮流方程对参数λ的偏导数，ejT表示除了第j个元素外其余元素全为0的行向量。τj的上标j和ej的下标j相对应，与ci的下标i无关，所以可以任取j算出所有负荷节点的ci值。n等于除平衡节点之外的节点数量。式 (3)中的τj与Abbott给出的鞍结型分岔点测试函数相同，因此此方法对应相同的理论基础[25]。已经证明[26]，在鞍结型分岔点处，τj/(w)i=0，因此对应ci=0。

(3)

1.2 指标的物理意义

由于视在功率与电压电流幅值的关系S=VI，其导数可以表示为式 (4)。

dS=IdV+VdI

(4)

因此

(5)

式中：dSi/dλ表示节点i的视在功率的导数，也是负荷功率的增长方向在节点i上的投影，可以视为已知常量，那么容易使得其满足dSi/dλ0。ZLi=Vi/Ii，表示负荷等值阻抗模。ZTHi=-dVi/dIi，在其他文献中被称为动态等值阻抗模[27]，与通常所采用的戴维南等值阻抗模不同。通常所采用的戴维南等值阻抗采用电压和电流相量变化量的比值，这里仅仅采用幅值变化量的比值。在物理意义上，式 (5)的形式类似于阻抗模值匹配指标，可以同样解释为在鞍结型分岔处，负荷等值阻抗模等于动态等值阻抗模。

1.3 保持指标具有连续特征的修正方法

与其他鞍结型分岔点逼近指标一样，在发电机节点类型发生变化时，指标ci可能会发生突变。

因此需要对式 (1)的定义式进行修正，使得该指标能够抵消潮流方程中节点类型变化带来的影响。该指标在突变点前后具有相似的变化斜率，所以只要消除突变量，即使在节点电压类型发生变化的情况下，保证指标随参数λ变化的光滑程度仍然较好。

(6)

图1 鞍结型分岔点逼近指标突变处理示意图Fig. 1 Schematic diagram of sudden change treatment of saddle-node bifurcation point approximation index

1.4 指标的线性化方法

(7)

其中参数βi可以通过最小方差的方式辨识，辨识方法如下。

步骤一：

A=[λ-E(λ),-c+E(c),-λ·c+E(λ·c)]

其中，λ=[λ1,λ2, …,λl]T；c=[ci1,ci2,…,cil]T；λ·c=[λ1·ci1,λ2·ci2,…,λk·cil]T；E(*)表示求平均数。

步骤二：

式中：Σi表示样本的方差；Aj,:表示A矩阵的第j行。

步骤三：求解式 (8)中的z向量。

(8)

其中第二个约束是为了使得βi与ci的符号异号，避免式 (7)中分母等于0。z(1)与z(3)分别表示z向量的第1个元素和第3个元素。步骤三是一个带约束的二次规划问题，可以利用现有的工具很容易求解。

步骤四：βi=z(1)

2 连续潮流与直接法的切换决策

2.1 节点类型转化识别的补逻辑

将发电机无功功率方程写成QGi(λ)-fqi(x)=0，其中QGi(λ)表示发电机i发出的无功功率，fqi(x)表示通过节点电压计算发电机节点i的无功功率的表达式。

(9)

式中：xj表示潮流方程中的变量，包括PQ节点的电压幅值和相角、PV节点的电压幅值。dxj/dλ由连续潮流方法中的预测环节得到。根据发电机无功的上限值和下限值，计算每个发电机节点的λ增量的估计值，如式 (10)所示。

(10)

(11)

式中：GPV表示发电机节点中PV节点的集合；GPQ表示发电机节点中PQ节点的集合；dVGi/dλ由连续潮流方法中的预测环节得到。然后取λ增量大于0中的最小值作为节点类型发生转化的估计点，如式 (12)和式 (13)所示。

(12)

(13)

但是当两个节点几乎同时发生类型转换时，提前预判节点类型转化难以识别其先后顺序，可能收敛到不满足发电机无功限制的解。为了避免这种情况下的多次缩小步长再次计算校正环节，本文增加以下两个补逻辑。

(1)当预判节点i将发生节点类型转化，通过计算收敛之后，节点j出现了越限，那么将节点j的类型转换之后，重新计算节点i的节点类型转化临界点。

(2)如果预判节点i和节点j发生节点类型转换时刻的λ值之差小于，且按照补逻辑1转换节点i的类型后仍然有新的发电机节点出现越界，那么就恢复节点i的节点类型，再以节点j作为类型转换的关键节点计算其节点类型转化的临界值。其中可以取一个较小的数，例如=0.001。

2.2 失稳类型的识别

失稳类型的正确识别依赖于节点类型的识别和鞍结型分岔点逼近指标的失稳判断。所以在计算流程上先识别节点类型转化，然后判断失稳类型。

每次计算得到节点类型转化临界点之后，通过节点类型转化前和转化后的切向量的符号变化可以判断该点是否为极限诱导分岔点[29，30]。如果是极限诱导分岔，则不需要进入鞍结型分岔计算的直接法，直接输出结果。

当鞍结型分岔点逼近指标判断鞍结点时，就不用再进入连续潮流的校正环节，转为采用鞍结型分岔点的直接计算方法，能够快速收敛到分岔点。判断的方法如下：

(1)在每一步连续潮流法的校正环节求解前，综合突变处理前的指标ci与母线电压幅值选择薄弱母线k。比较所有负荷节点的鞍结型分岔点逼近指标ci值，将其绝对值从小到大排序，记节点i的序号为Ord(ci)。相似地，将负荷电压母线的电压幅值也从小到大排序，记节点i地序号为Ord(Vi)，然后计算Ord(ci)+ Ord(Vi)，取结果最小的母线编号为k。

(14)

(15)

整体的具体决策过程见图2。为了提高算法的鲁棒性，保留了回到预测环节缩短步长再次计算校正环节的处理方式。

图2 总体计算流程Fig. 2 Overall calculation process

3 鞍结型分岔点计算的改进直接法

3.1 直流法的计算步骤

静态电压稳定的鞍结型分岔点计算等同于求解以下扩展潮流方程中的(x,λ)。

(16)

式中：rank(*)表示矩阵的秩；fx表示潮流方程的雅各比矩阵；g(x,λ)=0是鞍结型分岔点的判据，函数g也被称为测试函数[31]。

方程(16)的解即为鞍结型分岔点。利用牛顿法求解时需要形成方程(16)的增广雅各比矩阵J～，如式 (17)所示。fλ表示潮流方程对参数λ的偏导数。当f(x,λ)=f1(x)+f2(λ)时，fλ与x无关。gλ=0，gx表示g函数对代数变量x的偏导数。

(17)

函数g的表达式具有不同的等价形式。本文采用式 (18)形式的函数g，与式 (16)中的函数g等价[31，32]。

g=(d-cTA-1b)

(18)

式中：矩阵A表示fx矩阵中去掉第k行第k列后剩下的部分。向量b表示fx中第k列中除k行之外的元素。cT表示fx中的第k行中除第k列之外的元素。向量d表示fx矩阵中第k行第k列的元素。k是一个需要设定的参数，对直接法的收敛性有直接影响，本文选择k值的方法与鞍结型逼近指标预测时选择薄弱母线k的方法相同。

偏导数gx需要求解式 (18)中矩阵A的导数，也是f(x)在x点处展开的海森矩阵的一部分[33]。由于海森矩阵属于三维矩阵，其维度较高，不利于算法快速求解。

3.2 海森矩阵的快速计算方法

3.2.1 引理

引理1：一个矩阵，如果fx，R，T都是连续函数或者常数，那么一定存在唯一的向量u,W, 和实数g使得式 (19)和式 (20)成立，且fx矩阵的秩等于2n-r-1的充分必要条件是g=0[34]。

(19)

(20)

本文在引理1的基础上表示gx，并推导海森矩阵的隐式计算方法。将式(19)等号两边左乘(uT,g)，与式(20)联立得到：

g=uTfxW

(21)

取R=[01×n,1]T，TT=[01×n,1]，uT=(-cTA-1,1)P，W=P-1(-(-A-1b)T,1)T。其中矩阵P为一个置换矩阵，用于将雅各比矩阵fx中的第k行k列移动到最后一行一列。置换矩阵与矩阵fx相乘可以得到如式 (22)所示的分块形式，其中A∈R(2 n-r-1)×(2 n-r-1),b∈R(2 n-r-1)×1,c∈R(2 n-r-1)×1,d∈R。r表示PV节点(给定注入有功功率与电压幅值的节点)的数量。那么，将R，T，u和W的表达式代入式(21)中就可以得到式 (18)。说明在以上条件下，式 (18)中的函数g适用于引理1。

(22)

进一步通过式(19)和式(20)求解g对潮流方程变量x的导数的表达式。记g对任意变量的导数为g′，对式(21)等号两边同时求导得到式(23)。

g′=(uT)′(fxW)+uT(fx)′W+uTfxW′

(23)

根据式(19)和式(20)中fxW=Rg和uTfx=gTT，式(23)可以改写为以下形式。

g′=(uT)′(Rg)+uT(fx)′W+gTTW′

(24)

根据式(19)和式(20)中，TTW=1和uTR=1，得到

TTW′=-(T′)TW

(25)

(uT)′R=-uTR′

(26)

由于向量T向量R为常数向量，所以式(25)和(26)中的T′和R′都等于全0的向量。因此，将式(25)和式 GOTO(26)代入式(24)，得到g对任意参数的导数，如式(27)所示。

g′=uT(fx)′W

(27)

将式(27)应用到对所有代数变量求导，得到gx的表达式。

gx=uTfxxW

(28)

3.2.2 海森矩阵的隐式计算方法

式(28)中的fxx是潮流方程的海森矩阵，在大规模电力系统中，海森矩阵维度非常高。尽管对于电力系统的潮流方程，其海森矩阵高度稀疏，但稀疏矩阵计算技术很难在三维矩阵上直接利用。

根据方向导数的定义，沿着向量W的方向导数等于函数梯度与方向向量W的内积，即

(29)

式中：▽表示梯度算子，有▽f=fx；下标(i,j)表示矩阵中第i行第j列的元素；单个下标i或者j表示向量中第i个或者第j个元素。

同理可得

(30)

式中：(▽f)x=fxx，下标(:,j)表示矩阵的第j列。

那么，

(31)

以上分析过程说明，只需要把向量W当作x的一个增量，在这个增量方向上两次形成雅各比矩阵矩阵，就可以很方便的计算出gx。该计算过程计算量小，很容易在现有的潮流计算程序上实现。ε可以取一个较小的数值，理论上ε取得越小，gx的计算结果越精确，但是还需要考虑数值计算过程中的舍入误差，所以不宜取得过小。对于双精度浮点数，可以设置ε=10-7。

本文提出的海森矩阵隐式计算方法，显著提高了直接法的计算效率。

4 算例分析

在配置为i7-8700 处理器、8 G 内存的计算机上，用软件 Matlab 编制本文提出的计算程序，以算例case2736sp, case2869pagase, case6468rte[35], case9241pegase[36]为例进行分析，其网络规模分别为2 736条母线，2 869条母线，6 468条母线，和9 241条母线。本节选择了两种负荷功率增长方式，增长模式1表示全系统负荷功率增长的情况，假设全系统的负荷有功功率、无功功率等比例增长，全系统的发电机(平衡节点除外)发出的有功功率也等比例增长以分担负荷功率。增长模式2表示局部区域负荷功率增长的情况，假设编号最小的10个负荷节点的有功功率与无功功率等比例增长，全系统的发电机(平衡节点除外)等比例增长有功功率以分担负荷功率。大步长连续潮流采用拟弧长的参数化方法，直接法中采用隐式海森矩阵计算方法。

4.1 鞍结型分岔点逼近指标的性质

尽管节点类型转化会造成鞍结型分岔点逼近指标的突变，图3表明所提出的式 (6)的处理方式可以得到整体光滑性较好的指标变化曲线，鞍结型分岔点逼近指标的突变点的不连续特性得到了很好的修正，修正后的连续性是该指标在非线性变换后呈现线性特征的基础。

图3 鞍结型分岔点逼近指标(非线性变换前，突变修正前后)模式2Fig. 3 Saddle-node bifurcation point approximation index (before nonlinear transformation, whether jump corrected or not)

基于修正后的结果，通过参数β的辨识和式 (7)的非线性变换，得到的鞍结型分岔点逼近指标就具有了比较理想的线性趋势，如图4所示。

图4 鞍结型分岔点逼近指标(非线性变换后)Fig. 4 Saddle-node bifurcation point approximation index (after nonlinear transformation)

以上结果验证了所提出的鞍结型分岔点逼近指标的良好特性：在经过节点类型转化造成的突变修正和非线性变换之后，具有较好的线性特征。

4.2 补逻辑的有效性

考虑发电机节点类型转化，拟弧长的预测步长σ=1.0，计算不同算例、不同负荷增长模式下的静态电压失稳分岔点。计算过程中所触发的补逻辑次数如表 1所示。

表1 补逻辑触发次数比较Tab.1 Comparison of trigger times of complementary logic

表 1中SNB表示鞍结型分岔，LIB表示极限诱导分岔。结果表明，在绝大多数情况下，补逻辑1就足够有效地避免节点类型转化造成的缩小步长和再次计算校正环节。仅在case2736sp系统的负荷增长模式2的计算过程中触发了补逻辑2。

在算例case6468rte中，母线6 372与母线1 366几乎同时发生节点类型转化。虽然在预测步长比较大的情况下(即σ=1.0)通过式 (10)和式 (11)错误地识别到母线1 366先发生节点类型转化，从而造成母线6 372的无功越界。但是通过补逻辑1，将母线6 372转化为PQ节点后，计算结果正确的收敛到了预测步长σ取0.1时的正确解。

4.3 计算速度比较

为了比较步长大小对于鞍结点附近牛顿法计算收敛性以及计算速度的影响，不考虑发电机无功功率限制，比较以上4个算例在不同步长下的计算时间，如图5所示。

图5 计算时间比较Fig. 5 Comparison of computation time

图5中左侧线表示本文所采用的修正方法：通过鞍结型分岔点逼近指标评估运行点是否接近运行极限，然后通过直接法计算分岔点。右侧线表示采用拟弧长参数化方法的连续潮流计算方法。当连续潮流法计算过程中校正环节的牛顿法计算发散时，将步长缩小为原来的一半，再次计算校正环节[7]。

在连续潮流法的计算中，如果逼近临界点且步长过大，例如图5 中步长等于20的情况，因为牛顿法计算不收敛或者收敛慢，造成更多的校正环节重复计算以及更多的牛顿法迭代次数，反而降低了计算速度。

图5中的计算结果表明，在相同预测步长的情况下，本文方法比连续潮流法计算速度快，这说明通过本文提出的鞍结型分岔点逼近指标预测逼近鞍结型分岔的程度，判断当系统运行点接近负荷增长极限时，转为直接法计算鞍结型分岔具有计算速度上的优势。

4.4 计算精度比较

由于计算结果都是收敛后的潮流值，计算误差都在10-10以下，因此不论是本文方法还是连续潮流计算方法，所得到的计算结果都在PV曲线上。

图6比较了不同步长下的分岔点计算结果。可以看出：(1)本文的方法通过直接法计算分岔点，计算精度不受连续潮流步长的影响。而采用连续潮流法计算分岔点的结果会受到预测步长影响；(2)由于所得到的计算结果都在PV曲线上，且本文方法的计算结果都大于等于传统连续潮流的计算结果，说明本文的计算结果更加接近鞍结型分岔点。

图6 临界点计算精度的比较Fig. 6 Comparison of calculation accuracy of critical point

4.5 海森矩阵隐式计算效率

目前的稀疏矩阵技术针对二维矩阵之间或二维矩阵与向量之间的运算进行了优化，且本文的隐式形成海森矩阵能够比较方便地利用现有的稀疏矩阵技术。但是显式形成的海森矩阵属于三维矩阵，不方便直接利用现有稀疏矩阵技术，所以两者之间的计算时间难以进行有效的比较。因此，本文通过比较计算的非零元素个数与乘法次数来进行说明。

显式计算式 (28)中的海森矩阵fxxW，乘法计算的次数约等于fxx中非零元素的个数，因为向量W是稠密的；隐式的海森矩阵形成方法需要形成一次雅各比矩阵的计算量。根据fxx矩阵的元素组成，通过估计非零元素个数，从而估计形成过程中的乘法次数，如表 2所示。表格中的单位1等于导纳矩阵的非0元素个数，分析结果表明，在各个算例下，隐式形式海森矩阵的计算量都小于显式形成海森矩阵的计算量，而且隐式方法更容易在现有的潮流计算程序上实现。

表2 隐式形成海森矩阵与显式形成海森矩阵的计算量比较

5 结 论

本文提出了一种静态电压稳定临界点的快速计算方法。该方法具有以下3个特点：

(1)提出了一种线性性质较好的鞍结型分岔点逼近指标，利用该指标随负荷参数λ变化的线性变化优势，实现了对潮流解逼近鞍结型分岔点的有效估计。

(2)提出了连续潮流法与鞍结型分岔点计算直接法的切换策略，通过鞍结型分岔点逼近指标的判断，在逼近鞍结型分岔点时从连续潮流法切换到直接法，避免了连续潮流法难以快速准确地计算鞍结型分岔点的问题，同时为直接法提供了较好的计算初值，实现了两种计算方法的高效结合。

(3)提出了一种海森矩阵的隐式计算方法，具有容易实现、计算快的优点，提高了此类直接法的计算效率。在直角坐标系下，潮流方程的海森矩阵为常数矩阵，在本文的基础上可以进一步加速临界点的计算。

由于所提出的方法可以方便地处理节点类型的转化，在本文的基础上还可以考虑无功补偿设备的其他约束，以及增加新能源机组和直流输电系统的静态模型，这是本文下一步的工作。