Applied Mathematics, 2011, 2, 1539-1545
doi:10.4236/am.2011.212219 Published Online December 2011 (http://www.SciRP.org/journal/am)
Minimax Multivariate Control Chart Using a
Polynomial Function
Johnson Ademola Adewara1, Kayode Samuel Adekeye2, Osebekwin Ebenezer Asiribo3,
Samuel Babatope Adejuyigbe4
1Distance Learning Institu te, University of Lagos, Akoka, Nigeria
2Department o f Mat hematical Sci en ces, Redeemer’s University, Redemption City, Nigeria
3Department of Stat i s t i c s, University of Agriculture, Abeoku ta, Nigeria
4Department of Mechanical Engineering, University of Agriculture, Abeokuta, Nigeria
Received September 15, 2011; revised November 4, 2011; accepted Nov ember 12, 2011
Abstract
Minimax control chart uses the joint probability distribution of the maximum and minimum standardized
sample means to obtain the control limits for monitoring purpose. However, the derivation of the joint prob-
ability distribution needed to obtain the minimax control limits is complex. In this paper the multivariate
normal distribution is integrated numerically using Simpson’s one third rule to obtain a non-linear polyno-
mial (NLP) function. This NLP function is then substituted and solved numerically using Newton Raphson
method to obtain the control limits for the minimax control chart. The approach helps to overcome the prob-
lem of obtaining the joint probability distribution needed for estimating the control limits of both the maxi-
mum and the minimum statistic for monitoring multivariate process.
Keywords: Minimax, Non-Linear Polynomial, Process, Maximum and Minimum
1. Introduction
Multivariate statistical process control (MSPC) is parti-
cularly important in the industries where data are col-
lected on more than one variable. In practice, most of the
quality characteristics to be controlled and monitored are
not independent. The reason is that most of the variables
involved are interconnected, that is, they are correlated.
Hence, to monitor these interconnected or correlated va-
riables is not simple but rather complex, especially for
manufacturing processes. The use of multiple univariate
control charts does not deliver a useful solution in th is si-
tuation. The problems are that, the overall probability of
signaling a false “out-of-control” situation is not contro-
lled and more seriously the correlation among the vari-
ables are ignored.
In recent years, multivariate statistical process control
(MSPC) procedures have enjoyed wide application in in-
dustry. This has resulted from expanded capability to
monitor the key variables of a process with sensor and
measurement technolog y, and th e wid espread av ailab ility
of computers and statistical software programs that in-
corporate multivariate SPC cap ability.
Simultaneously, there have been many new technical
developments that have made multivariate SPC more
useful. For example, many authors have investigated me-
thods of monitoring multivariate co ntinuous data. [1] de-
veloped the multivariate T2 statistic for quality control
purposes. Multivariate generalizations of the CUSUM
procedure have been studied by [2] and [3] developed
and investigated multivariate exponenti a l ly wei ght ed mov-
ing averages to identify quality problems. The use of
multivariate exponentially weighted moving averages in
monitoring multivariate data have been enhanced by [4].
Monitoring principal components of multivariate data
has been studied by [5]. [6] discussed multivariate mini-
max control chart and he us ed the joint probab ility distri-
bution function of the minimum and maximum stan-
dardized sample means to derive the control limits to make
decision if the process is in or out of control. However,
the derivation of the joint probability d istribution needed
to obtain the minimax control limits as discussed by [6]
is complex. In this paper we propose a Non-Linear Poly-
nomial function (NLP) approach to multivariate minimax
control chart to monitor con tinuous data as an alternative
approach to the use of joint probability for both the ma-

J. A. ADEWARA ET AL. 1540
ximum ([]
p
Z
) and minimum ([1]
Z
) limits used by [6].
The minimax control limits derived by [6] is modified
and the multivariate normal distribution is integ rated nu-
merically using Simpson’s one third rule to obtain a non-
linear polynomial (NLP) function. This NLP function is
then substituted and solved numerically using Newton
Raphson method to obtain the control limits for the mi-
nimax control chart.
2. Non Linear Polynomial Function
Polynomia ls are popu lar in curve an d surface rep resenta-
tions and many critical problems arising in Computer
Aided Geometric Design such as surface integration, are
reduced to finding the zero set of a system of nonlinear
polynomial equations () 0fx (2.1)
where 12 (, )
n
f
ff f and each i
f
is a polynomial
of independent variables 12
(, ,)
l
X
XX X

. Several
root-finding algorithms for multivariate polynomial sys-
tems (2.1) have been used in practice. Newton type
methods, which are classified as local solution techni-
ques, have been applied to many problems since they are
quadratically convergent and produce accurate results.
They, however, require good initial approximations of
the roots of the systems, and fail to provide full assur-
ance that all roots have been found. These limitation s can
be overcome by global solution technique, which can be
categorized into three different types as proposed by [7].
The different types are algebraic and hybrid methods,
homotropy methods, and subdivision methods. Among
these techniques, the subdivision methods have been
widely used in practice because of their performance and
efficiency. The Interval Projected Polyhedral (IPP) algo-
rithm proposed by [7] and [8] is one example, and it has
been successfully applied to various problems. One par-
ticular interest is locating zeros of a univariate applica-
tion of polynomial [9].
It is a critical problem in diverse fields such as control
theory and many literature has been devoted to it (see e.g.
[10]).
Most of the root finding algorithms, however, experi-
ence difficulties in dealing with roots with high multi-
plicity such as performance deterioration and lack of ro-
bustness in numerical computation. For example, the IPP
algorithm, which belongs to the su bdivision class of me-
thods, slows down drastically and suf fers from prolifera-
tion of boxes that are assumed to enclose roots. More-
over, since a root with high multiplicity is unstable with
respect to small perturbation , round-off errors dur ing flo-
ating point arithmetic may change the topological aspect
in such a way that a cluster of roots could be formed
around the root.
Solving univariate polynomials with multiple roots is
an important but difficult task. [9] collated nine methods
to bound multiple roots of polynomials and compared
them rigorously. He also proposed a new hybrid algo-
rithm which gives numerically nearly optimal bounds for
multiple roots of univariate polynomials. Even though
these methods work well in most cases, it is not easy for
a user to control the size of th e bound of a root in general.
[11] used the Sturm sequences to compute all roots of a
univariate polynomial, but his approach relies on the di-
vision of polynomials to compute Sturm sequences. So,
it is not numerically robust unless exact arithmetic or
symbolic computation is used.
This paper focuses on the particular case where the


,1,,iQ   in
111
1
(, ,)
(, ,)
N
QQ N
xfs s
x
fs s








are
f
functions
polynomials.
Thus, we write ()
X
fs

where f is a nonlinear
function :
N
Q
f
CC and 1Q
f
f constitute the com-
ponent of f. The source separation problem consists of
recovering the sources 1
N
s
s from the observation
1,, Q
x
x for all i,()
i
f
Cs

, where stands for
the set of polynomials in variables 1
()Cs
,,
N
s
s and with
coefficients in C. This restriction is partly justified by th e
difficulty to tackle the nonlinear case because of its gene-
rality. In addition, polynomials constitute an important
class of nonlinear models which may represent accept-
able approximations of certain nonlinearities.
Finally, an important reason to deal with this model is
the following:
Consider the case where the multidimensional source
vector belongs to a finite set .
Although seemingly restrictive, this situation is highly
interesting since it occurs in digital communications,
where the emitted source sequences belong to a finite
alphabet depending on the modulation used.

(1)( )
,,na saa   
An important observation is that if s

 and A is
finite, all instantaneous mixtures of the sources can be
expressed as polyno mial mixtures. This follo ws immedi-
ately from the fact that any function on a finite set can be
interpolated by a polynomial in a way similar to La-
grange polynomial interpolation [12]. It follows that po-
lynomial mixtures constitute the general model of non-
linear mixtures in the case of sources belonging to a fi-
nite alphabet.
The Model
111
1
(, ,)
(, ,)
N
QQ N
xfs s
x
fs s








(2.2)
is a polynomial, and in order to be able to resort to alge-

J. A. ADEWARA ET AL.1541
braic techniques, we will restrict the separator to the
class of polynomial func tions in 1,,
Q
x
x, that is,


,
ii
g
Cx  .
2.1. Simpson One Third Rule
The Simpson’s 1/3rd rule is a numerical method for find-
ing the integral ()d
b
a
f
xx
 within some finite limits a
and b. Simpson’s 1/3rd rule approximates ()
f
x with a
polynomial of degree two p(x), i.e., a parabola between
the two limits a and b, and then finds the integral of that
bounded parabola, and is used to represent the approxi-
mate integral ()d
b
a
f
xx
. The integral of the approxi-
mated function is the area under the parabola bounded by
the points a and b by the positive side of the x axis. The
quadratic function has three points common to the
function f(x), as follows: The end points of the approxi-
mate quadratic function p(x) is the same as the function
f(x) at points a and b. p(x) takes the same value of the
function f(x) at point ()mab 2.
Thus three points are fixed each in equal interval
amb and a parabola is drawn through these three
points f(a), f(m), f(b). The area under the parabola
through these points bounded by a and b with the
positive side of the X axis is found and used as the
approximated integral value. The iterative formula below
can be used to find the integral of a function f(x) using
Simpson’s 1/3rd rule.
012
() () ()
3
SI h
I
fxfx fx


(2.3)
2.2. Newton Raphson Method
The Newton-Raphson method is based on the principle
that if the initial guess of the root of f (x) = 0 is at i
x
,
then if one draws the tangent to the curve at ()
i
f
x, the
point 1
()
i
x
 where the tangent crosses the x-axis is an
improved estimate of the root.
Using the definition of the slope of a function, at x = xi
1
() 0
() tan
ii
fx
fx
x
x




 (2.4)
Form Equation (2.4) we have
1()
()
i
ii i
f
x
xx
f
x

 (2.5)
Equation (2.4) is called the Newton-Raphson formula
for solving nonlinear equations of the form ()0fx

.
So starting with an initial guess, i
x
, one can find the
next guess, 1i
x
 by using Equation (2.5). One can re-
peat this process until one finds the root within a desir-
able tolerance.
Algorithm
The steps of the Newton-Raphson method to find the
root of an equation () 0fx


are:
Step 1. Evaluate ()
f
x

symbolically
Step 2. Use an initial guess of the root, i xi, to estimate
the new value of the root , 1i
x
 as 1()
()
i
ii i
f
x
xx
f
x


Step 3. Find the absolute relative approximate error


as
1
1
*100
ii
i
xx
x






Step 4. Compare the absolute relative approximate
error with the pre-specified relative error tolerance,
,If
s



 , then go to Step 2, or else stop the algo-
rithm. Also, check if the number of iterations has ex-
ceeded the maximum number of iterations allowed. If so,
one needs to terminate the algorithm and notify the user.
2.3. Minimax Control Chart
The Minimax control chart developed by Sepulveda [6],
and as discussed in [13] is similar to the charts propo sed
by [14] and [15]. The minimax control chart uses the mi-
nimum and maximum standardized sample means to
make the decision if the process should be considered in
control or out of control. However, the minimax chart
uses both lower and upper control limits on both the ma-
ximum and minimum standardized sample means. This
is facilitated by the development of the capability to de-
termine the value of the joint density function of the
maximum and min imu m stan dar dized s ampl e means . Th is
not only facilitates a method for setting the control limits,
but also allows for the comparison of the performance of
the minimax chart relative to other charts through com-
putation of the out-of-cont rol a verage run lengt h.
Minimax control chart is used to standardize all p
means and to monitor the maximum and the minimum of
those standardized sample means. To do this, the sample
average vector 2
1
(, ,, )
p
X
XX X
 is calculated and
its elements are standardized using the expression:
()nX
Z






 (2.6)
where

is the population mean and

is the stan-
dard deviation. The vector [], 1,2,,
i
Z
Zi p is
now defined as the standardized sample mean vector.
The maximum sample mean ()
p
Z
is defined as the
maximum of the elements of the vector Z, that is,
[] ( ),max
p
i
Z
Z

Also, the minimum standardized sam-

J. A. ADEWARA ET AL. 1542
)ple mean [1] (
Z
is defined as the minimum of the ele-
ments of the vector Z.
The control limits that was proposed by [6] is modi-
fied to solve for both upper and lower of maximum and
upper and lower of minimum as given in the expression
below:




:()d1
:()d
:1( )d
:1( )d
UCf Zz
f Zz
UCf Zz
f Zz







[]
[]
[1]
[1]
p
p
L
LCL
L
LCL






 
 





 










(2.7)
3. Results
The data for this research work were collected from the
production line of a manufacturing company that pro-
duces soft drinks. The samples were drawn from the lines
on each variable of the production. The data are secon-
dary and multivariate in nature. The data had five vari-
ables which are: X1 = Contents in ml, X2 = Brev brix, X3
= pressure, X4 = Gas volume (CO2) and X5 = Tempera-
ture. Thus,

12
,,,
n
X
XX X
. We assumed that the
variables are normally distributed since we are dealing
with continuous data. The multivariate normal distribu-
tion was integrated numerically using Simpson’s one
third rule. Simpson’s rule is a numerical method that
approximates the value of a definite integral by using
quadratic polynomials. This approach was applied to the
multivariate normal distribution to obtain a non-linear
polynomial (NLP) function. This (NLP) function over-
comes the problem of obtaining the joint prob ability dis-
tribution needed for the control limits of both the maxi-
mum ([]
p
Z
) and the minimum ([1]
Z
) statistic. This
method was used to determine the position of the five
control limits of the chart stated in Equation (2.4). In
other to obtain the control limits, an algorithm was deve-
loped and implemented on C language to fit the polyno-
mial function in the form Z = 0.0024x5 + 0.000005x4 –
0.0444x3 – 0.00006x2 + 0.3805x + 0.4988. Using the ob-
tained polynomial equation, the algorithm in 2.2.1 was
then used to obtain the control limits for both the mini-
mum a nd t he ma x imum st a t ist i c s. T he nu m e ric a l s olut ion
for the control limits using the developed algorithm is
presented in the Appendix. The control limits for mini-
mum and maximum statistics for the five variables under
consideration are presented in Table 1.
Using the obtained control limits in Table 1, the pro-
cess under study was tested for stability. To test for the
stability of the process, Equ ation (2.6) was used to trans-
form the data to obtain the minimum and the maximum
Table 1. The upper and lower control limit for both maxi-
mum and minimum statistics.
UCL[p] LCL[p] UCL[1] LCL[1]
2.4185 1.954 3.0306 2.7195
–2.9458 –1.3148 –3.877 –3.3942
Table 2. The maximum and the minimum values.
Maximum0.0187250.0214260.021987 0.022794 0.0357959
Minimum0.005012–0.00226–0.03127 –0.03031–0.00954
values for the five variables. The obtained minimum and
maximum values are presented in Table 2.
4. Discussion of Result
The values in Table 2 are arranged from the lowest to
the highest. Thus minimum of Z[p] is 0.018725, and ma-
ximu m of Z[p] is 0.0357959. Also the minimum of Z[1] is
–0.031271 and the maximum of Z[1] is 0.005012. Com-
paring these values with the control limit in Table 1, the
result shows that the minimum and maximum values
obtained are within the con trol limits. Hen ce, the produ c-
tion process under consideration can be adjudged as be-
ing stable.
5. Conclusions
Minimax multivariate control chart is another sensitive
multivariate control chart that has upper and lower con-
trol limits for both maximum and the minimum statistics
for monitoring a multivariate process. The paper has ad-
dressed the use of numerical solution for obtaining the
control limits of the minimax con trol chart as an alterna-
tive to the use of joint probabilit y distribution.
6. References
[1] H. Hotelling, “Multivariable Quality Control—Illustrated
by the Air Testing of Sample Bombsights,” In: C. Ei-
senhart, M. W. Hastay and W. A. Wallis, Eds., Tech-
niques of Statistical Analysis, McGraw Hill, New York,
1947, pp. 111-184.
[2] W. H. Woodall and M. M. Ncube, “Multivariate Cusum
Quality Control Procedures,” Technometrics, Vol. 27, No.
3, 1985, pp. 285-292. doi:10.2307/1269710
[3] C. A. Lowry, W. H. Woodall, C. W. Champ and S. E.
Rigdon, “A Multivariate Exponentially Weighted Moving
Average Control Chart,” Technometrics, Vol. 34, No. 1,
1992, p. 46. doi:10.2307/1269551
[4] G. C. Runger, J. B. Keats, D. C. Montgomery and R. D.
Scranton, “Improving the Performance of the Multivari-
ate Exponentially Weighted Moving Average Control
Chart,” Quality and Reliability International, Vol. 15, No.

J. A. ADEWARA ET AL. 1543
3, 1996, pp. 161-166.
doi:10.1002/(SICI)1099-1638(199905/06)15:3<161::AID
-QRE215>3.0.CO;2-V
[5] C. M. Mastrangelo, G. C. Runger and D. C. Montgomery,
“Statistical Process Monitoring with Principal Compo-
nents,” Quality and Reliability International, Vol. 12, No.
3, 1996, pp. 203-210.
doi:10.1002/(SICI)1099-1638(199605)12:3<203::AID-Q
RE12>3.0.CO;2-B
[6] A. Sepulveda, “The Minimax Control Chart for Multi-
variate Quality Control,” Dissertation, Department of In-
dustrial and Systems Engineering, Virginia Polytechnic
Institute and State University, Blacksburg, 1996.
[7] N. M. Patrikalakis and T. Maekawa, “Shape Interrogation
for Computer Aided Design and Manufacturing,” Springer-
Verlag, Heidelberg, 2002.
[8] E.C. Sherbrooke and N. M. Patrikalakis, “Computation of
the Solutions of Nonlinear Polynomial Systems,” Com-
puter Aided Geometric Design, Vol. 10, No. 5, 1993, pp.
379-405. doi:10.1016/0167-8396(93)90019-Y
[9] S. M. Rump, “Ten Methods To Bound Multiple Roots of
Polynomials,” Journal of Computational and Applied
Mathematics, Vol. 156, No. 2, 2003, pp. 403-432.
doi:10.1016/S0377-0427(03)00381-9
[10] .M. McNamee, “A Bibliography On Roots of Polynomi-
als,” Journal of Computational and Applied Mathematics,
Vol. 47, No. 3, 1993, pp. 391-394.
doi:10.1016/0377-0427(93)90064-I
[11] H. S. Wilf, “A Global Bisection Algorithm for Comput-
ing the Zeros of Polynomials in the Complex Plane,”
Journal of the Association for Computing Machinery, Vol.
25, No. 3, 1978, pp. 415-420.
doi:10.1145/322077.322084
[12] A. C. David, L. John and O. Donal, “Ideals, Varieties and
Algorithms: An Introduction to Computational Algebraic
Geometry and Commutative Algebra,” Springer, 2nd Edi-
tion, Berlin, 1996.
[13] J. Rehmert, “A Performance Analysis of the Minimax
Multivariate Quality Control Chart,” M.Sc. Dissertation,
Department of Industrial and Systems Engineering, Vir-
ginia Polytechnic Institute and State University, Blacks-
burg, 1997.
[14] A. J. Hayter and K. L. Tsui , “Identification and Quantifi-
cation in Multivariate Quality Control Problems,” Jour-
nal of Quality Technology, Vol. 26, 1994, pp. 197-208.
[15] N. H. Timm, “Multivariate Quality Control Using Finite
Intersection Tests,” Journal of Quality Technology, Vol.
28, 1996, pp. 233-243.

J. A. ADEWARA ET AL. 1544
Appendix: Numerical Solution
The numerical solutions for the control limits is
554 352
( )0.00245100.4446100.38050.4988ZxXX XXX


The multiple integral of Equation (2.6) is solved below:

554 352
665 4532
0.00245100.4446100.38050.4988 d
0.0024510 0.444610 0.38050.4988
65432
X
XXX X
XXXXX
x
X


  



665 4532
766554 3 2
0.0024510 0.444610 0.38050.4988d
65432
0.00245 100.4446 100.38050.4988
766554433212
XXXXX
x
XXXXXX











766554 3 2
867 65543
0.00245 100.4446 100.38050.4988d
766554433212
0.00245 100.4446 100.38050.4988
768657 546435324123
XXXXXX
x
XXXXXX










 

9687 565
0.00245 100.4446100.38050.4988
76896578546743563245 1234
4
X
XXX X



 
X
The possible solutions for the numerical algorithm of the multiple integral are given below.
7.5249
4.2219 4.8351
4.2219 4.8351
2.9458
i
i




35 24
3.1708 102.0783 100.95),:
x
x roots 
0.24197 2.5085
0.24197 2.5085
2.4185
8.0042 4.6104
8.0042 4.6104
i
i
i
i




799857 76
7.9366 102.9763105.2857 101.6667 10
x
xx

x


7.5651
4.223 4.7726
4.223 4.7726
1.3148
i
i




35 24
3.1708 102.0783 100.05),:
x
x roots


2
2
5.8498101.2355 ,
5.8498 101.2355
1.1954
8.0048 4.6049
8.0048 4.6049
i roots
i
i
i







J. A. ADEWARA ET AL. 1545
79 98 57
12*3.141592654(7.9366 102.9763 105.285710
x
xx




7.4507
4.2262 4.928
4.2262 4.928
3.877
i
i




76 3524
1.6667103.1708 102.0783 10))0.95),:
x
xx

 root
0.36978 3.1125
0.36978 3.1125
3.0306
8.0031 4.6195
8.0031 4.6195
i
i
i
i




79 98 57
1(2*3.141592654(7.9366 102.9763 105.2857 10
x
xx




74962
4.2229 4.8741
4.2229 4.8741
3.3942
i
i




76 3524
1.6667 103.1708 102.078310))0.05,:
x
xx roots

 
0.30276 2.8096
0.30276 2.8096
2.7195
8.0037 4.614
8.0037 4.614
i
i
i
i



