最適化問題は、マーケティング予算配分の最適化、配送ルートの最適化、スケジュール最適化など、何かを最適化する問題を扱うものです。
最適化問題には、登場する数式や最適解の条件などによって、線形計画問題や非線形計画問題、混合整数計画問題などがあります。
基本となるのが、線形計画問題です。線形式しか登場しません。前回、ベクトルや行列による数式表現のケースで簡単に説明しました。
実務では、非線形な関数を扱うことも多いです。
今回は、「Pyomoでの非線形計画問題を解く方法」についてお話しをします。
モデリングツールとソルバー
最適化問題を解くには、次の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
pipでもインストールできます。以下です。
pip install ipopt
おまけに、BonminとCouenneのインストール方法も記載しておきます。
以下は、Bonminのインストール方法です。LinuxとmacOSに対応しWindowsに対応していません。
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
データ（目的関数や制約式の係数など）
目的関数や制約式の係数などを、行列およびベクトル（正確には、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))
最適解の計算
先ず、モデリングです。
モデルのインスタンスを生成します。
以下、コードです。
model = pyo.ConcreteModel()
変数を定義します。2要素のベクトルxです。
以下、コードです。
# 変数の添字 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)
次の制約条件を定義します。
\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)
モデリングが終了したので、ソルバーを設定し最適化します。
以下、コードです。
# ソルバーの設定
opt = pyo.SolverFactory('ipopt')
# 最適化の実施
res = opt.solve(model)
結果を確認します。
以下、コードです。
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での非線形計画問題を解く方法」というお話しをしました。
次回以降から、具体的な事例を扱っていきます。
