The Performance of the MEWMA Control Chart When Parameters Are Estimated

Keywords: EWMA, Estimation Effect, Average Run Length, Statistical Process Control Average Run Length, EWMA, Estimation Effect, Multivariate Control charts, Statistical Process Control

In the present article, we study the effect of estimating the vector of means and variance-covariance matrix in the performance of the MEWMA control chart. We show through simulation that the performance of the MEWMA control chart is affected when the parameters are estimated compared to the known parameters case. The control limit that achieves the desired in control ARL when estimating the parameters is computed in several cases and formulas used to find approximately its values are given. Finally, an optimal design strategy for the MEWMA chart with estimated parameters is presented.

Introduction

Control charts are a valuable tool for detecting an out of control process. During the last decades there is an increasing interest for multivariate control chart techniques due to the fact that many processes involve several quality characteristics (variables) that are correlated. One of the most known control charts that belongs to this category is the multivariate exponentially weighted moving average (MEWMA) chart. The MEWMA control chart was introduced by [START_REF] Lowry | A multivariate exponentially weighted moving average control chart[END_REF] and it is used to detect fast small to moderate shifts in a process. [START_REF] Prabhu | Designing a multivariate EWMA control chart[END_REF] showed how to design this chart.

Generally speaking, the monitoring of a multivariate process includes two phases: Phase I and Phase II. In Phase I, one checks the stability of the process and estimates the unknown parameters. However, the main purpose in Phase II is to detect out-ofcontrol conditions as quickly as possible. The performance of a Phase I control chart is measured in terms of the probability of chart signals. On the other hand, the performance of a Phase II control chart is usually measured in terms of the parameters of the run length, where the run length is defined as the number of samples taken until the chart produces a signal. In order to use the MEWMA control chart during a Phase II application, the existing procedures for designing it are based on the assumption of known process mean vector µ µ µ µ and variance covariance matrix Σ Σ Σ Σ. In practice, however, µ µ µ µ and Σ Σ Σ Σ are usually unknown and we have to estimate these parameters using a historical data set (Phase I). Using parameters estimates with design procedures intended for known parameters can result in significantly deteriorated chart performance since these estimators add extra variability in the process.

Several authors have dealt with this problem for different univariate control charts, see e.g., [START_REF] Chen | The mean and standard deviation of the run length distribution of X charts when control limits are estimated[END_REF] and [START_REF] Jones | The performance of exponentially weighted moving average charts with estimated parameters[END_REF][START_REF] Jones | The run length distribution of the CUSUM with estimated parameters[END_REF]. [START_REF] Quesenberry | The effect of sample size on estimated limits for X and X control charts[END_REF] studied the effect of parameter estimation on the performance of the univariate Shewhart chart. This researcher suggested that at least 100 Phase I subgroups with sample size of 5 observations is needed to achieve expected statistical performance of the Shewhart chart. [START_REF] Jones | The performance of exponentially weighted moving average charts with estimated parameters[END_REF] showed that a number of subgroups more than 100 is needed to achieve expected statistical performance of the univariate EWMA T chart [START_REF] Nedumaran | On constructing 2 T control charts for online process monitoring[END_REF] and [START_REF] Champ | Properties of the 2 T Control Chart when Parameters are Estimated[END_REF]). None of the previous studies investigated the effect of parameter estimation on the performance of the MEWMA chart. This has been pointed by [START_REF] Jensen | Effects of parameter estimation on control chart properties: A literature review[END_REF] along with the conjecture that no attention has been given to the multivariate EWMA chart a fact that motivated this work. Our main focus in this paper is to investigate the effect of estimating µ µ µ µ and Σ Σ Σ Σ from an in-control Phase I data set on the performance of the MEWMA control chart used to monitor the process mean vector µ µ µ µ.

In Section 2, we present the MEWMA chart and its main properties. The effect of estimating the parameters on the performance of the proposed scheme is given in Section 3. In Section 4, we present the control limits that give the desired in-control average run length (ARL) when m samples of size n are used for estimating the parameters of the p-variate MEWMA chart. Section 5 presents the way to design the MEWMA control chart. Finally some comments and concluding remarks are given in Section 6.

The MEWMA control chart

Let i X , i= 1,2,… be ( 1 × p ) random vectors taken at regular time intervals and assumed to follow a p-variate normal distribution ) , (

0 0 Σ µ p N
, where 0 µ is the known p-dimensional vector of means and 0 Σ is the known variance-covariance matrix. Each vector records the values of p correlated quality characteristics that we want to monitor simultaneously and we further assume that these s ' i X , i=1, 2,… are independent. The random vector i X can be either the sample mean vector or individual observation vector at time i. The main concern is to detect shifts in the process mean vector µ from the target vector 0 µ . We assume without loss of generality that 0 µ = 0. The multivariate EWMA control chart statistic given by Lowry Z is usually set equal to the in-control mean vector of the process (i.e., 0 0 µ Z = ). A straightforward note is that if R=I then the MEWMA control chart is equivalent to the T² Chart.

et al. (1992) is 1 ) ( - - + = i i i Z R I RX Z , i= 1, 2, 3,… (1)
The MEWMA chart gives an out-of-control signal if

i i i Z Σ Z Z ' >h where i Z Σ is
the variance-covariance matrix of i Z and h is the control limit used to achieve the specified in-control ARL. The variance-covariance matrix of

i Z when r r r r p = = = = ... 2 1 is equal to ( ) Σ Σ Z ) 2 ( ] 1 1 [ 2 r r r i i - - - = .
The MEWMA chart statistic is usually constructed in terms of the asymptotic covariance matrix

Σ Σ Z ) 2 ( r r i - = .
(2)

When µ µ µ µ 0 and Σ Σ Σ Σ 0 are unknown, then m in-control Phase I samples, each of size n are used to estimate them. The in-control process mean vector µ µ µ µ 0 is estimated by

m m j j ∑ = = 1 X X , (3)
where j X is the j th sample mean vector. Also, the unknown variance-covariance

matrix Σ Σ Σ Σ is estimated by m m j j ∑ = = 1 S S , (4)
where j S is the within sample variance-covariance matrix.

For simplicity, we consider only charts with equal smoothing parameters and chart statistic calculated based on the asymptotic covariance matrix in Eq. ( 2). Lowry (1992) showed that the MEWMA chart with equal smoothing parameters is directionally invariant. [START_REF] Champ | Properties of the 2 T Control Chart when Parameters are Estimated[END_REF] showed that the in-control run length distribution of the directionally invariant T 2 chart statistic does not depend on the unknown process parameters. The in-control run length distribution of the T 2 chart statistic depends only on the number of variables p, the number of Phase I groups m, and the sample size within each group n. The out-of-control run length distribution depends on these constants as well as the statistical distance between the in-control and out-of-control mean vectors. This distance is usually referred to as the noncentrality parameter. Thus, the performance analysis of the T 2 chart can be given without knowledge of the in-control values of the process parameters or their estimates. [START_REF] Lowry | A multivariate exponentially weighted moving average control chart[END_REF] showed that the ARL performance of the MEWMA control chart with equal smoothing parameters (i.e., r r r r p = = = = ..

) depends only on the non-centrality parameter. This means that the MEWMA chart has the property of directional invariance. The use of unequal weights cancels out this property since in this case the ARL depends on the direction of the shift. It can be shown that the run length distribution of the directionally invariant MEWMA control chart with estimated parameters depends only on the values of p, m, n, and the smoothing parameter r. The proof is given in [START_REF] Champ | Properties of Multivariate Control Chart with Estimated Parameters[END_REF]. Rigdon (1995a,b) gave an integral and a double integral equation for the calculation of in-control and out-of-control ARLs of the MEWMA chart. [START_REF] Runger | A Markov chain model for the multivariate EWMA control chart[END_REF] used a Markov chain approximation to determine the run length performance of the MEWMA chart. Evaluating the run length performance of the MEWMA chart using the double integral equations or Markov chain approximation when the parameters are estimated from m in-control Phase I samples is intractable.

Thus, the results we present in this paper are rather based on extensive Monte Carlo simulation. [START_REF] Prabhu | Designing a multivariate EWMA control chart[END_REF] provided recommendations for the selection of parameters for a MEWMA chart. To detect quickly small sustained shifts in the process mean vector, it is usually recommended that one design the MEWMA chart using small values of the smoothing parameter r; see [START_REF] Lowry | A multivariate exponentially weighted moving average control chart[END_REF]. To achieve robustness to violations of the multivariate normality assumption, [START_REF] Stoumbos | Robustness to non-normality of the multivariate EWMA control chart[END_REF] recommended the use of smoothing parameters in the range from 0.02-0.05 for the MEWMA chart. However, as mentioned in [START_REF] Mahmoud | A multivariate adaptive exponentially weighted moving average control chart[END_REF], a very large sudden shift in the process mean vector may not be detected immediately if small values of r are used. When using a small value of r, the current observation receives a small weight r, and an observation far from the target value may not result in an immediate out-of-control signal. Moreover, the MEWMA charts with small values of r can build up an arbitrary large amount of inertia; see for example [START_REF] Woodall | The inertial properties of quality control charts[END_REF] and [START_REF] Mahmoud | A multivariate adaptive exponentially weighted moving average control chart[END_REF]. In this paper, we investigate the performance of the MEWMA chart with some different recommended values of the smoothing parameters r=0.05, 0.10, and 0.20 when the incontrol parameters are estimated.

Effect of estimated parameters on the performance of the MEWMA chart

When estimates are used instead of unknown parameters, the sampling distribution of the MEWMA chart statistics should account for the variability in the estimators. If the researcher does not take this variability into account, the in-control and out-of-control performance of the control charts can be strongly affected. The direct result of not taking this variability into account is the significant increase of the number of false alarms, especially when the sample size of the phase I data set is small. As the sample size of the phase I data set increases, the estimation error decreases. In this section we investigate using simulation the performance of the MEWMA chart with estimated parameters.

As mentioned in Section 2, when the mean vector and the covariance matrix are unknown they are usually estimated using an in-control reference sample consisting of m subgroups of size n from an in-control ) , (

0 0 Σ µ p N where 1 ≥ m and p n m > -) 1 (
. Then the estimators used are X µ = 0 ˆ and S Σ 0 = ˆ, where X and S are as defined in Eqs ( 3) and (4), respectively. [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF]

showed that if a random sample X 1 , X 2 , …, X n is taken from ) , ( 0 0 Σ µ p N
at sampling stage j, then the sample

F o r P e e r R e v i e w O n l y mean vector j X ~) , ( 0 0 n N p Σ µ
and the matrix

j n S ) 1 ( - ~) 1 , ( 0 - n Wishart p Σ
, where

) 1 , ( 0 - n Wishart p Σ
is the p-variate Wishart distribution with parameters 0 Σ and 1 n degrees of freedom. Accordingly, the in-control estimator of the process mean

vector X ~) , ( 0 0 nm N p Σ µ and the matrix [m(n-1) S ] ~)) 1 ( , ( 0 - n m Wishart p Σ .
As mentioned in Section 2, the in-control run length distribution of the MEWMA chart statistic depends only on p, m, n and r. Since the ARL performance of the MEWMA chart depends on the in-control parameters or their estimates only through the non-centrality parameter, we assume without loss of generality that 0 µ = 0 and I Σ = 0 . In our simulation studies we considered five different values for p (p=2, 3, 4, 5, or 6), nine values for m (m= 30, 40, 50, 70, 100, 150, 200, 300, or 500), four values for n (n= 3, 5, 10, or 15), and three values for r (r=0.05, 0.10, or 0.20). The values used for p, n, and m were chosen in order to represent small to large values of these parameters. The three selected values for the smoothing parameter r are the most recommended values in the literature.

In our simulation study, the procedure used to investigate the effect of p, m, n and r on the in-control ARL performance of the MEWMA chart can be described as follows:

1-Using 50,000 simulation runs, we first estimated the upper control limit h needed to produce an in-control ARL of 200 when different values of the smoothing parameter r is used in designing the MEWMA chart statistic with known parameters. The performance and design of a multivariate control charts based on an in-control ARL of 200 were studied by several researchers; see for example [START_REF] Lowry | A multivariate exponentially weighted moving average control chart[END_REF], [START_REF] Champ | Properties of the 2 T Control Chart when Parameters are Estimated[END_REF], and [START_REF] Mahmoud | A multivariate adaptive exponentially weighted moving average control chart[END_REF]. Table 1 gives the values of h according to different values of p and r. For example, when r=0.05 the upper control limit that gives the desired in-control ARL is h=7.36 in case of p=2, h=9.41 in case of p=3 and h=12.97 in case of p=5.

2-Then, we generated a random vector X from 3-Then we generated a random vector i X′ to represent the new Phase II process information observed at the i th sampling stage.

4-The MEWMA statistic defined in Eq. ( 1) was then calculated based on the estimated parameters X and S and compared to the upper control limit h calculated in Step 1.

5-Steps 3-4 are repeated until a signal is given. When the signal is given then we record the run length.

6-Steps 2-5 are repeated 50,000 times and the ARL is then calculated. The procedure IML in the statistical package SAS was used to generate the random vectors and matrices and to estimate the ARL values.

( [START_REF] Quesenberry | The effect of sample size on estimated limits for X and X control charts[END_REF] and [START_REF] Jones | The statistical design of the EWMA control charts with estimated parameters[END_REF]. Our simulation results show that the sufficient number of Phase I samples depends on the smoothing parameter r. These results show that larger number of Phase I samples are required when smaller smoothing constants are used. For example, while at least 2000 samples of size 5 is needed to achieve the expected statistical performance of the 2-variate MEWMA chart with estimated parameters when r=0.05, only 1400 and 1000 samples are sufficient to achieve it when r=0.10 and 0.20, respectively. As expected, monitoring a mean vector with a dimension more than 2 using the MEWMA chart with estimated parameters requires larger number of Phase I samples than those listed for p=2 to achieve the desired ARL performance. For example, when p=6 and r=0.20, more than 2500 samples of size 5 are needed so that the in-control ARL produced by the MEWMA chart is 200.

Corrected Limits for the MEWMA chart with estimated parameters

In some practical applications the data are plentiful and hence waiting until 1000 or more Phase I samples to accumulate is feasible. However, in most real-life on small number of samples, the statistical performance of the chart can be very poor.

On the other hand, waiting a long time to obtain the necessary data may result in the presence of undetected parameter shifts during this time. Therefore, several researchers have developed design procedures for some quality control charts that do not require the assumption of known parameters; see [START_REF] Quesenberry | The effect of sample size on estimated limits for X and X control charts[END_REF], [START_REF] Jones | The statistical design of the EWMA control charts with estimated parameters[END_REF], and [START_REF] Champ | Properties of the 2 T Control Chart when Parameters are Estimated[END_REF]. The main idea of these procedures was to use wider control limits for the chart in order to reflect the variability of the parameter estimates used in constructing the chart statistic. [START_REF] Champ | Properties of the 2 T Control Chart when Parameters are Estimated[END_REF] referred to these limits as the "corrected limits". In this section we develop using simulation the corrected limits of the MEWMA chart that give exactly an in-control ARL of 200 when m Phase I samples of size n are used in estimating the unknown parameters.

Tables 5-7 give the simulated corrected control limits for the p-variate MEWMA chart with parameters estimated from m samples of size n that give incontrol ARL of 200 when r=0.05, 0.10, and 0.20, respectively. These limits were obtained by the simulation algorithm (steps 2-6) described in Section 3. An additional binary search similar to that described in [START_REF] Champ | Properties of the 2 T Control Chart when Parameters are Estimated[END_REF] was used to obtain these limits. The main idea of the binary search was to adjust the limits until the desired ARL value was obtained.

(Insert Tables 5-7 about here)

The results in Tables 567show that the corrected limits depend on the size of m, n, p, and r. Smaller values m and/or n require wider corrected limits to achieve the desired in-control ARL. In practice, if m and/or n are different from the values given in Tables 567, the practitioner may use interpolation to determine an appropriate corrected control limit. Alternatively, we give estimated corrected limits for the pvariate MEWMA charts in Tables 8-10 when r=0.05, 0.10, and 0.20, respectively.

The linear equation models given in these tables are least squares estimates of the corrected limits using the logarithm of m (with base=10) as the independent variable in simple linear regression models. For example, when n=5, p=2, and r=0.05, Table 8 gives the estimate of the corrected limits as a linear function of the logarithm of m by: (5)

When n=5, p=2, and r=0.05, if m=210 samples are available, the estimated corrected limits calculated from Equation ( 5) is 8.14. Using 50,000 simulation runs, this estimated corrected limit produces an in-control ARL of 208. Each of the least squares estimates of the corrected limits listed in Tables 8-10 produced a coefficient of determination R 2 of at least 90%.

(Insert Tables 8-10 about here)

Optimal design of the MEWMA chart with estimated parameters

Although the corrected control limits given in Section 4 result in a chart with the desired in-control ARL, a subsequent increase in the out-of-control ARL of the MEWMA chart with estimated parameters over the chart with known parameters is expected. The optimal statistical design of an MEWMA chart requires not only determining the control limit h that satisfies the desired in-control ARL, but also the smoothing parameter r that minimizes the out-of-control ARL for a shift of a specified size measured in terms of the non-centrality parameter. As mentioned before, the out-of-control run length performance of the MEWMA chart does not depend on the parameters or their estimates. This performance rather depends on the values of m and n through the value of the corrected control limit, the smoothing parameter r, the number of monitored variables p, and the size of the shift measured in terms of the non-centrality parameter λ.

In this section, we compute the value of the smoothing parameter r and the corresponding corrected limit h that optimize the out-of-control performance of the p- 2. For each combination of r and h, a standardized shift in the non-centrality parameter of size n / δ is introduced and the out-of-control ARL is then estimated. This shift was made by changing the process mean vector from

µ =0 to 1 µ µ = such that 1 1 1 µ Σ µ - ′ = n / δ
. The optimal values of r and h are those that minimize the out-of-control ARL for the specified shift.

Tables 11 and12 give the optimal combination of r and h for different values of m and p when δ =0.5 and 1, respectively. The second column in both tables gives the optimal values of r in detecting shifts of various sizes for an in-control ARL of 200.

The following columns give the values of the upper control limit h that achieve the desired in-control ARL of 200 when m samples of size n=5 are used in estimating the unknown parameters and the out-of-control ARL corresponding to the specified shift.

The last two columns in these tables give the upper control limits and the out-ofcontrol ARLs when the MEWMA statistic is calculated based on known parameters.

(Insert Tables 11-12 about here)

Our simulation results show that the optimal value of r does not depend on the number m of Phase I samples or on whether the MEWMA chart statistic is calculated based on known or estimated parameters. Also, this optimal value does not vary much as p changes. The optimal values of r given in this paper are very close to the optimal values given in Table 2 in Lowry et al. (1992, p. 49).

The results in Tables 1112show that the out-of-control ARL performance improves significantly as the number of samples m increases. For example, when p=2, the out-of-control ARL for detecting a shift of size δ=0.5 decreases from 46.41 when m=30 to 29.72 when m=200. When known parameters are used in constructing the 2variate MEWMA chart statistic, the out-of-control ARL for detecting this size of shift is 26.50, as shown in Table 11.

Conclusions

In this paper we have studied the effect of estimating the parameters on the MEWMA control chart. Specifically, we have shown through simulation that the performance of the MEWMA chart is seriously affected if the vector of means and the covariance matrix are estimated based on a small number of Phase I samples. The value of the control limit is computed for several combinations of the number of samples and sample size so that the in-control performance of the MEWMA chart in the case of known and estimated parameters is equivalent. Finally, an optimal design of the MEWMA chart is presented, providing to the practitioners the suitable smoothing parameter for fast detection of an out of control situation for different shift sizes. .13 14.22 14.77 13.24 13.12 12.35 12.28 11.86 11.76 11.23 11.10 4 0.13 17.22 16.58 15.87 14.82 14.67 13.63 14.00 12.88 13.19 11.99 Page

  problem in the multivariate case has been studied only in the case of the 2

  Σ . The random vector X and the matrix S were the random information given in the parameter estimates stage (Phase I).

  number of Phase I samples is usually not available. As shown in Section 3, if the practitioner uses a control chart with parameters estimated based

  variate MEWMA chart for different values of m of size n=5. The MEWMA chart is usually recommended for detecting optimally small sustained shifts. In our study, two sizes of the process shift were considered; a standardized shift in the non-centrality parameter from λ=0 to n / δ λ = was considered, where δ =0.5 or 1. Three values were considered for p (p=2, 3, or 4) and 4 values for m (m=30, 50, 100, or 200). The procedure used to find the optimal design of the MEWMA chart with estimated parameters can be described as follows: given values of m and p, all the combinations of r and h that produce an in-control ARL of 200 were obtained based on Steps 2-6 in the procedure described in Section 3.

  All corrections suggested by the referee have been incorporated in the revised version. Page 7. The referee correctly pointed out that this proposition is wrong. We have deleted the proposition. Pages 8, 9, 10, 12, 13, 14. All corrections suggested by the referee have been incorporated in the revised version.For the referee's note in Page 9: The MEWMA chart produces an in-control ARL of around 197.2 when m=2500, p=6, r=0.2 and n=5.

  As for the effect of the smoothing parameter r used in designing the MEWMA chart, the results in Tables2-4show that as r increases the number of false alarms produced by the MEWMA chart decreases. For example, if m=30, n=3, p=2 and r increases from 0.05 to 0.10, the ARL increases from 85.82 to 94.68. If the same m, n, and p are used and r increases from 0.10 to 0.20, the in-control ARL increases from 94.68 to 111.18. Notice that the T 2 chart is a special case of the MEWMA chart when r=1. Our simulation study shows that if r increases from 0.20 to 1, the in-control ARL increases from 111.18 to 164.02.In order to reduce the variability in the sampling distribution of the estimates so that the ARL performance of a chart with estimated parameters is not significantly

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	30 to 500, the y
	ARL increases from 85.82 to 177.15. On the contrary, the number of false alarms
	produced by the MEWMA chart increases as p increases. Compared to the previous
	example, when m=30, n=3, r=0.05, and p=6, the in-control ARL produced by the
	MEWMA chart is reduced to 41.49.
	(Insert Tables 2-4 about here)

Insert Table 1 about here)

Tables 2-4 give the in-control ARL values of the MEWMA chart for different values of m, n, p, and r. As shown in these tables, the in-control performance of the MEWMA control chart with estimated parameters is strongly affected if the researcher uses the upper control limit h designed based on known parameters. For example, as shown in Table

2

, when the smoothing parameter r=0.05 and h= 7.36 are used in designing the MEWMA chart for monitoring bivariate quality characteristics (p=2) and the parameters are estimated using m=30 Phase I samples with n=3 observations, the in-control ARL produced by the chart is only 85.82. The in-control ARL of this chart when the parameters are known is 200. This example shows that the in-control run length distribution is seriously affected by ignoring the variability added by estimating the unknown parameters. As expected, as m and/or n increases, the number of false alarms produced by the MEWMA chart decreases, as shown in Tables

2-4

. For instance, in the previous example if m increases from affected, the number of Phase I samples m should increase to an appropriate level. For example,

[START_REF] Champ | Properties of the 2 T Control Chart when Parameters are Estimated[END_REF]

recommended a number of Phase I samples m>900/n when designing the 2-variate 2

T control chart with estimated parameters. Several researchers gave the appropriate number of samples required to achieve a desired incontrol ARL when a control chart with estimated parameters is used when the sample size n=5; see for example

Table 1 :

 1 The upper control limit value that produce an in-control ARL of 200 when known parameters are used in constructing the MEWMA chart statistic, given for different values of p and r.

	r

URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.

ca Communications in Statistics -Simulation and Computation

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	16

Table 2 :

 2 The ARL values when m Phase I samples with p variables, each of size n are used to estimate the unknown parameters and the UCL used is the one that gives incontrol ARL of 200 when the parameters are known. The smoothing parameter used in calculating the MEWMA chart statistic is r=0.05.

	m

Table 3 :

 3 The ARL values when m Phase I samples with p variables, each of size n are used to estimate the unknown parameters and the UCL used is the one that gives incontrol ARL of 200 when the parameters are known. The smoothing parameter used in calculating the MEWMA chart statistic is r=0.10.

	m

Table 4 :

 4 The ARL values when m Phase I samples with p variables, each of size n are used to estimate the unknown parameters and the UCL used is the one that gives incontrol ARL of 200 when the parameters are known. The smoothing parameter used in calculating the MEWMA chart statistic is r=0.20.

	m

Table 5 :

 5 The UCL values that produce in-control ARL of 200 when m Phase I samples with p variables, each of size n are used to estimate the unknown parameters.

	The smoothing parameter used in calculating the MEWMA chart statistic is r=0.05.
	p	n				m
		30	40	50	70	100	150	200	300	500
		3 10.25	9.82	9.45	8.98	8.62	8.28	8.08	7.85	7.66
	2	5 10.23	9.78	9.41	8.96	8.60	8.26	8.06	7.84	7.65
		10 10.20	9.74	9.37	8.96	8.60	8.26	8.06	7.83	7.65
		15 10.17	9.68	9.35	8.95	8.59	8.25	8.05	7.83	7.65
		3 14.17 13.25 12.70 11.95 11.30 10.75 10.45 10.13	9.86
	3	5 13.95 13.10 12.54 11.80 11.22 10.68 10.41 10.12	9.85
		10 13.73 12.95 12.40 11.68 11.21 10.65 10.39 10.11	9.85
		15 13.63 12.84 12.33 11.65 11.20 10.65 10.37 10.10	9.85
		3 18.10 16.63 15.80 14.70 13.85 13.10 12.65 12.23 11.85
	4	5 17.45 16.25 15.40 14.43 13.63 12.93 12.55 12.17 11.82
		10 17.00 15.97 15.26 14.35 13.57 12.90 12.53 12.15 11.80
		15 16.90 15.85 15.20 14.30 13.55 12.89 12.53 12.15 11.80
		3 21.97 20.00 18.73 17.40 16.23 15.27 14.70 14.10 13.70
	5	5 20.83 19.30 18.28 17.00 15.97 15.10 14.58 14.00 13.68
		10 20.31 18.88 18.10 16.80 15.85 14.98 14.53 13.98 13.65
		15 20.22 18.70 17.93 16.77 15.80 14.95 14.50 13.98 13.65
		3 25.95 23.40 21.65 20.00 18.61 17.33 16.72 16.05 15.48
	6	5 24.22 22.35 21.08 19.44 18.20 17.08 16.53 15.95 15.44
		10 23.45 21.75 20.65 19.15 18.00 16.97 16.44 15.93 15.43
		15 23.30 21.49 20.47 19.10 17.93 16.94 16.40 15.92 15.43

Table 6 :

 6 The UCL values that produce in-control ARL of 200 when m Phase I samples with p variables, each of size n are used to estimate the unknown parameters.

	The smoothing parameter used in calculating the MEWMA chart statistic is r=0.10.
	p	n				m
		30	40	50	70	100	150	200	300	500
		3 10.84 10.44 10.17	9.84	9.58	9.30	9.13	8.99	8.86
	2	5 10.82 10.40 10.14	9.81	9.55	9.27	9.11	8.97	8.85
		10 10.81 10.37 10.12	9.79	9.53	9.26	9.10	8.97	8.85
		15 10.80 10.36 10.11	9.77	9.52	9.26	9.10	8.97	8.85
		3 14.47 13.72 13.25 12.70 12.18 11.81 11.53 11.29 11.10
	3	5 14.12 13.45 13.08 12.54 12.09 11.76 11.47 11.27 11.09
		10 13.93 13.36 12.98 12.48 12.04 11.72 11.46 11.26 11.09
		15 13.87 13.32 12.93 12.42 12.02 11.69 11.46 11.25 11.09
		3 18.00 16.75 16.12 15.30 14.62 14.03 13.70 13.40 13.14
	4	5 17.32 16.38 15.80 15.07 14.40 13.92 13.65 13.35 13.11
		10 16.88 16.12 15.56 14.92 14.33 13.86 13.63 13.32 13.10
		15 16.79 16.01 15.47 14.82 14.28 13.83 13.62 13.31 13.09
		3 21.37 19.80 18.94 17.86 16.91 16.18 15.77 15.39 15.10
	5	5 20.24 19.16 18.35 17.40 16.63 15.99 15.65 15.35 15.06
		10 19.73 18.70 18.02 17.15 16.46 15.90 15.62 15.32 15.03
		15 19.58 18.55 17.85 17.08 16.40 15.87 15.60 15.30 15.00
		3 24.85 22.91 21.75 20.25 19.17 18.27 17.77 17.27 16.90
	6	5 23.25 21.78 20.82 19.71 18.72 17.96 17.59 17.18 16.82
		10 22.41 21.10 20.31 19.28 18.48 17.82 17.53 17.14 16.77
		15 22.23 20.96 20.19 19.22 18.43 17.76 17.50 17.11 16.74

Table 7 :

 7 The UCL values that produce in-control ARL of 200 when m Phase I samples with p variables, each of size n are used to estimate the unknown parameters.The smoothing parameter used in calculating the MEWMA chart statistic is r=0.20.

	m

Table 8 :

 8 The least square estimate of the UCL as a function of the logarithm of m. The estimate produces in-control ARL close to 200 when the smoothing parameter used in calculating the p-variate MEWMA chart statistic is r=0.05.

	n

Table 9 :

 9 The least square estimate of the UCL as a function of the logarithm of m. The estimate produces in-control ARL close to 200 when the smoothing parameter used in calculating the p-variate MEWMA chart statistic is r=0.10.

	n

Table 10 :

 10 The least square estimate of the UCL as a function of the logarithm of m.The estimate produces in-control ARL close to 200 when the smoothing parameter used in calculating the p-variate MEWMA chart statistic is r=0.20.

	n

URL: http://mc.manuscriptcentral.com/lssp E-mail: comstat@univmail.cis.mcmaster.

ca Communications in Statistics -Simulation and Computation

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	23

Table 11 :

 11 The optimal value of the smoothing parameter r and the corresponding UCL that give an in-control ARL of 200 when using m samples of size n=5 in estimating the parameters of the p-variate MEWMA chart used to detect a sustained shift in the non-centrality parameter from λ to .3746.41 9.55 38.05 8.86 32.31 8.37 29.72 7.72 26.50 3 0.05 13.95 54.43 12.54 44.46 11.22 36.93 10.41 33.57 9.37 29.78 4 0.05 17.45 59.31 15.40 48.60 13.63 41.36 12.55 36.63 11.20 32.23

	λ	+	5 . 0	/	n	.
	m

Table 12 :

 12 The optimal value of the smoothing parameter r and the corresponding UCL that give an in-control ARL of 200 when using m samples of size n=5 in estimating the parameters of the p-variate MEWMA chart used to detect a sustained shift in the non-centrality parameter from λ to

	λ	+	1	/	n	.
	m

Communications in Statistics -Simulation and Computation


Hongyan ; Li

Joern Meissner

Capacitated Dynamic Lot Sizing with Capacity Acquisition

Keywords: CAPACITATED LOT SIZING, CAPACITY PLANNING, HEURISTICS Approximation supply chain management, lot sizing, capacity, approximation, heuristics supply chain management, lot sizing, capacity, approximation, heuristics

One of the fundamental problems in operations management is determining the optimal investment in capacity. Capacity investment consumes resources and the decision, once made, is often irreversible. Moreover, the available capacity level affects the action space for production and inventory planning decisions directly. In this paper, we address the joint capacitated lot sizing and capacity acquisition problem. The firm can produce goods in each of the finite periods into which the production season is partitioned. Fixed as well as variable production costs are incurred for each production batch, along with inventory carrying costs. The production per period is limited by a capacity restriction. The underlying capacity must be purchased up front for the upcoming season and remains constant over the entire season. We assume that the capacity acquisition cost is smooth and convex. For this situation, we develop a model which combines the complexity of time-varying demand and cost functions and of scale economies arising from dynamic lot-sizing costs with the purchase cost of capacity. We propose a heuristic algorithm that runs in polynomial time to determine a good capacity level and corresponding lot sizing plan simultaneously. Numerical experiments show that our method is a good trade-off between solution quality and running time.

One of the fundamental problems in operations management is figuring out how to determine optimal investment in capacity. A firm's capacity determines its maximal potential production per time unit. Acquiring capacity is usually costly and time consuming, and once the investment is made, the cost is often partially or completely irreversible, as installed capacity is difficult to adjust in the short term. Moreover, the decision on how much capacity to acquire also strongly influences the action space for future operations planning. Obviously, acquisition of too much capacity wastes investment that could be used for other important operation activities such as new product development and marketing. Too little capacity means long waiting times, missed sales opportunities and lost revenue. Therefore, it is necessary to find an effective and comprehensive method to determine the proper capacity configuration for operations with specific planning horizons.

In this paper, we consider a single-production facility that produces a single product item to satisfy a known demand. Firms must determine optimal capacity and at the same time solve a capacitated lot-sizing problem. The major difference between our study and previous efforts to address capacitated lot-sizing problems, such as the well-known papers of [START_REF] Wagner | Dynamic version of the economic lot size model[END_REF] and [START_REF] Zangwill | Minimum concave cost flows in certain networks[END_REF], is that in our model, the capacity level is an internal decision. We consider capacity-acquisition, production, and inventory-holding costs and formulate the problem as a cost-minimizing Non-Linear Mixed Integer Programming (NLMIP) model. It belongs to a problem class with a quadratic constrains, which is generally N P-hard according to the classification given in [START_REF] Garey | Multi-resource investment strategies: operational hedging under demand uncertainty[END_REF]. In our case a pseudo-polynomial solution approach is available. We find that this obvious method is computationally unattractive, and we therefore develop a heuristic algorithm. Our numerical experiments show that our method results in substantial improvements in running time with only minor sacrifice in solution quality.

This study seeks to provide a building block for more complicated models involving lot sizing and capacity decisions, for example, multiple product or multiple stage capacity acquisition and lot sizing problem. Solutions of a pseudo-polynomial algorithm are used as benchmark to measure the performance of our heuristics algorithms. The remainder of this paper is organized as follows. We review the relevant literature in Section 2. Section 3 introduces the relevant notation and the basic model. In Section 4 we propose a heuristic to solve this problem. A computational study and numerical results are presented in Section 5. Finally, the conclusions and future research directions are given in Section 6. Two venues of research are relevant with our joint lot-sizing and capacity-acquisition problem including dynamic lot-sizing and capacity-investment studies. Here, we concentrate on the most representative problem or the one most closely relevant to our interest. The aim of capacity-acquisition decisions is to select the proper capacity that not only satisfies demand completely, but also reduces overcapacity. The research on capacity acquisition includes two major streams: the traditional mathematical programming models and the economic models.

The flexible capacity investment and management problems was addressed at a relatively early stage with mathematical programming methods. [START_REF] Fine | Optimal investment in production-flexible manufacturing capacity[END_REF] introduce a two-stage stochastic programming model and analyze the cost-flexibility trade-offs involved in the investment in product-flexible manufacturing capacity for a firm. They address the sensitivity of the firm's optimal capacity investment decision to the costs of capacity, demand distribution, and risk level. van Mieghem (1998) studies the optimal investment problem of flexible manufacturing capacity as a function of product prices, investment costs and demand uncertainty for a two-product production environment. He suggests finding the optimal capacity by solving a multi-dimensional news-vendor problem assuming continuous demand and capacity. [START_REF] Netessine | Flexible service capacity: optimal investment and the impact of demand correlation[END_REF] propose a one-period flexible-service capacity optimization and allocation model taking the capacity acquisition, usage, and shortage costs into account. While each paper mentioned above considered the multiple products and multiple resources problems with demand uncertainties, their focuses were limited to single-period models.

Apart from the studies on flexible capacity investment, many efforts have also been made to solve generalized capacity-investment problems. Harrison and van Mieghem (1999) develop a single-period planning model to incorporate both capacity investment and production decisions for a multiple-product manufacturing firm. This study yields a multi-dimensional descriptive model generated from the "news-vendor model", and gives qualitative insights into real-world capacity-planning and capital-budgeting practices. Nevertheless, the decisions on optimal capacity investment are highly generalized, and the production plan decisions were not explicitly presented. [START_REF] Van Mieghem | Newsvendor networks: inventory management and capacity investment with discretionary activities[END_REF] extend the work of Harrison and van Mieghem (1999) to include an operations environment with multiple products, production processes, storage facilities, and inventory management. Moreover, they investigate how the structural properties of a single period extend to a multi-period setting. They also improve previous studies by considering some inventory-management issues. Since the news-vendor model is developed and applied to capacity-decision problems, it has been an important analysis technique to model and solve complex capacityoptimization problems under uncertainty. [START_REF] Burnetas | Future capacity procurements under unknown demand and increasing costs[END_REF] propose a newsvendor-like characterization of the optimal production policy on capacity under unknown demand and increasing costs within a finite horizon with discrete time periods. The approach focuses on the trade-off between increasing production costs and the learning mechanism about demand, neglecting set-up costs and capacity-supplying limitations.

Although the research has produced masses of data, lot-sizing remains one of the most difficult problems in production planning. This subject has been studied extensively in the literature. More than 50 years ago, [START_REF] Wagner | Dynamic version of the economic lot size model[END_REF] developed a forward algorithm for a general dynamic version of the uncapacitated economic lot-sizing model. The zero-inventory ordering theorem is a key contribution in this paper for the uncapacitated cases. Although many alternative algorithms have been presented, the dynamic programming method remains the major approach to solving lot-sizing problems. More recent studies consider a dynamic lot-sizing model with general cost structure. [START_REF] Federgruen | A simple forward algorithm to solve general dynamic lot sizing models with n periods in o(n log n) or o(n) time[END_REF] present a simple forward algorithm which solves the general dynamic lot-sizing model in O(T log T ) time and in O(T ) under mild assumptions on the cost data. This is an key improvement over the previously recommended well-known shortest path algorithm solution in O(T 2 ) space. [START_REF] Wagelmans | Economic lot sizing: an o(n log n)algorithm that runs in linear time in the wagner-whitin case[END_REF] extend the range of allowable cost data to allow for coefficients that are unrestricted in sign. They develop an alternative algorithm to solve the resulting problem in O(T log T ) time.

The uncapacitated lot-sizing problem is however an ideal case and hardly applicable to real-world operations. Capacity constraints always heavily influence production-plan decision-making. Furthermore, the general capacitated lot-sizing problem is N P-hard, see [START_REF] Bitran | Computational complexity of the capacitated lot size problem[END_REF]. For the special case of a constant capacity restriction over the entire planning horizon, a number of efficient algorithms are capable of calculating an optimal production plan. For example, [START_REF] Florian | Deterministic production planning with concave costs and capacity constraints[END_REF] present an algorithm with the computational complexity O(T 4 ) for the capacitated lot-sizing problem and explored the important properties of an optimal production plan. The optimal plan consists of a sequence of optimal sub-plans. [START_REF] Baker | An algorithm for the dynamic lot-size problem with time-varying production capacity constraints[END_REF] discover some important properties of an optimal solution to the problem when the production and inventory-holding costs are constant. Some other studies tried to relax the strict cost-structure restrictions in the algorithms reviewed above. [START_REF] Kirca | An efficient algorithm for the capacitated single item dynamic lot size problem[END_REF] present a dynamic programming-based algorithm with the computational complexity of O(T 4 ) and [START_REF] Shaw | An algorithm for single-item capacitated economic lot sizing with piecewise linear production costs and general holding costs[END_REF] propose (2007), van Hoesel and[START_REF] Van Hoesel | An O(T 3 ) algorithm for the economic lot-sizing problem with constant capacities[END_REF] and [START_REF] Chen | A new dynamic programming algorithm for the single item capacitated dynamic lot size model[END_REF] etc. In addition, heuristic algorithms are also applied to solve dynamic lot-sizing problems for efficiency. For example, [START_REF] Kiran | A combined heuristic approach to dynamic lot sizing problems[END_REF] proposes a combined heuristic algorithm based on the performance analysis of Silver-Meal heuristics. [START_REF] Alfieri | Lp-based heuristics for the capacitated lot-sizing problem: the interaction of model formulation and solution algorithm[END_REF] consider the application of trivial LP-based rounding heuristics to the capacitated lot-sizing problem.

Given the large number of research results on lot sizing it is impossible to examine carefully all of them. We refer interested readers to the following review papers. [START_REF] Karimi | The capacitated lot sizing problem: a review of models and algorithms[END_REF] review single-level lot-sizing problems, their variants and solution approaches. The authors introduce factors affecting formulation and the complexity of production-planning problems, and introduce different variants of lot-sizing and scheduling problems. Both exact and heuristic approaches for the problem are discussed. [START_REF] Jans | Modeling industrial lot sizing problems: a review[END_REF] present an overview on recent developments on the deterministic dynamic lot-sizing, focusing on the modeling of various industrial extensions and not on the solution approaches. The authors first define several different basic lot-sizing problems, and propose some extensions of these problems.

However, these studies all address capacity-investment or production-planning problems separately. The implications of combining these problems are rarely discussed. As an exception, [START_REF] Atamturk | Capacity acquisition, subcontracting, and lot sizing[END_REF] study a problem on capacity acquisition, subcontracting, and lot sizing. This is the only study we have come across that is closely related to our studies. However, the authors only discuss some special cases of production and holding-cost structure. Moreover, the study still focus on solving a series of capacitated lot-sizing problems discretely, causing the computational complexity to increase exponentially with the number of planning periods and demands. Additionally, [START_REF] Ahmed | Dynamic capacity acquisition and assignment under uncertainty[END_REF] study a dynamic capacity-acquisition and assignment problem in a simplified operations setting to determine the resource capacity and allocation of the resources to tasks. This study actually proposes a capacity-expansion and planning approach without considering inventory carry-over costs and the determination of production plans.

In summary, while progress has been made on investigating capacity-acquisition and lot-sizing decisions, the research has as yet yielded few results that pertain to joint optimization of capacity acquisition and production decisions under a capacitated lotsizing cost structure, even for a single firm. In this section, we analyze the capacity-acquisition and lot-sizing problem. A firm has to determine the optimal capacity to purchase and set a corresponding lot-sizing plan simultaneously.

The firm produces an item or product that consumes a common resource during its production. The amount of the resource the firm purchased is assumed to be the capacity limit in a dynamic lot-sizing setting. An example of this might be the number of trucks to lease, the work force to hire and other supportive activities for production.

The firm has to purchase the capacity for the entire planning horizon and can then use the capacity over the planning horizon. The capacity must satisfy the demand constraints and the excess capacity will be disposed of without extra disposal costs.

The production plan will be considered in a planning horizon of T periods. If the firms face a natural sales season introducing a new model or variant in each season, a natural choice of T arises, e.g. T = 52 weeks in the automobile manufacturing industry operating with a weekly production and sales schedule. Otherwise T is chosen large enough to ensure that the firms' decisions pertaining to the initial periods of the planning horizon are not affected by this truncation of the planning process.

The firm has a demand stream during the planning horizon, known only to the firm itself and following some predictable seasonality pattern. Thus, let

d t
= the demand faced by firm in period t, t = 1, . . . , T

The firm produces goods via a production and distribution process that, in principle, allows for inventory replenishment in each period. As in standard dynamic lot-sizing problems, we assume that fixed as well as variable production costs are incurred as well as inventory carrying costs, which are proportional to each unit end-of-the-period inventory. We assume that all fixed order costs stay constant over the planning horizon, while all other cost parameters may fluctuate in arbitrary ways. We define the cost parameters and decision variables as follows.

Costs: The firm needs to acquire the capacity in question on a spot market prior to the season.

f = the
We name the acquisition cost as A(C) and assume it is smooth and convex. Such an assumption is reasonable, among other explanations, when the purchase of the firm influences the market price. Our proposed algorithm works with any function A(C)

that is a continuous and convex function. As a simple illustration, let the market price for the resource be p = Λ+θC where Λ and θ are positive constants. Hence the capacity acquisition cost is:

A(C) = p • C = C(Λ + θC) (1)
We make following assumptions that are fairly standard in the peer literature. The inventory at the beginning of the planning level and the end of the planning horizon is zero respectively. Demand shortage is not allowed, because, for the deterministic case, it is optimal to pursue the 100% service level. The setup times are not considered.

For deterministic setup times, it is easy to be included after optimal lot-sizing strategy is determined by moving the setup period forward. For the stochastic setup times, it would result in a completely different problem. Capacity can be acquired at any positive amount. While one could argue that these assumption are far from reality in many cases, we chose this setup for two reasons. Firstly, we want to investigate the principal relationship between capacity and cost in the lot-sizing context and to view our model as a building block for further extensions. Secondly, many extensions make the model actually easier to solve or are theoretically not interesting. We will discuss some possible extensions in the directions for future research in the conclusions of this article.

This gives rise to the following formulation of the problem P:

z = min T t=1 (a t x t + h t I t + f y t ) + A(C) (2)
I t = x t -d t + I t-1 , ∀ t = 1, . . . , T (3a)
x t ≤ Cy t , ∀ t = 1, . . . , T (3b)

I 0 = I T = 0 (3c) x t ≥ 0, I t ≥ 0, y t ∈ {0, 1}, C ≥ 0, ∀ t = 1, . . . , T . ( 3d
)
where the objective function (2) minimizes the production and inventory-holding costs as well as the acquisition costs of the capacity. Constraints on the problem are: Equation (3a) ensures that inventory is balanced; Production is restricted by (3b); Equation (3c) sets initial and final inventories to zero; and the bounds of the variables are restricted by (3d). Solving the model entails simultaneously determining the optimal capacity, order periods, and production amounts in each order period. Capacity is assumed to be a continuous variable, meaning that capacity can be acquired at any non-negative level. It would be possible to linearise the quadratic constraints (3b) using a big-M formulation by replacing each constraint by two separate ones. We have tried to solve the resulting programme, which has a quadratic objective function in the case of linear capacity acquisition price, using CPLEX and found this computationally unattractive. We report on this in §5.

4 The Heuristic

Basic idea of the heuristic

The simultaneous calculation of an optimal capacity and an optimal production plan as explained above is a mixed integer nonlinear programming problem. This problem class is generally N P-hard according to [START_REF] Poljak | Convex relaxations of (0, 1)-quadratic programming[END_REF] and [START_REF] Bussieck | Mixed integer nonlinear programming[END_REF]. However, the capacitated lot-sizing problem with constant capacity can be solved in polynomial time. For example, [START_REF] Florian | Deterministic production planning with concave costs and capacity constraints[END_REF] suggest

an O(T 4 ) algorithm, and alternative approaches are also suggested by van Hoesel and [START_REF] Van Hoesel | An O(T 3 ) algorithm for the economic lot-sizing problem with constant capacities[END_REF] and [START_REF] Chen | A new dynamic programming algorithm for the single item capacitated dynamic lot size model[END_REF] that run in O(T 3 ) time. Therefore, the problem P can be solved by discretizing the interval of potential values for the capacities and solving for each of those values. So it is not N P-hard in the strong sense, and can be solved in pseudo-polynomial time.

Solving problems with reasonable sizes by discretizing the solution space for the capacities with CPLEX, although theoretically satisfactory, has shown to result in large computational times that make such a methodology impractical (see §5 for details). Therefore, in this section, we develop a O(T 3 log T ) heuristic algorithm that improves the computational efficiency dramatically.

To facilitate the presentation of our algorithm, we use the following notation. We define D(t) = t j=1 d j to be the cumulative demand in the first t periods, t = 1, . . . , T ;

X(t) = t j=1
x j to be the cumulative production level in the first t periods, t = 1, . . . , T ;

H(i, j) = j-1
k=i h k to be the cost of holding a product from period i to period j, ∀ 1 ≤ i < j ≤ T ;

H(i) =

T k=i h k to be the cost of holding a product from period i to the end of the planning horizon;

C n min = to be the minimum capacity that allows a feasible solution with n setups; Θ(n) = {1, ℓ 2 , . . . , ℓ n } to be a setup strategy with the fixed setup number n, n = 1, . . . , T . The orders in periods 1, ℓ 2 , ..., ℓ n obey the assumption that the available capacity is at least C n min .

In analogy to the algorithm presented by [START_REF] Federgruen | Competition under time-varying demands and dynamic lot sizing costs[END_REF], who present an algorithm for a combined pricing and uncapacitated lot sizing problem, the heuristic developed here considers each possible number of setups n, n = 1, . . . , T separately and determines the best capacity and production plan. We solve the following problem:

π * (C) = K n (C) + A(C) (4) = min n min C (nf + F n (C) + C (Λ + θC)) (5)
where the function F n (C) represents the production and inventory cost for a fixed setup number n.

The algorithm consists of three major steps:

Step 1 For each setup numbers n = 1, • • • , T , construct an initial solution with the minimal capacity that allows a feasible solution;

Step 2 Under each setup number n, update lot size plan and calculate the cost savings with the incremental capacity. This allows us to determine the best capacity acquisition level and lot size plan for each n respectively.

Step 3 Determine the optimal capacity and lot size solution by comparing the total cost over the setup numbers. The algorithm details are presented in the following sub-sections. We also analyze the complexity of the algorithm in Section 4.5.

Construction of the initial solution

For each number of setups n, we first find the minimal capacity that allows a feasible solution to the problem. This minimal capacity C n min can be calculated as follows:

C n min = max D(T ) n , max t=1,...,T D(t) t ∀n = 1, . . . , T (6)
After determining this minimal capacity, we find the number of order periods necessary for each fixed setup number n = 1, 2, . . . , T . We start with a solution that places the orders as late as possible under the minimal feasible capacity C n min , and then we improve the solution by shifting the orders forward or backward if this is beneficial.

The procedure is fully described in Algorithm 1. While it does not yield the optimal solution in general, in the important case of no prevailing speculative cost motives and C n min being determined as the average demand per period, it does result in an optimal initial solution:

Proposition 1 Assume that there is no speculative cost motive, i.e. a(s) + H(s, t) ≥ a(t) for all 1 ≤ s < t ≤ T , and that C n min =

D(T )

n , then Algorithm 1 results in an optimal solution for the fixed setup number n.

Proof : Let the initial production strategy from the Algorithm 1 be Θ 0 = {ℓ 0 1 , ℓ 0 2 , . . . , ℓ 0 n }, and moreover, since

C n min = D(T )
n , the production quantity in each setup period has to be C n min in order to satisfy demands. The proposition will be proved if we show the minimal cost π * = π (C n min |Θ 0 ). Suppose that the strategy Θ 0 is not optimal given the condition described in Proposition 1, there exists another production strategy Θ = {ℓ 1 , ℓ 2 , . . . , ℓ n } which makes π (C n min |Θ) ≤ π (C n min |Θ 0 ). According to the algorithm, the setups ℓ 0 i , i = 1, . . . , n cannot be postponed in order to satisfy the feasibility of solution, thus, there must exist at least one i, so that ℓ 0 i-1 < ℓ i < ℓ 0 i . This means that a(ℓ i ) + H(ℓ i , ℓ 0 i ) ≤ a(ℓ 0 i ). It contradicts the assumption of no speculative cost motive, a(s) + H(s, t) ≥ a(t) for all 1 ≤ s < t ≤ T . Thus, Algorithm 1 results in an optimal solution.

Update with increased capacity

Having found an initial solution, we update it with increased capacity. We introduce the following additional notation: The list of potential saving opportunities Ω is created first, and then elements of potential savings Ω are converted to a list of active savings Φ that we pursue at a given capacity increase. Each time a saving opportunity is exhausted, we check whether another element can be brought from Ω to Φ. Once Ω is empty, stop the algorithm. Each element of Ω is a quadruplet of the form ℓ -, ℓ + , δ, ǫ , ℓ -represents the period in which production is to be decreased, ℓ + is the period in which production is to be increased, and δ is the potential cost saving per unit, and ǫ denotes the maximum number of units for which the savings opportunity can be exploited.

After finding the initial solution, we update the production and lot sizing plan while the capacity increases. For any given number of order periods, we examine the possibility of improving the solution by using the additional capacity that the company might acquire by comparing the cost of such a change between two adjacent order periods. The two options are: either a shift of production to a previous order period or a postponement to a later order period. The first case, shifting the production earlier, creates no problems and can be repeated until the decreasing order period reaches zero. A postponement is potentially problematic, but can be done either until the first decreasing period has reached zero production level or until a further decrease leads to an infeasible solution. The maximum decrease is given by:

ǫ = min    x ℓ i ,   i k=1 x ℓ k - ℓ i -1 k=1 d k      (7)
In Algorithm 2, under the fixed setup number n, we compare each pair of sequential setups in period ℓ i and ℓ i+1 , i = 1, . . . , n -1 to determine {ℓ -, ℓ + , δ, ǫ}, and adding it to Ω.

Based on the saving opportunities matrix generated from the Algorithm 2, we sort the potential savings candidates Ω. Next, the Algorithm 3 moves to realize the savings.

In order to keep the linear decrease of lot sizing cost, we consider the capacity increases in a variable step size that is the minimum value of ǫ in the active savings candidate list Φ. The value of the current capacity adding a step size will be a breakpoint of capacity increasing. Upon reaching one of the breakpoints, the savings opportunity has been if a l i + H(l i , l i+1 ) < a l i+1 then 5:

Insert new element in Ω: (l i+1 , l i , a l i+1 -a l i -H(l i , l i+1 ), x l i+1 )

6:

else 7:

if

a l i + H(l i , l i+1 ) > a l i+1 then 8: if X(l i ) -D(l i+1 -1) < x l i then 9:
Insert new element in Ω:

{l i , l i+1 , a(l i ) + H(l i , l i+1 ) -a l i+1 , x l i } 10: else 11: Insert new element in Ω: {l i , l i+1 , a l i + H(l i , l i+1 ) -a l i+1 , X(l i ) -D(l i+1 -1)} 12:
end if

13:
end if 14:

end if 15: end for exhausted and is removed from the calculation. We have two options: either we stop when one order period has reached zero, with the reasoning that we can reach a similar solution in a run with n -1 setups or, since there is no harm from the point of view of complexity, we can proceed until our list is empty.

According to the heuristic procedure described above, some structural properties of the lot-sizing function K n (C) are realized and we clarified them in Lemma 1 below.

An outline of proof is provided to help illustrate the algorithm and the Lemma.

Lemma 1 For a fixed setup number n, the lot-sizing cost function K n (C) is piecewiselinear, non-increasing and convex in capacity.

Proof : The lot-sizing cost function is K n (C) = F n (C) + nf . Since fixed setup cost is constant, if the total production and inventory cost function F n (C) is piecewise-linear decreasing in capacity. Given a production plan x ℓ 1 , . . . , x ℓ n , the production and inventory cost function is

F n (C) = n i=1   a ℓ i x ℓ i + ℓ i+1 j=ℓ i +1 h j X(j) -D(j)    . ( 8
)
In order to prove that F n (C) is piecewise-linear decreasing in capacity, the following three properties of the function need to be proved respectively (all discussion below is based on a fixed setup number n):

(1) F n (C) is non-increasing in capacity.

If capacity increases to be C ′ , the production plan {x ℓ 1 , . . . , x ℓ n } is still feasible, and the decision space is broader, therefore, we have at least (3) F n (C) is continuous and convex.

F n (C ′ ) ≤ F n (C).
In the heuristic algorithms, a new breakpoint of capacity increase is always calculated based on the capacity level of the previous iteration. In addition, the solution of the production plan of an iteration is always the initial solution of the next iteration. Therefore, we see that the cost function F n (C) is continuous. Moreover, since the slope of the function F n (C) results from picking various elements from Ω, the convexity is a direct result of our picking elements in decreasing order of their savings.

Finally, given the setup number n, the fixed setup cost is constant, and therefore, the conclusion holds.

On a more abstract level, the above result restates a well-known result from parametric analysis for linear programmes, see for example [START_REF] Dinkelbach | Sensitivitätsanalysen und parametrische Programmierung. Ökonometrie und Unternehmensforschung[END_REF]. However, the above provides some intuition for our algorithm. For an illustration of Lemma 1, see Figure 1 selected from a numerical example discussed in Section 5.

Calculation of the optimal capacity

Upon obtaining the piecewise-linear functions for each individual setup number n, we calculate the optimal capacity to acquire by finding the appropriate breakpoint. According to the Lemma 1, and finite possible setup numbers, the optimal solution is obtained by comparing the minimal costs of all possible setup numbers. The optimal While this procedure may not be optimal, our computational experiments show that, in many cases, our results are very close to optimality. In the numerical study Section, we compare the heuristics solution with the solution obtained by a full enumeration over the discretized decision space method and using CPLEX 11.0 to solve the individual instances.

At last, we would like to comment on some possible variants on the model and algorithm. The algorithm is robust if the acquisition cost A(C) is a step function. The only modification is to consider a series of sub-capacity intervals caused by the step function within the entire capacity range. However, the decreasing or concave capacity acquisition function will not lead to the results we achieved.

Complexity

To assess the complexity of the algorithm, first, the individual complexity of the major algorithm steps are described. Note that the steps taken to solve the problem have to be repeated T times, once for each potential setup n. Then in each iteration under

Numerical Example

In this section, we present computational examples for our heuristic. Using the heuristic algorithm we developed for the capacity acquisition and lot sizing problem, a numerical study is carried out to show the robust performance of the algorithm. It is assumed that the firm faces a planning horizon of T = 54 periods with varying seasonal demand. The demand behaves according to:

d t = β t * d (9)
We consider six different seasonality patterns {β t : t = 1, . . . , 54} as follows:

(I) Time-invariant demand functions: β t = 1 ; t = 1, . . . , 54

(II) Linear Growth: β t = 0.25 + 1.5

(VI) Cyclical Pattern:

β t =        0.25 + 0.75(t -1) , t = 1, . . . ,
where tmod6 denotes t modulo 6. The first pattern reflects a situation where demand The results indicate that the heuristic algorithm performs effectively and efficiently.

First, the gaps between the heuristic and the pseudo-optimal solutions are very small with an overall average gap of 2.71%; this is acceptable given the extremely short computational times of around one second. Second, under the different demand patterns, the average gap does not vary dramatically. For the constant demand pattern, the average gap is the least, about 1% or less. For the holiday demand scenarios, the average gaps are higher but remains below 5% (see column 11 in Table 2).

Conclusion

In this paper, we consider the capacitated lot-sizing problem with capacity acquisition. We develop an efficient heuristic that solves the capacity acquisition, production, inventory decisions simultaneously with a complexity of O(T 3 log T ). Our numerical study shows that our heuristic algorithm performs well while using substantially less time compared to pseudo-optimal approach where the potential capacity space is discretized, while losing only a modest amount of accuracy. Our algorithm can serve as a building block for more sophisticated problems involving capacity-acquisition and dynamic lot-sizing decisions. Given that the demand for a firm's product is often affected by its pricing decisions, our future research will consider price-dependent demands. Additionally, more realistic problem settings, for instance, multiple stage production or multiple product lot-sizing and capacityacquisition should also be investigated. However, we would like to note that many extensions that aim to reflect better that reality might make the mode easier to solve.

  algorithm for the capacitated lot-sizing problem with general holding costs and piecewise linear production costs. The algorithm of the latter reduces the computation time to O(T 2 d), where d is the average demand when production cost is linear. Many other contributions in this area includeFedergruen et al.

  i = 1, . . . , T , where ξ i = {0, 1}; Φ = a list of active savings generated from Ω; Γ = (ǫ min , Savings) to be the executive list to update the lot sizing plan in each iteration of computation.

  Build sorted list of potential savings opportunities Ω 1: Given: Set of order periods Θ(n) = {1, ℓ 2 , . . . , ℓ n } 2: new list Ω 3: for i = 1 : 1 : n -1 do 4:

  F n (C) is piecewise-linear in capacity. According to the initial solution from Algorithm 1 and Algorithm 2, search all cost saving opportunities and record them in array Ω which allows the capacity to vary in the range [C n min , D(T )]. Furthermore, by Algorithm 3, we deal with the saving opportunities array Ω. According to the heuristic procedure, the computation includes a finite number of iterations based on different capacity levels. In each iteration, we define and calculate an active cost-saving array Φ = {φ m , m = 1, 2, . . . , M}, where φ m = {ℓ -, ℓ + , δ m , ε m }. For the detailed steps please refer to the algorithm. From the array Φ, we determine a capacity increase quantity ∆C = min ǫ m with the unit cost saving M m=1 δ m . Thus, cost function F n (C) is linear non-increasing in capacity interval (C, C + ∆C]. Capacity level (C + ∆C) is a new breakpoint of capacity increase.

Figure 1 :

 1 Figure 1: An example of cost variation with the capacity increase under a fixed setup number

  n, finding the minimal feasible capacity can be done in O(T ); next, finding the initial solution takes O(T ) for the first phase and O(T ) for the final initial solution;Third, in the update procedure of initial solution with the capacity increase, finding potential savings again takes O(T ), since n pairs at most have to be evaluated. After obtaining the list of potential savings, this list has to be sorted once, which takes O(T log T ) using Quicksort or a similar algorithm. Finally, searching the list for potential savings to determine each breakpoint can be done in O(T ). Given that there are at most n breakpoints, this leaves us with a complexity of O(T ) to update the solution from the previous lower capacity level.Considering the relationships (paralleled or hierarchical) between the steps, the integrated algorithm complexity is as follows. For each setup, we find the minimum capacity, find the initial solution, optimize the initial solution, and finally compare the optimal solution for each setup number which cause a complexity of O(T3 ). Based on each initial solution under C n min , the improvement procedure including the determination of Ω with a complexity of O(T ), the sorting of potential savings adding another O(T log T ) and updating of solution adding O(T ) again. Doing this totally introduces a complexity O(T log T ). The final comparison of each solution of each setup number gives a complexity of O(T 2 ). Taking everything into account, we have a complexity of O(T 3 ) + O(T 3 log T ) + O(T 2 ). Without loss of the generality, the overall heuristic algorithm terminates in O(T 3 log T ).

  functions are time-invariant and the second (third) pattern one with linear growth (decline). The fourth and fifth patterns represent a planning horizon with a single season of peak demands either at the beginning or at the end of the planning horizon. The last pattern (VI) is cyclical with a cycle length of six periods, such that demands in the two middle periods of each cycle are 7 times their value in the first and last period, while β t = 1 in the remaining two periods of the cycle. Using all combinations of the 6 demand patterns and the 3 TBO levels, we generate a number (18) of test problems. An average demand d = 50 is assumed. A group of replicable examples are first provided to show the differences between optimal and heuristic solutions on capacity, costs and setup numbers. We use a t = 15 and h t = 5 and analyse three different setup cost levels considering the assumption of no speculative inventory in firms. In addition, in order to calculate the capacity acquisition cost, we choose constants Λ = 200 and θ = 1. We determine the fixed setup cost indirectly by first choosing the EOQ-cycle time "Time-between-Orders (TBO)" = 2f hd and determine the setup cost value f from this identity. Let average demand d = 50, and then

For

  example, if a only a discrete number of potential capacity choices are under consideration, the problem becomes as easy as solving a capacitated lot-sizing problem for each of these choices. If cost is piecewise linear, one would need to solve for each of the intervals, compare the solution to the start and end value of the interval and pick the best solution overall. While this study effectively solves the capacity-acquisition and lot-sizing problem, it is based on deterministic demand and constant capacity assumptions. A promising avenue of future research would be taking into consideration demand uncertainty and time-varying capacity. Under stochastic demands, the phenomenon of demand shortage as well as service level constraint should be addressed as well. It might also be fruitful to analyze the problem and develop algorithms based on nonlinear production and inventory cost, or time varying fixed setup costs. In addition, as a future direction of research it would be very interesting to look into the capacity-acquisition and lot-sizing problem in a competitive environment.

Table 1 :

 1 Comparison of heuristic and pseudo-optimal solutions for our standard test problems ] for lowTBO values, the interval [2,6] for mediumTBO values and [5, 10] for high TBO values. Additionally, a random factor ε t ∼ U [0.5, 1.5] is added to the demand, and thus, d t = β t * ( d) * ε t . Maintain the cost data same as described above and run the code iteratively (i.e. 10 iterations). There are multiple instances for each TBO and demand pattern combination generated and solved. We calculate the average gaps between the heuristic and optimal solutions. The results are presented in Table2.

	Test		Demand Pseudo-optimal Solution	Heuristic Solution
	Problem	TBO	Pattern	CostSetup Capacity	CostSetupCapacity	Gap
	(1)	(2)	(3)	(4) (5)	(6)	(7) (8) (9)	(10)
	1		DP1	80000	54	50	80000	54 50	0.00%
	2		DP2	86155	38	80	87324	37 83	1.36%
	3	Low	DP3	87529	36	88	89621	44 89	2.39%
	4		DP4	125981 25	166	127893 23 166	1.52%
	5		DP5	118319 28	122	126993 29 115	7.33%
	6		DP6	83491	37	76	85581	45 63	2.50%
	7		DP1	161625 27	100	161625 27 100	0.00%
	8		DP2	163045 22	125	164315 22 123	0.78%
	9	Medium	DP3	162275 23	120	164283 23 119	1.24%
	10		DP4	175731 17	166	179342 17 168	2.06%
	11		DP5	175556 19	149	178471 16 170	1.66%
	12		DP6	163456 23	121	163611 18 151	0.09%
	13		DP1	250500 18	150	250500 18 150	0.00%
	14		DP2	251456 16	171	252621 16 170	0.46%
	15	High	DP3	249664 15	183	250772 15 181	0.44%
	16		DP4	253550 14	195	260615 12 226	2.79%
	17		DP5	253191 14	196	255037 15 181	0.73%
	18		DP6	250685 18	155	251361 18 151	0.27%

Table 2 :

 2 Average gaps between the heuristic and pseudo-optimal solutions and corresponding CPU computation times

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk

end if 12: end for 13: for i = 2 : 1 : n do for j = l i-1 : 1 : l i -1 do 17:

if V > a j + H(j, l i ) -a l i then 18:

V = a j + H(j, l i ) -a l i ; if y t = 1 then 32:

end if 35: end for The heuristic algorithm is coded using Matlab 7.5. We compare the heuristic solutions with benchmark solution as shown in Table 1. The bechmark solutions are obtained by discretising the potential capacity space, not necessarily to integer values, and evaluate the cost function z and calling CPLEX 11.0 solver in Matlab environment.

We consider discrete capacity levels C min , C min + △, C min + 2△, • • • , C max , where △ is assumed to be integer such as 1, 2, • • • , and C min = C T min , C max is the minimum capacity level to allow the optimal uncapacitated lot sizing solutions. The problems with discretised capacity values which are a series of capacitated lot sizing problems can be solved by polynomial time algorithm and standard MIP solvers. Since the pseudopolynomial algorithm is not the focus in this paper, we simply use CPLEX to solve the individual capacitated lot sizing problems to optimality here. Upon obtaining the piecewise-linear functions for each discretized capacity level, the optimal capacity and lot sizing strategy are obtained when the total capacity acquisition and lot sizing cost reaches minimal. The problem instances are solved on a Pentium 4 PC with 1G RAM, and the solution is pseudo-optimal.

As mentioned before, for the special case of quadratic capacity acquisition cost function, the problem P could in theory be handled by Mixed Integer Quadratic Programming (MIQP) solver in CPLEX. According to ILOG (2007), even relatively small integer programming models still take enormous amounts of computing time to solve and it is a very common occurrence with MIPs that the programme runs out off memory. Indeed this is what we have experienced when we attempted to solve instances by the CPLEX MIQP solver. All instances run in excess computation time; CPLEX ran out off memory in the majority of cases. For those instances that reached optimality, the gap was nearly zero when compared to our discretisation, which needed far less computational time.

Comparing with the pseudo-optimal solutions, the heuristic solutions are reasonable since the heuristic algorithm also suggests similar setup numbers and capacity levels (see Table 1). However, the differences between the solutions are relevant with the demand pattern. For the constant demand pattern (DP1), the heuristic algorithm provides the optimal solution. Readers may notice that, for test problem 5 in Table 1, the gap between the optimal cost and heuristic cost is oddly big comparing with other test problems. Since one specific example is not representative, we need to find the average performance of the heuristics algorithm.

We add the randomness of the demand and fixed setup costs from the hypothetical examples above. The TBO value is generated from a uniform distribution on the interval