PuLPで彩色問題を解く

内容

下の東京都の白地図を、隣り合う自治体で色が同じにならないように4色で塗り分けます。
塗り方はPuLPという数理計画ソルバを用いて計算します。

前提

• 色は赤、緑、青、黄の4色。
• 地方自治体は23区 + 26市 + 奥多摩町、日の出町、瑞穂町、檜原村の計53。島嶼部は除きます。
• 隣接は調べながら自力で書き出しました（ここが一番つらかった。）。江東区が港区、品川区、大田区と接しているのは、埋立地なので、無視しても良かったのですが、平面グラフとしての矛盾は生じなかったので現実に即す形にしました。特に飛び地のみで隣接する自治体はなかったかと思います。

colors = ['R', 'G', 'B', 'Y']

areas = [
    'Adachi', 'Arakawa', 'Bunkyo', 'Chiyoda', 'Chuo', 'Edogawa',
    'Itabashi', 'Katsushika', 'Kita', 'Koto', 'Meguro', 'Minato',
    'Nakano', 'Nerima', 'Ota', 'Setagaya', 'Shibuya', 'Shinagawa',
    'Shinjuku', 'Suginami', 'Sumida', 'Taito', 'Toshima',
    'Akishima', 'Akiruno', 'Chofu', 'Fuchu', 'Fussa', 'Hachioji',
    'Hamura', 'Higashikurume', 'Higashimurayama', 'Higashiyamato',
    'Hino', 'Inagi', 'Kiyose', 'Kodaira', 'Koganei', 'Kokubunji',
    'Komae', 'Kunitachi', 'Machida', 'Mitaka', 'Musashimurayama',
    'Musashino', 'Nishitokyo', 'Ome', 'Tachikawa', 'Tama',
    'Hinode', 'Mizuho', 'Okutama', 'Hinohara'
]

neighborings = [
    [0, 1], [0, 7], [0, 8], [0, 20], [1, 2], [1, 8], [1, 20], [1, 21],
    [2, 3], [2, 8], [2, 18], [2, 21], [2, 22], [3, 4], [3, 11], [3, 18],
    [3, 21], [4, 9], [4, 11], [4, 20], [4, 21], [5, 7], [5, 9], [5, 20],
    [6, 8], [6, 13], [6, 22], [7, 20], [8, 22], [9, 11], [9, 14], [9, 17],
    [9, 20], [10, 14], [10, 15], [10, 16], [10, 17], [11, 16], [11, 17],
    [11, 18], [12, 13],  [12, 16], [12, 18], [12, 19], [12, 22], [13, 19],
    [13, 22], [13, 44], [13, 45], [14, 15], [14, 17], [15, 16], [15, 19],
    [15, 25], [15, 39], [15, 42], [16, 17], [16, 18], [16, 19], [18, 22],
    [19 ,42], [19, 44], [20, 21], [23, 27], [23, 28], [23, 33], [23, 47],
    [24, 27], [24, 28], [24, 29], [24, 46], [24, 49], [24, 51], [24, 52],
    [25, 26], [25, 34], [25, 37], [25, 39], [25, 42], [26, 33], [26, 34],
    [26, 37], [26, 38], [26, 40], [26, 48], [27, 28], [27, 29], [27, 43],
    [27, 47], [27, 50], [28, 33], [28, 41], [28, 48], [28, 52], [29, 46],
    [29, 50], [30, 31], [30, 35], [30, 36], [30, 45], [31, 32], [31, 35],
    [31, 36], [32, 36], [32, 43], [32, 47], [33, 40], [33, 47], [33, 48],
    [34, 48], [36, 37], [36, 38], [36, 45], [36, 47], [37, 38], [37, 42],
    [37, 44], [37, 45], [38, 40], [38, 47], [40, 47], [41, 48], [42, 44],
    [43, 47], [43, 50], [44, 45], [46, 49], [46, 50], [46, 51], [51, 52]
]

考え方

自分で思いついたわけではなく、https://www.jstage.jst.go.jp/article/bjsiam/23/2/23_KJ00008726325/_pdf を参考にしただけですが、以下の通り考えることで、整数計画問題と捉えることができます。
$N$個の領域に対して使える色が$M$色あるとする。
まず、それぞれの領域 $i\in\mathbb{N}\ (1\le i\le N)$に関して、どの色が使われるかを4つのbool値(1のとき、その色で領域が塗られている。)で管理します。

x_{ic} \in \{0,1\}\ (1\le c\le M)\\ 
\sum_{c=1}^M x_{ic}= 1\\

また $i,j\ (i\ne j)$ が隣接しているとき、

x_{ic} + x_{jc}\le1\ (1\le c\le M)

が成り立っている必要があります。これらは全て一次方程式、一次不等式でかけるので、目的関数を伴わない整数計画問題となります。

コード

PuLPの使い方は、他記事を読んでいただくとして、先ほどの考え方をコードにすると以下のようになります。

import pulp

# colors, areas, neighboringsの宣言部分

problem = pulp.LpProblem('Four Colors', pulp.LpMinimize)

var = {
    area: pulp.LpVariable.dicts(area, colors, cat = pulp.LpBinary) for area in areas
}

for area in areas:
    problem += sum(var[area].values()) == 1

for pair in neighborings:
    for color in colors:
        problem += var[areas[pair[0]]][color] + var[areas[pair[1]]][color] <= 1

result = problem.solve()

if pulp.LpStatus[result] == 'Optimal':
    output = { color: [] for color in colors }
    for area, v in var.items():
        for color, flg in v.items():
            if flg.value() == 1:
                output[color].append(area)
                break
    print(output)
else:
    print('解なし')

結果

https://n.freemap.jp/ で色塗りしました。

補足

3色の場合解なしでした。