CN102790389A - 一种暂态稳定裕度概率分布的计算方法 - Google Patents

Abstract

本发明公开了一种暂态稳定裕度概率分布的计算方法，将分布在电力系统中的各个节点功率的概率数字特征，变换成失稳模式对应的两群驱动功率和负荷功率的概率数字特征；根据扩展等面积准则EEAC，将系统稳定裕度按Taylor级数表示成以两群驱动功率和负荷功率为变量、以稳定裕度灵敏度为系数的线性化表达式，以求取稳定裕度的概率数字特征；再利用Gram-Charlier级数计算出系统暂态稳定裕度的概率分布。可直接用于指导电力系统规划、设计、运行和控制等各个环节中安全稳定量化分析与优化决策的具体实践。

Description

一种暂态稳定裕度概率分布的计算方法

技术领域

本发明属于电力系统及其自动化技术领域，具体涉及一种应用在节点功率注入不确定条件下的暂态稳定裕度概率分布的计算方法。

背景技术

电力工业是国民经济发展中最重要的基础能源产业，电力系统的安全稳定高效运行具有重大意义。大规模新能源发电的并网运行和需求侧响应机制是现代电网的重要特征，但风能、太阳能发电的随机性、间歇性和波动性特点，以及用电侧灵活多变的需求方式，却给电力系统的安全稳定特性带来了深刻的变化。

稳定分析是系统安全性评估的不可或缺的一部分。传统的稳定分析方法局限于用确定性的观点对待系统工况和扰动场景，凭经验根据时域仿真得到的系统受扰轨迹进行稳定性的定性判断，因而无法量化主要随机因素对稳定性的影响。针对传统稳定性分析准则所面临的挑战，国内外的发展趋势是采用风险的观点审视电力系统的规划、设计、运行和控制等各个环节，而概率暂态稳定分析是实现安全性风险评估技术的基础。

电力系统概率暂态稳定分析的任务是根据影响稳定性的主要随机因素的统计特性来确定系统的暂态稳定性概率指标。已有研究成果中，蒙特卡罗仿真法基于对随机变量的统计试验和随机模拟，稳定域法采用系统工况点处于稳定域外的比例来反映系统失稳概率。这些方法只给出系统稳定或失稳的概率，却无法评估系统稳定或失稳的程度，更无法提供决策所需要的灵敏度信息，而且，由于计算费时或方法复杂等原因，难以处理不确定变量个数较多的情况。

间歇性电源和用电负荷分布在电力系统广域空间中，构成了大量的节点功率注入不确定性因素。鉴于在针对确定性系统工况和扰动场景的安全稳定评估技术领域，扩展等面积准则(EEAC)是惟一得到严格证明并投入工程应用的稳定性量化分析工具，可给出稳定裕度、稳定模式、控制措施灵敏度等信息，因此，有必要将EEAC向处理不确定因素的方向扩展，提出新的节点功率注入不确定条件下暂态稳定裕度概率分布的计算方法。

发明内容

针对现有技术的不足，本发明提供在一种节点功率注入不确定条件下的暂态稳定裕度概率分布的计算方法，为含有分布在广域空间中多个随机变量的电力系统稳定裕度概率分布计算提供了高效率的解决方法。不确定是指节点功率注入量随机，但遵循概率分布。

本发明提供的一种暂态稳定裕度概率分布的计算方法，其改进之处在于，所述方法包括如下步骤：

（1）确定初始状态下受扰系统的失稳模式和稳定裕度初始值η₀，并由用户确定采用的计算模式，包括基于纯解析法的暂态稳定裕度概率分布计算和综合时域仿真与解析法的暂态稳定裕度概率分布计算；

所述基于纯解析法的暂态稳定裕度概率分布计算包括如下步骤：

（2）计算所述受扰系统稳定裕度对临界群S驱动功率的灵敏度系数I：令临界群S的驱动功率P_ms和余下群A的驱动功率P_ma的初始值分别等于各自群内机组驱动功率均值的总和，通过静态EEAC（SEEAC）算法求取稳定裕度对临界群驱动功率的灵敏度系数I；所述灵敏度系数I包括一阶灵敏度系数

和二阶灵敏度系数

其中临界群和余下群中发电机机组通常有多个，临界群S的驱动功率P_ms和余下群A的驱动功率P_ma的初始值分别等于各自群内不同机组的驱动功率均值求和。

（3）计算临界群S驱动功率的概率数字特征，即各阶半不变量γ_υ(P_ms)：根据临界群内各发电机组驱动功率的概率分布，先求取每台发电机驱动功率P_mk的各阶半不变量γ_υ(P_mk)，再按临界群S中发电机组个数对所述各阶半不变量γ_υ(P_mk)求和，求取P_ms的各阶半不变量γ_υ(P_ms)，为

${\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ms}}\right)=\underset{k=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{mk}}\right)$

其中，k为机组编号，υ表示半不变量的阶数，Ns为临界群中发电机组个数。

（4）求取受扰系统稳定裕度的各阶中心矩β_υ(η)：利用稳定裕度η对临界群S驱动功率P_ms的灵敏度系数I和临界群S驱动功率的各阶半不变量γ_υ(P_ms)，计算稳定裕度的各阶半不变量γ_υ(η)，进而通过半不变量和中心矩的标准数学转换公式，求取稳定裕度的各阶中心矩β_υ(η)。

（5）计算稳定裕度的累积分布函数，完成理想两群失稳模式的暂态稳定裕度概率分布计算。将稳定裕度的初始值η₀及其各阶中心矩β_υ(η)代入Gram-Charlier级数计算公式，获得受扰系统稳定裕度η的累积分布函数，完成了基于纯解析法、针对理想两群失稳模式的暂态稳定裕度概率分布的快速计算。其中，理想两群是指互补群群内各发电机转子角轨迹间的间隙在整个动态过程中保持不变。

所述综合时域仿真与解析法的暂态稳定裕度概率分布计算包括如下步骤：

（6）确定节点功率变化范围内的失稳模式变化情况：根据系统中具有功率不确定性的发电机和负荷功率变化范围的上下限，计算出对应步骤（1）初始分群的临界群驱动功率P_ms、余下群驱动功率P_ma、临界群负荷功率P_ls、余下群负荷功率P_la变化范围的上下限，分别按照增加临界群驱动功率和余下群负荷功率，同时减少余下群驱动功率和临界群负荷功率，以及减少临界群驱动功率和余下群负荷功率，同时增加余下群驱动功率和临界群负荷功率的两个方向，依据人工设定的功率变化步长及其在各机组和负荷之间的分配规则，调整并形成一组系统稳态运行工况，再针对每一个系统稳态运行工况，进行针对设定故障下的详细时域仿真，识别失稳模式变化，并记录每一个新的失稳模式i所对应的稳定裕度η′_(i)，以及临界群的驱动功率P_ms0(i)、临界群的负荷功率P_ls0(i)、余下群驱动功率P_ma0(i)和余下群的负荷功率P_la0(i)。节点功率变化范围是指所计算系统中节点功率实际变化范围，由各节点功率概率分布的上下限确定，初始分群是指系统中各节点功率等于其概率分布均值时，针对设定故障，利用详细时域仿真和扩展等面积准则(EEAC)得到的系统分群情况。

（7）针对每个失稳模式，计算稳定裕度的灵敏度系数II：对每个失稳模式，采用摄动法，通过时域仿真和EEAC准则求取稳定裕度对所述临界群驱动功率P_ms的灵敏度系数

临界群负荷功率P_ls的灵敏度系数

余下群驱动功率P_ma的灵敏度系数和余下群负荷功率P_la的灵敏度系数

（8）针对每个失稳模式，求取该失稳模式下稳定裕度线性变化对应的各节点功率等于其概率分布均值时的稳定裕度初始值。若失稳模式的个数为1，则稳定裕度初始值η′₀₍₁₎＝η₀；若失稳模式有变化，则根据稳定裕度的线性化展开式，结合步骤（7）得到的灵敏度系数II和步骤（6）得到的失稳模式变化后稳定裕度η′_(i)，利用数学反推获得每个失稳模式下的稳定裕度初始值η′_0(i)。

（9）针对每个失稳模式，计算临界群驱动功率P_ms的各阶半不变量、临界群负荷功率P_ls的各阶半不变量、余下群驱动功率P_ma的各阶半不变量和余下群负荷功率P_la的各阶半不变量：先根据系统内所有机组驱动功率和负荷功率的概率分布，求取单个节点的机组驱动功率和负荷功率的各阶半不变量γ_υ，再根据每个失稳模式的分群情况，分别按照临界群及余下群内机组和负荷节点的个数分类求和，获得P_ms、P_ls、P_ma和P_la的各阶半不变量。

（10）针对每个失稳模式，求取稳定裕度η的各阶中心矩：利用步骤（7）得到的稳定裕度对所述临界群驱动功率P_ms的灵敏度系数II、临界群负荷功率P_ls的灵敏度系数II、余下群驱动功率P_ma的灵敏度系数II和余下群负荷功率P_la的灵敏度系数II，以及步骤（9）得到的所述临界群驱动功率P_ms的各阶半不变量、临界群负荷功率P_ls的各阶半不变量、余下群驱动功率P_ma的各阶半不变量和余下群负荷功率P_la的各阶半不变量，计算稳定裕度η的各阶半不变量，进而通过半不变量和中心矩的标准数学转换公式，求取稳定裕度η的各阶中心矩β_υ(η)。

（11）求取稳定裕度的累积分布函数：针对每个失稳模式，将稳定裕度的各阶中心矩β_υ(η)和稳定裕度初始值η′_0(i)代入Gram-Charlier级数计算公式，获得稳定裕度η的累积分布函数。若在节点功率变化范围内失稳模式发生变化，则按失稳模式进行分段，通过分段函数综合每个失稳模式对应的稳定裕度概率分布，完成了综合采用时域仿真与解析方法的暂态稳定裕度概率分布计算。

其中，所述步骤（4）中，在计算稳定裕度η的方差D(η)（即2阶中心距）时，按式(1)计入临界群驱动功率P_ms方差E(Δ²P_ms)的影响，而在计算稳定裕度的高阶（大于等于3阶）中心矩时则忽略P_ms方差的影响；其中

$D\left(\mathrm{\η}\right)=S_{P_{\mathrm{ms}}}^{2}E\left({\mathrm{\Δ}}^2{P}_{\mathrm{ms}}\right)+\frac{1}{4}{\left[S_{P_{\mathrm{ms}}}^{\left(2\right)}E\left({\mathrm{\Δ}}^2{P}_{\mathrm{ms}}\right)\right]}^2---\left(1\right)$

其中，所述步骤（4）中，按式(2)求取稳定裕度η对临界群驱动功率P_ms的一阶灵敏度系数

按式(3)求取稳定裕度的各阶半不变量γ_υ(η)；

$S_{P_{\mathrm{ms}}}=\frac{\∂\mathrm{\η}}{{\∂P}_{\mathrm{ms}}}+\frac{\∂\mathrm{\η}}{{\∂P}_{\mathrm{ma}}}\frac{{\mathrm{dP}}_{\mathrm{ma}}}{{\mathrm{dP}}_{\mathrm{ms}}}---\left(2\right)$

${\mathrm{\γ}}_{\mathrm{\υ}}\left(\mathrm{\η}\right)=S_{P_{\mathrm{ms}}}^{\mathrm{\υ}}{\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ms}}\right)---\left(3\right)$

其中，所述步骤（10）中，按式(4)根据时域仿真和EEAC获得的灵敏度系数

和求取稳定裕度η的各阶半不变量γ_υ(η)；

${\mathrm{\γ}}_{\mathrm{\υ}}\left(\mathrm{\η}\right)=\left[\begin{array}{cccc}{S}_{P_{\mathrm{ms}}}^{\mathrm{\υ}}& S_{P_{\mathrm{ls}}}^{\mathrm{\υ}}& S_{P_{\mathrm{ma}}}^{\mathrm{\υ}}& S_{P_{\mathrm{la}}}^{\mathrm{\υ}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ms}}\right)\\ {\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ls}}\right)\\ {\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ma}}\right)\\ {\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{la}}\right)\end{array}\right]---\left(4\right)$

其中，所述步骤（8)中，对每个失稳模式i，根据稳定裕度的线性化展开关系，按式(5)进行数学反推，求取各节点功率等于其概率分布均值时的稳定裕度初始值η′_0(i)。

${\mathrm{\η}}_{\left(i\right)}^{\′}={\mathrm{\η}}_{0\left(i\right)}^{\′}+\left[\begin{array}{cccc}{S}_{P_{\mathrm{ms}\left(i\right)}}& S_{P_{\mathrm{ls}\left(i\right)}}& S_{P_{\mathrm{ma}\left(i\right)}}& S_{P_{\mathrm{la}\left(i\right)}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ΔP}}_{\mathrm{ms}\left(i\right)}\\ {\mathrm{\ΔP}}_{\mathrm{ls}\left(i\right)}\\ {\mathrm{\ΔP}}_{\mathrm{ma}\left(i\right)}\\ {\mathrm{\ΔP}}_{\mathrm{la}\left(i\right)}\end{array}\right]---\left(5\right)$

式中 ${\mathrm{\ΔP}}_{\mathrm{ms}\left(i\right)}=P_{\mathrm{ms}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{ms}},$ ${\mathrm{\ΔP}}_{\mathrm{ls}\left(i\right)}=P_{\mathrm{ls}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{ls}},$ ${\mathrm{\ΔP}}_{\mathrm{ma}\left(i\right)}=P_{\mathrm{ma}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{ma}},$ ${\mathrm{\ΔP}}_{\mathrm{la}\left(i\right)}=P_{\mathrm{la}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{la}},$

和

分别为各节点功率等于其概率分布均值时的临界群驱动功率、临界群负荷功率、余下群驱动功率和余下群负荷功率。

与现有技术比，本发明的有益效果为：

1）本发明将分布在电力系统中的各个节点功率的概率数字特征，变换成失稳模式对应的两群驱动功率和负荷功率的概率数字特征，再与系统稳定裕度进行数学关联，为含有分布在广域空间中多个随机变量的电力系统稳定裕度概率分布计算提供了高效率的解决方法。

2）本发明将基于扩展等面积准则（EEAC）的稳定裕度线性化展开式与Gram-Charlier级数有机结合，先求解系统稳定裕度的概率数字特征，再计算其概率分布，避免了大量的时域仿真计算，提高了概率暂态稳定评估的效率。

3）本发明提出的纯解析方法能够计及系统内发电机机组驱动功率的不确定性，在理想两群失稳模式下计算精度高，具有速度快的特点；综合采用时域仿真与解析方法的混合算法能计及机组驱动功率和负荷功率的不确定性，通过时域仿真求取各项灵敏度系数，适用于理想两群、非理想两群和多个失稳模式的情况，具有精度高的优点。

4）本发明为系统运行工况不确定性条件下的概率稳定裕度计算提供了解决方案，可直接用于指导电力系统规划、设计、运行和控制等各个环节中安全稳定量化分析与优化决策的具体实践。

附图说明

图1为本发明提供的一种暂态稳定裕度概率分布的计算方法流程图。

图2为本发明提供的3机9节点系统接线图。

图3为本发明提供的单一失稳模式下暂态稳定裕度概率分布计算结果示意图。

图4为本发明提供的系统稳定裕度随临界群驱动功率变化的曲线图。

图5为本发明提供的用分段函数描述的系统稳定裕度概率分布图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步的详细说明。

本发明的主要过程是：将分布在电力系统中的各个节点功率的概率数字特征，变换成失稳模式对应的两群驱动功率和负荷功率的概率数字特征；根据扩展等面积准则EEAC，将系统稳定裕度按Taylor级数表示成以两群驱动功率和负荷功率为变量、以稳定裕度灵敏度为系数的线性化表达式，以求取稳定裕度的概率数字特征；再利用Gram-Charlier级数计算出系统暂态稳定裕度的概率分布。本发明提供了两种计算模式，其中纯解析方法能够计及系统内发电机机组驱动功率的不确定性，在理想两群失稳模式下计算精度高，具有速度快的特点；综合采用时域仿真与解析方法的混合算法能计及机组驱动功率和负荷功率的不确定性，通过时域仿真求取各项灵敏度系数，并用分段函数的方法解决由失稳模式变化导致稳定裕度或灵敏度系数突变的描述问题，适用于理想两群、非理想两群和多个失稳模式的情况，具有精度高的优点。

本实施例的流程图如图1所示，本实施将本发明方法应用于IEEE标准3机9节点系统，其系统接线图如图2所示，在已知发电机G2驱动功率概率分布的前提下，计算系统在线路8-9三相短路故障下的稳定裕度概率分布。为便于说明问题，假设系统中仅发电机G2的驱动功率具有随机性。

（1）令发电机G2的驱动功率等于其驱动功率概率分布均值，运用EEAC对该稳态运行工况进行详细时域仿真，得到系统的失稳模式(发电机G2为临界群，发电机G1,G3为余下群)，以及系统稳定裕度初始值η₀，本实施例分别用基于纯解析法的暂态稳定裕度概率分布计算和综合时域仿真与解析法的暂态稳定裕度概率分布计算两种计算模式进行计算说明。

基于纯解析法的暂态稳定裕度概率分布计算包括如下步骤：

（2）由SEEAC中灵敏度系数求解公式，求出系统稳定裕度对临界群驱动功率的一阶灵敏度系数

和二阶灵敏度系数

（3）通过发电机G2驱动功率的概率分布，求取G2驱动功率的各阶半不变量γ_υ(P_m2)。由于该失稳模式下临界群仅含有G2一台发电机，因此临界群驱动功率的各阶半不变量即等于G2驱动功率的各阶半不变量γ_υ(P_ms)＝γ_υ(P_m2)；

（4）按照式(1)计算稳定裕度的方差，即2阶中心距，式中E(Δ²P_ms)为临界群驱动功率的方差，在本例中即为发电机G2驱动功率的方差；概率分布的数值计算表明，在Gram-Charlier级数中，越是高阶的半不变量，对最终累积分布函数的影响就越小，为了减少计算量，忽略6阶及以上的各阶半不变量，按照式(3)计算系统稳定裕度的3，4，5阶半不变量，然后通过半不变量和中心距的数学转换公式求出稳定裕度的3，4，5阶中心距；

（5）将稳定裕度的初始值η₀及其各阶中心矩β_υ(η)代入Gram-Charlier级数计算公式，获得稳定裕度η的累积分布函数，完成基于纯解析法的暂态稳定裕度概率分布的快速计算。结果如附图3所示，横坐标为系统稳定裕度，纵坐标为概率。由图可知系统稳定裕度不大于横坐标所对应稳定裕度值的概率，例如A点表示系统在设定故障下稳定裕度不大于40%的概率为0.7。

综合时域仿真与解析法的暂态稳定裕度概率分布计算包括如下步骤：

（6）由于本实施例中仅有发电机G2的驱动功率具有不确定性，其余发电机驱动功率及负荷功率是确定的，因此只需要根据发电机G2驱动功率的概率分布即可计算出步骤（1）初始分群的临界群驱动功率变化范围的上下限。依据人工设定的功率步长，分别按照逐渐增加G2的驱动功率和逐渐减少G2的驱动功率这两个方向调整运行方式，形成新的运行工况，并对每一个新的运行工况进行详细时域仿真。

得到附图4所示的系统稳定裕度随临界群驱动功率的变化曲线，图中η_m和η′_m均在时域仿真时获得。由图可知：当发电机G2驱动功率在0到P′_ms之间时，系统失稳模式的临界群仅包含发电机G2；当发电机G2驱动功率大于P′_ms时，系统失稳模式发生变化，临界群中包含发电机G2和G3。此时，系统分群情况发生变化：由临界群仅包含发电机G2、余下群包含发电机G1和发电机G3，变成了临界群中包含发电机G2和G3、余下群仅包含发电机G1。图中曲线不连续的原因是，当失稳模式发生变化后，系统的临界群驱动功率变为G2和G3驱动功率之和；

（7）分别针对上述2个失稳模式，通过时域仿真摄动求取稳定裕度对临界群驱动功率P_ms的一阶灵敏度系数

和

由于本实施例中负荷功率及余下群驱动功率是确定的，因此不必要计算和

（8）针对失稳模式1(临界群包含发电机G2，余下群包含发电机G1和G3)，根据步骤1），稳定裕度初始值η′₀₍₁₎＝η₀；针对失稳模式2(临界群包含发电机G2和G3,余下群包含发电机G1)，根据式(5)求解稳定裕度初始值η′₀₍₂₎；

（9）针对每个失稳模式，计算临界群驱动功率的各阶半不变量。第一个失稳模式下，γ_υ(P_ms)＝γ_υ(P_m2)；第二个失稳模式下，由于G3驱动功率恒定，所以γ_υ(P_ms)＝γ_υ(P_m2)依旧成立；

（10）针对每个失稳模式，根据式(4)求取稳定裕度的各阶半不变量，继而根据半不变量和中心距的数学转换关系，求取稳定裕度的各阶中心距；

（11）针对每个失稳模式，将稳定裕度的各阶中心矩和稳定裕度初始值代入Gram-Charlier级数计算公式，分别得到两个失稳模式下的稳定裕度概率分布曲线。由于系统失稳模式变化，临界群包含机组不同，需要结合两个失稳模式下的稳定裕度概率分布，通过分段函数综合描述系统的稳定裕度概率分布，如图5所示。至此，完成了综合采用时域仿真与解析方法的暂态稳定裕度概率分布计算。

本实施例的字母含义说明：

η：表示受扰系统的稳定裕度，其中η₀为稳定裕度初始值；

S：表示临界群；A：表示余下群；

P_ms：表示临界群的驱动功率，其中P_mk表示临界群内每台发电机的驱动功率；

分别表示稳定裕度对临界群驱动功率的一阶灵敏度系数和二阶灵敏度系数；

P_ls：表示临界群的负荷功率；

P_ma：表示余下群的驱动功率；P_la：表示余下群的负荷功率；

表示第i个失稳模式下稳定裕度对临界群驱动功率的一阶灵敏度系数；

表示第i个失稳模式下稳定裕度对临界群负荷功率的一阶灵敏度系数；

表示第i个失稳模式下稳定裕度对余下群驱动功率的一阶灵敏度系数；

表示第i个失稳模式下稳定裕度对余下群负荷功率的一阶灵敏度系数；

γ_υ(P_ms)：表示临界群驱动功率的各阶半不变量，其中γ_υ(P_mk)表示每台发电机驱动功率的各阶半不变量；

β_υ(η)：表示受扰系统的稳定裕度的各阶中心矩；

γ_υ(η)：表示受扰系统的稳定裕度的各阶半不变量；

η′_(i)：表示所述步骤(6)中，随着节点功率的变化，系统由第i-1个失稳模式变化成第i个失稳模式后，新失稳模式对应运行工况下的系统稳定裕度；

η′₀：表示仅有一个失稳模式时，系统中各发电机驱动功率、负荷功率等于其概率分布均值时的系统稳定裕度，亦称为该失稳模式下的稳定裕度初始值；η′_0(i)表示有多个失稳模式时，按照第i个失稳模式对系统进行分群，系统中各发电机驱动功率、负荷功率等于其概率分布均值时的系统稳定裕度，在本发明中亦称为第i个失稳模式下的稳定裕度初始值；

P_ms0(i)：表示随着节点功率的变化，系统由第i-1个失稳模式变化为第i个失稳模式后，新失稳模式对应运行工况下的系统临界群驱动功率；

P_ma0(i)：表示随着节点功率的变化，系统由第i-1个失稳模式变化为第i个失稳模式后，新失稳模式对应运行工况下的系统临界群负荷功率；

P_ls0(i)：表示随着节点功率的变化，系统由第i-1个失稳模式变化为第i个失稳模式后，新失稳模式对应运行工况下的系统余下群驱动功率；

P_la0(i)：表示随着节点功率的变化，系统由第i-1个失稳模式变化为第i个失稳模式后，新失稳模式对应运行工况下的系统余下群负荷功率；

γ_υ：表示第ν阶半不变量；

D(η)：表示受扰系统稳定裕度的方差（2阶中心矩）；

E(Δ²P_ms)：表示临界群驱动功率的方差；

分别表示系统稳定裕度对临界群驱动功率、临界群负荷功率、余下群驱动功率和余下群负荷功率的灵敏度系数的ν次方；

分别表示第i个失稳模式下，系统稳定裕度对临界群驱动功率、临界群负荷功率、余下群驱动功率和余下群负荷功率的灵敏度系数；

ΔP_ms(i)、ΔP_ls(i)、ΔP_ma(i)、ΔP_la(i)：表示第i个失稳模式下，P_ms0(i)P_ls0(i)、P_ma0(i)和P_la0(i)与系统中各节点功率等于其概率分布均值时的临界群驱动功率

临界群负荷功率

余下群驱动功率

和余下群负荷功率

的偏差；

P′_ms：表示发电机G2驱动功率变化导致系统失稳模式改变时驱动功率的临界值，P′_ms的具体取值要根据实施例的情况确定。

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求范围当中。

Claims (17)

1.一种暂态稳定裕度概率分布的计算方法，其特征在于，所述方法包括如下步骤：

（1）确定初始状态下受扰系统的失稳模式和稳定裕度初始值η₀，并由用户确定采用的计算模式，包括基于纯解析法的暂态稳定裕度概率分布计算和综合时域仿真与解析法的暂态稳定裕度概率分布计算；

所述基于纯解析法的暂态稳定裕度概率分布计算包括如下步骤：

（2）计算所述受扰系统稳定裕度对临界群S驱动功率的灵敏度系数I；

（3）计算所述临界群S驱动功率的概率数字特征；

（4）计算所述受扰系统稳定裕度的各阶中心矩；

（5）计算所述受扰系统稳定裕度的累积分布函数；

所述综合时域仿真和解析法的暂态稳定裕度概率分布计算包括如下步骤：

（6）计算确定节点功率变化范围内的失稳模式变化情况；

（7）对所述各失稳模式计算稳定裕度对临界群驱动功率P_ms、临界群负荷功率P_ls、余下群驱动功率P_ma、余下群负荷功率P_la的灵敏度系数II；

（8）对所述各失稳模式计算该失稳模式下稳定裕度线性变化对应的各节点功率等于其概率分布均值时的稳定裕度初始值η′_0(i)；

（9）对所述失稳模式计算临界群驱动功率的概率数字特征、临界群负荷功率的概率数字特征、余下群驱动功率的概率数字特征和余下群负荷功率的概率数字特征；

（10）对所述各失稳模式计算稳定裕度的各阶中心矩；

（11）通过分段函数求取稳定裕度的累积分布函数。

2.如权利要求1所述的计算方法，其特征在于，所述概率数字特征即各阶半不变量γ_υ。

3.如权利要求1所述的计算方法，其特征在于，所述步骤（1）是令所述受扰系统中各节点功率等于其概率分布均值，对所述受扰系统中的设定故障，利用详细时域仿真和扩展等面积准则（EEAC），得到所述受扰系统的失稳模式和稳定裕度初始值η₀。

4.如权利要求1所述的计算方法，其特征在于，所述步骤（2）是令所述临界群S的驱动功率P_ms和余下群A的驱动功率P_ma的初始值分别等于各自群内机组驱动功率均值的总和，通过静态EEAC（SEEAC）算法求取所述受扰系统的稳定裕度对所述临界群S驱动功率P_ms的灵敏度系数I；所述灵敏度系数I包括一阶灵敏度系数

和二阶灵敏度系数

5.如权利要求1所述的计算方法，其特征在于，所述步骤（3）是根据所述临界群S内每台发电机组驱动功率P_mk的概率分布，求取每台发电机驱动功率P_mk的各阶半不变量γ_υ(P_mk)，按所述临界群S中发电机组个数对所述各阶半不变量γ_υ(P_mk)求和；求取P_ms的各阶半不变量γ_υ(P_ms)，为

${\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ms}}\right)=\underset{k=1}{\overset{\mathrm{Ns}}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{mk}}\right)$

其中，k为机组编号，υ表示半不变量的阶数，Ns为临界群机组个数。

6.如权利要求1所述的计算方法，其特征在于，所述步骤（4）是用所述灵敏度系数I和所述各阶半不变量γ_υ(P_ms)，计算所述受扰系统稳定裕度的各阶半不变量γ_υ(η)，通过半不变量和中心矩的标准数学转换公式，求取所述受扰系统稳定裕度的各阶中心矩β_υ(η)。

7.如权利要求1所述的计算方法，其特征在于，所述步骤（5）是将所述受扰系统稳定裕度的初始值η₀及其各阶中心矩β_υ(η)代入Gram-Charlier级数计算公式，得到所述受扰系统稳定裕度η的累积分布函数。

8.如权利要求1所述的计算方法，其特征在于，所述步骤（6）是根据所述受扰系统中具有功率不确定性的发电机和负荷功率变化范围的上下限，计算出对应步骤（1）初始状态下受扰系统失稳模式的临界群驱动功率P_ms、余下群驱动功率P_ma、临界群负荷功率P_ls、余下群负荷功率P_la变化范围的上下限，分别按照增加所述临界群S驱动功率和余下群A负荷功率，同时减少余下群A驱动功率和临界群S负荷功率，以及减少所述临界群S驱动功率和余下群A负荷功率，同时增加余下群A驱动功率和临界群S负荷功率两个方向，依据设定的功率变化步长及其在各机组和负荷之间的分配规则，生成一组系统稳态运行工况，再针对每一个系统稳态运行工况，进行针对设定故障下的时域仿真，识别失稳模式变化，记录每一个新的失稳模式i所对应的稳定裕度η′_(i)，以及临界群S的驱动功率P_ms0(i)、临界群S的负荷功率P_ls0(i)、余下群A的驱动功率P_ma0(i)和余下群A的负荷功率P_la0(i)。

9.如权利要求1所述的计算方法，其特征在于，所述步骤（7）是对每个失稳模式，采用摄动法，通过时域仿真和EEAC求取稳定裕度对所述临界群驱动功率P_ms的灵敏度系数

临界群负荷功率P_ls的灵敏度系数

余下群驱动功率P_ma的灵敏度系数

和余下群负荷功率P_la的灵敏度系数

10.如权利要求1所述的计算方法，其特征在于，所述步骤（8）中若所述受扰系统的失稳模式的个数为1，则稳定裕度初始值η₀₍₁₎=η₀；若所述受扰系统的失稳模式有变化，则根据稳定裕度的线性化展开式，结合步骤（7）得到的灵敏度系数II和步骤（6）得到的失稳模式变化后稳定裕度η′_(i)，利用数学反推获得每个失稳模式下的稳定裕度初始值η′_0(i)。

11.如权利要求1所述的计算方法，其特征在于，所述步骤（9）是根据所述受扰系统内发电机组驱动功率和负荷功率的概率分布，求取单个节点的发电机组驱动功率和负荷功率的各阶半不变量γ_υ，再根据每个失稳模式的分群情况，分别按照所述临界群S及余下群A内机组和负荷节点的个数分类求和，对应得到所述临界群驱动功率P_ms的各阶半不变量、临界群负荷功率P_ls的各阶半不变量、余下群驱动功率P_ma的各阶半不变量和余下群负荷功率P_la的各阶半不变量。

12.如权利要求1所述的计算方法，其特征在于，所述步骤（10）是用步骤（7）得到所述灵敏度系数II和步骤（9）得到的所述临界群驱动功率P_ms的各阶半不变量、临界群负荷功率P_ls的各阶半不变量、余下群驱动功率P_ma的各阶半不变量和余下群负荷功率P_la的各阶半不变量，计算所述稳定裕度η的各阶半不变量，进而通过半不变量和中心矩的标准数学转换公式，求取稳定裕度η的各阶中心矩β_υ(η)。

13.如权利要求1所述的计算方法，其特征在于，所述步骤（11）是针对每个失稳模式，将其稳定裕度的各阶中心矩β_υ(η)和稳定裕度初始值η′_0(i)代入Gram-Charlier级数计算公式，获得各失稳模式下稳定裕度η的累积分布函数，并利用分段函数综合描述。

14.如权利要求1所述的计算方法，其特征在于，所述步骤（4）中，在计算稳定裕度η的方差D(η)时，表达式如下：

$D\left(\mathrm{\η}\right)=S_{P_{\mathrm{ms}}}^{2}E\left({\mathrm{\Δ}}^2{P}_{\mathrm{ms}}\right)+\frac{1}{4}{\left[S_{P_{\mathrm{ms}}}^{\left(2\right)}E\left({\mathrm{\Δ}}^2{P}_{\mathrm{ms}}\right)\right]}^2---\left(1\right)$

15.如权利要求1所述的计算方法，其特征在于，所述步骤（4）中，按式(2)求取稳定裕度η对临界群驱动功率P_ms的一阶灵敏度系数按式(3)求取稳定裕度的各阶半不变量γ_υ(η)；

$S_{P_{\mathrm{ms}}}=\frac{\∂\mathrm{\η}}{{\∂P}_{\mathrm{ms}}}+\frac{\∂\mathrm{\η}}{{\∂P}_{\mathrm{ma}}}\frac{{\mathrm{dP}}_{\mathrm{ma}}}{{\mathrm{dP}}_{\mathrm{ms}}}---\left(2\right)$

${\mathrm{\γ}}_{\mathrm{\υ}}\left(\mathrm{\η}\right)=S_{P_{\mathrm{ms}}}^{\mathrm{\υ}}{\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ms}}\right)---\left(3\right)$

16.如权利要求1所述的计算方法，其特征在于，所述步骤（10）中，按式(4)根据时域仿真和EEAC获得的灵敏度系数

灵敏度系数

灵敏度系数

和灵敏度系数求取稳定裕度η的各阶半不变量γ_υ(η)；

${\mathrm{\γ}}_{\mathrm{\υ}}\left(\mathrm{\η}\right)=\left[\begin{array}{cccc}{S}_{P_{\mathrm{ms}}}^{\mathrm{\υ}}& S_{P_{\mathrm{ls}}}^{\mathrm{\υ}}& S_{P_{\mathrm{ma}}}^{\mathrm{\υ}}& S_{P_{\mathrm{la}}}^{\mathrm{\υ}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ms}}\right)\\ {\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ls}}\right)\\ {\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{ma}}\right)\\ {\mathrm{\γ}}_{\mathrm{\υ}}\left(P_{\mathrm{la}}\right)\end{array}\right]---\left(4\right)$

17.如权利要求1所述的计算方法，其特征在于，所述步骤（8)中，对每个失稳模式i，根据稳定裕度的线性化展开关系，按式(5)进行数学反推，求取各节点功率等于其概率分布均值时的稳定裕度初始值η′_0(i)。

${\mathrm{\η}}_{\left(i\right)}^{\′}={\mathrm{\η}}_{0\left(i\right)}^{\′}+\left[\begin{array}{cccc}{S}_{P_{\mathrm{ms}}\left(i\right)}& S_{P_{\mathrm{ls}\left(i\right)}}& S_{P_{\mathrm{ma}}\left(i\right)}& S_{P_{\mathrm{la}\left(i\right)}}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ΔP}}_{\mathrm{ms}\left(i\right)}\\ {\mathrm{\ΔP}}_{\mathrm{ls}\left(i\right)}\\ {\mathrm{\ΔP}}_{\mathrm{ma}\left(i\right)}\\ {\mathrm{\ΔP}}_{\mathrm{la}\left(i\right)}\end{array}\right]---\left(5\right)$

式中 ${\mathrm{\ΔP}}_{\mathrm{ms}\left(i\right)}=P_{\mathrm{ms}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{ms}},$ ${\mathrm{\ΔP}}_{\mathrm{ls}\left(i\right)}=P_{\mathrm{ls}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{ls}},$ ${\mathrm{\ΔP}}_{\mathrm{ma}\left(i\right)}=P_{\mathrm{ma}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{ma}},$ ${\mathrm{\ΔP}}_{\mathrm{la}\left(i\right)}$ $P_{\mathrm{la}0\left(i\right)}-{\stackrel{\‾}{P}}_{\mathrm{la}},$

和分别为各节点功率等于其概率分布均值时的临界群驱动功率、临界群负荷功率、余下群驱动功率和余下群负荷功率。

Images