二次规划：有效集法（Active-Set）

最新推荐文章于 2024-09-11 07:46:28 发布
最新推荐文章于 2024-09-11 07:46:28 发布
阅读量1.1w 收藏 34
点赞数 3
分类专栏： 最优化算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/zc02051126/article/details/8588644
版权

参考参考了《Practical Optimization Methods With Mathematic Applications》中的8.4节中介绍的有效集法（Active-Set），有效集法只能优化中等规模的最优化问题，如果是大规模的最优化问题，应该采用其它算法，代码如下：（代码执行环境VS2010）

QP.h
#ifndef _QP_
#define _QP_
#pragma comment (lib,"./bin/blasd.lib")
#pragma comment (lib,"./bin/clapackd.lib")
#pragma comment (lib,"./bin/libf2cd.lib")
#pragma comment (lib,"./bin/tmglibd.lib")

#include <fstream>
#include <assert.h>
#include <vector>
#include <algorithm>
#include <iostream>

#include "Lapack.h"
using namespace std;

/*
Copy right 
用拉格朗日乘子发求解等式约束的二次规划问题
求解问题为：min: 0.5*trans(x)*H*x + trans(c)*x; st: A*x = b

参考了《Practical Optimization Methods With Mathematic Applications》中的8.4节中介绍的有效集法（Active-Set）
*/
class QuadEquLag //等式约束问题
{
private:
	//
	double** rTmp;
public:
	double* optimizer(double** H,int& m, int& n, double* c, double** A, double* b);//返回最优解
	void calculateQ(double** reverseH, int& row, int& col, double** A,double** Q);
	void calculateR(double** reverseH, int& row, int& col, double** A,double** R);
	void twoMatrixMultiply(double** a, int& aRow, int& aCol, char* aType, double** b, int& bRow, int& bCol, char* bType,double** out);
	void inverseMatrix(double** a, int& aRow, int& aCol); //计算矩阵a的逆矩阵
	
	void testMatrixInverse();	
};

class QuadProg : public QuadEquLag //不等式约束问题
{
private:

public:
	//double* optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0);
	double optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0); //有效集算法只能处理中等规模的数据
	int rank(double** a, int& row, int& col, double tol=1e-10); //用SVD分解计算矩阵的秩
	void solveLinearEquations(double** a, int& row, int& col, double*);
	bool greaterThanZero(double* in, int& n);
	bool equalZero(double* in, int& n);
	void updateActiveSet(vector<int>& activeSet, double* mu);
	double findAfak(vector<int>& activeSet, double* dk, double** A , double* b,double* xk, int& minIndex, int& m, int& n);
	void updataXk(double* xk,double& afak, double* dk, int& n);
};

#endif

 

#include "QP.h"

void QuadEquLag::testMatrixInverse()
{
	int N = 500;
	double* A = new double[9];
	
	int LDA = 500;
	int* IPIV = new int[N];
	int LWORK = N;
	double* WORK = new double[N];
	int INFO ;
	ifstream readin;
	readin.open("./data/m.txt",ios_base::in);
	int k = 0;
	if(readin.is_open())
	{
		for(int i = 0; i < N; i++)
		{
			for(int j = 0; j < N; j++)
			{
				double tmp = 0;
				readin>>tmp;
				A[k++] = tmp;
			}
		}
	}
	int m = N; int n = m; 
	dgetrf_(&m, &n, A, &m, IPIV, &INFO);
	dgetri_(&N, A, &LDA, IPIV, WORK, &LWORK, &INFO);
}


QP.CPP
double* QuadEquLag::optimizer(double** H,int& m, int& n, double* c, double** A, double* b)
{
	//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
	double* res = new double[n];

	inverseMatrix(H, n, n);
	double** R = new double*[n];
	for(int i = 0; i < n; i++)
	{
		R[i] = new double[m];
	}
	
	calculateR(H, m, n, A,R);
	rTmp = R;

	double** Q = new double*[n];
	for(int i = 0; i < n; i++)
	{
		Q[i] = new double[n];
	}
	calculateQ(H, m, n, A,Q);
	
	double* Qc = new double[n];
	for(int i = 0; i < n; i++)
	{
		double tmp1 = 0;
		double tmp2 = 0;
		for(int k1 = 0; k1 < n; k1++)
		{
			tmp1 = tmp1 + Q[i][k1]*c[k1];
		}

		for(int k2 = 0; k2 < m; k2++)
		{
			tmp2 = tmp2 + R[i][k2]*b[k2];
		}

		res[i] = -tmp1 + tmp2;
	}

	for(int i = 0; i < n; i++)
	{
		delete[] R[i];
	}
	delete[] R;
	for(int i = 0; i < n; i++)
	{
		delete[] Q[i];
	}
	delete[] Q;
	return res;
}
void QuadEquLag::calculateQ(double** reverseH, int& m, int& n, double** A,double** Q)
{
	//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
	================计算H逆乘A转置=====================
	//double** rHtA = new double*[n]; //n*m维
	//for(int i = 0; i < n; i++)
	//{
	//	rHtA[i] = new double[m];
	//}
	//twoMatrixMultiply(reverseH,n,n,"N",A,m,n,"T",rHtA);
	=====================================================
	===============计算A乘H逆乘A转置=====================
	//double** ArH = new double*[m]; //m*n维
	//for(int i = 0; i < m; i++)
	//{
	//	ArH[i] = new double[n];
	//}
	//twoMatrixMultiply(A,m,n,"N",reverseH,m,n,"N",ArH);
	//double** ArHtA = new double*[m];
	//for(int i = 0; i < m; i++)
	//{
	//	ArHtA[i] = new double[m];
	//}
	//twoMatrixMultiply(ArH,m,n,"N",A,m,n,"T",ArHtA);
	//inverseMatrix(ArHtA, m, m);
	释放内存空间
	//for(int i = 0; i < m; i++)
	//{
	//	delete[] ArH[i];
	//}
	//delete[] ArH;
	//double** tmp1 = new double*[n];
	//for(int i = 0; i < n; i++)
	//{
	//	tmp1[i] = new double[m];
	//}
	//calculateR(reverseH, m, n, A,tmp1);
	//=====================================================
	//=====================计算A乘H逆======================
	double** AtH = new double*[m];
	for(int i = 0; i < m; i++)
	{
		AtH[i] = new double[n];
	}
	twoMatrixMultiply(A,m,n,"N", reverseH,n,n,"N",AtH);
	//=====================================================
	释放内存
	//for(int i = 0; i < n; i++)
	//{
	//	delete[] rHtA[i];
	//}
	//delete[] rHtA;
	//for(int i = 0; i < m; i++)
	//{
	//	delete[] ArHtA[i];
	//}
	//delete[] ArHtA;

	double** tmp2 = new double*[n];
	for(int i = 0; i < n; i++)
	{
		tmp2[i] = new double[n];
	}
	twoMatrixMultiply(rTmp,n,m,"N", AtH,m,n,"N",tmp2);
	//释放内存空间
	/*for(int i =0; i < n; i++)
	{
		delete[] tmp1[i];
	}
	delete[] tmp1;*/
	for(int i =0; i < m; i++)
	{
		delete[] AtH[i];
	}
	delete[] AtH;
	//两个矩阵相减
	for(int i = 0; i < n; i++)
	{
		for(int j = 0; j < n; j++)
		{
			Q[i][j] = reverseH[i][j] - tmp2[i][j];
		}
	}

	//=================内存清理==================
	for(int i =0; i < n; i++)
	{
		delete[] tmp2[i];
	}
	delete[] tmp2;
	//===========================================
}

void QuadEquLag::calculateR(double** reverseH, int& m, int& n, double** A,double** R)
{
	//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
	//================计算H逆乘A转置=====================
	double** rHtA = new double*[n]; //n*m维
	for(int i = 0; i < n; i++)
	{
		rHtA[i] = new double[m];
	}
	twoMatrixMultiply(reverseH,n,n,"N",A,m,n,"T",rHtA);
	//=====================================================
	//===============计算A乘H逆乘A转置=====================
	double** ArH = new double*[m]; //m*n维
	for(int i = 0; i < m; i++)
	{
		ArH[i] = new double[n];
	}
	twoMatrixMultiply(A,m,n,"N",reverseH,n,n,"N",ArH);
	double** ArHtA = new double*[m];
	for(int i = 0; i < m; i++)
	{
		ArHtA[i] = new double[m];
	}
	twoMatrixMultiply(ArH,m,n,"N",A,m,n,"T",ArHtA);
	inverseMatrix(ArHtA, m, m);
	//释放内存空间
	for(int i = 0; i < m; i++)
	{
		delete[] ArH[i];
	}
	delete[] ArH;

	twoMatrixMultiply(rHtA,n,m,"N",ArHtA,m,m,"N",R);

	//释放内存
	for(int i = 0; i < n; i++)
	{
		delete[] rHtA[i];
	}
	delete[] rHtA;
	for(int i = 0; i < m; i++)
	{
		delete[] ArHtA[i];
	}
	delete[] ArHtA;
}
void QuadEquLag::twoMatrixMultiply(double** a, int& aRow, int& aCol, char* aType, double** b, int& bRow, int& bCol, char* bType,double** out)
{
	if(strcmp(aType,"N") == 0 && strcmp(bType,"N") == 0)//a*b
	{
		for(int i = 0; i < aRow; i++)
		{
			for(int j = 0; j < bCol; j++)
			{
				double tmp = 0;
				for(int k = 0; k < aCol; k++)
				{
					tmp = tmp + a[i][k]*b[k][j];
				}
				out[i][j] = tmp;
			}
		}
	}
	else if(strcmp(aType,"N") == 0 && strcmp(bType,"T") == 0)//a*trans(b)
	{
		for(int i = 0; i < aRow; i++)
		{
			for(int j = 0; j < bRow; j++)
			{
				double tmp = 0.0;
				for(int k = 0; k < bCol; k++)
				{
					tmp = tmp + a[i][k]*b[j][k];
				}
				out[i][j] = tmp;
			}
		}
	}
	else if(strcmp(aType,"T") == 0 && strcmp(bType,"N") == 0)//trans(a)*b
	{

	}
	else if(strcmp(aType,"T") == 0 && strcmp(bType,"T") == 0)//trans(a)*trans(b)
	{

	}

}
void QuadEquLag::inverseMatrix(double** a, int& aRow, int& aCol)
{
	int N = aRow;
	double* A = new double[aRow * aCol];
	
	int LDA = aCol;
	int* IPIV = new int[N];
	int LWORK = N;
	double* WORK = new double[N];
	int INFO ;
	
	int k = 0;
	
	for(int i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			A[k++] = a[j][i];	
		}
	}
	int m = N; int n = m; 
	dgetrf_(&m, &n, A, &m, IPIV, &INFO);
	if(INFO != 0)
	{
		printf("dgetrf is wrong");
		assert(0);
	}
	INFO = 1;
	dgetri_(&N, A, &LDA, IPIV, WORK, &LWORK, &INFO);
	if(INFO != 0)
	{
		printf("dgetri is wrong");
		assert(0);
	}
	k = 0;
	for(int i = 0; i < N; i++)
	{
		for(int j = 0; j < N; j++)
		{
			a[j][i] = A[k++];	
		}
	}
	delete[] A;
}

//double* QuadProg::optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0)
//{
//	//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
//
//	
//
//	//vector<int> J0 ;
//	//
//	1)计算有效集J0，构造矩阵Aa0，设k=0
//	//for(int i = 0; i < m; i++)
//	//{
//	//		
//	//	double tmp = 0.0;
//	//	for(int j = 0; j < n; j++)
//	//	{
//	//		tmp = tmp + A[i][j] * x0[j];
//	//	}
//	//	if(tmp == b[i])
//	//	{
//	//		J0.push_back(i);
//	//	}
//	//}
//	//int nn=0;
//	//while(true)
//	//{
//	//	//2)计算gk，并判断rank[trans(Aak) gk] 时候等于 rank(trans(Aak))
//	//	double* gk = new double[n];
//	//	for(int i = 0; i < n; i++)
//	//	{
//	//		double tmp = 0.0;
//	//		for(int j = 0; j < n; j++)
//	//		{
//	//			tmp = tmp + H[i][j] * x0[j];
//	//		}
//	//		gk[i] = c[i] + tmp;
//	//	}
//	//	double** Aakgk = new double*[n];
//	//	for(int i = 0; i < n; i++)
//	//	{
//	//		Aakgk[i] = new double[J0.size() + 1];
//	//	}
//	//	
//	//	for(unsigned int i = 0; i < J0.size(); i++)
//	//	{
//	//		for(int j = 0; j < n; j++)
//	//		{
//	//			Aakgk[j][i] = A[J0[i]][j];
//	//		}
//	//	}
//	//	
//	//	int last = J0.size();
//	//	for(int j = 0; j < n; j++)
//	//	{
//	//		Aakgk[j][last] = gk[j];
//	//	}
//	//	
//	//	int col = J0.size() + 1;
//	//	int rankAakgk = rank(Aakgk,n,col);
//	//	col--;
//	//	int rankAak = rank(Aakgk, n, col);
//
//	//	if(rankAakgk == rankAak)
//	//	{
//	//	//3)如果gk中的值全大于零，则迭代停止，找到最优值，gk相当于罚参数
//	//		solveLinearEquations(Aakgk, n, col, gk);//gk作为输入参数和输出参数，最终的方程组的解存储在此
//	//	
//	//		for(int i = 0; i < n; i++)
//	//		{
//	//			delete[] Aakgk[i];
//	//		}
//	//		delete[] Aakgk;
//	//	
//	//		/*if(greaterThanZero(gk,n) == true)
//	//		{
//	//			delete[] gk;
//	//			break;
//	//		}*/
//	//		if(greaterThanZero(gk,n) == false)
//	//		{
//	//			updateActiveSet(J0, gk );//更新有效集，将gk中的负值对应的有效集去掉，如果gk中全部为负值，去掉最小的负值所对应的有效集
//	//			/*delete[] gk;*/
//	//		}
//	//	}
//	//
//	//	//4)求解dk
//
//	//	double** Aak = new double*[n];
//	//	int mTmp ;
//	//	double* bTmp = NULL;
//	//	if(J0.size() == 0)
//	//	{
//	//		for(int i = 0; i < n; i++)
//	//		{
//	//			Aak[i] = new double[1];
//	//			Aak[i][0] = 0.0;
//	//		}
//	//		mTmp = 1;
//	//		bTmp = new double[1];
//	//		for(unsigned int i = 0; i < 1; i++)
//	//		{
//	//			bTmp[i] = 0.0;
//	//		}
//	//	}
//	//	else
//	//	{
//	//		for(int i = 0; i < n; i++)
//	//		{
//	//			Aak[i] = new double[J0.size()];
//	//		}
//	//		for(unsigned int i = 0; i < J0.size(); i++)
//	//		{
//	//			for(int j = 0; j < n; j++)
//	//			{
//	//				Aak[j][i] = A[J0[i]][j];
//	//			}
//	//		}
//	//	
//	//		mTmp =(unsigned) J0.size();
//	//		bTmp = new double[J0.size()];
//	//		for(unsigned int i = 0; i < J0.size(); i++)
//	//		{
//	//			bTmp[i] = 0.0;
//	//		}
//	//	}
//	//	cout<<nn++<<endl;
//
//	//	double* dk = QuadEquLag::optimizer(H, mTmp, n, gk, Aak, bTmp);
//	//	if(equalZero(dk,n) == true && greaterThanZero(gk,n) == true)
//	//	{
//	//		delete[] gk;
//	//		break;
//	//	}
//	//	{
//	//		delete[] gk;
//	//	}
//	//	for(int i = 0; i < n; i++)
//	//	{
//	//		delete[] Aak[i];
//	//	}
//	//	delete[] Aak;
//	//	delete[] bTmp;
//	//	//5)求解afak，并设置x(k+1) = x(k) + afak*dk;
//	//	double* xk = x0;
//	//	int minIndex;
//	//	double afak = findAfak(J0, dk,A, b, xk,minIndex, m,n);
//	//	updataXk(xk,afak,dk,n);
//
//	//	//6)
//	//	if(afak < 1)
//	//	{
//	//		J0.push_back(minIndex);
//	//	}
//	//	sort(J0.begin(), J0.end());
//	//}
//	
//	return x0;
//}

double QuadProg::optimizer(double** H,int& m, int& n, double* c, double** A, double* b, double* x0)
{
	//H为n维对称阵，A为m*n矩阵，c为n维列向量，b为m维列向量
	vector<int> J0 ;
	bool come5 = true;
	//0)计算有效集J0，构造矩阵Aa0，设k=0
	for(int i = 0; i < m; i++)
	{
			
		double tmp = 0.0;
		for(int j = 0; j < n; j++)
		{
			tmp = tmp + A[i][j] * x0[j];
		}
		if(tmp == b[i])
		{
			J0.push_back(i);
		}
	}
	int nn=0;
	double* gk = new double[n];
	int nm = 0;
	while(true)
	{
		//1)计算gk，并判断rank[trans(Aak) gk] 时候等于 rank(trans(Aak))
		cout << nm++ <<endl;
		if(come5 == true)
		{
			come5 = false;

			
			for(int i = 0; i < n; i++)
			{
				double tmp = 0.0;
				for(int j = 0; j < n; j++)
				{
					tmp = tmp + H[i][j] * x0[j];
				}
				gk[i] = c[i] + tmp;
			}
		}
		//2)求解线性方程组[Q A'; A 0]*[dk v]' = [-gk 0]
		double* gkk = new double[n+J0.size()];
		double** QA = new double*[(n+J0.size())];
		for(unsigned int i = 0; i < n+J0.size(); i++)
		{
			QA[i] = new double[n+J0.size()];
		}
		for(int i = 0; i < n; i++)
		{
			for(int j = 0; j < n; j++)
			{
				QA[i][j] = H[i][j];
			}
		}
		for(unsigned int i = 0; i < J0.size(); i++)
		{
			for(int j = 0; j < n; j++)
			{
				QA[i+n][j] = A[J0[i]][j];
			}
		}
		for(unsigned int i = 0; i < J0.size(); i++)
		{
			for(int j = 0; j < n; j++)
			{
				QA[j][i+n] = A[J0[i]][j];
			}
		}
		for(unsigned int i = 0; i < J0.size(); i++)
		{
			for(unsigned int j = 0; j < J0.size(); j++)
			{
				QA[i+n][j+n] = 0.0;
			}
		}
		for(unsigned int i = 0; i < n+J0.size(); i++)
		{
			if(i < n)
				gkk[i] = -gk[i];
			else
				gkk[i] = 0.0;
		}
		int row = n + J0.size();
		int col = n + J0.size();

		double* gkkk = new double[n+J0.size()];
		for(unsigned int i = 0; i < n + J0.size(); i++)
		{
			gkkk[i] = gkk[i];
		}
		solveLinearEquations(QA, row, col, gkkk);
		for(unsigned int i = 0; i < n+J0.size(); i++)
		{
			delete[] QA[i];
		}
		delete[] QA;

		//3)如果dk=0,转到6）
		if(equalZero(gkkk,n) != true)
		{
		//4)计算步长
			int minIndex;
			double afak =findAfak(J0, gkkk, A ,  b, x0, minIndex,  m, n);
			if(afak < 1)
			{
				J0.push_back(minIndex);
			}
			sort(J0.begin(), J0.end());
		//5)更新xk
			updataXk(x0,afak,gkkk,n);
			come5 = true;
			delete[] gkk;
			delete[] gkkk;
			continue;
		}
		else if(equalZero(gkkk,n) == true)
		{
		//6）
			double* v = new double[J0.size()];
			for(unsigned int i = 0; i < J0.size(); i++)
			{
				v[i] = gkkk[n+i];
			}
			int num = J0.size();
			if(greaterThanZero(v,num) == true)
			{
				break;
			}
			else
			{
				updateActiveSet(J0, v);
			}
			delete[] v;
			
			delete[] gkkk;
			continue;

		}
	}
	delete[] gk;
	
	double obj = 0;
	double* o = new double[n];
	for(int i = 0; i < n; i++)
	{
		double tmp = 0.0;
		for(int j = 0; j < n; j++)
		{
			tmp = tmp + x0[j]*H[i][j];
		}
		o[i] = tmp;
	}
	double cx = 0.0;
	for(int i = 0; i < n; i++)
	{
		obj = obj + o[i] * x0[i] ;
		cx = cx + c[i] * x0[i];
	}
	delete[] o;

	return obj * 0.5 + cx;
		
		;
}

bool QuadProg::equalZero(double* in, int& n)
{
	for(int i = 0; i < n; i++)
	{
		if(in[i] > 0.0000000001)
			return false;
	}
	return true;
}

void QuadProg::updataXk(double* xk,double& afak, double* dk, int& n)
{
	for(int i = 0; i < n; i++)
	{
		xk[i] = xk[i] + afak * dk[i];
	}
}
double QuadProg::findAfak(vector<int>& activeSet, double* dk, double** A, double* b, double* xk, int& minIndex, int& m, int& n)
{
	
	double minVal = 0.0;
	vector<double> minVals;
	vector<int> minIndexs;
	
	for(int i = 0; i < m; i++)
	{
		vector<int>::iterator iter = find(activeSet.begin(),activeSet.end(),i);
		if(iter == activeSet.end()) //i不在有效集中
		{
			double aixk = 0.0;
			double aidk = 0.0;
			for(int j = 0; j < n; j++)
			{
				aidk = aidk + A[i][j] * dk[j];
				aixk = aixk + A[i][j] * xk[j];
			}
			if(aidk > 0.0)
			{
				double rating = (-aixk + b[i])/(aidk);
				minVals.push_back(rating);
				minIndexs.push_back(i);
			}
		}
	}
	minVal = minVals[0];
	minIndex = minIndexs[0];
	for(unsigned int i = 1; i< minVals.size(); i++)
	{
		if(minVals[i] < minVal)
		{
			minVal = minVals[i];
			minIndex = minIndexs[i];
		}
	}
	if(1 < minVal)
		return 1.0;
	else
		return minVal;
}
void QuadProg::updateActiveSet(vector<int>& activeSet, double* mu)
{
	int m =(int) activeSet.size();
	int minIndex = 0;
	double minVal = mu[0];
	for(int i = 1; i < m; i++)
	{
		if(mu[i] < minVal)
		{
			minVal = mu[i];
			minIndex = i;
		}
	}

	if(minVal < 0.0)
	{
		activeSet.erase(activeSet.begin() + minIndex);
	}
	
}
bool QuadProg::greaterThanZero(double* in, int& n)
{
	for(int i = 0; i < n; i++)
	{
		if(in[i] < 0.0)
			return false;
	}
	return true;
}
int QuadProg::rank(double** a,int& row, int& col, double tol)
{
	/*s = svd(A);rank = sum(s > tol);*/
	//===========SVD=============
	char jobu[]="A";
	char jobvt[]="A";
	int m = row;
	int n = col;
	double *b=new double[m*n];
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<m;j++)
		{
			b[i*m+j]=a[j][i];
		}
	}
	int lda=m;

	int minVal = min(m,n);
	double *s=new double[minVal];
	
	int ldu=m;
	double *u=new double[ldu*m];
	
	int ldvt=n;
	double *vt=new double[ldvt*n];
	int lwork=max(1,5*min(m,n));
	double *work=new double[lwork];
	
	int info;

	dgesvd_(jobu, jobvt, &m, &n, b, &lda, s, u, &ldu, vt, &ldvt, work, &lwork,&info);
	//===========================
	delete[] b;
	delete[] u;
	delete[] vt;
	int rank = 0;
	for(int i = 0; i < minVal; i++)
	{
		if(s[i] > tol)
			rank++;
	}

	delete[] s;
	return rank;
}

void QuadProg::solveLinearEquations(double** a, int& row, int& col, double* b)
{
	int N = row;
	int NRHS = 1;
	double* A = new double[(row) * (col)];

	int index = 0;
	for(int i = 0; i < row; i++)
	{
		for(int j = 0; j < col; j++)
		{
			A[index++] = a[j][i];
		}
	}

	int LDA = row;
	int* IPIV = new int[N];
	double* B = b;
	int LDB = max(1, N);
	int INFO;
	dgesv_(&N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO );

	delete[] A;
	delete[] IPIV;
}

测试主程序为：

#include "src\Optimizer\QP.h"
#include <fstream>
using namespace std;

int main(int argv, char** argc)
{
	double** H = new double*[400];
	double** A = new double*[400];
	for(int i = 0; i < 400; i++)
	{
		A[i] = new double[400];
		H[i] = new double[400];
	}

	ifstream readin;
	readin.open("./data/H.txt");
	if(readin.is_open())
	{
		for(int i = 0; i < 400; i++)
		{
			for(int j = 0; j < 400; j++)
			{
				double tmp;
				readin >> tmp;
				H[i][j] = tmp;
			}
		}
	}
	readin.close();

	double* f = new double[400];
	readin.open("./data/f.txt");
	if(readin.is_open())
	{
		for(int i = 0; i < 400; i++)
		{
			double tmp;
			readin >> tmp;
			f[i] = tmp;
		}
	}
	readin.close();

	for(int i = 0; i < 400; i++)
	{
		for(int j = 0; j < 400; j++)
		{
			A[i][j] = 0.0;
		}
		A[i][i] = 1.0;
	}

	double* b = new double[400];
	readin.open("./data/b.txt");
	if(readin.is_open())
	{
		for(int i = 0; i < 400; i++)
		{
			double tmp;
			readin >> tmp;
			b[i] = tmp;
		}
	}
	b[399] = 1000000000;
	readin.close();
	double* x0 = new double[400];
	for(int i = 0; i < 400; i++)
	{
		x0[i] = 0.5;
	}
	
	int m = 400;
	int n = 400;
	QuadProg* qu = new QuadProg();
	double obj = qu->optimizer( H,m,n,  f, A, b, x0);
	return 1;
}
本程序中用到了Lapack库和测试数据（Lapack是进行矩阵计算的库，Matlab中也采用了它）， 可以到此处下载，漏掉的Lapack.h文件见这里： 点击打开链接
文献学习(part16)--Oracle Based Active Set Algorithm for Scalable Elastic Net Subspace Clustering
小山羊的学习日志
12-02 401
学习笔记，仅供参考，有错必纠 文章目录Oracle Based Active Set Algorithm for Scalable Elastic Net Subspace ClusteringAbstractIntroductionPrior WorkPaper ContributionsElastic Net: Geometry and a New AlgorithmGeometric structure of the elastic net solutionA new active-set algo
一般凸二次规划的有效集算法 Matlab程序
jin的专栏
04-19 1万+
这是09年5月份很用心写的一个 运筹学中的一个程序，使得我对Matlab脚本语言有了一个深入的了解！ 但这篇很早已经在百度空间中发表，所以这次呢 作为在CSDN的第一篇博文，我选择它是有其实意义的，那就是要在 Linux 和 Socket 方面的研究将坚持同样的心态。 运筹学是大家所熟悉的，其实也就是问题的最优化求解，然而如何从程序角度去诠释复杂的算法，是一件非常有意思的事情。 由于程序中有大量的注释，所以这里就不再敖述，有问题可以留下评论 我们接着讨论。
凸二次规划的有效集.rar
04-25
最优化算法中的凸二次规划的有效集法，非常好用的程序，保证可以运行，包括四个M文件，有两个作用一样，只不过一个是直接执行，一个是放入命令窗口执行
二次规划及其MATLAB实现
最新发布
weidl001的博客
09-11 1100
二次规划（Quadratic Programming, QP）是一类重要的优化问题，其目标函数为二次函数，约束条件为线性不等式或等式。二次规划问题在工程、经济、金融等领域有广泛应用，如投资组合优化、人脸表情动画的权重求解、机械设计中的最优控制等。二次规划问题中的目标函数是二次型，约束条件为线性，因此问题的解法需要考虑二次型的性质（如正定性、半正定性等）。考虑一个典型的投资组合优化问题，目标是最小化投资组合的风险（即方差），并满足一定的收益约束。函数求解二次规划问题，输出的结果是投资组合的最优权重分配。
求解二次规划问题的拉格朗日及有效集方法(包含Matlab代码)
04-13
本资源主要内容：二次规划师非线性优化中的一种特殊情形，它的目标函数是二次实函数，约束函数都是线性函数。由于二次规划比较简单，便于求解（仅次于线性规划），并且一些非线性优化问题可以转化为求解一些列的二次规划问题，因此二次规划的求解方法较早引起人们的重视，称为求解非线性优化的一个重要途径。二次规划的算法较多，本文仅介绍求解等式约束凸二尺规划的拉格朗日方法以及求解一般约束凸二次规划的有效集方法。 本资源包含：《求解二次规划问题的拉格朗日及有效集方法》文档以及文档所用到的所有Matlab代码，非常适合初学者学习和研究！
有效集法介绍（Active Set Method）
青梅煮酒
10-29 8294
单纯性法（Simplex Method）是“线性规划之父”George Dantzig 最著名的成果，也是求解线性规划最有力的算法之一。而这一算法在求解二次规划（Quadratic Programming, QP）时的升级版就是有效集法（Active Set Method, ASM）。这两种算法的特点都是迭代点会循着约束边界前进，直到达到问题的最优点。本文对用于求解 QP 命题的 Primal A...
有效集法（Active Set），内点法（Interior-Point）及序列二次规划（SQP）法
热门推荐
更适合青年研究者的资源库！公众号：杰哥的无人驾驶便利店
12-01 1万+
有效集法（Active Set） 单纯性法（Simplex Method）是“线性规划之父”George Dantzig 最著名的成果，也是求解线性规划最有力的算法之一。而这一算法在求解二次规划（Quadratic Programming, QP）时的升级版就是有效集法（Active Set Method, ASM）。这两种算法的特点都是迭代点会循着约束边界前进，直到达到问题的最优点。 我们首先来讲一下约束部分： 假设 （就是这个向量一共有d个标量组成），则 的几何意象就是d维空间中的 d
Matlab_Active-set算法和实例
11-07
Active-set方法是一种在优化问题中解决凸二次规划（Quadratic Programming, QP）的有效算法。它主要应用于处理约束优化问题，尤其在那些包含线性不等式约束的情况下。在Matlab环境中，利用Active-set算法能够高效地...
Active_set_method2:求解二次规划问题的一种有用的活动集方法-matlab开发
05-31
总的来说，`Active_set_method2.m.zip`这个压缩包中的代码提供了一个使用MATLAB面向对象编程实现的二次规划问题的活动集方法。通过理解这个算法的原理和MATLAB的实现，我们可以更好地解决实际问题，并为优化问题的...
NX二次开发UF-ASSEM-set-active-arrangement 函数介绍
06-15
NX二次开发UF_ASSEM_set_active_arrangement 函数介绍，Ufun提供了一系列丰富的 API 函数，可以帮助用户实现自动化、定制化和扩展 NX 软件的功能。无论您是从事机械设计、制造、模具设计、逆向工程、CAE 分析等领域...
activeset_bc:具有约束的严格凸二次程序的快速活动集算法。-matlab开发
05-28
具有约束的严格凸二次程序的快速活动集算法。 [x，lam，ws，oc，cc] = activeset_bc（H，g，LB，UB，tol，maxiter）， 解决形式的严格凸二次程序最小0.5x'* H * x + g'* x，st LB <= x <= UB， 其中H = H'> 0...
有效集算法的几种MATLAB代码分享
09-14
二次规划是非线性优化中的一种特殊情形，它的目标函数是二次实函数，约束函数都是线性函数。 由于二次规划比较简单，便于求解（仅次于线性规划），并且一些非线性优化问题可以转化为求解一些列的二次规划问题，因此二次规划的求解方法较早引起人们的重视，称为求解非线性优化的一个重要途径。二次规划的算法较多， 本文打包了上传几种求解等式约束凸二尺规划的拉格朗日方法以及求解一般约束凸二次规划的有效集方法。 本资源包含：《求解二次规划问题的拉格朗日及有效集方法》文档以及文档所用到的所有Matlab代码，非常适合初学者学习和研究！
qp有效集算法
02-28
采用有效集策略的二次规划算法实现，c++语言。成功应用于预测控制在线控制律的求解
凸二次规划的详细解释以及求解步骤
05-15
凸二次规划有效集解法的解释以及求解步骤，解释了凸集集合并且给出了图示，希望能给读者很好的解释
算法大全(matlab算法集)
01-06
各种各样的计算方法在matlab中的实现
不等式约束二次规划——有效集法
w_w_y的博客
05-09 6913
不等式约束二次规划——有效集法
论文阅读笔记（5）：Oracle Based Active Set Algorithm for Scalable Elastic Net Subspace Clustering
weixin_44876302的博客
07-18 616
论文阅读笔记（5）：Oracle Based Active Set Algorithm for Scalable Elastic Net Subspace Clustering，基于Oracle的可伸缩弹性网络子空间聚类主动集算法摘要 注：16年CVPR，老文章但是个人觉得写得很好而且系统 摘要 ...
有效集 matlab代码
weixin_33871366的博客
06-28 2121
%有效集 function activeset H=[2 -1; -1 4]; c=[-1 -10]'; Ae=[ ]; be=[ ]; Ai=[-3 -2; 1 0; 0 1]; bi=[-6 0 0]'; x0=[0 0]'; [x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) function ...
数学建模 | 优化 | SQP、Active-Set、信赖域反射算法三种算法比较
m0_74083602的博客
08-26 628
信赖域反射算法可以理解为“稳扎稳打”的方法。它不会贸然做大步的改变，而是设定一个“信赖域”（相当于一步之内能走的最大距离），在这个区域内寻找最优解。如果当前的解不够好，它会调整信赖域的大小，并继续优化。有效集法有点像“挑挑拣拣”的策略。它先猜测哪些约束在最优解处会“紧贴”或“激活”，然后只针对这些“激活”的约束来求解问题。如果猜得不对，就调整并重复这个过程。SQP 可以理解为一种“逐步优化”的方法。它将一个复杂的非线性优化问题拆分成一系列更简单的二次规划问题，然后逐步求解这些简单问题，最终逼近最优解。
数值优化：二次规划与 active set 方法
"该资源是一份关于数值优化的讲义，涵盖了二次规划、有效集方法以及序列二次规划法。由Che-Rung Lee撰写，日期为2011年4月27日。主要内容包括二次规划的通用形式、拉格朗日乘子法、Karush-Kuhn-Tucker (KKT) 条件...
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值