零件的参数设计摘要:本题目对零件的参数这一问题,综合考虑重新设计零件的参数(包括标定值和容差),并与原设计进行比较,得出最优化的数学模型,并对模型进行求解,最后用计算机模拟对模型的最优解进行检验。
由题意知粒子分离器的参数y 由零件参数1234567,,,,,,x x x x x x x 的参数决定,参数i x 的容差等级决定了产品的成本,y 偏离0y 的值决定了产品的损失,问题就是寻找零件的最优标定值和最优等级搭配,使得批量生产时的总费用最少。
一、 问题的重述:一件产品由若干零件组装而成,标志产品性能的某个参数取决于这些零件的参数。
零件参数包括标定值和容差两部分。
进行成批生产时,标定值表示一批零件该参数的平均值,容差则给出了参数偏离其标定值的容许范围。
若将零件参数视为随机变量,则标定值代表期望值,在生产部门无特殊要求时,容差通常规定为均方差的3倍。
进行零件参数设计,就是要确定其标定值和容差。
这时要考虑两方面因素:一是当各零件组装成产品时,如果产品参数偏离预先设定的目标值,就会造成质量损失,偏离越大,损失越大;二是零件容差的大小决定了其制造成本,容差设计得越小,成本越高。
试通过如下的具体问题给出一般的零件参数设计方法。
粒子分离器某参数(记作y )由7个零件的参数(记作x 1,x 2,...,x 7)决定,经验公式为:7616.1242356.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元。
零件参数的标定值有一定的容许范围;容差分为A、B、C三个等级,用与标定值的相对值表示,A等为±1%,B等为±5%,C等为±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 所造成的损失和零件成本,先列出总费用的数学模型表达如下:72311000*(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 ,则()()71i ii iff x y f x x x =∂==+∆∂∑将标定值xi (i=1,2,3,…..7)带入经验公式得()i y f x =得71ii ify y y x x =∂∆=-=∆∂∑ 由于在加工零件时,在标定值知道的情况下,加工误差服从正太分布,即()2~0,i x N σ∆ 且∆xi 相互独立,由正态分布性质可知()2~0,y y N σ∆,()2~,y y N y σ ,由误差传递公式得22277211i yi i i i i i i f f x x x x σσσ==⎛⎫⎛⎫⎛⎫∂∂== ⎪ ⎪ ⎪∂∂⎝⎭⎝⎭⎝⎭∑∑由于容差均为方差的3倍,容差与标定值的比值为容差等级,则30.01,0.05,0.1i i x σ⎛⎫= ⎪⎝⎭, y 的分布密度函数为()()221y 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.82 1.21.6P y d y y d y ψψ=+⎰⎰y 偏离00.3y ±的概率,即废品的概率为()()()()1.23 1.8P 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 造成的损失和零件成本,设计最优零件参数的模型建立如下目标函数7237min 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=resultfval=inf;ticB(1)=2;B(5)=3;for b=2:3B(2)=b;for c=1:3B(3)=c;for d=1:3B(4)=d;for f=1:3B(6)=f;for g=1:2B(7)=g;[fv,x]=getcost(B);if fv<fvalXmin=x;Bmin=B;fval=fv;end;end;end;end;end;end;f=fval,Xmin,Bmin,p=getP(Xmin,Bmin)tocsimulation(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;elsef=f+20;end;switch (B(3))case 1f=f+200;case 2f=f+50;case 3f=f+20;end;switch (B(4))case 1f=f+500;case 2;f=f+100;case 3f=f+50;end;switch (BC(6))case 1f=f+100;case 2f=f+25;case 3f=f+10;end;if(B(7)==1)f=f+100elsef=f+25;end;f=f+p(2)*1000+p(3)*9000;function f=getfcY(MU,B)f=0;B=int32(B);for i=1:7if B(i)==1sigma(i)=MU(i)*0.01/3;end;if B(i)==2sigma(i)=MU(i)*0.01/3;end;if B(i)==3sigma(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:itera = abs(Yfun(getparaX(MU,B))-1.5);if a<0.1fF(1)=f(1)+1;end;if a<0.3&a>=0.1f(2)=f(2)+1;end;if a>=0.3f(3)=f(3)+1end;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:7if B(i) = = 1sigma0(i) = MU(i)*0.01;end;if B(i) = = 2sigma0(i) = MU(i)*0.05;end;if B(i) = = 3sigma0(i) = MU(i)*0.1;end;f(i)=normrnd(MU(I),sigma(i)/3);end;。