优化方法MATLAB编程 大连理工大学.docx

1、优化方法MATLAB编程 大连理工大学优化方法上机大作业 学 院： 姓 名： 学 号： 指导老师： 肖 现 涛 第一题源程序如下：function zy_x = di1ti(x)%di1ti是用来求解优化作业第一题的函数。x0=x; yimuxulong=0.000001;g0=g(x0);s0=-g0;A=2*ones(100,100);k=0;while k100 lanmed=-(g0)*s0/(s0*A*s0); x=x0+lanmed*s0; g=g(x); k=k+1; if norm(g) x=zeros(100,1); di1ti(x)After 1 iterations,ob

2、tain the optimal solution. The optimal solution is -0.250000.The optimal x is ans.ans =0.005*ones(100,1).第二题1. 最速下降法。源程序如下：function zy_x=di2titidu(x)%该函数用来解大作业第二题。x0=x; yimuxulong=1e-5; k=0;g0=g(x0); s0=-g0;while k=0 if norm(g0)=0&i(f(x0)+c*lanmed*g0*s0) lanmed=lanmed/2; i=i+1; else break; end end x

3、=x0+lanmed*s0; x0=x; g0=g(x); s0=-g0; k=k+1; endendzy_x=x;zyj=f(x);fprintf(after %d iterations,obtain the optimal solution.nnThe optimal solution is %f.nn The optimal x is ans.n,k,zyj); function f=f(x)x1=1 0 0 0*x;x2=0 1 0 0*x;x3=0 0 1 0*x;x4=0 0 0 1*x;f=(x1-1)2+(x3-1)2+100*(x2-x12)2+100*(x4-x32)2;f

4、unction g=g(x)x1=1 0 0 0*x;x2=0 1 0 0*x;x3=0 0 1 0*x;x4=0 0 0 1*x;g=2*(x1-1)-400*x1*(x2-x12);200*(x2-x12);2*(x3-1)-400*x3*(x4-x32);200*(x4-x32); x=-1.2 1 -1.2 1; di2titidu(x)after 5945 iterations,obtain the optimal solution.The optimal solution is 0.000000. The optimal x is ans.ans = 1.0000 1.0000 1

5、.00001.00002. 牛顿法源程序如下：function zy_x=di2tinewton(x)%该函数用来解大作业第二题。x0=x; yimuxulong=1e-5; k=0;g0=g(x0); h0=h(x0);s0=-inv(h0)*g0;while k=0&k1000 if norm(g0) x=-1.2 1 -1.2 1; di2tinewton(x)after 6 iterations,obtain the optimal solution.The optimal solution is 0.000000. The optimal x is ans.ans = 1.0000

6、1.0000 1.00001.0000可以看出，用Newton法经过6次迭代就能求出最优解。3. BFGS法源程序如下：function zy_x=di2tiBFGS(x)%该函数用来解大作业第二题。x0=x; yimuxulong=1e-5; k=0;g0=g(x0); H0=eye(4);s0=-H0*g0;while k=0&k100 if norm(g0)=0&i(f(x0)+c*lanmed*g0*s0) lanmed=lanmed/2; i=i+1; else break; end end x=x0+lanmed*s0; dete_x=x-x0; dete_g=g(x)-g0; m

7、iu=1+dete_g*H0*dete_g/(dete_x*dete_g); H=H0+(miu*dete_x*dete_x-H0*dete_g*dete_x-dete_x*dete_g*H0)/(dete_x*dete_g); s=-H*g(x); x0=x; s0=s; H0=H; g0=g(x); k=k+1; endendzy_x=x;zyj=f(x);fprintf(after %d iterations,obtain the optimal solution.nnThe optimal solution is %f.nn The optimal x is ans.n,k,zyj);

8、 function f=f(x)x1=1 0 0 0*x;x2=0 1 0 0*x;x3=0 0 1 0*x;x4=0 0 0 1*x;f=(x1-1)2+(x3-1)2+100*(x2-x12)2+100*(x4-x32)2;function g=g(x)x1=1 0 0 0*x;x2=0 1 0 0*x;x3=0 0 1 0*x;x4=0 0 0 1*x;g=2*(x1-1)-400*x1*(x2-x12);200*(x2-x12);2*(x3-1)-400*x3*(x4-x32);200*(x4-x32);代入初始值，计算结果如下： x=-1.2 1 -1.2 1; di2tiBFGS(

9、x)after 53 iterations,obtain the optimal solution.The optimal solution is 0.000000. The optimal x is ans.ans = 1.0000 1.0000 1.00001.0000第三题1. 惩罚函数法源程序如下：function zy_x=di3ti(x)%该函数用来解大作业第三题。x0=x; M=100; c=4; m=1;while m0 g0=g(x0,M); yimuxulong=1e-5;k=0;s0=-inv(H(x0,M)*g0; while k=0 if norm(g0)yimuxu

10、long break; else x=x0+s0; %牛顿法； x0=x; g0=g(x,M); s0=-inv(H(x0,M)*g0; k=k+1; end end if max(abs(h(x),g1(x),g2(x),g3(x)0.5 break; else M=M*c; m=m+1; endendzy_x=x;zyj=f(x); fprintf(after %d iterations,obtain the optimal solution.nnThe optimal solution is %f.nn The optimal x is ans.n,m,zyj); function F=

11、F(x,M) x1=1 0*x;x2=0 1*x;F=4*x1-x22-12+M*(h2+g12+g22+g32);function g=g(x,M) x1=1 0*x;x2=0 1*x;g=4+M*(-4*(25-x12-x22)*x1+2*(10*x1-x12+10*x2-x22-34)*(10-2*x1)+2*x1);-2*x2+M*(-4*(25-x12-x22)*x2+2*(10*x1-x12+10*x2-x22-34)*(10-2*x2)+2*x2);function H=H(x,M)x1=1 0*x;x2=0 1*x;H=M*(24*x12-120*x1+8*x22-40*x2+

12、238),M*(16*x1*x2-40*x1-40*x2+200);M*(16*x1*x2-40*x1-40*x2+200),-2+M*(24*x22-120*x2+8*x12-40*x1+238);function f=f(x)x1=1 0*x;x2=0 1*x;f=4*x1-x22-12;function h=h(x)x1=1 0*x;x2=0 1*x;h=25-x12-x22;function g1=g1(x)x1=1 0*x;x2=0 1*x;g=10*x1-x12+10*x2-x22-34;if g=0 g2=0;else g2=x1;endfunction g3=g3(x)x1=1

13、 0*x;x2=0 1*x;if x2=0 g3=0;else g3=x2;end代入任意初始值，运算结果如下。 x=rand(2,1); di3ticf(x)after 1 iterations,obtain the optimal solution.The optimal solution is -31.490552. The optimal x is ans.ans = 1.00244.84772. 乘子法源程序如下：function x,mu,lambda,output=multphr(fun,hf,gf,dfun,dhf,dgf,x0)%功能: 用乘子法解一般约束问题: min f(

14、x), s.t. h(x)=0, g(x).=0%输入: x0是初始点, fun, dfun分别是目标函数及其梯度；% hf, dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；% gf, dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；%输出: x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量; % output是结构变量, 输出近似极小值f, 迭代次数, 内迭代次数等maxk=500;c=2.0;eta=2.0;theta=0.8;k=0;ink=0;epsilon=0.00001;x=x0;he=feval(hf,x);gi=feva

15、l(gf,x);n=length(x);l=length(he);m=length(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while(btakepsilon&kepsilon if k=2&btaktheta*btaold c=eta*c; end for i=1:l mu(i)=mu(i)-c*he(i); end for i=1:m lambda(i)=max(0,lambda(i)-c*gi(i); end k=k+1; btaold=btak; x0=x; endendf=feval(fun,x);output.fv

16、al=f;output.iter=k;%增广拉格朗日函数function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c)f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for i=1:l psi=psi-he(i)*mu(i); s1=s1+he(i)2;endpsi=psi+0.5*c*s1;s2=0;for i=1:m s3=max(0,lambda(i)-c*gi(i); s2=s2+s32-lambda(i)2;endpsi=

17、psi+s2/(2*c);%不等式约束函数文件g1.mfunction gi=g1(x)gi=10*x(1)-x(1)2+10*x(2)-x(2)2-34;%目标函数的梯度文件df1.mfunction g=df1(x)g=4, -2*x(2);%等式约束（向量）函数的Jacobi矩阵（转置）文件dh1.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);function x,val,k=bfgs(fun,gfun,

18、x0,varargin)maxk=500;rho=0.55;sigma=0.4;epsilon=0.00001;k=0;n=length(x0);Bk=eye(n);while(kmaxk) gk=feval(gfun,x0,varargin:); if(norm(gk)epsilon) break; end dk=-Bkgk; m=0;mk=0; while(m20) newf=feval(fun,x0+rhom*dk,varargin:); oldf=feval(fun,x0,varargin:); if(newf0) Bk=Bk-(Bk*sk*sk*Bk)/(sk*Bk*sk)+(yk*yk)/(yk*sk); end k=k+1;x0=x;endval=feval(fun,x0,varargin:);结果x=2 2;x,mu,lambda,output=multphr(fun,hf,gf1,df,dh,dg,x0)x = 1.0013 4.8987mu = 0.7701lambda = 0 0 0.9434output = fval: -31.9923iter: 4