單純形法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]<<"\t";
        }
    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
/////////////////////////////////////