MATLB机械优化设计程序XXXX.docx
- 文档编号:25395302
- 上传时间:2023-06-08
- 格式:DOCX
- 页数:47
- 大小:87.43KB
MATLB机械优化设计程序XXXX.docx
《MATLB机械优化设计程序XXXX.docx》由会员分享,可在线阅读,更多相关《MATLB机械优化设计程序XXXX.docx(47页珍藏版)》请在冰豆网上搜索。
MATLB机械优化设计程序XXXX
MATLB机械优化设计程序XXXX
%例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);
iff4 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 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-4给定初始点 ,用阻尼牛顿法求函数 的极小点 sym
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLB 机械 优化 设计 程序 XXXX