Pedro Terán
email: teran@unizar.es

On a uniform law of large numbers for random sets and subdifferentials of random functions

Keywords: Clarke generalized gradient, Uniform Strong Law of Large Numbers, Upper semicontinuous multifunction, Subdifferential

In this paper, we strengthen the convergence in a uniform law of large numbers for random upper semicontinuous multifunctions of Shapiro and Xu. The proof is based on an abstract law of large numbers in a metric space endowed with a convex combination operation. Convergence in the Hausdorff metric is obtained, whereas the original result presented a weakened form of convergence of excess functionals. As a consequence, another law of large numbers for subdifferentials of random functions is improved as well.

Introduction

Shapiro and Xu [START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF] have recently presented a uniform law of large numbers for multifunctions, motivated by consistency analysis of stationary points but the non-uniqueness of such points requires replacing the ordinary law of large numbers for random variables by a law of large numbers for random sets (or multifunctions) given by subdifferentials of the random functions above.

The most salient feature of their result is that the LLN is needed to be uniform over compact sets but the subdifferential (Clarke generalized gradient) at a point does not vary continuously, it is only upper semicontinuous in the usual sense of multifunctions. This leads Shapiro and Xu to working with a one-sided version of the Hausdorff metric (convergence of excess functions). In [START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF]Theorem 2], they essentially prove that, if X is a random multifunction, upper semicontinuous, then almost surely the average n -1 X i (x, •) is asymptotically covered by the unions of expectations of X(y, •) where y is near x, and the expectation of X(x, •) is asymptotically covered by the unions of averages at points near x, uniformly as x ranges over a compact metric space.

Thus two (asymmetric) convergences are proven, each being a relaxed form of half what it takes to prove Hausdorff convergence. This is to account for X being only u.s.c. Of course, if X is assumed to be continuous then one can do without those relaxations and convergence in the Hausdorff metric is reached [START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF]Theorem 4].

They remark [6, Remark 2] that continuity might be needed to obtain convergence in the Hausdorff metric. That is eminently reasonable, yet we set out to prove that it is not correct. In fact, the uniform law of large numbers can hold without relaxing Hausdorff convergence if X is only upper semicontinuous, and we will show that even this is not strictly necessary. This leads to an important sharpening of their basic result.

The uniform law of large numbers will be obtained as an application of an abstract law of large numbers in metric spaces due to Molchanov and the present author [START_REF] Terán | The law of large numbers in a metric space with a convex combination operation[END_REF]. That paper presents a set of five axiomatic conditions under which an otherwise arbitrary convex combination operation admits

A c c e p t e d m a n u s c r i p t

the construction of an expectation operator satisfying the strong law of large numbers.

In this paper, the carrier space will be that of upper semicontinuous multifunctions on a compact space. Proving the law of large numbers reduces to showing that the natural definition of convex combinations fulfils those five conditions.

Convex combination spaces

Let (E, d) be a metric space endowed with a convex combination operation [•, •] yielding a point of E for each tuple of weights (λ 1 , . . . , λ n ) and points (u 1 , . . . , u n ):

(λ 1 , . . . , λ n , u 1 , . . . , u n ) → [λ 1 , u 1 ; . . . ; λ n , u n ],
where n ≥ 2, λ i > 0,

λ i = 1, u i ∈ E.
(i) (Commutativity) For every permutation σ of {1, . . . , n}, [λ i , u i ] n i=1 = [λ σ(i) , u σ(i) ] n i=1 ; (ii) (Associativity) [λ i , u i ] n+2 i=1 = [λ 1 , u i ; . . . ; λ n , u n ; λ n+1 +λ n+2 , [ λ n+j λ n+1 +λ n+2 ; u n+j ] 2 j=1 ]; (iii) (Continuity) If u, v ∈ E and λ (k) → λ ∈ (0, 1) as k → ∞, then [λ (k) , u; 1 -λ (k) , v] → [λ, u; 1 -λ, v] ; (iv) (Negative curvature) For all u 1 , u 2 , v 1 , v 2 ∈ E and λ ∈ (0, 1), d([λ, u 1 ; 1 -λ, u 2 ], [λ, v 1 ; 1 -λ, v 2 ]) ≤ λd(u 1 , v 1 ) + (1 -λ)d(u 2 , v 2 ) ;
(v) (Convexification) For each u ∈ E, there exists lim n→∞ [n -1 , u] n i=1 , which will be denoted by Ku (or K E u if there is ambiguity).

A c c e p t e d m a n u s c r i p t

These axioms are substantially weaker than those defining a normed linear space. It is proven in [START_REF] Terán | The law of large numbers in a metric space with a convex combination operation[END_REF]Sections 4 and 5] that they suffice to construct an expectation operator E satisfying the law of large numbers. E is defined on the space L 1 E of all Borel functions X such that Ed(X, u) is finite for some arbitrary u ∈ E. Examples, besides the usual operations in a Banach space, include their uplift to compact sets, convolution of measures, some geometric definitions of average in metric spaces, and others (see [START_REF] Terán | The law of large numbers in a metric space with a convex combination operation[END_REF]Section 9]). The aim was to find a unifying framework valid for various non-Banach spaces.

Theorem 1. Let E be a separable complete convex combination space. Then E admits an expectation operator E such that

(a) If X = r j=1 I Ω j u j is a simple function, then EX = [P (Ω j ), Ku j ] r j=1 ; (b) E extends by continuity to all X ∈ L 1 E ; (c) If X ∈ L 1 E and {X i } i∈N are pairwise i.i.d. as X, then [n -1 , X i ] n i=1 → EX almost surely.

Preliminaries

The space R d is endowed with the Euclidean norm. Its closed unit ball will be denoted by B, and its open unit ball by B o .

For any closed set A ⊂ R d , we denote by K(A) the family of all nonempty compact subsets of A. The excess of A over a set C is defined to be e(A, C)

= inf{ε > 0 | A ⊂ C + εB}.
The Hausdorff metric between them is

d H (A, C) = max{e(A, C), e(C, A)}.
Thus the excess is just a one-sided version of the Hausdorff metric. The norm of A is defined to be

A = d H (A, {0}) = sup x∈A |x|.
Note that the space K(R d ) is an example of convex combination space, having the convex hull mapping co as its convexification operator. In this case, the expectation operator is the convex hull of the Aumann expectation, see [START_REF] Artstein | A strong law of large numbers for random compact sets[END_REF].

A c c e p t e d m a n u s c r i p t

An important property of d H is

d H (co A, co C) ≤ d H (A, C).
In fact, the convexification operator of any convex combination space is always non-expansive [7, Proposition 3.6].

A random closed set is an Effros measurable mapping X from a measurable space (Ω, A) to F , i.e. the sets

{ω ∈ Ω | X(ω) ∩ U = Ø} are measurable for each open U ⊂ E. Functions from (Ω, A) to (K(R d ), d H )
are random closed sets if and only if they are Borel.

The abbreviation p.i.i.d. stands for 'pairwise independent and identically distributed'.

The space of u.s.c. multifunctions

Let (X , ρ) denote a compact metric space, and let U be the space of upper semicontinuous multifunctions on X with non-empty compact values, namely u ∈ U if for every ε > 0 and x ∈ R d there is some δ > 0 such that ρ(x, y) < δ implies u(y) ⊂ u(x) + εB o . Equivalently, the excess of u(y) over u(x), e(u(y), u(x)) = inf{ε > 0 | u(y) ⊂ u(x) + εB} goes to 0 as y → x.

We define naturally the convex combination of several u.s.c. multifunctions u i to take on the values

[λ i , u i ] n i=1 (x) = n i=1 λ i u i (x).
The space U is endowed with the uniform metric

H(u, v) = sup x∈X d H (u(x), v(x)).
Given an underlying measurable space (Ω, A), a random u.s.c. multifunction is defined in [START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF] to mean a mapping X : Ω × X → K(R d ) such that X(•, x) is measurable for each x ∈ X . It can be equivalently regarded as a random

A c c e p t e d m a n u s c r i p t

element of U with its cylindrical σ-algebra, and we will not really distinguish between both descriptions. For instance, in this context a 'Borel function' will be a random element of U with the Borel σ-algebra generated by the metric H. Notice that, even if X is compact metric and K(R d ) is separable, (U, H) is not separable in general. In order to apply Theorem 1, we will need to impose the assumption that the range (modulo a null set) of the random multifunction is separable. This is equivalent to demanding that the random multifunctions can be approximated by simple U-valued mappings; under the Continuum Hypothesis or even weaker set-theoretical assumptions, every Borel U-valued function has that property 1 . However, this is not a restriction to our main aim of showing that LLN convergence in the Hausdorff metric can be achieved, with a certain generality, for non-continuous multifunctions.

Theorem 2. The space (U, H) is a convex combination space where the convexification operator K U is given by

(K U u)(x) = co u(x). Proof. Set R(u) = x∈X u(x).
We begin by proving that R(u) is compact. (a) R(u) is bounded. By the upper semicontinuity of u, for each x ∈ X there exists some δ(x) > 0 such that

ρ(x, y) < δ(x) =⇒ u(y) ⊂ u(x) + B o . Since {x + δ(x)B o } x∈X is an open cover of X , it has a finite subcover {x i + δ(x i )B o } k i=1 . Consequently, for each y ∈ X , u(y) ⊂ k i=1 [u(x i ) + εB o ],
1 Being a family of functions from X to K(R d ), the space (U, H) has separability character bounded above by the cardinality of P(X ×K(R d )). Since both factors are separable metric spaces, that is in turn bounded by 2 c . Under the hypothesis that c is a measure-free cardinal (much weaker than the Continuum Hypothesis), 2 c is so also and [4, Proposition 438D] or [3, Appendix! C, Theorem 3.2] ensure that every U-valued Borel function has essentially separable range.

A c c e p t e d m a n u s c r i p t

so that R(u) is included in a finite union of bounded sets. (b) R(u) is closed. Take z n → z with {z n } n ⊂ R(u). Then there exist x n ∈ X such that z n ∈ u(x n ) for each n. Because X is compact metric, {x n } n has a subsequence (still denoted x n ) converging to some x ∈ X . Since u is upper semicontinuous, for any given ε > 0 and all sufficiently large n we have

z n ∈ u(x n ) ⊂ u(x) + εB.
But the latter set is closed, so z ∈ u(x) + εB. The arbitrariness of ε yields z ∈ u(x) ⊂ R(u).

We proceed now to checking properties (i) through (v). Properties (i) and (ii) are immediate. The compactness of R(u) for u ∈ U is the key tool in order to prove properties (iii) and (v).

As to (iii), we have to show that

sup x∈X d H (λ (k) u(x) + (1 -λ (k) )v(x), λu(x) + (1 -λ)v(x)) → 0 whenever λ (k) → λ ∈ (0, 1). To that purpose, it is enough to prove sup A∈K(R(u)),C∈K(R(v)) d H (λ (k) A + (1 -λ (k) )C, λA + (1 -λ)C) → 0.
Reasoning by contradiction, assume the contrary. Then there exists a subsequence {k } k , some ε > 0 and some sets

A k ∈ K(R(u)), C k ∈ K(R(v)) such that ε < d H (λ (k ) A k + (1 -λ (k ) )C k , λA k + (1 -λ)C k ) for all k.
By the compactness of R(u) and R(v), the family

K(R(u)) × K(R(v)) is d H -compact [5, Propositions 1-4-1 and 1-4-4]. Then we take a further subsequence {k } k such that A k → A, C k → C for some A ∈ K(R(u)), C ∈ K(R(v)). By the triangle inequality, ε < d H (λ (k ) A k + (1 -λ (k ) )C k , λ (k ) A + (1 -λ (k ) )C)+ d H (λ (k ) A + (1 -λ (k ) )C, λA + (1 -λ)C)+ d H (λA + (1 -λ)C, λA k + (1 -λ)C k ) = (I) + (II) + (III). But (I) ≤ λ (k ) d H (A k , A) + (1 -λ (k ) )d H (C k , C) → 0,

A c c e p t e d m a n u s c r i p t

similarly (III) → 0 and finally

(II) ≤ d H (λ (k ) A, λA) + d H ((1 -λ (k ) )C, (1 -λ)C) ≤ |λ (k ) -λ| • A + |λ (k ) -λ| • C → 0.
This is a contradiction, so (iii) is proven. Property (iv), negative curvature of U, follows easily from that of K(R d ):

H(λu 1 + (1 -λ)u 2 , λv 1 + (1 -λ)v 2 ) = sup x∈X d H (λu 1 (x) + (1 -λ)u 2 (x), λv 1 (x) + (1 -λ)v 2 (x)) ≤ sup x∈X (λd H (u 1 (x), v 1 (x)) + (1 -λ)d H (u 2 (x), v 2 (x))) ≤ ≤ λH(u 1 , v 1 ) + (1 -λ)H(u 2 , v 2 ).
Finally, in order to establish property (v), the existence of K U , we need to prove

H(n -1 n i=1 u, K U u) → 0.
Similar to (iii), it suffices to show that sup

A∈K(R(u)) d H (n -1 n i=1 A, co A) → 0.
Assume the contrary; a similar compactness argument yields a value ε > 0, a subsequence {n } n and sets

A n , A ∈ K(R(u)) such that d H (A n , A) → 0 and ε < d H ((n ) -1 n i=1 A n , co A n ) ≤ d H ((n ) -1 n i=1 A n , (n ) -1 n i=1 A) + d H ((n ) -1 n i=1 A, co A) +d H (co A, co A n ) = (IV ) + (V ) + (V I).

A c c e p t e d m a n u s c r i p t

Then (IV ) ≤ n • (n ) -1 d H (A n , A) → 0,
(V ) goes to 0 by [START_REF] Artstein | Convexification in limit laws of random sets in Banach spaces[END_REF] and

(V I) ≤ d H (A, A n ) → 0.
A contradiction obtains, so (v) is proven.

We will also need the following lemma about U.

Lemma 3. The metric space (U, H) is complete.

Proof. Since K(R d ) X , with the uniform metric, is complete, it suffices to prove that U is closed in that space. The proof is rather standard. Let

{u n } n ⊂ U H-converge to some u ∈ K(R d ) X
, and fix ε > 0. Then, d H (u n (y), u(y)) < ε/3 for some n ∈ N and all y ∈ X . Let x ∈ X . By the upper semicontinuity of u n there exists δ > 0 such that ρ(x, y) < ε implies

u n (x) ⊂ u n (y) + (ε/3)B o . We deduce that u(x) ⊂ u n (x) + ε 3 B o ⊂ u n (y) + 2ε 3 B o ⊂ u(y) + εB o .

Laws of large numbers

It is not hard now to derive the LLN for random u.s.c. multifunctions. We adopt the following notation from [START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF]: if X n are p.i.i.d. as X, Proof. For each n ∈ N there exists a null set N n ⊂ Ω such that the image of Ω\N n by X n is H-separable. The H-closure of n X n (Ω\N n ), denoted V, is then a complete (by Lemma 3) separable convex combination space. Theorems 2 and 1 yield now the following. Let E be the expectation operator in the space V. Then H(n -1 n i=1 X i , EX) → 0 almost surely, under the assumption that X ∈ L 1 V , namely EH(X, u) < ∞ for some (arbitrary) u ∈ U.

S n (x) = n -1 n i=1 X i (x, •), E(x) = E co X(x, •).
But taking u to be the null function, we obtain

H(X(ω), 0) = sup x∈X X(x, ω) ≤ κ(ω)
so the existence of κ implies that X ∈ L 1 V indeed holds. Moreover, one checks routinely, using approximations by simple functions, that E(x) = (EX)(x) for every x ∈ X , so convergence

H(n -1 n i=1 X i , EX) → 0 is the same thing as sup x∈X d H (S n (x), E(x)) → 0.
This should be compared to [START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF]Theorem 2], where it is proven (without the separability restriction on the range of X) that, for any r > 0, In this case the excess is the same as the Hausdorff metric, but observe that every u.s.c. singleton-valued function u identifies with the continuous function ũ given by u(x) = {ũ(x)}. But we have just shown the same conclusion by requiring only compact range: now x u(x) is just the range of ũ.

The application of this result to Clarke generalized gradients of random functions mirrors that in [START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF]. Since it follows by the same techniques therein, proofs are omitted. The reader is referred to [6, Section 3], on which the following is heavily based, for details.

Let {ξ n : Ω → R d } n be a p.i.i.d. sequence supported on a closed set Ξ. Consider a Carathéodory function F : R k × Ξ → R, such that the expected value function f (x) := EF (x, ξ) is finite, and the sample average function fn (x) = n -1 n i=1 F (x, ξ i ). Take X to be in K(R d ) and set F ξ = F (•, ξ). We assume that F ξ is Lipschitzian on a neighbourhood of X , and denote its Clarke generalized gradient at x ∈ X by ∂F ξ (x).

With this notation, we have the following variant of [6, Theorem 5] with strengthened convergence. Note that we are applying Theorem 4 with ∂F ξ i (ω) (x) in the role of X i (x, ω).

1 A c c e p t e d m a n u s c r i p tf

 1 in non-convex non-smooth stochastic optimization. Stationary points of the problem minx Ef (x, ξ) (x, ξ i ),

Theorem 4 .

 4 Let X be a random u.s.c. multifunction whose range is essentially H-separable, and let {X n } n be p.i.i.d. as X. If there exists an integrable function κ such thatX(x, ω) ≤ κ(ω) for all (x, ω) ∈ X × Ω, d H (S n (x), E(x)) → 0 almost surely.

A c c e p t e d m a n u s c r i p t

  ,y)≤r S n (y) → 0.Clearly, those statements become weaker as r increases. The convergence we prove is equivalent to taking r = 0 above.Shapiro and Xu [6, Theorem 4] get to prove convergence in the case r = 0 by assuming additionally that X is almost surely continuous. But then its range is separable, because the whole space of continuous functions from X to the metric space K(R d ) is separable. Thus Theorem 4 includes[START_REF] Shapiro | Uniform laws of large numbers for set-valued mappings and subdifferentials of random functions[END_REF] Theorem 4]. Remark 1. From the proof of Theorem 2 one can perceive that upper semicontinuity is used only to ensure that x∈X u(x) is compact for u ∈ U, therefore it can actually be replaced by the latter condition, which is not linked to continuity anymore. Thus continuity properties seem not to be really essential to this problem, contrary to what might seem from Shapiro and Xu's approach. Remark 2. On the other hand, [6, Corollary 1] considers the point-valued case in which X(x, •) is almost surely a singleton, proving sup x∈X e(S n (x), E(x)) → 0 a.s.

Theorem 5 .Remark 3 .

 53 With the definitions above, if the following conditions are met: (i) There exists an integrable functionκ : Ξ → R such that |F (x, ξ) -F (y, ξ)| ≤ κ(ξ)|x -y| for all x, y ∈ X , ξ ∈ Ξ, F (•,ξ) is regular at every x ∈ X , for almost every ξ ∈ Ξ, (iii) The range of ξ → ∂F ξ is essentially separable; then sup x∈X d H (∂ fn (x), ∂f(x)) → 0 almost surely. This covers [6, Theorem 6] as well.

  The notation [λ 1 , u 1 ; . . . ; λ n , u n ] and its shorthand [λ i , u i ] n i=1 convey the usual meaning of λ 1 u 1 + . . . + λ n u n and

	n i=1 λ i u i but E is not assumed to have an addition.
	Terán and Molchanov [7] call E a convex combination space if the following
	axioms are satisfied:

Acknowledgements: I am grateful to Alexander Shapiro and Huifu Xu for their kind remarks. Research has been partially funded by the Spanish Ministerio de Educación y Ciencia, under its research grants MTM2005-02254 and TSI2005-02511, and the Gobierno de Aragón, under its research grant PM2004-052. This paper was conceived during a 6-hour stay at the Reina Sofía airport in Tenerife.


Simulation Optimization of Pull Control Policies for Serial Manufacturing Lines and Assembly Manufacturing Systems with Genetic Algorithms

Keywords: ASSEMBLY SYSTEMS, GENETIC ALGORITHMS ASSEMBLY SYSTEMS, GENETIC ALGORITHMS JIT manufacturing systems, pull production control policies, genetic algorithm

Several efficient pull production control policies for serial lines that implement the lean/JIT manufacturing philosophy can be found in the production management literature. A recent development that is less well-studied than the serial line case is the application of pull type policies to assembly systems where manufacturing operations take place both sequentially and in parallel. Systems of this type contain assembly stations where two or more parts from lower hierarchically manufacturing stations merge in order to produce a single part of the subsequent stage. In this paper we extend the application of the Base Stock, Kanban, CONWIP, CONWIP/Kanban Hybrid and finally the Extended Kanban production control policies to assembly systems that produce final products of a single type. Discrete-event simulation is utilized in order to evaluate the performance of serial lines and assembly systems. It is essential to determine the best control parameters for each policy when operating in the same environment. The approach that we propose and probe for the problem of control parameter selection is that of a genetic algorithm with resampling, a technique used for the optimization of stochastic objective functions. Finally, we report our findings from numerical experiments conducted for two serial line simulation scenarios and two assembly system simulation scenarios.

Introduction

Effective production control, which translates to the control of the material flow that takes place within a manufacturing system in a way that resolves the trade-off between a high service rate and a low holding costs level plays a key role to the competitiveness of any manufacturing system. The recognition of the importance of effective production control by the research community gave rise to a significant number of sometimes quite disparate production control mechanisms. However, many production control mechanisms for flow-shop systems fall under two fundamental categories: push type and pull type. Production control schemes that utilize future demand forecasts to control the manufacturing process are known as push type production control mechanisms. On the other hand, pull type control policies react to actual occurrences of demand only. Various definitions of push and pull (as well as a clarification of the meaning of these two terms) can be found in [START_REF] Pyke | Push and pull in manufacturing and distributions systems[END_REF]. Pull type control policies implement the JIT manufacturing philosophy and have attracted considerable attention over the past years as they are widely considered to outperform push type mechanisms. According to the JIT manufacturing philosophy a manufacturing system should maintain the minimum levels of safety stocks that are required in order to meet the target service level while in the same time has the ability to react rapidly to incoming orders. Two fundamental pull control policies are Kanban and Base Stock [START_REF] Buzacott | Stochastic Models of Manufacturing Systems[END_REF]. Generalized Kanban (see [START_REF] Buzacott | A general approach for coordinating production in multiple cell manufacturing systems[END_REF] for example) and Extended Kanban [START_REF] Dallery | Extended kanban control system: combining kanban and base stock[END_REF] are two policies that were developed in an attempt to combine the merits of Kanban and Base Stock. [START_REF] Spearman | CONWIP: a pull alternative to Kanban[END_REF] introduced CONWIP, a policy that can be seen as a special case of the Kanban control policy. Applying CONWIP control in the manufacturing system as a whole and Kanban control in all of the system's manufacturing stages except the last one resulted in a new hybrid control policy known as CONWIP/Kanban Hybrid (see [START_REF] Paternina -Arboleda | Intelligent dynamic control policies for serial production lines[END_REF].

The bibliography regarding serial manufacturing lines controlled by the host of the above mentioned policies is quite voluminous. However, many important pull type control policies can be also applied to manufacturing systems containing assembly stations. An assembly station is a manufacturing cell where two or more parts from upstream manufacturing stations merge in order to produce a single part of the subsequent stage. Up to now the research carried out for systems of this type is relatively limited. Two variants of the Kanban control policy for assembly systems are analyzed in Di Mascolo and [START_REF] Di Mascolo | Performance evaluation of kanban controlled assembly systems[END_REF] and [START_REF] Matta | Analysis of assembly systems controlled with kanbans[END_REF]. The Extended Kanban for assembly systems was presented in [START_REF] Chaouiya | The extended kanban control system for production coordination of assembly manufacturing systems[END_REF]. [START_REF] Takahashi | Comparing CONWIP, synchronized CONWIP and Kanban in complex supply chains[END_REF] study assembly systems controlled by Kanban, CONWIP and a variation called Synchronized CONWIP using simulation. In this paper we extend the application of the Base Stock, Kanban, CONWIP, CONWIP/Kanban Hybrid and finally the Extended Kanban production control policies to assembly systems that produce final products of a single type. The release of component parts in an assembly station is simultaneous. In this case, a component part that it is stored in its corresponding output buffer and has the controller authorization to enter the assembly station will have to wait until all the other components become available too. This required synchronization intensifies the queueing effects that appear between consecutive manufacturing/assembly stations.

The system performance of serial manufacturing lines and assembly manufacturing systems can be evaluated by discrete-event simulation [START_REF] Law | Simulation Modelling and Analysis[END_REF] or via queuing-model techniques (for example see [START_REF] Manitz | Queueing-model based analysis of assembly lines with finite buffers and general service times[END_REF] and [START_REF] Matta | Analysis of assembly systems controlled with kanbans[END_REF]). In this paper we use simulation to evaluate the major performance measures of the system such as average throughput, average WIP and average level of backordered demands. Pull production control policies are heuristics that are characterized by a small number of control parameters that assume integer values. The performance of a system operating under a certain pull type control policy is significantly affected by parameter selection. Therefore, it is obvious that control parameter selection is a task of grave importance when it comes to the design of a pull type manufacturing system. It is reasonable that we would like to set the control parameters of each policy when operating in the same environment to their best values in order to gain insight regarding their behavior and proceed to performance comparisons. The approach that we propose and probe for the problem of control parameter optimization is that of simulation optimization, a practice that is rapidly becoming one of the most wide-spread tools in the field of complex systems design/analysis. Indicatively, we refer the reader to [START_REF] Swisher | A survey of simulation optimization techniques and procedures[END_REF] for a survey on simulation optimization bibliography. In this paper, we interface the simulation model of the system we wish to examine with a genetic algorithm. Related work can be found in [START_REF] Dengiz | Simulation optimization using tabu search[END_REF] as well as in [START_REF] Spinellis | A simulated annealing approach for buffer allocation in reliable production lines[END_REF] where simulation is used in conjunction with tabu search and simulated annealing respectively, in order to determine the optimum design parameters of manufacturing/inventory systems. [START_REF] Finke | Multiple machine JIT scheduling: a tabu search approach[END_REF] also use tabu search to solve a variation of the classic flow shop system where machines are organized in tiers. The choice of a genetic algorithm as the optimization method was made on the basis of the prevalent, but not as well-studied in the literature as one might expect conception, that genetic algorithms are inherently well-suited for "noisy" environments. Due to the fact that the objective function that we wish to maximize is stochastic as its values are generated by the simulation model we use resampling, i.e., performing multiple evaluations of the same parameter vector and using the mean of these evaluations as the fitness measurement of this individual, a practice discussed by [START_REF] Fitzpatrick | Genetic algorithms in noisy environments[END_REF] and [START_REF] Hammel | Evolution strategies on noisy functions. How to improve convergence properties[END_REF].

The remaining material of this paper is structured as follows. Sections 2 -2.1 give a brief description of six important pull production control policies for serial manufacturing lines. Sections 3-3.5 are devoted to the presentation of the Kanban, Base Stock, CONWIP, CONWIP/Kanban Hybrid and Extended Kanban control policies for assembly manufacturing systems. In sections 4-4.1 we discuss the main aspects of the simulation optimization methodology that we followed and namely, the construction of the simulation model, the formal definition of the parameter optimization problem and issues concerning the genetic algorithm that was used. We report our findings from the simulation experiments that we conducted for two serial line and two assembly system cases in sections 5-5.2. Finally, in section 6 we state our concluding remarks and point to possible directions for future research.

System description of a serial manufacturing line

A typical manufacturing line produces a single product type and consists of a number of manufacturing stations with intermediate buffers. Raw materials undergo operations such as machining, forming etc., and are gradually converted into final products to be delivered to customers. Raw materials are assumed to be available to the system at any time, so the first machine is never starved. Customer demands that arrive to the system request the release of one finished part from the finished goods buffer. Whenever a demand arrives to the system it is immediately satisfied from the finished parts inventory. If there are no parts available in the last buffer the demand is backordered. Manufacturing facilities have the ability to work on only one part at a time during a manufacturing cycle. As soon as a stage i part is manufactured, it is placed in the output buffer of that stage. A control policy determines when a stage i finished part is to be released from the output buffer of that stage to manufacturing facility is available at that time. There is no delay in material handling between stages. All machines have random production time, time between failures and repair time. Time intervals between demand arrivals are also stochastic. Each pull production control policy coordinates the release of parts from one stage to another in a different way. In other words, the timing with which the information of a demand arrival is transmitted from the last production stage to all of the upstream stages is policy -dependent. The following sections briefly explain the way that the Kanban, Base Stock, CONWIP, CONWIP/Kanban Hybrid, Extended Kanban and Generalized Kanban control policies for serial lines operate. A serial line consisting of three manufacturing stations is depicted in Figure 1.

Base Stock control policy

A Base Stock manufacturing line is completely described by N parameters, the base stock levels i B of each production station,

N i ,..., 2 , 1 =
, where N is the number of the system's manufacturing stations. The i B parameters correspond to the number of parts that exist in the system's buffers at the time the system is in its initial state that is before any demands have arrived to the system. This control policy operates as follows.

When a demand arrives to the system it is immediately transmitted to every manufacturing station, authorizing it to start working on a new stage i part. Base Stock has the advantage of reacting rapidly to incoming demand, with the drawbacks of providing no upper bounds for the system's buffer levels, as well as the relatively loose coordination of the manufacturing stages.

Kanban control policy

A Kanban manufacturing line's control parameters are the production authorizations i K of each station, N i ,..., 2 , 1 = . The Kanban policy offers very tight synchronization between the various production stations of the system. The i K parameter corresponds to the maximum number of parts that are allowed in station i (manufacturing facility -output buffer). Production station i is authorized to start working on a new part as soon as a finished station i part is released from its output buffer. The information of a demand arrival is transmitted from the last manufacturing station to the first one station -by -station. If there is a point where a part is not available this transmission is interrupted. Due to this mechanism of back-propagating the demand arrival information, Kanban systems have relatively long response times to costomer orders.

CONWIP control policy

CONWIP is an abbreviation for CONstand Work In Process. According to this policy the total number of parts that exist in the system (Work In Process) can never exceed a certain level, which is the single control parameter of the CONWIP policy. All machines in a CONWIP line are authorized to produce whenever they have this ability (they are operational and have a raw part to work on) except the first one. The first machine of the system is authorized to start working on a new part as soon as a unit from the finished parts buffer is released to a customer.

Generalized Kanban and Extended Kanban control policy

These two control policies combine the merits of Base Stock and Kanban as they react rapidly to the arrival of demands and effectively control the WIP at the same time. They are described by two parameters per stage, the base stocks i B and the production authorizations ) ( ,

i i i B K K >
, that are borrowed from the Base Stock and Kanban policies respectively. The finite number of production authorizations guarantees that the system's inventories will not exceed the pre -defined levels, but the stage coordination here is not as tight as in Kanban.

A station can be granted a production authorization even if a part is not released from its output buffer. For a detailed description of the way Generalized Kanban and Extended Kanban operate the reader is referred to [START_REF] Liberopoulos | A unified framework for pull control mechanisms in multi -stage manufacturing systems[END_REF], [START_REF] Dallery | Extended kanban control system: combining kanban and base stock[END_REF] and [START_REF] Buzacott | A general approach for coordinating production in multiple cell manufacturing systems[END_REF].

CONWIP/Kanban Hybrid

A CONWIP/Kanban Hybrid system, as implied, operates under a combination of the CONWIP and Kanban control policies. Departure of a finished part from the system authorizes the first stage to allow a new raw part to enter the system. All stages except the last one have a finite number of production authorizations 1

,..., 2 , 1 , - = N i K i
. Stage production authorizations i K and the total WIP that is allowed in the system are CONWIP/Kanban Hybrid's control parameters.

System description of a manufacturing system with assembly stations

In this research paper we examine manufacturing systems with assembly stations that produce a single type of final products. A manufacturing station has only one immediate upstream system station and in its turn, feeds only one subsequent manufacturing or assembly station. On the other hand, an assembly station has two or more immediate upstream stations that supply it with the component parts that will merge in order to produce a single part. An assembly station also gives its finished parts to only one immediate downstream station. We limit ourselves to the case where component parts are released in an assembly station simultaneously (see [START_REF] Chaouiya | The extended kanban control system for production coordination of assembly manufacturing systems[END_REF] and Di Mascolo and [START_REF] Di Mascolo | Performance evaluation of kanban controlled assembly systems[END_REF] for another variant). This implies that a component part that is stored in output buffer i and has the authorization to enter an immediately downsteam assembly station will have to wait until all the other components become available too. Since this part has not exited buffer i it is incorporated in the estimation of buffer's inventory level. An unlimited supply of raw materials is assumed for this type of system too. Customer demands are satisfied instantaneously from the finished goods stock. If there are no parts available in the last buffer the demand is backordered. Manufacturing facilities and assembly stations have the ability to work on only one part at a time during a manufacturing cycle. As soon as a stage i part is finished, it is placed in the output buffer of that stage. A control policy determines when a station i finished part is authorized to be released from the output buffer of that station to the next manufacturing/assembly facility. All stations in the system have random operation time, time between failures and repair time. Time intervals between demand arrivals are also stochastic. The following sections describe the Kanban, Base Stock, CONWIP, CONWIP/Kanban Hybrid and Extended Kanban control policies for assembly manufacturing systems. The systems are depicted as queuing networks with synchronization stations, a notion borrowed by [START_REF] Liberopoulos | A unified framework for pull control mechanisms in multi -stage manufacturing systems[END_REF]. In our opinion this is the most illustrative and easy to comprehend modeling framework. A synchronization station with four input queues and three outgoing customers is shown in Figure 2. The perpendicular line in Figure 2 symbolizes a server with instant service times that is fed by queues a, b, c and d. As soon as there is at least one customer in each of the input queues, these customers enter the server and instantaneously, one customer of type e, one of type f and one of type g exit the server. For additional details the reader is referred to [START_REF] Liberopoulos | A unified framework for pull control mechanisms in multi -stage manufacturing systems[END_REF]. For simplicity reasons, all systems in the following sections consist of two manufacturing stations and one assembly station, however everything which holds for the dynamics of each control policy for this small system can be easily extended for larger systems.

Base Stock control policy for assembly systems

Figure 3 represents the model of an assembly system controlled by the Base Stock policy with two manufacturing stations and one assembly station. In our model

3 , 2 , 1 , = i O i that contain 3 , 2 , 1 , = i B i stage i finished parts. The integers 3 , 2 , 1 , = i B i
that will be hereafter called as the base stocks of the stations are the system's only control parameters.

Kanban control policy for assembly systems

A Kanban system with two manufacturing stations in parallel that feed an assembly station is illustrated in Figure 4.

3 , 2 , 1 , = i S i
, stands for the system's manufacturing/assembly station i while queues represent the input buffers of stations i =1, 2, 3. We remind that input buffers are not physical areas of storage (see section 3.1). When the system is in its initial state, that is before any demands have arrived to the system, all queues in the model are empty except of the raw material buffers

3 , 2 , 1 , = i OK i that contain 3 , 2 , 1 , = i K i
stage i finished parts. Each station in the system has a fixed number of production authorizations or kanbans which are the system's control parameters. In other words, the inventory in all of the system's output buffers

3 , 2 , 1 , = i OK i cannot exceed the predefined levels. The kanbans of the stations 3 , 2 , 1 , = i S i are equal to the numbers 3 , 2 , 1 , = i K i
of parts that exist in the output buffers initially. The material and information flows are not independent, as they are in the Base Stock policy. The arrival of a demand to the system "pulls" a part from the output buffer of station 3 DK . This way the information of a customer arrival is propagated backwards through the release of kanbans. If at some point of the system, there is an empty output buffer, then this propagation is ceased. It should have become clear by now that the Kanban control policy for assembly systems offers very tight coordination between the system's manufacturing/assembly stations.

CONWIP control policy for assembly systems

CONWIP/Kanban Hybrid for assembly systems

A CONWIP/Kanban system with two manufacturing stations in parallel that feed an assembly station is illustrated in Figure 6

Extended Kanban control policy for assembly systems

O . Queues , 3 , 2 , 1 , = i FK i contain , 3 , 2 , 1 , = -i B K i i
kanban authorizations and all other queues in the system are vacant. The control parameters of the Extended Kanban system are the total station i production authorizations i K and the base stocks

i B , . 3 , 2 , 1 = i As one may see, in the extreme case of , 3 , 2 , 1 , = = i B K i i
the system is reduced to simple Kanban assembly system. The information of a customer arrival is transmitted instantly to all stations, a property borrowed from the Base Stock policy. The Extended Kanban control policy has a fixed number of kanbans in order to constrain the WIP, but in the same time offers more loose coordination between manufacturing/assembly stations than the original Kanban policy. By examining Figure 7, it becomes clear that a station can be granted an authorization to start working on a new part even if no parts have exited its ouput buffer.

Simulation optimization of serial manufacturing lines and assembly manufacturing systems

We use discrete-event simulation to evaluate the three major performance measures of serial manufacturing lines and assembly systems which are average throughput (number of final products delivered to customers per time unit), average WIP (WorkInProcess) and average level of backordered demands. Pull production control policies are simple heuristics that are characterized by a small number of control parameters that assume integer values. We seek to set the control parameters of each policy when operating in the same environment to the values that maximize a properly defined objective function in order to proceed to performance comparisons. To serve this purpose we make use of simulation optimization which, according to Hall and Bowden (1998), is "... the practice of linking an optimization method with a simulation model to determine appropriate settings of certain input parameters so as to maximize the performance of the simulated system". In this paper we interface the simulation model of the serial lines and assembly systems we wish to examine with a genetic algorithm.

Simulation of manufacturing serial lines and assembly systems

Discrete -event simulation, [START_REF] Law | Simulation Modelling and Analysis[END_REF], is the modeling tool that we used in order to simulate the operation of serial lines and assembly systems. According to this methodology, the simulation program is constructed on the basis of the discrete events that are possible to take place in a given system. In the case of an assembly system or a manufacturing serial line, we have the following events.

1 Since we examine pull type production control policies the event that appears first in our model is a demand arrival. The discrete -event simulation algorithm of a manufacturing serial line or an assembly system that operates under a pull production control policy is presented in Listing 3.1. Initialization. Manufacturing system's topology, service times, failure times, demand times etc., duration of the simulation sim T , t = 0, b t = 0, where t is the current simulation time and b t the previous simulation event's time of appearance.

WHILE ( t < sim T )
Find minimum time of occurrence from the next events listing Event execution. The system's state variables are adjusted with the use of its state equations Controller call. Determine whether the station i will be operational or idle during the next simulation cycle.

Where necessary, update the next events time of occurrence Listing 3.1. Manufacturing system discrete -event simulation algorithm

Optimization problem -Objective function

The mathematical formulation of the parameter optimization problem for serial lines and assembly systems controlled by pull production control policies is given below.

Let x = [ 1 x 2 x .... n x ], ∈ i
x Z, be the control parameter vector of some pull production control policy. For example, x could be the station i production authorizations (kanbans) in a Kanban system, the initial buffer levels in a Base Stock system etc. The objective is to find the control parameter values x that maximize the expected value of the stochastic objective function ( F x, ω), where ( F x, ω) is the derived profit from the system's operation per time unit.

Maximize:

( [F E x, ω)]
(1) subject to:

∈ i x Z, n i ,...., 2 , 1 =
ω is used to denote the stochastic nature of F. The evaluation of function F is the result of a simulation experiment. The value of F is calculated according to 2.

∑ = - - = N i i i B b H h pTH F 1 (2)
F is the derived profit from the system's operation per time unit, p the profit from a finished part's sale to a customer, TH the system's average throughput which is calculated by dividing the number of parts the system's last station produced by the simulation time sim T . h i stands for the cost of storing one part in output buffer i per time unit while i H is the average inventory level in output buffer i. The cost of one unsatisfied demand per time unit is denoted as b and finally, B is the system's average number of backordered demands. Repeating a simulation for the same value of x would result in different values of F due to the stochasticity of the simulation model that is denoted by ω.

Genetic algorithm

Genetic algorithms are stochastic iterative search techniques that do not require any kind of information on the objective function's structural properties as the search direction is determined solely by means of objective function evaluation. A genetic algorithm generates a population of individuals-solutions at each iteration. The next population is created by using computations that involve stochastic selections. Due to the fact that the objective function that we wish to maximize is stochastic we use resampling, a practice that is discussed in [START_REF] Beyer | Evolutionary algorithms in noisy environments: theoretical issues and guidelines for practice[END_REF], [START_REF] Fitzpatrick | Genetic algorithms in noisy environments[END_REF] and [START_REF] Hammel | Evolution strategies on noisy functions. How to improve convergence properties[END_REF]. For a given

( [F E x, ω)]= ∑ = m i m 1 1 ( F x, ω) i , x = const (3)
This technique is perhaps the most straight-forward when it comes to evaluation of stochastic functions but it can be very expensive computationally as the calls to the simulation model increase by a factor of m. On the other hand, the higher the value of m is, the more accurate the evaluation of ( [F E x, ω)] becomes. Another factor that affects the quality of the estimate of ( [F E x, ω)] is the simulation time sim T for which we run the model. We determine the values of m and sim T by trial and error in order to have a fairly accurate estimate of ( [F E x, ω)] without exhaustively long simulation experiments. This is a very delicate point as estimates that are not sufficiently accurate may cause the algorithm to incorrectly interpret improvements caused by the stochasticity of the simulation model as meliorations caused by a change in the parameter set or fail to detect better solutions that cause minor improvements in the objective function. The fact that we are dealing with a stochastic objective function must be taken into consideration when setting the stopping criterion too. Our algorithm terminates when a number of g successive generations do not yield more than a p% improvement on the objective function. It is essential to determine roughly through experimentation the range of ( [F E x, ω)]s that a given x will produce in order to set parameter p appropriatelly. This way we ensure that the algorithm will not be trapped in an endless loop and eventually terminate at some point. In cases where we suspect that a more thorough search of the decision space is needed, as a complementary stopping criterion, we prevent the algorithm from stopping before it has produced a minimum number of generations. The remaining parameters of the genetic algorithm are the standard ones of population size s, crossover probability p cross and mutation probability p mut .

Results -Discussion

Serial manufacturing line results

We examined a four -stage manufacturing line with equal operation times. Machines operate with service rates which are normally distributed random variables with mean 1.0 parts/time unit and s.d. 0.1 (

) 1 . 0 , 0 . 1 ( ~N R p ).
Repair to failure times are exponentially distributed with mean 1000.0 time units. Failures are operation dependent. Repair times are also assumed exponential with a MTTR of 10.0 time units. We focused our attention on two cases regarding the demand arrival process; one where the system has to satisfy a moderate demand for finished goods and one where the demand arrival rate is higher and close to the maximum throughput of the given system. In the first case, times between two successive customer arrivals are exponential random variables with mean 1.66 time units (0.6 customers per time unit), while in the second case the mean time between customer arrivals is equal to 1.25 time units (0.8 customers per time unit). Simulation time is set to 30000.0 time units. The sale of one finished part to the customers results in a 120,0 units profit. Inventory costs for storing one part per time unit in buffer i are shown in Table 1.

Table 1 should be placed about here

Linear backorder costs are incurred at a rate of 10.0 whenever there are backorders. The authors conducted a series of experiments in an attempt to come up with the most suitable genetic algorithm parameter set for this specific problem. The genetic algorithm's parameters that appeared to yield the best results were: population size = 40, crossover probability = 0.5, mutation probability 0.05, parameter p = 5% and g = 40 (see section 4.3). Each individual-parameter vector was evaluated 30 times, (m=30). In some cases we used the complementary stopping criterion in order to achieve a more thorough search of the solution space with the minimum number of generations to be produced by the genetic algorithm set to 100. The best parameters of each policy for the first simulation case are presented in Table 2. We denote base stock in buffer i by writing 3 presents the average WIP distribution through out the system and the backordered demands average together with the corresponding expected values of F. These figures were calculated by averaging the output of 30 executions of each policy's simulation model, where every model was allowed to run for 300,000.0 time units, that is ten times the simulation time used for fitness evaluation in the genetic algorithm.

Table 2 should be placed about here Table 3 should be placed about here By reading the results presented in Table 2 we observe that the genetic algorithm converged to relatively low values for the control parameters of all pull production control policies. Especially for the case of the Kanban control policy, the i K s were set to the lowest possible values these parameters can assume (if a station in a Kanban system has zero production authorizations then this station would never produce a single part). Note that by increasing the number of kanbans, base stocks or CONWIP type authorizations, the average throughput of the serial line increases but so does the WIP along with the related costs. The major cost factors in our serial line simulation scenario are the storage costs, and primarily the inventory level in the system's last buffer, as the average of backordered demands is weighted moderately. Given the production rates of the machines and the rate of incoming orders, we could infer the system can easily handle customer arrivals without having to maintain large safety stocks and this is depicted in the parameter sets that the algorithm converged to. The characteristic attribute of the Kanban policy which is tight coordination between manufacturing stations is evident in the first row of Table 3. We observe that the WIP is distributed practically evenly between buffers 1 to 4. However the average WIP in buffers 1, 2 and 3 are the highest among all control policies. The Base Stock control policy achieves a significantly leaner operation in expense of the lower service level. In a system that operates under the CONWIP policy all the machines except the first one have the perpetual authorization to produce whenever they can. There is a constant flow of parts through the system that ends in the finished goods buffer where parts are accumulated, and that is why the average finished products inventory is the highest of all control policies. As a consequence, parts are not stalled in intermediate buffers due to e.g. lack of production authorization which in turn results in very low average WIP levels in the first three buffers as it is evident when observing Table 3. The existence of additional station i production authorizations in the CONWIP/Kanban Hybrid policy does not alter the system's performance in a dramatic way, as one may observe in Table 3. The application of the Generalized Kanban control policy resulted in the highest profit rate in our simulation cases, while Extended Kanban ranked second. Both Generalized Kanban and Extended Kanban offer effective WIP control and maintenance of a high service rate. These two control policies have a finite number of production authorizations in order to prevent the system's Work In Process inventories from growing unboundedly but the stage coordination is not as tight as in Kanban. This is a result of the

n i B K i i ,... 2 , 1 , = -
, initially free kanban cards that Generalized Kanban and Extended Kanban make use of. Examination of the best parameters table (Table 2) shows that this technique allows the allocation of a greater number of production authorizations than simple Kanban that speed up the system's reaction to demand occurrences, while the base stocks approximate those used in the Base Stock system. The best parameters of each policy found by the genetic algorithm for the second simulation case, (demand arrival rate=0.8 arrivals/time unit), are presented in Table 4. Figure 9 shows the fitness of the best solution found by the genetic algorithm in each iteration for every control policy versus the number of generations. 5 presents the average WIP distribution through out the system and the backordered demands average together with the corresponding expected values of F. Again, these figures were calculated by averaging the output of 30 executions of each policy's simulation model, where every model was allowed to run for 300,000.0 time units. As we expected, the genetic algorithm selected higher values for the control policies' parameters, due to the need for increased safety stocks of intermediate and final products, caused by the higher demand rate. What we are primarily interested in, is to determine whether the increase of the demand rate modifies the ranking of the six production coordination mechanisms reported for the first simulation scenario. The Base Stock policy still manages to outperform the Kanban policy. An interesting point is that this happens for entirely different reasons in the two simulation cases. In the first, low demand rate case, the Base Stock policy outperforms Kanban due to its lower average inventories caused by the zero installation stock in the first three buffers while in the second, high demand rate case, the same happens but this time due to its approximately 43% lower level of backorders. A possible explanation for this is that the selection of the objective function's weights (holding cost of one stage i part per time unit, cost of one backorder per time unit) favors the Base Stock policy over the Kanban policy. The Generalized and Extended Kanban systems rank in the first two positions in both simulation cases giving strong indications that they can excel in various demand rate environments. Perhaps the most interesting observation that stems from this second series of optimization experiments is the striking improvement in the performance of the CONWIP and CONWIP/Kanban Hybrid control policies. More specifically, the CONWIP policy ranked last in the low demand simulation case, while in the high demand case, it exhibits a performance comparable to those of the two-parameter-per stage policies. The inherent characteristic of the CONWIP policy, which is uninterrupted flow of materials through the system and WIP accumulation in the finished goods buffer, is what renders it highly sensitive to demand fluctuations for these two simulation cases. However, notice that the CONWIP and the Hybrid policies would probably not perform that well in the high demand rate case if the serial line was not balanced, as it is in our study. For examble, consider a situation where the bottleneck of the system is located, i.e. in the second machine due to fact that this machine is more susceptible to failure than the other machines. In that case, the WIP would tend to accumulate in the first buffer, leading to excessive storage costs.

Assembly system results

We examined an assembly system that consists of four manufacturing stations and three assembly stations. The topology of the system is shown in Figure 10. ). Repair to failure times are exponentially distributed with mean 1000.0 time units. Failures are operation dependent. Repair times are also assumed exponential with a MTTR of 10.0 time units. Similarily to the serial line case, we examined two scenarios regarding the demand; in the first scenario times between two successive customer arrivals are exponential random variables with mean 1.66 time units (0.6 customers per time unit), while in the second scenario the mean time between customer arrivals is equal to 1.25 time units (0.8 customers per time unit). Simulation time is set to 40000.0 time units. The sale of one finished part to the customers results in a 70.0 units profit. Inventory costs for storing one part per time unit in buffer i are shown in Table 6.

Table 6 should be placed about here Linear backorder costs are incurred at a rate of 10.0 whenever there are backorders. Significant emphasis is given to customer service level and this is done through the inventory/backordered demands costs scheme. The cost per unsatisfied demand has a much larger impact on the objective function's value, than the cost of a part stored in buffer i. The parameters of the genetic algorithm were: population size = 30, crossover probability = 0.5, mutation probability 0.05, parameter p = 5% and g = 40 (see section 4.3). Each individual-parameter vector was evaluated 30 times, (m=30). In some cases we used the complementary stopping criterion in order to achieve a more thorough search of the solution space with the minimum number of generations to be produced by the genetic algorithm set to 100. The best parameters for each policy for the low demand rate case are presented in Table 7. We denote base stock in buffer i by writing i 8 presents the average WIP distribution through out the assembly system and the backordered demands average together with the corresponding expected values of F. These performance measures were calculated by averaging the output of 30 executions of each policy's simulation model, where every model was allowed to run for 300,000.0 time units.

Table 7 should be placed about here Table 8 should be placed about here In the assembly system of our simulation scenario component parts are released in an assembly station simultaneously (see sections 3.1-3.5). This means that a component part that is stored in its corresponding output buffer and has the authorization to enter the assembly station will have to wait until all the other components become available too. The simultaneous release of component parts in the assembly station can be treated as a part of the control policy or as a limitation posed due to technical reasons, i.e., the assembly operation cannot start until all components are gathered. Whichever the case may be, this required synchronization intensifies the queuing effects that appear between consecutive manufacturing/assembly stations. The Kanban and the Base Stock control policies rank in the two last places by the expected profit ] [F E criterion. The average number of backordered demands that the Kanban policy achieves is the highest among all competing policies. Evidently, the mechanism of authorizing production in the system's workstations through the release of kanban cards has a negative effect on the system's response time to incoming customer orders. On the other hand, the Base Stock policy yields the lower average of backorders, in the expense of relatively poor WIP control. CONWIP is one of the simplest control policies but this does not prevent it from excelling in the task of controlling an assembly system. WIP is constrained effectively in all storage facilities except the last one, something which characterizes this policy, while in the same time, incoming orders are satisfied almost instantly (average of backorders=0.86). The CONWIP/Kanban control policy for assembly systems performs worse than the CONWIP policy regarding the WIP in the output buffers of the manufacturing stations 1, 2, 3 and 4. Far betters results are obtained for output buffers 5 and 6, where the WIP is kept at the lowest levels among all competing policies throughout the simulation period. The Extended Kanban control policy yielded the highest expected profit rate in our assembly system simulation scenario, however given the stochasticity of the objective function we could argue that its performance is practically the same with that of the CONWIP mechanism. The statistics of this policy's performance displayed in Table 8 indicate that the existence of

n i B K i i ,... 2 , 1 , = -
, initially free kanban cards in the Extended Kanban assembly system allow the unhampered flow of information within the system's bounds which gives it the ability to react rapidly to demand occurrences without the need for excessive safety stocks. The best parameters for each policy for the second simulation case, (demand rate =0.8) are presented in Table 9. Figure 12 shows the fitness of the best solution in each iteration found by the genetic algorithm for every control policy versus the number of generations. 10 presents the average WIP distribution through out the assembly system and the backordered demands average together with the corresponding expected values of F. These performance measures were calculated by averaging the output of 30 executions of each policy's simulation model, where every model was allowed to run for 300,000.0 time units. Unlike the serial line system examined in section 5.1 the ranking of the control policies is not altered significantly with the exception of the CONWIP and CONWIP/Kanban hybrid control policies which swap positions. However, the difference in the values of the objective function is marginal and could probably be attributed to the simulation itself. The reasons that we did not observe a change in the ranking of the policies can be found in the different inventory/backorder cost schemes which were used for the serial line and the assembly system. We remind that in the serial line case the major cost factor was the average inventory level of the last buffer while in the assembly system the most heavily weighted cost factor was the average number of backorders. It is interesting to highlight the differences between the parameter sets that the genetic algorithm selected for the two demand rate cases. With the exception of the CONWIP/Kanban hybrid policy, in all policies, the base stock in the last buffer is significantly lower in the moderate demand case than in the high demand case. The differences between the base stock levels of the two demand rate cases are substantially smaller as we move upstream. This indicates that the throughput, and as a consequence the number of backorders, is primarily affected by the initial inventory of the last buffer and that the impact of inventory in the system's throughput diminishes in upstream buffers.

Conclusion and future research

Effective production control is a topic that has attracted the research community's interest since many years as it is of great importance to any manufacturing system's competitiveness. In this paper we briefly present six important pull type production control policies that are used in manufacturing serial lines. Then we extend the use of five of them and namely Kanban, Base Stock, CONWIP, CONWIP/Kanban Hybrid and Extended Kanban to assembly systems with simultaneous release of components in the assembly stations. Each one of these control policies is completely described by a number of control parameters. Pull type control policies are simple heuristics so there are no significant difficulties in their implementation except that of the control parameter selection, a multi -dimensional optimization problem. As a measure of both the serial lines' and the assembly systems' performance we used a linear function of the system's average throughput, the average WIP and the average level of backorders. Simulation optimization is our approach to the parameter optimization problem. The optimization method that we used was a genetic algorithm that uses resampling, a choice that was dictated by the fact that the objective function we wanted to maximize was stochastic. We could argue that a genetic algorithm linked to a simulation model proves to be a valuable tool for the design of complex systems like the ones treated in this paper despite the fact that a considerable amount of fine-tuning is needed in order to determine appropriate settings for the genetic algorithm. We examined two serial manufacturing line simulation scenarios and two assembly system simulation scenarios. In all cases we used stochastic production times, times between failures, times to repair and demand. Based on the control parameter sets that we obtained using the genetic algorithm, we presented performance comparisons between the pull type control policies we presented. We could argue that the main qualitative attributes of all of the pull type control policies for serial lines seem to apply to systems that contain assembly stations too. The two-parameter-per stage policies, (Generalized and Extended Kanban), seem to be the more efficient ones regardless of the cost function weights, the demand rate and the existence or not of parallel machines in the system. The uncoupling of the information flow from the flow of parts in the system, offers the ability to opt for less base stock than the Kanban policy and achieves better response to demand occurrences. Our simulated data indicate that, when fine-tuned, the Generalized and Extended Kanban mechanisms can outperform less sophisticated mechanisms like the Kanban and the Base Stock. However, the increased number of parameters that need to be tuned is an important disadvantage that cannot be overlooked. A real-world manufacturing system can easily consist of, e.g. 20 machines, which translates to 40 parameters that need to be optimized. This poses a significant computational burden especially if high precision results are required. The CONWIP and CONWIP/Kanban Hybrid policies can approximate the performance of the Generalized and Extended Kanban mechanisms for certain settings as shown in sections 5.1 and 5.2, while being computationally cheaper. On the other hand, these two policies can perform poorly in situations where the demand rate is moderate and the cost of final product inventory is high. Pull type control policies as efficient as they can be, still remain heuristic, sub -optimal control mechanisms that do not receive any feedback in terms of some kind of reward signal on the produce/do not produce decisions they make. Here is where the Artificial Intelligence field comes into scope. Development of new and efficient, AI based production controllers is a challenging but most promising research topic, and some work has already been done to this direction. Another possible alternative would be the analysis and simulation of different kinds of manufacturing systems such as a multi -product manufacturing lines, pull control mechanisms that utilize Advance Demand Information etc.

Policy B 1 / K 1 B 2 / K 2 B 3 / K 3 B 4 / K 4 C 1 1 1 1 Kanban 1 1 1 1 - 0 0 0 2 Base Stock ∞ ∞ ∞ ∞ - 1 1 2 0 CONWIP ∞ ∞ ∞ ∞ 4 2 1 1 0 CONWIP/Kanban Hybrid 2 1 1 ∞ 4 0 0 1 1 Extended Kanban

  facility. The finished stage i part is released to the next stage provided that the stage 1 + i

Figure 1

 1 Figure 1 should be placed about here Figure 1. Three station manufacturing line

Figure 2 Figure 2 .

 22 Figure 2 should be placed about here Figure 2. Synchronization station with four input queues and three outgoing customers

Figure 3 D

 3 Figure 3 should be placed about here Figure 3. Base Stock assembly system with two manufacturing stations and one assembly station

Figure 4

 4 Figure 4 should be placed about here Figure 4. Kanban assembly system with two manufacturing stations and one assembly station

Figure 7

 7 Figure 7 should be placed about here Figure 7. Extended Kanban assembly system with two manufacturing stations and one assembly station

  vector x we measure ( F x, ω) m times and use the mean of these measures as the expected value of ( F x, ω).

iB

  and station i production authorizations by i K . CONWIP type production authorizations are represented by the symbol C. Figure8shows the fitness of the best solution found by the genetic algorithm in each iteration for every control policy versus the number of generations.

Figure 8

 8 Figure 8 should be placed about here Figure 8. GA max E[F] plots -serial manufacturing lines -demand rate = 0.6

Figure 9

 9 Figure 9 should be placed about here Figure 9. GA max E[F] plots -serial manufacturing lines -demand rate = 0.8

Figure 10

 10 Figure 10 should be placed about here Figure 10. Assembly system topology

  type production authorizations are represented by the symbol C. Figure11shows the fitness of the best solution in each iteration found by the genetic algorithm for every control policy versus the number of generations.

Figure 11

 11 Figure 11 should be placed about here Figure 11. GA max E[F] plots -assembly systems -demand rate = 0.6

Figure 12

 12 Figure 12 should be placed about here Figure 12. GA max E[F] plots -assembly systems -demand rate = 0.8

  Figures

  Figure 1. Three station manufacturing line

Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .PageFigure 9 .FigureFigure 11 .

 234567911 Figure 2. Synchronization station with four input queues and three outgoing customers

O

  Input buffers are not physical spaces in our model, they are merely used for illustrative reasons, and to be more specific, to flag a change in the state of finished station to input buffer i I it means that it has been authorized by the control policy to be processed in station i as soon as this becomes possible, i.e. the station might be occupied by another part or it could be under repair. Through the use of these virtual buffers we explicitly show the timing with which a station is given the order to start working on a new part regardless of the fact that this station could or could not carry out this order at that time. Material flow is symbolized by continuous lines while the flow of information is shown by dashed lines. Initially, all queues in the model are empty except of the raw material buffers 0

	whereas queue 4 D contains the customer demands for final products. Finally,	I i	, = i	, 1	, 3 , 2	stands for the input
	buffer of station i . 1 -i	parts. When a part
	"moves" from output buffer	i	1 -
							O and
	'
	3 , 2 represent the output buffers of the corresponding , 1 , = i S i , are the systems v 3 , 2 , 1 , = i , 3 O i , 2 e manufacturing/assembly stations while queues stations. By 0 O and ' , 1 , = i o r P e r R e i e w O n l D F y
							5

0

O we denote the raw materials storing facilities. By assumption these two queues cannot be totally empty at any time. Queues i contain the demands for a unit of a station i finished part 0 O , and the output buffers

  Figure5presents the model of CONWIP assembly system with two manufacturing stations and one assembly station.finished station i parts (base stocks) and all of the remaining queues in the system are empty. The system's control parameters are the maximum number C of parts allowed in the finished products buffer 3 O CONWIP assembly system with two manufacturing stations and one assembly station All stations in a CONWIP assembly system are authorized to produce whenever they have this ability except of the ones that are fed by the raw materials buffers. In Figure5these stations are DC is sent one CONWIP type production authorization as soon as one product exits buffer 3O . The WIP in assembly systems controlled by the CONWIP policy tends to accumulate in the finished products buffer.

			S i	, = i	, 1	3 , 2	, is the system's manufacturing/assembly station i while queues	O i	, = i	, 1	3 , 2	represent the
	output buffers of the corresponding stations. We denote the raw materials buffers that, by assumption, can
	never be depleted as 0 O and	' O . CONWIP type production authorizations are placed in queues 0	1 DC and
	2 DC . Queue 4 D contains the customer unsatisfied demands and	I i	, = i	, 1	, 3 , 2	stands for the input buffer of
	station i . The initial state of the CONWIP assembly system is the following: the output buffers contain
	B i	, = i	, 1	, 3 , 2
	(CONWIP type production authorizations), and the base stocks	B i	, = i	, 1	, 3 , 2	.
											Figure 5 should be placed about here
			Figure 5. 2 , 1 , = i S i	. These stations are
	allowed to start processing a new part by the moment that they are granted a CONWIP type production
	authorization (queues	1 DC and	2 DC ). Each one of the queues	1 DC and	2

  An Extended Kanban assembly system with simultaneous release of components (see Chaouyia et al. for another variant) in the assembly stations is depicted in Figure7. Queues =1,2,3. We remind that input buffers are not physical areas of storage (see section 3.1).Dashed lines symbolize information flow while solid lines represent material flow. (Here kanban authorizations/cards are treated as physical entities.). In the system's initial state, i

									S i	, = i	, 1	3 , 2	, are the systems
	manufacturing/assembly stations while queues	OK i	, = i	, 1	3 , 2	represent the output buffers of the corresponding
	stations 1,2,3 and queues 0 O and	' O are the raw material buffers. Queues 0	FK i	, = i	, 1	, 3 , 2	contain station i
	production authorizations (kanbans). Queues	D i	, = i	, 1	, 3 , 2	contain demands for station i finished products,
	queue 4 D holds the customer demands for final products and finally queues	I i	, = i	, 1	, 3 , 2	represent the input
	buffers of stations i B station i parts are stored in
	queues	OK i	, = i	, 1	3 , 2	and raw materials are placed in queues 0 O and	0	'

Table 4

 4

	should be placed about here

  Table 5 should be placed about here

Table 9

 9

	should be placed about here
	Table 10 should be placed about here

Table 1 .

 1 Inventory costs -serial line simulation cases 1-2.

	Tables
		Buffer 1	Buffer 2	Buffer 3	Buffer 4
	Inventory cost	10.0	10.0	15.0	20.0

Table 2 .

 2 Best control parameters -Serial line simulation case 1. (demand arrival rate = 0.6 )

Table 3 .

 3 Average inventories -backorders average -expected profit -serial line simulation case 1. (demand arrival rate = 0.6). The cells contain the mean of these performance measures and the standard deviation in parenthesis.

	Policy	Buffer 1	Buffer 2	Buffer 3	Buffer 4	Backorders average	[F E	]
	Kanban	0.39 (0.0)	0.37 (0.0)	0.36 (0.0)	0.34 (0.0)	1.14 (0.06) 40.87 (0.64)
	Base Stock 0.17 (0.02)	0.16 (0.02) 0.16 (0.02)	0.26 (0.0)	1.79 (0.04) 43.13 (0.83)
	CONWIP	0.05 (0.0)	0.04 (0.0)	0.04 (0.0)	1.25 (0.01)	0.9 (0.05) 36.33 (0.46)
	CONWIP/
	Kanban	0.09 (0.0)	0.04 (0.0)	0.02 (0.0)	1.24 (0.01) 0.94 (0.07) 36.12 (0.61)
	Hybrid
	Extended
	Kanban	0.1 (0.0)	0.08 (0.0)	0.17 (0.0)	0.17 (0.0)	1.81 (0.04) 46.01 (0.45)
	Generalize
	d Kanban	0.03 (0.0)	0.04 (0.0)	0.16 (0.0)	0.17 (0.0)	1.9 (0.05)	46.55 (0.48)

Table 4 .

 4 Best control parameters -Serial line simulation case 2. (demand arrival rate = 0.8 )

	Policy	B 1 / K 1 B 2 / K 2 B 3 / K 3 B 4 / K 4	C
	Kanban	1 1	2 2	1 1	1 1	-
	Base Stock	0 ∞	2 ∞	1 ∞	2 ∞	-
	CONWIP	0	1	2	2	5
		∞	∞	∞	∞
	CONWIP/Kanban Hybrid	1 1	2 2	2 2	1 ∞	6
	Extended Kanban	0 2	2 2	1 2	1 6	-
	Generalized Kanban	0 6	1 5	1 6	2 5	-

Table 5 .

 5 Average inventories -backorders average -expected profit -serial line simulation case 2. (demand arrival rate = 0.8). The cells contain the mean of the performance measures and the standard deviation in parenthesis.

	Policy	Buffer 1	Buffer 2	Buffer 3	Buffer 4	Backorders average	[F E	]
	Kanban	0.18 (0.0)	0.88 (0.01)	0.16 (0.01)	0.12 (0.0)	5.9 (0.32)	21.67 (2.91)
	Base Stock 0.71 (0.06)	0.97 (0.06)	0.80 (0.07)	0.45 (0.0)	3.38 (0.15)	24.23 (2.75)
	CONWIP	0.19 (0.0)	0.17 (0.00)	0.16 (0.0)	0.81 (0.01) 3.83 (0.28)	35.35 (2.67)
	CONWIP/
	Kanban Hybrid	0.06 (0.0)	0.11 (0.0)	0.11 (0.0)	1.16 (0.02) 4.11 (0.27)	28.28 (2.14)
	Extended Kanban	0.30 (0.01)	0.51 (0.0)	0.31 (0.0)	0.13 (0.0)	4.48 (0.30)	35.74 (2.93)
	Generalize d Kanban	0.3 (0.01)	0.34 (0.01) 0.34 (0.01)	0.32 (0.0)	4.18 (0.22)	36.22 (2.32)

Table 6 .

 6 Inventory costs -assembly system simulation cases 1-2.

		Buffer 1 Buffer 2 Buffer 3 Buffer 4 Buffer 5 Buffer 6 Buffer 7
	Inventory cost	0.5	0.5	0.5	0.5	1.0	1.0	2.0

Table 7 .

 7 Best control parameters -Assembly system simulation case 1. (demand arrival rate = 0.6). B 1 / K 1 B 2 / K 2 B 3 /K 3 B 4 / K 4 B 5 / K 5 B 6 /K 6 B 7 /K 7 2

	Policy
									-
	Base stock	0 ∞	0 ∞	0 ∞	0 ∞	2 ∞	2 ∞	5 ∞	-
	Conwip	2	2	2	2	0	0	2	4
		∞	∞	∞	∞	∞	∞	∞
	Hybrid	4	4	4	4	1	1	0	5
		4	4	4	4	1	1	∞
	Extended Kanban	0 2	0 2	0 2	0 2	0 4	0 4	4 4	-

Table 8 .

 8 Average inventories -backorders average -expected profit -assembly system simulation case 1. (demand arrival rate = 0.6). The cells contain the mean of the performance measures and the standard deviation in parenthesis.

	Policy	Buffer1	Buffer 2	Buffer 3	Buffer 4	Buffer 5	Buffer 6	Buffer 7	Backorders average	[F E	]
	Kanban	0.39 (0.0)	0.39 (0.0)	0.39 (0.0)	0.39 (0.0)	0.38 (0.0)	0.38 (0.0)	1.0 (0.0)	0.94 (0.06)	29.0 (0.6)
	Base Stock	0.24 (0.03)	0.24 (0.03)	0.24 (0.03)	0.24 (0.02)	1.19 (0.04)	1.18 (0.04)	3.64 (0.02)	0.45 (0.03)	27.38 (0.42)
	CONWIP	0.06 (0.0)	0.06 (0.0)	0.07 (0.0)	0.06 (0.0)	0.12 (0.0)	0.12 (0.01)	1.71 (0.01)	0.86 (0.06)	29.56 (0.6)
	Hybrid	0.24 (0.01)	0.24 (0.01)	0.24 (0.01)	0.24 (0.01)	0.03 (0.0)	0.03 (0.0)	2.53 (0.01)	0.69 (0.04)	29.48 (0.38)
	Extended Kanban	0.11 (0.0)	0.1 (0.0)	0.1 (0.0)	0.1 (0.0)	0.19 (0.03)	0.19 (0.04)	1.71 (0.02)	0.81 (0.06)	29.74 (0.51)

Table 9 .

 9 Best control parameters -Assembly system simulation case 2. (demand arrival rate = 0.8).

	Policy	B 1 / K 1 B 2 / K 2 B 3 /K 3 B 4 / K 4 B 5 / K 5 B 6 /K 6	B 7 /K 7	C
	Kanban	2 2	2 2	2 2	2 2	4 4	4 4	9 9	-
	Base stock	1 ∞	1 ∞	1 ∞	1 ∞	6 ∞	6 ∞	15 ∞	-
	Conwip	4	4	4	4	1	1	8	13
		∞	∞	∞	∞	∞	∞	∞
	Hybrid	5	5	5	5	8	8	0	13
		5	5	5	5	8	8	∞
	Extended Kanban	0 5	0 5	0 5	0 5	1 9	1 9	10 11	-

Table 10 .

 10 Average inventories -backorders average -expected profit -assembly system simulation case 2. (demand arrival rate = 0.8). The cells contain the mean of the performance measures and the standard deviation in parenthesis.

	Policy	Buffer1	Buffer 2	Buffer 3	Buffer 4	Buffer 5	Buffer 6	Buffer 7	Backorders average	[F E	]
	Kanban	1.13 (0.0)	1.13 (0.0)	1.13 (0.0)	1.13 (0.0)	2.81 (0.02)	2.81 (0.02)	5.5 (0.06)	1.8 (0.2)	19.12 (1.92)
	Base Stock	1.11 (0.1)	1.12 (0.1)	1.1 (0.09)	1.1 (0.08	5.81 (0.1)	5.83 (0.1)	11.2 (0.06)	0.59 (0.06)	13.82 (0.76)
	CONWIP 0.64 (0.04)	0.63 (0.04)	0.63 (0.05)	0.63 (0.05)	1.14 (0.06)	1.16 (0.06)	7.2 (0.08)	1.56 (0.14)	22.42 (1.42)
	Hybrid	0.59 (0.03)	0.59 (0.02)	0.6 (0.02)	0.6 (0.03)	0.93 (0.05)	0.94 (0.05)	7.16 (0.07)	1.61 (0.19)	22.52 (1.9)
	Extended Kanban	0.61 (0.04)	0.61 (0.04)	0.61 (0.02)	0.61 (0.04)	1.08 (0.06)	1.08 (0.05)	5.46 (0.06)	1.86 (0.28)	22.86 (1.89)

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