大连理工优化方法大作业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)&(q=0)dk=d;x=x+d*s;flag=1;end k=k+1; if(p=0)&(q0)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=0.0) d=-g; end end if(norm(g)eps) break; end m=0;mk=0; while(m20) if(feval(fun,x0+rhom*d) 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(k3)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)&(q=0) dk=d;flag=1;end if(p=0)&(q 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)&(q=0)dk=d;flag=1;end if(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)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)&(q=0) dk=d;flag=1;end if(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=(-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)&(q=0) dk=d;flag=1;end if(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.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,:)*xbi(i)+epsilon), index(i)=0; endendwhile(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(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;endendendk=k+1;else exitflag=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;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.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(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.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