基于低阶SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法与流程

本发明涉及电力系统技术领域，特别是涉及基于低阶sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法，sod-ps-ii-r为“solutionoperatordiscretization-pesudospectrumiiandrotation”的英文缩写，中文含义为：解算子伪谱离散化降阶和旋转。

背景技术

随着科技的发展，尤其是电力科学与技术的发展，互联电力系统的规模逐渐增大，区间低频振荡问题越来越显著。传统的安装电力系统稳定器(powersystemstabilizer,pss)的方式，其反馈控制信号来源于当地，故不能有效抑制互联电力系统的区间振荡。

广域测量系统(wide-areameasurementsystem，wams)的出现为大规模互联电力系统稳定分析与控制的发展注入新的动力。基于wams提供的广域信息的互联电网低频振荡控制，通过引入有效反映区间振荡模式的广域反馈信号，能够获得较好的阻尼控制性能，其为阻尼互联电网中的区域间低频振荡，进而提高系统的输电能力提供了新的控制手段，具有良好并且广泛的应用前景。

广域信号在由不同通信介质(如光纤、电话线、数字微波、卫星等)组成的wams通信网络中传输和处理时，通常存在几十到几百毫秒间变化的通信延时。时滞是导致电力系统控制律失效、运行状况恶化和系统失稳的一个重要因素。因此，利用广域测量信息进行电力系统闭环控制时，必须计及时滞对电力系统的影响。

在现代电力系统中，小干扰稳定的关注点在于机电振荡问题。以状态空间模型为基础的特征值分析法是研究机电振荡的强有力工具。目前，研究人员已经提出了许多计算大规模电力系统关键特征值子集的有效方法，主要包括基于降阶系统的选择模式分析法，aesops算法和s-矩阵法，幂法、牛顿法、rayleigh商迭代法等序贯法以及同时迭代法、arnoldi算法、重新分解的双重迭代和jacobi-davidson方法等子空间迭代法。这些方法在计算部分特征值时均利用了增广状态矩阵的稀疏性，其中多数方法都是通过对原系统进行谱变换从而改变特征谱的分布，然后求取系统的特征值，最后再通过反变换得到原系统的关键特征值。但是，以上提到的方法均未考虑时滞对电力系统的影响。

中国发明专利基于padé近似的时滞电力系统特征值计算与稳定性判别方法.201210271783.8:[p].利用pade近似多项式逼近时滞环节，来计算系统最右侧的关键特征值，据此判断系统的时滞稳定性。中国发明专利基于eigd的大规模时滞电力系统特征值计算方法.201510055743.3.中国,201510055743.3[p].提出了一种基于显示igd(explicitigd,eigd)的大规模时滞电力系统特征值计算。利用计算得到的系统最右侧的关键特征值，来判断系统在固定时滞下的稳定性。

这些时滞稳定性判别方法，均需要通过多次扫描[0.1,2.5]hz低频振荡频率范围内、靠近虚轴的关键特征值，才能判断系统的时滞稳定性。

另外，现有方法的矩阵维数较大导致进行电力系统特征值计算时计算量大，计算时间长。

技术实现要素：

为了解决现有技术的不足，本发明提供了基于低阶sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法，低阶sod-ps-ii-r算法只需计算解算子离散化近似矩阵的模值最大的设定个数个特征值，通过一次计算，便能得到大规模时滞电力系统的机电振荡模式。

基于低阶sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法，包括：

建立时滞电力系统模型，根据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的对应关系，将计算时滞电力系统模型的特征值转化成计算时滞电力系统模型的解算子的特征值，从而将计算时滞电力系统机电振荡模式的问题转化为计算解算子的特征值问题；

将时滞电力系统模型的状态变量分为与时滞无关项和时滞相关项，改写时滞电力系统的数学模型，实现对其解算子维数的降阶；

通过伪谱方法对时滞电力系统模型的解算子进行离散化，得到解算子的离散化矩阵；

利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，计算解算子离散化近似矩阵模值大于1的特征值；

采用序贯法或子空间法，从模值大于1的特征值中找到解算子离散化近似矩阵的模值最大的设定个数个特征值；

将模值最大的设定个数个特征值依次经过谱映射、坐标反旋转和牛顿校验处理，最终得到时滞电力系统模型的特征值λ，特征值λ即对应时滞电力系统的机电振荡模式。

进一步优选的技术方案，所述时滞电力系统模型如下：

式中：为电力系统的状态变量向量，t为当前时刻，m为系统时滞环节个数，0<τ1<τ2<…<τi…<τm为时滞环节的时滞常数，τm为最大时滞，为系统状态矩阵，为稠密矩阵，为系统时滞状态矩阵，为稀疏矩阵，为t时刻系统状态变量导数的增量，δx(t)为t时刻系统状态变量的增量，δx(t-τi)为t-τi时刻系统状态变量的增量，δx(0)为系统状态变量的初始值即初始条件，并简写为

进一步优选的技术方案，所述时滞电力系统模型的线性化系统的特征方程为：

式中：λ为特征值，v为特征值对应的右特征向量。

进一步优选的技术方案，将时滞电力系统模型的状态变量分为与时滞无关项和时滞相关项则式(1)的第一式可改写为：

式中，n为系统状态变量个数。为t时刻系统与时滞相关状态变量导数的增量，为t时刻系统与时滞无关状态变量导数的增量，δx1(t)为t时刻系统与时滞相关状态变量的增量，δx2(t)为t时刻系统与时滞无关状态变量的增量，δx1(t-τi)为t-τi时刻系统与时滞相关状态变量的增量，δx2(t-τi)为t-τi时刻系统与时滞无关状态变量的增量；

故，式(3)可写作

式中，δx1(t)为系统与时滞相关状态变量的初始值即初始条件，可以简写为δx2(t)为系统与时滞相关状态变量的初始值即初始条件，可以简写为

进一步优选的技术方案，解算子t(h):x→x定义为将空间x上的初值条件转移到h时刻之后时滞电力系统解分段的线性算子；其中，h为转移步长，0≤h≤τm；

其中，s为积分变量，θ为变量，和分别为0和θ+h时刻时滞电力系统的状态；

由谱映射定理可知，解算子t(h)的特征值μ与时滞系统的特征值λ之间具有如下关系：

式中：σ(t(h))表示解算子的谱，\表示集合差运算。

进一步优选的技术方案，与解算子t(h)对应的、标准基形式的离散化矩阵表示如下：

式中：

式(8)中，n表示本发明提供的方法所取子区间个数，由h＝τmax/n知其决定了步长h，m为每个子区间采用第二类切比雪夫多项式的m+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数。矩阵下半块第一列为零矩阵，其余为零矩阵，为n2维单位矩阵；

式(9)中，为常数矩阵，其前m行元素完全由拉格朗日系数决定，最后一行元素均为常数；为n1维单位矩阵，为n2维单位矩阵。adj,j＝1,…,m-1，为降阶的时滞状态矩阵。式(10)中，为常数矩阵，为常数矩阵，为常数矩阵，其前m行元素完全由拉格朗日系数决定，最后一行元素均为常数。为n1维单位矩阵，为n2维单位矩阵。adi,i＝1,…,m，为降阶的时滞状态矩阵。adi降阶后的时滞状态矩。

进一步优选的技术方案，利用坐标旋转预处理方法，将坐标轴逆时针旋转θ角度，设坐标旋转后时滞电力系统模型的特征值为λ′，于是式(2)中的λ用λ′e^-jθ代替，可以得到坐标旋转后的特征方程：

式中：

λ′＝λe^jθ

τ′i＝τie^-jθ,i＝1,...,m

在坐标旋转后的特征方程中，取：

τ′i＝τie^-jθ≈τi,i＝1,...,m

从而，得到近似特征方程如下：

式中：λ″为λ′的近似值。

进一步优选的技术方案，解算子旋转后离散化矩阵的近似矩阵可表示为：

式中：

其中，和分别为经过旋转-放大预处理后的rm和bm,n。a′11、a′12、a′21、a′22分别为经过旋转-放大预处理后的a11、a12、a21、a22。

显然，当旋转角度θ＝0°时，t″m,n就退化为tm,n；

设t″m,n的特征值为μ″，其与λ″之间满足如下关系：

进一步优选的技术方案，解算子旋转后离散化矩阵的近似矩阵t″m,n的维数为(m+1)n1+(nm+1)n2，采用迭代特征值算法计算其模值最大的设定个数个特征值；

假设在第k次迭代时，需要计算t″m,n与向量的乘积，具体步骤如下：

步骤(1)：由于t″m,n具有的逻辑结构，可知v^k的第(m+1)n1+n2+1～(m+1)n1+1+((n'-1)m+1)n2个分量等于z^k的(m+1)n+1～(m+1)n1+(n'm+1)n2第个分量；

步骤(2)：将向量v^k分为和并按列压缩为矩阵和相应地，有：其中vec(·)为将矩阵压缩为列向量的运算；

步骤(3)：计算r为中间向量；

步骤(4)：计算

进一步优选的技术方案，所述步骤(3)的步骤如下：

将b″m,n表达式代入，可得：

由上述分析可知，要计算r^k，首先要计算v1^kp1^t并将其压缩为一个(m+1)n1维的列向量，然后计算对其进行求和，然后将其压缩为一个(m+1)n2维的列向量并乘以系数e^jθ，再计算v1^kp1^t并将其压缩为一个(m+1)n2维的列向量，然后将两个列向量相加得到一个一维的列向量，与v1^kp1^t压缩形成的(m+1)n1维的列向量一起组成r^k。

进一步优选的技术方案，所述步骤(4)的步骤如下：

采用迭代算法来计算求解过程中，涉及矩阵向量乘法运算b＝r″my，其中

首先，将向量y先分为和并分别按列压缩为矩阵和进而，得到b＝r″my的稀疏实现步骤如下：

进一步优选的技术方案，在计算得到解算子离散化近似矩阵的模值最大的设定个数个特征值μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，λ的计算公式如下：

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

第一、本发明提出的低阶sod-ps-ii-r算法充分考虑了电力系统的实际规模，以及通信延时对系统的影响，在用于计算实际系统机电振荡模式对应的关键特征值时具有很好的效果。

第二、本发明提出的低阶sod-ps-ii-r算法因施加旋转，较低阶sod-ps-ii算法的收敛速度更快。

第三、本发明的核心创新点在于考虑了时滞电力系统状态变量与时滞的相关性，使时滞电力系统模型的维数大大降低，据此得到低阶的时滞系统解算子离散化矩阵，从而提高计算解算子离散化矩阵特征值的效率。

附图说明

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

图1为时滞电力系统示意图。

图2为基于低阶sod-ps-ii-r时滞电力系统机电振荡模式计算方法的流程图。

具体实施方式

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

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

如图1所示：在实际大规模电力系统的建模过程中引入通信时滞。时滞电力系统包含无时滞电力系统、广域反馈时滞、广域阻尼控制器和广域输出时滞四部分，各部分之间的连接关系如图所示。图1中，无时滞电力系统的输出为yf，ydf为考虑时滞后的广域反馈信号并作为阻尼控制器的输入，yc为广域阻尼控制器的输出，ydc为无时滞电力系统的控制输入。

将时滞电力系统模型在其稳态运行点附近进行线性化，即可得到用于时滞电力系统小干扰稳定分析的线性化模型。

根据谱映射定理，所述时滞电力系统位于虚轴附近、频率在[0,2]hz范围内的特征值λ，被映射为解算子t(h)的特征值，并密集地分布在单位圆的附近。

时滞电力系统的特征值对应的运行模式主要可分为三类，机电振荡模式，缓慢模式，迅速模式。其中机电振荡模式对应的特征值所在区域，是本算法求解的对象。

时滞系统位于左半复平面的特征值被映射为解算子模值小于1的特征值，而时滞系统位于右半复平面的特征值被映射为解算子模值大于1的特征值。因此，利用解算子的特征值，就可以判断原时滞系统的稳定性。如果解算子至少存在一个模值大于1的特征值，则可以判断原时滞系统是不稳定的，如果解算子所有特征值的模值均小于1，则原时滞系统是渐进稳定的。

本申请的整体技术构思为：将所述时滞电力系统的解算子t(h)进行伪谱离散化，通过计算有限维离散化矩阵tm,n模值最大的设定个数个特征值，便可将其转化为时滞电力系统机电振荡模式对应的关键特征值。

为了提高所述时滞电力系统特征值算法的收敛性，需要增加tm,n的特征值之间的相对距离。利用坐标旋转预处理方法可将时滞电力系统的阻尼比小于给定常数ζ的关键特征值λ变换为对应的tm,n模值大于1的特征值。

由于所述时滞电力系统解算子t(h)的离散化点数必须为整数，要求坐标旋转后的最大时滞为实数。当旋转角度不为0时，得到旋转后的近似特征方程和对应的解算子离散化近似矩阵t″m,n。

对于所述时滞电力系统离散化近似矩阵t″m,n的维数为(m+1)n1+(n'm+1)n2，当系统规模较大时，需要采用稀疏特征值算法(如隐式重启动arnoldi算法)，来计算其模值最大的设定个数个特征值μ″。在计算得到μ″之后，首先通过谱映射关系得到与之对应的经过坐标旋转后时滞电力系统的特征值λ″。然后，再通过坐标反旋转得到时滞电力系统的特征值λ″e^-jθ。

本申请的一种典型的实施方式中，如图2所示，基于低阶sod-ps-ii-r的时滞电力系统机电振荡模式对应的特征值计算方法，包括如下步骤：

s1：建立时滞电力系统模型；

s2：将时滞电力系统状态变量分为与时滞无关项和时滞相关项，改写时滞电力系统的数学模型，实现对其解算子维数的降阶；

s3：通过伪谱离散化，得到解算子t(h)的有限维离散化矩阵tm,n；

s4：利用坐标旋转预处理方法，将时滞电力系统的阻尼比小于给定常数ζ的关键特征值λ，变换为tm,n旋转后的近似矩阵t″m,n模值大于1的部分特征值μ″；

s5：采用序贯法或子空间法(如隐式重启动arnoldi算法)来计算步骤s3得到的解算子离散化近似矩阵t″m,n的模值最大的设定个数个特征值μ″；

s6：在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验得到时滞电力系统的特征值λ。

至此，时滞电力系统的机电振荡模式对应的设定个数个关键特征值已被计算得到。

所述步骤s1中，考虑广域通信时滞影响后，可用如下的一组时滞微分方程组来描述电力系统为：

式中：为电力系统的状态变量向量，t为当前时刻，m为系统时滞环节个数。0<τ1<τ2<…<τi…<τm为时滞环节的时滞常数，τm为最大时滞。为系统状态矩阵，为稠密矩阵。为系统时滞状态矩阵，为稀疏矩阵。为t时刻系统状态变量导数的增量，δx(t)为t时刻系统状态变量的增量，δx(t-τi)为t-τi时刻系统状态变量的增量。δx(0)为系统状态变量的初始值(即初始条件)，可以简写为

上式表示的线性化系统的特征方程为：

式中：λ和v分别为特征值和相应的右特征向量。

本申请中，根据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的关系，将计算时滞电力系统模型的特征值转化成计算解算子的特征值；从而将计算时滞电力系统机电振荡模式的问题转化为计算解算子的特征值问题。

所述步骤s2中，状态变量分为与时滞无关项和时滞相关项则时滞电力系统的微分方程可改写为：

式中，n为系统状态变量个数。为t时刻系统与时滞相关状态变量导数的增量，为t时刻系统与时滞无关状态变量导数的增量。δx1(t)为t时刻系统与时滞相关状态变量的增量，δx2(t)为t时刻系统与时滞无关状态变量的增量。δx1(t-τi)为t-τi时刻系统与时滞相关状态变量的增量，δx2(t-τi)为t-τi时刻系统与时滞无关状态变量的增量。

故，描述时滞电力系统的微分方程可写作

式中，δx1(0)为系统与时滞相关状态变量的初始值(即初始条件)，并简写为δx2(0)为系统与时滞相关状态变量的初始值(即初始条件)，并简写为

步骤s3中，解算子t(h):x→x定义为将空间x上的初值条件转移到h(转移步长，0≤h≤τm)时刻之后系统解分段的线性算子。

由谱映射定理可知，解算子t(h)的特征值μ与时滞系统的特征值λ之间具有如下关系。

式中：σ(t(h))表示解算子的谱，\表示集合差运算。

时滞电力系统位于复平面左半侧的特征值被映射到解算子单位圆之内，而时滞电力系统位于复平面右半侧的特征值被映射到单位圆之外。因此，利用解算子的特征值，即可判断原时滞电力系统的稳定性。若解算子有模值大于1的特征值，则原时滞电力系统是不稳定的，若解算子所有特征值的模值均小于1，则原时滞电力系统是渐进稳定的。

解算子t(h)是描述x→x映射的无穷维线性算子。为了计算解算子的特征值并依此得到原时滞系统机电振荡模式对应的特征值，首先采用伪谱方法(pesudospectral,ps)对t(h)进行离散化，得到一个与解算子对应的、有限维的近似矩阵，进而计算近似矩阵的特征值并得到原时滞系统机电振荡模式对应的特征值。与解算子t(h)对应的、标准基形式的离散化矩阵可表示如下：

式中：

式中，n表示本发明提供的方法所取子区间个数，由h＝τmax/n知其决定了步长h，m为每个子区间采用第二类切比雪夫多项式的m+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数。矩阵下半块第一列为零矩阵，其余为零矩阵。为n2维单位矩阵。

式中，为常数矩阵。其前m行元素完全由拉格朗日系数决定，最后一行元素为常数。为n1维单位矩阵，为n2维单位矩阵。adj,j＝1,…,m-1，为降阶的时滞状态矩阵。

式中，为常数矩阵，为常数矩阵，为常数矩阵。其前m行元素完全由拉格朗日系数决定，最后一行元素均为常数。为n1维单位矩阵，为n2维单位矩阵。adi,i＝1,…,m，为降阶的时滞状态矩阵。

步骤s4中，利用坐标旋转预处理方法，将坐标轴逆时针旋转θ角度，将时滞电力系统的阻尼比小于给定常数ζ(＝sinθ)的关键特征值λ，变换为tm,n模值大于1的部分特征值μ。旋转后的虚轴对应着原坐标系中阻尼比为ζ的等阻尼比线。

设坐标旋转后时滞电力系统的特征值为λ′，于是将谱映射关系中的λ用λ′e^-jθ代替，便得到坐标旋转后的特征方程：

式中：

λ′＝λe^jθ

τ′i＝τie^-jθ,i＝1,...,m

由于n为整数，由关系h＝τmax/n可知，坐标旋转后的最大时滞τ′max须为实数。然而，当θ≠0时，由τ′i表达式可知τ′i(i＝1,…,m)必为复数。显然，这与τ′max为实数是相矛盾的。因此，在坐标旋转后的特征方程中，取：

τ′i＝τie^-jθ≈τi,i＝1,...,m

从而，得到近似特征方程如下：

式中：λ″为λ′的近似值。

上式对应的解算子离散化矩阵可表示为：

式中：

显然，当旋转角度θ＝0°时，t″m,n就退化为tm,n。

设t″m,n的特征值为μ″，其与λ″之间满足如下关系：

步骤s5中，矩阵t″m,n的维数为(m+1)n1+(n'm+1)n2。对于大规模电力系统，矩阵t″m,n的维数将十分巨大。因此，在应用解算子对应的离散化矩阵t″m,n求解大规模电力系统的时滞特征值时，必须采用迭代特征值算法(序贯法或子空间法)计算其模值最大的设定个数个特征值。

假设在第k次迭代时，需要计算t″m,n与向量的乘积，具体步骤如下：

(1)：由于t′m具有的特殊逻辑结构，可知v^k的第(m+1)n1+n2+1～(m+1)n1+1+((n'-1)m+1)n2个分量等于z^k的(m+1)n+1～(m+1)n1+(n'm+1)n2第个分量。

(2)：将向量v^k分为和并按列压缩为矩阵和相应地，有：其中vec(·)为将矩阵压缩为列向量的运算。

(3)：计算r为中间向量。

(4)：计算

步骤(3)的步骤如下：

将b″m,n表达式代入，可得：

由上述分析可知，要计算r^k，首先要计算v1^kp1^t并将其压缩为一个(m+1)n1维的列向量，然后计算对其进行求和，然后将其压缩为一个(m+1)n2维的列向量并乘以系数e^jθ，再计算v1^kp1^t并将其压缩为一个(m+1)n2维的列向量，然后将两个列向量相加得到一个一维的列向量，与v1^kp1^t压缩形成的(m+1)n1维的列向量一起组成r^k。

步骤(4)的步骤如下：

易知没有显示表达。因而，这里采用迭代算法来计算求解过程中，涉及矩阵向量乘法运算b＝r″my，其中

首先，将向量y先分为和并分别按列压缩为矩阵和进而，可以得到b＝r″my的稀疏实现步骤如下：

步骤s6中，在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，λ的计算公式如下：

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