短期负荷预测是在考虑气温因素和日期类型影响的前提下，对未来24小时、几天或几周的用电需求进行估算和预测^[1].

常见的负荷预测模型包括：极限学习机^[2]、支持向量机^[3]、卡尔曼滤波^[4]和随机森林^[5]等. 近年来，随着深度学习的发展，深度学习模型如卷积神经网络^[6]、循环神经网络（recurrent neural network，RNN）^[7]、深度置信网络^[8]等逐渐被应用到电力负荷预测中. 其中，RNN在语音识别、自然语言处理中应用广泛. 该网络通过循环的方式处理数据，可以在进行数据预测时更好地拟合数据在时间序列上的特性^[9]. 基于负荷序列的时序特性，RNN被许多学者引入负荷预测中^[10]. Tang等^[11]使用K-Means将数据集分类至同一个聚类中，再使用EEMD算法将负荷序列分解为较稳定的分量，针对每个分量建立双向循环神经网络（bidirectional recurrent neural network，BiRNN）模型进行预测. 然而RNN在计算距离较远时涉及雅可比矩阵的累乘，易出现梯度膨胀和梯度弥散^[12]. 长短时记忆网络（long short-term memory，LSTM）在网络的结构上相比RNN做出了改善，避免了RNN网络中的梯度膨胀和梯度弥散出现. Kong等^[13]在公开获取的真实住宅智能电表数据集中对LSTM以及多种同类模型进行测试，最终LSTM在居民家庭短期负荷预测方面优于其他同类模型.

门控循环单元（gated recurrent unit，GRU）是改良的LSTM，不容易使网络过拟合^[14-15]. 近年来，许多学者将GRU应用于负荷预测，并且在各自的实验数据集中都取得了良好的预测效果^[16-17]. Afrasiabi等^[18]使用卷积神经网络与GRU的组合模型预测风速和居民负荷，组合模型的预测效果优于对比模型. Tang等^[19]采用基于LSTM和GRU的BiRNN模型预测短期电力负荷，并在2个数据集上进行验证.

经验模态分解（empirical mode decomposition，EMD）可以将复杂的信号分解为相对平稳的模态分量，以减少不同分量之间的相互干扰^[20]. 集合经验模态分解（ensemble empirical mode decomposition，EEMD）^[21-23]在原始信号中引入正态分布的白噪声，显著改善EMD模式混合问题. Tang等^[24]提出EEMD和最小二乘支持向量回归（least squares support vector regression，LSSVR）的混合模型，并且在1993年3月至2010年1月的中国核能消费月度数据集中取得了良好的预测精度.

上述方法忽略了模型在迭代训练的过程中产生的误差信息. 本文提出包含误差修正环节的组合预测模型。模型使用叠式双向门控循环单元（stacked bidirectional gated recurrent unit, SBiGRU）学习负荷数据的主要特征，误差信息体现在SBiGRU模型预测产生的误差序列中；使用完整自适应噪声集成经验模态分解（complete ensemble empirical mode decomposition with adaptive noise, CEEMDAN）-SBiGRU模型对误差序列进行拟合；将SBiGRU模型与CEEMDAN-SBiGRU模型相结合得到本文组合预测模型.

GRU在设计过程中参考LSTM的网络结构，将LSTM中的遗忘门与输入门集合为更新门，同时融合细胞状态和更新状态. 与LSTM相比，GRU减少了一个门控单元，因此模型参数相对更少，运算效率显著提高. GRU的网络结构如图1所示. GRU通过t时刻的输入 $ {{{\boldsymbol{x}}}}_{{t}} $与t−1时刻的隐藏状态 $ {{{\boldsymbol{h}}}}_{{t−1}} $控制重置门 $ {{{\boldsymbol{r}}}}_{{t}} $与更新门 $ {{{\boldsymbol{z}}}}_{{t}} $， $\sigma $为激活函数Sigmoid； $ {{{\boldsymbol{h}}}}_{{t−1}} $经过重置门之后与 $ {{{\boldsymbol{x}}}}_{{t}} $相加并通过激活函数tanh将数据放缩至[−1，1]后得到 $ {{{\boldsymbol{h}}}}_{{t}} $.

(1) $ {{{\boldsymbol{z}}}}_{{t}}=\sigma \left({{{\boldsymbol{W}}}}_{{\rm{{z}}}}\left[{{{\boldsymbol{h}}}}_{{t-1}}{,}{{{\boldsymbol{x}}}}_{{t}}\right]\right){,} $

(2) $ {{{\boldsymbol{r}}}}_{{t}}=\sigma \left({{{\boldsymbol{W}}}}_{{{\rm{r}}}}\left[{{{\boldsymbol{h}}}}_{t {-1}}{,}{{{\boldsymbol{x}}}}_{{t}}\right]\right){,} $

(3) $ {\widetilde{{\boldsymbol{h}}}}_{{t}}={{\rm{tanh}}}\left({{\boldsymbol{W}}}\left[{{{\boldsymbol{r}}}}_{{t}}{{{\boldsymbol{h}}}}_{t {-1}}{,}{{{\boldsymbol{x}}}}_{{t}}\right]\right){,} $

(4) $ {{{\boldsymbol{h}}}}_{{t}}=\left({{\boldsymbol{1}}-}{{{\boldsymbol{z}}}}_{{t}}\right){{{\boldsymbol{h}}}}_{t{-1}}+{{{\boldsymbol{z}}}}_{{t}}{\widetilde {{\boldsymbol{h}}}}_{{t}}{.} $

式中：W_z、W_r、W分别为重置门、更新门、隐藏层网络状态的权值矩阵。

双向门控循环单元（bidirectional gated circulation unit，BiGRU）由前向GRU网络与后向GRU网络构成，结构如图2所示. BiGRU神经网络基于整个时间序列进行预测，它将隐藏层划分为前向与后向2个对立的部分，分别读取过去与未来时刻的信息. 第1层前向GRU计算当前时刻序列的顺序信息，第2层后向GRU计算相同时刻序列的逆序信息，最后BiGRU在t时刻的隐藏层状态由前向隐藏层状态和后向隐藏层状态加权得到. 其计算过程如下：

(5) $ {\overrightarrow{{\boldsymbol{h}}}}_{{t}}={{\rm{GRU}}}\left({{{\boldsymbol{x}}}}_{{t}}{,}{\overrightarrow{{\boldsymbol{h}}}}_{t\text{+1}}\right){,} $

(6) $ {\overleftarrow{{{\boldsymbol{h}}}}}_{{t}}={{\rm{GRU}}}\left({{{\boldsymbol{x}}}}_{{t}}{,}{\overleftarrow{{{\boldsymbol{h}}}}}_{t\text{-1}}\right){,} $

(7) $ {{{\boldsymbol{h}}}}_{{t}}={{{\boldsymbol{W}}}}_{{t}}{\overrightarrow{{\boldsymbol{h}}}}_{{t}}+{{{\boldsymbol{V}}}}_{{t}}{\overleftarrow{{{\boldsymbol{h}}}}}_{{t}}+{{{\boldsymbol{b}}}}_{{t}}{.} $

式中： ${\overrightarrow{{\boldsymbol{h}}}}_{{t}}$为t时刻前向隐藏层状态， ${\overleftarrow{{{\boldsymbol{h}}}}}_{{t}}$为t时刻后向隐藏层状态，W_t、V_t为 ${\overrightarrow{{\boldsymbol{h}}}}_{{t}}$与 ${\overleftarrow{{{\boldsymbol{h}}}}}_{{t}}$所对应的权重， ${{{\boldsymbol{b}}}}_{{t}}$为t时刻的隐藏层的偏置.

SBiGRU的网络结构是由多个BiGRU堆叠而成的。第n层BiGRU在t时刻的输入为 ${{x}}_{{t}}^{{n}}$，输出为 ${{x}}_{{t}}^{n\text{+1}}$，同时第n层的输出 ${{x}}_{{t}}^{n\text{+1}}$又作为n+1层的输入，这种多层网络结构设计增强了模型的数据泛化能力.

CEEMDAN在EEMD分解的过程中添加自适应的白噪声并且计算唯一的余量信号以获取本征模态函数（intrinsic mode function, IMF）分量，可以在较少的实验次数中重构信号序列，消除自适应白噪声在EEMD分解过程中所产生的误差.

定义 $ {{{\boldsymbol{E}}}}_{\text{i}}{()} $为通过EMD获得的第i个模态分量的计算算子，CEEMDAN算法的计算步骤如下. 1）在k次实验中，对原始信号 $ {{{\boldsymbol{x}}}}_{{t}}{+}\delta _{\text{0}}{\boldsymbol{\omega}} _{{j}} $进行分解，其中 $ \delta _{\text{0}} $为高斯白噪声的标准差， $ {\boldsymbol{\omega}} _{{j}} $为高斯白噪声. 通过EMD分解得到第1个IMF分量并获取唯一的余量信号 $ {{{\boldsymbol{r}}}}_{\text{1}}\left({t}\right) $：

(8) $ {{{\bf{IMF}}}}_{\text{1}}\left({t}\right)=\frac{\text{1}}{{k}}\sum\limits_{i\text{=1}}^{{k}}{{{\bf{IMF}}}}_{{1}}^{{i}}(t){,} $

(9) $ {{{\boldsymbol{r}}}}_{\text{1}}\left({t}\right)={{\boldsymbol{x}}}\left({t}\right)-{{{\bf{IMF}}}}_{\text{1}}\left({t}\right){.} $

2）继续获取第2个IMF分量：

(10) $ {{{\bf{IMF}}}}_{\text{2}}\left({t}\right)=\frac{\text{1}}{{k}}\sum\limits_{i\text{=1}}^{{k}}{{{\boldsymbol{E}}}}_{\text{1}}\left({{{\boldsymbol{r}}}}_{\text{1}}\left({t}\right)+\delta _{\text{0}}{{{\boldsymbol{E}}}}_{\text{1}}\left({\boldsymbol{\omega}} _{{j}}\left({t}\right)\right)\right){.} $

3）重复以上步骤，计算第n个余量信号：

(11) $ {{{\boldsymbol{r}}}}_{{n}}\left({t}\right)={{{\boldsymbol{r}}}}_{n\text{-1}}\left({t}\right)-{{{\bf{IMF}}}}_{{k}}\left({t}\right){.} $

则第n+1个IMF分量为

(12) $ {{{\bf{IMF}}}}_{n\text{+1}}\left({t}\right)=\frac{\text{1}}{{k}}\sum\limits_{i\text{=1}}^{{k}}{{{\boldsymbol{E}}}}_{\text{1}}\left({{{\boldsymbol{r}}}}_{{n}}\left({t}\right)+\delta _{{n}}{{{\boldsymbol{E}}}}_{{n}}\left({\boldsymbol{\omega}} _{{j}}\left({t}\right)\right)\right){.} $

4）重复步骤3），直到余量信号呈单调趋势，分解停止原始信号 $ {{\boldsymbol{x}}}\left({t}\right) $被分解为

(13) $ {{\boldsymbol{x}}}\left({t}\right)=\sum\limits_{i\text{=1}}^{{N}}{{{\bf{IMF}}}}_{{i}}{+{\boldsymbol{r}}}\left({t}\right){.} $

式中：N为最终模态分量的个数， $ {{\boldsymbol{r}}}\left({t}\right) $为最终单调的余量信号.

1）负荷预测. 使用训练集对第1阶段SBiGRU模型进行有监督训练，待模型收敛后预测训练集与测试集，得到训练集预测值 $ {{P}}_{\text{btrain}} $与测试集预测值 $ {{P}}_{\text{b}} $，SBiGRU训练过程中产生的误差序列 $ {{P}}_{\text{e}} $计算公式为

(14) $ {{P}}_{\text{e}}={{P}}_{\text{btrain}}-{{T}}_{\text{train}}{.} $

式中： $ {{T}}_{\text{train}} $为训练集的真实值.

2）误差预测. 使用CEEMDAN算法将 $ {{P}}_{\text{e}} $分解为数个IMF分量与1个残余分量，对每项分量分别建立SBiGRU模型进行学习与预测，累加预测结果得到误差的预测值 $ {{P}}_{\text{e}}^{'} $.

3）计算组合模型的预测值. 组合模型的预测结果 $ {P} $为测试集的预测值与误差的预测值之和，即：

(15) $ {P}={{P}}_{\text{b}}+{{P}}_{\text{e}}^{'}{.} $

组合模型的预测流程如图3所示.

用2016年加拿大阿尔伯塔省电力负荷数据验证组合模型的预测精度. 数据采样周期为1 h，每天包括24个采样点. 同时收集该区域的温度信息.如 图4所示为数据集中部分负荷序列L随时间t的变化曲线. 不同季节所考虑的气象因素应有所差别，本次实验只收集到了负荷数据和温度信息，为了避免这种差别对预测精度的影响，选取2016年7月至10月（夏季）共计2 928个样本作为研究对象.

由图4可以看出，负荷序列在双休日期间明显低于工作日，这说明日期类型会在一定程度上影响负荷值. 因此，在负荷预测中应考虑日期类型.

计算日期类型、气温与负荷序列之间的相关性. 日期类型用虚拟编码模式表示，以[0，1]表示[工作日，双休日].

计算得到气温与负荷序列的Pearson相关系数为0.768 2，两者具有较强的正相关性，因此须将气温作为负荷序的影响因素. 同时，计算得到日期类型与负荷序列的Pearson相关系数为0.182 73，两者呈弱相关性，结合图4的分析结果，双休日期间的负荷值明显低于工作日负荷值，因此将日期类型作为负荷序列的影响因素.

使用MaxMin的方法对数据进行归一化处理，将所有数据的值域归一化到[0，1]，将归一化的数据传入到模型中，再将模型的输出结果进行反归一化，计算的公式为：

(16) $ {{x}}_{{i}}^{'}=\frac{{{x}}_{{i}}-{{x}}_{\text{min}}}{{{x}}_{\text{max}}-{{x}}_{\text{min}}}{.} $

式中： $ {{x}}_{{i}}^{'} $、 $ {{x}}_{{i}} $分别为归一化后的值、数据值， $ {{x}}_{\text{max}} $、 $ {{x}}_{\text{min}} $分别为最大值、最小值.

使用2016年7月至9月共计92天的数据作为训练集，使用10月1~30日的数据作为测试集，以历史24个时刻点的数据预测未来24个时刻点的负荷数据.

将92天的训练集转化为三维矩阵，三维矩阵的维度信息为

(17) $ \left[{D}_{1}, {D}_{2},{ \cdots ,}{D}_{92}\right] , \left[{L}_{1}, {L}_{2},{ \cdots, }{L}_{24}\right] , \left[{X}_{1}, {X}_{2}, {X}_{3}\right]. $

式中 $：\left[{D}_{1},{D}_{2},{ \cdots },{D}_{92}\right]$为数据集长度92天， $[{L}_{1},{L}_{2},{ \cdots }, $ $ {L}_{24}]$为每天所包含的24个时刻的数据， $\left[{X}_{1},{X}_{2},{X}_{3}\right]$为当前时刻的负荷数据、气温、日期类型.

同理，测试集的维度信息为

(18) $ \left\{{D}_{1},{D}_{2}, { \cdots },{D}_{30}\right\}, \left\{{L}_{1}, {L}_{2},{ \cdots },{L}_{24}\right\} , \left\{{X}_{1},{X}_{2},{X}_{3}\right\}. $

模型的训练过程采用Adam优化算法. Adam结合了RMSprop善于处理非平稳目标与Adagrad善于处理稀疏梯度的优点，其计算公式为

(19) $ {{m}}_{{t}}=\mu {{m}}_{t\text{-1}}+\left(\text{1}{-\mu}\right){{g}}_{{t}}{,} $

(20) $ {{n}}_{{t}}={v}{{n}}_{t\text{-1}}+\left(\text{1}{-v}\right){{g}}_{{t}}^{\text{2}}{,} $

(21) $ {\widehat{m}}_{{t}}=\frac{{{m}}_{{t}}}{\text{1}-{\mu }^{{t}}}{,} $

(22) $ {\widehat{n}}_{{t}}=\frac{{{m}}_{{t}}}{\text{1-}{{v}}^{{t}}}{,} $

(23) $ \Delta \theta _{{t}}=-\frac{{\widehat{m}}_{{t}}}{\sqrt{{\widehat{n}}_{{t}}}+\epsilon }\eta {.} $

式中：g_t为梯度； $ \epsilon $为平滑指数，主要用于防止分母为零； $ \mu $、 $ {v} $为动量因子；m_t、n_t为梯度的一阶矩阵估计，是对期望 ${E}\left|{{g}}_{{t}}\right|$、 ${E}\left|{{g}}_{{t}}^{\text{2}}\right|$的估计； ${\widehat{m}}_{{t}}$、 ${\widehat{n}}_{{t}}$分别为m_t、n_t的修正，可对期望进行近似的无偏估计；Δθ_t为更新的目标参数；η为学习率.

模型采用平均绝对百分误差 $ {{E}}_{\text{MAPE}} $、均方根误差 $ {{E}}_{\text{RMSE}} $、最大绝对百分误差 $ {{E}}_{\text{Max}} $作为评价指标，计算公式分别为

(24) $ {{E}}_{\text{MAPE}}=\frac{\text{1}}{{n}}\sum\limits_{{i}{=1}}^{{n}}\left|\frac{{{x}}_{{i}}-{\widehat{{x}}}_{{i}}}{{{x}}_{{i}}}\right| \times \text{100} {\text{%}} {,} $

(25) $ {{E}}_{\text{RMSE}}=\sqrt{\frac{\text{1}}{{n}}\sum\limits_{i\text{=1}}^{{n}}{\left({{x}}_{{i}}-{\widehat{x}}_{{i}}\right)}^{\text{2}}}{,} $

(26) $ {{E}}_{\text{Max}}={{\rm{max}}}\left({{E}}_{\text{MAPE}}\right){.} $

式中： ${\widehat{x}}_{{i}}$为神经网络在i时刻的预测值.

实验在Ubuntu 18.04系统下利用Pytorch1.10完成，实验PC的处理器为Inter Core i5-8300H，内存为16 GB，显卡为NVIDIA GeForce GTX 1050Ti.

为了科学地设计SBiGRU模型结构，使用测试集对不同层数的SBiGRU模型进行实验分析. 以单层BiGRU，2+1层、3+1层、4+1层SBiGRU为实验样本，结合文献[5]，将迭代次数设置为较大经验值2 000，评价指标采用平均绝对百分误差与训练时长. 实验结果如表1所示. 表中，T为模型的训练时长.

根据表1确定第1阶段SBiGRU模型结构为2+1层，该模型由2层BiGRU与1层全连接层组合而成. 模型结构如图5所示. 以第一个BiGRU层的输入为例，{91，24，3}={输入数据集长度为91天，每日包含的24个时刻点，特征维度}.

模型训练方式：通过历史24个时刻点的负荷数据、日期类型和温度信息预测未来24个时刻点的负荷数据，以7月1日至9月29日的数据作为网络的输入，以7月2日至9月30日的负荷数据作为网络的输出，以Adam算法对网络参数进行优化. 为了获取训练过程中产生的误差样本，待SBiGRU模型收敛后，将训练集输入至训练完毕的模型中，得到如图6所示训练集的预测结果 $ {{P}}_{\rm{btrain}} $，根据式（14）求得如图7所示第1阶段SBiGRU模型在预测过程中产生的误差序列 $ {P}_{\text{e}} $.

观察图7可以发现，误差序列 $ {{P}}_{\text{e}} $周期性较差并且包含噪声与大量高频分量，这会增大误差序列的预测难度. 因此第2阶段使用CEEMDAN算法将 $ {{P}}_{\text{e}} $分解为8个IMF分量与1个残余分量，各项IMF分量的负荷L随时间t的变化曲线如图8所示.

观察图8可以发现，IMF₁~IMF₄分量的频率较高且周期性不明显，可以看作误差序列的高频分量. IMF₅~IMF₇分量的周期性明显，可以看作误差序列的周期分量. IMF₈分量的频率较低且周期性不明显，可以看作误差序列的低频分量. IMF₉为误差序列的残余分量.

将CEEMDAN算法分解得到的8个IMF分量与1个残余分量分别建立第2阶段的SBiGRU模型，其结构设计过程与第1阶段的SBiGRU结构类似，通过实验确定网络结构为2+1层，此时的网络的输入为误差分量，其特征维数为1，因此，第2阶段的SBiGRU输入层维度信息为{91，24，1}，此时设置各项SBiGRU模型的输入为7月1日至9月29的误差分量，输出为7月2日至9月30日的误差分量，同样以Adam算法对网络参数进行迭代优化.

将各个IMF分量 ${c}_{{i}}\left({t}\right)$以及残差分量 ${{r}}_{{n}}\left({t}\right)$的预测值累加得到误差的预测值，重构后的误差序列计算式为

(27) $ {p}\left({t}\right)=\sum\limits_{{i=1}}^{{n}}{c}_{{i}}\left({t}\right)+{{r}}_{{n}}\left({t}\right){.} $

实验通过10月1~29日的数据预测10月2~30日的负荷数据.

待第1阶段SBiGRU模型与第2阶段CEEMDAN-SBiGRU模型收敛后，通过第1阶段SBiGRU模型预测得到10月2~30日负荷序列的主要特征 $ {{P}}_{\text{b}} $，通过CEEMDAN-SBiGRU模型预测得到10月2~30日的误差预测值 $ {{P}}_{\text{e}}^{'} $，根据式（15）求得最终的预测结果P. 选取SBiGRU- SBiGRU、SBiGRU、BiRNN^[19]、传统负荷预测方法SVR做对比实验，其中SBiGRU- SBiGRU是以SBiGRU为误差修正模型的组合预测模型. 为了更清楚地观察预测效果，仅展示10月3~5日的预测结果，如图9所示. 由图可知，与SBiGRU- SBiGRU、SBiGRU、BiRNN和SVR相比，本文提出的组合模型具有更好的预测性能，能够更加有效的捕捉负荷序列的复杂非线性变化，因此可以更加准确与稳定的对负荷序列进行短期预测.

如图10所示为5种模型在预测集中的实际值A与预测值P的相关分布图，R为Pearson相关系数，用来评价T与P的相关性. 若图中相关分布点越趋近标准虚线y=x，表示模型预测值越接近真实值，相应的Pearson相关系数也越趋近1. 由图可知，模型对负荷序列的拟合程度最高，与其他模型相比具有最好的预测精度；SBiGRU- SBiGRU模型也取得了不错的预测精度，但整体上还是比本文模型的预测精度低；SBiGRU模型没有进行误差修正，因此预测精度比以上2种组合模型低，因为该模型采用多层的网络结构设计，所以预测精度高于BiRNN模型；机器学习模型SVR既无法深度挖掘负荷序列的特征，也不能提取负荷序列过去、未来时刻的时序信息，在实验结果中表现最差.

5种模型的预测结果误差对比如表2所示. 由表可知，SVR模型在预测集中表现最差，在测试集中平均准确率为96.90%；BiRNN模型能够学习负荷序列的时序信息，预测精度相对SVR模型得到了改善，为97.90%；SBiGRU模型在预测过程中考虑了过去、未来时刻的数据信息，并且以2+1层的网络结构提升模型的数据挖掘能力，因此预测精度高于BiRNN模型，测试集中平均准确率为98.21%；SBiGRU-SBiGRU模型在第1阶段中使用SBiGRU预测负荷的主要特征，在第2阶段中再次使用SBiGRU作为误差修正模型，提高了预测精度，日均准确率为98.45%；本文模型在第2阶段中使用CEEMDAN算法对误差序列进行分解，可以在不同的时间尺度上对误差序列进行预测，因此日均准确率在5种模型中最高，为98.86%.

如图11所示为5种模型的误差直方图，纵坐标为对应误差值的数量N，若P_e越集中于0附近，则代表模型的性能越好. 从图中可以看出，本文模型的误差直方图在0附近最为集中且最高，说明这种组合预测模型的性能在5种模型中表现最优.

如图12所示为5种模型的绝对误差箱线图. 由图可知，本文模型的箱线图的箱体范围最窄且上边缘最低，进一步说明该模型能够将预测误差控制在较小的范围内.

（1）本文模型、SBiGRU- SBiGRU组合模型的评价指标均优于SBiGRU、BiRNN、SVR单一模型，说明组合预测模型相比单一预测模型具有更好的预测精度.

（2）本文模型将误差序列分解为数个更为平稳的分量，并且对各分量建立SBiGRU预测模型，能够很好地拟合误差序列，在5种模型中具有最好的预测精度.

（3）在误差修正阶段，本文模型采用SBiGRU模型拟合每项误差分量. 下一步计划根据各误差分量的特点，采用不同的模型对误差分量进行预测，以提高误差分量的预测精度.