改进单纯形算法.txt

运筹学单纯形算法

//模型:
//	max   f＝CX
//	s.t.  AX=b (X>=0, b>=0)
//dim_m：矩阵A的行数
//dim_n：矩阵A的列数
//ptr_C_Array：C向量数组，数组长度为dim_n;
//ptr_b_Array：b向量数组，长度为dim_m;
//ptr_A_Matrix：矩阵A数组，长度为[dim_m*dim_n];
//ptr_X_Array：解向量输出数组，长度为dim_n;
//value：返回值类型，0为唯一解，1为无界解，2为无穷多解。
//函数返回值为f

/*例子
	单纯形法m*n (m= 3, n=5 )
	目标函数C＝[2,1,0,0,0]
	向量b＝[15,24,5]
	矩阵A=
	[ 0,5,1,0,0
	  6,2,0,1,0
	  1,1,0,0,1]
*/

double LG_WorkOut(double * ptr_C_Array,double * ptr_b_Array,
				double  * ptr_A_Matrix,int dim_m,int dim_n,double * & ptr_X_Array,int & value)
{
	//G矩阵（（m＋1）×（m＋1））初始化
		double * ptr_Matrix=new double[(dim_m+1)*(dim_m+1)];
	int icount=0,jcount=0,index=0;
	for(icount=0;icount<dim_m;icount++)
		for(jcount=0;jcount<dim_m;jcount++)
			if(icount==jcount)	ptr_Matrix[icount*(dim_m+1)+jcount]=1;
			else ptr_Matrix[icount*(dim_m+1)+jcount]=0;
	for(icount=0;icount<dim_m+1;icount++)	ptr_Matrix[dim_m*(dim_m+1)+icount]=0;
	for(icount=0;icount<dim_m;icount++)	ptr_Matrix[icount*(dim_m+1)+dim_m]=ptr_b_Array[icount];

	//初始化记录rj的数组
	CArray<Value_R *,Value_R *> V_Value_R;
	Value_R * ptrValue_R=new Value_R[dim_n-dim_m];
	for(icount=0;icount<dim_n-dim_m;icount++)
	{
		ptrValue_R[icount].index=icount;
		ptrValue_R[icount].value=ptr_C_Array[icount];
		V_Value_R.Add(ptrValue_R+icount);
	}

	//初始化Base_Array数组
	int * ptrBaseArray=new int[dim_m];
	for( icount=0; icount < dim_m; icount++ )
		ptrBaseArray[ icount ] = dim_n - dim_m + icount;

	DWORD TimeTick=::GetTickCount();
	do{
		DWORD time_Last=::GetTickCount();
		if( time_Last - TimeTick > 10000 )
		{
			delete ptr_Matrix;
			delete ptrBaseArray;
			delete ptrValue_R;
			value=3;
			return 0;
		}
		int find_l=-1;//当前在A矩阵中要进行转轴的列下标
		int find_index=-1;
		for(icount=0; icount < dim_n-dim_m; icount++)
		{
			if( V_Value_R.GetAt(icount)->value > 0 )
			{
				if( find_l == -1 )
				{
					find_l = V_Value_R.GetAt(icount)->index;
					find_index = icount;
				}
				else if( V_Value_R.GetAt(icount)->value > V_Value_R.GetAt(find_index)->value )
				{
					find_l = V_Value_R.GetAt(icount)->index;
					find_index = icount;
				}
			}
		}
		if( find_l == -1 )//没有找到 > 0 的检验数
		{
			int ZeroNum=0;
			for( icount=0; icount < dim_n-dim_m; icount++) 
				if( V_Value_R.GetAt( icount )->value == 0 ) ZeroNum++;
			if( ZeroNum ) value = 2;
			break;
		}
		
		int find_k=-1; //当前在A矩阵中要进行转轴的行下标
		double * ptr_yt=new double[dim_m+1];  //构造yt数组，判断其是否均<＝0
		for(icount=0;icount<dim_m;icount++)
		{
			double yt=0;
			for(jcount=0;jcount<dim_m;jcount++)
				yt+=ptr_Matrix[icount*(dim_m+1)+jcount]*ptr_A_Matrix[jcount*dim_n+find_l];
			ptr_yt[icount]=yt;
		}
		ptr_yt[dim_m]=V_Value_R.GetAt(find_index)->value;

		for(icount=0;icount<dim_m;icount++)
		{
			if(ptr_yt[icount]>0) 
			{
				if(find_k == -1)	find_k = icount;
				else
				{
					double comp_index=ptr_Matrix[icount*(dim_m+1)+dim_m]/ptr_yt[icount];
					double comp_find_k=ptr_Matrix[find_k*(dim_m+1)+dim_m]/ptr_yt[find_k];
					if( comp_index < comp_find_k )	find_k = icount;
				}
			}
		}
		if(find_k==-1)//当前规划在可行域内无下界。
		{
			delete ptr_Matrix;
			delete ptrBaseArray;
			delete ptrValue_R;
			delete ptr_yt;
			value=1;
			return 0;
		}
		
		//ptrBaseArray
		int changeHang=ptrBaseArray[find_k];//得到当前出基列的列号
		ptrBaseArray[find_k]=find_l;//添加入基列的列号到出基列的位置
		V_Value_R.GetAt(find_index)->index=changeHang;//将出基列的列号置于原来在检验数数组中入基列列号的位置

		//find_l,find_k分别是矩阵A的列下标和行下标
		//重新购置G矩阵（(m＋1)×(m＋1)）
		double * ptrTmpMatrix=new double[(dim_m+1)*(dim_m+1)];
		for(icount=0;icount<dim_m+1;icount++)
			for(jcount=0;jcount<dim_m+1;jcount++)
				ptrTmpMatrix[icount*(dim_m+1)+jcount]=ptr_Matrix[icount*(dim_m+1)+jcount];

		for(icount=0;icount<dim_m+1;icount++)
			for(jcount=0;jcount<dim_m+1;jcount++)
				if(icount==find_k)	
					ptr_Matrix[icount*(dim_m+1)+jcount] = ptrTmpMatrix[find_k*(dim_m+1)+jcount]/ptr_yt[find_k];
				else ptr_Matrix[icount*(dim_m+1)+jcount] -= ptrTmpMatrix[find_k*(dim_m+1)+jcount]*ptr_yt[icount]/ptr_yt[find_k];
		
		delete ptrTmpMatrix;
		delete ptr_yt;
	
		//更新V_Value_R数组
		for(icount=0;icount<dim_n-dim_m;icount++)
		{
			double value_and=0;
			int ValueIndex=V_Value_R.GetAt(icount)->index;
			for(jcount=0;jcount<dim_m;jcount++)
				value_and+=ptr_Matrix[dim_m*(dim_m+1)+jcount]*ptr_A_Matrix[jcount*dim_n+ValueIndex];
			V_Value_R.GetAt(icount)->value=ptr_C_Array[ValueIndex]+value_and;
		}
	}while(1);

	for(icount=0;icount<dim_m;icount++)
			ptr_X_Array[ ptrBaseArray[icount] ] = ptr_Matrix[icount*(dim_m+1)+dim_m];

	double WorkOut=0;
	for(icount=0;icount<dim_n;icount++)
		WorkOut+=ptr_C_Array[icount]*ptr_X_Array[icount];

	delete ptr_Matrix;
	delete ptrBaseArray;
	delete ptrValue_R;

	return WorkOut;	
}