优化方法MATLAB编程 大连理工大学.docx
《优化方法MATLAB编程 大连理工大学.docx》由会员分享,可在线阅读,更多相关《优化方法MATLAB编程 大连理工大学.docx(19页珍藏版)》请在冰豆网上搜索。
优化方法MATLAB编程大连理工大学
优化方法上机大作业
学院:
姓名:
学号:
指导老师:
肖现涛
第一题
源程序如下:
functionzy_x=di1ti(x)
%di1ti是用来求解优化作业第一题的函数。
x0=x;yimuxulong=0.000001;
g0=g(x0);s0=-g0;
A=2*ones(100,100);
k=0;
whilek<100
lanmed=-(g0)'*s0/(s0'*A*s0);
x=x0+lanmed*s0;
g=g(x);
k=k+1;
ifnorm(g)zy_x=x;
fprintf('After%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolutionis\n%f.\n\nTheoptimal"x"is"ans".',k,f(x))
break;
end
miu=norm(g)^2/norm(g0)^2;
s=-g+miu*s0;
g0=g;s0=s;x0=x;
end
functionf=f(x)
f=(x'*ones(100,1))^2-x'*ones(100,1);
functiong=g(x)
g=(2*x'*ones(100,1))*ones(100,1)-ones(100,1);
代入x0,运行结果如下:
>>x=zeros(100,1);
>>di1ti(x)
After1iterations,obtaintheoptimalsolution.
Theoptimalsolutionis
-0.250000.
Theoptimal"x"is"ans".
ans=0.005*ones(100,1).
第二题
1.最速下降法。
源程序如下:
functionzy_x=di2titidu(x)
%该函数用来解大作业第二题。
x0=x;yimuxulong=1e-5;k=0;
g0=g(x0);s0=-g0;
whilek>=0
ifnorm(g0)break;
else
lanmed=10;c=0.1;i=0;
whilei>=0&i<100
x=x0+lanmed*s0;
iff(x)>(f(x0)+c*lanmed*g0'*s0)
lanmed=lanmed/2;
i=i+1;
else
break;
end
end
x=x0+lanmed*s0;
x0=x;
g0=g(x);
s0=-g0;
k=k+1;
end
end
zy_x=x;
zyj=f(x);
fprintf('after%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolutionis%f.\n\nTheoptimal"x"is"ans".\n',k,zyj);
functionf=f(x)
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2;
functiong=g(x)
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x4-x3^2)];
>>x=[-1.21-1.21]';
>>di2titidu(x)
after5945iterations,obtaintheoptimalsolution.
Theoptimalsolutionis0.000000.
Theoptimal"x"is"ans".
ans=
1.0000
1.0000
1.0000
1.0000
2.牛顿法
源程序如下:
functionzy_x=di2tinewton(x)
%该函数用来解大作业第二题。
x0=x;yimuxulong=1e-5;k=0;
g0=g(x0);h0=h(x0);s0=-inv(h0)*g0;
whilek>=0&k<1000
ifnorm(g0)break;
else
x=x0+s0;
x0=x;
g0=g(x);
h0=h(x);
s0=-inv(h0)*g0;
k=k+1;
end
end
zy_x=x;
zyj=f(x);
fprintf('after%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolutionis%f.\n\nTheoptimal"x"is"ans".\n',k,zyj);
functionf=f(x)
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2;
functiong=g(x)
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x4-x3^2)];
functionh=h(x)
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
h=[2+1200*x1^2-400*x2-400*x100;-400*x120000;002+1200*x3^2-400*x4-400*x3;00-400*x3200];
代入初始值,运行结果如下:
>>x=[-1.21-1.21]';
>>di2tinewton(x)
after6iterations,obtaintheoptimalsolution.
Theoptimalsolutionis0.000000.
Theoptimal"x"is"ans".
ans=
1.0000
1.0000
1.0000
1.0000
可以看出,用Newton法经过6次迭代就能求出最优解。
3.BFGS法
源程序如下:
functionzy_x=di2tiBFGS(x)
%该函数用来解大作业第二题。
x0=x;yimuxulong=1e-5;k=0;
g0=g(x0);H0=eye(4);s0=-H0*g0;
whilek>=0&k<100
ifnorm(g0)break;
else
lanmed=10;c=0.1;i=0;
whilei>=0&i<100
x=x0+lanmed*s0;
iff(x)>(f(x0)+c*lanmed*g0'*s0)
lanmed=lanmed/2;
i=i+1;
else
break;
end
end
x=x0+lanmed*s0;
dete_x=x-x0;
dete_g=g(x)-g0;
miu=1+dete_g'*H0*dete_g/(dete_x'*dete_g);
H=H0+(miu*dete_x*dete_x'-H0*dete_g*dete_x'-dete_x*dete_g'*H0)/(dete_x'*dete_g);
s=-H*g(x);
x0=x;
s0=s;
H0=H;
g0=g(x);
k=k+1;
end
end
zy_x=x;
zyj=f(x);
fprintf('after%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolutionis%f.\n\nTheoptimal"x"is"ans".\n',k,zyj);
functionf=f(x)
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
f=(x1-1)^2+(x3-1)^2+100*(x2-x1^2)^2+100*(x4-x3^2)^2;
functiong=g(x)
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
g=[2*(x1-1)-400*x1*(x2-x1^2);200*(x2-x1^2);2*(x3-1)-400*x3*(x4-x3^2);200*(x4-x3^2)];
代入初始值,计算结果如下:
>>x=[-1.21-1.21]';
>>di2tiBFGS(x)
after53iterations,obtaintheoptimalsolution.
Theoptimalsolutionis0.000000.
Theoptimal"x"is"ans".
ans=
1.0000
1.0000
1.0000
1.0000
第三题
1.惩罚函数法
源程序如下:
functionzy_x=di3ti(x)
%该函数用来解大作业第三题。
x0=x;M=100;c=4;m=1;
whilem>0
g0=g(x0,M);yimuxulong=1e-5;k=0;s0=-inv(H(x0,M))*g0;
whilek>=0
ifnorm(g0)break;
else
x=x0+s0;%牛顿法;
x0=x;
g0=g(x,M);
s0=-inv(H(x0,M))*g0;
k=k+1;
end
end
ifmax([abs(h(x)),g1(x),g2(x),g3(x)])<0.5
break;
else
M=M*c;
m=m+1;
end
end
zy_x=x;
zyj=f(x);
fprintf('after%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolutionis%f.\n\nTheoptimal"x"is"ans".\n',m,zyj);
functionF=F(x,M)
x1=[10]*x;
x2=[01]*x;
F=4*x1-x2^2-12+M*(h^2+g1^2+g2^2+g3^2);
functiong=g(x,M)
x1=[10]*x;
x2=[01]*x;
g=[4+M*(-4*(25-x1^2-x2^2)*x1+2*(10*x1-x1^2+10*x2-x2^2-34)*(10-2*x1)+2*x1);-2*x2+M*(-4*(25-x1^2-x2^2)*x2+2*(10*x1-x1^2+10*x2-x2^2-34)*(10-2*x2)+2*x2)];
functionH=H(x,M)
x1=[10]*x;
x2=[01]*x;
H=[M*(24*x1^2-120*x1+8*x2^2-40*x2+238),M*(16*x1*x2-40*x1-40*x2+200);M*(16*x1*x2-40*x1-40*x2+200),-2+M*(24*x2^2-120*x2+8*x1^2-40*x1+238)];
functionf=f(x)
x1=[10]*x;
x2=[01]*x;
f=4*x1-x2^2-12;
functionh=h(x)
x1=[10]*x;
x2=[01]*x;
h=25-x1^2-x2^2;
functiong1=g1(x)
x1=[10]*x;
x2=[01]*x;
g=10*x1-x1^2+10*x2-x2^2-34;
ifg<0
g1=g;
else
g1=0;
end
functiong2=g2(x)
x1=[10]*x;
x2=[01]*x;
ifx1>=0
g2=0;
else
g2=x1;
end
functiong3=g3(x)
x1=[10]*x;
x2=[01]*x;
ifx2>=0
g3=0;
else
g3=x2;
end
代入任意初始值,运算结果如下。
>>x=rand(2,1);
>>di3ticf(x)
after1iterations,obtaintheoptimalsolution.
Theoptimalsolutionis-31.490552.
Theoptimal"x"is"ans".
ans=
1.0024
4.8477
2.乘子法
源程序如下:
function[x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)
%功能:
用乘子法解一般约束问题:
minf(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;
fori=1:
l
btak=btak+he(i)^2;
end
%更新乘子向量
fori=1:
m
temp=min(gi(i),lambda(i)/c);
btak=btak+temp^2;
end
btak=sqrt(btak);
ifbtak>epsilon
ifk>=2&&btak>theta*btaold
c=eta*c;
end
fori=1:
l
mu(i)=mu(i)-c*he(i);
end
fori=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;
%增广拉格朗日函数
functionpsi=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;
fori=1:
l
psi=psi-he(i)*mu(i);
s1=s1+he(i)^2;
end
psi=psi+0.5*c*s1;
s2=0;
fori=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
functiongi=g1(x)
gi=10*x
(1)-x
(1)^2+10*x
(2)-x
(2)^2-34;
%目标函数的梯度文件df1.m
functiong=df1(x)
g=[4,-2*x
(2)]';
%等式约束(向量)函数的Jacobi矩阵(转置)文件dh1.m
functiondhe=dh1(x)
dhe=[-2*x
(1),-2*x
(2)]'
%不等式约束(向量)函数的Jacobi矩阵(转置)文件dg1.m
functiondgi=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(kgk=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(newfmk=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=[22]';
[x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0)
x=
1.0013
4.8987
mu=
0.7701
lambda=
0
0
0.9434
output=
fval:
-31.9923
iter:
4