0 引言

近年来，风电作为快速增长的发电产业已引起能源企业越来越多的关注.风能的不确定性及变化因素使得风场输出功率很难控制，因此电力系统中大规模风电的并网给电网的稳定性及供电的可靠性带来了严峻的挑战^[1-2].精确可靠的短期风电功率预测有利于风电机组和常规发电机组的调度安排及运营成本优化.

目前，风电功率预测的方法主要有物理方法、统计方法等^[3].统计方法旨在通过对过去历史数据的统计分析，找出输入输出数据之间的联系，在短期预测中应用非常广泛.传统的统计方法有ARIMA(autoregressive integrated moving average model)、卡尔曼滤波等^[4-5]，但这些基于线性回归的预测方法难以描述风电功率或风速时间序列的内在非线性特性.更先进的一类现代统计方法，如神经网络^[6]、SVM^[7]等，具有很强的非线性逼近能力，在短期和超短期预测中能获得很好的预测结果.

在SVM的基础上，MKL方法由文[8]提出，其思想是通过不同核函数的线性凸组合来降低核函数及其参数选择对建模或分类精度的影响，通过增强模型的自动选择能力来提高模型的泛化性，避免了SVM等单核学习方法对参数的依赖性^[9-10].文[11]将基于SimpleMKL算法的MKL方法应用于多分类的变压器故障诊断中，取得了很好的诊断效果.文[12]将MKL方法应用于若干混沌时间序列基准问题的预测中，取得了很好的预测效果.文[13]提出一种MKL-LSSVM方法，将其应用于永磁同步电机的混沌建模问题中，具有很好的建模效果.因此，在SVM成功应用于风电功率预测的基础上，将MKL方法应用于短期风电功率预测是一种可行的选择.

经验小波变换(EWT)^[14]是一种自适应小波分解方法，其核心思想是通过构造合适的正交小波滤波器组对原始信号的傅里叶谱进行自适应划分，以提取具有紧支撑傅里叶频谱的模态信号分量，与经验模态分解(empirical mode decomposition，EMD)^[15]方法相比，能避免产生模态混叠及虚假模态的现象.文[15]提出了一种基于集成EMD—近似熵和回声状态网络的风电功率组合预测模型，其预测精度优于最小二乘支持向量机(least squares support veotor maohine，LSSVM).文[16-17]将EWT应用于机械故障诊断，实验结果表明EWT方法明显优于EMD方法，诊断效果显著.文[18]将EMD与SVM相结合，应用于短期风电功率预测中，与单一预测方法相比，这种组合预测方法的预测效果更为显著.

考虑到EWT与MKL方法的各自优点，本文提出一类基于不同算法实现的EWT-MKL组合预测方法.先利用EWT构造正交小波滤波器组，以提取各个模态分量，针对各个模态分量形成的子序列构建不同的MKL预测模型，由SimpleMKL、MKL-wrapper、MKL-chunking三种不同的算法完成模型的训练过程，最终对预测结果进行叠加.为验证本文方法的有效性，将其应用于不同季节的加拿大亚伯达省风电场的短期风电功率单步或多步直接预测中，并以美国NREL实验室所提供的实测风速数据集为实例，应用于短期风电功率单步或多步间接预测中，在同等条件下，还与基于RBF核的SVM方法或WSVM方法进行对比.

1 EWT方法

EWT基于包含在所处理信号谱信息中探测到的傅里叶支撑，自适应构建出小波滤波器组.通过该小波滤波器组从信号中提取具有紧支撑傅里叶谱的不同模态.

待分析信号f(t)进行傅里叶变换，记为，假定傅里叶支撑区间[0，π]被分割为L个连续区段，ω_l位于每个分割区段的边界上，即每个区段为[ω_l-1，ω_l](l=1，…，L-1)，以ω_l为中心，定义过渡带变量τ_l，过渡带的宽度T_l=2τ_l.而经验小波定义为每一个分割区间上的带通滤波器，由Littlewood-Paley和Meyer小波的构造可获得，EWT算法的实现步骤为：

Step 1    寻找由傅里叶谱幅值的h个极大值点{P_i}(∀i∈[1，h])构成的集合，按降序排列(P₁≥P₂≥…≥P_h)并归一化在区间[0, 1]上，设定阈值为χ=P_h+δ(P₁-P_h)，其中δ为相对幅值比，所估计出的小波数目L需保证P_i≥χ即可.每个区段的边界点ω_l大于阈值χ两个连续极大值点的中心且有ω₀=0与ω_L=π.

Step 2    选择τ_l=γω_l(0 < γ < 1)，γ为参数且有：

(1)

使得{φ₁(t)，{ψ_l(t)}_l=1^L}构成一个L²(·)空间上的经验小波紧框架集合.其中，φ_l(t)是经验尺度函数，ψ_l(t)是经验小波函数.

Step 3    φ_l(t)、ψ_l(t)的傅氏变换、计算为

(2)

(3)

其中函数β(x)的构造为

(4)

它满足性质：

Step 4    同经典小波变换，由式(3)、式(4)可知，细节系数W_f^e(l，t)及近似系数W_f^e(0，t)的计算为

(5)

(6)

其中，F^-1(·)表示傅里叶逆变换，表示复共轭.

由EWT的算法实现可知，信号的重构由式(7)给出：

(7)

信号分量f_n的表达为

(8)

(9)

其中，*为卷积符号.

2 MKL方法

SVM方法在复杂的非线性系统建模、控制中取得了成功应用，与其它核学习方法一样，其建模效果依赖于核函数及其参数的正确选择，往往需要交叉验证的重复步骤以获得合适的参数.与SVM方法相比，MKL学习方法呈现出一定的灵活性，可以使得核参数在较宽范围内进行选择，尤其是当数据具有不同来源时，可以根据其特点恰当地选择不同类型的核函数.

本文所使用的MKL方法中，其核函数k(x_i，x_j)考虑的是不同基础核函数的凸线性组合形式，即：

(10)

其中，d_m为权值系数，d_m≥0，，M为核函数的总数目；基础核函数k_m(x_i，x_j)=〈Φ(x_i)，Φ_m(x_j)〉，Φ_m(x)为x在特征空间上的映射.若给定训练数据(x_i，y_i)_i=1^N，其中x_i∈R^D，在原空间中可构建优化问题^[19]：

(11)

其中，ε为线性不敏感损失变量，ξ、ξ^*为松弛变量，C为惩罚因子，w_m是Φ_m(·)的权系数向量且w_m=d_mw′_m.

若f(x)为实值函数，ε不敏感损失函数定义为L^ε(x，y，f)=|y-f(x)|_ε=max(0，|y-f(x)|_ε-ε)，式(11)的优化问题提供了在分块内部具有非稀疏解、在m个分块水平上寻找稀疏解的可能性.

在文[19]的基础上，令变量，引入拉格朗日乘子α、，其中，α=[α₁，…，α_N]^T，，并经过适当的变形，可得形如式(12)的对偶优化问题：

(12)

其中，

(13)

式(12)的最优解在拉格朗日泛函的鞍点上取得，即有：

(14)

令式(14)中的L(·)对r的偏导数为0，得到.因此，式(14)可简化为

对目标函数L的优化问题求解等价于L关于α和求最小，而关于核函数的权值向量d求最大，即有：

(15)

定义变量，若α^*、是式(15)的最优解，则为最小且对于所有的α、在上述约束条件下均成立，因此式(15)的鞍点问题求解可等价于半无限线性规划问题求解，即：

(16)

由于α、的可能取值，式(16)有无穷多的约束.为获取式(16)的近似解，可使用“违背正则化最大约束”的收敛准则，即算法的停止条件是：

(17)

其中，η_MKL为MKL的精度参数，d^t和θ^t为(t-1)次迭代时的最优解，α^t和相应于下一次迭代时违背其约束条件新得到的最大值.

求解式(16)可使用基于列生成技术的MKL-wrapper算法及基于块分解技术的MKL-chunking算法.

2.1 MKL-wrapper算法

MKL-wrapper算法将优化问题分解为一个内部与外部子问题，通过应用内部问题的解作为输入，求解外部问题而交替进行，反之亦然，迭代至算法收敛，其外部循环构成一个“限制主问题”，即对于固定的α、值，基于已有的线性规划(LP)求解器获取最优的d值，而在内部循环中使用经典的SVM求解器识别不满足的约束条件，这在实际处理上尤为简单.

MKL-wrapper算法的实现步骤为：

Step 1    参数初始化.设置ε_MKL，t=1，S⁰=1，θ¹=-∞，d_m¹=1/M，m=1，2，…，M.

Step 2    计算，.

Step 3    计算.

Step 4    由式(16)，求解：

Step 5    由式(17)的准则可知，若|1-S^t/θ^t|≤η_MKL，则算法停止；否则，t=t+1，返回Step 2.

2.2 MKL-chunking算法

应用SVM优化问题的特殊结构，使用标准的优化软件包进行训练时，数据集包含的样本一般不超过数千个，而MKL-chunking算法使用了一种chunking分解技术，克服了此种限制且能够同时优化α、及d.分解的关键技术在于：少部分的优化变量处于工作集，其余的变量处于“冻结”状态.因此，MKLchunking算法使得求解系列具有“不变尺度”的SVM对偶优化子问题成为可行.

算法实现中，为了选择优化的下一组变量，需要临时计算SVM^light优化软件的输出，即：

(18)

若核函数的权值d变化，则储存的值无效，需要重新计算.为了避免全部重新计算，定义一个M×N的矩阵，其元素将保存每个核函数的输出，即：

(19)

因此，当d变化时，能被有效地重新计算，即有：

(20)

假设优化子问题的规模为Q，给定ε_MKL值.注意到各次迭代过程中，LP均会额外增加一个约束项，为算法描述的简洁，省略了“非活动”约束项的移除步骤. MKL-chunking算法的具体实现步骤为：

Step 1    设置η_MKL，t=1，g_m，i=0，，，d_m¹=1/M，m=1，…，M，i=1，…，N.

Step 2    由和α、选择Q个子优化变量，即i₁，…，i_Q，置中间变量α^old=α，.

Step 3    关于所选择的变量，求解SVM的对偶优化问题，更新α、，计算：

Step 4    计算，其中，

Step 5    若|1-S^t/θ^t|≥η_MKL，则：

否则：

Step 6    计算，i=1，2，…，N，由式(17)的准则可知，若|1-S^t/θ^t|≤η_MKL，则算法停止；否则，t=t+1，返回Step 2.

2.3 SimpleMKL算法

考虑到3.1节中基于半无限线性规划的MKL方法，其迭代求解算法往往需要多次迭代才能收敛于合理的解，文[20]给出了一种通过加权L₂范数正则化公式且采用具有附加的权值约束使得稀疏的核组合成为可行的新算法，除了对核组合的权值进行训练之外，还需求解标准的SVM优化问题，称之为SimpleMKL.

该算法实现过程中，线性核组合的稀疏性由核函数权值的L₁范数约束所控制，目标函数的求解可看作为一个光滑的凸优化问题，其算法简单，可基于SVM目标值应用梯度下降算法求解.在式(11)的基础上，SimpleMKL中的约束优化问题为

(21)

J(d)在对偶空间的形式为

(22)

其中，α、为拉格朗日乘子，根据强对偶定理可知SimpleMKL中的约束优化问题等价于其对偶问题，即：

(23)

考虑到J(d)的可微性，有：

(24)

其中，α^*、是α、的最优值.

算法执行时，一旦计算出J(d)的梯度，核权值系数d将通过梯度下降算法更新并确保d的约束满足.算法的停止准则为对偶间隙或最大的迭代次数，或者是2次连续迭代步数之间d的变化小于某设定值.更详细的算法描述参见文[20].

3 预测实例

本节的实例中，按照时间序列建模的方式，将EWT-MKL方法应用于不同地区的短期风电功率预测实例中，所建立的预测模型为

(25)

其中，Δ表示预测步长；x_t是历史时间序列值构成的多维输入向量，即x_t=[y_t-1，…，y_t-D]，D为嵌入维数；N为训练集大小；b^*为偏置量.

为评价本文所提出的组合预测方法的有效性，还与ε-SVM、SimpleMKL、MKL-wrapper等单一预测方法进行比较.实验均在Linux环境下，操作系统为Ubuntu12.04，基于Intel(R)双核处理器，主频为2.13 GHz，4 G内存的计算机上运行且SimpleMKL用Matlab语言，MKL-wrapper、MKL-chunking用python语言完成.

本节实验不同实例中，MKL方法均选择3个多项式核函数、6个RBF核函数作为基函数，形式为

(26)

(27)

其中，分别选取多项式的阶次p=1，2，3，选择核函数宽度2δ²∈{0.01，0.1，1，10，50，100}.设置ε=0.01，η_MKL=1e-4，正则化参数C=200.

评价指标分别采用平均绝对误差E_MAE、均方根误差E_RMSE、MAE的标准差σ_Std^[22]，其中σ_Std表示为

(28)

另一方面，为全面评价算法的应用效果，实例1为风电功率直接预测的方式，实例2为通过风速预测，再经“风速—功率”特性曲线转换间接预测风电功率.

“风速—功率”特性曲线表示风速与功率的对应关系，其定量关系为^[21]

(29)

其中，为预测风速，v_cut-in为“切入”风速，v_norm为“额定”风速，v_cut-out为“切出”风速，C(·)是由“风速—功率”特性曲线获得的曲线函数，p_norm为额定功率.

实际应用中，可结合所测量的历史风速与风电功率数据，获得与式(29)等价的“风速—功率”特性曲线，其实现步骤为：

Step 1    将风速区间[v_cut-in，v_norm]按0.2 m/s等间隔分割，计算第i个区间相应的风电功率均值μ_P(i)及标准差σ_P(i)，[μ_P(i)-σ_P(i)，μ_P(i)+σ_P(i)]外的样本点被舍弃.

Step 2    抛弃散点后，第i个风速区间所对应的风电功率在整个范围内被均分为10个功率区间，计算功率出现的频数F(i，j)(第j个功率区间的样本数)，并将其归一化，即，并按降序排列.

Step 3    设定阈值ρ_T，选择符合条件的前k个区间，则第i个风速区间内的风电功率为

(30)

其中p(i，j)为第i个风速区间中第j个功率子区间内的平均功率值.

因此，最终的“风速—风功”曲线由整个风速区间上的风电功率通过“滑窗”法依次获得.显然，短期风电功率间接预测的误差包括通过风速预测误差及“风速—功率”转换所形成的误差两个部分.

3.1 风电功率直接预测实例

采用加拿大亚伯达省电网公司(the Alberta electric system operator，AESO)提供的某风电场在2011年全年的实际风电功率数据集^[22]，该风电数据由SCADA系统采集并传送给AESO，其采样间隔为10 min.实验中将原始数据的连续3个采样点进行平均，作为一个新采样点，即采样间隔为30 min.选取4月24日至5月15日之间的全部数据作为训练，而以5月16日的数据作为测试数据集，应用EWT-MKL方法分别构建提前30 min、60 min、90 min的单步或多步预测模型来验证本文方法的有效性.

首先用EWT方法对训练数据集的原始风电功率时间序列进行变换，选择相对幅值比δ=0.5，可计算出小波的数目L=5.因此，EWT的输出由1个尺度函数与5个小波函数组成的带通滤波器组所构成. 图 1给出了经EWT变换后所提取的6个模态分量信号.

其次，针对风电功率的各分量信号时间序列(F₀~F₅)，可分别构建提前30 min、60 min、90 min的MKL预测模型，其中，F₀分量预测模型的嵌入维为8维，其它分量预测模型的嵌入维为10维.

图 2给出了针对F₀序列构建直接单步预测模型时，应用SimpleMKL、MKL-wrapper、MKL-chunking三种不同算法的MKL方法得到的核函数权值的分布对比图，其中横坐标所表示的d的索引号中，前3个对应于多项式核，余下的对应于高斯核.从图 2可以看出，无论采用何种算法实现的MKL模型，其核权值皆具有相似的稀疏解.

最后，根据各子序列预测模型的预测输出值，对预测结果进行叠加即可得到最终输出. 图 3给出了在测试集上基于不同算法实现的EWT-MKL方法进行提前30 min预测的结果比较. 图 4给出了其预测绝对值误差结果对比，由图 3、图 4可以看出，不同的EWT-MKL方法均能较为准确地预测实际风电功率值，预测误差波动都相对较小，呈现出较好的预测效果.

为进一步衡量EWT-MKL方法的预测效果，表 1~表 3列出了对测试数据集的数据进行提前30 min、60 min和90 min的预测结果数值对比及3种MKL算法收敛时间对比，SVM的时间是通过5折交叉验证寻优所用的时间.

由表 1~表 3可以看出，与SVM方法相比，MKL方法的预测精度略有提升，但提升的效果并不显著. SVM方法的RBF核参数及正则化参数是通过网格寻优的5折交叉验证方法所获得，而MKL方法避免了参数寻优的繁琐选择阶段，EWT-MKL组合预测方法则显示出了较高预测精度，与单一预测方法相比，提升效果较为明显.此外，通过算法训练的收敛时间对比可知，MKL-chunking算法的训练时间最快且优于SVM的训练时间，其次为SimpleMKL算法.

最后，随机选择2011年的春、夏、秋、冬四个季节中的4个测试周以E_MAE为评价指标进行验证，即选择2月的第2周、5月的第3周、8月的第2周、11月的第1周，进一步验证本文组合预测方法的有效性，不同预测模型对测试周的各天分别进行提前30 min、60 min、90 min预测，取7天的预测指标E_MAE平均值作为该测试周的预测指标E_MAE，由4个测试周的E_MAE的平均值得到全年的E_MAE，具体数值如表 4~表 6所示.从表 4~表 6的结果可以看出，无论是提前30 min还是提前60 min、90 min的短期风电功率多步预测，EWT-MKL方法均呈现出较高的预测精度，从全年预测结果的平均来看，EWT-MKL-chunking方法的预测效果相对较好.

从本节的实验可知，风电功率预测值和实际值都会存在一定程度的偏差.一般情况下，通过对训练数据进行交叉验证的方式，可调整模型的参数，若误差指标偏离较大，则需从重新考虑训练模型的嵌入维数及组合核函数的个数及类型，通过交叉验证实验的较优结果进行选取，从而使得预测模型的指标满足需求.然后，可将训练好的模型用于实际预测中，以起到制定合理的日发电计划调度，降低旋转备用容量及优化运营损失的作用.

3.2 风电功率间接预测实例

本小节的实例实验采用美国国家可再生能源实验室(the national renewable energy laboratory，NREL)所提供的西部风电数据集^[23]，从中选取丹佛市以西10英里位置，68个网格点，即680个Vestas v90 3-M型号的风机，所提供的时间间隔为10 min的平均风速与风电功率值，将数据集的连续2个数据采样点进行平均，作为一个新的数据，即数据采样的时间间隔为20 min.测试数据集包含的风速数据从2004年12月1日至2004年12月7日，之前的14 d所包含的风速数据作为训练数据集.

与文[21]的处理一致，通过分析历史风速、风电功率数据，参数设置为v_cut-in=2.64 m/s，v_cut-out=33.17 m/s，v_norm=20.44 m/s，按照等价的“风速—风功”特性曲线实现步骤，最终的“风速—功率”特性曲线由整个风速区间上的风电功率通过“滑窗”法依次获得，如图 5所示. 图 5中的黑色实线部分表示由实测风速、风电功率数据拟合出的风速—风电功率特性.

应用EWT-MKL组合预测方法分别构建提前60 min、120 min、180 min的风速预测模型，结合图 5的“风速—功率”特性曲线，以实现对风电功率的间接预测.首先用EWT方法对训练数据集的原始风电功率时间序列进行变换，选择相对幅值比δ=0.5，可计算出小波的数目L=4.因此，EWT的输出由1个尺度函数与4个小波函数组成的带通滤波器组所构成. 图 6给出了经EWT变换后所提取的5个模态分量信号.

其次，针对风速的各模态信号分量，以时间序列建模方式可分别构建提前60 min、120 min、180 min的MKL预测模型，各子序列模型的嵌入维均为6维.最后，对各分量预测模型的输出进行叠加即可得到最终的风速预测值，由图 5及式(29)计算得到对应风电功率预测值. 图 7、图 8给出了在测试集上(12月1日至7日)基于不同算法实现的EWT-MKL方法进行提前60 min间接预测时，转换后的风电功率多步预测结果对比及各点预测误差对比，图 9给出了相应的预测误差分布直方图(误差密度分布图)比较.由图 7~图 9的结果可以明显看出，不同的EWT-MKL方法均能较为准确地预测实际风电功率值，预测误差波动都相对较小，呈现出较好的预测效果.

为进一步定量比较不同的EWT-MKL方法的预测效果，在测试数据集上进行提前60 min、120 min和180 min的风速及转换后的风电功率多步预测误差的数值比较分别由表 7~表 9列出.由表 7~表 9可以看出，与SVM方法相比，MKL方法的预测精度略有提升，而EWT-MKL组合预测方法表现出了较高预测性能.此外，本文的EWT-MKL组合预测方法的预测精度也明显优于文[21]的结果.文[21]应用了基于WSVM方法的风电功率间接预测方法.

为了更直观地表明不同方法在不同预测步长下的预测性能，在测试集上分别进行提前20 min~180 min的单步与多步间接预测，其功率预测值与真实值的误差E_MAE随预测步长变化的情况由图 10给出.由图 10的结果比较可以看出，不同EWT-MKL方法均具有较小的预测误差且EWT-MKL-chunking方法的结果相对最好.

4 结论

针对短期风电功率预测，提出基于EWT-MKL的组合预测模型，基于SimpleMKL、MKL-wrapper、MKL-chunking三种算法实现MKL模型的训练；首先将其应用于不同季节的加拿大亚伯达省风电场的短期风电功率单步或多步直接预测中；其次，以美国NREL实验室所提供的实测风速数据集为实例，应用于短期风电功率单步或多步间接预测中.在同等条件下，还与SVM方法或WSVM方法进行了比较，结论如下：

1) SVM方法往往由于参数的不同初始化设置，使得预测性能变差，MKL方法则避免了SVM方法选择适当的核函数及相应“超参数”的困难，提升了模型的泛化能力.

2) EWT方法是一种自适应小波分解方法，鉴于其能避免EMD方法分解时产生的模态混叠及虚假模态缺点，将其与MKL方法结合并应用于短期风电功率预测中是一种可行且有效地选择.

3) 结合EWT和MKL二者优点的EWT-MKL组合预测方法，既保持了MKL方法较少依赖于模型参数及核函数的优点，又进一步提高了短期风电功率的预测能力.