首页 > 代码库 > 单纯形法C++实现
单纯形法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.
整个过程代码如下所示：
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 5163 //1 14 6 0 0 0 0 0164 //1 1 1 1 0 0 0 4165 //1 0 0 0 1 0 0 2166 //0 0 1 0 0 1 0 3167 //0 3 1 0 0 0 1 6168 /////////////////////////////////////
单纯形法C++实现