数模-零件的参数设计

零件的参数设计摘要：本题目对零件的参数这一问题，综合考虑重新设计零件的参数（包括标定值和容差），并与原设计进行比较，得出最优化的数学模型，并对模型进行求解，最后用计算机模拟对模型的最优解进行检验。

由题意知粒子分离器的参数y 由零件参数1234567,,,,,,x x x x x x x 的参数决定，参数i x 的容差等级决定了产品的成本，y 偏离0y 的值决定了产品的损失，问题就是寻找零件的最优标定值和最优等级搭配，使得批量生产时的总费用最少。

一、 问题的重述：一件产品由若干零件组装而成，标志产品性能的某个参数取决于这些零件的参数。

零件参数包括标定值和容差两部分。

进行成批生产时，标定值表示一批零件该参数的平均值，容差则给出了参数偏离其标定值的容许范围。

若将零件参数视为随机变量，则标定值代表期望值，在生产部门无特殊要求时，容差通常规定为均方差的３倍。

进行零件参数设计，就是要确定其标定值和容差。

这时要考虑两方面因素：一是当各零件组装成产品时，如果产品参数偏离预先设定的目标值，就会造成质量损失，偏离越大，损失越大；二是零件容差的大小决定了其制造成本，容差设计得越小，成本越高。

试通过如下的具体问题给出一般的零件参数设计方法。

粒子分离器某参数（记作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元。

零件参数的标定值有一定的容许范围；容差分为Ａ、Ｂ、Ｃ三个等级，用与标定值的相对值表示，Ａ等为±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 所造成的损失和零件成本，先列出总费用的数学模型表达如下：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;。