一种非线性谐波振荡的近似解析求解与稳定域分析方法与流程

1.本发明属于电力系统电能质量及稳定性分析领域，具体涉及一种非线性谐波振荡的近似解析求解与稳定域分析方法。


背景技术：

2.在我国分布式光伏发电、分散式风电等新能源发电装机容量快速增长背景下，通过构网型灵活性资源如燃气轮机发电机组、电动汽车储能电站、变频负荷等的灵活性调节构网能力，可有效提高新能源电力消纳比例，提升分布式新能源利用水平和配电系统运行效率，并保障配电网安全稳定供能。
3.然而，当非线性电力电子装置随着分布式光伏发电、分散式风电广泛接入配电系统时，在多个电力电子设备的暂动态累积影响下，配电系统电压和电流谐波表现为交互影响，其交互动态的非线性特性使得谐波等电能质量治理和稳定问题分析较为困难。


技术实现要素：

4.为了避免配电系统可靠运行的潜在风险和评估各种谐波的非线性交互影响，本发明提出了一种非线性谐波振荡的近似解析求解与稳定域分析方法。
5.为达到上述目的，本发明采用的技术方案如下：
6.一种非线性谐波振荡的近似解析求解与稳定域分析方法，针对含有构网型灵活性调节资源的电力电子化配用电网，建立具有普适性的非线性动力学模型，该模型的自然频率表现为周期变化，进一步通过分析自然频率的变化特征，依次分析ω＝ω0、ω＝2ω0、ω＝3ω0的情况下δ的近似解析解，ω为系统时变自然频率，ω0为系统特定自然频率，δ为非线性动力学系统轨迹相位角，并通过多阶正弦型和余弦型马蒂厄及其特征函数获得近似稳定边界，最后分析线性阻尼对近似解析解与稳定边界的影响，具体步骤包括：
7.步骤1：建立构网型电源幅相运动自然频率周期变化的动力学模型；
8.步骤2：求解非线性模型ω＝ω0和ωv＝2ω0情况下近似解析解和稳定边界，ωv为考虑构网型电源幅相幅值变化的自然频率；
9.步骤3：求解非线性模型ω＝2ω0和ωv＝2ω0情况下近似解析解和稳定边界；
10.步骤4：求解非线性模型ω＝3ω0和ωv＝2ω0情况下近似解析解和稳定边界；
11.步骤5：分析线性阻尼对非线性振荡稳定的影响。
12.进一步地，所述的步骤1中，建立构网型电源幅相运动自然频率周期变化的动力学模型，即采用一元即相位θ的二阶微分方程表达，幅值的变化通过二阶微分方程周期性变化的自然频率来表示，因此表现相位和幅值交互影响的非线性动力学模型为：
[0013][0014]
式中，θ为幅相运动中的相位变化值，γ为该幅相运动的阻尼，即ωn为幅值运动轨迹的自然频率，表现为准周期性变化，ωv为周
期性变化的频率，ε为小参量，t为时间，ω
′
为非周期性变化的频率，为二阶幅相运动方程相位的二阶导数，为二阶幅相运动方程相位的一阶导数；
[0015]
将坐标变换即θ＝e-γt
δ代入式(1)中，δ为系统变换后二阶幅相运动方程的相位，t为时间，γ为系统的阻尼，并采用马蒂厄方程形式，将式(1)幅值相位表示为如下非线性方程：
[0016][0017]
为系统变换后二阶幅相运动方程相位的二阶导数。
[0018]
进一步地，所述的步骤2中，通过两个线性无关的特解δ0＝cosω0t和δ0＝sinω0t，δ0为二阶幅相运动方程的特解，分别求解出该非线性方程的近似解析解如下：
[0019][0020]
ε满足非线性系统自然频率约束：
[0021][0022]
以及将δ0＝sinω0t代入非线性方程中计算得到：
[0023][0024]
相应地，ε满足非线性系统自然频率约束如下：
[0025][0026]
其中，为系统瞬时变化的自然频率。
[0027]
进一步地，所述的步骤3具体包括：近似解析解表现为余弦和正弦形式为：
[0028][0029][0030]
其系统稳定边界满足系统摄动量ε与系统自然频率之间的关系为：
[0031][0032][0033]
进一步地，所述的步骤4具体包括：近似解析解表现为余弦和正弦形式为：
[0034][0035][0036]
其系统稳定边界满足系统摄动量ε与系统自然频率之间的关系为：
[0037][0038][0039]
进一步地，所述的步骤5中具体包括：得到线性阻尼条件下的非线性稳定边界表达式如下：
[0040][0041][0042][0043][0044][0045]
上述式中，γ为阻尼系数；
[0046]
通过绘制式(29)-(33)的ε～ω
′
坐标曲线，绘制出非线性振荡系统的估计稳定区域。
[0047]
有益效果：
[0048]
相比传统基于事后曲线数据分析谐波振荡的方法，本发明从揭示非线性振荡模型机理切入，建立非线性振荡方程的近似解析求解方法，包括能够直接获得近似解析解、非线性二倍频振荡机理，其近似解析方法求解振荡轨迹速度更快，其对比如表1所示。
[0049]
表1显示了针对非线性二阶系统50秒、100秒、200秒动态特性的计算时长，由表可见：传统基于事后曲线数据进行傅里叶变换分析方法，所需计算时长分别为20.85秒、43.22秒、84.24秒；而采用本发明解析方法计算时长分别为7.26秒、12.15秒和23.18秒，加速比分别为65.18％、71.89％和72.48％。
[0050]
表1
[0051]
仿真时长事后数据分析方法本发明解析方法加速比50秒20.85秒7.26秒65.18％100秒43.22秒12.15秒71.89％200秒84.24秒23.18秒72.48％
[0052]
进一步地，本发明能够直接获得非线性振荡稳定边界。
附图说明
[0053]
图1为小扰动下面向谐波交互影响的谐波动力学相位轨迹示意图；
[0054]
图2为大扰动下面向谐波交互影响的谐波动力学相位轨迹示意图；
[0055]
图3为谐波动力学非线性自然频率曲线；
[0056]
图4为谐波动力学非线性交互的稳定域示意图。
具体实施方式
[0057]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0058]
下面结合附图与具体实施方式进一步说明本发明。
[0059]
本发明的一种非线性谐波振荡的近似解析求解与稳定域分析方法，包括如下步骤：
[0060]
步骤1：建立构网型幅相运动自然频率周期变化的动力学模型。
[0061]
为方便描述，幅相运动通常采用一元即相位二阶微分方程表达，而幅值的变化可通过二阶微分方程周期性变化的自然频率来表示，因此表现相位和幅值交互影响的非线性动力学模型可表达为：
[0062][0063]
式中，θ为幅相运动中的相位变化值，γ为该幅相运动的阻尼，即ωn为幅相运动轨迹的自然频率，表现为准周期性变化，ωv为周期性变化的频率，ε为小参量，t为时间，ω
′
为非周期性变化的频率，为二阶幅相运动方程相位的二阶导数，为二阶幅相运动方程相位的一阶导数。
[0064]
采用坐标变换即θ＝e-γt
δ，δ为系统变换后二阶幅相运动方程的相位，t为时间，γ为系统的阻尼，并将该变换代入式(1)中可得：
[0065][0066]
其中，为系统变换后二阶幅相运动方程相位的二阶导数。
[0067]
进一步地，定义ω
′
2-γ2＝ω2+εσ1+ε2σ2，σ1与σ2为扰动量，因此采用马蒂厄(mathieu)方程形式，可将式(2)幅值相位表示为如下非线性方程：
[0068][0069]
上式中，通常情况下，假定等值模型ωv＝2ω0值固定不变，随着系统结构参数的变化，ω≈ω0、ω≈2ω0、ω≈3ω0或者取ω0更高倍数的值。高频ω即(ω＞＞ωv)分量对系统性稳定性影响不大(或者系统运行状态通常处于稳定区域内，或者系统稳定域相对低频量较大)。同时为降低超低频ω即(ω＜＜ωv)对系统稳定分析计算的求解难度，研究中忽略超低频ω对系统稳定的影响。
[0070]
采用林茨泰德-庞加莱摄动方法，将式(3)中的非线性解δ展开成ε的二次级数形式：
[0071]
δ＝δ0+εδ1+ε2δ2ꢀꢀꢀ
(4)
[0072]
其中，δ0为基频项，δ1为谐波二阶项，δ2为谐波三阶项。
[0073]
将式(4)代入非线性方程(3)中，令等式两边ε的各次级数的系数相等，即可求出各阶近似的线性微分方程组，如下：
[0074]
[0075][0076][0077]
如下步骤将分别讨论ω＝ω0、ω＝2ω0、ω＝3ω0的情况下，δ的近似解析解及其近似稳定边界。
[0078]
步骤2：求解ω＝ω0和ωv＝2ω0情况下非线性方程近似解析解和稳定边界。
[0079]
当ω＝ω0时，方程(5a)具有两个线性无关的特解δ0＝cosω0t和δ0＝sinω0t。首先将δ0＝cosω0t代入方程(5b)中得到：
[0080][0081]
为避免久期项，令求解(6)可得：
[0082][0083]
进一步地，将δ0＝cosω0t和式(7)代入方程(5c)中得到：
[0084][0085]
为满足周期解条件，可令进而求解式(8)得到：
[0086][0087]
联合式(5a)、式(7)和式(9)，可求得该非线性方程的近似解析解如下：
[0088][0089]
由上式可见，当ε取值较大时，δ将出现宽频振荡。另外，ε满足非线性系统自然频率约束：
[0090][0091]
其中，为系统瞬时变化的自然频率。
[0092]
另一方面，将δ0＝sinω0t代入非线性方程中计算可得：
[0093][0094]
相应地，ε满足非线性系统自然频率约束如下：
[0095][0096]
步骤3：求解ω＝2ω0和ωv＝2ω0情况下非线性方程近似解析解和稳定边界。
[0097]
当ω＝2ω0时，通过式(5a)能够求解出线性无关的特解δ0＝cos2ω0t和δ0＝sin2ω0t，在此基础上，先将cos2ω0t代入(5b)，得到：
[0098]
[0099]
为避免久期项，可令σ1＝0，即可求解出：
[0100][0101]
进一步地，将δ0＝cos2ω0t和式(15)代入(5c)中，可得到δ2的表达式如下
[0102][0103]
为简化计算，假设δ2的解满足周期解条件，则令进一步求解式(16)得：
[0104][0105]
结合δ0＝cos2ω0t、式(15)和式(17)，可得到式(3)的近似解析解，即更新式(10)如下
[0106][0107]
则式(18)近似解析解表现出的振荡特性将取决于摄动量ε，因此需寻找摄动量ε与系统自然频率之间的关系，进而分析该非线性系统的结构稳定性：
[0108][0109]
若改用特解δ0＝sin2ω0t，则可计算出式(3)的近似解析解如下：
[0110][0111]
同时摄动量ε与系统自然频率之间的关系如下：
[0112][0113]
步骤4：求解ω＝3ω0和ωv＝2ω0情况下非线性方程近似解析解和稳定边界。
[0114]
当ω＝3ω0时，通过式(3a)能够求解出线性无关的特解δ0＝cos3ω0t和δ0＝sin3ω0t，在此基础上，先将cos3ω0t代入(3b)，为避免久期项，令σ1＝0，则可得到：
[0115][0116]
进一步地，将δ0＝cos3ω0t和式(22)中δ1带入(3c)，重新整合得到：
[0117][0118]
为避免久期项，令则可进一步得到：
[0119][0120]
则结合δ0＝cos3ω0t、式(22)和(24)，能够求解非线性方程的近似解析解如下：
[0121][0122]
对应地，稳定边界可采用非线性动力学系统自然频率ω2与ε的关系获得：
[0123][0124]
另一方面，当ω＝3ω0时，将特解δ0＝sin3ω0t带入式(3a)-(3c)中，进而获得近似解析解如下：
[0125][0126]
其稳定边界可表达为：
[0127][0128]
步骤5：线性阻尼对非线性振荡稳定的影响。
[0129]
结合ω
′
2-γ2＝ω2+εσ1+ε2σ2，并将式(11)和(13)，式(19)和(21)，式(26)和(28)带入即采用ω
x
替代ω，可得到线性阻尼条件下的非线性稳定边界表达式如下：
[0130][0131][0132][0133][0134][0135]
上述式中，γ为阻尼系数。通过绘制式(29)-(33)的ε～ω
′
坐标曲线，即可绘制出非线性振荡系统的估计稳定区域。
[0136]
下面结合附图进一步说明本发明。
[0137]
图1为小扰动下面向谐波交互影响的谐波动力学相位轨迹示意图；该图显示了非线性二阶动力学系统在小扰动下相角随时间的变化曲线，实线采用的摄动参数ε＝13.05，而虚线采用的ε＝43.05。当采用的摄动参数越大，其自然频率中的二倍频项系数越大，即二倍频项对动力学系统的动态特性影响越大，突破了阻尼带来的抑制作用，进而动力学曲线显示了振荡周期失稳现象。
[0138]
图2为大扰动下面向谐波交互影响的谐波动力学相位轨迹示意图；该图显示了非线性二阶动力学系统在大扰动下相角随时间的变化曲线，实线采用的摄动参数ε＝13.05，而虚线采用的ε＝43.05。当采用的摄动参数越大，其自然频率中的二倍频项系数越大，即二倍频项对动力学系统的动态特性影响越大，突破了阻尼带来的抑制作用，进而动力学曲线显示了振荡周期失稳现象。相比图1中曲线，大扰动下振荡幅度更大，即大扰动带来的激励作用更大。
[0139]
图3为谐波动力学非线性自然频率曲线；该图中，低摄动参数即ε＝13.05，自然频率振荡幅度较小，而高摄动参数ε＝43.05下自然频率振荡幅度较大，由此可见在本发明中，高摄动参数对系统非线性振荡的影响较大。
[0140]
图4为谐波动力学非线性交互的稳定域示意图。图中，黑色区域表示稳定区域，纵
坐标代表摄动参数ε，横坐标代表自然频率的平方；由图可见：摄动参数ε越大，系统越接近不稳定区域，且自然频率分别在ω＝ω0、ω＝2ω0、ω＝3ω0的情况下，系统稳定程度越差，即较小的摄动参数ε也会引起系统不稳定(系统不稳定区域在ε
→
0时接近横坐标)。
[0141]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。