做光伏功率预测的朋友应该都有过这种体会:白天的功率曲线像锯齿一样上下跳动,阴雨天又突然断崖式下跌,辐照度稍微有一点扰动,预测误差就被放大到没法看。传统单一模型在这种非平稳、强波动的序列面前,基本只能做到“大趋势对、细节全错”,一遇到天气突变就彻底失灵。我最近在MATLAB里完整实现了一套基于非线性二次分解的Ridge-RF-LSBoost组合预测模型,简单说就是先把原始光伏功率序列做两层分解,再用岭回归、随机森林和最小二乘提升三种模型分工预测,最后集成输出。整套流程跑下来,15分钟级光伏功率预测的效果比单一模型和普通分解模型稳定太多。
这篇内容适合三类人看:正在做光伏功率预测相关课题、需要写论文实验对比的研究生;做新能源功率预测工程落地、被精度指标反复折磨的工程师;以及第一次接触VMD、CEEMDAN这类信号分解方法,想在MATLAB里快速跑通一个完整案例的初学者。我会把从数据准备到模型集成的全流程拆开讲,包括参数为什么这么设、代码为什么这么写、哪些坑我踩过之后再也不碰,尽量做到你照着能复现。
1. 为什么需要“非线性二次分解+Ridge-RF-LSBoost”这种组合
1.1 光伏功率序列到底难在哪
光伏功率预测本质上是一个时间序列回归问题,但它比一般的销量预测、负荷预测更棘手。光伏功率受气象条件支配,而气象本身就是混沌系统,辐照度、温度、湿度、云量任何一个变量抖动,都会直接反映在功率曲线上。更麻烦的是,光伏功率序列有非常明显的多尺度特征:分钟级的高频脉动、小时级的云层过境波动、天级的日出日落周期、季节级的太阳高度角变化。
这些不同尺度的成分叠在一起,让原始功率序列呈现出强非平稳性、强波动性和复杂的非线性能量分布。你直接拿原始序列丢给随机森林或者LSTM,模型会非常痛苦——它既要学到缓慢变化的趋势,又要捕捉短促的突变,还要区分哪些波动是噪声、哪些是真实气象信号。结果往往是顾此失彼。
我见过不少人在这个问题上绕弯路:有人先对原始序列做一次EEMD分解,再用单一的LSTM逐分量建模,效果比直接建模好一些,但一到天气剧烈变化的日子,误差还是明显偏高。原因很简单:一次分解只能把序列划分到某个粒度,最高频的那个分量里仍然混杂着大量不规则扰动,后面模型照样要从“混沌”里学规律。
1.2 模型组合背后的分工逻辑
这套方案的核心思路是“分而治之,异构互补”。先用两层分解把复杂序列拆成多个相对平稳的子序列,然后针对不同子序列的特征,用不同模型分别预测,最后叠加输出。你可以把它类比成一个公司团队:有人擅长做长期战略规划,有人适合处理突发事件,有人专门扣细节,让一个人把所有活干完,效果一定不如专业分工。
三个模型的选择不是随手拍的,每一个都有明确的分工:
- Ridge(岭回归)本质是带L2正则化的线性回归,它的优势是稳定、抗共线性、不会过拟合,适合预测分解后那些变化平缓、趋势主导的低频分量,比如日周期趋势项和残差项。
- RF(随机森林)是Bagging类集成模型,通过多棵决策树并行投票来拟合非线性关系,对输入特征不敏感、对异常值有天然抵抗,适合处理中等频段、有一定非线性的分量。
- LSBoost(最小二乘提升)是Boosting类模型,每棵树都在拟合前一棵树的残差,擅长把前两个模型漏掉的细节一点点找补回来,适合处理最高频的波动分量。
这三个模型放在一起,正好覆盖了“线性趋势—非线性过渡—高频细节”三个层次。后面我会详细讲怎么给它们分配分量、怎么调参,以及集成时有哪些容易翻车的细节。
2. 核心细节拆解:二次分解怎么做、模型参数怎么定
2.1 第一层分解:VMD参数选择的现场经验
第一层分解我推荐使用VMD(变分模态分解),它在频带划分上比EMD/EEMD更干净,没有模态混叠那么严重,而且数学基础扎实,论文里也好解释。VMD的调用方式很简单,关键在参数设置。
% 第一层分解:VMD % 输入P为原始光伏功率序列(列向量) K = 6; % 模态数量 alpha = 2000; % 惩罚因子 tau = 0; % 噪声容忍度,设为0表示严格保真 DC = 0; % 第一个模态是否为直流分量 init = 1; % 中心频率初始化方式 tol = 1e-7; % 收敛容差 [u, ~, omega] = VMD(P, alpha, tau, K, DC, init, tol); % u是K行N列的矩阵,每一行是一个模态分量 % omega是每个模态的中心频率K是最要命的参数。K太小会欠分解,低频趋势和高频细节还绞在一起;K太大会过分解,把一个完整物理过程硬切成好几段,导致模态混叠和虚假分量。我试过一组数据上K从3取到10,K=6的时候各模态的中心频率分布最均匀,没有两个模态的中心频率贴得过近。通行的做法是观察VMD输出的中心频率序列:如果相邻两个模态的中心频率相差太小,说明K取大了;如果最高频模态的中心频率仍然很高且波动剧烈,说明K取小了,还得继续拆。
alpha惩罚因子控制的是模态带宽的紧凑程度。alpha越大,每个模态的频带越窄,分得越干净,但过大会把真实信号的高频部分切掉。我一般先试2000,再对比200和2000的结果,选重构误差最小的那组参数。重构误差就是sum(abs(sum(u,1) - P')),这个值越小,说明分解保真度越高。
2.2 第二层分解:为什么只分解“最复杂”的那一个分量
很多人做“二次分解”时有个误区:对全部模态再进行一次分解,最后得到十几个、二十几个分量,训练时间爆炸,而且很多子分量毫无物理意义,纯粹是把噪声拆成了更碎的噪声。
这套方案里的二次分解只针对“最复杂的那一个模态”。怎么定义最复杂?用样本熵(Sample Entropy),它衡量的是时间序列里新模式出现的概率,数值越大说明序列越复杂、越不可预测。低频趋势项的样本熵很低,高频波动项的样本熵很高。我们把VMD分解出的K个模态分别计算样本熵,挑出样本熵最大的那个模态,再做一次CEEMDAN分解。
% 计算样本熵(简化版,实际要写完整循环) % 数据长度N,嵌入维度m=2,容差r=0.1*std(seq) function SampEn = SampEn(seq, m, r) N = length(seq); % 构造模板向量 for i = 1 : N - m templateVec(i,:) = seq(i : i + m - 1); end % 统计匹配数量并计算条件概率比的对数 B = 0; A = 0; for i = 1 : N - m for j = i + 1 : N - m if max(abs(templateVec(i,:) - templateVec(j,:))) <= r B = B + 1; end end end for i = 1 : N - m - 1 for j = i + 1 : N - m - 1 if max(abs(seq(i:i+m) - seq(j:j+m))) <= r A = A + 1; end end end SampEn = -log(A / B); end % 计算每个模态的样本熵 for k = 1 : K SampEn_values(k) = SampEn(u(k,:), 2, 0.1 * std(u(k,:))); end [~, idx_complex] = max(SampEn_values); complex_component = u(idx_complex, :);之后对这个最复杂模态做CEEMDAN二次分解。CEEMDAN解决了EEMD的残留噪声问题,每次迭代都会自适应调整加噪幅值,分解结果更稳定,适合处理已经被VMD“切薄”了的高频分量。
% 第二层分解:CEEMDAN % 常用参数:Nstd为加噪标准差比例,NR为实现次数 Nstd = 0.2; NR = 500; MaxIter = 5000; [sub_modes, ~] = ceemdan(complex_component, Nstd, NR, MaxIter);有人会问,为什么CEEMDAN不做在原始序列上,而要做在VMD的某个模态上?因为VMD负责频带划分,CEEMDAN负责把复杂频带里的非线性细节再分离一层,两者是互补的。这个“先分频带、后拆复杂分量”的两级结构,就是标题里“非线性二次分解”的落地含义。经过两层分解后,总的子分量数量是K-1加上二次分解的子模态数量,每个子分量的平稳性都大幅提升,后面建模型就轻松多了。
2.3 三种模型的核心参数与适用性
给每个子分量分配模型时,我的原则是看样本熵:样本熵低的低频趋势分量给Ridge,样本熵中等的中频分量给RF,样本熵最高的高频子分量给LSBoost。这个规则不是死的,但实际跑下来命中率很高。
Ridge的核心参数是lambda正则化系数。lambda越大,模型权重被压得越狠,模型越简单。低频分量本身规律性强,lambda不用太大,我常用1e-3到1之间的网格搜索,配合交叉验证选最优。
% Ridge回归 lambda = 0.01; % 网格搜索确定 B_ridge = ridge(y_train_low, X_train_low, lambda, 0); y_pred_low = [ones(size(X_test_low,1),1), X_test_low] * B_ridge;注意MATLAB的ridge函数不会自动加截距列,返回的B是第一列截距、后面是各特征系数,预测时得手动拼一列1,这个顺序搞错会导致结果全乱。
RF的核心参数是树的数量和最小叶子节点数。树太少拟合不足,树太多训练时间线性增长但精度提升有限,通常100到300棵足够。最小叶子节点控制树的复杂度,设得大能防过拟合,但太小又欠拟合。对于中等频段分量,我常用MinLeafSize=5,效果比较平衡。
% 随机森林 numTrees = 200; Mdl_rf = TreeBagger(numTrees, X_train_mid, y_train_mid, ... 'Method', 'regression', 'MinLeafSize', 5); y_pred_mid = cellfun(@str2double, Mdl_rf.predict(X_test_mid));TreeBagger预测输出是cell数组,必须转成数值,这个我在老版本的MATLAB上被坑过好几次,新版输出格式有所变化,最好whos一下确认类型。
LSBoost的核心参数是学习率和迭代次数。学习率决定每一步走多大,调小能提高精度但需要更多轮数;迭代次数太多会过拟合。我通常配LearnRate=0.1、NumLearningCycles=200起步,再根据验证集曲线调整。还有一个参数容易被忽略:MaxNumSplits控制每棵树的分裂次数,LSBoost默认用浅树(桩树),如果分量结构不那么复杂,可以用默认值。
% 最小二乘提升 Mdl_lsboost = fitrensemble(X_train_high, y_train_high, ... 'Method', 'LSBoost', ... 'NumLearningCycles', 200, ... 'LearnRate', 0.1); y_pred_high = predict(Mdl_lsboost, X_test_high);表格对比一下三种模型的核心特性,方便后续按需选取:
| 模型 | 核心机制 | 擅长分量 | 关键参数 | 典型陷阱 |
|---|---|---|---|---|
| Ridge | 带L2正则的线性回归 | 低频趋势、平稳残差 | lambda | 忘记手动拼截距列 |
| RF | Bagging决策树集成 | 中频非线性分量 | 树数、最小叶子 | 预测输出为cell类型 |
| LSBoost | Boosting残差拟合 | 高频复杂分量 | 学习率、迭代轮数 | 学习率过大导致过拟合 |
3. MATLAB完整实操流程与代码实现
3.1 环境与数据准备
这一套代码需要MATLAB R2019b以上版本,最好安装统计和机器学习工具箱(Statistics and Machine Learning Toolbox)。VMD代码需要从网上下载VMD工具箱,CEEMDAN也有公开的MATLAB实现,这两个工具箱网上很好找,下载后把对应文件夹加入路径即可。
数据层面,我这边用的是某公开光伏数据集,采样间隔15分钟,时间跨度一年,包含有功功率(单位kW)和气象要素时序数据(辐照度、温度、湿度、风速)。拿到数据后第一件事是画一张原始功率曲线,眼睛扫一遍:看有没有超过装机容量的异常值、有没有夜间非零值、有没有连续缺失段。异常值用filloutliers替换或者直接剔除,夜间非零值要结合辐照度阈值判断是真实阴天杂散辐射还是传感器噪声。数据质量不过关的话,后面模型再花哨也是垃圾进垃圾出。
% 数据读取 data = readtable('pv_data.csv'); P_raw = data.Power; % 原始功率 G = data.Irradiance; % 辐照度 T = data.Temperature; % 温度 H = data.Humidity; % 湿度 % 异常值处理 P_raw = filloutliers(P_raw, 'linear', 'movmedian', 96);3.2 特征构建与数据划分
分解完成得到子分量后,需要对每个子分量构建特征矩阵。光伏功率预测最自然的特征就是历史功率值(滞后项)和当前气象要素。滞后项的数量一般取采样点的24倍对应6小时(15分钟一个点),但也可以根据自相关函数截尾位置确定。滞后项太多会让特征维度爆炸,太少又吃不到周期信息,我常用24到48个滞后点。
每个子分量的特征矩阵这样拼:当前时刻的辐照度、温度、湿度、风速,加上该子分量的前24个滞后值。
% 构建特征矩阵(以某个子分量为示意) numLag = 24; X = []; y = []; for t = numLag + 1 : length(component) features = [G(t), T(t), H(t), component(t-numLag : t-1)]; X = [X; features]; y = [y; component(t)]; end数据划分必须按时间顺序切分,绝对不能随机打乱。时间序列一旦打乱,训练集和测试集就混进了未来信息,模型精度虚高,工程上一部署就露馅。我通常前70%训练、后30%测试,如果数据量充足,可以再用其中一部分做滚动验证。这里有个常见错误是归一化:必须先只统计训练集的均值和标准差,再用这套统计量去归一化测试集,禁止在全数据集上一次性归一化。否则测试集信息通过归一化参数泄漏到训练过程,结果严重失真。
% 用训练集统计量归一化测试集,防止未来信息泄漏 [Xn_train, PS] = mapminmax(X_train', 0, 1); Xn_test = mapminmax('apply', X_test', PS); Xn_train = Xn_train'; Xn_test = Xn_test';mapminmax默认按行处理,所以要转置成“特征×样本”格式,这行代码的顺序最容易被忽略。PS结构体里保存着训练集的min和max,测试集通过apply复用同一套映射,这就是正确做法的核心。
3.3 全流程集成的代码骨架
把前面几部分串起来,整个流程的骨架是这样:
% ===================== 完整流程骨架 ===================== % 步骤1:第一层分解 VMD [u, ~, omega] = VMD(P_raw, alpha, tau, K, DC, init, tol); % 步骤2:计算样本熵,定位最复杂模态 SampEn_values = getSampEnAll(u); [~, idx_complex] = max(SampEn_values); % 步骤3:第二层分解 CEEMDAN,作用于最复杂模态 [sub_modes, ~] = ceemdan(u(idx_complex, :), Nstd, NR, MaxIter); % 步骤4:合并所有子分量 all_modes = [u(1:idx_complex-1, :); u(idx_complex+1:end, :); sub_modes]; % 步骤5:对每个子分量构建特征矩阵并训练对应模型 for i = 1 : size(all_modes, 1) [X_i, y_i] = buildFeatures(all_modes(i, :), G, T, H, numLag); [X_train_i, y_train_i, X_test_i, y_test_i] = splitSequence(X_i, y_i, ratio); % 根据样本熵选择模型 if sampEn_i < th1 pred_i = trainRidge(X_train_i, y_train_i, X_test_i); elseif sampEn_i < th2 pred_i = trainRF(X_train_i, y_train_i, X_test_i); else pred_i = trainLSBoost(X_train_i, y_train_i, X_test_i); end % 拼接所有分量预测 pred_components(i, :) = pred_i; end % 步骤6:叠加输出最终预测 final_pred = sum(pred_components, 1);集成预测时,大家最容易忽略“分量预测对齐”的问题。每个子分量的特征构建是从numLag+1开始的,所以预测序列天然比原始序列短numLag个点。在全流程合成时,要么先预测再加回延迟段,要么统一计算并裁剪到公共长度。我习惯先在buildFeatures里就把每个子分量都延展一个NaN段,保证所有分量长度与原始数据一致,最后直接相加不报错。
3.4 误差评价与结果可视化
评价指标是论文和工程验收都得过的一道关。我通常输出四个:RMSE(均方根误差)、MAE(平均绝对误差)、MAPE(平均绝对百分比误差)、R2(决定系数)。RMSE对大误差敏感,能反映“最差的时候有多差”;MAPE对量纲敏感,光伏功率在夜间接近零值时,MAPE会瞬间爆炸,所以计算时要用辐照度阈值过滤掉夜间点,或者对分母加一个min保护。
% 评价指标计算 RMSE = sqrt(mean((y_true - y_pred).^2)); MAE = mean(abs(y_true - y_pred)); MAPE = mean(abs((y_true - y_pred) ./ max(y_true, 1e-3))) * 100; SS_res = sum((y_true - y_pred).^2); SS_tot = sum((y_true - mean(y_true)).^2); R2 = 1 - SS_res / SS_tot;可视化建议画“原始值-预测值”对比曲线加误差柱状图,如果验证集含有一段晴天和一段阴雨天,对比效果会非常直观。我经常看到有人只报指标不画图,在论文评审或工程汇报时特别吃亏——图一放出来,模型哪里好哪里坏一眼就清楚了。阴天段的误差通常是晴天的两到三倍,这个现象本身就可以作为后续改进的切入点。
4. 常见问题排查与避坑实录
这一节把我在运行和调试这套模型时遇到的高频问题整理成一个速查表,每一个都是真金白银踩出来的。
| 现象 | 可能原因 | 排查与解决 |
|---|---|---|
| 预测值偶尔出现负数 | 分解分量叠加后出现负值,或模型外推越界 | 对最终功率做max(0, pred),并裁剪到装机容量上限 |
| K取大后各模态中心频率相互靠拢 | VMD过分解,模态混叠 | 降低K,观察中心频率相邻差值应呈明显梯度 |
| CEEMDAN分解后子模态仍有噪声 | NR实现次数不足 | 将NR提高到500以上,或增大归一化标准差参数 |
| Ridge预测低频项总是偏低 | lambda过大,模型被压得太简单 | 网格搜索lambda,观察训练/验证误差拐点 |
| RF训练很慢但精度不增 | 树数破千,收益递减 | 先固定300棵树,用MaxNumSplits控深度 |
| LSBoost过拟合明显 | 学习率太大、迭代轮数过多或树过深 | LearnRate降到0.05,轮数减半,限制每棵树深度 |
| 各项指标训练集极高、测试集崩掉 | 特征包含未来信息或归一化泄漏 | 检查是否全局归一化;检查滞后特征是否误用未来时刻值 |
| 夜间MAPE异常高 | 功率接近0,分母太小 | 用辐照度>50W/m^2过滤白天点再算MAPE |
还有一个坑是MATLAB的fitrensemble在LSBoost下默认用树桩,也就是每棵树只能分裂一次。如果分量结构稍微复杂一点,树桩就能拟合得很好;但如果验证损失一直降不下来但也不过拟合,可以去调整MaxNumSplits,我一般试4到8的范围。
关于负功率的修正,我多说一句。光伏功率预测和别的回归不同,输出有明确的物理边界:白天最大不超过装机容量,夜间无辐照度时就应该接近零。模型预测输出落在边界之外时,直接置零或裁剪是合理的,不是“作弊”,而是把先验物理知识注入后处理。这个步骤很多论文里不写,工程上却非常有用。
5. 后续扩展方向与个人经验总结
这套“两层分解+异构模型集成”的思路,并不只适用于光伏功率预测。我后来把它迁移到短期风电功率预测和电力负荷预测上,把输入信号换成对应序列,特征矩阵做相应调整,整体框架几乎原封不动就能复用。风电功率的波动比光伏更随机,分解层数可能需要调大一些,但架构性的东西是一致的。
如果你想继续往深度做,有三个方向值得考虑。
第一是天气类型聚类。光伏功率在不同天气类型下呈现出完全不同的统计特征——晴天几乎是一条平滑钟形曲线,多云天是随机锯齿状,雨天是低位平坦波动。直接把所有天气混在一起建模,等于强迫一个模型同时学会三种完全不同的场景。可以先用K-means或DTW聚类把历史样本按天气类型分组,每组分别训练一套分解参数和模型参数,预测时先判定天气类型再选择对应模型,精度还能再上一个台阶。
第二是辐照度预报误差的兼容。代码里把辐照度作为特征输入,但在真正的短期预测中,未来时刻的辐照度只能来自数值天气预报或卫星反演,预报复盘是“以确定的输入预测不确定的未来”和“以不确定的输入预测不确定的未来”,这两件事的难度完全不同。在做实验对比时,要明确区分你是用了实测辐照度还是预报辐照度,否则工程上会翻车。
第三是分量预测策略的进化。这里的Ridge、RF、LSBoost可以被替换成GRU、Transformer或者轻量级深度学习模型,尤其是那些样本熵仍然偏高的高频子分量,序列结构复杂但长度较短,小规模GRU往往能比LSBoost吃得更透。
我个人在实际操作中最深的体会是:这套模型的创新点不在某个单一模型,而在于“分层分解+异构模型专业分工”的架构。你没必要跟风堆叠多个大模型,先把分解和分工做好了,三个基础模型也能跑出很能打的结果。再往后,你可以逐步把其中表现最差的那个分量换成更先进的模型,这种迭代方式风险低、可控性强,每一步的精度提升都能落到指标上。
最后再分享一个小技巧:每次调参后,把完整结果存成结构体,包括分解参数、模型参数、预测结果、评价指标和画图数据。连续对比时直接在结构体层面操作,不用每次重新跑全流程。这个习惯省了我大量时间,也让我在复盘模型效果时能立刻定位是哪一层分解、哪一个模型在拖后腿。