线性规划入门·单纯形

许久之前开始学的算法了…今天终于学完了。主要参考的材料是董克凡的$2016$集训队论文和$\rm{Candy?}$的代码整理,在此处致以敬意。


$\rm{0x01\quad Preface}$

首先明确,线性规划主要解决的问题是:

$$
\text{最大化}\quad \sum c_ix_i \quad(i = 1, 2,3 \cdots n) \\
\text{满足约束} \quad \sum \limits_{j=1}^{n} a_{i,j}x_{j} \leq b_i\quad (i = 1,2,3\cdots m) \\ \qquad \qquad \qquad \quad x_j \geq 0\quad (j=1,2,3\cdots n)
$$
不失一般性的,我们定义最大化的函数为目标函数($\rm{Aim-Func}$) ,定义约束函数的集合为约束函数集($\rm{Constraint-Set}$)。那么朴素的线性规划可以看求一组向量{$x_1,x_2\cdots x_n$},使之既可以做约束函数的因变量,又满足其目标函数的值为$max$。

同时,规定所有的$x$均满足$x>0$。

当然,存在一种更加赏心悦目的矩阵表示方式:
$$
\text{最大化}\quad \boldsymbol{c^Tx} \\
\text{满足约束} \quad A\boldsymbol{x\leq b} \\\
\it{\qquad \qquad \quad} \boldsymbol {x} \geq 0
$$
其中$\boldsymbol{c,x,b}$均为一维向量,$A$为系数矩阵。

那么我们在高中数学必修五里面运用的智障做法是,通过每个约束确定一个凸包,再用目标函数不断平移以求得与凸包的切点/切线来达到最大值,但是这样的做法通常不具有一般性。所以需要引入一种更常用的方法来解决这类问题。

$\rm{0x02\quad}$ 松弛型矩阵与$\rm{Pivot}$操作

我们考虑一种更加友好的线性规划方式——松弛型矩阵,即将原来的矩阵添加几个无实际作用的新变量$x_i$,使之换一个样子:
$$
\text{最大化}\quad \sum c_ix_i \quad(i = 1, 2,3 \cdots n) \\
\text{满足约束} \quad x_{n+i} = b_i - \sum \limits_{j=1}^{n} a_{i,j}x_{j} \quad (i = 1,2,3\cdots m) \\ \qquad\quad x_j \geq 0\quad (j=1,2,3\cdots n+m)
$$
两种表示是等价的,但是我们更倾向于松弛型这种简洁的表述方式。

同时我们规定以下:

$1.$基变量:在松弛型约束中,等式左边的变量。

$2.$非基变量:在松弛型约束中,等式右边的变量。

那我们定义一次转轴操作$(pivot)$将一个基变量换进等式右边,换出一个非基变量的过程

假设我们在第$i$组约束中,有一个变量$x_k~(k>n)$是基变量,我们要换出一个非基变量$x_p~(1\leq p\leq n)$,那么就会由:
$$
x_k= b_i -\sum \limits_{j=1}^{n} a_{i,j}x_{j}
$$
变成
$$
x_p = \frac{b_i - \sum \limits_{j\neq p}a_{i,j}x_j -x_k}{a_{i,p}}
$$
而单纯性做法的本质就是通过不断转轴,实现目标函数不断变大。

看上去似乎比较抽象?我们考虑对一次转轴操作,我们需要在转轴之后将原来目标函数里的非基变量代入,即用$(b_i - \sum _{j\neq p}a_{i,j}x_j -x_k)/a_{i,p}$这个东西去替换掉$x_p$,那么其中的常数项$\frac{b_i}{a_{i,p}}$就可以作为目标函数中的一个常数值,当$x_i~(i=1,2,3\cdots n)$均为$0$时,目标函数的值即为此。那么转轴操作就是通过这样的操作使得目标函数里的常数值不断增大,达到最优解。

注意,转轴之后的得到的$x_p$的表达式,不仅要代入目标函数,也要代入其余的约束。

$\rm{0x03}\quad \rm{Simplex}$

那么接下来,我们考虑单纯性做法的完整过程。

不失一般性的,我们假设所有$b_i\geq 0$

首先,我们对于转轴操作结束,回代一次之后,会发现目标函数中肯定会有至少一项系数变为负值,即转进来的前·基变量$x_k$,那么增大$x_k$一定会让结果目标函数变小。所以我们可以断言,当目标函数里的所有变量系数均为负值时,目标函数的最优值就会是现在目标函数中的常数值——因为我们在前文已经假定$x_i\geq 0$了。

同样,我们每次转轴操作需要保持原来线性规划的不变性,换句话说就是我们每次转轴时都需要找一个对与某个非基变量限制最紧的约束,将其换出。原因是我们考虑目标函数中的$x_i$,当其系数$\geq 0$时,$\rm{Aim }\propto x_i$,所以我们需要找一个最紧的约束遏止$x_i$的增长(即使我们不想)。

那么伪代码如下,$A,\boldsymbol{b,c}$的定义一开始已给出:


Simplex(A, b, c){
$\qquad$ initialization(A,b,c) ;
$\qquad$ while $∃e$ that $c_e>0${
$\qquad$$\qquad$find the index $l$ that $A_{l,e}$ $> 0$ and minimizes $b_l/A_{l,e}$
$\qquad$$\qquad$if $∀l, A_{l,e} ≤ 0$
$\qquad$$\qquad$$\qquad$return Unbounded ;
$\qquad$$\qquad$else
$\qquad$$\qquad$$\qquad$pivot(A, b, c, l, e);
}


值得注意的是,如果在寻找完最紧约束后,存在一个$A_{l,e} \leq 0$,那么我们在这组约束里无论怎么增大$x_e$都会使得目标函数增大,所以直接return Unbounded

此时我们还需要进行额外的操作,即初始化。因为我们转轴操作的目的是使目标函数增大,所以我们需要保证所有的$b_i\geq 0$,但一开始给定的$b_i$很可能不满足这一约束,所以我们需要进行一次初始化。那么对于初始化,我们的思想是做另一次线性规划。鬼知道那些发明这玩意儿的神仙怎么构造出的这个线性规划。

我们考虑一个这样的线性规划:
$$
\text{最大化}\quad -x_0 \\
\text{满足约束} \quad x_{n+i} = b_i - \sum \limits_{j=1}^{n} a_{i,j}x_{j} +x_0\quad (i = 1,2,3\cdots m)\\
\it{\qquad \qquad \qquad \quad x_j} \geq 0\quad (j=0,1,2,3\cdots n)
$$
首先,无论怎样,先思考他会不会对原来的线性规划产生什么影响——答案是:不会。因为我们考虑当$-x_0$最大时,一定是$x_0=0$的时候。所以对于约束中的$x_{n+i}$,大小不会变化。

那么我们考虑这种初始化方式的正确性。我们考虑每次都将$x_0$作为换入变量(即我们一开始将$x_0$视作非基变量),那么$pivot$之后就会有
$$
x_0 = -b_i+x_{i+n}+\sum\limits_{j=1}^{n}a_{i,j}x_j
$$
由于我们一定找的时限制最小的,那么$b_i<0\Longrightarrow-b_i >0$,这组约束合法。

对于其余的约束,我们考虑代入之后,会有:
$$
x_{n+k}= b_k-b_i+\sum\limits_{j=1}^{n}(a_{k,j}-a_{i,j})x_j+x_{i+n}
$$
由于$b_i$最小,那么$b_k\geq b_i \Longrightarrow b_k-b_i\geq 0$。任务完成。

那么就可以愉快地上代码辣(≧▽≦)/!题目来自于$\rm{UOJ179}$

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
49
50
51
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>

#define MAXN 1010

using namespace std ;
const double eps = 1e-8, INF = 1e15 ;
int N, M, T, Id[MAXN], i, j ; double A[MAXN][MAXN], Get[MAXN] ;

void Pivot(int l, int e){
swap(Id[N + l], Id[e]) ;//交换两个变量
double t = A[l][e] ; A[l][e] = 1.0 ;
for(i = 0 ; i <= N ; ++ i) A[l][i] /= t ;//类似于高斯消元的第一步,把这一项挪到等式右边,所以要先都除以系数
for(i = 0 ; i <= M ; ++ i)
if(i != l && fabs(A[i][e]) > eps)
for(t = A[i][e], A[i][e] = j = 0 ; j <= N ; ++ j)
A[i][j] -= A[l][j] * t ;//此处则是代换。而由于代数式中肯定不会出现我们的A_{i,e}(毕竟是换到了等式右边),所以对于这一项的系数要特殊处理成0
}
bool Init(){
while(true){
int e = 0, l = 0 ; double t = -eps ;
for (i = 1 ; i <= M ; ++ i) if (A[i][0] < t) l = i, t = A[i][0] ; if(!l) return 1 ;//找系数最小的负值项
for (j = 1 ; j <= N ; ++ j) if (A[l][j] < -eps && (!e || (rand() & 1))) e = j ;
if (!e) return puts("Infeasible"), false ; Pivot(l, e) ;
}//不合法的线性规划,因为我们假设全部的系数都为正,且我们前面已经假定了全部x_i为正,那么不可能满足标准型里面的小于b_l——此时b_l为负值。
return 1 ;
}
bool Simplex(){
while(1){
int l = 0, e = 0 ; double MAX = INF ;
for(j = 1 ; j <= N ; ++ j) if (A[0][j] > eps) { e = j ; break ; } if(!e) break ;
for(i = 1 ; i <= M ; ++ i) if (A[i][e] > eps && A[i][0] / A[i][e] < MAX) MAX = A[i][0] / A[i][e], l = i ;//选个最紧的约束
if (!l) return puts("Unbounded"), false ; Pivot(l,e) ;
}
return 1 ;
}
int main(){
cin >> N >> M >> T ; srand(19260817) ;
for (i = 1 ; i <= N ; ++ i) Id[i] = i ;
for (i = 1 ; i <= N ; ++ i) cin >> A[0][i] ;
for (i = 1 ; i <= M ; cin >> A[i ++][0])
for (j = 1 ; j <= N ; ++ j) cin >> A[i][j] ;
if (!(Init() && Simplex())) return 0 ; printf("%.8lf\n", -A[0][0]) ;
if(T){
for(i = 1 ; i <= M ; ++ i) Get[Id[N + i]] = A[i][0] ;
for(i = 1 ; i <= N ; ++ i) printf("%.8lf ",Get[i]) ;
}
return 0 ;
}

然鹅事实上这份代码只有$97pts$…好像剩下3分就几乎没有人得过,除了kcz和std……

至于中间的srand,是由于我们随便$pivot$一个$e$就好,于是就听长者的

$\rm{0x04}\quad Afterwords$

一直想学,一直没有机会学。

说起来确实有点儿小激动。第一次写单纯型的代码是前不久的一次周末,早上点起来之后头昏脑涨地扑向Candy?的博客,学了起来。在那之后才发现原来线性规划并不常考,甚至说,不考。但是我挺喜欢这种感觉,感觉自己仿佛比出题人会的还多

$\rm{Reference}$

-------------本文结束感谢您的阅读-------------