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
例4-