最优化方法报告二.docx
《最优化方法报告二.docx》由会员分享，可在线阅读，更多相关《最优化方法报告二.docx（17页珍藏版）》请在冰豆网上搜索。
最优化方法报告二
最优化方法实验报告
姓名：
班级：
学号：
一.从外罚，内罚，乘子法中选一编程
乘子法（PHR算法）
1.算法原理
乘子法是Powell和Hestenes于1969年针对等式约束优化问题同时独立提出的一种优化算法，后于1973年经Rockfellar推广到求解不等式约束优化问题。
其基本思想是：
从原问题的拉格朗日函数出发，再加上适当的罚函数，从而将原问题转化为求解一系列的无约束优化子问题。
PHR算法考虑同时带有等式和不等式约束的优化问题的乘子法
minf（x）
s.t.hi（x）=0，i=1，……，l
gi（x）≥0，i=1，……，m
其基本思想是：
把解等式约束优化问题的乘子法推广到不等式约束优化问题，及先引进辅助变量把不等式约束化为等式约束，然后再利用最优化性条件消去辅助变量。
2.算法步骤
步0选取初始值。
给定x0∈R2,μ1∈Rl,λ1∈Rm,σ1>0,0≤ε<<1,ν∈（0,1）,η>1.令k:
=1.
步1求解子问题。
以xk-1为初始点，求解无约束子问题
minψ（x,μk,λk,σk）,
的极小点xk，其中
ψ（x,μ,λ,σ）=f（x）-
+
+
步2检验终止条件。
若βk≤ε，其中
，则停止迭代，输出xk作为原问题的近似极小点；否则，转步3.
步3更新罚参数。
若
，令
；否则，
。
步4更新乘子向量。
计算
步5令k:
=k+1，转步1.
3.matlab编程
function[x,mu,lambda,output]=multphr（fun,hf,gf,dfun,dhf,dgf,x0）
%功能:
用乘子法解一般约束问题:
minf（x）,s.t.h（x）=0,g（x）>=0
%输入:
x0是初始点,fun,dfun分别是目标函数及其梯度；
%hf,dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；
%gf,dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；
%输出:
x是近似最优点，mu,lambda分别是相应于等式约束和不
%等式约束的乘子向量；output是结构变量，输出近似极小值f,迭
%代次数
maxk=500;%最大迭代次数
sigma=2.0;%罚因子
eta=2.0;theta=0.8;%PHR算法中的实参数
k=0;ink=0;%k,ink分别是外迭代和内迭代次数
epsilon=1e-5;%终止误差值
x=x0;he=feval（hf,x）;gi=feval（gf,x）;
n=length（x）;l=length（he）;m=length（gi）;
%选取乘子向量的初始值
mu=0.1*ones（l,1）;lambda=0.1*ones（m,1）;
btak=10;btaold=10;%用来检验终止条件的两个值
while（btak>epsilon&k%调用BFGS算法程序求解无约束子问题
[x,ival,ik]=bfgs（'mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma）;
ink=ink+ik;
he=feval（hf,x）;gi=feval（gf,x）;
btak=0.0;
for（i=1:
l）,btak=btak+he（i）^2;end
fori=1:
m
temp=min（gi（i）,lambda（i）/sigma）;
btak=btak+temp^2;
end
btak=sqrt（btak）;
ifbtak>epsilon
if（k>=2&btak>theta*btaold）
sigma=eta*sigma;
end
%更新乘子向量
for（i=1:
l）,mu（i）=mu（i）-sigma*he（i）;end
for（i=1:
m）
lambda（i）=max（0.0,lambda（i）-sigma*gi（i））;
end
end
k=k+1;
btaold=btak;
x0=x;
end
f=feval（fun,x）;
output.fval=f;
output.iter=k;
output.inner_iter=ink;
output.bta=btak;
%xstar=[0.5*（sqrt（7）-1）;0.25*（sqrt（7）+1）];
%err1=norm（x-xstar）
4.算例运行情况
利用程序求解约束优化问题
取初始点为x0=（3,3）T，该问题有精确解
.
解：
目标函数文件f1.m，
functionf=f1（x）
f=（x
（1）-2.0）^2+（x
（2）-1.0）^2;
等式约束函数文件h1.m
functionhe=h1（x）
he=x
（1）-2.0*x
（2）+1.0;
不等式约束函数文件g1.m
functiongi=g1（x）
gi=-0.25*x
（1）^2-x
（2）^2+1;
目标函数的梯度文件df1.m，
functiong=df1（x）
g=[2.0*（x
（1）-2.0）,2.0*（x
（2）-1.0）]';
等式约束（向量）函数的Jacobi矩阵（转置）文件dh1.m，
functiondhe=dh1（x）
dhe=[1.0,-2.0]';
不等式约束（向量）函数的Jacobi矩阵（转置）文件dg1.m
functiondgi=dg1（x）
dgi=[-0.5*x
（1）,-2.0*x
（2）]';
在matlab命令窗口输入：
x0=[3,3]';
[x,mu,lambda,output]=multphr（'f1','h1','g1','df1','dh1','dg1',x0）
结果：
x=
0.8229
0.9114
mu=
-1.5945
lambda=
1.8465
output=
fval:
1.3934
iter:
23
inner_iter:
82
bta:
9.5419e-006
5.存在的问题与改进思路
每次迭代都需要更新罚参数和乘子向量，导致运算量变大，但这是为了消除外罚函数法中增广目标函数的病态性质，所以只能尝试引进其他迭代方法
二．凸二次规划的有效集方法编程
1.算法原理
2.算法步骤
步0选取初值。
给定初始可行点
，令k:
=0.
步1解子问题。
确定相应的有效集
。
求解子问题
可以得到极小点dk和拉格朗日乘子向量λk。
若
，转步3；否则，转步2.
步2检验终止准则。
计算拉格朗日乘子
其中
令
若（λk）t≥0，则xk是全局极小点，停算；否则，若（λk）t<0，则令Sk:
=Sk\{t},转步1.
步3确定步长αk。
令
其中
令xk+1:
=xk+αkdk.
步4若αk=1，则令Sk:
=Sk；否则，若αk<1,则令Sk:
=Sk∪{jk},其中jk满足
步5令k:
=k+1,转步1.
3.matlab编程
function[x,lamk,exitflag,output]=qpact（H,c,Ae,be,Ai,bi,x0）
%功能:
用有效集方法解一般约束二次规划问题:
%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）;
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（tm1tm=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;
%%%%%%%%求解子问题%%%%%%%%%%%%%%%
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=zeros（m,1）;
end
4.算例运行情况
利用程序求解下列二次规划问题：
s.t.-2x1-x2≥-3,
x1-x2≥-1,
-x1-2x2≥-2,
x1≥0,x2≥0.
该问题有精确解
最优值f（x*）=
.
解：
在matlab命令窗口输入下列命令：
>>H=[1-1;-12];
>>c=[-6-2]';
>>Ai=[-2-1;1-1;-1-2;10;01];
>>bi=[-3-1-200]';
>>x0=[00]';
>>[x,lambda,exitflag,output]=qpact（H,c,[],[],Ai,bi,x0）
运行结果：
x=
1.3333
0.3333
lambda=
2.4444
0.1111
exitflag=
0
output=
fval:
-8.1111
iter:
7
5.存在的问题与改进思路
问题：
有效集算法是一个非常有效的求解凸二次规划的方法，美中不足的是在调整Jk的过程中，当
时，则以
作为搜索方向，虽然保证了dk为下降方向，但不能保证为最速下降方向。
改进思路：
可以将投影梯度算法的方法引进来，将
在xk点的投影方向作为搜索方向
三．某一算法的改进
有效集算法是一个非常有效的求解凸二次规划的方法，美中不足的是在调整Jk的过程中，当
时，则以
作为搜索方向，虽然保证了dk为下降方向，但不能保证为最速下降方向。
可以将投影梯度算法的方法引进来，将
在xk点的投影方向作为搜索方向.从而对有效集作如下修改：
设xk是问题
的一个可行解，且有效指标集为Jk（xk）=E∪I（xk）,求解相应的等式约束问题
假设所得的解为
且相应的lagrange乘子向量λk，首先检验
（此算法改进摘抄于豆丁网论文>毕业论文>二次规划的改进有效集算法，，，向老师忏悔，没有学透，无法自主开发有效的算法改进）
升级会员 