基于迁移模糊系统的短期电力负荷预测建模

Short-Term Power Load Forecasting Modeling Based on Transfer Fuzzy System

摘要: 针对电力负荷数据缺失导致预测精度降低的问题，本文提出基于TSK迁移模糊系统(TSK-TFS)结合变分模态分解(VMD)、迁移成分分析(TCA)和改进斑马优化算法(IZOA)的短期电力负荷预测模型(IZOA-VMD-TSK-TFS-TCA)。首先利用VMD将电力负荷数据分解为若干子序列，并利用TCA将与电力负荷相关的因素降维；其次对斑马优化算法进行改进，利用改进后的斑马优化算法(IZOA)对TSK-TFS的参数寻优，并利用减法聚类算法得到聚类个数，把源域中的数据输入TSK模糊系统训练得到前件参数和后件参数并保留，继承参数并利用一部分目标域数据训练得到后件参数；最后根据得到的后件参数并经过计算得到测试集(另一部分目标域数据)若干子序列的预测值，将各个子序列的预测值叠加得到短期电力负荷的预测值。仿真实验结果表明，本文提出的IZOA-VMD-TSK-TFS-TCA短期电力负荷预测模型具有较高的预测精度，经过统计检验也证实了该模型具有较优的预测性能。

Abstract: To address the issue of reduced prediction accuracy caused by missing power load data, this paper proposes a short-term power load prediction model (IZOA-VMD-TSK-TFS-TCA) based on TSK Transfer Fuzzy System (TSK-TFS) combined with Variational Mode Decomposition (VMD), Transfer Component Analysis (TCA) and Improved Zebra Optimization Algorithm (IZOA). Firstly, VMD is used to decompose the power load data into several subsequences and TCA is used to reduce the dimensionality of factors related to power load. Secondly, IZOA is used to optimize the parameters of TSK-TFS and Subtractive Clustering Algorithm is used to obtain the number of clusters. The data in the source domain is input into the TSK fuzzy system to obtain the predecessor parameters and successor parameters and retain them. The parameters are inherited and the successor parameters are obtained using some of the target domain data. Finally, the predicted values of several subsequences of the test set (another portion of the target domain data) are obtained according to the obtained parameters of the subsequent components and calculated, and the predicted values of each subsequence are superimposed to obtain the predicted values of short-term power load. The experimental results show that the IZOA-VMD-TSK-TFS-TCA short-term power load prediction model has higher prediction accuracy. Statistical tests also confirmed that the model had superior predictive performance.

1. 引言

电力系统的管理和运行依赖于电力负荷的预测 [1] 。对电力负荷准确的预测有利于提高电力系统的稳定性，节约资源、降低成本、进行电力管理和实施经济调度计划 [2] 。

电力负荷预测按照时间尺度可分为长期、中期、短期、超短期四类 [3] 。按照预测方法分类，电力负荷预测可以分为数学统计方法和机器学习方法 [4] 。基于数学统计的方法有多元线性回归 [5] 、卡尔曼滤波法 [6] 、时间序列法 [7] 等。基于机器学习的方法主要有支持向量机(Support Vector Machine, SVM) [8] 、极限学习机(Extreme Learning Machine, ELM) [9] 、长短期神经网络(Long Short-Term Memory, LSTM) [10] 和门控循环单元(Gated Recurrent Unit, GRU) [11] 等。数学统计和机器学习方法已经广泛应用于电力负荷预测，由于机器学习生成的模型是黑箱，很难对其进行解释，而模糊系统具有可解释性和强非线性逼近能力，因此本文利用模糊系统对短期电力负荷进行预测。

I-型模糊集合(Type-l Fuzzy Set, T1-FS)最早 [12] 由Zadeh教授于1965年提出，模糊集合论的提出使计算机突破了无法处理模糊概念的禁锢。Zadeh教授 [13] 于1975年提出II-型模糊集合(Type-2 Fuzzy Set, T2-FS)的概念，增强了集合的模糊性，进而提高了处理不确定性问题的能力。在模糊集的基础上提出了模糊系统，模糊系统具有很强的非线性逼近能力和可解释性 [14] 。文献 [15] 利用模糊系统并结合相对湿度和温度数据对短期电力负荷进行预测，仿真实验表明利用模糊系统建立短期电力负荷预测模型预测效果较好。

由于模糊系统中的可调参数会影响其预测性能，利用优化算法可以提高模糊系统的预测效果，斑马优化算法(Zebra Optimization Algorithm, ZOA) [16] 具有较强的寻优能力，且其收敛速度较快。利用斑马优化算法对模糊系统中的可调参数寻优，可以提高模糊系统预测的准确性。

为降低电力负荷序列非线性和非平稳性对模型预测性能的影响，采取信号分解的方法可以有效地降低序列的非平稳性。针对负荷数据的非线性性和非平稳性，赵一鸣等 [17] 首先利用经验模态分解(Empirical Mode Decomposition, EMD)将电力负荷序列分解为若干个子序列，其次利用改进后的粒子群算法(Improved Particle Swarm Optimization, IPSO)优化LSTM的参数，最后将各个序列预测值叠加得到电力负荷的预测值。曹广华等 [18] 利用集合经验模态分解将原始数据分解为若干变量，其次利用樽海鞘群优化算法对反向传播神经网络(Back Propagation, BP)和LSTM的参数寻优，最后分解得到的若干子序列分为两组，分别代入BP和LSTM模型中得到电力负荷的预测值。张未等 [19] 利用变分模态分解(Variational Mode Decomposition, VMD)将负荷数据分解，将各个序列利用LSTM预测，再将各个分量联合天气等特征数据利用轻量级梯度提升机(Light Gradient Boosting Machine, LightGBM)进行负荷预测，仿真实验结果表明经过VMD分解后得到的预测结果更精确。

针对电力负荷数据缺失的问题，迁移学习能够从源域数据中得到有效的信息，从而将其迁移至目标域数据中，有效地解决了目标域信息缺失的问题。针对目标域数据不足导致的模型训练不充分，欧阳福莲等 [20] 首先利用多尺度卷积神经网络(Convolutional Neural Network, CNN)提取特征，其次将特征作为双向长短期网络(Bidirectional Long Short-Term Memory, Bi-LSTM)的输入，并引入注意力机制调节权重，利用源域数据训练模型并保留每层最优的训练权重，最后利用目标域数据对模型进行微调并得到预测值。孙晓燕等 [21] 首先根据源域数据利用聚类算法和GRU提取特征和分类模型，其次利用该模型对目标域中小样本进行提取特征并分类，融合样本相似度和时间遗忘因子，最后根据融合特征和迁移学习策略对目标域中的电力负荷进行预测，仿真实验证实了迁移学习的有效性；姜建国等 [22] 利用参数迁移方法结合1DCNN-LSTM对短期负荷进行预测，仿真实验表明在数据缺失的情况下，迁移学习可以有效地提升模型的预测精度。

针对电力负荷序列的非线性性和数据存在缺失的问题，本文首先利用VMD分解将电力负荷数据分解为若干子序列以降低原序列对模型预测性能的影响，将天气等因素利用迁移成分分析的方法进行降维并将源域和目标域之间的距离减小，其次利用改进的斑马优化算法对0阶L2型TSK迁移模糊系统 [23] 的参数进行寻优并利用减法聚类算法确定聚类数，最后将寻优得到的参数和聚类数输入迁移模糊系统中，并将输入代入此模型训练并进行预测，将预测得到的子序列重组得到预测值。本文以欧洲的两个城市的电力负荷数据集为基础，构建IZOA-VMD-TSK-TFS-TCA短期电力负荷预测模型，经过仿真实验分析，证明本文提出的模型具有较高的精度。

2. 相关方法介绍

2.1. 0阶(Takagi-Sugeuo-Kang) TSK模糊系统

TSK模糊系统具有良好的非线性逼近能力，在0阶TSK模糊系统 [24] 中，利用“If-Then”规则定义模糊系统的规则，第 $i$ 条规则为：

$\begin{array}{l}\text{IF}{x}_{1}is{A}_1^k\wedge x_{2}is{A}_2^k\wedge \cdots \wedge x_{d}is{A}_d^k\\ \text{THEN}{y}^k=f^k\left(x\right)=p_0^k,k=1,\cdots ,K\end{array}$ (1)

经过去模糊化处理后得到0阶TSK模糊系统的输出为：

$y^0={\displaystyle \underset{k=1}{\overset{K}{\sum }}\frac{{\mu }^k\left(x\right)}{{\displaystyle \underset{k\text{'}=1}{\overset{K}{\sum }}{\mu }^{k\text{'}}\left(x\right)}}}{f}^k\left(x\right)={\displaystyle \underset{k=1}{\overset{K}{\sum }}{\tilde{\mu }}^k\left(x\right)}{f}^k\left(x\right)={\displaystyle \underset{k=1}{\overset{K}{\sum }}{\tilde{\mu }}^k\left(x\right)p_0^k}$ (2)

其中

${\mu }^k\left(x\right)={\displaystyle \underset{i=1}{\overset{d}{\prod }}{\mu }_{A_i^k}\left(x_i\right)}$ (3)

使用高斯函数作为隶属度函数，则式(3)中 ${\mu }_{A_i}^k\left(x_i\right)$ 可表示为：

${\mu }_{A_i^k}\left(x_i\right)=\exp \left(\frac{-{\left(x_i-c_i^k\right)}^2}{2{\delta }_i^k}\right)$ ， (4)

其中参数 $c_i^k$ 和 ${\delta }_i^k$ 若采用FCM算法获取，则：

$c_i^k=\frac{{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\mu }_{jk}{x}_{ji}}}{{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\mu }_{jk}}}$ ，(5)

${\delta }_i^k=h{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\mu }_{jk}{\left(x_{ji}-c_i^k\right)}^2/{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\mu }_{jk}}}$ ， (6)

其中 ${\mu }_{jk}$ 为输入向量 $x_j={\left(x_{j1},x_{j2},\cdots ,x_{jd}\right)}^T$ 隶属于第 $k$ 类的隶属度， $h$ 为尺度参数。若令

$x_g={\left({\tilde{\mu }}^1\left(x\right),{\tilde{\mu }}^2\left(x\right),\cdots ,{\tilde{\mu }}^K\left(x\right)\right)}^T$ ， (7)

$p_g={\left(p_0^1,p_0^2,\cdots ,p_0^K\right)}^T$ ， (8)

则式(2)可表示为：

$y^0=p_{{}_g}^T{x}_g$ 。 (9)

通过以上计算，可将0阶TSK型模糊系统的参数学习问题最终转化为线性回归问题。由文献 [15] 可知，0阶L2型TSK模糊系统的优化目标函数为：

$\begin{array}{l}\underset{p_{g,}{\xi }^{+},{\xi }^{-},\epsilon }{{\displaystyle \min }}L\left(p_{g,}{\xi }^{+},{\xi }^{-},\epsilon \right)=\frac{1}{N\tau }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[{\left({\xi }_i^{+}\right)}^2+{\left({\xi }_i^{-}\right)}^2\right]}+\frac{1}{2}\left(p_{{}_g}^T{p}_g\right)+\frac{2}{\tau }\cdot \epsilon \\ s.t.\{\begin{array}{c}{y}_i-p_{{}_g}^T{x}_{gi}<\epsilon +{\xi }_i^{+}\\ p_{{}_g}^T{x}_{gi}-y_i<\epsilon +{\xi }_i^{-}\end{array},\forall i\end{array}$ ，(10)

其中 ${\xi }_i^{+}$ 和 ${\xi }_i^{-}$ 为L2范式惩罚项，利用最优化理论，式(10)的对偶问题可表示为：

$\begin{array}{l}\max -{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left({\lambda }_i^{+}-{\lambda }_i^{-}\right)}}\left({\lambda }_j^{+}-{\lambda }_j^{-}\right)\cdot x_{gi}^T{x}_{gi}-\frac{N\tau }{2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left({\lambda }_i^{+}\right)}^2}-\frac{N\tau }{2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left({\lambda }_i^{-}\right)}^2}+{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\lambda }_i^{+}}\cdot y_i\cdot \tau -{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\lambda }_i^{-}}\cdot y_i\cdot \tau \\ s.t.{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left({\lambda }_i^{+}+{\lambda }_i^{-}\right)}=1,{\lambda }_i^{+}\ge 0,{\lambda }_i^{-}\ge 0,\forall i\end{array}$ (11)

令

$\tilde{\alpha }={\left({\tilde{\alpha }}_1,\cdots ,{\tilde{\alpha }}_{2N}\right)}^T={\left({\left({\lambda }^{+}\right)}^T,{\left({\lambda }^{-}\right)}^T\right)}^T$ ， (12)

$z_i=\{\begin{array}{ll}{x}_{gi},\hfill & i=1,2,\cdots ,N,\hfill \\ -x_{g\left(i-N\right)},\hfill & i=N+1,\cdots ,2N.\hfill \end{array}$ ， (13)

$\beta ={\left(\tau y^T,-\tau y^T\right)}^T,y=\left(y_1,\cdots ,y_N\right)$ ， (14)

通过式(12)~(14)，可将式(11)简化为标准二次规划的形式：

$\begin{array}{l}\underset{\tilde{\alpha }}{{\displaystyle \text{argmax}}}-{\tilde{\alpha }}^T\tilde{K}\tilde{\alpha }+{\tilde{\alpha }}^T\beta \\ s.t.{\tilde{\alpha }}^{T}1=1,{\alpha }_i\ge 0,\forall i\end{array}$ ， (15)

其中 $\tilde{K}={\left[{\tilde{k}}_{ij}\right]}_{2N\times 2N},{\tilde{k}}_{ij}=z_i^T{z}_i+\frac{N\tau }{2}{\delta }_{ij}$ 为核矩阵， ${\delta }_{ij}=\{\begin{array}{c}1,i=j\\ 0,i\ne j\end{array}$ 。

通过对(15)的求解得到对偶问题的解，则原问题中各个参数的解为：

$p_g=\frac{2}{\tau }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left({\lambda }_i^{+}-{\lambda }_i^{-}\right)}{x}_{gi}$ ， (16)

${\xi }_i^{+}=N{\lambda }_i^{+}$ ， (17)

${\xi }_i^{-}=N{\lambda }_i^{-}$ ， (18)

$\epsilon ={\displaystyle \underset{i=1}{\overset{N}{\sum }}\left({\lambda }_i^{+}-{\lambda }_i^{-}\right)}{y}_i-\frac{N}{2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[{\left({\lambda }_i^{+}\right)}^2-{\left({\lambda }_i^{-}\right)}^2\right]}-\frac{1}{\tau }{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left({\lambda }_i^{+}-{\lambda }_i^{-}\right)}\left({\lambda }_j^{+}-{\lambda }_j^{-}\right)}\cdot x_{{}_{gi}}^T{x}_{gi}$ ， (19)

式(16)即为0阶L2型TSK模糊系统后件参数的最优解。

2.2. 斑马优化算法及改进

2.2.1. 斑马优化算法

斑马优化算法(ZOA) [16] 通过模拟斑马的觅食行为和针对捕食者的防御行为进行寻优，斑马优化算法寻优能力强且收敛速度较快，ZOA主要分为3个阶段，首先初始化寻优空间的斑马种群 $X={\left\{X_i\right\}}_{i=1,\cdots ,N}={\left\{x_{i,j}\right\}}_{j=1,\cdots ,m}$ ，

$x_{i,j}=l{b}_j+r\left(u{b}_j-l{b}_j\right)$ ， (20)

其中 $X_i$ 为第 $i$ 个斑马的值， $x_{i,j}$ 为第 $i$ 个斑马的第 $j$ 维度值，第 $i$ 个斑马对应的目标函数值为 $F_i$ ，斑马 $l{b}_j$ 为寻优下界， $u{b}_j$ 为寻优上界， $r$ 为[0, 1]之间的随机数。其次是斑马觅食阶段，位置更新公式为：

$x_{i,j}^{new,P1}=x_{i,j}+r\left(P{Z}_j-I\times x_{i,j}\right)$ ， (21)

$X_i=\{\begin{array}{ll}{X}_i^{new,P1},\hfill & F_i^{new,P1}<F_i,\hfill \\ X_i,\hfill & else.\hfill \end{array}$ ， (22)

其中 $X_i^{new,P1}$ 为第 $i$ 个斑马在觅食阶段的值， $x_{{}_{i,j}}^{new,P1}$ 为第 $i$ 个斑马在觅食阶段第 $j$ 维度值， $r$ 为[0, 1]之间的随机数， $I$ 为属于集合{1, 2}的随机值， $P{Z}_j$ 为先锋斑马第 $j$ 维度值， $F_i^{new,P1}$ 为第 $i$ 个斑马在觅食阶段的目标函数值，最后是针对捕食者的防御阶段，斑马的防御策略根据捕食者的分为两种，若捕食者为狮子，斑马选择逃跑战略即 $P_S\le 0.5$ 对应的函数式；若捕食者为其他的动物，则斑马采取攻击捕食者的策略即 $else$ 对应的函数式。斑马根据(23)式和式(24)更新位置，

$x_{i,j}^{new,P2}=\{\begin{array}{ll}{x}_{i,j}+R\left(2r-1\right)\left(1-\frac{t}{T}\right)x_{i,j},\hfill & P_s\le 0.5,\hfill \\ x_{i,j}+r\left(A{Z}_j-I\cdot x_{i,j}\right),\hfill & else.\hfill \end{array}$ ， (23)

$X_i=\{\begin{array}{l}{X}_i^{new,P2},F_i^{new,P2}<F_i\hfill \\ X_i,else.\hfill \end{array}$ ， (24)

其中 $X_i^{new,P2}$ 为第 $i$ 个斑马在防御阶段的值， $x_{{}_{i,j}}^{new,P2}$ 为第 $i$ 个斑马在防御阶段第 $j$ 维度值， $A{Z}_j$ 为被攻击斑马的第 $j$ 维度值， $F_i^{new,P2}$ 为第 $i$ 个斑马在防御阶段的目标函数值， $t$ 为迭代次数， $T$ 为最大迭代次数， $R$ 为常数0.01， $P_s\in [0,1]$ 为在两种策略中的选择一种策略的概率。

2.2.2. 改进的斑马优化算法

为了更好地得到优化算法的寻优结果，本文引入莱维飞行对斑马优化算法的种群初始化，在防御阶段利用正态分布数生成两种策略的切换概率。

1) 初始种群会很大程度地影响优化算法的寻优速度和结果，利用莱维飞行初始化斑马种群的寻优空间，使得寻优空间更具有随机性，增加全局搜索能力，即：

$s=\frac{u}{|v{|}^{1/\beta }}$ ， (25)

其中 $u\sim N\left(0,{\sigma }^2\right),v\sim N\left(0,1\right),\sigma ={\left\{\frac{\Gamma \left(1+\beta \right)\sin \left(\frac{\pi \beta }{2}\right)}{\beta \Gamma \left(\frac{1+\beta }{2}\right)2^{\frac{\beta -1}{2}}}\right\}}^{\frac{1}{\beta }}$ ， $\beta =1.5$ ，则初始化种群为：

$x_{i,j}=l{b}_j+\left(s_{i,j}-\left[s_{i,j}\right]\right)\left(u{b}_j-l{b}_j\right)$ ， (26)

其中 $s_{i,j}$ 为莱维飞行生成的步长， $[\cdot ]$ 为取整符号。

2) 第二个改进是防御阶段中两种策略的切换概率，由于两种策略分别是针对狮子与其他捕食者，而遇见其他捕食者的概率会更大，因此设置 $P_S\sim N\left(0,1\right)且P_S=\{\begin{array}{c}1,P_S\ge 1\\ 0,P_S\le 0\\ P_{S,}else\end{array}$ ，且将式(23)改为：

$x_{i,j}^{new,P2}=\{\begin{array}{ll}{x}_{i,j}+R\left(2r-1\right)\left(1-\frac{t}{T}\right)x_{i,j},\hfill & P_s>0.5,\hfill \\ x_{i,j}+r\left(A{Z}_j-I\cdot x_{i,j}\right),\hfill & else.\hfill \end{array}$ (27)

2.3. 变分模态分解

2.3.1. 变分模态分解

变分模态分解(VMD) [25] 是一种自适应、完全递归的一种信号分解的方法，且变分模态的分解数可以指定，构造和求解变分问题是变分模态分解的核心，首先构造变分问题：

$\{\begin{array}{l}\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{\min }\left\{{\displaystyle \underset{k=1}{\overset{K}{\sum }}|\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\cdot u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\Vert {|}_2^2}\right\}\hfill \\ s.t.{\displaystyle \underset{k=1}{\overset{K}{\sum }}{u}_k=f\left(t\right)}\hfill \end{array}$ ， (28)

其中 $K$ 为模态分解数， $u_k$ 为各模态函数， ${\omega }_k$ 为各模态中心频率， $f\left(t\right)$ 为原始信号，其次对变分问题求解，引入拉格朗日乘子 $\lambda$ 和二次惩罚项得到求解公式：

$L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \right)=\alpha {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\cdot u_k\left(t\right)\right]e^{-j{\omega }_{k}t}\Vert }_2^2}+{\Vert f\left(t\right)-{\displaystyle \underset{k=1}{\overset{K}{\sum }}{u}_k\left(t\right)}\Vert }_2^2+\langle \lambda \left(t\right),f\left(t\right)-{\displaystyle \underset{k=1}{\overset{K}{\sum }}{u}_k\left(t\right)}\rangle$ ， (29)

其中 $\alpha$ 为惩罚项，对变分问题求解上式的鞍点可得：

${\hat{u}}_k^{n+1}\left(\omega \right)=\frac{\hat{f}\left(\omega \right)-{\displaystyle \underset{i\ne k}{\sum }{\hat{u}}_i\left(\omega \right)+\hat{\lambda }\left(\omega \right)/2}}{1+2\alpha {\left(\omega -{\omega }_k\right)}^2}$ ， (30)

${\omega }_k^{n+1}=\frac{{\int }_0^{\infty }\omega {\left|{\hat{u}}_k\left(\omega \right)\right|}^{2}d\omega }{{\int }_0^{\infty }{\left|{\hat{u}}_k\left(\omega \right)\right|}^{2}d\omega }$ ， (31)

${\hat{\lambda }}^{n+1}\left(\omega \right)={\hat{\lambda }}^n\left(\omega \right)+\tau \left(\hat{f}\left(\omega \right)-{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\hat{u}}_k^{n+1}\left(\omega \right)}\right)$ ， (32)

当满足下式时，则求解完成，

${\displaystyle \underset{k=1}{\overset{K}{\sum }}{\Vert {\hat{u}}_k^{n+1}-{\hat{u}}_k^n\Vert }_2^2/{\Vert {\hat{u}}_k^n\Vert }_2^2}<\epsilon$ . (33)

2.3.2. 最小样本熵

由于变分模态分解的分解个数可以随机选定，样本熵可以测量时间序列的复杂度，具有一致性，为了得到更有效的分解成分，利用最小样本熵确定分解个数，最小样本熵的计算步骤如下：

step1：将时间序列 $X=\left\{x_1,x_2,\cdots ,x_N\right\}$ 构成 $l$ 维矢量，即：

$X_i=\left\{x_i,x_{i+1},\cdots ,x_{i+l-1}\right\},i=1,2,\cdots ,N-l+1.$ (34)

step2：定义距离 $dist\left(X_i,X_j\right)$ 为：

$dist\left(X_i,X_j\right)=\underset{k\in \left(0,l-1\right)}{\max }\left|x_{i+k}-x_{j+k}\right|$ ， (35)

step3：给定阈值 $r$ ，统计 $dist\left(X_i,X_j\right)<r$ 的数量并计算矢量总个数 $N-l$ 的比值，即：

$B_i^l\left(r\right)=\frac{num\left\{dist\left(X_i,X_j\right)<r\right\}}{N-l}$ ， (36)

step4：对所有由式(36)得到的结果求平均，即：

$B^l\left(r\right)=\frac{1}{N-l+1}{\displaystyle \underset{i=1}{\overset{N-l+1}{\sum }}{B}_i^l\left(r\right)}$ ， (37)

则样本熵的估计值为：

$Sampen\left(l,r,N\right)=-\ln \left(\frac{B^{l+1}\left(r\right)}{B^l\left(r\right)}\right)$ ， (38)

则对于VMD分解不同的分解数 $k$ ，计算最小样本熵为：

$Sampe{n}_{\min }=\underset{k}{\min }Sampen\left(l,r,N\right)$ ， (39)

利用上式即可找到最小样本熵趋于稳定对应分解数，即为VMD的分解模态数。

2.4. 迁移成分分析

迁移成分分析(TCA) [26] 在不改变数据结构和域间的分布变化的情况下，能够减小源域数据和目标域数据之间差异，使得不同域之间的差异在映射后的子空间中降低，且可以将数据进行降维。假设存在一个特征映射 $\varphi$ ，使得 $P\left(\varphi \left(X_s\right)\right)\approx P\left(\varphi \left(X_t\right)\right)$ ，迁移成分分析将最大均值差异(Maximize Mean Discrepancy, MMD)作为度量准则，使不同域之间的分布差异最小化，最大均值差异的计算公式为：

$D\left(X_s,X_t\right)={\Vert \frac{1}{n_1}{\displaystyle \underset{i=1}{\overset{n_1}{\sum }}\varphi \left(x_i\right)-}\frac{1}{n_2}{\displaystyle \underset{j=1}{\overset{n_2}{\sum }}\varphi \left(x_j\right)}\Vert }_H^2$ ， (40)

其中 $X_s$ 为源域数据， $X_t$ 为目标域数据， $n_1$ 和 $n_2$ 分别为源域和目标域的样本个数， $\left|\right|\cdot |{|}_H$ 为再生核希尔伯特空间范数，利用核矩阵 $K$ 和度量矩阵 $L$ 可将最大均值差异计算公式转化为：

$D\left(X_s,X_t\right)=tr\left(KL\right)$ ， (41)

其中

$K=\left[\begin{array}{cc}{K}_{s,s}& K_{s,t}\\ K_{t,s}& K_{t,t}\end{array}\right]\in R^{\left(n_1+n_2\right)\times \left(n_1+n_2\right)}$ ， (42)

$L_{i,j}=\{\begin{array}{ll}\frac{1}{n_1{}^2},\hfill & x_i,x_j\in D_s,\hfill \\ \frac{1}{n_2{}^2},\hfill & x_i,x_j\in D_t,\hfill \\ -\frac{1}{n_1{n}_2},\hfill & else.\hfill \end{array}$ (43)

对式(42)中的矩阵 $K$ 规格减小，得到 $\tilde{K}$ ：

$\tilde{K}=\left(K{K}^{-1/2}\tilde{W}\right)\left({\tilde{W}}^T{K}^{-1/2}K\right)=KW{W}^{T}K$ ， (44)

其中 $\tilde{W}\in R^{\left(n_1+n_2\right)\times m}$ ， $W=K^{-1/2}\tilde{W}$ ，将距离公式转化为 $D\left(X_s,X_t\right)=tr\left(W^{T}KLKW\right)$ ，将其代入式(40)中，则TCA的目标函数为：

$\begin{array}{l}\underset{W}{\min }tr\left(W^{T}KLKW\right)+\mu tr\left(W^{T}W\right)\\ s.t.W^{T}KHKW=I_m\end{array}$ ， (45)

其中 $H=I_{n_1+n_2}-1/\left(n_1+n_2\right){11}^T$ ， $I$ 为单位矩阵，1为 $n_1+n_2$ 行的全1列向量， $\mu$ 为权衡因子， $tr\left(W^{T}W\right)$ 为正则项。

3. 基于IZOA-VMD-TSK-TFS-TCA的预测模型建立

3.1. 0阶L2型TSK迁移模糊系统

利用迁移学习策略可以有效地解决由于数据缺失导致的模型预测性能降低问题，TSK模糊系统具有较好的非线性逼近性且具有可解释性，将迁移策略嵌入至0阶L2型TSK模糊系统中，可以提高模型预测的准确性，0阶L2型TSK迁移模糊系统根据文献 [22] 可知引入迁移机制的目标函数为：

$\begin{array}{l}\underset{p_g,{\xi }^{+},{\xi }^{-},\epsilon }{{\displaystyle \min }}\tilde{L}\left(p_g,{\xi }^{+},{\xi }^{-},\epsilon \right)={\Theta }_c\left(p_g,{\xi }^{+},{\xi }^{-},\epsilon \right)+\lambda {\Theta }_h\left(p_g,p_{g_0}\right)\\ s.t.\{\begin{array}{c}{y}_i-p_{{}_g}^T{x}_{gi}<\epsilon +{\xi }_i^{+}\\ p_{{}_g}^T{x}_{gi}-y_i<\epsilon +{\xi }_i^{-}\end{array},\forall i\end{array}$ ， (46)

${\Theta }_c\left(p_g,{\xi }^{+},{\xi }^{-},\epsilon \right)=\frac{1}{N\tau }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left[{\left({\xi }_i^{+}\right)}^2+{\left({\xi }_i^{-}\right)}^2\right]}+\frac{1}{2}\left(p_{{}_g}^T{p}_g\right)+\frac{2}{\tau }\cdot \epsilon$ ， (47)

${\Theta }_h\left(p_g,p_{g_0}\right)={\left(p_g-p_{g_0}\right)}^T\left(p_g-p_{g_0}\right)$ ， (48)

其中 ${\Theta }_c$ 为利用当前场景数据对模型训练学习， ${\Theta }_h$ 为利用历史场景迁移学习， $\lambda$ 为平衡当前场景和历史场景的参数。根据文献 [22] 求解上述问题的对偶问题即可得到原问题的解，即：

$p_g=\frac{2\lambda }{1+2\lambda }{p}_{g_0}+\frac{1}{1+2\lambda }{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left({\alpha }_i-{\alpha }_i^{\ast }\right)}{x}_{gi}$ ， (49)

其中 ${\alpha }_i$ 和 ${\alpha }_{{}_i}^{\ast }$ 为对偶问题的解， $p_{g_0}$ 为迁移原始场景中的后件参数。

3.2. IZOA-VMD-TSK-TFS-TCA模型

本文提出的IZOA-VMD-TSK-TFS-TCA模型具体步骤如下，令源域数据 $D_h=\left\{X_h;Y_h\right\}$ ， $X_h=\left[X_{h1},X_{h2}\cdots ,X_{h6}\right]$ ，目标域数据为 $D_m=\left\{X_m;Y_m\right\}$ ， $X_m=\left[X_{m1},X_{m2}\cdots ,X_{m6}\right]$ 。

step1：首先将短期电力负荷原始信号 $Y=\left[Y_h,Y_m\right]$ 进行变分模态分解为 $K$ 个若干子序列得到 $Y^{\prime }=\left[{Y^{\prime }}_{h1},{Y^{\prime }}_{m1};{Y^{\prime }}_{h2},{Y^{\prime }}_{m2}\cdots ;{Y^{\prime }}_{hK},{Y^{\prime }}_{mK}\right]$ ，选取3个天气因素和3个日期因素共6维即 $X=\left[X_h,X_m\right]$ ；

step2：利用迁移成分分析将原始数据中源域和目标域的天气因素和日期因素进行降维并将两个域之间的距离减小得到 $X^{\prime }=\left[{X^{\prime }}_h,{X^{\prime }}_m\right]$ 为3维数据；

step3：利用 $Y^{\prime }$ 和 $X^{\prime }$ 分别利用减法聚类算法得到聚类个数 $M=\left[M_1,M_2\cdots ,M_K\right]$ ；

step4：将 $Y^{\prime }$ 和 $X^{\prime }$ 划分为训练集和测试集；

step5：根据step4中的训练集数据利用改进后的斑马优化算法对0阶L2型TSK迁移模糊系统的参数进行寻优；

step6：将最优参数代入0阶L2型TSK模糊系统中，利用step4中训练集对TSK模糊系统训练，保留后件参数 $p_{g_0}$ ；

step7：继承前件参数和后件参数并对step4中测试集进行预测；

step8：将得到的各个预测值序列重组叠加，即可得到测试集上短期电力负荷的预测值。

IZOA-VMD-TSK-TFS-TCA模型流程图如图1所示，图1展示分为4个部分，第一部分为对数据进行处理，即利用VMD对短期电力负荷数据进行分解，利用TCA对天气和日期因素降维，且降低源域和目标域之间的差异；第二部分根据第一部分处理得到的数据采用减法聚类得到聚类数，即是TSK模糊系统的前件参数 $M$ ；第三部分为利用改进后的斑马优化算法对模糊迁移系统中的参数寻优；第四部分为利用第二和三部分得到的参数，根据第一部分得到的数据，采用迁移学习的策略对测试集上的短期电力负荷预测。

Figure 1. IZOA-VMD-TSK-TFS-TCA model flowchart

图1. IZOA-VMD-TSK-TFS-TCA模型流程图

4. 仿真实验

4.1. 数据描述

本文的数据来源于欧洲国家的电力负荷数据(https://data.open-power-system-data.org/)，数据记录每小时电力负荷，本文选取爱沙尼亚(EE)和拉托维亚(LV)的短期电力负荷数据，数据集选取拉托维亚2018年4月1日至2018年5月31日和爱沙尼亚2018年4月4日至2018年4月18日的短期电力负荷数据，其中爱沙尼亚在时间段4月9日12时至4月10日5时存在数据缺失。拉托维亚的原始数据如图2和图3所示，将爱沙尼亚作为源域，拉托维亚作为目标域，选取的天气因素包括气温、直接辐射和散射辐射，日期因素包括第几个小时、一周中的第几天、是否为休息日(休息日为1，非休息日为0)。划分训练集和测试集，训练集为源域的数据和目标域中2018年4月4日至2018年4月11日期间的数据，测试集为爱沙尼亚2018年4月12日至2018年4月18日。爱沙尼亚电力负荷原始数据如图2所示，拉托维亚电力负荷原始数据如图3所示。

Figure 2. LV power load raw data

图2. LV电力负荷原始数据

Figure 3. EE power load raw data

图3. EE电力负荷原始数据

表1展示了爱沙尼亚(源域数据)和拉托维亚(目标域数据)的最大值、最小值、均值、标准差、偏度和峰度。

Table 1. Statistical characteristics of raw data

表1. 原始数据的统计特性

4.2. VMD分解电力负荷序列分析

利用样本熵可以找到不同的分解个数对应的最小样本熵，再提取得到的最小样本熵的最小值，最小值对应的分解个数即为最优的分解个数，不同的分解个数对应的最小样本熵的最小值如图4所示：

Figure 4. The minimum sample entropy corresponding to different decomposition numbers

图4. 不同分解个数对应的最小样本熵

从图4可知，当分解个数达到7时，对应的最小样本熵趋于稳定，因此选取分解序列的个数为7。VMD电力负荷分解序列如图5所示。

Figure 5. VMD power load decomposition sequence diagram

图5. VMD电力负荷分解序列图

4.3. IZOA测试

为了检验改进的斑马优化算法的收敛性，利用文献 [16] 中的单峰目标函数(F1~F7)、多峰目标函数(F8~F13)和固定维的多峰目标函数(F14~F23)作测试。为了验证IZOA的性能，将其与ZOA和PSO作对比，设定算法的参数如下，

1) 设定IZOA和ZOA的种群数为30，最大迭代数为500；

2) 设定PSO的种群数为30，最大迭代数为500，学习因子设置为1.5和2，选定30次测试结果的均值和标准差作为评价指标，测试结果如表2所示：

Table 2. Test results of IZOA, ZOA, and PSO on the benchmark function

表2. IZOA、ZOA和PSO在基准函数上测试结果

从表3中可以看到IZOA在求解目标函数F1~F4、F7、F12、F14和F19~F22具有较好的效果；在F6、F9~F11上与ZOA效果相同而PSO效果较差；在F16~F17上均值相同，标准差较大；在F5上为标准差较小，而均值较大；在F13上IZOA的均值较小，标准差较大；F8、F15、F18和F23上IZOA得到的均值和方差更大。

4.4. 模型评价指标

为了验证模型的预测性能，利用以下3个统计量评价模型，分别为平均绝对百分误差 $\left(MAPE\right)$ 、平均误差 $\left(MAE\right)$ 以及均方根误差 $\left(RMSE\right)$ ：

$MAPE=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\frac{|y_i-{\hat{y}}_i|}{y_i}}\times 100\%$ ， (50)

$MAE=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left|y_i-{\hat{y}}_i\right|}$ ， (51)

$RMSE=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(y_i-{\hat{y}}_i\right)}^2}}{N}}$ ， (52)

其中 $y_i$ 为第 $i$ 个时刻的电力负荷的真实值， $y_i$ 为第 $i$ 个时刻的电力负荷的预测值， $N$ 为测试集中样本的总个数，以上3个指标的值越小，则表明模型的预测性能越优。

4.5. 预测结果分析

首先利用减法聚类确定迁移模糊系统的聚类数，再将数据做归一化处理，使得数据分布在[0, 1]，归一化处理的公式为

$x_{ij}=\frac{{x^{\prime }}_{ij}-x_{\min }}{x_{\max }-x_{\min }}$ ， (53)

原始数据 $X=\left\{{x^{\prime }}_{ij}\right\}$ 其中 $y_{\max }$ 为1， $y_{\min }$ 为0， $x_{\max }$ 和 $x_{\min }$ 分别为 ${x^{\prime }}_{ij}$ 所在列的最大值和最小值， $x_{ij}$ 为归一化后的数据，设置斑马优化算法的初始值为4，最大迭代次数为20；迁移成分分析的正则化参数设为0.6，将天气等因素降至3维，选定高斯核函数且带宽设置为0.4，为了验证本文提出模型的预测精度，利用以下几种模型作对比，分别为：

1) 迁移和迁移成分分析的情况下利用ZOA优化算法结合VMD分解、EEMD分解和不分解模型；利用PSO算法结合VMD模型；只利用VMD分解模型，用以上模型作对比；

2) 迁移情况下利用VMD分解、EEMD分解和不分解作对比；

3) 不迁移的情况下利用VMD分解和不分解作对比；

4) 融合VMD分解、LSTM和TCA利用源域和目标域中所有数据对测试集数据进行预测，融合VMD和LSTM以及单独用LSTM利用目标域数据对测试集数据进行预测。

1)中的模型预测值与真实值的对比如图6所示，可以看出在类型1的模型中IZOA-VMD-TSK-TFS-TCA短期电力负荷模型预测性能更优。2)中模型预测值与真实值的对比如图7所示，3)中模型预测值与真实值的对比如图8所示，4)中模型预测值与真实值的对比如图9所示，在不进行迁移学习，直接采用有数据缺失的目标域数据集对测试集预测，预测效果并不理想，采用迁移学习策略而未对数据进行迁移成分分析，预测效果也不理想。可以看出迁移学习策略和TCA在目标域数据缺失的情况下可以有效地提高预测精度。图10~13通过箱线图的形式展示了各个模型误差值，利用箱线图可以得到误差的四分位数和数据集中程度。图10中IZOA-VMD-TSK-TFS-TCA误差的中位数更接近于0，表明平均水平更接近于0，且比其余模型的误差数据更集中在0附近。

Figure 6. Short-term power load forecasting chart based on decomposition, migration, and TCA combination strategies

图6. 分解、迁移和TCA组合策略的短期电力负荷预测图

Figure 7. Short-term power load forecasting chart based on decomposition and migration combination strategies

图7. 分解、迁移组合策略的短期电力负荷预测图

Figure 8. Short term electricity load forecasting chart without migration

图8. 不进行迁移的短期电力负荷预测图

Figure 9. Short term power load prediction chart of LSTM model

图9. LSTM模型的短期电力负荷预测图

Figure 10. Error analysis of short-term power load using decomposition, migration, and TCA combination strategies

图10. 分解、迁移和TCA组合策略短期电力负荷的误差

Figure 11. Error analysis of short-term power load based on decomposition and migration combination strategies

图11. 分解、迁移组合策略的短期电力负荷的误差

Figure 12. Error in short-term electricity load without migration

图12. 不进行迁移的短期电力负荷的误差

Figure 13. Short term power load prediction chart of LSTM model

图13. LSTM模型的短期电力负荷预测图

表3展示了不同模型做预测的结果，可以看出进行迁移并对数据进行VMD并利用改进斑马优化算法和迁移成分分析时，预测效果最好，ZOA-TSK-TFS与ZOA-TSK-FS模型比较，采用迁移学习策略时MAPE降低了0.23%，ZOA-TSK-TFS与ZOA-TSK-TFS-TCA相比，采用TCA使得MAPE减少1.35%，ZOA-TSK-TFS-TCA与ZOA-VMD-TSK-TFS-TCA相比，对电力负荷序列分解后模型的MAPE减少了1.25%，对斑马优化算法改进后与ZOA-VMD-TSK-TFS-TCA相比，MAPE减少了0.11%。经过迁移策略、VMD、改进后的斑马优化算法和迁移成分分析可以有效的提高模型的预测精度，使得IZOA-VMD-TSK-TFS-TCA模型对短期电力负荷的预测更准确。

Table 3. Statistical indicator values for each model

表3. 各个模型的统计指标值

4.6. 模型的检验

为了定量的分析本文提出的IZOA-VMD-TSK-TFS-TCA短期电力负荷模型与其余模型的预测性能的差异，利用指标改善率( $I{R}_j$ )，其定义为

$I{R}_1=\frac{MAP{E}_i-MAP{E}_0}{MAP{E}_i}\times 100\%,i=1,2,\cdots ,13$ ， (54)

$I{R}_2=\frac{MA{E}_i-MA{E}_0}{MA{E}_i}\times 100\%,i=1,2,\cdots ,13$ ， (55)

$I{R}_3=\frac{RMS{E}_i-RMS{E}_0}{RMS{E}_i}\times 100\%,i=1,2,\cdots ,13$ ， (56)

其中 $MAP{E}_i$ 、 $MA{E}_i$ 和 $RMS{E}_i$ 为本文作为对比模型的统计指标值， $MAP{E}_0$ 、 $MA{E}_0$ 和 $RMS{E}_0$ 为本文提出的IZOA-VMD-TSK-TFS-TCA模型的统计指标值， $IR$ 结果如表4所示，

Table 4. Sensitivity of each model to IZOA-VMD-TSK-TFS-TCA

表4. 各个模型与IZOA-VMD-TSK-TFS-TCA的敏感度

从表4中可以看出 $I{R}_1$ 、 $I{R}_2$ 和 $I{R}_3$ 的范围分别是[0.0686 0.7189]、[0.0843 0.7011]和[0.09120.7097]，本文提出的模型较ZOA-LSTM可以较大改善模型的 $MAPE$ 、 $MAE$ 和 $RMSE$ ；而较ZOA-VMD-TSK-TFS-TCA本文提出的模型则对三个统计指标改善较小，从表4可以看出本文提出的模型与对比模型相较，均在一定的程度上改善了预测模型的三个统计指标值。

5. 结束语

本文针对短期电力负荷序列的非线性性和数据缺失导致预测精度下降的问题，提出IZOA-VMD-TSK-TFS-TCA短期电力负荷预测模型，通过实验结果可得到以下结论：

1) 本文利用VMD对短期电力负荷数据进行分解，降低了短期电力负荷的非线性和非平稳性，同时采用TCA对天气和日期等因素的数据进行降维且减小源域和目标域之间的距离。对ZOA算法进行改进，并利用基准函数进行测试，结果表明IZOA算法的性能优于ZOA算法和PSO算法。

2) 利用IZOA算法对TSK-TFS的可调参数进行寻优，避免盲目设置参数使得模型的预测性能下降。针对目标域数据缺失的问题，本文利用迁移学习的策略，利用TSK-TFS模型将源域的有效信息迁移至目标域中，解决了由于目标域中数据缺失导致的预测精度降低的问题。

3) 将IZOA-VMD-TSK-TFS-TCA模型与ZOA-VMD-TSK-TFS-TCA、ZOA-EEMD-TSK-TFS-TCA、PSO-VMD-TSK-TFS-TCA、VMD-TSK-TFS-TCA、ZOA-TSK-TFS-TCA、ZOA-VMD-TSK-TFS、ZOA-EEMD-TSK-TFS、ZOA-TSK-TFS、ZOA-VMD-TSK-FS和ZOA-TSK-FS、ZOA-VMD-LSTM-TCA、ZOA-VMD-LSTM和ZOA-LSTM模型的实验结果作对比，仿真结果表明，本文提出的模型预测精度更高。同时利用改善率对模型进行检验，结果表明IZOA-VMD-TSK-TFS-TCA短期电力负荷模型的预测性能更优。

基金项目

东华理工大学校级项目(DHYC-202338)。