单纯形法只有两个约束条件_线性规划之单纯形法【超详解+图解】

使用单纯型法来求解线性规划，输入单纯型法的松弛形式，是一个大矩阵，第一行为目标函数的系数，且最后一个数字为当前轴值下的 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)

充分条件：

任何最大流、最小费用最大流的线性规划都是全幺模矩阵