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

1、大连理工优化方法大作业MATLAB编程fun ctio n x,dk,k=fjqx(x,s)flag=0;a=0;b=0;k=0;d=1;while (flag=0)p,q=getpq(x,d,s);if (P=0) &( q=0)dk=d;x=x+d*s;flag=1;endk=k+1;if (p=0)&(q0)a=d;d=min2*d,(d+b)/2;endend%定义求函数值的函数 fun ，当输入为 x0= (x1 ， x2 )时，输出为 f function f=fun(x)f=(x(2)-x(1)A2)A2+(1-x(1)F2;function gf=gfun(x)gf=-4*x(

2、1)*(x (2) -x(1)A2)+2*(x(1)-1),2*(x(2)-x(1)A2);function p,q=getpq(x,d,s)p=fun(x)-fun(x+d*s)+0.20*d*gfun(x)*s;q=gfun(x+d*s)*s-0.60*gfun(x)*s;结果：x=0,1;s=-1,1;x,dk,k=fjqx(x,s)x =-0.0000 1.0000dk =1.1102e-016k =54习题2取初始= (0.0. 0,0)rl用兵柜梯皮法求解下面无约東优化 问题：min f (x) = x孑2xX2 十 2x孑 + x孑 H-爲X2天3 十 2xj + 3|X2 *3

3、,其中步长g的选取可利用习題1戎精确一维披索.注:通过比习题验证共範梯度法求辉门无二次西数极小点至多需 要“次迭代.fun ctio n f= fun( X )%所求问题目标函数f=X(1)A2-2*X(1)*X (2)+2*X(2)A2+X(3)A2+ X(4) A2-X( 2)*X(3)+2*X(1)+3*X(2)-X(3);end function g= gfun( X )%所求问题目标函数梯度g=2*X(1)-2*X(2)+2,-2*X(1)+4*X(2)-X(3)+3,2*X (3) -X (2)-1,2*X(4);end function x,val,k = frcg( fun,g

4、fun,xO )%功能：用FR共轭梯度法求无约束问题最小值%输入：x0是初始点，fun和gfun分别是目标函数和梯度%输出：x、val分别是最优点和最优值， k是迭代次数maxk=5000; %最大迭代次数rho=0.5;sigma=0.4;k=0;eps=10e-6;n=length(x0);while (k=0.0)d=-g;endendif (norm(g)eps)break ;endm=0;mk=0;while (m20)if(feval(fu n,xO+rhoAm*d) x0=0,0,0,0; x,val,k = frcg( fun,gfun,x0 )x =-4.0000 -3.00

5、00 -1.0000 0val =-8.0000k =或者function x,f,k=second(x)k=0;dk=dfun(x);g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);while (norm(g1)=0.02)if (k=3)k=0;g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);else if (k3)u=( norm(g1)A2)/( norm(gO)A2);s=-g1+u*s;k=k+1;g0=g1;dk=dfun(x);x=x+dk*s;g1=gfun(x);endendf=fun(x);endfunction f=f

6、un(x)f=x(1F2-2*x(1)*x (2)+2*x (2)A2+x(3)A2+x(4)A2-x (2) *x (3)+2*x(1)+3*x(2)-x(3);function gf=gfun(x)gf=2*x(1)-2*x(2)+2,-2*x(1)+4*x(2)-x(3)+3,2*x(3)-x(2)-1,2*x(4);function p,q=con(x,d)ss=-gfun(x);p=fun(x)-fun(x+d*ss)+0.2*d*gfun(x)*(ss);q=gfun(x+d*ss)*(ss)-0.6*gfun(x)*(ss);function dk=dfun(x)flag=0;a

7、=0;d=1;while (flag=0)p,q=con(x,d);if (p=0)&(q=0)dk=d;flag=1;endif (p=0)&(q x,f,k=second(x)-1.0090 0x =-4.0147 -3.0132f = -7.9999k = 1取初始点3 = （0）二考虑下面无约東优化问题：min f（x）二冷 + 2x2 + exp（xf + 天孑），其中歩长Qk的选取可別用习题1或精确一维搜索搜索方向 为一 HNW取垃=b取皿=R2f防）9耳丈啟为BFG5公式亠通过此习题体会上述三种算法的收敛速度.fun ctio n f,x,k=third_1(x) k=0;g=g

8、fu n(x);while (norm(g)=0.001) s=-g;dk=dfu n( x,s);x=x+dk*s;k=k+1;g=gfu n(x);f=fun( x);endfun ctio n f=fun(x)f=x(1)+2*x(2)A2+exp(x(1)A2+x(2)A2);fun ctio n gf=gfu n(x)gf=1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2); function j_1,j_2=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s); j_2=g

9、fun(x+d*s)*(s)-0.5*gfun(x)*(s); function dk=dfun(x,s) % 获取步长 flag=0;a=0;d=1;while (flag=0)p,q=con(x,d,s);if (p=0)&(q=0)dk=d;flag=1;endif (p=0)&(q=0.001)s=(-H*g);dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)A2+exp(x(1F2+x(2)A2);function gf=gfun(x) gf=1+2*x(1)*exp(x(

10、1F2+x(2)A2),4*x(2)+2*x(2)*(x(1F2+x(2)A2);function ggf=ggfun(x)ggf=(4*x(1)A2+2)*exp(x(1)A2+x (2) A2),4*x(1)*x (2) *exp(x(1)A2+x(2)A2); 4*x(1)*x(2)*exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2);function j_1,j_2=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s);j_2=gfun(x+d*s)*(s)-0.5*gfun(x)*(s

11、);function dk=dfun(x,s) % 步长获取flag=0;a=0;d=1;b=10000;while (flag=0)p,q=con(x,d,s);if (p=0)&(q=0)dk=d;flag=1;endif (p=0)&(q=(d+b)/2d=(d+b)/2;else d=2*d;endendEnd结果：x=0,1;f,x,k=third_2(x)f =0.7729x = -0.4193 0.0001k =8(2) 程序：function f,x,k=third_3(x)k=0;X=cell(2);g=cell(2);X1=x;H=eye(2);g1=gfun(X1);s=

12、(-H*g1);dk=dfun(X1,s);X2=X1+dk*s;g2=gfun(X2);while (norm(g2)=0.001)dx=X2-X1;dg=g2-g1;v=dx/(dx*dg)-(H*dg)/(dg*H*dg);h1=H*dg*dg*H/(dg*H*dg);h2=dx*dx/(dx*dx);h3=dg*H*dg*v*v;H=H-h1+h2+h3;k=k+1;X1=X2;g1=gfun(X1);s=(-H*g1);dk=dfun(X1,s);X2=X1+dk*s;g2=gfun(X2);norm(g2);f=fun(x);x=X2;endfunction f=fun(x)f=x

13、(1)+2*x(2)A2+exp(x(1F2+x(2)A2);function gf=gfun(x)gf=1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2);function ggf=ggfun(x)ggf=(4*x(1)A2+2)*exp(x(1)A2+x(2)A2),4*x(1)*x(2)*exp(x(1)A2+x(2)A2);4*x(1)*x(2)* exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2);function p,q=con(x,d,s)p=fun(x)-fun

14、(x+d*s)+0.1*d*gfun(x)*(s);q=gfun(x+d*s)*(s)-0.5*gfun(x)*(s);function dk=dfun(x,s)flag=0;a=0;d=1;b=10000;while (flag=0)p,q=con(x,d,s);if (p=0)&(q=0)dk=d;flag=1;endif (p=0)&(q=(d+b)/2d=(d+b)/2;else d=2*d;endendend结果：x=0,1;f,x,k=third_3(x)f =0.77290.0000x = -0.4195k=6习题4*U用有效集法求解下面勺勺二次规划问题:mmSi.(XI 一 I

15、)2 + (x2 一 2.5)2X1 - 2X2 + 2 0-Xi 2(2 + 6 0-Xi + 2X2 + 2 0xi,x2 0function callqpactH=2 0; 0 2;c=-2 -5;Ae= ; be=;Ai=1 -2; -1 -2; -1 2;1 0;0 1;bi=-2 -6 -2 0 0;x0=0 0;x,lambda,exitflag,output=qpact(H,c,Ae,be,Ai,bi,xO)fun ctio n x,lamk,exitflag,output=qpact(H,c,Ae,be,Ai,bi,x0)epsilo n=1.0e-9; err=1.0e-6

16、;k=0; x=x0; n=len gth(x); kmax=1.0e3;n e=le ngth(be); ni=le ngth(bi); lamk=zeros( ne+n i,1);in dex=ones(n i,1);for (i=1:ni)if(Ai(i,:)*xbi(i)+epsil on), i ndex(i)=0; end while (k0), Aee=Ae; endfor (j=1:ni)if (index(j)0), Aee=Aee; Ai(j,:); endendgk=H*x+c;m1,n1 = size(Aee);dk,lamk=qsubp(H,gk,Aee,zeros(

17、m1,1);if (norm(dk)ne)y,jk=min(lamk(ne+1:length(lamk);endif (y=0)exitflag=0;elseexitflag=1;for (i=1:ni)if (index(i)&(ne+sum(index(1:i)=jk)index(i)=0; break ;endendk=k+1;elseexitflag=1;alpha=1.0; tm=1.0;for (i=1:ni)if (index(i)=0)&(Ai(i,:)*dk0)tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if (tm1tm)tm=tm1; ti=i;

18、endendendalpha=min(alpha,tm);x=x+alpha*dk;if (tm0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae)*rb;x=ginvH*(Ae*lambda-c);elsex=-ginvH*c;lambda=0;end结果callqpactx =1.40001.7000lambda =0.8000exitflag =0output =fval: -6.4500iter: 7function x,mu,lambda,output=multphr(fu n, hf,gf,dfu n, dhf,dgf,xO)%功能：用乘子法解

19、一般约束问题 ：min f(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;i nk=0;epsilo n=0.00001;x=xO;he=feval(hf,

20、x);gi=feval(gf,x);n=len gth(x);l=le ngth(he);m=le ngth(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while (btakepsilon&kepsilonif k=2&btaktheta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endlambda(i)=max(0,lambda(i)-c*gi(i);endk=k+1;btaold=btak;x0=x;endendf=feval(fun,x);output.fval=f;out

21、put.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:lpsi=psi-he(i)*mu(i);s仁 s1+he(y2;psi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i);s2=s2+s3A2-lambda(i)A2;endpsi=psi+s2/(2*c);%

22、 不等式约束函数文件 g1.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-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,

23、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 ;enddk=-Bkgk;m=0;mk=0;while (m20)n ewf=feval(fu n, x0+rhoAm*dk,vararg in :);oldf=feval(fun,x0,varargin:);if(newf0)Bk=Bk-(Bk*sk*sk*Bk)/(sk*Bk*sk)+(yk

24、*yk)/(yk*sk);endk=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.00134.8987mu =0.7701lambda =000.9434output =fval: -31.9923iter: 4习题6利用序列二次规划方法求解习题5中的约束优化问题:min 4xi 一好一 12s.t. 25 - x? x孑=Q10xt 一 召 + 10旳-xj - 34 0 X1,X2 0f=3,1,1；A=2,1,1;1,-1,-1;b=2;-1;lb=0,0,0； x=li nprog(f,A,b,zeros(3),0,0,0,lb)结果：Optimization terminated.0.00000.50000.5000