电力需求预测挑战赛｜Datawhale AI夏令营第二期｜代码及笔记分享

最新分数：215.14567（13/1605）

赛题简介

赛事链接：电力需求预测挑战赛（科大讯飞xDatawhale）
给定多个房屋对应电力消耗历史N天的相关序列数据等信息，预测房屋对应电力的消耗。主要特征包括房屋id，日标识dt（1为数据集最近1天），房屋类型type，预测目标电力消耗target。本赛题是一个典型的多变量时间序列问题。当下时间序列预测的方法主要有三种：
1.传统的时间序列预测方法，ARIMA和指数平滑法
2.基于机器学习的方法，lightgbm、xgboost、catboost
3.基于深度学习的方法，RNN、LSTM、Wavenet

数据读取

导入必要的库，pandas、numpy、scipy等用于数据读取和计算，sklearn、lightgbm、xgboost、catboost等用于模型训练，tqdm用于处理进程的可视化展示，joblib用于模型的保存和加载，以及其他涉及系统变量的库。
训练集共有2877305行数据，时间范围为11-506不等。
测试集共有58320行数据，时间范围为1-10。

探索性数据分析

基本信息
定义data_stats函数展示数据集的基本信息情况，包括属性个数、最大属性占比、特征类型。

def data_stats(data):
    stats = []
    for col in data.columns:
        stats.append((col, data[col].nunique(),
                      round(data[col].value_counts(normalize = True, dropna = False).values[0] * 100, 3),
                      data[col].dtype))
        stats_df = pd.DataFrame(stats, columns = ['特征','属性个数','最大属性占比','特征类型'])
    return stats_df

data_stats(train)

	特征	属性个数	最大属性占比	特征类型
0	id	5832	0.017	object
1	dt	496	0.203	int64
2	type	19	20.722	int64
3	target	181304	0.004	float64

data_stats(test)

	特征	属性个数	最大属性占比	特征类型
0	id	5832	0.017	object
1	dt	10	10.000	int64
2	type	19	20.645	int64

不同type类型对应target的柱状图
可以发现不同type对应的电力需求存在比较大的差异，因此type应该是模型训练中的一个重要特征，可以根据type构建多个特征。

type_target_df = train.groupby('type')['target'].mean().reset_index()
plt.figure(figsize=(8, 4))
plt.bar(type_target_df['type'], type_target_df['target'], color=['#5ba2e6', '#f89732'])
plt.xlabel('Type')
plt.ylabel('Average Target Value')
plt.title('Bar Chart of Target by Type')
plt.show()

不同id的target变化趋势折线图
可以看到不同id的电力需求变化趋势存在着很大的差异，因此对于id也可以尝试进行标签编码作为类别特征用于模型训练。

unique_ids = train['id'].unique()
fig, axes = plt.subplots(4, 4, figsize=(20, 16), sharex=True, sharey=True)
axes = axes.flatten()

for i, unique_id in enumerate(unique_ids[:16]): 
    specific_id_df = train[train['id'] == unique_id]
    axes[i].plot(specific_id_df['dt'], specific_id_df['target'], linestyle='-', color='#5ba2e6')
    axes[i].set_title(f"ID: {unique_id}")
    axes[i].set_xlabel('DateTime')
    axes[i].set_ylabel('Target Value')

plt.tight_layout()
plt.show()

特征工程

特征构造

历史平移特征

提取具有相似性的信息，将t-1、t-2、…、t-n个时间单位的值用作第t个时间单位的特征，时间区间为一个月左右，也可以根据数据的周期性来确定。

for i in tqdm(list(range(10,36))+[2*30,3*30,4*30,5*30,6*30,7*30,8*30,9*30,10*30,11*30,12*30]):
        data[f'target_shift{i}'] = data.groupby('id')['target'].shift(i)
    
data['3_mean_target'] = (data['target_shift10'] + data['target_shift11'] + data['target_shift12']) / 3
data['5_mean_target'] = (data['3_mean_target']*3 + data['target_shift13'] + data['target_shift14']) / 5
data['7_mean_target'] = (data['5_mean_target']*5 + data['target_shift15'] + data['target_shift16']) / 7

历史平移+差分特征

可以帮助获取相邻阶段的增长差异，描述数据的涨减变化情况。

for i in tqdm(list(range(10,36))+[2*30,3*30,4*30,5*30,6*30,7*30,8*30,9*30,10*30,11*30,12*30]):
    for d in range(1,4):
        data[f'target_shift{i}_diff{d}'] = data.groupby('id')[f'target_shift{i}'].diff(d)

窗口统计特征

窗口统计可以构建不同的窗口大小，然后基于窗口范围进统计均值、最大值、最小值、中位数、方差的信息，可以反映最近阶段数据的变化情况。

for win in tqdm([7,10,14,21,28,30,35,50,70,90,120,150,180,240,360]):
    data[f'target_win{win}_mean'] = data.groupby('id')['target'].rolling(window=win, min_periods=3, closed='left').mean().values
    data[f'target_win{win}_max'] = data.groupby('id')['target'].rolling(window=win, min_periods=3, closed='left').max().values
    data[f'target_win{win}_min'] = data.groupby('id')['target'].rolling(window=win, min_periods=3, closed='left').min().values
    data[f'target_win{win}_std'] = data.groupby('id')['target'].rolling(window=win, min_periods=3, closed='left').std().values

历史平移+窗口统计特征

通过历史平移获取上个阶段的信息。

for i in [10]:
    for win in tqdm([7,10,14,21,28,30,35,50,70,90,120,150,180,240,360]):
        data[f'target_shift{i}_win{win}_mean'] = data.groupby('id')[f'target_shift{i}'].rolling(window=win, min_periods=3, closed='left').mean().values
        data[f'target_shift{i}_win{win}_max'] = data.groupby('id')[f'target_shift{i}'].rolling(window=win, min_periods=3, closed='left').max().values
        data[f'target_shift{i}_win{win}_min'] = data.groupby('id')[f'target_shift{i}'].rolling(window=win, min_periods=3, closed='left').min().values
        data[f'target_shift{i}_win{win}_sum'] = data.groupby('id')[f'target_shift{i}'].rolling(window=win, min_periods=3, closed='left').sum().values
        data[f'target_shift{i}_win{win}_std'] = data.groupby('id')[f'target_shift{i}'].rolling(window=win, min_periods=3, closed='left').std().values

自动化特征构造

OpenFE是一个自动化特征生成的框架，可以高效的构建有效的新特征，显著地提升GBDT（例如LightGBM，XGBoost等）和各种SOTA神经网络（例如Transformer，AutoInt，TabNet等）在表格类数据的效果。

from openfe import OpenFE, tree_to_formula, transform, TwoStageFeatureSelector
ofe = OpenFE()
ofe_features = ofe.fit(data=train_x, label=train_y,
                       metric='rmse', categorical_features=['type'],
                       n_data_blocks=8, feature_boosting=True,
                       n_jobs=8, task='regression',verbose=False)
fs = TwoStageFeatureSelector(metric='rmse', categorical_features=['type'],n_jobs=8)
features = fs.fit(data=train_x, label=train_y)

根据上述的自动化特征构造和选择，可以构建如下的特征：

data['GroupByThenRank_id_type'] = data.groupby('type')['id'].rank(ascending=True, pct=True).values
data['GroupByThenRank_dt_type'] = data.groupby('type')['dt'].rank(ascending=True, pct=True).values
data['log_dt'] = np.log(np.abs(data.dt.replace(0, np.nan)))

相关性分析

相关性分析被广泛应用于评估特征与目标变量之间的相关性。通过计算相关系数或其他相关性度量，我们可以筛选出与目标变量相关性较高的特征。设定合适的阈值，我们可以选择保留相关性高于阈值的特征，从而减少特征空间的维度，提高模型的训练效率。

train_se = data_fe[data_fe.target.notnull()].reset_index(drop=True)
features = train_se.columns[1:]
corr = []
for feat in features:
    corr.append(abs(train_se[[feat, 'target']].fillna(0).corr().values[0][1]))

se = pd.Series(corr, index=features).sort_values(ascending=False)

模型训练

# 进行数据切分
train = data_fe[data_fe.target.notnull()].reset_index(drop=True)
test = data_fe[data_fe.target.isnull()].reset_index(drop=True)

# 确定输入特征
train_cols = [f for f in data_fe.columns if f not in ['id','target']]

交叉验证

关于时间序列数据的交叉验证，能够直接使用且不打乱时间顺序的有TimeSeriesSplit、K-Fold方法，此外还有MonteCarloCV、GroupTimeSeriesSplit两种自定义的方法。
TimeSeriesSplit：时间序列被分成几个大小相等的折叠，然后每一次折首先被用来测试一个模型，然后重新训练它，除了第一折只用于训练。
K-Fold：将原始数据集分成K个子集，称为折（fold）。 然后，模型在K个不同的训练集上进行K次训练和验证，在每次训练中，其中一个折被作为验证集，而剩下的K-1个折被用作训练集。在某些情况下，K-fold交叉验证对时间序列是有用的。
MonteCarloCV：训练集的大小在每次迭代过程中都是固定的，验证原点是随机选择的，这个原点标志着训练集的结束和验证的开始。
GroupTimeSeriesSplit：首先计算一下unique的个数，根据group的分组进行相关的按照时间的组合，每次平移一个group，其中前几个group对应的数据为训练集，而紧接着时间后的一个group的数据为验证集。
在本数据集中，由于整个数据的时间序列不是连续的，只是每个id对应的时间序列是连续的，因此保留观测的时间顺序反而不是很有效，甚至会降低分数，在这里依旧适用K-Fold的方法，使用shuffle=True来打乱顺序。

from sklearn.model_selection import StratifiedKFold, KFold, GroupKFold
import lightgbm as lgb
import xgboost as xgb
from catboost import CatBoostRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error
def cv_model(clf, train_x, train_y, test_x, clf_name, seed = 2024):
    '''
    clf：调用模型
    train_x：训练数据
    train_y：训练数据对应标签
    test_x：测试数据
    clf_name：选择使用模型名
    seed：随机种子
    '''
    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]
        
        if clf_name == "lgb":
            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': 'mae',
                'min_child_weight': 6,
                'num_leaves': 2 ** 6,
                'lambda_l2': 10,
                'feature_fraction': 0.8,
                'bagging_fraction': 0.8,
                'bagging_freq': 4,
                'learning_rate': 0.1,
                'seed': 2024,
                'nthread' : -1
            }
            callback = [lgb.early_stopping(stopping_rounds=100, verbose=100)]
            model = clf.train(params, train_matrix, 1000, valid_sets=[train_matrix, valid_matrix],
                              categorical_feature=['type'], callbacks = callback)
            val_pred = model.predict(val_x, num_iteration=model.best_iteration)
            test_pred = model.predict(test_x, num_iteration=model.best_iteration)
            joblib.dump(model,"lgbmodel_157_5.pkl")
        
        if clf_name == "xgb":
            xgb_params = {
              'booster': 'gbtree', 
              'objective': 'reg:squarederror',
              'eval_metric': 'mae',
              'max_depth': 5,
              'lambda': 10,
              'subsample': 0.7,
              'colsample_bytree': 0.7,
              'colsample_bylevel': 0.7,
              'eta': 0.1,
              'tree_method': 'hist',
              'seed': 2024,
              'nthread': -1
              }
            train_matrix = clf.DMatrix(trn_x , label=trn_y)
            valid_matrix = clf.DMatrix(val_x , label=val_y)
            test_matrix = clf.DMatrix(test_x)
            
            watchlist = [(train_matrix, 'train'),(valid_matrix, 'eval')]
            
            model = clf.train(xgb_params, train_matrix, num_boost_round=1000, evals=watchlist, verbose_eval=100, early_stopping_rounds=100)
            val_pred  = model.predict(valid_matrix)
            test_pred = model.predict(test_matrix)
            joblib.dump(model,"xgbmodel_157_5.pkl")
            
        if clf_name == "cat":
            params = {'learning_rate': 0.1, 'depth': 5, 'bootstrap_type':'Bernoulli',
                      'thread_count':-1, 'od_type': 'Iter', 'od_wait': 100, 
                      'random_seed': 2024, 'allow_writing_files': False}
            
            model = clf(iterations=1000, **params)
            model.fit(trn_x, trn_y, eval_set=(val_x, val_y),
                      use_best_model=True, 
                      metric_period=100,
                      cat_features=['type'],
                      verbose=1)
            
            val_pred  = model.predict(val_x)
            test_pred = model.predict(test_x)
            joblib.dump(model,"catmodel_165_5.pkl")
        
        oof[valid_index] = val_pred
        test_predict += test_pred / kf.n_splits
        
        score = mean_absolute_error(val_y, val_pred)
        cv_scores.append(score)
        print(cv_scores)
        
    return oof, test_predict

# 选择lightgbm模型
lgb_oof, lgb_test = cv_model(lgb, train[train_cols], train['target'], test[train_cols], 'lgb')
# 选择xgboost模型
xgb_oof, xgb_test = cv_model(xgb, train[train_cols], train['target'], test[train_cols], 'xgb')
# 选择catboost模型
cat_oof, cat_test = cv_model(CatBoostRegressor, train[train_cols], train['target'], test[train_cols], 'cat')

# 进行取平均融合
final_test = (lgb_test + xgb_test + cat_test) / 3

模型融合

利用多个基学习器学习原数据，然后将这几个基学习学习到的数据交给第二层模型进行拟合。说白了就是将第一层模型的输出作为第二层模型的输入。

from sklearn.linear_model import Ridge
def stack_model(oof_1, oof_2, oof_3, predictions_1, predictions_2, predictions_3, y):
    '''
    输入的oof_1, oof_2, oof_3可以对应lgb_oof，xgb_oof，cat_oof
    predictions_1, predictions_2, predictions_3对应lgb_test，xgb_test，cat_test
    '''
    train_stack = pd.concat([oof_1, oof_2, oof_3], axis=1)
    test_stack = pd.concat([predictions_1, predictions_2, predictions_3], axis=1)
    
    oof = np.zeros((train_stack.shape[0],))
    predictions = np.zeros((test_stack.shape[0],))
    scores = []
    
    from sklearn.model_selection import RepeatedKFold
    folds = RepeatedKFold(n_splits=5, n_repeats=2, random_state=2021)
    
    for fold_, (trn_idx, val_idx) in enumerate(folds.split(train_stack, train_stack)): 
        print("fold n°{}".format(fold_+1))
        trn_data, trn_y = train_stack.loc[trn_idx], y[trn_idx]
        val_data, val_y = train_stack.loc[val_idx], y[val_idx]
        
        clf = Ridge(random_state=2021)
        clf.fit(trn_data, trn_y)

        oof[val_idx] = clf.predict(val_data)
        predictions += clf.predict(test_stack) / (5 * 2)
        
        score_single = mean_absolute_error(val_y, oof[val_idx])
        scores.append(score_single)
        print(f'{fold_+1}/{5}', score_single)
    print('mean: ',np.mean(scores))
   
    return oof, predictions

stack_oof, stack_pred = stack_model(pd.DataFrame(lgb_oof), pd.DataFrame(xgb_oof), pd.DataFrame(cat_oof), 
                                    pd.DataFrame(lgb_test), pd.DataFrame(xgb_test), pd.DataFrame(cat_test), train['target'])

Optuna超参数调优

def cat_optuna(trial,data=data,target=target):

    param = {
        'loss_function': 'RMSE',
        'task_type': 'GPU',
        'random_state': 2024,
        'n_estimators':  1000,
        'od_type': 'Iter',
        'od_wait': 50,
        'use_best_model': True,
        'l2_leaf_reg': trial.suggest_float('l2_leaf_reg', 1e-3, 10.0),
        'learning_rate': trial.suggest_float('learning_rate', 0.001, 1.0),
        'max_depth': trial.suggest_int('max_depth', 1, 16),
        'bootstrap_type': trial.suggest_categorical('bootstrap_type', ['Bayesian', 'Bernoulli', 'MVS']),
        'min_data_in_leaf': trial.suggest_int('min_data_in_leaf', 1, 100),
    }

    if param["bootstrap_type"] == "Bayesian":
        param["bagging_temperature"] = trial.suggest_float("bagging_temperature", 0, 10)
    elif param["bootstrap_type"] == "Bernoulli":
        param["subsample"] = trial.suggest_float("subsample", 0.1, 1, log=True)

    kf = KFold(n_splits=5, shuffle=True, random_state=2024)
    catcv_scores=[]
    for i, (train_index, valid_index) in enumerate(kf.split(data, target)):
        print('************************************ {} ************************************'.format(str(i+1)))
        trn_x, trn_y, val_x, val_y = data.iloc[train_index], target[train_index], data.iloc[valid_index], target[valid_index]
        
        model = CatBoostRegressor(**param)
        model.fit(trn_x, trn_y, eval_set=(val_x, val_y), cat_features=['type','id'],
                  verbose=0, early_stopping_rounds=50)
        val_preds = model.predict(val_x)
        rmse = mean_squared_error(val_y, val_preds)
        catcv_scores.append(rmse)

    return sum(catcv_scores) / len(catcv_scores)

optuna.visualization.plot_optimization_history(study)

optuna.visualization.plot_param_importances(study)

optuna.visualization.plot_edf(study)

参考资料

1.机器学习算法竞赛实战
2.时间序列数据的特征工程总结
3.OpenFE
4.时间序列的蒙特卡罗交叉验证
5.9个时间序列交叉验证方法的介绍和对比
6.时间序列特有的交叉验证方法GroupTimeSeriesSplit