图1 直驱风机并网系统结构
Fig.1 The structural diagram of PMSG connected to grid
摘要 电力电子化电力系统的宽频振荡问题严重危害电网的安全稳定运行，及时确定宽频振荡的传播路径和振荡源位置，对振荡的抑制至关重要。该文从因果关系的角度提出一种基于Copula传递熵的设备级和网络级宽频振荡传播路径分析及振荡源定位方法，通过分析多变量之间的因果传递方向和因果强度系数，利用有向加权图构建宽频振荡因果网络，在此基础上确定振荡源位置，并通过保留因果强度系数最大的支路，确定振荡的主要传播路径。采用所提方法一方面能够从设备级分析控制器内部各状态变量的因果关联性；另一方面能够从网络级分析电网中各节点振荡数据的因果传递性，从而确定宽频振荡在控制器内部和电网中的传播路径及振荡源位置。最后，在强迫振荡、风机与弱电网交互、发电机轴系振荡中的仿真算例表明，所提方法适用于多种机理的宽频振荡，能够同时实现设备级和网络级宽频振荡传播路径分析和振荡源定位。
关键词：宽频振荡 振荡传播 振荡源定位 Copula传递熵 因果分析
在“双碳”目标下，风电等新能源的大规模并网[1]、柔性直流/交流输电工程的建设[2]和大功率直流负荷的接入[3]等为电力系统引入了大量电力电子设备。而多样化电力电子设备在系统中的动态交互作用，可能会引发电力系统宽频振荡，其在电网中的传播可能造成电网联锁事故的发生，严重危害电力系统的安全稳定运行[4-5]。及时确定宽频振荡的传播路径和振荡源位置，对于振荡的抑制至关重要。
根据是否需要系统数学模型及参数，宽频振荡分析方法可以分为基于模型的分析方法和基于数据的分析方法。特征值分析法是最常用的基于模型的分析方法，通过建立系统的状态空间模型，并求解其特征根，进而对宽频振荡进行分析。文献[6]基于特征向量的物理含义提出一种确定振荡传播路径的分析方法，该方法能够确定参与振荡的主要设备及其控制环节。文献[7]利用闭环传递函数对直驱风电场内部及直驱风电场和弱交流系统之间的次同步分量通路进行研究，并分析了控制器参数对通路阻尼特性的影响。此外，还可以采用阻抗法，将系统等效为一个串并联的阻抗网络，并对其进行分析。文献[8]通过建立直驱风机并网系统的阻抗模型，推导了变流器控制参数对系统稳定性的影响。然而，对于含有大量非线性环节的电力电子化电力系统，基于模型的分析方法存在参数未知、模型阶数过高等局限性[9]，难以在实际系统中广泛应用[10]。
随着电力系统量测技术的发展和大量同步相量测量单元的安装[11]，基于量测数据对宽频振荡进行分析研究成为一个新思路，此类方法能够有效地摆脱对系统模型的依赖，因此受到广泛关注。常用的基于数据的方法包括阻抗测量法[12]、能量法[13-14]和人工智能方法[15-16]。文献[12]提出了确定次/超同步振荡源的阻抗判据，利用阻抗测量法实现了振荡源的定位。然而，阻抗判据受到阻抗耦合现象的影响，可能会导致分析结果不准确[17]，此外，在线测量的阻抗精确性容易受到扰动情况、系统运行状态、数据质量、辨识算法等多方面的影响[18]，从而影响定位结果。文献[13]分析了次同步振荡中暂态能量流和阻尼的关系，通过计算网络中的振荡能量流确定了次同步振荡源的位置。文献[14]提出了一种基于小波耗散能量谱的强迫振荡源定位方法，可利用系统的实测数据定位强迫振荡源。能量法在系统中能够准确定位的前提是无损网络和恒功率负荷[19]，而这样的前提在实际系统中通常难以满足。此外，基于能量法的母线级定位判据最初是以同步发电机组为核心展开的[13]，对不同机理宽频振荡的适应性需要开展更深入的研究。人工智能近年来迅速发展，文献[15]采用深度迁移学习和对抗式迁移学习方法，基于瞬时的振荡数据，实现了风电场中次同步振荡源的定位。文献[16]采用稀疏贝叶斯学习方法对发电机等效阻尼系数的正负进行评估，实现了负阻尼振荡源的定位。基于人工智能的方法对数据的规模和质量有较高的要求，然而目前实际系统中可获取的宽频振荡样本十分有限，因此人工智能的宽频振荡分析方法如何在实际系统中可靠应用还有待进一步研究。
随着信息理论的发展[20]，基于数理统计的方法为宽频振荡的分析提供了新的思路。此类方法对数据间的关联进行研究，发掘出振荡的内在规律。文献[21]利用互信息建立了次同步振荡各影响因素的评价模型，分析得到了振荡的关键影响因素。文献[22]利用偏有向相干分析方法，量化评估了双馈风机内部各能量元件的因果关系，得到了振荡在各元件间的传递方向和振荡源位置。目前基于数理统计的宽频振荡分析方法研究较少，考虑到新型电力系统中宽频振荡数据为多变量的强非线性和因果性的时间序列，如何从基于数理统计的因果关系角度分析宽频振荡的传播特性、实现振荡源定位值得深入研究。
本文首先从传递函数的角度对宽频振荡在设备内部和电网中的传播特性进行了理论分析，在此基础上考虑宽频振荡的非线性，采用适用于多变量非线性复杂因果分析的Copula传递熵[23-25]构建因果网络，一方面从设备级分析控制器内部各状态变量的因果关联性，另一方面从网络级分析电网中各节点振荡数据的因果传递性，在此基础上确定振荡源位置，并通过保留因果强度系数最大的支路，从而确定宽频振荡在控制器内部和电网中的主要传播路径。本文提出的基于Copula传递熵的宽频振荡传播路径分析和振荡源定位方法无模型假设，对数据分布没有要求，且具有较为明确的物理意义。针对不同频率（次/超同步、中高频）和多种机理（强迫振荡、风机与弱电网交互、发电机轴系振荡）的振荡算例表明，所提方法适用于多种机理的宽频振荡，且可以同时实现设备级和网络级宽频振荡传播路径和振荡源位置的确定。
以图1所示直驱风机并网系统为例对设备级的宽频振荡传播特性进行分析。其主要包括永磁同步机（Permanent Magnet Synchronous Generator, PMSG）及其传动系统、机侧换流器（Machine Side Converter, MSC）及其控制系统和网侧换流器（Grid Side Converter, GSC）及其控制系统。其中，Pin、Pe分别为机侧输出功率、网侧输入功率，Cdc为直驱风机直流电容，Udc为直流电容电压，utd和utq分别为GSC端口电压d、q轴分量，ugd和ugq分别为并网点电压d、q轴分量，igd和igq分别为GSC输出电流d、q轴分量，Kpv和Kiv分别为电压外环的比例、积分参数，Kp、Ki分别为锁相环的比例、积分参数，Kpi、Kii分别为电流内环的比例、积分参数，上标“c”为控制系统中的变量，上标“*”表示该变量的参考值，Lf为滤波电抗。
图1 直驱风机并网系统结构
Fig.1 The structural diagram of PMSG connected to grid
由系统的拓扑结构可以得到小信号方程为
（2）
式中，“
”表示小信号量；
为电网频率。
经过拉氏变换，式（2）可以写为
（3）
在合适的控制方式下，直驱风机的机侧部分不参与其网侧部分的振荡，因此本文只考虑网侧换流器，由其控制框图可得
（4）
（5）
其中，控制系统中的变量和交流系统中的变量关系为
（6）
（7）
联立式（5）～式（7），得到
（8）
其中，锁相环输出相位扰动为
（9）
令Hdc(s)、Hu(s)、Hi(s)和Hpll(s)分别为直流电容环节、直流电压外环、电流内环和锁相环的传递函数，其表达式为
（10）
联立式（3）～式（5）、式（8）～式（10），可以得到直驱风机并网系统的传递函数框图，如图2所示。
图2 直驱风机内部传递函数框图
Fig.2 Internal transfer function block diagram of PMSG
图2传递函数框图反映了直驱风电场内部各状态变量之间的传递关系，当系统存在频率为ωd的振荡模态时，令两变量间的传递函数的模值为该模态在两变量间的传播强度，见表1。
表1 各状态变量间的振荡传播强度
Tab.1 Oscillations intensity between state variables
传播强度 ①② ②③ ③⑤ ④⑥ ⑤⑦ ⑥⑧ ⑦⑨ ⑧⑩
通过以上分析可以实现宽频振荡在直驱风机内部传播特性的定量计算，对于状态变量
和
，到的传播强度越大，表明对的响应越大，状态变量到越是主导传播路径。因此，在知道振荡源和风机模型及其详细参数的条件下，可以根据状态变量间传播强度的大小确定振荡传播的主导路径。
对于一个包含风机、同步发电机和柔直等环节的复杂电力系统，联立其风机动态方程、同步发电机方程、柔直控制方程和系统网络方程，可以得到整个系统的状态空间方程，即
（11）
式中，A为状态变量矩阵；B为输入矩阵；C为输出矩阵；Δx为系统的状态变量向量；ΔiR为注入电网的扰动向量；Δy为系统各母线的输出向量。对上述状态空间方程进行拉氏变换，可得
（12）
联立得到传递函数矩阵，即
（13）
当风机与系统交互产生宽频振荡，振荡分量会反映在其并网点的电压、电流中，并会在系统各个节点得到响应。可利用整个并网系统的多输入多输出传递函数，从频域的角度来定量分析宽频振荡的传播特性。
定义振荡传播后在系统各节点输出变量幅值的变化程度为该节点的振荡传播强度[26]，即
（14）
式中，Wij(s)为传递函数矩阵的元素，表示第i个节点输出变量对输入扰动ΔiR的响应；║∆iR║2表示输入扰动向量的2范数；
为输入扰动向量的元素。hi(s)越大，表示扰动对该节点输出变量幅值的影响越大，则该节点的振荡传播强度越大，扰动位置到该节点的位置即为振荡传播的主要路径。因此，若系统的状态空间方程及参数已知，则可以通过上述过程实现宽频振荡的网络级传播特性分析。
通过第1节的分析可知，从传递函数的角度对宽频振荡传播特性进行理论计算仅适用于系统模型和参数已知且振荡源已经确定的情况，难以对模型未知的实际系统振荡的传播路径和振荡源位置进行确定。因此本节从数理统计的角度提出了一种基于Copula传递熵的设备级和网络级宽频振荡传播路径分析及振荡源定位方法，在系统模型未知的情况下能够实现振荡传播路径分析及振荡源定位。
在不同的理论中宽频振荡的因有不同的定义，从阻抗的角度来说，宽频振荡的因是等效电阻为负的元件，从能量的角度来说，宽频振荡的因是耗散能量为负的元件，从功率的角度来说，宽频振荡的原因是该频率下消耗功率为负的元件。以上因果关系均基于一定假设下物理模型的推导，这些方法在高比例新能源和高比例电力电子设备系统中的应用也存在一些挑战。因此，本文从信息论的角度对宽频振荡的因果关系进行定义：如果某个振荡变量的历史状态对其他变量的当前状态存在很大的影响，则说明该变量的动态变化为其余变量动态变化的“因”，其他变量为“果”。
从网络级角度来看，电力系统中设备与电网在公共连接点（Point of Common Coupling, PCC）交互引发振荡，此时系统中其他节点的电压电流根据电网运行条件和网络在该振荡频率下的参数发生相应改变，此时PCC的振荡可以看作是因，网络其他节点的振荡可以看作是果。从设备级角度来看，当电力电子设备内部某个控制环节发生故障或者参数设置不合理时，可能会导致设备内部各状态变量的振荡，此时该控制环节输出的状态变量可以看作为“因”，其他状态变量可以看作为“果”。
为了对宽频振荡中因果分析的适用性进行验证，本节首先从电力系统线性模型入手，分析了传递函数和格兰杰因果[27]（Granger-causality, G-causality）的对应关系，并采用传递熵的方法将因果分析进一步推广至非线性的宽频振荡问题分析中，具体分析过程如下。
对于电力系统线性模型，其状态空间方程的离散形式为
（15）
式中，D为直接传递矩阵。
用滞后算子L表示
，将xt和yt相互表示，则式（15）可写为
（16）
（17）
在格兰杰因果关系检验中，对于X和Y两个变量，假定被解释变量和解释变量每一个变量的预测信息全部包含在这些变量的时间序列中，检验通过估计如下回归实现：
（18）
（19）
式中，Xt、Xt-j分别为变量X在t和t-j时刻的值；Yt、Yt-j分别为变量Y在t和t-j时刻的值；ext、eyt分别为变量X和Y的随机误差；aj、bj、cj、dj为各参数矩阵的元素。
若式（18）中X的稀疏估计值在统计上整体显著不为0，同时式（19）中Y的稀疏估计值在统计上整体显著为0，则称X是引起Y变化的原因；反之亦然。当变量X和Y之间满足式（16）、式（17）时，则可以通过式（18）、式（19）的格兰杰因果关系检验对其进行因果关系的分析。
然而，格兰杰因果仅适用于线性因果关系，因此含有大量非线性环节的电力电子化电力系统宽频振荡问题需要一种新的因果分析方法。传递熵是分析非线性系统中因果关系的重要工具，它检测方向性和动态信息，但不采用任何特定的函数形式来描述不同变量之间的相互关系。因此，本节采用适用于多变量非线性复杂因果分析的Copula传递熵来分析宽频振荡的传播问题。
传递熵用来度量随机过程之间定向信息的传递。假设
、
为两组平稳的马尔可夫时间序列，T为时间序列的长度，X到Y的传递熵
定义为
（20）
其中
（21）
（22）
（23）
式中，
和t分别为X和Y的滞后时间；k和l分别为X和Y的嵌入维度；
/
为联合/条件概率密度函数。如果在未来的某一时刻，过程Y的状态完全由本身的历史状态决定，则此时
。传递熵也可以采用熵的形式进行表达，即
（24）
式中，
为变量间的联合熵；
为变量的边缘熵。式（24）表明，传递熵可以排除X和Y受到共同的历史状态和外部环境的影响，只反映X的历史状态对Y的当前状态的影响。
实际上，由随机变量联合熵、边缘熵和Copula熵的定义，可以推导出
（25）
式（25）表明联合熵能够分解为两部分：两个边缘熵和一个Copula熵。如图3所示，图形总面积为联合熵，两个外侧圆环分别为两变量X和Y的边缘熵，中心部分即为度量X和Y相关性的Copula熵。
>
图3 边缘熵、Copula熵与联合熵关系示意图
Fig.3 The relationship among marginal entropy, Copula entropy and joint entropy
因此，根据以上关系，将Copula熵引入传递熵的表达式中，即
式（26）表明传递熵可直接表示为三项 Copula 熵的组合，并不依赖边缘分布和联合分布的具体形式。由文献[24]可知，Copula熵的定义为
（27）
式中，
为变量X和Y的Copula密度函数，由Copula函数分别对X和Y求导求得，可以表示为
（28）
式中，
为变量X和Y的Copula函数，根据数据集，采用秩统计[28]的方法进行参数估计。
由于Copula传递熵具有不对称性，即
，因此可以通过比较两者的大小来判断因果关系的方向性，并定义其中较大的数值为因果强度系数。若传递熵＞
，则X是Y的原因，此时因果关系的传播方向是
，因果强度系数为的数值；反之，若传递熵＜，则Y是X的原因，此时因果关系的传播方向是
，因果强度系数为的数值。
宽频振荡在电力系统中的因果传播关系可以用有向加权图来表示，即
（29）
式中，V为网络节点集合，在设备级振荡分析中，由控制器的各响应状态变量组成，在网络级振荡分析中，由各母线的功率信号组成；W为各边的边权集合，权值为各节点间的因果强度系数，通过对经过预处理后的振荡数据计算Copula传递熵获得，即
（30）
式中，Wij为节点vi与节点vj之间的因果强度系数。根据Wij可以绘制出宽频振荡因果网络。
在宽频振荡因果关系网络的基础上，计算网络中每个节点vi的出度Dout(vi)，其计算公式为
（31）
某个节点的出度Dout(vi)越大，表明该节点的历史状态对其他节点之间的当前状态影响越大，说明该节点的动态变化为其余节点动态变化的因，该节点所在状态变量或母线即为振荡源位置。
进一步，在振荡源位置确定后，从其所在状态变量或母线出发，保留因果强度系数最大的支路，即为该振荡模态的主要传播路径。
如图4所示，本文所提出的基于Copula传递熵的设备级和网络级宽频振荡传播路径分析及振荡源定位方法具体步骤如下：
图4 基于Copula传递熵的设备级和网络级宽频振荡传播路径分析及振荡源定位方法流程
Fig.4 Flow chart of device and network level wideband oscillations propagation path analysis and oscillations source localization method based on Copula transfer entropy
1）数据预处理
由于待分析的振荡数据包含系统中不同电压等级的各种电气量，数据大小的数量级相差很大，因此首先对振荡数据进行归一化处理。采用零均值归一化的方法，其公式为
（32）
式中，n为该时间序列的长度；
为该组数据的均值。由于宽频振荡的因果分析只需要信号中的振荡周期项，因此需要对振荡信号进行分解。采用经验模态分解方法[29]（Empirical Mode Decomposition，EMD），可以将原始信号分解成为一系列固有模态分量和一个残余分量，去除其中的残余分量。
2）基于Copula传递熵的因果分析
对于每个模态，计算各变量间的Copula传递熵，并确定因果传递方向。若传递熵＞，则此模态的因果传播方向是，因果强度系数为的数值；反之，若传递熵＜，则此模态的因果传播方向是，因果强度系数为的数值。
3）宽频振荡因果网络的构建
利用有向加权图构建宽频振荡因果网络，其中，网络各节点为状态变量或母线功率信号，各节点间的边权为对应变量间的因果强度系数。
4）振荡源位置的确定
在宽频振荡因果网络的基础上，根据式（31）计算网络中每个节点的出度，出度最大的节点被确定为振荡源位置。
5）振荡传播路径的确定
在振荡源位置确定后，从其所在状态变量或母线出发，每个节点保留因果强度系数最大的支路，即为该振荡模态的主要传播路径。
为了验证本文所提宽频振荡因果分析方法的准确性和可行性，在Matlab/Simulink中搭建直驱风机并网系统和含风电场的四机两区系统，分别对控制器内部各环节和复杂系统各母线之间的振荡因果关系进行分析说明。本节中Copula传递熵的计算采用开源Python工具包Copent。其中，嵌入维度的取值取决于时间序列的马尔可夫阶数[30]，在暂态稳定分析中，可以认为电力系统下一刻所处的状态只与系统现在的状态相关[31]，因此电力系统的时间序列皆为一阶马尔可夫过程，两变量的嵌入维度k和l的取值均为1。由于振荡信号呈现周期性，且不同信号之间存在一定的时延，因此计算不同信号的Copula传递熵时，滞后时间不是一个定值，而是通过选取Copula传递熵最大时对应的滞后时间来确定的。文献[32]已经证明，设变量X的滞后时间τ为0，则当Y的滞后时间t等于X和Y的实际时延δ时，两者之间的Copula传递熵最大，最能反映两变量间的因果关系，即
（33）
因此，Y的滞后时间t取值为两变量的实际时延δ。
Case 1采用的是直驱风机并网系统，其拓扑如图1所示。为了激发出该系统的振荡模态，当系统稳定运行时，调整网侧控制器电压外环参数，并在5 s时，在GSC 端口电压d轴分量utd上施加5 V的电压阶跃扰动，持续时间为0.2 s，得到14.8 Hz的次同步振荡模态，扰动清除后的GSC输出有功功率波形如图5所示。
选取扰动清除后1 s内的数据进行分析，经过EMD模态分解和归一化处理之后，部分状态变量之间的Copula传递熵计算结果见表2。
图5 GSC输出有功功率的振荡波形
Fig.5 Oscillation waveform of GSC output active power
表2 部分状态变量间的Copula传递熵
Tab.2 Copula transfer entropy between partial state variables
传递熵()传递熵() ①②0.5540.101 ②③2.3560.019 ③⑤1.0780.214 ④⑥1.1040.351 ⑤⑦1.2320.223 ⑥⑧1.2170.218 ⑦⑨1.0250.305 ⑧⑩1.1150.342
由表2可知，状态变量
作为原因，作为结果时的Copula传递熵明显大于作为原因，作为结果时的Copula传递熵，因此表2中的状态变量为因，状态变量为果。
接着，根据表中数据构建各状态变量间的加权有向图，如图6所示，并采用2.3节的出度公式计算得到图中各节点的出度，见表3。
图6 宽频振荡因果关系的加权有向图
Fig.6 Weighted directed graph of causal relationships in wideband oscillations
表3 各状态变量的出度Dout计算结果
Tab.3 Calculation results of Dout of each state variable
参数数值 ①ΔPe6.554 ②ΔUdc14.252 ③5.467 ④9.802 ⑤6.692 ⑥6.279 ⑦4.979 ⑧5.199
由表3可知，状态变量
的出度最大，因此为振荡源所在位置，其为因时所对应的传递函数为网侧控制器电压外环环节，与实际的振荡产生原因一致，由此验证了所提方法定位设备内部振荡源的有效性。
直驱风机内部振荡的传播路径如图7所示。从状态变量出发，如图7中红色虚线所示，系统存在多条振荡传播路径。保留因果强度系数最大的支路，图7中虚线圆圈内为该支路上的状态变量，红色实线即为直驱风机内部该振荡模态的主要传播路径。
图7 直驱风机内部振荡的传播路径
Fig.7 Propagation path of oscillations inside PMSG
从图7可以看出，该次同步振荡模态主要与直驱风电机的直流电容环节、网侧直流电压外环和电流内环有关，这些环节形成了该模态的内部传播路径。
含风电场的四机两区系统的拓扑结构如图8所示，共有12条母线，风电场W1由200台永磁同步直驱风机构成，总容量为300 MW，经由变压器从母线5连接至母线11。经过分析得到系统存在三个与风机相关的振荡模态，分别为14.5 Hz的次同步振荡、80 Hz的超同步振荡和111 Hz的中高频振荡。本节将针对不同频率的振荡（次/超同步、中高频）和不同机理的振荡（强迫振荡、风机与弱电网交互、发电机轴系振荡）分别进行验证。
图8 含风电场、不含轴系的四机两区模型
Fig.8 Four-machine two-area system model with wind farm and without generator shaft system
考虑到实际同步相量测量装置，本系统的信号采样频率设置为2 400 Hz，数据长度为2 400个采样点，从扰动注入时刻开始采样，采集与母线相连的各支路功率作为母线的振荡信号，共计1 s。
3.2.1 强迫振荡
系统稳定运行状态下，分别将14.5 Hz和80 Hz的正弦电压扰动源加在风机并网的PCC处，其幅值为5%的母线额定电压，大小为34.5 V，持续时间为1 s，采样得到各母线的功率振荡信号，对其进行模态提取与归一化处理，分别得到振荡频率为14.5 Hz和80 Hz的振荡信号，如图9和图10所示。
图9 预处理后的次同步振荡信号
Fig.9 Preprocessed subsynchronous oscillation signals
图10 预处理后的超同步振荡信号
Fig.10 Preprocessed supersynchronous oscillation signals
根据系统的拓扑关系，对于含不同支路的母线，根据本文所提方法，分别对与其相邻的各支路振荡信号进行Copula传递熵的计算。以次同步振荡模态为例，如图11所示，其中紫色实线箭头表示左侧母线为因，右侧母线为果，绿色虚线箭头表示右侧母线为因，左侧母线为果，图中数值为Copula传递熵的数值。
图11 次同步振荡分量Copula传递熵的计算结果示意图
Fig.11 Schematic diagram of the calculation results of Copula transfer entropy for subsynchronous oscillations
在此基础上，通过比较紫色实线箭头和绿色虚线箭头上数值的大小判断因果传播的方向，得到次同步振荡传播的因果路径图，采用同样的步骤得到超同步振荡传播的因果路径图，如图12所示。其中，红色实线箭头的方向代表次同步振荡传播的方向，蓝色虚线箭头的方向代表超同步振荡传播的方向，其上的Copula传递熵即为振荡传播的因果强度系数。
图12 次、超同步振荡模态的传播路径
Fig.12 Propagation path of subsynchronous and supersynchronous oscillation mode
由图12可知，振荡传播路径的起点皆为母线5，说明振荡源位于母线5处，由于此强迫振荡模态的扰动源施加在风电场出口，即为母线5处，因此验证了所提方法对强迫振荡源定位的准确性，且对不同频率的振荡模态具有同样的适用性。
在此基础上，利用式（14）分别计算次、超同步振荡模态各节点的振荡传播强度
，结果见表4、表5，其中，母线a为因，母线b为果。
由表4、表5可知，振荡传播强度与因果强度系数存在一致的趋势，因此验证了本文方法分析强迫振荡传播特性的准确性。
3.2.2 风机与弱电网交互引发的自然振荡
系统稳定运行状态下，在风机并网的PCC处施加一个幅值为5 V，持续时间为0.2 s的阶跃扰动引发振荡，采样得到各母线的自然振荡信号，对其进行模态提取与归一化处理，得到振荡频率为111 Hz的12组时序振荡信号，如图13所示。同样采用本文所提方法计算Copula传递熵，作出振荡传播的因果路径如图14所示。
表4 次同步模态的振荡传播强度和因果强度系数
Tab.4 Oscillation propagation intensity and causal intensity coefficient of subsynchronous oscillation mode
母线a母线b因果强度系数 母线5母线111.0082.380 母线11母线41.0052.297 母线11母线100.9952.043 母线11母线121.0542.385 母线7母线20.8560.966 母线7母线60.9731.503
表5 超同步模态的振荡传播强度和因果强度系数
Tab.5 Oscillation propagation intensity and causal intensity coefficient of supersynchronous oscillation mode
母线a母线b因果强度系数 母线5母线111.0962.104 母线11母线41.0131.197 母线11母线101.0641.711 母线11母线121.0881.840 母线7母线20.9950.962 母线7母线60.9130.942
图13 预处理后的自然振荡信号
Fig.13 Preprocessed natural oscillation signals
由图14可知，振荡源同样在母线5处，由于此振荡模态由风机引发，即为母线5处，因此验证了所提方法对自然振荡源定位的准确性。
在此基础上，同样计算得到各节点的振荡传播强度，结果见表6。
由表6可知，振荡传播强度与因果强度系数存在一致的趋势，因此验证了本文方法分析自然振荡传播特性的准确性。
图14 自然振荡模态的传播路径
Fig.14 Propagation path of natural oscillation mode
表6 自然振荡模态的振荡传播强度和因果强度系数
Tab.6 Oscillation propagation intensity and causal intensity coefficient of natural oscillation mode
母线a母线b因果强度系数 母线5母线111.1142.126 母线11母线41.0010.757 母线11母线101.1061.331 母线11母线121.0981.609 母线7母线21.1021.608 母线7母线60.9961.114
3.2.3 发电机轴系振荡
进一步，为了验证本方法在多个不同机理振荡源情况下的有效性，如图15所示，在四机两区系统的3号发电机上增加了轴系模型，从而获取次同步振荡模态[33]。轴系采用六质量块模型，参数设置选用IEEE次同步谐振第一标准模型中的默认设置[34]。
图15 含风电场和轴系的四机两区模型
Fig.15 Four-machine two-area system model with wind farm and generator shaft system
在此基础上，在风机并网的PCC处施加一个幅值为5 V，持续时间为0.2 s的阶跃扰动引发振荡，采样得到各母线的功率振荡信号，对其进行模态分解与归一化处理，得到频率为80 Hz的风机与弱电网交互引发的振荡模态和频率为12 Hz的发电机轴系引发的振荡模态，预处理后的时序振荡信号如图16、图17所示。
采用本文所提方法分别对两个振荡模态计算Copula传递熵，作出各自的振荡传播因果路径如图18所示。其中，红色实线箭头的方向代表风机与弱电网交互振荡传播的方向，蓝色虚线箭头的方向代表发电机轴系振荡传播的方向，其上的Copula传递熵即为振荡传播的因果强度系数。
图16 风机与弱电网交互引发的振荡模态波形
Fig.16 Oscillation mode waveforms caused by interaction between wind turbines and weak grid
图17 发电机轴系引发的振荡模态波形
Fig.17 Oscillation mode waveforms caused by generator shaft system
图18 多种机理的振荡模态传播路径
Fig.18 Propagation path of oscillation modes of different mechanisms
由图18可知，80 Hz的振荡源在5号母线，即风机出口处，12 Hz的次同步振荡源在3号母线处，与轴系模型所加位置一致。因此所提方法能够实现多个不同机理振荡源情况下的传播路径确定及振荡源定位。
在此基础上，对于两个模态分别计算得到各节点的振荡传播强度，结果见表7、表8。其中，模态1为风机与弱电网交互引发的振荡模态，模态2为发电机轴系引发的振荡模态。
表7 模态1的振荡传播强度和因果强度系数
Tab.7 Oscillation propagation intensity and causal intensity coefficient of oscillation mode 1
母线a母线b因果强度系数 母线5母线111.0030.647 母线11母线41.0721.254 母线11母线100.9760.759 母线11母线120.8100.678 母线7母线21.0070.984 母线7母线60.8130.619
表8 模态2的振荡传播强度和因果强度系数
Tab.8 Oscillation propagation intensity and causal intensity coefficient of oscillation mode 2
母线a母线b因果强度系数 母线3母线121.0021.075 母线11母线41.0851.210 母线11母线51.1101.373 母线11母线101.1051.206 母线7母线20.9930.676 母线7母线60.8420.298
由表7、表8可知，振荡传播强度与因果强度系数存在一致的趋势，因此验证了本文方法分析多个不同机理振荡传播特性的准确性。
通过以上分析，对于不同机理的振荡，由传递函数计算得到的振荡传播强度与由量测数据计算得到的因果强度系数存在一致的趋势，因此验证了所提方法分析宽频振荡传播特性的准确性。然而，振荡传播强度的计算需要系统模型和参数已知，且其只能反映振荡传播的大小，不能反映方向，因此无法确定振荡源。而本文所提方法基于数理统计，摆脱了对系统模型的依赖，且能够同时确定振荡源位置和振荡传播路径，对于宽频振荡的分析具有重要意义。
本文基于Copula传递熵提出一种设备级和网络级宽频振荡传播路径分析及振荡源定位方法，所提方法首先对振荡数据进行预处理，再利用Copula传递熵分别从设备级和网络级两个层面对振荡数据间的因果关系进行定量计算，最后通过构建因果网络确定宽频振荡的传播路径和振荡源的所在位置。与传统的振荡分析方法相比，所提方法不仅摆脱了对模型的依赖，不受样本数量的限制，且物理意义明确。直驱风机并网模型和含风电场的四机两区模型的仿真验证了所提方法的有效性，得到的主要结论如下：
1）本文采用Copula传递熵的方法对振荡数据间的因果关系进行分析，该方法无模型假设，对数据分布没有要求，且具有较为明确的物理意义。
2）本文所提的基于Copula传递熵的宽频振荡传播路径分析及振荡源定位方法，能够从设备级分析控制器内部各状态变量的因果关联性，从而确定宽频振荡在控制器内部的传播路径和振荡源位置，为研究振荡传播的内在机理提供支撑。
3）本文所提的基于Copula传递熵的宽频振荡传播路径分析及振荡源定位方法，能够从网络级分析电网中各节点振荡数据的因果传递性，从而确定宽频振荡在电网中的传播路径和振荡源位置，便于采取振荡抑制的针对性措施。
4）本文所提的基于Copula传递熵的设备级和网络级宽频振荡传播路径分析及振荡源定位方法，基于振荡传播的因果关系，未对振荡机理进行假设，因此适用于多种机理振荡的传播路径分析和振荡源定位。
参考文献
[1] 樊肖杰, 迟永宁, 马士聪, 等. 大规模海上风电接入电网关键技术与技术标准的研究及应用[J]. 电网技术, 2022, 46(8): 2859-2870.
Fan Xiaojie, Chi Yongning, Ma Shicong, et al. Research and application of key technologies and technical standards for large-scale offshore wind farms connecting to power grid[J]. Power System Technology, 2022, 46(8): 2859-2870.
[2] 饶宏, 周月宾, 李巍巍, 等. 柔性直流输电技术的工程应用和发展展望[J]. 电力系统自动化, 2023, 47(1): 1-11.
Rao Hong, Zhou Yuebin, Li Weiwei, et al. Engineering application and development prospect of VSC-HVDC transmission technology[J]. Automation of Electric Power Systems, 2023, 47(1): 1-11.
[3] Wang Xiongfei, Blaabjerg F. Harmonic stability in power electronic-based power systems: concept, modeling, and analysis[J]. IEEE Transactions on Smart Grid, 2018, 10(3): 2858-2870.
[4] 孙华东, 徐式蕴, 许涛, 等. 电力系统安全稳定性的定义与分类探析[J]. 中国电机工程学报, 2022, 42(21): 7796-7809.
Sun Huadong, Xu Shiyun, Xu Tao, et al. Research on definition and classification of power system security and stability[J]. Proceedings of the CSEE, 2022, 42(21): 7796-7809.
[5] Liu Huakun, Xie Xiaorong, He Jingbo, et al. Subsynchronous interaction between direct-drive PMSG based wind farms and weak AC networks[J]. IEEE Transactions on Power Systems, 2017, 32(6): 4708-4720.
[6] 黄碧月. 直驱风电并网的电力系统动态特性分析与次同步振荡研究[D]. 武汉: 华中科技大学, 2021.
Huang Biyue. Analysis on dynamic characteristics and study on subsynchronous oscillation of power system with PMSG integration[D]. Wuhan: Huazhong University of Science and Technology, 2021.
[7] 高本锋, 王义, 曾四鸣, 等. 直驱风电场并入弱交流电网的次同步分量通路及阻尼特性分析[J]. 中国电机工程学报, 2022, 42(14): 5089-5103.
Gao Benfeng, Wang Yi, Zeng Siming, et al. Analysis of sub-synchronous component path and damping characteristics of D-PMSG-based wind farm incorporated into weak AC grid[J]. Proceedings of the CSEE, 2022, 42(14): 5089-5103.
[8] 王利超, 于永军, 张明远, 等. 直驱风电机组阻抗建模及次同步振荡影响因素分析[J]. 电力工程技术, 2020, 39(1): 170-177.
Wang Lichao, Yu Yongjun, Zhang Mingyuan, et al. Impedance model and analysis of subsynchronous oscillation influence factors for grid-connected full-converter wind turbines[J]. Electric Power Engineering Technology, 2020, 39(1): 170-177.
[9] 马宁宁, 谢小荣, 贺静波, 等. 高比例新能源和电力电子设备电力系统的宽频振荡研究综述[J]. 中国电机工程学报, 2020, 40(15): 4720-4732.
Ma Ningning, Xie Xiaorong, He Jingbo, et al. Review of wide-band oscillation in renewable and power electronics highly integrated power systems[J]. Proceedings of the CSEE, 2020, 40(15): 4720-4732.
[10] Gu Kanghui, Wu Feng, Zhang Xiaoping. Sub-synchronous interactions in power systems with wind turbines: a review[J]. IET Renewable Power Generation, 2019, 13(1): 4-15.
[11] 朱茂林, 刘灏, 毕天姝. 考虑风电场量测相关性的双馈风力发电机鲁棒动态状态估计[J]. 电工技术学报, 2023, 38(3): 726-740.
Zhu Maolin, Liu Hao, Bi Tianshu. Robust dynamic state estimation of doubly-fed induction generator considering measurement correlation in wind farms[J]. Transactions of China Electrotechnical Society, 2023, 38(3): 726-740.
[12] Xie Xiaorong, Zhan Ying, Shair J, et al. Identifying the source of subsynchronous control interaction via wide-area monitoring of sub/super-synchronous power flows[J]. IEEE Transactions on Power Delivery, 2020, 35(5): 2177-2185.
[13] 陈磊, 王文婕, 王茂海, 等. 利用暂态能量流的次同步强迫振荡扰动源定位及阻尼评估[J]. 电力系统自动化, 2016, 40(19): 1-8.
Chen Lei, Wang Wenjie, Wang Maohai, et al. Disturbance source location of subsynchronous forced oscillation and damping evaluation using transient energy flow[J]. Automation of Electric Power Systems, 2016, 40(19): 1-8.
[14] 姜涛, 高浛, 李筱静, 等. 基于小波耗散能量谱的电力系统强迫振荡源定位[J]. 电工技术学报, 2023, 38(7): 1737-1750.
Jiang Tao, Gao Han, Li Xiaojing, et al. Forced oscillation source location in power system using wavelet dissipation energy spectrum[J]. Transactions of China Electrotechnical Society, 2023, 38(7): 1737-1750.
[15] 陈剑, 杜文娟, 王海风. 基于对抗式迁移学习的含柔性高压直流输电的风电系统次同步振荡源定位[J]. 电工技术学报, 2021, 36(22): 4703-4715.
Chen Jian, Du Wenjuan, Wang Haifeng. Location method of subsynchronous oscillation source in wind power system with VSC-HVDC based on adversarial transfer learning[J]. Transactions of China Electro-technical Society, 2021, 36(22): 4703-4715.
[16] Ping Zuowei, Li Xiuting, He Wei, et al. Sparse learning of network-reduced models for locating low frequency oscillations in power systems[J]. Applied Energy, 2020, 262: 114541.
[17] 谢小荣, 刘华坤, 贺静波, 等. 新能源发电并网系统的小信号阻抗/导纳网络建模方法[J]. 电力系统自动化, 2017, 41(12): 26-32.
Xie Xiaorong, Liu Huakun, He Jingbo, et al. Small-signal impedance/admittance network modeling for grid-connected renewable energy generation systems[J]. Automation of Electric Power Systems, 2017, 41(12): 26-32.
[18] Gong Hong, Yang Dongsheng, Wang Xiongfei. Impact analysis and mitigation of synchronization dynamics for DQ impedance measurement[J]. IEEE Transactions on Power Electronics, 2019, 34(9): 8797-8807.
[19] Maslennikov S, Litvinov E. ISO new England experience in locating the source of oscillations online[J]. IEEE Transactions on Power Systems, 2021, 36(1): 495-503.
[20] 周一辰, 孙佳辉, 王书祥, 等. 基于高维模型表达方法的新能源电力系统小干扰失稳风险评估[J]. 电力系统自动化, 2022, 46(14): 73-82.
Zhou Yichen, Sun Jiahui, Wang Shuxiang, et al. Risk assessment of small-signal instability for renewable power system based on high-dimensional model representation method[J]. Automation of Electric Power Systems, 2022, 46(14): 73-82.
[21] 徐衍会, 刘慧, 成蕴丹. 次同步振荡在交直流电网中传播的关键影响因素[J]. 现代电力, 2023(4): 1-11.
Xu Yanhui, Liu Hui, Cheng Yundan. Key Influencing factors on propagation of sub-synchronous oscillations in AC and DC grids[J]. Modern Electric Power, 2023(4): 1-11.
[22] 张涌新. 双馈风电并网系统宽频振荡机理及振荡源定位研究[D]. 北京: 华北电力大学, 2020.
Zhang Yongxin. Research on the mechanism of wide-range-frequency oscillation and the location of oscillation source in the grid connected system of doubly fed wind power[D]. Beijing: North China Electric Power University, 2020.
[23] Ma Jian, Sun Zengqi. Mutual information is copula entropy[J]. Tsinghua Science & Technology, 2011, 16(1): 51-54.
[24] Ma Jian. Estimating transfer entropy via copula entropy[EB/OL]. arXiv, 2019: 1910.04375. https:// arxiv.org/abs/1910.04375.pdf.
[25] 唐锦, 张书怡, 吴秋伟, 等. 基于Copula函数与等概率逆变换的风电出力场景生成方法[J]. 电力工程技术, 2021, 40(6): 86-94.
Tang Jin, Zhang Shuyi, Wu Qiuwei, et al. Wind power output scenario generation method based on Copula function and equal probability inverse transformation[J]. Electric Power Engineering Technology, 2021, 40(6): 86-94.
[26] 李立, 张青蕾, 王康, 等. 风电并网系统次同步振荡传播路径频域定位方法[J]. 华北电力大学学报(自然科学版), 2023, 50(2): 54-62.
Li Li, Zhang Qinglei, Wang Kang, et al. Frequency domain location method of subsynchronous oscillation propagation path in wind power integrated system[J]. Journal of North China Electric Power University (Natural Science Edition), 2023, 50(2): 54-62.
[27] Granger C W J. Investigating causal relations by econometric models and cross-spectral methods[M]// Essays in Econometrics vol II: Collected Papers of Clive. New York: Cambridge University Press, 2008: 31-47.
[28] 耿雪青, 佘青山, 张启忠, 等. 基于Copula的多变量运动想象脑电信号因果分析方法[J]. 航天医学与医学工程, 2018, 31(1): 49-56.
Geng Xueqing, She Qingshan, Zhang Qizhong, et al. Multivariate causality analysis method of motor imagery EEG signals based on copula[J]. Space Medicine & Medical Engineering, 2018, 31(1): 49-56.
[29] 姜涛, 刘博涵, 李雪, 等. 基于自适应投影多元经验模态分解的电力系统强迫振荡源定位[J]. 电工技术学报, 2023, 38(13): 3527-3538.
Jiang Tao, Liu Bohan, Li Xue, et al. Forced oscillation location in power systems using adaptive projection intrinsically transformed multiple empirical mode decomposition[J]. Transactions of China Electrote-chnical Society, 2023, 38(13): 3527-3538.
[30] 王洪安. 上肢运动的肌间耦合分析方法研究[D]. 杭州: 杭州电子科技大学, 2021.
Wang Hongshan. Research on intermuscular coupling analysis method of upper limb movement[D]. Hangzhou: Hangzhou Dianzi University, 2021.
[31] 叶圣永, 王晓茹, 周曙, 等. 基于马尔可夫链蒙特卡罗方法的电力系统暂态稳定概率评估[J]. 电工技术学报, 2012, 27(6): 168-174.
Ye Shengyong, Wang Xiaoru, Zhou Shu, et al. Power system probabilistic transient stability assessment based on Markov chain Monte Carlo method[J]. Transactions of China Electrotechnical Society, 2012, 27(6): 168-174.
[32] Wibral M, Pampu N, Priesemann V, et al. Measuring information-transfer delays[J]. PLoS One, 2013, 8(2): e55809.
[33] 高本锋, 董涵枭, 卢亚军, 等. 直驱风电场并网对直流输电引起的火电机组轴系扭振影响机理分析[J].电工技术学报, 2024, 39(7): 1971-1984.
Gao Benfeng, Dong Hanxiao, Lu Yajun, et al. Mechanism analysis of the influence of direct drive wind farm integration on SSTI of thermal generator caused by LCC-HVDC[J]. Transactions of China Electrotechnical Society, 2024, 39(7): 1971-1984.
[34] 吴楠, 李尚盛, 查晓明. 改进的同步电机阻尼绕组电流观测器[J]. 电力系统保护与控制, 2015, 43(1): 67-72.
Wu Nan, Li Shangsheng, Zha Xiaoming. Improved synchronous generator damper current observer[J]. Power System Protection and Control, 2015, 43(1): 67-72.
Abstract Wideband oscillations in the power electronic power system seriously endangers the safe and stable operation of the power grid. Timely determination of the propagation path and source location of wideband oscillations is crucial for suppressing oscillations. Model based analysis methods have limitations such as unknown parameters and high complexity, restaining their wide application in practical power systems. With the development of information theory, the method based on mathematical statistics becomes a promising method for the analysis of wideband oscillations. These methods focus on the correlation between data items and discover the inherent laws of wideband oscillations. Considering that the data of wideband oscillations in the new-generation power system is a multivariable strongly nonlinear and causal time series, it is worth conducting in-depth research on how to analyze the propagation characteristics of wideband oscillations and achieve oscillations source localization from the perspective of causal relationships based on mathematical statistics. Therefore, this article proposes a device and network level wideband oscillations propagation path analysis and oscillations source localization method based on Copula transfer entropy from the perspective of causality.
Firstly, the zero-mean normalization method is used to preprocess the analyzed oscillations data. Secondly, the original signal can be decomposed into a series of intrinsic mode components and a residual component using empirical mode decomposition (EMD). The residual component can be removed as it is independent of oscillations. Thirdly, for each mode, calculate the Copula transfer entropy between each variable to determine the direction of causal transfer. Then define the numerical value of Copula transfer entropy as the causal strength coefficient. Fourthly, a wideband oscillations causal network is constructed using directed weighted graphs, where each node in the network is a state variable or bus power signal, and the edge weights between each node are the causal strength coefficients between the corresponding variables. Fifthly, calculate the out-degree of each node in the network, and the node with the highest out-degree is determined as the position of the oscillations source. Finally, starting from the state variable or bus where the oscillations source is located, retain the branch with the highest causal strength coefficient, which is the main propagation path of the oscillation mode. At device level, the causal relationship between state variables within the controller can be revealed, thereby determining the key state variable and oscillations propagation path. At network level, the causal relationship between the power of each bus in the power grid can be revealed, and then the bus where the oscillations source is located and the oscillations propagation path can be determined.
On one hand, a simulation example of PMSG connected to grid is used to analyze the causal relationship of oscillations in various links within the controller. The propagation path diagram of the internal oscillations of the PMSG is obtained using the method proposed in this article. The sub synchronous oscillation mode of the system is mainly related to the DC capacitor link, outer voltage loop, and inner current loop of the PMSG. These links form the internal propagation path of this oscillation mode.
On the other hand, a simulation example of a four-machine two-area system containing a wind farm is used to analyze the causal relationship of oscillations between the busbars of a complex system. Simulation examples in forced oscillations, interaction between wind turbines and weak grid, and generator shaft system oscillations demonstrate that the proposed method is suitable for various mechanisms of wideband oscillations, and can achieve wideband oscillations propagation path analysis and oscillations source localization from the network level.
The following conclusions can be drawn from the simulation analysis: (1) The method based on Copula transfer entropy avoids establishing a detailed model, and its simulation results are consistent with the oscillations propagation path obtained from theoretical derivation. (2) The method is suitable for various mechanisms of wideband oscillations such as forced oscillations, interaction between wind turbines and weak grid, and generator shaft system oscillations. (3) The method can simultaneously achieve wideband oscillations propagation path analysis and oscillations source localization at device and network level.
Keywords：Wideband oscillations, oscillations propagation, oscillations source localization, Copula transfer entropy, causal analysis
中图分类号：TM712
DOI: 10.19595/j.cnki.1000-6753.tces.230873
国家自然科学基金资助项目（52377084, 52007028）。
收稿日期 2023-06-07
改稿日期 2023-10-17
冯 双 女，1990年生，副教授，博士生导师，研究方向为电力电子化电力系统分析与控制、人工智能技术在电力系统中的应用。
E-mail：sfeng@seu.edu.cn（通信作者）
杨 浩 男，1998年生，硕士研究生，研究方向为基于量测数据的电力系统宽频振荡分析。
E-mail：yang_hao_@seu.edu.cn
（编辑 赫 蕾）