CN102832617A - 一种基于精确化模式判别的大电网暂态稳定分析方法 - Google Patents

Abstract

本发明提供一种基于精确化模式判别的大电网暂态稳定分析方法，包括如下步骤：步骤S1：获得电网动态响应数据，其中，所述相关数据包括SCADA、PMU和网络拓扑结构数据；步骤S2：对大电网的动态响应数据进行特征提取进行特征提取；步骤3：利用DHMM方法，进行实时模式判别。本发明方法考虑到了电气特征量的动态响应特性，并且多维特征分析同时进行模式识别，提高了人工智能类评估方法的识别的精度和速度，有助于大电网暂态稳定分析工作的开展。

Description

一种基于精确化模式判别的大电网暂态稳定分析方法

技术领域

本发明属于电力工程领域，尤其涉及一种基于精确化模式判别的大电网暂态稳定分析方法。

背景技术

目前，常见的大电网暂态稳定分析方法有以下几种：

（1）时域仿真法：

该方法包括步骤：1）求解根据电力系统中发电机、负荷等各种元件间的拓扑关系所得到描述系统暂态行为的微分-代数方程组，得到状态变量和代数变量随时间变化的轨迹；2）根据研究时间内各发电机转子之间相对角度的变化情况，判断系统的稳定性。时域仿真法可以根据需要选择不同详细程度的元件数学模型，且结果准确、可靠。

（2）直接法：

直接法主要有两大类：基于暂态能量函数（Transient Energy Function,TEF）的方法和基于扩展等面积准则（Extended Equal Area Criterion,EEAC）的方法。

基于暂态能量函数方法的思想源于单机无穷大系统的“等面积定则”，暂态能量函数包含暂态动能和暂态势能两部分，其中暂态动能是大扰动造成发电机失去同步的能量，暂态势能由位置能量、磁能和耗散能量等构成。故障切除后，系统能否保持稳定运行取决于故障后系统吸收暂态能量的能力，若故障后暂态动能完全转化为暂态势能，则系统是暂态稳定的；否则，系统是不稳定的。

基于扩展等面积准则的方法假定系统失稳为双机模式，将发电机按一定的规则分为各自同步的两群，进行双机等值，并最终转化为单机无穷大系统，然后用“等面积准则”判断系统的暂态稳定性。

直接法具有计算速度快、能够定量分析系统的稳定度等优点。

（3）人工智能法：

人工神经网络（ANN）、支持向量机（SVM）、决策树，贝叶斯分类器等算法及其使用不同特征提取方法的改进算法，都是目前用于暂态稳定分析的主流方法。人工智能法能够通过离线少量数据学习、快速建立辨识模型、用于实际数据辨识等优点。

但是上述方法还存在如下问题：

时域仿真法存在的缺点有：计算速度慢，不能适应大规模电力系统的在线应用；自身不包含暂态稳定的判据，只能得到发电机功角等信息，通过发电机转子之间最大摆开角是否大于某一指定的阈值，判断系统遭受大干扰后能否保持稳定运行，然而关于阈值的取值至今仍无标准，通常需要靠经验确定；此外，无法定量地分析暂态稳定的程度。

直接法也存在一些较大的缺陷，如模型适应性差，目前采用直接法的大电网暂态稳定分析大多基于简化模型的假设，不能计及各种复杂动态元件对系统稳定性的影响；只能分析系统的第一摆动稳定性，不适用于更复杂的失稳模式；分析结果的可靠性难以保证。

而人工智能法很少涉及大电网从受到扰动到建立新平衡的动态演变过程，而且要求的训练样本多，运算时间长。

发明内容

针对现有技术中大电网暂态稳定分析方法存在的计算速度慢、分析结果不够准确的问题，本发明提出了一种充分考虑了电气特征量的动态响应特性的、基于精确化模式判别的大电网暂态稳定分析方法，该方法能更快速、更精确的进行大电网暂态稳定分析。

为了解决上述技术问题，本发明采用如下技术方案：

一种基于精确化模式判别的大电网暂态稳定分析方法，包括步骤：

步骤S1、采集不同时刻的大电网的电气特征量，并根据电气特征量获得大电网的动态响应数据；

步骤S2、对大电网的动态响应数据进行特征提取；

步骤S3、利用离散隐马尔科夫模型对进行了特征提取的大电网动态响应数据进行实时模式判别，即可得到大电网暂态稳定分析结果。

上述步骤S1进一步包括以下子步骤：

步骤S1-1、采集大电网的电气特征量，并将所采集的大电网电气特征量进行整合；

步骤S1-2、采用相对灵敏度的方法筛选出相对灵敏度较高的电气特征量构成电气特征量子集，所述的相对灵敏度较高的电气特征量为相对灵敏度高于预设值m的电气特征量、或相对灵敏度最大的n个电气特征量、或相对灵敏度既高于预设值m的电气特征量中最大的n个电气特征量，所述的m、n为预先设置的自然数；

步骤S1-3、对电气特征量子集中的电气特征量进行归一化处理，即得到大电网的动态响应数据。

上述步骤S2利用PCA方法对大电网的动态响应数据进行特征提取。

上诉步骤S3进一步包括数据训练和大电网暂态稳定分析两个步骤，所述的数据训练进一步包括以下子步骤：

步骤S3a-1、对作为训练数据的电气特征量子集中的电气特征向量进行基于时间的离散，得到离散观测值样本，所述的电气特征量子集训练数据包括失稳的大电网电气特征量子集样本数据和稳定的大电网电气特征量子集样本数据；

步骤S3a-2、对离散观测值样本进行训练得到失稳的训练样本集和稳定的训练样本集；

所述的大电网暂态稳定分析进一步包括以下子步骤：

步骤S3b-1、对待分析稳定性的电气特征量子集中的电气特征向量进行基于时间的离散，得到离散观测值样本；

步骤S3b-2、将步骤S3b-1所得离散观测值样本与步骤S3a-2所得的失稳的训练样本集和稳定的训练样本集进行比对，从而得到大电网暂态稳定性的分析结果。

上述步骤S3a-1中是采用基于劳合社的算法对电气特征量子集中的电气特征向量进行基于时间的离散。

上述步骤S3a-2中采用离散隐马尔科夫模型对离散观测值样本进行训练，具体为：采用如下公式（1）和（2）对离散观测值样本进行样本训练：

${\stackrel{\‾}{a}}_{\mathrm{ij}}=\frac{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{t=1}{\overset{T_l-1}{\mathrm{\Σ}}}{{\mathrm{\α}}_t}^{*\left(l\right)}\left(i\right)a_{\mathrm{ij}}{b}_j\left(O_{t+1}^{\left(l\right)}\right){\mathrm{\β}}_{t+1}^{*\left(l\right)}\left(j\right)/P(O^{\left(l\right)}/\mathrm{\λ})}{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\underset{t=1}{\overset{T_l-1}{\mathrm{\Σ}}}{{\mathrm{\α}}_t}^{*\left(l\right)}\left(i\right){\mathrm{\β}}_t^{*\left(l\right)}\left(j\right)/P(O^{\left(l\right)}/\mathrm{\λ})}---\left(1\right)$

1≤i≤N,1≤j≤N

1≤i≤N,1≤k≤M

其中：

为训练后的状态转移概率矩阵A中第i行j列的元素，a_ij为训练前的状态转移概率矩阵A中第i行j列的元素；

为训练后的观测值概率矩阵B中第j行k列的元素；

α_t ^*(l)(i)、α_t ^*(l)(j)分别为改进后的第l个训练序列中第i、j行的元素在t时刻的前向变量， ${{\mathrm{\α}}_t}^{*\left(l\right)}\left(i\right)={\mathrm{\α}}_t^{\left(l\right)}\left(i\right)/{\mathrm{\Φ}}_1{\mathrm{\Φ}}_2\·\·\·{\mathrm{\Φ}}_t,$ ${{\mathrm{\α}}_t}^{*\left(l\right)}\left(j\right)={\mathrm{\α}}_t^{\left(l\right)}\left(i\right)/{\mathrm{\Φ}}_1{\mathrm{\Φ}}_2\·\·\·{\mathrm{\Φ}}_t,$ ${\mathrm{\Φ}}_m={\∑}_{i=1}^N{\mathrm{\α}}_m^{\left(l\right)}\left(i\right),$ m＝1,2,...,t，

分别为第l个训练序列中第i、j、m行的元素在t时刻的前向变量；

分别为改进后的第l个训练序列中第j行的元素在t、t+1时刻的后向变量， ${\mathrm{\β}}_t^{*\left(l\right)}\left(i\right)=\frac{{\widetilde{\mathrm{\β}}}_t^{\left(l\right)}\left(i\right)}{{\mathrm{\Φ}}_{t+1}},$ ${\mathrm{\Φ}}_{t+1}={\∑}_{i=1}^N{\mathrm{\α}}_{t+1}^{\left(l\right)}\left(i\right),$ ${\mathrm{\β}}_t^{*\left(l\right)}\left(i\right)={\∑}_{j=1}^N{a}_{\mathrm{ij}}{b}_j\left(O_{t+1}^{\left(l\right)}\right){\mathrm{\β}}_{t+1}^{*\left(l\right)}\left(j\right),$ t＝1,2,…,T-1，β_T ^*(l)(i)＝1，

为第l个训练序列中第i行的元素在t+1时刻的后向变量，β_T ^*(l)(i)为改进后的第l个训练序列中第j行的元素在T时刻的后向变量，T表示观测结束的时间；

O^(l)为第l个离散观测样本序列，为在t+1时刻的第l个离散观测样本序列，O_t为t时刻的观测值，V_k为选定进行样本训练的观测值；

为t+1时刻观测值

的概率密度函数；

P(O^(l)/λ)为在模型λ下、观测值O^(l)的条件概率；

t表示观测时刻，T_l是电气特征值序列在时间维度上的采样长度，N是马尔科夫链状态的数目，M是每个状态对应可能的观测值的个数，L为离散观测样本的维数。

上述步骤S3b-2是采用离散隐马尔科夫模型比较步骤S3b-1所得离散观测值样本与步骤S3a-2所得的失稳的训练样本集和稳定的训练样本集。

与现有技术相比，本发明具有以下优点和有益效果：

本发明方法充分考虑了电气特征量的动态响应特性，并且能够对多维特征分析同时进行模式识别，提高了人工智能类评估方法的识别的精度和速度，有助于大电网暂态稳定分析工作的开展。

附图说明

图1是本发明的大电网暂态稳定分析的精确化模式判别方法的流程图；

图2是5状态左右型DHMM；

图3是重庆电网示意图；

图4是重庆电网系统DHMM迭代曲线。

具体实施方式

本发明提供了一种基于精确化模式判别的大电网暂态稳定分析方法，本方法充分考虑到了大电网电气特征量的动态响应特性，在对大电网的动态响应数据进行特征提取后，利用DHMM方法对进行了特征提取的大电网动态响应数据进行实时模式判别，即科学而精确的进行了大电网暂态稳定分析。本发明方法可解决现有大电网暂态稳定分析方法存在的计算速度慢、分析结果不够准确的问题。

下面将对本发明方法做进一步说明图1，本发明包括如下步骤：

步骤S1、采集不同时刻的大电网的电气特征量，并根据不同时刻的电气特征量获得大电网的动态响应数据。

该步骤的具体实施方式为：采集大电网不同时刻的数据采集与监视控制系统（SupervisoryControlAnd DataAcquisition，SCADA）数据、同步相量测量单元（Phasor Measurement Unit，PMU）数据和网络拓扑结构数据，并将所采集的电气特征量数据进行整合，得到大电网的动态响应数据。

步骤S2、利用主成分分析（Principal ComponentAnalysis，PCA）方法，对大电网的动态响应数据进行特征提取。

主成分分析法(PCA)是多元统计学中一种常用的降维方法。它的基本思想是用较少的几个不相关的新变量，代替原有较多的相关变量，并且新变量是原有变量的线性组合。所选取的新变量即被称为主成分，选取的原则是尽可能保留原有变量中所包含的信息。从统计学的角度分析，一个变量所含有的信息可用其方差来表征。如果某个变量的方差越大，那么它所包含的信息量也就越大。

设有k个p维随机矢量，其中任一p维随机矢量X_i表示为X_i=(x_1i，x_2i，…x_pi)^T，X_i的主成分为Y_i，i＝1,2,…,n，且n≤p，该主成分满足如下条件：

1)

其中a_i为p×1维数字矢量；

2)a_i满足

且a_i的取值原则是使得Y_i的方差达到最大；

3)各主成分Y₁、Y₂、…、Y之间是相互独立的。

将上述主成分原则推广到含n个样本的样本集X_p×n，X_p×n＝(X₁,X₂,...,X_n)，X_i=(x_1i，x_2i，…x_pi)^T，i＝1,2,…,n，且n≤p，设样本集X_p×n的协方差矩阵为S_p×n，S_p×n＝(S₁,S₂,...,S_n)，S_i＝(s_1i，s_2i，...,s_pi)^T，由公式（1）求协方差矩阵S_p×n：

$S_{p\×n}=\frac{\underset{i=1,u=1}{\overset{i=n,u=p}{\mathrm{\Σ}}}(x_{\mathrm{iu}}-{\stackrel{\‾}{x}}_i){(x_{\mathrm{iu}}-\stackrel{\‾}{x})}^T}{n}---\left(1\right)$

其中，

x_iu为X_i中元素，i＝1,2,…,n，u＝1,2,...,p，

为X_i中各元素的算术平均值，即

设S_p×n的特征根分别为λ₁≥…≥λ_n≥0，则S_p×n的第i个主成分的方差就是λ_i。设上述特征根对应的正则正交特征向量为a₁、a₂、…、a_n，令a_i＝(a_1i，…,a_ni)^T，则

为样本集X_p×n的第i个样本主成分。

是主成分Y_k的贡献率，它表示主成分Y_k所保留的原样本集的离散程度的信息的比例，第一主成分的贡献率最大，表明它解释原始样本集的能力最强；

是前k个主成分Y₁、Y₂、...、Y_k的累积贡献率，它表示前k个主成分Y₁、Y₂、...、Y_k保留的原样本集的离散程度信息的比例。

通常取累积贡献率达到一个比较高的百分比，例如85%～95%，此时用主成分来代替原来的样本集，降维的效果较好，信息的损失量也较少。本发明中，大电网的动态响应数据即为样本集X_p×n，采用上述方法对电网的动态响应数据进行特征提取。

步骤S3、利用离散隐马尔科夫模型(Discrete Hidden Markov Model，DHMM)对进行了特征提取的大电网动态响应数据进行实时模式判别，即可得到大电网稳定状态的判断结果。

为了更加清楚本发明的基于精确化模式判别的大电网暂态稳定分析方法，下面将结合较佳实施例作更详细的阐述。

首先，将采集的SCADA数据、PMU数据和网络拓扑结构数据（即电气特征量）进行整合，得到大电网的动态响应数据。

从采集的电气特征量中选取59维电气特征量作为原始电气特征量。为了对比各电气特征量在不同扰动下的变化趋势，采用相对灵敏度的方法筛选出对电网变化灵敏度高的电气特征量子集，具体为：选定发生变化前的电气特征量作为基准值，记为W₀，发生变化后的电气特征量为W₁，则相对灵敏度W_x见式（2）：

$W_x=\frac{\mathrm{\ΔW}}{W_0}\×100\%=\frac{W_1-W_0}{W_0}\×100\%---\left(2\right)$

依据相对灵敏度的大小对59维电气特征量进行升序排列，选择相对灵敏度绝对值大于等于100%的电气特征量中的前32维电气特征量作为电气特征量子集，并对电气特征量子集中的特征量进行归一化处理。

本具体实施中电气特征量子集中的32维电气特征量如表1所示。

表1电气特征量子集

然后，利用PCA方法对相对灵敏度较高的32维电气特征量进行特征提取。

PCA方法是通过降维的思想来精简变量，将多个相关的原始变量指标转化为几个独立的综合指标。主成分是原始变量指标的线性组合，其转换的理念是让原始变量指标的线性组合的变异达到最大，这样就可以尽量减少信息的损失，仅利用几个重要的主成分就解释原始数据的大部分变异，简化问题，更好的揭示电力系统变量之间的规律，提高分析效率。

最后，利用DHMM方法对经特征提取后的大电网动态响应数据进行实时模式判别，参见图1，具体包括如下过程：

（a）基于劳合社的算法（Lloyd's algorithm）的标量量化编码；

在特征提取过程中获得的32维电气特征量是关于时间连续变化的特征向量，为了方便DHMM分类器的训练和识别，需要对上述特征向量进行“标量化”，即通过一定编码方式，实现关于时间连续变化的特征向量的离散化。

劳合社的算法是通信领域的成熟信源编码技术，其基本思想是通过试探方式，减少输入信号的量化失真以获得最佳的分区参数，对本发明而言，这里的输入信号是电气特征值。该算法能先根据信号幅值将电气特征量分为L-1个区间，将电气特征量按信号幅值的大小进行升序排列，然后将信号的L个相邻区间映射为L个离散值，幅值为x的信号值，其对应的离散值index(x)的定义见式（3）：

$\mathrm{index}\left(x\right)=\left\{\begin{array}{cc}1& x\&le;\mathrm{patition}\left(1\right)\\ i& \mathrm{patition}\left(i\right)<x\&le;\mathrm{patition}(i+1)\\ N& \mathrm{patition}\left(N\right)<x\end{array}\right.---\left(3\right)$

其中，

patition(i)为第i个分区的幅值，i为自然数，1≤i≤L。

（b）数据训练

DHMM是马尔科夫链的一种形态；它包括了两个随机过程：一是马尔科夫链，描述状态转移的统计规律；二是观测矩阵，描述状态与观测值的统计关系。因各个状态是隐含的，所以只能通过观测值来估计整个随机过程。

DHMM常用式（4）描述，简记为式（5）：

λ＝(N,M,π,A,B)（4）

λ＝(π,A,B)（5）

其中，N是马尔科夫链状态的数目；M是每个状态对应可能的观测值的个数；π是初始概率分布；A是状态转移概率矩阵；B是观测值概率矩阵。经过实用化推导，其训练公式如式（6）和（7）所示：

${\stackrel{\&OverBar;}{a}}_{\mathrm{ij}}=\frac{\underset{t=1}{\overset{T-1}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right)a_{\mathrm{ij}}{b}_j\left(O_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)/P(O/\mathrm{\&lambda;})}{\underset{t=1}{\overset{T_l-1}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_t\left(i\right){\mathrm{\&beta;}}_t\left(j\right)/P(O/\mathrm{\&lambda;})}(1\&le;i\&le;N,1\&le;j\&le;N)---\left(6\right)$

其中，

为训练后的状态转移概率矩阵A中第i行j列的元素；

a_ij为训练前的状态转移概率矩阵A中第i行j列的元素；

α_t(i)、α_t(j)分别为训练前的状态转移概率矩阵A中第i、j行的元素在t时刻的前向变量；

b_j(O_t+1)为t+1时刻观测值O_t+1的概率密度函数；

β_t(j)、β_t+1(j)为训练前的观测值概率矩阵B中第j行的元素在t、t+1时刻的后向变量；

为训练后的观测值概率矩阵B中第j行k列的元素，即训练后在j状态观测到k所对应观测值的概率；

P(O/λ)为在模型λ下、观测值为O的条件概率；

O_t为t时刻的观测值；

V_k为选定进行样本训练的观测值；

t表示观测时刻；

T表示观测结束的时间；

T_l是电气特征值序列在时间维度上的采样长度；

N是马尔科夫链状态的数目；

M是每个状态对应可能的观测值的个数。

上述步骤（a）所获得的32维离散观测样本index(x)不能直接用于DHMM分类器进行训练。因为DHMM算法的本质是求解极大似然条件概率，在多次迭代过程中，多个0到1之间的概率值反复做乘法运算，容易造成向下溢出。为了避免向下溢出，本发明对现有的DHMM训练方法进行了改进。

将L维离散观测样本index(x)序列表示为如式（8）所示：

O＝[O⁽¹⁾,O⁽²⁾,…O^(L)]（8）

其中，

表示第l个离散观测样本序列，即第l个经劳合社算法处理后的电气特征量序列，l为自然数，且l∈[1,L]，T_l是电气特征值序列在时间维度上的采样长度，本具体实施中，L=32。

为防止迭代时出现向下溢出问题，本发明对现有的DHMM训练公式进行改进，具体为：

1）对式（6）和式（7）中前向变量α_t(i)进行如下处理：

α₁(i)＝π_ib_i(o₁),1≤i≤N    （9）

${{\mathrm{\&alpha;}}_1}^{*}=\frac{{\mathrm{\&alpha;}}_1\left(i\right)}{{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&alpha;}}_1\left(i\right)}\&ap;\frac{{\mathrm{\&alpha;}}_1\left(i\right)}{{\mathrm{\&Phi;}}_1},1\&le;i\&le;N---\left(10\right)$

${\widetilde{\mathrm{\&alpha;}}}_{t+1}\left(j\right)=\left[{\&Sum;}_{i=1}^N{\mathrm{\&alpha;}}_t^{*}\left(i\right)a_{\mathrm{ij}}\right]b_j\left(o_{t+1}\right),1\&le;j\&le;N,t=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,T-1---\left(11\right)$

${\mathrm{\&alpha;}}_{t+1}^{*}\left(j\right)=\frac{{\widetilde{\mathrm{\&alpha;}}}_{t+1}\left(j\right)}{{\&Sum;}_{j=1}^N{\widetilde{\mathrm{\&alpha;}}}_{t+1}\left(j\right)}\&ap;\frac{{\widetilde{\mathrm{\&alpha;}}}_{t+1}\left(j\right)}{{\mathrm{\&Phi;}}_{t+1}},1\&le;j\&le;N,t=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,T-1---\left(12\right)$

其中，

α₁(i)为训练前的状态转移概率矩阵A中第i行的元素在1时刻的前向变量；即，训练前的状态转移概率矩阵A中第i行的元素初始时刻的前向变量；

π₁为状态转移概率矩阵A中第i行的元素的初始概率分布；

b_i(o_t)为训练前的状态转移概率矩阵A中第i行的元素的t时刻观测值o_t的概率密度函数；

α₁ ^*(i)是对α₁(i)进行改进后的初始时刻前向变量，即，改进后状态转移概率矩阵A中第i行的元素初始时刻的前向变量；

为对α_t+1(i)进行改进后的t+1时刻的前向变量，即，改进后状态转移概率矩阵A中第i行的元素在t+1时刻的前向变量；

a_ij为训练前的状态转移概率矩阵A中第i行j列的元素；

b_j(o_t+1)为t+1时刻观测值O_t+1的概率密度函数；

Ф_t定义为 ${\mathrm{\&Phi;}}_t={\&Sum;}_{i=1}^N{\mathrm{\&alpha;}}_t\left(i\right),$ 则， ${\mathrm{\&Phi;}}_1={\&Sum;}_{i=1}^N{\mathrm{\&alpha;}}_1\left(i\right),$ ${\mathrm{\&Phi;}}_{t+1}={\&Sum;}_{i=1}^N{\mathrm{\&alpha;}}_{t+1}\left(i\right);$

定义为 ${\widetilde{\mathrm{\&alpha;}}}_{t+1}\left(j\right)=\left[{\&Sum;}_{i=1}^N{\mathrm{\&alpha;}}_t^{*}\left(i\right)a_{\mathrm{ij}}\right]b_j\left(o_{t+1}\right);$

t表示观测时刻；

T表示观测结束的时间；

N是马尔科夫链状态的数目。

2）对式（6）和式（7）中后向变量β_t(j)进行如下处理：

β_T(i)＝1,1≤i≤N    （13）

β_T ^*(i)＝1,1≤i≤N    （14）

${\widetilde{\mathrm{\&beta;}}}_t\left(i\right)={\&Sum;}_{j=1}^N{a}_{\mathrm{ij}}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}^{*}\left(j\right),1\&le;i\&le;N,t=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,T-1---\left(15\right)$

${\mathrm{\&beta;}}_t^{*}\left(i\right)=\frac{{\widetilde{\mathrm{\&beta;}}}_t\left(i\right)}{{\mathrm{\&Phi;}}_{t+1}},1\&le;i\&le;N,t=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,T-1---\left(16\right)$

其中：

β_T(i)为训练前的观测值概率矩阵B中第i行的元素在T时刻的后向变量；

β_T ^*(i)是对β_T(i)进行改进后的T时刻的后向变量，即，改进后观测值概率矩阵B中第i行的元素在T时刻的后向变量；

是对β_t(i)进行改进后的t时刻的后向变量，即，改进后观测值概率矩阵B中第i行元素在t时刻的后向变量；

为自定义的变量，其定义为

a_ij为训练前的状态转移概率矩阵A中第i行j列的元素；

b_j(o_t+1)为t+1时刻观测值O_t+1的概率密度函数；

是对β_t+1(j)进行改进后的t+1时刻的后向变量，即，改进后观测值概率矩阵B中第j行元素在t+1时刻的后向变量；

Ф_t+1定义为 ${\mathrm{\&Phi;}}_{t+1}={\&Sum;}_{i=1}^N{\mathrm{\&alpha;}}_{t+1}\left(i\right);$

t表示观测时刻；

T表示观测结束的时间；

N是马尔科夫链状态的数目。

3）对式（6）和式（7）中条件概率P(O/λ)进行如下处理：

对α和β做了上述处理之后，为了保持原有公式（6）和（7）的计算结果不变，必须在条件概率P(O/λ)的计算公式中做相应的处理，以消除比例因子的影响。由公式（9）～（12）可以推出：

${\mathrm{\&alpha;}}_t^{*}\left(i\right)={\mathrm{\&alpha;}}_t\left(i\right)/{\mathrm{\&Phi;}}_1{\mathrm{\&Phi;}}_2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&Phi;}}_t---\left(17\right)$

${\mathrm{\&Phi;}}_t={\&Sum;}_{j=1}^N{\widetilde{\mathrm{\&alpha;}}}_t\left(j\right)={\&Sum;}_{j=1}^N\left[{\&Sum;}_{i=1}^N{\mathrm{\&alpha;}}_{t-1}^{*}\left(i\right)a_{\mathrm{ij}}\right]b_j\left(o_t\right)={\&Sum;}_{j=1}^N{\mathrm{\&alpha;}}_t\left(j\right)/{\mathrm{\&Phi;}}_1{\mathrm{\&Phi;}}_2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&Phi;}}_{t-1}---\left(18\right)$

因此

${\&Sum;}_{j=1}^N{\mathrm{\&alpha;}}_t\left(j\right)={\mathrm{\&Phi;}}_1{\mathrm{\&Phi;}}_2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&Phi;}}_t---\left(19\right)$

即

$P(O/\mathrm{\&lambda;})={\&Sum;}_{j=1}^N{\mathrm{\&alpha;}}_T\left(j\right)={\mathrm{\&Phi;}}_1{\mathrm{\&Phi;}}_2\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&Phi;}}_T---\left(20\right)$

对公式（20）取对数，可得：

$\mathrm{lgP}(O/\mathrm{\&lambda;})={\&Sum;}_{j=1}^N\lg {\mathrm{\&Phi;}}_t---\left(21\right)$

采用下列式（22）和（23）对离散观测样本index(x)进行样本训练。

${\stackrel{\&OverBar;}{a}}_{\mathrm{ij}}=\frac{\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\underset{t=1}{\overset{T_l-1}{\mathrm{\&Sigma;}}}{{\mathrm{\&alpha;}}_t}^{*\left(l\right)}\left(i\right)a_{\mathrm{ij}}{b}_j\left(O_{t+1}^{\left(l\right)}\right){\mathrm{\&beta;}}_{t+1}^{*\left(l\right)}\left(j\right)/P(O^{\left(l\right)}/\mathrm{\&lambda;})}{\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\underset{t=1}{\overset{T_l-1}{\mathrm{\&Sigma;}}}{{\mathrm{\&alpha;}}_t}^{*\left(l\right)}\left(i\right){\mathrm{\&beta;}}_t^{*\left(l\right)}\left(j\right)/P(O^{\left(l\right)}/\mathrm{\&lambda;})}---\left(22\right)$

1≤i≤N,1≤j≤N

1≤i≤N,1≤k≤M

其中，

为训练后的状态转移概率矩阵A中第i行j列的元素；

a_ij为训练前的状态转移概率矩阵A中第i行j列的元素；

为训练后的观测值概率矩阵B中第j行k列的元素，即训练后在j状态观测到k所对应观测值的概率；

α_t ^*(l)(i)、α_t ^*(l)(j)分别为改进后的第l个训练序列中第i、j行的元素在t时刻的前向变量；

分别为改进后的第l个训练序列中第j行的元素在t、t+1时刻的后向变量；

为t+1时刻观测值

的概率密度函数；

P(O^(l)/λ)为在模型λ下、观测值O^(l)的条件概率；

O_t为t时刻的观测值；

V_k为选定进行样本训练的观测值；

t表示观测时刻；

T_l是电气特征值序列在时间维度上的采样长度；

N是马尔科夫链状态的数目；

M是每个状态对应可能的观测值的个数；

L为离散观测样本的维数，在本具体实施中，L=32。

在进行基于DHMM的模式识别前，必须进行样本训练。用来进行样本训练的训练数据包括稳定的样本集和失稳的样本集，训练数据采用步骤S1和步骤S2方法获取并进行特征提取，并进行关于时间连续的离散化，然后按照公式（22）、（23）进行训练，具体过程为：

在DHMM分类器中，每一个独立的DHMM自动进行如下训练：

1）DHMM模型初始化概率分布π、状态转移概率矩阵A和观测值概率矩阵B，即可到第1次迭代的条件概率P(O|λ₁)；

2）使用Baum-Welch（BW）算法计算得到DHMM模型的新参数，即状态转移概率矩阵A和观测值该矩阵B；

3）使用DHMM模型的新参数，依据Viterbi算法计算第j次迭代的最优条件概率P(O |λ_j)；

4）将第j次迭代的条件概率P(O|λ_j)同第j-1次迭代的条件概率P(O|λ_j-1)比较，判断是否满足下面收敛条件|P(O|λ_j)-P(O|λ_j-1)|≤ε，其中，ε为设定的收敛门限值；

5）满足收敛条件，停止迭代，输出该模型参数；如果不收敛，继续迭代，直到收敛或者达到设定的最大迭代次数。

（c）DHMM模型初始化

电力系统从遭受扰动到暂态稳定过程按时间顺序可大致分为5个阶段，即正常状态、扰动初期、发生阶段、发展阶段和恢复阶段，按照时间顺序对应DHMM的5个马尔科夫（Markov）状态，即N＝5，状态转移概率矩阵A为5×5的方阵；观测值矩阵O是一个32维的特征值离散序列。该DHMM模型可视为“左右型”DHMM模型，如图2所示。

初始概率分布π和状态转移概率矩阵A初始化如式（10）和（11）所示，初始化观测值概率矩阵B满足在幅值上服从正态分布即可，因为在后续的多次迭代过程中会得到较好的训练修正值。

π＝[1 0 0 0 0]（10）

$A=\left[\begin{array}{ccccc}0.5& 0.5& 0& 0& 0\\ 0& 0.5& 0.5& 0& 0\\ 0& 0& 0.5& 0.5& 0\\ 0& 0& 0& 0.5& 0.5\\ 0& 0& 0& 0& 1\end{array}\right]---\left(11\right)$

（d）基于DHMM的模式识别过程

原始特征序列（即本具体实施中原始获取的59维电气特征量）经过基于相对灵敏度的筛选后，成为降维后的电气特征量子集（即表1中的32维电气特征量）；该电气特征量子集经过归一化和基于劳合社的算法的标量量化之后输入DHMM分类器，可以输入训练数据进行数据训练，训练数据包括稳定的样本集和失稳的样本集，也可以输入待辨识稳定性的数据进行稳定性辨识；DHMM分类器将输入的待辨识稳定性的数据与训练数据进行比对，并输出辨识结果。

DHMM分类器输出结果的判断准则是，计算得到稳定DHMM的概率大于失稳的概率，辨识结果为系统稳定；反之为系统失稳。

下面列举一个应用实施例来进一步说明本发明方法的优点和有益效果。

重庆电网实际系统的网架结构如图3所示，采用2010年实际运行的数据模型，在省级电网A与B之间的联络线：L1-L2和L12-L13为500kV重载线路，一旦发生短路故障，容易引起系统暂态失稳。以上述四条联络线为研究对象，为获得更多的样本，将基准负荷80%-130%细化递增步长，但要均匀分布；发电机有功出力和机端电压也是80%-130%细化递增步长，且均匀分布。上述四条联络线故障均得到500个数据样本，共计2000个数据样本，将其中1000个数据样本用作训练样本，另1000个数据样本用作待辨识稳定性的数据。采用DHMM模型对1000个训练样本进行训练收敛迭代曲线如图4所示。

经过35次迭代，趋于平稳，达到最大的条件概率。以该训练模型的参数，按照图1所示流程进行暂态稳定评估。为了进行比较，使用ANN算法对上述算例进行仿真，在训练之前采用PCA方法将32维特征向量压缩到10维，辨识结果如表2所示。

表2某实际电网系统DHMM与ANN算法的性能比较

由表2所示，DHMM和ANN算法都能达到相当高的辨识正确率，但DHMM需要的训练样本更少，收敛更加迅速。

Claims (7)

1.一种基于精确化模式判别的大电网暂态稳定分析方法，其特征在于，包括步骤：

步骤S1、采集不同时刻的大电网的电气特征量，并根据电气特征量获得大电网的动态响应数据；

步骤S2、对大电网的动态响应数据进行特征提取；

步骤S3、利用离散隐马尔科夫模型对进行了特征提取的大电网动态响应数据进行实时模式判别，即可得到大电网暂态稳定分析结果。

2.根据权利要求1所述的基于精确化模式判别的大电网暂态稳定分析方法，其特征在于：

所述的步骤S1进一步包括以下子步骤：

步骤S1-1、采集大电网的电气特征量，并将所采集的大电网电气特征量进行整合；

步骤S1-2、采用相对灵敏度的方法筛选出相对灵敏度较高的电气特征量构成电气特征量子集，所述的相对灵敏度较高的电气特征量为相对灵敏度高于预设值m的电气特征量、或相对灵敏度最大的n个电气特征量、或相对灵敏度既高于预设值m的电气特征量中最大的n个电气特征量，所述的m、n为预先设置的自然数；

步骤S1-3、对电气特征量子集中的电气特征量进行归一化处理，即得到大电网的动态响应数据。

3.根据权利要求1所述的基于精确化模式判别的大电网暂态稳定分析方法，其特征在于：

所述的步骤S2利用PCA方法对大电网的动态响应数据进行特征提取。

4.根据权利要求1所述的基于精确化模式判别的大电网暂态稳定分析方法，其特征在于：

所述的步骤S3进一步包括数据训练和大电网暂态稳定分析两个步骤，所述的数据训练进一步包括以下子步骤：

步骤S3a-1、对作为训练数据的电气特征量子集中的电气特征向量进行基于时间的离散，得到离散观测值样本，所述的电气特征量子集训练数据包括失稳的大电网电气特征量子集样本数据和稳定的大电网电气特征量子集样本数据；

步骤S3a-2、对离散观测值样本进行训练得到失稳的训练样本集和稳定的训练样本集；

所述的大电网暂态稳定分析进一步包括以下子步骤：

步骤S3b-1、对待分析稳定性的电气特征量子集中的电气特征向量进行基于时间的离散，得到离散观测值样本；

步骤S3b-2、将步骤S3b-1所得离散观测值样本与步骤S3a-2所得的失稳的训练样本集和稳定的训练样本集进行比对，从而得到大电网暂态稳定性的分析结果。

5.根据权利要求4所述的基于精确化模式判别的大电网暂态稳定分析方法，其特征在于：

所述的步骤S3a-1中是采用基于劳合社的算法对电气特征量子集中的电气特征向量进行 基于时间的离散。

6.根据权利要求4所述的基于精确化模式判别的大电网暂态稳定分析方法，其特征在于：

所述的步骤S3a-2中采用离散隐马尔科夫模型对离散观测值样本进行训练，具体为：采用如下公式（1）和（2）对离散观测值样本进行样本训练：

1≤i≤N,1≤j≤N

1≤i≤N,1≤k≤M

其中：

为训练后的状态转移概率矩阵A中第i行j列的元素，a_ij为训练前的状态转移概率矩阵A中第i行j列的元素；

为训练后的观测值概率矩阵B中第j行k列的元素；

α_t ^*(l)(i)、α_t ^*(l)(j)分别为改进后的第l个训练序列中第i、j行的元素在t时刻的前向变量， 

m＝1,2,...,t， 

分别为第l个训练序列中第i、j、m行的元素在t时刻的前向变量；

分别为改进后的第l个训练序列中第j行的元素在t、t+1时刻的后向变量， 

t＝1,2,…,T-1，β_T ^*(l)(i)＝1， 

为第l个训练序列中第i行的元素在t+1时刻的后向变量，β_T ^*(l)(i)为改进后的第l个训练序列中第j行的元素在T时刻的后向变量，T表示观测结束的时间；

O^(l)为第l个离散观测样本序列， 

为在t+1时刻的第l个离散观测样本序列，O_t为t时刻的观测值，V_k为选定进行样本训练的观测值； 

为t+1时刻观测值 

的概率密度函数；

P(O^(l)/λ)为在模型λ下、观测值O^(l)的条件概率；

t表示观测时刻，T_l是电气特征值序列在时间维度上的采样长度，N是马尔科夫链状态的数目，M是每个状态对应可能的观测值的个数，L为离散观测样本的维数。

7.根据权利要求4所述的基于精确化模式判别的大电网暂态稳定分析方法，其特征在于：

所述的步骤S3b-2是采用离散隐马尔科夫模型比较步骤S3b-1所得离散观测值样本与步骤S3a-2所得的失稳的训练样本集和稳定的训练样本集。 

Images