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

1、優(yōu)化方法上機(jī)大作業(yè)院 系：化工與環(huán)境生命學(xué)部 姓 名：李翔宇學(xué) 號(hào)： 31607007 指導(dǎo)教師 :肖現(xiàn)濤第一題:編寫(xiě)程序求解下述問(wèn)題min y(z) = (1 xi)2+ 100(x2 xj)2.初始點(diǎn)取/ = 0,精度取 s = le - 4,步長(zhǎng)由 Armijo 線捜索生成，方向 分別由下列方法生成：最速下降法BFGS 方法共輒梯度法1.最速下降法源程序如下：function x_star = ZSXJ (x0,eps) gk = grad(xO);res = norm(g k);k = 0;while res eps & k f0 + 0.0001*ak*slope ak =

2、ak/2;xk = x0 + ak*dk;f1 = fun(xk);endk = 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)A2 + 100*以以(2)咲咲(1)八八2)八八2;end function g = grad(x) g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)A2-x(2); g(2) = 20

3、0*(x(2)-x(1)A2);end運(yùn)行結(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 -The 9146-th iter, the residu

4、al is0.000100 xk =0.99990.9998MATLAB截屏: xO= Qj 01： esp=le-4: xk=ZSXJ txOj eps)一TheI-th it er?theresidualis13.372075一The 2-thit er jthe residualis12,07&876一Iht3-thit亡theresidualis 11. 054105一TheWh it er jtheresidualis 10. 421221一一Ih?5-th it ertheresidualis10. 020369Tl_ -:丄一一丄ll_ -_I; Ji_ 1J _一The

5、9142-th iter,theresidualis 0.000111- The9143-th iterjtheresidualis 0.000108一一The9144-th iterthe匸esidualis D. 000105一-The9145-th iterfth?residualis 0. 000102一Th?9143-th iter,theresidualis 0. 0001000.99990.99932. 牛頓法源程序如下：function x_star = NEWTON (x0,eps) gk = grad(x0);bk = grad2(xO)A(-1);res = norm(g

6、k);k = 0;while res eps & k xO=O,O;eps=1e-4; xk= NEWTON (x0,eps)-The 1-th iter, the residual is 447.213595-The 2-th iter, the residual is 0.000000 xk =1.00001.0000MATALB截屏; x0= 0,DY : esp=1; xkNEirONtaO, eps)The 1-th iter, the residual is 447,213595一The 2-th iter, the residual is0.OOQOQOxk =L 000

7、0|L 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 f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endk = 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)*(

8、y0)+(fa0*(fa0)/(fa0)*(y0);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)A2 + 100*債債(2)咲咲(1)八八2)八八2;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)A2-x(2);g(2) = 200*(x(2)-x(1)A2);end運(yùn)行結(jié)果： x0=0,0; esp=1e-4; xk= B

9、fgs(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-th iter, the residual is0.000099 xk =1.00011.0002x0= 0, 0J:esp=le-4: xk=Bfgs (xOj eps)一The1 thit erjtheresidualis 3

10、. 271712一The2 thit erjtheresidualis 2.381565一The3 thittheresidualis 3.448742一The4 thit er,the residualis 3, 162431 The 5 thit erjtheresidualis 2. 989084MATLAB截屏:The 1515-thiter, theresidualis 0.000108The 1516-thiter, theresidualis 0.000368The 1517-thiter, theresidualis 0. 0000991.00011.00024. 共軛梯度法源

11、程序如下：function x_star =Conj (x0,eps)gk = grad(x0);res = norm(gk);k = 0;dk = -gk;while res eps & k f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(xk);endd0=dk;g0=gk;k=k+1;xO=xk;gk=grad(xk);f=(norm(gk)/norm(gO)A2;res=norm(gk);dk=-gk+f*d0;fprintf(-The %d-th iter, the residual is %fn,k,res);endx

12、_star = xk;endfunction f=fun(x)f=(1-x(1)A2+100*(x(2)-x(1)A2)A2;endfunction g=grad(x)g=zeros(2,1);g(1)=400*x(1)A3-400*x(1)*x(2)+2*x(1)-2;g(2)=-200*x(1)A2+200*x(2);end運(yùn)行結(jié)果： xO=O,O; eps=1e-4; xk=Conj(xO,eps)-The 1-th iter, the residual is 3.271712-The 2-th iter, the residual is 1.380542-The 3-th iter,

13、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 is0.000134 -The 76-th iter, the residual is 0.000057 xk =0.99990.9999MATLAB截屏: ito=OT: esp=l*-4 : xk=Conj eps) Thel-thit ef

14、jtheresidualis 3. 271712Ihe2-thit erjtheresidualIS 1. 3805423-thit ertheresidualIS 4. 527780I he4-thit ertheresidualIS 0.8505&6I he5-thit erjtheresidualis 0.559005Tin尺一+Vi i十戶(hù)T十Har戶(hù) uT HIT-alH QRR7J.il丄丄比U 3111丄LBL,LILS丄 Mth UUUUZJ一The70-th iter,theresidualis 0. 001423一The71-th iter,theresiduali

15、s 0. 002841-The72-th itertheresidualis 0. 002945一The73-th itetjtheresidualis 0. 001532 Th? 74-th iter,theresidualis 0. 000403 The 75-thittheresidualis 0. 000134一Iht76-th iter,ther*sidualis 0.0000570.99990-9999第二題:編寫(xiě)程序利用增廣拉格朗日方法求解下述問(wèn)題初始點(diǎn)取= 0,精度取 s= lc-4.解：目標(biāo)函數(shù)文件fl.mfunction f=f1(x)f=4*x(1)-x(2)A2-12;

16、等式約束函數(shù)文件hl.mfunction he=h1(x)he=25-x(1)A2-x(2)A2;不等式約束函數(shù)文件gl.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-34;目標(biāo)函數(shù)的梯度文件dfl.mfunction g=df1(x)g = 4, 2.0*x(2);等式約束(向量)函數(shù)的Jacobi矩陣(轉(zhuǎn)置)文件dhl.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),

17、 10-2*x(2);然后在Matlab命令窗口輸入如下命令：x0=0,0;X,mu,lambda,output=multphr(f1,h1,g1,df1,dh1,dg1,x0);得到如下輸出：x =4.898717426488211.00128197198571算法編程利用程序調(diào)用格式function x, mu, lambda, out put =mul t phr (furij hf, gf, dfun, dhfj dgf, xO)paxk=500;si gma=2* 0; etaz2. 0; thetazO* 8;k=0; ink=0;epsilon=le-4;x=xO; he-fev

18、al (hffx); gi=feval (gf x);n=length(K); 1=1 eng th (he); nFlength(gi),rru=O l*ones(lj 1); lairbda=O l*ones(nA1); btak=10; btaold=10;hile(btakepsilon & kepsilon if(k=2 & btak theta*btaold) sigeta+sigma,end foriDui)zmu(i)-sigina*he(i); endfor (i=l：nOlambda(i)zmax(0.0, lairbda(i)-signia*gi(i);e

19、ndendk 二 k+1;btaold=btak;xO=x;endf=f eval (fun, x);output. fral=f; output. iter=k;output* inner_iter-ink;output,bta-btak;function x, val, k=bfgs (fun gfun, xO, varargin) maxk=500;rho=0. 55; sigmal=O. 4; epsi1onl = le5;k=0; n=length(xO);Bk=eye (n) ; %Bk=feval ( Hess , xO);while(kmaxk)gk=feval (gfun,

20、xO, varargin ： );if(norm(gk)epsilonl), break; enddk=-Bkgk;n=0; ink=O;while(m20)亠newf=feval (fun, xO+rhoirFdk, varacrgin ： ); oldf=feval (fun, xO, varargin ： );if(newf 0)Bk=Bk(Bk*sk*skJ*Bk)/(sk *Bk*sk)+(yk*ykJ)/(yk *sk); end k=k+l;xO=x;endval = f eval (fun, xO, v ar argin ： );第三題:下載安裝 CVXhhttp: CVX 編

21、寫(xiě)代碼求解下述問(wèn)題遲 1 +牝 一 1 0.工 2 0利用 CVX編寫(xiě)代碼求解下述問(wèn)題min3J：I X2 3 帀S.t.2淤十念 2 +3 W 2x + 2x2+ 3x3 0.1.解：將目標(biāo)函數(shù)改寫(xiě)為向量形式: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 =0cvx_end運(yùn)算結(jié)果：Calling SDPT3 4.0: 7 variables, 3 equality constraintsFor improved efficiency, S

22、DPT3 is solving the dual problem.nun1-22x1 4J*2subject tonum. of constraints = 3dim. of socp var = 4, num. of socp blk = 1dim. of linear var = 3*SDPT3: Infeasible path-following algorithms*version predcorr gam expon scale_dataNT 1 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime

23、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 11|1.000|0.987|4.3e-007|1.5e-001|1.6e+001| 9.043148e+000 -2.714056e-001| 0:0:01| chol 1 12|1.000|1.000|2.6e-007|7.6e-003|1.4e+000| 1.234938e+000 -5.011630e-002| 0:0:01| chol 1 13|1.000|1.000|2.4e-007|7.6e-004|3.0e-

24、001| 4.166959e-001 1.181563e-001| 0:0:01| chol4|0.892|0.877|6.4e-008|1.6e-004|5.2e-002| 2.773022e-001 2.265122e-001| 0:0:01| chol5|1.000|1.000|1.0e-008|7.6e-006|1.5e-002| 2.579468e-001 2.427203e-001| 0:0:01| chol6|0.905|0.904|3.1e-009|1.4e-006|2.3e-003| 2.511936e-001 2.488619e-001| 0:0:01| chol7|1.0

25、00|1.000|6.1e-009|7.7e-008|6.6e-004| 2.503336e-001 2.496718e-001| 0:0:01| chol8|0.903|0.903|1.8e-009|1.5e-008|1.0e-004| 2.500507e-001 2.499497e-001| 0:0:01| chol9|1.000|1.000|4.9e-010|3.5e-010|2.9e-005| 2.500143e-001 2.499857e-001| 0:0:01| chol10|0.904|0.904|5.7e-011|1.3e-010|4.4e-006| 2.500022e-001

26、 2.499978e-001| 0:0:01| chol11|1.000|1.000|5.2e-013|1.1e-011|1.2e-006| 2.500006e-001 2.499994e-001| 0:0:01| chol12|1.000|1.000|5.9e-013|1.0e-012|1.8e-007| 2.500001e-001 2.499999e-001| 0:0:01| chol13|1.000|1.000|1.7e-012|1.0e-012|4.2e-008| 2.500000e-001 2.500000e-001| 0:0:01| chol14|1.000|1.000|2.3e-

27、012|1.0e-012|7.3e-009| 2.500000e-001 2.500000e-001| 0:0:01|stop: max(relative gap, infeasibilities) 1.49e-008number of iterations = 14primal objective value = 2.50000004e-001dual objective value = 2.49999996e-001gap := trace(XZ) = 7.29e-009relative gap= 4.86e-009actual relative gap = 4.86e-009rel. p

28、rimal infeas (scaled problem) = 2.33e-012rel. dual = 1.00e-012rel. primal infeas (unscaled problem) = 0.00e+000rel. dual = 0.00e+000norm(X), norm(y), norm(Z) = 3.2e+000, 1.5e+000, 1.9e+000norm(A), norm(b), norm(C) = 3.9e+000, 4.2e+000, 2.6e+000Total CPU time (secs) = 0.99CPU time per iteration = 0.0

29、7termination code= 0DIMACS: 3.3e-012 0.0e+000 1.3e-012 0.0e+000 4.9e-009 4.9e-009Status: 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_beginvariable x(n) minimize( a*x)subject toC * x =0 cvx_end運(yùn)行結(jié)果：Calling SDPT3 4.0: 6 variables, 3 equality constraintsn

30、um. of constraints = 3dim. of linear var = 6*SDPT3: Infeasible path-following algorithms*version predcorr gam expon scale_dataNT 1 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime0|0.000|0.000|1.1e+001|5.1e+000|6.0e+002|-7.000000e+001 0.000000e+000| 0:0:00| chol 1 11|0.912|1.000

31、|9.4e-001|4.6e-002|6.5e+001|-5.606627e+000 -2.967567e+001| 0:0:00| chol 1 12|1.000|1.000|1.3e-007|4.6e-003|8.5e+000|-2.723981e+000 -1.113509e+001| 0:0:00| chol 1 13|1.000|0.961|2.3e-008|6.2e-004|1.8e+000|-4.348354e+000 -6.122853e+000| 0:0:00| chol 1 14|0.881|1.000|2.2e-008|4.6e-005|3.7e-001|-5.255152e+000 -5.622375e+000|