本程序实现的是单纯型法解线性规划问题,在检验向量e<=0的情况下,迭代得到最优解。
首先,要求你在data.txt文件中输入数据,输入方式,例:
3
1 1 1
2 3
3 1 1 >= 1
-1 4 1 >= 2
便可出现执行结果!
实现程序如下:
- #include<iostream>
- #include<fstream>
- #include<string>
- #include<cstdlib>
- #define INF -100000
- #define MAX 100000
- using namespace std;
- int N1, N2;
- //约束条件系数
- double **a;
- double *b; // 向量b
- double *e;//检验数向量
- double RHS = 0;//目标函数值
- int *x; // 存放基本可行解
- int m; // 添加的辅助变量的个数
- int JudgeSym(string symbol)
- {
- if(symbol==">=")
- return -1;
- if(symbol=="<=")
- return 1;
- if(symbol=="=")
- return 0;
- }
- // 找出初值解
- void Begin()
- {
- //进行列搜索, 得到初始解
- int index = 0;
- for(int j=1; j<N2+m; ++j){
- int cnt = 0;
- for(int i=1; i<N1; ++i)
- if(a[i][j]){
- ++cnt;
- index = i;
- }
- if(cnt==1&&x[index]==0){
- if(a[index][j] != 1){
- for(int k=1; k<N2; ++k)
- if(k!=j)
- a[index][k] /= a[index][j];
- b[index] /= a[index][j];
- a[index][j] = 1;
- }
- x[index] = j;
- }
- }
- cout<<"得到对偶问题的初始基本可行解为:"<<endl;
- for(int i=1; i<N1; ++i)
- cout<<"x("<<x[i]<<")"<<" = "<<b[i]/a[i][x[i]]<<endl;
- }
- void GetFirstForm()
- {
- double tmp = 0;
- for(int j=1; j<N1; ++j)
- if(e[x[j]]!=0){
- tmp = e[x[j]];
- for(int i=1; i<N2; ++i){
- /*May Have Error */
- e[i] -= tmp*a[j][i];
- }
- RHS -= tmp*b[j];
- }
- }
- int CheckRvec()
- {
- double min = MAX;
- int index = 0;
- for(int i=1; i<N1; ++i)
- if(b[i]<0&&b[i]<min){
- min =b[i];
- index = i;
- }
- if(index!=0)
- return index;
- cout<<"迭代得到一个最优解:RHS = "<<RHS<<endl;
- for(int i=1; i<N1; ++i)
- cout<<"x("<<x[i]<<")"<<" = "<<b[i]/a[i][x[i]]<<endl;
- return -1;
- }
- int MinEj(int Xiout)
- {
- double min = MAX;
- int index = -1;
- for(int j=1; j<N2+m; ++j)
- if(a[Xiout][j]<0&&e[j]/a[Xiout][j]<min){
- index = j;
- min = e[j]/a[Xiout][j];
- }
- if(index == -1)
- cout<<"a(r,j)>=0;无最优解!"<<endl;
- return index;
- }
- void InputandInit()
- {
- void Standardlize(int,int);
- ifstream infile("data.txt");
- infile>>N2;
- N2 += 1;
- //存储没有添加辅助变量的
- double *tpe = new double[N2];
- for(int i=1; i<N2; ++i){
- infile>>tpe[i];
- tpe[i] = -tpe[i];
- }/*修改,根据BeginOne()的返回值,判断是否进行e的再次录入*/
- infile>>N1>>N2;
- m = N1; /*最多添加的辅助变量的个数*/
- N1 += 1;
- N2 += 1;
- //正式进行目标函数的录入
- e = new double[N2+m];
- memset(e,0,sizeof(double)*(N2+m));
- memcpy(e,tpe,sizeof(double)*N2);
- delete tpe;/**/
- //为 系数矩阵申请内存
- a = new double*[N1];
- for(int i=1; i<N1; ++i){
- a[i] = new double[N2+m]; /*a : N1/N2+N1 */
- memset(a[i],0,sizeof(double)*(N2+m));
- }
- //为 右端向量申请内存
- b = new double[N1];
- memset(b,0,sizeof(double)*N1);
- m = 0; /*申请内存结束后,重新置m为零,使m记录实际添加的人工变量的个数*/
- //输入a系数矩阵
- for(int i=1; i<N1; ++i){
- for(int j=1; j<N2; ++j){
- infile>>a[i][j];
- }
- string Symbol;
- infile>>Symbol;
- infile>>b[i];
- int numSym = JudgeSym(Symbol);
- Standardlize(numSym,i);
- }
- }
- void Standardlize(int numSym, int line)
- {
- int row = N1 -1;
- int col = N2 -1;
- static int cnt = 1;
- switch(numSym){
- case -1:
- for(int j=1; j<N2; ++j)
- a[line][j] = -a[line][j];
- b[line] = -b[line];
- a[line][col+cnt] = 1;/*添加辅助变量*/
- ++cnt;
- ++m; /*辅助变量加1*/
- break;
- case 1:
- a[line][col+cnt] = 1;
- ++cnt;
- ++m; /*辅助变量加1*/
- break;
- default: return;
- }
- return;
- }
- void print()
- {
- static int cnt = 0;
- cout<<left<<showpos;
- cout.width(5);
- cout.precision(3);
- cout<<cnt<<"-Z ";
- for(int i=1; i<N2+m; ++i)
- cout<<e[i]<<" ";
- cout<<RHS<<endl;
- for(int i=1; i<N1; ++i){
- cout<<"x("<<x[i]<<") ";
- for(int j=1; j<N2+m; ++j)
- cout<<a[i][j]<<" ";
- cout<<b[i]<<endl;
- }
- cout<<"***********************"<<endl;
- ++cnt;
- }
- int main()
- {
- InputandInit();
- /* Until then; data input finished */
- x = new int[N1];
- memset(x,0,sizeof(int)*N1);
- Begin();
- GetFirstForm();
- //打印第一张单纯行表
- cout<<"The First Form Of Simplex Method:"<<endl;
- print();
- /*循环得出最优解*/
- int Xiout = 0, ej = 0; /*找到最小的e(j)/a(Xiout)(j)*/
- while((Xiout=CheckRvec())!=-1 &&(ej = MinEj(Xiout))!=-1){
- x[Xiout] = ej;
- for(int i=1; i<N2+m; ++i)
- if(ej!=i){
- a[Xiout][i] = a[Xiout][i]/a[Xiout][ej];
- }
- /*Error will be found*/
- b[Xiout]/=a[Xiout][ej];
- a[Xiout][ej] = 1;
- //化检验数向量
- double tmp = e[ej];
- for(int i=1; i<N2+m; ++i)
- e[i] -= tmp*a[Xiout][i];
- RHS -= tmp*b[Xiout];
- /*化右端向量*/
- for(int i=1; i<N1; ++i){
- if(i!=Xiout){
- tmp = a[i][ej];
- for(int j=1; j<N2+m; ++j)
- a[i][j] -= tmp*a[Xiout][j];
- b[i] -= tmp*b[Xiout];
- }
- }
- print();
- }
- system("pause");
- return 0;
- }
本文出自 “东方快翔” 博客,请务必保留此出处http://hustluy.blog.51cto.com/1792080/412999