MATLB机械优化设计程序XXXX

例2-1求22121()44fXxxx在点(1)3,2TX和点(2)2,0TX的梯度。%例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=24gmk1=4.4721igk1=0.44720.8944gfk2=00gmk2=0Warning:Dividebyzero.igk2=NaNNaN例2-2把函数2242121212112()44.54222fXxxxxxxxxx在点()[2.0,2.5]kTX展开泰勒二次近似式111212222.02.0401011311(,),62.0,2.52.52.5104222xxxxxxxx221212123272/2420210xxxxxx%例2-2Taylor展开symsx1x2f=4+4.5*x1-4*x2+x1^2+2*x2^2-2*x1*x2+x1^4-2*x1^2*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=[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展开式:')symsx1x2fkTL=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的表达式:22424+9/2x1-4x2+x1+2x2-2x1x2+x1-2x1x2函数f的一阶偏导数表达式:39/2+2x1-2x2+4x1-4x1x22-4+4x2-2x1-2x1函数f的二阶偏导数表达式：39/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展开式:2232-79/2x1+4x2+20x1-10x1x2+2x2例2-3求函数22212323312()252263fXxxxxxxxx的极值点和极值%例2-3求函数的极值symsx1x2x3f=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的表达式:2222x1+5x2+x3+2x2x3+2x1x3-6x2+3函数f的1阶偏导数：4x1+2x310x2+2x3-62x3+2x2+2x1函数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已知二维约束问题2212min()(3)fXxx受约束为2112()40gXxx22()0gXx*[2,0]TX例2-5MATLAB实现，用M文件判别函数的凸性：%例2-5判别函数的凸性symsx1x2f=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为凸函数');endM文件的运行结果如下函数f的表达式:2260-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点的梯度值%约束函数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=A\b;%解线性方程组，求拉格朗日乘子disp('拉格朗日乘子:')disp(lamda)iflamda=0disp('xk点是约束极小点')elsedisp('xk点不是约束极小点')enddisp('目标函数最小值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=41gradg2=[0,-1]约束函数g2的梯度[0,-1]gradg2k=0-1A=401-1b=20拉格朗日乘子：0.50000.5000xk点是约束极小点目标函数最小值minf(xk)minf=1例3-1利用进退法求42()25ftttt的极值区间，取初始点0，步长为0.1symstf=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==3esp=1.0e-6;endx1=x0;k=0;h=h0;while1x4=x1+h;%试探步k=k+1;f4=subs(f,findsym(f),x4);f1=subs(f,findsym(f),x1);iff4f1x2=x1;x1=x4;f2=f1;f1=f4;h=2*h;%加大步长elseifk==1h=-h;%反向搜索x2=x4;f2=f4;elsex3=x2;x2=x1;x1=x4;break;endendendminx=min(x1,x3);maxx=x1+x3-minx;formatshort;M函数文件的运行结果如下x1=0.3000x2=1.5000例3-2黄金分割法求一元函数2()1036fxxx的极小点，设初始搜索区间[,][10,10]ab作两步迭代运算symst;f=