优化方法作业第二版

1、优化方法上机大作业院 系：化工与环境生命学部 姓 名：李翔宇学 号： 31607007 指导教师 :肖现涛第一题:编写程序求解下述问题min y(z) = (1 xi)2 + 100(x2 xj)2.初始点取/ = 0,精度取s = le - 4,步长由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 = ak/2;xk

2、 = 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) = 200*(x(2)-x(1

3、)A2);end 运行结果： 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 residual is 0.

4、000100 xk =0.99990.9998MATLAB 截屏: xO= Qj 01 ： esp=le-4: xk=ZSXJ txOj eps)一 TheI-thit er?theresidualis13.372075一 The2-thit er jtheresidualis12,07&876一 Iht3-thit亡theresidualis11. 054105一 TheWhit er jtheresidualis10. 421221一一Ih?5-thit ertheresidualis10. 020369Tl_ -:丄一一丄 ll_ -I; Ji1J _一 The9142-thiter,t

5、heresidualis0.000111- The9143-thiterjtheresidualis0.000108一一 The9144-thiterthe匸esidualisD. 000105一-The9145-thiterfth?residualis0. 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 ep

6、s & 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.000000xk =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 is 0.OOQOQOxk =L 0000|L 00003. BFGS 方法源程序如

7、下：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)*(y0)+(fa0*(fa0)/(fa0)*(y0)

8、;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运行结果： x0=0,0; esp=1e-4; xk= Bfgs(x0,eps)-The 1-th iter, th

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

10、erjtheresidualis2.381565一 The3thittheresidualis3.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 = gr

11、ad(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_star = xk;endfunction f=fun(x)f=(1-x(1)A2+100

12、*(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运行结果： 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, the residual is 4.527780-The 4-th iter, the re

13、sidual 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.9999MATLAB 截屏: ito= ot : esp=l*-4 : xk=Conj eps) Thel-thit efjtheresidualis3. 271712Ihe2-thit erjthere

14、sidualIS1. 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丄Mth UUUUZJ一 The70-thiter,theresidualis0. 001423一 The71-thiter,theresidualis0. 002841-The72-thitertheresidualis0. 002945一 T

15、he73-thitetjtheresidualis0. 001532 Th?74-thiter,theresidualis0. 000403 The75-thittheresidualis0. 000134一 Iht76-thiter,ther*sidualis0.0000570.99990-9999第二题:编写程序利用增广拉格朗日方法求解下述问题初始点取= 0,精度取s= lc-4.解：目标函数文件fl.mfunction f=f1(x)f=4*x(1)-x(2)A2-12;等式约束函数文件hl.mfunction he=h1(x)he=25-x(1)A2-x(2)A2;不等式约束函数文件g

16、l.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-34;目标函数的梯度文件 dfl.mfunction g=df1(x)g = 4, 2.0*x(2);等式约束(向量)函数的 Jacobi矩阵(转置)文件 dhl.mfunction dhe=dh1(x)dhe = -2*x(1), -2*x(2);不等式约束(向量)函数的 Jacobi矩阵(转置)文件 dg1.mfunction dgi=dg1(x)dgi = 10-2*x(1), 10-2*x(2);然后在Matlab命令窗口输入如下命令：x0=0,0;X,mu,lambda,outp

17、ut=multphr(f1,h1,g1,df1,dh1,dg1,x0);得到如下输出：x =4.898717426488211.00128197198571算法编程利用程序调用格式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-feval (hff x); gi=feval (gf x);n=length(K); 1=1 eng th (

18、he); nFlength(gi),rru=O l*ones(lj 1); lairbda=O l*ones(nA 1); 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);endendk二k+1;btaold=btak;xO=x;endf=f eval (fun, x);output. f

19、ral=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, xO, varargin ： );if(norm(gk)epsilonl), break; enddk=-Bkgk;n=

20、0; ink=O;while(m20)亠newf=feval (fun, xO+rhoirFdk, varacrgin ： ); oldf=feval (fun, xO, varargin ： );i f (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 ： );第三题:下载安装 CVXh http: 利用CVX编写代码求解下述问题nun1-22x1 4j*2subject to迟1 + 牝 一 1 0.工2 0利用CV

21、X编写代码求解下述问题min 3j：i X2 3帀S.t.2淤十念2 +3 W 2x + 2x2 + 3x3 0.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 =0cvx_end运算结果：Calling SDPT3 4.0: 7 variables, 3 equality constraintsFor improved efficiency, SDPT3 is solving the dual problem

22、.num. 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 cputime0|0.000|0.000|8.0e-001|6.5e+000|3.1e+002| 1.000000e+001

23、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-001| 4.166959e-0011.181563e-001| 0:0:01| chol4|0.892|0.8

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

25、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.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.2

26、e-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:01| 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

27、(relative gap, infeasibilities) 1.49e-008number of iterations = 14primal objective value = 2.50000004e-001 dual objective 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.

28、 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.07termination code= 0DIMACS: 3.3e-012 0.0e+000 1.3e-012 0.0e+000 4.

29、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 运行结果：Calling SDPT3 4.0: 6 variables, 3 equality constraintsnum. of constraints = 3dim. of linear var = 6*SDPT3: Infeasible

30、path-following algorithms*version predcorr gam expon scale_dataNT 1 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 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.622