基于互信息贝叶斯网络的配电网拓扑鲁棒辨识算法

任鹏哲，刘友波，刘挺坚，何培东，张扬帆，邓舒予

（1. 四川大学电气工程学院，四川省成都市610065；2. 国网四川省电力公司计量中心，四川省成都市610045）

0 引言

近年来，分布式发电和新型负荷接入配电网的规模不断扩大，其波动性、随机性和不确定性对配电网运行方式产生了较大影响。为保证安全、稳定和经济运行，电网的运行方式需要面临较为频繁的变化［1］。同时，为了实现配电网可观、可测和可控性，需要研究实时跟踪配电网运行状态的技术，而作为配电网状态估计的基础，电网拓扑的识别显得尤为重要。

电网拓扑识别的传统方法有树搜索法［2］、矩阵法［3-4］、图划分法［5］、转移潮流法［6］、图论法［7-8］、不确定性推理法［9］等。上述方法依托数据采集与监控（supervisory control and data acquisition，SCADA）系统，根据采集的线路潮流、节点电压或设备开关状态进行拓扑识别。SCADA 系统数据稳定性差，遥信数据经常误报和漏报，遥测数据较为单一，容错能力不足。此外，SCADA 系统数据缺乏精准时标，采集数据有可能来自不同时刻，导致信息时标错位，给识别带来干扰。微型同步相量测量装置（micro phasor measurement unit，μPMU）［10］的出现与应用为解决上述问题提供了新思路。μPMU 数据采样率和精准度高，采样数据带有时标且能保持时序同步，可靠性和容错性好。已有部分学者基于μPMU 提出拓扑辨识方法：文献［11-12］通过列举可能拓扑，将测量的电压相角与各拓扑下的估计值进行内积匹配，实现拓扑辨识；文献［13］基于历史拓扑形成拓扑库，根据拓扑变化前后节点电压变化趋势识别拓扑；文献［14］提出一种μPMU 投票算法，每个μPMU 根据已有数据估计各拓扑置信度，以综合置信度最高的拓扑作为识别结果。上述基于μPMU 的识别算法也存在2 个问题：①通过逐个比较潮流估计值和测量值确定拓扑的算法，在规模较大或含有较多可行拓扑的网络中难以保证时效性；②基于实时潮流计算的算法鲁棒性较弱，在μPMU 故障、负荷和光伏等关键数据缺失时容易产生错误结果。

贝叶斯网络（Bayesian network，BN）［15］是一种可以直观反映变量因果关系的推理模型。该模型鲁棒性强，在输入信息异常或缺失时仍能保证良好的推理效果。但BN 在处理连续型变量时，需要人为设置变量区间数目［16］，难以适应较多连续性变量的情况。 为此，引入最大互信息系数（maximal information coefficient，MIC）方法的网格划分［17］衡量BN 中连续型节点的区间数目，有效减少了网络训练复杂度，保证了拟合效果。基于拉丁超立方抽样（Latin hypercube sampling，LHS）［18］生成光伏-负荷数据，实现场景在样本空间的均匀分布，光伏和负荷之间无明显相关性。在此前提下，为其设置相同的区间数目，进一步简化了BN 的训练过程，并有效提高了拓扑识别的泛化能力。算例结果表明，所提方法在时效性及鲁棒性方面表现出较强优势。

1 互信息BN

1.1 MIC 方法

MIC 方法是一种基于互信息和网格划分的关联评价方法，能够捕获变量间不同类型的联系并衡量其程度。计算变量x 和y 的MIC 步骤如下。

步骤1：在给定的网格分辨率下，对x 和y 组成的二维空间按不同网格划分方案进行划分，并求出不同方案的最大互信息值I，表达式为：

式中：P（x）和P（y）分别为变量x 和y 的边缘分布概率；P（x，y）为x 和y 的联合分布概率；X 和Y 分别为x 和y 的集合。

步骤2：将求得的最大互信息值进行归一化处理，归一化后的取值范围为［0，1］。

步骤3：改变网格分辨率并重新执行步骤1 和2，以不同分辨率下互信息最大值作为MIC 的值IMIC，表达式为：

式中：a 和b 分别为x 和y 方向上的网格划分个数；B为网格数上限。

经上述处理后，可以获得x 和y 的最大互信息值及对应的网格划分情况。本文就是利用该网格划分来衡量BN 连续型节点的区间划分。

1.2 BN 模型

BN 是一种将图论与概率学结合在一起的概率推理模型，在结构上是一个有向无环图。附录A 图A1 展示了BN 的模型，可用G=（V，E）来表示。其中，V 为节点集，V={vi}，vi为第i 个节点，有因果关系的2 个节点可以用有向线段连接，E 为有向线段集。

对于BN 中任意随机变量，其联合概率可由各自局部条件概率分布相乘得出，即

考虑条件独立性假设，即给定某节点的父节点集，该节点独立于它的所有非后代节点，式（3）可以表达为：

式中：k 为节点数；π（vi）为vi的父节点集。

在得到网络结构后，各节点的关系随之确定。利用BN 进行推理时，通过输入已知证据，即可更新各节点的后验概率，得到所要推知事件的概率。

1.3 BN 学习

BN 学习分为结构学习和参数学习。结构学习是指确定BN 的拓扑结构；参数学习是指确定网络各节点的条件概率和联合分布概率。

1.3.1 BN 结构学习

常用的BN 结构学习算法［19］有K2 学习算法和爬山搜索算法。这2 种算法均基于评分搜索函数，通过寻找评分值最高的网络结构来获得最佳的网络结构，其中，评分搜索函数S 代表了网络结构与数据集的拟合程度，其表达式为：

式中：n 为网络节点数；qi为节点xi父节点取值组合的数量；ri为节点xi取值数目；mijk为样本集中节点xi为第k 个取值，其父节点为第j 种取值组合的样本数量；mij为mijk在k ∈[1，ri]内的和。

与爬山搜索算法相比，K2 学习算法可以通过预先给出的节点顺序。在寻找BN 最优拓扑的过程中融入先验信息，但在未知因果关系的情况下通过主观定义节点顺序会降低结构学习的准确性，故选择爬山算法进行BN 的结构学习。首先，以随机生成的方式产生大量的初始BN 结构。然后，以评分搜索函数为依据，搜索评分最高的网络结构作为最优结构。

1.3.2 参数学习

参数学习［20］的主要方法有最大似然估计法、贝叶斯估计法和最大后验估计法。由于本文的训练样本为较完备的数据集，故采用完全基于数据、不需考虑网络参数先验分布的最大似然估计法，最大似然估计的表达式为：

式中：D 为样本集；θ 为参数；P（D|θ）为θ 的条件概率函数；L（θ|D）为θ 的似然度函数。

令θ 在定义域内变动，寻找使L（θ|D）达到最大的θ 值称为参数的最大似然估计，即

式中：θ*为令L（θ|D）取最大值时的θ 值。

2 基于BN 的配电网拓扑识别

目前，配电网的拓扑分析大多是在实时量测数据的基础上，通过电力系统状态估计进行的。数据的测量错误或缺失将会直接影响拓扑识别的正确性。此外，基于实时潮流计算模型的拓扑识别算法需要依次对可能的拓扑进行匹配估计，其效率不高。为此，提出一种基于BN 的配电网拓扑识别方法。其核心思想是通过LHS 生成在样本空间内分布均匀的光伏-负荷场景。通过离线的潮流计算模型获取各可行拓扑在不同光伏-负荷场景下的状态量（装有μPMU 节点的电压信息），利用BN 拟合拓扑光伏-负荷状况与电压数据的映射关系，形成基于BN 的拓扑识别器。当需要估计当前拓扑时，首先，通过μPMU 获取在时间上同步的测量点电压信息。然后，由SCADA 系统或预测系统［21］采集实时的光伏-负荷状况。最后，将上述采集到的数据作为BN 的先验信息，并设置拓扑节点为输出节点即可获得各个拓扑的概率，按照式（8）寻找可能概率最高的拓扑作为识别结果。

式中：Pi为BN 推理后给出的第i 个拓扑的可能概率。

具体的算法流程如图1 所示。

2.1 BN 训练样本的生成

2.1.1 配电网拓扑结构

本文的主要研究对象为辐射状结构的中压配电网，其结构如附录A 图A2 所示。实际中，辐射状配电网的运行拓扑存在着相关约束和规则，使得其数量保持在一定范围内。在该前提下，利用拓扑库的形式对历史拓扑进行存储，可有效提高拓扑识别效率和效果。按照文献［22］所提约束，共生成N 个可行拓扑进行算法的验证。

2.1.2 LHS 和样本的生成

LHS 是一种分层抽样技术，其采样均匀度高，能保证采样结果在样本空间的均匀分布。附录A图A3 展示了利用LHS 和简单随机抽样（simple random sampling，SRS）在二维空间生成样本的分布情况，可以看出LHS 的样本分布更加均匀，基于LHS 生成样本可有效提高拓扑识别器的泛化能力。

利用LHS 生成BN 训练样本的步骤如下。

步骤1：利用LHS 生成M 个负荷和光伏出力在一定范围内波动的光伏-负荷场景。这M 个场景与拓扑库中的N 个可行拓扑进行组合，共形成M×N个样本。

步骤2：根据样本设置光伏、负荷和拓扑状态等变量，经潮流计算后记录设置μPMU 的节点电压幅值和相位，形成含有电压信息的样本集，作为贝叶斯拓扑识别模型的训练集和测试集。

步骤3：考虑到实际中采集的数据包含随机噪声等干扰，为了提高贝叶斯拓扑识别模型的抗噪声能力，在训练集样本中加入带有误差的噪声样本。误差e 服从正态分布，其概率密度函数f (e)的表达式为：

式中：μ 为e 的平均数；σ 为e 的标准差。

2.2 贝叶斯拓扑推理网络的搭建

根据2.1.2 节生成的训练样本，贝叶斯拓扑推理网络的构建步骤如下。

步骤1：对拓扑库中的拓扑结构进行编号，按各支路的开断状态形成二进制编码，以二进制编码的大小顺序设置各拓扑的编号。

步骤2：样本数据预处理。样本中光伏、负荷和电压等变量属于连续型数据，若不经处理直接进行BN 训练，一方面会导致网络节点的属性过多影响训练效率；另一方面会造成训练过拟合，影响网络的鲁棒性，故需要对上述数据进行离散处理。目前，对BN 连续型节点的处理方法大多通过人为划分，并根据实际效果迭代调整，主观性强。连续型节点较多的情况是一个维度极高的组合问题，难以适应本文基于多维连续数据进行拓扑识别的场景。为此，引入MIC 方法确定连续型节点区间划分数目。该方法仅需调整光伏-负荷节点的区间数量，电压信息变量的区间划分可通过寻找具有MIC 方法的网格划分得出，在简化训练的同时保证了BN 的拟合效果。

由于MIC 方法是通过网格划分量度各变量的相关性，用MIC 网格划分连续型数据的区间可以尽可能保留变量间相关性。获得离散区间数量后，由ChiMerge 法［16］对各变量进行离散化处理。ChiMerge 法基于卡方分析，将相邻区间中卡方值最小的区域进行合并，直至区间数量收敛至设定值，卡方值χ2的表达式为：

式中：m 为每次比较的区间数目；n1为样本簇个数；Aij为第i 类区间内第j 簇样本的个数；Eij为Aij的数学期望。

步骤3：经前2 步获取多维变量（负荷、光伏数据、电压信息和拓扑结构）训练样本后，采用1.3.1 节的爬山搜索算法获得BN 的拓扑结构，用最大似然估计法学习BN 的参数，从而建立起能够刻画随机变量相关性的BN 模型。

3 算例分析

3.1 样本说明

在改造的IEEE 33 节点系统和IEEE 123 节点系统进行验证，改造后的系统结构如附录A 图A4所示。其中，由于μPMU 的布点策略不在本文的研究范围内，μPMU 安放位置参考文献［23-24］；负荷设置为综合负荷模型（恒阻抗负荷为20%，恒电流负荷为30%，恒功率负荷为50%），假设区域内负荷保持一致性，通过LHS 实现系统总负荷的随机波动，各节点负荷波动同总负荷保持一致；光伏机组设置为PQ 模型，功率因数为0.9，主要接在辐射状结构的末端节点或连接有较多支路的分支节点上，提高距电源较远区域的电压水平［25］，具体接入情况如附录A 表A1 所示。

按照2.1 节的拓扑生成规则和约束，分别在IEEE 33 节点系统和IEEE 123 节点系统中生成100 种可行拓扑进行识别效果验证。设置光伏出力和负荷在80%～120%波动，生成500 种运行场景。运行场景的数量说明如附录B 图B1 所示。将以上数据与可行拓扑组合，共形成50 000 个样本，经潮流计算生成包含节点电压信息的50 000 个样本数据。其中，45 000 个样本加入5%、10%和20%的随机误差作为BN 的训练集，1 000 个样本作为调试离散区间数量的调参样本，剩余4 000 个样本作为测试集。

在连续型数据离散过程中，光伏-负荷节点初始区间数设置为5，通过调参样本的识别效果寻找最优区间数目。其中，由于负荷和光伏数据通过LHS 生成，波动范围相同，在区间内均服从均匀分布，各变量之间无明显的相关性。在此前提下，可为负荷和光伏节点设置相同的离散区间数，以进一步减少离散过程的复杂度。离散区间的寻优过程如附录B 图B2 所示。在样本处理完毕后，按照1.3 节的BN 结构学习方法获取网络结构。以IEEE 33 节点系统为例，附录B 图B3 展示了获得的贝叶斯拓扑推理网络结构。图中，1 为拓扑节点，2 为负荷节点，3至6 为光伏节点，7 至30 为电压信息（幅值和相位）节点。在获得BN 结构后，通过极大似然估计法学习网络参数。至此，基于BN 的拓扑推理模型构建完毕。设置负荷、光伏和节点电压信息节点为输入节点，拓扑节点为输出节点，输入实时采集的负荷、光伏和节点电压等信息即可获得各个拓扑的概率，以概率最高的拓扑作为最终识别结果。

3.2 仿真验证

3.2.1 识别正确率测试

实际中，测量数据会受到外界噪声干扰。为了模拟真实情况下的识别效果，在μPMU 给出的电压数据和SCADA 系统/预测系统给出的光伏和负荷数据加入服从正态分布的噪声e1和e2，共设置了6 种噪声组合：①e1=0.5%，e2=1%；②e1=1%，e2=5%；③e1=5%，e2=10%；④e1=10%，e2=15%；⑤e1=15%，e2=20%；⑥e1=20%，e2=30%。其中，组合⑤和⑥为模拟测量装置不稳定的情况，可认为是数据时标未对齐的情形。分别用本文方法（BN）、估计匹配法（estimation matching method，EMM）［11-12］和投票法（vote method，VM）［14］进行测试，IEEE 33 节点系统下的识别效果如图2 所示。IEEE 123 节点系统结构和针对方法适应性的仿真结果如附录A表A2 和图A5 至图A7 所示，各项测试的总结果如附录C 表C1 至表C4 所示。

由图2 可知，3 种方法在测量误差较小的情况下均具有较高的识别率，但在模拟量测装置不稳定、误差较大时，均容易产生错误辨识。

3.2.2 识别时效性测试

为了验证所提方法在时效性的优势，分别在IEEE 33 节点系统和123 节点系统中生成100、200、300、400 和500 种可行拓扑进行测试。仿真环境为MATLAB2019b，处理器型号为Core i5-9400F，主频为2.90 GHz。3 种方法的时效性对比如图3所示。

图2 IEEE 33 节点系统下的识别效果对比Fig.2 Comparison of identification effects in IEEE 33-bus system

图3 方法时效性测试Fig.3 Timeliness test of method

由图3 可知，随着配电网可行拓扑数量增长，VM 和EMM 的识别时间呈线性增长，而BN 用时无明显变化。随配电网规模增大、测量数据和估计数据增多，3 种方法的识别时间均有所增长，但BN 仍保持在较低的范围内，有着明显的速度优势。

3.2.3 识别鲁棒性测试

为了验证所提方法在μPMU 故障情况下相较于其他方法的优势，设τ 为μPMU 的故障比例，在误差组合为e1=1%、e2=5%的系统中设置如下场景：τ=0、τ=10%、τ=20%、τ=30%、τ=40% 和τ=50%。在测试集中随机抽取100 个样本进行拓扑识别鲁棒性的验证。其中，为了保证测试的科学性，每种故障场景下进行10 次随机位置μPMU 故障的识别（故障数量为0 时为1 组）。2 种方法的测试样本相同，并取识别正确率的均值作为最终结果，测试结果如图4 所示。

图4 μPMU 部分故障下的拓扑识别Fig.4 Topology identification under partial fault of μPMU

由图4 可知，BN 在μPMU 存在故障的情况下仍能保证相对准确的识别结果，而EMM 和VM 效果受μPMU 故障数量的影响较大。这是由于2 种基于实时潮流计算的方法对数据敏感度较高，在可用数据量多时效果较好，但当测量装置故障或可用数据较少时，效果就大打折扣。本文所提方法基于BN的推理是一个通过先验信息（采集数据）不断更新后验信息（未知数据）的过程，通过仅有的证据可以实现对全局节点的推理，具有较强的鲁棒性。

除了上述μPMU 故障情况外，在光伏和负荷等影响潮流计算的关键数据存在缺失时，BN 仍能保证较好的效果。这是由于配电网的节点电压等信息与网络的注入功率本身存在着非线性关系［26］，而基于互信息的数据离散有效保证了BN 对光伏、负荷和节点电压信息相关性的挖掘与分析。即使这些信息不完全，所提方法也能通过有限的数据给出较为准确的结果。图5 展示了噪声为e1=1%、e2=5%以及光伏或负荷数据存在缺失时的识别效果。测试样本与μ PMU 故障测试相同，每种缺失率下设置10 组随机节点数据缺失情况（100% 缺失率下为1 组），取识别正确率均值作为最终结果。

图5 光伏/负荷数据部分缺失下的拓扑识别Fig.5 Topology identification under partial data missing of photovoltaic/load data

由图5 可知，BN 在光伏或负荷数据缺失率由0升至100%的过程中识别正确率有所降低，但即使缺失在100% 的情况下，识别正确率仍能保持在75%附近。EMM 和VM 的识别效果随缺失率的提高而大幅下降，在缺失率为20%时识别正确率低于30%，缺失率为60%时识别正确率低于10%，说明所提方法在关键数据缺失的情况下仍能保证对拓扑结构的准确识别，具有较强的鲁棒性。

4 结语

本文由离线方式计算不同拓扑在光伏-负荷场景下的潮流数据，并根据光伏、负荷、电压信息及拓扑的相关性样本搭建了BN。

1）与实时估计匹配方法相比，BN 具有更强的鲁棒性，在部分μPMU 存在故障、负荷或分布式发电等关键数据丢失的情况下仍能保证较高的识别精度，在实际中有较强的应用价值。

2）基于离线学习的BN 具有极高的时效性，且识别时间不随配电网规模的增大和可行拓扑的数目增大而明显增长，在大规模配电网拓扑识别应用中具有较强的技术优势。

所提出的BN 可以实现对配电网未知拓扑状态的快速和准确辨识。在下一阶段，将考虑负荷与分布式发电的关联性问题［27-28］，结合区域内光照和风速等信息对光伏-负荷和拓扑变化的潜在关系作进一步分析，并着重提升拓扑识别器的抗噪声能力。此外，将考虑μPMU 的优化配置情况，在保证方法效果的前提下尽可能减少放置的μPMU 数量。

附录见本刊网络版（http：//www.aeps-info.com/aeps/ch/index.aspx），扫英文摘要后二维码可以阅读网络全文。