大连理工优化方法大作业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';

结果：

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;

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)

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;

end

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

dk=d;

flag=1;

end

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

a=d;

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

End

结果：x=[0,0,0,0];

>> [x,f,k]=second(x)

x =-4.0147 -3.0132 -1.0090 0

f = -7.9999

k = 1

function [f,x,k]=third_1(x)

k=0;

g=gfun(x);

while(norm(g)>=0.001)

s=-g;

dk=dfun(x,s);

x=x+dk*s;

k=k+1;

g=gfun(x);

f=fun(x);

end

function f=fun(x)

f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);

function gf=gfun(x)

gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function [j_1,j_2]=con(x,d,s)

j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';

j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';

function dk=dfun(x,s)%获取步长

flag=0;

a=0;

while(flag==0)

[p,q]=con(x,d,s);

if (p<0)

b=d;

d=(d+a)/2;

end

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

dk=d;

flag=1;

end

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

a=d;

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

end

end

结果：

x=[0,1];

[f,x,k]=third_1(x)

f =0.7729

x = -0.4196 0.0001

k =8

（1）程序：

function [f,x,k]=third_2(x)

k=0;

H=inv(ggfun(x));

g=gfun(x);

while(norm(g)>=0.001)

s=(-H*g')';

dk=dfun(x,s);

x=x+dk*s;

k=k+1;

g=gfun(x);

f=fun(x);

end

function f=fun(x)

f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);

function gf=gfun(x)

gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)

ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^ 2);

4*x(1)*x(2)*exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2)]

;

function [j_1,j_2]=con(x,d,s)

j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)'; j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s)% 步长获取

flag=0;

a=0;

d=1;

b=10000;

while(flag==0)

[p,q]=con(x,d,s);

if (p<0)

b=d;

d=(d+a)/2;

end

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

dk=d;

flag=1;

end

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

a=d;

if 2*d>=(d+b)/2

d=(d+b)/2;

else d=2*d;

end

end

End

结果：

x=[0,1];

[f,x,k]=third_2(x)

f =0.7729

x = -0.4193 0.0001

k =8

（2）程序：

function [f,x,k]=third_3(x)

k=0;

X=cell(2);

g=cell(2);

X{1}=x;

H=eye(2);

g{1}=gfun(X{1});

s=(-H*g{1}')';

dk=dfun(X{1},s);

X{2}=X{1}+dk*s;

g{2}=gfun(X{2});

while(norm(g{2})>=0.001)

dx=X{2}-X{1};

dg=g{2}-g{1};

v=dx/(dx*dg')-(H*dg')'/(dg*H*dg');

h1=H*dg'*dg*H/(dg*H*dg');

h2=dx'*dx/(dx*dx');

h3=dg*H*dg'*v'*v;

H=H-h1+h2+h3;

k=k+1;

X{1}=X{2};

g{1}=gfun(X{1});

s=(-H*g{1}')';

dk=dfun(X{1},s);

X{2}=X{1}+dk*s;

g{2}=gfun(X{2});

norm(g{2});

f=fun(x);

x=X{2};

end

function f=fun(x)

f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);

function gf=gfun(x)

gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)

ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^ 2);4*x(1)*x(2)*exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^ 2);

function [p,q]=con(x,d,s)

p=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';

q=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';

function dk=dfun(x,s)

flag=0;

a=0;

d=1;

b=10000;

while(flag==0)

[p,q]=con(x,d,s);

if (p<0)

b=d;

d=(d+a)/2;

end

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

dk=d;

end

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

a=d;

if 2*d>=(d+b)/2

d=(d+b)/2;

else d=2*d;

end

end

end

结果：

x=[0,1];

[f,x,k]=third_3(x)

f =0.7729

x = -0.4195 0.0000

k=6

function callqpact

H=[2 0; 0 2];

c=[-2 -5]';

Ae=[ ]; be=[ ];

Ai=[1 -2; -1 -2; -1 2;1 0;0 1];

bi=[-2 -6 -2 0 0]';

x0=[0 0]';

[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)

function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilon=1.0e-9; err=1.0e-6;

k=0; x=x0; n=length(x); kmax=1.0e3;

ne=length(be); ni=length(bi); lamk=zeros(ne+ni,1);

index=ones(ni,1);

for (i=1:ni)

if(Ai(i,:)*x>bi(i)+epsilon), index(i)=0; end

end

while(k<=kmax)

if(ne>0), Aee=Ae; end

for(j=1:ni)

if(index(j)>0), Aee=[Aee; Ai(j,:)]; end end

gk=H*x+c;

[m1,n1] = size(Aee);

[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1)); if(norm(dk)<=err)

y=0.0;

if(length(lamk)>ne)

[y,jk]=min(lamk(ne+1:length(lamk))); end

if(y>=0)

exitflag=0;

else

exitflag=1;

for(i=1:ni)

if(index(i)&(ne+sum(index(1:i)))==jk) index(i)=0; break;

end

end

end

k=k+1;

else

exitflag=1;

alpha=1.0; tm=1.0;

for(i=1:ni)

if((index(i)==0)&(Ai(i,:)*dk<0))

tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);

if(tm1

tm=tm1; ti=i;

end

end

end

alpha=min(alpha,tm);

x=x+alpha*dk;

if(tm<1), index(ti)=1; end

end

if(exitflag==0), break; end

k=k+1;

end

output.fval=0.5*x'*H*x+c'*x;

output.iter=k;

function [x,lambda]=qsubp(H,c,Ae,be)

ginvH=pinv(H);

[m,n]=size(Ae);

if(m>0)

rb=Ae*ginvH*c + be;

lambda=pinv(Ae*ginvH*Ae')*rb;

x=ginvH*(Ae'*lambda-c);

else

x=-ginvH*c;

lambda=0;

end

结果

>>callqpact

x =

1.4000

1.7000

lambda =

0.8000

exitflag =

output =

fval: -6.4500

iter: 7

function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)

%功能: 用乘子法解一般约束问题: min f(x), s.t. h(x)=0, g(x).=0

%输入: x0是初始点, fun, dfun分别是目标函数及其梯度；

% hf, dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；

% gf, dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；

%输出: x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量; % output是结构变量, 输出近似极小值f, 迭代次数, 内迭代次数等

maxk=500;

c=2.0;

eta=2.0;theta=0.8;

k=0;ink=0;

epsilon=0.00001;

x=x0;he=feval(hf,x);gi=feval(gf,x);

n=length(x);l=length(he);m=length(gi);

mu=zeros(l,1);lambda=zeros(m,1);

btak=10;btaold=10;

while(btak>epsilon&&k

%调用BFGS算法程序求解无约束子问题

[x,ival,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c); ink=ink+ik;

he=feval(hf,x);gi=feval(gf,x);

btak=0;

for i=1:l

btak=btak+he(i)^2;

end

%更新乘子向量

for i=1:m

temp=min(gi(i),lambda(i)/c);

btak=btak+temp^2;

end

btak=sqrt(btak);

if btak>epsilon

if k>=2&&btak>theta*btaold

c=eta*c;

end

for i=1:l

mu(i)=mu(i)-c*he(i);

end

for i=1:m

lambda(i)=max(0,lambda(i)-c*gi(i));

end

k=k+1;

btaold=btak;

x0=x;

end

end

f=feval(fun,x);

output.fval=f;

output.iter=k;

%增广拉格朗日函数

function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c)

f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);

l=length(he);m=length(gi);

psi=f;s1=0;

for i=1:l

psi=psi-he(i)*mu(i);

s1=s1+he(i)^2;

end

psi=psi+0.5*c*s1;

s2=0;

for i=1:m

s3=max(0,lambda(i)-c*gi(i));

s2=s2+s3^2-lambda(i)^2;

end

psi=psi+s2/(2*c);

%不等式约束函数文件g1.m

function gi=g1(x)

gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;

%目标函数的梯度文件df1.m

function g=df1(x)

g=[4, -2*x(2)]';

%等式约束（向量）函数的Jacobi矩阵（转置）文件dh1.m function dhe=dh1(x)

dhe=[-2*x(1), -2*x(2)]'

%不等式约束（向量）函数的Jacobi矩阵（转置）文件dg1.m

function dgi=dg1(x)

dgi=[10-2*x(1), 10-2*x(2)]';

function [x,val,k]=bfgs(fun,gfun,x0,varargin)

maxk=500;

rho=0.55;sigma=0.4;epsilon=0.00001;

k=0;n=length(x0);

Bk=eye(n);

while(k

gk=feval(gfun,x0,varargin{:});

if(norm(gk)

break;

end

dk=-Bk\gk;

m=0;mk=0;

while(m<20)

newf=feval(fun,x0+rho^m*dk,varargin{:});

oldf=feval(fun,x0,varargin{:});

if(newf

mk=m;

break;

end

m=m+1;

end

x=x0+rho^mk*dk;

sk=x-x0;

yk=feval(gfun,x,varargin{:})-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,varargin{:});

结果

x=[2 2]';

[x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0) x =

1.0013

4.8987

mu =

0.7701

lambda =

0.9434

output =

fval: -31.9923

iter: 4

f=[3,1,1];

A=[2,1,1;1,-1,-1];

b=[2;-1];

lb=[0,0,0];

x=linprog(f,A,b,zeros(3),[0,0,0]',lb) 结果：

Optimization terminated.

x =

0.0000

0.5000

0.5000