# 线性规划之单纯形法【超详解+图解】

## 2.线性规划的一般形式

$\begin{array}{l}&space;\max&space;z&space;=&space;{x_1}&space;+&space;{x_2}\\&space;s.t.\left\{&space;{\begin{array}{*{20}{c}}&space;{2{x_1}&space;+&space;{x_2}&space;\le&space;12}\\&space;{{x_1}&space;+&space;2{x_2}&space;\le&space;9}\\&space;{{x_1},{x_2}&space;\ge&space;0}&space;\end{array}}&space;\right.&space;\end{array}$

# 4.线性规划的标准形式

$\begin{array}{l}&space;\max&space;z&space;=&space;\sum\limits_{j&space;=&space;1}^n&space;{{c_j}{x_j}}&space;\\&space;s.t.\left\{&space;{\begin{array}{*{20}{c}}&space;{\sum\limits_{j&space;=&space;1}^n&space;{{a_{ij}}{x_j}&space;=&space;{b_j},&space;i&space;=&space;1,2,...,m}&space;}\\&space;{xj&space;\ge&space;0,&space;j&space;=&space;1,2,...,n}&space;\end{array}}&space;\right.&space;\end{array}$

$\begin{array}{l}&space;\max&space;z&space;=&space;CX\\&space;AX&space;=&space;b\\&space;X&space;\ge&space;0\\&space;A&space;=&space;\left[&space;{\begin{array}{*{20}{c}}&space;{{a_{11}}}&{{a_{12}}}&&space;\cdots&space;&{{a_{1n}}}\\&space;\vdots&space;&&space;\vdots&space;&&space;\ddots&space;&&space;\vdots&space;\\&space;{{a_{m1}}}&{{a_{m2}}}&&space;\cdots&space;&{{a_{mn}}}&space;\end{array}}&space;\right]&space;\end{array}$

1）目标函数要求max

2）约束条件均为等式

3）决策变量为非负约束

1）若目标函数为最小化，可以通过取负，求最大化

2）约束不等式为小于等于不等式，可以在左端加入非负松弛变量，转变为等式，比如：

${x_1}&space;+&space;2{x_2}&space;\le&space;9&space;\Rightarrow&space;\left\{&space;{\begin{array}{*{20}{c}}&space;{{x_1}&space;+&space;2{x_2}&space;+&space;{x_3}&space;=&space;9}\\&space;{{x_3}&space;\ge&space;0}&space;\end{array}}&space;\right.$

同理，约束不等式为大于等于不等式时，可以在左端减去一个非负松弛变量，变为等式。

3）若存在取值无约束的变量，可转变为两个非负变量的差，比如：

$-&space;\infty&space;\le&space;{x_k}&space;\le&space;+&space;\infty&space;\Rightarrow&space;\left\{&space;{\begin{array}{*{20}{c}}&space;{{x_k}&space;=&space;{x_m}&space;-&space;{x_n}}\\&space;{{x_m},{x_n}&space;\ge&space;0}&space;\end{array}}&space;\right.$

$\begin{array}{l}&space;\max&space;z&space;=&space;{x_1}&space;+&space;{x_2}\\&space;s.t.\left\{&space;{\begin{array}{*{20}{c}}&space;{2{x_1}&space;+&space;{x_2}&space;+&space;{x_3}&space;=&space;12}\\&space;{{x_1}&space;+&space;2{x_2}&space;+&space;{x_4}&space;=&space;9}\\&space;{{x_1},{x_2},{x_3},{x_4}&space;\ge&space;0}&space;\end{array}}&space;\right.&space;\end{array}$

# 5.单纯形法

## 5.1几何意义

$\inline&space;x_i'&space;=&space;{C_i}&space;+&space;\sum\limits_{j&space;=&space;m&space;+&space;1}^n&space;{{m_{ij}}x_j'}&space;(i&space;=&space;1,2,...,m)$

$\inline&space;x_i'&space;=&space;{C_i}$

$\inline&space;\begin{array}{l}&space;\max&space;z&space;=&space;{x_1}&space;+&space;{x_2}\\&space;s.t.\left\{&space;{\begin{array}{*{20}{c}}&space;{2{x_1}&space;+&space;{x_2}&space;+&space;{x_3}&space;=&space;12}\\&space;{{x_1}&space;+&space;2{x_2}&space;+&space;{x_4}&space;=&space;9}\\&space;{{x_1},{x_2},{x_3},{x_4}&space;\ge&space;0}&space;\end{array}}&space;\right.&space;\end{array}$

$\inline&space;\left[&space;{\begin{array}{*{20}{c}}&space;{\rm{X}}&{{x_1}}&{{x_2}}&{{x_3}}&{{x_4}}&b\\&space;{}&2&1&1&0&{12}\\&space;{}&1&2&0&1&9\\&space;{\rm{C}}&1&1&0&0&z&space;\end{array}}&space;\right]&space;\to&space;\left[&space;{\begin{array}{*{20}{c}}&space;{\rm{X}}&{{x_1}}&{{x_2}}&{{x_3}}&{{x_4}}&b\\&space;{}&{{\textstyle{3&space;\over&space;2}}}&0&1&{&space;-&space;{\textstyle{1&space;\over&space;2}}}&{{\textstyle{{15}&space;\over&space;2}}}\\&space;{}&{{\textstyle{1&space;\over&space;2}}}&1&0&{{\textstyle{1&space;\over&space;2}}}&{{\textstyle{9&space;\over&space;2}}}\\&space;{\rm{C}}&{{\textstyle{1&space;\over&space;2}}}&0&0&{&space;-&space;{\textstyle{1&space;\over&space;2}}}&{z&space;-&space;{\textstyle{9&space;\over&space;2}}}&space;\end{array}}&space;\right]$

X1=0表示可行解在x轴上；X4=0表示可行解在x1+2x2=9的直线上。那么，求得的可行解即表示这两条直线的交点，也是可行域的顶点，如图所示：

所以，通过选择不同的基变量，可以获得不同的可行域的顶点。

## 5.2如何判断最优

$\inline&space;x_i'&space;=&space;{C_i}&space;+&space;\sum\limits_{j&space;=&space;m&space;+&space;1}^n&space;{{m_{ij}}x_j'}&space;(i&space;=&space;1,2,...,m)$

$\inline&space;z&space;=&space;{z_0}&space;+&space;\sum\limits_{j&space;=&space;m&space;+&space;1}^n&space;{{\sigma&space;_j}x_j'}$

## 5.4如何选择被替换的基变量

$\inline&space;\left[&space;{\begin{array}{*{20}{c}}&space;{\rm{X}}&{{x_1}}&{{x_2}}&{{x_3}}&{{x_4}}&b\\&space;{}&2&1&1&0&{12}\\&space;{}&1&2&0&1&9\\&space;{\rm{C}}&1&1&0&0&z&space;\end{array}}&space;\right]&space;\to&space;\left[&space;{\begin{array}{*{20}{c}}&space;{\rm{X}}&{{x_1}}&{{x_2}}&{{x_3}}&{{x_4}}&b\\&space;{}&{{\textstyle{3&space;\over&space;2}}}&0&1&{&space;-&space;{\textstyle{1&space;\over&space;2}}}&{{\textstyle{{15}&space;\over&space;2}}}\\&space;{}&{{\textstyle{1&space;\over&space;2}}}&1&0&{{\textstyle{1&space;\over&space;2}}}&{{\textstyle{9&space;\over&space;2}}}\\&space;{\rm{C}}&{{\textstyle{1&space;\over&space;2}}}&0&0&{&space;-&space;{\textstyle{1&space;\over&space;2}}}&{z&space;-&space;{\textstyle{9&space;\over&space;2}}}&space;\end{array}}&space;\right]$

## 5.5终止条件

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 * * * *

 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 * * * *

  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 {
10     int x = 0, y = 0;
11     double cmax = -INT_MAX;
12     vector<double> C = Matrix[0];
13     vector<double> B;
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<int>::iterator it = P.begin() ; it != P.end() ; it++)
47     {
48         if( Matrix[x][*it] != 0 )
49         {
50             //cout<<"erase "<<*it<<endl;
51             P.erase(*it);
52             break;
53         }
54     }
55     P.insert(y);
57     return true;
58 }
59
60 void pnt()
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<<Matrix[i][j]<<"\t";
67         }
68     cout<<endl;
69     }
70     cout<<"result z:"<<-Matrix[0][cn-1]<<endl;
71 }
72
73 void Gaussian(pair<size_t, size_t> 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 void solve()
96 {
97     pair<size_t, size_t> t;
98     while(1)
99     {
100
101         pnt();
102         if( Pivot(t) == 0 )
103         {
104             return;
105         }
106         cout<<t.first<<" "<<t.second<<endl;
107         for( set<int>::iterator it = P.begin(); it != P.end()  ; it++ )
108         {
109             cout<<*it<<" ";
110         }
111         cout<<endl;
112         Gaussian(t);
113     }
114 }
115
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         vector<double> vectmp;
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 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的向量

2.松弛型

xn+i=bi-sigma{aijxj}>=0

3.替入变量 xe(非基变量)

替出变量 xl(基本变量)

4.可行解

基本解：所有非基变量设为0

基本可行解

5.单纯形法的过程中B和N不断交换，在n维空间中不断走

“相当于不等式上的高斯消元”

【代码实现】：

pivot是转动操作

（一个优化 a[i][e]为0的约束没必要带入了

simplex是主过程

c[e]>eps

【对偶原理】：

1.原始线性规划 对偶线性规划

2.对于

min c x^T
Ax = b
x >= 0

x^T表示转置

det(A)表示A的行列式
A_i表示把A的第i列换为b之后的矩阵

1. 必要不充分：只含1，-1，0。因为单个元素可以看作行列式……
2. 充分必要：对它进行高斯消元的主元操作……（好像叫转轴？啊反正就是消别人的未知数……），得来的还是全幺模矩阵……这个是因为那个啥啥幺模矩阵组成了一个乘法群？用这个性质我们可以不用double。
3. 您可以手工屠矩阵用定义证它是全幺模！
4. 如果数学太差，您也可以写一个O(4^n * n^3)的强程序证它是全幺模！

orzorzorz

 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 M=10005,N=1005,INF=1e9;
9 const double eps=1e-6;
11     char c=getchar();int x=0,f=1;
12     while(c<'0'||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 int n,m;
18 double a[M][N],b[M],c[N],v;
19 void pivot(int l,int e){
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 double simplex(){
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) return v;
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 int main(){
53     for(int i=1;i<=m;i++){
55         for(int j=s;j<=t;j++) a[i][j]=1;
57     }
58     printf("%d",(int)(simplex()+0.5));
59 }

BZOJ 3550

BZOJ 3112

BZOJ 3265

BZOJ 1061

BZOJ 1061

POJ 1275

# 代码实现：

UOJ#179. 线性规划

 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;
11     char c=getchar();int x=0,f=1;
12     while(c<'0'||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 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){
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 bool init(){
30     while(true){
31         int e=0,l=0;
32         for(int i=1;i<=m;i++) if(a[i][0]<-eps && (!l||(rand()&1))) l=i;
33         if(!l) break;
34         for(int j=1;j<=n;j++) if(a[l][j]<-eps && (!e||(rand()&1))) e=j;
35         if(!e) {puts("Infeasible");return false;}
36         Pivot(l,e);
37     }
38     return true;
39 }
40 bool simplex(){
41     while(true){
42         int l=0,e=0; double mn=INF;
43         for(int j=1;j<=n;j++)
44             if(a[0][j]>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]<mn)
47             mn=a[i][0]/a[i][e],l=i;
48         if(!l) {puts("Unbounded");return false;}
49         Pivot(l,e);
50     }
51     return true;
52 }
53 int main(){
54     freopen("in","r",stdin);
55     srand(317);
58     for(int i=1;i<=m;i++){
61     }
62     for(int i=1;i<=n;i++) id[i]=i;
63     if(init() && simplex()){
64         printf("%.8lf\n",-a[0][0]);
65         if(type){
66             for(int i=1;i<=m;i++) ans[id[n+i]]=a[i][0];
67             for(int i=1;i<=n;i++) printf("%.8lf ",ans[i]);
68         }
69     }
70 }

# 全幺模矩阵(totally unimodular matrix)

posted @ 2017-06-30 11:00  Angel_Kitty  阅读(121138)  评论(33编辑  收藏  举报