运筹学练习Python精解——线性规划

练习1

考虑下面线性规划，其中$c_1$的值不能确定，使用图解法找出所对应所有不同$c_1(-\infty <c_1<\infty)$值的$(x_1,x_2)$最优解。

$$\max \quad Z = c_1 x_1 + x_2 \\ \text{s.t.} \quad \begin{cases} x_1 + x_2 \leq 6 \\ x_1 + 2x_2 \leq 10 \\ x_1 \geq 0 \\ x_2 \geq 0 \end{cases}$$

1.1 约束条件的交点

• 求交点

$x_1 + x_2 = 6x$ 和 $x_1 + 2x_2 = 10$ 的交点：

$$\begin{cases} x_1 + x_2 = 6 \\ x_1 + 2x_2 = 10 \end{cases}$$

所以交点是(2,4)。
$x_1 + x_2 = 6$ 和 $x_1 = 0$的交点：

$$\begin{cases} x_1 = 0 \\ 0 + x_2 = 6 \implies x_2 = 6 \end{cases}$$

但 (0,6) 不满足约束$x_1 + 2x_2 \leq 10$，因此它不是可行解。
$x_1 + 2x_2 = 10$ 和 $x_1 = 0$的交点：

$$\begin{cases} x_1 = 0 \\ 0 + 2x_2 = 10 \implies x_2 = 5 \end{cases}$$

所以交点是(0,5)。
$x_2 = 0$ 和 $x_1 + x_2 = 6$ 的交点：

$$\begin{cases} x_2 = 0 \\ x_1 + 0 = 6 \implies x_1 = 6 \end{cases}$$

所以交点是(6,0)。
$x_2 = 0$ 和 $x_1 + 2x_2 = 10$ 的交点：

$$\begin{cases} x_2 = 0 \\ x_1 + 0 = 10 \implies x_1 = 10 \end{cases}$$

但 (10,0) 不满足约束 $x_1 + x_2 \leq 6$，因此它不是可行解。

最终的可行区域顶点为：(2,4)；0,5)；(6,0)，得到目标函数 $Z = c_1 x_1 + x_2$在各个顶点的值为：在 (2,4)：$Z = 2c_1 + 4$ ；在 (0,5)：Z = 5 ；在 (6,0)：$Z = 6c_1$。

1.2 比较顶点的目标函数值

• 比较 (2,4) 和 (0,5) 的目标函数值
  我们要确定$c_1$ 的取值范围使得在这两个点上目标函数值相等：

$$2c_1 + 4 = 5 \implies 2c_1 = 1 \implies c_1 = 0.5$$

当 $c_1 < 0.5$ ，(0,5)的目标函数值大于 (2,4) ；
当 $c_1 > 0.5$，(2,4) 的目标函数值大于 (0,5) ；
当 $c_1 = 0.5$，在 (0,5) 和 (2,4) 之间的任意点都是最优解。

• 比较 (2,4) 和 (6,0) 的目标函数值

我们要确定$c_1$ 的取值范围使得在这两个点上目标函数值相等：

$$2c_1 + 4 = 6c_1 \implies 4 = 4c_1 \implies c_1 = 1$$

当 $c_1 < 1$， (2,4) 的目标函数值大于 (6,0) ；
当 $c_1 > 1$，(6,0) 的目标函数值大于 (2,4) ；
当 $c_1 = 1$，在 (2,4) 和 (6,0) 之间的任意点都是最优解。

1.3 汇总结果

通过上面的比较，我们可以总结：

$$(x_1, x_2) = \begin{cases} (0, 5) & \text{if } c_1 < 0.5 \\ \text{任意点在 } (0, 5) \text{ 和 } (2, 4) \text{之间} & \text{if } c_1 = 0.5 \\ (2, 4) & \text{if } 0.5 < c_1 < 1 \\ \text{任意点在 } (2, 4) \text{ 和 } (6, 0) \text{之间} & \text{if } c_1 = 1 \\ (6, 0) & \text{if } c_1 > 1 \end{cases}$$

这样，我们就得到了在不同 $c_1$ 取值范围下的最优解。

1.4 Python展示

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog

# Define the constraints
x1 = np.linspace(0, 10, 400)
x2_1 = 6 - x1
x2_2 = (10 - x1) / 2

# Plot the feasible region
plt.figure(figsize=(8, 6))
plt.plot(x1, x2_1, label=r'$x_1 + x_2 \leq 6$')
plt.plot(x1, x2_2, label=r'$x_1 + 2x_2 \leq 10$')
plt.xlim((0, 10))
plt.ylim((0, 6))

# Fill feasible region
plt.fill_between(x1, np.minimum(x2_1, x2_2), 0, where=(np.minimum(x2_1, x2_2)>0), color='grey', alpha=0.5)

# Labels and title
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='black', lw=1)
plt.legend()
plt.title('Feasible Region')

# Solve for different values of c1
c1_values = [-1, 0, 1, 2] # You can change these values to any desired range
solutions = []

for c1 in c1_values:
    c = [-c1, -1]  # Coefficients in objective function for maximization problem
    A = [[1, 1], [1, 2]]
    b = [6, 10]
    x0_bounds = (0, None)
    x1_bounds = (0, None)
    res = linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds], method='highs')
    solutions.append((res.x[0], res.x[1], c1))

    # Plot the optimal points
    plt.plot(res.x[0], res.x[1], 'o', label=f'Optimal for $c_1 = {c1}$')

plt.legend()
plt.show()

for sol in solutions:
    print(f"For c1 = {sol[2]}: Optimal solution (x1, x2) = ({sol[0]:.2f}, {sol[1]:.2f})")

练习2

某企业有三个车间生产同一种产品。每件产品由四个零件 1 和三个零件 2 组成。两个零件需耗用两种原材料 A 和 B。已知这两种原材料的供应量分别为 300kg 和 500kg。由于三个车间拥有的设备及工艺条件不同，每个工班原材料耗用量和零件产量也不同。见下表(三个车间每班用料和生产情况),问三个车间就各开多少个班，才能使该产品的配套数达到最大？

车间	A 材料	B 材料	零件1	零件2
一车间	8	6	7	5
二车间	5	9	6	9
三车间	3	8	8	4

2.1数学模型

有三个车间，每个车间的每班原材料耗用量和零件产量不同。我们的目标是确定每个车间开多少个班，以使得产品的配套数（即最终的产品数）最大。每个产品由4个零件1和3个零件2组成。我们有两个原材料A和B，供应量分别为300kg和500kg。

• 决策变量：
  • $x_1$ 表示一车间的班数
  • $x_2$​ 表示二车间的班数
  • $x_3$​ 表示三车间的班数
• 目标函数： 我们要最大化产品的配套数。每个产品需要4个零件1和3个零件2。假设最终产品的配套数为P，那么我们可以定义一个辅助变量来表示该目标：
  $$P = \max \left( \frac{7x_1 + 6x_2 + 8x_3}{4}, \frac{5x_1 + 9x_2 + 4x_3}{3} \right)$$

• 约束条件：
  • 原材料A的供应量不能超过300kg： $8x_1 + 5x_2 + 3x_3 \leq 300$
  • 原材料B的供应量不能超过500kg： $6x_1 + 9x_2 + 8x_3 \leq 500$
  • 零件1和零件2必须满足产品配套需求： $4P \leq 7x_1 + 6x_2 + 8x_3 \quad​ 3P \leq 5x_1 + 9x_2 + 4x_3$
  • 非负约束： $x_1 \geq 0, \quad x_2 \geq 0, \quad x_3 \geq 0$
• 求解方法： 使用线性规划求解工具，如Python的SciPy库或其他优化工具，求解这个线性规划问题。

$$\begin{align} \text{max} \quad & P \\ \text{subject to} \quad & 7x_1 + 6x_2 + 8x_3 - 4P \geq 0 \\ & 5x_1 + 9x_2 + 4x_3 - 3P \geq 0 \\ & 8x_1 + 5x_2 + 3x_3 \leq 300 \\ & 6x_1 + 9x_2 + 8x_3 \leq 500 \\ & x, x_2, x_3, P \geq 0 \end{align}$$

2.2 Python求解

我们使用Python的SciPy库来求解这个线性规划问题。

import pulp

# 创建一个线性规划问题
lp = pulp.LpProblem("Maximize_P", pulp.LpMaximize)

# 定义变量，x1, x2, x3为非负连续变量，P为非负整数变量
x1 = pulp.LpVariable('x1', lowBound=0, cat='Continuous')
x2 = pulp.LpVariable('x2', lowBound=0, cat='Continuous')
x3 = pulp.LpVariable('x3', lowBound=0, cat='Continuous')
P = pulp.LpVariable('P', lowBound=0, cat='Integer')

# 添加目标函数
lp += P

# 添加约束条件
lp += 7*x1 + 6*x2 + 8*x3 - 4*P >= 0
lp += 5*x1 + 9*x2 + 4*x3 - 3*P >= 0
lp += 8*x1 + 5*x2 + 3*x3 <= 300
lp += 6*x1 + 9*x2 + 8*x3 <= 500

# 求解问题
status = lp.solve()

# 输出结果
if status == pulp.LpStatusOptimal:
    print("Optimal integer value of P:", pulp.value(P))
    print("x1:", pulp.value(x1))
    print("x2:", pulp.value(x2))
    print("x3:", pulp.value(x3))
else:
    print("Failed to find an optimal solution.")

运行上述代码将会得到各车间的最优班次安排和最大产品配套数。

Optimal integer value of P: 116.0
x1: 12.571429
x2: 16.190476
x3: 34.857143

练习3

某厂生产甲、乙两种产品，需要A、B两种原料，生产消耗等参数如下表（表中的消耗系数为千克/件）。

产品原料	甲	乙	可用量(千克)	原料成本(元/千克)
A	2	4	160	1.0
B	3	2	180	2.0
销售价（元）	13	16

（1）请构造数学模型使该厂利润最大，并求解。
（2）原料A、B的影子价格各为多少。
（3）工厂可在市场上买到原料A。工厂是否应该购买该原料以扩大生产？

3.1 数学模型

决策变量

• $x$ 为生产甲产品的数量
• $y$ 为生产乙产品的数量

约束条件：

• $2x + 4y \leq 160$（原料A的限制）
• $3x + 2y \leq 180$ （原料B的限制）
• $x \geq，y \geq 0$ (可行性条件)

目标函数（利润最大化）： $$\text{max } Z = 13x + 16y - (2x + 3x\times 2 + 4y + 2y\times 2) = 5x + 8y$$

其中，销售价格减去成本得到的利润为 $5x + 8y$。

3.2 Python求解

from scipy.optimize import linprog

# 目标函数的系数（注意我们求解的是最小化问题，所以取相反数）
c = [-5, -8]

# 不等式左侧系数矩阵
A = [
    [2, 4],
    [3, 2]
]

# 不等式右侧常数
b = [160, 180]

# 求解线性规划问题
result = linprog(c, A_ub=A, b_ub=b, bounds=(0, None), method='highs')

# 输出结果
x, y = result.x
profit = -result.fun

print(f"生产甲产品的数量: {x:.2f}")
print(f"生产乙产品的数量: {y:.2f}")
print(f"最大利润: {profit:.2f}")

生产甲产品的数量: 50.00
生产乙产品的数量: 15.00
最大利润: 370.00

# 提取影子价格
# 使用 result.dual_ub 来获取影子价格
shadow_prices = result.ineqlin.marginals
shadow_price_A = -shadow_prices[0]
shadow_price_B = -shadow_prices[1]

print(f"原料A的影子价格: {shadow_price_A:.2f}")
print(f"原料B的影子价格: {shadow_price_B:.2f}")

# 是否应购买原料A
if shadow_price_A > 1.0:
    print("应购买原料A以扩大生产。")
else:
    print("不应购买原料A。")

原料A的影子价格: 1.75
原料B的影子价格: 0.50
应购买原料A以扩大生产。

练习4

纽约市在平时的一天里警察巡逻的需求如下：

时段	时间	巡逻警察数（最低限）
1	凌晨 12 点 - 凌晨 4 点	6
2	凌晨 4 点 - 早上 8 点	4
3	早上 8 点 - 正午 12 点	14
4	正午 12 点 - 下午 4 点	8
5	下午 4 点 - 晚上 8 点	12
6	晚上 8 点 - 凌晨 12 点	16

每个警察的班次为连续工作8小时，我们需要确定这6个班次各安排多少巡警可以使得满足人员需求的条件下使用的警力最少。

4.1 数学模型

设每个时段起始时的巡逻警数量分别为 $x_1, x_2, x_3, x_4, x_5, x_6$，其中：

• $x_1$ 表示午夜12点到早上8点的巡逻警数量
• $x_2$ 表示凌晨4点到中午12点的巡逻警数量
• $x_3$ 表示早上8点到下午4点的巡逻警数量
• $x_4$ 表示中午12点到晚上8点的巡逻警数量
• $x_5$ 表示下午4点到凌晨12点的巡逻警数量
• $x_6$ 表示晚上8点到凌晨4点的巡逻警数量

根据需求，我们得到以下约束条件：

$x_1 + x_6 \geq 6$
$x_1 + x_2 \geq 4$
$x_2 + x_3 \geq 14$
$x_3 + x_4 \geq 8$
$x_4 + x_5 \geq 12$
$x_5 + x_6 \geq 16$

目标是最小化巡逻警总数，即最小化 $x_1 + x_2 + x_3 + x_4 + x_5 + x_6$

4.2 Python 求解

我们可以使用 scipy.optimize 库中的 linprog 函数来求解这个线性规划问题。

import numpy as np
from scipy.optimize import linprog

# 系数矩阵
A = np.array([
    [-1,  0,  0,  0,  0, -1],
    [-1, -1,  0,  0,  0,  0],
    [ 0, -1, -1,  0,  0,  0],
    [ 0,  0, -1, -1,  0,  0],
    [ 0,  0,  0, -1, -1,  0],
    [ 0,  0,  0,  0, -1, -1]
])

# 需求矩阵
b = np.array([-6, -4, -14, -8, -12, -16])

# 目标函数系数
c = np.array([1, 1, 1, 1, 1, 1])

# 求解线性规划问题
res = linprog(c, A_ub=A, b_ub=b, bounds=(0, None), method='highs')

# 输出结果
print(f"最小巡逻警数: {res.fun}")
print(f"各时段巡逻警数: {res.x}")

最小巡逻警数: 32.0
各时段巡逻警数: [ 2.  2. 12.  0. 12.  4.]

练习5

某工厂在今后4个月内需租用仓库堆放物资。已知各月份所需仓库面积如下表；仓库租用费用随和同期定，期限越长折扣越大，具体数据见下表。租用仓库的合同每月初都可办理，每份合同将具体规定租用面积数和期限。因此该厂可根据需要在任何一个月初办理租用合同，每次可签一份或多份和合同。工厂需决定如何租用可使总的所付租用费用最小。请建立此问题的线性规划模型。

月份	1	2	3	4	合同租用期限（个月）	1	2	3	4
所需仓库面积（百平米）	15	10	20	12	租用费用（元/百平米）	2800	4500	6000	7300

5.1 建模三要素

• 决策变量
  我们定义决策变量 $x_{i,j}$，表示在第$i$个月租用$j$个月的仓库面积数（百平米）。

• 目标函数
  最小化总租用费用。

• 约束条件
  满足各月份的仓库面积需求;租用面积必须非负。

5.2 数学模型

我们将模型用数学公式表达如下：

• 目标函数：

$$\text{Min} \quad Z = \sum_{i=1}^{4} \sum_{j=1}^{4} c_{j} x_{i,j}$$

其中，$c_j$ 是租用j个月的费用。

• 约束条件：

对于每个月的需求约束：

$$\begin{aligned} x_{1,1} + x_{1,2} + x_{1,3} + x_{1,4} & \geq 15 \\ x_{2,1} + x_{1,2} + x_{1,3} + x_{2,4} & \geq 10 \\ x_{3,1} + x_{2,2} + x_{1,3} + x_{3,4} & \geq 20 \\ x_{4,1} + x_{3,2} + x_{2,3} + x_{1,4} & \geq 12 \\\end{aligned}$$

• 非负约束：

$$x_{i,j} \geq 0 \quad \forall i, j$$

5.3 Python求解

import pulp

# 定义模型
model = pulp.LpProblem("Minimize_Warehouse_Rental_Cost", pulp.LpMinimize)

# 定义决策变量
x = pulp.LpVariable.dicts("x", [(i, j) for i in range(1, 5) for j in range(1, 5)], lowBound=0, cat='Continuous')

# 定义租用费用（元/百平米）
cost = {1: 2800, 2: 4500, 3: 6000, 4: 7300}

# 添加目标函数
model += pulp.lpSum(cost[j] * x[i, j] for i in range(1, 5) for j in range(1, 5))

# 添加约束条件
# 月份需求
demand = {1: 15, 2: 10, 3: 20, 4: 12}
for i in range(1, 5):
    model += pulp.lpSum(x[k, j] for j in range(1, 5) for k in range(max(1, i - j + 1), i + 1)) >= demand[i]

# 求解模型
model.solve()

# 输出结果
if model.status == pulp.LpStatusOptimal:
    print("Optimal solution found.")
    for i in range(1, 5):
        for j in range(1, 5):
            print(f"x[{i},{j}] = {pulp.value(x[i, j])}")
    print(f"Total cost: {pulp.value(model.objective)}")
else:
    print("Failed to find an optimal solution.")

Optimal solution found.
x[1,1] = 3.0   x[1,2] = 0.0   x[1,3] = 0.0
x[1,4] = 12.0  x[2,1] = 0.0   x[2,2] = 0.0
x[2,3] = 0.0   x[2,4] = 0.0   x[3,1] = 8.0
x[3,2] = 0.0   x[3,3] = 0.0   x[3,4] = 0.0
x[4,1] = 0.0   x[4,2] = 0.0   x[4,3] = 0.0
x[4,4] = 0.0
Total cost: 118400.0

练习6

某厂生产A、B两种产品，需经过金工和装配两个车间加工，有关数据如下表所示。产品B无论批次大小，每件产品生产成本为400元。产品A的生产成本是分段函数，第1件到70件，每件成本为200元，从第71件开始，每件成本为190元。试建立线性规划整数模型使得该厂生产产品的总利润最大。

工时定额（小时/件）	A	B	总有效工时
金工	4	2	480小时
装配	2	5	500元
售价（元/件）	300	520

6.1 数学模型

• 决策变量：
  $x_A$为生产产品A的数量
  $x_B$为生产产品B的数量

• 目标函数

我们需要最大化总利润。根据题意：产品A的售价为300元/件，前70件的生产成本为200元/件，从第71件起的生产成本为190元/件；产品B的售价为520元/件，生产成本为400元/件。
前70件产品A的利润为： $300x_A - 200x_A = 100x_A$
超过70件产品A的利润为：$300(x_A - 70) - 190(x_A - 70) = 110(x_A - 70)$
产品B的利润为：$520x_B - 400x_B = 120x_B$
所以，目标函数是： $$\text{max } Z = 100 \cdot \min(x_A, 70) + 110 \cdot \max(x_A - 70, 0) + 120x_B$$

• 约束条件：
  金工车间时间限制： $4x_A + 3x_B \leq 480$
  装配车间时间限制： $2x_A + 5x_B \leq 500$
  非负约束： $x_A \geq 0, x_B \geq 0$

6.2 Python求解

import pulp

# 创建一个LP问题实例
prob = pulp.LpProblem("Maximize_Profit", pulp.LpMaximize)

# 定义决策变量
x_A = pulp.LpVariable('x_A', lowBound=0, cat='Integer')
x_B = pulp.LpVariable('x_B', lowBound=0, cat='Integer')

# 分段线性函数的辅助变量
x_A1 = pulp.LpVariable('x_A1', lowBound=0, upBound=70, cat='Integer')
x_A2 = pulp.LpVariable('x_A2', lowBound=0, cat='Integer')

# 定义目标函数
# 前70件产品A的利润
profit_A_1 = 100 * x_A1
# 超过70件产品A的利润
profit_A_2 = 110 * x_A2
# 产品B的利润
profit_B = 120 * x_B

# 总利润
prob += profit_A_1 + profit_A_2 + profit_B

# 添加约束条件
prob += 4 * x_A + 3 * x_B <= 480, "工车间时间限制"
prob += 2 * x_A + 5 * x_B <= 500, "装配车间时间限制"
prob += x_A1 + x_A2 == x_A, "分段线性函数约束"

# 求解问题
status = prob.solve()

# 输出结果
print(f"生产产品A的数量: {pulp.value(x_A):.0f}")
print(f"生产产品B的数量: {pulp.value(x_B):.0f}")
print(f"最大利润: {pulp.value(prob.objective):.2f}")

生产产品A的数量: 64
生产产品B的数量: 74
最大利润: 15920.00

练习7

考虑下面问题，其中$k>0$但其值不能确定。若其最优解为(2,3)，试确定$k$的取值范围。

$$\begin{aligned} &\text{Max} && Z = x_1 + 2x_2 \\ &\text{subject to} && -x_1 + x_2 \leq 2 \\ & && x_2 \leq 3 \\ & && kx_1 + x_2 \leq 2k + 3 \\ & && x_1 \geq 0, \, x_2 \geq 0 \\ \end{aligned}$$

7.1 确定可行域的各个顶点

首先，我们找到可行域的边界，通过求解各个约束的交点来确定可行域的顶点。

• 约束1和约束2的交点：

  $$\begin{cases} -x_1 + x_2 = 2 \\ x_2 = 3 \end{cases}$$

  解得交点为 (1,3)

• 约束1和约束3的交点：

  $$\begin{cases} -x_1 + x_2 = 2 \\ kx_1 + x_2 = 2k + 3 \end{cases}$$

  联立解得：

$$\begin{cases} x_2 = x_1 + 2 \\ x_2 = -kx_1 + 2k + 3 \end{cases}$$

整理得交点为:

$$\left(\frac{2k + 1}{1 + k}, \frac{2k + 1}{1 + k} + 2\right)$$

由于$k>0$，所以$\frac{2k + 1}{1 + k}+2 >3$，所以该点不是可行域的顶点。

• 约束2和约束3的交点：

  $$\begin{cases} x_2 = 3 \\ kx_1 + x_2 = 2k + 3 \end{cases}$$

  解得交点为 (2,3)

• 约束1和坐标轴的交点：

  • 当 $x_1$ = 0:

    $2x_2 = 2$

    交点为 (0,2)

  • 当$x_2$ = 0:

    $x_1 = 0$

    交点为 (0,0)

• 约束3和坐标轴的交点：

  • 当 $x_1$= 0:

    $x_2 \leq 3$

    交点为 (0,3)（但是被 $x_2 \leq 3$ 覆盖）。

  • 当 $x_2 = 0$:

  $$kx_1 \leq 2k + 3 \implies x_1 \leq \frac{2k + 3}{k}$$

  这个点只在 $k \neq 0$ 时有意义。

7.2 计算顶点的目标函数值

顶点 (0, 0) 的目标函数值： Z = 0
顶点 (0, 2) 的目标函数值： Z = 4
顶点 (1, 3) 的目标函数值： Z = 7
顶点 (2, 3) 的目标函数值： Z = 8
顶点 $\left(\frac{2k + 3}{k}, 0\right)$ 的目标函数值： $Z = \frac{2k + 3}{ k}$

7.3 $k$ 的取值范围

满足最优解为 (2, 3)的 $k$ 的取值范围

$$\frac{2k + 3}{k} \le 8$$

得到$k$的取值范围是 $k \geq \frac{1}{2}$。

7.4 Python程序

import numpy as np
import matplotlib.pyplot as plt

# 定义约束条件的直线方程
def constraints(k):
    x1 = np.linspace(0, 5, 400)
    x2_1 = x1 + 2
    x2_2 = 3 * np.ones_like(x1)
    x2_3 = 2*k + 3 - k*x1
    
    return x1, x2_1, x2_2, x2_3

# 计算顶点
def compute_vertices(k):
    vertices = [
        (0, 0),
        (0, 2),
        (1, 3),
        (2, 3),
        ((2*k + 3) / k, 0)
    ]
    return vertices

# 绘图
def plot_feasible_region(k):
    x1, x2_1, x2_2, x2_3 = constraints(k)
    
    plt.figure(figsize=(12, 10))
    
    # 绘制约束条件
    plt.plot(x1, x2_1, label='$-x_1 + x_2 = 2$', color='blue', linewidth=2)
    plt.plot(x1, x2_2, label='$x_2 = 3$', linestyle='--', color='green', linewidth=2)
    plt.plot(x1, x2_3, label=f'${k}x_1 + x_2 = {2*k + 3}$', color='red', linewidth=2)
    
    # 填充可行域
    plt.fill_between(x1, 0, np.minimum(np.minimum(x2_1, x2_2), x2_3), color='gray', alpha=0.3)
    
    # 标注顶点
    vertices = compute_vertices(k)
    for vert in vertices:
        if vert[0] <= 5 and vert[1] <= 5:
            plt.scatter(*vert, color='red', s=100)  # s 参数增加点的大小
            plt.text(vert[0], vert[1], f'({vert[0]:.2f}, {vert[1]:.2f})', fontsize=14, fontweight='bold')
    
    # 标注最优解
    optimal_solution = (2, 3)
    plt.scatter(*optimal_solution, color='orange', s=100)
    plt.text(optimal_solution[0], optimal_solution[1], f'Optimal ({optimal_solution[0]}, {optimal_solution[1]})', fontsize=14, fontweight='bold', color='orange')
    
    plt.xlim(0, 5)
    plt.ylim(0, 5)
    plt.xlabel('$x_1$', fontsize=16, fontweight='bold')
    plt.ylabel('$x_2$', fontsize=16, fontweight='bold')
    plt.legend(fontsize=14)
    plt.title(f'Feasible Region for $k = {k}$', fontsize=18, fontweight='bold')
    plt.grid(True)
    plt.show()

# 计算k=1/4时的可行域
k = 1/4
plot_feasible_region(k)

$k$=1/4	$k$=1

练习8

某公司用两种原油（ A 和 B ）混合加工成两种汽油（甲和乙）。甲、乙两种汽油含原油A的最低比例分别为 50％和 60％，每吨售价分别为 4800 元和 5600 元。该公司现有原油 A 和 B 的库存量分别为 500 吨和 1000 吨，还可以从市场上买到不超过 1500吨的原油 A 。原油 A 的市场价为：购买量不超过 500 吨时的单价为 10000 元/吨；购买量超过 500 吨单不超过 1000 吨时，超过 500 吨的部分 8000 元/吨；购买量超过 1000 吨时，超过 1000 吨的部分 6000 元/吨。该公司应如何安排原油的采购和加工。

8.1 数学模型

• 定义决策变量

$x_A$: 从市场购买的原油A的数量（吨）
$y_1$: 生产的甲种汽油数量（吨）
$y_2$: 生产的乙种汽油数量（吨）
$y_{1A}$: 用于生产甲种汽油的原油A数量（吨）
$y_{1B}$: 用于生产甲种汽油的原油B数量（吨）
$y_{2A}$: 用于生产乙种汽油的原油A数量（吨）
$y_{2B}$: 用于生产乙种汽油的原油B数量（吨）

• 建立约束条件
  库存和采购约束
  原油A的库存和购买总量不能超过可用量： $y_{1A} + y_{2A} \leq 500 + x_A$
  原油B的使用不能超过库存量： $y_{1B} + y_{2B} \leq 1000$
  市场上原油A的采购量限制： $0 \leq x_A \leq 1500$
  加工比例约束
  甲种汽油含原油A的比例至少为50%： $y_{1A} \geq 0.5 y_1$
  乙种汽油含原油A的比例至少为60%： $y_{2A} \geq 0.6 y_2$
  生产总量约束
  甲种汽油总量： $y_1 = y_{1A} + y_{1B}$
  乙种汽油总量： $y_2 = y_{2A} + y_{2B}$

• 目标函数

目标是最大化利润。利润 = 销售收入 - 采购成本。 销售收入为：

$$4800y_1 + 5600y_2$$

采购成本分段计算：

$$\text{采购成本} = \begin{cases} 10000 \times x_A & \text{if } x_A \leq 500 \\ 10000 \times 500 + 8000 \times (x_A - 500) & \text{if } 500 < x_A \leq 1000 \\ 10000 \times 500 + 8000 \times 500 + 6000 \times (x_A - 1000) & \text{if } x_A > 1000 \end{cases}$$

8.2 Python程序

import pulp

# 创建线性规划问题
model = pulp.LpProblem("Oil_Production", pulp.LpMaximize)

# 定义决策变量
x_A = pulp.LpVariable('x_A', lowBound=0, upBound=1500, cat='Continuous')
y_1 = pulp.LpVariable('y_1', lowBound=0, cat='Continuous')
y_2 = pulp.LpVariable('y_2', lowBound=0, cat='Continuous')
y_1A = pulp.LpVariable('y_1A', lowBound=0, cat='Continuous')
y_1B = pulp.LpVariable('y_1B', lowBound=0, cat='Continuous')
y_2A = pulp.LpVariable('y_2A', lowBound=0, cat='Continuous')
y_2B = pulp.LpVariable('y_2B', lowBound=0, cat='Continuous')

# 库存和采购约束
model += y_1A + y_2A <= 500 + x_A
model += y_1B + y_2B <= 1000

# 加工比例约束
model += y_1A >= 0.5 * y_1
model += y_2A >= 0.6 * y_2

# 生产总量约束
model += y_1 == y_1A + y_1B
model += y_2 == y_2A + y_2B

# 定义辅助变量表示分段采购量
x_A1 = pulp.LpVariable('x_A1', lowBound=0, upBound=500, cat='Continuous')
x_A2 = pulp.LpVariable('x_A2', lowBound=0, upBound=500, cat='Continuous')
x_A3 = pulp.LpVariable('x_A3', lowBound=0, upBound=500, cat='Continuous')

# x_A 分段
model += x_A == x_A1 + x_A2 + x_A3
model += x_A2 <= 500
model += x_A3 <= 500

# 目标函数中的采购成本
purchase_cost = (10000 * x_A1 + 
                 8000 * x_A2 + 
                 6000 * x_A3)

# 目标函数：最大化利润
model += 4800 * y_1 + 5600 * y_2 - purchase_cost

# 求解问题
model.solve()

# 输出结果（保留两位小数）
print(f"Optimal value of x_A: {pulp.value(x_A):.2f}")
print(f"Optimal value of y_1: {pulp.value(y_1):.2f}")
print(f"Optimal value of y_2: {pulp.value(y_2):.2f}")
print(f"Optimal value of y_1A: {pulp.value(y_1A):.2f}")
print(f"Optimal value of y_1B: {pulp.value(y_1B):.2f}")
print(f"Optimal value of y_2A: {pulp.value(y_2A):.2f}")
print(f"Optimal value of y_2B: {pulp.value(y_2B):.2f}")
print(f"Maximum Profit: {pulp.value(model.objective):.2f}")

Optimal value of x_A: 1000.00
Optimal value of y_1: 0.00
Optimal value of y_2: 2500.00
Optimal value of y_1A: 0.00
Optimal value of y_1B: 0.00
Optimal value of y_2A: 1500.00
Optimal value of y_2B: 1000.00
Maximum Profit: 7000000.00

练习9

某大学计算机实验室聘用4名大学生(序号1,2,3,4)和两名研究生(序号5,6)值班答疑。一直每人从周一到周五最多可安排的值班时间及每小时值班报酬如下表所示。

学生代号	报酬（元/小时）	周一	周二	周三	周四	周五
1	10.0	6	0	6	0	7
2	10.0	0	6	0	6	0
3	9.9	4	8	3	0	5
4	9.8	5	5	6	0	4
5	10.8	3	0	4	8	0
6	11.3	0	6	0	6	3

该实验室开放时间为上午8:00至晚上10:00，开放时间内须有且仅需一名学生值班。又规定每名大学生每周值班不少于8小时；研究生每周不少于7小时。建立该实验室总支付报酬为最小的数学模型。

• 符号定义

  • $x_{ijk}$：表示学生 $i$ 在周 $j$ 的第 $k$ 天值班的小时数。
  • $i$ 代表学生（1 到 6）
  • $j$ 代表周几（1 到 5，分别代表周一到周五）
  • $k$ 代表天（1 到 5）

• 目标函数
  最小化总报酬，即所有学生在所有工作日的值班时间与对应的报酬之积的总和。

$$\text{Min} \sum_{i=1}^{6} \sum_{j=1}^{5} \sum_{k=1}^{5} c_i \times x_{ijk}$$

其中，$c_i$ 为学生 $i$ 的报酬。

• 约束条件

  • 每天的值班时间不超过表中最大可安排的值班时间：
  $$x_{ijk} \leq \text{max-hours}_{ij}$$
  • 每天必须有且仅有一名学生值班：
  $$\sum_{i=1}^{6} x_{ijk} = 1$$
  • 每名大学生每周值班不少于 8 小时：
  $$\sum_{j=1}^{5} \sum_{k=1}^{5} x_{ijk} \geq 8 \quad \text{for } i \in \{1, 2, 3, 4\}$$
  • 每名研究生每周值班不少于 7 小时：
    $\sum_{j=1}^{5} \sum_{k=1}^{5} x_{ijk} \geq 7 \quad \text{for } i \in {5, 6}$
  • 每名学生每天的值班时间不能超过其最大可安排的值班时间：
  $$x_{ijk} \leq \text{max-hours}_{ij}$$
  • 非负约束：
  $$x_{ijk} \geq 0$$

• 最终模型

$$\text{Min} \sum_{i=1}^{6} \sum_{j=1}^{5} \sum_{k=1}^{5} c_i \times x_{ijk}$$

- 每天的值班时间不超过表中最大可安排的值班时间。
- 每天必须有且仅有一名学生值班。
- 每名大学生每周值班不少于 8 小时。
- 每名研究生每周值班不少于 7 小时。
- 每名学生每天的值班时间不能超过其最大可安排的值班时间。
- 非负约束

通过求解上述整数线性规划模型，可以得到满足所有约束条件且总报酬最小的实验室值班安排方案。

import pulp
import numpy as np

# 数据输入
students = [1, 2, 3, 4, 5, 6]
days = ['周一', '周二', '周三', '周四', '周五']
max_hours = {
    (1, '周一'): 6, (1, '周二'): 0, (1, '周三'): 6, (1, '周四'): 0, (1, '周五'): 7,
    (2, '周一'): 0, (2, '周二'): 6, (2, '周三'): 0, (2, '周四'): 6, (2, '周五'): 0,
    (3, '周一'): 4, (3, '周二'): 8, (3, '周三'): 3, (3, '周四'): 0, (3, '周五'): 5,
    (4, '周一'): 5, (4, '周二'): 5, (4, '周三'): 6, (4, '周四'): 0, (4, '周五'): 4,
    (5, '周一'): 3, (5, '周二'): 0, (5, '周三'): 4, (5, '周四'): 8, (5, '周五'): 0,
    (6, '周一'): 0, (6, '周二'): 6, (6, '周三'): 0, (6, '周四'): 6, (6, '周五'): 3
}
hourly_rate = {1: 10.0, 2: 10.0, 3: 9.9, 4: 9.8, 5: 10.8, 6: 11.3}

# 定义问题
prob = pulp.LpProblem("Minimize_Cost", pulp.LpMinimize)

# 定义决策变量
x = pulp.LpVariable.dicts("x", [(i, d) for i in students for d in days], lowBound=0, cat='Integer')

# 目标函数
prob += pulp.lpSum([hourly_rate[i] * x[i, d] for i in students for d in days])

# 约束条件
# 每天值班时间约束
for d in days:
    prob += pulp.lpSum([x[i, d] for i in students]) == 14

# 每名学生每天值班时间不能超过最大值
for i in students:
    for d in days:
        prob += x[i, d] <= max_hours[i, d]

# 每名大学生每周值班不少于 8 小时
for i in [1, 2, 3, 4]:
    prob += pulp.lpSum([x[i, d] for d in days]) >= 8

# 每名研究生每周值班不少于 7 小时
for i in [5, 6]:
    prob += pulp.lpSum([x[i, d] for d in days]) >= 7

# 求解问题
prob.solve()

# 输出结果矩阵
result_matrix = np.zeros((len(students), len(days)))

for i in students:
    for d_idx, d in enumerate(days):
        result_matrix[i-1][d_idx] = x[i, d].varValue

# 格式化输出结果矩阵，使其保留两位小数
formatted_result_matrix = np.around(result_matrix, 2)

print("排班结果矩阵（每个元素表示学生在对应天的值班小时数）：")
print(formatted_result_matrix)

# 输出总成本
total_cost = round(pulp.value(prob.objective), 2)
print(f"Total Cost: {total_cost}")

排班结果矩阵（每个元素表示学生在对应天的值班小时数）：
[[2. 0. 3. 0. 4.]
 [0. 2. 0. 6. 0.]
 [4. 7. 3. 0. 5.]
 [5. 5. 6. 0. 4.]
 [3. 0. 2. 2. 0.]
 [0. 0. 0. 6. 1.]]
Total Cost: 708.8

练习10

10.1 问题提出

现有数额为 $M$ 的资金用作投资，投资对象为 4 种资产 $s_i (i = 1, ..., 4)$。在这一时期内：
购买 $s_i$ 的收益率为 $r_i$，风险损失率为 $q_i$
购买 $s_i$ 需要付 $p_i$ 的交易费
假定同期银行存款利率是 $r = 5%$，既无交易费又无风险

$$s_i$$	$$r_i %$$	$$q_i %$$	$$p_i %$$
$s_1$	28	2.5	1
$s_2$	21	1.5	2
$s_3$	23	5.5	4.5
$s_4$	25	2.6	6.5

10.2 基本假设与符号规定

除了上述定义的 $s_i, r_i, q_i, p_i (i = 1, ..., 4)$ 外，另令 $s_0$ 指存入银行， $r_0 = 5\%, p_0 = q_0 = 0$
$x_i$ 表示投资 $s_i$ 的资金，
定义 $M = \sum x_i = 1$
定义总体风险为投资的所有 $s_i$ 中风险最大的一个，用 $R$ 表示： $R = \max \{ q_i x_i \}$
$Q$ 表示总体收益： $Q = \sum_{i=0}^{n=4} (r_i - p_i) x_i$
$r_i, q_i, p_i$ 在投资期间内不发生改变，且各个不同的投资项目之间相互独立

10.3 模型建立

该问题是一个多目标规划模型，目标是提高收益并降低风险（ $\max Q, \min R$ ）：

$$\begin{cases} \max;\sum_{i=0}^{n=4}(r_i-p_i)x_i \\ \min\;\max\{q_ix_i\} \end{cases}$$

约束条件为：

$$\begin{cases} \sum_{i=0}^{n=4}(1+p_i)x_i=M=1 \\ x_i\geq0 \end{cases}$$

针对多目标模型（目标1、目标2），我们通常可以通过固定目标1优化目标2，固定目标2优化目标1，或对目标1、2设置相对权重来转换为单目标模型

• 模型一：固定风险，优化收益

固定最大可承受的风险为 $a$，即需满足 $\frac{R}{M} \leq a$：

$$\max \sum_{i=0}^{n=4} (r_i - p_i) x_i \\ \text{s.t.} \quad \frac{(q_i x_i)}{M} \leq a \\ \sum_{i=0}^{n=4} (1 + p_i) x_i = M$$

• 模型二：固定收益，优化风险

固定最小可承受的收益为 $k$，即需满足 $Q \geq k$：

$$\min\;\max\{q_ix_i\}\\ \text{s.t.} \quad \sum_{i=0}^{n=4}(r_i-p_i)x_i \geq k \\ $$ \sum_{i=0}^{n=4}(1+p_i)x_i = M$$

• 模型三：设置收益与风险的相对权重

对风险、收益分别赋予权重 $s, (1 - s)$：

$$\min\;s(\max\{q_ix_i\})-(1-s)\sum_{i=0}^{n=4}(r_i-p_i)x_i \\ s.t. \sum_{i=0}^{n=4}(1+p_i)x_i = M$$

10.4 模型求解（以模型一为例）

$$\begin{aligned} &\text{min} \quad f = (-0.05, -0.27, -0.19, -0.185, -0.185) \begin{pmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \\ x_4 \end{pmatrix}^T \\ &\text{s.t.} \quad \begin{cases} x_0 + 1.01x_1 + 1.02x_2 + 1.045x_3 + 1.065x_4 = 1 \\ 0.025x_1 \leq a \\ 0.015x_2 \leq a \\ 0.055x_3 \leq a \\ 0.026x_4 \leq a \\ x_i \geq 0 \quad (i = 0, 1, \ldots, 4) \end{cases} \end{aligned}$$

由于 $a$ 是任意给定的风险度，因此不妨从 $a = 0, \Delta = 0.001$ 开始，循环探求不同 $a$ 对于最大收益的影响：

import numpy as np
from scipy.optimize import linprog
import matplotlib.pyplot as plt

a = 0
a_values = []
Q_values = []

while a < 0.05:
    c = [-0.05, -0.27, -0.19, -0.185, -0.185]
    A = np.zeros((4, 5))
    A[:, 1:] = np.diag([0.025, 0.015, 0.055, 0.026])
    b = a * np.ones(4)
    Aeq = np.array([[1, 1.01, 1.02, 1.045, 1.065]])
    beq = 1
    bounds = [(0, None)] * 5

    res = linprog(c, A_ub=A, b_ub=b, A_eq=Aeq, b_eq=beq, bounds=bounds, method='highs')
    
    if res.success:
        Q = -res.fun
        a_values.append(a)
        Q_values.append(Q)
    
    a += 0.001

plt.plot(a_values, Q_values, '*k')
plt.xlabel('a')
plt.ylabel('Q')
plt.show()

由图可知最大收益与风险承受能力的关系，特别的：
$a_1 = 0.006$是第一个拐点，其左侧最大收益随可承受的风险增加而快速增长，右侧的增长则较为缓慢。
$a_2 = 0.025$ 是第二个拐点，其右侧最大收益不再增长。
因此，可以判断$a = 0.006$是一个比较好的风险承受能力，当选择$a = 0.006$ 时，求解线性规划可得：

$$Q = 0.2019, x_0 = 0, x_1 = 0.24, x_2 = 0.4, x_3 = 0.1091, x_4 = 0.2212$$