JuMP による数理最適化超入門

●はじめに

「数理最適化 (Mathematical Optimization)」あるいは「数理計画法 (Mathematical Programing)」は、変数に関する不等式や等式で表される制約の条件下で、目的の関数を最小 (あるいは最大) にする変数の値を求める問題です。数理最適化は、問題を解くプログラム (最適化ソルバー) が商用非商用とわず開発されていて、それを使って解くのが一般的です。もちろん商用のソルバーのほうが高性能なのですが、近年ソルバーの性能は著しく向上していて、非商用 (無料) のソルバーでもかなりの問題を解くことができるようです。

ソルバーは数理モデル (変数、制約条件、目的関数など) を入力とし、それを解いた結果 (最適値) を出力します。数理モデルの記述には標準的なフォーマットがあるようですが、それらを自分で書くのは大変なので、モデリング言語を使うのが一般的です。たとえば、Python のモジュール PuLP を使うと、Python で変数や数式を記述することができます。PuLP については、拙作のページ Python3 Programming : PuLP による数理最適化超入門で説明しています。

Julia で数理最適化を行う際は、JuMP.jl という非常に優秀な共通モデリングライブラリを使用するのが標準です。JuMP で利用できるオープンソースのソルバーはいくつかありますが、今回は HiGHS (ハイズ) を使うことにしましょう。HiGHS は従来のオープンソースソルバー (CBC や GLPK など) を大きく上回る性能を持つと言われていて、現在は世界中で標準的なツールとして採用されています。

本ページでは JuMP と HiGHS を使って数理最適化の問題にチャレンジしてみようと思います。なお、M.Hiroi は JuMP を使うのは初めてなので、何か勘違いや間違いがあると思います。お気づきの点がありましたら、ご指摘いただけると助かります。とりあえず、拙作のページ PuLP による数理最適化超入門で取り上げた問題を JuMP で解くことができればいいな、と思っています。たいしたことはできませんが、よろしければお付き合いくださいませ。

●インストール

Julia のパッケージモード (REPL で ] キーを押す) で次のコマンドを実行します。

(@v1.12) pkg> add JuMP HiGHS

これで JuMP と HiGHS がインストールされます。インストール終了後 st で状態を表示すると、M.Hiroi の環境では次のように表示されます。

(@v1.12) pkg> st
Status `~/.julia/environments/v1.12/Project.toml`
  [336ed68f] CSV v0.10.16
  [a93c6f00] DataFrames v1.8.2
  [31c24e10] Distributions v0.25.126
  [87dc4568] HiGHS v1.23.0
  [4076af6c] JuMP v1.30.1
  [91a5bcdd] Plots v1.41.6
  [2913bbd2] StatsBase v0.34.11
  [f3b207a7] StatsPlots v0.15.8
  [bd369af6] Tables v1.12.1

インストールが終了したら、Backspace または Ctrl + C を押してください。通常の julia> プロンプトに戻ります。

●JuMP の基礎知識

JuMP の基本的な使い方は簡単です。手順は次のようになります。

model = Model(solver_constructor) で数理モデルを格納するオブジェクトを生成する

solver_constructor : HiGHS.Optimizer などのソルバーを指定する

@variable(model, name, type) で変数を定義する
- name : 変数名, 範囲 (上限・下限) 指定もできる
- type : 型 (デフォルトは実数)
- 例 : @variable(model, x >= 0, Int) （0以上の整数変数 \(x\)）
@objective(model, sense, expression) で目的関数 expression を設定する

sense : 最大化 (Max) または最小化 (Min) を指定する
expression には @variable で生成した変数を使用する

@constraint(model, [name], expression) で制約条件 expression を追加する

条件に名前 name を付けることで、後から参照や変更が可能
expression には @variable で生成した変数を使用する
制約条件は左辺式演算子右辺式の形式で指定する
使用できる演算子は <=. ==, >= の 3 つだけ (それ以外の演算子は NG)
制約条件はいくつでも指定することができる
expression には以下の関数が使用できる

sum(list), list の要素の合計値を表す
LinearAlgebra.dot(list1, list2), list1 と list2 の内積を表す
list の大きさを N とすると、sum(list1[i] * list2[i] for i = 1 : N) のほうが高速とのこと

println(model) で数理モデルを表示できる
optimize!(model) で構築した数理モデルをソルバーに渡し、最適化計算を実行する

termination_status(model), ソルバーの終了状態 (ステータス) を返す
最適解が見つかった場合は OPTIMAL、解が存在しない場合は INFEASIBLE などが返ってくる
objective_value(model), 最適化された目的関数の値 (最大値または最小値) を取得する
value(variable), 得られた変数の具体的な値 (最適解)

このほかに、数式に名前を付けるマクロ @expression や、配列に変数や数式を格納することも簡単にできます。

@expression(model, name, expr)
model 内で再利用可能な数式を作成する (name で参照可)
@variable(model, array[1:N] >= 0)
配列 array (大きさ N) に下限値 0 の変数を指定する
多次元配列でも OK
@constraint(model, [i = 1 : N], array[i] <= M)
N 個の制約を一括生成する
第 2 引数で name[ ... ] とすると、制約に名前 name を付けることもできる

それでは実際に JuMP を使って線形計画法を解いてみましょう。問題は参考文献『Ｃ言語による最新アルゴリズム事典』の線形計画法の例題をお借りしました。

目的関数:
    z = x + y + 1
制約条件:
    3 * x + 5 * y <= 15
    2 * x + y >= 4
    x - y == 1
    x >= 0
    y >= 0

上記の条件で z の最小値を求めます。プログラムは次のようになります。

リスト : 線形計画法 (test01.jl)

using JuMP, HiGHS                 

# モデルの生成
model = Model(HiGHS.Optimizer)    

# 変数の定義 (下限を 0 に設定)
@variable(model, x >= 0)          
@variable(model, y >= 0)          

# 目的関数
@objective(model, Min, x + y + 1) 

# 制約条件
@constraint(model, 3x + 5y <= 15) 
@constraint(model, 2x + y >= 4)   
@constraint(model, x - y == 1)    

# 実行
optimize!(model)

# 結果の確認
println("最適化ステータス: ", termination_status(model))
println("目的関数の最大値: ", objective_value(model))
println("最適解 x = ", value(x))
println("最適解 y = ", value(y))

@variable で変数 x, y を定義します。x > = 0, y > = 0 は変数の値が非負であることを表しています。もちろん、@constraint で条件を追加してもかまいません。あとは x と y を使って目的関数と制約条件を追加していくだけです。数式の書き方は Julia と同じなので難しいところはないでしょう。最後に、optimize!(model) でソルバーを実行して、結果を表示します。

実行結果は次のようになります。

$ julia test01.jl
Running HiGHS 1.14.0 (git hash: 7df0786de3): Copyright (c) 2026 under Apache 2.0 license terms
Using BLAS: blastrampoline
LP has 3 rows; 2 cols; 6 nonzeros
Coefficient ranges:
  Matrix  [1e+00, 5e+00]
  Cost    [1e+00, 1e+00]
  Bound   [0e+00, 0e+00]
  RHS     [1e+00, 2e+01]
Presolving model
0 rows, 0 cols, 0 nonzeros 0s
0 rows, 0 cols, 0 nonzeros 0s
Presolve reductions: rows 0(-3); columns 0(-2); nonzeros 0(-6) - Reduced to empty
Performed postsolve
Solving the original LP from the solution after postsolve

Model status        : Optimal
Objective value     :  3.3333333333e+00
P-D objective error :  1.1584935909e-16
HiGHS run time      :          0.01
最適化ステータス: OPTIMAL
目的関数の最大値: 3.333333333333334
最適解 x = 1.6666666666666667
最適解 y = 0.6666666666666667

実行環境 : Julia 1.12.6, JuMP 1.14.0, Ubunts 24.04 LTS (WSL2, Windows 11), intel CORE i5-1235U 1.30GHz

数式を入力するだけで、簡単に解を求めることができました。なお、HiGHS のログを非表示にする場合は、optimize!(model) を呼び出す前に、関数 set_silent(model) を実行してください。また、println(model) で数理モデルを表示することができます。

$ julia test01.jl
Min x + y + 1
Subject to
 x - y = 1
 2 x + y ≥ 4
 3 x + 5 y ≤ 15
 x ≥ 0
 y ≥ 0

最適化ステータス: OPTIMAL
目的関数の最大値: 3.333333333333334
最適解 x = 1.6666666666666667
最適解 y = 0.6666666666666667

●線形計画法の典型的な例題

次は、線形計画法の典型的な例題を 3 つ取り上げます。

[栄養問題]

下表に示す 3 種類の食品 (a, b, c) を使って、2 種類の栄養素 (x, y) の摂取量を満たす一番安い組み合わせを求めてください。

表 : 食品の栄養素
	食品a	食品b	食品c	摂取量
栄養素x	3	1	2	15
栄養素y	1	2	4	10
単価	4	2	5

食品の購入量を変数 a, b, c で表すことにすると、目的関数は次式のようになります。

目的関数:
    z = 4 * a + 2 * b + 5 * c

この式が最小となる各変数の値を求めます。一日に摂取すべき栄養素には下限値があるので、これが制約条件となります。

制約条件:
    3 * a + 1 * b + 2 * c >= 15,   栄養素 x の制約条件
    1 * a + 2 * b + 4 * c >= 10,   栄養素 y の制約条件
    a, b, c >= 0,                  非負条件

あとはこの式を JuMP でプログラムするだけです。

リスト : 栄養問題 (sample01.jl)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 変数の定義
@variable(model, a >= 0)
@variable(model, b >= 0)
@variable(model, c >= 0)

# 目的関数
@objective(model, Min, 4a + 2b + 5c)

# 制約条件
@constraint(model, 3a + b + 2c >= 15)
@constraint(model, a + 2b + 4c >= 10)
println(model)

# 実行
set_silent(model)
optimize!(model)

# 結果の表示
println("結果: ", termination_status(model))
println("最適解 a = ", value(a))
println("最適解 b = ", value(b))
println("最適解 c = ", value(c))
println("最小値 = ", objective_value(model))

プログラムは簡単なので説明は割愛させていただきます。実行結果を示します。

$ julia sample01.jl
Min 4 a + 2 b + 5 c
Subject to
 3 a + b + 2 c ≥ 15
 a + 2 b + 4 c ≥ 10
 a ≥ 0
 b ≥ 0
 c ≥ 0

結果: OPTIMAL
最適解 a = 4.0
最適解 b = 3.0
最適解 c = 0.0
最小値 = 22.0

食品 a と b だけを購入したほうが安くなりました。もし、食品 b の単価が 2 から 3 に値上がりしたとすると、結果は次のようになります。

Min 4 a + 3 b + 5 c
Subject to
 3 a + b + 2 c ≥ 15
 a + 2 b + 4 c ≥ 10
 a ≥ 0
 b ≥ 0
 c ≥ 0

結果: OPTIMAL
最適解 a = 4.0
最適解 b = 0.0
最適解 c = 1.5
最小値 = 23.5

この場合は食品 a と c だけを購入したほうが安くなります。ですが、経費は 22.0 から 23.5 に増えてしまいました。

[生産計画問題]

4 つの原料 (w, x, y, z) を使って 4 つの製品 (a, b, c, d) を生産します。製品 1 kg を作るのに必要な原材料、一日に使用できる原料の量 (kg)、製品の利益 (万円) が下表のように定義されています。一日の利益が最大となるように、各製品の生産量を決めてください。

表 : 原料の使用料と利益
	原料w	原料x	原料y	原料z	利益
製品a	2	1	0	0	5
製品b	0	2	1	0	3
製品c	0	0	1	2	2
製品d	1	0	0	2	4
使用量	4	8	6	10

各製品の生産量を変数 a, b, c, d で表すことにすると、目的関数は次式のようになります。

目的関数:
    g = 5 * a + 3 * b + 2 * c + 4 * d

この式が最大となる各変数の値を求めます。一日に使用できる原料は上限値があるので、これが制約条件となります。

制約条件:
    2 * a + 1 * d <= 4,   原料 w の制約
    1 * a + 2 * b <= 8,   原料 x の制約
    1 * b + 1 * c <= 6,   原料 y の制約
    2 * c + 2 * d <= 10,  原料 z の制約
    a, b, c, d >= 0,      非負条件

あとはこの式を JuMP でプログラムするだけです。

リスト : 生産計画問題 (sample02.jl)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 変数の定義
@variable(model, a >= 0)
@variable(model, b >= 0)
@variable(model, c >= 0)
@variable(model, d >= 0)

# 目的関数
@objective(model, Max, 5a + 3b + 2c + 4d)

# 制約条件
@constraint(model, 2a + d <= 4)
@constraint(model, a + 2b <= 8)
@constraint(model, b + c <= 6)
@constraint(model, 2c + 2d <= 10)

# 実行
println(model)
set_silent(model)
optimize!(model)

# 結果の表示
println("結果: ", termination_status(model))
println("最適解 a = ", value(a))
println("最適解 b = ", value(b))
println("最適解 c = ", value(c))
println("最適解 d = ", value(d))
println("最小値 = ", objective_value(model))

実行結果を示します。

$ julia sample02.jl
Max 5 a + 3 b + 2 c + 4 d
Subject to
 2 a + d ≤ 4
 a + 2 b ≤ 8
 b + c ≤ 6
 2 c + 2 d ≤ 10
 a ≥ 0
 b ≥ 0
 c ≥ 0
 d ≥ 0

結果: OPTIMAL
最適解 a = 0.0
最適解 b = 4.0
最適解 c = 1.0
最適解 d = 4.0
最小値 = 30.0

製品 a の生産量が 0 になってしまいました。そこで、変数の制約を a, b, c, d >= 1 に変更して実行してみましょう。

Max 5 a + 3 b + 2 c + 4 d
Subject to
 2 a + d ≤ 4
 a + 2 b ≤ 8
 b + c ≤ 6
 2 c + 2 d ≤ 10
 a ≥ 1
 b ≥ 1
 c ≥ 1
 d ≥ 1

結果: OPTIMAL
最適解 a = 1.0
最適解 b = 3.5
最適解 c = 2.5
最適解 d = 2.0
最小値 = 28.5

すべての製品を生産しますが、利益 g は 30.0 から 28.5 に低下しました。製品 a の利益が高いと違った結果になるでしょう。興味のある方は試してみてください。

[輸送問題]

工場 (x, y) から商品を店 (a, b, c) に配送します。供給量、需要量、輸送コストが下表で与えられているとき、総輸送コストが最小となる配送の仕方を求めてください。

表 : 輸送コスト
	店 a	店 b	店 c	供給量
工場 x	10	6	16	8
工場 y	8	8	4	16
需要量	12	4	8

工場 x から店に配送する商品の量を変数 xa, xb, xc とし、工場 y から店に配送する商品の量を変数 ya, yb, yc とします。目的関数は次式のようになります。

目的関数:
    z = 10 * xa + 6 * xb + 16 * xc + 8 * ya + 8 * yb + 4 * yc

この式が最小となる各変数の値を求めます。工場の供給量と店の需要には制限があるので、これが制約条件となります。

制約条件:
    xa + xb + xc == 8,  工場 x の供給量
    ya + yb + yc == 16, 工場 y の供給量
    xa + ya == 12,      店 a の需要
    xb + yb == 4,       店 b の需要
    xc + yc == 8,       店 c の需要
    xa, xb, xc >= 0,    非負条件
    ya, yb, yc >= 0

あとはこの式を JuMP でプログラムするだけです。

リスト : 輸送問題

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 変数の定義
@variable(model, xa >= 0)
@variable(model, xb >= 0)
@variable(model, xc >= 0)
@variable(model, ya >= 0)
@variable(model, yb >= 0)
@variable(model, yc >= 0)

# 目的関数
@objective(model, Min, 10xa + 6xb + 16xc + 8ya + 8yb + 4yc)

# 制約条件
@constraint(model, xa + xb + xc == 8)
@constraint(model, ya + yb + yc == 16)
@constraint(model, xa + ya == 12)
@constraint(model, xb + yb == 4)
@constraint(model, xc + yc == 8)

# 実行
println(model)
set_silent(model)
optimize!(model)

# 結果の表示
println("結果: ", termination_status(model))
println("最適解 xa = ", value(xa))
println("最適解 xb = ", value(xb))
println("最適解 xc = ", value(xc))
println("最適解 ya = ", value(ya))
println("最適解 yb = ", value(yb))
println("最適解 yc = ", value(yc))
println("最小値 = ", objective_value(model))

実行結果を示します。

$ julia sample03.jl
Min 10 xa + 6 xb + 16 xc + 8 ya + 8 yb + 4 yc
Subject to
 xa + xb + xc = 8
 ya + yb + yc = 16
 xa + ya = 12
 xb + yb = 4
 xc + yc = 8
 xa ≥ 0
 xb ≥ 0
 xc ≥ 0
 ya ≥ 0
 yb ≥ 0
 yc ≥ 0

結果: OPTIMAL
最適解 xa = 4.0
最適解 xb = 4.0
最適解 xc = 0.0
最適解 ya = 8.0
最適解 yb = 0.0
最適解 yc = 8.0
最小値 = 160.0

この例では、工場 x, y ともに半分ずつの量を店に配送するのが一番安くなりました。輸送コストが変化すると、当然ですが違った結果になります。興味のある方はいろいろ試してみてください。

●鶴亀算

JuMP を使って連立方程式を解くこともできます。簡単な例として「鶴亀算」を解いてみましょう。

鶴と亀、合わせて 100 匹いる。足の合計が 272 本のとき、鶴と亀はそれぞれ何匹ずついるか。
鶴と亀とトンボが合わせて 10 匹いる。足の合計が 38 本で羽の合計が 14 枚であるとき、鶴と亀とトンボはそれぞれ何匹ずついるか。(トンボの足は 6 本で羽は 4 枚)
鶏と犬とタコ、合わせて 24 匹が台所にいる。足の合計が 102 本のとき、鶏、犬、タコはそれぞれ何匹ずついるか。

問題 1 は次の連立方程式を解けば求めることができます。

x + y = 100
2x + 4y = 272

リスト : 鶴亀算 問題１ (turukame.jl)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 変数の定義
@variable(model, x, Int)
@variable(model, y, Int)

# 目的関数 (省略可)
@objective(model, Min, x)

# 制約
@constraint(model, x + y == 100)
@constraint(model, 2x + 4y == 272)

# 実行
println(model)
set_silent(model)
optimize!(model)

# 結果の表示
println("結果: ", termination_status(model))
println("目的関数の値: ", objective_value(model))
println("最適解 x = ", value(x))
println("最適解 y = ", value(y))

目的関数は x の最小値を求めていますが、解は 1 通りしかないので、x の最大値でも y の最小値でもかまいません。実行結果を示します。

$ julia turukame.jl
Min x
Subject to
 x + y = 100
 2 x + 4 y = 272
 x integer
 y integer

結果: OPTIMAL
目的関数の値: 64.0
最適解 x = 64.0
最適解 y = 36.0

なお、JuMP の場合、目的関数を省略することができます。その場合、ソルバーは「目的関数の値を最大化・最小化する」のではなく、「すべての制約条件を満たす解 (実行可能解) を 1 つ見つける」という動作になります。解が複数存在する場合、ソルバーは最初に見つけた解を返します。どの解が最初に見つかるかは、ソルバーのアルゴリズムに依存します。このとき objective_value(model) の値は 0.0 になります。実際に試してみましょう。

$ julia turukame.jl
Feasibility
Subject to
 x + y = 100
 2 x + 4 y = 272
 x integer
 y integer

結果: OPTIMAL
目的関数の値: 0.0
最適解 x = 64.0
最適解 y = 36.0

問題 2 は次の連立方程式を解けば求めることができます。

 x +  y +  z = 10
2x + 4y + 6z = 38
2x +      4z = 14

プログラムと実行結果を示します。

リスト : 鶴亀算 問題２ (turukame2.jl)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 変数の定義
@variable(model, x, Int)
@variable(model, y, Int)
@variable(model, z, Int)

# 制約
@constraint(model, x + y + z == 10)
@constraint(model, 2x + 4y + 6z == 38)
@constraint(model, 2x + 4z == 14)

# 実行
println(model)
set_silent(model)
optimize!(model)

# 結果の表示
println("結果: ", termination_status(model))
println("最適解 x = ", value(x))
println("最適解 y = ", value(y))
println("最適解 z = ", value(z))

$ julia turukame2.jl
Feasibility
Subject to
 x + y + z = 10
 2 x + 4 y + 6 z = 38
 2 x + 4 z = 14
 x integer
 y integer
 z integer

結果: OPTIMAL
最適解 x = 3.0
最適解 y = 5.0
最適解 z = 2.0

問題 3 は答えを 1 つに決めることはできませんが、適切な目的関数と制約条件を設定すれば、それを満たす解を求めることができます。

プログラムは次のようになります。

リスト : 鶴亀算 問題３ (turukame3.jl)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 変数の定義
@variable(model, x >= 1, Int)
@variable(model, y >= 1, Int)
@variable(model, z >= 1, Int)

# 目的関数
@objective(model, Max, z)

# 制約
@constraint(model, x + y + z == 24)
@constraint(model, 2x + 4y + 8z == 102)

# 実行
println(model)
set_silent(model)
optimize!(model)

# 結果の表示
println("結果: ", termination_status(model))
println("目的関数の値: ", objective_value(model))
println("最適解 x = ", value(x))
println("最適解 y = ", value(y))
println("最適解 z = ", value(z))

問題３の場合、変数の値を正、目的関数を z とし、z が最大となる解を求めます。結果は次のようになりました。

$ julia turukame3.jl
Max z
Subject to
 x + y + z = 24
 2 x + 4 y + 8 z = 102
 x ≥ 1
 y ≥ 1
 z ≥ 1
 x integer
 y integer
 z integer

結果: OPTIMAL
目的関数の値: 8.0
最適解 x = 13.0
最適解 y = 3.0
最適解 z = 8.0

●部分和問題

部分和問題 (subset-sum problem) は、要素が数値の集合 S において、要素の総和が M となる部分集合があるか判定する問題です。たとえば、集合 {2, 3, 5, 8} の場合、総和が 10 となる部分集合は {2, 3, 5} と {2, 8} がありますが、14 となる部分集合はありません。部分集合の総数は、要素数を n とすると 2ⁿ 個になるので、n が大きくなるとナイーブな方法では時間がかかってしまいます。

実際には、動的計画法を使って現実的な時間で部分和問題を解くことができるといわれています。最適化ソルバーでも高速に解くことができるか試してみましょう。

部分和問題は集合 X の要素を x_i とし、その係数を a_iすると、次の等式を満たすか判定する問題になります。

\( \displaystyle \sum_{i=1}^{n} a_i x_i = M, \quad (a_i = 0 \ \mathrm{or} \ 1) \)

係数 a_i の値が 0 か 1 かを決める問題になります。このような問題を「0-1 整数計画問題」といいます。JuMP を使うと、とても簡単にプログラムすることができます。次のリストを見てください。

リスト : 部分和問題 (subsetsum.jl)

using JuMP, HiGHS

nums = [  1,   2,   3,   5,   8,   13,   21,   34,   55,    89,
        144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946]

function solver(xs, n)
  size = length(xs)
  model = Model(HiGHS.Optimizer)

  # 変数
  @variable(model, vs[1:size], Bin)
  # 式の定義
  @expression(model, expr, sum(xs[i] * vs[i] for i = 1 : size))
  # 制約
  @constraint(model, expr == n)
  # 実行
  println(model)
  set_silent(model)
  optimize!(model)

  println("結果 ", termination_status(model))
  println("値 ", objective_value(model))
  println([xs[i] for i = 1 : size if value(vs[i]) == 1])
end

関数 solver の引数 xs が数の集合で、引数 n が求める部分集合の総和を表します。JuMP の変数は xs の要素と同じ数だけ生成して、配列 vs に格納します。つまり、xs[i] に対応する変数が vs[i] になります。JuMP の場合、マクロ @variable で要素が変数の配列を簡単に定義することができます。このとき、データ型に Bin を指定して変数の値を 0 or 1 にします。

次に、@expression で内積を求める式を定義します。これは sum() にジェネレータ式 xs[i] * vs[i] for i = 1 : size を渡すだけです。この式に名前 expr を付けます。すると、制約 @constraint の条件式は expr == n と表すことができます。解があれば結果は Optimal になり、解がない場合は Infeasible になります。

リスト nums はフィボナッチ数列になっています。要素の総和を M とすると、1 から M までの整数は、要素を組み合わせて必ず作ることができます。これはフィボナッチ数列の面白い特徴です。これをテストに使います。それでは実際に試してみましょう。

julia> include("subsetsum.jl")
... 略 ...

julia> solver(nums, sum(nums) - 1)
Feasibility
Min vs[1] + 2 vs[2] + 3 vs[3] + 5 vs[4] + 8 vs[5] + 13 vs[6] + 21 vs[7] + 34 vs[8] + 55 vs[9] + 89 vs[10] + 
144 vs[11] + 233 vs[12] + 377 vs[13] + 610 vs[14] + 987 vs[15] + 1597 vs[16] + 2584 vs[17] + 4181 vs[18] + 
6765 vs[19] + 10946 vs[20]
Subject to
 vs[1] + 2 vs[2] + 3 vs[3] + 5 vs[4] + 8 vs[5] + 13 vs[6] + 21 vs[7] + 34 vs[8] + 55 vs[9] + 89 vs[10] + 
144 vs[11] + 233 vs[12] + 377 vs[13] + 610 vs[14] + 987 vs[15] + 1597 vs[16] + 2584 vs[17] + 4181 vs[18] + 
6765 vs[19] + 10946 vs[20] = 28654
 vs[1] binary
 vs[2] binary
 vs[3] binary
 vs[4] binary
 vs[5] binary
 vs[6] binary
 vs[7] binary
 vs[8] binary
 vs[9] binary
 vs[10] binary
 vs[11] binary
 vs[12] binary
 vs[13] binary
 vs[14] binary
 vs[15] binary
 vs[16] binary
 vs[17] binary
 vs[18] binary
 vs[19] binary
 vs[20] binary

結果 OPTIMAL
値 0.0
[2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946]

julia> solver(nums, sum(nums) + 1)
... 略 ...

結果 INFEASIBLE
値 0.0
Int64[]

変数の個数が少ないので一瞬で答えを求めることができます。ところで、解が無くてもそれに近い値がほしい場合があるでしょう。この場合、制約条件を expr <= n とし、目的関数で最大値を探索すれば、n 以下でもっとも近い値を求めることができます。プログラムは次のようになります。

リスト : 部分和問題 (2)

function solver1(xs, n)
  size = length(xs)
  model = Model(HiGHS.Optimizer)

  # 変数
  @variable(model, vs[1:size], Bin)
  # 式の定義
  @expression(model, expr, sum(xs[i] * vs[i] for i = 1 : size))
  # 目的関数
  @objective(model, Max, expr)
  # 制約
  @constraint(model, expr <= n)
  # 実行
  println(model)
  set_silent(model)
  optimize!(model)

  println("結果 ", termination_status(model))
  println("値 ", objective_value(model))
  println([xs[i] for i = 1 : size if value(vs[i]) == 1])
end

それでは実際に試してみましょう。

julia> solver1(nums, sum(nums) - 1)
... 略 ...

結果 OPTIMAL
値 28654.0
[2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946]

julia> solver1(nums, 30000)
... 略 ...

結果 OPTIMAL
値 28655.0
[1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946]

引数 n と求めた解が等しければ、n と等しくなる部分集合を求めることができました。そうでなければ、n と等しい部分集合はありませんが、それに近い値と部分集合を得ることができます。

●ナップザック問題

次は「ナップザック問題」を取り上げます。ナップザック (knapsack) とは辞書を引いてみると、ランドセルのような背中にせおう四角形の袋や箱のことを意味します。ここでは物を入れる袋と簡単に考えてください。

ここで、ナップザックの中に品物を詰め込むことを考えてみます。一つのナップザックと複数の品物が与えられたとき、袋に詰めた品物の合計金額が最大になるような選び方を求めることが「ナップザック問題」です。ナップザック問題にはバリエーションがあって、同じ品物をいくつも入れて良い場合と、一つしか入れてはいけない場合があります。後者の場合を「0-1 ナップザック問題」といいます。

ナップザック問題は、部分和問題と同様に NP 問題になります。これは厳密に解を求めようとすると、全ての場合について総当たりで調べるしか方法がなく、データ数が多くなると時間がかかるため、現実的な時間では解答を出すことができないというものです。品物の詰め方が難問の一つ、といわれてもピンとこないと思いますが、ナップザック問題は品物の種類が増えるにしたがって、その組み合わせ方が爆発的に増えるのです。

ところが、幸いなことにナップザック問題は実用的には解決済みの問題と考えられています。とくに有名なのが「動的計画法」を用いた解法です。ナップザックと品物の大きさを整数値に限定すれば、動的計画法を用いることで厳密解を求めることが可能です。興味のある方は拙作のページ「Algorithms with Python: 動的計画法」をお読みください。

まず最初に、0-1 ナップザック問題を解いてみましょう。0-1 ナップザック問題を数式で表すと次のようになります。

W はナップザックの大きさで、V の最大値を求める

\(\begin{array}{l} \displaystyle \sum_{i=1}^{n} a_i w_i \leq W \\ \displaystyle \sum_{i=1}^{n} a_i p_i = V \end{array}\)

\(w_i\) : 重さ
\(p_i\) : 金額
\(a_i\) : 0 or 1
\(n\) : 品物の個数

基本的には、部分和問題と同じようなプログラムになりますが、ナップザックに入る範囲で、最大の金額になるような入れ方を求める必要があります。それでは実際に簡単な問題を解いてみましょう。

[問題１]

下表に示す品物をサイズ 15 のナップザックに入れるとき、金額が最大となる入れ方を求めてください。

品物	金額	サイズ
A	4	3
B	5	4
C	6	5
D	8	7
E	10	9

出典 : Coprisによる制約プログラミング入門, (田村直之さん)

プログラムは次のようになります。

リスト : 問題１の解法 (knap1.jl)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 品物
goods = ['a', 'b', 'c', 'd', 'e']
# 金額
price = [4, 5, 6, 8, 10]
# 重さ
weight = [3, 4, 5, 7, 9]
# 大きさ
size = length(goods)

# 変数の定義
@variable(model, xs[1:size], Bin)
# 目的関数
@objective(model, Max, sum(price[i] * xs[i] for i = 1 : size))
# 制約条件
@constraint(model, sum(weight[i] * xs[i] for i = 1 : size) <= 15)

# 実行
println(model)
set_silent(model)
optimize!(model)
# 結果
println("結果 ", termination_status(model))
println("値 ", objective_value(model))
println(value.(xs))

$ julia knap1.jl
Max 4 xs[1] + 5 xs[2] + 6 xs[3] + 8 xs[4] + 10 xs[5]
Subject to
 3 xs[1] + 4 xs[2] + 5 xs[3] + 7 xs[4] + 9 xs[5] ≤ 15
 xs[1] binary
 xs[2] binary
 xs[3] binary
 xs[4] binary
 xs[5] binary

結果 OPTIMAL
値 18.0
[1.0, 0.0, 1.0, 1.0, 0.0]

金額の最大値は 18 で、選択した品物は A, C, D の 3 つです。

次は同じ品物をいくつ選んでもよい問題を解いてみましょう。この場合、変数は Bin ではなく Int になります。

[問題２]

下表に示す品物をサイズ 10 のナップザックに入れるとき、金額が最大となる入れ方を求めてください。なお、同じ品物を何個入れてもかまいません。

品物	金額	サイズ
A	6	4
B	4	3
C	1	1

リスト : 問題２の解法 (knap2.jl)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 品物
goods = ['a', 'b', 'c']
# 金額
price = [6, 4, 1]
# 重さ
weight = [4, 3, 1]
# 大きさ
size = length(goods)

# 変数の定義
@variable(model, xs[1:size] >= 0, Int)
# 目的関数
@objective(model, Max, sum(price[i] * xs[i] for i = 1 : size))
# 制約条件
@constraint(model, sum(weight[i] * xs[i] for i = 1 : size) <= 10)

# 実行
println(model)
set_silent(model)
optimize!(model)
# 結果
println("結果 ", termination_status(model))
println("値 ", objective_value(model))
println(value.(xs))

$ julia knap2.jl
Max 6 xs[1] + 4 xs[2] + xs[3]
Subject to
 4 xs[1] + 3 xs[2] + xs[3] ≤ 10
 xs[1] ≥ 0
 xs[2] ≥ 0
 xs[3] ≥ 0
 xs[1] integer
 xs[2] integer
 xs[3] integer

結果 OPTIMAL
値 14.0
[2.0, 0.0, 2.0]

金額の最大値は 14 で、A を 2 個、C を 2 個選びます。なお、A を 1 個、B を 2 個選んでも金額は 14 になります。

もう一つ、簡単な例を示しましょう。

[問題３]

下表に示す品物をサイズ w のナップザックに入れるとき、金額が最大となる入れ方を求めてください。なお、同じ品物を何個入れてもかまいません。

品物	金額	サイズ
A	91	3
B	120	4
C	610	20
D	930	30

リスト : 問題３の解法 (knap3.jl)

using JuMP, HiGHS

function solver(w)
  model = Model(HiGHS.Optimizer)

  # 品物
  goods = ['a', 'b', 'c', 'd']
  # 金額
  price = [91, 120, 610, 930]
  # 重さ
  weight = [3, 4, 20, 30]
  # 大きさ
  size = length(goods)

  # 変数の定義
  @variable(model, xs[1:size] >= 0, Int)
  # 目的関数
  @objective(model, Max, sum(price[i] * xs[i] for i = 1 : size))
  # 制約条件
  @constraint(model, sum(weight[i] * xs[i] for i = 1 : size) <= w)

  # 実行
  println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  println("値 ", objective_value(model))
  println(value.(xs))
end

julia> include("knap3.jl")
solver (generic function with 1 method)

julia> solver(500)
Max 91 xs[1] + 120 xs[2] + 610 xs[3] + 930 xs[4]
Subject to
 3 xs[1] + 4 xs[2] + 20 xs[3] + 30 xs[4] ≤ 500
 xs[1] ≥ 0
 xs[2] ≥ 0
 xs[3] ≥ 0
 xs[4] ≥ 0
 xs[1] integer
 xs[2] integer
 xs[3] integer
 xs[4] integer

結果 OPTIMAL
値 15490.0
[0.0, 0.0, 1.0, 16.0]

julia> solver(1000)
Max 91 xs[1] + 120 xs[2] + 610 xs[3] + 930 xs[4]
Subject to
 3 xs[1] + 4 xs[2] + 20 xs[3] + 30 xs[4] ≤ 1000
 xs[1] ≥ 0
 xs[2] ≥ 0
 xs[3] ≥ 0
 xs[4] ≥ 0
 xs[1] integer
 xs[2] integer
 xs[3] integer
 xs[4] integer

結果 OPTIMAL
値 30992.0
[2.0, 1.0, 0.0, 33.0]

サイズ w が大きな値でも高速に解を求めることができました。興味のある方はいろいろ試してみてください。

Julia Programming