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

时间 2019-11-19

标签线性规划单纯超详解+图解繁體版

原文原文链接

1.做用

单纯形法是解决线性规划问题的一个有效的算法。线性规划就是在一组线性约束条件下，求解目标函数最优解的问题。php

2.线性规划的通常形式

在约束条件下，寻找目标函数z的最大值。ios

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

3.线性规划的可行域

知足线性规划问题约束条件的全部点组成的集合就是线性规划的可行域。若可行域有界（如下主要考虑有界可行域），线性规划问题的目标函数最优解必然在可行域的顶点上达到最优。c++

单纯形法就是经过设置不一样的基向量，通过矩阵的线性变换，求得基可行解（可行域顶点），并判断该解是否最优，不然继续设置另外一组基向量，重复执行以上步骤，直到找到最优解。因此，单纯形法的求解过程是一个循环迭代的过程。算法

图1 可行域数据结构

4.线性规划的标准形式

在说明单纯形法的原理以前，须要明白线性规划的标准形式。由于单纯形算法是经过线性规划的标准形来求解的。通常，规定线性规划的标准形式为：函数

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

写成矩阵形式：学习

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

标准形的形式为：优化

1）目标函数要求maxspa

2）约束条件均为等式.net

3）决策变量为非负约束

普通线性规划化为标准形：

1）若目标函数为最小化，能够经过取负，求最大化

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

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

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

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

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

本文最开始的线性规划问题转化为标准形为：

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

5.单纯形法

5.1几何意义

在标准形中，有m个约束条件（不包括非负约束），n个决策变量，且（n>=m）。首先，选取m个基变量，基变量对应约束系数矩阵的列向量线性无关。经过矩阵的线性变换，基变量可由非基变量表示：

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

若是令非基变量等于0，可求得基变量的值：

$x_i' = {C_i}$

若是为可行解的话，Ci大于0。那么它的几何意义是什么呢？

仍是经过上述具体的线性规划问题来讲明。

若是选择x二、x3为基变量，那么令x一、x4等于0，能够去求解基变量x二、x3的值。对系数矩阵作行变换，以下所示，x2=9/2，x3=15/2

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

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

图2

因此，经过选择不一样的基变量，能够得到不一样的可行域的顶点。

5.2如何判断最优

如前所述，基变量可由非基变量表示：

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

目标函数z也能够彻底由非基变量表示：

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

当达到最优解时，全部的 ${{\sigma _j}}$ 应小于等于0。当存在j， ${{\sigma _j}}$ >0时，当前解不是最优解，为何？

当前的目标函数值为z0，其中全部的非基变量值均取0。由以前分析可知， ${x_j'}$ =0表明可行域的某个边界，是 ${x_j'}$ 的最小值。若是可行解逐步离开这个边界， ${x_j'}$ 会变大，由于 ${{\sigma _j}}$ >0，显然目标函数的取值也会变大，因此当前解不是最优解。咱们须要寻找新的基变量。

5.3如何选择新的基变量

若是存在多个 ${{\sigma _j}}$ >0，选择最大的 ${{\sigma _j}}$ >0对应的变量做为基变量，这表示目标函数随着 ${x_j'}$ 的增长，增加的最快。

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

假如咱们选择非基变量做为下一轮的基变量，那么被替换基变量 ${x_j'}$ 在下一轮中做为非基变量，等于0。选择 ${x_j'}$ 的原则：替换后应该尽可能使值最大（由于上面已分析过，目标函数会随着的增大而增大）。

继续经过上面的例子来讲明：

从最后一行能够看到，x1的系数为1/2>0，因此选x二、x3为基变量并无是目标函数达到最优。下一轮选取x1做为基变量，替换x二、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

咱们能够获得其松弛形式：

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 <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);
 56     //cout<<"add "<<y<<endl;
 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的向量

最大化 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<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;
10 inline int read(){
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(){
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

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;
10 inline int read(){
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);
56     n=read();m=read();type=read();
57     for(int i=1;i<=n;i++) a[0][i]=read();
58     for(int i=1;i<=m;i++){
59         for(int j=1;j<=n;j++) a[i][j]=read();
60         a[i][0]=read();
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)

充分条件：

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