隐藏
线性规划单纯形学习笔记 | Bill Yang's Blog

路终会有尽头,但视野总能看到更远的地方。

0%

线性规划单纯形学习笔记

若对线性规划很熟悉的同学可以直接跳过这一段阅读松弛型线性规划。

线性规划

定义

一般线性规划

在一般线性规划问题中,我们希望优化一组线性不等式约束的线性函数。
已知一组实数$a_1,a_2,\ldots,a_n$和一组变量$x_1,x_2,\ldots,x_n$。我们给出定义在这些变量上的线性函数$f$:

如果$b$是一个实数而$f$是一个线性函数,则等式$f(x_1,x_2,\ldots,x_n)=b$是线性等式。
$f(x_1,x_2,\ldots,x_n)\le b,f(x_1,x_2,\ldots,x_n)\ge b$是线性不等式。
我们一般使用线性约束来表示线性等式或线性不等式。
在线性规划中,不允许严格的不等式。

在满足线性约束的条件下,我们往往希望最大化一个线性函数$g$:

我们称线性函数$g$为目标函数

整数线性规划

如果我们在一个线性规划中加入额外的要求,所有的变量都取整数值,那么就得到了一个整数线性规划
仅找出一个可行解就是NP Hard的,因此还没有已知的整数线性规划的多项式算法。
据说若系数矩阵是全幺模矩阵,最优解一定是整数。(我不是很懂,留着后面讲)

标准型与松弛型

标准型

在标准型中,我们已知$n$个实数$c_j$,$m$个实数$b_i$,以及$mn$个实数$a_{ij}$,其中$i=1,2,\ldots,m;j=1,2,\ldots,n$。我们希望找到$n$个实数$x_i$,
最大化:

满足约束条件:

我们称最大化的函数为目标函数,下面的$n+m$个不等式为约束
其中$n$个约束为非负约束。一个任意的线性规划不必有非负约束,但标准型需要。
但是这样的形式不优美,我们将其整理成矩阵形式会更方便。
构造一个$m\times n$的矩阵$A=(a_{ij})$,一个$m$维向量$b=(b_i)$,一个$n$维向量$c=(c_j)$,以及一个$n$维向量$x=(x_j)$,我们可以将标准式整理为矩阵形式:
最大化:

满足约束:

其中目标函数$c^Tx$表示向量的内积,$Ax$是一个矩阵向量乘积,而$x\ge0$表示每个$x$非负。
我们可以用一个元组$(A,b,c)$来表示一个标准型的线性规划。

转换线性规划为标准型

已知一个线性规划满足若干线性约束,要求最小化或最大化它,我们总可以将其转化为标准型。
一个线性规划不是标准型可能有如下$4$种原因。

  • 目标函数可能是最小化,而不是最大化。
  • 可能有变量不具有非负约束。
  • 可能有等式约束,即有一个等号而不是小于等于号。
  • 可能有不等式约束,但不是小于等于号,而是一个大于等于号。

我们逐一将以上线性规划转化为等价的标准型。

  1. 目标函数可能是最小化,而不是最大化。
    将目标函数取负即可。
  2. 可能有变量不具有非负约束。
    若有$x_j$不满足非负约束,那么我们将每个$x_j$出现的位置均用$x_j’-x_j’’$代替,而$x_j’,x_j’’$满足非负约束。
  3. 可能有等式约束,即有一个等号而不是小于等于号。
    用两个不等式约束$f(x_1,x_2,\ldots,x_n)\le b,f(x_1,x_2,\ldots,x_n)\ge b$来替换等式约束$f(x_1,x_2,\ldots,x_n)=b$。
  4. 可能有不等式约束,但不是小于等于号,而是一个大于等于号。
    不等式两边同时取负即可。

举算法导论上的一个例子:



转换线性规划为松弛型

为了利用单纯形算法高效地求解线性规划,我们更喜欢将其表示为某些约束是等式约束的形式。即,只有非负约束是不等式约束,其它约束都是等式约束。
设$\sum_{j=1}^na_{ij}x_j\le b_i$是一个不等式约束。
我们引入新变量$s$,并重写不等式为两个约束:

我们称$s$是一个松弛变量
将一个标准型线性规划转化为松弛型是容易的。

我们将上述例子中的标准型转化为松弛型:

我们称等式左边的变量($x_4,x_5,x_6$)为基本变量,等式右边的变量($x_1,x_2,x_3$)为非基本变量

这样,我们可以用矩阵形式表示松弛型:
我们将常数项写在最前面,并将所有常数项提取出来形成矩阵$a$。
如上述例子写成矩阵形式即为:

不难发现其实矩阵$a$是由标准型的矩阵$A$,向量$b,c$拼接而成。
注意不同之处,矩阵$a$中的矩阵$A’$与标准型的$A$所有对应项均为相反数,这样可以方便我们编程的处理。

单纯形算法

单纯形算法是求解线性规划的经典方法。虽然不是多项式算法,但与椭球法内点法这些多项式算法相比,实用价值更大,且在实际运行中单纯形算法往往非常快速。

单纯形算法与高斯消元法有些类似,每次迭代将系统重写为等价形式,在经过一定次数的迭代后使它的解很容易得到。我们可以把单纯形算法视作不等式上的高斯消元法。

我们从一个例子讲起,考虑下面标准型的线性规划:
最大化:

满足约束:

为了利用单纯形算法,我们必须将标准型转化为松弛型。
我们将上述标准型改写为松弛型:

我们首先构造出一组基本解:将等式右边的非基本变量设为$0$,然后计算等式左边基本变量的值。
在这个例子中,基本解是$(x_1,x_2,\ldots,x_6)=(0,0,0,30,24,36)$,其目标值$z=(3\cdot0)+(1\cdot0)+(2\cdot0)=0$。

接下来我们进行迭代,每次迭代我们将一个基本变量与非基本变量互换并代入约束条件重写松弛型线性规划。
我们先选择一个在目标函数中系数为正的非基本变量$x_e$,而且尽可能增加$x_e$的值而不违反任何约束。变量$x_e$称为基本变量,其他某个变量$x_l$称为非基本变量。其他基本变量和目标函数的值也可能改变。

我们考虑增加上述例子中$x_1$的值。
当增加$x_1$的值时,$x_4,x_5,x_6$的值均减小。因为对每个变量都有非负约束,所以我们不能允许它们中任意一个变成负值。
如果$x_1$增加到$30$以上,$x_4$变成负值,而当$x_1$分半增加到$12,9$以上,$x_5,x_6$分别变成负值。
第三个约束$x_6=36-4x_1-x_2-2x_3 \ge0$是最紧的约束,它限制了我们可以将$x_1$的值增加多少。
因此,我们互换$x_1$和$x_6$,解方程得到:

为了重写右边包含$x_1$的所有等式,我们用上述等式来取代$x_1$,如重写第一个约束可以得到:

类似的,我们可以将上述的系统改写:

我们称此操作为一个转动(pivot)
一次转动选取一个非基本变量$x_e$(称为替入变量)和一个基本变量$x_l$(称为替出变量),然后互换两者的角色。

当前的基本解将非基本变量的值设为$0$,于是得到$(9,0,0,21,6,0)$,目标值为$27$。

我们继续寻找一个可增加值得新变量。我们不想增加$x_6$,因为当它的值增加时,目标值减小。
我们可以尝试添加$x_2,x_3$,设选择$x_3$。
第三个约束限制最紧,要求$x_3$只能增加$\frac32$,因此重写第三个约束,然后将得到的新等式$x_3=\frac32-\frac{3x_2}8-\frac{x_5}4+\frac{x_6}8$替换进上一次的系统中得到新的系统:


这个系统存在关联的基本解$(\frac{33}4,0,\frac32,\frac{69}4,0,0)$,目标值为$\frac{111}4$。现在增加目标值的唯一方法是增加$z_2$。这三个约束分别给出了上界$132,4,\infty$。我们将$x_2$增加到$4$,它变成非基本变量。然后为$x_2$解等式,并代入其他等式得到新系统:

此时我们得到解$(8,4,0,18,0,0)$,目标值为$28$是最优的。
线性规划的系数不必是整数,过渡解也不一定是整数,最终解也不一定是整数,这个例子存在整数解纯属巧合。

单纯形算法实现

  1. 每次选择编号第一个正的非基本变量,将其作为替入变量。若找不到替入变量,说明迭代已结束。

    1
    2
    3
    4
    5
    6
    for(int i=1; i<=m; i++)
    if(dcmp(a[0][i])>0) {
    in=i;
    break;
    }
    if(!in)return a[0][0];
  2. 选择一个对替入变量限制最紧的基本变量作为替出变量。若无法找到说明该线性规划的解无界。

    1
    2
    3
    4
    5
    6
    for(int i=1; i<=n; i++)
    if(dcmp(a[i][in])<0&&a[i][0]/-a[i][in]<Min) {
    Min=a[i][0]/-a[i][in];
    out=i;
    }
    if(!out)throw ; //unbounded
  3. 转动

    1. 将替出变量所在等式重写

      1
      2
      3
      for(int i=0; i<=m; i++) //重置out约束
      if(i!=in)a[out][i]/=-a[out][in];
      a[out][in]=1/a[out][in];
    2. 用重写的等式重新计算其他约束

      1
      2
      3
      4
      5
      6
      for(int i=0; i<=n; i++) { //重新计算其他约束
      if(i==out||dcmp(a[i][in])==0)continue;
      double t=a[i][in];
      a[i][in]=0;
      for(int j=0; j<=m; j++)a[i][j]+=t*a[out][j];
      }

结束!代码就是这么简单!
当然单纯形还有更简洁的写法(可以在$30$行内写完),但是我并不提倡这种写法,因为它不好理解。
在我的代码中间要注意几个细节:

  • 在上述的描述中矩阵是$m$行$n$列,而在程序中为了方便(我)实现使用的是$n$行$m$列(在下面的模板的$init$中进行了转化)。
  • $a$矩阵的系数是标准型矩阵的相反数,也就是和上述松弛型的矩阵形式相同。
  • 本模板没有寻找初始解的操作,我们在之后解释。
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
41
42
43
44
45
46
47
48
const int maxn=10005,maxm=1005;
const double eps=1e-6;

int dcmp(double x) {
if(fabs(x)<=eps)return 0;
return x>eps?1:-1;
}

int n,m;

struct Simplex {
int n,m;
double a[maxn][maxm];
void init(int m,int n) { //a矩阵m行n列
this->n=m;
this->m=n;
}
void pivot(int in,int out) {
for(int i=0; i<=m; i++) //重置out约束
if(i!=in)a[out][i]/=-a[out][in];
a[out][in]=1/a[out][in];
for(int i=0; i<=n; i++) { //重新计算其他约束
if(i==out||dcmp(a[i][in])==0)continue;
double t=a[i][in];
a[i][in]=0;
for(int j=0; j<=m; j++)a[i][j]+=t*a[out][j];
}
}
double Solve() {
while(true) {
int in=0,out=0;
double Min=1e18;
for(int i=1; i<=m; i++)
if(dcmp(a[0][i])>0) {
in=i;
break;
}
if(!in)return a[0][0];
for(int i=1; i<=n; i++)
if(dcmp(a[i][in])<0&&a[i][0]/-a[i][in]<Min) {
Min=a[i][0]/-a[i][in];
out=i;
}
if(!out)throw ; //unbounded
pivot(in,out);
}
}
} fst;

时间复杂度

引理

如果单纯形算法在$C_{n+m}^m$次迭代内不终止,那么它是循环的。
证明:基本变量集合$B$唯一确定了一个松弛型。总共有$n+m$个变量,且$\left|B\right|=m$,所以至多有$C_{n+m}^m$种选择$B$的方式。所以,如果单纯形运行超过$C_{n+m}^m$次迭代,它必然循环。

循环是可能发生的,但很罕见。

Bland规则

如果单纯形算法的寻找替入替出变量总是选择具有最小下标的变量来打破目标值不变的局面,那么单纯形算法必然停止。

因此我们可以使用Bland规则所描述的方法执行单纯形算法(模板已采用),那么单纯形算法的时间复杂度上界是$O(C_{n+m}^m nm)$,然而单纯形算法的时间效率远远达不到其上界。
因此若题目不刻意卡单纯形算法,只要能用网络流做的题用单纯形算法不会超时。

对偶

我们在上面留下了一个疑问,如何寻找一个可行的初始解。
确实可以用辅助线性规划来寻找初始解。
然而在算法竞赛中一般并不会用到辅助线性规划寻找初始解,而是使用对偶的方法强行使 非基本变量为$0$ 是一组可行的初始解。
如果确实要用到寻找初始解请往下翻。

若有一个线性规划形式如下:
最小化目标函数:

满足约束:

此时我们发现形式与标准式不符,我们通过前文所讲的调整使其变为标准式:

最大化目标函数:

满足约束:

我们发现对于限制条件$-Ax\le -B$,非基本变量为$0$的初始解并不能满足。

此时我们可以通过将问题对偶,使该线性规划变成如下形式:

最大化目标函数:

满足约束:

我们惊讶的发现问题竟然变成了标准式,且满足非基本变量为$0$的初始解。

没错这就是对偶,我们可以保证经过对偶后的问题与原问题具有相同的解。

具体的证明就是一个配凑系数的过程,此处就不证明了。

将问题对偶后我们就可以用松弛型解决线性规划了。

当然我们有更方便的方法:
我们可以先不考虑对偶后的问题,直接强行将其转化为松弛型线性规划,列出简化后的矩阵$a$,并将其转置得到矩阵$a^T$,通过单纯形解决松弛型$a^T$即可。

寻找初始解

可以使用算导上的辅助线性规划的做法。
也可以使用candy?说明的随机初始化
在所有常数项$\lt0$的约束中随机选一个作为替出变量,再随机选一个系数$\lt0$的作为替入变量,然后$pivot$一下后常数项就变正了。

注意这样的寻找初始解方法可能会死循环,这就是玄学了。

1
2
3
4
5
6
7
8
9
10
11
12
13
bool find() {
while(true) {
int in=0,out=0;
for(int i=1; i<=n; i++)
if(dcmp(a[i][0])<0&&(!out||(rand()&1)))out=i;
if(!out)break;
for(int j=1; j<=m; j++)
if(dcmp(a[out][j])>0&&(!in||(rand()&1)))in=j;
if(!in)return false; //无解
pivot(in,out);
}
return true;
}

全幺模矩阵

若系数矩阵$A$是全幺模矩阵,那么可以保证所有过程均为整数,此时我们可以去除系数为$0$的项进行优化,并且不必使用double。
此时线性规划的解与整数线性规划的解相同。

如何证明一个矩阵是全幺模矩阵?
首先全幺模矩阵中的元素$\in \lbrace -1,0,1\rbrace$。

充分条件1

能将矩阵的行划分为两个不相交集合$N,M$满足:
每列最多两个非$0$元素。
若$1$列上有$2$个异号元素,那么它们分别属于$N,M$。
若$1$列上有$2$个同号元素,那么它们同时属于$N$或$M$。

充分条件2

能表示成网络流。
当然能用网络流了就没必要写单纯形了。(往往网络流比单纯形更快,加了对幺模矩阵的优化除外)

充分条件3

有一种列的排列能使得每行的$1$是连续的。

充分条件4

对于每列单调不增,且任意$2\times2$的子矩阵行列式的值为$-1,0,1$。

性质1

单位矩阵是全幺模矩阵

性质2

全幺模矩阵的乘积是全幺模矩阵

性质3

全幺模矩阵的逆是全幺模矩阵

性质4

进行线性规划的所有中间变量和解均为整数。

参考资料

姥爷们赏瓶冰阔落吧~