Python 説明変数間が従属関係にある制約条件下での多目的最適化「PyCaret×Optuna」

　本記事では、説明変数間が従属関係にある制約条件下での多目的最適化を実施して、下図のようなパレート解を取得する雛形コードを載せました。

f:id:HK29:20210530192629p:plain

下図はその実施例です。説明変数3つ「PTRATIO」「INDUS」「DIS」の和が固定値25という制約下で乱数により、多目的最適化を実施した結果です。

f:id:HK29:20210530190704p:plain

■本プログラム

#!/usr/bin/env python
# coding: utf-8

# In[1]:


# PyCaret用の回帰モデルを読み込む
from pycaret.regression import *

load_model_y1 = load_model('best_model_PRICE')
load_model_y2 = load_model('best_model_RM')
print(load_model_y1)
print(load_model_y2)


# In[2]:


import os
import optuna
import numpy as np
import matplotlib.pyplot as plt

outfile_name = 'history_data'

if os.path.isfile(outfile_name + '.csv'):
    os.remove(outfile_name + '.csv')

# 目的関数の定義。複数の目的変数を戻り値とする
def objective(trial):
    # 指定範囲の乱数の生成
    lstat = trial.suggest_float('LSTAT', 1.7, 38) # 変数lstatを上下限1.7~38の範囲で浮動小数点
    ptratio = trial.suggest_float('PTRATIO', 12.6, 22)
    indus = trial.suggest_float('INDUS', 0.4, 24)
    dis = trial.suggest_float('DIS', 1.1, 12)
    age = trial.suggest_float('AGE', 3, 100)
    chas = trial.suggest_categorical('CHAS', ['0', '1']) # カテゴリ変数の場合
    
    ##### 説明変数間の制約条件の設定
    # 各パラメータをひとつのパラメータセットに合成する
    np_list = [ptratio, indus, dis]
    param_arr = np.vstack(np_list).T
    # パラメータA,B,Cの合計がSになるように変換する。
    S = 25
    arr = param_arr * S / param_arr.sum(axis=1)[:, np.newaxis]
    # アンパック代入
    ptratio, indus, dis = arr[0]
    #####
    
    # 予測するデータセットの作成（データフレーム型）
    df = pd.DataFrame([(lstat, ptratio, indus, dis, age, chas)],
                     columns= ['LSTAT', 'PTRATIO', 'INDUS', 'DIS', 'AGE', 'CHAS'])
    #print(df)
    
    # 予測する
    predict_y1 = predict_model(load_model_y1, data=df).values.tolist()[0][-1]
    predict_y2 = predict_model(load_model_y2, data=df).values.tolist()[0][-1]
    
    # pandasデータフレームに結合する
    df['predict_y1'] = predict_y1
    df['predict_y2'] = predict_y2    
    
    # 履歴をcsvファイルに出力する
    df.to_csv(outfile_name + '.csv', mode='a', index = False, header = False)
    
    return predict_y1, predict_y2
    
# 最適化の条件設定
study = optuna.multi_objective.create_study(
    directions=["minimize", "maximize"], # "minimize" "maximize"
    sampler=optuna.multi_objective.samplers.NSGAIIMultiObjectiveSampler(seed = 1)
    #sampler=optuna.multi_objective.samplers.RandomMultiObjectiveSampler(seed = 1)
)
# 最適化の実行
study.optimize(objective, n_trials=200)

# 作成した履歴csvファイルにヘッダーを追記する
DF = pd.read_csv(outfile_name + '.csv')
DF.columns = ['LSTAT', 'PTRATIO', 'INDUS', 'DIS', 'AGE', 'CHAS',
              'predict_y1', 'predict_y2']
DF.index = np.arange(1, len(DF) + 1) #インデックス番号を1から振り直す
DF.to_csv(outfile_name + '.csv', mode='w', header=True)

# パレート解の取得。get_pareto_front_trials()メソッドを使用
trials = {str(trial.values): trial for trial in study.get_pareto_front_trials()}
trials = list(trials.values())
trials.sort(key=lambda t: t.values)
# パレート解であった水準Noの抽出
pareto_no_list = []
for i, trial in enumerate(trials, start=1):
    print('pareto -> ', trial.number, trial.values[0], trial.values[1])
    pareto_no_list.append(trial.number)
# パレート解のみを抽出（indexを指定して抽出）
DF_pareto = DF[DF.index.isin(pareto_no_list)]
DF_pareto.to_csv('pareto_data.csv', index = True)
# パレート解以外を抽出
#DF_not_pareto = DF.drop(index = pareto_no_list)
#DF_not_pareto.to_csv('not_pareto_data.csv', index = True)

# パレート解を図示
plt.rcParams["font.size"] = 22
plt.figure(figsize=(10, 6))
plt.title("multiobjective optimization")
plt.xlabel("PRICE (predict_y1)")
plt.ylabel("RM (predict_y2)")
plt.grid()
plt.scatter(DF['predict_y1'], DF['predict_y2'], s=100, c='blue', label='all trials')
plt.scatter(DF_pareto['predict_y1'], DF_pareto['predict_y2'], s=100, c='red', label='pareto front')
# オプション
ax = plt.gca() # get current axes 現在の軸設定データを取得する
x_min, x_max = ax.get_xlim()
ax.set_xlim(0, x_max)
y_min, y_max = ax.get_ylim()
ax.set_ylim(0, y_max)
ax.legend(bbox_to_anchor=(1, 0.95)) # 凡例の表示と位置の指定
plt.tight_layout() # 全体を整える
#plt.savefig("pareto_graph.png")
#plt.close()


# In[ ]: