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

Julia の最適化パッケージ JuMP の例を、引き続き紹介します。
過去記事 → 「最適化における Julia : [JuMP事始め](http://qiita.com/tenfu2tea/items/499ff93897f25a0256c0)」, 「[同: JuMP で N-Queen を解く](http://qiita.com/tenfu2tea/items/b1ad78cc9f0942fdbc93)」

# ゲーム理論を解く
## Python - PuLP
ネタ元: Qiita記事 「[組合せ最適でゲーム理論を解く](http://qiita.com/Tsutomu-KKE@github/items/eeb19fc086640562513e)」 @Tsutomu-KKE@github さん

```Python
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

```Julia
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
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
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記事 「[組合せ最適化を使おう - ナップサック問題](http://qiita.com/Tsutomu-KKE@github/items/bfbf4c185ed7004b5721#%E3%83%8A%E3%83%83%E3%83%97%E3%82%B5%E3%83%83%E3%82%AF%E5%95%8F%E9%A1%8C-1) 」@Tsutomu-KKE@github さん

```Python
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

```Julia
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
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
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 のみで定式化できる例しか取り組んでいません。