Demo entry 6660936

4

   

Submitted by 4 on Nov 18, 2017 at 03:57
Language: C. Code size: 6.8 kB.

#include<stdio.h>
#include<math.h>
float f(float x1,float x2);				//惩罚函数//
float g(float k1,float k2,float r);
void powell(float r0);					//鲍威尔法子函数//
float golden(float a1,float a4,float r0,float k1,float k2,float q1,float q2);	//0.618法子函数//
float search(float k1,float k2,float q1,float q2,float r0);  	//一维搜索法子函数//

float v1(float x1,float x2,float i);
float v2(float x1,float x2,float i);
float m(float x1,float x2,float i);

float g1(float x1,float x2);    //约束函数如下//
float g2(float x1,float x2);
float g3(float x1,float x2);
float g4(float x1,float x2);
float g5(float x1,float x2);
float g6(float x1,float x2);
float g7(float x1,float x2);

float x0[2]={4,4};
void main(void)              //以惩罚函数法为主程序//
{
	float r0,c,e1,e2,sb;
		printf("请输入惩罚因子r0=");
		scanf("%f",&r0);
		printf("请输入缩减系数c=");
		scanf("%f",&c);
		printf("请输入收敛精度(e1和e2(逗号分隔)):");
		scanf("%f,%f",&e1,&e2);
	float x00[2],y,l,sum;
	int i,k;
	sum=0.0;
	for(i=0;i<2;i++)
	x00[i]=x0[i];

	for(i=0;;i++)
	{
		powell(r0);  //调用了鲍威尔法子函数//
		sum=0.0;
		for(k=0;k<2;k++);
		sum=sum+(x00[k]-x0[k])*(x00[k]-x0[k]);
		l=sqrt(sum);
		sb=fabs((g(x0[0],x0[1],r0)-g(x00[0],x00[1],r0))/g(x00[0],x00[1],r0));
		printf("sb=%f,l=%f\n",sb,l);
		if(fabs(l)<e1&&(f(x0[0],x0[1])-f(x00[0],x00[1]))/f(x00[0],x00[1])<e2)
		break;
		for(k=0;k<2;k++);
		x00[k]=x0[k];
		r0=c*r0;
	};

	y=f(x0[0],x0[1]);
	printf("由惩罚函数求得最优解如下:\n");
	printf("两个设计变量:x0[0]=%f,x0[1]=%f\n",x0[0],x0[1]);
	printf("目标函数极小值:%f\n", y);
}

void powell(float r0)            //鲍威尔法子函数//
{
	float e[2],x1[4][2],y[4],del[2],d[3][2],a,mdel,l,sum,s,h;  
	int m,j,i,k,z;
		s=0.00001;
		mdel=0.0;
		sum=0.0;
	for(k=0;k<2;k++)
	{
		for(z=0;z<2;z++)
		e[z]=0.0;
		e[k]=1;                                             
		for(i=0;i<2;i++)
		d[k][i]=e[i];
	}                                       
	for(j=0;;j++)
	{
		printf("x[0]=%f,x0[1]=%f\n",x0[0],x0[1]);                                             
		for(k=0;k<2;k++)
		x1[0][k]=x0[k];
	for(k=0;k<2;k++)
	{                                                  
	   a=search(x1[k][0],x1[k][1],d[k][0],d[k][1],r0);  //调用一维搜索法子函数找最优步长//  
		for(i=0;i<2;i++)
		x1[k+1][i]=x1[k][i]+a*d[k][i];
	};                                                        
	for(k=0;k<2;k++)
		d[2][k]=x1[2][k]-x1[0][k];
	for(k=0;k<2;k++)
		x1[3][k]=x1[2][k]+d[2][k];
	for(k=0;k<=3;k++)
		y[k]=g(x1[k][0],x1[k][1],r0);
	for(k=0;k<2;k++)
	{
		del[k]=y[k]-y[k+1];                                  
	if(del[k]>mdel)
	{
		mdel=del[k];                                        
		m=k;
	};                                                         
}                                                          
if(y[3]<y[0]&&(y[0]-2*y[2]+y[3])*(y[0]-y[2]-mdel)*(y[0]-y[2]-mdel)<0.5*mdel*(y[0]-y[3])*(y[0]-y[3]))
{
	a=search(x1[2][0],x1[2][1],d[2][0],d[2][1],r0);             
	for(k=0;k<2;k++)
		x0[k]=x1[2][k]+a*d[2][k];
	for(k=m;k<2;k++)
	{
		for(i=0;i<2;i++)                                      
		d[m][i]=d[m+1][i];    
	};                                                   
}                                                         
else
{
	if(y[3]>y[2])                                          
	{
		for(k=0;k<2;k++)                                        
		x0[k]=x1[2][k];
	}                                          
	else
	{
		for(k=0;k<2;k++)                                      
		x0[k]=x1[3][k];
	};                                       
};                                                    
for(k=0;k<2;k++);
	sum=sum+(x1[2][k]-x1[0][k])*(x1[2][k]-x1[0][k]);
	l=sqrt(sum);
	printf("l=%f",l);
if(fabs(l)<s||(g(x1[2][0],x1[2][1],r0)-g(x1[0][0],x1[0][1],r0))/g(x1[0][0],x1[0][1],r0)<s)
	break;
	sum=0.0;
};
	h=f(x0[0],x0[1]);
	printf("%f\n",h);                                               
}                                                           
float search(float k1,float k2,float q1,float q2,float r0)   //一维搜索法子函数//
{
	float a1,a2,a3,h,c1,c2,c3,a,h0,x1,x2,x3;
        a=0;
        h0=0.0001;
	a1=a;h=h0;
        x1=k1+a1*q1;
        x2=k2+a1*q2;
        c1=g(x1,x2,r0);
	a2=a1+h;
        x1=k1+a2*q1;
        x2=k2+a2*q2;
	c2=g(x1,x2,r0);
	if(c2>c1)
	{
		h=0.0-h;
		a3=a1;
		c3=c1;
		a1=a2;
		a2=a3;
		c1=c2;
		c2=c3;
		a3=a2+h;
		x1=k1+a3*q1;
		x2=k2+a3*q2;
		c3=g(x1,x2,r0);
	}
	else
	{
		a3=a2+h;
        x1=k1+a3*q1;
        x2=k2+a3*q2;
		c3=g(x1,x2,r0);
	}
	while(c3<c2)
	{
		a1=a2;
		a2=a3;
		c1=c2;
		c2=c3;
		a3=a2+h;
		x1=k1+a3*q1;
		x2=k2+a3*q2;
		c3=g(x1,x2,r0);}
		a=golden(a1,a3,r0,k1,k2,q1,q2);  //一维搜索法子函数中调用0.618法搜索最优步长//
			return a;
	}
float golden(float a1,float a4,float r0,float k1,float k2,float q1,float q2)   //0.618法子函数//
{
   float a2,a3,y2,y3,i,e,k,x1,x2;
		   k=0.618;e=0.0001;
		   a2=a4-k*(a4-a1);
		   x1=k1+a2*q1;
		   x2=k2+a2*q2;
		   y2=g(x1,x2,r0);
		   a3=a1+k*(a4-a1);
		   x1=k1+a3*q1;
		   x2=k2+a3*q2;
		   y3=g(x1,x2,r0);
   for(i=0;i<=10000;i++)
     {
        if(y2>=y3)
        { 
			a1=a2;a2=a3;y2=y3;
			a3=a1+k*(a4-a1);
			x1=k1+a3*q1;
		    x2=k2+a3*q2;
			y3=g(x1,x2,r0);
		}
        else
        { 
			a4=a3;a3=a2;y3=y2;
            a2=a4-k*(a4-a1);
            x1=k1+a2*q1;
            x2=k2+a2*q2;
            y2=g(x1,x2,r0);
		}
        if(fabs((a3-a2)/a3)<e&&fabs((y3-y2)/y3)<e)
			break;
        else
			continue;
      }
	return (a2+a3)/2;
}
float f(float x1,float x2)    //目标函数//
{
	float y,sum,i;
	sum=0.0;
	for(i=3.14/60;i<=3.14/2;i=i+3.14/60)
		sum=sum+(v1(x1,x2,i)-v2(x1,x2,i))*(v1(x1,x2,i)-v2(x1,x2,i))*(3.14/60);
	return sum;
}
float v1(float x1,float x2,float i)
{
	float y,r;
	r=m(x1,x2,i);
	y=3.14-acos((r*r+x2*x2-x1*x1)/(2*r*x2))-acos((r*r+24)/(10*r));
	return y;
}
float v2(float x1,float x2,float i)
{
	float y;
	y=acos(((1+x1)*(1+x1)-x2*x2-25)/(10*x2))+(2/(3*3.14))*i*i;
	return y;
}
float m(float x1,float x2,float i)
{
	float y;
	y=sqrt(26-10*cos(i+acos(((1+x1)*(1+x1)-x2*x2+25)/(10*(1+x1)))));
	return y;
}
float g(float k1,float k2,float r)
{
	float y;
	y=f(k1,k2)-r*(1/g1(k1,k2)+1/g2(k1,k2)+1/g3(k1,k2)+1/g4(k1,k2)+1/g5(k1,k2)+1/g6(k1,k2)+1/g7(k1,k2));
	return y;
}

float g1(float x1,float x2)
{
	float y;
	y=-x1;
	return y;
}
float g2(float x1,float x2)
{
	float y;
	y=-x2;
	return y;
}
float g3(float x1,float x2)
{
	float y;
	y=6-x1-x2;
	return y;
}
float g4(float x1,float x2)
{
	float y;
	y=x1-x2-4;
	return y;
}
float g5(float x1,float x2)
{
	float y;
	y=x2-x1-4;
	return y;
}
float g6(float x1,float x2)
{	
	float y;
	y=x1*x1+x2*x2-1.414*x1*x2-16;
	return y;
}
float g7(float x1,float x2)
{	
	float y;
	y=36-x1*x1-x2*x2-1.414*x1*x2;
	return y;
}

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).