感動した！ min-max max-min の最適化表現 手法

# 初めに
本記事は「数理最適化 Advent Calendar 2020」の第1日目となります。

https://qiita.com/advent-calendar/2020/mo

まだ空きがありますので、基本的なことでもいいので知識をアウトプットして共有していきましょう！

# 数理最適化とは


From [WikiPedia]
[WikiPedia]: https://ja.wikipedia.org/wiki/%E6%95%B0%E7%90%86%E6%9C%80%E9%81%A9%E5%8C%96
> 数学の計算機科学やオペレーションズリサーチの分野における数理最適化（すうりさいてきか、英: mathematical optimization）とは、（ある条件に関して）最もよい元を、利用可能な集合から選択することをいう。


難しいことはおいておいて、数理最適化の世界では『もっともよい解』を判断する基準 や 『ある条件』を数式で表す。

それぞれの数式を以下のように呼ぶ。

- 『もっともよい解』を判断する基準を表す数式 = 目的関数
- 『ある条件』を表す数式 = 制約式

# 目的関数について

目的関数とは、最適にしたいことを表す式です。定義的には何でもOK何ですが、解を得る、特にコンピュータで解くには2次以上になるとOSSで解を得ることは難しいため、頑張って1次の式で表現します。

さて、目的関数ですが色々な実問題を解いていると、いくつかの課題の最小値を最大化したい（または最大値を最小化したい）というケースが出てきます。疑似的に数式で表すと、

$$ \max  {\min f(x)}$$

のような形です。普通に解こうと考えると、これは難しい問題になりますが、あるテクニックを使うときれいに解け、この手法に感動したため本解説を作成しました。

結論から言うと、別の変数$z$を定義し、このzを最大化することを目的と置きます。
また、同時に$z \le f(x)$の制約式を加えます。これにより、変数$z$が各$f(x)$の最小値であることが保証されます。

これだけではわかりづらいので、次に一例をあげます。

# ある数学の未解決問題

1m×1mのエリア内になるべく距離を距離を離して複数の円を置くとき、どのような配置にすればよいでしょうか？
古くからある課題で、わたしは秋山 仁先生が何かの雑誌に簡単な記事で知りました。（たしか雑誌では花壇の中に花を植えるみたいな感じだった）
この問題ですが、問題定義が有名なためか、数学の問題としてはかなりポピュラーで「点がXこの時の最適解はXXX」という形でそれぞれの点の個数ごとに証明がされているようです。
類似問題だけを取り扱ったマニアックなページがいかにあります。

http://www.packomania.com/

さて、この問題の目的関数は、「各点間の最小距離を最大化すること」です。すなわち、2点間の距離を算出する$f(x_i,x_j)$があるとして、目的関数は先ほど出た

$$max(min f(x_i,x_j))$$

の形となります。


# 線形計画に落とす

2点の間の距離を出すためには、二つの点$P_1,P_2$の座標をそれぞれ$(x_1,y_1),(x_2,y_2)$とすると

$$ \sqrt{(x_1-x_2)^2+(y_1-y_2)^2} $$

となります。距離を扱う以上、普通に計算すると2次の最適化となり、解を得ることが難しくなります。
これをどうにかするのは今回の本旨ではないので、点の位置をグリッドのどこかとなるように決めます。つまり、問題を離散化して組み合わせ最適化とすれば、距離は事前に計算することができ、線形の最適化問題として表現できるようになります。
※別の難しさは発生しますが

例えば、10×10のマス目に2点を配置する場合は対角に配置するのがよさそうです。
では、3点ではどうでしょう？おそらく、1点は角に配置され、残り2点は正三角形を構成するように配置するのがBestではないでしょうか？

| 2点の配置例 | 3点の配置例 |
|:-:|:-:|
| ![twopoints.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/98c9e330-6ed9-96d0-db75-c8050dd1065d.png) |  ![threepoints.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/5892c448-b1d7-c069-cfe4-4edabc3f92c4.png)|

# 定式化
定式化のまとめから記載すると、目的関数と制約式は以下となります。
目的関数： $z$ の最大化
制約式：
$ 制約式 (1) : \sum x_{(i,j)} = N $
$ 制約式 (2) : z \le  (1-x_{(i,j)}) * M + (1-x_{(k,l)}) * M  + dist((i,j),(k,l))$
$(ただし、dist((i,j),(k,l)) は点(i,j)と、点(k,l)の間の距離、Mは十分大きな距離)$

## 目的関数
配置した点間の最短距離をzと置くと、目的関数はzの最大化となります。

## 制約式
配置する点の数を$N$とするとします。また、10×10に区切られたマスのうち、点を配置するマスであれば、1、そうでなければ0となる変数$x_{(i,j)}$を定義します。

この時、配置する点が$N$のため、すべての変数$x_{(i,j)}$を合計するとNとなります。
すなわち

$ 制約式 : (1) \sum x_{(i,j)} = N $

次に配置した距離間がzより大きくなるという制約式を加えます。これがキモです。
$z \le 点間の距離 (ただし、すべての点の組み合わせにおいて)$
です。

そうすると、$z$はすべての点間の距離の最小値となりますが、目的関数が$z$の最大化のため、$max(min f(x_i,x_j))$を実現することができます。

$x_{(i,j)}$を用いてどう表現するか悩みましたが、以下で表現できました。

$ 制約式 (2) : z \le  (1-x_{(i,j)}) * M + (1-x_{(k,l)}) * M  + dist((i,j),(k,l))$

ここで、distは点$(i,j)$と、点$(k,l)$の間の距離、Mは十分大きな距離とします。（1×1の対角の2点を配置した場合、とりうる距離が最大$\sqrt2$となります。その、$\sqrt2$より大きい値とします）

この時、点$(i,j)$と、点$(k,l)$の両方が選ばれた場合（両方が1の場合）、上記制約式は$z \le dist((i,j),(k,l))$となり、意味を持ちます。
点$(i,j)$と、点$(k,l)$の少なくとも一方が0の場合は、$z \le 大きな数$となり、常に満たされる不等式となります。

#実装
Zは実数、xは離散の変数（0，1）となるので、MIP(Mixed Integer Linear Problem）となります。

```python:離散点配置問題
from pulp import LpVariable,LpProblem,LpBinary,LpContinuous,LpMaximize,LpStatus,value,PULP_CBC_CMD
from math import sqrt
from itertools import product

######定数定義
#格子の数
SIZE = 10
#ばらまくPointの数
N = 5

#2点間の距離を算出する
def dist(p1,p2):
    return(sqrt((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)/SIZE)

######変数定義
#どこにPointを置くか
x = LpVariable.dict(name="x",indexs=(range(SIZE),range(SIZE)),lowBound=0,upBound=1,cat=LpBinary)
#今回のキモ 多目的のための変数(どんなに良くても対角の√2以下)
z = LpVariable("z",lowBound=0,upBound=1.5,cat=LpContinuous)

#問題定義
p = LpProblem(name="SpreadingPointsProblem",sense=LpMaximize)

#zの最大化問題として定義
p += z

#N個のPointが配置される
tmp_sub = 0
for i,j in product(range(SIZE),range(SIZE)):
    tmp_sub += x[(i,j)]
p += tmp_sub == N

#確認した2点に同時に点が配置されているとしたばあい、その距離はz以上
BIGM = 100.0
for i,j in product(range(SIZE),range(SIZE)):
    for k,l in product(range(SIZE),range(SIZE)):
        if i==k and j==l:
            continue
        #x[(i,j)] x[(k,l)]が両方1のときのみ意味ある不等式となる。
        p += z <= (1-x[(i,j)])*BIGM + (1-x[(k,l)])*BIGM + dist((i,j),(k,l))

solver = PULP_CBC_CMD(msg=0,mip=True,threads=15)
p.solve(solver)
print(LpStatus[p.status])
print(f"Best Distance = {value(p.objective)}")

#結果出力
for j in range(SIZE):
    print(",".join(str(int(x[(j,i)].value())) for i in range(SIZE)))
```


``` shell:5点での実行結果
Optimal
Best Distance = 0.56568542
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,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
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,1

```

| 5点の最適配置 | 6点の最適配置 |7点の最適配置 |
|:-:|:-:|:-:|
| ![fivepoints.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/bd4d9a1d-8cb5-66c7-4a2f-c9c8c9cddcbd.png) | ![sixpoints.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/769ac3a8-1468-1c5b-3939-2758c5f1a5cf.png)  |![sevenpoints.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/5989851b-79f0-db02-a169-b65d701b6ac0.png)|

#最後に

Advent Calendar2日目はtaka_horibeさんの「1から主双対内点法を理解する話」です！
ありがたい！

https://qiita.com/advent-calendar/2020/mo