1、優(yōu)化方法上機(jī)大作業(yè)院 系：化工與環(huán)境生命學(xué)部 姓 名：李翔宇學(xué) 號(hào)： 31607007 指導(dǎo)教師 :肖現(xiàn)濤第一題:編寫程序求解下述問(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<=10000 dk = -gk;

2、ak =1; f0 = fun(xO);f1 = fun(xO+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.0001*ak*slope ak = 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 + 10

3、0*以(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) = 200*(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 resi

4、dual 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 residual is 0.000100 xk =0.99990.9998MATLAB 截屏:>> xO= Qj 01 ：>> esp=le-4:>> xk=ZSXJ txOj eps)一 TheI-thit er?theresidualis13.372075一 The2-thit er jtheresi

5、dualis12,07&876一 Iht3-thit亡theresidualis11. 054105一 TheWhit er jtheresidualis10. 421221一一Ih?5-thit ertheresidualis10. 020369Tl_ -:丄一一丄 ll_ -I; Ji'1J _一 The9142-thiter,theresidualis0.000111- The9143-thiterjtheresidualis0.000108一一 The9144-thiter>the匸esidualisD. 000105一-The9145-thiterfth?res

6、idualis0. 000102一 Th?9143-thiter,theresidualis0. 0001000.99990.99932. 牛頓法源程序如下：function x_star = NEWTON (x0,eps) gk = grad(x0);bk = grad2(xO)A(-1);res = norm(gk);k = 0;while res > eps && k<=1000dk=-bk*gk; xk=x0+dk;k = k+1;x0 = xk;gk = grad(xk);bk = grad2(xk)A(-1);res = norm(gk);fprintf

7、('-The %d-th iter, the residual is %fn',k,res); endx_star = xk;endfunction f = fun(x)f = (1-x(1)A2 + 100*(x(2)-x(1)A2)A2;end function g = grad2(x)g = zeros(2,2); g(1,1)=2+400*(3*x(1)A2-x(2);g(1,2)=-400*x(1);g(2,1)=-400*x(1);g(2,2)=200;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400

8、*x(1)*(x(1)A2-x(2);g(2) = 200*(x(2)-x(1)A2);end運(yùn)行結(jié)果：>> 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.000000xk =1.00001.0000MATALB 截屏;» x0= 0, DY :>> esp=1;» xkNEirONtaO, eps)The 1-th iter, the residua

9、l is 447,213595 一The 2-th iter, the residual is 0.OOQOQOxk =L 0000|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<=1000dk = -Hk*gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 +

10、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)'*(y0)+(fa0*(fa0)')/(fa0)'*(y0);res = norm(gk);fprintf('-The %d-th iter, the residual is %fn',k,r

11、es); 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= Bfgs(x0,eps)-The 1-th iter, the residual is 3.271712-The 2-th iter

12、, 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 is 0.000099 xk =1.00011.0002MATLAB截屏:x0= 0, 0J :esp=le-4: xk=Bfgs (xOj eps)一 The1thit erjtheresidualis3. 271712一 The2thit erjtheresidualis2.381565一 The3thitth

13、eresidualis3.448742一The4thit er,theresidualis3, 162431 The5thit erjtheresidualis2. 989084The 1515-thiter,theresidualis0.000108The 1516-thiter,theresidualis0.000368The 1517-thiter,theresidualis0. 0000991.00011.00024. 共軛梯度法源程序如下：function x_star =Conj (x0,eps)gk = grad(x0);res = norm(gk);k = 0;dk = -gk

14、;while res > eps && k<=1000ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk); while f1 > 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

15、is %fn',k,res);endx_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

16、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.000057 xk =0.99990.9999M

17、ATLAB 截屏:» ito= ot :>> esp=l*-4 :>> xk=Conj eps) Thel-thit efjtheresidualis3. 271712Ihe2-thit erjtheresidualIS1. 3805423-thit ertheresidualIS4. 527780I he4-thit ertheresidualIS0.8505&6I he5-thit erjtheresidualis0.559005Tin 戶尺一+ Vii十戶t十Har 戶 uT Hit-alH QRR7J.il丄丄比U 3111丄 LBL,LILS

18、丄Mth UUUUZJ一 The70-thiter,theresidualis0. 001423一 The71-thiter,theresidualis0. 002841-The72-thitertheresidualis0. 002945一 The73-thitetjtheresidualis0. 001532 Th?74-thiter,theresidualis0. 000403 The75-thittheresidualis0. 000134一 Iht76-thiter,ther*sidualis0.0000570.99990-9999第二題:編寫程序利用增廣拉格朗日方法求解下述問(wèn)題初始

19、點(diǎn)取= 0,精度取s= lc-4.解：目標(biāo)函數(shù)文件fl.mfunction f=f1(x)f=4*x(1)-x(2)A2-12;等式約束函數(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*

20、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('f1','h1','g1','df1','dh1','dg1',x0);得到如下輸出：x =4.898717426488211.00128197198571算法編程利用程序調(diào)用格式function x,

21、 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-feval (hff x); gi=feval (gf> x);n=length(K); 1=1 eng th (he); nFlength(gi),rru=O> l*ones(lj 1); lairbda=O» l*ones(nA 1); btak=10; btaold=10;hile(b

22、tak>epsilon & k<maxk)爲(wèi) ival, ik=bfgsC 叩si"," dirpsi", xO* fun, hf, gf, dfun, dhf, dgf, mf 1 ambda, sigma); ink二ink+ik;he=feval (h£ x); gi=feval (gf, x);btak-0,0;丹for (1=1:1) btak=btak+he(i) "2; endfor i=l:mteitp=min(gi (i), lairbda(i)/sigma); btak-btak+temp"2;

23、end btakzsqrt(btak); if btak>epsilon 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);endendk二k+1;btaold=btak;xO=x;endf=f eval (fun, x);output. fral=f; output. iter=k;output* inner_iter-ink;output,bta-

24、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(k<maxk)gk=feval (gfun, xO, varargin ： );if(norm(gk)<epsilonl), break; enddk=-Bkgk;n»=0; ink=O;while(m<20)亠newf=feval (f

25、un, xO+rhoirFdk, varacrgin ： ); oldf=feval (fun, xO, varargin ： );i f (newf <ol df+si gir)al*rho m*gk，*dk) rnk=m; break;endm=irrH;endx= xO+rho ink*dk;sk=xxO; yk=feval (gfun x, varargin ： ) -gk;i f (yk *sk>0)Bk=Bk(Bk*sk*skJ *Bk)/(sk *Bk*sk)+(yk*ykJ )/(yk *sk); end k=k+l;xO=x;endval = f eval (fu

26、n, xO, v ar argin ： );第三題:下載安裝 CVXh http: 利用CVX編寫代碼求解下述問(wèn)題nun1-22x1 4j*2subject to遲1 + 牝 一 1 < 0. xi > 0.工2 > 0利用CVX編寫代碼求解下述問(wèn)題min 3j：i X2 3帀S.t.2淤十念2 +3 W 2x + 2x2 + 3x3 < 5 2叭 + 2x2 + 巧 W 6 37 > 0.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)mi

27、nimize( x'*a*x-b*x)subject toc * x <= 1x>=0cvx_end運(yùn)算結(jié)果：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 a

28、lgorithms*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|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.

29、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-001| 4.166959e-0011.181563e-001| 0:0:01| chol4|0.892|0.877|6.4e-008|1.6e-004|5.2e-002| 2.773022e-0012.265122e-001| 0:0:01| chol5|1.000|1.000|1.0e-008|7.6e-006|1.5e-002| 2.579468e-

30、0012.427203e-001| 0:0:01| chol6|0.905|0.904|3.1e-009|1.4e-006|2.3e-003| 2.511936e-0012.488619e-001| 0:0:01| chol7|1.000|1.000|6.1e-009|7.7e-008|6.6e-004| 2.503336e-0012.496718e-001| 0:0:01| chol8|0.903|0.903|1.8e-009|1.5e-008|1.0e-004| 2.500507e-0012.499497e-001| 0:0:01| chol9|1.000|1.000|4.9e-010|3

31、.5e-010|2.9e-005| 2.500143e-0012.499857e-001| 0:0:01| chol10|0.904|0.904|5.7e-011|1.3e-010|4.4e-006| 2.500022e-0012.499978e-001| 0:0:01| chol11|1.000|1.000|5.2e-013|1.1e-011|1.2e-006| 2.500006e-0012.499994e-001| 0:0:01| chol12|1.000|1.000|5.9e-013|1.0e-012|1.8e-007| 2.500001e-0012.499999e-001| 0:0:0

32、1| chol13|1.000|1.000|1.7e-012|1.0e-012|4.2e-008| 2.500000e-0012.500000e-001| 0:0:01| chol14|1.000|1.000|2.3e-012|1.0e-012|7.3e-009| 2.500000e-0012.500000e-001| 0:0:01|stop: max(relative gap, infeasibilities) < 1.49e-008number of iterations = 14primal objective value = 2.50000004e-001 dual object

33、ive value = 2.49999996e-001 gap := trace(XZ) = 7.29e-009relative gap= 4.86e-009actual relative gap = 4.86e-009rel. primal infeas (scaled problem) = 2.33e-012rel. dual """ = 1.00e-012rel. primal infeas (unscaled problem) = 0.00e+000rel. dual """ = 0.00e+000norm(X), norm(

34、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.07termination 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

35、 -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 <= bx>=0 cvx_end 運(yùn)行結(jié)果：Calling SDPT3 4.0: 6 variables, 3 equality constraintsnum. of constraints = 3dim. of linear var = 6*SDPT3: Infeasible path-following algorithms*version predcorr gam expon scale_dataNT 1

36、0.000 1 0it 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 11|0.912|1.000|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 1