simplex.cpp

vc下实现的最优化程序代码

// Simplex.cpp: implementation of the CSimplex class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "REOPT.h"
#include "Simplex.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CSimplex::CSimplex()
{
	m_IsLimitary=true;
}

CSimplex::~CSimplex()
{

}

void CSimplex::Compute_B_N1_N2()
{
	int i,j;
	for(i=0;i<m_QuatCunt;i++)
	{
		for(j=0;j<m_VarB;j++)
		{
			m_B(i,j)=m_A(i,m_NumB[j]);
		}
		for(j=0;j<m_VarL;j++)
		{
			m_N1(i,j)=m_A(i,m_NumN1[j]);
		}
		for(j=0;j<m_VarU;j++)
		{
			m_N1(i,j)=m_A(i,m_NumN2[j]);
		}
	}
}

void CSimplex::Compute_Cb_Cn1_Cn2()
{
	int i,j;
	for(i=0;i<m_TarCunt;i++)
	{
		for(j=0;j<m_VarB;j++)
		{
			m_Cb(i,j)=m_C(i,m_NumB[j]);
		}
		for(j=0;j<m_VarL;j++)
		{
			m_Cn1(i,j)=m_C(i,m_NumN1[j]);
		}
		for(j=0;j<m_VarU;j++)
		{
			m_Cn2(i,j)=m_C(i,m_NumN2[j]);
		}
	}
}

float CSimplex::ComputeQ(LorU J)
{
	assert(J==L||J==U);

	CArrayMatrix ta,tb;
	float Q,Q1,Q2,min;
	Q=Q1=Q2=9999999;
	min=-1;

	ta=~m_B*m_N1;
	Q=m_U(m_k,0)-m_L(m_k,0);
	min=Q;
	if(J==L)
	{
		if(ta(m_k,0)>0)
		{
			tb=~m_B*m_b-m_L;
			Q1=tb.FindExtremumInVector(m_r,false)/ta(m_k,0);
		}
		else if(ta(m_k,0)<0)
		{
			tb=m_U-~m_B*m_b;
			Q2=tb.FindExtremumInVector(m_r,false)/ta(m_k,0)*(-1);
		}
		else
		{
			m_IsLimitary=false;
		}
	}
	else if(J==U)
	{
		if(ta(m_k,0)<0)
		{
			tb=~m_B*m_b-m_L;
			Q1=tb.FindExtremumInVector(m_r,false)/ta(m_k,0)*(-1);
		}
		else if(ta(m_k,0)>0)
		{
			tb=m_U-~m_B*m_b;
			Q2=tb.FindExtremumInVector(m_r,false)/ta(m_k,0);
		}
		else
		{
			m_IsLimitary=false;
		}
	}
	if(m_IsLimitary)
	{
		min=Q1<min?Q1:min;
		min=Q2<min?Q2:min;
	}
	return min;
}

void CSimplex::SetVariableCount(int m,int l,int n,int p,int q,int u)
{
    m_VarB=m;//基变量数
    m_VarL=l;//L集合变量数
	m_VarU=n;//N集合变量数
	m_VarCunt=m+l+n;
	m_QuatCunt=p;//约束方程数
	m_TarCunt=q;
	m_UnQuatCunt=u;
}

float CSimplex::Compute_lmd()
{
	CArrayMatrix temp1,temp2;
	int k1,k2;
	float lmd1,lmd2,lmd;
	
	k1=k2=-1;
	lmd1=lmd2=-0.1;

	temp1.SetMatrixSize(m_VarL);
	temp1=m_Cb*(~m_B)*m_N1-m_Cn1;
	lmd1=temp1.FindExtremumInVector(k1,true);
	temp2.SetMatrixSize(m_VarU);
	temp2=(m_Cb*(~m_B)*m_N2-m_Cn2)*(-1);
	lmd2=temp2.FindExtremumInVector(k2,true);

	if(lmd1>=lmd2)
	{
		m_k=k1;
		lmd=lmd1;
	}
	else
	{
		m_k=k2;
		lmd=lmd2;
	}
	return lmd;
}

void CSimplex::InitData(int m,int l,int n,int p,int q,int u)
{
    SetVariableCount(m,l,n,p,q,u);
	SetScript();
	SetMatrixSize();
	Compute_B_N1_N2();
	Compute_Cb_Cn1_Cn2();
	m_f=Compute_f0();
	m_lmd=-1;
	m_k=m_r=-1;
	m_Kind=bothnot;

	for(int i=0;i<m_VarCunt;i++)
		m_X(i,0)=0;
	SetMatrix_Xb();
	SetMatrix_Big();
}

void CSimplex::SetMatrixSize()
{
	m_A.SetMatrixSize(m_VarCunt*m_QuatCunt);
	m_B.SetMatrixSize(m_VarB*m_QuatCunt);
	m_N1.SetMatrixSize(m_VarL*m_QuatCunt);
	m_N2.SetMatrixSize(m_VarU*m_QuatCunt);

	m_C.SetMatrixSize(m_VarCunt*m_TarCunt);
	m_Cb.SetMatrixSize(m_VarB*m_TarCunt);
	m_Cn1.SetMatrixSize(m_VarL*m_TarCunt);
	m_Cn2.SetMatrixSize(m_VarU*m_TarCunt);

	m_X.SetMatrixSize(m_VarCunt);
	m_Xb.SetMatrixSize(m_VarB);
	m_Xn1.SetMatrixSize(m_VarL);
	m_Xn2.SetMatrixSize(m_VarU);

	m_b.SetMatrixSize(m_QuatCunt);

	m_L.SetMatrixSize(m_UnQuatCunt);
	m_U.SetMatrixSize(m_UnQuatCunt);
}

void CSimplex::SetScript()
{
	m_NumB.SetSize(m_VarB);
	m_NumN1.SetSize(m_VarL);
	m_NumN2.SetSize(m_VarU);

	int i=0;
	for(i=0;i<m_VarB;i++)
		m_NumB[i]=m_VarCunt-i;
	for(i=0;i<m_VarL;i++)
		m_NumN1[i]=m_VarCunt-i-m_VarB;
	for(i=0;i<m_VarU;i++)
		m_NumN2[i]=m_VarCunt-i-m_VarB-m_VarL;
}



void CSimplex::JudgeScript(int k)
{
	int i=0;
	for (i=0;i<m_VarL;i++)
	{
		if(k==m_NumN1[i])
		{
			m_Kind=L;
			break;
		}
	}
	
	for (i=0;i<m_VarU;i++)
	{
		if(k==m_NumN2[i])
		{
			m_Kind=U;
			break;
		}
	}
}

void CSimplex::ExecuteSLP()
{
	float epson=0.0000001;
//	InitData();
	
	while(m_lmd=Compute_lmd()>=epson)//得到m_k,lmd>0
	{
		JudgeScript(m_k);//确定m_k属于哪个集合,
		m_Q=ComputeQ(m_Kind);//得到m_r
		if(m_IsLimitary)
		{
			ReviseData();
		}
		else
		{
			AfxMessageBox("无有限最优解！");
			goto end;
		}
	}

	PrintResult();//当lmd<0时，得到最优解，打印结果
end:
	;
}

void CSimplex::ReviseData()
{
	assert(m_IsLimitary);
	if(m_Q==m_U(m_k,0)-m_L(m_k,0))
	{
		if(m_Kind==L)
			m_X(m_k,0)=m_U(m_k,0);
		else
			m_X(m_k,0)=m_L(m_k,0);
	}
	else
	{
		assert(ReSetBase());//换基；
		GaussMethod();//左端高斯消去
		ReviseRight();//修正右端
	}

	ReviseMatrix_X();
	ReviseSetScript();
	Compute_B_N1_N2();
	Compute_Cb_Cn1_Cn2();
	SetMatrix_Xb();
}

void CSimplex::GaussMethod()
{
	int i,j;

	for(i=0;i<m_QuatCunt+1;i++)
	{
		if(i==m_r)
			continue;
		else
		{
			for(j=0;j<m_VarCunt;i++)
			{
				m_Big(i,j)=-m_Big(m_r,m_k)*m_Big(m_r,j)+m_Big(i,j);
			}
		}
	}
}

void CSimplex::ReviseRight()
{
	if(m_Kind==L)
	{
		m_X(m_k,0)=m_L(m_k,0)+m_Q;
		CArrayMatrix Ak,temp1,temp2;
		Ak.SetMatrixSize(m_QuatCunt);

		temp1=~m_B*m_N1;
		for(int i=0;i<m_QuatCunt;i++)
		{
			Ak(i,0)=temp1(i,m_k);
		}

		m_Xb=~m_B*m_b-Ak*m_Q;
		temp2=(m_Cb*(~m_B)*m_N1-m_Cn1);
		m_f=m_f-m_Q*temp2(m_k,0);
	}
	else
	{
		m_X(m_k,0)=m_U(m_k,0)-m_Q;
		CArrayMatrix Ak,temp1,temp2;
		Ak.SetMatrixSize(m_QuatCunt);

		temp1=~m_B*m_N1;
		for(int i=0;i<m_QuatCunt;i++)
		{
			Ak(i,0)=temp1(i,m_k);
		}

		m_Xb=~m_B*m_b+Ak*m_Q;
		temp2=(m_Cb*(~m_B)*m_N2-m_Cn2);
		m_f=m_f+m_Q*temp2(m_k,0);
	}
}
void CSimplex::PrintResult()
{
	
}

BOOL CSimplex::ReSetBase()
{
	bool IsSuccess=false;
	int i=0;

	for(i=0;i<m_VarB;i++)
	{
		if(m_NumB[i]==m_r)
		{
			m_NumB[i]=m_k;
			IsSuccess=true;
			break;
		}
	}

	return IsSuccess;
}

BOOL CSimplex::ReviseSetScript()
{
	BOOL IsSuccess=false;
	int i,j;

	if(m_Kind==L)
	{
		for(i=0;i<m_VarL;i++)
		{
			if(m_NumN1[i]==m_k)
			{
				for(j=i;j<m_VarL-1;j++)
				{
					m_NumN1[j]=m_NumN1[j+1];
				}
				m_VarL=m_VarL-1;
				m_VarU=m_VarU+1;
				m_NumN2[m_VarU-1]=m_k;
				m_NumN1.SetSize(m_VarL);
				m_NumN2.SetSize(m_VarU);
				IsSuccess=true;
				break;
			}
		}
	}
	else
	{
		for(i=0;i<m_VarU;i++)
		{
			if(m_NumN2[i]==m_k)
			{
				for(j=i;j<m_VarL-1;j++)
				{
					m_NumN2[j]=m_NumN2[j+1];
				}
				m_VarL=m_VarL+1;
				m_VarU=m_VarU-1;
				m_NumN1[m_VarL-1]=m_k;
				m_NumN1.SetSize(m_VarL);
				m_NumN2.SetSize(m_VarU);
				IsSuccess=true;
				break;
			}
		}
	}

	return IsSuccess;
}

float CSimplex::Compute_f0()
{
	float f=0.0;
	CArrayMatrix temp;
	temp=m_Cb*~m_B*m_b-(m_Cb*~m_B*m_N1-m_Cn1)*m_Xn1-(m_Cb*~m_B*m_N2-m_Cn2)*m_Xn2;
	f=temp(0,0);
	return f;
}

void CSimplex::SetMatrix_Xb()
{
	int i=0;
	for(i=0;i<m_VarB;i++)
	{
		m_Xb(i,0)=m_X(m_NumB[i],0);
	}
}


void CSimplex::ReviseMatrix_X()
{
	int i=0;
	for(;i<m_VarB;i++)
	{
		m_X(m_NumB[i],0)=m_Xb(i,0);
	}
}

void CSimplex::SetMatrix_Big()
{
	m_Big.SetMatrixSize((m_QuatCunt+m_TarCunt)*(m_VarCunt));
	int i,j;

	for(i=0;i<m_QuatCunt;i++)
	{
		for(j=0;j<m_VarCunt;j++)
			m_Big(i,j)=m_A(i,j);
	}

	for(i=m_QuatCunt;i<m_QuatCunt+m_TarCunt;i++)
	{
		for(j=0;j<m_VarCunt;j++)
			m_Big(i,j)=m_C(i-m_QuatCunt,j);
	}
}