OPTIMIZATION OF CONTINUOUSFLOW TRANSFER LINES WITH DELAY USING IPA Iyad Mourania, Sophie Hennequina, Xiaolan Xiea, b a
INRIA/MACSI team and LGIPM, ISGMPBat. A. Ile du Saulcy, 57045 Metz Cedex 1 France {mourani, hennequin}@loria.fr b Ecole des Mines de SaintEtienne 158 cours Fauriel, 42023 SaintEtienne cedex 2 France
[email protected]
Abstract: This paper addresses the optimization of failureprone transfer lines with delays and constant demand. For this purpose, we propose a new continuous flow model. Machines are subject to timedependent 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 simulationbased 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: continuousflow model, transfer lines, transfer delays, simulationbased 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 leadtimes 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 workinprocess. 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 continuousflow singlestage singleproduct 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 ndimensional 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 singleproduct 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 operationdependant failures. Xie (2002a) addressed the performance evaluation and optimization of failureprone discreteevent system by using a fluidstochasticevent graph model, which is a decisionfree 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 continuousflow 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 timedependent 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
MK1
τK1
BK1
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 Mi1 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 Mi1 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 intransit 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). : interarrival 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 Bi1 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 timesdependant 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. SIMULATIONBASED 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 kth 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 23. 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 gradientbased optimization algorithm → to minimize J N (θ ) . It is expected that this optimal →
solution converges to true optimal solution of J (θ ) as N increases. The gradientbased optimization procedure is as follows: → ⎛ ⎞ ∂ hin +1 = ⎜ hin + sn J N (θ n ) ⎟ ∂ h i ⎝ ⎠
+
6. NUMERICAL RESULTS In this section, we consider 3machine 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 stepsize at iteration n. c+ 5
Two stepsizes are used. At the beginning, we use the Armijo stepsize (Armijo, 1966; Polak, 1997), i.e. the maximum stepsize sn = γ m such that →
→
J N (θ n +1 ) − J N (θ n ) ≤ −α γ
m
MTBFi MTTRi 100 20
Ui 4
D 0.5
The simulationbased optimization is performed for N = 1000 repairs of machine M3.
→
 ∇J N (θ n ) 2
with α = 0.5 and γ = 0.5. Armijo stepsize (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 stepsize (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 stepsize. It reduces it if no → improvement of J N (θ ) is observed after a certain number of iterations. SimulationOptimization 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 24. Step 3. Update h using subgradient method with either Armijo stepsize or standard subgradient stepsize. Step 4. Repeat Steps 23 till convergence.
Fig. 4. Cost function versus the buffer levels for 3machine 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 simulationbased 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 simulationbased 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 stepsize)
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 stepsize)
2
7. CONCLUSIONS In this paper, we have considered continuousflow 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 timedependent 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 simulationbased optimization algorithm with samplepath 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 operationdependant failures. REFERENCES Armijo, L. (1966), “Minimization of functions having Lipschitz continuous firstpartial derivatives”, Pacific Journal of Mathematics, vol. 16, pp. 13. Fu, M. and X. Xie (2002), “Derivative estimation for buffer capacity of continuous transfer lines subject to operationdependent failures”, Discrete Event Dynamic Systems: Theory and Applications, vol. 12, pp. 447469.
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. 557580. 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. 37323737. Mourani, I., S. Hennequin and X. Xie (2005b), “Simulationbased Optimization of a SingleStage FailureProne 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. 109138. Polak, E. (1997), Optimization algorithms and consistent approximations, SpringerVerlag, 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. 14071422. Xie, X. (2002a), “Fluidstochasticevent graphs for evaluation and optimization of discreteevent system with failures”, IEEE Transactions on Robotics and Automation, vol. 18, no. 3, pp. 360367. Xie, X. (2002b), “Evaluation and optimization of twostage continuous transfer lines subject to timedependent failures”, Discrete Event Dynamic Systems: Theory and Applications, vol. 12, pp. 109122.