LLMで数理最適化のモデリング

OptiMUS,Claude,ChatGPTで数理最適化の定式化とコーディングを試みる

問題クッキーとケーキ

以下の問題を OptiMUS とClaude とChatGPTに聞いてみた．

クッキーとケーキを作るのに、次のように薄力粉とバターを使います。薄力粉は1200g、バターは900gしかありません。クッキーとケーキの合計重量を最大化するには、それぞれ何個作ればよいでしょうか？

クッキー1個は110gで、1gあたり薄力粉0.3gとバター0.3gが必要
ケーキ1個は200gで、1gあたり薄力粉0.4gとバター0.2gが必要

参照： Pythonで学ぶ数理最適化による問題解決入門（翔泳社）ビープラウド著

OptiMUSによるコード

OptiMUSは， https://optimus-solver.vercel.app/ で試せるマルチ・エージェント型のLLMである．

間違い． 1gあたりを読み違えている．

# https://optimus-solver.vercel.app/
# import json
# import numpy as np

# import gurobipy as gp

# with open("dummy_data.json", "r") as f:
#     data = json.load(f)


# PastryFlour = data["PastryFlour"]
# CookieWeight = data["CookieWeight"]
# CakeWeight = data["CakeWeight"]
# FlourPerGramCookie = data["FlourPerGramCookie"]
# ButterPerGramCake = data["ButterPerGramCake"]
# ButterPerGramCookie = data["ButterPerGramCookie"]
# FlourPerGramCake = data["FlourPerGramCake"]
# Butter = data["Butter"]

# # Define model
# model = gp.Model('model')


# # ====== Define variables ====== 
# ButterUsedForCookies = model.addVar(vtype=gp.GRB.CONTINUOUS, name="ButterUsedForCookies")
# FlourUsedForCakes = model.addVar(vtype=gp.GRB.CONTINUOUS, name="FlourUsedForCakes")
# ButterUsedForCakes = model.addVar(vtype=gp.GRB.CONTINUOUS, name="ButterUsedForCakes")
# CakesProduced = model.addVar(vtype=gp.GRB.INTEGER, name="CakesProduced")
# CookiesProduced = model.addVar(vtype=gp.GRB.INTEGER, name="CookiesProduced")

# # ====== Define constraints ====== 

# # Define the cookie production variable as an integer
# CookiesProduced = model.addVar(vtype=gp.GRB.INTEGER, name="CookiesProduced")

# # Define the integer variable for the number of cakes produced
# CakesProduced = model.addVar(vtype=gp.GRB.INTEGER, name="CakesProduced")

# # Add non-negative integer constraint for the number of cookies produced
# CookiesProduced = model.addVar(vtype=gp.GRB.INTEGER, name="CookiesProduced")
# model.addConstr(CookiesProduced >= 0, name="non_negative_cookies_produced")

# # Add non-negativity constraint for the number of cakes produced
# model.addConstr(CakesProduced >= 0, name="non_negativity_cakes")

# # Add constraint for the total flour usage for cookies and cakes
# model.addConstr(FlourPerGramCookie * CookieWeight * CookiesProduced + FlourPerGramCake * CakeWeight * CakesProduced <= PastryFlour, name="flour_capacity")

# # Add butter usage constraint
# model.addConstr(CakesProduced * CakeWeight * ButterPerGramCake + CookiesProduced * CookieWeight * ButterPerGramCookie <= Butter, name="butter_limit")

# # Constraint: The total flour used for cookies must not exceed the total available pastry flour
# model.addConstr(CookiesProduced * CookieWeight * FlourPerGramCookie <= PastryFlour, name="PastryFlourLimit")

# # Add constraint for butter used in cookies production
# ButterUsedForCookies = model.addVar(vtype=gp.GRB.CONTINUOUS, name="ButterUsedForCookies")
# CookiesProduced = model.addVar(vtype=gp.GRB.INTEGER, name="CookiesProduced")
# CookieWeight = 110
# ButterPerGramCookie = 0.3

# model.addConstr(ButterUsedForCookies == CookiesProduced * CookieWeight * ButterPerGramCookie, name="butter_usage")

# # Add flour usage constraint for cakes
# model.addConstr(FlourUsedForCakes == CakeWeight * FlourPerGramCake * CakesProduced, name="flour_usage_cakes")

# # Add constraint for total butter used in cakes
# model.addConstr(ButterUsedForCakes == CakeWeight * ButterPerGramCake * CakesProduced, name="butter_usage_cakes_constraint")

# # ====== Define objective ====== 

# # Set objective
# model.setObjective(CookieWeight * CookiesProduced + CakeWeight * CakesProduced, gp.GRB.MAXIMIZE)

# # Optimize model
# model.optimize()


# # Get model status
# status = model.status


# # Get solver information
# solving_info = {}
# if status == gp.GRB.OPTIMAL:
#     solving_info["status"] = model.status
#     solving_info["objective_value"] = model.objVal
#     solving_info["variables"] = [
#         {
#             "symbol": var.VarName,
#             "value": var.X,
#         }
#         for var in model.getVars()
#     ]
#     solving_info["runtime"] = model.Runtime
#     solving_info["iteration_count"] = model.IterCount
    



# # Get objective value
# obj_val = model.objVal

Claudeによるコード

同様の間違い．古いChatGPTもほぼ同様の間違えをする．

#Claudeによる解答 
#以下の問題を数理最適化で定式化してください．
#クッキーとケーキを作るのに、次のように薄力粉とバターを使います。薄力粉は1200g、バターは900gしかありません。クッキーとケーキの合計重量を最大化するには、それぞれ何個作ればよいでしょうか？ - クッキー1個は110gで、1gあたり薄力粉0.3gとバター0.3gが必要 - ケーキ1個は200gで、1gあたり薄力粉0.4gとバター0.2gが必要
#>>>
# この問題は線形計画問題として定式化できます。
# 【決定変数】

# x: クッキーの個数
# y: ケーキの個数

# 【目的関数】
# 最大化したい目的は、クッキーとケーキの合計重量です。
# 目的関数 = 110x + 200y
# 【制約条件】

# 薄力粉の使用量制約
# 0.3x + 0.4y ≤ 1200
# バターの使用量制約
# 0.3x + 0.2y ≤ 900
# 変数の非負条件
# x ≥ 0, y ≥ 0

# したがって、この問題は以下の線形計画問題として定式化できます。
# 【目的関数】
# 最大化: 110x + 200y
# 【制約条件】
# 0.3x + 0.4y ≤ 1200
# 0.3x + 0.2y ≤ 900
# x ≥ 0
# y ≥ 0
# この線形計画問題を解くことで、クッキーとケーキの合計重量を最大化する最適な個数を求めることができます。

from gurobipy import Model, GRB

# モデルの作成
model = Model("cookie_cake")

# 決定変数の定義
x = model.addVar(vtype=GRB.INTEGER, name="cookies")
y = model.addVar(vtype=GRB.INTEGER, name="cakes")

# 目的関数の設定
model.setObjective(110*x + 200*y, GRB.MAXIMIZE)

# 制約条件の設定
model.addConstr(0.3*x + 0.4*y <= 1200, "flour_constraint")
model.addConstr(0.3*x + 0.2*y <= 900, "butter_constraint")
model.addConstr(x >= 0, "x_non_negative")
model.addConstr(y >= 0, "y_non_negative")

# モデルの最適化
model.optimize()

# 最適解の出力
if model.status == GRB.OPTIMAL:
   print(f"最適解:")
   print(f"クッキーの個数: {int(x.x)}")
   print(f"ケーキの個数: {int(y.x)}")
   print(f"合計重量: {110*int(x.x) + 200*int(y.x)}")
else:
   print("最適解が見つかりませんでした。")

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 4 rows, 2 columns and 6 nonzeros
Model fingerprint: 0x6e1f517d
Variable types: 0 continuous, 2 integer (0 binary)
Coefficient statistics:
  Matrix range     [2e-01, 1e+00]
  Objective range  [1e+02, 2e+02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [9e+02, 1e+03]
Found heuristic solution: objective 330000.00000
Presolve removed 2 rows and 0 columns
Presolve time: 0.00s
Presolved: 2 rows, 2 columns, 4 nonzeros
Variable types: 0 continuous, 2 integer (0 binary)
Found heuristic solution: objective 330090.00000

Root relaxation: objective 6.000000e+05, 0 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0    600000.00000 600000.000  0.00%     -    0s

Explored 1 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 16 (of 16 available processors)

Solution count 3: 600000 330090 330000 

Optimal solution found (tolerance 1.00e-04)
Best objective 6.000000000000e+05, best bound 6.000000000000e+05, gap 0.0000%
最適解:
クッキーの個数: 0
ケーキの個数: 3000
合計重量: 600000

ChatGPT4o によるコード

この問題を線形計画問題として定式化するためには、まず変数を定義し、制約と目的関数を設定します。

変数：

\(x\): クッキーの個数
\(y\): ケーキの個数

クッキー1個は110gで、1gあたり薄力粉0.3gとバター0.3gが必要であるため、クッキー1個には以下の量の材料が必要です：

薄力粉：110g × 0.3 = 33g
バター：110g × 0.3 = 33g

ケーキ1個は200gで、1gあたり薄力粉0.4gとバター0.2gが必要であるため、ケーキ1個には以下の量の材料が必要です：

薄力粉：200g × 0.4 = 80g
バター：200g × 0.2 = 40g

したがって、以下のように定式化できます：

目的関数（クッキーとケーキの合計重量を最大化する）：

\[\text{maximize } 110x + 200y\]

制約条件：

\[33x + 80y \leq 1200 \quad (\text{薄力粉の制約})\] \[33x + 40y \leq 900 \quad (\text{バターの制約})\] \[x \geq 0\] \[y \geq 0\]

これで線形計画問題として定式化されました。次に、この問題を解くために単体法やその他の線形計画ソルバーを使用することができます。

手計算による解の概要

制約条件から計算を行います。

まず、制約を式の形式に揃えます：

\(33x + 80y \leq 1200\)
\(33x + 40y \leq 900\)

これらをグラフにプロットし、交点を求めます。二つの不等式の交点を計算します。

交点の計算：

\(33x + 80y = 1200\)
\(33x + 40y = 900\)

これらの式から\(y\) を消去します：

\[33x + 80y = 1200\] \[33x + 40y = 900\]

上の式から下の式を引くと： \[40y = 300\] \[y = 7.5\]

この値をどちらかの元の式に代入して \(x\) を求めます。ここでは 2 番目の式を使います：

\[33x + 40(7.5) = 900\] \[33x + 300 = 900\] \[33x = 600\] \[x = 18.18\]

交点の座標は (18.18, 7.5) です。この点が存在するかどうかをチェックします。

次に、この点での目的関数の値を計算します：

\[110(18.18) + 200(7.5)\]

これが最大値かどうかは、制約の範囲内の他の端点もチェックします：

\(x = 0\) の場合： \(33(0) + 80y \leq 1200 \rightarrow y = 15\)
\(y = 0\) の場合： \(33x + 80(0) \leq 1200 \rightarrow x = 36.36\)

これらのポイントに対しても目的関数の値を計算します。

以上の計算結果から、具体的な最適解が得られるでしょう。適切なソルバーを使うと、これを自動的に行うことができます。

コードを要求すると正解を出す．

from gurobipy import Model, GRB

# モデルの作成
model = Model("Cookie_and_Cake")

# 変数の追加
x = model.addVar(name="Cookies", vtype=GRB.CONTINUOUS)
y = model.addVar(name="Cakes", vtype=GRB.CONTINUOUS)

# 目的関数の設定
model.setObjective(110 * x + 200 * y, GRB.MAXIMIZE)

# 制約条件の追加
model.addConstr(33 * x + 80 * y <= 1200, "Flour")
model.addConstr(33 * x + 40 * y <= 900, "Butter")

# モデルの最適化
model.optimize()

# 最適解の取得
if model.status == GRB.OPTIMAL:
    print(f"Optimal number of cookies: {x.X}")
    print(f"Optimal number of cakes: {y.X}")
    print(f"Maximum total weight: {model.ObjVal}")
else:
    print("No optimal solution found")

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 2 rows, 2 columns and 4 nonzeros
Model fingerprint: 0x48bc11cf
Coefficient statistics:
  Matrix range     [3e+01, 8e+01]
  Objective range  [1e+02, 2e+02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [9e+02, 1e+03]
Presolve time: 0.00s
Presolved: 2 rows, 2 columns, 4 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    2.6250000e+31   1.968750e+30   2.625000e+01      0s
       2    3.5000000e+03   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.01 seconds (0.00 work units)
Optimal objective  3.500000000e+03
Optimal number of cookies: 18.181818181818183
Optimal number of cakes: 7.5
Maximum total weight: 3500.0

嘘つき島パズル

ある島には正直族と嘘つき族と呼ばれる2 種類の人たちが仲良く住んでいる．正直族は必ず本当のことを言い，嘘つき族は必ず嘘をつく．

あなたは，この島の人たちが正直族か嘘つき族なのかの調査を依頼された．

最初の家の旦那に聞いたところ「夫婦は両方とも嘘つき族だよ」という答えだった．

上の情報をもとに夫婦が嘘つき族かどうかを判別するための問題を解くための整数計画の定式化を示してください．

最初は数理最適化を使わずに論理のみで解くが，プロンプトで指示を追加すると正解を出す．

import gurobipy as gp
from gurobipy import GRB

# モデルの作成
model = gp.Model("truth_and_lie")

# 変数の作成
x_H = model.addVar(vtype=GRB.BINARY, name="x_H")
x_W = model.addVar(vtype=GRB.BINARY, name="x_W")

# 制約条件の追加
# 旦那が嘘つきである場合
model.addConstr((x_H == 0) >> (x_H + x_W >= 1), "c1")  # x_H = 0 なら x_H + x_W >= 1
# 旦那が正直である場合
model.addConstr((x_H == 1) >> (x_H + x_W == 0), "c2")  # x_H = 1 なら x_H + x_W == 0

# 目的関数の設定（この場合はダミー、制約を満たす解を見つけることが目的）
model.setObjective(0, GRB.MINIMIZE)

# 最適化の実行
model.optimize()

# 結果の表示
if model.status == GRB.OPTIMAL:
    print(f"x_H (旦那): {x_H.X}")
    print(f"x_W (妻): {x_W.X}")
else:
    print("最適解が見つかりませんでした。")

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 0 rows, 2 columns and 0 nonzeros
Model fingerprint: 0x2f9bab4b
Model has 2 general constraints
Variable types: 0 continuous, 2 integer (2 binary)
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  Objective range  [0e+00, 0e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [0e+00, 0e+00]
  GenCon rhs range [1e+00, 1e+00]
  GenCon coe range [1e+00, 1e+00]
Presolve removed 0 rows and 2 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 16 available processors)

Solution count 1: 0 

Optimal solution found (tolerance 1.00e-04)
Best objective 0.000000000000e+00, best bound 0.000000000000e+00, gap 0.0000%
x_H (旦那): 0.0
x_W (妻): 1.0

OR Tools (cp-sat)を指定すると，以下のコードを出す（正解）．

from ortools.sat.python import cp_model

# モデルの作成
model = cp_model.CpModel()

# 変数の作成
x_H = model.NewBoolVar('x_H')
x_W = model.NewBoolVar('x_W')

# 制約条件の追加
# 旦那が嘘つきである場合
model.Add(x_H == 0).OnlyEnforceIf(x_H.Not())
model.Add(x_H + x_W >= 1).OnlyEnforceIf(x_H.Not())

# 旦那が正直である場合
model.Add(x_H == 1).OnlyEnforceIf(x_H)
model.Add(x_H + x_W == 0).OnlyEnforceIf(x_H)

# ソルバーの作成と解の探索
solver = cp_model.CpSolver()
status = solver.Solve(model)

# 結果の表示
if status == cp_model.OPTIMAL:
    print(f"x_H (旦那): {solver.Value(x_H)}")
    print(f"x_W (妻): {solver.Value(x_W)}")
else:
    print("最適解が見つかりませんでした。")

x_H (旦那): 0
x_W (妻): 1

OptiMUSの例題

https://github.com/teshnizi/OptiMUS で配布しているデータを入力して問題を解いてみる．

Linear and Convex Optimization: A Mathematical Approach, Michael H. Veatch
Model Building in Mathematical Programming, H. Paul Williams

データは流通情報工学演習のページで配布している．

データセットを解凍すると， 3つの本の練習問題ごとに問題ファイルが出てくる．各フォルダには，以下のテキストファイルがある．

description.txt: 問題の概要を記述したテキストファイル
data.json: サンプルデータのJSONファイル
obj.txt: 最適目的関数値
output.json: 出力例（最適解）

生産計画

ある会社が，3つの製品A,B,Cを製造している。製造工程には、「切断」と「プレス」という2つの工程がある．すべての製品はこの2つの瞬間を通過しなければならない。

1日8時間使用可能な切削部門の生産能力は，製品Aは1時間に2000個、製品Bは1時間に1600個，製品Cは1時間に1100個である．

1日8時間使用可能なプレス部門の生産能力は，製品Aは1時間に1000個、製品Bは1時間に1500個、製品Cは1時間に2400個である．

製品切り替えの段取り時間（費用）はかからないものとする．

製品A,B,Cの1個あたりの利益をそれぞれ12,9,8（万円）としたとき，利益を最大化する製造計画を立案せよ．

この問題は、製品A、B、Cの製造量を決定して、切断とプレスの生産能力制約を満たしつつ、総利益を最大化することを目的とした線形計画問題です。Gurobiを用いて解くためのPythonコードを示します。

問題の定式化

変数

\(x_A\) : 製品Aの製造量
\(x_B\) : 製品Bの製造量
\(x_C\) : 製品Cの製造量

制約

切断部門の生産能力制約 \[\frac{x_A}{2000} + \frac{x_B}{1600} + \frac{x_C}{1100} \leq 8\]
プレス部門の生産能力制約 \[\frac{x_A}{1000} + \frac{x_B}{1500} + \frac{x_C}{2400} \leq 8\]

目的関数

総利益の最大化
\[12x_A + 9x_B + 8x_C\]

コードの説明

パラメータ:
- profits: 各製品の利益
- cutting_capacity: 切断部門の生産能力
- pressing_capacity: プレス部門の生産能力
- working_hours: 1日の労働時間
変数の追加:
- x['A'], x['B'], x['C']: 製品A, B, Cの製造量
目的関数の設定:
- 各製品の利益の合計を最大化
制約の追加:
- 切断部門の制約
- プレス部門の制約
モデルの最適化:
- model.optimize(): 最適解を求める

このコードを実行すると、Gurobiが最適な製造量を計算し、最大の利益をもたらす製造計画を示します。

from gurobipy import Model, GRB, quicksum

# モデルの作成
model = Model("Production_Planning")

# パラメータ
profits = {'A': 12, 'B': 9, 'C': 8}
cutting_capacity = {'A': 2000, 'B': 1600, 'C': 1100}
pressing_capacity = {'A': 1000, 'B': 1500, 'C': 2400}
working_hours = 8

# 変数の追加
x = model.addVars(['A', 'B', 'C'], vtype=GRB.CONTINUOUS, name="x")

# 目的関数の設定（利益の最大化）
model.setObjective(quicksum(profits[i] * x[i] for i in ['A', 'B', 'C']), GRB.MAXIMIZE)

# 切断部門の制約
model.addConstr(quicksum(x[i] / cutting_capacity[i] for i in ['A', 'B', 'C']) <= working_hours, "Cutting")

# プレス部門の制約
model.addConstr(quicksum(x[i] / pressing_capacity[i] for i in ['A', 'B', 'C']) <= working_hours, "Pressing")

# モデルの最適化
model.optimize()

# 結果の表示
if model.status == GRB.OPTIMAL:
    print("Optimal solution found:")
    for i in ['A', 'B', 'C']:
        print(f"  Product {i}: {x[i].X} units")
    print(f"Total profit: {model.ObjVal} (in 10,000 yen)")
else:
    print("No optimal solution found")

Set parameter Username
Academic license - for non-commercial use only - expires 2024-11-13
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 2 rows, 3 columns and 6 nonzeros
Model fingerprint: 0x37e2d5ea
Coefficient statistics:
  Matrix range     [4e-04, 1e-03]
  Objective range  [8e+00, 1e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [8e+00, 8e+00]
Presolve time: 0.01s
Presolved: 2 rows, 3 columns, 6 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    9.2800000e+32   3.149576e+30   9.280000e+02      0s
       2    1.1312432e+05   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.01 seconds (0.00 work units)
Optimal objective  1.131243243e+05
Optimal solution found:
  Product A: 5621.621621621622 units
  Product B: 0.0 units
  Product C: 5708.108108108108 units
Total profit: 113124.32432432433 (in 10,000 yen)

スポーツ用品の生産

linear_and_convex_optimization (problem 1)

Wild Sports produces \(M\) different products using \(N\) different raw materials.
The company has \(available_{i}\) of raw material \(i\) available.
Product requires units of material per unit produced.
Product has a selling price of per unit.
Product has a production cost of per unit produced.
Product has a forecasted maximum demand of for the month.

INPUT FORMAT:

{
    "available": [available_{i} for i in 1, ..., N],
    "requirements": [[req_{i,j} for i in 1, ..., N] for j in 1, ..., M],
    "prices": [price_{j} for j in 1, ..., M],
    "costs": [cost_{j} for j in 1, ..., M],
    "demands": [demand_{j} for j in 1, ..., M]
}

OBJECTIVE: How many of each product should the sports produce in order to maximize its profit for the month?

OUTPUT INFO:

the amount of each product produced .
shows the total profit of the company for the month.

OUTPUT FORMAT:

{
    "amount": [amount_{j} for j in 1, ..., M],
    "total_profit": total_profit
}

以下の問題を数理最適化で定式化します。

問題の定式化

変数

\(x_j\) : 製品 \(j\) の生産量（単位）

パラメータ

\(\text{available}_i\) : 原材料\(i\) の利用可能量
\(\text{req}_{i,j}\) : 製品\(j\) の生産に必要な原材料\(i\) の単位数
\(\text{price}_j\) : 製品\(j\) の販売価格（単位）
\(\text{cost}_j\) : 製品\(j\) の生産コスト（単位）
\(\text{demand}_j\) : 製品\(j\) の予測最大需要量

目的関数

会社の月間利益を最大化する。

\[ \text{Maximize } \sum_{j=1}^M (\text{price}_j - \text{cost}_j) x_j \]

制約条件

原材料の利用可能量制約:

\[ \sum_{j=1}^M \text{req}_{i,j} x_j \leq \text{available}_i \quad \forall i = 1, 2, \ldots, N \]

予測最大需要量制約: \[ x_j \leq \text{demand}_j \quad \forall j = 1, 2, \ldots, M \]
非負制約: \[ x_j \geq 0 \quad \forall j = 1, 2, \ldots, M \]

入力形式

{
    "available": [available_i for i in 1, ..., N],
    "requirements": [[req_i_j for i in 1, ..., N] for j in 1, ..., M],
    "prices": [price_j for j in 1, ..., M],
    "costs": [cost_j for j in 1, ..., M],
    "demands": [demand_j for j in 1, ..., M]
}

出力形式

{
    "amount": [amount_j for j in 1, ..., M],
    "total_profit": total_profit
}

このモデルでは、各製品の生産量を決定し、利益を最大化するための生産計画を最適化します。結果として、生産量の配分と総利益を出力します。

from gurobipy import Model, GRB

# 入力データの例
available = [100, 150, 200]  # 利用可能な原材料の例
requirements = [
    [2, 1, 3],  # 製品1に必要な各原材料の単位
    [1, 2, 1],  # 製品2に必要な各原材料の単位
    [3, 2, 2]   # 製品3に必要な各原材料の単位
]
prices = [20, 30, 40]  # 製品の販売価格の例
costs = [10, 15, 20]   # 製品の生産コストの例
demands = [50, 60, 70] # 製品の最大需要量の例

# 定数
N = len(available)  # 原材料の数
M = len(prices)     # 製品の数

# モデルの作成
model = Model("WildSports")

# 変数の作成
x = model.addVars(M, name="x", lb=0, vtype=GRB.CONTINUOUS)

# 目的関数の設定
profit = sum((prices[j] - costs[j]) * x[j] for j in range(M))
model.setObjective(profit, GRB.MAXIMIZE)

# 制約条件の設定
# 原材料の利用可能量制約
for i in range(N):
    model.addConstr(sum(requirements[j][i] * x[j] for j in range(M)) <= available[i], name=f"raw_material_{i}")

# 需要制約
for j in range(M):
    model.addConstr(x[j] <= demands[j], name=f"demand_{j}")

# 最適化の実行
model.optimize()

# 結果の出力
amount = [x[j].X for j in range(M)]
total_profit = model.ObjVal

# 出力形式に合わせて結果を表示
output = {
    "amount": amount,
    "total_profit": total_profit
}
print(output)

Set parameter Username
Academic license - for non-commercial use only - expires 2024-11-13
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 6 rows, 3 columns and 12 nonzeros
Model fingerprint: 0x57c4962d
Coefficient statistics:
  Matrix range     [1e+00, 3e+00]
  Objective range  [1e+01, 2e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e+01, 2e+02]
Presolve removed 3 rows and 0 columns
Presolve time: 0.01s
Presolved: 3 rows, 3 columns, 9 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    1.5000000e+03   1.625000e+01   0.000000e+00      0s
       1    1.1666667e+03   0.000000e+00   0.000000e+00      0s

Solved in 1 iterations and 0.01 seconds (0.00 work units)
Optimal objective  1.166666667e+03
{'amount': [0.0, 60.0, 13.333333333333334], 'total_profit': 1166.6666666666667}

食品製造問題1

model_building_in_mathematical_programming (problem 1)

定式化の詳細

パラメータ

\(I\): 油の種類数
\(M\): 月数
\(\text{price}_{i,m}\): 月\(m\) に油\(i\) を購入する価格
\(\text{sell\_price}\): 最終製品の販売価格（トン当たり）
\(\text{is\_vegetable}_{i}\): 油\(i\) が植物油であれば1、非植物油であれば0
\(\text{max\_vegetable\_refining\_per\_month}\): 月に精製できる植物油の最大トン数
\(\text{max\_non\_vegetable\_refining\_per\_month}\): 月に精製できる非植物油の最大トン数
\(\text{storage\_size}\): 各油の最大貯蔵量（トン）
\(\text{storage\_cost}\): 貯蔵コスト（トンあたり月額）
\(\text{max\_hardness}\): 最終製品の最大硬度
\(\text{min\_hardness}\): 最終製品の最小硬度
\(\text{hardness}_{i}\): 油\(i\) の硬度
\(\text{init\_amount}\): 各油の初期在庫量（トン）

変数

\(\text{buyquantity}_{i,m}\): 月\(m\) に油\(i\) を購入する量（トン）
\(\text{refine}_{i,m}\): 月\(m\) に油\(i\) を精製する量（トン）
\(\text{storage}_{i,m}\): 月\(m\) の終わりに油\(i\) を貯蔵している量（トン）

目的関数

利益の最大化を目指します。利益は総収入から総コスト（購入コストと貯蔵コスト）を引いたものです。

\[ \text{Maximize} \quad \sum_{m=1}^{M} \sum_{i=1}^{I} (\text{sell\_price} \cdot \text{refine}_{i,m} - \text{price}_{i,m} \cdot \text{buyquantity}_{i,m} - \text{storage\_cost} \cdot \text{storage}_{i,m}) \]

制約条件

在庫バランス：

\[ \text{storage}_{i,m} = \text{storage}_{i,m-1} + \text{buyquantity}_{i,m} - \text{refine}_{i,m} \quad \forall i, \forall m \]

精製制約：

\[ \sum_{i: \text{is\_vegetable}_i = 1} \text{refine}_{i,m} \leq \text{max\_vegetable\_refining\_per\_month} \quad \forall m \]

\[ \sum_{i: \text{is\_vegetable}_i = 0} \text{refine}_{i,m} \leq \text{max\_non\_vegetable\_refining\_per\_month} \quad \forall m \]

貯蔵制約：

\[\text{storage}_{i,m} \leq \text{storage\_size} \quad \forall i, \forall m\]

硬度制約： \[ \text{min\_hardness} \leq \frac{\sum_{i=1}^{I} \text{hardness}_i \cdot \text{refine}_{i,m}}{\sum_{i=1}^{I} \text{refine}_{i,m}} \leq \text{max\_hardness} \quad \forall m \]
初期在庫制約： \[ \text{storage}_{i,0} = \text{init\_amount} \quad \forall i \]
最終在庫制約：

\[ \text{storage}_{i,M} = \text{init\_amount} \quad \forall i \]

説明

データの読み込み：
- 提供されたデータを読み込みます。
変数の定義：
- buyquantity：各油を各月に購入する量。
- refine：各

結果

入力データの月と製品を逆に入力しているため，月数と製品数が異なる場合にはエラー（同じ場合には動いてしまうので，取りにくいバグ）
非線形の制約になるので， Gurobiでは解けずに実行エラー

修正プロンプトと返答

修正プロンプト

硬度制約が線形表現での除算で表されているのでエラーするようです．線形式に変換してください．

硬度制約を線形形式に変換するために、以下のように線形化します。これは、分母を掛け算にして除算を避ける方法です。

\[ \text{min\_hardness} \cdot \sum_{i=1}^{I} \text{refine}_{i,m} \leq \sum_{i=1}^{I} \text{hardness}_{i} \cdot \text{refine}_{i,m} \leq \text{max\_hardness} \cdot \sum_{i=1}^{I} \text{refine}_{i,m} \quad \forall m \]

これで硬度制約が線形形式で表現され、最適化問題を解くことができます。

硬度制約を線形形式に変換しました。これは、分母の総精製量を両辺に掛けることで除算を避けています。
それ以外の部分は元のコードと同様です。

import gurobipy as gp
from gurobipy import GRB

def optimize_food_production(data):
    I = len(data['buy_price'][0])  # Number of oil types
    M = len(data['buy_price'])  # Number of months
    
    # Parameters
    buy_price = data['buy_price']
    sell_price = data['sell_price']
    is_vegetable = data['is_vegetable']
    max_vegetable_refining_per_month = data['max_vegetable_refining_per_month']
    max_non_vegetable_refining_per_month = data['max_non_vegetable_refining_per_month']
    storage_size = data['storage_size']
    storage_cost = data['storage_cost']
    max_hardness = data['max_hardness']
    min_hardness = data['min_hardness']
    hardness = data['hardness']
    init_amount = data['init_amount']
    
    # Create a new model
    model = gp.Model("food_production")
    
    # Decision variables
    buy = model.addVars(I, M, name="buy", vtype=GRB.CONTINUOUS)
    refine = model.addVars(I, M, name="refine", vtype=GRB.CONTINUOUS)
    storage = model.addVars(I, M, name="storage", vtype=GRB.CONTINUOUS)
    
    # Auxiliary variables for total refined amount and total hardness
    total_refined = model.addVars(M, name="total_refined", vtype=GRB.CONTINUOUS)
    total_hardness = model.addVars(M, name="total_hardness", vtype=GRB.CONTINUOUS)
    
    # Objective function: Maximize profit
    revenue = gp.quicksum(refine[i, m] * sell_price for i in range(I) for m in range(M))
    purchase_cost = gp.quicksum(buy[i, m] * buy_price[m][i] for i in range(I) for m in range(M))
    storage_cost_total = gp.quicksum(storage[i, m] * storage_cost for i in range(I) for m in range(M))
    model.setObjective(revenue - purchase_cost - storage_cost_total, GRB.MAXIMIZE)
    
    # Constraints
    
    # Initial storage
    for i in range(I):
        model.addConstr(storage[i, 0] == init_amount)
    
    # Storage capacity constraints
    for m in range(M):
        for i in range(I):
            model.addConstr(storage[i, m] <= storage_size)
    
    # Refining constraints
    for m in range(M):
        model.addConstr(gp.quicksum(refine[i, m] for i in range(I) if is_vegetable[i] == 1) <= max_vegetable_refining_per_month)
        model.addConstr(gp.quicksum(refine[i, m] for i in range(I) if is_vegetable[i] == 0) <= max_non_vegetable_refining_per_month)
    
    # Total refined amount and hardness constraints
    for m in range(M):
        model.addConstr(total_refined[m] == gp.quicksum(refine[i, m] for i in range(I)))
        model.addConstr(total_hardness[m] == gp.quicksum(refine[i, m] * hardness[i] for i in range(I)))
        model.addConstr(total_hardness[m] >= total_refined[m] * min_hardness)
        model.addConstr(total_hardness[m] <= total_refined[m] * max_hardness)
    
    # Storage dynamics
    for i in range(I):
        for m in range(1, M):
            model.addConstr(storage[i, m] == storage[i, m-1] + buy[i, m] - refine[i, m])
    
    # Final storage must be equal to initial amount
    for i in range(I):
        model.addConstr(storage[i, M-1] == init_amount)
    
    # Optimize model
    model.optimize()
    
    # Prepare the result
    result = {
        "buy": [[buy[i, m].x for i in range(I)] for m in range(M)],
        "refine": [[refine[i, m].x for i in range(I)] for m in range(M)],
        "storage": [[storage[i, m].x for i in range(I)] for m in range(M)]
    }
    
    return result

# Example input
{
    "buy_price":[
        [110, 120, 130, 110, 115],
        [130, 130, 110, 90, 115],
        [110, 140, 130, 100, 95],
        [120, 110, 120, 120, 125],
        [100, 120, 150, 110, 105],
        [90, 100, 140, 80, 135]
    ],
    "sell_price": 150,
    "is_vegetable": [1, 1, 0, 0, 0],
    "max_vegetable_refining_per_month": 200,
    "max_non_vegetable_refining_per_month": 250,
    "storage_size": 1000,
    "storage_cost": 5,
    "min_hardness": 3,
    "max_hardness": 6,
    "hardness": [8.8, 6.1, 2.0, 4.2, 5.0],
    "init_amount": 500
}

# Call the optimization function
result = optimize_food_production(data)
print(result)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 32 rows, 24 columns and 62 nonzeros
Model fingerprint: 0x4a01e00d
Coefficient statistics:
  Matrix range     [1e+00, 1e+01]
  Objective range  [5e+00, 2e+02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [2e+02, 1e+03]
Presolve removed 32 rows and 24 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    3.1650000e+05   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.01 seconds (0.00 work units)
Optimal objective  3.165000000e+05
{'buy': [[0.0, 0.0], [1000.0, 200.0], [0.0, 600.0]], 'refine': [[500.0, 400.0], [500.0, 400.0], [500.0, 400.0]], 'storage': [[200.0, 200.0], [700.0, 0.0], [200.0, 200.0]]}

食品製造問題2

model_building_in_mathematical_programming (problem 2)

この問題をGurobiで解くために、混合整数最適化問題として定式化します。以下のGurobiモデルを用いて、入力データに基づいて最適な購入、精製、および保管の方針を求めます。

モデルの定式化

追加パラメータ

\(\text{min\_usage}\)：オイルの最小使用量
\(\text{dependencies}_{i}\)：オイル \(i\) を使用する際に必要な依存関係

追加変数

\(\text{use}_{i,m}\)：月 \(m\) にオイル \(i\) を使用するかどうか（0または1のバイナリ変数）

制約条件

各オイルの月ごとの保管量は、保管容量を超えないこと
各月の植物油の精製量は最大精製量を超えないこと
各月の非植物油の精製量は最大精製量を超えないこと
製品の硬度は、指定された範囲内に収まること
各オイルの使用量が、使用する場合には最小使用量以上であること
各オイルの使用に対する依存関係を満たすこと
各月に使用されるオイルの種類数が3種類以下であること
最終月の保管量が初期保管量と同じであること

import gurobipy as gp
from gurobipy import GRB

def optimize_food_production(data):
    I = len(data['buy_price'][0])  # Number of oil types
    M = len(data['buy_price'])  # Number of months
    
    # Parameters
    buy_price = data['buy_price']
    sell_price = data['sell_price']
    is_vegetable = data['is_vegetable']
    max_vegetable_refining_per_month = data['max_vegetable_refining_per_month']
    max_non_vegetable_refining_per_month = data['max_non_vegetable_refining_per_month']
    storage_size = data['storage_size']
    storage_cost = data['storage_cost']
    max_hardness = data['max_hardness']
    min_hardness = data['min_hardness']
    hardness = data['hardness']
    init_amount = data['init_amount']
    min_usage = data['min_usage']
    dependencies = data['dependencies']
    
    # Create a new model
    model = gp.Model("food_production")
    
    # Decision variables
    buy = model.addVars(I, M, name="buy", vtype=GRB.CONTINUOUS)
    refine = model.addVars(I, M, name="refine", vtype=GRB.CONTINUOUS)
    storage = model.addVars(I, M, name="storage", vtype=GRB.CONTINUOUS)
    use = model.addVars(I, M, name="use", vtype=GRB.BINARY)
    
    # Auxiliary variables for total refined amount and total hardness
    total_refined = model.addVars(M, name="total_refined", vtype=GRB.CONTINUOUS)
    total_hardness = model.addVars(M, name="total_hardness", vtype=GRB.CONTINUOUS)
    
    # Objective function: Maximize profit
    revenue = gp.quicksum(refine[i, m] * sell_price for i in range(I) for m in range(M))
    purchase_cost = gp.quicksum(buy[i, m] * buy_price[m][i] for i in range(I) for m in range(M))
    storage_cost_total = gp.quicksum(storage[i, m] * storage_cost for i in range(I) for m in range(M))
    model.setObjective(revenue - purchase_cost - storage_cost_total, GRB.MAXIMIZE)
    
    # Constraints
    
    # Initial storage
    for i in range(I):
        model.addConstr(storage[i, 0] == init_amount)
    
    # Storage capacity constraints
    for m in range(M):
        for i in range(I):
            model.addConstr(storage[i, m] <= storage_size)
    
    # Refining constraints
    for m in range(M):
        model.addConstr(gp.quicksum(refine[i, m] for i in range(I) if is_vegetable[i] == 1) <= max_vegetable_refining_per_month)
        model.addConstr(gp.quicksum(refine[i, m] for i in range(I) if is_vegetable[i] == 0) <= max_non_vegetable_refining_per_month)
    
    # Total refined amount and hardness constraints
    for m in range(M):
        model.addConstr(total_refined[m] == gp.quicksum(refine[i, m] for i in range(I)))
        model.addConstr(total_hardness[m] == gp.quicksum(refine[i, m] * hardness[i] for i in range(I)))
        model.addConstr(total_hardness[m] >= total_refined[m] * min_hardness)
        model.addConstr(total_hardness[m] <= total_refined[m] * max_hardness)
    
    # Minimum usage constraints
    for i in range(I):
        for m in range(M):
            model.addConstr(refine[i, m] >= min_usage * use[i, m])
    
    # Dependency constraints
    for i in range(I):
        for j in range(I):
            if dependencies[i][j] == 1:
                for m in range(M):
                    model.addConstr(use[i, m] <= use[j, m])
    
    # Limit number of oils used in each month
    for m in range(M):
        model.addConstr(gp.quicksum(use[i, m] for i in range(I)) <= 3)
    
    # Storage dynamics
    for i in range(I):
        for m in range(1, M):
            model.addConstr(storage[i, m] == storage[i, m-1] + buy[i, m] - refine[i, m])
    
    # Final storage must be equal to initial amount
    for i in range(I):
        model.addConstr(storage[i, M-1] == init_amount)
    
    # Optimize model
    model.optimize()
    
    # Prepare the result
    result = {
        "buy": [[buy[i, m].x for i in range(I)] for m in range(M)],
        "refine": [[refine[i, m].x for i in range(I)] for m in range(M)],
        "storage": [[storage[i, m].x for i in range(I)] for m in range(M)]
    }
    
    return result

# Example input
{
    "buy_price":[
        [110, 120, 130, 110, 115],
        [130, 130, 110, 90, 115],
        [110, 140, 130, 100, 95],
        [120, 110, 120, 120, 125],
        [100, 120, 150, 110, 105],
        [90, 100, 140, 80, 135]
    ],
    "sell_price": 150,
    "is_vegetable": [1, 1, 0, 0, 0],
    "max_vegetable_refining_per_month": 200,
    "max_non_vegetable_refining_per_month": 250,
    "storage_size": 1000,
    "storage_cost": 5,
    "min_hardness": 3,
    "max_hardness": 6,
    "hardness": [8.8, 6.1, 2.0, 4.2, 5.0],
    "init_amount": 500,
    "min_usage": 20,
    "dependencies": [
        [0, 0, 0, 0, 1],
        [0, 0, 0, 0, 1],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0]
    ]
}

# Call the optimization function
result = optimize_food_production(data)
print(result)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 47 rows, 30 columns and 92 nonzeros
Model fingerprint: 0x42107314
Variable types: 24 continuous, 6 integer (6 binary)
Coefficient statistics:
  Matrix range     [1e+00, 5e+01]
  Objective range  [5e+00, 2e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [3e+00, 1e+03]
Presolve removed 47 rows and 30 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 16 available processors)

Solution count 1: 316500 

Optimal solution found (tolerance 1.00e-04)
Best objective 3.165000000000e+05, best bound 3.165000000000e+05, gap 0.0000%
{'buy': [[0.0, 0.0], [1000.0, 200.0], [0.0, 600.0]], 'refine': [[500.0, 400.0], [500.0, 400.0], [500.0, 400.0]], 'storage': [[200.0, 200.0], [700.0, 0.0], [200.0, 200.0]]}

工場計画問題1

model_building_in_mathematical_programming (problem 3)

エンジニアリング工場は複数の製品を機械で製造し、工場が保有する機械 \(m\) の数は \(num_{m}\)
各製品 \(k\) は利益 \(profit_{k}\) をもたらす（単位あたりの販売価格から原材料費を差し引いたものとして定義）
製品 \(k\) が機械 \(m\) で必要とする単位生産時間（時間）は \(time_{k, m}\)
現在の月（1月）およびその後の数か月間、特定の機械はメンテナンスのために停止する予定
月 \(i\) において、停止する機械 \(m\) の数は \(maintain_{i, m}\)
各製品に対して各月ごとにマーケティング上の制限が存在する
月 \(i\) における製品 \(k\) の制限は \(limit_{k, i}\)
各製品を最大100単位まで保管することが可能であり、保管コストは月あたり単位 \(store_price\)
現在は在庫が存在しないが、今後数か月の終わりには各製品の在庫を \(keep_quantity\) 保持する必要がある
工場は週6日、1日 \(n\_workhours\) 時間稼働している
シーケンシング問題を考慮する必要はない
便宜上、1か月は24日であると仮定できる

入力形式:

{
    "num_machines": [num_{m} for m in 1,...,M],
    "profit": [profit_{k} for k in 1,...K],
    "time": [[time_{k, m} for m in 1,...,M] for k in 1,...,K],
    "maintain": [[maintain_{i, m} for m in 1,...,M] for i in 1,...,I],
    "limit": [[limit_{k, i} for i in 1,...,I] for k in 1,...,K],
    "store_price": store_price,
    "keep_quantity": keep_quantity,
    "n_workhours": n_workhours
}

目的:

会社が利益を最大化するために追求すべき販売、保管および製造の方針は何か

出力情報:

\(sell\) は各製品の各月における販売数量を表す
\(manufacture\) は各月における製造された製品の数量を表す
\(storage\) は各製品の各月における保管量を表す

出力形式:

{
    "sell": [[sell_{k, i} for k in 1,...,K] for i in 1,...I],
    "manufacture": [[manufacture_{k, i} for k in 1,...,K] for i in 1,...I],
    "storage": [[storage_{k, i} for k in 1,...,K] for i in 1,...,I]
}

定式化

この問題を線形計画問題（LP問題）として定式化するために、以下の変数、パラメータ、制約、目的関数を定義します。

変数

\(\text{sell}_{k,i}\): 月 \(i\) に製品 \(k\) を販売する数量
\(\text{manufacture}_{k,i}\): 月 \(i\) に製品 \(k\) を製造する数量
\(\text{storage}_{k,i}\): 月 \(i\) に製品 \(k\) の在庫数

パラメータ

\(\text{num\_machines}[m]\): 機械 \(m\) の総数
\(\text{profit}[k]\): 製品 \(k\) の単位利益
\(\text{time}[k][m]\): 製品 \(k\) の製造に機械 \(m\) が必要な時間 (hours/unit)
\(\text{maintain}[i][m]\): 月 \(i\) にメンテナンス中の機械 \(m\) の数
\(\text{limit}[k][i]\): 月 \(i\) の製品 \(k\) の販売上限
\(\text{store\_price}\): 単位在庫コスト
\(\text{keep\_quantity}\): 月末の各製品の保持在庫数
\(\text{n\_workhours}\): 1日の作業時間
\(\text{n\_days\_per\_month} = 24\): 1ヶ月の日数

目的関数

最大化すべき目的関数は、総利益から総在庫コストを引いたものです。

\[ \text{Maximize} \quad \sum_{i=1}^{I} \sum_{k=1}^{K} \left( \text{profit}[k] \cdot \text{sell}_{k,i} - \text{store\_price} \cdot \text{storage}_{k,i} \right) \]

制約

製品の在庫量のバランス: 各月の在庫は、前月の在庫に当月の製造量を加え、販売量を引いたものと等しくなければなりません。 \[ \text{storage}_{k,i} = \text{storage}_{k,i-1} + \text{manufacture}_{k,i} - \text{sell}_{k,i} \quad \forall k, \forall i \] （ただし、\(\text{storage}_{k,0} = 0\) とする）
機械の利用可能時間: 各機械の利用可能時間は、作業時間とメンテナンス時間に基づいて制限されます。 \[ \sum_{k=1}^{K} \text{time}[k][m] \cdot \text{manufacture}_{k,i} \leq (\text{num\_machines}[m] - \text{maintain}[i][m]) \cdot \text{n\_workhours} \cdot \text{n\_days\_per\_month} \quad \forall m, \forall i \]
販売上限: 各製品の販売数量は、その月の販売上限を超えてはなりません。 \[ \text{sell}_{k,i} \leq \text{limit}[k][i] \quad \forall k, \forall i \]
在庫上限: 各製品の在庫は100を超えてはなりません。 \[ \text{storage}_{k,i} \leq 100 \quad \forall k, \forall i \]
期末在庫: 各製品の最終月の在庫は \(\text{keep\_quantity}\) である必要があります。 \[ \text{storage}_{k,I} = \text{keep\_quantity} \quad \forall k \]

OptiMUSのサンプルデータにバグがあった．正しい入力データを準備することも最低限重要である．

import gurobipy as gp
from gurobipy import GRB

def optimize_factory_production(data):
    num_machines = data['num_machines']
    profit = data['profit']
    time = data['time']
    maintain = data['maintain']
    limit = data['limit']
    store_price = data['store_price']
    keep_quantity = data['keep_quantity']
    n_workhours = data['n_workhours']
    
    M = len(num_machines)  # Number of machine types
    K = len(profit)  # Number of product types
    I = len(limit[0])  # Number of months

    days_per_month = 24  # Assuming each month has 24 working days
    total_workhours = n_workhours * 6 * days_per_month

    # Create a new model
    model = gp.Model("factory_production")
    
    # Decision variables
    sell = model.addVars(K, I, name="sell", vtype=GRB.CONTINUOUS)
    manufacture = model.addVars(K, I, name="manufacture", vtype=GRB.CONTINUOUS)
    storage = model.addVars(K, I, name="storage", vtype=GRB.CONTINUOUS)
    
    # Objective function: Maximize profit
    revenue = gp.quicksum(sell[k, i] * profit[k] for k in range(K) for i in range(I))
    storage_cost = gp.quicksum(storage[k, i] * store_price for k in range(K) for i in range(I))
    model.setObjective(revenue - storage_cost, GRB.MAXIMIZE)
    
    # Constraints
    
    # Initial storage (no initial stock)
    for k in range(K):
        model.addConstr(storage[k, 0] == 0)
    
    # Storage capacity constraints
    for i in range(I):
        for k in range(K):
            model.addConstr(storage[k, i] <= 100)
    
    # End-of-month storage requirements
    for k in range(K):
        model.addConstr(storage[k, I-1] == keep_quantity)
    
    # Market limitations
    for i in range(I):
        for k in range(K):
            model.addConstr(sell[k, i] <= limit[k][i])
    
    # Manufacturing time constraints
    for i in range(I):
        for m in range(M):
            #データの行列の入れ間違え； maintain[i][m] を変更して対処
            available_hours = total_workhours * (num_machines[m] - maintain[ｍ][i])
            model.addConstr(gp.quicksum(manufacture[k, i] * time[k][m] for k in range(K)) <= available_hours)
    
    # Storage dynamics
    for k in range(K):
        for i in range(1, I):
            model.addConstr(storage[k, i] == storage[k, i-1] + manufacture[k, i] - sell[k, i])
    
    # Ensure manufactured and sold quantities are non-negative
    for k in range(K):
        for i in range(I):
            model.addConstr(manufacture[k, i] >= 0)
            model.addConstr(sell[k, i] >= 0)
    
    # Optimize model
    model.optimize()
    
    # Prepare the result
    result = {
        "sell": [[sell[k, i].x for i in range(I)] for k in range(K)],
        "manufacture": [[manufacture[k, i].x for i in range(I)] for k in range(K)],
        "storage": [[storage[k, i].x for i in range(I)] for k in range(K)]
    }
    
    return result

# Example input
data = {
    "num_machines": [4, 2, 3, 1, 1],
    "profit": [10, 6, 8, 4, 11, 9, 3],
    "time": [
        [0.5, 0.1, 0.2, 0.05, 0.0],
        [0.7, 0.2, 0.0, 0.03, 0.0],
        [0.0, 0.0, 0.8, 0.0, 0.01],
        [0.0, 0.3, 0.0, 0.07, 0.0],
        [0.3, 0.0, 0.0, 0.1, 0.05],
        [0.2, 0.6, 0.0, 0.0, 0.0],  #データのミス（この行がない）
        [0.5, 0.0, 0.6, 0.08, 0.05]
    ],
    "maintain": [ #行と列が反対であるが， プログラムを変更して対処
        [1, 0, 0, 0, 1, 0],
        [0, 0, 0, 1, 1, 0],
        [0, 2, 0, 0, 0, 1],
        [0, 0, 1, 0, 0, 0],
        [0, 0, 0, 0, 0, 1]
    ],
    "limit": [
        [500,600,300,200,0,500],
        [1000,500,600,300,100,500],
        [300,200,0,400,500,100],
        [300,0,0,500,100,300],
        [800,400,500,200,1000,1100],
        [200,300,400,0,300,500],
        [100,150,100,100,0,60]
    ],
    "store_price": 0.5,
    "keep_quantity": 100,
    "n_workhours": 8.0
}

# Call the optimization function
result = optimize_factory_production(data)
print(result)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 247 rows, 126 columns and 442 nonzeros
Model fingerprint: 0x3336220a
Coefficient statistics:
  Matrix range     [1e-02, 1e+00]
  Objective range  [5e-01, 1e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [6e+01, 5e+03]
Presolve removed 247 rows and 126 columns
Presolve time: 0.01s
Presolve: All rows and columns removed
Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    9.3450000e+04   0.000000e+00   1.000000e+00      0s
Extra simplex iterations after uncrush: 2
       2    9.3450000e+04   0.000000e+00   0.000000e+00      0s

Solved in 2 iterations and 0.01 seconds (0.00 work units)
Optimal objective  9.345000000e+04
{'sell': [[500.0, 600.0, 100.0, 200.0, 0.0, 500.0], [1000.0, 500.0, 100.0, 300.0, 100.0, 500.0], [300.0, 200.0, 0.0, 400.0, 500.0, 0.0], [300.0, 0.0, 0.0, 500.0, 100.0, 300.0], [800.0, 400.0, 100.0, 200.0, 1000.0, 0.0], [200.0, 300.0, 400.0, 0.0, 300.0, 500.0], [100.0, 150.0, 100.0, 100.0, 0.0, 0.0]], 'manufacture': [[0.0, 700.0, 0.0, 200.0, 0.0, 600.0], [0.0, 600.0, 0.0, 300.0, 100.0, 600.0], [0.0, 200.0, 0.0, 400.0, 600.0, 0.0], [0.0, 0.0, 0.0, 500.0, 100.0, 400.0], [0.0, 500.0, 0.0, 200.0, 1100.0, 0.0], [0.0, 300.0, 400.0, 0.0, 300.0, 600.0], [0.0, 250.0, 0.0, 100.0, 100.0, 0.0]], 'storage': [[0.0, 100.0, 0.0, 0.0, 0.0, 100.0], [0.0, 100.0, 0.0, 0.0, 0.0, 100.0], [0.0, 0.0, 0.0, 0.0, 100.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 100.0, 0.0, 0.0, 100.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 100.0, 0.0, 0.0, 100.0, 100.0]]}

工場計画問題2

model_building_in_mathematical_programming (problem 4)

上の問題にメンテナンスの月を選択するように変更したモデル．

修正ポイント

メンテナンスデータの修正: 各機械タイプごとに必要なメンテナンス月数 down を考慮。
決定変数の修正: maintain 変数を追加し、各月にメンテナンス中の機械数を表すように修正。
制約の追加: メンテナンス制約を追加し、各機械タイプごとに必要なメンテナンス月数を確保する。
総作業時間の修正: 一月あたりの総作業時間を二交代制として修正。

このコードを使用して、工場の製造、販売、保管、およびメンテナンスの方針を最適化することができます。

import gurobipy as gp
from gurobipy import GRB

def optimize_factory_production(data):
    num_machines = data['num_machines']
    profit = data['profit']
    time = data['time']
    down = data['down']
    limit = data['limit']
    store_price = data['store_price']
    keep_quantity = data['keep_quantity']
    n_workhours = data['n_workhours']
    
    M = len(num_machines)  # Number of machine types
    K = len(profit)  # Number of product types
    I = len(limit[0])  # Number of months

    days_per_month = 24  # Assuming each month has 24 working days
    total_workhours = n_workhours * 2 * 6 * days_per_month

    # Create a new model
    model = gp.Model("factory_production")
    
    # Decision variables
    sell = model.addVars(K, I, name="sell", vtype=GRB.CONTINUOUS)
    manufacture = model.addVars(K, I, name="manufacture", vtype=GRB.CONTINUOUS)
    storage = model.addVars(K, I, name="storage", vtype=GRB.CONTINUOUS)
    maintain = model.addVars(M, I, name="maintain", vtype=GRB.INTEGER)

    # Objective function: Maximize profit
    revenue = gp.quicksum(sell[k, i] * profit[k] for k in range(K) for i in range(I))
    storage_cost = gp.quicksum(storage[k, i] * store_price for k in range(K) for i in range(I))
    model.setObjective(revenue - storage_cost, GRB.MAXIMIZE)
    
    # Constraints
    
    # Initial storage (no initial stock)
    for k in range(K):
        model.addConstr(storage[k, 0] == 0)
    
    # Storage capacity constraints
    for i in range(I):
        for k in range(K):
            model.addConstr(storage[k, i] <= 100)
    
    # End-of-month storage requirements
    for k in range(K):
        model.addConstr(storage[k, I-1] == keep_quantity)
    
    # Market limitations
    for i in range(I):
        for k in range(K):
            model.addConstr(sell[k, i] <= limit[k][i])
    
    # Manufacturing time constraints considering maintenance
    for i in range(I):
        for m in range(M):
            available_hours = total_workhours * (num_machines[m] - maintain[m, i])
            model.addConstr(gp.quicksum(manufacture[k, i] * time[k][m] for k in range(K)) <= available_hours)
    
    # Maintenance constraints
    for m in range(M):
        model.addConstr(gp.quicksum(maintain[m, i] for i in range(I)) == down[m])

    # Storage dynamics
    for k in range(K):
        for i in range(1, I):
            model.addConstr(storage[k, i] == storage[k, i-1] + manufacture[k, i] - sell[k, i])
    
    # Ensure manufactured and sold quantities are non-negative
    for k in range(K):
        for i in range(I):
            model.addConstr(manufacture[k, i] >= 0)
            model.addConstr(sell[k, i] >= 0)
    
    # Optimize model
    model.optimize()
    
    # Prepare the result
    result = {
        "sell": [[sell[k, i].x for i in range(I)] for k in range(K)],
        "manufacture": [[manufacture[k, i].x for i in range(I)] for k in range(K)],
        "storage": [[storage[k, i].x for i in range(I)] for k in range(K)],
        "maintain": [[maintain[m, i].x for i in range(I)] for m in range(M)]
    }
    
    return result

data = {
    "num_machines": [4, 2, 3, 1, 1],
    "profit": [10, 6, 8, 4, 11, 9, 3],
    "time": [
        [0.5, 0.1, 0.2, 0.05, 0.0],
        [0.7, 0.2, 0.0, 0.03, 0.0],
        [0.0, 0.0, 0.8, 0.0, 0.01],
        [0.0, 0.3, 0.0, 0.07, 0.0],
        [0.3, 0.0, 0.0, 0.1, 0.05],
        [0.2, 0.6, 0.0, 0.0, 0.0],  #データのミス（この行がない）
        [0.5, 0.0, 0.6, 0.08, 0.05]
    ],
    "maintain": [
        [1, 0, 0, 0, 1, 0],
        [0, 0, 0, 1, 1, 0],
        [0, 2, 0, 0, 0, 1],
        [0, 0, 1, 0, 0, 0],
        [0, 0, 0, 0, 0, 1]
    ],
    "limit": [
        [500,600,300,200,0,500],
        [1000,500,600,300,100,500],
        [300,200,0,400,500,100],
        [300,0,0,500,100,300],
        [800,400,500,200,1000,1100],
        [200,300,400,0,300,500],
        [100,150,100,100,0,60]
    ],
    "down": [2, 2, 3, 1, 1],
    "store_price": 0.5,
    "keep_quantity": 100,
    "n_workhours": 8.0
}

# Call the optimization function
result = optimize_factory_production(data)
print(result)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 252 rows, 156 columns and 502 nonzeros
Model fingerprint: 0x205416d6
Variable types: 126 continuous, 30 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e-02, 2e+03]
  Objective range  [5e-01, 1e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 9e+03]
Found heuristic solution: objective 25150.000000
Presolve removed 252 rows and 156 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 16 available processors)

Solution count 2: 116280 25150 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.162800000000e+05, best bound 1.162800000000e+05, gap 0.0000%
{'sell': [[500.0, 600.0, 300.0, 200.0, 0.0, 500.0], [1000.0, 500.0, 600.0, 300.0, 100.0, 500.0], [300.0, 200.0, 0.0, 400.0, 500.0, 100.0], [300.0, 0.0, 0.0, 500.0, 100.0, 300.0], [800.0, 400.0, 500.0, 200.0, 1000.0, 1100.0], [200.0, 300.0, 400.0, 0.0, 300.0, 500.0], [100.0, 150.0, 100.0, 100.0, 0.0, 60.0]], 'manufacture': [[0.0, 600.0, 300.0, 200.0, 0.0, 600.0], [0.0, 500.0, 600.0, 300.0, 100.0, 600.0], [0.0, 200.0, 0.0, 400.0, 500.0, 200.0], [0.0, 0.0, 0.0, 500.0, 100.0, 400.0], [0.0, 400.0, 500.0, 200.0, 1000.0, 1200.0], [0.0, 300.0, 400.0, 0.0, 300.0, 600.0], [0.0, 150.0, 100.0, 100.0, 0.0, 160.0]], 'storage': [[0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0], [0.0, 0.0, 0.0, 0.0, 0.0, 100.0]], 'maintain': [[2.0, -0.0, -0.0, -0.0, -0.0, -0.0], [2.0, -0.0, -0.0, -0.0, -0.0, -0.0], [3.0, -0.0, -0.0, -0.0, -0.0, -0.0], [1.0, -0.0, -0.0, -0.0, -0.0, -0.0], [1.0, -0.0, -0.0, -0.0, -0.0, -0.0]]}

発電最適化

introduction_to_linear_optimization (problem 3)

問題の定式化

変数

\(x_{\text{coal}, t}\): 年 \(t\) に追加する石炭火力発電容量 (MW)
\(x_{\text{nuke}, t}\): 年 \(t\) に追加する原子力発電容量 (MW)

パラメータ

\(\text{demand}_t\): 年 \(t\) の需要 (MW)
\(\text{oil\_cap}_t\): 年 \(t\) に利用可能な石油火力発電容量 (MW)
\(\text{coal\_cost}\): 石炭火力発電の1 MWあたりの資本コスト
\(\text{nuke\_cost}\): 原子力発電の1 MWあたりの資本コスト
\(\text{max\_nuke}\): 原子力発電容量の総容量に対する最大割合 (%)
\(\text{coal\_life}\): 石炭火力発電所の寿命 (年)
\(\text{nuke\_life}\): 原子力発電所の寿命 (年)

目的関数

総コストを最小化する： \[ \text{minimize} \quad \sum_{t=1}^{T} (\text{coal\_cost} \cdot x_{\text{coal}, t} + \text{nuke\_cost} \cdot x_{\text{nuke}, t}) \]

制約条件

供給が需要を満たす制約： \[ \sum_{j=\max(1, t-\text{coal\_life}+1)}^t x_{\text{coal}, j} + \sum_{j=\max(1, t-\text{nuke\_life}+1)}^t x_{\text{nuke}, j} + \text{oil\_cap}_t \geq \text{demand}_t \quad \forall t \]
原子力発電容量の上限制約： \[ \sum_{j=\max(1, t-\text{nuke\_life}+1)}^t x_{\text{nuke}, j} \leq \frac{\text{max\_nuke}}{100} \left( \sum_{j=\max(1, t-\text{coal\_life}+1)}^t x_{\text{coal}, j} + \sum_{j=\max(1, t-\text{nuke\_life}+1)}^t x_{\text{nuke}, j} + \text{oil\_cap}_t \right) \quad \forall t \]
非負制約： \[ x_{\text{coal}, t} \geq 0 \quad \forall t \] \[ x_{\text{nuke}, t} \geq 0 \quad \forall t \]

出力フォーマット

{
    "coal_cap_added": [coal_t for t in 1, ..., T],
    "nuke_cap_added": [nuke_t for t in 1, ..., T],
    "total_cost": total_cost
}

from gurobipy import Model, GRB, quicksum

def electricity_capacity_planning(demand, oil_cap, coal_cost, nuke_cost, max_nuke, coal_life, nuke_life):
    T = len(demand)
    
    # モデルの作成
    model = Model("electricity_capacity_planning")
    
    # 変数の作成
    coal_cap_added = model.addVars(T, name="coal_cap_added", lb=0)
    nuke_cap_added = model.addVars(T, name="nuke_cap_added", lb=0)
    
    # 目的関数の設定
    model.setObjective(quicksum(coal_cost * coal_cap_added[t] + nuke_cost * nuke_cap_added[t] for t in range(T)), GRB.MINIMIZE)
    
    # 供給が需要を満たす制約
    for t in range(T):
        model.addConstr(quicksum(coal_cap_added[j] for j in range(max(0, t - coal_life + 1), t + 1)) +
                        quicksum(nuke_cap_added[j] for j in range(max(0, t - nuke_life + 1), t + 1)) +
                        oil_cap[t] >= demand[t], name=f"demand_constraint_{t}")
    
    # 原子力発電容量の上限制約
    for t in range(T):
        nuke_capacity = quicksum(nuke_cap_added[j] for j in range(max(0, t - nuke_life + 1), t + 1))
        total_capacity = quicksum(coal_cap_added[j] for j in range(max(0, t - coal_life + 1), t + 1)) + \
                         quicksum(nuke_cap_added[j] for j in range(max(0, t - nuke_life + 1), t + 1)) + \
                         oil_cap[t]
        model.addConstr(nuke_capacity <= (max_nuke / 100) * total_capacity, name=f"nuke_constraint_{t}")
    
    # 最適化の実行
    model.optimize()
    
    if model.status == GRB.OPTIMAL:
        coal_cap_added_solution = [coal_cap_added[t].x for t in range(T)]
        nuke_cap_added_solution = [nuke_cap_added[t].x for t in range(T)]
        total_cost = model.objVal
        return {
            "coal_cap_added": coal_cap_added_solution,
            "nuke_cap_added": nuke_cap_added_solution,
            "total_cost": total_cost
        }
    else:
        return None

# 入力データの例
input_data = {
  "demand": [10, 15, 17, 20, 13, 19, 10, 25, 30, 33, 30, 35],
  "oil_cap": [15, 14, 13, 10, 10, 7, 5, 5, 5, 5, 5, 5],
  "coal_cost": 10,
  "nuke_cost": 5,
  "max_nuke": 20,
  "coal_life": 5,
  "nuke_life": 10
}

# 結果の計算
output_data = electricity_capacity_planning(**input_data)
print(output_data)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 24 rows, 24 columns and 250 nonzeros
Model fingerprint: 0xbe87bf6e
Coefficient statistics:
  Matrix range     [2e-01, 1e+00]
  Objective range  [5e+00, 1e+01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [1e+00, 3e+01]
Presolve removed 1 rows and 0 columns
Presolve time: 0.00s
Presolved: 23 rows, 24 columns, 248 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    0.0000000e+00   1.630000e+02   0.000000e+00      0s
      12    3.5300000e+02   0.000000e+00   0.000000e+00      0s

Solved in 12 iterations and 0.01 seconds (0.00 work units)
Optimal objective  3.530000000e+02
{'coal_cap_added': [0.0, 1.0, 7.800000000000001, 0.0, 0.0, 0.0, 0.0, 15.0, 4.0, 2.3999999999999995, 0.0, 1.6000000000000005], 'nuke_cap_added': [0.0, 0.0, 0.0, 1.1999999999999997, 0.0, 2.0, 0.0, 1.8000000000000003, 1.0, 0.6000000000000005, 0.0, 0.39999999999999947], 'total_cost': 353.0}

人員計画

model_building_in_mathematical_programming (problem 5)

A company is undergoing a number of changes that will affect its manpower requirements in future years.
Owing to the installation of new machinery, fewer unskilled but more skilled and semi-skilled workers will be required.
In addition to this, a downturn in trade is expected in the next year, which will reduce the need for workers in all categories.
The estimated manpower requirements for manpower in year is .
The current number of manpower is .
The company wishes to decide its policy with regard to recruitment, retraining, redundancy and short-time working in the next years.
There is a natural wastage of labour.
A fairly large number of workers leave during their first year. After this, the rate of leaving is much smaller.
Taking this into account, the wastage rates of manpower with less than one year’s service is .
The wastage rates of manpower with more than one year’s service is .
There has been no recent recruitment and all the workers in the current labour force have been employed for more than one year.
It is possible to recruit a limited number of workers from outside.
In any one year, the maximum number that can be recruited in manpower is .
The redundancy payment to manpower is .
It is possible to employ up to more workers over the whole company than are needed, but the extra costs per employee per year is for manpower .
Up to workers in each category of skill can be put on short-time working.
The cost of this (per employee per year) is for manpower .
An employee on short-time working meets the production requirements of half a full-time employee.

INPUT FORMAT:

{
    "requirement": [[requirement_{k, i} for i in 1,...,I] for k in 1,...,K],
    "strength": [strength_{k} for k in 1,...,K],
    "lessonewaste": [lessonewaste_{k} for k in 1,...,K],
    "moreonewaste": [lessonewaste_{k} for k in 1,...,K],
    "recruit": [recruit_{k} for k in 1,...,K],
    "costredundancy": [costredundancy_{k} for k in 1,...,K],
    "num_overman": num_overman,
    "costoverman": [costoverman_{k} for k in 1,...,K],
    "num_shortwork": num_shortwork
    "costshort": [costshort_{k} for k in 1,...,K],
}

OBJECTIVE: The company’s declared objective is to minimize redundancy. How should they operate in order to do this?

OUTPUT INFO:

represents the number of each recruit manpower to employ each year
represents the number of each overmanning manpower to employ each year
represents the number of each short-time manpower to employ each year

OUTPUT FORMAT:

{
    "recruit": [[recruit_{k, i} for i in 1,...,I] for k in 1,...,K],
    "overmanning": [[overmanning_{k, i} for i in 1,...,I] for k in 1,...,K],
    "short": [[short_{k, i} for i in 1,...,I] for k in 1,...,K]
}

この問題を線形計画問題（LP問題）として定式化するために、以下の変数、パラメータ、制約、目的関数を定義します。

定式化

変数

\(\text{recruit}_{k,i}\): 年 \(i\) に職種 \(k\) の新規雇用者数
\(\text{overmanning}_{k,i}\): 年 \(i\) に職種 \(k\) の過剰雇用者数
\(\text{short}_{k,i}\): 年 \(i\) に職種 \(k\) の短時間勤務者数
\(\text{redundancy}_{k,i}\): 年 \(i\) に職種 \(k\) の余剰人員数

パラメータ

\(\text{requirement}[k][i]\): 年 \(i\) に必要な職種 \(k\) の人員数
\(\text{strength}[k]\): 現在の職種 \(k\) の人員数
\(\text{lessonewaste}[k]\): 職種 \(k\) の1年未満の自然離職率
\(\text{moreonewaste}[k]\): 職種 \(k\) の1年以上の自然離職率
\(\text{recruit\_limit}[k]\): 年 \(i\) に雇用できる職種 \(k\) の最大人数
\(\text{costredundancy}[k]\): 職種 \(k\) の余剰人員1人あたりのコスト
\(\text{num\_overman}\): 会社全体で許可される過剰雇用者数の上限
\(\text{costoverman}[k]\): 職種 \(k\) の過剰雇用者1人あたりのコスト
\(\text{num\_shortwork}\): 各職種で短時間勤務が許可される人数の上限
\(\text{costshort}[k]\): 職種 \(k\) の短時間勤務者1人あたりのコスト

目的関数

余剰人員コストを最小化するための目的関数:

\[ \text{Minimize} \quad \sum_{i=1}^{I} \sum_{k=1}^{K} \left( \text{costredundancy}[k] \cdot \text{redundancy}_{k,i} \right) \]

制約

人員バランス: 各年の人員は前年の人員、採用者、短時間勤務者、過剰雇用者、余剰人員、および自然離職を考慮してバランスを取る必要があります。

\[ \text{strength}_{k, i} = \text{strength}_{k, i-1} \cdot (1 - \text{moreonewaste}[k]) + \text{recruit}_{k, i} \cdot (1 - \text{lessonewaste}[k]) - \text{redundancy}_{k, i} \]
労働力要求: 各年の労働力要求は、短時間勤務者を考慮してバランスを取る必要があります。

\[ \text{strength}_{k, i} + \text{overmanning}_{k, i} + 0.5 \cdot \text{short}_{k, i} \geq \text{requirement}_{k, i} \]
新規採用の上限: 各年における職種 \(k\) の新規採用者数は、その職種の採用上限を超えてはなりません。

\[ \text{recruit}_{k, i} \leq \text{recruit\_limit}[k] \]
過剰雇用の上限: 会社全体での過剰雇用者数は、指定された上限を超えてはなりません。

\[ \sum_{k=1}^{K} \text{overmanning}_{k, i} \leq \text{num\_overman} \]
短時間勤務の上限: 各職種における短時間勤務者数は、その職種の上限を超えてはなりません。

\[ \text{short}_{k, i} \leq \text{num\_shortwork} \]

入力フォーマット

{
    "requirement": [[requirement_{k, i} for i in 1,...,I] for k in 1,...,K],
    "strength": [strength_{k} for k in 1,...,K],
    "lessonewaste": [lessonewaste_{k} for k in 1,...,K],
    "moreonewaste": [moreonewaste_{k} for k in 1,...,K],
    "recruit": [recruit_{k} for k in 1,...,K],
    "costredundancy": [costredundancy_{k} for k in 1,...,K],
    "num_overman": num_overman,
    "costoverman": [costoverman_{k} for k in 1,...,K],
    "num_shortwork": num_shortwork,
    "costshort": [costshort_{k} for k in 1,...,K]
}

出力フォーマット

{
    "recruit": [[recruit_{k, i} for i in 1,...,I] for k in 1,...,K],
    "overmanning": [[overmanning_{k, i} for i in 1,...,I] for k in 1,...,K],
    "short": [[short_{k, i} for i in 1,...,I] for k in 1,...,K]
}

問題文が不正確で，ダウングレード（上位のスキルをもった労働者を下位のスキルの労働者として使用すること）の変数が考慮されていない．また，問題文で述べられているretrainingに対する変数が抜けている．これは，スキルを上位のものにするために必要である．したがって生成された定式化もコードも不完全なものとなり，結果は実行不能になる．

from gurobipy import Model, GRB, quicksum

# 入力データ
data = {
  "requirement": [
    [1000, 1400, 1000],
    [500, 2000, 1500],
    [0, 2500, 2000]
  ],
  "strength": [2000, 1500, 1000],
  "lessonewaste": [0.25, 0.2, 0.1],
  "moreonewaste": [0.1, 0.05, 0.05],
  "recruit": [500, 800, 500],
  "costredundancy": [200, 500, 500],
  "num_overman": 150,
  "costoverman": [1500, 2000, 3000],
  "num_shortwork": 50,
  "costshort": [500, 400, 400]
}

requirement = data["requirement"]
strength = data["strength"]
lessonewaste = data["lessonewaste"]
moreonewaste = data["moreonewaste"]
recruit_limit = data["recruit"]
costredundancy = data["costredundancy"]
num_overman = data["num_overman"]
costoverman = data["costoverman"]
num_shortwork = data["num_shortwork"]
costshort = data["costshort"]

# 定数
K = len(strength)
I = len(requirement[0])

# モデルの作成
model = Model('ManpowerPlanning')

# 変数の作成
recruit = model.addVars(K, I, name="recruit", vtype=GRB.CONTINUOUS, lb=0)
overmanning = model.addVars(K, I, name="overmanning", vtype=GRB.CONTINUOUS, lb=0)
short = model.addVars(K, I, name="short", vtype=GRB.CONTINUOUS, lb=0)
redundancy = model.addVars(K, I, name="redundancy", vtype=GRB.CONTINUOUS, lb=0)
strength_vars = model.addVars(K, I+1, name="strength", vtype=GRB.CONTINUOUS, lb=0)

# 初期人員数
for k in range(K):
    model.addConstr(strength_vars[k, 0] == strength[k])

# 人員バランスの制約
for k in range(K):
    for i in range(1, I+1):
        if i == 1:
            model.addConstr(
                strength_vars[k, i] == strength_vars[k, i-1] * (1 - moreonewaste[k]) + recruit[k, i-1] * (1 - lessonewaste[k]) - redundancy[k, i-1]
            )
        else:
            model.addConstr(
                strength_vars[k, i] == strength_vars[k, i-1] * (1 - moreonewaste[k]) + recruit[k, i-1] * (1 - lessonewaste[k]) - redundancy[k, i-1]
            )

# 労働力要求の制約
for k in range(K):
    for i in range(I):
        model.addConstr(
            strength_vars[k, i+1] - overmanning[k, i] + 0.5 * short[k, i] == requirement[k][i]
        )

# 新規採用の上限
for k in range(K):
    for i in range(I):
        model.addConstr(
            recruit[k, i] <= recruit_limit[k]
        )

# 過剰雇用の上限
for i in range(I):
    model.addConstr(
        quicksum(overmanning[k, i] for k in range(K)) <= num_overman
    )

# 短時間勤務の上限
for k in range(K):
    for i in range(I):
        model.addConstr(
            short[k, i] <= num_shortwork
        )

# 目的関数の設定
model.setObjective(
    quicksum(costredundancy[k] * redundancy[k, i] for k in range(K) for i in range(I)), GRB.MINIMIZE
)

# 最適化の実行
model.optimize()

# 結果の表示
if model.status == GRB.OPTIMAL:
    recruit_sol = model.getAttr('x', recruit)
    overmanning_sol = model.getAttr('x', overmanning)
    short_sol = model.getAttr('x', short)

    output = {
        "recruit": [[recruit_sol[k, i] for i in range(I)] for k in range(K)],
        "overmanning": [[overmanning_sol[k, i] for i in range(I)] for k in range(K)],
        "short": [[short_sol[k, i] for i in range(I)] for k in range(K)],
    }

    print(output)
else:
    print("Optimal solution was not found.")

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 42 rows, 48 columns and 93 nonzeros
Model fingerprint: 0x9064da2e
Coefficient statistics:
  Matrix range     [5e-01, 1e+00]
  Objective range  [2e+02, 5e+02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e+01, 2e+03]
Presolve removed 21 rows and 3 columns
Presolve time: 0.00s

Solved in 0 iterations and 0.00 seconds (0.00 work units)
Infeasible model
Optimal solution was not found.

産業連関表

model_building_in_mathematical_programming (problem 11)

An economy consists of several industries. （産業連関表；本の9節）

Each unit produced by one of the industries (a unit will be taken as £1’s worth of value of production) requires inputs from possibly its own industry as well as other industries.
Industry requires input from industy .
Industry requires units of manpower.
There is a time lag in the economy so that the output in year t + 1 requires an input in year t.
Output from an industry may also be used to build productive capacity for itself or other industries in future years.
To give unit increases (capacity for £1’s worth of extra production) in productive capacity of industry , units of industry and units of manpower are needed.
Input from an industry in year t results in a (permanent) increase in productive capacity in year t + 2.
Stocks of goods may be held from year to year.
At present (year 0), the stocks of industry is , and the capacity of industry is .
An exogenous consumption requirement of units from industry has to be satisfied in every year except year 0.

INPUT FORMAT:

{
    "inputone": [[inputone_{k, j} for j in 1,...,K] for k in 1,...,K],
    "manpowerone": [manpowerone_{k} for k in 1,...,K],
    "inputtwo": [[inputtwo_{k, j} for j in 1,...,K] for k in 1,...,K],
    "manpowertwo": [manpowertwo_{k} for k in 1,...,K],
    "stock": [stock_{k} for k in 1,...,K],
    "capacity": [capacity_{k} for k in 1,...,K],
    "demand": [demand_{k} for k in 1,...,K]
}

OBJECTIVE: Maximize the total manpower requirement over five years

OUTPUT INFO:

represents the amount of units produced by each industry in each year
represents the amount of units used to build productive capacity for each industry in each year
represents the amount of stock of each industry held in each year

###OUTPUT FORMAT:

{
    "produce": [[produce_{k, t} for t in 1,...,T] for k in 1,...,K],
    "buildcapa": [[buildcapa_{k, t} for t in 1,...,T] for k in 1,...,K],
    "stockhold": [[stock_{k, t} for t in 1,...,T] for for k in 1,...,K]
}

{
    "inputone": [[0.1, 0.1, 0.2],
                 [0.5, 0.1, 0.1],
                 [0.5, 0.2, 0.2]],
    "manpowerone": [0.6, 0.3, 0.2],
    "inputtwo": [[0.0, 0.1, 0.2],
                 [0.7, 0.1, 0.1],
                 [0.9, 0.2, 0.2]],
    "manpowertwo": [0.4, 0.2, 0.1],
    "stock": [150, 80, 100],
    "capacity": [300, 350, 280],
    "demand": [6e+07, 6e+07, 3e+07]
}

この問題を解決するために、次のような数理最適化モデルを定式化します。以下のコードは、Gurobiを使用して実装します。

モデルの定式化

変数
- produce[k, t]: 年度\(t\) に産業\(k\) が生産する単位量。
- buildcapa[k, t]: 年度\(t\) に産業\(k\) が将来の生産能力のために使用する単位量。
- stockhold[k, t]: 年度\(t\) に産業\(k\) が保有する在庫量。
制約
- 生産量と在庫の関係：年度\(t\) の生産量と在庫は、需要と生産能力を満たさなければならない。
- 生産能力の更新：産業\(k\) の生産能力は、前年の生産能力と新たに追加された生産能力の和。
- インプットのバランス：各産業の生産には、他の産業からのインプットが必要。
- 労働力の使用：各産業の生産には一定の労働力が必要。
目的関数
- 総労働力の最大化。

説明

データの読み込み:
- 提供されたデータを読み込みます。
変数の作成:
- produce: 各産業が各年に生産する量。
- buildcapa: 各産業が将来の生産能力を構築するために使用する量。
- stockhold: 各産業が各年に保有する在庫量。
制約の設定:
- 初期在庫を設定し、各年の生産と在庫のバランスを制約します。
- 各産業の生産能力の制約を設定します。
- 各産業のインプットバランスの制約を設定します。
- 各産業の労働力使用の制約を設定します。
目的関数の設定:
- 総労働力の最大化を目的とします。
最適化の実行:
- モデルを最適化し、結果を表示します。

このコードを実行することで、各産業が各年にどれだけ生産し、どれだけの能力を構築し、どれだけの在庫を保有するべきかが求められます。

from gurobipy import Model, GRB, quicksum

# 入力データ
data = {
    "inputone": [[0.1, 0.1, 0.2],
                 [0.5, 0.1, 0.1],
                 [0.5, 0.2, 0.2]],
    "manpowerone": [0.6, 0.3, 0.2],
    "inputtwo": [[0.0, 0.1, 0.2],
                 [0.7, 0.1, 0.1],
                 [0.9, 0.2, 0.2]],
    "manpowertwo": [0.4, 0.2, 0.1],
    "stock": [150, 80, 100],
    "capacity": [300, 350, 280],
    "demand": [60, 60, 30]
}

inputone = data["inputone"]
manpowerone = data["manpowerone"]
inputtwo = data["inputtwo"]
manpowertwo = data["manpowertwo"]
stock = data["stock"]
capacity = data["capacity"]
demand = data["demand"]

K = len(inputone)  # 産業の数
T = 5  # 年数

# モデルの作成
model = Model("EconomyOptimization")

# 変数の作成
produce = model.addVars(K, T, name="produce", lb=0, vtype=GRB.CONTINUOUS)
buildcapa = model.addVars(K, T, name="buildcapa", lb=0, vtype=GRB.CONTINUOUS)
stockhold = model.addVars(K, T+1, name="stockhold", lb=0, vtype=GRB.CONTINUOUS)

# 初期在庫を設定
for k in range(K):
    model.addConstr(stockhold[k, 0] == stock[k])

# 制約条件の設定
for t in range(T):
    for k in range(K):
        # 在庫と生産量のバランス
        if t == 0:
            model.addConstr(produce[k, t] + stockhold[k, t] == demand[k] + stockhold[k, t+1])
        else:
            model.addConstr(produce[k, t] + stockhold[k, t] == demand[k] + stockhold[k, t+1] + produce[k, t-1])

        # 生産能力の制約
        model.addConstr(produce[k, t] <= capacity[k])
        
        # インプットのバランス
        model.addConstr(quicksum(inputone[k][j] * produce[j, t] for j in range(K)) + quicksum(inputtwo[k][j] * buildcapa[j, t] for j in range(K)) 
                        <= stockhold[k, t] + produce[k, t])

        # 労働力の使用
        model.addConstr(quicksum(manpowerone[j] * produce[j, t] for j in range(K))
                        + quicksum(manpowertwo[j] * buildcapa[j, t] for j in range(K)) <= 470)  # 470を労働力の上限と仮定

# 目的関数の設定
total_manpower = quicksum(manpowerone[k] * produce[k, t] for k in range(K) for t in range(T)) + quicksum(manpowertwo[k] * buildcapa[k, t] for k in range(K) for t in range(T))
model.setObjective(total_manpower, GRB.MAXIMIZE)

# 最適化の実行
model.optimize()

# 結果の出力
produce_result = [[produce[k, t].X for t in range(T)] for k in range(K)]
buildcapa_result = [[buildcapa[k, t].X for t in range(T)] for k in range(K)]
stockhold_result = [[stockhold[k, t].X for t in range(T+1)] for k in range(K)]

# 出力形式に合わせて結果を表示
output = {
    "produce": produce_result,
    "buildcapa": buildcapa_result,
    "stockhold": stockhold_result
}

print(output)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 63 rows, 48 columns and 265 nonzeros
Model fingerprint: 0x813e3bca
Coefficient statistics:
  Matrix range     [1e-01, 1e+00]
  Objective range  [1e-01, 6e-01]
  Bounds range     [0e+00, 0e+00]
  RHS range        [3e+01, 5e+02]
Presolve removed 43 rows and 18 columns
Presolve time: 0.00s
Presolved: 20 rows, 30 columns, 127 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    9.9512024e+03   6.678242e+03   0.000000e+00      0s
      20    2.3250000e+03   0.000000e+00   0.000000e+00      0s

Solved in 20 iterations and 0.01 seconds (0.00 work units)
Optimal objective  2.325000000e+03
{'produce': [[300.0, 0.0, 300.0, 90.0, 300.0], [350.0, 128.75, 350.0, 222.5, 350.0], [280.0, 280.0, 280.0, 280.0, 280.0]], 'buildcapa': [[0.0, 0.0, 322.5, 0.0, 270.0], [520.0, 1876.875, 0.0, 1466.25, 105.0], [0.0, 0.0, 0.0, 0.0, 0.0]], 'stockhold': [[150.0, 390.0, 30.0, 270.0, 0.0, 150.0], [80.0, 370.0, 88.74999999999997, 250.0, 62.5, 130.0], [100.0, 350.0, 320.0, 290.0, 260.0, 230.0]]}

起動停止問題

model_building_in_mathematical_programming (problem 16)

A number of power stations are committed to meeting the following electricity load demands over a day.
Period of a day has demand .
There are several types of generating unit available: the number of -th type units is .
Type generator has to work between a minimum level and a maximum level .
There is an hourly cost of running each type generator at minimum level .
In addition, there is an extra hourly cost for each megawatt at which a type unit is operated above the minimum level.
Starting up a type generator also involves a cost .

INPUT FORMAT:

{
    "demand": [demand_{t} for t in 1,...,T],
    "num": [num_{k} for in 1,...,K],
    "minlevel": [minlevel_{k} for in 1,...,K],
    "maxlevel": [minlevel_{k} for in 1,...,K],
    "runcost": [runcost_{k} for in 1,...,K],
    "extracost": [extracost_{k} for in 1,...,K],
    "startcost": [startcost_{k} for in 1,...,K]
}

OBJECTIVE: Which generators should be working in which periods of the day to minimize total cost?

OUTPUT INFO:

represents the number of generators of each type that are on in each period

OUTPUT FORMAT:

{
    "numon": [[numon_{k, t} for t in 1,...,T] for k in 1,...,K]
}

以下の問題を混合整数計画問題 (MIP) として定式化します。

問題の定式化

変数

\(P_{k,t}\): 時刻\(t\) におけるタイプ\(k\) の発電機の発電量 (連続変数)
\(U_{k,t}\): 時刻\(t\) におけるタイプ\(k\) の発電機の稼働状態 (バイナリ変数, 1=稼働, 0=停止)
\(S_{k,t}\): 時刻\(t\) におけるタイプ\(k\) の発電機の起動状態 (バイナリ変数, 1=起動, 0=非起動)

パラメータ

\(\text{demand}_t\): 時刻\(t\) の需要 (MW)
\(\text{num}_k\): タイプ\(k\) の発電機の台数
\(\text{minlevel}_k\): タイプ\(k\) の発電機の最低発電量 (MW)
\(\text{maxlevel}_k\): タイプ\(k\) の発電機の最大発電量 (MW)
\(\text{runcost}_k\): タイプ\(k\) の発電機が最低発電量で稼働するためのコスト (単位時間あたり)
\(\text{extracost}_k\): タイプ\(k\) の発電機が最低発電量を超えて稼働する場合の追加コスト (MWあたり)
\(\text{startcost}_k\): タイプ\(k\) の発電機の起動コスト

目的関数

総コストを最小化する。

\[ \text{Minimize} \sum_{t=1}^T \sum_{k=1}^K \left( \text{runcost}_k \cdot U_{k,t} + \text{extracost}_k \cdot (P_{k,t} - \text{minlevel}_k \cdot U_{k,t}) + \text{startcost}_k \cdot S_{k,t} \right) \]

制約条件

発電量制約: \[ \sum_{k=1}^K P_{k,t} \geq \text{demand}_t \quad \forall t \]
発電機の出力制約: \[ \text{minlevel}_k \cdot U_{k,t} \leq P_{k,t} \leq \text{maxlevel}_k \cdot U_{k,t} \quad \forall k, \forall t \]
発電機の台数制約: \[ \sum_{t=1}^T U_{k,t} \leq \text{num}_k \quad \forall k \]
発電機の起動状態制約: \[ S_{k,t} \geq U_{k,t} - U_{k,t-1} \quad \forall k, \forall t > 1 \]

\[ S_{k,1} \geq U_{k,1} \quad \forall k \]

入力形式

{
    "demand": [demand_t for t in 1,...,T],
    "num": [num_k for k in 1,...,K],
    "minlevel": [minlevel_k for k in 1,...,K],
    "maxlevel": [maxlevel_k for k in 1,...,K],
    "runcost": [runcost_k for k in 1,...,K],
    "extracost": [extracost_k for k in 1,...,K],
    "startcost": [startcost_k for k in 1,...,K]
}

出力形式

{
    "numon": [[numon_{k, t} for t in 1,...,T] for k in 1,...,K]
}

問題を読み違えて実行不能になる．　

#間違い！
from gurobipy import Model, GRB

# 入力データ
data = {
    "demand": [15000, 30000, 25000, 40000, 27000],
    "num": [12, 10, 5],
    "minlevel": [850, 1250, 1500],
    "maxlevel": [2000, 1750, 4000],
    "runcost": [1000, 2600, 3000],
    "extracost": [2.0, 1.3, 3.0],
    "startcost": [2000, 1000, 500]
}

demand = data["demand"]
num = data["num"]
minlevel = data["minlevel"]
maxlevel = data["maxlevel"]
runcost = data["runcost"]
extracost = data["extracost"]
startcost = data["startcost"]

T = len(demand)  # 時間の数
K = len(num)     # 発電機のタイプの数

# モデルの作成
model = Model("PowerGeneration")

# 変数の作成
P = model.addVars(K, T, name="P", lb=0, vtype=GRB.CONTINUOUS)
U = model.addVars(K, T, name="U", vtype=GRB.INTEGER)
S = model.addVars(K, T, name="S", vtype=GRB.BINARY)

# 目的関数の設定
cost = sum(runcost[k] * U[k, t] + extracost[k] * (P[k, t] - minlevel[k] * U[k, t]) + startcost[k] * S[k, t]
           for k in range(K) for t in range(T))
model.setObjective(cost, GRB.MINIMIZE)

# 制約条件の設定
# 需要制約
for t in range(T):
    model.addConstr(sum(P[k, t] for k in range(K)) >= demand[t], name=f"demand_{t}")

# 発電量の制約
for k in range(K):
    for t in range(T):
        model.addConstr(minlevel[k] * U[k, t] <= P[k, t], name=f"minlevel_{k}_{t}")
        model.addConstr(P[k, t] <= maxlevel[k] * U[k, t], name=f"maxlevel_{k}_{t}")

# 発電機の台数制約
for k in range(K):
    model.addConstr(sum(U[k, t] for t in range(T)) <= num[k], name=f"num_{k}")

# 起動状態の制約
for k in range(K):
    for t in range(1, T):
        model.addConstr(S[k, t] >= U[k, t] - U[k, t-1], name=f"start_{k}_{t}")
    model.addConstr(S[k, 1] >= U[k, 1], name=f"start_{k}_1")

# 最適化の実行
model.optimize()

# 結果の出力
numon = [[int(U[k, t].X) for t in range(T)] for k in range(K)]

# 出力形式に合わせて結果を表示
output = {
    "numon": numon
}
print(output)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 53 rows, 45 columns and 132 nonzeros
Model fingerprint: 0xd809e7b2
Variable types: 15 continuous, 30 integer (15 binary)
Coefficient statistics:
  Matrix range     [1e+00, 4e+03]
  Objective range  [1e+00, 2e+03]
  Bounds range     [1e+00, 1e+00]
  RHS range        [5e+00, 4e+04]
Presolve removed 0 rows and 3 columns
Presolve time: 0.00s

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 16 available processors)

Solution count 0

Model is infeasible or unbounded
Best objective -, best bound -, gap -

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[17], line 63
     60 model.optimize()
     62 # 結果の出力
---> 63 numon = [[int(U[k, t].X) for t in range(T)] for k in range(K)]
     65 # 出力形式に合わせて結果を表示
     66 output = {
     67     "numon": numon
     68 }

Cell In[17], line 63, in <listcomp>(.0)
     60 model.optimize()
     62 # 結果の出力
---> 63 numon = [[int(U[k, t].X) for t in range(T)] for k in range(K)]
     65 # 出力形式に合わせて結果を表示
     66 output = {
     67     "numon": numon
     68 }

Cell In[17], line 63, in <listcomp>(.0)
     60 model.optimize()
     62 # 結果の出力
---> 63 numon = [[int(U[k, t].X) for t in range(T)] for k in range(K)]
     65 # 出力形式に合わせて結果を表示
     66 output = {
     67     "numon": numon
     68 }

File src/gurobipy/var.pxi:125, in gurobipy.Var.__getattr__()

File src/gurobipy/var.pxi:153, in gurobipy.Var.getAttr()

File src/gurobipy/attrutil.pxi:100, in gurobipy.__getattr()

AttributeError: Unable to retrieve attribute 'X'

修正プロンプト

発電機は３つのタイプがあり，各タイプの台数上限が"num": [12, 10, 5] に保管されています．
つまり，0-1変数SとUは，各々の発電機に対して定義しなければなりません．コードをそのように修正してください．

一応実行可能にはなるが，まだおかしなコードを出力する．

from gurobipy import Model, GRB

# 入力データ
data = {
    "demand": [15000, 30000, 25000, 40000, 27000],
    "num": [12, 10, 5],
    "minlevel": [850, 1250, 1500],
    "maxlevel": [2000, 1750, 4000],
    "runcost": [1000, 2600, 3000],
    "extracost": [2.0, 1.3, 3.0],
    "startcost": [2000, 1000, 500]
}

demand = data["demand"]
num = data["num"]
minlevel = data["minlevel"]
maxlevel = data["maxlevel"]
runcost = data["runcost"]
extracost = data["extracost"]
startcost = data["startcost"]

T = len(demand)  # 時間の数
K = len(num)     # 発電機のタイプの数

# モデルの作成
model = Model("PowerGeneration")

# 変数の作成
P = model.addVars(K, max(num), T, name="P", lb=0, vtype=GRB.CONTINUOUS)
U = model.addVars(K, max(num), T, name="U", vtype=GRB.BINARY)
S = model.addVars(K, max(num), T, name="S", vtype=GRB.BINARY)

# 目的関数の設定
cost = sum(runcost[k] * U[k, n, t] + extracost[k] * (P[k, n, t] - minlevel[k] * U[k, n, t]) + startcost[k] * S[k, n, t]
           for k in range(K) for n in range(num[k]) for t in range(T))
model.setObjective(cost, GRB.MINIMIZE)

# 制約条件の設定
# 需要制約
for t in range(T):
    model.addConstr(sum(P[k, n, t] for k in range(K) for n in range(num[k])) >= demand[t], name=f"demand_{t}")

# 発電量の制約
for k in range(K):
    for n in range(num[k]):
        for t in range(T):
            model.addConstr(minlevel[k] * U[k, n, t] <= P[k, n, t], name=f"minlevel_{k}_{n}_{t}")
            model.addConstr(P[k, n, t] <= maxlevel[k] * U[k, n, t], name=f"maxlevel_{k}_{n}_{t}")

# 起動状態の制約
for k in range(K):
    for n in range(num[k]):
        for t in range(1, T):
            model.addConstr(S[k, n, t] >= U[k, n, t] - U[k, n, t-1], name=f"start_{k}_{n}_{t}")
        model.addConstr(S[k, n, 1] >= U[k, n, 1], name=f"start_{k}_{n}_1")

# 最適化の実行
model.optimize()

# 結果の出力
numon = [[int(U[k, n, t].X) for t in range(T)] for k in range(K) for n in range(num[k])]

# 出力形式に合わせて結果を表示
output = {
    "numon": numon
}
print(output)

Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Xeon(R) W-2140B CPU @ 3.20GHz
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads

Optimize a model with 410 rows, 540 columns and 1053 nonzeros
Model fingerprint: 0xd9b74fb5
Variable types: 180 continuous, 360 integer (360 binary)
Coefficient statistics:
  Matrix range     [1e+00, 4e+03]
  Objective range  [1e+00, 2e+03]
  Bounds range     [1e+00, 1e+00]
  RHS range        [2e+04, 4e+04]
Found heuristic solution: objective 328050.00000
Presolve removed 109 rows and 243 columns
Presolve time: 0.03s
Presolved: 301 rows, 297 columns, 783 nonzeros
Found heuristic solution: objective 319244.00000
Variable types: 0 continuous, 297 integer (189 binary)
Found heuristic solution: objective 319238.00000

Root relaxation: objective 2.558429e+05, 216 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 255842.857    0    6 319238.000 255842.857  19.9%     -    0s
H    0     0                    257011.50000 255842.857  0.45%     -    0s
H    0     0                    256500.00000 255842.857  0.26%     -    0s
     0     0 256056.598    0   17 256500.000 256056.598  0.17%     -    0s
     0     0 256056.598    0    6 256500.000 256056.598  0.17%     -    0s
     0     0 256056.598    0   28 256500.000 256056.598  0.17%     -    0s
     0     0 256250.351    0   24 256500.000 256250.351  0.10%     -    0s
     0     0 256254.098    0   15 256500.000 256254.098  0.10%     -    0s
     0     0 256254.098    0   27 256500.000 256254.098  0.10%     -    0s
     0     0 256254.098    0   24 256500.000 256254.098  0.10%     -    0s
     0     0 256254.098    0   30 256500.000 256254.098  0.10%     -    0s
     0     0 256254.098    0   32 256500.000 256254.098  0.10%     -    0s
     0     0 256254.098    0    8 256500.000 256254.098  0.10%     -    0s
     0     0 256260.832    0    8 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0   12 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0   24 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0    8 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0   12 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0    8 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0    8 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0    2 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0   15 256500.000 256260.832  0.09%     -    0s
     0     0 256260.832    0   18 256500.000 256260.832  0.09%     -    0s
     0     0 256291.667    0    5 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0    5 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0   19 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0    7 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0   19 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0    7 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0    7 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0    3 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0   12 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0    5 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0    8 256500.000 256291.667  0.08%     -    0s
     0     0 256291.667    0   20 256500.000 256291.667  0.08%     -    0s

Cutting planes:
  Gomory: 2
  MIR: 4

Explored 1 nodes (814 simplex iterations) in 0.24 seconds (0.04 work units)
Thread count was 16 (of 16 available processors)

Solution count 7: 256500 256500 256500 ... 328050

Optimal solution found (tolerance 1.00e-04)
Best objective 2.565000000000e+05, best bound 2.565000000000e+05, gap 0.0000%
{'numon': [[1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 1, 1, 1, 1], [1, 0, 0, 0, 0], [0, 1, 1, 1, 1], [0, 1, 1, 1, 1], [0, 1, 1, 1, 1], [0, 1, 1, 1, 1], [0, 1, 1, 1, 1], [0, 1, 1, 1, 1], [0, 1, 1, 1, 1], [1, 1, 1, 1, 1], [0, 0, 0, 1, 1], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 0], [0, 0, 0, 0, 0]]}

結論

ChatGPT4oは，きちんとした日本語で書けば定式化とコードを出してくれる（ことが多い）．

うまくいかない場合には，修正のプロンプトを正しく書けば，（不格好だが）直してくれる（こともある）．

実用化のためには，データ不足や誤った入力に対して，正しい情報を引き出すための返答をするエージェントなどの開発が必要になる．

問題 クッキーとケーキ

OptiMUSによるコード

Claudeによるコード

ChatGPT4o によるコード

手計算による解の概要

嘘つき島パズル

OptiMUSの例題

生産計画

問題の定式化

変数

制約

目的関数

コードの説明

スポーツ用品の生産

INPUT FORMAT:

OUTPUT INFO:

OUTPUT FORMAT:

問題の定式化

変数

パラメータ

目的関数

制約条件

入力形式

出力形式

食品製造問題1

定式化の詳細

パラメータ

変数

目的関数

制約条件

説明

結果

修正プロンプトと返答

食品製造問題2

モデルの定式化

追加パラメータ

追加変数

制約条件

工場計画問題1

入力形式:

目的:

出力情報:

出力形式:

定式化

変数

パラメータ

目的関数

制約

工場計画問題2

修正ポイント

発電最適化

問題の定式化

変数

パラメータ

目的関数

制約条件

出力フォーマット

人員計画

INPUT FORMAT:

OUTPUT INFO:

OUTPUT FORMAT:

定式化

変数

パラメータ

目的関数

制約

入力フォーマット

出力フォーマット

産業連関表

INPUT FORMAT:

OUTPUT INFO:

モデルの定式化

説明

起動停止問題

INPUT FORMAT:

OUTPUT INFO:

OUTPUT FORMAT:

問題の定式化

変数

パラメータ

問題クッキーとケーキ