线性规划之单纯形法

最后更新于：2023年2月24日下午

在 Codeforces 有人问了我两个优化问题，一个是非线性规划（具体说是凸优化问题）。一般来说非线性规划没有什么具体的算法，但是凸优化，可以转化成凸包，然后转换成线段（如果是二维的）上的最值问题就搞定了。另一个是线性规划（所有线性规划其实也是特殊的凸优化），线性规划有著名的单纯形算法，7 年前就学过，但一直没写过代码。趁这次机会重新学习一下单纯形算法，并给出代码。英文版解答

参考教材：运筹学第三版清华大学出版社，另外 OI-wiki 上讲的简洁清晰但是不够全面。

Gauss 消元法

之前一直不写这个模板的原因：可能无解，可能唯一解，可能无穷多个解，double 有判断，很烦。

求解 \(Ax = b\)，如果无解就输出空向量，否则输出（某一个）答案向量，无穷解的话随便输出一个。

浮点数版

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

std::vector<double> Gauss(std::vector<std::vector<double>> A, std::vector<double> b) {
  int n = A.size(), m = A[0].size();
  std::vector<double> x(m);
  std::vector<int> p(m);
  std::iota(p.begin(), p.end(), 0);
  const double eps = 1e-12;
  auto findNonZero = [&](int i) { // 实际上找最大的比较好
    for (int row = i; row < n; ++row) if (fabs(A[row][i]) > eps) return row;
    return n;
  };
  auto triangleGauss = [&](int sz) { // A[i][i] = 1
    std::vector<double> x(sz);
    for (int i = sz - 1; i >=0; --i) {
      x[i] = b[i];
      for (int row = 0; row < i; ++row) b[row] -= A[row][i] * x[i];
    }
    x.resize(A[0].size());
    return x;
  };
  int sz = n;
  for (int i = 0, row; i < n; ++i) {
    while (i < m) {
      row = findNonZero(i);
      if (row != n) break;
      for (int j = 0; j < n; ++j) A[j][i] = A[j][m - 1];
      std::swap(p[i], p[--m]);
    }
    if (i == m) {
      for (int row = m; row < n; ++row) if (fabs(b[row])) {
        // std::cout << "\nNo answer\n";
        return std::vector<double>();
      }
      sz = i;
      break;
    }
    if (row != i) {
      std::swap(A[row], A[i]);
      std::swap(b[row], b[i]);
    }
    b[i] /= A[i][i];
    for (int j = m - 1; j >= i; --j) A[i][j] /= A[i][i];
    for (int row = i + 1; row < n; ++row) {
      b[row] -= A[row][i] * b[i];
      for (int j = m - 1; j >= i; --j) {
        A[row][j] -= A[row][i] * A[i][j];
      }
    }
  }
  // if (sz != A[0].size()) std::cout << "\nInfinite answer\n";
  auto xt = triangleGauss(sz);
  for (int t = 0; t < A[0].size(); ++t) x[p[t]] = xt[t];
  return x;
}

int main() {
  // freopen("in","r",stdin);
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int cas = 1;
  std::cin >> cas;
  while (cas--) {
    int n, m;
    std::cin >> n >> m;
    std::vector<std::vector<double>> a(n, std::vector<double>(m, 0));
    for (auto &x : a) {
      for (auto &i : x) std::cin >> i;
    }
    std::vector<double> b(n);
    for (auto &x : b) std::cin >> x;
    auto x = Gauss(a, b);
    for (auto t : x) std::cout << t << "\n";
  }
  return 0;
}

有限域版

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

std::vector<LL> Gauss(std::vector<std::vector<LL>> A, std::vector<LL> b) {
  int n = A.size(), m = A[0].size();
  std::vector<LL> x(m);
  std::vector<int> p(m);
  std::iota(p.begin(), p.end(), 0);
  const LL M = 998244353;
  std::function<LL(LL)> inv = [&](LL a) -> LL {
    return a == 1 ? 1 : (M - M / a) * inv(M % a) % M;
  };
  auto sub = [](LL &x, LL y) {
    (x -= y) < 0 && (x += M);
  };
  auto triangleGauss = [&](int sz) { // A[i][i] = 1
    std::vector<LL> x(sz);
    for (int i = sz - 1; i >=0; --i) {
      x[i] = (b[i] + M) % M;
      for (int row = 0; row < i; ++row) sub(b[row], A[row][i] * x[i] % M);
    }
    x.resize(A[0].size());
    return x;
  };
  auto findNonZero = [&](int i) {
    for (int row = i; row < n; ++row) if (A[row][i]) return row;
    return n;
  };
  int sz = n;
  for (int i = 0, row; i < n; ++i) {
    while (i < m) {
      row = findNonZero(i);
      if (row != n) break;
      for (int j = 0; j < n; ++j) A[j][i] = A[j][m - 1];
      std::swap(p[i], p[--m]);
    }
    if (i == m) {
      for (int row = m; row < n; ++row) if (b[row]) {
        // std::cout << "\nNo answer\n";
        return std::vector<LL>();
      }
      sz = i;
      break;
    }
    if (row != i) {
      std::swap(A[i], A[row]);
      std::swap(b[i], b[row]);
    }
    LL inva = inv(A[i][i]);
    (b[i] *= inva) %= M;
    for (int j = m - 1; j >= i; --j) (A[i][j] *= inva) %= M;
    for (int row = i + 1; row < n; ++row) {
      sub(b[row], A[row][i] * b[i] % M);
      for (int j = m - 1; j >= i; --j) {
        sub(A[row][j], A[row][i] * A[i][j] % M);
      }
    }
  }
  // if (sz != A[0].size()) std::cout << "\nInfinite answer\n";
  auto xt = triangleGauss(sz);
  for (int t = 0; t < A[0].size(); ++t) x[p[t]] = xt[t];
  return x;
}

int main() {
  // freopen("in","r",stdin);
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int cas = 1;
  std::cin >> cas;
  while (cas--) {
    int n, m;
    std::cin >> n >> m;
    std::vector<std::vector<LL>> a(n, std::vector<LL>(m, 0));
    for (auto &x : a) {
      for (auto &i : x) std::cin >> i;
    }
    std::vector<LL> b(n);
    for (auto &x : b) std::cin >> x;
    auto x = Gauss(a, b);
    for (auto t : x) std::cout << t << "\n";
  }
  return 0;
}

上面做法是先化成上三角再求，其实也可以先直接化成对角的，两种都挺好，就不过不改了。

单纯形法

首先无论是什么样的线性规划问题，都先化成标准形式：\(\max z = \sum c_i x_i\)，其中 \(Ax = b, x_i \geq 0\)（这里输入的 \(b\) 必然都是非负，否则显然无可行解），并且保证 \(A\) 的左边是一个单位阵。简单的说通过 “大 M 法” 使得 (\(b_1, \cdots, b_n, 0, \cdots, 0\)) 是可以可行解。

求极小值通过 \(c\) 取负号解决
\(\geq b\) 通过加一个变量变成等号
\(\leq b\) 通过减一个变量变成等号

#include <bits/stdc++.h>
#define watch(x) std::cout << (#x) << " is " << (x) << std::endl
using LL = long long;
using pii = std::pair<int, int>;
using pll = std::pair<LL, LL>;

using VD = std::vector<double>;
const double eps = 1e-10;
const double inf = 1e10;
// make sure that A = (I, A') and b >= 0, compute max cx
VD simplexCore(VD c, std::vector<VD> A, VD b) {
  int n = A.size(), m = c.size();
  std::vector<int> p(m);
  std::iota(p.begin(), p.end(), 0);
  for (int i = 0; i < n; ++i) A[i].emplace_back(b[i]);
  c.emplace_back(0);
  A.emplace_back(c);
  for (int j = n; j <= m; ++j) {
    for (int i = 0; i < n; ++i) {
      A[n][j] -= A[n][i] * A[i][j];
    }
  }
  auto check = [&]() -> bool {
    for (int j = n; j < m; ++j) if (A[n][j] > eps) {
      bool flag = false;
      for (int i = 0; i < n; ++i) if (A[i][j] > eps) {
        flag = true;
        break;
      }
      if (!flag) return false;
    }
    return true;
  };
  while (1) {
    int ch = std::max_element(A[n].begin() + n, A[n].begin() + m) - A[n].begin(), hc;
    if (A[n][ch] < eps) break;
    assert(check()); // otherwise unbounded, no max solution
    double theta = DBL_MAX;
    for (int i = 0; i < n; ++i) if (A[i][ch] > eps && A[i].back() / A[i][ch] < theta) {
      theta = A[i].back() / A[i][ch];
      hc = i;
    }
    std::swap(p[ch], p[hc]);
    double tmp = 1 / A[hc][ch];
    for (int j = n; j <= m; ++j) A[hc][j] *= tmp;
    for (int i = 0; i <= n; ++i) if (i != hc) {
      for (int j = n; j <= m; ++j) if (j != ch) {
        A[i][j] -= A[i][ch] * A[hc][j];
      }
    }
    for (int i = 0; i <= n; ++i) A[i][ch] *= -tmp;
    A[hc][ch] = tmp;
  }
  VD x(m);
  for (int i = 0; i < n; ++i) x[p[i]] = A[i].back();
  watch(-A.back().back()); // max_val
  return x; // point Corresponds to max_val
}
// compute max cx, with Aqx = bq and Alq x <= blq, end of 0 can be ommit in A and Aq
VD simplex(VD c, std::vector<VD> Aq, VD bq, std::vector<VD> Alq, VD blq) {
  assert(Aq.size() == bq.size());
  assert(Alq.size() == blq.size());
  int n = Aq.size() + Alq.size();
  int m = c.size();
  for (int i = 0; i < bq.size(); ++i) if (bq[i] < -eps) {
    for (auto &x : Aq[i]) x = -x;
    bq[i] = -bq[i];
  }
  for (int i = 0; i < blq.size(); ++i) if (blq[i] < -eps) {
    for (auto &x : Alq[i]) x = -x;
    ++m;
  }
  std::vector<VD> A(n, VD(n + m));
  VD f(n + m), b(n);
  int now = n + c.size();
  for (int i = 0; i < n; ++i) A[i][i] = 1;
  for (int i = 0; i < Aq.size(); ++i) {
    for (int j = 0; j < Aq[i].size(); ++j) A[i][n + j] = Aq[i][j];
    b[i] = bq[i];
    f[i] = -inf;
  }
  for (int i = 0; i < Alq.size(); ++i) {
    for (int j = 0; j < Alq[i].size(); ++j) A[i + Aq.size()][n + j] = Alq[i][j];
    if (blq[i] < -eps) {
      A[i + Aq.size()][now++] = -1;
      f[i + Aq.size()] = -inf;
    }
    b[i + Aq.size()] = fabs(blq[i]);
  }
  for (int i = 0; i < c.size(); ++i) f[n + i] = c[i];
  auto x = simplexCore(f, A, b);
  return VD(x.begin() + n, x.begin() + n + c.size());
}
int main() {
  // freopen("in","r",stdin);
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int cas = 1;
  std::cin >> cas;
  while (cas--) {
    int nlq, nq, m;
    std::cin >> nlq >> nq >> m;
    std::vector<std::vector<double>> Alq(nlq, std::vector<double>(m, 0));
    std::vector<std::vector<double>> Aq(nq, std::vector<double>(m, 0));
    std::vector<double> blq(nlq), bq(nq), c(m);
    for (auto &x : c) std::cin >> x;
    for (auto &x : Alq) for (auto &i : x) std::cin >> i;
    for (auto &x : blq) std::cin >> x;
    for (auto &x : Aq) for (auto &i : x) std::cin >> i;
    for (auto &x : bq) std::cin >> x;
    auto x = simplex(c, Aq, bq, Alq, blq);
    for (auto t : x) std::cout << t << "\n";
  }
  return 0;
}
/* 输入样例：（和 matlab linprog 样例一致）
1
6 1 2
1 0.33333333

1 1
1 0.25
1 -1
-0.25 -1
-1 -1
-1 1
2 1 2 1 -1 2

1 0.25
0.5
*/

SagMath 线性规划文档以及 Matlab linprog 文档

对偶理论

例题：605C，利用对偶理论之后，此题两种做法，一种是利用三分法搞定，当然了要特别注意精度问题，另一种利用半平面的交（可以查看 HDU 模板）。

以下内容以后有空再补吧

整数规划

混合型规划

#algorithm

图论 cpp 模板上一篇