CPLEXサンプル steelmill（製鉄所:在庫突き合わせ問題）を読み解く

CPLEXのサンプルのsteelmill（製鉄所:在庫突き合わせ問題）を読み解きます。
この問題は、CSPLibという制約プログラミング用のテスト問題のライブラリーで提供されている[038: Steel Mill Slab Design](https://www.csplib.org/Problems/prob038/)が元ネタになっています。

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

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

# 問題

製鉄所ではスラブという半製品材料の鋼片を圧延加工して、製品である鋼板を作っていきます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/fd99d2f2-328f-caf6-c637-5969ca4a09e9.png)

#### 参考:鉄の出来るまで（製鋼から熱間圧延まで） 

https://www.chiba-muse.or.jp/SCIENCE/16jfe/tetsu2.html

スラブには大きいものから小さいものまで重さによる種類があります。また、一つのスラブから2色まで色を付けることができます。一方で鋼板の注文もそれぞれ必要な重さと色があります。どのサイズのスラブからどの注文の鋼板を切り出すと、最も無駄なくスラブをつかえるかを決定します。
製鉄所ではスラブという半製品材料の鋼片を圧延加工して、製品である鋼板を作っていきます。スラブには大きいものから軽いものまで重さによる種類があります。また、一つのスラブから2色まで色を付けることができます。
一方で鋼板の注文もそれぞれ必要な重さと色があります。どのサイズのスラブからどの注文の鋼板を切り出すと、最も無駄なくスラブをつかえるかを決定します。同じ重さのスラブを何枚使ってもかまいません。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/3353a5a6-910d-57bb-c61e-398e899f30f8.png)



## 利用データ
1.  最大の利用可能なスラブ数: 10
1.  スラブあたりの色数上限:2
2.  スラブの重さの種類

|スラブの重さの種類|
|----:|
|0|
|10|
|15|
|25|



4.各注文の必要量と色番号


|注文番号|重さ|色番号|
|----:|----:|----|
|1|22|5緑|
|2|9|3黄|
|3|9|4黄緑|
|4|8|5緑|
|5|8|7青|
|6|6|3黄|
|7|5|6青緑|
|8|3|0ピンク|
|9|3|2オレンジ|
|10|3|3黄|
|11|3|1赤|
|12|3|5緑|


## 決定変数
0.  各注文をどのスラブに割り当てるか
0.  使用する各スラブの容量


## 目的関数
「使用する各容量のスラブから注文分する部分を除いた」ロス分の合計量


## 制約
0.  すべての注文が各スラブに割り当てられる
0.  スラブあたり色数は2色以内

## 問題の種類
制約プログラミング

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

## 目的関数
ロスの合計:3

## 決定変数
0.  どの注文をどのスラブに割り当てるか：色のついたセル
0.  各スラブの使用量

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/92a96aa9-1bd4-f6ca-3197-b5c923a36c6e.png)



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

## 利用データ
まずパラメーターの定義を行っています。
```java:パラメーター.mod
/// パラメーター
tuple TParams {
  int maxNbSlabs; //最大の利用可能なスラブ数
  int maxColorsPerSlab; //スラブあたりの色数上限 
}

/// 単一行のテーブル形式.
TParams Params = ...;
```
これらの定義はsteelmill.datファイルの中で行っています。

```java:パラメーター.dat
Params =
<
  10, // 最大の利用可能なスラブ数
   2   // スラブあたりの色数上限
>;
```


次にスラブ関係のデータの定義を行っています。

```java:スラブデータ.mod
// 利用可能なスラブの総数。 
// 理論的にはこれに制限はありませんが、実用的な最適化モデルを作成するために、合理的な上限を課します。
int nbSlabs = Params.maxNbSlabs;
range Slabs = 1..nbSlabs;

// スラブの重さの種類
{int} slabWeights = ...;
// 最も重いスラブ
int maxSlabWeight = max (sw in slabWeights) sw;
```
スラブの重さの種類はsteelmill.datファイルの中で行っています。
サイズがないスラブも定義しています。スラブは最大10個まで利用可能ですが、あまったスラブはサイズ0で定義したいためです。

```java:スラブデータ.dat
// 利用可能なスラブの総数
// サイズ 0 のスラブはスラブがないことを表します
slabWeights = { 0, 10, 15, 25 };
```

maxSlabWeightはslabWeightsの最大値ですので、問題ブラウザで見ると以下のようになっています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/ac69b5ad-535d-456a-f78b-a01663cc7504.png)

次に注文関係のデータの定義を行っています。

```java:注文データ.mod
// 注文コイルの重さと色
tuple TOrder {
   key int index;
   int weight;
   int color;
}
{TOrder} Orders = ...;

int nbOrders = card(Orders);
range orderIndexRange = 1..nbOrders;
```
スラブの重さの種類はsteelmill.datファイルの中で行っています。

```java:注文データ.dat
Orders =
{
    <1, 22, 5>
  , <2,  9, 3>
  , <3,  9, 4>
  , <4,  8, 5>
  , <5,  8, 7>
  , <6,  6, 3>
  , <7,  5, 6>
  , <8,  3, 0>
  , <9,  3, 2>
  , <10, 3, 3>
  , <11, 3, 1>
  , <12, 3, 5>
};

```

Ordersを問題ブラウザで見ると以下のようになっています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/18bc6894-3fb7-9920-87e5-77a2aed6bdcf.png)

[card関数](https://www.ibm.com/docs/ja/icos/20.1.0?topic=functions-card)はセット内の項目数を返します。ですので、nbOrdersは以下のように12になります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/b30beed4-5ab1-c948-78f7-d4a3b18b849f.png)

[unionキーワード](https://www.ibm.com/docs/ja/icos/20.1.0?topic=keywords-union)は重複を取り除いてリストを作ります。

```java:注文の中で使われている全種類の色
{int} allcolors = union(o in Orders) { o.color};
```

Ordersからcolorの列だけを取り出して重複を取り除いたリストを作っています

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/05967ea3-b182-47db-0df3-dfe775f7dc39.png)

[item関数](https://www.ibm.com/docs/ja/icos/20.1.0?topic=functions-card)はセット内のタプルへのアクセスを行います。
```java:各注文の重さの配列
// CPはパック制約のために配列が必要です
int orderWeights[ i in orderIndexRange ] = item(Orders, i-1).weight;
```
Ordersからweightの列だけを取り出してリストを作っています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/e71d1d57-de53-cde9-6de2-0cb3414047ae.png)


## 決定変数
次に決定変数の定義を行っています。
dvarが決定変数を意味します。
「各注文をどのスラブに割り当てるか」も
「使用する各スラブの容量」は1次元の配列で定義されます。

```java:決定変数
// 各注文をどのスラブに割り当てるか
dvar int productionSlab[Orders] in Slabs;

// 使用する各スラブの容量
dvar int slabUse[Slabs] in 0..maxSlabWeight;
```

これらの決定変数は、求められた解では以下のようになりました。productionSlabで注文数分で12個あり、値はスラブ数分の1から10の値をとります。slabUseはスラブ数分で10個あり、値は0から25（最大のスラブ容量）までの値をとります。決定変数の合計数で22個になります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/796f4bf5-79af-f061-0592-9f2754429f22.png)




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

関数の定義に先立って、「必要量に対するロス量のリスト」lossを作ります。

[minキーワード](https://www.ibm.com/docs/ja/icos/20.1.0?topic=keywords-min)は関連式の集合の最小値を計算するキーワードです。

```java:必要量に対するロス量のリスト
// 損失は、注文量を製造するのに必要な量に対して、用意できるスラブの容量が大きい場合に発生します
int loss[use in 0..maxSlabWeight] =
  min (sw in slabWeights : sw >= use) (sw - use); 
```
lossは0から25までのリストです。0トンから25トンまでの注文量を作った場合に発生するロスを計算しています。例えば1トンの注文量を作ることができる一番小さいスラブは10トンです。この時のロスは10トン-1トン＝9トン発生します。このようにして0トンから25トンまでの注文量に対するロスをリスト化しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/c9b9fcb0-5b9d-264b-be50-f48f09aa0d13.png)


次に合計損失を計算します。

[sumキーワード](https://www.ibm.com/docs/ja/icos/20.1.0?topic=keywords-sum)は関連式の集合の合計値を計算するキーワードです。


```java:合計損失
dexpr int totalLoss = sum(s in Slabs) loss[slabUse[s]];

/////中略//////

//全てのスラブでのロスの合計を最小化する
minimize totalLoss ;
```
slabUse[1-10]からスラブの使用量を取り出し、それをlossリストに入れて、ロス量を取り出しています。
例えばSlab2は14トンですので、loss[14]で1トンのロスが出ます。
これらの合計値を算出してこれを最小化することを目的関数としています。


図で表すと以下のようなイメージです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/c33b3b17-d46d-57be-8a9d-f2539e144123.png)


## 探索フェーズ指定
制約プログラミング(CP)では[探索フェーズ](https://www.ibm.com/docs/ja/icos/20.1.0?topic=phases-what-is-search-phase)をガイドできます。ここでは決定変数はslabUseとproductionSlabがありますが、
「setSearchPhases」をつかって、productionSlabを優先して決定する決定変数としています。
問題の性質を人間が考えれば、productionSlabが決まれば自動的にslabUseが決まるので、productionSlabを優先して決めた方がいいと考えるのですが、ソルバーには業務イメージはありません。このようなヒントを与えることでより効率的な求解ができることがあります。これは必須ではないのでこの部分をコメントアウトしても動作はします。

```java:探索フェーズの定義
execute {
  var f = cp.factory;
  cp.setSearchPhases(f.searchPhase(productionSlab));
}
```

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

---略----
}

```

### 制約1: すべての注文は容量内のスラブに割り当てられる

[pack関数](https://www.ibm.com/docs/ja/icos/20.1.0?topic=functions-pack)はコンテナーのセットの積載量を維持する制約です。

CPオプティマイザー関数の[pack](https://www.ibm.com/docs/ja/icos/20.1.0?topic=f-pack)の説明の方がわかりやすかったです。
各コンテナーの容量が順守されるようにコンテナーに全積み荷を割り当てる制約です。

Pack梱包制約
(dvar int[ ]：load コンテナーの容量 ←各スラブの利用量
,dvar int[ ]:where 積み荷のコンテナの場所　←各注文に対する各スラブの割り当て
,int[] :size 積み荷のサイズ)　←各注文の必要量



```java:制約1: すべての注文は容量内のスラブに割り当てられる
  packCt:
    // すべての注文は容量内のスラブに割り当てられる
    pack(slabUse, productionSlab, orderWeights);

```



図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/d112836d-f852-b828-1b6b-2d6d771cb554.png)



### 制約2: スラブあたり色数はParams.maxColorsPerSlab以内
「開設した倉庫のみを利用」の制約は以下のように定義されます。

forall (s in Slabs)でスラブ数分の制約を作っています。
colorCt:で制約にラベルを付けています。


```java:制約2:スラブあたり色数はParams.maxColorsPerSlab以内
  forall (s in Slabs) {
    colorCt:
      // スラブあたり色数はParams.maxColorsPerSlab以内
      sum (c in allcolors) (
        or (o in Orders : o.color == c) (productionSlab[o] == s)
      ) <= Params.maxColorsPerSlab; 
  }
```
例として、Slab3で考えてみます。
o.color == 1赤の注文は1,8です。注文1のcolor=5青、注文8のcolor=0ピンクですから、どちらも1赤ではないので、OR（FALSE,FALSE)です。
一方でo.color == 5緑の注文は1,4,12です。注文1のcolor=5青ですが、注文4と注文12のcolor= 5緑ですから、OR(FALSE,TRUE,TRUE)なので、TRUEです。

このようにすべて見ていくと5青と2オレンジはTRUEです。これをsum()すると2です。つまり2色という意味です。<= Params.maxColorsPerSlab(=2)ですからすべてのスラブについて2色以下になるという制約になっています。

or (o in Orders : o.color == 0) (productionSlab[o] == 3)==FALSE
or (o in Orders : o.color == 1) (productionSlab[o] == 3)==FALSE
or (o in Orders : o.color == 2) (productionSlab[o] == 3)==TRUE
or (o in Orders : o.color == 3) (productionSlab[o] == 3)==FALSE
or (o in Orders : o.color == 4) (productionSlab[o] == 3)==FALSE
or (o in Orders : o.color == 5) (productionSlab[o] == 3)==TRUE
or (o in Orders : o.color == 6) (productionSlab[o] == 3)==FALSE
or (o in Orders : o.color == 7) (productionSlab[o] == 3)==FALSE




図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/29cc0a14-bc7b-8253-a61a-a7a5fc31b0ca.png)

## 結果出力
ここまででモデルとしては完成していますが、決定変数をわかりやすく加工して表示します。

```java:スラブ毎の注文番号
{int} fromSlab[s in Slabs] = 
{ o.index | o in Orders : productionSlab[o] == s };
```
内包表記になっています。フィルター「productionSlab[o] == s」の条件にあうOrdersのindexを抜き出しています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/0f2db4cd-ca07-b11a-7548-09e228899f6c.png)

```java:スラブ毎の色種
{int} slabColors[s in Slabs] = { o.color | o in Orders : o.index in fromSlab[s] };
```

内包表記になっています。フィルター「o.index in fromSlab[s]」の条件にあうOrdersのcolorを抜き出しています。


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/24a37f21-07f7-884a-92d3-09f42dae30fe.png)


最後に[execute](https://www.ibm.com/docs/ja/icos/20.1.0?topic=keywords-execute)でILOG Scriptを実行し、[write](https://www.ibm.com/docs/ja/icos/20.1.0?topic=functions-write)でスクリプトログに書き出しています。
 if (Opl.card(fromSlab[s]) > 0) でスラブの容量が0より大きいものだけを出力しています。つまりSlab1-5が表示されます。
slabColors[s]で書き出したcolorsは「colors = {5 0}」となり中カッコが残っています。、
「for (o in fromSlab[s]) { write(" " + o);}」で書き出したOrdersは「Orders = 1 8」とO中カッコがとれて表示されています。

```java:結果出力
execute {
  for (s in Slabs) {
    if (Opl.card(fromSlab[s]) > 0) {
      write("Slab " + s + ", Loss = " + loss[slabUse[s]]
          + ", colors =" + slabColors[s] + ", Orders =");
      for (o in fromSlab[s]) {
        write(" " + o);
      }
      writeln();
    }
  }
}
```
COSでみると以下のように表示されています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/440c2fd7-8a69-15c7-3063-fa246512390b.png)


# 完成OPL

製品サンプルとほぼ同じ内容ですが、日本語コメントを入れたり順序を入れ替えたり、添え字を入れ替えたりして読みやすくしています。

```Java:steelmilljp.mod

using CP;

/// パラメーター
tuple TParams {
  int maxNbSlabs; //最大の利用可能なスラブ数
  int maxColorsPerSlab; //スラブあたりの色数上限 
}

/// 単一行のテーブル形式.
TParams Params = ...;

// 利用可能なスラブの総数。 
// 理論的にはこれに制限はありませんが、実用的な最適化モデルを作成するために、合理的な上限を課します。
int nbSlabs = Params.maxNbSlabs;
range Slabs = 1..nbSlabs;

// スラブの重さの種類
{int} slabWeights = ...;
// 最も重いスラブ
int maxSlabWeight = max (sw in slabWeights) sw;

// 注文コイルの重さと色
tuple TOrder {
   key int index;
   int weight;
   int color;
}
{TOrder} Orders = ...;

int nbOrders = card(Orders);
range orderIndexRange = 1..nbOrders;

//注文の中で使われている全種類の色
{int} allcolors = union(o in Orders) { o.color};

// 各注文の重さの配列
// CPはパック制約のために配列が必要です
int orderWeights[ i in orderIndexRange ] = item(Orders, i-1).weight;


////決定変数
// 各注文をどのスラブに割り当てるか
dvar int productionSlab[Orders] in Slabs;

// 使用する各スラブの容量
dvar int slabUse[Slabs] in 0..maxSlabWeight;

//　必要量に対するロス量のリスト
// 損失は、注文量を製造するのに必要な量に対して、用意できるスラブの容量が大きい場合に発生します
int loss[use in 0..maxSlabWeight] =
  min (sw in slabWeights : sw >= use) (sw - use); 

// 合計損失
dexpr int totalLoss = sum(s in Slabs) loss[slabUse[s]];
//dexpr int totalUsed = sum(s in Slabs) slabUse[s];

// 探索フェーズの定義
execute {
  var f = cp.factory;
  cp.setSearchPhases(f.searchPhase(productionSlab));
}

//全てのスラブでのロスの合計を最小化する
minimize totalLoss ;

subject to {
  packCt:
    // すべての注文は容量内のスラブに割り当てられる
    pack(slabUse, productionSlab, orderWeights);

  forall (s in Slabs) {
    colorCt:
      // スラブあたり色数はParams.maxColorsPerSlab以内
      sum (c in allcolors) (
        or (o in Orders : o.color == c) (productionSlab[o] == s)
      ) <= Params.maxColorsPerSlab; 
  }
}

//// 結果出力
//　スラブ毎の注文番号
{int} fromSlab[s in Slabs] = 
{ o.index | o in Orders : productionSlab[o] == s };

//　スラブ毎の色種
{int} slabColors[s in Slabs] = { o.color | o in Orders : o.index in fromSlab[s] };

// 結果出力
execute {
  for (s in Slabs) {
    if (Opl.card(fromSlab[s]) > 0) {
      write("Slab " + s + ", Loss = " + loss[slabUse[s]]
          + ", colors =" + slabColors[s] + ", Orders =");
      for (o in fromSlab[s]) {
        write(" " + o);
      }
      writeln();
    }
  }
}

```

データはわざとロスが出るデータにしています。

```Java:steelmilljp.dat
Params =
<
  10, // 最大の利用可能なスラブ数
   2   // スラブあたりの色数上限
>;

Orders =
{
    <1, 22, 5>
  , <2,  9, 3>
  , <3,  9, 4>
  , <4,  8, 5>
  , <5,  8, 7>
  , <6,  6, 3>
  , <7,  5, 6>
  , <8,  3, 0>
  , <9,  3, 2>
  , <10, 3, 3>
  , <11, 3, 1>
  , <12, 3, 5>
};

// 利用可能なスラブの総数
// サイズ 0 のスラブはスラブがないことを表します
slabWeights = { 0, 10, 15, 25 };
```

# 参考
#### 在庫問題 - IBM Documentation 

https://www.ibm.com/docs/ja/icos/20.1.0?topic=problem-inventory

#### csplib038: Steel Mill Slab Design

https://www.csplib.org/Problems/prob038/

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

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