最適化問題をPythonのPyomoライブラリーで解こうその3：非線形計画問題

最適化問題は、マーケティング予算配分の最適化、配送ルートの最適化、スケジュール最適化など、何かを最適化する問題を扱うものです。

最適化問題には、登場する数式や最適解の条件などによって、線形計画問題や非線形計画問題、混合整数計画問題などがあります。

基本となるのが、線形計画問題です。線形式しか登場しません。前回、ベクトルや行列による数式表現のケースで簡単に説明しました。

実務では、非線形な関数を扱うことも多いです。

今回は、「Pyomoでの非線形計画問題を解く方法」についてお話しをします。

Contents [hide]

モデリングツールとソルバー
ソルバーのインストール
今回、最適化する問題
必要なモジュールの読み込み
データ（目的関数や制約式の係数など）
最適解の計算
まとめ

モデリングツールとソルバー

最適化問題を解くには、次の2つのツールが必要になります。

モデリングツール；最適化問題を記述するためのツール
ソルバー：記述された最適化問題を計算するためのツール

Pyomoは、モデリングツールです。

そのため、ソルバーを別途準備する必要があります。

線形計画問題を解くための無料のソルバーの1つが、GLPK（GNU Linear Programming Kit)です。前回利用しました。

今回は、非線形計画問題を解くためのソルバーが必要になります。

IPOPT
Bonmin
Couenne

Bonminは、Windowsには対応していません。Couenneは、Pythonからインストールできず別途インストールする必要があります。

ということで、今回はIPOPT （Interior Point OPTimizer）というソルバーを使います。

IPOPTは主双対内点法を利用したソルバーで、非線形最適化問題に対応しています。

幾つか問題があります。大きいところでは、以下の3点です。

大域的最適解に弱い（凸計画問題であれば問題ない）
離散に弱い（混合整数非線形計画問題用のソルバーでないため）
大規模問題に弱い（高速ではなく中速）

ちなみに、BonminとCouenneは混合整数非線形計画問題に対応しています。

ソルバーのインストール

IPOPTをcondaでインストールするときは、以下です。

conda install -c conda-forge ipopt

conda install -c conda-forge ipopt

pipでもインストールできます。以下です。

pip install ipopt

pip install ipopt

おまけに、BonminとCouenneのインストール方法も記載しておきます。

以下は、Bonminのインストール方法です。LinuxとmacOSに対応しWindowsに対応していません。

conda install -c conda-forge coinbonmin

conda install -c conda-forge coinbonmin

Windows版のBonminをインストールする場合は、以下のサイトからダウンロードしインストールして下さい。Couenneもここからダウンロードできます。

Download COIN-OR binary
https://www.coin-or.org/download/binary/

今回、最適化する問題

今回、最適化する問題です。

\displaystyle<br />\begin{aligned}<br />& \underset{x} \text{min} && f(x_1,x_2) = 4 x_1^2 + 0.25 x_2^2 \\<br />& \text{subject to} && x_1 x_2 \ge 1 \\<br />&&& x_1 \ge 0.5\\<br />&&& x_2 \ge 2\\<br />\end{aligned}<br />

必要なモジュールの読み込み

先ずは、利用するモジュールを読み込みます。

以下、コードです。

import numpy as np

import pyomo.environ as pyo

from pyomo.opt import SolverFactory

import numpy as np import pyomo.environ as pyo from pyomo.opt import SolverFactory

import numpy as np
import pyomo.environ as pyo
from pyomo.opt import SolverFactory

データ（目的関数や制約式の係数など）

目的関数や制約式の係数などを、行列およびベクトル（正確には、NumPyの配列（array））として設定し、それを辞書(dictionary)へ変換します。

以下、コードです。

# 変数の数

I = 2

# データ（係数など）

## 目的変数の係数

a_data = np.array([4, 0.25])

## ベクトル（各制約式の下限）

lower_data = np.array([0.5, 2])

# 辞書(dictionary)へ変換

## 目的変数の係数

a = dict((i, a_data[i-1]) for i in range(1,I+1))

## ベクトル（各制約式の下限）

lower = dict((i, lower_data[i-1]) for i in range(1,I+1))

# # 変数の数 # I = 2 # # データ（係数など） # ## 目的変数の係数 a_data = np.array([4, 0.25]) ## ベクトル（各制約式の下限） lower_data = np.array([0.5, 2]) # # 辞書(dictionary)へ変換 # ## 目的変数の係数 ## 目的変数の係数 a = dict((i, a_data[i-1]) for i in range(1,I+1)) ## ベクトル（各制約式の下限） lower = dict((i, lower_data[i-1]) for i in range(1,I+1))

#
# 変数の数
#

I = 2

#
# データ（係数など）
#

## 目的変数の係数
a_data = np.array([4, 0.25])

## ベクトル（各制約式の下限）
lower_data = np.array([0.5, 2])

#
# 辞書(dictionary)へ変換
#

## 目的変数の係数

## 目的変数の係数
a = dict((i, a_data[i-1]) for i in range(1,I+1))

## ベクトル（各制約式の下限）
lower = dict((i, lower_data[i-1]) for i in range(1,I+1))

最適解の計算

先ず、モデリングです。

モデルのインスタンスを生成します。

以下、コードです。

model = pyo.ConcreteModel()

model = pyo.ConcreteModel()

変数を定義します。2要素のベクトル $x$ です。

以下、コードです。

# 変数の添字

model.I = pyo.Set(initialize=range(1, I+1))

# 変数の定義

model.x = pyo.Var(model.I)

# 変数の添字 model.I = pyo.Set(initialize=range(1, I+1)) # 変数の定義 model.x = pyo.Var(model.I)

# 変数の添字
model.I = pyo.Set(initialize=range(1, I+1))

# 変数の定義
model.x = pyo.Var(model.I)

次の目的関数を定義します。

\displaystyle<br />\begin{aligned}<br />&& f(x_1,x_2) = 4x_1^2 + 0.25x_2^2<br />\end{aligned}<br />

目的関数を数式として表現し関数として定義した後に、その関数を目的関数に設定します。

以下、コードです。

# 目的関数の数式の定義

def ObjRule(model):

return sum(a[i] * (model.x[i])**2 for i in model.I)

# 目的関数として設定

model.obj = pyo.Objective(rule = ObjRule, sense = pyo.minimize)

# 目的関数の数式の定義 def ObjRule(model): return sum(a[i] * (model.x[i])**2 for i in model.I) # 目的関数として設定 model.obj = pyo.Objective(rule = ObjRule, sense = pyo.minimize)

# 目的関数の数式の定義
def ObjRule(model):
    return sum(a[i] * (model.x[i])**2 for i in model.I)

# 目的関数として設定
model.obj = pyo.Objective(rule = ObjRule, sense = pyo.minimize)

次の制約条件を定義します。

\displaystyle<br />\begin{aligned}<br />&&& x_1 x_2 \ge 1 \\<br />&&& x_1 \ge 0.5, x_2 \ge 2\\<br />\end{aligned}<br />

制約条件を数式として表現し関数として定義した後に、その関数を制約条件に設定します。

以下、コードです。

# 制約1

def Construle1(model):

return model.x[1]*model.x[2] >= 1

model.eq1 = pyo.Constraint(rule = Construle1)

# 制約2

def Construle2(model, i):

return model.x[i] >= lower[i]

model.eq2 = pyo.Constraint(model.I, rule = Construle2)

# 制約1 def Construle1(model): return model.x[1]*model.x[2] >= 1 model.eq1 = pyo.Constraint(rule = Construle1) # 制約2 def Construle2(model, i): return model.x[i] >= lower[i] model.eq2 = pyo.Constraint(model.I, rule = Construle2)

# 制約1
def Construle1(model):
    return model.x[1]*model.x[2] >= 1

model.eq1 = pyo.Constraint(rule = Construle1)

# 制約2
def Construle2(model, i):
    return model.x[i] >= lower[i]

model.eq2 = pyo.Constraint(model.I, rule = Construle2)

モデリングが終了したので、ソルバーを設定し最適化します。

以下、コードです。

# ソルバーの設定

opt = pyo.SolverFactory('ipopt')

# 最適化の実施

res = opt.solve(model)

# ソルバーの設定 opt = pyo.SolverFactory('ipopt') # 最適化の実施 res = opt.solve(model)

# ソルバーの設定
opt = pyo.SolverFactory('ipopt')

# 最適化の実施
res = opt.solve(model)

結果を確認します。

以下、コードです。

print(model.display())

print('\n')

print('optimum value = ', model.obj())

print("x = ", model.x[:]())

print(model.display()) print('\n') print('optimum value = ', model.obj()) print("x = ", model.x[:]())

print(model.display())
print('\n')
print('optimum value = ', model.obj())
print("x = ", model.x[:]())

以下、実行結果です。

最適解は……

\displaystyle<br />\begin{aligned}<br />&&& x_1 = 0.49999999558007313\\<br />&&& x_2 = 2.0000000022840774\\<br />\end{aligned}<br />

最適値（最小値）は……

\displaystyle<br />\begin{aligned}<br />&& f(x_1,x_2) = 1.99999998460437<br />\end{aligned}<br />

まとめ

今回は、「Pyomoでの非線形計画問題を解く方法」というお話しをしました。

次回以降から、具体的な事例を扱っていきます。