Economic and Economic-Statistical Design of Multivariate Bayesian Control Chart with Variable Sampling Intervals

Due to the proper performance of Bayesian control chart in detecting process shifts, it recently has become the subject of interest. It has been proved that on Bayesian and traditional control charts, the economic and statistical performances of the variable sampling interval (VSI) scheme are superior to those of the fixed ratio sampling (FRS) strategy in detecting small to moderate shifts. This paper studies the VSI multivariate Bayesian control chart based on economic and economic-statistical designs. Since finding the distribution of Bayesian statistic is t complicated, we apply Monte Carlo method and we employ artificial bee colony (ABC) algorithm to obtain the optimal design parameters (sample size, sampling intervals, warning limit and control limit). In the end, this case study is compared with VSI Hotelling’s T2 control chart and it is shown that this approach is more desirable statistically and economically.


Introduction
The source of shifts in process means are random and assignable causes. Control charts are used to monitor a process in order to detect assignable causes. The traditional control chart used for monitoring the mean of a single variable is X-bar and for multiple variables is T 2 -Hotelling control chart which was proposed by Hotelling (1947). Girshick and Rubin (1952) considered a Bayes approach to a quality control model. Also Taylor (1965 and1967) showed that the traditional control techniques were not optimal. Calabrese (1995) presented Bayesian Process Control for Attributes. Tagaras and Nikolaidis (2002) evaluated the economic performance of various adaptive control schemes to derive conclusions about their relative effectiveness. They used the analysis that concentrated on Bayesian control charts for monitoring the process mean in finite production run. Recently, other researchers such as Maikis (2009) and Tavakoli et al. (2014 have attended the Bayesian control charts. The traditional sampling strategy in control charts is the fixed ratio sampling (FRS) scheme wherein samples of fixed size are obtained at fixed intervals between successive samples. Although the FRS control charts have a good performance in detecting large shifts in the process mean, it is proved that its efficiency to detect small and moderate shifts in the process mean is week. To improve this weakness, variable sampling interval (VSI) is the one procedure. In this scheme which was first proposed by Reynolds et al. (1988), sampling intervals are varied such that the next sample should be taken after a short time if a sample statistic shows an indication of trouble, and the next sample will be taken after long time if there is no such indication. Seif and Sadeghifar (2015) applied this method in special case of Hotelling's T 2 control chart. When the strategy of the process designer to obtain design parameters (sample size, sampling intervals, warning limit and control limit) is maximizing the time between false alarms and minimizing the time to detect merely an off-target process, control chart is designed statistically. Unlike statistical designs (SD) in which only the statistical properties are taken into consideration, economic designs (ED) of control charts, is known as a pioneer in minimizing the expected cost per time unit without statistical considerations. Duncan (1956 and1971), Lorenzen and Vance (1986), Banerjee and Rahim (1988) and Torabian et al. (2010) are researchers that applied economic design on their control charts. Woodall (1986) expressed that the EDs have a high type I error probability which can lead to irrelevant process adjustment or a loss of trust in the control process. To overcome this weakness, Saniga (1989) proposed economicstatistical design (ESD). He included the statistical constraints upon the ED and showed that not only its economic properties are as good efficient as EDs, but also its statistical properties are acceptable. These constraints due to the designer needs, can be Type I error, power of the chart, adjusted average time to signal (AATS) and average number of false alarm (ANF). Many authors such as Zhang and Berardi (1997), Prabhu et al. (1997), Yang and Rahim (2005), Yeong and Yanjing (2015) and Faraz et al. (2009) considered this scheme in their researches.
In this paper, we proposed VSI multivariate Bayesian control chart and. According to this sampling method, the performance of considered control improved. Since finding the distribution of Bayesian statistic is complicated, we utilized Monte Carlo method. The optimum design parameters are obtained by using artificial bee colony (ABC) algorithm that is newer, faster and more accurate than Genetic algorithm which is commonly used. We also compared VSI and FRS Multivariate Bayesian charts based on ESD and showed that the performance of VSI scheme is more satisfying. Additionally, we compared the economic performance of VSI and FRS charts for various input parameters separately. This comparison shows that using economic statistical design can control statistical features while not incurring significant costs. The paper is organized as follows. After the introduction to the subject matter, in section 2 we define the VSI multivariate Bayesian control chart. In the next section, the cost model proposed by Costa and Rahim (2001) is used to build economic and economic-statistical designs of VSI multivariate Bayesian control chart. Optimization method and ABC algorithm are briefly reviewed in section 4. In section 5, the design optimum parameters for economic and economic-statistical design of VSI multivariate Bayesian control chart are obtained and for a better comparison, the values of VSI Hotelling's T 2 approach are also presented. Finally, we drew some conclusions.

VSI Multivariate Bayesian Control Chart
Assume that we want to monitor and maintain the output of a process which is described by the p quality characteristic. The considered assumptions are as follows: 1. In the process quality control, the VSI Bayesian control scheme is employed to monitor p related quality characteristics mean. 2. The p quality characteristics are multivariate normally distributed with mean vector and covariance matrix . 3. The process starts in the in-control state with mean 0 and the length of time that the process stays in this state is exponentially distributed with mean 1  ⁄ . 4. By single assignable cause the process mean may shift from 0 to 1 = 0 + . 5. The covariance matrix  is constant over time and drifting process from this assumption is not a subject of this paper. 6. The process is not self-corrective. That is, when the process shifts to an out-of-control state, it returns to the in-control condition only by the management intervention upon appropriate corrective actions. 7. After the shift, the process mean remains off target until the assignable cause is eliminated. 8. During the search for an assignable cause, the process is shut down. Consider that X is a 1 × random vector wherein jth element is the jth quality characteristic. We assume that it is multivariate normally distributed, denoted by ( 0 , ), where 0 is the 1 × mean vector and  is the × covariance matrix of X. We assume that the process starts in the in-control state, the length of time that process stays in the in-control state has an exponential distribution with mean 1  ⁄ and by one assignable cause may shift the mean of the quality characteristic (process mean) from 0 to 1 = 0 + where = [ 1 , 2 , … ] and ≠ 0. Since in practice, 0 and  are generally unknown, therefore in commissioning phase it is necessary to estimate them from m initial samples, each of size n. ̿ and ̅ are the sample estimates of 0 and  respectively, where and is a p dimensional random vector of the quality characteristic in which ijth element is the jth element of ith sample.
To monitor and control p correlated characteristics of the process outputs, in the implementation phase, we should take a random sample of size n from the production at per time t.
is the kth sample at time t, which is dimensional random vector. The statistic of the Hotelling's T 2 in this phase is as follow: where ̅̅̅ is the mean of s which is independent of the estimators of the parameters. In this case and in this phase, if > 1, 2 has a non-central F distribution with parameters p, ( − 1) − + 1 and ′  −1 . It is 2 is distributed as F distribution with p and ( − 1) − + 1 degrees of freedom.
We define the statistic of Bayesian control chart as the posterior probability that the process is in-control. It is denoted by  at the tth stage of sampling. Thus where 0 signifies all samples in the commissioning phase and U 1 , … , U t are values which obtain by all observation in the implementation phase till time t as follows: = ( − 1) − + 1 ( + 1)( − 1) 2 Theorem 1. As respect to the deffinition of  , it is equal as follow: where ( | = 0 ) is density function of given the process is in-control Also ( | = 1 ) is density functions of given the process is out-of-control .
proof: It is considered X (t) = {X 1 , X 2 , … , X i , … X t } as a collection of samples till time t. To simplicity, the joint density function of X (t) and μ t is considered at the point (x (t) , μ 0 ) by f(x (t) , μ t = μ 0 ).
The first term on the right hand of equation 2 is equal as follws: Similarly, the second term on the right hand of above equation is equal because, it is obvious that ( = 1 | −1 = 1 ) = 1. Therefore, As was expressed and equations 3 and 4 and due to the fact that ( = 0 | ( ) = ( ) ) = ( ( ) , = 0 ) , the theorem will be proofed. □ It is obvious that 0 <  < 1 and production continues as long as  > , where is a control limit. Also, ( | = We know that in the traditional sampling strategy (FRS) in control charts, samples are taken of fixed size with a fixed time interval between samples. Reynolds et al. (1988) showed that the variable sampling interval (VSI) is effective to improve the performance of FRS strategy. They split the chart into three regions: the central (safe) region, the warning region, and the action region. In this literature it was considered that if the current point falls in the safe region, there is a strong indication that the process is operating properly, so it is possible to wait longer than the usual time for the next sample (maximum h). If the sample point falls in the warning region, it is an indication that the process needs adjustment, so the next sample should be taken after a shorter period of time than usual (minimum h). Finally, a sample point falling in the action region is an indication that the process is operating in an out-of-control state. Thus, the production process should be stopped to search. Therefore, the parameters of a VSI control chart that should be determined come as , ℎ 1 , ℎ 2 , and where ℎ 1 > ℎ 2 . These regions and limits are shown in the following figure 1. The chosen procedure of sampling interval is as follows: 1. If  falls in d1 region the next sample is taken after ℎ 1 unit of time 2. If  falls in d2 region the next sample is taken after ℎ 2 unit of time 3. If  falls in d3 region, the chart will show that the process is out-of-control. However if = μ 0 , then the signal is a false alarm.

Economic and Economic-statistical design
In this paper, the cost model proposed by Costa and Rahim (2001) has been used to build our objective function. Also we considered certain assumptions about the process which we have studied. The assumptions that have been summarized below are relatively standard and most economic models incorporate them to some degree.

The cost function
According to this model the process cycle consists of the following phases, as illustrated in Figure 2: 1.
In control period, 2.
Time to detect an assignable cause, 4.
Time to repair an assignable cause.
These four steps form the basis of the renewal reward process that can be used to calculate the quality cycle cost per hour for a specified set of design parameters.

Figure 2: A quality cycle
Assume that ATC is the average time from the start of the production until the first signal after the process shifts and ANF is the average number of false alarms. Therefore, the expected length of a production cycle is given by where T 0 is the average time taken for a false alarm when the process is in control and T 1 is the average time to find and remove the assignable cause.
The expected net profit from a production cycle is given by where V 0 and V 1 are the average profit per hour earned when the process is operating in control and out-of-control, respectively, a 3 is the average cost to detect and remove the assignable cause, a 3 ′ is the average consequence cost of a false alarm, a 2 is the cost per an inspected item and ANI is the average number of inspected items per cycle. Now, based on Costa and Rahim (2001), the loss function E(L), or the average production cycle cost per hour is given by The goal of the ED VSI scheme is to determine the five design parameters (sample size n, sampling intervals h1 and h2, warning limit w and control limit k) which minimize ( ) without any statistical limitation. But based on ESD some statistical constraints are applied on ( ) as follow: Since the Bayesian statistic is complicated and finding its distribution is too difficult, we utilized Monte Carlo method to obtain the AATS, ANF and ANI. We know that in the traditional schemes, since their statistics have a specified distribution, computing these parameters and transition probability matrix is possible and easy. For example in X-bar chart, when quality characteristic has a normal distribution, the chart statistic, i.e. X ̅ , is normally distributed too and we can apply Markov chain approach to compute mentioned parameters and matrix.
According to the definition of ( ), it is sufficient to compute three parameters ATC, ANF and ANI. Also it should be noted that: Now, let be the time from the start of production until the first signal after the process shifts. Thus = ( ). If we simulate production cycle N times, then for large values of N, a good approximation of ATC is ∑ =1 , where is the value of in the ith cycle. Other parameters can be computed in the same way.

Optimization method and ABC approach
As it was expressed the purpose of the ED VSI scheme is to find the five design parameters (n, h1, h2, wand k) which minimize ( ) without any statistical limitation. But based on ESD some statistical constraints have been applied on ( ) as follows: which we fixed 0 = 7 and 0 = 0.5. Therefore, based on this design, the optimization problem can be written as follow: Also it should be noted that the process and cost parameters (input parameters) are 2 , 3 ′ , 3 , 0 , 1 , 0 , 1 ,  . Clearly, it is desirable to minimize false alarms and detection of process shifts as quickly as possible. Thus it is desirable to have small ANF and AATS values. This nonlinear constrained optimization problem is solved through the artificial bee colony (ABC) algorithm to obtain the optimum design parameters of the Bayesian control chart. In 2005, Karaboga proposed an artificial bee colony (ABC), which is based on a particular intelligent behavior of honeybee swarms. ABC is developed based on inspecting the behaviors of real bees in finding nectar and sharing the information of food sources with the bees in the hive. Agents in ABC are the Employed bees, the Onlooker bees and the Scout bees. The employed bees stays on a food source and provides the neighborhood of the source in its memory. The onlooker bees gets the information of food sources from the employed bees in the hive and select one of the food sources to gather the nectar. The Scouts is responsible for finding new food, the new nectar and sources. Procedures of ABC are as follows: The probability of selecting a nectar source is: where is the probability of selecting the ith employed bee, S is the number of employed bees, is the position of the ith employed bee and ( ) is fitness value. The new position is calculated as follow: where is the position of the onlooker bee, t is the iteration number, is the randomly chosen employed bee, j is

Numerical comparison
To compare the performance of the VSI multivariate Bayesian control chart with other schemes, we coded the optimization and simulation algorithms of this scheme and the VSI Hotelling's T 2 in Matlab. The different values of input parameters that were proposed by Costa and Rahim (2001), is brought in Table 1. For each case of Table 1 and based on ED and ESD, the optimal design parameters and resulting loss for VSI multivariate Bayesian and VSI Hotelling's T 2 is computed, separately. Also, for better and exhaustive comparison the noted values of FRS multivariate Bayesian scheme based on ED and ESD which were shown by , are presented. The values of the optimal design parameters and resulting loss have been shown in the tables as follows: -The VSI multivariate Bayesian and FRS multivariate Bayesian approach based on ESD in Table 2, -The VSI multivariate Bayesian and FRS multivariate Bayesian approach based on ED in Table 3, -The VSI Hotelling's T 2 approach based on ESD in Table 4 and -The VSI Hotelling's T 2 approach based on ED in Table 5. In this research, all of the optimum values have been obtained with p=2.     All the cases of Table 1 are pertinent to some kinds of times and costs of process except rows 7, 10 and 11 which are related to the different kinds of shifts. In the row 1 and 6 the large shift is considered and in the rows 10, 11 and 13 small shifts are taken into consideration. Also in other rows moderate shifts are attended. Since the process times and costs are very different in each process production, it is more important that we pay attention to the performance and behavior of control charts in various shifts.
The results of VSI multivariate Bayesian control chart based on ESD and ED showed that this control chart based on ED, compared to ESD has a lower cost in some rows of Table 1 but under statistical properties in all the cases of shifts, costs and times of process (was shown in Table 1) is poor (Figurs 3 and 4). These results prove our hypothesis about the features of ESDs and the weakness of EDs which in EDs, the statistical properties despite of their importance in decision about production processes, are not considered. Also, these values show that even the cost under ESD in the rows 2, 6, 8, 9, 10 and 12 of Table 1 is better than ED and in the row 11, almost is as good as ED. It means that, the performance of ESD in small process shifts (rows 10, 11 and 13), when the length of time that the process stays in the in-control state is short (row 12) and in some cases of times and costs (rows 2, 6, 8 and 9) is better than ED policy. As the importance of statistical properties and due to the economic characteristic in ESD, it can be said that the performance of VSI Bayesian control chart based on ESD is more desirable than ED.  Table 2 show that the VSI compared to FRS multivariate Bayesian control chart based on ESD has a better statistical performance in all kinds of process shifts, costs and times except when the process has a large shift (row 7). Although economic features in VSI strategy in the cases which the process has small shifts (rows 10, 11 and 13) and rows 2, 6 and 12 has a more desirable performance, in other cases of moderate shifts and different times and costs of process a weaker performance rather than FRS policy is observable ( Figures 5 and 6). The values which are shown in Table 3 express that VSI policy in multivariate Bayesian control chart is economically effective when small shifts in process occur. In other case the performance of FRS strategy in multivariate Bayesian control chart based on ED is more favorable than VSI model. Therefore, the performance of VSI multivariate Bayesian control chart compared to FRS is economically and statistically more effective in small shifts of process.  Tables 4 and 5. As it was shown in these tables, the efficiency of VSI Hotelling's T 2 approach -statistically and economically-is more feeble than VSI Bayesian scheme. All the statistical and economic properties reported in tables are greater than those of VSI multivariate Bayesian approach while the smaller values of these characteristics are more desirable. This holds in the ED and ESD. Therefore, the performance of VSI multivariate Bayesian control chart compared to VSI Hotelling's T 2 approach is a lot more admirable under any circumstances (Figures 7 -9).

An illustrative example
This section illustrates ED and ESD of the VSI Bayesian control chart with an industrial example which is considered by Chen (2006). Consider two quality characteristics of soft drink fabrication process are the pressure inside the soft drink bottle and the gas volume presented in the drink. The in-control process mean and variance are unknown but they are estimated by 50 initial samples. The estimated values of μ 0 and ∑ are ̿ = (6.93, 3.86) and ̅ = [ 0.06 −0.04 −0.04 0.8 ], respectively. According to the previous runs, the shift will happen after 100 hours which it means  = 0.01. In this process, it has been determined that the cost for each inspected item is $5 ( 2 = 5). The average amount of time wasted searching for the assignable cause is about 2.5 hours (T0=2.5) and the average time to find and remove the assignable cause (which is given as T1) is equal to 1. The process is subject to the magnitude of shift one (d=1). The earned profit of the process is $500 (V0=500) for in-control period and about $100 (V1=100) for out-ofcontrol period. The average consequence cost of a false alarm is $500 ( 3 ′ = 500) and the average cost to detect and remove the assignable cause is $500 ( 3 ). Determining both ED and ESD separately, optimal design parameters and economic and statistical properties for VSI Hotelling's T 2 and VSI Multivariate Bayesian control charts will be obtained based on above mentioned values. The results are shown in Table 6 and 7 respectively.  Based on the results, it is clear that the performance of VSI Bayesian control chart is more desirable than X-bar control chart, economically and statistically; however, the x-bar chart is commonly used in the industry.
In this example, using Bayesian control chart based on ED, the operation saving per hour is

LOSS-ED
LOSS of Hotelling's T2 approach -ED LOSS of Bayesian approach -ED Also, using Bayesian control chart based on ESD, the signal of process shift will be received 0.862 hours faster than X-bar chart: − − = 0.862.

Conclusion
In this paper, we first proposed an economic and economic-statistical designs of VSI multivariate Bayesian control chart and then determined optimal design parameters based on them. In this research Mont Carlo method was used since finding the distribution of Bayesian statistic seemed difficult.
The optimal values of economic and statistical properties of VSI multivariate Bayesian control chart based on ESD and ED show that this scheme compered to VSI Hotelling's T 2 control chart based on ED and ESD has a more favorable performance in all cases and conditions of process shifts, costs and times. Also, these results proved our hypothesis about the features of ESDs and the weakness of EDs. The statistical characteristics were unfavorable in ED and economic property in ESD almost was as good as ED. So generally it can be said, the performance of VSI multivariate Bayesian control chart based on ESD was more desirable than ED model. Finally, the results show that features of VSI multivariate Bayesian control chart based on ESD is more desirable than ESD of FRS multivariate Bayesian control chart in all conditions. Also, the performance of VSI multivariate Bayesian scheme under ED was economically more effective than ED of FRS multivariate Bayesian approach when the process has small shifts.