GoogleのOR-Toolsで数理最適化モデルの練習問題を解く(4)ナンバープレースを解く

#はじめに

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」の練習問題を解く第四回目の記事です。

第一回目の記事は以下です。まずはこちらをご覧下さい。
>GoogleのOR-Toolsで数理最適化モデルの練習問題を解く(1)一番易しいマス埋め問題
>https://qiita.com/ttlabo/private/7e8c93f387814f99931f

#ナンバープレース(ナンプレ)を解く

参考テキストに出てくる練習問題の3番目です。
以下の問題をさっそく行ってみましょう。

>:speech_balloon: 問題

```ruby:7.3.問題
下図のナンバープレース(ナンプレ)を解け。但し、条件は以下とする。
 ①9x9の全マスに、1～9の数字を必ず埋める。
 ②どの1行、どの1列、どの3x3においても、同じ数字は1回だけしか使えない。
 ③数字が埋まっている場所では、その数字を使う。
```

![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/c97223ed-57c7-7043-1696-8e3b2955ba6a.jpg)


>:question: **考え方**

■数理モデルを考える
図2のような、9x9x9の729個の箱を考えます。
軸方向を行、列、セル数字とします。セル数字とは、ナンプレのそれぞれの枠に入る数字という意味です。

![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/93707c8b-f81c-50c5-64b5-02a76cebcf3c.jpg)

1つの箱は、セル数字が選ばれているか選ばれていないかのどちらかの状態を持つとします。

例えば、1行1列の下から1番目の箱は、ナンプレの1行目1列目セルにおいて数字1が選ばれている場合は1、選ばれていない場合は0を持つとします。
これにならい、1行目1列目セルの下から2番目の箱は数字2が選ばれているかいないか、同様に一番上の箱は数字9が選ばれているかいないかを表します。
これらが9行目9列目まで同じように連続します。
一般化して、x行y列セル数字zの箱が1の場合はナンプレのx行目y列目のマスの数字がzであることを意味します。

このような箱を想定し、定式化を行います。

>:a: **解答**

この729個の一つ一つの箱が0または1になります。
制約条件は以下です。

(1)ナンプレの1つのマスに数字は1つだけ入る。つまり1x1x9(行x列xセル数字)の9箱の合計値が1である。
(2)1行のマスに同じ数字は1つだけ入る。9x1x1の9箱の合計値が1
(3)1列のマスに同じ数字は1つだけ入る。1x9x1の9箱の合計値が1
(4)3x3マスに同じ数字は1つだけ入る。3x3x1の9箱の合計値が1
　※3x3マスとは、問題文付属図の太線で示されている9つのマスです
(5)数字が初めから指定されているマス(問題文付属図)では、その数字とする。1x1x1の1箱の合計値が1

プログラムを検討します。
プログラムの記載内容については、基本的にGoogleのOR-Toolsのガイドに沿っています。(https://developers.google.com/optimization)

プログラムの冒頭におまじないを書きます。

```ruby:7.4.renshu.py
from __future__ import print_function
from ortools.linear_solver import pywraplp
import numpy as np
```

混合整数計画ソルバーで解きますので、以下に宣言します。

```ruby:7.4.renshu.py
# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
```

問題に提示されているナンプレの9x9マスを以下の2次元配列preに定義します。

```ruby:7.4.renshu.py
pre = [[0,0,6,0,0,0,0,0,1],
       [0,7,0,0,6,0,0,5,0],
       [8,0,0,1,0,3,2,0,0],
       [0,0,5,0,4,0,8,0,0],
       [0,4,0,7,0,2,0,9,0],
       [0,0,8,0,1,0,7,0,0],
       [0,0,1,2,0,5,0,0,3],
       [0,6,0,0,7,0,0,8,0],
       [2,0,0,0,0,0,4,0,0]]
```

9x9x9個の箱を3次元の配列変数diceと宣言します。

```ruby:7.4.renshu.py
dice = {}
```

9x9x9個の箱(サイコロ箱)について、それぞれの箱を変数で定義します。
行を１～9まで変数xでループさせます。
列を１～9まで変数yでループさせます。
セル数字方向に対して変数zでループさせます。

```ruby:7.4.renshu.py
# サイコロ箱の変数定義
for x in range(9):
    for y in range(9):
        for z in range(9):
            dice[x,y,z] = solver.BoolVar("dice%i%i%i" % (x,y,z))
```


以下、制約条件を定義します。

**制約条件(1)**
`ナンプレの1つのマスに数字は1つだけ入る。つまり1x1x9(行x列xセル数字)の9箱の合計値が1である。`
合計値を取得する際には、solver.Sum()を使用します。合計の対象となる値はdiceのz方向の1マスから9マスまでで、内包表記でzを１～9までループさせています(実際にはインデックス値で0～8までループ)。

```ruby:7.4.renshu.py
for x in range(9):
    for y in range(9):
        solver.Add(solver.Sum([dice[x,y,z] for z in range (9)]) == 1)
```

**制約条件(2)**
`1行のマスに同じ数字は1つだけ入る。9x1x1の9箱の合計値が1`

```ruby:7.4.renshu.py
for y in range(9):
    for z in range(9):
        solver.Add(solver.Sum([dice[x,y,z] for x in range (9)]) == 1)
```

**制約条件(3)**
`(3)1列のマスに同じ数字は1つだけ入る。1x9x1の9箱の合計値が1`

```ruby:7.4.renshu.py
for x in range(9):
    for z in range(9):
        solver.Add(solver.Sum([dice[x,y,z] for y in range (9)]) == 1)
```

**制約条件(4)**
`(4)3x3マスに同じ数字は1つだけ入る。3x3x1の9箱の合計値が1`

```ruby:7.4.renshu.py
# z軸方向にループ
for z in range(9):
    # 3x3マスの中心セルがとるx座標をループ
    for m1 in [1,4,7]:
        # 3x3マスの中心セルがとるy座標をループ
        for m2 in [1,4,7]:
            solver.Add(solver.Sum([dice[x,y,z] for x in [m1 - 1,m1,m1 + 1] for y in [m2 - 1,m2,m2 + 1]]) == 1)
```

**制約条件(5)**
`(5)数字が初めから指定されているマス(問題文付属図)では、その数字とする。1x1x1の1箱の合計値が1`

```ruby:7.4.renshu.py
for x in range(9):
    for y in range(9):
        if pre[x][y] != 0:
            z = pre[x][y]
            solver.Add(dice[x,y, z - 1] == 1)
```

制約条件は以上です。
解を求めます。
今回は、最大値化または最小値化する必要はありません。

```ruby:7.4.renshu.py
# 求解実行
status = solver.Solve()
```

結果を検証します。

```ruby:7.4.renshu.py
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())

    # 表示用の配列pを作る
    p = []
    for x in range(9):
        for y in range(9):
            for z in range(9):
                if dice[x,y,z].solution_value() > 0.5:
                    p.append(z + 1)

    print(np.reshape(p, (9, 9)))                    
    print("Time = ", solver.WallTime(), " milliseconds")
else:
    print('The problem does not have an optimal solution.')
```

今回の実行では、以下のような解が得られました。
>Solution:
>ok
>Objective value = 0.0
>[[5 3 6 8 2 7 9 4 1]
> [1 7 2 9 6 4 3 5 8]
> [8 9 4 1 5 3 2 6 7]
> [7 1 5 3 4 9 8 2 6]
> [6 4 3 7 8 2 1 9 5]
> [9 2 8 5 1 6 7 3 4]
> [4 8 1 2 9 5 6 7 3]
> [3 6 9 4 7 1 5 8 2]
> [2 5 7 6 3 8 4 1 9]]
>Time =  477  milliseconds

上記でナンプレが解けました。
これをマスに記載したものが以下の図です。
ピンク色の部分は問題出題時に番号が予め決められていたセルです。

![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/0c19b46d-c976-4083-d5ec-be7bea6ab844.jpeg)


以下にプログラムの全体を記載します。

```ruby:7.4.renshu.py
from __future__ import print_function
from ortools.linear_solver import pywraplp
import numpy as np

# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

pre = [[0,0,6,0,0,0,0,0,1],
       [0,7,0,0,6,0,0,5,0],
       [8,0,0,1,0,3,2,0,0],
       [0,0,5,0,4,0,8,0,0],
       [0,4,0,7,0,2,0,9,0],
       [0,0,8,0,1,0,7,0,0],
       [0,0,1,2,0,5,0,0,3],
       [0,6,0,0,7,0,0,8,0],
       [2,0,0,0,0,0,4,0,0]]


# サイコロ箱の変数定義
dice = {}
for x in range(9):
    for y in range(9):
        for z in range(9):
            dice[x,y,z] = solver.BoolVar("dice%i%i%i" % (x,y,z))

# 制約条件1 - 1マスに数字は一つ(z方向にループした際に1が入る箱は9個のうち1個)
for x in range(9):
    for y in range(9):
        solver.Add(solver.Sum([dice[x,y,z] for z in range (9)]) == 1)

# 制約条件2 - 縦(y方向)に同じ数字は一つ(y方向にループした際に1が入る箱は9個のうち1個)
for x in range(9):
    for z in range(9):
        solver.Add(solver.Sum([dice[x,y,z] for y in range (9)]) == 1)

# 制約条件3 - 横(x方向)に同じ数字は一つ(x方向にループした際に1が入る箱は9個のうち1個)
for y in range(9):
    for z in range(9):
        solver.Add(solver.Sum([dice[x,y,z] for x in range (9)]) == 1)

# 制約条件4 - 3x3マスに同じ数字は一つ
# z軸方向にループ
for z in range(9):
    # 3x3マスの中心セルがとるx座標をループ
    for m1 in [1,4,7]:
        # 3x3マスの中心セルがとるy座標をループ
        for m2 in [1,4,7]:
            solver.Add(solver.Sum([dice[x,y,z] for x in [m1 - 1,m1,m1 + 1] for y in [m2 - 1,m2,m2 + 1]]) == 1)

# 制約条件5 - 数字が入っているマスは、その数字
for x in range(9):
    for y in range(9):
        if pre[x][y] != 0:
            z = pre[x][y]
            #print('z=',z)
            solver.Add(dice[x,y, z - 1] == 1)

# 求解実行
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())

    # 表示用の配列pを作る
    p = []
    for x in range(9):
        for y in range(9):
            for z in range(9):
                if dice[x,y,z].solution_value() > 0.5:
                    p.append(z + 1)

    print(np.reshape(p, (9, 9)))                    
    print("Time = ", solver.WallTime(), " milliseconds")
else:
    print('The problem does not have an optimal solution.')
```

#出展

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」に記載の練習問題を参考にさせて頂いております。

■参考テキスト
「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」
　斉藤努[著]
　近代科学社[出版]

![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/df8e6661-9dd6-559b-5ae8-0bf79f73d1a5.jpeg)

# ご意見など

ご意見、間違い訂正などございましたらお寄せ下さい。