求解二次规划问题的拉格朗日及有效集方法.doc

求解二次规划问题的拉格朗日及有效集方法

学 院：数学与统计学院
班 级：硕2041班
姓 名：王彭
学 号：3112054028
指导教师：阮小娥
同 组 人：钱东东
——最优化方法课程实验报告
求解二次规划问题的拉格朗
日及有效集方法

求解二次规划问题的拉格朗日
及有效集方法
摘要
二次规划师非线性优化中的一种特殊情形，它的目标函数是二次实函数，约束函数都是线性函数。由于二次规划比较简单，便于求解（仅次于线性规划），并且一些非线性优化问题可以转化为求解一些列的二次规划问题，因此二次规划的求解方法较早引起人们的重视，称为求解非线性优化的一个重要途径。二次规划的算法较多，本文仅介绍求解等式约束凸二尺规划的拉格朗日方法以及求解一般约束凸二次规划的有效集方法。

关键字：二次规划，拉格朗日方法，有效集方法。
【目录】
摘要 - 1 -
1 等式约束凸二次规划的解法 - 3 -
1.1 问题描述 - 3 -
1.2 拉格朗日方法求解等式约束二次规划问题 - 3 -
1.2.1 拉格朗日方法的推导 - 3 -
1.2.2 拉格朗日方法的应用 - 4 -
2 一般凸二次规划问题的解法 - 5 -
2.1 问题描述 - 5 -
2.2 有效集法求解一般凸二次规划问题 - 6 -
2.2.1 有效集方法的理论推导 - 6 -
2.2.2 有效集方法的算法步骤 - 9 -
2.2.3 有效集方法的应用 - 10 -
3 总结与体会 - 11 -
4 附录 - 11 -
4.1 拉格朗日方法的matlab程序 - 11 -
4.2 有效集方法的Matlab程序 - 11 -

1 等式约束凸二次规划的解法
1.1 问题描述
我们考虑如下的二次规划问题
(1.1)
其中对称正定，行满秩，，。
1.2 拉格朗日方法求解等式约束二次规划问题
1.2.1 拉格朗日方法的推导
首先写出拉格朗日函数：
，(1.2)
令
，
得到方程组

将上述方程组写成分块矩阵形式：
(1.3)
我们称伤处方程组的系数矩阵

为拉格朗日矩阵。
下面的定理给出了线性方程组(1.1)有唯一解的充分条件。
定理1 设对称正定，行满秩。若在问题(1.1)的解处满足二阶充分条件，即

则线性方程组(1.4)的系数矩阵非奇异，即方程组（1.4）有唯一解。其中，方程组(1.4)为（1.1）对应的齐次方程组：
（1.4）.
下面，我们来推导方程(1.3)的求解公式。根据定理1，拉格朗日矩阵必然是非奇异的，故可设其逆为
.
由恒等式

可得

于是由上述四个等式得到矩阵的表达式

因此，由(1.3)可得解得表达式

其中，分别由(1.5),(1.6),(1.7)给出。
下面给出和的另一种等价表达式。设是问题(1.1)的任一可行点，即满足。而在此点处目标函数的梯度为，利用和，可将(1.8)改写为

1.2.2 拉格朗日方法的应用
（1）拉格朗日方法的Matlab程序见附录。
（2）利用拉格朗日方法求解下列问题：

解 容易写出

利用Matlab程序求解该问题可以结果如下：

2 一般凸二次规划问题的解法
2.1 问题描述
考虑一般二次规划

其中是阶对称阵。记，下面的定理给出了问题(2.1)的一个最优性充要条件。
定理2 是二次规划问题(2.1)的局部极小点当且仅当
（1）存在，使得

(2)记

则对于任意的，均有.
容易发现，问题(2.1)是凸二次规划的充要条件是半正定。此时，定理2的第二部分自然满足。注意到凸优化问题的局部极小点也是全局极小点的性质，我们有下面的定理：
定理3 是凸二次规划的全局极小点的充要条件是满足条件，即存在，使得

2.2 有效集法求解一般凸二次规划问题
2.2.1 有效集方法的理论推导
首先引入下面的定理，它是有效集方法理论基础。
定理4 设是一般凸二次规划问题的全局极小点，且在处的有效集为，则也是下列等式约束凸二次规划

的全局极小点。
从上述定理可以发现，有效集方法的最大难点是事先一般不知道有效集，因此只有想办法构造一个集合序列去逼近它，即从初始点出发，计算有效集，解对应的等式约束子问题。重复这一做法，得到有效集序列使之，以获得原问题的最优解。
基于上述定理，我们分4步来介绍有效集方法的算法原理和实施步骤。
第一步，形成子问题并求出搜索方向.设是问题(2.1)的一个可行点，据此确定相应的有效集，其中求解相应的子问题

上述问题等价于

其中设求出问题(2.4)的全局极小点为是对应的拉格朗日乘子。
第二步，进行线性搜索确定步长因子.假设，分两种情形讨论。
(1) 若是问题(2.1)的可行点，即

则令
(2) 若不是问题(2.1)的可行点，则通过线性搜索求出下降最好的可行点。注意到目标函数是凸二次函数，那么这一点应该在可行域的边界上达到。因此只要求出满足可行条件的最大步长即可。
当时，对于任意的，都有和，此时，不受限制。当时，即第个约束是严格的不等式约束，此时要求满足，即

注意到上式右端非正，故当时，上式恒成立。而当时，由上式可解得

故有

合并第(1)(2)可得
.
第三步，修正.当时，有效集不变，即.而当时，
，
故，因此在处增加了一个有效约束，即.
第四步，考虑的情形。此时是问题(2.3)的全局极小点。若这是对应的不等式约束的拉格朗日乘子均为非负，则也是问题(2.1)的全局极小点，迭代终止。否则，如果对应的不等式约束的拉格朗日乘子有负的分量，那么需要重新寻找一个下降可行方向。
设，现在要求一个下降可行方向，满足且，为简便计算，按下述方式选取：

即

另一方面，注意到是子问题(2.3)的全局极小点，故有

即

其中

从而，由(2.6)知

于是有

上式表明，由(2.6)确定的是一个下降可行方向。因此，令，则修正后的子问题

的全局极小点必然是原问题的一个下降可行方向。
2.2.2 有效集方法的算法步骤
经过上面的分析和推导，我们现在可写出有效集方法的算法步骤：
（1）选取初值。给定初始可行点，令.
（2）解子问题。确定相应的有效集.求解子问题

得极小点和拉格朗日乘子向量.若转步骤（4）；否则，转步骤（3）。
（3）检验终止准则。计算拉格朗日乘子

其中

令

若，则是全局极小点，停算。否则，若，则令，转步骤（2）。
（4）确定步长.令，其中

令
（5）若，则令；否则，若，则令，其中满足

（6）令，转步骤(1)。
2.2.3 有效集方法的应用
（1）有效集法的Matlab程序见附录。
（2）用有效集方法求解下列二次规划问题：

解 首先确定有关数据：

利用Matlab计算可得结果如下：

3 总结与体会
通过本次实验，笔者对求解等式约束凸二次规划问题的拉格朗日方法及一般约束条件下凸二次规划问题的有效集方法有了较深的认识和了解。
感谢阮老师的悉心教诲和指导，使得笔者对最优化课程中的理论推导、算法步骤及编程都比较熟悉，对以后的科研工作有很好的指导和借鉴意义。
4 附录
4.1 拉格朗日方法的matlab程序
(1) 拉格朗日算法函数

%本程序用拉格朗日方法求解灯饰约束条件的二次规划问题。
function [x,lam,fval]=qlag(H,A,b,c)
%功能：用拉格朗日方法求解等式约束二次规划：
% min f(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 -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)

4.2 有效集方法的Matlab程序
（1）有效集方法函数

%本程序主要适用于求解一般约束条件下的凸二次规划问题
function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)
%功能：用有效集方法解一般约束二次规划问题
%min f(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);
for i=1:ni
if Ai(i,:)*x>bi(i)+epsilon
index(i)=0;
end
end
%算法主程序
while k<=kmax
%求解子问题
Aee=[];
if ne>0
Aee=Ae;
end
for j=1:ni
if index(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));
if norm(dk)<=err
y=0.0;
if length(lamk)>ne
[y,jk]=min(lamk(ne+1:length(lamk)));
end
if y>=0
exitflag=0;
else
exitflag=1;
for i=1:ni
if index(i)&&(ne+sum(index(1:i)))==jk
index(i)=0;
break;
end
end
end
k=k+1;
else
exitflag=1;
%求步长
alpha=1.0;tm=1.0;
for i=1:ni
if (index(i)==0)&&(Ai(i,:)*dk<0)
tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);
if tm1<tm
tm=tm1;
ti=i;
end
end
end
alpha=min(alpha,tm);
x=x+alpha*dk;
%修正有效集
if tm<1
index(ti)=1;
end
end
if exitflag==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);
if m>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）有效集方法求解一般约束凸二次规划问题主程序

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)
- 13 -