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

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


* [scipyの公式レファレンスはこちら](https://docs.scipy.org/doc/scipy-0.18.1/reference/optimize.html)


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


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

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

詳しい理論については以下のリンクから参照してみてください。
* [株式会社NTTデータ数理システム](http://www.msi.co.jp/nuopt/docs/v18/manual/html/0B-03-03.html)
* [制約を伴う最適化問題の逐次二次計画接近法による一解法](https://www.jstage.jst.go.jp/article/ieejeiss1972/102/9/102_9_207/_article/-char/ja/)


# 例題1【2変数関数の最小値を求める】
まず、以下の等式制約条件式と不等式制約条件式における2変数関数の最小値とその時の解、${x_1,x_2}$を求めてみます。

```math
\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}
```

## ライブラリのインポート

```python:ライブラリのインポート
from scipy.optimize import minimize
import numpy
```
まず、最適化に必要なライブラリをインポートしておきます。
`scipy`がインストールされていない場合は、コマンドプロンプトかAnaconda promptから`pip install scipy`コマンドを打つことでインストールすることができます。
## 目的関数と制約条件式の設定
```math
\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}
```
```python:目的関数と制約条件式の設定
#　目的関数
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]
```
等式制約条件式および不等式制約条件式については全て左辺に移項した式で定義する必要があります。

## 逐次二次計画法の適用
```python
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)$の値|

## 解の取り出し
また、得られた解についてはこのように変数に代入することで取り出せます。
```python:
x1,x2 = result.x[0],result.x[1]
print(f"x1 = {x1}",f"x2 = {x2}")
```

## 例題１のコードまとめ
<details><summary>ここをクリック</summary>

```python:逐次二次計画法の適用
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)
```
</details>


# 例題２【4変数関数の最小値を求める】
もう1問解いてみます。
解探索の初期点は、
$(x_1,x_2,x_3,x_4)=(1,5,5,1)$とします。

```math
\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}
```
## 逐次二次計画法の実装
```python:実装
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])
```
こちらも求めることができました。

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


# 参考文献
>* [しっかり学ぶ数理最適化 モデルからアルゴリズムまで](https://www.amazon.co.jp/%E3%81%97%E3%81%A3%E3%81%8B%E3%82%8A%E5%AD%A6%E3%81%B6%E6%95%B0%E7%90%86%E6%9C%80%E9%81%A9%E5%8C%96-%E3%83%A2%E3%83%87%E3%83%AB%E3%81%8B%E3%82%89%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0%E3%81%BE%E3%81%A7-KS%E6%83%85%E5%A0%B1%E7%A7%91%E5%AD%A6%E5%B0%82%E9%96%80%E6%9B%B8-%E6%A2%85%E8%B0%B7-%E4%BF%8A%E6%B2%BB/dp/4065212707/ref=sr_1_1?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&crid=1DX9M8YFIUPPB&keywords=%E6%A2%85%E8%B0%B7&qid=1649313566&sprefix=%E6%A2%85%E8%B0%B7%2Caps%2C234&sr=8-1)
>* [Pythonによる数理最適化入門 (実践Pythonライブラリー)](https://www.amazon.co.jp/Python%E3%81%AB%E3%82%88%E3%82%8B%E6%95%B0%E7%90%86%E6%9C%80%E9%81%A9%E5%8C%96%E5%85%A5%E9%96%80-%E5%AE%9F%E8%B7%B5Python%E3%83%A9%E3%82%A4%E3%83%96%E3%83%A9%E3%83%AA%E3%83%BC-%E4%B9%85%E4%BF%9D-%E5%B9%B9%E9%9B%84/dp/4254128959/ref=sr_1_3?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&crid=22UVL29AXJYU2&keywords=%E4%B9%85%E4%BF%9D%E5%B9%B9%E9%9B%84&qid=1649313664&sprefix=%E4%B9%85%E4%BF%9D%E5%B9%B9%E9%9B%84%2Caps%2C182&sr=8-3)
>* [非線形システムの最適化](https://www.amazon.co.jp/%E9%9D%9E%E7%B7%9A%E5%BD%A2%E3%82%B7%E3%82%B9%E3%83%86%E3%83%A0%E3%81%AE%E6%9C%80%E9%81%A9%E5%8C%96%E2%80%95%E4%B8%80%E7%9B%AE%E7%9A%84%E3%81%8B%E3%82%89%E5%A4%9A%E7%9B%AE%E7%9A%84%E3%81%B8-%E5%9D%82%E5%92%8C-%E6%AD%A3%E6%95%8F/dp/4627912501/ref=sr_1_12?crid=3EK0MSON2HCNV&keywords=%E5%9D%82%E5%92%8C%E6%AD%A3%E6%95%8F&qid=1649314324&sprefix=%E5%9D%82%E5%92%8C%2Caps%2C217&sr=8-12)
>* [SciPy Beginner's Guide for Optimization](https://www.youtube.com/watch?v=cXHvC_FGx24)