最長しりとり問題を混合整数線形計画法として定式化しMATLABで解く

与えられた単語から最も長いしりとりを作る問題（以下、「最長しりとり問題」）を考えます。ここでは、最長しりとり問題をMILP（混合整数計画法）として定式化します。また、MATLABのOptimization ToolboxでMILPを実装し、小田急の駅（70個）とアイドルマスター全員（325人）について解いていきます。

## 最長しりとり問題の例
与えられた単語から、なるべく長くしりとりが続くように、選択・並び替えをします。
```plaintext
1:ごりら
2:ぱせり
3:りんご
4:らっぱ
5:ぱん
```
が与えられたとしたら、最長しりとりの一例として
$$りんご\to ごりら\to らっぱ\to ぱん$$
などが挙げられます。

## 単語リスト
ひらがなで書かれた単語リストをテキストファイル形式で用意し、読み込みます。実際の例は本記事下端の付録に示します。単語数を$N$とします。
```matlab
%% 単語を読み込む
words_string = readlines("---.txt");
N = length(words_string);
heads = zeros(N,1);
tails = zeros(N,1);
for i=1:N
    word_char = char(words_string(i));
    heads(i) = word_char(1);
    tails(i) = word_char(end);
end
disp(N)
```

## ノードと隣接行列
ここでは、しりとりを$N+2$個のノードからなる有向グラフとして表すこととします。ノードをそれぞれ次のように割り当てます。
* $1$から$N$：与えられた単語
* $S:=N+1$ : しりとりの開始を表すもの（「開始ノード$S$」）
* $E:=N+2$：しりとりの終了を表すもの（「終了ノード$E$」）

<details><summary>ノードの例</summary>

次のように$N=5$個の単語が与えられるとします。
```plaintext
1:ごりら
2:ぱせり
3:りんご
4:らっぱ
5:ぱん
```
この時、次のしりとりは
$$\texttt{START}\longrightarrow りんご\to ごりら\to らっぱ\to ぱん\longrightarrow \texttt{END}$$
次のようなノードの移動として表されます。
$$6 \longrightarrow 3\to1\to4\to5 \longrightarrow 7$$

---
</details>

次に、隣接行列$G\in\mathbb R^{ (N+2)\times(N+2) }$を作ります。隣接行列は、しりとりでの単語の接続の可否を表します。具体的には、
$$G_{ij}=\begin{cases}
1& i\ne j\text{かつ単語}i\text{の後に単語}j\text{を続けられる} \\\\
0& \text{それ以外}
\end{cases}$$
と定めます。
ただし、開始ノードの後には全ての単語を続けられるので、
$$G_{S j} = 1, \quad\forall j\in\\{1,\dots,N\\}$$
また、全ての単語は終了ノードへ続けられるので、
$$G_{i E} = 1, \quad\forall i\in\\{1,\dots,N\\}$$
とします。

<details><summary>隣接行列の例</summary>

```plaintext
1:ごりら
2:ぱせり
3:りんご
4:らっぱ
5:ぱん
```
この時、隣接行列は
$$\begin{bmatrix}
0&0&0&🦍📯&0 &0&1\\\\
0&0&🌿🍎&0&0 &0&1\\\\
🍎🦍&0&0&0&0 &0&1\\\\
0&📯🌿‬&0&0&📯🍞 &0&1\\\\
0&0&0&0&0 &0&1\\\\
1&1&1&1&1 &0&0\\\\
0&0&0&0&0 &0&0
\end{bmatrix}$$
です（絵文字が入っているところは全て$1$）。

---
</details>


```matlab
%% 隣接行列G(i,j) = 単語iの後に単語jを続けられれば1、そうでなければ0
% N+1 : 開始ノード
% N+2 : 終了ノード
START_INDEX = N+1;
END_INDEX = N+2;

G = zeros(N+2);
for i=1:N
for j=1:N
    if i==j
        continue
    end
    if tails(i)==heads(j)
        G(i,j) = 1;
    end
end
end
G(START_INDEX,1:N) = 1;%開始ノードの後には全ての単語を続けられる
G(1:N,END_INDEX) = 1;%すべての単語の後に終了ノードを続けられる
n_connectable = sum(G,"all")
```

## 最適化問題の構築
隣接行列$G$のうち、成分が1のインデックス$(i,j)$の集合を$M$とします。$M$は有向グラフのエッジの集合であると共に、しりとりにおいて接続可能な単語の組を表していることになります。
最長しりとりは、この$M$の中からうまく単語の組を取り出し、なるべく長く繋げていく問題です。以下より、最適化問題を示します。

ちなみに、上のコードの`n_connectable`$=|M|$です。

### 変数
$|M|$個の整数の変数$x_{ij}\in\\{0,1\\},~(i,j)\in M$を用意し、
$$ x_{ij} = \begin{cases}
1 & 単語の組(i,j)を使う。つまり、単語iの後に単語jを続けることとする。 \\\\
0 & 使わない
\end{cases}$$
という意味とします。これら$|M|$個の変数を最適化により決定することで、しりとりを作ります。

<details><summary>変数の例</summary>

```plaintext
1:ごりら
2:ぱせり
3:りんご
4:らっぱ
5:ぱん
```
という単語リストがあったとします。この時、「ごりら→らっぱ」という単語の組を使うのなら$x_{14}=1$、使わないのなら$x_{14}=0$となります。

---
</details>

```matlab
%% 最適化変数x = 単語接続を表す。単語iの後に単語jを続ければ1、そうでなければ0。
x_name = @(i,j) "i"+i+"j"+j;

% {m}:m番目の組の{i; j; xの名前; x}
idx_list = cell(n_connectable, 4);
% フィールド名:xの名前
x_struct = struct;
% {i}:単語iから続く組(i,*)のリスト
wherei_list = cell(N+2,1);
% {j}:単語jへ続く組(*,j)のリスト
wherej_list = cell(N+2,1);

m = 1;
for i=1:N+2
for j=1:N+2
if G(i,j) == 1
    name = x_name(i,j);
    x = optimvar( ...
        name, ...
        1, ...
        Type="integer", ...
        LowerBound=0,UpperBound=1);

    idx_list{m,1} = i;
    idx_list{m,2} = j;
    idx_list{m,3} = name;
    idx_list{m,4} = x;
    m = m + 1;
    x_struct.(name) = x;
    wherei_list{i} = [wherei_list{i};x];
    wherej_list{j} = [wherej_list{j};x];
end
end
end
```
後から変数を探しやすくするために、ここで色んな変数を作っておきます💦

### 問題と目的関数
なるべく長いしりとりを作るということは、なるべく多くの単語の組を使うということなので、問題と目的関数は次のようになります。
$$\max_{ x_{ij} ~,~ (i,j)\in M } \sum_{(i,j)\in M}x_{ij}$$

```matlab
prob = optimproblem(ObjectiveSense="maximize");
% 目的関数
prob.Objective = sum(struct2array(x_struct))
```
$x_{ij}$がそれぞれ整数であり、また、制約（以下より説明します）が全て線形制約となっています。したがって、この問題は混合整数線形計画問題となります。

### 制約1. 開始ノードから続く単語は1つ
開始ノードからは任意の単語へ続くことができますが、その回数は1回だけです。これは次のように表されます。
$$ \sum_{(i,j)\in M, i=S} x_{ij} =1 $$
```matlab
from_start_x_list = wherei_list{START_INDEX};
prob.Constraints.Start = sum(from_start_x_list)==1
```

### 制約2. 終了ノードへ続く単語は1つ
任意の単語から終了ノードへ続くことができますが、その回数は1回だけです。これは次のように表されます。
$$ \sum_{(i,j)\in M, j=E} x_{ij} =1$$
```matlab
to_end_x_list = wherej_list{END_INDEX};
prob.Constraints.End = sum(to_end_x_list)==1
```

### 制約3. 単語iから続く単語の総数は1以下
しりとりは、有向グラフ$G$上を一筆書きで辿っていくのと同じです。ですので、同じ単語から2回以上出ていってはいけません。
$$ \sum_{(i,j)\in M, i=i^\ast} x_{ij} \le 1 
,\quad\forall i^\ast\in\\{1,\dots,N\\} $$
```matlab
for i=1:N
from_i_x_list = wherei_list{i};
if ~isempty(from_i_x_list)
    prob.Constraints.("From"+i) ...
        = sum(from_i_x_list) <= 1;
end
end
```

### 制約4. 単語jへ続く単語の総数は1以下
制約3と同様、同じ単語へ2回以上入っていってはいけません。
$$ \sum_{(i,j)\in M, j=j^\ast} x_{ij} \le 1 
,\quad\forall j^\ast\in\\{1,\dots,N\\} $$
```matlab
for j=1:N
to_j_x_list = wherej_list{j};
if ~isempty(to_j_x_list)
    prob.Constraints.("To"+j) ...
        = sum(to_j_x_list) <= 1;
end
end
```

### 制約5. 単語の接続について
単語$k$へ続く単語が有るのなら、単語$k$から続く単語も有るべきです。逆に、
単語$k$へ続く単語が無いのなら、単語$k$から続く単語も無いべきです。

制約3と4より、単語$k$へ続く総回数と単語$k$から続く総回数は$\\{0,1\\}$のいずれかなので、この制約は次のように表されます。
$$ \sum_{(i,j)\in M ~,~ j=k} x_{ij}=
\sum_{(i,j)\in M ~,~ i=k} x_{ij}
~,\quad k\in\\{1,\dots,N\\}
$$
左辺が任意の単語から単語$k$へ続く総回数、右辺が単語$k$から任意の単語へ続く総回数を表しています。

```matlab
f = waitbar(0);
for k=1:N
    waitbar(k/N, f, ""+k+"/"+N+"");

    to_k_x_list = wherej_list{k};
    to_k_count = sum(to_k_x_list);
    from_k_x_list = wherei_list{k};
    from_k_count = sum(from_k_x_list);

    prob.Constraints.("Connection"+k) ...
        = to_k_count == from_k_count;
end
close(f);
```
このコードはちょっと重いので、ループでプログレスバーを表すようにしています。

## 実験１：小田急の駅
まずは小田急の駅（$N=70$個）を与えてみます。接続可能な単語の組数は$|M|=225$となりました。なので、225個の変数の最適化問題となります。
```matlab
%% 求解する
options_intlinprog = optimoptions("intlinprog", ...
    Display="iter", ...
    RelativeGapTolerance=1e-6, ...
    IntegerTolerance=1e-6);
sol = solve(prob, Options=options_intlinprog);
```
なんか解けちゃいました。MATLAB先輩、さすがです。。。
```matlab
%% 結果をまとめる
% R(i,j) = 単語iの後に単語jを続けていれば1、そうでなければ0。
R = zeros(N+2);
for m = 1:n_connectable
    i = idx_list{m,1};
    j = idx_list{m,2};
    name = idx_list{m,3};
    x = sol.(name);
    R(i,j) = x;
end

% 結果を見る
for i0 = find(R(START_INDEX,:))
i = i0;
i_history = [i];
while true
    disp(i+":"+words_string(i))
    i = find(R(i,:),1);
    if isempty(i)
        disp("これに続く単語が存在しない")
        break
    end
    if find(i_history == i)
        disp("無限ループになっている")
        break
    end
    if i == END_INDEX
        disp("終了ノードに到達した つながった単語の数:"+length(i_history))
        break
    end
    i_history = [i_history i];
end
end
```
このコードを実行してしりとりの結果を見ると・・・？
```plaintext
66:ふじさわほんまち
12:ちとせふなばし
23:しんゆりがおか
42:かいせい
20:いくた
26:たまがわがくえんまえ
32:えびな
終了ノードに到達した つながった単語の数:7
```
素晴らしい。なんかいい感じのしりとりになっています。Rustで実装した、深さ優先探索を使った最長しりとりソルバー（別記事でまた書こうと思います）で確認したところ、このしりとりは確かに最長であることが分かりました。

## 実験２：アイドルマスター全員
次は、アイドルマスターのキャラクター全員（$N=325$）でやってみます。接続可能な単語の組数は$|M|=3525$でした🫠🫠🫠多すぎる

以下に結果を示します！
```plaintext
243:ふじい とも
297:ももせ りお
69:おとくら ゆうき
86:きた ひなこ
110:こうさか うみ
282:みふね みゆ
312:ゆめみ りあむ
287:むかい たくみ
284:みやお みや
301:やお ふぇいふぇい
44:いむら せつな
199:なかの ゆか
80:かわしま みずき
97:きよすみ くろう
46:うじいえ むつみ
273:みうら あずさ
127:さぎさわ ふみか
77:かみじょう はるな
207:なんば えみ
283:みむら かなこ
119:ころん くりす
161:すずみや せいか
74:かたぎり さなえ
51:えびはら なほ
256:ほんだ みお
57:おおさき あまな
198:ながとみ はすみ
277:みずの みどり
318:りゅうざき かおる
終了ノードに到達した つながった単語の数:29
```
29人繋げることができました。また、無事に終了ノードまで至れています。

ですが、あまり性能が出ていませんね・・・というのも、アイドルマスターでは少なくとも100人以上しりとりを繋げられることが、Rustの方のプログラムで確認できているのです。
局所的最適解があれば（あるのかな🤔）それにハマっていたり、変数があまりにも多すぎてソルバーの性能の限界に来ていたり、そもそもの定式化の質が低いことが考えられます。

## おわりに
最長しりとり問題をMILP（混合整数線形計画法）として定式化し、MATLABに解いてもらいました。
単語数が少ない時には高精度の結果が得られることを確認できた一方で、単語数が多い時は精度に難ありとなりました。もっと良いソルバー（MATLABより良いソルバーなんてないと思うけど）や、もっと良い制約条件等の組み方をして、精度が上がるか検証する必要があると思います。

これとは別に、最長しりとり問題を深さ優先探索で求めるプログラムをRustで作っています。これについては後日また別の記事で出そうと思います。

あと、個人的に、最適化問題は数学で表すよりもMATLABコードで表したほうが説明的で分かりやすいですね。

## 付録：単語リストテキストファイル
```plaintext:odakyu.txt
しんじゅく
みなみしんじゅく
さんぐうばし
よよぎはちまん
よよぎうえはら
ひがしきたざわ
しもきたざわ
せたがやだいた
うめがおか
ごうとくじ
きょうどう
ちとせふなばし
そしがやおおくら
せいじょうがくえんまえ
きたみ
こまえ
いずみたまがわ
のぼりと
むこうがおかゆうえん
いくた
よみうりらんどまえ
ゆりがおか
しんゆりがおか
かきお
つるかわ
たまがわがくえんまえ
まちだ
さがみおおの
おだきゅうさがみはら
そうぶだいまえ
ざま
えびな
あつぎ
ほんあつぎ
あいこういしだ
いせはら
つるまきおんせん
とうかいだいがくまえ
はだの
しぶさわ
しんまつだ
かいせい
かやま
とみず
ほたるだ
あしがら
おだわら
さつきだい
くりひら
くろかわ
はるひの
おだきゅうながやま
おだきゅうたませんたあ
からきだ
ひがしりんかん
ちゅうおうりんかん
みなみりんかん
つるま
やまと
さくらがおか
こうざしぶや
ちょうご
しょうなんだい
むつあいにちだいまえ
ぜんぎょう
ふじさわほんまち
ふじさわ
ほんくげぬま
くげぬまかいがん
かたせえのしま
```

アイドルマスターについては、次のページをごにょごにょやって名前を取り出してきました。

https://idollist.idolmaster-official.jp/search