B Rivet

H Cecotti

E Maby

J Mattout

Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface

Keywords: Brain-Computer Interface, P300, EEG, BLDA, Sensor selection

A challenge in designing a Brain-Computer Interface (BCI) is the choice of the channels, e.g. the most relevant sensors. Although a setup with many sensors can be more efficient for the detection of Event-Related Potential (ERP) like the P300, it is relevant to consider only a low number of sensors for a commercial or clinical BCI application. Indeed, a reduced number of sensors can naturally increase the user comfort by reducing the time required for the installation of the EEG (electroencephalogram) cap and can decrease the price of the device. In this study, the influence of spatial filtering during the process of sensor selection is addressed. Three sensor selection approaches based on a recursive backward elimination are compared. Two of them maximize the Signal to Signal-plus-Noise Ratio (SSNR) for the different sensor subsets while the third one maximizes the differences between the averaged P300 waveform and the non P300 waveform. We show that the locations of the most relevant sensors subsets for the detection of the P300 are highly dependent on the use of spatial filtering. Applied on data from 20 healthy subjects, this study proves that subsets obtained where sensors are suppressed in relation to their individual SSNR are less efficient than when sensors are suppressed in relation to their contribution once the different selected sensors are combined for enhancing the signal. In other words, it highlights the difference between estimating the P300 projection on the scalp and evaluating the more efficient sensor subsets for a P300-BCI. Finally, this study explores the issue of channel commonality across subjects. The results support the conclusion that spatial filters during sensor selection procedure allow selecting better sensors for a visual P300 Brain-Computer Interface.

Introduction

A Brain-computer interface (BCI) enables a user to communicate through the direct and real-time measurements of brain activity [START_REF] Allison | Brain-computer interface systems: progress and prospects[END_REF]. For patients suffering from severe motor disabilities, e.g. locked-in patients, a BCI may often represent the only mean of communication between the patients and their family. BCIs may actually provide the only communication pathway for patients who are unable to communicate via conventional means [START_REF] Birbaumer | Brain-computer interfaces: communication and restoration of movement in paralysis[END_REF]. Although a remaining challenge is to get an efficient working BCI, other challenges aim at reducing the cost of a BCI and to bring BCI technology to mass market, to healthy persons for new BCI applications [START_REF] Guger | How many people are able to operate an EEG-based brain-computer interface (BCI)?[END_REF][START_REF] Guger | How many people are able to control a P300-based brain.computer interface (BCI[END_REF]. This necessary step can require the reduction of the number of electrodes and the adaptation of their locations for each user [START_REF] Sannelli | On optimal channel configurations for SMR-based brain-computer interfaces[END_REF]. Indeed, reducing the number of sensors yields more comfort for the user, decreases installation time duration and may substantially reduce the financial cost of the BCI setup since the cost of an EEG cap and an amplifier vary in relation to the number of channels. Besides, the reduction of the number of sensors can also reduce the power consumption for wireless EEG caps [START_REF] Shih | Sensor selection for energy-efficient ambulatory medical monitoring[END_REF]. This point is relevant for BCI applications like video games or for other applications where the number of sensors is limited. Finally sensor selection may facilitate the classifier work by selecting a reduced and more relevant set of input features.

Several strategies are possible for selecting a subset of sensors. For instance, one can select sensors based on prior knowledge from neuroscience studies that can provide the source of the potentials. In such a case, the subset is fixed and may jeopardize the performance in some subjects as the optimal sensor subset is subject-dependent [START_REF] Hoffmann | An efficient P300-based brain-computer interface for disabled subjects[END_REF].

To obtain optimal performance, it is mandatory to clearly identify subject-specific optimal sensor subsets. For a given set of N sensors, there are 2 N -1 different possible subsets, without the empty subset. To find the optimal subset, three main searching Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface3 approaches can be considered: exhaustive, random or sequential search. The exhaustive search is usually not possible: the search space is often exponentially prohibitive. The sequential and random search does not guarantee optimality as it gives up completeness. However, it is easily implementable and can provide suboptimal solutions given some evaluation criterion. Several variations are described in the literature like the greedy hillclimbing approach, forward selection, backward elimination or bi-directional selection.

For instance, recursive feature elimination was used for sensor selection in BCI based on motor imagery [START_REF] Lal | Support vector channel selection in BCI[END_REF][START_REF] Schroder | Robust EEG channel selection across subjects for braincomputer interfaces[END_REF]. In this paper, we consider a recursive backward elimination strategy to determine which sensor should be removed in the widely used P300 speller paradigm [START_REF] Donchin | Assessing the speed of a P300-based braincomputer interface[END_REF]. The impact of spatial filtering during the selection of the more relevant sensor subsets is then evaluated with the P300 speller accuracy. Moreover an analysis of commonality of selected channels across subjects is studied. The remainder of the paper is organized as follows: the sensor selection strategy and the different methods are described in the second section; the materials and the experimental protocol are presented in the third section; and finally, results are compared and discussed in the last section.

Methods

Data modeling: xDAWN framework

The proposed method to select adaptively a relevant subset of sensors is intimately related to the special structure shared by the data due to the oddball paradigm used in a visual P300 BCI. Indeed, the xDAWN framework [START_REF] Rivet | xDAWN algorithm to enhance evoked potentials: application to brain-computer interface[END_REF] provides an algebraic model of the recorded signals X ∈ R Nt×Ns (where N t is the number of time samples in the recorded EEG signal and N s the number of sensors) by considering that X is composed of three additive terms: (i) the P300 responses A 1 related to target stimuli, (ii) responses

Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface4

A 2 common to every stimuli, and (iii) the residual noise H. As a consequence, one can write X as

X = D 1 A 1 + D 2 A 2 + H, (1)
where A 1 ∈ R N 1 ×Ns and A 2 ∈ R N 2 ×Ns are the patterns representing one target evoked potential (i.e. the P300 response) whose length is N 1 samples and one common evoked potential whose length is N 2 samples, respectively.

D 1 ∈ R Nt×N 1 (resp. D 2 ∈ R Nt×N 2 )
is a Toeplitz matrix whose first column entries are set to zero except for those that correspond to a target (resp. every) onset which are equal to one. D 1 A 1 and D 2 A 2 thus represent the superposition of all targets and common responses, respectively. This parametrization takes into account the possible overlapping between target potentials and/or between common potentials and it allows to estimate these responses in the least mean square sense by

    Â1 Â2     = ([D 1 ; D 2 ] T [D 1 ; D 2 ]) -1 [D 1 ; D 2 ] T X (2)
where [D

(e i , X) = e T i ÂT 1 D T 1 D 1 Â1 e i e T i X T Xe i , ( 3
)
where e i ∈ R Ns is a zero vector but the i-th element equal to 1. Y e i is thus the i-th column of a matrix Y .

The purpose of applying spatial filter u ∈ R Ns is to create a virtual sensor so that the resulting signal Xu provides a larger SSNR ρ(u, X) than SSNR achieved by actual sensors max i ρ(e i , X). xDAWN framework [START_REF] Rivet | xDAWN algorithm to enhance evoked potentials: application to brain-computer interface[END_REF] also provides a way to estimate

N f spatial filters u i , i ∈ {1, • • • , N f }, that maximize the SSNR (3).
Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface5

Sensor selection

The chosen method for adaptively selecting a relevant subset of sensors is based on a recursive backward elimination. Starting with all sensors, it consists in alternatively testing each sensor for its significance and in removing the least relevant one at each iteration step. An irrelevant sensor is a sensor whose removal barely impairs the performance or selection criterion. In this work, we eliminate two sensors at a time, leaving us with the most significant remaining subset. Elimination goes on until every sensor has been eliminated. At the end of the selection process, sensors can be ranked according to their revealed significance. A relevant sensor will be eliminated at the end of the iterative process while a useless one will be eliminated along the very first iterations.

The cost function that determines the relevance of a sensor subset is based on the evaluation of the SSNR. The cost function of the sensor subset I can be estimated as a function of the individual SSNR ρ(e i , X) of each actual sensor

s a (I) = i∈I ρ(e i , X). ( 4
)
This solution is denoted AS as actual sensor. Moreover, the SSNR can also be estimated after some spatial filtering. In this case, the score function is estimated globally as a function of the grouped SSNR from the virtual sensors

s v (I) = N f k=1 ρ(u k , XE I ), (5)
where E I is the matrix whose columns are the concatenation of vectors e i with i ∈ I.

This technique is referred as V S as virtual sensor.

As a reference, the score function based on the largest difference between the averaged P300 trials waveform and the averaged non P300 trials waveform is finally considered and this method is denoted (Ref). In this case the score function of the

Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface6

AS: SSNR V S: Spatial Filters (SF) + SSNR Ref: Largest difference between P300 and Non P300 trials Table 1. The three criteria for backward selecting the most relevant sensors.

sensor subset I is estimated as

s r (I) = i∈I (x P 300 i -x N P 300 i ) T (x P 300 i -x N P 300 i ), (6)
where x P 300 i (resp. x N P 300 i

) is the vector of the averaged P300 (resp. non P300) trials waveform at channel i.

Finally, the backward elimination procedure is summarized in Algorithm 1 where the considered score functions are presented in Table 1. ρ(e i , X), and method V S that selects the most relevant sensors in relation to their grouped SSNR ρ(u, XE I ).

Materials

The EEG signal was recorded on 20 healthy subjects (mean age= 26, sd= 5.7; 7 females) with the OpenViBE framework [START_REF] Maby | The OpenViBE P300-Speller scenario: a thorough online evaluation[END_REF]. Each subject filled a consent form prior to the experiment. Subjects were wearing a standard EEG cap with 32 electrodes subsampled from the 10-10 system [START_REF] Sharbrough | Guidelines for standard electrode position nomenclature[END_REF]. The EEG activity was recorded continuously from 32 active electrodes (actiCap, Brain Products GmbH, Munich). The electrodes for the ground and the reference were placed on the forehead and the nose, respectively.

For the evaluation of the spatial filters impact during sensor selection, two sessions were considered. The P300 speller was presented to the subjects as described in [START_REF] Farwell | Talking off the top of your head: toward a mental prosthesis utilizing event-related brain potentials[END_REF],

with a square matrix of 36 symbols, where one symbol is a letter or a digit. The EEG signal was sampled at 100Hz and was bandpass filtered between 1Hz and 12.5Hz with a Butterworth filter (order=4). Finally the features vector is obtained by a downsampling to 25Hz for speeding up the computation time and reducing the number of features. Such decimation parameter has been successfully applied in other studies about the P300 speller [START_REF] Cecotti | Convolutional neural networks for P300 detection with application to brain-computer interfaces[END_REF][START_REF] Rakotomamonjy | BCI competition iii : Dataset ii -ensemble of SVMs for BCI P300 speller[END_REF]. For each sensor and for each trial, the signals were normalized so that they had a zero mean and a standard deviation equal to one. A Bayesian linear discriminant analysis (BLDA) classifier is considered for the detection of the P300 [START_REF] Hoffmann | An efficient P300-based brain-computer interface for disabled subjects[END_REF][START_REF] Mackay | Bayesian interpolation[END_REF]. ‡ Note that in this study, the different ISI are not exploited as a parameter.

Results

The P300 speller was evaluated with different sensor subsets as defined in the previous section. For each selected sensor subset, the signal was first enhanced by using spatial filters before the classification stage performed by a BLDA for the detection of the P300.

For the classification, a sequence of 0.6s from a visual stimulus has been considered for each trial. Therefore, the P300 response is expected to be centered on the chosen time sequence. For each subject, the number of P300 trials and non P300 trials are 1200 and 6000, respectively, for the test session. Based on a previous study [START_REF] Rivet | Sensor selection for P300 speller brain computer interface[END_REF], we only consider a subset of four selected sensors by the different methods (AS, V S and Ref) since it has been shown that four sensors is a good tradeoff between the number of sensors and the classification accuracy. The impact of spatial filters during the selection process is firstly assessed (Section 4.1). The commonality of channels across subjects is then studied (Section 4.2).

Impact of spatial filters during sensor selection

Figure 1 presents the recognition rate of the P300 speller, with 10 repetitions, for the subsets of four sensors obtained with the three methods (AS, V S, Ref ), and for the whole set of 32 sensors (All). The average recognition rate is 80.4%, 76.4%, 90.3% and 95.8% for the subsets obtained with Ref , AS, V S and the whole set (All), respectively.

Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface9

The performance of the P300 speller is always better (excepted for subject 2) with the four sensors selected with method V S (i.e. with spatial filters) compared to AS or to Ref . Two-tailed t-tests (5% significance level) indicate that the method V S provides better sensors than AS (t (19) = 4.42, sd = 8.4, p < 5e -4 ) and than Ref (t (19) = 4.10, sd = 6.49, p < 1e -3 ). A pairwise two-tailed t-test has also been used to compare the speller accuracy with sensors obtained with V S and the whole set. As expected, the result of the t-test indicates a rejection of the null hypothesis at the 5% significance level. It shows that the whole set of 32 sensors (All) provides a better performance than the 4 sensors obtained with V S (t (19) = 2.62, sd = 5.36, p < 2e -2 ).

The position of the best sensors selected with the methods AS and V S is depicted for every subject in Figure 2. The grey color represents the sensor relevance, i.e. a dark grey means the location is relevant. For both methods, the sensor subsets are not homogeneously located across subjects, highlighting the variability across subjects.

In a general way, V S provides sparser sensors than AS. Nevertheless, as depicted in Figure 3, the grand average of the evoked response obtained after the spatial filtering step of the four sensors selected by AS and V S show that the V S method provides a better defined evoked response, i.e. a less large P300 wave, than the AS method. Note that the N200 component is also visible before the P300 peak for both AS and V S. 5) depicts the averaged P300 and non P300 waves across all subjects (resp. for Subject 20) and all sessions. This latter subject is a good example to provide evidence that electrodes with a relative low amplitude can be relevant in a sensor selection procedure. Indeed, one may directly point out the projection of the P300 over P Z as the amplitude of the P300 is higher than with other sensors. However, this electrode was not selected in the best sensor subsets as presented in Figure 1, (S20 (AS) and S20 (VS)). Instead, electrodes closer to the bilateral occipito-temporal parts Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface10 were selected.

S1 (AS) S1 (V S) S2 (AS) S2 (V S) S3 (AS) S3 (V S) S4 (AS) S4 (V S) S5 (AS) S5 (V S) S6 (AS) S6 (V S) S7 (AS) S7 (V S) S8 (AS) S8 (V S) S9 (AS) S9 (V S) S10 (AS) S10 (V S) S11 (AS) S11 (V S) S12 (AS) S12 (V S) S13 (AS) S13 (V S) S14 (AS) S14 (V S) S15 (AS) S15 (V S) S16 (AS) S16 (V S) S17 (AS) S17 (V S) S18 (AS) S18 (V S) S19 (AS) S19 (V S) S20 (AS) S20 (V S)

Commonality of channels across subjects

In Figure 6, an analysis of the differences and similarities of selected sensors for each subject is provided. The mean, the standard deviation (S.D.) and the kurtosis across the 20 subjects with four selected sensors are depicted in Fig. 6 AS and V S: mainly in the occipital and parietal areas. However, it is challenging to extract some common sensors that may be useful for every subject. The S.D. is higher in the occipital and parietal areas, suggesting that the order of these sensors can vary a lot while keeping a high importance in the sensor selection. Finally, the figures of the last column aim at depicting the sensors that could be relevant for few subjects. The kurtosis of the sensor rank suggests that although the majority of the relevant sensors are in the occipital and parietal area, some sensors in the frontal area remain important for some subjects to enhance the P300.

Even if this last observation proves the necessity to customize the sensor location for each subject so that the performance are maximized, we compare two criteria to select a common subset of four sensors for all the subjects (Table 2 and Figure 7). The

Conclusion and discussion

Three ways for the analysis of the best sensor location for P300 detection have been proposed. In every solution, the method is based on a recursive backward elimination.

In This study allowed comparing different ways for analyzing the location of the P300 on the scalp. The different experiments have shown that it is possible to achieve relatively good performance with only four sensors, with an accuracy of about 90% with 10 epochs. For the different subsets of four sensors, which are personalized to each subject, two sensors are common to half of the subjects with AS (P O 9 , O 1 ). With V S, only P 8 is common to more than half of the subjects (Table 2). Some sensors are located on the occipital area, confirming previous works suggesting that occipital sites have also an important role [START_REF] Bianchi | Which physiological components are more suitable for visual ERP based brain-computer interface? a preliminary MEG/EEG study[END_REF][START_REF] Cecotti | Convolutional neural networks for P300 detection with application to brain-computer interfaces[END_REF][START_REF] Krusienski | Toward enhanced P300 speller performance[END_REF].

Although the P300 response has been discovered for over 40 years [START_REF] Sutton | Evoked potential correlates of stimulus uncertainty[END_REF], its full understanding remains challenging and elusive [START_REF] Polich | Updating P300: An integrative theory of P3a and P3b[END_REF]. The best locations for the sensors depend on the person and highlight both the complexity underlying the P300 process and the need of a personalized/adaptive P300-BCI with a limited number of sensors. Several studies in the literature have located the sources of the P300. Maeno et al. [START_REF] Maeno | Source modeling of the P300 eventrelated response using magnetoencephalography and electroencephalography measurements[END_REF] point out the P300 sources on the bilateral occipito-temporal gyrus for visual stimuli as this area lies in the ventral pathway for vision. In [START_REF] Yamaguchii | Anterior and posterior association cortex contributions to the somatosensory P300[END_REF], Yamaguchii et al. indicate that association cortex in the temporal-parietal junction is critical for generating the scalp-recorded target and novelty P300s, whereas dorsolateral frontal cortex contributes preferentially to novelty P300 generation. These sources are projected on the scalp and can explain the topographic maps presented in Figure 2. The inter-individual differences may be ascribed to the peculiar convolution of individual cortex with heterogeneous equivalent resulting dipole orientation [START_REF] Tarkka | Source localization of P300 from oddball, single stimulus, and omitted-stimulus paradigms[END_REF]. For the P300 speller paradigm, the type of P300 detected corresponds mainly to the P3b, which is related to context updating operations and subsequent memory storage. As the P3b is generated via parieto-temporal cortical structures, it could explain the high rank of P 8 and P 7 for most of the subjects. Future

Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface17 works will deal with the sensor selection in relation

Algorithm 1

 1 Backward elimination to sensor selection. Initialize I = {1, • • • , N s } for i = 1 to N s do for all j ∈ I do Compute score function s . (I \ j) of I \ j = {i ∈ I|i = j} end for Estimate less relevant sensor ĵ = arg max j s . (I \ j) Remove less relevant sensor: I ← I \ ĵ end for The objectives of the experiments are to evaluate and to compare the impact of the spatial filters during the selection process, i.e. the difference between method AS which eliminates the individual less relevant sensors in relation to their individual SSNR

  The first session was dedicated to the estimation of different subsets of sensors and for training the classifier. The second session was for testing the P300 speller. In both sessions, the number of repetitions was fixed to 10, i.e. each symbol is spelt 10 consecutive times. The first and second sessions contained 50 and 60 symbols, respectively. The Interval Stimulus Interval (ISI) was 170ms for the training session. For the test session, four groups of 15 characters had an ISI of 110, 150, 190 and 230ms, respectively ‡. The average ISI in the test session was like the training session, i.e. 170ms.

Figure 1 .

 1 Figure 1. Impact of the sensor subsets on the recognition rate of the P300 speller.

Figure 4 (

 4 Figure 4 (resp. Figure 5) depicts the averaged P300 and non P300 waves across all

Figure

  Figure 4 (resp. Figure 5) depicts the averaged P300 and non P300 waves across all

Figure 2 .

 2 Figure 2. Position of the 4 best sensors for the criteria AS and VS. A dark/light grey represents a relevant/non relevant sensor location.

Figure 3 .

 3 Figure 3. Averaged P300 and non P300 trials after spatial filtering of 4 sensors selected with AS and V S, during 600ms.

Figure 4 .

 4 Figure 4. Averaged ERP on the target and non-target across the different electrodes. The bold and thin lines represent the ERP on the target and non target, respectively.

Figure 5 .

 5 Figure 5. Averaged ERP on the target and non-target across the different electrodes for Subject 20, during 600ms. The vertical dash line corresponds to 300ms after a visual stimulus. The bold/thin line represents the averaged P300/non P300 response.

Mean

  first one (denoted 'Voting') is simply obtained by a voting procedure: the sensors are ranked in relation to their frequency of appearance in the four best sensors for the 20 Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface13

Figure 6 .

 6 Figure 6. Global rank for every sensor. The rank represents the step during the backward elimination when the sensors has been removed. (a) Voting selection

  (a) and by a joint selection (b). Each cell represents the electrode position in the international extended 10-10 system. For the voting procedure (a), the frequency this electrode has been selected as one of the four best sensors across the 20 subjects is also indicated.subjects. The second procedure (denoted 'Joint') refers to the backward elimination procedure (Algorithm 1) where the subject specific score function s . (I \ j) is replaced by the average of the 20 individual score functions related to each subject.

  the first case, sensors are suppressed in relation to their individual SSNR. In the second case, sensors are suppressed in relation to their contribution once the different selected sensors are combined for enhancing the signal. In the third solution, sensors are suppressed in relation to the largest differences between P300 and Non P300 trials waveforms. The results indicate that the sensors, which are contained in the best sensor subset are different than the sensors that provide individually the best SSNR. In other words, it proves that given a low number of sensors, the sensors contained in the best subset are not individually the best sensors. This observation shows that the sensor subset should be consider directly as a whole and not as the aggregation of individual sensors in spite of the backward elimination procedure. Sensors should complement each other in the subset to enhance the P300 detection. This highlights the difference between two problems: evaluating the best sensor subsets for a P300-BCI, and estimating the best P300 projection on the scalp.Finally, the sensor selection method is not based on the detection of the P300 but only on the estimation of the SSNR, thus avoiding different stages like classifiers training and cross validation. Spatial filtering is not only a necessary pre-processing Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface16 step for signal detection, it is also a relevant pre-processing step for sensor selection.

  1 ; D 2 ] is a matrix of size N t × (N 1 + N 2 ) obtained by concatenation of D 1 and

	D 2 . xDAWN model (1) and estimation (2) are useful to define the Signal to Signal plus
	Noise Ratio (SSNR) of the i-th sensor by
	ρ

Table 2 .

 2 ; 11 O 1 ; 10 P O 10 ; 9 P 8 ; 8 P 7 ; 7 P Z ; 5 P 3 ; 4 T P 10 ; 3 P 4 ; 3 Ref O 1 P 8 P 3 P O 9 P 7 O 2 P Z P 4 P O 10 O Z Top 10 sensors for AS and V S across 20 subjects selecting by a voting procedure

	Top	1	2		3	4	5		6	7	8	9	10
	AS P O 9 O 2 ; 3
	V S	P 8 ; 12	P 7 ; 8		P 3 ; 8	O Z ; 8 P Z ; 7 P O 9 ; 5 O 1 ; 4	T 8 ; 3	O 2 ; 3 P O 10 ; 3
	Ref O 1 ; 12 P O 9 ; 8	O 2 ; 8	P 3 ; 6 P 8 ; 6 P 4 ; 5 O (b) Joint selection P 7 ; 5
		Top	1	2	3	4	5	6	7	8	9	10
		AS	P 8 P 7 P O 9 O 1	P 3 P O 10 O 2 P 4	O z	P z
		V S P O 9 P 3 C Z	P 8 O Z	T 8	F 3 O 1	P 4	P Z

Z ; 4 P Z ; 4 P O 10 ; 4

  Comparison of classification accuracy achieved for all the 20 subjects by several selection procedures (Ref , AS and V S). 'Ind.', 'Joint' and 'Vot.' refer to a subject specific selection, to the joint selection and to the voting selection, respectively. All refers to the classification accuracy achieved with all the 32 sensors. On each box, the central mark is the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to the most extreme data points not considered outliers, and outliers are plotted individually.As already mentioned, there is a large diversity across subjects about the location of the most relevant sensors. Indeed, the largest value of the frequency of appearance among the four most relevant sensors is only 11 or 12 depending of the selection criterion: this means that this sensor is only common to about half of the subjects. It is also worth noting that the voting and joint selection procedures provide quite different selected sensors. Actually, for the V S method, if P 8 and P 3 are among the four selected sensors for both voting and joint selection methods, P O 9 is only ranked as the sixth most relevant method by the voting selection whereas it is the most relevant sensor for the joint selection. Even more surprising, C Z (resp. P 7 ) is in the four most relevant sensor for the joint (resp. voting) method while it is not among the 10 most relevant sensors for the voting (resp. joint) selection. Furthermore, as shown by Figure7, among the joint and voting selection, the joint and voting selections with virtual sensor (Joint V S and V ot. V S) provide the best classification accuracy. Indeed, with these four sensors common to all the 20 subjects selected by Joint V S (resp. V ot. V S), the classification accuracy is larger than 60% for 17 (resp. 16) of them and 16 (resp. 14) subjects achieved even more than 75% of classification accuracy. Moreover eleven subjects achieved an Impact of spatial filters during sensor selection in a visual P300 Brain-Computer Interface15 impressive classification accuracy larger than 90% with only four sensors selected by V ot. V S. Finally, this latter result highlights one more time the fact that using spatial filters during the sensor selection is essential to provide accurate subsets of channels classification accuracy for both Joint V S and V ot. V S is larger than with AS or

			100
	Classification	accuracy [%]	40 60 80
			20
			0
			Ind. Ref Ind. AS Ind. VS Joint Ref Joint AS Joint VS Vot. Ref Vot. AS Vot. VS	All
	Figure 7. Ref methods.

Acknowledgment

The authors thank the anonymous reviewers for their pertinent comments which allow to drastically improve the quality of this paper. This work has been supported by French National Research Agency (ANR) through TecSan program (project RoBIK ANR-09-TECS-013) and through DEFIS program (project Co-Adapt ANR-09-EMER-002).


Christoph H Glock
email: christoph.glock@uni-wuerzburg.de

Batch sizing with controllable production rates in a multistage production system

Keywords: INVENTORY CONTROL, BATCH SIZING, PRODUCTION variable production rate, multi-stage lot sizing inventory control, production, variable production rate, batch sizing, multi-stage lot-sizing

Batch sizing with controllable production rates in a multi-stage production system

Christoph H.

Introduction

The control of inventories in production and supply chain systems has frequently been the subject of interest for researchers in recent years. Reducing excessive inventory in a supply chain may increase the flexibility of the chain and simultaneously reduce the costs associated with producing, storing, and transporting products (see for example [START_REF] Lee | The Bullwhip Effect in Supply Chains[END_REF]). While early works on inventory control mainly focused on production systems with a single production and a single consumption stage, more recent publications have analysed interdependencies between multiple stages of a production system (see e.g. [START_REF] Bogaschewsky | Optimizing multi-stage production with constant lot size and varying number of unequal sized batches[END_REF][START_REF] Chiu | A comprehensive review of lot streaming[END_REF][START_REF] Glock | Supply chain coordination via integrated inventory models: a review[END_REF]).

In a recent paper, Glock (2010a) studied the impact of variable production rates on the buildup of inventory in a two-stage production system and showed that deviating from the 'design production rate' of the system may reduce the system's total costs. This is important, as varying the production rate gives production planners flexibility in smoothing material flows and in avoiding the accumulation of inventory at bottleneck stations. The findings of Glock (2010a), however, are not applicable to multi-stage production systems because of the interdependencies between consecutive stages that are neglected in inventory models which consider only two stages. This paper extends the model of Glock (2010a) to a multi-stage production system where the production rates at the stages may be varied between minimum and maximum values. This paper relates two streams of research, which are: (1) multi-stage production systems and

(2) variable production rates in inventory models. Studies that fall along the first stream analysed the impact of different transportation strategies on the system's total costs. For example, [START_REF] Szendrovits | Manufacturing cycle time determination for a multi-stage economic production quantity model[END_REF] presented a multi-stage production system and assumed that batches of equal sizes are transferred between subsequent stages. He showed that the overlapping of production cycles at subsequent production stages may reduce inventory in the system and lower total costs. [START_REF] Goyal | Note on manufacturing cycle time determination for a multi-stage economic production quantity model[END_REF] extended the model of Szendrovits by considering a fixed transportation cost to ship a batch to the next stage, and proposed a solution procedure.

Assuming that the number of batches may differ between production stages, Goyal (1977a) and [START_REF] Szendrovits | Optimizing Multi-Stage Production with Constant Lot Size and Varying Numbers of Batches[END_REF] developed solutions for models with equal-sized batch shipments. Drezner et al. (1984) showed that equal-sized batch shipments reduce the system's total costs. [START_REF] Bogaschewsky | Optimizing multi-stage production with constant lot size and varying number of unequal sized batches[END_REF] The studies of the second stream of research assumed that the mean output per unit of time of a production system may be influenced either by varying the rate of performance of the production equipment or by inserting idle times between task elements (see for example [START_REF] Buzacott | One-and Two-Stage Scheduling of Two Products with Distributed Inserted Idle Time: The Benefits of a Controllable Production Rate[END_REF][START_REF] Schweitzer | Optimizing processing rates for flexible manufacturing systems[END_REF]. In the second case, the time which is effectively spent on manufacturing is reduced. In this context, some studies assumed that the production rate may be changed during the production process, while other studies assumed that the production rate has to be fixed before a production run starts. In the first case, which is commonly referred to as the 'flexible case', the setup costs are only incurred when switching to producing a different product and not when adjusting the speed of the production equipment. In the second case, which is termed the 'rigid case', changing the setup of a machine during a production run is not possible either because it is technically impossible or associated with prohibitively high costs (Glock 2010a). This paper combines both research streams by studying the effects of variable production rates on a multi-stage inventory system. The impacts of the 'rigid case' and the 'flexible case' of the production system of concern are investigated for equal and unequal batch policies. The remainder of the paper is organised as follows: In the next section, the paper describes the problem and outlines the assumptions and definitions that will be used in the remaining sections of the paper. Accordingly, we develop formal models for the rigid and the flexible cases for equal-and unequal-sized batch shipment policies and propose solution procedures for the models. Section 4 contains a numerical study and section 5 concludes the paper.

Problem description

In the following section, we study the case where a single product is manufactured in a serial production system consisting of multiple stages and assume that batch shipments are transferred between subsequent stages. Production systems as the one described in this paper can be found in a variety of application areas, for example in the automotive industry or the electronic industry (see for example [START_REF] Kuo | Bottlenecks in serial production lines: A system-theoretic approach[END_REF][START_REF] Cochran | Optimum junction point location and inventory levels in serial hybrid push/pull production systems[END_REF][START_REF] Chiang | DT-Bottlenecks in Serial Production Lines: Theory and Application[END_REF]. As in Glock (2010a), we differentiate between two types of batch transfer policies that are discussed in [START_REF] Szendrovits | Manufacturing cycle time determination for a multi-stage economic production quantity model[END_REF] and [START_REF] Goyal | Determination of Optimum Production Quantity for a Two-Stage Production System[END_REF]. The first policy suggests transferring batches of equal sizes between subsequent stages, which allows for an overlap between production and consumption stages thus reducing the manufacturing cycle time and the inventory levels in the system. The second policy suggests transferring batches of increasing or decreasing sizes between subsequent stages. The increase or decrease in the batch size is constant and equal to the ratio of the production rate to the demand rate. Both policies have been studied in multi-stage production systems (see for example [START_REF] Szendrovits | Manufacturing cycle time determination for a multi-stage economic production quantity model[END_REF][START_REF] Bogaschewsky | Optimizing multi-stage production with constant lot size and varying number of unequal sized batches[END_REF] and will be investigated for variable production rates in the following section.

The following assumptions and notations are used to formulate the problem:

Assumptions:

1. The model's input parameters are deterministic and constant over time.

2. A uniform lot of size Q is manufactured at each stage.

3. The production rate at each stage can be varied between given limits p s,min and p s,max .

4. The minimum production rate at each stage exceeds the demand rate and p s,min > d.

5. The unit production costs are a function of the production rate.

6. The unit production cost function is convex in p s and p s 0 is the 'design production rate' that minimises unit production costs at stage s.

7. Batches are only sent to the subsequent stage when the inventory level at this stage reaches zero.

8. The number of batches is identical for each stage.

9. We explicitly focus on a time weighted cost function to enhance the practical applicability of our model. As [START_REF] Higgins | Manufacturing Planning and Control[END_REF] note, production plans are often made with planning horizons rolling through time as time progresses. Thus, companies who wish to implement our model as a heuristic planning tool may specify a planning period and repeat the planning process after the planning time has elapsed. Note that this assumption is not uncommon in the literature, cf. for example [START_REF] Bogaschewsky | Optimizing multi-stage production with constant lot size and varying number of unequal sized batches[END_REF] and Glock (2010a). However, normalizing the planning period to 1 time unit results in the common money-per-unit-time cost function.

10. Shortages are not allowed.

Model formulation and solution

The rigid case

a) Equal-sized batch shipments

The first scenario we analyse is the case where equal-sized batch shipments are transported to the subsequent stage and where the production rate may be varied only before the start of a production run. The corresponding inventory time plots are illustrated exemplarily in figure 1 (cf. the bold lines. Note that the dashed lines represent the case of a reduced production rate).

It is shown that a lot of size Q is manufactured at stage s and that m batches of size Q/m are transported to the subsequent stage whenever the inventory level in front of the stage drops to zero. Stage s+1 consumes the products delivered by stage s and produces a product which is then stored and successively shipped batch-wise to stage s+2, which is the final stage in the example. Note that our model is valid for an arbitrary number of stages and is not restricted to three stages as in figure 1. In the example, it can be seen that due to p s+1 > p s+2 , the second wherefore it has to be kept in stock for t w,s+1,1 = t p,s+2 -t p,s+1 time units. In contrast, due to p s+1 > p s , the production time of a batch at stage s exceeds the consumption time of a batch at stage s+1, wherefore batches 1,...,m-1 have to be kept in stock for t w,s,i time units before they are sent to the subsequent stage in order to build up enough inventory to assure an uninterrupted supply of materials. Further, the example illustrates that a change in the production rate at stage s+1 may influence inventory build-up at stage s+1 as well as at the preceding and the succeeding stage (cf. the dashed lines in figure 1).

----------Figure 1 ----------As described in Glock (2010a), the time weighted inventory at stage s+1 consists of a 'regular' inventory, which corresponds to the triangles shown in the second plot in figure 1, and inventory due to waiting times. The 'regular' inventory can be expressed as follows:

(1)

        + = +1 2 1 1 2 s s reg s p p m Q I
Inventory due to waiting times can be calculated as:

(2a) ( )         - - =         - = + - = + ∑ s s m i s s wait s p p m m Q p p i m Q I 1 1 1 2 1 1 1 2 1 1 1 2 2
for p s > p s+1 and (2b)

( )         - - =         - = + - = + ∑ 1 2 1 1 1 2 2 1 1 1 2 1 1 s s m i s s wait s p p m m Q p p i m Q I for p s+1 > p s
From (2a) and (2b) a general expression can be formulated as follows:

(3)

( ) s s wait s p p m m Q I 1 1 1 2 1 2 - - = +
The time weighted inventory is given as the sum of (1) and (3). The inventory carrying costs in the planning period for all S stages thus equal:

(4)

( ) ∑ = + +         - - +         + = S s s s s s s h p p m p p m DQ IC 1 1 1 1 1 1 1 1 2
Apart from the inventory carrying costs, the system encounters production costs which vary with a varying production rate. In the following, we assume that the unit production cost function is of the form specified by [START_REF] Eiamkanchanalai | Production lot sizing with variable production rate and explicit idle capacity cost[END_REF] and used in Glock (2010a), wherefore we can conclude that: The total costs consist of the sum of (4) and ( 5) and the setup and transportation costs in the planning period, which amount to (A s +mT s )D/Q. Thus, it follows that:

(6) ( ) ( ) ( ) ∑ ∑ ∑ = = = + + + + +         - - +         + = S s s s S s s s S s s s s s s p Dc Q D mT A h p p m p p m DQ TC 1 1 1 1 1 1 1 1 1 1 2
It can be easily shown that ( 6) is convex in Q for a given p-vector and given values for m.

Thus, the optimal solution for Q is given as: 6), the total cost function reduces to:

(7) ( ) ( ) ∑ ∑ = + + =         - - +         + + = S s s s s s s S s s s h p p m p p mT A m * Q 1 1 1 1 1 1 1 1 1 2 Substituting (7) in (
(8) ( ) ( ) ( ) ∑ ∑ ∑ = = = + + + +         - - +         + = S s s s S s S s s s s s s s s p Dc mT A h p p m p p m D TC 1 1 1 1 1 1 1 1 1 1 2
As could be shown by [START_REF] Szendrovits | Optimizing Multi-Stage Production with Constant Lot Size and Varying Numbers of Batches[END_REF], ( 8) is quasi-convex in m. Thus, for a given p-vector, an optimal solution for m may be calculated with the help of the following optimality condition:

(9) TC(m*-1) ≥ TC(m*) ≤ TC(m*+1)
The cost function shown in (8) was analysed by Glock (2010a) for the case S = 2 and p s > p s+1 , who derived rules of how to vary the production rate in order to reduce total costs. However, his results may only partially be transferred to the case of a multi-stage production system due to the interdependencies that arise between multiple production stages. The basic problem associated with varying production rates in a serial production system is illustrated in figure 2, which considers three production stages and two buffer stocks which are located between each two successive stages. If we consider the case where p s-1 > p s < p s+1 , inventory accumulates in the buffer stock in front of stage s since stage s is not able to consume the products delivered by the preceding stage as fast as stage s-1 produces them, and in the buffer stock in front of stage s+1 since a certain amount of inventory has to be build up prior to the start of the production cycle at stage s+1 in order to assure an uninterrupted supply of materials. If we now assume that the production rate is increased at stage s such that p s-1 > p s > p s+1 , inventory is reduced in the first buffer stock and partially transferred to the second (depending on the former and actual ratio of p s /p s+1 ). Thus, it becomes obvious that varying the production rate at any particular stage may influence upstream and downstream inventory, which is helpful in smoothing material flows and in exploiting different inventory carrying charges at different stages in the production system. The basic problem associated with optimizing (8) in p s is that (8) is not necessarily convex in p s for given values of m. However, we may state the following theorems to derive a solution:

Theorem 1: For m = 1, producing with p s < p s 0 cannot be optimal.

Proof: see appendix 1.

Theorem 2: For m = 2, it is optimal to produce with p s = p s 0 in case p s-1 < p s > p s+1 and to increase p s in all other instances.

Proof: see appendix 1.

Theorem 3: For m > 2,

• producing with p s < p s 0 cannot be optimal in case p

s-1 > p s > p s+1 and m/(m-2) > h s /h s-1 . If p s-1 > p s > p s+1 and m/(m-2) < h s /h s-1
, producing with p s > p s 0 cannot be optimal.

• producing with p s > p s 0 cannot be optimal in case p s-1 < p s > p s+1 .

• producing with p s < p s 0 cannot be optimal in case p s-1 > p s < p s+1 .

• producing with p s < p s 0 cannot be optimal in case p s-1 < p s < p s+1 and m/(m-2) > h s-1 /h s . If p s-1 > p s > p s+1 and m/(m-2) < h s-1 /h s , producing with p s > p s 0 cannot be optimal.

Proof: see appendix 1.

To calculate a good, but not necessarily optimal value for p s , we can apply a one-dimensional search algorithm. For example, we could conduct a grid search and, starting from p min , increase p stepwise until p max is reached and select the best solution found during the search process as the final solution (see [START_REF] Gill | Practical Optimization[END_REF] for other alternatives). The following procedure can be used to find a solution for m*, p s *, and Q*:

Step Step 3: Search p s ∈ [p s,min ; p s 0 ] or p s ∈ [p s 0 ; p s,max ] that minimises TC according to the conditions specified above. If TC > TC*, Goto Step 5.

Step 4: Set TC* = TC, m* = m, p s * = p s ∀s, and m = m+1. Goto Step 3.

Step 5: Find Q* from (6).

b) Unequal-sized batch shipments

In the case where unequal-sized batches are shipped to the subsequent stage, the size of the jth batch at stage s can be calculated as follows (see Goyal 1977b):

(10) q j,s = q 1,s λ s

j-1 = q 1,s (p s /p s+1 ) j-1 with ∑ ∑ = - = = = m i i s s , m i s , i q q Q 1 1 1 1 λ
It is obvious that subsequent batches increase in size in case p s > p s+1 and decrease if p s < p s+1 .

The corresponding inventory time plots that result if t p,i,s = t p,i-1,s+1 are illustrated in figure 3 (cf. the bold lines. The dashed lines again illustrate the case of a reduced production rate). As can be seen, there is no inventory due to waiting times if subsequent batches increase or decrease by λ s = p s /p s+1 .

----------Figure 3 ----------The time weighted inventory at stage s can be calculated as follows:

(11) ( ) ( )( ) ( )(
)

1 1 1 1 1 1 2 1 1 2 2 2 2 1 2 1 2 1 1 2 1 1 1 2 2 2 1 1 1 1 - + - +         + =         + = =         + + +         + +         + = + = - + + + + ∑ m s s s m s s s m i i s s s s , s s , m s s , m s , m s s , s s , s , s s , s s , s , s p p Q p p q p q p q q ... p q p q q p q p q q TWI λ λ λ λ λ
The inventory carrying costs can now be derived if the time weighted inventory is summed up over all S stages and if the resulting expression is multiplied with the number of lots in the planning period (D/Q) and the unit inventory carrying charges per unit of time h s . Further considering setup costs, transportation costs, and production costs leads to the total cost function:

(12) ( )( ) ( )( ) ( ) ( ) ∑ ∑ ∑ = = = + + + + - + - +         + = S s s s S s s s S s s m s s s m s s s p Dc Q D mT A h p p DQ TC 1 1 1 1 1 1 1 1 1 1 2 λ λ λ λ
Since ( 12) is convex in Q for given values for m and p s , a solution for the optimal order quantity (and with (10) for the optimal size of the first batch) can be calculated with the help of differential calculus. It follows: 13) in ( 12), the total cost function reduces to:

F o r P e e r R e v i e w O n l y (13) ( ) ( )( ) ( )( ) ∑ ∑ = + = - + - +         + + = S s s m s s s m s s s S s s s h p p mT A * Q 1 1 1 1 1 1 1 1 1 2 λ λ λ λ Substituting (
(14) ( )( ) ( )( ) ( ) ( ) ∑ ∑ ∑ = = = + + + - + - +         + = S s s s S s s s S s s m s s s m s s s p Dc mT A h p p D TC 1 1 1 1 1 1 1 1 1 1 2 λ λ λ λ
Since it may be shown that ( 14) is convex in m for given values of p s , an optimal solution may be derived with the help of the following optimality condition:

(15) TC(m*-1) ≥ TC(m*) ≤ TC(m*+1)
However, ( 14) is not necessarily convex in p s . To derive a solution for the production rate, we state the following theorem:

Theorem 4: For m = 1, producing with p s < p s 0 cannot be optimal.

Proof: follows directly from theorem 1.

Further, simulation studies indicated that ( 14) is quasi-convex in p s . To calculate a good, but not necessarily optimal value for p s , we may thus again apply a one-dimensional search algorithm over the intervals [p s,min ; p s,max ] or [p s 0 ; p s,max ], respectively. The following procedure can be used to find a solution for m*, p s *, and Q* (or q 1,s *, respectively):

Step 1: Set m = 1 and search p s ∈ [p s 0 ; p s,max ] that minimises TC. Set TC* = TC, m* = m, and

p s * = p s . Set m = m+1.
Step 2: Search p s ∈ [p s,min ; p s,max ] that minimises TC. If TC > TC*, Goto Step 4.

Step 3: Set TC* = TC, m* = m, p s * = p s , and m = m+1. Goto Step 2.

Step 4: Find Q* from (13) and q 1,s * from (10).

The flexible case a) Equal-sized batch shipments

In the previous section, it was assumed that the production rate has to be fixed before a production run starts. In the following, we relax this assumption and suppose that each batch may be produced with a separate production rate. In compliance with earlier work, we assume that the costs of changing the production rate are negligible. If each batch is manufactured with a separate production rate, the 'regular' inventory given in

(1) has to be reformulated as follows:

(16)

∑ = +         + = m i s , i s , i reg s p p m Q I 1 1 2 2 1 1 2
Similarly, inventory due to waiting times can be calculated as:

(17a) ∑ ∑ ∑ ∑ ∑ ∑ = = - = + = = - = +         - =         - = m i i j s , j i j s , j m i i j s , j i j s , j wait s p p m Q mp Q mp Q m Q I 2 2 1 1 1 2 2 2 2 1 1 1 1 1 for p s > p s+1 and (17b) ∑ ∑ ∑ ∑ ∑ ∑ = - - = + = = - - = + =         - =         - = m i m i j s , j m i j s , j m i m i j s , j m i j s , j wait s p p m Q mp Q mp Q m Q I 2 1 1 1 2 2 2 1 1 1 1 1 for p s < p s+1 and
A general expression can be formulated from ( 17a) and (17b) as follows:

(18)                 -         = ∑ ∑ ∑∑ ∑∑ ∑∑ = - - = + = = = = = - = + m m i j s , j m i j s , j m m i j s , j m i j s , j wait s p , p p , p m Q I 2 i 1 1 1 2 i 2 2 i 2 i 1 1 1 2 2 1 1 Min 1 1 Max
Note that ( 18) and ( 3) are identical for p j,s = p k,s ∀j,k. The time weighted inventory is again given as the sum of ( 16) and ( 18). The inventory carrying costs in the planning period thus equal:

(19) s S s s h m DQ IC ∑ = = 1 2 δ where         -         +         + = ∑ ∑ ∑∑ ∑∑ ∑∑ ∑ = - - = + = = = = = - = + = + m m i j s , j m i j s , j m m i j s , j m i j s , j m i s , i s , i s p , p p , p p p 2 i 1 1 1 2 i 2 2 i 2 i 1 1 1 1 1 1 1 Min 1 1 Max 1 1 2 1 δ
Considering production costs and setup and transportation costs leads to the total cost function: 20) is convex in Q for given values for m and a given p s -vector, the optimal solution for Q is given as follows: 22) is not necessarily convex in p i,s and m, we state the following theorems to derive a good solution:

(20) ( ) ( ) ∑∑ ∑ ∑ = = = = + + + = S s m i s , i s , i S s s s S s s s p c m D Q D mT A h m DQ TC 1 1 1 1 2 δ Since (
Theorem 5: In case p s-1 > p s > p s+1
• it is optimal to increase p 1,s .

• it is optimal to increase p m,s if h s-1 > h s and to reduce p m,s if h s < h s-1 .

• it is optimal to increase p 2,s ,...,p m-1,s if h s-1 > h s and to reduce p 2,s ,...,p m-1,s if h s-1 < h s .

Further, it is optimal to vary p j,s more than p k,s for j < k.

Proof: see appendix 2.

Theorem 6: In case p s-1 > p s < p s+1
• it is optimal to increase p j,s ∀j.

• it is optimal to vary p j,s stronger than p k,s for j < k if h s-1 > h s and vice versa.

Proof: see appendix 3.

Theorem 7: In case p s-1 < p s > p s+1 • it is optimal to increase p 1,s if h s > h s-1 and to reduce p 1,s if h s < h s-1 .
• it is optimal to increase p m,s if h s-1 > h s and to reduce p m,s if h s-1 < h s .

• it is optimal to reduce p 2,s ,...,p m-1,s . Further, it is optimal to reduce p j,s stronger than p k,s for j < k if h s-1 > h s and vice versa.

Proof: see appendix 4.

Theorem 8: In case p s-1 < p s < p s+1 • it is optimal to increase p 1,s if h s > h s-1 and to reduce p 1,s if h s < h s-1 .
• it is optimal to increase p m,s .

• it is optimal to increase p 2,s ,...,p m-1,s in case h s > h s-1 and vice versa. Further, it is optimal to vary p j,s stronger than p k,s for j > k.

Proof: see appendix 5.

To calculate a good, but not necessarily optimal value for p i,s , we can again apply a line search algorithm over the intervals identified above. Further, we calculate a solution for m with the help of the following condition: The following procedure can be applied to find a solution for m*, Q*, and the p s -vector: Step 2: Search p j,s ∈ [p s,min ; p s 0 ] or p j,s ∈ [p s 0 ; p s,max ] that minimises TC according to the conditions specified above. If TC > TC*, Goto Step 4.

Step 3: Set TC* = TC, m* = m, p j,s * = p j,s ∀j,s and m = m+1. Goto Step 2.

Step 4: Find Q* from (21).

b) Unequal-sized batch shipments

In the case where unequal-sized batches are transported to the subsequent stage, the size of the jth batch depends on the size of the first shipment and the production rate of batches 2 to j (cf. figure 3). Thus, it follows that (24)

∏ = + = j i s , i s , i s , s , j p p q q 2 1 1 and ∑∏ ∑ = = + = = = m j j i s , i s , i s , m j s , j p p q q Q 1 2 1 1 1
The time weighted inventory at stage s can be calculated as follows (cf. Glock 2010a):

(25)         +         = =         + + +         + +         + = + = = + + + + ∑ ∏ 1 1 2 2 1 2 1 1 1 2 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 s , j s , j m j j i s , i s , i s , s , m s , m s , m s , m s , m s , s , s , s , s , s , s , s , s , s , s p p p p q p q p q q ... p q p q q p q p q q TWI
The inventory costs may again be calculated by summing up the time weighted inventory over all S stages and by multiplying the resulting expression with the number of lots in the planning period (D/Q) and the unit inventory carrying charges per unit of time h s . Further considering setup, transportation and production costs leads to the total cost function: 26) is convex in Q for given values for m and p i,s , the optimal solution for Q is given as: 27) in ( 26) leads to: 28) is not necessarily convex in p i,s and m, we state the following theorem to derive a good solution:

(26) ( ) ( ) ∑∑ ∑∏ ∏ ∑ ∑ ∑∏ ∑ ∏ = = = = + = + = = = = + + = = + + + +                 +         = S s m j m j j i s , i s , i j i s , i s , i s , j s , j S s s s S s m j j i s , i s , i s s , j s , j m j j i s , i s , i p p p p p c D Q D mT A p p h p p p p DQ TC 1 1 1 2 1 2 1 1 1 2 1 2 1 1 1 2 2 1 1 1 2 Since (
= = + =                 +         + = S s m j j i s , i s , i s s , j s , j m j j i s , i s , i S s s s p p h p p p p mT A * Q 1 2 1 2 1 1 1 2 2 1 1 1 1 2 Substituting (
(28) ( ) ( )                 + +                 +         = ∑∑ ∑∏ ∏ ∑ ∑ ∑∏ ∑ ∏ = = = = + = + = = = = + + = = + S s m j m j j i s , i s , i j i s , i s , i s , j s , j S s S s s s m j j i s , i s , i s s , j s , j m j j i s , i s , i p p p p p c mT A p p h p p p p D TC 1 1 1 2 1 2 1 1 1 2 1 2 1 1 1 2 2 1 1 1 2 Since (
Theorem 9: Producing the first batch with p 1,s < p s 0 cannot be optimal.

Proof: this has been shown in Glock (2010a).

For batches 2 to m, a variation in p i,s may both lead to an increase or to a decrease in inventory. This is due to the fact that the variation of the production rate of a particular batch may both influence the size of subsequent batches as well as the transfer of inventory between subsequent buffer stocks. Thus, complex interdependencies arise between the production rates which are used to manufacture the batches at subsequent stages and the inventory carrying cost charges in the buffer stocks between the stages. Therefore, it is difficult to derive further propositions which are valid in general.

To calculate a good, but not necessarily optimal value for p i,s , we can again apply a line search algorithm over the intervals [p s 0 ; p s,max ] or [p s,min ; p s,max ], respectively. Further, we calculate a solution for m with the help of the following condition: The following procedure can be applied to find a solution for m*, Q*, and the p s -vector:

Step Step 2: Search p j,s ∈ [p s,min ; p s,max ] ∀s and j ∈ [2,…,m] that minimises TC. If TC > TC*, Goto

Step 4.

Step 3: Set TC* = TC, m* = m, p j,s * = p j,s ∀s and j∈[2,…,m], and m = m+1. Goto Step 2.

Step 4: Find Q* from (12).

Numerical Examples

To illustrate the effect of variable production rates on a multi-stage production system with batch shipments, we solved a set of test problems shown in table 1 where we considered three production stages and one final stage representing demand. In addition to the data given below, we assumed D = 1000 and the demand rate for stage S+1 to be 100. Each of these test problems was solved using the solution procedures described in section 3. Tables 2 to 6 summarise the results for the cases of constant, rigid and flexible production rates, respectively.

----------Table 1 ----------As shown in table 2, shipping batches of unequal sizes to the subsequent stage reduces total costs, while shipping batches of equal sizes leads to an increase in total costs. Further, it can be seen that shipping unequal-sized batches reduces inventory in all three buffer stocks (note that I s denotes the buffer stock between stages s and s+1), which is due to the fact that smaller initial shipments enable stages 2,...,S to initiate production earlier, which reduces inventory in the system.

The results shown in tables 3 to 6 illustrate that allowing the production rate to be varied reduces total costs. Further, it becomes obvious that the flexible case leads to lower total costs than the rigid case, and that transporting unequal-sized batches always leads to lower costs than the case where batches of equal size are shipped to the subsequent stage. Although the last result is well known from the literature (see for example [START_REF] Goyal | Determination of Optimum Production Quantity for a Two-Stage Production System[END_REF], it has to be considered that implementing unequal-sized batch shipments may lead to difficulties, for example due to problems in standardising transportation equipment. In such a case, transferring batches of unequal size from one stage to the next may lead to additional costs, which need to be considered when comparing both alternatives, especially when the cost difference between both models is not significant.

---------- ----------If we consider the rigid case first, the total costs of problem one may be reduced by increasing the production rate of the second stage and by reducing the production rate of stages 1 and 3 for both the cases of equal-and unequal-sized batch shipments. A comparison with the results

shown in table 2 indicates that this strategy shifts inventory from the first two buffer stocks to the last one, which is subject to lower inventory carrying charges in the present case. Further, it can be seen that the overall inventory in the system is reduced by moving the production rates closer to each other. Tables 5 and6 show that varying the production rates with every shipment leads to a further reduction in total costs, although inventory is not necessarily reduced. This is due to the fact that it is beneficial to increase the lot size, which may lead to higher inventory carrying costs, but which reduces the number of lots in the planning period and consequently transportation and setup costs.

Problem two illustrates the effect of a variation of the inventory carrying charges on the selection of the production rate. An increase in h s leads to higher inventory carrying costs, which enhances the advantage that results from varying the production rate. Further, as keeping inventory in the third buffer stock is now more expensive than keeping inventory in the first two buffers, it is beneficial to vary the production rates to such an extent that inventory is reduced in the last buffer stock. As can be seen in tables 3 to 6, the production rate of stage 3 is reduced further, which prevents that inventory is transferred too fast from buffer 2 to buffer 3. In addition, the system produces with a smaller lot size, which reduces overall inventory in the system.

If the transportation costs per shipment T s are increased, the number of shipments m is reduced, which, ceteris paribus, leads to higher inventory in the system (cf. problem 3). To compensate this increase, the production rate of the first batches is increased, which reduces inventory. This reaction is obvious, since it has been shown that waiting time-related inventory is especially problematic in case m is high, but that the weight of the production rate of a particular batch is reduced as m decreases (cf. appendices 2 to 5). A similar effect results from varying the setup costs A s , which also impacts the number of batch shipments in the planning period and the necessity of varying p s (cf. problem 4).

----------Table 3 ----------

The impact that results from varying the shape of the unit production cost function is illustrated in problems five and six. If the function is compressed, increasing or decreasing the production rate leads to a sharp increase in production costs. Thus, the higher the slope of the unit production cost function, the lower the variation in p i for both the rigid and the flexible case. Finally, problem seven illustrates the effect of a variation in the interval of feasible production rates. It is obvious that if the interval is extended, the system gains additional flexibility, which may reduce total costs.

----------Table 4 ----------

---------- Table 5 ---------- ---------- Table 6 ----------

Conclusion

This paper extended the work of Glock (2010a) to consider a multi-stage production system.

This system was too investigated for equal and unequal sized shipments between successive stages for production rates that have to be fixed prior to the start of a production run or that may be varied during the production process. The results of the paper showed that • allowing production rates to be varied reduces inventory in the system and results in lower total costs. The case where production rates may be varied for each batch produced leads to lower total costs than the case where the production rates have to be fixed before the start of a production run. Further, unequal-sized batch shipments lead to lower total costs than equal-sized batch shipments.

• system inventory may be reduced by bringing the production rates of subsequent stages of a serial production system closer to each other. This enables the system to process products smoothly without having to build up large buffer stocks in front of or in the back-end of bottleneck stations.

• varying production rates helps to shift inventory between buffer stocks, which enables the production planner to exploit different inventory carrying charges in the system. • the question whether the production rate of a production facility should be increased or reduced depends on (1) the number of batch shipments, (2) the production rates of the preceding and subsequent stage and (3) the inventory carrying charges applicable to the buffer stocks located in front of and in the back-end of the respective production stage.

Our research is especially important for companies facing high in-process inventories or bottleneck-problems. As could be shown, varying the production rates in a serial production system may significantly influence waiting time-related inventory and may thus help to control queues in front of subsequent production stages. This is consistent with the OPT-Philosophy (see for example [START_REF] Goldratt | Computerized shop floor scheduling[END_REF], which aims on maximizing the throughput of the production system as a whole, instead of the throughput of a single machine, and the JITphilosophy, which aims on reducing in-process inventories.

One important aspect that has not been addressed in this paper is that varying the production rate may impact the quality of the product under consideration. Thus, with an increasing or a decreasing production rate, the production process may produce a higher or a lower percentage of defective products that need to be reworked or scrapped. Although the production cost function formulated in (5) may reflect this relation in part, more complex interdependencies between the production rate of a production facility and the quality of the products produced on this facility may prevail in practice. Consequently, it might be necessary to consider a separate quality cost function to account for scrap, rework or shortages as a result of a variation in the production rates. One approach could be to use an exponentially distributed random variable, which describes the point in time at which a production process goes out of control, and to establish a connection between the production rate of the production facility and the mean of the variable (see for example [START_REF] Rosenblatt | Simultaneous determination of production cycle and inspection schedules in a production system[END_REF][START_REF] Zhu | Pushing Quality Improvement Along the Supply Chain[END_REF]). Thus, with an increase or a decrease of the production rate, the mean of the random variable would be changed, resulting in a lower or higher probability of the production process going out of control.

In order to increase the scope of our analysis, the model presented in this article could further be extended to include costs of changing the production rate or a different unit production cost function (see for example [START_REF] Khouja | The economic production lot size model under volume flexibility[END_REF][START_REF] Moon | Determination of optimal production rates on a single facility with dependent mold lifespan[END_REF]. Further, it would be interesting to permit a different number of batch shipments at each production stage, which has been shown to reduce total costs [START_REF] Bogaschewsky | Optimizing multi-stage production with constant lot size and varying number of unequal sized batches[END_REF]. If we consider three subsequent stages, we may differentiate between the four cases illustrated in table 7.

Appendices

---------- Table 7 ----------While component II is always positive for cases 1 to 4, component I adopts a positive value for m = 1 and a negative value for m > 2. It is obvious that in case of m = 1, all four cost functions may be reduced by increasing p s , while a reduction in p s leads to an increase in inventory carrying costs. Due to the fact that a variation in p s always leads to an increase in production costs, it cannot be optimal to produce with a production rate lower than p s,0 in case

m = 1.
In contrast, if m = 2, component I reduces to zero in all four cost functions (cf. table 7). Thus, inventory carrying costs may be reduced by reducing component II. As can be seen, component II is independent of p s in case 2, wherefore we can conclude that if p s-1 < p s > p s+1 , it is not beneficial to produce with another production rate than p s,0 , which would result in increased production costs, but not in lower inventory carrying costs. This effect is due to the fact that in case p s-1 < p s > p s+1 and m = 2, a variation in p s leads to a decrease (increase) in I reg that is exactly compensated by an increase (decrease) in I wait (it can be seen in figure 1 by comparing the bold and dashed lines that a decrease in p s leads to a higher regular inventory and a lower inventory due to waiting times). This has been shown by Glock (2010a) in the context of a two-stage production system.

For the other three cases, i.e. p s-1 > p s > p s+1 , p s-1 > p s < p s+1 , and p s-1 < p s < p s+1 , inventory carrying costs may be reduced by increasing p s . This effect can be explained by analysing the variation of I reg and I wait for the case that p s is varied to p s '. It can be seen that (A-3) is positive for m = 2 and p s ' > p s , wherefore we may conclude that an increase in p s leads to a reduction in the regular inventory that compensates the increase in waiting time-related inventory. Thus, since c p,s (p s ) increases with an increasing or a decreasing production rate, it is optimal to produce with p s > p s,0 . Finally, for m > 2, it can be seen that in case p s-1 > p s > p s+1 , the second bracket term in (A-3) is positive for m/(m-2) > h s /h s-1 . In this case, it is beneficial to increase p s , which would lead to a positive value for the first bracket term and consequently lower inventory carrying costs.

(A-3) ( ) ( ) ( )( ) ( ) ( ) ( )              < < + -         - < > +         - > < - +         - > > + -         - = ∆ + ∆ + - - + - - + - - + - - 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 2 2 1 1 2 2 1
In contrast, if m/(m-2) < h s /h s-1 , it is beneficial to reduce p s .

For the case p s-1 < p s > p s+1 , it can be seen that the third bracket term in (A-3) is always negative for m > 2. Thus, inventory carrying costs can be reduced by reducing p s , which leads to a positive value for ∆I reg + ∆I wait in total. If p s-1 > p s < p s+1 holds, (A-3) illustrates that the sum of ∆I reg + ∆I wait is independent of m, and that the change in inventory which results from varying p s to p s ' may be increased by enhancing the production rate at stage s.

Finally, in case p s-1 < p s < p s+1 , it can be seen that the second bracket term in (A-3) is positive for m/(m-2) > h s-1 /h s . In this case, it is beneficial to increase p s , which would lead to a positive value for the first bracket term and consequently lower inventory carrying costs. To the contrary, if m/(m-2) < h s /h s-1 , it is beneficial to reduce p s .

Appendix 2

In case p s-1 > p s > p s+1 , inventory costs for three subsequent stages as given in (19) reduce to:

(A-4) ( ) ( ) ( )     + - + - - + + - +     + +     + - + - + - -         + - - - + - - - - = + - = - - = - - - - - = + ∑ ∑ ∑ ∑ 1 III 1 1 1 1 1 II 1 1 1 1 1 1 2 1 I 1 2 1 1 2 1 1 1 1 1 2 1 2 1 2 1 2 2 1 2 3 2 2 3 s , m s s , m s s s , m s s , s s , s s s , s i j s , j s i j s , i s s s , j s s i j s , j s s , i s m i s , i s p h p h h p h p m h p h h m p h p h p h h p h h p h p h p h m DQ 4 3 42 1 4 4 3 4 4 2 1 4 4 4 4 8 4 4 4 4 7 6
As can be seen, only component II depends on p 1,s , i.e. the production rate of the first batch produced on stage s. Since component II is always positive, inventory carrying costs may be reduced by increasing p 1,s . This result is obvious, as an increase in the production rate of the first batch leads to a faster inventory build-up and enables stage s+1 to initiate consumption earlier, which reduces inventory in the system. Looking at the production rate of the last batch, it can be seen in (A-4) that only component III depends on p m,s . Due to the fact that III is positive for h s-1 > h s and negative for the opposite case, inventory carrying costs may be reduced by increasing p m,s if h s-1 > h s and vice versa. This result may be explained by the fact that a variation in p m,s leads to a shift of inventory between the buffer stocks in front of stages s and s+1, which enables the system to take advantage of different inventory carrying charges along the supply chain.

Finally, considering component I, it can be seen that an increase in p j,s for j = 2,...,m-1 leads to a reduction in inventory if h s-1 > h s and vice versa for the same reasons as explained above.

Further, it becomes obvious that the sum function in the component leads to a higher weight for the production rate of earlier batches, which is due to the fact that the production rate of batch j with j > 1 influences the waiting times of all subsequent batch shipments. Thus, if p j,s is reduced, the reduction for p j,s is larger than for p k,s with j < k, and vice versa.

Appendix 3

In case p s-1 > p s < p s+1 , inventory costs for three subsequent stages as given in ( 19) reduce to: In this case, component II and III, which are always positive, depend on p 1,s and p m,s , respectively. Thus, inventory carrying costs may be reduced by increasing p 1,s and p m,s , which leads to a lower value for components II and III. Further, it can be seen that component I is strictly positive and that it may be reduced by increasing p j,s for j = 2,...,m-1. Thereby, the sum function in component I leads to a higher weight for the production rate of earlier batches if h s-1 > h s and to a higher weight for the production rate of later batches if h s-1 < h s . This may be explained by the fact that a faster production rate for later batches enables the system to initiate consumption earlier, which reduces inventory in front of stage s+1. In contrast, if earlier batches are produced with a higher production rate, waiting times are reduced for batches stored in the buffer stock in front of stage s. Thus, by varying the production rates of subsequent batches differently, the system may take advantage of different inventory carrying costs in the system. As a consequence, we may conclude that it is optimal to increase p j,s stronger than p k,s for j < k if h s-1 > h s and vice versa. In case p s-1 < p s > p s+1 , inventory costs for three subsequent stages as given in (19) reduce to: As can be seen in (A-6), component II is positive in case h s > h s-1 and negative for the opposite case. Thus, if h s > h s-1 , it is optimal to increase p 1,s , while it is optimal to reduce p 1,s if h s < h s-1 . The opposite pertains for component III, which is positive for h s-1 > h s and negative for h s-1 < h s . Thus, we conclude that it is optimal to increase p m,s if h s-1 > h s and to reduce p m,s if h s-1 < h s .

(A-5) ( ) ( )     + - + + - - + - +     + +     - - + + + +         - - + - - - + - - - - = + + = - = - - - = = - - - - ∑ ∑ ∑ ∑ ∑ 1 III 1 1 1 1 1 II 1 1 1 1 1 1 1 1 I 2 1 1 1 1 2 2 1 1 1 1 2 1 2 1 2 2 1 2 2 2 s ,
(A-6) ( ) ( )     + - + - + - + - +     + +     + + - + -         + + + - - - + - - - - = + - = = - - - = + - - - = - - ∑ ∑ ∑ ∑ ∑ 1 III 1 1 1 1 1 II 1 1 1 1 1 1 2 1 I 1 2 1 1 1 2 1 1 1 1 1 1 2 1 2 1 2 2 1 2 2 3 2 s ,
As to the production rate for batches 2 to m-1, it can be seen that component I is strictly negative, wherefore inventory carrying costs may be reduced by reducing p j,s for j = 2,...,m-1.

As explained above, the sum function in component I leads to a higher weight for the production rate of earlier batches if h s-1 > h s and to a higher weight for the production rate of later batches if h s-1 < h s . As a consequence, we may conclude that it is optimal to reduce p j,s stronger than p k,s for j < k if h s-1 > h s and vice versa. Component II of (A-7) is positive in case h s > h s-1 and negative for the opposite case. Thus, if h s > h s-1 , it is optimal to increase p 1,s , while it is optimal to reduce p 1,s if h s < h s-1 . In contrast, component III is always positive, wherefore it is inventory carrying costs may be reduced by increasing p m,s .

Finally, component I is positive for h s > h s-1 and negative for the opposite case. Thus, in case h s > h s-1 , it is optimal to increase p j,s and to reduce p j,s if h s < h s-1 . Again, the sum function in component I leads to a higher weight for the production rate of later batches, wherefore it is beneficial to vary p j,s stronger than p k,s for j > k. .34, 233.36, 237.22, 240.84, 244.22, 247.42, 250.47} {250, 247.82, 244.43, 240.89, 237.16, 233.24, 250} {298.80, 270, 270, 270, 270, 270, 301.19} 181 .38, 247.29, 248.07, 248.85, 249.62} {206.49, 206.49, 206.49, 206.49, 206.49} {299.87, 296.91, 296.63, 296.35, 300 .04, 249.74, 249.81, 249.89, 249.96} {200.68, 200.68, 200.68, 200.68, 200.68} {299.99, 299.70, 299.68, 299.65, 300 .21, 233.10, 236.32, 240.16, 244.60, 249.01, 251.99} {251.99, 248.40, 244.30, 242.37, 234.96, 233.41, 255.07} {297.52, 250, 250, 250, 250, 250, 299 .05, 250.52, 249.31, 249.12, 249.91, 251.84} {212.22, 221.82, 228.77, 232.45, 233.27, 231.97} {320, 270, 270, 270, 270, 287 .30, 250.80, 250.53, 251.45} {212.41, 219.83, 224.13, 223.98} {320, 270, 270, 271.14} 170.53 166.08 515.85 15019.4 3 1 2 3 5 357.87 {254.13, 251.65, 250.75, 251.15, 252.92} {215.14, 224.85, 231.38, 233.92, 233.62} {320, 270, 270, 270, 286 .22, 249.59, 249.26, 249.21, 249.43, 249.91} {201.43, 203.39, 205.29, 207.11, 208.78, 209.80} {320, 273.12, 282.57, 293, 296.68, 298 .03, 249.96, 249.93, 249.94, 249.99} {200.16, 200.42, 200.66, 200.89, 201.09} {303.72, 298.19, 299.19, 299.65, 299 , 250.08, 248.68, 248.54, 249.53, 251.71} {211.85, 220.15, 225.64, 229.21, 231.26, 231} {350, 271.87, 250, 250, 264.69, 287

  finally developed a model where subsequent batch shipments are of unequal size. The complexity of the model required developing algorithms for finding the optimal solution.

  p s ) = unit production cost function of stage s A s = setup costs per setup at stage s T s = transportation costs per shipment at stage s d = demand rate in units per unit time D = total demand in the planning period h s = inventory carrying charges per unit per unit of time at stage s s = size of the jth batch produced at stage s λ s = the proportional increase/decrease in successive shipments within a batch production run at stage s m = number of batch shipments in a production cycle p s = production rate in units per unit of time at stage s with p S+1 = d p s 0 = design production rate at stage s a s,i = cost parameters of the unit production cost function Q = production lot size S = number of stages Definitions: I s reg = 'regular' inventory at stage s I s wait = inventory due to waiting times at stage s IC = inventory carrying costs PC = production costs TC = total costs in the planning period TWI = time weighted inventory |a| = the absolute value of a, i.e.

  at stage s+1 before the first batch is completely used up at stage s+2,

( 5 )

 5 PC s = Dc s (p s ) = D(a s,0 p s 2 -a s,1 p s +a s,2 )

  Figure2 ----------

  1: Set m = 1 and search p s ∈ [p s 0 ; p s,max ] ∀s that minimises TC. Set TC* = TC, m* = m, and p s * = p s ∀s. Step 2: Set m = 2 and p s = p s 0 ∀s where p s-1 < p s > p s-1 . Otherwise, search p s ∈ [p s 0 ; p s,max ] that minimises TC. If TC > TC*, Goto Step 5. Set TC* = TC, m* = m, p s * = p s ∀s, and m = m+1

( 23 )

 23 TC(m*-1) ≥ TC(m*) ≤ TC(m*+1)

  Set m = 1 and search p 1,s ∈ [p s 0 ; p s,max ] or p 1,s ∈ [p s,min ; p s 0 ], respectively, that minimises TC. Set TC* = TC, m* = m, and p 1,s * = p 1,s ∀s. Set m = m+1.

( 21 )

 21 TC(m*-1) ≥ TC(m*) ≤ TC(m*+1)

  1: Set m = 1 and search p 1,s ∈ [p s 0 ; p s,max ] that minimises TC. Set TC* = TC, m* = m, and p 1,s * = p 1,s . Set m = m+1.

Figure 1 :

 1 Figure 1: Inventory-time plots for the case of equal-sized batch shipments Figure 2: A three-stage production system with two buffer stocks Figure 3: Inventory-time plots for the case of unequal-sized batch shipments

Table 2

 2

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

Table 1 # s h s A s T s {a 0,s ,a 1,s ,a 2,s } {p s,min , p s,max }

 1

		1	3	250	20	{1/6000, 1/12, 11}	{230, 300}
	1	2	3	275	30	{1/7000, 2/35, 6.5}	{170, 250}
		3	2	200	25	{1/9000, 1/15, 10.5}	{270, 320}
		1	3	250	20	{1/6000, 1/12, 11}	{230, 300}
	2	2	5	275	30	{1/7000, 2/35, 6.5}	{170, 250}
		3	10	200	25	{1/9000, 1/15, 10.5}	{270, 320}
		1	3	250	40	{1/6000, 1/12, 11}	{230, 300}
	3	2	3	275	60	{1/7000, 2/35, 6.5}	{170, 250}
		3	2	200	50	{1/9000, 1/15, 10.5}	{270, 320}
		1	3	125	20	{1/6000, 1/12, 11}	{230, 300}
	4	2	3	137.5	30	{1/7000, 2/35, 6.5}	{170, 250}
		3	2	100	25	{1/9000, 1/15, 10.5}	{270, 320}
		1	3	250	20	{1/600, 10/12, 105.17}	{230, 300}
	5	2	3	275	30	{1/700, 20/35, 58.06}	{170, 250}
		3	2	200	25	{1/900, 10/15, 100.54}	{270, 320}
		1	3	250	20	{1/60, 100/12, 1042.67}	{230, 300}
	6	2	3	275	30	{1/70, 200/35, 572.34}	{170, 250}
		3	2	200	25	{1/90, 100/15, 1000.54}	{270, 320}
		1	3	250	20	{1/6000, 1/12, 11}	{200, 320}
	7	2	3	275	30	{1/7000, 2/35, 6.5}	{150, 270}
		3	2	200	25	{1/9000, 1/15, 10.5}	{250, 350}

Table 3

 3

	# m*	Q*	{p 1 ,...,p S }*	I 1	I 2	I 3	TC*
	1 7	332.88 {244.30, 244.30, 270}	194.66 240.98 1224.09 9764.98
	2 4	143.19 {244.62, 235.99, 270}	157.03 170.79 583.35 16476.0
	3 5	337.01 {245.84, 248.84, 270}	274.17 310.98 1310.61 11025.5
	4 4	209.50 {241.95, 238.63, 285.58} 222.49 255.55 864.08 8440.81
	5 5	263.2 {248.85, 206.48, 298}	320.05 372.38 1051.03 10885.8
	6 5	258.42 {249.88, 200.68, 299.81} 334.93 386.78 1037.53 10942.3
	7 7	336.81 {244.48, 244.48, 259.96} 196.81 226.10 1221.3 9757.59

Table 4

 4

	# m*	Q*	{p 1 ,...,p S }*	I 1	I 2	I 3	TC*
	1 6	318.542 {249.71, 228.04, 275.30} 228.04 234.32 1018.85 9157.69
	2 4	157.43 {250.96, 220.61, 270}	171.08 170.24 514.63 15051.3
	3 4	319.32 {252.71, 228.28, 284.91} 337.09 333.94 1068.16 10308.6
	4 4	224.06 {251.48, 222.43, 288.70} 241.74 241.43 753.64 7869.0
	5 5	294.24 {249.73, 204.40, 297.75} 282.35 306.84 985.46 9970.88
	6 5	291.83 {249.97, 200.47, 299.77} 287.24 315.56 980.46 9997.18
	7 6	318.52 {249.71, 228.04, 275.30} 228.04 234.32 1018.85 9157.69

Table 5

 5

	# s m* Q*	{p 1,s ,...,p m,s }*	I s	TC*
		1	{253
	1	2	7 333.67
		3

Table 6

 6

	# s m*	Q*	{p 1,s ,...,p m,s }*	I s	TC*
		1		{253
	1	2	6 331.13
		3

Table 1 :

 1 Test problems used for computational experimentation

Table 2 :

 2 Computational results for the case of fixed production rates

Table 3 :

 3 Computational results for the rigid case with equal-sized batch shipments

Table 4 :

 4 Computational results for the rigid case with unequal-sized batch shipments

Table 5 :

 5 Computational results for the flexible case with equal-sized batch shipments

Table 6 :

 6 Computational results for the flexible case with unequal-sized batch shipments

Table 7 :

 7 Inventory carrying cost for alternative ratios of the production rates

Page 32 of 35 http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.uk International Journal of Production Research

	Page 33 of 35	t w,s,2		t w,s,1
	Q
	m
					Time
	Inventory at stage s+1		t p,s	t w,s,2 = waiting time of batch 2 t w,s,1 = waiting time of batch 1
	Q	F o
	m	r
			P
			e r e	Time
	Inventory at stage s+1 Q	t w,s+1,1	R e v t w,s+1,2 i w e	t w,s+1,1 = waiting time of batch 2 t w,s+1,2 = waiting time of batch 3
	m			n O
		t p,s+1	t p,s+2		l y Time
	Inventory at stage s+2
	Q
	m
					Time

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

uk International Journal of Production Research

					Page 34 of 35
	stage s-1	buffer	stage s	buffer	stage s+1
	p s-1		p s		p s+1

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

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

uk International Journal of Production Research

	Page 35 of 35
				Time
	Inventory at stage s+1
	F o
	r
	P
	e r e		Time
	Inventory at stage s+1	w e i v e R
	q 1,s+1	O n
	t p,1,s+1 t p,2,s+1	t p,3,s+1	l y	Time
	Inventory at stage s+2
				Time

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

uk International Journal of Production Research

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

Acknowledgement

The author is grateful to the anonymous referees for their valuable comments and suggestions on an earlier version of this paper. The author further wishes to thank the editorial team and the reviewers for the fast reviewing process.

Appendix 5

In case p s-1 < p s < p s+1 , inventory costs for three subsequent stages as given in (19) reduce to: