线性规划 单纯型算法

(Updated: )
(7 mins to read)

n个变量$x_1,…,x_n$​，m个线性约束，形如$\sum \limits_{i=1}^na_{ij}x_i \leq b_j$​,要求确定这n个变量的值，$\max z =\sum\limits_{i=1}^n c_ix_i$单纯型可以被卡到指数级，但一般情况下速度较快，相比于MCMF，代码更短，而且更直观以下用矩阵形式表达（规定小写为列向量，大写为矩阵）

标准形式

最大化$c^Tx$，满足约束$Ax\leq b,x_i \geq 0$

松弛形式

最大化$c^Tx$，满足约束$Ax=b,x_i \geq 0$只要对每个标准形式的约束条件增加一个变量$x_{n+i} \geq 0$即可

对偶问题

最小化$b^Tx$，满足约束$A^Tx \leq c,x_i \geq 0$目标函数的系数c与约束条件右边的系数b交换，并将约束系数矩阵转置，最大化与最小化互换，改变不等式号

$a_{00}=-z,a_{0i}=c_i,a_{i0}=b_i$m行n列，表示n个变量，m个约束，单纯型使用的是松弛形式，只不过没有显示表示出来当约束系数$a_{ij}$​为0，1，-1（全幺模矩阵）时用以下代码

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
int n, m;
int a[N][M], q[M];
const int INF = 0x3f3f3f3f;
void cal(int l, int e)
{
  a[l][e] = 1; int p = 0;
  for(int j=0; j<=n; j++) if(a[l][j]) q[++p] = j;
  for(int i=0, t=0; i<=m; i++)
    if(i!=l && (t=a[i][e]))
    {
      a[i][e] = 0;
      for(int j=1; j<=p; j++) a[i][q[j]] -= t*a[l][q[j]];
    }
}
int work()
{
  while(true)
  {
    int l = 0, e = 0, mn = INF;
    for(int j=1; j<=n; j++)
      if(a[0][j]>0) { e = j; break; }
    if(!e) break;
    for(int i=1; i<=m; i++)
      if(a[i][e]>0 && a[i][0]<mn) mn = a[i][0], l = i;
    assert(l);
    cal(l, e);
  }
  return -a[0][0];
}

一般情况下，原问题或者对偶问题中至少有一个零解是可行解，所以不需要init，否则需要先通过init找到一组初始解，再进行simplex优化

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
52
53
54
55
56
57
58
const double eps = 1e-8, INF = 1e15;
int n, m, type, id[M*2], q[M];
double a[N][M], ans[M];
void pivot(int l, int e)
{
  swap(id[l+n], id[e]);
  double t = a[l][e]; a[l][e] = 1.0; int p = 0;
  for(int j=0; j<=n; j++) a[l][j] /= t;
  for(int j=0; j<=n; j++) if(fabs(a[l][j])>eps) q[++p] = j;
  for(int i=0; i<=m; i++)
    if(i!=l && fabs(a[i][e])>eps)
    {
      double t = a[i][e]; a[i][e] = 0.0;
      for(int j=1; j<=p; j++) a[i][q[j]] -= t*a[l][q[j]];
    }
}
bool init()
{
  while(true)
  {
    int l = 0, e = 0;
    for(int i=1; i<=m; i++)
      if(a[i][0]<-eps && (!l || (rand()&1))) l = i;
    if(!l) break;
    for(int j=1; j<=n; j++)
      if(a[l][j]<-eps && (!e || (rand()&1))) e = j;
    if(!e) return puts("Infeasible"), false; //无解
    pivot(l, e);
  }
  return true;
}
bool simplex()
{
  while(true)
  {
    int l = 0, e = 0; double w = eps;
    for(int j=1; j<=n; j++) //选择最大的那个（更快找到最优解），或者选择第一个大于0的（可以防止死循环）
      if(a[0][j]>w) w = a[0][j], e = j;
    if(!e) break;
    w = INF;
    for(int i=1; i<=m; i++)
      if(a[i][e]>eps && a[i][0]/a[i][e]<w)
        w = a[i][0]/a[i][e], l = i;
    if(!l) return puts("Unbounded"), false;//解无界
    pivot(l, e);
  }
  return true;
}
void solve()
{
  iota(id+1, id+n+1, 1);
  if(init() && simplex())
  {
    printf("%.8f\n", -a[0][0]);
    for(int i=1; i<=m; i++) ans[id[n+i]] = a[i][0];
    if(type) for(int i=1; i<=n; i++) printf("%.8f ", ans[i]);
  }
}

线性规划转费用流：列出标准形式，添加松弛变量转化为松弛型每个式子与前一个式子相减，如果每个变量只出现两次，且系数一正一负，就可以转化为费用流对于每个等式新建一个点，如果等式中右边常数b为正数，从源点s向该点连边（b,0），否则该点向汇点连边（-b,0）。对每个变量（包括添加的用来松弛的变量）从它系数为正的等式代表的点向它系数为负的等式代表的点连边（INF，c），c为目标函数中的系数