lpprodlg.cpp

运筹学单纯形算法

double CLPProDlg::LP_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)
{
	int * ptrBaseArray=new int[dim_m];
	for( int tmpcount=0; tmpcount < dim_m; tmpcount++ )	ptrBaseArray[ tmpcount ] = dim_n - dim_m + tmpcount;
	double * ptrTestArray=new double[dim_n];
	for( tmpcount=0; tmpcount < dim_n; tmpcount++ )	ptrTestArray[ tmpcount ] = ptr_C_Array[tmpcount];

	DWORD TimeTick=::GetTickCount();
	do{
		DWORD time_Last=::GetTickCount();
		if( time_Last - TimeTick > 10000 )
		{
			delete ptrTestArray;
			delete ptrBaseArray;
			value=3;
			return 0;
		}
		int  find_k = -1;
		for(int index=0;index<dim_n;index++)
		{
			if( ptrTestArray[index] > 0 ){
				if( find_k == -1 )	find_k = index;
				else if( ptrTestArray[index] > ptrTestArray[find_k] )	find_k=index;
			}
		}
		if( find_k == -1 ) 
		{
			int ZeroNum=0;
			for( int index=0; index<dim_n; index++)	if(ptrTestArray[index]==0) ZeroNum++;
			if( ZeroNum != dim_m ) value=2;
			break;
		}
		
		int find_l=-1;
		for(index=0;index<dim_m;index++)
		{
			//由于已经对 Pj<=0 进行了检验,所以必存在某一index
			//使得Matrix[index][find_k]>0
			if( ptr_A_Matrix[ index * dim_n + find_k ] >0 )
			{
				if( find_l == -1 )
					find_l =index;
				else 
				{
					double comp_index=ptr_b_Array[index]/ptr_A_Matrix[index*dim_n+find_k];
					double comp_find_l=ptr_b_Array[find_l]/ptr_A_Matrix[find_l*dim_n+find_k];
					if( comp_index < comp_find_l )	find_l=index;
				}
			}
		}
		if(find_l==-1)//存在Pj<=0;为无界解
		{
			delete ptrTestArray;
			delete ptrBaseArray;
			value=1;
			return 0;
		}
		//find_l,find_k都已经找到
		double * ptr_Matrix_tmp=new double[dim_m*dim_n];
		int icount=0,jcount=0;
		for(icount=0;icount<dim_m;icount++)
			for(jcount=0;jcount<dim_n;jcount++)
				ptr_Matrix_tmp[icount*dim_n+jcount]=ptr_A_Matrix[icount*dim_n+jcount];
		
		double find_bl=ptr_b_Array[find_l];
		double find_l_k=ptr_Matrix_tmp[find_l*dim_n+find_k];

		for(icount=0; icount<dim_m; icount++)
		{
			if( find_l==icount )
				ptr_b_Array[icount]/=find_l_k;
			else
			{
				double tmp=ptr_Matrix_tmp[icount*dim_n+find_k]*find_bl;
				tmp/=find_l_k;
				ptr_b_Array[icount]-=tmp;
			}

			for(int jcount=0; jcount<dim_n; jcount++)
			{
				if( find_l==icount )
					ptr_A_Matrix[icount*dim_n+jcount] /= find_l_k;
				else
				{
					double tmp = ptr_Matrix_tmp[icount*dim_n+find_k]*ptr_Matrix_tmp[find_l*dim_n+jcount];
					tmp /= find_l_k;
					ptr_A_Matrix[icount*dim_n+jcount]-=tmp;
				}
			}
		}
		double find_ck=ptrTestArray[find_k];
		for(index=0;index<dim_n;index++)
		{
			if(index==ptrBaseArray[find_l]) 
			{
				ptrTestArray[index]=find_ck/find_l_k;
				ptrTestArray[index]*=-1;
			}
			else
			{
				double tmp=ptr_Matrix_tmp[find_l*dim_n+index]/find_l_k;
				tmp*=find_ck;
				ptrTestArray[index]-=tmp;
			}
		}
		for(index=0;index<dim_m;index++)
		{
			if(index==find_l)
			{
				ptrBaseArray[index]=find_k;
				break;
			}
		}
		delete ptr_Matrix_tmp;
	}while(1);
	
	for(int icount=0;icount<dim_m;icount++)
			ptr_X_Array[ ptrBaseArray[icount] ] = ptr_b_Array[icount];

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

	delete ptrTestArray;
	delete ptrBaseArray;

	return WorkOut;
}

double CLPProDlg::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;	
}

double CLPProDlg::LD_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)
{
	int * ptrBaseArray=new int[dim_m];
	for( int tmpcount=0; tmpcount < dim_m; tmpcount++ )	ptrBaseArray[ tmpcount ] = dim_n - dim_m + tmpcount;

	double * ptrTestArray=new double[dim_n];
	for( tmpcount=0; tmpcount < dim_n; tmpcount++ )	ptrTestArray[ tmpcount ] = ptr_C_Array[tmpcount];

	DWORD TimeTick=::GetTickCount();
	do{
		DWORD time_Last=::GetTickCount();
		if( time_Last - TimeTick > 10000 )
		{
			delete ptrTestArray;
			delete ptrBaseArray;
			value=3;
			return 0;
		}
		int  find_k = -1;//行数
		for(int index=0;index<dim_m;index++)
		{
			if( ptr_b_Array[index] < 0 ){
				if( find_k == -1 )	find_k = index;
				else if( ptr_b_Array[index] < ptr_b_Array[find_k] ) find_k=index;
			}
		}
		if( find_k == -1 ) 
		{
			int ZeroNum=0;
			for( int index=0; index<dim_n; index++)	if(ptrTestArray[index]==0) ZeroNum++;
			if( ZeroNum != dim_m ) value=2;
			break;
		}

		int find_l=-1;//列数
		for(index=0;index<dim_n;index++)
		{
			if( ptr_A_Matrix[ find_k * dim_n + index ] < 0 )
			{
				if( find_l == -1 )
					find_l =index;
				else 
				{
					double comp_index=ptrTestArray[index]/ptr_A_Matrix[find_k*dim_n+index];
					double comp_find_l=ptrTestArray[find_l]/ptr_A_Matrix[find_k*dim_n+find_l];
					if( comp_index < comp_find_l )	find_l=index;
				}
			}
		}
		if(find_l==-1)//存在Pj<=0;无可行解
		{
			delete ptrTestArray;
			delete ptrBaseArray;
			value=1;
			return 0;
		}

		double * ptr_Matrix_tmp=new double[dim_m*dim_n];
		int icount=0,jcount=0;
		for(icount=0;icount<dim_m;icount++)
			for(jcount=0;jcount<dim_n;jcount++)
				ptr_Matrix_tmp[icount*dim_n+jcount]=ptr_A_Matrix[icount*dim_n+jcount];
		
		double find_bk=ptr_b_Array[find_k];
		double find_k_l=ptr_Matrix_tmp[find_k*dim_n+find_l];

		for(icount=0; icount<dim_m; icount++)
		{
			if( find_k==icount )
				ptr_b_Array[icount]/=find_k_l;
			else
			{
				double tmp=ptr_Matrix_tmp[icount*dim_n+find_l]*find_bk;
				tmp/=find_k_l;
				ptr_b_Array[icount]-=tmp;
			}

			for(int jcount=0; jcount<dim_n; jcount++)
			{
				if( find_k==icount )
					ptr_A_Matrix[icount*dim_n+jcount] /= find_k_l;
				else
				{
					double tmp = ptr_Matrix_tmp[icount*dim_n+find_l]*ptr_Matrix_tmp[find_k*dim_n+jcount];
					tmp /= find_k_l;
					ptr_A_Matrix[icount*dim_n+jcount]-=tmp;
				}
			}
		}
		double find_ck=ptrTestArray[find_l];
		for(index=0;index<dim_n;index++)
		{
			if(index==ptrBaseArray[find_k]) 
			{
				ptrTestArray[index]=find_ck/find_k_l;
				ptrTestArray[index]*=-1;
			}
			else
			{
				double tmp=ptr_Matrix_tmp[find_k*dim_n+index]/find_k_l;
				tmp*=find_ck;
				ptrTestArray[index]-=tmp;
			}
		}
		for(index=0;index<dim_m;index++)
		{
			if(index==find_k)
			{
				ptrBaseArray[index]=find_l;
				break;
			}
		}
		delete ptr_Matrix_tmp;
		
	}while(1);

	for(int icount=0;icount<dim_m;icount++)
			ptr_X_Array[ ptrBaseArray[icount] ] = ptr_b_Array[icount];

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

	delete ptrTestArray;
	delete ptrBaseArray;

	return WorkOut;
}