大连理工优化方法大作业MATLAB编程

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)

b=d;

d=(d+a)/2;

end

if(p>=0)&&(q>=0)

dk=d;

x=x+d*s;

flag=1;

end

k=k+1;

if(p>=0)&&(q<0)

a=d;

d=min{2*d,(d+b)/2};

end

end

%定义求函数值的函数fun，当输入为x0=（x1，x2）时，输出为f function 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.0000

dk =1.1102e-016

k =54

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

end

function 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)];

end

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)

break;

end

m=0;mk=0;

while(m<20)

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

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,x0);

end

结果：

>> x0=[0,0,0,0];

>> [ x,val,k ] = frcg( 'fun','gfun',x0 )

x =

-4.0000 -3.0000 -1.0000 0

val =

-8.0000

k =

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(k<3)

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

end

end

f=fun(x);

end

function 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)

b=d;

d=(d+a)/2;