单纯形法C++实现

作者：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 = 0, y = 0;
    double cmax = -INT_MAX;
    vector<double> C = Matrix[0];
    vector<double> B;

    for( size_t i = 0 ; i < bn ; i++ )
    {
        B.push_back(Matrix[i][cn-1]);
    }
    
    for( size_t i = 0 ; i < C.size(); i++ )//在非轴元素中找最大的c
    {
        if( cmax < C[i] && P.find(i) == P.end())
        {
            cmax = C[i];
            y = i;
        }
    }
    if( cmax < 0 )
    {
        return 0;
    }

    double bmin = INT_MAX;
    for( size_t i = 1 ; i < bn ; i++ )
    {
        double tmp = B[i]/Matrix[i][y];
       if( Matrix[i][y] != 0 && 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] != 0 )
        {
            //cout<<"erase "<<*it<<endl;
            P.erase(*it);
            break;
        }
    }
    P.insert(y);
    //cout<<"add "<<y<<endl;
    return true;
}

void pnt()
{
    for( size_t i = 0 ; i < Matrix.size() ; i++ )
    {
        for( size_t j = 0 ; j < Matrix[0].size() ; j++ )
        {
            cout<<Matrix[i][j]<<"	";
        }
    cout<<endl;
    }
    cout<<"result z:"<<-Matrix[0][cn-1]<<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 = 0 ; i < cn ; i++ )//主行归一化
    {
        Matrix[x][i] /= norm;
    }
    for( size_t i = 0 ; i < bn && i != x; i++ )
    {
        if( Matrix[i][y] != 0)
        {
            double tmpnorm = Matrix[i][y];
            for( size_t j = 0 ; j < cn ; j++ )
            {
                Matrix[i][j] = Matrix[i][j] - tmpnorm * Matrix[x][j];
            }
        }
    }
}

void solve()
{
    pair<size_t, size_t> t;
    while(1)
    {

        pnt();
        if( Pivot(t) == 0 )
        {
            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 = 0 ; i < bn ; i++ )
    {
        vector<double> vectmp;
        for( size_t j = 0 ; j < cn ; j++)
        {
            double tmp = 0;
            fin>>tmp;
            vectmp.push_back(tmp);
        }
        Matrix.push_back(vectmp);
    }

    for( size_t i = 0 ; i < bn-1 ; i++ )
    {
        P.insert(cn-i-2);
    }
    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
/////////////////////////////////////