2023年优化方法上机大作业.doc

优化措施 上机大作业

机械工程与材料能源学部
能源与动力学院

能源与环境工程

联络方式：

x0=[0;1]T;%初始值
s0=[-1;1]T;%初始搜索方向
c1=0.1;c2=0.5;a=0;b=inf;d=1;n=0;
x1=x0+d*s0;
g0=[(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));(x0(2)-x0(1)^2)];
g1=[(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));(x1(2)-x1(1)^2)];
f1=(x1(2)-x1(1)^2)^2+(1-x1(1))^2;
f0=(x0(2)-x0(1)^2)^2+(1-x0(1))^2;
while((f0-f1<-c1*d*g0'*s0)||(g1'*s0<c2*g0'*s0))
if ((f0-f1)<(-c1*d*g0'*s0))
b=d;d=(d+a)/2;
x1=x0+d*s0;
g0=[(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));(x0(2)-x0(1)^2)];
g1=[(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));(x1(2)-x1(1)^2)];
f1=(x1(2)-x1(1)^2)^2+(1-x1(1))^2;
f0=(x0(2)-x0(1)^2)^2+(1-x0(1))^2;
elseif (((g1')*s0)<(c2*(g0')*s0))
a=d;
if(2*d<=(d+b)/2)
d=2*d;
else
d=(d+b)/2;
end
x1=x0+d*s0;
g0=[(x0(2)-x0(1)^2)*x0(1)-2*(1-x0(1));(x0(2)-x0(1)^2)];
g1=[(x1(2)-x1(1)^2)*x1(1)-2*(1-x1(1));(x1(2)-x1(1)^2)];
f1=(x1(2)-x1(1)^2)^2+(1-x1(1))^2;
f0=(x0(2)-x0(1)^2)^2+(1-x0(1))^2;
end
end
x1

f1=(x1(2)-x1(1)^2)^2+(1-x1(1))^2

x1 =

-0.0000
1.0000

d =

1.1102e-016

f1 =

2

function f = fun( x )
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
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 = fun( x )
%UNTITLED4 Summary of this function goes here
% Detailed explanation goes here
g=[2 -2 0 0;-2 4 -1 0;0 -1 2 0;0 0 0 2]*x+[2;3;-1;0];

end
x0=[0;0;0;0]; %初始值
eps=1.0e-4; %精度
g0=gfun(x0);
s0=-g0;
n=0;
syms d1;
while norm(g0)>eps
if n<3
g=gfun(x0+d1*s0);
d= double(solve(s0'*g));
x1=x0+d*s0;
g1=gfun(x1);
if norm(g1)<eps
n=n+1;
x0=x1;
break
else
s0=-g1+(norm(g1)^2/norm(g0)^2)*s0;
x0=x1;
g0=g1;
end
elseif n=3
x0=x1;
g0=gfun(x0);
s0=-g0;
n=0;
end
n=n+1;
end
x0
n
fun(x0)

x0 =

-4
-3
-1
0

n =

3

ans =

-8

function f= fun3_1(x )
%FUN3 Summary of this function goes here
% Detailed explanation goes here
f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);
end

function g= gfun3_1(x)
%GFUN3_1 Summary of this function goes here
% Detailed explanation goes here
g=[1+2*x(1)*exp(x(1)^2+x(2)^2);4*x(2)+2*x(2)*exp(x(1)^2+x(2)^2)];
end
(1) 最速下降法
x0=[0;1];%初始值
eps=1.0e-5;%精度
n=0;
g0=gfun3_1(x0);
syms d1;
while norm(g0)>=eps
s0=-g0;
g=gfun3_1(x0+d1*s0);
d= double(solve(s0'*g));
x1=x0+d*s0;
g1=gfun3_1(x1);
if( norm(g1)<eps)
n=n+1;
x0=x1;
break;
else
x0=x1;
g0=gfun3_1(x0);
end
n=n+1;
end
f0=fun3_1(x0)
x0
n

f0 =

0.7729

x0 =

-0.4194
0.0000

n =

6

（2）牛顿法
function g2 = hesse(x)
%HESSE Summary of this function goes here
% Detailed explanation goes here
g2=[2*exp(x(1)^2+x(2)^2)+4*x(1)^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+2*exp(x(1)^2+x(2)^2)+4*x(2)^2*exp(x(1)^2+x(2)^2)];

end
x0=[0;1];%初始值
eps=1.0e-5;%精度
n=0;
g0=gfun3_1(x0);
g20=hesse(x0);
while norm(g0)>=eps
d=-g20\g0;
x1=x0+d;
g1=gfun3_1(x1);
if( norm(g1)<eps)
n=n+1;
x0=x1;
break;
else
x0=x1;
g0=gfun3_1(x0);
end
n=n+1;
end
f0=fun3_1(x0)
x0
n

f0 =

0.7729

x0 =

-0.4194
0.0000

n =

35

（3）BFGS

x0=[0;1];%初始值
eps=1.0e-5;%精度
n=0;
g0=gfun3_1(x0);
syms d1;
h0=eye(2);
while norm(g0)>=eps
s0=-h0*g0;
g=gfun3_1(x0+d1*s0);
d= double(solve(s0'*g));
x1=x0+d*s0;
g1=gfun3_1(x1);
if( norm(g1)<eps)
n=n+1;
x0=x1;
break;
else
h0=h0-(h0*(g1-g0)*(g1-g0)'*h0)/((g1-g0)'*h0*(g1-g0))+...
((x1-x0)*(x1-x0)')/((x1-x0)'*(g1-g0))+...
((g1-g0)'*h0*(g1-g0))*((x1-x0)*(x1-x0)');
x0=x1;
g0=gfun3_1(x0);
end
n=n+1;
end
f0=fun3_1(x0)
x0
n

f0 =

0.7729

x0 =

-0.4194
0.0000

n =

4

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)
Aee=[ ];
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)
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 = zeros(m,1);
end
callqpact.m文献
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)
运行成果：
callqpact
x =
1.4000
1.7000
lambda =
0.8000
exitflag =
0
output =

fval: -6.4500
iter: 7

Function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)
function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma)
function he=h1(x)
he=-x(1)^2-x(2)^2+25.0;
function f=f1(x)
f=4*x(1)-x(2)^2-12;
function gi=g1(x)
gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;
function dhe = dh1(x)
dhe = [-1*x(1), -1*x(2)]';
function dgi = dg1(x)
dgi = [10-2*x(1), 10-2*x(2)]';
function g=df1(x)
g = [4, -2.0*x(2)]';
function [x,val,k]=bfgs(fun,gfun,x0,varargin)
maxk=500;
sigma=2.0;
eta=2.0; theta=0.8;
k=0; ink=0;
epsilon=1e-5;
x=x0; he=feval(hf,x); gi=feval(gf,x);
n=length(x); l=length(he); m=length(gi);

mu=0.1*ones(l,1); lambda=0.1*ones(m,1);
btak=10; btaold=10;
while(btak>epsilon & k<maxk)
[x,ival,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma);
ink=ink+ik;
he=feval(hf,x); gi=feval(gf,x);
btak=0.0;
for (i=1:l), btak=btak+he(i)^2; end
for i=1:m
temp=min(gi(i),lambda(i)/sigma);
btak=btak+temp^2;
end
btak=sqrt(btak);
if btak>epsilon
if(k>=2 & btak > theta*btaold)
sigma=eta*sigma;
end
for (i=1:l), mu(i)=mu(i)-sigma*he(i); end
for (i=1:m)
lambda(i)=max(0.0,lambda(i)-sigma*gi(i));
end
end
k=k+1;
btaold=btak;
x0=x;
end
f=feval(fun,x);
output.fval=f;
output.iter=k;
output.inner_iter=ink;
output.bta=btak;
f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x);
l=length(he); m=length(gi);
psi=f; s1=0.0;
for(i=1:l)
psi=psi-he(i)*mu(i);
s1=s1+he(i)^2;
end
psi=psi+0.5*sigma*s1;
s2=0.0;
for(i=1:m)
s3=max(0.0, lambda(i) - sigma*gi(i));
s2=s2+s3^2-lambda(i)^2;
end
psi=psi+s2/(2.0*sigma);
maxk=500;
rho=0.55; sigma1=0.4; epsilon1=1e-5;
k=0; n=length(x0);
Bk=eye(n);
while(k<maxk)
gk=feval(gfun,x0,varargin{:});
if(norm(gk)<epsilon1), 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<oldf+sigma1*rho^m*gk'*dk)
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{:});

function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma)
dpsi=feval(dfun,x);
he=feval(hf,x); gi=feval(gf,x);
dhe=feval(dhf,x); dgi=feval(dgf,x);
l=length(he); m=length(gi);
for(i=1:l)
dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);
end
for(i=1:m)
dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);
end

>> x0=[1,1]';
[x,mu,lambda,output]=multphr('f1','h1','g1','df1','dh1','dg1',x0)

x =
1.0013
4.8987
mu =
2.0312
lambda =

0.7545
output =

fval: -31.9923
iter: 5
inner_iter: 58
bta: 4.3187e-07

function [x,mu,lam,val,k]=sqpm(x0,mu0,lam0)
maxk=100;
n=length(x0); l=length(mu0); m=length(lam0);
rho=0.5; eta=0.1; B0=eye(n);
x=x0; mu=mu0; lam=lam0;
Bk=B0; sigma=0.8;
epsilon1=1e-6; epsilon2=1e-5;
[hk,gk]=cons(x); dfk=df1(x);
[Ae,Ai]=dcons(x); Ak=[Ae; Ai];
k=0;
while(k<maxk)
[dk,mu,lam]=qpsubp(dfk,Bk,Ae,hk,Ai,gk);
mp1=norm(hk,1)+norm(max(-gk,0),1);
if(norm(dk,1)<epsilon1)&(mp1<epsilon2)
break;
end
deta=0.05;
tau=max(norm(mu,inf),norm(lam,inf));
if(sigma*(tau+deta)<1)
sigma=sigma;
else
sigma=1.0/(tau+2*deta);
end
im=0;
while(im<=20)
if(phi1(x+rho^im*dk,sigma)-phi1(x,sigma)<eta*rho^im*dphi1(x,sigma,dk))
mk=im;
break;
end
im=im+1;
if(im==20), mk=10; end
end
alpha=rho^mk; x1=x+alpha*dk;
[hk,gk]=cons(x1); dfk=df1(x1);
[Ae,Ai]=dcons(x1); Ak=[Ae; Ai];
lamu=pinv(Ak)'*dfk;
if(l>0&m>0)
mu=lamu(1:l); lam=lamu(l+1:l+m);
end
if(l==0), mu=[]; lam=lamu; end
if(m==0), mu=lamu; lam=[]; end
sk=alpha*dk;
yk=dlax(x1,mu,lam)-dlax(x,mu,lam);
if(sk'*yk>0.2*sk'*Bk*sk)
theta=1;
else
theta=0.8*sk'*Bk*sk/(sk'*Bk*sk-sk'*yk);
end
zk=theta*yk+(1-theta)*Bk*sk;
Bk=Bk+zk*zk'/(sk'*zk)-(Bk*sk)*(Bk*sk)'/(sk'*Bk*sk);
x=x1;
k=k+1;
end
val=f1(x);
function p=phi1(x,sigma)
f=f1(x); [h,g]=cons(x); gn=max(-g,0);
l0=length(h); m0=length(g);
if(l0==0), p=f+1.0/sigma*norm(gn,1); end
if(m0==0), p=f+1.0/sigma*norm(h,1); end
if(l0>0&m0>0)
p=f+1.0/sigma*(norm(h,1)+norm(gn,1));
end
function dp=dphi1(x,sigma,d)
df=df1(x); [h,g]=cons(x); gn=max(-g,0);
l0=length(h); m0=length(g);
if(l0==0), dp=df'*d-1.0/sigma*norm(gn,1); end
if(m0==0), dp=df'*d-1.0/sigma*norm(h,1); end
if(l0>0&m0>0)
dp=df'*d-1.0/sigma*(norm(h,1)+norm(gn,1));
end
function l=la(x,mu,lam)
f=f1(x);
[h,g]=cons(x);
l0=lemgth(h); m0=length(g);
if(l0==0), l=f-lam*g; end
if(m0==0), l=f-mu'*h; end
if(l0>0&m0>0)
l=f-mu'*h-lam'*g;
end
function dl=dlax(x,mu,lam)
df=df1(x);
[Ae,Ai]=dcons(x);
[m1,m2]=size(Ai); [l1,l2]=size(Ae);
if(l1==0), dl=df-Ai'*lam; end
if(m1==0), dl=df-Ae'*mu; end
if(l1>0&m1>0), dl=df-Ae'*mu-Ai'*lam; end
function f=f1(x)
f=4*x(1)-x(2)^2-12;
function df=df1(x)
df=[4, -2*x(2)]';
function [h,g]=cons(x)
h=[25-x(1)^2-x(2)^2];
g=[x(1);x(2)];
function [dh,dg]=dcons(x)
dh=[-2*x(1), -2*x(2)];
dg=[1 0; 0 1];
Qpsubp.m文献
function [d,mu,lam,val,k]=qpsubp(dfk,Bk,Ae,hk,Ai,gk)
n=length(dfk); l=length(hk); m=length(gk);
gamma=0.05; epsilon=1.0e-6; rho=0.5; sigma=0.2;
ep0=0.05; mu0=0.05*zeros(l,1); lam0=0.05*zeros(m,1);
d0=ones(n,1); u0=[ep0;zeros(n+l+m,1)];
z0=[ep0; d0; mu0;lam0,];
k=0;
z=z0; ep=ep0; d=d0; mu=mu0; lam=lam0;
while (k<=150)
dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);
if(norm(dh)<epsilon)
break;
end
A=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);
b=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)*u0-dh;
dz=A\b;
if(l>0&m>0)
de=dz(1); dd=dz(2:n+1); du=dz(n+2:n+l+1); dl=dz(n+l+2:n+l+m+1);
end
if(l==0)
de=dz(1); dd=dz(2:n+1); dl=dz(n+2:n+m+1);
end
if(m==0)
de=dz(1); dd=dz(2:n+1); du=dz(n+2:n+l+1);
end
i=0;
while (i<=20)
if(l>0&m>0) dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);
end
if(l==0) dh1=dah(ep+rho^i*de,d+rho^i*dd,mu,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);
end
if(m==0) dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam,dfk,Bk,Ae,hk,Ai,gk);
end
if(norm(dh1)<=(1-sigma*(1-gamma*ep0)*rho^i)*norm(dh))
mk=i;
break;
end
i=i+1;
if(i==20), mk=10; end
end
alpha=rho^mk;
if(l>0&m>0)
ep=ep+alpha*de; d=d+alpha*dd;
mu=mu+alpha*du; lam=lam+alpha*dl;
end
if(l==0)
ep=ep+alpha*de; d=d+alpha*dd;
lam=lam+alpha*dl;
end
if(m==0)
ep=ep+alpha*de; d=d+alpha*dd;
mu=mu+alpha*du;
end
k=k+1;
end
val=0.5*d'*Bk*d+dfk'*d;
function p=phi(ep,a,b)
p=a+b-sqrt(a^2+b^2+2*ep^2);
function dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk)
n=length(dfk); l=length(hk); m=length(gk);
dh=zeros(n+l+m+1,1);
dh(1)=ep;
if(l>0&m>0)
dh(2:n+1)=Bk*d-Ae'*mu-Ai'*lam+dfk;
dh(n+2:n+l+1)=hk+Ae*d;
for(i=1:m)
dh(n+l+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);
end
end
if(l==0)
dh(2:n+1)=Bk*d-Ai'*lam+dfk;
for(i=1:m)
dh(n+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);
end
end
if(m==0)
dh(2:n+1)=Bk*d-Ae'*mu+dfk;
dh(n+2:n+l+1)=hk+Ae*d;
end
dh=dh(:);
function bet=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)
dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);
bet=gamma*norm(dh)*min(1,norm(dh));
function [dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk)
m=length(gk);
dd1=zeros(m,m); dd2=zeros(m,m); v1=zeros(m,1);
for(i=1:m)
fm=sqrt(lam(i)^2+(gk(i)+Ai(i,:)*d)^2+2*ep^2);
dd1(i,i)=1-lam(i)/fm;
dd2(i,i)=1-(gk(i)+Ai(i,:)*d)/fm;
v1(i)=-2*ep/fm;
end
function A=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk)
n=length(dfk); l=length(hk); m=length(gk);
A=zeros(n+l+m+1,n+l+m+1);
[dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk);
if(l>0&m>0)
A=[1, zeros(1,n), zeros(1,l), zeros(1,m);
zeros(n,1), Bk, -Ae', - Ai';
zeros(l,1), Ae, zeros(l