CBC(Pulp, python-mip)ソルバーのログの読み方

Pulp[^1]やpython-mip[^2]からフリーで動かせるCBC(COIN-OR Brand-and-Cut)ソルバーのログの読み方をまとめてみました。
(個人的なメモとしてまとめたので間違いなどがあればご指摘お願いします。)

内容は以下のpdfのほぼ翻訳になります。

- CBC : COIN-OR Branch-and-Cut (A Short Guide to the CBC Command Line Interface)  
http://www.decom.ufop.br/haroldo/files/cbcCommandLine.pdf

# はじめに
上のpdfではメッセージ7までは`air03`と呼ばれるSet Partitioning問題のベンチマークを計算した時のログを使って説明されています。
メッセージ7以降は`air04'を用いて説明が行われています。

http://miplib2017.zib.de/instance_details_air03.html
http://miplib2017.zib.de/instance_details_air04.html

上のリンクから`mps`ファイルが入手できます。
python-mipの`read`を使って問題を読み込み計算することが出来ます。
(Pulpでは`read`関数が準備されていない...。`write`はありますが。)

# ログの読み方
## メッセージ1~3
![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/315019/d57ad9cb-70fe-b605-f0f3-db371056956f.png)

### メッセージ1

```
Continuous objective value is 338864 - 0.05 seconds
```

モデルの線型緩和解を解いた時の目的項の値(コスト)と計算時間が表示されています。
この例ではコストが`338864`の線型緩和解が`0.05`秒で得られています。
この値が元の問題の上界下界になります。
(最小化の場合は下界、最大化の場合は上界になります。以降では最小化問題を考えます。)
何も特別な操作を行っていないので最もラフに見積もった下界になります。

### メッセージ2

```
Cgl0004I processed model has 120 rows, 8456 columns (8456 integer) and 71651 elements
```

`Cgl`と記載された行では**前処理**が行われています。
`Cgl`はCut Generation Library[^3]の略のようです。

前処理で削減した後の制約条件(`rows`)、変数(`columns`)、非ゼロ要素の数(`elements`)が最後の行に表示されています。
前処理を行うと以下のように問題のサイズが小さくなっています。

|   |  制約条件  | 変数| 非ゼロ要素 |
| ---- | ---- | ---- | ---- |
|  前処理前  |  823  | 8904| 72965 |
|  前処理後  |  120  | 8456| 71651 |

前処理で問題のサイズが小さくなればなるほどMIPソルバーにとって解きやすい定式化になっているそうです。

### メッセージ3

```
Cbc0038I Pass 1: suminf. 8.33333 (22) obj.  341524 iterations 106
Cbc0038I Pass 2: suminf. 8.33333 (22) obj.  341524 iterations 3
Cbc0038I Pass 3: suminf. 8.33333 (22) obj.  341524 iterations 70
Cbc0038I Pass 4: suminf. 7.20000 (20) obj.  342390 iterations 75
Cbc0038I Pass 5: suminf. 1.50000 (3)  obj.  343697 iterations 45
Cbc0038I Pass 6: suminf. 1.50000 (3)  obj.  343697 iterations 55
       ...
Cbc0038I Pass 12: suminf. 0.00000 (0) obj.  362176 iterations 144
Cbc0038I Solution found of 362176
```

`Cbc0038I`が始まったことを表しています。
`Cbc0038I`では**Feasibility Pump (M. Fischetti, Glover, & Lodi, 2005)[^4],[^5]**と呼ばれる手法を用いて実行可能解となる初期解(暫定解)を探索しているようです。
(`I`はInitialの略でしょうか。)
12のPassが終了した後に`Cbc0038I Solution found of 362176`と表示されており、コストが`362176`の暫定解が得られたことが分かります。
メッセージ1に表示されていた下界が`338864`だったので、少なくとも`[338864, 362176]`の範囲に最適解(厳密解)があることが分かります。
暫定解のコストと下界との差(ギャップ)が小さいほどCBCの性能は良くなります。

## メッセージ4~6
![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/315019/4be9f4a6-0efe-97e2-cbde-c1bbcbf96132.png)

### メッセージ4

```
Cbc0012I Integer solution of 340160 found by DiveCoefficient after 14 iterations and 0 n odes (0.97 seconds)
```

**DiveCoefficient[^5],[^6]**と呼ばれる前処理の結果が表示されています。
(私は詳細を知りません。)
その結果さらに暫定解が更新されてギャップが`[338864, 340160]`の区間になっています。

### メッセージ5

```
Cbc0013I At root node, 5 cuts changed objective from 338864.25 to 340160 in 2 passes
```

カット系のアルゴリズムが実行されているようです。
特に**双対問題**の上界下界を改善するために、得られた緩和解の小数値を取り除くカットをいくつか試しているみたいです。
この例では5種類のカットが実行されたようです。
詳しい内容がログの下の方に表示されています。

```
Cbc0014I Cut generator 0 (Probing) - 0 row cuts average 0.0 elements, 160 column cuts (160 active) in 0.840 seconds - new frequency is 1
Cbc0014I Cut generator 1 (Gomory) - 4 row cuts average 1257.2 elements, 0 column cuts (3 active) in 0.020 seconds - new frequency is -100
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.020 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 10 row cuts average 3.7 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is 1
Cbc0014I Cut generator 6 (TwoMirCuts) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.030 seconds - new frequency is -100
```

この問題では0,1,3番目のカットで十分のようです。
(Gomory Cut[^7]とClique Cutが特に有効であることが分かります。)
双対問題の上界が`340160`であることが判定できました。
暫定解のコストも`340160`であるので最適解であることが分かります。

### メッセージ6

```
Result - Finished objective 340160 after 0 nodes and 11 iterations - took 4.75 seconds (total time 5.06)
Total time 5.14
```

最終的なコストの値と計算時間が表示されています。
今回の問題はroot nodeで十分良い実行可能解や双対問題の上界が得られているのでこれ以上分枝カット法の枝刈りは行われていません。


## メッセージ7, 8

![fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/315019/9d60857d-b17c-b422-5568-b086ae9332ea.png)

これらのメッセージはより難しい問題(`air04`)を計算したときに表示されるログになります。
`air04`では分枝限定法の前にこの探索が行われているようです。
分子限定法の木構造の探索で探索したノードの数、探索がまだ終わっていないノードの数がそれぞれ表示されています。
新しい整数が得られるたびにそのコストや用いた手法に関する情報が表示されています。(メッセージ8)

以上がログの読み方になります。
参考までにpython-mipを使って`air03`, `air04`を計算した時のログを以下にUpしました。
https://github.com/Noriaki416/Qiita/tree/master/CBC_log

`air0#_cbc_log.txt`がログになります。

# チューニング
ログを読むことでどのプロセスが問題を解く上でネックになっているのかが分かります。
その原因を特定してチューニングを行うことで求解性能を改善できるはずです。
冒頭のpdfのp6以降に詳細なチューニングパラメータが記載されています。

Pulpを使ってチューニングする場合は`options`で指定すればいいみたいです。
下のブログで詳しく解説されています。

https://inarizuuuushi.hatenablog.com/entry/2019/03/07/090000

> optionsではそれ以外のCBCオプションを設定できます.
> 例えばmaxsolオプションを使用する場合は options = ['maxsol 1']とします.
> maxsolオプションは, 暫定解が何個見つかった時に計算を終了するかを指定できます. 最適解でなくても実行可能解が１つ欲しい場合はmaxsol 1とすればより速い時間で求めることができます.

`maxsol`以外のパラメータについてもGAMS[^8]のサイトに記載されています。
色々弄ってみると違いが見えて面白いのかもしれません。

以上になります。
ご指摘などありましたら遠慮なくお願いします。

[^1]: https://pypi.org/project/PuLP/
[^2]: https://www.python-mip.com
[^3]: https://projects.coin-or.org/Cgl
[^4]: http://www.dei.unipd.it/~fisch/papers/feasibility_pump.pdf
[^5]: https://www.zib.de/berthold/primal_heuristics.pdf (SCIPのpdfですがFeasibility PumpやDiving~に関する説明が記述されています。)
[^6]: https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiZ9Y-RhoDtAhUDEqYKHZJvD58QFjABegQIBhAC&url=https%3A%2F%2Fopus4.kobv.de%2Fopus4-zib%2Ffiles%2F1029%2FBerthold_Primal_Heuristics_For_Mixed_Integer_Programs.pdf&usg=AOvVaw3F8ik1VHxnn_3GdP_rk89N
[^7]: http://dopal.cs.uec.ac.jp/okamotoy/lect/2013/opt/handout06.pdf
[^8]: https://www.gams.com/latest/docs/S_CBC.html#CBC_general