最適化
Julia
パズル
組合せ最適化
JuMP

最適化における Julia - JuMP で色々解く

More than 1 year has passed since last update.

Julia の最適化パッケージ JuMP の例を、引き続き紹介します。
過去記事 → 「最適化における Julia : JuMP事始め」, 「同: JuMP で N-Queen を解く

ゲーム理論を解く

Python - PuLP

ネタ元: Qiita記事 「組合せ最適でゲーム理論を解く@Tsutomu-KKE@github さん

from pulp import *
from ortoolpy import addvar, addvars

a = [[0, -1, 1], [4, 0, -1], [-1, 1, 0]]
m = LpProblem(sense=LpMaximize) # 数理モデル
xyz = addvars(3) # 変数 x,y,z
w = addvar(lowBound=None) # 変数 w
m += w # 目的関数
m += lpSum(xyz) == 1 # 制約条件
for i in range(3):
    m += lpDot(a[i], xyz) >= w # 制約条件
m.solve() # 求解
print(value(w), [value(v) for v in xyz])
>>>
0.16666667 [0.16666667, 0.33333333, 0.5]

Julia - JuMP

using JuMP
using Clp

a = [[0, -1, 1], [4, 0, -1], [-1, 1, 0]]
m = Model( solver=ClpSolver() ) # 数理モデル
@variable(m, x[1:3] >= 0.0) # 変数 x[1],x[2],x[3]
@variable(m, w)             # 変数 w
@objective(m, Max, w)       # 目的関数
@constraint(m, sum(x[i] for i=1:3) == 1.0) # 制約条件
@constraint(m, [i=1:3], dot(a[i],x) >= w)  # 制約条件
status=solve(m)

■ 2017/05/23 修正--- Model() を作成するときに solver=ClpSolver()を明示する様にしました 修正終わり---

ちなみに、制約条件をコマンドライン (REPL)から入力すると、以下の数式が表示されます。

julia> @constraint(m, sum(x[i] for i=1:3) == 1.0)
x[1] + x[2] + x[3] = 1

julia> @constraint(m, [i=1:3], dot(a[i],x) >= w)
3-element Array{JuMP.ConstraintRef,1}:
 -x[2] + x[3] - w  0 
 4 x[1] - x[3] - w  0
 -x[1] + x[2] - w  0 

最適化の結果

julia> status
:Optimal

julia> getobjectivevalue(m)
0.16666666666666669

julia> getvalue(x)
3-element Array{Float64,1}:
 0.166667
 0.333333
 0.5     

ナップサック問題

Python - PuLP

ネタ元: Qiita記事 「組合せ最適化を使おう - ナップサック問題@Tsutomu-KKE@github さん

from pulp import *
size = [21, 11, 15, 9, 34, 25, 41, 52]
weight = [22, 12, 16, 10, 35, 26, 42, 53]
capacity = 100
r = range(len(size))
m = LpProblem(sense=LpMaximize) # 数理モデル
x = [LpVariable('x%d'%i, cat=LpBinary) for i in r] # 変数
m += lpDot(weight, x) # 目的関数
m += lpDot(size, x) <= capacity # 制約
m.solve()
print((value(m.objective), [i for i in r if value(x[i]) > 0.5]))
>>>
(105.0, [0, 1, 3, 4, 5])

Julia - JuMP

using JuMP
using Cbc

size = [21, 11, 15, 9, 34, 25, 41, 52]
weight = [22, 12, 16, 10, 35, 26, 42, 53]
capacity = 100

n=length(size)
m=Model(solver=CbcSolver() )             # モデル
@variable(m, x[1:n], Bin)                # 変数
@objective(m, Max,dot(weight,x))         # 目的関数
@constraint(m, dot(size,x) <= capacity)  # 制約条件
status = solve(m)

■ 2017/05/23 修正--- solverClpSolver() から CbcSolver() に変更しました 修正終わり---

目的関数、制約条件の数式 (REPL, コマンドライン)

julia> @objective(m, Max,dot(weight,x))         # 目的関数
22 x[1] + 12 x[2] + 16 x[3] + 10 x[4] + 35 x[5] + 26 x[6] + 42 x[7] + 53 x[8]

julia> @constraint(m, dot(size,x) <= capacity)  # 制約条件
21 x[1] + 11 x[2] + 15 x[3] + 9 x[4] + 34 x[5] + 25 x[6] + 41 x[7] + 52 x[8]  100

最適化の結果

julia> status
:Optimal

julia> getobjectivevalue(m)
105.0

julia> getvalue(x)
8-element Array{Float64,1}:
 1.0
 1.0
 0.0
 1.0
 1.0
 1.0
 0.0
 0.0

以下続くかな

 いまのところ、JuMP のみで定式化できる例しか取り組んでいません。