优化方法MATLAB编程 大连理工大学Word格式.docx
《优化方法MATLAB编程 大连理工大学Word格式.docx》由会员分享，可在线阅读，更多相关《优化方法MATLAB编程 大连理工大学Word格式.docx（19页珍藏版）》请在冰点文库上搜索。
.'
k,f（x））
break;
end
miu=norm（g）^2/norm（g0）^2;
s=-g+miu*s0;
g0=g;
s0=s;
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"
.
ans=0.005*ones（100,1）.
第二题
1.最速下降法。
functionzy_x=di2titidu（x）
%该函数用来解大作业第二题。
yimuxulong=1e-5;
k=0;
s0=-g0;
whilek>
=0
ifnorm（g0）<
else
lanmed=10;
c=0.1;
i=0;
whilei>
=0&
i<
iff（x）>
（f（x0）+c*lanmed*g0'
*s0）
lanmed=lanmed/2;
i=i+1;
x0=x;
g0=g（x）;
zy_x=x;
zyj=f（x）;
fprintf（'
after%diterations,obtaintheoptimalsolution.\n\nTheoptimalsolutionis%f.\n\nTheoptimal"
.\n'
k,zyj）;
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;
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"
ans=
1.0000
1.0000
2.牛顿法
functionzy_x=di2tinewton（x）
h0=h（x0）;
s0=-inv（h0）*g0;
k<
1000
else
x=x0+s0;
h0=h（x）;
s0=-inv（h0）*g0;
functionh=h（x）
h=[2+1200*x1^2-400*x2-400*x100;
-400*x120000;
002+1200*x3^2-400*x4-400*x3;
00-400*x3200];
代入初始值，运行结果如下：
di2tinewton（x）
after6iterations,obtaintheoptimalsolution.
可以看出，用Newton法经过6次迭代就能求出最优解。
3.BFGS法
functionzy_x=di2tiBFGS（x）
H0=eye（4）;
s0=-H0*g0;
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'
s=-H*g（x）;
H0=H;
代入初始值，计算结果如下：
di2tiBFGS（x）
after53iterations,obtaintheoptimalsolution.
第三题
1.惩罚函数法
functionzy_x=di3ti（x）
%该函数用来解大作业第三题。
M=100;
c=4;
m=1;
whilem>
g0=g（x0,M）;
s0=-inv（H（x0,M））*g0;
whilek>
%牛顿法；
g0=g（x,M）;
s0=-inv（H（x0,M））*g0;
ifmax（[abs（h（x））,g1（x）,g2（x）,g3（x）]）<
0.5
M=M*c;
m=m+1;
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）
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）
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）];
f=4*x1-x2^2-12;
h=25-x1^2-x2^2;
functiong1=g1（x）
g=10*x1-x1^2+10*x2-x2^2-34;
ifg<
g1=g;
else
g1=0;
functiong2=g2（x）
ifx1>
g2=0;
g2=x1;
functiong3=g3（x）
ifx2>
g3=0;
g3=x2;
代入任意初始值，运算结果如下。
x=rand（2,1）;
di3ticf（x）
after1iterations,obtaintheoptimalsolution.
Theoptimalsolutionis-31.490552.
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;
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&
&
maxk）
%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs（'
mpsi'
'
dmpsi'
x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c）;
ink=ink+ik;
he=feval（hf,x）;
btak=0;
fori=1:
l
btak=btak+he（i）^2;
%更新乘子向量
m
temp=min（gi（i）,lambda（i）/c）;
btak=btak+temp^2;
btak=sqrt（btak）;
ifbtak>
epsilon
ifk>
=2&
btak>
theta*btaold
c=eta*c;
mu（i）=mu（i）-c*he（i）;
lambda（i）=max（0,lambda（i）-c*gi（i））;
btaold=btak;
f=feval（fun,x）;
output.fval=f;
output.iter=k;
%增广拉格朗日函数
functionpsi=mpsi（x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c）
psi=f;
s1=0;
fori=1:
psi=psi-he（i）*mu（i）;
s1=s1+he（i）^2;
psi=psi+0.5*c*s1;
s2=0;
s3=max（0,lambda（i）-c*gi（i））;
s2=s2+s3^2-lambda（i）^2;
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）
rho=0.55;
sigma=0.4;
n=length（x0）;
Bk=eye（n）;
while（k<
gk=feval（gfun,x0,varargin{:
}）;
if（norm（gk）<
epsilon）
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+sigma*rho^m*gk'
*dk）
mk=m;
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）;
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.9434
output=
fval:
-31.9923
iter:
4
升级会员 