大连理工大学优化方法上机作业

1、大连理工大学优化方法上机大作业优化方法上机大作业 学 院：电子信息与电气工程学部 姓 名： 学 号： 指导老师： 上机大作业（一）%目标函数function f=fun(x)f=100*(x(2)-x(1)2)2+(1-x(1)2;end%目标函数梯度function gf=gfun(x)gf=-400*x(1)*(x(2)-x(1)2)-2*(1-x(1);200*(x(2)-x(1)2);End%目标函数Hess矩阵function He=Hess(x)He=1200*x(1)2-400*x(2)+2,-400*x(1); -400*x(1), 200;end%线搜索步长function

2、mk=armijo(xk,dk)beta=0.5; sigma=0.2;m=0; maxm=20;while (m<=maxm) if(fun(xk+betam*dk)<=fun(xk)+sigma*betam*gfun(xk)'*dk) mk=m; break; end m=m+1;endalpha=betamknewxk=xk+alpha*dkfk=fun(xk)newfk=fun(newxk)%最速下降法function k,x,val=grad(fun,gfun,x0,epsilon)%功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标

3、函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值maxk=5000; %最大迭代次数beta=0.5; sigma=0.4;k=0;while(k<maxk) gk=feval(gfun,x0); %计算梯度 dk=-gk; %计算搜索方向 if(norm(gk)<epsilon), break;end %检验终止准则 m=0;mk=0; while(m<20) %用Armijo搜索步长 if(feval(fun,x0+betam*dk)<=feval(fun,x0)+sigma*betam*gk'

4、*dk) mk=m;break; end m=m+1; end x0=x0+betamk*dk; k=k+1;endx=x0;val=feval(fun,x0);>> x0=0;0;>> k,x,val=grad('fun','gfun',x0,1e-4)迭代次数：k = 1033x = 0.9999 0.9998val = 1.2390e-008%牛顿法x0=0;0;ep=1e-4;maxk=10;k=0;while(k<maxk) gk=gfun(x0); if(norm(gk)<ep) x=x0 miny=fun(x)

5、k0=k break; else H=inv(Hess(x0); x0=x0-H*gk; k=k+1; endendx = 1.0000 1.0000miny = 4.9304e-030迭代次数k0 = 2%BFGS方法function k,x,val=bfgs(fun,gfun,x0,varargin)%功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值N=1000;epsilon=1e-4;beta=0.55;sigma=0.4;n=length(x

6、0);Bk=eye(n);k=0;while(k<N) gk=feval(gfun,x0,varargin:); if(norm(gk)<epsilon), break;end dk=-Bkgk; m=0;mk=0; while(m<20) newf=feval(fun,x0+betam*dk,varargin:); oldf=feval(fun,x0,varargin:); if(newf<=oldf+sigma*betam*gk'*dk) mk=m;break; end m=m+1; end x=x0+betamk*dk; sk=x-x0; yk=feval

7、(gfun,x,varargin:)-gk; if(yk'*sk>0) Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk); end k=k+1; x0=x;endval=feval(fun,x0,varargin:);>> x0=0;0;>> k,x,val=bfgs('fun','gfun',x0)k = 20x = 1.0000 1.0000val = 2.2005e-011%共轭梯度法function k,x,val=frcg(fun,

8、gfun,x0,epsilon,N)if nargin<5,N=1000;endif nargin<4, epsilon=1e-4;endbeta=0.6;sigma=0.4;n=length(x0);k=0;while(k<N) gk=feval(gfun,x0); itern=k-(n+1)*floor(k/(n+1); itern=itern+1; if(itern=1) dk=-gk; else betak=(gk'*gk)/(g0'*g0); dk=-gk+betak*d0; gd=gk'*dk; if(gd>=0),dk=-gk;en

9、d end if(norm(gk)<epsilon),break;end m=0;mk=0; while(m<20) if(feval(fun,x0+betam*dk)<=feval(fun,x0)+sigma*betam*gk'*dk) mk=m;break; end m=m+1; end x=x0+betam*dk; g0=gk; d0=dk; x0=x;k=k+1;endval=feval(fun,x); >> x0=0;0;k,x,val=frcg('fun','gfun',x0,1e-4,1000)k = 122x

10、 = 1.0001 1.0002val = 7.2372e-009上机大作业（二）%目标函数function f_x=fun(x)f_x=4*x(1)-x(2)2-12;%等式约束条件function he=hf(x)he=25-x(1)2-x(2)2;end%不等式约束条件function gi_x=gi(x,i)switch i case 1 gi_x=10*x(1)-x(1)2+10*x(2)-x(2)2-34; case 2 gi_x=x(1); case 3 gi_x=x(2); otherwiseend%求目标函数的梯度function L_grad=grad(x,lambda,c

11、igma)d_f=4;2*x(2);d_g(:,1)=-2*x(1);-2*x(2);d_g(:,2)=10-2*x(1);10-2*x(2);d_g(:,3)=1;0;d_g(:,4)=0;1;L_grad=d_f+(lambda(1)+cigma*hf(x)*d_g(:,1);for i=1:3 if lambda(i+1)+cigma*gi(x,i)<0 L_grad=L_grad+(lambda(i+1)+cigma*gi(x,i)*d_g(:,i+1); continue endend%增广拉格朗日函数function LA=lag(x,lambda,cee)LA=fun(x)

12、+lambda(1)*hf(x)+0.5*cee*hf(x)2;for i=1:3 LA=LA+1/(2*cee)*(min(0,lambda(i+1)+cee*gi(x,i)2-lambda(i+1)2);endfunction xk=BFGS(x0,eps,lambda,cigma)gk=grad(x0,lambda,cigma);res_B=norm(gk);k_B=0;a_=1e-4;rho=0.5;c=1e-4;length_x=length(x0);I=eye(length_x);Hk=I;while res_B>eps&&k_B<=10000 dk=-

13、Hk*gk; m=0; while m<=5000 if lag(x0+a_*rhom*dk,lambda,cigma)-lag(x0,lambda,cigma)<=c*a_*rhom*gk'*dk mk=m; break; end m=m+1; end ak=a_*rhomk; xk=x0+ak*dk; delta=xk-x0; y=grad(xk,lambda,cigma)-gk; Hk=(I-(delta*y')/(delta'*y)*Hk*(I-(y*delta')/(delta'*y)+(delta*delta')/(del

14、ta'*y); k_B=k_B+1;x0=xk;gk=y+gk;res_B=norm(gk);end%增广拉格朗日法function val_min=ALM(x0,eps)lambda=zeros(4,1);cigma=5;alpha=10;k=1; res=abs(hf(x0),0,0,0;for i=1:3 res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma);endres=max(res);while res>eps&&k<1000xk=BFGS(x0,eps,lambda,cigma);lambda(1)=

15、lambda(1)+cigma*hf(xk);for i=1:3 lambda(i+1)=lambda(i+1)+min(0,lambda(i+1)+gi(x0,1);endk=k+1;cigma=alpha*cigma;x0=xk;res=norm(hf(x0),0,0,0;for i=1:3 res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma);endres=max(res);endval_min=fun(xk);fprintf('k=%dn',k);fprintf('fmin=%.4fn',val_min);fp

16、rintf('x=%.4f;%.4fn',xk(1),xk(2);>> x0=0;0;>> val_min=ALM(x0,1e-4)k=10fmin=-31.4003x=1.0984;4.8779val_min = -31.4003上机大作业（三）A=1 1;-1 0;0 -1;n=2;b=1;0;0;G=0.5 0;0 2;c=2 4;cvx_solver sdpt3cvx_beginvariable x(n)minimize (x'*G*x-c*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value