MATLB机械优化设计程序XXXX.docx

《MATLB机械优化设计程序XXXX.docx》由会员分享，可在线阅读，更多相关《MATLB机械优化设计程序XXXX.docx（32页珍藏版）》请在三一办公上搜索。

1、例2-1求 在点和点的梯度。%例2-1梯度的计算syms x1 x2 %定义符号变量f=x12+x22-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(

2、1),Xzuobiao2(2) %计算Xzuobiao1点的梯度值gmk2=norm(gfk2) %计算Xzuobiao1点的梯度模igk2=gfk2/gmk2 %计算Xzuobiao1点的梯度单位向量gradf = 2*x1-4, 2*x2gfk1 = 2 4gmk1 = 4.4721igk1 = 0.4472 0.8944gfk2 = 0 0gmk2 = 0Warning: Divide by zero.igk2 = NaN NaN例2-2把函数在点展开泰勒二次近似式%例2-2 Taylor展开 syms x1 x2f=4+4.5*x1-4*x2+x12+2*x22-2*x1*x2+x14

3、-2*x12*x2disp(函数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=dx1

4、x1 dx1x2;dx2x1 dx2x2;pretty(simplify(H);%计算xk点的值x1=2.0; x2=2.5;disp(函数在xk点的函数值:)fk=subs(f)disp(函数在xk点的一节偏导数矩阵:)dk=subs(dx1 dx2)disp(函数xk点的海色矩阵:)HK=subs(dx1x1 dx1x2;dx2x1 dx2x2)disp(函数在xk点的二阶Taylor展开式:)syms x1 x2fkTL=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/

5、2*x1-4*x2+x12+2*x22-2*x1*x2+x14-2*x12*x2 函数f的表达式: 2 2 4 2 4 + 9/2 x1 - 4 x2 + x1 + 2 x2 - 2 x1 x2 + x1 - 2 x1 x2函数f的一阶偏导数表达式: 3 9/2 + 2 x1 - 2 x2 + 4 x1 - 4 x1 x2 2 -4 + 4 x2 - 2 x1 - 2 x1函数f的二阶偏导数表达式： 3 9/2 + 2 x1 - 2 x2 + 4 x1 - 4 x1 x2 2 2 + 12 x1 - 4 x2 -2 - 4 x1 -2 - 4 x1 4 函数在xk点的函数值:fk = 5.50

6、00函数在xk点的一节偏导数矩阵:dk = 15.5000 -6.0000函数xk点的海色矩阵:HK = 40 -10 -10 4函数在xk点的二阶Taylor展开式: 2 2 32 - 79/2 x1 + 4 x2 + 20 x1 - 10 x1 x2 + 2 x2例2-3求函数的极值点和极值%例2-3 求函数的极值syms x1 x2 x3f=2*x12+5*x22+x32+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);d

7、sx3=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

8、的2阶偏导数，构成海色矩阵disp(函数f的2阶偏导数矩阵)H=dsx1x1 dsx1x2 dsx1x3;dsx2x1 dsx2x2 dsx2x3;dsx3x1 dsx3x2 dsx3x3%计算海色矩阵的正定性D,p=chol(subs(H);if p=0;disp(海色矩阵为正定，函数f有极小点：);end%计算极值存在的必要条件，求极值点坐标x1,x2,x3=solve(dsx1,dsx2,dsx3,x1,x2,x3);disp(极值点坐标:)fprintf(1,x1=%3.4fn,subs(x1);fprintf(1,x2=%3.4fn,subs(x2);fprintf(1,x3=%3.

9、4fn,subs(x3);disp(在极值点，函数f数值：)fmb=subs(f)M文件的运行结果如下函数f的表达式: 2 2 2 2 x1 + 5 x2 + x3 + 2 x2 x3 + 2 x1 x3 - 6 x2 + 3函数f的1阶偏导数： 4 x1 + 2 x3 10 x2 + 2 x3 - 6 2 x3 + 2 x2 + 2 x1函数f的2阶偏导数矩阵 H = 4, 0, 2 0, 10, 2 2, 2, 2海色矩阵为正定，函数f有极小点：极值点坐标:x1=1.0000x2=1.0000x3=-2.0000在极值点，函数f数值： fmb = 0例2-5 已知二维约束问题受约束为 例2

10、-5 MATLAB实现，用M文件判别函数的凸性：%例2-5判别函数的凸性syms x1 x2f=60-10*x1-4*x2+x12+x22-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

11、的2阶偏导数，构成海色矩阵disp(函数f的2阶偏导数矩阵)H=dsx1x1 dsx1x2;dsx2x1 dsx2x2%计算函数矩阵的正定性d,p=chol(subs(H);if p=0;disp(海色矩阵为正定，函数f为凸函数);endM文件的运行结果如下函数f的表达式: 2 2 60 - 10 x1 - 4 x2 + x1 + x2 - x1 x2函数f的1阶偏导数: -10 + 2 x1 - x2 -4 + 2 x2 - x1函数f的2阶偏导数矩阵H = 2, -1 -1, 2海色矩阵为正定，函数f为凸函数%例2-6K-T条件syms x1 x2 v %定义目标函数和约束函数的符号变量%

12、目标函数和约束函数f=(x1-3)2+x22; %目标函数f的表达式g1=x12+x2-4;g2=-x2;g3=-x1;v=x1,x2;%计算xk点的约束函数值x1=2;x2=0; %xk点的坐标值disp(xk点约束函数数值:)g=subs(g1 g2 g3)disp(根据g1=0和g2=0，判断g1和g2为起作用约束:)%计算xk的梯度%目标函数的梯度gradf=jacobian(f); %计算目标函数的梯度disp(目标函数的梯度)disp(gradf) %显示目标函数的梯度gradfk=subs(subs(gradf,x1),x2)%显示目标函数xk点的梯度值%约束函数g1gradg1

13、=jacobian(g1);disp(约束函数g1的梯度:)disp(gradg1)gradg1k=subs(subs(gradg1,x1),x2)%约束函数g2gradg2=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=Ab; %解线性方程组，求

14、拉格朗日乘子disp(拉格朗日乘子:)disp(lamda)if lamda=0 disp(xk点是约束极小点)else disp(xk点不是约束极小点)enddisp(目标函数最小值minf(xk)minf=subs(f) %显示目标函数的最小值M文件的运行结果如下xk点约束函数数值：g = 0 0 -2根据g1=0和g2=0,判断g1和g2为起作用约束：目标函数的梯度 2*x1-6, 2*x2 gradfk = -2 0 约束函数g1的梯度 2*x1, 1gradg1k = 4 1gradg2 = 0, -1 约束函数g2的梯度 0, -1gradg2k = 0 -1A = 4 0 1 -

15、1b = 2 0拉格朗日乘子： 0.5000 0.5000xk点是约束极小点目标函数最小值minf(xk)minf = 1例3-1 利用进退法求的极值区间，取初始点0，步长为0.1syms tf=t4-t2-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;format long;if nargin=3 esp=1.0e-6;endx1=x0;k=0;h=h0;wh

16、ile 1 x4=x1+h; %试探步k=k+1; f4=subs(f,findsym(f),x4); f1=subs(f,findsym(f),x1); if f4eps & kfu 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 if k=100000 disp(找不到最优点！);x=NaN; minf=NaN; return; end x=(a+b)/2; minf=subs(f,findsym(f

17、),x); format short;M函数文件的运行结果如下：x = 5.0000fx =11.0000例3-3利用二次插值法求函数的最优解，设初始搜索区间为syms t;f=sin(t);x,fx=minPWX(f,4,5)二次插值法一维搜索函数文件minPWX如下：function x,minf=minPWX(f,a,b,eps)%目标函数：f;%初始收缩区间左端点：a;%初始收缩区间左端点：b;%精度：eps;%目标函数取最小值时的自变量：x;%目标函数的最小值：minfformat long;if nargin=3 eps=1.0e-6;endt0=(a+b)/2;k=0;tol=1

18、; while toleps fa=subs(f,findsym(f),a); %搜索区间左端点函数值 fb=subs(f,findsym(f),b); %搜索区间右端点函数值 ft0=subs(f,findsym(f),t0); %内插点函数值 tu=fa*(b2-t02)+fb*(t02-a2)+ft0*(a2-b2); td=fa*(b-t0)+fb*(t0-a)+ft0*(a-b);t1=tu/2/td; %插值多项式的极小点 ft1=subs(f,findsym(f),t1); %插值多项式的极小值 tol=abs(t1-t0); if ft1=ft0 if t1=t0 b=t0;

19、%更新搜索区间右端点 t0=t1; %更新内插点 else a=t0; %更新搜索区间左端点 t0=t1; %更新内插点 end k=k+1; else if t1eps 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); format short;M函数文件的运行结果如下：x =4fx = -156例3-5用fminbnd求函数在区间上的

20、极小值x,fval,exitflag,output=fminbnd(x4-x2+x-1,-2,1)所得结果为x = -0.8846fval = -2.0548exitflag =1output = iterations: 11 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 1x112 char从输出结果可以看出，fminbnd用了黄金分割算法和抛物线算法来求本例的极小值，exitflag =1表明成功求得函数的极小值，迭代次数11次，要查看结果精度，可以接着在命令窗口中输入.ou

21、tput.messageans =Optimization terminated:the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004说明求得的结果的精度为1.0e-4，如果想提高精度，options参数来指定，在命令窗口输入 opt=optimset(Tolx,1.0e-6); format long; x,fval,exitflag,output=fminbnd(x4-x2+x-1,-2,1,opt)所得结果为x = -0.88464616447476fval = -2.05

22、478406218540exitflag =1output = iterations: 11 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 1x112 char这样求得的结果就有了1.0e-6的精度。 例4-2 利用梯度法求目标函数的极小值，设初始点为，收敛精度syms t s;f=t2+s2-t*s-10*t-4*s+60;x,mf=minFD(f,0 0,t s)梯度法函数文件minFD如下：function x,minf = minFD(f,x0,var,eps)%目标函数

23、:f;%初始点:x0;%自变量向量:var;%精度:eps;%目标函数取最小值时的自变量值:x;%目标函数的最小值:minfformat long;if nargin = 3 eps = 1.0e-6;endsyms l;tol = 1;gradf = - jacobian(f,var);while toleps 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

24、*v; x0 = x1;endx = x1;minf = subs(subs(f,x(1),x(2);format short;M函数文件的运行结果如下：x = 8.0000 6.0000mf =8.0000例4-3函数的初始点取，用牛顿法求他的极值点syms t s;f=t2-4*s2;x,mf=minNT(f,1 1,t s) 牛顿法函数文件minNT如下 function x,minf = minNT(f,x0,var,eps%目标函数：f;%初始点：x0;%自变量向量var;%精度：eps;%目标函数取最小时的自变量值：x;%目标函数最小值：minf;format long;if na

25、rgin = 3 eps = 1.0e-6;endtol = 1;x0 = transpose(x0);gradf = jacobian(f,var); %梯度方向jacf = jacobian(gradf,var); %雅克比矩阵 while toleps 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;endx = x1;minf = subs(f,var,x);format

26、short;M函数文件的运行结果如下：x = 0 0mf = 0例4-4给定初始点，用阻尼牛顿法求函数的极小点syms t s z;f=(t-s+z)2+(-t+s+z)2+(t+s+z)2x,mf=minMNT(f,0.5 1 0.5,t s z)阻尼牛顿法函数文件minMNT如下function x,minf = minMNT(f,x0,var,eps)format long;if nargin = 3 eps = 1.0e-6;endtol = 1;x0 = transpose(x0);syms l;gradf = jacobian(f,var);jacf = jacobian(grad

27、f,var);while toleps v = subs(gradf,var,x0); tol = norm(v); pv = subs(jacf,var,x0); p = -inv(pv)*transpose(v); y = x0 + l*p; yf = subs(f,var,y);a,b = minJT(yf,0,0.1); %进退法求单峰区间 xm = minHJ(yf,a,b); %黄金分割法进行一维搜素 x1 = x0 + xm*p; x0 = x1;endx = x1;minf = subs(f,var,x);format short;M函数文件的运行结果如下：x = 1.0e-0

28、15 *-0.3468 -0.6936 -0.3468 mf =2.4053e-030例4-5 给定初始点用共轭梯度法求解syms t s;f=t2+4*s2;x,mf=minGETD(f,1 1,t s)共轭梯度法函数文件minGETD如下function x, minf = minGETD (f,x0,var,eps)format long;if nargin=3 eps=1.0e-6;endx0 = transpose(x0);n = length(var);syms l;gradf = jacobian(f,var);v0=subs(gradf,var,x0);p = -transpo

29、se(v0);k = 0;while 1 v = subs(gradf,var,x0); tol = norm(v); if tol=eps x = x0; break; end y = x0 + l*p; yf = subs(f,var,y); a,b = minJT(yf,0,0.1); %进退法确定单峰区间 xm = minPWX(yf,a,b); %二次插值一维搜素 x1 = x0 + xm*p; vk = subs(gradf,var,x1); tol = norm(vk);if tol=eps x = x1; break; end if k+1=n x0 = x1; continue;else lamda = dot(vk,vk)/dot(v,v); p = -transpose(vk) + lamda*p; k = k+1;