- 製造業出身のデータサイエンティストがお送りする記事
- 今回はphysboを使って多目的最適化問題をベイズ最適化で解いてみました。
##はじめに
過去に多目的最適化問題については、記事を何個か書いておりますので参考にして頂けますと幸いです。
- 多目的最適化問題のNSGA-Ⅱを勉強したの整理しました
- 多目的最適化問題のNSGA-Ⅱを実装してみた
- 多目的最適化問題のNSGA-Ⅲを勉強したの整理しました
- 多目的最適化問題のNSGA-Ⅲを実装してみた
- pymooを使って多目的最適化のNSGA-Ⅱを実装してみた
##使用するライブラリー（physbo）
今回は最適化ライブラリphysboを使って実装したいと思います。
##physboの実装
今回は実装サンプルを試してみました。
# 必要なライブラリーのインストール
import numpy as np
import matplotlib.pyplot as plt
import physbo
import itertools
%matplotlib inline
次にテスト関数を作成します。
def vlmop2_minus(x):
n = x.shape[1]
y1 = 1 - np.exp(-1 * np.sum((x - 1/np.sqrt(n)) ** 2, axis = 1))
y2 = 1 - np.exp(-1 * np.sum((x + 1/np.sqrt(n)) ** 2, axis = 1))
return np.c_[-y1, -y2]
探索範囲を指定します。
a = np.linspace(-2, 2, 101)
test_X = np.array(list(itertools.product(a, a)))
シミュレータを定義します。
class simulator(object):
def __init__(self, X):
self.t = vlmop2_minus(X)
def __call__( self, action):
return self.t[action]
simu = simulator(test_X)
今回の設定した関数は、1つは右上にピークがあり、1つは左下にピークがある関数です。
可視化すると以下のような感じです。
# plot objective 1
plt.figure()
plt.imshow(simu.t[:,0].reshape((101,101)),
vmin=-1.0,
vmax=0.0,
origin="lower",
extent=[-2.0, 2.0, -2.0, 2.0]
)
plt.title("objective 1")
plt.colorbar()
plt.plot([1.0/np.sqrt(2.0)],
[1.0/np.sqrt(2.0)],
'*'
)
plt.show()
# plot objective 2
plt.figure()
plt.imshow(simu.t[:,1].reshape((101,101)),
vmin=-1.0,
vmax=0.0,
origin="lower",
extent=[-2.0, 2.0, -2.0, 2.0]
)
plt.title("objective 2")
plt.colorbar()
plt.plot([-1.0/np.sqrt(2.0)],
[-1.0/np.sqrt(2.0)],
'*'
)
plt.show()
次に最適化を実行していきます。
- policyのセット
- 多目的最適化用のphysbo.search.discrete_multi.policyを利用
- num_objectivesに目的関数の数を指定
policy = physbo.search.discrete_multi.policy(test_X=test_X,
num_objectives=2
)
policy.set_seed(0)
ベイズ最適化は、下記3パターンの方法があるそうです。
- HVPI (Hypervolume-based Probability of Improvement)
- EHVI (Expected Hyper-Volume Improvement)
- TS (Thompson Sampling)
今回はHVPIを使用してみたいと思います。
HVPI (Hypervolume-based Probability of Improvement)は、多次元の目的関数空間における非劣解領域 (non-dominated region) の改善確率をscoreとして算出する手法だそうです。
policy = physbo.search.discrete_multi.policy(test_X=test_X,
num_objectives=2
)
policy.set_seed(0)
policy.random_search(max_num_probes=10,
simulator=simu
)
res_HVPI = policy.bayes_search(max_num_probes=40,
simulator=simu,
score='HVPI',
interval=10
)
パレート解の取得は以下のコードでできます。
front, front_num = res_HVPI.export_pareto_front()
front, front_num
結果はこんな感じで返ってきます。
各目的変数の値と対象のインデックス番号が得られます。
(array([[-9.82978321e-01, -3.32414921e-04],
[-9.71056940e-01, -3.45637262e-02],
[-9.50786244e-01, -6.76228653e-02],
[-9.34458810e-01, -1.15668709e-01],
[-9.13824628e-01, -1.71906452e-01],
[-8.88492700e-01, -2.36876025e-01],
[-8.68576210e-01, -2.82707709e-01],
[-8.57549341e-01, -3.05692550e-01],
[-8.33710080e-01, -3.53626320e-01],
[-8.20907701e-01, -3.78330099e-01],
[-7.92933832e-01, -4.26778524e-01],
[-7.77704642e-01, -4.50447880e-01],
[-7.62117827e-01, -4.74823143e-01],
[-7.45438101e-01, -4.98117246e-01],
[-7.28459156e-01, -5.21910482e-01],
[-7.14031659e-01, -5.50368010e-01],
[-6.92015440e-01, -5.67552511e-01],
[-6.72522879e-01, -5.89370083e-01],
[-6.52909089e-01, -6.11332656e-01],
[-6.32120559e-01, -6.32120559e-01],
[-6.11332656e-01, -6.52909089e-01],
[-5.89370083e-01, -6.72522879e-01],
[-5.70311337e-01, -6.93980247e-01],
[-5.44575729e-01, -7.10347737e-01],
[-5.21910482e-01, -7.28459156e-01],
[-4.98117246e-01, -7.45438101e-01],
[-4.74823143e-01, -7.62117827e-01],
[-4.26778524e-01, -7.92933832e-01],
[-4.02089723e-01, -8.07119688e-01],
[-3.90152305e-01, -8.24313473e-01],
[-3.53626320e-01, -8.33710080e-01],
[-3.05692550e-01, -8.57549341e-01],
[-2.58961961e-01, -8.78749507e-01],
[-2.16642598e-01, -8.97780597e-01],
[-1.71906452e-01, -9.13824628e-01],
[-1.15668709e-01, -9.34458810e-01],
[-6.76228653e-02, -9.50786244e-01],
[-3.51124618e-02, -9.67608404e-01],
[-1.46847591e-03, -9.78680443e-01]]),
array([21, 15, 35, 26, 32, 18, 43, 27, 40, 16, 31, 49, 23, 47, 33, 12, 34,
45, 22, 44, 28, 41, 13, 42, 29, 48, 19, 30, 46, 10, 37, 24, 38, 17,
36, 25, 39, 11, 20]))
パレート解の可視化は以下のような感じでできます。
def plot_pareto_front(res):
front, front_num = res.export_pareto_front()
dominated = [i for i in range(res.num_runs) if i not in front_num]
points = res.fx[dominated, :]
plt.figure(figsize=(7, 7))
plt.scatter(res.fx[dominated,0], res.fx[dominated,1], c = "blue")
plt.scatter(front[:, 0], front[:, 1], c = "red")
plt.title('Pareto front')
plt.xlabel('Objective 1')
plt.ylabel('Objective 2')
plt.xlim([-1.0, 0.0])
plt.ylim([-1.0, 0.0])
plot_pareto_front(res_HVPI)
ちなみに全探索もできます。
policy = physbo.search.discrete_multi.policy(test_X=test_X,
num_objectives=2
)
policy.set_seed(0)
N = test_X.shape[0]
res_all = policy.random_search(max_num_probes=N,
simulator=simu
)
plot_pareto_front(res_all)
##さいごに
最後まで読んで頂き、ありがとうございました。
今回は、physboを使って多目的最適化問題をベイズ最適化を使って解きました。
今回は単純な問題だったので時間はかかりませんでしたが、複雑な問題となった場合、実行時間がどこまでかかるのかは気になりました。
時間がある時に試してみようと思います。
訂正要望がありましたら、ご連絡頂けますと幸いです。