DataWhale市场博弈和价格预测

news2024/11/16 4:37:26

文章目录

  • 市场博弈和价格预测
    • 任务
    • 数据
    • 出清价格的形成步骤
    • 代码运行
  • Task1：跑通baseline
    • 使用ABM估计市场出清价格
    • 转换耗煤量为机组报价
    • 保存结果为submit.csv
  • Task2：时间序列挖掘+ABM构建学习
    • 数据基本信息
    • 统计指标分析
    • 数据分布
    • 分时统计特征
      • 不同小时的总负荷和电价
    • 负电价与高电价形成原因分析
      • 负电价
      • 高电价
    • 双变量分析
      • 相关系数
    • 总结
    • ABM构建学习与实战
      • 无法避免的贫富差距——玻尔兹曼财富模型
        • 构建环境
        • 定义agent
  • Task3：进阶上分方向讨论
    • 时间序列挖掘
      • 特征工程
        • 其他时序特征
      • 模型训练
      • 模型融合

市场博弈和价格预测

比赛链接：http://competition.sais.com.cn/competitionDetail/532232/format

任务

本次比赛提供某市场的出清价格以及市场参与主体相关参数信息，期待选手使用Agent-based model根据历史市场需求、外部环境以及主体基本信息对主体博弈行为进行刻画，预测未来情景下的市场出清价格。
本质是一个回归问题，需要预测2023年7月1日到2024年4月18日每15分钟的市场出清价格，虽然可以纯粹依靠时间序列模型完成，但比赛方也提到期待选手使用ABM模型建模获取市场出清价格，因此可以在两方面进行尝试上分（ps：如果感觉ABM较难，可以专注于时间序列挖掘）。

数据

electricity price.csv：电力市场的市场出清价格，市场需求等信息。训练集范围为2021年12月1日到2023年7月1日，共计55392个点；测试集范围为2023年7月1日到2024年4月18日，共计28228个点

  1. Day/Time：交易时间，中国电力现货市场15分钟结算一次，一天共96个交易点
  2. demand：区域内电力总负荷（总需求），单位为MW
  3. clearing price (CNY/MWh)：市场出清电价，单位为元/MW·h

unit.csv：存放市场供给者（各发电机组）的参数信息。机组数据包含549个不同的火电机组

  1. unit ID：每个机组唯一的ID
  2. Capacity（MW）：机组的额定容量（额定功率），越高机组的发电能力越强
  3. utilization hour (h) ：电厂的年平均运行小时数，需要注意多个机组可能共同属于一个电厂，有相同的值
  4. coal consumption (g coal/KWh)：每发一度电需要耗费多少煤炭，为成本参数
  5. power consumption rate：电厂单位时间内耗电量与发电量的百分比，例如单位时间耗电量为500度电，发电量为10000度电，利用率就是500/10000=5%。

出清价格的形成步骤

  1. 所有发电机组申报自己卖出的电价和电量
  2. 市场根据机组报价，从低到高排序，依次从低价开始成交
  3. 当成交的容量和大于等于总需求时，达到市场出清（供需平衡），这时候最后一个达成交易的机组报价为市场出清价格
    假如市场总需求为3000MW，四个机组依次报价报量如下：
  • 机组1：报价150元，容量500MW
  • 机组2：报价200元，容量500MW
  • 机组3：报价250元，容量2000MW
  • 机组4：报价400元，容量2500MW

按价格从低到高，机组1先成交，随后是机组2、机组3，此时容量和等于总需求，不再有额外的电力需求，市场达到出清状态，市场出清价格为250元。
这里我们还可以发现一些现象：

  • 机组4由于过高的价格竞标失败，导致没有卖出任何电力。现实中机组4不得不停掉一些发电机。但煤电机组频繁关停会影响设备寿命和性能，还可能延误后续高价时段的正常发电，最终只能亏本卖出。
  • 机组1实际上能卖更多钱，如果他能预料到机组2和3报价为200和250，他就可以报价300元，保证最后一个出清的是他，而不是只能卖150元（但它无法报特别高的价，因为电力是公共物品，受政府管控，设置了价格上下限）

因此报价本身就是一个充满策略博弈的行为，就像在股市中总想做买入价格最低，卖出价格最高的那个人。

代码运行

以下baseline代码均在kaggle平台运行

Task1：跑通baseline

学习链接：https://linklearner.com/activity/12/2/2

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

这里我们做了几件事：

  1. 合并市场数据的日期和时间，便于利用pandas高效的时间序列处理能力（例如我们能直接用大小于号取一段时间序列）
  2. 处理市场数据以符合最终的提交格式，保存到本地便于后续直接使用

submit.csv提交格式

day,time,clearing price （CNY/MWh）
2024/4/1 , 0:15 , 352.334 
2024/4/1 , 0:30 , 355.536
  1. 设置市场数据的列顺序，并去除合并前的时间列

在进行进一步编写代码前，有效的预处理能大量减少后续迭代方案的耗时

# 读取市场数据
electricity_price = pd.read_csv('/kaggle/input/20240730a/electricity price.csv')
# 读取市场主体（各发电机组）数据
unit = pd.read_csv('/kaggle/input/20240730a/unit.csv')

"""
准备示例提交数据sample_submit
1. electricity_price["clearing price (CNY/MWh)"].isna()找到出清价格为缺失值的行，即要预测的目标
2. 去除demand列，符合最后的提交格式 
"""

sample_submit = electricity_price[electricity_price["clearing price (CNY/MWh)"].isna()].drop(columns="demand")
sample_submit.to_csv("sample_submit.csv", index=False)

# 将day和time列合并成timestamp列，便于提取时间戳特征
electricity_price["timestamp"] = pd.to_datetime(
    electricity_price["day"] + " " + electricity_price["time"].str.replace("24:00:00", "00:00"))

# 处理24:00:00的情况，即表示第二天的00:00:00
mask = electricity_price['timestamp'].dt.time == pd.Timestamp('00:00:00').time()

# 需要将这些行的日期部分加一天
electricity_price.loc[mask, 'timestamp'] += pd.Timedelta(days=1)

# 设置列的顺序，同时去除day和time列
electricity_price = electricity_price[["timestamp", "demand", "clearing price (CNY/MWh)"]]

我们可以查看一下市场数据，这里包含约55000个出清电价和电力负荷数据，市场每15分钟生成一次出清电价

  1. timestamp：时间戳
  2. demand：区域内电力总负荷（总需求），单位为MW
  3. clearing price (CNY/MWh)：市场出清电价，单位为元/MW·h

其中MW为功率单位，MW·h为能量单位，表示1MW的电器运行一小时产生的能量（例如1度电就是1KW·h）

electricity_price.head()  # 显示前5行数据
timestampdemandclearing price (CNY/MWh)
02021-12-01 00:15:0040334.18350.80
12021-12-01 00:30:0040523.15350.80
22021-12-01 00:45:0040374.74350.80
32021-12-01 01:00:0040111.55350.80
42021-12-01 01:15:0040067.50348.93

机组数据包含549个不同的火电机组

  1. unit ID：每个机组唯一的ID
  2. Capacity（MW）：机组的额定容量（额定功率），越高机组的发电能力越强
  3. utilization hour (h) ：电厂的年平均运行小时数，需要注意多个机组可能共同属于一个电厂，有相同的utilization hour (h)
  4. coal consumption (g coal/KWh)：每发一度电需要耗费多少煤炭，为成本参数
  5. power consumption rate：电厂单位时间内耗电量与发电量的百分比，例如单位时间耗电量为500度电，发电量为10000度电，利用率就是500/10000=5%。
unit.head()
unit IDCapacity（MW）utilization hour (h)coal consumption (g coal/KWh)power consumption rate (%)
01110.02069.12266.076.91
12160.05509.22292.706.91
23160.03562.79293.356.91
34160.05684.12284.886.91
45220.02231.35323.088.54

使用ABM估计市场出清价格

Agent-Based Modeling (ABM) 是一种模拟系统的复杂行为和动态的计算方法。它通过定义系统中的个体(称为“代理(Agent)”)及其相互作用规则来进行建模。每个代理具有独特的属性和行为。ABM 模型通过模拟代理之间的互动，观察宏观层面的模式和现象，帮助理解复杂系统的行为和演化。

在这里，市场的主体是各个发电机组，他们在电力现货市场中根据市场信息/其他机组信息进行报价，并以获得最大利润为目标。

这里我们用边际成本定价法作为baseline，即每个机组都以其自身每生产一度电的成本进行报价，而不考虑市场信息和其他机组信息。

类比来说就是只会以成本价卖出商品（事实上这就是完全竞争市场的长期均衡状态）

随后，按报价从低到高排序，依次接受报价，直到累计的功率超过了总需求，市场达到出清状态，并以最后一个满足的报价作为市场出清价格的估计。

如果感到无法理解，Task1链接文档下准备了一些背景知识，可以参考阅读。

sorted_unit = unit.sort_values("coal consumption (g coal/KWh)")  # 按照一度电的耗煤量（近似为边际成本）降序排序
sorted_unit.head()
unit IDCapacity（MW）utilization hour (h)coal consumption (g coal/KWh)power consumption rate (%)
11411560.03318.063.00.00
23823925.06486.073.05.33
15015130.03531.074.04.07
14814930.03531.074.03.61
14915030.03531.074.08.72
# 预先计算 sorted_unit 的累积和
sorted_unit['cumulative_capacity'] = sorted_unit['Capacity（MW）'].cumsum()

prices = []

# 找到最后一个满足总需求的机组报价
for demand in electricity_price["demand"]:
    price = sorted_unit[sorted_unit['cumulative_capacity'] >= demand]["coal consumption (g coal/KWh)"].iloc[0]
    prices.append(price)

print(len(prices))
prices[:5]


83520

[279.0, 279.0, 279.0, 279.0, 279.0]

转换耗煤量为机组报价

由于这里耗煤量单位是g,虽然和成本正相关但不能等同于报价（事实上应该考虑煤价和其他固定资产开销成本）。由于我们不知道如何将耗煤量转为报价，因此我们简单地使用线性回归，拟合估计价格到真实市场出清价格。

例如A集合[1,2,3]可通过乘以2转为B集合[2,4,6]，假如我们无法知道计算方法，可以用一个函数拟合这样的映射关系。在这里A集合就是耗煤量，B集合为实际报价

model = LinearRegression()
# 55392为训练集的长度
train_length = 55392
prices = np.array(prices).reshape(-1, 1)
X = prices[:train_length]
y = electricity_price["clearing price (CNY/MWh)"].iloc[:train_length].values.reshape(-1, 1)
model.fit(X, y)

可以看出我们的拟合方程式是：P=11.26*N-2763.16，其中N为耗煤量，P为机组报价

model.coef_, model.intercept_

(array([[11.26542666]]), array([-2763.16464251]))

使用边际成本定价的机组报价为出清价格估计

y_pred = model.predict(prices[train_length:])
y_pred = y_pred.flatten()  # 2维矩阵转为1维
y_pred[:5]

array([407.60234502, 401.74432316, 401.74432316, 391.1548221 ,
       401.74432316])

保存结果为submit.csv

成功后就可以提交到官网啦：http://competition.sais.com.cn/competitionDetail/532232/mySubmissions

sample_submit["clearing price (CNY/MWh)"] = y_pred
sample_submit.head()
daytimeclearing price (CNY/MWh)
553922023/7/10:15407.602345
553932023/7/10:30401.744323
553942023/7/10:45401.744323
553952023/7/11:00391.154822
553962023/7/11:15401.744323
sample_submit.to_csv("submit1.csv", index=False)

Task2：时间序列挖掘+ABM构建学习

学习链接：https://datawhaler.feishu.cn/wiki/VPC6wDI8Si31ZCkKknVcqWftnEf

import pandas as pd
import seaborn as sns
import matplotlib.pylab as plt
import warnings
warnings.filterwarnings('ignore')
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = ['WenQuanYi Micro Hei']
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

数据基本信息

# 这里读取的数据文件其实就是task1中经过时间戳处理的数据
electricity_price = pd.read_csv('/kaggle/input/20240730a/electricity_price_parsed.csv', parse_dates=["timestamp"], index_col=0)
electricity_price.columns = ["demand", "price"]
electricity_price.head()
demandprice
timestamp
2021-12-01 00:15:0040334.18350.80
2021-12-01 00:30:0040523.15350.80
2021-12-01 00:45:0040374.74350.80
2021-12-01 01:00:0040111.55350.80
2021-12-01 01:15:0040067.50348.93
# 创建测试集掩码，标记出所有价格为 NaN 的数据行
test_mask = electricity_price["price"].isna()

# 创建训练集掩码，标记出所有价格不为 NaN 的数据行,其中~代表布尔取反，即True和False互换。 `[True,False]`和 `~[False,True]`一致
train_mask = ~test_mask

# 打印训练集的范围和总长度
# 训练集范围从训练集中最早的日期到最新的日期
print(f"训练集范围：{electricity_price[train_mask].index.min()} --> {electricity_price[train_mask].index.max()}\t总长度{len(electricity_price[train_mask])}")

# 打印测试集的范围和总长度
# 测试集范围从测试集中最早的日期到最新的日期
print(f"测试集范围：{electricity_price[test_mask].index.min()} --> {electricity_price[test_mask].index.max()}\t总长度{len(electricity_price[test_mask])}")


训练集范围：2021-12-01 00:15:00 --> 2023-07-01 00:00:00	总长度55392
测试集范围：2023-07-01 00:15:00 --> 2024-04-19 00:00:00	总长度28128

electricity_price.info() 

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 83520 entries, 2021-12-01 00:15:00 to 2024-04-19 00:00:00
Data columns (total 2 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   demand  83520 non-null  float64
 1   price   55392 non-null  float64
dtypes: float64(2)
memory usage: 1.9 MB

进一步观察数据，发现数据中的出清价格在一小时内都是相同的，但需求会在一小时中变动。

可以写一个简单函数确认一下

electricity_price.head(12)
demandprice
timestamp
2021-12-01 00:15:0040334.18350.80
2021-12-01 00:30:0040523.15350.80
2021-12-01 00:45:0040374.74350.80
2021-12-01 01:00:0040111.55350.80
2021-12-01 01:15:0040067.50348.93
2021-12-01 01:30:0039562.02348.93
2021-12-01 01:45:0039629.90348.93
2021-12-01 02:00:0039121.31348.93
2021-12-01 02:15:0039124.08347.01
2021-12-01 02:30:0038730.58347.01
2021-12-01 02:45:0038711.77347.01
2021-12-01 03:00:0038139.21347.01
def check_repeated(data, repeat_count=4):    
    """
    检查给定数据序列中是否存在元素不断重复的情况。

    参数:
    data (list): 要检查的序列数据。
    repeat_count (int): 每个元素应重复的次数。默认值为4。

    返回:
    None
    """
    # 以步长 repeat_count 遍历 data 的索引
    for i in range(0, len(data), repeat_count):
        # 从索引 i 开始，取长度为 repeat_count 的子序列
        subsequence = data[i:i + repeat_count]
        
        # 如果子序列的独特元素数量不等于 1，则表示不是同一元素重复
        if len(set(subsequence)) != 1:
            print(f"序列数据不是元素不断重复 {repeat_count} 次")
            return  # 发现不满足条件的情况后，直接返回
        
    # 如果遍历完所有子序列，未发现不满足条件的情况，则输出满足条件的信息
    print(f"序列数据是元素不断重复 {repeat_count} 次")

# 调用函数并传入特定的电价数据
check_repeated(electricity_price[train_mask]["price"])

序列数据是元素不断重复 4 次

从图中也能确认，在一小时中出清价格一致。

# 使用loc方法选择指定日期的数据，绘制价格图表
electricity_price.loc["2023-01-03"].plot(y="price", figsize=(18, 5), marker='o')

# 设置图表的标题
plt.title("2023/1/3 Price")

Text(0.5, 1.0, '2023/1/3 Price')

在这里插入图片描述

统计指标分析

首先我们从最基本的单变量数值关系入手，探究其统计指标的表现。具体来说可以关注：

  • 中心趋势
    • 均值：df[feature].mean()
    • 中位数 df[feature].median()
    • 最大值 df[feature].max()
    • 最小值 df[feature].min()
    • 众数 df[feature].mode()
  • 变异程度
    • 标准差 df[feature].std()
    • 极差 df[feature].apply(lambda x: x.max() - x.min())
    • 四分位数 df[feature].quantile([0.25, 0.5, 0.75])
    • 变异系数 df[feature].std()/df[feature].mean()
    • 偏度和峰度 df[feature].skew(), df[feature].kurtosis()
  • 变化率：df[feature].pct_change() 计算公式为 x 2 − x 1 x 1 \dfrac{x_2-x_1}{x_1} x1​x2​−x1​​
electricity_price.describe()
demandprice
count83520.00000055392.000000
mean36595.327516360.002312
std9588.519396203.745125
min6993.310000-85.000000
25%30214.722500264.142500
50%36704.025000401.695000
75%42868.577500471.962500
max75501.3300001296.000000

数据分布

从总负荷和出清价格的数值分布来看，可以发现以下几个特点：

  • 总负荷基本服从正态分布，符合现实中的耗电规律
  • 出清价格中有大量的异常负价格
  • 1-100之间的低价格较多，表明这些时候火电
  • 有一定的异常高价（800以上）

这是一个很反常的现象，启发我们去探究负价格形成的原因。如果能找出引起低价/高价的关键信息，对模型提升会有较大帮助

# sns.displot(...) 函数用于绘制数据的分布图
ax = sns.displot(
    electricity_price,  # 输入数据，包含绘制所需的列
    x="demand",         # 指定绘图的列，这里是 "demand" 列
    aspect=1.5,         # 图形的宽高比，1.5 表示宽度是高度的 1.5 倍
    height=5,           # 图形的高度设置为 5 英寸
    kde=True            # 启用核密度估计（KDE），用于绘制数据的平滑概率密度曲线
)

# 设置图形的标题
ax.set(title="demand distripution")

在这里插入图片描述

ax=sns.displot(electricity_price,x="price",aspect=1.5,height=5,kde=True);
ax.set(title="clearing price distripution")

在这里插入图片描述

分时统计特征

不同小时的总负荷和电价

既然是时间序列，我们会考虑到指标值在不同时间上的不同表现，去发现时间上的特征趋势。首先分析不同年份下，一天不同小时的总负荷和电价。我们可以发现，数据中用电存在两个高峰期和一个低谷期

  • 早高峰：一般在6-9点
  • 晚高峰：一般在16-21点
  • 低谷期：10-15点

这似乎很反常，理论上10-15点都在进行大量的工业和商业活动，为什么反而电价和总负荷会更低呢，这里需要了解一个概念：“鸭子曲线”

鸭子曲线：由于火电和光伏发电互为替代品，当一天太阳出来后，太阳能逐渐开始替代火电，并在14点达到最大，进而导致火电受光伏发电竞争而降价。而在傍晚时太阳落山，光伏机组迅速减小发电，此时火电开始集中发电，价格迅速上升，形成了一天中典型的“两高峰，一低谷”的态势。

这样我们又增加了几个先验信息：

  • 考虑其他新能源（尤其是光伏）的影响对预测价格意义重大
  • 光伏受天气、季节等因素影响，说明还需要借助外部天气数据辅助预测。
  • 随着时间推移和中国碳中和的发展，光伏必定会在更大程度上替代火电，因此可以猜测2024年的火电价格会进一步下降。
# 提取小时信息，并创建一个新列 "hour"
electricity_price["hour"] = electricity_price.index.hour
# 提取月份信息，并创建一个新列 "month"
electricity_price["month"] = electricity_price.index.month
# 提取日期信息，并创建一个新列 "day"
electricity_price["day"] = electricity_price.index.day
# 提取星期几的信息（0 = 周一, 6 = 周日），并创建一个新列 "weekday"
electricity_price["weekday"] = electricity_price.index.weekday
# 提取年份信息，并创建一个新列 "year"
electricity_price["year"] = electricity_price.index.year

# 创建一个包含两个子图的绘图区域，图形大小为 20x9 英寸，并启用约束布局
fig, ax = plt.subplots(
    1,                       # 子图的行数为 1
    2,                       # 子图的列数为 2
    figsize=(20, 9),         # 设置整个图形的大小为 20x9 英寸
    constrained_layout=True  # 启用约束布局以自动调整子图位置和大小
)

# 绘制不同年份的分时电价线图
sns.lineplot(
    electricity_price.groupby(["hour", "year"])["price"].mean().reset_index(),  # 按小时和年份分组，计算每组的平均电价，并重置索引
    x="hour",                  # x 轴数据为小时
    y="price",                 # y 轴数据为平均电价
    ax=ax[0],                  # 将图绘制到第一个子图 (ax[0])
    marker="o",                # 数据点标记为圆圈
    hue="year",                # 根据年份不同设置不同的颜色
    palette="Set1"             # 使用 "Set1" 调色板
)
ax[0].set_title("Prices for each hour of the day in different years")  # 设置第一个子图的标题
ax[0].set_xticks(range(24))         # 设置 x 轴刻度为 0 到 23（小时范围）

# 绘制不同年份的分时总负荷线图
sns.lineplot(
    electricity_price.groupby(["hour", "year"])["demand"].mean().reset_index(),  # 按小时和年份分组，计算每组的平均总负荷，并重置索引
    x="hour",                  # x 轴数据为小时
    y="demand",                # y 轴数据为平均总负荷
    ax=ax[1],                  # 将图绘制到第二个子图 (ax[1])
    marker="o",                # 数据点标记为圆圈
    hue="year",                # 根据年份不同设置不同的颜色
    palette="Set1"             # 使用 "Set1" 调色板
)
ax[1].set_title("Demand for each hour of the day in different years")  # 设置第二个子图的标题
ax[1].set_xticks(range(24));         # 设置 x 轴刻度为 0 到 23（小时范围）

在这里插入图片描述

还可以进一步用透视表和热力图来分析不同月份下的电价变化趋势，从中我们可以发现几个趋势

  • 从1月到到5月，出清电价峰谷不断扩大，价格持续走低
  • 6-8月出清电价开始回升，并在8月达到顶峰
  • 9-12月出清电价开始迅速下降，但在11-12月略有回暖

可以根据地理知识做出猜测：

  • 从冬至到夏至，日出时间不断提早，光伏发电更早抢占火电市场空间，导致火电更早在早高峰期间跌价。同时1-5月是大风期，风电也进一步驱使火电降价。
  • 6-8月为小风期，风电减弱。同时雨带迁移到北方，光伏受天气影响发电量大幅减少，火电迎来高价区（尤其是8月降水最多，火电价格最高）
  • 9-12月为大风期，风电加强。日出开始推迟，火电在低谷期的价格逐渐回升。

这启发我们除了构造月份特征，还能构造不同时期的指示变量特征（如1-5月、6-8月）

# 创建一个透视表，计算不同月份和时间下的电价
pivot = pd.pivot_table(
    electricity_price,          # 输入的 DataFrame
    values="price",             # 透视表中要填充的值，这里是 "price"
    index="month",              # 设置行索引为月份
    columns="hour"              # 设置列索引为小时
)

# 将透视表中的数据类型转换为整数
pivot = pivot.astype(int)
pivot
hour0123456789...14151617181920212223
month
1275238214196201233295395490405...240352511656599529515458390340
2315277251228231264311391438347...177269403530586545522468408378
3385345314295303359437430356222...116215353466542486472423395416
4442405372355368429426363245127...119198325433503497493463435466
5468439415403413469443369251148...161251352422468475479471460469
6470453447441447452417338256205...267338433492509509501487490491
7385355324293284304304276245269...342389436476513542562557541465
8440390344306289315340318273293...443516594633648669648620627542
9513503475449456504498375234118...164317484590628572504462463493
10474468457442442459484412255125...198350515608584513479458454468
11363368376377371351332313285232...268355450517442375351317322355
12356328305290292319351418453373...249374558682602514476429382385

12 rows × 24 columns

# 创建一个图形，大小为 17x10 英寸
plt.figure(figsize=(17, 10))

# 绘制热图，显示不同月份和时间下的电价
sns.heatmap(
    pivot,                      # 透视表数据
    cmap="coolwarm",            # 使用 "coolwarm" 调色板，显示热图的颜色
    linewidths=0.5,             # 设置单元格之间的分隔线宽度为 0.5
    annot=True,                 # 启用单元格值的注释
    fmt=".0f",                  # 注释的格式为整数
    annot_kws={"size": 12,      # 注释文本的大小设置为 12
               "weight": "bold",  # 注释文本的字体加粗
               "color": "black"}  # 注释文本的颜色为黑色
)

# 设置图形的标题
plt.title("Prices under different months and hours")

Text(0.5, 1.0, 'Prices under different months and hours')

在这里插入图片描述

负电价与高电价形成原因分析

先前提到，数据中有大量负电价和一部分的异常高价，在这里我们探究其形成原因，并思考如何预测这种模式。

负电价

从不同小时下负电价的出现次数来看，低谷期的负电价较为明显，可能是受市场竞争导致电价中标失败，只能亏本售出。

minus_mask = electricity_price["price"] < 0

plt.figure(figsize=(12, 7))
ax = sns.countplot(electricity_price[minus_mask], x="hour")
ax.set(title="Hourly distribution of frequency occurrence of negative electricity prices")

[Text(0.5, 1.0, 'Hourly distribution of frequency occurrence of negative electricity prices')]

在这里插入图片描述

plt.figure(figsize=(12, 7))
ax = sns.countplot(electricity_price[minus_mask], x="month")
ax.set(title="Monthly distribution of frequency occurrence of negative electricity prices")

[Text(0.5, 1.0, 'Monthly distribution of frequency occurrence of negative electricity prices')]

在这里插入图片描述

通过按负电价出现频数排序，可以发现最典型的一次集中的负电价发生在2022年的五一假期中，连续出现了7次的负电价。猜测是因为火电厂发电需要人为监控（例如开关机，设置出力，而风电和光伏放着自己发电就好），假期期间大多数员工放假，从而导致火力发电量下降，出现负电价。

# 统计一下负电价出现最多的日期
result = (
    # 选择满足 minus_mask 条件的数据
    electricity_price[minus_mask]
    # 按月份和日期分组，计算每个组合的记录数量
    .groupby(["month", "day"])["price"]
    .size()                      # 计算每个分组的大小（即每个分组的记录数）
    .reset_index()              # 重置索引，使 groupby 结果成为 DataFrame，并保留分组字段为列
    .sort_values("price",       # 按 "price" 列排序
                 ascending=False)        # 降序排序
    .head(15)                   # 选择排序后的前 15 行
)
result # 下面price这一列其实是为负电价时时刻的频数
monthdayprice
945288
202176
935168
985668
724568
764964
7741060
1913157
4122656
3421956
171112956
10551652
6332752
812049
2921049
# 创建一个图形，大小为 20x8 英寸
plt.figure(figsize=(20, 8))

# 绘制从 2022年5月1日到2022年5月9日的电价数据的折线图
ax = sns.lineplot(
    electricity_price.loc["2022-05-01":"2022-05-09"]["price"],  # 从 DataFrame 中选择时间范围内的电价数据
    color="black"  # 设置折线的颜色为黑色
)

# 在每一天的 10:00 到 15:00 之间添加黄色半透明的高亮区域
for i in range(1, 10):
    ax.axvspan(
        f"2022-05-0{i} 10:00:00",  # 高亮区域的开始时间
        f"2022-05-0{i} 15:00:00",  # 高亮区域的结束时间
        color='yellow',            # 高亮区域的颜色为黄色
        alpha=0.2                  # 高亮区域的透明度设置为 0.2（0 完全透明，1 完全不透明）
    )

# 添加一条红色的虚线，y 值为 0，用于显示负电价的参考线
plt.axhline(
    y=0,                      # y 轴的值为 0
    color="red",              # 虚线的颜色为红色
    linestyle="--"           # 虚线的线型为虚线
)

# 设置 x 轴刻度标签的旋转角度为 45 度，并且水平对齐方式为右对齐
plt.setp(
    ax.get_xticklabels(),     # 获取 x 轴刻度标签
    rotation=45,              # 设置标签的旋转角度为 45 度
    ha='right'                # 设置标签的水平对齐方式为右对齐
)

# 设置图形的标题
plt.title("A period of typical negative electricity price trend: May 1, 2022 - May 9, 2022")

Text(0.5, 1.0, 'A period of typical negative electricity price trend: May 1, 2022 - May 9, 2022')

在这里插入图片描述

我们还可以进一步探究假期对电价的影响，例如2022年的大年三十是1月31日，由于调休，1月28日-1月30日上班，1月31日-2月6日放假。果然在春节假期中也出现了大量的负电价。
这说明节假日是判断负电价的重要特征

小测试：尝试画出2023年春节放假期间，1月21日-1月27日中的电价波动，并对比2023年的1月28日-2月6日电价情况，这是否更加证明假期特征的重要性？

# 假期对负电价的影响 2022年1月27日 - 2022年2月7日
plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2022-01-27":"2022-02-07"]["price"], color="black")
plt.axhline(y=0, color="red", linestyle="--")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.title("Impact of holidays on negative electricity prices January 27, 2022 - February 7, 2022")

Text(0.5, 1.0, 'Impact of holidays on negative electricity prices January 27, 2022 - February 7, 2022')

在这里插入图片描述

# 假期对负电价的影响 2023年1月21日 - 2022年2月6日
plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2023-01-21":"2023-02-06"]["price"], color="black")
plt.axhline(y=0, color="red", linestyle="--")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.title("Impact of holidays on negative electricity prices January 21, 2023 - February 6, 2023")

Text(0.5, 1.0, 'Impact of holidays on negative electricity prices January 21, 2023 - February 6, 2023')

在这里插入图片描述

高电价

和负电价相反，高电价主要集中在日落后，此时光伏发电下降，火电有较大的竞价空间。

# 计算电价的上限阈值，使用 3 标准差原则来检测异常值
upper_threshold = (
    electricity_price["price"].mean() +      # 电价的均值
    3 * electricity_price["price"].std()      # 加上 3 倍的电价标准差
)

# 创建一个布尔型掩码，用于标识电价高于上限阈值的异常值
high_abnormal_mask = (
    electricity_price["price"] > upper_threshold  # 判断电价是否大于计算得到的上限阈值
)

# 高电价频数出现的小时分布
plt.figure(figsize=(12, 7))
ax = sns.countplot(electricity_price[high_abnormal_mask], x="hour")
ax.set(title="The hourly distribution of frequency occurrences of high electricity prices")

[Text(0.5, 1.0, 'The hourly distribution of frequency occurrences of high electricity prices')]

在这里插入图片描述

# 高电价频数出现的月分布
plt.figure(figsize=(12,7))
ax=sns.countplot(electricity_price[high_abnormal_mask],x="month")
ax.set(title="Monthly distribution of frequency occurrence of high electricity prices");

在这里插入图片描述

通过对比明显高价的2022年8月3日-8月6日之间的数据，可以发现其在总需求上并不突出，甚至比后续日期中总需求更高的时间段价格还高，表明其受到了外部因素的影响

# 从电价数据中筛选出高于上限阈值的异常值
(
electricity_price[high_abnormal_mask]
    # 按月份和日期分组，并计算每个组合的异常值记录数量
    .groupby(["month", "day"])["price"]
    .size()                      # 计算每个分组的大小（即每个分组的记录数）
    .reset_index()              # 重置索引，将分组字段转换为 DataFrame 的列
    .sort_values("price",       # 按 "price" 列排序
        ascending=False)        # 降序排序
    .head(15)                   # 选择排序后的前 15 行
)
monthdayprice
2581536
228536
208320
218420
1742916
29122916
27122712
185611
51108
102108
154177
0124
3164
4184
2154
plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2022-08-01":"2022-08-08"]["demand"], color="black")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.axvspan("2022-08-03", "2022-08-06", color="yellow", alpha=0.2)
plt.title("2022/8/1 - 2022/8/8/ Demond")

Text(0.5, 1.0, '2022/8/1 - 2022/8/8/ Demond')

在这里插入图片描述

plt.figure(figsize=(20, 8))
ax = sns.lineplot(electricity_price.loc["2022-08-01":"2022-08-08"]["price"], color="black")
plt.axhline(y=upper_threshold, color="red", linestyle="--")
plt.setp(ax.get_xticklabels(), rotation=45, ha='right')
plt.title("A period of high electricity price from August 1, 2022 to August 8, 2022 from 3-6")

Text(0.5, 1.0, 'A period of high electricity price from August 1, 2022 to August 8, 2022 from 3-6')

在这里插入图片描述

双变量分析

本题中特征均为连续型变量，因此我们可以采用画散点图，计算相关系数的方法来观察两个变量的关系。

相关系数

在本案例中，demand和price的相关系数为0.34，可以看出二者有较高的线性相关性。

并且我们还能发现，在高需求下，价格变动趋缓

electricity_price[["demand","price"]].corr()
demandprice
demand1.0000000.343577
price0.3435771.000000
# 创建一个图形，大小为 8x10 英寸
plt.figure(figsize=(8, 10))

# 绘制回归图（散点图及拟合线）
sns.regplot(
    data=electricity_price.loc["2022"],  # 选择2022年的数据
    x="demand",                         # x 轴的变量为 "demand"
    y="price",                          # y 轴的变量为 "price"
    scatter_kws={                       # 设置散点图的样式
        "s": 0.5,                      # 散点的大小设置为 0.5
        "alpha": 0.6,                  # 散点的透明度设置为 0.6
        "color": "black"               # 散点的颜色设置为黑色
    },
    color="red",                        # 拟合线的颜色设置为红色
    lowess=True                          # 启用局部加权回归（Lowess）以拟合数据
)

<Axes: xlabel='demand', ylabel='price'>

在这里插入图片描述

总结

综合以上，我们可以总结出几点：

  1. 气象状况对出清价格有较大影响
  2. 节假日对出清价格有较大影响，易于出现负值
  3. 总负荷与出清价格线性关系很高，但总体呈现分段线性的特征
  4. 不同月份/小时下的出清价格受市场竞争影响较大
  5. 碳中和不断发展，火电价格有总体下降的趋势

本案例旨在提供探索性数据分析的一般步骤。可以自行尝试使用seaborn，plotly等python库做进一步的可视化，挖掘序列的更多信息，从而指导后续的特征构造。

ABM构建学习与实战

首先安装本次要用到的mesa库，其是一个方便的构建ABM的开源框架。
pip install mesa

无法避免的贫富差距——玻尔兹曼财富模型

首先我们思考一个场景，假设有这样一个世界，所有人生而平等，即持有相同的1单位资本。现在他们按照一个规则行动：如果自身持有资本大于0，将1单位资本随机给一个其他的人。一个社会按这种规则不断演化，会产生什么结果呢，试着去模拟一下。

构建环境

首先来认识一下Mesa中的环境，其定义了Agent交互的场所、时间规则等。我们可以认为每当环境做出一次迭代时，所有Agent都被唤醒，各自按照其规则执行一步。

import numpy as np
import pandas as pd
import seaborn as sns
from tqdm import tqdm
import matplotlib.pylab as plt
plt.style.use('ggplot')
# 设置matplotlib显示中文
plt.rcParams['font.sans-serif'] = ['SimHei'] 
plt.rcParams['axes.unicode_minus'] = False

# Mesa相关库
from mesa.model import Model
from mesa.time import RandomActivation
from mesa.agent import Agent
from mesa.datacollection import DataCollector

Mesa中的环境必须继承自mesa.Model，且必须实现step()方法

  • 初始化时需要定义一些环境设置
    • super().__init__()，调用父类的初始化方法，即继承一些有用的属性和方法（用到了Python面向对象编程的知识）
    • 智能体所处地图：离散方格、连续方格、网络、甚至是下载的RPG地图…… 这里我们没有设置地图。
    • 智能体调度器：描述以何种顺序触发智能体的行为。一般智能体的行动分为异步和同步，本案例定义了一个RandomActivation，会先打乱所有智能体的顺序，再按顺序依次让智能体执行自己的动作。
    • （辅助工具）DataCollector：主要用于记录实验中的一些数据，例如每个Agent的某个属性；也可以传入一个函数，接受Agent/Model进行某种处理来获取数据。
  • step：让环境迭代一次，即模拟所有Agent交互一次。可以在这里加入一些数据记录函数（比如存储每个Agent的财富，发生某种动作的Agent占比等）
def compute_gini(model): # 衡量财富或收入分配不平等程度的指标。Gini 系数的值范围从 0（完全平等）到 1（完全不平等）。
    agent_wealths = [agent.wealth for agent in model.schedule.agents] # 从模型的所有代理中提取财富信息，形成一个财富列表
    x = sorted(agent_wealths) # 将财富列表按升序排序。
    N = len(model.agents) # 代理的总数
    B = sum(xi * (N - i) for i, xi in enumerate(x)) / (N * sum(x)) # 计算加权总财富，用于后续的 Gini 系数计算。
    return 1 + (1 / N) - 2 * B # gini系数



class WealthModel(Model):
    """
    模拟财富分布的环境类，继承自 Model 类。
    """

    def __init__(self):
        """
        初始化 WealthEnv 环境。
        调用父类的初始化方法，并创建一个随机激活调度器（RandomActivation）。
        调度器每个时间步随机选择代理并激活它们，从而打乱了激活的顺序。例如有四个Agent，可能激活顺序是1,3,2,4
        """
        super().__init__()
        self.schedule = RandomActivation(self)  # 创建一个随机激活调度器，用于管理代理的激活顺序
        self.datacollector = DataCollector(model_reporters={"gini":compute_gini},agent_reporters={"Wealth": "wealth"})

    def step(self) -> None:
        """
        进行一步模拟。调用调度器的 step 方法，激活所有代理进行一步操作。
        """
        self.datacollector.collect(self) # 在每一步模拟中，收集模型和代理的数据。
        self.schedule.step()  # 调用调度器的 step 方法，激活所有代理进行一步操作

    def add_agents(self, agents: list[Agent]):
        """
        添加代理到环境中。

        参数:
        agents (list[Agent]): 要添加的代理列表。
        """
        for agent in agents:
            self.schedule.add(agent)  # 将每个代理添加到调度器中
定义agent

回想之前提到的规则：如果自身持有资本不为0，将1单位资本随机给一个其他的人。可以大致构想，假设我是这个Agent，可能会按下面的流程与其他Agent互动：

  1. 检查自己的财富状况，如果为0，则什么也不做
  2. 如果财富大于0，从除自己外的所有Agent中随机抽一个Agent
  3. 将对方财富加1，自身财富减1（也可以尝试其他数字）
  4. 停止行动，等待下一轮环境迭代

mesa的Agent必须有unique_id，model两个属性：

  • unique_id：用于区分Agent
  • model：先前定义好的环境

此外Agent必须实现step方法，表示让Agent根据自己规则行动一步

class WealthAgent(Agent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 1

    def step(self):
        if self.wealth > 0: 
            # 如果财富大于0，从除自己外的所有Agent中随机抽一个Agent
            # self.random.choice等价于random库的choice
            # self.model.schedule.agents能获得环境中的agent列表
            other_agent = self.random.choice(self.model.schedule.agents)
            other_agent.wealth += 1
            self.wealth -= 1

目前我们就已经完成了ABM的构建了。下面我们用1000个Agent运行200个时间步，看看最终的总体财富分布如何：

steps = 200

env = WealthModel()
agents = [WealthAgent(i, env) for i in range(1000)] # 定义Agent列表，里面每个元素为交互的主体
env.add_agents(agents) # # 将Agent放入环境中

# 初始每个Agent拥有的财富是一样的
# [agent.wealth for agent in agents] # 都是1

# env.step() # 环境迭代一步
# 迭代一步后Agent各自交互，财富发生变化
# [agent.wealth for agent in agents] # 有的变为0，有的增加

for i in tqdm(range(steps)):
    env.step()

100%|██████████| 200/200 [01:47<00:00,  1.86it/s]

可以发现最终财富趋向于两级分化：80%的人持有20%的财富，而20%的人持有80%的财富。

同时我们可以监控基尼系数（衡量财富不平等的统计指标，0为绝对平等，1为绝对不平等）。发现基尼系数迅速达到并稳定在了0.6左右，表示有相当高的贫富差距。

可以通过数学证明，最终的财富分布为物理学中的玻尔兹曼分布，试想这是为什么？

使用env.datacollector.get_agent_vars_dataframe()获取监控数据Dataframe。第一列固定为模拟的时间步，第二列固定为每个Agent的ID，后续列为监控的数值

df = env.datacollector.get_agent_vars_dataframe()
df
Wealth
StepAgentID
001
11
21
31
41
.........
1996230
3240
1254
00
1260

200000 rows × 1 columns

last_wealth=df.xs(steps-1,level="Step") # 多级索引的提取方法，xs指定获取哪个层级下，为指定值的行
last_wealth
Wealth
AgentID
1362
8702
9802
1712
4770
......
6230
3240
1254
00
1260

1000 rows × 1 columns

ax=sns.histplot(last_wealth, discrete=True)
ax.set(title="The final wealth distribution",xlabel="Wealth")


[Text(0.5, 1.0, 'The final wealth distribution'), Text(0.5, 0, 'Wealth')]

在这里插入图片描述

model_df=env.datacollector.get_model_vars_dataframe()
model_df
gini
00.000000
10.517086
20.571684
30.600080
40.600286
......
1950.626884
1960.640286
1970.620418
1980.634310
1990.637590

200 rows × 1 columns

plt.figure(figsize=(12,7))
ax=sns.lineplot(model_df)
ax.set(title='Change in Gini coefficient',ylabel="Gini coefficient",xlabel="Step");

在这里插入图片描述

当然，我们还能尝试调整模型/智能体的参数，并观察还会有什么有趣的结果：

  • 富二代模型：现在不是所有智能体生而平等了，有些智能体天生就有更多的财富，倾向于奢侈生活（每次给出2单位资本），但也因为地位而更容易获得资本（一次增加2单位资本）。

现在我们换成普通智能体一开始持有50元，富二代则是250元，普通智能体初始个数为280，富二代为20，运行5000步：
可以发现现在由于普通智能体都有一定资本了，富二代可以靠吸取底层财富保持自己相对稳定的财富水平，而普通智能体资产很难与他们的财富有交集。当然一些富二代也会因为过度挥霍而沦为普通智能体。

class WealthAgent(Agent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 50

    def step(self):
        if self.wealth > 0:
            other_agent = self.random.choice(self.model.schedule.agents)
            if type(other_agent) == RichSecondAgent:
                if self.wealth >= 2:
                    self.wealth -= 2
                    other_agent.wealth += 2
                    return
            other_agent.wealth += 1
            self.wealth -= 1


class RichSecondAgent(Agent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)
        self.wealth = 250

    def step(self):
        if self.wealth > 0:
            other_agent = self.random.choice(self.model.schedule.agents)
            other_agent.wealth += 2
            self.wealth -= 2

from tqdm import tqdm
steps = 5000

env = WealthModel()
agents = [WealthAgent(i, env) for i in range(280)]
agents+=[RichSecondAgent(i,env)  for i in range(280,300)]
env.add_agents(agents)

for i in tqdm(range(steps)):
    env.step()

100%|██████████| 5000/5000 [05:37<00:00, 14.84it/s]

agent_df=env.datacollector.get_agent_vars_dataframe()
wealths=agent_df[agent_df.index.get_level_values("AgentID").isin([2,3,4,280,285,297])]
plt.figure(figsize=(20,9))
sns.lineplot(wealths,hue="AgentID",x="Step",y="Wealth",palette="Set1");

在这里插入图片描述

当然还能尝试其他模型，看看会涌现出什么特性

  • 负债累累模型：现在允许智能体负债，但一旦负债超过100单位，智能体将出局。
  • 劫富济贫模型：每10个时间步，最高与最低的智能体财富将会被三七分。
  • 求锦鲤模型：最高财富的智能体将会在每100步中进行一次抽奖，将自身财富10%随机给一个其他智能体，获奖的智能体此后将自身花费增加2单位

是不是感觉很神奇，仅仅一个规则就模拟出了财富分布不均、社会阶层固化等现象。事实上如果参数和交互更加复杂，ABM还能更为真实地模拟和解释生活中的其他经济学现象。

如果对Mesa的其他模型感兴趣，还可以参考Mesa Example做进一步的复现。

Task3：进阶上分方向讨论

学习链接：https://datawhaler.feishu.cn/wiki/DafPwzGSriBXd7kaMdtcngQlnzc

本节设置了三个主要的上分方向,可以自由尝试：

  • 时间序列挖掘（较容易）
  • ABM报价策略优化（难度中等偏上）
  • 强化学习（难度高）

时间序列挖掘

本baseline主要展示如何用时间序列方法解决该比赛，主要特点为：

  1. 使用了一些基本的时间序列特征构造方法
  2. 以线性回归模型和树模型LightGBM为例训练模型
  3. ensemble方法：模型预测结果融合
  4. 基于先验信息调整预测结果
import numpy as np
import pandas as pd
import seaborn as sns
import os
from tqdm import tqdm
import pandas as pd
import seaborn as sns
import matplotlib.pylab as plt
import warnings
warnings.filterwarnings('ignore')
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = ["WenQuanYi Micro Hei",'SimHei']
plt.rcParams['axes.unicode_minus'] = False

train_length=55392

electricity_price=pd.read_csv('/kaggle/input/20240730a/electricity_price_parsed.csv',parse_dates=["timestamp"],index_col=0)
electricity_price.columns=["demand","price"]
sample_submit=pd.read_csv('/kaggle/working/sample_submit.csv')

train_data=electricity_price.copy(deep=True)
electricity_price.head()
demandprice
timestamp
2021-12-01 00:15:0040334.18350.80
2021-12-01 00:30:0040523.15350.80
2021-12-01 00:45:0040374.74350.80
2021-12-01 01:00:0040111.55350.80
2021-12-01 01:15:0040067.50348.93

特征工程

# 将电价数据的时间索引信息添加到训练数据中

# 提取时间索引的小时信息，并添加到训练数据中，创建 "hour" 列
train_data["hour"] = electricity_price.index.hour

# 提取时间索引的日期信息，并添加到训练数据中，创建 "day" 列
train_data["day"] = electricity_price.index.day

# 提取时间索引的月份信息，并添加到训练数据中，创建 "month" 列
train_data["month"] = electricity_price.index.month

# 提取时间索引的年份信息，并添加到训练数据中，创建 "year" 列
train_data["year"] = electricity_price.index.year

# 提取时间索引的星期信息，并添加到训练数据中，创建 "weekday" 列
train_data["weekday"] = electricity_price.index.weekday

# 根据月份信息，判断是否为风季（1-5月和9-12月），创建布尔型 "is_windy_season" 列
train_data["is_windy_season"] = electricity_price.index.month.isin([1, 2, 3, 4, 5, 9, 10, 11, 12])

# 根据小时信息，判断是否为低谷时段（10-15点），创建布尔型 "is_valley" 列
train_data["is_valley"] = electricity_price.index.hour.isin([10, 11, 12, 13, 14, 15])

# 提取时间索引的季度信息，并添加到训练数据中，创建 "quarter" 列
train_data["quarter"] = electricity_price.index.quarter

# 对时间特征进行独热编码（One-Hot Encoding），删除第一列以避免多重共线性
train_data = pd.get_dummies(
    data=train_data,        # 需要进行独热编码的 DataFrame
    columns=["hour", "day", "month", "year", "weekday"],  # 需要独热编码的列
    drop_first=True         # 删除第一列以避免多重共线性
)

基于Task2中的EDA，我们发现节假日中电力价格更容易为负值，因此我们构造典型的春节和劳动节的节假日特征。当然你也可以尝试构造其他的节假日特征（端午节、清明节、国庆节……）

可以使用python的holidays库来构造节假日特征，但考虑到其给出的日期范围不够准确，所以这里进行了手动构造

def generate_holiday_dates(start_dates, duration):
    holidays = []  # 初始化一个空列表，用于存储节假日日期
    for start_date in start_dates:  # 遍历每个节假日的开始日期
        # 生成从 start_date 开始的日期范围，持续时间为 duration 天
        holidays.extend(pd.date_range(start=start_date, periods=duration).tolist())
    return holidays  # 返回所有节假日日期的列表

spring_festival_start_dates = ["2022-01-31", "2023-01-21", "2024-02-10"]
labor_start_dates = ["2022-04-30", "2023-04-29"]
spring_festivals = generate_holiday_dates(spring_festival_start_dates, 7)
labor = generate_holiday_dates(labor_start_dates, 5)

# 判断训练数据的索引是否在春节日期列表中，生成布尔型列 "is_spring_festival"
train_data["is_spring_festival"] = train_data.index.isin(spring_festivals)

# 判断训练数据的索引是否在劳动节日期列表中，生成布尔型列 "is_labor"
train_data["is_labor"] = train_data.index.isin(labor)

可以发现一段时间内，总负荷如果有下降趋势，则下一个出清价格也可能下降。因此联想到构造基于总需求的窗口特征。

def cal_range(x):
    """
    计算极差（最大值和最小值之差）。

    参数：
    x (pd.Series): 输入的时间序列数据。

    返回：
    float: 极差值。

    示例：
    >>> import pandas as pd
    >>> x = pd.Series([1, 2, 3, 4, 5])
    >>> cal_range(x)
    4
    """
    return x.max() - x.min()


def increase_num(x):
    """
    计算序列中发生增长的次数。

    参数：
    x (pd.Series): 输入的时间序列数据。

    返回：
    int: 序列中增长的次数。

    示例：
    >>> x = pd.Series([1, 2, 3, 2, 4])
    >>> increase_num(x)
    3
    """
    return (x.diff() > 0).sum()


def decrease_num(x):
    """
    计算序列中发生下降的次数。

    参数：
    x (pd.Series): 输入的时间序列数据。

    返回：
    int: 序列中下降的次数。

    示例：
    >>> x = pd.Series([1, 2, 1, 3, 2])
    >>> decrease_num(x)
    2
    """
    return (x.diff() < 0).sum()


def increase_mean(x):
    """
    计算序列中上升部分的均值。

    参数：
    x (pd.Series): 输入的时间序列数据。

    返回：
    float: 序列中上升部分的均值。

    示例：
    >>> x = pd.Series([1, 2, 3, 2, 4])
    >>> diff = x.diff()
    >>> diff
    0    NaN
    1    1.0
    2    1.0
    3   -1.0
    4    2.0
    dtype: float64
    >>> increase_mean(x)
    1.33
    """
    diff = x.diff()
    return diff[diff > 0].mean()


def decrease_mean(x):
    """
    计算序列中下降的均值（取绝对值）。

    参数：
    x (pd.Series): 输入的时间序列数据。

    返回：
    float: 序列中下降的均值（绝对值）。

    示例：
    >>> import pandas as pd
    >>> x = pd.Series([4, 3, 5, 2, 6])
    >>> decrease_mean(x)
    2.0
    """
    diff = x.diff()
    return diff[diff < 0].abs().mean()


def increase_std(x):
    """
    计算序列中上升部分的标准差。

    参数：
    x (pd.Series): 输入的时间序列数据。

    返回：
    float: 序列中上升部分的标准差。

    示例：
    >>> import pandas as pd
    >>> x = pd.Series([1, 2, 3, 2, 4])
    >>> increase_std(x)
    0.5773502691896257
    """
    diff = x.diff()
    return diff[diff > 0].std()


def decrease_std(x):
    """
    计算序列中下降部分的标准差。

    参数：
    x (pd.Series): 输入的时间序列数据。

    返回：
    float: 序列中下降部分的标准差。

    示例：
    >>> import pandas as pd
    >>> x = pd.Series([4, 3, 5, 2, 6])
    >>> decrease_std(x)
    1.4142135623730951
    """
    diff = x.diff()
    return diff[diff < 0].std()


from tqdm import tqdm  # 导入 tqdm 库用于显示进度条

# 定义滚动窗口大小的列表
window_sizes = [4, 12, 24]

# 遍历每个窗口大小
with tqdm(window_sizes) as pbar:
    for window_size in pbar:
        # 定义要应用的聚合函数列表
        functions = ["mean", "std", "min", "max", cal_range, increase_num,
                     decrease_num, increase_mean, decrease_mean, increase_std, decrease_std]

        # 遍历每个聚合函数
        for func in functions:
            # 获取函数名称，如果是字符串则直接使用，否则使用函数的 __name__ 属性
            func_name = func if type(func) == str else func.__name__

            # 生成新列名，格式为 demand_rolling_{window_size}_{func_name}
            column_name = f"demand_rolling_{window_size}_{func_name}"

            # 计算滚动窗口的聚合值，并将结果添加到 train_data 中
            train_data[column_name] = train_data["demand"].rolling(
                window=window_size,        # 滚动窗口大小
                min_periods=window_size//2,  # 最小观测值数
                closed="left"         # 滚动窗口在左侧闭合
            ).agg(func)              # 应用聚合函数

            pbar.set_postfix({"window_size": window_size, "func": func_name})


100%|██████████| 3/3 [09:09<00:00, 183.25s/it, window_size=24, func=decrease_std]

我们还可以加入一些其他的时序特征，例如demand的滞后特征，差分特征，百分比特征等。通常我们会尝试构造尽可能多的特征，随后在特征筛选中留下对结果有显著影响的特征。

其他时序特征

  1. shift 方法

    • 描述: 将序列中的值向前或向后移动指定的位数，填充移位后产生的缺失值。
    • 语法: Series.shift(periods=1, freq=None, axis=0, fill_value=None)
    • 参数:
      • periods：移位的时期数，正数表示向后移位，负数表示向前移位。
      • freq：时间间隔或时间偏移量，可选。
      • axis：要移动的轴，默认为0。
      • fill_value：移位后填充缺失值的值，可选。
    • 示例:
      import pandas as pd
demand = pd.Series([10, 20, 30, 40])
demand_shift_1 = demand.shift(1)
print(demand_shift_1)
# 输出：
# 0     NaN
# 1    10.0
# 2    20.0
# 3    30.0
# dtype: float64

  2. diff 方法

    • 描述: 计算序列中相邻值之间的差异，结果是当前值减去前一个值。
    • 语法: Series.diff(periods=1, axis=0)
    • 参数:
      • periods：要计算差异的时期数，默认为1。
      • axis：要计算差异的轴，默认为0。
    • 示例:
      import pandas as pd
demand = pd.Series([10, 20, 30, 40])
demand_diff_1 = demand.diff(1)
print(demand_diff_1)
# 输出：
# 0    NaN
# 1    10.0
# 2    10.0
# 3    10.0
# dtype: float64

  3. pct_change 方法

    • 描述: 计算序列中相邻值的百分比变化，结果是当前值减去前一个值再除以前一个值。
    • 语法: Series.pct_change(periods=1, fill_method='pad', limit=None, freq=None)
    • 参数:
      • periods：要计算百分比变化的时期数，默认为1。
      • fill_method：用于填充缺失值的方法（如 ‘pad’ 或 ‘bfill’），默认为 ‘pad’。
      • limit：填充时的最大数量，可选。
      • freq：时间间隔或时间偏移量，可选。
    • 示例:
      import pandas as pd
demand = pd.Series([10, 20, 30, 40])
demand_pct_1 = demand.pct_change(1)
print(demand_pct_1)
# 输出：
# 0         NaN
# 1    1.000000
# 2    0.500000
# 3    0.333333
# dtype: float64

# 添加新的特征列：demand_shift_1，表示将 demand 列中的值向后移动一位
# shift(1) 的结果是当前行的值等于前一行的值，第一行的值为 NaN
train_data["demand_shift_1"] = train_data["demand"].shift(1)

# 添加新的特征列：demand_diff_1，表示 demand 列中相邻值的差
# diff(1) 的结果是当前行的值减去前一行的值，第一行的值为 NaN
train_data["demand_diff_1"] = train_data["demand"].diff(1)

# 添加新的特征列：demand_pct_1，表示 demand 列中相邻值的百分比变化
# pct_change(1) 的结果是当前行的值减去前一行的值再除以前一行的值，第一行的值为 NaN
train_data["demand_pct_1"] = train_data["demand"].pct_change(1)

模型训练

一般数据科学比赛中，线性模型和树模型是比较常用的两种回归模型，在这里我们将线性回归和LightGBM作为Baseline的两个模型。
可以尝试其他模型例如：

  • 基于深度学习的方法：LSTM、NBeats、Transformer……
  • 基于时序模型的方法：ARIMA、Prophet、VARMAX……
  • 其他线性模型：Ridge、ElasticNet、Lasso……
  • 其他树模型：CatBoost、XGBoost、LightGBM

这里我们使用bfill()和ffill()主要因为获取时序特征会形成缺失值，比如我们将数据平移一层，那么第一个值由于没有前面的值，就会变为NaN，即缺失值。当然你可以使用其他缺失值填充方式（KNN插值、线性插值、模型预测缺失值）

本例中我们使用最简单的基于时序的填充方法。

# 从 train_data 中创建训练集和测试集特征数据 (X) 和目标数据 (y)

# 创建训练集特征数据 X_train
# 1. 从 train_data 中选择前 train_length 行，去除 "price" 列
# 2. 使用 bfill 方法向后填充缺失值
# 3. 使用 ffill 方法向前填充缺失值
X_train = train_data.iloc[:train_length].drop(columns=["price"]).bfill().ffill()

# 创建测试集特征数据 X_test
X_test = train_data.iloc[train_length:].drop(columns=["price"]).bfill().ffill()

# 创建训练集目标数据 y_train
y_train = train_data.iloc[:train_length][["price"]]

from sklearn.linear_model import LinearRegression
from lightgbm import LGBMRegressor

# 创建 LGBMRegressor 模型对象，设置参数
# num_leaves：树的叶子数，控制模型复杂度
# n_estimators：迭代次数，即树的数量
# verbose：日志信息级别，-1 表示不输出日志信息
lgb_model = LGBMRegressor(num_leaves=2**5-1, n_estimators=300, verbose=-1)

# 创建线性回归模型对象
linear_model = LinearRegression()

# 使用训练集数据训练 LGBMRegressor 模型
# X_train：训练集特征数据
# y_train：训练集目标数据
lgb_model.fit(X_train, y_train)

# 使用训练集数据中的 "demand" 特征训练线性回归模型
# X_train[["demand"]]：训练集特征数据中仅包含 "demand" 列
# y_train：训练集目标数据
linear_model.fit(X_train[["demand"]], y_train)

# 使用训练好的 LGBMRegressor 模型预测测试集特征数据
# X_test：测试集特征数据
# 返回预测的目标值
lgb_pred = lgb_model.predict(X_test)

# 使用训练好的线性回归模型预测测试集特征数据中的 "demand" 列
# X_test[["demand"]]：测试集特征数据中仅包含 "demand" 列
# 返回预测的目标值，并将结果展平为一维数组
linear_pred = linear_model.predict(X_test[["demand"]]).flatten()

模型融合

考虑到demand对价格的线性程度高，但直接线性回归无法捕捉到时序中的一些非线性特征，因此我们融合树模型和线性模型的结果，由于我们没有其他先验信息，选择直接取均值。

此外还记得我们在Task2中提到，由于鸭子曲线，火电厂在低谷期时的价格会越来越低，并且整体火电价格也会下降，因此我们可以对最终预测结果进行一个整体上的缩放，保证符合未来的趋势。

# 简单求均值
y_pred = (lgb_pred+linear_pred)/2
y_pred *= 0.95  # 进行少量修正

sample_submit["clearing price (CNY/MWh)"] = y_pred
sample_submit.head()
daytimeclearing price (CNY/MWh)
02023/7/10:15391.704814
12023/7/10:30386.056237
22023/7/10:45379.044964
32023/7/11:00376.880418
42023/7/11:15382.486849
sample_submit.to_csv("submit3.csv", index=False)

上述结果提交后分数为：11302.0459。

# 定义每4行计算均值并替代原值的函数
def replace_with_group_mean(df, column_name, group_size):
    df[column_name] = df[column_name].groupby(df.index // group_size).transform('mean')
    return df

# 使用函数处理 'clearing price (CNY/MWh)' 列
sample_submit = replace_with_group_mean(sample_submit, 'clearing price (CNY/MWh)', 4)
sample_submit.head()
daytimeclearing price (CNY/MWh)
02023/7/10:15383.421608
12023/7/10:30383.421608
22023/7/10:45383.421608
32023/7/11:00383.421608
42023/7/11:15368.875564
sample_submit.to_csv("submit3.1.csv", index=False)

上述结果提交后分数为：11152.1809。

下面进行调参看是否有所提高。下面使用了GPU。

from sklearn.model_selection import GridSearchCV

param_grid = {
    'num_leaves': [31, 63, 127],
    'max_depth': [-1, 10, 20],
    'learning_rate': [0.1, 0.01, 0.001],
    'n_estimators': [100, 200, 300],
    'subsample': [0.8, 1.0],
    'colsample_bytree': [0.8, 1.0]
}

grid_search = GridSearchCV(estimator=LGBMRegressor(device='gpu',verbose=-1),
                           param_grid=param_grid,
                           scoring='neg_mean_squared_error',
                           cv=5)
grid_search.fit(X_train, y_train)
lgb_model = grid_search.best_estimator_
lgb_model.fit(X_train, y_train)
lgb_pred = lgb_model.predict(X_test)

X_train['demand_squared'] = X_train['demand'] ** 2
X_test['demand_squared'] = X_test['demand'] ** 2
linear_model = LinearRegression()
linear_model.fit(X_train[["demand",'demand_squared']], y_train)
linear_pred = linear_model.predict(X_test[["demand",'demand_squared']]).flatten()

# 简单求均值
y_pred = (lgb_pred+linear_pred)/2
y_pred *= 0.95  # 进行少量修正

sample_submit["clearing price (CNY/MWh)"] = y_pred

# 使用函数处理 'clearing price (CNY/MWh)' 列
sample_submit = replace_with_group_mean(sample_submit, 'clearing price (CNY/MWh)', 4)
sample_submit.head()
daytimeclearing price (CNY/MWh)
02023/7/10:15382.140680
12023/7/10:30382.140680
22023/7/10:45382.140680
32023/7/11:00382.140680
42023/7/11:15357.452726
sample_submit.to_csv("submit3.2.csv", index=False)

上述结果提交后分数为：11210.5533。

综上，我这里最好的分数是：
在这里插入图片描述

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/1982050.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

错误代码0x80004005：全面解析与解决策略

错误代码0x80004005：全面解析与解决策略

在Windows操作系统中&#xff0c;错误代码0x80004005是一个相当常见的提示&#xff0c;它通常表示一个未指定的错误或一个无法找到的错误。这个错误可能出现在多种场景中&#xff0c;包括软件安装、系统更新、网络连接、文件操作等。由于其普遍性和模糊性&#xff0c;理解和解决…
阅读更多...
搭建链路追踪系统 Jaeger

搭建链路追踪系统 Jaeger

搭建链路追踪系统 Jaeger 架构 架构图 架构设计 在设计以 Jaeger 为中心的链路追踪系统时&#xff0c;我们采用微服务架构模式。 Jaeger 作为核心追踪系统&#xff0c;负责收集、存储和展示分布式追踪数据。Elasticsearch 作为后端存储&#xff0c;提供数据的持久化和高效…
阅读更多...
网鼎杯之unfinish

网鼎杯之unfinish

题目来源&#xff1a;网鼎杯 2018 题目描述&#xff1a;SQL 问题&#xff1a;给一个页面&#xff0c;如何判断其有漏洞 通过awvs扫描得知存在 register.php 注册页面&#xff0c;并且注册界面存在SQL盲注漏洞。 题目提示SQL&#xff0c;在注册页面测试无果&#xff0c;且发现输…
阅读更多...
解决Ubuntu报“无法解析域名cn.archive.ubuntu.com“问题

解决Ubuntu报“无法解析域名cn.archive.ubuntu.com“问题

今天在Ubuntu系统上&#xff0c;使用sudo apt update命令&#xff0c;进行更新时&#xff0c;弹出"无法解析域名 cn.archive.ubuntu.com"问题&#xff0c;如图(1)所示&#xff1a; 图(1) 弹出"无法解析域名 cn.archive.ubuntu.com" 错误 出现这种现象的原因…
阅读更多...
C++并发编程(一)：线程基础

C++并发编程(一)：线程基础

简介 本文学习的是 b 站 up 恋恋风辰的并发编程教学视频做的一些笔记补充。 教程视频链接如下&#xff1a;线程基础&#xff1a;视频教程 文档链接如下&#xff1a;线程基础&#xff1a;笔记文档 理论上直接看 up 提供的笔记文档即可&#xff0c;我这里主要是记录一些我自己…
阅读更多...
免费qq号码估价的工具和软件

免费qq号码估价的工具和软件

目前有多种 qq 号码估价的工具和软件。例如&#xff0c;晒号网的 QQ 估价器可以根据 QQ 号码等级、QQ 号码资深度、QQ 号码年限、活跃时间等进行准确的 QQ 号码估价。此外&#xff0c;还有其他一些相关的估价软件和平台&#xff0c;如 QQ 号码估价 2.0 全新玩法&#xff0c;利用…
阅读更多...
复现opendrivelab的“点云预测”项目

复现opendrivelab的“点云预测”项目

本文的主要工作就是复现下述论文中的算法。 该论文全称&#xff1a;Visual Point Cloud Forecasting 论文内容在此不做过多介绍&#xff0c;直接上项目。 一、准备工作 首先通读readme.md文件的内容&#xff0c;了解所需要的相关依赖和数据等内容。 一定要多读几遍&#xf…
阅读更多...
C# udp通信测试助手

C# udp通信测试助手

1、UI界面和最终实现功能测试 2、代码 using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Net; using System.Net.Sockets; using System.Text; using System.Threa…
阅读更多...
【C++BFS算法】2998. 使 X 和 Y 相等的最少操作次数

【C++BFS算法】2998. 使 X 和 Y 相等的最少操作次数

本文涉及知识点 CBFS算法 LeetCode2998. 使 X 和 Y 相等的最少操作次数 给你两个正整数 x 和 y 。 一次操作中&#xff0c;你可以执行以下四种操作之一&#xff1a; 如果 x 是 11 的倍数&#xff0c;将 x 除以 11 。 如果 x 是 5 的倍数&#xff0c;将 x 除以 5 。 将 x 减 …
阅读更多...
Linux项目实战——业务系统（后附百问网课程视频链接）

Linux项目实战——业务系统（后附百问网课程视频链接）

处理配置文件 一、main函数流程 初始化显示系统——>初始化输入系统——>初始化文字系统——>初始化页面系统——>业务系统 二、业务系统框架 1. 读取配置文件 2. 根据配置文件生成按钮、界面 3. 读取输入事件 4. 根据输入事件找到按钮 5. 调用按键的OnPressed函…
阅读更多...
小程序开发怎么申请HTTPS证书？

小程序开发怎么申请HTTPS证书？

小程序开发中申请SSL证书的流程包括选择信任可靠的SSL证书服务商、申请SSL证书、配置服务器、测试验证以等。以下将具体分析这个过程&#xff1a; 1、选择信任可靠的SSL证书服务商&#xff1a; 选择信任可靠的证书服务商&#xff0c;致命的证书服务商可以提供可靠的产品和完善…
阅读更多...
【已解决】页面操作系统功能，诡异报错500nginx错误

【已解决】页面操作系统功能，诡异报错500nginx错误

【已解决】页面操作系统功能&#xff0c;诡异报错500nginx错误&#xff0c;后台没有任何报错信息 不知道啥原因 清理了浏览器缓存 也没有效果 还有一个表现情况&#xff0c;同样的操作&#xff0c;有时可以又是不行 因为报错ng的代理问题&#xff0c;检查了ng配置 后续经过同…
阅读更多...
Unity强化工程 之 SpriteEditer SingleMode

Unity强化工程 之 SpriteEditer SingleMode

本文仅作笔记学习和分享&#xff0c;不用做任何商业用途 本文包括但不限于unity官方手册&#xff0c;unity唐老狮等教程知识&#xff0c;如有不足还请斧正 因为unity不只是3d需要&#xff0c;还有2d游戏需要大量编辑处理图片素材&#xff0c;所以需要了解Sprite&#xff08;精灵…
阅读更多...
Cesium初探-实体

Cesium初探-实体

在 Cesium 中&#xff0c;"实体"&#xff08;Entity&#xff09;是一个核心概念&#xff0c;它代表了可以在场景中渲染的任何东西&#xff0c;从简单的点、线、多边形到复杂的模型和图像。实体可以用来表示地理空间数据&#xff0c;如地标、轨迹、卫星轨道等。 以下…
阅读更多...
【数据结构】栈的概念、结构和实现详解

【数据结构】栈的概念、结构和实现详解

本文来介绍一下数据结构中的栈&#xff0c;以及如何用C语言去实现。 1. 栈的概念及结构 栈&#xff1a;一种特殊的线性表&#xff0c;它只允许在固定的一端进行插入和删除元素的操作。 进行数据插入和删除操作的一端称为栈顶&#xff0c;另一端称为栈底。 栈中元素遵循后进先出…
阅读更多...
PXE实现linux系统批量自动安装

PXE实现linux系统批量自动安装

实验环境&#xff1a; 1、RHEL7主机 2、主机图形化 3、配置网络可用 4、关闭vmware dhcp功能 一、实验环境准备 1、主机图形化 在安装安装RHEL7系统时&#xff0c;选择图形化安装&#xff0c;如果没有选择&#xff0c;可以在后面通过命令安装&#xff0c;如下&#xff1…
阅读更多...
Pinia状态管理库

Pinia状态管理库

为了跨组件传递JWT令牌&#xff0c;我们就会利用Pinia状态管理库&#xff0c;它允许跨组件或页面共享状态。 使用Pinia步骤&#xff1a; 安装pinia&#xff1a;cnpm install pinia 在vue应用实例中使用pinia 在src/stores/token.js中定义store 在组件中使用store 1.在main.js文…
阅读更多...
sql注入靶场sqli-labs常见sql注入漏洞详解

sql注入靶场sqli-labs常见sql注入漏洞详解

目录 sqli-labs-less1 1.less1普通解法 1.在url里面填写不同的值&#xff0c;返回的内容也不同&#xff0c;证明&#xff0c;数值是进入数据库进行比对了的&#xff08;可以被注入&#xff09; 2.判断最终执行的sql语句的后面还有内容吗&#xff0c;并且能够判断是字符型的拼接…
阅读更多...
MySQL增删改查（基础）

MySQL增删改查（基础）

1、. 新增&#xff08;Create&#xff09; 语法&#xff1a; INSERT [INTO] table_name[(column [, column] ...)] VALUES (value_list) [, (value_list)] ... 例子&#xff1a; -- 创建一张学生表 DROP TABLE IF EXISTS student; CREATE TABLE student (id INT,sn INT com…
阅读更多...
电子琴——Arduino

电子琴——Arduino

音调有7个音调&#xff0c;分别是哆来咪发索莱西&#xff1b;如果用蜂鸣器来发出这七个音调就要分别设置这七个音调对应频率。 电子琴实现需要物品有&#xff0c;arduino开发板一个&#xff0c;按键7个&#xff0c;蜂鸣器1个&#xff0c;杜邦线若干 重点讲一下按键原理 按键开…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023