单纯形解决线形规划问题的代码(c)
总算完成了,下面是代码。 又:因为运筹学并没有学过,这次凑巧帮管院一个同学写单纯形的代码, 所以就有了这个。是从昨天开始看书的,如果有错误,还请见谅------觉得那个运筹学还真不容易呢,虽比不上咱们的离散:)~~~~~~~~
测试结果: Successful: X[3]=77.142860 X[1]=7.142854 X[2]=28.571430 max Z=585.714294 Press any key to continue
/* *data.txt */
5 2 1 0 0 2 3 0 1 0 1 5 0 0 1
170 100 150
10 18 0 0 0
1 0 0 0 1 0 0 0 1
1 0 0 0 1 0 0 0 1
0 0 0
2 3 4
/* 本题的例子是:
max Z=10*x1+18*x2; s.t. 5*x1+2*x2<=170 2*x1+3*x2<=100 x1+5*x2<=150 x1,x2>=0
本文件的数据是上式经过加入松弛变量后的标准形式的数据
*/
/* *fmain.cpp */ #include <stdio.h> #include <stdlib.h> #define M 3 #define N 5
/* 解决形如 max z=cx s.t. Ax=b, x>=0 的线形规划问题
_B[M][M]是B[M][M]的逆阵 */ float A[M][N],b[M], c[N],B[M][M], cn[M],X[M],Y[M], Kao[N],P[M], E[M][M],_B[M][M]; /* 解序号 */ int Xcount[M]; float sita; /* 进基变量,出基变量 */ int xin,xout;
int init(FILE *pf); void computeE(); void compute_B(); void computeXY(); int computeKao(); int computePsita(); void computeBcn(); int succ(); int fall();
void main() { FILE *pf; float result=0.0; if(!init(pf)) { printf("Error:can't open data file,please check...\n"); exit(-1); } computeXY(); xin=computeKao(); while(!succ()) { xout=computePsita(); if(fall()) { printf("Fall...already exit...\n"); exit(0); } computeBcn(); computeE(); compute_B(); computeXY(); xin=computeKao(); } printf("Successful:\n"); for(int i=0;i<M;i++) printf("X[%d]=%f\t",Xcount[i]+1,X[i]);
for(int i=0;i<M;i++) result+=c[Xcount[i]]*X[i]; printf("\nmax Z=%f\n",result); } /* 用data.txt文件初始化,包括: A,b,c,B,_B,cn,Xcount */ int init(FILE *pf) { pf=fopen("data.txt","r"); if(pf==NULL){ fclose(pf); return 0; } /*input A[M][N]*/ for(int i=0;i<M;i++) for(int j=0;j<N;j++){ fscanf(pf,"%f",&A[i][j]); }
/*input b[M]*/ for(int i=0;i<M;i++){ fscanf(pf,"%f",&b[i]); }
/*input c[N]*/ for(int i=0;i<N;i++){ fscanf(pf,"%f",&c[i]); }
/*input B[M][M]*/ for(int i=0;i<M;i++) for(int j=0;j<M;j++) fscanf(pf,"%f",&B[i][j]);
/*input _B[M][M]*/ for(int i=0;i<M;i++) for(int j=0;j<M;j++){ fscanf(pf,"%f",&_B[i][j]); }
/*input cn[M]*/ for(int i=0;i<M;i++) fscanf(pf,"%f",&cn[i]);
/* input Xcount[M] */ for(int i=0;i<M;i++){ fscanf(pf,"%d",&Xcount[i]); }
/*其它*/ for(int i=0;i<M;i++) P[i]=X[i]=Y[i]=0;
for(int i=0;i<N;i++) Kao[i]=0;
fclose(pf); return 1; } /* 计算初等变换矩阵E */ void computeE() { for(int i=0;i<M;i++) for(int j=0;j<M;j++) { E[i][j]=0; if(i==j) E[i][j]=1; }
for(int i=0;i<M;i++) E[i][xout]=-P[i]/P[xout]; E[xout][xout]=1/P[xout]; } /* 计算B逆矩阵_B */ void compute_B() { float __B[M][M]; for(int i=0;i<M;i++) for(int j=0;j<M;j++) __B[i][j]=0;
for(int i=0;i<M;i++) for(int j=0;j<M;j++) for(int k=0;k<M;k++) __B[i][j]+=E[i][k]*_B[k][j];
for(int i=0;i<M;i++) for(int j=0;j<M;j++){ _B[i][j]=__B[i][j]; }
}
/* 计算 X=_B*b ,Y=_B*cn */ void computeXY() { for(int i=0;i<M;i++) P[i]=X[i]=Y[i]=0;
for(int i=0;i<M;i++){ for(int j=0;j<M;j++) X[i]+=_B[i][j]*b[j]; } for(int i=0;i<M;i++){ for(int j=0;j<M;j++) Y[i]+=cn[j]*_B[j][i]; }
} /* 计算检验数 kao=c-cn*_B*A xin<-return */ int computeKao() { float midval=0.0; int k; for(int i=0;i<N;i++){ for(int j=0;j<M;j++) midval+=Y[j]*A[j][i]; Kao[i]=c[i]-midval; midval=0; } midval=Kao[0]; k=0; for(int j=0;j<N;j++) if(Kao[j]>midval){ midval=Kao[j]; k=j; }
return k; } /*计算进基列向量P,判断出基变量 xout<-return */ int computePsita() { int k=0; for(int i=0;i<M;i++){ for(int j=0;j<M;j++) P[i]+=_B[i][j]*A[j][xin]; }
sita=X[0]/P[0]; for(int j=0;j<M;j++) if(sita>X[j]/P[j]) { sita=X[j]/P[j]; k=j; }
return k; } /* 计算矩阵B,cn */ void computeBcn() { for(int i=0;i<M;i++) B[i][xout]=A[i][xin]; cn[xout]=c[xin];
Xcount[xout]=xin; } /* 判断是否成功 */ int succ() { int k=1; for(int i=0;i<N;i++) k=k&&(Kao[i]<=0.0);
return k; } /* 算法失败,退出 */ int fall() { int k=1; for(int i=0;i<N;i++) k=k&&(P[i]<=0.0); return k; }

|