優(yōu)化方法作業(yè)第二版

1、優(yōu)化方法上機(jī)大作業(yè) 院 系：化工與環(huán)境生命學(xué)部姓 名：李翔宇學(xué) 號：31607007指導(dǎo)教師:肖現(xiàn)濤第一題：1. 最速下降法源程序如下：function x_star = ZSXJ(x0,eps) gk = grad(x0); res = norm(gk); k = 0; while res > eps && k<=10000 dk = -gk; ak =1; f0 = fun(x0); f1 = fun(x0+ak*dk); slope = dot(gk,dk); while f1 > f0 + 0.0001*ak*slope ak = ak/2; xk =

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

3、= 200*(x(2)-x(1)2); end 運行結(jié)果：>> x0=0,0'>> esp=1e-4;>> xk=ZSXJ(x0,eps)-The 1-th iter, the residual is 13.372079-The 2-th iter, the residual is 12.079876-The 3-th iter, the residual is 11.054105-The 9144-th iter, the residual is 0.000105-The 9145-th iter, the residual is 0.000102

4、-The 9146-th iter, the residual is 0.000100xk = 0.99990.9998MATLAB截屏：2. 牛頓法源程序如下：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

5、); 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; end function 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);

6、g(1)=2*(x(1)-1)+400*x(1)*(x(1)2-x(2); g(2) = 200*(x(2)-x(1)2); end 運行結(jié)果：>> 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.00001.0000MATALB截屏;3. BFGS方法源程序如下：function x_star = Bfgs(x0,eps) g0 = grad(x0); gk=g

7、0; 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)-f

8、a0*(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*

9、(x(1)-1)+400*x(1)*(x(1)2-x(2); g(2) = 200*(x(2)-x(1)2); end 運行結(jié)果：>> x0=0,0'>> esp=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 1516-th iter, the residual is 0.000368-The 1517-

10、th iter, the residual is 0.000099xk = 1.0001 1.0002MATLAB截屏：4. 共軛梯度法源程序如下：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 =

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

12、(1)*x(2)+2*x(1)-2;g(2)=-200*x(1)2+200*x(2);end 運行結(jié)果：>> 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 resid

13、ual 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.9999MATLAB截屏：第二題：解 ：目標(biāo)函數(shù)文件 f1.mfunction f=f1(x)f=4*x(1)-x(2)2-12;等式約束函數(shù)文件 h1.mfunction he=h1(x)he=25-x(1)2-x(2)2;不等式約束函數(shù)文件 g1.mfunction gi=g1(x)gi=1

14、0*x(1)-x(1)2+10*x(2)-x(2)2-34;目標(biāo)函數(shù)的梯度文件 df1.mfunction g=df1(x)g = 4, 2.0*x(2)'等式約束（向量）函數(shù)的Jacobi矩陣（轉(zhuǎn)置）文件dh1.mfunction dhe=dh1(x)dhe = -2*x(1), -2*x(2)'不等式約束（向量）函數(shù)的Jacobi矩陣（轉(zhuǎn)置）文件dg1.mfunction dgi=dg1(x)dgi = 10-2*x(1), 10-2*x(2)'然后在 Matlab 命令窗口輸入如下命令:x0=0,0;x,mu,lambda,output=multphr('

15、f1','h1','g1','df1','dh1','dg1',x0);得到如下輸出:x =4.89871742648821算法編程利用程序調(diào)用格式第三題：1.解：將目標(biāo)函數(shù)改寫為向量形式：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運算結(jié)果:Calling SDPT3 4.0: 7 vari

16、ables, 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 pinfe

17、as 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| c

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

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

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

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

22、alue = 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)

23、, 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. 程序

24、代碼：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運行結(jié)果：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 e

25、xpon 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-00

26、3|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| 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.