最近研究で最適化ソルバを使う機会が増えたので,そのメモ
(英語を読むことが苦手なので,毎回英語で説明を読むくらいなら日本語で書いておくという気持ち)
間違えているところがあれば指摘していただけると幸いです.
今回の記事では線形計画問題(LP)の解き方しか書かないが,CPLEXでは整数計画問題(IP), 混合整数計画問題(MIP), 凸二次計画問題(CQP), 凸二次制約つき凸二次計画問題(QCQP)など,様々な問題を解くことができる.
前提条件として,CPLEXはインストール済み,Pathも通してあるとする.
学生はAcademic版が無料で手に入る?
たとえば,以下の問題を解く時を考える:
\min 3 x_1 + 24 x_2 + 13 x_3 \\
\text{ subject to }
\begin{cases}
110 x_1 + 205x_2 + 160 x_3 \geq 2000 \\
4 x_1 + 32 x_2 + 13 x_3 \geq 55 \\
x1, x2, x3 \geq 0
\end{cases}
$\min$の後ろを「目的関数」,subject to の後ろを「制約式」という.
目的関数が線形,制約式も線形なので,これは線形計画問題(Linear Programming, LP)とよばれる問題の一種.
余談: この最適化問題が解をもつならば,双対定理より,双対ギャップ 0 で双対問題も解をもつことが知られている.
まずは,この問題をLPの標準形に直す.
LPの標準形は,以下の形で表される:
\min \mathbf{c} \cdot \mathbf{x} \\
\text{subject to }
\begin{cases}
\mathrm{A} \mathbf{x} \leq \mathbf{b}, \\
\mathbf{x} \geq \mathbf{0}
\end{cases}
ここで,$\mathrm{A} \in \mathbb{R}^{m \times n}, \mathbf{b} \in \mathbb{R}^{m}$は制約を行列の形で表したもので,不等号は各要素について成り立つという意味で用いられる:
\forall i \in \{ 1, 2, .., m \},\\
(\mathrm{A} \mathbf{x})_i \leq b_i
もともと解きたかった問題をこの形に直すと,
\min \begin{bmatrix} 3 & 24 & 13 \end{bmatrix}
\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} \\
\text{subject to }
\begin{cases}
\begin{bmatrix}
-110 & -205 & -160 \\
-4 & -32 & -13
\end{bmatrix}
\begin{bmatrix}
x_1 \\ x_2 \\ x_3
\end{bmatrix}
\leq
\begin{bmatrix}
-2000 \\ -55
\end{bmatrix} \\
x_1, x_2, x_3 \geq 0
\end{cases}
この式のうち,
c = \begin{bmatrix} 3 & 24 & 13 \end{bmatrix},
\mathbf{x} = \begin{bmatrix} x_1 & x_2 & x_3 \end{bmatrix},
\mathrm{A} = \begin{bmatrix}
-110 & -205 & -160 \\
-4 & -32 & -13
\end{bmatrix},
\mathbf{b} = \begin{bmatrix} -2000 \\ -55 \end{bmatrix}
と置けば,LPの標準形になることがわかる.
この問題をIBMの最適化ソルバCPLEXを用いて解く.
CPLEXを使う上で$ \mathrm{A} \mathbf{x} \leq \mathbf{b} $の形にする必要はないが,個人的にこの形が好きなのでこの形にした
import cplex
# ときたい最適化問題の定数部分を行列,ベクトルに変形する
# min c x
# s.t. A x ≤ b
# x ≥ 0
c = [3, 24, 13]
x = ["x{}".format(i) for i in range(1, 4)] # ["x1", "x2", "x3"]
A = [[-110, -205, -160], \
[-4, -32, -13]]
sense = ['L', 'L'] # Ax ≤ b の不等号部分の配列
# 'L' : (Ax)_i ≤ b_i のときに書く
# 'G' : (Ax)_i ≥ b_i のときに書く
# 'E' : (Ax)_i = b_i のときに書く
b = [-2000, -55]
lb = [0, 0, 0] #x1, x2, x3 ≥ 0 という意味
# インスタンスを作る
# Problem(問題)の頭文字をとって P とした
P = cplex.Cplex()
P.objective.set_sense(P.objective.sense.minimize) # 最小化問題か最大化問題かを表す文
# 最大化ならP.objective.sense.maximize
P.set_problem_type(P.problem_type.LP) # 線形計画問題であるという文
P.variables.add(obj=c, lb=lb, names=x) # 目的関数,変数を追加
for i in range(len(A)):
expr = cplex.SparsePair(ind=x, val=A[i])
P.linear_constraints.add(lin_expr=expr, senses=sense[i], rhs=[b[I]]) # 線形制約を一つずつ追加
P.solve() # 最適化問題を解く
P.solution.get_values() # 最適値を表示
P.solution.get_objective_value() # 最適解 [x1, x2, x3] を表示
P.solution.get_dual_values() # 双対問題の最適解を表示
CPLEXではデフォルトで途中計算がコンソールに表示されるが,これをなくすには
P.set_results_stream(results_file=None)
と書けば良い.
おわりに
気が向いたり,要望があればQPを解くとき等の書き方ものせます.