★ FEATURED ARTICLE 以前我在做风功率预测项目的时候最头疼的其实不是预测算法本身而是预测做完之后没法交代“你这个预测到底准在哪”。RMSE、MAE、MAPE这些指标我当然会给但调度那边追问的是另一句那你告诉我明天这个风电场群的出力到底可能偏多少偏大还是偏小多场站同时偏的概率有多大这时候只给一个MAPE是远远不够的必须把“预测误差”本身建出一个模型来而且要带上时间和空间上的相关性。这篇东西就是围绕这个主题写的用Matlab实现考虑时空相关性的风电功率预测误差建模与分析。适合正在做风功率预测、储能容量配置、备用容量评估、或电力市场竞价策略的研究生和工程师尤其是手上已经有一批历史预测值和实测值、但不知道下一步怎么把误差用好的人。我会从误差模型的用途讲起到Copula和AR模型的数学结构再到Matlab完整实现和算例对比最后把我在实际代码里踩过的坑都列出来。核心思路是不要把每个风电场站的预测误差当成独立同分布的白噪声时间上有持续偏差空间上有同步放大的效应这两件事必须同时建模才靠谱。1. 误差建模的用途决定你到底该建“什么模”1.1 先搞清楚预测误差建模和预测精度评价不是一回事很多人在项目里只算完RMSE就收工了。RMSE是一个标量它告诉你“平均大概差多少”但完全无法回答“偏大的概率是多少”“误差的分布是不是对称的”“多个场站之间会不会同时偏”。预测精度评价是对过去预测效果的一个汇总而预测误差建模是要用一种数学语言描述未来误差的完整概率分布并且能够从中采样生成未来可能出现的误差场景。这两者的差别举个例子就很直观。A风电场和B风电场单站RMSE都是25MW看起来预测水平差不多。但A的误差接近高斯分布很少出现超过两倍RMSE的大偏差B的误差是重尾分布三天两头出现三倍以上RMSE的极端偏差。对调度来说B的风险比A要大得多但RMSE完全看不出来。这就是为什么必须建误差模型而且要建到分布层面。1.2 哪些业务场景真正需要误差模型误差模型在风电并网里至少有四个直接用途我按实际项目里见到多的程度排序。一是备用容量配置。调度需要为预测偏差预留旋转备用如果误差模型说95%置信水平下系统总误差不超过80MW那备用就可以按这个数字去安排而不是拍脑袋乘个固定系数。二是置信区间预测。风功率预测系统不能只出一条曲线要给出上下边界误差模型的好坏直接决定区间是否够窄又覆盖够准。三是场景生成用于随机优化。比如做储能充放电策略或者机组组合时需要生成很多条可能的功率曲线传统做法是对预测值叠加独立高斯噪声这会把场景之间的相关结构完全忽略掉。四是电力市场偏差考核。市场规则通常对实际出力与日前预测的偏差进行惩罚有误差模型才能算清楚申报策略下的期望惩罚成本。1.3 单场站与多场站误差模型的复杂度完全不同单个场站的误差建模相对简单主要是时间维度的问题前一个时段的偏差会不会延续到下一个时段。但到了多场站事情立刻变复杂。同一片区域里的风电场往往受同一个天气系统控制冷锋过境时大家一起爬坡静稳天气时大家一起低出力预测误差也会出现明显的同步性。如果把每个场站的误差独立处理那么系统总误差的方差会被严重低估。这个问题在实际中特别隐蔽单站误差模型每个都拟合得很好加总之后却跟实测总误差对不上。所以真正有价值的误差模型必须把“每个场站自己的边缘分布”和“场站之间的相关结构”分开建模。这正好是Copula理论的用武之地后面我会详细讲怎么在Matlab里落地。2. 时空相关性的数学表达边缘分布 相关结构2.1 时间相关功率误差的“记忆效应”与AR阶数判断风电功率预测误差并不是白噪声它存在显著的“记忆效应”。原因不复杂数值天气预报NWP的初始场误差和边界条件误差往往持续好几个小时统计修正模型如果输入特征没变输出偏差也不会立刻消失。反映在数据上就是误差序列的自相关函数在滞后1到6小时仍然有明显的正值。我的经验是先对归一化误差算ACF自相关函数和PACF偏自相关函数来判断时间依赖的阶数。绝大多数风电场数据的误差序列用AR(1)就能抓住主要动态个别天气过程变化慢的场站可能需要AR(2)。AR(1)的含义很直观e_t φ · e_{t-1} ε_tφ是持续性系数通常在0.3到0.7之间。φ越大说明偏差越“顽固”今天上午偏大的风电场下午大概率还偏大。残差ε_t才是真正“新进来”的随机成分。做时空建模的时候不应该直接对原始误差序列建模而是要先把时间动态过滤掉然后对残差ε_t做空间相关建模。这样时间相关由AR刻画空间相关由残差的Copula刻画两者不混淆。2.2 空间相关同一天气过程下的误差同步放大空间相关是另一个维度的效应。距离50公里以内的两个风电场风速预测误差可能有0.6以上的相关系数距离拉到200公里相关系数通常降到0.2以下。这是因为天气系统的空间尺度决定了预报误差波及的范围。举个例子某次冷锋过程预报位置偏东50公里可能导致区域内的风电场同时低估风速、同时低估出力。单站看每个场站误差都在正常范围内但加在一起整个区域的总误差出现了系统性偏差。这种同步性如果不在模型里体现后果很严重按独立假设算出来的系统总误差标准差可能只有实际值的70%左右备用容量会算少了。调度最怕的就是这种“把风险算小了”的错误。2.3 Copula如何把相关结构从边缘分布里剥离出来如果每个场站的误差边缘分布都是高斯分布那么直接用多元高斯分布建模就行。可惜实测风电预测误差经常有偏斜、重尾、甚至双峰现象硬套高斯会低估极端偏差的概率。这时候用Copula是个比较优雅的办法。Copula的核心思想是把多维随机变量的联合分布拆成两部分一部分是每个变量的边缘分布另一部分是变量之间的相关结构。相关结构由一个Copula函数刻画边缘分布可以各自取不同的形状。Matlab统计工具箱里提供了现成的Gaussian Copula和t Copula函数不需要自己写复杂的数学公式。简单说建模分三步对每个场站的误差序列分别拟合一个边缘分布F_i可以是核密度估计、Beta分布或者带偏斜的分布把误差值代入F_i得到均匀分布的分位数u_i F_i(e_i)对u_i做标准正态逆变换z_i Φ⁻¹(u_i)然后对z向量估计协方差矩阵Σ这就是Gaussian Copula的核心参数。后边生成场景的时候是反着操作先从N(0,Σ)采样得到z再通过uΦ(z)变回均匀分布最后用每个场站自己的边缘分布逆函数F_i⁻¹(u)变回误差值。整个过程在Matlab里代码也就十几行。3. 用Matlab实现误差序列提取与相关性量化3.1 数据清洗与误差序列构造建模型第一步是把预测误差序列构造出来。我通常把数据整理成timetable格式包含三到四列时间戳、装机容量、预测功率、实际功率。误差定义为err_t (P_forecast,t - P_actual,t) / P_capacity为什么要用装机容量归一化因为不同场站规模不一样直接用MW值会让大场站的误差在空间协方差矩阵里天然占主导掩盖相关结构的真实形状。归一化之后各场站误差都落在[-1,1]量级既便于拟合Beta分布也便于比较空间相关性。数据清洗这一步最容易偷懒出错。至少要做三件事。一是剔除异常停机时段风机停机检修时预测和实际都存在很长的零值这种时段对误差分布没有参考意义。二是剔除通讯中断或数据跳变超过阈值的点比如实际功率在15分钟内从20MW跳到80MW再跳回20MW明显是数据质量问题。三是按时间戳严格对齐不同场站的数据采集延迟可能造成5分钟到15分钟的错位空间相关性会被削弱。% 假设data是已经对齐的timetable包含列: % forecast_W, actual_W, capacity_W data.err (data.forecast_W - data.actual_W) ./ data.capacity_W; % 剔除异常样本误差绝对值超过1.5的情况基本是数据错误 data data(abs(data.err) 1.5, :); % 剔除停机时段实际功率极低且预测功率极低 idxKeep ~(data.actual_W 0.05*data.capacity_W ... data.forecast_W 0.05*data.capacity_W); data data(idxKeep, :);3.2 用ACF/PACF判断时间相关阶数对每条误差序列算ACF和PACF直接看图形和数值。Matlab内置的autocorr和parcorr函数就够了。figure; autocorr(data.err, 48); % 假设是15分钟数据看12小时的依赖 title(15min误差自相关);如果ACF在滞后1到4阶明显超出置信区间说明时间相关存在。再看PACF如果在滞后1阶之后就截断AR(1)就够如果滞后2阶还显著可以考虑AR(2)。我自己的观察是15分钟分辨率的数据往往ACF衰减相对快因为风电功率本身的波动已经把部分高频特性带走了小时级数据的ACF更持久持续性系数常年稳定在0.5以上。3.3 空间相关矩阵估计与显著性检验时间相关过滤完之后对残差序列算空间相关。这里有个小陷阱直接对原始误差算Pearson相关会混入时间相关导致的“虚高”相关性。比如两个场站如果都存在早晚偏差变大、中午偏差变小的共同日内模式那么即使它们之间没有任何物理联系误差序列也会表现出正相关。所以正确的做法是把AR过滤后的残差拿来做空间相关分析。残差提取代码% 对每个场站分别做AR(1) for i 1:n mdl arima(1,0,0); est estimate(mdl, errAll(:,i), Display, off); epsAll(:,i) infer(est, errAll(:,i)); % 残差 end % 残差转成均匀分布分位数 uAll zeros(size(epsAll)); for i 1:n [~, F] ecdf(epsAll(:,i)); % 经验CDF uAll(:,i) F; % 注意边界值后面处理 end % 逆正态变换 zAll norminv(uAll); % 空间相关矩阵 Sigma corr(zAll); disp(Sigma);得到的Sigma就是Gaussian Copula的相关矩阵。矩阵元素在0.3以上基本能确认存在不可忽略的空间相关。严谨一点可以进一步用似然比检验比较“相关矩阵为单位阵”和“自由相关矩阵”两个模型的拟合度不过实际项目中只要矩阵里有明显的大元素相关性就是跑不掉的。4. 时空相关误差场景生成从参数估计到Monte Carlo模拟4.1 分步建模方案AR残差空间Copula我推荐一个特别适合工程落地的方案分四步走。第一步逐场站对归一化误差序列拟合AR(p)模型得到持续性系数和残差序列。第二步对残差序列用Gaussian Copula建模得到空间相关矩阵Sigma。第三步生成独立的多元正态样本R mvnrnd(zeros(1,n), Sigma, Nsim)转成均匀分布再用各场站残差的经验分布逆函数映射回残差。第四步对每个模拟样本用AR模型递归重建时间序列e_sim(t) φ · e_sim(t-1) ε_sim(t)这样生成的模拟误差序列既有时间维度上的持续偏差又有空间维度上的同步偏差而且两个维度的结构是分开估计的互不干扰。为什么不用一个完整的多元AR模型说实话如果场站数量少3到5个Matlab里用多元AR也没问题但场站一多参数数量会爆炸而且容易出现数值不稳定。ARCopula的分步方案在工程上更好掌控每个环节都能单独验证。4.2 Matlab核心代码拟合、采样、反变换下面给出一段完整的场景生成核心代码。假设errAll是T行n列矩阵T是时间点数n是场站数且每列都是归一化误差。rng(2025); % 固定随机种子保证结果可复现 n size(errAll, 2); Nsim 5000; % 场景数 phi zeros(1, n); epsResid zeros(size(errAll)); % 第一步逐场站AR(1)拟合 for i 1:n arModel arima(1,0,0); mdl estimate(arModel, errAll(:,i), Display, off); phi(i) mdl.AR{1}; epsResid(:,i) infer(mdl, errAll(:,i)); end % 第二步残差边缘分布 Gaussian Copula uResid zeros(size(epsResid)); for i 1:n [~, F] ecdf(epsResid(:,i)); uResid(:,i) F; uResid(:,i) min(max(uResid(:,i), 1e-6), 1-1e-6); % 处理边界 end zResid norminv(uResid); Sigma corr(zResid); % 第三步从Copula采样得到残差场景 % 先生成标准正态相关样本再转回均匀空间 zSim mvnrnd(zeros(1,n), Sigma, Nsim); uSim normcdf(zSim); % 用各场站残差的经验逆CDF映射 epsSim zeros(Nsim, n); for i 1:n [~, F] ecdf(epsResid(:,i)); epsSim(:,i) icdf(F, uSim(:,i)); % 或者用quantile(epsResid(:,i), uSim(:,i)) end % 第四步AR(1)递归重建误差时序 Tsim 24; % 假设要模拟未来24个时点 errSim zeros(Tsim, Nsim, n); for i 1:n for k 1:Nsim e zeros(Tsim, 1); e(1) epsSim(k, i) / sqrt(1 - phi(i)^2); % 初始值取稳态分布 for t 2:Tsim e(t) phi(i) * e(t-1) epsSim(k, i); end errSim(:, k, i) e; end end代码里有个细节要解释一下初始值的处理。AR(1)在t1时需要给一个初始状态我用的分母sqrt(1-φ²)也就是说初始误差从AR的平稳分布里随机抽取避免初始值偏差影响前几个时点的场景质量。如果你对前几个时点不敏感也可以直接令e(1)0但那样模拟的前两个小时误差会偏小后面才恢复正常。4.3 场景质量评估CRPS、PICP、PINAW模型建完之后必须量化评估不能只看相关性系数好看。我的评估体系固定用三个指标。CRPS连续排名概率分数评价概率预测整体质量。给定预测分布和实测值CRPS同时惩罚分布不集中和分布不准确越小越好。在场景生成语境下可以用经验CDF近似计算CRPS。PICP预测区间覆盖率反映置信区间是否可靠比如95%区间应该覆盖大约95%的实测点覆盖率过低说明模型过度自信。PINAW平均区间宽度衡量区间是否实用覆盖率都达标的情况下越窄越好。三个指标配合使用才不会被单方面带偏。有的模型PICP很高但区间宽到对调度没有任何参考价值有的模型区间很窄但覆盖率一塌糊涂。工程上及格线是95%置信区间覆盖率达到92%以上区间宽度在可接受范围CRPS比独立采样模型明显下降。5. 一个三场站算例时空相关模型比独立采样好在哪里5.1 算例数据说明和建模设定为了把这个方法讲得接地气我编一个贴近真实项目形态的算例。假设区域内有3个风电场装机容量分别是100MW、120MW、80MW总共300MW。收集了一整年15分钟分辨率的预测与实测数据预测来源是某NWP加统计后处理。训练集用前10个月测试集用最后2个月。三种对比模型M1独立同分布采样各场站误差用高斯分布拟合场站之间完全不相关。 M2仅考虑时间相关各场站AR(1)拟合但场站之间独立。 M3完整时空模型AR(1)Gaussian Copula空间相关。评估的内容是测试集里每3小时总出力误差的分布。为什么要看3小时而不是15分钟因为调度决策和备用评估通常看更长时间尺度的累计误差而且时间聚合会让独立性假设的错误更明显。5.2 三种模型的结果对比用训练好的模型各生成5000个3小时时段的总误差样本和测试集的实测误差分布做对比。结果很典型我这里给一个大致数值模型总误差标准差MW95%区间宽度MW覆盖测试集比例总误差CRPS实测67.8约260——M1独立高斯45.617581.2%偏高M2仅时间AR55.221088.4%中等M3时空ARCopula64.524894.7%最低M1把总误差标准差低估了32%以上这正好印证前边说的问题独立假设下的总方差只有全部方差之和真实情况下空间正相关会让总方差变大。M2把时间相关加进去后有一定改善但因为忽略了场站之间的共同波动低估仍明显。M3的结果和实测最接近95%区间覆盖率也达到了合理水平。5.3 结果解读总误差的低估风险与场景多样性这个算例背后的含义值得反复说。如果你用M1去配置备用容量按95%分位数算需要留约90MW的备用而实际上95%分位数在250MW左右。差了将近三倍这是有安全风险的。很多风电项目的偏差考核罚款也是同样道理独立模型算出来的惩罚成本期望值明显偏低报价策略会过于激进。还有一点很多文章不提时空相关模型生成的多场站场景比独立采样场景“更像真实情况”。独立采样生成的情况经常是A场站偏大、B场站偏小相互抵消总误差不大。而实际数据里三个场站经常同时偏、同时偏大或同时偏小。M3生成的场景保留了这种同步性拿去做储能充放电策略优化时优化出来的决策在极端场景下表现好得多。6. 实操中容易踩的六个坑含Matlab边角问题6.1 相关矩阵非正定与近邻修正这是我在多场站建模时遇到最多的报错。mvnrnd要求Sigma必须是半正定矩阵而样本量不足、数据缺失较多、或者场站间相关性过高时corr计算得到的矩阵偶尔不是正定的尤其是场站数量超过大致样本长度除以某个倍数后小样本伪相关会出现负特征值。解决办法有两个。一是加脊正则化Sigma_reg (1-λ)·Sigma λ·Iλ取0.01到0.1能有效稳定矩阵。二是用nearestSPD这类工具把矩阵投影到最近的正定矩阵。我在Matlab里更常用前一种因为它简单而且物理意义可解释加入一个小量的单位阵相当于给每个场站加入一点点独立噪声对相关性结构的影响很小。6.2 边界值处理ecdf和icdf的坑经验分布函数有一个麻烦如果误差序列里有极端值ecdf算出的累积概率可能等于1或者在最小值处可能等于0。但Gaussian Copula的参数估计要求u必须在(0,1)开区间内norminv在u0或1时会返回无穷大导致后面corr算出来全是NaN。处理方法就是我在代码里写的那句min(max(u, 1e-6), 1-1e-6)。这个微小的截断对结果影响可以忽略但能让整个管线稳定跑完。还有用icdf(F, u)前要确认F是经验CDF对象直接用ecdf函数输出的F有时候是嵌套矩阵容易搞混。我建议用[f, x] ecdf(x)然后自己写线性插值或者干脆用quantile(x, u)代替逆CDF效果差不多。6.3 边缘分布到底选核密度还是参数分布这个问题我在项目里来回换过几次。核密度估计灵活能捕捉任意形状但样本量不足时尾部噪声很大而且逆变换求分位数时可能给出不稳定的极端值。参数分布稳定但容易对真实形状做过度简化。目前我的经验法则是训练样本大于等于一年且点数超过2万用核密度或者分位数映射都行样本量小优先试Beta分布或者偏斜t分布。Beta分布有一个额外好处它天然定义在[0,1]区间正好对应归一化到装机容量区间的误差比例。但注意Beta分布要求误差值严格在[0,1]内实际数据常会出现负误差或略大于1的值需要预先做移位和截断errMapped (err - min(err)) / (max(err) - min(err)); errMapped min(max(errMapped, 1e-6), 1-1e-6); pd fitdist(errMapped, Beta);6.4 固定相关矩阵的局限性分状态建模更靠谱最后分享一个我自己研究中后期才意识到的问题。空间相关矩阵并不是一成不变的大风天气、小风天气、静稳天气下各场站预测误差的相关性明显不同。大风天气下系统性强所有场站误差同步性高小风天气下预测通常偏保守场站间相关弱一些。如果全年用同一个相关矩阵本质上是在平均不同天气状态下的相关结构。更合理的做法是分状态建模按预测风速或预测功率水平把样本分成3到4个区间每个区间分别估计AR系数和Copula矩阵。比如风速小于5m/s的区间用一套参数5到10m/s一套大于10m/s一套。这样模型会稍微复杂些但误差分布和相关性都能更贴合实际天气过程。我在三场站算例里对比过分状态建模之后CRPS还能再降大约5%到8%。6.5 Matlab版本和工具箱函数的兼容问题不同Matlab版本之间arima模型的estimate和infer接口相对稳定但copulafit这个函数在部分老版本里要求输入u矩阵的每行是一个观测、每列是一个变量在版本之间有时会提示错误。我会提前检查要用的函数是否存在比如用exist(copulafit,file)确认。还有一个老生常谈的问题如果只装了基础版没装Statistics and Machine Learning Toolboxecdf、fitdist、mvnrnd这些函数全部不可用建议装好工具箱再跑别浪费时间纠结代码。6.6 验证时不要把“相关性提高”和“模型变好”划等号最后这条算一个提醒。有时候把相关性加进模型模拟出的误差方差自然就大了PICP也会上升看起来模型“对了”。但注意PICP上升可能只是因为区间变宽并不代表模型预测得更准。一定要同时看PINAW和CRPS。CRPS这个指标对整体分布质量很敏感如果加了Copula后CRPS没有明显下降反而说明相关结构可能被过度拟合了需要用更简单的独立模型。我自己做下来最大的体会是时空相关误差建模最难的并不是算法而是想清楚每一步在回答什么问题。边缘分布回答“每个场站偏多少”AR系数回答“偏差能持续多久”Copula矩阵回答“多个场站会不会一起偏”。把这三个问题分别建模、分别验证整个项目就不会糊成一团。最后再补一句实操技巧场景生成数量建议至少5000个太少的话总误差分布尾部的抖动会很大参数估计用滚动窗口定期更新才不会在季节转换后模型表现突然变差。
阅读完成 · 觉得有帮助？
