线性规划（算法导论）

花了将近四天的时间初略的把<<算法导论>>线性规划的一章的伪代码实现了，调试的时间略长，当中遇到很多的问题，只有自己动手才会体会的更深。虽然调试出来了，但代码还很不精简，以后有时间来优化一下。

       

#include<iostream>
#include<algorithm>
#include<numeric>
#include<vector>
#include<set>
//#include<math.h>
using namespace std;
#define N 4 //the count of the non basic variables
#define M 3 //the count of the inequations

bool flag=0;

class MODULE
{
public:
	MODULE();
	MODULE & operator=(MODULE &m);
	float ma[M+1][N+M+1];
	float mb[M+1];
	float mc[N+M+1];//目标函数非基本变量对应的系数
	int l[M+1];//每行对应的基本变量下标
	int r[N+1];//非基本变量下标
	float v;
};
MODULE MM;
MODULE::MODULE()
{
	for(int i=1;i<=M;i++)
		for(int j=1;j<=N+M;j++)
			ma[i][j]=0;
	for(int i=1;i<=N+M;i++){
		mc[i]=0;
		if(i<=N)
			r[i]=i;
	}
	for(int i=1;i<=M;i++){
		mb[i]=0;
		l[i]=(i+N);
	}
		
	v=0;
}
MODULE & MODULE::operator=(MODULE &m){
	for(int i=1;i<=M;i++)
		for(int j=0;j<=N+M;j++)
			ma[i][j]=m.ma[i][j];
	for(int i=0;i<=N+M;i++)
		mc[i]=m.mc[i];
		
	for(int i=1;i<=M;i++){
		mb[i]=m.mb[i];l[i]=m.l[i];
	}
	for(int i=0;i<=N;i++)
		r[i]=m.r[i];
	
		
	v=m.v;
	return *this;
}


int minL(float a[],int n)
{
	int s=a[1];int m=1;
	for(int i=1;i<=n;i++)
		if(a[i]<s)
		{s=a[i];m=i;}
	return m;
}

int find(int i,int a[],int n){
	for(int j=1;j<=n;j++)
		if(i==a[j])
			return j;
	return 0;
}

MODULE PIVOT(MODULE m,int l,int e){//将l行对应的基本变量换为非基本变量，将xe换为基本变量
	m.mb[l]=m.mb[l]/m.ma[l][e]; //更新l行的参数
	for(int j=0;j<=N;j++){
		if(m.r[j]==e) continue;
		else{
			m.ma[l][m.r[j]]=m.ma[l][m.r[j]]/m.ma[l][e];
		}
	}
	m.ma[l][m.l[l]]=1/m.ma[l][e];
	m.ma[l][e]=0;
   //compute the coefficients of the remaining constraints
	for(int i=1;i<=M;i++){
		if(i==l) continue;
		else{
			m.mb[i]=m.mb[i]-m.ma[i][e]*m.mb[l];
		}
		for(int j=0;j<=N;j++){
		if(m.r[j]==e) continue;
		else{
			m.ma[i][m.r[j]]=m.ma[i][m.r[j]]-m.ma[i][e]*m.ma[l][m.r[j]];
		}
	}
		m.ma[i][m.l[l]]=0-m.ma[i][e]*m.ma[l][m.l[l]];
		m.ma[i][e]=0;
	}
	//compute the objective function.
	m.v=m.v+m.mc[e]*m.mb[l];
	for(int j=0;j<=N;j++){
		if(m.r[j]==e) continue;
		else{
			m.mc[m.r[j]]=m.mc[m.r[j]]-m.mc[e]*m.ma[l][m.r[j]];
		}
	}
	m.mc[m.l[l]]=0-m.mc[e]*m.ma[l][m.l[l]];
	m.mc[e]=0;
	//compute new sets of basic and nonbasic variables.
	
	for(int j=0;j<=N;j++)
	if(m.r[j]==e) 
		m.r[j]=m.l[l];
	m.l[l]=e;
	return m;
}

MODULE initial_simplex(MODULE m1)
{
	float ori_c[N+1];
	for(int i=1;i<=N;i++)
		ori_c[i]=m1.mc[i];

	//MODULE m1=MM;
	
	int k=minL(m1.mb,M);
	if(m1.mb[k]>=0)
		return m1;
	for(int i=1;i<=M;i++){
		m1.ma[i][0]=-1;
		//m1.v=-x[0];
	}
	for(int i=1;i<=N;i++)//set the objective function to -x0.
		m1.mc[i]=0;
	m1.mc[0]=-1;
	m1.r[0]=0;
	int l=k;
	m1=PIVOT(m1,l,0);
	//simplex  line2----11
	float s=1000;
	int min_e=1000;
	for(int j=0;j<=N;j++){
	for(int i=0;i<=N;i++)
		if(m1.mc[m1.r[i]]>0){
			if(s>m1.mc[m1.r[i]]){
				s=m1.mc[m1.r[i]];min_e=m1.r[i];
		}
		}
		if(s==1000) break;
		s=1000;
		float deta_i=10000000;int l=1;
		for(int i=1;i<=M;i++){
			if(m1.ma[i][min_e]>0.000001){
				if(deta_i>m1.mb[i]/m1.ma[i][min_e]){
					deta_i=m1.mb[i]/m1.ma[i][min_e];
					l=i;}
			}
		}
			m1=PIVOT(m1,l,min_e);
	}

	
 if(!find(0,m1.l,M)){ //if the basic solutions sets /x0=0
	
	for(int i=0;i<=N;i++)//the final slack form with x0 removed and the original objective function restored
		m1.mc[i]=0;
	m1.v=0;
	
	for(int i=1;i<=N;i++){
		int k1=find(i,m1.l,M); //i在不在基本变量里面，若在找到对应的行标k1
		if(k1){
			for(int j=0;j<=N;j++)
				m1.mc[m1.r[j]]=m1.mc[m1.r[j]]-ori_c[i]*m1.ma[k1][m1.r[j]]; //计算非基本变量的系数
			m1.v=m1.v+m1.mb[k1]*ori_c[i];
		}
		else
			m1.mc[i]=m1.mc[i]+ori_c[i];
	}
	for(int i=1;i<=M;i++)//取消了x0的非基本变量
		m1.ma[i][0]=0;
			m1.mc[0]=0;
	return m1;
 }
 else{//若基本解不能使x0==0，退出

	 cout<<"infeasible"<<endl;exit(0);}

}

void simplex(MODULE m1)
{
	MODULE m2=initial_simplex(m1);//通过initial_simplex返回基本解可行的松弛型
	if(!flag){
	float s=1000;
	int min_e=1000;
	for(int j=0;j<=N;j++){
	for(int i=0;i<=N;i++)
		if(m2.mc[m2.r[i]]>0.000001){//每次找非基本变量前面系数最小的非基本变量换入
			if(s>m2.mc[m2.r[i]]){
				s=m2.mc[m2.r[i]];min_e=m2.r[i];
		}
		}
		if(s==1000) break;
		s=1000;
		float deta_i=10000000;int l=1;
		for(int i=1;i<=M;i++){//找使限制非基本变量最紧的行对应的基本变量作为换出变量
			if(m2.ma[i][min_e]>0){
				if(deta_i>m2.mb[i]/m2.ma[i][min_e]){
					deta_i=m2.mb[i]/m2.ma[i][min_e];
					l=i;}
			}
		}
			m2=PIVOT(m2,l,min_e);
	}
	
	cout<<"result="<<m2.v<<endl;
	}

	else
		cout<<"result ubound"<<endl;
}


int main()
{
	
	for(int i=1;i<=N;i++)
		cin>>MM.mc[i];
	for(int i=N+1;i<=N+M;i++)
		MM.mc[i]=0;

	for(int j=1;j<=M;j++)
		cin>>MM.mb[j];
	for(int i=1;i<=M;i++)
		MM.ma[i][0]=0;
	for(int i=1;i<=M;i++){
		for(int j=1;j<=N;j++)
			cin>>MM.ma[i][j];
		for(int j=N+1;j<=N+M;j++)
			MM.ma[i][j]=0;
	}
	simplex(MM);

	return 0;
}	

举例：z=2x+y

            x-y>=0;x+y-1<=0;y>=-1 

 转换为松弛型：令x=x1-x2 ; y=x3-x4  得到：

           z=2x1-2x2+x3-x4

         -x1+x2+x3-x4<=0;x1-x2+x3-x4<=1;-x3+x4<=1; x1,x2,x3,x4>=0 (N=4,M=3)

      

例2：当b[]数组中含有负值时 

z=2x1-x2

2x1-x2<=2;x1-5x2<=-4;x1,x2>=0 (N=2,M=2)