- 閱讀本文你將會知道
- 線性規劃簡介
- 線性規劃的標準形
- 一般型轉標準型
- < 與 ≤
- 線性規劃的鬆弛形
- 標準型轉鬆弛形
- 單純形演算法
- 基本可行解
- 如何判斷最優
- 旋轉操作
- 如何透過旋轉更新解？
- 退化與布蘭德規則
- 基本不可行解
- 單純形演算法的幾何意義
- 單純形演算法的時間複雜度分析
- 線性規劃問題有更優的做法嗎？
- 對偶定理
- 全么模矩陣
- 例題
- 「UOJ #179」線性規劃
- 「NOI2008」志願者招募
- 「ABC231H」Minimum Coloring
- 「SHOI2004」最小生成樹
- 「CF1430G」Yet Another DAG Problem
- 參考資料
閱讀本文你將會知道
- 線性規劃與單純形演算法
- 單純形演算法時間複雜度分析
- 單純形 C++ 程式碼實現
- 單純形演算法在演算法題目中的運用
線性規劃簡介
首先引入機床廠問題：
某機床廠生產甲、乙兩種機床，每臺銷售後的利潤分別為 \(4000\) 元與 \(3000\) 元。 生產甲機床需用 \(A\)、\(B\) 機器加工，加工時間分別為每臺 \(2\) 小時和 \(1\) 小時；生產乙機床需用 \(A\)、\(B\)、\(C\) 三種機器加工，加工時間為每臺各一小時。若每天可用於加工的機器時數分別為 \(A\) 機器 \(10\) 小時、\(B\) 機器 \(8\) 小時和 \(C\) 機器 \(7\) 小時，問該廠每天應生產甲、乙機床各幾臺，才能使該天的總利潤最大？
上述問題的數學模型：設該廠生產 \(x_1\) 臺甲機床和 \(x_2\) 臺乙機床時總利潤最大，則 \(x_1\)、\(x_2\) 應滿足：
這裡變數 \(x_1\)、\(x_2\) 稱之為決策變數，上面的最大式被稱為問題的目標函式，下面約束條件中的幾個不等式是問題的約束條件，記為 \(\text{s.t.}\)（即 subject to）。由於上面的目標函式及約束條件均為線性函式，故被稱為線性規劃問題。
線性規劃的標準形
線性規劃的標準形式如下：
一般型轉標準型
- 如果題目要求最小值，那麼可以把目標函式的係數均乘上 \(-1\) 轉換為最大值；
- 如果約束為 \(ax=b\)，那麼轉換為 \(ax\le b\) 和 \(ax\ge b\) 兩個約束；
- 如果約束為 \(ax\ge b\)，該約束的係數均乘上 \(-1\) 就可以變成 \(-ax\le-b\)；
- 如果對 \(x\) 的值域沒有要求，那麼可以用 \(x'-x''\) 替代 \(x\)，其中 \(x',x''\ge0\)；
- 其它情況讀者可以自行嘗試。
< 與 ≤
那麼有人可能要問了，為什麼不能將 \(x<1\) 的一般型轉換為標準型呢？
舉個例子：
發現是無解的。
線性規劃的鬆弛形
因為不等式處理起來不如等式方便，所以我們定義鬆弛型：
標準型轉鬆弛形
考慮加入輔助變數 \(x_{n+1}\sim x_{n+m}\)，如下：
那麼就有：
以機床廠問題為例：
可以被轉換成：
其中，\(x_3\sim x_5\) 是輔助變數。我們給它們起一個名字，叫做基本變數。同理，\(x_1\sim x_2\) 就是非基本變數。
也就是說，\(x_{1}\sim x_n\) 為非基本變數，\(x_{n+1}\sim x_{n+m}\) 為基本變數。
單純形演算法
我們接下來就會介紹單純形演算法——一個求解線性規劃的經典方法。
它有三步：
- 找到一個初始的基本可行解；
- 不斷的進行旋轉（pivot）操作；
- 重複 2 直到結果最優。
我們以上面的線性規劃鬆弛型為例子（此後，我們會省略變數均大於等於 \(0\) 這一要求，請讀者自行理解）：
基本可行解
此時我們要進行第一步，找到一個基本可行解。
所謂基本可行解，就是找到一組基本變數的取值，使得其滿足題目要求。
我們假設非基本變數取值都是 \(0\)，以此得到基本變數的取值，然後容易發現，\(x_1=0,x_2=0,x_3=10,x_4=8,x_5=7,-z=0\)。
發現它們都是滿足 \(x_1,x_2,x_3,x_4,x_5\ge0\) 的要求的。所以這是基本可行解。
當然，也存在不合法的初始情況，稱為基本不可行解，這一點我們後面再說。
如何判斷最優
我們知道，目標函式可以由非基本變數得出：
當達到最優解時，我們發現 \(c_j\) 必定都小於等於 \(0\)。
為什麼，因為若 \(c_j\) 大於 \(0\)，那麼 \(z\) 的值則會與 \(x_j\) 成正比例關係，若 \(x_j\) 變大，\(z\) 也會變大，所以當前並不是最優解。
我們發現，當前的 \(z=4x_1+3x_2\)，都是大於 \(0\) 的係數，不是最優值。
所以，當我們發現 \(z\) 的係數全部都小於等於 \(0\) 時，此時答案最優！
如何修改目標函式的係數呢？我們可以透過旋轉操作。
旋轉操作
所謂旋轉操作，就是將非基本變數與基本變數交換的操作。
可能有點拗口，舉個例子：
將 \(x_1\) 與 \(x_3\) 交換，我們就用 \(x_3\) 來表示 \(x_1\)：
然後把 \(x_1\) 當作基本變數，將 \(x_3\) 當作非基本變數，把 \(x_1\) 代入其它式子：
因為交換了，我們把 \(x_3\) 移動到 \(x_1\) 原本那一列，看看和原式子有什麼不同：
比較一下：
我們發現，對於將基本變數 \(x_3\) 與非基本變數 \(x_1\) 交換的操作，原本的 \(x_3\) 這一行（第二行）除去 \(x_1\) 自己，其它的係數都除以了 \(x_1\) 的係數 \(2\)，而 \(x_1\) 自己則是因為與 \(x_3\) 交換變成了 \(x_3\) 的係數（也就是 \(1\)）之後才除以了 \(x_1\) 原本的係數。
仔細觀察其它列的變換，我們發現它們都減去了自己原本 \(x_1\) 係數倍的 \(x_3\) 這一行的值，比如對於第一行的第一項，它減去了原本 \(x_1\) 的係數（也就是 \(4\)）乘以第二行的第一項 \(5\)（注意是更新後的）也就是減去了 \(20\)；比如對於第一行的第二項，它自己先是變成了 \(x_3\) 的係數 \(0\)，再減去了原本 \(x_1\) 的係數（也就是 \(4\)）乘以第二行的第二項 \(\frac{1}{2}\)（注意也是更新後的）也就是減去了 \(2\)。
沒有理解沒有關係，我們換個角度：
我們以矩陣的形式觀察，首先可以將這原本的鬆弛型式子的係數轉換為矩陣（定義矩陣下標從 \(0\) 開始，矩陣右邊多出來的三列就是基本變數的係數，我們將基本變數標紅）：
我們進行換元操作時把 \(x_1\) 與 \(x_3\) 交換，\(x_1\) 變成基本變數，\(x_3\) 變成非基本變數：
同時我們對第二行進行處理以保證基本變數係數均為 \(1\)：
我們此時需要將其他行原本的 \(x_1\) 係數化為 \(0\)，實際上就是將每一項都減去這一行原本 \(x_1\) 係數倍的第二行，同時因為其他行的 \(x_3\) 係數都是 \(0\)，所以需要替換：
得到：
發現實際上三個基本變數都是 \(1\) 或者 \(0\)，所以我們將它們省略，得到一個更簡略的矩陣：
同理的變換方法，交換系數（左邊是包含那三個基本變數的過程，右邊是省略三個基本變數（注意，是基本變數，不一定是 \(x_3,x_4,x_5\)）的過程，對照著看更清晰）：
將其它行的 \(x_1\) 都消掉：
在程式碼實現上，明顯省略三個基本變數會更加好寫。
對於將基本變數 \(x_{n+l}\) 與非基本變數 \(x_e\) 交換，有以下程式碼：
void pivot(int l, int e) {
double t = a[l][e];
a[l][e] = 1; // 變成 x_n+l 的係數 1
for(int i = 0; i <= n; i ++) a[l][i] /= t;
for(int i = 0; i <= m; i ++) if(i != l && abs(a[i][e]) > eps) {
t = a[i][e]; a[i][e] = 0; // 變成 x_n+l 在這一行的係數 0
for(int j = 0; j <= n; j ++) {
a[i][j] -= a[l][j] * t;
}
}
}
如何透過旋轉更新解？
我們的目標是將第一行目標函式的係數都變成小於等於 \(0\) 的，而我們發現將在目標函式中係數大於 \(0\) 的非基本變數 \(x_1\) 與在這列係數同樣大於 \(0\) 的非基本變數 \(x_3\) 交換時，可以使得這一個目標函式的係數變得小於等於 \(0\)。
原理是在交換時 a[i][e] = 0 又會 a[i][j] -= a[l][j] * t，所以當 \(i=0,j=e\) 時，這一個係數將會變得小於等於 \(0\)。
所以我們考慮找到一個 \(a_{0,e}>0(1\le e\le n)\)，如果沒有，那麼就滿足目標函式的係數全部都小於等於 \(0\)，當前就是最優解。
否則我們就再一次找到一個 \(a_{l,e}>0(1\le l\le m)\)，如果沒有，說明當前 \(e\) 無法變成非正係數，就說明解可以無窮大（unbounded）。
然後，我們旋轉 \(l\) 與 \(e\)。
重複這個操作就可以了。
退化與布蘭德規則
在進行上述過程中，我們可能會進入一個死迴圈，目標值不變，我們稱當前遇到了退化（degeneracy），退化可能會導致死迴圈。而應對它的方法就是布蘭德規則（bland），我們可以根據如下法則選擇 \(l\) 與 \(e\)：
- 在選擇 \(e\) 時，選擇下標最小的那個；
- 在選擇 \(l\) 時，選擇約束最緊的那個（也就是 \(\frac{a_{l,0}}{a_{l,e}}\) 最小的那個，它限制了取值範圍）。
根據布蘭德規則，我們可以寫出如下程式碼：
void simplex() {
while(1) {
ll l = 0, e = 0;
ld mn = inf;
for(ll i = 1; i <= n; i ++) {
if(a[0][i] > eps) {
e = i;
break;
}
}
if(!e) break;
for(ll i = 1; i <= m; i ++) {
if(a[i][e] > eps && a[i][0] / a[i][e] < mn) {
l = i;
mn = a[i][0] / a[i][e];
}
}
if(!l) {
printf("Unbounded");
exit(0);
}
pivot(l, e);
}
}
基本不可行解
我們發現，在求解基本可行解中，如果出現某一個基本變數小於 \(0\) 的情況，是不可行的。
比如：
那麼轉換為標準型：
轉換為鬆弛型：
我們發現，當 \(x_1=0,x_2=0\) 時，\(x_3=-10,x_4=-3\)，它是不滿足 \(x_3\ge0,x_4\ge0\) 的要求的，所以我們認為它不合法。
我們觀察以下基本不可行解有什麼特點，我們將上面的鬆弛型整理一下：
我們發現，當我們將這些非基本變數都取為 \(0\) 時，因為 \(a_{1,0}=-10\)，所以導致 \(x_3=-10\)。同理因為 \(a_{2,0}=-8\)，所以 \(x_4=-8\)。
也就是說，一個初始解可行當且僅當對於任意的 \(i\) 都滿足 \(a_{i,0}\ge0\)。
生成一個初始可行解有幾種方法：一種是建立一個輔助線性規劃（auxiliary linear program），但在演算法競賽中我們常用的是第二種方法，隨機選擇法：
我們先隨機找到一個 \(a_{l,0}<0\)，我們希望能將其變為大於等於 \(0\) 的。如果沒有這樣的 \(l\)，說明當前已經是初始可行解。
否則再隨機找到一個 \(a_{l,e}<0\)，根據上面的經驗，我們將 \(l\) 與 \(e\) 旋轉之後就可以把 \(a_{l,0}\) 變為大於等於 \(0\) 的。如果不存在這樣的 \(e\) 說明當前這個位置不能變成大於等於 \(0\) 的，那麼無解（infeasible）。
void init() {
while(1) {
ll l = 0, e = 0;
for(ll i = 1; i <= m; i ++) {
if(a[i][0] < -eps && (!l || rnd() % 2 == 0)) {
l = i;
}
}
if(!l) break;
for(ll i = 1; i <= n; i ++) {
if(a[l][i] < -eps && (!e || rnd() % 2 == 0)) {
e = i;
}
}
if(!e) {
printf("Infeasible");
exit(0);
}
pivot(l, e);
}
}
單純形演算法的幾何意義
我們換一個角度來看線性規劃：透過圖的視角。考慮以下問題：
我們可以畫出以下圖表：
容易發現，最優解必定在頂點上，不需要考慮內部點。（因為圍成的可行域一定是凸的）
同時，因為可行域是凸的，我們可以求每個頂點的高度，找出其中最高的一個，肯定就是最優點。
但是我們還有個更簡單的方法：
先找到一個頂點，然後從這個頂點，沿著某條邊線，走到下一個頂點，直到最優。方向的選擇可以有很多種，最多使用的是比較短視的方法：沿著最陡峭的那一條，追求當前步上升最快。
因為可行域它是凸的，就保證了要麼解無窮大要麼只有一個極值。
當我們進行一次旋轉操作時，相當於沿著一條邊移動到另一個頂點上，所以進行若干次操作後必定可以移動到最值上。
我們設移動的方向是 \(\lambda\)，距離是 \(\theta\)，我們進行旋轉操作相當於走 \(x'=x-\theta\lambda\)。
對應到單純形中的矩陣，例如還是機床廠問題：
我們的方向 \(\lambda\) 就相當與某一個非基本變數的那一列，比如此處我們選擇第一列，可以發現 \(4(-z)+-x_1+0x_2+2x_3+1x_4+0x_5=0\)，那麼就是 \(\lambda=(4,-1,0,2,1,0)\)。
我們走多少呢？走過多會超出區域，過少會達不到頂點。可以發現最多隻能是 \(\theta=5\)，也就是限制最緊的那一個：\(10\div2=5\)，我們就可以得到解：
\(x'=(0,0,0,10,8,7)-(4,-1,0,2,1,0)\times5=(-20,5,0,0,3,7)\)，發現是對應的：
同時，如果可行域不包括原點，那麼我們也是需要建立一個初始可行解的，否則我們就沒有在頂點上跑，矛盾。
單純形演算法的時間複雜度分析
我們發現，旋轉一次的時間複雜度為 \(O(nm)\)。假設旋轉 \(k\) 次，那麼時間複雜度就是 \(O(knm)\)。
在很長一段時間內，人們認為單純形是多項式時間複雜度的。直到 V. Klee and G. L. Minty[1972] 構造了一個例子，我們稱它為 Klee–Minty 問題：
其中，\(0<\delta\le\frac{1}{3}\)。
該問題的時間複雜度是質數級別的。該問題的可行域是頂點被擾動了的單位超立方體（unit hypercube），如果選擇全零作為初始可行解進行單純形演算法，在幾何意義下，單純形演算法將會遍歷每一個頂點，進行 \(2^n-1\) 次轉動操作，才可以得到最優解。
但是，Borgwardt (1982) 證明單純形演算法的平均複雜度是多項式時間的；Haimovich (1983) 證明了迭代次數的數學期望實際上是線性的；Spielman and Teng (2004) 引入了平滑型複雜度理論（smoothed analysis）。Spielman & Teng 定理斷言：線上性規劃問題上加入隨機高斯擾動，單純形演算法期望用多項式步數求解。
所以我們在最開始時進行若干次隨機擾動可以使得單純形演算法期望可以在多項式時間複雜度內求解。一般來說，旋轉的次數是在 \(2(n+m)\) 左右的。
線性規劃問題有更優的做法嗎？
答案是有的。儘管單純形是指數時間複雜度，但是 L. G. Khachiyan 提出的橢球法與 N. Karmarkar 提出的內點法具有多項式時間複雜度。這裡由於篇幅問題不展開。
而在一般的演算法競賽中，單純形更為常用且表現更好。
對偶定理
當我們的線性規劃是求最小值的同時約束也都是大於等於時，我們不僅可以透過將係數都乘以 \(-1\) 以轉換為一般型，也可以透過對偶定理實現。
我們稱原問題為 LP，對偶問題為 DP。
原問題有：
對偶問題有：
簡單來說，我們將單純形矩陣進行矩陣轉置後再進行樸素單純形即可。還是以機床廠問題為例：
它的 DP 是：
即原矩陣的轉置。
| 對偶問題（下）\原問題（右） | 最優解 | 無界解 | 無可行解 |
|---|---|---|---|
| 最優解 | \(x\) | \ | \ |
| 無界解 | \ | \ | \(x\) |
| 無可行解 | \ | \(x\) | 無法判斷 |
所以部分問題可以透過轉換為對偶問題以省略初始找可行解的過程。
全么模矩陣
若矩陣滿足任意一個子方陣的行列式為 \(0,-1,1\)，那麼我們稱這個矩陣為全么模矩陣（totally unimodular matrix）。
若矩陣是全么模矩陣，該線性規劃最優解為整數。
我們觀察後可以發現，如果一個規劃問題它的可行域多面體的所有頂點都是整數點的話（例如圖中 \(P\) 就滿足這個條件，而 \(P1\) 和 \(P2\) 都不滿足這個條件），那就可以滿足線性規劃最優解為整數了。
證明全么模矩陣可行域的頂點都在整點上：一個頂點無非是將一些線性無關的不等式改成等式後的線性方程組的解。如果矩陣 \(A\) 是全么模的，若滿足 \(A_S\) 是 \(A\) 的非奇異方陣，那麼就有 \(\det A_S=\pm1\)（因為全么模矩陣保證了 \(\det A_S=0\text{ or }\pm1\)，而非奇異方陣保證了 \(\det A_S\neq0\)）。若 \(\det A_S^{(i)}|b_S\) 是將方陣 \(A_S\) 的第 \(i\) 列替換為 \(b\) 得到的矩陣，那麼依克萊姆法則（Cramer's Rule）可知：
因為 \(A_S^{(i)}|b_S\) 方陣中的每一個數均為整數，所以 \(\det A_S^{(i)}|b_S\) 為整數。同時因為 \(\det A_S=\pm1\)，所以 \(x_i\) 必定為整數。以此類推，所以頂點均為整點。
以下命題中的矩陣是全么模矩陣：
- 無向二分圖的關聯矩陣是全么模矩陣（無向二分圖的關聯矩陣為行表示結點，列表示邊，如果結點和邊關聯，則單元格值為 \(1\)，否則為 \(0\)）；
- 有向圖的關聯矩陣是全么模矩陣（有向圖的關聯矩陣為行表示結點，列表示邊，每條邊與入點的單元格值為 \(1\)，與出點的單元格值為 \(-1\)，否則為 \(0\)）；
- 任何最大流、最小費用最大流的線性規劃都是全么模矩陣。
具體證明可以檢視參考資料。
證明一個矩陣是全么模其實還有一個騷操作：我們隨機造若干資料，如果這些資料的結果都是整數，那麼這個矩陣八成是全么模的。
例題
「UOJ #179」線性規劃
求解一個 \(n\) 個變數與 \(m\) 條約束的標準型線性規劃。
\(n,m\le20\)。
如果你實現樸素的單純形，你就會發現沒有透過 hack 資料，提交記錄。實際上我們可以在開頭進行若干次的隨機擾動以實現期望線性。
可惜的是依舊沒過（我太菜了），這次是被卡精度了，實現一下高精度浮點數或許能過。
#include <bits/stdc++.h>
using namespace std;
#define N 100
#define ll long long
#define ld long double
#define eps 1e-8
#define inf 1e15
ll n, m, t;
ld a[N][N];
ll r[N];
ld ans[N];
std::mt19937 rnd(time(0));
void pivot(ll l, ll e) {
swap(r[n + l], r[e]);
ld t = a[l][e];
a[l][e] = 1;
for(ll i = 0; i <= n; i ++) a[l][i] /= t;
for(ll i = 0; i <= m; i ++) if(i != l && abs(a[i][e]) > eps) {
t = a[i][e]; a[i][e] = 0;
for(ll j = 0; j <= n; j ++) {
a[i][j] -= a[l][j] * t;
}
}
}
void noise() {
for(ll i = 1; i <= 10; i ++) {
ll l = 0, e = 0;
for(ll i = 1; i <= m; i ++) {
if(abs(a[i][0]) > eps && (!l || rnd() % 2 == 0)) {
l = i;
}
}
if(!l) break;
for(ll i = 1; i <= n; i ++) {
if(abs(a[l][i]) > eps && (!e || rnd() % 2 == 0)) {
e = i;
}
}
if(!e) continue;
pivot(l, e);
}
}
void init() {
while(1) {
ll l = 0, e = 0;
for(ll i = 1; i <= m; i ++) {
if(a[i][0] < -eps && (!l || rnd() % 2 == 0)) {
l = i;
}
}
if(!l) break;
for(ll i = 1; i <= n; i ++) {
if(a[l][i] < -eps && (!e || rnd() % 2 == 0)) {
e = i;
}
}
if(!e) {
printf("Infeasible");
exit(0);
}
pivot(l, e);
}
}
void simplex() {
while(1) {
ll l = 0, e = 0;
ld mn = inf;
for(ll i = 1; i <= n; i ++) {
if(a[0][i] > eps) {
e = i;
break;
}
}
if(!e) break;
for(ll i = 1; i <= m; i ++) {
if(a[i][e] > eps && a[i][0] / a[i][e] < mn) {
l = i;
mn = a[i][0] / a[i][e];
}
}
if(!l) {
printf("Unbounded");
exit(0);
}
pivot(l, e);
}
}
int main() {
scanf("%lld %lld %lld", &n, &m, &t);
for(ll i = 1; i <= n; i ++) {
scanf("%Lf", &a[0][i]);
}
for(ll i = 1; i <= m; i ++) {
for(ll j = 1; j <= n; j ++) {
scanf("%Lf", &a[i][j]);
}
scanf("%Lf", &a[i][0]);
}
for(ll i = 1; i <= n; i ++) r[i] = i;
noise();
init();
simplex();
if(abs(a[0][0]) < eps) printf("0\n");
else printf("%.10Lf\n", -a[0][0]);
if(t) {
for(ll i = 1; i <= m; i ++) ans[r[n + i]] = a[i][0];
for(ll i = 1; i <= n; i ++) printf("%.10Lf ", ans[i]);
}
}
「NOI2008」志願者招募
\(n\) 天，每天需要 \(a_i\) 個人。有 \(m\) 種人，每種人可以從 \(s_i\) 工作到 \(t_i\)，費用為 \(c_i\)，求最小費用。
\(1\le n\le1000,1\le m\le10000\)。
單純形板題，設變數 \(x_1\sim x_m\) 表示每種志願者招募多少個，第 \(i\) 條約束是第 \(i\) 天可以工作志願者的和大於等於 \(a_i\)。最小化 \(x_ic_i\) 的和。
發現既是求最小值又是大於等於的約束，所以可以透過對偶實現。
我們發現志願者都要求是正數的，所以需要證明這個矩陣是全么模矩陣。
但是我們發現矩陣比較特殊，都是 \(0\) 和 \(1\)，且每一列的 \(1\) 有且僅有連續的一段。
我們可以將這個矩陣的每一列乘上 \(-1\) 再累加到後一列，我們發現這些操作對矩陣的行列式是沒有影響的：
此時該矩陣屬於有向圖的關聯矩陣，是全么模矩陣，所以可以得到最優整數解。
同理，我們可以在開頭進行若干次的隨機擾動以實現期望線性。
#include <bits/stdc++.h>
using namespace std;
#define ld double
#define N 1010
#define M 10010
const ld eps = 1e-8, inf = 1e9;
int n, m;
ld a[M][N];
std::mt19937 rnd(time(0));
void pivot(int l, int e) {
ld t = a[l][e];
a[l][e] = 1;
for(int i = 0; i <= n; i ++) {
a[l][i] /= t;
}
for(int i = 0; i <= m; i ++) if(i != l && abs(a[i][e]) > eps){
t = a[i][e];
a[i][e] = 0;
for(int j = 0; j <= n; j ++) {
a[i][j] -= a[l][j] * t;
}
}
}
void noise() {
for(int i = 1; i <= 100; i ++) {
int l = 0, e = 0;
for(int i = 1; i <= m; i ++) {
if((a[i][0] > eps || a[i][0] < -eps) && (!l || rnd() % 2 == 0)) {
l = i;
}
}
if(!l) break;
for(int i = 1; i <= n; i ++) {
if((a[l][i] > eps || a[l][i] < -eps) && (!e || rnd() % 2 == 0)) {
e = i;
}
}
if(!e) continue;
pivot(l, e);
}
}
void init() {
while(1) {
int l = 0, e = 0;
for(int i = 1; i <= m; i ++) {
if(a[i][0] < -eps && (!l || rnd() % 2)) l = i;
}
if(!l) break;
for(int i = 1; i <= n; i ++) {
if(a[l][i] < -eps && (!e || rnd() % 2)) e = i;
}
pivot(l, e);
}
}
void simplex() {
while(1) {
int l = 0, e = 0;
ld mn = inf;
for(int i = 1; i <= n; i ++) {
if(a[0][i] > eps) {
e = i;
break;
}
}
if(!e) break;
for(int i = 1; i <= m; i ++) {
if(a[i][e] > eps && a[i][0] / a[i][e] < mn) {
l = i;
mn = a[i][0] / a[i][e];
}
}
pivot(l, e);
}
}
int main() {
scanf("%d %d", &n, &m);
for(int i = 1; i <= n; i ++) {
scanf("%lf", &a[0][i]);
}
for(int i = 1; i <= m; i ++) {
int s, t;
scanf("%d %d %lf", &s, &t, &a[i][0]);
for(int j = s; j <= t; j ++) {
a[i][j] = 1.0;
}
}
noise();
init();
simplex();
printf("%.lf", -a[0][0]);
}
「ABC231H」Minimum Coloring
一個 \(H \times W\) 的網格圖，初始所有點都是白色的。
有 \(N\) 個點可以被改變成黑色，這 \(N\) 個點的座標是 \(a_i,b_i\)，改變顏色的代價是 \(c_i\)。
你需要找到最小代價使得每行每列都至少有一個黑色節點。
資料保證有解。
\(1\leq N,H,W\leq10^3\)
首先，我們設 \(x_i=0/1\) 表示第 \(i\) 個塗不塗。然後為了保證每一行都有，前 \(h\) 個約束為第 \(i\) 行所有的 \(x\) 加起來大於等於 \(1\)；為了保證每一列都有，後 \(w\) 個約束為第 \(i-h\) 列所有的 \(x\) 加起來大於等於 \(1\)。目標就是要求 \(\sum x_ic_i\) 儘可能小。
對偶一下就可以了。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define N 2010
double eps = 1e-8, inf = 1e15;
ll h, w, n, m;
double a[N][N];
void pivot(ll l, ll e) {
double t = a[l][e];
a[l][e] = 1;
for(ll i = 0; i <= m; i ++) {
a[l][e] /= t;
}
for(ll i = 0; i <= n; i ++) if(i != l && abs(a[l][e]) > eps) {
t = a[i][e], a[i][e] = 0;
for(ll j = 0; j <= m; j ++) {
a[i][j] -= t * a[l][j];
}
}
}
void simplex() {
while(1) {
ll l = 0, e = 0;
double mn = inf;
for(ll i = 1; i <= m; i ++) if(a[0][i] > eps) {
e = i;
break;
}
if(!e) break;
for(ll i = 1; i <= n; i ++) if(a[i][e] > eps && a[i][0] / a[i][e] < mn) {
l = i;
mn = a[i][0] / a[i][e];
}
pivot(l, e);
}
}
int main() {
scanf("%lld %lld %lld", &h, &w, &n);
m = h + w;
for(ll i = 1; i <= n; i ++) {
ll x, y, c;
scanf("%lld %lld %lld", &x, &y, &c);
a[i][0] = c;
a[i][x] = 1;
a[i][h + y] = 1;
}
for(ll i = 1; i <= m; i ++) {
a[0][i] = 1;
}
simplex();
printf("%.lf", -a[0][0]);
}
「SHOI2004」最小生成樹
給定一個 \(n\) 點 \(m\) 邊的簡單圖，每條邊將邊權修改為 \(w_i'\gets w_i+c\) 具有代價 \(|c|\)，給定簡單圖上的一棵生成樹 \(T\)，要求最小的代價修改簡單圖上的每一個邊權使得這顆生成樹 \(T\) 變為最小生成樹。
\(1\le n\le 50,1\le m\le 1500\)
對於一條非樹邊 \(j\)，它肯定跟若干條樹邊構成了一個環，那麼這個非樹邊權值一定要大於等於環上的所有邊。設其中一條是 \(i\)，我們有一個貪心的策略——減小 \(i\) 的邊權，增加 \(j\) 的邊權。
所以有：\(w_i-x_i\le w_j+x_j\)。
也就是 \(w_i-w_j\le x_i+x_j\)。目標是最小化 \(x\) 的和。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define N 60
#define M 1510
const double eps = 1e-8, inf = 1e9;
mt19937 rnd(time(0));
ll n, m;
bool vis[N];
ll head[N], nxt[M * 2], to[M * 2], cnt;
ll road[N][N];
ll U[M], V[M], W[M];
void addEdge(ll u, ll v) {
cnt ++;
to[cnt] = v;
nxt[cnt] = head[u];
head[u] = cnt;
}
ll fa[N], dep[N];
void dfs(ll u) {
for(ll i = head[u]; i; i = nxt[i]) {
ll v = to[i];
if(v == fa[u]) continue;
fa[v] = u;
dep[v] = dep[u] + 1;
dfs(v);
}
}
double a[M][M];
ll tot;
void pivot(ll l, ll e) {
double t = a[l][e]; a[l][e] = 1;
for(ll i = 0; i <= tot; i ++) {
a[l][i] /= t;
}
for(ll i = 0; i <= m; i ++) if(i != l && abs(a[i][e]) > eps) {
t = a[i][e]; a[i][e] = 0;
for(ll j = 0; j <= tot; j ++) {
a[i][j] -= t * a[l][j];
}
}
}
void init() {
while(1) {
ll l = 0, e = 0;
for(ll i = 1; i <= m; i ++) {
if(a[i][0] < -eps && (!l || rnd() % 2 == 0)) {
l = i;
}
}
if(!l) break;
for(ll i = 1; i <= tot; i ++) {
if(a[l][i] < -eps && (!e || rnd() % 2 == 0)) {
e = i;
}
}
pivot(l, e);
}
}
void simplex() {
while(1) {
ll l = 0, e = 0;
double mn = inf;
for(ll i = 1; i <= tot; i ++) {
if(a[0][i] > eps) {
e = i;
break;
}
}
if(!e) break;
for(ll i = 1; i <= m; i ++) {
if(a[i][e] > eps && a[i][0] / a[i][e] < mn) {
l = i;
mn = a[i][0] / a[i][e];
}
}
pivot(l, e);
}
}
int main() {
scanf("%lld %lld", &n, &m);
for(ll i = 1; i <= m; i ++) {
scanf("%lld %lld %lld", &U[i], &V[i], &W[i]);
road[U[i]][V[i]] = road[V[i]][U[i]] = i;
}
for(ll i = 1; i < n; i ++) {
ll u, v;
scanf("%lld %lld", &u, &v);
addEdge(u, v);
addEdge(v, u);
}
dfs(1);
for(ll i = 1; i <= m; i ++) {
ll u = U[i], v = V[i];
a[i][0] = 1;
while(u != v) {
if(dep[u] < dep[v]) swap(u, v);
ll x = u;
u = fa[u];
if(W[road[x][u]] > W[i]) {
// W[road[x][u]] - x[road[x][u]] <= W[i] + x[i]
// x[i] + x[road[x][u]] >= W[road[x][u]] - W[i]
a[0][++ tot] = W[road[x][u]] - W[i];
a[i][tot] = 1;
a[road[x][u]][tot] = 1;
}
}
}
init();
simplex();
if(a[0][0] > -eps && a[0][0] < eps) printf("0");
else printf("%.lf", -a[0][0]);
}
「CF1430G」Yet Another DAG Problem
給定一個 \(n\) 點 \(m\) 邊的有向無環圖，每條邊都有 \(w_i\) 的權重，給每個點分配權值 \(a_i\)，對於每條連線 \((u,v)\) 的邊，定義其權值為 \(b_i=a_u-a_v\)，要求：
\(b_i>0\)
\(\sum w_ib_i\) 最小
請輸出一種分配方案。
\(1\le n\le 18,1\le m\le n(n-1)/2\)
這裡 \(b_i>0\) 看似不是線性規劃，實際上因為 \(b\) 為整數，所以同等於 \(b_i\ge 1\)。每條邊可以被描述為限制 \(1a_u+-1a_v\ge1\)。\(\min \sum w_ib_i\) 可以被拆分為 \(\min\sum(w_i a_u+-w_i a_v)\)，然後單純形做就好了。
因為要求每一個點的取值，所以我不用對偶，直接將每一項係數乘以 \(-1\) 也可以達到同樣的效果。
有向圖的關聯矩陣是全么模矩陣，所以最優解是整數的。
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define db long double
#define N 20
#define M 400
const db inf = 1e9, eps = 1e-8;
ll n, m;
db a[M][N], ans[N];
ll r[N + M];
mt19937 rnd(114514191);
void pivot(ll l, ll e) {
swap(r[n + l], r[e]);
db t = a[l][e];
a[l][e] = 1;
for(ll i = 0; i <= n; i ++) {
a[l][i] /= t;
}
for(ll i = 0; i <= m; i ++) if(i != l && abs(a[i][e]) > eps) {
t = a[i][e]; a[i][e] = 0;
for(ll j = 0; j <= n; j ++) {
a[i][j] -= t * a[l][j];
}
}
}
void init() {
while(1) {
ll l = 0, e = 0;
for(ll i = 1; i <= m; i ++) {
if(a[i][0] < -eps && (!l || rnd() % 2)) {
l = i;
}
}
if(!l) break;
for(ll i = 1; i <= n; i ++) {
if(a[l][i] < -eps && (!e || rnd() % 2)) {
e = i;
}
}
pivot(l, e);
}
}
void simplex() {
while(1) {
ll l = 0, e = 0;
db mn = inf;
for(ll i = 1; i <= n; i ++) {
if(a[0][i] > eps) {
e = i;
break;
}
}
if(!e) break;
for(ll i = 1; i <= m; i ++) {
if(a[i][e] > eps && a[i][0] / a[i][e] < mn) {
l = i;
mn = a[i][0] / a[i][e];
}
}
pivot(l, e);
}
}
int main() {
scanf("%lld %lld", &n, &m);
for(ll i = 1; i <= m; i ++) {
ll u, v, w;
scanf("%lld %lld %lld", &u, &v, &w);
a[i][u] = -1;
a[i][v] = 1;
a[i][0] = -1;
a[0][u] -= w;
a[0][v] += w;
}
for(ll i = 1; i <= n; i ++) {
r[i] = i;
}
init();
simplex();
for(ll i = 1; i <= m; i ++) {
ans[r[n + i]] = a[i][0];
}
for(ll i = 1; i <= n; i ++) {
printf("%.Lf ", ans[i]);
}
}
參考資料
- 線性規劃-單純形演算法詳解 | 細語呢喃
- 數學規劃（3）單純形法的進一步討論
- 線性規劃 | Daltao's blog!
轉載時請附上鍊接：https://www.cnblogs.com/znpdco/p/18147653