CN104680426A - 基于伊藤微分的时滞电力系统随机稳定性分析方法及系统 - Google Patents

Abstract

本发明涉及电力系统稳定分析领域，尤其涉及一种基于伊藤微分的时滞电力系统随机稳定性分析方法及系统，包括：首先，建立考虑随机因素的Lyapunov-Krasovskii目标泛函，并借助伊藤微分公式求解目标泛函的弱无穷小生成算子，然后，在该算子中加入由Newton-Leibniz公式构造的自由权项以降低保守性，在此基础上，建立一组线性矩阵不等式，并求解考虑不同随机激励强度的电力系统在过负荷、弱阻尼以及高励磁倍数等恶劣工况下，所能承受的最大时滞。IEEE16机68节点系统的仿真结果验证了本发明的有效性及可行性，同时保守性分析结果表明本发明具有较低的保守性。

Description

基于伊藤微分的时滞电力系统随机稳定性分析方法及系统

技术领域

本发明涉及电力系统稳定分析领域，尤其涉及一种基于伊藤微分的时滞电力系统随机稳定性分析方法及系统。

背景技术

广域反馈控制器为电力系统稳定分析带来新的契机，但同时也引入时滞的问题，而时滞的存在必然会削弱广域控制器的性能，以致出现负阻尼的情况，因此，迫切需要对系统的时滞稳定性进行深入研究。

在时滞稳定性分析方面，已有不少理论成果，其稳定性结论大致可分为两大类：1)时滞独立稳定性判据。若一个稳定判据中不含时滞h_t，则该稳定判据称为时滞独立稳定性判据。时滞独立稳定性判据主要利用Rekasius变换处理时滞系统的特征方程，以求解系统时滞稳定上限，但该方法需要在时滞空间中搜寻系统的关键特征值，计算量较大，很难适用于大规模系统的计算。时滞独立稳定性判据允许系统的时滞不确定或未知，但是，当时滞较小时，时滞独立稳定条件将比时滞相关稳定条件的保守性更高。2)时滞依赖稳定性。若一个稳定判据包含h_t，则该稳定判据称为时滞依赖稳定判据。时滞依赖稳定判据主要基于Lyapunov-Krasovskii泛函。该类判据充分利用了时滞变化等信息，能够考虑系统存在不确定性以及时变时滞等情况，具有较低的保守性。然而，上述成果均未考虑接入电力系统的随机因素，随着新能源发电的快速发展以及电动汽车等各种随机负荷的不断投入使用，电力系统将无法利用确定性的微分方程理论分析系统的时滞稳定性。

例如文献[1]：D.Yue and Q.-L.Han.Delay-dependent exponential stability ofstochastic systems with time-varying delay,nonlinearity,and Markovian switehing.IEEE Automatic Control,2005,50(2):217-222和文献[2]：陈云.随机时滞系统的分析与综合[D].杭州：浙江大学，2008。文献[1]研究了随机系统在考虑时变时滞、马尔可夫切换以及非线性时均方意义下的时滞依赖稳定性问题；文献[2]讨论了具有时变时滞的非线性系统的鲁棒指数稳定性，然后推广到随机时滞系统中讨论其时滞相关指数稳定性问题。文献[1]中的方法可以分析含有非线性以及马尔科夫切换的随机时滞稳定性问题，文献[2]通过定义一个辅助的向量并引入适当的自由权矩阵，无需考虑矩阵不等式约束条件，但是文献[1]和文献[2]均存在保守性较大的问题。

发明内容

为了解决传统时滞电力系统分析未考虑接入电力系统的随机因素的问题，本发明提出了一种基于伊藤微分的时滞电力系统随机稳定性分析方法，包括：

步骤1、采集时滞电力系统网络结构参数；

步骤2、将采集到的时滞电力系统网络结构参数构建得到系统的状态矩阵和时滞矩阵，并与系统的状态向量一起根据时滞系统方程来构造考虑随机因素的Lyapunov-Krasovskii泛函，并借助伊藤微分公式求解目标泛函的弱无穷小算子；

步骤3、利用Newton-Leibniz公式建立与自由权矩阵L和M有关的等式，并将该等式与必要的松散项加入到Lyapunov-Krasovskii泛函的弱无穷小算子之中，得到时滞系统稳定判据；

步骤4、将时滞系统稳定判据等价变换成符合广义特征值最优法求解的标准形式，并将经Schur降阶后的状态矩阵和时滞矩阵以及统计得到的随机激励强度值带入标准形式，求解得到随机时滞系统的时滞稳定上限。

所述步骤2中的Lyapunov-Krasovskii泛函为：

$\begin{array}{c}V(x_t,t)=x^T\left(t\right)\mathrm{Px}\left(t\right)+{\∫}_{t-h_t}^t{x}^T\left(s\right)\mathrm{Qx}\left(s\right)\mathrm{ds}\\ +{\∫}_{t-\stackrel{\‾}{h}}^t{x}^T\left(s\right)\mathrm{Kx}\left(s\right)\mathrm{ds}\\ +{\∫}_{-\stackrel{\‾}{h}}^0{\∫}_{t+\mathrm{\θ}}^t{y}^T\left(s\right)\mathrm{Ry}\left(s\right)\mathrm{dsd\θ}\\ +{\∫}_{-\stackrel{\‾}{h}}^0{\∫}_{t+\mathrm{\θ}}^t{f}^T\left(s\right)\mathrm{Zf}\left(s\right)\mathrm{dsd\θ}\end{array}---\left(3\right)$

其中，x(t)∈Rⁿ是系统的状态向量，时滞h_t满足为时滞稳定上限，为时滞h_t的一阶导数，μ为时滞的最大变化率，向量城y(t)∈R_n满足：y(t)dt＝dx(t),P、Q、R、K、Z为具有适当维数的正定矩阵，是待求变量，为时滞稳定上限，f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵；s,θ为替换变量t的多重积分变量。

所述步骤2中Lyapunov-Krasovskii泛函的弱无穷小算子为：

$\begin{array}{c}\mathrm{\ζV}(t,x\left(t\right))=2x^T\left(t\right)\mathrm{Pf}\left(t\right)+x^T\left(t\right)(Q+K)x\left(t\right)\\ -(1-{\stackrel{\·}{h}}_t)x^T(t-h_t)\mathrm{Qx}(t-h_t)\\ -x^T(t-\stackrel{\‾}{h})\mathrm{Kx}(x-\stackrel{\‾}{h})\\ +\stackrel{\‾}{h}{y}^T\left(t\right)\mathrm{Ry}\left(t\right)-{\∫}_{t-\stackrel{\‾}{h}}^t{y}^T\left(t\right)\mathrm{Ry}\left(t\right)\mathrm{ds}\\ +\stackrel{\‾}{h}{f}^T\left(t\right)\mathrm{Zf}\left(t\right)-{\∫}_{t-\stackrel{\‾}{h}}^t{f}^T\left(t\right)\mathrm{Zf}\left(t\right)\mathrm{ds}\\ +\mathrm{tr}\left\{g^T\left(x\left(t\right)\right)\mathrm{Pg}\left(x\left(t\right)\right)\right\}\end{array}---\left(5\right)$

其中，x(t)∈Rⁿ是系统的状态向量，时滞h_t满足为时滞稳定上限，为时滞h_t的一阶导数，μ为时滞的最大变化率，向量城y(t)∈R_n满足：y(t)dt＝dx(t),P、Q、R、K、Z为具有适当维数的正定矩阵，是待求变量，f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵；s,θ为替换变量t的多重积分变量。

所述步骤3中改进自由权矩阵L和M有关的等式为：

$2{\mathrm{\ξ}}^T\left(t\right)L[x\left(t\right)-x(t-\stackrel{\‾}{h})-{\∫}_{t-\stackrel{\‾}{h}}^t\mathrm{dx}\left(t\right)]=0$

2y^T(t)M[(Ax+A_dx_t-ht-y(t))dt+g(x(t))dω(t)]＝0

其中， $\mathrm{\ξ}\left(t\right)={\left[\begin{array}{ccccc}{x}^T& x^T(t-h_t)& x^T(t-\stackrel{\‾}{h})& y^T\left(t\right)& f^T\left(t\right)]\end{array}\right]}^T,$ x(t)∈Rⁿ是状态向量，向量城y(t)∈R_n满足：y(t)dt＝dx(t),为时滞稳定上限，μ为时滞的最大变化率，时滞h_t满足时滞d(t)与一阶导数满足条件：f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵。

所述步骤3中必要的松散项为：

2ζ^T(t)S[Ax(t)-A_dx(t-h_t)-f(t)]＝0

其中，f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵，h_t为时滞。

所述步骤3中时滞系统随机稳定判据如下：

考虑随机因素的时滞电力系统模型表示为：

其中，x(t)∈Rⁿ是系统的状态向量，时滞h_t满足矩阵A,A_d为已知实矩阵，W(t)是定义在完备概率空间(Ω,F,P)上的Wiener过程，且满足ε(dW(t))＝0,ε(dW²(t))＝dt；

定义向量城y(t)∈R_n，使其满足：

对于式(1)所示的随机时滞电力系统，给定标量h>0和μ，若存在P＝P^T>0，Q＝Q^T≥0，R＝R^T≥0，K＝K^T>0，Z＝Z^T≥0，以及适维矩阵L、M和S，使得如下线性矩阵不等式(14)成立，则随机时滞系统(1)均方稳定；

$\left[\begin{array}{cc}\mathrm{\&Omega;}& L\\ *& -T_1-Z/\stackrel{\&OverBar;}{h}\end{array}\right]<0---\left(14\right)$

其中：

$\mathrm{\&Omega;}={\mathrm{\&Omega;}}_1+{\mathrm{\&Omega;}}_2+{\mathrm{\&Omega;}}_2^T$

Ω₂＝[L+SA SA_d -L 0 -S]

${\mathrm{\&Omega;}}_1=\left[\begin{array}{ccccc}\begin{array}{c}Q+K+G_1^{T}P{G}_1\\ -R/\stackrel{\&OverBar;}{h}\end{array}& R/\stackrel{\&OverBar;}{h}& 0& A^{T}M& P\\ *& \begin{array}{c}{G}_2^{T}P{G}_2-(1-\mathrm{\&mu;})Q\\ -R/\stackrel{\&OverBar;}{h}\end{array}& 0& A_d^{T}M& 0\\ *& *& -K& 0& 0\\ *& *& *& \stackrel{\&OverBar;}{h}R-M-M^T& 0\\ *& *& *& *& \stackrel{\&OverBar;}{h}Z\end{array}\right]$

其中，P、Q、R、K、Z、T₁为具有适当维数的正定矩阵，是待求变量；L、M、S为具有适当维数的实矩阵，是待求变量；Ω、Ω₁、Ω₂为具有适当维数的矩阵，G₁、G₂为已知变量。

所述时滞电力系统网络结构参数包括：线路阻抗值、母线电压幅值和相位、发电机有功与无功出力、有功负荷与无功负荷、发电机电抗后电势与发电机励磁电势、发电机阻尼系数、发电机定子时间常数与发电机励磁回路时间常数、发电机励磁回路放大系数、发电机励磁电势与发电机机端电压参考值、发电机稳态电抗、发电机暂态电抗、发电机次暂态电抗、发电机功角与转速。

所述系统状态向量包括：发电机功角、发电机角速度、发电机电抗后电势、励磁系统电压调节器输出电压、励磁系统输出电压、励磁系统励磁反馈电压、原动机与调速器的ΔP_m,Δμ,Δx_m，m＝1表示为汽轮机，m＝2表示为水轮机、电力系统稳定器的Δy₁,Δy₂,Δy₃，其中，ΔP_m为机械功率增量，Δμ为汽门开度增量，Δx₁为硬反馈量增量，Δx₂为总反馈量增量；对于PSS系统,Δy₁为复位状态变量增量，Δy₂为第一次超前滞后补偿状态变量增量，Δy₃为第二次超前滞后补偿状态变量增量。

所述随机激励强度值通过统计得到，类别包括：负荷的随机波动、原动机扭矩的随机振动、控制回路的测量噪声、互联大电网中功率角的随机小幅振荡、风力发电中风速的随机变化。

一种基于伊藤微分的时滞电力系统随机稳定性分析系统，包括：数据采集模块、时滞上限求解模块和结果输出模块；其中，数据采集模块通过时滞上限求解模块与结果输出模块相连

所述数据采集模块用于采集时滞电力系统网络结构参数，并将采集数据发送至时滞上限求解模块；

所述时滞上限求解模块的功能是：

首先将采集到的时滞电力系统网络结构参数构建得到系统的状态矩阵和时滞矩阵，并与系统的状态向量一起根据时滞系统方程来构造考虑随机因素的Lyapunov-Krasovskii泛函，并借助伊藤微分公式求解目标泛函的弱无穷小算子；

然后利用Newton-Leibniz公式建立与自由权矩阵L和M有关的等式，并将该等式与必要的松散项加入到Lyapunov-Krasovskii泛函的弱无穷小算子之中，得到时滞系统稳定判据；

最后将时滞系统稳定判据等价变换成符合广义特征值最优法求解的标准形式，并将经Schur降阶后的状态矩阵和时滞矩阵以及统计得到的随机激励强度值带入标准形式，求解得到随机时滞系统的时滞稳定上限；

所述结果输出模块用于输出时滞稳定上限结果。

本发明的有益效果在于：本发明利用改进伊藤微分建立时滞系统随机稳定判据，借助广义特征值法对系统时滞上限进行求解。仿真结果表明，该方法能够有效降低时滞稳定上限计算过程中的保守性，具有很好的正确性和有效性。

附图说明

图1为考虑电力系统随机特性的时滞稳定上限计算系统结构；

图2为IEEE16机68节点系统结构图；

图3为时，16机系统不同时滞下发电机G1-G13相对功角动态响应；

图4为时，16机系统不同时滞下发电机G8-G16相对功角动态响应；

图5为时，16机系统不同时滞下发电机G1-G16相对功角动态响应；

图6为时，16机系统不同时滞下发电机G6-G15相对功角动态响应；

图7为时，16机系统不同时滞下发电机G9-G15相对功角动态响应；

图8为时，16机系统不同时滞下发电机G4-G13相对功角动态响应

具体实施方式

下面结合附图，对优选实施例作详细说明。

如图1所示，一种基于伊藤微分的时滞电力系统随机稳定性分析系统，包括：数据采集模块、时滞上限求解模块和结果输出模块，其中，数据采集模块通过时滞上限求解模块与结果输出模块相连。

所述数据采集模块用于采集时滞电力系统网络结构参数，并将采集数据发送至时滞上限求解模块；

所述时滞上限求解模块的功能是：

首先将采集到的时滞电力系统网络结构参数构建得到系统的状态矩阵和时滞矩阵，并与系统的状态向量一起根据时滞系统方程来构造考虑随机因素的Lyapunov-Krasovskii泛函，并借助伊藤微分公式求解目标泛函的弱无穷小算子；

然后利用Newton-Leibniz公式建立与自由权矩阵L和M有关的等式，并将该等式与必要的松散项加入到Lyapunov-Krasovskii泛函的弱无穷小算子之中，得到时滞系统稳定判据；

最后将时滞系统稳定判据等价变换成符合广义特征值最优法求解的标准形式，并将经Schur降阶后的状态矩阵和时滞矩阵以及统计得到的随机激励强度值带入标准形式，求解得到随机时滞系统的时滞稳定上限；

所述结果输出模块用于输出时滞稳定上限结果。

本发明提出了一种基于伊藤微分的时滞电力系统随机稳定性分析方法，包括：

1.采集时滞电力系统网络结构参数，包括：线路阻抗值、母线电压幅值和相位、发电机有功与无功出力、有功负荷与无功负荷、发电机电抗后电势与发电机励磁电势、发电机阻尼系数、发电机定子时间常数与发电机励磁回路时间常数、发电机励磁回路放大系数、发电机励磁电势与发电机机端电压参考值、发电机稳态电抗、发电机暂态电抗、发电机次暂态电抗、发电机功角与转速；

2.时滞上限的求解

考虑随机因素的时滞电力系统模型可表示为：

其中，x(t)∈Rⁿ是系统的状态向量，时滞h_t满足矩阵A,A_d为已知实矩阵，W(t)是定义在完备概率空间(Ω,F,P)上的Wiener过程，且满足ε(dW(t))＝0,ε(dW²(t))＝dt。

所述系统状态向量包括：发电机功角、发电机角速度、发电机电抗后电势、励磁系统电压调节器输出电压、励磁系统输出电压、励磁系统励磁反馈电压、原动机与调速器的ΔP_m,Δμ,Δx_m，m＝1表示为汽轮机，m＝2表示为水轮机、电力系统稳定器的Δy₁,Δy₂,Δy₃，其中，ΔP_m为机械功率增量，Δμ为汽门开度增量，Δx₁为硬反馈量增量，Δx₂为总反馈量增量；对于PSS系统,Δy₁为复位状态变量增量，Δy₂为第一次超前滞后补偿状态变量增量，Δy₃为第二次超前滞后补偿状态变量增量。

所述完备概率空间(Ω,F,P)中，Ω为样本空间,F为Ω子集的σ代数,P为F上的概率测度。本发明中，W(t)是发电机角速度的完备概率空间。

然后定义向量城y(t)∈R_n，使其满足：

构造如下形式的Lyapunov-Krasovskii泛函：

$\begin{array}{c}V(x_t,t)=x^T\left(t\right)\mathrm{Px}\left(t\right)+{\&Integral;}_{t-h_t}^t{x}^T\left(s\right)\mathrm{Qx}\left(s\right)\mathrm{ds}\\ +{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t{x}^T\left(s\right)\mathrm{Kx}\left(s\right)\mathrm{ds}\\ +{\&Integral;}_{-\stackrel{\&OverBar;}{h}}^0{\&Integral;}_{t+\mathrm{\&theta;}}^t{y}^T\left(s\right)\mathrm{Ry}\left(s\right)\mathrm{dsd\&theta;}\\ +{\&Integral;}_{-\stackrel{\&OverBar;}{h}}^0{\&Integral;}_{t+\mathrm{\&theta;}}^t{f}^T\left(s\right)\mathrm{Zf}\left(s\right)\mathrm{dsd\&theta;}\end{array}---\left(3\right)$

由伊藤微分公式可知，V(x_t,t)沿着系统(3)的随机微分为：

dV(t,x(t))＝ζVdt+2x^T(t)Pgdω(t)           (4)

式(4)中，弱无穷小生成算子为：

$\begin{array}{c}\mathrm{\&zeta;V}(t,x\left(t\right))=2x^T\left(t\right)\mathrm{Pf}\left(t\right)+x^T\left(t\right)(Q+K)x\left(t\right)\\ -(1-{\stackrel{\&CenterDot;}{h}}_t)x^T(t-h_t)\mathrm{Qx}(t-h_t)\\ -x^T(t-\stackrel{\&OverBar;}{h})\mathrm{Kx}(x-\stackrel{\&OverBar;}{h})\\ +\stackrel{\&OverBar;}{h}{y}^T\left(t\right)\mathrm{Ry}\left(t\right)-{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t{y}^T\left(t\right)\mathrm{Ry}\left(t\right)\mathrm{ds}\\ +\stackrel{\&OverBar;}{h}{f}^T\left(t\right)\mathrm{Zf}\left(t\right)-{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t{f}^T\left(t\right)\mathrm{Zf}\left(t\right)\mathrm{ds}\\ +\mathrm{tr}\left\{g^T\left(x\left(t\right)\right)\mathrm{Pg}\left(x\left(t\right)\right)\right\}\end{array}---\left(5\right)$

由Newton-Leibniz公式可知，对于任意的适维矩阵L和M，下列式子成立：

$2{\mathrm{\&xi;}}^T\left(t\right)L[x\left(t\right)-x(t-\stackrel{\&OverBar;}{h})-{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t\mathrm{dx}\left(t\right)]=0---\left(6\right)$

2y^T(t)M[(Ax+A_dx_t-ht-y(t))dt+gdω(t)]＝0       (7)

其中， $\mathrm{\&xi;}\left(t\right)={\left[\begin{array}{ccccc}{x}^T& x^T(t-h_t)& x^T(t-\stackrel{\&OverBar;}{h})& y^T\left(t\right)& f^T\left(t\right)]\end{array}\right]}^T.$

同时，由f(t)＝Ax(t)+A_dx(t-h_t)，可知对于任意的适维矩阵S，以下等式恒成立：

2ζ^T(t)S[Ax(t)-A_dx(t-h_t)-f(t)]＝0      (8)

考虑如下不等式：

±2x^Ty≤x^TP_- ¹x+y^TPy           (9)

其中，x、y是具有适当维数的任意向量，T为任意对阵正定矩阵。

进一步分别利用Jensen不等式和式(9)，可得以下不等式分别成立：

$-{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t{y}_s^{T}R{y}_s\mathrm{ds}\&le;\left[\begin{array}{c}x\\ x_{t-h_t}\end{array}\right]\left[\begin{array}{cc}-R/\stackrel{\&OverBar;}{h}& R/\stackrel{\&OverBar;}{h}\\ R/\stackrel{\&OverBar;}{h}& -R/\stackrel{\&OverBar;}{h}\end{array}\right]\left[\begin{array}{cc}x& x_{t-h_t}\end{array}\right]---\left(10\right)$

$\begin{array}{c}-2{\mathrm{\&xi;}}^T\left(t\right)N{\&Integral;}_{t-h_t}^{t}d{x}_s=-2{\mathrm{\&xi;}}^T\left(t\right)N{\&Integral;}_{t-h_t}^{t}f\left(s\right)\mathrm{ds}\\ -2{\mathrm{\&xi;}}^T\left(t\right)N{\&Integral;}_{t-h_t}^{t}g\left(s\right)\mathrm{d\&omega;}\left(s\right)\\ \&le;-2{\mathrm{\&zeta;}}^T\left(t\right)N{\&Integral;}_{t-h_t}^{t}f\left(s\right)\mathrm{ds}\\ +{\mathrm{\&xi;}}^T\left(t\right)N{T}_1^{-1}{N}^T\mathrm{\&xi;}\left(t\right)\\ +{\left({\&Integral;}_{t-h_t}^{t}g\left(x\right)\mathrm{d\&omega;}\left(s\right)\right)}^T{T}_1{\&Integral;}_{t-h_t}^{t}g\left(s\right)\mathrm{d\&omega;}\left(s\right)\end{array}---\left(11\right)$

综合以上分析，将式(6)、(7)、(8)、(10)和(11)代入式(5)，整理可得：

$\begin{array}{c}\mathrm{\&zeta;V}(t,x\left(t\right))\&le;{\mathrm{\&xi;}}^T\left(t\right)(\mathrm{\&Omega;}+L{T_1}^{-1}{L}^T+\stackrel{\&OverBar;}{h}\&CenterDot;L{Z}^{-1}{L}^T)\mathrm{\&xi;}\left(t\right)\\ -{\&Integral;}_{t-h_t}^t[{\mathrm{\&xi;}}^T\left(t\right)L+f^T\left(t\right)Z]Z^{-1}[L^T\mathrm{\&xi;}\left(t\right)+\mathrm{Zf}\left(t\right)]\mathrm{ds}\end{array}---\left(12\right)$

由于Z>0，因此保证下式成立：

$\mathrm{\&Omega;}+L{T}_2^{-1}{L}^T+\stackrel{\&OverBar;}{h}\&CenterDot;L{Z}^{-1}{L}^{\mathrm{TT}}<0---\left(13\right)$

若对于任意的ξ(t)≠0，均有Ε{ζV(t,x(t))}<0成立，则随机系统(1)在概率空间上渐近均方稳定。基于Schur补定理和式(13)，时滞系统随机稳定判据如下：

对于式(1)所示的随机时滞电力系统，给定标量和μ，若存在P＝P^T>0，Q＝Q^T≥0，R＝R^T≥0，K＝K^T>0，Z＝Z^T≥0，以及适维矩阵L、M和S，使得如下线性矩阵不等式(14)成立，则随机时滞系统(1)均方稳定。

$\left[\begin{array}{cc}\mathrm{\&Omega;}& L\\ *& -T_1-Z/\stackrel{\&OverBar;}{h}\end{array}\right]<0---\left(14\right)$

其中：

$\mathrm{\&Omega;}={\mathrm{\&Omega;}}_1+{\mathrm{\&Omega;}}_2+{\mathrm{\&Omega;}}_2^T$

Ω₂＝[L+SA SA_d-L 0 -S]

${\mathrm{\&Omega;}}_1=\left[\begin{array}{ccccc}\begin{array}{c}Q+K+G_1^{T}P{G}_1\\ -R/\stackrel{\&OverBar;}{h}\end{array}& R/\stackrel{\&OverBar;}{h}& 0& A^{T}M& P\\ *& \begin{array}{c}{G}_2^{T}P{G}_2-(1-\mathrm{\&mu;})Q\\ -R/\stackrel{\&OverBar;}{h}\end{array}& 0& A_d^{T}M& 0\\ *& *& -K& 0& 0\\ *& *& *& \stackrel{\&OverBar;}{h}R-M-M^T& 0\\ *& *& *& *& \stackrel{\&OverBar;}{h}Z\end{array}\right]$

其中，式(14)中变量的定义：P、Q、R、K、Z、N、T₁、T₂为具有适当维数的正定矩阵，是待求变量；L、M、S为具有适当维数的实矩阵，是待求变量。

式(14)表征的矩阵不等式仅能判定系统是否稳定，而无法获取系统时滞稳定上限等信息，考虑到时滞稳定上限的求解是一个具有线性不等式约束的凸优化问题，具有广义特征值(GEVP)的形式。因此，本发明提出利用gevp方法计算系统的时滞稳定上限。由于式(14)不是标准的gevp形式，需要进行必要的处理。根据Schur补性质，式(14)可变形为：

$\left[\begin{array}{cccc}\mathrm{\&Omega;}& L& 0& 0\\ *& -T_1-Z/\stackrel{\&OverBar;}{h}& R& 0\\ *& *& -R/\stackrel{\&OverBar;}{h}& Z\\ *& *& *& -Z/\stackrel{\&OverBar;}{h}\end{array}\right]<0---\left(15\right)$

令 $-Z/\stackrel{\&OverBar;}{h}<-Z_1,-R/\stackrel{\&OverBar;}{h}<-R_1$ 以及 $d=1/\stackrel{\&OverBar;}{h},$ 得：

$\left[\begin{array}{ccc}{T}_1& 0& 0\\ *& T_2& 0\\ *& *& T_3\end{array}\right]<d\&CenterDot;\left[\begin{array}{ccc}{Z}_1& 0& 0\\ *& Z_2& 0\\ *& *& K\end{array}\right]---\left(16\right)$

在式(16)中，将-Z₁,-R₁分别代替得到：

$\left[\begin{array}{cccc}\mathrm{\&Omega;}& L& 0& 0\\ *& -T_1-Z_1& R& 0\\ *& *& -R_1& Z\\ *& *& *& -Z\end{array}\right]<0---\left(17\right)$

由式(15)-(17)，时滞稳定上限问题即可转化为如下优化问题：

$\begin{array}{c}\min \mathrm{\&tau;}\\ s.t.\left(16\right),\left(17\right)\end{array}---\left(18\right)$

其中，τ为待求的实数，是未知变量；通过求解式(18)，以式(16)和式(17)为约束的最小τ，最终利用可求出随机系统时滞稳定上限。

3.输出时滞稳定上限结果数值。

实施例

基于MATLAB仿真软件搭建的图2所示的IEEE16机68节点系统考查考虑电力系统随机特性的时滞随机稳定性分析的有效性和通用性。该系统可分为5大区域。区域1、2和3分别为等值系统，区域4为纽约系统，区域5为新英格兰系统，区域4和5之间有三条联络线，分别为线路1-2、线路8-9和线路1-27，其中线路1-2和线路8-9均为双回线。发电机采用经典模型，励磁采用IEEE-DC2型励磁，负荷模型为50％的恒有功负荷，50％的恒无功负荷。

1过负荷运行状态

固定其他参数，16机系统负荷整体大小处于1.2倍额定负荷，随机激励强度分别取0.03和0.06。将16机系统经Schur降阶后的状态矩阵A和时滞矩阵Ad以及随机激励强度值代入式(16)和式(17)，求得随机系统在和时的时滞稳定上限分别为h_0.03＝111.5ms和h_0.06＝98.7ms。为验证本文所提方法的可行性及有效性，基于H2/H∞控制方法对IEEE16机68节点系统进行阻尼控制器设计，当时，将时滞分别设置为0ms，50ms，111.5ms和150ms，观察发电机1-13之间的相对功角差在不同时滞下的动态响应曲线，如图3所示，当时，将时滞分别设置为0ms，50ms，98.7ms和130ms，观察发电机8-16之间的相对功角差在不同时滞下的动态响应曲线，如图4所示。

由图3和图4可知，当时滞小于时滞稳定上限时，随机电力系统处于稳定状态，且系统能够在10s内迅速阻尼区间振荡；当时滞等于时滞稳定上限时，随机系统处于弱阻尼状态，发电机的相对功角差出现一定摆动；当时滞大于时滞稳定上限时，随机激励下的电力系统由弱阻尼状态发展为不稳定状态，功角差曲线将等幅甚至增幅振荡。

利用prony算法分析发电机1-13和发电机8-16在不同时滞下功角差曲线的阻尼比，结果如表1和表2所示。由表1和表2可知，随机激励为0.03，时滞达到111.5ms时，功角差曲线的阻尼比为12.35％，随机激励为0.06，时滞达到98.7ms时，功角差曲线的阻尼比为12.35％。在这两种随机激励下，虽然系统仍存在一定的阻尼，但其阻尼比已经降至13％以下，阻尼控制器将无法满足控制要求。对比随机系统在和下的时滞稳定上限，发现随机电力系统在随机激励强度由0.03增加到0.06时，时滞稳定上限值减小，说明随机激励的增大导致系统稳定性变差。同时，表1和表2的分析结果表明本方法在系统处于过负荷时，能够有效求解电力系统在随机激励下的时滞稳定上限。

表1时,16机系统G1与G13功角差各时滞时间下的阻尼比

表2时,16机系统G8与G16功角差各时滞时间下的阻尼比

2弱阻尼运行状态

固定其他参数，16台发电机阻尼系数D取额定阻尼系数的50％，随机激励强度分别取0.03和0.06。利用降阶后的状态矩阵A和时滞矩阵Ad以及随机激励强度值，求得随机系统在随机激励强度为0.02和0.06时，时滞稳定上限分别为h_0.03＝97.1ms和h_0.06＝83.6ms。当时，观察发电机1-16之间的相对功角差在不同时滞下的动态响应曲线，如图5所示。当时，在不同时滞下，发电机6-15之间的相对功角差动态响应曲线，如图6所示。

由图5和图6可知，当时滞小于时滞稳定上限时，发电机相对功角差在10s后基本保持为恒定值，随机电力系统处于稳定状态，而且系统能够在10s内迅速阻尼区间振荡；当时滞等于时滞稳定上限时，随机系统处于弱阻尼状态，发电机的相对功角差出现一定的摆动；当时滞大于时滞稳定上限时，发电机功角差曲线出现等幅甚至增幅振荡，随机激励下的电力系统从弱阻尼状态发展为不稳定状态。

表3时,16机系统G1与G16功角差各时滞时间下的阻尼比

表4时,16机系统G6与G15功角差各时滞时间下的阻尼比

表3和表4分别为发电机1-16和发电机6-15在不同时滞下功角差曲线的阻尼比。由表3和表4可知，时滞小于时滞稳定上限时，随机电力系统的阻尼比较大，处于正常运行状态，但是当时滞达到时滞稳定上限时，阻尼比已经降至12％以下，分别为10.89％和11.46％，此时随机激励下的系统仍存在一定的阻尼，但阻尼控制器无法满足控制要求。通过分析图5、图6、表3以及表4可知，本文提出的随机电力系统时滞稳定性分析方法在弱阻尼运行工况下可行，且能够有效求解弱阻尼电力系统在随机激励下的时滞稳定上限。

3高励磁放大倍数运行状态

固定其他参数，16台发电机励磁放大倍数Ka取额定励磁放大倍数的110％，随机激励强度分别取0.03和0.06。利用降阶后的状态矩阵A和时滞矩阵Ad以及随机激励强度值，求得随机系统在不同随机激励强度下的时滞稳定上限分别为h_0.03＝126.3ms和h_0.06＝99.4ms。当时，不同时滞下发电机9-15之间的相对功角差动态响应曲线如图7所示，当时，观察发电机4-13之间的相对功角差动态响应曲线，如图8所示。

由图7和图8可以看出，当时滞小于时滞稳定上限时，系统可在10s内迅速阻尼区间振荡，10s后发电机相对功角差为恒定值，随机电力系统处于稳定状态；当时滞等于时滞稳定上限时，发电机的相对功角差在恒定值附近波动，随机系统处于弱阻尼状态；当时滞大于时滞稳定上限时，随机激励下的电力系统从弱阻尼状态发展为不稳定状态，相对功角差曲线出现等幅甚至增幅振荡。

表5时,16机系统G9与G15功角差各时滞时间下的阻尼比

表5和表6分别给出了发电机9-15和发电机4-13在不同时滞下功角差曲线的阻尼比。由表5和表6可知，随机激励为0.03，时滞达到126.3ms时，功角差曲线的阻尼比为8.88％，随机激励为0.06，时滞达到99.4ms时，功角差曲线的阻尼比为10.84％，阻尼比已经降至11％以下，随机激励下的电力系统处于弱阻尼状态。通过分析图7、图8、表5以及表6可知，本文提出的基于伊藤微分分析随机电力系统时滞稳定性的方法在高励磁放大倍数运行工况下可行，且能够有效求解电力系统在随机激励下的时滞稳定上限。

表6时,16机系统G4与G13功角差各时滞时间下的阻尼比

4保守性分析

考虑到公平性，基于如下系统，验证本文所提方法的低保守性，系统参数如下：

$\begin{array}{c}{A}_1=\left[\begin{array}{cc}-2& 0\\ -1& -1\end{array}\right],A_{d1}=\left[\begin{array}{cc}-1& 0\\ -0.5& -1\end{array}\right]\\ G_1=G_2=\sqrt{0.1}*\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\end{array}---\left(19\right)$

对于不同的μ，利用LMI算法求出不同方法对应的时滞稳定上限如表7所示。对比表中结果，可知本方法具有较低的保守性。

表7不同方法所能得到的时滞稳定上限

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (10)

1.一种基于伊藤微分的时滞电力系统随机稳定性分析方法，其特征在于，包括：

步骤1、采集时滞电力系统网络结构参数；

步骤2、将采集到的时滞电力系统网络结构参数构建得到系统的状态矩阵和时滞矩阵，并与系统的状态向量一起根据时滞系统方程来构造考虑随机因素的Lyapunov-Krasovskii泛函，并借助伊藤微分公式求解目标泛函的弱无穷小算子；

步骤3、利用Newton-Leibniz公式建立与自由权矩阵L和M有关的等式，并将该等式与必要的松散项加入到Lyapunov-Krasovskii泛函的弱无穷小算子之中，得到时滞系统稳定判据；

步骤4、将时滞系统稳定判据等价变换成符合广义特征值最优法求解的标准形式，并将经Schur降阶后的状态矩阵和时滞矩阵以及统计得到的随机激励强度值带入标准形式，求解得到随机时滞系统的时滞稳定上限。

2.根据权利要求1所述方法，其特征在于，所述步骤2中的Lyapunov-Krasovskii泛函为：

$\begin{array}{c}V(x_t,t)=x^T\left(t\right)\mathrm{Px}\left(t\right)+{\&Integral;}_{t-h_t}^t{x}^T\left(s\right)\mathrm{Qx}\left(s\right)\mathrm{ds}\\ +{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t{x}^T\left(s\right)\mathrm{Kx}\left(s\right)\mathrm{ds}\\ +{\&Integral;}_{-\stackrel{\&OverBar;}{h}}^0{\&Integral;}_{t+\mathrm{\&theta;}}^t{y}^T\left(s\right)\mathrm{Ry}\left(s\right)\mathrm{dsd\&theta;}\\ +{\&Integral;}_{-\stackrel{\&OverBar;}{h}}^0{\&Integral;}_{t+\mathrm{\&theta;}}^t{f}^T\left(s\right)\mathrm{Zf}\left(s\right)\mathrm{dsd\&theta;}\end{array}---\left(3\right)$

其中，x(t)∈Rⁿ是系统的状态向量，时滞h_t满足 为时滞稳定上限，为时滞h_t的一阶导数，μ为时滞的最大变化率，向量城y(t)∈R_n满足： P、Q、R、K、Z为具有适当维数的正定矩阵，是待求变量，为时滞稳定上限，f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵；s,θ为替换变量t的多重积分变量。

3.根据权利要求1所述方法，其特征在于，所述步骤2中Lyapunov-Krasovskii泛函的弱无穷小算子为：

$\begin{array}{c}\mathrm{\&zeta;V}(t,x\left(t\right))=2x^T\left(t\right)\mathrm{Pf}\left(t\right)+x^T\left(t\right)(Q+K)x\left(t\right)\\ -(1-{\stackrel{.}{h}}_t)x^T(t-h_t)\mathrm{Qx}(t-h_t)\\ -x^T(t-\stackrel{\&OverBar;}{h})\mathrm{Kx}(t-\stackrel{\&OverBar;}{h})\\ +\stackrel{\&OverBar;}{h}{y}^T\left(t\right)\mathrm{Ry}\left(t\right)-{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t{y}^T\left(t\right)\mathrm{Ry}\left(t\right)\mathrm{ds}\\ +\stackrel{\&OverBar;}{h}{f}^T\left(t\right)\mathrm{Zf}\left(t\right)-{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t{f}^T\left(t\right)\mathrm{Zf}\left(t\right)\mathrm{ds}\\ +\mathrm{tr}\left\{g^T\left(x\left(t\right)\right)\mathrm{Pg}\left(x\left(t\right)\right)\right\}\end{array}---\left(5\right)$

其中，x(t)∈Rⁿ是系统的状态向量，时滞h_t满足 为时滞稳定上限，为时滞h_t的一阶导数，μ为时滞的最大变化率，向量城y(t)∈R_n满足： P、Q、R、K、Z为具有适当维数的正定矩阵，是待求变量，f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵；s,θ为替换变量t的多重积分变量。

4.根据权利要求1所述方法，其特征在于，所述步骤3中改进自由权矩阵L和M有关的等式为：

$2{\mathrm{\&xi;}}^T\left(t\right)L[x\left(t\right)-x(t-\stackrel{\&OverBar;}{h})-{\&Integral;}_{t-\stackrel{\&OverBar;}{h}}^t\mathrm{dx}\left(t\right)]=0$

${2y}^T\left(t\right)M[(\mathrm{Ax}+A_d{x}_{t-h_t}-y\left(t\right))\mathrm{dt}+g\left(x\left(t\right)\right)\mathrm{d\&omega;}\left(t\right)]=0$

其中， $\mathrm{\&xi;}\left(t\right)={\left[\begin{array}{ccccc}{x}^T& x^T(t-h_t)& x^T(t-\stackrel{\&OverBar;}{h})& y^T\left(t\right)& f^T\left(t\right)\end{array}\right]}^T,$ (t)∈Rⁿ是状态向量，向量城y(t)∈R_n满足： 为时滞稳定上限，μ为时滞的最大变化率，时滞h_t满足时滞d(t)与一阶导数满足条件：f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵。

5.根据权利要求1所述方法，其特征在于，所述步骤3中必要的松散项为：

2ζ^T(t)S[Ax(t)-A_dx(t-h_t)-f(t)]＝0

其中，f(t)＝Ax(t)+A_dx(t-h_t)，A为状态矩阵，A_d为时滞矩阵，h_t为时滞。

6.根据权利要求1所述方法，其特征在于，所述步骤3中时滞系统随机稳定判据如下：

考虑随机因素的时滞电力系统模型表示为：

其中，x(t)∈Rⁿ是系统的状态向量，时滞h_t满足矩阵A,A_d为已知实矩阵，W(t)是定义在完备概率空间(Ω,F,P)上的Wiener过程，且满足ε(dW(t))＝0,ε(dW²(t))＝dt；

定义向量城y(t)∈R_n，使其满足：

对于式(1)所示的随机时滞电力系统，给定标量和μ，若存在P＝P^T>0，Q＝Q^T≥0，R＝R^T≥0，K＝K^T>0，Z＝Z^T≥0，以及适维矩阵L、M和S，使得如下线性矩阵不等式(14)成立，则随机时滞系统(1)均方稳定；

$\left[\begin{array}{cc}\mathrm{\&Omega;}& L\\ *& -T_1-Z/\stackrel{\&OverBar;}{h}\end{array}\right]<0---\left(14\right)$

其中：

$\mathrm{\&Omega;}={\mathrm{\&Omega;}}_1+{\mathrm{\&Omega;}}_2+{\mathrm{\&Omega;}}_2^T$

Ω₂＝[L+SA SA_d -L 0 -S]

${\mathrm{\&Omega;}}_1=\left[\begin{array}{ccccc}Q+K+G_1^{T}P{G}_1& R/\stackrel{\&OverBar;}{h}& 0& A^{T}M& P\\ -R/\stackrel{\&OverBar;}{h}& & & & \\ *& G_2^{T}P{G}_2-(1-\mathrm{\&mu;})Q& 0& A_d^{T}M& 0\\ & -R/\stackrel{\&OverBar;}{h}& & & \\ *& *& -K& 0& 0\\ *& *& *& \stackrel{\&OverBar;}{h}R-M-M^T& 0\\ *& *& *& *& \stackrel{\&OverBar;}{h}Z\end{array}\right]$

其中，P、Q、R、K、Z、T₁为具有适当维数的正定矩阵，是待求变量；L、M、S为具有适当维数的实矩阵，是待求变量；Ω、Ω₁、Ω₂为具有适当维数的矩阵，G₁、G₂为已知变量。

7.根据权利要求1所述方法，其特征在于，所述时滞电力系统网络结构参数包括：线路阻抗值、母线电压幅值和相位、发电机有功与无功出力、有功负荷与无功负荷、发电机电抗后电势与发电机励磁电势、发电机阻尼系数、发电机定子时间常数与发电机励磁回路时间常数、发电机励磁回路放大系数、发电机励磁电势与发电机机端电压参考值、发电机稳态电抗、发电机暂态电抗、发电机次暂态电抗、发电机功角与转速。

8.根据权利要求1所述方法，其特征在于，所述系统状态向量包括：发电机功角、发电机角速度、发电机电抗后电势、励磁系统电压调节器输出电压、励磁系统输出电压、励磁系统励磁反馈电压、原动机与调速器的ΔP_m,Δμ,Δx_m，m＝1表示为汽轮机，m＝2表示为水轮机、电力系统稳定器的Δy₁,Δy₂,Δy₃，其中，ΔP_m为机械功率增量，Δμ为汽门开度增量，Δx₁为硬反馈量增量，Δx₂为总反馈量增量；对于PSS系统,Δy₁为复位状态变量增量，Δy₂为第一次超前滞后补偿状态变量增量，Δy₃为第二次超前滞后补偿状态变量增量。

9.根据权利要求1所述方法，其特征在于，所述随机激励强度值通过统计得到，种类包括：负荷的随机波动、原动机扭矩的随机振动、控制回路的测量噪声、互联大电网中功率角的随机小幅振荡、风力发电中风速的随机变化。

10.一种基于伊藤微分的时滞电力系统随机稳定性分析系统，其特征在于，包括：数据采集模块、时滞上限求解模块和结果输出模块；其中，数据采集模块通过时滞上限求解模块与结果输出模块相连；

所述数据采集模块用于采集时滞电力系统网络结构参数，并将采集数据发送至时滞上限求解模块；

所述时滞上限求解模块的功能是：

首先将采集到的时滞电力系统网络结构参数构建得到系统的状态矩阵和时滞矩阵，并与系统的状态向量一起根据时滞系统方程来构造考虑随机因素的Lyapunov-Krasovskii泛函，并借助伊藤微分公式求解目标泛函的弱无穷小算子；

然后利用Newton-Leibniz公式建立与自由权矩阵L和M有关的等式，并将该等式与必要的松散项加入到Lyapunov-Krasovskii泛函的弱无穷小算子之中，得到时滞系统稳定判据；

最后将时滞系统稳定判据等价变换成符合广义特征值最优法求解的标准形式，并将经Schur降阶后的状态矩阵和时滞矩阵以及统计得到的随机激励强度值带入标准形式，求解得到随机时滞系统的时滞稳定上限；

所述结果输出模块用于输出时滞稳定上限结果。