Cyber-Physical Coordinated Assessments with Incorporations of Transient Stability Analysis

随着能源互联网建设和电网数字化转型加快，电力信息物理融合系统（CPS）正加速电力运行调度网、全域安全防护网和用户终端数据网之间的协同互动，促进建设以“电力大数据服务公共社会管理和经济发展”为目标的新型数字产业，发展以高比例分布式能源为主体的现代化综合能源系统^[1-2]。作为关键支撑技术，自2006年首次由美国国家科学院提出以来，CPS的安全性定义和研究边界一直存在分歧。一般而言，电力CPS的安全问题是包含通信控制安全和工程运行安全的综合安全问题^[3]，其研究本质是基于通用的数学模型描述信息系统的通道中断、信息错位、延迟以及错误等故障，在通信智能化技术高度集中的CPS模型中，对物理侧能源系统的工控可靠性影响^[4]。

国内外电力CPS安全研究学者在信息安全脆弱性分析^[5]，基于网络攻击事件的电力运行故障分析^[6]，电力CPS离散事件的控制建模^[7]，电力CPS的半实物仿真^[8]等领域开展了大量研究并取得多项成果。但上述部分研究的重点大都落在信息侧或物理侧的单边安全问题上，未能提出改善CPS综合安全性的量化方案，其中部分结果缺少电力工程安全环境的计算校验，安全结论还需验证^[9]。电力CPS的安全问题本质可类比“木桶理论”，其上限应由整体系统的安全短板决定，即需要考虑信息物理双重脆弱性的协同作用^[6]。CPS环境中任何信息安全漏洞或工程设计缺陷都有可能被恶意利用，其危害会被放大，直接影响到电力系统综合安全。筛选信息安全事件及其潜在诱发的工程安全事件，有助于快速定位CPS系统内部的安全“短板”，帮助安全研究人员在规划设计阶段提前做好预案^[10]。

为了更好解决电力CPS综合安全问题，学者前期工作^[10-11]介绍了基于电力监控系统搭建入侵攻击模型，给出潜在入侵风险的概率建模，并根据变电站中断故障响应给出基于稳态潮流分析的全枚举方案。但上述工作缺少CPS联合仿真模型校验，缺少信息安全事件对工程物理安全的直接表征。Liu在CPS融合仿真测试中明确了中间人攻击（MIM）、拒绝服务攻击（DoS）等网络攻击诱发的非传统电力工程故障的复杂特性^[12]，在此基础上，学者基于RTDS搭建了广域测量网络（WAMS）架构的半实物仿真系统，将控制设备和测量装置的通讯协议整合到仿真系统中，指出基于CPS模型的控制理论可有效限制电力系统级联故障^[13-14]。上述研究采用的OPNET和RTDS/RT-LAB混合仿真平台可以实现多组高度复杂模型的实时动态仿真^[15]，为相关理论及应用提供了坚实的技术支撑，但由于其投资成本高，开发周期长，不支持后期开发，目前暂未普及^[16-17]。

为提高仿真平台的适用性和可拓展性，学者基于软件定义网络（SDN）理论，将电力控制采集元件抽象为独立单元，自定义攻击场景，验证了IEC 61850协议的网络攻击特征^[18]。为保证研究结果易于复现，安全结论便于推广，文章通过开源软件搭建了电力变电站CPS联合仿真的初级测试平台，主要贡献包括：（1）基于实际工程的变电站监控系统工程接线，通过网络拓扑仿真器Mininet搭建了变电站过程层网络的虚拟化框架^[19-20]，在虚拟化的通信主机中调用数据包分析进程Scapy设置通信系统故障^[21]；（2）通过调用外部的PSD-BPA暂/稳计算核心，脱离了传统的“故障卡”限制，实现了计算文件中故障信息的自动更新，实现了基于Mininet-BPA的CPS轻量级联合仿真，并为安稳分析工程师提供在线实时计算多组电力故障的可行思路；（3）根据典型的网络安全事件，通过实际的工程算例验证了基于电力CPS联合仿真框架的多站点过载级联故障的快速检测方案，并给出了对应的安稳指标。

基于SDN理论，网络模拟器Mininet提供了一种集成路由器、交换机、链路、主机等多个对象联合建模的整体虚拟化方案^[19]。本章以电力监控系统为基础，介绍基于Mininet-Scapy框架的轻量化信息网络拓扑仿真方案。图1左侧给出了电力变电站母线母联间隔中，过程层网络拓扑的CPS简化模型，图1的右下框图为对应的信息模块及接线，右上框图为变电站过程层网络的Mininet虚拟化架构。在Linux环境中，通过仿真器Mininet虚拟化的各个设备和主机在Linux中被表示为各个Bash进程，共享核心、主机文件系统以及进程标识号。但进程在网络命名空间中会被分配专用的网络接口、路由、防火墙规则等，使得各个Bash进程可以运行在独立且完整的网络环境中，运行的系统级代码和功能也互不影响。由此，文章认为基于Mininet创建的虚拟网络可以较为准确地体现真实物理主机的网络特征。

Figure 1.  The model of power substation CPS and the corresponding Mininet-based virtualization framework

图1中基于Mininet模拟器的信息仿真框架里除了与真实系统对应的信息主机和交换机节点外，另有控制器（Controller）节点。这是由于Mininet是基于SDN开发的项目，后者典型特征是将信息网络的各控制功能抽象并实例化为一个单独的功能模块，用户可通过设置控制器参数来配置各个模块的详细功能，以此测试用户定制化的通信或路由协议。文章不涉及此项功能，选择默认交换机模型，即控制器参数设置为远程控制器方案“Controller”设置为“Controller-Remote”，即不额外定义本地的控制器模块。图1中控制器通道为虚线，其余通信链路为实线，以示区别。

为了提高仿真测试的准确性和适用性，根据IEC 61850标准，文章提出的变电站信息仿真框架中还包括了SV和GOOSE通信协议，以此区分合并单元和智能终端的通信内容。文章选择了开放、交互式的数据包开发程序Scapy作为构建和收发GOOSE（SV）数据包的功能模块，通过配置其内置函数参数，可自定义收发主机的目的IP、MAC地址以及网卡信息等。根据前文介绍，基于Mininet虚拟化部署的主机节点可独立调用虚拟机的核心进程和应用，由此各主机调用、运行Scapy的数据包收发功能不会互相影响。

图2给出了基于Mininet-Scapy联合仿真的变电站信息网络虚拟化流程。主要功能伪代码见附录A。注意到，代码中调用了专业的网络仿真平台NS3作为基础仿真环境，这是因为NS3支持更多通信功能和协议，且能为Mininet提供外部接口，让Mininet主机运行NS3内部的通信协议，提高仿真准确性。例如，附录A中的链路仿真模块“CSMALink”便是基于NS3开发。详细的Mininet嵌入NS3融合开发技术不属于文章范畴，此处不展开讨论。

Figure 2.  Simulation solutions of the information system based on Mininet-Scapy framework

由变电站信息网络虚拟化框架可知，Mininet各个主机可以查看其通信及控制详情，直接调用Wireshark程序即可查询所有网卡上的数据流量，进而得到网络的时延信息。本章将继续介绍基于Mininet-Scapy信息通信框架的系统运行安全分析模块，并根据线路的过载结果给出系统安全风险的量化方案。根据工程安全的分析思路，文章基于“最坏”考虑给出安全事件$ \mathit{X} $假设：变电站$ i $同时出现电力故障和信息故障。例如，攻击者劫持变电站通信网络且持续潜伏，直到变电站出现全站失压故障，并同时发动分布式拒绝服务攻击（DDoS），恶意占用合并单元和智能终端的传输通道，以拖延测控保护装置接收故障信号，从而扩大故障影响。

便于阐述，本节令网络攻击造成的时延为${\tau }_{{\rm{x}}}$，系统变电站集合${\boldsymbol{N}}$，线路集合${\boldsymbol{L}}$，借助图论形式${{{\boldsymbol{G}}}}({\boldsymbol{N}},{\boldsymbol{L}})$来表达系统的电网拓扑，记安稳计算程序为$ \mathbf{c}\mathbf{a}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{l}\mathbf{a}\mathbf{t}\mathbf{e}(\bullet ) $，其结果${\boldsymbol{r}}$为高阶矩阵，包括线路最高功率、额定功率、系统功角、系统频率以及系统电压的计算结果，分别记为${P}_{\mathrm{m}\mathrm{a}\mathrm{x}}，{P}_{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}，\delta \text{´}，\mu \text{´}$以及$\sigma \text{´}$，可由.read()函数展开。算法1给出了传输线过载引起的级联故障的工程安全计算方法，如表1所示。

算法1第4行给出了系统的过载指标$ {\mathrm{\gamma }}_{\mathrm{p}} $，其中$\kappa $为过载稳定阈值。$ m $表示当发生第$ n $组全站失压故障时，过载线路集合${{{\boldsymbol{M}}}}$的索引，${P}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{{\rm{m}}}}$和${P}_{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e},\mathrm{{\rm{m}}}}$分别表示线路$ m $出现的最高功率和额定功率。此外，算法还设置了功角稳定性指标${\gamma }_{\text{δ}}$，频率稳定性指标${\gamma }_{\text{μ}}$以及电压稳定性指标${\mathrm{\gamma }}_{{\text{σ}}}$，初始值设为0。算法1还选取了最大系统功角、最低系统频率和最低系统电压，即${\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}、{\mu }_{\mathrm{m}\mathrm{i}\mathrm{n}}$和$ {\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}} $，作为表征各分量严重性的稳定判据。例如，在计算周期内，$ tϵ(0,{t}_{\mathrm{d}}) $，若系统的功角差超过$ {\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}} $，则给出判断：系统在下个计算周期内极有可能出现失步运行。相较于线路过载，失步运行的故障影响更为严重，于是将其对应的功角严重性指标${\mathrm{\gamma }}_{\text{δ}}$设置为1.0。同样地，指标${\mathrm{\gamma }}_{\text{μ}}$和${\mathrm{\gamma }}_{\text{σ}}$设置为1，分别表示仿真中算法检测到系统在该计算周期内或下个周期将发生严重的频率降低和电压崩溃故障。算法给出${v}_{n}$作为第$ n $组故障的综合安全指标，可知，其数值越高，故障越严重，即系统稳定性越低。功能实现的伪代码见附录B。

由附录B可知，文章选择的安全计算程序$ \mathbf{c}\mathbf{a}\mathbf{l}\mathbf{c}\mathbf{u}\mathbf{l}\mathbf{a}\mathbf{t}\mathbf{e}\left(\bullet \right) $包括了工程计算中常用的电力系统安全分析软件PSD-BPA中内置的稳态/暂态计算核心“bpa_pfnt.exe”和“bpa_swnt.exe”，主要用于工程计算中的稳态潮流计算和暂态故障计算。此外，在update_swi()函数中，通过Python脚本实现了故障卡的连续、自动编辑和更新，无需手动编辑“.LSD”故障卡。由电力工程安全的稳控策略表可知，直接切除过载线路容易扩大故障影响，过载保护设备动作较为谨慎。由此，在算法实现中，暂态分析函数“Xcade_swi()”中只选取了过载最严重的线路作为典型的过载故障信息添加到故障卡更新函数“update_swi()”中，用以区分初始故障。

根据前文第1章和第2章的介绍可知，文章基于Mininet虚拟化框架联合电力安全计算软件PSD-BPA，得到了变电站信息物理联合测试方案。本小节将详细介绍仿真的测试环境及相关参数。根据南方电网系统某220 kV保信子站的工程接线图，算例选取了变电中典型的保护间隔作为研究对象，给出了联合仿真方案，如图3的右上框图所示。文章的算例仿真主要基于VirtualBox-Ubuntu搭建的Linux环境。算例在已嵌入Mininet模块的NS3环境中启用“./waf”环境变量的编译模式，并在该模式中直接运行Mininet脚本即可让虚拟主机接入NS3环境。图3中左上框图给出了变电站通信网络虚拟化框架及相关界面展示：配置NS3-Mininet的联合仿真通道，通过Python编辑Mininet的网络仿真脚本，通过“net.addSwitch()”函数配置1个主交换机和5个间隔交换机，通过“net.addHost()”函数配置了20个通信节点，并根据图3的间隔信息，通过Mininet中内置的链路配置函数“CSMLink()”配置了34条链路。

Figure 3.  The CPS-based co-simulation scheme of power system analysis on delayed-operations of protective relay induced by information network failures in the 220 kV power substation

算例基于第2章中的“最坏”假设，选取Scapy作为收发流量测试包的操作工具，通过内置的“sendp()”函数定义发包规则、网卡和间隔时间来模拟满足IEC协议的GOOSE和SV报文通信^[22]。在目的地址的主机上调用Wireshark抓取特定网卡的报文，解析并处理报文。最后，算例将变电站信息网络仿真结果导入至算法1，并调用PSD-BPA的计算核心实现信息物理联合仿真。仿真中主要的编译语言为Python 3.7，暂稳计算核心版本为5.7.1。

1）变电站通信网络链路的时延测试结果

由3.1节环境设置可知，算例通过Mininet添加了20个通信节点（host），记为${h}_{1}，{h}_{2}，\cdots ，{h}_{20}$。其中，${h}_{2}，{h}_{3}，{h}_{4}，{h}_{5}$分别表示图3中的母联测控主机保护装置，母线保护装置，变压器保护装置以及线路保护装置；${h}_{6}，{h}_{7}，{h}_{8}$依次表示母线母联间隔内1套合并单元和2套智能终端的3个主机配置；以此类推，${h}_{9}，{h}_，\cdots，{h}_{20}$表示其余4个间隔内的主机配置。表2给出了未受攻击时的系统延时的测试结果。结果显示，20个通信节点互相通信时，部分节点和链路的负载较高，某一条链路会出现极高时延，可以通过配置冗余链路或有针对性配置低时延设备来改善。但在网络入侵攻击的假设场景下，为了体现入侵攻击对系统时延的影响，算例选取15 ms作为系统设备的默认时延设置，以便消除由通信链路单一或部分节点负载高等客观因素对网络安全结论的影响。

算例在Mininet中搭建了220 kV变电站的网络虚拟化模型。表3给出了变电站保护装置在不同测试场景下的平均时延结果，其中保护1～保护5分别为母联保护、110 kV变压器保护、35 kV变压器保护以及厂站进线与出现的线路保护。各组保护系统内配置了1套合并单元和2套智能终端，根据附录中的设置，初始带宽设置为1 Gbps。为区分不同的设备，保护1～保护3的固定时延设置为5 ms～20 ms，保护4和保护5的固定时延设置为40 ms～50 ms。注意到，时延仿真测试中某条链路出现了较大时延，这是因为此处采用的是基础通信模型，网络拓扑未做优化处理。${\mathit{X}}_{1}，{\mathit{X}}_{2}$分别表示保护1和2子系统中的智能终端发生了信息安全事件$ \mathit{X} $；$ {\mathit{X}}_{3} $表示保护2和保护3子系统中同时出现了信息安全事件。

由信息仿真结果可知，出现信息故障$ {\mathit{X}}_{3} $后，通信网络中约64%的链路时延增加了10倍或以上，其中2号主变的MU主机与110 kV变压器保护主机的通信时延，从正常条件下的5.16 ms增加至520.11 ms。此时信息故障将极大地阻碍合并单元向保护装置传递测量信号，存在严重运行安全隐患。

2）基于海南电网夏大运行方式数据的仿真结果

算例选取2027年海南电网夏大运行方式数据作为系统工程安全的基础数据。注意到，电网中接入78个变电站站点，包括466个母线节点，涵盖23个供电区域，接入功率11.6 GW，挂有11.3 GW负荷。算例假设各站点的信息物理拓扑与配置均同图3一致，并选择表1中场景$ {\mathit{X}}_{3} $作为信息故障模型。可知，通信系统的平均时延为288.5 ms，保护设备接收到真实测量信号的时间将延缓至少12～13周波。过载稳定裕度$\kappa $取温度在15～40 ℃时载流量的变化区间^[23]，$\kappa $$\in \left[\mathrm{1,1.35}\right]$。

根据上述计算条件，首先对各站点的母线出线做N−1三相短路故障校验，共计检测线路1400条，其中无法满足计算条件的线路200条，即在预设的计算周期500 ms内，电力系统的安稳计算无法在约束条件内求解可行运行方案。在严重的信息故障前提下，工程安全故障范围可能会被恶意扩散，此时电力系统极有可能不再满足基本N−1准则，电网可靠性降低。

基于算法1，图4给出了海南电网各个变电站站点在信息故障场景$\mathit{X}_{3}$下发生全站失压故障的联合仿真结果。其中红色和蓝色图例分别表示发生初始故障后，有/否考虑线路过载特性的暂态稳定结果。黑色表示暂稳安全指标为1.0的厂站。例如，“望楼站”为220 kV厂站，挂有负荷179.8 MW，当该站同时出现信息故障和全站失压故障时，注意到，计算周期内暂态分析仍然收敛，系统的整体过载指标为0.785，且未检到其他安全故障。但在考虑过载保护动作的场景中，系统在发生初始故障后会切除部分严重过载的线路，此时系统的功角差将超过计算极限，在计算周期内无法得到可行解，预计将出现严重的失步运行故障，此时该站的稳定性指标为1.0。同样地，“八所站”为500 kV变电站，挂有负荷215 MW，在发生初始故障时，计算得到系统的过载指标为0.73；当考虑过载保护动作后，系统的功角稳定性和频率稳定性指标均为1.0，预计将出现严重的失步运行和频率崩溃故障。

Figure 4.  The screening results of cyber-physical coordinated evaluations based on the peak-load dataset of Hainan grid

注意到，在严重的信息故障条件下，大部分厂站的过载保护动作将降低系统运行的安全阈值，其中“望楼”“炼化”等11个站点出现了包括潜在的失步运行、电压崩溃故障以及系统频率降低甚至频率崩溃等严重运行故障。但对于“文昌气电”站、“东方电厂”等厂站，过载保护的积极动作提升了其运行系统的安全阈值。下一步的工作将研究详细的安稳控制策略对CPS工程安全问题的影响机制。

文章沿用了前文的网络攻击假设，在站控层网络的入侵风险评估模型的基础上新增了过程层网络的信息物理联合仿真，基于继电保护装置的时延测试验证了入侵风险的评估方法。文章提出了可快速筛选潜在“高危风险节点或设备”的CPS联合仿真方法，为电力CPS安全工程师提供了一种支持自定义开发及调试系统安全的开源方法。基于实际的工程算例，文章成功识别出高危风险站点和线路，量化了安稳控制设备的信息故障对电力能源系统的影响，同时验证了过载保护等传统电力保护设备对CPS安全事件的影响机制，例如，对负荷节点的和发电机节点应配置不同的过载保护控制策略便可提升系统的稳定裕度。此外，文章在联合仿真中脱离了传统的“故障卡”限制，实现了计算文件中故障信息的自动更新，在提高计算效率的同时也为安稳分析工程师在线实时计算多组电力故障提供了可行思路。

未来能源领域的工作重点是搭建以新能源为主体的高效能源体系。更多智能化、集成化的控制和感知设备将促进电力运行网络与信息网络的数据融合，电力CPS的安全场景愈发复杂。后续的研究将继续优化联合仿真功能代码，提高数据交互效率，同时探究安稳控制设备的信息安全事件对高比例新能源广泛接入的新型电力系统的影响。

该部分主要介绍Linux环境中基于Mininet各子模块开发的变电站信息仿真平台。

# 输入参数:

# H,主机节点list

# S,交换机节点list

# L,系统链路list

# DataRate,带宽

# Delay,基础时延

from mininet.net import Mininet, Controller, CLI,ns3

from mininet.ns3 import CSMALink

def SubstationNets_Initialize(H, S, L, DataRate, Delay):

net=Minine t(controller=‘Controller-Remote’, autoSetMacs=True)

net. addController('Default_Controller')

for h in H:

h=net. addHost(h. name, ip, defaultRoute=None)

h. setIP(h. IP, intf=h. MAC)

h.setMAC(h.MAC, intf=h.MAC)

for s in S:

s.addSwitch(s.name)

for l in L:

CSMALink(l,DataRate,Delay)

net.start()

mininet.ns3.start() // 激活NS3环境

CLI(net)

If__name__==‘__main__’:

SubstationNets_Initialize(H,S,L,DataRate,Delay)

# 计算参数:

# tao_x,信息仿真的时延结果

# sub, 系统各变电站的母线节点汇总list

# 初始化delta_max, miu_min, sigma_min, kappa, t_0, t_d,

import os

def calculate():

os.run(‘bpa_pfnt.exe’)

os.run(‘bpa_swnt.exe’)

if__name__==‘__main__’:

gamma_delta, gamma_miu, gamma_sigma=0

for n in sub:

update_swi(n,tao_x) # update_swi()修改.swi文件，模拟

r=calculate() # 系统切除全站失压故障

p_list=r.read() # p_list储存系统过载线路信息

for p in p_list:

if p._overloads_time in [tao_x,t_d] # 调取线路过载信息，并判断其是否

Xcade_swi(p_list) # 在程序的计算周期内，Xcade_swi()

r_Xcade=calculate() # 可在原.swi卡中添加过载线路信息

p_list,_delta, _miu, _sigma=r_Xcade.read()

if _delta>delta_max: gamma_delta=1

if _miu<miu_min: gamma_miu=1

if _sigma<sigma_min:gamma_miu=1

if gamma_delta==1 or gamma_miu==1 or gamma_miu==1: pass

else

p_rate, p_max = r.read(p_list)

gamma_p = sum(p_max)/kappa*sum(p_rate)