大连理工优化方法大作业MATLAB编程.docx
- 文档编号:9960345
- 上传时间:2023-02-07
- 格式:DOCX
- 页数:23
- 大小:217.63KB
大连理工优化方法大作业MATLAB编程.docx
《大连理工优化方法大作业MATLAB编程.docx》由会员分享,可在线阅读,更多相关《大连理工优化方法大作业MATLAB编程.docx(23页珍藏版)》请在冰豆网上搜索。
大连理工优化方法大作业MATLAB编程
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=min{2*d,(d+b)/2};
end
end
%定义求函数值的函数fun,当输入为x0=(x1,x2)时,输出为f
functionf=fun(x)
f=(x
(2)-x
(1)^2)^2+(1-x
(1))^2;
functiongf=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.00001.0000
dk=1.1102e-016
k=54
functionf=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);
end
functiong=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,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 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) break; end m=0;mk=0; while(m<20) if(feval(fun,x0+rho^m*d) mk=m;break; end m=m+1; end x0=x0+rho^mk*d; val=feval(fun,x0); g0=g;d0=d; k=k+1; end x=x0; val=feval(fun,x0); end 结果: >>x0=[0,0,0,0]; >>[x,val,k]=frcg('fun','gfun',x0) x= -4.0000-3.0000-1.00000 val= -8.0000 k= 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); elseif(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); end end f=fun(x); end functionf=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); functiongf=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)'; functiondk=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=min{2*d,(d+b)/2}; end End 结果: x=[0,0,0,0]; >>[x,f,k]=second(x) x=-4.0147-3.0132-1.00900 f=-7.9999 k=1 function[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); end functionf=fun(x) f=x (1)+2*x (2)^2+exp(x (1)^2+x (2)^2); functiongf=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)'; functiondk=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=min{2*d,(d+b)/2}; end end 结果: x=[0,1]; [f,x,k]=third_1(x) f=0.7729 x=-0.41960.0001 k=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); end functionf=fun(x) f=x (1)+2*x (2)^2+exp(x (1)^2+x (2)^2); functiongf=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)]; functionggf=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)'; functiondk=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; if2*d>=(d+b)/2 d=(d+b)/2; elsed=2*d; end end End 结果: x=[0,1]; [f,x,k]=third_2(x) f=0.7729 x=-0.41930.0001 k=8 (2)程序: function[f,x,k]=third_3(x) k=0; X=cell (2); g=cell (2); X{1}=x; H=eye (2); g{1}=gfun(X{1}); s=(-H*g{1}')'; dk=dfun(X{1},s); X{2}=X{1}+dk*s; g{2}=gfun(X{2}); while(norm(g{2})>=0.001) dx=X{2}-X{1}; dg=g{2}-g{1}; 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; X{1}=X{2}; g{1}=gfun(X{1}); s=(-H*g{1}')'; dk=dfun(X{1},s); X{2}=X{1}+dk*s; g{2}=gfun(X{2}); norm(g{2}); f=fun(x); x=X{2}; end functionf=fun(x) f=x (1)+2*x (2)^2+exp(x (1)^2+x (2)^2); functiongf=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)]; functionggf=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)'; functiondk=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; if2*d>=(d+b)/2 d=(d+b)/2; elsed=2*d; end end end 结果: x=[0,1]; [f,x,k]=third_3(x) f=0.7729 x=-0.41950.0000 k=6 functioncallqpact H=[20;02]; c=[-2-5]'; Ae=[];be=[]; Ai=[1-2;-1-2;-12;10;01]; bi=[-2-6-200]'; x0=[00]'; [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;end end while(k<=kmax) Aee=[]; if(ne>0),Aee=Ae;end for(j=1: ni) if(index(j)>0),Aee=[Aee;Ai(j,: )];end end gk=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))); end if(y>=0) exitflag=0; else exitflag=1; for(i=1: ni) if(index(i)&(ne+sum(index(1: i)))==jk) index(i)=0;break; end end end k=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=tm1;ti=i; end end end alpha=min(alpha,tm); x=x+alpha*dk; if(tm<1),index(ti)=1;end end if(exitflag==0),break;end k=k+1; end output.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); else x=-ginvH*c; lambda=0; end 结果 >>callqpact x= 1.4000 1.7000 lambda= 0.8000 exitflag= 0 output= fval: -6.4500 iter: 7 function[x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0) %功能: 用乘子法解一般约束问题: minf(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 %调用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; fori=1: l btak=btak+he(i)^2; end %更新乘子向量 fori=1: m temp=min(gi(i),lambda(i)/c); btak=btak+temp^2; end btak=sqrt(btak); ifbtak>epsilon ifk>=2&&btak>theta*btaold c=eta*c; end fori=1: l mu(i)=mu(i)-c*he(i); end fori=1: m lambda(i)=max(0,lambda(i)-c*gi(i)); end k=k+1; btaold=btak; x0=x; end end f=feval(fun,x); output.fval=f; output.iter=k; %增广拉格朗日函数 functionpsi=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; fori=1: l psi=psi-he(i)*mu(i); s1=s1+he(i)^2; end psi=psi+0.5*c*s1; s2=0; fori=1: m s3=max(0,lambda(i)-c*gi(i)); s2=s2+s3^2-lambda(i)^2; end psi=psi+s2/(2*c); %不等式约束函数文件g1.m functiongi=g1(x) gi=10*x (1)-x (1)^2+10*x (2)-x (2)^2-34; %目标函数的梯度文件df1.m functiong=df1(x) g=[4,-2*x (2)]'; %等式约束(向量)函数的Jacobi矩阵(转置)文件dh1.m functiondhe=dh1(x) dhe=[-2*x (1),-2*x (2)]' %不等式约束(向量)函数的Jacobi矩阵(转置)文件dg1.m functiondgi=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 gk=feval(gfun,x0,varargin{: }); if(norm(gk) break; end dk=-Bk\gk; m=0;mk=0; while(m<20) newf=feval(fun,x0+rho^m*dk,varargin{: }); oldf=feval(fun,x0,varargin{: }); if(newf mk=m; break; end m=m+1; end x=x0+rho^mk*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; end val=feval(fun,x0,varargin{: }); 结果 x=[22]'; [x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0) x= 1.0013 4.8987 mu= 0.7701 lambda= 0 0 0.9434 output= fval: -31.9923 iter: 4 f=[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) 结果: Optimizationterminated. x= 0.0000 0.5000 0.5000
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工 优化 方法 作业 MATLAB 编程