CN102801158B

CN102801158B - Method for calculating time-lag electric power system eigenvalue and discriminating stability based on Pade approximation

Info

Publication number: CN102801158B
Application number: CN201210271783.8A
Authority: CN
Inventors: 牛新生; 叶华; 王春义; 贾善杰; 刘晓明
Original assignee: Economic and Technological Research Institute of State Grid Shandong Electric Power Co Ltd
Current assignee: Economic and Technological Research Institute of State Grid Shandong Electric Power Co Ltd
Priority date: 2012-07-31
Filing date: 2012-07-31
Publication date: 2014-07-09
Anticipated expiration: 2032-07-31
Also published as: CN102801158A

Abstract

The invention discloses a time-delay power system eigenvalue calculation and stability discrimination method based on Pade approximation. It utilizes Pade approximation to approximate the time-delay of wide-area feedback to a rational polynomial. The connection between the controllers is used to establish a linearized model of the time-delay power system. Finally, some characteristic roots of the time-delay system are directly obtained according to the system state matrix, and then the time-delay stability of the system is judged. The eigenvalue calculation results of the four-machine two-area example system, this method can accurately solve some eigenvalues and eigenvectors related to the dynamic components in the time-delay system, and can conveniently and accurately judge the time-delay stability of the system, and correctly solve The number of eigenvalues and calculation precision related to the delay link are related to the order of the rational polynomial. In addition, this method has the advantages of small calculation amount and short calculation time.

Description

Eigenvalue Calculation and Stability Discrimination Method of Time-delay Power System Based on Pade Approximation

技术领域 technical field

本发明涉及一种时滞电力系统特征值计算与稳定性判别方法，尤其涉及一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法。The invention relates to a time-delay power system eigenvalue calculation and stability discrimination method, in particular to a time-delay power system eigenvalue calculation and stability discrimination method based on Pade approximation.

背景技术 Background technique

随着电力系统规模的扩大，区域间的低频振荡正成为限制电网传输能力的瓶颈，而现有的阻尼控制器（主要是按照单机无穷大系统设计的电力系统稳定器）并不能很好地解决这一问题。根本原因在于：(1)不能直接利用相对功角和相对角速度构成闭环控制。虽然采用相对功角和相对角速度来实现阻尼控制是最直接和有效的，但长期以来缺少必要的测量方法，只能采用其他间接变量来代替，导致控制效果不佳。(2)限于本地局部信息。采用本地测量构成反馈控制，不能很好地反映区间振荡模式，导致控制系统虽然能阻尼本地振荡模式，但难以有效地抑制区间振荡模式。With the expansion of the scale of the power system, the low-frequency oscillation between regions is becoming a bottleneck limiting the transmission capacity of the power grid, and the existing damping controllers (mainly power system stabilizers designed according to the single-machine infinite system) cannot solve this problem well. a question. The fundamental reasons are: (1) The relative power angle and relative angular velocity cannot be directly used to form a closed-loop control. Although using relative power angle and relative angular velocity to realize damping control is the most direct and effective, there has been a lack of necessary measurement methods for a long time, and other indirect variables can only be used instead, resulting in poor control effect. (2) Limited to local local information. Using local measurement to form feedback control cannot reflect the interval oscillation mode well, so that although the control system can damp the local oscillation mode, it is difficult to effectively suppress the interval oscillation mode.

随着信息处理和通信技术的迅猛发展，同步相量技术和广域测量系统（Wide AreaMeasurement System,WAMS），给电力系统的监测、分析和控制提供了新的手段，为互联电网阻尼控制带来了新的契机。同步相量测量单元可同步采集表征电网运行状态的几乎所有的变量，最为关键的是它能测量发电机的内电势、转子角、角速度，母线电压相位等与低频振荡密切相关的量。在高速通信网络（如电力数据宽带网）的支持下，各相量测量单元采集的带时标的数据能以较小的延时传递到数据中心站，完成同步处理和分析，构成广域测量系统。With the rapid development of information processing and communication technology, synchrophasor technology and Wide Area Measurement System (Wide Area Measurement System, WAMS) provide new means for the monitoring, analysis and control of the power system, and bring great benefits to the damping control of the interconnected grid. a new opportunity. The synchrophasor measurement unit can synchronously collect almost all variables that characterize the operating state of the power grid. The most important thing is that it can measure the internal potential of the generator, rotor angle, angular velocity, bus voltage phase and other quantities closely related to low-frequency oscillation. With the support of high-speed communication network (such as power data broadband network), the time-scaled data collected by each phasor measurement unit can be transmitted to the data center station with a small delay, and complete synchronous processing and analysis to form a wide-area measurement system .

随着WAMS的不断成熟与完善，将其应用于大电网的闭环控制必然是未来电力系统控制发展的方向之一。WAMS可在一定的延时内获取机组间的相对功角和角速度，并向分散布置的阻尼控制器提供全局信息，使得克服现有的阻尼控制器只能利用本地测量构成反馈控制这一固有缺陷，并有效抑制本地和区间两种模式的低频振荡成为可能。With the continuous maturity and improvement of WAMS, applying it to the closed-loop control of the large power grid must be one of the future development directions of power system control. WAMS can obtain the relative power angle and angular velocity between units within a certain delay, and provide global information to the distributed damping controller, so as to overcome the inherent defect that the existing damping controller can only use local measurement to form feedback control , and it becomes possible to effectively suppress low-frequency oscillations in both local and interval modes.

然而，广域信息在通信网络中传输和处理，存在较大的时滞。时滞是导致系统控制律失效、运行状况恶化和系统失稳的一种重要诱因[Wu H X,Tsakalis K S,Heydt G T.Evaluation oftime delay effects to wide-area power system stabilizer design[J].IEEE Transactions on PowerSystems,2004,19(4):1935-1941.]，利用广域信息进行电力系统的闭环控制时，必须计及时滞的影响。However, there is a large time lag in the transmission and processing of wide-area information in the communication network. Time delay is an important cause of system control law failure, deterioration of operating conditions and system instability[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 PowerSystems,2004,19(4):1935-1941.], when using wide-area information for closed-loop control of power systems, the influence of time delay must be taken into account.

在考虑时滞影响后的电力系统线性化微分-代数方程对应的特征方程中，标量时滞被转化为指数项，因而，特征方程是一个超越方程，为此，在求取电力系统的时滞稳定裕度时，通常采用函数变换的方法,如Rekasius变换[Rifat S,Nejat O.A novel stability study on multipletime-delay systems(MTDS)using the root clustering paradigm[C].Proceedings of the AmericanControl Conference,Boston,MA,2004,5422-5427.]、Lambert-W函数[Yi S,Nelson P W,Ulsoy AG.Time-delay systems:analysis and control suing the Lambert W function[M].Singapore:WorldScientific Publishing Company,2010]、SCF方法[Chen J,Gu G，Nett C N.A new method forcomputing 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.JiangQuanyuan,Zhang Pengxiang,Cao Yijia.Wide-area FACTS damping control in consideration offeedback signal's time delays[J].Proceedings of the CSEE,2006,26(7):82-88.]，此类方法通常与系统模型降阶相结合以降低计算量，但仍存在固有的保守性缺点。In the characteristic equation corresponding to the linearized differential-algebraic equation of the power system after considering the influence of time delay, the scalar time delay is transformed into an exponential term, therefore, the characteristic equation is a transcendental equation. Therefore, when calculating the time delay of the power system For the stability margin, the method of function transformation is usually adopted, such as Rekasius transformation [Rifat S, Nejat O.A novel stability study on multipletime-delay systems (MTDS) using the root clustering paradigm [C]. Proceedings of the American Control Conference, Boston, MA ,2004,5422-5427.], Lambert-W function [Yi S, Nelson P W, Ulsoy AG. Time-delay systems: analysis and control suing the Lambert W function [M]. Singapore: World Scientific Publishing Company, 2010], The SCF method [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.] transforms the transcendental items, To avoid the difficulty of directly solving the characteristic equation, this type of method has the disadvantages of complex transformation and large amount of calculation. In addition, the adequacy criterion of time-delay-dependent stability based on Linear Matrix Inequality (LMI) [Yu Li. Robust Control - Linear Matrix Inequality Output Method [M]. Beijing: Tsinghua University Press, 2002. ], is also widely used in time-delay stability analysis and time-delay upper limit solution [Jiang Quanyuan, Zhang Pengxiang, Cao Jia. Wide-area FACTS damping control considering the influence of feedback signal time-delay[J]. Chinese Journal of Electrical Engineering, 2006,26 (7):82-88.JiangQuanyuan, Zhang Pengxiang,Cao Yijia.Wide-area FACTS damping control in consideration offeedback signal's time delays[J].Proceedings of the CSEE,2006,26(7):82-88.], Such methods are usually combined with system model reduction to reduce the amount of calculation, but still have inherent conservative shortcomings.

在广域阻尼控制研究中，Pade近似是一种常用的时滞环节处理方法。通过有理多项式来逼近时滞环节，进而可以方便地利用经典和现代控制理论设计广域阻尼控制律[Wu H X,Ni H,Heyde G T.The impact of time delay on control design in power systems[C].Proceedings of IEEEPower Engineering Society Winter Meeting,New York,USA,2002,1511-1516.王成山,石颉.考虑时间延迟影响的电力系统稳定器设计[J].中国电机工程学报,2007,27(10):1-6.WangChengshan,Shi Jie.PSS designing with consideration of time delay impact[J].Proceedings of theCSEE,2007,27(10):1-6.石颉,王成山.考虑广域信息时延影响的H∞阻尼控制器[J].中国电机工程学报,2008,28(1):30-34.Shi Jie,Wang Chengshan.Design of H∞controller for dampinginterarea 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 ofpower systems based on WAMS with signal transmission delays[J].Proceedings of the CSEE,2010,30(19):37-43.]，并通过时域仿真验证控制器的有效性。In the study of wide-area damping control, the Pade approximation is a commonly used method for dealing with time-delay links. The time-delay link is approximated by rational polynomials, and then the wide-area damping control law can be conveniently designed using classical and modern control theories[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. Wang Chengshan, Shi Jie. Power System Stabilizer Design Considering Time Delay Effect [J]. Chinese Journal of Electrical Engineering, 2007, 27( 10):1-6.WangChengshan, Shi Jie.PSS designing with consideration of time delay impact[J].Proceedings of the CSEE,2007,27(10):1-6.Shi Jie,Wang Chengshan.When considering wide-area information H∞ damping controller for damping interarea oscillations with consideration of the delay of remote signal [J].Proceedings of the CSEE,2008,28(1):30-34. Hu Zhijian, Zhao Yishu. Robust Stability Control of Interconnected Power System Considering Time Delay of Wide-area Measurement System[J]. Chinese Journal of Electrical Engineering ,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. ], and verify the effectiveness of the controller by time domain simulation.

特征值分析方法是小扰动稳定性分析的基本而最有效的方法。文献[贾宏杰,陈建华,余晓丹.时滞环节对电力系统小扰动稳定性的影响[J].电力系统自动化,2006,30(10):5-8,17.The eigenvalue analysis method is the basic and most effective method for small disturbance stability analysis. Literature[Jia Hongjie, Chen Jianhua, Yu Xiaodan. Influence of time-delay link on small disturbance stability of power system[J].Automation of Electric Power System,2006,30(10):5-8,17.

Jia Hongjie,Chen Jianhua,Yu Xiaodan.Impact of time delay on power system small signalstability[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 anddesign 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 determiningpower system delay margins with considering two practical constraints[J].Automation of ElectricPower Systems,2008,32(9):7-10,19.]提出了一种直接、有效搜索复平面上特定边界上（如虚轴、特征值实部或阻尼比等于给定常数）时滞系统的关键特征值的方法。该方法能够得到的精确的特征值，但搜索过程的计算量较大。Jia Hongjie, Chen Jianhua, Yu Xiaodan.Impact of time delay on power system small signalstability[J].Automation of Electric Power Systems,2006,20(5):5-8,17.] and [Yuan Ye, Cheng Lin, Sun Yuanzhang, et al. Time-delay influence analysis and time-delay compensation design of wide-area damping control[J]. Electric Power System Automation, 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.] By calculating the eigenvalues of the single-machine infinite system and the equivalent two-machine system respectively To analyze the influence of time delay on the small disturbance stability of power system, but the detailed eigenvalue calculation method is not given. Literature [Jia Hongjie, Yu Xiaodan. Time-delay stability margin of power system under two kinds of practical constraints [J]. Electric Power System Automation, 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 ElectricPower Systems,2008,32(9):7-10,19.] proposes a direct and effective search for specific boundaries on the complex plane (such as imaginary axis, feature method to value the real part or damping ratio equal to a given constant) for the key eigenvalues of a time-delay system. This method can obtain accurate eigenvalues, but the calculation of the search process is relatively large.

特征值分析方法已经形成了比较成熟和完善的理论，并在电力工业的实践中获得广泛应用。如果能够提出时滞系统的特征值计算方法，进而沿用经典的特征值分析的思路和理论框架来分析时滞电力系统的小扰动稳定性，无论对于完善和丰富基于特征值的小扰动稳定性分析理论，还是促进广域阻尼控制的工程应用，都将具有重要的意义和价值。基于这种思想，本发明提出了一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，通过建模和嵌入时滞环节的状态空间表达，进而可以直接利用常规或稀疏特征值方法求得系统的部分特征值，进而判别系统的时滞稳定性。针对四机两区算例系统，通过与离散化特征值求解方法[Engelborghs K,Roose D.On stability of LMS methods and characteristic roots of delaydifferential equations[J].SIAM Journal on Numerical Analysis,,2003,40(2):629-650.]计算结果的对比，验证了本发明方法的正确性和有效性。The eigenvalue analysis method has formed a relatively mature and perfect theory, and has been widely used in the practice of the electric power industry. If the eigenvalue calculation method of the time-delay system can be proposed, and then the classic eigenvalue analysis ideas and theoretical framework are used to analyze the small-disturbance stability of the time-delay power system, whether it is necessary to improve and enrich the small-disturbance stability analysis based on eigenvalues Both the theory and the engineering application to promote wide-area damping control will be of great significance and value. Based on this idea, the present invention proposes a Pade approximation-based method for calculating the eigenvalues of time-delay power systems and determining stability. By modeling and embedding the state space expression of time-delay links, conventional or sparse eigenvalues can be directly used The method obtains some eigenvalues of the system, and then judges the time-delay stability of the system. For the example system of four machines and two areas, by solving the method of discrete eigenvalues [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.] The comparison of calculation results has verified the correctness and effectiveness of the method of the present invention.

广域测量信号在传输和处理过程中产生的时滞，使电力系统成为一个时滞系统。时滞电力系统各部分之间的连接关系如图2所示。The time delay generated during the transmission and processing of the wide-area measurement signal makes the power system a time-delay system. The connection relationship between various parts of the time-delay power system is shown in Figure 2.

无时滞电力系统模型Delay-free power system model

设无广域阻尼控制器时描述电力系统的微分-代数方程组为：Assuming that there is no wide-area damping controller, the differential-algebraic equations describing the power system are:

$\{\begin{matrix} \overset{\cdot &Center Dot;}{x x} = = f f ((x x,, y the y)) \\ 00 = = g g ((x x,, y the y)) \end{matrix} - - - - - - ((11))$

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

为系统状态变量的微分。In the formula, f is the differential equation describing the component dynamics, g is the network equation, x is the system state variable, y is the system algebraic variable (node voltage),

is the differential of the system state variable.

无时滞电力系统的输出为u，其通过广域反馈作为阻尼控制器的输入。设f_f为描述u与(x,y)之间联系的函数，u＝f_f(x,y)。y_c为广域阻尼控制器的输出，并作为无时滞电力系统的控制输入。在稳态运行点(x₀,y₀)对式(1)和u进行线性化，可得：The output of the time-delay-free power system is u, which is used as the input of the damping controller through wide-area feedback. Let f _f be a function describing the relation between u and (x, y), u=f _f (x, y). y _c is the output of the wide-area damping controller and serves as the control input of the time-delay-free power system. Linearize equation (1) and u at the steady-state operating point (x ₀ , y ₀ ), and get:

$\{\begin{matrix} Δ Δ \overset{\cdot &Center Dot;}{x x} = = AΔx AΔx + + BΔy BΔy + + EΔ EΔ {y the y}_{c c} \\ 00 = = CΔx CΔx + + DΔy DΔy \\ Δu Δu = = {K K}_{11} Δx Δx + + {K K}_{22} Δy Δy \end{matrix} - - - - - - ((22))$

式中，A、B、C、D分别为为微分方程f和网络方程g相对于系统状态变量x和代数变量y的偏导数，即 $A = \frac{&PartialD; f}{&PartialD; x},$ $B = \frac{&PartialD; f}{&PartialD; y},$ $C = \frac{&PartialD; g}{&PartialD; x},$ $D = \frac{&PartialD; g}{&PartialD; y},$ $K_{1} = \frac{&PartialD; f_{f}}{&PartialD; x},$ $K_{2} = \frac{&PartialD; f_{f}}{&PartialD; y}, .$ E是系统的输入矩阵，E中的非零元素表征了附加广域阻尼控制器与被附加控制设备（如发电机励磁调节器、FACTS设备等）之间的连接关系。E中非零元素所在的行，对应着控制设备中放大环节的输出变量在无时滞电力系统状态变量x中的位置，非零元素值为放大环节的放大倍数与时间常数的比值。In the formula, A, B, C, and D are the partial derivatives of the differential equation f and the network equation g relative to the system state variable x and the algebraic variable y, namely $A = \frac{&PartialD; f}{&PartialD; x},$ $B = \frac{&PartialD; f}{&PartialD; the y},$ $C = \frac{&PartialD; g}{&PartialD; x},$ $D. = \frac{&PartialD; g}{&PartialD; the y},$ $K_{1} = \frac{&PartialD; f_{f}}{&PartialD; x},$ $K_{2} = \frac{&PartialD; f_{f}}{&PartialD; the y}, .$ E is the input matrix of the system, and the non-zero elements in E represent the connection relationship between the additional wide-area damping controller and the additional control equipment (such as generator excitation regulator, FACTS equipment, etc.). The row of the non-zero element in E corresponds to the position of the output variable of the amplification link in the control device in the state variable x of the non-delay power system, and the value of the non-zero element is the ratio of the amplification factor of the amplification link to the time constant.

广域阻尼控制器状态空间表达State-space expression of wide-area damping controller

广域阻尼控制器的动态及输出可由如下微分-代数方程组表示：The dynamics and output of the wide-area damping controller can be expressed by the following differential-algebraic equations:

$\{\begin{matrix} {\overset{\cdot &Center Dot;}{x x}}_{c c} = = {f f}_{c c} (({x x}_{c c},, {y the y}_{d d})) \\ {y the y}_{c c} = = {g g}_{c c} (({x x}_{c c},, {y the y}_{d d})) \end{matrix} - - - - - - ((33))$

式中，f_c为微分方程，g_c为代数方程，x_c为广域阻尼控制器中的状态变量，y_d为广域阻尼控制器的输入。y_d=ue^-τt，u为广域阻尼控制器的输入，τ＝[τ₁,…,τ_i,…,τ_m]^T为所有广域阻尼控制器的时滞形成的向量，τ_i>0为第i个时滞环节的时滞常数，i=1,2,…,m，m为正整数，表示系统中时滞环节的总个数。where f _c is a differential equation, g _c is an algebraic equation, x _c is the state variable in the wide-area damping controller, and y _d is the input of the wide-area damping controller. y _d =ue ^-τt , u is the input of the wide-area damping controller, τ=[τ ₁ ,…,τ _i ,…,τ _m ] ^T is the vector formed by time-delays of all wide-area damping controllers, τ _i >0 is the time-delay constant of the i-th time-delay link, i=1,2,...,m, m is a positive integer, indicating the total number of time-delay links in the system.

方程(3)对应的线性化方程为：The linearization equation corresponding to equation (3) is:

$\{\begin{matrix} Δ Δ {\overset{\cdot &Center Dot;}{x x}}_{c c} = = {A A}_{c c} Δ Δ {x x}_{c c} + + {B B}_{c c} Δ Δ {y the y}_{d d} \\ Δ Δ {y the y}_{c c} = = {C C}_{c c} Δ Δ {x x}_{c c} + + {D D.}_{c c} Δ Δ {y the y}_{d d} \end{matrix} - - - - - - ((44))$

式中，A_c、B_c、C_c、D_c分别表示微分方程f_c和代数方程g_c相对于状态变量x_c和代数变量y_d的偏导数，即 $A_{c} = \frac{&PartialD; f_{c}}{&PartialD; x_{c}},$ $B_{c} = \frac{&PartialD; f_{c}}{&PartialD; y_{d}},$ $C_{c} = \frac{{&PartialD; g}_{c}}{{&PartialD; x}_{c}},$ $D_{c} = \frac{{&PartialD; g}_{c}}{{&PartialD; y}_{c}} .$ In the formula, A _c , B _c , C _c , and D _c represent the partial derivatives of the differential equation f _c and the algebraic equation g _c with respect to the state variable x _c and the algebraic variable y _d respectively, namely $A_{c} = \frac{&PartialD; f_{c}}{&PartialD; x_{c}},$ $B_{c} = \frac{&PartialD; f_{c}}{&PartialD; {the y}_{d}},$ $C_{c} = \frac{{&PartialD; g}_{c}}{{&PartialD; x}_{c}},$ ${D.}_{c} = \frac{{&PartialD; g}_{c}}{{&PartialD; the y}_{c}} .$

时滞电力系统模型Time-delay power system model

将描述系统动态元件的微分方程f和描述广域阻尼控制器动态的微分方程f_c写在一起，形成矩阵方程向量f′，即

相应地，将二者中的状态变量x和x_c写在一起，则形成状态变量向量x′，即

于是，考虑广域反馈时滞τ后，电力系统可由如下时滞微分-代数方程描述：Write the differential equation f describing the dynamic components of the system and the differential equation f _c describing the dynamics of the wide-area damping controller together to form a matrix equation vector f′, namely

Correspondingly, write the state variables x and x _c in the two together to form the state variable vector x′, that is

Then, after considering the wide-area feedback delay τ, the power system can be described by the following delay differential-algebraic equation:

$\{\begin{matrix} {\overset{\cdot &Center Dot;}{x x}}^{' '} = = {f f}^{' '} (({x x}^{' '} y the y,, {x x}_{τ τ 11}^{' '},, {y the y}_{τ τ 11},, {x x}_{τ τ 22}^{' '},, {y the y}_{τ τ 22},, . . . . . .,, {x x}_{τm τm}^{' '},, {y the y}_{τm τm})) \\ 00 = = g g (({x x}^{' '},, y the y)) \\ 00 = = g g (({{x x}_{τ τ 11}^{' '},, y the y}_{τ τ 11})) \\ 00 = = g g (({x x}_{τ τ 22}^{' '},, {y the y}_{τ τ 22})) \\ . . \\ . . \\ . . \\ 00 = = g g (({x x}_{τm τm}^{' '},, {y the y}_{τm τm})) \end{matrix} - - - - - - ((55))$

式中，[x′_τi,y_τi]＝[x′(t-τ_i),y(t-τ_i)]为时滞状态变量和代数变量，i为正整数。In the formula, [x′ _τi ,y _τi ]=[x′(t-τ _i ), y(t-τ _i )] is the time-delay state variable and algebraic variable, and i is a positive integer.

在稳态运行点(x′₀,y₀)处对(5)进行线性化，可得全系统的线性化模型：Linearize (5) at the steady-state operating point (x′ ₀ ,y ₀ ), and the linearized model of the whole system can be obtained:

$\{\begin{matrix} Δ Δ {\overset{\cdot &Center Dot;}{x x}}^{' '} = = {A A}_{00}^{' '} Δ Δ {x x}^{' '} {B B}_{00}^{' '} Δy Δy + + {Σ Σ}_{i i = = 11}^{m m} (({A A}_{τi τi}^{' '} Δ Δ {x x}_{τi τi}^{' '} + + {B B}_{τi τi}^{' '} Δ Δ {y the y}_{τi τi})) \\ 00 = = {C C}_{00}^{' '} Δ Δ {x x}^{' '} + + {D D.}_{00}^{' '} Δy Δ y \\ 00 = = {C C}_{τ τ 11}^{' '} Δ Δ {x x}_{τ τ 11}^{' '} + + {D D.}_{τ τ 11}^{' '} Δ Δ {y the y}_{τ τ 11} \\ . . \\ . . \\ . . \\ 00 = = {C C}_{τm τm}^{' '} Δ Δ {x x}_{τm τm}^{' '} + + {D D.}_{τm τm}^{' '} Δ Δ {y the y}_{τm τm} \end{matrix}- - - - - - - ((66))$

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

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

In the formula: Α ₀ ′, B′ ₀ , C′ ₀ , D′ ₀ represent the partial derivatives of differential equation f′ and algebraic equation g to state variable x′ and algebraic variable y respectively,

Α′ _τi , B′ _τi , C′ _τi , D′ _τi represent the partial derivatives of the differential equation f′ and the algebraic equation g to the time-delay state variable x′ _τi and the algebraic variable y _τi respectively

时滞电力系统特征值求解问题Problems of Solving Eigenvalues of Time-delayed Power Systems

当D₀′和D′_τi(i=1,2,…,m)非奇异时，消去代数变量，方程(6)可简化为：When D ₀ ′ and D′ _τi (i=1,2,…,m) are non-singular, the algebraic variables are eliminated, and equation (6) can be simplified as:

$Δ Δ {\overset{\cdot &Center Dot;}{x x}}^{' '} = = {\overset{~ ~}{A A}}_{00}^{' '} Δ Δ {x x}^{' '} + + {Σ Σ}_{i i = = 11}^{m m} {\overset{~ ~}{A A}}_{τi τi}^{' '} Δ Δ {x x}_{τi τi}^{' '} - - - - - - ((77))$

式中，In the formula,

$\begin{matrix} {\overset{~ ~}{A A}}_{00}^{' '} = = {A A}_{00}^{' '} - - {B B}_{00}^{' '} {(({D D.}_{00}^{' '}))}^{- - 11} {C C}_{00}^{' '} \\ {\overset{~ ~}{A A}}_{τi τi}^{' '} = = {A A}_{τi τi}^{' '} - - {B B}_{τi τi}^{' '} {(({D D.}_{τi τi}^{' '}))}^{- - 11} {C C}_{τi τi}^{' '} \end{matrix} - - - - - - ((88))$

式(7)表示的线性化系统的特征方程为：The characteristic equation of the linearized system represented by formula (7) is:

$det det ((λI λ I - - {\overset{~ ~}{A A}}_{00}^{' '} - - {Σ Σ}_{i i = = 11}^{m m} {\overset{~ ~}{A A}}_{τi τi}^{' '} {e e}^{- - λ λ {τ τ}_{i i}})) = = 00 - - - - - - ((99))$

式中，λ为系统的特征值。In the formula, λ is the eigenvalue of the system.

时滞动力系统的稳定性理论指出[廖晓昕.动力系统的稳定性理论和应用[M].北京:国防工业出版社,2000.]，如果线性化系统(7)的全部特征值都具有负实部，则时滞系统(5)在稳态运行点(x′₀,y₀)处是小扰动稳定的；反之，若至少存在一个具有正实部的特征值，则系统在该点处是小扰动不稳定的。The stability theory of time-delay dynamical system points out [Liao Xiaoxin. Stability Theory and Application of Dynamical System [M]. Beijing: National Defense Industry Press, 2000.], if all the eigenvalues of the linearized system (7) have negative real part, then the time-delay system (5) is small-disturbance stable at the steady-state operating point (x′ ₀ , y ₀ ); on the contrary, if there is at least one eigenvalue with a positive real part, the system at this point is Small perturbations are unstable.

然而，由于时滞及指数项

的存在，线性化系统的特征方程(9)为超越方程，其有无穷多个解。因此，通过直接求解该方程得到系统的关键特征值并判别系统的时滞稳定性变得非常困难。However, due to the time lag and the exponential term

The existence of , the characteristic equation (9) of the linearized system is a transcendental equation, which has infinitely many solutions. Therefore, it becomes very difficult to obtain the key eigenvalues of the system and judge the time-delay stability of the system by directly solving the equation.

中国专利201010123345.8，提出了一种直接、有效搜索复平面上特定边界上（如虚轴、特征值实部或阻尼比等于给定常数）时滞系统的关键特征值的方法。该方法能够得到的精确的特征值，但搜索过程的计算量较大。Chinese patent 201010123345.8 proposes a method to directly and efficiently search for key eigenvalues of time-delay systems on specific boundaries (such as imaginary axes, real parts of eigenvalues, or damping ratios equal to a given constant) on the complex plane. This method can obtain accurate eigenvalues, but the calculation of the search process is relatively large.

中国专利200810151217.7、200910070255.4、200910070254.X均是基于线性矩阵不等式（Linear Matrix Inequlity，LMI）方法来判别时滞系统的小干扰稳定性。这类方法在判别系统的稳定性时存在固有的保守性。Chinese patents 200810151217.7, 200910070255.4, and 200910070254.X are all based on the linear matrix inequality (Linear Matrix Inequlity, LMI) method to determine the small disturbance stability of the time-delay system. Such methods are inherently conservative in judging the stability of the system.

发明内容 Contents of the invention

本发明的目的就是为了解决上述问题，提供一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，它不需要迭代和搜索就可以准确地求解与时滞系统中动态元件相关的部分特征值和特征向量，利用计算得到的特征值可以直接判断系统的小干扰稳定性，具有计算量小，计算时间少，不存在任何保守性的优点。The purpose of the present invention is exactly in order to solve above-mentioned problem, provides a kind of time-delay power system eigenvalue calculation and stability discrimination method based on Pade approximation, it does not need iteration and search and just can accurately solve the dynamic element relevant in time-delay system Some eigenvalues and eigenvectors can be used to directly judge the small-disturbance stability of the system by using the calculated eigenvalues, which has the advantages of small calculation amount, less calculation time, and no conservatism.

为了实现上述目的，本发明采用如下技术方案：In order to achieve the above object, the present invention adopts the following technical solutions:

一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，具体步骤如下：A Pade approximation based eigenvalue calculation and stability discrimination method for time-delayed power systems, the specific steps are as follows:

步骤一：指定广域反馈时滞τ_i；Step 1: specify wide-area feedback delay τ _i ;

步骤二：进行Pade近似，得到时滞环节的近似有理多项式；Step 2: Perform Pade approximation to obtain the approximate rational polynomial of the delay link;

步骤三：将有理多项式转化为状态空间表达式；Step 3: Transform rational polynomials into state space expressions;

步骤四：进行平衡化处理；Step 4: Perform balance processing;

步骤五：将时滞环节与无时滞电力系统、广域阻尼控制器进行连接，得到时滞电力系统的线性化模型；Step 5: Connect the time-delay link with the time-delay-free power system and the wide-area damping controller to obtain the linearized model of the time-delay power system;

步骤六：利用QR法或稀疏特征值算法求解系统的特征值；Step 6: Use the QR method or the sparse eigenvalue algorithm to solve the eigenvalues of the system;

步骤七：直接判别系统的时滞稳定性。Step seven: directly judge the time-delay stability of the system.

所述步骤二中，在拉普拉斯域中，第i个时滞环节可表示为

t是时间，s是频率。Pade近似是一种利用[l,k]阶有理多项式逼近

的方法：In the second step, in the Laplace domain, the i-th time-delay link can be expressed as

t is time and s is frequency. The Pade approximation is a rational polynomial approximation using [l,k] order

Methods:

${e e}^{- - {τ τ}_{i i} s the s} \approx \approx R R ((s the s)) = = \frac{{N N}_{l l} ((s the s))}{{N N}_{k k} ((s the s))} = = \frac{{a a}_{00} + + {a a}_{11} {τ τ}_{i i} s the s + + . . . . . . + + {a a}_{l l} {(({τ τ}_{i i} s the s))}^{l l}}{} - - - - - - ((1010))$

式中，l和k为正整数。系数a_i和b_i为实数，可以由下式求出：In the formula, l and k are positive integers. The coefficients a _i and b _i are real numbers and can be obtained by the following formula:

${a a}_{j j} = = {((- - 11))}^{j j} \frac{((l l + + k k - - j j))!! l l!!}{j j!! ((l l - - j j))!!} - - - - - - ((1111))$

${b b}_{j j} = = \frac{((l l + + k k - - j j))!! k k!!}{j j!! ((k k - - j j))!!} - - - - - - ((1212))$

所述步骤二在Pade近似中，阶数l和k越大，有理多项式R(s)越接近于

通常情况下，取l=k；时滞越小，频域中R(s)与

相位一致的区间越大，即频带越宽。In the Pade approximation, the step 2, the larger the order l and k, the closer the rational polynomial R(s) is to

Usually, l=k is taken; the smaller the time lag, the R(s) and

The larger the phase coincident interval, that is, the wider the frequency band.

所述步骤三中，将式(10)的分子和分母同时除以

进而利用传递函数实现问题的一般性方法[刘豹.现代控制理论[M].北京:机械工业出版社,2000.]，将基于Pade近似得到的时滞环节的有理多项式逼近R(s)转换为等价的状态空间表达，对于第i个时滞环节，以变量增量形式表示的能控标准型为：In described step 3, the numerator and the denominator of formula (10) are divided by

Then use the transfer function to realize the general method of the problem [Liu Bao. Modern Control Theory [M]. Beijing: Machinery Industry Press, 2000.], transform the rational polynomial approximation R(s) based on the time-delay link obtained by Pade approximation is an equivalent state space expression, for the i-th time-delay link, the controllable standard form expressed in the form of variable increment is:

$\{\begin{matrix} Δ Δ {\overset{\cdot \cdot}{x x}}_{di di} = = {\overset{~ ~}{A A}}_{di di} Δ Δ {x x}_{di di} + + {\overset{~ ~}{B B}}_{d d} Δ Δ {u u}_{i i} \\ Δ Δ {y the y}_{di di} = = {\overset{~ ~}{C C}}_{di di} Δ Δ {x x}_{di di} + + {\overset{~ ~}{D D.}}_{di di} Δ Δ {u u}_{i i} \end{matrix} - - - - - - ((1313))$

式（13）中，系数矩阵

的具体表达式为：In formula (13), the coefficient matrix

The specific expression is:

k为Pade近似阶数。k is the approximate order of Pade.

在较小时滞τ_i和较高阶数的情况下，由Pade近似得到的传递函数R(s)中的系数之间，以及由传递函数转换得到的状态空间表达式的系数矩阵A_di和C_di的非零元素之间，在数量级上相差很大，采用线性变换使得时滞环节的状态空间表达式的系数矩阵更加平衡：In the case of small delay τ _i and higher order, between the coefficients in the transfer function R(s) approximated by Pade, and the coefficient matrices A _di and C of the state space expression obtained by the transfer function transformation There is a large difference in magnitude between the non-zero elements of _di , and the linear transformation makes the coefficient matrix of the state-space expression of the time-delay link more balanced:

${A A}_{di di} = = {T T}^{- - 11} {\overset{~ ~}{A A}}_{di di} T T,, {B B}_{di di} = = {T T}^{- - 11} {\overset{~ ~}{B B}}_{di di},, {C C}_{di di} = = {\overset{~ ~}{C C}}_{di di} T T,, {D D.}_{di di} = = {\overset{~ ~}{D D.}}_{di di} - - - - - - ((1515))$

式中，T为对角变换矩阵；由此可知，平衡化处理后的系数矩阵仍然符合能控标准型；In the formula, T is a diagonal transformation matrix; it can be seen that the coefficient matrix after balancing processing still conforms to the controllable standard form;

所述步骤五：将步骤三中得到的时滞环节的状态空间表达式，与无时滞电力系统的线性化模型、广域阻尼控制器的线性化模型相连接，从而建立包含时滞环节的电力系统小扰动稳定性分析的线性化模型：Said step five: connect the state space expression of the time-delay link obtained in step three with the linearization model of the power system without time-delay and the linearization model of the wide-area damping controller, thereby establishing a time-delay link A linearized model for small-disturbance stability analysis of power systems:

$\{\begin{matrix} Δ Δ {\overset{\cdot &Center Dot;}{x x}}^{' '' '} = = {A A}^{' '' '} Δ Δ {x x}^{' '' '} + + {B B}^{' '' '} Δy Δ y \\ 00 = = {C C}^{' '' '} Δ Δ {x x}^{' '' '} + + {D D.}^{' '' '} Δy Δy \end{matrix} - - - - - - ((1616))$

式中，x″为系统状态变量，y为系统代数变量（节点电压），A"、B"、C"、D"为系数矩阵；In the formula, x" is the system state variable, y is the system algebraic variable (node voltage), and A", B", C", D" are coefficient matrices;

假设第i个时滞环节的状态变量x_d和与之相连的广域阻尼控制器的状态变量x_c排列在无时滞电力系统状态变量x之后，即 $Δ x^{''} = {[\begin{matrix} {Δx}^{T} & . & . & . & {Δx}_{di}^{T} & {Δx}_{ci}^{T} & . & . & . \end{matrix}]}^{T};$ 设无时滞电力系统的系数矩阵为A，B，C，D；广域阻尼控制器的系数矩阵为A_c，B_c，C_c，D_c，广域反馈信号的系数矩阵分别为K_1i和K_2i；于是，式(16)中系数矩阵A"、B"、C"、D"具体可表示为：Assume that the state variable x _d of the i-th time-delay link and the state variable x _c of the wide-area damping controller connected to it are arranged after the state variable x of the time-delay-free power system, that is $Δ x^{''} = {[\begin{matrix} {Δx}^{T} & . & . & . & {Δx}_{di}^{T} & {Δx}_{ci}^{T} & . & . & . \end{matrix}]}^{T};$ Suppose the coefficient matrix of the time-delay-free power system is A, B, C, D; the coefficient matrix of the wide-area damping controller is A _c , B _c , C _c , D _c , and the coefficient matrix of the wide-area feedback signal is K _1i and K _2i ; then, the coefficient matrices A", B", C", and D" in formula (16) can be specifically expressed as:

式(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 signalstability analysis[J].Proceedings of the CSEE,2005,25(2):17-21.]；E为稀疏矩阵，Aci、Bci、Cci、Dci与广域阻尼控制器的具体结构有关；In formula (17), the coefficient matrices A _di , B _di , C _di , and D _di satisfy the controllable standard type; A is a block diagonal matrix, B and C are block sparse matrices, and D is 2×2 block sparse Matrix[Du Zhengchun, Liu Wei, Fang Wanliang, et al. Sparse implementation of a key eigenvalue calculation in small disturbance stability analysis[J]. Chinese Journal of Electrical Engineering, 2005,25(2):17-21.Du Zhengchun, Liu Wei, Fang Wanliang, et al.A sparse method for the calculation of critical eigenvalue in small signalstability analysis[J].Proceedings of the CSEE,2005,25(2):17-21.]; E is a sparse matrix, Aci , Bci, Cci, Dci are related to the specific structure of the wide-area damping controller;

所述式(17)中各系数矩阵具有如下特点：Each coefficient matrix in the formula (17) has the following characteristics:

(1)若将广域反馈信号中涉及的动态元件、时滞环节、广域阻尼控制器状态变量的系数矩阵作为一个子块，则A"仍为分块对角阵；B"仅在时滞环节和广域阻尼控制器状态变量所在行上增加了零矩阵或少量几列非零元素，B"仍为稀疏阵；(1) If the dynamic components involved in the wide-area feedback signal, the time-delay link, and the coefficient matrix of the state variables of the wide-area damping controller are used as a sub-block, then A" is still a block diagonal matrix; B" only when A zero matrix or a few columns of non-zero elements are added to the row where the hysteresis link and the state variable of the wide-area damping controller are located, and B" is still a sparse matrix;

(2)C"仅在时滞环节和广域阻尼控制器状态变量所在列上增加了零矩阵，D"与无时滞电力系统相应的系数矩阵完全一样。(2) C" only adds a zero matrix on the column where the time-delay link and the state variables of the wide-area damping controller are located, and D" is exactly the same as the corresponding coefficient matrix of the power system without time-delay.

综上可知，包含时滞环节的电力系统的线性化方程系数矩阵A"、B"、C"、D"，与无时滞电力系统在线性化方程系数矩阵A、B、C、D，具有完全相同的稀疏结构。In summary, the linearized equation coefficient matrices A", B", C", and D" of the power system including time-delay links, and the linearized equation coefficient matrices A, B, C, and D of the power system without time-delay, have The exact same sparse structure.

所述步骤六：当D"非奇异时，消去代数变量y，式(16)可简化为：Described step 6: when D " non-singularity, eliminate algebraic variable y, formula (16) can be simplified as:

$Δ Δ {\overset{\cdot &Center Dot;}{x x}}^{' '' '} = = {\overset{~ ~}{A A}}^{' '' '} Δ Δ {x x}^{' '' '} - - - - - - ((1818))$

式中，

为基于Pade近似得到的时滞系统的状态矩阵。In the formula,

is the state matrix of the time-delay system based on the Pade approximation.

通过计算

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

Part of the eigenvalues of the time-delay power system can be obtained. In addition, since the coefficient matrix of the power system with time-delay links and the conventional non-delay power system has the same sparse characteristics, the sparse processing technology used in the calculation of some key eigenvalues of large-scale non-delay power systems[Du Zhengchun, Liu Wei, Fang Wanliang, et al. A Sparse Realization of Calculation of Key Eigenvalues in Small Disturbance Stability Analysis[J]. Proceedings of the Chinese Society for Electrical Engineering, 2005, 25(2): 17-21. DuZhengchun, 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.] and calculation method[Du Zhengchun, Liu Wei, Fang Wanliang, et al. .Calculation of key eigenvalues in small disturbance stability analysis based on Jacobi-Davidson method[J].Proceedings of the Chinese Society for Electrical Engineering,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.], still applicable to the calculation of time-delay power system Eigenvalues.

本发明的有益效果：利用Pade近似将广域反馈时滞逼近为一个有理多项式，通过与无时滞电力系统和广域阻尼控制器之间的连接，建立时滞电力系统的的线性化模型，最后根据系统状态矩阵直接求出时滞系统的部分特征根；该方法能够较准确地求解与时滞系统中动态元件相关的部分特征值和特征向量，正确求解与时滞环节相关的特征值个数和计算精度，与有理多项式的阶数有关。Beneficial effects of the present invention: use Pade approximation to approximate the wide-area feedback time-delay as a rational polynomial, and establish a linearized model of the time-delay power system through the connection between the time-delay-free power system and the wide-area damping controller, Finally, according to the system state matrix, some eigenvalues of the time-delay system are directly obtained; this method can accurately solve some eigenvalues and eigenvectors related to the dynamic components in the time-delay system, and correctly solve the eigenvalues related to the time-delay link. The number and calculation precision are related to the order of the rational polynomial.

利用Pade近似有理多项式来逼近时滞环节，进而计算得到含有广域通信时滞电力系统的部分特征值，从而直接判断电力系统的小干扰稳定性。The Pade approximate rational polynomial is used to approximate the time-delay link, and then some eigenvalues of the power system with wide-area communication time-delay are calculated, so as to directly judge the small-disturbance stability of the power system.

附图说明 Description of drawings

图1为本发明的整体流程图；Fig. 1 is the overall flowchart of the present invention;

图2为时滞电力系统示意图；Figure 2 is a schematic diagram of a time-delay power system;

图3为两区四机系统图；Figure 3 is a system diagram of two zones and four machines;

图4为实部大于-50的部分特征值；Figure 4 shows some eigenvalues whose real part is greater than -50;

图5为实部大于-10的部分特征值。Figure 5 shows some eigenvalues whose real part is greater than -10.

其中，1.无时滞电力系统，2.e^-sτ，3.广域阻尼控制器。Among them, 1. Delay-free power system, 2. e ^-sτ , 3. Wide-area damping controller.

具体实施方式 Detailed ways

下面结合附图与实施例对本发明作进一步说明。The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

如图1所示，一种基于Pade近似的时滞电力系统特征值计算与稳定性判别方法，具体步骤如下：As shown in Fig. 1, a Pade approximation-based eigenvalue calculation and stability discrimination method for time-delayed power systems, the specific steps are as follows:

步骤一：指定广域反馈时滞τ_i和Pade近似阶数k；Step 1: specify wide-area feedback delay τ _i and Pade approximation order k;

步骤四：进行平衡化处理；Step 4: Perform balance processing;

所述步骤二中，在拉普拉斯域中，第i个时滞环节

可表示为

Pade近似是一种利用[l,k]阶有理多项式逼近

的方法：In the second step, in the Laplace domain, the i-th time-delay link

can be expressed as

The Pade approximation is a rational polynomial approximation using [l,k] order

Methods:

${e e}^{- - {τ τ}_{i i} s the s} \approx \approx R R ((s the s)) = = \frac{{N N}_{l l} ((s the s))}{{N N}_{k k} ((s the s))} = = \frac{{a a}_{00} + + {a a}_{11} {τ τ}_{i i} s the s + + . . . . . . + + {a a}_{l l} {(({τ τ}_{i i} s the s))}^{l l}}{} - - - - - - ((1919))$

式中，系数ai和bi可以由下式求出：In the formula, the coefficients ai and bi can be obtained by the following formula:

${a a}_{j j} = = {((- - 11))}^{j j} \frac{((l l + + k k - - j j))!! l l!!}{j j!! ((l l - - j j))!!} - - - - - - ((2020))$

${b b}_{j j} = = \frac{((l l + + k k - - j j))!! k k!!}{j j!! ((k k - - j j))!!} - - - - - - ((21 twenty one))$

所述步骤二在Pade近似中，阶数l和k越大，有理多项式R(s)越接近于通常情况下，取l=k；时滞越小，频域中R(s)与

相位一致的区间越大，即频带越宽。In the Pade approximation, the step 2, the larger the order l and k, the closer the rational polynomial R(s) is to Usually, l=k is taken; the smaller the time lag, the R(s) and

所述步骤三中，将式(10)的分子和分母同时除以

$\{\begin{matrix} Δ Δ {\overset{\cdot \cdot}{x x}}_{di di} = = {\overset{~ ~}{A A}}_{di di} Δ Δ {x x}_{di di} + + {\overset{~ ~}{B B}}_{d d} Δ Δ {u u}_{i i} \\ Δ Δ {y the y}_{di di} = = {\overset{~ ~}{C C}}_{di di} Δ Δ {x x}_{di di} + + {\overset{~ ~}{D D.}}_{di di} Δ Δ {u u}_{i i} \end{matrix} - - - - - - ((22 twenty two))$

式（13）中，In formula (13),

${A A}_{di di} = = {T T}^{- - 11} {\overset{~ ~}{A A}}_{di di} T T,, {B B}_{di di} = = {T T}^{- - 11} {\overset{~ ~}{B B}}_{di di},, {C C}_{di di} = = {\overset{~ ~}{C C}}_{di di} T T,, {D D.}_{di di} = = {\overset{~ ~}{D D.}}_{di di} - - - - - - ((24 twenty four))$

所述步骤五：将步骤三中得到的时滞环节的状态空间表达式，与无时滞电力系统的线性化模型、广域阻尼控制器的线性化模型相连接，从而建立包含时滞环节的电力系统小扰动稳定性分析的线性化模型：Said step five: connect the state space expression of the time-delay link obtained in step three with the linearized model of the power system without time-delay and the linearized model of the wide-area damping controller, thereby establishing a time-delay link A linearized model for small-disturbance stability analysis of power systems:

$\{\begin{matrix} Δ Δ {\overset{\cdot \cdot}{x x}}^{' '' '} = = {A A}^{' '' '} Δ Δ {x x}^{' '' '} + + {B B}^{' '' '} Δy Δ y \\ 00 = = {C C}^{' '' '} Δ Δ {x x}^{' '' '} + + {D D.}^{' '' '} Δy Δy \end{matrix} - - - - - - ((2525))$

(1)若将广域反馈信号中涉及的动态元件、时滞环节、广域阻尼控制器状态变量的系数矩阵作为一个子块，则A"仍为分块对角阵；B"仅在时滞环节和广域阻尼控制器状态变量所在行上增加了零矩阵或少量几列非零元素，B"仍为稀疏阵；(1) If the dynamic components involved in the wide-area feedback signal, the time-delay link, and the coefficient matrix of the state variable of the wide-area damping controller are used as a sub-block, then A" is still a block diagonal matrix; B" only when A zero matrix or a few columns of non-zero elements are added to the row where the hysteresis link and the state variable of the wide-area damping controller are located, and B" is still a sparse matrix;

$Δ Δ {\overset{\cdot &Center Dot;}{x x}}^{' '' '} = = {\overset{~ ~}{A A}}^{' '' '} Δ Δ {x x}^{' '' '} - - - - - - ((2727))$

式中，

为基于Pade近似得到的时滞系统的状态矩阵。In the formula,

is the state matrix of the time-delay system based on the Pade approximation.

通过计算

广域测量信号在传输和处理过程中产生的时滞，使电力系统成为一个时滞系统。时滞电力系统各部分之间的连接关系如图2所示，无时滞电力系统与e^-sτ连接，e^-sτ与广域阻尼控制器连接。The time delay generated during the transmission and processing of the wide-area measurement signal makes the power system a time-delay system. The connection relationship between various parts of the time-delay power system is shown in Figure 2. The time-delay-free power system is connected to e ^-sτ , and e ^-sτ is connected to the wide-area damping controller.

如图3所示的两区四机电力系统为例，验证本发明提出的时滞电力系统特征值计算方法的正确性和有效性,系统的参数详见[Kunder P.Power System Stabilityand Control[M].NewYork:McGraw-Hill,1994.]。Take the two-zone four-machine power system as shown in Figure 3 as an example to verify the correctness and effectiveness of the time-delay power system eigenvalue calculation method proposed by the present invention. The parameters of the system are detailed in [Kunder P.Power System Stability and Control[M ]. New York: McGraw-Hill, 1994.].

在发电机G1和G3安装本地电力系统稳定器（Power System Stabilizer,PSS）的基础上，考虑在G1上装设以G1和G3相对转速偏差Δω13为反馈信号的广域PSS，进一步提高系统的阻尼水平。On the basis of installing local power system stabilizers (Power System Stabilizer, PSS) on generators G1 and G3, it is considered to install a wide-area PSS on G1 with the relative speed deviation Δω13 of G1 and G3 as the feedback signal to further improve the damping level of the system .

本发明以基于离散化的时滞系统特征值求解方法计算得到的部分特征值作为精确解，并作为验证本发明提出的基于Pade近似的特征值计算方法有效性和精确性的基准。The present invention uses part of the eigenvalues calculated based on the discretization-based eigenvalue solving method of the time-delay system as the exact solution, and serves as a benchmark for verifying the effectiveness and accuracy of the Pade approximation-based eigenvalue calculation method proposed by the present invention.

本发明借助Matlab工具箱DDE-BIFTOOL中p_stabil函数计算得到的算例系统的部分特征值的精确解，其中，牛顿法迭代的收敛精度取为1e-8。The present invention uses the p_stabil function in the Matlab toolbox DDE-BIFTOOL to calculate the exact solution of some eigenvalues of the example system, wherein the convergence accuracy of Newton's method iteration is taken as 1e-8.

针对算例系统，在不同的反馈时滞τ和不同的Pade近似有理多项式的阶数k下，利用本发明提出的基于Pade近似的特征值计算方法进行了大量的计算，并与基于离散化的特征值计算方法的结果进行了对比和分析，下面仅以部分计算结果来说明本发明方法的正确性和有效性。For the example system, under different feedback delays τ and different order k of Pade approximation rational polynomials, a large number of calculations are performed using the Pade approximation-based eigenvalue calculation method proposed by the present invention, and combined with the discretization-based The results of the eigenvalue calculation method have been compared and analyzed, and the correctness and effectiveness of the method of the present invention will be illustrated with only part of the calculation results below.

当反馈时滞τ=0.30s时，利用基于离散化的特征值计算方法、基于5阶和20阶Pade近似的特征值计算方法，得到的实部大于-50的部分特征值，如图4所示。图4右下角矩形框内的、实部大于-10的部分特征值的详细分布如图5所示，由图5可知：When the feedback delay τ=0.30s, using the discretization-based eigenvalue calculation method and the eigenvalue calculation method based on the 5th order and 20th order Pade approximation, some eigenvalues whose real part is greater than -50 are obtained, as shown in Figure 4 Show. The detailed distribution of some eigenvalues in the rectangular box in the lower right corner of Figure 4, whose real part is greater than -10, is shown in Figure 5. From Figure 5, we can see that:

(1)当k=5和20时，利用基于Pade近似的特征值计算方法，能够较为准确地计算得到实部大于-10的复平面内的部分特征值，通过进一步分析可知，其与描述系统动态元件的状态变量对应。(1) When k=5 and 20, the eigenvalue calculation method based on the Pade approximation can be used to calculate more accurately some eigenvalues in the complex plane whose real part is greater than -10. Through further analysis, it can be seen that it is consistent with the description system Corresponding to the state variable of the dynamic element.

(2)对于k=5，在实部小于-10的复平面内，基于Pade近似的特征值计算方法，会得到若干个错误的特征值，如：-15.06761008±23.47922927i、-21.53755703±15.83642620i、-36.68588408、-49.61254513；此外，还有多个特征值未被计算出来，通过进一步分析可知，遗漏和错误求解的特征值，主要与有理多项式的各阶变量对应。(2) For k=5, in the complex plane whose real part is less than -10, based on the Pade approximation eigenvalue calculation method, several wrong eigenvalues will be obtained, such as: -15.06761008±23.47922927i, -21.53755703±15.83642620i , -36.68588408, -49.61254513; In addition, there are many eigenvalues that have not been calculated. Through further analysis, it can be seen that the eigenvalues that are missing and wrongly solved mainly correspond to the variables of each order of the rational polynomial.

(3)只有当有理多项式的阶数显著增大到k=20时，才能保证本发明方法在实部大于-41的复平面内不出现漏根、错根的情况，与有理多项式变量对应的部分特征值也能较为准确地求得。当阶数增加到30时，在实部大于-50的复平面内，本发明与基于离散化的特征值求解方法得到特征值完全相同。(3) only when the order of rational polynomial increases significantly to k=20, just can guarantee that the inventive method does not occur the situation of leakage root, wrong root in the complex plane that real part is greater than-41, and the corresponding of rational polynomial variable Some eigenvalues can also be obtained more accurately. When the order increases to 30, in the complex plane whose real part is greater than -50, the present invention is exactly the same as the eigenvalue obtained by the discretization-based eigenvalue solving method.

表1τ=0.3s时，基于离散化和Pade近似计算得到的特征值Table 1 When τ=0.3s, the eigenvalues calculated based on discretization and Pade approximation

(4)在增加Pade近似有理多项式阶数以准确求取与有理多项式变量对应的部分特征值的同时，与系统动态元件相对应的部分特征值的精度也显著地得到提高。如表1所示，例如，当k=5时，本发明方法与基于离散化方法得到的特征值之间的最大绝对误差max(abs(ΔRe(λ),ΔIm(λ)))=1.042e-4，当k增加到7时，最大绝对误差相应地减小到4.673e-8，特征值计算精度得到显著的提高。(4) While increasing the order of Pade approximate rational polynomials to accurately obtain some eigenvalues corresponding to rational polynomial variables, the accuracy of some eigenvalues corresponding to system dynamic components is also significantly improved. As shown in Table 1, for example, when k=5, the maximum absolute error max(abs(ΔRe(λ),ΔIm(λ)))=1.042e between the method of the present invention and the eigenvalue obtained based on the discretization method -4, when k increases to 7, the maximum absolute error is correspondingly reduced to 4.673e-8, and the calculation accuracy of eigenvalues is significantly improved.

(5)由计算得到的系统的特征值可知，当τ=0.3s时，系统是小干扰稳定的。(5) According to the calculated eigenvalues of the system, when τ=0.3s, the system is stable with little disturbance.

为了深入考察Pade近似有理多项式的阶数k对特征值计算精度的影响，在不同时滞值下，基于不同阶数Pade近似的特征值计算方法得到的特征值，与基于离散化方法得到的特征值之间的最大绝对误差如表2所示。不难看出，在给定精度要求下，随着时滞的增大，Pade近似有理多项式阶数k须相应地增大。由表2中可知，对于4机2区域系统，τ在[0.05，0.5]范围内变化时，6阶Pade近似就能使本发明方法计算得到的特征值达到1e-5的计算精度。In order to further investigate the influence of the order k of the Pade approximation rational polynomial on the calculation accuracy of the eigenvalues, under different delay values, the eigenvalues obtained by the eigenvalue calculation method based on the Pade approximation of different orders are compared with the eigenvalues obtained by the discretization method. The maximum absolute error between the values is shown in Table 2. It is not difficult to see that, under a given precision requirement, as the delay increases, the order k of the Pade approximation rational polynomial must increase accordingly. It can be seen from Table 2 that for a 4-machine 2-area system, when τ varies in the range of [0.05, 0.5], the 6th-order Pade approximation can make the eigenvalues calculated by the method of the present invention reach the calculation accuracy of 1e-5.

表2不同时滞和Pade近似阶数下，特征值的计算精度Table 2 Calculation accuracy of eigenvalues under different delays and Pade approximation orders

当反馈时滞τ=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模态分量的模值较大。When the feedback time lag τ=0.30, the exact value of the right eigenvector corresponding to the system partial eigenvalue obtained based on the discretized eigenvalue calculation method is shown in rows 2 to 4 of Table 3. Through further modal analysis, the three electromechanical oscillation modes and corresponding modes of the system can be obtained, as shown in rows 2 to 4 of Table 3. Among them, λ1 represents the oscillation of G1 and G2 relative to G3 and G4; λ2 represents the oscillation of G1 relative to G2 and G3 relative to G4, but the modulus value corresponding to the modal components of G1 and G2 is relatively large; λ3 represents the relative Due to the oscillation of G2 and G3 relative to G4, but the modulus values corresponding to the modal components of G3 and G4 are relatively large.

利用5阶和7阶Pade近似的计算得到的系统部分特征值和相应的右特征向量中与各发电机转速相应的分量，如表3的第5~10行所示。通过与模态的精确值的对比，可知本发明可以准确计算得到振荡模态的幅值和相位。k=5时，模态幅值的最大绝对误差为4e-8，相位的最大绝对误差为1e-2°；k=7时，幅值的最大绝度误差为1e-8；模态相位的最大绝对误差为9e-6°。The partial eigenvalues of the system obtained by using the 5th-order and 7th-order Pade approximation and the corresponding components in the corresponding right eigenvectors corresponding to the rotation speed of each generator are shown in rows 5 to 10 of Table 3. Through comparison with the precise value of the mode, it can be seen that the present invention can accurately calculate the amplitude and phase of the oscillation mode. When k=5, the maximum absolute error of modal amplitude is 4e-8, and the maximum absolute error of phase is 1e-2°; when k=7, the maximum absolute error of amplitude is 1e-8; The maximum absolute error is 9e-6°.

表3τ=0.3s时，Pade近似计算得到的模态及其精确值Table 3 When τ=0.3s, the modes and their exact values obtained by Pade’s approximate calculation

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。Although the specific implementation of the present invention has been described above in conjunction with the accompanying drawings, it does not limit the protection scope of the present invention. Those skilled in the art should understand that on the basis of the technical solution of the present invention, those skilled in the art do not need to pay creative work Various modifications or variations that can be made are still within the protection scope of the present invention.

Claims

1. calculate and a stability distinguishing method based on the approximate time-lag power system characteristic value of Pade, it is characterized in that, concrete steps are as follows:

Step 1: specify wide area Feedback Delays τ _iwith the approximate exponent number k of Pade;

Step 2: carry out Pade and be similar to, obtain the approximate rational polynominal of Time Delay; In described step 2, in Laplace domain, i Time Delay

be expressed as

t is the time, and s is frequency; The value of i is positive integer; Pade is approximate is that a kind of [l, k] rank rational polynominal of utilizing is approached

method:

e^{- τ_{i} s} \approx R (s) = \frac{N_{l} (s)}{N_{k} (s)} = \frac{a_{0} + a_{1} τ_{i} s + . . . + a_{l} {(τ_{i} s)}^{l}}{b_{0} + b_{1} τ_{i} s + . . . + b_{k} {(τ_{i} s)}^{k}} - - - (10)

In formula, l and k are positive integer; Coefficient a _mand b _jobtained by following formula:

a_{m} = {(- 1)}^{m} \frac{(l + k - m)! l!}{m! (l - m)!} - - - (11);

b_{j} = \frac{(l + k - j)! k!}{j! (k - j)!} - - - (12);

Wherein, 0≤m≤l, 0≤j≤k;

Step 3: rational polynominal is converted into state-space expression;

Step 4: carry out equilibrating processing;

Step 5: by Time Delay be connected without time-lag power system, wide area damping control, obtain the inearized model of time-lag power system;

Step 6: the characteristic value of utilizing QR method or sparse eigenvalue Algorithm for Solving system;

Step 7: the directly time lag stability of judgement system.

2. as claimed in claim 1 a kind ofly calculate and stability distinguishing method based on the approximate time-lag power system characteristic value of Pade, is characterized in that, described step 2 is in Pade is approximate, and exponent number l and k are larger, rational polynominal R (s) more close to

under normal circumstances, get l=k; Time lag is less, R in frequency domain (s) with

the interval that phase place is consistent is larger, and frequency band is wider.

3. as claimed in claim 1 a kind ofly calculate and stability distinguishing method based on the approximate time-lag power system characteristic value of Pade, is characterized in that, in described step 3, by the molecule of formula (10) and denominator simultaneously divided by

and then utilize the general approach of transfer function problem of implementation, rational polynominal based on the approximate Time Delay obtaining of Pade is approached to R (s) and be converted to state space expression of equal value, for i Time Delay, represent take variable incremental form can control standard type as:

\{\begin{matrix} Δ {\dot{x}}_{di} = {\tilde{A}}_{di} Δ x_{di} + {\tilde{B}}_{di} Δ u_{i} \\ Δ y_{di} = {\tilde{C}}_{di} Δ x_{di} + {\tilde{D}}_{di} Δ u_{i} \end{matrix} - - - (13)

In formula (13), coefficient matrix

expression be:

4. as claimed in claim 1 a kind of based on the Pade approximate calculating of time-lag power system characteristic value and stability distinguishing method, it is characterized in that, described step 4, in the situation that counting compared with Small Time Lag τ i and higher-order, between the coefficient in the approximate transfer function R (s) obtaining of Pade, and the coefficient matrices A of the state-space expression being converted to by transfer function _diand C _dinonzero element between, on the order of magnitude, differ greatly, adopt linear transformation to make the coefficient matrix balance more of the state-space expression of Time Delay:

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

In formula, T is diagonal transformation matrix; Hence one can see that, and equilibrating coefficient matrix after treatment still meets can control standard type.

5. as claimed in claim 1 a kind of based on the Pade approximate calculating of time-lag power system characteristic value and stability distinguishing method, it is characterized in that, described step 5: by the state-space expression of the Time Delay obtaining in step 3, be connected with the inearized model without time-lag power system, the inearized model of wide area damping control, thereby set up the inearized model of the Study of Power System Small Disturbance Stability analysis that comprises Time Delay:

\{\begin{matrix} Δ {\dot{x}}^{'} = A^{''} Δ x^{'} + B^{'} Δy \\ 0 = C^{'} Δ x^{'} + D^{'} Δy \end{matrix} - - - (16)

" be system state variables, y is system algebraically variable, and A ", B ", C ", D " are coefficient matrix in formula, x;

The state variable x of the wide area damping control of supposing the state variable xd of i Time Delay and be attached thereto _cbe arranged in without after time-lag power system state variable x,

Δ x^{'} = {[\begin{matrix} Δ x^{T} & \cdot \cdot \cdot & Δ x_{di}^{T} & Δ x_{ci}^{T} & \cdot \cdot \cdot \end{matrix}]}^{T};

If be A without the coefficient matrix of time-lag power system, B, C, D; The coefficient matrix of wide area damping control is A _c, B _c, C _c, D _c, the coefficient matrix of wide area feedback signal is respectively K _1iand K _2i; So coefficient matrices A in formula (16) ", B ", C ", D " are specifically expressed as:

In formula (17), coefficient matrices A _di, B _di, C _di, D _disatisfiedly can control standard type; A is a point Block diagonal matrix, and B, C are piecemeal sparse matrix, and D is 2 × 2 piecemeal sparse matrixes; E is sparse matrix, A _ci, B _ci, C _ci, D _cirelevant with the concrete structure of wide area damping control.

6. as claimed in claim 5 a kind ofly calculate and stability distinguishing method based on the approximate time-lag power system characteristic value of Pade, is characterized in that, in described step 5 formula (17), each coefficient matrix has following features:

(1) " be still a point Block diagonal matrix if using the coefficient matrix of the dynamic element relating in wide area feedback signal, Time Delay, wide area damping control state variable as sub-block, an A; " only on Time Delay and wide area damping control state variable are expert at, increased null matrix or a small amount of several row nonzero element, B " is still Sparse Array to B;

(2) C " has only increased null matrix, D " with just the same without the corresponding coefficient matrix of time-lag power system in Time Delay and wide area damping control state variable column;

In summary, the lienarized equation coefficient matrices A of the electric power system that comprises Time Delay ", B ", C ", D ", with without time-lag power system at linearisation equation coefficient matrix A, B, C, D, there is identical sparsity structure.

7. as claimed in claim 1 a kind ofly calculate and stability distinguishing method based on the approximate time-lag power system characteristic value of Pade, is characterized in that described step 6: as D " when nonsingular, cancellation algebraically variable y, formula (16) is reduced to:

Δ {\dot{x}}^{'} = {\tilde{A}}^{'} Δ x^{'} - - - (18)

In formula,

for the state matrix based on the approximate time-lag system obtaining of Pade;

By calculating

characteristic value, must arrive the partial feature value of time-lag power system; In addition, because the electric power system that comprises Time Delay has identical sparse characteristic with conventional without the coefficient matrix of time-lag power system, the sparse treatment technology and the computational methods that in the time calculating on a large scale without time-lag power system part critical eigenvalue, use, stand good in the partial feature value of calculating time-lag power system.