Matlab最优化编程例子

题目：分别用最速下降法、FR共轭梯度法、DFP法和BFGS法求解问题：

22minf(x)?x1?2x1x2?4x2?x1?3x2

取初始点x(1)?(1,1)T，通过Matlab编程实现求解过程。 公用函数如下：

1、function f= fun( X ) %所求问题目标函数

f=X(1)^2-2*X(1)*X(2)+4*X(2)^2+X(1)-3*X(2); end

2、function g= gfun( X ) %所求问题目标函数梯度

g=[2*X(1)-2*X(2)+1,-2*X(1)+8*X(2)-3]; end

3、function He = Hess( X ) %所求问题目标函数Hesse矩阵 n=length(X); He=zeros(n,n); He=[2,-2; -2,4];

End

解法一：最速下降法

function [ x,val,k ] = grad( fun,gfun,x0 ) %功能：用最速下降法求无约束问题最小值

%输入：x0是初始点，fun和gfun分别是目标函数和梯度 %输出：x、val分别是最优点和最优值，k是迭代次数 maxk=5000;%最大迭代次数 rho=0.5;sigma=0.4; k=0;eps=10e-6; while(k

g=feval(gfun,x0);%计算梯度 d=-g;%计算搜索方向 if(norm(d)

if(feval(fun,x0+rho^m*d)

x0=x0+rho^mk*d;

k=k+1; end x=x0;

val=feval(fun,x0); end

解法二：FR共轭梯度法

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)

while(m<20)

if(feval(fun,x0+rho^m*d)

x0=x0+rho^mk*d; val=feval(fun,x0); g0=g;d0=d; k=k+1; end x=x0;

val=feval(fun,x0); end

解法三：DFP法

function [ x,val,k ] = dfp( fun,gfun,x0 ) %功能：用DFP法求无约束问题最小值

%输入：x0是初始点，fun和gfun分别是目标函数和梯度 %输出：x、val分别是最优点和最优值，k是迭代次数 maxk=5000;%最大迭代次数 rho=0.5;sigma=0.4; k=0;eps=10e-6; n=length(x0);

Hk=inv(feval('Hess',x0)); while(k

gk=feval(gfun,x0); if(norm(gk)

dk=-Hk*gk'; dk=dk'; m=0;mk=0; while(m<20)

if(feval(fun,x0+rho^m*dk)

?P校正

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); end

解法四：BFGS法

function [ x,val,k ] = bfgs( fun,gfun,x0 ) %功能：用BFGS法求无约束问题最小值

%输入：x0是初始点，fun和gfun分别是目标函数和梯度 %输出：x、val分别是最优点和最优值，k是迭代次数 maxk=5000;%最大迭代次数 rho=0.5;sigma=0.4; k=0;eps=10e-6; n=length(x0); Bk=eye(n); while(k

gk=feval(gfun,x0); if(norm(gk)

dk=-Bk*gk'; m=0;mk=0; while(m<20)

new=sigma*rho^m*gk*dk; old=feval(fun,x0);

if(feval(fun,x0+rho^m*dk')

?GS校正

x=x0+rho^mk*dk'; sk=x-x0;

yk=feval(gfun,x)-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); end