混合随机激励下的多机电力系统稳定性分析

2023-12-23段长宏马思想

振动与冲击 2023年24期

段长宏, 李升, 马思想

(1.南京工程学院电力工程学院,南京 211167; 2.中广核新能源南通有限公司,江苏南通 210019)

近年来,电力系统的随机性明显增加,在“双碳”的目标下,新型电力系统将会呈现“双高”与“双随机”等特点,电力系统中的随机性必不可忽视[1-2]。随着科学的发展,太阳能、风能等新能源发电与电动汽车并网给电力系统带来不断的随机扰动,现代互联电网规模增加,由负荷、故障带来的随机扰动也更加普遍,传统的电力系统确定性分析法已不再适用,因此就必须从随机稳定性方面去分析随机激励对电力系统的影响[3]。

电力系统的稳定包括功角、电压和频率的稳定。功角稳定又分为小扰动与大扰动稳定。其中,当小扰动足够小的时候,可以直接采用Lyapunov第一法进行分析。对于大扰动稳定分析主要用数值方法与能量函数法。文献[4]用扩展面积法将多机系统等效为两群系统,建立起随机微分方程,基于拟哈密顿随机平均法(quasi-Hamiltonian system and stochastic averaging,QHSSA)得出多机系统电力系统能量函数,分析随机电力系统暂态稳定性,并与蒙特卡洛模拟法相比有较好的一致性。关于频率的稳定,文献[5]提出以域内概率作为随机激励对系统频率动态安全性影响的量化评估指标,并通过后向Kolmogorov方程,提出了快速求解域内概率的半解析法。随机扰动又在随机理论中称为“随机激励”,初值与参数在电力系统运行中是固定不变的,外部激励则是时变的。文献[6]基于单机系统,构建了带有随机激励功率项的随机微分方程,定量分析不同的随机激励强度,激励步长与仿真计算步长对系统功角的影响。文献[7]基于随机微积分理论,提出了一种新的随机暂态稳定性评估框架。文献[8]在确定型非线性系统的基础上引入随机激励,利用EM数值算法进行迭代,并通过均值均方稳定性证明了随机小扰动下的系统功角与转速的稳定性。文献[9]提出随机动力学模型的线性化随机激励强度最优阈值模型,得出了线性化最优阈值的条件。文献[10]提出改进EEAC法,将多机系统的随机微分方程模型通过同调分群等值为单机无穷大模型,推导出含伊藤积分的加减速面积,求出不同随机激励强度的极限切除时间。文献[11]分析考虑在随机激励影响下单机无穷大的系数矩阵特征值的3种不同情况,并证明了这3种情况下均值均差的界与随机激励强度与系统相应参数之间的联系。

前述文献都是将外部激励近似为单一的高斯型激励,但电力系统的随机性是复杂的。本文将连续且平稳的随机激励视为Gauss过程,将负荷冲击,机组投切等离散型偶发随机事件近似为复合泊松过程,利用确定型非线性系统的基本思路和方法,在此基础上增加高斯与泊松混合激励而确立为带跳随机微分方程,利用Heun数值方法进行求解,并通过随机稳定性分析方法证明高斯与泊松混合激励下的电力系统的均值、均方稳定性,给出了均值、均方的界与跳跃激励强度和随机激励强度及系统参数之间的关系,讨论相同随机激励强度下不同的跳跃激励强度与相同跳跃激励强度下不同随机激励强度对多机电力系统的功角的影响。最后,通过多机系统(4机11节点系统)的功角和转速响应曲线与系统均值、均方的界进行分析比较,验证了本文提出的理论的合理性。

1 带跳随机微分方程模型的建立

Gauss白噪声一般用来近似平稳、连续的外部随机激励,但电力系统中发生的功率突变、机组投切等冲击突发事件时用Poisson白噪声来近似更加的合适。考虑Gauss白噪声与Poisson白噪声扰动的n维向量动态系统,其带跳随机微分运动方程[12]如下

dX(t)=F[X(t),t]dt+G[X(t),t]dB(t)+
H[X(t),t]dC(t),X(t0)=X0

(1)

式(1)的积分形式为

(2)

式中:X(t0),B(t)与C(t)在正文中,是相互独立的,X(t)=[X1(t),X2(t),…,Xn(t)]T为n维矢量随机状态变量,在{B(t),t∈T}随机过程中,在时间段t1

(3)

复合Poisson过程形式导数具有Poisson白噪声特征即dC(t)=WH(t)dt。{N(t),t∈T}为服从参数λ的泊松过程,表示0～T时间内发生事件的次数,任意s,t≥0,N(t+s)-N(t)～P(λt)。C(t)为独立增量过程,{Yk,k=1,2,…,n}为一独立同分布的随机变量,且与{N(t),t≥0}为单位阶跃函数。Yk与tk分别为随机事件发生时的脉冲强度与发生事件的时刻。

2 Heun数值算法

随机微分方程的解析解很难得出,一般都是在特殊情形下才会有解析解,如今的随机微分方程的研究大部分使用数值解法。如今的随机微分方程比较成熟的数值解法有EM(Euler-Maruyama)法、Milstein法、Heun法、和Runge-Kutta法[14],EM法是其中最简单的,也是使用最多的,但其收敛阶较低,数值稳定性较低,本文采用收敛阶较高的Heun数值算法是基于梯形公式对EM算法进行改进。对式(1)进行时间离散化可近似为

(4)

Heun数值迭代格式

(5)

图1 标准Gauss白噪声与Poisson白噪声路径Fig.1 Standard Gauss white noise and Poisson white noise paths

图2 布朗运动路径与包含复合泊松的布朗运动路径Fig.2 Brownian motion path and Brownian motion path with compound Poisson

3 混合激励下多机系统的随机稳定性分析

在确定性的情况下,第i台发电机的转子运动方程为

(6)

(7)

(8)

(9)

δij=δi-δj

(10)

式中:Mi为第i台发电机的惯性;常数δij为发电机功角;ωi为发电机转速;ω0为初始角速度;t为时间;Di为阻尼系数;Pei为电磁功率;n为发电机得总台数;δij0为第i台发电机与第j台发电机之间的初值功角差;Gii为第i台发电机的自电导;Yij为第i台与第j台发电机之间的互导纳;αij第i台与第j台发电机之间阻抗角的余角;Pmi为机械功率数值等于Pei的稳态值[15]。

在式(6)的基础上增加随机激励(由Gauss过程驱动)和泊松激励(由复合Poisson过程驱动),则运动方程变为

(11)

(12)

式中:σi为第i台发电机受到高斯激励WG的强度;γi为第i台发电机受到泊松激励的强度。式(12)则为受高斯和泊松混合激励的发电机转子运动方程。

对式(9)减去式(8)可得

(13)

(14)

(15)

其中

式(15)可以写成

(16)

式(16)的解为

(17)

证明式(16)的均值稳定性

(E‖X‖)2≤E[‖X‖2]=E[XT(t)X(t)]

(18)

将式(17)代入式(18)得

根据范数的定义与随机过程的数字特征推导(具体推导过程见附录A)可得

(19)

由上述关于系统式(16)的均值、均方稳定性的证明可以判断,当满足

(20)

(21)

4 仿真分析

选用4机11节点系统来验证电力系统在高斯与泊松混合激励下的稳定性。系统参数为M1=13 s,M2=13 s,M3=12.5 s,M4=12.5 s,D1=0.2,D2=0.25,D3=0.3,D4=0.35。其他具体参数参考文献[17]。通过计算可得λA=-0.008 79,C0=2.608 25。首先取系统跳跃事件发生的概率λ=0.02,即每秒发生跳跃事件的概率为0.02,则150 s内发生3次跳跃事件, 0～150 s内发生3次跳跃的过程,如图3所示。为便于之后的研究,统一取随机激励强度σ=0.01与跳跃激励强度γ=0.01,将以上参数代入式(21)时,可得均方的界E[XT(t)X(t)]<1.461 64×10-4。在此激励强度下的功角曲线与转速曲线,如图4、图5所示。通过功角与转速的仿真曲线可以看出,虽然系统受到高斯与泊松混合激励的影响,但由于强度较小,同时转子具有惯性,功角波动不超过0.5°,转速波动不超0.005,验证系统具有均值均方稳定性。

图3 跳跃事件发生时刻与次数Fig.3 Occurrence time and number of jumping events

图4 四机系统在σi=0.01,γi=0.01的功角曲线Fig.4 Power-angle curve of the four-machine system at σi=0.01,γi=0.01

图5 四机系统在σi=0.01,γi=0.01的转速曲线Fig.5 Rotation speed curve of the four-machine system at σi=0.01,γi=0.01

下面定量分析随机激励相等的情况下不同跳跃激励强度对系统功角的影响,为了避免随机激励对系统的影响,统一取随机激励σi=0.01。分别取8组不同的跳跃激励强度来计算均值、均方的界限,具体结果如表1所示。一号机在相同随机激励强度不同跳跃激励强度的功角响应曲线图,如图6所示。由图6可知,在跳跃激励强度取表1第1组数据时,即γ=[0.01,0.02,0.03,0.04],跳跃激励强度较小,功角状态曲线波动较小,基本上可以运行在扰动前的平衡状态。图6中点线“-·”是在跳跃激励强度取第4组数据时即γ=[0.01,0.15,0.30,0.45]时的功角状态曲线图,可以看出由于跳跃激励强度的增加,在75 s,90 s与120 s发生跳跃事件时,系统功角的峰值出现了比较明显的增大,随着跳跃激励强度的增加,在发生跳跃事件时功角曲线波动也逐渐增加。而在跳跃激励强度γ=[0.01,0.7,1.4,2.1]时,仿真结果如图6实线“-”所示,在跳跃激励强度较大时,发生跳跃事件后功角响应曲线波动幅度出现了非常明显的增大且出现了严重的波动,即在跳跃激励强度增大时,系统均值、均方的上界增大,系统无法稳定运行在扰动前的平衡状态,最后可得系统失稳。

表1 相同随机激励强度下不同跳跃激励强度情况时系统式(16)的均值均方的界Tab.1 The bound of mean square of system (16) under the same random excitation intensity and different jumping excitation intensity

图6 1号机分别取表1中3组数据时的功角曲线对比Fig.6 Comparison of power angle curves of machine 1 when three groups of data in Tab.1 are respectively taken

以下分析对单个发电机跳跃激励强度发生变化时对系统状态量的影响,取γ=[0.01,0.01,0.01,0.01]与γ=[0.02,0.01,0.01,0.01],仅改变1号机的跳跃激励强度,将上述跳跃激励强度代入式(21)可得在单个发电机激励强度增加时,式(16)的均值、均方临界值也会增加。1号机本身跳跃激励强度增加时的功角曲线对比图与3号机在仅1号机跳跃激励强度发生改变时的功角曲线对比图,如图7、图8所示。由仿真可以验证,随机激励强度不改变时两种跳跃激励强度的功角曲线在发生跳跃事件之前保持一致,在任意一台发电机跳跃激励强度发生改变时,本身及其他发电机功角响应曲线波动幅度都会增大。

图8 3号机在1号机跳跃激励强度γ=0.01与γ=0.02下的功角曲线Fig.8 Power-angle curve curves of machine No.3 under the jumping excitation intensity γ=0.01 and γ=0.02 of machine No.1

下面分析相同的跳跃激励强度不同的随机激励强度下的攻角响应曲线,4台发电机跳跃激励强度统一取0.01,同理取8组不同的随机激励强度来计算系统式(16)的均值、均方稳定的界限,结果如表2所示。

表2 相同跳跃激励强度下不同随机激励强度情况时系统式(16)的均值均方的界Tab.2 Bounds of mean square of system (16) for different random excitation intensities under the same jump excitation intensities

图9为相同跳跃激励强度下不同随机激励强度的1号机功角响应曲线图,可以看出在取第1组随机激励强度时,仿真曲线如图9中虚线“--”所示,功角响应曲线波动幅度较小,基本上稳定在平衡状态。随着随机激励强度的增加,仿真如图9中点线“-·”与实线“-”所示,功角响应曲线逐步向下远离扰动前的平衡点,并逐渐失稳。

图9 1号机分别取表1中3组数据时的功角曲线对比Fig.9 Comparison of power angle curves of machine 1 when three groups of data in table 2 are respectively taken

在跳跃激励强度保持不变的同时,仅改变1号机的随机激励强度,即取σ=[0.01,0.01,0.01,0.01]与σ=[0.02,0.01,0.01,0.01]两组数据进行仿真,结果如图10～图11所示。在75 s第一次跳跃事件发生前,两种激励强度下的1号机与3号机功角曲线相差较小,但随着跳跃事件的发生,即扰动的叠加,使两种激励强度的功角曲线开始逐渐偏离,相差变大。同时,由仿真也可得出,在仅一台发电机随机激励强度发生改变时,其他发电机的功角响应曲线变化与此台发电机功角响应曲线变化相似。