CPLEXサンプル 工場計画1を読み解く

CPLEXのサンプルの[工場計画I](https://www.ibm.com/docs/ja/icos/20.1.0?topic=library-manufacturing#usroplsamples.uss_opl_modlib.1019740__usroplsamples.uss_opl_modlib.1019805)Factory Planning（CPLEX OPL CPLEXサンプル MIP）を読み解きます。

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

* 実行環境
    * 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) の Factory Planning Iを題材としています。多期間製品混合です。
最大の粗利を上げられる組合せで、どの製品を毎月何個製造し、何個販売し、何個在庫するかを決定します。
例えば、以下のような制約があります。
各生産プロセスには各月の生産キャパシティーがあります。また各製品の需要には各月の最大需要数があります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/706e4293-3ef2-88ec-08cc-593a452e7820.png)


工場計画 1 の例の特性は次のとおりです。

業種: 製造
機能: OPL プロジェクト
技法: 混合整数計画（線形計画法)
複雑さ: 基本
キーワード: 多期間、混合
プロジェクトの場所: examples/opl/models/FactoryPlanning/factoryPlanning1.mod
モデル名: factoryPlanning1.mod

* サンプル

●コード

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

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

●サンプルデータ

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



## 利用データ
1.  製品種ごとの利益

| |利益|
|:----|----:|
|Prod1|10|
|Prod2|6|
|Prod3|8|
|Prod4|4|
|Prod5|11|
|Prod6|9|
|Prod7|3|


２.製品種ごとの必要プロセス時間
| |Prod1|Prod2|Prod3|Prod4|Prod5|Prod6|Prod7|
|----:|----:|----:|----:|----:|----:|----:|----:|
|研磨(Grind)|0.5|0.7|0|0|0.3|0.2|0.5|
|Vドリル(VDrill)|0.1|0.2|0|0.3|0|0.6|0|
|HDドリル(HDrill)|0.2|0|0.8|0|0|0|0.6|
|くり抜き(Bore)|0.05|0.03|0|0.07|0.1|0|0.08|
|かんながけ(Plane)|0|0|0.01|0|0.05|0|0.05|

３.  各月の製品種毎の最大需要数

| |Prod1|Prod2|Prod3|Prod4|Prod5|Prod6|Prod7|
|----:|----:|----:|----:|----:|----:|----:|----:|
|1|500|1000|300|300|800|200|100|
|2|600|500|200|0|400|300|150|
|3|300|600|0|0|500|400|100|
|4|200|300|400|500|200|0|100|
|5|0|100|500|100|1000|300|0|
|6|500|500|100|300|1100|500|60|



４.各プロセスの毎月の稼働可能台数

|月|1|2|3|4|5|6|
|:----|:----|:----|:----|:----|:----|:----|
|研磨(Grind)|3|4|4|4|3|4|
|Vドリル(VDrill)|2|2|2|1|1|2|
|HDドリル(HDrill)|3|1|3|3|3|2|
|くり抜き(Bore)|1|1|0|1|1|1|
|かんながけ(Plane)|1|1|1|1|1|0|

ひと月の稼働時間/台:384

つまり、一月の研磨は3台*384H＝1,152H行うことが可能です。

４．在庫関連の情報
１製品の在庫コスト/月：0.5
初期在庫数： 0
最終在庫数: 50
最大在庫可能数:100

## 決定変数
0. 各製品の毎月の製造数
0. 各製品の毎月の販売数
0. 各製品の毎月の在庫数


## 目的関数
粗利益＝全製品利益-全製品の在庫コスト


## 制約
0. 在庫数が最大在庫可能数以下
0. 販売数が最大需要数以下
0. 毎月の各プロセスの製造可能数以下
0. 前月在庫数＋製造数と販売数＋当月在庫が等しい
0. 初期在庫数と最終在庫数は固定値


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

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

## 目的関数
粗利益の合計:93,711.5

## 決定変数
１. 各製品の毎月の製造数

|月|1|2|3|4|5|6|
|----:|----:|----:|----:|----:|----:|----:|
|Prod1|500|700|0|200|0|550|
|Prod2|888|600|0|300|100|550|
|Prod3|383|117|0|400|600|0|
|Prod4|300|0|0|500|100|350|
|Prod5|800|500|0|200|1100|0|
|Prod6|200|300|400|0|300|550|
|Prod7|0|250|0|100|100|0|


２. 各製品の毎月の販売数

|月|1|2|3|4|5|6|
|----:|----:|----:|----:|----:|----:|----:|
|Prod1|500|600|100|200|0|500|
|Prod2|888|500|100|300|100|500|
|Prod3|300|200|0|400|500|50|
|Prod4|300|0|0|500|100|300|
|Prod5|800|400|100|200|1000|50|
|Prod6|200|300|400|0|300|500|
|Prod7|0|150|100|100|0|50|

３．各製品の毎月の在庫数

|月|0|1|2|3|4|5|6|
|----:|----:|----:|----:|----:|----:|----:|----:|
|Prod1|0|0|100|0|0|0|50|
|Prod2|0|0|100|0|0|0|50|
|Prod3|0|83|0|0|0|100|50|
|Prod4|0|0|0|0|0|0|50|
|Prod5|0|0|100|0|0|100|50|
|Prod6|0|0|0|0|0|0|50|
|Prod7|0|0|100|0|0|100|50|


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

## 利用データ

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

```java:製品名とプロセス名の定義(mod)
// 製品名
{string} Prod = ...;
// プロセス名
{string} Process = ...;
```


```java:製品名とプロセス名の定義(dat)
// 製品名
Prod = {Prod1, Prod2, Prod3, Prod4, Prod5, Prod6, Prod7};
//　プロセス名
Process = {Grind, VDrill, HDrill, Bore, Plane};
```


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

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

```java:計画を立てる月数(dat)
//計画を立てる月数
NbMonths = 6;
```


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

```java:製品関係(mod)
// 製品種ごとの利益
float ProfitProd[Prod] = ...;
// 製品種ごとの必要プロセス時間[プロセス,月]
float ProcessProd[Process][Prod] = ...;
// 各月の製品種毎の最大需要数[月,製品]
int MarketProd[Months][Prod] = ...;
```
```java:製品関係(dat)
// 製品種ごとの利益
// profitProd[j] is profit per unit for product j
ProfitProd = [10 6 8 4 11 9 3];

// 製品種ごとの必要プロセス時間[プロセス,月]
// processProd[i,j] gives hours of process i required by product j
ProcessProd = [[0.5  0.7  0.0  0.0  0.3  0.2  0.5 ]
               [0.1  0.2  0.0  0.3  0.0  0.6  0.0 ]
               [0.2  0.0  0.8  0.0  0.0  0.0  0.6 ]
               [0.05 0.03 0.0  0.07 0.1  0.0  0.08]
               [0.0  0.0  0.01 0.0  0.05 0.0  0.05]];

// 各月の製品種毎の最大需要数[月,製品]
// marketProd[i,j] gives marketing limitation on product j for month i
MarketProd = [[ 500 1000  300  300  800  200  100]
              [ 600  500  200    0  400  300  150]
              [ 300  600    0    0  500  400  100]
              [ 200  300  400  500  200    0  100]
              [   0  100  500  100 1000  300    0]
              [ 500  500  100  300 1100  500   60]]; 
```
問題ブラウザで見ると以下のようになっています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/37807c00-e4ad-c36c-1369-20de135dfcb5.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/ac896d1c-ae05-c823-523a-aaf3a9757fef.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/5cd3dc61-70f8-4cae-6836-f5a8fccf382b.png)


プロセスの稼働可能台数がNumProcessの定義です。
processProdの定義は時間単位なので、NumProcessにHoursMonthを掛ける必要があります。


```java:プロセスの稼働可能時間(mod)
// ひと月の稼働時間/台
float HoursMonth = ...;
// 各プロセスの毎月の稼働可能台数[プロセス,月]
float NumProcess[Process][Months] = ...;
```
```java:プロセスの稼働可能時間(dat)
// ひと月の稼働時間/台
HoursMonth = 384; // 2 eight hour shifts per day, 24 working days per month;
// 各プロセスの毎月の稼働可能台数[プロセス,月]
NumProcess = [[3 4 4 4 3 4]
              [2 2 2 1 1 2]
              [3 1 3 3 3 2]
              [1 1 0 1 1 1]
              [1 1 1 1 1 0]];
```

問題ブラウザでNumProcessを見ると以下のようになっています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/2d2b6836-90a9-8c0e-22e6-e0d3750366d9.png)


在庫関係の定義です。

```java:在庫関係(mod)
// １製品の在庫コスト/月
float CostHold = ...;
// 初期在庫数
float StartHold = ...;
// 最終在庫数
float EndHold = ...;
// 最大在庫可能数
int MaxHold = ...;
```
```java:在庫関係(dat)
// １製品の在庫コスト/月
CostHold  = 0.5;
// 初期在庫数
StartHold = 0;
// 最終在庫数
EndHold   = 50;
// 最大在庫可能数
MaxHold   = 100;
```

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

Holdのx軸はMonthsではなくR0になっています。0月が初期在庫になります。
Holdは「制約1:在庫数が最大在庫可能数以下」を一緒に定義しています。
Sellは「制約2:販売数が最大需要数以下」を一緒に定義しています。

```java:決定変数
//// 決定変数
// 製造数[製品,月]
dvar int+ Make[Prod][Months];
// 在庫数[製品,初期在庫+月]。制約1:在庫数が最大在庫可能数以下
dvar int+ Hold[Prod][R0] in 0..MaxHold;
// 販売数[製品,月]。制約2:販売数が最大需要数以下
dvar int+ Sell[j in Prod][m in Months] in 0..MarketProd[m,j];
```

これらの決定変数は、求められた解では以下のようになりました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/64f36fe4-20bf-b2bf-3c64-25d83baac1df.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/97356772-4344-e55e-2ae4-24b038add2d6.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/8600ec04-ac64-5ef9-22dd-54c286ccb0c4.png)



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

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


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


```java:決定式のコスト（Cost）、利益（Profit）、粗利益(GrossProfit)
// 利益＝製品ごとの利益*販売数
dexpr float Profit = 
  sum (j in Prod, m in Months) ProfitProd[j] * Sell[j][m];
// コスト＝在庫コスト*在庫量
dexpr float Cost = 
  sum (j in Prod, m in Months) CostHold * Hold[j][m];
// 粗利=利益-コスト
dexpr float GrossProfit = Profit - Cost;
```

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/931709cb-4892-f6b1-cb08-75221a049664.png)


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

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

---略----
}

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

以下の2つの制約は決定変数の定義の中で定義済みです。
制約1:在庫数が最大在庫可能数以下
制約2:販売数が最大需要数以下

### 制約3: 毎月の各プロセスの製造可能数以下
forall(m in Months, i in Process)で月数＊プロセス数分の制約を作っています。
左辺「製品種ごとの必要プロセス時間（ProcessProd）」 * 「製造数（Make）」で「各プロセスごとの月々の稼働時間」を計算しています。
右辺「各プロセスの毎月の稼働可能台数(NumProcess) 」* 「ひと月の稼働時間/台（HoursMonth）」
で「毎月の各プロセスの製造可能時間」を計算しています。

```java:制約3: 毎月の各プロセスの製造可能数以下
forall(m in Months, i in Process)
  ct1Capacity: sum(j in Prod) ProcessProd[i][j] * Make[j][m]
         <= NumProcess[i][m] * HoursMonth;
```
問題ブラウザで見ると以下のようになります。スラックで余った稼働時間が確認できます。3月のBoreはスラック0ですので、ここがボトルネックになっていることがわかります。3月のBoreの稼働台数を増やすことで生産量を増やせる可能性があるといえます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/06bbfb59-7c9f-e6fb-8d14-4c671cd72594.png)


図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/3b717e95-ba04-b89f-6a3d-912b562d796c.png)


### 制約4:前月在庫数＋製造数と販売数＋当月在庫が等しい
forall(j in Prod, m in Months)で製品＊月数分の制約を作っています。
前月在庫数＋製造数と販売数＋当月在庫が一致するという流量保存の制約です。

```java:制約4:前月在庫数＋製造数と販売数＋当月在庫が等しい
forall(j in Prod, m in Months)
  ct2InvBal: Hold[j][m-1] + Make[j][m] == Sell[j][m] + Hold[j][m];  
}  
```
問題ブラウザで見ると以下のようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/af745275-23a7-cff5-c5b7-68c59d18b7c5.png)


図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/cafa2f92-696e-c9ed-7d2e-fb5a13e1cddf.png)




### 制約5:初期在庫数と最終在庫数は固定値
forall(j in Prod)で製品数分の制約を作っています。
初期在庫(Hold[j][0])を0(StartHold)に固定している。
最終月在庫(Hold[j][NbMonths])を50(EndHold)に固定している。


```java:制約5:初期在庫数と最終在庫数は固定値
forall(j in Prod) {
  ct3StartInv: Hold[j][0] == StartHold;
  ct4EndInv: Hold[j][NbMonths] == EndHold;
}
```
問題ブラウザで見ると以下のようになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/17879257-4af8-805f-9602-1e061ce4c0ff.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/c0fc2ce6-3f90-3a75-dea6-6fc3f9da1574.png)



図にすると以下のようなイメージになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/a6d5f93f-f486-5bdc-4d4b-97a2194f83fa.png)

## 結果の表示
解の決定変数を表示します。
plan[m][j] は、各製品 j を毎月 m にどれだけ製造、販売、保有するかを出力します。

```java:制約（mod）
execute DISPLAY {
   //plan[m][j] は、各製品 j を毎月 m にどれだけ製造、販売、保有するかを記述します。
   for(var m in Months)
      for(var j in Prod)
         writeln("plan[",m,"][",j,"] = <Make:",Make[j][m],", Sell:",Sell[j][m],", Hold:",Hold[j][m],">");
}
```
スクリプト・ログのビューに出力されます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/223763/ba0364b8-8266-8750-c527-0ba75d26a98a.png)



# 完成OPL

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

```Java:factoryPlanning1hkwd.mod
// 工場計画I
// https://www.ibm.com/docs/ja/icos/20.1.0?topic=library-manufacturing#usroplsamples.uss_opl_modlib.1019740__usroplsamples.uss_opl_modlib.1019805


//// データ
// 製品名
{string} Prod = ...;
// プロセス名
{string} Process = ...;

// 計画を立てる月数
int NbMonths = ...;
range Months = 1..NbMonths;
//初期在庫＋月。在庫の制約を最初の月とそれ以降で分けないため
range R0 = 0..NbMonths;

// 製品種ごとの利益
float ProfitProd[Prod] = ...;
// 製品種ごとの必要プロセス時間[プロセス,月]
float ProcessProd[Process][Prod] = ...;
// 各月の製品種毎の最大需要数[月,製品]
int MarketProd[Months][Prod] = ...;

// ひと月の稼働時間/台
float HoursMonth = ...;
// 各プロセスの毎月の稼働可能台数[プロセス,月]
float NumProcess[Process][Months] = ...;

// １製品の在庫コスト/月
float CostHold = ...;
// 初期在庫数
float StartHold = ...;
// 最終在庫数
float EndHold = ...;
// 最大在庫可能数
int MaxHold = ...;



//// 決定変数
// 製造数[製品,月]
//dvar float+ Make[Prod][Months];
dvar int+ Make[Prod][Months];
// 在庫数[製品,初期在庫+月]。制約1:在庫数が最大在庫可能数以下
//dvar float+ Hold[Prod][R0] in 0..MaxHold;
dvar int+ Hold[Prod][R0] in 0..MaxHold;
// 販売数[製品,月]。制約2:販売数が最大需要数以下
//dvar float+ Sell[j in Prod][m in Months] in 0..MarketProd[m,j];
dvar int+ Sell[j in Prod][m in Months] in 0..MarketProd[m,j];

//// 決定関数
// 利益＝製品ごとの利益*販売数
dexpr float Profit = 
  sum (j in Prod, m in Months) ProfitProd[j] * Sell[j][m];
// コスト＝在庫コスト*在庫量
dexpr float Cost = 
  sum (j in Prod, m in Months) CostHold * Hold[j][m];
// 粗利=利益-コスト
dexpr float GrossProfit = Profit - Cost;

  
//// 目的関数
// 粗利の最大化
maximize GrossProfit;
        
//// 制約
subject to {
  // Limits on process capacity
  // 制約3:毎月の各プロセスの製造可能数以下
  forall(m in Months, i in Process)
    ct1Capacity: sum(j in Prod) ProcessProd[i][j] * Make[j][m]
           <= NumProcess[i][m] * HoursMonth;

  // Inventory balance
  // 制約4:前月在庫数＋製造数と販売数＋当月在庫が等しい
  forall(j in Prod, m in Months)
    ct2InvBal: Hold[j][m-1] + Make[j][m] == Sell[j][m] + Hold[j][m];

  // Starting and ending inventories are fixed
  // 制約5:初期在庫数と最終在庫数は固定値
  forall(j in Prod) {
    ct3StartInv: Hold[j][0] == StartHold;    
    ct4EndInv: Hold[j][NbMonths] == EndHold;
  }
}


// 結果の表示
execute DISPLAY {
   //plan[m][j] describes how much to make, sell, and hold of each product j in each month m
   // plan[m][j] は、各製品 j を毎月 m にどれだけ製造、販売、保有するかを記述します。
   for(var m in Months)
      for(var j in Prod)
         writeln("plan[",m,"][",j,"] = <Make:",Make[j][m],", Sell:",Sell[j][m],", Hold:",Hold[j][m],">");
}
```

データは以下です。

```Java:factoryPlanning1hkwd.dat
//計画を立てる月数
NbMonths = 6;

// 製品名
Prod = {Prod1, Prod2, Prod3, Prod4, Prod5, Prod6, Prod7};
// プロセス名
Process = {Grind, VDrill, HDrill, Bore, Plane};


// 製品種ごとの利益
// profitProd[j] is profit per unit for product j
ProfitProd = [10 6 8 4 11 9 3];

// 製品種ごとの必要プロセス時間[プロセス,月]
// processProd[i,j] gives hours of process i required by product j
ProcessProd = [[0.5  0.7  0.0  0.0  0.3  0.2  0.5 ]
               [0.1  0.2  0.0  0.3  0.0  0.6  0.0 ]
               [0.2  0.0  0.8  0.0  0.0  0.0  0.6 ]
               [0.05 0.03 0.0  0.07 0.1  0.0  0.08]
               [0.0  0.0  0.01 0.0  0.05 0.0  0.05]];

// 各月の製品種毎の最大需要数[月,製品]
// marketProd[i,j] gives marketing limitation on product j for month i
MarketProd = [[ 500 1000  300  300  800  200  100]
              [ 600  500  200    0  400  300  150]
              [ 300  600    0    0  500  400  100]
              [ 200  300  400  500  200    0  100]
              [   0  100  500  100 1000  300    0]
              [ 500  500  100  300 1100  500   60]]; 


// Process capacity
// ひと月の稼働時間/台
HoursMonth = 384; // 2 eight hour shifts per day, 24 working days per month;

// 各プロセスの毎月の稼働可能台数[プロセス,月]
NumProcess = [[3 4 4 4 3 4]
              [2 2 2 1 1 2]
              [3 1 3 3 3 2]
              [1 1 0 1 1 1]
              [1 1 1 1 1 0]];

         
// １製品の在庫コスト/月
CostHold  = 0.5;
// 初期在庫数
StartHold = 0;
// 最終在庫数
EndHold   = 50;
// 最大在庫可能数
MaxHold   = 100;
```

# 参考
#### 製造：工場計画I
https://www.ibm.com/docs/ja/icos/20.1.0?topic=library-manufacturing#usroplsamples.uss_opl_modlib.1019740__usroplsamples.uss_opl_modlib.1019805


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

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