1
Help us understand the problem. What are the problem?

posted at

Scipyで多変数関数の最小値を求める(逐次二次計画法の利用)

はじめに

pythonのscipyモジュールを使用して多変数関数の最小値を求める方法について紹介します。
scipyの公式レファレンスには多くの方法が記載されていますが、今回は特に逐次二次計画法(SLSQP)について取り上げます。

実行環境

OS : Windows10
環境 : Anaconda(エディタ:Jupyter Notebook)
言語 : Python3
使用ライブラリ : scpiy,numpy

逐次二次計画法(SLSQP)

逐次二次計画法 (Sequenitial Least Quadratic Programming) とは、制約付き非線形計画問題を解く手法の1つです。

最適化したい目的関数に対して、非線形な等式制約条件式と不等式制約条件式が複数ある場合に適用することができます。

詳しい理論については以下のリンクから参照してみてください。

例題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という変数に計算結果が格納されます。

例題1の実行結果

     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}")

例題1のコードまとめ

ここをクリック
逐次二次計画法の適用
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)

例題2【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)

例題2の実行結果

     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])

こちらも求めることができました。

解探索を始める初期点によって結果は変わりますので、色々試してみてください。

参考文献

Register as a new user and use Qiita more conveniently

  1. You can follow users and tags
  2. you can stock useful information
  3. You can make editorial suggestions for articles
What you can do with signing up
1
Help us understand the problem. What are the problem?