基于低阶EIGD的时滞电力系统稳定性分析方法与流程

本发明涉及电力系统稳定性分析技术领域，特别是涉及基于低阶显式无穷小生成元离散化(explicitlyinfinitesimalgeneratordiscretization,eigd)的时滞电力系统稳定性分析方法。

背景技术：

现代电力系统的运行和控制无时无刻不依赖于一个可靠的信息系统。基于计算机技术、通信技术和传感技术的电力信息系统与电力一次系统紧密而有机地结合在一起。电力系统在本质上是一个信息物理融合的动力系统，即cpps(cyber-physicalpowersystem)。广域测量系统(wide-areameasurementsystem，wams)的出现给大规模互联电力系统稳定分析与控制的发展带来新的契机，其高级应用包括：电力系统动态监测与状态估计、参数辨识、稳定性监测与评估、低频振荡辨识和广域阻尼控制、故障定位与广域保护等。

20多年以来，国内外学者针对基于广域测量信息的电力系统分析与控制开展了大量研究。作为wams的一项高级应用，基于广域测量信息的广域阻尼控制，通过引入有效反映区间低频振荡模式的广域反馈信号，如发电机相对转速和功角、联络线功率等，能够显著增强对制约大规模互联电网输电能力的区间低频振荡问题的控制能力，进而为提高系统的输电能力提供了新的控制手段，具有良好而又广泛的应用前景。广域测量信号在采集、路由、传输和处理过程中存在数十到几百毫秒的时延。广域阻尼控制回路中的通信时滞，对控制器的性能产生影响并为电力系统带来运行风险。

基于谱离散化的特征值分析方法(spectraldiscretization-basedeigen-analysismethods,sdms)可有效分析时滞对cpps的小干扰稳定性的影响。其原理是利用谱算子将时滞cpps(delayedcpps,dcpps)的谱映射为状态转移方程和抽象柯西方程的谱，特征方程中不存在超越的指数项。为表述方便，后续将“时滞信息物理融合的电力系统”简称为“时滞电力系统”。因此，大规模的多时滞电力系统的关键特征值可以通过传统特征分析法求解。

例如，中国发明专利申请号为201510055743.3，专利名称为：基于eigd的大规模时滞电力系统特征值计算方法，提出了一种基于显式无穷小生成元离散化(explicitinfinitesimalgeneratordiscretization,eigd)大规模时滞电力系统特征值计算；中国发明专利申请号为201510829129.8，专利名称为：基于sod-ps-r算法的时滞电力系统机电振荡模式计算方法，提出了一种基于解算子伪谱离散化(sod-ps)的大规模时滞电力系统特征值计算方法。

上述sdms方法的共有弊端是谱算子的离散化矩阵维数较大，导致特征值计算的计算量大，计算时间长。为解决这个问题，可以利用克罗内克积(kroneckerproudct)的性质减少特征值计算中的矩阵-向量乘积量。但是从本质上来说，这种方式并没有减低谱离散化矩阵的矩阵维数。显然，如果能够降低离散化矩阵的阶数，sdms方法的效率能够显著提高。

技术实现要素：

为了解决现有技术的不足，本发明提供了基于低阶eigd的时滞电力系统稳定性分析方法，能够降低离散化矩阵的阶数，提高计算效率。

基于低阶eigd的时滞电力系统稳定性分析方法，包括：

建立时滞电力系统的动态模型并线性化为一组时滞微分方程；

将时滞微分方程中时滞电力系统的状态变量重新排序，划分为与时滞无关的状态变量和与时滞有关的状态变量，在对过去时刻的时滞电力系统状态进行离散化时，剔除与时滞无关的状态变量的离散化，其相应的离散化矩阵的维数减少；

利用无穷小生成元将时滞微分方程转换为常微分方程，从而将时滞电力系统的特征方程等价转换为常规特征方程，常规特征方程通过无穷小生成元描述，进而将时滞电力系统的特征值计算转化为计算无穷小生成元的特征值；

由于无穷小生成元的特征值有无穷多个，需要对无穷小生成元进行离散化，首先，通过时滞区间上的多个离散点将时滞区间划分为离散函数空间，从而将连续函数离散化为分块向量，然后通过拉格朗日插值多项式近似离散点处的函数值，排列并改写这些离散点处的状态方程得到无穷小生成元的离散化矩阵，从而将无限维的特征值问题转化为有限维的特征值问题；

依据时滞电力系统的特征值与无穷小生成元特征值之间一一对应的谱映射关系，将计算时滞电力系统最右端的特征值近似为无穷小生成元离散化矩阵最右端的特征值；

将无穷小生成元离散化矩阵经过位移-逆变换，无穷小生成元最右端的特征值变换为模值最大的特征值，计算无穷小生成元离散化矩阵经过位移-逆变换后的矩阵模值最大的特征值，并将其反变换为系统最右端特征值的近似值；

对近似特征值进行修正，得到时滞电力系统的精确特征值，利用得到的精确特征值判断电力系统的稳定性。

进一步优选的技术方案，所述建立时滞电力系统的动态模型并线性化为一组时滞微分方程，具体为：将时滞电力系统模型在平衡点处进行泰勒展开，得到线性化后的系统动态模型，即通过一组时滞微分方程描述时滞电力系统的动态模型。

进一步优选的技术方案，所述通过时滞微分方程表示时滞电力系统模型，具体为：

式中：x∈r^n×1为电力系统的状态变量向量；n为系统状态变量总数；t为当前运行时刻；m为时滞个数；τ1,…,τm为时滞常数，τmax表示最大的时滞；δx(t)为t时刻系统状态变量的增量；δx(t-τi)为t-τi时刻系统状态变量的增量；为t时刻系统状态变量增量的导数；δx(0)为系统状态变量的初始值即初始状态，并简写为和分别表示稠密的系统状态矩阵和稀疏的时滞状态矩阵。

在后续分析中，和只具有理论分析的意义，为充分利用其稀疏性，计算时采用其增广形式，具体可以表示为：

式中：a0，b0，为高度稀疏的系统增广状态矩阵，l为系统代数变量总数。

式(1)所示系统的特征方程为：

式中：λ和分别表示dcpps的特征值及其对应的右特征向量。

进一步优选的技术方案，将时滞电力系统的状态变量重新排序，依据其时滞状态矩阵相应元素是否为0将其划分为与时滞无关的状态和与时滞相关的状态其中n1+n2＝n。

具体来说，如果在所有时滞状态矩阵中，某个状态变量对应的行、列元素均为0，则将该变量当做与时滞无关的状态，剩余的变量则为与时滞无关的状态。相应地，公式(1)的第一式可改写为：

式中，和分别表示照状态变量的顺序和划分对系统状态矩阵的改写并分块，均为稠密矩阵，和分别表示依照状态变量的顺序和划分对时滞状态矩阵进行改写并分块：

其中，011，012和021分别表示n1×n1、n1×n2和n2×n1维的零矩阵，表示定义，即将adi,22定义为adi；

故，式(4)写作

相应地，系统的特征方程可改写为：

由于式(7)中存在指数项系统的特征值有无穷个。

进一步优选的技术方案，设状态空间是由区间[-τmax,0]到n维实数空间映射的连续函数构成的巴拿赫(banach)空间，并赋有上确界范数巴拿赫空间是无穷维空间，解算子x→x用来表征了时滞电力系统初始状态和不同时刻的状态之间的关系。

具体地，其将θ(-τmax≤θ≤0)时刻时滞电力系统初始状态转移到θ+h(h≥0)时刻系统状态δx(h+θ)(＝δxh(θ))的线性有界算子，无穷小生成元x→x定义为解算子在h＝0时刻的导数，其中，c表示复数泛函空间。

将公式(1)的第一式代入无穷小生成元可得

由式(8)可以看出，时滞微分方程式(1)被转化为常规齐次线性微分方程，其系数即为无穷小生成元其状态为δxt，无穷小生成元的特征方程和特征函数分别为：

这表明，时滞电力系统的特征值就是无穷小生成元的特征值，即：

式中：σ(·)表示矩阵或算子的全部特征值集合，即谱，从式(11)可以看出，时滞电力系统的特征值计算可直接转化为无穷小生成元的特征值计算。

进一步优选的技术方案，给定任意正整数n，区间[-τmax,0]上n+1个离散点形成的集合ωn定义为ωn:＝{θn,j,j＝0,1,...,n}，且有-τmax≤θn,n≤…≤θn,1≤θn,0＝0，设表示集合ωn上定义的离散函数空间。

由于系统的初始状态是x上的连续函数，也即可以表示为通过离散函数空间可将离散化为分块向量其中离散函数是连续函数在离散点θn,j处函数值的近似，即设lnφ表示唯一存在的次数不超过n的拉格朗日插值多项式，且满足

式(1)对应的特征方程表示为

设配置多项式lnφ可由阶数等于或小于n的切比雪夫多项式的一组基表示，即：

式中：为常数向量，j＝0,…,n。其中，

tj(·)为第j阶第一类切比雪夫多项式，即tj(·)＝cos(jarccos(·))，arccos表示反余弦函数。

可通过以下方程来降阶表示dcpps的状态方程

式中，uj(·)为第j阶第二类切比雪夫多项式。

无穷小生成元的伪谱离散化近似矩阵xn→xn为：

式中，为伴随矩阵，为块上三角矩阵。πn和σn可具体表示为：

式中：和分别表示维数为n1×n1和n2×n2的单位矩阵。表示一个高度稀疏的伴随常数矩阵。rj(j＝1,…,n+1)由系统状态矩阵经过简单运算得到，即：

矩阵-与向量v相乘时，将γ′改写为的形式，充分利用kronecker积的性质减少运算量。其中，表示拉格朗日系数向量。

进一步优选的技术方案，采用隐式arnoldi算法计算系统最右端的特征值时，需要计算位移-逆矩阵模值最大(并依次递减)的部分关键特征值。(σ′n)^-1和σn具有完全相同的分块稀疏结构，而且的各个非零子块都可以被显式地表达出来，见式(19)：

γ′＝[in-r′1…-r′n](20)

式(19)中，表示维数为nn2×nn2的单位矩阵，式(21)-(23)中，s表示位移点。

采用牛顿迭代法对得到的系统近似特征值进行修正。

利用得到的系统准确特征值λ判定系统的稳定性：若存在λ的实部大于零，即re(λ)>0，系统不稳定。若所有re(λ)<0，系统渐进稳定。若存在re(λ)＝0，则系统临界稳定。其中，re(·)表示特征值的实部。

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

第一、本发明提出的基于低阶eigd的特征值计算方法用于系统的稳定性分析，充分考虑了时滞的影响。

第二、本发明通过减少与时滞无关状态变量的离散化，大大减少了显式无穷小生成元离散化矩阵的维数，从而在保证精算精度的同时，大大降低了dcpps关键特征值计算的计算量。尤其是对大规模dcpps，本发明提出的特征值计算效率与无时滞电力系统计算效率相当。从本质上提高了原有eigd方法的计算效率。

第三、本发明充分利用离散化矩阵和系统增广矩阵的稀疏性，减少本发明的特征值计算方法中矩阵-向量乘积的计算量。

附图说明

构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

图1为本发明的方法流程图；

图2为山东电网示意图；

表1为本发明提出的方法与eigd方法在选取不同的参数时的计算效率对比。

具体实施方式

应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

如图1所示：基于lowe-ordereigd的时滞电力系统稳定性分析方法，包括如下步骤：

步骤(1)：建立dcpps动态模型并对其进行线性化，系统状态可通过一组时滞微分方程描述，从而得到系统的状态矩阵和时滞状态矩阵，并依次给出dcpps的特征方程。

步骤(2)：将时滞电力系统的状态变量重新排序，并划分为与时滞无关的状态和与时滞有关的状态。随之，系统状态矩阵和时滞状态矩阵重新排列并分块。此时，时滞电力系统的状态方程分为两部分。由于存在指数项，系统的特征值有无穷个。

步骤(3)：利用无穷小生成元将时滞微分方程转换为常微分方程，从而将时滞电力系统的特征方程等价地转换为常规特征方程，进而将时滞电力系统的特征值计算转化为计算无穷小生成元的特征值。

步骤(4)：通过时滞区间上的多个离散点将时滞区间划分为离散函数空间，从而将连续函数离散化为分块向量。然后通过拉格朗日插值多项式近似离散点处的函数值，排列并改写这些离散点处的状态方程得到无穷小生成元的离散化矩阵，从而将无限维的特征值问题转化为有限维的特征值问题。

步骤(5)：依据时滞电力系统的特征值与无穷小生成元特征值之间的谱映射关系，将计算时滞电力系统最右端的特征值近似为无穷小生成元离散化矩阵最右端的特征值。

步骤(6)：采用隐式arnoldi算法计算步骤(4)得到的无穷小生成元离散化矩阵经过位移-逆变换后的矩阵模值最大的特征值，并将其反变换为系统最右端特征值的近似值。

步骤(7)：采用牛顿迭代法对步骤(6)中得到的系统近似特征值进行修正，从而得到时滞电力系统的精确特征值和特征向量。

步骤(8)：利用步骤(7)中得到的系统准确特征值判定系统的稳定性。

所述步骤(1)中，线性化的时滞电力系统模型可通过以下时滞微分方程描述：

式中：x∈r^n×1为电力系统的状态变量向量，n为系统状态变量总数。t为当前运行时刻。τ1,…,τm为时滞常数，τmax表示最大的时滞。δx(t)为t时刻系统状态变量的增量；δx(t-τi)为t-τi时刻系统状态变量的增量；为t时刻系统状态变量增量的导数。δx(0)为系统状态变量的初始值(即初始条件)，并简写为和分别表示稠密的系统状态矩阵和稀疏的时滞状态矩阵，其增广形式为

式中：为高度稀疏的系统增广状态矩阵。l为系统代数变量总数。

式(1)所示系统的特征方程为：

式中：λ和分别表示dcpps的特征值及其对应的右特征向量。

所述步骤(2)中，将时滞电力系统的状态变量重新排序，并将其划分为与时滞无关的状态和与时滞相关的状态其中n1+n2＝n。相应地，步骤(1)中公式(1)的第一式可改写为：

式中，和分别表示依照状态变量的顺序和划分对时滞状态矩阵进行改写并分块。

其中，011，012和021分别表示n1×n1、n1×n2和n2×n1维的零矩阵。表示定义，即将adi,22定义为adi。

故，式(4)可写作

相应地，系统的特征方程可改写为：

由于式(7)中存在指数项系统的特征值有无穷个。

所述步骤(3)中，设状态空间是由区间[-τmax,0]到n维实数空间映射的连续函数构成的巴拿赫(banach)空间，并赋有上确界范数巴拿赫空间是无穷维空间。解算子x→x用来表征了时滞电力系统初始状态和不同时刻的状态之间的关系。具体地，其将θ(-τmax≤θ≤0)时刻时滞电力系统初始状态转移到θ+h(h≥0)时刻系统状态δx(h+θ)(＝δxh(θ))的线性有界算子。无穷小生成元x→x定义为解算子在h＝0时刻的导数。将步骤(1)中公式(1)的第一式代入无穷小生成元可得

由式(8)可以看出，时滞微分方程式(1)被转化为常规齐次线性微分方程，其系数即为无穷小生成元其状态为δxt。无穷小生成元的特征方程和特征函数分别为：

这表明，时滞电力系统的特征值就是无穷小生成元的特征值，即：

式中：σ(·)表示矩阵或算子的全部特征值集合，即谱。从式(11)可以看出，时滞电力系统的特征值计算可直接转化为无穷小生成元的特征值计算。

所述步骤(4)中，给定任意正整数n，区间[-τmax,0]上n+1个离散点形成的集合ωn定义为ωn:＝{θn,j,j＝0,1,...,n}，且有-τmax≤θn,n≤…≤θn,1≤θn,0＝0。设表示集合ωn上定义的离散函数空间。因此，x上的任意连续函数被离散化为分块向量其中离散函数是连续函数在离散点θn,j处函数值的近似，即设lnφ表示唯一存在的次数不超过n的拉格朗日插值多项式，且满足

步骤(1)中式(1)对应的特征方程可表示为

设配置多项式lnφ可由阶数等于或小于n的切比雪夫多项式的一组基表示，即：

式中：为常数向量，j＝0,…,n。其中，tj(·)为第j阶第一类切比雪夫多项式，即tj(·)＝cos(jarccos(·))，arccos表示反余弦函数。

可通过以下方程来降阶表示dcpps的状态方程

式中，uj(·)为第j阶第二类切比雪夫多项式。

所述步骤(4)中，无穷小生成元的伪谱离散化近似矩阵xn→xn为：

式中，为伴随矩阵，为块上三角矩阵。πn和σn可具体表示为：

式中：和分别表示维数为n1×n1和n2×n2的单位矩阵。表示一个高度稀疏的伴随常数矩阵。rj(j＝1,…,n+1)由系统状态矩阵经过简单运算得到，即：

所述步骤(4)中，矩阵-与向量v相乘时，将γ′改写为的形式，充分利用kronecker积的性质减少运算量。其中，表示拉格朗日系数向量。

所述步骤(6)中，需要对采用隐式arnoldi算法系统最右端的特征值时，需要计算位移-逆矩阵模值最大(递减)的部分关键特征值。(σ′n)^-1和σn具有完全相同的分块稀疏结构，而且的各个非零子块都可以被显式地表达出来，见式(19)。

γ′＝[in-r′1…-r′n](20)

式(19)中，表示维数为nn2×nn2的单位矩阵。式(21)-(23)中，s表示位移点。

所述步骤(7)中，采用【中国发明专利基于eigd的大规模时滞电力系统特征值计算方法,201510055743.3[p]】中的牛顿迭代法对步骤(6)中得到的系统近似特征值进行修正。

所述步骤(8)中，利用步骤(7)中得到的系统准确特征值λ判定系统的稳定性：若存在λ的实部大于零，即re(λ)>0，系统不稳定。若所有re(λ)<0，系统渐进稳定。若存在re(λ)＝0，则系统临界稳定。其中，re(·)表示特征值的实部。

下面给出更为详细的试验例子，如图2所示：利用山东电网系统，来验证本发明提出的基于低阶eigd方法的时滞电力系统稳定性分析方法的高效性。所有的分析均在matlab中和在intel3.4ghz8gbram台式计算机上进行。

某水平年下山东电网的规模如下：母线516个，变压器支路和输电线路936条，同步发电机114台，负荷299个。由特征值分析可知，系统存在一个频率为0.78hz、阻尼比6.97％的区间低频振荡模式，表现为山东聊城地区的发电机组相对于东部沿海地区发电机组的振荡。为了进一步提高该区间振荡模式的阻尼，考虑在聊城厂#1和#2机组上装设广域lqr。广域反馈信号分别为威海厂#3机组相对于聊城厂#1和#2机组的转速差和功角差，控制增益均为40和0.125。假设两个控制回路的综合时滞分别为90ms和100ms。系统状态变量和代数变量的维数分别为n＝1128和l＝5765。

为验证本发明提出的方法的高效性，与原有eigd方法和无时滞电力系统进行对比，表1给出了3种方法在选取不同参数时的计算效率。为计算系统的关键特征值，选取离散点数n＝10或20，位移点s＝7j和13j，在每个位移点计算的特征值个数r＝20或50。在该系统中，原有eigd方法得到的离散化矩阵约为无时滞电力系统维数的n倍。本发明提出的方法矩阵维数只比无时滞电力系统的矩阵维数略大，其迭代次数与原eigd方法大致相同，计算时间约为原eigd方法的1/16～1/8(0.08～0.12)。当系统规模增大时，其节省的计算时间更大。在同一个系统下，减少的计算时间与随n和r的增加而增加。与无时滞电力系统相比，本发明提出的方法单次计算时间略高，总计算时间仍在相同数量级上。

表1

本发明通过对系统的状态变量重新排序，剔除与时滞无关的状态变量的离散化，减少离散化矩阵的维数，从而减少dcpps的关键特征值计算的计算量。通过本发明提出的lowe-ordereigd的特征值计算方法计算dcpps的关键特征值时，尤其是针对大规模dcpps，其计算量与常规无时滞电力系统的计算量相当。与原有的基于eigd的特征值计算方法相比，其计算效率可提高10倍左右。

以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。