1. 为什么光伏功率预测要从“点”走向“概率”在很长一段时间里我做光伏功率预测都习惯性只关注“明天中午12点能发多少电”。这个数值用起来很直观调度看一个数报告里写一个数软件界面大屏上也显示一个数。但真正和电力调度的人聊过之后才知道确定性预测在强波动天气下几乎没有决策价值——你说12点出50兆瓦结果实际只有20兆瓦那调度是按50兆瓦去安排备用容量还是按20兆瓦如果他按50兆瓦安排电网频率就是风险按20兆瓦安排又可能因为弃光或备用不足产生经济损失。所以做光伏功率预测不能只给一条期望曲线还应该告诉调度预测值有多大的不确定性值落在哪个区间内的概率是多少哪些时刻容易发生大幅波动。这就是概率预测。概率预测的落地形态大致有三种预测区间、分位数曲线、以及随机场景集合。要做这些单纯靠常用的神经网络或者支持向量回归很难直接给出联合概率分布尤其是当我们还希望捕捉“多个时间点之间”、“多个光伏电站之间”的时空相关性时问题就变得更复杂。我最后采用的方案是把单调广义学习系统MBLS和Copula理论组合起来MBLS负责做确定性预测并且通过单调性约束保证输出符合光伏功率的物理规律Copula则把这些输出在不同时刻、不同站点之间的相依结构连接起来形成完整的时空联合概率分布。这套方案的核心价值在于“拆分”先用一个灵活且易于训练的回归模型拿到边缘分布再单独用Copula去建模相关性结构。这样各模块可以分别调试也避免了直接用深度生成模型带来的训练不稳定问题。需要说明的是这整套流程是一套可复现的Matlab实现不需要深度学习框架也不需要GPU一台普通工作站就能跑完离线训练和在线预测。这篇文章适合谁看如果你正在做光伏功率预测、风电功率预测或者任何带时空依赖性的新能源出力预测如果你手里已经有了一批历史功率数据和气象数据却不知道怎样才能输出带置信区间的预测结果如果你听说过Copula但对怎么和机器学习模型结合没有直观概念——那这篇笔记应该能帮你少走不少弯路。我会尽量把模型设计、数据构造、Matlab代码、评估方法都讲清楚也会把踩过的坑一并列出来。2. 模型骨架MBLS与Copula如何协作2.1 广义学习系统的底子是随机特征映射加岭回归聊MBLS之前先得聊聊它的前身——广义学习系统Broad Learning SystemBLS。BLS是深度网络之外的一种增量式学习结构思路非常朴素先用一系列随机生成的权重把输入映射到高维特征空间然后再通过岭回归求输出权重矩阵。这个流程可以拆成三步。第一步对输入X做特征映射生成“特征节点”。常见的映射方式是对X做线性变换后加激活函数比如Z activation(X * W_e b_e)其中W_e是随机生成的映射矩阵。第二步把特征节点再经过一次随机投影和激活生成“增强节点”H activation(Z * W_h b_h)第三步把Z和H拼在一起形成一个广义特征矩阵A [Z | H]然后求解岭回归权重W_out (A^T A λ I)^{-1} A^T Y这里的λ是正则化系数I是单位矩阵。整个过程可以理解为把输入数据随机映射到高维空间再用线性回归拟合目标但因为特征维数很高所以拟合能力并不弱而且训练速度远快于传统BP神经网络。BLS最吸引人的地方是它不需要反向传播没有层层梯度传播带来的误差积累。特征节点和增强节点都是随机映射只有最后一层岭回归是可学习的所以训练过程非常稳定耗时也短。但它的缺点也很明显——随机映射没有任何约束输出不一定满足我们对物理世界的先验认识。2.2 引入单调约束让预测结果更符合物理常识光伏功率和辐照度之间理论上存在一个明显的单调递增关系辐照度越高功率通常越大。虽然由于组件温度、遮挡分布、逆变器限功率等因素实际数据中会出现局部波动但从整体物理规律来看预测模型的输出不应该出现“辐照度上升、预测功率下降”这种反直觉的现象。常规BLS不做任何单调性约束拟合出来的曲线在训练数据稀疏的区域经常会出现上下振荡。为了抑制这种振荡就有研究者提出了单调广义学习系统Monotonic Broad Learning SystemMBLS它通过在损失函数上增加惩罚项或者对网络权重施加非负/非正的符号限制保证模型输出关于某些指定输入变量单调。具体实现上我在代码里做的是“符号约束法”针对需要单调的输入维度要求第一层映射矩阵W_e中对应列的方向满足某种一致性从而让输出对该维度保持单调方向。比如我们希望功率关于辐照度单调递增那么从输入到特征节点的映射权重中辐照度对应的那一列在取激活函数之前就要保证与最终输出权重组合后的整体效果为正。这种方法的好处是不引入额外的对抗训练或复杂的约束优化只需要在岭回归求解时把约束转换成二次规划问题或者简单地在梯度更新后做投影。Matlab里有现成的优化工具箱可以搞定也可以自己写一个迭代投影算法。实际用下来MBLS比普通BLS在高辐照度突变区间更稳定预测曲线很少出现“锯齿状”波动。当然单调约束不是万能的。如果数据里本身存在“温度很高导致组件转换效率下降、功率反而降低”的效应而你把温度也列为单调递增变量那就违背物理事实了。所以实际建模时我只对“辐照度”这一个变量施加单调递增约束其他变量比如温度、湿度、风速都保持自由拟合让数据自己说话。2.3 Copula把“分时刻的边缘分布”缝成“时空联合分布”当我们用MBLS得到每个时刻的确定性预测y_hat(t)下一步就是描述不确定性。假设预测误差服从某种分布我们可以用历史残差来估计每个时刻的预测误差分布例如用核密度估计或者参数化分布拟合得到时刻t的累积分布函数F_t。但问题在于光伏功率在连续时刻之间是强相关的。中午12点的误差如果是因为云团遮挡造成的那么12点05分的误差大概率也偏大。如果我们简单地把每个时刻的误差分布独立抽样再生成一条“完整功率曲线”那这条曲线会非常毛糙一会儿高一会儿低完全不像真实的功率曲线。Copula的作用就是把这些边缘分布按照真实的相依结构“粘”在一起。Copula函数的本质是一个多元联合分布函数它的边缘分布都服从均匀分布U(0,1)因此可以用来连接任意边缘分布。描述相依结构时我们不关心每个变量本身是正态还是偏态只需要知道变量之间排序上的相关性——比如“12点功率高的时候12点05分功率往往也高”这种同步趋势。实际建模时我用的是阿基米德Copula族中的Gumbel Copula因为它可以刻画上尾相关性也就是在功率高值区间出现“抱团”现象。光伏功率在晴天时容易接近额定值此时多个时刻会同时维持在较高水平这种同步性用Gumbel Copula描述比较贴合数据特征。你也可以试试Clayton Copula下尾相关强或者Frank Copula对称相关具体哪个好要靠AIC/BIC去选。Copula的建模对象不一定是原始误差也可以是对概率积分变换后的残差序列。标准做法是把历史残差用经验累积分布函数做变换得到近似均匀分布的样本u_t然后对u_t拟合Copula参数。预测阶段再从Copula联合分布中抽样得到对应的u_1, u_2, ..., u_T再通过逆累计分布函数转成每个时刻的预测误差场景叠加MBLS的预测均值就得到一条完整的光伏功率时空场景。2.4 整体训练与预测流程把MBLS和Copula串起来整个流程分成离线和在线两个部分。离线阶段第一整理训练数据构造特征矩阵X和历史功率矩阵Y第二训练MBLS模型得到每个时刻的确定性预测第三计算训练集残差估计每个预测时刻的残差边缘分布比如用核密度估计KDE第四把残差转换为均匀分布样本拟合Copula参数。在线预测阶段第一输入当前的气象预报和时空特征用MBLS预测出未来的确定性功率曲线第二用训练好的Copula抽样生成M条时间相关性一致的误差场景第三将误差场景叠加到确定性预测上得到M条光伏功率场景第四对场景做统计输出各个时刻的预测区间、分位数曲线、期望功率等。这个过程非常模块化。如果你只想快速验证Copula部分可以不训练MBLS先用简单的线性回归替代确定性模型如果你只想看MBLS的效果也可以暂时不接Copula输出点预测后直接做残差分位数区间。模块化的好处是每个环节都能单独测试出了问题能快速定位。3. 数据准备时空特征的构造和边缘分布估计3.1 我们到底要让模型“看”什么数据光伏功率预测的输入特征一般分为几类历史功率数据、数值天气预报NWP数据、实测气象数据、以及时间标记时角、日期等。我在这个项目里用到的核心特征有历史功率序列当前时刻前1到3小时内的功率值步长15分钟作为时序记忆。辐照度预测序列来自数值天气预报的全球水平辐照度GHI预测值这是我们要施加单调约束的变量。温度预测环境温度会影响组件效率一般作为自由变量进入模型。时间编码包含一天中的小时、一年中的日期能帮助模型捕捉季节性规律。我会把它们编码成正弦/余弦两个分量避免“0点和23点”的距离被算法误判。除了当前预测时刻的特征之外时空预测模型还要考虑相邻预测时刻的特征。我构造特征的方式是“滑窗展开”如果要预测未来4个小时共16个时刻的功率我就把从预测起始时刻之前12个时刻以及之后若干时刻的气象特征全都拼在一起形成一个大的特征向量。这个做法看起来很笨但在BLS这种“高维特征线性输出”的结构里非常有效——它不需要像LSTM那样按时间步展开直接把上下文信息喂给模型即可。数据清洗方面需要注意三点剔除非发电时段比如夜间功率恒为0的记录如果保留会让模型学到“夜里功率为0”这种与实际光照无关的错误规律剔除逆变器限功率时段这类时段功率被削平真实辐照度远高于功率对应值极容易造成训练残差异常剔除数据跳变超过物理阈值的异常点。3.2 边缘分布估计别直接假设正态刚开始做概率预测的时候我图省事直接假设预测残差服从正态分布用均值加减1.96倍标准差构造95%置信区间。结果在晴天场景下看着还行一到多云天气区间明显偏窄因为真实残差分布是右偏的——少数时刻由于云团遮挡功率骤降产生很大的负残差。后来我改用核密度估计KDE来估计每个预测时刻的残差边缘分布。KDE的思想很简单对一组残差样本用平滑核函数在每个样本点上叠加一个“小山峰”得到连续的概率密度估计。Matlab自带的ksdensity函数可以直接调用带宽bandwidth建议用“Silverman”规则或者自己交叉验证选择。这里有一个容易被忽略的细节边缘分布应该是“每个预测时刻”或者“每组相似气象条件”单独估计而不是把所有时刻的残差混在一起估计一个总分布。因为正午时段的预测误差通常比清晨和傍晚更大不同季节也不同。我采用的做法是按照预测时刻比如未来第1个15分钟、第2个15分钟……分别收集历史残差然后对每个提前期估计一个边缘分布。这样下来一共有16个边缘分布如果预测16个时刻。刚开始觉得麻烦但评估CRPS和区间覆盖率时效果提升非常明显。3.3 站点/时间点之间的相关性到底怎么刻画“时空”二字包含两个维度时间上的连续相关以及空间上多个站点之间的相关。如果你的项目只针对单个光伏电站那重点就是时间维度如果你要做区域光伏群那空间维度就必须考虑——不同站点共享同一个云图和气象过程它们之间的功率不可能完全独立。我在处理单电站多时刻预测时先把16个预测时刻的残差变换成均匀分布样本u_1, ..., u_16然后用一个16维的Copula来刻画它们之间的相依结构。16维的维数并不算高Gumbel Copula或者高斯Copula都能处理。如果是多电站场景我倾向于“先降维再Copula”对M个电站、T个时刻的数据可以把M×T个变量看作一个高维随机向量但直接拟合所有变量的Copula很容易因为维数过高而失败。我的做法是先用主成分分析PCA或者因子分析把空间维度压缩到几个主成分再对主成分序列拟合一个较低维的Copula。抽样时生成主成分的随机场景再映射回原始变量空间。还有一种更简单但同样实用的路子用“高斯Copula 相关矩阵”来捕捉所有时空变量的相关性。高斯Copula的密度表达式简单参数就是相关矩阵Matlab里计算相关矩阵非常方便。缺点是它无法刻画尾部相关性但如果你主要关心的是区间预测而不是极端事件高斯Copula已经够用。4. Matlab实现的关键代码与参数选择4.1 MBLS主体函数签名与核心计算我自己写的MBLS训练函数长这样function model trainMBLS(X, Y, s, c, lambda, mono_idx) % X: 输入特征矩阵N行D列 % Y: 目标矩阵N行T列T为预测提前期数 % s: 特征节点个数 % c: 增强节点个数 % lambda: 岭回归正则化系数 % mono_idx: 需要单调递增的变量所在列索引第一步是生成随机映射权重。为了保证可复现我会用rng固定随机种子rng(42); We randn(size(X,2), s) / sqrt(size(X,2)); be rand(1, s) * 0.1; Z tanh(X * We be);然后生成增强节点Wh randn(s, c) / sqrt(s); bh rand(1, c) * 0.1; H tanh(Z * Wh bh); A [Z, H];如果你需要在某些维度上施加单调递增约束最简单的做法是构造A之后把正则化目标变成带线性不等式约束的二次规划。Matlab里有命令quadprog可以求解但写成公式比较繁琐。我在这里分享一个折中方案对mono_idx对应的权重进行梯度投影即在岭回归迭代过程中每次更新后把对应权重的负值清零。实现岭回归的解析解W_out (A * A lambda * eye(size(A,2))) \ (A * Y);如果特征维数特别高比如超过几万用解析解求逆会偏慢。这时可以改用最小二乘迭代法或者用随机特征映射配合批量梯度下降。但一般情况下BLS的增强节点数量控制在几百到几千解析解完全够用。4.2 Copula参数估计和场景抽样代码Matlab里对阿基米德Copula的支持不算直接但可以借助copulafit和copularnd函数来处理高斯Copula和t-Copula。如果你的数据适合Gumbel Copula那需要自己实现参数估计。我在这里给出一个高斯Copula的完整流程。先把残差变换成均匀分布样本。假设残差矩阵residual是N行T列U zeros(size(residual)); for t 1:T U(:,t) ksdensity(residual(:,t), residual(:,t), function, cdf); end因为ksdensity返回的是经验CDF所以U的每一列都近似服从U(0,1)。下一步拟合高斯Copula[Rho, nu] copulafit(Gaussian, U);如果要生成1000条场景可以用Usim copularnd(Gaussian, Rho, 1000);然后对Usim的每一列使用之前KDE估计的逆CDF把均匀分布样本逆变换回残差空间。Matlab里没有直接的icdf函数用于KDE所以需要事先存好CDF的插值点F_t ksdensity(residual(:,t), residual(:,t), function, cdf); edge_t residual(:,t); resid_sim_t interp1(F_t, edge_t, Usim(:,t), linear, extrap);最后把场景叠加到MBLS的确定性预测上scenarios repmat(Y_pred, 1, 1000) resid_sim; % 注意维度布局这里我踩过一个坑interp1在插值到区间边界外的时候可能会产生非常大的外推值因为U(0,1)的样本可能落在历史残差的最小值或最大值之外。解决办法是不要让样本完全均匀分布在0或1附近或者在ksdensity估计CDF时对边界做截断处理例如把CDF值限制在[0.001, 0.999]之间这样逆变换不会追求无穷远。4.3 训练中的超参数选择经验MBLS的超参数主要有三个特征节点数s、增强节点数c、岭回归正则化系数lambda。我在实践中的取值规律如下s决定特征映射的宽度一般取100~1000。如果输入特征维数是几十s取300左右就有不错的拟合能力如果特征维数上百可以适当增加到500以上。c决定非线性增强程度通常取s的0.5~1倍。太多会导致过拟合太少则欠拟合。lambda控制正则化强度我一般用交叉验证网格搜索候选值是10^{-3}到10^{-6}之间。BLS和深度网络不同它对lambda很敏感太大会把权重压得过小导致欠拟合太小又会放大噪声。实际测试中lambda1e-4到1e-5效果最好。另外特征映射中的激活函数我选用tanh比ReLU稳定因为ReLU在随机初始化的情况下容易出现大量神经元“死亡”导致特征矩阵稀疏。如果你发现训练误差很大先检查激活函数和随机权重尺度。权重初始化时除以sqrt(输入维度)这个细节很重要我一开始没注意直接用randn结果特征值相差几个数量级岭回归直接算出了NaN。5. 评估概率预测CRPS、区间覆盖和PIT5.1 定量指标CRPS与Pinball Loss点预测用RMSE和MAE评估概率预测则需要专门的指标。最常用的是连续排名概率分数CRPS和分位数损失Pinball Loss。CRPS可以理解为概率预测和真实观测之间距离的度量。如果真实值是y预测的累积分布函数是F那么CRPS定义为CRPS ∫ (F(z) - I(y ≤ z))² dz在Matlab中如果已经得到了预测场景我通常直接用场景集合构造经验分布函数然后数值积分计算CRPS。也可以用场景的排序来计算近似值将M个预测场景从小到大排序然后用离散近似公式求解CRPS。CRPS越低说明概率预测既准又锐利。Pinball Loss则是针对分位数预测设计的。假设我们的预测分位数是q_τ(t)实际观测是y(t)则损失为loss (y(t) - q_τ(t)) * τ, 当 y(t) q_τ(t) loss (q_τ(t) - y(t)) * (1 - τ), 当 y(t) q_τ(t)如果预测的是多个分位数例如5%、10%、……、95%把所有分位数的Pinball Loss平均一下就能得到综合评价指标。5.2 PIT直方图检验概率校准CRPS很低不一定代表预测分布校准良好。我见过有的模型CRPS指标还可以但预测区间严重偏向一侧导致保守但不校准。此时PITProbability Integral Transform直方图是一个很好的可视化诊断工具。PIT的计算方法是对于每个测试样本计算观测值在预测累计分布函数下的位置PIT F(y_obs)。如果概率预测是完美的校准那么PIT应该服从均匀分布U(0,1)。直方图如果呈现U形说明预测分布过窄极端值出现频率高于预期如果呈现倒U形说明预测分布过宽如果明显偏斜说明存在系统性偏差。我在评估时通常会画三张图PIT直方图、95%预测区间的时间序列图、CRPS随提前期的变化曲线。前两张帮助判断定性表现第三张用于和别的模型做定量对比。5.3 场景可视化怎么看生成多条场景曲线之后不要只把它当作“一堆折线”。我习惯把场景集合分成几个典型类别晴天场景、多云场景、波动场景。通过观察不同场景之间的相关性来判断Copula是否真正捕捉到了时间相依性。比如晴天场景下所有场景曲线应当像“一束紧扎的面条”互相之间非常接近多云场景下场景曲线会出现明显分化但每条曲线本身应该有相对平滑的走势而不是逐个时刻随机跳变。如果你发现场景曲线在时间上“抖动”得像白噪声那大概率是Copula没用好可能相依结构估计偏低或者你漏掉了时间维度上的相关性实际是把它当成了独立抽样。6. 我踩过的坑和排查思路6.1 Copula对边缘分布很敏感第一次跑通完整流程时我直接用Matlab的ecdf函数来估计CDF然后把残差变换成U样本。结果拟合出来的Copula相关矩阵过于“乐观”生成的场景几乎全挤在一起。后来才发现ecdf在样本量少的时候经验CDF是阶梯状的U样本在边界处会出现大量0或1严重干扰Copula参数估计。改用ksdensity的function,cdf选项之后U样本变得平滑相关矩阵也恢复正常。这件事给我的教训是边缘分布估计和Copula参数估计是一对“耦合环节”。边缘分布如果不够平滑后面全白搭。6.2 MBLS的单调性约束导致特征矩阵病态我在给辐照度增加单调约束时一开始用的是“权重投影法”——每次迭代后把对应权重中的负值强行置零。结果训练出的模型在训练集上误差很小但测试集上性能明显变差原因是强行置零破坏了原有特征空间的表示能力导致特征矩阵近似奇异。后来我改为在损失函数上增加单调性惩罚项即在原有岭回归目标里加入一项惩罚使得目标函数虽然鼓励拟合数据但也惩罚单调方向上的符号错误。这样有效缓解了病态问题。如果你也遇到这个问题可以考虑把惩罚系数从0.01开始调不要一上来就给太大权重。6.3 数据泄露归一化的时候把手伸到了未来这个错我犯得很典型。我在构造特征时把整段时间序列标准化用的是整个数据集的均值和标准差包括测试集。这样训练时模型“偷看”了未来数据的分布导致在线预测时性能大幅下降。正确的做法是只用训练时期的数据统计量来做标准化测试数据在评估时也用训练集的统计量进行变换。这一点在时间序列任务里必须刻在脑门上。6.4 多站点场景下的维度爆炸我把单电站方案直接扩展到10个站点、16个提前期总共160个变量尝试直接拟合高维Copula结果参数估计矩阵大得离谱而且优化极不稳定。后来用PCA降维后把主成分个数控制在10个再把主成分的场景映射回原空间效果好了很多。如果你也想用高维Copula可以考虑“藤结构”Vine Copula它能把高维相关结构拆成多个二维条件Copula灵活度更高。但Matlab里原生支持不多需要额外工具箱或者自己实现。对大多数光伏预测项目来说PCA高斯Copula已经能解决90%的问题。7. 个人体会与下一步扩展做这个项目最大的感受是“拆分问题”比“堆模型”更重要。光伏功率预测如果直接端到端生成概率场景训练和调参的复杂度会非常高而且很难解释。把问题拆成“确定性预测边缘分布Copula相关结构”三个独立环节后每个环节都可以用成熟的模块处理出现问题时定位也极快。如果你也想在自己的数据上复现这套方案我建议从单站点、16个提前期开始先用线性回归或分位数回归替代MBLS跑通整个Copula流程。然后再把确定性模型换成MBLS对比一下CRPS变化。这样即使MBLS部分调试不顺利Copula部分也已经验证过了不会把两个模块的问题混在一起。下一步我想尝试的方向有两个。第一个是使用时变的Copula参数让相关结构能随季节和天气类型变化另一个是把MBLS的单调约束扩展到多个变量比如对“辐照度”和“温度”分别施加不同方向的单调约束但需要更加精细的约束设计避免互相冲突。这些内容还在试验阶段等跑出可靠的结果再来分享。