南京信息工程大学数学与统计学院，江苏 南京

收稿日期：2021年5月21日；录用日期：2021年6月9日；发布日期：2021年6月25日

摘要

本文针对多元正态总体均值向量偏移问题，将双抽样策略与可变样本容量控制图相结合，构建了双抽样可变样本容量T²控制图(DSVSS T²图)。根据过程的实际运行状态，构建了马尔可夫链转移概率矩阵。综合考虑过程正常和异常两种状态下的运行成本、误报警损失、抽样成本及维修费用，结合Lorenzen-Vance经济模型，提出DSVSS T²控制图的经济设计模型并分析了均值向量变化对决策变量的影响。

关键词

双抽样，可变样本容量，马尔可夫链，经济优化

Economic Design Optimization of DSVSS T² Control Chart

Peng Chen, Bin Zhang^*, Xiu Chen

College of Mathematics and Statistics, Nanjing University of Information Science & Technology, Nanjing Jiangsu

Received: May 21^st, 2021; accepted: Jun. 9^th, 2021; published: Jun. 25^th, 2021

ABSTRACT

Aiming at the problem of mean vector shift of multivariate normal population, this paper combines double sampling strategy with variable sample size control chart, and constructs double sampling variable sample size control chart (DSVSS T² chart). According to the actual running state of the process, the Markov chain transition probability matrix is constructed. Considering the operation cost, false alarm loss, sampling cost and maintenance cost of normal and abnormal process, based on Lorenzen-Vance economic model, an economic design model of DSVSS T² control chart is proposed, and the influence of mean vector change on decision variables is analyzed.

Keywords:Double Sampling, Variable Sampling Interval, Markov Chain, Economic Design

Copyright © 2021 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

对于一个复杂产品，通常包含两个或两个以上质量特性，且质量特性之间存在相关关系，采用单变量控制图对每个质量特性分别进行监控往往并不适合 [1]。同时监控多个质量特性的控制方案能够准确反映生产状态，及时发现过程异常。Hotelling (1947)第一次提出多变量T²控制图，可以有效监测多质量特性的生产过程中发生的偏移 [2]。

传统T²控制图的参数(抽样间隔h，样本容量n，控制限k)都是确定的，这类控制图统称为静态FSR控制图。静态FSR控制图对过程较小的偏移不敏感，难以满足精益的生产环境需求。基于此，Reynolds (1988)首次提出可变抽样间隔(VSI)控制图 [3]。随后，变采样间隔(VSI)、变样本容量(VSS)、变样本量和采样间隔(VSSI)、变参数(VP)等四种自适应监控方案得到了广泛研究。Aparisi (2001)将自适应方案加入多变量控制图，研究了可变采样间隔(VSI) T²控制图 [4]。Faraz (2009)提出了可变样本容量(VSS) T²控制图 [5]。Mahadik and Shirke (2011)构建了联合可变样本容量和可变抽样区间(VSSI)的T²控制图 [6]。Saniga (2015)研究了变参数(VP) T²控制图，并将其性能与其他自适应T²图进行了比较 [7]。

Daudin (1992)将双抽样思想加入到控制图中，即通过两阶段抽样来监控过程，并与传统的控制图进行对比，发现其表现出更好的统计特性，能够减少抽样数量，提高监测效率 [8]。Carot (2002)将双抽样和可变样本间隔相结合，提出DSVSI $\bar{X}$ 控制图 [9]。PH Lee (2012)构建了DSVSI $\bar{X}$ 控制图的经济优化方法 [10]。Mehdi Katebi (2020)将双抽样和自适应方案同时应用于多元控制图，提出了DSVSSI T²控制图，将所提出的T²控制图与现有的标准T²控制图和其他自适应多变量T²图进行比较，发现DSVSSI T²图能更快地检测平均向量中小到中等的偏移 [10]。

自适应控制图可变抽样间隔在实际生产中往往难以操控，使用DSVSS控制图，可以降低质量检测人员的操控难度。在过程出现异常时，样本量的扩大能够获得更多的过程信息。万强(2020)基于计件数据，提出双抽样可变样本间隔np图的经济设计方法 [11]。本文将双抽样可变样本容量的思想应用于多质量特性计量数据的过程监控中，提出多质量特性的DSVSS T²控制图，构建了DSVSS T²控制图的经济优化模型。

2. 双抽样可变样本容量T²控制图

本文首先将多质量特性的抽样过程分为两个阶段，包含三种不同的检测水平。第一阶段采用两种检测水平，即小样本检测水平和正常样本检测水平；第二阶段采用大样本检测水平。两个阶段如图1所示。

图1. DSVSS T²控制图两个阶段

现假定生产过程有p个需要监控的质量特性， $X_1,X_2,\cdots ,X_p$，质量特性向量 $X=\left(X_1,X_2,\cdots ,X_p\right)$ 服从p维正态分布，均值向量 $\mu =\left({\mu }_1,{\mu }_2,\cdots ,{\mu }_p\right)$，协方差为 $\Sigma$。抽取容量为n的样本，计算每个样本的样本均值 $\overline{{{\rm X}}_i}=\left(\overline{X_1},\overline{X_2},\cdots ,\overline{X_p}\right)$，计算样本统计量T²， $T^2=n{\left(\overline{X_i}-{\mu }_{\text{0}}\right)}^{\prime }\Sigma \left(\overline{X_i}-{\mu }_{\text{0}}\right)$，其中 ${\mu }_0$ 为过程受控状态下的均值向量，若计算的T²超出控制限，则过程被标记为失控。

DSVSS T²控制图具有 $\left(n_1,n_2,n_3\right)$，$\left({\omega }_1,{\omega }_2\right)$，$\left(k_1,k_2\right)$，h八个设计参数，其中n₁，n₂分别表示一阶段抽样的小样本容量和正常样本容量，n₃表示二阶段抽样的样本容量， $n_1<n_2<n_3$ ； ${\omega }_{\text{1}}$ 为扩大抽样样本容量的警告限，ω₂为进入二阶段抽样的警告限， ${\omega }_{\text{1}}<{\omega }_{\text{2}}$ ；k₁，k₂分别表示两个阶段的控制限， $k_2<k_1$。h为固定抽样间隔。为方便描述，现记： $I_1=\left(0,{\omega }_1\right)$，$I_2=\left({\omega }_1,{\omega }_2\right)$，$I_3=\left({\omega }_2,k_1\right)$，$I_4=\left(k_1,\infty \right)$，$I_5=\left(0,k_2\right)$，$I_6=\left(k_2,\infty \right)$。操作步骤如下：

步骤一：抽取样本量为 $n_1$ 的小样本，计算样本统计量 $T_{11}^2$，进入步骤二。

步骤二：1) 样本统计量 $T_{11}^2$ 落在 $I_1$，则过程受控，返回步骤一。

2) 样本统计量 $T_{11}^2$ 落在 $I_2$，则过程受控，进入步骤三。

3) 样本统计量 $T_{11}^2$ 落在 $I_3$，则过程受控，进入步骤五。

4) 样本统计量 $T_{11}^2$ 落在 $I_4$，则认为过程失控，应立即寻找失控原因。

步骤三：抽取样本量为 $n_2$ 的大样本，计算样本统计量 $T_{12}^2$，进入步骤四。

步骤四：1) 样本统计量 $T_{12}^2$ 落在 $I_1$，则过程受控，返回步骤一。

2) 样本统计量 $T_{12}^2$ 落在 $I_2$，则过程受控，返回步骤三。

3) 样本统计量 $T_{12}^2$ 落在 $I_3$，则过程受控，进入步骤五。

4) 样本统计量 $T_{12}^2$ 落在 $I_4$，则认为过程失控，应立即寻找失控原因。

步骤五：进入二次抽样阶段，立即抽取样本量为 $n_3$ 的样本，计算样本统计量 $T_2^2$，进入步骤六。

步骤六：1) 样本统计量 $T_2^2$ 落在 $I_5$，则过程受控，返回步骤三。

2) 样本统计量 $T_2^2$ 落在 $I_6$，则认为过程失控，应立即寻找失控原因。

其中 $T_{11}{}^2=n_1{\left(\overline{{{\rm X}}_1}-{\mu }_{\text{0}}\right)}^{\prime }\Sigma \left(\overline{{{\rm X}}_1}-{\mu }_{\text{0}}\right)$，$\overline{X_1}={\displaystyle \underset{i=1}{\overset{n_1}{\sum }}{X}_i}$ ； $T_{12}{}^2=n_2{\left(\overline{X_2}-{\mu }_{\text{0}}\right)}^{\prime }\Sigma \left(\overline{X_2}-{\mu }_{\text{0}}\right)$，$\overline{X_2}={\displaystyle \underset{i=1}{\overset{n_2}{\sum }}{X}_i}$ ； $T_2{}^2=n_3{\left(\overline{X_3}-{\mu }_{\text{0}}\right)}^{\prime }\Sigma \left(\overline{X_2}-{\mu }_{\text{0}}\right)$，$\overline{X_3}={\displaystyle \underset{i=1}{\overset{n_3}{\sum }}{X}_i}$。

3. 转移概率矩阵

马尔可夫链作为一种随机过程，在控制图的相关领域中，常用来预测控制图各个状态的转移概率，利用状态转移矩阵去计算控制图的统计指标，从而对控制图进行统计和经济优化。基于本文的双抽样可变样本容量T²控制图，为简单起见，假设 ${\mu }_0$ 和 $\Sigma$ 已知，则T²服从自由度为p的中心卡方分布。若过程出现异常，均值发生偏移，T²服从非中心参数为 $\eta$ 的非中心卡方分布。现将第一阶段的样本统计量 $T_{11}^2$ 和 $T_{12}^2$ 统称为 $T_1^2$，则可将过程分为以下四种状态：

状态一： $T_1^2\in I_1$

状态二： $T_1^2\in I_2$

状态三： $T_1^2\in I_3\cap T_2^2\in I_5$

状态四： $T_1^2\in I_4$，$T_1^2\in I_3\cap T_2^2\in I_6$

其中状态四为吸收态，马尔可夫链转移矩阵P如下：

$P=\left(\begin{array}{cccc}{P}_{11}& P_{12}& P_{13}& P_{14}\\ P_{21}& P_{22}& P_{23}& P_{24}\\ P_{31}& P_{32}& P_{33}& P_{34}\\ 0& 0& 0& 1\end{array}\right)$

其中 $P_{ij}$ 为状态i转移到状态j的概率，因此

$P_{11}=P\left(T_1^2\in I_1|n_1;h\right)=P\left(T_1^2\le w_1\right)=F\left(w_1,p,\eta =n_1{\delta }^2\right)$

$P_{1\text{2}}=P\left(T_1^2\in I_{\text{2}}|n_1;h\right)=P\left(w_1\le T_1^2\le w_{\text{2}}\right)=F\left(w_{\text{2}},p,\eta =n_1{\delta }^2\right)-F\left(w_{\text{1}},p,\eta =n_1{\delta }^2\right)$

$\begin{array}{c}{P}_{1\text{3}}=P\left(T_1^2\in I_{\text{3}}\text{\&}{T}_{\text{2}}^2\in I_{\text{5}}|n_1;h\right)=P\left(w_{\text{2}}\le T_1^2\le k_{\text{1}}\right)\times P\left(T_2^2\le k_{\text{2}}\right)\\ =\left[F\left(k_{\text{1}},p,\eta =n_1{\delta }^2\right)-F\left(w_{\text{2}},p,\eta =n_1{\delta }^2\right)\right]\times F\left(k_{\text{2}},p,\eta =n_3{\delta }^2\right)\end{array}$

$P_{14}=1-{\displaystyle \underset{i=1}{\overset{3}{\sum }}{P}_{1i}}$

$P_{\text{2}1}=P_{\text{3}1}=P\left(T_1^2\in I_1|n_{\text{2}};h\right)=P\left(T_1^2\le w_1\right)=F\left(w_1,p,\eta =n_{\text{2}}{\delta }^2\right)$

$P_{\text{22}}=P_{\text{32}}=P\left(T_1^2\in I_{\text{2}}|n_{\text{2}};h\right)=P\left(w_1\le T_1^2\le w_{\text{2}}\right)=F\left(w_{\text{2}},p,\eta =n_{\text{2}}{\delta }^2\right)-F\left(w_{\text{1}},p,\eta =n_{\text{2}}{\delta }^2\right)$

$\begin{array}{c}{P}_{\text{23}}=P_{\text{33}}=P\left(T_1^2\in I_{\text{3}}\text{\&}{T}_{\text{2}}^2\in I_{\text{5}}|n_{\text{2}};h\right)=P\left(w_{\text{2}}\le T_1^2\le k_{\text{1}}\right)\times P\left(T_2^2\le k_{\text{2}}\right)\\ =\left[F\left(k_{\text{1}},p,\eta =n_{\text{2}}{\delta }^2\right)-F\left(w_{\text{2}},p,\eta =n_{\text{2}}{\delta }^2\right)\right]\times F\left(k_{\text{2}},p,\eta =n_3{\delta }^2\right)\end{array}$

$P_{\text{2}4}=P_{\text{3}4}=1-{\displaystyle \underset{i=1}{\overset{3}{\sum }}{P}_{\text{2}i}}=1-{\displaystyle \underset{i=1}{\overset{3}{\sum }}{P}_{\text{3}i}}$

其中 $F\left(x,p,\eta \right)$ 为非中心卡方分布的累计概率分布函数，p为自由度，即需要监控的质量特征的个数， $\eta$ 为非中心参数，计算公式为 $\eta =n{\delta }^{\text{2}}$，其中 $\delta$ 为多变量均值向量漂移量。

4. 统计指标计算

统计指标常用来测评控制图的性能表现。主要有受控状态下过程平均运行链长 $AR{L}_0$,平均周期时间 $AT{S}_{\text{0}}$ 以及一个控制周期的平均抽样数 $AN{I}_{\text{0}}$ 和失控状态下平均周期时间 $AT{S}_1$ 以及一个控制周期的平均抽样数 $AN{I}_1$，计算公式分别如下：

$AR{L}_0=b^{\prime }{\left(I-Q_0\right)}^{-1}l$

$AT{S}_0=b^{\prime }{\left(I-Q_0\right)}^{-1}{h}_0$

$AN{I}_0=b^{\prime }{\left(I-Q_0\right)}^{-1}N$

$AT{S}_1=b^{\prime }{\left(I-Q_1\right)}^{-1}{h}_0$

$AN{I}_1=b^{\prime }{\left(I-Q_1\right)}^{-1}N$

其中 $b^{\prime }=\left(b_{\text{1}},b_2,b_3\right)$ 为初始状态概率矩阵，为避免在过程开始就出现异常，一般设 $b^{\prime }=\left(\text{0},\text{1},\text{0}\right)$ ；I为3 × 3单位矩阵； $Q_{\text{0}}$ 为多变量均值不发生漂移时 $\left(\delta =\text{0}\right)$ 转移概率矩阵P删除吸收态后的3 × 3矩阵， $Q_{\text{1}}$ 为多变量均值发生漂移时 $\left(\delta >\text{0}\right)$ 转移概率矩阵P删除吸收态后的3 × 3矩阵； $l^{\prime }=\left(1,1,1\right)$ ； $h$ 为取样间隔，由于该控制图采样间隔固定，则取样间隔向量 ${h^{\prime }}_{\text{0}}=\left(h,h,h\right)$ ；N为取样个数向量,则 $N^{\prime }=\left(n_1,n_2,n_3\right)$。

5. 经济模型构建

为方便该控制图经济成本的计算，现假设如下：

· 控制图从受控状态开始运行；

· 控制图受控时间服从均值为 $\text{1/}\lambda$ 的指数分布；

· 控制图发生警报后，便停止运行，查找异常原因并完全维修；

· 二次抽样间隔为0。

本文采用Lorenzen-Vance经济模型 [11]，首先计算控制图的单个过程周期时间 $E\left(T\right)$。其中包括受控时间x，其服从均值为 $\text{1/}\lambda$ 的指数分布；失控时间AATS，由 $AT{S}_1$ 减去过程漂移与之前最后一次抽样之间

的时间期望值 $\tau$，因为受控时间x服从均值为 $\text{1/}\lambda$ 的指数分布， $\tau =E\left(x|0<x<h\right)=\frac{1-{\text{e}}^{-\lambda h}\left(\lambda h+1\right)}{\lambda \left(1-{\text{e}}^{-\lambda h}\right)}$，$AATS=AT{S}_1-\tau$ ；查找过程异常并完全维修的时间 $\epsilon$，是一个常量，则有

$E\left(T\right)=\frac{\text{1}}{\lambda }+AATS+\epsilon$

将控制图从开始运行至报警记为一个运行周期，则控制图在一个运行周期内的成本，主要由五个部分组成：

1) 受控状态下的过程成本

记受控状态下单位时间成本为 $c_{\text{1}}$，因为受控时间服从均值为 $\text{1/}\lambda$ 的指数分布，则受控状态下的过程成本为 $c_1/\lambda$。

2) 失控状态下的过程成本

记失控状态下的单位时间成本为 $c_{\text{2}}$，失控时间为AATS，则失控状态下的过程成本为 $c_{\text{2}}\times AATS$。

3) 误报警成本

记寻找一个误报警的成本为 $c_3$，受控状态下的期望抽样次数为s，受控下的平均抽样次数等于过程开始到过程均值发生偏移前最后一次抽样的时间长度除以平均抽样间隔，由于抽样间隔固定，则

$s=\frac{1/\lambda -\tau }{E\left(h\right)}=\frac{1/\lambda -\tau }{h}$，误报警的概率为 $1/AR{L}_0$，则误报警成本为 $\frac{c_{\text{3}}s}{AR{L}_0}$。

4) 抽样和检测成本

记一次抽样检测的固定成本为 $c_{\text{4}}$，则单个过程周期抽样检测的固定成本为 $c_{\text{4}}\left(s+AATS/h\right)$ ；单个样品检测成本为 $c_{\text{5}}$，受控状态下，平均周期时间和平均抽样数分别为 $AT{S}_{\text{0}}$ 和 $AN{I}_{\text{0}}$，则受控状态下的平均

单位时间抽样数为 $\frac{AN{I}_{\text{0}}}{AT{S}_{\text{0}}}$，当受控时间服从均值为 $\text{1/}\lambda$ 的指数分布时，受控状态下的样本抽样与检测成本为 $c_{\text{5}}\times 1/\lambda \times \frac{AN{I}_{\text{0}}}{AT{S}_{\text{0}}}$ ；失控状态下，平均周期时间和平均抽样数分别为 $AT{S}_{\text{1}}$ 和 $AN{I}_{\text{1}}$，则失控状态下的平均单位时间抽样数为 $\frac{AN{I}_{\text{1}}}{AT{S}_{\text{1}}}$，失控时间为AATS，失控状态下的样本抽样与检测成本为 $c_{\text{5}}\times AATS\times \frac{AN{I}_{\text{1}}}{AT{S}_{\text{1}}}$ ；则抽样和检测的总成本为 $c_4\left(s+\frac{AATS}{h}\right)+c_{\text{5}}\left(1/\lambda \times \frac{AN{I}_{\text{0}}}{AT{S}_{\text{0}}}+AATS\times \frac{AN{I}_{\text{1}}}{AT{S}_{\text{1}}}\right)$。

5) 查找异常原因并完全维修的成本

记寻找异常原因并完全维修所需成本为 $c_6$，为一个常量。

综上，一个运行周期内的总成本 $E\left(C\right)$ 为：

$E\left(C\right)=\frac{c_{\text{1}}}{\lambda }+c_{\text{2}}\times AATS+\frac{c_{\text{3}}s}{AR{L}_0}+c_4\left(s+\frac{AATS}{h}\right)+c_{\text{5}}\left(1/\lambda \times \frac{AN{I}_{\text{0}}}{AT{S}_{\text{0}}}+AATS\times \frac{AN{I}_{\text{1}}}{AT{S}_{\text{1}}}\right)+c_{\text{6}}$

则单位时间的运行成本函数 $E\left(A\right)=\frac{E\left(C\right)}{E\left(T\right)}$，该控制图的经济优化模型的目标函数为：

$\begin{array}{l}\min E\left(A\right)\\ s.t.0<n_1<n_2<n_3\le n_{\max }\\ n_1,n_2,n_3\in N^{+}\\ 0<k_2<k_1\\ 0<w_l<w_n<k_1\end{array}$

优化目标是确定控制图参数 $\left(w_1,w_2,k_1,k_2,n_{\text{1}},n_2,n_3,h\right)$ 使得该控制图单位时间的运行成本最小。

6. 数值计算

假定某生产过程需要同时监控的质量特性有3个，即 $p=3$。系统异常运行时仅导致均值发生偏移，多变量均值向量偏移量 $\delta =0.9$。系统受控状态下的运行时间服从均值为 $1/\lambda =100$ 的指数分布。本文的模型参数初始值参照Lorenzen和Vance [12] 在其文献中的设定值，即受控状态下单位时间成本 $c_1=114.24$，失控状态下的单位时间成本 $c_2=949.2$，寻找一个误报警的成本 $c_3=977.4$，一次抽样检测的固定成本 $c_4=2$，单个样品成本 $c_5=4.22$，寻找异常原因并完全维修所需成本 $c_6=1086$，查找过程异常并完全维

修的时间 $\epsilon =5/60$。将上述数据代入运行成本函数 $E\left(A\right)=\frac{E\left(C\right)}{E\left(T\right)}$，利用MATLAB遗传算法工具箱，得到8个决策变量 $\left(n_1,n_2,n_3,h,w_1,w_2,k_1,k_2\right)$ 和单位时间运行成本 $E\left(A\right)$ 的最优解如表1所示。

表1. 最优决策变量及目标函数值

7. 灵敏度分析

为了研究均值向量漂移对决策变量及目标函数的影响，在上节基础之上进行灵敏度分析，如表2所示。

表2. 均值向量漂移量δ对决策变量的影响

由表2可知随着漂移量δ的增大，第一阶段抽样的样本容量n₁，n₂不变，第二阶段抽样的样本容量n₃减小，抽样间隔h减小，警戒限w₁，w₂增加，控制限k₁，k₂增加，单位时间运行成本 $E\left(A\right)$ 减小，说明对于DSVSS T²控制图而言，当生产过程均值向量发生微小偏移时单位时间成本较高，这是因为尽管自适应控制图对微小偏移的反应足够及时，能有效避免过程异常时的生产损失，但当微小偏移发生时，需要抽取更多的样本去获得生产过程的信息，抽样成本将会增加，单位时间成本将增加。

8. 结论

本文提出的DSVSS T²控制图可以同时监控多个关键质量特性，将双抽样(DS)策略和可变样本容量(VSS)联合，提高了控制图在检测平均向量中位移的灵敏度且易于操控。利用马尔可夫链方法导出控制图统计指标的计算公式并构建单位时间成本模型，提出双抽样可变样本容量(DSVSS) T²控制图的优化设计方法。灵敏度分析表明随着均值向量漂移量递增，一阶段抽样样本量不变，二阶段抽样样本量递减，控制图的抽样间隔递减，警戒限和控制限递增，单位时间成本递减。

文章引用

陈 鹏,张 斌,陈 秀. DSVSS T²控制图经济设计模型
Economic Design Optimization of DSVSS T² Control Chart[J]. 应用数学进展, 2021, 10(06): 2186-2193. https://doi.org/10.12677/AAM.2021.106228

参考文献