运筹学练习Python精解——运输问题

练习1

某公司在如下3个地方生产商品，并运送到另外7个地点进行销售，请问该公司如何配送运输的成本最低？

产地\销地	FRA	DET	LAN	WIN	STL	FRE	LAF	供应量
GARY	inf	14	11	inf	16	inf	8	1400
CLEV	27	8	12	9	26	inf	17	2600
PITT	24	inf	inf	13	28	99	inf	2900
需求量	900	1200	600	400	1700	1100	1000	6900

import numpy as np
from scipy.optimize import linprog

# 定义成本矩阵，inf表示无限大，表示不可运输
cost_matrix = [
    [np.inf, 14, 11, np.inf, 16, np.inf, 8],
    [27, 8, 12, 9, 26, np.inf, 17],
    [24, np.inf, np.inf, 13, 28, 99, np.inf]
]

# 将inf替换为一个足够大的数字，以避免数值问题
large_number = 1e6
cost_matrix = np.where(np.isinf(cost_matrix), large_number, cost_matrix)

# 将成本矩阵展平成一维数组
c = cost_matrix.flatten()

# 定义供应量和需求量
supply = [1400, 2600, 2900]
demand = [900, 1200, 600, 400, 1700, 1100, 1000]

# 构建约束矩阵
num_supply = len(supply)
num_demand = len(demand)

A_eq = []
b_eq = []

# 供应量约束
for i in range(num_supply):
    constraint = [0] * (num_supply * num_demand)
    for j in range(num_demand):
        constraint[i * num_demand + j] = 1
    A_eq.append(constraint)
    b_eq.append(supply[i])

# 需求量约束
for j in range(num_demand):
    constraint = [0] * (num_supply * num_demand)
    for i in range(num_supply):
        constraint[i * num_demand + j] = 1
    A_eq.append(constraint)
    b_eq.append(demand[j])

A_eq = np.array(A_eq)
b_eq = np.array(b_eq)

# 定义变量边界
x_bounds = [(0, None)] * (num_supply * num_demand)

# 求解线性规划问题
result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=x_bounds, method='highs')

# 输出结果
if result.success:
    print("Optimal value:", result.fun)
    x_values = result.x.reshape((num_supply, num_demand))
    print("x values:")
    print(x_values)
else:
    print("No solution found")

Optimal value: 200500.0
x values:
[[   0.    0.    0.    0.  800.    0.  600.]
 [   0. 1200.  600.  400.    0.    0.  400.]
 [ 900.    0.    0.    0.  900. 1100.    0.]]

练习2

如下表的运输问题中总需要量超过总供应量(方框中的数字是单位运费)。假定对销地$B_1$、$B_2$和$B_3$未满足需要量的单位罚款成本是5、3和2，试建立该问题的数学模型，并探讨能否将其转变为产销平衡运输问题。

产地\销地	B1	B2	B3	供 应 量
A1	5	1	7	10
A2	6	4	6	80
A3	3	2	5	15
需求量	75	20	50	145\105

2.1 数学模型

设 $x_{ij}$ 为从产地 $A_i$ 运输到销地 $B_jBj​ 的货物量。

• 目标函数

最小化总运输成本和罚款成本： $$\min \sum_{i=1}^{3} \sum_{j=1}^{3} c_{ij} x_{ij} + \sum_{j=1}^{3} p_j y_j$$ 其中，
$c_{ij}$​ 是从 $A_i$​ 到 $B_j$ 的单位运费；$p_j$​ 是销地 $B_j$​ 的单位罚款成本；$y_j$​ 是销地 $B_j$​ 的未满足需要量。

• 约束条件
  供应量约束： $\sum\_{j=1}^{3} x_{ij} \leq s_i \quad \forall i = 1, 2, 3$
  需求量约束（包括未满足的部分）： $\sum_{i=1}^{3} x_{ij} + y_j = d_j \quad \forall j = 1, 2, 3$
  非负性约束： $x_{ij} \geq 0 \quad \forall i = 1, 2, 3；y_j \geq 0 \quad \forall j = 1, 2, 3$
• 参数

单位运价$c_{ij}:$

$$\begin{matrix} & B_1 & B_2 & B_3 \\ A_1 & 5 & 1 & 7 \\ A_2 & 6 & 4 & 6 \\ A_3 & 3 & 2 & 5 \end{matrix}$$

供应量 $s_i$​：
$$s = [10, 80, 15]$$
需求量$d_j$：

$$d = [75, 20, 50]$$

单位罚款成本$p_j$​：
$$p = [5, 3, 2]$$

2.2 转化为平衡型运输问题

为将其转化为平衡型运输问题，定义一个虚拟产地 $A_4$​，其供应量为未满足的总需求量，即 $\sum d_j - \sum s_i = 145 - 105 = 40$。
新增的运费矩阵行对应的是虚拟产地 $A_4$​​，其到各个销地 $B_j$​ 的运费为对应的罚款成本。

新的运费矩阵为：

$$\begin{matrix} & B_1 & B_2 & B_3 \\ A_1 & 5 & 1 & 7 \\ A_2 & 6 & 4 & 6 \\ A_3 & 3 & 2 & 5 \\ A_4 & 5 & 3 & 2 \\ \end{matrix}$$

新的供应量向量$s$ 为：

$$s = [10, 80, 15, 40]$$

需求量向量保持不变：

$$d = [75, 20, 50]$$

2.3 Python计算程序

import numpy as np
from scipy.optimize import linprog

# 定义参数
c = [5, 1, 7, 6, 4, 6, 3, 2, 5, 5, 3, 2]
s = [10, 80, 15, 40]
d = [75, 20, 50]

# 构建目标函数系数向量
c = np.array(c)

# 构建约束矩阵
A_eq = [
    [1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1],
    [1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0],
    [0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0],
    [0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1]
]

A_eq = np.array(A_eq)
b_eq = np.array(s + d)

# 定义边界
x_bounds = [(0, None) for _ in range(len(c))]

# 求解线性规划问题
result = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=x_bounds, method='highs')

# 输出结果
if result.success:
    print("Optimal value:", result.fun)
    print("x values:", result.x.reshape(4, 3))
else:
    print("No solution found")

Optimal value: 595.0
x values: 
[[ 0. 10.  0.]
 [60. 10. 10.]
 [15.  0.  0.]
 [ 0.  0. 40.]]

练习3

在下表不平衡运输问题中(方框中的数字是单位运费)，若产地 有单位物质未运出，就要发生存储成本。假定在产地$A1$、$A2$和$A3$的单位存储成本是5、4和3，又假定产地 的供应量必须全部运出，试建立该问题的数学模型，并探讨能否将其转变为产销平衡运输问题。

产地\销地	B1	B2	B3	供 应 量
A1	1	2	1	20
A2	0	4	5	40
A3	2	3	3	30
需求量	30	20	20	70\90

3.1 数学模型

供应地$A_1, A_2, A_3$ 的供应量分别是 20, 40, 30；需求地$B_1, B_2, B_3$的需求量分别是 30, 20, 20；存储成本分别是 ￥s_1 = 5, s_2 = 4, s_3 = 3 $。运输成本矩阵：

$$C = \begin{bmatrix} 1 & 2 & 1 \\ 0 & 4 & 5 \\ 2 & 3 & 3 \end{bmatrix}$$

3.2 转化为平衡型运输问题

• 供应和需求总量
  供应总量 20 + 40 + 30 = 90 万吨；需求总量 30 + 20 + 20 = 70 万吨

• 引入假设需求节点 $B_4$
  假设需求节点的需求量是 20，对应的运输成本为各供应地的存储成本

• 更新后的运输成本矩阵

$$\begin{bmatrix} 1 & 2 & 1 & 5 \\ 0 & 4 & 5 & 4 \\ 2 & 3 & 3 & 3 \end{bmatrix}$$

• 新的需求量

$$\begin{bmatrix} 30, 20, 20, 20 \end{bmatrix}$$

使用上述新的运输成本矩阵和需求量构建平衡型运输问题。

3.3 Python求解

import numpy as np
from scipy.optimize import linprog

# 供应量
supply = np.array([20, 40, 30])

# 需求量
demand = np.array([30, 20, 20, 20])

# 运输成本矩阵
cost = np.array([
    [1, 2, 1, 5],
    [0, 4, 5, 4],
    [2, 3, 3, 3]
])

# 将成本矩阵展平
c = cost.flatten()

# 创建约束矩阵
A_eq = []
for i in range(len(supply)):
    row = [0] * len(c)
    for j in range(len(demand)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

for j in range(len(demand)):
    row = [0] * len(c)
    for i in range(len(supply)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

A_eq = np.array(A_eq)

# 供应和需求量
b_eq = np.concatenate([supply, demand])

# 求解线性规划问题
result = linprog(c, A_eq=A_eq, b_eq=b_eq, method='highs')

# 打印结果
if result.success:
    print("最优解:")
    x = result.x.reshape(cost.shape)
    for i in range(len(supply)):
        for j in range(len(demand)):
            print(f"从 A{i+1} 到 B{j+1} 运输量: {x[i, j]}")
    print(f"最小总运输成本: {result.fun}")
else:
    print("没有找到最优解")

练习4

某厂月底安排某一产品在下月四周的生产计划。估计每件产品在第一周与第二周的生产成本为150元，后两周的生产成本为170元，各周产品需求量分别为700件、800件、1000件和1200件。工厂每周至多生产产品900件，在第二周、第三周可加班生产。加班生产时每周可增产三百件，但生产成本每件需增加30元。过剩产品的储存费为每周15元。试安排生产，使总成本最小，建立运输模型。

需求总量：700 + 800 + 1000 + 1200 = 3700
最高生产量：900 + 900 + 300 + 900 + 300 + 900 = 4200
需求缺少：4200 - 3700 = 500

周别	B1	B2	B3	B4	供应量
A1	150	165	180	195	900
A2	M	150	165	180	900
A3	M	180	195	210	300
A4	M	M	170	185	900
A5	M	M	200	215	300
A6	M	M	M	170	900
需求量	700	800	1000	1200	3700\4200

4.1 转化为平衡型运输问题

周别	B1	B2	B3	B4	B5	供应量
A1	150	165	180	195	0	900
A2	M	150	165	180	0	900
A3	M	180	195	210	0	300
A4	M	M	170	185	0	900
A5	M	M	200	215	0	300
A6	M	M	M	170	0	900
需求量	700	800	1000	1200	500	4200\4200

4.2 Python程序

import numpy as np
from scipy.optimize import linprog
import pandas as pd

# 定义非常大的数字来表示不可行路径
M = 9999

# 供应量
supply = np.array([900, 900, 300, 900, 300, 900])

# 需求量
demand = np.array([700, 800, 1000, 1200, 500])

# 运输成本矩阵
cost = np.array([
    [150, 165, 180, 195, 0],
    [M, 150, 165, 180, 0],
    [M, 180, 195, 210, 0],
    [M, M, 170, 185, 0],
    [M, M, 200, 215, 0],
    [M, M, M, 170, 0]
])

# 将成本矩阵展平
c = cost.flatten()

# 创建约束矩阵
A_eq = []
for i in range(len(supply)):
    row = [0] * len(c)
    for j in range(len(demand)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

for j in range(len(demand)):
    row = [0] * len(c)
    for i in range(len(supply)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

A_eq = np.array(A_eq)

# 供应和需求量
b_eq = np.concatenate([supply, demand])

# 求解线性规划问题
result = linprog(c, A_eq=A_eq, b_eq=b_eq, method='highs')

# 打印结果
if result.success:
    # 重塑结果
    x = result.x.reshape(cost.shape)
    
    # 创建一个DataFrame来显示结果
    df = pd.DataFrame(x, columns=['B1', 'B2', 'B3', 'B4', 'B5'], index=['A1', 'A2', 'A3', 'A4', 'A5', 'A6'])
    df.loc['需求量'] = demand
    df['供应量'] = np.append(supply, [0])

    print("最优解:")
    print(df)
    print(f"\n最小总运输成本: {result.fun}")
else:
    print("没有找到最优解")

最优解:
        B1     B2      B3      B4     B5  供应量
A1   700.0    0.0   200.0     0.0    0.0  900
A2     0.0  700.0   200.0     0.0    0.0  900
A3     0.0  100.0     0.0     0.0  200.0  300
A4     0.0    0.0   600.0   300.0    0.0  900
A5     0.0    0.0     0.0     0.0  300.0  300
A6     0.0    0.0     0.0   900.0    0.0  900
需求量  700.0  800.0  1000.0  1200.0  500.0  

最小总运输成本: 607500.0

练习5

有甲、乙、丙3个城市，每年分别需要煤炭320、250、350万吨，由A、B两个煤矿负责供应。已知煤矿年产量A为400万吨；B为450万吨，从两煤矿至各城市煤炭运价(元/吨)见下表。由于需求大于产量，经协商平均，甲城市必要时可少供0～30万吨，乙城市需求量须全部满足，丙城市需求量不少于270万吨。试求将甲、乙两矿煤炭全部分配出去，满足上述条件又使总运费为最低的调运方案。

产地\销地	甲	乙	丙	供 应 量
A	15	18	22	400
B	21	25	16	450
需求量	320	250	350

5.1 转化为平衡型运输问题

产地\销地	甲	甲'	乙	丙	丙'	供 应 量
A	15	15	18	22	22	400
B	21	21	25	16	16	450
需求量	290	30	250	270	80	920\850

产地\销地	甲	甲'	乙	丙	丙'	供 应 量
A	15	15	18	22	22	400
B	21	21	25	16	16	450
C	0	0	0	0	0	70
需求量	290	30	250	270	80	920\920

5.2 Python程序

import numpy as np
from scipy.optimize import linprog
import pandas as pd

# 供应量
supply = np.array([400, 450, 70])

# 需求量
demand = np.array([290, 30, 250, 270, 80])

# 运输成本矩阵
cost = np.array([
    [15, 15, 18, 22, 22],
    [21, 21, 25, 16, 16],
    [0, 0, 0, 0, 0]
])

# 将成本矩阵展平
c = cost.flatten()

# 创建约束矩阵
A_eq = []
for i in range(len(supply)):
    row = [0] * len(c)
    for j in range(len(demand)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

for j in range(len(demand)):
    row = [0] * len(c)
    for i in range(len(supply)):
        row[i * len(demand) + j] = 1
    A_eq.append(row)

A_eq = np.array(A_eq)

# 供应和需求量
b_eq = np.concatenate([supply, demand])

# 求解线性规划问题
result = linprog(c, A_eq=A_eq, b_eq=b_eq, method='highs')

# 打印结果
if result.success:
    # 重塑结果
    x = result.x.reshape(cost.shape)
    
    # 创建一个DataFrame来显示结果
    df = pd.DataFrame(x, columns=['甲', '甲\'', '乙', '丙', '丙\''], index=['A', 'B', 'C'])
    df.loc['需求量'] = demand
    df['供应量'] = np.append(supply, [0])

    print("最优解:")
    print(df)
    print(f"\n最小总运输成本: {result.fun}")
else:
    print("没有找到最优解")

最优解:
         甲    甲'      乙      丙    丙'  供应量
A    190.0  30.0  180.0    0.0   0.0  400
B    100.0   0.0    0.0  270.0  80.0  450
C      0.0   0.0   70.0    0.0   0.0   70
需求量  290.0  30.0  250.0  270.0  80.0    0

最小总运输成本: 14240.0

练习6

某餐馆承办宴会, 每晚连续举行, 共举行五次。宴会上需用特殊的餐巾, 根据参加的人数, 预计每晚的需要量为：第一天 1000 条，第二天700条，第三天 800 条，第四工 1200 条，第五工 1500 条，五天之后，所有的餐巾作废。宴会中用过的餐巾经过洗涤处理后可以重复使用，这样可以降低使用成本。已知每条新餐巾需要1元的费用 送洗时可选择两种方式:快洗仅需要一天时间, 每条洗涤费用为 0.2 元, 慢洗需要两天时间, 每条洗涤费用 0.1 元。问：如何安排,可使总费用最低?

6.1 数学模型

设 $x_j$ ——第 $j$ 天使用新毛巾的数量; $y_{i j}$——第 $i$ 天送第 $j$ 天使用快洗餐巾的数量; $\mathrm{z}_{\mathrm{ij}}$ 一第 $\mathrm{i}$ 天送第 $\mathrm{j}$ 天使用慢洗餐巾的数量;

$$\begin{aligned} & \operatorname{Min} \quad \mathrm{z}=\sum \mathrm{x}_{\mathrm{j}}+\sum \sum 0.2 \mathrm{y}_{\mathrm{ij}}+\sum \sum 0.1 \mathrm{z}_{\mathrm{ij}} \\ & \text { 第一天: } \mathrm{x}_1=1000 \\ & \text { 第二天: } x_2+y_{12}=700 \\ & \text { 第三天: } x_3+z_{13}+y_{23}=800 \\ & \text { 第四天: } x_4+z_{14}+z_{24}+y_{34}=1200 \\ & \text { 新购餐巾: } x_1+x_2+x_3+x_4+x_5 \leq 5200 \\ & \text { 第一天送洗: } \mathrm{y}_{12}+\mathrm{z}_{13}+\mathrm{z}_{14}+\mathrm{z}_{15} \leq 1000 \\ & \text { 第二天送洗: } \mathrm{y}_{23}+\mathrm{z}_{24}+\mathrm{z}_{25} \leq 700 \\ & \text { 第三天送洗: } \mathrm{y}_{34}+\mathrm{z}_{35} \leq 800 \\ & \text { 第五天: } x_5+z_{15}+z_{25}+z_{35}+y_{45}=1500 \\ & \text { 第四天送洗: } \mathrm{y}_{45} \leq 1200 \\ & x_j \geq 0, \quad y_{i j} \geq 0, \quad z_{i j} \geq 0, \quad (i=1, \cdots, 4 ; \quad j=1, \cdots, 5) \\ \end{aligned}$$

6.2 转化为平衡型运输问题

供应需求	I	II	III	IV	V	VI	产量
新购	1	1	1	1	1	0	5200
第一天	M	0.2	0.1	0.1	0.1	0	1000
第二天	M	M	0.2	0.1	0.1	0	700
第三天	M	M	M	0.2	0.1	0	800
第四天	M	M	M	0.2	0.1	0	1200
销量	1000	700	800	1200	1500	3700

6.3 Python程序

import numpy as np
from scipy.optimize import linprog
import pandas as pd

# 定义问题的参数
cost = [1, 1, 1, 1, 1, 0, 0.2, 0.2, 0.2, 0.2, 0.1, 0.1, 0.1, 0.1]  # 目标函数系数

# 约束矩阵
A_eq = [
    [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # 第一天
    [0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],  # 第二天
    [0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0],  # 第三天
    [0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0],  # 第四天
    [0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1]   # 第五天
]

b_eq = [1000, 700, 800, 1200, 1500]

A_ub = [
    [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # 新购餐巾
    [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0],  # 第一天送洗
    [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0],  # 第二天送洗
    [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0],  # 第三天送洗
    [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0]   # 第四天送洗
]

b_ub = [5200, 1000, 700, 800, 1200]

# 解决线性规划问题
res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, A_ub=A_ub, b_ub=b_ub, method='highs')

# 提取结果
solution = res.x
optimal_value = res.fun

# 打印最优值
print("最优值:", optimal_value)

# 打印最优解
print("最优解:", solution)

# 创建结果表格
columns = ['I', 'II', 'III', 'IV', 'V', 'VI', '产量']
rows = ['新购', '第一天', '第二天', '第三天', '第四天', '销量']

data = [
    [solution[0], solution[1], solution[2], solution[3], solution[4], solution[5], 5200],
    ['M', solution[6], solution[7], solution[8], solution[9], 0, 1000],
    ['M', 'M', solution[10], solution[11], solution[12], 0, 700],
    ['M', 'M', 'M', solution[13], 0, 800],
    ['M', 'M', 'M', 'M', 0, 1200],
    [1000, 700, 800, 1200, 1500, 3700, '']
]

df = pd.DataFrame(data, columns=columns, index=rows)

# 打印结果表格
print("\n最优解表格:")
print(df)

最优值: 2210.0
最优解表格:
          I     II    III      IV       V      VI    产量
新购   1000.0    0.0  800.0    -0.0     0.0     0.0  5200
第一天       M  700.0   -0.0     0.0     0.0     0.0  1000
第二天       M      M    0.0  1200.0  1500.0     0.0   700
第三天       M      M      M     0.0     0.0   800.0  None
第四天       M      M      M       M     0.0  1200.0  None
销量     1000    700    800    1200  1500.0  3700.0  

练习7

某工厂按合同规定必须于当年的每个季度末分别提供 $10、15、25、20$台同一规格的柴油机。已知该厂的生产能力及生产每台柴油机的成本如表示。又如果生产出来的柴油机当季不交货, 每台每积压一个季度需要存储维护费用 0.15 万元。要求在完成合同的情况下, 做出使全年生产费用最小的决策。

季度	生产能力 (台)	单位成本 (万元/台)
I	25	10.8
II	35	11.1
III	30	11.0
IV	10	11.3

7.1 数学模型

设 $x_{i j}$ 第 $i$ 季度生产, 用于第 $\mathrm{j}$ 季度交货的数量。
$\quad \min \quad z=\sum_{i=1}^4 \sum_{ j=1}^4 c_{i j}$
供应:
$\quad \begin{array}{r}\mathrm{x}_{11}+\mathrm{x}_{12}+\mathrm{x}_{13}+\mathrm{x}_{14} \leq 25 \\ \mathrm{x}_{22}+\mathrm{x}_{23}+\mathrm{x}_{24} \leq 35\\ \mathrm{x}_{33}+\mathrm{x}_{34} \leq 30\\ \mathrm{x}_{44} \leq 10 \end{array}$
需求:
$\begin{cases}x_{11} & =10 \\ x_{12}+x_{22} & =15 \\ x_{13}+x_{23}+x_{33} & =25 \\ x_{14}+x_{24}+x_{34}+x_{44} & =20\end{cases}$

单位费用表(万元)

供应\需求	I	II	III	IV	供应量
I	10.8	10.95	11.10	11.25	25
II	M	11.10	11.25	11.40	35
III	M	M	11.00	11.15	30
IV	M	M	M	11.30	10
需求量	10	15	25	20

7.2 转化为平衡型运输问题

供应\需求	I	II	III	IV	虚拟需求点	供应量
I	10.8	10.95	11.10	11.25	0	25
II	M	11.10	11.25	11.40	0	35
III	M	M	11.00	11.15	0	30
IV	M	M	M	11.30	0	10
需求量	10	15	25	20	30

7.3 Python程序

from scipy.optimize import linprog
import numpy as np

# 成本矩阵
costs = np.array([
    [10.8, 10.95, 11.10, 11.25, 0],
    [1e6, 11.10, 11.25, 11.40, 0],
    [1e6, 1e6, 11.00, 11.15, 0],
    [1e6, 1e6, 1e6, 11.30, 0]
])

# 供应量
supply = np.array([25, 35, 30, 10])

# 需求量
demand = np.array([10, 15, 25, 20, 30])

# 创建目标函数的系数列表
c = costs.flatten()

# 创建左边的不等式约束矩阵和右边的约束向量
A_eq = []
b_eq = []

# 供应量约束
for i in range(len(supply)):
    row = [0] * costs.size
    row[i*costs.shape[1]:(i+1)*costs.shape[1]] = [1] * costs.shape[1]
    A_eq.append(row)
    b_eq.append(supply[i])

# 需求量约束
for j in range(len(demand)):
    row = [0] * costs.size
    row[j::costs.shape[1]] = [1] * costs.shape[0]
    A_eq.append(row)
    b_eq.append(demand[j])

# 转换为numpy数组
A_eq = np.array(A_eq)
b_eq = np.array(b_eq)

# 使用linprog求解
result = linprog(c, A_eq=A_eq, b_eq=b_eq, method='highs')

# 提取结果
x = result.x.reshape(costs.shape)
total_cost = result.fun

# 输出最优解和最优值
print("最优解：")
print(x)
print("\n最优值：")
print(total_cost)

最优解：
[[10. 10.  0.  5.  0.]
 [ 0.  5.  0.  0. 30.]
 [ 0.  0. 25.  5.  0.]
 [ 0.  0.  0. 10.  0.]]

最优值：
773.0

练习8

布里斯班机场扩建运输方案

• 问题描述
  布里斯班机场扩建需要将大量泥土从4个供应地运往7个需求地。下表展示了不同供应地和需求地间的距离（单位：百米）以及每个供应地的可供泥土量（单位：立方米）。

供应地/需求地	延展地	干塘	道路	停车场	消防站	工业园区	周边道路	可得数量
Apron	26	12	10	18	11	8	20	660
Term	28	14	12	20	13	10	22	301
Cargo	20	26	20	2	6	22	18	271
Access	26	10	4	16	24	14	21	99
需求量	247	394	265	105	90	85	145	1331\ 1331

• 分析
  属于典型运输问题，去连续变量，非负，LP问题，模型构造基本没有难度，但为了完整起见，写得详细一点：
  • I：供应地编号集合，$i = 1,2,3,4$
  • J：需求地编号集合$j = 1,2,\ldots,7$

在MATLAB和Python中分别建立模型，因笔者对两种语言特征的编写倾向，前者以编号形式，后者个人的见解是采用字典更灵活。

• 决策变量
  $x_{i,j}$：从供应地$i$运至需求地$j$的泥土立方米数$i = 1,2,\ldots,4; j = 1,2,\ldots,7$

• 约束条件

  • 任意4个供应地之一，向7个位置供应泥土的总量不超过自身产量；
  • 任意7个位置自4个供应地接收的泥土总和不能小于该位置的项目需求。

• 目标函数
  最小化总的运距与运量二者的乘积。

• 已知参数

  • $P_{i,j}$：从4个供应地到7个需求地的距离（dist）
  • $S_i$：供应地的最多泥土供应量（supply）
  • $D_j$：需求地所需泥土的需求数量（demand）

• 完整模型

$$\begin{aligned} \min \quad & \sum_{i=1}^4 \sum_{j=1}^7 x_{i,j} P_{i,j} \\ s.t.: \quad & \left\{ \begin{aligned} & \sum_{j=1}^7 x_{i,j} \le S_i, \quad \forall i \\ & \sum_{i=1}^4 x_{i,j} \ge D_j, \quad \forall j \\ & x_{i,j} \ge 0, \quad \forall(i = 1,2,3,4; j = 1,2,\ldots,7) \end{aligned} \right. \end{aligned}$$

import cvxpy as cp
import numpy as np
import pandas as pd

# Data Initialization
supply_nodes = ["Apron", "Term", "Cargo", "Access"]
demand_nodes = ["延展地", "干塘", "道路", "停车场", "消防站", "工业园区", "周边道路"]

# Cost matrix (transportation cost from each supply node to each demand node)
cost_matrix = np.array([
    [26, 12, 10, 18, 11, 8, 20],
    [28, 14, 12, 20, 13, 10, 22],
    [20, 26, 20, 2, 6, 22, 18],
    [26, 10, 4, 16, 24, 14, 21]
])

# Supply and demand quantities
supply_quantities = np.array([660, 301, 271, 99])
demand_quantities = np.array([247, 394, 265, 105, 90, 85, 145])  # Adjust last value to match supply

# Initialize optimization variables
x = cp.Variable((len(supply_nodes), len(demand_nodes)), nonneg=True)

# Objective function: Minimize transportation cost
objective = cp.Minimize(cp.sum(cp.multiply(cost_matrix, x)))

# Constraints
constraints = [
    cp.sum(x, axis=1) <= supply_quantities,  # Supply constraints
    cp.sum(x, axis=0) == demand_quantities   # Demand constraints
]

# Define the problem
prob = cp.Problem(objective, constraints)

# Solve the problem
result = prob.solve()

if result is not None and x.value is not None:
    # Rounding solution to nearest integers
    solution = np.round(x.value).astype(int)  # Round to the nearest integer and convert to int

    # Create a DataFrame for the solution
    solution_df = pd.DataFrame(solution, index=supply_nodes, columns=demand_nodes)

    print("Optimal Value (Total Transportation Cost):", result)
    print("Solution (Transportation Plan):")
    print(solution_df)

    # Add the supply quantities to the solution DataFrame
    solution_df['可得数量'] = supply_quantities

    # Add the demand quantities as a new row using pd.concat
    demand_row = pd.DataFrame(demand_quantities.reshape(1, -1), index=['需求量'], columns=demand_nodes)
    solution_df = pd.concat([solution_df, demand_row])

    # Output the result in a similar table structure
    print("\nResult Table:")
    print(solution_df)
else:
    print("No feasible solution found.")
#总距离：17592

Solution (Transportation Plan):
        延展地   干塘   道路  停车场  消防站  工业园区  周边道路
Apron    49  283  106    0   63    59   101
Term     32  111   60    0   27    26    44
Cargo   166    0    0  105    0     0     0
Access    0    0   99    0    0     0     0

运输方案写成了从4个供应地(行)向7个需求地(列)的运量形式，例如从供应地1(Apron)运往需求地3(道路)的运输量为 $x_{1,3}=106$。