# Building Cost Function 3D Benchmarks to Improve the Economic Statistical Design ofX

Authors: Santiago-Omar Caballero-Morales
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Shewhart or control charts are important Statistical Process Control (SPC) techniques
used for prompt detection of failures in a manufacturing process and minimization of production costs which are modelled with nonlinear functions (cost functions). Heuristic methods have been used to find the chart’s parameters integrated within the cost function that best
comply with economic and statistical restrictions. However heuristic estimation is highly
dependent on the size of the search space, the set of initial solutions, and the exploration
operators. In this paper the 3D analysis of the cost function is presented to more accurately
identify the search space associated with each parameter ofX
¯ control charts and to improve estimation. The parameters estimated with this approach were more accurate than those
estimated with Hooke and Jeeves (HJ) and Genetic Algorithms (GAs) under different failure
distributions. The results presented in this work can be used as a benchmark to evaluate
and improve the performance of other heuristic methods.


---

## Body

## 1. Introduction

Shewhart or control charts are tools used to determine whether or not a production process is in a state of statistical control and thus if the entities being produced are within quality requirements. These requirements are set by Upper and Lower Control Limits (UCL, LCL). If the quality attribute or feature from sampled entities (e.g., weight, length) is not within these limits, then the process is in an out-of-control state (nonconforming entities are being produced). In this case finding and correcting the “assignable cause” of this state (e.g., a failure) must be performed [1].A control chart requires three main parameters:n, the size of the sample; h, the length of the interval between samples; and k, the coefficient of the control limits [1]. These parameters are estimated based on economic and statistical restrictions because there are costs and times associated with sampling and searching assignable causes. Considering arbitrary values for these parameters without the economic and statistical factors may lead to unnecessary expenses in the manufacturing process.(i)
Small intervals between samples (smallh) would increase sampling frequency and the process’ cycle time with the associated production costs. Depending on the nature of the entities being produced, unnecessary sampling may lead to product loss.
(ii)
Close control limits (smallk) would increase the frequency of (possibly false) failure alarms and the associated costs of interrupting the manufacturing process. A false alarm may lead to rejection of products that do not necessarily have low quality (rejection of conforming entities).The Economic Statistical Design (ESD) of control charts considers the statistical requirements such as the probabilities of error Types I (α, detecting an out-of-control state when the process is fine) and II (β, not detecting an out-of-control state when the process is not fine) for the estimation of n, h, and k. These parameters are integrated within a cost function which includes the time and money costs associated with sampling and searching/repairing assignable causes [2–10]. Finding the most suitable values of n, h, and k that minimize the cost function considering the economic and statistical restrictions is not an easy task. This is because of the number of variables and the complexity of the cost function model. The ESD of control charts is considered a combinatorial optimization problem [11].Among the approaches, methods, or algorithms to find a suitable solution for the ESD of control charts, Hooke and Jeeves (HJ) [4, 12], Genetic Algorithms (GAs) [2, 5, 6, 11, 13, 14], Combinatorial Methods (CB) [7], Tabu Search (TS) [15], and Teaching-Learning Based Optimization (TLBO) [16] have been proposed. Heuristic methods however do not guarantee finding the optimal values for n, h, and k that minimize the multivariable cost function. This is because these methods or algorithms select the better solutions in every iteration step knowing only the solutions in the previous iteration steps [17]. Particularly, the search space defined by the initial set of solutions of these methods may be different; thus convergence to feasible solutions may be achieved at different times.Achieving efficient algorithms is important because the effectiveness of the ESD and cost minimization depends on the accuracy of determination of then, h, and k parameters [16]. In the absence of the optimal results the assessment of these algorithms is performed by comparing them with other algorithms for solving the same problem or case study. Execution time (speed), computational and space efficiency, accuracy, and precision are criteria for evaluating the performance of algorithms [18]. In order to evaluate the performance of algorithms considering these criteria the following approaches have been proposed.(i)
Development of benchmarks: these may consist of sets of optimal results, or sets of specifically designed problems to test the abilities (search strategies, operators) of the algorithms (e.g., benchmark problems) [17, 19].
(ii)
Algorithm visualization: it enables the externalization of the space exploration strategies of the algorithm for analysis and improvements [20].
(iii)
Characterization of the type of problems for which the algorithms have good performance [19].The contribution of this paper is the characterization and visualization of the cost function for the development of benchmarks for the ESD ofX
¯ control charts. The analysis and 3D visualization of the cost function are presented to provide more specific intervals for n, h, and k to identify the most suitable search space for the global optimum and also to provide an insight of the behavior of the cost function for different failure mechanisms. This information is important to provide a benchmark to evaluate and improve the performance of algorithms for the ESD of  X
¯ control charts.The cost function presented by Rahim and Banerjee [21] was used for this work as it considers general failure distribution in contrast to the well known cost function of Lorenzen and Vance [22] which considers an exponential distribution. The 3D visualization was achieved with support of reductions of the nonlinear factors associated with the cost function such as the failure probability distributions to make computations faster. In order to test the outcome of the 3D visualization a direct search was performed on the specific 3D search space identified with this approach. The parameters estimated within this space were statistically more accurate than those obtained with the HJ and GAs heuristic methods.The structure of this paper is as follows: in Section2 the details of the cost function and the failure probability distributions are presented. Then, in Section 3 the details of the method used to create the 3D visualization models are presented while in Section 4 sets of these visualizations (libraries) are presented. In Section 5 the application of the 3D visualizations to improve estimation of parameters is discussed. A comparison of performance with other methods as HJ and GA is presented. Finally in Section 6 the conclusions and future work are discussed.


## 2. Cost Function

The cost function models the interaction between the costs, times, failure patterns, and restrictions associated with a production cycle which represents the interval from the starting production time (in-control state) to the time when a failure caused by an assignable cause occurs. This cycle includes the time required to detect and repair the assignable cause [21].The objective of the ESD is to minimize the average cost per unit of time of a production cycle which is expressed as(1)
min
⁡
Z
=
E
C
E
T
,
where E
(
C
) is the Expected Cost per Cycle and E
(
T
) the Expected Cycle Length [4, 21, 23]. In the following sections the expressions for E
(
C
) and E
(
T
) with general failure distributions (which represent the failure patterns) are presented.


### 2.1. Expected Cost per Cycle and Cycle Length

The cost function considered in this work is an extension of the function presented by Lorenzen and Vance [22] which has been used and adapted by many researchers in the SPC field [24–30]. The ratio of (2) and (3) represents this cost function as presented by Rahim and Banerjee [21] for processes with general failure distribution and constant sampling intervals. In Table 1 the elements of (2) and (3) are described:
(2)
E
C
=
D
0
h
1
-
F
h
F
h
+
D
1
h
1
+
β
1
-
β
+
Y
α
1
-
F
(
h
)
F
(
h
)
+
(
a
+
b
n
)
1
F
(
h
)
+
β
1
-
β
+
τ
(
D
0
-
D
1
)
+
W
(3)
E
(
T
)
=
h
1
F
(
h
)
+
β
1
-
β
+
α
Z
0
1
-
F
(
h
)
F
(
h
)
+
Z
1
.Table 1
Definitions of elements of the cost function.
F
(
h
)
Probability of failure of a process at the end of an intervalh.
α
Probability of the error Type I: false positive detection of the out-of-control state.
β
Probability of the error Type II: false negative detection of the out-of-control state.
τ
Point of time withinh where the process changes from the in-control state to the out-of-control state.
Z
0
Time associated with finding an assignable cause of failure when there is a false positive detection of the out-of-control state.
Z
1
Time associated with finding an assignable cause of failure and restoring the process to an in-control state.
a
+
b
n
Sampling cost: constant (a) and variable (b) cost.
D
0
In-control state production cost per hour.
D
1
Out-of-control state production cost per hour:D
1
>
D
0.
Y
Cost of a false positive detection of the out-of-control state.
W
Cost associated with finding an assignable cause of failure and restoring the process to an in-control state.Whilen  and  h  are explicitly presented in the cost function, k is implicitly presented in α and β which are defined by the following expressions for the X
¯ control chart assuming normality [26, 31]:
(4)
α
=
2
1
-
ϕ
(
k
)
,
(5)
β
=
ϕ
(
k
-
δ
n
)
-
ϕ
(
-
k
-
δ
n
)
,
where δ is the magnitude of the change of the mean of the process to be detected. The mean of the process changes from μ to μ
1
=
μ
+
δ
σ where σ is the standard deviation. ϕ
(
v
) denotes the distribution function of a standard normal variable v.


### 2.2. Failure Probability Distributions

Failure density functions (f
(
t
)) are important elements for reliability prediction [32]. These functions are used to model the failure rate patterns of a process and must be considered for the estimation of production costs. For any density function the associated cumulative probability distribution function (F
(
t
)) represents the probability that a unit, randomly taken from a population, will fail at most in time t [33].If instead of taking one unitn units are taken at the end of a time interval h, F
(
t
)
→
F
(
h
) where F
(
h
) represents the probability of the process failing (changing to an out-of-control state) at the end of the sampling interval h [34]. Considering a failure density function f
(
t
), the failure probability distribution F
(
h
) can be estimated by the following expression:
(6)
F
(
h
)
=
∫
t
=
0
t
=
h
‍
f
(
t
)
d
t
.In this work the exponential and Weibull functions are considered. The exponential function is commonly used for the estimation of failure patterns and is defined by the following expression:(7)
f
E
(
t
;
λ
)
=
λ
e
-
λ
t
,
where λ is the main parameter of the density function and represents the known number of failures per unit of time.The Weibull density function is used for processes where causes of failures are assumed to be of different nature [35]. This function is defined by the following expression [36]:
(8)
f
W
(
t
;
c
,
m
)
=
m
c
-
m
t
m
-
1
e
-
t
/
c
m
,
where c represents the time when the process is likely to fail, identifying in this way the life expectancy of the process. c is known as thescale parameter and m as theform/shape parameter. When m
=
1 the Weibull density function is approximated to the exponential function with λ
=
1
/
c. Failure patterns which conform to a Weibull probability distribution with m
>
1 are considered to be age-related failures while those conforming to a Weibull distribution with m
=
1 are considered to be random failures [35].In order to compute the cost function defined by the ratio of (2) and (3), F
(
h
) must be computed for (7) and (8). Particularly for the Weibull function, the estimation of the element τ involves some complex integrals. In Section 3.3 the methodology performed for the computation of F
(
h
) and τ is presented.


## 2.1. Expected Cost per Cycle and Cycle Length

The cost function considered in this work is an extension of the function presented by Lorenzen and Vance [22] which has been used and adapted by many researchers in the SPC field [24–30]. The ratio of (2) and (3) represents this cost function as presented by Rahim and Banerjee [21] for processes with general failure distribution and constant sampling intervals. In Table 1 the elements of (2) and (3) are described:
(2)
E
C
=
D
0
h
1
-
F
h
F
h
+
D
1
h
1
+
β
1
-
β
+
Y
α
1
-
F
(
h
)
F
(
h
)
+
(
a
+
b
n
)
1
F
(
h
)
+
β
1
-
β
+
τ
(
D
0
-
D
1
)
+
W
(3)
E
(
T
)
=
h
1
F
(
h
)
+
β
1
-
β
+
α
Z
0
1
-
F
(
h
)
F
(
h
)
+
Z
1
.Table 1
Definitions of elements of the cost function.
F
(
h
)
Probability of failure of a process at the end of an intervalh.
α
Probability of the error Type I: false positive detection of the out-of-control state.
β
Probability of the error Type II: false negative detection of the out-of-control state.
τ
Point of time withinh where the process changes from the in-control state to the out-of-control state.
Z
0
Time associated with finding an assignable cause of failure when there is a false positive detection of the out-of-control state.
Z
1
Time associated with finding an assignable cause of failure and restoring the process to an in-control state.
a
+
b
n
Sampling cost: constant (a) and variable (b) cost.
D
0
In-control state production cost per hour.
D
1
Out-of-control state production cost per hour:D
1
>
D
0.
Y
Cost of a false positive detection of the out-of-control state.
W
Cost associated with finding an assignable cause of failure and restoring the process to an in-control state.Whilen  and  h  are explicitly presented in the cost function, k is implicitly presented in α and β which are defined by the following expressions for the X
¯ control chart assuming normality [26, 31]:
(4)
α
=
2
1
-
ϕ
(
k
)
,
(5)
β
=
ϕ
(
k
-
δ
n
)
-
ϕ
(
-
k
-
δ
n
)
,
where δ is the magnitude of the change of the mean of the process to be detected. The mean of the process changes from μ to μ
1
=
μ
+
δ
σ where σ is the standard deviation. ϕ
(
v
) denotes the distribution function of a standard normal variable v.


## 2.2. Failure Probability Distributions

Failure density functions (f
(
t
)) are important elements for reliability prediction [32]. These functions are used to model the failure rate patterns of a process and must be considered for the estimation of production costs. For any density function the associated cumulative probability distribution function (F
(
t
)) represents the probability that a unit, randomly taken from a population, will fail at most in time t [33].If instead of taking one unitn units are taken at the end of a time interval h, F
(
t
)
→
F
(
h
) where F
(
h
) represents the probability of the process failing (changing to an out-of-control state) at the end of the sampling interval h [34]. Considering a failure density function f
(
t
), the failure probability distribution F
(
h
) can be estimated by the following expression:
(6)
F
(
h
)
=
∫
t
=
0
t
=
h
‍
f
(
t
)
d
t
.In this work the exponential and Weibull functions are considered. The exponential function is commonly used for the estimation of failure patterns and is defined by the following expression:(7)
f
E
(
t
;
λ
)
=
λ
e
-
λ
t
,
where λ is the main parameter of the density function and represents the known number of failures per unit of time.The Weibull density function is used for processes where causes of failures are assumed to be of different nature [35]. This function is defined by the following expression [36]:
(8)
f
W
(
t
;
c
,
m
)
=
m
c
-
m
t
m
-
1
e
-
t
/
c
m
,
where c represents the time when the process is likely to fail, identifying in this way the life expectancy of the process. c is known as thescale parameter and m as theform/shape parameter. When m
=
1 the Weibull density function is approximated to the exponential function with λ
=
1
/
c. Failure patterns which conform to a Weibull probability distribution with m
>
1 are considered to be age-related failures while those conforming to a Weibull distribution with m
=
1 are considered to be random failures [35].In order to compute the cost function defined by the ratio of (2) and (3), F
(
h
) must be computed for (7) and (8). Particularly for the Weibull function, the estimation of the element τ involves some complex integrals. In Section 3.3 the methodology performed for the computation of F
(
h
) and τ is presented.


## 3. Building the 3D Benchmark

The ESD problem consists in finding the set of parametersn, h, and k that minimizes the cost function Z
=
E
(
C
)
/
E
(
T
)
=
f
(
n
,
h
,
k
). In this case there are four variables: three to be estimated (
n
,
h
,
k
) and one which is dependent on these variables (the associated Z value). Because of this situation a direct 3D visualization is not possible. Hence, the 3D benchmark consists of the integration of different 2D and 3D models. In this section the steps followed to build these model sets are presented.In order to perform the computations required to build the 2D and 3D models for the 3D benchmarks the following equipment was used: DELL Inspiron 15R with Intel Core i7-4500U CPU at 2.40 GHz with 8 GB in RAM. Programming was performed with MATLAB R2008a (version 7.6.0.324).


### 3.1. Finding Optimal Intervals forn, k, and δ

An initial link betweenn, k, α, and β is defined by (4) and (5) where β
≤
β
^ and α
≤
α
^. β
^ and α
^ are maximum error levels and commonly β
^
≤
0.20 and α
^
≤
0.20. At this point there is no influence from other elements or costs from Z.The first step is to find an optimal interval forn. For this purpose the approximation presented in [37] for n with a specific value for β is used:
(9)
1
-
β
^
=
ϕ
(
δ
n
-
k
)
,
(10)
ϕ
-
1
(
1
-
β
^
)
=
δ
n
-
k
,
(11)
n
=
ϕ
-
1
(
1
-
β
^
)
+
k
δ
2
.Whilen depends on k and β, α is dependent only on k. From (4)  k is expressed as
(12)
k
=
ϕ
-
1
(
1
-
α
2
)
.By replacing (12) into (11) the estimation of n now depends on α and β. For any α and β the interval for n can be obtained with the following expression:
(13)
n
=
ϕ
-
1
(
1
-
β
)
+
ϕ
-
1
(
1
-
α
/
2
)
δ
2
.The optimal interval forn now depends on required intervals for α and β. In this case the interval from 0.01 (minimum) to 0.20 (maximum) was considered for both error probabilities (e.g., α and β
≤
0.20).The second step is to estimate an optimal interval fork. This is performed directly from (13) as defined by (12). In Table 2 the k values estimated for α
≤
0.20 are presented where an inverse correlation between both parameters is observed.Table 2
k values associated with α
≤
0.20 assuming normality for the X
- control char.
α
k
α
k
0.0020
3.0902
0.1040
1.6258
0.0080
2.6521
0.1100
1.5982
0.0140
2.4573
0.1160
1.5718
0.0200
2.3263
0.1220
1.5464
0.0260
2.2262
0.1280
1.5220
0.0320
2.1444
0.1340
1.4985
0.0380
2.0749
0.1400
1.4758
0.0440
2.0141
0.1460
1.4538
0.0500
1.9600
0.1520
1.4325
0.0560
1.9110
0.1580
1.4118
0.0620
1.8663
0.1640
1.3917
0.0680
1.8250
0.1700
1.3722
0.0740
1.7866
0.1760
1.3532
0.0800
1.7507
0.1820
1.3346
0.0860
1.7169
0.1880
1.3165
0.0920
1.6849
0.1940
1.2988
0.0980
1.6546
0.2000
1.2816In Figure1 the 3D visualizations of (13) considering the optimal intervals for k and β and three δ levels are presented.3D visualization ofn considering different values for k, β, and δ.
(a)
δ
=
0.10
(b)
δ
=
0.50
(c)
δ
=
0.90Forα
=
(
0.00,0.20
) from (12) it is estimated that k
=
(
1.2816,3.0902
) (inverse interval). While the overall 3D pattern does not change through Figures 1(a)–1(c), the difference is mainly observed on the magnitude of n which is highly dependent on δ. In Figure 2 the 2D visualization of the intervals (minimum and maximum values) for n considering different values for δ across β and α
≤
0.20 is presented. The intervals presented in Figures 1 and 2 can be used for any cost function for X
¯ control charts assuming normality.Figure 2
2D visualization of the intervals forn considering different values for δ.Although in practiceδ is considered to be equal to 0.5, the overall effect of δ on n has not been presented in graphical form as in Figure 2. For this work three scenarios are considered for δ based on the levels defined in [2, 4, 31, 38]: δ
=
(
0.25,0.50,1.00
). This leads to the following optimal intervals for n based on data presented in Figure 2: n = (73–570, 19–143, 5–36). In Table 3 the associated intervals for n, k, and δ are summarized.Table 3
Intervals defined forn, k, and δ assuming normality for the X
- control chart.
Parameter
Lower bound
Upper bound
n (δ
=
0.25)
73
570
n (δ
=
0.50)
19
143
n (δ
=
1.00)
5
36
k
1.2816
3.0902


### 3.2. Setting the Constant Elements ofZ

The cost function defined by (2) and (3) has many elements besides the n, h, k, and F
(
h
) parameters which were described in Table 1. In the literature this cost function has been studied with a set of standard values for the constant elements a, b, Y, W, Z
0, Z
1, D
0, and D
1 for specific δ values [2, 4, 31, 38]. In Table 4 these sets are presented for each level of δ.Table 4
Standard costs and times forZ [2].
δ
Z
0 (hours)
Z
1 (hours)
D
0 ($)
D
1 ($)
W ($)
Y ($)
a ($)
b ($)
0.25
0.025
0.10
25.00
475.00
550.00
250.00
10.00
2.11
0.50
0.250
1.00
50.00
950.00
1100.00
500.00
20.00
4.22
1.00
0.500
10.00
100.00
1900.00
2200.00
1000.00
40.00
8.44


### 3.3. Estimating the Probability Distributions andτ

F
(
h
) for the exponential and Weibull density functions is obtained from (6) as follows:
(14)
F
E
(
h
)
=
∫
0
h
‍
f
E
(
t
)
d
t
=
1
-
e
-
λ
h
.
F
W
(
h
)
=
∫
0
h
‍
f
W
(
t
)
d
t
=
-
c
-
m
h
c
-
m
h
m
(
-
1
+
e
-
(
h
/
c
)
m
)
.An important element of (2) is τ which is defined as the point of time within h where the process changes from the in-control state to the out-of-control state. This element is estimated as
(15)
τ
=
∫
0
h
‍
t
f
(
t
)
d
t
∫
0
h
‍
f
(
t
)
d
t
=
∫
0
h
‍
t
f
(
t
)
d
t
F
(
h
)
.Considering the exponential distributionτ is expressed as
(16)
τ
E
=
∫
0
h
‍
t
f
E
(
t
)
d
t
F
E
(
h
)
=
-
(
-
1
+
e
-
λ
h
+
e
-
λ
h
λ
h
)
/
λ
1
-
e
-
λ
h
.Equation (16) is widely presented in the literature [5, 11, 26, 30]; however for the Weibull distribution an expression for τ
W for any m and c has not been found. In this work an expression for τ
W was obtained with support of the symbolic toolbox of MATLAB and additional algebraic reductions:
(17)
τ
W
=
∫
0
h
‍
t
f
W
(
t
)
d
t
F
W
(
h
)
=
R
-
c
-
m
h
/
c
-
m
h
m
(
-
1
+
e
-
(
h
/
c
)
m
)
.In (17)  R is a complex term which is defined as follows:
(18)
R
=
c
-
m
h
m
+
1
m
h
c
m
-
(
1
/
2
)
/
m
e
-
(
1
/
2
)
h
/
c
m
S
,
(19)
S
=
W
t
1
h
c
-
m
m
2
+
2
h
c
-
m
m
2
+
h
c
-
m
m
+
W
t
2
h
c
-
2
m
4
m
2
+
4
m
+
1
×
m
+
1
2
m
+
1
3
m
+
1
-
1
.The elementsW
t
1 and W
t
2 in (19) are the Whittaker functions M
a
1
,
b
x
(
z
) and M
a
2
,
b
x
(
z
) which are defined in terms of Kummer’s confluent hypergeometric function M (also known as the 1
F
1 function) [39] as follows:
(20)
W
t
1
=
M
a
1
,
b
x
(
z
)
=
e
-
z
/
2
z
1
/
2
+
b
x
M
1
2
+
b
x
-
a
1
,
1
+
2
b
x
,
z
,
W
t
2
=
M
a
2
,
b
x
(
z
)
=
e
-
z
/
2
z
1
/
2
+
b
x
M
1
2
+
b
x
-
a
2
,
1
+
2
b
x
,
z
,
where
(21)
a
1
=
(
m
+
1
)
m
-
(
1
/
2
)
(
2
m
+
1
)
m
,
a
2
=
(
m
+
1
)
m
-
(
1
/
2
)
(
2
m
+
1
)
m
+
1
,
b
x
=
(
1
/
2
)
(
2
m
+
1
)
m
+
1
2
,
z
=
h
c
m
.The hypergeometric functions from (20) were computed with the MATLAB function hypergeom [40, 41] as follows:
(22)
M
1
2
+
b
x
-
a
1
,
1
+
2
b
x
,
z
=
hypergeom
1
2
+
b
x
-
a
1
,
1
+
2
b
x
,
z
,
M
(
1
2
+
b
x
-
a
2
,
1
+
2
b
x
,
z
)
=
hypergeom
1
2
+
b
x
-
a
2
,
1
+
2
b
x
,
z
.


### 3.4. Estimating the 3D Visualization ofZ

Although then and k values defined within the intervals presented in Table 3 would produce complying α and β values for the statistical restrictions α, β
≤
0.20, some values may actually produce undefined elements in Z. Also, the suitable interval for h has yet to be identified. Thus, a second set of intervals within those presented in Table 3 must be identified and this is the objective of the 3D visualization of Z.In Figure3 the structure and pseudocode of the 3D plotting algorithm is presented. Initially the cost and time parameters for Z must be provided in addition to the type and associated parameters of the failure probability distribution (F
(
h
)) and τ. Then sets of evaluation vectors n
e, h
e, and k
e are created to specify the values of n, h, and k to be evaluated in Z. These values must be within the associated intervals presented in Table 3. For n the first element of n
e is n
l which is the lower bound for n while the last element is n
u which is the upper bound for n. As presented in Table 3 these bounds are dependent on δ.Figure 3
Algorithm for 3D visualization.Each set ofn, h, and k values are defined in such a way that vectors n
e, h
e, and k
e have the same number of elements. Then each set of n, h, and k values is evaluated in Z and the 3D matrix Cost is created where Cost
(
i
,
j
,
l
) contains the Z value for the ith element of n
e, the jth element of h
e, and the lth element of k
e. Then this 3D space is divided into gaps to visualize the patterns of the cost magnitudes in Cost. In general five gaps were considered and these were identified with a specific color. The colored points of these gaps were plotted into the 3D space by using the plot
3 command in MATLAB.In Figure4 the 3D visualization of Z with δ = 0.25 and exponential distribution for F
(
h
) and τ are presented. In this case Gap1 concentrates the minimum Cost values within 1.0% of the distance between the total minimum and maximum costs in Z. This gap is delimited by the following intervals: n
∈
(
73,113
), h
∈
(
1.59,5.0
), and k
∈
(
1.2816,1.805
).Example of 3D visualization forZ with δ
=
0.25 and exponential failure distribution with λ
=
0.25.
(a)
2D visualization
(b)
3D visualizationSomething very important to mention about this technique is that the total minimum found in this process is not the optimal (or near optimal) minimum of the wholeZ function. This is because n
e, h
e, and k
e do not contain all continuous values within the search space given the condition that n
e, h
e, and k
e must have the same size. Hence just some points of the search space are actually plotted.In this case the 3D visualization serves as a guide to identify the regions where the optimal parameters are likely to be found by a more specific search algorithm instead of searching across the whole solution space. Hence, the 3D visualization can support the explorative performance of the search process which then can be focused on the exploitative search.For more specific values forα (where k is given by Table 2 and (12)) and β the minimum values may be present in a different gap. In such case Figures 4(a) and 4(b) can be used in the following way to obtain an estimate for the parameters’ intervals:(i)
by having a specifick value for a required α a lower and upper bound for h can be obtained;
(ii)
then, an interval forn can be estimated based on the region delimited by the associated bounds of k and h. Each gap’s color can be used to guide the delimitation of the region.In Section4 sets of 3D visualizations with different settings of δ and failure distributions are presented. Then in Section 5 examples about how to use these 3D visualizations to improve estimation of parameters are presented and discussed.


## 3.1. Finding Optimal Intervals forn, k, and δ

An initial link betweenn, k, α, and β is defined by (4) and (5) where β
≤
β
^ and α
≤
α
^. β
^ and α
^ are maximum error levels and commonly β
^
≤
0.20 and α
^
≤
0.20. At this point there is no influence from other elements or costs from Z.The first step is to find an optimal interval forn. For this purpose the approximation presented in [37] for n with a specific value for β is used:
(9)
1
-
β
^
=
ϕ
(
δ
n
-
k
)
,
(10)
ϕ
-
1
(
1
-
β
^
)
=
δ
n
-
k
,
(11)
n
=
ϕ
-
1
(
1
-
β
^
)
+
k
δ
2
.Whilen depends on k and β, α is dependent only on k. From (4)  k is expressed as
(12)
k
=
ϕ
-
1
(
1
-
α
2
)
.By replacing (12) into (11) the estimation of n now depends on α and β. For any α and β the interval for n can be obtained with the following expression:
(13)
n
=
ϕ
-
1
(
1
-
β
)
+
ϕ
-
1
(
1
-
α
/
2
)
δ
2
.The optimal interval forn now depends on required intervals for α and β. In this case the interval from 0.01 (minimum) to 0.20 (maximum) was considered for both error probabilities (e.g., α and β
≤
0.20).The second step is to estimate an optimal interval fork. This is performed directly from (13) as defined by (12). In Table 2 the k values estimated for α
≤
0.20 are presented where an inverse correlation between both parameters is observed.Table 2
k values associated with α
≤
0.20 assuming normality for the X
- control char.
α
k
α
k
0.0020
3.0902
0.1040
1.6258
0.0080
2.6521
0.1100
1.5982
0.0140
2.4573
0.1160
1.5718
0.0200
2.3263
0.1220
1.5464
0.0260
2.2262
0.1280
1.5220
0.0320
2.1444
0.1340
1.4985
0.0380
2.0749
0.1400
1.4758
0.0440
2.0141
0.1460
1.4538
0.0500
1.9600
0.1520
1.4325
0.0560
1.9110
0.1580
1.4118
0.0620
1.8663
0.1640
1.3917
0.0680
1.8250
0.1700
1.3722
0.0740
1.7866
0.1760
1.3532
0.0800
1.7507
0.1820
1.3346
0.0860
1.7169
0.1880
1.3165
0.0920
1.6849
0.1940
1.2988
0.0980
1.6546
0.2000
1.2816In Figure1 the 3D visualizations of (13) considering the optimal intervals for k and β and three δ levels are presented.3D visualization ofn considering different values for k, β, and δ.
(a)
δ
=
0.10
(b)
δ
=
0.50
(c)
δ
=
0.90Forα
=
(
0.00,0.20
) from (12) it is estimated that k
=
(
1.2816,3.0902
) (inverse interval). While the overall 3D pattern does not change through Figures 1(a)–1(c), the difference is mainly observed on the magnitude of n which is highly dependent on δ. In Figure 2 the 2D visualization of the intervals (minimum and maximum values) for n considering different values for δ across β and α
≤
0.20 is presented. The intervals presented in Figures 1 and 2 can be used for any cost function for X
¯ control charts assuming normality.Figure 2
2D visualization of the intervals forn considering different values for δ.Although in practiceδ is considered to be equal to 0.5, the overall effect of δ on n has not been presented in graphical form as in Figure 2. For this work three scenarios are considered for δ based on the levels defined in [2, 4, 31, 38]: δ
=
(
0.25,0.50,1.00
). This leads to the following optimal intervals for n based on data presented in Figure 2: n = (73–570, 19–143, 5–36). In Table 3 the associated intervals for n, k, and δ are summarized.Table 3
Intervals defined forn, k, and δ assuming normality for the X
- control chart.
Parameter
Lower bound
Upper bound
n (δ
=
0.25)
73
570
n (δ
=
0.50)
19
143
n (δ
=
1.00)
5
36
k
1.2816
3.0902


## 3.2. Setting the Constant Elements ofZ

The cost function defined by (2) and (3) has many elements besides the n, h, k, and F
(
h
) parameters which were described in Table 1. In the literature this cost function has been studied with a set of standard values for the constant elements a, b, Y, W, Z
0, Z
1, D
0, and D
1 for specific δ values [2, 4, 31, 38]. In Table 4 these sets are presented for each level of δ.Table 4
Standard costs and times forZ [2].
δ
Z
0 (hours)
Z
1 (hours)
D
0 ($)
D
1 ($)
W ($)
Y ($)
a ($)
b ($)
0.25
0.025
0.10
25.00
475.00
550.00
250.00
10.00
2.11
0.50
0.250
1.00
50.00
950.00
1100.00
500.00
20.00
4.22
1.00
0.500
10.00
100.00
1900.00
2200.00
1000.00
40.00
8.44


## 3.3. Estimating the Probability Distributions andτ

F
(
h
) for the exponential and Weibull density functions is obtained from (6) as follows:
(14)
F
E
(
h
)
=
∫
0
h
‍
f
E
(
t
)
d
t
=
1
-
e
-
λ
h
.
F
W
(
h
)
=
∫
0
h
‍
f
W
(
t
)
d
t
=
-
c
-
m
h
c
-
m
h
m
(
-
1
+
e
-
(
h
/
c
)
m
)
.An important element of (2) is τ which is defined as the point of time within h where the process changes from the in-control state to the out-of-control state. This element is estimated as
(15)
τ
=
∫
0
h
‍
t
f
(
t
)
d
t
∫
0
h
‍
f
(
t
)
d
t
=
∫
0
h
‍
t
f
(
t
)
d
t
F
(
h
)
.Considering the exponential distributionτ is expressed as
(16)
τ
E
=
∫
0
h
‍
t
f
E
(
t
)
d
t
F
E
(
h
)
=
-
(
-
1
+
e
-
λ
h
+
e
-
λ
h
λ
h
)
/
λ
1
-
e
-
λ
h
.Equation (16) is widely presented in the literature [5, 11, 26, 30]; however for the Weibull distribution an expression for τ
W for any m and c has not been found. In this work an expression for τ
W was obtained with support of the symbolic toolbox of MATLAB and additional algebraic reductions:
(17)
τ
W
=
∫
0
h
‍
t
f
W
(
t
)
d
t
F
W
(
h
)
=
R
-
c
-
m
h
/
c
-
m
h
m
(
-
1
+
e
-
(
h
/
c
)
m
)
.In (17)  R is a complex term which is defined as follows:
(18)
R
=
c
-
m
h
m
+
1
m
h
c
m
-
(
1
/
2
)
/
m
e
-
(
1
/
2
)
h
/
c
m
S
,
(19)
S
=
W
t
1
h
c
-
m
m
2
+
2
h
c
-
m
m
2
+
h
c
-
m
m
+
W
t
2
h
c
-
2
m
4
m
2
+
4
m
+
1
×
m
+
1
2
m
+
1
3
m
+
1
-
1
.The elementsW
t
1 and W
t
2 in (19) are the Whittaker functions M
a
1
,
b
x
(
z
) and M
a
2
,
b
x
(
z
) which are defined in terms of Kummer’s confluent hypergeometric function M (also known as the 1
F
1 function) [39] as follows:
(20)
W
t
1
=
M
a
1
,
b
x
(
z
)
=
e
-
z
/
2
z
1
/
2
+
b
x
M
1
2
+
b
x
-
a
1
,
1
+
2
b
x
,
z
,
W
t
2
=
M
a
2
,
b
x
(
z
)
=
e
-
z
/
2
z
1
/
2
+
b
x
M
1
2
+
b
x
-
a
2
,
1
+
2
b
x
,
z
,
where
(21)
a
1
=
(
m
+
1
)
m
-
(
1
/
2
)
(
2
m
+
1
)
m
,
a
2
=
(
m
+
1
)
m
-
(
1
/
2
)
(
2
m
+
1
)
m
+
1
,
b
x
=
(
1
/
2
)
(
2
m
+
1
)
m
+
1
2
,
z
=
h
c
m
.The hypergeometric functions from (20) were computed with the MATLAB function hypergeom [40, 41] as follows:
(22)
M
1
2
+
b
x
-
a
1
,
1
+
2
b
x
,
z
=
hypergeom
1
2
+
b
x
-
a
1
,
1
+
2
b
x
,
z
,
M
(
1
2
+
b
x
-
a
2
,
1
+
2
b
x
,
z
)
=
hypergeom
1
2
+
b
x
-
a
2
,
1
+
2
b
x
,
z
.


## 3.4. Estimating the 3D Visualization ofZ

Although then and k values defined within the intervals presented in Table 3 would produce complying α and β values for the statistical restrictions α, β
≤
0.20, some values may actually produce undefined elements in Z. Also, the suitable interval for h has yet to be identified. Thus, a second set of intervals within those presented in Table 3 must be identified and this is the objective of the 3D visualization of Z.In Figure3 the structure and pseudocode of the 3D plotting algorithm is presented. Initially the cost and time parameters for Z must be provided in addition to the type and associated parameters of the failure probability distribution (F
(
h
)) and τ. Then sets of evaluation vectors n
e, h
e, and k
e are created to specify the values of n, h, and k to be evaluated in Z. These values must be within the associated intervals presented in Table 3. For n the first element of n
e is n
l which is the lower bound for n while the last element is n
u which is the upper bound for n. As presented in Table 3 these bounds are dependent on δ.Figure 3
Algorithm for 3D visualization.Each set ofn, h, and k values are defined in such a way that vectors n
e, h
e, and k
e have the same number of elements. Then each set of n, h, and k values is evaluated in Z and the 3D matrix Cost is created where Cost
(
i
,
j
,
l
) contains the Z value for the ith element of n
e, the jth element of h
e, and the lth element of k
e. Then this 3D space is divided into gaps to visualize the patterns of the cost magnitudes in Cost. In general five gaps were considered and these were identified with a specific color. The colored points of these gaps were plotted into the 3D space by using the plot
3 command in MATLAB.In Figure4 the 3D visualization of Z with δ = 0.25 and exponential distribution for F
(
h
) and τ are presented. In this case Gap1 concentrates the minimum Cost values within 1.0% of the distance between the total minimum and maximum costs in Z. This gap is delimited by the following intervals: n
∈
(
73,113
), h
∈
(
1.59,5.0
), and k
∈
(
1.2816,1.805
).Example of 3D visualization forZ with δ
=
0.25 and exponential failure distribution with λ
=
0.25.
(a)
2D visualization
(b)
3D visualizationSomething very important to mention about this technique is that the total minimum found in this process is not the optimal (or near optimal) minimum of the wholeZ function. This is because n
e, h
e, and k
e do not contain all continuous values within the search space given the condition that n
e, h
e, and k
e must have the same size. Hence just some points of the search space are actually plotted.In this case the 3D visualization serves as a guide to identify the regions where the optimal parameters are likely to be found by a more specific search algorithm instead of searching across the whole solution space. Hence, the 3D visualization can support the explorative performance of the search process which then can be focused on the exploitative search.For more specific values forα (where k is given by Table 2 and (12)) and β the minimum values may be present in a different gap. In such case Figures 4(a) and 4(b) can be used in the following way to obtain an estimate for the parameters’ intervals:(i)
by having a specifick value for a required α a lower and upper bound for h can be obtained;
(ii)
then, an interval forn can be estimated based on the region delimited by the associated bounds of k and h. Each gap’s color can be used to guide the delimitation of the region.In Section4 sets of 3D visualizations with different settings of δ and failure distributions are presented. Then in Section 5 examples about how to use these 3D visualizations to improve estimation of parameters are presented and discussed.


## 4. Libraries of 3D Visualizations

In Table5 the parameters considered for the 3D visualizations with different failure distributions are presented. In Figure 5 the selected 3D visualizations for the exponential distribution are presented while in Figures 6, 7, and 8 the visualizations for the Weibull distribution are presented.Table 5
Standard parameters forF
(
h
).
F
(
h
)
Level 1
Level 2
Level 3
Exponential
λ
=
0.05
λ
=
0.25
λ
=
0.50
Weibull
m
=
0.5, 1
/
c
=
0.05
m
=
0.5, 1
/
c
=
0.25
m
=
0.5, 1
/
c
=
0.50
Weibull
m
=
2.0, 1
/
c
=
0.05
m
=
2.0, 1
/
c
=
0.25
m
=
2.0, 1
/
c
=
0.50
Weibull
m
=
3.0, 1
/
c
=
0.05
m
=
3.0, 1
/
c
=
0.25
m
=
3.0, 1
/
c
=
0.50Library of 3D visualizations forZ with exponential failure distribution: λ
=
0.05, λ
=
0.25, and λ
=
0.50.
(a)
δ
=
0.25, λ
=
0.05
(b)
δ
=
0.25, λ
=
0.25
(c)
δ
=
0.25, λ
=
0.50
(d)
δ
=
0.50, λ
=
0.05
(e)
δ
=
0.50, λ
=
0.25
(f)
δ
=
0.50, λ
=
0.50
(g)
δ
=
1.00, λ
=
0.05
(h)
δ
=
1.00, λ
=
0.25
(i)
δ
=
1.00, λ
=
0.50Library of 3D visualizations forZ with δ
=
0.25 and Weibull failure distribution: 1
/
c
=
0.05, 1
/
c
=
0.25, and 1
/
c
=
0.50.
(a)
m
=
0.5, 1
/
c
=
0.05
(b)
m
=
0.5, 1
/
c
=
0.25
(c)
m
=
0.5, 1
/
c
=
0.50
(d)
m
=
2.0, 1
/
c
=
0.05
(e)
m
=
2.0, 1
/
c
=
0.25
(f)
m
=
2.0, 1
/
c
=
0.50
(g)
m
=
3.0, 1
/
c
=
0.05
(h)
m
=
3.0, 1
/
c
=
0.25
(i)
m
=
3.0, 1
/
c
=
0.50Library of 3D visualizations forZ with δ
=
0.50 and Weibull failure distribution: 1
/
c
=
0.05, 1
/
c
=
0.25, and 1
/
c
=
0.50.
(a)
m
=
0.5, 1
/
c
=
0.05
(b)
m
=
0.5, 1
/
c
=
0.25
(c)
m
=
0.5, 1
/
c
=
0.50
(d)
m
=
2.0, 1
/
c
=
0.05
(e)
m
=
2.0, 1
/
c
=
0.25
(f)
m
=
2.0, 1
/
c
=
0.50
(g)
m
=
3.0, 1
/
c
=
0.05
(h)
m
=
3.0, 1
/
c
=
0.25
(i)
m
=
3.0, 1
/
c
=
0.50Library of 3D visualizations forZ with δ
=
1.00 and Weibull failure distribution: 1
/
c
=
0.05, 1
/
c
=
0.25, and 1
/
c
=
0.50.
(a)
m
=
0.5, 1
/
c
=
0.05
(b)
m
=
0.5, 1
/
c
=
0.25
(c)
m
=
0.5, 1
/
c
=
0.50
(d)
m
=
2.0, 1
/
c
=
0.05
(e)
m
=
2.0, 1
/
c
=
0.25
(f)
m
=
2.0, 1
/
c
=
0.50
(g)
m
=
3.0, 1
/
c
=
0.05
(h)
m
=
3.0, 1
/
c
=
0.25
(i)
m
=
3.0, 1
/
c
=
0.50


## 5. Experiments and Results

Experiments were performed in order to evaluate the performance of the 3D visualizations to improve estimation of the control chart’s parameters. For these experimentsβ
^
≤
0.20 and α
^
≤
0.15 were considered. A direct search method across the regions identified with the 3D visualizations was performed to evaluate the performance of this technique. Hooke and Jeeves (HJ) and Genetic Algorithms (GA) were considered for comparison purposes. For this, the schemes presented in [12] for HJ and in [6] for GA were used.


### 5.1. Example of Estimation Process

In this case the ESD of anX
¯ control chart with δ
=
0.25 and Weibull failure distribution is considered. The Weibull distribution has the following parameters: m
=
0.50 and 1
/
c
=
0.25. The 3D model presented in Figure 6(b) is used for this case. As presented in Figure 9 two main views (front and left side) of the 3D model are used to estimate the intervals for each parameter.Figure 9
Example of estimation of intervals for each parameter using the 3D visualization.The first step is to identify the gap where the minimumZ is likely to be found according to α
^
≤
0.15. For this purpose the lower bound for k is estimated from the restriction α
^
≤
0.15 which from (12) implies k
≥
1.4395. As observed in the front view (
1
) of the 3D model the region delimited by Gap1 is within the space defined by k
≥
1.4395. Hence Gap1 is the region where the overall minimum of Z is likely to be found for this example.From the same front view(
1
) the upper bound for k is estimated and the interval k
∈
(
1.4395,2.115
) is identified. Using the same front view the lower and upper bound for h are estimated and the interval h
∈
(
2.5,5.0
) is identified. Finally with the left-side view (
2
) the lower and upper bound for n are estimated and the interval n
∈
(
73,193
) is identified.From this point an estimate ofn, h, and k can be given in terms of intervals. However to provide a specific value for Z and the set of parameters a search method must be performed. In this case a direct search method over all possible values within the (more specific) region delimited by these intervals can be applied leading to a more exploitative search.


### 5.2.  Results for Different Cases

HJ and GA were used to find the optimal parameters for each case presented in Figures5, 6, 7, and 8. For the 3D visualization method (3D Vis) the intervals identified for each case are presented in Table 6. The parameters estimated with HJ, GA, and 3D Vis are presented in Tables 7, 8, 9, and 10 for the exponential and Weibull distributions.Table 6
Gap1 intervals consideringα
^
≤
0.15.
F
(
h
)
δ
Parameters
n
l
n
u
h
l
h
u
k
l
k
u
λ
=
0.05
73
153
2.08
5.00
1.43
2.19
Exponential
0.25
λ
=
0.25
73
113
1.46
5.00
1.43
1.81
λ
=
0.50
73
153
1.67
5.00
1.43
1.88
λ
=
0.05
19
39
1.67
4.38
1.43
2.26
Exponential
0.50
λ
=
0.25
19
29
1.04
1.88
1.43
1.88
λ
=
0.50
19
23
0.83
1.46
1.43
1.51
λ
=
0.05
8
12
1.29
1.94
1.81
2.68
Exponential
1.00
λ
=
0.25
9
15
0.64
0.80
1.81
2.45
λ
=
0.50
9
15
0.48
0.64
1.81
2.45
Weibull
m
=
0.50
1
/
c
=
0.05
73
133
2.92
5.00
1.44
2.19
0.25
1
/
c
=
0.25
73
193
2.50
5.00
1.44
2.12
1
/
c
=
0.50
73
173
2.72
5.00
1.44
2.04
Weibull
m
=
0.50
1
/
c
=
0.05
19
64
3.13
5.00
1.44
2.57
0.50
1
/
c
=
0.25
24
54
2.29
5.00
1.44
1.96
1
/
c
=
0.50
24
44
1.88
3.96
1.44
1.66
Weibull
m
=
0.50
1
/
c
=
0.05
9
20
1.94
5.00
1.69
2.74
1.00
1
/
c
=
0.25
8
12
0.48
0.64
1.98
2.39
1
/
c
=
0.50
8
12
0.27
0.48
1.87
2.39
Weibull
m
=
2.00
1
/
c
=
0.05
73
113
1.67
2.92
1.44
1.66
0.25
1
/
c
=
0.25
73
93
1.25
1.67
1.44
1.51
1
/
c
=
0.50
73
113
3.33
5.00
1.44
1.81
Weibull
m
=
2.00
1
/
c
=
0.05
19
34
1.04
1.67
1.44
2.04
0.50
1
/
c
=
0.25
19
29
0.63
0.84
1.44
1.81
1
/
c
=
0.50
19
29
0.42
0.84
1.44
1.81
Weibull
m
=
2.00
1
/
c
=
0.05
8
13
0.97
1.61
1.81
2.80
1.00
1
/
c
=
0.25
8
20
0.81
1.29
1.57
2.74
1
/
c
=
0.50
8
20
0.64
0.97
1.52
2.68
Weibull
m
=
3.00
1
/
c
=
0.05
73
93
1.04
1.88
1.44
1.51
0.25
1
/
c
=
0.25
73
85
0.83
1.25
1.44
1.45
1
/
c
=
0.50
73
85
0.63
1.04
1.44
1.45
Weibull
m
=
3.00
1
/
c
=
0.05
19
34
1.04
1.25
1.44
2.04
0.50
1
/
c
=
0.25
19
29
0.63
0.84
1.44
1.81
1
/
c
=
0.50
19
29
0.42
0.84
1.44
1.81
Weibull
m
=
3.00
1
/
c
=
0.05
8
20
0.97
1.61
1.69
2.80
1.00
1
/
c
=
0.25
8
23
0.81
1.29
1.44
2.80
1
/
c
=
0.50
8
23
0.81
1.13
1.44
2.86Table 7
Accuracies of the 3D Vis method in comparison with HJ and GA: exponential failure distribution.
Method
δ
λ
n
h
k
α
β
Z
0.05
84
3.9674
1.4442
0.1487
0.1984
159.84
0.25
0.25
84
2.7926
1.4471
0.1479
0.1992
337.63
0.50
87
3.8406
1.4855
0.1374
0.1986
448.12
0.05
27
2.7447
1.5725
0.1158
0.1525
229.47
HJ
0.50
0.25
24
1.4118
1.4728
0.1408
0.1643
464.81
0.50
22
1.1061
1.4715
0.1412
0.1911
609.54
0.05
12
1.8757
2.1508
0.0315
0.0945
262.03
1.00
0.25
12
0.8472
2.1153
0.0344
0.0887
304.87
0.50
11
0.5785
2.0766
0.0378
0.1075
300.81
0.05
89
3.8265
1.5060
0.1321
0.1969
161.97
0.25
0.25
93
2.9533
1.5691
0.1166
0.1999
343.15
0.50
86
2.8405
1.4692
0.1418
0.1978
450.08
0.05
27
2.8190
1.5358
0.1246
0.1440
229.50
GA
0.50
0.25
24
1.4427
1.4396
0.1500
0.1562
464.62
0.50
24
1.1698
1.4395
0.1500
0.1562
609.44
0.05
10
1.6132
2.2231
0.0262
0.1738
263.88
1.00
0.25
10
0.6060
2.2469
0.0246
0.1800
308.33
0.50
9
0.5188
2.0447
0.0409
0.1697
302.14
0.05
84
4.3500
1.4400
0.1499
0.1972
159.40
0.25
0.25
84
2.8000
1.4400
0.1499
0.1972
337.43
0.50
84
3.7400
1.4400
0.1499
0.1972
446.46
0.05
27
2.7700
1.5700
0.1164
0.1519
229.46
3D Vis
0.50
0.25
23
1.4100
1.4400
0.1499
0.1690
464.58
0.50
21
1.0900
1.4400
0.1499
0.1972
608.95
0.05
12
1.8500
2.1700
0.0300
0.0978
262.01
1.00
0.25
11
0.8000
2.1000
0.0357
0.1119
304.87
0.50
11
0.5700
2.0700
0.0385
0.1063
300.79Table 8
Accuracies of the 3D Vis method in comparison with HJ and GA: Weibull failure distribution withm
=
0
.
50.
Method
δ
1
/
c
n
h
k
α
β
Z
0.05
88
3.8947
1.4938
0.1352
0.1972
133.53
0.25
0.25
88
3.4561
1.4929
0.1355
0.1970
284.13
0.50
88
3.4257
1.5007
0.1334
0.1991
397.12
0.05
29
3.7472
1.6544
0.0980
0.1496
190.11
HJ
0.50
0.25
38
3.8785
1.4805
0.1387
0.0546
432.06
0.50
36
2.8071
1.5394
0.1237
0.0721
609.37
0.05
13
3.0000
2.1228
0.0338
0.0691
259.11
1.00
0.25
10
0.4732
2.1968
0.0280
0.1672
320.68
0.50
9
0.2767
2.1578
0.0309
0.1998
302.82
0.05
99
3.8213
1.6451
0.0999
0.1998
138.62
0.25
0.25
90
3.9306
1.5235
0.1276
0.1981
280.60
0.50
91
3.5717
1.5365
0.1244
0.1981
396.60
0.05
29
3.9993
1.5870
0.1125
0.1344
188.41
GA
0.50
0.25
33
3.2033
1.6592
0.0971
0.1125
436.16
0.50
25
2.4127
1.4527
0.1463
0.1474
611.79
0.05
10
2.4603
2.1888
0.0286
0.1651
263.23
1.00
0.25
8
0.5002
1.9841
0.0472
0.1992
321.98
0.50
8
0.5100
1.9667
0.0492
0.1944
314.96
0.05
85
5.0000
1.4525
0.1464
0.1969
123.00
0.25
0.25
95
5.0000
1.4395
0.1500
0.1593
273.06
0.50
86
4.9900
1.4395
0.1500
0.1896
383.55
0.05
32
4.9850
1.5795
0.1142
0.1058
184.04
3D Vis
0.50
0.25
36
3.5700
1.4395
0.1500
0.0593
431.36
0.50
33
2.7350
1.4395
0.1500
0.0760
607.52
0.05
13
3.0000
2.1300
0.0332
0.0700
259.11
1.00
0.25
10
0.5100
2.1700
0.0300
0.1605
320.53
0.50
9
0.2700
2.1400
0.0324
0.1949
302.71Table 9
Accuracies of the 3D Vis method in comparison with HJ and GA: Weibull failure distribution withm
=
2
.
00.
Method
δ
1
/
c
n
h
k
α
β
Z
0.05
91
2.0000
1.5368
0.1243
0.1982
229.99
0.25
0.25
87
1.1926
1.4809
0.1386
0.1973
400.34
0.50
88
3.4406
1.4962
0.1346
0.1979
516.35
0.05
28
1.3000
1.7938
0.0729
0.1971
270.59
HJ
0.50
0.25
28
0.7561
1.7929
0.0730
0.1969
444.33
0.50
26
0.6257
1.7007
0.0890
0.1980
541.21
0.05
11
1.3000
2.2544
0.0242
0.1441
263.23
1.00
0.25
13
1.0000
2.0805
0.0375
0.0636
288.67
0.50
12
0.8071
2.0394
0.0414
0.0771
288.84
0.05
90
2.0769
1.5194
0.1287
0.1970
229.31
0.25
0.25
99
1.2073
1.6432
0.1003
0.1992
416.15
0.50
93
3.8578
1.5666
0.1172
0.1992
514.82
0.05
24
1.4386
1.5782
0.1145
0.1918
270.00
GA
0.50
0.25
23
0.8823
1.5463
0.1220
0.1972
441.89
0.50
24
0.6575
1.4507
0.1469
0.1589
546.09
0.05
8
1.2560
1.9792
0.0478
0.1979
266.44
1.00
0.25
8
0.8372
1.9746
0.0483
0.1966
294.58
0.50
11
0.7471
1.9763
0.0481
0.1971
296.11
0.05
84
1.9900
1.4495
0.1472
0.1999
224.99
0.25
0.25
84
1.2500
1.4495
0.1472
0.1999
396.69
0.50
84
5.0000
1.4395
0.1500
0.1971
502.14
0.05
24
1.3400
1.5995
0.1097
0.1976
268.72
3D Vis
0.50
0.25
22
0.7750
1.4995
0.1337
0.1988
438.42
0.50
21
0.6200
1.4495
0.1472
0.1999
536.13
0.05
11
1.3400
2.2000
0.0278
0.1321
263.11
1.00
0.25
12
0.9667
2.0730
0.0382
0.0821
288.62
0.50
13
0.8333
2.0467
0.0407
0.0595
288.75Table 10
Accuracies of the 3D Vis method in comparison with HJ and GA: Weibull failure distribution withm
=
3
.
00.
Method
δ
1
/
c
n
h
k
α
β
Z
0.05
94
1.4708
1.5785
0.1145
0.1989
265.10
0.25
0.25
94
0.9899
1.5750
0.1153
0.1980
408.36
0.50
89
3.8030
1.5051
0.1323
0.1967
525.15
0.05
30
1.0628
1.8879
0.0590
0.1975
280.22
HJ
0.50
0.25
26
0.6986
1.6610
0.0967
0.1871
400.82
0.50
28
0.6228
1.8031
0.0714
0.1997
471.59
0.05
11
1.1775
2.2431
0.0249
0.1415
260.99
1.00
0.25
13
1.1000
1.9935
0.0462
0.0535
277.85
0.50
13
0.9540
1.9798
0.0477
0.0520
278.10
0.05
84
1.6472
1.4423
0.1492
0.1978
259.88
0.25
0.25
94
0.9069
1.5788
0.1144
0.1990
407.81
0.50
89
0.7631
1.5082
0.1315
0.1975
483.76
0.05
31
1.0235
1.9160
0.0554
0.1927
283.05
GA
0.50
0.25
30
0.6695
1.8805
0.0600
0.1954
406.84
0.50
24
0.6120
1.4459
0.1482
0.1578
478.75
0.05
8
1.1491
1.9850
0.0471
0.1995
263.71
1.00
0.25
8
0.8793
1.9745
0.0483
0.1966
285.30
0.50
8
0.7760
1.9754
0.0482
0.1968
287.46
0.05
84
1.4300
1.4495
0.1472
0.1999
255.19
0.25
0.25
84
0.9200
1.4395
0.1500
0.1971
392.14
0.50
84
0.7650
1.4395
0.1500
0.1971
474.57
0.05
25
1.0800
1.6495
0.0990
0.1975
275.71
3D Vis
0.50
0.25
23
0.6950
1.5495
0.1213
0.1981
396.94
0.50
23
0.5800
1.5495
0.1213
0.1981
463.77
0.05
11
1.2300
2.2100
0.0271
0.1342
260.82
1.00
0.25
13
1.0667
2.0395
0.0414
0.0587
277.70
0.50
13
0.9467
1.9595
0.0501
0.0499
278.07A factorial analysis was performed on the data presented in Tables7 to 10 to assess the effect of the proposed method on the task of minimizing Z. In Figure 10 the “Main Effects” plots on the mean of Z across the following factors and levels are presented.(i)
The failure probability distribution (F
(
h
)): exponential (Exp), Weibull (Weib) with m
=
0.5, m
=
2.0, and m
=
3.0.
(ii)
The magnitude of the change of the mean of the process to be detected (δ): 0.25, 0.50, and 1.00.
(iii)
The estimation method for the control chart’s parameters: HJ, GA, and 3D Vis.
(iv)
The failure rate of the distributions (determined byλ for the exponential distribution and 1
/
c for the Weibull distribution): 0.05, 0.25, 0.50.Figure 10
Main effects plots.The mean value ofZ changes according to these factors observing a positive correlation between Z and the failure rate given by λ and 1
/
c. For the estimation method, GA presents the higher Z with 347.705 while HJ presents an average of 346.050. In contrast, the parameters estimated with the 3D Vis method provided an average of 341.147. This difference was tested for statistical significance by means of a “2-sample t-test” performed with Minitab Ver. 16.2.3.0.In Figures11 and 12 the results of the statistical tests for the hyphoteses H
0
:
μ
HJ
-
μ
3
D
=
0 and H
0
:
μ
GA
-
μ
3
D
=
0 are presented. In both cases the paired t-test concluded that the results obtained with GA and HJ were significantly different from the results obtained with the proposed 3D method thus rejecting H
0.Figure 11
Pairedt-test for H
0
:
μ
HJ
-
μ
3
D
=
0.Figure 12
Pairedt-test for H
0
:
μ
GA
-
μ
3
D
=
0.Because the results presented in Tables7 to 10 for the 3D Vis method were obtained with a direct search within the specific regions delimited by Gap1 it can be assumed that these are optimal (or at least near optimal) results. As presented in Figures 11 and 12 these results were statistically better than those obtained with well established heuristics as HJ and GA. Thus, these results can be used as benchmark to evaluate the performance of other algorithms or methods.


## 5.1. Example of Estimation Process

In this case the ESD of anX
¯ control chart with δ
=
0.25 and Weibull failure distribution is considered. The Weibull distribution has the following parameters: m
=
0.50 and 1
/
c
=
0.25. The 3D model presented in Figure 6(b) is used for this case. As presented in Figure 9 two main views (front and left side) of the 3D model are used to estimate the intervals for each parameter.Figure 9
Example of estimation of intervals for each parameter using the 3D visualization.The first step is to identify the gap where the minimumZ is likely to be found according to α
^
≤
0.15. For this purpose the lower bound for k is estimated from the restriction α
^
≤
0.15 which from (12) implies k
≥
1.4395. As observed in the front view (
1
) of the 3D model the region delimited by Gap1 is within the space defined by k
≥
1.4395. Hence Gap1 is the region where the overall minimum of Z is likely to be found for this example.From the same front view(
1
) the upper bound for k is estimated and the interval k
∈
(
1.4395,2.115
) is identified. Using the same front view the lower and upper bound for h are estimated and the interval h
∈
(
2.5,5.0
) is identified. Finally with the left-side view (
2
) the lower and upper bound for n are estimated and the interval n
∈
(
73,193
) is identified.From this point an estimate ofn, h, and k can be given in terms of intervals. However to provide a specific value for Z and the set of parameters a search method must be performed. In this case a direct search method over all possible values within the (more specific) region delimited by these intervals can be applied leading to a more exploitative search.


## 5.2.  Results for Different Cases

HJ and GA were used to find the optimal parameters for each case presented in Figures5, 6, 7, and 8. For the 3D visualization method (3D Vis) the intervals identified for each case are presented in Table 6. The parameters estimated with HJ, GA, and 3D Vis are presented in Tables 7, 8, 9, and 10 for the exponential and Weibull distributions.Table 6
Gap1 intervals consideringα
^
≤
0.15.
F
(
h
)
δ
Parameters
n
l
n
u
h
l
h
u
k
l
k
u
λ
=
0.05
73
153
2.08
5.00
1.43
2.19
Exponential
0.25
λ
=
0.25
73
113
1.46
5.00
1.43
1.81
λ
=
0.50
73
153
1.67
5.00
1.43
1.88
λ
=
0.05
19
39
1.67
4.38
1.43
2.26
Exponential
0.50
λ
=
0.25
19
29
1.04
1.88
1.43
1.88
λ
=
0.50
19
23
0.83
1.46
1.43
1.51
λ
=
0.05
8
12
1.29
1.94
1.81
2.68
Exponential
1.00
λ
=
0.25
9
15
0.64
0.80
1.81
2.45
λ
=
0.50
9
15
0.48
0.64
1.81
2.45
Weibull
m
=
0.50
1
/
c
=
0.05
73
133
2.92
5.00
1.44
2.19
0.25
1
/
c
=
0.25
73
193
2.50
5.00
1.44
2.12
1
/
c
=
0.50
73
173
2.72
5.00
1.44
2.04
Weibull
m
=
0.50
1
/
c
=
0.05
19
64
3.13
5.00
1.44
2.57
0.50
1
/
c
=
0.25
24
54
2.29
5.00
1.44
1.96
1
/
c
=
0.50
24
44
1.88
3.96
1.44
1.66
Weibull
m
=
0.50
1
/
c
=
0.05
9
20
1.94
5.00
1.69
2.74
1.00
1
/
c
=
0.25
8
12
0.48
0.64
1.98
2.39
1
/
c
=
0.50
8
12
0.27
0.48
1.87
2.39
Weibull
m
=
2.00
1
/
c
=
0.05
73
113
1.67
2.92
1.44
1.66
0.25
1
/
c
=
0.25
73
93
1.25
1.67
1.44
1.51
1
/
c
=
0.50
73
113
3.33
5.00
1.44
1.81
Weibull
m
=
2.00
1
/
c
=
0.05
19
34
1.04
1.67
1.44
2.04
0.50
1
/
c
=
0.25
19
29
0.63
0.84
1.44
1.81
1
/
c
=
0.50
19
29
0.42
0.84
1.44
1.81
Weibull
m
=
2.00
1
/
c
=
0.05
8
13
0.97
1.61
1.81
2.80
1.00
1
/
c
=
0.25
8
20
0.81
1.29
1.57
2.74
1
/
c
=
0.50
8
20
0.64
0.97
1.52
2.68
Weibull
m
=
3.00
1
/
c
=
0.05
73
93
1.04
1.88
1.44
1.51
0.25
1
/
c
=
0.25
73
85
0.83
1.25
1.44
1.45
1
/
c
=
0.50
73
85
0.63
1.04
1.44
1.45
Weibull
m
=
3.00
1
/
c
=
0.05
19
34
1.04
1.25
1.44
2.04
0.50
1
/
c
=
0.25
19
29
0.63
0.84
1.44
1.81
1
/
c
=
0.50
19
29
0.42
0.84
1.44
1.81
Weibull
m
=
3.00
1
/
c
=
0.05
8
20
0.97
1.61
1.69
2.80
1.00
1
/
c
=
0.25
8
23
0.81
1.29
1.44
2.80
1
/
c
=
0.50
8
23
0.81
1.13
1.44
2.86Table 7
Accuracies of the 3D Vis method in comparison with HJ and GA: exponential failure distribution.
Method
δ
λ
n
h
k
α
β
Z
0.05
84
3.9674
1.4442
0.1487
0.1984
159.84
0.25
0.25
84
2.7926
1.4471
0.1479
0.1992
337.63
0.50
87
3.8406
1.4855
0.1374
0.1986
448.12
0.05
27
2.7447
1.5725
0.1158
0.1525
229.47
HJ
0.50
0.25
24
1.4118
1.4728
0.1408
0.1643
464.81
0.50
22
1.1061
1.4715
0.1412
0.1911
609.54
0.05
12
1.8757
2.1508
0.0315
0.0945
262.03
1.00
0.25
12
0.8472
2.1153
0.0344
0.0887
304.87
0.50
11
0.5785
2.0766
0.0378
0.1075
300.81
0.05
89
3.8265
1.5060
0.1321
0.1969
161.97
0.25
0.25
93
2.9533
1.5691
0.1166
0.1999
343.15
0.50
86
2.8405
1.4692
0.1418
0.1978
450.08
0.05
27
2.8190
1.5358
0.1246
0.1440
229.50
GA
0.50
0.25
24
1.4427
1.4396
0.1500
0.1562
464.62
0.50
24
1.1698
1.4395
0.1500
0.1562
609.44
0.05
10
1.6132
2.2231
0.0262
0.1738
263.88
1.00
0.25
10
0.6060
2.2469
0.0246
0.1800
308.33
0.50
9
0.5188
2.0447
0.0409
0.1697
302.14
0.05
84
4.3500
1.4400
0.1499
0.1972
159.40
0.25
0.25
84
2.8000
1.4400
0.1499
0.1972
337.43
0.50
84
3.7400
1.4400
0.1499
0.1972
446.46
0.05
27
2.7700
1.5700
0.1164
0.1519
229.46
3D Vis
0.50
0.25
23
1.4100
1.4400
0.1499
0.1690
464.58
0.50
21
1.0900
1.4400
0.1499
0.1972
608.95
0.05
12
1.8500
2.1700
0.0300
0.0978
262.01
1.00
0.25
11
0.8000
2.1000
0.0357
0.1119
304.87
0.50
11
0.5700
2.0700
0.0385
0.1063
300.79Table 8
Accuracies of the 3D Vis method in comparison with HJ and GA: Weibull failure distribution withm
=
0
.
50.
Method
δ
1
/
c
n
h
k
α
β
Z
0.05
88
3.8947
1.4938
0.1352
0.1972
133.53
0.25
0.25
88
3.4561
1.4929
0.1355
0.1970
284.13
0.50
88
3.4257
1.5007
0.1334
0.1991
397.12
0.05
29
3.7472
1.6544
0.0980
0.1496
190.11
HJ
0.50
0.25
38
3.8785
1.4805
0.1387
0.0546
432.06
0.50
36
2.8071
1.5394
0.1237
0.0721
609.37
0.05
13
3.0000
2.1228
0.0338
0.0691
259.11
1.00
0.25
10
0.4732
2.1968
0.0280
0.1672
320.68
0.50
9
0.2767
2.1578
0.0309
0.1998
302.82
0.05
99
3.8213
1.6451
0.0999
0.1998
138.62
0.25
0.25
90
3.9306
1.5235
0.1276
0.1981
280.60
0.50
91
3.5717
1.5365
0.1244
0.1981
396.60
0.05
29
3.9993
1.5870
0.1125
0.1344
188.41
GA
0.50
0.25
33
3.2033
1.6592
0.0971
0.1125
436.16
0.50
25
2.4127
1.4527
0.1463
0.1474
611.79
0.05
10
2.4603
2.1888
0.0286
0.1651
263.23
1.00
0.25
8
0.5002
1.9841
0.0472
0.1992
321.98
0.50
8
0.5100
1.9667
0.0492
0.1944
314.96
0.05
85
5.0000
1.4525
0.1464
0.1969
123.00
0.25
0.25
95
5.0000
1.4395
0.1500
0.1593
273.06
0.50
86
4.9900
1.4395
0.1500
0.1896
383.55
0.05
32
4.9850
1.5795
0.1142
0.1058
184.04
3D Vis
0.50
0.25
36
3.5700
1.4395
0.1500
0.0593
431.36
0.50
33
2.7350
1.4395
0.1500
0.0760
607.52
0.05
13
3.0000
2.1300
0.0332
0.0700
259.11
1.00
0.25
10
0.5100
2.1700
0.0300
0.1605
320.53
0.50
9
0.2700
2.1400
0.0324
0.1949
302.71Table 9
Accuracies of the 3D Vis method in comparison with HJ and GA: Weibull failure distribution withm
=
2
.
00.
Method
δ
1
/
c
n
h
k
α
β
Z
0.05
91
2.0000
1.5368
0.1243
0.1982
229.99
0.25
0.25
87
1.1926
1.4809
0.1386
0.1973
400.34
0.50
88
3.4406
1.4962
0.1346
0.1979
516.35
0.05
28
1.3000
1.7938
0.0729
0.1971
270.59
HJ
0.50
0.25
28
0.7561
1.7929
0.0730
0.1969
444.33
0.50
26
0.6257
1.7007
0.0890
0.1980
541.21
0.05
11
1.3000
2.2544
0.0242
0.1441
263.23
1.00
0.25
13
1.0000
2.0805
0.0375
0.0636
288.67
0.50
12
0.8071
2.0394
0.0414
0.0771
288.84
0.05
90
2.0769
1.5194
0.1287
0.1970
229.31
0.25
0.25
99
1.2073
1.6432
0.1003
0.1992
416.15
0.50
93
3.8578
1.5666
0.1172
0.1992
514.82
0.05
24
1.4386
1.5782
0.1145
0.1918
270.00
GA
0.50
0.25
23
0.8823
1.5463
0.1220
0.1972
441.89
0.50
24
0.6575
1.4507
0.1469
0.1589
546.09
0.05
8
1.2560
1.9792
0.0478
0.1979
266.44
1.00
0.25
8
0.8372
1.9746
0.0483
0.1966
294.58
0.50
11
0.7471
1.9763
0.0481
0.1971
296.11
0.05
84
1.9900
1.4495
0.1472
0.1999
224.99
0.25
0.25
84
1.2500
1.4495
0.1472
0.1999
396.69
0.50
84
5.0000
1.4395
0.1500
0.1971
502.14
0.05
24
1.3400
1.5995
0.1097
0.1976
268.72
3D Vis
0.50
0.25
22
0.7750
1.4995
0.1337
0.1988
438.42
0.50
21
0.6200
1.4495
0.1472
0.1999
536.13
0.05
11
1.3400
2.2000
0.0278
0.1321
263.11
1.00
0.25
12
0.9667
2.0730
0.0382
0.0821
288.62
0.50
13
0.8333
2.0467
0.0407
0.0595
288.75Table 10
Accuracies of the 3D Vis method in comparison with HJ and GA: Weibull failure distribution withm
=
3
.
00.
Method
δ
1
/
c
n
h
k
α
β
Z
0.05
94
1.4708
1.5785
0.1145
0.1989
265.10
0.25
0.25
94
0.9899
1.5750
0.1153
0.1980
408.36
0.50
89
3.8030
1.5051
0.1323
0.1967
525.15
0.05
30
1.0628
1.8879
0.0590
0.1975
280.22
HJ
0.50
0.25
26
0.6986
1.6610
0.0967
0.1871
400.82
0.50
28
0.6228
1.8031
0.0714
0.1997
471.59
0.05
11
1.1775
2.2431
0.0249
0.1415
260.99
1.00
0.25
13
1.1000
1.9935
0.0462
0.0535
277.85
0.50
13
0.9540
1.9798
0.0477
0.0520
278.10
0.05
84
1.6472
1.4423
0.1492
0.1978
259.88
0.25
0.25
94
0.9069
1.5788
0.1144
0.1990
407.81
0.50
89
0.7631
1.5082
0.1315
0.1975
483.76
0.05
31
1.0235
1.9160
0.0554
0.1927
283.05
GA
0.50
0.25
30
0.6695
1.8805
0.0600
0.1954
406.84
0.50
24
0.6120
1.4459
0.1482
0.1578
478.75
0.05
8
1.1491
1.9850
0.0471
0.1995
263.71
1.00
0.25
8
0.8793
1.9745
0.0483
0.1966
285.30
0.50
8
0.7760
1.9754
0.0482
0.1968
287.46
0.05
84
1.4300
1.4495
0.1472
0.1999
255.19
0.25
0.25
84
0.9200
1.4395
0.1500
0.1971
392.14
0.50
84
0.7650
1.4395
0.1500
0.1971
474.57
0.05
25
1.0800
1.6495
0.0990
0.1975
275.71
3D Vis
0.50
0.25
23
0.6950
1.5495
0.1213
0.1981
396.94
0.50
23
0.5800
1.5495
0.1213
0.1981
463.77
0.05
11
1.2300
2.2100
0.0271
0.1342
260.82
1.00
0.25
13
1.0667
2.0395
0.0414
0.0587
277.70
0.50
13
0.9467
1.9595
0.0501
0.0499
278.07A factorial analysis was performed on the data presented in Tables7 to 10 to assess the effect of the proposed method on the task of minimizing Z. In Figure 10 the “Main Effects” plots on the mean of Z across the following factors and levels are presented.(i)
The failure probability distribution (F
(
h
)): exponential (Exp), Weibull (Weib) with m
=
0.5, m
=
2.0, and m
=
3.0.
(ii)
The magnitude of the change of the mean of the process to be detected (δ): 0.25, 0.50, and 1.00.
(iii)
The estimation method for the control chart’s parameters: HJ, GA, and 3D Vis.
(iv)
The failure rate of the distributions (determined byλ for the exponential distribution and 1
/
c for the Weibull distribution): 0.05, 0.25, 0.50.Figure 10
Main effects plots.The mean value ofZ changes according to these factors observing a positive correlation between Z and the failure rate given by λ and 1
/
c. For the estimation method, GA presents the higher Z with 347.705 while HJ presents an average of 346.050. In contrast, the parameters estimated with the 3D Vis method provided an average of 341.147. This difference was tested for statistical significance by means of a “2-sample t-test” performed with Minitab Ver. 16.2.3.0.In Figures11 and 12 the results of the statistical tests for the hyphoteses H
0
:
μ
HJ
-
μ
3
D
=
0 and H
0
:
μ
GA
-
μ
3
D
=
0 are presented. In both cases the paired t-test concluded that the results obtained with GA and HJ were significantly different from the results obtained with the proposed 3D method thus rejecting H
0.Figure 11
Pairedt-test for H
0
:
μ
HJ
-
μ
3
D
=
0.Figure 12
Pairedt-test for H
0
:
μ
GA
-
μ
3
D
=
0.Because the results presented in Tables7 to 10 for the 3D Vis method were obtained with a direct search within the specific regions delimited by Gap1 it can be assumed that these are optimal (or at least near optimal) results. As presented in Figures 11 and 12 these results were statistically better than those obtained with well established heuristics as HJ and GA. Thus, these results can be used as benchmark to evaluate the performance of other algorithms or methods.


## 6. Conclusions and Future Work

In this work an important cost function considered for the ESD ofX
¯ control charts was analyzed to provide benchmark (optimal) values for the estimation of the parameters n, h, and k. This was also performed to provide the means for evaluation of heuristic algorithms for the estimation of the chart’s parameters. Particularly the analyses presented for n and α-β-δ ((13), Figure 2) and k and α ((12), Table 2) are valid for the general ESD when normality is assumed regardless of the cost function considered.Among the applications of the present work the following can be mentioned.(i)
The 3D plotting algorithm presented in Figure3 can be adapted to visualize other multivariable cost functions and (18) and (19) can be used to explore the patterns of cost functions with different Weibull parameters.
(ii)
The libraries presented in Figures5 to 8 provide the means to further analyze and/or understand the effect of the failure distribution and associated parameters on the minimum and maximum values of the cost function Z.
(iii)
If alternative values for the constant parameters of the cost function (Table4) do not change the behavior of the 3D visualization (e.g., the gap’s regions) the libraries can be used to get a graphical and fast estimation of the chart’s parameters as described in Section 5.1.
(iv)
The intervals identified with the 3D visualizations can be used to guide the search strategies of heuristic algorithms.The concept of 3D visualization is important to understand the dynamics of different processes and thus it should be explored for other cases of ESD. Future directions of this work are the following.(i)
To perform the 3D visualizations of other cost functions and failure distributions for other control charts as EWMA [11, 30], p [13], and S
/
S
2 [42, 43].
(ii)
To estimate the valid intervals forn, h, and k considering nonnormality [2, 44] and also to analyze the effect of the assumption of nonnormality on the behavior of the optimal regions of the cost function.
(iii)
To perform the 3D visualization of cost functions with variable sampling intervals.
(iv)
To develop the mathematical foundations to analyze the effect of alternative constant parameters on the behavior of the cost function. In such case, the number of variables in the cost function will be increased and thus alternativeN-D visualization models have to be developed.
(v)
To develop the mathematical model of the gaps’ patterns observed in the 3D visualizations to reduce the complexity of the cost function.



---

*Source: 879456-2014-11-09.xml*