1 .等式约束二次规划问题1)。
另外,求解对称矩阵、等式约束二次规划问题一般采用直接消元法、正交分解法、拉格朗日法。
用拉格朗日方法构造Eq.1)的拉格朗日乘子函数
在KKT的条件下得到:
写成行列的形式
系数矩阵称为Lagrange矩阵,对称不定。 如果Lagrange矩阵的逆存在
得到的最佳解.如果存在就能得到
如果不存在Lagrange矩阵的逆,就可以求出伪逆。
2 .等式和不等式约束二次规划问题
由于具有等式约束的二次规划问题已经有了成熟的解法,要求求解具有不等式约束的二次规划问题,自然有将不等式约束转化为等式约束的思路。 有效集合法Active Set Method )是求解具有不等式约束的二次规划问题的典型算法。首先看下问题的几何表述
可行域:目标问题的可行域是包含边界和内部的空间内的凸多面体。 例如,上图中的绿色多边形区域。 当等值面:矩阵g归一化时,目标函数的等值面是空间中一族同心超椭圆体,所有椭圆体相似且同向各轴同向)。 例如,上图中的同心椭圆曲线族。 越是内层椭圆体,目标函数的取值越好低)。 在优化问题:几何表示下,我们的目标是在可行范围内找到到达同心椭圆曲线族最内层的点。 例如,在上图中,p是无约束最优解,q是约束最优解。 添加约束后,最优解一定会出现在某不等式约束的边界上。
给定有效集:任一可行解,其有效集合是满足等号成立的限制条件的集合。 当然,有效集一定包含所有等式约束,也包括不等式约束的子集。 请注意,这里等号是成立的。 也就是说,如果当前解使不等式约束的等号成立,则此不等式约束包含在有效集中。 将最优有效集:最优解的有效集改写为最优有效集3358www.Sina.com/最优有效集的约束全部等式约束,舍弃其他约束,用Lagrange乘法直接求解即可。 虽然约束的数目有限,但是由于最佳有效集是其子集,所以最佳有效集存在有限的多个可能性。 于是自然想到了以下有名的暴力解读法——贫困法。
每次从整体约束中选取一个子集包括所有等式约束和部分不等式约束),称为此次尝试的。,求解工作集一个工作集合对应的子优化问题然后,检查该解是否为原始问题的可行解,即是否也满足工作集以外的约束。 是的,如果是,我会记录下来。 否则就扔掉。 调查所有可能的子集后,在所有记录的解中选择目标函数值最佳的即可。
穷举法可以在有限步骤内确实得到目标问题的解,但其运算量往往超出想象,特别是标量和不等式约束条件个数较多时,要求在以下两个方面进行改善。
利用3358www.Sina.com/凸二次规划的特点,建立了一组识别规则实际为KKT条件),可以直接判断可行解是否为最优解,所以算法在尝试最优解后立即停止,剩余可能性http://www.com 代替随机选择下一个要尝试的工作集,通过求解子以上两个改进,形成了效率大大提高的有效集法。
工作集对应的子优化问题。
在改进最优解的识别方式:有效集合法中,作为反复的出发点需要初始可行解。 由于可行解与目标函数无关,得到初始可行解的方法与线性规划完全一致。 典型的方法有两阶段法和大m法。 因此,我们假设已经得到了初始可行解x0,得到了有效集合,作为W0迭代求解。 迭代过程主要由两部分组成,一部分是寻找当前有效集的最优解,另一部分是更新有效集得到最优集。 重复k k=1,2,3,)步骤后,假设现在的解为xk,工作集合为wk,根据xk在可行领域的位置,可以分为4种情况。 如果33558www.Sina.com/xk满足KKT条件,则成为最佳解,算法停止。 33558www.Sina.com/xk在当前工作集等式约束集)下是最优的,但如果工作集中至少存在一个约束,则它原本是不等式约束,当xk移动到不等式约束的可行侧时,可以进一步降低目标函数值不移动迭代点即xk 1=xk脱离工作集的约束)恢复不等式约束)恢复不等式约束33558www.Sina.com/xk在当前工作集等式约束集)下仍可改进,但步骤已满出原问题的可行域) )将正好到达原可行域边界的位置作为下一个迭代点xk 1,同时将与冲突边界对应的不等式约束添加到工作集中。 情况3的例子是下图的x3 改进尝试各种可能性的顺序:
若 xk在当前工作集(等式约束集)下还可以改进,并且步长能够走满(即改进到最优仍未出原始可行域),那么就取这个最优点作为下一个迭代点 xk+1,并保持工作集不变。情形4的例子请参见下图中 x1和 x4。
下面是使用有效集求解带不等式约束优化问题的经典例子。
G = eye2,2)*2;d = [-2 -5]’;A = [-1 2; 1 2; 1 -2; -1 0; 0 -1];b = [2;6;2;0;0];Aeq = [-1 1];beq = [0.03];MaxIterations = 100;Tolerance = 1.0e-6;x0 = [0;1]; workset0= [3;5];x = active_set_methodx0,G,d, A, b, Aeq, beq, MaxIterations, Tolerance)%options = optimoptions@quadprog,’Algorithm’,’active-set’,’MaxIterations’,500);y = quadprogG, d, A, b, Aeq, beq, [], [], x0, options)function x = active_set_methodx0,G,d, Ai, bi, Aeq, beq, MaxIterations, Tolerance)%min 1/2*x’Gx+d’x%st. A*x <= b% Aeq*x = b ne = sizeAeq,1); ni = sizeAi,1); workset = zerosni,1); x = x0;%给定初始解 for i = 1:ni %将满足初始解的不等式都加入有效集 if absAii,:)*x0 -bii)) <= Tolerance workseti) = 1; end end for k = 1:MaxIterations %迭代求解 a = [Aeq; Aiworkset == 1,:)]; b_ = [beq; biworkset == 1,:)]; [x1, lambda]=kktG, d, a, b_); s = x1 – x; if norms,2) < Tolerance %得到的解是本次有效集中的最优解 x = x1; min_lambda = 0; if lengthlambda) > ne [min_lambda,index] = minlambdane+1:end)); end if min_lambda >= 0 %满足KKT条件,得到最优解 break; else %寻找新的有效集 for i=1:ni if workseti)&&sumworkset1:i))==index workseti)=0; break; end end end else %得到的解不是本次有效集中的最优解,继续优化 [alpha, index] = Alphaworkset, x, s, Ai, bi); x = x + alpha*s; if alpha <1 worksetindex) = 1; end end end function [s, lambda]=kktG, g, Ae, be) ws = sizeAe,1); kkt_A = [G, Ae’; Ae, zerosws,ws)]; kkt_B = [-g;be]; slambda = pinvkkt_A)*kkt_B; %kkt_A不一定可逆,需要求伪逆 dim = sizeG,2); s = slambda1:dim); lambda = slambdadim+1:end); end function [alpha, index] = Alphaworkset, x, s, Ai, bi) outset = findworkset == 0); Atp = Aioutset,:)*s; Atp_index = Atp>0; out_plus_ind = outsetAtp_index); alphatset = biout_plus_ind) – Aiout_plus_ind,:)*x)./AtpAtp_index); [alpha, index] = min[alphatset;1]); if alpha < 1 index = out_plus_indindex); end endend
参考:
https://zhuanlan.zhihu.com/p/29525367
https://zhuanlan.zhihu.com/p/26514613
https://blog.csdn.net/dymodi/article/details/50358278
https://zhuanlan.zhihu.com/p/50906541