容量制限なし施設配置問題

# 数理最適化とは

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）とは、（ある条件に関して）最もよい元を、利用可能な集合から選択することをいう。


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

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

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

----------------------

# 問題定義

> 顧客の集合(需要地点)$D$と施設の配置可能地点の集合$F$が与えられ、各施設$i \in F$を開設するための固定コストを$f_i$、各顧客$j \in D$と核施設$i$の間に1単位の需要が移動するときにかかる輸送コストを$c_{ij}$とする。このとき、施設開設固定コストと輸送コストの総和が最小となるように開設する施設を選択せよ

<div align="center">
<img src=https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/121950/95141788-640f-6b38-9dd5-6065f0d6693c.png alt="属性" title="問題イメージ">
</div>

----------------------

# 定式化
## 目的関数

目的はコストの最小化。この問題ではコストは需要の移動にかかるコスト、施設の開設コストの2種があり、その和となる。

$$ 目的関数(最小化)： \sum_{i \in F} f_i y_i + \sum_{i \in F} \sum_{j \in D} c_{ij} x_{ij}$$

### 需要の移動にかかるコスト
施設$i$で満たす顧客$j$の需要の割合を$x_{ij}$とする。
この時、輸送にかかるコストは$c_{ij} x_{ij}$であり、すべての施設と顧客の組み合わせを考えると
$$ \sum_{i \in F} \sum_{j \in D} c_{ij} x_{ij} $$

### 施設の開設にかかるコスト
施設$i$を開設するかどうかを$y_i$で表すと、施設開設にかかる総コストは
$$ \sum_{i \in F} f_i y_i $$

----------------------

# 制約式
開設していない施設($y_i=0$)からは需要を満たせない、開設した施設($y_i=1$)からは需要を満たす(割合なので最大 $x_{ij}=1$)ことができるので、
$$ x_{ij} \le y_i \hspace{50pt} \forall i \in F, \forall j \in D$$

顧客$j$ごとに、需要の割合を足して1になる必要があるため、
$$ \sum_{i \in F} x_{ij} = 1 \hspace{50pt} \forall j \in D $$

$x_{ij}$は割合のため、
$$  x_{ij} \ge 0 \hspace{50pt} \forall i \in F, \forall j \in D$$

$y_i$は施設を作る場合1,作らない場合0の2値なので

$$ y_i \in \\{ 0,1 \\} $$

----------------------
#実装

```python:Pulpを用いた実装
from pulp import LpProblem,lpSum,LpBinary,LpContinuous,LpVariable,LpMinimize
import random
import numpy as np
from itertools import product

random.seed(1)

#顧客の集合(顧客はjで表す)
D_SIZE = 30

#施設候補場所(施設はiで表す)
F_SIZE = 5

###問題設定
#施設建設にかかるコスト
f = np.random.randint(low=100,high=200,size=F_SIZE)
#各施設から顧客への輸送コスト
c = np.random.randint(low=10,high=100,size=(F_SIZE,D_SIZE))

p = LpProblem(name="FLP",sense=LpMinimize)

#施設iから顧客jに運ぶ製品の割合（どれだけ需要を満たすかの割合）
x = LpVariable.dict('x',indexs=(range(F_SIZE),range(D_SIZE)),lowBound=0,upBound=1,cat=LpContinuous)

#施設iを開設するか(1)しないか(0)を表す変数y
y = LpVariable.dict('y',indexs=(range(F_SIZE)),lowBound=0,upBound=1,cat=LpBinary)

#目的関数定義
p += lpSum([f[i]*y[i] for i in range(F_SIZE)]) + lpSum([c[(i,j)]*x[(i,j)] for i,j in product(range(F_SIZE),range(D_SIZE))])

#施設のありなしによって、製品をその施設から出荷できるかどうか
for i,j in product(range(F_SIZE),range(D_SIZE)):
    p += x[(i,j)] <= y[i]

#割合の合計は1
for j in range(D_SIZE):
    p += lpSum([x[(i,j)] for i in range(F_SIZE)]) == 1.0

#解く
p.solve()

#最適化できたら解を出力
if p.status == 1:
    print(", ".join(f"y[{i}]={int(y[i].value())}" for i in range(F_SIZE)))
    for j in range(D_SIZE):
        print(", ".join(f"x[{i:02},{j:02}]={x[(i,j)].value()}" for i in range(F_SIZE)))
```
----------------------
#実行結果

最適化処理は0.01秒と出ています。
結果を見て気が付きましたが、もし施設を建てることになるとコストが最も安い経路ですべての需要を満たすことになります。
もう少し供給量や輸送量などの上限がないと、あまり最適化問題としての面白みがない気がします。

```shell:出力
Welcome to the CBC MILP Solver 
Version: 2.9.0 
Build Date: Feb 12 2015 

command line - XXXX
/tmp/304c5513b4f14eaaa8719c3635811fc8-pulp.mps ratio None allow None threads None presolve on strong None gomory on knapsack on probing on branch printingOptions all solution /tmp/304c5513b4f14eaaa8719c3635811fc8-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 185 COLUMNS
At line 801 RHS
At line 982 BOUNDS
At line 1138 ENDATA
Problem MODEL has 180 rows, 155 columns and 450 elements
Coin0008I MODEL read with 0 errors
String of None is illegal for double parameter ratioGap value remains 0
String of None is illegal for double parameter allowableGap value remains 0
String of None is illegal for integer parameter threads value remains 0
String of None is illegal for integer parameter strongBranching value remains 5
Option for gomoryCuts changed from ifmove to on
Option for knapsackCuts changed from ifmove to on
Continuous objective value is 1275 - 0.00 seconds
Cgl0004I processed model has 180 rows, 155 columns (5 integer (5 of which binary)) and 450 elements
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 1275
Cbc0038I Relaxing continuous gives 1275
Cbc0038I Before mini branch and bound, 5 integers at bound fixed and 150 continuous
Cbc0038I Mini branch and bound did not improve solution (0.00 seconds)
Cbc0038I After 0.00 seconds - Feasibility pump exiting with objective of 1275 - took 0.00 seconds
Cbc0012I Integer solution of 1275 found by feasibility pump after 0 iterations and 0 nodes (0.00 seconds)
Cbc0001I Search completed - best objective 1275, took 0 iterations and 0 nodes (0.00 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 1275 to 1275
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                1275.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.01   (Wallclock seconds):       0.01

y[0]=1, y[1]=1, y[2]=1, y[3]=0, y[4]=0
x[00,00]=0.0, x[01,00]=0.0, x[02,00]=1.0, x[03,00]=0.0, x[04,00]=0.0
x[00,01]=0.0, x[01,01]=1.0, x[02,01]=0.0, x[03,01]=0.0, x[04,01]=0.0
x[00,02]=0.0, x[01,02]=1.0, x[02,02]=0.0, x[03,02]=0.0, x[04,02]=0.0
x[00,03]=0.0, x[01,03]=0.0, x[02,03]=1.0, x[03,03]=0.0, x[04,03]=0.0
x[00,04]=0.0, x[01,04]=0.0, x[02,04]=1.0, x[03,04]=0.0, x[04,04]=0.0
x[00,05]=0.0, x[01,05]=0.0, x[02,05]=1.0, x[03,05]=0.0, x[04,05]=0.0
x[00,06]=1.0, x[01,06]=0.0, x[02,06]=0.0, x[03,06]=0.0, x[04,06]=0.0
x[00,07]=1.0, x[01,07]=0.0, x[02,07]=0.0, x[03,07]=0.0, x[04,07]=0.0
x[00,08]=0.0, x[01,08]=1.0, x[02,08]=0.0, x[03,08]=0.0, x[04,08]=0.0
x[00,09]=0.0, x[01,09]=0.0, x[02,09]=1.0, x[03,09]=0.0, x[04,09]=0.0
x[00,10]=1.0, x[01,10]=0.0, x[02,10]=0.0, x[03,10]=0.0, x[04,10]=0.0
x[00,11]=1.0, x[01,11]=0.0, x[02,11]=0.0, x[03,11]=0.0, x[04,11]=0.0
x[00,12]=1.0, x[01,12]=0.0, x[02,12]=0.0, x[03,12]=0.0, x[04,12]=0.0
x[00,13]=0.0, x[01,13]=1.0, x[02,13]=0.0, x[03,13]=0.0, x[04,13]=0.0
x[00,14]=1.0, x[01,14]=0.0, x[02,14]=0.0, x[03,14]=0.0, x[04,14]=0.0
x[00,15]=0.0, x[01,15]=1.0, x[02,15]=0.0, x[03,15]=0.0, x[04,15]=0.0
x[00,16]=1.0, x[01,16]=0.0, x[02,16]=0.0, x[03,16]=0.0, x[04,16]=0.0
x[00,17]=1.0, x[01,17]=0.0, x[02,17]=0.0, x[03,17]=0.0, x[04,17]=0.0
x[00,18]=1.0, x[01,18]=0.0, x[02,18]=0.0, x[03,18]=0.0, x[04,18]=0.0
x[00,19]=0.0, x[01,19]=0.0, x[02,19]=1.0, x[03,19]=0.0, x[04,19]=0.0
x[00,20]=0.0, x[01,20]=0.0, x[02,20]=1.0, x[03,20]=0.0, x[04,20]=0.0
x[00,21]=0.0, x[01,21]=0.0, x[02,21]=1.0, x[03,21]=0.0, x[04,21]=0.0
x[00,22]=1.0, x[01,22]=0.0, x[02,22]=0.0, x[03,22]=0.0, x[04,22]=0.0
x[00,23]=0.0, x[01,23]=1.0, x[02,23]=0.0, x[03,23]=0.0, x[04,23]=0.0
x[00,24]=0.0, x[01,24]=1.0, x[02,24]=0.0, x[03,24]=0.0, x[04,24]=0.0
x[00,25]=0.0, x[01,25]=1.0, x[02,25]=0.0, x[03,25]=0.0, x[04,25]=0.0
x[00,26]=0.0, x[01,26]=1.0, x[02,26]=0.0, x[03,26]=0.0, x[04,26]=0.0
x[00,27]=1.0, x[01,27]=0.0, x[02,27]=0.0, x[03,27]=0.0, x[04,27]=0.0
x[00,28]=0.0, x[01,28]=1.0, x[02,28]=0.0, x[03,28]=0.0, x[04,28]=0.0
x[00,29]=0.0, x[01,29]=1.0, x[02,29]=0.0, x[03,29]=0.0, x[04,29]=0.0
```

----------------------

# おまけ

施設配置問題には色々なバリエーションがあるが、多くは大規模な整数計画問題として定式化されNP困難なため厳密解を求めるのは難しい。

## バリエーション例

- 移動距離の総和を最小にするミニサム問題
- 集合被覆問題として解く、カバリング型の問題など

## 適用分野

- 工場、倉庫、配送センター、新店舗などの立地選択や統廃合
- 消防署、病院などの緊急施設、廃棄物処理場施設の立地選択