数学建模 零件参数的优化设计 ()

合集下载
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

零件参数的
优化设计
摘 要
本文建立了一个非线性多变量优化模型。

已知粒子分离器的参数y 由零件参数)72,1( =i x i 决定,参数i x 的容差等级决定了产品的成本。

总费用就包括y 偏离y 0造成的损失和零件成本。

问题是要寻找零件的标定值和容差等级的最佳搭配,使得批量生产中总费用最小。

我们将问题的解决分成了两个步骤:1.预先给定容差等级组合,在确定容差等级的情况下,寻找最佳标定值。

2.采用穷举法遍历所有容差等级组合,寻找最佳组合,使得在某个标定值下,总费用最小。

在第二步中,由于容差等级组合固定为108种,所以只要在第一步的基础上,遍历所有容差等级组合即可。

但是,这就要求,在第一步的求解中,需要一个最佳的模型使得求解效率尽可能的要高,只有这样才能尽量节省计算时间。

经过对模型以及matlab 代码的综合优化,最终程序运行时间仅为3.995秒。

最终计算出的各个零件的标定值为:
i x ={0.0750,0.3750,0.1250,0.1200,1.2919,15.9904,0.5625},
等级为:B B C C B B B d ,,,,,,=
一台粒子分离器的总费用为:421.7878元
与原结果相比较,总费用由3074.8(元/个)降低到421.7878(元/个),降幅为86.28%,结果是令人满意的。

为了检验结果的正确性,我们用计算机产生随机数的方式对模型的最优解进行模拟检验,模拟结果与模型求解的结果基本吻合。

最后,我们还对模型进行了误差分析,给出了改进方向,使得模型更容易推广。

关键字:零件参数 非线性规划 期望 方差
一、问题重述
一件产品由若干零件组装而成,标志产品性能的某个参数取决于这些零件的参数。

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

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

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

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

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

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

粒子分离器某参数(记作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 的目标值(记作y 0)为1.50。

当y 偏离y 0+0.1时,产品为次品,质量损失为1,000元;当y 偏离y 0+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 偏离y 0造成的损失和零件成本,重新设计零件参数(包括标定值和容差),并与原设计比较,总费用降低了多少?
二、模型假设
1、将各零件参数视为随机变量,且各自服从正态分布;
2、假设组成离子分离器的各零件互不影响,即各零件参数互相独立;
3、假设小概率事件不可能发生,即认为各零件参数只可能出现在容许范围内;
4、在大批量生产过程中,整批零件都处于同一等级,。

本题可认为1000各零件都为A 等、B 等或C 等;
5、生产过程中出质量损失外无其他形式的损失;
6、在质量损失计算过程中,认为所有函数都是连续可导的。

三、符号说明
i x :第i 类零件参数的标定值(i=1,2……7);
i x ∆:第i 类零件参数的实际值相对目标值的偏差(i=1,2……7);
i r :第i 类零件参数的容差(i=1,2,……7);
i σ:第i 类零件参数的方差(i=1,2,……7);
i i b a ,:标定值i x 的上、下限;
y :离子分离器某参数的实际值;
0y :离子分离器该参数的目标值;
y :离子分离器某参数的均值;
y ∆:离子分离器某参数的实际值y 相对平均值y 的偏差;
y σ:离子分离器某参数的方差;
1P :一批产品中正品的概率;
2P :一批产品中次品的概率;
3P :一批产品中废品的概率;
W :一批产品的总费用(包括损失和成本费)

ij C :第i 类零件对应容差等级为j 的成本(j=A,B,C ) 单位:元/个。

四、问题分析
该问题是一定约束条件下的最优化问题,经分析题意,拟建立以总费用为目标函数的非线性规划模型。

总费用由损失费和成本费两部分组成,零件成本由简单的线性代数式决定,而损失费涉及概率分布的非线性函数。

要求出损失费,就必须知道一批产品的次品率和废品率,结合各类零件都服从),(2
i i x N ,可假设y 也服从正态分布,联想正态分布的性质——当各变量均服从正态分布时,其线性组合也服从正态分布。

题中所给经验公式为一复杂的非线性的公式,无法直接对其分析处理,所以需借助泰勒公式将其展开并作相应处理使其线性化。

而对于零件成本,需先确定容差等级才能求得成本费。

由容差等级和各类零件的标定值
i x 便可知道给类零件的容差i r 。

最后,便将问题转化为i x 、i r 关于总目标函数的
最优解的问题上。

在进行零件参数设计时,如果零件设计不妥,造成产品参数偏离预先设定值,
就会造成质量损失,且偏差越大,损失也越大;零件容差的大小决定了其制造成本,容差设计得越小(即精度越高)零件成本越高。

合理的设计方案应既省费用又能满足产品的预先设定值,设计方向应该如下:
(1)设计的零件参数,要保证由零件组装成的产品参数符合该产品的预先设定值,即使有偏离也应是在满足设计最优下的容许范围。

(2)零件参数(包括标定值和容差等级)的设计应使总费用最小为优。

此外分析零件的成本及产品的质量损失不难发现,质量损失对费用的影响远大于零件成本对费用的影响,因而设计零件参数时,主要考虑提高产品质量来达到减少费用的目的。

五、模型建立
为了确定原设计中标定值(x i i (,,,)=127 的期望值)及已给的容差对产品性能参数影响而导致的总损失W ,即确定y 偏离目标值y 0所造成的损失和零件成本,先列出总损失的数学模型表达如下: )90001000(1000327
1P P C W i ij ++⨯=∑=
当然,为了确定总损失W ,必须知道1P 、2P 、3P (即正品、次品及废品的概率)。

为此,将经验公式用泰勒公式在)72,1( ==i x X i 处展开并略去二次以上高次项后来研究y 的概率分布,设y x f =)(,则

=∆∂∂+==7
1)()(i i i
i x x f
x f y X f 将标定值)72,1( =i x i 带入经验公式即得 )(i x f y = 所以 i i i
x x f
y y y ∆∂∂=-=∆∑
=7
1 由于在加工零件时,在标定值知道的情况下,加工误差服从正态分布,即 )(2
,0N ~i i x σ∆ 且i x ∆相互独立,由正态分布性质可知
),0(~2
y N y σ∆ ),(~2
y y N y σ
由误差传递公式得 227
12
7
12
)()()(i i i i i
i i i y
x x x f x f σσσ∑∑==∂∂=∂∂= (1)
由于容差为均方差的3倍,容差与标定值的比值为容差等级,则
3
1
.0,305.0,301.0=
i
i
x σ y 的分布密度函数为
()
2
2
21)(y y y e
y y
σσπϕ-=
-
y 偏离1.00±y 的概率,即次品的概率为
⎰⎰+=8
.16
.14
.12
.12)()()()(y d y y d y P ϕϕ (2)
y 偏离3.00±y 的概率,即废品的概率为
⎰⎰+∞

-+=8
.12
.13)()()()(y d y y d y P ϕϕ (3)
由于y 偏离0y 越远,损失越大,所以在y σ固定时,调整y 使之等于目标值0y 可降低损失。

取0y y y -=∆即0y y =,则 )1
.0(
2y
P σφ= )3
.0(
3y
P σφ=
)(t φ为标准正态分布函数。

综合考虑y 偏离y 0造成的损失和零件成本,设计最优零件参数的模型建立如下:
目标函数
min )90001000(1000327
1P P C W i ij ++⨯=∑=
s.t. )72,1( =≤≤i a x b i
i i )72,1()(0 ==i x f y i
六、模型求解
初略分析
对于原给定的设计方案,利用matlab编程计算(见附录),计算结果如下:
由于按原设计方案设计的产品正品率过低,损失费过高,显然设计不够合理。

进一步分析发现,参数均值y=1.7256偏离目标值
y=1.5太远,致使损失过大。

尽管原设计方案保证了正本最低,但由于零件参数的精度过低,导致正品率也过低。

所以我们应综合考虑成本费和损失费。

模型的实现过程:
本模型通过matlab进行求解,我们通过理论模型求解和随机模拟的求解过程如下:在给定容差等级的情况下,利用matlab中求解非线性规划的函数fmincon,通过多次迭代求解,最终求得一组最优解。

最初,我们设定的fmincon 函数的目标函数就是总费用,约束条件为各个标定值的容许范围,以及各零件标
y,即1.5。

然而,在迭代过程中我们发现,求解定值带入产品参数表达式应为
过程十分慢,在给定容差等级的确定的情况下,计算最优标定值需要将近400秒,如果在此基础上对108种容错等级进行穷举查找最优组合,将需要大概12小时。

显然这是不合理的。

因此,我们在仔细对matlab实现代码研究发现,求解过程之所以慢,是因为代码中存在多次调用求偏导和积分的函数,在fmincon 的多次迭代中,耗费大量时间。

所以,为了提高求解速度,我们首先利用matlab 中diff函数对产品参数中的各个表达式进行求偏导,然后得到多个带参表达式,利用int函数对y的概率密度函数进行积分,分别得到出现次品和废品概率的表达式,然后将这些表达式写进程序里,这样在求解过程中就不需要在每一次迭代中都要求偏导和积分了,修改后的程序运行时间大大减少。

程序流程图
离子分离器参数均值y =1.5 离子分离器参数方差y σ=0.0689
模型检验
对设计方案进行动态模拟,由于每种零件参数均服从正态分布,用正态分布
x的计随机数发生器在每种零件参数允许范围内产生1000个随机数参与真实值
i
算随机模拟N次后结果如下:
正品率次品率废品率成本费损失费总费用0.8570 0.1430 0.0000 275 143 418
σ=0.0689画出y的概率分布图,再对x随机取样画根据最优解的y=1.5,
y
出y的概率分布图(见图6.1),由图可知:两组数据所画概率分布图的拟合度相当高,进一步确保了模型的正确性。

图6.1概率分布图对比图
通过以上数据,与原设计方案所得结果相比较,总费用由3074.8(元/个)降低到421.7878(元/个),降幅为86.28%,结果是令人满意的。

七、误差分析
1、在建模过程中,通过泰勒公式将)
y=展开并略去二次及以上项使线
(X
f
性化,不可避免地产生了截断误差,所以展开后的式子只是原经验公式的近似关
系式。

但在一般情况下,线性化和求总和在实用上具有足够的精度,所以由于函数线性化而略去的高次项可以忽略不计。

在函数关系式较复杂的情况下,将其线性化更具有明显的优势。

2、本模型忽略了小概率事件发生的可能,认为零件的参数只可能出现在允 范围内,即[]i i i i x x σσ3,3+-。

现实中,小概率事件仍有发生的可能性,但在大批量生产中,小概率事件的发生对最终结果没有影响,所以可以忽略。

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

八、模型的评价及推广
1.优点
(1)建模过程中,采用泰勒公式将经验公式简化,并假设各零件参数都服从满足大量数据的正态分布,使得整个模型的建立及求解得到大大简化。

(2)本模型运用概率统计与优化知识对零件参数进行优化设计。

通过建立一个反映设计要求的数学模型,利用MATLAB 软件,经过编程来实现对设计方案参数的调整,将总费用由3074.8(元/个)降低到421.7878(元/个),降幅达到86.28%,结果还是令人十分满意的。

(3)本模型在程序运算的过程中,做了适当处理,将每次循环本该由计算机求偏导和积分的提前人为处理,将求偏导和积分后的算式写入程序中,这样大大节约了运算时间,将运行时间由几个小时缩短为3.0995s 。

2.缺点
(1)本模型在模型的求解过程中,对一些可接受范围内的误差直接进行了忽略,因而对于结果的精确性还是会有一定的影响。

(2)本模型是建立在一些假设中的,所有实用性受到了限制,在实际生产中,如果可以把更多的一些因素考虑进去应该会更好。

在已假定的条件下,本模型的优化结果是好的。

3推广
此模型有较强的应用价值。

工程中往往因为某个零件的选取不当,而影响产品的参数,使可靠性降低,造成了极大的经济损失。

所以需综合考虑零件成本和质量,以求获得最大的经济效益。

本模型具有广泛的适用性,很容易加以推广。

模型中的设计变量
x i i (,,,)=127 可以推广到n 个的情形,即设计变量x i n R i n (,,,)=∈12 ,其中设计空间R n
是一个n 维空间。

本模不仅适用于粒子分离器参数的设计,而且也可用于类似的机构、零部件、工艺设备等的基本参数的设计问题;容差等级同样可推广应用。

参考文献
【1】
韩之俊,姚平中,《概率与统计》,国防工业出版社,1985 【2】
陈宝林,《最优化理论与算法》,清华大学出版社,1989 【3】
裘宗燕,《数学软件系统的应用及程序设计》,北京大学出版社,1994 【4】
许波,《Matlab 工程数学应用》,清华大学出版社,2001
附录:matlab代码:
function f=result
%穷举108种容错等级组合求解全局最优解
fval=inf;
tic
%Bmin=[2 3 3 3 3 3 2];
%Xmin
B(1)=2;
B(5)=3;
for i=2:3
B(2)=i;
for j=1:3
B(3)=j;
for t=1:3
B(4)=t;
for g=1:3
B(6)=g;
for m=1:2
B(7)=m;
[fv,x]=getcost(B);
if fv<fval
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=simulation(MU,B)
%用随机法和计算的结果进行模拟比较
for i=1:10000
y(i)=Yfun(getparaX(MU,B));
end;
[f,xi] = ksdensity(y);
plot(xi,f); % 画经验概率密度曲线
hold on;
y0=Yfun(MU);
fc=getfcY(MU,B);
%{
x = normrnd(y0,fc,1,10000);
[f1,xj] = ksdensity(x);
plot(xj,f1,'r');
%}
x0=min(y):0.01:max(y);
y=((2*pi)^0.5*fc)^(-1)*exp(-(x0-y0).^2/2/fc^2);
plot(x0,y,'r');
%{
x=min(y):0.01:max(y);
yg=gaussmf(x,[fc,y0]);
plot(x,yg,'r');
%}
title('对照图');
gtext('注:蓝线为对x随机取样求得的y分布');
gtext('红线为根据模型计算出的y分布');
xlabel('y');
ylabel('y的概率密度');
hold off;
function [f,x]=getcost(B)
%在给定容差等级的情况下求最优的标定值,使得Y的均值为y0的情况下,方差最小MU=[0.1 0.3 0.1 0.1 1.5 16 0.75];%给定初始的标定值
options=optimset('LargeScale','off','Display','off');%,'Tolx',1.0000e-032);
[x,fval]=fmincon('getfcY',MU,[],[],[],[],[],[],'mycon',options,B);
x,B,f=cost(x,B)
function [c,ceq]=mycon(MU,B)
%求最优标定值时的约束条件
%c为不等式约束
%ceq为等式约束
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)
%当标定值为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;
f=f+50;
switch (B(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)
%对于所给的标定值和容差求Y的方差
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.05/3;
end;
if B(i)==3
sigma(i)=MU(i)*0.1/3;
end;
end;
x1=MU(1);x2=MU(2);x3=MU(3);x4=MU(4);x5=MU(5);x6=MU(6);x7=MU(7);
%求Y对各变量的偏导的评分与对应的方差乘积之和
f=(pd1(x1,x2,x3,x4,x5,x6,x7)*sigma(1))^2;f=f+(pd2(x1,x2,x3,x4,x5,x6,x7)*sigma(2))^2 ;
f=f+(pd3(x1,x2,x3,x4,x5,x6,x7)*sigma(3))^2;f=f+(pd4(x1,x2,x3,x4,x5,x6,x7)*sigma(4)) ^2;
f=f+(pd5(x1,x2,x3,x4,x5,x6,x7)*sigma(5))^2;f=f+(pd6(x1,x2,x3,x4,x5,x6,x7)*sigma(6)) ^2;
f=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)
%Y对x1的偏导
f=8721/50/x5*(x3/(x2-x1))^(17/20)*((1-131/50*(1-9/25/...
(x4/x2)^(14/25))^(3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)+...
148257/1000*x1/x5/(x3/(x2-x1))^(3/20)*((1-131/50*(1-9/25/(x4/x2)...
^(14/25))^(3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)*x3/(x2-x1)^2;
function f=pd2(x1,x2,x3,x4,x5,x6,x7)
%Y对x2的偏导
f=-148257/1000*x1/x5/(x3/(x2-x1))^...
(3/20)*((1-131/50*(1-9/25/(x4/x2)^(14/25))^(...
3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)*x3/(x2-x1)^...
2+8721/100*x1/x5*(x3/(x2-x1))^(17/20)/((1-131/50*(1-9/25/...
(x4/x2)^(14/25))^(3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)*(24759/31250*...
(1-9/25/(x4/x2)^(14/25))^(1/2)/(x4/x2)^(2/5)*x4/x2^2+3799/1250*(1-9/25/...
(x4/x2)^(14/25))^(3/2)*(x4/x2)^(4/25)*x4/x2^2)/x6/x7;
function f=pd3(x1,x2,x3,x4,x5,x6,x7)
%Y对x3的偏导
f=148257/1000*x1/x5/(x3/(x2-x1))^(3/20)*((1-131/50*...
(1-9/25/(x4/x2)^(14/25))^(3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)/(x2-x1);
function f=pd4(x1,x2,x3,x4,x5,x6,x7)
%Y对x4的偏导
f=8721/100*x1/x5*(x3/(x2-x1))^(17/20)/((1-131/50*(1-9/25/(x4/x2)^(14/25))^...
(3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)*(-24759/31250*(1-9/25/(x4/x2)^(14/25))^...
(1/2)/(x4/x2)^(2/5)/x2-3799/1250*(1-9/25/(x4/x2)^(14/25))^(3/2)*(x4/x2)^(4/25)/x2)/ x6/x7;
function f=pd5(x1,x2,x3,x4,x5,x6,x7)
%Y对x5的偏导
f=-8721/50*x1/x5^2*(x3/(x2-x1))^(17/20)*((1-131/50*(1-9/25/(x4/x2)^(14/25))^(3/2)*( x4/x2)^(29/25))/x6/x7)^(1/2);
function f=pd6(x1,x2,x3,x4,x5,x6,x7)
%Y对x6的偏导
f=-8721/100*x1/x5*(x3/(x2-x1))^(17/20)/((1-131/50*(1-9/25/(x4/x2)^(14/25))^...
(3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)*(1-131/50*(1-9/25/(x4/x2)^(14/25))^(3/2)*(x4/x2 )^(29/25))/x6^2/x7;
function f=pd7(x1,x2,x3,x4,x5,x6,x7)
%Y对x7的偏导
f=-8721/100*x1/x5*(x3/(x2-x1))^(17/20)/((1-131/50*(1-9/25/...
(x4/x2)^(14/25))^(3/2)*(x4/x2)^(29/25))/x6/x7)^(1/2)*...
(1-131/50*(1-9/25/(x4/x2)^(14/25))^(3/2)*(x4/x2)^(29/25))/x6/x7^2;
function f=getP(MU,B)
%当标定值为MU,容差等级为B时,求正品、次品、废品的概率
yb=Yfun(MU);
fc=getfcY(MU,B);
%syms x0 u a0;yy=subs(((2*pi)^0.5*a0)^(-1)*exp(-(x0-u)^2/2/a0^2),'u',yb);
%yy=subs(yy,'a0',fc);
%y0=1.5;
f(2)=jf1(yb,fc);
f(3)=jf2(yb,fc);
%f(1)=0;
%f(2)=(cdf('normal',y0+0.3,yb,fc) -cdf('normal',y0+0.1,yb,fc))*2 ;
%f(3)=2*cdf('normal',y0-0.3,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)
%Y的表达式
f=174.42*x(1)/x(5)*(x(3)/(x(2)-x(1)))^0.85*...
((1-2.62*(1-0.36*(x(4)/x(2))^(-0.56))^1.5*...
(x(4)/x(2))^1.16)/x(6)/x(7))^0.5;
function f=geteveryP(MU,B,iter)
%利用标定值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
f(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;
f
function f=getparaX(MU,B)
%利用标定值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),sigma0(i)/3);
% while( ~(f(i)>(MU(i)-sigma0(i)) & f(i)<(MU(i)+sigma0(i)))) %f(i)=normrnd(MU(i),sigma0(i)/3);
%end; end;。

相关文档
最新文档