单纯形法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

我们可以构造单纯性表，其中最后一行打星的列为轴值。

在单纯性表中，我们发现非轴值的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出轴的过程。变换后的单纯性表为：

继续计算，我们得到：

此时我们发现，所有非轴的x的系数全部小于零，即增大任何非轴的x值并不能使得目标函数最大，从而得到最优解32.

整个过程代码如下所示：

  1 #include <cstdio>
  2 #include <climits>
  3 #include <cstdlib>
  4 #include <iostream>
  5 #include <cstring>
  6 #include  <vector>
  7 #include  <fstream>
  8 #include <set>
  9 using namespace std;
 10 vector<vector<double> > Matrix;
 11 double Z;
 12 set<int> P;
 13 size_t cn, bn;
 14
 15 bool Pivot(pair<size_t, size_t> &p)//返回0表示所有的非轴元素都小于0
 16 {
 17     int x = 0, y = 0;
 18     double cmax = -INT_MAX;
 19     vector<double> C = Matrix[0];
 20     vector<double> B;
 21
 22     for( size_t i = 0 ; i < bn ; i++ )
 23     {
 24         B.push_back(Matrix[i][cn-1]);
 25     }
 26
 27     for( size_t i = 0 ; i < C.size(); i++ )//在非轴元素中找最大的c
 28     {
 29         if( cmax < C[i] && P.find(i) == P.end())
 30         {
 31             cmax = C[i];
 32             y = i;
 33         }
 34     }
 35     if( cmax < 0 )
 36     {
 37         return 0;
 38     }
 39
 40     double bmin = INT_MAX;
 41     for( size_t i = 1 ; i < bn ; i++ )
 42     {
 43         double tmp = B[i]/Matrix[i][y];
 44        if( Matrix[i][y] != 0 && bmin > tmp )
 45        {
 46            bmin = tmp;
 47            x = i;
 48        }
 49     }
 50
 51     p = make_pair(x, y);
 52
 53     for( set<int>::iterator it = P.begin() ; it != P.end() ; it++)
 54     {
 55         if( Matrix[x][*it] != 0 )
 56         {
 57             //cout<<"erase "<<*it<<endl;
 58             P.erase(*it);
 59             break;
 60         }
 61     }
 62     P.insert(y);
 63     //cout<<"add "<<y<<endl;
 64     return true;
 65 }
 66
 67 void pnt()
 68 {
 69     for( size_t i = 0 ; i < Matrix.size() ; i++ )
 70     {
 71         for( size_t j = 0 ; j < Matrix[0].size() ; j++ )
 72         {
 73             cout<<Matrix[i][j]<<"\t";
 74         }
 75     cout<<endl;
 76     }
 77     cout<<"result z:"<<-Matrix[0][cn-1]<<endl;
 78 }
 79
 80 void Gaussian(pair<size_t, size_t> p)//行变换
 81 {
 82     size_t  x = p.first;
 83     size_t y = p.second;
 84     double norm = Matrix[x][y];
 85     for( size_t i = 0 ; i < cn ; i++ )//主行归一化
 86     {
 87         Matrix[x][i] /= norm;
 88     }
 89     for( size_t i = 0 ; i < bn && i != x; i++ )
 90     {
 91         if( Matrix[i][y] != 0)
 92         {
 93             double tmpnorm = Matrix[i][y];
 94             for( size_t j = 0 ; j < cn ; j++ )
 95             {
 96                 Matrix[i][j] = Matrix[i][j] - tmpnorm * Matrix[x][j];
 97             }
 98         }
 99     }
100 }
101
102 void solve()
103 {
104     pair<size_t, size_t> t;
105     while(1)
106     {
107
108         pnt();
109         if( Pivot(t) == 0 )
110         {
111             return;
112         }
113         cout<<t.first<<" "<<t.second<<endl;
114         for( set<int>::iterator it = P.begin(); it != P.end()  ; it++ )
115         {
116             cout<<*it<<" ";
117         }
118         cout<<endl;
119         Gaussian(t);
120     }
121 }
122
123 int main(int argc, char *argv[])
124 {
125     ifstream fin;
126     fin.open("./test");
127     fin>>cn>>bn;
128     for( size_t i = 0 ; i < bn ; i++ )
129     {
130         vector<double> vectmp;
131         for( size_t j = 0 ; j < cn ; j++)
132         {
133             double tmp = 0;
134             fin>>tmp;
135             vectmp.push_back(tmp);
136         }
137         Matrix.push_back(vectmp);
138     }
139
140     for( size_t i = 0 ; i < bn-1 ; i++ )
141     {
142         P.insert(cn-i-2);
143     }
144     solve();
145 }
146 /////////////////////////////////////
147 //glpk input:
148 ///* Variables */
149 //var x1 >= 0;
150 //var x2 >= 0;
151 //var x3 >= 0;
152 ///* Object function */
153 //maximize z: x1 + 14*x2 + 6*x3;
154 ///* Constrains */
155 //s.t. con1: x1 + x2 + x3 <= 4;
156 //s.t. con2: x1  <= 2;
157 //s.t. con3: x3  <= 3;
158 //s.t. con4: 3*x2 + x3  <= 6;
159 //end;
160 /////////////////////////////////////
161 //myinput:
162 //8 5
163 //1 14 6 0 0 0 0 0
164 //1 1 1 1 0 0 0 4
165 //1 0 0 0 1 0 0 2
166 //0 0 1 0 0 1 0 3
167 //0 3 1 0 0 0 1 6
168 /////////////////////////////////////