QUBO++を活用した最適化モデルの実装解説
中長期生産計画の最適化
はじめに
製造業では、将来需要に対応するために「どの製品をどのラインで生産するか」「いつ設備投資を行うか」「どれだけ能力増強すべきか」といった意思決定が重要になります。
しかし、実際では以下のような制約が複雑に絡み合います。
- 製品ごとに生産可能なラインが異なる
- 設備投資によって製造可能製品を増やせる
- 生産能力には上限がある
- ライン能力にも上限がある
- 投資コストを抑えつつ利益を最大化したい
- 複数年にわたる計画を同時に検討したい
本記事では、これらの課題を QUBO（Quadratic Unconstrained Binary Optimization） として定式化し、Pythonライブラリ pyqbpp （QUBO++のpythonライブラリ）を用いて解く実装例を詳しく解説します。
背景
生産計画最適化の難しさ
製造業の中長期生産計画では以下のような問題が発生します。
例えば、
- ある製品の需要が増加する
- しかしこの製品は現在ライン（国内工場）でしか生産できない
- 現在ラインの能力には限界がある
- 別ラインへ生産移管するには設備投資が必要
- 製品の需要とラインの生産能力を応じて、各ラインに分けて生産する
このような状況では、「どこに・いつ・いくら投資すると最大利益になるか」を評価する必要があります。
さらに、
- 投資は複数年にまたがる
- 投資効果は翌年以降も継続する
- 生産ラインは複数存在する
- 製品数も多い
- 複数の投資候補が存在する
ため、組合せ数が急激に増加します。
この種の問題は設備投資計画、生産能力計画（Capacity Planning）、サプライチェーン設計などで頻繁に現れます。
従来はMILP（Mixed Integer Linear Programming）で解くことが一般的ですが、本実装ではQUBOへ変換し、QUBO++で扱える形に定式化しています。
中長期生産計画の最適化
問題定義・データ
問題定義
計画期間中において、
- 各製品の需要を満たす
- 製品生産能力制約を守る
- ライン能力制約を守る
- 必要に応じて設備投資を行う
- 総利益を最大化する
投資計画と生産計画を同時に決定します。
データ概要
具体的な問題では、6製品（P01~P06）、2生産ライン（L01, L02）、3年間（Y26, Y27, Y28）の計画期間を対象として、生産能力増強および設備投資を含めた最適化を行います。ここでは最適化モデルで利用される入力データを整理して説明します。
- 各製品は異なる需要量、生産能力、利益率、設備投資条件を持っています。
- 製品によっては特定ラインでしか生産できず、設備投資によって生産可能ラインを拡張できます。
- 各年度で能力増強や設備投資の意思決定を行います。
現在の製造可能性
| 製品 | L01 | L02 |
|---|---|---|
| P01 | 1 | 0 |
| P02 | 1 | 0 |
| P03 | 1 | 0 |
| P04 | 1 | 1 |
| P05 | 0 | 1 |
| P06 | 0 | 1 |
ある製品を対象ラインで現在製造可能かどうかを表します。
1：現在生産可能、0：生産不可を表します。
製造可能化投資コスト（単位：千万円）
| 製品 | L01 | L02 |
|---|---|---|
| P01 | 0 | 20 |
| P02 | 0 | 50 |
| P03 | 0 | 100 |
| P04 | 0 | 0 |
| P05 | 300 | 0 |
| P06 | 1000 | 0 |
現在生産できないラインで製造可能にするための投資額です。
※ 0：現在製造可能で投資が必要ないことです。
現在の製品別生産能力（単位：万個）
| 製品 | L01 | L02 |
|---|---|---|
| P01 | 15 | 0 |
| P02 | 10 | 0 |
| P03 | 5 | 0 |
| P04 | 8 | 13 |
| P05 | 0 | 13 |
| P06 | 0 | 0 |
各製品が各ラインで生産可能な最大数量を表します。
現在のライン総能力（単位：万個）
| ライン | 能力 |
|---|---|
| L01 | 15 |
| L02 | 20 |
ライン全体で処理可能な生産数量の上限です。製品個別の能力が十分あっても、このライン能力を超えて生産することはできません。
製品能力増強オプション（単位：万個）
| オプション | 増強量 |
|---|---|
| Option 0 | 0 |
| Option 1 | 2 |
| Option 2 | 4 |
各年度で製品ごとに選択可能な能力増強オプションです。投資を実施するとその効果は以降の年度にも継続します。※ 0：製品能力増強しないことです。
製品能力増強コスト（単位: 千万円）
| 製品 | L01 (+2万個) | L01 (+4万個) | L02 (+2万個) | L02 (+4万個) |
|---|---|---|---|---|
| P01 | 10 | 15 | 20 | 30 |
| P02 | 20 | 40 | 20 | 40 |
| P03 | 10 | 10 | 10 | 10 |
| P04 | 20 | 60 | 20 | 60 |
| P05 | 10 | 50 | 10 | 50 |
| P06 | 20 | 40 | 30 | 60 |
製品やラインによって増強コストが異なるため、どの製品へ投資するかが利益に大きく影響します。
ライン能力増強オプション（単位：万個）
| オプション | 増強量 |
|---|---|
| Option 0 | 0 |
| Option 1 | 4 |
| Option 2 | 8 |
各年度で生産ラインごとに選択可能な能力増強オプションです。
※ 0：ライン能力増強しないことです。
ライン能力増強コスト（単位: 千万円）
| ライン | +4万個 | +8万個 |
|---|---|---|
| L01 | 1,000 | 1,000 |
| L02 | 1,500 | 1,500 |
ライン能力増強は製品能力増強と比較して非常に高額な投資であり、慎重な判断が必要です。
年度別需要予測（単位：万個）
| 年度 | P01 | P02 | P03 | P04 | P05 | P06 |
|---|---|---|---|---|---|---|
| Y26 | 18 | 5 | 2 | 7 | 3 | 0 |
| Y27 | 0 | 8 | 3 | 7 | 3 | 0 |
| Y28 | 0 | 10 | 2 | 8 | 4 | 4 |
需要パターンを見ると、製品 P02 は継続的に増加傾向にあり、製品 P06 は Y28 の年度から新たな需要が発生しています。そのため、将来需要を見越した設備投資計画が重要になります。
製品別・ライン別利益（単位：万円）
| 製品 | L01 | L02 |
|---|---|---|
| P01 | 100 | 150 |
| P02 | 80 | 90 |
| P03 | 100 | 110 |
| P04 | 120 | 121 |
| P05 | 200 | 162 |
| P06 | 50 | 60 |
各製品を各ラインで1単位生産した場合の利益を表します。このような利益差が存在するため、単純に需要を満たすだけではなく、設備投資を行って生産ラインを変更した方が全体利益が高くなる場合があります。
数理モデル構成
制約式と目的式
-
c_expr：制約式 -
o_expr：目的式
c_expr = qbpp.Expr()
o_expr = qbpp.Expr()
製品能力増強選択変数 (prod_capup_sel_v)
- 「どの年度に、どの製品を、どのラインで、どれだけ能力増強するか」を表す意思決定変数（バイナリ変数）です。
- 4次元：年度 × 製品 × ライン × 能力増強オプション
-
prod_capup_sel_c制約は各年度・各製品・各ラインについて、能力増強オプションを必ず1つだけ選択する
# 変数: 年度・製品・ラインごとの能力増強オプション選択
prod_capup_sel_v = qbpp.var(
"p_up_sel_v", shape=(n_years, n_prods, n_lines, n_prod_capups)
)
# 制約: 各年度・製品・ラインごとに能力増強オプションを1つだけ選択する
prod_capup_sel_c = qbpp.constrain(qbpp.vector_sum(prod_capup_sel_v), equal=1)
c_expr += qbpp.sum(prod_capup_sel_c)
ライン能力増強選択変数（line_capup_sel_v）
- 「どの年度に、どの生産ラインへ、どれだけライン能力増強を行うか」を表す意思決定変数（バイナリ変数）です。
- 3次元：年度 × ライン × ライン能力増強オプション
-
line_capup_sel_c制約は各年度・各ラインに対して、能力増強オプションを必ず1つだけ選択することを保証しています
# 変数: 年度・ラインごとのライン能力増強オプション選択変数
line_capup_sel_v = qbpp.var("l_up_sel_v", shape=(n_years, n_lines, n_line_capups))
# 制約: 各年度・製品・ラインの組み合わせごとに能力増強オプションを必ず1つ選択する
line_capup_sel_c = qbpp.constrain(qbpp.vector_sum(line_capup_sel_v), equal=1)
c_expr += qbpp.sum(line_capup_sel_c)
製造可能性管理（make_v）
-
make_vは各年度における製造可能状態を表すバイナリ変数です。 - 3次元：年度 × 製品 × ライン
# 変数：特定の年度において、製品が当該ラインで生産可能かどうかを示す。
make_v = qbpp.var("make_v", shape=(n_years, n_prods, n_lines))
-
make_change_c制約は生産可能状態の単調増加制約で、一度生産可能になったら、その後は生産不可に戻らないことを保証します。
# 制約: 製造可能状態は 0→1 のみ変化可能で、1→0 は不可
# make_v[yi-1][pi][li] <= make[yi][pi][li]
make_change_c = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
if yi == 0:
# 初年度は初期状態と比較
make_change_c[yi][pi][li] = qbpp.constrain(
make_v[yi][pi][li] - curr_make_data[pi][li], between=(0, 1)
)
else:
# 2年目以降は前年度状態と比較
make_change_c[yi][pi][li] = qbpp.constrain(
make_v[yi][pi][li] - make_v[yi - 1][pi][li], between=(0, 1)
)
c_expr += qbpp.sum(make_change_c)
- 製品能力増強が選択された場合、その製品が生産可能状態になることを保証する制約を、以下のように定義できます。
-
prod_capup_sel_v[yi][pi][li][pci] = 1の場合、制約式prod_capup_sel_v[yi][pi][li][pci] × (1 - make_v[yi][pi][li]) = 0より、
make_v[yi][pi][li] = 1であることが保証されます。この制約式にはバイナリ変数同士の積が含まれています。 -
多くのQUBOライブラリでは、
constrain(expr, equal=0)をペナルティ項
$P(\mathrm{expr})^2$
に変換して扱います。そのため、制約式のペナルティ項では、内部的には4次の多項式（HUBO: Higher-Order Binary Optimization）として表現される場合があります。QUBO++はHUBOを直接扱えるため、このような制約もシンプルに記述できます。 -
また、この制約を
make_change_cと組み合わせることで、「製品能力増強が行わ れた製品は生産可能状態となる」という条件を満たすことができます。
-
# 制約: 製品能力増強が選択された場合、その製品は生産可能状態になる
#################################
## QUBO++でHUBOを使用する場合 ##
make_c = qbpp.expr(shape=(n_years, n_prods, n_lines, n_prod_capups))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
for pci in range(1, n_prod_capups):
make_c[yi][pi][li][pci] = qbpp.constrain(
prod_capup_sel_v[yi][pi][li][pci] * (1 - make_v[yi][pi][li]),
equal=0,
)
c_expr += qbpp.sum(make_c)
## HUBO使用終了 ##
#################################
- 別の方法として、HUBOを用いずに同等の制約を表現することもできます。
-
make_cum_vは、「対象年度までに製造可能化投資が実施されたか」を累積的に表す変数です。 -
make_c制約は、「累積投資が存在する場合は生産可能状態にする」ための制約です。 -
QUBOでは if 文を直接表現できません。そのため、
make_cum_v[yi][pi][li] > 0 → make_v[yi][pi][li] = 1という論理条件を数式として表現する必要があります。Big-M法は、そのためによく用いられる手法の一つです。
############################### ## QUBO++でHUBOを使用しない ## # 変数: 現年度まで累積した生産可能状態 make_cum_v = qbpp.expr(shape=(n_years, n_prods, n_lines)) for yi in range(n_years): for pi in range(n_prods): for li in range(n_lines): make_cum_v[yi][pi][li] += curr_make_data[pi][li] # 対象年度までの能力増強投資を累積 for i in range(yi + 1): for pci in range(1, n_prod_capups): make_cum_v[yi][pi][li] += prod_capup_sel_v[i][pi][li][pci] # 制約: 製品・ラインへの投資が行われたら生産可能状態になる make_c = qbpp.expr(shape=(n_years, n_prods, n_lines)) for yi in range(n_years): for pi in range(n_prods): for li in range(n_lines): # make_cum_v >= 1 のとき make_v = 1、それ以外は 0 M = n_years + 2 # upper bound value for make_acc_v make_c[yi][pi][li] = qbpp.constrain( (M * make_v[yi][pi][li]) - make_cum_v[yi][pi][li], between=(0, M) ) c_expr += qbpp.sum(make_c) ## HUBO未使用版終了 ## ################################# -
make_cum_v[yi][pi][li] > 0の場合、制約式(M × make_v[yi][pi][li]) - make_cum_v[yi][pi][li] ≥ 0を満たすためには、
make_v[yi][pi][li] = 1である必要があります。例えば、
make_cum_v[yi][pi][li]の最大値が 4 の場合 (３年間に生産強化を投資する)、M = 5 と設定すると、-
make_v = 0のとき0 - make_cum_v < 0となり制約違反になります。 -
make_v = 1のとき5 - make_cum_v ≥ 0が常に成立します。
したがって、
make_cum_v > 0の場合にはmake_v = 1が強制されます。なお、make_cum_v = 0の場合は、make_v = 0またはmake_v = 1のいずれも取り得ます。このときの値は、現在の製造可能状態を表すcurr_make_dataに依存します。 -
-
生産能力と投資額の計算
- 「設備投資を実施した結果、各年度でどれだけの生産能力が利用可能になるか」と「そのためにどれだけの投資コストが発生するか」を計算します。
-
prod_cap_v：各年度・各製品・各ラインにおいて利用可能な生産能力を表す変数です。過去年度に実施した能力増強を累積しています。 -
prod_inv_v：各年度に発生する製品別投資額を表しています。投資額の構成
投資額は以下2種類から構成されます。- 製造可能化投資： 生産できないラインを生産可能にするための投資
- 能力増強投資：既存設備の能力を増やす投資（選択した増強オプションに応じて投資額を計算します）
-
line_cap_v：各年度におけるライン全体の生産能力を表します。製品能力と同様に、過去の増強投資を累積します。 -
line_inv_v：各年度・ライン能力増強に必要な投資額を表します。選択された能力増強オプションに対応するコストを加算しています。
# 変数: 承認済み能力増強後の製品能力
prod_cap_v = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
prod_cap_v[yi][pi][li] += curr_prod_cap_data[pi][li]
for i in range(yi + 1):
for pci in range(n_prod_capups):
prod_cap_v[yi][pi][li] += (
prod_capup_sel_v[i][pi][li][pci] * prod_capup_data[pci]
)
# 変数: 製品・ライン能力増強の投資額
prod_inv_v = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
if yi == 0:
invest_new_cost = (
make_v[yi][pi][li] - curr_make_data[pi][li]
) * new_make_cost_data[pi][li]
else:
invest_new_cost = (
make_v[yi][pi][li] - make_v[yi - 1][pi][li]
) * new_make_cost_data[pi][li]
prod_inv_v[yi][pi][li] += invest_new_cost
for pci in range(n_prod_capups):
prod_inv_v[yi][pi][li] += (
prod_capup_sel_v[yi][pi][li][pci]
* prod_capup_cost_data[pi][li][pci]
)
# 変数: ライン能力増強後の能力
line_cap_v = qbpp.expr(shape=(n_years, n_lines))
for yi in range(n_years):
for li in range(n_lines):
line_cap_v[yi][li] += curr_line_cap_data[li]
for i in range(yi + 1):
for lci in range(n_line_capups):
line_cap_v[yi][li] += (
line_capup_sel_v[i][li][lci] * line_capup_data[lci]
)
# 変数: ライン能力増強投資額
line_inv_v = qbpp.expr(shape=(n_years, n_lines))
for yi in range(n_years):
for li in range(n_lines):
for lci in range(n_line_capups):
line_inv_v[yi][li] += (
line_capup_sel_v[yi][li][lci] * line_capup_cost_data[li][lci]
)
製品別総生産能力制約（prod_cap_sum_c）
- 各年度・各製品について、全ラインの生産能力の合計が需要量を満たしていることを保証する制約を定義しています。
- 製品ごとの需要を満たすために、十分な生産能力が確保されているかを確認する重要な制約です。
# 制約: 全ラインの製品能力合計が需要量を満たすこと
prod_cap_sum_v = qbpp.expr(shape=(n_years, n_prods))
prod_cap_sum_c = qbpp.expr(shape=(n_years, n_prods))
for yi in range(n_years):
for pi in range(n_prods):
ub_prod_cap_sum = 0
for li in range(n_lines):
ub_prod_cap_sum += curr_prod_cap_data[pi][li] + (yi + 1) * ub_prod_capup
prod_cap_sum_v[yi][pi] += prod_cap_v[yi][pi][li]
prod_cap_sum_c[yi][pi] = qbpp.constrain(
prod_cap_sum_v[yi][pi] - demand_data[yi][pi], between=(0, ub_prod_cap_sum)
)
c_expr += qbpp.sum(prod_cap_sum_c)
生産量変数と需要充足制約（prod_qty_v / prod_qty_sum_v / prod_qty_sum_c）
-
prod_qty_v：各年度において、各製品を各ラインで実際にどれだけ生産するかを表す意思決定変数です。- 本モデルの中で最終的な生産計画を表現する最も重要な変数の一つです。 -
prod_qty_sum_v：各製品について全ラインの生産量を合計した値を表します。 - 「各製品の総生産量が需要量と完全に一致する」ことを保証しています。
- 生産能力があるだけでは不十分で、実際に需要量分を生産する必要があります。そのため、本制約は生産計画を成立させる上で非常に重要な役割を持っています。
# 変数: 年度ごとの製品・ライン別生産量
prod_qty_v = qbpp.var(
"prod_qty", shape=(n_years, n_prods, n_lines), between=(0, ub_demand)
)
# 変数: 年度ごとの製品総生産量
prod_qty_sum_v = qbpp.expr(shape=(n_years, n_prods))
# 制約: 総生産量は需要量と一致すること
prod_qty_sum_c = qbpp.expr(shape=(n_years, n_prods))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
prod_qty_sum_v[yi][pi] += prod_qty_v[yi][pi][li]
prod_qty_sum_c[yi][pi] = qbpp.constrain(
prod_qty_sum_v[yi][pi] - demand_data[yi][pi], equal=0
)
c_expr += qbpp.sum(prod_qty_sum_c)
製品能力制約とライン能力制約（prod_cap_c / line_load_c）
- 生産計画では需要を満たすことが重要ですが、それ以上に重要なのは、
- 製品別の設備能力を超えないこと
- ライン全体の処理能力を超えないこと
# 制約: 製品・ラインごとの生産量は利用可能能力を超えないこと
prod_cap_c = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
ub_prod_cap_to_year = curr_prod_cap_data[pi][li] + (yi + 1) * ub_prod_capup
prod_cap_c[yi][pi][li] = qbpp.constrain(
prod_cap_v[yi][pi][li] - prod_qty_v[yi][pi][li],
between=(0, ub_prod_cap_to_year),
)
c_expr += qbpp.sum(prod_cap_c)
# 変数: 年度・ラインごとのラインの総負荷
line_load_v = qbpp.expr(shape=(n_years, n_lines))
# 制約: ラインの総負荷はライン能力を超えないこと
line_load_c = qbpp.expr(shape=(n_years, n_lines))
for yi in range(n_years):
for li in range(n_lines):
max_line_cap_to_year = curr_line_cap_data[li] + (yi + 1) * ub_line_capup
for pi in range(n_prods):
line_load_v[yi][li] += prod_qty_v[yi][pi][li]
line_load_c[yi][li] = qbpp.constrain(
line_cap_v[yi][li] - line_load_v[yi][li], between=(0, max_line_cap_to_year)
)
c_expr += qbpp.sum(line_load_c)
売上・投資額・利益の計算（revenue_v / invest_v / profit_v）
-
revenue_v：各年度において生産活動から得られる総売上を表します。計算方法は、各製品・各ラインについて、「生産量 × 単位利益」を計算し、合計しています。 -
invest_v：各年度で発生する総設備投資額を表します。投資額の構成は製品関連投資とライン関連投資を合計されます。 -
profit_v：各年度の最終利益を表します。「利益 = 売上 − 投資額」
# 変数: 各年度の売上
revenue_v = qbpp.expr(shape=(n_years))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
revenue_v[yi] += prod_qty_v[yi][pi][li] * unit_profit_data[pi][li] * 10
# 変数: 各年度の投資額
invest_v = qbpp.expr(shape=(n_years))
for yi in range(n_years):
# 製品能力増強投資
for pi in range(n_prods):
for li in range(n_lines):
invest_v[yi] += prod_inv_v[yi][pi][li]
# ライン能力増強投資
for li in range(n_lines):
invest_v[yi] += line_inv_v[yi][li]
# 変数: 各年度の利益
profit_v = qbpp.expr(shape=(n_years))
for yi in range(n_years):
profit_v[yi] += revenue_v[yi] - invest_v[yi]
目的関数（Objective Function）の定義
-
total_revenue_v：計画期間全体で得られる総売上を表します。 -
total_invest_v：計画期間全体で発生する設備投資額を表します。 - 「総利益 = 総売上 − 総投資額」を最大化するのは目的ですが、QUBOソルバーは基本的に最小化問題として問題を解きます。そのため、最大化問題を最小化問題へ変換する必要があります。
# 変数: 総売上利益
total_revenue_v = qbpp.expr()
for yi in range(n_years):
total_revenue_v += revenue_v[yi]
# 変数: 総投資額
total_invest_v = qbpp.expr()
for yi in range(n_years):
total_invest_v += invest_v[yi]
# 目的関数: 総利益最大化
o_expr += (-1) * total_revenue_v + total_invest_v
最終最適化モデルの構築（f）
- これまで作成した「制約条件」と「目的関数」を1つのQUBOモデルへ統合しています。
- QUBO（Quadratic Unconstrained Binary Optimization）では、本来の数理最適化問題に存在する制約条件を、ペナルティ項として目的関数へ組み込む必要があります。そのため、この部分はモデル全体の最終組み立てを行う重要な処理です。
- ペナルティとは制約違反を許してしまうため、制約違反に対して非常に大きな罰則を与える必要があります。
# 最終のQUBO式
PENALTY = 100000
f = c_expr * PENALTY + o_expr
f.simplify_as_binary()
最適化実行
- 構築したQUBOモデルをソルバーへ渡し、制約を満たしながら最も利益の高い生産・投資計画を探索しています。
- QUBO++の「EasySolver」または「ABS3Solver」で600秒の時間に最適解を探索します。
# solver = qbpp.EasySolver(f)
solver = qbpp.ABS3Solver(f)
sol = solver.search(time_limit=600, enable_default_callback=1)
出力結果
- ソルバーの計算ログ（モデルのエネルギー）
TTS = 0.000s Energy = 570800000
TTS = 0.000s Energy = 514799200
TTS = 0.000s Energy = 469998232
TTS = 0.000s Energy = 426797264
...
- 製品別結果：需要量、生産量、生産能力、能力増強量、投資額
- ライン別結果：ライン負荷、ライン能力、ライン増強量、投資額
- 年次業績
-------------------------------
年度: Y26
製品: P01, 需要量: 18, 生産量: 18, 生産能力: 21
ライン: L01, 生産/能力: 14/15, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 4/4, 能力増強量: 4 (新規), 投資額: 50
製品: P02, 需要量: 5, 生産量: 5, 生産能力: 14
ライン: L01, 生産/能力: 1/10, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 4/4, 能力増強量: 4 (新規), 投資額: 90
製品: P03, 需要量: 2, 生産量: 2, 生産能力: 9
ライン: L01, 生産/能力: 0/5, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 2/4, 能力増強量: 4 (新規), 投資額: 110
製品: P04, 需要量: 7, 生産量: 7, 生産能力: 21
ライン: L01, 生産/能力: 0/8, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 7/13, 能力増強量: 0, 投資額: 0
製品: P05, 需要量: 3, 生産量: 3, 生産能力: 17
ライン: L01, 生産/能力: 3/4, 能力増強量: 4 (新規), 投資額: 350
ライン: L02, 生産/能力: 0/13, 能力増強量: 0, 投資額: 0
製品: P06, 需要量: 0, 生産量: 0, 生産能力: 2
ライン: L01, 生産/能力: 0/0, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 0/2, 能力増強量: 2, 投資額: 30
---------
ライン: L01, 負荷/能力: 18/23, 能力増強量: 8, 投資額: 1000
ライン: L02, 負荷/能力: 17/20, 能力増強量: 0, 投資額: 0
-------------------------------
年度: Y27
製品: P01, 需要量: 0, 生産量: 0, 生産能力: 21
ライン: L01, 生産/能力: 0/15, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 0/4, 能力増強量: 0, 投資額: 0
製品: P02, 需要量: 8, 生産量: 8, 生産能力: 18
ライン: L01, 生産/能力: 0/10, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 8/8, 能力増強量: 4, 投資額: 40
製品: P03, 需要量: 3, 生産量: 3, 生産能力: 9
ライン: L01, 生産/能力: 0/5, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 3/4, 能力増強量: 0, 投資額: 0
製品: P04, 需要量: 7, 生産量: 7, 生産能力: 21
ライン: L01, 生産/能力: 0/8, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 7/13, 能力増強量: 0, 投資額: 0
製品: P05, 需要量: 3, 生産量: 3, 生産能力: 17
ライン: L01, 生産/能力: 3/4, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 0/13, 能力増強量: 0, 投資額: 0
製品: P06, 需要量: 0, 生産量: 0, 生産能力: 4
ライン: L01, 生産/能力: 0/0, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 0/4, 能力増強量: 2, 投資額: 30
---------
ライン: L01, 負荷/能力: 3/23, 能力増強量: 0, 投資額: 0
ライン: L02, 負荷/能力: 18/20, 能力増強量: 0, 投資額: 0
-------------------------------
年度: Y28
製品: P01, 需要量: 0, 生産量: 0, 生産能力: 21
ライン: L01, 生産/能力: 0/15, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 0/4, 能力増強量: 0, 投資額: 0
製品: P02, 需要量: 10, 生産量: 10, 生産能力: 20
ライン: L01, 生産/能力: 0/10, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 10/10, 能力増強量: 2, 投資額: 20
製品: P03, 需要量: 2, 生産量: 2, 生産能力: 9
ライン: L01, 生産/能力: 0/5, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 2/4, 能力増強量: 0, 投資額: 0
製品: P04, 需要量: 8, 生産量: 8, 生産能力: 21
ライン: L01, 生産/能力: 4/8, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 4/13, 能力増強量: 0, 投資額: 0
製品: P05, 需要量: 4, 生産量: 4, 生産能力: 17
ライン: L01, 生産/能力: 4/4, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 0/13, 能力増強量: 0, 投資額: 0
製品: P06, 需要量: 4, 生産量: 4, 生産能力: 4
ライン: L01, 生産/能力: 0/0, 能力増強量: 0, 投資額: 0
ライン: L02, 生産/能力: 4/4, 能力増強量: 0, 投資額: 0
---------
ライン: L01, 負荷/能力: 8/23, 能力増強量: 0, 投資額: 0
ライン: L02, 負荷/能力: 20/20, 能力増強量: 0, 投資額: 0
=============================
年度: Y26, 売上利益額: 41070, 投資額: 1630, 利益額: 39440
年度: Y27, 売上利益額: 24970, 投資額: 70, 利益額: 24900
年度: Y28, 売上利益額: 31240, 投資額: 20, 利益額: 31220
総売上利益額: 97280
総投資額: 1720
総利益額: 95560
結論
- 本プログラムでは、生産能力計画と設備投資計画を統合的に最適化するQUBOモデルを実装したものです。
- 従来は別々に検討されることの多かった「製造可能化投資」「製品能力増強」「ライン能力増強」「生産量配分」を単一モデルに統合し、総利益最大化という観点で最適解を探索しています。
- 特に重要なのは、投資効果を年度間で累積管理している点です。これにより短期的な利益だけでなく、中長期的な設備戦略や生産ネットワーク設計まで評価できるモデルになっています。
- さらにQUBO形式で実装されているため、量子アニーリングマシンやイジングマシンへの適用も容易であり、製造業における大規模な設備投資最適化問題への応用が期待されます。
- 今後は、在庫制約、リードタイム、複数工場、設備停止期間、投資予算制約などを追加することで、より実務に近い生産計画モデルへ発展させることが可能です。
参考
全プログラム
Show code
import pyqbpp as qbpp
##############################
### 問題データ ###
##############################
# 製品一覧
prod_list = ["P01", "P02", "P03", "P04", "P05", "P06"]
# 生産ライン一覧
line_list = ["L01", "L02"]
# 計画年度一覧
year_list = ["Y26", "Y27", "Y28"]
# 現在の製品・ライン適合マトリクス
# 1 = 当該ラインで生産可能, 0 = 生産不可
curr_make_data = [[1, 0], [1, 0], [1, 0], [1, 1], [0, 1], [0, 1]]
# 製品を生産可能にするための投資コスト（単位: 千万円）
# 0 = 既に当該ラインで生産可能
new_make_cost_data = [
[0, 20],
[0, 50],
[0, 100],
[0, 0],
[300, 0],
[1000, 0],
]
# 各製品・各ラインの現在の生産能力（単位: 万個）
curr_prod_cap_data = [
[15, 0],
[10, 0],
[5, 0],
[8, 13],
[0, 13],
[0, 0],
]
# 各ラインの現在の総生産能力（単位: 万個）
curr_line_cap_data = [15, 20]
# 製品能力増強オプション（単位: 万個）
# 0 = 能力増強なし
prod_capup_data = [0, 2, 4]
# 製品能力増強オプションごとのコスト（単位: 千万円）
prod_capup_cost_data = [
[[0, 10, 15], [0, 20, 30]],
[[0, 20, 40], [0, 20, 40]],
[[0, 10, 10], [0, 10, 10]],
[[0, 20, 60], [0, 20, 60]],
[[0, 10, 50], [0, 10, 50]],
[[0, 20, 40], [0, 30, 60]],
]
# ライン能力増強オプション（単位: 万個）
# 0 = 能力増強なし
line_capup_data = [0, 4, 8]
# ライン能力増強コスト（単位: 千万円）
line_capup_cost_data = [[0, 1000, 1000], [0, 1500, 1500]]
# 各年度・各製品の需要量（単位: 万個）
demand_data = [
[18, 5, 2, 7, 3, 0],
[0, 8, 3, 7, 3, 0],
[0, 10, 2, 8, 4, 4],
]
# 製品ごと・ラインごとの単位利益（単位: 万円）
unit_profit_data = [[100, 150], [80, 90], [100, 110], [120, 121], [200, 162], [50, 60]]
# 年間で選択可能な製品能力増強量の最大値
ub_prod_capup = max(prod_capup_data)
# 年間で選択可能なライン能力増強量の最大値
ub_line_capup = max(line_capup_data)
# 全年度・全製品における最大需要量
ub_demand = max(x for row in demand_data for x in row)
n_prods = len(prod_list)
n_lines = len(line_list)
n_prod_capups = len(prod_capup_data)
n_line_capups = len(line_capup_data)
n_years = len(year_list)
##############################
### QUBO問題の作成 ###
##############################
# 制約式全体を保持する式
c_expr = qbpp.Expr()
# 目的関数全体を保持する式
o_expr = qbpp.Expr()
# 変数: 年度・製品・ラインごとの能力増強オプション選択
prod_capup_sel_v = qbpp.var(
"p_up_sel_v", shape=(n_years, n_prods, n_lines, n_prod_capups)
)
# 制約: 各年度・製品・ラインごとに能力増強オプションを1つだけ選択する
prod_capup_sel_c = qbpp.constrain(qbpp.vector_sum(prod_capup_sel_v), equal=1)
c_expr += qbpp.sum(prod_capup_sel_c)
# 変数: 年度・ラインごとのライン能力増強オプション選択変数
line_capup_sel_v = qbpp.var("l_up_sel_v", shape=(n_years, n_lines, n_line_capups))
# 制約: 各年度・製品・ラインの組み合わせごとに能力増強オプションを必ず1つ選択する
line_capup_sel_c = qbpp.constrain(qbpp.vector_sum(line_capup_sel_v), equal=1)
c_expr += qbpp.sum(line_capup_sel_c)
# 変数: 指定年度において製品を当該ラインで生産可能かを表す変数
make_v = qbpp.var("make_v", shape=(n_years, n_prods, n_lines))
# 制約: 製造可能状態は 0→1 のみ変化可能で、1→0 は不可
# make_v[yi-1][pi][li] <= make[yi][pi][li]
make_change_c = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
if yi == 0:
# 初年度は初期状態と比較
make_change_c[yi][pi][li] = qbpp.constrain(
make_v[yi][pi][li] - curr_make_data[pi][li], between=(0, 1)
)
else:
# 2年目以降は前年度状態と比較
make_change_c[yi][pi][li] = qbpp.constrain(
make_v[yi][pi][li] - make_v[yi - 1][pi][li], between=(0, 1)
)
c_expr += qbpp.sum(make_change_c)
# 制約: 製品能力増強が選択された場合、その製品は生産可能状態になる
#################################
## QUBO++でHUBOを使用する場合 ##
make_c = qbpp.expr(shape=(n_years, n_prods, n_lines, n_prod_capups))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
for pci in range(1, n_prod_capups):
make_c[yi][pi][li][pci] = qbpp.constrain(
prod_capup_sel_v[yi][pi][li][pci] * (1 - make_v[yi][pi][li]),
equal=0,
)
c_expr += qbpp.sum(make_c)
## HUBO使用終了 ##
#################################
# ###############################
# ## QUBO++でHUBOを使用しない ##
# # 変数: 現年度まで累積した生産可能状態
# make_cum_v = qbpp.expr(shape=(n_years, n_prods, n_lines))
# for yi in range(n_years):
# for pi in range(n_prods):
# for li in range(n_lines):
# make_cum_v[yi][pi][li] += curr_make_data[pi][li]
# # 対象年度までの能力増強投資を累積
# for i in range(yi + 1):
# for pci in range(1, n_prod_capups):
# make_cum_v[yi][pi][li] += prod_capup_sel_v[i][pi][li][pci]
# # 制約: 製品・ラインへの投資が行われたら生産可能状態になる
# make_c = qbpp.expr(shape=(n_years, n_prods, n_lines))
# for yi in range(n_years):
# for pi in range(n_prods):
# for li in range(n_lines):
# # make_cum_v >= 1 のとき make_v = 1、それ以外は 0
# M = n_years + 2 # upper bound value for make_acc_v
# make_c[yi][pi][li] = qbpp.constrain(
# (M * make_v[yi][pi][li]) - make_cum_v[yi][pi][li], between=(0, M)
# )
# c_expr += qbpp.sum(make_c)
# ## HUBO未使用版終了 ##
# #################################
# 変数: 承認済み能力増強後の製品能力
prod_cap_v = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
prod_cap_v[yi][pi][li] += curr_prod_cap_data[pi][li]
for i in range(yi + 1):
for pci in range(n_prod_capups):
prod_cap_v[yi][pi][li] += (
prod_capup_sel_v[i][pi][li][pci] * prod_capup_data[pci]
)
# 変数: 製品・ライン能力増強の投資額
prod_inv_v = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
if yi == 0:
invest_new_cost = (
make_v[yi][pi][li] - curr_make_data[pi][li]
) * new_make_cost_data[pi][li]
else:
invest_new_cost = (
make_v[yi][pi][li] - make_v[yi - 1][pi][li]
) * new_make_cost_data[pi][li]
prod_inv_v[yi][pi][li] += invest_new_cost
for pci in range(n_prod_capups):
prod_inv_v[yi][pi][li] += (
prod_capup_sel_v[yi][pi][li][pci]
* prod_capup_cost_data[pi][li][pci]
)
# 変数: ライン能力増強後の能力
line_cap_v = qbpp.expr(shape=(n_years, n_lines))
for yi in range(n_years):
for li in range(n_lines):
line_cap_v[yi][li] += curr_line_cap_data[li]
for i in range(yi + 1):
for lci in range(n_line_capups):
line_cap_v[yi][li] += (
line_capup_sel_v[i][li][lci] * line_capup_data[lci]
)
# 変数: ライン能力増強投資額
line_inv_v = qbpp.expr(shape=(n_years, n_lines))
for yi in range(n_years):
for li in range(n_lines):
for lci in range(n_line_capups):
line_inv_v[yi][li] += (
line_capup_sel_v[yi][li][lci] * line_capup_cost_data[li][lci]
)
# 制約: 全ラインの製品能力合計が需要量を満たすこと
prod_cap_sum_v = qbpp.expr(shape=(n_years, n_prods))
prod_cap_sum_c = qbpp.expr(shape=(n_years, n_prods))
for yi in range(n_years):
for pi in range(n_prods):
ub_prod_cap_sum = 0
for li in range(n_lines):
ub_prod_cap_sum += curr_prod_cap_data[pi][li] + (yi + 1) * ub_prod_capup
prod_cap_sum_v[yi][pi] += prod_cap_v[yi][pi][li]
prod_cap_sum_c[yi][pi] = qbpp.constrain(
prod_cap_sum_v[yi][pi] - demand_data[yi][pi], between=(0, ub_prod_cap_sum)
)
c_expr += qbpp.sum(prod_cap_sum_c)
# 変数: 年度ごとの製品・ライン別生産量
prod_qty_v = qbpp.var(
"prod_qty", shape=(n_years, n_prods, n_lines), between=(0, ub_demand)
)
# 変数: 年度ごとの製品総生産量
prod_qty_sum_v = qbpp.expr(shape=(n_years, n_prods))
# 制約: 総生産量は需要量と一致すること
prod_qty_sum_c = qbpp.expr(shape=(n_years, n_prods))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
prod_qty_sum_v[yi][pi] += prod_qty_v[yi][pi][li]
prod_qty_sum_c[yi][pi] = qbpp.constrain(
prod_qty_sum_v[yi][pi] - demand_data[yi][pi], equal=0
)
c_expr += qbpp.sum(prod_qty_sum_c)
# 制約: 製品・ラインごとの生産量は利用可能能力を超えないこと
prod_cap_c = qbpp.expr(shape=(n_years, n_prods, n_lines))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
ub_prod_cap_to_year = curr_prod_cap_data[pi][li] + (yi + 1) * ub_prod_capup
prod_cap_c[yi][pi][li] = qbpp.constrain(
prod_cap_v[yi][pi][li] - prod_qty_v[yi][pi][li],
between=(0, ub_prod_cap_to_year),
)
c_expr += qbpp.sum(prod_cap_c)
# 変数: 年度・ラインごとのラインの総負荷
line_load_v = qbpp.expr(shape=(n_years, n_lines))
# 制約: ラインの総負荷はライン能力を超えないこと
line_load_c = qbpp.expr(shape=(n_years, n_lines))
for yi in range(n_years):
for li in range(n_lines):
max_line_cap_to_year = curr_line_cap_data[li] + (yi + 1) * ub_line_capup
for pi in range(n_prods):
line_load_v[yi][li] += prod_qty_v[yi][pi][li]
line_load_c[yi][li] = qbpp.constrain(
line_cap_v[yi][li] - line_load_v[yi][li], between=(0, max_line_cap_to_year)
)
c_expr += qbpp.sum(line_load_c)
# 変数: 各年度の売上利益
revenue_v = qbpp.expr(shape=(n_years))
for yi in range(n_years):
for pi in range(n_prods):
for li in range(n_lines):
revenue_v[yi] += prod_qty_v[yi][pi][li] * unit_profit_data[pi][li] * 10
# 変数: 各年度の投資額
invest_v = qbpp.expr(shape=(n_years))
for yi in range(n_years):
# 製品能力増強投資
for pi in range(n_prods):
for li in range(n_lines):
invest_v[yi] += prod_inv_v[yi][pi][li]
# ライン能力増強投資
for li in range(n_lines):
invest_v[yi] += line_inv_v[yi][li]
# 変数: 各年度の利益
profit_v = qbpp.expr(shape=(n_years))
for yi in range(n_years):
profit_v[yi] += revenue_v[yi] - invest_v[yi]
# 変数: 総売上利益
total_revenue_v = qbpp.expr()
for yi in range(n_years):
total_revenue_v += revenue_v[yi]
# 変数: 総投資額
total_invest_v = qbpp.expr()
for yi in range(n_years):
total_invest_v += invest_v[yi]
# 目的関数: 総利益最大化
o_expr += (-1) * total_revenue_v + total_invest_v
# 最終のQUBO式
PENALTY = 100000
f = c_expr * PENALTY + o_expr
f.simplify_as_binary()
###############################
### Solverで求解 ###
###############################
solver = qbpp.EasySolver(f)
# solver = qbpp.ABS3Solver(f)
sol = solver.search(time_limit=600, enable_default_callback=1)
###############################
### 解の表示 ###
###############################
for yi in range(n_years):
print("-------------------------------")
print(f"年度: {year_list[yi]}")
for pi in range(n_prods):
print(
f" 製品: {prod_list[pi]}, "
+ f"需要量: {demand_data[yi][pi]}, "
+ f"生産量: {sol(prod_qty_sum_v[yi][pi])}, "
+ f"生産能力: {sol(prod_cap_sum_v[yi][pi])}"
)
for li in range(n_lines):
capup = 0
for pci in range(n_prod_capups):
if sol(prod_capup_sel_v[yi][pi][li][pci]) == 1:
capup = prod_capup_data[pci]
if yi == 0:
last_makeable = curr_make_data[pi][li]
else:
last_makeable = sol(make_v[yi - 1][pi][li])
is_new_make = last_makeable != sol(make_v[yi][pi][li])
if last_makeable != sol(make_v[yi][pi][li]):
new_text = " (新規)"
else:
new_text = ""
print(
f" ライン: {line_list[li]}, "
+ f"生産/能力: {sol(prod_qty_v[yi][pi][li])}/{sol(prod_cap_v[yi][pi][li])}, "
+ f"能力増強量: {capup}{new_text}, "
+ f"投資額: {sol(prod_inv_v[yi][pi][li])}"
)
print("---------")
for li in range(n_lines):
l_capup = 0
for lci in range(n_line_capups):
if sol(line_capup_sel_v[yi][li][lci]) == 1:
l_capup = line_capup_data[lci]
print(
f" ライン: {line_list[li]}, "
+ f"負荷/能力: {sol(line_load_v[yi][li])}/{sol(line_cap_v[yi][li])}, "
+ f"能力増強量: {l_capup}, "
+ f"投資額: {sol(line_inv_v[yi][li])}"
)
print("=============================")
for yi in range(n_years):
print(
f"年度: {year_list[yi]}, 売上利益額: {sol(revenue_v[yi])}, 投資額: {sol(invest_v[yi])}, 利益額: {sol(profit_v[yi])}"
)
total_revenue = sol(total_revenue_v)
total_invest = sol(total_invest_v)
total_profit = total_revenue - total_invest
print(f"総売上利益額: {total_revenue}")
print(f"総投資額: {total_invest}")
print(f"総利益額: {total_profit}")