基于受扰轨迹动态特性的数值积分提前终止方法与流程

本发明涉及一种基于受扰轨迹动态特性的数值积分提前终止方法，属于电力系统及其自动化技术领域。

背景技术：

电力系统暂态稳定分析，离不开数值积分。基于小步长逐步积分求得的受扰轨迹可以精确地反映复杂系统的非线性和非自治因素对暂态过程的影响。

但数值积分法本身并不能提供失稳的充要条件。因此，若仅获得受扰轨迹，电力系统专家只能依靠直觉和经验来判断所得到的受扰轨迹是否稳定，但却不知道该轨迹离开临界条件有多远。换言之，只能通过检查受扰轨迹是否超过经验限值来定性地判断该系统是否失稳。为了避免将任何稳定轨迹误判为失稳，必须将该限值取得很大。

扩展等面积准则(eeac)通过保稳的互补群惯量中心-相对运动(cccoi-rm)变换，将rⁿ多机系统的稳定性充要条件保留在某一个r¹等值映象系统中，从而将原高维系统稳定性的定性或定量分析问题严格地转换为多个映象平面轨迹的数据挖掘问题。eeac算法从受扰轨迹，而不是原始的数学模型出发，适用于任意可积的高维、强非自治、强非线性的运动系统，当然，其前提是在求取受扰轨迹时采用了相应的数学模型。

可见，数值积分仍是eeac算法的基础和前提。换言之，由数值积分给出多长积分区间的受扰轨迹，基于eeac算法亦仅能准确表征这段受扰轨迹的量化稳定裕度。

数值积分的时间长度，都依靠专家经验的事前确定。针对某一固定的积分区间，例如3s：对于大多研究算例，若受扰轨迹在这段积分区间内没有明显发散，可被定性地判定为暂态稳定；而对于若干其他算例，如缓慢发展的多摆失稳情形，则很可能被漏检。换言之，一个确定的积分长度对大多研究算例而言可能已经太长，但对某些算例可能却远远不够。为了尽可能地保证暂稳准确判定，就必须将积分区间设置得较大(如10s)，这样不论对失稳算例还是稳定算例，计算量都将大大增加，却仍旧未能免除误判的可能。

技术实现要素：

本发明的目的在于克服现有技术中的不足，提供一种基于受扰轨迹动态特性的数值积分提前终止方法，以解决现有技术中导致的上述多项缺陷或缺陷之一。

为达到上述目的，本发明是采用下述技术方案实现的：基于受扰轨迹动态特性的数值积分提前终止方法，所述方法包括：根据受扰轨迹信息，结合包括：失稳模式易变度、受扰轨迹非同调度、最大转子角间隙、同向相邻摆次转子角间隙信息相似度及振荡衰减度在内的动态特性，判断是否提前终止数值积分。

进一步的，所述方法具体包括如下步骤：

对测试算例全集中的每个算例分别执行下述步骤：

步骤a：求取算例的失稳模式易变度k；

步骤b：采用小步长逐步积分方法向前积分一步，直至到达主导映象首摆最远点或动态鞍点；

步骤c：若采用小步长逐步积分已达主导映象首摆最远点，则根据下述判定条件判定是否提前终止积分：

1)失稳模式易变度k值小于第一阈值ε1；

2)主导映象机械功率pm大于电磁功率分量pc；

若同时满足上述两个判定条件则提前终止积分；否则，执行下一步；

步骤d：若采用小步长逐步积分已达主导映象首摆动态鞍点，则根据下述判定条件判定是否提前终止积分：

1)失稳模式易变度k值小于第一阈值ε1；

若满足上述判定条件则提前终止积分；否则，执行下一步；

步骤e：基于首摆暂态过程受扰轨迹，求取首摆非同调度μ1；

步骤f：若采用小步长逐步积分已达主导映象首摆最远点，则根据下述判定条件判定是否提前终止积分：

1)失稳模式易变度k值小于第二阈值ε2；

2)主导映象机械功率pm大于电磁功率分量pc；

3)首摆非同调度μ1值小于第三阈值ε3；

若同时满足上述三个判定条件则提前终止积分；否则，执行下一步；

步骤g：若采用小步长逐步积分已达主导映象首摆动态鞍点，则根据下述判定条件判定是否提前终止积分：

1)失稳模式易变度k值小于第二阈值ε2；

2)首摆非同调度μ1值小于第三阈值ε3；

若同时满足上述两个判定条件则提前终止积分；否则，执行下一步；

步骤h：采用小步长逐步积分方法向前积分一步，直至最大转子角间隙大于预设的第四阈值ε4，或到达主导映象第三摆最远点或动态鞍点；

步骤i：若最大转子角间隙已大于第四阈值ε4，则提前终止积分；否则，执行下一步；其中：δij表示机组i转子角同机组j转子角的差值，即：δij＝δi-δj；i表示领前群中机组i；j表示余下群中机组j；s表示领前群；a表示余下群；

步骤j：若采用小步长逐步积分已达主导映象第三摆最远点或动态鞍点，求取前三摆非同调度μ1-3；若前三摆非同调度μ1-3值小于第五阈值ε5，则提前终止积分；否则，执行下一步；

步骤k：采用小步长逐步积分方法向前积分一步，直至最大转子角间隙大于预设的第四阈值ε4，或到达主导映象下一摆最远点；

步骤l：若最大转子角间隙已大于第四阈值ε4，则提前终止积分；否则，执行下一步；

步骤m：若采用小步长逐步积分已达主导映象下一摆最远点，从该摆次始，向后溯洄四个摆次，若已小于第四摆，执行步骤k，否则执行下一步；

步骤n：根据步骤m中向后溯洄四个摆次的暂态过程，求取正向及反向相邻摆次转子角间隙信息相似度，若相似度小于第六阈值ε6，则提前终止积分；否则，执行下一步；

步骤o：若从该摆次始，向后溯洄六个摆次，若已小于第四摆，则执行步骤k，否则执行下一步；

步骤p：根据步骤o中向后溯洄的六个摆次的暂态过程，求取连续两个正向及反向相邻摆次振荡衰减度，若振荡衰减度小于第七阈值ε7，则提前终止积分；否则，执行步骤k。

进一步的，所述失稳模式易变度k采用大步长泰勒级数展开快速求取。

进一步的，首摆非同调度μ1的求取方法包括如下步骤：

μ1＝max{μs.d,μs.p,μa.d,μa.p}(1)

式(1)中，μs.d、μs.p、μa.d、μa.p分别为首摆故障中领前群、故障后领前群、故障中余下群、故障后余下群非同调度指标，计算公式分别如式(2)～(5)所示：

式(2)～(5)中，m为惯量，ω为角速度；下标s、a分别代表领前群等值机、余下群等值机；下标d、p、τ分别代表故障中、故障后、故障清除瞬间。

进一步的，第三阈值ε3的取值为：

式(6)中，k为研究系统所含机组数。

进一步的，前三摆非同调度μ1-3的求取方法如下：

μ1-3＝max{μ1,μ2,μ3}(7)

式(7)中，μ2、μ3分别为第2摆、第3摆非同调度。

进一步的，第2摆非同调度μ2的计算方法如下：

μ2＝max{μs.2,μa.2}(8)

式(8)中，μs.2、μa.2分别为第二摆领前群、余下群非同调度指标，计算方法分别如式(9)、(10)所示：

式(9)、(10)中，m为惯量，ω为角速度；下标s、a分别代表领前群等值机、余下群等值机；下标fep1、dcp2分别代表到达主导映象的首摆最远点和第二摆动态中心点时刻。

进一步的，第五阈值ε5的取值为：

式(11)中，k为研究系统所含机组数。

进一步的，正向相邻摆次转子角间隙信息相似度的计算方法如下：

式(12)中，d2l-1、d2l+1分别为到达主导映象第(2l-1)、(2l+1)摆最远点时刻各机转子角间隙值的标准差。

进一步的，连续两个正向相邻摆次振荡衰减度的计算方法如下：

式(13)中，δ为主导映象系统等值转子角，下标fep2l-1代表到达主导映象第(2l-1)摆最远点时刻；特别地，fep0代表稳态平衡点时刻。

综上，本发明提供的基于受扰轨迹动态特性的数值积分提前终止方法，从稳定性机理的多个视角出发执行数值积分终止的提前判别，由此在确保定性分析的正确性和定量分析精度的前提下进一步大大减少了计算量。对于解决高比例可再生能源入网情形下的暂态稳定分析所面临的强不确定性、强非线性及维数灾等问题具有重大的理论和工程意义。

附图说明

图1是根据本发明实施例的流程图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

应当理解：

当且仅当多机电力系统采用经典模型，而且动态轨迹具有理想两群特性时，其主导映象是定常哈密顿系统。理想的两群特性是指整个动态过程中，同群各发电机的加速度一直保持一致。此时，该系统的失稳模式只有首摆失稳和第2摆失稳两种形式，既不可能发生多摆失稳，也不可能有不同的失稳映象。

不论是互补群各群内部的非同调因素，还是多机运动方程中的非哈密顿因素，对于r¹映象来说都可被视为在相应的单自由度哈密顿系统上的扰动。设系统在有限的观察时段内稳定，如果要据此而得到“系统在该时段后仍然保持长期稳定”的断言，则必须保证主导映象在观察期以后不再有强时变因素。

主导映象当前的稳定裕度越大，且时变因素越弱，则观察期以后保持稳定的可能性越大。显然，多机动态过程的两群特性越接近于理想，则时变因素就越弱。

由此，量化时变因素的途径之一即为衡量多机动态过程的两群特性，换言之，即两群内各机组的非同调特性。

不过，群内非同调特性并非衡量时变因素强弱的唯一途径。稳定性机理分析和大量算例仿真一致表明：对于时变程度较弱的算例，算例失稳模式不易随扰动参数(如故障清除时间τ值)变化；对于时变程度较强的算例，算例失稳模式易随扰动参数(如故障清除时间τ值)变化。由此，失稳模式易变度指标亦可用于指导数值积分的提前终止。

至此，从时变因素影响稳定裕度随积分区间敏感程度的因果分析的角度出发，可利用失稳模式易变度、受扰轨迹非同调度两指标指导数值积分的提前终止。

另一方面，随着暂稳态势的推演，在后继摆动期间，若各机转子角间隙信息保持相对一致或主导映象系统呈振荡衰减情形，则亦可得到“系统在该时段后仍然保持长期稳定”的结论。由此，同向相邻摆次转子角间隙信息相似度及振荡衰减度指标亦可用于指导数值积分的提前终止。

上述四个指标中，失稳模式易变度指标由大步长泰勒级数展开快速求得，其余三个指标均为利用数值积分求取受扰轨迹的中间结果。由此，该方法所需计算量的绝大部分仍为受扰轨迹的求取，执行提前判稳所增加的计算量小到几乎无法统计。这样，对各算例执行暂稳量化分析所需的积分区间依赖于各自的动态特性，而不是盲目地取为定值。

本实施例即根据研究算例时变程度影响稳定裕度对积分区间敏感程度的因果分析，提供了一种基于受扰轨迹动态特性的数值积分提前终止方法，为含时变因素强弱程度不同的各算例匹配各不相同的积分区间，以期进一步兼顾暂稳分析精度与速度。具体而言，本实施例的步骤如图1所示：

图1中步骤1描述的是暂态稳定详细分析流程启动后，取出测试算例全集中的某个算例，利用大步长泰勒级数展开算法快速求取该算例失稳模式易变度k。具体说来，针对该算例，设置其故障清除时间为0.8s，并由大步长泰勒级数展开获得该算例故障中、故障后各机转子角轨迹；选取n个时刻作为故障清除后各机组暂态稳定态势发展较为/完全明显的n个时间断面，分别观察每一时间断面的各机组转子角，获取各时间断面最大、次最大转子角间隙，并由此求得表征各时间断面各机组转子角间隙信息的指标k1、k2、…、kn；以各时间断面各机组转子角间隙信息指标k1、k2、…、kn中的最大值k表征该算例失稳模式易变度，即k＝max{k1,k2,l,kn}。

图1中步骤2描述的是应用利用小步长逐步积分技术向前积分一步，直至到达主导映象首摆最远点(fep)或动态鞍点(dsp)。

图1中步骤3描述的是决定是否提前终止数值积分的判定规则：若基于小步长逐步积分已达主导映象首摆fep，同时由步骤1反映的失稳模式易变度k值小于阈值ε1；进一步地，首摆暂态过程中、主导映象机械功率pm大于电磁功率分量pc，则判定可提前终止积分，执行步骤18，否则执行步骤4。

图1中步骤4描述的是决定是否提前终止数值积分的判定规则：若基于小步长逐步积分已达主导映象首摆dsp，同时由步骤1反映的失稳模式易变度k值小于阈值ε1，则判定可提前终止积分，执行步骤18，否则执行步骤5。

步骤3和步骤4中，ε1为静态阈值，它是根据不同实际系统的大量典型算例，以可靠性为首要原则优化求得的，针对不同系统不同工况中的各算例具有强壮性，在不同系统、工况及算例下都不变。阈值ε1的值为0.095。

图1中步骤5披露了一种反映首摆暂态过程中受扰轨迹非同调度的计算方法：通过分别求取首摆故障中领前群、故障后领前群、故障中余下群、故障后余下群非同调度指标μs.d、μs.p、μa.d、μa.p，按式(1)取四者中的最大值来表征首摆非同调度μ1：

μ1＝max{μs.d,μs.p,μa.d,μa.p}(1)

指标μs.d、μs.p、μa.d、μa.p则按式(2)～(5)求取：

式(2)～(5)中，m为惯量，ω为角速度；i表示领前群中机组i；j表示余下群中机组j；s表示领前群；a表示余下群；下标s、a分别代表领前群等值机、余下群等值机；下标d、p、τ分别代表故障中、故障后、故障清除瞬间。

图1中步骤6描述的是决定是否提前终止数值积分的判定规则：若基于小步长逐步积分已达主导映象首摆fep，同时由步骤1反映的失稳模式易变度k值小于阈值ε2；此外首摆暂态过程中、主导映象机械功率pm大于电磁功率分量pc；进一步地，由步骤5反映的首摆非同调度μ1值小于阈值ε3，则判定可提前终止积分，执行步骤18，否则执行步骤7。

图1中步骤7描述的是决定是否提前终止数值积分的判定规则：若基于小步长逐步积分已达主导映象首摆dsp，同时由步骤1反映的失稳模式易变度k值小于阈值ε2；此外由步骤5反映的首摆非同调度μ1值小于阈值ε3，则判定可提前终止积分，执行步骤18，否则执行步骤8。

步骤6和步骤7中ε2为静态阈值、ε3为动态阈值，它们都是根据不同实际系统的大量典型算例，以可靠性为首要原则优化求得的，针对不同系统不同工况中的各算例具有强壮性，在不同系统、工况及算例下都不变。阈值ε1、ε3分别按式(6)、(7)取值：

ε2＝0.448(6)

式(7)中，k为研究系统所含机组数。

图1中步骤8描述的是利用小步长逐步积分技术向前积分一步，直至最大转子角间隙大于某阈值，或到达主导映象第3摆fep或dsp。

图1中步骤9描述的是决定是否提前终止数值积分的判定规则：若最大转子角间隙已大于阈值ε4，则判定可提前终止积分，执行步骤18，否则执行步骤10。

该步骤中，阈值ε4为静态阈值，它是根据不同实际系统的大量典型算例，以可靠性为首要原则优化求得的，针对不同系统不同工况中的各算例具有强壮性，在不同系统、工况及算例下都不变。阈值ε4的值为500°。

图1中步骤10披露了一种反映前三摆暂态过程中受扰轨迹非同调度的计算方法：通过分别求取首摆、第2摆、第3摆非同调度μ1、μ2、μ3，按式(8)取三者中的最大值来表征前三摆非同调度μ1-3：

μ1-3＝max{μ1,μ2,μ3}(8)

按式(1)求取指标μ1，按式(9)求取指标μ2：

μ2＝max{μs.2,μa.2}(9)

式(9)中，μs.2、μa.2分别为第2摆领前群、余下群非同调度指标，可分别按式(10)、(11)定义之：

式(10)、(11)中，m为惯量，ω为角速度；下标i、j分别代表领前群、余下群某机组i、j；下标fep1、dcp2分别代表到达主导映象的首摆fep和第二摆dcp时刻。

指标μ3的求取则可依μ2的求取步骤执行。

图1中步骤11描述的是决定是否提前终止数值积分的判定规则：若由步骤10反映的前三摆非同调度μ1-3值小于阈值ε5，则判定可提前终止积分，执行步骤18，否则执行步骤12。

该步骤中，阈值ε5为动态阈值，它是根据不同实际系统的大量典型算例，以可靠性为首要原则优化求得的，针对不同系统不同工况中的各算例具有强壮性，在不同系统、工况及算例下都不变。阈值ε5的按式(12)取值。

式(12)中，k为研究系统所含机组数。

图1中步骤12描述的是利用小步长逐步积分技术向前积分一步，直至最大转子角间隙大于某阈值，或到达主导映象下一摆fep。

图1中步骤13描述的是决定是否提前终止数值积分的判定规则：若最大转子角间隙已大于阈值ε4，则判定可提前终止积分，执行步骤18，否则执行步骤14。

该步骤中，阈值ε4为静态阈值，它同步骤9中阈值ε4一致。

图1中步骤14描述的是，若基于小步长逐步积分已达主导映象下一摆fep，从该摆次始，向后溯洄四个摆次，若已小于第4摆，执行步骤12，否则执行步骤15。

图1中步骤15披露了一种正向及反向相邻摆次转子角间隙信息相似度(s)的计算方法，以及决定是否提前终止数值积分的判定规则：

通过比较主导映象正向某两个相邻摆次到达各自fep时刻、各机转子角间隙值的标准差的差异，按式(13)求取正向相邻摆次转子角间隙信息相似度(s)：

式(13)中，d2l-1、d2l+1分别为到达主导映象第2l-1、2l+1摆fep时刻，各机转子角间隙值的标准差。

反向相邻摆次转子角间隙信息相似度(s)可依此类推。

若正向/反向相邻摆次转子角间隙信息相似度s值小于阈值ε6，则判定可提前终止积分，执行步骤18，否则执行步骤16。

该步骤中，阈值ε6为静态阈值，它是根据不同实际系统的大量典型算例，以可靠性为首要原则优化求得的，针对不同系统不同工况中的各算例具有强壮性，在不同系统、工况及算例下都不变。阈值ε6的值为0.100。

图1中步骤16描述的是，若基于小步长逐步积分已达主导映象下一摆fep，从该摆次始，向后溯洄六个摆次，若已小于第4摆，执行步骤12，否则执行步骤17。

图1中步骤17披露了一种连续两个正向及反向相邻摆次振荡衰减度(e)的计算方法，以及决定是否提前终止数值积分的判定规则：

通过比较主导映象正向三个相邻摆次的等值转子角行程间的两两差异，按式(14)求取连续两个正向相邻摆次振荡衰减度(e)：

式(14)中，δ为主导映象系统等值转子角，下标fep2l-1代表到达主导映象第2l-1摆fep时刻；特别地，fep0代表初始(稳态平衡点)时刻。

连续两个反向相邻摆次振荡衰减度(e)的定义可依此类推。

若连续两个正向及反向相邻摆次振荡衰减度e值小于阈值ε7，则判定可提前终止积分，执行步骤18，否则执行步骤12。

该步骤中，阈值ε7为静态阈值，它是根据不同实际系统的大量典型算例，以可靠性为首要原则优化求得的，针对不同系统不同工况中的各算例具有强壮性，在不同系统、工况及算例下都不变。阈值ε7的值为0.800。

图1中步骤18描述的是，如果测试算例全集中每个算例都已完成数值积分，结束数值积分提前终止流程，否则取下一算例执行步骤1。

作为本实施例的具体计算，以ieee10g39n系统，以及海南(2009年数据)、山东(2004及2012年数据，分别记为山东a及山东b)、江西(2011年数据)、浙江(2012及2013年数据，分别记为浙江a及浙江b)、河南(2011年数据)、新疆(2012年数据)这9个系统为测试系统。

除ieee10g39n系统仅考虑经典模型外，其他各工程系统均计及发电机详细模型、励磁器、调速器等的影响。

测试算例为线路三相永久短路故障，故障地点大多选择研究系统的关键节点，故障清除时间则取大多取研究算例的(准)临界清除时间，并特别注意将发生多摆失稳、产生病态dsp(vdsp)、失先稳再同步情形的算例纳入测试算例中。

各系统算例数及各算例暂稳情形分别如表1、2所示，执行暂稳分析时、将观察时间设定为10s：

表1

表2

表3显示这357个测试算例经本发明所述技术方案提前终止积分情况：

表3

表3显示：

对于超过60％的测试算例，仅需基于主导映象首摆信息即可提前终止积分，可提前终止积分者，均为暂态稳定的算例；

对于超过80％的测试算例，仅需基于主导映象前三摆信息即可提前终止积分，可提前终止积分者，既包括暂态稳定的算例，也包括前三摆失稳的算例。

在剩下的不到20％的算例中，包括了如下算例：

13个多摆失稳的算例，在检测到转子角间隙值大于阈值后可终止积分。

1个先失稳再同步的算例，在判别其达到平稳振荡后(12摆次)可终止积分。

50个暂态稳定的算例，在判别其达到平稳振荡或振荡衰减状态后可终止积分。对于这50个算例，进一步的统计结果如表4所示：

表4

表4中，4个算例即使基于30s积分区间仍然无法提前终止积分，进一步仿真发现：对于这些算例，需至少积分至60s才可望满足平稳振荡条件以终止积分。而对于其他46个已能提前终止积分者，所采用方法的统计结果如表5所示：

表5

对于17个满足振荡衰减条件可提前终止积分者，亦可通过满足平稳振荡条件来提前终止积分，但积分时间需增加；对于29个满足平稳振荡条件可提前终止积分者，则无法利用基于振荡衰减的判据来提前终止积分。

应用本发明所述技术方案实际执行时，超过80％的测试算例仅需基于前三摆次信息即可实现提前终止积分；对基于前三摆次信息仍不能提前终止积分者，每向前积分一个摆次，均依次判别是否满足平稳振荡或振荡衰减条件，来决定可否提前终止积分，以提高效率。

进一步地，通过对测试算例进行校核，对于可提前终止数值积分的算例，分别比较根据已得轨迹执行数据挖掘所获的量化裕度信息，和根据固定积分区间(如10s)所得积分轨迹的量化裕度信息间的差异，可知其最大相对差值不超过0.01％。

综上，本发明可为各算例根据上述动态信息自动匹配恰当的积分区间，相较原基于固定积分区间的暂态稳定分析方法，可在确保算例分析精度不变的前提下，大大提高分析速度；对进一步兼顾在线暂态稳定分析精度和速度，解决高比例可再生能源入网情形下的暂态稳定分析可能面临的强不确定性、强非线性及维数灾问题，具有重大的理论和工程意义。

虽然本发明已以较佳实施例公开如上，但实施例并不是用来限定本发明的。在不脱离本发明之精神和范围内，所做的任何等效变化或润饰，同样属于本发明之保护范围。因此本发明的保护范围应当以本申请的权利要求所界定的内容为标准。