GoogleのOR-Toolsで数理最適化モデルの練習問題を解く(3)生産最適化問題

#はじめに

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」の練習問題を解く第三回目の記事です。

第一回目の記事は以下です。まずはこちらをご覧下さい。
>GoogleのOR-Toolsで数理最適化モデルの練習問題を解く(1)一番易しいマス埋め問題
>https://qiita.com/ttlabo/private/7e8c93f387814f99931f

#生産最適化問題

参考テキストに出てくる練習問題の3番目です。
以下の問題をさっそく行ってみましょう。

>:speech_balloon: 問題

```ruby:7.3.問題
原料に限りがある状態で製品を生産して、利益を最大化したい。
製品ごとの原料価格と利益、原料ごとの在庫は以下の通り。
```
![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/1874e1fd-04ef-1705-4be6-1aa5a9e5f00e.jpeg)

>:question: **考え方**

この問題は、典型問題と言われる数理最適化の問題のうち「ロジスティクス・ネットワーク設計問題」の一つです。
制約条件と目的関数を考えます。

**制約条件**
　問題文には「原料に限りがある」と書かれています。これが制約条件を考える手がかりです。

**目的関数**
　利益を最大化することが目的関数のヒントになります。

>:a: **解答**

制約条件を考えます。
製品1を生産する量をx、製品2を生産する量をyとします。
x,yは正の整数値を取ります。

**制約条件**
原料1は在庫が40あります。製品1と製品2に使用する原料1は40以下となります。
`1 * x + 2 * y <= 40`

次に、原料2は在庫が80あります。製品1と製品2に使用する原料2は80以下となります。
`4 * x + 4 * y <= 80`

同様に原料3について、以下が成り立ちます。
`3 * x + 1 * y <= 50`

**目的関数**
利益を考えます。
製品1から得られる利益は5.0、製品2から得られる利益は4.0です。
これらの利益の合計が最大になることが目的関数です。
`利益 = x * 5.0 + y * 4.0`

プログラムを検討します。
プログラムの記載内容については、基本的にGoogleのOR-Toolsのガイドに沿っています。
(https://developers.google.com/optimization)


プログラムの冒頭におまじないを書きます。

```ruby:7.3.renshu.py
from __future__ import print_function
from ortools.linear_solver import pywraplp
```

混合整数計画ソルバーで解きますので、以下に宣言します。

```ruby:7.3.renshu.py
# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
                             pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
```

x,yはマイナス値を取らない条件を付けます。
x,yの最大値を100としておきます。

```ruby:7.3.renshu.py
# x,yについての非負条件(マイナス値でないこと)
x = solver.IntVar(0.0, 100, 'x')
y = solver.IntVar(0.0, 100, 'y')
```

3つの制約条件を定義します。

```ruby:7.3.renshu.py
solver.Add(1 * x + 2 * y <= 40)
solver.Add(4 * x + 4 * y <= 80)
solver.Add(3 * x + 1 * y <= 50)
```

本問題は最大化問題なので、solver.Maximizeを使用します。

```ruby:7.3.renshu.py
#利益を最大化する
solver.Maximize(x * 5.0 + y * 4.0)
```

求解を実行します。

```7.3.renshu.py
# 求解実行
status = solver.Solve()
```

結果を確認します。

```7.3.renshu.py
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    if x.solution_value() > 0.5:
        print('x=',x.solution_value())
        print('y=',y.solution_value())
    print("Time = ", solver.WallTime(), " milliseconds")
else:
    print('The problem does not have an optimal solution.')
```

最適化計算結果は以下となりました。

>Solution:
>ok
>Objective value = 95.0
>x= 15.0
>y= 5.0
>Time =  369  milliseconds

製品1を15個、製品2を5個生産すればよいことがわかります。
最大化される利益は95.0です。

プログラム全体です。

```7.3.renshu.py
from __future__ import print_function
from ortools.linear_solver import pywraplp

# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

x = solver.IntVar(0.0, 100, 'x')
y = solver.IntVar(0.0, 100, 'y')

# 制約条件
solver.Add(1 * x + 2 * y <= 40)
solver.Add(4 * x + 4 * y <= 80)
solver.Add(3 * x + 1 * y <= 50)


# 利益を最大化する
solver.Maximize(x * 5.0 + y * 4.0)

# 求解実行
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    if x.solution_value() > 0.5:
        print('x=',x.solution_value())
        print('y=',y.solution_value())
    print("Time = ", solver.WallTime(), " milliseconds")
else:
    print('The problem does not have an optimal solution.')
```

#出展

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」に記載の練習問題を参考にさせて頂いております。

■参考テキスト
「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」
　斉藤努[著]
　近代科学社[出版]

![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/df8e6661-9dd6-559b-5ae8-0bf79f73d1a5.jpeg)

# ご意見など

ご意見、間違い訂正などございましたらお寄せ下さい。