电工杯 2021 B 范文：短期负荷预测（56 天×24 h 合成数据 · 傅里叶+温度二次特征回归 + 残差 AR(1) 滚动修正 + 三模型递进评价）

一、摘要

针对电工杯 2021 年 B 题「负荷预测」，本文以一套 56 天×24 h（共 1344 点）的区域电网逐时负荷与气温合成数据为对象，建立并逐级验证了「朴素基线 → 气象特征回归 → 残差时序修正」的三层递进预测体系。数据合成遵循真实电网负荷的三重结构：日内双峰形态（早峰 10:30、晚峰 19 时）、温度二次敏感（低温采暖与高温制冷双重抬升）与周历效应（周六、周日负荷系数分别降为 0.88 与 0.85）。预测模型方面：模型一采用「星期×小时」同期均值基线，测试集 MAPE 高达 16.62%——它只记得「形状」而看不见温度，暴露了纯历法基线在气温波动下的失效；模型二以日周期傅里叶对（sin/cos ωh、sin/cos 2ωh）刻画日内形状、以温度与制冷度²/采暖度²刻画气象弹性、以周末哑元刻画周历效应，共 9 维特征做多元线性回归（正规方程 + 列主元高斯消元纯标准库求解），MAPE 压缩至 4.19%、RMSE 30.30 MW，较基线改善 74.8%，且回归系数量级与物理直觉一一对应（周末哑元 −74.26 MW、采暖度² 3.30 > 制冷度² 1.10）；模型三针对回归残差 lag-1 自相关高达 0.648 的强正自相关结构，建立 AR(1) 滚动修正（ $\hat y_t=X_t\hat\beta+\varphi\,\hat e_{t-1}$ ， $\varphi=0.6483$ ），在提前 1 h 滚动场景下把 MAPE 进一步压到 2.99%（较回归再降 28.6%），残差 lag-1 自相关由 0.648 白化至 0.443。逐日检验显示最优模型全周 MAPE 稳定于 2.57%~3.43%，无单日失效。全文纯标准库实现、固定种子（SEED=20260904）可复现，并确定性补齐本题缺失的练习数据集 data/dgcup2021b.csv（1344 行），正文、配图、附录与真源四路数字一致。

二、问题重述

短期负荷预测是电网调度计划（机组组合、联络线计划、需求响应安排）的入口环节：预测偏高则机组空转浪费燃料，偏低则被迫高价购电甚至拉闸。赛题要求基于历史用电数据与气象信息建立时间序列模型，完成高精度预测并给出科学评价。其实质可拆为四个子问题：

结构识别：负荷序列由哪些成分构成？日周期形状、气象弹性、周历节律各占多少？
模型构建：如何把领域知识翻译成可估计的数学结构，使参数可解释、可复现？
精度提升：在日尺度模型之上，如何利用「近期误差蕴含近期信息」的时序相关性做超短期修正？
科学评价：MAPE、RMSE、MAE 各自度量什么？残差是否还残留可利用结构？模型在整周内是否稳定？

本文按此四问组织：先合成并解剖数据（第四、五节），再逐级建立三层模型并用统一测试窗对比（第六节），最后做残差诊断与逐日稳健性检验（第七、八节）。

三、模型假设与符号

H1（合成数据）：官方数据未公开，本文按真实电网负荷结构合成 56 天逐时数据并固定种子，方法链对真实数据同样适用。
H2（准静态气象弹性）：负荷对气温的响应在小时尺度上即时完成，无热惯性延迟；弹性用二次型（制冷度 $C=\max(T-18,0)$ 、采暖度 $H=\max(14-T,0)$ ）表达。
H3（周历平稳）：周内节律在样本窗内稳定，周六/周日效应可用常数哑元刻画。
H4（加性结构）：负荷 = 形状项 + 气象项 + 周历项 + 随机项，各成分可加且随机项均值为零。
H5（滚动信息可用）：模型三工作在提前 1 h 滚动场景，即预测第 $t$ 小时时， $t-1$ 小时的真实负荷已可观测。

主要符号：第 $d$ 天第 $h$ 小时气温 $T_{d,h}$ （°C）、负荷 $y_{d,h}$ （MW）；特征向量 $x=(1,\sin\omega h,\cos\omega h,\sin 2\omega h,\cos 2\omega h,D_{\mathrm{wk}},T,C^2,H^2)$ ， $\omega=2\pi/24$ ；回归系数 $\hat\beta$ ；残差 $e=y-X\hat\beta$ ；AR(1) 系数 $\varphi=\sum_t e_{t-1}e_t/\sum_t e_{t-1}^2$ ；评价指标 MAPE $=\frac1n\sum|y-\hat y|/y\times100\%$ 、RMSE、MAE。

四、数据合成与负荷结构探索（问题一）

4.1 建模框架

图1 建模框架

图 1 给出全链路：数据合成 → 特征工程 → 正规方程回归 → 残差 AR 建模 → 滚动预测与评价。整条链路只用 Python 标准库（math/random），任何环境可逐位复现。

4.2 序列形态与测试窗

图2 最后两周逐时负荷

图 2 展示最后两周的逐时负荷：日内「早峰 + 更高晚峰」的双峰形态清晰可辨，工作日与周末之间可见约一成的台阶式下移。蓝色为训练期末段（第 43~~49 天），红色为测试周（第 50~~56 天，共 168 点）——所有模型的对比都在这同一窗口上进行，杜绝「用测试集调参」的数据泄漏。全样本温度范围 −0.7~~22.1 °C、负荷 365.7~~1122.4 MW（均值 542.0 MW），峰谷差近 3 倍，凸显高精度预测的调度价值。

4.3 温度弹性与周历节律

图3 负荷-温度散点

图4 平均日内曲线

图 3 的负荷—温度散点呈明显的「两翼上翘」：以 14~18 °C 为舒适区，低于 14 °C 时采暖负荷按温度缺口的平方抬升（低温翼更陡），高于 18 °C 时制冷负荷同样二次放大——这正是 H2 中二次型弹性设定的事实依据。图 4 按周几分组的平均日内曲线显示：晚峰（19 时）为主峰、早峰（10:30 附近）为次峰；周六、周日整体下移（系数 0.88 / 0.85，即约 12% 与 15%），且周末早峰更平缓——「晚峰为主、周末塌腰」是后续特征设计直接引用的结构先验。

五、三层递进预测模型（问题二）

5.1 模型一：同期均值基线

最朴素的可用基线是「历史同期条件均值」：对测试周每对（星期几 $w$ ，小时 $h$ ），取训练集中同组合的负荷均值 $\hat y=\bar y_{w,h}$ 。它保留了日内形状与周历节律，但完全不含气象输入——训练窗内气温从 6 °C 缓变到 15 °C 再回落，同期均值把不同温度状态的日子硬性平均，误差自然失控。测试集 MAPE 16.62%、RMSE 178.28 MW：这个数字本身就是结论的一部分——负荷预测的第一杠杆不是更复杂的算法，而是把气象信息写进特征。

5.2 模型二：傅里叶 + 温度二次 + 周末哑元的线性回归

按第四节的结构先验构造 9 维特征并做最小二乘：

$y = \beta_0+\beta_1\sin\omega h+\beta_2\cos\omega h+\beta_3\sin2\omega h+\beta_4\cos2\omega h+\beta_5D_{\mathrm{wk}}+\beta_6 T+\beta_7 C^2+\beta_8 H^2+\varepsilon$

其中一阶傅里叶对刻画「日单循环」、二阶对补充晚峰更尖的非对称形状； $C^2,H^2$ 即图 3 的两翼。求解用正规方程 $(X'X)\hat\beta=X'y$ ，列主元高斯消元实现（9 维方程组毫秒级完成）。表 1 给出全部系数及其解读。

表 1 回归系数与物理解读

特征	系数	解读
常数项	451.81	舒适区基础负荷（MW），与图 4 谷值吻合
sin(ωh) / cos(ωh)	−31.79 / −20.01	日单循环基波
sin(2ωh) / cos(2ωh)	−36.44 / −11.40	双峰非对称修正（晚峰更尖）
周末哑元	−74.26	周末整体减载，量级与图 4 台阶一致
温度 T	2.33	舒适区内的弱线性漂移
制冷度² C²	1.10	高温制冷翼弹性
采暖度² H²	3.30	低温采暖翼弹性约为制冷翼 3 倍

每个系数量级都能在探索性图上找到对应物——这不是巧合，而是「结构先验 → 特征设计 → 参数可解释」闭环的直接证据。

5.3 模型三：残差 AR(1) 滚动修正

回归残差 $e_t$ 的 lag-1 自相关高达 0.648：今天的误差与上一小时的误差强烈同号，意味着回归面「追不上」负荷的日内漂移（如天气实况与合成气候的逐时偏离）。据此建立 AR(1) 修正：

$\hat y_t^{(3)} = X_t\hat\beta + \varphi\,\hat e_{t-1},\qquad \varphi=\frac{\sum_t e_{t-1}e_t}{\sum_t e_{t-1}^2}=0.6483$

$\varphi$ 在训练残差上估计；测试时逐时滚动——每预测完一小时，立即用该小时的真实负荷刷新残差估计 $\hat e_t$ 。该模型工作在提前 1 h 超短期场景（与日前计划的模型一/二口径不同，正文明确区分），其价值在于把「刚刚发生的偏差」外推一步。

六、预测精度与残差诊断（问题三）

图5 测试周前3天预测对比

图6 三模型MAPE对比

图 5 叠画测试周前 3 天的实际负荷与两条预测曲线：同期均值（红）系统性偏离——它按「历史平均温度的日子」出牌，而测试周偏冷，采暖负荷被整体低估；回归（蓝）则紧贴实际曲线，峰谷位置与幅度均对齐。表 2 汇总统一测试窗上的三模型指标。

表 2 三模型测试集指标对比（168 点）

模型	MAPE	RMSE (MW)	MAE (MW)	场景
一：同期均值	16.62%	178.28	127.45	日前计划
二：线性回归	4.19%	30.30	25.88	日前计划
三：回归+AR(1)	2.99%	22.10	18.39	提前 1 h 滚动

递进关系一目了然：特征工程把 MAPE 从 16.62% 砍到 4.19%（改善 74.8%），残差时序修正再压 28.6% 到 2.99%；RMSE 的改善（178.28 → 30.30 → 22.10 MW）说明大误差样本被系统性消灭而非仅均值对齐。

图7 残差自相关对比

图 7 从诊断角度复核模型三的机制：回归残差的 |自相关| 在 lag 1~3 显著非零（0.65/0.45/0.30 一带），呈典型 AR 结构；AR(1) 修正后 lag-1 降至 0.443、高阶滞后全面回落——残差被显著白化，模型三不是「碰巧更准」，而是吃掉了残差中真实存在的可预测成分。

七、逐日稳健性检验

图8 逐日MAPE

表 3 最优模型（模型三）逐日 MAPE

测试日	D50	D51	D52	D53	D54	D55	D56
MAPE (%)	3.43	3.05	3.19	3.28	2.57	2.64	2.76

七个测试日 MAPE 全部落于 2.57%~3.43%，极差不足 0.9 个百分点，无单日失效或周末特异性劣化——模型在整周尺度上稳定，可支撑调度计划的例行编制。

八、模型检验、优缺点与拓展

检验：① 一致性——正文全部数字由真源脚本一次运行产出，附录独立执行结果逐位一致（四路一致）；② 无泄漏——三层模型共享同一训练/测试切分，AR(1) 的 $\varphi$ 仅在训练残差上估计；③ 可解释性——表 1 每个系数均可在图 3/图 4 上找到物理对应。

主要优点：① 三层递进把「基线—特征—时序」的每一层增益干净分离（74.8% 与 28.6%），归因清晰，便于向评审与调度员解释；② 纯标准库、固定种子，零环境依赖；③ 残差诊断闭环（先测自相关、再建 AR、后验白化），方法论完整。

局限与改进方向：① 线性 + 二次型的气象响应在极端温度日（超出训练温度域）存在外推风险，可引入样条或分段回归；② H2 忽略建筑热惯性，真实系统宜加入温度滞后项 $T_{t-1},T_{t-2}$ ；③ 周历效应用常数哑元，节假日（调休、重大活动）需单独日历特征；④ 进一步可对比 BP 神经网络 / SVM 等非线性基学习器，或对多气象站点、负荷场景做组合预测与概率区间预报。

九、结论

本文围绕电工杯 2021 B 题建立了「结构探索 → 特征工程 → 残差修正」的短期负荷预测完整方案，在 56 天×24 h 合成数据（1344 点，固定种子）上的核心结论：① 负荷由日内双峰、温度二次弹性与周历节律三重结构叠加，负荷—温度散点两翼上翘、周末下移 12%15%；② 纯历法同期均值基线 MAPE 高达 16.62%，证明气象特征是负荷预测的第一杠杆；③ 傅里叶 + 温度二次 + 周末哑元的 9 维线性回归把 MAPE 压到 4.19%（RMSE 30.30 MW），系数全部物理可解释；④ 针对残差 lag-1 自相关 0.648 的强 AR 结构，AR(1) 滚动修正进一步把提前 1 h 预测的 MAPE 压至 2.99%（再降 28.6%），残差显著白化，逐日 MAPE 稳定于 2.57%3.43%。方法论层面，「每一层增益单独计量」的递进评价框架可直接迁移到真实电网的预测工程与竞赛写作中。

参考文献

[1] 王成山, 王赛一. 电力系统负荷预测方法与实践[M]. 北京: 中国电力出版社, 2018.

[2] 康重庆, 夏清, 刘梅. 电力系统负荷预测[M]. 北京: 高等教育出版社, 2007.

[3] Box G E P, Jenkins G M, Reinsel G C, et al. Time Series Analysis: Forecasting and Control[M]. 5th ed. Hoboken: Wiley, 2015.

[4] Hong T, Fan S. Probabilistic electric load forecasting: A tutorial review[J]. International Journal of Forecasting, 2016, 32(3): 914–938.

[5] 中国电机工程学会. 全国大学生电工数学建模竞赛历年赛题汇编[M]. 北京: 中国电力出版社, 2020.

附录：核心 Python 实现

以下代码为真源 tools/gen_dgcup2021b.py 的等价核心摘录（固定种子、纯标准库），可在本文件所在目录直接运行，输出与正文全部数字一致。

# -*- coding: utf-8 -*-
"""短期负荷预测：同期均值 / 傅里叶+温度回归 / 残差 AR(1) 修正核心。"""
import math
import random

SEED = 20260904
NDAYS = 56
NH = 24
NTRAIN = 49                       # 前 49 天训练，后 7 天测试


def gen_data():
    rng = random.Random(SEED)
    rows = []
    for d in range(NDAYS):
        wd = (d + 3) % 7
        day_mean = 6.0 + 9.0 * math.sin(math.pi * d / 55.0)
        for h in range(NH):
            t = day_mean + 5.5 * math.sin(2 * math.pi * (h - 9) / 24.0) \
                + rng.gauss(0.0, 0.8)
            base = 430.0 \
                + 120.0 * math.exp(-((h - 10.5) / 2.4) ** 2) \
                + 150.0 * math.exp(-((h - 19.0) / 2.6) ** 2)
            cdd = max(t - 18.0, 0.0)
            hdd = max(14.0 - t, 0.0)
            week = 1.0 if wd <= 4 else (0.88 if wd == 5 else 0.85)
            load = (base + 6.5 * cdd * cdd + 3.2 * hdd * hdd) * week \
                + rng.gauss(0.0, 7.0)
            rows.append(dict(d=d, h=h, wd=wd,
                             temp=round(t, 2), load=round(load, 2)))
    return rows


def solve_lin(A, b):
    """列主元高斯消元解 Ax=b。"""
    n = len(A)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for c in range(n):
        piv = max(range(c, n), key=lambda r: abs(M[r][c]))
        if abs(M[piv][c]) < 1e-12:
            raise ValueError("singular")
        M[c], M[piv] = M[piv], M[c]
        for r in range(c + 1, n):
            f = M[r][c] / M[c][c]
            for k in range(c, n + 1):
                M[r][k] -= f * M[c][k]
    x = [0.0] * n
    for r in range(n - 1, -1, -1):
        s = M[r][n] - sum(M[r][k] * x[k] for k in range(r + 1, n))
        x[r] = s / M[r][r]
    return x


def ols(X, y):
    p = len(X[0])
    XtX = [[sum(X[i][a] * X[i][bb] for i in range(len(X)))
            for bb in range(p)] for a in range(p)]
    Xty = [sum(X[i][a] * y[i] for i in range(len(X))) for a in range(p)]
    return solve_lin(XtX, Xty)


def features(temp, h, wd):
    hr = 2 * math.pi * h / 24.0
    weekend = 1.0 if wd >= 5 else 0.0
    cdd = max(temp - 18.0, 0.0)
    hdd = max(14.0 - temp, 0.0)
    return [1.0, math.sin(hr), math.cos(hr),
            math.sin(2 * hr), math.cos(2 * hr),
            weekend, temp, cdd * cdd, hdd * hdd]


FEANAMES = ["常数", "sin(w*h)", "cos(w*h)", "sin(2w*h)", "cos(2w*h)",
            "周末哑元", "温度", "制冷度2", "采暖度2"]


def metrics(y, yhat):
    n = len(y)
    mape = sum(abs(y[i] - yhat[i]) / y[i] for i in range(n)) / n * 100.0
    rmse = math.sqrt(sum((y[i] - yhat[i]) ** 2 for i in range(n)) / n)
    mae = sum(abs(y[i] - yhat[i]) for i in range(n)) / n
    return mape, rmse, mae


def autocorr(e, lag):
    n = len(e)
    m = sum(e) / n
    num = sum((e[i] - m) * (e[i + lag] - m) for i in range(n - lag))
    den = sum((v - m) ** 2 for v in e)
    return num / den if den > 0 else 0.0


def run():
    rows = gen_data()
    tr = rows[:NTRAIN * NH]
    te = rows[NTRAIN * NH:]
    temps = [r["temp"] for r in rows]
    loads = [r["load"] for r in rows]

    print("== 数据 ==")
    print("样本 %d 天×%d 时=%d 点 训练 %d 点 测试 %d 点"
          % (NDAYS, NH, len(rows), len(tr), len(te)))
    print("温度 %.1f~%.1f C 负荷 %.1f~%.1f MW 均值 %.1f MW"
          % (min(temps), max(temps), min(loads), max(loads),
             sum(loads) / len(loads)))

    cond = {}
    for r in tr:
        cond.setdefault((r["wd"], r["h"]), []).append(r["load"])
    condm = {k: sum(v) / len(v) for k, v in cond.items()}
    pred1 = [condm[(r["wd"], r["h"])] for r in te]
    m1 = metrics([r["load"] for r in te], pred1)
    print("== 模型一 同期均值基线 ==")
    print("MAPE %.2f%% RMSE %.2f MW MAE %.2f MW" % m1)

    Xt = [features(r["temp"], r["h"], r["wd"]) for r in tr]
    yt = [r["load"] for r in tr]
    beta = ols(Xt, yt)
    Xe = [features(r["temp"], r["h"], r["wd"]) for r in te]
    ye = [r["load"] for r in te]
    pred2 = [sum(b * x for b, x in zip(beta, row)) for row in Xe]
    m2 = metrics(ye, pred2)
    print("== 模型二 多元线性回归 ==")
    print("MAPE %.2f%% RMSE %.2f MW MAE %.2f MW" % m2)
    print("回归系数:", " ".join("%s=%.2f" % (n, b)
                                for n, b in zip(FEANAMES, beta)))

    fit = [sum(b * x for b, x in zip(beta, row)) for row in Xt]
    etr = [yt[i] - fit[i] for i in range(len(yt))]
    phi = sum(etr[i - 1] * etr[i] for i in range(1, len(etr))) \
        / sum(v * v for v in etr[:-1])
    print("残差 AR(1): phi=%.4f 训练残差 lag-1 自相关 %.4f"
          % (phi, autocorr(etr, 1)))

    pred3 = []
    e_prev = etr[-1]
    for i, r in enumerate(te):
        base_p = pred2[i]
        pred3.append(base_p + phi * e_prev)
        e_prev = ye[i] - base_p
    m3 = metrics(ye, pred3)
    print("== 模型三 回归+AR(1) 滚动修正 ==")
    print("MAPE %.2f%% RMSE %.2f MW MAE %.2f MW" % m3)
    imp_b = (m1[0] - m2[0]) / m1[0] * 100.0
    imp_a = (m2[0] - m3[0]) / m2[0] * 100.0
    print("回归较基线 MAPE 降低 %.1f%% AR 修正再降 %.1f%%"
          % (imp_b, imp_a))
    print("AR 后残差 lag-1 自相关 %.4f"
          % autocorr([ye[i] - pred3[i] for i in range(len(ye))], 1))

    daily = []
    for d0 in range(NTRAIN, NDAYS):
        seg = [(i, r) for i, r in enumerate(te) if r["d"] == d0]
        mp = sum(abs(ye[i] - pred3[i]) / ye[i] for i, _ in seg) / len(seg) \
            * 100.0
        daily.append(mp)
    print("逐日 MAPE(模型三):", " ".join("%.2f" % v for v in daily))
    return dict(rows=rows, tr=tr, te=te, beta=beta, phi=phi,
                m1=m1, m2=m2, m3=m3, pred1=pred1, pred2=pred2,
                pred3=pred3, daily=daily, etr=etr)


if __name__ == "__main__":
    run()