MATLB机械优化设计程序.docx
《MATLB机械优化设计程序.docx》由会员分享，可在线阅读，更多相关《MATLB机械优化设计程序.docx（40页珍藏版）》请在冰豆网上搜索。
MATLB机械优化设计程序
例2-1求
在点
和点
的梯度。
%例2-1梯度的计算
symsx1x2%定义符号变量
f=x1^2+x2^2-4*x1+4;%定义二维目标函数
gradf=jacobian（f）%计算函数梯度
Xzuobiao1=[3,2];Xzuobiao2=[2,0];%定义Xzuobiao点坐标
gfk1=subs（subs（gradf,Xzuobiao1
（1））,Xzuobiao1
（2））%计算Xzuobiao1点的梯度值
gmk1=norm（gfk1）%计算Xzuobiao1点的梯度模
igk1=gfk1/gmk1%计算Xzuobiao1点的梯度单位向量
gfk2=subs（subs（gradf,Xzuobiao2
（1））,Xzuobiao2
（2））%计算Xzuobiao1点的梯度值
gmk2=norm（gfk2）%计算Xzuobiao1点的梯度模
igk2=gfk2/gmk2%计算Xzuobiao1点的梯度单位向量
gradf=
[2*x1-4,2*x2]
gfk1=
24
gmk1=
4.4721
igk1=
0.44720.8944
gfk2=
00
gmk2=
0
Warning:
Dividebyzero.
igk2=
NaNNaN
例2-2把函数
在点
展开泰勒二次近似式
%例2-2Taylor展开
symsx1x2
f=4+4.5*x1-4*x2+x1^2+2*x2^2-2*x1*x2+x1^4-2*x1^2*x2
disp（'函数f的表达式:
'）
pretty（simplify（f））;
%计算函数的一阶偏导数
dx1=diff（f,x1）;
dx2=diff（f,x2）;
disp（'函数f的一阶偏导数表达式:
'）
pretty（simplify（dx1））;
pretty（simplify（dx2））;
%计算函数的二阶偏导数
dx1x1=diff（f,x1,2）;
dx1x2=diff（dx1,x2）;
dx2x1=diff（dx2,x1）;
dx2x2=diff（f,x2,2）;
%根据函数f的二阶偏导数，构成Hessian矩阵
disp（'函数f的二阶偏导数表达式：
'）
pretty（simplify（dx1））;
H=[dx1x1dx1x2;dx2x1dx2x2];
pretty（simplify（H））;
%计算xk点的值
x1=2.0;x2=2.5;
disp（'函数在xk点的函数值:
'）
fk=subs（f）
disp（'函数在xk点的一节偏导数矩阵:
'）
dk=subs（[dx1dx2]）
disp（'函数xk点的海色矩阵:
'）
HK=subs（[dx1x1dx1x2;dx2x1dx2x2]）
disp（'函数在xk点的二阶Taylor展开式:
'）
symsx1x2
fkTL=fk+dk*[x1-2.0;x2-2.5]+0.5*[x1-2.0,x2-2.5]*HK*[x1-2.0;x2-2.5];
pretty（simplify（fkTL））;
f=
4+9/2*x1-4*x2+x1^2+2*x2^2-2*x1*x2+x1^4-2*x1^2*x2
函数f的表达式:
2242
4+9/2x1-4x2+x1+2x2-2x1x2+x1-2x1x2
函数f的一阶偏导数表达式:
3
9/2+2x1-2x2+4x1-4x1x2
2
-4+4x2-2x1-2x1
函数f的二阶偏导数表达式：
3
9/2+2x1-2x2+4x1-4x1x2
[2]
[2+12x1-4x2-2-4x1]
[]
[-2-4x14]
函数在xk点的函数值:
fk=
5.5000
函数在xk点的一节偏导数矩阵:
dk=
15.5000-6.0000
函数xk点的海色矩阵:
HK=
40-10
-104
函数在xk点的二阶Taylor展开式:
22
32-79/2x1+4x2+20x1-10x1x2+2x2
例2-3求函数
的极值点和极值
%例2-3求函数的极值
symsx1x2x3
f=2*x1^2+5*x2^2+x3^2+2*x2*x3+2*x1*x3-6*x2+3;
disp（'函数f的表达式:
'）
pretty（simplify（f））;
latex（f）;
%计算函数的1阶偏导数
dsx1=diff（f,x1）;
dsx2=diff（f,x2）;
dsx3=diff（f,x3）;
disp（'函数f的1阶偏导数：
'）
pretty（simplify（dsx1））;
pretty（simplify（dsx2））;
pretty（simplify（dsx3））;
%计算函数的2阶偏导数
dsx1x1=diff（f,x1,2）;
dsx1x2=diff（dsx1,x2）;
dsx1x3=diff（dsx1,x3）;
dsx2x1=diff（dsx2,x1）;
dsx2x2=diff（f,x2,2）;
dsx2x3=diff（dsx2,x3）;
dsx3x1=diff（dsx3,x1）;
dsx3x2=diff（dsx3,x2）;
dsx3x3=diff（f,x3,2）;
%根据函数f的2阶偏导数，构成海色矩阵
disp（'函数f的2阶偏导数矩阵'）
H=[dsx1x1dsx1x2dsx1x3;dsx2x1dsx2x2dsx2x3;dsx3x1dsx3x2dsx3x3]
%计算海色矩阵的正定性
[D,p]=chol（subs（H））;
ifp==0;disp（'海色矩阵为正定，函数f有极小点：
'）;end
%计算极值存在的必要条件，求极值点坐标
[x1,x2,x3]=solve（dsx1,dsx2,dsx3,'x1,x2,x3'）;
disp（'极值点坐标:
'）
fprintf（1,'x1=%3.4f\n',subs（x1））;
fprintf（1,'x2=%3.4f\n',subs（x2））;
fprintf（1,'x3=%3.4f\n',subs（x3））;
disp（'在极值点，函数f数值：
'）
fmb=subs（f）
M文件的运行结果如下
函数f的表达式:
222
2x1+5x2+x3+2x2x3+2x1x3-6x2+3
函数f的1阶偏导数：
4x1+2x3
10x2+2x3-6
2x3+2x2+2x1
函数f的2阶偏导数矩阵
H=
[4,0,2]
[0,10,2]
[2,2,2]
海色矩阵为正定，函数f有极小点：
极值点坐标:
x1=1.0000
x2=1.0000
x3=-2.0000
在极值点，函数f数值：
fmb=0
例2-5已知二维约束问题
受约束为
例2-5MATLAB实现，用M文件判别函数的凸性：
%例2-5判别函数的凸性
symsx1x2
f=60-10*x1-4*x2+x1^2+x2^2-x1*x2;
disp（'函数f的表达式:
'）
pretty（simplify（f））;
dsx1=diff（f,x1）;
dsx2=diff（f,x2）;
disp（'函数f的1阶偏导数:
'）
pretty（simplify（dsx1））;
pretty（simplify（dsx2））;
%计算函数的2阶偏导数
dsx1x1=diff（f,x1,2）;
dsx1x2=diff（dsx1,x2）;
dsx2x1=diff（dsx2,x1）;
dsx2x2=diff（f,x2,2）;
%根据函数f的2阶偏导数，构成海色矩阵
disp（'函数f的2阶偏导数矩阵'）
H=[dsx1x1dsx1x2;dsx2x1dsx2x2]
%计算函数矩阵的正定性
[d,p]=chol（subs（H））;
ifp==0;disp（'海色矩阵为正定，函数f为凸函数'）;end
M文件的运行结果如下
函数f的表达式:
22
60-10x1-4x2+x1+x2-x1x2
函数f的1阶偏导数:
-10+2x1-x2
-4+2x2-x1
函数f的2阶偏导数矩阵
H=
[2,-1]
[-1,2]
海色矩阵为正定，函数f为凸函数
%例2-6K-T条件
symsx1x2v%定义目标函数和约束函数的符号变量
%目标函数和约束函数
f=（x1-3）^2+x2^2;%目标函数f的表达式
g1=x1^2+x2-4;
g2=-x2;
g3=-x1;
v=[x1,x2];
%计算xk点的约束函数值
x1=2;x2=0;%xk点的坐标值
disp（'xk点约束函数数值:
'）
g=subs（[g1g2g3]）
disp（'根据g1=0和g2=0，判断g1和g2为起作用约束:
'）
%计算xk的梯度
%目标函数的梯度
gradf=jacobian（f）;%计算目标函数的梯度
disp（'目标函数的梯度'）
disp（gradf）%显示目标函数的梯度
gradfk=subs（subs（gradf,x1）,x2）%显示目标函数xk点的梯度值
%约束函数g1
gradg1=jacobian（g1）;
disp（'约束函数g1的梯度:
'）
disp（gradg1）
gradg1k=subs（subs（gradg1,x1）,x2）
%约束函数g2
gradg2=jacobian（g2,v）
disp（'约束函数g2的梯度：
'）
disp（gradg2）
gradg2k=subs（subs（gradg2,x1）,x2）
%根据kt条件建立线性方程组
A=[gradg1k
（1）,gradg2k
（1）;gradg1k
（2）,gradg2k
（2）]%建立k-T条件线性方程组的系数矩阵
b=[-gradfk
（1）;-gradfk
（2）]%建立K-T条件线性方程组的常数向量
lamda=A\b;%解线性方程组，求拉格朗日乘子
disp（'拉格朗日乘子:
'）
disp（lamda）
iflamda>=0
disp（'xk点是约束极小点'）
else
disp（'xk点不是约束极小点'）
end
disp（'目标函数最小值minf（xk）'）
minf=subs（f）%显示目标函数的最小值
M文件的运行结果如下
xk点约束函数数值：
g=00-2根据g1=0和g2=0,判断g1和g2为起作用约束：
目标函数的梯度
[2*x1-6,2*x2]
gradfk=-20约束函数g1的梯度
[2*x1,1]
gradg1k=41
gradg2=[0,-1]约束函数g2的梯度
[0,-1]
gradg2k=0-1
A=40
1-1
b=2
0
拉格朗日乘子：
0.5000
0.5000
xk点是约束极小点
目标函数最小值minf（xk）
minf=1
例3-1利用进退法求
的极值区间，取初始点0，步长为0.1
symst
f=t^4-t^2-2*t+5;
[x1,x2]=minJT（f,0,0.1）
进退法确定搜索区间函数文件minJT如下：
function[minx,maxx]=minJT（f,x0,h0,eps）
%目标函数：
f;
%初始点：
x0;
%初始步长：
h0;
%精度：
esp;
%区间左端点：
minx;
%区间右端点：
maxx;
formatlong;
ifnargin==3
esp=1.0e-6;
end
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）;
iff4x2=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
end
minx=min（x1,x3）;
maxx=x1+x3-minx;
formatshort;
M函数文件的运行结果如下
x1=
0.3000
x2=
1.5000
例3-2黄金分割法求一元函数
的极小点，设初始搜索区间
作两步迭代运算
symst;
f=t^2-10*t+36;
[x,fx]=minHJ（f,-10,10）
黄金分割法一维搜索函数文件minHJ如下：
function[x,minf]=minHJ（f,a,b,eps）
%目标函数：
f;
%搜素区间左端点：
a;
%搜索区间右端点：
b;
%精度：
eps;
%目标函数取最小值时自变量值：
x;
%目标函数的最小值：
minf;
formatlong;
ifnargin==3
eps=1.0e-6;
end
l=a+0.382*（b-a）;%试探点
u=a+0.618*（b-a）;%试探点
k=1;
tol=b-a;
whiletol>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）;%缩小搜索区间
else
b=u;%改变区间右端点
u=l;
l=a+0.382*（b-a）;%缩小搜索区间
end
k=k+1;
tol=abs（b-a）;
end
ifk==100000
disp（'找不到最优点！
'）;
x=NaN;
minf=NaN;
return;
end
x=（a+b）/2;
minf=subs（f,findsym（f）,x）;
formatshort;
M函数文件的运行结果如下：
x=5.0000
fx=11.0000
例3-3利用二次插值法求函数
的最优解，设初始搜索区间为
symst;
f=sin（t）;
[x,fx]=minPWX（f,4,5）
二次插值法一维搜索函数文件minPWX如下：
function[x,minf]=minPWX（f,a,b,eps）
%目标函数：
f;
%初始收缩区间左端点：
a;
%初始收缩区间左端点：
b;
%精度：
eps;
%目标函数取最小值时的自变量：
x;
%目标函数的最小值：
minf
formatlong;
ifnargin==3
eps=1.0e-6;
end
t0=（a+b）/2;
k=0;
tol=1;
whiletol>eps
fa=subs（f,findsym（f）,a）;%搜索区间左端点函数值
fb=subs（f,findsym（f）,b）;%搜索区间右端点函数值
ft0=subs（f,findsym（f）,t0）;%内插点函数值
tu=fa*（b^2-t0^2）+fb*（t0^2-a^2）+ft0*（a^2-b^2）;
td=fa*（b-t0）+fb*（t0-a）+ft0*（a-b）;
t1=tu/2/td;%插值多项式的极小点
ft1=subs（f,findsym（f）,t1）;%插值多项式的极小值
tol=abs（t1-t0）;
ifft1<=ft0
ift1<=t0
b=t0;%更新搜索区间右端点
t0=t1;%更新内插点
else
a=t0;%更新搜索区间左端点
t0=t1;%更新内插点
end
k=k+1;
else
ift1<=t0
a=t1;
else
b=t1;
end
k=k+1;
end
end
x=t1;
minf=subs（f,findsym（f）,x）;
formatshort;
M函数文件的运行结果如下：
x=4.7124
fx=-1.0000
例3-4试用牛顿法求函数
的极小点
symst;
f=t^4-4*t^3-6*t^2-16*t+4;
[x,fx]=minNewton（f,3）
牛顿法一维搜索函数文件minNewton如下
function[x,minf]=minNewton（f,x0,eps）
%目标函数：
f;
%初始点：
x0;
%精度：
eps;
%目标函数取最小值时的自变量：
x;
%目标函数的最小值：
minf
formatlong
ifnargin==2
eps=1.0e-6;
end
df=diff（f）;%一阶导数
d2f=diff（df）;%二阶导数
k=0;
tol=1;
whiletol>eps
dfx=subs（df,findsym（df）,x0）;%一阶导数值
d2fx=subs（d2f,findsym（d2f）,x0）;%二阶导数值
x1=x0-dfx/d2fx;
k=k+1;
tol=abs（dfx）;
x0=x1;
end
x=x1;
minf=subs（f,findsym（f）,x）;
formatshort;
M函数文件的运行结果如下：
x=4
fx=-156
例3-5用fminbnd求函数
在区间
上的极小值
[x,fval,exitflag,output]=fminbnd（'x^4-x^2+x-1',-2,1）
所得结果为
x=-0.8846
fval=-2.0548
exitflag=1
output=iterations:
11
funcCount:
14
algorithm:
'goldensectionsearch,parabolicinterpolation'
message:
[1x112char]
从输出结果可以看出，fminbnd用了黄金分割算法和抛物线算法来求本例的极小值，exitflag=1表明成功求得函数的极小值，迭代次数11次，要查看结果精度，可以接着在命令窗口中输入.
output.message
ans=Optimizationterminated:
thecurrentxsatisfiestheterminationcriteriausingOPTIONS.TolXof
1.000000e-004
说明求得的结果的精度为1.0e-4，如果想提高精度，options参数来指定，在命令窗口输入
>>opt=optimset（'Tolx',1.0e-6）;
>>formatlong;
>>[x,fval,exitflag,output]=fminbnd（'x^4-x^2+x-1',-2,1,opt）
所得结果为
x=-0.88464616447476
fval=-2.05478406218540
exitflag=1
output=iterations:
11
funcCount:
14
algorithm:
'goldensectionsearch,parabolicinterpolation'
message:
[1x112char]
这样求得的结果就有了1.0e-6的精度。
例4-2利用梯度法求目标函数
的极小值，设初始点为
，收敛精度
symsts;
f=t^2+s^2-t*s-10*t-4*s+60;
[x,mf]=minFD（f,[00],[ts]）
梯度法函数文件minFD如下：
function[x,minf]=minFD（f,x0,var,eps）
%目标函数:
f;
%初始点:
x0;
%自变量向量:
var;
%精度:
eps;
%目标函数取最小值时的自变量值:
x;
%目标函数的最小值:
minf
formatlong;
ifnargin==3
eps=1.0e-6;
end
symsl;
tol=1;
gradf=-jacobian（f,var）;
whiletol>eps
v=subs（gradf,var,x0）;
tol=norm（v）;
y=x0+l*v;
yf=subs（f,var,y）;
[a,b]=minJT（yf,0,0.1）;
xm=minHJ（yf,a,b）;%用黄金分割法进行一维搜索
x1=x0+xm*v;
x0=x1;
end
x=x1;
minf=subs（subs（f,x
（1））,x
（2））;
formatshort;
M函数文件的运行结果如下：
x=8.00006.0000
mf=8.0000
例4-3函数
的初始点取
，用牛顿法求他的极值点
symsts;
f=t^2-4*s^2;
[x,mf]=minNT（f,[11],[ts]）
牛顿法函数文件minNT如下
function[x,minf]=minNT（f,x0,var,eps
%目标函数：
f;
%初始点：
x0;
%自变量向量var;
%精度：
eps;
%目标函数取最小时的自变量值：
x;
%目标函数最小值：
minf;
formatlong;
ifnargin==3
eps=1.0e-6;
end
tol=1;
x0=transpose（x0）;
gradf=jacobian（f,var）;%梯度方向
jacf=jacobian（gradf,var）;%雅克比矩阵
whiletol>eps
v=subs（gradf,var,x0）;
tol=norm（v）;
pv=subs（jacf,var,x0）;
p=-inv（pv）*transpose（v）;%搜索方向
p=double（p）;
x1=x0+p;
x0=x1;
end
x=x1;
minf=subs（f,var,x）;
formatshort;
M函数文件的运行结果如下：
x=0
0
mf=0
升级会员 