OPTIMIZATION OF CONTINUOUS-FLOW TRANSFER LINES WITH DELAY USING IPA

OPTIMIZATION OF CONTINUOUS-FLOW TRANSFER LINES WITH DELAY USING IPA

OPTIMIZATION OF CONTINUOUS-FLOW TRANSFER LINES WITH DELAY USING IPA Iyad Mourania, Sophie Hennequina, Xiaolan Xiea, b a INRIA/MACSI team and LGIPM, I...

437KB Sizes 0 Downloads 13 Views

OPTIMIZATION OF CONTINUOUS-FLOW TRANSFER LINES WITH DELAY USING IPA Iyad Mourania, Sophie Hennequina, Xiaolan Xiea, b a

INRIA/MACSI team and LGIPM, ISGMP-Bat. A. Ile du Saulcy, 57045 Metz Cedex 1 France {mourani, hennequin}@loria.fr b Ecole des Mines de Saint-Etienne 158 cours Fauriel, 42023 Saint-Etienne cedex 2 France [email protected]

Abstract: This paper addresses the optimization of failure-prone transfer lines with delays and constant demand. For this purpose, we propose a new continuous flow model. Machines are subject to time-dependent failures, times to failure and times to repair are random variable with general distribution. Material flowing out a machine waits a period of time called delay for material transfer before arriving at its downstream buffer. We present a simulation-based optimization method for determining optimal buffer levels in order to minimize the long run average cost. The optimization algorithm is based on the Infinitesimal Perturbation Analysis for estimation of gradients along the simulation. Copyright © 2006 IFAC Keywords: continuous-flow model, transfer lines, transfer delays, simulation-based optimization, perturbation analysis.

1. INTRODUCTION Continuous flow models have been widely used for optimal control and design of manufacturing system. Techniques for optimization of system parameters are simpler than those for discrete parameter optimization. Although continuous flow models offer an interesting way to reduce the complexity inherent in traditional discrete parts modeling, existing continuous flow models neglect some important characteristics of manufacturing systems such as production lead-times and transportation delays. In existing continuous flow models, material flows continuously and instantaneously from machines to machines. However, many manufacturing processes have significant delays in the material flow, such delays occur in oven processes (e.g. semiconductor diffusion), drying processes and testing. These delays usually have great impact on performance measures such as customer response time and work-in-process. Unfortunately, most existing continuous flow models do not take into account these delays. To our knowledge, there are only three exceptions. i) Van

Ryzin, et al. (1991) explicitly considered the impact of delays for optimal flow control of job shops composed of a set of flexible unreliable machines in order to minimize the discount and long run average cost. A heuristic control policy for a flow shop with delay is derived using theoretical arguments and approximations. ii) Mourani, et al. (2005a) extended the model of Van Ryzin, et al. (1991) and proposed a continuous Petri net model with delays for performance modeling and optimization of transfer lines. The new model is more realistic than classical continuous flow models yet keeps the simplicity and analyticity of continuous flow models. iii) Mourani, et al. (2005b) addressed the optimization of a continuous-flow single-stage single-product manufacturing system with delay and constant demand. They found by using Infinitesimal Perturbation Analysis (IPA) technique, the optimal buffer level, which minimizes the long run average cost. The IPA techniques have been widely considered for control and optimization of discrete event systems. In the pioneer work, motivated by buffer storage

optimization problem in a production line, Ho, et al. (1979) developed an efficient technique called perturbation analysis (PA). It enables one to compute the sensitivity of a performance measure with respect to the system parameters by a single simulation run. Ho and Cao (1991) developed an IPA technique for the efficient computation of n-dimensional gradient vector of performance measure, J(θ), of a Discrete Event Dynamic System with respect to its parameter θ (such as buffer size, inflow rate, service rate) using only one statistical experiment of the system. This is opposed to the traditional method of estimating sensitivity information such as dJ(θ)/dθ. IPA calculates directly the sample derivative dL(θ, ξ)/dθ using information on the nominal trajectory (θ, ξ) alone. The basic idea is the following: if the perturbation introduced into the trajectory (θ, ξ) are sufficient small, then we can assume that the event sequence of the perturbed trajectory (θ + ∆θ, ξ) remains unchanged from the nominal, i.e., the two trajectories are deterministically similar in the order of their event sequence. In this case the derivative dL(θ, ξ) /dθ can be calculated easily, see also (Glasserman, 1990).

2. CONTINUOUS FLOW MODEL WITH DELAY This paper considers continuous flow models of single-product transfer line composed of K machines (M1, M2, .., MK) and K buffers (B1, B2,…, BK) of finite capacity. Material flows continuously from outside the system to the first machine, then waits a period of time called delay τ1 for material transfer, before arriving to the first buffer B1, then the second machine M2, then waits a second period of time called delay τ2 before arriving to the second buffer B2, and so forth to the last machine, where it leaves the system (see Fig.1). We assume that τi > 0, ∀ 0 < i < K. Hence a delay is considered between a machine Mi and its downstream buffer Bi. That means that the parts produced on Mi would not arrive immediately in the buffer Bi. So, a delay will occur in the delivering of the parts to the downstream machine. These delays are introduced to explicitly account for processing times and transportation times. D

Fu and Xie (2002), estimated the derivatives of the throughput rate with respect to buffer capacity for continuous flow models of a transfer line comprising two machines separated by a buffer of finite capacity and subject to operation-dependant failures. Xie (2002a) addressed the performance evaluation and optimization of failure-prone discrete-event system by using a fluid-stochastic-event graph model, which is a decision-free Petri net. In addition, Xie (2002b) considered continuous flow transfer lines composed of two machines subject to timedependent failures and separated by a buffer of finite capacity. He established a set of evolution equations that determines the continuous state variables, i.e., cumulative productions and the buffer level, at epochs of discrete events. Based on these evolution equations, he proved the concavity of the throughput rates of the machines and derived gradient estimators and proposed a single sample path optimization algorithm. In this paper we consider a continuous-flow model of transfer line composed of K machines and K buffers of finite capacity. Material transfer between a machine and its downstream buffer takes a finite amount of time called delay. Machines are subject to time-dependent failures. Time to failure and times to repair are generally distributed. The goal of this paper is to find, by using IPA technique, the vector of optimal buffer level, which minimizes the long run average cost. The rest of the paper is organized as follows. Section 2 presents the continuous flow model with delays of transfer lines. Section 3 presents the IPA technique for gradient estimation. The cost function and its derivatives are addressed in Section 4. A simulationbased optimization algorithm is presented in Section 5. Numerical results are given in Section 6. We conclude the paper in Section 7.

M1

τ1

B1

M2

τ2

MK-1

τK-1

BK-1

MK

τK

BK

Fig. 1. Production line with a delay between Mi and Bi

We assume that the input buffer of the first machine is never empty (first machine is never starved) and the output buffer of the last machine is never full (the last machine is never blocked). The following notations at time t will be used throughout the paper (with i = 1, 2,…, K): αi(t)

: the failure state of machine Mi with αi(t) = 1 if Mi is up and αi(t) = 0 otherwise. Let α(t) = [α1(t), α2(t),…., αK(t)]. : the remaining life time (until failure of ri(t) machine Mi if αi(t) = 1 or repair of machine Mi if αi(t) = 0)). : the age (since last failure or repair) of ai(t) machine Mi in state αi(t). Let a(t) = [a1(t), a2(t),……, aK(t)]. : the buffer capacity of the buffer Bi. It is hi extended to include both materials in buffer Bi and materials that have been produced by machine Mi-1 and are in transit to buffer Bi. Physical meaning of this buffer capacity is either the number of kanbans in a JIT system or the physical buffer capacity plus the capacity of convey connecting machine Mi-1 and buffer Bi. xi(t) ∈ : the quantity of material in the buffer Bi. Let x(t) = [x1(t), x2(t), …, xK(t)]. [0, hi] yi(t) : the inventory position of the buffer Bi, i.e. the sum of the inventory level xi(t) and the parts in-transit to Bi (see Fig 2). Let y(t) = [y1(t), y2(t), …, yK(t)]. : buffer capacity of the buffer Bi, i.e. the hi upper limit of yi(t). : the kth event. ek : the epoch of event ek with t0 = 0. tk ui(t) ∈ : the production rate of the machine Mi with Ui as the maximal production rate of [0, Ui]

udi(t)

∆k TBFi,k

TTRi,k D

Mi. : the rate of the flow entering buffer Bi taking into account the delay τi for material transfer, with udi(t) = ui(t - τi) (see Fig 2). : inter-arrival times of events (time between ek and ek+1 i.e., ∆k = tk+1 - tk). : kth time to failure of machine Mi. For each i, TBFi,k are independent and identically distributed (i.i.d.) random variables of general distribution. : kth time to repair of machine Mi. For each i, TTRi,k are i.i.d. random variables of general distribution. : constant demand rate.

at its upstream buffer Bi-1 and availability of free kanbans Hi – yi of stage i. The dynamics of the system are given by the following equations: • The state S(t) of the system at time t, ∀ t ∈[tk, tk+1) yi(t) = yi(tk) + [(ui(tk) - ui+1(tk)) × (t - tk)], xi(t) = xi(tk) + [(ui(tk - τi) - ui+1(tk)) × (t - tk)], ri(t) = ri(tk) - (t - tk). •

⎧ri (tk ), ⎪ xi (tk ) ⎪ , ⎪⎪ ui +1 (tk ) − ui (tk − τ i ) ∆ k = tk +1 − tk = ⎨ ⎪ hi − yi (tk ) , ⎪ ui (tk ) − ui +1 (tk ) ⎪ ⎪⎩t Ni ( tk −τ i ) + τ i − tk

Machines are subject to times-dependant failures. This means that a machine can fail even it is not working on a part. Each machine can be either operational (αi(t) = 1) or down (αi(t) = 0). When it is operational, it can be either working or idle (if it is starved or blocked). Consider that an event ek causes a change in the production rate ui of the machine Mi at time tk. The effect of this change on the downstream buffer Bi will arrive after τi time units (see Fig. 2). Mi ui(t)

τi

udi= ui(t-τi)

Bi

Mi+1 xi(t)

yi(t)

ui+1(t)

Fig. 2. Delay between Mi and Bi

Different events are possible: the failure of a machine Mi, the repair of Mi, buffer full of Bi (i.e. yi(t) reaches Hi), buffer empty of Bi (i.e. xi(t) becomes 0) denoted respectively, Fi, Ri, BFi, BEi. Furthermore, any change in the production rate ui (with i = 1, 2, …., K) of a machine Mi is coupled with an event called (DLi). Thus if there is a change in a production rate ui at time tk, so udi will change to the same rate at time (tk + τi). So ek ∈{ Fi, Ri, BFi, BEi, DLi}. For the implementation purpose of the simulation, a FIFO queue is used to record all active events (DLi). Whenever the production rate of a machine change at time t, we add into the queue a couple containing the new value of ui and t + τi, i.e. the corresponding time of the change of udi. The control policy is defined as follows: if ⎧0, ⎪U , if ⎪⎪ i u i (t ) = ⎨min(ud i −1 (t ),U i ), if ⎪min(u (t ),U ), if i +1 i ⎪ ⎪⎩min(u i +1 (t ), ud i −1 (t ),U i ), if

α i (t ) = 0, xi −1 > 0, y i < hi , α i (t ) = 1, xi −1 = 0, y i < hi , α i (t ) = 1, xi −1 > 0, y i = hi , α i (t ) = 1, xi −1 = 0, y i = hi , α i (t ) = 1.

where ud0(t) = ∞, uK+1(t) = D. In words, the control policy is a kanban policy which ensures that the material in each stage including material in buffer Bi and material in transit does not exceed a given number of kanbans. The production rate of each machine Mi is also limited by availability of material

Next event epoch ∆k = tk+1 - tk can be determined as follows: if ek +1 = Ri or Fi if ek +1 = BEi if ek +1 = BFi if ek +1 = DLi

where Ni(t) is the index k of the event ek corresponding to the last change of ui prior to time t. •

Next state S(tk+1):

⎧TBFi or TTRi , if ek +1 = Ri or Fi ri (tk +1 ) = ⎨ otherwise ⎩ri (tk ) − ∆ k ⎧0 , xi (tk +1 ) = ⎨ ⎩ xi (tk ) + ((ui (tk − τi ) − ui +1 (tk )) × ∆ k )

⎧hi , yi (tk +1 ) = ⎨ ⎩ yi (tk ) + ((ui (tk ) − ui +1 (tk )) × ∆ k )

if ek +1 = BEi otherwise

if ek +1 = BFi otherwise

tk+1 = tk + ∆k 3. IPA For the purpose of IPA technique, we derive firstly the next event epoch ∆k, and then the next state S(tk+1) with respect to the vector of optimal buffer level which is given by: ⎡h1 ⎤ ⎢h ⎥ ⎢ 2⎥. → → θ = h = ⎢. ⎥ ⎢ ⎥ ⎢. ⎥ ⎢hn ⎥ ⎣ ⎦

The initial conditions proposed in our case are: xi(0) = yi(0) = hi/2 and αi(0) = 1. • The derivatives of the next event epoch ∆k with →



respect to θ = h is given as follows: ⎧ ∂ri (tk ) if ek +1 = Ri or Fi ⎪ → , ⎪ ∂θ ⎪ ∂xi (tk ) 1 , if ek +1 = BEi ⎪ → τ u ( t ) − u ( t − ) i k i ∂∆ k ⎪ i +1 k ∂θ =⎨ → ∂h ∂y (t ) 1 ∂θ ⎪ ( i − i →k ), if ek +1 = BFi ⎪ ui (tk ) − ui +1 (tk ) ∂ θ→ ∂θ ⎪ ⎪ ∂t Ni (tk −τ i ) ∂τ i ∂tk + →− → if ek +1 = DLi → ⎪ ⎪⎩ ∂ θ ∂θ ∂θ



The derivatives of the next state S(tk+1) with →



respect to θ = h is given as follows: ∂ri (tk +1 ) →

∂θ

if ek +1 = Ri or Fi

⎧0, ⎪ = ⎨ ∂ri (tk ) ∂∆ k − → ⎪ → ∂θ ⎩ ∂θ



∂θ



∂J N (θ )

otherwise



∂θ

if ek +1 = BEi ⎧0, ⎪ ∂∆ = ⎨ ∂xi (tk ) + (u i (tk − τ i ) − ui +1 (tk )) →k ⎪ → ∂θ ⎩ ∂θ

∂xi (tk +1 )

otherwise

⎧ ∂hi if ek +1 = BFi ⎪ →, ∂yi (tk +1 ) ⎪ ∂ θ =⎨ → ⎪ ∂yi (tk ) + (u (t ) − u (t )) ∂∆ k ∂θ otherwise i k i +1 k → ⎪ → ∂θ ⎩ ∂θ

∂t k +1 →

∂θ



=

∂t k →

∂θ

+

The derivative of the cost function with respect to the vector of optimal buffer level can be written as follows:





∂xi (t0 ) →

∂θ

∂yi (t0 ) →

∂θ

∂ t0 →

∂θ

∑∑ k



∂θ

i







Gik( θ ) and its derivative ∂Gik( θ )/∂ θ . Thus we will discuss all the possible cases needed for exact computing of the cost function and its derivative. In the following and for the sake of simplicity xi,k, xi,k+1 are used to denote respectively xi(tk) and xi(tk+1) used before.

∂θ







∂Gik (θ )

Case 1: xi,k ≥ 0, xi,k+1 ≥ 0

∂∆ k .

Gik (θ ) = (



respect to θ = h is given as follows: ∂θ

∂t N / ∂ θ → 1 J (θ ) + tN tN

For calculating this derivative, we need to evaluate the cost function between each two successive events

The derivatives of the state S(t0) at time t0 with ∂ri (t 0 )



=−



∂Gik (θ ) →

∂θ

=0

) × ∆k × c+

∂∆ c + ∂xi ,k +1 ∂xi ,k (( → + → )∆ k + →k ( xi ,k +1 + xi ,k )) 2 ∂θ ∂θ ∂θ

=

xi ,k +1 + xi ,k



Gik (θ ) = −(

if i = j if i ≠ j

⎧0.5, =⎨ ⎩0,

2

Case 2: xi,k ≤ 0, xi,k+1 ≤ 0

if i = j if i ≠ j

⎧0.5, =⎨ ⎩0,

xi ,k +1 + xi ,k



∂Gik (θ ) →

∂θ

=0.

=−

2

) × ∆k × c−

∂∆ c − ∂xi ,k +1 ∂xi ,k (( → + → ) ∆ k + →k ( xi ,k +1 + xi ,k )) 2 ∂θ ∂θ ∂θ

Case 3: xi,k ≥ 0, xi,k+1 ≤ 0

4. COST FUNCTION AND ITS DERIVATIVES

This case is illustrated in Fig. 3. The inventory trajectory crosses 0 at time t*, i.e. xi(t*) = 0 with:

xi(t)

The average cost function is the long run average inventory holding and backlogging cost defined as the limiting function of: →

J N (θ ) =

1 tN



0

Note that cost structure is very common in inventory control literature except that the inventory holding cost is assumed to be identical for all stages. Results of this paper can be easily extended to take into account stage specific inventory holding cost and cost related to inventory in transit.

1 tN

N −1

∑∑ G k =0

ik



(θ )

i



tk +1

∫ g ( x (s)) ds . i

tk

t * − tk =

i

tk+1

t

xi ,k u i +1,k − ud i ,k

Thus, →

Gik (θ ) = →

∂Gik (θ ) →

∂θ

=

+ 2 − 2 1 (c xi ,k + c xi ,k +1 ) 2 u i +1,k − ud i , k

∂xi ,k ∂x 1 ( c + × xi , k × → + c − × xi ,k +1 × i ,→k +1 ) ui +1,k − ud i ,k ∂θ ∂θ

Case 4: xk ≤ 0, xk+1 ≥ 0 →



∂Gik (θ ) →

∂θ

with

c

-

Fig. 3. Cost function for the case where xi,k ≥ 0, xi,k+1 ≤ 0

Gik (θ ) = −

It can be rewritten as follows:

t*

xi(tk+1)

i

where c+, c- denote the unit inventory holding cost and backlogging cost with c+ > 0, c- > 0.

Gik (θ ) =

tk

∫ ∑ gi ( xi (s)) ds

+ ⎪⎧c x if x ≥ 0 g i ( x) = ⎨ − ⎪⎩−c x if x < 0



c+

tN

where N is the constant number of events, gi(xi(s)) corresponds to the inventory cost which is given by:

J N (θ ) =

xi(tk)

J (θ )

=

− 2 + 2 − 2 + 2 1 (−c xi ,k + c xi ,k +1 ) 1 (c xi ,k + c xi ,k +1 ) = 2 ui+1,k − udi ,k 2 udi ,k − ui+1,k

∂xi ,k ∂x 1 (c − × xi , k × → + c + × xi , k +1 × i ,→k +1 ) ud i ,k − u i +1,k ∂θ ∂θ

5. SIMULATION-BASED OPTIMIZATION ALGORITHM In this section, we present an algorithm for → minimizing the average cost J (θ ) . The basic idea is to approximate the optimization of long run average cost by the optimization of a sample path function. In this paper, we exploit the following property. If (i) the total number of events N is defined as a function of failure/repair events such as k-th repair of machine MK and if (ii) common times to failure

The simulation of the continuous flow model is a discrete event algorithm as described in Section 2 and progresses from event to event. At each event, state variables, clock variables and their gradients are updated according to equations of Sections 2-3. n → Cumulative cost G (θ ) and its gradient

∑∑ k =0



n

∑∑ k =0

i

∂Gik (θ )

ik

i

are also computed according to



∂θ

equations of Section 4. At the end of the simulation, →





JN( θ ) and ∂JN( θ )/∂ θ are determined using results of Section 4.



(TBF) and times to repair (TTR) are used for all θ , → then the sample function J N (θ ) is a continuous function. We then use a gradient-based optimization algorithm → to minimize J N (θ ) . It is expected that this optimal →

solution converges to true optimal solution of J (θ ) as N increases. The gradient-based optimization procedure is as follows: → ⎛ ⎞ ∂ hin +1 = ⎜ hin + sn J N (θ n ) ⎟ ∂ h i ⎝ ⎠

+

6. NUMERICAL RESULTS In this section, we consider 3-machine production lines with 3 identical machines except the delays with τ1 = 2.0 time units, τ2= 2.2 time units, τ3= 2.5 time units. Times between failure and times to repair are exponentially distributed with rate λi and µi respectively, i.e., mean time between failure MTBFi = 1/λi and mean time to repair MTTRi = 1/µi. Initial conditions are xi(0) = yi(0) = hi/2 and

αi(0) = 1. The parameters are summarized in

Table 1.

Table 1 Simulation data

where sn is the step-size at iteration n. c+ 5

Two step-sizes are used. At the beginning, we use the Armijo step-size (Armijo, 1966; Polak, 1997), i.e. the maximum step-size sn = γ m such that →



J N (θ n +1 ) − J N (θ n ) ≤ −α γ

m

MTBFi MTTRi 100 20

Ui 4

D 0.5

The simulation-based optimization is performed for N = 1000 repairs of machine M3.



|| ∇J N (θ n ) ||2

with α = 0.5 and γ = 0.5. Armijo step-size (Armijo, 1966; Polak, 1997) allows fast descent but its determination is time consuming. Consequently, → n +1

c250

→ n

when the cost improvement J N (θ ) − J N (θ ) is smaller than a given percentage, we switch to standard gradient optimization procedure with the following step-size (Nedić and Bertsekas, 2002):

We

first

run

the

simulation

for

evaluating



exhaustively the cost function J (θ ) by varying h1, h2 and h3 along integer points. The minimal cost is obtained with h1 = 7. The resulting cost function for h1 = 7 is plotted in Fig. 4.



sn = η

J N (θ n ) →

|| ∇J N (θ n ) ||2

At the beginning, η is chosen such that sn is equal to the last Armijo step-size. It reduces it if no → improvement of J N (θ ) is observed after a certain number of iterations. Simulation-Optimization Algorithm: Step 1. Initialization: h ← h0 Step 2. Simulate the continuous flow model with parameter vector h and with common random variable stream till a given event fulfilling above →





condition (i) to determine JN( θ ) and ∂JN( θ )/∂ θ using results of Sections 2-4. Step 3. Update h using sub-gradient method with either Armijo step-size or standard sub-gradient stepsize. Step 4. Repeat Steps 2-3 till convergence.

Fig. 4. Cost function versus the buffer levels for 3-machine line

The cost function is convex and there exist a minimum cost value equal to (306.34 monetary unit) which corresponds to the following vector of buffer level: (h1 = 7 parts, h2 = 14 parts, h3 = 32 parts).

The simulation-based optimization algorithm is then run to optimize the transfer line model. The cost function converges rapidly towards the optimal value of Fig. 4 (306.23 monetary unit), and the vector of optimal buffer level corresponding to that optimal cost value is: (h1 = 6.49 parts, h2 = 15.63 parts, h3 = 32.36 parts) as given in Table 2. Table 2 Results of simulation-based optimization N° of iteration

Cost value

Vector of buffer level (h1(parts), h2(parts), h3(parts))

(monetary unit)

0 1

387.88

(20.00, 20.00, 20.00)

327.28

(11.69, 18.99, 31.36)

(with Armijo step-size)

309.58

(7.61, 16.98, 30.92)

3 4 5 6 7 8 9 10 11

307.84 307.23 306.63 306.41 306.31 306.27 306.25 306.24 306.23

(5.77, 16.01, 31.51) (7.07, 15.98, 32.11) (6.74, 15.78, 32.22) (6.61, 15.7, 32.29) (6.55, 15.66, 32.32) (6.52, 15.64, 32.34) (6.51, 15.64, 32.33) (6.50, 15.63, 32.35) (6.49, 15.63, 32.36)

(with Armijo step-size)

2

7. CONCLUSIONS In this paper, we have considered continuous-flow transfer lines composed of K machines and K buffers with delays (such as transportation delays) and constant demand. Material waits a period of time called delay for material transfer before arriving from machine to its downstream buffer. Machines are subject to time-dependent failures i.e., the machine can fail even it is not working on a part. Times to failure and times to repair are random variable with general distribution. A simulation-based optimization algorithm with sample-path gradients has been developed to determine the vector of optimal buffer level, which minimizes the long run average cost. Future research is necessary to investigate the convexity of the cost function obtained by simulation. In addition, it would be interesting to study the case of transfer line where the machines are subject to operation-dependant failures. REFERENCES Armijo, L. (1966), “Minimization of functions having Lipschitz continuous first-partial derivatives”, Pacific Journal of Mathematics, vol. 16, pp. 1-3. Fu, M. and X. Xie (2002), “Derivative estimation for buffer capacity of continuous transfer lines subject to operation-dependent failures”, Discrete Event Dynamic Systems: Theory and Applications, vol. 12, pp. 447-469.

Glasserman, P. (1990), Gradient Estimation via Perturbation Analysis, Kluwer Academic Publisher. Ho, Y.C., A. Eyler and T.T. Chien (1979), “A gradient technique for general buffer storage design in a serial production line”, International Journal of Production Research, vol. 17, pp. 557-580. Ho, Y. and X.R. Cao (1991), Perturbation analysis of discrete event dynamic systems. Boston, MA: Kluwer Academic Publishers. Mourani, I., S. Hennequin and X. Xie (2005a), “Continuous Petri nets with delays for performance evaluation of transfer lines”, Proceedings of IEEE International Conference on Robotics and Automation (ICRA2005), Barcelona, Spain, pp. 3732-3737. Mourani, I., S. Hennequin and X. Xie (2005b), “Simulation-based Optimization of a SingleStage Failure-Prone Manufacturing System with Transportation Delay” Proceedings of International Conference on Industrial Engineering and Systems Management(IESM), Marrakech, Morocco. Nedić, A. and D.P. Bertsekas (2002) “Incremental subgradient methods for nondifferentiable optimization”, SIAM on Optimization, vol. 12, No. 1, pp. 109-138. Polak, E. (1997), Optimization algorithms and consistent approximations, Springer-Verlag, New York. Van Ryzin, G.J., S.X.C. Lou and S.B. Gershwin (1991), “Scheduling job shops with delay”, Int. J. Prod. Res., vol. 29, no. 7, pp. 1407-1422. Xie, X. (2002a), “Fluid-stochastic-event graphs for evaluation and optimization of discrete-event system with failures”, IEEE Transactions on Robotics and Automation, vol. 18, no. 3, pp. 360-367. Xie, X. (2002b), “Evaluation and optimization of two-stage continuous transfer lines subject to time-dependent failures”, Discrete Event Dynamic Systems: Theory and Applications, vol. 12, pp. 109-122.