線性規劃之單純形講解
1.作用
單純形法是解決線性規劃問題的一個有效的 ofollow,noindex">演算法 。線性規劃就是在一組線性約束條件下，求解目標函式最優解的問題。
2. 線性規劃的一般形式
在約束條件下，尋找目標函式z的最大值。
3.線性規劃的可行域
滿足線性規劃問題約束條件的所有點組成的集合就是線性規劃的可行域。若可行域有界（以下主要考慮有界可行域），線性規劃問題的目標函式最優解必然在可行域的頂點上達到最優。
單純形法就是通過設定不同的基向量，經過矩陣的線性變換，求得基可行解（可行域頂點），並判斷該解是否最優，否則繼續設定另一組基向量，重複執行以上步驟，直到找到最優解。所以，單純形法的求解過程是一個迴圈迭代的過程。
圖1 可行域
4.線性規劃的標準形式
在說明單純形法的原理之前，需要明白線性規劃的標準形式。因為單純形演算法是通過線性規劃的標準形來求解的。一般，規定線性規劃的標準形式為：
寫成矩陣形式：
標準形的形式為：
1）目標函式要求max
2）約束條件均為等式
3）決策變數為非負約束
普通線性規劃化為標準形：
1）若目標函式為最小化，可以通過取負，求最大化
2）約束不等式為小於等於不等式，可以在左端加入非負鬆弛變數，轉變為等式，比如：
同理，約束不等式為大於等於不等式時，可以在左端減去一個非負鬆弛變數，變為等式。
3）若存在取值無約束的變數，可轉變為兩個非負變數的差，比如：
本文最開始的線性規劃問題轉化為標準形為：
5.單純形法
5.1幾何意義
在標準形中，有m個約束條件（不包括非負約束），n個決策變數，且（n>=m）。首先，選取m個基變數 
Line"/> ，基變數對應約束係數矩陣的列向量線性無關。通過矩陣的線性變換，基變數可由非基變量表示：
如果令非基變數等於0，可求得基變數的值 ：
如果為可行解的話，Ci大於0。那麼它的幾何意義是什麼呢？
還是通過上述具體的線性規劃問題來說明。
如果選擇x2、x3為基變數，那麼令x1、x4等於0，可以去求解基變數x2、x3的值。對係數矩陣做行變換，如下所示，x2=9/2，x3=15/2
X1=0表示可行解在x軸上；X4=0表示可行解在x1+2x2=9的直線上。那麼，求得的可行解即表示這兩條直線的交點，也是可行域的頂點，
所以，通過選擇不同的基變數，可以獲得不同的可行域的頂點。
5.2如何判斷最優
如前所述，基變數可由非基變量表示：
目標函式z也可以完全由非基變量表示：
當達到最優解時，所有的 
應小於等於0。當存在 j ， >0時，當前解不是最優解，為什麼？
當前的目標函式值為z0，其中所有的非基變數值均取0。由之前分析可知， 
=0代表可行域的某個邊界，是 的最小值。如果可行解逐步離開這個邊界， 會變大，因為 >0，顯然目標函式的取值也會變大，所以當前解不是最優解。我們需要尋找新的基變數。
5.3如何選擇新的基變數
如果存在多個 >0，選擇最大的 >0對應的變數作為基變數，這表示目標函式隨著 的增加，增長的最快。
5.4如何選擇被替換的基變數
假如我們選擇非基變數 
作為下一輪的基變數，那麼被替換基變數 在下一輪中作為非基變數，等於0。選擇 的原則：替換後應該儘量使 值最大（因為上面已分析過，目標函式會隨著 的增大而增大）。
繼續通過上面的例子來說明：
從最後一行可以看到，x1的係數為1/2>0，所以選x2、x3為基變數並沒有是目標函式達到最優。下一輪選取x1作為基變數，替換x2、x3中的某個變數。
第一行是符號
第二行：若x1替換x3作為基變數，x3=0時，x1=(15/2)/(3/2)=5
第三行：若x1替換x2作為基變數，x2=0時，x1=(9/2)/(1/2)=9
顯然，應該把x2作為非基變數。
5.5終止條件
當目標函式用非基變數的線性組合表示時，所有的係數均不大於0，則表示目標函式達到最優。
如果，有一個非基變數的係數為0，其他的均小於0，表示目標函式的最優解有無窮多個。這是因為目標函式的梯度與某一邊界正交，在這個邊界上，目標函式的取值均相等，且為最優。
使用單純型法來求解線性規劃，輸入單純型法的鬆弛形式，是一個大矩陣，第一行為目標函式的係數，且最後一個數字為當前軸值下的 z 值。下面每一行代表一個約束，數字代表係數每行最後一個數字代表 b 值。
演算法和使用單純性表求解線性規劃相同。
對於線性規劃問題：
Max x1 + 14* x2 + 6*x3
s . t . x1 + x2 + x3 <= 4
x1 <= 2
x3 <= 3
3*x2 + x3 <= 6
x1,x2,x3 >= 0
我們可以得到其鬆弛形式：
Maxx1 + 14*x2 + 6*x3
s.t. x1 + x2 + x3 + x4 = 4
x1 + x5 = 2
x3 + x6 = 3
3*x2 + x3 + x7 = 6
x1 , x2 , x3 , x4 , x5 , x6 , x7 ≥ 0
我們可以構造單純性表，其中最後一行打星的列為軸值。
| x1 | x2 | x3 | x4 | x5 | x6 | x7 | b |
| c1=1 | c2=14 | c3=6 | c4=0 | c5=0 | c6=0 | c7=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，再上面的例子中我們選擇c2作為新的軸，加入軸集合中，那麼誰該出軸呢？
其實我們由於每個x都大於零，對於x2它的增加是有所限制的，如果x2過大，由於其他的限制條件，就會使得其他的x小於零，於是我們應該讓x2一直增大，直到有一個其他的x剛好等於0為止，那麼這個x就被換出軸。
我們可以發現，對於約束方程1，即第一行約束，x2最大可以為4（4/1），對於約束方程4，x2最大可以為3(6/3），因此x2最大隻能為他們之間最小的那個，這樣才能保證每個x都大於零。因此使用第4行，來對各行進行高斯行變換，使得二列第四行中的每個x都變成零，也包括c2。這樣我們就完成了把x2入軸，x7出軸的過程。變換後的單純性表為：
| x1 | x2 | x3 | x4 | x5 | x6 | x7 | b |
| c1=1 | c2=0 | c3=1.33 | c4=0 | c5=0 | c6=0 | c7=-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 |
| * | * | * | * |
繼續計算，我們得到：
| x1 | x2 | x3 | x4 | x5 | x6 | x7 | b |
| c1=-1 | c2=0 | c3=0 | c4=0 | c5=-2 | c6=0 | c7=0 | -z=-32 |
| 1.5 | 0 | 1 | 1.5 | 0 | 0 | -0.5 | 3 |
| 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的係數全部小於零，即增大任何非軸的x值並不能使得目標函式最大，從而得到最優解32.
整個過程程式碼如下所示：
1 #include <bits/stdc++.h> 2 using namespace std; 3 vector<vector<double> > Matrix; 4 double Z; 5 set<int> P; 6 size_t cn, bn; 7 8 bool Pivot(pair<size_t, size_t> &p)//返回0表示所有的非軸元素都小於0 9 { 10int x = 0, y = 0; 11double cmax = -INT_MAX; 12vector<double> C = Matrix[0]; 13vector<double> B; 14 15for( size_t i = 0 ; i < bn ; i++ ) 16{ 17B.push_back(Matrix[i][cn-1]); 18} 19 20for( size_t i = 0 ; i < C.size(); i++ )//在非軸元素中找最大的c 21{ 22if( cmax < C[i] && P.find(i) == P.end()) 23{ 24cmax = C[i]; 25y = i; 26} 27} 28if( cmax < 0 ) 29{ 30return 0; 31} 32 33double bmin = INT_MAX; 34for( size_t i = 1 ; i < bn ; i++ ) 35{ 36double tmp = B[i]/Matrix[i][y]; 37if( Matrix[i][y] != 0 && bmin > tmp ) 38{ 39bmin = tmp; 40x = i; 41} 42} 43 44p = make_pair(x, y); 45 46for( set<int>::iterator it = P.begin() ; it != P.end() ; it++) 47{ 48if( Matrix[x][*it] != 0 ) 49{ 50//cout<<"erase "<<*it<<endl; 51P.erase(*it); 52break; 53} 54} 55P.insert(y); 56//cout<<"add "<<y<<endl; 57return true; 58 } 59 60 void pnt() 61 { 62for( size_t i = 0 ; i < Matrix.size() ; i++ ) 63{ 64for( size_t j = 0 ; j < Matrix[0].size() ; j++ ) 65{ 66cout<<Matrix[i][j]<<"\t"; 67} 68cout<<endl; 69} 70cout<<"result z:"<<-Matrix[0][cn-1]<<endl; 71 } 72 73 void Gaussian(pair<size_t, size_t> p)//行變換 74 { 75size_tx = p.first; 76size_t y = p.second; 77double norm = Matrix[x][y]; 78for( size_t i = 0 ; i < cn ; i++ )//主行歸一化 79{ 80Matrix[x][i] /= norm; 81} 82for( size_t i = 0 ; i < bn && i != x; i++ ) 83{ 84if( Matrix[i][y] != 0) 85{ 86double tmpnorm = Matrix[i][y]; 87for( size_t j = 0 ; j < cn ; j++ ) 88{ 89Matrix[i][j] = Matrix[i][j] - tmpnorm * Matrix[x][j]; 90} 91} 92} 93 } 94 95 void solve() 96 { 97pair<size_t, size_t> t; 98while(1) 99{ 100 101pnt(); 102if( Pivot(t) == 0 ) 103{ 104return; 105} 106cout<<t.first<<" "<<t.second<<endl; 107for( set<int>::iterator it = P.begin(); it != P.end(); it++ ) 108{ 109cout<<*it<<" "; 110} 111cout<<endl; 112Gaussian(t); 113} 114 } 115 116 int main(int argc, char *argv[]) 117 { 118//ifstream fin; 119//fin.open("./test"); 120cin>>cn>>bn; 121for( size_t i = 0 ; i < bn ; i++ ) 122{ 123vector<double> vectmp; 124for( size_t j = 0 ; j < cn ; j++) 125{ 126double tmp = 0; 127cin>>tmp; 128vectmp.push_back(tmp); 129} 130Matrix.push_back(vectmp); 131} 132 133for( size_t i = 0 ; i < bn-1 ; i++ ) 134{ 135P.insert(cn-i-2); 136} 137solve(); 138 } 139 ///////////////////////////////////// 140 //glpk input: 141 ///* Variables */ 142 //var x1 >= 0; 143 //var x2 >= 0; 144 //var x3 >= 0; 145 ///* Object function */ 146 //maximize z: x1 + 14*x2 + 6*x3; 147 ///* Constrains */ 148 //s.t. con1: x1 + x2 + x3 <= 4; 149 //s.t. con2: x1<= 2; 150 //s.t. con3: x3<= 3; 151 //s.t. con4: 3*x2 + x3<= 6; 152 //end; 153 ///////////////////////////////////// 154 //myinput: 155 /* 156 8 5 157 1 14 6 0 0 0 0 0 158 1 1 1 1 0 0 0 4 159 1 0 0 0 1 0 0 2 160 0 0 1 0 0 1 0 3 161 0 3 1 0 0 0 1 6 162 */ 163 /////////////////////////////////////
【理論羅列】：
1.標準型
m個約束 n個變數用x向量表示 A是一個m*n的矩陣 c是一個n的向量 b是一個m的向量
最大化 cx
滿足約束 Ax<=b x>0
2.鬆弛型
基本變數 B |B|=m 一個約束對應一個 表示鬆弛量 叫做鬆弛變數(基本變數)
非基變數 N |N|=n
xn+i=bi-sigma{aijxj}>=0
3.替入變數 xe(非基變數)
替出變數 xl(基本變數)
4.可行解
基本解：所有非基變數設為0
基本可行解
5.單純形法的過程中B和N不斷交換，在n維空間中不斷走
“相當於不等式上的高斯消元”
【程式碼實現】：
pivot是轉動操作
基本思想就是改寫l這個約束為xe作為基本變數，然後把這個新xe的值帶到其他約束和目標函式中，就消去xe了
改寫和帶入時要修改b和a 目標函式則是 c和v
轉動時l和e並沒有像演算法導論上一樣a矩陣用了兩行分別是a[l][]和a[e][]（這樣佔用記憶體大），而是用了同一行，這樣a矩陣的行數=|B|，列數=|N|
也就是說，約束條件只用m個，儘管B和N不斷交換，但同一時間還是隻有m個約束(基本變數)n個非基變數
注意改寫成鬆弛型後a矩陣實際係數為負
（一個優化 a[i][e]為0的約束沒必要帶入了
simplex是主過程
基本思想是找到一個c[e]>0的，然後找對這個e限制最緊的l，轉動這組l e
注意精度控制eps
c[e]>eps
還有找l的時候a[i][e]>eps才行
【對偶原理】：
1.原始線性規劃 對偶線性規劃
2.對於
最大化 cx
滿足約束 Ax<=b x>0
對偶問題為
最小化 bx
滿足約束 ATx>=c x>0 （AT為A的轉置）
可以轉化很多問題來避免初始解不可行
說從前有個線性規劃
min c x^T
Ax = b
x >= 0
這裡面A是矩陣，x、b、c都是向量
x^T表示轉置
啊……我們假設以上出現的所有元素都是整數……
那麼Ax = b要是恰好方程數等於未知數數，且解出來恰好為正數是不是就超開心？ （假設是線性無關的）
根據那個啥法則，x_i = det(A_i) / det(A)
det(A)表示A的行列式
A_i表示把A的第i列換為b之後的矩陣
如果det(A_i)恰好是det(A)的倍數那不就超開心？這樣
但是現實是殘酷的，往往這傢伙會除出來小數，然後整數規劃就坑爹了。
我們現在還是假定“恰好方程數等於未知數數，且解出來恰好為正數”
我們再加個限制：det(A) = 1或-1
就233了吧。
解一定是整數了。
於是可以順利變為整數規劃。我們把det(A) = 1, -1的矩陣稱為么模矩陣。
但是現實是殘酷的，“恰好方程數等於未知數數，且解出來恰好為正數”往往不滿足。
但是其實沒關係。由於每個方程都對應著一個平面，所以解的空間是單純形，最優解一定會出現在頂點上。
何為頂點？就是平面的交點。
何為平面？一共m + n個：Ax = b是m個方程，x = 0是n個方程。（本來是x >= 0，我們只靠慮切割空間的平面……）
要是頂點都是整點不是超開心？等價於從這m + n個方程中任取n個方程把它解掉得到的解是整數解。
通過前面的分析我們知道，如果恰好選取的這n個方程的係數矩陣的行列式值為1，-1就超開心了。當然要是行列式值為0對應著無解或無窮多解的情況，它又不是頂點管它做甚……
考察係數矩陣
一個是A，好大好大
另一個是x = 0的係數，易知就是單位矩陣I
你從I中選哪幾行……由於行列式的性質……一行*k加到另一行上去行列式的值不變，那麼對應的未知數就會被幹掉……
所以等價於……
從A中任意選取是一個子方陣，它的行列式的值只可能為1, -1, 0。
這樣的矩陣我們稱為全么模矩陣。
OrzOrzOrz
附上一道題BZOJ 1061
1 /* 2 線性規劃單純形 3 1.原始線性規劃 對偶線性規劃對於最大化 cx滿足約束 Ax<=b,x>0; 4 2.對偶問題為最小化 bx滿足約束 ATx>=c,x>0（AT為A的轉置）; 5 可以轉化很多問題來避免初始解不可行 6 */ 7 #include<bits/stdc++.h> 8 using namespace std; 9 #define clr(a,b) memset(a,b,sizeof a) 10 #define lowbit(x) x&-x 11 typedef long long ll; 12 const int INF=0x3f3f3f3f; 13 inline int read() 14 { 15int x=0,f=1;char ch=getchar(); 16while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();} 17while(ch>='0'&&ch<='9') {x=x*10+ch-'0';ch=getchar();} 18return x*f; 19 } 20 const int M=10005; 21 const int N=1005; 22 const double eps=1e-6; 23 int n,m; 24 double a[M][N],b[M],c[N],v; 25 inline void pivot(int l,int e)//矩陣的轉秩 26 { 27b[l]/=a[l][e]; 28for(int j=1;j<=n;++j) 29{ 30if(j!=e) a[l][j]/=a[l][e]; 31} 32a[l][e]=1/a[l][e]; 33for(int i=1;i<=m;++i) 34{ 35if(i!=l&&fabs(a[i][e])>0) 36{ 37b[i]-=a[i][e]*b[l]; 38for(int j=1;j<=n;++j) 39{ 40if(j!=e) a[i][j]-=a[i][e]*a[l][j]; 41} 42a[i][e]=-a[i][e]*a[l][e]; 43} 44} 45v+=c[e]*b[l]; 46for(int j=1;j<=n;++j) 47{ 48if(j!=e) c[j]-=c[e]*a[l][j]; 49} 50c[e]=-c[e]*a[l][e]; 51 } 52 53 inline double simplex() 54 { 55while(1) 56{ 57int e=0,l=0; 58for(e=1;e<=n;++e) 59{ 60if(c[e]>eps) break; 61} 62if(e==n+1) return v; 63double mn=INF; 64for(int i=1;i<=m;++i) 65{ 66if(a[i][e]>eps&&mn>b[i]/a[i][e]) mn=b[i]/a[i][e],l=i; 67} 68if(mn==INF) return INF; 69pivot(l,e); 70} 71 } 72 73 int main() 74 { 75n=read(),m=read(); 76for(int i=1;i<=n;++i) c[i]=read(); 77for(int i=1;i<=m;++i) 78{ 79int s,t; 80s=read(),t=read(); 81for(int j=s;j<=t;++j) a[i][j]=1; 82b[i]=read(); 83} 84printf("%d\n",(int)(simplex()+0.5)); 85return 0; 86 } View Code
附上幾道題的題號，練習學習一下：
BZOJ 3550
BZOJ 3112
BZOJ 3265
BZOJ 1061
POJ 1275
標準型:
轉化為標準型:
鬆弛型:
單純型演算法
轉軸：
初始化：
程式碼實現：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #include<cmath> 6 using namespace std; 7 typedef long long ll; 8 const int N=25; 9 const double eps=1e-8,INF=1e15; 10 inline int read(){ 11char c=getchar();int x=0,f=1; 12while(c<'0'||c>'9'){if(c=='-')f=-1; c=getchar();} 13while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();} 14return x*f; 15 } 16 int n,m,type; 17 double a[N][N],ans[N]; 18 int id[N<<1]; 19 int q[N]; 20 void Pivot(int l,int e){ 21swap(id[n+l],id[e]); 22double t=a[l][e]; a[l][e]=1; 23for(int j=0;j<=n;j++) a[l][j]/=t; 24for(int i=0;i<=m;i++) if(i!=l && abs(a[i][e])>eps){ 25t=a[i][e]; a[i][e]=0; 26for(int j=0;j<=n;j++) a[i][j]-=a[l][j]*t; 27} 28 } 29 bool init(){ 30while(true){ 31int e=0,l=0; 32for(int i=1;i<=m;i++) if(a[i][0]<-eps && (!l||(rand()&1))) l=i; 33if(!l) break; 34for(int j=1;j<=n;j++) if(a[l][j]<-eps && (!e||(rand()&1))) e=j; 35if(!e) {puts("Infeasible");return false;} 36Pivot(l,e); 37} 38return true; 39 } 40 bool simplex(){ 41while(true){ 42int l=0,e=0; double mn=INF; 43for(int j=1;j<=n;j++) 44if(a[0][j]>eps) {e=j;break;} 45if(!e) break; 46for(int i=1;i<=m;i++) if(a[i][e]>eps && a[i][0]/a[i][e]<mn) 47mn=a[i][0]/a[i][e],l=i; 48if(!l) {puts("Unbounded");return false;} 49Pivot(l,e); 50} 51return true; 52 } 53 int main(){ 54freopen("in","r",stdin); 55srand(317); 56n=read();m=read();type=read(); 57for(int i=1;i<=n;i++) a[0][i]=read(); 58for(int i=1;i<=m;i++){ 59for(int j=1;j<=n;j++) a[i][j]=read(); 60a[i][0]=read(); 61} 62for(int i=1;i<=n;i++) id[i]=i; 63if(init() && simplex()){ 64printf("%.8lf\n",-a[0][0]); 65if(type){ 66for(int i=1;i<=m;i++) ans[id[n+i]]=a[i][0]; 67for(int i=1;i<=n;i++) printf("%.8lf ",ans[i]); 68} 69} 70 }
對偶性：
全么模矩陣(totally unimodular matrix)
充分條件：
