最优化理论 附录C Matlab算法示例

CMatlab程序Matlab程序.§C.1一维搜索的算例确定初始搜索区间的进退法（3.1）程序示例一function[a,b,c]=JintuifaI(phi,xk,dk,alpha0,h0)%该函数是一维搜索过程中确定搜索区间的进退法算法之一phi目标函数xk代表当前迭代点dk当前搜索方向alpha0初始点h0初始步长a搜索区间的左端点b搜索区间的右端点eps=1.0e-6;alpha1=alpha0;phi1=phi(xk+alpha1*dk);k=0;h=h0;while1k=k+1;tt=alpha1+h;ifphi_tmp<phi1alpha2=alpha1;phi2=phi1;alpha1=tt;h=2*h;k=k+1;elseifk==1alpha2=tt;phi2=phi_tmp;h=-h;elsea=min(tt,alpha2);b=max(tt,alpha2);c=alpha1;breakendendend上面的函数可以通过如下主程序调用.clear;clc;n=2;x=ones(n,1);xx=ones(n,1);dd=floor(4*rand(n,1)+1);alpha=0;h=0.1;[a,b,c]=JintuifaI(f,xx,dd,alpha,h)示例二function[a,b,c]=JintuifaII(f,n,xk,dk,alpha0,h0)%该函数是一维搜索过程中确定搜索区间的进退法算法之二eps=1.0e-6;alpha1=alpha0;phi1=eval(ff)k=0;h=h0;while1alpha_tmp=alpha1+h;k=k+1;phi_tmp=eval(ff);ifphi_tmp<phi1alpha2=alpha1;phi2=phi1;alpha1=alpha_tmp;phi1=phi_tmp;h=2*h;k=k+1;elseifk==1phi2=phi_tmp;h=-h;elsea=min(alpha_tmp,alpha2)c=alpha1breakendendendfunction f=3;fori=1:nf=f+x(i)^2;end调用这个函数的主程序为clear;clc;n=2;xx=ones(n,1)dd=floor(4*rand(n,1)+1)t0=0;h=0.1;[a,b,c]=JintuifaII('fun31',n,xx,dd,t0,h);示例三function[a,b,c]=JintuifaIII(f,n,xk,dk,alpha0,h0)%该函数是一维搜索过程中确定搜索区间的进退法算法之三eps=1.0e-6;symst;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f);alpha1=alpha0;%phi1=subs(phi,'t',alpha1);t=alpha1;phi1=eval(phi);k=0;h=h0;while1alpha_tmp=alpha1+h;k=k+1;% t=alpha_tmp;phi_tmp=eval(phi);ifphi_tmp<phi1alpha2=alpha1;phi2=phi1;alpha1=alpha_tmp;phi1=phi_tmp;h=2*h;k=k+1;elseifk==1phi2=phi_tmp;h=-h;elsea=min(alpha_tmp,alpha2);c=alpha1;breakendendend调用这个函数的主程序如下.clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h0=0.5;[a,b,c]=JintuifaIII(f,n,xx,dd,t0,h0)0.618法（3.2）程序function[alpha_min,phi_min,kk]=Hjfengefa(f,n,xk,dk,a,b)%该函数是一维搜索的黄金分割算法eps=1.0e-6;l=a+0.382*(b-a);u=a+0.618*(b-a);k=1;ff=[f,'(xk+l*dk,n)'];phi1=eval(ff);phi2=eval(ff);delta=eps;while1ifphi1>phi2ifb-l<=deltaalpha_min=u;phi_min=phi2;break;elsea=l;l=u;phi1=phi2;u=a+0.618*(b-a);phi2=eval(ff);k=k+1;endelseifu-a<=deltaalpha_min=l;phi_min=phi1;break;elseb=u;u=l;phi2=phi1;l=a+0.382*(b-a);phi1=eval(ff);k=k+1;endendendalpha_min=(a+b)/2;ff=[f,'(xk+alpha_min*dk,n)'];kk=k;调用这个函数的主程序为clear;clc;n=2;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h=0.1;[a,b,c]=JintuifaII('fun31',n,xx,dd,t0,h);[alpha_min,phi_min,kk]=Hjfengefa('fun31',n,xx,dd,a,b)Fabonacci法程序function[alpha_min,phi_min,kk]=Fabonacci(f,nn,xk,dk,a,b)%该函数是一维搜索的斐波那契算法eps=1.0e-6;delta=eps;N=(b-a)/delta;F=ones(2,1);n=2;whileF(n)<Nn=n+1;F(n)=F(n-1)+F(n-2);endl=a+(F(n-2)/F(n))*(b-a);u=a+(F(n-1)/F(n))*(b-a);ff=[f,'(xk+l*dk,nn)'];phi1=eval(ff);phi2=eval(ff);if phi1>phi2a=l;l=u;phi1=phi2;u=a+(F(k-1)/F(k))*(b-a);phi2=eval(ff);elseb=u;u=l;phi2=phi1;l=a+(F(k-2)/F(k))*(b-a);phi1=eval(ff);endendalpha_min=(a+b)/2;phi_min=eval(ff);kk=n-1;调用这个函数的主程序为clear;clc;n=2;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h=0.1;[a,b,c]=JintuifaII('fun31',n,xx,dd,t0,h);[alpha_min,phi_min,kk]=Fabonacci('fun31',n,xx,dd,a,b)一维搜索的二分法（3.3）程序function[alpha_min,phi_min,kk]=Erfenfa(f,n,xk,dk,a,b)%该函数是精确一维搜索的二分法算法eps=1.0e-6;symst;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f);g=diff(phi);k=0;while1alpha=(a+b)/2;t=alpha;tmp=eval(g);iftmp>0b=alpha;elsea=alpha;endk=k+1;if(delta<eps)breakendendalpha_min=(a+b)/2;t=alpha_min;phi_min=eval(phi);kk=k;调用这个函数的主程序为clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h0=0.5;[a,b,c]=JintuifaIII(f,n,xx,dd,t0,h0)[alpha_min,phi_min,kk]=Erfenfa(f,n,xx,dd,a,b)§C.1.5一维搜索的插值法程序§C.1.5.1一点二次插值法的程序function[alpha_min,phi_min,kk]=Nchazhi(f,n,xk,dk,c)%该函数是精确一维搜索的一点二次插值（牛顿经典插值）算法，适合于凸函数eps=1.0e-6;symst;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f);g=diff(phi);h=diff(g);k=0;t=c;alpha=c;while1ga=eval(g);ha=eval(h);ifabs(ga)<epsbreakendalpha=alpha-ga/ha;t=alpha;k=k+1;endalpha_min=alpha;phi_min=eval(phi);kk=k;调用这个函数的主程序如下.clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h0=0.5;[a,b,c]=JintuifaIII(f,n,xx,dd,t0,h0)[alpha_min,phi_min,kk]=Nchazhi(f,n,xx,dd,c)§C.1.5.2两点二次插值法的程序（3.4）function[alpha_min,phi_min,kk]=LiangdianI(f,n,xk,dk,alpha,h)%Ieps=1.0e-6;symst;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f);g=diff(phi);t=alpha;phi1=eval(phi);g1=eval(g);rho=0.5;alpha1=alpha;k=0;while1ifg1>0h=-h;endwhile1alpha2=alpha1+h;t=alpha2;phi2=eval(phi);ifphi2>phi1+g1*(alpha2-alpha1)breakendh=2*h;endt=alpha_tmp;phi_tmp=eval(phi);g_tmp=eval(g);ifabs(g_tmp)<epsbreakendalpha1=alpha_tmp;phi1=phi_tmp;g1=g_tmp;h=rho*h;k=k+1;endalpha_min=alpha_tmp;phi_min=phi_tmp;kk=k;调用该函数的主程序为clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h0=0.5;[alpha_min,phi_min,kk]=LiangdianI(f,n,xx,dd,t0,h0)示例二function[alpha_min,phi_min,kk]=LiangdianI(f,n,xk,dk,alpha,h)%IIeps=1.0e-6;symst;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f);g=diff(phi);alpha1=alpha;t=alpha;phi1=eval(phi);g1=eval(g);ifg1>0h=-h;endwhile1alpha2=alpha1+h;t=alpha2;g2=eval(g);ifg1*g2<=0breakendh=2*h;endk=0;while1ifdelta>0alpha_tmp=alpha1-g1*delta;elsealpha_tmp=(alpha1+alpha2)/2;endt=alpha_tmp;g_tmp=eval(g);ifabs(g_tmp)<epsbreakendifg1*g_tmp<0alpha2=alpha_tmp;g2=g_tmp;elsealpha1=alpha_tmp;g1=g_tmp;endk=k+1;endalpha_min=alpha_tmp;phi_min=eval(phi);kk=k;调用该函数的主程序为clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h0=0.5;[alpha_min,phi_min,kk]=LiangdianII(f,n,xx,dd,t0,h0)§C.1.5.3三点二次插值法（3.5）程序function[alpha_min,phi_min,kk]=SandianCZF(f,n,xk,dk,alpha1,alpha3,alpha2)%该函数是精确一维搜索的三点二次插值算法eps=1.0e-6;phi1=eval(ff);phi2=eval(ff);phi3=eval(ff);kk=0;while1ifabs(a)<epsalpha_min=alpha2;breakendif(b-alpha1)*(b-alpha3)>0alpha_min=alpha2;breakendifabs(b-alpha2)<epsalpha_min=b;breakendphib=eval(ff);if(b-alpha1)*(b-alpha2)>0ifphib<phi2alpha1=alpha2;phi1=phi2;alpha2=b;phi2=phib;elsealpha3=b;phi3=phib;endelseifphib<phi2alpha3=alpha2;phi3=phi2;alpha2=b;phi2=phib;elsealpha1=b;phi1=phib;endendkk=kk+1;endphi_min=eval(ff);调用该函数的主程序为clear;clc;n=2;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h=0.1;[a,b,c]=JintuifaII('fun31',n,xx,dd,t0,h);[alpha_min,phi_min,kk]=SandianCZF('fun31',n,xx,dd,a,b,c)不精确一维搜索（3.6）程序function[alpha_min,phi_min,kk]=ArmijoGoldstein(f,n,xk,dk,alpha,h)%该函数是不精确一维搜索的Armijo-Goldstein算法eps=1.0e-6;rho=0.1;tt=10;symst;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f);g=diff(phi);a=alpha;t=a;phi1=eval(phi);g1=eval(g);ifg1>0h=-h;enda_tmp=a+h;flag=0;kk=0;while1t=a_tmp;phi_tmp=eval(phi);tmp1=phi_tmp-phi1-rho*g1*(a_tmp-alpha);iftmp1>0b=a_tmp;flag=1;elseiftmp2<0a=a_tmpelsebreak;endifflag==0a_tmp=a+tt*h;elsea_tmp=(b+a)/2;endkk=kk+1;endt=alpha_min;phi_min=eval(phi);调用这个函数的主程序为clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h0=0.1;[alpha_min,phi_min,kk]=ArmijoGoldstein(f,n,xx,dd,t0,h0)不精确一维搜索（3.7）程序function[alpha_min,phi_min,kk]=WolfePowell(f,n,xk,dk,alpha,h)%该函数是不精确一维搜索的Wolfe-Powell算法eps=1.0e-6;rho=0.1;sigma=0.5;symsts;s=t;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f)g=diff(phi)a=alpha;t=a;phi0=eval(phi);g0=eval(g);rr=0;ifg0>0t=alpha-s;phi=eval(phi);g=diff(phi);t=0;phi0=eval(phi);g0=eval(g);rr=1;endkk=0;flag=0;a_tmp=a+h;while1t=a_tmp;phi_tmp=eval(phi);iftmp1>0b=a_tmp;flag=1;a_tmp=(b+a)/2;elseg_tmp=subs(g);% tmp2=sigma*abs(g0)-abs(g_tmp);iftmp2<0a=a_tmp;ifflag==0h=2*h;a_tmp=a+h;elsea_tmp=(b+a)/2;endelsebreak;endendkk=kk+1;endt=alpha_min;phi_min=eval(phi);ifrr==1t=alpha-alpha_min;end调用这个函数的主程序为clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);t0=0;h0=0.5;[alpha_min,phi_min,kk]=WolfePowell(f,n,xx,dd,t0,h0)§C.1.8不精确一维搜索的后退方法程序function[alpha_min,phi_min,kk]=Armijoonly(f,n,xk,dk)%该函数是不精确一维搜索的仅用简单Armijo准则的后退算法eps=1.0e-6;rho=0.1;rr=0.5;;symst;x0=xk+t*dk;fori=1:neval(['x'num2str(i)'=x0('num2str(i)');'])endphi=eval(f);g=diff(phi);t=0;phi0=eval(phi);g0=eval(g);t=1;ifg0>0t=-1;endkk=0;while1phi_tmp=eval(phi);iftmp>0t=rr*t;elsealpha_min=t;phi_min=phi_tmp;break;endkk=kk+1;end调用这个函数的主程序为clear;clc;n=2;fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;xx=ones(n,1);dd=floor(4*rand(n,1)+1);[alpha_min,phi_min,kk]=Armijoonly(f,n,xx,dd)§C.2利用导数的无约束最优化算例§C.2.1最速下降法function[x_min,f_min,kk]=Zuisuxiajiangfa(f,n,x0)%该函数是求解无约束最优化问题的最速下降算法eps=1.0e-6;fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;t0=0;h=0.1;x=x0;while1fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx=eval(gg);fx=eval(ff);ifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;enddk=-gx';% x=x+alpha*dk;k=k+1;end这里调用的函数定义为function fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;g=jacobian(f,x);调用这个函数的主程序为clear;clc;n=2;xx=3*ones(n,1);[x_min,f_min,kk]=Zuisuxiajiangfa('fun4',n,xx)§C.2.2牛顿法function[x_min,f_min,kk]=Newtonfa(f,n,x0)%该函数是求解无约束最优化问题的阻尼牛顿算法eps=1.0e-6;eps1=0.01;fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;x=x0;while1fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx=eval(gg);fx=eval(ff);hx=eval(hh);ifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;enddk=-gx';ifabs(det(hx))>epsiftmp<eps1dk=-gx';endendx=x+alpha*dk;k=k+1;end这里调用的函数定义为function fori=1:nx(i)=sym(['x',num2str(i)]);endf=4*x(1)^2+x(2)^2-2*x(1)*x(2)+3*x(2)+7;g=jacobian(f,x);h=jacobian(g,x);调用这个函数的主程序为clear;clc;n=2;xx=3*ones(n,1);[x_min,f_min,kk]=Newtonfa('fun4',n,xx)§C.2.3共轭梯度法function[x_min,f_min,kk]=Gongetidufa(f,n,x0)%该函数是求解无约束最优化问题的共轭梯度算法eps=1.0e-6;fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;t0=0;h=0.1;x=x0;while1fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx=eval(gg);fx=eval(ff);ifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;endifk>0beta=gx*gx'/(gx_old*gx_oldFletcher-Reeves公式% beta=gx*(gx-gx_old)'/(gx_old*gx_oldPolak-Ribiere-Polyak公式% beta=gx*(gx-gx_old)'/((gx-gx_old)*dk_oldCrowder-Wolfe公式% beta=-gx*gx'/(gx_old*dk_oldDixon公式% beta=gx*gx'/((gx-gx_old)*dk_oldDai-Yuan公式dk=-gx'+beta*dk_old;elsedk=-gx';end% alpha=WolfePowell(ff,n,x,dk,t0,h);x=x+alpha*dk;gx_old=gx;dk_old=dk;k=k+1;end调用这个函数的主程序为clear;clc;n=2;xx=3*ones(n,1);[x_min,f_min,kk]=Gongetidufa('fun4',n,xx)§C.2.4拟牛顿法function[x_min,f_min,kk]=NiNewtonfa(f,n,x0)%该函数是求解无约束最优化问题的拟牛顿算法eps=1.0e-6;eps1=0.01;H=eye(n,n);fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;t0=0;h=0.1;x=x0;while1fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx=eval(gg);fx=eval(ff);ifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;endifk>0ss=x-x_old;bb=H*yy;aa=ss-bb;% A=aa*aa'/(yy'*aa);% H=H+A对称秩一校正B=ss*ss'/(ss'*yy)-bb*bb'/(bb'*yy);H=H+B;DFP（Davidon-Fletcher-Powell）校正公式% C=(aa*ss'+ss*aa')/(ss'*yy)-(aa'*yy)*ss*ss'/(ss'*yy)^2;% H=H+C;BFGS（Broyden-Flecther-Goldforb-Shanno）校正公式enddk=-H*gx';gx_old=gx;x_old=x;%alpha=Armijoonly(ff,n,x,dk);% alpha=WolfePowell(ff,n,x,dk,t0,h);x=x+alpha*dk;k=k+1;end调用这个函数的主程序为clear;clc;n=2;xx=3*ones(n,1);[x_min,f_min,kk]=NiNewtonfa('fun4',n,xx)§C.2.5信赖域算法§C.2.5.1信赖域方法子问题的折线步方法程序functionx_new=Zhexianbufa(A,b,h,x)%mins^tAs+b^ts,s.t.||s||<=h的折线步算法eps=1.0e-6;tmp1=b'*A*b;tmp2=norm(b);dc=-b/tmp2;iftmp1<0ac=h;elseac=min(h,tmp2^3/tmp1);endxc=x+ac*dc;x_new=xc;iftmp1>0xn=x-inv(A)*b;ifnorm(xn-x)<=hx_new=xn;returnelseifh>ac+epsx_tmp=xn-xc;nt=x_tmp'*x_tmp;nn=(xc-x)'*(xc-x);ll=(sqrt(delta)-nnt)/nt;x_new=(1-ll)*xc+ll*xn;endend调用这个函数的主程序为clear;clc;A=[3,1;1,2];b=ones(2,1);x=[2,2]';h=0.1;x_new=Zhexianbufa(A,b,h,x)§C.2.5.2信赖域方法（6.1）的程序function[x_min,f_min,kk]=Xinlaiyufangfa(f,n,x0)%该函数是求解无约束最优化问题的信赖域算法eps=1.0e-6;fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;x=x0;fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx=eval(gg);fx=eval(ff);hx=eval(hh);h=norm(gx);while1ifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;endx=Zhexianbufa(hx,gx',h,x);fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx1=eval(gg);fx1=eval(ff);hx1=eval(hh);xx=x-x0;nn=norm(xx);r=(fx-fx1)/(fx-q);ifr<0.25h=nn/4;ifr<=0x=x0;endelseifr>0.75ifabs(nn-h)<epsh=2*h;endendx0=x;fx=fx1;gx=gx1;hx=hx1;k=k+1;end调用这个函数的主程序为clear;clc;n=2;xx=3*ones(n,1);§C.2.5.3Levenberg-Marquardt方法（6.2）的程序function[x_min,f_min,kk]=xinlaiyuLMfa(f,n,x0)%Levenberg-Marquardt算法eps=1.0e-6;mu=0.01;fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;x=x0;fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx=eval(gg);fx=eval(ff);hx=eval(hh);while1ifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;endwhile1G=hx+mu*eye(n,n);R=chol(G);ifabs(det(R))>epsbreak;endendP=inv(R);gxx=gx*P;s=-P*gxx';x=x0+s;fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endgx1=eval(gg);fx1=eval(ff);hx1=eval(hh);nn=norm(s);q=s'*hx*s/2+gx*s+fx;r=(fx-fx1)/(fx-q);if r<0.25ifr<=0x=x0;endelseif r>0.75mu=mu/2;x0=x;fx=fx1;gx=gx1;hx=hx1;k=k+1;endend调用这个函数的主程序为clear;clc;n=2;xx=3*ones(n,1);[x_min,f_min,kk]=xinlaiyuLMfa('fun4',n,xx)§C.2.6最小二乘问题的算法§C.2.6.1线性最小二乘问题的正交分解法（算法7.1）的程序functionx_new=Llsmethod(A,b)%该函数是求解线性最小二乘问题的算法eps=1.0e-6;[m,n]=size(A);B=[Ab];[Q,R]=qr(B);RA=R;RA(:,n+1)=[];R(all(RA==0,2),:)=[];RA(all(RA==0,2),:)=[];Rb=R(:,n+1);x_new=inv(RA)*Rb;调用这个函数的主程序为clear;clc;A=[1,2,3,1;4,-2,2,4;3,-1,4,5;4,3,-4,-5;7,-6,5,4;2,-1,-3,6]b=[1,2,3,4,5,6]'x_new=Llsmethod(A,b)§C.2.6.2高斯牛顿法（算法7.2）的程序function[x_min,f_min,kk]=Zuixiaoerchengfa(f,n,x0)%该函数是求解非线性最小二乘问题的高斯牛顿算法eps=1.0e-6;fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;x=x0;while1fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endA=eval(JJ);b=-eval(gg);fx=eval(ff);gx=A'*bifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;enddk=Llsmethod(A,b);x=x+dk;k=k+1;end这里调用的函数定义为function n=2;fori=1:nx(i)=sym(['x',num2str(i)],'real');endT=50:5:120;R=[34.78,28.61,23.65,19.63,16.37,13.72,11.54,9.744,8.261,7.03,6.005,5.147,4.427,3.82,3.307];f=0;fori=1:15% res(i)=R(i)-x(1)*exp(x(2)/(T(i)+x(3)));res(i)=R(i)-x(1)+(T(i)+x(2))^2;f=f+res(i)^2;endres=res';J=jacobian(res);调用这个函数的主程序为clear;clc;n=2;xx=[1,1]';[x_min,f_min,kk]=Zuixiaoerchengfa('fun7',n,xx)§C.2.6.3求解Levenberg-Marquardt（7.3）程序function[x_min,f_min,kk]=LMmethodforLSP(f,n,x0)%Levenberg-Marquardt算法eps=1.0e-6;eta1=0.1;eta2=0.75;gamma1=0.5;gamma2=2;rr=0.5;mu=0.1;h_max=200;fori=1:nx(i)=sym(['x',num2str(i)]);endk=0;x=x0;h=1;fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endhx=eval(JJ)'*eval(JJ);gx=eval(JJ)'*eval(gg);fx=eval(ff);while1ifnorm(gx)<epsx_min=x;f_min=fx;kk=k;break;endA=hx+mu*eye(n,n);x=Zhexianbufa(A,gx,h,x);s=x-x0;fori=1:neval(['x'num2str(i)'=x('num2str(i)');'])endhx1=eval(JJ)'*eval(JJ);fx1=eval(ff);q=s'*hx*s/2+gx'*s+fx;iffx-q>epsr=(fx-fx1)/(fx-q);elser=1;endifr<eta1mu=mu*rrx=x0;elseifr>eta2x0=x;fx=fx1;gx=gx1;hx=hx1;k=k+1;endend调用这个函数的主程序为clear;clc;xx=[1,1]';[x_min,f_min,kk]=LMmethodforLSP('fun7',2,xx)§C.3不用导数的无约束最优化算例§C.3.1单纯形调优法程序function[x_min,f_min,kk]=DCtiaoyoufa(f,x0,l)%该函数是不用导数的最优化方法之单纯形调优算法x0初始点，x0须为列向量eps=1.0e-6;alpha=1;beta=0.5;gamma=2;delta=0.5;n=length(x0);A=x0*ones(1,n)+l*eye(n);A(:,n+1)=x0;fori=1:n+1fx(i)=feval(f,A(:,i));endkk=0;while1x_bar=zeros(n,1);fori=1:nx_bar=x_bar+A(:,z(i));endx_bar=(1/n)*x_barx_bar为非最大值点的重心x_r=x_bar+alpha*(x_bar-A(:,z(n+1x_r为反射点fxr=feval(f,x_r);iffxr<y(1)x_e=x_bar+gamma*(x_bar-A(:,z(n+1x_e为延伸点iffeval(f,x_e)<y(1)A(:,z(n+1))=x_e;fx(z(n+1))=feval(f,x_e);elseA(:,z(n+1))=x_r;fx(z(n+1))=fxr;endelseiffxr<y(n)A(:,z(n+1))=x_r;fx(z(n+1))=fxr;elseiffxr<y(n+1)A(:,z(n+1))=x_r;fx(z(n+1))=fxr;endx_c=x_bar-beta*(x_bar-A(:,z(n+1x_r为收缩点iffeval(f,x_c)<fx(z(n+1))A(:,z(n+1))=x_c;fx(z(n+1))=feval(f,x_c);elsefori=2:n+1A(:,z(i))=A(:,z(1))+delta*(A(:,z(i))-A(:,z(1)));%减小棱长fx(z(i))=feval(f,A(:,z(i)));endendendenderr1=0;fori=1:n+1err1=err1+(y(i)-f_bar)^2;enderr1=sqrt(err1/(n+1));iferr1<epsbreakendfxx=y(z(1))；xx=A(:,z(1))；kk=kk+1;endx_min=A(:,z(1));f_min=y(z(1));这个函数中的输入函数为function ff=fun8(x)G=[2,1;1,2];r=-3*ones(2,1);调用它的主程序为clear;clc;n=2;a0=-2*ones(n,1);h0=1;[x_min,f_min,kk]=DCtiaoyoufa(@fun8,a0,h0)§C.3.2模式搜索法程序function[x_min,f_min,kk]=MSsousuofa(f,x0,l)%该函数是不用导数的最优化方法之模式搜索算法eps=1.0e-6;beta=0.5;[n,m]=size(x0);A=eye(n);kk=0;x=x0;while1y=x;fori=1:niff(y+l(i)*A(:,i))<f(y)y=y+l(i)*A(:,i);elseiff(y-l(i)*A(:,i))<f(y)y=y-l(i)*A(:,i);elsel(i)=beta*l(i);endend %探测移动iff(y)<f(x)tmp=x;x=y;y=2*x-tmp; %模式移动elseify~=xy=x;elsea=max(l)ifa<epsbreak;elsel=beta*lendendkk=kk+1;endx_min=x;f_min=f(x);调用这个函数的主程序为clear;clc;n=2;x=ones(n,1);a0=5-10*rand(n,1)t=2;h0=t*ones(n,1);[x_min,f_min,kk]=MSsousuofa(f,a0,h0)§C.3.3Rosenbork法程序function[x_min,f_min,kk]=Rosenbork(f,x0,l)%Rosenbork算法eps=1.0e-6;beta=0.5;gamma=2;[n,m]=size(x0);A=eye(n);kk=0;x=x0;flag=0;while1y=x;while1z=y;fori=1:niff(y+l(i)*A(:,i))<f(y)y=y+l(i)*A(:,i);l(i)=gamma*l(i);elsel(i)=-beta*l(i);endendiff(y)<f(z)continuef(y)<f(x)x_old=x;x=y;ifnorm(x-x_old)<epsflag=1;endbreakelsea=max(abs(l));ifa<epsflag=1;break;endendendifflag==1breakendforj=1:nifabs(v(j))<epsB(:,j)=A(:,j);elseB(:,j)=zeros(n,1);fori=j:nB(:,j)=B(:,j)+v(i)*A(:,i);endendendA(:,1)=B(:,1)/norm(B(:,1));forj=2:nA(:,j)=B(:,j);fori=1:j-1A(:,j)=A(:,j)-B(:,j)'*A(:,i)*A(:,i);endendkk=kk+1;endx_min=x;f_min=f(x);调用这个函数的主程序为clear;clc;n=2;x=ones(n,1);a0=5-10*rand(n,1)t=2;h0=t*ones(n,1);[x_min,f_min,kk]=Rosenbork(f,a0,h0)§C.3.4Powell法程序示例一function[x_min,f_min,kk]=Powell_I(f,nn,x0,h0)%该函数是不用导数的最优化方法之Powell原始算法eps=1.0e-6;[n,m]=size(x0);A=eye(n);kk=0;x=x0;while1y=x;fori=1:nalpha0=0;y=y+alpha_min*A(:,i);endalpha0=0;y=y+alpha_min*d;x_old=x;x=y;ifnorm(x-x_old)<epsbreakelsefori=1:n-1A(:,i)=A(:,i+1);endA(:,n)=d;endkk=kk+1;endx_min=x;f_min=eval(ff);这里被调用的函数定义为function ff=fun8(x,n)G=[2,1;1,2];r=-3*ones(2,1);调用这个函数的主程序为clear;clc;n=2;xx=2*ones(n,1);h=0.1;[x_min,f_min,kk]=Powell_I('fun8',n,xx,h)示例二function[x_min,f_min,kk]=Powell_II(f,nn,x0,h0)%该函数是不用导数的最优化方法之Powell修正算法eps=1.0e-6;[n,m]=size(x0);A=eye(n);kk=0;x=x0;fx=[f,'(x,nn)'];y=x;xx=x;while1fori=1:nalpha0=0;y=y+alpha_min*A(:,i);x=y;fx=[f,'(x,nn)'];endfori=1:nz(i)=ff(i)-ff(i+1);end[zz,ii]=max(z);d=y-xx;alpha0=0;alpha_min=Hjfengefa(f,nn,y,d,a,b);y=y+alpha_min*d;alpha_bar=1+alpha_min;x=y;fxx=eval(fx);ifnorm(x-xx)<epsbreakelseifabs(alpha_bar)>tmpfori=ii:n-1AA(:,i)=AA(:,i+1);endAA(:,n)=d;endxx=x;ff