优化设计上机指导书2012

上机指导书

适用专业: 机械工程 课程代码: 8201581 编写单位: 机械工程学院 编 写 人: 杨昌明 审 核 人: 审 批 人:

《优化设计及有限元分析》上机指导书v2011

目 录

实验一 （实验代码1）„„„„„„„„„„„„„„„„„„„„„„„„„„ 1 实验二 （实验代码2）„„„„„„„„„„„„„„„„„„„„„„„„„„ 9 注释 „„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„ 14 主要参考文献

15 - 1 -

„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„ 《优化设计及有限元分析》上机指导书v2011

MATLAB Optimization Toolbox应用基础

专题Ⅰ 单变量函数最小值问题（一维问题）

一、上机目的：

1、掌握一维搜索方法的基本概念、算法等知识点； 2、熟悉Matlab软件环境及有关语句和函数的使用方法。

二、原理和方法：

函数fminbnd是用来寻找单变量函数在固定区间内的最小值点及最小值。MATLAB各工具箱的每个函数，包括优化工具箱函数，都有多种不同的调用格式。其调用格式可使用在线帮助系统查询。

函数fminbnd最常用的调用格式为： [x, fval] = fminbnd(fun, x1, x2)

%返回函数fun在区间x1三、上机内容与步骤：

1、熟悉MATLAB软件环境，了解用户界面、掌握帮助系统使用等（重点）。

2、如求解函数f(x)=sin(x)在区间[0, 2]内的最小值及最小值点，使用如下命令即可：

>> [x, fval] = fminbnd(@sin, 0, 2*pi)

%符号“@”表明目标函数为MATLAB自定义的正弦函数sin.m。 x =

4.7124 %最小值点赋值给x。

fval = -1.0000 %最小值赋值给fval。

3、如果目标函数并不是MATLAB自定义的函数，需要首先定义一个函数M文件，以该文件的名称来代替上述命令行中的“sin”部分。

如希望求出f(x)x3cosxxlogxex在[0，1]内的最小值及最小值点，首先编写目标函

数的M文件opt1.m，如下：

function f = myfun(x)

%编写函数M文件时，上句代码中的“myfun”并不一定与文件名相同。 f = (x^3+cos(x)+x*log(x))/exp(x)

然后调用fminbnd函数，在命令窗口中输入： >> [x, fval]=fminbnd(@opt1, 0, 1) 

x =

0.5223 fval =

0.3974

%注意此处使用“opt1”而不是“myfun”。

4、如果目标函数的表达式较为简单，也可以使用表达式的字符串而不必编写目标函数。如上例也可以使用[x, fval]=fminbnd('(x^3+cos(x)+x*log(x))/exp(x)',0,1)实现。

注：在MATLAB软件中，function为关键字，%为注释符号，>>为命令窗口中的提示符。

- 2 -

《优化设计及有限元分析》上机指导书v2011

专题Ⅱ 线性方程组的应用 一、上机目的：

1、了解线性规划问题及可行解、最优解的概念 ； 2、掌握Matlab软件关于求解线性规划的语句和方法。 二、原理和方法：

在工程实践中，很多重要的实际问题都是线性的（至少能够用线性函数很好的近似表示），所以我们一般把这些问题化为线性的目标函数和约束条件进行分析，通常将目标函数和约束都是线性表达式的规划问题称为线性规划。

其一般形式是：

矩阵形式表示为：

mins.t.fcxAxb,x0Tminfc1x1c2x2cnxna11x1a12x2a1nxnb1a21x1a22x2a2nxnb2axaxaxbm22mnnmm11xi0(i1,2,,n)s.t.线性规划的可行解是满足约束条件的解；

线性规划的最优解是使目标函数达到最优的可行解。

一般求解线性规划的常用方法是单纯形法和改进的单纯形法，这类方法的基本思路是先求得一个可行解，检验是否为最优解；若不是，可用迭代的方法找到另一个更优的可行解，经过有限次迭代后，可以找到可行解中的最优解或者判定无最优解。 三、上机内容与步骤：

在Matlab优化工具箱中，linprog函数是使用单纯形法求解下述线性规划问题的函数。

命令格式为：

[x,fval]linprog(c,A,b,aeq,beq,vlb,vub)[x,fval]linprog(c,A,b,aeq,beq,vlb,vub,x0)mins.t.fcxAxb,aeqxbeq;vlbxvubT其中：A为不等式约束条件矩阵，b, c分别为目标函数的系数向量和不等式约束条件中最右边的数值向量；也可设置解向量的上界vlb和下界vub，即解向量必须满足

- 3 -

《优化设计及有限元分析》上机指导书v2011

vlb<=x<=vub；还可预先设置初始解向量x0。

如没有不等式约束，而只有等式约束条件时，A=[ ],b=[ ]; aeq与beq分别为等式约束矩阵；

输出的结果：x表示最优解向量；fval表示最优值。

某一线性规划问题，如下：

c = [-3;1;1]; A = [1 -2 1; 4 -1 -2]; b = [11;-3]; aeq = [2 0 -1]; beq = -1; vlb = [0;0;0];

[x,fval] = linprog(c, A, b, aeq, beq, vlb)

则在命令窗口中输入：optest，可得到计算结果：x1 = 4,x2 = 1,x3 = 9，fval = 2。

- 4 -

maxf3x1x2x3x12x2x3114x1x22x332x1x31x0,i1,2,3is.t.考虑到linprog函数只解决形式如下的线性规划。

mins.t.fcxAxb,aeqxbeq;x0T所以先要将线性规划变换为如下形式：

minf3x1x2x32x1x31x12x2x3114x1x22x33x0,i1,2,3is.t.建立M文件如下（如以文件名optest保存）：

《优化设计及有限元分析》上机指导书v2011

专题Ⅲ 无约束规划计算方法 一、上机目的

1、了解无约束规划问题的求解原理与方法 ； 2、会用Matlab软件求解无约束规划问题。 二、原理和方法

无约束规划问题的解法一般按目标函数的形式分为两大类： 一类是一元函数的一维搜索法，如黄金分割法、插值法等； 另一类是求解多元函数的下降迭代法。 三、上机内容与步骤

在Matlab软件中，求解无约束规划的常用命令是： x=fminunc(‘fun’, x0)

其中，fun函数应预先定义到M文件中。

问题如下：

初始点x0 = [-2 4] 首先建立函数文件fun02.m

function f = fun02(x)

f = 3/2*x(1)^2 + 1/2*x(2)^2 - x(1)*x(2) - 2*x(1) 在命令窗口输入：

x0=[-2;4]

x=fminunc(@fun02,x0)

计算结果：极小值为-1，最优点：x1=1,x2=1。

专题Ⅳ 约束非线性规划计算方法 一、上机目的

1、了解约束非线性规划问题的求解原理与方法； 2、会用Matlab软件求解约束非线性规划问题。 二、原理和方法

对于约束非线性规划，随着目标函数和约束条件的不同，解法也不同，一般来说，有两类方法：

（1）、将约束问题转化为无约束问题的求解方法； （2）、用线性规划来逼近非线性规划； 三、上机内容与步骤

约束非线性规划的一般形式为：

- 5 -

minf(x)32x1212x2x1x22x12 《优化设计及有限元分析》上机指导书v2011

s.tminf(x)xAxb,aeq*xbeq(线性约束)）g(x)0,ceq(x)0（非线性约束lbxub其中，f(x)为多元实值函数；g(x)为向量函数，并且f(x)，g(x)中至少有一个函数是非线性函数的（否则成为线性规划问题）。

在Matlab优化工具箱中，fmincon函数是用SQP算法来解决一般的约束非线性规划的函数，其命令格式为：

x = fmincon(‘fun’, x0, A, b)

x = fmincon(‘fun’, x0, A, b, Aeq, beq) x = fmincon(‘fun’, x0, A, b, Aeq, beq, lb, ub) x = fmincon(‘fun’, x0, A, b, Aeq, beq, lb,ub, nonlcon)

对于如下约束非线性规划问题：

min(x11)(x22)(x33)(x44)2222x1x2x3x45s.t.3x13x22x3x410x0i初值为[1;1;1;1]建立一个M文件 fun05.m，代码如下：

function y=fun05(x)

y=(x(1)-1)^2+(x(2)-2)^2+(x(3)-3)^2+(x(4)-4)^2; 存储为fun05.m文件. 在命令窗口输入：

x0 = [1;1;1;1]; A= [1 1 1 1;3 3 2 1]; B = [5;10]; Aeq = []; Beq = []; Lb = [0;0;0;0];

[x, g] = fmincon(@fun05,x0,A,B,Aeq,Beq,Lb) 计算结果为： X = 0.0000 0.6667 1.6665 2.6668

g = 6.3333

- 6 -

《优化设计及有限元分析》上机指导书v2011

实验一 齿轮传动系统优化设计

一、上机目的和任务

1、 机械优化设计是一门实践性很强的课程，通过实际上机计算达到以下目的： 1）加深对机械优化设计方法的基本理论和算法的理解。 2）培养学生独立编制计算机程序的能力。 3）掌握常用优化方法程序的使用方法。

4）培养学生灵活运用优化设计方法解决工程实际问题的能力。

二、上机内容

某二级斜齿圆柱齿轮减速器（如图1所示），高速级输入功率 P1 =6.2KW，转速n1=1450r/min;总传动比i=31.5，齿轮宽度系数a=0.4，齿轮材料和热处理：大齿轮45号钢正火187~207HB，小齿轮45号钢调质228~255HB，工作寿命10年以上。要求按照总中心距 a最小来确定齿轮传动方案。

图1 传动系统简图

三、上机步骤

该优化问题属于约束非线性规划问题，采用MATLAB优化工具箱进行求解。 求约束非线性规划问题的MATLAB函数是fmincon，约束非线性规划问题的数学模型表示为：

minf(X) s.t. AX  b(线性不等式约束) AeqX = beq（线性等式约束）

C（X）0（非线性不等式约束）

Ceq（X）=0（非线性等式约束） Lb X Ub（边界约束）

函数fmincon的使用格式为：

[ x, fval, exitflag, output, hession ]=

fmincon( @fun, x0, A, b, Aeq, beq, Lb, Ub, ’NIc’, options，P1, P2, …)

- 7 -

《优化设计及有限元分析》上机指导书v2011

其中，输出参数有：

x是返回目标函数的最优解；

fval是返回目标函数在最优解x点的函数值； exitflag是返回算法的终止标志；

output是返回优化算法信息的一个数据结构； grad是返回目标函数在最优解x点的梯度；

hession是返回目标函数在最优解x点的hession矩阵值。

输入参数有： fun是调用目标函数的函数文件名；

x0是初始点；

线性不等式约束约束条件的系数矩阵A和常数向量b； 线性等式约束条件的系数矩阵Aeq和常数向量beq；

设计变量X的下界向量Lb和上界向量Ub；

’NIc’是定义非线性约束条件的函数名； options是设置优化选项参数；

P1,P2等是传递给fun的附加参数。

参数A, b, Aeq, beq, Lb, Ub, options如果没有定义，可用空矩阵符号“[ ]”代替。

(一) 建立优化设计数学模型

将涉及总中心距a齿轮传动方案的6个独立参数作为设计变量

TT X = [mn1,mn2,z1,z3,i1,][x1,x2,x3,x4,x5,x6]

式中，mn1,mn2分别为高速级和低速级齿轮副的模数；

z1,z3分别为高速级和低速级小齿轮齿数；

i1为高速级传动比；

为齿轮副螺旋角。

减速器总中心距a最小的目标函数为：

() minfXx1x3(1x5)x2x4(131.x5512coxs6

)性能约束包括：齿轮接触强度条件、齿根弯曲强度条件、高速级大齿轮与低速轴不干涉条件等。根据齿轮材料与热处理规范，得到齿面许用接触应力[H]=518.75MPa，齿根许用弯曲应力

[F]1,3=153.5MPa和[F]2,4=141.6MPa.。根据传递功率和转速，在齿轮强度计算条件中代入有关数据：

高速轴转矩T1 =41690Nmm，中间轴转矩 T2 =40440 Nmm，高速轴和低速轴载荷系数 K1=1.225和

K2=1.204。

边界约束包括：根据传递功率和转速估计高速级和低速级齿轮副模数的范围；综合考虑传动平稳、

- 8 -

《优化设计及有限元分析》上机指导书v2011

轴向力不能太大、轴齿轮的分度圆直径不能太小与两级传动的大齿轮浸油深度大致相近等因素，估计两级传动大齿轮的齿数范围、高速级传动比范围和齿轮副螺旋角范围等。

因此，建立如下17个不等式约束条件：

g1(X)cosx63.079102336x1x3x50 （高速级齿轮接触强度条件）

433g2(X)x5cosx61.07110g3(X)cosx69.939102225x2x40 （低速级齿轮接触强度条件）

3233(1x5)x1x30 （高速级大齿轮弯曲强度条件）

4g4(X)x5cosx61.07610(31.5x5)x2x40（低速级大齿轮弯曲强度条件）

32 g5(X)x5[2(x150)cosx6x1x3x5]x2x4(31.5x5)0 （大齿轮与轴不干涉条件）

g6(X)2x10 （高速级齿轮副模数的下限） g7(X)x150 （高速级齿轮副模数的上限） g8(X)3.5x20 （低速级齿轮副模数的下限） g9(X)x260 （低速级齿轮副模数的上限） g10(X)14x30 （高速级小齿轮齿数的下限） g11(X)x3220 （高速级小齿轮齿数的上限） g12(X)16x40 (低速级齿轮数的下限) g13(X)x4220 （低速级齿轮齿数的上限) g14(X)5.8x50 (高速级传动比的下限) g15(X)x570 (高速级传动比的上限) g16(X)8x60 (齿轮副螺旋角的下限) g17(X)x6150 （齿轮副螺旋角的上限)

（二）编制优化设计的M文件（主程序设计、目标函数子程序设计、非线性约束子程序设计）

% 两级斜齿轮传动中心距优化设计

% 1-----减速器中心距优化设计主程序 % 设计变量的初始值

x0=[2；4；18；20；6.4；10]； % 设计变量的下界与上界 lb=[2；3.5；14；16；5.8；8]； ub=[5；6；22；22；7；15]；

- 9 -

《优化设计及有限元分析》上机指导书v2011

% 使用多维约束优化命令fmincon（调用目标函数jsqyh_f和非线性约束函数jsqyh_g） % 无线性不等式约束，设计变量的系数矩阵a=[]和常数项向量b=[]不需定义 % 没有等式约束，则参数：系数矩阵Aeq=[]和常数项向量beq=[]不需定义

[x，fn]= fmincon（@jsqyh_f,x0,[ ],[ ],[ ],[ ],lb,ub,@jsqyh_g）;

disp ＇******** 两级斜齿轮传动中心距优化设计最优解 ********＇ fprintf(1，＇ 高速级齿轮副模数 Mn1 = %3.4f mm\\n＇,x(1))

fprintf(1，＇ 低速级齿轮副模数 Mn2 = %3.4f mm\\n＇,x(2)) fprintf(1，＇ 高速级小齿轮齿数 z1 = %3.4f \\n＇,x(3))

fprintf(1，＇ 低速级小齿轮齿数 z3 = %3.4f \\n＇, x(4)) fprintf(1，＇ 高速级齿轮副传动比 i1 = %3.4f \\n＇,x(5)) fprintf(1，＇ 齿轮副螺旋角 beta = %3.4f 度\\n＇,x(6)) fprintf(1，＇ 减速器总中心距 a12 = %3.4f mm\\n＇,fn)

% 调用多维约束化非线性约束函数（jsqyh_g）计算最优点x*的性能约束函数值

g=jsqyh_g(x);

disp ＇ ---------- 最优点的性能约束函数值 -----------＇ fprintf(1，＇ 高速级齿轮副接触疲劳强度约束函数值 g1= %3.4f\\n＇,g(1)) fprintf(1，＇ 低速级齿轮副接触疲劳强度约束函数值 g2= %3.4f\\n＇,g(2)) fprintf(1，＇ 高速级大齿轮齿根弯曲强度约束函数值 g3= %3.4f\\n＇,g(3)) fprintf(1，＇ 低速级大齿轮齿根弯曲强度约束函数值 g4= %3.4f\\n＇,g(4)) fprintf(1，＇ 大齿轮齿顶与轴不干涉几何约束函数值 g5= %3.4f\\n＇,g(5))

% 2------两级斜齿轮减速器总中心距的目标函数（jsqyh_f） function f=jsqyh_f(x); hd=pi/180;

a1=x(1)*x(3)*(1+x(5)); a2=x(2)*x(4)*(1+31.5/x(5)); cb=2*cos(x(6)*hd); f=(a1+a2)/cb;

% 3--------两级斜齿轮减速器优化设计的非线性不等式的约束函数（jsqyh_g） function [g,ceq]=jsqyh_g(x); hd=pi/180;

g(1)=cos(x(6)*hd)^3-3.079e-6*x(1)^3*x(3)^3*x(5); g(2)=x(5)^2*cos(x(6)*hd)^3-1.701e-4*x(2)^3*x(4)^3;

g(3)=cos(x(6)*hd)^2-9.939e-5*(1+x(5))*x(1)^3*x(3)^2;

g(4)=x(5)^2.*cos(x(6)*hd)^2-1.076e-4*(31.5+x(5)*x(2)^3*x(4)^2;

g(5)=x(5)*(2*(x(1)+50)*cos(x(6)*hd)+x(1)*x(2)*x(3)-x(2)*x(4)*(31.5+x(5)) ceq=[];

M文件的运行结果为：

********** 两级斜齿轮传动中心距优化设计最优解 ********** 高速级齿轮副模数 Mn1 = 2.0471 mm

- 10 -

《优化设计及有限元分析》上机指导书v2011

低速级齿轮副模数 Mn2 = 3.6059 mm 高速级小齿轮齿数 z1 = 18.5067 低速极小齿轮齿数 z3 = 16.0000 高速级齿轮副传动比 i1 = 5.8000

齿轮副螺旋角 beta = 8.0000 度 减速器总中心距 a12 = 317.4186 mm ————————最优点的性能约束函数值————————

高速机齿轮副接触疲劳强度约束函数值 g1= 0.0000 低速机齿轮副接触疲劳强度约束函数值 g2= 0.0000 高速级大齿轮齿根弯曲强度约束函数值 g3= -1.0052 低速级大齿轮齿根弯曲强度约束函数值 g4= -15.1854

大齿轮齿顶与轴不干涉几何约束函数值 g5= -761.7929

（三） 优化结果处理

经检验，最优点位于性能约束g1(X)、g2(X) 和边界约束 g6(X)、 g12(X)、g14(X)、g16(X) 的交集上。

mn1、mn2高速级和低速级齿轮副模数按照规范圆整为标准值mn1 =2mm和mn2 =4mm；高速级

小齿轮齿数圆整为整数 z1=19；根据高速机传动比 i1=5.8，则高速级大齿轮齿数为z2 =i1z1 =110；根据低速级传动比 i2减速器总中心距：

amn1(z1z2)mn2(z3z4)2cos2(19110)4(1687)2cos8ii131.55.8 =5.413，则高速级大齿轮齿数为 z4i2z=87。 3mm338.292mm

如果将减速器总中心距圆整为 a=340mm，则齿轮副螺旋角调整为 9.838。

四、上机报告要求

上机结束后，学生应完成一份完整的上机报告（手写稿），主要内容包括： 1、 上机目的及上机内容。 2、 优化方法的基本原理。 3、 优化程序代码。 4、 计算结果及分析。

- 11 -