Kaggle参戦記：エストニアの電力需給予測にいっちょかみで参加

こんばんは！「人間の文章かAIの文章化を見分けるコンペ（LLM detect）」は残念な結果に終わってしまいましたが、その後も懲りずに新たなコンペにチャレンジしています。

今度は、「Enefit」という、エストニアの電力会社が主催する「ソーラーパネルを設置しているプロシューマーの電力需給予測」のコンペに参加してみました。LLM detectが終わってからの参加だったので、その時点ですでに「締め切り8日間前」と、かなり切羽詰まった状況だったのですが、例によってパブリックに共有されているノートブックを参考にキャッチアップ＆性能向上して、パブリックLBでは2731チーム中232位(上位8%)で終了。なんとか銅メダル圏内に位置しています。前回同様、また大規模Shake-downの可能性もあるため、ぬか喜びはできませんが...

現在は、実データを基にしたプライベートLBスコアを算出中です。4月末までのデータを対象に、プライベートLBスコアが決まりますので、最終結果が出るのはまだしばらく先です。また結果が出たらご報告します。

Enefitのコンペを通じて学んだこと

今回のコンペでは、主に下記の2つのノートブックを参考にさせていただきました。

• 「事前にトレーニングされたモデルを使った予測」
• 「CatBoostを使ったベースライン&クロスバリデーション」

上記のノートブックを通じて、
(1) 事前トレーニングされたモデルを使った予測
(2) ノートブックのモジュール化（特徴量抽出ロジックの共通化）
(3) Optunaを使ったハイパーパラメータの最適化
のテクニックを学ぶことができました。

(1) 事前トレーニングされたモデルを使った予測

「事前にトレーニングされたモデルを使った予測」のノートブックでは、前半のセルでデータを読み込み、予測に必要な特徴量を生成するロジックが定義されています。そして、「Loading Trained Models」のセクションで、おもむろに事前学習されたモデルを読み込んでいます。

.py

c1 = load('/kaggle/input/enefit-trained-model/voting_regressor_consumption_model.joblib')
p1 = load('/kaggle/input/enefit-trained-model/voting_regressor_production_model.joblib')

このコンペでは、「電力消費量」と「電力供給量」を予測するタスクになっているため、「電力消費量予測用のモデル」がc1、「電力供給量予測用のモデル」がp1になっています。そして、予測のための関数が下記のように定義されています。

.py

def predict(df_features,model_consumption=c1,model_production=p1):
    predictions = np.zeros(len(df_features))

    mask = df_features["is_consumption"] == 1
    predictions[mask.values] = model_consumption.predict(
            df_features[mask]
    ).clip(0)

    mask = df_features["is_consumption"] == 0
    predictions[mask.values] = model_production.predict(
            df_features[mask]
    ).clip(0)

    return predictions

特徴量のis_consumptionが1なら、「電力消費量を予測」するためにモデルc1を使い、is_consumptionが0なら、「電力供給量を予測」するためにモデルp1が使われていますね。Kaggleのノートブックの処理時間の制限は「9時間以内」であることが多いようですが、このようにして、時間のかかる学習処理を切り離して、「事前学習したモデルを使って予測するだけのノートブック」を作ってsubmitもできるわけですね。

自分で事前学習させたモデルを読み込ませるには、いくつか方法があります。1つは、モデルを学習してファイルに出力するようなノートブックを作成し、そのノートブックを「Add Data」⇒「Your Notebooks」でInput Dataとして取り込めるようにする方法。この方法だと、モデルを学習させて、新しいモデルができたら、常に新しいモデルを使うような設定にすることも可能です。

2つ目の方法は、モデルを学習して出力したファイルを、一度ダウンロードして、予測用のノートブックにアップロードする方法。こちらは、「Add Data」ボタンの横の「↑（Upload Data）」ボタンでアップロードできます。私の場合、いろんなパラメータを変えながら、様々なモデルを作って、組み合わせたかったので、後者の方法を採用しました。

元のノートブックでは、is_consumptionが1か0で、消費量と供給量を個別に予測するような作りになっていましたが、他にもis_businessという特徴量があって、「一般家庭」か「企業」かを区別するセグメントもあったので、私は、is_consumption × is_businessの組み合わせで4種類のモデルを学習させるように改造しました。

(2) ノートブックのモジュール化

(1)事前トレーニングされたモデルを使った予測ができるようになると、学習に使った特徴量抽出ロジックと、予測に使う特徴量抽出ロジックを共通化したい、という気持ちが高まってきました。共通化していないと、学習に使う特徴量生成ロジックを変更すると、それと同じように予測側の特徴量生成ロジックも変更しなければならず、同じコードを2重管理する必要が出てきます。

これを解決するためには、特徴量抽出ロジックをライブラリ化することですが、これも比較的簡単にできることが分かりました。Kaggleノートブックの「File」メニューに「Set as Utility Script」を選択すると、ノートブック形式ではなく、.py形式で保存が可能になります。

例えば、特徴量抽出ロジックを「feature_generation.py」という名前で保存したとします。すると、学習用・予測用のノートブックの「File」メニューの「Add Utility Script」から「feature_generation.py」を選択して取り込めるようになります。取り込んだら、最初の方のセルで、

.py

import feature_generation as fg

として、importすれば、自作の特徴量抽出ロジックを、学習用・予測用に共通して使えるようになります。

ここで一点ご注意。共通化された特徴量抽出ロジックを修正／更新すると、学習・予測用のノートブックで、「更新されたライブラリを読み込みますか？」というポップアップが出ます。しかし、そこで「OK」ボタンを押すだけだと、実行時に修正／更新された内容が反映されません。Kaggleのノートブックに限らず、JupyterやVSCでノートブックを実行する時によくある話ですが、ライブラリの修正後の内容を反映して実行するには、一度、カーネルをリセットしてimportしなおす必要があります。私はこれを忘れていて、ちょっとハマりました。

ちなみに、ノートブックの最初にautoreloadを使うと、カーネルをリセットしなくても、importしなおしてくれるらしいです。

.py

%load_ext autoreload
%autoreload 2

(3)Optunaを使ったハイパーパラメータの最適化

これは、知っている人にとっては今更感がありますが、私は今回初めてOptunaを使って、その手軽さに感動しましたので、改めてご紹介します。

「CatBoostを使ったベースライン&クロスバリデーション」での使い方が、とても分かりやすかったです。

キモは、下記の部分。

.py

# loss function
def objective(trial):
    config = {
        'n_estimators': trial.suggest_int('n_estimators', 50, 250),
        'learning_rate': trial.suggest_float('learning_rate', 5e-3, 0.75,log=True),
        'depth': trial.suggest_int('depth', 2, 10, log=True),
        'l2_leaf_reg': trial.suggest_float('l2_leaf_reg',1e-8,100,log=True),
        'model_size_reg': trial.suggest_float('model_size_reg',1e-8,100,log=True),
        'colsample_bylevel': trial.suggest_float("colsample_bylevel", 0.1, 1),
        'subsample': trial.suggest_float("subsample", 0.5, 1)
    }
    
    cv = TimeSeriesSplit(n_splits=3)
    cv_mae = [None]*3
    for i, (train_index, test_index) in enumerate(cv.split(timesteps)):
        cv_mae[i] = fit_and_test_fold(config, X, y, timesteps[train_index], timesteps[test_index])
        
    # saving the individual fold holdout metrics 
    # uncomment this line if you don't want this
    trial.set_user_attr('split_mae', cv_mae)
        
    return np.mean(cv_mae)


sampler = optuna.samplers.TPESampler(
    n_startup_trials=10, seed=1234
)

study = optuna.create_study(
    directions=['minimize'],sampler=sampler,study_name='catboost'
)

# maximum of 50 trials or 2 hr wall clock time
study.optimize(objective, n_trials=50, timeout= 7200) 
_ = joblib.dump(study, 'catboost_enefit_hyperopt.pkl')

下の行から順に見ていきましょう。

.py

# maximum of 50 trials or 2 hr wall clock time
study.optimize(objective, n_trials=50, timeout= 7200) 
_ = joblib.dump(study, 'catboost_enefit_hyperopt.pkl')

studyというオブジェクトを作って、50回または7200秒（2時間）、objectiveの値を最適化しなさい、という命令があります。で、最適化が終わったら、その結果をpklファイルに保存しています。

それでは、objectiveって何でしょうか。こんな感じで定義されていますね。

.py

def objective(trial):
    config = {
        'n_estimators': trial.suggest_int('n_estimators', 50, 250),
        'learning_rate': trial.suggest_float('learning_rate', 5e-3, 0.75,log=True),
        'depth': trial.suggest_int('depth', 2, 10, log=True),
        'l2_leaf_reg': trial.suggest_float('l2_leaf_reg',1e-8,100,log=True),
        'model_size_reg': trial.suggest_float('model_size_reg',1e-8,100,log=True),
        'colsample_bylevel': trial.suggest_float("colsample_bylevel", 0.1, 1),
        'subsample': trial.suggest_float("subsample", 0.5, 1)
    }
    
    cv = TimeSeriesSplit(n_splits=3)
    cv_mae = [None]*3
    for i, (train_index, test_index) in enumerate(cv.split(timesteps)):
        cv_mae[i] = fit_and_test_fold(config, X, y, timesteps[train_index], timesteps[test_index])
        
    # saving the individual fold holdout metrics 
    # uncomment this line if you don't want this
    trial.set_user_attr('split_mae', cv_mae)
        
    return np.mean(cv_mae)

configは、このソースコードの場合、CatBoostRegressorに設定するハイパーパラメータです。trialというオブジェクトを使って、「ここら辺に最適解がありそうかな？」と思う範囲でパラメータを振ってもらうように指定しています。で、時系列の予測なのでTimeSeriesSplitで3通りの異なる期間でクロスバリデーションして、結果をMAE（Mean Absolute Error）で評価し、3期間のMAEの平均値を最小化するように最適化を行う、という処理になっています。

あとは、最適化が終わったら、下記のコードで、平均MAEが小さい順のパラメータの一覧を確認できます。

.py

# get results
results = study.trials_dataframe(attrs=('number','value','duration', 'params'))
results = results.rename(columns={'value':'mae'})
results['duration'] = results['duration']/np.timedelta64(1, 's')
results = results.sort_values(by='mae',ascending=True)
results.to_csv(f'mae_catboost.csv',index=False)
results.head(10)

結果はこんな感じ。22回目の試行で平均MAEが約68で最小値をとって、その時のパラメータが、これこれこうだった。というのが読み取れますね。

そして、最適化されたベストなパラメータで、再度学習しなおして、出来上がったモデルをファイルに保存。これで、「最適なパラメータで事前学習されたモデル」のいっちょあがりです。

.py

model = fit_model(X, y, n_jobs=4, config=study.best_params, verbose=20)
model.save_model('catboost_energy_pred.cbm',format='cbm')

私のノートブックでは、下記のようなちょっとした改造をほどこして使いました。

• 最適化にかける回数／時間の制約を50回／2時間 ⇒ 1万回／9時間に目いっぱい拡大
• 平均MAEではなく、最大MAEを最小化するように最適化（「大きく外さない」という発想）
• CatBoostのParameter Tuningを参照し、n_estimators、depthなどの範囲を再調整

最適なパラメータを見つけるためのTPESamplerというアルゴリズムがなかなか素晴らしく、バックグラウンドで走らせておくと、みるみる最適値に近づいていってくれるので、見ていて気持ちが良かったです。TPEのようなベイズ最適化に基づく最適化の仕組みは、このページのぐりぐり動くグラフが、イメージしやすかったです。

また、ここでは、CatBoostのパラメータの最適化の例について説明しましたが、他にも異なるモデル同士のアンサンブルの際の重みを最適化するのにも使いました。
特徴量Xで学習したmodelsが予測した値pred1と、特徴量X2で学習したmodels2が予測した値pred2を線形結合してpredを求める時、

pred = α × (w × pred1 + (1 - w) × pred2) + β

w、α、βを最適化するためのソースコードは下記の通り。configの中身がalpha、beta、wに代わっている。

.py

def test_fold(config:Dict, models, X, X2, y, year_month_train, year_month_test) -> float:
    first_dates_month = pd.to_datetime(X[['year', 'month']].assign(day=1))
    
    train_index = first_dates_month.isin(year_month_train)
    test_index = first_dates_month.isin(year_month_test)
    
    X_test = X[test_index]
    y_test = y[test_index]
    X2_test = X2[test_index]
    
    alpha = config["alpha"]
    beta = config["beta"]
    w = config["w"]
    y_test_pred1 = predict(X_test, models)
    y_test_pred2 = predict(X2_test, models2)
    
    y_test_pred = (alpha * (w * y_test_pred1 + (1 - w) * y_test_pred2) + beta).clip(0)
    return mean_absolute_error(y_test, y_test_pred)


# loss function
def objective(trial):
    config = {
        'alpha': trial.suggest_float('alpha', 0.8, 1.2),
        'beta': trial.suggest_float('beta', -10, 10),
        'w' : trial.suggest_float('w', 0, 1)
    }

    cv = TimeSeriesSplit(n_splits=3)
    cv_mae = [None] * 3
    for i, (train_index, test_index) in enumerate(cv.split(timesteps)):
        cv_mae[i] = test_fold(config, models, X, X2, y, timesteps[train_index], timesteps[test_index])
        print(f"CV-{i}: mae =", cv_mae[i])

    # saving the individual fold holdout metrics 
    # uncomment this line if you don't want this
    ave_mae = np.mean(cv_mae)
    max_mae = np.max(cv_mae)
    
    trial.set_user_attr('max_mae', max_mae)
    trial.set_user_attr('ave_mae', ave_mae)
    
    return max_mae


sampler = optuna.samplers.TPESampler(
    n_startup_trials=10, seed=777
)

study = optuna.create_study(
    directions=['minimize'],
    sampler=sampler,
    study_name='model ensemble'
)

# maximum of 10000 trials or 9 hr wall clock time
study.optimize(objective, n_trials=10000, timeout= 9 * 3600) 
_ = joblib.dump(study, 'ensemble_enefit_hyperopt.pkl')