200622130049.cpp

单纯形算法描述 在地学分析课上学到的一点东西 用ＶＣ写的

#include <stdio.h>
#include <stdlib.h>
#define MAX 100
#define STP 100
int stop=1; //迭代记数变量
int status; //iterative迭代返回值：1唯一最优，0无界解，-1无穷多最优解 -2迭代超过限制次数
int step=1; //目前阶段
double a[MAX][MAX],b[MAX],c[MAX],temp_c[MAX],max=0; //方程组相关系数
int num_x; //变量个数 
int num_st; //约束方程数
int num_ar=0; //人工变量个数
int arti[MAX]; //人工变量下标
int base[MAX]; //基变量下标
int ma_mi; //1为求最大值，2为求最小值
void create(); //建立方程组
void iterative(); //单纯型法迭代
void output(); //输出结果
void banner(); //打印程序标题
void exchange(double c[MAX],double temp_c[MAX]); //交换两阶段价值系数
void show(); //输出方程组
void main() {
int i,j,k;
banner();
create();
//保存原价值系数，转换为第一阶段价值系数
for(i=1;i<=num_x;i++) {
  k=0;
  for(j=1;j<=num_ar;j++) if(i==arti[j]) k=1;
  if(k==1) temp_c[i]=-1;
  else temp_c[i]=0;
}
exchange(c,temp_c);
printf("\n\n第一阶段问题为：\n\n");
show();
step++;
printf("\n\n按回车开始第一阶段迭代");
getchar(); 
getchar();
iterative();
if(status==-2) {
  puts("迭代超过限制次数强行终止！\n");
  puts("\n按回车结束");
  getchar();
  exit(0);
}
output();
if(max!=0) {
  puts("\n\n原问题无可行解。\n");
  puts("\n按回车结束");
  getchar();
  exit(0);
}
//转换为第二阶段价值系数
exchange(c,temp_c);
//把人工变量列全设为0
for(i=1;i<=num_ar;i++) {
  c[arti[i]]=0;
  for(j=1;j<=num_st;j++) a[j][arti[i]]=0;
}
puts("\n\n第二阶段问题为：\n\n");
show();
puts("\n\n按回车开始第二阶段迭代");
getchar(); 
iterative();
switch(status) {
case 1:
  output();
  puts("\n\n原问题有唯一最优解。\n");
  puts("\n按回车结束");
  getchar();
  exit(0);
case 0:
  puts("\n\n原问题为无界解。\n");
  puts("\n按回车结束");
  getchar();
  exit(0);
case -1:
  output();
  puts("\n\n原问题有无穷多最优解。\n");
  puts("\n按回车结束");
  getchar();
  exit(0);
case -2:
  puts("迭代超过限制次数强行终止！\n");
  puts("\n按回车结束");
  getchar();
  exit(0);
}//switch
}
void banner() {

printf("\t\t         单纯型法解线性规划问题\n");

printf("\n");
}
void show() {
//对方程组以自然的格式输出，系数为零的x不显示
//为1的不显示系数1，-1系数只显示负号
int i,j,k;
switch(step) {
case 1:
  printf("min z= ");
  printf("x[%d]",arti[1]);
  for(i=2;i<=num_ar;i++) printf(" + x[%d]",arti[i]);
  break;
case 2:
printf("max z= ");
printf("%lg x[%d]",c[1],1);
for(i=2;i<=num_x;i++) {
  if(c[i]==1) printf(" + x[%d]",i);
  else if(c[i]==-1) printf(" - x[%d]",i);
  else if(c[i]>=0) printf(" +%lg x[%d]",c[i],i);
  else printf(" %lg x[%d]",c[i],i);
}
break;
}
printf("\nst:\n");
for(i=1;i<=num_st;i++) {
  k=0;
  for(j=1;j<=num_x;j++) { 
   if(a[i][j]!=0) {
    if(a[i][j]==1&&k!=0) printf(" + x[%d]",j);
    else if(a[i][j]==1&&k==0) printf("  x[%d]",j);
    else if(a[i][j]==-1) printf(" - x[%d]",j);
    else if(a[i][j]>=0&&k!=0) printf(" +%lg x[%d]",a[i][j],j);
    else if(a[i][j]>=0&&k==0) printf(" %lg x[%d]",a[i][j],j);
    else printf(" %lg x[%d]",a[i][j],j);
    k=1;
   }
  }
  printf(" == %lg\n",b[i]); 
}
printf(" x[1]~x[%d]>=0",num_x);
}
void exchange(double c[MAX],double temp_c[MAX]) 
{
int i;
double temp[MAX];
for(i=1;i<=num_x;i++) {
  temp[i]=temp_c[i];
  temp_c[i]=c[i];
  c[i]=temp[i];
}
}
void create() {
//输入方程组系数，每个方程输完后回显确认
int i,j,k,re_st[MAX],tnum_x,num_addv=0,num_ba=0;
char confirm;

while(1) {
  printf("请选择：1、max，2、min：(1/2)");
  scanf("%d",&ma_mi);
  if(ma_mi!=1&&ma_mi!=2) printf("输入错误，重新选择。");
  else break;
}

while(1) {
  printf("指定变量个数：");
  scanf("%d",&num_x);
  printf("输入价值系数c1-c%d:\n",num_x);
  for(i=1;i<=num_x;i++) {   
   printf("c%d=",i);
   scanf("%lf",&c[i]);
  }
  if(ma_mi==1) printf("max z= ");
  else printf("min z= ");
  printf("%lg x[%d]",c[1],1);
  for(i=2;i<=num_x;i++) {
   if(c[i]>=0) printf(" +%lg x[%d]",c[i],i);
   else printf("  %lg x[%d]",c[i],i);
  }
printf("\n正确吗？：(y/n)");
getchar();
confirm=getchar();
if (confirm=='y') break;
else if(confirm=='n') continue;
}
printf("输入约束方程组个数：");
scanf("%d",&num_st);
for(i=1;i<=num_st;i++) {
  printf("st.%d:\n",i);
  while(1) {
   printf("请选择：1、==，2、>=，3、<= ：(1/2/3)");
   scanf("%d",&re_st[i]);
   if(re_st[i]!=1&&re_st[i]!=2&&re_st[i]!=3) printf("输入错误，请重新选择。");
   else break;
  }
  printf("输入系数：\n");
  for(j=1;j<=num_x;j++) {  
   printf("a%d=",j);
   scanf("%lf",&a[i][j]);
  }
  printf("输入约束条件右端项：\nb%d=",i);
  scanf("%lf",&b[i]);
  
  printf("st.%i:\n",i);
  printf("%lg x[%d]",a[i][1],1);
  for(j=2;j<=num_x;j++) {
   if(a[i][j]>=0) printf(" +%lg x[%d]",a[i][j],j);
   else printf(" %lg x[%d]",a[i][j],j);
  }
  switch(re_st[i]) {
   case 1: printf(" == %lg",b[i]); break;
   case 2: printf(" >= %lg",b[i]); break;
   case 3: printf(" <= %lg",b[i]); break;
  }
  while(1) {
   printf("\n正确吗？(y/n)");
   getchar();
   confirm=getchar();
   if (confirm=='y') break;
   else if(confirm=='n') {i-=1; break;}
  }
}
//显示输入的方程组
printf("\n原问题为：\n\n");
if(ma_mi==1) printf("max z= ");
else printf("min z= ");
printf("%lg x[%d]",c[1],1);
for(i=2;i<=num_x;i++) {
  if(c[i]==1) printf(" + x[%d]",i);
  else if(c[i]==-1) printf(" - x[%d]",i);
  else if(c[i]>=0) printf(" +%lg x[%d]",c[i],i);
  else printf(" %lg x[%d]",c[i],i);
}
printf("\nst:\n");
for(i=1;i<=num_st;i++) { 
  k=0;
  for(j=1;j<=num_x;j++) {
   if(a[i][j]!=0) {
    if(a[i][j]==1&&k!=0) printf(" + x[%d]",j);
    else if(a[i][j]==1&&k==0) printf("  x[%d]",j);
    else if(a[i][j]==-1) printf(" - x[%d]",j);
    else if(a[i][j]>=0&&k!=0) printf(" +%lg x[%d]",a[i][j],j);
    else if(a[i][j]>=0&&k==0) printf(" %lg x[%d]",a[i][j],j);
    else printf(" %lg x[%d]",a[i][j],j);
    k=1;
   }
  }
  switch(re_st[i]) {
   case 1: 
    printf(" == %lg\n",b[i]); 
    break;
   case 2: 
    printf(" >= %lg\n",b[i]); 
    break;
   case 3: 
    printf(" <= %lg\n",b[i]); 
    break;
  }
}
printf(" x[1]~x[%d]>=0\n",num_x);
tnum_x=num_x;
for(i=1;i<=num_st;i++) {
  switch(re_st[i]) {
  case 1:
  case 3:
   num_x+=1;
   break;
  case 2:
   num_x+=2;
   break;
  }
}
//化为标准形式
if(ma_mi==2) for(i=1;i<=tnum_x;i++) c[i]*=-1; //求最小值时，系数变相反数
for(i=1;i<=num_st;i++) {
  switch(re_st[i]) {
   case 1:
    num_addv++;
    num_ba++;
    num_ar++;
    c[tnum_x+num_addv]=0; 
    base[num_ba]=arti[num_ar]=tnum_x+num_addv;
    for(j=tnum_x+1;j<=num_x;j++) 
     if(j==tnum_x+num_addv) a[i][tnum_x+num_addv]=1;
     else a[i][j]=0;
    break;
   case 2:
    num_addv++;
    c[tnum_x+num_addv]=0;
    num_addv++;
    num_ba++;
    num_ar++;
    c[tnum_x+num_addv]=0; 
    base[num_ba]=arti[num_ar]=tnum_x+num_addv;
    for(j=tnum_x+1;j<=num_x;j++) 
     if(j==tnum_x+num_addv-1) a[i][tnum_x+num_addv-1]=-1;
     else if(j==tnum_x+num_addv) a[i][tnum_x+num_addv]=1;
     else a[i][j]=0;
    break;
   case 3:
    num_addv++;
    num_ba++;
    c[tnum_x+num_addv]=0;
    base[num_ba]=tnum_x+num_addv;
    for(j=tnum_x+1;j<=num_x;j++) 
     if(j==tnum_x+num_addv) a[i][tnum_x+num_addv]=1;
     else a[i][j]=0;
    break;
  }//switch
}//增加松弛变量、剩余变量、人工变量、确定基变量
//显示标准化后的方程组
printf("\n化为标准形式后：\n\n");
if(ma_mi==1) printf("max z= ");
else printf("max z'= ");
printf("%lg x[%d]",c[1],1);
for(i=2;i<=num_x;i++) { 
  k=0;
  for(j=1;j<=num_ar;j++)
  if(i==arti[j]) k=1;
  if(k==1) printf(" -M x[%d]",i);
  else if(c[i]==1) printf(" + x[%d]",i);
  else if(c[i]==-1) printf(" - x[%d]",i);
  else if(c[i]>=0) printf(" +%lg x[%d]",c[i],i);
  else printf(" %lg x[%d]",c[i],i);
}
printf("\nst:\n");
for(i=1;i<=num_st;i++) {
  k=0;
  for(j=1;j<=num_x;j++) {
   if(a[i][j]!=0) {
    if(a[i][j]==1&&k!=0) printf(" + x[%d]",j);
    else if(a[i][j]==1&&k==0) printf("  x[%d]",j);
    else if(a[i][j]==-1) printf(" - x[%d]",j);
    else if(a[i][j]>=0&&k!=0) printf(" +%lg x[%d]",a[i][j],j);
    else if(a[i][j]>=0&&k==0) printf(" %lg x[%d]",a[i][j],j);
    else printf(" %lg x[%d]",a[i][j],j);
    k=1;
   }
  }
  printf(" == %lg\n",b[i]); 
}
printf(" x[1]~x[%d]>=0",num_x);
}
void iterative() {
int i,j,k,k_a,k_f,l; //k_a,k_f值为0或1，记录当前下标在arti[]或base[]里的搜索结果
int base_elem;
int base_out,base_in;
double sigma[MAX],temp;
double value_be; //高斯消元里保存主元素值
printf("\n\n第%d次迭代：\n\n",stop); 
for(i=1;i<=num_st;i++) {
  printf("c%d=%lg\t",base[i],c[base[i]]);
  printf("b%d=%lg\t",i,b[i]);
  switch(step) {
   case 1:
    for(j=1;j<=num_x;j++)
    {
     printf("a[%d][%d]=%lg\t",i,j,a[i][j]);
    }
    printf("\n");
    break;
   case 2:
    for(j=1;j<=num_x;j++) {
     k_a=0;
     for(l=1;l<=num_ar;l++) if(j==arti[l])k_a=1;
     if(k_a!=1) printf("a[%d][%d]=%lg\t",i,j,a[i][j]);
    }
    printf("\n");
    break;
  }
}
//求检验数sigma
for(i=1;i<=num_x;i++) {
  sigma[i]=c[i];
  for(j=1;j<=num_st;j++) sigma[i]-=c[base[j]]*a[j][i];
  for(j=1;j<=num_st;j++) if(i==base[j]) sigma[i]=0;
  switch(step) {
   case 1:
    printf("sigma[%d]=%lg\t",i,sigma[i]);
    break;
   case 2:
    k_a=0;
    for(l=1;l<=num_ar;l++) if(i==arti[l]) k_a=1;
    if(k_a!=1) printf("sigma[%d]=%lg\t",i,sigma[i]);
    break;
  }
}
putchar('\n');
//检验检验数sigma是否全小于等于0
k=0;
for(i=1;i<=num_x;i++) {
  if(sigma[i]>0) 
  k=1;
}
if(k==0) {
  //sigma是全小于等于0时，检查是否为无穷多最优解
  for(i=1;i<=num_x;i++) { 
   k_f=k_a=0;
   for(j=1;j<=num_ar;j++)
   if(i==arti[j]) k_a=1;
   if(sigma[i]==0&&k_a!=1) {
    for(j=1;j<=num_st;j++) if(i==base[j]) k_f=1;
   if(k_f==0) {status=-1; return;}
   }
  }
  status=1; 
  return;
}
//检查是否为无界解
for(i=1;i<=num_x;i++) { 
  k_f=0;
  if(sigma[i]>0) {
   for(j=1;j<=num_st;j++) if(a[j][i]>0) k_f=1;
   if(k_f!=1) {status=0; return;}
  }
}
//确定换入变量
for(i=1;i<=num_x;i++) { 
  k=0;
  for(j=1;j<=num_st;j++) if(i==base[j]) k=1;
  if(k==0&&sigma[i]>0) temp=sigma[i]-1;
}//temp赋初值
for(i=1;i<=num_x;i++) {
  k=0;
  for(j=1;j<=num_st;j++) if(i==base[j]) k=1;
  if(k==0)
   if(sigma[i]>temp&&sigma[i]>0) {
    base_in=i;
    temp=sigma[i];
   }
}
//确定换出变量
for(i=1;i<=num_st;i++) 
  if(a[i][base_in]>0) {
   temp=b[i]/a[i][base_in]+1;
   break;
  }//temp赋初值
for(i=1;i<=num_st;i++) {
  if(b[i]/a[i][base_in]<=temp&&a[i][base_in]>0) {
   for(j=1;j<=num_ar;j++)
    if(base[i]==arti[j]) { 
     base_out=base[i];
     base_elem=i;
     temp=b[i]/a[i][base_in];
     break;
    }
  }//人工变量优先换出
  if(b[i]/a[i][base_in]<temp&&a[i][base_in]>0) {
   base_out=base[i];
   base_elem=i;
   temp=b[i]/a[i][base_in];
  }
}
printf(" 基变量：");
for(i=1;i<=num_st;i++) printf("x[%d] ",base[i]);
printf("换入变量：x[%d] 换出变量：x[%d]",base_in,base_out);
//基变量变换，进行新方程初始化后迭代
for(i=1;i<=num_st;i++) {
  if(base[i]==base_out) base[i]=base_in;
}
//初始化主元素行系数
value_be=a[base_elem][base_in];
b[base_elem]/=value_be;
for(i=1;i<=num_x;i++) a[base_elem][i]/=value_be;
for(i=1;i<=num_st;i++) {
  if(i!=base_elem) {
   b[i]-=b[base_elem]*a[i][base_in];
   value_be=a[i][base_in];
   for(j=1;j<=num_x;j++) a[i][j]-=a[base_elem][j]*value_be;
  }
}
stop++;
if(stop>STP) {status=-2; return;}
iterative();
}
void output() {
int i,j;
double X[MAX];
printf("\n结果如下：\n");
printf("\nX=(");
for(i=1;i<=num_x;i++) {
  for(j=1;j<=num_st;j++)
  if(i==base[j]) {X[i]=b[j];break;}
  else X[i]=0;
  printf("%lg ",X[i]);
}
printf(")");
for(i=1;i<=num_x;i++) max+=c[i]*X[i];
if(ma_mi==1) printf("\nMax z= %lf\n",max);
else printf("\nMin z= %lf\n",-max);
}