两阶段生产计划问题

在生产计划问题中，决策者必须决定如何购买原材料、劳动力和其他资源，以生产最终产品来最大化利润。

在本案例研究中，一家公司（GTC）生产扳手和钳子，受限于钢材的供应、机器能力（成型和组装）、劳动力以及市场需求。GTC希望确定购买多少钢材。问题的复杂性在于，当前可用的组装能力和产品对盈利的贡献是未知的，但将在下一期开始时得知。

因此，在这个时期，GTC必须：

• 确定需要购买多少钢材。

在下个周期开始时，GTC 将首先确定可用的装配产能以及每单位扳手和钳子的收入，然后 GTC 将决定

• 需要生产多少扳手和钳子。

不确定性表现为四种可能的情景之一，每种情景的概率相等。

我们首先导入 PuLP 包。

import pulp

接下来，我们将读取数据。这里，我们将数据读取为向量。在实际使用中，这可能是从数据库中读取的。首先，是不随场景变化的数据元素。这些元素各有两种值，一种对应扳手，另一种对应钳子。

products = ["wrenches", "pliers"]
price = [130, 100]
steel = [1.5, 1]
molding = [1, 1]
assembly = [0.3, 0.5]
capsteel = 27
capmolding = 21
LB = [0, 0]
capacity_ub = [15, 16]
steelprice = 58

下一组参数是与四种场景相对应的参数。

scenarios = [0, 1, 2, 3]
pscenario = [0.25, 0.25, 0.25, 0.25]
wrenchearnings = [160, 160, 90, 90]
plierearnings = [100, 100, 100, 100]
capassembly = [8, 10, 8, 10]

接下来，我们将创建代表产品和场景组合的列表。这些列表稍后将用于创建参数的字典。

production = [(j, i) for j in scenarios for i in products]
pricescenario = [[wrenchearnings[j], plierearnings[j]] for j in scenarios]
priceitems = [item for sublist in pricescenario for item in sublist]

接下来，我们使用 dict(zip(…)) 将这些列表转换为字典。这样做是为了我们可以通过有意义的名称来引用参数。

price_dict = dict(zip(production, priceitems))
capacity_dict = dict(zip(products, capacity_ub * 4))
steel_dict = dict(zip(products, steel))
molding_dict = dict(zip(products, molding))
assembly_dict = dict(zip(products, assembly))

要定义我们的决策变量，我们使用函数 pulp.LpVariable.dicts()，该函数创建带有相关索引值的字典。

production_vars = pulp.LpVariable.dicts(
    "production", (scenarios, products), lowBound=0, cat="Continuous"
)
steelpurchase = pulp.LpVariable("steelpurchase", lowBound=0, cat="Continuous")

我们创建 LpProblem 然后制定目标函数。注意这是一个最大化问题，因为目标是最大化净收入。

gemstoneprob = pulp.LpProblem("The Gemstone Tool Problem", pulp.LpMaximize)

目标函数使用 pulp.lpSum() 函数指定。注意，它是通过 += 添加到问题中的。

gemstoneprob += (
    pulp.lpSum(
        [
            pscenario[j] * (price_dict[(j, i)] * production_vars[j][i])
            for (j, i) in production
        ]
        - steelpurchase * steelprice
    ),
    "Total cost",
)

然后我们加入约束条件。这里的约束条件是基于场景和产品设置的，并使用 for i in list: 表示法指定。在每个约束条件中，求和表达式使用列表推导式表示。注意，约束条件与目标函数不同，因为每个约束条件都以逻辑比较（通常是 <= 或 >=，但也可以是 ==）结束。最后，这里，文件为每个约束条件命名，名称包括约束条件适用的特定场景或产品。

for j in scenarios:
    gemstoneprob += pulp.lpSum(
        [steel_dict[i] * production_vars[j][i] for i in products]
    ) - steelpurchase <= 0, ("Steel capacity" + str(j))
    gemstoneprob += pulp.lpSum(
        [molding_dict[i] * production_vars[j][i] for i in products]
    ) <= capmolding, ("molding capacity" + str(j))
    gemstoneprob += pulp.lpSum(
        [assembly_dict[i] * production_vars[j][i] for i in products]
    ) <= capassembly[j], ("assembly capacity" + str(j))
    for i in products:
        gemstoneprob += production_vars[j][i] <= capacity_dict[i], (
            "capacity " + str(i) + str(j)
        )

完整的文件可以在这里找到 Two_stage_Stochastic_GemstoneTools.py

r"""
Author: Louis Luangkesorn <lugerpitt@gmail.com> 2019
https://github.com/lluang

Title: Gemstone Optimization problem

Problem taken from Data, Models, and Decisions by Bertsimas and Freund, 4th Edition
DMD 7.2

## 2 stage problem

- **Scenarios:** $s \in S = (1, 2, 3, 4)$
- **Probability scenario occuring:** $p^s$
- **Cost of steel:** $c$
- **Total steel:** $cap_{steel}$
- **Total molding and assembly hours:** $cap_{molding}, cap_{assembly}^s$
- **Wrench and plier earnings by scenario:** $w^s, p^s$
- **Max demand wrenches and pliers:** $UB_w, UB_p$
- **Decision variables**
  - $(W_{t+1}^s, P_{t+1}^s)$
- **Objective**   $Max \sum_s (p^s * (w^s W_{t+1}^s + p^s P_{t+1}^s) - c$
- **Steel Constraint:** $1.5W_{t+1}^1 + P_{t+1}^1 - C \le 0$
- **Molding Constraint:** $W_{t+1}^1 + P_{t+1}^1 \le cap_{molding}$
- **Assembly Constraint:** $0.3 W_{t+1}^1 + 0.5 P_{t+1}^1  \le cap_{molding}^s$
- **Demand Limit W:** $W \le UB_w$
- **Demand Limit P:** $P \le UB_p$
- **Nonnegativity:** $W, P \ge 0$
"""

import pulp

# parameters
products = ["wrenches", "pliers"]
price = [130, 100]
steel = [1.5, 1]
molding = [1, 1]
assembly = [0.3, 0.5]
capsteel = 27
capmolding = 21
LB = [0, 0]
capacity_ub = [15, 16]
steelprice = 58
scenarios = [0, 1, 2, 3]
pscenario = [0.25, 0.25, 0.25, 0.25]
wrenchearnings = [160, 160, 90, 90]
plierearnings = [100, 100, 100, 100]
capassembly = [8, 10, 8, 10]

production = [(j, i) for j in scenarios for i in products]
pricescenario = [[wrenchearnings[j], plierearnings[j]] for j in scenarios]
priceitems = [item for sublist in pricescenario for item in sublist]

# create dictionaries for the parameters
price_dict = dict(zip(production, priceitems))
capacity_dict = dict(zip(products, capacity_ub * 4))
steel_dict = dict(zip(products, steel))
molding_dict = dict(zip(products, molding))
assembly_dict = dict(zip(products, assembly))

# Create variables and parameters as dictionaries
production_vars = pulp.LpVariable.dicts(
    "production", (scenarios, products), lowBound=0, cat="Continuous"
)
steelpurchase = pulp.LpVariable("steelpurchase", lowBound=0, cat="Continuous")

# Create the 'gemstoneprob' variable to specify
gemstoneprob = pulp.LpProblem("The Gemstone Tool Problem", pulp.LpMaximize)

# The objective function is added to 'gemstoneprob' first
gemstoneprob += (
    pulp.lpSum(
        [
            pscenario[j] * (price_dict[(j, i)] * production_vars[j][i])
            for (j, i) in production
        ]
        - steelpurchase * steelprice
    ),
    "Total cost",
)

for j in scenarios:
    gemstoneprob += pulp.lpSum(
        [steel_dict[i] * production_vars[j][i] for i in products]
    ) - steelpurchase <= 0, ("Steel capacity" + str(j))
    gemstoneprob += pulp.lpSum(
        [molding_dict[i] * production_vars[j][i] for i in products]
    ) <= capmolding, ("molding capacity" + str(j))
    gemstoneprob += pulp.lpSum(
        [assembly_dict[i] * production_vars[j][i] for i in products]
    ) <= capassembly[j], ("assembly capacity" + str(j))
    for i in products:
        gemstoneprob += production_vars[j][i] <= capacity_dict[i], (
            "capacity " + str(i) + str(j)
        )

# Print problem
print(gemstoneprob)

# The problem data is written to an .lp file
gemstoneprob.writeLP("gemstoneprob.lp")
# The problem is solved using PuLP's choice of Solver
gemstoneprob.solve()
# The status of the solution is printed to the screen
print("Status:", pulp.LpStatus[gemstoneprob.status])

# OUTPUT

# Each of the variables is printed with it's resolved optimum value
for v in gemstoneprob.variables():
    print(v.name, "=", v.varValue)
production = [v.varValue for v in gemstoneprob.variables()]

# The optimised objective function value is printed to the console
print("Total price = ", pulp.value(gemstoneprob.objective))