1. 核心算法：GM(1,1) 灰色预测模型究竟怎么工作的？

代码做了什么：

grey_prediction 函数实现了 GM(1,1) 算法。它接收一串历史负荷数据，通过内部的矩阵运算，输出未来指定步长（默认2小时）的预测数据。

深度解读（为什么这么写 & 原理剖析）：

• 累加生成序列 (1-AGO)： 现实中的原始数据往往是杂乱无章的。代码中 x1 = np.cumsum(x0) 这一步是灰色系统的核心思想——“累加”。通过将前 个数据累加起来，原本无规律的序列会呈现出明显的指数增长规律。

• 构建数据矩阵与最小二乘法： 算法的本质是拟合一个微分方程。代码通过构建矩阵 $B$ 和向量 ，并使用最小二乘法求解参数（发展系数）和（灰色作用量）。这里的核心数学表达使用了线性代数中的经典公式：

  

  对应代码：u = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)。

• 累减还原 (IAGO)： 预测出累加序列的值后，必须使用 x1_hat_k_plus_1 - x1_hat_k 进行“差分”还原，才能得到我们真正需要的原始负荷预测值。

• 优势与适用场景： 为什么这里用灰色预测而不是火热的深度学习（如 LSTM）？因为 GM(1,1) 极度契合小样本预测（只需4个以上数据即可），且计算资源消耗极低。对于这种只有单日 24 小时数据的场景，灰色预测比神经网络更稳健。

• 异常处理的艺术： 注意代码中包含了 try...except np.linalg.LinAlgError:。在实际开发中，数据可能会导致矩阵不可逆。加上这层保护，程序就不会直接崩溃，提升了代码的健壮性。

2. 数据可视化：如何画出优雅且连续的预测曲线？

代码做了什么：

see_load 函数使用 matplotlib 将实际负荷和预测负荷绘制在同一张折线图上。

深度解读（避坑指南 & 技巧分享）：

• 尺寸对齐的“潜规则”： 在用 matplotlib 绘图时，很多新手经常会遇到 ValueError: All arrays must be of the same length 或者维度索引报错。这是因为 X 轴坐标和 Y 轴数据的长度不一致。

• 曲线平滑拼接技巧： 如果只画预测点，图上的实际曲线和预测曲线会断开，视觉效果很差。代码中这行 y_pred = [daily_load[-1]] + predicted_load 是一个非常实用的高阶技巧。它强行把实际负荷的最后一个点作为预测曲线的起点，配合 x_pred 的坐标延续，完美实现了两条曲线的无缝拼接。

• 中文字体配置： plt.rcParams['font.sans-serif'] = ["Fangsong"] 这一句直接解决了图表中文显示为方块的痛点，是写 Python 数据分析脚本的标配起手式。

3. 业务逻辑与代码优化：峰谷平电价计算

代码做了什么：

get_time_type 划分峰谷平时段，classify_load 统计各时段总电量，classify_load_bill 计算最终账单。

深度解读（设计思路 & 进阶替代方案）：

• 逻辑拆分的模块化思维： 这段代码没有把所有东西塞进一个函数，而是按照“打标签 -> 算电量 -> 算钱”的逻辑拆分成了三个独立函数。这使得代码具有极高的可读性和复用性。

• 替代方案对比： 目前代码使用 for 循环搭配繁琐的 if-elif 来判断时间区间。对于新手来说，这非常直观易懂。

  • 进阶提示： 但如果在工业级的大数据量场景中，这种纯 Python 循环效率较低。后续可以引导读者思考如何使用 Pandas 库的 pd.cut() 或者 numpy.select() 进行向量化操作，这会让代码行数锐减且运行速度成倍提升。

4. 程序的入口：Main 函数的流程控制

深度解读：

if __name__ == '__main__': 搭配 main() 函数是 Python 非常好的编程习惯。将耗时或阻塞进程的操作（如 plt.show() 弹窗）放在最后一步，确保了控制台的计算结果能第一时间打印出来，提升了用户体验。

详细代码如下

import matplotlib.pyplot as plt
import numpy as np

# 1. 设置全局字体为仿宋（解决中文显示问题）
plt.rcParams['font.sans-serif'] = ["Fangsong"]
# 可选：解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False


def grey_prediction(data, predict_steps=2):
    """
    使用 GM(1,1) 灰色预测模型预测后续数据
    :param data: 原始数据列表 (list or numpy array)
    :param predict_steps: 需要预测的步数
    :return: 预测的数据列表
    """
    x0 = np.array(data)
    n = len(x0)

    # 1. 累加生成序列 (1-AGO)
    x1 = np.cumsum(x0)

    # 2. 计算紧邻均值生成序列 (z1)
    z1 = 0.5 * (x1[:-1] + x1[1:])

    # 3. 构造数据矩阵 B 和数据向量 Y
    B = np.zeros((n - 1, 2))
    B[:, 0] = -z1
    B[:, 1] = 1
    Y = x0[1:].reshape((n - 1, 1))

    # 4. 计算参数 a 和 b (最小二乘法)
    try:
        u = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
    except np.linalg.LinAlgError:
        print("矩阵不可逆，无法进行灰色预测计算。")
        return []

    a, b = u[0][0], u[1][0]

    # 5. 构建预测模型并进行预测
    predictions = []
    for k in range(n, n + predict_steps):
        # 预测累加序列的值
        x1_hat_k = (x0[0] - b / a) * np.exp(-a * (k - 1)) + b / a
        x1_hat_k_plus_1 = (x0[0] - b / a) * np.exp(-a * k) + b / a

        # 6. 累减还原 (IAGO) 得到原始序列的预测值
        x0_hat_k_plus_1 = x1_hat_k_plus_1 - x1_hat_k
        predictions.append(round(x0_hat_k_plus_1, 2))  # 保留两位小数

    return predictions


def see_load(daily_load, predicted_load=None):
    fig, ax = plt.subplots()
    """fig表示生成的一系列绘图构成的图形。ax用来定义图像文件"""

    # 绘制实际负荷
    ax.plot(range(len(daily_load)), daily_load, label="实际负荷")

    # 如果有预测数据，则绘制预测曲线
    if predicted_load:
        # X轴坐标延续实际负荷的长度
        x_pred = range(len(daily_load) - 1, len(daily_load) + len(predicted_load))
        # Y轴数据：拼接最后一个实际点和所有预测点，使线条连续
        y_pred = [daily_load[-1]] + predicted_load
        # 绘制红色虚线作为预测曲线
        ax.plot(x_pred, y_pred, 'r--', label="预测负荷 (后2小时)")
        ax.legend()  # 显示图例

    ax.set_title("电厂的24小时负荷及预测情况", fontsize=18)
    ax.set_xlabel('时段/小时', fontsize=18)  # 设置X轴的命名与字体大小
    ax.set_ylabel('逐时负荷/MW', fontsize=18)  # 设置y轴的命名与字体大小
    plt.show()  # 展示图像


def get_time_type(daily_load):
    time_type = []
    for i in range(len(daily_load)):
        if (i >= 8 and i < 11) or (i >= 18 and i < 23):
            time_type.append("peak")
        elif (i >= 7 and i < 8) or (i >= 11 and i <= 18) or (i == 23):
            time_type.append("flat")
        else:
            time_type.append("valley")
    return time_type


def classify_load(daily_load):
    time_types_list = get_time_type(daily_load)
    peak_all_load, flat_all_load, valley_all_load = 0, 0, 0
    for i in range(len(time_types_list)):
        if time_types_list[i] == "peak":
            peak_all_load += daily_load[i]
        elif time_types_list[i] == "flat":
            flat_all_load += daily_load[i]
        else:
            valley_all_load += daily_load[i]
    print(f"峰时段用电总量是: {peak_all_load}（MW）")
    print(f"平时段用电总量是: {flat_all_load}（MW）")
    print(f"谷时段用电总量是: {valley_all_load}（MW）")
    return peak_all_load, flat_all_load, valley_all_load


def classify_load_bill(peak_load, flat_load, valley_load, time_prices):
    """计算各个时段的电费"""
    peak_load_bill = peak_load * time_prices["peak"]  # 峰时段电费
    flat_load_bill = flat_load * time_prices["flat"]  # 平时段电费
    valley_load_bill = valley_load * time_prices["valley"]  # 谷时段电费
    load_bill_all = peak_load_bill + flat_load_bill + valley_load_bill
    print("===========================================")
    print(f"峰时段的电费是: {peak_load_bill:.2f}元")
    print(f"平时段的电费是: {flat_load_bill:.2f}元")
    print(f"谷时段的电费是: {valley_load_bill:.2f}元")
    print(f"当天的电费是: {load_bill_all:.2f}元")


def load_max(daily_load):
    """定位当日最大负荷及对应时刻"""
    idx = daily_load.index(max(daily_load))
    Type_time = get_time_type(daily_load)
    print("===========================================")
    print(f"当日最大负荷为{daily_load[idx]}MW, 对应第{idx}时段, 属于{Type_time[idx]}时刻")


def calc_electric_cost(daily_load, time_prices):
    """计算各时段总用电量"""
    peak_load, flat_load, valley_load = classify_load(daily_load)
    classify_load_bill(peak_load, flat_load, valley_load, time_prices)


def main():
    """假设每个时段消耗的负荷数量"""
    daily_load = [50, 48, 45, 42, 40, 46, 55, 65, 80, 90, 85, 75, 70, 68, 66, 69, 72, 85, 95, 100, 98, 92, 80, 60]
    time_prices = {"peak": 0.85, "flat": 0.55, "valley": 0.25}  # 存入分段信息的字典

    # 1. 执行灰色预测
    print("===========================================")
    print("正在计算 GM(1,1) 灰色预测模型...")
    predicted_load = grey_prediction(daily_load, predict_steps=2)
    print(f"预测未来2小时的负荷数据为: {predicted_load[0]}MW, {predicted_load[1]}MW")
    print("===========================================")

    # 2. 核心分析模块
    calc_electric_cost(daily_load, time_prices)  # 定义计算电费模块
    load_max(daily_load)

    # 3. 数据可视化 (将其放在最后，防止阻塞命令行打印)
    see_load(daily_load, predicted_load)


if __name__ == '__main__':
    main()