二次规划

2012-2013（1）专业课程实践论文二次规划王霄，0818180103，R数学08-1班于灏，0818180105，R数学08-1班夏之秋，0818180110，R数学08-1班一、算法理论二次规划是非线性优化中的一种特殊情形,它的目标函数是二次实函数,约束函数都是线性函数.由于二次规划比较简单,便于求解(仅次于线性规划),并且一些非线性优化问题可以转化为求解一系列的二次规划问题,因此二次规划的求解方法较早引起人们的重视,成为求解非线性优化的一个重要途径.二次规划的算法较多,本论文仅介绍求解一般约束凸二次规划的有效集方法.考虑一般二次规划},,,1{,0},,,1{,0..,21minmlIibxalEibxatsxcHxxiTiiTiTT有效集方法的最大难点是事先一般不知道有效集*xS,因此只有想办法构造一个集合序列去逼近它.即从初始点0x出发,计算有效集0xS,解对应的等式约束子问题.重复这一做法,得到有效集序列kxS,,1,0k,使之*xSxSk,以获得原问题的最优解.我们分六步来介绍有效集方法的算法原理和实施步骤.第一步选定初始值.给定初始可行点nRx0，令0k。第二步求解子问题.确定相应的有效集kkxIES，求解子问题kTiTkTkSidatsdgHdddq,0..21min得到极小点kd和拉格朗日乘子向量k。若0kd转入步三；否则转步二。第三步检验终止准则.计算拉格朗日乘子kkkgB其中kikkkkkkkSiaAHAATHABcHxg,,,111，令iktkmin。若0tk，则kx是全局极小点，停算。否则若0tk，则令tSSkk\，转步一。第四步确定步长k.令,1mink，其中kTikTiikdaxabmin，其中0kTida。令kkkkdxx1。第五步若1k，则令kkSS1，否则，若1k，则令kkkjSS1，其中kj满足kTjkTjjkdaxabkkk。第六步令1kk，转步一。二、算法框图开始选定初始值，给定初始可行点，令0k确定相应的有效集kkxIES得到极小点kd和拉格朗日乘子向量k若0kd计算拉格朗日乘子kkkgB，令iktkmin若0tk令tSSkk\kx是全局极小点结束三、算法程序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,:)*xbi(i)+epsilon),index(i)=0;endend%算法主程序while(k=kmax)%求解子问题Aee=[];if(ne0),Aee=Ae;endfor(j=1:ni)if(index(j)0),Aee=[Aee;Ai(j,:)];endendgk=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)));endif(y=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk)index(i)=0;break;endendendk=k+1;elseexitflag=1;%求步长alpha=1.0;tm=1.0;for(i=1:ni)if((index(i)==0)&(Ai(i,:)*dk0))tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if(tm1tm)tm=tm1;ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;%修正有效集if(tm1),index(ti)=1;endendif(exitflag==0),break;endk=k+1;endoutput.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(m0)rb=Ae*ginvH*c+be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;end四、算法实现例1.求解二次规划问题,0,0,22,1,32..,2621min2121212121222121xxxxxxxxtsxxxxxxxf解：在Matlab命令窗口输入下列命令:H=[1-1;-12];c=[-6-2]';Ae=[];be=[];Ai=[-2-1;1-1;-1-2;10;01];bi=[-3-1-200]';x0=[00]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)的计算结果为x=1.33330.3333lambda=2.44440.1111exitflag=0output=fval:-8.1111iter:7该问题有精确解Tx31,34*，最优值918*xf。例2.求解二次规划问题.0,0,623..,102min212121222121xxxxtsxxxxxxxf解：在Matlab命令窗口输入下列命令:H=[2-1;-14];c=[-1-10]';Ae=[];be=[];Ai=[-3-2;10;01];bi=[-600]';x0=[00]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)的计算结果为x=0.50002.2500lambda=0.7500exitflag=0output=fval:-13.7500iter:8