优化方法-2016.docx

优化方法上机大作业上机大作业：编写程序求解下述问题 min x f(x) = (1x1)2 + 100(x2 x12)2. 初始点取 x0 = 0, 精度取=1e4,步长由 Armijo 线搜索生成, 方向分别由下列方法生成: 1 最速下降法2 牛顿法3 BFGS 方法 4 共轭梯度法1. 最速下降法源程序如下：function x_star = steepest(x0,eps) gk = grad(x0); res = norm(gk); k = 0; while res eps k=1000 dk = -gk; ak =1; f0 = fun(x0); f1 = fun(x0+ak*dk); slope = dot(gk,dk); while f1 f0 + 0.1*ak*slope ak = ak/2; xk = x0 + ak*dk; f1 = fun(xk); end k = k+1; x0 = xk; gk = grad(xk);res = norm(gk); fprintf(-The %d-th iter, the residual is %fn,k,res); end x_star = xk; end function f = fun(x) f = (1-x(1)2 + 100*(x(2)-x(1)2)2; endfunction g = grad(x) g = zeros(2,1); g(1)=2*(x(1)-1)+400*x(1)*(x(1)2-x(2); g(2) = 200*(x(2)-x(1)2); end 运行结果： x0=0,0;eps=1e-4eps = 1.0000e-004 xk=steepest(x0,eps)-The 1-th iter, the residual is 3.271712-The 2-th iter, the residual is 2.504194-The 3-th iter, the residual is 2.073282-The 998-th iter, the residual is 0.449447-The 999-th iter, the residual is 0.449447-The 1000-th iter, the residual is 0.449447-The 1001-th iter, the residual is 0.449447xk = 0.63690.40382. 牛顿法源程序如下：function x_star = newton(x0,eps) gk = grad(x0); bk = grad2(x0)(-1); res = norm(gk); k = 0; while res eps k=1000 dk=-bk*gk; xk=x0+dk; k = k+1; x0 = xk; gk = grad(xk); bk = grad2(xk)(-1); res = norm(gk); fprintf(-The %d-th iter, the residual is %fn,k,res); end x_star = xk; end function f = fun(x) f = (1-x(1)2 + 100*(x(2)-x(1)2)2; endfunction g = grad2(x) g = zeros(2,2); g(1,1)=2+400*(3*x(1)2-x(2); g(1,2)=-400*x(1); g(2,1)=-400*x(1); g(2,2)=200; end function g = grad(x) g = zeros(2,1); g(1)=2*(x(1)-1)+400*x(1)*(x(1)2-x(2); g(2) = 200*(x(2)-x(1)2); end 运行结果： x0=0,0;eps=1e-4; xk=newton(x0,eps)-The 1-th iter, the residual is 447.213595-The 2-th iter, the residual is 0.000000xk = 1.0000 1.00003. BFGS方法源程序如下：function x_star = bfgs(x0,eps) g0 = grad(x0); gk=g0; res = norm(gk); Hk=eye(2); k = 0; while res eps k=1000 dk = -Hk*gk; ak =1; f0 = fun(x0); f1 = fun(x0+ak*dk); slope = dot(gk,dk); while f1 f0 + 0.1*ak*slope ak = ak/2; xk = x0 + ak*dk; f1 = fun(xk); end k = k+1; fa0=xk-x0; x0 = xk; g0=gk;gk = grad(xk);y0=gk-g0;Hk=(eye(2)-fa0*(y0)/(fa0)*(y0)*(eye(2)-(y0)*(fa0)/(fa0)*(y0)+(fa0*(fa0)/(fa0)*(y0);res = norm(gk); fprintf(-The %d-th iter, the residual is %fn,k,res); end x_star = xk; endfunction f=fun(x) f=(1-x(1)2 + 100*(x(2)-x(1)2)2; endfunction g = grad(x) g = zeros(2,1); g(1)=2*(x(1)-1)+400*x(1)*(x(1)2-x(2); g(2) = 200*(x(2)-x(1)2); end 运行结果： x0=0,0; eps=1e-4; xk=bfgs(x0,eps)-The 1-th iter, the residual is 3.271712-The 2-th iter, the residual is 2.381565-The 3-th iter, the residual is 3.448742-The 998-th iter, the residual is 0.004690-The 999-th iter, the residual is 0.008407-The 1000-th iter, the residual is 0.005314-The 1001-th iter, the residual is 0.010880xk = 0.99550.99114. 共轭梯度法源程序如下：function x_star =conj(x0,eps) gk = grad(x0);res = norm(gk); k = 0; dk = -gk; while res eps k=1000 ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk); while f1 f0 + 0.1*ak*slope ak = ak/2; xk = x0 + ak*dk; f1 = fun(xk); end d0=dk;g0=gk; k=k+1; x0=xk;gk=grad(xk);f=(norm(gk)/norm(g0)2; res=norm(gk);dk=-gk+f*d0;fprintf(-The %d-th iter, the residual is %fn,k,res); end x_star = xk; end function f=fun(x)f=(1-x(1)2+100*(x(2)-x(1)2)2;endfunction g=grad(x)g=zeros(2,1);g(1)=400*x(1)3-400*x(1)*x(2)+2*x(1)-2;g(2)=-200*x(1)2+200*x(2);end 运行结果： x0=0,0; eps=1e-4; xk=Conj(x0,eps)-The 1-th iter, the residual is 3.271712-The 2-th iter, the residual is 1.380542-The 3-th iter, the residual is 4.527780-The 4-th iter, the residual is 0.850596-The 73-th iter, the residual is 0.001532-The 74-th iter, the residual is 0.000402-The 75-th iter, the residual is 0.000134-The 76-th iter, the residual is 0.000057xk = 0.9999 0.9999上机大作业：编写程序利用增广拉格朗日方法求解下述问题 min 4x1 x2 2 12s.t. 25x12 x22 = 0 10x1 x12 + 10x2 x22 34 0x 1,x2 0初始点取 x0 = 0, 精度取= 1e4. 主程序：function x,mu,lamda,output=main(fun,hf,gf,dfun,dhf,dgf,x0)maxk=2000;theta=0.8; eta=2.0;k=0; ink=0;ep=1e-4; sigma=0.4;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); lamda=0.1*ones(m,1);betak=10; betaold=10; while(betakep kmaxk)ik,x,val=bfgs(lagrang,dlagrang,x0,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma);ink=ink+ik;he=feval(hf,x); gi=feval(gf,x);betak=sqrt(norm(he,2)2+norm(min(gi,lamda/sigma),2)2);if betakepmu=mu-sigma*he;lamda=max(0.0,lamda-sigma*gi);if(k=2 betaktheta*betaold)sigma=eta*sigma;endendk=k+1;betaold=betak;x0=x;endf=feval(fun,x);output.fval=f;output.iter=k;output.inner_iter=ink;output.beta=betak;增广拉格朗日函数function lag=lagrang(x,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma)f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x);l=length(he); m=length(gi);lag=f; s1=0.0;for(i=1:l)lag=lag-he(i)*mu(i);s1=s1+he(i)2;endlag=lag+0.5*sigma*s1;s2=0.0;for(i=1:m)s3=max(0.0,lamda(i)-sigma*gi(i);s2=s2+s32-lamda(i)2;endlag=lag+s2/(2.0*sigma);增广拉格朗日梯度函数function dlag=dlagrang(x,fun,hf,gf,dfun,dhf,dgf,mu,lamda,sigma)dlag=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)dlag=dlag+(sigma*he(i)-mu(i)*dhe(:,i);endfor(i=1:m)dlag=dlag+(sigma*gi(i)-lamda(i)*dgi(:,i);end线搜索程序 基于BFGS方法function k,x,val=bfgs(fun,gfun,x0,varargin)Max=1000;ep=1.e-4;beta=0.55; sigma1=0.4;n=length(x0); Bk=eye(n);k=0;while(kMax)gk=feval(gfun,x0,varargin:);if(norm(gk)ep), break; enddk=-Bkgk; m=0; mk=0;while(m20) newf=feval(fun,x0+betam*dk,varargin:);oldf=feval(fun,x0,varargin:);if(newf=oldf+sigma1*betam*gk*dk)mk=m; break;endm=m+1;endx=x0+betamk*dk;sk=x-x0;yk=feval(gfun,x,varargin:)-gk;if(yk*sk0)Bk=Bk-(Bk*sk*sk*Bk)/(sk*Bk*sk)+(yk*yk)/(yk*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin:);目标函数文件function f=f1(x)f=4*x(1)-x(2)2-12;等式约束文件function he=h1(x)he=25-x(1)2-x(2)2;不等式约束function gi=g1(x)gi=zeros(3,1);gi(1)=10*x(1)-x(1)2+10*x(2)-x(2)2-34;gi(2)=x(1);gi(3)=x(2);目标函数梯度文件function g=df1(x)g=4;-2*x(1);等式函数梯度文件function dhe=dh1(x)dhe=-2*x(1);-2*x(2);不等式函数梯度文件function dgi=dg1(x)dgi=10-2*x(1),1,0；10-2*x(2)，0,1;输入数据X0=0;0x,mu,lamda,output=main(f1,h1,g1,df1,dh1,dg1,x0)输出数据x = 1.0013 4.8987mu = 0.0158lamda = 0.5542 0 0output = fval: -31.9924 iter: 5 inner_iter: 33 beta: 8.4937e-005上机大作业：1.解：将目标函数改写为向量形式：x*a*x-b*x程序代码：n=2;a=0.5,0;0,1;b=2 4;c=1 1;cvx_beginvariable x(n)minimize( x*a*x-b*x)subject toc * x = 1x=0cvx_end运算结果:Calling SDPT3 4.0: 7 variables, 3 equality constraints For improved efficiency, SDPT3 is solving the dual problem.- num. of constraints = 3 dim. of socp var = 4, num. of socp blk = 1 dim. of linear var = 3* SDPT3: Infeasible path-following algorithms* version predcorr gam expon scale_data NT 1 0.000 1 0 it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime- 0|0.000|0.000|8.0e-001|6.5e+000|3.1e+002| 1.000000e+001 0.000000e+000| 0:0:00| chol 1 1 1|1.000|0.987|4.3e-007|1.5e-001|1.6e+001| 9.043148e+000 -2.714056e-001| 0:0:01| chol 1 1 2|1.000|1.000|2.6e-007|7.6e-003|1.4e+000| 1.234938e+000 -5.011630e-002| 0:0:01| chol 1 1 3|1.000|1.000|2.4e-007|7.6e-004|3.0e-001| 4.166959e-001 1.181563e-001| 0:0:01| chol 1 1 4|0.892|0.877|6.4e-008|1.6e-004|5.2e-002| 2.773022e-001 2.265122e-001| 0:0:01| chol 1 1 5|1.000|1.000|1.0e-008|7.6e-006|1.5e-002| 2.579468e-001 2.427203e-001| 0:0:01| chol 1 1 6|0.905|0.904|3.1e-009|1.4e-006|2.3e-003| 2.511936e-001 2.488619e-001| 0:0:01| chol 1 1 7|1.000|1.000|6.1e-009|7.7e-008|6.6e-004| 2.503336e-001 2.496718e-001| 0:0:01| chol 1 1 8|0.903|0.903|1.8e-009|1.5e-008|1.0e-004| 2.500507e-001 2.499497e-001| 0:0:01| chol 1 1 9|1.000|1.000|4.9e-010|3.5e-010|2.9e-005| 2.500143e-001 2.499857e-001| 0:0:01| chol 1 1 10|0.904|0.904|5.7e-011|1.3e-010|4.4e-006| 2.500022e-001 2.499978e-001| 0:0:01| chol 2 2 11|1.000|1.000|5.2e-013|1.1e-011|1.2e-006| 2.500006e-001 2.499994e-001| 0:0:01| chol 2 2 12|1.000|1.000|5.9e-013|1.0e-012|1.8e-007| 2.500001e-001 2.499999e-001| 0:0:01| chol 2 2 13|1.000|1.000|1.7e-012|1.0e-012|4.2e-008| 2.500000e-001 2.500000e-001| 0:0:01| chol 2 2 14|1.000|1.000|2.3e-012|1.0e-012|7.3e-009| 2.500000e-001 2.500000e-001| 0:0:01| stop: max(relative gap, infeasibilities) 1.49e-008- number of iterations = 14 primal objective value = 2.50000004e-001 dual objective value = 2.49999996e-001 gap := trace(XZ) = 7.29e-009 relative gap = 4.86e-009 actual relative gap = 4.86e-009 rel. primal infeas (scaled problem) = 2.33e-012 rel. dual = 1.00e-012 rel. primal infeas (unscaled problem) = 0.00e+000 rel. dual = 0.00e+000 norm(X), norm(y), norm(Z) = 3.2e+000, 1.5e+000, 1.9e+000 norm(A), norm(b), norm(C) = 3.9e+000, 4.2e+000, 2.6e+000 Total CPU time (secs) = 0.99 CPU time per iteration = 0.07 termination code = 0 DIMACS: 3.3e-012 0.0e+000 1.3e-012 0.0e+000 4.9e-009 4.9e-009- -Status: SolvedOptimal value (cvx_optval): -32. 解：将目标函数改写为向量形式：程序代码：n=3;a=-3 -1 -3;b=2;5;6;C=2 1 1;1 2 3;2 2 1;cvx_begin variable x(n) minimize( a*x) subject to C * x = b x=0cvx_end运行结果：Calling SDPT3 4.0: 6 variables, 3 equality constraints- num. of constraints = 3 dim. of linear var = 6* SDPT3: Infeasible path-following algorithms* version predcorr gam expon scale_data NT 1 0.000 1 0 it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime- 0|0.000|0.000|1.1e+001|5.1e+000|6.0e+002|-7.000000e+001 0.000000e+000| 0:0:00| chol 1 1 1|0.912|1.000|9.4e-001|4.6e-002|6.5e+001|-5.606627e+000 -2.967567e+001| 0:0:00| chol 1 1 2|1.000|1.000|1.3e-007|4.6e-003|8.5e+000|-2.723981e+000 -1.113509e+001| 0:0:00| chol 1 1 3|1.000|0.961|2.3e-008|6.2e-004|1.8e+000|-4.348354e+000 -6.122853e+000| 0:0:00| chol 1 1 4|0.881|1.000|2.2e-008|4.6e-005|3.7e-001|-5.255152e+000 -5.622375e+000| 0:0:00| c