基于微分求解模型的大规模电力系统病态潮流计算

林 珑 李建杰

(1. 重庆水利电力职业技术学院， 重庆 402160; 2. 重庆大学 输配电装备及系统安全与新技术国家重点实验室， 重庆 400030)

随着智能电网的进一步发展，我国电网规模不断扩大，结构日益复杂，区域间耦合日益紧密[1-3]．面对以上变化，各级调度应具备全网统一的网络分析能力，传统的潮流计算方法面临挑战[4]．而国民经济的迅速发展使得电网常常重负荷运行，此时电力系统潮流方程呈现病态特征，需要特殊的计算方法．

病态潮流计算相关研究由来已久．最早出现，也是迄今为止应用最广泛的方法是最优乘子法[5-9]，在牛顿-拉夫逊法的基础上进行鲁棒性变换，添加乘子并优选其数值，从而搜索病态潮流解．另一种病态潮流的求解算法是非线性规划法[10]，该方法将潮流计算中非线性方程组求解的问题转化为一组无约束的非线性规划问题，优势在于总是能够给出潮流有解或无解的解答，但计算量过大且结果依赖于优化求解方法的选择．此外，文献[11]构造了基于同伦法的病态潮流算法，该方法发挥同伦方法的大范围收敛性，使得病态潮流易于收敛．文献[12]则采用隐式Cholesky分解方法进行大规模病态潮流计算．文献[13]对多种病态潮流算法进行了论述以及比较分析．

以上方法能够完成特定应用场景下的病态潮流计算，并已被证明在一定规模的电力系统中是有效的．但其有效性受系统规模限制的特性明显，在系统节点数目超过1000后，方法的收敛性受到影响，同时计算耗时呈指数特性上升，对于大规模系统效率极低，无法满足全网一体化分析需求．因此有必要研究适用于大规模电力系统的新型病态潮流计算方法．

本文针对病态潮流方程的数学特性，采用连续牛顿法将非线性方程组转化为常微分方程组，并进行雅克比矩阵的相应处理，构建潮流计算的微分求解模型．在此基础上，采用一定的数值积分方法，能够准确高效地完成大规模病态电力系统的潮流求解．结合实际系统算例，其正确性得以验证．

1 数学原理简述

连续牛顿法是一种将非线性方程组求解问题转化为常微分方程组积分求解问题的数学方法，最初是通过将牛顿-拉夫逊法和显示欧拉法求解步骤进行类比和替换形成的[14]．现对其原理简述如下：对于如式(1)的非线性方程组

f(x)=0

(1)

牛顿-拉夫逊迭代法[15]在定点进行一阶泰勒展开，得到用方程雅克比矩阵表示的基本迭代格式如式(2)所示．式中x为待求解变量，J(xn)为在xn处的雅克比矩阵．

xn+1=xn+Δxn

Δxn=-J(xn)-1*f(xn)

(2)

另一方面，对于如式(3)的常微分方程组，

(3)

显式欧拉法[16]给出的积分公式为：

xn+1=xn+Δxn

Δxn=Δt*g(xn)

(4)

其中Δt为积分步长，一般有Δt=1．对比式(2)和式(4)可发现，如果做如式(5)的假设，则两式完全相同．

g(xn)=-J(xn)-1*f(xn)

(5)

由式(1)，(3)和(5)，可得到下式：

(6)

式(6)即为CNM的基本数学模型．结合一定的数值积分方法，该方程可以求解．设积分求解得到的结果为x*，易得：

g(x*)=-J(x*)-1*f(x*)=0

(7)

只要J(x*)非奇异，必有f(x*)=0，即x*为原非线性方程组的解．

应注意的是，式(6)是否具有稳定平衡点也需要验证，否则不能保证数值积分过程可以求解得到稳定结果．设g(x)=0的解为x=x*，则x*一定为原方程的一个平衡点，但该平衡点不一定稳定．由相关数学知识可知，x*为稳定平衡点的条件为g(x)在x=x*处的梯度矩阵的所有特征值均在左半平面(即右端函数雅克比矩阵负定)．对式(6)进行梯度矩阵的计算，

(8)

明显是负定矩阵，式(6)具有稳定平衡点．

2 病态潮流微分求解模型

2.1 针对潮流计算构建的微分方程

电力系统潮流方程表征电力系统中的功率平衡关系，由电力网络参数与节点电气参数构成．对于N节点电力系统，设其非平衡节点数目为n=N-1，其中PV节点数目为r，PQ节点数目为n-r，则在极坐标下其潮流方程形式如式(9)所示．其中ΔPi、ΔQi为节点i的有功和无功功率不平衡量，在系统稳态时其值应为0；PGi、PLi为节点i的发电机有功出力和负荷有功注入，QGi、QLi为节点i的发电机无功出力和负荷无功注入；Vi、Vj为节点i和j的电压幅值，θij=θi-θj为节点i和j的相角差值；Gij和Bij为网络节点电导矩阵和电纳矩阵中第i行j列的元素．

(9)

潮流方程本质为一组复杂的非线性方程组，可用式(1)表达，其中f为不平衡功率计算对应的映射，x为潮流计算中的未知状态变量，一般为节点电压幅值和相角．利用前述CNM数学方法，可将其转化为式(6)所示的常微分方程组进行求解，能够得到正确的潮流解．

当电力系统重负荷运行时，潮流方程呈现病态，其显著特征之一即为潮流雅克比矩阵接近奇异，对应于式(6)即为J(x)-1的计算效率很低，对于大规模系统该求逆过程更是极为困难，整体求解遇到瓶颈．为充分利用CNM方法对于基态和病态潮流的一致收敛性，避免由于雅克比矩阵求逆造成的收敛性和效率问题，本文对CNM方法进行了定雅克比矩阵处理，即在数值计算过程中采用固定的雅克比矩阵，其因子表或逆矩阵均只需要形成一次．采用定雅克比矩阵的CNM方法得到的微分方程如式(10)所示，其中J0为固定的雅克比矩阵，一般由初始点x0处计算得到．

(10)

采用定雅克比矩阵处理后，微分求解模型的正确性需要验证．首先，由于J0非奇异，对式(10)进行积分求解的结果一定为原方程的解．至于式(10)稳定平衡点的存在性问题，同样计算梯度矩阵，可得：

(11)

无法直接判断其是否负定．本文研究发现选择平凡点处的J0(如初始求解点处)，采用本文算法计算实际电力系统病态潮流时矩阵-J0-1*J(x*)的特征根均在左半平面，满足具有稳定平衡点的条件．后文算例测试结果也印证了其正确性．

2.2 数值积分方法应用

对式(10)所示微分方程进行求解时，需采用一定的数值积分方法．尽管任何有效的数值积分方法均可应用于通过CNM建立的微分方程的求解[17]，但不同方法的收敛性和收敛速度有很大差异．由第1节中相关内容可知，若采用一阶的显示欧拉法，则本文算法计算过程与普通的牛顿-拉夫逊法基本相同；为完成大规模系统病态潮流的求解，需选用更高阶数、更精确的数值积分方法，这也正体现了采用本文算法与牛顿法计算的区别．

考虑到计算需求，本文主要采用龙格库塔积分方法．该方法分为不同阶数的显式和隐式求解方法，具有不同特性[18]．为保证算法的最优性能，研究中将3种不同阶数的显式龙格库塔法RK1、RK2、RK4分别应用于实际系统算例的求解中，通过收敛性和计算效率的对比进行优选，最终选择RK2作为具体实现时的积分方法．其积分计算式如式(12)所示，其中xi为第i步积分时由电压幅值和相角组成的状态变量，ki(1)、ki(2)为计算的中间向量，ti为虚拟时间变量序列，Δt为固定步长．

(12)

2.3 大规模病态潮流求解流程

图1 大规模病态潮流求解流程

3 算例分析

本文应用了不同规模的电力系统算例，为方便说明，将n节点系统算例简称为casen．文中小规模算例case39，case118，case300，均为IEEE标准算例；大规模系统算例case2383，case3120，case9241以及万节点以上算例case13659则为来自波兰电网与欧洲电网的真实算例．为全面分析算法对于不同规模、不同病态水平的系统潮流求解时的性能，本文采用发电机有功和负荷有功、无功同比例增长的方法构造病态潮流算例，即如式(13)所示．其中λ为负荷增长因子，PL、PL0分别为病态算例和标准算例的所有负荷有功功率，无功功率同理；PG、PG0分别为病态算例和标准算例中所有发电机的有功出力．

PL=λPL0，QL=λQL0，PG=λPG0

(13)

为提高算法效率，本文采用C/C++语言进行算法编程实现，仿真平台为Microsoft Visual Studio 2012；同时采用Matlab的MATPOWER计算库采用牛顿-拉夫逊法进行潮流计算结果正确性对比，仿真平台为Matlab R2015b．计算机配置为Intel®Xeon(R) CPU E3-1230 v5，8GB内存．计算中收敛精度均为10-4．

不失一般性，以case2383为例，本文算法及牛顿-拉夫逊法迭代次数随负荷因子的变化见表1．在负荷水平比较低时，N-R法收敛较快，本文算法由于固定雅可比矩阵的处理而收敛次数较多但稳定于10次左右，此时还未进入病态；随着负荷水平不断提高，收敛越来越困难，在λ=2.0时N-R法已经无法收敛，本文算法仍能够求解，计算的残差曲线如图2所示；λ=2.1时本文算法也无法收敛．

表1 不同负荷水平下算法收敛性(迭代次数)对比

图2 潮流病态求解残差曲线

以最优乘子法计算相应病态潮流的结果作为参考值，评估本文算法的电压和相角误差，进一步验证算法准确性．同样添加变雅克比矩阵的测试结果，针对不同规模的病态潮流算例进行验证，其结果见表2，表中误差值为各个节点的电压(标幺值)和相角(弧度)与参考值作差的绝对值中的最大值．结果显示，电压和相角误差均不超过5×10-3，在工程可以允许的范围内．

表2 病态潮流求解正确性

*注：采用最优乘子法对case3120与case13659的相应病态潮流进行计算时无法收敛，因此以变雅克比法作为参考值．

为进一步验证建立的微分求解模型的合理性，结合2.1节中相应理论分析，针对不同规模的病态算例，计算如式(11)所示梯度矩阵在复平面上最靠右的特征值的实部a，结果见表3．可以看到，所有特征根均位于负半平面，微分方程组具有稳定平衡点．

表3 梯度矩阵特征根分析

以上算例分析结果能够验证本文算法应用于病态潮流计算时的正确性．对于大规模系统而言，算法计算效率同样需要关注．数值积分方法的步长对于常微分方程组求解的收敛速度有较大影响，对于本文方法同样如此．在应用算法时为达到效率的最优化，应选取使得迭代次数最少的积分步长，但并未有研究能够对积分步长的选取进行理论指导，因此本文对不同步长时算法的收敛速度进行了实际测试与比较．求解时在Δt=0～1.5之间设置间隔为0.01的150个点，对Δt值在1.0附近波动时的算法迭代次数进行测试，结果如图3所示．迭代次数上限设置为100次，达到该次数即视为不收敛．可以看出对于不同规模系统具有相似规律，在积分步长过小或过大时无法收敛，最优积分步长在0.8～1.0之间．

图3 不同积分步长时算法收敛速度

为进一步提高整体效率，本文对于因子分解及前代回代计算等步骤均采用了Intel MKL[19]计算库中的相应函数，不同规模系统病态潮流计算的整体耗时结果见表4．可见算法对于上千节点大规模系统的病态潮流计算均有较高效率，对于超过万节点的系统也在1s内能够给出病态潮流计算结果．

表4 算法整体效率

4 结 论

本文完成了基于微分求解模型的大规模电力系统病态潮流计算，主要结论如下：1)基于连续牛顿法，非线性方程组的求解可以转化为常微分方程组的数值积分求解．2)建立电力潮流的微分求解模型，并进行相应简化处理，使其适用于大规模电力系统的病态潮流求解．3)将数值积分方法应用于求解中，得到整体积分计算式，并根据电力网络特点构建整体求解算法流程．4)结合具体系统算例进行分析，通过对N-R法以及最优乘子法潮流计算的结果对比，本文算法求解病态潮流时的有效性和准确性得以证明；同时结果显示算法对于大规模系统有较高效率，是可以应用于实际场景的实用性算法．