CPLEXサンプル warehouseを読み解く

CPLEXのサンプルのwarehouseを読み解きます。


<サンプル導入ディレクトリ>opl\examples\opl\warehouse\warehouse.mod
にあります。

* 実行環境
    * CPLEX 22.1.1
    * Windows 11 64bit

# 問題

既存の店舗に供給するための倉庫を開設する場所をいくつか検討しています。それぞれの倉庫には、開発固定費と対応できる店舗数を示す最大容量があります。各店舗は1つの倉庫からしか供給を受けず、店舗への供給コストは選択した倉庫によって異なります。最もコストを下げられるように、どの倉庫を開設し、どの店舗にどの倉庫から供給するかを決定します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/32a62287-4ff0-c1b3-044d-9fc6639d1791.png)

## 利用データ
1.   一つの倉庫の開設に必要な固定費: 30
2.  各倉庫から供給できる店舗数上限

|倉庫|ボン|ボルドー|ロンドン|パリ|ローマ|
|----:|----:|----:|----:|----:|----:|
|容量|1 |4 |2 |1 |3 |


3.各倉庫から各店舗への輸送費

|倉庫|ボン|ボルドー|ロンドン|パリ|ローマ|
|----|----:|----:|----:|----:|----:|
|店舗 1 |20 |24 |11 |25 |30 |
|店舗 2 |28 |27 |82 |83 |74 |
|店舗 3 |74 |97 |71 |96 |70 |
|店舗 4 |2 |55 |73 |69 |61 |
|店舗 5 |46 |96 |59 |83 |4 |
|店舗 6 |42 |22 |29 |67 |59 |
|店舗 7 |1 |5 |73 |59 |56 |
|店舗 8 |10 |73 |13 |43 |96 |
|店舗 9 |93 |35 |63 |85 |46 |
|店舗 10 |47 |65 |55 |71 |95|

## 目的関数
「倉庫開設固定費＋各倉庫から各店舗への輸送費＝総コスト」を最小化

## 決定変数
0.  各倉庫を開設するかどうか
0.  どの倉庫からどの店舗に供給するか

## 制約
0.  各店舗は一つの倉庫から必ず供給を受ける
0.  開設した倉庫からのみ供給する
0.  各倉庫から供給できる店舗数には上限がある

## 問題の種類
整数計画

# 解
先に解を紹介します。

## 目的関数
合計費用：383

## 決定変数
0.  各倉庫を開設するかどうか：青いセルで表示した倉庫
0.  どの倉庫からどの店舗に供給するか:緑のセルの倉庫と店舗の組み合わせ

この解で面白いのは以下のような点になります。
- 「ボン」は「店舗7」への費用が1で最低だが、「店舗4」を優先した方が合計では低くなる
- 「店舗1」は「ローマ」よりも「パリ」から供給した方がコストは安いが、「パリ」の倉庫の開設費用を抑えた方がトータル費用は抑えることができる

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/18ae54a7-1d45-84a2-4020-21c7e926e43b.png)

# OPLコード解説
OPLコードの解説を行っていきます。

## 変数定義
まず配列やレンジの変数の定義を行っています。
```opl:配列定義.mod
{string} Warehouses = ...;
int NbStores = ...;
range Stores = 0..NbStores-1;
```
これらの定義はwarehouse.datファイルの中で行っています。

```dat:配列定義.dat
Warehouses = {Bonn,Bordeaux,London,Paris,Rome};
NbStores = 10;
```

## 利用データ
次に利用データの定義を行っています。
各店舗と倉庫間の輸送費は2次元の配列で定義されます。

```opl:利用データ.mod
//開設固定費
int Fixed = ...;
//各倉庫の許容量
int Capacity[Warehouses] = ...;
//各店舗と倉庫間の輸送費
int SupplyCost[Stores][Warehouses] = ...;
```
これらの定義はwarehouse.datファイルの中で行っています。

```dat:利用データ.dat
Fixed = 30;
Capacity = [1,4,2,1,3];
SupplyCost = [ 
       [ 20, 24, 11, 25, 30 ], 
       [ 28, 27, 82, 83, 74 ],
       [ 74, 97, 71, 96, 70 ],
       [  2, 55, 73, 69, 61 ],
       [ 46, 96, 59, 83,  4 ],
       [ 42, 22, 29, 67, 59 ],
       [  1,  5, 73, 59, 56 ],
       [ 10, 73, 13, 43, 96 ],
       [ 93, 35, 63, 85, 46 ],
       [ 47, 65, 55, 71, 95 ] ];
```
CapacityとSupplyCostは配列であるため、問題ブラウザで見ると以下のようになっています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/be18f1d2-bbae-2992-7917-7131ae6c572f.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/4a87289b-fc0d-cf24-5a6e-e7a461788112.png)

## 決定変数
次に決定変数の定義を行っています。
dvarが決定変数を意味します。
「倉庫の開設の有無」は1次元の配列、
「倉庫と店舗の供給の有無」は2次元の配列で定義されます。

```opl:決定変数.mod
////決定変数
//倉庫の開設の有無
dvar boolean Open[Warehouses];
//倉庫と店舗の供給の有無
dvar boolean Supply[Stores][Warehouses];
```

これらの決定変数は、求められた解では以下のようになりました。どちらもBooleanなので「0,1」のフラグで表されます。Openで5個、Supplyで50個のフラグがあり、合計で55個になります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/243535a0-4a74-b166-3177-2b32817f7613.png)


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/f9d6f8e1-3b08-d36a-b685-35cd406da188.png)




## 目的関数
次に目的関数の定義を行っています。

まず、minimizeなので最小化問題になります。

sumは総和であるΣを意味します。
sum( w in Warehouses ) Open[w] *Fixed は、開設したすべての倉庫数＊開設固定費を意味します。

sum( w in Warehouses , s in Stores ) のようにsumの中に二つの配列がある場合は、ΣΣが二つ並んでいることを意味します。
sum( w in Warehouses , s in Stores ) SupplyCost[s][w] * Supply[s][w]は「供給コスト」と「供給の有無」を総当たりで掛け算した合計ということになります。つまり供給の総コストになります。

```python:目的関数.mod
////目的関数
minimize
  //開設コスト
  sum( w in Warehouses ) Open[w] * Fixed 
  //輸送コスト
  + sum( w in Warehouses , s in Stores ) SupplyCost[s][w] * Supply[s][w];
```


数式で表すと以下です。

```math
minimize \sum_{}^{w} Open_w*Fixed+\sum_{}^{w}\sum_{}^{s}SupplyCost_{sw}*{Supply}_{sw}
```

図で表すと以下のようなイメージです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/acbd621e-9d89-6780-7f35-7c776bbd2111.png)



## 制約
次に制約の定義を行っています。
subject toの中カッコで囲みます。
```python:制約.mod
subject to{

---略----
}

```

### 制約1: 各店舗は一つの倉庫

「各店舗は一つの倉庫」の制約は以下のように定義されます。
すべての店舗で、利用する倉庫をかならず一つ割り当てるという意味です。

forall( s in Stores )で店舗数分の制約を作っています。
ctEachStoreHasOneWarehouse:で制約にラベルを付けています。
sum( w in  Warehouses ) Supply[s][w]は、ある店舗に供給するすべての倉庫数の合計をだしています。
店舗は必ず一つの倉庫から納入するので「==1」の制約になります。

```python:制約1: 各店舗は一つの倉庫.mod
  //各店舗は一つの倉庫
  forall( s in Stores )
    ctEachStoreHasOneWarehouse:
      sum( w in  Warehouses ) 
        Supply[s][w] == 1;
```

数式で表すと以下です。

```math
\sum_{}^{w}Supply_{sw}=1　 \{すべての店舗sに対して\}\\
```

店舗数分の10個の制約が作られています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/27901abe-f171-8992-fde8-eeedf7cd4466.png)


図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/f3786435-a614-933f-98dc-572ee24496b7.png)


### 制約2: 開設した倉庫のみを利用
「開設した倉庫のみを利用」の制約は以下のように定義されます。

forall( w in Warehouses, s in Stores )で店舗数×倉庫数分の制約を作っています。
ctUseOpenWarehouses:で制約にラベルを付けています。
供給の有無を示すSupply[s][w]には「0」か「1」が入ります。 開設するかどうかのOpen[w]にも「0」か「1」が入ります。開設していない（Open[w]=0）場合にはその倉庫からの供給は受けられないので その倉庫からすべての店舗への供給のSupply[s]はすべて0でなければなりません。倉庫が開設される（Open[w]=1）場合も、店舗が供給を受けるとは限らないので「＝」ではなく「<=」になります。

```python:制約2: 開設した倉庫のみを利用.mod
  //開設した倉庫のみを利用
  forall( w in Warehouses, s in Stores )
    ctUseOpenWarehouses:
      Supply[s][w] <= Open[w];
```

数式で表すと以下です。

```math
Supply_{sw} <= Open_w \{すべての倉庫wとすべての店舗sに対して\}\\
```

店舗数×倉庫数分で50個の制約が作られています。
「スラック」とは制約式の余裕量を示しています。開設した倉庫のフラグとSupplyのフラグの差です。開設した倉庫には１が立ち、Supplyも1になったら、その差であるスラックは0になります。倉庫Bonnと店舗3のスラックは供給されたので０になります。

また一方で、倉庫Parisは開設されなかったので、すべてのスラックは全て０です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/388107c1-c937-656b-0a64-ec90cfd6ffef.png)

図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/9fe5a7d3-6f9e-8240-225b-73bd4e0f82b3.png)



### 制約3: 倉庫許容量以下の供給
「倉庫許容量以下の供給」の制約は以下のように定義されます。各倉庫の許容量以下の供給であることを示しています。

forall( w in Warehouses )で倉庫数分の制約を作っています。
ctMaxUseOfWarehouse:で制約にラベルを付けています。

sum( s in Stores ) Supply[s][w]はある倉庫に供給する店舗すべての供給量です。これがその倉庫の許容量Capacity[w]以下である必要があります。


```python:制約3: 倉庫許容量以下の供給.mod
  //倉庫許容量以下の供給
  forall( w in Warehouses )
    ctMaxUseOfWarehouse:         
      sum( s in Stores ) 
        Supply[s][w] <= Capacity[w];
```

数式で表すと以下です。

```math
\sum_{}^{s}Supply_{sw}<=Capacity_w　 \{すべての倉庫wに対して\}\\
```

倉庫数分の5個の制約が作られています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/f003e3a8-bacb-e487-5a51-455a33c2f18d.png)

図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/8e632dd8-8350-a834-09b2-07ef72d45c84.png)

# 完成OPL

製品サンプルと同じ内容ですが、日本語コメントを入れたり順序を入れ替えています。

```Java:warehouse.mod

//配列定義
{string} Warehouses = ...;
int NbStores = ...;
range Stores = 0..NbStores-1;

////利用データ
//開設固定費
int Fixed = ...;
//各倉庫の許容量
int Capacity[Warehouses] = ...;
//各店舗と倉庫間の輸送費
int SupplyCost[Stores][Warehouses] = ...;

////決定変数
//倉庫の開設の有無
dvar boolean Open[Warehouses];
//倉庫と店舗の供給の有無
dvar boolean Supply[Stores][Warehouses];


////目的関数
minimize
  //開設コスト
  sum( w in Warehouses ) Open[w] * Fixed 
  //輸送コスト
  + sum( w in Warehouses , s in Stores ) SupplyCost[s][w] * Supply[s][w];
    
////制約
subject to{
  //各店舗は一つの倉庫
  forall( s in Stores )
    ctEachStoreHasOneWarehouse:
      sum( w in  Warehouses ) Supply[s][w] == 1;
  //開設した倉庫のみを利用
  forall( w in Warehouses, s in Stores )
    ctUseOpenWarehouses:
      Supply[s][w] <= Open[w];
  //倉庫許容量以下の供給
  forall( w in Warehouses )
    ctMaxUseOfWarehouse:         
      sum( s in Stores ) Supply[s][w] <= Capacity[w];
}

{int} Storesof[w in Warehouses] = { s | s in Stores : Supply[s][w] == 1 };
execute DISPLAY_RESULTS{
  writeln("Open=",Open);
  writeln("Storesof=",Storesof);
}
```

# おまけ

制約2と3は一緒に書いてしまう方法もあります。
Capacity[w]*Open[w]としてしまう方法です。この場合Parisのように開設しない場合はCapacityが結局0になります。こちらの方が直感的にわかりやすいかと思いました。


```Java:warehouse.mod
////制約
subject to{
  //各店舗は一つの倉庫
  forall( s in Stores )
    ctEachStoreHasOneWarehouse:
      sum( w in  Warehouses ) Supply[s][w] == 1;
  //開設した倉庫のみを利用
  //forall( w in Warehouses, s in Stores )
  // ctUseOpenWarehouses:
  //    Supply[s][w] <= Open[w];

  //開設した倉庫許容量以下の供給
  forall( w in Warehouses )
    ctMaxUseOfWarehouse:         
      sum( s in Stores ) Supply[s][w] <= Capacity[w]*Open[w];
}
```

# 参考
#### Warehouse location - IBM Documentation 

https://www.ibm.com/docs/ja/icos/20.1.0?topic=programming-warehouse-location

CPLEXサンプルを読み解く記事一覧

https://qiita.com/kawada2017/items/bdb85f66e5f1de6039d6