基于VMD和主客观赋权的次同步振荡能量主导因素辨识

徐衍会, 刘慧,成蕴丹,孙冠群,蔡德福,王尔玺

(1. 华北电力大学电气与电子工程学院,北京市102206;2. 国网湖北省电力有限公司电力科学研究院,武汉市 430077)

0 引 言

随着风电等可再生能源场站在电网的广泛应用,含高比例新能源电网将成为电力系统的典型形态。风电并网系统次同步振荡事件在世界范围内广泛发生,对电力系统安全稳定运行造成了严重威胁[1-4]。在次同步振荡特性分析中,现有研究大多采用建立状态空间模型求解特征根[5-7],或建立阻抗模型分析等效阻抗特性[8-10]。但状态空间法及阻抗模型都需要对所研究系统进行详细建模,且需要获取各元件相应的参数,同时系统模型阶数较高,计算耗时较长,多用于离线分析。宽频量测技术的发展,使数据驱动的次同步振荡分析成为可能[11-13]。因此,有学者基于电力系统量测数据提出了基于端口能量的能量函数分析方法[14-15],其构造过程无需计及风机内部的动态过程,有效避免对电力电子环节复杂特性的建模过程,充分利用风机端口可量测电气量,降低了能量函数的构造难度。因此,本文从风电场并网系统量测数据出发,利用端口能量的能量拟合函数表征次同步振荡能量特性,并对能量特性的影响因素建立评估模型辨识其主导因素。

基于量测数据的次同步振荡分析方法有Prony法[16-17]、快速傅里叶分解法[18-19]以及模态分解法等。Prony法采用不同振幅、相位、衰减因子以及频率的指数函数拟合时域信号,然而Prony法拟合结果阶数较高,计算效率较低,同时计算精度不高。快速傅里叶分解法将时域信号转换至频域,提取整数次基波分量,但由于实际电力系统中难以做到精准同步采样,易产生栅栏效应和频谱泄露等问题[20-21]。模态分解法可以分为小波变换、希尔伯特黄变换(Hilbert-Huang transform, HHT)两类[19]。小波变换在傅里叶变换的基础上,用小波函数取代三角函数,解决了傅里叶变换中窗口大小不能随频率变化的困难,但由于小波变换中引入的小波基没有标准的选取原则,不合适的小波基将会产生较大误差。HHT模态分解法由经验模态分解(empirical mode decomposition, EMD)和希尔伯特变换(Hilbert transform, HT)两部分构成[22]。HHT首先将时域信号进行EMD分解,再将分解后的本征模态函数(intrinsic mode fuction, IMF)分量进行HT变换,适用于分析复杂非线性时域信号的时频属性。然而EMD在信号分解过程中会引入虚假分量,产生模态混叠问题[23],使得IMF分量失去单一尺度特征。因此,本文引入变分模态分解(variational mode decomposition, VMD)方法对时域信号进行模态分解,将原始信号分解为一系列不同频率的模态分量,相比于经验模态分解方法解决了模态混叠弊端[24],分解得到准确的次同步模态分量用以提取次同步振荡能量函数中的能量特性。

能量函数法从能量的角度对电力系统振荡进行了分析研究,取得了良好的效果[25-27]。文献[28]首先将能量流法推广至次同步振荡分析,验证了暂态能量流对次同步振荡的适用性。文献[29-30]通过建立状态空间模型,分析了各因素对能流功率的影响。但实际复杂电力系统建模困难,同时不同运行工况、不同控制策略、不同外部环境下新能源发电机组次同步振荡特性具有很大差异,难以建立各结构形态下的振荡分析模型。因此,为研究不同场景下次同步振荡能量特性,本文从风电场运行时序数据出发实现模型驱动向数据驱动的转变,建立各影响因素对于能量特性的评估模型,进而辨识次同步振荡能量特性的主导因素。

基于上述分析,本文提出一种数据驱动的次同步振荡能量特性提取及影响因素辨识方法。首先,基于VMD算法避免经验模态分解方法中的模态混叠问题,提取准确的次同步模态时域信号;其次,推导次同步振荡下端口能量的能量函数表达式,利用次同步模态分量进行能量函数计算;最后利用主客观赋权方法建立能量特性的影响因素评估模型,并辨识影响次同步振荡能量特性的主导因素。通过理想信号以及仿真模型试验,本文所提方法能够较为准确地提取模态信息,进而辨识能量特性的主导因素,有效解决实际工程中次同步振荡问题。

1 变分模态分解

1.1 变分模态分解理论

变分模态分解是信号处理中一种自适应的完全非递归的模态变分和信号分析方法。VMD基于经典维纳滤波理论求解变分问题,得到各中心频率及带宽限制,进而得到各中心频率在频域中对应的有效成分。

VMD可以在搜索求解的迭代过程中自适应对各模态中心频率及其对应带宽进行匹配,使得每个模态的估计带宽值最小,实现信号频谱剖分及不同分量分离,最终获得变分问题的最优解。同时,VMD方法可以克服模态分量混叠、过包络及欠包络等经典模态分解中的问题,可以从非平稳性序列中分解得到各频率分量尺度下相对平稳的子序列,相较于经典模态分解方法更适合处理复杂度高、非线性强的时间序列,有助于分析非线性电力系统中波动性较强的电气量。

VMD的模态分解过程可以概括为变分问题的构造以及变分问题的求解。求解过程中要求各模态中心频率对应带宽之和最小,同时要求所有模态分量之和为原始信号。VMD相较于一般模态分解方法重新定义了更严格的有限带宽的IMF,并认为每个IMF分量都是调幅调频信号,其具有特定中心频率、有限带宽,如式(1)所示。

uk(t)=Ak(t)cos[φk(t)]

(1)

式中:φk(t)为信号uk(t)对应相位;相位导数dφk(t)/dt为信号uk(t)的瞬时频率;Ak(t)为信号uk(t)的包络线幅值。

泛函是一种从任意向量空间到标量的映射,求解泛函极值的问题称为变分问题。求解VMD变分问题的极值要求各模态分量中心频率的带宽之和最小,约束条件为所有模态分量之和等于原始信号。VMD约束变分模型如下:

(2)

式中:uk为将原始信号f(t)分解为K个中心频率不同的IMF分量,k=1,2,…,K;wk为各IMF模态下的中心频率;*为卷积运算符;∂t为对函数求时间的偏导数;δ(t)为狄拉克函数。

为解决上述约束变分模型的最优化问题,需要利用二次惩罚项和拉格朗日乘子法,将约束变分问题转变为非约束变分问题,引入增广Lagrangian函数,得到拓展的Lagrange方程,如式(3)所示。

(3)

式中:τ为惩罚因子,可以降低高斯噪声的干扰;λ(t)为拉格朗日乘子;<·>为内积运算符号。

为求解非约束优化问题,引入交替方向乘子(alternating direction method of multipliers, ADMM)迭代算法结合Parseval/Plancherel、傅里叶等距变换等方法优化得出各模态分量及其中心频率。进而迭代搜寻增广Lagrange方程的鞍点,得到交替寻优迭代的更新泛函uk、wk和λ的表达式如下:

(4)

(5)

(6)

(7)

VMD算法的具体过程可以概括为:

2)n=n+1;

1.2 VMD算法辨识理想信号

为验证VMD算法的数据特征提取能力及参数优化的有效性,构造模态差距较大的理想信号,并加入14.09 dB高斯白噪声对VMD模态分解性能进行验证。理想信号如式(8)所示。

f(t)=207e-0.48tcos(2π×25.4t-π/6)+
20e0.29tcos(2π×15.5t+π/4)+
396e-0.32tcos(2π×41.2t)

(8)

对理想信号f(t)进行VMD分解,设置模态分量数K=3,VMD分解后各分量对应时域波形以及频谱图如图1所示。

图1 理想信号经VMD分解后的IMF分量Fig.1 The IMF components of the ideal signal decomposed by VMD

为了进一步验证VMD算法分解结果的准确性,分别用EMD、集合经验模态分解(ensemble empirical mode decomposition,EEMD)方法对信号分解后进行主要振荡模态频率提取,如表1所示。

表1 不同模态分解方法对主要振荡模态提取结果Table 1 Frequency extraction results from different modal decompo-sition methods

从表1中可以看出VMD分解结果能够很好处理模态混叠的信号,较为完整地提取所需频率(15.532 Hz)的时域波形,同时基本完全剥离了非主导模态(25.398 Hz、41.197 Hz)时域分量。为对模态分解结果进行定量对比,引入信噪比(signal to noise ratio,SNR)、均方根误差(root mean square error,RMSE)等评价指标。

(9)

(10)

对比不同分解方法精度评价指标可以得出VMD方法模态分解后信噪比得到明显增加,均方根误差得到明显下降,证明VMD更适用于处理多模态时域信号,并能较为准确地提取所需频率的时域波形。

2 次同步振荡能量特性提取

为研究风机次同步振荡能量特性,需要构造能够表征风机次同步振荡的能量函数。由于实际电力系统中可以量测得到端口电气量,故而本文采用基于端口能量的能量函数定义。从节点i经支路Lij流向节点j的能量函数如式(11)所示。

(11)

双馈风机(doubly fed induction generator, DFIG)并网正常运行时,双馈风机端口电压三相对称。若系统发生次同步振荡,则电压电流中会产生次同步频率的扰动分量。以a相为例,端口电压电流可以表示为:

ua=U0cos(ω0t+δ0)+Uercos(ωert+δer)

(12)

ia=I0cos(ω0t+δ0)+Iercos(ωert+θer)

(13)

式中:U0和I0分别为工频电压和工频电流有效值;Uer和Ier分别为次同步频率电压和次同步频率电流有效值;ω0=2πf0、ωer分别为工频角频率和次同步频率角频率,f0=50 Hz;δ0、θ0分别为工频电压和工频电流的初相位;δer、θer分别为次同步频率电压和次同步频率电流的初相位。

将三相端口电压经Park变换转换至dq坐标系后,端口电压电流可以表示为:

(14)

(15)

式中:ω-=ω0-ωer;ud、uq为dq坐标系下端口电压;id、iq为dq坐标系下端口电流。

将式(14)、(15)代入式(11)可得节点i能量函数:

(16)

同理,对于直驱风电场并入弱电网系统,也可以推导出端口能量表达式,与双馈风机不同,直驱风机(permanent magnet sgnchronous wind turbine generator, PMSG)端口含有超同步分量,其表达式为:

(17)

式中:Usub和Isub分别为次同步频率电压和次同步频率电流有效值;Usup、Isup分别为超同步频率电压和超同步频率电流有效值;δsub、θsub分别为次同步电压和次同步电流的初相位;δsup、θsup分别为超同步频率电压和超同步频率电流的初相位。

从式(16)和式(17)可以看出系统发生次同步振荡后,能量函数由元件能量变化项和元件消耗或产生的能量项两部分构成。其中元件能量变化项为周期分量;元件消耗或产生的能量项为一次函数形式的非周期分量,其一次项系数可以代表元件的阻尼特性,称为能流功率。

能量函数可由实际电力系统中采集的电气量计算得到,然而实际量测信号中包含噪声分量等干扰信号,当干扰信号分量作用较大时,能量函数法计算得到的误差较大,对风机次同步振荡能量函数的特征提取造成困难。此外,从式(16)、(17)可以看出能流功率仅与次、超同步模态电气量相关。因此,将采集的电气量进行模态分解处理,提取次同步模态时域波形,可以消除其他各频段分量对能流功率计算的影响,实现对次同步振荡能量特性的准确提取。

在模态分解方法中,变分模态分解相较于经典模态分解方法可以避免模态混叠现象,对采样和噪声具有更强的鲁棒性,可以得到更准确的次同步频率电压和电流波形。因此本文引入变分模态分解方法提取时域波形的次同步模态时域波形,获得更准确的能流功率,并且避免实际系统中采用的实时采集电气量进行大量能量函数的积分计算,大幅提升能量函数的计算效率,进而为风机并网次同步振荡能量特性影响因素辨识奠定基础。

3 次同步振荡能量特性影响因素辨识

文献[30]考虑风速变化对风电机组能流功率的影响建立机理模型,然而,影响能量特性的因素众多,能流功率不仅与风速相关,而且也与风机参数、系统当前运行状况等相关,各影响因素对能流功率的影响程度尚不相同,故需建立评估模型,筛选与能量特性强相关的影响因素。

基于上述分析,本文采用主客观赋权法对次同步振荡能量特性影响因素进行综合评估,其中,各因素的主观权重由层次分析法(analytic hierar-chy process, AHP)获得,客观权重由多元线性回归分析获得。

3.1 基于层次分析法的主观权重

本文采用基于比例标度构造法的AHP判断矩阵来进行专家打分,采用该方法构造的判断矩阵进行分析不需要进行一致性检验,从而提升AHP方法的实用性以及可靠性,AHP具体步骤如下:

1)根据专家意见得到m个评价指标的重要性排序,表示为:x1≥x2≥…≥xm,其中xj表示主观重要性排名第j个的评价指标。

2)根据AHP的比例标度得到相邻两个指标之间的标度值tj,tj取值如表3所示。

表3 AHP方法的标度值Table 3 Scale values of AHP method

3)基于标度值的传递,得到m个影响因素的主观判断矩阵G,如式(18)所示:

(18)

gjk为矩阵G中第j行第k列元素,j,k=1,2, … ,m,表示指标j相较于指标k的重要程度,矩阵G中的元素满足:gjj=1,gkj=1/gjk,gjk=gjwgwk, 基于矩阵G中的元素,第j个影响因素的主观权重为ω′j计算公式如式(19)所示:

(19)

3.2 基于多元线性回归法的客观权重

通过变量回归拟合得到各影响因素与能量函数的关联关系,故可以使用多元线性回归方法量化表征各因素与能量特性的线性相关性以获取客观权重,多元线性回归模型如式(20)所示:

Y=β0+β1y1+β2y2+…+βmym+ε

(20)

式中:Y为代表能量特性的因变量;y1,y2, …,ym为m个与能量特性相关的自变量;β0为常数项系数,β1,β2,…,βm为自变量对应的回归系数,可以采用最小二乘法求解其拟合值;ε为均值为0的随机变量,用来代表随机因素对因变量的影响。对回归系数β1,β2, …,βm取绝对值后进行归一化处理,即可获得第j个影响因素的客观权重ωj″,如式(21)所示:

(21)

3.3 主客观赋权法综合权重系数计算

若设置h个系统运行工况,m个系统能量特性影响因素,则样本矩阵R′如式(22)所示:

(22)

式中:Rij表示第i个工况下第j个影响因素的值,考虑到不同影响因素的量纲不同,故首先对样本矩阵各元素做标准化处理,处理方法如式(23):

(23)

式中:rij为标准化处理后的样本值,取值范围为[0,1],Rjmax=max{R1j,R2j,…,Rhj},Rjmin=min{R1j,R2j,…,Rhj},则可构成标准化矩阵R,如式(24)所示:

(24)

采用线性加权的方法综合考虑主观权重和客观权重,得到综合权重表达式为:

ω=αω′+βω″

(25)

式中:ω为综合权重系数;ω′为主观权重;ω″为客观权重;α为主观权重系数;β为客观权重系数,α和β满足:α+β=1且α,β≥0。

为了充分体现客观权重和主观权重且不偏颇于某一方,构造如式(26)所示优化模型:

(26)

(27)

根据式(27)得到主客观影响因素权重系数,代入式(24)得到最终的综合权重{ω1,ω2, … ,ωm},本文所研究系统次同步振荡能量特性主导因素辨识方法总体流程如图2所示,通过综合权重可以定量评价各影响因素对系统能量特性的影响程度。

图2 次同步振荡能量主导因素辨识方法流程图Fig.2 Flow chart of the identification of dominant factors of sub-synchronous oscillation energy

4 算例研究

本文基于VMD进行次同步振荡能量特性提取,结合主客观赋权方法进行能量函数拟合并辨识其主导因素。为验证上述方法有效性,在电磁暂态仿真平台PSCAD/EMTDC中搭建双馈风电场经串补并入无穷大电网模型以及直驱风电场并入弱电网模型进行时域仿真,产生能量特性影响因素辨识研究数据集。

4.1 双馈风电场经串补并入无穷大电网

4.1.1 VMD模态分解辨识

在PSCAD/EMTDC中搭建如图3所示的双馈风机并入含串联补偿系统,等值双馈风机由30台相同的额定功率为1.5 MW的双馈风机组成,假设风电场所有风电机组的运行状态相同,使用单台双馈风机代表所有风电机组,其等值参数可以由单台风电机组参数获得,单台风机参数如表4所示。双馈风机分别由0.69 kV/35 kV、35 kV/220 kV两个升压变压器经串补线路接入交流电网。

表4 单台双馈风机具体参数Table 4 Specific parameters of single DFIG

图3 双馈风电场并入串补电网系统等值模型Fig.3 Equivalence model of a doubly-fed wind farm integrated into a series-complementary grid system

对如图3所示系统进行仿真,获取数据,设置10 s时串补电容并入电网引发次同步振荡。为获取仿真数据集进行波形辨识,改变交流线路串补电容大小、风电场风速、风电场控制参数、风电场额定容量等构建数据集,各变量变化范围如表5所示。

表5 双馈风机变量设置Table 5 Variable settings for DFIG

三相电压经VMD分解前后信号如图4所示,得到IMF分量及其对应频谱,如图5所示。

图4 次同步振荡下三相电压时域波形Fig.4 Three-phase voltage time domain waveform under subsyn-chronous oscillation

图5 VMD分解后IMF分量Fig.5 IMF components after VMD decomposition

对数据集中每组的电压电流量进行VMD模态分解,得到次同步模态时域波形,由式(16)拟合计算得到每组数据对应的能流功率,用来表征双馈风机并网系统次同步振荡的能量特性,所用计算机处理器为11th Gen Intel(R) Core(TM) i7-11800H @ 2.30 GHz,本文所用数据的采样频率为200 Hz,时间段为100 ms,通过对时序数据进行测试,次同步振荡能量特性辨识时间为0.131 s,如果使用性能更好的计算机将会进一步缩短能量特性提取所需时间。

4.1.2 影响次同步振荡能量特性变量选择(双馈风机)

以双馈风机并入串补电网引发次同步振荡为例,振荡受到电网外部运行条件以及风机自身运行状态的影响。电网外部运行条件中,串补电容的大小对次同步振荡具有较大影响,其电容值越大,发生次同步振荡风险越大;同时风电场出口有功功率、无功功率以及电压可以反映次同步振荡波动情况。风机自身运行状态中,风速对次同步振荡具有较大影响,风速越小,发生次同步振荡风险越大;其次RSC控制参数中转子侧电流内环比例系数对振荡影响较大[31],比例系数越大,发生次同步振荡风险越大;此外发电机转速也在一定程度上反映了次同步振荡波动程度。因此,选择串补电容值x1、风速x2、风电场出口侧有功功率x3,出口侧无功功率x4、出口侧电压x5、发电机转速x6、RSC电流内环比例系数x7作为次同步振荡能量特性的自变量。按照运行工况、控制参数、外部环境分类,结果如表6所示。同时选择能够表征风机并网次同步振荡能量特性的能流功率y作为因变量进行分析,并探究各自变量与能量特性之间的相关关系。

表6 能量特性的自变量分类(双馈风机)Table 6 Classification of independent variables for energy characteristics(DFIG)

4.1.3 主客观赋权法分析(双馈风机)

首先,基于AHP计算主观权重,根据专家意见对7个影响因素进行重要性排序:x1≥x7≥x3≥x2≥x4≥x5≥x6,设置标度值t1=1.2,t2=1.4,t3=1.2,t4=1.4,t5=1.2,t6=1.2,得到m个影响因素的主观判断矩阵G,如表7所示,则得到的主观权重ω′如表8所示。

表7 各因素的主观判断矩阵GTable 7 Subjective judgment matrix G for each factor

表8 各因素的主观权重ω′Table 8 The subjective weights ω′ of the factors

其次,建立多元回归模型获取客观权重ω″,以各影响因素x1～x7作为自变量,能量函数的能流功率作为因变量进行回归分析。分别进行多元线性回归和逐步回归,得到回归模型,如式(28)所示,回归模型显著性检验量如表9所示。

表9 回归模型显著性指标Table 9 Regression model significance metric

y=-0.233 1+0.616 2x1+1.379 4x2-1.871 3x3+
0.024 0x4-0.204 6x5+0.209 8x6-0.612 5x7

(28)

表中F统计量对应的p值1.637×10-42远小于显著性水平0.05对应的p值,因此认为能流功率与各影响因素之间存在显著的线性回归关系。由回归系数可计算得到客观权重ω″,各影响因素的客观权重如表10所示。

表10 各因素的客观权重ω″Table 10 The subjective weights ω″ of the factors

根据式(27)计算可得到综合权重系数α=0.407 3,β=0.592 7,故可由式(24)最终计算得到综合权重ω,如表11所示。由表11可知,所研究的7个影响因素中,串补电容值x1、风速x2、风电场出口侧有功功率x3、RSC电流内环比例系数x7与能流功率有较强的相关性,而风电场出口侧无功功率x4、风电场出口侧电压x5、发电机转速x6与能流功率关系较弱,故在对能流功率的进一步研究中,应当着重分析辨识出的4个强相关影响因素,以更有效地在工程实际中解决次同步振荡问题。

表11 各因素的综合权重ωTable 11 The combined weights ω of the factors

4.2 直驱风电场并入弱电网

在PSCAD/EMTDC中搭建如图6所示的直驱风电场并入弱电网系统,等值直驱风机由20台相同的额定功率为5 MW的直驱风机组成,假设风电场所有风电机组的运行状态相同,使用单台直驱风机代表所有风电机组,其等值参数可以由单台风电机组参数获得,具体参数如表12所示。

表12 单台直驱风机具体参数Table 12 Specific parameters of single PMSG

图6 直驱风电场并入弱电网系统等值模型Fig.6 Equivalent model of direct driven wind farm connected to weak grid system

4.2.1 影响次同步振荡能量特性变量选择(直驱风机)

选择短路比x1、风速x2、GSC电流内环比例系数x3,锁相环比例系数x4、风电场出口侧有功功率x5、风电场出口侧无功功率x6、风电场出口侧电压x7作为该系统次同步振荡能量特性的自变量,以能流功率y作为因变量进行变量拟合,并探究各自变量与能量特性之间的相关关系,按照运行工况、控制参数、外部环境分类,结果如表13所示。

表13 能量特性的自变量分类(直驱风机)Table 13 Classification of independent variables for energy characteristics(PMSG)

4.2.2 主客观赋权法分析(直驱风机)

通过本文所建的相关性评估模型,分别获取主观权重ω′和客观权重ω″,然后结合权重系数α=0.117 5,β=0.882 5,由式(25)最终计算得到综合权重ω,各权重数值如表14所示。

表14 各影响因素的权重(直驱风机)Table 14 The combined weights of the factors(PMSG)

由表14可知,在直驱风电场并入弱电网系统,短路比x1、GSC电流内环比例系数x3、锁相环比例系数x4、风速x2与能流功率有较强的相关性,而风电场出口侧有功功率x5、风电场出口侧无功功率x6、风电场出口侧电压x7与能流功率关系较弱。

5 结 论

本文提出了一种基于VMD和主客观赋权方法的次同步振荡能量特性提取及主导因素辨识方法,并通过PSCAD搭建的双馈风机并网系统进行验证,结论如下:

1)本文基于VMD对系统可量测电气量进行模态分解,相较于经典模态分解方法避免了模态混叠的问题,进而准确提取出可以表征次同步振荡能量特性的能流功率。

2)通过基于层次分析和多元线性回归的主客观赋权法,优化得到综合权重,进而辨识出次同步振荡能量特性的主导因素。分析表明,在双馈风机经串补并网系统中,串补度、风速、风电场出口有功功率以及RSC电流内环比例系数是影响次同步振荡能量特性的主导因素;在直驱风电场并入弱电网系统中,短路比、GSC电流内环比例系数、锁相环比例系数以及风速是影响次同步振荡能量特性的主导因素。

3)本文所提方法不依赖于系统模型,不拘束于特定的系统运行方式和机组控制策略,仅利用实际电力系统中可量测的电气量,可以获得次同步振荡能量特性的主导因素,对基于能量的次同步振荡分析与抑制起到支撑作用。