CN101777767B - 一种时滞电力系统稳定的判别方法 - Google Patents

Abstract

本发明属于电力系统技术领域，涉及一种时滞电力系统稳定的判别方法：首先利用参数变换技术，将大范围时滞稳定域边界的求解，转换为对有限参数空间的搜索和平移过程，大大减少了时滞稳定域边界的计算量；进一步，利用四点插值法，通过对稳定域边界上临界点实施微扰以确定稳定域的构成和边界局部性质。本发明仅需要在有限范围内追踪一个复矩阵的特征值轨迹，即可求解系统在大范围时滞空间中的稳定域及其边界，计算效率较高。

Description

一种时滞电力系统稳定的判别方法

技术领域

本发明属于电力系统技术领域，涉及一种时滞判稳方法。

背景技术

自然界中，系统未来的发展趋势既取决于当前状态，也与过去状态有关，这类现象称为时滞。在广域环境下电力系统的量测时滞十分明显，因此研究时滞环节对系统稳定分析和控制器设计的影响意义重大^[1-5]。

在时滞空间中，电力系统全部稳定运行点构成了系统的稳定域，在进行广域控制器设计时，需保证系统运行点位于稳定域内。已有一些方法探讨时滞空间稳定域的求解，[6]采用频域方法直接求解时滞稳定域边界上的临界点；[7-10]利用Lyapunov方法，通过构造不同的Lyapunov函数来实现求解，但其计算结果存在保守性；[11-12]采用Rekasius变换，通过在(-∞，+∞)内求解一个增广多项式的零点来确定时滞稳定临界点，[13]则通过Lambert-W函数和根轨迹追踪来求解时滞稳定临界点，[11-13]方法的计算量都较大；[14]利用复矩阵变换，通过在有限范围追踪一个复矩阵的特征谱实现时滞稳定临界点的求解，[15]在[14]方法基础上加以改进，使之可以考虑(超)低频振荡现象。

上述方法仅在电力系统时滞较小时有效。而在实际电力系统中，如果通信系统出现了拥塞，控制回路的量测时滞将可能变得很大，此时若广域控制器设计不合理，可能会对系统稳定运行产生危害。因此，需要研究电力系统在大范围时滞空间中稳定域的拓扑构成和边界性质，此时上述方法都难以胜任。

参考文献

1.廖晓昕，动力系统的稳定性理论和应用[M]，国防工业出版社，北京，2000.Liao Xiaoxin，Theory and application of stability for dynamic systems[M]，National Defense IndustryPress，Beijing，2000.

2.Wu Hongxia，Tsakalis K S，Heydt G T，Evaluation of time delay effects to wide-area power systemstabilizer design[J]，IEEE Trans.on Power Systems，2004，19(4)：1935-1941.

3.胡志祥，谢小荣，童陆园，广域阻尼控制延迟特性分析及多项式拟合补偿[J]，电力系统自动化，2005，29(20)：29-34.

Hu Zhixiang，Xie Xiaorong，Tong Luyuan，Characteristic analysis and polynomial fitting basedcompensation of the time delays in wide-area damping control system[J]，Automation of Electric PowerSystems，2005，29(20)：29-34.

4.江全元，邹振宇，曹一家等，考虑时滞影响的电力系统分析和时滞控制研究进展[J]，电力系统自动化，2005，29(3)：2-7.

Jiang Quanyuan，Zou Zhenyu，Cao Yijia，et al，Overview of power systems stability analysis andwide-area control in consideration of time delay[J]，Automation of Electric Power Systems，2005，29(3)：2-7.

5.Yu Xiaodan，Jia Hongjie，Zhao Jinli，A LMI based approach to power system stability analysis with timedelay，Proc of 2008 Tencon2008，Nov.18-21，2008，Hyderabad，India.Vol.1，pp.1-6.

6.E Fridman and M Gil，Stability of linear systems with time-varying delays：a direct frequency domainapproach[J]，Journal of Computational and Applied Mathematics，2007，200(1)：61-66.

7.P Park，A delay-dependent stability criterion for systems with uncertain time-invariant delays[J]，IEEETrans.Automatic Control，1999，44(4)：876-877.

8.Xu Shengyuan and L James，Improved delay-dependent stability criteria for time-delay systems[J]，IEEETrans on Automatic Control，2005，50(3)：384-387.

9.He Yong，Wu Min and She Jinhua，Delay-dependent stability criteria for linear systems with multipletime delays[J]，IEE Proc-Control Theory and Applications，2006，153(4)：447-452.

10.Liu Hailin and Chen Guohua，Delay-dependent stability for neural networks with time-varying delay，Chaos，Solitons & Fractals，2007，33(1)：171-177.

11.Rifat S，Nej at O，A novel stability study on multiple time-delay systems(MTDS)using the rootclustering paradigm[C]，Proc.of the American Control Conference，Boston，MA，2004.6.30-7.2，(6)：5422-5427.

12.贾宏杰，尚蕊和张宝贵，电力系统时滞稳定裕度求解方法[J]，电力系统自动化，2007，31(2)：5-11.Jia Hongjie，Shang Rui and Zhang Baogui，Computation of delay stability margin of electric powersystems[J]，Automation of Electric Power Systems，2007，31(5)：5-11.

13.Cheng Yicheng and Hwang Chyi，Use of the Lambert W function for time-domain analysis of feedbackfractional delay systems[J]，IEE Proceedings-Control Theory and Applications，2006，153(2)：167-174.

14.Jia hongjie and Yu Xiaodan，A simple method for power system stability analysis with multiple timedelays[C]，Proc.of IEEE 2008 General Meeting，20-24July，Pittsburgh，PA USA.pp.4596157

15.贾宏杰和余晓丹，两种实际约束下的电力系统时滞稳定裕度[J]，电力系统自动化，2008，32(9)：7-10，19.Jia Hongjie and Yu Xiaodan，Method of determining power system delay margins with considering twopractic al constraints[J]，Automation of Electric Power Systems，2008，32(9)：7-10，19.

16.J.K.Hale，Theory of Functional Differential Equations，New York，Springer-Verlag，1977.

发明内容

为解决这一问题，本发明提供了一种简便方法，仅需要在有限范围内追踪一个复矩阵的特征值轨迹，即可求解系统在大范围时滞空间中的稳定域及其边界，计算效率较高。为此，本发明采用如下的技术方案：

一种时滞电力系统稳定的判别方法，其特征在于，包括下列步骤：

第一步：将含有时滞环节的时滞电力系统在平衡点(x₀，y₀)处进行线性化：

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{x}=A_0\mathrm{\Δx}+A_1\mathrm{\Δ}{x}_{\mathrm{\τ}1}+A_2\mathrm{\Δ}{x}_{\mathrm{\τ}2}+\·\·\·+A_k\mathrm{\Δ}{x}_{\mathrm{\τk}}\\ +B_0\mathrm{\Δy}+B_1\mathrm{\Δ}{y}_{\mathrm{\τ}1}+B_2\mathrm{\Δ}{y}_{\mathrm{\τ}2}+\·\·\·B_k{\mathrm{\Δ}}_{\mathrm{\τk}}\\ 0=C_0\mathrm{\Δx}+D_0\mathrm{\Δy}\\ \begin{array}{cc}0=C_i\mathrm{\Δ}{x}_{\mathrm{\τi}}+D_i\mathrm{\Δ}{y}_{\mathrm{\τi}}& i=\mathrm{1,2},\·\·\·k\end{array}\end{array}\right.$

式中：x∈Rⁿ，y∈R^m分别为系统的状态变量和代数变量；Δ表示各变量的增量；设(x_τi，y_τi)＝(x(t-τ_i)，y(t-τ_i))为系统时滞状态变量和时滞代数变量；τ₁，τ₂，…，τ_k为时滞常数；且有

$A_0={\frac{\∂f}{\∂x}|}_p,B_0={\frac{\∂f}{\∂y}|}_p,C_0={\frac{\∂g}{\∂x}|}_p,D_0={\frac{\∂g}{\∂y}|}_p$

$A_i={\frac{\∂f}{{\∂x}_{\mathrm{\τi}}}|}_p,B_i={\frac{\∂f}{{\∂y}_{\mathrm{\τi}}}|}_p,C_i={\frac{\∂g}{{\∂x}_{\mathrm{\τi}}}|}_p,D_i={\frac{\∂g}{{\∂y}_{\mathrm{\τi}}}|}_p$

分别为系统方程对状态变量和时滞变量的导数，

当矩阵D₀，D₁，D₂…，D_k非奇异时，方程(1)可简化为：

$\mathrm{\Δ}\stackrel{\·}{x}={\tilde{A}}_0\mathrm{\Δx}+{\tilde{A}}_1\mathrm{\Δ}{x}_{\mathrm{\τ}1}+{\tilde{A}}_2\mathrm{\Δ}{x}_{\mathrm{\τ}2}+\·\·\·+{\tilde{A}}_k\mathrm{\Δ}{x}_{\mathrm{\τk}}$

${\tilde{A}}_i=A_i-B_i{D}_i^{-1}{C}_i$ i＝0，1…，k

给定实际电力系统中的真实时滞长度(τ₀，τ₁，……τ_m)，

并令：ε_τ＝(ε₁，ε₂，...，ε_k)∈R^k，ε_i＞0，λ^c＝j·ω^c，ω^c＞0，

其中λ^c代表实部为零的系统特征值，ω^c代表改特征值的虚部，则时滞系统的特征方程可表示为： $\mathrm{\Γ}\left({\mathrm{\ω}}^c\right)=\det (\mathrm{\λI}-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-j{\mathrm{\ϵ}}_i^c})=0;$

第二步：设

并定义函数R_p(·)用于表示矩阵M(ε_τ)拥有正实部特征值的个数：m＝R_p(M(ε_τ))；

第三步：将每一个参数ε_i对应的搜索空间[0，2π)进行n等分，确定搜索步长h＝2π/n；

第四步：，其中，循环进行M(ε_τ)矩阵特征值的跟踪，并确定系统时滞空间的临界点，方法如下：设ε₁，ε₂，…，ε_k为用于循环的变量，分别用于循环的第1，2，…，k层，且都从2π/n以步长h递增至2π，在循环体内部顺序执行如下四个步骤的操作：①ε_τ＝(ε₁，ε₂，…，ε_k)；②计算此时M(ε_r)的特征值α(ε_τ)；③计算此时M(ε_r)正实部特征值数m；④若存在相邻特征值穿越虚轴则计算λ^c和对应的

第五步：进一步，通过下式即得到在时滞τ_i方向上(-∞，+∞)范围内所有与λ^c相关的系统稳定临界点集合：

r_i的最大取值由下式来确定：

函数R_d(·)为取整函数；

第六步：并利用下式决定时滞稳定区域的边界曲线：

i＝1，2，…，k；

第七步：令时滞τ_i微增Δτ_i，i＝1，2，...，k，利用如下插值公式求解此时对应的系统关键特征值：

y₁＝λ^c

$y_2={\mathrm{\λ}}^c+d\·\exp \left(\frac{2\mathrm{\π}}{3}j\right)$

$y_3={\mathrm{\λ}}^c+d\·\exp (-\frac{2\mathrm{\π}}{3}j)$

y₄＝λ^c+d

$x_m=\det (y_m\·I-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-y_m({\mathrm{\τ}}_i+\mathrm{\Δ}{\mathrm{\τ}}_i)}),$ m＝1，2，3，4其中：d的取值为一个微小的增量，由拉格朗日插值公式，求得微扰后系统的关键特征值为： ${\mathrm{\λ}}_i^{+}=\underset{m=1}{\overset{4}{\mathrm{\Σ}}}{y}_m\underset{j\≠m}{\underset{1\≤j\≤4}{\mathrm{\Π}}}\frac{-x_j}{x_m-x_j};$

第八步：令τ_i微减Δτ_i，利用第七步的插值公式及

m＝1，2，3，4，求解此时的关键特征值λ_i ^-；

第九步：利用λ^c、λ_i ⁺和λ_i ^-在复平面的分布情况，可得系统关键特征值沿τ_i方向的变化规律；

第十步：对每一维时滞参数τ_i，i＝1，2，3，...，k，均采用上述过程判断对应的系统关键特征值随τ_i的变化规律，并采用下列原则精确判断系统沿各个方向稳定区域与不稳定区域的分布以及其边界的相关性质：如果关键特征值越过了虚轴来到了左半平面，且此时系统已经没有位于右半平面的特征值，则系统是稳定的；如果关键特征值越过了虚轴来到了右半平面，则系统是不稳定的；

第十一步：判断时滞长度数据(τ₁，……τ_k)是否位于稳定域内，若位于，则判定(τ₁，……τ_k)在k维空间的稳定区域中，否则，则判定(τ₁，……τ_m)不在k维空间的稳定区域中。

针对传统时滞稳定域分析方法，只能适用于时滞较小情况的不足，本发明给出的时滞电力系统稳定的判别方法，是一种适用于大范围时滞稳定域分析的有效方法。它首先利用参数变换，通过对一边长为2π的k-维空间的超正方体进行临界点搜索；然后通过平移和反变换，实现电力系统大范围时滞稳定域边界的求解；最后采用四点插值法，通过对稳定域边界点实施微扰来判断稳定域的构成和边界性质。

附图说明

图1临界点搜索原理。

图2验证系统在(ε₁，ε₂)空间中的稳定域边界。

图3验证系统在(τ₁，τ₂)空间中的稳定域边界。

图4验证系统大范围时滞空间中的稳定域。

图5WSCC-3机9节点系统在(ε₁，ε₂)空间的稳定域边界。

图6WSCC-3机9节点系统在(τ₁，τ₂)空间的稳定域边界。

图7WSCC-3机9节点系统在(τ₁，τ₂)空间的稳定域

图8WSCC-3机9节点系统时滞较小时的稳定区域。

具体实施方式

下面从电力系统时滞模型、本发明所依据的稳定判据及其证明、本发明的时滞系统判稳方法以及实施例几个方面对本发明做进一步详述。

1.电力系统时滞模型

将含有时滞环节的时滞电力系统在平衡点(x₀，y₀)处进行线性化：

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{x}=A_0\mathrm{\Δx}+A_1\mathrm{\Δ}{x}_{\mathrm{\τ}1}+A_2\mathrm{\Δ}{x}_{\mathrm{\τ}2}+\·\·\·+A_k\mathrm{\Δ}{x}_{\mathrm{\τk}}\\ +B_0\mathrm{\Δy}+B_1\mathrm{\Δ}{y}_{\mathrm{\τ}1}+B_2\mathrm{\Δ}{y}_{\mathrm{\τ}2}+\·\·\·B_k{\mathrm{\Δ}}_{\mathrm{\τk}}\\ 0=C_0\mathrm{\Δx}+D_0\mathrm{\Δy}\\ \begin{array}{cc}0=C_i\mathrm{\Δ}{x}_{\mathrm{\τi}}+D_i\mathrm{\Δ}{y}_{\mathrm{\τi}}& i=\mathrm{1,2},\·\·\·k\end{array}\end{array}\right.---\left(1\right)$

式中：x∈Rⁿ，y∈R^m分别为系统的状态变量和代数变量；Δ表示个变量的增量；(x_τi，y_τi)＝(x(t-τ_i)，y(t-τ_i))为系统时滞状态变量和时滞代数变量；τ₁，τ₂，…，τ_k为时滞常数；且有

$A_0={\frac{\∂f}{\∂x}|}_p,B_0={\frac{\∂f}{\∂y}|}_p,C_0={\frac{\∂g}{\∂x}|}_p,D_0={\frac{\∂g}{\∂y}|}_p---\left(2\right)$

$A_i={\frac{\∂f}{{\∂x}_{\mathrm{\τi}}}|}_p,B_i={\frac{\∂f}{{\∂y}_{\mathrm{\τi}}}|}_p,C_i={\frac{\∂g}{{\∂x}_{\mathrm{\τi}}}|}_p,D_i={\frac{\∂g}{{\∂y}_{\mathrm{\τi}}}|}_p---\left(3\right)$

分别为系统方程对状态变量和时滞变量的导数。

当矩阵D₀，D₁，D₂…，D_k非奇异时，方程(1)可简化为：

$\mathrm{\Δ}\stackrel{\·}{x}={\tilde{A}}_0\mathrm{\Δx}+{\tilde{A}}_1\mathrm{\Δ}{x}_{\mathrm{\τ}1}+{\tilde{A}}_2\mathrm{\Δ}{x}_{\mathrm{\τ}2}+\·\·\·+{\tilde{A}}_k\mathrm{\Δ}{x}_{\mathrm{\τk}}---\left(4\right)$

${\tilde{A}}_i=A_i-B_i{D}_i^{-1}{C}_i$ i＝0，1…，k    (5)

进一步地，时滞系统的特征方程可表示为：

$\mathrm{\Γ}\left(\mathrm{\λ}\right)=\det (\mathrm{\λI}-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-\mathrm{\λ}{\mathrm{\τ}}_i})=0.---\left(6\right)$

2本发明的时滞电力系统稳定的判别方法

给定实际电力系统中的真实时滞长度(τ0，τ1，……τm)，按照下列方法求解该系统大范围时滞稳定域。

a.时滞稳定域边界求解的预备工作

令：

$M\left({\mathrm{\ϵ}}_{\mathrm{\τ}}\right)={\tilde{A}}_0+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-j\·{\mathrm{\ϵ}}_i}---\left(7\right)$

其中：ε_τ＝(ε₁，ε₂，...，ε_k)∈R^k，ε_i＞0。

令：

λ^c＝j·ω^c，ω^c＞0 ${\mathrm{\ϵ}}_i^c={\mathrm{\ω}}^c\·{\mathrm{\τ}}_i---\left(8\right)$

其中λ^c代表实部为零的系统特征值，ω^c代表改特征值的虚部，则式(6)此时将变为：

$\mathrm{\Γ}\left({\mathrm{\ω}}^c\right)=\det (\mathrm{\λI}-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-j{\mathrm{\ω}}^c\·{\mathrm{\τ}}_i})=\det (\mathrm{\λI}-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-j{\mathrm{\ϵ}}_i^c})=0---\left(9\right)$

进一步，通过下式即可得到在时滞τ_i方向上(-∞，+∞)范围内所有与τ^c相关的系统稳定临界点集合：

$T_i^c=(\frac{{\mathrm{\ϵ}}_i^c}{{\mathrm{\ω}}^c},\frac{{\mathrm{\ϵ}}_i^c\±r_i\·2\mathrm{\π}}{{\mathrm{\ω}}^c})---\left(10\right)$

在实际系统中，考虑到τ_i取值大于0，则式(10)中只取正号即可。对于时滞τ_i，当我们关心其范围为[0，τ_i ^max]时，在利用式(10)式确定系统的临界点时，r_i的最大取值由下式来确定：

$r_{i,\max }=R_d\left(\frac{{\mathrm{\τ}}_i^{\max }\·{\mathrm{\ω}}^c-{\mathrm{\ϵ}}_i^c}{2\mathrm{\π}}\right)+1---\left(11\right)$

函数R_d(·)为取整函数。

b.时滞稳定域边界快速求解

第一步：利用式(7)和(8)，对原系统特征方程(6)进行变换，形成矩阵M(ε_τ)。

为便于描述，定义函数R_p(·)用于表示矩阵M(ε_τ)拥有正实部特征值的个数：

m＝R_p(M(ε_τ))    (12)

第二步：将每一个参数ε_i对应的搜索空间[0，2π)进行n等分，确定搜索步长：

h＝2π/n    (13)

第三步：循环进行M(ε_τ)矩阵特征值的跟踪，并确定系统时滞空间的临界点。设ε₁，ε₂，…，ε_k为用于循环的变量，分别用于循环的第1，2，…，k层，且都从2π/n以步长h递增至2π，在循环体内部执行如下操作：

ε_τ＝(ε₁，ε₂，…，ε_k)

计算此时M(ε_r)的特征值α(ε_τ)

计算此时M(ε_r)正实部特征值数m

若存在相邻特征值穿越虚轴则计算λ^c和对应的

第四步：根据所关心的时滞取值范围，利用(11)式计算各个时滞方向上的r_i的最大取值；进一步，利用(10)式进行临界点的平移和反变换。

第五步：输出计算结果，并利用下式决定时滞稳定区域的边界曲线：

${\mathrm{\τ}}_i^c=\frac{{\mathrm{\ϵ}}_i^c}{{\mathrm{\λ}}^c},$ i＝1，2，…，k    (14)

c.稳定区域判别

第一步：令时滞τ_i微增Δτ_i，i＝1，2，…，k，利用如下插值公式求解此时对应的系统关键特征值：

y₁＝λ^c    (15)

$y_2={\mathrm{\λ}}^c+d\·\exp \left(\frac{2\mathrm{\π}}{3}j\right)---\left(16\right)$

$y_3={\mathrm{\λ}}^c+d\·\exp (-\frac{2\mathrm{\π}}{3}j)---\left(17\right)$

y₄＝λ^c+d    (18)

$x_m=\det (y_m\·I-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-y_m({\mathrm{\τ}}_i+\mathrm{\Δ}{\mathrm{\τ}}_i)}),m=\mathrm{1,2,3,4}---\left(19\right)$

其中：d为一微小增量，取值越小越好，本发明中取0.001即可。由拉格朗日插值公式，可求得微扰后系统的关键特征值为：

${\mathrm{\λ}}_i^{+}=\underset{m=1}{\overset{4}{\mathrm{\Σ}}}{y}_m\underset{j\≠m}{\underset{1\≤j\≤4}{\mathrm{\Π}}}\frac{-x_j}{x_m-x_j}---\left(20\right)$

第二步：令τ_i微减Δτ_i，同样利用(15)-(20)求解此时的关键特征值λ_i ^-，此时式(19)将改为：

$x_m=\det (y_m\·I-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-y_m({\mathrm{\τ}}_i-\mathrm{\Δ}{\mathrm{\τ}}_i)})---\left(21\right)$

m＝1，2，3，4

第三步：利用λ^c、λ_i ⁺和λ_i ^-在复平面的分布情况，可得系统关键特征值沿τ_i方向的变化规律。

第四步：对每一维时滞参数τ_i，i＝1，2，3，...，k，均采用上述过程判断对应的系统关键特征值随τ_i的变化规律，并采用下列原则精确判断系统沿各个方向稳定区域与不稳定区域的分布以及其边界的相关性质：

如果关键特征值越过了虚轴来到了左半平面，且此时系统已经没有位于右半平面的特征值，则系统是稳定的；如果关键特征值越过了虚轴来到了由半平面，则系统是不稳定的。

判断时滞长度数据(τ₁，……τ_k)是否位于稳定域内，若位于，则判定(τ₁，……τ_k)在k维空间的稳定区域中，否则，则判定(τ₁，……τ_m)不在k维空间的稳定区域中。

3算例分析

首先利用[11]所提供的一典型时滞系统验证本发明方法的有效性，进而利用该方法分析WSCC-3机9节点系统大范围时滞稳定域的规律。

典型算例

对于[11]所提供的验证算例，经推导可得该系统的特征方程如下：

$A+B{e}^{-{\mathrm{\τ}}_1\mathrm{\λ}}+C{e}^{-{\mathrm{\τ}}_2\mathrm{\λ}}+8e^{-({\mathrm{\τ}}_1+{\mathrm{\τ}}_2)\mathrm{\λ}}=0---\left(22a\right)$

其中：

A＝s²+7.1s+21.1425    (22b)

B＝6s+14.8            (22c)

C＝2s+7.3             (22d)

●时滞参数空间稳定域边界的求解

利用式(10)和(11)，对式(22)进行参数变换，结果如下：

$A+B{e}^{-{\mathrm{\ϵ}}_{1}j}+C{e}^{-{\mathrm{\ϵ}}_{2}j}+8e^{-({\mathrm{\ϵ}}_1+{\mathrm{\ϵ}}_2)j}=0---\left(23\right)$

利用本发明所提供方法，在(ε₁，ε₂)空间中的[0，2π)范围内计算系统临界点，所得结果示于图2(对应于图中阴影部分)。可以看到，系统的时滞稳定域边界(临界点集合)在(ε₁，ε₂)空间中构成一条封闭曲线。

假设我们关心的时滞参数范围为：0≤τ₁≤5.0，0≤τ₂≤5.0，即则利用式(11)，可求得：

r_1，max＝r_2，max＝5    (24)

注意，此时在利用式(11)时，应取稳定域边界上所有临界点中虚部ω^c最大的参与计算。

利用式(10)和式(2)的结果，将(ε₁，ε₂)空间中稳定域边界的计算结果进行平移，然后反变换到(τ₁，τ₂)空间，所得结果示于图3。由变换式(14)可知，(ε₁，ε₂)空间中的一条封闭曲线变换到(τ₁，τ₂)空间后仍将是一条封闭曲线。考虑到稳定域边界上不同临界点对应的ω^c取值不同，图2中的封闭曲线(临界点集合)在变换到(τ₁，τ₂)空间后，将变得不规则。

●稳定区域的判别

由图2和图3可知，无论是在(ε₁，ε₂)空间还是在(τ₁，τ₂)空间，系统稳定域的边界都由一族封闭曲线共同构成。对于每一条封闭曲线，其内外区域的稳定性刚好相反；同时任意两条封闭曲线均可利用式(10)通过平移后实现互换，因此，利用任意一条封闭曲线判断其内外区域的稳定性，所得结果可推广应用到其它所有的边界封闭曲线。

本发明以图3中左下角封闭曲线为对象，取其上一点P加以研究，该点坐标：

对应的关键特征值为：λ^c＝±j3.3224：

1)对P点沿τ₁方向进行微扰分析

令τ₁ ^c增加Δτ＝d＝0.001(下同)，利用本发明的方法，可求得此时的关键特征值为：λ⁺＝0.0006±j3.2194；进一步，令τ₁ ^c减少Δτ，可求得此时的关键特征值为：λ^-＝-0.0006±j3.3255。

2)对P点沿τ₂方向进行微扰分析

令τ₂ ^c增加Δτ，利用本发明的方法，求得此时的关键特征值为：λ⁺＝-0.0003±j3.3224；进一步，令τ₂ ^c减少Δτ，对应的关键特征值变为：λ^-＝0.0004±j3.3225。

从上面分析可知，该封闭曲线的内部为不稳定区域，外部为稳定区域。将此结果在(τ₁，τ₂)空间中的其他封闭曲线进行推广，可得系统在0≤τ₁≤5.0和0≤τ₂≤5.0范围内稳定域的计算结果，如图4所示。图中灰色区域为不稳定区域，白色区域为稳定区域，所得结果与[11]中的结果完全相同。由于[11]方法需要在(-∞，+∞)范围内搜索系统所有的临界点，计算量巨大，与之相比，本发明仅需在一个边长为2π的正方形空间中进行临界点搜索，因此具有更高的计算效率。

WSCC-3机9节点算例

WSCC-3机9节点系统，考虑发电机G2和G3的控制回路中同时存在时滞，其中G2回路时滞记为τ₁，G3回路时滞记为τ₂。系统模型推导参见[14]，分析时，取负荷水平2.0p.u，P_m2＝P_m3＝1.0p.u，V_ref2＝V_ref3＝1.03p.u.，系统其他参数取值同[14]。

在平衡点处，系统的时滞模型可表示为式(6)所示形式，系数矩阵

的取值见附录。假设我们关心的时滞参数取值范围为：0≤τ₁≤7.5和0≤τ₂≤7.5，下面利用本发明方法进行求解。

●时滞参数空间稳定域边界的求解

利用式(7)和(9)，对WSCC-3机9节点系统时滞方程进行变换可得：

进一步，利用本发明所给方法，可求得在(ε₁，ε₂)空间中的稳定域边界，同时利用式(11)，根据所关心的时滞参数范围，可求得：

r_1，max＝r_2，max＝6    (25)

通过平移操作后，所得最终结果示于图5。从中不难看出，WSCC-3机9节点系统在(ε₁，ε₂)空间中的稳定域边界也由周期性的封闭曲线构成，但靠近ε₁和ε₂轴的封闭曲线，均与两坐标轴相交。进一步，通过式(14)进行反变换可得系统在(τ₁，τ₂)空间中的稳定域边界，结果示于图6。

●稳定区域的判别

与上述的典型算例类似，只需利用2.3节方法，对图6中任一封闭曲线上的任何一点进行微扰分析，即可判断系统的稳定区域。本发明利用图6中的Q点加以判断，其分析过程与3.1节类似，为简单起见在此略去。所得结果如图7所示，其中灰色区域为不稳定区域，白色区域则为稳定区域。对图7结果进行深入分析可发现如下规律：

1)将图7中左下角邻近原点的稳定区域(对应于时滞参数取值较小情况)放大后示于图8，通过对比分析可知，其与[12，14]的结果完全相同。但[12，14]方法仅适用于时滞较小情况，而本发明方法则在大范围时滞空间中适用，且计算效率更高。

2)WSCC-3机9节点系统的稳定域由一些互不连通的稳定分区构成。电力系统传统控制设备的量测时滞较小，因此只需关注图8所示的小时滞稳定分区即可。但当量测环节出现阻滞，导致量测时滞变大时，则需要保证控制设备运行在图7所示的稳定区域范围之内。

3)随着时滞τ₁，τ₂数值的增大，WSCC-3机9节点系统的稳定分区不断变小直至完全消失，这表明控制回路时滞是造成系统不稳定的一个重要原因，时滞越大对系统的稳定运行越不利。

下面是WSCC-3机9节点系统算例数据。

${\tilde{A}}_0=\left[\begin{array}{cccccccccc}0& 377& 0& 0& 0& 0& 0& 0& 0& 0\\ -0.1492& -0.0039& -0.0337& -0.1127& 0& 0.1006& 0& 0.1116& 0.0607& 0\\ -0.0231& 0& -0.2559& 0.0516& 0.1667& 0.1471& 0& 0.4597& 0.0130& 0\\ -1.8676& 0& 0.2386& -5.0160& 0& 0.9221& 0& 0.2634& 0.7509& 0\\ 0& 0& -2308.6264& 959.2936& -50.0000& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 377& 0& 0& 0\\ 0.2148& 0& 0.1946& 0.1231& 0& -0.3616& -0.0083& -0.0870& -0.2987& 0\\ 0.1345& 0& 0.3502& 0.0159& 0& -0.0141& 0& -0.1416& 0.0259& 0.1250\\ 2.4094& 0& 0.4318& 1.8511& 0& -5.6875& 0& -0.2429& -14.2197& 0\\ 0& 0& 0& 0& 0& 0& 0& -2358.7906& 828.3157& -50.0000\end{array}\right]---\left(F1\right)$

${\tilde{A}}_{\mathrm{\τ}1}=\left[\begin{array}{cccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -2358.3506& 0& 2271.9084& -593.7330& 0& -173.7266& 0& -942.8067& 86.9951& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]---\left(F2\right)$

${\tilde{A}}_{\mathrm{\τ}2}=\left[\begin{array}{cccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -250.7616& 0& -814.1500& 13.5902& 0& -123.8625& 0& 32.4877& -343.2048& 0\end{array}\right]---\left(F3\right)$

Claims (1)

1.一种时滞电力系统稳定的判别方法，其特征在于，包括下列步骤：

第一步：将含有时滞环节的时滞电力系统在平衡点(x₀，y₀)处进行线性化，得到方程(1)：

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{x}=A_0\mathrm{\Δx}+A_1{\mathrm{\Δx}}_{\mathrm{\τ}1}+A_2{\mathrm{\Δx}}_{\mathrm{\τ}2}+...+A_k{\mathrm{\Δx}}_{\mathrm{\τk}}\\ +B_0\mathrm{\Δy}+B_1{\mathrm{\Δy}}_{\mathrm{\τ}1}+B_2{\mathrm{\Δy}}_{\mathrm{\τ}2}+...+B_k{\mathrm{\Δ}}_{\mathrm{\τk}}\\ \begin{array}{cc}0=C_0\mathrm{\Δx}+D_0\mathrm{\Δy}& \\ 0=C_i{\mathrm{\Δx}}_{\mathrm{\τi}}+D_i{\mathrm{\Δy}}_{\mathrm{\τi}}& i=\mathrm{1,2},...k\end{array}\end{array}\right.$

式中：x∈Rⁿ，y∈R^m分别为系统的状态变量和代数变量；Δ表示各变量的增量；

设(x_τi，y_τi)＝(x(t-τ_i)，y(t-τ_i))为系统时滞状态变量和时滞代数变量；τ₁，τ₂，…，τ_k为时滞常数；

且有

$A_0=\frac{\∂f}{\∂x}{|}_p,B_0=\frac{\∂f}{\∂y}{|}_p,$ $C_0=\frac{\∂g}{\∂x}{|}_p,D_0=\frac{\∂g}{\∂y}{|}_p$

$A_i=\frac{\∂f}{\∂x_{\mathrm{\τi}}}{|}_p,B_i=\frac{\∂f}{\∂y_{\mathrm{\τi}}}{|}_p,$ $C_i=\frac{\∂g}{\∂x_{\mathrm{\τi}}}{|}_p,D_i=\frac{\∂g}{\∂y_{\mathrm{\τi}}}{|}_p$

分别为系统方程对状态变量和时滞变量的导数，

当矩阵D₀，D₁，D₂…，D_k非奇异时，方程(1)简化为：

$\mathrm{\Δ}\stackrel{\·}{x}={\tilde{A}}_0\mathrm{\Δx}+{\tilde{A}}_1{\mathrm{\Δx}}_{\mathrm{\τ}1}+{\tilde{A}}_2{\mathrm{\Δx}}_{\mathrm{\τ}2}+...+{\tilde{A}}_k{\mathrm{\Δx}}_{\mathrm{\τk}}$

${\tilde{A}}_0=A_0-B_0{D}_0^{-1}{C}_0$

${\tilde{A}}_i=A_i-B_i{D}_i^{-1}{C}_i$

给定实际电力系统中的真实时滞长度(τ₀，τ₁，……τ_m)，

并令：ε_τ＝(ε₁，ε₂，…，ε_k)∈R^k，ε_i＞0，λ^c＝j·ω^c，ω^c＞0，

其中λ^c代表实部为零的系统特征值，ω^c代表该特征值的虚部，则时滞系统的特征方程表示为： $\mathrm{\Γ}\left({\mathrm{\ω}}^c\right)=\det (\mathrm{\λI}-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-j{\mathrm{\ϵ}}_i^c})=0;$

第二步：设

并定义函数R_p(·)用于表示矩阵M(ε_τ)拥有正实部特征值的个数：m＝R_p(M(ε_τ))；

第三步：将每一个参数ε_i对应的搜索空间[0，2π)进行n等分，确定搜索步长h＝2π/n；

第四步：其中，循环进行M(ε_τ)矩阵特征值的跟踪，并确定系统时滞空间的临界点，方法如下：设ε₁，ε₂，…，ε_k为用于循环的变量，分别用于循环的第1，2，…，k层，且都从2π/n以步长h递增至2π，在循环体内部顺序执行如下四个步骤的操作：①ε_τ＝(ε₁，ε₂，…，ε_k)；②计算此时M(ε_r)的特征值α(ε_τ)；③计算此时M(ε_r)正实部特征值数m；④若存在相邻特征值穿越虚轴则计算λ^c和对应的

第五步：进一步，通过下式即得到在时滞τ_i方向上(-∞，+∞)范围内所有与λ^c相关的系统稳定临界点集合：

r_i的最大取值由下式来确定：

$r_{i,\max }=R_d\left(\frac{{\mathrm{\τ}}_i^{\max }\·{\mathrm{\ω}}^c-{\mathrm{\ϵ}}_i^c}{2\mathrm{\π}}\right)+1,$ 函数R_d(·)为取整函数；

第六步：并利用下式决定时滞稳定区域的边界曲线：

第七步：令时滞τ_i微增Δτ_i，利用如下插值公式求解此时对应的系统关键特征值：

y₁＝λ^c

$y_2={\mathrm{\λ}}^c+d\·\exp \left(\frac{2\mathrm{\π}}{3}j\right)$

$y_3={\mathrm{\λ}}^c+d\·\exp (-\frac{2\mathrm{\π}}{3}j)$

y₄＝λ^c+d

$x_m=\det (y_m\·I-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-y_m({\mathrm{\τ}}_i+{\mathrm{\Δ\τ}}_i)}),m=\mathrm{1,2,3,4}$ 其中：d的取值为一个微小的增量，由拉格朗日插值公式，求得微扰后系统的关键特征值为： ${\mathrm{\λ}}_i^{+}=\underset{m=1}{\overset{4}{\mathrm{\Σ}}}{y}_m\underset{\underset{j\≠m}{1\≤j\≤4}}{\mathrm{\Π}}\frac{-x_j}{x_m-x_j};$

第八步：令τ_i微减Δτ_i，利用第七步的插值公式及 $x_m=\det (y_m\·I-{\tilde{A}}_0-\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\tilde{A}}_i{e}^{-y_m({\mathrm{\τ}}_i-{\mathrm{\Δ\τ}}_i)}),$ m＝1，2，3，4，求解此时的关键特征值

第九步：利用λ^c、

和

在复平面的分布情况，得到系统关键特征值沿τ_i方向的变化规律；

第十步：对每一维时滞参数τ_i，均采用上述过程判断对应的系统关键特征值随τ_i的变化规律，并采用下列原则精确判断系统沿各个方向稳定区域与不稳定区域的分布以及其边界的相关性质：如果关键特征值越过了虚轴来到了左半平面，且此时系统已经没有位于右半平面的特征值，则系统是稳定的；如果关键特征值越过了虚轴来到了右半平面，则系统是不稳定的；

第十一步：判断时滞长度数据(τ₁，……τ_k)是否位于稳定域内，若位于，则判定(τ₁，……τ_k)在k维空间的稳定区域中，否则，则判定(τ₁，……τ_k)不在k维空间的稳定区域中。

Images