二次规划问题和MATLAB函数quadprog的使用

二次规划(Quadratic Programming,QP)问题的一般形式为：

其中，，为纯量，为阶对称矩阵。易知二次规划的Hesse矩阵等于。如果为半正定矩阵，则称此规划为凸二次规划，否则为非凸规划。对于凸二次规划，目标函数q(x)是一个凸函数。如果有至少一个向量x满足约束而且q(x)在可行域有下界，二次规划问题就有一个全局最小值x。 如果G是正定矩阵，则称此规划为严格凸二次规划，此时全局最小值就是唯一的。如果G=0，二次规划问题就变成线性规划问题。根据优化理论，一个点x 成为全局最小值的必要条件是满足Karush-Kuhn-Tucker(KKT)条件。当q(x)是凸函数时，KKT条件也是充分条件。

---

根据约束条件的不同，二次规划可分为等式约束二次规划问题和不等式约束二次规划问题。等式约束二次规划问题即只含有等式约束，常见的解法有直接消去法、广义消去法、拉格朗日(Lagrange)法；对于不等式约束二次规划问题，其基本思想是把不等式约束转化为等式约束再求解，常见解法有有效集(active set)方法，有效集方法在每步迭代中把有效约束作为等式约束，然后可以用拉格朗日法求解，重复直到求得最优解。

        很多学者专门研究各类二次规划的求解方法，如文献[4][5]，对于非数学专业的的人来讲更重要的是怎么把二次规划当作一个工具去使用它。

【例】求如下二次规划问题。

【分析】首先应该把目标函数表示成如下矩阵形式：

                                                              

这里细说一下如何写成矩阵形式。(根据规划问题求解H矩阵可以用hessian函数直接求得，后面会说)

        首先，向量x是很容易写出的，因为f(x)包含两个变量x1和x2，因此

                                                                 

       其次，向量f只与两个变量x1和x2的一次项有关，所以，因此

                                                              

最后，矩阵H只与两个变量x1和x2的二次项有关，所以​​​​​​​，这里要注意的是不同于二次型，这里有个系数1/2，所以矩阵H的元素是二次型中的矩阵元素大小的两倍。给出一个规律：设矩阵H第i行第j列的元素大小为H(i,j)，二次项xixj的系数为a(i,j)，则

                                                          

本例中，，这是由于x1的平方项（即x1x1）系数为1/2，所以第1行第1列的元素为1=2*(1/2)，x2的平方项（即x2x2）系数为1，所以第2行第2列的元素为2=2*1，x1x2项（即x2x1）的系数为-1，所以第1行第2列和第2行第1列的元素均为-1。

        目标函数搞定之后，下面来看约束条件部分，约束条件应该写成如下形式：

                                                          

本例中约束条件只有不等式约束，因此Aeq和beq为空，对于A和b很容易就可以得出来：

                                               

 而约束条件中对变量x1和x2只给出下限，没有给上限，因此ub为空，

得到了所有的参数，将参数输入MATLAB，编程如下：（代码是直接在Command Window中一行一行录入的，所以每行前面有符号“>>”）

>> H = [1 -1; -1 2];
>> g= [-2; -6];
>> A = [1 1; -1 2; 2 1];
>> b = [2; 2; 3];
>> lb = [0; 0];
>> [x,fval,exitflag,output,lambda] = quadprog(H,g,A,b,[],[],lb)

输出结果，此处只展示x;x = 0.6667 1.3333

求取H的简单方法：

H = hessian(f,[x,y]);其中f为求解方程，[x,y]为变量，

注意需要转换为double类型：H = double(H);

求解一次项系数 F 
% fexp = expand(f); 
fcol = collect(f,{‘x’,’y’}); % x^2 - 4*x + y^2 - 6*y + 18 
disp(fcol);