Python数学建模算法与应用（司守奎）–第4章随笔

第 4 章  线性规划和整数规划模型

4.1 线性规划（Linear Programming，LP）

4.1.1 线性规划模型

线性规划模型的一般形式

，

或简写为

其向量表示形式为

其矩阵表示形式为

其中， 为目标函数的系数向量，又称为价值向量；

 为决策向量；

 为约束方程组的系数矩阵；

  为  的列向量，又称约束方程组的系数向量；

 为约束方程组的常数向量。

4.1.2 模型求解及应用

可以使用 Python 的 cvxpy 库，用于求解凸优化问题。http://www.vcxpy.org/

例 4.2 求解线性规划模型

# 程序文件 ex4_2.py
import cvxpy as cp
from numpy import array
c = array([70, 50, 60])                      # 定义目标向量
a = array([[2, 4, 3],[3, 1, 5], [7, 3, 5]])  # 定义约束矩阵
b = array([150, 160, 200])                   # 定义约束条件的右边向量
x = cp.Variable(3, pos=True)                 # 定义 3 个决策变量
obj = cp.Maximize(c@x)                       # 构造目标函数
cons = [a@x<=b]                              # 构造约束条件
prob = cp.Problem(obj, cons)
prob.solve(solver='GLPK_MI')                 # 求解问题
print('最优解为：', x.value)
print('最优值为：', prob.value)

例 4.3 某部分今后 5 年考虑以下投资项目，已知：

项目A，从第一年到第四年每年年初需投资，并于次年末收回本利 115%。

项目B，从第三年初需投资，到第五年末能回收本利125%，最大投资额不超过4万元。

项目C，第二年初需投资，到第五年末能收回本利140%，最大投资额不超过3万元。

项目D，五年内每年初购买，于当年末归还，利息收益6%。

部门现有资金10万元，试问如何确定这些项目每年投资额，使第五年末本利总额最大？

建立线性规划模型

# 程序文件 ex4_3.py
import cvxpy as cp
x = cp.Variable((5, 4), pos=True)
obj = cp.Maximize(1.15*x[3, 0]+1.40*x[1, 2]+1.25*x[2, 1]+1.06*x[4, 3])
cons = [x[0, 0]+x[0, 3] == 100000,
       x[1, 0]+x[1, 2]+x[1, 3] == 1.06*x[0, 3],
       x[2, 0]+x[2, 1]+x[2, 3] == 1.15*x[0, 0]+1.06*x[1, 3],
       x[3, 0]+x[3, 3] == 1.15*x[1, 0]+1.06*x[2, 3],
       x[4, 3] == 1.15*x[2, 0]+1.06*x[3, 3],
       x[2, 1]<=40000,x[1,2]<=30000]
prob = cp.Problem(obj, cons)
prob.solve(solver='GLPK_MI')
print('最优解为：', x.value)
print('最优值为：', prob.value)

例 4.4 捷运公司在下一年度的 1~4 月拟租用仓库堆放物资。已知各月所需仓库面积如表 4.2 所示，仓库租借费用随合同同期而定，期限越长，折扣越大。合同每月初均可办理，规定租用面积和期限。每次办理可签一份合同，也可签若干份租用面积和期限不同的合同，试确定该公司签订租借合同的最优决策，使所付租借费用最小。

表 4.2 所需仓库面积和租借仓库费用数据

月份	1	2	3	4
所需仓库面积 / 100	15	10	20	12
合同租借期限 / 月	1	2	3	4
合同期内的租费 / 元	2800	4500	6000	7300

建立线性规划模型

# 程序文件 ex4_4.py
import cvxpy as cp
x = cp.Variable((4, 4), pos=True)
obj = cp.Minimize(2800*sum(x[:,0])+4500*sum(x[:3,1])+6000*sum(x[:2,2])+7300*x[0,3])
cons = [sum(x[1,:])>=15,
       sum(x[0,1:])+sum(x[2,:3])>=10,
       sum(x[0,2:])+sum(x[1,1:3])+sum(x[2,:2])>=20,
       x[0,3]+x[1,2]+x[2,1]+x[3,0]>=12]
prob = cp.Problem(obj, cons)
prob.solve(solver='GLPK_MI')
print('最优解为：\n',x.value)
print('最优值为：',prob.value)

例 4.5 计算 6 个产地 8 个销地的最小费用运输问题，单位商品运价如表 4.3 所示

表 4.3 单位商品运价表

									产量
	6	4	6	7	4	2	5	9	60
	4	9	5	3	8	5	8	2	55
	5	2	1	9	7	4	3	3	51
	7	6	7	3	9	2	7	1	43
	2	3	9	5	7	2	6	5	41
	5	5	2	2	8	1	4	3	52
销量	35	37	22	32	41	32	43	38

建立线性规划模型

# 程序文件 ex4_5_1.py
import numpy as np
import cvxpy as cp
import pandas as pd
c = np.genfromtxt('data4_5_1.txt', dtype=float, max_rows=6, usecols=range(8)) # 读前 6 行前 8 列数据
e = np.genfromtxt('data4_5_1.txt', dtype=float, max_rows=6, usecols=8)        # 读最后一列数据
d = np.genfromtxt('data4_5_1.txt', dtype=float, skip_header=6)                # 读最后一行数据
x = cp.Variable((6, 8), pos=True)
obj = cp.Minimize(cp.sum(cp.multiply(c, x)))
con = [cp.sum(x,axis=0)==d,
      cp.sum(x,axis=1)<=e]
prob = cp.Problem(obj, con)
prob.solve(solver='GLPK_MI')
print('最优解为：\n', x.value)
print('最优值为：', prob.value)
xd = pd.DataFrame(x.value)
xd.to_excel('data4_5_2.xlsx')                                                 # 数据写到 excel 文件，便于做表使用

# 通过 excel 文件传递数据
# 程序文件 ex4_5_2.py
import cvxpy as cp
import pandas as pd
data = pd.read_excel('data4_5_3.xlsx', header=None)
data = data.values
c = data[:-1, :-1]
d = data[-1, :-1]
e = data[:-1, -1]
x = cp.Variable((6, 8), pos=True)
obj = cp.Minimize(cp.sum(cp.multiply(c, x)))
con = [cp.sum(x,axis=0)==d,
      cp.sum(x,axis=1)<=e]
prob = cp.Problem(obj, con)
prob.solve(solver='GLPK_MI')
print('最优解为：\n', x.value)
print('最优值为：', prob.value)
xd = pd.DataFrame(x.value)
xd.to_excel('data4_5_4.xlsx')

4.2 整数规划（Integer Programming，IP）

4.2.1 整数线性规划模型

 从决策变量的取值范围，可分为

（1）纯整数规划

（2）混合整数规划：决策变量一部分必须取整，另一部分可以不取整

（3）0-1整数规划：决策变量智能取 0 或 1。

例 4.6 背包问题

一个旅行者外出旅行，携带一背包，装一些最有用的东西，共有  件物品供选择。已知每件物品的“使用价值”  和重量 ，要求

（1）最多携带物品的质量为 。

（2）每件物品要么不带，要么只能整件携带。

这是决策问题中比较经典的 0-1 规划问题，要么带要么不带，是一个二值逻辑问题。

例 4.7 指派问题

某单位有  项任务，正好需要  个人去完成，假设分配每个人只能完成一项任务。设  表示分配第  个人去完成第  项任务的费用（时间等），问应如何指派，完成任务的总费用最小？

建立 0-1 整数规划模型：

例 4.8 旅行商问题（货郎担问题）

有一推销员，从城市  出发，要遍历城市  各一次，最后返回  。已知从  到   的旅费为 ，试问应按怎样的次序访问这些城市，使得总旅费最少？

建立 0-1 整数规划模型：

4.2.2 整数线性规划模型的求解

例 4.9 为了生产需要，某工厂的一条生产线需要每天 24h 不间断运转，但是每天不同时段所需要的工人最低数量不同，具体数据如表 4.5 所示。已知每名工人的连续工作时间为 8h。则该工厂应该为生产线配备多少名工人，才能保证生产线的正常运转？

表 4.5 工人数量需求表

班次	1	2	3	4	5	6
时间段	0:00-4:00	4:00-8:00	8:00-12:00	12:00-16:00	16:00-20:00	20:00-24:00
工人数量	35	40	50	45	55	30

# 程序文件 e4_9_1.py
import cvxpy as cp
x = cp.Variable(6, integer=True)
obj = cp.Minimize(sum(x))
cons = [x[0]+x[5]>=35,x[0]+x[1]>=40,
       x[1]+x[2]>=50,x[2]+x[3]>=45,
       x[3]+x[4]>=55,x[4]+x[5]>=30,
       x>=0]
prob = cp.Problem(obj, cons)
prob.solve(solver='GLPK_MI')
print('最优值为：', prob.value)
print('最优解为：', x.value)

# 求余运算
# 程序文件 ex4_9_2.py
import cvxpy as cp
import numpy as np
a = np.array([35, 40, 50, 45, 55, 30])
x = cp.Variable(6, integer=True)
obj = cp.Minimize(sum(x))
cons = [x>=0]
for i in range(6):
    cons.append(x[(i-1)%6]+x[i]>=a[i])
prob = cp.Problem(obj, cons)
prob.solve(solver='GLPK_MI')
print('最优值为：', prob.value)
print('最优解为：', x.value)

例 4.10 某连锁超市经营企业为了扩大规模，新租用 5 个门店，经过装修后再营业。现有 4 家装饰公司分别对 5 个门店的装修费用进行报价，具体数据如表 4.6 所列。为保证这些质量，规定每个装修公司最多承担两个店面的装修任务。为节省装修费用，该企业该如何分配装修任务？

表 4.6 装修费用表

门店	1	2	3	4	5
	15	13.8	12.5	11	14.3
	14.5	14	13.2	10.5	15
	13.8	13	12.8	11.3	14.6
	14.7	13.6	13	11.6	14

建立 0-1 整数规划模型：

# 程序文件 ex4_10.py
import cvxpy as cp
import numpy as np
c = np.loadtxt('data4_10.txt')
x = cp.Variable((4, 5), integer=True)           # 定义决策变量
obj = cp.Minimize(cp.sum(cp.multiply(c, x)))    # 构造目标函数
cons = [0<=x, x<=1, cp.sum(x, axis=0)==1,       # 构造约束条件
        cp.sum(x, axis=1)<=2]
prob = cp.Problem(obj, cons)
prob.solve(solver='GLPK_MI')                    # 求解问题
print('最优解为：\n', x.value)
print('最优值为：', prob.value)

例 4.11 已知 10 个商业网点的坐标如表 4.7 所列，现要在 10 个网点中选择适当位置设置供应站，要求供应站只能覆盖 10km 之内的网点，且每个供应站最多供应 5 个网点，如何设置才能使供应站的数目最小，并求最小供应站的个数。

表 4.7 商业网点 x 坐标和 y 坐标数据

	9.4888	8.7928	11.5960	11.5643	5.6756	9.8497	9.1756	13.1385	15.4663	15.5464
	5.6718	10.3868	3.9294	4.4325	9.9658	17.6632	6.1517	11.8569	8.8721	15.5868

建立 0-1 整数规划模型：

# 程序文件 ex4_11.py
import cvxpy as cp
import numpy as np
a = np.loadtxt('data4_11.txt')
d = np.zeros((10, 10))
for i in range(10):
    for j in range(10):
        d[i, j] = np.linalg.norm(a[:, i]-a[:, j])
x = cp.Variable(10, integer=True)
y = cp.Variable((10, 10), integer=True)
obj = cp.Minimize(sum(x))
cons = [sum(y)>=1, cp.sum(y, axis=1)<=5,
        x>=0, x<=1, y>=0, y<=1]
for i in range(10):
    cons.append(x[i]==y[i, j])
    for j in range(10):
        cons.append(d[i, j]*y[i, j]<=10*x[i])
        cons.append(x[i]>=y[i, j])
prob = cp.Problem(obj, cons)
prob.solve(solver='GLPK_MI')
print('最优值为：', prob.value)
print('最优解为：\n', x.value)
print('----------\n', y.value)

4.3 投资的收益与风险

例 4.12 （选自 1998 年全国大学生数学建模竞赛 A 题）市场上有  种资产（如股票、债券、……） 供投资者选择，某公司有数额为  的一笔相当大的资金可用作一个时期的投资。公司财务分析人员对这  种 资产进行评估，估算出在这一时期内购买资产  的平均收益率为 ，并预测出购买  的风险损失率为 。考虑到投资越分散，总的风险越小，公司确定，当用这笔资金购买若干种资产时，总体风险可用所投资的  中最大的一个风险来度量。

购买  要付交易费，费率为 ，并且当购买额不超过给定值  时，交易费按购买  计算（不买无须付费）。另外，假设同期银行存款利率是 ，且既无交易费又无风险。

已知  时的相关数据如表 4.8 所列。

表 4.8 4 种资产的相关数据

				元
	28	2.5	1	103
	21	1.5	2	198
	23	5.5	4.5	52
	25	2.6	6.5	40

试给该公司设计一种投资组合方案，即用给定的资金 ，有选择地购买若干种资产或存银行生息，使净收益尽可能大，而总体风险尽可能小。

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

模型二：固定盈利水平，极小化风险

模型三：两个目标函数加权求和

模型一求解

# 程序文件 ex4_12_1.py
import cvxpy as cp 
import pylab as plt

b = plt.array([0.025, 0.015, 0.055, 0.026])
c = plt.array([0.05, 0.27, 0.19, 0.185, 0.185])
x = cp.Variable(5, pos=True)
aeq = plt.array([1, 1.01, 1.02, 1.045, 1.065])
obj = cp.Maximize(c @ x)
a = 0; aa = []; Q = []; X = []; M = 10000;
while a < 0.05:
    con = [aeq @ x == M, cp.multiply(b, x[1:]) <= a*M]
    prob = cp.Problem(obj, con)
    prob.solve(solver='GLPK_MI')
    aa.append(a)
    Q.append(prob.value)
    X.append(x.value)
    a = a+0.001
plt.rc('text', usetex=True)
plt.rc('font', size=15)
plt.plot(aa, Q, 'r*')
plt.xlabel('$a$')
plt.ylabel('$Q$', rotation=0)
plt.show()

模型三求解

# 程序文件 ex4_12_2.py
import cvxpy as cp
import numpy as np
import pylab as plt

plt.rc('font', family='SimHei')
plt.rc('font', size=15)
x = cp.Variable(6, pos=True)
r = np.array([0.05, 0.28, 0.21, 0.23, 0.25])
p = np.array([0, 0.01, 0.02, 0.045, 0.065])
q = np.array([0, 0.025, 0.015, 0.055, 0.026])

def LP(w):
    V = []                # 风险初始化
    Q = []                # 收益初始化
    X = []                # 最优解的初始化
    con = [(1+p) @ x[:-1] == 10000, cp.multiply(q[1:], x[1:5]) <= x[5]]
    for i in range(len(w)):
        obj = cp.Minimize(w[i]*x[5]-(1-w[i])*((r-p) @ x[:-1]))
        prob = cp.Problem(obj, con)
        prob.solve(solver='GLPK_MI')
        xx = x.value      # 提出所有决策变量的取值
        V.append(max(q*xx[:-1]))
        Q.append((r-p) @ xx[:-1])
        X.append(xx)
    print('w=', w)
    print('V=', np.round(V, 2))
    print('Q=', np.round(Q, 2))
    plt.figure()
    plt.plot(V, Q, '*-')
    plt.grid('on')
    plt.xlabel('风险（元）')
    plt.ylabel('收益（元）')
    return X
    
w1 = np.arange(0, 1.1, 0.1)
LP(w1)
print('---------------')
w2 = np.array([0.766, 0.767, 0.810, 0.811, 0.824, 0.825, 0.962, 0.963, 1.0])
X = LP(w2)
print(X[-3])
plt.show()

4.4 比赛项目排序问题

例 4.13 （选自 2005 年电工杯数学建模竞赛 B 题）

在各种运动比赛中，为了使比赛公平、公正、合理地举行，一个基本要求是：在比赛项目排序过程中，尽可能使每个运动员不连续参加两项比赛，以便运动员恢复体力，发挥正常水平。

表 4.11 所列为某个小型运动会的比赛报名表。有 14 个比赛项目，40 名运动员参加比赛。表中第 1 行表示 14 个比赛项目，第 1 列表示 40 名运动员，表中 “#” 位置表示运动员参加此项比赛。建立此问题的数学模型，要求合理安排比赛项目顺序，使连续参加两项比赛的运动员人次尽可能地少。

表 4.11 某小型运动会的比赛报名表

1

2

3

4

5

6

7

8

9

10

11

12

13

14

1

#

#

#

#

2

#

#

#

3

#

#

#

4

#

#

#

5

#

#

#

6

#

#

7

#

#

8

#

#

9

#

#

#

#

10

#

#

#

#

11

#

#

#

#

12

#

#

13

#

#

#

14

#

#

#

15

#

#

#

16

#

#

#

17

#

#

18

#

#

19

#

#

20

#

#

21

#

#

22

#

#

23

#

#

24

#

#

#

#

25

#

#

#

26

#

#

27

#

#

28

#

#

29

#

#

#

30

#

#

31

#

#

#

32

#

#

33

#

#

34

#

#

#

#

35

#

#

#

36

#

#

37

#

#

#

38

#

#

#

#

39

#

#

#

#

40

#

#

#

#

TSP 模型可表示为

# 程序文件 ex4_13.py
import numpy as np
import pandas as pd
import cvxpy as cp

a = pd.read_excel("data4_13_1.xlsx", header=None)
a = a.values
a[np.isnan(a)]=0                             # 把空格对应的不确定值替换为0
m,n = a.shape
w = np.ones((n+1, n+1))*10000000             # 邻接矩阵初始化
for i in range(n):
    for j in range(n):
        if i != j:
            w[i, j] = sum(a[:, i]*a[:, j])
for i in range(n):
    w[i, n] = 0
    w[n, i] = 0
wd = pd.DataFrame(w)
wd.to_excel('data4_13_2.xlsx')               # 把邻接矩阵保存到 Excel 文件
x = cp.Variable((n+1, n+1), integer=True)
u = cp.Variable(n+1, integer=True)
obj = cp.Minimize(cp.sum(cp.multiply(w, x)))
con = [cp.sum(x, axis=0) == 1, cp.sum(x, axis=1) == 1,
      x>=0, x<=1, u[0]==0, u[1:]>=1, u[1:]<=n]
for i in range(n+1):
    for j in range(1, n+1):
        con.append(u[i]-u[j]+(n+1)*x[i, j]<=n)
prob = cp.Problem(obj, con)
prob.solve(solver='GLPK_MI')
print('最优值为：', prob.value)
print('最优解为：\n', x.value)
i,j = np.nonzero(x.value)
print('xij=1对应的行列位置如下：')
print(i+1)
print(j+1)

最优值为： 2.0
最优解为：
 [[0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
xij=1对应的行列位置如下：
[ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15]
[ 5 14  6  9 11  2  3  1  8 13  7 10 15 12  4]

TSP 回路为：