求解二次规划问题的拉格朗]日及有效集方法——最优化方法课程实验报告学 院:数学与统计学院班 级:硕2041班姓 名：王彭学 号：28指导教师：阮小娥同组人：钱东东求解二次规划问题的拉格朗日及有效集方法摘要二次规划师非线性优化中的一种特殊情形，它的目标函数是二次实函数，约 束函数都是线性函数由于二次规划比较简单，便于求解（仅次于线性规划）， 并且一些非线性优化问题可以转化为求解一些列的二次规划问题，因此二次规划 的求解方法较早引起人们的重视，称为求解非线性优化的一个重要途径二次规 划的算法较多，本文仅介绍求解等式约束凸二尺规划的拉格朗日方法以及求解一 般约束凸二次规划的有效集方法关键字：二次规划，拉格朗日方法，有效集方法目录】摘要 错误!未定义书签1等式约束凸二次规划的解法 错误!未定义书签问题描述 错误!未定义书签拉格朗日方法求解等式约束二次规划问题 错误!未定义书签拉格朗日方法的推导 错误!未定义书签拉格朗日方法的应用 错误!未定义书签2 一般凸二次规划问题的解法 错误!未定义书签问题描述 错误!未定义书签有效集法求解一般凸二次规划问题 错误!未定义书签有效集方法的理论推导 错误!未定义书签。
有效集方法的算法步骤 错误!未定义书签有效集方法的应用 错误!未定义书签3总结与体会 错误!未定义书签4附录 错误!未定义书签拉格朗日方法的matlab程序 错误!未定义书签有效集方法的Matlab程序 错误!未定义书签1等式约束凸二次规划的解法问题描述min — xtHx + ctx,我们考虑如下的二次规划问题1 2s.t. Ax = b其中H e Rnxn对称正定，A e Rmxn行满秩，c, X G Rn，b e Rm拉格朗日方法求解等式约束二次规划问题拉格朗日方法的推导首先写出拉格朗日函数：L(x,人)=xtHx + ctx - Xt (Ax - b)V L(x, X) = 0,V/(x, X) = 0得到方程组Hx - ATX = -c 一 Ax = -b.将上述方程组写成分块矩阵形式:-At-a 0 _||_X「|_- b我们称伤处方程组的系数矩阵H-A-At0为拉格朗日矩阵卜面的定理给出了线性方程组有唯一解的充分条件定理1设H e Rmxn对称正定，A e Rmxn行满秩若在问题的解X*处满足二 阶充分条件，即dTHd > 0, Vd e Rn,d0, Ad = 0,则线性方程组的系数矩阵非奇异，即方程组（）有唯一解。
其中，方程组为（） 对应的齐次方程组：[H-At「d「-A0 _Lv J=0（）.&下面，我们来推导方程的求解公式根据定理1，拉格朗日矩阵必然是非奇 异的，故可设其逆为一 H-At=一 G-Bt --A0JL- Bc由恒等式\H-At ]]G-Bt -=-1n0 一nxm-A0 JL- bC JL0mxnIm可得HG + At B = I-HBt - AtC = 0nxm-AG = 0mxnABt = Im于是由上述四个等式得到矩阵G,B,C的表达式G = H-1 - H-1 At (AH-1 At )-i AH-1, (1.5)B = (AH -1 At ) -1 AH -1, (1.6)(1.7)C = -( AH -1 At ) -1.因此，由可得解得表达式X=一 Bt[-c]="-Gc + BTb—力-Bc_- b jBc - Cb(1.8)其中，G,B,C分别由,,给出下面给出X和厂的另一种等价表达式设xk是问题的任一可行点，即xk满 足Axk = b而在此点处目标函数的梯度为gk = W （七）=叫+c，利用xk和gk，气一 Gg(1.9)可将改写为Bg拉格朗日方法的应用（1） 拉格朗日方法的Matlab程序见附录。
2） 利用拉格朗日方法求解下列问题：min x2 + 2x2 + x2 — 2x x + x , st. x + x + x = 4, 2x - x + x = 2.解容易写出一 2-20「「11「「4「H =-240,c =0,A =,b =2-11」2_ 002」1_1」1—」 —,利用Matlab程序求解该问题可以结果如下:1.9000909090909091.9545454545454550. 136363636363637Lam =2.636363636363636-1.363636363036363£val =3.9772727272727282一般凸二次规划问题的解法问题描述考虑一般二次规划I 1 Tmin —xtHx + ctx,2s.t. aTx一b = 0,i e E = {1,•••,/},, (2.1)i iaTx-b > 0,i e I = {l +1,…,m}ii其中H是n阶对称阵记I(x*) = {i I aTx* -b = 0,i e I}，下面的定理给出了问题的一个最优性充要条件定理2 x*是二次规划问题的局部极小点当且仅当(1)存在人* e Rm，使得‘Hx+ c-& a-& a = 0,ieE ielaTx* - b = 0, i e E,aTx* - b > 0, i e I,X* > 0, i e I; X* = 0, i e I \ I ( x *)ii⑵记S = {d e Rn \{0}ldTa. = 0,ie E;dTa. > 0,i e I(x*);dTa. = 0,i e I(x*)且人* > 0}.则对于任意的d e S，均有dTHd > 0.容易发现，问题是凸二次规划的充要条件是H半正定。
此时，定理2的第二 部分自然满足注意到凸优化问题的局部极小点也是全局极小点的性质，我们有卜面的定理:定理3 x*是凸二次规划的全局极小点的充要条件是x*满足KT条件，即存在人* e Rm，使得"+ c-&a =0,ieE ieIaTx* - b = 0, i e E,i iaTx* - b > 0,i e I,iiX* > 0,i e I;X* = 0,i e I \ I(x*).ii有效集法求解一般凸二次规划问题有效集方法的理论推导首先引入下面的定理，它是有效集方法理论基础定理4设x*是一般凸二次规划问题的全局极小点，且在x*处的有效集为S(x*)= E I(x*)，则x*也是下列等式约束凸二次规划(2.2)I 1 丁口 , T min — xt Hx + ctx, 2st. aTx-b = 0,i e S(x*).i的全局极小点从上述定理可以发现，有效集方法的最大难点是事先一般不知道有效集S(x*)，因此只有想办法构造一个集合序列去逼近它，即从初始点x0出发，计算 有效集S(x0)，解对应的等式约束子问题重复这一做法，得到有效集序列{S(七)},k = 0,1,…，使之S(七)T s3*)，以获得原问题的最优解。
基于上述定理，我们分4步来介绍有效集方法的算法原理和实施步骤第一步，形成子问题并求出搜索方向dk.设七是问题的一个可行点，据此确定相应的有效集sk其中 I(x ) = {i I aTXk - b = 0, i e I}.求解相应的aT x子问题2 xtHx + ctx, st. aTx - b = 0, i e S .min(2.3)上述问题等价于min q (d) = -2drHd + grd,s.t. aTd = 0, i e S .(2.4)其中 k=Gxk + C-设求出问题的全局极小点为dk,人k是对应的拉格朗日乘子第二步，进行线性搜索确定步长因子ak .假设dk丰0分两种情形讨论k + dk是问题的可行点，即aT (x + d ) - b = 0, i e E 及 aT (x + d ) - b > 0, i e I.i k k i i k k i则令a = 1, x = x + d .(2) 若xk + dk不是问题的可行点，则通过线性搜索求出下降最好的可行点 注意到目标函数是凸二次函数，那么这一点应该在可行域的边界上达到因此只 要求出满足可行条件的最大步长a即可k当i e S时，对于任意的a > 0，都有aTd = 0和aT (x +a d ) = aTx = b， k k i k i k k k i k i此时，ak > 0不受限制。
当i冬S/寸，即第i个约束是严格的不等式约束，此时要 求a满足aT闵+a kdk) > b ,艮口a a’d > b - aTx , i e S .k i k i i k k注意到上式右端非正，故当aTd k > 0时，上式恒成立而当aTd k <0时，由上式可解得a < b 一件.k aTd故有a =a=min ] 一ai Xk | a’d < 0 >. k 〔 E ，k 〕合并第(1)(2)可得a = min{1,a k} (2.5).第三步，修正S. .当a k=1时，有效集不变，即SkH：= Sk .而当a k < 1时，故 a’ (x +a d ) = bik k k k ik—b - a’xak ―以k 'kaTdk k， i k k因此在xk+(处增加了一个有效约束，即Sk+i := Sk {ik}.第四步，考虑dk = 0的情形此时xk是问题的全局极小点若这是对应的不 等式约束的拉格朗日乘子均为非负，则xk也是问题的全局极小点，迭代终止 否则，如果对应的不等式约束的拉格朗日乘子有负的分量，那么需要重新寻找一个下降可行方向设七广0, e I亳)，现在要求一个下降可行方向 dk，满足grdk < 0且a’dk = 0, Vj e E;a^d^ > 0, Vj e I(x^)，为简便计算，按下述方式选取d^ :a\ (x + d ) > b ,a (x + d ) = b , Vj e S , j。
j ,j k k j k k| ardk > 0,队=0：月j e Sk, j 壬 j, a*)另一方面，注意到七是子问题的全局极小点，故有"A?人:其中A广("les :"广Les -k k从而，gTd = XtAtd .由知k k k k kATd = Z (aTd )e = (aT d )e ,jeS^ J h h于是有gTd = Xt (aT d )e, =Mk (m d ) < 0.上式表明，由确定的d是一个下降可行方向因此，令S ' = S \{j }，则修正 k k k k后的子问题min q (d) = 2dTHd + gTd,s.t. aTd = 0, i e S '的全局极小点必然是原问题的一个下降可行方向有效集方法的算法步骤经过上面的分析和推导，我们现在可写出有效集方法的算法步骤:(1) 选取初值给定初始可行点x0 e Rn，令k := 0 .(2) 解子问题确定相应的有效集S广E I(xk).求解子问题min q (d) = 2dTHd + gTd,s.t. aTd = 0, i e S ,得极小点dk和拉格朗日乘子向量人•若 dk0转步骤（4）；否则，转步骤（3）。
3）检验终止准则计算拉格朗日乘子其中gk=Hx + c,Bk = (AH -i A.t ) -i AH -i,=(ai )ieSk(X ) = min {(X ) }.k t ieI (xR k 1若(Xk)则xk是全局极小点，停算否则，若（X ）则令S := Sk \{t}，转步骤（2）4）确定步长a .令a = min{1,ak}，其中a = min] ~~七 ' | aTd < 01.k ieSk [ aTdk令x ：= x +a d .（5）若a = 1则令S := S ；kk+1否则，若a上＜ 1，则令Sk i = Sk{/J，其中匕满足— b - aT xk aT djk k|（6）令k := k +1，转步骤（1）有效集方法的应用（1）有效集法的Matlab程序见附录2)用有效集方法求解下列二次规划问题:min f (x) = x 2 — xx + 2x 2 — x 一 10x ,1 1 2 2 1 2< s.t. — 3x 一 2x 2 -6气 > 0, x2 > 0.解首先确定有关数据：「- 3—2—62—1-—1 -H =,c =,Ae = [], be = [], Ai =10, bi =0—14—101——1」 —1010利用Matlab计算可得结果如下:X =0. 50002. 250CT=lambda =output -0.750Cfva.1: -13.7500iter: 83总结与体会通过本次实验，笔者对求解等式约束凸二次规划问题的拉格朗日方法及一般 约束条件下凸二次规划问题的有效集方法有了较深的认识和了解。
感谢阮老师的悉心教诲和指导，使得笔者对最优化课程中的理论推导、算法 步骤及编程都比较熟悉，对以后的科研工作有很好的指导和借鉴意义4附录拉格朗日方法的matlab程序(1) 拉格朗日算法函数%本程序用拉格朗日方法求解灯饰约束条件的二次规划问题function [x,lam,fval]=qlag(H,A,b,c)%功能：用拉格朗日方法求解等式约束二次规划：—% min f(x)=*x'Hx+c'x， . Ax=b%输入：H,c分别是目标函数的矩阵和向量，A,b分别是约束条件中的矩阵和向量%输出：(x,lam)是KT点，fval是最优值IH=inv(H);AHA=A*IH*A';IAHA=inv(AHA);AIH=A*IH;G=IH-AIH'*IAHA*AIH;B=IAHA*AIH;C=-IAHA;x=B'*b-G*c;lam=B*c-C*b;fval=*x'*H*x+c'*x;⑵拉格朗日算法求解等式约束凸二次规划问题主程序：H=[2 -2 0;-2 4 0;0 0 2];c=[0 0 1]';A=[1 1 1;2 -1 1];b=[4 2]';[x,lam,fval]=qlag(H,A,b,c)有效集方法的Matlab程序(1)有效集方法函数%本程序主要适用于求解一般约束条件下的凸二次规划问题function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)%功能：用有效集方法解一般约束二次规划问题%min f(x)=*x'*H*x+c'*x,% . a'_i*x-b_i=0,(i=1,...,l),】% a'_i*x-b_i>=0,(i=l+1,...，m)%输入：x0是初始点，H,c分别是目标函数二次矩阵和向量；%Ae=(a_1,...,a_l)',be=(b_1,...,b_l);%Ai=(a_{l+1},...,a_m),bi=(b_{l+1},...,b_m)'.%输出：x是最优解，lambda是对应的乘子向量；output是结构变量% 输出极小值f(x),迭代次数k等信息，exitflag是算法终止类型%初始化epsilon=;err=;k=0;x=x0;n=length(x);kmax=;ne=length(be);ni=length(bi);lamk=zeros(ne+ni,1);(index=ones(ni,1);for i=1:niif Ai(i,:)*x>bi(i)+epsilonindex(i)=0;end%算法主程序while k<=kmax%求解子问题Aee=[];(if ne>0Aee=Ae;endfor j=1:niif index(j)>0Aee=[Aee;Ai(j,:)];endendgk=H*x+c;[m1,n1]=size(Aee);{[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1));if norm(dk)<=erry=;if length(lamk)>ne[y,jk]=min(lamk(ne+1:length(lamk)));endif y>=0exitflag=0;elseexitflag=1;for i=1:niif index(i)&&(ne+sum(index(1:i)))==jkindex(i)=0;break;endendendk=k+1;elseexitflag=1;%求步长alpha=;tm=;for i=1:niif (index(i)==0)&&(Ai(i,:)*dk<0)tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if tm10rb=Ae*ginvH*c+be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;(3) 有效集方法求解一般约束凸二次规划问题主程序clc;clear;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)。