CN105205502A - 一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法 - Google Patents

Abstract

本发明公开了一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，包括：找出电压降落的时间点，在电压降落的时间点对应的扰动时刻，进行负荷动特性提取与分类；判断负荷类别间的变化是否具有马氏性，将所有数据按时间平均分段，对每段数据基于最大似然的思想建立马氏链的概率转移矩阵；判断数字特征是否改变，如果未改变，则转入步骤五，否则依据矩阵相应的数字特征对时间段的负荷数据进行聚类，对有改变了数据特征的各时间段负荷数据求概率转移矩阵；利用马尔可夫蒙特卡罗模拟，描述负荷变化情况；利用隐马尔科夫模型处理反映负荷类别转化的序列。本方法对马尔柯夫链蒙特卡罗模拟进行了改进，有效降低了矩阵迭代后陷入稳态的可能性。

Description

一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法

技术领域

本发明涉及一种基于马尔柯夫蒙特卡罗的负荷特性分类综合方法。

背景技术

负荷建模是电力系统建模中的一个基础和关键问题。建立能够反映负荷特性准确的负荷模型一直是挑战性的难点问题。负荷建模的最大困难在于负荷的随机性、时变性,它包括负荷大小的改变和负荷组成成分的改变，虽然如此，但是负荷特性存在一定的规律性。

为了在把握规律的基础上，解决负荷构成的随机性、时变性问题，对负荷动特性进行分类与综合。负荷动特性分类与综合是将不同时间采集的同一变电站的动态负荷特性数据中负荷构成成分相似的归为一类，并通过综合的方法建立每一类的负荷模型；且综合得到的模型能反映该类负荷的动特性特征，同时具有对不同强度扰动的内插外推能力。

发明内容

为解决现有技术存在的不足，本发明公开了一种针对电力负荷随机性、时变性的分类特点，提供了一种精确度高的针对电力负荷随机性、时变性的分类综合方法，克服了现有负荷建模方法无法考虑负荷的随机性、时变性问题。

为实现上述目的，本发明的具体方案如下：

一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，包括以下步骤：

步骤一：找出电压降落的时间点，在电压降落的时间点对应的扰动时刻，进行负荷动特性提取与分类；

步骤二：判断负荷类别间的变化是否具有马氏性，如果是，则进行步骤三，否则，结束；

步骤三：将所有数据按时间平均分段，对每段数据基于最大似然的思想建立马氏链的概率转移矩阵；

步骤四：判断数字特征即概率转移矩阵的稳态分布作为一个向量是否改变，当向量的中心距超过规定门阈值时，认为数字特征已经改变，依据概率转移矩阵相应的数字特征对时间段的负荷数据进行聚类，对有改变了数据特征的各时间段负荷数据求概率转移矩阵，然后转入步骤五，如果数字特征未改变，则直接转入步骤五；

步骤五：利用概率转移矩阵生成累积状态转移率矩阵，用马尔柯夫蒙特卡罗仿真(MCMC)对负荷状态在整个时段内的变化进行内插与外推，得到完整的负荷数据序列；

步骤六：利用隐马尔科夫模型处理步骤五的结果即反映负荷类别转化的序列，得到更为精确的负荷类别变化数据。

进一步的，所述步骤四中，数字特征取为各类别的稳态分布：

$J={\[L\left(\underset{m\→\mathrm{\∞}}{\lim }{P}^m\right)-{\mathrm{\π}}^{*}\]}^2,\underset{j=1}{\overset{6}{\Σ}}{p}_{ij}=1,p_{ij}>0$

式中，P为概率转移矩阵，p_ij为概率转移矩阵的第i行第j列的元素，函数的返回值为矩阵的任意行矢量，π^*表示第一步求取的最优稳态分布，求J的最小值。

对于这样一个非线性优化问题，应用遗传算法，容易得到近似极值，其中稳态分布的存在性可转化为最小化P^m的列矢量与列矢量均值的差的二范数。

进一步的，所述步骤四中，依据矩阵相应的数字特征对时间段的负荷数据进行聚类：设x为聚类向量，将目标函数设为样本到所属类别中心的距离，采用模糊C均值聚类(fuzzyc-meansalgorithm,FCM)，用迭代的方法求取最小目标函数。

FCM算法把样本空间x＝{x₁,x₂,…,x_n}分为c类(2≤c≤n)，任一样本点x_i不会严格划分为某一类。FCM用模糊划分，使用隶属度u_ij(0≤u_ij≤1)来确定样本点x_i属于第j(0≤j≤c)类的程度。如样本空间x的一个模糊子集所对应的隶属度矩阵是一个模糊隶属度矩阵，用U＝{U_ij}表示。隶属度矩阵U具有如下性质：

$\left\{\begin{array}{c}0\&le;u_{ij}\&le;1\\ \underset{j=1}{\overset{c}{\&Sigma;}}{u}_{ij}=1\\ 0<\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}<n\end{array}\right.---\left(1\right)$

式中，u_ij为隶属度。

FCM算法就是在式(1)的约束条件下使目标函数J最小化，即：

$J(U,c_1,...,c_c)=\underset{i=1}{\overset{c}{\&Sigma;}}\underset{j}{\overset{n}{\&Sigma;}}{u}_{ij}^m{d}_{ij}^2,---\left(2\right)$

式中，m∈[1,∞)是模糊加权系数；c_j是c类中第j类的聚类中心；是x_i到c_j的欧拉距离。

其步骤如下：

(1)给定聚类数c，模糊加权系数m，迭代停止阈值ε；设迭代次数k＝0，最大迭代次数k_max；用值在[0,1]间的随机数初始化隶属矩阵U，使其满足式(1)中的约束条件；

(2)由式(2)计算聚类中心

$c_j=\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}^m{x}_i/\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}$

式中，u_ij为隶属度，x_i为样本；

(3)由聚类中心c_j更新隶属度矩阵U^(k+1)，即

$u_{ij}^{(k+1)}={\&lsqb;\underset{i=1}{\overset{c}{\&Sigma;}}{(d_{ij}^{(k+1)}/d_{il}^{(k+1)})}^{2/(m-1)}\&rsqb;}^{-1},(l=1,2,...c)$

式中，d_ij为样本和各类中心的距离，上标(k+1)表示第k次迭代的相应数值；

(4)给定收敛的判别精度ε>0，若||U^(k+1)-U^(k)||≤ε，停止迭代；否则置k＝k+1，并返回步骤(2)；

(5)得到x的一个最优模糊C划分隶属度矩阵U＝{U_ij}和聚类中心c＝{c_ij}。

进一步的，所述步骤五中，在进行内插和外推时，用概率转移矩阵生成累积状态转移率矩阵：

P_cum中元素的取值如下：

$P_{cum,il}=\left\{\begin{array}{c}0,l=0\\ \underset{j<l}{\&Sigma;}{p}_{ij},0<l\&le;N+1\end{array}\right..$

其中，p_ij为概率转移矩阵P的第i行第j列的元素。换言之，经过本步骤的处理，累积状态转移率矩阵的P_cum,il为P的所有第i行前j列元素之和。

进一步的，所述步骤五中，利用马尔柯夫蒙特卡罗仿真MCMC法生成给定时间长度负荷类别变化时间序列(给定时间上限)步骤如下：

1)随机产生一个在区间[0,1]内的数u；

2)将u与Pcum的第p行元素(元素行号与当前所处状态号相同)进行比较，假设u落在累积状态转移率矩阵P_cum的P_cum,pq到P_cum,p(q+1)这个范围内，则认为负荷状态的序列的下一时刻状态为q；

3)如果生成的序列已经满足时间长度的要求，则转步骤1)继续，否则，令当前状态变为q，转步骤1)步继续，直到达到时间要求。

进一步的，所述步骤五中，上述步骤在实现时，每次模拟只对相邻故障时间间隔内的负荷状态进行模拟。

进一步的，所述步骤五中，因为故障时间点随机出现，各时间点之间的时间间隔不统一，相邻两个时间点之间可能还有其他的负荷状态，所以要统一步长为常数△t，累积状态转移率矩阵的处理相应也要进行改进：

转移概率矩阵中，P₁₂表示经过一个步长的时间后，从1状态转化为2状态的概率，也可以表示，现在为状态2，一个步长的时间之前的状态为1的概率，所以

式中，P_12,k为已知△t时刻，负荷状态为1，则第(k+1)△t时刻，负荷处于第2种状态的概率。

P_12,kt＝P_21,k为已知第(k+1)△t时刻，负荷状态为1，则△t时刻，负荷处于第2种状态的概率。

问题可以细化为假设已知△t、10△t两个时刻的负荷状态分别为a与b，已知马尔柯夫链P，采用蒙特卡洛仿真内插确定[△t,10△t]时间段内各时间点负荷状态。以确定5△t时刻的负荷情况为例，P⁴、(P⁵)^T分别为首、末两端与5△t时刻负荷状态的转换概率分布。

根据马尔可夫性，△t时刻的情况与5△t时刻有关，而5△t时刻的情况与△t时刻有关，即一个马尔柯夫过程的正向和反向都具有马尔柯夫性，所以，5△t时刻的情况应该结合△t、10△t两个时刻的负荷状态进行确定。

对P⁴进行处理，令

P_cum1中最后一列元素为1，元素的取值如下：

$P_{cum1,kl}=\left\{\begin{array}{c}0,l=0\\ \underset{j<l}{\&Sigma;}{p}_{kj,4},0<l\&le;N+1\end{array}\right.$

对(P⁵)^T进行处理，令

P_cum2中元素的取值如下：

$P_{cum2,kl}=\left\{\begin{array}{c}0,l=0\\ \underset{j<l}{\&Sigma;}{p}_{kj,5t},0<l\&le;N+1\end{array}\right.$

再将P_cum2各行归一化。

从△t、10△t两个时刻对5△t时刻的情况进行蒙特卡罗仿真，将P_cum＝(P_cum1+P_cum2)/2，然后利用前文所述的利用MCMC法生成给定时间长度负荷类别变化时间序列的处理方法进行仿真。

进一步的，所述步骤六中，隐马尔科夫模型HMM是一个双重随机过程，它有两个组成部分：

a马尔可夫链----描述类别的转移，用转移概率描述；

b一般随机过程----描述状态与观察序列间的关系，用观察值概率描述。

两个随机过程的关系如图1所示，显而易见，两者之间是嵌套的关系。

一个隐马尔科夫模型HMM常用五元数组来描述：λ＝(Q,O,A,B,π)。其中：状态的有限集合Q＝{q₁,q₂,…,q_n,}，n为马氏链的状态数，即对应本文中的各种负荷类别；观察值的有限集合O＝{o₁,o₂,…,o_m}，m为每种类别对应的可能的观察值数目；状态转移概率矩阵A；观察值概率矩阵B；初始概率分布向量π。

用Baum-Welch算法，基于HMM模型识别负荷特性的变化。步骤如下：

1)建立初始模型即待训练模型λ₀；

2)基于λ₀以及观察值序列O，训练新模型λ；

3)如果logP(O|λ)-logP(O|λ₀)<D，说明训练已经达到预期效果，算法结束，式中，P(O|λ)表示在实际模型为λ的情况下，观察序列为O的概率，D为设定门阈值；

4)否则，令λ₀＝λ，继续步骤2)工作。

训练过程前，需要定义两个辅助变量。第一个辅助变量为ξ_t(i,j)，表示给定模型λ和观察序列条件下，从i到j的转移概率：

$\begin{array}{c}{\mathrm{\&xi;}}_t\left(i,j\right)=P\left(q_t=i,q_{t+1}=j|X,\mathrm{\&lambda;}\right)\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right)a_{ij}{b}_j\left(O_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}{\underset{i=1}{\overset{N}{\&Sigma;}}\underset{j=1}{\overset{N}{\&Sigma;}}{\mathrm{\&alpha;}}_t\left(i\right)a_{ij}{b}_j\left(x_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}\end{array}$

式中，α_t(i)＝P(o₁,o₂,…,o_t,x_t＝i)，a_ij＝P(x_t+1＝j|x_t＝i)，b_j(x_t+1)＝P(o_t+1|x_t+1＝j)，β_t+1(j)＝P(o_t+1,o_t+2,…,o_T|q_t＝i,λ)。

第二个辅助变量为为t时刻处于状态S_i的概率，为整个过程中从状态S_i转出次数的预期，为从S_i跳转到S_j次数的预期。

训练过程：为当t＝1时，处于q_i的概率。

${\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_{ij}=\frac{\underset{t=1}{\overset{T-1}{\&Sigma;}}{\mathrm{\&xi;}}_t(i,j)}{\underset{t=1}{\overset{T-1}{\&Sigma;}}{\mathrm{\&gamma;}}_t\left(i\right)},1\&le;i\&le;N,1\&le;j\&le;N$

${\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{jk}=\frac{\underset{t=1,o_t=k}{\overset{T}{\&Sigma;}}{\mathrm{\&gamma;}}_t\left(j\right)}{\underset{t=1}{\overset{T}{\&Sigma;}}\mathrm{\&gamma;}\left(j\right)},1\&le;j\&le;N,1\&le;k\&le;N$

分别为状态转移概率矩阵A；观察值概率矩阵B中的元素的估计值。当收敛或达到预置次数时，即可得到状态转移概率矩阵A作为本发明的最终结果。

本发明的有益效果：

1.马氏链的概率转移矩阵更能体现负荷特性各类别之间的转化规律。

2.可精确确定负荷类别的故障时间点是随机出现的，各时间点之间的时间间隔不统一，相邻两个时间点之间可能还有其他的负荷类别，基于最大似然思想建立马氏链的概率转移矩阵，体现负荷类别转化时的真实情况，误差较小。

3.数字特征是对马氏链特性的体现，当各数字特征都明显改变的时候，反映的正是马氏链的变化。通过数字特征对时间段进行聚类，针对不同类的时段分别建立概率矩阵，有效减小了误差。

4.本方法对马尔柯夫链蒙特卡罗模拟进行了改进，有效降低了矩阵迭代后陷入稳态的可能性。

5.本发明之电力负荷时变性分类综合方法适用于针对同一变电站在一段时间内负荷变化情况的描述，能够满足误差要求，极具工程实用价值。

附图说明

图1两个随机过程的关系示意图。

图2本发明针对电力负荷时变性的分类综合方法流程方框图。

具体实施方式：

下面结合附图对本发明进行详细说明：

如图2所示，执行步骤01，开始；

其次，执行步骤02，计算统计量此统计量服从自由度为(n-1)²的χ²分布。

设m_ij表示X⁽⁰⁾(t)从状态i经过一步转移到状态j的频数,并将转移频数矩阵的各列之和除以各行各列的总和,所得到的值记为P'_j：

${P^{\&prime;}}_j=\frac{\underset{j=1}{\overset{m}{\&Sigma;}}{m}_{ij}}{\underset{i=1}{\overset{m}{\&Sigma;}}\underset{j=1}{\overset{m}{\&Sigma;}}{m}_{ij}}$

记

${P^{\&prime;}}_{ij}=\frac{m_{ij}}{\underset{j=1}{\overset{m}{\&Sigma;}}{m}_{ij}}$

则统计量

${\mathrm{\&chi;}}^2=2\underset{j=1}{\overset{m}{\&Sigma;}}\underset{k=1}{\overset{m}{\&Sigma;}}{N}_{jk}\left|ln({p^{\&prime;}}_{jk}/{p^{\&prime;}}_k)\right|$

给定置信度α，查表可得到值，如果则认为该序列符合马氏性，否则不是马氏链。

再次，执行步骤03，采取基于PSO的非线性规划求取概率转移矩阵。

本文采取基于粒子群模拟退火算法(particleswarmoptimizationsimulatedannealingalgorithm,PSOSA)的非线性规划求取概率转移矩阵。其步骤如下

1)约束条件为当i与k都是1到N之间的整数时，

$\left\{\begin{array}{c}\underset{k=1}{\overset{N}{\&Sigma;}}{P}_{ik}=1\\ P_{ik}\&GreaterEqual;0\end{array}\right.$

其中，P_ik为转移矩阵各元素。

2)确定适应度函数。适应度函数的输入量是马尔可夫矩阵，适应度函数的输出量是概率之和。某相邻两个故障时刻的负荷类别分别为甲乙，若按理论计算，即计及步长地按马氏链的规律转化，得到相应矩阵，将此矩阵中表明从甲到乙的概率计入计数器，以此类推，将各次转化的概率叠加，即得到适应度函数的输出量，问题也转化为求此适应度函数的最大值。

$\min f=min(-\underset{i=1}{\overset{m-1}{\&Sigma;}}{p}_{s_i{s}_{i+1}}^{\left(t_{i+1}-t_i\right)})$

其中，m为检测时间点数，t_i为第i个检测点，s_i为第i个检测点处的负荷类别，为中第s_i行第s_i+1列元素。

3)初始化粒子群,对于每一个粒子计算目标函数值,确定自身最优位置p_i和全体最优位置p_g。

4)执行一次SA搜索,在温度T下从初始位置x₀开始,对解产生随机扰动生成新解并用Metropolis准则判定是否接受,重复进行L_k(马尔科夫链长度)次直到得到一个新解x'.执行退火操作降低温度，T_k+1＝K*T_k。

5)执行一次PSO搜索，粒子群向个体最优与全局最优移动，并更新个体最优与全局最优。

$\begin{array}{c}{v}_i^d={\mathrm{wv}}_i^d+c_1{r}_1\left(p_i^d-x_i^d\right)+c_2{r}_2\left(p_g^d-x_i^d\right)\\ x_i^d=x_i^d+{\mathrm{\&alpha;v}}_i^d\end{array}$

其中，w为惯性因子，x为当前粒子所在位置，v为例子飞翔速度。

6)若f(x')<f(p_g)，随机在全体粒子中选取一个粒子i令x_i＝x',返回5)；若f(x')>f(p_g),则令x₀＝p_g，返回4)；若f(x')＝f(p_g)且模拟退火算法和粒子群算法都满足终止条件,结束计算。

最终得到x的最优值，即概率转移矩阵

然后，执行步骤04，求各类别的稳态分布

$J={\&lsqb;L\left(\underset{m\&RightArrow;\mathrm{\&infin;}}{\lim }{P}^m\right)-{\mathrm{\&pi;}}^{*}\&rsqb;}^2,\underset{j=1}{\overset{6}{\&Sigma;}}{p}_{ij}=1,p_{ij}>0$

式中，P为概率转移矩阵，p_ij为概率转移矩阵的第i行第j列的元素，函数的返回值为矩阵的任意行矢量，π^*表示第一步求取的最优稳态分布。求J的最小值。

对于这样一个非线性优化问题，应用遗传算法，容易得到近似极值，其中稳态分布的存在性可转化为最小化P^m的列矢量与列矢量均值的差的二范数。

依据矩阵相应的数字特征对时间段的负荷数据进行聚类：设x为聚类向量，将目标函数设为样本到所属类别中心的距离，采用模糊C均值聚类(fuzzyc-meansalgorithm,FCM)，用迭代的方法求取最小目标函数。

FCM算法把样本空间x＝{x₁,x₂,…,x_n}分为c类(2≤c≤n)，任一样本点x_i不会严格划分为某一类。FCM用模糊划分，使用隶属度u_ij(0≤u_ij≤1)来确定样本点x_i属于第j(0≤j≤c)类的程度。如样本空间x的一个模糊子集所对应的隶属度矩阵是一个模糊隶属度矩阵，用U＝{U_ij}表示。隶属度矩阵U具有如下性质：

$\left\{\begin{array}{c}0\&le;u_{ij}\&le;1\\ \underset{j=1}{\overset{c}{\&Sigma;}}{u}_{ij}=1\\ 0<\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}<n\end{array}\right.---\left(1\right)$

式中，u_ij为隶属度。

FCM算法就是在式(1)的约束条件下使目标函数J最小化，即：

$J(U,c_1,...,c_c)=\underset{i=1}{\overset{c}{\&Sigma;}}\underset{j}{\overset{n}{\&Sigma;}}{u}_{ij}^m{d}_{ij}^2,---\left(2\right)$

式中，m∈[1,∞)是模糊加权系数；c_j是c类中第j类的聚类中心；是x_i到c_j的欧拉距离。

其步骤如下：

(1)给定聚类数c，模糊加权系数m，迭代停止阈值ε；设迭代次数k＝0，最大迭代次数k_max；用值在[0,1]间的随机数初始化隶属矩阵U，使其满足式(1)中的约束条件；

(2)由式(2)计算聚类中心

$c_j=\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}^m{x}_i/\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}$

式中，u_ij为隶属度，x_i为样本。

(3)由聚类中心c_j更新隶属度矩阵U^(k+1)，即

$u_{ij}^{(k+1)}={\&lsqb;\underset{i=1}{\overset{c}{\&Sigma;}}{(d_{ij}^{(k+1)}/d_{il}^{(k+1)})}^{2/(m-1)}\&rsqb;}^{-1},(l=1,2,...c)$

式中，d_ij为样本和各类中心的距离，上标(k+1)表示第k次迭代的相应数值，。

(4)给定收敛的判别精度ε>0，若||U^(k+1)-U^(k)||≤ε，停止迭代；否则置k＝k+1，并返回步骤(2)；

(5)得到x的一个最优模糊C划分隶属度矩阵U＝{U_ij}和聚类中心c＝{c_ij}。

接着，执行步骤05，对系统进行内插和外推。

用概率转移矩阵生成累积状态转移率矩阵。

P_cum中元素的取值如下：

$P_{cum,kl}=\left\{\begin{array}{c}0,l=0\\ \underset{j<l}{\&Sigma;}{p}_{kj},0<l\&le;N+1\end{array}\right.$

利用MCMC法生成给定时间长度负荷类别变化时间序列(给定时间上限)步骤如下：

1)随机产生一个在区间[0,1]内的数u；

2)将u与Pcum的第p行元素(元素行号与当前所处状态号相同)进行比较，假设u落在累积状态转移率矩阵P_cum的P_cum,pq到P_cum,p(q+1)这个范围内，则认为负荷状态的序列的下一时刻状态为q；

3)如果生成的序列已经满足时间长度的要求，则转步骤1)继续，否则，令当前状态变为q，转步骤1)步继续，直到达到时间要求；

4)每次模拟只对相邻故障时间间隔内的负荷状态进行模拟。

但是，本文所解决的问题有其特殊性，因为故障时间点随机出现，各时间点之间的时间间隔不统一，相邻两个时间点之间可能还有其他的负荷状态，所以要统一步长为常数△t，累积状态转移率矩阵的处理相应也要进行改进处理。

转移概率矩阵中，P₁₂表示经过一个步长的时间后，从1状态转化为2状态的概率，也可以表示，现在为状态2，一个步长的时间之前的状态为1的概率，所以

式中，P_12,k为已知△t时刻，负荷状态为1，则第(k+1)△t时刻，负荷处于第2种状态的概率。：

P_12,kt＝P_21,k为已知第(k+1)△t时刻，负荷状态为1，则△t时刻，负荷处于第2种状态的概率。。

问题可以细化为假设已知△t、10△t两个时刻的负荷状态分别为a与b，已知马尔柯夫链P，采用蒙特卡洛仿真内插确定[△t,10△t]时间段内各时间点负荷状态。以确定5△t时刻的负荷情况为例，P⁴、(P⁵)^T分别为首、末两端与5△t时刻负荷状态的转换概率分布。

根据马尔可夫性，△t时刻的情况与5△t时刻有关，而5△t时刻的情况与△t时刻有关，即一个马尔柯夫过程的正向和反向都具有马尔柯夫性，所以，5△t时刻的情况应该结合△t、10△t两个时刻的负荷状态进行确定。

对P⁴进行处理，令

P_cum1中最后一列元素为1，元素的取值如下：

$P_{cum1,kl}=\left\{\begin{array}{c}0,l=0\\ \underset{j<l}{\&Sigma;}{p}_{kj,4},0<l\&le;N+1\end{array}\right.$

对(P⁵)^T进行处理，令

P_cum2中元素的取值如下：

$P_{cum2,kl}=\left\{\begin{array}{c}0,l=0\\ \underset{j<l}{\&Sigma;}{p}_{kj,5t},0<l\&le;N+1\end{array}\right.$

再将P_cum2各行归一化。

从△t、10△t两个时刻对5△t时刻的情况进行蒙特卡罗仿真，将P_cum＝(P_cum1+P_cum2)/2，然后利用前文所述的利用MCMC法生成给定时间长度负荷类别变化时间序列的处理方法进行仿真。

最后，执行步骤06，用Baum-Welch算法，基于HMM模型识别负荷特性的变化。步骤如下：

1)初始模型(待训练模型)λ₀，

2)基于λ₀以及观察值序列O，训练新模型λ；

3)如果logP(O|λ)-logP(O|λ₀)<D，说明训练已经达到预期效果，算法结束，式中，P(O|λ)表示在实际模型为λ的情况下，观察序列为O的概率，D为设定门阈值；

4)否则，令λ₀＝λ，继续第2步工作。

给定模型λ和观察序列条件下，从i到j的转移概率定义为ξ_t(i,j)。

$\begin{array}{c}{\mathrm{\&xi;}}_t\left(i,j\right)=P\left(q_t=i,q_{t+1}=j|X,\mathrm{\&lambda;}\right)\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right)a_{ij}{b}_j\left(o_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}{\underset{i=1}{\overset{N}{\&Sigma;}}\underset{j=1}{\overset{N}{\&Sigma;}}{\mathrm{\&alpha;}}_t\left(i\right)a_{ij}{b}_j\left(x_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}\end{array}$

式(6-1)中，α_t(i)＝P(o₁,o₂,…,o_t,x_t＝i)，a_ij＝P(x_t+1＝j|x_t＝i)，b_j(x_t+1)＝P(o_t+1|x_t+1＝j)，β_t+1(j)＝P(o_t+1,o_t+2,…,o_T|q_t＝i,λ)。

第二个辅助变量为为t时刻处于状态S_i的概率，为整个过程中从状态S_i转出次数的预期，为从S_i跳转到S_j次数的预期。

训练过程：为当t＝1时，处于q_i的概率。

${\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_{ij}=\frac{\underset{t=1}{\overset{T-1}{\&Sigma;}}{\mathrm{\&xi;}}_t(i,j)}{\underset{t=1}{\overset{T-1}{\&Sigma;}}{\mathrm{\&gamma;}}_t\left(i\right)},1\&le;i\&le;N,1\&le;j\&le;N$

${\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{jk}=\frac{\underset{t=1,o_t=k}{\overset{T}{\&Sigma;}}{\mathrm{\&gamma;}}_t\left(j\right)}{\underset{t=1}{\overset{T}{\&Sigma;}}\mathrm{\&gamma;}\left(j\right)},1\&le;j\&le;N,1\&le;k\&le;N$

与分别为状态转移概率矩阵A；观察值概率矩阵B中的元素的估计值。当收敛或达到预置次数时，即可得到状态转移概率矩阵A作为本发明的最终结果。

其中，步骤3、步骤4和步骤5中，建立概率转移矩阵并以此进行内插与外推均在MATLAB仿真平台中完成。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (10)

1.一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，包括以下步骤：

步骤一：找出电压降落的时间点，在电压降落的时间点对应的扰动时刻，进行负荷动特性提取与分类；

步骤二：判断负荷类别间的变化是否具有马氏性，如果是，则进行步骤三，否则，结束；

步骤三：将所有数据按时间平均分段，对每段数据基于最大似然的思想建立马氏链的概率转移矩阵；

步骤四：判断数字特征即概率转移矩阵的稳态分布作为一个向量是否改变，当向量的中心距超过规定门阈值时，认为数字特征已经改变，依据概率转移矩阵相应的数字特征对时间段的负荷数据进行聚类，对有改变了数据特征的各时间段负荷数据求概率转移矩阵，然后转入步骤五，如果数字特征未改变，则直接转入步骤五；

步骤五：利用概率转移矩阵生成累积状态转移率矩阵，用马尔柯夫蒙特卡罗仿真对负荷状态在整个时段内的变化进行内插与外推，得到完整的负荷数据序列；

步骤六：利用隐马尔科夫模型处理步骤五的结果即反映负荷类别转化的序列，得到更为精确的负荷类别变化数据。

2.如权利要求1所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，所述步骤四中，数字特征取为各类别的稳态分布：

$J={\&lsqb;L\left(\underset{m\&RightArrow;\mathrm{\&infin;}}{\lim }{P}^m\right)-{\mathrm{\&pi;}}^{*}\&rsqb;}^2,\underset{j=1}{\overset{6}{\&Sigma;}}{p}_{ij}=1,p_{ij}>0$

式中，P为概率转移矩阵，p_ij为概率转移矩阵的第i行第j列的元素，函数的返回值为矩阵的任意行矢量，π^*表示第一步求取的最优稳态分布，求J的最小值。

3.如权利要求1所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，所述步骤四中，依据矩阵相应的数字特征对时间段的负荷数据进行聚类：设x为聚类向量，将目标函数设为样本到所属类别中心的距离，采用模糊C均值聚类，用迭代的方法求取最小目标函数。

4.如权利要求3所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，用迭代的方法求取最小目标函数：

FCM算法把样本空间x＝{x₁,x₂,…,x_n}分为c类，2≤c≤n，任一样本点x_i不会严格划分为某一类，FCM用模糊划分，使用隶属度u_ij(0≤u_ij≤1)来确定样本点x_i属于第j(0≤j≤c)类的程度。如样本空间x的一个模糊子集所对应的隶属度矩阵是一个模糊隶属度矩阵，用U＝{U_ij}表示，隶属度矩阵U具有如下性质：

$\left\{\begin{array}{c}0\&le;u_{ij}\&le;1\\ \underset{j=1}{\overset{c}{\&Sigma;}}{u}_{ij}=1\\ 0<\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}<n\end{array}\right.---\left(1\right)$

式中，u_ij为隶属度；

FCM算法就是在式(1)的约束条件下使目标函数J最小化，即：

$J(U,c_1,...,c_c)=\underset{i=1}{\overset{c}{\&Sigma;}}\underset{j}{\overset{n}{\&Sigma;}}{u}_{ij}^m{d}_{ij}^2,---\left(2\right)$

式中，m∈[1,∞)是模糊加权系数；c_j是c类中第j类的聚类中心；是x_i到c_j的欧拉距离。

5.如权利要求4所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，FCM算法就是在式(1)的约束条件下使目标函数J最小化，其步骤如下：

(1)给定聚类数c，模糊加权系数m，迭代停止阈值ε；设迭代次数k＝0，最大迭代次数k_max；用值在[0,1]间的随机数初始化隶属矩阵U，使其满足式(1)中的约束条件；

(2)由式(2)计算聚类中心

$c_j=\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}^m{x}_i/\underset{i=1}{\overset{n}{\&Sigma;}}{u}_{ij}$

式中，u_ij为隶属度，x_i为样本；

(3)由聚类中心c_j更新隶属度矩阵U^(k+1)，即

$u_{ij}^{(k+1)}={\&lsqb;\underset{i=1}{\overset{c}{\&Sigma;}}{(d_{ij}^{(k+1)}/d_{il}^{(k+1)})}^{2/(m-1)}\&rsqb;}^{-1},(l=1,2,...c)$

式中，d_ij为样本和各类中心的距离，上标(k+1)表示第k次迭代的相应数值；

(4)给定收敛的判别精度ε>0，若||U^(k+1)-U^(k)||≤ε，停止迭代；否则置k＝k+1，并返回步骤(2)；

(5)得到x的一个最优模糊C划分隶属度矩阵U＝{U_ij}和聚类中心c＝{c_ij}。

6.如权利要求1所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，所述步骤五中，在进行内插和外推时，用概率转移矩阵生成累积状态转移率矩阵：

P_cum中元素的取值如下：

$P_{cum,il}=\left\{\begin{array}{c}0,l=0\\ \underset{j<l}{\&Sigma;}{p}_{ij},0<l\&le;N+1\end{array}\right.$

其中，p_ij为概率转移矩阵P的第i行第j列的元素，经过本步骤的处理，累积状态转移率矩阵的P_cum,il为P的所有第i行前j列元素之和。

7.如权利要求1所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，所述步骤五中，利用马尔柯夫蒙特卡罗仿真MCMC法生成给定时间长度负荷类别变化时间序列步骤如下：

1)随机产生一个在区间[0,1]内的数u；

2)将u与Pcum的第p行元素进行比较，假设u落在累积状态转移率矩阵P_cum的P_cum,pq到P_cum,p(q+1)这个范围内，则认为负荷状态的序列的下一时刻状态为q；

3)如果生成的序列已经满足时间长度的要求，则转步骤1)继续，否则，令当前状态变为q，转步骤1)继续，直到达到时间要求。

8.如权利要求1所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，所述步骤五中，因为故障时间点随机出现，各时间点之间的时间间隔不统一，相邻两个时间点之间可能还有其他的负荷状态，所以要统一步长为常数△t，累积状态转移率矩阵的处理相应也要进行改进：

转移概率矩阵中，P₁₂表示经过一个步长的时间后，从1状态转化为2状态的概率，也可以表示，现在为状态2，一个步长的时间之前的状态为1的概率，所以

式中，P_12,k为已知△t时刻，负荷状态为1，则第(k+1)△t时刻，负荷处于第2种状态的概率。

P_12,kt＝P_21,k为已知第(k+1)△t时刻，负荷状态为1，则△t时刻，负荷处于第2种状态的概率。

9.如权利要求1所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，所述步骤六中，隐马尔科夫模型HMM是一个双重随机过程，它有两个组成部分：

a马尔可夫链----描述类别的转移，用转移概率描述；

b一般随机过程----描述状态与观察序列间的关系，用观察值概率描述；

一个隐马尔科夫模型HMM常用五元数组来描述：λ＝(Q,O,A,B,π)，其中：状态的有限集合Q＝{q₁,q₂,…,q_n,}，n为马氏链的状态数；观察值的有限集合O＝{o₁,o₂,…,o_m}，m为每种类别对应的可能的观察值数目；状态转移概率矩阵A；观察值概率矩阵B；初始概率分布向量π。

10.如权利要求9所述的一种基于马尔柯夫蒙特卡罗的负荷特性综合分类方法，其特征是，用Baum-Welch算法，基于HMM模型识别负荷特性的变化，步骤如下：

1)建立初始模型即待训练模型λ₀；

2)基于λ₀以及观察值序列O，训练新模型λ；

3)如果logP(O|λ)-logP(O|λ₀)<D，说明训练已经达到预期效果，算法结束，式中，P(O|λ)表示在实际模型为λ的情况下，观察序列为O的概率，D为设定门阈值；

4)否则，令λ₀＝λ，继续步骤2)工作；

训练过程前，需要定义两个辅助变量，第一个辅助变量为ξ_t(i,j)，表示给定模型λ和观察序列条件下，从i到j的转移概率：

$\begin{array}{c}{\mathrm{\&xi;}}_t(i,j)=P(q_t=i,q_{t+1}=j|X,\mathrm{\&lambda;})\\ =\frac{{\mathrm{\&alpha;}}_t\left(i\right)a_{ij}{b}_j\left(O_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}{\underset{i=1}{\overset{N}{\&Sigma;}}\underset{j=1}{\overset{N}{\&Sigma;}}{\mathrm{\&alpha;}}_t\left(i\right)a_{ij}{b}_j\left(x_{t+1}\right){\mathrm{\&beta;}}_{t+1}\left(j\right)}\end{array}$

式中，α_t(i)＝P(o₁,o₂,…,o_t,x_t＝i)，a_ij＝P(x_t+1＝j|x_t＝i)，b_j(x_t+1)＝P(o_t+1|x_t+1＝j)，β_t+1(j)＝P(o_t+1,o_t+2,…,o_T|q_t＝i,λ)；

第二个辅助变量为为t时刻处于状态S_i的概率，为整个过程中从状态S_i转出次数的预期，为从S_i跳转到S_j次数的预期；

训练过程：为当t＝1时，处于q_i的概率；

${\stackrel{\&OverBar;}{\mathrm{\&alpha;}}}_{ij}=\frac{\underset{t=1}{\overset{T-1}{\&Sigma;}}{\mathrm{\&xi;}}_t(i,j)}{\underset{t=1}{\overset{T-1}{\&Sigma;}}{\mathrm{\&gamma;}}_t\left(i\right)},1\&le;i\&le;N,1\&le;j\&le;N$

${\stackrel{\&OverBar;}{\mathrm{\&beta;}}}_{jk}=\frac{\underset{t=1,o_t=k}{\overset{T}{\&Sigma;}}{\mathrm{\&gamma;}}_t\left(j\right)}{\underset{t=1}{\overset{T}{\&Sigma;}}\mathrm{\&gamma;}\left(j\right)},1\&le;j\&le;N,1\&le;k\&le;N$

分别为状态转移概率矩阵A、观察值概率矩阵B中的元素的估计值，当收敛或达到预置次数时，即可得到状态转移概率矩阵A作为最终结果。