resrict.cpp

约束最优化 最优化试验程序 自己编写 结果经过验证

/*************************************************************/
/************************p2003050221**************************/
/*************************************************************/
/*************************************************************/
#include<iostream.h>
double DOIT(double m1[2],double m2[2]){
	//DOIT
	double sum=0;
	for(int i=0;i<2;i++)
		sum+=m1[i]*m2[i];
	return sum;
}//DOIT
void calg(double x[2],double g[2]){
	//calg
	g[0]=2*x[0]-2;
	g[1]=2*x[1]-4;
}//calg
void QP(int I[3],double G[2][2],double x[2],double A[3][2],double d[2],double namta[3]){
	//QP问题的求解
	//LU=========方法
	int i,j,k,q;
    int count=0;
	double g[2];
	double tempx[5];
	double L[6][6];
	for(i=0;i<3;i++)
		if(I[i])count++;
	if(I[0]&&I[2]){
		for(i=0;i<2;i++)
			A[1][i]=A[2][i];
	}
	calg(x,g);
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
		L[i][j]=G[i][j];
	for(i=0;i<2;i++)
		for(j=2;j<2+count;j++)
			L[i][j]=-A[j-2][i];
	for(i=2;i<2+count;i++)
		for(j=0;j<2;j++)
			L[i][j]=A[i-2][j];
	for(i=2;i<count+2;i++)
		for(j=2;j<=count+2;j++)
			L[i][j]=0;
	for(i=0;i<2;i++)
		L[i][count+2]=-g[i];
	cout<<"求解新的QP问题增广矩阵如下======";
		for(i=0;i<2+count;i++){
		cout<<'\n';
		for(j=0;j<3+count;j++)
			cout<<L[i][j]<<'\t';}
		cout<<'\n';
		/*==lu分解过程=*/
	double temp1;
	for(k=0;k<2+count;k++)
	{
		for(j=k;j<=2+count;j++){
			temp1=0;
		for(q=0;q<k-1;q++)
			temp1+=L[k][q]*L[q][j];
		L[k][j]=L[k][j]-temp1;
		}
		for(i=k+1;i<2+count;i++){
			temp1=0;
		for(q=0;q<k-1;q++)
			temp1+=L[k][q]*L[q][k];
		L[i][k]=(L[i][k]-temp1)/L[k][k];
		}
	}
//	for(i=0;i<2+count;i++){
//		cout<<'\n';
//		for(j=0;j<3+count;j++)
//			cout<<L[i][j]<<'\t';}
	//====回带过程=======
	tempx[count+1]=L[1+count][2+count]/L[1+count][1+count];
	double sum=0;
	for(i=count;i>=0;i--){
		sum=0;
		for(j=i+1;j<count+2;j++)
			sum+=L[i][j]*tempx[j];
		tempx[i]=(L[i][2+count]-sum)/L[i][i];
	}
	cout<<"====*****************解如下:";
	for(i=0;i<count+2;i++)
		cout<<tempx[i]<<'\t';
	for(i=0;i<2;i++)
	   d[i]=tempx[i];
	for(j=2;j<2+count;j++)
		namta[j-2]=tempx[j];
	cout<<'\n';
//	for(i=0;i<2;i++)
//		cout<<d[i]<<'\t';
//	for(i=0;i<2;i++)
//		cout<<namta[i]<<'\t';
}
int if_d(double d[2]){
	for(int i=0;i<2;i++)
		if(d[i])return 1;
		return 0;
}
double min(double namta[3],int &temp_I){

	if(namta[0]<=namta[1]){temp_I=0;return namta[0];}
	else{temp_I=1;return namta[1];}
	
}
void compute_alpha(int I[3],double A[3][2],double d[2],double b[3],double x[2],double &alpha){
	//compute_alpha
	double temp1,temp2;
	int k;
	double temp[3];
	for(int i=0;i<3;i++)
        if(I[i]==0){
			temp2=DOIT(A[i],d);
			if(temp2<0){
				temp1=DOIT(A[i],x);
				temp[i]=(b[i]-temp1)/temp2;
				k=i;
			}
		}
		if(temp[k]<1){alpha=temp[k];I[k]=1;}
		else alpha=1;
}

void main()
{
	/*step1*/
	int i;
	int I[3]={1,1,0};//有效集,1为是,0为否
	int k=1;
	double x[2]={0,0};//initial point
	double G[2][2]={2,0,0,2};//HESSEN
	double A[3][2]={1,0,0,1,-1,-1};//s.t.
	double b[3]={0,0,-1};
	double d[2],namta[3];
	double namta1;
	int temp_I;
	double alpha;
step2:;	/*step2*/
	QP(I,G,x,A,d,namta);
//		for(i=0;i<2;i++)
//		cout<<d[i]<<'\t';
//	for(i=0;i<2;i++)
//		cout<<namta[i]<<'\t';
	if(if_d(d))goto step4;
	else goto step3;
step3:  
	namta1=min(namta,temp_I);
	cout<<"namata为"<<namta1<<'\n';
	if(namta1>=0)goto end;
	else I[temp_I]=0;
	goto step2;
step4: 
	compute_alpha(I,A,d,b,x,alpha);
//	for(i=0;i<3;i++)
//		cout<<I[i];
	cout<<"alpha为"<<alpha<<'\n';
	for(i=0;i<2;i++)
		x[i]=x[i]+alpha*d[i];
//step6:
	k++;goto step2;
end:
	cout<<"x0:"<<x[0]<<'\t'<<"x1:"<<x[1]<<endl;
	cout<<"k:"<<k;
}