單純形算法

定義

單純形法是解決線性規劃問題的一個有效的算法。

線性規劃就是在一組線性約束條件下，求解目標函數最優解的問題。

線性規劃的一般形式

在約束條件下，尋找目標函數 \(z\) 的最大值：

\[ \max \ z = x_1 + x_2 \]

\[ s.t \begin{cases} 2x_1 + x_2 \leq 12 \\ x_1 + 2x_2 \leq 9 \\ x_1, x_2 \geq 0 \end{cases} \]

線性規劃的可行域

滿足線性規劃問題約束條件的所有點組成的集合就是線性規劃的可行域。若可行域有界（以下主要考慮有界可行域），線性規劃問題的目標函數最優解必然在可行域的頂點上達到最優。

單純形法就是通過設置不同的基向量，經過矩陣的線性變換，求得基可行解（可行域頂點），並判斷該解是否最優，否則繼續設置另一組基向量，重複執行以上步驟，直到找到最優解。所以，單純形法的求解過程是一個循環迭代的過程。

kexingyu

圖1 可行域

線性規劃的標準形式

在説明單純形法的原理之前，需要明白線性規劃的標準形式。因為單純形算法是通過線性規劃的標準形來求解的。一般，規定線性規劃的標準形式為：

\[ \max \ z = \sum_{j = 1}^{n}c_jx_j \]

\[ s.t \begin{cases} \displaystyle \sum_{j = 1}^{n}a_{ij}x_j = b_i, i = 1,2,\dots,m\\ x_j \geq 0 , j = 1,2,\dots,n \\ \end{cases} \]

寫成矩陣形式：

\[ \max \ z = CX \]

\[ AX = b \]

\[ X \geq 0 \]

\[ A = \begin{bmatrix} a_{11} & a_{12} & \ldots & a_{1n}\\ a_{21} & a_{22} & \ldots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \ldots & a_{mn} \end{bmatrix} \]

標準形的形式為：

目標函數要求 \(\max\)
約束條件均為等式
決策變量為非負約束

普通線性規劃化為標準形：

若目標函數為最小化，可以通過取負，求最大化
約束不等式為小於等於不等式，可以在左端加入非負變量，轉變為等式，比如：

\[ x_1 + 2x_2 \leq 9 \implies \begin{cases} x_1 + 2x_2 + x_3 = 9 \\ x_3 \geq 0 \end{cases} \]

同理，約束不等式為大於等於不等式時，可以在左端減去一個非負鬆弛變量，變為等式。
若存在取值無約束的變量，可轉變為兩個非負變量的差，比如：

\[ -\infty \leq x_k \leq +\infty \implies \begin{cases} x_k = x_m - x_n \\ x_m,x_n \geq 0 \end{cases} \]

本文最開始的線性規劃問題轉化為標準形為：

\[ \max \ z = x_1 + x_2 \]

\[ s.t \begin{cases} 2x_1 + x_2 + x_3 = 12 \\ x_1 + 2x_2 + x_4 = 9 \\ x_1, x_2, x_3, x_4 \geq 0 \end{cases} \]

單純形法

幾何意義

在標準形中，有 \(m\) 個約束條件（不包括非負約束），\(n\) 個決策變量，且 \(n \geq m\)。首先選取 \(m\) 個基變量 \(x_j^{'}(j = 1, 2, \ldots, m )\)，基變量對應約束係數矩陣的列向量線性無關。通過矩陣的線性變換，基變量可由非基變量表示：

\[ x_i^{'} = C_i + \sum_{j = m + 1}^{n}m_{ij}x_j^{'}(i = 1, 2, \ldots , m) \]

如果令非基變量等於 \(0\)，可求得基變量的值：

\[ x_i^{'} = C_i \]

如果為可行解的話，\(C_i\) 大於 \(0\)。那麼它的幾何意義是什麼呢？

還是通過上述具體的線性規劃問題來説明：

\[ \max \ z = x_1 + x_2 \]

\[ s.t \begin{cases} 2x_1 + x_2 + x_3 = 12 \\ x_1 + 2x_2 + x_4 = 9 \\ x_1, x_2, x_3, x_4 \geq 0 \end{cases} \]

如果選擇 \(x_2\)、\(x_3\) 為基變量，那麼令 \(x_1\)、\(x_4\) 等於 \(0\)，可以去求解基變量 \(x_2\)、\(x_3\) 的值。對係數矩陣做行變換，如下所示，\(x_2=\dfrac{9}{2}\)，\(x_3=\dfrac{15}{2}\)。

\[ \begin{bmatrix} {\mathrm{X}} & {x_{1}} & {x_{2}} & {x_{3}} & {x_{4}} & {b} \\ {} & {2} & {1} & {1} & {0} & {12} \\ {} & {1} & {2} & {0} & {1} & {9} \\ {\mathrm{C}} & {1} & {1} & {0} & {0} & {z} \end{bmatrix}\to\begin{bmatrix} {\mathrm{X}} & {x_{1}} & {x_{2}} & {x_{3}} & {x_{4}} & {b} \\ {} & {\frac{3}{2}} & {0} & {1} & {-\frac{1}{2}} & {\frac{15}{2}} \\ {} & {\frac{1}{2}} & {1} & {0} & {\frac{1}{2}} & {\frac{9}{2}} \\ {\mathrm{C}} & {\frac{1}{2}} & {0} & {0} & {-\frac{1}{2}} & {z-\frac{9}{2}} \end{bmatrix} \]

\(X_1=0\) 表示可行解在 \(y\) 軸上；\(X_4=0\) 表示可行解在 \(x_1+2x_2=9\) 的直線上。那麼，求得的可行解即表示這兩條直線的交點，也是可行域的頂點，如圖所示：

kexingyu_point

圖2

所以，通過選擇不同的基變量，可以獲得不同的可行域的頂點。

如何判斷最優

如前所述，基變量可由非基變量表示：

\[ x_i^{'} = C_i + \sum_{j = m + 1}^{n}m_{ij}x_j^{'}(i = 1, 2, \ldots , m) \]

目標函數 \(z\) 也可以完全由非基變量表示：

\[ z = z_0 + \sum_{j = m + 1}^{n} \sigma_j x_j^{'} \]

當達到最優解時，所有的 \(\sigma_j\) 應小於等於 \(0\)，當存在 \(j\)，\(\sigma_j > 0\) 時，當前解不是最優解，為什麼？

當前的目標函數值為 \(z_0\)，其中所有的非基變量值均取 \(0\)。由之前分析可知，\(x_j^{'} = 0\) 代表可行域的某個邊界，是 \(x_j^{'}\) 的最小值。如果可行解逐步離開這個邊界，\(x_j^{'}\) 會變大，因為 \(\sigma_j > 0\)，顯然目標函數的取值也會變大，所以當前解不是最優解。我們需要尋找新的基變量。

如何選擇新的基變量

如果存在多個 \(\sigma_j > 0\)，選擇最大的 \(\sigma_j > 0\) 對應的變量作為基變量，這表示目標函數隨着 \(x_j^{'}\) 的增加，增長的最快。

如何選擇被替換的基變量

假如我們選擇非基變量 \(x_s^{'}\) 作為下一輪的基變量，那麼被替換基變量 \(x_j^{'}\) 在下一輪中作為非基變量，等於 \(0\)。選擇 \(x_j^{'}\) 的原則：替換後應該儘量使 \(x_s^{'}\) 值最大（因為上面已分析過，目標函數會隨着 \(x_s^{'}\) 的增大而增大），但要保證替換基變量後的解仍是可行解，因此應該選擇最緊的限制。

繼續通過上面的例子來説明：

\[ \begin{bmatrix} {\mathrm{X}} & {x_{1}} & {x_{2}} & {x_{3}} & {x_{4}} & {b} \\ {} & {2} & {1} & {1} & {0} & {12} \\ {} & {1} & {2} & {0} & {1} & {9} \\ {\mathrm{C}} & {1} & {1} & {0} & {0} & {z} \end{bmatrix}\to\begin{bmatrix} {\mathrm{X}} & {x_{1}} & {x_{2}} & {x_{3}} & {x_{4}} & {b} \\ {} & {\frac{3}{2}} & {0} & {1} & {-\frac{1}{2}} & {\frac{15}{2}} \\ {} & {\frac{1}{2}} & {1} & {0} & {\frac{1}{2}} & {\frac{9}{2}} \\ {\mathrm{C}} & {\frac{1}{2}} & {0} & {0} & {-\frac{1}{2}} & {z-\frac{9}{2}} \end{bmatrix} \]

從最後一行可以看到，\(x_1\) 的係數為 \(\dfrac{1}{2}>0\)，所以選 \(x_2\)、\(x_3\) 為基變量並沒有是目標函數達到最優。下一輪選取 \(x_1\) 作為基變量，替換 \(x_2\)、\(x_3\) 中的某個變量。

第一行是符號

第二行：若 \(x_1\) 替換 \(x_3\) 作為基變量，\(x_3=0\) 時，\(x_1=\dfrac{15/2}{3/2}=5\)

第三行：若 \(x_1\) 替換 \(x_2\) 作為基變量，\(x_2=0\) 時，\(x_1=\dfrac{9/2}{1/2}=9\)

儘管替換 \(x_2\) 後，\(x_1\) 的值更大，但將它代入 \(x_3\) 後會發現 \(x_3\) 的值為負，不滿足約束。從幾何的角度來看，選擇 \(x_2\) 和 \(x_4\) 作為非基變量，得到的解是直線 \(x_2=0\) 和 \(x_1 + 2x_2 = 9\) 的交點，它在可行域外。因此應該選擇 \(x_3\) 作為非基變量。

終止條件

當目標函數用非基變量的線性組合表示時，所有的係數均不大於 \(0\)，則表示目標函數達到最優。

如果，有一個非基變量的係數為 \(0\)，其他的均小於 \(0\)，表示目標函數的最優解有無窮多個。這是因為目標函數的梯度與某一邊界正交，在這個邊界上，目標函數的取值均相等，且為最優。

使用單純形法來求解線性規劃，輸入單純形法的鬆弛形式，是一個大矩陣，第一行為目標函數的係數，且最後一個數字為當前軸值下的 \(z\) 值。下面每一行代表一個約束，數字代表係數每行最後一個數字代表 \(b\) 值。

算法和使用單純性表求解線性規劃相同。

對於線性規劃問題：

\[ \max \ x_1 + 14x_2 + 6x_3 \]

\[ s.t \begin{cases} x_1 + x_2 + x_3 \leq 4 \\ x_1 \leq 2 \\ x_3 \leq 3 \\ 3x_2 + x_3 \leq 6 \\ x_1, x_2, x_3 \geq 0 \end{cases} \]

我們可以得到其鬆弛形式：

\[ \max \ x_1 + 14x_2 + 6x_3 \]

\[ s.t \begin{cases} x_1 + x_2 + x_3 + x_4 = 4 \\ x_1 + x_5 = 2 \\ x_3 + x_6 = 3 \\ 3x_2 + x_3 + x_7 = 6 \\ x_1, x_2, x_3, x_4, x_5, x_6, x_7 \geq 0 \end{cases} \]

我們可以構造單純形表，其中最後一行打星的列為軸值。

\(x_1\)	\(x_2\)	\(x_3\)	\(x_4\)	\(x_5\)	\(x_6\)	\(x_7\)	\(b\)
\(c_1=1\)	\(c_2=14\)	\(c_3=6\)	\(c_4=0\)	\(c_5=0\)	\(c_6=0\)	\(c_7=0\)	\(-z=0\)
\(1\)	\(1\)	\(1\)	\(1\)	\(0\)	\(0\)	\(0\)	\(4\)
\(1\)	\(0\)	\(0\)	\(0\)	\(1\)	\(0\)	\(0\)	\(2\)
\(0\)	\(0\)	\(1\)	\(0\)	\(0\)	\(1\)	\(0\)	\(3\)
\(0\)	\(3\)	\(1\)	\(0\)	\(0\)	\(0\)	\(1\)	\(6\)
			\(*\)	\(*\)	\(*\)	\(*\)

在單純形表中，我們發現非軸值的 \(x\) 上的係數大於零，因此可以通過增加這些個 \(x\) 的值，來使目標函數增加。我們可以貪心的選擇最大的 \(c\)，在上面的例子中我們選擇 \(c_2\) 作為新的軸，加入軸集合中，那麼誰該出軸呢？

其實我們由於每個 \(x\) 都大於零，對於 \(x_2\) 它的增加是有所限制的，如果 \(x_2\) 過大，由於其他的限制條件，就會使得其他的 \(x\) 小於零，於是我們應該讓 \(x_2\) 一直增大，直到有一個其他的 \(x\) 剛好等於 \(0\) 為止，那麼這個 \(x\) 就被換出軸。

我們可以發現，對於約束方程 \(1\)，即第一行約束，\(x_2\) 最大可以為 \(4\)（\(\dfrac{4}{1}\)），對於約束方程 \(4\)，\(x_2\) 最大可以為 \(2\)（\(\dfrac{6}{3}\)），因此 \(x_2\) 最大隻能為他們之間最小的那個，這樣才能保證每個 \(x\) 都大於零。因此使用第 \(4\) 行，來對各行進行高斯行變換，使得第二列第四行中的每個 \(x\) 都變成零，也包括 \(c_2\)。這樣我們就完成了把 \(x_2\) 入軸，\(x_7\) 出軸的過程。變換後的單純形表為：

\(x_1\)	\(x_2\)	\(x_3\)	\(x_4\)	\(x_5\)	\(x_6\)	\(x_7\)	\(b\)
\(c_1=1\)	\(c_2=0\)	\(c_3=1.33\)	\(c_4=0\)	\(c_5=0\)	\(c_6=0\)	\(c_7=-4.67\)	\(-z=-28\)
\(1\)	\(0\)	\(0.67\)	\(1\)	\(0\)	\(0\)	\(-0.33\)	\(2\)
\(1\)	\(0\)	\(0\)	\(0\)	\(1\)	\(0\)	\(0\)	\(2\)
\(0\)	\(0\)	\(1\)	\(0\)	\(0\)	\(1\)	\(0\)	\(3\)
\(0\)	\(1\)	\(0.33\)	\(0\)	\(0\)	\(0\)	\(0.33\)	\(2\)
	\(*\)		\(*\)	\(*\)	\(*\)

繼續計算，我們得到：

\(x_1\)	\(x_2\)	\(x_3\)	\(x_4\)	\(x_5\)	\(x_6\)	\(x_7\)	\(b\)
\(c_1=-1\)	\(c_2=0\)	\(c_3=0\)	\(c_4=-2\)	\(c_5=0\)	\(c_6=0\)	\(c_7=-4\)	\(-z=-32\)
\(1.5\)	\(0\)	\(1\)	\(1.5\)	\(0\)	\(0\)	\(-0.5\)	\(3\)
\(1\)	\(0\)	\(0\)	\(0\)	\(1\)	\(0\)	\(0\)	\(2\)
\(-1.5\)	\(0\)	\(0\)	\(-1.5\)	\(0\)	\(1\)	\(0.5\)	\(0\)
\(-0.5\)	\(1\)	\(0\)	\(-0.5\)	\(0\)	\(0\)	\(0.5\)	\(1\)
	\(*\)	\(*\)		\(*\)	\(*\)

此時我們發現，所有非軸的 \(x\) 的係數全部小於零，即增大任何非軸的 \(x\) 值並不能使得目標函數最大，從而得到最優解 \(32\)。

整個過程代碼如下所示：

#include <bits/stdc++.h>
using namespace std;
vector<vector<double> > Matrix;
double Z;
set<int> P;
size_t cn, bn;

bool Pivot(pair<size_t, size_t> &p) {  // 返回0表示所有的非軸元素都小於0
  int x = 0, y = 0;
  double cmax = -INT_MAX;
  vector<double> C = Matrix[0];
  vector<double> B;

  for (size_t i = 0; i < bn; i++) {
    B.push_back(Matrix[i][cn - 1]);
  }

  for (size_t i = 0; i < C.size(); i++) {  // 在非軸元素中找最大的c
    if (cmax < C[i] && P.find(i) == P.end()) {
      cmax = C[i];
      y = i;
    }
  }
  if (cmax < 0) {
    return 0;
  }

  double bmin = INT_MAX;
  for (size_t i = 1; i < bn; i++) {
    double tmp = B[i] / Matrix[i][y];
    if (Matrix[i][y] != 0 && bmin > tmp && tmp > 0) {
      bmin = tmp;
      x = i;
    }
  }

  p = make_pair(x, y);

  for (set<int>::iterator it = P.begin(); it != P.end(); it++) {
    if (Matrix[x][*it] != 0) {
      // cout<<"erase "<<*it<<endl;
      P.erase(*it);
      break;
    }
  }
  P.insert(y);
  // cout<<"add "<<y<<endl;
  return true;
}

void pnt() {
  for (size_t i = 0; i < Matrix.size(); i++) {
    for (size_t j = 0; j < Matrix[0].size(); j++) {
      cout << Matrix[i][j] << "\t";
    }
    cout << endl;
  }
  cout << "result z:" << -Matrix[0][cn - 1] << endl;
}

void Gaussian(pair<size_t, size_t> p) {  // 行變換
  size_t x = p.first;
  size_t y = p.second;
  double norm = Matrix[x][y];
  for (size_t i = 0; i < cn; i++) {  // 主行歸一化
    Matrix[x][i] /= norm;
  }
  for (size_t i = 0; i < bn; i++) {
    if (i != x && Matrix[i][y] != 0) {
      double tmpnorm = Matrix[i][y];
      for (size_t j = 0; j < cn; j++) {
        Matrix[i][j] = Matrix[i][j] - tmpnorm * Matrix[x][j];
      }
    }
  }
}

void solve() {
  pair<size_t, size_t> t;
  while (1) {
    pnt();
    if (Pivot(t) == 0) {
      return;
    }
    cout << t.first << " " << t.second << endl;
    for (set<int>::iterator it = P.begin(); it != P.end(); it++) {
      cout << *it << " ";
    }
    cout << endl;
    Gaussian(t);
  }
}

int main(int argc, char *argv[]) {
  // ifstream fin;
  // fin.open("./test");
  cin >> cn >> bn;
  for (size_t i = 0; i < bn; i++) {
    vector<double> vectmp;
    for (size_t j = 0; j < cn; j++) {
      double tmp = 0;
      cin >> tmp;
      vectmp.push_back(tmp);
    }
    Matrix.push_back(vectmp);
  }

  for (size_t i = 0; i < bn - 1; i++) {
    P.insert(cn - i - 2);
  }
  solve();
}

/////////////////////////////////////
// glpk input:
///* Variables */
// var x1 >= 0;
// var x2 >= 0;
// var x3 >= 0;
///* Object function */
// maximize z: x1 + 14*x2 + 6*x3;
///* Constrains */
// s.t. con1: x1 + x2 + x3 <= 4;
// s.t. con2: x1  <= 2;
// s.t. con3: x3  <= 3;
// s.t. con4: 3*x2 + x3  <= 6;
// end;
/////////////////////////////////////
// myinput:
/*
8 5
1 14 6 0 0 0 0 0
1 1 1 1 0 0 0 4
1 0 0 0 1 0 0 2
0 0 1 0 0 1 0 3
0 3 1 0 0 0 1 6
*/
/////////////////////////////////////

理論羅列

標準型

\(m+n\) 個約束 \(n\) 個變量用 \(x\) 向量表示，\(A\) 是一個 \(m\times n\) 的矩陣，\(c\) 是一個 \(n\) 的向量，\(b\) 是一個 \(m\) 的向量，最大化 \(cx\) 滿足約束 \(Ax \leq b,x > 0\)。

最大化 \(\sum_{j=1}^nc_jx_j\) 滿足如下約束條件：

\[ \sum_{j = 1}^na_{ij}x_j \leq b_i,i = 1,2,\ldots,m \]

\[ x_j \geq 0,j = 1,2,\ldots,n \]

\(n\) 個變量，\(m+n\) 個約束，構造 \(m \times n\) 的矩陣 \(A\)，\(m\) 維向量 \(b\)，\(n\) 維向量 \(c\)

最大化 \(C^Tx\) 滿足如下約束條件：

\[ Ax \leq b \]

\[ x \geq 0 \]

轉換為標準型

若目標函數要求取最小值，那麼可以對其取相反數變成取最大值。對於限制條件 \(f(x_1, x_2, \ldots ,x_n) = b\)，可以用兩個不等式 \(f(x_1, x_2, \ldots, x_n) \leq b,-f(x_1,x_2,\ldots,x_n) \leq -b\) 描述，對於限制條件 \(f(x_1,x_2,\ldots,x_n) \geq b\)，可以用不等式 \(-f(x_1,x_2,\ldots,x_n) \leq -b\) 描述。對於無限制的變量 \(x\)，可以將其拆為兩個非負變量 \(x_0,x_1\)，使得 \(x = x_0 - x_1\)。

鬆弛型

基本變量 \(B\)，\(|B|=m\)，一個約束對應一個，表示鬆弛量，叫做鬆弛變量（基本變量）

非基變量 \(N\)，\(|N|=n\)，\(x_n + i = b_i - \sum a_{ij}x_j \geq 0\)

鬆弛變量 \(x_{n+i}\)

\[ \sum_{j = 1}^na_{ij}x_j \leq b_i \rightarrow x_{n + i} = b_i - \sum_{j = 1}^{n}a_{ij}x_j, x_{n+i} \geq 0 \]

等式左側為基本變量，右側為非基本變量。

變量

替入變量 \(x_e\)（非基變量）
替出變量 \(x_l\)（基本變量）

可行解

基本解：所有非基變量設為 \(0\)，基本變量為右側的常數
基本可行解：所有 \(b_i \geq 0\)

注：單純形法的過程中 \(B\) 和 \(N\) 不斷交換，在 \(n\) 維空間中不斷走，「相當於不等式上的高斯消元」。

轉軸

選取一個非基本變量 \(x_e\) 為替入變量，基本變量 \(x_l\) 為替出變量，將其互換，為了防止循環，根據 Bland 規則，選擇下標最小的變量。

Bland 規則 可以參看：最優化方法

初始化

在所有 \(b_i < 0\) 的約束中隨機選一個作為 \(x_l\)，再隨機選一個 \(a_{le} < 0\) 作為 \(x_e\)，然後 \(pivot(l,e)\) 後 \(b_i\) 就變正了。

算法實現

每個約束定義了 \(n\) 維空間中的一個半空間（超平面），交集形成的可行域是一個凸區域稱為單純形。目標函數是一個超平面，最優解在凸區域定點處取得。通過不斷的轉軸操作，在 \(n\) 維凸區域的頂點上不斷移動（轉軸），使得基本解的目標值不斷變大，最終達到最優解。

以下問題可以轉換為單純形：

最短路

最大流

最小費用最大流

多商品流

基本思想就是改寫 \(l\) 這個約束為 \(x_e\) 作為基本變量，然後把這個新 \(x_e\) 的值帶到其他約束和目標函數中，就消去 \(x_e\) 了。改寫和帶入時要修改 \(b\) 和 \(a\)，目標函數則是 \(c\) 和 \(v\)。

轉動時，\(l\) 和 \(e\) 並沒有像算法導論上一樣，\(a\) 矩陣用了兩行分別是 \(a_{l, \square}\) 和 \(a_{e, \square}\)（這樣佔用內存大），而是用了同一行，這樣 \(a\) 矩陣的行數 \(=|B|\)，列數 \(=|N|\)。

也就是説，約束條件只用 \(m\) 個，儘管 \(B\) 和 \(N\) 不斷交換，但同一時間還是隻有 \(m\) 個約束（基本變量），\(n\) 個非基變量，注意改寫成鬆弛型後 \(a\) 矩陣實際係數為負。（一個優化為 \(a_{i,e}\) 的約束沒必要帶入了）。

simplex 是主過程，基本思想是找到一個 \(c_e>0\) 的，然後找對這個 \(e\) 限制最緊的 \(l\)，轉動這組 \(l,e\)，注意精度控制 \(\epsilon\)，\(c_e>\epsilon\)，還有找 \(l\) 的時候 \(a_{i,e}>\epsilon\) 才行。

例題「NOI2008」志願者招募

題目大意：長度為 \(n\) 的序列，第 \(i\) 位至少 \(b_i\)，\(m\) 種區間使 \([l_i,r_i] + 1\) 代價為 \(a_i\)。

原始問題 \(m\) 個變量，\(n\) 個約束，當 \(l_j \leq i \leq r_j\)，\(a_{ij} = 1\)。

對偶問題 \(n\) 個變量，\(m\) 個約束

\[ \max \ \sum_{i=1}^n b_iy_i \]

\[ s.t \ \sum_{l_i \leq j \leq r_i}y_j \leq c_i, y_i \geq 0 \]

把對應出的係數矩陣代入到單純形算法就可以求出最優解了。

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int M = 10005, N = 1005, INF = 1e9;

int read() {  // 快读
  char c = getchar();
  int x = 0, f = 1;
  while (c < '0' || c > '9') {
    if (c == '-') f = -1;
    c = getchar();
  }
  while (c >= '0' && c <= '9') {
    x = x * 10 + c - '0';
    c = getchar();
  }
  return x * f;
}

int n, m;
double a[M][N], b[M], c[N], v;

void pivot(int l, int e) {  // 转轴操作函数
  b[l] /= a[l][e];
  for (int j = 1; j <= n; j++)
    if (j != e) a[l][j] /= a[l][e];
  a[l][e] = 1 / a[l][e];

  for (int i = 1; i <= m; i++)
    if (i != l && fabs(a[i][e]) > 0) {
      b[i] -= a[i][e] * b[l];
      for (int j = 1; j <= n; j++)
        if (j != e) a[i][j] -= a[i][e] * a[l][j];
      a[i][e] = -a[i][e] * a[l][e];
    }

  v += c[e] * b[l];
  for (int j = 1; j <= n; j++)
    if (j != e) c[j] -= c[e] * a[l][j];
  c[e] = -c[e] * a[l][e];

  // swap(B[l],N[e])
}

double simplex() {
  while (true) {
    int e = 0, l = 0;
    for (e = 1; e <= n; e++)
      if (c[e] > (double)0) break;
    if (e == n + 1) return v;  // 此时v即为最优解
    double mn = INF;
    for (int i = 1; i <= m; i++) {
      if (a[i][e] > (double)0 && mn > b[i] / a[i][e]) {
        mn = b[i] / a[i][e];  // 找对这个e限制最紧的l
        l = i;
      }
    }
    if (mn == INF) return INF;  // unbounded
    pivot(l, e);                // 转动l,e
  }
}

int main() {
  n = read();
  m = read();
  for (int i = 1; i <= n; i++) c[i] = read();
  for (int i = 1; i <= m; i++) {
    int s = read(), t = read();
    for (int j = s; j <= t; j++) a[i][j] = 1;  // 表示第i种志愿者在j时间可以服务
    b[i] = read();
  }
  printf("%d", (int)(simplex() + 0.5));
}

對偶原理

最大化與最小化互換，常數與目標函數互換，改變不等號，變量與約束對應。

\[ \max \ c^Tx: Ax \leq b, x \geq 0 \]

\[ \min \ b^Ty: A^Ty \geq c, y \geq 0 \]

\(d_{uv}\) 表示 \(u,v\) 是否匹配

\[ \max \ \sum_{(u,v) \in E}c_{uv}d_{uv} \]

\[ s.t \begin{cases} \sum_{(v) \in Y} d_{uv} \leq 1, u \in X \\ \sum_{(u) \in X} d_{uv} \leq 1, v \in Y \\ d_{u,v} \in \{0,1\} \end{cases} \]

令 \(p_u,p_v\) 為兩類約束對偶之後的變量

\[ \min \ \sum_{u \in x} p_u + \sum_{v \in Y} p_v \]

\[ s.t \begin{cases} p_u + p_v \geq c_{uv} \\ u \in X, v \in Y\\ p_u, p_v \geq 0 \end{cases} \]

全幺模矩陣（Totally Unimodular Matrix）

充分條件：

僅有 \(-1,0,1\) 構成
每列至多兩個非零數
行可分為兩個集合：
- 一列包含兩個同號非零數，兩行不在同一個集合
- 一列包含兩個異號非零數，兩行在同一個集合

線性規劃中 \(A\) 為全幺模矩陣，則單純形法過程中所有係數 \(\in -1,0,1\)，可以去除係數為 \(0\) 的項進行優化！

注：任何最大流、最小費用最大流的線性規劃都是全幺模矩陣

更多詳細的解釋參看：https://www.cnblogs.com/ECJTUACM-873284962/p/7097864.html

習題練習

UOJ#179. 線性規劃

參考資料

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用