机械优化设计课本中编程实例

燕山大学

机械优化设计论文

专

班

学

姓

业： 12机械工程 级： 工学部1班 号： 名：

2012年 12月 05日

摘 要:

机械优化设计是将最优化原理和计算技术应用于设计领域,为工程设计提供一种重要的科学设计方法。机械优化设计包括建立优化设计问题的数学模型和选择恰当的优化方法与程序两方面的内容。由于机械优化设计是应用数学方法寻求机械设计的最优方案，所以首先要根据实际的机械设计问题建立相应的数学模型，即用数学形式来描述实际设计问题。在建立数学模型时，需要用专业知识确定设计的限制条件和所追求的目标，确立各设计变量之间的相互关系等。机械优化设计问题的数学模型可以是解析式，实验数据或经验公式。虽然它们给出的形式不同，但都是反应设计变量之间的数量关系的。MATLAB是目前国际上最流行的科学与工程计算的软件工具, 它具有强大的数值分析、矩阵运算、信号处理、图形显示、模拟仿真和最优化设计等功能。本文用MATLAB来解决机械设计中的几个常见的问题。

关键词:MATLAB ;优化;机械设计;软件

1 引 言

近年来发展起来的计算机辅助设计，在引入优化设计方法后，使得在设计过程中既能够不断选择设计参数并评选出最优设计方案，又可以加快设计速度，缩短设计周期。在科学技术发展要求机械产品更新日益缩短的今天，把优化设计方法与计算机辅助设计结合起来，使设计过程完全自动化，已成为设计方法的一个重要趋势。

2 采用MATLAB软件进行优化设计

2.1.问题描述:

求ft28t36的最优解

2.1.1规划模型的建立:

目标函数

ft210t36

约束条件

无约束

2.1.2对应的程序:

clc

clear

syms t

f=t^2-10*t+36;

x1=0;

h=2;

f1=subs(f,x1);

x2=x1+h;

f2=subs(f,x2);

f3=f2-1;

t=1;

if (f1-f2)>0

while f3

f3=subs(f,x2+t*h);

t=t+1;

end

x3=x2+(t-1)*h;

else

f3=f2;f2=f1;t=1;f1=f2-1;x3=x2;x2=x1;

while f1

f1=subs(f,x2-t*h);

t=t+1;

end

x1=x2-(t-1)*h;

end

a=x1;

b=x3;

e=1e-05;

k=0.618;

a1=b-k*(b-a);

a2=a+k*(b-a);

f1=subs(f,a1);

f2=subs(f,a2);

c=(b-a)/2;

while c>e

if f1>f2

a=a1;a1=a2;

a2=a+k*(b-a);

f1=f2;

f2=subs(f,a2);

t_min=a2;

f_min=f2;

else

b=a2;a2=a1;

a1=b-k*(b-a);

f2=f1;

f1=subs(f,a1);

t_min=a1;

f_min=f1;

end

c=(b-a)/2;

end

t_min

f_min

2.1.3运行结果

2.1.4 程序基本原理

此程序运用的是黄金分割法，黄金分割法适用于〔a，b〕区间上的任何单谷函数求极小值问题。黄金分割法是建立在区间消去法原理基础上的试探方法，即在搜索区间〔a，b〕内插入两点并计算其函数值，将区间分成三段。应用函数的单谷性质，通过函数值大小的比较，删去其中一段，使搜索区间无限缩小，从而得到极小点的数值近似解。

2.1.5 结 论

黄金分割法可以通过最少的试验次数，找到“最佳点”。它是一个一维的搜索方法，并且比较精确，配合计算机MATLAB软件更能实现求解的快速性和准确性。

2.2.问题描述:

求minf(x)(x112x21)2(49x149x284x12324x2618)2的最优解

2.2.1规划模型的建立:

目标函数

minf(x)(x112x21)2(49x149x284x12324x2618)2 22

初始条件:x01,1 fx03330769 T

约束条件

无约束

2.2.2对应的程序:

syms t1 t2

f=(t1^2+12*t2-1)^2+(49*t1+49*t2+84*t1+2324*t2-681)^2;x0=[1;1]

b=symvar(f);

n=size(b,2);g=cell(n,1);G=cell(n,n);H0=eye(n,n);

for i=1:n

g{i,1}=diff(f,b(i));

end

for i=1:n %产生海塞矩阵

for j=1:n

G{j,i}=diff(g{i},b(j));

end

end

g0=subs(g,b,x0); %产生在x0点的梯度

d0=H0*g0;

a=buchang(x0,d0,f);

x1=x0-a*d0; %产生第一个点

xcha=x1-x0; %计算两点之间的距离

dis=mo(xcha);k=0;

while dis>1e-5 %收敛条件

g1=subs(g,b,x1); %判断在迭代点是否梯度为0

if g1==0

Gy=subs(G,b,x1); %观察海赛矩阵

break

end

y0=g1-g0;

s0=x1-x0;

H1=H0+(s0*s0')/(s0'*y0)-(H0*y0*y0'*H0)/(y0'*H0*y0);

d1=H1*g1;

a=buchang(x1,d1,f);

x2=x1-a*d1;

g0=g1;x0=x1;x1=x2;H0=H1;

xcha=x1-x0;dis=mo(xcha);

x1 %输出每次的迭代点

subs(f,b,x1) %迭代点的函数值

k=k+1 %显示迭代次数

end

2.2.3运行结果

2.1.4 程序基本原理

变量的尺度变换是放大或缩小各个坐标。通过尺度变换可以把函数的偏心程度降低到最低到最低限制。尺度变换技巧能显著地改进几乎所有极小化方法的收敛性质。用最速下降法求f(x1,x2)x125x2的极小值时，需要进行10次迭代才能达到最小值点x*0,0。但是，若变换 T22

y1x1

y25x2

即把x2的尺度放大5倍，就可以将等值线为椭圆的函数f(x1,x2)变换成等值线为圆的函数(y1,y2)y1y2，从而消除了函数的偏心，用最速下降法只需一次迭代即可求得最小值点。

在Newton法中，基本迭代公式

Xk1XktkPk， 22

其中，tk1，

Pk[2f(Xk)]1f(Xk)，

于是有

Xk1XkGk1gk,k0,1,2··· （1）

其中X0是初始点，gk和Gk分别是目标函数f(X)在点Xk的梯度和Hesse矩阵.

为了消除这个迭代公式中的Hesse逆矩阵Gk1，可用某种近似矩阵HkH(Xk)来替换它，即构造一个矩阵序列{Hk}去逼近Hesse逆矩阵序列{Gk1}

此时式（1）变为

Xk1XkHkgk

事实上，式中PkHkgk无非是确定了第k次迭代的搜索方向，为了取得更大的灵活性，我们考虑更一般的的迭代公式

Xk1XktkHkgk (2)

其中步长因子tk通过从Xk出发沿PkHkgk作直线搜索来确定.式（2）是代表很长的一类迭代公式.例如，当HkI（单位矩阵）时，它变为最速下降法的迭代公式.为使Hk确实与Gk1近似并且有容易计算的特点，必须对Hk附加某些条件：

第一，为保证迭代公式具有下降性质，要求{Hk}中的每一个矩阵都是对称 正定的.

理由是，为使搜索方向PkHkgk是下降方向，只要

TTgkPkgkHkgk0

成立即可，即

TgkHkgk0

成立.当Hk对称正定时，此公式必然成立，从而保证式（2）具有下降性质.

第二，要求Hk之间的迭代具有简单形式.显然，

Hk1HkEk (3)

是最简单的形式了.其中Ek称为校正矩阵，式（3）称为校正公式.

第三，必须满足拟Newton条件.即：

Hk1(gk1gk)(Xk1Xk) (4)

为了书写方便也记

ykgk1gk

SkXk1Xk

于是拟Newton条件可写为

Hk1ykSk (5)

有式（3）和（5）知，Ek必须满足

(HkEk)ykSk

或EkykSkHkyk (6)

2.2.5 结 论

通过对本题的求解结果x*0.3363,0.2681 f(x*)5.4318与精确解T

 f(x*)5.9225的比较虽然DFP变尺度法具有综合了梯x*0.28581,0.27936T

度法、牛顿法的优点而又避弃它们各自的缺点，只需计算一阶偏导数，无需计算二阶偏导数及其逆矩阵，对目标函数的初始点选择均无严格要求，收敛速度快的特点,但是还是会有一定的误差。

2.3.问题描述:

求minf(x)(x12)2(x21)2的最优解

2.3.1规划模型的建立:

目标函数

minf(x)(x12)2(x21)2

约束条件

g1x1x20

g1x1x220

初始条件

x03,3 fx05 T22

2.3.2对应的程序:

function y=fun1(x)

y=(x(1)-2).^2+(x(2)-1).^2;

end

function [c,ceq]=gt1(x)

c=x(1).^2-x(2);

ceq=[];

end

A=[1,1];b=2;

x0=[1;2];lb=[];ub=[];

[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun1,x0,A,b,[],[],lb,ub,@gt1);

2.3.3运行结果

2.3.4 程序基本原理

fmincon运行首先检查有无梯度表达提供，如有则选则大规模算法

（subspace trust region），由此涉及到Hessian阵的近似计算，由于已提供了梯度的公式，则Hessian阵可以直接通过有限差分计算。但是如果用户直接提供了Hessian计算公式，则直接计算。 如果没有梯度表达式提供，fmincon选则SQP算法，算法中Hessian阵可以通过BFGS迭代，初始Hessian阵任给。注意BFGS公式中q项是需要计算目标函数梯度得到的。所以Hessian矩阵的近似计算是需要用到有限差分法。

2.3.5 结论

由以上实例可知,采用基于MATLAB软件的计算方法解决优化设计问题,求解程序简单、编程量小、求解效率高,能满足求解离散优化问题的要求,保证得到符合工程要求的最优解,从而可以获得很好的优化结果。

3 总结与体会

利用MATLAB软件可以求解线性规划、无约束规划、二次规划和有约束非线性规划等优化设计问题。该软件编程工作量小，初始参数输入简单，符合工程设计语言，有着很大的优越性。通过给出的优化实例可以看出，应用该软件求解机械优化设计问题非常方便。通过学习机械优化设计让我了解了许多优化算法。真正体会到了科学优化设计的强大与简洁明了,与传统的优化设计方法相比较,大大提高了设计的效率和质量。现代化设计工作是借助电子计算机，应用一些精确

度较高的力学数值分析方法，优化软件进行分析计算，找最优设计方案，实现理论设计代替经验设计。在进行程序求解的过程中，因为是初学Matlab软件，对很多问题的关键点不能够掌握，非线性约束如何书写，函数格式如何书写，变量未定义等等或大或小的问题，但是在一步步排除错误、重新编写程序的过程中，渐渐的对Mtalab熟悉起来，懂得了一些优化方法的简单计算过程和原理，省去了繁琐复杂的优化计算过程。

燕山大学

机械优化设计论文

专

班

学

姓

业： 12机械工程 级： 工学部1班 号： 名：

2012年 12月 05日

摘 要:

机械优化设计是将最优化原理和计算技术应用于设计领域,为工程设计提供一种重要的科学设计方法。机械优化设计包括建立优化设计问题的数学模型和选择恰当的优化方法与程序两方面的内容。由于机械优化设计是应用数学方法寻求机械设计的最优方案，所以首先要根据实际的机械设计问题建立相应的数学模型，即用数学形式来描述实际设计问题。在建立数学模型时，需要用专业知识确定设计的限制条件和所追求的目标，确立各设计变量之间的相互关系等。机械优化设计问题的数学模型可以是解析式，实验数据或经验公式。虽然它们给出的形式不同，但都是反应设计变量之间的数量关系的。MATLAB是目前国际上最流行的科学与工程计算的软件工具, 它具有强大的数值分析、矩阵运算、信号处理、图形显示、模拟仿真和最优化设计等功能。本文用MATLAB来解决机械设计中的几个常见的问题。

关键词:MATLAB ;优化;机械设计;软件

1 引 言

近年来发展起来的计算机辅助设计，在引入优化设计方法后，使得在设计过程中既能够不断选择设计参数并评选出最优设计方案，又可以加快设计速度，缩短设计周期。在科学技术发展要求机械产品更新日益缩短的今天，把优化设计方法与计算机辅助设计结合起来，使设计过程完全自动化，已成为设计方法的一个重要趋势。

2 采用MATLAB软件进行优化设计

2.1.问题描述:

求ft28t36的最优解

2.1.1规划模型的建立:

目标函数

ft210t36

约束条件

无约束

2.1.2对应的程序:

clc

clear

syms t

f=t^2-10*t+36;

x1=0;

h=2;

f1=subs(f,x1);

x2=x1+h;

f2=subs(f,x2);

f3=f2-1;

t=1;

if (f1-f2)>0

while f3

f3=subs(f,x2+t*h);

t=t+1;

end

x3=x2+(t-1)*h;

else

f3=f2;f2=f1;t=1;f1=f2-1;x3=x2;x2=x1;

while f1

f1=subs(f,x2-t*h);

t=t+1;

end

x1=x2-(t-1)*h;

end

a=x1;

b=x3;

e=1e-05;

k=0.618;

a1=b-k*(b-a);

a2=a+k*(b-a);

f1=subs(f,a1);

f2=subs(f,a2);

c=(b-a)/2;

while c>e

if f1>f2

a=a1;a1=a2;

a2=a+k*(b-a);

f1=f2;

f2=subs(f,a2);

t_min=a2;

f_min=f2;

else

b=a2;a2=a1;

a1=b-k*(b-a);

f2=f1;

f1=subs(f,a1);

t_min=a1;

f_min=f1;

end

c=(b-a)/2;

end

t_min

f_min

2.1.3运行结果

2.1.4 程序基本原理

此程序运用的是黄金分割法，黄金分割法适用于〔a，b〕区间上的任何单谷函数求极小值问题。黄金分割法是建立在区间消去法原理基础上的试探方法，即在搜索区间〔a，b〕内插入两点并计算其函数值，将区间分成三段。应用函数的单谷性质，通过函数值大小的比较，删去其中一段，使搜索区间无限缩小，从而得到极小点的数值近似解。

2.1.5 结 论

黄金分割法可以通过最少的试验次数，找到“最佳点”。它是一个一维的搜索方法，并且比较精确，配合计算机MATLAB软件更能实现求解的快速性和准确性。

2.2.问题描述:

求minf(x)(x112x21)2(49x149x284x12324x2618)2的最优解

2.2.1规划模型的建立:

目标函数

minf(x)(x112x21)2(49x149x284x12324x2618)2 22

初始条件:x01,1 fx03330769 T

约束条件

无约束

2.2.2对应的程序:

syms t1 t2

f=(t1^2+12*t2-1)^2+(49*t1+49*t2+84*t1+2324*t2-681)^2;x0=[1;1]

b=symvar(f);

n=size(b,2);g=cell(n,1);G=cell(n,n);H0=eye(n,n);

for i=1:n

g{i,1}=diff(f,b(i));

end

for i=1:n %产生海塞矩阵

for j=1:n

G{j,i}=diff(g{i},b(j));

end

end

g0=subs(g,b,x0); %产生在x0点的梯度

d0=H0*g0;

a=buchang(x0,d0,f);

x1=x0-a*d0; %产生第一个点

xcha=x1-x0; %计算两点之间的距离

dis=mo(xcha);k=0;

while dis>1e-5 %收敛条件

g1=subs(g,b,x1); %判断在迭代点是否梯度为0

if g1==0

Gy=subs(G,b,x1); %观察海赛矩阵

break

end

y0=g1-g0;

s0=x1-x0;

H1=H0+(s0*s0')/(s0'*y0)-(H0*y0*y0'*H0)/(y0'*H0*y0);

d1=H1*g1;

a=buchang(x1,d1,f);

x2=x1-a*d1;

g0=g1;x0=x1;x1=x2;H0=H1;

xcha=x1-x0;dis=mo(xcha);

x1 %输出每次的迭代点

subs(f,b,x1) %迭代点的函数值

k=k+1 %显示迭代次数

end

2.2.3运行结果

2.1.4 程序基本原理

变量的尺度变换是放大或缩小各个坐标。通过尺度变换可以把函数的偏心程度降低到最低到最低限制。尺度变换技巧能显著地改进几乎所有极小化方法的收敛性质。用最速下降法求f(x1,x2)x125x2的极小值时，需要进行10次迭代才能达到最小值点x*0,0。但是，若变换 T22

y1x1

y25x2

即把x2的尺度放大5倍，就可以将等值线为椭圆的函数f(x1,x2)变换成等值线为圆的函数(y1,y2)y1y2，从而消除了函数的偏心，用最速下降法只需一次迭代即可求得最小值点。

在Newton法中，基本迭代公式

Xk1XktkPk， 22

其中，tk1，

Pk[2f(Xk)]1f(Xk)，

于是有

Xk1XkGk1gk,k0,1,2··· （1）

其中X0是初始点，gk和Gk分别是目标函数f(X)在点Xk的梯度和Hesse矩阵.

为了消除这个迭代公式中的Hesse逆矩阵Gk1，可用某种近似矩阵HkH(Xk)来替换它，即构造一个矩阵序列{Hk}去逼近Hesse逆矩阵序列{Gk1}

此时式（1）变为

Xk1XkHkgk

事实上，式中PkHkgk无非是确定了第k次迭代的搜索方向，为了取得更大的灵活性，我们考虑更一般的的迭代公式

Xk1XktkHkgk (2)

其中步长因子tk通过从Xk出发沿PkHkgk作直线搜索来确定.式（2）是代表很长的一类迭代公式.例如，当HkI（单位矩阵）时，它变为最速下降法的迭代公式.为使Hk确实与Gk1近似并且有容易计算的特点，必须对Hk附加某些条件：

第一，为保证迭代公式具有下降性质，要求{Hk}中的每一个矩阵都是对称 正定的.

理由是，为使搜索方向PkHkgk是下降方向，只要

TTgkPkgkHkgk0

成立即可，即

TgkHkgk0

成立.当Hk对称正定时，此公式必然成立，从而保证式（2）具有下降性质.

第二，要求Hk之间的迭代具有简单形式.显然，

Hk1HkEk (3)

是最简单的形式了.其中Ek称为校正矩阵，式（3）称为校正公式.

第三，必须满足拟Newton条件.即：

Hk1(gk1gk)(Xk1Xk) (4)

为了书写方便也记

ykgk1gk

SkXk1Xk

于是拟Newton条件可写为

Hk1ykSk (5)

有式（3）和（5）知，Ek必须满足

(HkEk)ykSk

或EkykSkHkyk (6)

2.2.5 结 论

通过对本题的求解结果x*0.3363,0.2681 f(x*)5.4318与精确解T

 f(x*)5.9225的比较虽然DFP变尺度法具有综合了梯x*0.28581,0.27936T

度法、牛顿法的优点而又避弃它们各自的缺点，只需计算一阶偏导数，无需计算二阶偏导数及其逆矩阵，对目标函数的初始点选择均无严格要求，收敛速度快的特点,但是还是会有一定的误差。

2.3.问题描述:

求minf(x)(x12)2(x21)2的最优解

2.3.1规划模型的建立:

目标函数

minf(x)(x12)2(x21)2

约束条件

g1x1x20

g1x1x220

初始条件

x03,3 fx05 T22

2.3.2对应的程序:

function y=fun1(x)

y=(x(1)-2).^2+(x(2)-1).^2;

end

function [c,ceq]=gt1(x)

c=x(1).^2-x(2);

ceq=[];

end

A=[1,1];b=2;

x0=[1;2];lb=[];ub=[];

[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun1,x0,A,b,[],[],lb,ub,@gt1);

2.3.3运行结果

2.3.4 程序基本原理

fmincon运行首先检查有无梯度表达提供，如有则选则大规模算法

（subspace trust region），由此涉及到Hessian阵的近似计算，由于已提供了梯度的公式，则Hessian阵可以直接通过有限差分计算。但是如果用户直接提供了Hessian计算公式，则直接计算。 如果没有梯度表达式提供，fmincon选则SQP算法，算法中Hessian阵可以通过BFGS迭代，初始Hessian阵任给。注意BFGS公式中q项是需要计算目标函数梯度得到的。所以Hessian矩阵的近似计算是需要用到有限差分法。

2.3.5 结论

由以上实例可知,采用基于MATLAB软件的计算方法解决优化设计问题,求解程序简单、编程量小、求解效率高,能满足求解离散优化问题的要求,保证得到符合工程要求的最优解,从而可以获得很好的优化结果。

3 总结与体会

利用MATLAB软件可以求解线性规划、无约束规划、二次规划和有约束非线性规划等优化设计问题。该软件编程工作量小，初始参数输入简单，符合工程设计语言，有着很大的优越性。通过给出的优化实例可以看出，应用该软件求解机械优化设计问题非常方便。通过学习机械优化设计让我了解了许多优化算法。真正体会到了科学优化设计的强大与简洁明了,与传统的优化设计方法相比较,大大提高了设计的效率和质量。现代化设计工作是借助电子计算机，应用一些精确

度较高的力学数值分析方法，优化软件进行分析计算，找最优设计方案，实现理论设计代替经验设计。在进行程序求解的过程中，因为是初学Matlab软件，对很多问题的关键点不能够掌握，非线性约束如何书写，函数格式如何书写，变量未定义等等或大或小的问题，但是在一步步排除错误、重新编写程序的过程中，渐渐的对Mtalab熟悉起来，懂得了一些优化方法的简单计算过程和原理，省去了繁琐复杂的优化计算过程。