CN104615882B - 基于eigd的大规模时滞电力***特征值计算方法 - Google Patents
基于eigd的大规模时滞电力***特征值计算方法 Download PDFInfo
- Publication number
- CN104615882B CN104615882B CN201510055743.3A CN201510055743A CN104615882B CN 104615882 B CN104615882 B CN 104615882B CN 201510055743 A CN201510055743 A CN 201510055743A CN 104615882 B CN104615882 B CN 104615882B
- Authority
- CN
- China
- Prior art keywords
- matrix
- formula
- time
- power system
- vector
- 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.)
- Active
Links
Landscapes
- Supply And Distribution Of Alternating Current (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了基于EIGD的大规模时滞电力***特征值计算方法,步骤如下:建立时滞电力***模型,进行时滞电力***模型的公式转化,转换为抽象柯西问题,进而将计算时滞电力***模型的特征值转化成计算转化后公式的无穷小生成元的特征值;对无穷小生成元进行离散化,得到无穷小生成元的近似矩阵;近似矩阵通过位移处理,得到位移处理后的近似矩阵,位移处理后的近似矩阵经过逆变换,得到逆矩阵;从而将要求计算的部分特征值变为模值最大的部分特征值;采用Arnoldi算法来计算得到的逆矩阵的模值最大的部分特征值;从而得到近似矩阵的特征值为;利用牛顿迭代法进行校验,通过计算得到时滞电力***精确的特征值和特征向量。
Description
技术领域
本发明涉及一种基于EIGD的大规模时滞电力***特征值计算方法,EIGD的中文解释:显示无穷小生成元离散化,EIGD的英文全称:Explicit Infinitesimal GeneratorDiscretization。
背景技术
广域测量***(Wide-Area Measurement System,WAMS)的出现给大规模互联电力***稳定分析与控制的发展带来新的契机。基于WAMS提供的广域信息的互联电网低频振荡控制,通过引入有效反映区间振荡模式的广域反馈信号,能够获得较好的阻尼控制性能,其为解决互联电网中的区域间低频振荡问题,进而提高***的输电能力提供了新的控制手段,具有良好而又广泛的应用前景。
广域信号在由不同通信介质(如光纤、电话线、数字微波、卫星等)组成的WAMS通信网络中传输和处理时,存在几十到几百毫秒间变化的通信延时。时滞是导致***控制律失效、运行状况恶化和***失稳的一种重要诱因[Wu H X,Tsakalis K S,Heydt GT.Evaluation of time delay effects to wide-area power system stabilizerdesign.IEEE Trans.Power Syst.,2004,19(4):1935-1941.]。因此,利用广域测量信息进行电力***闭环控制时,必须计及时滞的影响。
目前,考虑广域通信时滞的电力***小干扰稳定性分析方法,总体上可分为时域方法和频率方法两大类。
(1)时域方法
时域方法,主要基于时滞依赖的Lyapunov-Krasovskii稳定定理和Razumikhin定理,并利用线性矩阵不等式(Linear Matrix Inequality,LMI)技术来设计鲁棒控制器,在保证***稳定性的同时,能够得到***能够承受的时滞上限。时域方法的缺点主要有以下两点:(1)基于Lyapunov的时滞依赖的稳定性判据,仅为***稳定的充分条件,存在固有的保守性;(2)当应用于大规模电力***时,需要与模型降阶技术相结合以降低计算量。但是,***模型的精确性对稳定性的判别存在一定的影响。
(2)频域方法
●函数变换方法
在频域范围内,当考虑时滞影响后,电力***线性化微分-代数方程对应的特征方程中标量时滞被转化为指数项。因而,特征方程是一个超越方程。为此,在求取电力***的时滞稳定裕度时,通常采用函数变换的方法(如Rekasius变换[Rifat S,Nejat O.A novelstability study on multiple time-delay systems(MTDS)using the root clusteringparadigm.Proc of the American Control Conference,Boston,MA,2004,5422-5427.]、Lambert-W函数[Yi S,Nelson P W,Ulsoy A G.Time-delay systems:analysis andcontrol using the Lambert W function.Singapore:World Scientific PublishingCompany,2010:]、Pade近似,等)对超越项进行变换,避免直接求解特征方程的困难。
函数变换方法存在如下不足:(1)Rekasius变换方法,可以将超越的特征方程改写为一个代数方程,从而能够计算得到虚轴上所有特征值。然而,该方法不能得到复平面上其他特征值的分布。(2)当***存在对称时滞且***矩阵能够同时被三角化时,时滞***的特征谱可以利用Lambert-W函数显示表示。可见,该方法存在较强的前提和假设条件。(3)Pade近似是一种常用的时滞环节处理方法。通过有理多项式来逼近时滞环节,进而可以方便地利用经典和现代控制理论设计广域阻尼控制律,并通过时域仿真验证控制器的有效性。然而,Pade近似的准确性对阻尼控制器设计的影响却未被深入分析和研究。
●谱离散化方法
近年来,在数值分析和计算数学领域,谱离散化方法大量研究用于计算时滞***的部分特征值。文献[Milano F,Anghel M.Impact of time delays on power systemstability.IEEE Trans.Circuit Syst.-I,2012,59(4):889-900.]利用无穷小生成元离散化方法计算励磁***中发电机端电压量测环节存在单一时滞情况下,***的部分特征值。文献[Liang H,Choi B J,Zhuang W H,等.Stability enhancement of decentralizedInverter control through wireless communications in microgrids.IEEETrans.Smart Grid,2013,4(1):321-331.]利用无穷小生成元离散化方法分析在廉价通信方案下,当微电网中逆变器控制环节存在时滞时,***的小干扰稳定性。
谱离散化方法存在如下不足:(1)无穷小生成元的近似矩阵,没有清晰、解析的表达形式,难以利用稀疏特征值方法计算得到***的部分特征值;(2)没有利用***线性化增广矩阵的稀疏特性,计算量较大,不能用于大规模电力***时滞特征值的计算。
发明内容
本发明的目的就是为了解决上述问题,提供一种基于EIGD的大规模时滞电力***特征值计算方法,它具有计算量小,内存要求低,计算速度快,计算效率高的优点,能够用于考虑广域时滞影响的大规模电力***的特征值计算和时滞稳定性判别。
为了实现上述目的,本发明采用如下技术方案:
基于EIGD的大规模时滞电力***特征值计算方法,包括如下步骤:
步骤(1):建立时滞电力***模型,依据时滞电力***模型的状态变量在设定区间段内是连续可微的,进行时滞电力***模型的公式转化,转换为抽象柯西问题,进而将计算时滞电力***模型的特征值转化成计算转化后公式的无穷小生成元的特征值;为了得到无穷小生成元的特征值,对无穷小生成元进行离散化,从而得到无穷小生成元的近似矩阵
步骤(2):近似矩阵通过位移处理,得到位移处理后的近似矩阵位移处理后的近似矩阵经过逆变换,得到逆矩阵从而将复平面上距离某个设定点最近的、设定数量个数的特征值,变换为设定数量个数的模值最大特征值;
步骤(3):采用隐式重启动Arnoldi算法来计算步骤(2)得到的逆矩阵的模值最大特征值λ";从而得到近似矩阵的特征值为λ;
步骤(4):利用牛顿迭代法对步骤(3)得到的特征值为λ进行校验,通过计算得到时滞电力***精确的特征值和特征向量。
所述步骤(1)的时滞电力***模型如下:
式中:t为当前时刻;Δx(t)∈Rn×1为增量形式的电力***状态变量向量,n为***状态变量总数;表示电力***状态变量向量的导数;τi>0为第i个时滞环节的时滞常数,i=1,2,…,m;其中,最大的时滞表示为τmax;为***状态矩阵,为一个稠密矩阵;(i=1,…,m)为***时滞状态矩阵,为稀疏矩阵;
稠密矩阵和稀疏矩阵具体表示为:
式中:Αi∈Rn×n,Bi∈Rn×l,Ci∈Rl×n,Di∈Rl×l为线性化***矩阵,l为***代数变量总数。
所述步骤(1)的抽象柯西问题:
式中:u(t)为一个连续可微函数,且u(t)=x(t+θ),θ∈[-τmax,0]。
所述步骤(1)的无穷小生成元的近似矩阵的计算公式为:
式中:ΠN∈R(N+1)n×(N+1)n为伴随矩阵,ΣN∈R(N+1)n×(N+1)n为块上三角矩阵;
伴随矩阵ΠN和块上三角矩阵ΣN具体表示为:
式中:表示克罗内克积;Rj(j=0,1,…,N)由***状态矩阵经过代数运算得到,即:
所述步骤(2)的位移处理之前需要做的步骤如下:
式(1)表示的线性化***的特征方程为:
式中:λ为特征值,v为特征值对应的右特征向量;
式(2)的等价增广形式为:
式中:w∈Rl×1为中间和辅助向量;
设In为n阶单位矩阵,则A'(λ)和B'(λ)具体表示为:
所述步骤(2)的位移处理步骤如下:
给定位移点s,将式(2)中的λ用λ'+s代替,则得到位移后的特征方程:
式中:
λ′=λ-s
所述步骤(2)中:
设表示经过位移处理后的近似矩阵;逆矩阵的计算过程为
式中:(Σ′N)-1显示地表示如下:
(Σ′N)-1=[Σ′0;[0Nn,INn]]
Σ′0=(R0′)-1[In,R1′,…,R′N]
式中:0Nn∈RNn×n为零矩阵,INn∈RNn×n为单位矩阵;[;]和[,]分别表示在列和行的方向上进行向量或矩阵连接;A'(s)和B'(s)将式(4)和(5)中的λ用s代替后得到。
所述步骤(3)的步骤如下:
设第k个Krylov向量表示为qk,则矩阵向量乘积的流程如下:
步骤(3-1):利用稀疏矩阵乘法,计算得到
步骤(3-2):从列的方向上,将r压缩为矩阵其中
步骤(3-3):赋值z=v0;
步骤(3-4):赋值j=1;
步骤(3-5):计算
步骤(3-6):计算
步骤(3-7):求和z=z+p+w;
步骤(3-8):将j增加1,重复步骤(3-5)-(3-7),直至j=N;
步骤(3-9):计算qk+1(1:n)=(R0′)-1z;
步骤(3-10):qk+1((n+1):(N+1)n)=r((n+1):(N+1)n)。
在步骤(3-10)中,qk+1(i1:i2)和r(i1:i2)分别表示抽取qk+1和r的第i1到i2个元素形成的列向量。
设利用隐式重启动Arnoldi算法计算得到的特征值为λ",则稀疏近似矩阵对应的特征值为λ=1/λ"+s=λ′+s。此外,与λ对应的Krylov向量的前n个分量是特征向量v的良好的估计和近似。
所述步骤(4)的步骤为:
将式(3)线性化,得到牛顿迭代的修正方程式如下:
式中:λ(k)为第k次迭代时特征值,v(k)为第k次迭代时特征向量,w(k)为第k次迭代时中间辅助向量的估计值;Δλ(k)为第k次迭代时特征值的修正量;Δv(k)为第k次迭代时特征向量的修正量;Δw(k)为第k次迭代时中间辅助向量的估计值的修正量;f(k)表示第k次迭代时增广形式特征方程的不平衡量。
为了保证特征向量v(k)的唯一性,假设其第一个分量并在牛顿迭代过程中保持不变,相应地,对v(k)中的其他元素进行规格化,为了避免式(10)表示的修正方程降阶并保留其稀疏特性,将Δv(k)的第一个分量用Δλ(k)代替。相应地,用g'(k)代替J'(k)的第一列,则得到修正后的雅可比矩阵为J"(k)=[g'(k),J'(k)(1:(n+l),2:(n+l))];
经过上述处理后,修正方程式(10)变为:
式中:J(k)为将J"(k)的第一列用单位向量e1代替后得到的矩阵,其与J'(k)具有相同的稀疏结构;g(k)为g'(k)与e1之差;
利用Sherman-Morrisony公式,得式(11)的解:
牛顿迭代收敛以后,v′(k)的第一个元素就是待求的时滞电力***的精确特征值λ。
所述步骤(4)如果计算得到的全部特征值λ都具有负实部,则时滞电力***在稳态运行点处是小扰动稳定的;反之,若至少存在一个具有正实部的特征值,则***在该点处是小扰动不稳定的。
本发明的有益效果:
1由上述流程看出,步骤(3-5)和步骤(3-9)的计算量最大、耗时最长。从而可以推知,整个流程的计算量近似正比于集合ΩN中的离散点数N。
计算化无穷小生成元的近似矩阵和它的逆矩阵所具有的显示、稀疏特性,使得可以充分利用***线性化增广系数矩阵Αi,Bi,Ci,Di(i=0,1,…,m)的稀疏特性,以稀疏的方式实现步骤(3-5)和(3-9),从而大大减小计算量和内存要求,提高了计算效率,能够用于考虑广域时滞影响的大规模电力***特征值计算和时滞稳定性判别。
2能够高效、准确地计算得到,用于判别电力***时滞稳定性判别所需的、复平面上最右侧的部分特征值。此外,由于不需要计算出***的所有特征值,就可以实现电力***时滞稳定性的判别。总起来说,该方法具有很高的计算效率。
附图说明
图1为本发明的流程图;
图2为四机两区域算例***;
图3为N=50时,和的特征值;
图4为N=50时,和的特征值。
具体实施方式
下面结合附图与实施例对本发明作进一步说明。
1.时滞电力***模型
考虑广域通信时滞影响后,电力***可以如下的一组时滞微分方程组来描述:
式中:x∈Rn×1为电力***的状态变量向量,n为***状态变量总数。t为当前时刻,τi>0为第i个时滞环节的时滞常数,i=1,2,…,m。其中,最大的时滞表示为τmax。为***状态矩阵,为一个稠密矩阵;为***时滞状态矩阵,为稀疏矩阵。它们可具体表示为:
式中:Αi∈Rn×n,Bi∈Rn×l,Ci∈Rl×n,Di∈Rl×l为线性化***矩阵,l为***代数变量总数。
式(1)表示的线性化***的特征方程为:
式中:λ和v分别为特征值和相应的右特征向量。
式(2)的等价增广形式为:
式中:w∈Rl×1为中间和辅助向量。设In为n阶单位矩阵,则A'(λ)和B'(λ)可具体表示为:
2.基于显示无穷小生成元离散化的时滞特征值计算方法------基本理论
由于时滞电力***的状态变量x(t)在区间[-τmax,∞)上是连续可微的,可以证明式(1)可以转化为下列抽象柯西问题:
式中:u(t)为一个连续可微函数,且u(t)=x(t+θ),θ∈[-τmax,0]。
由谱映射定理可知,式(1)表示的时滞***的特征值等于无穷小生成元的特征值,因此可以通过计算的特征值,来得到时滞***的特征值。然而,的特征值本质上是一个无穷维问题。因此,需要对进行离散化,从而得到其近似矩阵
给定正整数N,区间[-τmax,0]上N+1个不同离散点形成的集合ΩN定义为ΩN:={θN,j,j=1,…,N+1},且有-τmax≤θN,j≤…≤θN,N+1=0。集合ΩN中的N个非零离散点,选取为经过位移和归一化处理后的第N阶第一类切比雪夫多项式的零点。在XN上,式(6)表示的抽象柯西问题可以被离散化,进而可以得到无穷小生成元的近似矩阵
式中:ΠN∈R(N+1)n×(N+1)n为伴随矩阵,ΣN∈R(N+1)n×(N+1)n为块上三角矩阵。ΠN和ΣN可具体表示为:
式中:表示克罗内克积。Rj(j=0,1,…,N)由***状态矩阵经过代数运算得到,即:
3.基于显示无穷小生成元离散化的时滞特征值计算方法------稀疏实现
3.1位移逆变换
电力***机电振荡模式对应的特征值,通常位于复平面虚轴附近,决定着电力***的小干扰稳定性。与其他特征值相比,这部分特征值的模值较小。因此,需要利用位移逆变换将近似矩阵从而将要求计算的部分特征值变为模值最大的部分特征值。
给定位移点s,将式(2)中的λ用λ'+s代替,则可以得到位移后的特征方程:
式中:
λ′=λ-s
设表示经过位移处理后的近似矩阵。其逆矩阵为:
式中:(Σ′N)-1可以显示地表示如下:
(Σ′N)-1=[Σ′0;[0Nn,INn]]
Σ′0=(R0′)-1[In,R1′,…,R′N]
式中:0Nn∈RNn×n为零矩阵,INn∈RNn×n为单位矩阵。[;]和[,]分别表示在列和行的方向上进行向量或矩阵连接。A'(s)和B'(s)可以将式(4)和(5)中的λ用s代替后得到。
3.2矩阵向量乘积的稀疏实现
采用隐式重启动Arnoldi算法来计算经过位移逆变换后近似矩阵模值最大的部分特征值。在Arnoldi算法中,耗时最长、计算量最大的操作就是在形成Krylov向量过程中的矩阵向量乘积运算。设第k个Krylov向量表示为qk,则矩阵向量乘积的流程如下:
步骤1:利用稀疏矩阵乘法,计算得到
步骤2:从列的方向上,将r压缩为矩阵其中
步骤3:赋值z=v0;
步骤4:赋值j=1;
步骤5:计算
步骤6:计算
步骤7:求和z=z+p+w;
步骤8:将j增加1,然后重复步骤5-步骤7,直至j=N;
步骤9:计算qk+1(1:n)=(R0′)-1z;
步骤10:qk+1((n+1):(N+1)n)=r((n+1):(N+1)n)。
在步骤10中,qk+1(i1:i2)和r(i1:i2)分别表示抽取qk+1和r的第i1到i2个元素形成的列向量。
由上述流程可以看出,步骤5和步骤9的计算量最大、耗时最长。从而可以推知,整个流程的计算量近似正比于集合ΩN中的离散点数N。
利用***线性化增广系数矩阵的稀疏特性,步骤5和9可以以稀疏的方式实现,从而大大减小计算量,提高计算效率。
设利用隐式重启动Arnoldi算法计算得到的特征值为λ",则稀疏近似矩阵对应的特征值为λ=1/λ"+s=λ′+s。此外,与λ对应的Krylov向量的前n个分量是特征向量v的良好的估计和近似。
3.3牛顿校验
为了改善计算得到的特征值和特征向量的精度,可以利用牛顿法,计算得到时滞电力***精确的特征值和相应的特征向量。
将式(3)线性化,得到牛顿迭代的修正方程式如下:
式中:λ(k)、v(k)和w(k)分别为第k次迭代时特征值、特征向量和中间辅助向量的估计值。前缀Δ表示它们各自的修正量。f(k)表示第k次迭代时增广形式特征方程的不平衡量。
为了保证特征向量v(k)的唯一性,假设其第一个分量并在牛顿迭代过程中保持不变。相应地,对v(k)中的其他元素进行规格化,其值可能大于1。为了避免式(10)表示的修正方程降阶并保留其稀疏特性,将Δv(k)的第一个分量用Δλ(k)代替。相应地,用g'(k)代替J'(k)的第一列,则得到修正后的雅可比矩阵为J"(k)=[g'(k),J'(k)(1:(n+l),2:(n+l))]。经过上述处理后,修正方程式(10)变为:
式中:J(k)为将J"(k)的第一列用单位向量e1代替后得到的矩阵,其与J'(k)具有相同的稀疏结构。g(k)为g'(k)与e1之差。
利用Sherman-Morrisony公式,可得式(11)的解:
牛顿迭代收敛以后,v′(k)的第一个元素就是待求的时滞电力***的精确特征值λ。
如果计算得到的全部特征值λ都具有负实部,则时滞电力***在稳态运行点处是小扰动稳定的;反之,若至少存在一个具有正实部的特征值,则***在该点处是小扰动不稳定的。
4.算例验证
利用四机两区算例***,来验证本文提出的基于显示无穷小生成元离散化的大规模时滞电力***的稀疏特征值计算方法的有效性。所有的分析均在Matlab中和在InterCore i5 4×3.4GHz 8GB RAM个人计算机上进行。
四机两区算例***如图2所示。所有发电机采用高增益晶闸管励磁***,并装设以转速偏差为输入的电力***稳定器(PSS)。在此基础上,在G1上装设以G1和G3相对转速偏差Δω13为反馈信号的广域PSS,以进一步提高对区域间低频振荡模式的阻尼。。假设广域PSS的反馈时滞和输入时滞分别为τ1=150ms和τ2=100ms。***状态变量和代数变量的维数分别为n=56和l=22。
4.1稀疏近似矩阵的逼近能力分析
为了验证本发明提出的基于显示无穷小生成元离散化的大规模时滞电力***的稀疏特征值计算方法的有效性,首先分析无穷小生成元离散化得到的稀疏近似矩阵逼近无穷小生成元的能力和准确性。在四机两区***中,当N=50时,的维数为(N+1)n=2856。首先,利用QR算法计算得到的全部特征值,如图3所示。接着,将这些特征值作为牛顿法的初始值,进行校验。其中一部分特征值能够收敛为无穷小生成元的特征值,如图3所示。
下面对和的特征值进行比较和分析。由图3可知,在Re(λ)∈[-48,0]、Im(λ)∈[0,820]的区域内,以及点(-100,0)附近,和的特征值能较好地吻合。图3中,剩余特征值为的特征值,经过牛顿迭代后,它们不能收敛为的相应的特征值。但是,它们均位于已经收敛的那部分特征值的左侧。通过上述分析可知,在一定程度上能较准确地逼近但这对于电力***小干扰稳定性分析而言已经足够。
4.2方法的准确性分析
利用隐式重启动Arnoldi算法计算稀疏近似矩阵的部分特征值。由于四机两区算例***的规模较小,可以通过为算法选择合适的参数,从而一次性得到实部大于某个数值、***最右侧的部分特征值。当N=50时,位移点选择为s=0.05,计算得到其附近的r=200个特征值,如4所示。由图可知,Arnoldi算法能够准确计算得到频率近似矩阵最大的特征值非常接近无穷小生成元的第19个特征值λ19。
下面对本发明所提出的方法和文献[Breda D,Maset S,VermiglioR.Pseudospectral differencing methods for characteristic roots of delaydifferential equations.SIAM Journal on Scientific Computing,2005,27(2):482-495.]提出的基于隐式无穷小生成元离散化的特征值计算方法(Implicit InfinitesimalGenerator Discretization,IIGD方法)进行对比。利用IIGD方法计算得到s=0.05附近200个特征值,如图4所示。与本发明提出的基于显示无穷小生成元离散化的特征值计算方法(EIGD)的结果相比,差别在于IIGD方法能够准确的求得的频率最大的特征值为λ20,比EIGD求得的特征值多一个。
1.1计算时间分析
四种情况下,即N=20,40,50和60时,利用QR计算近似矩阵的全部特征值,以及利用非稀疏IGD和稀疏IGD方法分别计算近似矩阵和在位移点s=0.05附近r=200个特征值。计算表1所示。
表1 QR、IIGD和EIGD算法所用的计算时间比较
N=20 | N=40 | N=50 | N=60 | |
QR | 1.97s | 13.71s | 27.06s | 47.70s |
IIGD | 10.14s | 25.84s | 34.96s | 45.81s |
EIGD | 6.24s | 10.61s | 12.69s | 15.21s |
首先,当N=20时,虽然QR算法计算了稀疏近似矩阵的全部特征值,但是由于其在***规模较小的情况下非常鲁棒且收敛收敛速度快,其比EIGD算法计算200个特征值所用的时间少得多。当N增加时,QR算法的计算量和计算时间急剧增加。相比较而言,EIGD算法的计算时间增加比较平缓。这是因为,在EIGD方法中,矩阵向量乘积的计算量与N大体上呈线性关系。其次,当N=40,50和60时,EIGD方法所需的计算时间比QR和IIGD方法少。
上述虽然结合附图对本发明的具体实施方式进行了描述,但并非对本发明保护范围的限制,所属领域技术人员应该明白,在本发明的技术方案的基础上,本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。
Claims (10)
1.基于EIGD的大规模时滞电力***特征值计算方法,其特征是,包括如下步骤:
步骤(1):建立时滞电力***模型,依据时滞电力***模型的状态变量在设定区间段内是连续可微的,进行时滞电力***模型的公式转化,转换为抽象柯西问题,进而将计算时滞电力***模型的特征值转化成计算转化后公式的无穷小生成元的特征值;为了得到无穷小生成元的特征值,对无穷小生成元进行离散化,从而得到无穷小生成元的近似矩阵
步骤(2):近似矩阵通过位移处理,得到位移处理后的近似矩阵位移处理后的近似矩阵经过逆变换,得到逆矩阵从而将复平面上距离某个设定点最近的、设定数量个数的特征值,变换为设定数量个数的模值最大特征值;
步骤(3):采用隐式重启动Arnoldi算法来计算步骤(2)得到的逆矩阵的模值最大特征值λ";从而得到近似矩阵的特征值为λ;
步骤(4):利用牛顿迭代法对步骤(3)得到的特征值λ进行校验,通过计算得到时滞电力***精确的特征值和特征向量。
2.如权利要求1所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(1)的时滞电力***模型如下:
式中:t为当前时刻;Δx(t)∈Rn×1为增量形式的电力***状态变量向量,n为***状态变量总数;表示电力***状态变量向量的导数;τi>0为第i个时滞环节的时滞常数,i=1,2,…,m;其中,最大的时滞表示为τmax;为***状态矩阵,为一个稠密矩阵;i=1,…,m为***时滞状态矩阵,为稀疏矩阵;
稠密矩阵和稀疏矩阵具体表示为:
式中:Αi∈Rn×n,Bi∈Rn×l,Ci∈Rl×n,Di∈Rl×l为线性化***矩阵,l为***代数变量总数。
3.如权利要求2所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(1)的抽象柯西问题:
式中:u(t)为一个连续可微函数,且u(t)=x(t+θ),θ∈[-τmax,0]。
4.如权利要求3所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(1)的无穷小生成元的近似矩阵的计算公式为:
式中:ΠN∈R(N+1)n×(N+1)n为伴随矩阵,ΣN∈R(N+1)n×(N+1)n为块上三角矩阵;
伴随矩阵ΠN和块上三角矩阵ΣN具体表示为:
式中:表示克罗内克积;Rj由***状态矩阵和***时滞状态矩阵经过代数运算得到,i=1,…,m,j=0,1,…,N;即:
5.如权利要求4所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(2)的位移处理之前需要做的步骤如下:
式(1)表示的线性化***的特征方程为:
式中:λ为特征值,v为特征值对应的右特征向量;
式(2)的等价增广形式为:
式中:w∈Rl×1为中间和辅助向量;
设In为n阶单位矩阵,则A'(λ)和B'(λ)具体表示为:
6.如权利要求5所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(2)的位移处理步骤如下:
给定位移点s,将式(2)中的λ用λ'+s代替,则得到位移后的特征方程:
式中:
λ′=λ-s
7.如权利要求6所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(2)中:
设表示经过位移处理后的近似矩阵;逆矩阵的计算过程为
式中:(Σ′N)-1显示地表示如下:
(Σ′N)-1=[Σ′0;[0Nn,INn]]
Σ′0=(R′0)-1[In,R′1,…,R′N]
式中:0Nn∈RNn×n为零矩阵,INn∈RNn×n为单位矩阵;[;]和[,]分别表示在列和行的方向上进行向量或矩阵连接;A'(s)和B'(s)将式(4)和(5)中的λ用s代替后得到。
8.如权利要求7所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(3)的步骤如下:
设第k个Krylov向量表示为qk,则矩阵向量乘积的流程如下:
步骤(3-1):利用稀疏矩阵乘法,计算得到
步骤(3-2):从列的方向上,将r压缩为矩阵其中
步骤(3-3):赋值z=v0;
步骤(3-4):赋值j=1;
步骤(3-5):计算
步骤(3-6):计算
步骤(3-7):求和z=z+p+w;
步骤(3-8):将j增加1,重复步骤(3-5)-(3-7),直至j=N;
步骤(3-9):计算qk+1(1:n)=(R′0)-1z;
步骤(3-10):qk+1((n+1):(N+1)n)=r((n+1):(N+1)n);
在步骤(3-10)中,qk+1(i1:i2)和r(i1:i2)分别表示抽取qk+1和r的第i1到i2个元素形成的列向量;
设利用隐式重启动Arnoldi算法计算得到的特征值为λ",则稀疏近似矩阵对应的特征值为λ=1/λ"+s=λ′+s;此外,与λ对应的Krylov向量的前n个分量是特征向量v的良好的估计和近似。
9.如权利要求8所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(4)的步骤为:
将式(3)线性化,得到牛顿迭代的修正方程式如下:
式中:λ(k)为第k次迭代时特征值,v(k)为第k次迭代时特征向量,w(k)为第k次迭代时中间辅助向量的估计值;Δλ(k)为第k次迭代时特征值的修正量;Δv(k)为第k次迭代时特征向量的修正量;Δw(k)为第k次迭代时中间辅助向量的估计值的修正量;f(k)表示第k次迭代时增广形式特征方程的不平衡量;
为了保证特征向量v(k)的唯一性,假设其第一个分量在牛顿迭代过程中保持不变,相应地,对v(k)中的其他元素进行规格化,为了避免式(10)表示的修正方程降阶并保留其稀疏特性,将Δv(k)的第一个分量用Δλ(k)代替;相应地,用g'(k)代替J'(k)的第一列,则得到修正后的雅可比矩阵为J"(k)=[g'(k),J'(k)(1:(n+l),2:(n+l))];
经过上述处理后,修正方程式(10)变为:
式中:J(k)为将J"(k)的第一列用单位向量e1代替后得到的矩阵,其与J'(k)具有相同的稀疏结构;g(k)为g'(k)与e1之差;
利用Sherman-Morrisony公式,得式(11)的解:
牛顿迭代收敛以后,v′(k)的第一个元素就是待求的时滞电力***的精确特征值λ。
10.如权利要求1所述的基于EIGD的大规模时滞电力***特征值计算方法,其特征是,所述步骤(4)如果计算得到的全部特征值λ都具有负实部,则时滞电力***在稳态运行点处是小扰动稳定的;反之,若至少存在一个具有正实部的特征值,则***在该点处是小扰动不稳定的。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510055743.3A CN104615882B (zh) | 2015-02-03 | 2015-02-03 | 基于eigd的大规模时滞电力***特征值计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510055743.3A CN104615882B (zh) | 2015-02-03 | 2015-02-03 | 基于eigd的大规模时滞电力***特征值计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104615882A CN104615882A (zh) | 2015-05-13 |
CN104615882B true CN104615882B (zh) | 2017-06-27 |
Family
ID=53150323
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510055743.3A Active CN104615882B (zh) | 2015-02-03 | 2015-02-03 | 基于eigd的大规模时滞电力***特征值计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104615882B (zh) |
Families Citing this family (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104809666B (zh) * | 2015-05-20 | 2017-12-29 | 华北电力大学 | 一种搜索圆半径可控的电力***小干扰特征值计算方法 |
CN105468909B (zh) * | 2015-11-24 | 2018-01-30 | 山东大学 | 基于sod‑ps‑r算法的时滞电力***机电振荡模式计算方法 |
CN109685400B (zh) * | 2018-02-24 | 2020-07-31 | 山东大学 | 基于时间积分igd的时滞电力***稳定性判别方法 |
CN108647906B (zh) * | 2018-05-25 | 2020-05-08 | 山东大学 | 基于低阶eigd的时滞电力***稳定性分析方法 |
CN108879669B (zh) * | 2018-07-13 | 2021-03-23 | 山东大学 | 基于低阶iigd算法的时滞电力***特征值分析方法 |
CN109033022B (zh) * | 2018-07-13 | 2020-07-31 | 山东大学 | 基于低阶sod-lms算法的时滞电力***特征值计算方法 |
CN108808703B (zh) * | 2018-07-13 | 2020-07-31 | 山东大学 | 基于低阶igd-irk的时滞电力***小干扰稳定性分析方法 |
CN108923421B (zh) * | 2018-07-13 | 2020-04-21 | 山东大学 | 基于低阶sod-irk算法的时滞电力***关键特征值计算方法 |
CN109583065B (zh) * | 2018-11-20 | 2022-10-28 | 西安交通大学 | 基于部分变量离散化的时滞电力***特征值计算方法 |
CN109615209B (zh) * | 2018-12-05 | 2021-08-03 | 山东大学 | 一种时滞电力***特征值计算方法及*** |
CN112365115A (zh) * | 2020-08-26 | 2021-02-12 | 天津大学 | 一种电网信息能源***稳定性评估方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101777767B (zh) * | 2010-03-15 | 2012-07-25 | 天津大学 | 一种时滞电力***稳定的判别方法 |
CN102801158B (zh) * | 2012-07-31 | 2014-07-09 | 国网山东省电力公司经济技术研究院 | 基于Pade近似的时滞电力***特征值计算与稳定性判别方法 |
CN103838965B (zh) * | 2014-02-26 | 2017-01-04 | 华北电力大学 | 基于广义特征值的时滞稳定上限计算***及其计算方法 |
-
2015
- 2015-02-03 CN CN201510055743.3A patent/CN104615882B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN104615882A (zh) | 2015-05-13 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104615882B (zh) | 基于eigd的大规模时滞电力***特征值计算方法 | |
Khazraj et al. | A performance comparison between extended Kalman Filter and unscented Kalman Filter in power system dynamic state estimation | |
CN106532711B (zh) | 随迭代和节点类型改变雅可比矩阵的牛顿法潮流计算方法 | |
Liu et al. | Solving power system differential algebraic equations using differential transformation | |
CN102801158B (zh) | 基于Pade近似的时滞电力***特征值计算与稳定性判别方法 | |
CN105449665B (zh) | 基于sod‑ps的时滞电力***稳定性判别方法 | |
Zhuravlev et al. | Calculation of electric drives with electric machines of unconventional design | |
CN105468909B (zh) | 基于sod‑ps‑r算法的时滞电力***机电振荡模式计算方法 | |
CN105701568B (zh) | 一种启发式的配电网状态估计量测位置快速优化方法 | |
US11569682B2 (en) | System and method for a fast power network simulator | |
Zhang et al. | Towards highly efficient state estimation with nonlinear measurements in distribution systems | |
CN111860950B (zh) | 一种虚拟电厂功率可调能力的概率评估方法 | |
CN105184027A (zh) | 一种基于交互式多模型算法的电力负荷建模方法 | |
CN107069696A (zh) | 一种电力***状态估计的并行计算方法 | |
CN108647906A (zh) | 基于低阶eigd的时滞电力***稳定性分析方法 | |
CN108054757A (zh) | 一种内嵌无功和电压的n-1闭环安全校核方法 | |
CN109617080A (zh) | 基于改进的雅可比矩阵的直角坐标牛顿法潮流计算方法 | |
CN107834593A (zh) | 一种下垂控制型孤岛微电网静态电压稳定概率评估方法 | |
CN109698505B (zh) | 大电网静态电压稳定在线防控的调控量化映射计算方法 | |
CN109494748A (zh) | 基于节点类型和修正的雅可比矩阵的牛顿法潮流计算方法 | |
CN106410811B (zh) | 首次迭代小阻抗支路端点改变雅可比矩阵的潮流计算方法 | |
Wang et al. | Probabilistic static voltage stability analysis considering the correlation of wind power | |
US20220244687A1 (en) | Stability region determining method for distributed power cyber-physical system with multiple time delays | |
CN108808706A (zh) | 基于sod-ps-ii-r算法的时滞电力***机电振荡模式计算方法 | |
CN108242808A (zh) | 基于igd-lms的时滞电力***稳定性判别方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |