CPLEXサンプル 配合 1を読み解く

CPLEXのサンプルの[配合 I](https://www.ibm.com/docs/ja/icos/20.1.0?topic=library-manufacturing#d114357e164) 食品製造 1（CPLEX OPL CPLEXサンプル MIP）を読み解きます。

<サンプル導入ディレクトリ>opl\examples/opl/models/BlendingMultiPeriod
にあります。

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

# 問題

このモデルでは、[「Model Building in Mathematical Programming」(H.P. Williams 著)](https://www.amazon.co.jp/Model-Building-Mathematical-Programming-Williams/dp/1118443330) の Food manufacture 1を題材としています。多期間に渡る配合の問題です。

6種類の粗油を毎月何トン仕入れ、何トン保管するか、その組合せでマーガリンを毎月何トン製造するかを決定します。最大の粗利が出るようにします。
6種類の粗油の組合せによってマーガリンの硬度が変わりますが、3－6の硬度で作る必要があります。
また、毎月各粗油の仕入れ値は変動します。また保管には固定費がかかります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/ab08925c-1aca-31b1-348a-aedc02955aec.png)


配合 1 の例の特性は次のとおりです。

業種: 製造
機能: OPL モデル
技法: 線形計画法
複雑さ: 基本
キーワード: 配合、製造、食品、多期間
プロジェクトの場所: examples/opl/models/BlendingMultiPeriod
モデル名: blending1.mod

* サンプル

●コード

https://github.com/hkwd/230823cplexsample/blob/master/240115FactoryPlanning/factoryPlanning1hkwd.mod

この記事でのサンプルは、オリジナルのサンプルとほぼ同じ内容ですが、日本語コメントを入れたり順序を入れ替えたり、添え字を入れ替えたりして読みやすくしています

●サンプルデータ

https://github.com/hkwd/230823cplexsample/blob/master/240115FactoryPlanning/factoryPlanning1hkwd.dat



## 利用データ
1.  各粗油の仕入値/月


| |1|2|3|4|5|6|
|:----|----:|----:|----:|----:|----:|----:|
|Veg1|110|130|110|120|100|90|
|Veg2|120|130|140|110|120|100|
|Oil1|130|110|130|120|150|140|
|Oil2|110|90|100|120|110|80|
|Oil3|115|115|95|125|105|135|



２.各粗油の硬度
|粗油|硬度|
|----:|----:|
|Veg1|8.8|
|Veg2|6.1|
|Oil1|2|
|Oil2|4.2|
|Oil3|5|


３. その他の情報
|データ|値|
|:----|----:|
|製品利益|150|
|植物油の最大生産量/月|200|
|非植物油の最大生産量/月|250|
|保管コスト/月|5|
|粗油の初期保管量と最終月保管量|500|
|各原料油の最大保管量|1000|


## 決定変数
0. 各粗油の各月の購入量
0. 各粗油の各月の保管量
0. 各粗油の各月の利用量
0. 各月の製品生産量


## 目的関数
全製品利益＝製品生産量＊製品利益
全製品コスト＝各粗油の仕入値＊購入量＋保管量＊保管コスト）

総利益＝全製品利益－全製品コスト


## 制約
0. 在庫量が最大保管量以下
0. 植物油と非植物の月あたりの最大使用量以下
0. 1トンあたりの硬度は3以上6以下
0. 利用量と生産量が等しい
0. 前月在庫量＋購入量と利用量＋当月在庫量が等しい
0. 初期在庫量と最終在庫量は固定値


## 問題の種類
混合整数計画

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

## 目的関数
総利益の合計:107,840

## 決定変数
１. 各粗油の各月の購入量

| 月|1|2|3|4|5|6|
|:----|----:|----:|----:|----:|----:|----:|
|Veg1|0|0|0|0|0|659|
|Veg2|0|0|0|0|0|541|
|Oil1|0|0|0|0|0|0|
|Oil2|0|750|0|0|0|750|
|Oil3|0|0|0|0|0|0|



２. 各粗油の各月の保管量

| 月|初期在庫|1|2|3|4|5|6|
|:----|----:|----:|----:|----:|----:|----:|----:|
|Veg1|500|477|477|318|159|0|500|
|Veg2|500|323|123|82|41|0|500|
|Oil1|500|500|500|500|500|500|500|
|Oil2|500|250|750|500|250|0|500|
|Oil3|500|500|500|500|500|500|500|


３．各粗油の各月の利用量

|月|1|2|3|4|5|6|
|:----|----:|----:|----:|----:|----:|----:|
|Veg1|23|0|159|159|159|159|
|Veg2|177|200|41|41|41|41|
|Oil1|0|0|0|0|0|0|
|Oil2|250|250|250|250|250|250|
|Oil3|0|0|0|0|0|0|

４．各月の製品生産量

|月|生産量|
|----:|----:|
|1|450|
|2|450|
|3|450|
|4|450|
|5|450|
|6|450|

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

## 利用データ

まず製品名とプロセス名のデータの定義を行っています。
データの中身の定義はblending1hkwd.datファイルの中で行っています。

```java:粗油名の定義(mod)
// 各粗油名
{string} Raw = ...;
```


```java:粗油名の定義(dat)
// 各粗油名
Raw = {Veg1, Veg2, Oil1, Oil2, Oil3};
```


生産月数として1月－6月を定義しています。
Rは1月－6月に加えて、0月も定義しています。これは在庫テーブルの月の情報です。在庫に関しては制約を1月と2月以降で別のものを定義しなくてよいようにというテクニカルな理由と、1月以前に在庫があった場合にも対応できるようにするという意図があります。

```java:生産月数(mod)
//　生産月数
int NbMonths = ...;
range Months = 1..NbMonths;
// 初期在庫＋月。在庫の制約を最初の月とそれ以降で分けないため
range R = 0..NbMonths;
```

```java:生産月数(dat)
// 生産月数
NbMonths = 6;
```


粗油と製品関係のデータを定義しています。

```java:粗油と製品関係(mod)
// 各粗油の仕入値/月
float CostRaw[Raw][Months] = ...;
//　各粗油の硬度
float HardRaw[Raw] = ...;
// 製品利益
float ProfitProd = ...;
// 植物油かそれ以外かの識別フラグ
float IsVeg[Raw] = ...;
float IsOil[Raw] = ...;
// 植物油の最大生産量/月
float MaxVeg = ...;
// 非植物油の最大生産量/月
float MaxOil = ...;
```
```java:粗油と製品関係(dat)
// 各粗油の仕入値/月
CostRaw = [[110 130 110 120 100  90]
           [120 130 140 110 120 100]
           [130 110 130 120 150 140]
           [110  90 100 120 110  80]
           [115 115  95 125 105 135]];
//　各粗油の硬度
HardRaw = [8.8 6.1 2 4.2 5];
// 植物油かそれ以外かの識別フラグ
IsVeg = [1 1 0 0 0];
IsOil = [0 0 1 1 1];
// 製品利益
ProfitProd = 150;
// 植物油の最大生産量/月
MaxVeg = 200;
// 非植物油の最大生産量/月
MaxOil = 250;
```
問題ブラウザで見ると以下のようになっています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/ca5b3e9c-8998-583e-9697-e81c78bbdb1a.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/f51bc4da-e12d-fa93-1a50-a65949960b59.png)
IsVegとIsOilは、植物油と非植物油を見分けるためのフラグになっています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/d6ced1c5-5e94-611c-8ad9-0c44907d4c2c.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/fdbf52a1-e45c-30cb-9278-19e4647bf606.png)

保管系データの定義です。

```java:保管系データ(mod)
// 保管コスト/月
float CostStore = ...;
// 粗油の初期保管量と最終月保管量
float StartEndStore = ...;
// 各原料油の最大保管量
int MaxStore = ...;
```
```java:保管系データ(dat)
// 保管コスト/月
CostStore = 5;
// 粗油の初期保管量と最終月保管量
StartEndStore = 500;
// 各原料油の最大保管量
MaxStore = 1000;
```

## 決定変数
次に決定変数の定義を行っています。
dvarが決定変数を意味します。

Storeのy軸はMonthsではなくRになっています。0月が初期在庫になります。
Storeは「制約1:在庫数が最大保管量以下」を一緒に定義しています。

```java:決定変数
//// 決定変数
// 購入量[各粗油,月]
dvar int+ Buy[Raw][Months];
// 保管量[各粗油,月]　制約１：最大保管量以下
dvar int+ Store[Raw][R] in 0..MaxStore;
// 利用量[各粗油,月]
dvar int+ Use[Raw][Months];
// 製品生産量[月]
dvar int+ p[Months];
```

これらの決定変数は、求められた解では以下のようになりました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/692688b4-1753-3151-3fe6-7ab1a2bfb9df.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/247e646f-3dc7-c8a3-6fb7-8e9bf2ab306c.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/94001c4e-ff29-1e61-4b49-a7bc5cd6888b.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/08fa3e0d-d622-082a-0962-b4683e232bc2.png)




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

関数の定義に先立って、決定式のコスト（Cost）、利益（Profit）、粗利益(GrossProfit)の決定関数を作ります。


[sumキーワード](https://www.ibm.com/docs/ja/icos/20.1.0?topic=keywords-sum)は関連式の集合の合計値を計算するキーワードです。
各粗油の仕入値(CostRaw)は粗油毎月毎に異なります。
製品利益(ProfitProd)在庫コスト(CostStore)は粗油と毎月共通です。


```java:決定式のコスト（Cost）、利益（Profit）、粗利益
// 利益＝製品利益＊製品生産量
dexpr float Profit = 
  sum (m in Months) ProfitProd * p[m];
// コスト＝各粗油の仕入値＊購入量＋保管コスト＊保管量
dexpr float Cost =
  sum(j in Raw, m in Months)
         (CostRaw[j][m] * Buy[j][m] + CostStore * Store[j][m]);
// 粗利=利益-コスト
dexpr float GrossProfit = Profit - Cost;
```

図にすると以下のようになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/2679a15f-9a45-e460-a021-35a5ad693e07.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/978e6f51-5a40-a292-63fe-24af8ab2ba09.png)




決定関数として定義した粗利益を最大化することが目的関数になります。

```java:目的関数 粗利益の最大化
maximize GrossProfit;
```
## 制約
次に制約の定義を行っています。
subject toの中カッコで囲みます。
```java:制約（mod）
subject to{

---略----
}

```
ctXXX:で制約にラベルを付けています。

以下の制約は決定変数の定義の中で定義済みです。
制約1:在庫量が最大保管量以下


### 制約2: 植物油と非植物油は月あたりの最大使用量以下
forall(m in Months)で月数分の制約を作っています。
「植物油の識別フラグ（IsVeg）」 * 「使用量（Use）」で「植物油の使用量」を計算し、それがMaxVeg以下である条件になっています。
非植物油についても同様の条件をつくっています。



```java:制約2:植物油と非植物油は月あたりの最大使用量以下
ctMaxUseVeg: sum(j in Raw) IsVeg[j] * Use[j][m] <= MaxVeg;
ctMaxUseOil: sum(j in Raw) IsOil[j] * Use[j][m] <= MaxOil;
```
問題ブラウザで見ると以下のようになります。スラックは0ですので、全月で植物油と非植物油の使用量限界まで使っていることがわかります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/a57a83fd-c518-5969-fbcc-7037b170b7f5.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/3d88ce0e-dfd4-207e-e943-218198260be2.png)




図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/7369ea06-6328-e112-7cb3-ead2aac9502b.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/dd44b801-33dd-b82c-d383-4c12927987c5.png)



### 制約3:1トンあたりの硬度は3以上6以下
forall(m in Months)で月数分の制約を作っています。

右辺では、「各粗油の硬度（HardRaw）」 * 「使用量（Use）」で製造したマーガリンの総硬度を出しています。

左辺は3 *「生産量（p）」以上、6 * 「生産量（p）」以下の硬度であることを条件にしています。
```java:制約3:1トンあたりの硬度は3以上6以下
ctHard1: sum(j in Raw) HardRaw[j] * Use[j][m]   >= 3 * p[m];
ctHard2: sum(j in Raw) HardRaw[j] * Use[j][m]   <= 6 * p[m];
```
問題ブラウザで見ると以下のようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/184f7165-b3d5-f49e-bb78-6a497782af4e.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/6311abe2-87d8-a9d4-6138-a9477252cb9d.png)



図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/023104f7-41cc-0006-121a-c33db900b26a.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/e148f13f-2dbd-7a56-ee82-d7f8f28e4bd0.png)


これは本当は以下のように、左辺を「p[m]」で割ってあげた方1トン当たりの硬度になるので、人間の可読性は高まります。しかしながら、決定変数を分母に割り算をすると線形の問題にならなくなるので、「p[m]」は右辺で掛け算にしています。
```java:制約3:1トンあたりの硬度は3以上6以下
ctHard1: sum(j in Raw) HardRaw[j] * Use[j][m]/p[m]  >= 3 * ;
ctHard2: sum(j in Raw) HardRaw[j] * Use[j][m]/p[m]   <= 6 * ;
```

#### 参考:決定変数を分母にした割合で制約を作る
https://qiita.com/spssfun2017/items/4d148487dd6b9a48cf88



### 制約4:利用量と生産量が等しい
forall(m in Months)で月数分の制約を作っています。
すべての粗油の「使用量（Use）」==「生産量（p）」という条件です。

```java:制約4:利用量と生産量が等しい
ctMatBal: sum(j in Raw) Use[j][m] == p[m];
```

問題ブラウザで見ると以下のようになります。
5月はスラックが0ではありませんが、e-14であり誤差です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/112471be-84e4-fa76-34b7-bf97c7e490f6.png)


図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/106a8a43-dce3-959a-5468-16eca36c111f.png)

### 制約5:前月在庫量＋購入量と利用量＋当月在庫量が等しい
forall(j in Raw, m in Months)で粗油分＊月数分の制約を作っています。
前月在庫量＋購入量と利用量＋当月在庫量が一致するという流量保存の制約です。

```java:制約5:前月在庫量＋購入量と利用量＋当月在庫量が等しい
forall(j in Raw, m in Months)
  ctInvVal: Store[j][m-1] + Buy[j][m] == Use[j][m] + Store[j][m];
```
問題ブラウザで見ると以下のようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/80cd256f-06a4-7e31-1455-456a03160d6a.png)

図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/18548815-75e7-38df-9414-d30d7fd29d0b.png)

### 制約6:初期在庫量と最終在庫量は固定値
forall(j in Raw) で粗油数分の制約を作っています。
初期在庫(Store[j][0])を500(StartEndStore)に固定している。
最終月在庫(Store[j][NbMonths])を500(StartEndStore)に固定している。


```java:制約6:初期在庫量と最終在庫量は固定値
forall(j in Raw) {
  ctStartInv: Store[j][0] == StartEndStore;
  ctEndInv: Store[j][NbMonths] == StartEndStore;
}
```
問題ブラウザで見ると以下のようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/0baec4cc-a28c-7caa-a506-114afbfd12b9.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/8562d21a-ae8a-f643-a742-30067f35aa85.png)



図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/e037fb58-433a-91be-005c-9ba0309ee6a0.png)


## 結果の表示
解の決定変数を表示します。
plan[m][j] は、各粗油 j を毎月 m にどれだけ購入、使用、保管するかを出力します。

```java:結果の表示（mod）
execute DISPLAY {
   for (var m in Months)
      for (var j in Raw)
         writeln("plan[",m,"][",j,"] = <buy:",Buy[j][m],",use:",Use[j][m],",store:",Store[j][m],">");
```
スクリプト・ログのビューに出力されます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/930a5ead-9b11-b888-6ad0-942e86391a0f.png)



# 完成OPL

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

```Java:blending1hkwd.mod
// 食品製造 1

// 各粗油名
{string} Raw = ...;

//　生産月数
int NbMonths = ...;
range Months = 1..NbMonths;
// 初期在庫＋月。在庫の制約を最初の月とそれ以降で分けないため
range R = 0..NbMonths;

// 各粗油の仕入値/月
float CostRaw[Raw][Months] = ...;
//　各粗油の硬度
float HardRaw[Raw] = ...;
// 製品利益
float ProfitProd = ...;
// 植物油かそれ以外かの識別フラグ
float IsVeg[Raw] = ...;
float IsOil[Raw] = ...;
// 植物油の最大生産量/月
float MaxVeg = ...;
// 非植物油の最大生産量/月
float MaxOil = ...;


// 保管コスト/月
float CostStore = ...;
// 粗油の初期保管量と最終月保管量
float StartEndStore = ...;
// 各原料油の最大保管量
int MaxStore = ...;


//// 決定変数
// 購入量[各粗油,月]
dvar int+ Buy[Raw][Months];
// 保管量[各粗油,月]　制約１：最大保管量以下
dvar int+ Store[Raw][R] in 0..MaxStore;
// 利用量[各粗油,月]
dvar int+ Use[Raw][Months];
// 製品生産量[月]
dvar int+ p[Months];


// 利益＝製品利益＊製品生産量
dexpr float Profit = 
  sum (m in Months) ProfitProd * p[m];
// コスト＝各粗油の仕入値＊購入量＋保管コスト＊保管量
dexpr float Cost =
  sum(j in Raw, m in Months)
         (CostRaw[j][m] * Buy[j][m] + CostStore * Store[j][m]);
// 粗利=利益-コスト
dexpr float GrossProfit = Profit - Cost;

//　目的関数 粗利益の最大化
maximize GrossProfit;


//　制約
subject to {

  forall(m in Months) {
      
    // Maximum usage per month   
    // 制約2:植物油と非植物油の月あたりの最大使用量以下
    ctMaxUseVeg: sum(j in Raw) IsVeg[j] * Use[j][m] <= MaxVeg;
    ctMaxUseOil: sum(j in Raw) IsOil[j] * Use[j][m] <= MaxOil;
     
    // Hardness constraints
    // 制約3:1トンあたりの硬度は3以上6以下
    //ctHard1: sum(j in Raw) HardRaw[j] * Use[j][m] - 6 * p[m] <= 0;
    //ctHard2: sum(j in Raw) HardRaw[j] * Use[j][m] - 3 * p[m] >= 0;
    ctHard1: sum(j in Raw) HardRaw[j] * Use[j][m]   >= 3 * p[m];
    ctHard2: sum(j in Raw) HardRaw[j] * Use[j][m]   <= 6 * p[m];
    //ctHard1: sum(j in Raw) HardRaw[j] * Use[j][m] / (p[m]*1.0) >= 3;
    //ctHard2: sum(j in Raw) HardRaw[j] * Use[j][m] / (p[m]*1.0) <= 6;

    // Material balance constraints
    // 制約4:利用量と生産量が等しい
    //ctMatBal: sum(j in Raw) Use[j][m] - p[m] == 0;
    ctMatBal: sum(j in Raw) Use[j][m] == p[m];
  }

  // Inventory balance
  // 制約5:前月在庫量＋購入量と利用量＋当月在庫量が等しい
  forall(j in Raw, m in Months)
    ctInvVal: Store[j][m-1] + Buy[j][m] == Use[j][m] + Store[j][m];

  // Starting and ending inventories are fixed
  // 制約6:初期在庫量と最終在庫量は固定値
  forall(j in Raw) {
    ctStartInv: Store[j][0] == StartEndStore;
    ctEndInv: Store[j][NbMonths] == StartEndStore;
  }
}


//Display the plan for each month and each raw material
//plan[m][j] = <Buy[j][m], Use[j][m], Store[j][m]>
execute DISPLAY {
   for (var m in Months)
      for (var j in Raw)
         writeln("plan[",m,"][",j,"] = <buy:",Buy[j][m],",use:",Use[j][m],",store:",Store[j][m],">");
}
```

データは以下です。

```Java:blending1hkwd.dat
// 各粗油名
Raw = {Veg1, Veg2, Oil1, Oil2, Oil3};
// 生産月数
NbMonths = 6;

// 各粗油の仕入値/月
CostRaw = [[110 130 110 120 100  90]
           [120 130 140 110 120 100]
           [130 110 130 120 150 140]
           [110  90 100 120 110  80]
           [115 115  95 125 105 135]];
//　各粗油の硬度
HardRaw = [8.8 6.1 2 4.2 5];
// 植物油かそれ以外かの識別フラグ
IsVeg = [1 1 0 0 0];
IsOil = [0 0 1 1 1];
// 製品利益
ProfitProd = 150;
// 植物油の最大生産量/月
MaxVeg = 200;
// 非植物油の最大生産量/月
MaxOil = 250;

// 保管コスト/月
CostStore = 5;
// 粗油の初期保管量と最終月保管量
StartEndStore = 500;
// 各原料油の最大保管量
MaxStore = 1000;
```

# 参考
#### 製造：配合 I
https://www.ibm.com/docs/ja/icos/20.1.0?topic=library-manufacturing#d114357e164

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

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