从酸奶厂到宣传片：用Python+PuLP搞定线性规划与排序问题（附完整代码）

当你第一次翻开数学建模教材，看到那些抽象的资源分配、任务调度问题时，是否曾感到无从下手？本文将以两个典型问题为例——酸奶厂生产优化和宣传片制作调度，手把手带你用Python的PuLP库和经典排序算法，将纸上模型转化为可执行的代码解决方案。我们会从问题理解、模型构建到代码实现完整走一遍，并提供可直接运行的Jupyter Notebook代码片段。

1. 酸奶厂生产优化：线性规划实战

假设你是一家乳制品厂的生产经理，需要决定A、B、C三种酸奶的产量以最大化利润。工厂有以下资源限制：

• 原料甲：400千克
• 原料乙：50千克
• 设备工时：1200小时

每种酸奶的单位消耗和利润如下：

产品	设备工时	原料甲(kg)	原料乙(kg)	利润(元)
A	3	1	0.2	5.5
B	5	1	0.3	10
C	4	1	0.15	6

1.1 建立数学模型

这是一个典型的线性规划问题，目标函数和约束条件可表示为：

最大化利润 ：

5.5*A + 10*B + 6*C

约束条件 ：

3*A + 5*B + 4*C ≤ 1200  (设备工时)
1*A + 1*B + 1*C ≤ 400   (原料甲)
0.2*A + 0.3*B + 0.15*C ≤ 50 (原料乙)
A, B, C ≥ 0

1.2 Python实现：PuLP库应用

from pulp import *

# 创建问题实例
prob = LpProblem("Yogurt_Production", LpMaximize)

# 定义决策变量
A = LpVariable("A", lowBound=0, cat='Continuous')
B = LpVariable("B", lowBound=0, cat='Continuous')
C = LpVariable("C", lowBound=0, cat='Continuous')

# 定义目标函数
prob += 5.5*A + 10*B + 6*C, "Total Profit"

# 添加约束条件
prob += 3*A + 5*B + 4*C <= 1200, "Machine Time"
prob += 1*A + 1*B + 1*C <= 400, "Material X"
prob += 0.2*A + 0.3*B + 0.15*C <= 50, "Material Y"

# 求解问题
prob.solve()

# 输出结果
print(f"Status: {LpStatus[prob.status]}")
print(f"Optimal Production Plan:")
print(f"A: {value(A):.2f} units")
print(f"B: {value(B):.2f} units")
print(f"C: {value(C):.2f} units")
print(f"Total Profit: {value(prob.objective):.2f} yuan")

提示：PuLP默认使用CBC求解器，对于小型问题足够高效。如需处理更大规模问题，可考虑商业求解器如Gurobi或CPLEX。

1.3 结果分析与可视化

运行上述代码后，我们得到了最优生产方案。为了更直观地理解资源分配，可以绘制资源使用情况的堆叠条形图：

import matplotlib.pyplot as plt

resources = ['Machine Time', 'Material X', 'Material Y']
total = [1200, 400, 50]
used = [
    3*value(A) + 5*value(B) + 4*value(C),
    1*value(A) + 1*value(B) + 1*value(C),
    0.2*value(A) + 0.3*value(B) + 0.15*value(C)
]

plt.figure(figsize=(10, 6))
bars = plt.bar(resources, used)
plt.bar(resources, [t-u for t,u in zip(total, used)], bottom=used, color='lightgray')
plt.ylabel('Hours/Kilograms')
plt.title('Resource Utilization')
for bar in bars:
    height = bar.get_height()
    plt.text(bar.get_x() + bar.get_width()/2., height,
             f'{height/total[bars.index(bar)]*100:.1f}%',
             ha='center', va='bottom')
plt.show()

2. 宣传片制作调度：带惩罚的排序问题

现在考虑第二个问题：一家公司需要制作三个宣传片，每个宣传片有制作周期、截止日期和延迟罚金。目标是安排制作顺序使总罚金最小。

给定数据：

宣传片	制作周期(Pi)	截止日期(Di)	每日罚金(Ci)
1	5	10	2
2	8	15	3
3	4	12	1

2.1 问题建模

这是一个典型的单机调度问题，我们需要：

1. 确定宣传片的制作顺序
2. 计算每个宣传片的完成时间
3. 计算总延迟罚金

数学模型的关键变量：

• 排序变量：表示宣传片的制作顺序
• 完成时间：前一个宣传片完成后开始下一个
• 延迟时间：完成时间减去截止日期（若为正）

2.2 Python实现：穷举法求解

对于小规模问题（n≤10），穷举所有排列是可行的：

from itertools import permutations

# 输入数据
projects = {
    1: {'duration': 5, 'deadline': 10, 'penalty': 2},
    2: {'duration': 8, 'deadline': 15, 'penalty': 3},
    3: {'duration': 4, 'deadline': 12, 'penalty': 1}
}

# 计算给定顺序的总罚金
def calculate_penalty(order):
    time = 0
    total_penalty = 0
    for project in order:
        time += projects[project]['duration']
        delay = max(0, time - projects[project]['deadline'])
        total_penalty += delay * projects[project]['penalty']
    return total_penalty

# 穷举所有可能的顺序
best_order = None
min_penalty = float('inf')
for order in permutations([1, 2, 3]):
    current_penalty = calculate_penalty(order)
    if current_penalty < min_penalty:
        min_penalty = current_penalty
        best_order = order

print(f"Optimal Order: {best_order}")
print(f"Minimum Total Penalty: {min_penalty}")

2.3 更高效的解法：EDD与SPT规则

对于大规模问题，穷举法不可行。我们可以考虑启发式规则：

# 最早截止日期优先(EDD)
edd_order = sorted(projects.keys(), key=lambda x: projects[x]['deadline'])
edd_penalty = calculate_penalty(edd_order)

# 最短处理时间优先(SPT)
spt_order = sorted(projects.keys(), key=lambda x: projects[x]['duration'])
spt_penalty = calculate_penalty(spt_order)

print(f"EDD Order: {edd_order}, Penalty: {edd_penalty}")
print(f"SPT Order: {spt_order}, Penalty: {spt_penalty}")

注意：实际应用中，可能需要结合多种规则或使用更高级的算法如分支定界法。

3. 进阶应用：利润随产量变化的非线性规划

回到酸奶厂问题，现在假设产品A的利润是其产量的函数：p = 6 - 0.01x。这变成了一个非线性规划问题。

3.1 数学模型调整

新的目标函数：

(6 - 0.01*A)*A + 10*B + 6*C

约束条件保持不变。

3.2 Python实现：SciPy优化

from scipy.optimize import minimize

def objective(x):
    A, B, C = x
    return -((6 - 0.01*A)*A + 10*B + 6*C)  # 负号因为minimize

# 约束条件
constraints = [
    {'type': 'ineq', 'fun': lambda x: 1200 - (3*x[0] + 5*x[1] + 4*x[2])},
    {'type': 'ineq', 'fun': lambda x: 400 - (x[0] + x[1] + x[2])},
    {'type': 'ineq', 'fun': lambda x: 50 - (0.2*x[0] + 0.3*x[1] + 0.15*x[2])}
]

# 变量边界
bounds = [(0, None), (0, None), (0, None)]

# 初始猜测
x0 = [100, 100, 100]

# 求解
solution = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints)

print(f"Optimal Solution: A={solution.x[0]:.2f}, B={solution.x[1]:.2f}, C={solution.x[2]:.2f}")
print(f"Maximum Profit: {-solution.fun:.2f}")

4. 工程实践中的注意事项

在实际项目中应用这些技术时，有几个关键点需要考虑：

1. 数据准备与验证 ：

   • 确保所有输入参数准确无误
   • 对约束条件进行敏感性分析
   • 检查模型假设是否合理（如线性关系假设）

2. 性能考量 ：

   • 对于大规模问题，选择合适的求解器
   • 考虑使用并行计算或分布式优化技术
   • 可能需要简化模型以提高求解效率

3. 结果解释与应用 ：

   • 理解解的稳定性（微小变化是否导致解大幅变动）
   • 考虑实施中的实际限制（如最小生产批量）
   • 准备应对意外情况的备选方案

以下是一个简单的敏感性分析示例，展示设备工时变化对最优利润的影响：

import numpy as np
import matplotlib.pyplot as plt

machine_hours = np.linspace(1000, 1500, 20)
profits = []

for hours in machine_hours:
    prob = LpProblem("Sensitivity_Analysis", LpMaximize)
    A = LpVariable("A", lowBound=0)
    B = LpVariable("B", lowBound=0)
    C = LpVariable("C", lowBound=0)
    prob += 5.5*A + 10*B + 6*C
    prob += 3*A + 5*B + 4*C <= hours
    prob += 1*A + 1*B + 1*C <= 400
    prob += 0.2*A + 0.3*B + 0.15*C <= 50
    prob.solve()
    profits.append(value(prob.objective))

plt.plot(machine_hours, profits)
plt.xlabel('Available Machine Hours')
plt.ylabel('Maximum Profit')
plt.title('Sensitivity Analysis: Machine Hours vs Profit')
plt.grid(True)
plt.show()