介绍
已知:
- 一个 m × n m×n m×n大小的实矩阵 A A A。
- 大小为 m m m的实向量 b b b。
- 大小为 n n n的实向量 c c c。
在线性规划的目标中,我们的目标是找到一个大小为 n n n的实向量 x x x,使得① A x ≤ b Ax≤b Ax≤b,② c T x c^Tx cTx最大。( c T c^T cT是 c c c的转置矩阵)
分析:
- 对于 A x ≤ b Ax≤b Ax≤b,左侧可以看做是 x x x中个元素的线性组合。
- 右侧则是对左侧每一组线性组合的约束。
- ≤ b ≤b ≤b说明这个约束是通过指定上界来进行约束。
- 实际上也可以用 = = =或 ≥ ≥ ≥来约束。
- 在这样的约束下,所追求的是目标函数 f = c T x f=c^Tx f=cTx的最大值(也可能最小值)。
- 而目标函数同样也是 x x x中各元素的线性组合。
对于具体的问题来说,我们需要分析问题,从中找到已知条件: A 、 b 、 c A、b、c A、b、c,然后就可以求解。
求解
求解算法本文不作介绍,而是直接利用Python中的cvxpy库进行求解。(库的安装,参考官网。cvxpy使用指南。)
案例1
从文库中随意找到一个线性规划的案例。该案例提供了答案,可以用来验证。
案例截图如下:
首先从案例中提取出已知条件中的 A 、 b 、 c A、b、c A、b、c。
其中 A 、 b A、b A、b与条件约束有关,可能有多组。因为cvxpy支持多组约束的计算,所以为了问题的清晰表述,将分成了三组约束:
- 相等约束
根据 x 1 + x 4 = 400 ; x 2 + x 5 = 600 ; x 3 + x 6 = 500 x_1+x_4=400;x_2+x_5=600;x_3+x_6=500 x1+x4=400;x2+x5=600;x3+x6=500得出:
# 约束1
A1 = np.array([[1, 0, 0, 1, 0, 0],
[0, 1, 0, 0, 1, 0],
[0, 0, 1, 0, 0, 1],
])
b1 = np.array([400, 600, 500])
- 小于等于约束
根据 0.4 x 1 + 1.1 x 2 + x 3 ≤ 800 ; 0.4 x 1 + 1.1 x 2 + x 3 ≤ 800 ; 0.5 x 4 + 1.2 x 5 + 1.3 x 6 ≤ 900 0.4x_1+1.1x_2+x_3≤800;0.4x_1+1.1x_2+x_3≤800;0.5x_4+1.2x_5+1.3x_6≤900 0.4x1+1.1x2+x3≤800;0.4x1+1.1x2+x3≤800;0.5x4+1.2x5+1.3x6≤900得出:
# 约束2
A2 = np.array([[0.4, 1.1, 1, 0, 0, 0],
[0, 0, 0, 0.5, 1.2, 1.3],
])
b2 = np.array([800, 900])
- 大于等于约束
根据 x i ≥ 0 , i = 1 , 2 , . . . , 6 x_i≥0,i=1,2,...,6 xi≥0,i=1,2,...,6得出:
# 约束3
A3 = np.array([[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 1],
])
b3 = np.array([0, 0, 0, 0, 0, 0])
根据 m i n z = 13 x 1 + 9 x 2 + 10 x 3 + 11 x 4 + 12 x 5 + 8 x 6 minz = 13x_1+9x_2+10x_3+11x_4+12x_5+8x_6 minz=13x1+9x2+10x3+11x4+12x5+8x6,提取出目标函数结构中的 c c c:
# 用于目标函数
c = np.array([13, 9, 10, 11, 12, 8])
最后的求解代码:
import cvxpy as cp
import numpy as np
# ...把前面三个约束、向量c的代码加上
# 定义自变量
n = 6 # 有6个自变量:x1,x2,x3,x4,x5,x6
x = cp.Variable(n)
# 定义问题,把三个约束条件添加上
prob = cp.Problem(cp.Minimize(c.T @ x),
[A1 @ x == b1, A2 @ x <= b2, A3 @ x >= b3])
# 解决问题
prob.solve()
# 输出结果
print("\n目标函数的最小值", prob.value)
print("所求的x矩阵")
print(x.value)
# 对x向量各元素取整数后再输出
for item in x.value:
print(round(item))
最后向量 x x x的整数解(与文库中提供的答案结果一致):
0
600
0
400
0
500
目标函数的最小值:13800.00000202733
案例2
案例来源:该博客中的第一个案例(营业员上班时间安排的问题)。
import cvxpy as cp
import numpy as np
A = np.array([[1, 0, 0, 1, 1, 1, 1],
[1, 1, 0, 0, 1, 1, 1],
[1, 1, 1, 0, 0, 1, 1],
[1, 1, 1, 1, 0, 0, 1],
[1, 1, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 1],
[1, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 1],
])
b = np.array([300, 300, 350, 400, 480, 600, 500, 0, 0, 0, 0, 0, 0, 0])
c = np.array([1, 1, 1, 1, 1, 1, 1])
n = 7
x = cp.Variable(n)
prob = cp.Problem(cp.Minimize(c.T @ x),
[A @ x >= b])
prob.solve()
# Print result.
print("\nThe optimal value is", prob.value)
print("A solution x is")
for item in x.value:
print(round(item))
print("A dual solution is")
print(prob.constraints[0].dual_value)