时变电力系统稳定性分析系统及方法

摘要

一种时变电力系统稳定性分析系统及方法，所述系统包括数据采集模块，风速拟合预测模块、时变系统建立与降阶模块、稳定判据求解模块和结果输出模块，数据采集模块用于采集数据；风速拟合预测模块对短时间内条件风速趋势进行拟合预测，并将风速区间化，根据拟合模型计算各区间的条件特征风速概率密度矩阵；时变系统建立与降阶模块建立考虑风速随机特征的连续马尔科夫时变电力系统模型，并对系统进行降阶；稳定判据求解模块利用线性矩阵不等式中的可行性问题求解方法求解系统是否稳定。通过本发明的时变电力系统稳定性分析系统及方法，能够有效解决风机接入下风速随机特征导致的时变电力系统稳定性不易判别的难题。

说明书

技术领域

本发明涉及电力系统分析和控制技术领域，特别是涉及到时变电力系统的 稳定性分析技术。

背景技术

风电的发展是世界各国未来能源战略的重要组成部分，但是风机出力的波 动往往导致系统运行工况变化，这为电力系统的安全稳定与经济运行带来了巨 大风险。因此，迫切需要对含风速随机特征的时变电力系统稳定性进行深入分 析。

现有技术中，有关风电接入下电力系统的稳定性分析技术主要分为两类， 一类是基于风机自身的稳定性分析，另一类是基于电网侧单一运行工况下的稳 定性分析。基于风机自身的稳定性分析以不同风机接入下风功率变化对风机自 身稳定性的影响为研究对象，分析不同风机的稳定性差异，然而，这类分析并 未对风机接入后网侧电力系统的稳定性进行分析。基于电网侧单一运行工况下 的稳定性分析以风电接入下风功率变化时网侧电力系统振荡模式为分析对象， 分析网侧电力系统稳定性，然而，该方法仅针对电力系统的单一运行工况进行 分析，并未对风机出力波动较大导致系统运行工况发生大范围变化的情况进行 分析。因此，现有技术中，对于风电接入下电力系统的分析不全面，准确性低。

发明内容

鉴于此，本发明的目的在于克服现有技术的困难，避免仅对风力发电机自 身进行稳定性分析，因而忽视了风力发电机接入后网侧电力系统的稳定性；也 避免了以风电接入下风功率变化时网侧电力系统振荡模式为分析目标，因此忽 视了风力发电机出力波动较大导致系统运行工况发生大范围变化的问题。为此 本发明提出一种考虑风速随机特征的时变电力系统稳定性分析系统及方法。

本发明从电力系统实际状况出发，能够有效解决风机接入下风速随机特征 导致的电力系统时变稳定性不易判别的难题。首先利用双参数韦伯分布模型对 短时间内条件风速趋势进行拟合预测，再将风速区间化，计算各区间的条件特 征风速概率密度矩阵，根据条件特征风速确定系统的运行工况。然后，根据系 统网架信息建立考虑风速随机特征的连续马尔科夫时变电力系统模型，并降阶 系统。其次，本发明构建含连续马尔科夫时变电力系统模型的李雅普诺夫泛函， 并在该泛函的弱无穷小算子中运用邓肯引理，推导满足干扰衰减度的鲁棒随机 稳定性的线性矩阵不等式，最后将稳定性问题转化为可行性问题求解系统是否 稳定。

为了实现此目的，本发明采取的技术方案为如下。

一种时变电力系统稳定性分析系统，所述系统包括依次连接的数据采集模 块，风速拟合预测模块、时变系统建立与降阶模块、稳定判据求解模块和结果 输出模块，数据采集模块还连接至时变系统建立与降阶模块；

所述数据采集模块用于采集风电场风速与风机出力信息、网络结构参数、 系统内发电机频率、功角，并将采集数据发送至风速拟合预测模块和时变系统 建立与降阶模块；

所述风速拟合预测模块根据数据采集模块发送的信息与风电场历史实际数 据，对短时间内条件风速趋势进行拟合预测，并将风速区间化，根据拟合模型 计算各区间的条件特征风速概率密度矩阵；

所述时变系统建立与降阶模块将每个区间作为时变系统的一个子工况，根 据风机风速与出力对应关系以及数据采集模块的采集数据建立考虑风速随机特 征的连续马尔科夫时变电力系统模型，并对系统进行降阶；

所述的稳定判据求解模块用于形成时变电力系统鲁棒随机稳定判据，利用 线性矩阵不等式中的可行性问题求解方法求解系统是否稳定；

所述的结果输出模块用于输出系统稳定性判别结果。

一种时变电力系统稳定性分析方法，所述方法包括：

A、采集风电场风速与风机出力信息、网络结构参数、系统内发电机频率、 功角；

B、根据采集的数据与风电场历史实际数据，对短时间内条件风速趋势进行 拟合预测，并将风速区间化，根据拟合模型计算各区间的条件特征风速概率密 度矩阵；

C、将每个区间作为时变系统的一个子工况，根据风机风速与出力对应关系 以及采集的数据建立考虑风速随机特征的连续马尔科夫时变电力系统模型，并 对系统进行降阶；

D、形成时变电力系统鲁棒随机稳定判据，利用线性矩阵不等式中的可行性 问题求解方法求解系统是否稳定；

E、输出时变电力系统稳定分析结果。

步骤B中，所述的条件特征风速概率密度为某一条件特征风速下，预定时 间间隔后的风速概率密度。

另外步骤B中，对短时间内条件风速趋势进行拟合预测，并将风速区间化， 根据拟合模型计算各区间的条件特征风速概率密度矩阵包括：

B1、利用双参数韦伯分布模型拟合时变系统各条件风速下的变化趋势；

B2、按照一定的区间长度将风速在切入风速与切出风速之间进行区间化， 利用步骤B1中各条件风速下的拟合模型确定各区间的条件特征风速概率密度矩 阵。

其中，所述双参数韦伯分布模型为：

$g\left(v_i\right|v_j)=\left(\frac{k}{c}\right){\left(\frac{v_i}{c}\right)}^{k-1}\exp [-{\left(\frac{v_i}{c}\right)}^k],$

其中，v_i为当前时刻风速；

v_j为条件风速；

k为形状系数；

c为尺度系数；

g(v_i|v_j)为条件风速为v_j时当前风速为v_i的概率密度。

步骤C中，根据风机风速与出力对应关系以及采集的数据建立考虑风速随 机特征的连续马尔科夫时变电力系统模型，并对系统进行降阶包括：

C1、将每个区间作为时变系统的一个子工况，根据风机风速与出力的对应 关系确定各区间特征风速对应的风机输出功率，并利用该风机输出功率以及采 集的数据建立考虑风速随机特征的连续马尔科夫时变电力系统模型；

C2、在保证所关心频带输入输出特性不变的前提下对连续马尔科夫时变电 力系统模型降阶。

其中，所述风机风速与出力对应关系为：

$P\left(v_i\right)=\left(\begin{array}{cc}0& 0\le v_i<v_{\mathrm{ci}}{v}_i\ge v_{\mathrm{co}}\\ \frac{v_i-v_{\mathrm{ci}}}{v_R-v_{\mathrm{ci}}}{P}_R& v_{\mathrm{ci}}<v_i<v_R\\ P_R& v_R\le v_i<v_{\mathrm{co}}\end{array}\right),$

其中，v_ci为风机切入风速；

v_co为风机切出风速；

v_R为风机额定风速；

P_R为风机额定输出功率。

特别地，步骤C2中所述对系统进行降阶包括：

C21、将系统状态方程划分为：

$\left(\begin{array}{c}{\stackrel{·}{X}}_1\\ {\stackrel{·}{X}}_2\end{array}\right)=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right),$

其中，X₁＝[Δω^T,Δδ^T]为保留变量，Δω表示发电机功角转速的变化量，Δδ 表示发电机功角变化量，X₂为待消去的其他变量；

C22、消去X₂后为X₁，其中I为单位阵，p为 微分算子。

步骤D中，形成时变电力系统鲁棒随机稳定判据，利用线性矩阵不等式中 的可行性问题求解方法求解系统是否稳定包括：

D1、构建含连续马尔科夫时变电力系统模型的李雅普诺夫泛函，并在该泛 函的弱无穷小算子中适用邓肯引理，获得满足干扰衰减度γ的鲁棒随机稳定性的 线性矩阵不等式，得到时变电力系统鲁棒随机稳定判据；

D2、利用线性矩阵不等式中的可行性问题求解方法确定系统是否稳定。

其中，所述鲁棒随机稳定性判据具体可以表示为：给定正常数γ>0，如果 存在一组正定对称矩阵P_i>0，i∈S，使得如下一组矩阵不等式成立：

$\left(\begin{array}{cc}{J}_i& C_i^T{L}_i+P_i{G}_i\\ *& L_i^T{L}_i-{\gamma }^{2}I\end{array}\right)<0,$

$J_i=A_i^T{P}_i+P_i{A}_i+\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+C_i^T{C}_i,$

则当u(t)≡0时，时变系统鲁棒随机稳定，且满足扰动衰减度γ，即：

${\left|\right|z\left|\right|}_2\le \gamma {[{\left|\right|\omega \left|\right|}_2^2+x_0^{T}P\left(s_0\right)x_0]}^{\frac{1}{2}}.$

通过采用本发明的时变电力系统稳定性分析系统及方法，能够有效解决风 机接入下风速随机特征导致的时变电力系统稳定性不易判别的难题。另外，本 发明的时变电力系统稳定性分析系统及方法无需获取系统运行轨迹即可对时变 电力系统进行稳定性判别，与现有技术中的时域仿真法相比，降低了计算量， 提高了判别效率。

附图说明

图1为本发明实施方式中时变电力系统稳定性分析系统的结构示意图。

图2为风机功率曲线示意图。

图3为条件风速模式信息表。

图4为风速模式下条件风速概率密度的韦伯分布拟合系数表。

图5为本发明一个应用示例的IEEE4机11节点系统结构示意图。

图6为本发明一个应用示例的4机系统单工况变化情况下，发电机2-3相对 功角动态响应示意图。

图7为本发明一个应用示例的4机系统多工况变化情况下，发电机2-3相对 功角动态响应示意图

图8为本发明另一个应用示例的IEEE16机68节点系统结构示意图。

图9为本发明另一个应用示例的16机系统单工况变化情况下，发电机1-8 相对功角动态响应示意图。

图10为本发明另一个应用示例的16机系统多工况变化情况下，发电机1-8 相对功角动态响应示意图。

图11为本发明另一个应用示例的16机系统各运行模式最小阻尼比。

具体实施方式

下面结合附图，对本发明作详细说明。

以下公开详细的示范实施例。然而，此处公开的具体结构和功能细节仅仅 是出于描述示范实施例的目的。

然而，应该理解，本发明不局限于公开的具体示范实施例，而是覆盖落入 本公开范围内的所有修改、等同物和替换物。在对全部附图的描述中，相同的 附图标记表示相同的元件。

同时应该理解，如在此所用的术语“和/或”包括一个或多个相关的列出 项的任意和所有组合。另外应该理解，当部件或单元被称为“连接”或“耦接” 到另一部件或单元时，它可以直接连接或耦接到其他部件或单元，或者也可以 存在中间部件或单元。此外，用来描述部件或单元之间关系的其他词语应该按 照相同的方式理解(例如，“之间”对“直接之间”、“相邻”对“直接相邻” 等)。

为了介绍本发明的技术方案，首先说明本发明的技术原理。

风机出力受风速随机特征影响呈现多种运行方式，若恰当的对风速的随机 特征进行分析，模拟含风机出力影响的电力系统运行状态，则能对考虑风速随 机特征的时变电力系统稳定性进行有效分析。

定义条件风速概率密度g(v_i|v_j)为某一条件特征风速下，预定时间间隔之后 (例如采样间隔通常选为10min，所述预定时间间隔为这样一个采样间隔或几个 采样间隔之后)的风速概率密度。采用双参数韦伯分布拟合风速条件概率密度 曲线。

双参数韦伯分布模型具体表达式为：

$g\left(v_i\right|v_j)=\left(\frac{k}{c}\right){\left(\frac{v_i}{c}\right)}^{k-1}\exp [-{\left(\frac{v_i}{c}\right)}^k]---\left(1\right)$

其中，v_i为当前时刻风速；

v_j为条件风速；

k为形状系数；

c为尺度系数；

g(v_i|v_j)为条件风速为v_j时当前风速为v_i的概率密度。

风速的概率分布函数为：

$G\left(v_i\right)=1-\exp [-{\left(\frac{v_i}{c}\right)}^k],---\left(2\right)$

忽略电气损耗以及风电场尾流等因素，风机出力与风速之间的近似函数关 系可以表示为：

$P\left(v_i\right)=\left(\begin{array}{cc}0& 0\le v_i<v_{\mathrm{ci}}{v}_i\ge v_{\mathrm{co}}\\ \frac{v_i-v_{\mathrm{ci}}}{v_R-v_{\mathrm{ci}}}{P}_R& v_{\mathrm{ci}}<v_i<v_R\\ P_R& v_R\le v_i<v_{\mathrm{co}}\end{array}\right)---\left(3\right)$

其中，v_ci为风机切入风速，

v_co为风机切出风速，

v_R为风机额定风速，

P_R为风机额定输出功率，

P(v_i)为风速为v_i的风机额定输出功率。

由式(1)-(3)得到风机输出功率的分段概率函数为：

$F\left(P\right)=\left(\begin{array}{cc}G\left(v_{\mathrm{ci}}\right)+1-G\left(v_{\mathrm{co}}\right)& P=0\\ G\left(v_i\right)-G\left(v_{\mathrm{ci}}\right)& 0<P<P_R\\ G\left(v_{\mathrm{co}}\right)-G\left(v_R\right)& P=P_R\end{array}\right)---\left(4\right)$

另外，时变电力系统的模型可用连续马尔科夫系统表示为：

$\left(\begin{array}{c}{x}_t=(A\left(s_t\right)+\mathrm{\Delta A}(s_t,t))x_t+B\left(s_t\right)u_t+G\left(s_t\right){\omega }_t\\ z_t=C\left(s_t\right)x_t+D\left(s_t\right)u_t+L\left(s_t\right){\omega }_t\end{array}\right)---\left(5\right)$

其中，x_t∈Rⁿ是状态向量，

u_t∈R^p是控制输入向量，

t是时间，

ω_t为白噪声，满足：$E\left\{d{\omega }_t\right\}=0,E\left\{{\mathrm{d\omega }}_t^2\right\}=\mathrm{Wdt}(W>0),$

{s(t),t≥0}是在有限空间S＝{1,2,…,l}中取值的连续马尔科夫过程，对应各 个风速模式下的系统运行工况，描述系统工况随风机出力变化的演化过程。其 状态转移概率密度可以表示为：

$p_{\mathrm{ij}}=P({\eta }_{t+\Delta }=j|{\eta }_t=i)=\left(\begin{array}{cc}{\pi }_{\mathrm{ij}}\Delta +o\left(\Delta \right),& i\ne j\\ 1+{\pi }_{\mathrm{ii}}\Delta +o\left(\Delta \right),& i=j\end{array}\right)---\left(6\right)$

其中，$\underset{\Delta \to \infty }{\lim }\frac{o\left(\Delta \right)}{\Delta }=0,(\Delta >0),$

η_t为时刻t时的运行工况，Δ为时刻t的变化量，

o(Δ)是时刻t的变化量Δ的高阶无穷小量，

矩阵π表示马尔科夫转移概率密度矩阵，π_ij是系统运行工况在t时刻处于i， 而在t+Δ时刻处于j的转换概率密度，并且有：

${\pi }_{\mathrm{ii}}=-\underset{i\ne j}{\Sigma }{\pi }_{\mathrm{ij}}({\pi }_{\mathrm{ij}}\ge 0,i\ne j)---\left(7\right)$

对于每一个s_t＝i∈S，记A(s_t)、B(s_t)、G(s_t)、C(s_t)、D(s_t)、L(s_t)分别表 示为A_i、B_i、G_i、C_i、D_i、L_i，并且不确定参数满足匹配性条件：

ΔA(s_t,t)＝ΔA(i,t)＝H_iF(i,t)M_i         (8)

其中，H_i和M_i为已知矩阵，I为单位矩阵，实矩阵F(i,t)反映了系统不确 定参数的结构信息，满足：

F^T(i,t)F(i,t)≤I         (9)

接下来根据系统的状态方程说明稳定性分析过程。

对于式(5)表示的时变电力系统的模型，当u(t)≡0以及所有的初始条件 x₀∈Rⁿ和s₀∈S成立时，若满足：

${\int }_0^{\infty }E\left\{{\left|\right|x(t,x_0,s_0)\left|\right|}^2\right\}\mathrm{dt}<\infty ---\left(10\right)$

那么，形如式(5)的时变电力系统随机稳定。

另外，鲁棒随机稳定性判据具体可以表示为：

给定正常数γ>0，如果存在一组正定对称矩阵P_i>0，i∈S，使得如下一组 矩阵不等式成立：

$\left(\begin{array}{cc}{J}_i& C_i^T{L}_i+P_i{G}_i\\ *& L_i^T{L}_i-{\gamma }^{2}I\end{array}\right)<0---\left(11\right)$

$J_i=A_i^T{P}_i+P_i{A}_i+\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+C_i^T{C}_i---\left(12\right)$

则当u(t)≡0时，时变系统鲁棒随机稳定，且满足扰动衰减度γ，即：

${\left|\right|z\left|\right|}_2\le \gamma {[{\left|\right|\omega \left|\right|}_2^2+x_0^{T}P\left(s_0\right)x_0]}^{\frac{1}{2}}---\left(13\right)$

利用Schur定理，由式(11)可知：

$A_i^T{P}_i+P_i{A}_i+\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+C_i^T{C}_i<0---\left(14\right)$

由于对i∈S，故：

$A_i^T{P}_i+P_i{A}_i+\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j<0---\left(15\right)$

因此，只需确定式(15)成立的条件下，系统(5)鲁棒随机稳定，即可说明鲁棒 随机稳定判据成立。

因此，构造李雅普诺夫泛函为：

$V(x_t,s_t=i)=V(x_t,i)=x_t^T{P}_i{x}_i(P_i>0,i\in S)---\left(16\right)$

当ω(t)＝0时，由泛函V(x_t,i)的弱无穷小算子可知：

$\mathrm{lV}(x_t,i)=\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}V(x_t,i)+x_t^T{A}_i^T\frac{\partial }{\partial x}V(x_t,i)=x_t^T[\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+A_i^T{P}_j+P_j{A}_i]x_t<0---\left(17\right)$

即存在一组正定矩阵Q_i>0(i∈S)，使得：

$\mathrm{lV}(x_t,i)=-x_t^T{Q}_i{x}_t,$

从而有：

$\mathrm{lV}(x_t,i)\le -\underset{i\in S}{\min }{\lambda }_{\min }\left[Q_i\right]x_t^T{x}_t<0---\left(18\right)$

由邓肯引理，有：

$\left(\begin{array}{c}E\left[V(x_t,i)\right]-E\left[V(x_0,s_0)\right]=E\left[{\int }_0^t\mathrm{lV}(x_r,s_r)\mathrm{dr}\right]\\ \le -\underset{i\in S}{\min }{\lambda }_{\min }\left[Q_i\right]E\left[{\int }_0^t{x}_r^T{x}_r\mathrm{dr}\right|(x_0,s_0)]\end{array}\right)---\left(19\right)$

由上式可知：

$\left(\begin{array}{c}\underset{i\in S}{\min }{\lambda }_{\min }\left[Q_i\right]E\left[{\int }_0^t{x}_r^T{x}_r\mathrm{dr}\right|(x_0,s_0)]\\ \le E\left[V(x_0,s_0)\right]-E\left[V(x_t,i)\right]\le E\left[V(x_0,s_0)\right]\end{array}\right)---\left(20\right)$

因此，有：

$E\left[{\int }_0^t{x}_r^T{x}_r\mathrm{dr}\right|(x_0,s_0)]\le \frac{E\left[V(x_0,s_0)\right]}{\underset{i\in S}{\min }{\lambda }_{\min }\left[Q_i\right]}---\left(21\right)$

当t→∞时，有：

$E\left[{\int }_0^{\infty }{x}_r^T{x}_r\mathrm{dr}\right|(x_0,s_0)]\le \frac{E\left[V(x_0,s_0)\right]}{\underset{i\in S}{\min }{\lambda }_{\min }\left[Q_i\right]}<\infty ---\left(22\right)$

因此，由式(15)-(22)知式(5)所表示的时变电力系统随机稳定。

当ω(t)≠0时，泛函V(x_t,i)的弱无穷小算子可表示为：

$\mathrm{lV}(x_t,i)=x_t^T[\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+A_i^T{P}_j+P_j{A}_i]x_t+2x_t^T{P}_i{G}_i{\omega }_t---\left(23\right)$

定义指标函数为：

$J_T=E\left[{\int }_0^T(z_t^T{z}_t-{\gamma }^2{\omega }_t^T{\omega }_t)\mathrm{dt}\right]---\left(24\right)$

由于：

$\left(\begin{array}{c}{z}_t^T{z}_t-{\gamma }^2{\omega }_t^T{\omega }_t+\mathrm{lV}(x_t,i)=x_t^T{C}_i^T{C}_i{x}_t-{\gamma }^2{\omega }_t^T{\omega }_t+x_t^T{C}_i^T{L}_i{\omega }_t+{\omega }_t^T{L}_i^T{C}_i{x}_t\\ +{\omega }_t^T{L}_i^T{L}_i{\omega }_t+x_t^T[\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+A_i^T{P}_j+P_j{A}_i]x_t+{2\omega }_t^T{G}_i^T{P}_i{x}_t+x_t^T{P}_i{G}_i{\omega }_t\\ =\left(\begin{array}{cc}{x}_t^T& {\omega }_t^T\end{array}\right)\left(\begin{array}{cc}{J}_i& C_i^T{L}_i+P_i{G}_i\\ *& L_i^T{L}_i-{\gamma }^{2}I\end{array}\right)\left(\begin{array}{c}{x}_t\\ {\omega }_t\end{array}\right)\end{array}\right)$

因此，

$J_T=E\left[{\int }_0^T(z_t^T{z}_t-{\gamma }^2{\omega }_t^T{\omega }_t+\mathrm{lV}(x_t,i))\mathrm{dt}\right]-E\left[{\int }_0^T\left(\mathrm{lV}(x_t,i)\right)\mathrm{dt}\right]---\left(25\right)$

由邓肯引理，有：

$E\left[{\int }_0^T\left(\mathrm{lV}(x_t,i)\right)\mathrm{dt}\right]=E\left[V(x_T,s_T)\right]-E\left[V(x_0,s_0)\right],$

从而有

$J_T=E\left\{{\int }_0^T\left(\begin{array}{cc}{x}_t^T& {\omega }_t^T\end{array}\right)\left(\begin{array}{cc}{J}_i& C_i^T{L}_i+P_i{G}_i\\ *& L_i^T{L}_i-{\gamma }^{2}I\end{array}\right)\left(\begin{array}{c}{x}_t\\ {\omega }_t\end{array}\right)\mathrm{dt}\right\}-E\left[V(x_T,s_T)\right]+E\left[V(x_0,s_0)\right]---\left(26\right)$

由于E[V(x_T,s_T)]≥0以及式(15)，有：

J_T≤E[V(x₀,s₀)]，

当T→∞时，有：

${\left|\right|z\left|\right|}_2\le \gamma {[{\left|\right|\omega \left|\right|}_2^2+x_0^{T}P\left(s_0\right)x_0]}^{\frac{1}{2}},$

即由式(23)-(26)知，式(11)使(5)表示的时变电力系统具有扰动衰减度γ。

最终，系统(5)所示时变电力系统的稳定性可以采用如下优化问题表示：

$\underset{P_i,\gamma }{\min }t$               (27)

s.t.(11)(12)

利用可行性问题求解方法对LMI进行求解，若存在t＜0以及满足判据所示 正定矩阵P_i>0，则使(5)表示的时变电力系统鲁棒随机稳定，且满足扰动衰减度 γ；否则，时变电力系统不稳定。

因此，如图1所示，本发明包括一种时变电力系统稳定性分析系统，所述 系统包括依次连接的数据采集模块，风速拟合预测模块、时变系统建立与降阶 模块、稳定判据求解模块和结果输出模块，数据采集模块还连接至时变系统建 立与降阶模块；

所述数据采集模块用于采集风电场风速与风机出力信息、网络结构参数、 系统内发电机频率、功角，并将采集数据发送至风速拟合预测模块和时变系统 建立与降阶模块；

所述风速拟合预测模块根据数据采集模块发送的信息与风电场历史实际数 据，对短时间内条件风速趋势进行拟合预测，并将风速区间化，根据拟合模型 计算各区间的条件特征风速概率密度矩阵；

所述时变系统建立与降阶模块将每个区间作为时变系统的一个子工况，根 据风机风速与出力对应关系以及数据采集模块的采集数据建立考虑风速随机特 征的连续马尔科夫时变电力系统模型，并对系统进行降阶；

所述的稳定判据求解模块用于形成时变电力系统鲁棒随机稳定判据，利用 线性矩阵不等式中的可行性问题求解方法求解系统是否稳定；

所述的结果输出模块用于输出系统稳定性判别结果。

通过本发明的时变电力系统稳定性分析系统，能够有效解决风机接入下风 速随机特征导致的电力系统时变稳定性不易判别的难题。首先利用双参数韦伯 分布模型对短时间内条件风速趋势进行拟合预测，再将风速区间化，计算各区 间的条件特征风速概率密度矩阵，根据条件特征风速确定系统的运行工况。然 后，根据系统网架信息建立考虑风速随机特征的连续马尔科夫时变电力系统模 型，并降阶系统。

与之相应，本发明还公开了一种时变电力系统稳定性分析方法，所述方法 包括：

A、采集风电场风速与风机出力信息、网络结构参数、系统内发电机频率、 功角；

B、根据采集的数据与风电场历史实际数据，对短时间内条件风速趋势进行 拟合预测，并将风速区间化，根据拟合模型计算各区间的条件特征风速概率密 度矩阵；

C、将每个区间作为时变系统的一个子工况，根据风机风速与出力对应关系 以及采集的数据建立考虑风速随机特征的连续马尔科夫时变电力系统模型，并 对系统进行降阶；

D、形成时变电力系统鲁棒随机稳定判据，利用线性矩阵不等式中的可行性 问题求解方法求解系统是否稳定；

E、输出时变电力系统稳定分析结果。

在一个具体实施方式中，步骤B中所述的条件特征风速概率密度为某一条 件特征风速下，预定时间间隔后(例如采样间隔通常选为10min，所述预定时间 间隔为这样一个采样间隔或几个采样间隔之后)的风速概率密度。

特别地，步骤B中对短时间内条件风速趋势进行拟合预测，并将风速区间 化，根据拟合模型计算各区间的条件特征风速概率密度矩阵包括：

B1、利用双参数韦伯分布模型拟合预测时变系统各条件风速下的变化趋势；

B2、按照一定的区间长度将风速在切入风速与切出风速之间进行区间化， 利用步骤B1中各条件风速下的拟合模型确定各区间的条件特征风速概率密度矩 阵。

更具体地，一个具体实施方式中，所述双参数韦伯分布模型为：

$g\left(v_i\right|v_j)=\left(\frac{k}{c}\right){\left(\frac{v_i}{c}\right)}^{k-1}\exp [-{\left(\frac{v_i}{c}\right)}^k],$

其中，v_i为当前时刻风速；

v_j为条件风速；

k为形状系数；

c为尺度系数；

g(v_i|v_j)为条件风速为v_j时当前风速为v_i的概率密度。

另外，一个实施方式的步骤C中，根据风机风速与出力对应关系以及采集 的数据建立考虑风速随机特征的连续马尔科夫时变电力系统模型，并对系统进 行降阶包括：

C1、将每个区间作为时变系统的一个子工况，根据风机风速与出力的对应 关系确定各区间特征风速对应的风机输出功率，并利用该风机输出功率以及采 集的数据建立考虑风速随机特征的连续马尔科夫时变电力系统模型；

特别地，所述风机风速与出力对应关系为：

$P\left(v_i\right)=\left(\begin{array}{cc}0& 0\le v_i<v_{\mathrm{ci}}{v}_i\ge v_{\mathrm{co}}\\ \frac{v_i-v_{\mathrm{ci}}}{v_R-v_{\mathrm{ci}}}{P}_R& v_{\mathrm{ci}}<v_i<v_R\\ P_R& v_R\le v_i<v_{\mathrm{co}}\end{array}\right),$

其中，v_ci为风机切入风速；

v_co为风机切出风速；

v_R为风机额定风速；

P_R为风机额定输出功率。

另外，连续马尔科夫时变电力系统模型表达式具体为：

$\left(\begin{array}{c}{x}_t=(A\left(s_t\right)+\mathrm{\Delta A}(s_t,t))x_t+B\left(s_t\right)u_t+G\left(s_t\right){\omega }_t\\ z_t=C\left(s_t\right)x_t+D\left(s_t\right)u_t+L\left(s_t\right){\omega }_t\end{array}\right)$

其中，x_t∈Rⁿ是状态向量，

u_t∈R^p是控制输入向量，

t是时间，

ω_t为白噪声，满足：$E\left\{d{\omega }_t\right\}=0,E\left\{{\mathrm{d\omega }}_t^2\right\}=\mathrm{Wdt}(W>0),$

{s(t),t≥0}是在有限空间S＝{1,2,…,l}中取值的连续马尔科夫过程，对应各 个风速模式下的系统运行工况，描述系统工况随风机出力变化的演化过程。其 状态转移概率密度可以表示为：

$p_{\mathrm{ij}}=P({\eta }_{t+\Delta }=j|{\eta }_t=i)=\left(\begin{array}{cc}{\pi }_{\mathrm{ij}}\Delta +o\left(\Delta \right),& i\ne j\\ 1+{\pi }_{\mathrm{ii}}\Delta +o\left(\Delta \right),& i=j\end{array}\right)$

其中，$\underset{\Delta \to \infty }{\lim }\frac{o\left(\Delta \right)}{\Delta }=0,(\Delta >0),$

η_t为时刻t时的运行工况，Δ为时刻t的变化量，

o(Δ)是时刻t的变化量Δ的高阶无穷小量，

矩阵π表示马尔科夫转移概率密度矩阵，π_ij是系统运行工况在t时刻处于i， 而在t+Δ时刻处于j的转换概率密度，并且：

${\pi }_{\mathrm{ii}}=-\underset{i\ne j}{\Sigma }{\pi }_{\mathrm{ij}}({\pi }_{\mathrm{ij}}\ge 0,i\ne j),$

对于每一个s_t＝i∈S，记A(s_t)、B(s_t)、G(s_t)、C(s_t)、D(s_t)、L(s_t)分别表 示为A_i、B_i、G_i、C_i、D_i、L_i，并且不确定参数满足匹配性条件：

ΔA(s_t,t)＝ΔA(i,t)＝H_iF(i,t)M_i，

其中，H_i和M_i为已知矩阵，I为单位矩阵，实矩阵F(i,t)反映了系统不确 定参数的结构信息，满足：

F^T(i,t)F(i,t)≤I。

C2、在保证所关心频带输入输出特性不变的前提下对连续马尔科夫时变电 力系统模型降阶。

特别地，步骤C2中所述对系统进行降阶包括：

C21、将系统状态方程划分为：

$\left(\begin{array}{c}{\stackrel{·}{X}}_1\\ {\stackrel{·}{X}}_2\end{array}\right)=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)\left(\begin{array}{c}{X}_1\\ X_2\end{array}\right),$

其中，X₁＝[Δω^T,Δδ^T]为保留变量，Δω表示发电机功角转速的变化量，Δδ表 示发电机功角变化量，X₂为其他变量；

C22、消去X₂后为其中I为单位阵，p为 微分算子。

可以将上式改写为其中A_r(p)为运算形式的降阶系统系数 阵。

由上三式可以得到两个重要性质：

(1)如果p＝λ₁(i＝1,2,…,N)为第一式相应系统特征根，即|λ_iI-A|＝0,则 p＝λ_i也为第二式或第三式形式上降阶的系统特征根，即亦有|λ_iI-A_r(λ_i)|＝0， 特征根不发生变化，系统模式不变。

(2)对于原系统，λ_i的特征向量u_i，有Au_i＝λ_iu_i。设降阶系统特征根λ_i相应 的特征向量为u_ri，即A_r(λ_i)u_ri＝λ_iu_ri，则u_ri和u_i中保留变量X_r相对应的元素相 等，即特征向量的相应元素不变。因此，在X_r保留变量处去观察同一模式λ_i的 振荡时，相对幅值即相位不变，或者说模态不变。

因此，所关心的频带输入输出特性被完整的保留下来，同时实现了降阶的 效果。

另外在本发明一具体实施方式中，步骤D中，形成时变电力系统鲁棒随机 稳定判据，利用线性矩阵不等式中的可行性问题求解方法求解系统是否稳定包 括：

D1、构建含连续马尔科夫时变电力系统模型的李雅普诺夫泛函，并在该泛 函的弱无穷小算子中运用邓肯引理，推导满足干扰衰减度γ的鲁棒随机稳定性的 线性矩阵不等式，得到时变电力系统鲁棒随机稳定判据；

D2、利用线性矩阵不等式中的可行性问题求解方法求解系统是否稳定。

特别地，所述含连续马尔科夫时变电力系统模型的李雅普诺夫泛函表达式 为：

$V(x_t,s_t=i)=V(x_t,i)=x_t^T{P}_i{x}_i(P_i>0,i\in S),$

其中P_i为正定矩阵。

另外，所述的弱无穷小算子具体表达式为：

$\mathrm{lV}(x_t,i)=\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}V(x_t,i)+x_t^T{A}_i^T\frac{\partial }{\partial x}V(x_t,i)=x_t^T[\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+A_i^T{P}_j+P_j{A}_i]x_t<0,$

而所述的邓肯引理表述为：设函数f(x)的定义域为I：如果对于属于I内某 个区间上的任意两个自变量的值x₁、x₂，当x₁大于x₂时都有f(x₁)大于f(x₂).那 么就说f(x)在这个区间上是增函数。如果对于属于I内某个区间上的任意两个 自变量的值x₁、x₂，当x₁小于x₂时都有f(x₁)小于f(x₂)，那么就是f(x)在这个 区间上是减函数。

另外，所述鲁棒随机稳定性判据具体可以表示为：给定正常数γ>0，如果 存在一组正定对称矩阵P_i>0，i∈S，使得如下一组矩阵不等式成立：

$\left(\begin{array}{cc}{J}_i& C_i^T{L}_i+P_i{G}_i\\ *& L_i^T{L}_i-{\gamma }^{2}I\end{array}\right)<0,$

$J_i=A_i^T{P}_i+P_i{A}_i+\underset{j=1}{\overset{l}{\Sigma }}{\pi }_{\mathrm{ij}}{P}_j+C_i^T{C}_i,$

则当u(t)≡0时，时变系统鲁棒随机稳定，且满足扰动衰减度γ，即：

${\left|\right|z\left|\right|}_2\le \gamma {[{\left|\right|\omega \left|\right|}_2^2+x_0^{T}P\left(s_0\right)x_0]}^{\frac{1}{2}}.$

在本发明具体实施方式中，所述的线性矩阵不等式中可行性问题的求解方 法具体可以表述为：

对于给定的线性矩阵不等式系统A(x)＜B(x)，通过求解下式所示辅助凸优 化问题，寻找到全局最小标量值t_min，如果t_min＜0，则存在x∈Rⁿ或矩阵X∈Rⁿ满足下式：

min t

s.t.A(x)-B(x)≤tI。

以下通过两个具体示例来说明本发明时变电力系统稳定性分析系统及方法 的技术效果。

示例1

以某风电场1个月内实测风速数据为例进行分析，其中切入风速为3m/s， 额定风速为12m/s，切出风速为22m/s，风机输出功率曲线如图2所示。以1m/s 为区间宽度将风速划分为10个模式，构成连续马尔科夫模型状态集合，记作 S＝{1,2,…,10}，各模式信息在图3中给出。

根据图3中模式划分对各风速模式下条件风速概率密度曲线进行韦伯分布 拟合，拟合结果如图4所示。根据图4及式(1)求得各条件风速概率密度矩阵如 式(28)所示。

搭建如图5所示的IEEE4机11节点系统，该系统包含区域1与区域2两个 区域，其中发电机1和发电机2位于区域1，发电机3和发电机4位于区域2， 两个区域的联络线7-8与8-9均为双回线。发电机采用6阶详细模型，励磁系统 采用快速励磁，基准模型下的负荷采用50％恒阻抗和50％恒电流模型。分析时 将区域2中发电机G4以等容量的功率源代替用来模拟风机输出功率的工况变化 情况。

$\pi =\left(\begin{array}{cccccccccc}-0.2414& 0.2413& 0.0001& 0& 0& 0& 0& 0& 0& 0\\ 0.2061& -0.4255& 0.2184& 0.001& 0& 0& 0& 0& 0& 0\\ 0.0469& 0.2205& -0.5102& 0.2388& 0.004& 0& 0& 0& 0& 0\\ 0.0079& 0.0520& 0.2102& -0.5481& 0.2694& 0.0086& 0& 0& 0& 0\\ 0.0012& 0.0103& 0.0553& 0.202& -0.5736& 0.2901& 0.0147& 0& 0& 0\\ 0.0001& 0.0013& 0.01& 0.0523& 0.1953& -0.5678& 0.2975& 0.0113& 0& 0\\ 0& 0.0001& 0.0011& 0.0086& 0.0483& 0.1979& -0.5319& 0.2728& 0.0031& 0.0001\\ 0& 0& 0.0001& 0.0009& 0.0073& 0.045& 0.2057& -0.4779& 0.2189& 0.0003\\ 0& 0& 0& 0& 0.0005& 0.0046& 0.0357& 0.204& -0.3848& 0.14\\ 0& 0& 0& 0& 0& 0& 0& 0& 0.0049& -0.0049\end{array}\right)---\left(28\right)$

首先，利用模态分析法得到四机系统相应风机出力模式下的状态矩阵并利 用SMA方法降阶，将降阶后的系统状态矩阵与式(28)所示状态转移概率密度矩 阵代入式(27)，扰动衰减度γ取0.1，利用可行性问题求解方法求得该LMI成立， 且存在正定矩阵P_i，4机系统鲁棒随机稳定。其中，式(29)-(30)给出了部分正定 矩阵P_i的取值。

$P_1=\left(\begin{array}{cccccc}1.062& -0.945& 0.083& 2.552& 1.408& 10.81\\ -0.945& 1.162& -0.145& -2.065& 1.058& -13.30\\ 0.083& -0.145& 1.837& -12.01& 12.54& 0.772\\ 2.552& -2.065& -12.01& 4184.3& -938.6& -1333.5\\ 1.407& 1.058& 12.53& -938.6& 4012.6& -1766.2\\ 10.81& -13.31& 0.772& -1333.5& -1766.2& 6190.9\end{array}\right)---\left(29\right)$

$P_3=\left(\begin{array}{cccccc}1.117& -0.891& -0.012& 1.807& 0.635& 3.632\\ -0.891& 1.182& -0.171& -1.736& 0.886& -4.96\\ -0.012& -0.171& 2.047& -4.119& 4.671& 0.562\\ 1.807& -1.736& -4.119& 3588.8& -1091.1& -1054.9\\ 0.635& 0.886& 4.671& -1091.1& 4152.4& -1684.2\\ 3.632& -4.961& 0.562& -1054.9& -1684.2& 6430.9\end{array}\right)---\left(30\right)$

由上述求解结果可知，该4机系统满足鲁棒随机稳定判据所示条件，即该 系统在考虑风速随机特征时为鲁棒随机稳定的系统。图6和图7分别就系统实 际运行中可能出现的运行工况变化情况进行时域分析，其中，图6为系统仅发 生一次工况变化情况下，发电机G2与G3间的相对功角动态响应示意图，图 7为系统发生多工况变化情况时发电机G2与G3间的相对功角动态响应示意 图。图中系统实际分析结果与本发明时变电力系统稳定性分析方法的判断结果 一致。由此可知，本文提出的时变电力系统稳定性分析方法可以在不获取运行 轨迹的前提下准确判别系统稳定性，计算量少，具有较好的有效性和准确性。

示例2

搭建如图8所示的IEEE16机68节点系统，进一步考查考虑风速随机特 征的时变电力系统稳定性分析方法的有效性和通用性。该系统可分为5大区 域，其中区域1、2和3为等值系统，区域4为纽约系统，区域5为新英格兰 系统，将区域3中发电机G16以等容量的功率源代替用来模拟风机输出功率 的变化情况。发电机采用6阶详细模型，励磁采用IEEE-DC1型励磁,负荷模 型采用WECC负荷模型，80％的恒有功负荷，80％的恒无功阻抗负荷，20％的 动态负荷。

首先利用SMA方法分别对各工况下系统状态矩阵降阶，并将降阶后的系 统状态矩阵与式(28)所示概率密度转移矩阵代入式(27)，扰动衰减度γ取0.1， 利用可行性问题求解方法求得该LMI成立。但是，由于所得矩阵P_i中P₆存在 负数特征根-15.3827，即矩阵P_i存在非正定矩阵，因此，该16机系统在考虑风 机出力变化时不满足鲁棒随机稳定条件。

图9和图10分别就系统实际运行中可能出现的运行工况变化情况进行时 域分析，其中，图9为系统仅发生一次工况变化情况时发电机G1与G8间的 相对功角动态响应示意图，图10为多次工况变化情况时发电机G1与G8间的 相对功角动态响应示意图。从图中可知系统在单次与多次工况变化情况时均存 在不稳定情况，这与本方法的判稳结果一致。

分别对系统各运行工况的状态矩阵进行分析，各状态矩阵最小阻尼比如图 11所示。由图11可知，各运行模式下状态矩阵特征值均具有负实部，此时， 若利用特征值分析法单独分析任一运行工况稳定性，分析结果均为稳定。由此 可知，在考虑风速随机特征的时变电力系统中，系统从一个稳定工况切换到另 一个稳定工况时可能发生失稳，该现象是传统针对单一工况的稳定性分析方法 无法分析的，本文提出的时变电力系统稳定性分析方法能够有效针对这一现象 进行稳定性判别分析，方法简便，具有良好的准确性和有效性。

需要说明的是，上述实施方式仅为本发明较佳的实施方案，不能将其理解 为对本发明保护范围的限制，在未脱离本发明构思前提下，对本发明所做的任 何微小变化与修饰均属于本发明的保护范围。