简介面向电力系统调度、风电功率预测及时间序列分析学习者提供基于MATLAB实现的ARMA模型示例用于对风速或风电功率进行短期预测。风速序列通常具有自相关性和随机波动ARMA模型可有效刻画此类动态特征为电网调度提供参考。压缩包内共1个文件为MATLAB脚本文件包体仅2KB结构紧凑适合快速阅读、调试并迁移到自己的风电场数据集中。脚本覆盖从数据预处理、自相关与偏自相关分析定阶到使用arima函数完成参数估计、残差检验再到输出预测结果的关键流程完整展现了ARMA(p,q)模型的建模步骤。已有470人学习可帮助读者理解自回归与移动平均模型的机理掌握风电场数据时序建模的基本方法并作为后续扩展ARIMA或组合预测模型的起点代码也保留了较清晰的注释阅读时能快速对应每个建模环节降低从数学原理到实际运行之间的转换门槛。1. 别急着estimateARMA1.zip里的风速预测为什么值得先跑通我手头这份ARMA1.zip里只有两个文件ARMA1.zip和ARMA1.m。解压后其实就是一个ARMA1.m用MATLAB实现对风电场风速或风电功率的预测。很多人拿到这类资源第一反应是改改路径跑一下出图就完事。但风功率预测比普通时序预测更敏感风速序列有明显的日周期性、湍流噪声偶尔还有测风仪结冰造成的跳零。如果直接把原始数据送进arima函数结果往往残差里还藏着周期项。这篇笔记就把我在风电数据上跑ARMA的完整流程写出来从平稳性检验到定阶、参数估计、预测最后是滚动预测时容易踩的坑。适合刚接触MATLAB时间序列工具箱的工程师也适合想从SVM/神经网络切回经典模型对比基线的人。2. 风速数据清洗与平稳性检验先让ARMA的前提成立2.1 为什么我建议用10分钟平均风速而不是瞬时值ARMA模型的数学基础是弱平稳时间序列均值恒定、协方差只依赖时间间隔。实测风速每秒都在波动瞬时值序列几乎都是非平稳的而且带有严重的异方差。风电场常用的SCADA数据通常提供1秒或10秒采样值我一般先重采样成10分钟平均风速。原因有二一是10分钟平均基本能反映湍流尺度以外的变化二是风功率曲线IEC标准也是按10分钟平均定义的。重采样用retime如果数据是timetable格式或者直接reshape加mean。注意不要在重采样时丢掉NaNretime默认会保留后续要处理。% 假设data是分钟级风速时间戳在t列风速在v列 data readmatrix(wind_data.csv); % 第一列时间序号第二列风速m/s v_raw data(:,2); v_10min movmean(v_raw, 10, Endpoints, fill); % 简单10点滑动平均实际按时间窗这里movmean是滑动平均参数10表示窗口长度Endpoints,fill保证首尾不收缩。严格按10分钟间隔应该对每600秒的数据块求平均但快速验证时滑动窗口更方便。处理完风速后风电功率可以做同样的重采样但功率序列还会受风机控制策略影响经常有截至风速以上的平台区ARMA对这类平台区的拟合会偏保守。2.2 缺失值、异常值和野点哪些要插值哪些要剔除测风塔数据最常见的两个问题一个是通信中断导致的整段缺失另一个是叶片结冰或仪器故障导致的长时间零值或尖峰值。整段缺失超过2小时的我选择删除而不是插值因为线性插值会人为引入低频趋势影响ADF检验判断。单个离散野点比如风速在0.5秒内突变超过5m/s用filloutliers处理方法选linear。注意filloutliers默认用中位数绝对偏差对风数据可能过敏感我会先用isoutlier看一遍。v_clean filloutliers(v_10min, linear, ThresholdFactor, 5); % ThresholdFactor越大越不敏感 v_clean fillmissing(v_clean, movmedian, 6); % 对剩余缺失用窗口长度6的中位数填补参数说明ThresholdFactor是基于MAD的倍率风电数据取4到6比较合适取3会误删正常阵风。fillmissing的movmedian窗口我一般取6等效补1小时内的缺口超过1小时的缺口还是建议直接滤掉该段。处理完之后要做一次平稳性检验不要看折线图差不多就往下走。数据异常类型判别方式处理方式适用场景瞬时尖峰相邻点差5m/sfilloutliers(linear)测风仪误差长时间零值连续超过2小时为0整段剔除叶片结冰随机缺失少于6个连续NaNfillmissing(movmedian)通信丢包长段缺失超过2小时连续NaN剔除并断开通道故障这个策略不是固定的。长段缺失如果后续预测需要可以用ACF估计出的前一天同一时段数据做占位但那样会引入对称偏差所以我在短期样本里宁可少一段也不能掺水分。表格里的判别方式是我在代码里实际用的条件超过这个阈值就该人工review。2.3 ADF检验和差分阶数写脚本而不是肉眼判断ADF检验原假设是序列存在单位根即非平稳。MATLAB的adftest返回h和pValueh1表示拒绝原假设序列平稳。我习惯写一个小循环对不同差分阶数分别检验maxD 2; results zeros(maxD1, 2); for d 0:maxD if d 0 y v_clean; else y diff(v_clean, d); end y y(~isnan(y)); % 差分后首部出现NaN [h, p] adftest(y, Alpha, 0.05); results(d1,:) [h, p]; end这个循环输出d0、1、2对应的判断结果。多数情况原始10分钟平均风速的p值在0.1到0.5之间不能拒绝单位根说明非平稳一阶差分后p值通常小于0.01h1。如果一阶差分后仍不平稳先不要急着二阶差分大概率是存在周期性应该考虑SARIMA或者加入傅里叶项。风功率序列比风速更容易出现“团状波动”二阶差分可能过度差分手流失真信息所以结果表里如果d2才平稳我会警惕数据里是否含着另一台机组的启停残影。提示adftest只测单位根不测季节性。如果风速数据有明显日周期d1后ACF图还会在滞后96如果是10分钟数据每24小时为144个点附近有显著峰值这时候要回到数据层面做季节差分而不是无限增加差分阶数。3. MATLAB里用arima()定阶ACF/PACF和信息准则怎么配合3.1 自相关函数和偏自相关函数读图ARMA(p,q)模型写作 y_t c φ1 y_{t-1} ... φp y_{t-p} ε_t θ1 ε_{t-1} ... θq ε_{t-q}。识别阶数最经典的工具是ACF和PACF。对于差分后的平稳序列画autocorr和parcorry diff(v_clean); % 一阶差分后的序列 y y(~isnan(y)); figure subplot(2,1,1) autocorr(y, NumLags, 60) title(ACF - First Difference) subplot(2,1,2) parcorr(y, NumLags, 60) title(PACF - First Difference)NumLags取60对于10分钟数据相当于10小时能看出日内趋势。读图规律是ACF拖尾、PACF截尾则选AR(p)ACF截尾、PACF拖尾则选MA(q)两边都拖尾则考虑ARMA(p,q)。风速序列很少出现干净的截尾更多是两边拖尾所以直接用ARMA组合而非纯AR或纯MA。另外ACF在48、144位置有峰说明还有周期成分需要先用lag算子或者季节差分不能直接套ARMA。3.2 AIC/BIC自动定阶遍历p、q而不是凭经验读图只能缩小范围。我通常把p和q都限定在0到5遍历后按AIC、BIC排序。MATLAB里可以对arima模型用estimate再取aic但估计161个组合耗时有点久。更快的办法是直接用armax不这里用arima但可以先用aictestMATLAB没有直接函数写循环即可。maxPQ 5; aicVals zeros(maxPQ1, maxPQ1); bicVals zeros(maxPQ1, maxPQ1); for p 0:maxPQ for q 0:maxPQ if p 0 q 0, continue; end Mdl arima(p,1,q); % d1直接对原始风速序列建模 try [EstMdl, ~, ~] estimate(Mdl, v_clean, Display, off); aicVals(p1,q1) EstMdl.Summary.AIC; bicVals(p1,q1) EstMdl.Summary.BIC; catch aicVals(p1,q1) inf; bicVals(p1,q1) inf; end end end [minAIC, idx] min(aicVals(:)); [pBest, qBest] ind2sub(size(aicVals), idx); pBest pBest - 1; qBest qBest - 1;这段代码直接对原始风速做差分阶数d1的ARMA估计遍历p、q。estimate会自动处理NaN不v_clean不应该有NaN。try是因为高阶模型可能估计失败返回inf。ind2sub把线性索引恢复成p和q注意MATLAB的p,q从0开始所以减1。这样选出来的最优阶数可以直接用于预测。不过我要提醒一点summary中的AIC字段在不同MATLAB版本有差异如果EstMdl.Summary.AIC取不到就改用aicbic(LogL, numParams, numObs)。我遇到过2020a之前的版本没有Summary对象所以备错方式要写好。3.3 为什么BIC选出的阶数常常比AIC低AIC更注重拟合优度BIC对阶数惩罚更重。在风电数据上AIC容易选出p5,q5的过度参数化模型样本外预测前几步波动剧烈。我一般以BIC为第一参考然后看估计出的系数是否显著。下面是我跑过的一组对比数据能说明问题候选模型AICBICLjung-Box p值结论ARMA(2,1)10234.110210.70.203可选ARMA(5,4)10211.510232.90.612过拟合风险ARMA(1,1)10302.810295.40.011残差有自相关这个表是示意但规律是真的ARMA(5,4)的AIC最低BIC却更高因为多出来的滞后项没有实际贡献。ARMA(1,1)的Ljung-Box p值0.011残差里还有信息不可取。折中后选ARMA(2,1)BIC较小残差检验通过。注意不要看到AIC最低就直接用要把候选模型的残差图叠起来看尤其看预测步数4-8步的误差。4. 参数估计与残差白噪声检验模型有没有学成4.1 estimate的基本参数和使用习惯选好了p、q用estimate估计。我一般会打开Display同时保留结果对象。Mdl arima(pBest, 1, qBest); Mdl.Constant 0; % 风速差分后均值接近0通常不设常数项 [EstMdl, estParams, LogL] estimate(Mdl, v_clean, Display, params);参数说明Mdl.Constant 0可以减小参数方差但如果原始风速有明显上升趋势还是保留常数项。estParams里包含估计值、标准误和t统计量。注意如果Mdl没有指定AR和MA滞后默认是AR(1)和MA(1)不是AR(p)和MA(q)所以要先明确arima(p,1,q)。还有一个细节v_clean的末尾不能有NaNestimate虽然会跳过NaN但差分后的起始点也会少一个导致样本量突然减少影响AIC比较时的自由度。所以在前一章数据清洗时我会把所有结果全部整理成不含NaN的均匀序列。4.2 怎么读估计结果AR系数、MA系数和常数项假设输出Constant -0.0023 0.0011 -2.091 AR{1} 0.7821 0.0321 24.37 AR{2} -0.1023 0.0445 -2.30 MA{1} 0.4213 0.0632 6.66 Variance 0.5834 0.0116 50.31AR{1}接近0.8说明差分序列有较强的惯性MA{1}接近0.4说明白噪声冲击的影响会持续一段时间。这些系数都在±1之间实质根在单位圆内说明模型稳定。MATLAB对AR多项式要求所有根在单位圆外如果估计结果里有根的模小于1infer会报警。这里可以跑一下[res, logL] infer(EstMdl, v_clean);infer返回残差序列注意它是在给定模型和数据条件下计算一步预测残差比resid更适合后续诊断。logL是似然函数值用于比较非嵌套模型。res前几个值可能是NaN因为AR项需要前序观测。我一般会丢掉前10个NaN再做检验因为lbqtest不会自动忽略NaN传进去会报错。4.3 残差白噪声检验Ljung-Box的lags到底选多少残差如果还存在自相关说明模型没提取干净。常用lbqtestresValid res(~isnan(res)); [h, p, stat] lbqtest(resValid, Lags, [5, 10, 20, 40]);返回值h是0/1向量对应每个lags。Lags选多个值不要只看一个。风速序列的滞后5和滞后40含义完全不同滞后5代表半小时内的短记忆滞后40代表约7小时前的周期性残留。我要求所有h0p0.05。如果滞后40对应p很小通常说明每天的峰值时段没有被ARMA捕捉考虑引入时间变量或改用SARIMA。还有一种常见问题是残差方差不是常数风电功率在大风天和小风天方差差异明显lbqtest检验不出来可以画autocorr(res.^2)看平方残差是否相关。若相关则说明有ARCH效应需要GARCH类模型而不是直接上ARMA-GARCH组合。这里要强调一下AR和MA阶数不能无限提高来压残差那会把噪声也学进模型。4.4 残差正态性与预测置信区间的关系lbqtest通过后很多人不管正态性。但forecast生成的置信区间是基于残差正态假设的。检查qqplot或jbtest。风速残差通常略带厚尾置信区间会略乐观问题不大。但若残差分布明显偏态说明模型倾向于低估上升速度。可以看skewness(res)若为负则预测风速往往偏大。这时我会在预测输出后加一点保守修正比如把预测下限再压低5%。不过这是经验做法不是统计学结论。实际上在风功率场景里我更关心残差的方差是否随风速增大而增大所以我会按风速分箱计算残差方差如果高风速箱的方差是低风速箱的3倍以上就不要再简单使用同方差置信区间了。5. 风速/风电功率预测与误差评估forecast()怎么用5.1 用估计好的模型做多步预测forecast函数可以给出点预测和预测误差方差但要注意需要提供前导数据。如果直接用forecast(EstMdl, 48, Y0, v_clean)它会输出未来48步预测。这里的48步是10分钟即8小时。实际风电调度需要未来24小时甚至72小时预测这时ARMA的预测误差很快变大因为ARMA本质是短记忆模型。因此我会把预测步数限制在6-12步再多就开始退化。numSteps 12; % 2小时预测 [YF, YMSE] forecast(EstMdl, numSteps, Y0, v_clean); Ypred YF; % 直接就是原尺度风速因为模型内部处理了差分 Yerr sqrt(YMSE);实际上forecast返回的YF是原始数据尺度还是差分后尺度当模型有d0时forecast会返回原始尺度这是MATLAB中arima的方便之处内部做差分的逆变换。但要小心Y0的长度至少等于模型的最大AR滞后阶数如果不够forecast会报错或自动截断。YMSE是预测误差方差不是置信区间半宽。95%置信限为YF ± 1.96*sqrt(YMSE)不要直接拿YMSE当误差棒。5.2 滚动预测还是直接多步预测调度业务里的选择在风电场我通常固定模型每10分钟滚动预测未来2小时。直接多步预测的误差在12步之后会快速膨胀因为ARMA的冲击衰减很快。滚动方式更贴近业务YpredAll zeros(numSteps,1); inputData v_clean; for k 1:numSteps if k 1 inputData [inputData; YpredAll(k-1)]; % 把上一步预测值当作观测 end [YpredAll(k), ~] forecast(EstMdl, 1, Y0, inputData); end这个循环里把前一步预测当作输入会累积误差但能模拟在线更新。另一种做法是每一次都用最新观测值重新估计模型但每10分钟一次estimate太慢常见做法是每小时重新估计一次参数中间几次用固定模型滚动。注意如果inputData过长每次循环都要重新截取最近500点否则forecast内部会处理整个历史序列影响计算速度。我一般这样写inputData v_clean(max(1,end-500):end); % 只留最近500个点5.3 误差指标MAE、RMSE和“风速≤6m/s时的偏差点”不能只看RMSE。风电功率预测更关心拐点和低风速段。功率是风速的三次方风速误差5%对应功率误差约16%。因此我用三段式评估风速区间样本比例RMSE(m/s)对应功率误差(MW)0-4 m/s20%0.430.084-12 m/s55%0.610.4212 m/s25%0.780.35这里数值是示意但思路是评估必须分区间。低于切入风速时ARMA预测的绝对值误差很小但相对误差很大大于额定风速时风速变化对功率影响钝化。如果有功率曲线最好把风速预测转换成功率后再算RMSE不要直接对功率建模因为功率序列的分布更像尖峰厚尾ARMA假设不成立。5.4 预测结果可视化画图对比实测和预测figure plot(1:numSteps, Yreal, k-o, LineWidth, 1.2); hold on; plot(1:numSteps, Ypred, r--, LineWidth, 1.2); upper Ypred 1.96*Yerr; lower Ypred - 1.96*Yerr; fill([1:numSteps flip(1:numSteps)], [upper flip(lower)], ... r, FaceAlpha, 0.2, EdgeColor, none); legend({Actual,Forecast,95% PI}); xlabel(Step (10 min)); ylabel(Wind speed (m/s));这里fill的用法要注意flip让上界和下界围成多边形FaceAlpha控制透明。注意upper和lower如果维度不对需要转置。预报区间在风电调度里比点预测更有用如果区间宽度大于1.5m/s调度员会倾向于安排备用容量。我通常在图上再画一条功率曲线对应的功率轴但这里保持风速轴更清晰。6. 滚动预测技巧避免模型失效的三个检查点6.1 每次rolling前先检查数据尾部的质量滚动预测最容易踩坑的是测风塔在预测开始时刻出现一段NaNforecast会忽略尾部NaN导致输入长度减少但步数不变预测值直接偏掉。我的检查句式assert(sum(isnan(v_clean(end-30:end))) 0, 尾部有NaN需要先平滑填补)如果尾部有缺失用前一个10分钟的观测值填充不要用全段均值。因为ARMA中最近的历史值权重最大用全段均值会把最近的风速突变拉平。6.2 固定窗口重估参数ARMA参数会随季节变化。夏天热力湍流强AR系数的惯性低于冬天。我会维护一个每小时重估的定时器% 伪代码实际在回调函数中 if mod(minute(now), 60) 0 EstMdl estimate(Mdl, v_clean(end-720:end), Display,off); end这里窗口取720个点即5天。窗口太短参数波动大太长跟不上季节变化。参数重估后要重新做一次lbqtest如果p值突然小于0.05说明数据进入新的状态。注意mod用minute(now)不是标准函数实际建议用mod(posixtime(datetime(now)), 3600)来控制否则在某些MATLAB版本里now返回日期序列号minute只能取当前分钟同一个小时的不同分钟都会进入导致重估频繁。这里我给的是逻辑示意部署时要知道这个坑。6.3 模型退化指标预测误差连续上升时的退出策略只靠RMSE判断退化会有滞后。我习惯记录每步预测的误差相对前一天的同一时刻的变化率超过阈值就把模型回退到上一步参数。这个阈值取什么风速预测中MAE连续6个点超过历史P90就触发重估。不要等全天结束。如果重估后依然不佳说明ARMA已经不合适应该切换支持向量回归或LSTM做对照。ARMA作为强基线在数据稳定时往往和深度学习差距不大但计算成本低一个数量级这个优势在线上部署时很值钱。if sum(maeWindow prctile(maeHistory, 90)) 6 reestimate true; end这个阈值需要现场根据业务调。比如冬夏两季的P90不同我会把maeHistory按月份分组避免跨季节比较。重估时使用卡尔曼滤波器可能更平滑但estimate在每小时一次的计算成本完全可接受。把这个阈值记录在配置文件里现场值班人员可以直接调整不需要改代码。本文还有配套的精品资源点击获取
MATLAB ARMA风速预测：从平稳性检验到滚动预测全流程
简介面向电力系统调度、风电功率预测及时间序列分析学习者提供基于MATLAB实现的ARMA模型示例用于对风速或风电功率进行短期预测。风速序列通常具有自相关性和随机波动ARMA模型可有效刻画此类动态特征为电网调度提供参考。压缩包内共1个文件为MATLAB脚本文件包体仅2KB结构紧凑适合快速阅读、调试并迁移到自己的风电场数据集中。脚本覆盖从数据预处理、自相关与偏自相关分析定阶到使用arima函数完成参数估计、残差检验再到输出预测结果的关键流程完整展现了ARMA(p,q)模型的建模步骤。已有470人学习可帮助读者理解自回归与移动平均模型的机理掌握风电场数据时序建模的基本方法并作为后续扩展ARIMA或组合预测模型的起点代码也保留了较清晰的注释阅读时能快速对应每个建模环节降低从数学原理到实际运行之间的转换门槛。1. 别急着estimateARMA1.zip里的风速预测为什么值得先跑通我手头这份ARMA1.zip里只有两个文件ARMA1.zip和ARMA1.m。解压后其实就是一个ARMA1.m用MATLAB实现对风电场风速或风电功率的预测。很多人拿到这类资源第一反应是改改路径跑一下出图就完事。但风功率预测比普通时序预测更敏感风速序列有明显的日周期性、湍流噪声偶尔还有测风仪结冰造成的跳零。如果直接把原始数据送进arima函数结果往往残差里还藏着周期项。这篇笔记就把我在风电数据上跑ARMA的完整流程写出来从平稳性检验到定阶、参数估计、预测最后是滚动预测时容易踩的坑。适合刚接触MATLAB时间序列工具箱的工程师也适合想从SVM/神经网络切回经典模型对比基线的人。2. 风速数据清洗与平稳性检验先让ARMA的前提成立2.1 为什么我建议用10分钟平均风速而不是瞬时值ARMA模型的数学基础是弱平稳时间序列均值恒定、协方差只依赖时间间隔。实测风速每秒都在波动瞬时值序列几乎都是非平稳的而且带有严重的异方差。风电场常用的SCADA数据通常提供1秒或10秒采样值我一般先重采样成10分钟平均风速。原因有二一是10分钟平均基本能反映湍流尺度以外的变化二是风功率曲线IEC标准也是按10分钟平均定义的。重采样用retime如果数据是timetable格式或者直接reshape加mean。注意不要在重采样时丢掉NaNretime默认会保留后续要处理。% 假设data是分钟级风速时间戳在t列风速在v列 data readmatrix(wind_data.csv); % 第一列时间序号第二列风速m/s v_raw data(:,2); v_10min movmean(v_raw, 10, Endpoints, fill); % 简单10点滑动平均实际按时间窗这里movmean是滑动平均参数10表示窗口长度Endpoints,fill保证首尾不收缩。严格按10分钟间隔应该对每600秒的数据块求平均但快速验证时滑动窗口更方便。处理完风速后风电功率可以做同样的重采样但功率序列还会受风机控制策略影响经常有截至风速以上的平台区ARMA对这类平台区的拟合会偏保守。2.2 缺失值、异常值和野点哪些要插值哪些要剔除测风塔数据最常见的两个问题一个是通信中断导致的整段缺失另一个是叶片结冰或仪器故障导致的长时间零值或尖峰值。整段缺失超过2小时的我选择删除而不是插值因为线性插值会人为引入低频趋势影响ADF检验判断。单个离散野点比如风速在0.5秒内突变超过5m/s用filloutliers处理方法选linear。注意filloutliers默认用中位数绝对偏差对风数据可能过敏感我会先用isoutlier看一遍。v_clean filloutliers(v_10min, linear, ThresholdFactor, 5); % ThresholdFactor越大越不敏感 v_clean fillmissing(v_clean, movmedian, 6); % 对剩余缺失用窗口长度6的中位数填补参数说明ThresholdFactor是基于MAD的倍率风电数据取4到6比较合适取3会误删正常阵风。fillmissing的movmedian窗口我一般取6等效补1小时内的缺口超过1小时的缺口还是建议直接滤掉该段。处理完之后要做一次平稳性检验不要看折线图差不多就往下走。数据异常类型判别方式处理方式适用场景瞬时尖峰相邻点差5m/sfilloutliers(linear)测风仪误差长时间零值连续超过2小时为0整段剔除叶片结冰随机缺失少于6个连续NaNfillmissing(movmedian)通信丢包长段缺失超过2小时连续NaN剔除并断开通道故障这个策略不是固定的。长段缺失如果后续预测需要可以用ACF估计出的前一天同一时段数据做占位但那样会引入对称偏差所以我在短期样本里宁可少一段也不能掺水分。表格里的判别方式是我在代码里实际用的条件超过这个阈值就该人工review。2.3 ADF检验和差分阶数写脚本而不是肉眼判断ADF检验原假设是序列存在单位根即非平稳。MATLAB的adftest返回h和pValueh1表示拒绝原假设序列平稳。我习惯写一个小循环对不同差分阶数分别检验maxD 2; results zeros(maxD1, 2); for d 0:maxD if d 0 y v_clean; else y diff(v_clean, d); end y y(~isnan(y)); % 差分后首部出现NaN [h, p] adftest(y, Alpha, 0.05); results(d1,:) [h, p]; end这个循环输出d0、1、2对应的判断结果。多数情况原始10分钟平均风速的p值在0.1到0.5之间不能拒绝单位根说明非平稳一阶差分后p值通常小于0.01h1。如果一阶差分后仍不平稳先不要急着二阶差分大概率是存在周期性应该考虑SARIMA或者加入傅里叶项。风功率序列比风速更容易出现“团状波动”二阶差分可能过度差分手流失真信息所以结果表里如果d2才平稳我会警惕数据里是否含着另一台机组的启停残影。提示adftest只测单位根不测季节性。如果风速数据有明显日周期d1后ACF图还会在滞后96如果是10分钟数据每24小时为144个点附近有显著峰值这时候要回到数据层面做季节差分而不是无限增加差分阶数。3. MATLAB里用arima()定阶ACF/PACF和信息准则怎么配合3.1 自相关函数和偏自相关函数读图ARMA(p,q)模型写作 y_t c φ1 y_{t-1} ... φp y_{t-p} ε_t θ1 ε_{t-1} ... θq ε_{t-q}。识别阶数最经典的工具是ACF和PACF。对于差分后的平稳序列画autocorr和parcorry diff(v_clean); % 一阶差分后的序列 y y(~isnan(y)); figure subplot(2,1,1) autocorr(y, NumLags, 60) title(ACF - First Difference) subplot(2,1,2) parcorr(y, NumLags, 60) title(PACF - First Difference)NumLags取60对于10分钟数据相当于10小时能看出日内趋势。读图规律是ACF拖尾、PACF截尾则选AR(p)ACF截尾、PACF拖尾则选MA(q)两边都拖尾则考虑ARMA(p,q)。风速序列很少出现干净的截尾更多是两边拖尾所以直接用ARMA组合而非纯AR或纯MA。另外ACF在48、144位置有峰说明还有周期成分需要先用lag算子或者季节差分不能直接套ARMA。3.2 AIC/BIC自动定阶遍历p、q而不是凭经验读图只能缩小范围。我通常把p和q都限定在0到5遍历后按AIC、BIC排序。MATLAB里可以对arima模型用estimate再取aic但估计161个组合耗时有点久。更快的办法是直接用armax不这里用arima但可以先用aictestMATLAB没有直接函数写循环即可。maxPQ 5; aicVals zeros(maxPQ1, maxPQ1); bicVals zeros(maxPQ1, maxPQ1); for p 0:maxPQ for q 0:maxPQ if p 0 q 0, continue; end Mdl arima(p,1,q); % d1直接对原始风速序列建模 try [EstMdl, ~, ~] estimate(Mdl, v_clean, Display, off); aicVals(p1,q1) EstMdl.Summary.AIC; bicVals(p1,q1) EstMdl.Summary.BIC; catch aicVals(p1,q1) inf; bicVals(p1,q1) inf; end end end [minAIC, idx] min(aicVals(:)); [pBest, qBest] ind2sub(size(aicVals), idx); pBest pBest - 1; qBest qBest - 1;这段代码直接对原始风速做差分阶数d1的ARMA估计遍历p、q。estimate会自动处理NaN不v_clean不应该有NaN。try是因为高阶模型可能估计失败返回inf。ind2sub把线性索引恢复成p和q注意MATLAB的p,q从0开始所以减1。这样选出来的最优阶数可以直接用于预测。不过我要提醒一点summary中的AIC字段在不同MATLAB版本有差异如果EstMdl.Summary.AIC取不到就改用aicbic(LogL, numParams, numObs)。我遇到过2020a之前的版本没有Summary对象所以备错方式要写好。3.3 为什么BIC选出的阶数常常比AIC低AIC更注重拟合优度BIC对阶数惩罚更重。在风电数据上AIC容易选出p5,q5的过度参数化模型样本外预测前几步波动剧烈。我一般以BIC为第一参考然后看估计出的系数是否显著。下面是我跑过的一组对比数据能说明问题候选模型AICBICLjung-Box p值结论ARMA(2,1)10234.110210.70.203可选ARMA(5,4)10211.510232.90.612过拟合风险ARMA(1,1)10302.810295.40.011残差有自相关这个表是示意但规律是真的ARMA(5,4)的AIC最低BIC却更高因为多出来的滞后项没有实际贡献。ARMA(1,1)的Ljung-Box p值0.011残差里还有信息不可取。折中后选ARMA(2,1)BIC较小残差检验通过。注意不要看到AIC最低就直接用要把候选模型的残差图叠起来看尤其看预测步数4-8步的误差。4. 参数估计与残差白噪声检验模型有没有学成4.1 estimate的基本参数和使用习惯选好了p、q用estimate估计。我一般会打开Display同时保留结果对象。Mdl arima(pBest, 1, qBest); Mdl.Constant 0; % 风速差分后均值接近0通常不设常数项 [EstMdl, estParams, LogL] estimate(Mdl, v_clean, Display, params);参数说明Mdl.Constant 0可以减小参数方差但如果原始风速有明显上升趋势还是保留常数项。estParams里包含估计值、标准误和t统计量。注意如果Mdl没有指定AR和MA滞后默认是AR(1)和MA(1)不是AR(p)和MA(q)所以要先明确arima(p,1,q)。还有一个细节v_clean的末尾不能有NaNestimate虽然会跳过NaN但差分后的起始点也会少一个导致样本量突然减少影响AIC比较时的自由度。所以在前一章数据清洗时我会把所有结果全部整理成不含NaN的均匀序列。4.2 怎么读估计结果AR系数、MA系数和常数项假设输出Constant -0.0023 0.0011 -2.091 AR{1} 0.7821 0.0321 24.37 AR{2} -0.1023 0.0445 -2.30 MA{1} 0.4213 0.0632 6.66 Variance 0.5834 0.0116 50.31AR{1}接近0.8说明差分序列有较强的惯性MA{1}接近0.4说明白噪声冲击的影响会持续一段时间。这些系数都在±1之间实质根在单位圆内说明模型稳定。MATLAB对AR多项式要求所有根在单位圆外如果估计结果里有根的模小于1infer会报警。这里可以跑一下[res, logL] infer(EstMdl, v_clean);infer返回残差序列注意它是在给定模型和数据条件下计算一步预测残差比resid更适合后续诊断。logL是似然函数值用于比较非嵌套模型。res前几个值可能是NaN因为AR项需要前序观测。我一般会丢掉前10个NaN再做检验因为lbqtest不会自动忽略NaN传进去会报错。4.3 残差白噪声检验Ljung-Box的lags到底选多少残差如果还存在自相关说明模型没提取干净。常用lbqtestresValid res(~isnan(res)); [h, p, stat] lbqtest(resValid, Lags, [5, 10, 20, 40]);返回值h是0/1向量对应每个lags。Lags选多个值不要只看一个。风速序列的滞后5和滞后40含义完全不同滞后5代表半小时内的短记忆滞后40代表约7小时前的周期性残留。我要求所有h0p0.05。如果滞后40对应p很小通常说明每天的峰值时段没有被ARMA捕捉考虑引入时间变量或改用SARIMA。还有一种常见问题是残差方差不是常数风电功率在大风天和小风天方差差异明显lbqtest检验不出来可以画autocorr(res.^2)看平方残差是否相关。若相关则说明有ARCH效应需要GARCH类模型而不是直接上ARMA-GARCH组合。这里要强调一下AR和MA阶数不能无限提高来压残差那会把噪声也学进模型。4.4 残差正态性与预测置信区间的关系lbqtest通过后很多人不管正态性。但forecast生成的置信区间是基于残差正态假设的。检查qqplot或jbtest。风速残差通常略带厚尾置信区间会略乐观问题不大。但若残差分布明显偏态说明模型倾向于低估上升速度。可以看skewness(res)若为负则预测风速往往偏大。这时我会在预测输出后加一点保守修正比如把预测下限再压低5%。不过这是经验做法不是统计学结论。实际上在风功率场景里我更关心残差的方差是否随风速增大而增大所以我会按风速分箱计算残差方差如果高风速箱的方差是低风速箱的3倍以上就不要再简单使用同方差置信区间了。5. 风速/风电功率预测与误差评估forecast()怎么用5.1 用估计好的模型做多步预测forecast函数可以给出点预测和预测误差方差但要注意需要提供前导数据。如果直接用forecast(EstMdl, 48, Y0, v_clean)它会输出未来48步预测。这里的48步是10分钟即8小时。实际风电调度需要未来24小时甚至72小时预测这时ARMA的预测误差很快变大因为ARMA本质是短记忆模型。因此我会把预测步数限制在6-12步再多就开始退化。numSteps 12; % 2小时预测 [YF, YMSE] forecast(EstMdl, numSteps, Y0, v_clean); Ypred YF; % 直接就是原尺度风速因为模型内部处理了差分 Yerr sqrt(YMSE);实际上forecast返回的YF是原始数据尺度还是差分后尺度当模型有d0时forecast会返回原始尺度这是MATLAB中arima的方便之处内部做差分的逆变换。但要小心Y0的长度至少等于模型的最大AR滞后阶数如果不够forecast会报错或自动截断。YMSE是预测误差方差不是置信区间半宽。95%置信限为YF ± 1.96*sqrt(YMSE)不要直接拿YMSE当误差棒。5.2 滚动预测还是直接多步预测调度业务里的选择在风电场我通常固定模型每10分钟滚动预测未来2小时。直接多步预测的误差在12步之后会快速膨胀因为ARMA的冲击衰减很快。滚动方式更贴近业务YpredAll zeros(numSteps,1); inputData v_clean; for k 1:numSteps if k 1 inputData [inputData; YpredAll(k-1)]; % 把上一步预测值当作观测 end [YpredAll(k), ~] forecast(EstMdl, 1, Y0, inputData); end这个循环里把前一步预测当作输入会累积误差但能模拟在线更新。另一种做法是每一次都用最新观测值重新估计模型但每10分钟一次estimate太慢常见做法是每小时重新估计一次参数中间几次用固定模型滚动。注意如果inputData过长每次循环都要重新截取最近500点否则forecast内部会处理整个历史序列影响计算速度。我一般这样写inputData v_clean(max(1,end-500):end); % 只留最近500个点5.3 误差指标MAE、RMSE和“风速≤6m/s时的偏差点”不能只看RMSE。风电功率预测更关心拐点和低风速段。功率是风速的三次方风速误差5%对应功率误差约16%。因此我用三段式评估风速区间样本比例RMSE(m/s)对应功率误差(MW)0-4 m/s20%0.430.084-12 m/s55%0.610.4212 m/s25%0.780.35这里数值是示意但思路是评估必须分区间。低于切入风速时ARMA预测的绝对值误差很小但相对误差很大大于额定风速时风速变化对功率影响钝化。如果有功率曲线最好把风速预测转换成功率后再算RMSE不要直接对功率建模因为功率序列的分布更像尖峰厚尾ARMA假设不成立。5.4 预测结果可视化画图对比实测和预测figure plot(1:numSteps, Yreal, k-o, LineWidth, 1.2); hold on; plot(1:numSteps, Ypred, r--, LineWidth, 1.2); upper Ypred 1.96*Yerr; lower Ypred - 1.96*Yerr; fill([1:numSteps flip(1:numSteps)], [upper flip(lower)], ... r, FaceAlpha, 0.2, EdgeColor, none); legend({Actual,Forecast,95% PI}); xlabel(Step (10 min)); ylabel(Wind speed (m/s));这里fill的用法要注意flip让上界和下界围成多边形FaceAlpha控制透明。注意upper和lower如果维度不对需要转置。预报区间在风电调度里比点预测更有用如果区间宽度大于1.5m/s调度员会倾向于安排备用容量。我通常在图上再画一条功率曲线对应的功率轴但这里保持风速轴更清晰。6. 滚动预测技巧避免模型失效的三个检查点6.1 每次rolling前先检查数据尾部的质量滚动预测最容易踩坑的是测风塔在预测开始时刻出现一段NaNforecast会忽略尾部NaN导致输入长度减少但步数不变预测值直接偏掉。我的检查句式assert(sum(isnan(v_clean(end-30:end))) 0, 尾部有NaN需要先平滑填补)如果尾部有缺失用前一个10分钟的观测值填充不要用全段均值。因为ARMA中最近的历史值权重最大用全段均值会把最近的风速突变拉平。6.2 固定窗口重估参数ARMA参数会随季节变化。夏天热力湍流强AR系数的惯性低于冬天。我会维护一个每小时重估的定时器% 伪代码实际在回调函数中 if mod(minute(now), 60) 0 EstMdl estimate(Mdl, v_clean(end-720:end), Display,off); end这里窗口取720个点即5天。窗口太短参数波动大太长跟不上季节变化。参数重估后要重新做一次lbqtest如果p值突然小于0.05说明数据进入新的状态。注意mod用minute(now)不是标准函数实际建议用mod(posixtime(datetime(now)), 3600)来控制否则在某些MATLAB版本里now返回日期序列号minute只能取当前分钟同一个小时的不同分钟都会进入导致重估频繁。这里我给的是逻辑示意部署时要知道这个坑。6.3 模型退化指标预测误差连续上升时的退出策略只靠RMSE判断退化会有滞后。我习惯记录每步预测的误差相对前一天的同一时刻的变化率超过阈值就把模型回退到上一步参数。这个阈值取什么风速预测中MAE连续6个点超过历史P90就触发重估。不要等全天结束。如果重估后依然不佳说明ARMA已经不合适应该切换支持向量回归或LSTM做对照。ARMA作为强基线在数据稳定时往往和深度学习差距不大但计算成本低一个数量级这个优势在线上部署时很值钱。if sum(maeWindow prctile(maeHistory, 90)) 6 reestimate true; end这个阈值需要现场根据业务调。比如冬夏两季的P90不同我会把maeHistory按月份分组避免跨季节比较。重估时使用卡尔曼滤波器可能更平滑但estimate在每小时一次的计算成本完全可接受。把这个阈值记录在配置文件里现场值班人员可以直接调整不需要改代码。本文还有配套的精品资源点击获取