📄 greatm.c
字号:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// 大M单纯形法+参数b灵敏度测试
// Author: lubangjian
// Date: 2003/12/07
#define SUCCESS 1
#define FAIL 0
#define GREATM 0xffffffff//大M值
/*分配内存*/
void *_malloc(char *position,int memsize)
{
void *mem;
if ((mem=malloc(memsize))==NULL)
{
printf("%s:内存分配失败\n",position);
return NULL;
}
return mem;
}
//功能说明:取进基变量
//入参说明:sigma--判别数列表
// p--限制条件参数矩阵
// row--矩阵列数
// col--矩阵行数
//返 回 值:>0--进基变量列号,-1--找不到进基变量,0--当前已为最优解
int get_ent_para(double *sigma,double *p,int row,int col)
{
int colcount,rowcount,colmark;
double mino;
//取最小判别数值
for (colcount=1,mino=sigma[0],colmark=0;colcount<col;colcount++)
if (mino>sigma[colcount])
{
mino=sigma[colcount];
colmark=colcount;
}
if ((mino>=0.0)||(mino==-0.0))//当前的判别数均大于等于0,即当前为最优解
return -2;
//检查是否进基变量系数矩阵所在列是否全小于0
for (rowcount=0;rowcount<row;rowcount++)
{
if (p[rowcount*col+colmark]>0.0)
break;
}
if (rowcount>=row)//当前不可能有最优解
return -1;
return colmark;//返回进基变量所在列号
}
//功能说明:计算离基变量值,并由此确定离基变量
//入参说明:entp--进基变量参数
// p--限制条件参数矩阵
// b--基解
// sitar--离基变量判别数存放区
// row--矩阵列数
// col--矩阵行数
int get_out_para(int entp,double *p,double *b,double *sitar,int row,int col)
{
int rowcount;
double minsitar;
int rowmark;
//计算离基变量判别数值
for (rowcount=0;rowcount<row;rowcount++)
{
if (p[rowcount*col+entp]<=0)
{
sitar[rowcount]=-1.0;//表示该项无效
continue;
}
sitar[rowcount]=b[rowcount]/p[rowcount*col+entp];
}
//取离基变量
for (rowcount=1,minsitar=sitar[0],rowmark=0;rowcount<row;rowcount++)
{
if (sitar[rowcount]<=-1.0)
{
continue;
}
if ((minsitar>sitar[rowcount])||(minsitar<=-1.0))
{
rowmark=rowcount;
minsitar=sitar[rowcount];
}
}
return rowmark;//返因离基变量所在行号
}
//功能说明:输出最优化结果
int output(int *B,double *b,double *cb,int row,int col,int *greatmlist,double *solve,int printres)
{
int rowcount,colcount,greatMcount;
double funcval;
char *position="output";
for (colcount=0;colcount<col;colcount++)
solve[colcount]=0.0;
for (rowcount=0,funcval=0;rowcount<row;rowcount++)
{
if (b[rowcount]<0)
{
if (printres)
printf("无最优解!\n");
return FAIL;
}
solve[B[rowcount]]=b[rowcount];//最优解
funcval+=cb[rowcount]*b[rowcount];
}
for (colcount=0;colcount<col;colcount++)
if ((greatmlist[colcount])&&(solve[colcount]>0.0))
{
if (printres)
printf("无最优解!\n");
return FAIL;
}
if (printres)//需要输出结果
{
printf("最优解:(");
for (colcount=0;colcount<col;colcount++)
{
if (greatmlist[colcount]==0)
{
if (colcount>0)//此处默认为第一个变量不可能为大M
printf(",");
printf("%f",solve[colcount]);
}
}
printf(")\n函数值:%f\n",funcval);
}
return SUCCESS;
}
//功能说明:输出中间结果
int output_mid(double *cj,
double *b,
int *B,
double *p,
int row,
int col,
double *sigma,
double *cb,
double *cjo,
double *sitar)
{
int rowcount,colcount;
printf("\tcj \t|");
for (colcount=0;colcount<col;colcount++)
printf(" %f",cj[colcount]);
printf("\n");
for (rowcount=0;rowcount<row;rowcount++)
{
printf("%f|P%d|%f|",cb[rowcount],B[rowcount]+1,b[rowcount]);
for (colcount=0;colcount<col;colcount++)
printf(" %f",p[rowcount*col+colcount]);
printf(" |%f\n",sitar[rowcount]);
}
printf("\n\tsigma \t|");
for (colcount=0;colcount<col;colcount++)
printf(" %f",sigma[colcount]);
printf("\n");
return 0;
}
#define Free_Alocated_mem if (sigma) free(sigma);\
if (cb) free(cb);\
if (cjo) free(cjo);\
if (sitar) free(sitar);\
if (greatmlist) free(greatmlist);\
if (b) free(b);\
if (B) free(B);\
if (p) free(p)
//功能说明:大M单纯形法
//入参说明:cj--原形参数列表
// ent_b--初始基解
// ent_B--基变量参数
// ent_p--限制条件参数矩阵
// row--矩阵列数
// col--矩阵行数
// printres--是否打印结果
// needmid--是否需要中间结果
// solve--解
//返 回 值:SUCCESS--成功,FAIL--失败
int greatM(double *cj,double *ent_b,int *ent_B,double *ent_p,int row,int col,int printres,int needmid,double *solve)
{
double *sigma=NULL;//判别数
double *cb=NULL;
double *cjo=NULL;
int entp,outp;//进基变量,离基变量
double *sitar=NULL;//离期变量判别虽数
double tempo,tempb;
int res;
int *greatmlist=NULL;
int colcount,rowcount;
double greatm=GREATM;
char *position="greatM";
double *b=NULL,*p=NULL;
int *B=NULL;
if (!(b=(double *)_malloc(position,sizeof(double)*row)))
return FAIL;
if (!(B=(int *)_malloc(position,sizeof(int)*row)))
{
Free_Alocated_mem;
return FAIL;
}
if (!(p=(double *)_malloc(position,sizeof(double)*row*col)))
{
Free_Alocated_mem;
return FAIL;
}
if (!(sigma=(double *)_malloc(position,sizeof(double)*col)))//为判别数分配空间
{
Free_Alocated_mem;
return FAIL;
}
if (!(cb=(double *)_malloc(position,sizeof(double)*row)))//为CB分配空间
{
Free_Alocated_mem;
return FAIL;
}
if (!(sitar=(double *)_malloc(position,sizeof(double)*row)))//为离基变量判别数分配空间
{
Free_Alocated_mem;
return FAIL;
}
if (!(greatmlist=(int *)_malloc(position,sizeof(int)*col)))//为可能的大M空间分配内存
{
Free_Alocated_mem;
return FAIL;
}
memcpy(b,ent_b,sizeof(double)*row);
memcpy(B,ent_B,sizeof(int)*row);
memcpy(p,ent_p,sizeof(double)*row*col);
for (rowcount=0;rowcount<row;rowcount++)
cb[rowcount]=cj[B[rowcount]];
//取大M列表
for (colcount=0;colcount<col;colcount++)
{
if (cj[colcount]>=greatm)
greatmlist[colcount]=1;
else
greatmlist[colcount]=0;
}
//计算初始判别数
for (colcount=0;colcount<col;colcount++)
{
tempo=0.0;
for (rowcount=0;rowcount<row;rowcount++)//计算一列的判别数
tempo+=p[rowcount*col+colcount]*cj[B[rowcount]];
sigma[colcount]=cj[colcount]-tempo;
}
while(1)//逐次迭代取最优解
{
//检查基解
for (rowcount=0;rowcount<row;rowcount++)
if (b[rowcount]<0)
{
if (printres)
printf("无最优解!\n");
return FAIL;
}
res=get_ent_para(sigma,p,row,col);//取进基变量
if (res==-1)//查找进基变量失败
return FAIL;
if (res==-2)//当前为最优解
{
res=output(B,b,cb,row,col,greatmlist,solve,printres);
Free_Alocated_mem;
return res;
}
entp=res;//取进基变量
outp=get_out_para(entp,p,b,sitar,row,col);//计算离基变判别数,取得离基变量
if (needmid)//需要输出中间结果
{
printf("entp=%d,outp=%d\n",entp+1,outp+1);
output_mid(cj,b,B,p,row,col,sigma,cb,cjo,sitar);
}
B[outp]=entp;//更新B值
cb[outp]=cj[entp];//更新cb值
//计算非离基行参数矩阵值
tempb=b[outp]/p[outp*col+entp];
for (rowcount=0;rowcount<row;rowcount++)
{
if (rowcount==outp)//当前为离基变量所在行
{
b[rowcount]/=p[outp*col+entp];/*计算与离基变量在同行的b值*/
continue;
}
b[rowcount]-=tempb*p[rowcount*col+entp];/*计算与离基变量不同行b值*/
tempo=p[rowcount*col+entp]/p[outp*col+entp];
for (colcount=0;colcount<col;colcount++)
p[rowcount*col+colcount]-=p[outp*col+colcount]*tempo;
}
//计算系数矩阵离基量所在行的值和新的判别数
tempo=sigma[entp]/p[outp*col+entp];
for (colcount=0;colcount<col;colcount++)
{
sigma[colcount]-=p[outp*col+colcount]*tempo;//新的判别数值
if (colcount!=entp)
p[outp*col+colcount]/=p[outp*col+entp];//与进基变量不在同列的系数矩阵值
}
p[outp*col+entp]=1.0;//
}
Free_Alocated_mem;
return FAIL;
}
#define delicacy 15
#define REALCJ 5
#define DOWN 0
#define UP 1
#define STRATDELICACY 100000 //起始精度
#define get_an_area(direction) from=b2[bcount];\
iscontinue=1;\
while(iscontinue)\
{\
if (direction==UP)\
b2[bcount]+=delicacyval;\
else\
b2[bcount]-=delicacyval;\
if (greatM(cj,b2,B,*p,3,7,0,0,tempsolve)==SUCCESS)/*当前取到最优值*/\
{\
for (rescount=0;rescount<REALCJ;rescount++)\
if (((solve[rescount]>0.0)&&(tempsolve[rescount]<=0.0))\
||((solve[rescount]<=0)&&(tempsolve[rescount]>0.0)))/*当前基解已发生变化*/\
{\
iscontinue=0;/*停止查找*/\
break;\
}\
if (iscontinue)\
from=b2[bcount];/*调整起始点*/\
}\
else\
{\
iscontinue=0;\
}\
if ((direction==UP)&&(b2[bcount]>=unlimit))/*认为无上限*/\
iscontinue=0;\
if ((direction==DOWN)&&(b2[bcount]<=-unlimit))/*认为无下限*/\
iscontinue=0;\
}
//调用例子
int main(int argc,char **argv)
{
int needmid=0;//默认为不需中间结果
double cj[7]={-3,1,1,0,0,GREATM,GREATM},solve[7],tempsolve[7];
double b[3]={11,3,1},bi,step,stepcount,from,to,delicacyval=1.0,b2[3],limit[2];
int B[3]={4,5,6},bcount,delicacycount,rescount,iscontinue;
double p[3][7]={{1,-2,1,0,1,0,0},{-4,1,2,-1,0,1,0},{-2,0,1,0,0,0,1}};
double unlimit=GREATM;
if (argc>1)
needmid=atoi(argv[1]);
if (greatM(cj,b,B,*p,3,7,1,needmid,solve)==SUCCESS)//当前取得最优解
{
printf("约束右边值b灵敏度分析(当测试绝对值>=%f时,不再测试!)\n",unlimit);
for (bcount=0;bcount<3;bcount++)//对b逐个进行灵敏度分析
{
limit[0]=limit[1]=b[bcount];
delicacyval=STRATDELICACY;//从精度1开始分析
for (delicacycount=0;delicacycount<delicacy;delicacycount++)
{
memcpy(b2,b,sizeof(double)*3);
b2[bcount]=limit[1];
get_an_area(UP);//取上限
limit[1]=from;
memcpy(b2,b,sizeof(double)*3);
b2[bcount]=limit[0];
get_an_area(DOWN);//取下限
limit[0]=from;
delicacyval/=10;//精度提高10
printf("有效范围:(%.10f,%.10f),测试精度:%10.10f\n",limit[0],limit[1],delicacyval);
}
printf("\n最终范围:");
if (limit[0]<=-unlimit)
printf("b%d:(无下限,",bcount+1);
else
printf("b%d:(%.10f,",bcount+1,limit[0]);
if (limit[1]>=unlimit)
printf("无上限),测试精度:%10.10f\n\n",delicacyval);
else
printf("%.10f),测试精度:%10.10f\n\n",limit[1],delicacyval);
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -