线性规划问题
首先引入如下的问题:
假设食物的各种营养成分、价格如下表:
Food
Energy(能量)
Protein(蛋白质)
Calcium(钙)
Price
Oatmeal(燕麦)
110
4
2
3
Whole milk(全奶)
160
8
285
9
Cherry pie(草莓派)
420
4
22
20
Pork with beans(猪肉)
260
14
80
19
要求我们买的食物中,至少要有2000的能量,55的蛋白质,800的钙,怎样买最省钱?
设买燕麦、全奶、草莓派、猪肉为x1,x2,x3,x4
于是我们可以写出如下的不等式组
其实这些不等式组就是线性规划方程(Linear programming formulation):
简单的说,线性规划就是在给定限制的情况下,求解目标。
可行域
来看一个算法导论中的例子,考虑如下的线性规划:
max x1+x2s.t 4×1–x2≤8 2×1+x2≤10 5×1–2×2≥−2 x1,x2≥0
我们可以画出下面的图:
看图a,灰色的区域就是这几个约束条件要求x1,x2所在的区域,而我们最后的解x1,x2也要在这里面。我们把这个区域称为可行域(feasible region)
图b可以直观的看出,最优解为8, 而 x1= 2 , x2=6
线性规划标准形式
线性规划的标准形式如下:
min cTxs.t Ax≤b x≥0
就是
求的是min
所有的约束为<=的形式
所有的变量均 >=0
如何变为标准形式?
原来是max, 直接*-1求min
若原来约束为=,转为 >= 和<=
约束原来为 >= 同样的*-1,就改变了<=
若有变量 xi < 0 ,那么用 x‘ – x”来替代,其中 x’>=0 x”>=0
线性规划松弛形式
松弛形式为:
min cTx
s.t. Ax=b
x≥0
就是通过引入变量把原来的 <= ,变为=的松弛形式.
如:
x1+x2≤2 x1+x2≥1×1,x2≥0
写为松弛形式就是
x1+x2+x3=2×1+x2+x4=1×1,x2,x3,x4≥0
<= vs <
为什么我们的线性规划的形式都是可以 <= 或者 >=的形式的?把等号去掉可以么?不可以
举个例子
max x
s.t. x≤1
max x
s.t. x<1
显然第二个是无解的。
单纯形算法的思想与例子
如何求解线性规划问题呢?
有一些工具如GLPK,Gurobi 等,不在本文的介绍范围内。
本文要介绍的是单纯形算法,它是求解线性规划的经典方法,虽然它的执行时间在最坏的情况下是非多项式的(指数时间复杂度),但是,在绝大部分情况下或者说实际运行过程中却是多项式时间。
它主要就三个步骤
找到一个初始的基本可行解
不断的进行旋转(pivot)操作
重复2直到结果不能改进为止
以下面的线性规划为例:
min −x1−14×2−6x3s.t. x1+x2+x3≤4 x1≤2 x3≤3 3×2+x3≤6 x1,x2,x3≥0
将其写为松弛的形式:
min −x1−14×2–6x3s.t. x1+x2+x3+x4=4 x1+x5=2 x3++x6=3 3×2+x3+x7=6 x1,x2,x3,x4,x5,x6,x7≥0
其实,就是等价于(仍然要求 x1,x2,x3,x4,x5,x6,x7 >=0):
z=−x1−14×2–6x3x4=4−x1–x2−x3x5=2–x1x6=3–x3x7=6–3×2–x3
在上述的等式的左边称为基本变量,而右边称为非基本变量。
现在来考虑基本解就是把等式右边的所有非基本变量设为0,然后计算左边基本变量的值。
这里,容易得到基本解为:(x1,x2….x7) = (0,0,0,4,2,3,6),而目标值z = 0,其实就是把基本变量xi设置为bi。
一般而言,基本解是可行的,我们称其为基本可行解。初始的基本解不可行的情况见后面的讨论,这里假设初始的基本解就是基本可行解,因此三个步骤中第一步完成了。
现在开始,来讨论上面的第二个步骤,就是旋转的操作。
我们每次选择一个在目标函数中的系数为负的非基本变量xe,然后尽可能的增加xe而不违反约束,并将xe用基本变量xl表示, 然后把xe变为基本变量,xl变为非基本变量。
这里,假设我们选择增加x1,那么在上述的等式(不包括目标函数z那行)中,第1个等式限制了x1 <=4(因为x4>=0),第2个等式有最严格的限制,它限制了x1 <=2,因此我们最多只能将x1增加到2,根据上面的第二个等式,我们有: x1 = 2 – x5,带入上面的等式就实现了xe和xl的替换:
z=−2−14×2–6×3+x5x4=2–x2−x3+x5x1=2–x5x6=3–x3x7=6–3×2–x3
这样其实就是一个转动(pivot)的过程,一次转动选取一个非基本变量(也叫替入变量)xe 和一个基本变量(也叫替出变量) xl ,然后替换二者的角色。执行一次转动的过程与之前所描述的线性规划是等价的。
同样的,将非基本变量设为0,于是得到:(x1,x2….x7) = (2,0,0,2,0,3,6), Z = -2,说明我们的目标减少到了-2
接下来是单纯形算法的第三步,就是不断的进行转动,直到无法进行改进为止,继续看看刚才的例子:
我们接着再执行一次转动,这次我们可以选择增大x2或者x3,而不能选择x5,因为增大x5之后,z也增大,而我们要求的是最小化z。假设选择了x2,那么第1个等式限制了x2 <=2 , 第4个等式限制了x2 <= 2,假设我们选择x4为替出变量,于是有: x2 = 2 – x3 – x4 + x5 ,带入得:
z=−30+8×3+14×4−13×5
x2=2−x3−x4+x5x1=2–x5x6=3–x3x7=2×3+3×4–3×5
此时,我们的基本解变为(x1,x2….x7) = (2,2,0,0,0,3,0), Z = -30
我们可以继续的选择增大x5,第4个等式具有最严格的限制(0 – 3×5 >=0),我们有x5 = 2/3 x3 + x4 – 1/3 x7
带入得
z=−30–23×3+x4+133x7x2=2–13×3−13x7x1=2–23×3–x4+13x7x6=3–x3x5=23×3+x4–13×7
此时,我们的基本解变为(x1,x2….x7) = (2,2,0,0,0,3,0), Z = -30,这时候并没有增加,但是下一步,我们可以选择增加 x3。第2个和第3个有最严格的限制,我们选第2个的话,得:x3 = 3 – 3/2 x1 – 3/2 x4 + 1/2 x7,然后老样子,继续带入:
z=−32+x1+2×4+4x7x2=1+12×1+12×4–12x7x3=3–32×1–32×4+12x7x6=32×1+32×4–12x7x5=2–x1
现在,已经没有可以继续增大的值了,停止转动,z=-32就是我们的解,而此时,基本解为:(x1,x2….x7) = (0,1,3,0,2,0,0),看看最开始的目标函数:z = -x1 -14×2 – 6×3 ,我们将x2=1,×3=3带入得,z=-32,说明我们经过一系列的旋转,最后得到了目标值。
退化(Degeneracy)
在旋转的过程中,可能会存在保持目标值不变的情况,这种现象称为退化。比如上面的例子中,两次等于-30.
可以说退化可能会导致循环(cycling)的情况,这是使得单纯形算法不会终止的唯一原因。还好上面的例子中,我们没有产生循环的情况,再次旋转,目标值继续降低。
《算法导论》是这样介绍退化产生循环的:
Degeneracy can prevent the simplex algorithm from terminating, because it can lead to a phenomenon known as cycling: the slack forms at two different iterations of SIMPLEX are identical. Because of degeneracy, SIMPLEX could choose a sequence of pivot operations that leave the objective value unchanged but repeat a slack form within the sequence. Since SIMPLEX is a deterministic algorithm, if it cycles, then it will cycle through the same series of slack forms forever, never terminating.
如何避免退化?一个方法就是使用Bland规则:
在选择替入变量和替出变量的时候,我们总是选择满足条件的下标最小值。
替入变量xe:目标条件中,系数为负数的第一个作为替入变量
替出变量xl:对所有的约束条件中,选择对xe约束最紧的第一个
在上面的例子中,我也是这么做的。^ ^
另一个方法是加入随机扰动。
无界(unbounded)的情况
有的线性规划问题是无界的,举个栗子对于下面的线性规划
min−x1–x2s.t. x1–x2≤1 −x1+x2≤1 x1,x2≥0
画出区域为:
显然可以不断的增大。让我们来看看单纯形算法是如何应对的:
上述的写成松弛形式为:
min −x1–x2s.t. x1–x2+x3=1 −x1+x2+x4=1 x1,x2,x3,x4≥0
也就是,
z=−x1–x2x3=1–x1+x2x4=1+x1–x2
选择x1 为替入变量,x3为替出变量,有:
z=−1–2×2+x3x1=1+x2–x3x4=2–x3
这时候我们只能选择x2 为替入变量,才能使得目标值变小,但是我们发现,对于x2没有任何的约束,也就是说,x2可以无限大,所以这是没有边界的情况。
这个情况是我们有一个替入变量,但是找不到一个替出变量导致的,这时候就是无界的情况了,写算法的时候注意判断一下即可。
从几何角度看单纯形算法
上面我们介绍单纯形算法的时候,是通过最直观的等式变换(就是旋转操作)介绍的。
我们知道,线性规划就是在可行域围成的多胞形中求解,现在从几何的视图来看看单纯形算法。
只需考虑顶点
让我再次召唤之前的图:
直观上看,最优解就在顶点上,不需要考虑内部点。
一个引入的证明
我们假设x(0) 是最优解,连接x(1)和x(0) 与 x(2)和x(3)相交于点x’
我们可以把x(0) 分解,x(0) = λ1 x(1) + (1 – λ1)x’ 其中λ1 = p / (p + q)
同样的把x‘ 分解,x’ = λ2 x(2) + (1 – λ2)x(3) 其中λ2 = r / (r + s)
因此有:x(0) = λ1 x(1) + (1 – λ1)λ2 x(2) + (1 – λ1) (1 – λ2)x(3),而λ1 + (1 – λ1)λ2 + (1 – λ1) (1 – λ2) = 1
设 cT x(1) 小于等于 cT x(2), cT x(3),因此有:
CTx(0)=λ1CTx(1)+(1–λ1)λ2CTx(2)+(1–λ1)(1–λ2)CTx(3) ≥λ1CTx(1)+(1–λ1)λ2CTx(1)+(1–λ1)(1–λ2)CTx(1) =CTx(1)
因此,x(1) 并不比x(0) 差。
小结
用几何的角度看待单纯形算法,主要有几点:
最优解可以在顶点上找到,不许考虑内部点
一次旋转就是一个顶点沿着一条边λ走θ倍到另一个顶点的过程
当然也需要注意初始化单纯形算法,比如之前的情况:
我们的顶点要在可行域才行,而不要跑到(0,0)去了。初始方法和之前的一样。
单纯形算法的调用(Python内置工具包)
python真的是非常强大。scipy包里面包含了很多科学计算相关的模块方法。
原题目:有2000元经费,需要采购单价为50元的若干桌子和单价为20元的若干椅子,你希望桌椅的总数尽可能的多,但要求椅子数量不少于桌子数量,且不多于桌子数量的1.5倍,那你需要怎样的一个采购方案呢?解:要采购x1张桌子,x2把椅子max z= x1 + x2s.t. x1 – x2 <= 01.5×1 >= x250×1 + 20×2 <= 2000x1, x2 >=0在python中此类线性规划问题可用lp solver解决scipy.optimize._linprog def linprog(c: int, A_ub: Optional[int] = None, b_ub: Optional[int] = None, A_eq: Optional[int] = None,b_eq: Optional[int] = None, bounds: Optional[Iterable] = None, method: Optional[str] = simplex, callback: Optional[Callable] = None,options: Optional[dict] = None) -> OptimizeResult矩阵A:就是约束条件的系数(等号左边的系数)矩阵B:就是约束条件的值(等号右边)矩阵C:目标函数的系数值from scipy import optimize as optimport numpy as np#参数#c是目标函数里变量的系数c=np.array([1,1])#a是不等式条件的变量系数a=np.array([[1,-1],[-1.5,1],[50,20]])#b是是不等式条件的常数项b=np.array([0,0,2000])#a1,b1是等式条件的变量系数和常数项,这个例子里无等式条件,不要这两项#a1=np.array([[1,1,1]])#b1=np.array([7])#限制lim1=(0,None)#(0,None)->(0,+无穷)lim2=(0,None)#调用函数ans=opt.linprog(-c,a,b,bounds=(lim1,lim2))#输出结果print(ans)#注意:我们这里的应用问题,椅子不能是0.5把,所以最后应该采购37把椅子
手动自己实现一个简单的单纯性算法
import numpy as npclass Simplex(object):
def __init__(self, obj, max_mode=False):
self.max_mode = max_mode # 默认是求min,如果是max目标函数要乘-1
self.mat = np.array([[0] + obj]) * (-1 if max_mode else 1) #矩阵先加入目标函数 def add_constraint(self, a, b):
self.mat = np.vstack([self.mat, [b] + a]) #矩阵加入约束 def solve(self):
m, n = self.mat.shape # 矩阵里有1行目标函数,m – 1行约束,应加入m-1个松弛变量
temp, B = np.vstack([np.zeros((1, m – 1)), np.eye(m – 1)]), list(range(n – 1, n + m – 1)) # temp是一个对角矩阵,B是个递增序列
mat = self.mat = np.hstack([self.mat, temp]) # 横向拼接
while mat[0, 1:].min() < 0: #判断目标函数里是否还有负系数项
col = np.where(mat[0, 1:] < 0)[0][0] + 1 # 在目标函数里找到第一个负系数的变量,找到替入变量
row = np.array([mat[i][0] / mat[i][col] if mat[i][col] > 0 else 0x7fffffff for i in
range(1, mat.shape[0])]).argmin() + 1 # 找到最严格约束的行,也就找到替出变量
if mat[row][col] <= 0: return None # 若无替出变量,原问题无界
mat[row] /= mat[row][col] #替入变量和替出变量互换
ids = np.arange(mat.shape[0]) != row
mat[ids] -= mat[row] * mat[ids, col:col + 1] # 对i!= row的每一行约束条件,做替入和替出变量的替换
B[row] = col #用B数组记录替入的替入变量
return mat[0][0] * (1 if self.max_mode else -1), {B[i]: mat[i, 0] for i in range(1, m) if B[i] < n} #返回目标值,对应x的解就是基本变量为对应的bi,非基本变量为0″””
minimize -x1 – 14×2 – 6×3
st
x1 + x2 + x3 <=4
x1 <= 2
x3 <= 3
3×2 + x3 <= 6
x1 ,x2 ,x3 >= 0
answer :-32
“””
t = Simplex([-1, -14, -6])
t.add_constraint([1, 1, 1], 4)
t.add_constraint([1, 0, 0], 2)
t.add_constraint([0, 0, 1], 3)
t.add_constraint([0, 3, 1], 6)
print(t.solve())
print(t.mat)
小结
给定一个线性规划L,就只有如下三种情形:
有一个有限目标值的最优解
不可行
无界
转载:细语呢喃 > 线性规划-单纯形算法详解本文在转载的过程中加入了一些自己代码的实现,想看原作者详情的请移步到:https://www.hrwhisper.me/introduction-to-simplex-algorithm/