提高电力系统潮流计算收敛性方法的比较分析

范翔，严正，赵文恺，范金燕，王虹富

（1．上海交通大学电力传输与功率变换控制教育部重点实验室，上海 200240；2.上海交通大学应用数学系，上海 200240；3．中国电力科学研究院，北京 100192）

随着输电规模的扩大，电网将形成互联大区域的复杂结构，系统运行工况变得多样化。电网方式安排人员在得到合理的系统运行工况前，会遇到许多潮流难以收敛的情况，而此时传统牛顿法计算失败，难以为后续的稳定分析及方式调整提供有效信息，因此，研究提高潮流计算收敛性的方法具有普遍的现实意义和重要的实用价值。

造成传统牛顿法或者PQ解耦法潮流计算难以收敛的主要原因是由于迭代过程中的雅可比矩阵奇异或接近奇异[1]。目前能够用于提高大规模电力系统潮流计算收敛性的典型方法主要有：最优乘子法[1-5]、张量法[6-9]和自适应LM方法[10]。文献[1]提出最优乘子法计算潮流，具有永不发散的特点，且与现有的潮流程序接口简单、易于实现，是一种通用提高潮流收敛性的方法，在各类难收敛潮流计算中效果较好；文献[6-7]将张量法引入潮流计算，考虑了潮流方程泰勒展开近似二阶项的影响，在部分区域重负荷情况下，张量法能够取得到较好解；文献[10]引入自适应LM方法计算潮流方程，该方法在迭代过程中雅可比矩阵始终非奇异，通过算例仿真验证了自适于应LM方法，减小了对初值选取的要求，扩大了收敛范围；文献[11]基于优化技术的潮流方法[11]，将求解潮流方程组转化成计算优化模型，使得问题复杂化。基于大范围收敛同伦法提高潮流计算收敛性，扩大了潮流方程组的收敛范围，但计算量大，未引入商业应用[12-14]。

本文研究在常规潮流基础上比较容易实现的3类方法：张量法、最优乘子法和自适应LM方法，及其数值过程。总结了各方法进行电力系统潮流计算的特点。基于数值过程，介绍了3种方法关键步的稀疏实现，比较了3种方法相对牛顿法单步迭代的额外计算开销。在1个标准系统和2个实际系统中进行仿真，仿真计算结果表明，相比其他2种方法，自适应LM方法应用于实际大系统需结合稀疏技术。

1 提高潮流收敛性的3种方法

令潮流方程为

式中，状态变量x=[V，θ]T。当潮流方程雅可比矩阵J（x）奇异时，常规牛顿法失效，需采用提高潮流收敛性的方法。

1.1 最优乘子法

所谓最优乘子法，即在每一步计算状态变量的修正量Δx之后，不直接用Δx去修正状态变量x，而是乘以一个标量乘子λ去修正，即

最优乘子法的核心在于乘子λ的计算，思想是在牛顿迭代方向上寻求最优步长，即

式中：λk为第k步迭代时的最优乘子；dk为第k次迭代的牛顿步。若在当前牛顿方向上搜寻不到最优步长，则λ→0，方法停止。若潮流结果收敛，则d→0，当方法收敛时λ→1。由此可见，引入最优乘子可有效避免因迭代步步幅过大而导致方法振荡或发散。最优乘子方法的具体计算流程参见文献[1-2]。

1.2 张量法

张量法潮流计算的实质为含二阶项潮流计算。通过对潮流方程的二阶展开项近似，计算出张量步，修正牛顿迭代步。

1.2.1 插值张量法

在当前迭代点xk处将式左侧展开至二阶，即

若潮流方程采用直角坐标形式，则张量Ak为常值三维矩阵。由文献[12]知，Ak可用前述迭代点表示为其中，sk=xi-xk，i∈[1，k-1]且使得{sk}线性无关；ak为已知向量。将Ak带入式（4）计算迭代步dk，得方程

（i=1，…，p）的二次方程组，即

式（6）称为张量方程组，若仅取前一次迭代点xk-1做插值，则式（6）仅有一个二次方程和一个待求参数β，回代入式（5），得

若式（7）不存在实数解，则采用正交变换将关于迭代步dk的最小二乘min‖F（xk+dk）‖2转换为关于张量参数β的最小二乘，计算出dk做为迭代步[10]，该方法依然要求雅可比矩阵非奇异。方法具体流程详见文献[10]。

1.2.2 直接张量法

设实际迭代步为d=dN+dT，其中dN为牛顿步，dT为张量步，且参考直接张量法原始文献[6]知，‖dN‖2≫‖dT‖2。将迭代步带入式（4），令其等于0，计算得到

1.2.3 数值特点

在良态迭代过程中，牛顿法模型准确，直接张量法的假设条件可以被满足。当潮流方程难以收敛，按照含二阶项潮流计算理论，牛顿法线性化模型的误差较大，牛顿步不准确，需引入张量步进行修正。根据张量修正步dT的计算表达式知，在雅可比矩阵趋近奇异时dT的幅值较大从而‖d‖2也较大，也会出现与常规牛顿法相同的振荡发散问题。且直接张量法缺乏收敛性理论证明。

两种张量法提高收敛性的特点为张量步对牛顿步的修正使得迭代步d变为

易知张量法收敛的充分条件是‖F+C‖2为雅可比矩阵最小奇异值σmin的高阶无穷小。由此可知，C所起补偿作用是在雅可比矩阵条件数较大或接近奇异时减小部分潮流方程的失配量，防止迭代步幅过大而发散。

综上，两类张量法的特点是：①引入二阶补偿项对潮流偏差量进行补偿，当迭代过程中雅可比矩阵接近奇异或条件数很大时，减小迭代步长；②基于插值的张量法由文献[9]奠定了数学基础，而直接张量法缺乏相关理论依据。

1.3 自适应LM方法

将式（1）左侧在当前迭代点处一阶展开为

式中，dk=xk+1-xk。

潮流方程的最小二乘模型为

当式（11）所得解x~满足G（x~）=0时，x~即为潮流方程的解。

将式（10）代入式（11）并引入步长约束，得到最初计算自适应LM方法迭代步的模型，即

式中，参数δ按照一定方式更新[16]。由文献[16]可得式（12）的解为

式中，μk为自适应LM方法阻尼因子，μk=αk‖F（xk）‖，其中αk为非负调节因子。自适应LM方法的基本计算流程[17]为

步骤1平启动，迭代步数k置1，设置α1和常数m及收敛精度ε，使得α1＞m，0＜p0＜p1＜p2＜1。

步骤2计算μk，按式（13）计算dk。

步骤3计算τk

选择是否接受dk

步骤4自适应因子的调整

步骤5用潮流收敛判据‖F（xk）‖2＜ε判别潮流收敛与否，不收敛则返回步骤2继续迭代，至最大迭代次数方法停止输出结果。

对比最优乘子法，当雅可比矩阵条件数很大或接近奇异时，牛顿步非下降方向，最优乘子强制为0，由于不能改变迭代方向，方法停止在潮流失配量较小的近似潮流解上。而式（13）中自适应LM方法通过改变阻尼因子μk来调整迭代步步幅和方向。因此在计算难收敛潮流工况时，自适应LM方法较最优乘子法有更好的适应性。另外，自适应LM方法在迭代过程中求逆矩阵JTJ+μI始终非奇异，因而自适应LM方法较具有更好的数值稳定性和收敛性。

综上所述，自适应LM方法在潮流计算有两个特点：①在计算过程中保持求逆矩阵JTJ+μI对称正定，不会遇到奇异问题，因此自适应LM方法具有较好的数值稳定性；②自适应LM方法通过改变阻尼因子μ的值，来同时改变迭代步长和迭代方向，对难收敛潮流具有更好的适应性。同时方法可收敛至潮流方程的最小二乘解。

2 各方法的稀疏实现及计算量比较

2.1 稀疏实现

张量法和最优乘子法的求逆矩阵结构与牛顿法雅可比矩阵相同，故可采用与常规牛顿法相同的矩阵LU分解（LU decomposition/三角分解）及前代回代提高计算效率。因此，张量法和最优乘子法均可计算大规模系统潮流。

由式知计算自适应LM方法迭代步dk的求逆矩阵JTJ+μI结构不同于牛顿法雅可比矩阵J，且较J矩阵稀疏度有所下降，但其仍为高度稀疏的矩阵，因此不能直接采用牛顿法的稀疏处理。

JTJ+μI在结构和数值上均对称，且为正定矩阵，可以采用稀疏Cholesky分解[18]，来求解稀疏线性方程组。由于矩阵JTJ+μI与系统导纳矩阵稀疏结构不同，所以不能沿用基于导纳矩阵的排序结果。此时采用近似最小度AMD（approx-imate minimum degree）排序算法[18]对JTJ+μI矩阵排序，减少Cholesky分解时的注入元个数。由于整个迭代过程中矩阵JTJ+μI的结构保持不变，故只需在首次迭代时对矩阵JTJ+μI进行排序即可。

2.2 计算量比较

基于插值张量法的每步迭代，当p=1时，求解式（6）时需求解一次大规模稀疏线性方程组，在求解式（7）时还需求解一次大规模线性方程组。由于单步迭代中求逆矩阵因子表固定，所以p=1时插值法的单步计算量约为普通牛顿法单步计算量的2倍。

对于直接张量法的每步迭代，计算dN和dT需各求解一次大规模稀疏线性方程组，因而单步计算量约为普通牛顿法的2倍。

最优乘子法仅在计算牛顿步时求解一次大规模稀疏线性方程组，因此最优乘子法单步迭代计算量约等于普通牛顿法。

自适应LM方法的单步迭代中，需要由式（13）来求解迭代步，而式（13）中的矩阵JTJ与状态估计中的系数矩阵有类似的结构，因此可以认为自适应LM方法单步计算量比普通牛顿法略有增加，本文的算例仿真部分验证了该结论。

综合第1节的数值机理分析，3种方法的单步计算量和收敛特点如表1所示。

表1 最优乘子法、张量法和自适应LM方法单步计算量及收敛特点汇总Tab.1 Computational costs per step and convergent characteristic of optimal multiplier method，tensor method and self-adapted LM method

3 算例仿真

仿真平台为主频3.0GHzCore2 CPU，4G内存的PC机，操作系统为Win7，采用Matlab2013a和C++的混合编程方式，C++实现与BPA的数据接口，Matlab实现各方法以及测试各算例。算例中的具体信息如表2所示。

表2 测试系统数据Tab.2 Data of test system s

比较几种方法的单步迭代时间，其中自适应LM方法调用Matlab函数chol和amd实现稀疏Cholesky分解及AMD排序。

算例1采用IEEE39节点系统构造重负荷难收敛潮流。

算例2采用来源于mat power4.1[19]波兰3375wp测试系统，采用mat power自带潮流求解器因迭代初始点雅可比矩阵奇异无法收敛。

算例3源于华东电网2004年夏高的BPA数据，加重系统部分区域负荷，至BPA潮流程序不收敛。模拟局部重负荷难收敛工况。

对以上3个算例分别采用张量法、最优乘子法和自适应LM方法计算各系统的潮流，验证并比较3种方法的收敛特性。

算例4采用国网21 479节点数据，用于比较文献[10]的LM方法和本文的自适应LM方法，在大规模系统潮流计算中的适应性。

对于重负荷系统，有可能出现PV-PQ节点的频繁转换，使得潮流计算的收敛域变窄[20]，在某种程度上恶化了潮流程序的收敛性。由于本文旨在对潮流计算的3种主要方法的收敛性和计算量进行分析和比较，因此，算例中不考虑系统中PV-PQ节点转换的情况。

在比较各方法的收敛特性之前，先通过华东04系统来比较各方法的单步计算量，结果如表3所示。

表3 单步迭代计算时间比较Tab.3 Comparisonin of computing time per iteration between each method s

由表3可知，各方法的单步计算时间与常规牛顿法的单步时间之比与表1中给出的结论基本一致。同时，也验证了用稀疏Cholesky分解和AMD排序能够显著提高自适应LM迭代步的计算效率。

3.1 算例1

采用IEEE39节点系统构造重负荷难收敛潮流。良态情况下各方法所得结果与BPA所得潮流结果比较，电压相角（rad）和幅值（p.u.）的绝对偏差为10-7，证明所写程序无误。

将非发电机节点的有功负荷扩大1.3倍，无功负荷扩大1.4倍，BPA牛顿法计算该潮流不收敛。采用张量法、自适应LM方法和最优乘子法计算潮流方程组迭代曲线如图1和图2所示。

自适应LM方法参数取值为：α1=0.1，m=10-5，p0=10-4，p1=0.75，p2=0.95。各方法的收敛准则均为潮流偏差‖F‖2≤10-5。

图1 两类张量法迭代曲线Fig.1 Iteration curves of two tensor methods

从图1可见，直接张量法在潮流偏差较大时陷入振荡，插值张量法补偿项取得的结果更好，但在得到10-2数量级的近似潮流解时，插值张量法陷入数值振荡。两种方法近似二阶补偿项C的补偿效果插值张量法优于直接张量法。

从第15次迭代可以看出，插值张量法二阶补偿项Cp的修正使得潮流偏差量减小为最小奇异值的高阶无穷小。而直接张量法二阶补偿项Cd的补偿性能不佳，导致步幅仍呈现增大趋向。第21次迭代时，补偿项不能补偿潮流偏差至合适值，此时最小奇异值为：0.03，因而产生振荡。张量法虽在迭代中得到最小潮流偏差量近似解，但因补偿项不能跟上雅可比矩阵奇异值减小的速度最终振荡，难以通过现有的方法将迭代停止在合适的潮流解上。

图2 最优乘子法和自适应LM方法迭代曲线Fig.2 Iteration curve of self-adaptive LM method and optimal multiplier method

从图2可以看出，初始迭代段自适应LM方法由于潮流偏差较大阻尼因子很大，因而迭代步-（JTJ+μI）-1JTF呈现出最速下降步-1/μJTF折线下降的特点。迭代至第5步时两种方法具有十分接近的潮流偏差量。第5、6步迭代时最优乘子为6.1×10-4和3.3×10-5；最优乘子法第5步迭代时雅可比矩阵的最小特征值为3.5×10-3，第6步迭代时雅可比矩阵的最小特征值为8.2×10-4；由此可见，雅可比矩阵趋近奇异，最优乘子趋近0，方法停滞。自适应LM方法在第5次迭代时求逆矩阵最小特征值为0.09，第6次迭代时求逆矩阵最小特征值为6.0×10-3，由此可验证自适应LM方法阻尼因子绕开了雅可比矩阵在第5次迭代接近奇异的区域，获得了潮流偏差量更小的解。

3.2 算例2

波兰3375wp潮流数据取自波兰2007年冬季晚高峰潮流断面。应用matpower自带潮流求解器计算该系统潮流发散。由于平启动迭代时雅可比矩阵趋近奇异，最优乘子法停滞，张量法均失效。系统的雅可比矩阵趋近奇异，最小特征值为0，最小奇异值为0。

由于首次迭代雅可比矩阵趋近奇异，最优乘子法停滞，张量法失效，只有自适应LM方法顺利计算出了该系统的潮流。自适应LM方法计算该系统潮流耗时1 s，迭代曲线如图3所示。

图3 自适应LM方法迭代功率失配量曲线Fig.3 Variation of power mismatches via self-adaptive LM method

从图3可以看出，经过10次迭代后，自适应LM方法平稳收敛至‖F‖2=7.2×10-8。

采用最优乘子法计算该系统时，首次迭代遇到雅可比矩阵趋近奇异，导致最优步长趋近0，方法停滞。

采用张量法计算时，由于雅可比矩阵趋近奇异，牛顿步无法计算，因此张量法失效。

采用自适应LM方法计算得到的电压如图4所示。

图4 电压幅值柱形图Fig.4 Voltage amplitude histogram

从图4可以看出，电压幅值在允许范围内。计算所得，最小相角为-37.29°，最大相角为3.21°。该潮流解未呈现出病态情况，表明自适应LM方法计算求得真实潮流解。

3.3 算例3

该算例采用2004年夏季华东电网实际潮流断面数据。原始数据为良态潮流数据，将SH区域的有功负荷扩大21.98%，无功负荷扩大60%；JS区域的有功负荷扩大7.95%，无功负荷扩大6.82%，以此模拟部分区域重负荷难收敛工况。此时BPA潮流程序计算发散，自适应LM方法、最优乘子法所得迭代曲线如图5所示。

图5 自适应LM方法及最优乘子法迭代曲线Fig.5 Iterative curves of self-adaptive LM method and optimal multiplier method

采用插值张量法计算该系统潮流，获得初始迭代点后，首次迭代雅可比矩阵奇异，由于初值非足够接近潮流解故不能获得牛顿步，即不能通过牛顿步求解张量修正步。采用直接张量法亦存在此问题。

取最优乘子法停滞时的雅可比矩阵Jop_end和功率偏差向量Fop_end，计算得到Jop_end的正交基Qop_end后发现同样，取自适应LM方法终止时的雅可比矩阵JLM_end和功率偏差向量FLM_end，计算得到FLM_end的正交基QLM_end后发现，因此，自适应LM方法得到了该系统潮流的精确最小二乘解。

3.4 算例4

本算例旨在比较文献[10]的方法与本文的自适应LM方法的收敛性。由于迭代过程中雅可比矩阵出现奇异，最优乘子法和两类张量法均失效。

[10]提出的LM方法，其阻尼因子为μ=0.001‖F‖2与功率偏差量呈线性关系。即潮流偏差大时阻尼因子大，潮流偏差小时阻尼因子小，具备提高潮流收敛性的LM方法阻尼因子的特点。但是该阻尼因子在当前迭代步不能调节，即不能获得当前迭代点的最优下降方向，从而增加迭代次数，计算速度慢。

本文的阻尼因子采用

式中，α为自适应因子，通过类似信赖域的方法进行调节，阻尼因子与功率偏差量呈非线性关系，能够获得当前迭代点合适的迭代方向和迭代步幅。通过国网21479节点算例可验证（见图6）。

图6 2种方法计算国网21479系统迭代曲线比较Fig.6 Comparison of two methods in terms of iterations for state grid-21479 system

从图6可以看出，采用自适应因子调节的阻尼因子更能适应大规模系统潮流计算。

4 结论

本文分析了张量法、最优乘子法、和自适应LM方法的数值特点，以此为基础比较了3种方法的收敛特点及单步计算量。并对1个标准系统以及2个实际系统进行仿真计算。

（1）在重负荷情况下，插值张量法利用二阶项能够较好地补偿功率偏差，对难收敛潮流具有较好地适应性；然而当雅可比矩阵趋近奇异时，易振荡；直接张量法补偿项效果不佳。

（2）当系统潮流呈现重负荷而难收敛时，最优乘子法最优乘子迅速减小至0，收敛到近似最小二乘，近似解处偏差量仍较大；当初始点难以确定甚至雅可比矩阵奇异时，最优乘子法失效。

（3）采用自适应LM方法提高潮流计算收敛性较最优乘子法和张量法有更好的鲁棒性和数值稳定性。由于阻尼因子的存在，自适应LM方法迭代过程中不会发生因雅可比矩阵奇异而导致的发散问题。自适应LM方法较小依赖初值，对初值难设定、重负荷等工况具有很好的适应性。

参考文献：

[1]Iwamoto S，Tamura Y. A load flow calculation method for ill-conditioned power systems[J]. IEEE Trans on Power Apparatus and Systems，1981，100（4）：1736-1743.

[2]王宪荣，包丽明（Wang Xianrong，Bao Liming）.极坐标系准最优乘子病态潮流解法研究（The study of ill-conditioned load flow using quasi-optimal factor in polar coordinates）[J].中国电机工程学报（Proceedings of the CSEE），1994，14（1）：40-45.

[3]Overbye T J.Computation of a practical method to restore power flow solvability[J].IEEE Trans on Power Systems，1995，10（1）：280-287.

[4]石飞，赵晋泉，王毅（ShiFei，Zhao Jinquan，Wang Yi）.计及发电机无功约束的最优乘子潮流计算方法比较（A comparison of Newton -Raphson optimal multipliers load flow methods considering generator reactive power limitation）[J].电力系统保护与控制（Power System Protection and Control），2009，37（2）：6-10.

[5]胡泽春，严正（Hu Zechun，Yan Zheng）.带最优乘子牛顿法在交直流系统潮流计算中的应用（Application of Newton load flow methods with optimal multiplier for AC/DC power system）[J].电力系统自动化（Automation of Electric Power Systems），2009，33（9）：26-31.

[6]Salgado R S，Zeitune A F.Power flow solutions through tensor methods[J].IET Generation，Transmission&Distribution，2009，3（5）：413-424.

[7]林济铿，吴鹏，袁龙，等（Lin Jikeng，Wu Peng，Yuan Long，et al）.基于张量法的电力系统潮流计算（Power flow computation based on tensor methods）[J].中国电机工程学报（Proceedings of the CSEE），2011，31（34）：113-119.

[8]Bouaricha A，Schnabel R B.Tensor methods for large sparse systems of nonlinear equations[J].Mathematical Programming，1998，82（3）：377-400.

[9]Schnabel RB，Frank PD.Tensor methods for nonlinear equations[J].SIAM Journal on Numerical Anal-ysis，1984，21（5）：815-843.

[10]Lagace PJ.Power flow methods for improving convergence[C]//38th Annual Conference of IEEE Industrial Electronics.Montreal，Canada，2012：1387-1392.

[11]覃智君，侯云鹤，吴复立（Qin Zhijun，Hou Yunhe，Wu Fulix）.大规模交直流系统潮流计算的实用化模型（Practi-cal model for large-scale AC-DC system power flow calcu-lation）[J].中国电机工程学报（Proceedings of the CSEE），2011，31（10）：95-101.

[12]周佃民，廖培金（Zhou Dianmin，Liao Peijin）.电力系统病态潮流的同伦方法求解（Homotopy method for illconditioned power system load flow calculation）[J].电力系统及其自动化学报（Proceedings of the CSU-EPSA），1999，11（5-6）：67-71.

[13]陈礼义，戴宏伟，张琦鹏（Chen Liyi，DaiHongwei，Zhang Qipeng）.一种大范围收敛的电力系统潮流算法——同伦延拓法（A large range convergent method for power system load calculation—Homotopy calculation method）[J].电力系统及其自动化学报（Proceedings of the CSUEPSA），1993，5（1）：67-74，30.

[14]孙小舟（Sun Xiaozhou）.使用“BPA电力系统分析程序”的体会（Experiences on using‘BPA power system analysis program’）[J].电力建设（Electric Power Construction），1998，19（11）：34-35.

[15]周国标，宋宝瑞，谢建利.数值计算[M].北京：高等教育出版社，2008.

[16]Dennis JE Jr.Numerical Methods for Unconstrained Optimization and Nonlinear Equations[M].北京：科学出版社，2009.

[17]Fan Jinyan，Yuan Yaxiang. On the quadratic convergence of the Levenberg-Marquardt method without nonsingularity assumption[J]. Computing，2005，74（1）：23-39.

[18]Davis T A.DirectMethods for Sparse Linear Systems[M].Philadelphia：SIAM，2006.

[19]Zimmerman R D，Murillo-Sanchez C E，Thomas R J.MATPOWER：Steady-state operations，planning，and analysis tools for power systems research and education[J].IEEE Transon Power Systems，2011，26（1）：12-19.

[20]赵晋泉，江晓东，张伯明（Zhao Jinquan，Jiang Xiaodong，Zhang Boming）.潮流计算中PV-PQ节点转换逻辑的研究（Study on PV-PQ bus type switching logic in power flow computation）[J].中国电机工程学报（Proceedings of the CSEE），2005，25（1）：54-59.