qp method.txt

测试无约束优化定理求解最小值,包括测试数据,在VC下可以实现,已经调试成功.

	//第一步：确定初始点，求其有效集
	I=new int[100*m];
	for(k=0;k<100;k++)
		for(j=0;j<m;j++)           //全部给初值为-1,表示有效集为空
		{
           I[m*k+j]=-1;
		}
	k=0;
	for(j=0;j<m;j++)                    //求X0的有效集
	{        
		    c=0;                        //c要初始化，否则出错
			//c=A[j]*X0[j]-B[j];
		    for(i=0;i<n;i++)
			{
				c+=A[i*m+j]*X0[i];
			}
			c=c-B[j];
			if(c>=-0.00000001&&c<=0.00000001)
			{
				I[j]=j;                //实际上是m*k+j,这里是第一次,所以省略
				set[k]++;
			}
	}
    p=multiply(G,X0,n,n,1);    //X0为初始点,只有一列
    gk=add(p,r,n); 
    delete []p; 
    d=new double[n]; 
	nami=new double[m]; 
	/*	for(j=0;j<m;j++)
		{
            cout<<I[j]<<" ";
		}*/
	
	while(flag!=1)
	{
		cout<<"第"<<k+1<<"次有效约束个数为: ";
		cout<<set[k]<<endl;
		if(set[k]==0)                          //无约束条件
		{ 
			    cout<<"解无约束二次规划:"<<endl;
				temp=tranverse(G,n); 
				p=multiply(temp,gk,n,n,1); 
				delete []temp; 
				for(i=0;i<n;i++) 
					d[i]=-p[i]; 
				delete []p; 
				for(i=0;i<m;i++) 
					nami[i]=0;
		} 
		else
		{ 
			ak=new double[n*set[k]]; 
			total=0;                      //控制ak的列
			for(i=0;i<m;i++) 
			{
				//cout<<I[m*k+i]<<" ";
				if(I[m*k+i]!=-1)
				{ 
					for(j=0;j<n;j++) 
					     ak[set[k]*j+total]=A[j*m+I[m*k+i]]; //第几个约束条件就代表取第几列放入ak
				    total++;
				}
			}
			/*for(i=0;i<n;i++)
			{
		        for(j=0;j<k;j++)
				{
                       cout<<ak[i*set[k]+j]<<" ";
				}
			}
            char ch=getchar();*/
				bk=new double[total]; 
			for(i=0;i<total;i++) 
				bk[i]=0;
			j=total;
			total=n+total; 
			matrixL=new double[total*total]; 
			matrixR=new double[total];
		    matrixL=ConstructMatrixL(G,ak,n,set[k]);
			for(i=0;i<total;i++)
			{	
				for(j=0;j<total;j++)
					cout<<matrixL[i*total+j]<<" ";
				cout<<endl;
			}
			cout<<endl;
            matrixR=ConstructMatrixR(gk,bk,n,set[k]);
			for(i=0;i<total;i++)
			{	
					cout<<matrixR[i]<<endl;
			}
            delete []ak; 
			delete []bk; 
			p=root(total,matrixL,matrixR); 
			delete []matrixL; 
            delete []matrixR;
			for(i=0;i<n;i++) 
				d[i]=p[i]; 
			j=0; 
			for(i=0;i<m;i++)                              //nami数组保存乘子向量。 
			{
				if(I[k*m+i]!=-1)
				{ 
					nami[i]=p[n+j]; 
				    j++;
				} 
			    else 
				nami[i]=0;
			} 
		}
		sign=1;
		for(i=0;i<n;i++) 
		{	
			if(d[i]<-0.00000001||d[i]>0.00000001)        //如果d[i]!=0
			{ 
					sign=0; 
					break;
			} 
		}
		if(sign==1)                                     //sign为1则d等于0.
		{                                                            
			cout<<"d全为0!!"<<endl; 
			flag=1; 
			x=0.0; 
			for(i=0;i<m;i++)                            //求最小的nami
			{
				if(nami[i]<-0.00000001)
				{ 
					flag=0; 
					if(nami[i]<x)
					{ 
						x=nami[i]; 
						y=i;
					}
				}
			}
				if(flag==1)                            //求解成功!
				{ 
					cout<<"求解成功!"<<endl;
					temp=new double[2];
					for(i=0;i<n;i++)
						temp[i]=X0[i];
					cout<<endl;
					return temp;                          //传入的变量不能作为返回值
					delete []d; 
					delete []gk;                          //有分歧
					//break;
				}
				else
				{
					cout<<"有效集个数减一个:";
					I[m*k+y]=-1;
					set[k]--;                        //有效集减一个,下一轮循环开始时给set[k+1]
                    for(i=0;i<m;i++)                                 //X0不变
					{
                       I[m*(k+1)+i]=I[m*k+i];
					}
				}
		} 
		else
		{ 
			a=Find_ak(I,A,B,X0,d,n,m,k); 
			a1=NumMul(d,a,n); 
			a0=add(X0,a1,n); 
			for(i=0;i<n;i++) 
				X0[i]=a0[i]; 
			delete []a0; 
			delete []a1; 
			delete []gk; 
			q=multiply(G,X0,n,n,1); 
			gk=add(q,r,n); 
			delete []q;
		}
		k=k+1;
		set[k]=set[k-1];
	}  
}
/************************************************************************************************/
void main()
{
	int i,j;
	double *p;
    //double q0[4][4]={1,1,2,1,4,5,8,7,6,3,2,4,-1,-2,-3,-4};
	double q0[16]={1,1,2,1,4,5,8,7,6,3,2,4,-1,-2,-3,-4};
	//double q1[4][3]={1,0,0,0,0,1,0,0,0,0,1,0};
	double q1[12]={1,0,0,0,0,1,0,0,0,0,1,0};
	double q2[16]={1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1};
    p=multiply(q0,q1,4,4,3);
	cout<<"乘法的结果是: ";
	for(i=0;i<12;i++)
	cout<<*(p+i)<<" ";
	delete []p;
	p=Div(q0,1,16);
	cout<<endl;
	cout<<"除法的结果是: ";
	for(i=0;i<16;i++)
	cout<<*(p+i)<<" ";
	delete []p;
	p=sub(q0,q2,16);
	cout<<endl;
	cout<<"减法的结果是: ";
	for(i=0;i<16;i++)
	cout<<*(p+i)<<" ";
    delete []p;
	p=add(q0,q2,16);
	cout<<endl;
	cout<<"加法的结果是: ";
	for(i=0;i<16;i++)
	cout<<*(p+i)<<" ";
    delete []p;
	cout<<endl;
/************************************以上测试矩阵运算********************************************/
    double matrix[4]={1,1,1,-1};
	double B[2]={3,2};
    //p=root(matrix,2);
	p=root(2,matrix,B);
	for(i=0;i<2;i++)
	cout<<*(p+i)<<" ";
	delete []p;
	cout<<endl;
/************************************以上测试方程组求解******************************************/
    double G[9]={1,2,3,4,5,6,7,8,9};
    double Ak[6]={1,2,3,4,5,6};
	p=ConstructMatrixL(G,Ak,3,2);
	for(i=0;i<5;i++)
	{	
		for(j=0;j<5;j++)
			cout<<*(p+i*5+j)<<" ";
		cout<<endl;
	}
	cout<<endl;
	delete []p;
	double r[3]={0,0,0};
    double Bk[2]={1,2};
    p=ConstructMatrixR(r,Bk,3,2);
	for(i=0;i<5;i++)
		cout<<p[i]<<" ";
	cout<<endl;
	delete []p;
/***********************************以上测试矩阵的合并*******************************************/
	double *TestG;
	double *TestA;
	double *TestB;
	double *Testr;
	double *matrixA;
	double *matrixB;
	int n,m;                              //n为G的维数，m为A的列数
	double *X;
	cout<<"下面是用定理求解最小值的测试部分: "<<endl;
    cout<<"请您输入二次规划式的二次型正定矩阵G的维数: ";
	cin>>n;
    TestG=new double[n*n];
    //cout<<"请您输入二次规划式的二次型正定矩阵G: ";
	for(i=0;i<n;i++)
	{
		cout<<"请您输入二次规划式的二次型正定矩阵G的第"<<i+1<<"行: ";
		for(j=0;j<n;j++)
		{
			cin>>TestG[i*n+j];
		}
	}
    cout<<"请您输入二次规划式的系数矩阵A的列数(行数为"<<n<<")";
	cin>>m;
    TestA=new double[n*m];
	for(i=0;i<n;i++)
	{
		cout<<"请您输入二次规划式的系数矩阵A的第"<<i+1<<"行: ";
		for(j=0;j<m;j++)
		{
			cin>>TestA[i*m+j];
		}
	}
    Testr=new double[n];
    cout<<"请您输入二次规划式的一次系数矩阵r: ";
    for(i=0;i<n;i++)
	{
		cin>>Testr[i];
	}
	TestB=new double[m];
    cout<<"请您输入二次规划式的系数矩阵B: ";
    for(i=0;i<m;i++)
	{
		cin>>TestB[i];
	}
    matrixA=new double[(n+m)*(n+m)];
    matrixA=ConstructMatrixL(TestG,TestA,n,m);
    matrixB=new double[n+m];
	matrixB=ConstructMatrixR(Testr,TestB,n,m);
    X=new double[n+m];
    X=root(n+m,matrixA,matrixB);
    for(i=0;i<n+m;i++)
	{
		cout<<X[i]<<" ";
	}
	cout<<endl;
    delete []TestG;
    delete []TestA;
	delete []TestB;
	delete []Testr;
	delete []matrixA;
	delete []matrixB;
	delete []X;
	cout<<endl;
/***********************************以上测试定理求解最小值***************************************/
	//double a[4]={1,3,4,7};
	double a[9]={1,3,5,2,4,6,5,6,8};
	p=new double[9];
	p=tranverse(a,3);
    for(i=0;i<3;i++) 
	{
		for(j=0;j<3;j++) 
           cout<<p[i*3+j]<<" "; 
            cout<<endl; 
	} 
	delete []p;
/***********************************以上测试矩阵求逆*********************************************/
    p=new double[2];
	/*double TG[4]={2,0,0,2};
    double TA[6]={1,-1,0,
		          1,0,-1};
    double TB[3]={1,0,0};
	double TX[2]={0,0};
    double Tr[2]={-2,-4};
	p=QP(TG,Tr,TA,TB,TX,2,3);
	double TG[4]={2,0,0,2};
    double TA[10]={ -1,1, 1,-1,0,
		            2,2,-2, 0,-1};
    double TB[5]={2,6,2,0,0};
	double TX[2]={2,0};
    double Tr[2]={-2,-5};*/
	cout<<"以下求解QP算法: "<<endl;
	cout<<"请您输入二次规划式的二次型正定矩阵G的维数: ";
	cin>>n;
	double *TG;
	TG=new double[n*n];
	for(i=0;i<n;i++)
	{
		cout<<"请您输入二次规划式的二次型正定矩阵G的第"<<i+1<<"行: ";
		for(j=0;j<n;j++)
		{
			cin>>TG[i*n+j];
		}
	}
    cout<<"请您输入二次规划式的系数矩阵A的列数(行数为"<<n<<")";
	cin>>m;
    double *TA;
	TA=new double[n*m];
	for(i=0;i<n;i++)
	{
		cout<<"请您输入二次规划式的系数矩阵A的第"<<i+1<<"行: ";
		for(j=0;j<m;j++)
		{
			cin>>TA[i*m+j];
		}
	}
    double *Tr;
	Tr=new double[n];
    cout<<"请您输入二次规划式的一次系数矩阵r: ";
    for(i=0;i<n;i++)
	{
		cin>>Tr[i];
	}
	double *TB;
	TB=new double[m];
    cout<<"请您输入二次规划式的系数矩阵B: ";
    for(i=0;i<m;i++)
	{
		cin>>TB[i];
	}
	double *TX;
	TX=new double[n];
    cout<<"请您输入初始解: ";
    for(i=0;i<n;i++)
	{
		cin>>TX[i];
	}
	cout<<endl;
	cout<<endl;
	p=QP(TG,Tr,TA,TB,TX,n,m);
	for(i=0;i<2;i++)
		cout<<p[i]<<" ";
    delete []p;
/***********************************以上测试QP算法***********************************************/
}