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 修正--- solver を ClpSolver() から 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 のみで定式化できる例しか取り組んでいません。