予測モデルを作って、説明変数から目的変数を予測することはできますが、反対に目的変数を最大化(最小化)する説明変数を探したいことがあります。サロゲート・モデル、配合調整、プロセス・インフォマティクスなどのパラメーター最適化をしたいケースです。
一方で通常最適化の技術では、定式関数であれば最適解を求めることができますが、関数の中身が不明な場合にはそれはできませんでした。
CPLEXのブラックボックス最適化は関数の中身が不明な場合でも最適解を探索することができるという機能になります。C言語ではバージョン20.1から、Pythonでは22.1からサポートされています。この記事ではPythonを使って説明します。
この機能をつかうことで予測モデルをブラックボックス関数として定義して目的変数を最大化(最小化)する説明変数を探すことができます。
- テスト環境
- cplex:22.1.1.0
- docplex:2.25.236
- python 3.10.9
- Windows 11 64 bit
データと問題
以下のような車のデータです。
mpg:燃費
engine:エンジン体積
horsepower:馬力
weight:重量
acceleration:加速性
このデータから馬力を最大化できるような、燃費、エンジン体積、重量、加速性を探していきます。
線形回帰モデルの作成
まずは「燃費、エンジン体積、重量、加速性」から「馬力」を予測する重回帰モデルを作ります。
from sklearn.linear_model import LinearRegression
lr = LinearRegression()
# 説明変数:燃費、エンジン体積、重量、加速性
X = df[['mpg', 'engine', 'weight', 'acceleration']] .to_numpy().tolist()
# 目的変数:馬力
Y = df['horsepower'].to_numpy().tolist()
# 線形回帰のモデルを作成
lr = LinearRegression()
lr.fit(X, Y)
print('coefficient = ', lr.coef_)
print('intercept = ', lr.intercept_)
結果は以下のようになりました。
coefficient = [-0.39658708 0.10224376 0.01800971 -4.72344381]
intercept = 113.6089240375573
得られた重回帰式は以下のようになります。
mpg * -0.39658708 + engine * 0.10224376 + weight * 0.01800971 + acceleration * -4.72344381 + 113.6089240375573
エンジン体積や重量は大きい方が、燃費や加速性は低い方が馬力が大きいということですね。
線形問題としてパラメーター(説明変数)を最適化する
線形回帰の場合は、式が明らかなので、ブラックボックス関数ではなくシンプルな線形最適化が可能です。
以下のコードになります。
# MPを利用
from docplex.mp.model import Model
m = Model(name='Linear')
# 決定変数=説明変数。上限下限は元データより定義
mpg = m.continuous_var(df['mpg'].min(), df['mpg'].max(), "mpg")
engine = m.continuous_var(df['engine'].min(), df['engine'].max(), "engine")
weight = m.continuous_var(df['weight'].min(), df['weight'].max(), "weight")
acceleration = m.continuous_var(df['acceleration'].min(), df['acceleration'].max(), "acceleration")
# 目的関数:目的変数の最大化
m.maximize(mpg * -0.39658708 + engine * 0.10224376 + weight * 0.01800971 + acceleration * -4.72344381 + 113.6089240375573)
# 求解:線形問題
s = m.solve()
m.print_solution()
簡単に主なところを解説します。
以下は、「燃費、エンジン体積、重量、加速性」という回帰モデルの説明変数を決定変数として定義しています。上限と下限は元データから取得しています。
mpg = m.continuous_var(df['mpg'].min(), df['mpg'].max(), "mpg")
engine = m.continuous_var(df['engine'].min(), df['engine'].max(), "engine")
weight = m.continuous_var(df['weight'].min(), df['weight'].max(), "weight")
acceleration = m.continuous_var(df['acceleration'].min(), df['acceleration'].max(), "acceleration")
以下は、先ほど得られた重回帰式をそのまま目的関数として定義しています。つまり、「燃費、エンジン体積、重量、加速性」を入力として回帰式で得られた「馬力」を最大化するという意味です。
m.maximize(mpg * -0.39658708 + engine * 0.10224376 + weight * 0.01800971 + acceleration * -4.72344381 + 113.6089240375573)
結果は以下のようになりました。
objective: 211.343
status: OPTIMAL_SOLUTION(2)
mpg=9.000
engine=455.000
weight=5140.000
acceleration=8.000
つまり
燃費=9.000
エンジン体積=455.000
重量=5140.000
加速性=8.000
の車だと、馬力=211.343が期待できるということになります。
実は単なる線形関数なので、負の係数を持つ説明変数は最低値を、正の係数を持つ説明変数は最大値を取っているだけです。
回帰モデルをブラックボックス関数としてパラメーター(説明変数)最適化する
次はついにブラックボックス関数を使ってみます。
回帰式をモデルのまま関数にして、これをCPLEXの目的関数として定義します。
まず回帰モデルをpythonの関数「predictHP」として定義します。
元々predictは説明変数の組合せを複数件入力して、複数件の目的変数が返ってくるので、配列の次元数を変換しています。
def predictHP(mpg, engine, weight, acceleration):
return lr.predict([[mpg, engine, weight, acceleration]])[0]
# テスト実行
predictHP(10, 10, 10, 10)
テスト結果は以下のようになりました。ちゃんとスコアリングができています。
63.611149954052415
次がメインのロジックです。
# CPを利用
from docplex.cp.model import *
import docplex.cp.solver.solver as solver
# ブラックボックス関数として線形回帰の関数を定義
bbf = CpoBlackboxFunction(predictHP)
mdl = CpoModel()
# 決定変数=説明変数。上限下限は元データより定義。CPは浮動小数点は扱えないので整数化
mpg = integer_var(df['mpg'].min(), df['mpg'].max(), "mpg")
engine = integer_var(int(df['engine'].min()), int(df['engine'].max()), "engine")
weight = integer_var(df['weight'].min(), df['weight'].max(), "weight")
acceleration = integer_var(int(df['acceleration'].min()), int(df['acceleration'].max()), "acceleration")
# 目的関数:目的変数の最大化。ブラックボックス関数を呼び出し
obj = bbf(mpg, engine, weight, acceleration)
mdl.maximize(obj)
# 求解:CP問題
msol = mdl.solve(Workers=1, TimeLimit=5)
msol.print_solution()
簡単に主なところを解説します。
まず、CPのエンジンを使っています。回帰式をそのまま埋め込んだ際には線形問題だったのでMPのエンジンを使っていました。ブラックボックス関数はCPエンジンの機能になります。
from docplex.cp.model import *
import docplex.cp.solver.solver as solver
以下では、先に定義したpython関数「predictHP」をCPのブラックボックス関数として定義しています。これで「predictHP」は「bbf」という関数として最適化式の中で利用可能になります。
bbf = CpoBlackboxFunction(predictHP)
以下では、「燃費、エンジン体積、重量、加速性」という回帰モデルの説明変数を決定変数として定義しています。上限と下限は元データから取得しています。
線形問題とよく似ています。しかし、線形問題ではcontinuous_varを使いましたが、CP問題ではcontinuous_varは使えないため、integer_varを使っています。
engineとaccelerationは小数点を含んでいるため、ここでは単純に整数に切り捨てています。(固定小数点であれば、あらかじめ10^nをかけておくことで、小数点以下の情報も活かすことはできます)
# 決定変数=説明変数。上限下限は元データより定義。CPは浮動小数点は扱えないので整数化
mpg = integer_var(df['mpg'].min(), df['mpg'].max(), "mpg")
engine = integer_var(int(df['engine'].min()), int(df['engine'].max()), "engine")
weight = integer_var(df['weight'].min(), df['weight'].max(), "weight")
acceleration = integer_var(int(df['acceleration'].min()), int(df['acceleration'].max()), "acceleration")
以下では、先ほど定義したブラックボックス関数をそのまま目的関数として定義しています。つまり、「燃費、エンジン体積、重量、加速性」を入力として、ブラックボックス関数である回帰モデルから得られた「馬力」を最大化するという意味です。
obj = bbf(mpg, engine, weight, acceleration)
mdl.maximize(obj)
以下は、求解をしています。線形問題と異なり、最適解を求めることには時間がかかりますので、時間制限TimeLimit=5をかけています。
msol = mdl.solve(Workers=1, TimeLimit=5)
msol.print_solution()
結果は以下のようになりました。
! --------------------------------------------------- CP Optimizer 22.1.1.0 --
! Maximization problem - 4 variables, 0 constraints, 1 blackbox
! TimeLimit = 5
! Workers = 1
! Initial process time : 0.01s (0.01s extraction + 0.00s propagation)
! . Log search space : 29.7 (before), 29.7 (after)
! . Memory usage : 275.4 kB (before), 275.4 kB (after)
! Using sequential search.
! ----------------------------------------------------------------------------
! Best Branches Non-fixed Branch decision
0 4 -
+ New bound is 1.797693e+308
* 108.2543 816 0.13s (gap > 10000%)
* 145.1055 824 0.13s (gap > 10000%)
145.1055 1000 1 23 != acceleration
* 189.7611 1053 0.15s (gap > 10000%)
* 191.9457 1290 0.16s (gap > 10000%)
* 198.4070 1700 0.23s (gap > 10000%)
198.4070 2000 1 1859 != weight
198.4070 3000 1 F 4547 = weight
* 198.6521 3329 0.49s (gap > 10000%)
198.6521 4000 1 F 3390 = weight
* 199.0487 4383 0.62s (gap > 10000%)
199.0487 5000 1 F 4017 = weight
* 199.4453 5501 0.75s (gap > 10000%)
199.4453 6000 1 40 != mpg
* 199.8418 6653 0.88s (gap > 10000%)
199.8418 7000 1 F 2335 = weight
* 200.2384 7839 1.01s (gap > 10000%)
200.2384 8000 1 239 != engine
200.2384 9000 1 F 21 = acceleration
! Time = 1.12s, Average fail depth = 8, Memory usage = 1.3 MB
! Current bound is 1.797693e+308 (gap > 10000%)
! Best Branches Non-fixed Branch decision
200.2384 10000 1 1903 != weight
200.2384 11000 1 F 2280 = weight
200.2384 12000 1 40 = mpg
200.2384 13000 1 F 2289 = weight
200.2384 14000 1 135 != engine
200.2384 15000 1 95 != engine
200.2384 16000 1 F 2286 = weight
200.2384 17000 1 F 2795 = weight
200.2384 18000 1 2149 != weight
200.2384 19000 1 2674 != weight
200.2384 20000 1 F 10 = acceleration
200.2384 21000 1 F 15 = mpg
200.2384 22000 1 4610 != weight
200.2384 23000 1 2645 != weight
* 211.3429 23778 3.18s (gap > 10000%)
211.3429 24000 1 27 != mpg
211.3429 25000 1 1617 != weight
211.3429 26000 1 8 = acceleration
211.3429 27000 1 3253 != weight
211.3429 28000 1 4489 != weight
! Time = 3.69s, Average fail depth = 8, Memory usage = 1.7 MB
! Current bound is 1.797693e+308 (gap > 10000%)
! Best Branches Non-fixed Branch decision
211.3429 29000 1 1690 != weight
211.3429 30000 1 F 1731 = weight
211.3429 31000 1 13 != acceleration
211.3429 32000 1 F 3215 = weight
211.3429 33000 1 83 != engine
211.3429 34000 1 F 5092 = weight
211.3429 35000 1 25 != mpg
211.3429 36000 1 F 2428 = weight
211.3429 37000 1 F 1710 = weight
211.3429 38000 1 F 9 != acceleration
211.3429 38351 1 F 3442 = weight
+ New bound is 1.572981e+308 (gap > 10000%)
! ----------------------------------------------------------------------------
! Search terminated by limit, 11 solutions found.
! Best objective : 211.3429 (gap > 10000%)
! Best bound : 1.572981e+308
! ----------------------------------------------------------------------------
! Number of branches : 38351
! Number of fails : 18816
! Number of b-box calls : 15954
! Total memory usage : 1.7 MB (1.7 MB CP Optimizer + 0.0 MB Concert)
! Time spent in solve : 5.01s (5.00s engine + 0.01s extraction)
! Search speed (br. / s) : 7671.7
! ----------------------------------------------------------------------------
-------------------------------------------------------------------------------
Model constraints: 0, variables: integer: 4, interval: 0, sequence: 0
Solve status: Feasible
Search status: SearchStopped, stop cause: SearchStoppedByLimit
Solve time: 5.0 sec
-------------------------------------------------------------------------------
Objective values: (211.3429232616431,), bounds: (1.5729814930045262e+308,), gaps: (7.44279e+305,)
Variables:
acceleration = 8
engine = 455
mpg = 9
weight = 5140
結果は
燃費=9
エンジン体積=455
重量=5140
加速性=8
の車だと、馬力=211.3429232616431が期待できるということになります。
うまいこと回帰式での結果と同じになりました。ランダムの要素があるので毎回同じになるとは限りませんが、求解時間を長くすれば安定はしやすくなります。
ログを見ると以下のようにありますので、11個の解をみつけて5秒で打ち切っています。
Search terminated by limit, 11 solutions found.
----中略-----
Search status: SearchStopped, stop cause: SearchStoppedByLimit
Solve time: 5.0 sec
もう少しログを見てみると、いろいろな説明変数の値を試して、少しずつ目的関数である「馬力」の値が改善されているのがわかります。
ブラックボックスなので、トライアンドエラーをしていることがわかります。
198.4070 2000 1 1859 != weight
198.4070 3000 1 F 4547 = weight
* 198.6521 3329 0.49s (gap > 10000%)
198.6521 4000 1 F 3390 = weight
* 199.0487 4383 0.62s (gap > 10000%)
199.0487 5000 1 F 4017 = weight
* 199.4453 5501 0.75s (gap > 10000%)
199.4453 6000 1 40 != mpg
* 199.8418 6653 0.88s (gap > 10000%)
199.8418 7000 1 F 2335 = weight
回帰式をブラックボックス関数として記述して実行
今行ったことをさらにわかりやすくするために、回帰式をブラックボックス関数として記述してみました。
# CPを利用
from docplex.cp.model import *
import docplex.cp.solver.solver as solver
# ブラックボックス関数として線形回帰の式を定義
bbf = CpoBlackboxFunction(lambda mpg, engine, weight, acceleration:
mpg * -0.39658708 + engine * 0.10224376 + weight * 0.01800971 + acceleration * -4.72344381 + 113.6089240375573)
mdl = CpoModel()
# 決定変数=説明変数。上限下限は元データより定義。CPは浮動小数点は扱えないので整数化
mpg = integer_var(9, 46, "mpg")
engine = integer_var(68, 455, "engine")
weight = integer_var(1613, 5140, "weight")
acceleration = integer_var(8, 25, "acceleration")
# 目的関数:目的変数の最大化。ブラックボックス関数を呼び出し
obj = bbf(mpg, engine, weight, acceleration)
mdl.maximize(obj)
# 求解:CP問題
msol = mdl.solve(Workers=1, TimeLimit=5)
msol.print_solution()
lamdaで無名関数で回帰式をそのままCpoBlackboxFunctionの中に定義しています。
bbf = CpoBlackboxFunction(lambda mpg, engine, weight, acceleration:
mpg * -0.39658708 + engine * 0.10224376 + weight * 0.01800971 + acceleration * -4.72344381 + 113.6089240375573)
上限と下限もベタ書きしてみました。この方が意味は理解しやすいと思います。
mpg = integer_var(9, 46, "mpg")
engine = integer_var(68, 455, "engine")
weight = integer_var(1613, 5140, "weight")
acceleration = integer_var(8, 25, "acceleration")
実行結果は変わらないので割愛します。ブラックボックス関数として定義した式の中身をCPエンジンは解釈せずにブラックボックスのまま実行しています。動きがよりイメージしやすいかと思います。
制約問題として式化できないランダムフォレストモデルをブラックボックス関数としてパラメーター(説明変数)最適化する
さてここまでは理解しやすいように、回帰モデル、回帰式をつかってブラックボックス関数を説明してきました。しかし、もちろん回帰モデルであれば、線形問題として解くのが正確かつ効率的で、ブラックボックス関数を使う必要は全くありません(というかそもそも他に制約がなければ、各説明変数の最小値や最大値を選べばいいだけです)。
しかし、機械学習では、回帰モデルのように簡単な式で書けないモデルがほとんどです。このような場合にブラックボックス関数は真価を発揮します。
ここでは式で書くことが困難なランダムフォレストのモデルをブラックボックス関数として利用していきます。
まず、ランダムフォレストのモデルをつくり、先ほどの回帰モデル同様に「predictHP」というpythonの関数を作ります。
from sklearn.ensemble import RandomForestRegressor
# 説明変数:燃費、エンジン体積、重量、加速性
X = df[['mpg', 'engine', 'weight', 'acceleration']] .to_numpy().tolist()
# 目的変数:馬力
Y = df['horsepower'].to_numpy().tolist()
# ランダムフォレストのモデルを作成
forest = RandomForestRegressor()
forest.fit(X, Y)
# ランダムフォレストのモデルを関数化
def predictHP(mpg, engine, weight, acceleration):
return forest.predict([[mpg, engine, weight, acceleration]])[0]
# テスト実行
predictHP(10, 10, 10, 10)
そして、最適化していきます。
実はコードの内容は回帰モデルと全く変わりません。「bbf = CpoBlackboxFunction(predictHP)」でブラックボックス関数としてランダムフォレストのモデルを関数化した「predictHP」を指定しているだけです。
# CPを利用
from docplex.cp.model import *
import docplex.cp.solver.solver as solver
# ブラックボックス関数としてランダムフォレストの関数を定義
bbf = CpoBlackboxFunction(predictHP)
mdl = CpoModel()
# 決定変数=説明変数。上限下限は元データより定義。CPは浮動小数点は扱えないので整数化
mpg = integer_var(df['mpg'].min(), df['mpg'].max(), "mpg")
engine = integer_var(int(df['engine'].min()), int(df['engine'].max()), "engine")
weight = integer_var(df['weight'].min(), df['weight'].max(), "weight")
acceleration = integer_var(int(df['acceleration'].min()), int(df['acceleration'].max()), "acceleration")
# 目的関数:目的変数の最大化。ブラックボックス関数を呼び出し
obj = bbf(mpg, engine, weight, acceleration)
mdl.maximize(obj)
# 求解:CP問題
msol = mdl.solve(Workers=1, TimeLimit=5)
msol.print_solution()
結果は以下のようになりました。
! --------------------------------------------------- CP Optimizer 22.1.1.0 --
! Maximization problem - 4 variables, 0 constraints, 1 blackbox
! TimeLimit = 5
! Workers = 1
! Initial process time : 0.01s (0.01s extraction + 0.00s propagation)
! . Log search space : 29.7 (before), 29.7 (after)
! . Memory usage : 275.4 kB (before), 275.4 kB (after)
! Using sequential search.
! ----------------------------------------------------------------------------
! Best Branches Non-fixed Branch decision
0 4 -
+ New bound is 1.797693e+308
* 74.39000 816 2.29s (gap > 10000%)
* 162.2500 820 2.29s (gap > 10000%)
* 182.8100 824 2.30s (gap > 10000%)
* 208.7800 836 2.31s (gap > 10000%)
208.7800 1000 1 F 1693 = weight
* 209.2800 1272 2.94s (gap > 10000%)
! ----------------------------------------------------------------------------
! Search terminated by limit, 5 solutions found.
! Best objective : 209.2800 (gap > 10000%)
! Best bound : 1.797693e+308
! ----------------------------------------------------------------------------
! Number of branches : 1992
! Number of fails : 966
! Number of b-box calls : 868
! Total memory usage : 1.1 MB (1.0 MB CP Optimizer + 0.0 MB Concert)
! Time spent in solve : 5.01s (5.00s engine + 0.01s extraction)
! Search speed (br. / s) : 398.6
! ----------------------------------------------------------------------------
-------------------------------------------------------------------------------
Model constraints: 0, variables: integer: 4, interval: 0, sequence: 0
Solve status: Feasible
Search status: SearchStopped, stop cause: SearchStoppedByLimit
Solve time: 5.0 sec
-------------------------------------------------------------------------------
Objective values: (209.28,), bounds: (1.7976931348623157e+308,), gaps: (8.58989e+305,)
Variables:
acceleration = 17
engine = 455
mpg = 46
weight = 5140
結果は
燃費=46
エンジン体積=455
重量=5140
加速性=17
の車だと、馬力=209.28が期待できるということになります。
ロジックのわからないランダムフォレスト・モデルの説明変数を最適化することができました!
なお、回帰モデルと比較してどちらの「馬力」がいい結果なのかは、そもそもモデルが違うので比較ができません。それよりもまずモデル自体の性能を別に評価して選択すべきです。
しかしながら、探索できたブランチ数が全く異っていることには注意が必要です。
- 回帰モデル:38,351
- ランダムフォレスト・モデル:1,992
これは一回あたりのスコアリングにかかる負荷が回帰モデルよりランダムフォレスト・モデルの方が重いからです。ブラックボックス関数では内部で何度も関数を呼び出し、その度にスコアリングが行われています。その差が表れています。ランダムフォレストの場合には、TimeLimitをもっと長くした方がよい結果が得られるかもしれません。
サンプル
サンプルコード
サンプルデータ
参考
Module docplex.cp.blackbox
CP オプティマイザーの新しいブラック・ボックス最適化機能 - IBM Documentation
CP Optimizer black-box expressions in Python - IBM Documentation
SPSS Modelerでパラメーター最適化ーサロゲート・モデル、配合調整、プロセス・インフォマティクスなどー #SPSS - Qiita