正文

单纯形法(最终版本)2005-09-24 16:05:00

【评论】 【打印】 【字体: 】 本文链接:http://blog.pfan.cn/elva6401/5212.html

分享到:

#include<iomanip.h>

#include <fstream.h>
ifstream fin("a.txt");
ofstream fout("b.txt");
#define Min -99999
#define Max 99999
int main()
{
float a[100][100];//方程的系数
float check[100];//检验数
float Csum;//Cj-Zj的和
float z[100];//目标函数
int m,n;//分别为变量数和方程个数
int XB[100];//基变量
int mark[100];//标记某个基是否是基变量
int in,out;//换入和换出基的位置
float inN,outN;//换入和换出基的值
float result;//最优值
int i,j;
int sum;
int posx,posy;
float temp,temp2;
fin>>m>>n;
for(i=0;i<m;i++) fin>>z[i];
for(i=0;i<n;i++)
for(j=0;j<=m;j++)
fin>>a[i][j];
//下面是找基的过程
for(i=0;i<100;i++)
XB[i]=mark[i]=0;
//先来找容易看出来的基
for(i=0;i<m;i++)
{
posy=0;
sum=0;
for(j=0;j<n;j++)
if(a[j][i]!=0)
{
if(sum==0) posy=j;
sum++;
}
if(sum==1 && mark[i]==0) //第i个变量只有第j个方程有系数且不是基
{
temp=a[posy][i];
for(j=0;j<=m;j++)//把它化为基
a[posy][j]/=temp;
XB[posy]=i;
mark[i]=1;
}
}

posy=0;
while(posy<n)
{
if(XB[posy]==0)//表示第posy个方程还没有找出基
{
posx=0;
while(mark[posx]==1 || a[posy][posx]==0)//对应的变量是基或等于0,换下一个变量
posx++;
temp=a[posy][posx];
for(i=0;i<=m;i++)//化为1
a[posy][i]/=temp;
for(i=0;i<n;i++)//其它的行化为0
{
 if(i!=posy)
 {
 temp2=a[i][posx];
 for(j=0;j<=m;j++)
 a[i][j]=a[i][j]-a[posy][j]*temp2;
 }
}
XB[posy]=posx;
mark[posx]=1;
}
posy++;
}

//把基全找出来了,放到了XB数组中,后面进行换基迭代

//算检验数
for(i=0;i<m;i++)
{
Csum=0;
for(j=0;j<n;j++)
Csum+=z[XB[j]]*a[j][i];
check[i]=z[i]-Csum;
}
//列出第1个表
/*fout<<"CB XB b ";
for(i=0;i<m;i++)
fout<<"X"<<i+1<<" ";
fout<<endl;*/
for(i=0;i<n;i++)
{
fout<<z[XB[i]]<<" ";
fout<<"X"<<XB[i]+1<<" ";
fout.setf(ios::fixed);
fout<<setprecision(3)<<a[i][m]<<" ";
//fout<<a[i][m]<<"  ";
for(j=0;j<m;j++)
fout<<a[i][j]<<"  ";
fout<<endl;
}
fout<<"Cj-Zj ";
for(i=0;i<m;i++)
fout<<check[i]<<" ";
fout<<endl;

 


//看一下b中是否有负数,要用对偶单纯形法
while(1)
{
 out=-1;
 outN=Max;
 for(i=0;i<n;i++)
 if(a[i][m]<outN)
 {
  out=i;
  outN=a[i][m];
 }
 if(outN>0) break;//b全部大于0
 in=-1;
 inN=Max;
 for(i=0;i<m;i++)
 if(a[out][i]<0 && check[i]<0)
 {
  if(check[i]/a[out][i]<inN)
  {
   in=i;
   inN=check[i]/a[out][i];
  }
 }
 if(inN==Max) {fout<<"\n无界解\n";break;}
temp=a[out][in];
 for(i=0;i<=m;i++)
  a[out][i]=a[out][i]/temp;
 for(i=0;i<n ;i++)
 if(i!=out)
 {
 temp=a[i][in];
 for(j=0;j<=m;j++)
  a[i][j]=a[i][j]-a[out][j]*temp;
 }
 temp=check[in];
 for(i=0;i<m;i++)
  check[i]=check[i]-a[out][i]*temp;
fout<<"入X"<<in+1<<" 出X"<<XB[out]+1<<endl;
for(i=0;i<n;i++)
{
fout<<z[XB[i]]<<" ";
fout<<"X"<<XB[i]+1<<" ";
fout<<a[i][m]<<"  ";
for(j=0;j<m;j++)
fout<<a[i][j]<<"  ";
fout<<endl;
}
fout<<"  Cj-Zj  ";
for(i=0;i<m;i++)
fout<<check[i]<<" ";
fout<<endl;
XB[out]=in;
}

 


while(1)
{
//经验是否达到最优
in=-1;
inN=Min;
for(i=0;i<m;i++)
if(check[i]>inN)//找最大的作为换入基
{
in=i;
inN=check[i];
}
if(inN<=0) break;//检验全比0小,达到最优.退出
//找换出基
out=-1;
outN=Max;
for(i=0;i<n;i++)
{
 if(a[i][in]>0)
 {
  if(a[i][m]/a[i][in]<outN)
  {
  outN=a[i][m]/a[i][in];
  out=i;
  }
 }
}
if(outN==Max) {fout<<"\n无界解\n";break;}//没有换出基,无解.退出
temp=a[out][in];
 for(i=0;i<=m;i++)
  a[out][i]=a[out][i]/temp;
 for(i=0;i<n ;i++)
 if(i!=out)
 {
 temp=a[i][in];
 for(j=0;j<=m;j++)
  a[i][j]=a[i][j]-a[out][j]*temp;
 }
 temp=check[in];
 for(i=0;i<m;i++)
  check[i]=check[i]-a[out][i]*temp;
fout<<"入X"<<in+1<<" 出X"<<XB[out]+1<<endl;
for(i=0;i<n;i++)
{
fout<<z[XB[i]]<<" ";
fout<<"X"<<XB[i]+1<<" ";
fout<<a[i][m]<<"  ";
for(j=0;j<m;j++)
fout<<a[i][j]<<"  ";
fout<<endl;
}
fout<<"  Cj-Zj  ";
for(i=0;i<m;i++)
fout<<check[i]<<" ";
fout<<endl;
XB[out]=in;
}
result=0;
for(i=0;i<n;i++)
{
//fout<<a[i][m]<<" "<<z[XB[i]]<<" ";
result+=a[i][m]*z[XB[i]];
}
float c[100];
for(i=0;i<m;i++)
c[i]=0;
for(i=0;i<n;i++)
c[XB[i]]=a[i][m];
for(i=0;i<m;i++)
fout<<"X"<<i+1<<"="<<c[i]<<" ";
fout<<"\n result="<<result<<endl;
}

阅读(5999) | 评论(6)


版权声明:编程爱好者网站为此博客服务提供商,如本文牵涉到版权问题,编程爱好者网站不承担相关责任,如有版权问题请直接与本文作者联系解决。谢谢!

评论

loading...
您需要登录后才能评论,请 登录 或者 注册