gongetidu.c

基于牛顿迭代法的二次规划程序

#define A 2
#define B 3
#define C 4
#define D 0
#include <math.h>

float x0[4],x1[4],g[3],d[4],dif[4],r;
int time[100],i,n;

/*目标函数*/
float ff(float x [4])
{ 
	float ff;
    ff = x[0]*x[0]+x[1]*x[1]+2*x[2]*x[2]+x[3]*x[3]-5*x[0]-5*x[1]-21*x[2]+7*x[3]+D ;
	return ff;
}

/*约束条件*/
float *pp(float x[4])
{
	float *p3;
	g[0] = -x[0]*x[0]-x[1]*x[1]-x[2]*x[2]-x[3]* x[3]-x[0]+ x[1]-x[2]+x[3]+A ;
	g[1] = -x[0]*x[0]-2* x[1]*x[1]- x[2]*x[2]-2* x[3]* x[3]+ x[0]+x[3]+B ;
	g[2] = -2*x[0]*x[0]- x[1]*x[1]- x[2]*x[2]+x[1]+ x[3]+C;
	p3=g;
	return(p3);
}


/*障碍函数*/
float gg(float x[4])
{
	float h[3],f,gg,*p0;
	int i;
	p0=pp(x);	
	for(i=0;i<3;i++)
		h[i]=*(p0+i);		
	f=ff(x);
	gg=f+r*(1/h[0]+1/h[1]+1/h[2]);
	return gg;
}



/*导数*/
float *diff(float x[4])
{
	float *di,*p1;
	int i;
	p1=pp(x);	
	for(i=0;i<3;i++)
		g[i]=*(p1+i);
	dif[0]=2*x[0]-5-r*((-2*x[0]-1)/(g[0]*g[0])+(-2*x[0]+1)/(g[1]*g[1])+(-4*x[0]-2)/(g[2]*g[2]));
	dif[1]=2*x[1]-5-r*((-2*x[1]+1)/(g[0]*g[0])+(-4*x[1])/(g[1]*g[1])+(-2*x[1]+1)/(g[2]*g[2]));
	dif[2]=4*x[2]-21-r*((-2*x[2]-1)/(g[0]*g[0])+(-2*x[2])/(g[1]*g[1])+(-2*x[2])/(g[2]*g[2]));
	dif[3]=2*x[3]+7-r*((-2*x[3]+1)/(g[0]*g[0])+(-4*x[3]+1)/(g[1]*g[1])+1/(g[2]*g[2]));
	di=dif;
	return(di);
}


/*求q*/
float findq(float x0[4],float x1[4])
{
	float h0[4],h1[4],q,*p2,a,b;
	int i;
	a=0;
	b=0;
	p2=diff(x0);
	for(i=0;i<4;i++)
		h0[i]=*(p2+i);
	p2=diff(x1);
	for(i=0;i<4;i++)
		h1[i]=*(p2+i);
	for(i=0;i<4;i++)
		a=a+h1[i]*h1[i];
	for(i=0;i<4;i++)
		b=b+h0[i]*h0[i];
	q=labs(a/b);
	return q;
}


/*主函数*/
main()
{
	float f,a,q,e,*p;
	int m,l;
	r=0.5;
	do
	{
		printf("Please input the initial values:\n");
		for(i=0;i<4;i++)
		{
			printf("x[%d]=",i);
			scanf("%f",&x0[i]);
		}
		p=pp(x0);
		if((*p)<0 || (*(p+1)<0) || (*(p+2)<0))
			printf("The values do not satisfy restriction\n");
	}
	while((*p)<0 || (*(p+1)<0) || (*(p+2)<0));


	p=diff(x0);
	for(i=0;i<4;i++)
		d[i]=-(*(p+i));
	e=d[0]*d[0]+d[1]*d[1]+d[2]*d[2]+d[3]*d[3];


	for(l=0;l<100;l++)
	{
		r=r*0.4;
		for(m=0;m<10000;m++)
		{
			a=0.0001;			
			for(i=0;i<4;i++)
				x1[i]=x0[i]+a*d[i];

			while((gg(x1)-gg(x0))>0.0001)
			{	
					a=a/2;
					for(i=0;i<4;i++)
						x1[i]=x0[i]+a*d[i];
			}

			p=pp(x1);
			if( (*p)<0 )
				break;
			q=findq(x0,x1);
			p=diff(x1);
			for(i=0;i<4;i++)
				d[i]=-(*(p+i))+q*d[i];
			e=d[0]*d[0]+d[1]*d[1]+d[2]*d[2]+d[3]*d[3];
			if(e<0.0001)
				break;
			for(i=0;i<4;i++)
				x0[i]=x1[i];
		}
		time[l]=m+1;	
		if(r*(g[0]+g[1]+g[2])<0.00001)
			break;
	}
	for(i=0;i<100;i++)
	{
		if(time[i]==0)
			break;
		n+=time[i];
	}
	printf("iterative: %d \n",n);
	printf("The values are:\n");
	for(i=0;i<4;i++)
		printf("x[%d]=%f\n",i,x0[i]);
	f=ff(x0);
	printf("The function value is:%f\n",f);
}