优化方法上机作业--2013

第一题（牛顿法和不精确一维搜索）牛顿法：syms x1 x2;f=(x2-(x1)*(x1))^2+(1-x1)^2;v=[x1,x2];df=jacobian(f,v);df=df.';G=jacobian(df,v);epson=1e-12;x0=[0,1]';g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});G1=subs(G,{x1,x2},{x 0(1,1),x0(2,1)});k=0;mul_count=0;sum_count=0;mul_count=mul_count+12;sum_count=sum_count+6;while(norm(g1)>epson)p=-G1\g1;x0=x0+p;g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});k=k+1;mul_count=mul_count+16;sum_count=sum_count+11;end;kx0mul_countsum_count运行结果：k =9x0 =11mul_count =156sum_count =105不精确一维搜索：fun1.m文件function f=fun1(x)f=(x(1)-1)^2+(x(2)-x(1)^2)^2;gfun1.m文件function gf = gfun1(x)gf=[2*x(1)*(x(1)^2-x(2))+2*(x(1)-1),2*(x(2)-x(1)^2)]';wolfepowell.m文件function [k,m,opt,x]=wolfepowell(xk,sk)max=1000;c1=0.1; c2=0.6;a=0;b=inf;dk=0.2;m=0;while(m<=max)if (fun1(xk)-fun1(xk+sk*dk)<-c1*dk*gfun1(xk)'*sk)b=dk;dk=(dk+a)/2;elseif(fun1(xk)-fun1(xk+sk*dk)>=-c1*dk*gfun1(xk)'*sk)&&(gfun1(xk+sk*dk)'*sk<c2*gf un1(xk)'*sk);a=dk;if 2*dk<(dk+b)/2dk=2*dk;elsedk=(dk+b)/2;endelsemk=m;break;endm=m+1;endm=mk;k=dk;x=xk+sk*dk;opt=fun1(x);运行结果：>>xk=[0,1]';sk=[1,-1]';[k,m,opt,x]= wolfepowell(xk,sk)k =0.8000m =2opt =1.0080x =0.80000.2000第二题（共轭梯度法）frcg.m文件function [x,val,k]=frcg(fun,gfun,x0)maxk=5000;rho=0.6;sigma=0.4;k=0; epsilon=1e-4;n=length(x0);while(k<maxk)g=feval(gfun,x0);itern=k-(n+1)*floor(k/(n+1));itern=itern+1;if(itern==1)d=-g;elsebeta=(g'*g)/(g0'*g0);d=-g+beta*d0; gd=g'*d;if(gd>=0.0)d=-g;endendif(norm(g)<epsilon), break; endm=0; mk=0;while(m<20)if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)mk=m; break;endm=m+1;endx0=x0+rho^mk*d;val=feval(fun,x0);g0=g; d0=d;k=k+1;endx=x0;val=feval(fun,x);fun.m文件function f=fun(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); gfun.m文件function gf=gfun(x)gf=[2*x(1)-2*x(2)+2, -2*x(1)+4*x(2)+3,2*x(3)-x(2)-1,2*x(4)]';运行结果：>>x0=[0 0 0 0]';[x,val,k]=frcg('fun','gfun',x0)x =-3.4064-2.6515-0.7032val =-7.8338k =5000第三题（BFGS法、最速下降法和牛顿法）fun3.m文件function f=fun3(x)f=(x(1)-1)^2+5*(x(2)-x(1)^2)^2;gfun3.m文件function gf = gfun3(x)gf=[20*x(1)*(x(1)^2-x(2))+2*(x(1)-1),10*(x(2)-x(1)^2)]'; Hessen.m文件function Hess=Hessen(x)Hess=[60*x(1)^2-20*x(2)+2,(-20)*x(1);(-20)*x(1),10];Steepest.m文件function [m,opt,x] = Steepest(xk)m=0;Eps=1.0e-4;max=1000;while(m<max)g=gfun3(xk);sk=-g;dk=wolfepowell(xk,sk);x=xk+dk*sk;g0=gfun3(x);if(norm(g0)<=Eps)break;elsexk=x;endm=m+1;endmk=m;x=xk;opt=fun3(xk);newton.m文件function [x,opt,mk] = newton(xk)m=0;Eps=1e-4;max=100;H=[1,0;0,1];while(m<max)g=gfun3(xk);M=Hess(xk);H=inv(M);sk=-H*g;dk=wolfepowell(xk,sk);x=xk+dk*sk;g0=gfun3(x);if(norm(g0)<=Eps)break;elsexk=x;endm=m+1;endmk=m;x=xk;opt=fun3(xk);bfgs.m文件function [x,opt,mk] = bfgs(xk)m=0;Eps=1e-4;max=100;H=[1,0;0,1];while(m<max)g=gfun3(xk);sk=-H*g;dk=wolfepowell(xk,sk);x=xk+dk*sk;g0=gfun3(x);if(norm(g0)<=Eps)break;elseH=H+((1+((g0-g)'*H*(g0-g))/((x-xk)'*(g0-g)))*(x-xk)*(x-xk)'-H*(g0-g)*(x-xk)'-(x-x k)*(g0-g)'*H)/((x-xk)'*(g0-g));xk=x;endm=m+1;endmk=m;x=xk;opt=fun3(xk);wolfepowell.m文件function [k,m,opt,x]=wolfepowell(xk,sk)max=1000;c1=0.1; c2=0.6;a=0;b=inf;dk=0.2;m=0;while(m<=max)if (fun3(xk)-fun3(xk+sk*dk)<-c1*dk*gfun3(xk)'*sk)b=dk;dk=(dk+a)/2;elseif(fun3(xk)-fun3(xk+sk*dk)>=-c1*dk*gfun3(xk)'*sk)&&(gfun3(xk+sk*dk)'*sk<c2*gf un3(xk)'*sk);a=dk;if 2*dk<(dk+b)/2dk=2*dk;elsedk=(dk+b)/2;endelsemk=m;break;endm=m+1;endm=mk;k=dk;x=xk+sk*dk;opt=fun3(x);输入参数：牛顿法运行结果：>> xk=[2,0]';[x,opt,mk] = newton(xk)x =1.00001.0000opt =3.2601e-10mk =12最速下降法运行结果：>>xk=[2,0]';[x,opt,mk] = Steepest(xk)x =0.99990.9998opt =8.2108e-09mk =354BFGS法运行结果：>>xk=[2,0]';[x,opt,mk] = bfgs(xk)x =0.99990.9999opt =4.7725e-09mk =16(0.9999, 0.9998) 最速下降法（迭代次数354）即得数值最优点为x= (1.0000, 1.0000) 牛顿法（迭代次数12）(0.9999, 0.9999) BFGS公式（迭代次数16）最优解近似为0第四题（乘子法）multphr.m文件function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)maxk=500;sigma=2.0;eta=2.0; theta=0.8;k=0; ink=0;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<maxk)[x,~,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; endfor (i=1:m);temp=min(gi(i),lambda(i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);if(btak>epsilon)if(k>=2&&btak> theta*btaold)sigma=eta*sigma;endfor (i=1:l), mu(i)=mu(i)-sigma*he(i); endfor (i=1:m)lambda(i)=max(0.0,lambda(i)-sigma*gi(i));endendk=k+1;btaold=btak;x0=x;endf=feval(fun,x);output.fval=f;output.iter=k;output.inner_iter=ink;output.bta=btak;mpsi.m文件function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x);l=length(he); m=length(gi);psi=f; s1=0.0;for(i=1:l)psi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*sigma*s1;s2=0.0;for(i=1:m)s3=max(0.0, lambda(i) - sigma*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2.0*sigma);h1.m文件function he=h1(x)he=-x(1)^2-x(2)^2+25.0;f1.m文件function f=f1(x)f=4*x(1)-x(2)^2-12;g1.m文件function gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;h1.m文件function dhe = dh1(x)dhe = [-1*x(1), -1*x(2)]';dg1.m文件function dgi = dg1(x)dgi = [10-2*x(1), 10-2*x(2)]';df1.m文件function g=df1(x)g = [4, -2.0*x(2)]';bfgs.m文件function [x,val,k]=bfgs(fun,gfun,x0,varargin)maxk=500;rho=0.55; sigma1=0.4; epsilon1=1e-5;k=0; n=length(x0);Bk=eye(n);while(k<maxk)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon1), break; enddk=-Bk\gk;m=0; mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma1*rho^m*gk'*dk)mk=m; break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0; yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1; x0=x;endval=feval(fun,x0,varargin{:});dmpsi.m文件function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) dpsi=feval(dfun,x);he=feval(hf,x); gi=feval(gf,x);dhe=feval(dhf,x); dgi=feval(dgf,x);l=length(he); m=length(gi);for(i=1:l)dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);endfor(i=1:m)dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);end运行结果：>>x0=[1,1]';[x,mu,lambda,output]=multphr('f1','h1','g1','df1','dh1','dg1',x0) x =1.00134.8987mu =2.0312lambda =0.7545output =fval: -31.9923iter: 5inner_iter: 58bta: 4.3187e-07第五题（有效集法）qpact.m文件function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) 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; endendwhile (k<=kmax)Aee=[ ];if(ne>0), 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,:)*dk<0))tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk); if(tm1<tm)tm=tm1; ti=i;endendendalpha=min(alpha,tm);x = x+alpha*dk;if(tm<1), index(ti)=1; endendif(exitflag==0), break; endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;qsubp.m文件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);elsex = -ginvH*c;lambda = zeros(m,1);endcallqpact.m文件function callqpactH=[2 0; 0 2];c=[-2 -5]';Ae=[ ]; be=[ ];Ai=[1 -2; -1 -2; -1 2;1 0;0 1];bi=[-2 -6 -2 0 0]';x0=[0 0]';[x, lambda, exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)运行结果：>>callqpactx =1.40001.7000lambda =0.8000exitflag =output =fval: -6.4500iter: 7第七题fun文件function f=fun(x)f=abs(x(1))+abs(x(2))+abs(x(3))+abs(x(4));主程序x0=[1;1;1;1];A=[]; b=[];Aeq=[1 0 0 0.5;0 1 0.2 0.3;0 0.1 1 0.2];beq=[-1;0.2;1]; VLB=[]; VUB=[];[x,fval]=fmincon('fun',x0,A,b,Aeq,beq,VLB,VUB)运行结果x =-1.0000-0.00001.00000.0000 fval =2.0000。