基于Pade近似的时滞电力系统特征值计算与稳定性判别方法

摘要

本发明公开了一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，它利用Pade近似将广域反馈时滞逼近为一个有理多项式，通过与无时滞电力系统和广域阻尼控制器之间的连接，建立时滞电力系统的的线性化模型，最后根据系统状态矩阵直接求出时滞系统的部分特征根，进而判别系统的时滞稳定性。四机两区域算例系统的特征值计算结果，该方法能够较准确地求解与时滞系统中动态元件相关的部分特征值和特征向量，可以方便准确地判别系统的时滞稳定性，正确求解与时滞环节相关的特征值个数和计算精度，与有理多项式的阶数有关。此外，该方法具有计算量小，计算时间短的优点。

说明书

技术领域

本发明涉及一种时滞电力系统特征值计算与稳定性判别方法，尤其涉及一种基于Pade近 似的时滞电力系统特征值计算与稳定性判别方法。

背景技术

随着电力系统规模的扩大，区域间的低频振荡正成为限制电网传输能力的瓶颈，而现有 的阻尼控制器（主要是按照单机无穷大系统设计的电力系统稳定器）并不能很好地解决这一 问题。根本原因在于：(1)不能直接利用相对功角和相对角速度构成闭环控制。虽然采用相 对功角和相对角速度来实现阻尼控制是最直接和有效的，但长期以来缺少必要的测量方法， 只能采用其他间接变量来代替，导致控制效果不佳。(2)限于本地局部信息。采用本地测量 构成反馈控制，不能很好地反映区间振荡模式，导致控制系统虽然能阻尼本地振荡模式，但 难以有效地抑制区间振荡模式。

随着信息处理和通信技术的迅猛发展，同步相量技术和广域测量系统（Wide Area  Measurement System,WAMS），给电力系统的监测、分析和控制提供了新的手段，为互联电网 阻尼控制带来了新的契机。同步相量测量单元可同步采集表征电网运行状态的几乎所有的变 量，最为关键的是它能测量发电机的内电势、转子角、角速度，母线电压相位等与低频振荡 密切相关的量。在高速通信网络（如电力数据宽带网）的支持下，各相量测量单元采集的带 时标的数据能以较小的延时传递到数据中心站，完成同步处理和分析，构成广域测量系统。

随着WAMS的不断成熟与完善，将其应用于大电网的闭环控制必然是未来电力系统控制 发展的方向之一。WAMS可在一定的延时内获取机组间的相对功角和角速度，并向分散布置 的阻尼控制器提供全局信息，使得克服现有的阻尼控制器只能利用本地测量构成反馈控制这 一固有缺陷，并有效抑制本地和区间两种模式的低频振荡成为可能。

然而，广域信息在通信网络中传输和处理，存在较大的时滞。时滞是导致系统控制律失 效、运行状况恶化和系统失稳的一种重要诱因[Wu H X,Tsakalis K S,Heydt G T.Evaluation of  time delay effects to wide-area power system stabilizer design[J].IEEE Transactions on Power  Systems,2004,19(4):1935-1941.]，利用广域信息进行电力系统的闭环控制时，必须计及时滞 的影响。

在考虑时滞影响后的电力系统线性化微分-代数方程对应的特征方程中，标量时滞被转化 为指数项，因而，特征方程是一个超越方程，为此，在求取电力系统的时滞稳定裕度时，通 常采用函数变换的方法,如Rekasius变换[Rifat S,Nejat O.A novel stability study on multiple  time-delay systems(MTDS)using the root clustering paradigm[C].Proceedings of the American  Control Conference,Boston,MA,2004,5422-5427.]、Lambert-W函数[Yi S,Nelson P W,Ulsoy A  G.Time-delay systems:analysis and control suing the Lambert W function[M].Singapore:World  Scientific Publishing Company,2010]、SCF方法[Chen J,Gu G，Nett C N.A new method for  computing delay margins for stability of linear delay systems[J].Systems and Control Letters,1995, 26(2):107-117.]对超越项进行变换，避免直接求解特征方程的困难，这类方法存在变换复杂和 计算量大的不足。此外，基于线性矩阵不等式（Linear Matrix Inequlity，LMI）的时滞依赖稳 定的充分性判据[俞立.鲁棒控制——线性矩阵不等式出力方法[M].北京:清华大学出版社, 2002.]，也被广泛用于时滞稳定性分析和时滞上限求解[江全元,张鹏翔,曹一家.计及反馈信 号时滞影响的广域FACTS阻尼控制[J].中国电机工程学报,2006,26(7):82-88.Jiang  Quanyuan,Zhang Pengxiang,Cao Yijia.Wide-area FACTS damping control in consideration of  feedback signal's time delays[J].Proceedings of the CSEE,2006,26(7):82-88.]，此类方法通常与 系统模型降阶相结合以降低计算量，但仍存在固有的保守性缺点。

在广域阻尼控制研究中，Pade近似是一种常用的时滞环节处理方法。通过有理多项式来 逼近时滞环节，进而可以方便地利用经典和现代控制理论设计广域阻尼控制律[Wu H X,Ni H, Heyde G T.The impact of time delay on control design in power systems[C].Proceedings of IEEE  Power Engineering Society Winter Meeting,New York,USA,2002,1511-1516.王成山,石颉.考 虑时间延迟影响的电力系统稳定器设计[J].中国电机工程学报,2007,27(10):1-6.Wang  Chengshan,Shi Jie.PSS designing with consideration of time delay impact[J].Proceedings of the  CSEE,2007,27(10):1-6.石颉,王成山.考虑广域信息时延影响的H∞阻尼控制器[J].中国电 机工程学报,2008,28(1):30-34.Shi Jie,Wang Chengshan.Design of H∞controller for damping  interarea oscillations with consideration of the delay of remote signal[J].Proceedings of the CSEE, 2008,28(1):30-34.胡志坚,赵义术.计及广域测量系统时滞的互联电力系统鲁棒稳定控制[J]. 中国电机工程学报,2010,30(19):37-43.Hu Zhijian,Zhao Yishu.Robust stability control of  power systems based on WAMS with signal transmission delays[J].Proceedings of the CSEE,2010, 30(19):37-43.]，并通过时域仿真验证控制器的有效性。

特征值分析方法是小扰动稳定性分析的基本而最有效的方法。文献[贾宏杰,陈建华,余 晓丹.时滞环节对电力系统小扰动稳定性的影响[J].电力系统自动化,2006,30(10):5-8,17.

Jia Hongjie,Chen Jianhua,Yu Xiaodan.Impact of time delay on power system small signal  stability[J].Automation of Electric Power Systems,2006,20(5):5-8,17.]和[袁野,程林,孙元章, 等.广域阻尼控制的时滞影响分析及其时滞补偿设计[J].电力系统自动化,2006,30(14):6-9. Yuan Ye,Cheng Lin,Sun Yuanzhang,et al.Effect of delay input on wide-area damping control and  design of compensation[J].Automation of Electric Power Systems,2006,30(14):6-9.]分别通过计 算单机无穷大系统和等值两机系统的特征值来分析时滞对电力系统小扰动稳定性的影响，但 未给出详细的特征值计算方法。文献[贾宏杰,余晓丹.2种实际约束下的电力系统时滞稳定裕 度[J].电力系统自动化,2008,32(9):7-10,19.Jia Hongjie,Yu Xiaodan.Method of determining  power system delay margins with considering two practical constraints[J].Automation of Electric  Power Systems,2008,32(9):7-10,19.]提出了一种直接、有效搜索复平面上特定边界上（如虚轴、 特征值实部或阻尼比等于给定常数）时滞系统的关键特征值的方法。该方法能够得到的精确 的特征值，但搜索过程的计算量较大。

特征值分析方法已经形成了比较成熟和完善的理论，并在电力工业的实践中获得广泛应 用。如果能够提出时滞系统的特征值计算方法，进而沿用经典的特征值分析的思路和理论框 架来分析时滞电力系统的小扰动稳定性，无论对于完善和丰富基于特征值的小扰动稳定性分 析理论，还是促进广域阻尼控制的工程应用，都将具有重要的意义和价值。基于这种思想， 本发明提出了一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，通过建模和 嵌入时滞环节的状态空间表达，进而可以直接利用常规或稀疏特征值方法求得系统的部分特 征值，进而判别系统的时滞稳定性。针对四机两区算例系统，通过与离散化特征值求解方法 [Engelborghs K,Roose D.On stability of LMS methods and characteristic roots of delay  differential equations[J].SIAM Journal on Numerical Analysis,,2003,40(2):629-650.]计算结果 的对比，验证了本发明方法的正确性和有效性。

广域测量信号在传输和处理过程中产生的时滞，使电力系统成为一个时滞系统。时滞电 力系统各部分之间的连接关系如图2所示。

无时滞电力系统模型

设无广域阻尼控制器时描述电力系统的微分-代数方程组为：

$\left(\begin{array}{c}\stackrel{·}{x}=f(x,y)\\ 0=g(x,y)\end{array}\right)---\left(1\right)$

式中，f为描述元件动态的微分方程，g为网络方程，x为系统状态变量，y为系统代数 变量（节点电压），为系统状态变量的微分。

无时滞电力系统的输出为u，其通过广域反馈作为阻尼控制器的输入。设f_f为描述u与 (x,y)之间联系的函数，u＝f_f(x,y)。y_c为广域阻尼控制器的输出，并作为无时滞电力系统的控 制输入。在稳态运行点(x₀,y₀)对式(1)和u进行线性化，可得：

$\left(\begin{array}{c}\Delta \stackrel{·}{x}=\mathrm{A\Delta x}+\mathrm{B\Delta y}+\mathrm{E\Delta }{y}_c\\ 0=\mathrm{C\Delta x}+\mathrm{D\Delta y}\\ \mathrm{\Delta u}=K_1\mathrm{\Delta x}+K_2\mathrm{\Delta y}\end{array}\right)---\left(2\right)$

式中，A、B、C、D分别为为微分方程f和网络方程g相对于系统状态变量x和代数变量y 的偏导数，即$A=\frac{\partial f}{\partial x},$$B=\frac{\partial f}{\partial y},$$C=\frac{\partial g}{\partial x},$$D=\frac{\partial g}{\partial y},$$K_1=\frac{\partial f_f}{\partial x},$$K_2=\frac{\partial f_f}{\partial y},.$E是系统的输入矩阵，E 中的非零元素表征了附加广域阻尼控制器与被附加控制设备（如发电机励磁调节器、FACTS 设备等）之间的连接关系。E中非零元素所在的行，对应着控制设备中放大环节的输出变量 在无时滞电力系统状态变量x中的位置，非零元素值为放大环节的放大倍数与时间常数的比 值。

广域阻尼控制器状态空间表达

广域阻尼控制器的动态及输出可由如下微分-代数方程组表示：

$\left(\begin{array}{c}{\stackrel{·}{x}}_c=f_c(x_c,y_d)\\ y_c=g_c(x_c,y_d)\end{array}\right)---\left(3\right)$

式中，f_c为微分方程，g_c为代数方程，x_c为广域阻尼控制器中的状态变量，y_d为广域阻尼 控制器的输入。y_d=ue^-τt，u为广域阻尼控制器的输入，τ＝[τ₁,…,τ_i,…,τ_m]^T为所有广域阻尼控 制器的时滞形成的向量，τ_i>0为第i个时滞环节的时滞常数，i=1,2,…,m，m为正整数，表示 系统中时滞环节的总个数。

方程(3)对应的线性化方程为：

$\left(\begin{array}{c}\Delta {\stackrel{·}{x}}_c=A_c\Delta x_c+B_c\Delta y_d\\ \Delta y_c=C_c\Delta x_c+D_c\Delta y_d\end{array}\right)---\left(4\right)$

式中，A_c、B_c、C_c、D_c分别表示微分方程f_c和代数方程g_c相对于状态变量x_c和代数变量y_d的偏导数，即$A_c=\frac{\partial f_c}{\partial x_c},$$B_c=\frac{\partial f_c}{\partial y_d},$$C_c=\frac{{\partial g}_c}{{\partial x}_c},$$D_c=\frac{{\partial g}_c}{{\partial y}_c}.$

时滞电力系统模型

将描述系统动态元件的微分方程f和描述广域阻尼控制器动态的微分方程f_c写在一 起，形成矩阵方程向量f′，即相应地，将二者中的状态变量x和x_c写在一起， 则形成状态变量向量x′，即于是，考虑广域反馈时滞τ后，电力系统可由如下时 滞微分-代数方程描述：

$\left(\begin{array}{c}{\stackrel{·}{x}}^{\prime }=f^{\prime }(x^{\prime }y,x_{\tau 1}^{\prime },y_{\tau 1},x_{\tau 2}^{\prime },y_{\tau 2},...,x_{\mathrm{\tau m}}^{\prime },y_{\mathrm{\tau m}})\\ 0=g(x^{\prime },y)\\ 0=g\left({x_{\tau 1}^{\prime },y}_{\tau 1}\right)\\ 0=g(x_{\tau 2}^{\prime },y_{\tau 2})\\ .\\ .\\ .\\ 0=g(x_{\mathrm{\tau m}}^{\prime },y_{\mathrm{\tau m}})\end{array}\right)---\left(5\right)$

式中，[x′_τi,y_τi]＝[x′(t-τ_i),y(t-τ_i)]为时滞状态变量和代数变量，i为正整数。

在稳态运行点(x′₀,y₀)处对(5)进行线性化，可得全系统的线性化模型：

$\left(\begin{array}{c}\Delta {\stackrel{·}{x}}^{\prime }=A_0^{\prime }\Delta x^{\prime }{B}_0^{\prime }\mathrm{\Delta y}+\underset{i=1}{\overset{m}{\Sigma }}(A_{\mathrm{\tau i}}^{\prime }\Delta x_{\mathrm{\tau i}}^{\prime }+B_{\mathrm{\tau i}}^{\prime }\Delta y_{\mathrm{\tau i}})\\ 0=C_0^{\prime }\Delta x^{\prime }+D_0^{\prime }\mathrm{\Delta y}\\ 0=C_{\tau 1}^{\prime }\Delta x_{\tau 1}^{\prime }+D_{\tau 1}^{\prime }\Delta y_{\tau 1}\\ .\\ .\\ .\\ 0=C_{\mathrm{\tau m}}^{\prime }\Delta x_{\mathrm{\tau m}}^{\prime }+D_{\mathrm{\tau m}}^{\prime }\Delta y_{\mathrm{\tau m}}\end{array}\right)---\left(6\right)$

式中：Α₀′、B′₀、C′₀、D′₀分别表示微分方程f′和代数方程g对状态变量x′和代数变量y 的偏导数，Α′_τi、B′_τi、C′_τi、D′_τi分别表示微分方程f′和 代数方程g对时滞状态状态变量x′_τi和代数变量y_τi的偏导数

时滞电力系统特征值求解问题

当D₀′和D′_τi(i=1,2,…,m)非奇异时，消去代数变量，方程(6)可简化为：

$\Delta {\stackrel{·}{x}}^{\prime }={\tilde{A}}_0^{\prime }\Delta x^{\prime }+\underset{i=1}{\overset{m}{\Sigma }}{\tilde{A}}_{\mathrm{\tau i}}^{\prime }\Delta x_{\mathrm{\tau i}}^{\prime }---\left(7\right)$

式中，

$\left(\begin{array}{c}{\tilde{A}}_0^{\prime }=A_0^{\prime }-B_0^{\prime }{\left(D_0^{\prime }\right)}^{-1}{C}_0^{\prime }\\ {\tilde{A}}_{\mathrm{\tau i}}^{\prime }=A_{\mathrm{\tau i}}^{\prime }-B_{\mathrm{\tau i}}^{\prime }{\left(D_{\mathrm{\tau i}}^{\prime }\right)}^{-1}{C}_{\mathrm{\tau i}}^{\prime }\end{array}\right)---\left(8\right)$

式(7)表示的线性化系统的特征方程为：

$\det (\mathrm{\lambda I}-{\tilde{A}}_0^{\prime }-\underset{i=1}{\overset{m}{\Sigma }}{\tilde{A}}_{\mathrm{\tau i}}^{\prime }{e}^{-\lambda {\tau }_i})=0---\left(9\right)$

式中，λ为系统的特征值。

时滞动力系统的稳定性理论指出[廖晓昕.动力系统的稳定性理论和应用[M].北京:国 防工业出版社,2000.]，如果线性化系统(7)的全部特征值都具有负实部，则时滞系统(5)在稳态 运行点(x′₀,y₀)处是小扰动稳定的；反之，若至少存在一个具有正实部的特征值，则系统在该 点处是小扰动不稳定的。

然而，由于时滞及指数项的存在，线性化系统的特征方程(9)为超越方程，其有无穷 多个解。因此，通过直接求解该方程得到系统的关键特征值并判别系统的时滞稳定性变得非 常困难。

中国专利201010123345.8，提出了一种直接、有效搜索复平面上特定边界上（如虚轴、 特征值实部或阻尼比等于给定常数）时滞系统的关键特征值的方法。该方法能够得到的精确 的特征值，但搜索过程的计算量较大。

中国专利200810151217.7、200910070255.4、200910070254.X均是基于线性矩阵不等式 （Linear Matrix Inequlity，LMI）方法来判别时滞系统的小干扰稳定性。这类方法在判别系统 的稳定性时存在固有的保守性。

发明内容

本发明的目的就是为了解决上述问题，提供一种基于Pade近似的时滞电力系统特征值计 算与稳定性判别方法，它不需要迭代和搜索就可以准确地求解与时滞系统中动态元件相关的 部分特征值和特征向量，利用计算得到的特征值可以直接判断系统的小干扰稳定性，具有计 算量小，计算时间少，不存在任何保守性的优点。

为了实现上述目的，本发明采用如下技术方案：

一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，具体步骤如下：

步骤一：指定广域反馈时滞τ_i；

步骤二：进行Pade近似，得到时滞环节的近似有理多项式；

步骤三：将有理多项式转化为状态空间表达式；

步骤四：进行平衡化处理；

步骤五：将时滞环节与无时滞电力系统、广域阻尼控制器进行连接，得到时滞电力系统 的线性化模型；

步骤六：利用QR法或稀疏特征值算法求解系统的特征值；

步骤七：直接判别系统的时滞稳定性。

所述步骤二中，在拉普拉斯域中，第i个时滞环节可表示为t是时间，s是频率。 Pade近似是一种利用[l,k]阶有理多项式逼近的方法：

$e^{-{\tau }_{i}s}\approx R\left(s\right)=\frac{N_l\left(s\right)}{N_k\left(s\right)}=\frac{a_0+a_1{\tau }_{i}s+...+a_l{\left({\tau }_{i}s\right)}^l}{b_0+b_1{\tau }_{i}s+...+b_k{\left({\tau }_{i}s\right)}^k}---\left(10\right)$

式中，l和k为正整数。系数a_i和b_i为实数，可以由下式求出：

$a_j={(-1)}^j\frac{(l+k-j)!l!}{j!(l-j)!}---\left(11\right)$

$b_j=\frac{(l+k-j)!k!}{j!(k-j)!}---\left(12\right)$

所述步骤二在Pade近似中，阶数l和k越大，有理多项式R(s)越接近于通常情况 下，取l=k；时滞越小，频域中R(s)与相位一致的区间越大，即频带越宽。

所述步骤三中，将式(10)的分子和分母同时除以进而利用传递函数实现问题的一 般性方法[刘豹.现代控制理论[M].北京:机械工业出版社,2000.]，将基于Pade近似得到的 时滞环节的有理多项式逼近R(s)转换为等价的状态空间表达，对于第i个时滞环节，以变量 增量形式表示的能控标准型为：

$\left(\begin{array}{c}\Delta {\stackrel{·}{x}}_{\mathrm{di}}={\tilde{A}}_{\mathrm{di}}\Delta x_{\mathrm{di}}+{\tilde{B}}_d\Delta u_i\\ \Delta y_{\mathrm{di}}={\tilde{C}}_{\mathrm{di}}\Delta x_{\mathrm{di}}+{\tilde{D}}_{\mathrm{di}}\Delta u_i\end{array}\right)---\left(13\right)$

式（13）中，系数矩阵的具体表达式为：

k为Pade近似阶数。

在较小时滞τ_i和较高阶数的情况下，由Pade近似得到的传递函数R(s)中的系数之间，以 及由传递函数转换得到的状态空间表达式的系数矩阵A_di和C_di的非零元素之间，在数量级上 相差很大，采用线性变换使得时滞环节的状态空间表达式的系数矩阵更加平衡：

$A_{\mathrm{di}}=T^{-1}{\tilde{A}}_{\mathrm{di}}T,B_{\mathrm{di}}=T^{-1}{\tilde{B}}_{\mathrm{di}},C_{\mathrm{di}}={\tilde{C}}_{\mathrm{di}}T,D_{\mathrm{di}}={\tilde{D}}_{\mathrm{di}}---\left(15\right)$

式中，T为对角变换矩阵；由此可知，平衡化处理后的系数矩阵仍然符合能控标准型；

所述步骤五：将步骤三中得到的时滞环节的状态空间表达式，与无时滞电力系统的线性 化模型、广域阻尼控制器的线性化模型相连接，从而建立包含时滞环节的电力系统小扰动稳 定性分析的线性化模型：

$\left(\begin{array}{c}\Delta {\stackrel{·}{x}}^{\prime \prime }=A^{\prime \prime }\Delta x^{\prime \prime }+B^{\prime \prime }\mathrm{\Delta y}\\ 0=C^{\prime \prime }\Delta x^{\prime \prime }+D^{\prime \prime }\mathrm{\Delta y}\end{array}\right)---\left(16\right)$

式中，x″为系统状态变量，y为系统代数变量（节点电压），A"、B"、C"、D"为系数矩 阵；

假设第i个时滞环节的状态变量x_d和与之相连的广域阻尼控制器的状态变量x_c排列在无 时滞电力系统状态变量x之后，即$\Delta x^{\prime \prime }={\left(\begin{array}{ccccccccc}{\mathrm{\Delta x}}^T& .& .& .& {\mathrm{\Delta x}}_{\mathrm{di}}^T& {\mathrm{\Delta x}}_{\mathrm{ci}}^T& .& .& .\end{array}\right)}^T;$设无时滞电力系统的系数矩 阵为A，B，C，D；广域阻尼控制器的系数矩阵为A_c，B_c，C_c，D_c，广域反馈信号的系数矩 阵分别为K_1i和K_2i；于是，式(16)中系数矩阵A"、B"、C"、D"具体可表示为：

式(17)中，系数矩阵A_di、B_di、C_di、D_di满足能控标准型；A为分块对角阵，B、C为分块 稀疏矩阵，D为2×2分块稀疏矩阵[杜正春,刘伟,方万良,等.小干扰稳定性分析中一种关 键特征值计算的稀疏实现[J].中国电机工程学报,2005,25(2):17-21.Du Zhengchun,Liu Wei, Fang Wanliang,et al.A sparse method for the calculation of critical eigenvalue in small signal  stability analysis[J].Proceedings of the CSEE,2005,25(2):17-21.]；E为稀疏矩阵，Aci、Bci、Cci、 Dci与广域阻尼控制器的具体结构有关；

所述式(17)中各系数矩阵具有如下特点：

(1)若将广域反馈信号中涉及的动态元件、时滞环节、广域阻尼控制器状态变量的系数 矩阵作为一个子块，则A"仍为分块对角阵；B"仅在时滞环节和广域阻尼控制器状态变量所在 行上增加了零矩阵或少量几列非零元素，B"仍为稀疏阵；

(2)C"仅在时滞环节和广域阻尼控制器状态变量所在列上增加了零矩阵，D"与无时滞电力 系统相应的系数矩阵完全一样。

综上可知，包含时滞环节的电力系统的线性化方程系数矩阵A"、B"、C"、D"，与无时 滞电力系统在线性化方程系数矩阵A、B、C、D，具有完全相同的稀疏结构。

所述步骤六：当D"非奇异时，消去代数变量y，式(16)可简化为：

$\Delta {\stackrel{·}{x}}^{\prime \prime }={\tilde{A}}^{\prime \prime }\Delta x^{\prime \prime }---\left(18\right)$

式中，为基于Pade近似得到的时滞系统的状态矩阵。

通过计算的特征值，就可以得到时滞电力系统的部分特征值。此外，由于包含时滞环 节的电力系统和常规无时滞电力系统的系数矩阵具有相同的稀疏特性，在计算大规模无时滞 电力系统部分关键特征值时所使用的稀疏处理技术[杜正春,刘伟,方万良,等.小干扰稳定 性分析中一种关键特征值计算的稀疏实现[J].中国电机工程学报,2005,25(2):17-21.Du  Zhengchun,Liu Wei,Fang Wanliang,et al.A sparse method for the calculation of critical  eigenvalue in small signal stability analysis[J].Proceedings of the CSEE,2005,25(2):17-21.]和计 算方法[杜正春,刘伟,方万良,等.基于Jacobi-Davidson方法的小干扰稳定性分析中的关键 特征值计算[J].中国电机工程学报,2005,25(14):19-24.Du Zhengchun,Liu Wei,Fang  Wanliang,et al.The application of the Jacobi-Davidson method to the calculation of critical  eigenvalues in the small signal stability analysis[J].Proceedings of the CSEE,2005,25(14): 19-24.]，仍然适用于计算时滞电力系统的部分特征值。

本发明的有益效果：利用Pade近似将广域反馈时滞逼近为一个有理多项式，通过与无时 滞电力系统和广域阻尼控制器之间的连接，建立时滞电力系统的的线性化模型，最后根据系 统状态矩阵直接求出时滞系统的部分特征根；该方法能够较准确地求解与时滞系统中动态元 件相关的部分特征值和特征向量，正确求解与时滞环节相关的特征值个数和计算精度，与有 理多项式的阶数有关。

利用Pade近似有理多项式来逼近时滞环节，进而计算得到含有广域通信时滞电力系统的 部分特征值，从而直接判断电力系统的小干扰稳定性。

附图说明

图1为本发明的整体流程图；

图2为时滞电力系统示意图；

图3为两区四机系统图；

图4为实部大于-50的部分特征值；

图5为实部大于-10的部分特征值。

其中，1.无时滞电力系统，2.e^-sτ，3.广域阻尼控制器。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

如图1所示，一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，具体步 骤如下：

步骤一：指定广域反馈时滞τ_i和Pade近似阶数k；

步骤二：进行Pade近似，得到时滞环节的近似有理多项式；

步骤三：将有理多项式转化为状态空间表达式；

步骤四：进行平衡化处理；

步骤五：将时滞环节与无时滞电力系统、广域阻尼控制器进行连接，得到时滞电力系统 的线性化模型；

步骤六：利用QR法或稀疏特征值算法求解系统的特征值；

步骤七：直接判别系统的时滞稳定性。

所述步骤二中，在拉普拉斯域中，第i个时滞环节可表示为Pade近似是一种利 用[l,k]阶有理多项式逼近的方法：

$e^{-{\tau }_{i}s}\approx R\left(s\right)=\frac{N_l\left(s\right)}{N_k\left(s\right)}=\frac{a_0+a_1{\tau }_{i}s+...+a_l{\left({\tau }_{i}s\right)}^l}{b_0+b_1{\tau }_{i}s+...+b_k{\left({\tau }_{i}s\right)}^k}---\left(19\right)$

式中，系数ai和bi可以由下式求出：

$a_j={(-1)}^j\frac{(l+k-j)!l!}{j!(l-j)!}---\left(20\right)$

$b_j=\frac{(l+k-j)!k!}{j!(k-j)!}---\left(21\right)$

所述步骤二在Pade近似中，阶数l和k越大，有理多项式R(s)越接近于通常情况 下，取l=k；时滞越小，频域中R(s)与相位一致的区间越大，即频带越宽。

所述步骤三中，将式(10)的分子和分母同时除以进而利用传递函数实现问题的一 般性方法[刘豹.现代控制理论[M].北京:机械工业出版社,2000.]，将基于Pade近似得到的 时滞环节的有理多项式逼近R(s)转换为等价的状态空间表达，对于第i个时滞环节，以变量 增量形式表示的能控标准型为：

$\left(\begin{array}{c}\Delta {\stackrel{·}{x}}_{\mathrm{di}}={\tilde{A}}_{\mathrm{di}}\Delta x_{\mathrm{di}}+{\tilde{B}}_d\Delta u_i\\ \Delta y_{\mathrm{di}}={\tilde{C}}_{\mathrm{di}}\Delta x_{\mathrm{di}}+{\tilde{D}}_{\mathrm{di}}\Delta u_i\end{array}\right)---\left(22\right)$

式（13）中，

在较小时滞τ_i和较高阶数的情况下，由Pade近似得到的传递函数R(s)中的系数之间，以 及由传递函数转换得到的状态空间表达式的系数矩阵A_di和C_di的非零元素之间，在数量级上 相差很大，采用线性变换使得时滞环节的状态空间表达式的系数矩阵更加平衡：

$A_{\mathrm{di}}=T^{-1}{\tilde{A}}_{\mathrm{di}}T,B_{\mathrm{di}}=T^{-1}{\tilde{B}}_{\mathrm{di}},C_{\mathrm{di}}={\tilde{C}}_{\mathrm{di}}T,D_{\mathrm{di}}={\tilde{D}}_{\mathrm{di}}---\left(24\right)$

式中，T为对角变换矩阵；由此可知，平衡化处理后的系数矩阵仍然符合能控标准型；

所述步骤五：将步骤三中得到的时滞环节的状态空间表达式，与无时滞电力系统的线性 化模型、广域阻尼控制器的线性化模型相连接，从而建立包含时滞环节的电力系统小扰动稳 定性分析的线性化模型：

$\left(\begin{array}{c}\Delta {\stackrel{·}{x}}^{\prime \prime }=A^{\prime \prime }\Delta x^{\prime \prime }+B^{\prime \prime }\mathrm{\Delta y}\\ 0=C^{\prime \prime }\Delta x^{\prime \prime }+D^{\prime \prime }\mathrm{\Delta y}\end{array}\right)---\left(25\right)$

式中，x″为系统状态变量，y为系统代数变量（节点电压），A"、B"、C"、D"为系数矩 阵；

假设第i个时滞环节的状态变量x_d和与之相连的广域阻尼控制器的状态变量x_c排列在无 时滞电力系统状态变量x之后，即$\Delta x^{\prime \prime }={\left(\begin{array}{ccccccccc}{\mathrm{\Delta x}}^T& .& .& .& {\mathrm{\Delta x}}_{\mathrm{di}}^T& {\mathrm{\Delta x}}_{\mathrm{ci}}^T& .& .& .\end{array}\right)}^T;$设无时滞电力系统的系数矩 阵为A，B，C，D；广域阻尼控制器的系数矩阵为A_c，B_c，C_c，D_c，广域反馈信号的系数矩 阵分别为K_1i和K_2i；于是，式(16)中系数矩阵A"、B"、C"、D"具体可表示为：

式(17)中，系数矩阵A_di、B_di、C_di、D_di满足能控标准型；A为分块对角阵，B、C为分块 稀疏矩阵，D为2×2分块稀疏矩阵[杜正春,刘伟,方万良,等.小干扰稳定性分析中一种关 键特征值计算的稀疏实现[J].中国电机工程学报,2005,25(2):17-21.Du Zhengchun,Liu Wei, Fang Wanliang,et al.A sparse method for the calculation of critical eigenvalue in small signal  stability analysis[J].Proceedings of the CSEE,2005,25(2):17-21.]；E为稀疏矩阵，Aci、Bci、Cci、 Dci与广域阻尼控制器的具体结构有关；

所述式(17)中各系数矩阵具有如下特点：

(1)若将广域反馈信号中涉及的动态元件、时滞环节、广域阻尼控制器状态变量的系数 矩阵作为一个子块，则A"仍为分块对角阵；B"仅在时滞环节和广域阻尼控制器状态变量所在 行上增加了零矩阵或少量几列非零元素，B"仍为稀疏阵；

(2)C"仅在时滞环节和广域阻尼控制器状态变量所在列上增加了零矩阵，D"与无时滞电力 系统相应的系数矩阵完全一样。

综上可知，包含时滞环节的电力系统的线性化方程系数矩阵A"、B"、C"、D"，与无时 滞电力系统在线性化方程系数矩阵A、B、C、D，具有完全相同的稀疏结构。

所述步骤六：当D"非奇异时，消去代数变量y，式(16)可简化为：

$\Delta {\stackrel{·}{x}}^{\prime \prime }={\tilde{A}}^{\prime \prime }\Delta x^{\prime \prime }---\left(27\right)$

式中，为基于Pade近似得到的时滞系统的状态矩阵。

通过计算的特征值，就可以得到时滞电力系统的部分特征值。此外，由于包含时滞环 节的电力系统和常规无时滞电力系统的系数矩阵具有相同的稀疏特性，在计算大规模无时滞 电力系统部分关键特征值时所使用的稀疏处理技术[杜正春,刘伟,方万良,等.小干扰稳定 性分析中一种关键特征值计算的稀疏实现[J].中国电机工程学报,2005,25(2):17-21.Du  Zhengchun,Liu Wei,Fang Wanliang,et al.A sparse method for the calculation of critical  eigenvalue in small signal stability analysis[J].Proceedings of the CSEE,2005,25(2):17-21.]和计 算方法[杜正春,刘伟,方万良,等.基于Jacobi-Davidson方法的小干扰稳定性分析中的关键 特征值计算[J].中国电机工程学报,2005,25(14):19-24.Du Zhengchun,Liu Wei,Fang  Wanliang,et al.The application of the Jacobi-Davidson method to the calculation of critical  eigenvalues in the small signal stability analysis[J].Proceedings of the CSEE,2005,25(14): 19-24.]，仍然适用于计算时滞电力系统的部分特征值。

广域测量信号在传输和处理过程中产生的时滞，使电力系统成为一个时滞系统。时滞电 力系统各部分之间的连接关系如图2所示，无时滞电力系统与e^-sτ连接，e^-sτ与广域阻尼控 制器连接。

如图3所示的两区四机电力系统为例，验证本发明提出的时滞电力系统特征值计算方法 的正确性和有效性,系统的参数详见[Kunder P.Power System Stabilityand Control[M].New  York:McGraw-Hill,1994.]。

在发电机G1和G3安装本地电力系统稳定器（Power System Stabilizer,PSS）的基础上， 考虑在G1上装设以G1和G3相对转速偏差Δω13为反馈信号的广域PSS，进一步提高系统 的阻尼水平。

本发明以基于离散化的时滞系统特征值求解方法计算得到的部分特征值作为精确解，并 作为验证本发明提出的基于Pade近似的特征值计算方法有效性和精确性的基准。

本发明借助Matlab工具箱DDE-BIFTOOL中p_stabil函数计算得到的算例系统的部分特 征值的精确解，其中，牛顿法迭代的收敛精度取为1e-8。

针对算例系统，在不同的反馈时滞τ和不同的Pade近似有理多项式的阶数k下，利用本 发明提出的基于Pade近似的特征值计算方法进行了大量的计算，并与基于离散化的特征值计 算方法的结果进行了对比和分析，下面仅以部分计算结果来说明本发明方法的正确性和有效 性。

当反馈时滞τ=0.30s时，利用基于离散化的特征值计算方法、基于5阶和20阶Pade近似 的特征值计算方法，得到的实部大于-50的部分特征值，如图4所示。图4右下角矩形框内的、 实部大于-10的部分特征值的详细分布如图5所示，由图5可知：

(1)当k=5和20时，利用基于Pade近似的特征值计算方法，能够较为准确地计算得到 实部大于-10的复平面内的部分特征值，通过进一步分析可知，其与描述系统动态元件的状态 变量对应。

(2)对于k=5，在实部小于-10的复平面内，基于Pade近似的特征值计算方法，会得到若 干个错误的特征值，如：-15.06761008±23.47922927i、-21.53755703±15.83642620i、 -36.68588408、-49.61254513；此外，还有多个特征值未被计算出来，通过进一步分析可知， 遗漏和错误求解的特征值，主要与有理多项式的各阶变量对应。

(3)只有当有理多项式的阶数显著增大到k=20时，才能保证本发明方法在实部大于-41 的复平面内不出现漏根、错根的情况，与有理多项式变量对应的部分特征值也能较为准确地 求得。当阶数增加到30时，在实部大于-50的复平面内，本发明与基于离散化的特征值求解 方法得到特征值完全相同。

表1τ=0.3s时，基于离散化和Pade近似计算得到的特征值

(4)在增加Pade近似有理多项式阶数以准确求取与有理多项式变量对应的部分特征值的 同时，与系统动态元件相对应的部分特征值的精度也显著地得到提高。如表1所示，例如， 当k=5时，本发明方法与基于离散化方法得到的特征值之间的最大绝对误差max(abs(ΔRe(λ), ΔIm(λ)))=1.042e-4，当k增加到7时，最大绝对误差相应地减小到4.673e-8，特征值计算精度 得到显著的提高。

(5)由计算得到的系统的特征值可知，当τ=0.3s时，系统是小干扰稳定的。

为了深入考察Pade近似有理多项式的阶数k对特征值计算精度的影响，在不同时滞值下， 基于不同阶数Pade近似的特征值计算方法得到的特征值，与基于离散化方法得到的特征值之 间的最大绝对误差如表2所示。不难看出，在给定精度要求下，随着时滞的增大，Pade近似 有理多项式阶数k须相应地增大。由表2中可知，对于4机2区域系统，τ在[0.05，0.5]范围 内变化时，6阶Pade近似就能使本发明方法计算得到的特征值达到1e-5的计算精度。

表2不同时滞和Pade近似阶数下，特征值的计算精度

当反馈时滞τ=0.30时，基于离散化的特征值计算方法得到的系统部分特征值对应的右特 征向量的准确值，如表3第2~4行所示。通过进一步的模态分析可得系统的三个机电振荡模 式及相应的模态，如表3第2~4行所示。其中，λ1表现为G1、G2相对于G3、G4的振荡； λ2表现为G1相对于G2、G3相对于G4的振荡，但对应G1、G2模态分量的模值较大；λ3 表现为G1相对于G2、G3相对于G4的振荡，但对应G3、G4模态分量的模值较大。

利用5阶和7阶Pade近似的计算得到的系统部分特征值和相应的右特征向量中与各发电 机转速相应的分量，如表3的第5~10行所示。通过与模态的精确值的对比，可知本发明可以 准确计算得到振荡模态的幅值和相位。k=5时，模态幅值的最大绝对误差为4e-8，相位的最 大绝对误差为1e-2°；k=7时，幅值的最大绝度误差为1e-8；模态相位的最大绝对误差为9e-6°。

表3τ=0.3s时，Pade近似计算得到的模态及其精确值

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限 制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付 出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。