单纯形法只有两个约束条件_线性规划之单纯形法【超详解+图解】
使用單純型法來求解線性規劃，輸入單純型法的松弛形式，是一個大矩陣，第一行為目標函數的系數，且最后一個數字為當前軸值下的 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
我們可以得到其松弛形式：
Max x1 + 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
2 using namespacestd;3 vector >Matrix;4 doubleZ;5 setP;6 size_t cn, bn;7
8 bool Pivot(pair &p)//返回0表示所有的非軸元素都小于0
9 {10 int x = 0, y = 0;11 double cmax = -INT_MAX;12 vector C = Matrix[0];13 vectorB;14
15 for( size_t i = 0 ; i < bn ; i++)16 {17 B.push_back(Matrix[i][cn-1]);18 }19
20 for( size_t i = 0 ; i < C.size(); i++ )//在非軸元素中找最大的c
21 {22 if( cmax < C[i] && P.find(i) ==P.end())23 {24 cmax =C[i];25 y =i;26 }27 }28 if( cmax < 0)29 {30 return 0;31 }32
33 double bmin =INT_MAX;34 for( size_t i = 1 ; i < bn ; i++)35 {36 double tmp = B[i]/Matrix[i][y];37 if( Matrix[i][y] != 0 && bmin >tmp )38 {39 bmin =tmp;40 x =i;41 }42 }43
44 p =make_pair(x, y);45
46 for( set::iterator it = P.begin() ; it != P.end() ; it++)47 {48 if( Matrix[x][*it] != 0)49 {50 //cout<
51 P.erase(*it);52 break;53 }54 }55 P.insert(y);56 //cout<
57 return true;58 }59
60 voidpnt()61 {62 for( size_t i = 0 ; i < Matrix.size() ; i++)63 {64 for( size_t j = 0 ; j < Matrix[0].size() ; j++)65 {66 cout<
73 void Gaussian(pair p)//行變換
74 {75 size_t x =p.first;76 size_t y =p.second;77 double norm =Matrix[x][y];78 for( size_t i = 0 ; i < cn ; i++ )//主行歸一化
79 {80 Matrix[x][i] /=norm;81 }82 for( size_t i = 0 ; i < bn && i != x; i++)83 {84 if( Matrix[i][y] != 0)85 {86 double tmpnorm =Matrix[i][y];87 for( size_t j = 0 ; j < cn ; j++)88 {89 Matrix[i][j] = Matrix[i][j] - tmpnorm *Matrix[x][j];90 }91 }92 }93 }94
95 voidsolve()96 {97 pairt;98 while(1)99 {100
101 pnt();102 if( Pivot(t) == 0)103 {104 return;105 }106 cout<::iterator it = P.begin(); it != P.end() ; it++)108 {109 cout<
116 int main(int argc, char *argv[])117 {118 //ifstream fin;119 //fin.open("./test");
120 cin>>cn>>bn;121 for( size_t i = 0 ; i < bn ; i++)122 {123 vectorvectmp;124 for( size_t j = 0 ; j < cn ; j++)125 {126 double tmp = 0;127 cin>>tmp;128 vectmp.push_back(tmp);129 }130 Matrix.push_back(vectmp);131 }132
133 for( size_t i = 0 ; i < bn-1 ; i++)134 {135 P.insert(cn-i-2);136 }137 solve();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 5157 1 14 6 0 0 0 0 0158 1 1 1 1 0 0 0 4159 1 0 0 0 1 0 0 2160 0 0 1 0 0 1 0 3161 0 3 1 0 0 0 1 6162 */
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。
這樣的矩陣我們稱為全幺模矩陣。
番外篇：
1. 必要不充分：只含1，-1，0。因為單個元素可以看作行列式……
2. 充分必要：對它進行高斯消元的主元操作……(好像叫轉軸？啊反正就是消別人的未知數……)，得來的還是全幺模矩陣……這個是因為那個啥啥幺模矩陣組成了一個乘法群？用這個性質我們可以不用double。
3. 您可以手工屠矩陣用定義證它是全幺模！
4. 如果數學太差，您也可以寫一個O(4^n * n^3)的強程序證它是全幺模！
orzorzorz
附上一道題BZOJ 1061
1 #include
2 #include
3 #include
4 #include
5 #include
6 using namespacestd;7 typedef long longll;8 const int M=10005,N=1005,INF=1e9;9 const double eps=1e-6;10 inline intread(){11 char c=getchar();int x=0,f=1;12 while(c'9'){if(c=='-')f=-1; c=getchar();}13 while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}14 return x*f;15 }16
17 intn,m;18 doublea[M][N],b[M],c[N],v;19 void pivot(int l,inte){20 b[l]/=a[l][e];21 for(int j=1;j<=n;j++) if(j!=e) a[l][j]/=a[l][e];22 a[l][e]=1/a[l][e];23
24 for(int i=1;i<=m;i++) if(i!=l&&fabs(a[i][e])>0){25 b[i]-=a[i][e]*b[l];26 for(int j=1;j<=n;j++) if(j!=e) a[i][j]-=a[i][e]*a[l][j];27 a[i][e]=-a[i][e]*a[l][e];28 }29
30 v+=c[e]*b[l];31 for(int j=1;j<=n;j++) if(j!=e) c[j]-=c[e]*a[l][j];32 c[e]=-c[e]*a[l][e];33
34 //swap(B[l],N[e])
35 }36
37 doublesimplex(){38 while(true){39 int e=0,l=0;40 for(e=1;e<=n;e++) if(c[e]>eps) break;41 if(e==n+1) returnv;42 double mn=INF;43 for(int i=1;i<=m;i++)44 if(a[i][e]>eps&&mn>b[i]/a[i][e]) mn=b[i]/a[i][e],l=i;45 if(mn==INF) return INF;//unbounded
46 pivot(l,e);47 }48 }49
50 intmain(){51 n=read();m=read();52 for(int i=1;i<=n;i++) c[i]=read();53 for(int i=1;i<=m;i++){54 int s=read(),t=read();55 for(int j=s;j<=t;j++) a[i][j]=1;56 b[i]=read();57 }58 printf("%d",(int)(simplex()+0.5));59 }
附上幾道題的題號，練習學習一下：
BZOJ 3550
BZOJ 3112
BZOJ 3265
BZOJ 1061
BZOJ 1061
POJ 1275
標準型:
轉化為標準型:
松弛型:
單純型算法
轉軸：
初始化：
代碼實現：
1 #include
2 #include
3 #include
4 #include
5 #include
6 using namespacestd;7 typedef long longll;8 const int N=25;9 const double eps=1e-8,INF=1e15;10 inline intread(){11 char c=getchar();int x=0,f=1;12 while(c'9'){if(c=='-')f=-1; c=getchar();}13 while(c>='0'&&c<='9'){x=x*10+c-'0'; c=getchar();}14 return x*f;15 }16 intn,m,type;17 doublea[N][N],ans[N];18 int id[N<<1];19 intq[N];20 void Pivot(int l,inte){21 swap(id[n+l],id[e]);22 double t=a[l][e]; a[l][e]=1;23 for(int j=0;j<=n;j++) a[l][j]/=t;24 for(int i=0;i<=m;i++) if(i!=l && abs(a[i][e])>eps){25 t=a[i][e]; a[i][e]=0;26 for(int j=0;j<=n;j++) a[i][j]-=a[l][j]*t;27 }28 }29 boolinit(){30 while(true){31 int e=0,l=0;32 for(int i=1;i<=m;i++) if(a[i][0]eps) {e=j;break;}45 if(!e) break;46 for(int i=1;i<=m;i++) if(a[i][e]>eps && a[i][0]/a[i][e]
對偶性：
全幺模矩陣(totally unimodular matrix)
充分條件：
任何最大流、最小費用最大流的線性規劃都是全幺模矩陣
總結
以上是生活随笔為你收集整理的单纯形法只有两个约束条件_线性规划之单纯形法【超详解+图解】的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Vue中使用echarts绘制地图，以及
- 下一篇： C语言 基础 1