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

●施設配置問題

「施設配置問題 (facility location problem)」は、工場や倉庫など重要な施設の配置を決定する問題です。施設配置問題には様々な種類の問題がありますが、今回は基本的な問題である p-median 問題と p-center 問題を取り上げます。どちらの問題も、施設の候補地の集合と、その利用者の集合が与えられていて、ある基準を満たすよう候補地の集合から p 個の施設を選択します。p-median 問題は、利用者の移動距離の総和が最小となる配置を求めます。p-center 問題は、利用者の最大の移動距離が最小となる配置を求めます。

●p-median 問題

簡単な例題として、施設の候補地の集合 F = {f1, f2, f3} と利用者の集合 D = {d1, d2, d3} において、施設 f を 2 つ選ぶ場合を考えてみましょう。下表に施設と利用者の距離と利用者の人数を示します。

    表 : 人数と距離

        : f1 : f2 : f3 
--------+----+----+-----
d1 (2)  : 11 :  3 :  6
d2 (10) :  1 :  7 :  9
d3 (4)  :  8 :  5 : 12

出典 : ロジスティクス工学 05 施設配置 (参考 URL 2)

施設 \(f_{i}\) の選択を表す変数 (Bin) を \(F_{i}\) とし、利用者 \(d_{i}\) が施設 \(f_{j}\) を使うことを表す変数 (Bin) を \(X_{ij}\) とします。\(d_{i}\) の人数を \(W_{i}, d_{i}\) と \(f_{j}\) の距離を \(C_{ij}\) とすると、目的関数と制約条件は次のようになります。

目的関数 (最小化):

\( \displaystyle \sum_i \sum_j W_{i} C_{ij} X_{ij} \)

制約条件:

1. 配置する施設は二つ
\( \displaystyle \sum_i F_{i} = 2 \)

2. 利用者が使用する施設は一つだけ
\( \displaystyle \sum_j X_{ij} = 1, \quad (i = 1, 2, 3) \)

3. \(F_{j}\) が 0 ならば \(X_{ij}\) も 0 になる
\( X_{ij} \leq F_{j}, \quad (i = 1, 2, 3, j = 1, 2, 3) \)

最後の制約条件は、施設 \(F_{j}\) を配置しないのであれば、その施設は利用できないことをことを表します。あとはこれを JuMP でプログラムするだけです。プログラムリストと実行結果を示します。

リスト : 施設配置問題 (p-median)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# コスト (距離)
#      f1  f2  f3
cs = [ 11   3   6;    # d1
        1   7   9;    # d2
        8   5  12 ]   # d3

# 人数
ws = [2, 10, 4]

# 変数
# fs[i]    : 施設 fi を配置するか否か
# xs[i, j] : di が fi を利用するか否か
@variable(model, fs[1:3], Bin)
@variable(model, xs[1:3, 1:3], Bin)

# 目的関数
@objective(model, Min, sum(ws[i] * sum(cs[i, j] * xs[i, j] for j = 1:3) for i = 1:3))

# 制約条件
@constraint(model, sum(fs) == 2)
for i = 1:3
  @constraint(model, sum(xs[i, j] for j = 1:3) == 1)
end
for i = 1:3
  for j = 1:3
    @constraint(model, xs[i, j] <= fs[j])
  end
end

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

結果 OPTIMAL
[1.0, 1.0, -0.0]
[0.0 1.0 0.0;
 1.0 0.0 0.0;
 0.0 1.0 0.0]
値 36.0

施設を配置するのは f1 と f2 で、d1 と d3 が f2 を利用し、d2 が f1 を利用します。このときの費用は 36 になりました。

次は、顧客 (1 つ) を平面上にランダムに n 個配置して、その中から施設となるものを p 個選んでみましょう。この場合、施設と顧客の距離の総和が最小となる p 個のグループに顧客を分けることになります。プログラムは次のようになります。

リスト : p-median 問題 (乱数での配置)

function solver_median(ps, p)
  size = length(ps)
  cs = make_matrix(ps)   # 隣接行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, xs[1:size, 1:size], Bin)
  @variable(model, fs[1:size], Bin)

  # 目的関数
  @objective(model, Min, sum(cs .* xs))

  # 制約条件
  @constraint(model, sum(fs) == p)
  for i = 1:size
    @constraint(model, sum(xs[i, j] for j = 1:size) == 1)
  end
  for i = 1:size
    @constraint(model, xs[i, i] >= fs[i])
    for j = 1:size
      @constraint(model, xs[i, j] <= fs[j])
    end
  end
  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  println(value.(fs))
  println(value.(xs))
  println("値 ", objective_value(model))
end

関数 make_matirx() は拙作のページ「巡回セールスマン問題」のプログラムを流用しています。制約条件では、次の制約式を追加しています。

\(X_{ii} \geq F_{i}\)

施設 \(F_{i}\) を選んだ場合、i 番目の顧客は自分自身を利用することになるので、\(X_{ii}\) の値は 1 になります。上記の制約条件を追加することで、fs[i] が 1 ならば xs[i, i] の値を 1 にすることができます。なお、この制約条件がなくてもプログラムは正常に動作しますが、制約を追加すると実行時間がほんのちょっとですが速くなるようです。

それでは実行してみましょう。最初は 12 個の中から施設を 3 つ選んでみましょう。

julia> solver_median(make_data(12), 3)
結果 OPTIMAL
[0.0,  0.0, 0.0, 1.0, 0.0, -0.0, 1.0, -0.0, 0.0, -0.0, 0.0, 1.0]
[ 0.0 -0.0  0.0  1.0 -0.0  0.0  0.0  0.0  0.0 -0.0  0.0  0.0;
 -0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0 -0.0  0.0  0.0  0.0;
  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -0.0  0.0  0.0  0.0  1.0;
 -0.0 -0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
  0.0  0.0  0.0  1.0  0.0 -0.0  0.0  0.0  0.0 -0.0  0.0  0.0;
  0.0  0.0  0.0  1.0 -0.0  0.0  0.0  0.0  0.0 -0.0  0.0  0.0;
  0.0  0.0  0.0  0.0 -0.0  0.0  1.0  0.0  0.0  0.0 -0.0  0.0;
  0.0  0.0 -0.0  0.0  0.0  0.0  0.0  0.0 -0.0  0.0  0.0  1.0;
 -0.0 -0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
  0.0  0.0  0.0  1.0 -0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0;
  0.0  0.0 -0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0;
  0.0  0.0 -0.0  0.0  0.0  0.0  0.0 -0.0  0.0  0.0  0.0  1.0]
値 1077.8907165041312

12 個の中から 4, 7, 12 番目を施設として選びます。4 番目の施設を利用するのが 7 個、7 番目の施設を利用するのが 2 個、12 番目の施設を利用するのが 3 個になり、最小値は 1077.9 となりました。どうやら正常に動作しているようです。

次は個数を増やして、実行時間を計測してみましょう。

julia> @time solver_median(make_data(60), 3)
結果 OPTIMAL
値 6374.914224617162
  0.214874 seconds (198.87 k allocations: 13.189 MiB)

julia> @time solver_median(make_data(80), 4)
結果 OPTIMAL
値 6841.956722259684
  0.380135 seconds (348.93 k allocations: 20.590 MiB)

julia> @time solver_median(make_data(100), 5)
結果 OPTIMAL
値 7817.760050808525
  0.661048 seconds (541.42 k allocations: 31.232 MiB, 3.75% gc time)

julia> @time solver_median(make_data(120), 6)
結果 OPTIMAL
値 8516.605447642038
  1.115110 seconds (776.36 k allocations: 52.238 MiB, 1.40% gc time)

julia> @time solver_median(make_data(140), 7)
結果 OPTIMAL
値 9090.820206226033
  2.384160 seconds (1.05 M allocations: 66.875 MiB, 1.44% gc time)

実行環境 : Julia 1.12.6, Ubunts 24.04 LTS (WSL2, Windows 11), intel CORE i5-1235U 1.30GHz
ソルバー : JuMP v1.30.1, HiGHS v1.24.0

顧客の個数が 100 程度であれば、制約式を強化しなくても高速に解くことができるようです。ところで、施設の個数の割合を変更すると、実行時間が変化するかもしれません。興味のある方はいろいろ試してみてください。

●p-center 問題

次は p-center 問題のプログラムを作りましょう。まずは最初に p-median と同じ例題を取り上げます。p-center 問題の場合、利用者の人数は関係なく、施設と利用者の距離が問題になります。

    表 : 距離

   : f1 : f2 : f3 
---+----+----+----
d1 : 11 :  3 :  6
d2 :  1 :  7 :  9
d3 :  8 :  5 : 12

出典 : 参考 URL 2『ロジスティクス工学 05 施設配置』

施設 \(f_{i}\) の選択を表す変数 (Bin) を \(F_{i}\) とし、利用者 \(d_{i}\) が施設 \(f_{j}\) を使うことを表す変数 (Bin) を \(X_{ij}\) とします。\(d_{i}\) と \(f_{j}\) の距離を \(C_{ij}\) とすると、目的関数と制約条件は次のようになります。

目的関数 (最小化): \(z\)

制約条件:

1. i と j の距離を z 以下にする
\( \displaystyle \sum_j C_{ij} X_{ij} \leq z \)

2. 配置する施設は二つ
\( \displaystyle \sum_i F_{i} = 2 \)

3. 利用者が使用する施設は一つだけ
\( \displaystyle \sum_j X_{ij} = 1, \quad (i = 0, 1, 2) \)

4. \(F_{j}\) が 0 ならば \(X_{ij}\) も 0 になる
\( X_{ij} \leq F_{j}, \quad (i = 0, 1, 2, j = 0, 1, 2) \)

ポイントは最大の移動距離を変数 z で表すところです。変数 z を目的関数にすることで、z を最小化することができます。そして、制約条件に最初の式を追加します。i -> j までの距離が z 以下である条件を追加することで、z が最小となるように \(X_{ij}\) が選択されます。

あとはこれを JuMP でプログラムするだけです。プログラムリストと実行結果を示します。

リスト : 施設配置問題 (p-center)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# コスト (距離)
#      f1  f2  f3
cs = [ 11   3   6;    # d1
        1   7   9;    # d2
        8   5  12 ]   # d3

# 変数
# fs[i]    : 施設 fi を配置するか否か
# xs[i, j] : di が fi を利用するか否か
#
@variable(model, fs[1:3], Bin)
@variable(model, xs[1:3, 1:3], Bin)
@variable(model, 0 <= z)

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

# 制約条件
@constraint(model, sum(fs) == 2)
for i = 1:3
  @constraint(model, sum(cs[i, j] * xs[i, j] for j = 1:3) <= z)
  @constraint(model, sum(xs[i, j] for j = 1:3) == 1)
end
for i = 1:3
  for j = 1:3
    @constraint(model, xs[i, j] <= fs[j])
  end
end

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

結果 OPTIMAL
[ 1.0, 1.0, -0.0]
[-0.0  1.0  -0.0; 
  1.0 -0.0  -0.0;
 -0.0  1.0  -0.0]
値 5.0

施設を配置するのは f1 と f2 で、d1 と d3 が f2 を利用し、d2 が f1 を利用します。このときの最大距離 z は d3 -> f2 の 5 になりました。d1 -> f2 は 3 で、d2 -> f1 は 1 なので、z よりも小さく条件を満たしていることがわかります。

次は、顧客 (1 つ) を平面上にランダムに n 個配置して、その中から施設となるものを p 個選んでみましょう。この場合、施設と顧客の最大距離が最小となる p 個のグループに顧客を分けることになります。プログラムは次のようになります。

リスト : p-center 問題 (乱数での配置)

function solver_center(ps, p)
  size = length(ps)
  cs = make_matrix(ps)   # 隣接行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, xs[1:size, 1:size], Bin)
  @variable(model, fs[1:size], Bin)
  @variable(model, 0 <= z)

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

  # 制約条件
  @constraint(model, sum(fs) == p)
  for i = 1:size
    @constraint(model, sum(cs[i, j] * xs[i, j] for j = 1:size) <= z)
  end
  for i = 1:size
    @constraint(model, sum(xs[i, j] for j = 1:size) == 1)
  end
  for i = 1:size
    @constraint(model, xs[i, i] >= fs[i])
    for j = 1:size
      @constraint(model, xs[i, j] <= fs[j])
    end
  end
  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  println(round.(Int, value.(fs)))
  println(round.(Int, value.(xs)))
  println("値 ", objective_value(model))
end

p-median の場合、制約条件 xs[i][i] >= fs[i] はなくても動作しましたが、p-center では必要になります。省略すると正常に動作しません。ご注意くださいませ。

それでは実行してみましょう。最初は 10 個の中から施設を 3 つ選んでみましょう。

julia> solver_center(make_data(10), 3)
結果 OPTIMAL
[1, 0, 0, 0, 0, 0, 1, 0, 0, 1]
[1 0 0 0 0 0 0 0 0 0;
 0 0 0 0 0 0 0 0 0 1;
 0 0 0 0 0 0 1 0 0 0;
 1 0 0 0 0 0 0 0 0 0;
 0 0 0 0 0 0 1 0 0 0;
 1 0 0 0 0 0 0 0 0 0;
 0 0 0 0 0 0 1 0 0 0;
 0 0 0 0 0 0 0 0 0 1;
 0 0 0 0 0 0 1 0 0 0;
 0 0 0 0 0 0 0 0 0 1]
値 154.18495387034275

10 個の中から 1, 7, 10 番目を施設として選びます。1 番目の施設を利用するのが 3 個、7 番目の施設を利用するのが 4 個、10 番目の施設を利用するのが 3 個になり、最小値は 154.2 となりました。どうやら正常に動作しているようです。

次は個数を増やして、実行時間を計測してみましょう。

julia> @time solver_center(make_data(20), 2)
結果 OPTIMAL
値 235.05318547086318
  0.464861 seconds (29.34 k allocations: 1.508 MiB)

julia> @time solver_center(make_data(40), 4)
結果 OPTIMAL
値 166.37307474468332
  1.612352 seconds (96.88 k allocations: 5.742 MiB)

julia> @time solver_center(make_data(60), 6)
結果 OPTIMAL
値 130.41855696180662
  4.781321 seconds (215.69 k allocations: 13.580 MiB)

julia> @time solver_center(make_data(80), 8)
結果 OPTIMAL
値 107.93516572461452
 11.528247 seconds (379.11 k allocations: 21.740 MiB)

julia> @time solver_center(make_data(100), 10)
結果 OPTIMAL
値 97.63708311906906
 14.890287 seconds (595.11 k allocations: 32.966 MiB)

p-center 問題は p-median 問題よりも実行時間がかかるようです。それでも、顧客が 100 個程度であれば 30 秒もかからずに解くことができました。興味のある方はいろいろ試してみてください。

●参考 URL

数理最適化入門（4）：施設配置の数理モデル [PDF], (田中健一さん)
ロジスティクス工学 05 施設配置 [PDF], (宮本裕一郎さん)

●プログラムリスト

#
# shisetu.jl : 施設配置問題 (p-median と p-center)
#
#              Copyright (C) 2026 Makoto Hiroi
#
using JuMP, HiGHS

# 乱数でデータを生成
function make_data(n)
  [(rand(10:490), rand(10:490)) for _ = 1:n]
end

# 距離を求める
function distance(p1, p2)
  dx = p1[1] - p2[1]
  dy = p1[2] - p2[2]
  sqrt(dx * dx + dy * dy)
end

# 隣接行列の生成
function make_matrix(ps)
  size = length(ps)
  xs = zeros(size, size)
  for i = 1:size
    for j = 1:size
      xs[i, j] = distance(ps[i], ps[j])
    end
  end
  xs
end

# 解を求める
function solver_median(ps, p)
  size = length(ps)
  cs = make_matrix(ps)   # 隣接行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, xs[1:size, 1:size], Bin)
  @variable(model, fs[1:size], Bin)

  # 目的関数
  @objective(model, Min, sum(cs .* xs))

  # 制約条件
  @constraint(model, sum(fs) == p)
  for i = 1:size
    @constraint(model, sum(xs[i, j] for j = 1:size) == 1)
  end
  for i = 1:size
    @constraint(model, xs[i, i] >= fs[i])
    for j = 1:size
      @constraint(model, xs[i, j] <= fs[j])
    end
  end
  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  # println(value.(fs))
  # println(value.(xs))
  println("値 ", objective_value(model))
end

function solver_center(ps, p)
  size = length(ps)
  cs = make_matrix(ps)   # 隣接行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, xs[1:size, 1:size], Bin)
  @variable(model, fs[1:size], Bin)
  @variable(model, 0 <= z)

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

  # 制約条件
  @constraint(model, sum(fs) == p)
  for i = 1:size
    @constraint(model, sum(cs[i, j] * xs[i, j] for j = 1:size) <= z)
  end
  for i = 1:size
    @constraint(model, sum(xs[i, j] for j = 1:size) == 1)
  end
  for i = 1:size
    @constraint(model, xs[i, i] >= fs[i])
    for j = 1:size
      @constraint(model, xs[i, j] <= fs[j])
    end
  end
  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  # println(round.(Int, value.(fs)))
  # println(round.(Int, value.(xs)))
  println("値 ", objective_value(model))
end

●施設配置問題 (2)

施設配置問題の続きです。今回は集合被覆問題と最大被覆問題を取り上げます。

●集合被覆問題

施設配置問題でいう「集合被覆問題 (set cover problem)」は、施設の候補地の集合 F と利用者の集合 D が与えられたとき、すべての利用者が施設を利用でき、なおかつ施設の個数を最小にする配置を求める問題です。数学での定義は参考 URL 3『集合被覆問題 - Wikipedia』をお読みください。

簡単な例題として、利用者 \(D_{i}\) が利用できる施設の候補地 \(F_{j}\) が下表のように与えられているとき、施設の個数が最小となる配置を求めてみましょう。

    表 : 利用可能な施設

    : f1 : f2 : f3 : f4 : f5 
----+----+----+----+----+----
 d1 :  1 :  0 :  0 :  0 :  1
 d2 :  1 :  1 :  0 :  0 :  0
 d3 :  0 :  1 :  1 :  0 :  0
 d4 :  0 :  0 :  1 :  1 :  0
 d5 :  0 :  0 :  0 :  1 :  1

    (出典 : 参考 URL 2 『ロジスティクス工学 05 施設配置』)

施設 \(f_{j}\) の選択を表す変数 (Bin) を \(F_{j}\) とし、上表の要素を \(W_{ij}\) で表すことにします。目的関数と制約条件は次のようになります。

目的関数 (最小化):

\( \displaystyle \sum_j F_{j} \)

制約条件:

\( \displaystyle \sum_j W_{ij} F_{j} \geq 1 \quad (i= 1, 2, 3, 4, 5) \)

目的関数は簡単ですね。制約条件ですが、利用者 \(D_{i}\) が利用できる施設は W の i 行目に定義されているので、その要素と F の要素の内積を計算して、その値が 1 以上であれば \(D_{i}\) は施設を利用することができます。

あとはこれを JuMP でプログラムするだけです。プログラムリストと実行結果を示します。

リスト : 集合被覆問題 (set cover problem)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 利用可能な施設
#      f1 f2 f3 f4 f5
ws = [  1  0  0  0  1;    # d1
        1  1  0  0  0;    # d2
        0  1  1  0  0;    # d3
        0  0  1  1  0;    # d4
        0  0  0  1  1 ]   # d5

# 変数
@variable(model, fs[1:5], Bin)

# 目的関数
@objective(model, Min, sum(fs))

# 制約条件
for i = 1:5
  @constraint(model, sum(ws[i, j] * fs[j] for j = 1:5) >= 1)
end

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

結果 OPTIMAL
[-0.0, 1.0, 1.0, -0.0, 1.0]
値 3.0

配置する施設は f2, f4, f5 の三つになりました。

次は、利用者を平面上にランダムに n 個配置し、半径 r 内にある施設を利用できる条件で、集合被覆問題を解いてみましょう。施設は利用者の中から選びます。プログラムは次のようになります。

リスト : 集合被覆問題 (乱数による配置)

# 隣接行列の生成
function make_matrix(ps, r)
  size = length(ps)
  xs = zeros(size, size)
  for i = 1:size
    p1 = ps[i]
    for j = 1:size
      if distance(p1, ps[j]) < r
        xs[i, j] = 1
      end
    end
  end
  xs
end

# 集合被覆問題
function solver_set(ps, r)
  size = length(ps)
  cs = make_matrix(ps, r)   # 行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, fs[1:size], Bin)

  # 目的関数
  @objective(model, Min, sum(fs))

  # 制約条件
  for i = 1:size
    @constraint(model, sum(cs[i, j] * fs[j] for j = 1:size) >= 1)
  end

  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  println(round.(Int, value.(fs)))
  println("z ", objective_value(model))
end

関数 make_matrix() は利用者 p1 と施設 p2 の距離を計算して、それが半径 r 以下であれば要素を 1 に、そうでなければ要素を 0 にします。これで利用者が利用できる施設を行列で表すことができます。

それでは実行してみましょう。半径は 100 としました。

julia> @time solver_set(make_data(25), 100)
結果 OPTIMAL
[0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0]
z 7.0
  0.015533 seconds (4.31 k allocations: 153.633 KiB, 0.75% compilation time)

julia> @time solver_set(make_data(50), 100)
結果 OPTIMAL
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 
 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1]
z 10.0
  0.012546 seconds (5.72 k allocations: 288.016 KiB)

julia> @time solver_set(make_data(75), 100)
結果 OPTIMAL
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]
z 9.0
  0.022003 seconds (7.37 k allocations: 448.359 KiB)

julia> @time solver_set(make_data(100), 100)
結果 OPTIMAL
[0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
z 10.0
  0.010761 seconds (9.05 k allocations: 636.578 KiB)

実行環境 : Julia 1.12.6, Ubunts 24.04 LTS (WSL2, Windows 11), intel CORE i5-1235U 1.30GHz
ソルバー : JuMP v1.30.1, HiGHS v1.24.0

実行時間は高速で、利用者の個数が 100 個でもすぐに解を求めることができました。もっと時間がかかると思っていたので、こんなに速いとは大変驚きました。最適化ソルバー (HiGHS) は優秀ですね。

●最大被覆問題

施設設置問題における「最大被覆問題 (max cover problem)」は、あらかじめ設置する施設の個数 (p) を決めておいて、できるだけ多くの利用者をカバーするような設置場所を求める問題です。簡単な例題として、設置する施設の個数が 2 つで、利用者 \(D_{i}\) が利用できる施設の候補地 \(F_{j}\) が下表のように与えられているとき、利用者の総和が最大となる設置場所を求めてみましょう。

        表 : 利用可能な施設

        : f1 : f2 : f3 : f4 : f5 
--------+----+----+----+----+----
 d1 (1) :  1 :  0 :  0 :  0 :  1
 d2 (4) :  1 :  1 :  0 :  0 :  0
 d3 (5) :  0 :  1 :  1 :  0 :  0
 d4 (3) :  0 :  0 :  1 :  1 :  0
 d5 (2) :  0 :  0 :  0 :  1 :  1

    カッコ内の数字は人数

(出典 : 参考 URL 2『ロジスティクス工学 05 施設配置』)

施設 \(f_{i}\) の選択を表す変数 (Bin) を \(F_{i}\), 利用者 \(d_{i}\) が施設を利用できるか否かを表す変数 (Bin) を \(X_{i}\), 利用者の人数を \(D_{i}\) とし、上表の要素を \(W_{ij}\) で表すことにします。目的関数と制約条件は次のようになります。

目的関数 (最大化):

\( \displaystyle \sum_i D_{i} X_{i} \)

制約条件:

\(\begin{array}{l} \displaystyle \sum_j F_{j} = 2 \\ \displaystyle \sum_j W_{ij} F_{j} \leq X_{i} \quad (i= 1, 2, 3, 4, 5) \end{array}\)

目的関数は簡単ですね。制約条件ですが、最初の式は設置する施設の数を表します。次の式は、\(X_{i}\) が 1 ならば施設を利用できるので、左辺式の値は 1 以上になることを表しています。\(X_{i}\) が 0 ならば施設を利用できないので、左辺式は 0 でもいいわけです。

あとはこれを JuMP でプログラムするだけです。プログラムリストと実行結果を示します。

リスト : 最大被覆問題 (max cover)

using JuMP, HiGHS

model = Model(HiGHS.Optimizer)

# 利用可能な施設
#      f1 f2 f3 f4 f5
ws = [  1  0  0  0  1;    # d1
        1  1  0  0  0;    # d2
        0  1  1  0  0;    # d3
        0  0  1  1  0;    # d4
        0  0  0  1  1 ]   # d5

# 利用者の人数
ds = [1, 4, 5, 3, 2]

# 変数
@variable(model, fs[1:5], Bin)
@variable(model, xs[1:5], Bin)

# 目的関数
@objective(model, Max, sum(ds .* xs))

# 制約条件
@constraint(model, sum(fs) == 2)
for i = 1:5
  @constraint(model, sum(ws[i, j] * fs[j] for j = 1:5) >= xs[i])
end

# 実行
#println(model)
set_silent(model)
optimize!(model)
# 結果
println("結果 ", termination_status(model))
println(round.(Int, value.(fs)))
println(round.(Int, value.(xs)))
println("z ", objective_value(model))

結果 OPTIMAL
[0, 1, 0, 1, 0]
[0, 1, 1, 1, 1]
z 14.0

設置する施設は f2 と f4 で、d2 から d5 までの計 14 人をカバーすることができます。残念ながら d1 の人は施設を利用することはできません。

次は、利用者 (一人) を平面上にランダムに n 個配置し、半径 r 内にある施設を利用できる条件で、最大被覆問題を解いてみましょう。施設は利用者の中から選びます。プログラムは次のようになります。

リスト : 最大被覆問題 (乱数による配置)

function solver_max(ps, r, p)
  size = length(ps)
  ws = make_matrix(ps, r)   # 隣接行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, fs[1:size], Bin)
  @variable(model, xs[1:size], Bin)

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

  # 制約条件
  @constraint(model, sum(fs) == p)
  for i = 1:size
    @constraint(model, sum(ws[i, j] * fs[j] for j = 1:size) >= xs[i])
  end

  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  # println(round.(Int, value.(fs)))
  z = objective_value(model)
  println(z, " ", z / size)
end

関数 make_matrix() は集合被覆問題で作成したものと同じです。プログラムは簡単なので説明は割愛いたしｈます。

それでは実行してみましょう。半径 r は 100 で、設置する施設数 p は 7 としました。

julia> @time solver_max(make_data(25), 100, 7)
結果 OPTIMAL
24.0 0.96
  0.006732 seconds (4.18 k allocations: 187.789 KiB)

julia> @time solver_max(make_data(50), 100, 7)
結果 OPTIMAL
49.0 0.98
  0.008640 seconds (5.72 k allocations: 347.625 KiB)

julia> @time solver_max(make_data(75), 100, 7)
結果 OPTIMAL
67.0 0.8933333333333333
  0.040559 seconds (7.30 k allocations: 533.492 KiB)

julia> @time solver_max(make_data(100), 100, 7)
結果 OPTIMAL
89.0 0.89
  0.019637 seconds (8.90 k allocations: 791.602 KiB)

集合被覆問題と同じように、利用者の個数が 100 個でもすぐに解を求めることができました。興味のある方は、施設の個数や半径を変更するなど、いろいろ試してみてください。

●参考 URL

数理最適化入門（4）：施設配置の数理モデル [PDF], (田中健一さん)
ロジスティクス工学 05 施設配置 [PDF], (宮本裕一郎さん)
集合被覆問題 - Wikipedia

●プログラムリスト

#
# cover.jl : 集合被覆問題と最大被覆問題 (乱数版)
#
#            Copyright (C) 2026 Makoto Hiroi
#
using JuMP, HiGHS

# 乱数でデータを生成
function make_data(n)
  [(rand(10:490), rand(10:490)) for _ = 1:n]
end

# 距離を求める
function distance(p1, p2)
  dx = p1[1] - p2[1]
  dy = p1[2] - p2[2]
  sqrt(dx * dx + dy * dy)
end

# 隣接行列の生成
function make_matrix(ps, r)
  size = length(ps)
  xs = zeros(size, size)
  for i = 1:size
    p1 = ps[i]
    for j = 1:size
      if distance(p1, ps[j]) < r
        xs[i, j] = 1
      end
    end
  end
  xs
end

# 集合被覆問題
function solver_set(ps, r)
  size = length(ps)
  cs = make_matrix(ps, r)   # 行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, fs[1:size], Bin)

  # 目的関数
  @objective(model, Min, sum(fs))

  # 制約条件
  for i = 1:size
    @constraint(model, sum(cs[i, j] * fs[j] for j = 1:size) >= 1)
  end

  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  println(round.(Int, value.(fs)))
  println("z ", objective_value(model))
end


# 最大被覆問題
function solver_max(ps, r, p)
  size = length(ps)
  ws = make_matrix(ps, r)   # 隣接行列

  model = Model(HiGHS.Optimizer)

  # 変数の生成
  @variable(model, fs[1:size], Bin)
  @variable(model, xs[1:size], Bin)

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

  # 制約条件
  @constraint(model, sum(fs) == p)
  for i = 1:size
    @constraint(model, sum(ws[i, j] * fs[j] for j = 1:size) >= xs[i])
  end

  # 実行
  #println(model)
  set_silent(model)
  optimize!(model)
  # 結果
  println("結果 ", termination_status(model))
  # println(round.(Int, value.(fs)))
  z = objective_value(model)
  println(z, " ", z / size)
end

Julia Programming