简介灰色模型实现的电力负荷预测源码包面向电力系统研究人员、数据分析初学者及需要做短期负荷预测的工程师解决小样本、非线性数据下的负荷趋势估计问题。资源共8个文件包括4个Matlab脚本m文件、2个Excel数据表、1张结果对比图与1份说明文档压缩包仅50KB轻量易用。已有272人学习下载。代码涵盖数据预处理、GM(1,1)模型构建、最小二乘参数求解、预测输出与残差检验等完整流程附带README说明和可视化图片可直接运行并对照理解灰色预测的核心步骤同时针对数据波动场景提供了模型改进与组合预测思路适合作为入门到进阶的实践模板。1. 灰色模型做电力负荷预测为什么数据越少越值得先试新台区投运、配电变压器扩容、工厂产线改造这类场景里历史用电数据往往只有几周的日峰值负荷。样本量这么小神经网络训练不起来统计回归又没有分布假设的支撑反倒是 Grey-Model灰色模型这类小样本方法用 4 个以上连续负荷点就能把短期趋势估出来。它在电力负荷预测里最常见的形态是 GM(1,1)对累加序列拟合一阶微分方程再把预测值累减还原成原始负荷。这种思路适合日电量、日最大负荷和未来 3~7 天的走势估算做电力系统规划、能效分析、运维调度的人都不需要依赖长历史曲线。下面从建模原理、级比检验、Python 实现到滚动窗口选参把一条能复现的路径讲完整。2. GM(1,1)的建模过程与可灰性判断先知道数据能不能用2.1 累加生成、背景值与最小二乘GM(1,1)建模的三步推导假设拿到了按时间顺序排列的负荷序列x(0) [x(0)(1), x(0)(2), …, x(0)(n)]。灰色模型不对原始序列直接拟合而是先做一次累加生成AGOx(1)(k) x(0)(1) x(0)(2) … x(0)(k)。累加后随机波动被摊薄。原始日负荷可能在 150~190 MW 之间来回跳累加序列却变成 156、318、476、646、814…… 这样一条接近指数形态的平滑曲线。GM(1,1) 中的第一个 1 表示一阶微分方程第二个 1 表示只有一个变量所以全称是一阶单变量灰色模型。对累加序列灰色模型采用一阶线性微分方程作为原型dx(1)/dt a·x(1) b。其中a称为发展系数b称为灰作用量。离散序列上算不出导数所以取相邻两个累加值的均值作为背景值z(1)(k) (x(1)(k) x(1)(k-1)) / 2把方程改写成x(0)(k) -a·z(1)(k) b。把 k2 到 n 的方程排成矩阵形式Y B·[a, b]^T。其中B矩阵第一列是-z第二列全为 1Y是原始序列第 2 到第 n 个点。用最小二乘解(B^T B)^(-1) B^T Y工程代码里通常写成伪逆。解出a、b后微分方程的解析解为x̂(1)(k1) (x(0)(1) - b/a)·e^(-a·k) b/a把k代入到n以后就得到累加域的预测值再用累减生成还原x̂(0)(k) x̂(1)(k) - x̂(1)(k-1)这里最容易出错的是指数上的符号e的指数是-a·k不是a·k。如果a是负数指数为正预测序列持续上升如果a为正数预测呈现衰减通常意味着数据本身处在下降段。2.2 级比检验与可容覆盖区间建模前先算这4个值不是任何负荷序列都能塞给 GM(1,1)。建模前先做级比检验。定义级比λ(k) x(0)(k-1) / x(0)(k)k 从 2 取到 n。指数型增长序列的相邻级比应该围绕某个常数小幅波动如果某个λ(k)明显偏离整体水平说明序列里有突变或周期跳动累加后偏离指数形态直接建模误差必然放大。理论上n 个样本的可容覆盖区间是(e^(-2/(n1)), e^(2/(n1)))。工程上很少手算直接查表更省事。样本数 n区间下限 e^(-2/(n1))区间上限 e^(2/(n1))60.75151.330770.77881.284080.80071.2488100.83361.1996120.85741.1663以[156, 162, 158, 170, 168, 176, 183]这 7 天日峰值负荷为例级比分别为 0.9630、1.0253、0.9294、1.0119、0.9545、0.9617全部落在 n7 的区间 (0.7788, 1.2840) 内说明该序列可以作为灰色模型的输入。检查时应重点看有没有个别级比越界如果越界点是检修日或错录的尖峰先对原始数据做清洗如果一半以上越界说明这段数据本身不适合灰色建模不要硬跑。级比计算用代码跑一遍更不容易出错import numpy as np seq np.array([156, 162, 158, 170, 168, 176, 183], dtypefloat) lmbda seq[:-1] / seq[1:] n len(seq) lower np.exp(-2 / (n 1)) upper np.exp(2 / (n 1)) print(级比序列:, np.round(lmbda, 4)) print(可容覆盖区间:, (round(lower, 4), round(upper, 4)))这里seq[:-1] / seq[1:]是在对相邻两项做比例计算得到的数组长度是n-1需要逐一对照上下限。注意如果在 pandas 里做务必先把索引重置为连续的 0..n-1否则切片会对不齐。提示级比检验建立在等间隔采样的前提上。日尺度负荷序列应保证逐日连续缺数日期要先插值不要用排序来代替时间顺序。2.3 发展系数a和灰作用量b在负荷场景中的实际含义估计出的a、b不只是中间变量它们的取值直接影响预测的可靠性。经验上|a|小于 0.3 时预测曲线平滑外推相对稳定0.3 到 0.5 之间适合做 1~3 步的短预测超过 0.8 时即使拟合阶段误差不大外推几步也会发散。负荷预测里如果算出很大的a先检查是否窗口太长或数据里有异常点而不是直接采用结果。b则承载原始序列的量纲和基准水平把单位从 MW 换成 kWb成倍变化a基本不变所以数据不需要像神经网络那样做 0-1 归一化保持原始单位和等间隔采样即可。3. 用Python实现Grey-Model电力负荷预测一套完整的GM(1,1)函数与检验3.1 最小可运行的GM(1,1)函数AGO、参数估计、累减还原一次完成无论网上找的 Grey-Model 包叫什么名字核心逻辑就集中在一个函数里累加生成、构造 B 矩阵、最小二乘求参、累减还原。下面这个实现只依赖 numpy可以直接复制运行import numpy as np def gm11(seq, pred_steps1): 灰色预测 GM(1,1) seq 原始一维负荷序列要求等间隔、全为正数长度 4 pred_steps 向后预测的步数 返回: fit 对历史数据的拟合值长度与 seq 相同 pred 未来 pred_steps 个预测值 a, b 模型参数 seq np.asarray(seq, dtypefloat) n len(seq) # 1. 一次累加生成 AGO摊平原始波动 x1 np.cumsum(seq) # 2. 背景值序列相邻累加值的均值 z (x1[1:] x1[:-1]) / 2.0 # 3. 构造方程组 Y B [a, b]^T B np.column_stack([-z, np.ones(n - 1)]) Y seq[1:].reshape(-1, 1) # 4. 最小二乘估计参数用 pinv 避免近奇异矩阵报错 theta np.linalg.pinv(B.T B) B.T Y a theta[0, 0] b theta[1, 0] # 5. 累加域预测 # x1_hat(k1) (x0(1) - b/a) * exp(-a*k) b/a k np.arange(1, n pred_steps) x1_hat np.r_[seq[0], (seq[0] - b / a) * np.exp(-a * k) b / a] # 6. 累减还原得到原始域的拟合值和预测值 fit np.r_[seq[0], np.diff(x1_hat[:n])] pred x1_hat[n:] - x1_hat[n - 1:-1] return fit, pred, a, b用 7 天日峰值负荷演示一次调用load [156, 162, 158, 170, 168, 176, 183] fit, pred, a, b gm11(load, pred_steps3) print(历史拟合:, np.round(fit, 2)) print(未来3天预测:, np.round(pred, 2)) print(a %.4f, b %.4f % (a, b))代码里有几点值得说明。第 4 步用pinv而不是inv当负荷序列接近线性时(B^T B)可能接近奇异用一般逆矩阵会抛出LinAlgError伪逆能稳定返回最小范数解。第 5 步生成的x1_hat长度为n pred_steps索引0对应原始第一个点索引1到n-1对应累加域拟合索引n以后对应未来累加值。第 6 步做差分时fit的第一个值保持为seq[0]预测值是x1_hat[n:] - x1_hat[n-1:-1]这就把累加域的差值还原成了原始负荷量纲。提示函数里没有做参数合法性校验。生产使用时先确认len(seq) 4、np.all(seq 0)否则后续计算没有意义。3.2 零值、负值和平移变换灰色模型数据预处理的3个规则灰色模型要求序列全为正数。负荷数据里的 0 值通常来自停电、仪表离线或表底未抄负值则可能是功率方向接反。出现 0 或负值级比计算会除零累加序列出现平台参数估计完全失真。处理规则单个 0 用前后两个时刻的均值补连续 0 用上一周同一天同时刻的负荷补实在补不齐就把这段数据截断用干净的子序列建模。如果序列整体是正数但个别点很小级比可能突破区间上限。常见做法是对序列做整体平移给每个点加一个常数偏移C建模预测后再把偏移减回去。简单实现如下def make_positive(seq): seq np.asarray(seq, dtypefloat) if np.all(seq 0): return seq, 0.0 offset -np.min(seq) 1.0 return seq offset, offset seq_shift, offset make_positive(load) fit_shift, pred_shift, a, b gm11(seq_shift, pred_steps3) fit fit_shift - offset pred pred_shift - offset平移变换只改变序列的位置不改变曲线的形状因此不影响发展系数a只影响b和还原后的预测值。要注意offset必须让平移后的最小值大于 0通常取-min(seq) 1.0就够用。3.3 后验差检验与精度等级C和P怎么看模型建完先做后验差检验不要只看拟合曲线贴得多近。检验指标有两个后验差比值C S2 / S1和小误差概率P。其中S1是原始序列的标准差S2是残差序列的标准差P表示残差与残差均值的偏差落在0.6745 * S1范围内的概率。计算代码如下def post_check(original, fitted): resid np.asarray(original) - np.asarray(fitted) S1 np.std(original, ddof1) S2 np.std(resid, ddof1) C S2 / S1 avg_resid np.mean(resid) P np.mean(np.abs(resid - avg_resid) 0.6745 * S1) return C, P C, P post_check(load, fit) print(后验差比值 C %.3f, 小误差概率 P %.3f % (C, P))精度等级的常见划分如下表。C 越小说明残差波动相对原始波动越弱模型捕捉趋势的能力越强P 越接近 1 说明残差越集中。精度等级后验差比值 C小误差概率 P工程结论一级优C 0.35P 0.95可以直接使用二级合格C 0.50P 0.80可用于短期预测三级勉强C 0.65P 0.70需要残差修正或换窗口四级不合格C 0.65P 0.70不建议继续使用如果检验结果落到三级以下优先去调整建模窗口长度而不是直接上残差修正。窗口问题不解决任何修正都是在错误基础上做补偿。4. 电力负荷预测的历史窗口选择滚动评估比拍脑袋定n更可靠4.1 输入长度n和预测步长7~14天窗口的取舍GM(1,1) 对输入窗口长度很敏感。同样的数据把 n 从 5 拉到 20预测值可能差出好几个百分点。原因在于累加序列会对趋势做放大窗口越长历史早期信息对指数拟合的影响越大模型越像是在平均历史而不是跟踪最近的走势。日粒度负荷我一般从 4 扫到 15周粒度负荷从 6 扫到 12秒级或小时级曲线不建议直接使用灰色模型先聚合成日尺度。窗口默认值常被写成 7理由是对齐一个完整周但灰色模型本身不感知周期。如果负荷呈现周一到周五高、周末低的双峰形态一条指数曲线拟合整个窗口预测点落在周末后的周一时前一个窗口末端是周日模型的趋势外推与实际方向相反。常见做法是分工作日和周末分别建模各用各的窗口。预测步长也要克制一次外推 1~3 步误差可控超过 5 步后指数外推会快速发散不如每收到一个真实值就重估一次参数。各数据粒度的大致参数范围可以参考下表数据粒度常用候选窗口 n推荐预测步长日峰值 / 日电量4~151~3周电量6~121~296点日内曲线不建议直接建模先聚合——4.2 用滚动窗口自动挑选nrolling_mape函数的3个关键边界与其凭经验定 n不如在已有历史上做一次滚动评估。做法是对每个候选 n从索引 n 开始依次取长度为 n 的训练段用 GM(1,1) 预测后续pred_steps个点计算实际值和预测值的 MAPE最后取平均。代码如下def rolling_mape(history, candidateslist(range(4, 16)), pred_steps1): 用滚动预测挑选合适的建模窗口长度 n history 连续历史负荷序列长度需大于 max(candidates) pred_steps candidates 候选的窗口长度列表 返回最优 n 及其平均 MAPE best_n None best_mape np.inf for n in candidates: mape_list [] for t in range(n, len(history) - pred_steps 1): train history[t - n:t] _, pred, _, _ gm11(train, pred_stepspred_steps) actual history[t:t pred_steps] mape_list.append(np.mean(np.abs(pred - actual) / actual)) cur_mape np.mean(mape_list) if cur_mape best_mape: best_mape cur_mape best_n n return best_n, best_mape三个边界条件决定了这段代码是否可靠。第一history的总长度必须大于max(candidates) pred_steps否则最内层循环遍历不到足够多的验证点。第二循环从t n开始保证第一个训练窗口恰好有 n 个点如果把range从 0 开始前面的若干步会用不完整的窗口建模型MAPE 会被低估。第三actual history[t:t pred_steps]与训练数据不重叠这是滚动验证的基本原则。实际使用时把候选 n 缩到 4 到 min(15, len(history)//2) 之间避免候选窗口超过历史一半。4.3 滚动更新而不是重复用旧模型每步都重新估计参数选定 n 之后生产中的预测方式应该是滚动更新。今天拿到真实负荷就把新值并入窗口尾部删掉窗口头部最旧的一个点重新调用gm11。代码上是这样n 7 load_history [156, 162, 158, 170, 168, 176, 183] # 模拟后续 5 天的滚动预测 for step in range(5): _, pred, a, b gm11(load_history[-n:], pred_steps1) next_load pred[0] # 有真实值就替换下一行没有真实值只能用预测值续算 load_history.append(next_load) load_history.pop(0) print(第 %d 次预测: %.2f MW, a%.4f % (step 1, next_load, a))注意load_history[-n:]每次取最近 n 个点窗口滚动时pop(0)会把最旧的点挤出。用预测值续算可以做完一个长序列的演示但生产环境切不可这样预测误差会逐日累积几天后预测值就和真实负荷脱节。应该由采集系统把真实值写入历史再触发下一次建模。4.4 强季节性和节假日冲击灰色模型撑不住的场景灰色模型对平滑趋势的拟合能力强对脉冲和相位变化几乎没有建模能力。夏季空调负荷脉冲、长假前后工厂停工复工会造成原始序列出现明显的突升突降。此时级比检验会大面积越界预测结果可能高出实际 20% 以上。我的做法是先做负荷结构分析把序列按工作日、周末、节假日分开每一类单独用灰色模型或者把原始负荷转成环比变化率对变化率序列建灰色模型再把变化率累乘还原成负荷。前一种做法保留原始量纲后一种更适合趋势频繁反转的场景。5. 残差修正与结果验证把电力负荷预测误差再压一档5.1 对残差序列再建GM(1,1)尾部修正的叠加方法当后验差检验达到三级但还达不到二级时残差修正是最直接的手段。把原始序列与拟合值的差当作一条新序列对它再跑一次 GM(1,1)用得到的残差预测值去修正原预测。由于残差序列包含正负值和 0建模前必须先平移。实现如下def residual_corrected(seq, pred_steps3): fit, pred, _, _ gm11(seq, pred_stepspred_steps) resid np.asarray(seq[:len(fit)]) - fit if len(resid) 6: return pred # 残差序列含负数时先平移 offset -np.min(resid) 1.0 resid_pos resid offset _, r_pred, _, _ gm11(resid_pos[-6:], pred_stepspred_steps) r_pred r_pred - offset return pred r_pred关键在resid_pos[-6:]这一段只用残差序列最近 6 个点做拟合而不是全量残差。残差本身的周期性很弱取太多点会把早期的旧模式强行外推叠加后反而破坏主模型的预测。修正后的结果要重新算后验差检验如果 C 没有降到二级以内说明残差中没有可建模的趋势这个修正不值得上线。5.2 先用持久性模型做基准判断灰色预测是否真的提升了在做残差修正之前先回答一个更基本的问题灰色模型比什么都别预测强多少。负荷预测里最常用的基线是持久性模型直接规定未来几小时的负荷等于当前值。日尺度上就是x̂(tm) x(t)。基线 MAPE 的计算代码很短def persistence_mape(history, pred_steps1): errs [] for t in range(len(history) - pred_steps): errs.append(abs(history[t] - history[t pred_steps]) / history[t pred_steps]) return np.mean(errs)把gm11的滚动 MAPE、persistence_mape的基线 MAPE、残差修正后的滚动 MAPE 放在同一张表里看。只有当灰色模型明显优于持久性模型并且残差修正还能再压低一个百分点以上这个方案才有工程意义。很多时候灰色模型拟合阶段误差很小滚动 MAPE 却比基线还高原因就是窗口选择不合理模型在用一个过时的趋势外推。5.3 验证路径三组输出一起看再决定是否启用上线前我会同时输出四样东西候选窗口 n、后验差比值 C、小误差概率 P、滚动 MAPE 与持久性基线 MAPE 的差值。判断顺序是先看 C 和 P 是否达到二级以上达不到就换窗口或清洗数据再看滚动 MAPE 是否比基线低低不了说明模型没有带来增量信息最后看残差修正是否把 MAPE 压低至少 1 个百分点没有就维持原模型。窗口长度用rolling_mape自动挑选生产预测时用 4.3 的滚动更新流程每收到一个真实值重估一次参数。先跑rolling_mape输出候选 n 和 MAPE把持久性基线和残差修正结果对齐到同一时间区间再决定生产配置里放哪组参数。本文还有配套的精品资源点击获取