CN108808706A - 基于sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 - Google Patents
基于sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 Download PDFInfo
- Publication number
- CN108808706A CN108808706A CN201810776545.XA CN201810776545A CN108808706A CN 108808706 A CN108808706 A CN 108808706A CN 201810776545 A CN201810776545 A CN 201810776545A CN 108808706 A CN108808706 A CN 108808706A
- Authority
- CN
- China
- Prior art keywords
- time
- power system
- matrix
- lag power
- lag
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J3/00—Circuit arrangements for ac mains or ac distribution networks
- H02J3/24—Arrangements for preventing or reducing oscillations of power in networks
-
- H—ELECTRICITY
- H02—GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
- H02J—CIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
- H02J2203/00—Indexing scheme relating to details of circuit arrangements for AC mains or AC distribution networks
- H02J2203/20—Simulating, e g planning, reliability check, modelling or computer assisted design [CAD]
Landscapes
- Engineering & Computer Science (AREA)
- Power Engineering (AREA)
- Supply And Distribution Of Alternating Current (AREA)
Abstract
本发明公开了基于SOD‑PS‑II‑R算法的时滞电力系统机电振荡模式计算方法，包括：建立时滞电力系统模型；通过伪谱方法对电力系统模型的解算子进行离散化，得到时滞电力系统模型的解算子的离散化矩阵；利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，计算旋转后的解算子离散化近似矩阵模值大于1的特征值；采用序贯法或子空间法，从上述模值大于1的特征值中找到解算子离散化近似矩阵的模值最大的设定数量的特征值；对模值最大的设定数量的特征值，依次经过谱映射、坐标反旋转和牛顿校验处理，最终得到时滞电力系统模型的特征值，其对应时滞电力系统的机电振荡模式。
Description
技术领域
本发明涉及电力系统技术领域，特别是涉及基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，SOD-PS-II-R为“Solution Operator Discretization-PesudoSpectrum II and Rotation”的英文缩写，中文含义为：解算子伪谱离散化和旋转。
背景技术
电力科学与技术的发展使得互联电力系统的规模逐渐增大，区间低频振荡问题愈发突出。传统采用安装电力系统稳定器(Power System Stabilizer,PSS)的方案来解决低频振荡问题，但是其反馈控制信号来源于通常较远，因此不能有效抑制互联电力系统的区间振荡。
广域测量系统(Wide-Area Measurement System，WAMS)的出现带给大规模互联电力系统稳定分析与控制新的发展机遇。基于WAMS提供的广域信息的互联电网低频振荡控制，通过引入有效反映互联电力系统区间振荡模式的广域反馈信号，能获得较好的阻尼控制性能，其为抑制互联电网中的区域间低频振荡，从而提高系统的输电能力提供了更好的控制手段，具有良好且广泛的应用前景。
广域信号在由不同通信介质(如光纤、电话线、数字微波、卫星等)组成的WAMS通信网络中传输和处理时，存在几十到几百毫秒间变化的通信延时。通信延时是引起系统控制律失效、运行状况恶化和系统失稳的一种重要原因。因此，在利用广域测量信息进行电力系统进行闭环控制时，必须计及通信延时的影响。
机电振荡问题是现代电力系统小干扰稳定的关注点。以状态空间模型为基础的特征值分析法是研究机电振荡的强有力工具。目前，研究人员已经提出了许多计算大规模电力系统关键特征值子集的有效方法，主要包括基于降阶系统的选择模式分析法，AESOPS算法和S-矩阵法，幂法、牛顿法、Rayleigh商迭代法等序贯法以及同时迭代法、Arnoldi算法、重新分解的双重迭代和Jacobi-Davidson方法等子空间迭代法。这些方法在计算电力系统部分特征值时都利用了增广状态矩阵的稀疏性，其中多数方法是通过对原系统进行谱变换从而改变特征谱的分布，然后求取系统的特征值，最后通过反变换得到原系统的关键特征值。但是，以上提到的方法均未考虑通信延时对电力系统的影响。
中国发明专利基于Padé近似的时滞电力系统特征值计算与稳定性判别方法.201210271783.8:[P].利用Pade近似多项式逼近时滞环节计算系统最右侧的关键特征值，从而判断系统的时滞稳定性。中国发明专利基于EIGD的大规模时滞电力系统特征值计算方法.201510055743.3.中国,201510055743.3[P].提出了一种基于显示IGD(Explicit IGD,EIGD)的大规模时滞电力系统特征值计算。根据计算得到的电力系统最右侧的关键特征值，可以判断系统在固定时滞下的稳定性。
上述这些时滞稳定性判别方法，均需要通过多次扫描[0.1,2.5]Hz低频振荡频率范围内、靠近虚轴的关键特征值，才能判断系统的时滞稳定性，操作步骤繁琐，计算量大。
发明内容
为了解决现有技术的不足，本发明提供了基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，用以得到大规模时滞电力系统的机电振荡模式。本发明SOD-PS-II-R算法只需计算一次便可得到解算子离散化近似矩阵的模值最大的设定数量个特征值，也即大规模时滞电力系统的机电振荡模式。
基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，包括：
建立时滞电力系统模型，依据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的对应关系，将计算时滞电力系统模型的特征值转化成计算解算子的特征值，从而将计算时滞电力系统机电振荡模式的问题转化为计算解算子的特征值问题；
通过伪谱方法对电力系统模型的解算子进行离散化，得到时滞电力系统模型的解算子的离散化矩阵；
利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，计算旋转后的解算子离散化近似矩阵模值大于1的特征值；
采用序贯法或子空间法，从上述模值大于1的特征值中找到解算子离散化近似矩阵的模值最大的设定数量的特征值；
对模值最大的设定数量的特征值，依次经过谱映射、坐标反旋转和牛顿校验处理，最终得到时滞电力系统模型的特征值，其对应时滞电力系统的机电振荡模式。
进一步优选的技术方案，所述时滞电力系统模型如下：
式中：为电力系统的状态变量向量，m为系统时滞环节的个数，t为当前时刻，0<τ1<τ2<…<τi…<τm为时滞环节的时滞常数，为t时刻系统状态变量导数的增量，Δx(t)为t时刻系统状态变量的增量，Δx(t-τi)为t-τi时刻系统状态变量的增量，Δx(0)为系统状态变量的初始值即初始条件，并简写为 为系统状态矩阵，为稠密矩阵，为系统时滞状态矩阵，为稀疏矩阵，i＝1,…,m。
进一步优选的技术方案，所述时滞电力系统模型的线性化系统的特征方程为：
式中：λ为特征值，v为特征值对应的右特征向量。
进一步优选的技术方案，解算子T(h):X→X定义为将空间X上的初值条件转移到h时刻之后时滞电力系统解分段的线性算子；其中，h为转移步长，0≤h≤τm；
其中，s为积分变量，θ为变量，和分别为0和θ+h时刻时滞电力系统的状态；
由谱映射定理可知，解算子T(h)的特征值μ与时滞电力系统的特征值λ之间具有如下关系：
式中：σ(T(h))表示解算子的谱，\表示集合差运算。
进一步优选的技术方案，与解算子T(h)对应的、标准基形式的离散化矩阵表示如下：
式中：
式(5)中，N表示本发明提供的方法所取子区间个数，由h＝τmax/N知其决定了步长h，M为每个子区间采用第二类切比雪夫多项式的M+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数。矩阵下半块最后一列为零矩阵，其余为零矩阵，In为n维单位矩阵；
式(6)中，为常数矩阵；其前M行元素完全由拉格朗日系数决定，最后一行元素为常数；In为n维单位矩阵；A0为系统状态矩阵；Aj，j＝1,…,m-1为系统时滞状态矩阵；
式(7)中，为常数矩阵，为常数矩阵。其前M行元素完全由拉格朗日系数决定，最后一行元素为常数；In为n维单位矩阵。
步长h＝τmax/N，N为所取子区间个数，M为采用第二类切比雪夫多项式的M+1个经过位移和归一化处理的极点对每个步长为h的子区间进行离散化的参数，Aj和Ai同为系统时滞状态矩阵，只是下标不同。LAj其前M行元素完全由拉格朗日系数决定，最后一行元素为常数。
进一步优选的技术方案，利用坐标旋转预处理方法，将坐标轴逆时针旋转θ角度，将时滞电力系统模型的阻尼比小于设定常数ζ的关键特征值λ，变换为TM,N模值大于1的特征值μ，旋转后的虚轴对应原坐标系中阻尼比为ζ的等阻尼比线，ζ＝sinθ；
所求电力系统特征值中小于给定阻尼比的部分特征值定义为关键特征值。
设坐标旋转后时滞电力系统模型的特征值为λ′，于是式(2)中的λ用λ′e-jθ代替，可以得到坐标旋转后的特征方程：
式中：
λ′＝λejθ
τ′i＝τie-jθ,i＝1,...,m
由于N为整数，由关系h＝τmax/N可知，坐标旋转后的最大时滞τ′max须为实数。然而，当θ≠0时，由τ′i表达式可知τ′i必为复数。因此，在坐标旋转后的特征方程中，取：
τ′i＝τie-jθ≈τi,i＝1,...,m
从而，得到近似特征方程如下：
式中：λ″为λ′的近似值。
进一步优选的技术方案，式(9)对应的解算子旋转后离散化矩阵的近似矩阵N′＝N/α，α为放大倍数，可表示为：
式中：
显然，当旋转角度θ＝0°时，T″M,N即为TM,N；
设T″M,N的特征值为μ″，其与λ″之间满足如下关系：
进一步优选的技术方案，解算子离散化近似矩阵T″M,N的维数为(N′M+1)n，采用迭代特征值算法计算其模值最大的设定个数个特征值；
具体计算过程如下：
假设在第k次迭代时，需要计算T″M,N与向量的乘积， 其中，n为系统状态变量个数，N为需要设定的子区间个数，N′＝N/α，α为放大倍数，M为每个子区间采用第二类切比雪夫多项式的M+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数。具体步骤如下：
步骤(1)：由于T″M,N具有的逻辑结构，可知vk的第n+1～((N′-1)M+1)n个分量等于的(M+1)n+1～(N′M+1)n第个分量；
步骤(2)：将向量vk按列压缩为一矩阵 相应地，有：vk＝vec(Vk)，其中vec(·)为将矩阵压缩为列向量的运算；
步骤(3)：计算rk为中间向量；r为中间向量，B″M,N为施加旋转-放大预处理技术后的BM,N；
步骤(4)：计算R″M为施加旋转-放大预处理技术后的RM；
进一步优选的技术方案，所述步骤(3)的步骤如下：
将B″M,N表达式代入，可得：
由上述分析可知，要计算rk，首先要计算然后再进行求和，然后将其压缩为一个(M+1)n维的列向量并乘以系数ejθ，再计算VkPT并将其压缩为一个(M+1)n维的列向量，最后将两个列向量相加。
进一步优选的技术方案，所述步骤(4)的步骤如下：
没有显示表达，因而，采用迭代算法来计算求解过程中，涉及矩阵向量乘法运算b＝R″My，其中y为列向量；
首先，将向量y按列压缩为一矩阵进而，得到b＝R″My的稀疏实现步骤如下：
可以稀疏实现。
进一步优选的技术方案，在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，λ的计算公式如下：
与现有技术相比，本发明的有益效果是：
第一、本发明提出的SOD-PS-II-R算法具有良好的适用性，能够计算实际大规模时滞电力系统机电振荡模式对应的特征值。
第二、本发明提出的SOD-PS-II-R算法在计算时滞电力系统机电振荡模式对应的特征值时，只需计算一次解算子离散化近似矩阵的模值最大的设定个数特征值。
第三，本发明提出的SOD-PS-II-R算法能够高效计算相差20倍以上的两时滞电力系统机电振荡模式特征值，补充了已有算法的不足。
第四、本发明的方法是基于SOD-PS-II-R的时滞电力系统机电振荡模式对应的特征值计算方法，核心创新点在用伪谱离散化将无限维的时滞电力系统转换为有限维的离散化矩阵，并对其加以旋转，然后利用矩阵向量相乘的稀疏实现以及子空间迭代算法得到时滞系统离散化矩阵的特征值，最后通过谱映射关系将离散化近似矩阵的特征值转换为时滞电力系统的特征值。
附图说明
构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。
图1为基于SOD-PS-II-R时滞电力系统机电振荡模式计算方法的流程图。
具体实施方式
应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
在实际大规模电力系统的建模过程中引入通信时滞。时滞电力系统包含无时滞电力系统、反馈时滞、阻尼控制器和输出时滞四部分。设无时滞电力系统的输出为yf，ydf为考虑时滞后的广域反馈信号并作为阻尼控制器的输入，yc为广域阻尼控制器的输出，ydc为无时滞电力系统的控制输入。
时滞电力系统的特征值对应的运行模式主要可分为三类，机电振荡模式，缓慢模式，迅速模式。其中机电振荡模式对应的特征值所在区域，是本算法求解的对象。
时滞电力系统位于复平面虚轴左侧的特征值被映射到解算子单位圆之内，即解算子特征值的模值小于1，而时滞电力系统位于复平面虚轴右侧的特征值被映射到解算子单位圆之外，即解算子特征值的模值大于1。因此，根据解算子的特征值，即可判断原时滞电力系统的稳定性。只要解算子存在模值大于1的特征值，便可判断原时滞电力系统是不稳定的，若解算子所有特征值的模值均小于1，则原时滞系统是渐进稳定的。
利用坐标旋转预处理方法，将坐标轴逆时针旋转θ角度，将时滞电力系统的阻尼比小于给定常数ζ(＝sinθ)的特征值λ，变换为TM,N模值大于1的部分特征值μ。旋转后的虚轴对应着原坐标系中阻尼比为ζ的虚线。
本申请的一种典型的实施方式中，如图1所示，基于SOD-PS-II-R的时滞电力系统机电振荡模式对应的特征值计算方法，包括如下步骤：
S1：建立时滞电力系统模型；
S2：通过伪谱离散化，得到解算子T(h)的有限维离散化矩阵TM,N；
S3：利用坐标旋转预处理方法，将时滞电力系统的阻尼比小于给定常数ζ的关键特征值λ，变换为TM,N旋转后的近似矩阵T″M,N模值大于1的部分特征值μ″；
S4：采用序贯法或子空间法(如隐式重启动Arnoldi算法)来计算步骤S3得到的解算子离散化近似矩阵T″M,N的模值最大的设定个数个特征值μ″；
S5：在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验得到时滞电力系统的特征值λ。
至此，已计算得到时滞电力系统的机电振荡模式对应的设定个数个关键特征值。
所述步骤S1中，考虑广域通信时滞影响后，电力系统可用如下的一组时滞微分方程组来描述：
式中：为电力系统的状态变量向量，m为系统时滞环节的个数，t为当前时刻。0<τ1<τ2<…<τi…<τm为时滞环节的时滞常数。为t时刻系统状态变量导数的增量，Δx(t)为t时刻系统状态变量的增量，Δx(t-τi)为t-τi时刻系统状态变量的增量。Δx(0)为系统状态变量的初始值(即初始条件)，并简写为 为系统状态矩阵，为稠密矩阵。 为系统时滞状态矩阵，为稀疏矩阵。
上式表示的线性化系统的特征方程为：
式中：λ和v分别为特征值和相应的右特征向量。
本申请依据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的关系，将计算时滞电力系统模型的特征值转化成计算解算子的特征值；从而将计算时滞电力系统机电振荡模式的问题转化为计算解算子的特征值问题。
所述步骤S2中，解算子T(h):X→X定义为将空间X上的初值条件转移到h(转移步长，0≤h≤τm)时刻之后系统解分段的线性算子。
由谱映射定理可知，解算子T(h)的特征值μ与时滞系统的特征值λ之间具有如下关系。
式中：σ(T(h))表示解算子的谱，\表示集合差运算。
时滞电力系统位于复平面虚轴左侧的特征值被映射到解算子单位圆之内，即解算子特征值的模值小于1，而时滞电力系统位于复平面虚轴右侧的特征值被映射到解算子单位圆之外，即解算子特征值的模值大于1。因此，根据解算子的特征值，即可判断原时滞电力系统的稳定性。只要解算子存在模值大于1的特征值，便可判断原时滞电力系统是不稳定的，若解算子所有特征值的模值均小于1，则原时滞系统是渐进稳定的。
解算子T(h)是描述X→X映射的无穷维线性算子。为了计算解算子的特征值并依此得到原时滞电力系统机电振荡模式对应的特征值，需要首先采用伪谱方法(Pesudospectral,PS)对解算子T(h)进行离散化，得到一个与解算子对应的、有限维的近似矩阵，进而计算近似矩阵的特征值并得到原时滞电力系统机电振荡模式对应的特征值。与解算子T(h)对应的、标准基形式的离散化矩阵TM,N可表示如下：
式中：
式中，N表示本发明提供的方法所取子区间个数，由h＝τmax/N知其决定了步长h，M为每个子区间采用第二类切比雪夫多项式的M+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数。TM,N矩阵下半块最后一列为零矩阵，其余为零矩阵。In为n维单位矩阵。
式中，为常数矩阵。其前M行元素完全由拉格朗日系数决定，最后一行元素为常数。In为n维单位矩阵；A0为系统状态矩阵；Aj，j＝1,…,m-1为系统时滞状态矩阵；
式中，为常数矩阵，为常数矩阵。其前M行元素完全由拉格朗日系数决定，最后一行元素为常数。
所述步骤S3中，利用坐标旋转预处理方法，将坐标轴逆时针旋转θ角度，将时滞电力系统的阻尼比小于给定常数ζ(＝sinθ)的关键特征值λ，变换为TM,N模值大于1的部分特征值μ。旋转后的虚轴对应着原坐标系中阻尼比为ζ的等阻尼比线。
设坐标旋转后时滞电力系统的特征值为λ′，于是将谱映射关系中的λ用λ′e-jθ代替，可以得到坐标旋转后的特征方程：
式中：
λ′＝λejθ
τ′i＝τie-jθ,i＝1,...,m
由于N为整数，由关系h＝τmax/N可知，坐标旋转后的最大时滞τ′max须为实数。然而，当θ≠0时，由τ′i表达式可知τ′i(i＝1,…,m)必为复数。因此，为解决此问题，在坐标旋转后的特征方程中，取：
τ′i＝τie-jθ≈τi,i＝1,...,m
从而，得到近似特征方程如下：
式中：λ″为λ′的近似值。
上式对应的解算子离散化矩阵T″M,N可表示为：
式中：
显然，当旋转角度θ＝0°时，T″M,N就退化为TM,N。
设T″M,N的特征值为μ″，其与λ″之间满足如下关系：
所述步骤S4中，矩阵T″M,N的维数为(N′M+1)n。对于大规模电力系统，矩阵T″M,N的维数将十分巨大。因此，在应用解算子对应的离散化矩阵T″M,N求解大规模电力系统的时滞特征值时，必须采用迭代特征值算法(序贯法或子空间法)计算其模值最大的设定个数个特征值。
假设在第k次迭代时，需要计算T″M,N与向量的乘积， 其中，n为系统状态变量个数，N为需要设定的子区间个数，N′＝N/α，α为放大倍数，M为每个子区间采用第二类切比雪夫多项式的M+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数。具体步骤如下：
(1)：由T″M,N具有的逻辑结构，可知vk的第n+1～((N′-1)M+1)n个分量等于的第(M+1)n+1～(N′M+1)n个分量。
(2)：将向量vk按列压缩为一矩阵 相应地，有：vk＝vec(Vk)，其中vec(·)为将矩阵压缩为列向量的运算。
(3)：计算rk为中间向量。
(4)：计算
步骤(3)的步骤如下：
将B″M,N表达式代入，可得：
由上述分析可知，要计算wk，首先要计算然后再进行求和，然后将其压缩为一个(M+1)n维的列向量并乘以系数ejθ，再计算VkPT并将其压缩为一个(M+1)n维的列向量，最后将两个列向量相加。
值得注意的是，可以稀疏实现，从而减小计算量，提高计算效率。
所述步骤(4)的步骤如下：
易知没有显示表达。因而，这里采用迭代算法来计算求解过程中，涉及矩阵向量乘法运算b＝R″My，其中
首先，将向量y按列压缩为一矩阵进而，可以得到b＝R″My的稀疏实现步骤如下：
值得注意的是，可以稀疏实现，从而减小计算量，提高计算效率。
所述步骤S5中，在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，λ的计算公式如下：
为了验证本发明提出的基于SOD-PS-II方法的时滞电力系统机电振荡模式计算方法的效果，以山东电网系统为基础，所有分析均在Matlab中和在Intel 3.4GHz 8GB RAM台式计算机上进行。
某水平年下山东电网的参数为：同步发电机114台，母线516个，变压器支路和输电线路936条，负荷299个。在聊城厂的两个机组上装设广域LQR，广域反馈信号分别为威海厂的一个机组对于聊城厂两个机组的转速差和功角差，控制增益均为40和0.125。系统状态变量的维数为n＝1128，系统代数变量的维数为l＝5765。
为验证本发明提出的方法的效果，与SOD-PS-R方法进行对比，表1给出了两种方法在选取不同时滞参数τ1和τ2时的计算时间。为得到系统的机电振荡模式的关键特征值，选取本发明提出的方法的离散点参数为M＝3，N＝3，SOD-PS-R方法的离散点参数为M＝3，N＝4，保证两种方法的矩阵维数相同，两种方法计算的特征值个数均为r＝100。随着两时滞相差倍数的增大，本发明提出的方法的计算速度相对SOD-PS-R算法的计算速度越快。
表1
以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。
Claims (10)
1.基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，包括：
建立时滞电力系统模型，依据时滞电力系统模型的特征值与时滞电力系统模型的解算子特征值之间的对应关系，将计算时滞电力系统模型的特征值转化成计算解算子的特征值，从而将计算时滞电力系统机电振荡模式的问题转化为计算解算子的特征值问题；
通过伪谱方法对电力系统模型的解算子进行离散化，得到时滞电力系统模型的解算子的离散化矩阵；
利用坐标旋转预处理方法，将解算子的离散化矩阵旋转并进行近似得到旋转后的解算子离散化近似矩阵，计算旋转后的解算子离散化近似矩阵模值大于1的特征值；
采用序贯法或子空间法，从上述模值大于1的特征值中找到解算子离散化近似矩阵的模值最大的设定数量的特征值；
对模值最大的设定数量的特征值，依次经过谱映射、坐标反旋转和牛顿校验处理，最终得到时滞电力系统模型的特征值，其对应时滞电力系统的机电振荡模式。
2.如权利要求1所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，所述时滞电力系统模型如下：
式中：为电力系统的状态变量向量，m为系统时滞环节的个数，t为当前时刻，0<τ1<τ2<…<τi…<τm为时滞环节的时滞常数，为t时刻系统状态变量导数的增量，Δx(t)为t时刻系统状态变量的增量，Δx(t-τi)为t-τi时刻系统状态变量的增量，Δx(0)为系统状态变量的初始值即初始条件，并简写为 为系统状态矩阵，为稠密矩阵，为系统时滞状态矩阵，为稀疏矩阵，i＝1,…,m；
所述时滞电力系统模型的线性化系统的特征方程为：
式中：λ为特征值，v为特征值对应的右特征向量。
3.如权利要求1所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，解算子T(h):X→X定义为将空间X上的初值条件转移到h时刻之后时滞电力系统解分段的线性算子；其中，h为转移步长，0≤h≤τm；
其中，s为积分变量，θ为变量，和分别为0和θ+h时刻时滞电力系统的状态；
由谱映射定理可知，解算子T(h)的特征值μ与时滞电力系统的特征值λ之间具有如下关系：
式中：σ(T(h))表示解算子的谱，\表示集合差运算。
4.如权利要求1所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，与解算子T(h)对应的、标准基形式的离散化矩阵TM,N表示如下：
式中：
式(5)中，N表示所取子区间个数，由h＝τmax/N知其决定了步长h，M为每个子区间采用第二类切比雪夫多项式的M+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数，矩阵下半块最后一列为零矩阵，其余为零矩阵，In为n维单位矩阵；
式(6)中，为常数矩阵，其前M行元素完全由拉格朗日系数决定，最后一行元素为常数；In为n维单位矩阵；A0为系统状态矩阵；Aj，j＝1,…,m-1为系统时滞状态矩阵；
式(7)中，为常数矩阵，为常数矩阵；其前M行元素完全由拉格朗日系数决定，最后一行元素为常数。
5.如权利要求1所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，利用坐标旋转预处理方法，将坐标轴逆时针旋转θ角度，将时滞电力系统模型的阻尼比小于设定常数ζ的关键特征值λ，变换为TM,N模值大于1的特征值μ，旋转后的虚轴对应原坐标系中阻尼比为ζ的等阻尼比线，ζ＝sinθ；
设坐标旋转后时滞电力系统模型的特征值为λ′，于是式(2)中的λ用λ′e-jθ代替，可以得到坐标旋转后的特征方程：
式中：
λ′＝λejθ
τi′＝τie-jθ,i＝1,...,m
由于N为整数，由关系h＝τmax/N可知，坐标旋转后的最大时滞τ′max须为实数，然而，当θ≠0时，由τi′表达式可知τi′必为复数，因此，在坐标旋转后的特征方程中，取：
τi′＝τie-jθ≈τi,i＝1,...,m
从而，得到近似特征方程如下：
式中：λ″为λ′的近似值。
6.如权利要求5所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，式(9)对应的解算子旋转后离散化矩阵的近似矩阵T″M,N可表示为：
式中：
显然，当旋转角度θ＝0°时，T″M,N即为TM,N；
设T″M,N的特征值为μ″，其与λ″之间满足如下关系：
7.如权利要求1所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，解算子离散化近似矩阵T″M,N的维数为(N′M+1)n，采用迭代特征值算法计算其模值最大的设定个数个特征值；
具体计算过程如下：
假设在第k次迭代时，需要计算T″M,N与向量的乘积， 其中，n为系统状态变量个数，N为需要设定的子区间个数，N′＝N/α，α为放大倍数，M为每个子区间采用第二类切比雪夫多项式的M+1个经过位移和归一化处理的极点对每个子区间进行离散化的参数，具体步骤如下：
步骤(1)：由于T″M,N具有的逻辑结构，可知vk的第n+1～((N′-1)M+1)n个分量等于zk的第(M+1)n+1～(N′M+1)n个分量；
步骤(2)：将向量vk按列压缩为一矩阵 相应地，有：vk＝vec(Vk)，其中vec(·)为将矩阵压缩为列向量的运算；
步骤(3)：计算rk为中间向量；
步骤(4)：计算
8.如权利要求7所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，所述步骤(3)的步骤如下：
将B″M,N表达式代入，可得：
由上述分析可知，要计算rk，首先要计算然后再进行求和，然后将其压缩为一个(M+1)n维的列向量并乘以系数ejθ，再计算VkPT并将其压缩为一个(M+1)n维的列向量，最后将两个列向量相加。
9.如权利要求7所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，所述步骤(4)的步骤如下：
没有显示表达，因而，采用迭代算法来计算求解过程中，涉及矩阵向量乘法运算b＝R″My，其中
首先，将向量y按列压缩为一矩阵进而，得到b＝R″My的稀疏实现步骤如下：
稀疏实现。
10.如权利要求1所述的基于SOD-PS-II-R算法的时滞电力系统机电振荡模式计算方法，其特征是，在计算得到μ″之后，依次经过谱映射、坐标反旋转和牛顿校验后得到时滞电力系统的特征值λ，在牛顿校验前，λ的计算公式如下：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810776545.XA CN108808706B (zh) | 2018-07-13 | 2018-07-13 | 基于sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810776545.XA CN108808706B (zh) | 2018-07-13 | 2018-07-13 | 基于sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108808706A true CN108808706A (zh) | 2018-11-13 |
| CN108808706B CN108808706B (zh) | 2020-05-08 |
Family
ID=64076755
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810776545.XA Active CN108808706B (zh) | 2018-07-13 | 2018-07-13 | 基于sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108808706B (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109615209A (zh) * | 2018-12-05 | 2019-04-12 | 山东大学 | 一种时滞电力系统特征值计算方法及系统 |
| CN110298847A (zh) * | 2019-06-27 | 2019-10-01 | 浙江工业大学 | 一种长时间背景收集的背景建模方法 |
| CN112165085A (zh) * | 2020-07-31 | 2021-01-01 | 山东大学 | 基于psod的时滞电力系统高效特征值分析方法及系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH1132439A (ja) * | 1997-07-09 | 1999-02-02 | Mitsubishi Heavy Ind Ltd | 系統連系用電力変換器 |
| CN105468909A (zh) * | 2015-11-24 | 2016-04-06 | 山东大学 | 基于sod-ps-r算法的时滞电力系统机电振荡模式计算方法 |
-
2018
- 2018-07-13 CN CN201810776545.XA patent/CN108808706B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH1132439A (ja) * | 1997-07-09 | 1999-02-02 | Mitsubishi Heavy Ind Ltd | 系統連系用電力変換器 |
| CN105468909A (zh) * | 2015-11-24 | 2016-04-06 | 山东大学 | 基于sod-ps-r算法的时滞电力系统机电振荡模式计算方法 |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109615209A (zh) * | 2018-12-05 | 2019-04-12 | 山东大学 | 一种时滞电力系统特征值计算方法及系统 |
| CN109615209B (zh) * | 2018-12-05 | 2021-08-03 | 山东大学 | 一种时滞电力系统特征值计算方法及系统 |
| CN110298847A (zh) * | 2019-06-27 | 2019-10-01 | 浙江工业大学 | 一种长时间背景收集的背景建模方法 |
| CN110298847B (zh) * | 2019-06-27 | 2021-06-04 | 浙江工业大学 | 一种长时间背景收集的背景建模方法 |
| CN112165085A (zh) * | 2020-07-31 | 2021-01-01 | 山东大学 | 基于psod的时滞电力系统高效特征值分析方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108808706B (zh) | 2020-05-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108808706A (zh) | 基于sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 | |
| CN105468909A (zh) | 基于sod-ps-r算法的时滞电力系统机电振荡模式计算方法 | |
| CN104615882B (zh) | 基于eigd的大规模时滞电力系统特征值计算方法 | |
| CN104298809B (zh) | 一种基于矩阵指数电磁暂态仿真的非线性建模求解方法 | |
| Yang et al. | Real-time FPGA-RTDS co-simulator for power systems | |
| CN105449665B (zh) | 基于sod‑ps的时滞电力系统稳定性判别方法 | |
| CN102609598B (zh) | 一种对大型电力系统进行电磁暂态仿真的方法 | |
| Ye et al. | Efficient eigen-analysis for large delayed cyber-physical power system using explicit infinitesimal generator discretization | |
| WO2012031398A1 (en) | Numerical method for simulating subsonic flows based on euler equations in lagrangian formulation | |
| CN106410848A (zh) | 一种电力电子多馈入电力系统小干扰稳定性评估方法 | |
| CN108647906A (zh) | 基于低阶eigd的时滞电力系统稳定性分析方法 | |
| CN112214905B (zh) | 一种电力系统宽频带建模分析与仿真方法及系统 | |
| Diao et al. | On parallelizing single dynamic simulation using HPC techniques and APIs of commercial software | |
| CN100470995C (zh) | 电力系统特征值的分布式计算方法 | |
| CN105978003A (zh) | 一种考虑时滞的电力系统附加广域阻尼控制器设计方法 | |
| CN108471112A (zh) | 一种输电线路的电磁暂态仿真方法及系统 | |
| Klopfer et al. | A diagonalized diagonal dominant alternating direction implicit (D3ADI) scheme and subiteration correction | |
| CN108321821A (zh) | 基于sod-irk的时滞电力系统稳定性判别方法 | |
| CN110162843A (zh) | 一种电网一次系统联合二次系统实时仿真建模方法与装置 | |
| CN104950261B (zh) | 电池的硬件在环仿真测试方法及系统 | |
| CN108242808A (zh) | 基于igd-lms的时滞电力系统稳定性判别方法 | |
| CN108808705A (zh) | 基于低阶sod-ps-ii-r算法的时滞电力系统机电振荡模式计算方法 | |
| Francescatto et al. | A semi‐coarsening strategy for unstructured multigrid based on agglomeration | |
| CN109685400A (zh) | 基于时间积分igd的时滞电力系统稳定性判别方法 | |
| CN108879669B (zh) | 基于低阶iigd算法的时滞电力系统特征值分析方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |