最优化方法源程序1.docx
- 文档编号:29830510
- 上传时间:2023-07-27
- 格式:DOCX
- 页数:31
- 大小:132.41KB
最优化方法源程序1.docx
《最优化方法源程序1.docx》由会员分享,可在线阅读,更多相关《最优化方法源程序1.docx(31页珍藏版)》请在冰豆网上搜索。
最优化方法源程序1
最优化方法源程序
一、线性规划问题
例解线性规划函数linprog()
线性规划规划矩阵形式:
格式为:
x=linprog(c,A,B,Aeq,Beq,LB,UB)
例解线性规划
即
%线性规划求解
clear;
c=[-100,-200]';
A=[1,1;1,0];
B=[500,200]';
Aeq=[2,6];
Beq=1200;
LB=[0,0]';UB=[];
[x,fval,exitflag,output]=linprog(c,A,B,Aeq,Beq,LB,UB)
执行结果:
Optimizationterminated.
x=
200.0000
133.3333
fval=
-4.6667e+004
exitflag=
1
output=
iterations:
4
algorithm:
'large-scale:
interiorpoint'
cgiterations:
0
message:
'Optimizationterminated.'
二、无约束一元非线性规划问题
例10.618法
function[s,phis,k,G,E]=golds(phi,a,b,delta,epsilon)
%功能:
0.618法精确线搜索
%输入:
phi是目标函数,a,b是搜索区间的两个端点
%delta,epsilon分别是自变量和函数值的容许误差
%输出:
s,phis分别是近似极小点和极小值,G是nx4矩阵,
%其第k行分别是a,p,q,b的第k次迭代值[ak,pk,qk,bk],
%E=[ds,dphi],分别是s和phis的误差限.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=(sqrt(5)-1)/2;h=b-a;phia=feval(phi,a);phib=feval(phi,b);
p=a+(1-t)*h;q=a+t*h;phip=feval(phi,p);phiq=feval(phi,q);
k=1;G(k,:
)=[a,p,q,b];
while(abs(phib-phia)>epsilon)|(h>delta)
if(phip b=q;phib=phiq;q=p;phiq=phip; h=b-a;p=a+(1-t)*h;phip=feval(phi,p); else a=p;phia=phip;p=q;phip=phiq; h=b-a;q=a+t*h;phiq=feval(phi,q); end k=k+1;G(k,: )=[a,p,q,b]; end ds=abs(b-a);dphi=abs(phib-phia); if(phip<=phiq) s=p;phis=phip; else s=q;phis=phiq; end E=[ds,dphi]; 程序调用: >>[s,phis,k,G,E]=golds(inline('2*x^2-x-1'),-1,1,0.16,3) s= 0.236067977499790 phis= -1.124611797498107 k= 7 G= -1.000000000000000-0.2360679774997900.2360679774997901.000000000000000 -0.2360679774997900.2360679774997900.5278640450004211.000000000000000 -0.2360679774997900.0557280900008410.2360679774997900.527864045000421 .0557******** .0557******** 0.1671842700025240.2360679774997900.2786404500042060.347524157501472 0.1671842700025240.2097567425069400.2360679774997900.278640450004206 E= 0.1114561800016830.012076339517143 例2最速降法 function[x,val,k]=grad(fun,gfun,x0) %功能: 用最速下降法求解无约束问题: minf(x) %输入: x0是初始点,fun,gfun分别是目标函数和梯度 %输出: x,val分别是近似最优点和最优值,k是迭代次数. maxk=5000;%最大迭代次数 rho=0.5;sigma=0.4; k=0;epsilon=1e-5; while(k g=feval(gfun,x0);%计算梯度 d=-g;%计算搜索方向 if(norm(d) m=0;mk=0; while(m<20)%Armijo搜索 if(feval(fun,x0+rho^m*d) mk=m;break; end m=m+1; end x0=x0+rho^mk*d; k=k+1; end x=x0; val=feval(fun,x0); functionf=fun(x) f=x (1)^2+4*x (2)^2; functiongf=gfun(x) gf=[2*x (1),8*x (2)]'; functionHe=Hess(x) n=length(x); He=zeros(n,n); He=[2,0; 0,8]; 程序调用: >>x0=[11]'; >>[x,val,k]=grad('fun','gfun',x0) x= 0 0 val= 0 k= 2 例MATLAB函数fminbnd() 用于一元函数无约束优化的局部最优解 常用格式如下: (1)x=fminbnd(fun,x1,x2) (2)x=fminbnd(fun,x1,x2,options) (3)[x,fval]=fminbnd(...) (4)[x,fval,exitflag]=fminbnd(...) (5)[x,fval,exitflag,output]=fminbnd(...) 其中(3)、(4)、(5)的等式右边可选用 (1)或 (2)的等式右边。 函数fminbnd的算法基于黄金分割法和二次插值法,它要求目标函数必须是连续函数,并可能只给出局部最优解. 1)fval是函数f(x)在解x处的值. 2)exitflag的值描述程序运行情况,如果exitflag的值大于0,则程序收敛于解;如果exitflag的值等于0,则函数计算达到最大次数;如果exitflag的值小于0,则程序未收敛到解. 3)output输出程序运行的某些信息.其中output.iterations是迭代次数;output.funccount是函数计算次数;output.algorithm程序所使用的算法;output.firstorderopt是一阶最优性的度量. 例求 在 中的最大值与最小值. f='2*exp(-x).*sin(x)'; fplot(f,[0,8]);%作图语句 [xmin,ymin]=fminbnd(f,0,8) f1='-2*exp(-x).*sin(x)';%标准函数为求最小故乘-1 [xmax,ymax]=fminbnd(f1,0,8) 运行结果: xmin= 3.9270 ymin= -0.0279 xmax= 0.7854 ymax= -0.6448(0.6448) P59例3.5 f='exp(-x)+x^2'; fplot(f,[0,1]);%作图语句 [xmin,ymin]=fminbnd(f,0,1) 运行结果: xmin= 0.3517 ymin= 0.8272 三、无约束多元非线性规划问题 例MATLAB函数函数fminunc() 用于多元函数无约束优化的局部最优解 命令格式为: (1)x=fminunc(fun,X0);或x=fminsearch(fun,X0) (2)x=fminunc(fun,X0,options); 或x=fminsearch(fun,X0,options) (3)[x,fval]=fminunc(...); 或[x,fval]=fminsearch(...) (4)[x,fval,exitflag]=fminunc(...); 或[x,fval,exitflag]=fminsearch (5)[x,fval,exitflag,output]=fminunc(...); 或[x,fval,exitflag,output]=fminsearch(...) fminsearch是用单纯形法寻优.fminunc的算法见以下几点说明: [1]fminunc为无约束优化提供了大型优化和中型优化算法。 由options中的参数LargeScale控制: LargeScale=’on’(默认值),使用大型算法 LargeScale=’off’(默认值),使用中型算法 [2]fminunc为中型优化算法的搜索方向提供了4种算法,由options中的参数HessUpdate控制: HessUpdate=’bfgs’(默认值),拟牛顿法的BFGS公式; HessUpdate=’dfp’,拟牛顿法的DFP公式; HessUpdate=’steepdesc’,最速下降法 [3]fminunc为中型优化算法的步长一维搜索提供了两种算法,由options中参数LineSearchType控制: LineSearchType=’quadcubic’(缺省值),混合的二次和三次多项式插值; LineSearchType=’cubicpoly’,三次多项式插 使用fminunc和fminsearch可能会得到局部最优解. P60例3.6 functionf=fun1(x) f=exp(x (1))*(4*x (1)^2+2*x (2)^2+4*x (1)*x (2)+2*x (2)+1); >>[x,fval]=fminunc('fun1',[-1,1]) x= 0.5000-1.0000 fval= 3.6609e-015 例香蕉函数 该问题有精确解 绘等高线图: [X,Y]=meshgrid(-2: .125: 2,-1: .125: 3); Z=100*(X.^2-Y).^2+(X-1).^2; conts=exp(3: 20); [C,h]=contour(X,Y,Z,conts); clabel(C,h)%等高线填标签 xlabel('x1'),ylabel('x2'),title('MinimizationoftheBananafunction') 例DFP算法求解无约束问题: 香蕉函数 function[x,val,k]=dfp(fun,gfun,x0) %功能: 用DFP算法求解无约束问题: minf(x) %输入: x0是初始点,fun,gfun分别是目标函数及其梯度 %输出: x,val分别是近似最优点和最优值,k是迭代次数. maxk=1e5;%给出最大迭代次数 rho=0.55;sigma=0.4;epsilon=1e-5; k=0;n=length(x0); Hk=eye(n); while(k gk=feval(gfun,x0);%计算梯度 if(norm(gk) dk=-Hk*gk;%解方程组,计算搜索方向 m=0;mk=0; while(m<20)%用Armijo搜索求步长 if(feval(fun,x0+rho^m*dk) mk=m;break; end m=m+1; end %DFP校正 x=x0+rho^mk*dk; sk=x-x0;yk=feval(gfun,x)-gk; if(sk'*yk>0) Hk=Hk-(Hk*yk*yk'*Hk)/(yk'*Hk*yk)+(sk*sk')/(sk'*yk); end k=k+1;x0=x; end val=feval(fun,x0); functionf=fun(x) f=100*(x (1)^2-x (2))^2+(x (1)-1)^2; functiongf=gfun(x) gf=[400*x (1)*(x (1)^2-x (2))+2*(x (1)-1),-200*(x (1)^2-x (2))]'; 程序调用: >>x0=[0,0]' >>[x,val,k]=dfp('fun','gfun',x0) x= 0.999999994334768 0.999999988038449 val= 7.192197151384610e-017 k= 29 例共轭梯度法求解无约束问题: 香蕉函数 function[x,val,k]=frcg(fun,gfun,x0) %功能: 用FR共轭梯度法求解无约束问题: minf(x) %输入: x0是初始点,fun,gfun分别是目标函数和梯度 %输出: x,val分别是近似最优点和最优值,k是迭代次数. maxk=5000;%最大迭代次数 rho=0.6;sigma=0.4; k=0;epsilon=1e-4; 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) m=0;mk=0; while(m<20)%Armijo搜索 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,x); functionf=fun(x) f=100*(x (1)^2-x (2))^2+(x (1)-1)^2; functiongf=gfun(x) gf=[400*x (1)*(x (1)^2-x (2))+2*(x (1)-1),-200*(x (1)^2-x (2))]'; 程序调用: >>x0=[-1,1]' >>[x,val,k]=frcg('fun','gfun',x0) x= 0.999903558672420 0.999806719806067 val= 9.317481519760974e-009 k= 143 例信赖域法求解无约束问题: 香蕉函数 function[xk,val,k]=trustm(x0) %功能: 牛顿型信赖域方法求解无约束优化问题minf(x) %输入: x0是初始迭代点 %输出: xk是近似极小点,val是近似极小值,k是迭代次数 n=length(x0);x=x0;dta=1; eta1=0.15;eta2=0.75;dtabar=2.0; tau1=0.5;tau2=2.0;epsilon=1e-6; k=0;Bk=Hess(x);%Bk=eye(n); while(k<150) gk=gfun(x); if(norm(gk) break; end [d,val,lam,ik]=trustq(gk,Bk,dta); deltaq=-qk(x,d); deltaf=fun(x)-fun(x+d); rk=deltaf/deltaq; if(rk<=eta1) dta=tau1*dta; elseif(rk>=eta2&norm(d)==dta) dta=min(tau2*dta,dtabar); else dta=dta; end end if(rk>eta1) x0=x;x=x+d; %sk=x-x0;yk=gfun(x)-gfun(x0); %vk=sqrt(yk'*Bk*yk)*(sk/(sk'*yk)-Bk*yk/(yk'*Bk*yk)); %Bk=Bk-Bk*yk*yk'*Bk/(yk'*Bk*yk)+sk*sk'/(sk'*yk)+vk*vk' %pause Bk=Hess(x); end k=k+1; end xk=x; val=fun(xk); %%%目标函数%%%%%%%%%%%%%%% functionf=fun(x) f=100*(x (1)^2-x (2))^2+(x (1)-1)^2; %%%子问题目标函数%%%%%%%%%%%%% functionqd=qk(x,d) gk=gfun(x);Bk=Hess(x); qd=gk'*d+0.5*d'*Bk*d; %%%目标函数的梯度%%%%%%%%%%%%%% functiongf=gfun(x) gf=[400*x (1)*(x (1)^2-x (2))+2*(x (1)-1),-200*(x (1)^2-x (2))]'; %%%目标函数的Hesse阵%%%%%%%%%%%%%% functionHe=Hess(x) n=length(x); He=zeros(n,n); He=[1200*x (1)^2-400*x (2)+2,-400*x (1); -400*x (1),200]; function[d,val,lam,k]=trustq(gk,Bk,dta) %功能: 求解信赖域子问题: minqk(d)=gk'*d+0.5*d'*Bk*d,s.t.||d||<=delta %输入: gk是xk处的梯度,Bk是第k次近似Hesse阵,dta是当前信赖域半径 %输出: d,val分别是子问题的最优点和最优值,lam是乘子值,k是迭代次数. n=length(gk);gamma=0.05; epsilon=1.0e-6;rho=0.6;sigma=0.2; mu0=0.05;lam0=0.05; d0=ones(n,1);z0=[mu0,lam0,d0']'; u0=[mu0,zeros(1,n+1)]'; k=0;%k为迭代次数 z=z0; mu=mu0;lam=lam0;d=d0; while(k<=150)%Step1ofthealgorithm dh=dah(mu,lam,d,gk,Bk,dta); if(norm(dh) break; end A=JacobiH(mu,lam,d,Bk,dta);b=beta(mu,lam,d,gk,Bk,dta,gamma)*u0-dh; B=inv(A);dz=B*b; dmu=dz (1);dlam=dz (2);dd=dz(3: n+2); m=0;mk=0; while(m<20) dhnew=dah(mu+rho^m*dmu,lam+rho^m*dlam,d+rho^m*dd,gk,Bk,dta); if(norm(dhnew)<=(1-sigma*(1-gamma*mu0)*rho^m)*dh) mk=m; break; end m=m+1; end alpha=rho^mk; mu=mu+alpha*dmu; lam=lam+alpha*dlam; d=d+alpha*dd; k=k+1; end val=gk'*d+0.5*d'*Bk*d; %%%%%%%%%%%%%%%%%%%%%%%%%% functionp=phi(mu,a,b) p=a+b-sqrt((a-b)^2+4*mu); %%%%%%%%%%%%%%%%%%%%%%%%%% functiondh=dah(mu,lam,d,gk,Bk,dta) n=length(d); dh (1)=mu;dh (2)=phi(mu,lam,dta^2-norm(d)^2); mh=(Bk+lam*eye(n))*d+gk; for(i=1: n) dh(2+i)=mh(i); end dh=dh(: ); %%%%%%%%%%%%%%%%%%%%%%%%%% functionbet=beta(mu,lam,d,gk,Bk,dta,gamma) dh=dah(mu,lam,d,gk,Bk,dta); bet=gamma*norm(dh)*min(1,norm(dh)); %%%%%%%%%%%%%%%%%%%%%%%%%% functionA=JacobiH(mu,lam,d,Bk,dta) n=length(d); A=zeros(n+2,n+2); pmu=-4*mu/sqrt((lam+norm(d)^2-dta^2)^2+4*mu^2); thetak=(lam+norm(d)^2-dta^2)/sqrt((lam+norm(d)^2-dta^2)^2+4*mu^2); A=[1,0,zeros(1,n); pmu,1-thetak,-2*(1+thetak)*d'; zeros(n,1),d,Bk+lam*eye(n)]; 程序调用: >>x0=[-1,1]' >>[x,val,k]=trustm(x0) x= 1.0000 1.0000 val= 4.3868e-017 k= 32 例应用MATLAB函数解无约束非线性规划 %method=1,采用BFGS法。 clf;clear x0=[-1.9,2];%赋初值。 %这里是学习的重点: OPTIONS是控制fminunc和fminsearch指令的重要参数, %用optimset('属性
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 优化 方法 源程序