最优化方法源程序1.doc

最优化方法源程序一、线性规划问题例 解线性规划函数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)执行结果：Optimization terminated.x = 200.0000 133.3333fval = -4.6667e+004exitflag = 1output = iterations: 4 algorithm: large-scale: interior point cgiterations: 0 message: Optimization terminated.二、无约束一元非线性规划问题例1 0.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)|(hdelta) if(phipphiq) 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; endds=abs(b-a); dphi=abs(phib-phia);if(phip s,phis,k,G,E=golds(inline(2*x2-x-1),-1,1,0.16,3)s = 0.236067977499790phis = -1.124611797498107k = 7G = -1.000000000000000 -0.236067977499790 0.236067977499790 1.000000000000000 -0.236067977499790 0.236067977499790 0.527864045000421 1.000000000000000 -0.236067977499790 0.055728090000841 0.236067977499790 0.527864045000421 0.055728090000841 0.236067977499790 0.347524157501472 0.527864045000421 0.055728090000841 0.167184270002524 0.236067977499790 0.347524157501472 0.167184270002524 0.236067977499790 0.278640450004206 0.347524157501472 0.167184270002524 0.209756742506940 0.236067977499790 0.278640450004206E = 0.111456180001683 0.012076339517143例2 最速降法function x,val,k=grad(fun,gfun,x0)% 功能: 用最速下降法求解无约束问题: min f(x)%输入: x0是初始点, fun, gfun分别是目标函数和梯度%输出: x, val分别是近似最优点和最优值, k是迭代次数.maxk=5000; %最大迭代次数rho=0.5;sigma=0.4;k=0; epsilon=1e-5;while(kmaxk) g=feval(gfun,x0); %计算梯度 d=-g; %计算搜索方向 if(norm(d)epsilon), break; end m=0; mk=0; while(m20) %Armijo搜索 if(feval(fun,x0+rhom*d) x0=1 1; x,val,k=grad(fun,gfun,x0)x = 0 0val = 0k = 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.9270ymin = -0.0279xmax = 0.7854ymax = -0.6448 (0.6448)P59 例3.5 f=exp(-x)+x2;fplot(f,0,1); %作图语句xmin,ymin=fminbnd (f,0,1)运行结果:xmin = 0.3517ymin =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 function f=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.0000fval = 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(Minimization of the Banana function)例 DFP算法求解无约束问题:香蕉函数function x,val,k=dfp(fun,gfun,x0)%功能: 用DFP算法求解无约束问题: min f(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(kmaxk) gk=feval(gfun,x0); %计算梯度 if(norm(gk)epsilon), break; end %检验终止准则 dk=-Hk*gk; %解方程组, 计算搜索方向 m=0; mk=0; while(m20) % 用Armijo搜索求步长 if(feval(fun,x0+rhom*dk)0) Hk=Hk-(Hk*yk*yk*Hk)/(yk*Hk*yk)+(sk*sk)/(sk*yk); end k=k+1; x0=x;endval=feval(fun,x0);function f=fun(x)f=100*(x(1)2-x(2)2+(x(1)-1)2;function gf=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.999999988038449val = 7.192197151384610e-017k = 29例 共轭梯度法求解无约束问题: 香蕉函数function x,val,k=frcg(fun,gfun,x0)% 功能: 用FR共轭梯度法求解无约束问题: min f(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=0.0) d=-g; end end if(norm(g)epsilon), break; end %检验终止条件 m=0; mk=0; while(m20) %Armijo搜索 if(feval(fun,x0+rhom*d) x0=-1,1 x,val,k= frcg (fun,gfun,x0)x = 0.999903558672420 0.999806719806067val = 9.317481519760974e-009k = 143例 信赖域法求解无约束问题: 香蕉函数function xk,val,k=trustm(x0)%功能: 牛顿型信赖域方法求解无约束优化问题 min f(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(k150) gk=gfun(x); if(norm(gk)epsilon) 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=eta2&norm(d)=dta) dta=min(tau2*dta,dtabar); else dta=dta; end end if(rketa1) 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;endxk=x;val=fun(xk);% 目标函数 %function f=fun(x) f=100*(x(1)2-x(2)2+(x(1)-1)2; % 子问题目标函数 %function qd=qk(x,d)gk=gfun(x); Bk=Hess(x);qd=gk*d+0.5*d*Bk*d;% 目标函数的梯度 %function gf=gfun(x)gf=400*x(1)*(x(1)2-x(2)+2*(x(1)-1), -200*(x(1)2-x(2);% 目标函数的Hesse阵 %function He=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)% 功能: 求解信赖域子问题: min qk(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) %Step1 of the algorithm dh=dah(mu,lam,d,gk,Bk,dta); if(norm(dh)epsilon) 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 (m20) dhnew=dah(mu+rhom*dmu,lam+rhom*dlam,d+rhom*dd,gk,Bk,dta); if(norm(dhnew)x0=-1,1 x,val,k= trustm(x0)x = 1.0000 1.0000val = 4.3868e-017k = 32例 应用MATLAB函数解无约束非线性规划% method1，采用BFGS法。 clf; clear x0=-1.9,2; %赋初值。% 这里是学习的重点： OPTIONS是控制fminunc和fminsearch指令的重要参数，%用optimset(属性，属性值,)指令改变设置，可以容易地控制算法。OPTIONS=optimset(LargeScale,off); %fminunc默认的大规模算法是“信赖域方法”，这是一种有效的算法；%将LargeScale的属性设置为off时,fminunc的默认中等规模的算法就是BFGS方法。 OPTIONS = optimset(OPTIONS,gradobj,on); %使用解析梯度。OPTIONS=optimset(OPTIONS,HessUpdate,bfgs); %此句可以省略%定义梯度函数GRAD=inline(100*(4*x(1)3-4*x(1)*x(2)+2*x(1)-2;100*(2*x(2)-2*x(1)2); f=inline(100*(x(2)-x(1)2)2+(1-x(1)2); %定义目标函数。x,fval,exitflag,output = fminunc(f,GRAD,x0,OPTIONS) % method2，采用DFP法。 OPTIONS=optimset(LargeScale,off); OPTIONS = optimset(OPTIONS,gradobj,on); OPTIONS=optimset(OPTIONS,HessUpdate,dfp); % 将HessUpdate属性设置为dfp就使fminunc指令采用DFP法。 GRAD=inline(100*(4*x(1)3-4*x(1)*x(2)+2*x(1)-2;100*(2*x(2)-2*x(1)2); f=inline(100*(x(2)-x(1)2)2+(1-x(1)2); x,fval,exitflag,output = fminunc(f,GRAD,x0,OPTIONS) % method3，采用最速下降法。 clf; clear x0=-1.9,2; %赋初值。 OPTIONS=optimset(LargeScale,off); OPTIONS = optimset(OPTIONS,gradobj,on); OPTIONS=optimset(OPTIONS,HessUpdate,steepdesc); %将HessUpdate属性设置为steepdesc就使fminunc指令采用最速下降法。 GRAD=inline(100*(4*x(1)3-4*x(1)*x(2)+2*x(1)-2;100*(2*x(2)-2*x(1)2); f=inline(100*(x(2)-x(1)2)2+(1-x(1)2); x,fval,exitflag,output = fminunc(f,GRAD,x0,OPTIONS); % method4，采用单纯形方法(Simplex Search)。 clf; clear x0=-1.9,2; %赋初值。 OPTIONS=optimset(LargeScale,off); OPTIONS = optimset(OPTIONS,gradobj,off); %该方法不使用导数，所以要设置gradobj属性为off。 f=inline(100*(x(2)-x(1)2)2+(1-x(1)2); % 套用本程序的格式定义目标函数。x,fval,exitflag,output = fminsearch(f,x0,OPTIONS)四、有约束非线性规划问题例 解二次规划的函数quadprog ()格式为:x,fval,exitflag,output = quadprog(H,c,A,b,Aeq,beq,lb,ub,x0,options)对应的二次规划矩阵形式： 例题4.11 解二次规划 即 H=2,0;0,2; c=6,0; A=-2,-1;-1,0;0,-1; b=-4,0,0;x,fval,exitflag,output = quadprog(H,c,A,b)执行结果：Optimization terminated.x = 1.0000 2.0000fval = 11.0000exitflag = 1output = iterations: 2 algorithm: medium-scale: active-set firstorderopt: cgiterations: message: Optimization terminated.或 OPTIONS=optimset(LargeScale,off);OPTIONS = optimset(OPTIONS,gradobj,off); x,fval,exitflag,output = quadprog(2,0;0,2, 6,0, -2,-1;-1,0;0,-1, -4,0,0, OPTIONS)或%二次规划求解clf; clear OPT=optimset;OPT.LargeScale=off;%不使用大规模问题算法H=2,0;0,2;c=6,0;A=-2,-1;-1,0;0,-1;b=-4,0,0;Aeq=;Beq=;x,fval,exitflag,output = quadprog(H,c,A,b,Aeq,Beq,OPT)例 函数fmincon ()有约束非线性规划问题为: 非线性规划求解的函数是fmincon,命令的基本格式如下：(1) x=fmincon(fun,X0,A,B)(2) x=fmincon(fun,X0,A,B,Aeq,Beq)(3) x=fmincon(fun,X0,A,B,Aeq,Beq,LB,UB)(4)x=fmincon(fun,X0,A,B,Aeq,Beq,LB,UB,nonlcon)(5)x=fmincon(fun,X0,A,b,Aeq,Beq,LB,UB,nonlcon,options) 其中x输出极值点，fun为M文件，X0迭代初值，LB，UB变量上下限，nonlcon 为给非线性约束写的M函数,options为参数说明。(6) x,fval= fmincon(.)注意：fmincon函数可能会给出局部最优解，这与初值X0的选取有关。例： 写成标准形式：s. t. 1）建立M文件nlfun.m定义目标函数：function f=nlfun(x);f=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);2）再建立M文件nlcon.m定义非线性约束：%非线性规划的非线性约束function g,ceq=nlcon(x)g=1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;ceq=x(1)2+x(2)-1;%等式约束3）建立主程序nlp_1.m：%非线性规划求解x0=1;-1;A=;B=;Aeq=;Beq=;%Aeq=;Beq=;LB=;UB=;OPT=optimset;OPT.LargeScale=off;%不使用大规模问题算法x,fval,exitflag,output=fmincon(nlfun,x0,A,B,Aeq,Beq,LB,UB,nlcon,OPT) 4）运行主程序结果:当x0=1;-1x = 1.358389899160045 -0.845222358215794fval = 13.718534108797590exitflag = 5output = iterations: 5 funcCount: 21 lssteplength: 0.500000000000000 stepsize: 0.001839539420487 algorithm: medium-scale: SQP, Quasi-Newton, line-search firstorderopt: 0.003406133987369 message: 1x172 char例4.12 求解约束非线性规划 运行主程序结果:x = -9.547405033629550 1.047405033629039fval = 0.023550379447862exitflag = 1output = iterations: 8 funcCount: 28 lssteplength: 1 stepsize: 4.250273253408122e-004 algorithm: medium-scale: SQP, Quasi-Newton, line-search firstorderopt: 8.512543350561175e-007 message: 1x144 char例4.13 求解约束非线性规划 主程序:%非线性规划求解x0=1;2;A=;B=;Aeq=1 1;Beq=3;LB=0 0;UB=2 2;OPT=optimset;OPT.LargeScale=off;%不使用大规模问题算法x,fval,exitflag,output=fmincon(nlfun_1,x0,A,B,Aeq,Beq,LB,UB,nlcon_1,OPT)目标函数:%非线性规划目标函数function f=nlfun_1(x)f=(4-x(2)*(x(1)-3)2;运行主程序结果:x = 2.000000000000000 1.000000000000000fval = 3exitflag = 1output = iterations: 2 funcCount: 9 lssteplength: 1 stepsize: 0 algorithm: medium-scale: SQP, Quasi-Newton, line-search firstorderopt: 2.220446049250313e-016 message: 1x144 charMat