はじめに
pythonのscipyモジュールを使用して多変数関数の最小値を求める方法について紹介します。
scipyの公式レファレンスには多くの方法が記載されていますが、今回は特に逐次二次計画法(SLSQP)について取り上げます。
実行環境
OS : Windows10
環境 : Anaconda(エディタ:Jupyter Notebook)
言語 : Python3
使用ライブラリ : scpiy,numpy
逐次二次計画法（SLSQP）
逐次二次計画法 (Sequenitial Least Quadratic Programming) とは、制約付き非線形計画問題を解く手法の１つです。
最適化したい目的関数に対して、非線形な等式制約条件式と不等式制約条件式が複数ある場合に適用することができます。
詳しい理論については以下のリンクから参照してみてください。
例題1【2変数関数の最小値を求める】
まず、以下の等式制約条件式と不等式制約条件式における2変数関数の最小値とその時の解、${x_1,x_2}$を求めてみます。
\text{minimize} \quad x_1^2+x_1\times x_2\\\
\\\
\text{subject to}
\begin{cases}
x_1^3+x_1\times x_2 =100\\\
\\\
x_1^2+4x_2 \geq 50 \\\
\\\
-100\leq x_1,x_2\leq 100
\end{cases}
ライブラリのインポート
from scipy.optimize import minimize
import numpy
まず、最適化に必要なライブラリをインポートしておきます。
scipyがインストールされていない場合は、コマンドプロンプトかAnaconda promptからpip install scipyコマンドを打つことでインストールすることができます。
目的関数と制約条件式の設定
\text{minimize} \quad x_1^2+x_1\times x_2\\\
\\\
\text{subject to}
\begin{cases}
x_1^3+x_1\times x_2 =100\\\
\\\
x_1^2+4x_2 \geq 50 \\\
\\\
-100\leq x_1,x_2\leq 100
\end{cases}
# 目的関数
def objective_fnc(x):
x1 = x[0]
x2 = x[1]
return x1**2 + x1*x2
# 等式制約条件
def equality_constraint(x):
x1 = x[0]
x2 = x[1]
return x1**3 + x1*x2 - 100
# 不等式制約条件
def inequality_constraint(x):
x1 = x[0]
x2 = x[1]
return x1**2 + 4*x2 - 50
bounds_x1 = (-100,100)
bounds_x2 = (-100,100)
bound = [bounds_x1,bounds_x2]
等式制約条件式および不等式制約条件式については全て左辺に移項した式で定義する必要があります。
逐次二次計画法の適用
constraint1 = {"type":"eq","fun":equality_constraint}
constraint2 = {"type":"ineq","fun":inequality_constraint}
constraint = [constraint1,constraint2]
x0=[5,5]
result=minimize(objective_fnc, x0, method="SLSQP", bounds=bound, constraints=constraint)
print(result)
解探索を始める初期点を$(x_1,x_2)=(5,5)$に設定して実行してみます。
また、今回はresultという変数に計算結果が格納されます。
例題１の実行結果
fun: 50.34449775048212
jac: array([16.49958134, 4.0409255 ])
message: 'Optimization terminated successfully'
nfev: 15
nit: 5
njev: 5
status: 0
success: True
x: array([4.04092525, 8.41773078])
最小値および最適解を求めることができました。
実行結果の主な内訳としては、以下の表のようになっています。
| 関数名 | 意味 |
|---|---|
| fun | 最適化された目的関数の値 |
| jac | ヤコビアンの値 |
| message | 探索終了原因の記述 |
| nfev, njev | 目的関数とそのヤコビアンの評価数 |
| nit | 最適解に至るまでの反復回数 |
| status, success | 終了状況 |
| x | 最適解$(x_1,x_2)$の値 |
解の取り出し
また、得られた解についてはこのように変数に代入することで取り出せます。
x1,x2 = result.x[0],result.x[1]
print(f"x1 = {x1}",f"x2 = {x2}")
例題１のコードまとめ
ここをクリック
from scipy.optimize import minimize
import numpy
# 最小化したい目的関数
def objective_fnc(x):
x1 = x[0]
x2 = x[1]
return x1**2 + x1*x2
# 等式制約条件
def equality_constraint(x):
x1 = x[0]
x2 = x[1]
return x1**3 + x1*x2 - 100
# 不等式制約条件
def inequality_constraint(x):
x1 = x[0]
x2 = x[1]
return x1**2 + x1*x2 - 50
bounds_x1 = (-100,100)
bounds_x2 = (-100,100)
bounds = [bounds_x1,bounds_x2]
constraint1 = {"type":"eq","fun":equality_constraint}
constraint2 = {"type":"ineq","fun":inequality_constraint}
constraint = [constraint1,constraint2]
x0=[5,5]
result = minimize(objective_fnc, x0, method = "SLSQP", bounds = bounds, constraints = constraint)
print(result)
例題２【4変数関数の最小値を求める】
もう1問解いてみます。
解探索の初期点は、
$(x_1,x_2,x_3,x_4)=(1,5,5,1)$とします。
\text{minimize} \quad x_1x_4(x_1+x_2+x_3)+x_3\\\
\\
\text{subject to}
\begin{cases}
x_1^2+x_2^2+x_3^2+x_4^2 =40\\\
\\
x_1x_2x_3x_4 \geq 25 \\\
\\
1\leq x_1,x_2,x_3,x_4\leq 5
\end{cases}
逐次二次計画法の実装
from scipy.optimize import minimize
import numpy
# 目的関数
def objective_fnc(x):
x1=x[0]
x2=x[1]
x3=x[2]
x4=x[3]
return x1*x4*(x1 + x2 + x3) + x3
# 等式制約条件
def equality_constraint(x):
x1=x[0]
x2=x[1]
x3=x[2]
x4=x[3]
return x1**2 + x2**2 + x3**2 + x4**2 -40
# 不等式制約条件
def inequality_constraint(x):
x1=x[0]
x2=x[1]
x3=x[2]
x4=x[3]
return x1*x2*x3*x4 - 25
bounds_x1=(1,5)
bounds_x2=(1,5)
bounds_x3=(1,5)
bounds_x4=(1,5)
bounds=[bounds_x1,bounds_x2,bounds_x3,bounds_x4]
constraint1={"type":"eq","fun":equality_constraint}
constraint2={"type":"ineq","fun":inequality_constraint}
constraint=[constraint1,constraint2]
x0=[1,5,5,1]
result=minimize(objective_fnc,x0,method="SLSQP",bounds=bounds,constraints=constraint)
print(result)
例題２の実行結果
fun: 17.01401724556073
jac: array([14.57227039, 1.37940764, 2.37940764, 9.56415081])
message: 'Optimization terminated successfully'
nfev: 25
nit: 5
njev: 5
status: 0
success: True
x: array([1. , 4.74299607, 3.82115466, 1.37940764])
こちらも求めることができました。
解探索を始める初期点によって結果は変わりますので、色々試してみてください。