张天昀的个人博客

单纯形算法模板

（2019年03月01日）

算法介绍

单纯形算法是问题求解中的一个（不是那么）重要的算法，用于求解线性规划问题。

代码模板

以下单纯形算法模板可以用于求线性规划问题是否有可行解、是否无解，并求出最优解。

（感谢李顶为同学指出这个模板早期版本中存在的bug。）

使用方法：将A[0][j]A[0][j]设置为目标函数的参数，整个函数为

f=∑j=1nA[0][j]f = \sum\limits_{j = 1}^{n} A[0][j]

第ii个约束条件为

∑j=1nA[i][j]⩽A[i][0]\sum\limits_{j = 1}^n A[i][j] \leqslant A[i][0]

运行simplex()的结果若为-1则无解，1则无界，0则代表有解。最优解为A[0][0]A[0][0]，v[j]v[j]为各个参数的值。

const double EPS = 1e-8;

int id[MAXN+MAXM] = {};
double v[MAXN] = {};
double a[MAXM][MAXN] = {};
 
int sgn(double x) {
  if (x < -EPS) return -1;
  return x > EPS ? 1 : 0;
}
 
void pivot(int r, int c) {
  swap(id[n + r], id[c]);
  double x = -a[r][c];
  a[r][c] = -1;
  for (int i = 0; i <= n; ++i) a[r][i] /= x;
  for (int i = 0; i <= m; ++i) {
    if (sgn(a[i][c]) && i != r) {
      x = a[i][c];
      a[i][c] = 0;
      for (int j = 0; j <= n; ++j) a[i][j] += x * a[r][j];
    }
  }
}
 
int simplex() {
  /* important: revert symbols of conditions */
  for (int i = 1; i <= m; ++i) {
    for (int j = 1; j <= n; ++j) {
      a[i][j] *= -1;
    }
  }
  for (int i = 1; i <= n; ++i) id[i] = i;
  /* initial-simplex */
  while (true) {
    int x = 0, y = 0;
    for (int i = 1; i <= m; ++i) {
      if (sgn(a[i][0]) < 0) { x = i; break; }
    }
    if (!x) break;
    for (int i = 1; i <= n; ++i) {
      if (sgn(a[x][i]) > 0) { y = i; break; }
    }
    if (!y) return -1; // infeasible
    pivot(x, y);
  }
  /* solve-simplex */
  while (true) {
    int x = 0, y = 0;
    for (int i = 1; i <= n; ++i) {
      if (sgn(a[0][i]) > 0) { x = i; break; }
    }
    if (!x) break; // finished
    double w = 0, t = 0; bool f = true;
    for (int i = 1; i <= m; ++i) {
      if (sgn(a[i][x]) < 0) {
        t = -a[i][0] / a[i][x];
        if (f || t < w) {
          w = t, y = i, f = false;
        }
      }
    }
    if (!y) { return 1; } // unbounded
    pivot(y, x);
  }
  for (int i = 1; i <= n; ++i) v[i] = 0;
  for (int i = n + 1; i <= n + m; ++i) v[id[i]] = a[i - n][0];
  return 0;
}

References