求解二次规划问题的拉格朗日及有效集方法.docx
《求解二次规划问题的拉格朗日及有效集方法.docx》由会员分享，可在线阅读，更多相关《求解二次规划问题的拉格朗日及有效集方法.docx（18页珍藏版）》请在冰豆网上搜索。
求解二次规划问题的拉格朗日及有效集方法
求解二次规划问题的拉格朗
日及有效集方法
最优化方法课程实验报告
院：
数学与统计学院
求解二次规划问题的拉格朗日
及有效集方法
摘要
二次规划师非线性优化中的一种特殊情形，它的目标函数是二次实函数，约束函数都是线性函数。
由于二次规划比较简单，便于求解（仅次于线性规划），并且一些非线性优化问题可以转化为求解一些列的二次规划问题，因此二次规划的求解方法较早引起人们的重视，称为求解非线性优化的一个重要途径。
二次规划的算法较多，本文仅介绍求解等式约束凸二尺规划的拉格朗日方法以及求解一般约束凸二次规划的有效集方法。
关键字：
二次规划，拉格朗日方法，有效集方法。
【目录】
摘要
等式约束凸二次规划的解法
1.1问题描述
1.2拉格朗日方法求解等式约束二次规划问题
1.2.1拉格朗日方法的推导
1.2.2拉格朗日方法的应用
一般凸二次规划问题的解法
2.1问题描述
2.2有效集法求解一般凸二次规划问题
2.2.1有效集方法的理论推导
2.2.2有效集方法的算法步骤
-11-
-11-
4.1拉格朗日方法的matlab程序
4.2有效集方法的Matlab程序
1等式约束凸二次规划的解法
1.1问题描述
我们考虑如下的二次规划问题
[■.1T口T
min-xHx+cx,z>八\2（1.1）is.t.Ax=b
其中对称正定，R^^行满秩，c,x亡Rn，b亡Rm
1.2拉格朗日方法求解等式约束二次规划问题
1.2.1拉格朗日方法的推导
首先写出拉格朗日函数:
1
L（x,A）=—xtHx+cTx—A（Ax-b），（1.2）
2
jL（x,k）=0,J丄（x,Q=0，
得到方程组
Hx-A冬=-c,
-Ax=-b.
将上述方程组写成分块矩阵形式:
〔H[-A
我们称伤处方程组的系数矩阵
-at-
0
足二阶充分条件，即
dTHdaOM5,d工O,Ad=0,
则线性方程组（1.4）的系数矩阵非奇异，即方程组（1.4）有唯一解。
其中，方程组（1.4）为（1.1）对应的齐次方程组：
下面，我们来推导方程是非奇异的，
故可设其逆为
由恒等式
可得
LHa
（1.3）
叮/0（1'4）-
的求解公式。
根据定理1,拉格朗日矩阵必然
-At
0
Tg-bt
/i-BC
〔H
-at「
〔G
-Bt"
r1n0nxm
[-A
0
[—B
C
|_0mxi1m”
HG+ATB=ln
—AG=0m河
—HBt-AtC=0恤
ABt=lm.
于是由上述四个等式得到矩阵G,B,C的表达式
11T1T11
G=H--HA（AHA）AH
4T4A.
B=（AHA）AH,
C=-（AH°AT）°
因此，由（1.3）可得解得表达式
ELF-叮〔"
0」B「
「-Gc+BTb]Cj[-b」[Bc-Cb」
其中，G,B,C分别由（1.5）,（1.6）,（1.7）给出。
（1.5）
（16）
（1.7）
（1.8）
F面给出X和几的另一种等价表达式。
设xk是问题（1.1）的任一可行点，即
Xk满足AXk=b。
而在此点处目标函数的梯度为
gk=Nf（Xk）=HXk中c，利用Xk
和gk，可将（1.8）改写为
[x[_[xk-GgkL^j"iBgk.
（1.9）
1.2.2拉格朗日方法的应用
（1）拉格朗日方法的Matlab程序见附录。
（2）利用拉格朗日方法求解下列问题：
min
s.t.
222
Xr+2x2+X3-2x1X2+X3,
x^i+X2+X3=4,
2X1—X2+X3
=2.
解容易写出
f
•2
-2
0-
f
•0-
H=
-2
4
0
c=
0
[
0
0
2
[
1
aJ1
-1
b邛
L2J
利用Matlab程序求解该问题可以结果如下:
k909090909090909
1.954545454545455
0.135********3637
Ian=
2.630303636363636
-L363636363636363
fval=
3,97；272727272：
2S
2一般凸二次规划问题的解法
2.1问题描述
考虑一般二次规划
min1xTH^cTx,
2
ajx—bi=0,i-E={1,…，l},,
{s.t.
!
ajx-bi>O,iJ={l+1,…,m}
（2.1）
其中H是n阶对称阵。
记l（x*）={ilaTx*-bi=0,i迂1}，
F面的定理给出了问题
（2.1）的一个最优性充要条件。
定理2x是二次规划问题（2.1）的局部极小点当且仅当
（1）存在ZeRm，使得
Hx+c—2；入iNj—送Zia^0,
任ig
T*
aiX—bi=0,i亡E,
T*
aiX-bi>0,i€I,
入*>0,i"M*=0,飞II（x*）
*,,*
=0,i€|（x）且M>0}.
⑵记
S={dRn{0}|dTa^0,^E;dTa^0,^I（x*）;dTai则对于任意的d亡S，均有dTHd>0.
容易发现，问题（2.1）是凸二次规划的充要条件是H半正定。
此时，定理2的第二部分自然满足。
注意到凸优化问题的局部极小点也是全局极小点的性质，我们有下面的定理：
定理3X*是凸二次规划的全局极小点的充要条件是X*满足KT条件，即存
在力迂Rm，使得
Hx+C—无几idi—2入iai=0,
任iS
aTX-bj=0,i忘E,
aTx-bj>0,i引，
人>0,i引；k*=0,i引I（x*）.
2.2有效集法求解一般凸二次规划问题
2.2.1有效集方法的理论推导
首先引入下面的定理，它是有效集方法理论基础。
定理4设X是一般凸二次规划问题的全局极小点，且在X处的有效集为
S^*）=eU|（x*），则x*也是下列等式约束凸二次规划
I.1T.T
Imin2XH^CX\（2.2）
[st.a：
X—b=0,i亡S（x）.
的全局极小点。
从上述定理可以发现，有效集方法的最大难点是事先一般不知道有效集
S（x*），因此只有想办法构造一个集合序列去逼近它，即从初始点X0出发，计算
-6-
有效集S（X0），解对应的等式约束子问题。
重复这一做法，得到有效集序列
{S（Xk）},k二。
1，…，使之S（Xk）TS（X），以获得原问题的最优解。
基于上述定理，我们分4步来介绍有效集方法的算法原理和实施步骤。
第一步，形成子问题并求出搜索方向dk.设Xk是问题（2.1）的一个可行点，据此确定相应的有效集Sk=EU|（xk），其中I（Xk）={i|aiTXk-bi=0,i忘I}.求解相应的子问题
\.1T+T
（2.3）
（2.4）
min-XHx+cX,
\T2
[s.t.aiX-bi=0,i忘Sk.
上述问题等价于
1TT
minqk（d）=-dHd+gkd,
2
s.t.ajd=0,i忘Sk.
其中X=Xk+4©-GXk+c.设求出问题（2.4）的全局极小点为dk,Xk是对应的拉格朗日乘子。
（1）若Xk+dk是冋题（2.1）的可行点,
则令叭Hxkdt=Xk+dk.
⑵若Xk+dk不是问题（2.1）的可行点，则通过线性搜索求出下降最好的可行
点。
注意到目标函数是凸二次函数，那么这一点应该在可行域的边界上达到。
因此只要求出满足可行条件的最大步长otk即可。
当i-Sk时，对于任意的叭>0，都有ajdk=0和aiT（XkEkdk）=aiTXk=bi，
此时，叭＞0不受限制。
当浮Sk时，即第i个约束是严格的不等式约束，此时要求ak满足aT（Xk+ctkdk）3bi，即
GkaTdk^b—aTXk,i忘
注意到上式右端非正，故当aTdk二0时，上式恒成立。
而当aTd^O时，由上式
-7-
可解得
aidk
故有
s=7k=min卩―；'兀laTdkIaidk
合并第
（1）
（2）可得
o-k=min{1,otk}（2.5）.
第三步，修正Sk.当ak=1时，有效集不变，即Sy=Sk.而当otkCl时,
.T
-bik-aikXk叫"k——aikdk
故aT（X^akdk^bi，因此在Xk十处增加了一个有效约束，即
kk
Sk十：
=SkU{ik}.
第四步，考虑dk=O的情形。
此时Xk是问题（2.3）的全局极小点。
若这是对应
的不等式约束的拉格朗日乘子均为非负，则Xk也是问题（2.1）的全局极小点，迭代
终止。
否则，如果对应的不等式约束的拉格朗日乘子有负的分量，那么需要重新
寻找一个下降可行方向。
设几jkv。
'jk引（Xk），现在要求一个下降可行方向dk，满足gidk<0且
aTdk=O,Vj-EeTdk>0,巧-I（Xk），为简便计算，按下述方式选取dk:
aTk（Xk+dk）Abjk,aT（Xk+dk）=bjNj-Sk,jHj
fa^dkAO,
bdk5Sk,jHjk,3
另一方面，注意到xk是子问题（2.3）的全局极小点，故有
Hxk+c-送几:
ai=0,
id
gk=Ak\,
其中
Ak=（ai）iSk，几k=（人）i旳.
从而，gldk^^Aldk.由（2.6）知
Aldk=£（aTdk）ej=（aTkdk）ejk，
jSk
于是有
gidk
=沐（ajkdk）ej=Xjk（ajkdk）<0.
上式表明，由（2.6）确定的dk疋
曰一个下降可行方向。
因此，令Sk'=Sk{jk}，贝M修
正后的子问题
4
minqk（d）=—dTHd+gTd,
2
s.t.aTd=0,i壬Sk'
的全局极小点必然是原问题的一个下降可行方向。
222有效集方法的算法步骤
经过上面的分析和推导，我们现在可写出有效集方法的算法步骤:
（1）选取初值。
给定初始可行点xo-Rn，令k=0.
（2）解子问题。
确定相应的有效集Sk=eU|（xk）.求解子问题
minqk（dH1dTH^gTd,
s.t.ajd=0,i忘Sk,
得极小点dk和拉格朗日乘子向量几k•右dkK0转步骤（4）;否则，转步骤（
（3）检验终止准则。
计算拉格朗日乘子
'“k=Bkgk,
其中
gk=HXk+c,Bk=（人H~AT）~AkHt宀=@）赳.
@k）t=m（in）{0k）i}.
若仏k）tX0，则Xk是全局极小点，停算。
否则，若（S）t"，则令Sk=Sk
{t},
转步骤
（2）。
（4）确定步长ak.令ak=min{1,otk}，其中
ak=minJbXk|aTdk<0，.
0Iaidki
令Xk+：
=Xk+akdk.
（5）若Gk=1，则令SkH1=Sk；否则，若ak<1,则令S"=SkU{jk}，其中jk满足
.T
-bjk-SXk
ak=T
ajkdk
（6）令k=k+1，转步骤
（1）。
223有效集方法的应用
有效集法的Matlab程序见附录。
用有效集方法求解下列二次规划问题:
22
f（x）hXj—x1X2+2x2—x1-10x2,
st.—3x1—2x2——6
为>0,X2二0.
解首先确定有关数据:
利用Matlab计算可得结果如下:
exitflag=
O'
0.5000
2.2500
3总结与体会
通过本次实验，笔者对求解等式约束凸二次规划问题的拉格朗日方法及一般约束条件下凸二次规划问题的有效集方法有了较深的认识和了解。
感谢阮老师的悉心教诲和指导，使得笔者对最优化课程中的理论推导、算法步骤及编程都比较熟悉，对以后的科研工作有很好的指导和借鉴意义。
4附录
4.1拉格朗日方法的matlab程序
（1）拉格朗日算法函数
%本程序用拉格朗日方法求解灯饰约束条件的二次规划问题。
function[x,lam,fval]=qlag（H,A,b,c）
%功能：
用拉格朗日方法求解等式约束二次规划：
%minf（x）=0.5*x'Hx+c'x，s.t.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=0.5*x'*H*x+c'*x;
（2）拉格朗日算法求解等式约束凸二次规划问题主程序:
H=[2-20;-240;002];c=[001]';
A=[111;2-11];
b=[42]';
[x,lam,fval]=qlag（H,A,b,c）
4.2有效集方法的Matlab程序
（1）有效集方法函数
%本程序主要适用于求解一般约束条件下的凸二次规划问题
function[x,lamk,exitflag,output]=qpact（H,c,Ae,be,Ai,bi,xO）
%功能：
用有效集方法解一般约束二次规划问题
%minf（x）=0.5*x'*H*x+c'*x,
%s.t.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=1.0e-9;err=1.0e-6;
k=0;x=x0;n=length（x）;kmax=1.0e3;
ne=length（be）;ni=length（bi）;lamk=zeros（ne+ni,1）;
index=ones（ni,1）;
fori=1:
ni
ifAi（i,:
）*x>bi（i）+epsilon
index（i）=0;
end
end
%算法主程序
whilek<=kmax
%求解子问题
Aee=[];
ifne>0
Aee=Ae;
end
forj=1:
ni
ifindex（j）>0
Aee=[Aee;Ai（j,:
）];
end
end
gk=H*x+c;
[m1,n1]=size（Aee）;
[dk,lamk]=qsubp（H,gk,Aee,zeros（m1,1））;
ifnorm（dk）v=eiT
y=0.0;
iflength（lamk）>ne
[y,jk]=min（lamk（ne+1:
length（lamk）））;
end
ify>=0
exitflag=0;
else
exitflag=1;
fori=1:
ni
ifindex（i）&&（ne+sum（index（1:
i）））==jkindex（i）=0;
break;
end
end
end
k=k+1;
else
exitflag=1;
%求步长
alpha=1.0;tm=1.0;
fori=1:
ni
if（index（i）==0）&&（Ai（i,:
）*dk<0）tm1=（bi（i）-Ai（i,:
）*x）/（Ai（i,:
）*dk）;
iftm1tm=tm1;
ti=i;
end
end
end
alpha=min（alpha,tm）;
x=x+aIpha*dk;
%修正有效集
iftm<1
index（ti）=1;
end
end
ifexitflag==0
break;
end
k=k+1;
end
output.fval=0.5*x'*H*x+c'*x;
output.iter=k;
（2）求解子问题函数
%求解子问题
function[x,lambda]=qsubp（H,c,Ae,be）ginvH=pinv（H）;
[m,n]=size（Ae）;
ifm>0
rb=Ae*ginvH*c+be;
lambda=pinv（Ae*ginvH*Ae'）*rb;x=ginvH*（Ae'*lambda-c）;
else
x=-ginvH*c;
lambda=0;
end
（3）有效集方法求解一般约束凸二次规划问题主程序
cic;
clear;
H=[2-1;-14];
c=[-1-10]';
Ae=[];be=[];
Ai=[-3-2;10;01];
bi=[-600]';
xO=[O0]';
[x,lambda,exitflag,output]=qpact（H,c,Ae,be,Ai,bi,xO）
升级会员 