优化方法编程大作业.pdf
优化方法编程大作业1 1：Matlab 代码：将下面代码输入 command 中即可：clcclearsyms x2 x1 fa=0;b=+inf;c1=0.1;c2=0.5;%rk 为搜索步长，xk 为变量，sk 为搜索方向,gk 为梯度值。rk=1;sk=-1;1;xk0=0;1;%xk0 初始值f=(x2-x12)2+(1-x1)2;%所求函数j=0;%循环的标志变量%循环求解最佳的步长while j=0;fk0=subs(f,x1,x2,xk0(1),xk0(2);%求 f 在 xk0 处的函数值 gk0=subs(jacobian(f,x1,x2),x1,x2,xk0(1),xk0(2);%求 f 在 xk0 处的梯度值 xk=xk0+rk*sk;fk=subs(f,x1,x2,xk(1),xk(2);%求 f 在 xk 处的函数值 gk=subs(jacobian(f,x1,x2),x1,x2,xk(1),xk(2);%求 f 在 xk 处的梯度值%搜索步长满足 wolfe-powell 第一个准则 if fk0-fk=-c1*rk*gk0*sk%搜索步长满足 wolfe-powell 第二个准则 if gk*sk=c2*gk0*sk%输出最佳的步长 rk=rk;return;%搜索步长不满足 wolfe-powell 准则，继续迭代 else j=j+1;a=rk;rk=min(2*rk,0.5*(rk+b);end else j=j+1;b=rk;rk=0.5*(a+rk);endendxkgkrk运行结果：j=53，迭代 54 次。xk=-0.0000 1.0000gk=-2.0000 2.0000rk=1.1102e-162首先定义三个函数 f（求函数值），g（求梯度值）和 wolfe（第一题的算法求步长）并存为 M 文件放在一个文件夹下：function f=f(x)f=x(1)2-2*x(1)*x(2)+2*x(2)2+x(3)2+x(4)2-x(2)*x(3)+2*x(1)+3*x(2)-x(3);endfunction g=g(yk0)syms y1 y2 y3 y4;y=y12-2*y1*y2+2*y22+y32+y42-y2*y3+2*y1+3*y2-y3;g=subs(jacobian(y,y1,y2,y3,y4),y1,y2,y3,y4,yk0(1),yk0(2),yk0(3),yk0(4);%求 gf 在 yk0 处的梯度值endfunction tk=wolfe(yk0,dk)a=0;b=+inf;c1=0.1;c2=0.5;%rk 为搜索步长，xk 为变量，sk 为搜索方向,gk 为梯度值。tk=1;%sk=-1;1;%xk0=0;1;%xk0 初始值%f=(x2-x12)2+(1-x1)2;%所求函数j=0;%循环的标志变量%循环求解最佳的步长while j=0;fk0=f(yk0);%求 f 在 xk0 处的函数值 gk0=g(yk0);%求 f 在 xk0 处的梯度值 yk=yk0+tk*dk;fk=f(yk);%求 f 在 xk 处的函数值 gk=g(yk);%求 f 在 xk 处的梯度值%搜索步长满足 wolfe-powell 第一个准则 if fk0-fk=-c1*tk*gk0*dk%搜索步长满足 wolfe-powell 第二个准则 if gk*dk=c2*gk0*dk%输出最佳的步长 tk=tk;return;%搜索步长不满足 wolfe-powell 准则，继续迭代 else j=j+1;a=tk;tk=min(2*tk,0.5*(tk+b);end else j=j+1;b=tk;tk=0.5*(a+tk);endendend将下面的代码输入 Matlab 的 command 中可得结果：xk0=0;0;0;0;k=0;n=4;gk0=g(xk0);%求 f 在 xk0 处的梯度值sk=-gk0;%sk 为 skwhile k=n-1 gk0=g(xk0);%求 f 在 xk0 处的梯度值 rk=wolfe(xk0,sk);xk=xk0+rk*sk;fk=f(xk);%求 f 在 xk 处的函数值 gk=g(xk);%求 f 在 xk 处的梯度值 xk0=xk;if(gk(1)2+gk(2)2+gk(3)2+gk(4)2)(1/2)=0;fk0=f(yk0);%求 f 在 xk0 处的函数值 gk0=g(yk0);%求 f 在 xk0 处的梯度值 yk=yk0+tk*dk;fk=f(yk);%求 f 在 xk 处的函数值 gk=g(yk);%求 f 在 xk 处的梯度值%搜索步长满足 wolfe-powell 第一个准则 if fk0-fk=-c1*tk*gk0*dk%搜索步长满足 wolfe-powell 第二个准则 if gk*dk=c2*gk0*dk%输出最佳的步长 tk=tk;return;%搜索步长不满足 wolfe-powell 准则，继续迭代 else j=j+1;a=tk;tk=min(2*tk,0.5*(tk+b);end else j=j+1;b=tk;tk=0.5*(a+tk);endendend（1）编写最速下降法函数minFD。function x,minf,k=minFD(x0)%功能：用变尺度求解无约束优化问题%输入：初始点 xk%输出：最优解 val，及最小点 x，迭代次数 k eps=1.0e-4;tol=1;k=0;while toleps v=-gfun(x0);tol=norm(v);step=buchang(x0,v);%一维搜索 xk=x0+step*v;x0=xk;k=k+1;endx=xk;minf=fun(x);在 command 中输入 xk=1,0;x,fv,k=minFD(xk)运行结果为：x=-0.4194 0fv=0.7729k=12（2）编写牛顿法求解极值minNT。function x,minf,k=minNT(x0)%功能：用变尺度求解无约束优化问题%输入：初始点 xk%输出：最优解 val，及最小点，迭代次数eps=1.0e-4;tol=1;k=0;while toleps v =gfun(x0);pv=Hess(x0);p=-inv(pv)*v;tol=norm(p);step=buchang(x0,p);%一维搜索 x1=x0+step*p;x0=x1;k=k+1;endx=x1;minf=fun(x);在 command 中输入 xk=1,0;x,fv,k=minNT(xk)运行结果为：x=-0.4194 0fv=0.7729k=6（3）编写 BFGS 法求解极值function x,minf,m=minBFGS(x0)%功能：BFGS 法求解无约束优化问题%输入：初始点 xk%输出：最优解 minf，及最小点 x，迭代次数 meps=1.0e-4;H=1,0;0,1;k=0;m=0;n=2;while 1 v =gfun(x0);p=-H*v;step=buchang(x0,p);%一维搜索 x1=x0+step*p;vk=gfun(x1);if norm(vk)xk=1,0;x,fv,k=minBFGS(xk)运行结果为：x=-0.4194 0fv=0.7729k=54首先定义拉格朗日法求解等式约束的函数QuadLagR。function xv,fv=QuadLagR(H,c,A,b)%正定矩阵:H%二次规划一次项系数向量：c%等式约束右端向量 b%目标函数取极小值时的自变量值：xv%目标函数极小值：fvinvH=inv(H);F=invH*transpose(A)*inv(A*invH*transpose(A)*A*invH-invH;D=inv(A*invH*transpose(A)*A*invH;xv=F*c+transpose(D)*b;fv=transpose(xv)*H*xv/2+transpose(c)*xv;有效集法求解函数：function xv,fv=ActivdeSet(H,c,A,b,x0,AcSet)%正定矩阵:H%二次规划一次项系数向量：c%等式约束右端向量 b%初始点：x0%初始约束集：AcSset%目标函数取极小值时的自变量值：xv%目标函数极小值：fvsz=size(A);xx=1:sz(1);nonAcSet=zeros(1,1);m=1;for i=1:sz(1)%寻找非约束指标集 if(isempty(find(AcSet=xx(i),1)nonAcSet(m)=i;m=m+1;else ;endendinvH=inv(H);bCon=1;while bCon cl=H*x0+c;Al=A(AcSet,:);bl=b(AcSet);xm=QuadLagR(H,cl,Al,zeros(length(AcSet),1);%用拉格朗日求解由约束指标集定义的等式约束二次规划 if xm=0%最优解为 0 trAl=transpose(Al);R=inv(Al*invH*trAl)*Al*invH;lamda=R*(H*xm+cl);%拉格朗日乘子 lmin,index=min(lamda);if lmin=0 xv=x0;%停止计算，得出最优解 fv=transpose(x0)*H*x0/2+transpose(c)*x0;bCon=0;return;else l=length(AcSet);%删除指标 AcSet（index）nonAcSet=nonAcSet AcSet(index);sa,sb=sort(nonAcSet);nonAcSet=sa;tmpAcS=AcSet(1:(index-1)AcSet(index+1):l);AcSet=tmpAcS;end else d=xm;mAlpha=inf;adinAcS=0;for i=1:length(nonAcSet)%确定 alpha if A(nonAcSet(i),:)*d 0 alpha=(b(nonAcSet(i)A(nonAcSet(i),:)*x0)/(A(nonAcSet(i),:)*d);if alpha mAlpha mAlpha=alpha;adinAcS=nonAcSet(i);end end end mXA=min(1,mAlpha);%确定 alpha-k x0=x0+mXA*d;if mXA=0-xv=x0;fv=transpose(x0)*H*x0/2+transpose(c)*x0;bCon=0;return;else l=length(AcSet);nonAcSet=nonAcSet AcSet(index);sa,sb=sort(nonAcSet);nonAcSet=sa;tmpAcS=AcSet(1:(index-1)AcSet(index+1):l);AcSet=tmpAcS;end end endend在 command 中输入 xv,fv=ActivdeSet(H,c,A,b,1;1,)运行结果：xv=0.8000 1.2000fv=-7.20005首先定义三组函数:1）定义目标函数及梯度函数function f=f1(x)f=4*x(1)-x(2)2-12;endfunction g=df1(x)g=4,-2.0*x(2);end2）定义约束条件函数及梯度函数约束：function he=h1(x)he=25-x(1)2-x(2)2;endfunction gi=g1(x)gi=10*x(1)-x(1)2-x(2)2+10*x(2)-34;end3）梯度：function dhe=dh1(x)dhe=-2*x(1),-2.0*x(2);endfunction dgi=dg1(x)dgi=10-2*x(1),10-2.0*x(2);end乘子法函数：function x,f,ink=multphr(fun,hf,gf,dfun,dhf,dgf,x0)%功能:用乘子法解一般约束问题:min f(x),s.t.h(x)=0,g(x)=0%输入:x0 是初始点,fun,dfun 分别是目标函数及其梯度；%hf,dhf 分别是等式约束（向量）函数及其Jacobi 矩阵的转置；%gf,dgf 分别是不等式约束（向量）函数及其Jacobi 矩阵的转置；%输出:x 是近似最优点，mu,lambda 分别是相应于等式约束和不%等式约束的乘子向量；output 是结构变量，输出近似极小值f,迭%代次数maxk=5000;%最大迭代次数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(btakepsilon&kepsilon 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;endf=feval(fun,x);command 中输入 x0=3,4;xv,fv,k=multphr(f1,h1,g1,df1,dh1,dg1,x0)运行结果：xv=1.0013 4.8987fv=-31.9923k=416 6目标函数：目标函数：function f=fun5(x)f=4*x(1)-x(2)2-12;end目标函数的梯度：目标函数的梯度：function g=dfun5(x)g=4;-2*x(2);end等式约束函数：等式约束函数：function he=h5(x)he=25-x(1)2-x(2)2;end等式约束函数的梯度：等式约束函数的梯度：function dhe=dh5(x)dhe=-2*x(1);-2*x(2);end第一个不等式约束的函数：第一个不等式约束的函数：function gi=g5_1(x)gi=10*x(1)-x(1)2+10*x(2)-x(2)2-34;end第一个不等式约束的函数的梯度：第一个不等式约束的函数的梯度：function dgi=dg5_1(x)dgi=10-2*x(1);10-2*x(2);end第二个不等式约束的函数：第二个不等式约束的函数：function gi=g5_2(x)gi=x(1);end第二个不等式约束的函数的梯度：第二个不等式约束的函数的梯度：function dgi=dg5_2(x)dgi=1;0;end第三个不等式约束的函数：第三个不等式约束的函数：function gi=g5_3(x)gi=x(2);end第三个不等式约束的函数的梯度：第三个不等式约束的函数的梯度：function dgi=dg5_3(x)dgi=0;1;end增广朗格朗日函数：增广朗格朗日函数：function psi=mpsi(x,fun5,h5,g5_1,g5_2,g5_3,mu,lambda1,lambda2,lambda3,sigma)%增广拉格朗日函数f=feval(fun5,x);he=feval(h5,x);g1=feval(g5_1,x);g2=feval(g5_2,x);g3=feval(g5_3,x);psi=f;s1=0.0;psi=psi+he*mu;s1=s1+he2;psi=psi+0.5*sigma*s1;s2=0.0;s3_1=min(0.0,lambda1+sigma*g1);s3_2=min(0.0,lambda2+sigma*g2);s3_3=min(0.0,lambda3+sigma*g3);s2=s2+(s3_12-lambda12)+(s3_22-lambda22)+(s3_32-lambda32);psi=psi+s2/(2.0*sigma);endpsi=psi+s2/(2.0*sigma);增广朗格朗日函数梯度：增广朗格朗日函数梯度：functiondpsi=dmpsi(x0,h5,g5_1,g5_2,g5_3,dfun5,dh5,dg5_1,dg5_2,dg5_3,mu,lambda1,lambda2,lambda3,sigma)%增广朗格朗日函数梯度dpsi=feval(dfun5,x0);he=feval(h5,x0);g1=feval(g5_1,x0);g2=feval(g5_2,x0);g3=feval(g5_3,x0);dhe=feval(dh5,x0);dg1=feval(dg5_1,x0);dg2=feval(dg5_2,x0);dg3=feval(dg5_3,x0);s3_1=min(0.0,lambda1+sigma*g1);s3_2=min(0.0,lambda2+sigma*g2);s3_3=min(0.0,lambda3+sigma*g3);dpsi=dpsi+(sigma*he+mu)*dhe;dpsi=dpsi+s3_1*dg1+s3_2*dg2+s3_3*dg3;endDFPDFP求步长：求步长：functionx,val,k=dfp(fun,gfun,x0,fun5,h5,g5_1,g5_2,g5_3,dfun5,dh5,dg5_1,dg5_2,dg5_3,mu,lambda1,lambda2,lambda3,sigma)maxk=500;%最大迭代次数rho=0.55;sigma1=0.4;epsilon1=1.0e-5;k=0;n=length(x0);Bk=eye(n);%Bk=feval(Hess,x0);while(kmaxk)gk=feval(gfun,x0,h5,g5_1,g5_2,g5_3,dfun5,dh5,dg5_1,dg5_2,dg5_3,mu,lambda1,lambda2,lambda3,sigma);%计算梯度 if(norm(gk)epsilon1)break;end%检验终止准则 dk=-Bk*gk;%解方程组，计算搜索方向 m=0;mk=0;while(m30)%用Armijo搜索步长newf=feval(fun,x0+rhom*dk,fun5,h5,g5_1,g5_2,g5_3,mu,lambda1,lambda2,lambda3,sigma);oldf=feval(fun,x0,fun5,h5,g5_1,g5_2,g5_3,mu,lambda1,lambda2,lambda3,sigma);if(newf0)Bk=Bk-(Bk*(yk*yk)*Bk)/(yk*Bk*yk)+(sk*sk)/(sk*yk);end k=k+1;x0=x;endval=feval(fun,x0,fun5,h5,g5_1,g5_2,g5_3,mu,lambda1,lambda2,lambda3,sigma);主程序：主程序：functionx,mu,lambda1,lambda2,lambda3,output=multphr(fun5,h5,g5_1,g5_2,g5_3,dfun5,dh5,dg5_1,dg5_2,dg5_3,x0)%x为最优解，mu为最终惩罚因子，lambda1、lambda2、lambda3分别为三个不等式的惩罚因子%x0为初值maxk=10000;%最大迭代次数sigma=2.0;%罚因子eta=2.0;theta=0.8;%PHR算法中的实参k=0;ink=0;%k,ink分别是外迭代法和内迭代次数epsilon=1e-5;%终止误差x=x0;he=feval(h5,x);g1=feval(g5_1,x);%选取乘子向量的初始值mu=0.1*ones(1,1);lambda1=0.1*ones(1,1);lambda2=0.1*ones(1,1);lambda3=0.1*ones(1,1);btak=10;btaold=10;%用来检验终止条件的两个值while(btakepsilon&kepsilon if(k=2&btak theta*btaold)sigma=eta*sigma;end%更新乘子问题 mu=mu+sigma*he;lambda1=min(0.0,lambda1+sigma*g1);lambda2=min(0.0,lambda2+sigma*g2);lambda3=min(0.0,lambda3+sigma*g3);end k=k+1;btaold=btak;x0=x;endf=feval(fun5,x);output.fval=f;%最优解output.iter=k;%外迭代次数output.inner_iter=ink;%内迭代次数output.bta=btak;%终止检验数计算结果：计算结果：输入：输入：x0=0;0;x,mu,lambda1,lambda2,lambda3,output=multphr(fun5,h5,g5_1,g5_2,g5_3,dfun5,dh5,dg5_1,dg5_2,dg5_3,x0)结果为：结果为：最优点：最优点：x=1.0013 4.8987等式罚因子等式罚因子：mu=-1.0156不等式罚因子不等式罚因子：lambda1=-0.7545lambda2=0lambda3=0output=最优值最优值 fval:-31.9923外迭代次数：外迭代次数：iter:5内迭代次数：内迭代次数：inner_iter:27终止检验数终止检验数：bta:2.8300e-0077 7在 command 中输入下面代码：a=2 1 1;1-1-1;b=2-1;c=3 1 2;x=linprog(c,a,b)运行结果：x=T 1.0e+12*(-0.7993,1.7216,-2.5207)
