最適化問題をPythonのPyomoライブラリーで解こうその2：線形計画問題（行列とベクトルを利用）

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

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

基本となるのが、線形計画問題です。線形式しか登場しません。前回、PythonのPyomoライブラリーで解く方法について簡単に説明しました。

前回の例で登場した最適化問題は小規模のものでしたが、現実は変数や数式の数は数百、数千、数万になります。

その場合、今回のように1つ1つ数式を定義していくことは現実的ではありません。

例えば、数式の係数や変数の上限や下限の制約を、外部からデータとして取得し、設定することが多いです。

今回は、「行列とベクトルを使った線形計画問題」というお話しをします。前回と同じ最適化問題を扱います。

前回の復習です。

Contents [hide]

前回の復習
　最適化の構成要素
　モデリングツールとソルバー
　最適化問題のモデリングし解くまでの流れ
最適化する問題（前回と同じ）
前回と異なるところ
必要なモジュールの読み込み
データ（目的関数や制約式の係数など）
最適解の計算
　変数に制約なし
　変数に非負制約あり
　変数に範囲制約あり
　変数に整数制約と範囲制約あり
まとめ

前回の復習

　最適化の構成要素

最適化には3つの構成要素があります。

目的関数
変数
制約条件

「ある制約条件のもと目的関数を最適化（最大化 or 最小化）する変数を見つける」のが最適化問題です。

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

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

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

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

GLPK（GNU Linear Programming Kit)は、線形計画問題を解くための無料のソルバーです。

　最適化問題のモデリングし解くまでの流れ

以下、最適化問題を解くまでの流れです。

モデリング
1. モデルのインスタンスの生成
2. 変数の定義
3. 目的関数の定義
4. 制約条件の定義
最適化
1. ソルバーの設定
2. 最適化の実施
結果確認

非常にシンプルです。

最適化する問題（前回と同じ）

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

\displaystyle<br />\begin{aligned}<br />& \underset{x} \text{min} && f(x_1,x_2) = 75 x_1 + 125 x_2 \\<br />& \text{subject to} && 6x_1 + 3x_2 \ge 38 \\<br />&&& 5x_1 + 21x_2 \ge 29\\<br />\end{aligned}<br />

ベクトルや行列で表現すると、以下のようになります。

\displaystyle<br />\begin{aligned}<br />& \underset{x} \text{min} && f(\bm{x}) = \bm{\alpha}^T \bm{x} \\<br />& \text{subject to} && \bm{A x} \ge \bm{\beta} \\<br />\end{aligned}<br />

\displaystyle<br />\bm{x}=\begin{bmatrix}<br />x_1 \\<br />x_2<br />\end{bmatrix}、<br />\bm{\alpha}=\begin{bmatrix}<br />75 \\<br />125<br />\end{bmatrix}、<br />\bm{A}=\begin{bmatrix}<br />6 & 3 \\<br />5 & 21<br />\end{bmatrix}、<br />\bm{\beta}=\begin{bmatrix}<br />38 \\<br />29<br />\end{bmatrix}<br />

さらに、変数に幾つかの制約をつけた場合、どうなるのかを見ていきます。

今回は、変数への制約の付け方を、次の3パターンで設定し見ていきます。

変数に非負制約あり
変数に範囲制約あり
変数に整数制約と範囲制約あり

前回と異なるところ

目的関数や制約式の係数などを、行列およびベクトル（正確には、NumPyの配列（array））として与えられたとして、前回と同じ最適化問題を解いていきます。

NumPyの配列（array）として与えられたデータを、辞書(dictionary)へ変換し、Pyomoでモデリングします。

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

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

以下、コードです。

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([75,125])

## 係数行列

t_data = np.array([[6, 3],

[5, 21]])

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

t0_data = np.array([38,29])

# 辞書(dictionary)へ変換

## 目的変数の係数

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

## 係数行列

t = dict(

((i, j), t_data[i-1][j-1])

for i in range(1, I+1)

for j in range(1, I+1)

)

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

t0 = dict((i, t0_data[i-1]) for i in range(1,I+1))

# # 変数の数 # I = 2 # # データ（係数など） # ## 目的変数の係数 a_data = np.array([75,125]) ## 係数行列 t_data = np.array([[6, 3], [5, 21]]) ## ベクトル（各制約式の下限） t0_data = np.array([38,29]) # # 辞書(dictionary)へ変換 # ## 目的変数の係数 a = dict((i, a_data[i-1]) for i in range(1, I+1)) ## 係数行列 t = dict( ((i, j), t_data[i-1][j-1]) for i in range(1, I+1) for j in range(1, I+1) ) ## ベクトル（各制約式の下限） t0 = dict((i, t0_data[i-1]) for i in range(1,I+1))

#
# 変数の数
#

I = 2

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

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

## 係数行列
t_data = np.array([[6, 3],
                   [5, 21]])

## ベクトル（各制約式の下限）
t0_data = np.array([38,29])

#
# 辞書(dictionary)へ変換
#

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

## 係数行列
t = dict(
    ((i, j), t_data[i-1][j-1])  
    for i in range(1, I+1) 
    for j in range(1, I+1)
)

## ベクトル（各制約式の下限）
t0 = dict((i, t0_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, within=pyo.Reals)

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

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

# 変数
model.x = pyo.Var(model.I, within=pyo.Reals)

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

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

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

以下、コードです。

# 目的関数の数式の定義

def ObjRule(model):

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

# 目的関数として設定

model.obj = pyo.Objective(rule = ObjRule,

sense = pyo.minimize)

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

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

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

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

\displaystyle<br />\begin{aligned}<br />&&& 6x_1 + 3x_2 \ge 38 \\<br />&&& 5x_1 + 21x_2 \ge 29\\<br />\end{aligned}<br />

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

以下、コードです。

# 制約条件の数式の定義

def Construle(model, i):

return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

# 制約条件として設定

model.eq = pyo.Constraint(model.I, rule = Construle)

# 制約条件の数式の定義 def Construle(model, i): return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i] # 制約条件として設定 model.eq = pyo.Constraint(model.I, rule = Construle)

# 制約条件の数式の定義
def Construle(model, i):
    return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

# 制約条件として設定
model.eq = pyo.Constraint(model.I, rule = Construle)

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

以下、コードです。

# 最適化

## ソルバーの設定

opt = pyo.SolverFactory('glpk')

## 最適化の実施

res = opt.solve(model)

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

# 最適化

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

## 最適化の実施
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 = 6.40540540540541\\<br />&&& x_2 = -0.144144144144144\\<br />\end{aligned}<br />

最適値（最小値）は……

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

　変数に非負制約あり

変数に非負制約をつけた場合です。

# モデルのインスタンスの生成

model = pyo.ConcreteModel()

# 変数の定義

# 添字

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

# 変数

model.x = pyo.Var(model.I, domain=pyo.NonNegativeReals)

# 目的関数

# 目的関数の数式の定義

def ObjRule(model):

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

# 目的関数として設定

model.obj = pyo.Objective(rule = ObjRule,

sense = pyo.minimize)

# 制約条件

# 制約条件の数式の定義

def Construle(model, i):

return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

# 制約条件として設定

model.eq1 = pyo.Constraint(model.I, rule = Construle)

# 最適化

# ソルバーの設定

opt = SolverFactory('glpk')

# 最適化の実施

res = opt.solve(model)

# 結果確認

print(model.display())

print('\n')

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

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

# # モデルのインスタンスの生成 # model = pyo.ConcreteModel() # # 変数の定義 # # 添字 model.I = pyo.Set(initialize=range(1, I+1)) # 変数 model.x = pyo.Var(model.I, domain=pyo.NonNegativeReals) # # 目的関数 # # 目的関数の数式の定義 def ObjRule(model): return sum(a[i] * model.x[i] for i in model.I) # 目的関数として設定 model.obj = pyo.Objective(rule = ObjRule, sense = pyo.minimize) # # 制約条件 # # 制約条件の数式の定義 def Construle(model, i): return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i] # 制約条件として設定 model.eq1 = pyo.Constraint(model.I, rule = Construle) # # 最適化 # # ソルバーの設定 opt = SolverFactory('glpk') # 最適化の実施 res = opt.solve(model) # # 結果確認 # print(model.display()) print('\n') print('optimum value = ', model.obj()) print("x = ", model.x[:]())

#
# モデルのインスタンスの生成
#

model = pyo.ConcreteModel()

#
# 変数の定義
#

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

# 変数
model.x = pyo.Var(model.I, domain=pyo.NonNegativeReals)

#
# 目的関数
#

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

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

#
# 制約条件
#

# 制約条件の数式の定義
def Construle(model, i):
    return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

# 制約条件として設定
model.eq1 = pyo.Constraint(model.I, rule = Construle)

#
# 最適化
#

# ソルバーの設定
opt = SolverFactory('glpk')

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

#
# 結果確認
#

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

以下、コードです。

以下、実行結果です。

　変数に範囲制約あり

変数の上限と下限のベクトルを設定します。こちらも、ベクトル（正確には、NumPyの配列（array））として設定し、それを辞書(dictionary)へ変換します。

以下、コードです。

# 上限と下限のデータの追加

# データ（係数など）

## 変数の上限

upper_data = np.array([5,10])

## 変数の下限

lower_data = np.array([-5,0])

# 辞書(dictionary)へ変換

## 変数の上限

upper = dict((i, upper_data[i-1]) for i in range(1,I+1))

## 変数の下限

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

# # 上限と下限のデータの追加 # # データ（係数など） ## 変数の上限 upper_data = np.array([5,10]) ## 変数の下限 lower_data = np.array([-5,0]) # 辞書(dictionary)へ変換 ## 変数の上限 upper = dict((i, upper_data[i-1]) for i in range(1,I+1)) ## 変数の下限 lower = dict((i, lower_data[i-1]) for i in range(1,I+1))

#
# 上限と下限のデータの追加
#

# データ（係数など）

## 変数の上限
upper_data = np.array([5,10])

## 変数の下限
lower_data = np.array([-5,0])

# 辞書(dictionary)へ変換

## 変数の上限
upper = dict((i, upper_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.I = pyo.Set(initialize=range(1, I+1))

# 変数

model.x = pyo.Var(model.I, domain=pyo.NonNegativeReals)

# 目的関数

# 目的関数の数式の定義

def ObjRule(model):

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

# 目的関数として設定

model.obj = pyo.Objective(rule = ObjRule,

sense = pyo.minimize)

# 制約条件

# 制約1

def Construle1(model, i):

return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

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

# 制約2

def Construle2(model, i):

return model.x[i] <= upper[i]

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

# 制約3

def Construle3(model, i):

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

model.eq3 = pyo.Constraint(model.I, rule = Construle3)

# 最適化

# ソルバーの設定

opt = SolverFactory('glpk')

# 最適化の実施

res = opt.solve(model)

# 結果確認

print(model.display())

print('\n')

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

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

# # モデルのインスタンスの生成 # model = pyo.ConcreteModel() # # 変数の定義 # # 添字 model.I = pyo.Set(initialize=range(1, I+1)) # 変数 model.x = pyo.Var(model.I, domain=pyo.NonNegativeReals) # # 目的関数 # # 目的関数の数式の定義 def ObjRule(model): return sum(a[i] * model.x[i] for i in model.I) # 目的関数として設定 model.obj = pyo.Objective(rule = ObjRule, sense = pyo.minimize) # # 制約条件 # # 制約1 def Construle1(model, i): return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i] model.eq1 = pyo.Constraint(model.I, rule = Construle1) # 制約2 def Construle2(model, i): return model.x[i] <= upper[i] model.eq2 = pyo.Constraint(model.I, rule = Construle2) # 制約3 def Construle3(model, i): return model.x[i] >= lower[i] model.eq3 = pyo.Constraint(model.I, rule = Construle3) # # 最適化 # # ソルバーの設定 opt = SolverFactory('glpk') # 最適化の実施 res = opt.solve(model) # # 結果確認 # print(model.display()) print('\n') print('optimum value = ', model.obj()) print("x = ", model.x[:]())

#
# モデルのインスタンスの生成
#

model = pyo.ConcreteModel()

#
# 変数の定義
#

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

# 変数
model.x = pyo.Var(model.I, domain=pyo.NonNegativeReals)

#
# 目的関数
#

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

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

#
# 制約条件
#

# 制約1
def Construle1(model, i):
    return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

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

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

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

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

model.eq3 = pyo.Constraint(model.I, rule = Construle3)

#
# 最適化
#

# ソルバーの設定
opt = SolverFactory('glpk')

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

#
# 結果確認
#

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

以下、実行結果です。

　変数に整数制約と範囲制約あり

変数に整数制約と範囲制約をつけた場合です。

以下、コードです。

# モデルのインスタンスの生成

model = pyo.ConcreteModel()

# 変数の定義

# 添字

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

# 変数

model.x = pyo.Var(model.I, domain=pyo.Integers)

# 目的関数

# 目的関数の数式の定義

def ObjRule(model):

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

# 目的関数として設定

model.obj = pyo.Objective(rule = ObjRule,

sense = pyo.minimize)

# 制約条件

# 制約1

def Construle1(model, i):

return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

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

# 制約2

def Construle2(model, i):

return model.x[i] <= upper[i]

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

# 制約3

def Construle3(model, i):

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

model.eq3 = pyo.Constraint(model.I, rule = Construle3)

# 最適化

# ソルバーの設定

opt = SolverFactory('glpk')

# 最適化の実施

res = opt.solve(model)

# 結果確認

print(model.display())

print('\n')

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

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

# # モデルのインスタンスの生成 # model = pyo.ConcreteModel() # # 変数の定義 # # 添字 model.I = pyo.Set(initialize=range(1, I+1)) # 変数 model.x = pyo.Var(model.I, domain=pyo.Integers) # # 目的関数 # # 目的関数の数式の定義 def ObjRule(model): return sum(a[i] * model.x[i] for i in model.I) # 目的関数として設定 model.obj = pyo.Objective(rule = ObjRule, sense = pyo.minimize) # # 制約条件 # # 制約1 def Construle1(model, i): return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i] model.eq1 = pyo.Constraint(model.I, rule = Construle1) # 制約2 def Construle2(model, i): return model.x[i] <= upper[i] model.eq2 = pyo.Constraint(model.I, rule = Construle2) # 制約3 def Construle3(model, i): return model.x[i] >= lower[i] model.eq3 = pyo.Constraint(model.I, rule = Construle3) # # 最適化 # # ソルバーの設定 opt = SolverFactory('glpk') # 最適化の実施 res = opt.solve(model) # # 結果確認 # print(model.display()) print('\n') print('optimum value = ', model.obj()) print("x = ", model.x[:]())

#
# モデルのインスタンスの生成
#

model = pyo.ConcreteModel()

#
# 変数の定義
#

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

# 変数
model.x = pyo.Var(model.I, domain=pyo.Integers)

#
# 目的関数
#

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

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

#
# 制約条件
#

# 制約1
def Construle1(model, i):
    return sum(t[i,j] * model.x[j] for j in model.I) >= t0[i]

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

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

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

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

model.eq3 = pyo.Constraint(model.I, rule = Construle3)

#
# 最適化
#

# ソルバーの設定
opt = SolverFactory('glpk')

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

#
# 結果確認
#

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

以下、実行結果です。