【风况预测评分规则-python实现】第五届全国工业互联网数据创新应用大赛-机组数据驱动的风电场短期风况预测

比赛官网

注：不保证验证代码百分百正确，从两种实现方式多次实验取得的结果相差不大来看，错误的可能性应该也不会太大

import numpy as np
import pandas as pd

""" 部分变量名来源（采用谷歌翻译，省略“风”）
风场: field
风机: machine
时段: period
时刻: time
风速: speed
风向: direction
"""

field = ['风场1', '风场2']
field_len = len(field)

machine = [f'x{i}' for i in list(range(26, 50+1)) + list(range(25, 49+1))]
machine_len = len(machine) // field_len

period = [f'{s}_{str(i).zfill(2)}' for s in ['春', '夏', '秋', '冬'] for i in range(1, 20+1)]
period_len = len(period)

time = [i * 30 for i in range(1, 20+1)]
time_len = len(time)

def generate_data():
    return pd.Dataframe(
        data={
            '风场': np.array(field * machine_len * period_len * time_len).reshape(-1, field_len).T.reshape(-1),
            '风机': np.array(machine * period_len * time_len).reshape(-1, field_len * machine_len).T.reshape(-1),
            '时段': np.array(np.array(period * time_len).reshape(-1, period_len).T.reshape(-1).tolist() * field_len * machine_len),
            '时刻': np.array(time * field_len * machine_len * period_len),
            '风速': np.around(np.random.random(field_len * machine_len * period_len * time_len), 4),
            '风向': np.around(np.random.random(field_len * machine_len * period_len * time_len), 4)
        }
    )

# 临时预测数据 = 生成数据 + 噪音
df_true = generate_data()
df_pred = df_true.copy()
df_pred[['风速', '风向']] = df_pred[['风速', '风向']] + np.random.random() * 0.2

%%time

# 用到的时候再算，自下向上实现

# 公式1
R = lambda : 100 / (1 + np.sum([E(0, i) for i in range(80)]) + np.sum([E(1, i) for i in range(80)]))

# 公式2
def E(f, i, p=machine_len):
    df_true_field = df_true[df_true['风场'] == field[f]]
    df_true_period = df_true_field[df_true_field['时段'] == period[i]]
    
    df_pred_field = df_pred[df_pred['风场'] == field[f]]
    df_pred_period = df_pred_field[df_pred_field['时段'] == period[i]]
    
    value = 0
    for j in range(p):
        v_true, d_true = np.vsplit(df_true_period[df_true_period['风机'] == machine[25*f:25*(f+1)][j]][['风速', '风向']].values.T, 2)
        v_pred, d_pred = np.vsplit(df_pred_period[df_pred_period['风机'] == machine[25*f:25*(f+1)][j]][['风速', '风向']].values.T, 2)

        value += 0.7 * m(v_pred, v_true) + 0.3 * n(v_true, field[f], d_pred, d_true)
    value /= p
    
    return value

# 公式 3
m = lambda v_pred, v_true: np.sum(np.abs(v_pred - v_true) * np.array([w(k) for k in range(1, 20+1)]))

# 公式 4
n = lambda v_true, field, d_pred, d_true: np.sum(a(v_true, field) * e(d_pred, d_true))

# 公式 5
e = lambda d_pred, d_true: np.minimum((d_pred - d_true) % 1, (d_true - d_pred) % 1)

# 公式 6
w = lambda k: 0.06 if k <= 10 else 0.04

# 公式 7
# 把 alpha 写成 a 是故意的
a = lambda v_true, field: (v_true <= 0.15 if field == '风场1' else 0.086 if field == '风场2' else 0) * np.array([w(k) for k in range(1, 20+1)])

print(R())

结果如下

84.60212382607445
Wall time: 6.14 s

%%time

# 先算再用

w = lambda k: 0.06 if k <= 10 else 0.04

# 两个风场分开计算

df_pred_1, df_pred_2 = [df_pred[df_pred['风场'] == key] for key in df_pred['风场'].unique()]
df_pred_1.index = df_pred_2.index = np.arange(40000)

df_true_1, df_true_2 = [df_true[df_true['风场'] == key] for key in df_true['风场'].unique()]
df_true_1.index = df_true_2.index = np.arange(40000)

# 分开 风速 和 风向 这两列

v_pred_1, d_pred_1 = np.vsplit(df_pred_1[['风速', '风向']].values.T, 2)
v_true_1, d_true_1 = np.vsplit(df_true_1[['风速', '风向']].values.T, 2)

v_pred_2, d_pred_2 = np.vsplit(df_pred_2[['风速', '风向']].values.T, 2)
v_true_2, d_true_2 = np.vsplit(df_true_2[['风速', '风向']].values.T, 2)

# --------------------主要思路--------------------
# 先计算全部的 mj 和 nj
# 根据 公式2 计算同一时段所有风机的平均预测误差
# 再计算同一个风场的所有时段的平均预测误差
# 最后根据 公式1 计算最终得分R
# --------------------------------------------------

# 以下是风场1

m1 = np.sum(np.abs(v_pred_1.reshape(25 * 80, 20) - v_true_1.reshape(25 * 80, 20)) * np.array([w(k) for k in range(1, 20+1)]), -1)

# 计算 a1(alpha_k 风向预测偏差的权重)时要注意阈值，风场1和风场2的阈值不一样
threshold = 0.15
a1 = (v_true_1.reshape(25 * 80, 20) <= threshold) * np.array([w(k) for k in range(1, 20+1)])
e1 = np.minimum((d_pred_1.reshape(25 * 80, 20) - d_true_1.reshape(25 * 80, 20)) % 1, (d_true_1.reshape(25 * 80, 20) - d_pred_1.reshape(25 * 80, 20)) % 1)
n1 = np.sum(a1 * e1, -1)

# reshape后，行 为风机，列 为时段，转置后，就能按公式2来计算了
E1 = np.sum((0.7 * m1 + 0.3 * n1).reshape(25, 80).T.sum(-1) / 25)

# --------------------------------------------------（过程和计算风场1时一样）

# 以下是风场2

m2 = np.sum(np.abs(v_pred_2.reshape(25 * 80, 20) - v_true_2.reshape(25 * 80, 20)) * np.array([w(k) for k in range(1, 20+1)]), -1)

threshold = 0.086
a2 = (v_true_2.reshape(25 * 80, 20) <= threshold) * np.array([w(k) for k in range(1, 20+1)])
e2 = np.minimum((d_pred_2.reshape(25 * 80, 20) - d_true_2.reshape(25 * 80, 20)) % 1, (d_true_2.reshape(25 * 80, 20) - d_pred_2.reshape(25 * 80, 20)) % 1)
n2 = np.sum(a2 * e2, -1)

E2 = np.sum((0.7 * m2 + 0.3 * n2).reshape(25, 80).T.sum(-1) / 25)

# --------------------------------------------------

R = 100 / (1 + E1 + E2)
print(R)

结果如下

84.60305404909366
Wall time: 46 ms

两种方法的计算结果相差不大，有一丢丢误差可能是计算过程中浮点精度的问题