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

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

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

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

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

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

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

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. 宣传片制作调度：带惩罚的排序问题

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

给定数据：

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()