optim.cpp

空调控制系统最优算法

// optim.cpp: implementation of the optim class.
//
//////////////////////////////////////////////////////////////////////

#include "optim.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
void Ahu_optim::Process_optim(float C[3],float d[3],float e[3],double Taiw,double Tchws,double Tchwo,double Tchwi,double dT,double Tfresh,double Taout,double Tao,double Troom,double Dfresh,double Daout,double Dao,double Droom,double arph,double Mwmin,double Mwmax,double Mamin,double Mamax,double x[2])
{
	int m=10;//惩罚因子
	double ee=0.001;//精确度
	int c=2;//m=m*c
	double t;//t1,t2比较后的最大值
	double f;//最终优化的函数值
	int g=0;//判断不等式约束的标志位

	double ff;//函数值
	//double x[2];//函数自变量,可作为函数process的返回值
	double p[]={1,2};

	double Ma;//空气质量流量
	double Mw;//冷冻水质量流量

	for(int k=0;;k++)
	{
		printf("第%d次迭代\n",k+1);
		g=Com_g(x,Mwmin,Mwmax,Mamin,Mamax);	
		ff=Powell(p,0.3,0.001,0.0001,2,x,m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
		printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff);
		t=Com_t(x,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
		printf("t=%f\n",t);
		if(t<ee)
			break;
		else
		{
			m=m*c;
			Ma=x[0];
			Mw=x[1];
			f=Ppump_nom*(d[0]+d[1]*(Mw/Mw_nom)+d[2]*(Mw/Mw_nom)*(Mw/Mw_nom)+d[3]*(Mw/Mw_nom)*(Mw/Mw_nom)*(Mw/Mw_nom))+Pfan_nom*(e[0]+e[1]*(Ma/Ma_nom)+e[2]*(Mw/Mw_nom)*(Mw/Mw_nom)+e[3]*(Mw/Mw_nom)*(Mw/Mw_nom)*(Mw/Mw_nom));
			printf("f=%f\n",f);
			printf("\n");
		}
	}
	getch();
}

/*构造目标函数*/
double Ahu_optim::Objf(double x[],int m,int g,float C[3],float d[3],float e[3],double Taiw,double Tchws,double Tchwo,double Tchwi,double dT,double Tfresh,double Taout,double Tao,double Troom,double Dfresh,double Daout,double Dao,double Droom,double arph,double Mwmin,double Mwmax,double Mamin,double Mamax)//目标函数
{
	double ff,dd;
	int b;
	int flag[4];//标志位，用于判断惩罚时用到哪个不等式约束

	double Ma,Mw;
	Ma=x[0];//空气质量流量
	Mw=x[1];//冷冻水质量流量
	
	double Hfresh,Haout,Hao,Hroom;//计算空气焓值
	Hfresh=Air_h(Dfresh,Tfresh);
	Haout=Air_h(Daout,Taout);
	Hao=Air_h(Dao,Tao);
	Hroom=Air_h(Droom,Troom);

	dd=(C[0]*pow(Ma,C[2])/(1+C[1]*pow(Ma/Mw,C[2]))*(Taiw-Tchws)-Mw*(Tchwo-(Tchwi+dT)))*(C[0]*pow(Ma,C[2])/(1+C[1]*pow(Ma/Mw,C[2]))*(Taiw-Tchws)-Mw*(Tchwo-(Tchwi+dT)))+(Ma*(arph*Hfresh+(1-arph)*Haout-Hao)-Ma*(Hroom-Hao))*(Ma*(arph*Hfresh+(1-arph)*Haout-Hao)-Ma*(Hroom-Hao));
	ff=Ppump_nom*(d[0]+d[1]*(Mw/Mw_nom)+d[2]*(Mw/Mw_nom)*(Mw/Mw_nom)+d[3]*(Mw/Mw_nom)*(Mw/Mw_nom)*(Mw/Mw_nom))+Pfan_nom*(e[0]+e[1]*(Ma/Ma_nom)+e[2]*(Mw/Mw_nom)*(Mw/Mw_nom)+e[3]*(Mw/Mw_nom)*(Mw/Mw_nom)*(Mw/Mw_nom))+m*dd;
	
	/*从com_g()函数得到的flag标志，通过和1异或后得到选用了哪个约束*/
	for(int i=0;i<4;i++)
	{
		b=g^1;
		if(g==g-1)
			flag[i]=1;//第i个不等式被用到
		else
			flag[i]=0;//第i个不等式没用到
		g=g>>1;
	}
	ff=ff+m*(flag[0]*(Mwmin-Mw)*(Mwmin-Mw)+flag[1]*(Mw-Mwmax)*(Mw-Mwmax)+flag[2]*(Mamin-Ma)*(Mamin-Ma)+flag[3]*(Ma-Mamax)*(Ma-Mamax));
	return(ff);
}



/*不等式约束条件函数型如：g=a1+a2*x[0]+a3*x[1]*/
/*采用二进制判断*/
int Ahu_optim::Com_g(double x[],double Mwmin,double Mwmax,double Mamin,double Mamax)
{
	//x[0]空气质量流量Ma
	//x[1]冷冻水质量流量Mw
	double g=0;
	int flag=0;
	for(int i=0;i<4;i++)
	{
		switch(i)
		{
			case 0:g=Mamin-x[0];
				if(g<=0)
					g=0;
				else
					flag=flag+1;//无约束优化包含不等式约束1
			case 1:g=-Mamax+x[0];
				if(g<=0)
					g=0;
				else
					flag=flag+2;//无约束优化包含不等式约束2
			case 2:g=Mwmin-x[1];
				if(g<=0)
					g=0;
				else
					flag=flag+4;//无约束优化包含不等式约束3
			case 3:g=-Mwmax+x[1];
				if(g<=0)
					g=0;
				else
					flag=flag+8;//无约束优化包含不等式约束4
		}
	}
	return flag;
}


/*计算t1,t2,max{t1,t2}*/
double Ahu_optim::Com_t(double x[],float C[3],float d[3],float e[3],double Taiw,double Tchws,double Tchwo,double Tchwi,double dT,double Tfresh,double Taout,double Tao,double Troom,double Dfresh,double Daout,double Dao,double Droom,double arph,double Mwmin,double Mwmax,double Mamin,double Mamax)
{
	double t1,t2,maxt1,maxt2;
	double temp;
	int i;

	double Ma=x[0];//空气质量流量
	double Mw=x[1];//冷冻水质量流量

	double Hfresh,Haout,Hao,Hroom;//计算空气焓值
	Hfresh=Air_h(Dfresh,Tfresh);
	Haout=Air_h(Daout,Taout);
	Hao=Air_h(Dao,Tao);
	Hroom=Air_h(Droom,Troom);

	for(i=0;i<2;i++)
	{
		switch(i)
		{
		case 0:t1=C[0]*pow(Ma,C[2])/(1+C[1]*pow(Ma/Mw,C[2]))*(Taiw-Tchws)-Mw*(Tchwo-(Tchwi+dT));
			if(t1<0)
				t1=(-1)*t1;				
			maxt1=t1;
			break;
		case 1:t1=Ma*(arph*Hfresh+(1-arph)*Haout-Hao)-Ma*(Hroom-Hao);
			if(t1<0)
				t1=(-1)*t1;	
			if(t1>maxt1)
				maxt1=t1;
			break;
		}
	}
	for(i=0;i<4;i++)
	{
		switch(i)
		{
		case 0:t2=Mamin-x[0];
			maxt2=t2;
			break;
		case 1:t2=-Mamax+x[0];
			if(t2>maxt2)
				maxt2=t2;
			break;
		case 2:t2=Mwmin-x[1];
			if(t2>maxt2)
				maxt2=t2;
			break;
		case 3:t2=-Mwmax+x[1];
			if(t2>maxt2)
				maxt2=t2;
			break;
		}
	}
	if(maxt1>=maxt2)
		temp=maxt1;
	else
		temp=maxt2;	
	return temp;
}


/*powell算法*/
double Ahu_optim::Gold(double a[],double b[],double eps,int n,double xx[],int m,int g,float C[3],float d[3],float e[3],double Taiw,double Tchws,double Tchwo,double Tchwi,double dT,double Tfresh,double Taout,double Tao,double Troom,double Dfresh,double Daout,double Dao,double Droom,double arph,double Mwmin,double Mwmax,double Mamin,double Mamax)//一维搜索0.618法
{
	int i;
	double f1,f2,*x[2],ff,q,w;
	for(i=0;i<2;i++)
		x[i]=(double *)malloc(n*sizeof(double));
	for(i=0;i<n;i++)
	{
		*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
		*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
	}
	f1=Objf(x[0],m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
	f2=Objf(x[1],m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
	do
	{
		if(f1>f2)
		{
			for(i=0;i<n;i++)
			{
				b[i]=*(x[0]+i);
				*(x[0]+i)=*(x[1]+i);
			}
			f1=f2;
			for(i=0;i<n;i++)
				*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
			f2=Objf(x[1],m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
		}
		else
		{
			for(i=0;i<n;i++)
			{
				a[i]=*(x[1]+i);
				*(x[1]+i)=*(x[0]+i);}
				f2=f1;
				for(i=0;i<n;i++)
					*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
				f1=Objf(x[0],m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
			}
			q=0;
			for(i=0;i<n;i++)
				q=q+(b[i]-a[i])*(b[i]-a[i]);
			w=sqrt(q);
	}while(w>eps);
	for(i=0;i<n;i++)
		xx[i]=0.5*(a[i]+b[i]);
	ff=Objf(xx,m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
	for(i=0;i<2;i++)
		free(x[i]);
	return(ff);
}


void Ahu_optim::Jtf(double x0[],double h0,double s[],int n,double a[],double b[],int m,int g,float C[3],float d[3],float e[3],double Taiw,double Tchws,double Tchwo,double Tchwi,double dT,double Tfresh,double Taout,double Tao,double Troom,double Dfresh,double Daout,double Dao,double Droom,double arph,double Mwmin,double Mwmax,double Mamin,double Mamax)
{
	int i;
	double *x[3],h,f1,f2,f3;
	for(i=0;i<3;i++)
		x[i]=(double *)malloc(n*sizeof(double));
	h=h0;
	for(i=0;i<n;i++)
		*(x[0]+i)=x0[i];
	f1=Objf(x[0],m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
	for(i=0;i<n;i++)
		*(x[1]+i)=*(x[0]+i)+h*s[i];
	f2=Objf(x[1],m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
	if(f2>=f1)
	{
		h=-h0;
		for(i=0;i<n;i++)
			*(x[2]+i)=*(x[0]+i);
		f3=f1;
		for(i=0;i<n;i++)
		{
			*(x[0]+i)=*(x[1]+i);
			*(x[1]+i)=*(x[2]+i);
		}
		f1=f2;
		f2=f3;
	}
	for(;;)
	{
		h=2*h;
		for(i=0;i<n;i++)
			*(x[2]+i)=*(x[1]+i)+h*s[i];
		f3=Objf(x[2],m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
		if(f2<f3) 
			break;
		else
		{
			for(i=0;i<n;i++)
			{
				*(x[0]+i)=*(x[1]+i);
				*(x[1]+i)=*(x[2]+i);
			}
			f1=f2;
			f2=f3;
		}
	}
	if(h<0)
		for(i=0;i<n;i++)
		{
			a[i]=*(x[2]+i);
			b[i]=*(x[0]+i);
		}
	else
		for(i=0;i<n;i++)
		{
			a[i]=*(x[0]+i);
			b[i]=*(x[2]+i);
		}
		for(i=0;i<3;i++)
			free(x[i]);
}


double Ahu_optim::Oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[],int m,int g,float C[3],float d[3],float e[3],double Taiw,double Tchws,double Tchwo,double Tchwi,double dT,double Tfresh,double Taout,double Tao,double Troom,double Dfresh,double Daout,double Dao,double Droom,double arph,double Mwmin,double Mwmax,double Mamin,double Mamax)
{
	double *a,*b,ff;
	a=(double *)malloc(n*sizeof(double));
	b=(double *)malloc(n*sizeof(double));
	Jtf(x0,h0,s,n,a,b,m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
	ff=Gold(a,b,epsg,n,x,m,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
	free(a);
	free(b);
	return (ff);
}

double Ahu_optim::Powell(double p[],double h0,double eps,double epsg,int n,double x[],int mm,int g,float C[3],float d[3],float e[3],double Taiw,double Tchws,double Tchwo,double Tchwi,double dT,double Tfresh,double Taout,double Tao,double Troom,double Dfresh,double Daout,double Dao,double Droom,double arph,double Mwmin,double Mwmax,double Mamin,double Mamax)
{
	int i,j,m;
	double *xx[4],*ss,*s;
	double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,dd;
	ss=(double *)malloc(n*(n+1)*sizeof(double));
	s=(double *)malloc(n*sizeof(double));
	for(i=0;i<n;i++)
	{
		for(j=0;j<=n;j++)
			*(ss+i*(n+1)+j)=0;
		*(ss+i*(n+1)+i)=1;
	}
	for(i=0;i<4;i++)
	xx[i]=(double *)malloc(n*sizeof(double));
	for(i=0;i<n;i++)
		*(xx[0]+i)=p[i];
	for(;;)
	{
		for(i=0;i<n;i++)
		{
			*(xx[1]+i)=*(xx[0]+i);
			x[i]=*(xx[1]+i);
		}
		f0=f1=Objf(x,mm,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
		dlt=-1;
		for(j=0;j<n;j++)
		{
			for(i=0;i<n;i++)
			{
				*(xx[0]+i)=x[i];
				*(s+i)=*(ss+i*(n+1)+j);
			}
			f=Oneoptim(xx[0],s,h0,epsg,n,x,mm,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
			df=f0-f;
			if(df>dlt)
			{
				dlt=df;
				m=j;
			}
		}
		sdx=0;
		for(i=0;i<n;i++)
			sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
		if(sdx<eps)
		{
			free(ss);
			free(s);
			for(i=0;i<4;i++)
				free(xx[i]);
			return(f);
		}
		for(i=0;i<n;i++)
			*(xx[2]+i)=x[i];
		f2=f;
		for(i=0;i<n;i++)
		{
			*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
			x[i]=*(xx[3]+i);
		}
		fx=Objf(x,mm,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
		f3=fx;
		q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
		dd=0.5*dlt*(f1-f3)*(f1-f3);
		if((f3<f1)||(q<dd))
		{
			if(f2<=f3)
				for(i=0;i<n;i++)
					*(xx[0]+i)=*(xx[2]+i);
			else
				for(i=0;i<n;i++)
					*(xx[0]+i)=*(xx[3]+i);
		}
		else
		{
			for(i=0;i<n;i++)
			{
				*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
				*(s+i)=*(ss+(i+1)*(n+1));
			}
			f=Oneoptim(xx[0],s,h0,epsg,n,x,mm,g,C,d,e,Taiw,Tchws,Tchwo,Tchwi,dT,Tfresh,Taout,Tao,Troom,Dfresh,Daout,Dao,Droom,arph,Mwmin,Mwmax,Mamin,Mamax);
			for(i=0;i<n;i++)
				*(xx[0]+i)=x[i];
			for(j=m+1;j<=n;j++)
				for(i=0;i<n;i++)
					*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
		}
	}
}


//d_含湿量
//t_摄氏温标
double Ahu_optim::Air_h(double d,double t)
{
	double h;
	h=(1.01+1.84*d)*t+2500*d;
	return h;
}