优化方法MATLAB编程大连理工大学.docx
《优化方法MATLAB编程大连理工大学.docx》由会员分享，可在线阅读，更多相关《优化方法MATLAB编程大连理工大学.docx（15页珍藏版）》请在冰豆网上搜索。
优化方法MATLAB编程大连理工大学
优化方法上机大作业
学院：
姓名：
学号：
指导老师：
肖现涛
第一题
q编写汁算机程序求解下列优化问XS,MATLABAC/C++均可,•将代码矜运行牆果打印'于考试前交到大黑1^81107室.
注;通过此习趣验证共辄梆度法求解门元二次曲款极小点至多需要卩次迭代+
源程序如下：
functionzy_x=dilti（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）vyimuxulong
zy_x=x;
fprintf（'After%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolutionis\n%f.\n\nTheoptimal"x"is"ans".',k,f（x））
break;
end
miu=norm（gF2/norm（gOF2;
s=-g+miu*s0;
gO=g;sO=s;xO=x;
end
functionf=f（x）
f=（x'*ones（100,1））A2-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）.
第二题
取初始点0=（-1.2,1,-1.2,l）rf＜虑下面无约束优化问題：
min（xi-I）2+（x3一I）2+100{x2-xf）2+100（x4一肩尸.
编写朮解上述问题的桂序，
«最速下降法，
O牛顿法T
oBFG5方法.
通过此习题体会上述三种算法的收蝕速度一
1.最速下降法。
源程序如下：
functionzy_x=di2titidu（x）
%该函数用来解大作业第二题x0=x;yimuxulong=1e-5;k=0;
g0=g（x0）;s0=-g0;
whilek>=0
ifnorm（gO）vyimuxulong
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\nTheoptimalsolution
is%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）A2+（x3-1）A2+100*（x2-x1A2）A2+100*（x4-x3A2）A2;
functiong=g（x）
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
g=[2*（x1-1）-400*x1*（x2-x1A2）;200*（x2-x1A2）;2*（x3-1）-400*x3*（x4-x3A2）;200*（x4-x3A
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）vyimuxulong
break;
else
x=x0+s0;
x0=x;
g0=g（x）;
h0=h（x）;
s0=-inv（h0）*g0;
k=k+1;
end
endzy_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）A2+（x3-1）A2+100*（x2-x1A2）A2+100*（x4-x3A2）A2;
functiong=g（x）
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
g=[2*（x1-1）-400*x1*（x2-x1A2）;200*（x2-x1A2）;2*（x3-1）-400*x3*（x4-x3A2）;200*（x4-x3A
2）]；functionh=h（x）
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
h=[2+1200*x1A2-400*x2-400*x100;-400*x120000;002+1200*x3A2-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;
gO=g（xO）;H0=eye（4）;s0=-H0*g0;
whilek>=0&k<100
ifnorm（g0）vyimuxulong
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=HO+（miu*dete_x*dete_x'-HO*dete_g*dete_x'-dete_x*dete_g'*HO）/（dete_x'*dete_g）;
s=-H*g（x）;
xO=x;
sO=s;
HO=H;
gO=g（x）;
k=k+1;
end
end
zy_x=x;
zyj=f（x）;
fprintf（'after%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolution
is%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）A2+（x3-1）A2+100*（x2-x1A2）A2+100*（x4-x3A2）A2;
functiong=g（x）
x1=[1000]*x;
x2=[0100]*x;
x3=[0010]*x;
x4=[0001]*x;
g=[2*（x1-1）-400*x1*（x2-x1A2）;200*（x2-x1A2）;2*（x3-1）-400*x3*（x4-x3A2）;200*（x4-x3A
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
第三题
习题3
编写求解下面的约束优化问题的惩罚函数法和乘子法程序:
min4xi—x孑一12
s.t.25—#—x孑=0
10xi—好+10x2—^2—34>0
Xi,X2>0
其中每决迭代需求解一个无约束极小化问題.
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）vyimuxulong
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\nTheoptimal
solutionis%f.\n\nTheoptimal"x"is"ans".\n',m,zyj）;
functionF=F（x,M）
x1=[10]*x;
x2=[01]*x;
F=4*x1-x2A2-12+M*（hA2+g1A2+g2A2+g3A2）;
functiong=g（x,M）
x1=[10]*x;
x2=[01]*x;
g=[4+M*（-4*（25-x1A2-x2A2）*x1+2*（10*x1-x1A2+10*x2-x2A2-34）*（10-2*x1）+2*x1）;-
2*x2+M*（-4*（25-x1A2-x2A2）*x2+2*（10*x1-xM2+10*x2-x2A2-34）*（10-2*x2）+2*x2）];
functionH=H（x,M）
x1=[10]*x;
x2=[01]*x;
H=[M*（24*x1A2-120*x1+8*x2A2-40*x2+238）,M*（16*x1*x2-40*x1-40*x2+200）;M*（16
*x1*x2-40*x1-40*x2+200）,-2+M*（24*x2A2-120*x2+8*x1A2-40*x1+238）];
functionf=f（x）
x1=[10]*x;
x2=[01]*x;
f=4*x1-x2A2-12;
functionh=h（x）
x1=[10]*x;
x2=[01]*x;
h=25-x1A2-x2A2;
functiong1=g1（x）
x1=[10]*x;
x2=[01]*x;
g=10*x1-x1A2+10*x2-x2A2-34;
ifg<0
g仁g;
else
g1=0；endfunctiong2=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&&kvmaxk）
%调用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）A2;
end
fori=1:
m
temp=min（gi（i）,lambda（i）/c）;btak=btak+tempA2;
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））;endk=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）A2;
end
psi=psi+0.5*c*s1;
s2=0;
fori=1:
m
s3=max（0,lambda（i）-c*gi（i））;s2=s2+s3A2-lambda（i）A2;
end
psi=psi+s2/（2*c）;
%不等式约束函数文件gl.m
functiongi=g1（x）gi=10*x
（1）-x
（1）A2+10*x
（2）-x
（2）八2-34;
%目标函数的梯度文件dfl.m
functiong=df1（x）
g=[4,-2*x
（2）]';
%等式约束（向量）函数的Jacobi矩阵（转置）文件dhl.m
functiondhe=dh1（x）
dhe=[-2*x
（1）,-2*x
（2）]'
%不等式约束（向量）函数的Jacobi矩阵（转置）文件dgl.mfunctiondgi=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（kvmaxk）
gk=feval（gfun,x0,varargin{:
}）;
if（norm（gk）break;
end
dk=-Bk\gk;
m=0;mk=0;
while（m<20）
newf=feval（fun,x0+rhoAm*dk,varargin{:
}）;
oldf=feval（fun,x0,varargin{:
}）;
if（newfmk=m;
break;
end
m=m+1;
end
x=x0+rhoAmk*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）;endk=k+1;x0=x;
endval=feval（fun,xO,varargin{:
}）;
结果
x=[22]';
[x,mu,lambda,output]=multphr（'fun','hf,'gf1','df,'dh','dg',xO）
x=
1.0013
4.8987
mu=
0.7701
lambda=
0
0
0.9434
output=
fval:
-31.9923
iter:
4
升级会员 