予測モデルを作って、説明変数から目的変数を予測することはできますが、反対に目的変数を最大化（最小化）する説明変数を探したいことがあります。サロゲート・モデル、配合調整、プロセス・インフォマティクスなどのパラメーター最適化をしたいケースです。

一方で通常最適化の技術では、定式関数であれば最適解を求めることができますが、関数の中身が不明な場合にはそれはできませんでした。

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エンジンの機能になります。

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をかけています。

ブラックボックス関数として線形回帰の関数を定義# 求解：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
 *      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)

上限と下限もベタ書きしてみました。この方が意味は理解しやすいと思います。

決定変数＝説明変数。上限下限は元データより定義。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")

実行結果は変わらないので割愛します。ブラックボックス関数として定義した式の中身を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