优化方法上机大作业.docx

1、优化方法上机大作业上机大作业：编写程序求解下述问题 min x f(x) = (1x1)2 + 100(x2 x12)2. 初始点取 x 0 = 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 xk = x0 + ak*dk; f1 = fun(xk); end k = k+1;x0 = xk; gk =

2、 grad(xk);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-4eps =1.0000e-004 xk

3、=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

4、.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 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.00001.00003.BFGS方法源程序如下：function x_star = bfgs(x

5、0,eps) g0 = grad(x0);gk=g0;res = norm(gk);Hk=eye(2);k = 0; while res eps 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

6、 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

7、 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

8、= norm(gk);k = 0; dk = -gk;while res eps 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)

9、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 r

10、esidual 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.99990.9999上机大作业：编写程序利用增广拉格朗日方法求解下述问题 min 4x1 x2 2 12s.t. 25x12 x22 = 0 10x1 x12 + 10x2 x22 34 0x 1,x2 0初始点取 x 0 = 0, 精度取 = 1e4. 主程序：function x,mu,lamda,

11、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 lamda=max(0.0,lamda-sigma*gi);if(k=2 endendk=k+1;bet

12、aold=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)

13、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=lengt

14、h(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(k0)Bk=Bk-(Bk*sk*sk*Bk)/(sk*Bk*sk)+(yk*yk)/(yk*sk);endk=k+

15、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);不等

16、式函数梯度文件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.00134.8987mu =0.0158lamda =0.554200output = fval: -31.9924iter: 5inner_iter: 33beta: 8.4937e-005上机大作业：1.解：将目标函数改写为向量形式：x*a*x-b*x程序代码：n=2;a=0.5,0;0,1;b=2 4;c=1 1;cvx_beginvariab

17、le x(n)minimize( x*a*x-b*x)subject toc * x =0cvx_end运算结果:Calling SDPT3 4.0: 7 variables, 3 equality constraintsFor improved efficiency, SDPT3 is solving the dual problem.-num. of constraints = 3dim. of socp var = 4, num. of socp blk = 1dim. of linear var = 3*SDPT3: Infeasible path-following algorith

18、ms*version predcorr gam expon scale_dataNT 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

19、|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|

20、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|

21、 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|

22、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) =0cvx_end运行结果：Calling SDPT

23、3 4.0: 6 variables, 3 equality constraints-num. of constraints = 3dim. of linear var = 6*SDPT3: Infeasible path-following algorithms*version predcorr gam expon scale_dataNT 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+00

24、1 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

25、 4|0.881|1.000|2.2e-008|4.6e-005|3.7e-001|-5.255152e+000 -5.622375e+000| 0:0:00| chol 1 1 5|0.995|0.962|1.6e-009|6.2e-006|1.5e-002|-5.394782e+000 -5.409213e+000| 0:0:00| chol 1 1 6|0.989|0.989|2.7e-010|5.2e-007|1.7e-004|-5.399940e+000 -5.400100e+000| 0:0:00| chol 1 1 7|0.989|0.989|5.3e-011|5.8e-009|

26、1.8e-006|-5.399999e+000 -5.400001e+000| 0:0:00| chol 1 1 8|1.000|0.994|2.8e-013|4.3e-011|2.7e-008|-5.400000e+000 -5.400000e+000| 0:0:00|stop: max(relative gap, infeasibilities) 1.49e-008-number of iterations = 8primal objective value = -5.39999999e+000dual objective value = -5.40000002e+000gap := tr

27、ace(XZ) = 2.66e-008relative gap = 2.26e-009actual relative gap = 2.21e-009rel. primal infeas (scaled problem) = 2.77e-013rel. dual “ “ “ = 4.31e-011rel. primal infeas (unscaled problem) = 0.00e+000rel. dual “ “ “ = 0.00e+000norm(X), norm(y), norm(Z) = 4.3e+000, 1.3e+000, 1.9e+000norm(A), norm(b), norm(C) = 6.7e+000, 9.1e+000, 5.4e+000Total CPU time (secs) = 0.11 CPU time per iteration = 0.01 termination code = 0DIMACS: 3.6e-013 0.0e+000 5.8e-011 0.0e+000 2.2e-009 2.3e-009-Status: SolvedOptimal value (cvx_optval): -5.4