作者：jostree 转载请注明出处 http://www.cnblogs.com/jostree/p/4156685.html
使用单纯型法来求解线性规划，输入单纯型法的松弛形式，是一个大矩阵，第一行为目标函数的系数，且最后一个数字为当前轴值下的 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.
整个过程代码如下所示：
#include <cstdio>
#include <climits>
#include <cstdlib>
#include <iostream>
#include <cstring>
#include <vector>
#include <fstream>
#include <set>
using namespace std;
vector<vector<double> > Matrix;
double Z;
set<int> P;
size_t cn, bn; bool Pivot(pair<size_t, size_t> &p)//返回0表示所有的非轴元素都小于0
{
int x = , y = ;
double cmax = -INT_MAX;
vector<double> C = Matrix[];
vector<double> B; for( size_t i = ; i < bn ; i++ )
{
B.push_back(Matrix[i][cn-]);
} for( size_t i = ; i < C.size(); i++ )//在非轴元素中找最大的c
{
if( cmax < C[i] && P.find(i) == P.end())
{
cmax = C[i];
y = i;
}
}
if( cmax < )
{
return ;
} double bmin = INT_MAX;
for( size_t i = ; i < bn ; i++ )
{
double tmp = B[i]/Matrix[i][y];
if( Matrix[i][y] != && bmin > tmp )
{
bmin = tmp;
x = i;
}
} p = make_pair(x, y); for( set<int>::iterator it = P.begin() ; it != P.end() ; it++)
{
if( Matrix[x][*it] != )
{
//cout<<"erase "<<*it<<endl;
P.erase(*it);
break;
}
}
P.insert(y);
//cout<<"add "<<y<<endl;
return true;
} void pnt()
{
for( size_t i = ; i < Matrix.size() ; i++ )
{
for( size_t j = ; j < Matrix[].size() ; j++ )
{
cout<<Matrix[i][j]<<"\t";
}
cout<<endl;
}
cout<<"result z:"<<-Matrix[][cn-]<<endl;
} void Gaussian(pair<size_t, size_t> p)//行变换
{
size_t x = p.first;
size_t y = p.second;
double norm = Matrix[x][y];
for( size_t i = ; i < cn ; i++ )//主行归一化
{
Matrix[x][i] /= norm;
}
for( size_t i = ; i < bn && i != x; i++ )
{
if( Matrix[i][y] != )
{
double tmpnorm = Matrix[i][y];
for( size_t j = ; j < cn ; j++ )
{
Matrix[i][j] = Matrix[i][j] - tmpnorm * Matrix[x][j];
}
}
}
} void solve()
{
pair<size_t, size_t> t;
while()
{ pnt();
if( Pivot(t) == )
{
return;
}
cout<<t.first<<" "<<t.second<<endl;
for( set<int>::iterator it = P.begin(); it != P.end() ; it++ )
{
cout<<*it<<" ";
}
cout<<endl;
Gaussian(t);
}
} int main(int argc, char *argv[])
{
ifstream fin;
fin.open("./test");
fin>>cn>>bn;
for( size_t i = ; i < bn ; i++ )
{
vector<double> vectmp;
for( size_t j = ; j < cn ; j++)
{
double tmp = ;
fin>>tmp;
vectmp.push_back(tmp);
}
Matrix.push_back(vectmp);
} for( size_t i = ; i < bn- ; i++ )
{
P.insert(cn-i-);
}
solve();
}
/////////////////////////////////////
//glpk input:
///* Variables */
//var x1 >= 0;
//var x2 >= 0;
//var x3 >= 0;
///* Object function */
//maximize z: x1 + 14*x2 + 6*x3;
///* Constrains */
//s.t. con1: x1 + x2 + x3 <= 4;
//s.t. con2: x1 <= 2;
//s.t. con3: x3 <= 3;
//s.t. con4: 3*x2 + x3 <= 6;
//end;
/////////////////////////////////////
//myinput:
//8 5
//1 14 6 0 0 0 0 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
/////////////////////////////////////