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

例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(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-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=dx1x1 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/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.5000函数在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);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=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.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-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的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 %定义目标函数和约束函数的符号变量%目标函数和约束函数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=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; %解线性方程组，求拉格朗日乘子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 -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;while 1 x4=x1+h; %试探步k=k+1; f4=subs(f,findsym(f),x4); f1=subs(f,findsym(f),x1); if f4<f1 x2=x1; x1=x4; f2=f1; f1=f4; h=2*h; %加大步长 else if k=1 h=-h; %反向搜索 x2=x4; f2=f4; else x3=x2; x2=x1;x1=x4; break; end endendminx=min(x1,x3);maxx=x1+x3-minx;format short;M函数文件的运行结果如下x1 = 0.3000x2 = 1.5000例3-2 黄金分割法求一元函数的极小点，设初始搜索区间作两步迭代运算syms t;f=t2-10*t+36;x,fx=minHJ(f,-10,10)黄金分割法一维搜索函数文件minHJ如下：function x,minf=minHJ(f,a,b,eps)%目标函数：f;%搜素区间左端点：a;%搜索区间右端点：b;%精度：eps;%目标函数取最小值时自变量值：x;%目标函数的最小值：minf;format long;if nargin=3 eps=1.0e-6;endl=a+0.382*(b-a); %试探点u=a+0.618*(b-a); %试探点k=1;tol=b-a;while tol>eps && k<100000 fl=subs(f,findsym(f),l); %试探点函数值 fu=subs(f,findsym(f),u); %试探点函数值 if fl>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 if k=100000 disp('找不到最优点！');x=NaN; minf=NaN; return; end x=(a+b)/2; minf=subs(f,findsym(f),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; while tol>eps 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; %更新搜索区间右端点 t0=t1; %更新内插点 else a=t0; %更新搜索区间左端点 t0=t1; %更新内插点 end k=k+1; else if t1<=t0 a=t1; else b=t1; end k=k+1; end end x=t1; minf=subs(f,findsym(f),x);format short;M函数文件的运行结果如下：x = 4.7124fx = -1.0000例3-4试用牛顿法求函数的极小点syms t;f=t4-4*t3-6*t2-16*t+4;x,fx=minNewton(f,3)牛顿法一维搜索函数文件minNewton如下function x,minf=minNewton(f,x0,eps)%目标函数：f;%初始点：x0;%精度：eps;%目标函数取最小值时的自变量：x;%目标函数的最小值：minfformat long if nargin=2 eps=1.0e-6; end df=diff(f); %一阶导数d2f=diff(df); %二阶导数 k=0; tol=1; while tol>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); format short;M函数文件的运行结果如下：x =4fx = -156例3-5用fminbnd求函数在区间上的极小值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次，要查看结果精度，可以接着在命令窗口中输入.output.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.05478406218540exitflag =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)%目标函数:f;%初始点:x0;%自变量向量:var;%精度:eps;%目标函数取最小值时的自变量值:x;%目标函数的最小值:minfformat long;if nargin = 3 eps = 1.0e-6;endsyms l;tol = 1;gradf = - jacobian(f,var);while tol>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;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 nargin = 3 eps = 1.0e-6;endtol = 1;x0 = transpose(x0);gradf = jacobian(f,var); %梯度方向jacf = jacobian(gradf,var); %雅克比矩阵 while tol>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;endx = x1;minf = subs(f,var,x);format 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(gradf,var);while tol>eps 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-015 *-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 = -transpose(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;