基于正态变换的变样本容量Shewhart比例控制图设计

周思阳,张 莹 ,柯楚贤

(1.交通运输部水运科学研究所,北京 100088;2.武汉理工大学 交通与物流工程学院,湖北 武汉 430063)

0 引言

控制图作为一种有效的统计过程控制(Statistical Process Control, SPC)工具,被广泛应用于科学研究和工业生产。在一些生产场景中,成分间的比例关系是一种关键质量特性,需要对其进行监控以保证产品质量的稳定性,这些场景大致分为两类:①生产过程中混合物的成分之间存在特定的比例关系,如食品生产中原料的配比[1]和港口配煤作业中原煤的配比[2];②生产过程中涉及物理或化学反应等特定操作前后产品成分的比例,如在电池回收过程中可回收部分的重量比例[3]。因此近年来,如何利用控制图监控两个变量之间的比例关系逐渐成为研究热点。

在比例控制图研究中,CELANO等[1]首先在样本容量n=1的情况下对两个正态变量之比进行监控,提出Shewhart-RZ 控制图,并对其统计性能进行了研究;随后,这项研究由CELANO等[4]扩展到样本容量n>1的情形;NGUYEN等[5]提出单边Shewhart-RZ控制图,并利用线性协变量模型研究了测量误差对其性能的影响;ZAIDI等[6]和TRAN等[7]将对两个变量间比例的监控扩展到对成分数据的监控,从而可以同时对混合物中的多种成分进行监控。当比例特性偏移较大时,Shewhart-RZ控制图能较快发现过程失控,但由于其仅能利用当前样本信息,当比例特性的偏移较小时,往往需要较长时间才能发现过程失控。研究中通常采用附加运行规则、与合格品链长(Conforming Run Length,CRL)控制图相结合、变抽样区间和变样本容量等策略中的一种或多种对Shewhart控制图进行优化,以改善其在偏移较小时的性能。其中针对Shewhart-RZ控制图,CELANO等[8]和NGUYEN等[9]分别将双边和单边Shewhart-RZ控制图与合格品链长控制图结合;TRAN等[10]和TRAN[11]采用不同附加运行规则对双边Shewhart-RZ控制图进行优化;NGUYEN等[12]采用变抽样区间策略对单边Shewhart-RZ控制图进行优化设计。目前采用变样本容量策略对单边Shewhart-RZ控制图进行优化的研究比较缺乏,有些研究采用两种记忆性控制图,即指数移动加权平均(Exponentially Weighted Moving Average, EWMA)和累计和(CUmulative SUM, CUSUM)控制图,对二元正态变量间的比例关系进行监控[13-17],计算结果表明其性能优于Shewhart-RZ控制图。

分析发现,动态控制图中的变样本容量策略能有效提高Shewhart型控制图的性能,并具有易于实施的良好特性。然而,目前比例控制图研究多集中于静态控制图领域,对动态控制图的研究较少,且均未考虑变样本容量策略在比例控制图中的应用。本文提出一种单边变样本容量Shewhart比例控制图,即VSS Shewhart-RZ*控制图,来监控二元正态变量间的比例关系,并比较其与单边Shewhart-RZ控制图的性能,结果表明相同情况下VSS Shewhart-RZ*控制图具有更小的ARL。

1 比例Z分布

假设W=(X,Y)T为二元正态随机向量,W～N(μW,∑W),μW为均值向量,∑W为协方差矩阵,具体计算公式如下:

(1)

(2)

式中ρ为变量X,Y之间的相关性系数。此时变量X,Y的变异系数可定义为γX=σX/μX,γY=σY/μY,同时记X,Y的标准差之比为ω=σX/σY。

令X,Y的比值为Z,Z=X/Y,由于无法得到比例Z的累计分布函数的准确表达式,CELANO等[4]对其近似估计为

(3)

式中:Φ(·)为标准正态分布的累积分布函数;A和B为参数z,γX,γY,ω,ρ的函数,

(4)

(5)

比例Z的概率分布函数可近似估计为

(6)

式中φ(·)为标准正态分布的概率分布函数。

比例Z累积分布函数的反函数可由式(3)推导得到

(7)

式中C1,C2,C3是参数γX,γY,ω,ρ的函数:

(8)

(9)

(10)

为了在生产过程中监控比例Z,在时刻i(i=1,2,…)抽取样本容量为n的样本{Wi,1,Wi,2,…,Wi,n},Wi,j=(Xi,j,Yi,j)T～N(μW,i,∑W,i)(j=1,…,n)为二元正态随机向量,其均值向量和协方差矩阵分别为

(11)

(12)

与CELANO等[4]的研究相同,对样本数据做以下假设:①不同抽样时刻下产品的规格可以发生变化,即当i≠k时,有μW,i≠μW,k,∑W,i≠∑W,k;②对于不同的样本,变量X,Y的变异系数γX,γY为已知常数,即对样本{Wi,1,Wi,2,…,Wi,n}可得σX,i=γX×μX,i,σY,i=γY×μY,i;③当过程处于受控状态时,变量X,Y的比例为已知常数z0=μX,i/μY,i(i=1,2,…)。此时观测统计量为

(13)

(14)

(15)

2 VSS Shewhart-RZ*控制图

2.1 VSS Shewhart-RZ*控制图的实施流程

由于比例Z的分布不对称,本文设计了两个单边VSS Shewhart-RZ*控制图分别监控比例向上或向下的偏移,以避免性能损失。在抽样时样本容量n有两种可能的取值,即小样本容量nS和大样本容量nL,nS

(1)在进行控制图优化设计时,需要计算n(i)=nS时的警戒限和控制限{(LWLS,LCLS),(UWLS,UCLS)},以及n(i)=nL时的警戒限和控制限{(LWLL,LCLL),(UWLL,UCLL)},因此需要花费更多的计算资源和时间对控制图的参数组合进行寻优。

(2)在实际使用过程中,使用者需要检查每个样本点所处的位置,以判断是否需要更换控制限和警戒限组合,增加了使用者的操作难度。

(16)

此时上单边VSS Shewhart-RZ*控制图的控制限和警戒限可固定为

UCL=z0×KU,

(17)

UWL=z0×WU。

(18)

与之类似,下单边VSS Shewhart-RZ*控制图的控制限和警戒限为

LCL=z0×KD,

(19)

LWL=z0×WD。

(20)

式中:KU,KD为控制限系数,WU,WD为警戒限系数,且满足WU≤KU和KD≤WD。

上单边VSS Shewhart-RZ*控制图的区域划分如图1a所示,工作机制如下:

下单边VSS Shewhart-RZ*控制图的区域划分如图1b所示,工作机制与上单边控制图的工作机制类似,不再赘述。

2.2 VSS Shewhart-RZ*控制图的性能指标

在变样本容量控制图研究中,常采用平均运行长度ARL和平均样本容量ASS作为控制图的性能指标。ARL指从过程开始运行到控制图发出过程失控信号时所经过的平均样本数,受控状态下的ARL0越大越好,此时控制图虚发警报的概率降低,可避免多余的停机检查;失控状态下的ARL1越小越好,此时控制图漏发警报的概率降低,可以快速检测到异常因素。ASS为样本容量的均值,ASS=E(n(i)),当控制图的性能达到需求时,ASS越小越好,更小的平均样本数表示抽样检测成本更小。比较VSS型控制图与固定样本容量的静态控制图时,使受控状态下的ARL0和ASS0保持一致,比较失控状态下的ARL1和ASS1,更小的ARL1或ASS1表示拥有更好的性能。

假设过程失控时,过程中存在异常因素使受控状态下的比例z0偏移至z1=τ×z0,其中τ>0,τ的取值表示偏移的大小和方向,τ>1对应比例向上的偏移,τ<1对应比例向下的偏移,过程处于受控状态时τ=1。同时考虑失控状态下,变量(X,Y)间的相关系数也可能发生偏移,即从受控状态下的ρ=ρ0偏移至ρ=ρ1。

采用马尔可夫链推导单边VSS Shewhart-RZ*控制图的性能指标。将单边VSS Shewhart-RZ*控制图建模为一个含有3种状态的马尔可夫链,其中前两种为受控状态,分别对应使用小样本容量n(i)=nS的情形和使用大样本容量n(i)=nL的情形,第3种为失控状态,在马尔可夫链模型中为吸收态。马尔可夫链的一步转移概率矩阵

(21)

式中:Q为2×2的转移概率矩阵,1=(1,1)T,r为2×1的向量,r=1-Q1,即矩阵P中每一行的和为1;pS(n(i))和pL(n(i))分别为从当前状态转移至小样本容量和大样本容量的概率,n(i)∈{nS,nL}表示当前状态。对于上单边VSS Shewhart-RZ*控制图,

(22)

(23)

对于下单边VSS Shewhart-RZ*控制图,

(24)

(25)

此时控制图的ARL为

ARL=qT(I-Q)-11。

(26)

式中:I为2×2的单位矩阵,q为初始概率向量,1=(1,1)T。当设定第1组样本为小样本容量,即n(1)=nS时,q=(1,0)T;当设定第1组样本为大样本容量,即n(1)=nL时,q=(0,1)T。

计算ASS时需要先将转移矩阵P转换成P*,当n(1)=nS时,

(27)

当n(1)=nL时

(28)

与矩阵P不同,P*没有吸收态,当马尔可夫链到达失控状态时,消除过程中的异常因素后立刻重启,使马尔可夫链恢复至初始状态。由P*定义的马尔可夫链有稳态概率向量π=(πS,πL,πOOC)T,可由下式求得:

(29)

其中稳态概率(πS,πL,πOOC)与样本容量(nS,nL,n(1))相联系。当n(1)=nS时,矩阵R可由矩阵P*先转置,再将对角线元素减1,并将第1行元素替换为1,得到

(30)

当n(1)=nL时,矩阵R可由矩阵P*先转置,对角线元素减1后,将第2行元素替换为1,得到

(31)

此时ASS可由下式计算:

ASS=(nS,nL,n(1))π。

(32)

2.3 VSS Shewhart-RZ*控制图的优化设计

VSS Shewhart-RZ*控制图的优化目标是在特定的比例偏移τ和相关系数ρ1下,确定最优的控制限系数KU(或KD)、警戒限系数WU(或WD)和样本容量组合(nS,nL),使受控状态下控制图的平均运行长度ARL0=L0,平均样本容量ASS0=n0,同时使失控状态下的ARL1最小。其中:L0与受控状态下控制图误发警报的概率α密切相关,L0=1/α;失控状态下的ARL1与漏发警报的概率β相关,ARL1=1/(1-β)。控制图的优化过程即在给定可接受的误发警报概率α下,使失控状态下控制图漏发警报的概率β最小。由于在实际应用中样本容量的最大值一般会有限制,与之前的研究[21-23]类似,本文设置样本容量可取的最大值为31,则nS,nL需要满足1≤nS

minARL1(nS,nL,WU,KU,γX,γY,ρ1,τ,z0)。

s.t.

ARL0(nS,nL,WU,KU,γX,γY,ρ0,τ=1,z0)=L0;

ASS0(nS,nL,WU,KU,γX,γY,ρ0,τ=1,z0)=n0;

WU≤KU;

1≤nS

nS,nL∈N+。

(33)

下单边VSS Shewhart-RZ*控制图的数学模型与上单边类似,仅警戒限和控制限之间大小关系的约束略有不同,需将其改写为KD≤WD。其中L0的取值一般需要事先给定,基于3σ原理常取L0=370.4,本文为了方便比较,与之前研究[4]保持一致,取ARL0=200。n0的取值为实际过程中受控状态下样本容量的期望值,为方便比较,本文取n0∈{5,15}。

由于非凸的混合整数非线性规划问题属于NP-hard问题,求解非常困难,目前缺少有效的求解算法,本文将整数变量分离,使问题化简为对多个非线性方程组来求解,具体过程为(如图2):

(1)确定γX,γY,ρ0,ρ1,τ,z0,n0,L0等参数的取值。

(2)令nS=1。

(3)令nL=n0+1。

(4)在当前参数下,采用MATLAB软件中的fslove函数求解由两个等式约束构成的非线性方程组,得到一组警戒限和控制限变量的取值(WU,KU)或(WD,KD)。

(5)根据当前(nS,nL,WU,KU)或(nS,nL,WD,KD)的取值计算ARL1指标。

(6)保持nS不变,令nL=nL+1。

(7)重复(4)～(6),直至nL=31。

(8)令nS=nS+1。

(9)重复(3)～(8),直至nS=n0-1。

3 数值分析及性能比较

在不同参数设置下求得单边VSS Shewhart-RZ*控制图的最优决策变量组合,计算ARL1性能指标,并与固定样本容量(Fixed Sample Size,FSS)的单边Shewhart-RZ控制图[5]、单边Synthetic-RZ控制图[9]和单边VSI Shewhart-RZ控制图[12]进行比较。同时由2.2节可知,n(1)取值会对控制图性能产生影响,目前变样本容量控制图文献中大多只研究了n(1)=nS的情况,对n(1)=nL的情形研究较少,因此将在不同n(1)取值下研究单边VSS Shewhart-RZ*控制图的性能。为方便对比,采用与文献[4]的参数组合,如n0∈{5,15},γX∈{0.01,0.2},γY∈{0.01,0.2},ρ∈{0.0,±0.4,±0.8},τ∈{0.9,0.95,0.98,0.99,1.01,1.02,1.05,1.1},考虑γX=γY,γX≠γY,ρ0=ρ1和ρ0≠ρ1等场景。

图3和图4所示为单边Shewhart-RZ控制图和单边VSS Shewhart-RZ*控制图(包括n(1)=nS和n(1)=nL两种情形,没有特别指明时默认n(1)=nS),在过程失控且X,Y的相关系数保持不变(ρ0=ρ1)时,ARL1的值随比例偏移变化的情况。图3主要展示变量X,Y的变异系数相等(γX=γY)时的情形,其中左边两列图为γX=γY=0.01的情形,右边两列图为γX=γY=0.2的情形。图4主要展示γX≠γY时的情形,其中左边两列图为(γX=γY)=(0.01,0.2)的情形,右边两列图为(γX=γY)=(0.2,0.01)的情形。单边VSS Shewhart-RZ*控制图与单边Synthetic-RZ控制图、单边VSI Shewahrt-RZ控制图的性能比较在图3中进行了展示,其中单边Synthetic-RZ控制图的性能指标为ARL,单边VSI Shewahrt-RZ控制图的性能指标为平均报警时间ATS。从图3和图4可得以下结论:

(1)单边VSS Shewhart-RZ*控制图性能受(γX,γY)和ρ影响。变量(X,Y)的变异系数(γX,γY)取值越小,控制图的性能越好。例如当ρ=-0.8,n0=5,τ=0.9时,若(γX,γY)=(0.01,0.01)则ARL1=1,若(γX,γY)=(0.2,0.2)则ARL1=19.8。ρ>0时控制图的性能优于ρ<0时的性能,即变量(X,Y)之间呈正相关时控制图的性能更好。例如当(γX,γY)=(0.2,0.2),n0=5,τ=0.9时,若ρ=-0.8则ARL1=19.8,若ρ=0.8则ARL1=2.6。

(2)对于相同百分比的偏移ΔZ=100×|τ-1|,当γX=γY时单边VSS Shewhart-RZ*控制图对下偏移的检测性能略好于对上偏移的检测性能。例如当(γX,γY)=(0.2,0.2),n0=5,ρ=-0.8,τ=0.99(或τ=1.01),即ΔZ=1%时,可得下单边控制图的ARL1=170.3(上单边控制图ARL1=170.6),如图3所示。当γX≠γY时,单边VSS Shewhart-RZ*控制图对上下偏移的性能趋势取决于变异系数(γX,γY)之间的大小关系。当γX<γY时,控制图对下偏移的检测性能更好;当γX>γY时,控制图对上偏移的检测性能更好。例如当(γX,γY)=(0.01,0.2),n0=5,ρ=-0.8,τ=0.99(或τ=1.01)时,ARL1=137.8(或ARL1=157.4);当(γX,γY)=(0.2,0.01),n0=5,ρ=-0.8,τ=0.99(或τ=1.01)时,ARL1=157(或ARL1=138.3),如图4所示。

图5所示为单边Shewhart-RZ控制图和单边VSS Shewhart-RZ*控制图,当过程失控X,Y的相关系数发生变化(ρ0≠ρ1,n0=5)时,ARL1随比例偏移变化的情况。受控状态下ρ0=±0.4,受异常因素影响偏移至ρ1=0.5×ρ0或ρ1=2×ρ0。从图5可以观察到以下趋势:增强负相关性会提高单边VSS Shewhart-RZ*控制图对比例偏移的敏感性,例如(γX,γY)=(0.01,0.2),ρ1=0.5×ρ0=-0.2,τ=0.99时,ARL1=148.1,如果ρ1=ρ0=-0.4则ARL1=136.6,如果ρ1=2×ρ0=-0.8则ARL1=117.5;增强正相关性会降低控制图对比例偏移的敏感性,例如(γX,γY)=(0.01,0.2),ρ1=0.5×ρ0=0.2,τ=0.99时,ARL1=123.3,如果ρ1=ρ0=0.4则ARL1=134.1,如果ρ1=2×ρ0=-0.8则ARL1=160.5。

在图3～图5中,通过对比可知,对于单边VSS Shewhart-RZ*控制图,第1次抽样采用大样本容量(n(1)=nL)的性能优于第1次抽样采用小样本容量(n(1)=nS),在(γX,γY)=(0.01,0.01)、偏移较小时性能优势尤为显著。n(1)=nL时,虽然单边VSS Shewhart-RZ*控制图拥有更小的ARL1,但是同时也拥有更大的ASS1,例如(γX,γY)=(0.01,0.01),n0=5,ρ0=ρ1=-0.8,τ=0.99,如果n(1)=nS则ARL1=4.7,ASS1=8.7,如果n(1)=nL则ARL1=1.7,ASS1=28.9。

在图3～图5中,通过对比可知,单边VSS Shewhart-RZ*控制图的性能优于单边FSS Shewhart-RZ控制图,说明变样本容量策略起到了很好的优化效果。当变量X,Y的变异系数均较小时,VSS策略对小偏移的优化效果更为明显;当变量X,Y的变异系数有一者较大时,VSS策略对大偏移的优化效果更为明显。在一些场景下,FSS Shewhart-RZ控制图与VSS Shewhart-RZ*控制图都能取得很好的效果,使ARL1=1,但是VSS Shewhart-RZ*控制图可以使ASS1

在图3中通过对比可知,在以下两种情况下,单边VSS Shewhart-RZ*控制图的性能优于单边Synthetic-RZ控制图:

(1)变量X,Y的变异系数较小,即(γX,γY)=(0.01,0.01)时,n(1)=nL的单边VSS Shewhart-RZ*控制图性能优于单边Synthetic-RZ控制图,且偏移较小时优势更为显著。例如(γX,γY)=(0.01,0.01),ρ0=ρ1=-0.8,n0=5,τ=1.01时,单边Synthetic-RZ控制图ARL1=5.2,n(1)=nL的单边VSS Shewhart-RZ*控制图ARL1=1.8。

(2)当变量X,Y的变异系数较大,即(γX,γY)=(0.2,0.2),比例偏移较大(τ>1.1或τ<0.9)且受控状态下的样本容量较小(n0=5)时,单边VSS Shewhart-RZ*控制图同样拥有更好的性能。例如(γX,γY)=(0.2,0.2),ρ0=ρ1=0,n0=5,τ=1.1时,单边Synthetic-RZ控制图ARL1=14.1,n(1)=nL的单边VSS Shewhart-RZ*控制图ARL1=8.4。

在VSS Shewhart-RZ*控制图与VSI Shewhart-RZ控制图的性能对比中,以平均报警时间ATS=E(h)×ARL作为性能指标。其中E(h)为平均抽样间隔,对于VSS Shewhart-RZ*控制图,其抽样间隔h为固定值,取h=1,则ATS=ARL。因为VSI Shewhart-RZ控制图在设计时满足受控状态下平均抽样区间E(h)0=1的约束,所以可以通过比较失控状态下VSI Shewhart-RZ控制图ATS1指标与VSS Shewhart-RZ*控制图ARL1指标的数值,来对比两者的性能。图3展示了抽样区间组合(hS,hL)=(0.1,1.1),τ∈{0.95,0.98,0.99,1.01,1.02,1.05}时,VSI Shewhart-RZ控制图失控状态下的ATS1[12]。可知在以下两种情况下,单边VSS Shewhart-RZ*控制图的性能优于VSI Shewhart-RZ控制图:

(1)当变量X,Y的变异系数较小,且相关系数ρ≤0、偏移和受控状态下的样本容量较小时,例如 (γX,γY)=(0.01,0.01),n0=5 ,ρ0=ρ1=-0.4 ,τ=1.01时,单边VSI Shewhart-RZ控制图ATS1=5.9,单边VSS Shewhart-RZ*控制图ARL1=4.0。

(2)当变量X,Y的变异系数较大,且偏移较大时,例如(γX,γY)=(0.2,0.2)n0=5,ρ0=ρ1=0,τ=0.95时,单边VSI Shewhart-RZ控制图ATS1=64.7,单边VSS Shewhart-RZ*控制图ARL1=53.5。

4 案例说明

表1 案例样本数据

针对τ=1.01的偏移,在ARL0=200,n0=5的条件下对采用的4种控制图进行优化设计,得到最优参数。对于单边Shewhart-RZ控制图计算可得UCL=1.015 4;对于单边Synthetic-RZ控制图,参考文献[9]计算可得合格品链长子图的控制限为H=4,单边Shewhart-RZ子图的控制限为UCL=1.010 7;对于单边VSI Shewhart-RZ控制图,参考文献[12]计算可得抽样间隔组合(hS,hf)=(0.1,1.1),警戒限UWL=1.007 5,控制限UCL=1.015 4;对于单边VSS Shewhart-RZ*控制图(n(1)=nS),根据第2章计算可得样本容量组合(nS,nL)=(3,20),警戒限UWL=1.159 0,控制限UCL=2.575 8。

根据样本数据及控制图参数绘制4种比例控制图如图6所示,其中Synthetic-RZ控制图仅绘制了单边Shewhart-RZ子图。由图6可知,在前10个受控样本下,所有控制图均未误发警报;当生产过程失控后,根据4种控制图各自的运行规则,单边Shewhart-RZ控制图在第14个样本点发出警报,报警时间tSh=4;单边Synthetic-RZ控制图在第12个样本点发出警报,报警时间tSyn=2;单边VSI Shewhart-RZ控制图在第14个样本点发出警报,报警时间tVSI=2.4;单边VSS Shewhart-RZ*控制图在第11～14样本点发出警报,报警时间tVSS=1。当过程处于失控状态时,4种比例控制图均在较短时间内发出了警报,说明了其有效性。其中VSS Shewhart-RZ*控制图在第11个点处及时发出警报,报警时间最短,且连续4点报警,具有很好的稳定性,效果最好。实际应用中,生产者应在VSS Shewhart-RZ*控制图发出警报后及时停产检查,在消除引起比例偏移的异常因素后再恢复生产,以免产生大量不合格品,保证产品质量的稳定性。

VSS Shewhart-RZ*控制图性能优秀,且易理解、易实施,可广泛应用于工程实践。首先其原理简单易理解,即检测统计量距离控制限较远,过程失控的可能性较小时,采用小样本容量抽样,以减少抽样成本;检测统计量距离控制限较近,过程失控的可能性较大时,采用大样本容量抽样,以更好地了解当前过程状态。其次,VSS策略在应用中仅需改动每次抽样的产品数量,便于操作,可广泛应用于人工抽样检测、自动化抽样检测等作业场景。

5 结束语