前段时间参加了华为的2017软件精英挑战赛,用到了单纯形算法求解线性规划问题,学习了正单纯形,对偶单纯形以及割平面法并用C语言实现了完整的simplex算法。
单纯形算法是用来求解线性规划问题的,其被用在了众多SMT Solver中,如Yices, Z3等,都是基于的单纯形算法,其算法效率在最坏情况下为指数级别,不过实现证明其效率在大多数情况下都是令人满意的。
标准形式
max{b+∑1≤k≤nckxk}
∑1≤k≤nA1xk≤b1
…
∑1≤k≤nAmxk≤bm
x1,...,xn≥0
对于 b1,...,bm≥0 的情况,可使用正单纯形算法求解,对于存在 bi<0 的情况,则需要使用对偶单纯形求解。这里求出来的解都是非整数解,若要得到目标式的整数最优解,在前面的基础上还需要使用割平面法,即引入Gomory(高莫雷)约束。
一些特殊情况的处理
1.对于不满足
xi≥0
的情况,可将
xi
化为两个非负数的差值,即
xi=xj−xk
, 其中
xj,xk≥0
2.如果限制式子为等式,可将其化为一个大于等于和小于等于,如
a+b=c⇔(a+b≤c∧a+b≥c)
3.如果要求的目标式
Z
的最小值,则可转化为先求
正单纯形算法
首先将标准形式写为如下形式(Slack Form)
Z=b+∑1≤k≤nckxk
y1=b1−∑1≤k≤nA1xk
…
ym=bm−∑1≤k≤nAmxk
这里引入了新的松弛变量
y1,...,ym≥0
,不难证明,新的形式跟原不等式等价。
xk
称为basic variable(基变量),
ym
称为non-basic variable(非基变量)。因为
b≥0
,所以所有的非基变量为0为一组可行解,我们就是从这个初始可行解触发,一步一步找到最优解。下面是具体的步骤:
1.查找
ck
且
ck>0
,若不存在(目标式已经无法再继续增大),跳转到第4步
2.找到
minbiAi
,即限制条件最强的一组式子
3.交换
xk
和
yi
,
xk
变为新的基变量,
yi
变为新的非基变量,将其他式子中的
xk
替换为新的值,然后回到第1步
4.将所有的非基变量取0值,即可算出所有的基变量的值以及目标式的最大值,算法结束
单纯形的物理意义
如上图所示,为7个限制不等式所构成的凸多面体,每个限制条件可看成正凸多面体的一个面,凸多面体的内部(包括表面)即为可行解区域。可以证明,最优解如果存在,则一定在凸多面体的顶点上。对于三维空间来说(更高维空间可类似推广),某个顶点比为三个或以上的平面的交点,这个交点满足对应相交平面的限制不等式。如点A,同时满足式子②③⑦。
算法中第3步中的交换因子(pivot),可以看成是从凸多面体的一个顶点走向另一个顶点。
对偶单纯形
在标准形式(Slack Form)中,如果
bm<0
,则不能直接用正单纯形求解,需要使用对偶单纯形,先将标准型化为正单纯形的标准形式,即
bm≥0
,然后再使用正单纯形进行求解。
对偶单纯形的求解步骤如下:
1.引入新的松弛变量
z≥0
,得到如下新的标准型
Z=b+∑1≤k≤nckxk−z
y1=b1−∑1≤k≤nA1xk+z
…
ym=bm−∑1≤k≤nAmxk+z
2.取
min{bi},0≤i≤m
,交换
z
与
3.为了使得加入
4.当
z
再次从基变量置换为非基变量后,此时
割平面法求整数最优解
以上求得的都是非整数最优解,即LP(Linear Programing),如果要求整数最优解(Interger Linear Programing),则需要在得到非整数最优解的基础上,对其进行切割,以得到整数解,这里采用割平面法。
设求出的最优解:
设
xi
不为整数,
xi=bi−∑1≤k≤nAkxk
,
xi
为基变量,
xk
为非基变量
将
bi
与
Ak
分离成一个整数与真分数之和,即
bi=[bi]+fi
,
Ak=[Ak]+fk
其中
[bi],[Ak]
为整数,
fi,fk
为真分数,则有
xi=[bi]+fi−∑1≤k≤n[Ak]xk−∑1≤k≤nfkxk
把整数项移到左边,分数项移到右边,得
xi−[bi]+∑1≤k≤n[Ak]xk=fi−∑1≤k≤nfkxk
因为等式两边都为整数,且
fi−∑1≤k≤nfkxk≤fi<1
则
fi−∑1≤k≤nfkxk≤0
, 这就是高莫雷(Gomory)约束。
对上式引入一个新的松弛变量
gi
得
gi=−fi+∑1≤k≤nfkxk
将其添加到原式中,使用对偶单纯形进行求解。如果算出解仍不为整数解,则选取新的Gomory约束,继续进行切割
用C实现的源码地址如下:
(https://github.com/zjc666/Simplex_Algorithm)