基于SVR与libsvm的电力负荷预测Matlab实现——时间序列预测实战
2026/9/16 5:07:43 网站建设 项目流程

简介:一套基于支持向量机回归的电力负荷预测模型,以Matlab为实现平台,主要面向电力系统时间序列预测的初学者,也适合需要在Matlab中快速完成回归模型的开发者。整体方案围绕单变量时间序列数据展开,集成了拟合优度、平均绝对误差、平均偏差误差等多维评价指标,代码质量较高;负荷数据以CSV文件提供,替换数据后即可直接运行,运行环境需在2020版本及以上。资源包为ZIP压缩包,共含八个文件,整体大小约985千字节,其中包括主程序源代码、适用于六十四位Windows系统的编译组件、示例负荷数据、算法参数说明文本,以及三张预测效果对比图片,能帮助使用者省去编译配置和参数查阅的流程。目前已有四百零八人浏览学习。借助该包可掌握从数据读取、模型训练、指标计算到结果可视化的完整回归预测流程,并可在现有框架内灵活修改数据或调整参数,开展新的电力负荷预测实验。

1. 先把电力负荷预测这件事说清楚

电力负荷预测本质上是把一个时间序列的未来值,用历史值“算”出来。传统回归模型处理这类问题,往往被两个问题卡住:一是负荷数据不是纯线性的,气温、节假日、生产作息会叠加出明显的非线性段;二是样本量不大,深度学习模型动辄几千上万的样本需求在这里并不现实。SVR(支持向量回归)正好落在两者之间,它用核函数把低维非线性映射到高维线性空间,在小样本下依然能控制泛化误差,所以在电力系统短期负荷预测里一直是性价比很高的基线模型。这套 Matlab 实现基于 libsvm 工具包完成 SVR 回归,输入一份 csv 格式的历史负荷序列,输出未来时段的预测值,并给出 R²、MAE、MBE 等评价指标。适合正在做时间序列预测课程设计、需要快速验证 SVR 效果的工程师,也适合想对比“传统统计模型 vs 机器学习模型”的算法岗同学。需要说明的是,它做的是单变量预测,也就是只用负荷自身的历史数据,不掺温度等外生变量,这也是很多论文里最基本的对比设置。下面直接从数据、代码、参数到坑位,完整过一遍。

2. 数据组织与单变量时间序列的样本构造

2.1 csv 里到底存了什么

data.csv 是这份资源的核心数据文件,用 Excel 或文本编辑器打开就能看到典型的两列结构:第一列是时间点,第二列是对应时刻的负荷值。采样间隔常见的是 15 分钟、半小时或 1 小时,具体间隔需要打开文件确认,因为后续构造特征时,lag的含义完全取决于这个间隔。用 Matlab 读取并预览数据:

% 读取 csv 数据 dataTable = readtable('data.csv'); head(dataTable, 10) % 预览前10行 % 取出负荷列,转为 double 数组 loadSeries = dataTable{:, 2};

readtable在 Matlab R2020 及以上版本里对 csv 的兼容性很好,不需要额外指定分隔符。这里把数据转成列向量是为了适配后面的样本矩阵构造。如果 csv 第一行是列名,dataTable{:, 2}会正确取出数值列;如果没有列名,第二列就是原始负荷。

2.2 为什么单变量预测要把序列切成“X 看过去,y 看未来”

单变量时间序列预测的第一步不是直接丢给 SVR,而是把一列连续序列转换成监督学习所需的样本对。假设我们用前p个时刻的负荷值预测第p+1个时刻,即自回归阶数为p,那么原始序列就会被滑窗切成若干组:

function [X, y] = createLagMatrix(series, p) n = length(series); % 样本数:前 p 个历史值至少留 p 个,所以是 n - p X = zeros(n - p, p); y = zeros(n - p, 1); for i = 1:n - p X(i, :) = series(i : i + p - 1)'; y(i) = series(i + p); end end

调用方式:

p = 7; % 用前7个时刻预测下一个时刻 [X, y] = createLagMatrix(loadSeries, p);

这段代码中,series(i : i + p - 1)是长度为p的滑动窗口,y(i)是窗口后紧邻的负荷值。参数p的选择非常关键:取得太小,模型看不到足够长的趋势;取得太大,特征维度增加而样本数减少,SVR 的训练时间线性上升,还可能引入噪声。常见做法是先画自相关图,看延迟几阶的相关系数明显不为零,再结合业务周期来定。负荷数据通常有日周期性,如果采样间隔是 1 小时,p = 24是合理的起点。

2.3 训练集 / 测试集划分要防“泄露”

时间序列划分不能随机打乱,否则未来信息进入训练集,评价指标会虚高。正确做法是按时序切分:

% 用前80%的数据训练,后20%测试 numSamples = length(y); trainNum = floor(numSamples * 0.8); X_train = X(1:trainNum, :); y_train = y(1:trainNum); X_test = X(trainNum+1:end, :); y_test = y(trainNum+1:end);

这里没有使用交叉验证,而是固定切分,原因在于时间序列相邻样本之间存在强相关性,K 折随机打乱会破坏时序依赖,导致验证结果失真。如果要更严谨,可以用时间序列交叉验证,即扩展窗口方式,但作为基线预测模型,固定切分已经足够说明问题。

3. SVR 原理与 libsvm 在 Matlab 中的接入

3.1 SVR 和普通回归的差别在哪

SVR 不是让所有样本点都落在回归线上,而是允许预测值与真实值之间存在一个大小为epsilon的误差管道。损失函数只惩罚超出管道的样本,这样模型不会为了迎合个别噪声点而剧烈扭曲。配合 RBF 核函数,SVR 在高维空间里拟合非线性关系,同时通过惩罚参数C控制模型复杂度。三个核心参数分别是:

参数作用取值建议
-s 3指定使用 epsilon-SVR 回归固定为 3
-t 2使用 RBF 径向基核函数固定为 2
-c惩罚系数,越大越容易过拟合1 到 100 之间搜索
-gRBF 核的 gamma 值,控制单个样本的影响半径0.001 到 1 之间搜索
-pepsilon 管道宽度,越大越平滑0.01 到 1 之间搜索

libsvm 在 Matlab 中的接口很直接:svmtrain负责训练,svmpredict负责预测。注意这个svmtrain不是 Matlab 自带的统计工具箱函数,而是 libsvm 编译产生的 mex 文件,两者函数名相同但行为不同,这一点后面排错部分会单独讲。

3.2 训练前必须先做归一化

SVR 对特征尺度敏感,RBF 核函数本身依赖样本间的欧氏距离,如果负荷值在几千到几万之间波动,而某个特征量级不同,距离计算会被大数值特征主导。常见做法是把训练特征和标签都缩放到 [0, 1] 区间:

% 特征归一化,使用训练集的 min/max [X_train_norm, ps_x] = mapminmax(X_train', 0, 1); X_train_norm = X_train_norm'; X_test_norm = mapminmax('apply', X_test', ps_x)'; % 标签归一化 [y_train_norm, ps_y] = mapminmax(y_train', 0, 1); y_train_norm = y_train_norm';

归一化时有一个典型错误:先用全部数据计算 min/max,再切分训练测试集。这会把测试集的统计信息泄露给训练过程,导致预测结果偏乐观。这里的做法是只用训练集计算映射参数ps_xps_y,测试集直接调用apply套用同一套映射。mapminmax默认按行处理,所以输入需要转置,输出后再转置回来。

归一化后的标签在预测阶段还需要还原:

% 预测并反归一化 [y_pred_norm, ~, ~] = svmpredict(y_train_norm, X_train_norm, model); % 训练集预测 [y_test_pred_norm, ~, ~] = svmpredict(y_test_norm, X_test_norm, model); % 测试集预测 % 反归一化到原始量纲 y_train_pred = mapminmax('reverse', y_train_pred_norm', ps_y)'; y_test_pred = mapminmax('reverse', y_test_pred_norm', ps_y)';

如果跳过反归一化,得到的预测值都在 [0, 1] 区间,无法和真实负荷直接对比,R² 算出来也会是错的。

3.3 svmtrain 的完整调用与参数解析

libsvm 在 Matlab 中训练模型的最小完整代码如下:

% 训练 SVR 模型 % -s 3 表示 epsilon-SVR % -t 2 表示 RBF 核 % -c 惩罚系数 % -g 核函数 gamma % -p epsilon 管道宽度 % -q 静默模式,不输出训练迭代信息 cmd = '-s 3 -t 2 -c 32 -g 0.125 -p 0.01 -q'; model = svmtrain(y_train_norm, X_train_norm, cmd);

svmtrain的参数顺序是固定的,先标签后特征,和 Matlab 自带fitrsvm的传参顺序不同。cmd是一个字符串,libsvm 用它解析参数,每个参数以空格分隔。-c 32表示惩罚系数为 32,-g 0.125是 gamma 值,这个组合在负荷预测里是比较常见的区间。如果不加-q,训练过程会输出迭代次数、支持向量数等信息,排错时建议去掉-q看完整输出。

3.4 用网格搜索找参数,而不是凭感觉

SVR 的参数组合对结果影响很大,手动试几个值很难得到理想效果。常见做法是网格搜索配合交叉验证,libsvm 自带的网格搜索脚本可以用,但这里给出一个更贴合 Matlab 流程的实现:

rng(42); % 固定随机种子,保证可复现 cList = [1, 10, 32, 100]; gList = [0.01, 0.05, 0.125, 0.5]; pList = [0.01, 0.1, 1]; bestMSE = inf; bestParams = ''; for c = cList for g = gList for p = pList cmd = sprintf('-s 3 -t 2 -c %f -g %f -p %f -q', c, g, p); % 5折交叉验证,libsvm 的 -v 参数返回分类准确率或回归MSE的均值 mse = svmtrain(y_train_norm, X_train_norm, [cmd ' -v 5']); if mse < bestMSE bestMSE = mse; bestParams = sprintf('-c %f -g %f -p %f', c, g, p); end end end end disp(['Best params: ', bestParams]); disp(['CV MSE: ', num2str(bestMSE)]);

-v 5表示 5 折交叉验证,此时svmtrain不再返回模型,而是返回交叉验证的平均 MSE。网格搜索的粒度决定模型上限,但搜索空间过大会让训练时间成倍增加。在实际项目中,我一般先粗搜确定数量级,再在小范围内细搜。上述代码把三组参数做笛卡尔积,最多训练 4×4×3=48 个模型,在负荷这类中等规模数据上完全可接受。

4. 完整训练流程与多指标评价

4.1 把整个流程串起来

在得到最优参数后,用全部训练数据重新训练最终模型,再对测试集预测。完整的流程可以组织成一个脚本:

%% 1. 读取数据 dataTable = readtable('data.csv'); loadSeries = dataTable{:, 2}; %% 2. 构造滞后特征 p = 24; % 小时级数据,用过去24小时预测下一小时 [X, y] = createLagMatrix(loadSeries, p); %% 3. 切分数据 numSamples = length(y); trainNum = floor(numSamples * 0.8); X_train = X(1:trainNum, :); y_train = y(1:trainNum); X_test = X(trainNum+1:end, :); y_test = y(trainNum+1:end); %% 4. 归一化 [X_train_norm, ps_x] = mapminmax(X_train', 0, 1); X_train_norm = X_train_norm'; X_test_norm = mapminmax('apply', X_test', ps_x)'; [y_train_norm, ps_y] = mapminmax(y_train', 0, 1); y_train_norm = y_train_norm'; %% 5. 最佳参数(由网格搜索得到) bestCmd = '-s 3 -t 2 -c 32 -g 0.125 -p 0.01 -q'; model = svmtrain(y_train_norm, X_train_norm, bestCmd); %% 6. 预测 [y_train_pred_norm, ~, ~] = svmpredict(y_train_norm, X_train_norm, model); [y_test_pred_norm, ~, ~] = svmpredict(y_test_norm, X_test_norm, model); %% 7. 反归一化 y_train_pred = mapminmax('reverse', y_train_pred_norm', ps_y)'; y_test_pred = mapminmax('reverse', y_test_pred_norm', ps_y)';

这份流程看起来平淡,但每一步都有隐藏细节。createLagMatrix函数要放在脚本同目录下,或直接写成子函数。mapminmax返回的ps_x是结构体,保存了映射的 min 和 max,测试集apply时不能重新计算。svmpredict的第一个参数虽然是实际标签,但在预测阶段真正使用的是第三个参数model,第一个参数只用来计算误差指标。

4.2 评价指标:不止 R² 一个

代码输出三个指标:R²、MAE、MBE。这三个指标各有偏重:

指标公式含义评价重点
1 - SS_res / SS_tot模型解释了真实方差的多少,越接近 1 越好
MAEmean(‖y_true - y_pred‖)平均绝对误差,反映误差的典型大小
MBEmean(y_pred - y_true)平均偏差,正值表示整体高估,负值表示整体低估

计算代码如下:

% 测试集指标计算 SS_res = sum((y_test - y_test_pred).^2); SS_tot = sum((y_test - mean(y_test)).^2); R2 = 1 - SS_res / SS_tot; MAE = mean(abs(y_test - y_test_pred)); MBE = mean(y_test_pred - y_test); fprintf('R2 = %.4f\n', R2); fprintf('MAE = %.4f\n', MAE); fprintf('MBE = %.4f\n', MBE);

R² 对异常值敏感,个别极端负荷点会让 SS_res 变大,R² 明显下降。MAE 更稳定。MBE 则是判断是否有系统性偏差的工具,如果 MBE 绝对值远大于 MAE,说明模型不只是随机误差,而是整体预测偏移,这时应该检查特征构造或归一化过程是否能反映负荷趋势。训练集的指标和测试集指标要同时看,如果训练集 R² 很高而测试集很低,多半是过拟合,需要降低-c或增大-p

4.3 画图看趋势:预测结果的可视化

数值指标之外,一定要画时序对比图,曲线重合度能直观暴露模型在峰值处的表现:

figure; plot(y_test, 'b-', 'LineWidth', 1.5); hold on; plot(y_test_pred, 'r--', 'LineWidth', 1.5); legend('真实值', '预测值', 'Location', 'best'); xlabel('样本点'); ylabel('负荷值 (MW)'); title('SVR 电力负荷预测结果对比'); grid on;

注意这里没有把测试集的时间索引完整保留,如果 data.csv 第一列是 datetime 类型,可以用dataTable{:, 1}截取对应测试区间的时刻作为 x 轴,图上时间语义更清晰。画图后重点观察峰值时刻,负荷预测的难点通常在早晚高峰,如果峰值处预测偏低,可以尝试提高-c让模型拟合更充分,同时适当缩小-p迫使回归线贴近样本。

5. 边界场景与常见误用

5.1 训练集 / 测试集划分的隐藏顺序问题

在 2.3 节中,切分数据前必须先构造滞后矩阵,再切分。如果先切分原始序列再分别构造滞后矩阵,会导致训练样本和测试样本的特征窗口相互重叠,测试集包含了训练集末尾的部分信息,验证结果虚高。这是一个在时间序列预测里特别容易踩的坑。

正确顺序是:完整序列 →createLagMatrix→ 切分Xy。滑动窗口在构造过程中跨越了潜在的训练测试边界,但由于样本是按时间排序的,测试集的第一个样本,其特征窗口完全位于训练集最后一段之后,不存在信息泄露。如果别人代码里把顺序写反了,指标再好看也不能作为模型真实性能的凭证。

5.2 libsvm 的 mex 文件与 Matlab 自带 svmtrain 冲突

项目里提供了svmtrain.mexw64svmpredict.mexw64,这两个文件的本质是编译好的 C 程序,Matlab 通过 mex 接口直接调用。如果你在代码里写svmtrain,有可能调用的是 Matlab 自带的统计和机器学习工具箱版svmtrain,而不是 libsvm 版。两者的参数格式完全不同:libsvm 用字符串cmd,自带工具箱用'KernelFunction', 'rbf'这样的键值对。而且在较新的 Matlab 版本中,自带版svmtrain已经被标记为不推荐,甚至可能被移除。

解决办法是确保 libsvm 文件夹在路径顶部,或者给脚本起始处加一句:

% 将 libsvm 所在文件夹置于路径最前 addpath('/your/path/to/libsvm/matlab');

如果要确认当前调用的到底是不是 libsvm 版本,可以用which svmtrain查看路径。如果显示的是...\toolbox\stats\...而不是...\libsvm\matlab\...,则需要addpath或重新编译。这个冲突在不同机器上表现不一致,而且错误信息往往是模糊的“参数数量不足”,排查时优先看which输出。

5.3 归一化时 NaN 和数据缺失处理

csv 数据在采集过程中可能因为采集设备故障或通信中断产生缺失值。如果直接让缺失值进入训练,svmtrain可能报错或结果异常。常见的处理方式有两种:

% 方法1:线性插值填补缺失 loadSeries = fillmissing(loadSeries, 'linear'); % 方法2:删除缺失时刻的样本(但需要同步删除特征窗口) nanIdx = isnan(loadSeries); loadSeries(nanIdx) = [];

方法一更适合负荷数据,因为负荷在相邻时刻高度相关,线性插值引入的偏差远小于删除样本带来的时序断裂。删除样本会导致时间间隔不连续,后续构造滞后矩阵时,跨过删除点的窗口实际上拼凑了不同日期的负荷,语义上不再连贯。如果缺失段太长,线性插值也不可靠,这时应该把这一段数据整段剔除,并且要意识到预测的连续性已经被打破。

5.4 多步预测 vs 单步预测的差异

这个模型是单变量单步预测——用历史 24 小时预测下一个小时。如果你想预测未来 24 小时,不能直接把模型跑 24 次迭代预测,因为每一轮预测的结果都会被当作真实值输入下一次预测,误差会逐步积累。常见做法是递归多步预测,但这会显著放大误差;另一种做法是多输出策略,一次预测未来 24 个值,但 libsvm 原生不直接支持多输出回归,需要通过多模型或向量回归技巧实现。如果只想做一个快速迭代版本,递归预测可以这么写:

% 递归预测未来 H 步 H = 24; currentWindow = X_test(1, :); % 初始特征取自测试集第一个样本 predSeq = zeros(H, 1); for step = 1:H currentNorm = mapminmax('apply', currentWindow', ps_x)'; predNorm = svmpredict(0, currentNorm, model); % 第一个参数任意占位 predStep = mapminmax('reverse', predNorm', ps_y)'; predSeq(step) = predStep; % 窗口滑动:去掉最旧的一个,加入最新预测值 currentWindow = [currentWindow(2:end), predStep]; end

这里每一步都把预测值拼到窗口末尾,再丢弃窗口最前面的值。递归预测的误差会随步长增加而增大,所以在评估时,通常只报告前几步的指标,或单独计算每一步的累积误差。这一段代码适合做对比实验,但不适合直接用于生产调度。

5.5 评估指标对比的标准口径

多指标评价时有一个容易被忽略的口径问题:R² 是基于全部测试样本计算的,但负荷序列有明显的日周期和季节性,如果测试集恰好包含一段平缓的夜间负荷,R² 可能会虚高。反过来,如果测试集包含多个高峰日,R² 会下降。更合理的方式是按天切分评估,计算每一天的平均指标,再取平均。这样得到的指标更能反映模型在不同场景下的稳定表现。下面是一个按天切分的简化版本:

% 假设按小时采样,每天24个点,测试集正好是整数天 numTestDays = floor(length(y_test) / 24); dailyR2 = zeros(numTestDays, 1); for d = 1:numTestDays idx = (d-1)*24+1 : d*24; ssRes = sum((y_test(idx) - y_test_pred(idx)).^2); ssTot = sum((y_test(idx) - mean(y_test(idx))).^2); dailyR2(d) = 1 - ssRes / ssTot; end meanDailyR2 = mean(dailyR2);

这段代码假设采样间隔为 1 小时,如果实际是 15 分钟,将 24 替换为 96 即可。按天计算指标后,还能观察哪几天模型表现差,通常是极端天气或节假日前后的负荷突变,这时仅靠单变量历史数据很难预测到位,需要引入外部变量或日历特征。

6. 进阶验证:残差分析与参数敏感度速查

6.1 残差序列自相关检验

SVR 预测的价值不仅在于误差小,还在于误差是否有规律。如果残差(真实值减预测值)表现出明显自相关,说明模型遗漏了某种时间结构,比如周期性成分没有进入特征窗口。一个快速检验残差自相关的方法是画自相关图:

residuals = y_test - y_test_pred; figure; autocorr(residuals, 40); title('预测残差自相关');

如果残差自相关在滞后 24 处显著超过置信区间,说明模型没有完全捕捉日周期性,这是特征构造的问题,而不是 SVR 的问题。这时可以尝试把历史窗口拉长到p = 48甚至p = 72,让模型看到更多历史信息。另一种做法是增加“当前时刻离上次高峰的间隔”作为特征,但这就属于特征工程范畴了。残差自相关检验的另一个作用是指出是否需要切换到带季节分量的模型,比如 SARIMA 或季节性分解后再预测。

6.2 参数敏感度速查表

共享一份快速参考表,用于在网格搜索前缩小范围,表中数据来自负荷预测任务中常见经验:

现象优先调整参数调整方向
训练集精度高,测试集精度差-c降低,比如从 100 降到 10
预测曲线过度平滑,峰值拍平-p降低,让误差管道收紧
训练集和测试集精度都差-g升高,比如从 0.01 升到 0.5
训练时间过长-c降低,同时考虑减小样本数
预测结果整体偏高或偏低检查 MBE可能是归一化边界效应或特征窗口末段值偏高

-g的影响比-c更敏感。gamma过小,RBF 核的作用半径很大,模型过于平滑;gamma过大,模型只看近距离样本,容易过拟合。在实际操作中,先用固定-p 0.01,对-c-g做二维搜索,确定大致区域后再加入-p的搜索,能显著节省时间。

6.3 结果复现的固定随机种子与版本注意项

libsvm 的训练过程本身不涉及随机数,但网格搜索如果使用交叉验证,svmtrain在内部可能对样本进行随机排序。虽然 libsvm 的交叉验证是按顺序切分的,不依赖随机种子,但为了复现结果,仍然推荐在脚本开头加rng(2024)。同时要注意 mex 文件依赖 Matlab 版本和操作系统位数,项目里提供的是mexw64,只能在 64 位 Windows 上的 Matlab 使用,如果换到 Linux 或 macOS,需要重新编译 libsvm。编译方法在 libsvm 的make.m里,Matlab 命令行直接运行:

% 在 libsvm/matlab 目录下执行 make;

编译前需要确保系统有可用的 C 编译器,Matlab 中运行mex -setup可以检查并选择编译器。如果使用的是 R2023b 及之后的版本,Matlab 默认使用 MinGW64 或 MSVC,make.m会自动挑选合适的配置。这个环节出错时,错误信息通常是找不到编译器或 SDK,先执行mex -setup解决环境问题,再重新运行make。编译成功后,svmtrain.mexw64会被重新生成,覆盖原文件。若不想覆盖原始文件,可以把编译产物放到另一个文件夹,然后通过addpath指向新位置。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询