0 回归-海上风电出力预测

汀沿河

已于 2024-04-17 11:13:54 修改

阅读量478

点赞数 2

分类专栏： # 2比赛常用的代码文章标签：回归数据挖掘人工智能

于 2024-04-17 11:09:42 首次发布

本文链接：https://blog.csdn.net/qq_28611929/article/details/137864170

版权

2比赛常用的代码专栏收录该内容

2 篇文章 0 订阅

订阅专栏

https://www.dcic-china.com/competitions/10098

分析一下：特征工程如何做。

时间特征：小时、分钟、一个星期中的第几天、一个月中的第几天。这些可以作为周期特征的标识。比如周六周日的人流会有很大的波动，这些如果不告诉模型它是很难学习到知识的。
业务特征：这方面需要查阅相关的知识点了。操作基本都是在对单个特征特殊处理f（x），两个特征之间做四则运算。同一业务特征做加减，不同领域特征做乘除。最好做出来的特征有实际的物理意义。
历史序列特征：滑动窗口、移动平均等等；我之前参加过一个做的特征工作是爆炸式的，也是惊讶了我，但是别人的结果是真的好。这玩意真有点迷，做尝试吧。
label处理。比如回归，如果能降低当前标签的量纲一定要做。可以与某个及其相关的特征做除法（减法），缩小变化，这样防止模型预测的结果不可控。

import numpy as np
import pandas as pd
import lightgbm as lgb
import xgboost as xgb
from catboost import CatBoostClassifier, CatBoostRegressor
from sklearn.model_selection import StratifiedKFold, KFold, GroupKFold
from sklearn.metrics import mean_squared_error, mean_absolute_error
import matplotlib.pyplot as plt
import tqdm
import sys
import os
import gc
import argparse
import warnings
warnings.filterwarnings('ignore')

# 读取数据
train_info = pd.read_csv('../data/first_data/A榜-训练集_海上风电预测_基本信息.csv', encoding='gbk')
train_df = pd.read_csv('../data/first_data/A榜-训练集_海上风电预测_气象变量及实际功率数据.csv', encoding='gbk')

test_info = pd.read_csv('../data/first_data/B榜-测试集_海上风电预测_基本信息.csv', encoding='gbk')
test_df = pd.read_csv('../data/first_data/B榜-测试集_海上风电预测_气象变量数据.csv', encoding='gbk')

submit_example = pd.read_csv('../data/first_data/submit_example.csv')

train_df = train_df.merge(train_info[['站点编号','装机容量(MW)']], on=['站点编号'], how='left')
test_df = test_df.merge(test_info[['站点编号','装机容量(MW)']], on=['站点编号'], how='left')

train_df['站点编号'] = train_df['站点编号'].apply(lambda x:int(x[1]))
test_df['站点编号'] = test_df['站点编号'].apply(lambda x:int(x[1]))

train_df.columns = ['stationId','time','airPressure','relativeHumidity','cloudiness','10mWindSpeed','10mWindDirection',
                 'temperature','irradiation','precipitation','100mWindSpeed','100mWindDirection','power','capacity']

test_df.columns = ['stationId','time','airPressure','relativeHumidity','cloudiness','10mWindSpeed','10mWindDirection',
                 'temperature','irradiation','precipitation','100mWindSpeed','100mWindDirection','capacity']

# 特征组合
train_df['100mWindSpeed/10mWindSpeed'] = train_df['100mWindSpeed'] / (train_df['10mWindSpeed'] + 0.0000001)
test_df['100mWindSpeed/10mWindSpeed'] = test_df['100mWindSpeed'] / (test_df['10mWindSpeed'] + 0.0000001)

train_df['100mWindDirection/10mWindDirection'] = train_df['100mWindDirection'] / (train_df['10mWindDirection'] + 0.0000001)
test_df['100mWindDirection/10mWindDirection'] = test_df['100mWindDirection'] / (test_df['10mWindDirection'] + 0.0000001)

train_df['10mWindDirection_new'] = train_df['10mWindDirection'] - 180
test_df['10mWindDirection_new'] = test_df['10mWindDirection'] - 180

# 差值
train_df['100mWindSpeed_10mWindSpeed'] = train_df['100mWindSpeed'] - train_df['10mWindSpeed'] 
test_df['100mWindSpeed_10mWindSpeed'] = test_df['100mWindSpeed'] - test_df['10mWindSpeed']

train_df['100mWindDirection_10mWindDirection'] = train_df['100mWindDirection'] - train_df['10mWindDirection']
test_df['100mWindDirection_10mWindDirection'] = test_df['100mWindDirection'] - test_df['10mWindDirection']

# 风切变指数
train_df['WindSpeed/WindDirectio'] = train_df['100mWindSpeed/10mWindSpeed'] / train_df['100mWindDirection/10mWindDirection']
test_df['WindSpeed/WindDirectio'] = test_df['100mWindSpeed/10mWindSpeed'] / test_df['100mWindDirection/10mWindDirection']

train_df['100mWindSpeed/10mWindSpeed_2'] = train_df['100mWindSpeed/10mWindSpeed'].apply(lambda x:np.log10(x)) / 10
test_df['100mWindSpeed/10mWindSpeed_2'] = test_df['100mWindSpeed/10mWindSpeed'].apply(lambda x:np.log10(x)) / 10

# 湿度/温度
train_df['relativeHumidity/temperature'] = train_df['relativeHumidity'] / (train_df['temperature'] + 0.0000001)
test_df['relativeHumidity/temperature'] = test_df['relativeHumidity'] / (test_df['temperature'] + 0.0000001)

# 辐射/温度
train_df['irradiation/temperature'] = train_df['irradiation'] / (train_df['temperature'] + 0.0000001)
test_df['irradiation/temperature'] = test_df['irradiation'] / (test_df['temperature'] + 0.0000001)

# 辐射/云量
train_df['irradiation/cloudiness'] = train_df['irradiation'] / (train_df['cloudiness'] + 0.0000001)
test_df['irradiation/cloudiness'] = test_df['irradiation'] / (test_df['cloudiness'] + 0.0000001)

# 是否降水
train_df['is_precipitation'] = train_df['precipitation'].apply(lambda x:1 if x>0 else 0)
test_df['is_precipitation'] = test_df['precipitation'].apply(lambda x:1 if x>0 else 0)

def get_time_feature(df, col):
    
    df_copy = df.copy()
    prefix = col + "_"
    df_copy[col] = df_copy[col].astype(str)
    
    df_copy[col] = pd.to_datetime(df_copy[col], format='%Y-%m-%d %H:%M')
    df_copy[prefix + 'month'] = df_copy[col].dt.month
    df_copy[prefix + 'day'] = df_copy[col].dt.day
    df_copy[prefix + 'hour'] = df_copy[col].dt.hour
    df_copy[prefix + 'minute'] = df_copy[col].dt.minute
    df_copy[prefix + 'weekofyear'] = df_copy[col].dt.weekofyear
    df_copy[prefix + 'dayofyear'] = df_copy[col].dt.dayofyear
    
    return df_copy   

train_df = get_time_feature(train_df, 'time')
test_df = get_time_feature(test_df, 'time')


# 合并训练数据和测试数据
train_df['is_test'] = 0
test_df['is_test'] = 1
df = pd.concat([train_df, test_df], axis=0).reset_index(drop=True)

# 构建特征
num_cols = ['airPressure','relativeHumidity','cloudiness','10mWindSpeed','10mWindDirection',
            'temperature','irradiation','precipitation','100mWindSpeed','100mWindDirection']

for col in tqdm.tqdm(num_cols):
    # 历史平移/差分特征
    for i in [1,2,3,4,5,6,7,15,30,50] + [1*96,2*96,3*96,4*96,5*96]:
        df[f'{col}_shift{i}'] = df.groupby('stationId')[col].shift(i)
        df[f'{col}_feture_shift{i}'] = df.groupby('stationId')[col].shift(-i)

        df[f'{col}_diff{i}'] = df[f'{col}_shift{i}'] - df[col]
        df[f'{col}_feture_diff{i}'] = df[f'{col}_feture_shift{i}'] - df[col]

        df[f'{col}_2diff{i}'] = df.groupby('stationId')[f'{col}_diff{i}'].diff(1)
        df[f'{col}_feture_2diff{i}'] = df.groupby('stationId')[f'{col}_feture_diff{i}'].diff(1)
    
    # 均值相关
    df[f'{col}_3mean'] = (df[f'{col}'] + df[f'{col}_feture_shift1'] + df[f'{col}_shift1'])/3
    df[f'{col}_5mean'] = (df[f'{col}_3mean']*3 + df[f'{col}_feture_shift2'] + df[f'{col}_shift2'])/5
    df[f'{col}_7mean'] = (df[f'{col}_5mean']*5 + df[f'{col}_feture_shift3'] + df[f'{col}_shift3'])/7
    df[f'{col}_9mean'] = (df[f'{col}_7mean']*7 + df[f'{col}_feture_shift4'] + df[f'{col}_shift4'])/9
    df[f'{col}_11mean'] = (df[f'{col}_9mean']*9 + df[f'{col}_feture_shift5'] + df[f'{col}_shift5'])/11
    
    df[f'{col}_shift_3_96_mean'] = (df[f'{col}_shift{1*96}'] + df[f'{col}_shift{2*96}'] + df[f'{col}_shift{3*96}'])/3
    df[f'{col}_shift_5_96_mean'] = (df[f'{col}_shift_3_96_mean']*3 + df[f'{col}_shift{4*96}'] + df[f'{col}_shift{5*96}'])/5
    df[f'{col}_future_shift_3_96_mean'] = (df[f'{col}_feture_shift{1*96}'] + df[f'{col}_feture_shift{2*96}'] + df[f'{col}_feture_shift{3*96}'])/3
    df[f'{col}_future_shift_5_96_mean'] = (df[f'{col}_future_shift_3_96_mean']*3 + df[f'{col}_feture_shift{4*96}'] + df[f'{col}_feture_shift{5*96}'])/3
    
    # 窗口统计
    for win in [3,5,7,14,28]:
        df[f'{col}_win{win}_mean'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').mean().values
        df[f'{col}_win{win}_max'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').max().values
        df[f'{col}_win{win}_min'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').min().values
        df[f'{col}_win{win}_std'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').std().values
        df[f'{col}_win{win}_skew'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').skew().values
        df[f'{col}_win{win}_kurt'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').kurt().values
        df[f'{col}_win{win}_median'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').median().values
        
        df = df.sort_values(['stationId','time'], ascending=False)
        
        df[f'{col}_future_win{win}_mean'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').mean().values
        df[f'{col}_future_win{win}_max'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').max().values
        df[f'{col}_future_win{win}_min'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').min().values
        df[f'{col}_future_win{win}_std'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').std().values
        df[f'{col}_future_win{win}_skew'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').skew().values
        df[f'{col}_future_win{win}_kurt'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').kurt().values
        df[f'{col}_future_win{win}_median'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').median().values
        
        df = df.sort_values(['stationId','time'], ascending=True)
        
        # 二阶特征
        df[f'{col}_win{win}_mean_loc_diff'] = df[col] - df[f'{col}_win{win}_mean']
        df[f'{col}_win{win}_max_loc_diff'] = df[col] - df[f'{col}_win{win}_max']
        df[f'{col}_win{win}_min_loc_diff'] = df[col] - df[f'{col}_win{win}_min']
        df[f'{col}_win{win}_median_loc_diff'] = df[col] - df[f'{col}_win{win}_median']
        
        df[f'{col}_future_win{win}_mean_loc_diff'] = df[col] - df[f'{col}_future_win{win}_mean']
        df[f'{col}_future_win{win}_max_loc_diff'] = df[col] - df[f'{col}_future_win{win}_max']
        df[f'{col}_future_win{win}_min_loc_diff'] = df[col] - df[f'{col}_future_win{win}_min']
        df[f'{col}_future_win{win}_median_loc_diff'] = df[col] - df[f'{col}_future_win{win}_median']
        
for col in ['is_precipitation']:
    for win in [4,8,12,20,48,96]:
        df[f'{col}_win{win}_mean'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').mean().values
        df[f'{col}_win{win}_sum'] = df.groupby('stationId')[col].rolling(window=win, min_periods=3, closed='left').sum().values


train_df = df[df.is_test==0].reset_index(drop=True)
test_df = df[df.is_test==1].reset_index(drop=True)
del df
gc.collect()

train_df = train_df[train_df['power']!='<NULL>'].reset_index(drop=True)
train_df['power'] = train_df['power'].astype(float)

cols = [f for f in test_df.columns if f not in ['time','power','is_test']] # capacity
def cv_model(clf, train_x, train_y, test_x, capacity, seed=2024):
    folds = 5
    kf = KFold(n_splits=folds, shuffle=True, random_state=seed)
    oof = np.zeros(train_x.shape[0])
    test_predict = np.zeros(test_x.shape[0])
    cv_scores = []
    
    for i, (train_index, valid_index) in enumerate(kf.split(train_x, train_y)):
        print('************************************ {} ************************************'.format(str(i+1)))
        trn_x, trn_y, val_x, val_y = train_x.iloc[train_index], train_y[train_index], train_x.iloc[valid_index], train_y[valid_index]
        
        # 转化目标,进行站点目标归一化
        trn_y = trn_y / capacity[train_index]
        val_y = val_y / capacity[valid_index]
        
        train_matrix = clf.Dataset(trn_x, label=trn_y)
        valid_matrix = clf.Dataset(val_x, label=val_y)
        params = {
            'boosting_type': 'gbdt',
            'objective': 'regression',
            'metric': 'rmse',
            'min_child_weight': 5,
            'num_leaves': 2 ** 8,
            'lambda_l2': 10,
            'feature_fraction': 0.8,
            'bagging_fraction': 0.8,
            'bagging_freq': 4,
            'learning_rate': 0.1,
            'seed': 2023,
            'nthread' : 16,
            'verbose' : -1,
        }
        model = clf.train(params, train_matrix, 3000, valid_sets=[train_matrix, valid_matrix],
                          categorical_feature=[], verbose_eval=500, early_stopping_rounds=200)
        val_pred = model.predict(val_x, num_iteration=model.best_iteration)
        test_pred = model.predict(test_x, num_iteration=model.best_iteration)
        
        oof[valid_index] = val_pred
        test_predict += test_pred / kf.n_splits
        
        score = 1/(1+np.sqrt(mean_squared_error(val_pred * capacity[valid_index], val_y * capacity[valid_index])))
        cv_scores.append(score)
        print(cv_scores)
        
        if i == 0:
            imp_df = pd.DataFrame()
            imp_df["feature"] = cols
            imp_df["importance_gain"] = model.feature_importance(importance_type='gain')
            imp_df["importance_split"] = model.feature_importance(importance_type='split')
            imp_df["mul"] = imp_df["importance_gain"]*imp_df["importance_split"]
            imp_df = imp_df.sort_values(by='mul',ascending=False)
            imp_df.to_csv('feature_importance.csv', index=False)
            print(imp_df[:30])
            
    return oof, test_predict

lgb_oof, lgb_test = cv_model(lgb, train_df[cols], train_df['power'], test_df[cols], train_df['capacity'])