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

精选优质文档-倾情为你奉上function 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) b=d;d=(d+a)/2; end if(p>=0)&&(q>=0)dk=d;x=x+d*s;flag=1;end k=k+1; if(p>=0)&&(q<0)a=d;d=min2*d,(d+b)/2;endend%定义求函数值的函数fun，当输入为x0=（x1，x2）时，输出为ffunction f=fun(x) f=(x(2)-x(1)2)2+(1-x(1)2;function gf=gfun(x) gf=-4*x(1)*(x(2)-x(1)2)+2*(x(1)-1),2*(x(2)-x(1)2);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 =54function f= fun( X )%所求问题目标函数f=X(1)2-2*X(1)*X(2)+2*X(2)2+X(3)2+ X(4)2- X(2)*X(3)+2*X(1)+3*X(2)-X(3);endfunction 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);endfunction x,val,k = frcg( fun,gfun,x0 )%功能：用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<maxk) g=feval(gfun,x0);%计算梯度 itern=k-(n+1)*floor(k/(n+1); itern=itern+1; %计算搜索方向 if(itern=1) d=-g; else beta=(g*g')/(g0*g0'); d=-g+beta*d0; gd=g'*d; if(gd>=0.0) d=-g; end end if(norm(g)<eps) break; end m=0;mk=0; while(m<20) if(feval(fun,x0+rhom*d)<feval(fun,x0)+sigma*rhom*g'*d) mk=m;break; end m=m+1; end x0=x0+rhomk*d; val=feval(fun,x0); g0=g;d0=d; k=k+1;endx=x0;val=feval(fun,x0);end结果：>> x0=0,0,0,0;>> x,val,k = frcg( 'fun','gfun',x0 )x = -4.0000 -3.0000 -1.0000 0val = -8.0000k =21或者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(k<3)u=(norm(g1)2)/(norm(g0)2); s=-g1+u*s;k=k+1;g0=g1;dk=dfun(x);x=x+dk*s;g1=gfun(x);endendf=fun(x);endfunction f=fun(x) f=x(1)2-2*x(1)*x(2)+2*x(2)2+x(3)2+x(4)2-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=0;d=1;while(flag=0)p,q=con(x,d);if (p<0) b=d;d=(d+a)/2; end if(p>=0)&&(q>=0) dk=d;flag=1;end if(p>=0)&&(q<0) a=d;d=min2*d,(d+b)/2;endEnd结果：x=0,0,0,0;>> x,f,k=second(x)x =-4.0147 -3.0132 -1.0090 0f = -7.9999k = 1function f,x,k=third_1(x) k=0; g=gfun(x);while(norm(g)>=0.001) s=-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)2+exp(x(1)2+x(2)2);function gf=gfun(x)gf=1+2*x(1)*exp(x(1)2+x(2)2),4*x(2)+2*x(2)*(x(1)2+x(2)2);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)'function dk=dfun(x,s)%获取步长flag=0;a=0;d=1;while(flag=0)p,q=con(x,d,s);if (p<0) b=d;d=(d+a)/2; end if(p>=0)&&(q>=0)dk=d;flag=1;end if(p>=0)&&(q<0) a=d;d=min2*d,(d+b)/2;endend结果：x=0,1; f,x,k=third_1(x) f =0.7729x = -0.4196 0.0001k =8（1） 程序：function f,x,k=third_2(x)k=0; H=inv(ggfun(x);g=gfun(x);while(norm(g)>=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)2+exp(x(1)2+x(2)2);function gf=gfun(x)gf=1+2*x(1)*exp(x(1)2+x(2)2),4*x(2)+2*x(2)*(x(1)2+x(2)2);function ggf=ggfun(x)ggf=(4*x(1)2+2)*exp(x(1)2+x(2)2),4*x(1)*x(2)*exp(x(1)2+x(2)2); 4*x(1)*x(2)*exp(x(1)2+x(2)2),4+(4*x(2)2+2)*exp(x(1)2+x(2)2);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)'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) b=d;d=(d+a)/2; end if(p>=0)&&(q>=0) dk=d;flag=1;end if(p>=0)&&(q<0) a=d;if 2*d>=(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=(-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(1)+2*x(2)2+exp(x(1)2+x(2)2);function gf=gfun(x)gf=1+2*x(1)*exp(x(1)2+x(2)2),4*x(2)+2*x(2)*(x(1)2+x(2)2);function ggf=ggfun(x)ggf=(4*x(1)2+2)*exp(x(1)2+x(2)2),4*x(1)*x(2)*exp(x(1)2+x(2)2);4*x(1)*x(2)*exp(x(1)2+x(2)2),4+(4*x(2)2+2)*exp(x(1)2+x(2)2);function p,q=con(x,d,s) p=fun(x)-fun(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) b=d;d=(d+a)/2; end if(p>=0)&&(q>=0) dk=d;flag=1;end if(p>=0)&&(q<0) a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendend结果：x=0,1;f,x,k=third_3(x)f =0.7729x = -0.4195 0.0000k=6function 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,x0)function x,lamk,exitflag,output=qpact(H,c,Ae,be,Ai,bi,x0)epsilon=1.0e-9; err=1.0e-6;k=0; x=x0; n=length(x); kmax=1.0e3;ne=length(be); ni=length(bi); lamk=zeros(ne+ni,1);index=ones(ni,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsilon), index(i)=0; endendwhile(k<=kmax)Aee=;if(ne>0), 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(m1,1);if(norm(dk)<=err)y=0.0;if(length(lamk)>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;endendendk=k+1;else exitflag=1;alpha=1.0; tm=1.0;for(i=1:ni)if(index(i)=0)&(Ai(i,:)*dk<0)tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if(tm1<tm)tm=tm1; ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1), index(ti)=1; endendif(exitflag=0), break; endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;function x,lambda=qsubp(H,c,Ae,be)ginvH=pinv(H);m,n=size(Ae);if(m>0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;end结果>>callqpactx = 1.4000 1.7000lambda = 0.8000exitflag = 0output = fval: -6.4500iter: 7function x,mu,lambda,output=multphr(fun,hf,gf,dfun,dhf,dgf,x0)%功能: 用乘子法解一般约束问题: 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;ink=0;epsilon=0.00001;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while(btak>epsilon&&k<maxk)%调用BFGS算法程序求解无约束子问题 x,ival,ik=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c); ink=ink+ik; he=feval(hf,x);gi=feval(gf,x); btak=0; for i=1:l btak=btak+he(i)2;end%更新乘子向量 for i=1:m temp=min(gi(i),lambda(i)/c); btak=btak+temp2; end btak=sqrt(btak); if btak>epsilon if k>=2&&btak>theta*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.fval=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=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,x0,varargin)maxk=500;rho=0.55;sigma=0.4;epsilon=0.00001;k=0;n=length(x0);Bk=eye(n);while(k<maxk) gk=feval(gfun,x0,varargin:); if(norm(gk)<epsilon) break; end dk=-Bkgk; m=0;mk=0; while(m<20) newf=feval(fun,x0+rhom*dk,varargin:); oldf=feval(fun,x0,varargin:); if(newf<oldf+sigma*rhom*gk'*dk) mk=m; break; end m=m+1; end x=x0+rhomk*dk; sk=x-x0; yk=feval(gfun,x,varargin:)-gk; if(yk'*sk>0) 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: 4f=3,1,1;A=2,1,1;1,-1,-1;b=2;-1;lb=0,0,0;x=linprog(f,A,b,zeros(3),0,0,0',lb)结果：Optimization terminated.x = 0.0000 0.5000 0.5000专心-专注-专业