优化方法-2016

1、优化方法上机大作业上机大作业I :编写程序求解下述问题 min x f(x) = (1 -xi)2 + 100(x 2 - xi2) 2.初 始点取x0 = 0,精度取e=1e-4,步长由Armijo 线搜索生成，方向 分别由下列方法生成:1 最速下降法2 牛顿法3 BFGS方法4 共轭梯度法1. 最速下降法源程序如下：function x_star = steepest(x0,eps)gk = grad(x0);res = norm(gk);k = 0;while res eps & k f0 + 0.1*ak*slopeak = ak/2;xk = x0 + ak*dk;f1 = fun(

2、xk);endk = k+1;x0 = xk;gk = grad(xk);res = norm(gk);fprintf(-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;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;eps=1e

3、-4eps =1.0000e-004 xk=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- -

4、The 1001-th iter, the residual is 0.449447xk =0.63690.40382. 牛顿法 源程序如下：function x_star = newton(x0,eps)gk = grad(x0);bk = grad2(x0)F(-1);res = norm(gk);k = 0;while res eps & k 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.00

5、001.00003. BFGST法源程序如下：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)+(f

6、a0*(fa0)/(fa0)*(y0);res = norm(gk);fprintf(-The %d-th iter, the residual is %fn,k,res);endx_star = xk;endfunction f=fun(x)f=(1-x(1)42 + 100*(x(2)-x(1)A2)A2;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运行结果： x0=0,0;- eps=1e-4;- xk=bfgs(x0,ep

7、s)- -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 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.

8、010880xk =0.99550.99114. 共轭梯度法 源程序如下：f unction 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;x0=xk;gk=grad(xk);f=(norm(gk)/norm(g0)F2;res=norm(gk);dk=-gk+f*d0;fprintf(-The %d-th iter, the

9、 residual 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运行结果： x0=0,0;- eps=1e-4;- xk=Conj(x0,eps)- -The 1-th iter, the residual is 3.271712- -The 2-th iter, the re

10、sidual 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.000057xk =0.99990.9999上机大作业n ：编写

11、程序利用增广拉格朗日方法求解下述问题min 4x 1 -x2 2 - 12s.t. 25 -x12 -x22 = 010X1 太12 + 10x 2 -X22 - 34 0x 1,x 2 0初始点取X0 = 0,精度取 = 1e-4.主程序：function x,mu,lamda,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);

12、 m=length(gi);mu=0.1*ones(l,1); lamda=0.1*ones(m,1);betak=10; betaold=10;while (betakep & kepmu=mu-sigma*he;lamda=max(0.0,lamda-sigma*gi);if(k=2 & betaktheta*betaold)sigma=eta*sigma;endend k=k+1;betaold=betak;x0=x;endf=feval(fun,x);output.fval=f;output.iter=k;output.inner_iter=ink;output.beta=betak;

13、增广拉格朗日函数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(i42;endlag=lag+0.5*sigma*s1;s2=0.0;for (i=1:m)s3=max(0.0,lamda(i)-sigma*gi(i);s2=s2+s3A2-lamda(iF2;end

14、lag=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=length(gi);for (i=1:l)dlag=dlag+(sigma*he(i)-mu(i)*dhe(:,i);endfor (i=1:m)dlag=dlag+(sigma*gi(i)-lamda(

15、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 (kMax)gk=feval(gfun,x0,varargin:);if (norm(gk)ep), break ; enddk=-Bkgk;m=0; mk=0;while (m20)newf=feval(fun,x0+betaAm*dk,varargin:);oldf=feval(fun,x0,vararg

16、in:);if (newf0)Bk=Bk-(Bk*sk*sk*Bk)/(sk*Bk*sk)+(yk*yk)/(yk*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin:);目标函数文件function f=f1(x)f=4*x(1)-x(2)A2-12;等式约束文件function he=h1(x)he=25-x(142-x(2)A2;不等式约束function gi=g1(x)gi=zeros(3,1);gi(1)=10*x(1)-x(1)A2+10*x(2)-x(2)A2-34;gi(2)=x(1);gi(3)=x(2);目标函数梯度文件funct

17、ion g=df1(x)g=4;-2*x(1);等式函数梯度文件function dhe=dh1(x)dhe=-2*x(1);-2*x(2);不等式函数梯度文件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上机大作业

18、m: 下载安装 CVX, http: /cvxrxom cvx/ 利用CVX编写代码求解下述问题mirir事 + 2 211 4 耳2subject to X + x2 i 0. g 0 利用CVX编写代码求解下述问题min -3工JC2 Bwf3S.t. 2应 + X2 + 2l L + 212 + 33 52叼 + 2口 61. 解：将目标函数改写为向量形式：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运算结果:Cal

19、ling 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 algorithms*version predcorr gam expon scale_dataNT 10.000 1 it pstep

20、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.000|1.000|2.6e-007|7.6e-003|1.4e+000| 1.234938e+000 -5.011630e-002|0:0:0

21、1| chol 1 13|1.000|1.000|2.4e-007|7.6e-004|3.0e-001| 4.166959e-001 1.181563e-001| 0:0:01| chol 1 14|0.892|0.877|6.4e-008|1.6e-004|5.2e-002| 2.773022e-001 2.265122e-001| 0:0:01| chol 1 15|1.000|1.000|1.0e-008|7.6e-006|1.5e-002| 2.579468e-001 2.427203e-001| 0:0:01|chol 1 16|0.905|0.904|3.1e-009|1.4e-0

22、06|2.3e-003| 2.511936e-001 2.488619e-001| 0:0:01|chol 1 17|1.000|1.000|6.1e-009|7.7e-008|6.6e-004| 2.503336e-001 2.496718e-001| 0:0:01| chol 1 18|0.903|0.903|1.8e-009|1.5e-008|1.0e-004| 2.500507e-001 2.499497e-001| 0:0:01|chol 1 19|1.000|1.000|4.9e-010|3.5e-010|2.9e-005| 2.500143e-001 2.499857e-001|

23、 0:0:01| chol 1 110|0.904|0.904|5.7e-011|1.3e-010|4.4e-006| 2.500022e-001 2.499978e-001| 0:0:01| chol 2 211|1.000|1.000|5.2e-013|1.1e-011|1.2e-006| 2.500006e-001 2.499994e-001| 0:0:01| chol 2 212|1.000|1.000|5.9e-013|1.0e-012|1.8e-007| 2.500001e-001 2.499999e-001| 0:0:01| chol 2 213|1.000|1.000|1.7e

24、-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 = 14primal objective value = 2.50000004e-001dual objective value = 2.49999996e-0

25、01gap := 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(y), norm(Z) = 3.2e+000, 1.5e+000, 1.9e+000norm(A), norm(b), norm(

26、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 -3;b=2;5;6;C=2 1 1;1 2 3;2 2 1;cvx_beginvariable x

27、(n)minimize( a*x)subject toC * x =0cvx_end运行结果：Calling SDPT3 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 10.000 1it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime

28、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.622375e+000|0:0:00| chol 1 15|0.995|0.962|1.6e-009|6.2e-006|1.5e-002|-5.394782e+000 -