数模-零件的参数设计

零件的参数设计

摘要：

本题目对零件的参数这一问题，综合考虑重新设计零件的参数（包括标定值和容差），并与原设计进行比较，得出最优化的数学模型，并对模型进行求解，最后用计算机模拟对模型的最优解进行检验。由题意知粒子分离器的参数y 由零件参数1234567,,,,,,x x x x x x x 的参数决定，参数i x 的容差等级决定了产品的成本，y 偏离0y 的值决定了产品的损失，问题就是寻找零件的最优标定值和最优等级搭配，使得批量生产时的总费用最少。

一、 问题的重述：

一件产品由若干零件组装而成，标志产品性能的某个参数取决于这些零件的参数。零件参数包括

标定值和容差两部分。进行成批生产时，标定值表示一批零件该参数的平均值，容差则给出了参数偏离其标定值的容许范围。若将零件参数视为随机变量，则标定值代表期望值，在生产部门无特殊要求时，容差通常规定为均方差的３倍。

进行零件参数设计，就是要确定其标定值和容差。这时要考虑两方面因素：一是当各零件组装成产品时，如果产品参数偏离预先设定的目标值，就会造成质量损失，偏离越大，损失越大；二是零件容差的大小决定了其制造成本，容差设计得越小，成本越高。 试通过如下的具体问题给出一般的零件参数设计方法。

粒子分离器某参数（记作y ）由7个零件的参数（记作x 1,x 2,...,x 7）决定，经验公式为：

7616

.1242

3

56

.02485.01235136.0162.2142.174x x x x x x x x x x x Y ???

? ?????

?

?????

???

??--????

? ??-????? ???=-

y 的目标值（记作0y ）为1.50。当y 偏离0y ±0.1时，产品为次品，质量损失为1,000元；当y 偏离0y ±0.3时，产品为废品，损失为9,000元。

零件参数的标定值有一定的容许范围；容差分为Ａ、Ｂ、Ｃ三个等级，用与标定值的相对值表示，Ａ等为±1%，Ｂ等为±5%，Ｃ等为±10%。7个零件参数标定值的容许范围，及不同容差等级零件的成本（元）如下表（符号／表示无此等级零件）：

现进行成批生产，每批产量1,000个。在原设计中，7个零件参数的标定值为：

x 1=0.1,x 2=0.3,x 3=0.1,x 4=0.1,x 5=1.5,x 6=16,x 7=0.75;容差均取最便宜的等级。

请你综合考虑y 偏离0y 造成的损失和零件成本，重新设计零件参数（包括标定值和容差），并与原设计比较，总费用降低了多少？

二、问题的假设

1、假设在加工零件时，在确定了标定值的情况下，零件的误差服从正太分布且各个零件的误差是

相互独立的。

2、假设制造零件的总费用只由零件的损失费用和成本组成，不必考虑其他外在因素。

3、假设题目所给的经验公式足够反映参数1234567,,,,,,x x x x x x x ，对参数y 的影响，而且经验公式有足够高的精度，即不考虑经验公式的误差。

三、符号说明

四、模型的建立

由题意可以知道，容差如果变大，则生产产品的的成本会降低，但同时y 偏离0y 的程度也增大，从而导致了损失的增加，由此我们要求出一个最优解，使得总费用最低。为了确定原设计中标定值（xi （i=1，2，3，….，7）的期望值）及已给的容差对产品性能参数影响而导致的总损失w ，即确定y 偏离目标值0y 所造成的损失和零件成本，先列出总费用的数学模型表达如下：

7

231

1000*(10009000)ij i W C P P ==++∑

为了确定总损失w ，必须知道123,,P P P （即正品、次品及废品的概率）。为此，用泰勒公式将经验公式在X=i x （i=1，2，3，…..7）处展开并略去高次项（原因：误差本来就在0.01级别，它的高阶无穷小完全可以忽略），后来研究y 的概率分布，设f （x ）=y ，

则

()()7

1i i

i i

f

f x y f x x x =?==+??∑

将标定值xi （i=1，2，3，…..7）带入经验公式得()i y f x =

得

7

1i

i i

f

y y y x x =??=-=??∑ 由于在加工零件时，在标定值知道的情况下，加工误差服从正太分布，即()2~0,i x N σ? 且?xi 相互独立，由正态分布性质可知

()2~0,y y N σ?，()

2

~,y y N y σ ，

由误差传递公式得

2

2

2

7

7

2

11i y

i i i i i i i f f x x x x σσσ==??????

??== ? ? ?????????

∑∑

由于容差均为方差的3倍，容差与标定值的比值为容差等级，则30.01,0.05,0.1i i x σ??

= ???

， y 的分布密度函数为

(

)()2

2

1y y y y e

σ

ψ--

=

产品为正品时y 的范围是[1.2 ]1.6

产品为次品时y 的范围是[1.2 ]1.4和[1.6 ]1.8， 产品为废品时y 的范围是(-∞ ]1.2和[1.8 )+∞

y 偏离00.1y ±的概率，即次品的概率为

()()()()1.4 1.8

2 1.2

1.6

P y d y y d y ψψ=+??

y 偏离00.3y ±的概率，即废品的概率为

()()()()1.23 1.8

P y d y y d y ψψ+∞

-∞

=+??

由于y 偏离0y 越远，损失越大，所以在y σ固定时，调整y 使之等于目标值0y 可降低损失。 取0y y y ?=- 即y =0y ， 则

20.1y P σ??

=Φ ? ?

??

，30.3y P σ??=Φ ? ??? ()t Φ为标准正太分布函数。综合考虑y 偏离0y 造成的损失和零件成本，设计最优零件参数的模型建立

如下目标函数

7

237

min 1000*(10009000)ij i W C P P ==++∑

五、模型的求解

初步分析，对于原给定的方案，利用matlab 编程计算（见附录），计算结果如下

由于按原设计方案设计的产品频率过低，损失费过高，显然设计不合理。进一步分析发现，参数均值

y =1.7256偏离目标值0y =1.5太远，致使损失过大。尽管原设计方案保证了成本最低，但由于零件参

数的精度过低，导致正品率也过低，损失较高。所以我们应综合考虑成本费和损失费。

模型的实现过程：本模型通过matlab 进行求解，我们通过理论模型求解和随机模拟的求解过程如下：在给定容差等级的情况下，利用matlab 中求解非线性规划的函数fmincon ，通过多次迭代求解，最终球的一组最优解。最初，我们设定的fmincon 函数目标函数就是总费用，约束条件为各个标定值的容许范围，以及各零件标定值带入产品参数表达式应为0y ，即1.5.然而，在迭代过程中我们发现，求解过程十分慢，因此，我们在仔细对matlab 实现代码进行研究发现，求解过程非常慢，为了提高速度，我们首先利用matlab 的diff 函数对产品参数中的各个表达式进行求偏导，然后得到多个带参表达式，利用int 函数对y 的概率密度函数进行积分，分别得到出现次品和废品的概率的表达式，然后将这些表达式写进程序里，这样在求解过程中就不需要在每一次迭代中都要求偏导和积分了，修改后的程序运行时间大大减少。

六、模型检验

对设计方案进行模型检验模拟，由于每种零件参数均服从正态分布，用正态分布随机发生器在每种零件参数允许的范围内产生1000个随机数参与真实值i x 的计算随机模拟多次后结果如下：

七、误差分析

1、在建模过程中，通过泰勒公式将()y f x =展开并略去高次项使线性化，不可避免地产生可截断误差，所以展开后的式子致使原经验公式的近似关系式。但在一般情况下，线性化和在求和在实用上具有足够的精度，所以由于函数线性化而略去的高次项可以忽略不计。在函数关系式叫复杂的情况下，将其线性化更具有明显的优势。

2、本模型忽略了小概率事件的发生的可能，认为零件的参数只可能出现在允许范围内，即

[]3,3i i i i x x σσ-+，现实中，小概率事件仍有能发生，但是在大批量生产中，小概率事件发生对最终结

果没有影响，所以可以忽略。

3、该模型对于质量损失的计算，将所有函数都看做连续函数，而这对于每个零件而言是不可能的，所以其中也会产生误差。

八、模型优缺点

优点：

1、 建模过程中，采用泰勒公式将经验公式简化，并假设各零件参数都满足大量数据的正态分布，使

得整个模型的建立及求解得到大大简化。

2、 本模型运用概率统计与优化知识对零件参数进行优化设计。通过建立一个反应设计要求的数学模

型，利用matlab 软件，经过编程来实现对设计方案参数的调整，将总费用由3074.8（元/个），结果还是令人十分满意的。 缺点：

1、 本模型在模型的求解过程中，对一些可接受范围内的误差直接进行了忽略，因而对于结果的精确

性还是会有影响。

2、 本模型时建立在一些假设中的，所有实用性受到了限制，在实际生产中，如果可以把更多的一些

因素考虑进去会更好。在已假设的条件下，本模型的优化结果还是好的。

附录：

function f=result

fval=inf;

tic

B(1)=2;

B(5)=3;

for b=2:3

B(2)=b;

for c=1:3

B(3)=c;

for d=1:3

B(4)=d;

for f=1:3

B(6)=f;

for g=1:2

B(7)=g;

[fv,x]=getcost(B);

if fv

Xmin=x;

Bmin=B;

fval=fv;

end;

end;

end;

end;

end;

end;

f=fval,Xmin,Bmin,p=getP(Xmin,Bmin)

toc

simulation(Xmin,Bmin);

function [f,x]=getcost(B)

MU=[0.1 0.3 0.1 0.1 1.5 16 0.75];

options=optimset('largescale','off');

[x,fval]=fmincon('getfcY',MU,[],[],[],[],[],[],'mycon',options,B); x,B,f=cost(x,B)

function [c,ceq]=mycon(MU,B)

c(1)=MU(1)-0.125;

c(2)0.075-MU(1);

c(3)=MU(2)-0.375;

c(4)=0.225-MU(2);

c(5)=MU(3)-0.125;

c(6)=0.075-MU(3);

c(7)=MU(4)-0.125;

c(8)=0.075-MU(4);

c(9)=MU(5)-1.875;

c(10)=1.125-MU(5);

c(11)=MU(6)-20;

c(12)=12-MU(6);

c(13)=MU(7)-0.935;

c(14)=0.5625-MU(7);

ceq(1)=Yfun(MU)-1.5;

function f=cost(MU,B)

f=25;

p=getP(MU,B;

if (B(2)=2)

f=f+50;

else

f=f+20;

end;

switch (B(3))

case 1

f=f+200;

case 2

f=f+50;

case 3

f=f+20;

end;

switch (B(4))

case 1

f=f+500;

case 2;

f=f+100;

case 3

f=f+50;

end;

switch (BC(6))

case 1

f=f+100;

case 2

f=f+25;

case 3

f=f+10;

end;

if(B(7)==1)

f=f+100

else

f=f+25;

end;

f=f+p(2)*1000+p(3)*9000;

function f=getfcY(MU,B)

f=0;

B=int32(B);

for i=1:7

if B(i)==1

sigma(i)=MU(i)*0.01/3;

end;

if B(i)==2

sigma(i)=MU(i)*0.01/3;

end;

if B(i)==3

sigma(i)=MU(i)*0.05/3;

end;

end;

x1=MU(1);x2=MU(2);x3=MU(3);x4=MU(4);x5=MU(5);x6=MU(6);x7=MU(7); f=(pd1(x1,x2,x3,x4,x5,x6,x7)*sigma(1))^2;

f=(pd2(x1,x2,x3,x4,x5,x6,x7)*sigma(2))^2;

f=(pd3(x1,x2,x3,x4,x5,x6,x7)*sigma(3))^2;

f=(pd4(x1,x2,x3,x4,x5,x6,x7)*sigma(4))^2;

f=(pd5(x1,x2,x3,x4,x5,x6,x7)*sigma(5))^2;

f=(pd6(x1,x2,x3,x4,x5,x6,x7)*sigma(6))^2;

f=(pd7(x1,x2,x3,x4,x5,x6,x7)*sigma(7))^2;

f=abs(f^0.5);

function f=pd1(x1,x2,x3,x4,x5,x6,x7)

f=8721/50/x(5)*(x(3)/(x(2)-x(1)))^(17/20)*((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*(x(4)/x(2))^(29/25))/x (6)/x(7))^(1/2)+148257/1000*x(1)/x(5)/(x(3)/(x(2)-x(1)))^(3/20)*((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)* (x(4)/x(2))^(29/25))/x(6)/x(7))^(1/2)*x(3)/(x(2)-x(1))^2;

function f=pd2(x1,x2,x3,x4,x5,x6,x7)

f=-148257/1000*x(1)/x(5)/(x(3)/(x(2)-x(1)))^...

(3/20)*((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))...

^(3/2)*(x(4)/x(2))^(29/25))/x(6)/x(7))^(1/2)*x(3)/(x(2)...

-x(1))^2+8721/100*x(1)/x(5)*(x(3)/(x(2)-x(1)))^(17/20)/...

((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*(x(4)/x(2))^...

(29/25))/x(6)/x(7))^(1/2)*(24759/31250*(1-9/25/(x(4)/x(2))^(14/25))^(1/2)/(x(4)/x(2))^...

(2/5)*x(4)/x(2)^2+3799/1250*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*(x(4)/x(2))^(4/25)*x(4)/x(2)^2)/x(6)/x(7);

function f=pd3(x1,x2,x3,x4,x5,x6,x7)

f=148257/1000*x(1)/x(5)/(x(3)/(x(2)-x(1)))^(3/20)*((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*(x(4)/x(2))^(2 9/25))/x(6)/x(7))^(1/2)/(x(2)-x(1));

function f=pd4(x1,x2,x3,x4,x5,x6,x7)

f=8721/100*x(1)/x(5)*(x(3)/(x(2)-x(1)))^(17/20)/...

((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*...

(x(4)/x(2))^(29/25))/x(6)/x(7))^(1/2)*(-24759/31250*...

(1-9/25/(x(4)/x(2))^(14/25))^(1/2)/(x(4)/x(2))^(2/5)/x(2)-...

3799/1250*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*(x(4)/x(2))^(4/25)/x(2))/x(6)/x(7);

function f=pd5(x1,x2,x3,x4,x5,x6,x7)

f=-8721/50*x(1)/x(5)^2*(x(3)/(x(2)-x(1)))^(17/20)*((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*(x(4)/x(2))^(2 9/25))/x(6)/x(7))^(1/2);

function f=pd6(x1,x2,x3,x4,x5,x6,x7)

f=d(6)=-8721/100*x(1)/x(5)*(x(3)/(x(2)-x(1)))^...

(17/20)/((1-131/50*(1-9/25/(x(4)/x(2))^(14/25))...

^(3/2)*(x(4)/x(2))^(29/25))/x(6)/x(7))^(1/2)*...

(1-131/50*(1-9/25/(x(4)/x(2))^(14/25))^(3/2)*(x(4)/x(2))^(29/25))/x(6)^2/x(7);

function f=pd7(x1,x2,x3,x4,x5,x6,x7)

f=-8721/100*x(1)/x(5)*(x(3)/(x(2)-x(1)))...

^(17/20)/((1-131/50*(1-9/25/(x(4)/x(2))^...

(14/25))^(3/2)*(x(4)/x(2))^(29/25))/...

x(6)/x(7))^(1/2)*(1-131/50*(1-9/25/(x(4)/x(2))...

^(14/25))^(3/2)*(x(4)/x(2))^(29/25))/x(6)/x(7)^2;

function f=getP(MU,B)

yb=Yfun(MU);

fc=getfcY(MU,B);

f(2)=jf1(yb,fc);

f(3)=jf2(yb,fc);

f(1)=1-f(2)-f(3);

f=double(f);

function f=jf1(u,a0)

f=--1125899906842624/5644425081792261*...

erf(1/10*2^(1/2)*(-9+5*u)/a0)*2^(1/2)*pi*(1/2)...

+1125899906842624/5644425081792261*erf(1/10)*2^(1/2)*...

(-8+5*u)/a0)*2^(1/2)*pi^(1/2)-1125899906842624/5644425081792261*...

erf(1/10*2^(1/2)*(-7+5*u)/a0)*2^(1/2)*pi*(1/2)+1125899906842624/5644425081792261*...

erf(1/10*2^(1/2)*(-6+5*u)/a0)*2^(1/2)*pi^(1/2);

function f = jf2(u,a0)

f = -1125899906842624/5644425081792261*erf...

(1/2*2^(1/2)*(-10+u)/a0)*2^(1/2)*pi^(1/2)+1125899906842624/...

5644425081792261*erf(1/10*2^(1/2)*(-9+5*u)/a0)*2^(1/2)*pi^(1/2)-...

1125899906842624/5644425081792261*erf(1/10*2^(1/2)*(-6+5*u)/a0)*2^...

(1/2)*pi^(1/2) + 1125899906842624/5644425081792261*erf(1/2*2^(1/2)*(10+u)/a0)*2^(1/2)*pi^(1/2)

function f=Yfun(x)

f=174.42*(x(1)/x(5))*(x(3)/(x(2)-x(1)))^0.85*sqrt((1-2.62*(1-0.36*(x(4)/...

x(2))^(-0.56))^1.5*(x(4)/x(2))^1.16)/(x(6)*x(7)));

function f = geteveryP(MU,B,iter)

f(1)=0;

f(2)=0;

f(3)=0;

for i = 1:iter

a = abs(Yfun(getparaX(MU,B))-1.5);

if a<0.1

fF(1)=f(1)+1;

end;

if a<0.3&a>=0.1

f(2)=f(2)+1;

end;

if a>=0.3

f(3)=f(3)+1

end;

end;

f(1)=f(1)/iter;

f(2)=f(2)/iter;

f(3)=f(3)/iter;

function f=getparaX(MU,B)

B=int32(B);

for i = 1:7

if B(i) = = 1

sigma0(i) = MU(i)*0.01;

end;

if B(i) = = 2

sigma0(i) = MU(i)*0.05;

end;

if B(i) = = 3

sigma0(i) = MU(i)*0.1;

end;

f(i)=normrnd(MU(I),sigma(i)/3);

end;