最优化方法实验报告(2)Word格式.doc
《最优化方法实验报告(2)Word格式.doc》由会员分享，可在线阅读，更多相关《最优化方法实验报告(2)Word格式.doc（23页珍藏版）》请在冰点文库上搜索。
d）令，转b）。
（二）牛顿法
牛顿法是基于多元函数的泰勒展开而来的，它将作为搜索方向，因此它的迭代公式可直接写出来：
用牛顿法求无约束问题的算法步骤如下：
a）给定初始点,精度，并令k=0；
b）若，停止，极小点为，否则转c）；
c）计算；
（三）共轭梯度法
共轭梯度法是利用目标函数梯度逐步产生共轭方向作为线搜索方向的方法，每次搜索方向都是在目标函数梯度的共轭方向，搜索步长通过一维极值算法确定。
a）给定初始点,精度；
c）；
d）用一维搜索方法求，使得
令，转e）；
e）若，停止，极小值为，否则转f）；
f）若转c），否则转g）；
g）令，
三、实验内容：
1．最速下降法的MATLAB实现
2．牛顿法的MATLAB实现
3．共轭梯度法的MATLAB实现
四、实验过程：
1．最速下降法的函数：
function[x,minf]=minFD（f,x0,var,eps）
%最速下降法主函数
ifnargin==3
eps=1.0e-6;
end
symsl;
tol=1;
gradf=-jacobian（f,var）;
whiletol>
eps
v=Funval（gradf,var,x0）;
tol=norm（v）;
y=x0+l*v;
yf=Funval（f,var,y）;
[a,b]=minJT（yf,0,0.1）;
%进退法求区间
xm=minHJ（yf,a,b）;
%黄金分割法
x1=x0+xm*v;
x0=x1;
x=x1;
minf=Funval（f,var,x）;
%进退法函数
function[minx,maxx]=minJT（f,x0,h0,eps）
x1=x0;
k=0;
h=h0;
while1
x4=x1+h;
k=k+1;
f4=subs（f,findsym（f）,x4）;
f1=subs（f,findsym（f）,x1）;
iff4<
f1
x2=x1;
x1=x4;
f2=f1;
f1=f4;
h=2*h;
else
ifk==1
h=-h;
x2=x4;
f2=f4;
else
x3=x2;
x2=x1;
x1=x4;
break;
end
end
minx=min（x1,x3）;
maxx=x1+x3-minx;
%黄金分割法函数
function[x,minf]=minHJ（f,a,b,eps）
l=a+0.382*（b-a）;
u=a+0.618*（b-a）;
k=1;
tol=b-a;
eps&
&
k<
100000
fl=subs（f,findsym（f）,l）;
fu=subs（f,findsym（f）,u）;
iffl>
fu
a=l;
l=u;
u=a+0.618*（b-a）;
b=u;
u=l;
l=a+0.382*（b-a）;
tol=abs（b-a）;
ifk==100000
disp（'
找不到最小值！
'
）;
x=NaN;
minf=NaN;
return;
x=（a+b）/2;
minf=subs（f,findsym（f）,x）;
2．牛顿法的函数：
function[x,minf]=minNT（f,x0,var,eps）
x0=transpose（x0）;
gradf=jacobian（f,var）;
jacf=jacobian（gradf,var）;
v=Funval（gradf,var,x0）;
pv=Funval（jacf,var,x0）;
p=-inv（pv）*transpose（v）;
p=double（p）;
x1=x0+p;
3.共轭梯度法的函数：
function[x,minf]=minGETD（f,x0,var,eps）
n=length（var）;
v0=Funval（gradf,var,x0）;
p=-transpose（v0）;
iftol<
=eps
x=x0;
break;
end
y=x0+l*p;
x1=x0+xm*p;
vk=Funval（gradf,var,x1）;
tol=norm（vk）;
x=x1;
ifk+1==n
x0=x1;
continue;
lamda=dot（vk,vk）/dot（v,v）;
p=-transpose（vk）+lamda*p;
k=k+1;
五、实验结果（总结/方案）
1、最速下降法：
用最速下降法求函数的极小值，初始点取
。
在commandwindow中输入：
>
symsts;
f=（t-4）^2+（s+2）^2+1;
[t,mf]=minFD（f,[1-3],[ts]）
输出结果：
x=4.0000-2.0000
mf=1
2、牛顿法：
用牛顿法求函数的极小值，其中初始点取为。
[t,mf]=minNT（f,[00],[ts]）
x=4
-2
mf=1
3、共轭梯度法：
用共轭梯度法求函数的极小值，其中初始值取。
f=（t-3）^2+s^2;
[t,mf]=minGETD（f,[-2,6],[ts]）
x=3.0000
0.0000
mf=2.00116e-037
实验四
约束最优化方法的MATLAB实现
2013年05月10日星期三实验成绩：
一、实验目的：
通过本次实验使学生较为熟练使用MATLAB软件，并能利用该软件进行约束最优化方法的计算。
二、实验内容：
1．罚函数法的MATLAB实现
2．可行方向法的MATLAB实现
三、实验背景：
（一）罚函数
外罚点函数是通过一系列罚因子，求罚函数的极小点来逼近原约束问题的最优点。
之所以称为外点罚函数法，是因为他是从可行域外部向约束边界逐步靠拢的。
用外点罚函数法求解线性约束优化问题的算法过程如下：
1）给定初始点，罚参数列及精度，置k=1；
2）构造罚函数；
3）用某种无约束非线性规划，以为初始点求解；
4）设最优解为，若满足某种终止条件，则停止迭代输出，否则令，转2）。
罚函数列的选法：
通常先选定一个初始常数和一个比例系数，则其余的可表示为。
终止条件可采用，其中。
（二）可行方向法
可行方向法是求解如下约束最优化问题的算法，其中约束条件为线性约束，A为约束系数矩阵，b为约束向量。
其基本思路是从可行点出发，沿着目标函数值减小的方向搜索求出新的可行点，如此迭代下去。
可行点是满足约束条件的点。
2、算法步骤
可行方向法的算法过程如下：
1）给定初始可行点，使其满足约束条件，置；
2）在处，将A，b分解成和，使，；
3）如果是空的，则令p=1，否则令;
4）计算，若则转6），否则转5）；
5）若是空的，则停止计算，输出，否则计算；
如果，则停止计算，输出，若包含负的分量，则选择一个负分量，去掉对应的行，转3）；
6）求一维约束问题，其中的计算方法如下：
，求出最优解，设为，令，置，转2）。
1、罚函数法的函数
function[x,minf]=minPF（f,x0,A,b,c1,p,var,eps）
formatlong;
ifnargin==7
eps=1.0e-4;
FE=0;
fori=1:
length（b）
FE=FE+（var*transpose（A（1,:
））-b（i））^2;
x1=transpose（x0）;
x2=inf;
M=c1*p;
FF=M*FE;
SumF=f+FF;
[x2,minf]=minNT（SumF,transpose（x1）,var）;
%牛顿法函数
ifnorm（x2-x1）<
x=x2;
c1=M;
x1=x2;
formatshort;
2、可行方向法的函数
function[x,minf]=minRosen（f,A,b,x0,var,eps）
ifnargin==5
sz=size（A）;
m=sz
（1）;
gf=jacobian（f,var）;
bConti=1;
whilebConti
k=0;
s=0;
A1=A;
A2=A;
b1=b;
b2=b;
fori=1:
m
dfun=A（i,:
）*x0-b（i）;
ifabs（dfun）<
0.000000001
k=k+1;
A1（k,:
）=A（i,:
b1（k,1）=b（i）;
s=s+1;
A2（s,:
b2（s,1）=b（i）;
ifk>
0
A1=A1（1:
k,:
b1=b1（1:
ifs>
A2=A2（1:
s,:
b2=b2（1:
while1
P=eye（n,n）;
ifk>
tM=transpose（A1）;
P=P-tM*inv（A1*tM）*A1;
gv=Funval（gf,var,x0）;
gv=transpose（gv）;
d=-P*gv;
ifd==0
ifk==0
x=x0;
bConti=0;
break;
else
w=inv（A1*tM）*A1*gv;
ifw>
=0
x=x0;
bConti=0;
break;
else
[u,index]=min（w）;
sA1=size（A1）;
ifsA1
（1）==1
k=0;
else
k=sA1
（2）;
A1=[A1（1:
（index-1）,:
A1（（index+1）:
sA1
（2）,:
）];
end
end
end
yl=x0+l*d;
tmpf=Funval（f,var,yl）;
bb=b2-A2*x0;
dd=A2*d;
ifdd>
=0
[tmpI,lm]=minJT（tmpf,0,0.1）;
lm=inf;
fori=1:
length（dd）
ifdd（i）<
ifbb（i）/dd（i）<
lm
lm=bb（i）/dd（i）;
[xm,minf]=minHJ（tmpf,0,lm,1.0e-14）;
tol=norm（xm*d）;
iftol<
eps
x0=x0+xm*d;
end
升级会员 