[发明专利]利用AP聚类-跳转持续改进MC的风电出力模拟生成方法

权利要求书

1.一种利用AP聚类-跳转持续改进MC的风电出力模拟生成方法，其特征是，它包括：对历史风电出力数据进行AP聚类、形成某一聚类类别下风电出力状态跳转序列、形成某一聚类类别下时间序列长度为一天的风电出力状态时间序列、在某一聚类类别下形成的初步风电出力时间序列上抽样并叠加风电出力波动量、根据类间转移概率矩阵模拟生成风电出力时间序列，具体内容为：

1)对历史风电出力数据进行AP聚类

近邻传播(Affinity Propagation，AP)聚类算法是一种基于“信息传递”的无监督聚类算法；首先建立样本集合的相似度矩阵S为式(1)，矩阵中的元素S_a,b表示样本x_a与x_b之间的相似度，采用负欧氏距离表示，为式(2)；S_a,a为聚类中心的参考度，作为能否成为聚类中心的评判标准，在初始时均设置为偏向参数，其直接影响最终聚类的个数，偏向系数的值越大，产生聚类个数越多；对单位时间段内既定采样间隔的实测历史风电出力，以其日均值和日最值差来构造能够反映对应各日风电出力特征的向量，并利用该向量对历史风电出力以天为单位进行AP聚类，其中：单位时间段一般是指月、季、年、多年；采样间隔通常是指1分钟、5分钟或10分钟；

S_a,b＝-||x_a－x_b||² (2)

式中：S为样本集合的相似度矩阵；S_a,b表示样本x_a与x_b之间的相似度；S_a,a为聚类中心的参考度；x_a和x_b为数据样本点；a＝1,2,……,n；b＝1,2,……,n；n为用于AP聚类的样本总数；||x_a-x_b||²表示x_a与x_b间的欧氏距离；

将样本集合的相似度矩阵S作为输入，然后计算每个样本点的吸引度r(x_a,x_k)与归属度a(x_a,x_k)；r(x_a,x_k)表示x_k点是否适合作为x_a点的聚类中心，即x_k对x_a的吸引度；a(x_a,x_k)表示x_a点是否选择x_k作为其聚类中心，即x_a对x_k的归属度；r(x_a,x_k)和a(x_a,x_k)越大，表明x_k越适合作为聚类中心；聚类过程中不断地更新吸引度矩阵和归属度矩阵，为避免迭代过程中发生震荡，引入阻尼系数，直到迭代产生稳定的聚类中心和类别归属情况；因此，近邻传播结束时，x_a的聚类中心确定为x_k；其中，x_k为式(3)：

x_k＝argmax{a(x_a,x_k)+r(x_a,x_k)} (3)

式中：argmax为求取函数值最大的自变量，即x_k为满足{a(x_a,x_k)+r(x_a,x_k)}取得最大值时的取值；a(x_a,x_k)为样本点x_a对样本点x_k的归属度；r(x_a,x_k)为样本点x_k对样本点x_a的吸引度；x_k和x_a为数据样本点；k＝1,2,……,n；a＝1,2,……,n；n为用于AP聚类的样本总数；

采用轮廓系数(Silhouette Coefficient，SC)作为聚类评价指标；轮廓系数越大表示聚类效果越好；轮廓系数SC(x_k)定义为式(4)-(6)：

式中：x_k，x_a，x_b为数据样本点；C(x_k)表示样本点x_k所属的聚类类别；C(x_a)表示样本点x_a所属的聚类类别；C(x_b)表示样本点x_b所属的聚类类别；轮廓系数SC(x_k)的值介于[-1,1]，反映了x_k作为聚类中心时，聚类中心与类内样本的平均距离I(x_k)是否显著区别于它到类外样本的平均距离O(x_k)；I(x_k)表示聚类中心与类内样本的平均距离；D(x_k,x_a)表示样本点x_k与x_a之间的距离；M表示与样本点x_k同属同一类别的数据点个数；O(x_k)表示聚类中心与类外样本的平均距离；D(x_k,x_b)表示样本点x_k与x_b之间的距离；G表示与样本点x_k不属于同一类别的数据点个数；

通过对历史风电出力数据以天为单位进行AP聚类，可以将出力水平与波动特性相似的日风电出力序列划归为同一类；

2)形成某一聚类类别下风电出力状态跳转序列

将历史风电出力数据按照(P′_min，P′_max)均匀划分成N个数据段，每个数据段对应一个状态，则每个状态所对应的风电出力范围大小为(P′_max-P′_min)/N；其中，P′_min为历史风电出力的最小值，P′_max为历史风电出力的最大值，N为状态数；针对传统马尔科夫链存在风电出力序列容易陷入某种状态而难以跳转，以致生成的风电出力时间序列可能出现在某一状态下持续时间过长的情况，通过采用状态跳转矩阵代替状态转移矩阵进行改进；

在由聚类重构序列产生状态转移矩阵时，当相邻两个数据之间的状态数跳变幅度大于所划分状态数的1/3时，认为这一转移不符合实际，将此处的状态转移矩阵值置零，进而得到改进的状态转移矩阵；

其次，将上述得到的状态转移矩阵对角线元素全部置零，并计算每个元素占该行所有元素之和的比例，将其作为新的概率值，所得到的概率矩阵即为状态跳转矩阵；状态跳转矩阵为式(7)-(8)；

p_ij＝P(E_t+1＝j|E_t＝i) (8)

式中：矩阵P为状态跳转矩阵；矩阵P的行对应风电出力当前的出力状态，列对应下一时刻的出力状态；N为状态数；i，j为风电出力状态；i＝1,2,……,N；j＝1,2,……,N；p_ij代表风电出力从状态i转移到状态j的概率；E_t、E_t+1为第t时刻和第t+1时刻风电出力值对应的状态；P(A|B)为条件概率函数，表示在B的条件下出现A的概率；

相应的累积状态跳转矩阵为式(9)：

Q中元素的取值为式(10)：

式中：矩阵Q为累积状态跳转矩阵；q_ul为累积状态跳转矩阵Q中第u行、第l列的元素；p_uj为状态跳转矩阵P中第u行、第j列的元素；N为状态数；u＝1,2,……,N；l＝1,2,……,N；j＝1,2,……,N；

通过上述累积状态跳转矩阵Q的形成过程，得到某一聚类类别下的累积状态跳转矩阵，即类内累积状态跳转矩阵；然后，根据类内累积状态跳转矩阵，抽样形成每一聚类类别下的状态跳转序列；在该聚类类别下，风电出力当前时刻的状态为E_t，下一时刻的状态为E_t+1，生成服从均匀分布的随机数ξ，若0ξ≤q_Et,1，则E_t+1＝1；若q_Et,lξ≤q_Et,l+1，则E_t+1＝l+1，其中，q_Et,1表示类内累积状态跳转矩阵Q的第E_t行、第1列的元素；q_Et,l表示类内累积状态跳转矩阵Q的第E_t行、第l列的元素；q_Et,l+1表示类内累积状态跳转矩阵Q的第E_t行、第l+1列的元素；

由于状态跳转矩阵的对角线元素为零，相应的类内累积状态跳转矩阵中对角线元素与其左边元素相等，即q_i,i-1等于q_ii，所以[q_i,i-1,q_ii]为空集，在利用马尔科夫链法生成风电出力状态序列时，ξ不会落在[q_i,i-1,q_ii]范围内；其中，q_ii为类内累积状态跳转矩阵Q中第i行、第i列的元素，i＝1,2,……,N；q_i,i-1为累积状态跳转矩阵Q中第i行、第i-1列的元素；

通过上述步骤，可以依据某一聚类类别下的类内累积状态跳转矩阵，抽样形成风电出力处在该聚类类别下时，相邻状态不相同的状态跳转序列；

3)形成某一聚类类别下时间序列长度为一天的风电出力状态时间序列

按照从风电出力进入到某个状态的第一个时刻开始记录，若其持续的时间为T，则记录该状态持续时间T一次的统计方法遍历历史风电出力时间序列，统计出风电出力在该状态下各个持续时间出现的次数，得到在该状态下持续时间的分布情况；

统计并绘制得到风电出力某一状态持续时间分布直方图后，采用一条平滑的曲线描述此直方图，能够避免对直方图抽样时存在的特殊性；在描述数据分布特性时，核密度估计(Kernel Density Estimation，KDE)不依赖于参数估计模型的选择，而且能够有效地避免直方图估计对于直方图组距和位置的依赖，所以采用KDE来描述风电出力状态持续时间分布特性；KDE表达式为式(11)：

式中：为风电出力状态持续时间的核密度估计函数；K(·)为核函数；h为带宽；x₁,x₂,……,x_W表示风电出力序列在某一状态下各个持续时间值的一组数据，数据数目为W个；x_w为风电出力状态持续时间的第w个样本值，w＝1,2,……,W；W表示风电出力序列在某一状态下出现了W个不同的持续时间；

在KDE中，带宽的选取决定了拟合曲线的平滑程度，带宽越大越平滑，但其拟合效果越差；选择高斯函数作为核函数，采用经验法则求取带宽，如式(12)：

h＝1.06σW^-1/5 (12)

式中：h为带宽；σ为风电出力状态持续时间正态分布标准差；W表示风电出力序列在某一状态下出现了W个不同的持续时间；

为得到某一聚类类别下时间序列长度为一天的风电出力的状态时间序列，在生成满足核密度估计函数的各状态持续时间的随机数集合后，遍历得到的该聚类类别下的风电出力状态跳转序列，并进行重复抽样：当前状态为E_t，则在状态E_t的持续时间随机数集合中任意抽取一个数值，作为状态E_t的持续时间，重复进行，直至生成的时间序列长度满足要求；

4)在某一聚类类别下形成的初步风电出力时间序列上抽样并叠加风电出力波动量

在形成某一聚类类别下的风电出力状态时间序列后，确定每一时刻的初步风电出力；风电出力当前时刻的状态为E_t，则可以随机生成该时刻的初步风电出力P′_t·cb的数值，P′_t·cb∈(P_Et.min，P_Et.max)，P_Et.min、P_Et.max分别为状态E_t对应的风电出力范围的最小值与最大值；

风电出力波动性是指相邻时刻的出力值存在差异，采用一阶差分波动量，即相邻2个单位时间内的出力差值来刻画波动特性；风电出力的一阶差分波动量表达式为式(13)：

ΔP′_t＝P′_t+1-P′_t (13)

式中：t表示风电出力时刻；ΔP′_t表示t时刻风电出力波动量；P′_t和P′_t+1分别表示t时刻和t+1时刻的风电出力值；

风电出力一阶差分波动量的概率密度分布拟合函数有：正态分布、t-location-scale分布、logistic分布；使用拟合精度更好的混合高斯分布刻画风电出力波动特性，其概率密度函数是若干个高斯概率密度函数的加权，当数据分布复杂时，混合高斯分布函数能够解决单高斯分布函数拟合精度的不足，其数学模型表达式为式(14)-(15)：

式中：f(x)为混合高斯分布函数；v＝1,2,……,V；V为混合高斯分布中单高斯分布的个数；a_v为混合系数；p_gs(x|b_v，c_v)为第v个单高斯分布的概率密度函数，其中，x表示服从单高斯分布的随机数；b_v和c_v分别代表第v个单高斯分布的均值和标准差；

采用最大期望(expectation-maximization，EM)算法求解混合高斯分布函数对应分布参数，EM算法分为两步：第1步称为期望(E)步，它基于初始值的参数或之前的迭代值来计算似然函数的期望；第2步是最大化(M)步，它将似然函数最大化，并转换为可获得的新参数值，EM算法是反复迭代的过程，直到上述两个步骤收敛；

在生成满足混合高斯分布函数的波动量随机数集合后，随机抽取波动量，并在形成的初步风电出力值上叠加波动量，如式(16)；遍历生成的初步风电出力时间序列，叠加波动分量完成，形成某一聚类类别下时间序列长度为一天的风电出力时间序列；

P′_t·fd＝P′_t·cb+ΔP′_t (16)

式中：P′_t·fd为某一聚类类别下风电出力时间序列t时刻的风电出力值；P′_t·cb为某一聚类类别下风电出力时间序列t时刻初步风电出力值；ΔP′_t表示t时刻风电出力波动量；

5)根据类间转移概率矩阵模拟生成风电出力时间序列

根据含风力发电的电力系统中长期优化调度的需要，确定模拟生成风电出力时间序列的时间长度，而中长期优化调度所需的时间长度一般为2天及以上；

分别统计相邻两天风电出力所属类别之间的转移次数，并计算相对应的转移概率得到类间转移概率矩阵，如式(17)-(18)，进而可以得到累积类间转移概率矩阵，如式(19)-(20)；

p_lj·rz＝P(H_e+1＝z|H_e＝r) (18)

式中：矩阵P_lj为类间转移概率矩阵；矩阵P_lj的行对应当日风电出力所属的聚类类别，列对应次日风电出力所属的聚类类别；L为AP聚类所产生的类别数量；r，z为风电出力所属的聚类类别；r＝1,2,……,L；z＝1,2,……,L；p_lj·rz代表风电出力所属的聚类类别由类别r转移到类别z的概率；H_e、H_e+1为当日和次日风电出力所属的聚类类别；P(A|B)为条件概率函数，表示在B的条件下出现A的概率；

式中：矩阵Q_lj为累积类间转移概率矩阵；q_lj·rg为累积类间转移概率矩阵Q_lj中第r行、第g列的元素；p_lj·rz为类间转移概率矩阵P_lj中第r行、第z列的元素；L为AP聚类所产生的类别数量；r＝1,2,……,L；z＝1,2,……,L；g＝1,2,……,L；

由当日风电出力所属的聚类类别，依据累积类间转移概率矩阵，抽样得到次日风电出力所属的AP聚类类别：当日风电出力所属的聚类类别为H_e，次日风电出力所属的聚类类别为H_e+1，生成服从均匀分布的随机数ε，若0ε≤q_He,1，则H_e+1＝1；若q_He,gε≤q_He,g+1，则H_e+1＝g+1，其中，q_He,1表示累积类间转移概率矩阵Q_lj的第H_e行、第1列的元素；q_He,g表示累积类间转移概率矩阵Q_lj的第H_e行、第g列的元素；q_He,g+1表示累积类间转移概率矩阵Q_lj的第H_e行、第g+1列的元素；

模拟生成既定采样间隔下，R天内的风电出力时间序列的步骤如下，R为正整数：

①判断第1天风电出力所属的AP聚类类别；

②通过执行步骤2)-4)，形成第1天的风电出力模拟生成时间序列，记为Ω₁，Ω₁＝{P′_1·fd(1),P′_2·fd(1),……,P′_t·fd(1)}，其中，Ω₁表示第1天的风电出力模拟生成序列；其中，P′_t·fd(1)表示第1天第t个时刻的模拟生成风电出力值，t表示风电出力时刻；

③由步骤5)计算得到风电出力累积类间转移概率矩阵，抽样确定第τ天风电出力所属的AP聚类类别，τ＝2,3,……,R；

④重复步骤2)-4)，得到第τ天的风电出力模拟生成时间序列，记为Ω_τ，Ω_τ＝{P′_1·fd(τ),P′_2·fd(τ),……,P′_t·fd(τ)}，其中，Ω_τ表示第τ天的风电出力模拟生成序列；P′_t·fd(τ)其中表示第τ天第t个时刻的模拟生成风电出力值；

⑤重复执行③-④，直至生成共R天的风电出力时间序列，记为Ω，Ω＝{Ω₁,……,Ω_τ,……,Ω_R}。