work1.cpp

本代码实现了系统工程中线性规划的单纯型解法和灵敏度分析

#include <stdio.h>
#include <stdlib.h>

int t=1; //总标志位，如果等于0说明无解

struct u
{
	int m;//行
	int n;//列
	float data;//数值
	struct u * r;
	struct u * l;
}* Basic;
struct line_h
{
	int s;      //行中总元素个数
	int m;      //行数
	int tag;    //1表示已经存在旋转主元，基变量个数不变，0则反之。
	int bac;    //一行中，基变量的下标数
	struct unit_p * u;   //指向此行数据节点的指针
	struct line_h * next;  //next指针
} * chart, * chart1;

struct unit_p
{
	int n;     //列数
	float data;  //数据
	struct unit_p * next; //next指针
};


void div (struct line_h * p,float a)    //把p指针所指的行元素都除以a
{
	struct unit_p *q;
	
	q=p->u;

	while (q!=NULL)
	{
		q->data/=a;
		q=q->next;
	}

}

void output_data()      // 输出最终答案
{
	struct unit_p *q;
	struct line_h *p, *p1;
	int c;
	int m;
	int bac=0,s_bac=0;

	FILE * fp;

	c=chart->s;
	p=chart;
	
	while(p!=NULL){m=p->m;p=p->next;}

	fp=fopen("out.txt","a");

	fprintf(fp,"\n当");

	while(m>0)
	{
		p=chart;

		while (p!=NULL)//找到最小的基变量
		{
			if (p->bac>s_bac)
				if (bac==0||bac>p->bac) {bac=p->bac;p1=p;}				
			p=p->next;
		}
		
		q=p1->u;
		while(q->next!=NULL) q=q->next;

		fprintf(fp,"x%d=%.3f,",bac,q->data);

		s_bac=bac;
		bac=0;
//		p1->bac=0;
		p1=NULL;
		q=NULL;

		m--;
	}
	
	q=chart->u;
	while(q->n!=c){q=q->next;}
	fprintf(fp,"其他为零时,最小值是 %.3f",-(q->data));

	fclose(fp);
}

void output_eror(int i)   //无解时输出错误信息
{
	FILE * fp;
	
	fp=fopen("out.txt","a");

	if (i==1) fprintf(fp,"\n没有可行解!");
	if (i==2) fprintf(fp,"\n目标函数无下界!");

	fclose(fp);
}

void add (struct line_h *p1,struct line_h *p2,int k)  //以p2为基准，单位化第k列
{
	struct unit_p *q1,*q2,*q3;
	float a;

	q1=p1->u;
	q2=p2->u;
	
	q3=q1;
	while(q3->n!=k) {q3=q3->next;}

	a=q3->data;

	while (q1!=NULL)
	{
		q1->data = q1->data - q2->data * a;
		q1=q1->next;
		q2=q2->next;
	}

}


void output_arry(int i)  //输出表中所有数据
{
	FILE * fp;
	struct line_h *p;
	struct unit_p *q;

	p=chart;

	fp=fopen("out.txt","a");

	if (fp==NULL) fp=fopen("out.txt","w");


	if (i==1) fprintf(fp,"初始矩阵-----------------------\n");
	if (i==2) fprintf(fp,"初始可行基---------------------\n");
	if (i==3) fprintf(fp,"换基---------------------------\n");

	
	while(p!=NULL)
	{
		q=p->u;
		while(q!=NULL)
		{
			fprintf(fp,"%8.3f",q->data);
			q=q->next ;
		}
		fprintf(fp,"\n");
		p=p->next;
	}

	fclose(fp);

}

int compute_q(int k)  //计算第k列的最小比值，并且返回所在的行数
{
	float a,b,min;
	int cont,c;
	struct line_h * p;
	struct unit_p * q;

	p=chart->next;
	min=-1;
	cont=0;

	while (p!=NULL)
	{
		c=p->s;
		q=p->u;

		while(q!=NULL)
		{
			if (q->n==k) a=q->data;
			if (q->n==c) b=q->data;
			q=q->next;
		}
		
		if (a>0)
		{
			if (min==-1) {min=b/a;cont=p->m;}
			else if (min>b/a)
				{
					min=b/a;
					cont=p->m;
				}
		}

		p=p->next;
	}
	
	return(cont);
}

void retort (int k,int l) //以l行k列为旋转主元，旋转矩阵。
{
	struct line_h * p,* p1;
	struct unit_p * q;
	float a;
	
	p=chart;
	while (p->m!=l) {p=p->next;}

	q=p->u;
	while (q->n!=k) {q=q->next;}

	if (q->data!=1)
	{
		a=q->data;
		div(p,a);
	}

	p1=chart;
    while (p1!=NULL)
	{
		if (p1!=p) add(p1,p,k);
		p1=p1->next;
	}
	
	p=chart->next;
	while(p->m!=l) p=p->next;
	p->bac=k;

}

void initial()    //求初始可行基
{
	int k,l;
	int s,m=0;
	int cont=0;
	struct line_h * p;
	

	s=chart->s;

	k=1;

	p=chart;
	while (p!=NULL){m++;p=p->next;}

	while (k<s&&cont<m)
	{
		l=compute_q(k);
		
		if (l!=0)
		{
			retort(k,l);
//			output_arry();
			
			p=chart;
			while(p->m!=l){p=p->next;}

			if (p->tag==0) {cont++;p->tag=1;}

		}

		k++;
	}



/*
	while (org->next!=NULL&&k<=s)
	{
		l=compute_q(k);

		if (l!=0)
		{

		retort(k,l);
		cont++;

		p=org;
		while (p->m!=l) {p=p->next;}

		p1=org;
		while (p1->m!=l-1) {p1=p1->next;}

		p1->next=p->next;
		p->next=NULL;
		
		if (fin==NULL) fin=p;
		else {p->next=fin;fin=p;}

		p=org;
		m=0;
		while(p!=NULL)
		{
			p->m=m;
			m++;
			p=p->next;
		}
		}  
		k++;

	}

	org->next=fin;

	m=0;
	while(org!=NULL)
	{
		org->m=m;
		m++;
		org=org->next; 
	}
*/
	output_arry(2);

	if (cont!=m-1) 
	{
//		printf("There is no available answer!");
		output_eror(1);
		t=0;
	}
}



void change()  //换基
{
	struct line_h * p;
	struct unit_p * q;
	int l,c;

	p=chart;
	q=p->u;
	c=p->s;

	while (q!=NULL&&q->n<c)
	{
		if (q->data<0) 
		{
			l=compute_q(q->n);
			if (l!=0) {retort(q->n,l);q=p->u;output_arry(3);}
			else 
			{
				output_eror(2);
				t=0;
				q=NULL;
			}
//			printf("There is no limit answer!");

		}
		else q=q->next;
	}



}





void input_data() //从文件中读入数据
{
	FILE *fp;
	float a;
	char c;
	struct line_h * p1,* p2;
	struct unit_p * q1,* q2,* q3;
	int m=0,n=1; //行初始0列初始1

	q1=NULL;
	q2=NULL;
	p1=NULL;
	p2=NULL;

	
	fp=fopen("dat.txt","r");

	do
	{
		fscanf(fp,"%f",&a);
		c=fgetc(fp);
		

		q1=(struct unit_p *)malloc(sizeof(struct unit_p));
		q1->data=a;
		q1->n=n++;
		q1->next=NULL;

		if (q2==NULL) {q2=q1;q3=q1;}
		else
		{
			q2->next=q1;
			q2=q1;
		}

		if (c=='\n'||c=='#')
		{
			p1=(struct line_h *)malloc(sizeof(struct line_h));
			p1->m=m++;
			p1->next=NULL;
			p1->u=q3;
			p1->s=n-1;
			p1->tag=0;
			p1->bac=0;

			if (p2==NULL) {p2=p1;chart=p1;}
			else
			{
				p2->next=p1;
				p2=p1;
			}

			n=1;
			q1=NULL;q2=NULL;q3=NULL;
			
		}




	}while(c!='#');



	fclose(fp);


	output_arry(1);


}

int IsBasicP(int i)
{
	struct line_h * p;
	p=chart->next;
	int j=0;
	
	while (p!=NULL)
	{
		if (p->bac==i) j=1;
		p=p->next;
	}
	return (j);
}

void analysis_c()//目标函数中价格系数的灵敏度分析
{
	int i,n;
	float c,a,cmax=0,cmin=0;
	struct line_h * p;
	struct unit_p * q, *q1;

	n=chart->s;
	
	
	FILE * fp;


	fp=fopen("out.txt","a");
	fprintf(fp,"\n\n目标函数中价格系数的灵敏度分析结果：");
	fclose(fp);

	for (i=1;i<n;i++)
	{
		if (IsBasicP(i))
		{
			p=chart->next;//找到基变量所对应的行；
			while (p->bac!=i) p=p->next;
			
			q=p->u;
			q1=chart->u;
			
			while(q->n<n)
			{
				if (q->data!=0&&q->data!=1) 
				{
					a=q1->data/q->data;
					if (a<0&&(a>cmax||cmax==0)) cmax=a;
					if (a>0&&(a<cmin||cmin==0)) cmin=a;
				}
				q1=q1->next;
				q=q->next;
			}
			
			fp=fopen("out.txt","a");

			fprintf(fp,"\n%.3f<=C%d<=%.3f",cmax,i,cmin);
			fclose(fp);

			cmax=0;cmin=0;

		}
		else
		{
			q=chart->u;
			while (q->n!=i) q=q->next; 
			
			if (q->data!=0) c=-1*q->data;
			else c=0;

			fp=fopen("out.txt","a");

			fprintf(fp,"\n%.3f<=C%d",c,i);
			fclose(fp);

		}
	}

}

void input_data1() //从文件中读入数据
{
	FILE *fp;
	float a;
	char c;
	struct line_h * p1,* p2;
	struct unit_p * q1,* q2,* q3;
	int m=0,n=1; //行初始0列初始1

	q1=NULL;
	q2=NULL;
	p1=NULL;
	p2=NULL;

	
	fp=fopen("dat.txt","r");

	do
	{
		fscanf(fp,"%f",&a);
		c=fgetc(fp);
		

		q1=(struct unit_p *)malloc(sizeof(struct unit_p));
		q1->data=a;
		q1->n=n++;
		q1->next=NULL;

		if (q2==NULL) {q2=q1;q3=q1;}
		else
		{
			q2->next=q1;
			q2=q1;
		}

		if (c=='\n'||c=='#')
		{
			p1=(struct line_h *)malloc(sizeof(struct line_h));
			p1->m=m++;
			p1->next=NULL;
			p1->u=q3;
			p1->s=n-1;
			p1->tag=0;
			p1->bac=0;

			if (p2==NULL) {p2=p1;chart1=p1;}
			else
			{
				p2->next=p1;
				p2=p1;
			}

			n=1;
			q1=NULL;q2=NULL;q3=NULL;
			
		}




	}while(c!='#');



	fclose(fp);
}


void analysis_b()
{

	
	
	input_data1();//重新读入数据到chart1

	float *a;
	int m=0;
	struct line_h* p;
	struct unit_p* q;
	int bac=0,s_bac=0;

	int i=0,j;

	FILE * fp;



	fp=fopen("out.txt","a");
	fprintf(fp,"\n\n资源系数灵敏度分析过程：");
	fclose(fp);

//计算基的个数
	p=chart1->next;
	while (p!=NULL) {m++;p=p->next;}

//读取基矩阵
	a=new float[m*m];
		
	for (j=0;j<m;j++)
	{
		p=chart->next;//依次查找基在的列
		while(p!=NULL)
		{
			if (p->bac>s_bac&&(bac==0||bac>p->bac)) bac=p->bac;
			p=p->next;
		}

		p=chart1->next;//读入数据到a
		while(p!=NULL)
		{
			q=p->u;
			while(q->n!=bac) q=q->next;
			a[i*m+j]=q->data;
			i++;
			p=p->next;
		}


		i=0;		
		
		s_bac=bac;
		bac=0;
	}
//输出基矩阵
	fp=fopen("out.txt","a");
	fprintf(fp,"\n此问题的基矩阵B为：");
	for (i=0;i<m;i++) 
	{
		fprintf(fp,"\n");
		for (j=0;j<m;j++)
			fprintf(fp,"%.3f ",a[i*m+j]);
	}
	fclose(fp);

//求B逆
	float *b;
	b=new float[m*m*2];
	
	for (i=0;i<m*m*2;i++)
	{
		if (i%(2*m)<m) b[i]=a[(i/(2*m))*m+(i%(2*m))];
		else if ((i%(2*m))-m==(i/(2*m))) b[i]=1;
		else b[i]=0;
	}

	float t;
	int k;
	for (i=0;i<m;i++)
	{
	
		t=b[i*2*m+i];
		for (j=0;j<2*m;j++)//把第i行按第i列单位化
			b[i*2*m+j]/=t;


		for (j=0;j<m;j++)
			if (j!=i)
			{
				t=b[j*2*m+i]/b[i*2*m+i];
				for (k=0;k<2*m;k++)
					b[j*2*m+k]-=b[i*2*m+k]*t;
			}
		
	}

	for (i=0;i<m;i++)
		for (j=0;j<m;j++)
			a[i*m+j]=b[i*2*m+j+m];

//输出B逆

		
	fp=fopen("out.txt","a");
	fprintf(fp,"\n所以B逆为：");
	for (i=0;i<m;i++) 
	{
		fprintf(fp,"\n");
		for (j=0;j<m;j++)
			fprintf(fp,"%.3f ",a[i*m+j]);
	}

	fclose(fp);
	
//求b的灵敏度
//读取x的解
	float *x;
	x=new float[m];

	bac=0;
	s_bac=0;
	for (i=0;i<m;i++)
	{	

		p=chart->next;
		while(p!=NULL)
		{
			if (p->bac>s_bac&&(bac>p->bac||bac==0)) {bac=p->bac;q=p->u;}
			p=p->next;
		}
		s_bac=bac;
		bac=0;
		
		while(q->next!=NULL) q=q->next;
		x[i]=q->data;

	}

	fp=fopen("out.txt","a");
	fprintf(fp,"\n所以资源系数b的灵敏度分析结果为：");
	fclose(fp);

	float bmin=0,bmax=0;

	for (j=0;j<m;j++)
	{
		for (i=0;i<m;i++)
		{
			if (a[i*m+j]>0&&(bmax==0||bmax<-1*x[i]/a[i*m+j])) bmax=-1*x[i]/a[i*m+j];
			if (a[i*m+j]<0&&(bmin==0||bmin>-1*x[i]/a[i*m+j])) bmin=-1*x[i]/a[i*m+j];
		}

		fp=fopen("out.txt","a");
		fprintf(fp,"\n%.3f<=b%d<=%.3f",bmax,j+1,bmin);
		fclose(fp);

		bmax=0;
		bmin=0;
	}


}

void main()
{
	input_data();

    if (t) initial();

	if (t) change();

	if (t) output_data();

	if (t) analysis_c();

	if (t) analysis_b();
}