0 引言
交直流电网互联规模不断扩大、可再生能源大规模并网，在改善电网大规模电能输送能力和优化电源结构的同时，也使电力系统运行中随机性和不确定性因素日益增加，给电力系统电压稳定性带来诸多挑战[1-5]。因此，研究各种随机性和不确定性因素影响下的电力系统电压稳定性对电力系统安全稳定运行具有极其重要的意义[6-9]。
静态电压稳定域（Static Voltage Stability Region,SVSR）是分析和评估电力系统电压稳定性的重要工具[10-13]，它已在研究强随机性和不确定性因素影响下的电力系统电压稳定性中发挥了巨大的作用，但构建SVSR关键在于对其边界的搜索。超平面近似作为目前最为通用的静态电压稳定域边界（Static Voltage Stability Region Boundary,SVSRB）构建方法，已在工程中得到了应用[14-16]。该方法核心思想是在功率注入空间内搜索系统所有可能运行方式下的鞍结分岔（Saddle Node Bifurcation,SNB）点，借助这些SNB点来近似SVSRB。因此，构建SVSRB关键是对SNB点的快速、准确搜索。
目前，电力系统SNB点的搜索方式主要有连续潮流（Continuation Power Flow,CPF）法[17-18]、最优潮流（Optional Power Flow,OPF）法[19]和直接法[20-21]。CPF法以系统当前运行点为始点，沿指定功率增长方向采用预测-校正方式追踪系统SNB点，该方法虽可获得较高精度的SNB点，但系统规模会严重制约该方法的计算效率。OPF法以系统负荷裕度最大为目标函数、潮流方程为约束条件，通过求解该优化模型获得电力系统的SNB点，但优化算法需要通过寻优以逐次逼近最优解，其求解过程相对耗时，且“维数灾”问题难以解决，不利于区域互联电力系统SNB点的快速搜索。直接法通过求解一组表征SNB点特性的非线性方程组来搜索SNB点，该方法对初值较为敏感且计算量较大，系统规模的扩大会导致该方法的计算量和存储空间急剧增加，故该方法也难以适用于大规模区域互联电力系统SNB点的快速搜索。为提高SVSRB的搜索效率，文献[22-23]分别提出了基于边界追踪的预测-校正算法和基于边界追踪的优化算法，以实现SVSRB的快速构建。但如何在实现区域互联电力系统SVSRB高精度构建前提下，进一步提高SVSRB的计算效率仍具有十分重要的研究价值。
相对于CPF法和OPF法，直接法通过求解一组表征SNB点特性的非线性方程组便可获取系统的SNB点，该方法计算原理简单、计算精度高，且在求解SNB点时涉及大量的线性方程组计算，虽然其方程组的计算量和计算复杂度较高，但易于通过GPU并行加速以提高SVSRB计算效率[24]。因此，本文以直接法为基础，提出一种基于CPU-GPU异构的电力系统静态电压稳定域边界并行计算方法。该方法依据SVSRB拓扑性质，采用沿边界追踪SVSRB的搜索方式，将已求SNB点处系统信息作为直接法搜索相邻SNB点的初值，以克服直接法对初值敏感的不足；再结合CPU-GPU异构平台实现SNB点的并行求解，以降低直接法计算量大、计算复杂度高的缺陷，从而提升SVSRB搜索效率；最后通过多组测试系统算例，对所提SVSRB并行计算方法的准确性和有效性进行分析、验证。
1 基于直接法的鞍结分岔点计算方法
由文献[25-27]可知，电力系统的SVSRB主要是由SNB点和极限诱导分岔（Limit Induced Bifurcation,LIB）点构成的，而LIB点可在SNB点搜索的基础上，利用参考文献[28]来判断是否到达，因此，SNB点的搜索至关重要。本文在构建SVSRB时重点关注SNB点的搜索，将SVSRB的搜索等效为搜索系统在各种可能运行方式下的SNB点。而直接法作为一种求解SNB点的常用方法，其原理简单，便于计算，已被用于电力系统SNB点的求解[20-21,28]。因此，本文采用直接法来搜索电力系统的SNB点。
1.1 直接法基本原理
直接法搜索电力系统SNB点的核心思想就是求解一组能够反映系统SNB点特性的非线性方程组，该非线性方程组的表达式为
式中，f(x,λ)=0为含参数的潮流方程，x为状态变量，l为负荷裕度；fx(x,λ)·r=0为潮流方程所对应Jacobi矩阵发生奇异的方程，fx(x,l)为潮流Jacobi矩阵，r为零特征值所对应的右特征向量，其维度与潮流方程一致；rp-1=0为规范化方程，rp代表r中的第p个元素。对于一给定电力网络，若系统总节点数为N，PQ节点共npq个，PV节点共npv个，令n=2npq+npv，则上述非线性方程组的维度为2n+1。
采用直接法求解电力系统SNB点的过程，实质就是采用牛拉法求解式（1）所示的非线性方程组，其求解过程中的修正方程组为
式中，fxxr∈Rn×n，其任意第l行、第d列的元素可表示为
，其第p个元素为1，其余元素均为0；fλ∈Rn为含参数的潮流方程对λ的偏导数。
由式（2）可求得修正量Δx、Δr和Δl，将所得修正量代入式（3）更新变量x、r和l。
设定收敛阈值ε，将式（3）更新后的x、l和r代入式（1），若满足收敛条件‖f(x,l)‖≤ε和‖fxr‖≤ε，则迭代完成，SNB点处x、r和l 求解完毕；否则，继续迭代，直至f(x,l)和fxr满足上述收敛条件为止。
1.2 修正方程组的降阶处理
由式（2）可知：修正方程组式（2）的维数是潮流方程维数的两倍，在迭代求解过程中，计算量较大，导致整个SNB点的求解耗时较高。为降低式（2）中修正量求解的复杂度和计算量，文献[20]将式（2）所示的2n+1维修正方程组通过换元变换拆分成式（4）～ 式（7）所示的四组同系数n+1维线性方程组。
各修正量与α、β、ξ和η关系为
式中，Δxp为薄弱节点电压幅值变化量，用于临界点处x、r和l的初值预测，其计算公式为
通过求解降阶后的线性方程组式（4）～式（7），并根据式（8）所示各修正量与α、β、ξ和η的关系，先求解修正量Δx、Δr和Δl，然后将Δx、Δr和Δl代入式（3），可对待求变量x、r和l进行修正，以提高修正量求解效率，减少数据所需存储空间，提高直接法搜索SNB点的计算效率。
2 基于直接法的SVSRB并行计算
采用第1节所述的直接法在功率注入空间内搜索大量SNB点，并借助这些SNB点即可近似系统SVSRB。然而，直接法对初值非常敏感，若采用传统的SVSRB搜索方法（简称“射线法”），不当的初值选取可能导致直接法求解结果不收敛或收敛耗时较长。同时，系统规模的扩大将导致直接法的计算量急剧增加，不利于大规模电力系统SVSRB的快速搜索。针对上述不足，本文提出一种基于边界追踪的直接法以实现SVSRB的并行计算。
2.1 基于边界追踪的直接法
边界追踪法的核心思想是利用SVSRB上相邻SNB点的系统运行信息变化较小这一特点，将已求SNB点的状态量x、右特征向量r和负荷裕度l作为下一相邻待求解SNB点的初值，以充分利用各SNB点信息，提高SVSRB构建效率。基于该思想，本文将直接法与边界追踪法相结合，提出一种基于边界追踪的SVSRB快速搜索的直接法，以有效利用各SNB点的系统运行信息，克服直接法初值选取困难的不足。由于基于边界追踪的直接法可以很好地保证求解SNB点时给定初值的精度，因此，有利于加快直接法计算SNB点的收敛速度，提高SVSRB的搜索效率。
图1给出了基于边界追踪的直接法搜索SVSRB的基本原理。首先给定初始功率增长方向d0，d0对应的功率增长方向角q0=arctan(ΔPj0/ΔPi0)，ΔPi0、ΔPj0分别为d0中对应于节点i、j的功率注入增长分量。以基态为起始点，根据近似二次曲线预估SNB点初值[20]，采用第1节所述的直接法计算功率增长方向d0下的SNB点0，SNB点0处的状态变量、负荷裕度和右特征向量分别为x0、l0和r0。改变d0中节点i、j的有功注入增长分量ΔPi0和ΔPj0，得到新功率增长方向d1为
图1 基于边界追踪的SVSRB搜索示意图
Fig.1 Exploring SVSRB using boundary tracking method
式中，
。首先将SNB点0的状态变量x0、负荷裕度l0和右特征向量r0作为求解d1下SNB点1的初值；然后，采用直接法计算SNB点1处的x1、l1和r1。同理，以SNB点1处的x1、l1和r1作为求解d2下SNB点2的初值，再采用直接法计算SNB点2处的x2、l2和r2，则可实现沿q 减小方向下由SNB点0至SNB点2的搜索。再沿q增大的方向，采用同样的求解方法即可实现由SNB点0至SNB点4′的快速搜索，进而快速构建二维有功注入空间的SVSRB。
2.2 直接法并行计算SVSRB的预处理
采用边界追踪的直接法很好地解决了直接法搜索电力系统SNB点时初值选取困难这一瓶颈。然而，直接法求解SNB点时，牛顿迭代修正方程组式（2）的维数是潮流方程组维数的两倍，计算量大、计算复杂度高，严重影响了SVSRB构建效率。尽管对修正方程组式（2）的降阶处理有效降低了修正方程组维数，但降阶处理使得需要求解的线性方程组数量增多，搜索SNB点的计算量依然很大，不利于SNB点的快速搜索。为此，本文基于CPU-GPU异构平台，将直接法求解式（4）～式（7）的四组n+1维线性方程组采用带Jacobi+ILU(0)两阶段预处理的稳定双共轭梯度法（BICGSTAB）进行求解，并借助GPU对该求解过程进行并行加速，以降低直接法搜索电力系统SNB点的计算量和复杂度，提高SVSRB搜索效率。
利用GPU的并行计算能力加速迭代法求解线性方程组是GPU的主要应用之一[29-31]，借助GPU的并行计算优势，并根据系数矩阵B既不对称也不正定，且具有高度稀疏性的特点，参考文献[24]，本文采用Jacobi预处理与ILU(0)预处理相结合的Jacobi+ILU(0)两阶段预处理方法，并结合BICGSTAB迭代法求解线性方程组式（4）～式（7），以提高修正量求解速度。
以式（4）为例，本节简要介绍Jacobi+ILU(0)预处理的BICGSTAB迭代法求解过程。将线性方程组式（4）简记为
式中，B为系数矩阵；α为待求向量；b为等式右侧列向量。对式（12）进行Jacobi+ILU(0)预处理，处理后的线性方程组等效为
式中，J为Jacobi预处理矩阵；L、U分别为ILU(0)分解得到的上、下三角稀疏矩阵。选取L为预处理矩阵，对式（13）进行变换，得到
2.3 基于直接法的SVSRB并行计算流程实现
采用基于BICGSTAB迭代法和Jacobi+ILU(0)两阶段预处理的线性方程组求解方法可以加快式（4）～式（7）的求解效率。为进一步提升直接法搜索SNB点的计算效率，本文采用CPU-GPU异构平台，将直接法搜索SNB点过程中计算量较大、计算耗时占比高的修正量求解部分由GPU完成，即借助GPU的并行加速技术对BICGSTAB迭代法进行加速求解，以提高修正量的求解效率，而其他逻辑性较强，计算量较小部分由CPU完成。具体而言，在CPU中，通过牛拉法对表征SNB点特性的非线性方程组式（1）线性化处理，得到相应的修正方程组式（2），并对式（2）进行降阶处理，得到如式（4）～式（7）所示的四组同系数线性方程组；然后将系数矩阵B、线性方程组（4）～式（7）右侧列向量、生成的预处理器J、L和U，以及其他相关参数从CPU移植到GPU中，在GPU中，采用带Jacobi+ILU(0)两阶段预处理的BICGSTAB迭代法对式（4）～式（7）进行并行求解，得到α、β、ξ和η；接着将α、β、ξ和η传输到CPU中，在CPU中，利用各修正量与α、β、ξ和η的关系，由式（8）得到修正量Δx、Δr和Δl，并按式（3）对待求量x、r和l进行修正，将修正后的x、l和r代入式（1），重复牛顿迭代过程，直到f(x,l)和fxr满足‖f(x,l)‖≤ε和‖fxr‖≤ε，则迭代完成，x、r和l求解完毕，至此，则可实现SNB点并行求解，从而提高SVSRB的搜索效率。
图2进一步给出了二维有功注入空间内所提基于CPU-GPU异构的静态电压稳定域边界并行计算的流程。
图2中采用CPU-GPU异构平台计算SNB点的具体过程如前所述，牛顿迭代求解SNB点的初值x0、l0和r0由SVSRB上相邻已求SNB点的系统运行信息给定。因SVSRB上相邻SNB点距离较近，系统运行信息变化较小，故采用基于边界追踪的直接法可以保证直接法搜索SNB时对初值的高质量选取，有利于加快直接法搜索SNB点的收敛速度，克服直接法初值不易选取的不足；而对于直接法计算量大、计算复杂度高这一缺陷，采用CPU-GPU异构平台，利用GPU并行技术对直接法搜索SNB点时计算量较大的四组同系数线性方程组式（4）～式（7）进行并行加速求解，大幅降低了直接法搜索SNB点的计算量与计算复杂度。上述两种处理方式的有机结合，有效解决了直接法初值敏感、计算量大这两大不足，实现了SVSRB上各SNB点的快速并行求解。
图2 基于CPU-GPU异构的SVSRB并行计算流程
Fig.2 Flow chart of the proposed CPU-GPU heterogeneous computing method
在搜索出SVSRB基础上，为进一步获取SVSRB近似解析表达式，以指导电网运行人员进行电压稳定性评估与控制，可参考文献[25]，再对SVSRB进行近似，得到SVSRB的近似解析表达式。
3 算例分析
针对所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法，本节首先以WECC3机9节点系统为例，验证本文所提方法的准确性；然后以欧洲13659节点系统为例，验证所提方法构建区域互联电力系统 SVSRB 的可行性；最后以case2737sop、case3120sp、case7092、case9241pegase等多个测试系统为例，验证所提方法的有效性。
本文算例计算硬件平台为 Dell PowerEdge C4130机架式服务器，CPU为Inter Xeon E5-2690，主频2.6GHz，内存256GB，GPU为Nvidia Tesla K80，总内存带宽480GB/s。软件平台为Matlab2018a，CUDA驱动版本为8.0。构建SVSRB时控制搜索方向的功率增长角Δq 为0.1rad，牛顿-拉夫逊法的收敛精度ε取10-2，最大迭代次数N取50，BICGSTAB迭代法的收敛阈值εb取10-2。
3.1 WECC3机9节点测试系统
本节首先以WECC3机9节点系统为例[32]，采用所提SVSRB并行计算方法，分别构建二维、三维有功功率注入空间内的SVSRB。
3.1.1 二维SVSRB构建
选择节点5、7的有功注入量为坐标轴，在以节点5、7形成的二维负荷有功注入空间内采用本文所提方法构建SVSRB，如图3所示。
图3 WECC-9系统的SVSRB
Fig.3 SVSRB of WECC-9 test system
由图3可知：首先，以基态为起始点，设置追踪第一个SNB点时所对应的功率增长方向d8=[0 0.326 0 0.170 0 0 0.180 0 0 0.200 0 0 0]T，d8对应的功率增长方向角q8=0.732 8rad，采用CPUGPU异构平台并行求解d8下的SNB点。具体过程如下：通过牛拉法对式（1）线性化处理并对牛顿迭代修正方程组式（2）进行降阶处理，得到式（4）～式（7）；然后，将系数矩阵B、线性方程组式（4）～式（7）右侧列向量、预处理器J、L和U，以及其他相关参数从CPU移植到GPU中，并利用GPU的并行加速技术，采用带两阶段预处理的BICGSTAB迭代法并行求解式（4）～式（7），得到中间变量α、β、ξ和η；再将求解的α、β、ξ和η传输到CPU中，利用各修正量与α、β、ξ和η的关系得到修正量Δx、Δr和Δl，并按式（3）对待求量x、r和l进行更新；重复上述过程，直到满足牛顿迭代收敛条件，则d8下的SNB点并行求解完毕，得到的初始SNB点如图3中SNB点8所示，其坐标为（2.388 0,2.149 2），最大负荷裕度l8=11.939 8(pu)。
然后，沿q 减小方向采用基于边界追踪的直接法并行搜索系统SNB点，详细过程如下：以SNB点8为起始点，令q7=q8-Δq，根据d8中节点5、7的有功增长分量ΔP5、ΔP7及q7得到新的功率增长方向d7=[0 0.326 0 0.170 0 0 0.159 1 0 0.217 0 0 0]T；以SNB点8的状态变量x8、负荷裕度l8及右特征向量r8作为直接法搜索d7下SNB点的初值，采用CPU-GPU异构平台实现SNB点的并行求解，得到d7对应的SNB点7，SNB点7的坐标为（2.661 2,1.951 8），最大负荷裕度l7=12.265 4(pu)。以此类推，将已求SNB点作为相邻待求解SNB点的起始点，其状态变量、负荷裕度和右特征向量作为直接法搜索相邻待求解SNB点的初值，并结合CPU-GPU异构平台实现SNB点的并行求解，依次得到SNB点6,5,4,…,0，其中SNB点0的功率增长方向角q0=-0.067 2rad，此时q0＜0，令q0=0，计算d0=[0 0.326 0 0.170 0 0 0 0 0.269 1 0 0]T下的SNB点0，并停止沿q 减小方向SNB点的搜索。
进一步，仍以SNB点8为起始点，类似上述过程，沿q增大的方向采用CPU-GPU异构的直接法依次搜索SNB点9,10,11,…,17，SNB点17的功率增长方向角为q17=1.632 8rad，此时有：q17＞π/2，令q17=π/2，计算d17=[0 0.326 0 0.170 0 0 0.269 1 0 0 0 0]T下的SNB点17，至此，系统在第一象限内沿q增大方向的SΝΒ点搜索结束。各SΝΒ点详细信息见表1。将图3搜索得到的SΝΒ点依次连接，即可获得由节点5、7形成的二维有功注入空间内的SVSRB。
表1 WECC-9系统中各SNB点信息
Tab.1 Detailed information of SNB of WECC-9 test system
图3进一步对比了本文所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法与CPF法、基于射线的直接法、基于边界追踪的直接法所构建的SVSRB。由图3可知：对于同一个功率增长方向，本文提出的SVSRB并行计算方法获得的SNB点和CPF、基于射线的直接法、基于边界追踪的直接法得到的SNB点近似重合，有效验证了所提的静态电压稳定域边界并行计算方法的准确性。
由于SNB点处潮流方程所对应的Jacobi矩阵发生奇异，Jacobi矩阵存在零特征值，故搜索得到的SNB点处Jacobi矩阵对应的最小特征值到0的距离可以作为衡量各种方法计算精度的标准。
图4分别对比了CPF、基于射线的直接法、基于边界追踪的直接法和本文所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法得到的SNB点处潮流Jacobi矩阵所对应的最小特征值。由图4中各SNB点处的最小特征值可以看出：本文所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法与基于射线的直接法、基于边界追踪的直接法计算得到的SNB点精度较CPF法搜索的SNB点精度高。图4中由CPF法所得的SNB点平均最小特征值为1.04×10-2，而本文所提方法与基于边界追踪的直接法、基于射线的直接法获得的最小特征值的平均值分别为1.21×10-3、1.20×10-3、2.01×10-3，该三种方法计算得到的SNB点精度大致相等，且比CPF的计算精度高出一个数量级，该对比结果表明：本文所提静态电压稳定域边界并行计算方法具有较高的计算精度。
图4 SVSRB的计算误差
Fig.4 Calculation errors of SVSRB
进一步，在图3搜索所得SVSRB基础上，参考文献[25]，对图3的SVSRB进行近似，近似结果如图5所示，对应的SVSRB局部近似边界H1、H2、H3、H4、H5、H6、H7、H8解析表达式如式（16）所示。
图5 WECC-9系统二维SVSRB的近似边界
Fig.5 Approximated 2-dimensional SVSRB of WECC-9 test case
3.1.2 计算效率分析
本节进一步对比了CPF法、基于射线的直接法、基于边界追踪的直接法与本文所提方法在上述场景下的计算效率。
表2给出了上述场景下不同方法构建二维SVSRB的计算时间。由表2可知，CPF在上述场景下的计算耗时远大于基于射线的直接法、基于边界追踪的直接法和本文所提出方法的计算耗时；基于边界追踪的直接法在上述场景中的计算耗时较基于射线的直接法、CPF法分别减少了24.8%和94.7%；本文所提方法的计算时间较基于边界追踪的直接法、基于射线的直接法耗时更长。出现上述结果的原因是：CPF法获得SNB点时涉及大量的预测-校正过程，该过程相对耗时，且随着SNB点数量的增多，CPF的耗时明显增加，因而基于CPF的SVSRB构建方法用时较长；而本文所提方法与基于射线的直接法、基于边界追踪的直接法均是对表征SNB点特性的非线性方程组进行求解，获得SNB点的过程相对简单，且所需的牛顿迭代次数远小于CPF法校正步的总迭代次数，因而计算时间较CPF更短；但因测试系统较小，CPU与GPU数据交互的耗时远大于GPU并行求解带来的加速性，因而本文所提方法较基于边界追踪的直接法、基于射线的直接法耗时更久。
表2 不同方法的计算时间对比
Tab.2 Computational times compared with different methods
由于直接法搜索电力系统SNB点的核心是利用牛拉法对式（1）进行求解，因此，构建SVSRB的计算耗时主要依赖于获取SNB点时的牛拉法总迭代次数。表3对比了本文所提方法与基于射线的直接法、基于边界追踪的直接法搜索SVSRB的迭代次数。
表3 不同方法的迭代次数
Tab.3 Total times of iterations compared with different methods
由表3的对比结果可知：基于边界追踪的直接法搜索SVSRB时所需的总迭代次数较少，因而搜索SVSRB时计算耗时较少；基于射线的直接法搜索SVSRB所需总迭代次数较多，计算耗时较基于边界追踪的直接法更长；而本文所提方法虽有较少的迭代次数，但构建SVSRB所需的计算时间比基于射线的直接法和基于边界追踪的直接法更久。其主要原因是：基于射线的直接法始终以基态为起始点，各SNB点的搜索相对独立，且需不断利用负荷参数二阶导数预估不同功率增长方向下的SNB点的初值，因此获得SNB点收敛解所需的牛顿迭代次数较多，收敛耗时较长；而基于边界追踪的直接法将已求SNB点作为相邻待求解SNB点的初始点，因初值更接近待求解SNB点的真实值，故获取SNB点时所需的迭代次数更少，构建SVSRB的计算耗时更短；本文所提方法是在基于边界追踪的直接法基础上利用CPU-GPU异构平台进行的SVSRB并行计算，因而牛顿总迭代次数较少，但本文采用CPU-GPU异构平台求解SNB点时，由于测试系统规模较小，CPU与GPU间数据交互耗时远大于GPU并行加速的计算耗时，因而构建SVSRB计算时间相对较长。
3.1.3 功率增长角Δq 的分析
由于SVSRB上获得各SNB点的功率增长方向受功率增长角Δθ的影响，Δθ过大或过小均会影响SVSRB的构建精度及构建速度。因此，本节进一步分析了Δθ对构建SVSRB的具体影响。
表4为不同Δθ时，在上述坐标空间中构建二维SVSRB所需计算耗时与牛拉法总迭代次数的对比。由表4可知：当Δθ＝0.05rad时，构建SVSRB所需的SNB点数量较多，故构建SVSRB时牛顿迭代的总迭代次数较多，计算耗时较长，但因相邻SNB点距离较近，因此平均迭代次数相对较少；当Δθ＝0.25rad时，构建SVSRB所需SNB点数量较少，故牛顿迭代法的总迭代次数相对较少，搜索SVSRB的计算耗时较短，但因相邻SNB点距离相对较远，因此牛拉法迭代时所需的平均迭代次数相对较多。
表4 不同Δθ 取值下的SVSRB计算效率对比
Tab.4 Computational efficiency comparisons of SVSRB with different Δθ
随着Δθ的不断增大，构建SVSRB所需SNB点数量不断减少，计算时间也相对减少，但SVSRB的构建精度与所搜索的SNB点数量密切相关，搜索的SNB点数量越多，构建的SVSRB越准确。因而，Δθ的增大在提高SVSRB构建效率的同时，也降低了SVSRB构建精度。图6给出了在节点5、7形成的二维负荷有功注入空间内，由不同Δθ下构建的SVSR局部边界。由图6可知，Δθ取0.1rad时，SVSR区域由A、B和C三部分组成，Δθ取0.25rad时，SVSR区域由A组成。显然，由区域A、B和C共同组成的区域更加接近真实的SVSR。Δθ取0.25rad时得到的SVSR区域较为保守，当系统运行在区域B或C时，会造成系统误报警，因此，Δθ的取值需在SVSRB的搜索精度和搜索效率之间达到平衡。对比表4和图6的计算结果可知：本文中Δθ取值为0.1rad是合理、可行的。
图6 不同Δθ构建的SVSRB
Fig.6 Constructed SVSRB compared with different Δθ
3.1.4 迭代步长的分析
本节进一步详细分析了牛拉法收敛精度ε和BICGSTAB迭代法收敛阈值εb的选取过程。
1）牛拉法收敛精度ε的选取
采用边界追踪法搜索SVSRB时，SVSRB上相邻SNB点处的系统运行信息变化较小，若牛拉法收敛精度过小，则在提高SNB点计算精度的同时，也会增加牛顿迭代次数，使得SNB点的计算耗时明显增加；而牛拉法收敛精度过大，则在减少牛顿迭代次数、提高SNB点计算效率的同时，可能无法保证SNB点的计算精度。
表5对比了牛拉法不同收敛精度对上述场景中SNB点8的计算精度和计算效率的影响。由表5可知，当牛拉法收敛精度为10-2时，直接法即可获得准确的负荷裕度，且相较于其他收敛精度，具有较少的迭代次数和较高的计算速度。
表5 不同ε下SNB点8的计算效率对比
Tab.5 Computational efficiency comparisons of 8th SNB with different ε
表6进一步对比了牛拉法收敛精度ε分别为10-2和10-8时SNB点8处系统各节点的电压幅值，并以收敛精度10-8为基准，计算了各节点电压幅值相对误差。由表6可知，ε=10-2时的各节点电压幅值相对误差均低于0.05%，进一步验证了所设牛拉法收敛精度的准确性。
表6 不同ε下SNB点8处各节点电压幅值对比
Tab.6 Comparison of voltage magnitude of 8th SNB with different ε
综上所述，本文设置的牛拉法收敛精度10-2可以很好地保证SNB点精度，且具有较高的计算效率，因此，本文中的牛拉法迭代步长选为10-2。
2）BICGSTAB迭代法收敛阈值εb的选取
BICGSTAB迭代法是以当前残差与预设阈值进行比较来判断何时停止迭代，故对阈值的修改即可实现对求解精度的控制[24]。本文中，BICGSTAB迭代法主要用于求解牛顿迭代内部的四个低维线性方程组式（4）～式（7），因此，其收敛阈值主要影响式（4）～式（7）的求解精度。然而，对式（4）～式（7）的求解并不是为了得到其精确解，而是为了使牛顿迭代可以向非线性方程组式（1）的真实解逼近，因此，BICGSTAB迭代法的收敛精度不需太高，能保证收敛即可。
通常，采用迭代求解法求解牛拉法内部的线性方程组时，迭代求解法的收敛精度应不低于牛拉法的收敛精度[31]，结合本文选取的牛拉法收敛精度，故BICGSTAB迭代法的收敛精度应不低于10-2。而将BICGSTAB迭代法的求解精度调低可以有效降低内部线性求解所需的迭代次数，提高计算效率，故本文中BICGSTAB迭代法的收敛阈值取为10-2。
3.1.5 三维SVSRB构建
本节在二维SVSRB构建的基础上，进一步借鉴文献[23]的方法，以节点5、7、9的有功注入为坐标轴，构建了如图7所示的三维SVSRB。由图7的三维SVSRB搜索结果可知：本文所提基于CPUGPU异构的静态电压稳定域边界并行计算方法在二维、三维有功注入空间均可以实现SVSRB的准确构建。
图7 WECC3机9节点系统的三维SVSRB
Fig.7 3-dimensional SVSRB of WECC 9 test system
进一步，借鉴文献[25]的安全域边界近似方法，对图7所示的三维SVSRB近似，结果如图8所示，相应的边界近似解析表达式详见附表1。
附表1 WECC-9节点系统的三维SVSRB近似边界解析表达式系数
App.Tab.1 Analytical expressions for the 3-dimensional SVSRB approximation hyperplanes of WECC-9 test system
图8 WECC-9系统的三维SVSRB近似边界
Fig.8 Approximated 3-dimensional SVSRB of WECC-9 test system
3.2 Case13659pegase测试系统
WECC3机9节点系统算例验证了所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法的准确性，本节进一步以欧洲电网13 659节点系统为例[33]，采用所提的静态电压稳定域边界并行计算方法分别构建case13659pegase系统的二维和三维SVSRB，以进一步验证所提方法构建大规模区域互联电力系统SVSRB的可行性。
分别选择节点4、5和节点4、6作为关键节点，采用本文所提的方法分别在节点4、5和4、6的二维有功注入空间构建SVSRB，结果如图9所示。进一步，以节点4、5、6的有功注入量为坐标轴构建三维SVSRB，结果如图10所示。
图9 case13659pegase系统的二维SVSRB
Fig.9 2-dimensional SVSRB of European 13659PEGASE test system
图9和图10结果表明：本文所提基于CPUGPU异构的静态电压稳定域边界并行计算方法可以用于高维电力系统SVSRB的构建，验证了本文所提方法的可行性。
进一步，对图9和图10所得的二维、三维SVSRB进行近似，结果如图11和图12所示。其中，图9a和图11a中，以节点4、5有功注入为坐标轴的二维SVSRB近似解析表达式为
图10 case13659pegase系统的三维SVSRB
Fig.10 3-dimensional SVSRB of European 13659PEGASE test system
图11 case13659pegase系统的二维SVSRB近似边界
Fig.11 Approximated 2-dimensional SVSRB of European 13659PEGASE test system
图12 case13659pegase系统的三维SVSRB近似边界
Fig.12 Approximated 3-dimensional SVSRB of European 13659PEGASE test system
图9b和图11b中，以节点4、6有功注入为坐标轴的二维SVSRB近似解析表达式为
图10和图12中，以节点4、5和6有功注入为坐标轴的SVSRB近似解析表达式详见附表2。
附表2 欧洲电网13659节点系统三维SVSRB近似边界解析表达式系数
App.Tab.2 Analytical expressions for the 3-dimensional SVSRB approximation hyperplanes of European 13659PEGASE test system
3.3 SVSRB并行计算效率分析
为验证所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法构建区域互联电力系统SVSRB的有效性，本节分别以case2737sop、case3120sp、case7092、case9241pegase、case11624和case13659pegase测试系统为例，对基于CPUGPU异构平台的BICGSTAB迭代法与Jacobi+ILU(0)预处理相结合的计算方法、基于CPU平台的BICGSTAB迭代法与Jacobi+ILU(0)预处理相结合的计算方法、基于CPU平台的LU分解算法以及基于射线的直接法搜索二维、三维SVSRB的计算耗时进行对比。上述四种算法搜索的二维和三维SVSRB计算耗时结果见表7和表8。
表7和表8中“—”表示测试系统的SVSRB构建失败。表中基于射线的直接法构建电力系统SVSRB失败的主要原因为：利用负荷参数二阶导数预估SNB点初值存在给定初值不稳定、不精确的现象，从而导致直接法计算部分测试系统SVSRB上的SNB点时结果不收敛。而基于CPU平台的LU分解算法构建电力系统SVSRB失败的主要原因是由于LU分解算法受功率增长方向的影响，不合理的功率增长方向会使LU分解算法求解线性方程组式（4）～式（7）出现错误解，从而导致牛顿迭代不收敛，SNB点求解失败。对比表7和表8结果可知：基于射线的直接法虽理论上可行，但实际应用时很难构建出系统的SVSRB；基于CPU平台的LU分解算法虽然可以快速搜索SVSRB，但该算法稳定性较差；而对式（4）～式（7）采用Jacobi+ILU(0)两阶段预处理的BICGSTAB迭代法与牛顿迭代法相结合的方式使得本文所提基于边界追踪的直接法对于上述测试系统无论在CPU平台还是CPU-GPU异构平台均可构建出系统完整的SVSRB，有效验证了本文所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法适用性更强、鲁棒性更高。其主要原因是：直接法对初值非常敏感，基于射线的直接法因以基态为起始点，基态与SNB点距离较远，直接法牛顿迭代时给定的初值不够合理，因而牛顿迭代时SNB点收敛较慢甚至不收敛，不能构建完整的SVSRB；而高度优化的LU分解算法可快速求解式（4）～式（7）并结合边界追踪法可实现SVSRB的高效构建，但基于CPU平台的LU分解算法对功率增长方向较为敏感，不当的功率增长方向会导致SNB点求解失败，无法构建完整的SVSRB；而本文所提并行计算方法因给定的牛顿迭代初值非常接近相邻待求解SNB点的真实值，SNB点初值精度较高，因此易获得SNB点收敛解，便于构建完整的SVSRB。
表7 不同测试系统的二维SVSRB计算时间对比
Tab.7 Computational times of 2-dimensional SVSRB for different bulk power grids
表8 不同测试系统的三维SVSRB计算时间对比
Tab.8 Computational times of 3-dimensional SVSRB for different bulk power grids
进一步分析表7和表8中CPU平台和CPUGPU异构平台下采用基于边界追踪的直接法搜索SVSRB的结果表明：对于case2737sop和case3120sp测试系统，本文所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法构建SVSRB的计算耗时相较基于CPU平台的基于边界追踪的直接法更长。其主要原因是：GPU作为CPU的外部设备，二者间的数据移植比较耗时[34]，利用GPU的并行技术加速对case2737sop和case3120sp测试系统线性方程组求解的同时，CPU和GPU之间的数据传输也会耗费一定时间，且对于 case2737sop和case3120sp测试系统而言，CPU和GPU间的数据交互耗时要多于GPU并行求解带来的加速性。当测试系统规模大于case3120sp时，采用所提基于CPUGPU异构的静态电压稳定域边界并行计算方法构建SVSRB的效率优势逐渐突显，对case13659pegase测试系统的二维SVSRB，所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法构建SVSRB的加速比是在CPU平台中构建SVSRB的2.25倍。
表9以CPU平台构建SVSRB的计算耗时为基准，给出了采用所提基于CPU-GPU异构平台的SVSRB并行计算方法在搜索二维和三维SVSRB时的相对计算效率。由于采用CPU-GPU异构平台的SVSRB搜索方法在case2737sop和case3120sp系统的计算耗时大于在CPU平台中的构建耗时，因而对应的表9中的计算效率提升量记为0。表9的计算耗时对比结果表明：对于大型区域互联电力系统，采用本文所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法可进一步加快SVSRB的构建效率，降低了构建SVSRB所需的计算耗时，验证了所提方法在构建大规模区域互联电力系统静态电压稳定域中的有效性。
表9 不同测试系统的SVSRB并行计算效率对比
Tab.9 Comparisons of computational efficiency for SVSRB in bulk power grids
上述对比结果表明：本文所提基于CPU-GPU异构的静态电压稳定域边界并行计算方法适用于区域互联电力系统SVSRB的构建，且在保证高精度的前提下具有较高的计算效率。
4 结论
本文提出一种基于CPU-GPU异构的电力系统静态电压稳定域边界并行计算方法，通过WECC3机9节点测试系统、波兰电网2737节点测试系统和3120节点测试系统、欧洲电网7092节点测试系统、9241节点测试系统、11624节点测试系统和13659节点测试系统算例对本文所提方法的准确性和有效性进行了分析验证，相关结论如下：
1）所提方法将直接法与边界追踪法相结合，成功克服了直接法初值选择的不足，且与CPU-GPU异构平台有效结合，解决了直接法计算量大、计算复杂度高的瓶颈。
2）所提方法可实现二维以及高维有功注入空间内SVSRB的准确构建,且采用CPU-GPU异构平台有效降低了SVSRB的构建耗时，提高了SVSRB的计算效率。
3）所提方法因给定初值稳定，且采用Jacobi+ILU(0)两阶段预处理的BICGSTAB迭代法求解线性方程组，具有良好的算法稳定性和适用性。
所提基于CPU-GPU异构的电力系统静态电压稳定域边界并行计算方法虽可实现由SNB点构成的SVSRB的快速构建，但对于实际电力系统，LIB点同样也可诱发系统电压失稳，因此，如何实现含SNB点和LIB点的SVSRB并行计算将是本文下一步的研究重点。
附 录
（续）
[1] 朱泽安,周修宁,王旭,等.基于稳暂态联合仿真模拟的区域多可再生能源系统评估决策[J].电工技术学报,2020,35(13):2780-2791.
Zhu Zean,Zhou Xiuning,Wang Xu,et al.Evaluation and decision-making of regional multi-renewable energy system based on steady-transient integrated simulation[J].Transactions of China Electrotechnical Society,2020,35(13):2780-2791.
[2] 张炎,丁明,韩平平,等.直流闭锁后风电送端系统暂态稳定及控制策略研究[J].电工技术学报,2020,35(17):3714-3726.
Zhang Yan,Ding Ming,Han Pingping,et al.Study on the transient stability and control schemes of the sending-end system involving wind power after UHVDC block[J]. Transactions of China Electrotechnical Society,2020,35(17):3714-3726.
[3] Hussien M G,Liu Yi,Xu Wei.Robust position observer for sensorless direct voltage control of standalone ship shaft brushless doubly-fed induction generators[J].CES Transactions on Electrical Machines and Systems,2019,3(4):363-376.
[4] 卢春宏,章玮.无电解电容永磁同步电机驱动系统母线电压解析分析[J].电机与控制学报,2019,23(10):15-22.
Lu Chunhong,Zhang Wei.Analytical analysis on PMSM control systems without electrolytic capacitor[J].Electric Machines and Control,2019,23(10):15-22.
[5] 张艳霞,尹佳鑫,蒙高鹏,等.一种基于广域测量信息的在线同调分群方法[J].电机与控制学报,2019,23(5):10-17.
Zhang Yanxia,Yin Jiaxin,Meng Gaopeng,et al.Method of online recognition of coherent generators based on wide area information [J].Electric Machines and Control,2019,23(5):10-17.
[6] 赖一峰.直流配电网的电压协同控制及稳定运行研究[J].电气技术,2019,20(7):42-47.
Lai Yifeng.Research on voltage synergy control and stable operation of DC distribution network[J].Electrical Engineering,2019,20(7):42-47.
[7] Jiang Tao,Bai Lingquan,Jia Hongjie,et al.Identification of voltage stability critical injection region in bulk power systems based on the relative gain of voltage coupling[J].IET Generation,Transmission &Distribution,2016,10(7):1495-1503.
[8] 乐健,廖小兵,李奔,等.基于扩展仿射模型的不确定性静态电压稳定性全局灵敏度分析[J].电工技术学报,2020,doi:10.19595/j.cnki.1000-6753.tces.200850.
Le Jian,Liao Xiaobing,Li Ben,et al.Global Sensitivity Analysis of Uncertain Static Voltage Stability Based on Extended Affine Model[J].Transactions of China Electrotechnical Society,2020,doi:10.19595/j.cnki.1000-6753.tces.200850.
[9] 陈中,严俊,朱政光,等.多馈入交直流混联系统解耦安全域的刻画及应用[J].电力系统自动化,2020,44(22):37-44.
Chen Zhong,Yan Jun,Zhu Zhengguang,et al.Characterization and application of decoupling security region in multi-infeed AC/DC hybrid system[J].Automation of Electric Power Systems,2020,44(22):37-44.
[10] 肖峻,肖居承,张黎元,等.配电网的严格与非严格安全边界[J].电工技术学报,2019,34(12):2637-2648.
Xiao Jun,Xiao Jucheng,Zhang Liyuan,et al.Strict and non-strict security boundary of distribution network[J].Transactions of China Electrotechnical Society,2019,34(12):2637-2648.
[11] 杨明,程凤璐,韩学山,等.电力系统实时调度的有效静态安全域法[J].中国电机工程学报,2015,35(6):1353-1362.
Yang Ming,Cheng Fenglu,Han Xueshan,et al.Realtime dispatch based on effective steady-state security regions of power system[J].Proceedings of the CSEE,2015,35(6):1353-1362.
[12] 姜涛,谭洪强,李雪,等.电力系统热安全域边界快速搜索的优化模型[J].中国电机工程学报,2019,39(22):6533-6546.
Jiang Tao,Tan Hongqiang,Li Xue,et al.A general optimization model for exploring static thermal security region boundary in bulk power systems[J].Proceedings of the CSEE,2019,39(22):6533-6546.
[13] 肖峻,曹严,张宝强,等.配电网安全域的全维直接观测[J].电工技术学报,2020,35(19):4171-4182.
Xiao Jun,Cao Yan,Zhang Baoqiang,et al.Fulldimensional direct observation of distribution system security region [J]. Transactions of China Electrotechnical Society,2020,35(19):4171-4182.
[14] 苗伟威,贾宏杰,董泽寅.基于有功负荷注入空间静态电压稳定域的最小切负荷算法[J].中国电机工程学报,2012,32(16):44-52.
Miao Weiwei,Jia Hongjie,Dong Zeyin.A minimum load shedding calculation algorithm based on static voltage stability region in active load power injection space[J].Proceedings of the CSEE,2012,32(16):44-52.
[15] 王成山,许晓菲,余贻鑫,等.电力系统电压稳定域的局部可视化描述及其应用[J].中国电机工程学报,2004,24(3):1-5.
Wang Chengshan,Xu Xiaofei,Yu Yixin,et al.Visualization of part of the static voltage stability region in power systems and its application[J].Proceedings of the CSEE,2004,24(3):1-5.
[16] 韩琪,余贻鑫,李慧玲,等.电力系统注入空间静态电压稳定域边界的实用表达式[J].中国电机工程学报,2005,25(5):8-14.
Han Qi,Yu Yixin,Li Huiling,et al.A practical boundary expression of static voltage stability region in injection space of power systems [J].Proceedings of the CSEE,2005,25(5):8-14.
[17] Adusumilli B S,Kumar B K.Modified affine arithmetic based continuation power flow analysis for voltage stability assessment under uncertainty[J].IET Generation,Transmission &Distribution,2018,12(18):4225-4232.
[18] Hao Sheng,Chiang H D.CDFLOW-a practical tool for tracing stationary behaviors of general distribution networks[J].IEEE Transactions on Power Systems,2014,29(3):1365-1371.
[19] 胡林麟.基于优化的电力系统潮流可行域边界计算与分析[D].天津:天津大学,2012.
[20] 江伟,王成山,余贻鑫,等.直接计算静态电压稳定临界点的新方法[J].中国电机工程学报,2006,26(10):1-6.
Jiang Wei,Wang Chenshan,Yu Yixin,et al.A new method for direct calculating the critical point of static voltage stability[J].Proceeding of the CSEE,2006,26(10):1-6.
[21] Jiang Tong,Wan Kaiyao,Feng Zhuocheng.Boundaryderivative direct method for computing saddle node bifurcation points in voltage stability analysis[J].International Journal of Electrical Power and Energy Systems,2019,112:199-208.
[22] 姜涛,张明宇,李雪,等.静态电压稳定域局部边界的快速搜索新方法[J].中国电机工程学报,2018,38(14):4126-4137,4318.
Jiang Tao,Zhang Mingyu,Li Xue,et al.A novel algorithm to explore static voltage stability region boundary in power systems[J].Proceedings of the CSEE,2018,38(14):4126-4137,4318.
[23] 姜涛,张明宇,崔晓丹,等.电力系统静态电压稳定域边界快速搜索的优化模型[J].电工技术学报,2018(1):4167-4179.
Jiang Tao,Zhang Mingyu,Cui Xiaodan,et al.A novel optimization model to explore static voltage stability region boundary in bulk power systems[J].Transactions of China Electrotechnical Society,2018(1):4167-4179.
[24] 李雪,刘烨,姜涛,等.电力系统负荷裕度的并行计算方法研究[J].中国电机工程学报,2019,39(17):5105-5117,5291.
Li Xue,Liu Ye,Jiang Tao,et al.Parallel calculation of Load margin in bulk power grid using CPU-GPU architecture[J].Proceedings of the CSEE,2019,39(17):5105-5117,5291.
[25] 姜涛,李晓辉,李雪,等.电力系统安全域边界通用搜索模型与近似方法[J].中国电机工程学报,2020,40(14):4411-4429.
Jiang Tao,Li Xiaohui,Li Xue,et al.General optimization model and piecewise approach for approximating security region boundary in bulk power systems[J].Proceedings of the CSEE,2020,40(14):4411-4429.
[26] Yang Tiankai,Yu Yixin.Static voltage security region-based coordinated voltage control in smart distribution grids[J].IEEE Transactions on Smart Grid,2018,9(6):5494-5502.
[27] 余贻鑫.电力系统安全域方法研究述评[J].天津大学学报,2008,41(6):635-646.
Yu Yixin.Review of study on methodology of security regions of power system[J].Journal of Tianjin University,2008,41(6):635-646.
[28] 陈昌,姜彤,万凯遥,等.直接计算静态电压稳定裕度的改进崩溃点法[J].电力自动化设备,2020,40(11):150-155.
Chen Chang,Jiang Tong,Wan Kaiyao,et al.Improved point of collapse method for direct calculation of static voltage stability margin[J].Electric Power Automation Equipment,2020,40(11):150-155.
[29] 宋晓喆,魏国,李雪,等.基于预处理BICGSTAB法的电力系统潮流并行计算方法[J].电力系统保护与控制,2020,48(20):18-28.
Song Xiaozhe,Wei Guo,Li Xue,et al.Parallel power flow computing in power grids based on a preconditioned BICGSTAB method[J]. Power Systems Protection and Control,2020,48(20):18-28.
[30] Li Xue,Li Fangxing,Yuan Haoyu,et al.GPU-based fast decoupled power flow with preconditioned iterative solver and inexact Newton method[J].IEEE Transactions on Power Systems,2017,32(4):2695-2703.
[31] Li Xue,Li Fangxing.GPU-based power flow analysis with Chebyshev preconditioner and conjugate gradient method[J].Electric Power Systems Research,2014,116(11):87-93.
[32] 贾宏杰,姜涛,姜懿郎,等.电力系统时滞稳定域临界点的快速搜索新方法[J].电力系统自动化,2013,37(6):30-36.
Jia Hongjie,Jiang Tao,Jiang Yilang,et al.A method to search critical point of time-delay stability region inpower systems[J].Automation of Electric Power Systems,2013,37(6):30-36.
[33] 姜涛,李晓辉,李雪,等.电力系统静态电压稳定域边界近似的空间切向量法[J].中国电机工程学报,2020,40(12):3729-3743.
Jiang Tao,Li Xiaohui,Li Xue,et al.Approximating voltage stability region boundary in power systems using tangent vector[J].Proceedings of the CSEE,2020,40(12):3729-3743.
[34] 叶盛.基于CPU-GPU-FPGA的异构计算系统及任务调度算法研究[D].西安:西安电子科技大学,2019.