Samuel Vercraene
email: samuel.vercraene@grenoble-inp.fr

Jean-Phillipe Gayon
email: jean-philippe.gayon@grenoble-inp.fr

A two-stage production/inventory system with admission control of product returns ⋆

Keywords: Multi-echelon, Reverse logistics, Production/Inventory system, Optimal Control, Remanufacturing

We consider a two-stage production-inventory system. Each stage is made of an exponential production server and an output inventory. The downstream inventory faces a Poisson demand and each inventory receives product returns that can be either accepted or rejected. The objective is to find a production/admission policy minimizing the expected discounted costs over an infinite horizon. We consider linear holding, backorder, production, admission and rejecting costs. We show that the optimal policy is characterized by state-dependent base-stock levels and state-dependent disposal thresholds. We also obtain monotonicity results for these switching curves. Interestingly, our results can be transposed to an hybrid singleechelon system with a production facility and a remanufacturing facility, both with exponential leadtimes. Finally, we design heuristics to control production and disposal and we carry out a numerical study to investigate their performances.

Introduction

Products are more and more returned in supply chains. Customers can return products a short time after purchase due to take-back commitments of the supplier. The proportion of returns is particularly important in electronic business where customers can not touch a product before purchasing it. Customers might also return used products a long time after purchase. This type of return has increased in recent years due to new regulations encouraging waste reduction, especially in Europe. Some industries also encourage it for economical and marketing reasons.

Returns constitute a reverse flow, from the customer to the supplier, which complicates inventory control. A recent review of literature on reverse logistics is provided by Ilgin and Gupta [START_REF] Ilgin | Environmentally conscious manufacturing and product recovery (ecmpro): A review of the state of the art[END_REF]. Inventory management constitutes a part of this review, and we can see that an abundant literature is devoted to the inventory control of single echelon systems. In several situations, the structure of the optimal policy is similar to the case without returns. For instance, Fleischmann and Kuik [START_REF] Fleischmann | On optimal inventory control with independent stochastic item returns[END_REF] consider a single inventory with stochastic demand and stochastic independent returns. To model the returns, they consider a demand that can be both positive or negative. They show average cost optimality of an (s, S) policy.

In this paper, we consider a two-stage production /inventory system with returns of products at each stage and a disposal option upon arrival of returns. Before presenting in detail our model, we briefly review the related literature of inventory control. In particular, we focus on models with a disposal option and on multi-echelon models.

Disposal option

We first review the inventory control literature considering a disposal option for returns. Heyman [START_REF] Heyman | Optimal disposal policies for a single-item inventory system with returns[END_REF] considers an inventory system with backorders, Poisson demand and Poisson returns. He assumes zero lead-times and linear costs for both procurement and remanufacturing. These strong assumptions imply that the optimal production policy is a make-to-order policy. The optimal disposal policy is a simple threshold policy: when the inventory level exceeds a certain disposal threshold R, every returned item is disposed upon arrival. He also provides an explicit expression for the optimal disposal threshold. Zerhouni et al. [START_REF] Zerhouni | Admission control of product returns in a reverse logistic context[END_REF] relax the assumption of zero manufacturing lead-time. Only one item at a time can be manufactured and the manufacturing leadtime is exponentially distributed. For both lost sale and backorder cases, the optimal policy is an (R, S) policy stating to produce when the net inventory is below S and to dispose returned items when the net inventory exceeds R. The ordering of the two optimal parameters depends on the ordering of the linear manufacturing and remanufacturing costs. In the backorder case, they obtain an explicit expression for the optimal base-stock level when the difference (S -R) is given. The disposal policy of [START_REF] Heyman | Optimal disposal policies for a single-item inventory system with returns[END_REF] can be seen as a special case of the (R, S) policy where S = 0.

Ouyang and Zhu [START_REF] Ouyang | An inventory control model with disposal for manufacturing/remanufacturing hybrid system[END_REF] relax the zero lead-time and fixed cost assumptions. The production lead-time is deterministic while the remanufacturing leadtime is neglected. They consider an (s p , Q p , s d ) where s p is the reorder point, Q p the order quantity and s d the disposal threshold (returns are remanufactured if and only if the stock level is below s d ). This policy was first introduced in the paper of van der Laan et al. [START_REF] Van Der Laan | Product remanufacturing and disposal: A numerical comparison of alternative control strategies[END_REF] which is presented in the next paragraph. Ouyang and Zhu derive analytical expressions for the average cost while van der Laan et al. use numerical methods. The optimal parameters can then be obtained by solving a difficult nonlinear integer programming. The (R, S) policy in [START_REF] Zerhouni | Admission control of product returns in a reverse logistic context[END_REF] can be seen as an (s p , Q p , s d ) with s p = S -1, Q p = 1 and s d = R.

Van der Laan et al. [START_REF] Van Der Laan | An (s,q) inventory model with remanufacturing and disposal[END_REF] extend the model of Heyman in a different direction. They model explicitly the inventory of remanufacturable products and the remanufacturing process. The repair lead-time is exponential and demand/return processes are again independent Poisson processes. They introduce the (s p , Q p , N ) policy where any return is disposed whenever the number of products waiting for repair equals N . The flow of repaired items is then the output process of an M/M/N/N queue. Analytical expressions are also computed for the average cost. Van der Laan et al. [START_REF] Van Der Laan | Product remanufacturing and disposal: A numerical comparison of alternative control strategies[END_REF] generalize the policies of [START_REF] Ouyang | An inventory control model with disposal for manufacturing/remanufacturing hybrid system[END_REF] and [START_REF] Van Der Laan | An (s,q) inventory model with remanufacturing and disposal[END_REF] with the (s p , Q p , s d , N ) policy. They compare numerically the optimal policy belonging to this class with the optimal policy belonging to (s p , Q p , s d ) (respectively (s p , Q p , N )) class of policies.

Van der Laan et Salomon [START_REF] Van Der Laan | Production planning and inventory control with remanufacturing and disposal[END_REF] consider a model with inter-occurrence times between demands (resp. returns) being Coxian-2 distributed. The demand process and return process are correlated. They compare an (s p , Q p , Q r , s d ) push-disposal policy with an (s p , Q p , s r , S r , s d ) pull-disposal policy to coordinate manufacturing and remanufacturing decisions. For the push-disposal policy, s p is the manufacturing reorder point and Q p the manufacturing order quantity. Returned products are remanufactured with batch size Q r and the remanufactured items are disposed if the inventory position reaches s d . For the pull-disposal policy, remanufacturing is initiated only when the inventory position is below s r and the remanufacturable inventory is above S r . Teunter and D. Vlachos [START_REF] Teunter | On the necessity of a disposal option for returned items that can be remanufactured[END_REF] complement the numerical study of the above model. Unlike [START_REF] Heyman | Optimal disposal policies for a single-item inventory system with returns[END_REF] and [START_REF] Zerhouni | Admission control of product returns in a reverse logistic context[END_REF], other papers do not guarantee that the optimal policy belongs to the class of policies under consideration.

Multi-echelon systems

In their seminal work, Clark and Scarf [START_REF] Clark | Optimal policies for a multi-echelon inventory problem[END_REF] studies a series inventory system with N stages, finite horizon, periodic review, linear holding and backorder cost, no setup cost and stochastic demand only on the downstream stage. They prove that a Base-stock echelon policy is optimal. These assumptions have been relaxed in several papers. For instance, Parker and Kapuscinski [START_REF] Parker | Optimal policies for a capacitated two-echelon inventory system[END_REF] add a constraint of capacity on the orders. In general, a base-stock echelon policy is not optimal when a capacity constraint is added. They show that a base-stock echelon policy is nearly optimal when the downstream echelon is not overloaded. In the contrary case, the Kanban policy is nearly optimal.

In production-inventory systems, replenishment is modelled in a different way than in classical inventory systems. Items are produced by servers one by one, or possibly by batches. Each unit, or batch, requires a random leadtime to be produced. Implicitly, in production-inventory systems, replenishments are capacitated. Veatch and Wein [START_REF] Veatch | Optimal control of a two-station tandem production/inventory system[END_REF] consider a two-stage system with one exponential server at each stage. Otherwise, their assumptions are similar to [START_REF] Clark | Optimal policies for a multi-echelon inventory problem[END_REF]. Again, as in [START_REF] Parker | Optimal policies for a capacitated two-echelon inventory system[END_REF], they show on numerical experiments that the optimal policy is not a basestock echelon policy. They investigate several policy (base-stock, Kanban, fixed buffer and conwip) and compare them to the optimal policy. Their conclusions are the same as those of [START_REF] Parker | Optimal policies for a capacitated two-echelon inventory system[END_REF] : the base-stock policy is generally the nearest optimal policy but when the downstream station is overloaded the Kanban policy is better. In another paper, Veatch and Wein [START_REF] Veatch | Monotone control of queueing networks[END_REF] characterize the form of the optimal policy for a series system with N stages and generalize their results to assembly systems, disassembly systems and to the routing problem. The optimal policy is a state-dependent base-stock policy. It is optimal to produce at stage i when the inventory level at stage i is smaller than a base-stock level, depending on inventory levels at all other stages. Dallery and Liberopoulos [START_REF] Dallery | Comparative modelling of multi-stage production-inventory control policies with lot sizing[END_REF] investigates a Kanban generalized policy. This policy is a mix between Kanban and base-stock echelon policy. They also provide a review of heuristic policies.

Fewer papers investigate the control of multi-echelon systems with product returns. DeCroix et al. [START_REF] Decroix | A series system with returns: Stationary analysis[END_REF] analyze a setting similar to [START_REF] Clark | Optimal policies for a multi-echelon inventory problem[END_REF] except that demand can be negative. They prove that a base-stock echelon is still optimal. They also propose a method to compute a near optimal policy, explain how to extend their model when returns occur at different stages and compare the base-stock echelon policy to fixed buffer policies. With the same type of modelling, DeCroix and Zipkin [START_REF] Decroix | Inventory management for an assembly system with product or component returns[END_REF] characterize the optimal policy for an assembly system and DeCroix [START_REF] Decroix | Optimal policy for a multiechelon inventory system with remanufacturing[END_REF] analyzes a simple series model with only one return in any of the stages and evaluates whether the base-stock model is optimal. His conclusion is that the base-stock echelon policy is optimal if the return is made at the upstream station. Vercraene et al. [START_REF] Vercraene | Two stages production/inventory control with products returns[END_REF] extend [START_REF] Veatch | Optimal control of a two-station tandem production/inventory system[END_REF] for the case with Poisson returns at each stage. The structure of the optimal policy is similar to [START_REF] Veatch | Optimal control of a two-station tandem production/inventory system[END_REF]. In an extensive numerical study, they compare the performances of base-stock, Kanban, fixed buffer and half-optimal policies.

Our contributions

To the best of our knowledge, we are not aware of papers dealing with a multi-stage inventory system with a disposal option upon arrival of returns . We fill this gap by considering an extension to [START_REF] Vercraene | Two stages production/inventory control with products returns[END_REF] where items can be rejected or accepted when they are returned.

In Section 2, we detail the problem formulation. More precisely, each stage is made of an exponential production server and an output inventory. The downstream inventory faces a Poisson demand and each inventory receives product returns, according to Poisson processes, that can be either accepted or rejected. The objective is to find a production/admission policy minimizing the expected discounted costs over an infinite horizon. We consider linear holding, backorder, production, admission and rejecting costs. In Section 3, we show that the optimal policy is characterized by four switching curves: a state-dependent base-stock level and a state-dependent admission threshold for each stage. We also derive monotonicity results for the switching curves.

In Section 4, we show that our results can be extended to the case with two types of returns at each stage: one type has a disposal option while the other one does not have a disposal option. Then we explain how our results can be transposed to an hybrid single-stage system with a production facility and a remanufacturing facility where both facilities have exponential replenishment leadtimes.

In Section 5, we design several simple heuristics and study the performances of some of them. The production servers can be controlled according to classic multi-echelon policies (e.g. Kanban, base-stock echelon, fixed buffer or generalized Kanban). At stage 1, due to the symmetry between production control and admission control, we suggest to use admission controls that are similar to the production controls at the upstream server. At stage 2, it is less clear that we can adapt the production policy to the admission control policy and we propose new control policies. After a brief explanation on how to compute the policies, we compare the performance of some of them in Section 6.

Problem formulation

We consider a 2-stage production/inventory system in series which satisfies end-customer demand (see Figure 1). Station M i produces items one by one. The production lead-time of station M i is exponentially distributed with rate µ i . Preemption is allowed and station M i can be started or stopped at any time. Produced items are stocked in a buffer B i just after M i . The end buffer M 2 sees customer demands arriving according to a Poisson process with rate λ. We assume that backorders are allowed. At time t, the on-hand inventory at B i is denoted by X i (t). Note that X 2 (t) can be negative due to backorders. When buffer B 1 is empty, the production is blocked at station M 2 . Returns of products occur at buffer B i according to an independent Poisson process with rate δ i . Another way to see these two return flows is to consider a situation where there is a single flow of returns for the whole system (rate δ 1 + δ 2 ) and, after an inspection, returned products are routed to the Finished Good (FG) inventory with probability δ 2 /δ 1 + δ 2 and to the Work-In-Process (WIP) inventory with probability δ 1 /(δ 1 + δ 2 ). To ensure the stability of queue lengths, we assume that the demand is smaller than the production capacity for each echelon: λ < µ 1 + δ 1 + δ 2 (echelon 1) and λ < µ 2 + δ 2 (echelon 2). Unlike the system without disposal options [START_REF] Vercraene | Two stages production/inventory control with products returns[END_REF], we don't have to assume that the total return rate is smaller than the demand rate. Upon arrival at buffer B i , a returned product is either accepted with cost C a i or disposed (rejected) with cost

C d i .
When a return is accepted in buffer B i , it can be used immediately as a new product (we neglect the remanufacturing lead-time). Let C r i = C a i -C d i be the relative remanufacturing cost. The unit production cost is

C p i . The system incurs in state (X 1 , X 2 ) a cost rate c(X 1 , X 2 ) = h 1 X 1 +h 2 X + 2 +bX - 2
where h i is the inventory holding cost per unit of time at buffer B i and b is the backorder cost per unit of time.

A production policy π specifies when to produce and when to accept returns, for each stage. The discounted expected cost over an infinite horizon of a policy π, with initial state (x 1 , x 2 ) and discount rate α > 0, is given by :

v π (x 1 , x 2 ) = E   +∞ 0 e -αt c(X 1 (t), X 2 (t))dt|X 1 (0) = x 1 , X 2 (0) = x 2 , π   + E   ∞ i=1 2 j=1 e -αT p j (i) C p j + e -αT a j (i) C a j + e -αT d j (i) C d j |X 1 (0) = x 1 , X 2 (0) = x 2 , π   (1)
where T p j (i), T a j (i) and T d j (i)) respectively represent the i th time instant at stage j when a product is produced, when a return item is accepted and when a return is disposed.

The objective is to minimize the expected discounted cost over an infinite horizon. This problem can be formulated as a continuous-time Markov Decision Process (MDP). Let v ⋆ be the optimal value function:

v ⋆ (x 1 , x 2 ) = min π v π (x 1 , x 2 ) = v π ⋆ (x 1 , x 2 ).
After uniformization of the MDP with rate τ = λ + µ 1 + µ 2 + δ 1 + δ 2 + α, we can transform the continuous-time MDP into a discrete time MDP [START_REF] Puterman | Marckov Decision Processes, Discrete stochastic, Dynamic programming[END_REF]. The optimal value function can be shown to satisfy the following optimality equations:

v ⋆ = T v ⋆ with T v(x 1 , x 2 ) = 1 τ            x 1 h 1 + x + 2 h 2 + x - 2 b + λv(x 1 , x 2 -1) +µ 1 min(v(x 1 , x 2 ), v(x 1 + 1, x 2 ) + C p 1 ) +µ 2 min(v(x 1 , x 2 ), v(x 1 -1, x 2 + 1) + C p 2 ) +δ 1 min(v(x 1 , x 2 ) + C d 1 , v(x 1 + 1, x 2 ) + C a 1 ) +δ 2 min(v(x 1 , x 2 ) + C d 2 , v(x 1 , x 2 + 1) + C a 2 )            (2)

Characterization of the Optimal Policy

In order to derive structural properties of the optimal policy, we will show that the optimal value function belongs to the following set of value functions V .

Definition 1 A value function v belongs to V if, for all x 1 , x 2 , we have

(a) v(x 1 + 1, x 2 ) -v(x 1 , x 2 ) ≤ v(x 1 + 1, x 2 + 1) -v(x 1 , x 2 + 1) (b) (i) v(x 1 + 1, x 2 ) -v(x 1 , x 2 + 1) ≤ v(x 1 + 2, x 2 ) + v(x 1 + 1, x 2 + 1) and (ii) v(x 1 , x 2 + 1) -v(x 1 + 1, x 2 ) ≤ v(x 1 , x 2 + 2) -v(x 1 + 1, x 2 + 1) (c) v(x 1 + 1, x 2 ) -v(x 1 , x 2 ) ≤ v(x 1 + 2, x 2 ) -v(x 1 + 1, x 2 ) (d) v(x 1 , x 2 + 1) -v(x 1 , x 2 ) ≤ v(x 1 , x 2 + 2) -v(x 1 , x 2 + 1)
In Koole [START_REF] Koole | Structural results for the control of queueing systems using event-based dynamic programming[END_REF], property (a) is called supermodularity and denoted Super(1, 2). Property (b) is called superconvexity and denoted SuperC [START_REF] Clark | Optimal policies for a multi-echelon inventory problem[END_REF][START_REF] Dallery | Comparative modelling of multi-stage production-inventory control policies with lot sizing[END_REF]. Finally (c) and (d) refer to the convexity of v in x 1 and x 2 and are denoted Conv(1) et Conv [START_REF] Dallery | Comparative modelling of multi-stage production-inventory control policies with lot sizing[END_REF]. Properties (a) and (b) imply together properties (c) and (d).

Lemma 2 If v ∈ V then T v ∈ V . Moreover v ⋆ belongs to V .
Proof: Due to page constraints, we provide only a sketch of the proof. As operator T is a contraction mapping, the fixed point theorem in a Banach space [START_REF] Puterman | Marckov Decision Processes, Discrete stochastic, Dynamic programming[END_REF] ensures that any sequence of value functions (v n ) defined as v n+1 = T v n will converge to the optimal value function v * , the unique solution of the optimality equations

v * = T v * .
The operator T is a linear combination of operators that propagate Super(1, 2) and SuperC(1, 2) (see [START_REF] Koole | Structural results for the control of queueing systems using event-based dynamic programming[END_REF]). As a result, if a value function v is in V , then the value function T v is also in V . If we take a null value function v 0 , it is clear that v 0 ∈ V . By recursion, we conclude that v * ∈ V . 2

We can now state our main theorem.

Theorem 3

The optimal policy is a state-dependent policy. There exists four switching curves s 1 (x 2 ), s 2 (x 1 ), r 1 (x 2 ) and r 2 (x 1 ) such that :

• Produce at M 1 if and only if x 1 < s 1 (x 2 ). Moreover s 1 (x 2 ) -1 ≤ s 1 (x 2 + 1) ≤ s 1 (x 2 ). • Produce at M 2 if and only if x 2 < s 2 (x 1 ). Moreover s 2 (x 1 ) -1 ≤ s 2 (x 1 + 1) ≤ s 2 (x 1 ). • Accept return at R 1 if and only if x 1 < r 1 (x 2 ). Moreover r 1 (x 2 ) -1 ≤ r 1 (x 2 + 1) ≤ r 1 (x 2 ).
• Accept return at R 2 if and only if x 2 < r 2 (x 1 ). Moreover r 2 (x 1 ) + 1 ≥ r 2 (x 1 + 1) ≥ r 2 (x 1 ). In addition, s 1 (x 2 ) and r 1 (x 2 ) can be ordered for each x 2 in the following way:

s 1 (x 2 ) ≤ r 1 (x 2 ) if C r 1 > C p 1 , s 1 (x 2 ) = r 1 (x 2 ) if C r 1 = C p 1 and s 1 (x 2 ) ≥ r 1 (x 2 ) if C r 1 > C p 1 .
Proof: As v * ∈ V , the four switching curves are well defined. Conv(1) ensures that we can define the production threshold

s 1 (x 2 ) = min[x 1 |v(x 1 + 1, x 2 ) -v(x 1 , x 2 ) + C p 1 > 0] and the disposal threshold r 1 (x 2 ) = min[x 1 |v(x 1 + 1, x 2 ) -v(x 1 , x 2 ) + C r 1 > 0]. SuperC(1,
2) ensures that we can define the production threshold

s 2 (x 1 ) = min[x 2 |v(x 1 -1, x 2 + 1) -v(x 1 , x 2 ) + C p 2 > 0]. Finally Conv(2) ensures that we can define the disposal threshold r 2 (x 1 ) = min[x 2 |v(x 1 , x 2 + 1) -v(x 1 , x 2 ) + C r 2 > 0].
The monotonicity results on the switching curves are also implied by the fact that v * ∈ V . For instance, Super(1, 2) ensures that s 1 (x 2 + 1) ≤ s 1 (x 2 ) and SuperC(1, 2) ensures that s 1 (x 2 ) -1 ≥ s 1 (x 2 + 1). The other monotonicity results are obtained in a similar way. 2

In order to illustrate Theorem 3, we provide a numerical example in Figure 2. It describes the switching curves in function of the state of the system. We can see four curves because there are four decision operators with minimization in T , so four decisions to take to control optimaly the system. The optimal policy is simple, we produce or accept returns in a facility if the system state is below the corresponding switching curve.

For a single-echelon system, the optimal policy is a simple (R, S) policy with S the base-stock level and R the disposal threshold [START_REF] Zerhouni | Intégration des flux inverses dans la gestion des stocks et de la production[END_REF]. x

2 x 1 s 1 (x 2 ) s 2 (x 1 ) r 1 (x 2 ) r 2 (x 1 )
Fig. 2. Optimal policy (µ1 = 0.5, µ2 = 0.7,

δ1 = 0.1, δ2 = 0.6, λ = 1, h1 = 1, h2 = 2, b = 4, C r 1 = -1, C r 2 = 5, C p 1 = C p 2 = 0)
4 Extension to forced returns Consider an extension of our previous model where there is a Poisson return flow without disposal option (rate δ ′ i ) in addition to the Poisson return flow with disposal option (Figure 3). Returns without a disposal option will be called forced returns and returns with a disposal option will be called optional returns. It is straightforward to extend Theorem 3 to this new problem. So the structure of the optimal policy is exactly the same than the problem with disposal option only.

Fig. 3. Extension with forced returns

There are now three necessary and sufficient conditions for the stability of the queue lengths. First, the demand rate must be larger than the sum of forced return rates:

λ > δ ′ 1 +δ ′ 2 .
Second, all forced returns at stage 1 must be treated by station M 2 : δ ′ 1 < µ 2 . Third, each echelon j must be able to satisfy the demand:

λ < µ j + 2 i=j (δ i +δ ′ i ) for j = 1, 2.
Interestingly, we can use this model with forced returns to model an hybrid production/remanufacturing system with positive production and remanufacturing leadtimes (Figure 4). It suffices to set δ 1 = δ ′ 2 = µ 1 = 0 in the model with forced returns to take another viewpoint (consider return as production and vice versa). The optional returns at stage 2 (R 2 ) can be seen as a controllable production facility with production rate δ 2 and preemption allowed. This manufacturing facility produces new products for the finished good inventory B 2 . In parallel, forced returns occur with rate δ ′ 1 . These returns are put in a remanufacturable inventory (B 1 ). Then M 2 is a controllable remanufacturing facility with rate µ 2 . When a product has been remanufactured, it is placed in the finished good inventory B 2 and can be sold as a new product. The cost structure include linear production, remanufacturing, return and holding costs. The problem is then to jointly control the manufacturing server and the remanufacturing server in order to minimize expected discounted costs. Since the optimal policy has a complex structure and is difficult to implement in industry (see Figure 2), we use, in this section, our theoretical results (see Theorem 3) to design simple and efficient policies. We propose the following controls for stations M 1 , M 2 and optional returns R 1 and R 2 .

Control of M 1

• Fixed buffer (F B s1 ) : produce when x 1 < s F B 1

• Base-stock echelon (BS s1 ) : produce when

x 1 + x 2 < s BS 1 • Kanban (KB s1 ) : produce when x 1 + x + 2 < s KB 1
We propose to use classic policies for multi-echelon systems (see e.g. [START_REF] Veatch | Optimal control of a two-station tandem production/inventory system[END_REF]): fixed buffer base-stock echelon and Kanban. The base-stock echelon (BS s1 ) states to produce when x 1 + x 2 < s BS 1 while the Kanban policy (KB s1 ) states to produce when x 1 + x + 2 < s KB 1 . We don't recommend to use a fixed buffer policy (F B s1 ) which states to produce when x 1 < s F B 1 . It is well known that it performs very poorly, most of the time. Control of M 2 :

• Fixed buffer (F B s2 ) : produce when x 2 < s F B

2

The main heuristic used in the literature to control M 2 is the fixed buffer policy F B s2 which states to produce when x 1 + x 2 < s F B 2 . It might be interesting to explore new heuristics to have a better approximation of s 2 (x 1 ) since this switching curve might depend a lot on x 1 . Control of R 1 :

• Fixed buffer (F B r1 ) : accept return when x 1 < r F B 1

• Base-stock echelon (BS r1 ) : accept return when

x 1 + x 2 < r BS 1 • Kanban (KB r1 ) : accept return when x 1 + x + 2 < r KB 1
Due to the equivalence between the control of M 1 and the control of optional returns R 1 , we suggest to use similar policies to control R 1 and M 1 : base-stock echelon (BS r1 ) which states to accept return when

x 1 + x 2 < r BS 1
and Kanban (KB r1 ) which states to accept return when

x 1 + x + 2 < r KB 1 . Control of R 2 :
• Fixed buffer (F B r2 ) : accept return when x 1 < r F B

2

• Base-stock return (BR r2 ) : accept return when

x 1 + x 2 < r BR 2
• Base-stock return with priority (BP r2 ) : accept return when x 1 + x 2 < r BP 2 or x 2 < 0 Fixed buffer (F B r2 ) which states to accept return when x 2 < r F B 2 is a very simple heuristic because it considers only local informations. However, we see in section 3 that optimal r 2 (x 1 ) decreases. We have an intuition about that : if the server M 2 is fast, stock B 1 and stock B 2 could be seen as one global stock. In the limit case (µ 2 → ∞) the system is like the single echelon of [START_REF] Zerhouni | Intégration des flux inverses dans la gestion des stocks et de la production[END_REF]. Due to this observation a new heuristic called base-stock return (BR r2 : accept return when x 1 + x 2 < r BR 2 ) seems relevant. But, an other limit observation have to be noted : if we have a lot of product in B 1 the system could be seen as a single echelon with an infinite stock B 1 , it assures that lim x1→∞ r 2 (x 1 ) = constant. From this observation we propose an heuristic called base-stock return with Priority (BP r2 ) which states to accept return when x 1 + x 2 < r BP 2 or x 2 < 0. The BP r2 heuristic could be extend to a policy (obviously better) with a variable level of priority x 2 < r ′BP

2

. This heuristic would be similar to the generalized Kanban [START_REF] Dallery | Comparative modelling of multi-stage production-inventory control policies with lot sizing[END_REF]. Here we chose 0 as priority level because we think it is the easiest to implement in industry.

Moreover we propose to try specific cases with s i = r i (i ∈ {1, 2}) and with the same heuristic for the control of M i and R i . We call these cases equal policies. They are obviously worse but easier to implement in industry. Optimal policies s 2 (x 1 ) (insceased) and r 2 (x 1 ) (decreased) δ' 1 =0.1 δ' 1 =0.9 In this paper, we prove the structure of the optimal policy of a two-stage production/inventory system with returns at each stage that could be disposed or not. Due to the complexity of the optimal policy, we propose heuristics policies to simplify their implementation in industry. At the end, we make an numerical study to test their efficiency. Our results are complementary to those in the literature but it would be interesting to extend the numerical study and investigate efficient heuristics to control the downstream server.

Fig. 1 .

 1 Fig. 1. Model

Fig. 4 .

 4 Fig. 4. Hybrid production/remanufacturing system

Fig. 5 .

 5 Fig. 5. Variation of µ2 for the set of parameters : µ1= 0.6, δ ′ 1 = 0, δ2 = 0.6, λ = 1, h1 = 1, h2 = 2, b = 4, C r 2 = 10, C p 1 = C p 2 = 0

Fig. 6 .

 6 Fig. 6. Variation of δ ′1 for the set of parameters :µ1 = 0, µ2 = 1.5, δ2 = 1, λ = 1, h1 = 1, h2 = 2, b = 4, C r 2 = 10, C p 1 = C p 2 = 0 7 Conclusion

Numerical study

In this numerical study, we focus on the average cost criterion. The optimal policy is computed by a value iteration algorithm [START_REF] Puterman | Marckov Decision Processes, Discrete stochastic, Dynamic programming[END_REF] based on the fact that any sequence v n+1 = T v n converges to the optimal value function. We truncate the state space in three directions and increase the state space until the average cost is no more sensitive to the truncation level. Heuristic policies are also computed by a value iteration algorithm where we adapt the operator T to the heuristic under consideration.

For the heuristic policies, we want to find the parameters s 1 , s 2 , r 1 , r 2 that minimize the average cost function C(s 1 , s 2 , r 1 , r 2 ) which is computed by the previous value iteration algorithm. This optimization problem is a non linear problem with integer variables that might be long to solve since evaluating a given policy might already take time. Therefore, we make the plausible assumption that the function C(s 1 , s 2 , r 1 , r 2 ) is unimodal. A function f is unimodal if for x, y and z on a line and y between x and z : f (x) is finite and f (x) ≤ f (y) implies f (y) ≤ f (z). The assumption of unimodularity has been validated on several instances. Based on this assumption, we can solve the system with the maximal gradient with constant step method. This method is very efficient here because we can start the optimization with an approximate value of s 1 , s 2 , r 1 and r 2 resulting from the calculation of the optimal policy. Note that if we want to valid that a cost C(s 1 , s 2 , r 1 , r 2 ) is a minimum we need to check that all adjacent points of (s 1 , s 2 , r 1 , r 2 ) are higher. With 4 dimensions we need to check 3 4 = 81 points, so the computation time can be long.

We limit our study to a system with δ 1 = δ ′ 2 = 0 to have only one flow which could be forced or/and controllable from production and flow which could be controllable from return. Moreover we make the classical assumption that h 1 = λ = 1. Under these hypotheses we compute several instances.

As expected, optimal policy respect the structure described in Section 3, the relative position between s 1 (x 2 ) and r 1 (x 2 ) is function of the sign of C r 1 -C p 1 and when M 2 is fast, we have r 2 (x 1 ) near curve x 1 + x 2 = constant. Note that we do not find a clear explanation of the relative position between s 2 (x 1 ) and r 2 (x 1 ) (see Figure 5). About heuristics, the first observation is that equal heuristics (i.e. heuristics with r π i =s π i ) are obviously worse than others heuristics and even often very bad with a cost gap over > 10%. Their implementation in industries is easier but the loss can be very important. The second observation is that, for control of M 1 , we obtain the same conclusion as [START_REF] Veatch | Optimal control of a two-station tandem production/inventory system[END_REF]. The Kanban policy on M 1 is better than others when M 2 is overloaded with regard to µ 1 . Ordinary the best heuristics for M 1 is the base stock echelon policy. The fixed buffer policy is numerically always the worse. The optimal thresholds plot in Figure 5 show this trend. The last observation is that the control of R 2 is more complicated. We give an more specific study about this control in the following paragraph.

We compute instances with heuristic BS s1 for the control of M 1 and F B s2 for the control of M 2 given Figures 5 and6. The first figure, shows the variation of production capacity in server M 2 . The rate µ 2 can vary on ]0.6, ∞[. Here we start from 0.65 because compute near limit of stability increased the state space and the computation time. F B r2 is the best heuristic when M 2 is overloaded (i.e. µ 2 near 0.6). But BR r2 and BP r2 are very efficient when utilization of M 2 is low (i.e. µ 2 >> λ). This phenomenon is expected because optimal curves r 2 (x 1 ) on the three other plots have this behavior. When µ 2 = 0.65, r 2 (x 1 ) is near horizontal so F B r2 threshold is a good approximation. But when µ 2 is bigger, r 2 (x 1 ) is near curve x 1 + x 2 = constant so BR r2 heuristic became the best. Note that s 2 (x 1 ) increases a lot when M 2 is overloaded. We do not have a clear intuition about this but a heuristic with threshold which increasing in x 1 could be tested in this case. Moreover, note that peaks and irregularity are due to discrete character of parameters s i and r 2 .

The BP r2 policy seems a good heuristic, but it is not always the case. For instance, if we have a lot of forced returns, then r 2 (x 1 ) became negative for all x 1 . So BP r2 became a bad approximation of optimal threshold because r BP 2 can not be negative. To solve this problem we could relax the constraint that priority level is 0 and look for the optimal priority level. Figure 6 shows the variation of forced return in stock B 1 . The return rate δ ′ 1 can varies on [0, 1[. This figure shows clearly the limit of BP r2 when optimal r 2 (x 1 ) < 0. Nevertheless for δ ′ 1 < 0.7 it is generally a good heuristic.


Stephane Chretien
email: stephane.chretien@univ-fcomte.fr

Jean-Marc Nicod

Laurent Philippe
email: laurent.philippe@femto-st.fr

Veronika Rehn-Sonigo
email: veronika.sonigo@univ-fcomte.fr

Lamiel Toch
email: lamiel.toch@gmail.com

Job Scheduling Using successive Linear Programming Approximations of a Sparse Model

In this paper we tackle the well-known problem of scheduling a collection of parallel jobs on a set of processors either in a cluster or in a multiprocessor computer. For the makespan objective, i.e., the completion time of the last job, this problem has been shown to be NP-Hard and several heuristics have already been proposed to minimize the execution time. We introduce a novel approach based on successive linear programming (LP) approximations of a sparse model. The idea is to relax an integer linear program and use ℓp norm-based operators to force the solver to find almost-integer solutions that can be assimilated to an integer solution. We consider the case where jobs are either rigid or moldable. A rigid parallel job is performed with a predefined number of processors while a moldable job can define the number of processors that it is using just before it starts its execution. We compare the scheduling approach with the classic Largest Task First list based algorithm and we show that our approach provides good results for small instances of the problem. The contributions of this paper are both the integration of mathematical methods in the scheduling world and the design of a promising approach which gives good results for scheduling problems with less than a hundred processors.

Introduction

Nowadays clusters of computers or large shared memory computers are widely used by many communities such as researchers, universities or industries to speed up their applications. Due to their cost these computing facilities are usually shared between several users and several parallel jobs must be run at the same time on the same platform. The problem of scheduling parallel jobs on clusters without knowing in advance the submission times of user jobs has been widely studied [START_REF] Ye | On-line scheduling of parallel jobs in a list[END_REF]. In this case the scheduling problem is said to be "on-line" [START_REF] Dror | On the definition of "on-line" in job scheduling problems[END_REF]. When all characteristics of the jobs are known in advance, the scheduling problem becomes "off-line" and it has been widely studied for sequential jobs [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF] and for parallel jobs [START_REF] Dutot | Bicriteria algorithm for scheduling jobs on cluster platforms[END_REF][START_REF] Dror | Job scheduling in multiprogrammed parallel systems[END_REF].

The "off-line" problem considered here depends on the job characteristics. In the literature one distinguishes three kinds of parallel jobs. Rigid jobs [START_REF] Lublin | The workload on parallel supercomputers: Modeling the characteristics of rigid jobs[END_REF] are performed with the number of processors originally required. Moldable jobs introduced by Turek et al. in [START_REF] Turek | Approximate algorithms scheduling parallelizable tasks[END_REF] may run with different numbers of processors but cannot change their allocation after their start. Malleable jobs [START_REF] Dutot | Scheduling on hierarchical clusters using malleable tasks[END_REF] can modify the number of allocated processors during their execution. The rigid job model can easily be used in most of the cases of parallel jobs. The two other models however need an interaction between the application and the scheduler to define the number of allocated processors. This is for instance the case of applications developed with the Bulk Synchronous Parallel (BSP) model introduced in [START_REF] Valiant | A bridging model for parallel computation[END_REF] that can be run as moldable jobs. Processor virtualization however could be a solution to transparently make standard parallel applications moldable as presented in [START_REF] Nicod | Using virtualization and job folding for batch scheduling[END_REF]. Applying virtualization to malleable jobs is probably more difficult as it would need to use virtual machine migration. For these reasons we focus on rigid and moldable jobs.

The problem of scheduling several parallel rigid and moldable jobs on homogeneous computing resources has been shown to be NP-Hard respectively in [START_REF] Dror | Job scheduling in multiprogrammed parallel systems[END_REF] and [START_REF] Dutot | Bicriteria algorithm for scheduling jobs on cluster platforms[END_REF]. Several previous works have already tackled the issue of providing heuristics that give efficient sub-optimal solutions. In [START_REF] Krim Amoura | Scheduling independent multiprocessor tasks[END_REF] static scheduling of rigid parallel jobs for minimizing the makespan is studied and in [START_REF] Foto | Scheduling to minimize the average completion time of dedicated tasks[END_REF] for minimizing the sum of the completion time of each job. In [START_REF] Dutot | Scheduling moldable BSP tasks[END_REF], Dutot et al. consider the problem of scheduling moldable jobs with the objective of minimizing the makespan. The authors present experimental results where the well-known Largest Task First (LTF) algorithm is the best for the makespan objective.

The contribution of this paper is a novel approach for scheduling a collection of rigid or moldable jobs using successive LP approximations based on the gradient operator. To the best of our knowledge there is no existing work using this promising approach based on the sparse recovery problem in statistics domain.

The remainder of the paper is organized as follows. In Section 2 we describe the problem and the model of moldable jobs. In Section 3 we present the sparsity promoting penalization as well as linear approximation principles. Then, in Section 4 we present how to adapt this method to our scheduling problem. In Section 5 we compare our technique with the algorithm developed by Dutot et al. in [START_REF] Dutot | Scheduling moldable BSP tasks[END_REF] and show experimental results to assess the performance of our approach, and finally we conclude and give future work directions in Section 6.

Framework

In this section we formally define the targeted framework and the problem. We consider the problem of scheduling a collection of n independent parallel jobs. We tackle both cases of rigid and moldable jobs.

The jobs are run on a homogeneous cluster of distributed computing nodes or on a shared memory multiprocessor or multicore computer. In a cluster each node is made up of identical processors which are in turn made up of identical cores. The scheduling policy used on most clusters does not pay any attention to the exact distribution of the cores allocated on the nodes provided that the job is parallel. For this reason, in this paper, we will only consider the number of allocated cores, assimilated to processors and called Processing Elements (PEs). The results can then be applied either on clusters or on multiprocessor-multicore computers. In the remainder of the paper m denotes the number of available PEs in the execution platform.

Rigid jobs are defined by an execution time and a static number of requested PEs, i.e., the job cannot be run on neither more nor less PEs than originally requested. Each rigid job i is defined by its number of requested PEs reqproc i and its duration reqtime i .

Moldable jobs can be run on a different number of PEs or cores but this number is fixed at the job execution start and cannot change during the execution. The considered moldable jobs respect the model defined in [START_REF] Dutot | Scheduling moldable BSP tasks[END_REF]. Let reqtime i be the duration of job i which requires at most reqproc i PEs. Let t i (n) be the duration of the job i if n PEs are allocated for job i. The relation between the duration of a job i and its number of allocated PEs is stated as:

∀i, ∀n ≤ reqproc i , t i (n) = reqproc i n reqtime i
Given this framework our objective is to minimize the makespan of the schedule. According to the α|β|γ (platform | application | optimized criterion) classification of scheduling problems given by Graham in [START_REF] Lewis | Optimization and approximation in deterministic sequencing and scheduling: a survey[END_REF], the above problem is denoted by P |parallel jobs|C max .

Sparsity promoting penalization with successive linearizations

The optimization method presented in the paper relays on two steps. First we formulate the problem as an integer linear program, then we relax it and apply the sparsity promoting penalization which tries to find almost integer solutions. As the sparsity promoting penalization implies to minimize a non linear objective function we use successive LP approximations to linearize it. In this section we detail the main steps of the method.

Sparsity promoting penalization

Recent works on the sparse recovery problem in statistics and signal processing have brought to general attention the fact that using non-differentiable penalties such as the ℓ p norm can be an efficient ersatz to combinatorial constraints in order to promote sparsity. This approach for constructing continuous relaxations to hard combinatorial problems is a key ingredient in e.g., the new field called Compressed Sensing which originated in the work of Candès, Romberg and Tao [START_REF] Candes | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF]. Donoho [START_REF] David | Compressed sensing[END_REF] showed that finding the sparsest solution to an underdetermined system of linear equations may sometimes be equivalent to finding the solution with smallest ℓ 1 -norm. This discovery lead to a intense research activity in the recent years focusing on finding weaker sufficient conditions on the linear system under which it is possible to prove this equivalence. It was found in particular that for matrices satisfying certain incoherence conditions (implying that the columns of the associated matrix are almost orthogonal), the equivalence between finding the sparsest and the least ℓ 1 norm solution holds for systems with a number of unknowns to the order of exponential of the number of equations. Other non-differentiable penalties have also been proposed in order to increase the performance of sparse recovery procedures. Candès, Wakin and Boyd proposed an iterative reweighted ℓ 1 procedure in [START_REF] Candes | Enhancing Sparsity by Reweighted L1 Minimization[END_REF]. In our setting, the standard ℓ 1 relaxation is not suitable. Indeed, as will be detailed in the sequel (e.g. equation 1 below), our constraints will always imply that the ℓ 1 is constant. A more appropriate sparsity promoting penalization in this case is the ℓ p -quasi-norm relaxation, for p ∈ (0, 1). This corresponds to minimizing x p := ( k x p k ) 1/p instead of x 1 , under the same design constraints. Such a non-convex relation was successfully implemented in, e.g. [START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF].

Linear and conic approximation

In physics and mathematics a function f is often approximated with a linear formulation at point x 0 , if f is differentiable at point x 0 . The gradient of a function with several parameters (f : R n → R), noted ∇f , is the vector whose components are equal to derivatives of f with respect to the parameters. Taylor's expansion gives

f (x + h) = f (x) + ∇f, h + o(h)
where x and h belong to R n , and represents the dot product. In cases such as x → x p , where f is non-differentiable, it is still possible linearize by using the appropriate generalization of the gradient, called the Clarke-subdifferential. In simple words, a non-differentiable function may have several tangents in a generalized sense and the Clarke-subdifferential, denoted by ∂f (x), is the set of all such generalized tangents. The nonsmooth counterpart to Taylor's expansion is given by

f (x + h) = f (x) + sup g∈∂f (x) g, h + o(h) .
In order to implement our ℓ p -based relaxation, we will implement successive linearizations on a standard linear programming solver.

Applying the method on the job scheduling problem

In this section we apply the method on the job scheduling problem. First it implies to define a sparse representation of the problem then we apply the two steps of sparsity promoting penalization and linear approximation.

Formulation as an integer linear program

In the defined framework a solution to the scheduling problem must provide at least the start time of the jobs for the rigid jobs as their duration and the number of used PEs are constants of the problem. For the moldable jobs, the duration depends upon the number of PEs that are allocated to the jobs. So the scheduled jobs are characterized by their start time and the number of allocated PEs and the duration of the job is determined as soon as this number of allocated PEs is determined. We call "configuration" of a job the number of allocated PEs. We call "position" of a job, its position determined by its start time in a discrete time scale. Finally, we call "slot" the couple (configuration, position). Let us create a list of slots (configurations, positions) for each job. The idea is to create a vector x i for each job i. Each component x i,j of the vector x i is a binary variable which indicates whether slot j of job i is chosen or not. Then we fix a time horizon T and we let a linear program find a solution. We iteratively reduce the time horizon T until the linear program cannot find a solution any more.

The following constant values are defined to formulate the problem:

• proc i,j : the number of PEs for the configuration j of job i

• nconf i : the number of all possible configurations for job i (for rigid jobs i, nconf i = 1)

• nslots i : the number of all possible slots for job i

• C i,s : the configuration index of job i used in the slot s of job i

• run i,s,t indicates whether in the slot s the job i is running at time t.

Then we define the binary variable x i,s which indicates whether slot s of job i is chosen or not. For each job i, we note x i the vector whose components are the values x i,s and we define a vector x which is equal to the concatenation of the n vectors x i of every job i,

1 ≤ i ≤ n.
The problem can be formulated as an integer linear program. Since we only have to determine whether a feasible solution exists or not for a given time horizon T , we only need the constraints to be respected. That is why we set all the coefficients of the variables in the objective function to 0. The problem is stated as: "find a feasible solution which respects the following linear constraints:"

∀1 ≤ i ≤ n, s=nsloti s=1 x i,s = 1 ( 1
)
∀1 ≤ t ≤ T , i=n i=1 s=nsloti s=1 x i,s × run i,s,t × proc i,Ci,s ≤ m (2)
Constraint 1 imposes the unicity of the chosen slot s on each job i. Constraint 2 means that at each time t, the set of all running jobs does not consume more than the m available PEs in the considered cluster.

Relaxation via sparsity promoting penalization

The solution of the Integer formulation of the problem cannot be found in polynomial execution time. So we make a relaxation of it. We transform all binary variables x i,s into rational variables and with 0 ≤ x i,s ≤ 1. Since vector x -the concatenation of x i vectors -indicates which slots are chosen, we are tempted to strongly enforce its sparsity. In fact, vector x must have exactly n "1" and many "0". Thus, we legitimately expect the binary constraints to be naturally recovered by imposing sufficient sparsity. Notice that the proposed constraints impose that the sum of the components of x is equal to one jobwise. Since the components are positive, this implies that the ℓ 1 norm is equal to one jobwise, which explains why minimizing the ℓ 1 norm for promoting sparsity is unfortunately useless in the present context. In order to overcome this difficulty, we chose to minimize the ℓ p norm non-convex function

f (x) = i x i p (3)
under constraints (1) and (2) for p ∈ ]0, 1[.

Successive LP approximation scheme

We now apply successive LP approximation schemes to linearize the problem. Let f i (x i ) = x i p , for all jobs i. Thus, f = i f i . We use the value of each variable computed during the previous iteration. We will use the following arbitrary choice g ∈ ∂f among all possible subgradients of f :

g i,j =      x p-1 i,j × f i (x i ) 1-p if x i,j = 0 0 otherwise. ( 4
|J | i f (x (iter) i ) + ∇f (x (iter) i ), xi -x (iter) i 8
x ← execute LP(J , m, T, proc, run) The method is implemented in Algorithm 1. It starts with any initial value e.g. the zero vector. First we compute a lower bound of the makespan at line 1, which is equal to the maximum between the duration of the longest job and i reqproci×reqtimei m

. The time horizon T is set to the makespan of the LTF list algorithm. If the linear program LP finds a satisfactory solution (line 9), it reduces the time horizon (line 12) until it cannot (line 23) before maxIter iterations. If it does not find a satisfactory solution with T = Listmakespan before maxIter iterations (line ), it increases the time horizon T (line 21). For a given time horizon T , it iteratively updates the objective function of the linear program (line 7) according to the subgradient-based Taylor approximation rule of the sparsity promoting penalization.

Improving the algorithm efficiency

During the experiment step of our work a problem appeared in the linear resolution. Satisfactory solutions for Algorithm 1 are only detected (at line 9) if all jobs i have their vector x i with exactly one "1" as the algorithm is designed to find exclusively exact solutions. In fact, for a given time horizon T , the successive linear approximations manage to find a schedule for most of the jobs of the collection but it let few jobs j of the collection with fuzzy schedules. That is to say, vectors x i contain exactly one "1" while vectors x j do not. In this case the algorithm often continues to iterate, even if x j is close to 1, until maxIter is reached without being able to find a solution. This leads to longer computing times for the algorithm while giving inefficient solutions.

So we modify Algorithm 1 and its detection criterion at line 9 as follows: when a valid rational schedule is found we keep the exact schedule for jobs i whose x i have exactly one "1" and we schedule the rest of the jobs for which the linear program gives fuzzy schedules with the LTF list algorithm. If a solution shorter than the time horizon T is found, then the f ound variable is set to true otherwise we continue to iterate.

Simulation and results

In this section we present the results obtained on the two versions of the algorithm and we compare them to the well-known Largest Task First algorithm. We assess both cases of rigid and moldable jobs. Notice that the problem we propose to solve is nonconvex and very high dimensional. Moreover, no theoretical guarantee for convergence of the proposed iterative procedure is available and it is well known that minimizing an ℓ p quasi-norm, 0 < p < 1, is NP-hard already. On the other hand, various non-convex ℓ p -based strategies have been successfully used for promoting sparsity in the literature. Despite the current lack of appropriate theoretical foundation, in most reported experiments the ℓ p -based approach managed to reach a local solution significantly superior to the ℓ 1 minimizer for, e.g., the Compressed Sensing reconstruction problem [START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF]. The goal of this section is to show that such a good performance can also be observed for the studied scheduling problem. Notice that we did not optimize the computational aspects of the problem, in particular, we made no use of the very special properties of the constraint matrix. This explains why the computing time is currently much higher than what could be obtained after a careful design of the algebraic aspects of our algorithms.

Experimental settings

Carrying out real experiments on clusters is difficult: experiments are not reproducible and may be long. Furthermore a cluster is expensive and meant to be used for calculations while experiments may monopolize it. For these reasons, we have developed a simulator of a homogeneous cluster based on a master/slave architecture. This simulator is also meant to check schedules obtained by the different algorithms. The simulator is implemented using SimGrid [START_REF] Casanova | Simgrid: A generic framework for large-scale distributed experiments[END_REF] and its MSG API. It takes a workload as input and it gives a schedule as output.

To simulate the job collection, we use synthetic workloads generated with uniform distributions. The parameters associated with a workload is the job granularity, the ratio of the duration of the longest job over the duration of the shortest one.

Assessing performance of Algorithm 1

In a first set of experiments the simulations have been run with a ℓ p norm where p = 0.1, maxIter is set to 15000 in the algorithm and the machine is made up of 64 PEs. We have scheduled a collection of 60 jobs and, for each number of jobs in the collection, we performed 40 experiments to compute an average value of the ratio of the makespan over the lower bound. The results where disappointing: they were far from the optimal and very time consuming.

So we ran another set of experiments with less jobs, ℓ p norm where p = 0.1, maxIter set to 15000. The machine is made up of 32 PEs and the number of PEs requested by each job is uniformly chosen between 1 and 8. The granularity is set to 25. For each number of jobs in the collection we perform 20 experiments, then we remove the best and the worse results in order to reduce the deviation, and we compute an average of the ratio of the makespan over the lower bound. Figure 1a shows the ratio of the makespan over the lower bound against the number of rigid jobs, while Figure 1b shows this ration for moldable jobs. In the figures the algorithms are noted succ. LP approx for our algorithm and LIST for the LTF implementation. The figures also show the standard deviation σ: the height of a vertical line is equal to 2σ.

We can note that for less than 25 jobs the successive LP approximation algorithm gives better results than the LTF algorithm and with more than 25 jobs the latter outperforms successive LP approximation. Note that after 40 jobs the performance ratio of LTF quickly tends toward 1.1 which means, on the one hand, that it probably finds most of the time the optimal solution and, on the other hand, that it is difficult to find better solutions. Moreover, with more than 25 jobs, the problem becomes so complex that the successive LP approximation algorithm must increase its T to find a solution. Under this threshold the maximum gain is about 15% for 16 moldable or rigid jobs. Furthermore the standard deviation of the experiments with our new approach is less than the standard deviation of LTF. We can easily understand that for 5 jobs the optimal is found due to the experimental settings: the number of PEs that each job requires is uniformly chosen between 1 and 8. As a consequence, all jobs may start at time 0. That also explains the peak with 15 jobs which do not necessarily start at time 0. We can also notice that for rigid jobs and moldable jobs the successive LP approximation algorithm has the same behavior, that is to say, when the number of jobs increases, the ratio of the makespan over the lower bound increases. As we can see in Figure 2a and Figure 2b, Algorithm 1 is very time consuming with both rigid jobs and moldable jobs compared to the LTF algorithm. Note however that for 15 jobs, in the case where gives the best results, the time taken by the LP approximation is not more than 1.5 minutes which is still reasonable. We assess the performance of the improved version in the following section.

Performance of the improved algorithm

To assess the improved algorithm, we performed experiments with two simulated machines made up of 64 and 128 PEs. The number of PEs requested by each job is uniformly chosen between 1 and 16 for the machine with 64 PEs, and between 1 and 32 for the machine with 128 PEs. Granularity is set to 25 for the machine with 64 PEs and to 10 for the machine with 128 PEs. We set p = 0.1 and maxIter = 200. For each number of jobs in the collection, we perform 40 experiments.

Figure 3a shows the ratio of the makespan over the lower bound against the number of rigid jobs in a cluster of 64 PEs, while in Figure 3b we consider scheduling moldable jobs. The performances of the new approach are better than LTF for moldable and rigid jobs, and better than the unmodified algorithm. Figure 4a and Figure 4b give the results for a machine with 128 PEs.

We can note that in the four cases the performance ratio between LTF and our approach is up to 20%. The results obtained with a 128 PEs machine show an improvement for the LTF algorithm compared to 64 case while the behavior of the new algorithm is quite similar. This is probably because our solution is very close to (if not at) the optimal solution and nothing more can be gained.

We have also recorded some statistics data after each execution of the linear program. On average with 64 PEs and 16 jobs almost 75% of jobs have exact schedules, while with 64 jobs 50% of them have exact schedules. We notice that when the number of jobs to schedule increases the number of exact schedules found by the linear program decreases. We get the same behaviour with 128 PEs: on average with 128 PEs and 16 jobs almost 80% of jobs have exact schedules, while with 128 jobs 60% of them have exact schedules.

Figures 5a and5b show the time spend by the two algorithms. We notice that the hybrid algorithm is less time consuming than the original algorithm but still consumes more time than the LTF algorithm.

Conclusion and Future Work

In this paper, we assess the use of successive linear programming approximations of a sparse model for job scheduling. This method is applied on clusters to schedule rigid and moldable jobs. Experimental results show that the pure successive LP approximation only gives good performances regarding the makespan for scheduling up to dozens jobs on a machine with dozens PEs. In contrast, a variant associated with LTF gives good results for bigger instances with up to a hundred jobs on machines with up to a hundred PEs. This variant is a good For future work we plan to implement the Split Bregman Method [START_REF] Goldstein | The split bregman method for l1-regularized problems[END_REF] to speed up the solving time and try other relaxations of the linear program. We also plan to use a multi-level scheduling approach for which we distinguish small jobs and large jobs. We then apply our method on different collections of jobs.

An important part of the simulations has been run thanks to the computing facilities of the Mésocentre de Calcul de Franche-Comté in Besançon, France.

4 proc← 5 ∀i 6 while

 456 compute the configurations (J ) ; run ← compute all possible slots (J , m, T ) ; iter ← 1 ; f ound ← f alse iter < maxIter and not f ound do 7 set the objective function of LP(J , m, T, proc, run) to

9 if 11 T 12 T ← T - 1 13f

 911121 ∀i, xi contains exactly one "1" then 10 sched ← convert into schedule (x, proc, run) ← makespan(sched) ound ← true 14 if incT = true then 15 end ← true 16 ∀i, k, x (iter) i,k ← x i,k ; iter ← iter + 1 17 if not f ound then 18 if T = listM akespan then 19

Figure 1 :

 1 Figure 1: Performance comparison 32 PEs

Figure 2 :

 2 Figure 2: Compute time for 32 PEs

Figure 3 :

 3 Figure 3: Performance of the algorithms with 64 PEs

Figure 4 :

 4 Figure 4: Performance of the algorithms with 128 PEs

Figure 5 :

 5 Figure 5: CPU Time consumed to compute a schedule with 64 PEs


K C Santosh
email: santosh.kc@loria.fr

Bart Lamiroy
email: bart.lamiroy@loria.fr

Laurent Wendling
email: laurent.wendling@parisdescartes.fr

Spatio-structural Symbol Description with Statistical Feature Add-on

In this paper, we present a method for symbol description based on both spatio-structural and statistical features computed on elementary visual parts, called 'vocabulary'. This extracted vocabulary is grouped by type (e.g., circle, corner ) and serves as a basis for an attributed relational graph where spatial relational descriptors formalise the links between the vertices, formed by these types, labelled with global shape descriptors. The obtained attributed relational graph description has interesting properties that allows it to be used efficiently for recognising structure and by comparing its attribute signatures. The method is experimentally validated in the context of electrical symbol recognition from wiring diagrams.

Introduction

Graphics recognition has an extremely rich state-of-the-art literature in symbol recognition and localisation. However, most methods are particularly suited for isolated line symbols, not for composed symbols connected to a complex environment [START_REF] Cordella | Symbol recognition in documents: a collection of techniques?[END_REF][START_REF] Lladós | Symbol Recognition: Current Advances and Perspectives[END_REF]. Considering the problem of symbol localisation in real documents, composed of individual parts and constrained by spatial relations for instance, one needs to be able to extract visual parts, characterise their shape description and formalise the possible links that exist between them. This integration of spatial relations and shape description of the extracted visual parts is going to be the core of this paper. The method is very much inspired by a real world industrial problem [START_REF] Tombre | Pattern recognition methods for querying and browsing technical documentation[END_REF][START_REF] Santosh | Inductive logic programming for symbol recognition[END_REF][START_REF] Tooley | Aircraft electrical and electronic systems: principles, operation and maintenance. Aircraft engineering principles and practice[END_REF]. Fig. 1 shows a few samples of the data related to it.

Global signal based descriptors [START_REF] Zhang | Review of shape representation and description techniques[END_REF] present a number of inconvenients in our context. They difficultly accommodate with connected or composite symbols and they are generally not well adapted for capturing small detail changes. In statistical approaches, signatures are simple with low computational cost. They are, unfortunately, primarily designed for applications where line symbols are isolated [START_REF] Tombre | Stable, robust and off-the-shelf methods for graphics recognition[END_REF]. Furthermore, discriminative power and robustness in general applications usually require optimal selection of features [START_REF] Barrat | A bayesian network for combining descriptors: application to symbol recognition[END_REF] or the fusion of different classifiers [START_REF] Terrades | Optimal classifier fusion in a nonbayesian probabilistic framework[END_REF]. Besides global signal based descriptors, another idea is to decompose the symbols into either vector based primitives like points, lines and arcs or into meaningful parts like circles, triangles and rectangles. These primitives are then used in structural descriptors like attributed relational graphs (ARG) [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF], region adjacency graphs (RAG) [START_REF] Lladós | Symbol recognition by error-tolerant subgraph matching between region adjacency graphs[END_REF], as well as deformable templates [START_REF] Valveny | A model for image generation and symbol recognition through the deformation of lineal shapes[END_REF]. In addition to the common drawback related to stability issues coming from segmentation and, error-prone raster-to-vector conversion, variability of the size of the underlying graph structures leads to computational complexity in matching. Structural approaches however, provide a powerful relational representation, conveying how parts are connected to each other, and are usually considered to be preserving generality and extensibility.

In this paper, we aim to combine the best of both structural and statistical approaches, and try to avoid the shortcomings of each of them. To do so, we decompose symbols by expressing their various parts in a fixed visual vocabulary, using spatial relations, graphs and signal based descriptors to describe the whole shape. Our symbol description is explained in Section 2. Symbols can be compared by computing matching scores based on vertex and edge alignment (cf. Section 3). In Section 4, we validate our method and compare it with the state-of-the-art. The paper concludes in Section 5. This paper is the extension of the previously published work [START_REF] Santosh | Symbol recognition using spatial relations[END_REF], where we have validated the use of spatial relations for symbol recognition, but which did not include the vertex signatures developed here.

Symbol Description

Expanding on previously published work [START_REF] Tombre | Pattern recognition methods for querying and browsing technical documentation[END_REF][START_REF] Santosh | Symbol recognition using spatial relations[END_REF], we use a set of well controlled elementary part detectors to define a visual vocabulary. In our case, they consist of: circles, corners, loose ends and thick (filled) components. More formally, we denote the type set as, T = {T thick , T circle , T corner , T extremity }. Such visual elementary parts are extracted with the help of image treatment analysis operations [START_REF] Rendek | The search for genericity in graphics recognition applications: Design issues of the qgar software system[END_REF]. While, in the general case, this vocabulary can be of any kind from any type of bag-of-features, related to what is visually pertinent in the application context under consideration, our current vocabulary is related to electrical symbols. It can be easily extended to adapt to other domains.

Rather than using the detected elements as a basis for expressing and computing spatial relations, we group them together according to their types, as shown in Fig. 2. Now, we represent whole symbol by a complete ARG as a 4tuple G = (V, E, F A , F E ) where V is the set of vertices, E ⊆ V × V is the set of graph edges, F A : V → A V is a function assigning attributes to the vertices and A V representing a set of vocabulary type set T as well as their global shape signatures S, and F E : E → ℜ E is a function assigning labels to the edges where ℜ represents spatial relations of the edge E. Following the symbol in Fig. 2, the resulting graph whose attribute type set is {T 1 , T 2 , T 3 }, can be expressed as

G = { V = {T 1 , T 2 , T 3 }, E = {(T 1 , T 2 ), (T 1 , T 3 ), (T 2 , T 3 )}, F A = {((T 1 , T circle ), S T1 ), ((T 2 , T corner ), S T2 ), ((T 3 , T extremity ), S T3 )}, F E = {((T 1 , T 2 ), ℜ 1,2 ), ((T 1 , T 3 ), ℜ 1,3 ), ((T 2 , T 3 ), ℜ 2,3 )}}.
Since this forms a complete graph, it is obvious that there exist r = t(t-1) 2 edges for t attribute types. However, because of the use of fixed and completely labelled attributes, we can avoid the NP-hardness of the matching problem (cf.

Section 3).

In what follows, we explain how edges are labelled with spatial relations, computed between the vertices, and how vertices are labelled with shape features. This results in symbols being represented by ARGs like the one depicted in Fig. 2.

Edge Signatures using Spatial Relations. The choice of spatial relation models depends on configuration of the studied objects. If the objects are far enough from each other, their relations can be approximated by their centres and their discretised angle [START_REF] Miyajima | Spatial Organization in 2D Segmented Images: Representation and Recognition of Primitive Spatial Relations[END_REF]. Otherwise, if they are neither too far nor too close, relations can be approximated by their minimum bounding rectangle (MBR) [START_REF] Papadias | Spatial relations, minimum bounding rectangles, and spatial data structures[END_REF] as long as the objects are regular in shape (and their respective MBR are nonoverlapping). Approaches like angle histograms [START_REF] Wang | Human-Based Spatial Relationship Generalization Through Neural/Fuzzy Approaches[END_REF], tend to be more capable of dealing with overlapping. However, since they consider all pixels, their computational cost increases dramatically. Our work is inspired by the concept of fuzzy relations that take degrees of truth which is more natural than using standard, all-or-none relations [START_REF] Freeman | The modelling of spatial relations[END_REF].

Our model can be explained as follows. Given two shapes A and B for which we want to compute the relative position, and given a reference point R p , we cover the surrounding space at regular radial intervals of Θ = 2π/m by using a radial line. We compute this unique reference point R p = C A +C B 2 from the centroids C A and C B of the shapes under consideration. The uniqueness of R p thus avoids possible errors related to the choice of either of them as a reference. As depicted in Fig. 3, the line rotates over a cycle, and intersecting with object

X (A or B), generates a boolean histogram H, H(X, R p ) = [I(R p , jΘ)] j=0,...,m-1 , where I(R p , θ i ) = 1 if line(R p , θ i ) ∩ X = ∅ 0 otherwise.
This is extended wlog to the sector defined by two successive angle values (θ i and θ i+1 ) and is normalised with respect to the total area of the studied object such that H(.) = 1. This histogram provides both spatial and structural information. Fig. 3 provides an illustration for an arbitrary object X.

The histogram can be made rotationally invariant by projecting the smallest angle made by R p and C X on the horizontal axis. Also, translation does not affect at all, since it uses C X . Scaling does not produce any difference in H as it is normalised.

Vertex Signature via Shape Features. Each vertex has a distinct vocabulary type containing different shape and size information. Since spatial relations only encode relative positioning and point distributions, and do not completely exploit global shape information in the way shape descriptors do, we study the pertinence of R-signature [START_REF] Tabbone | A new shape descriptor defined on the radon transform[END_REF], region based Zernike moments (ZM) [START_REF] Kim | A region-based shape descriptor using zernike moments[END_REF], generic fourier descriptors (GFD) [START_REF] Zhang | Shape-based image retrieval using generic fourier descriptor[END_REF] and shape context (SC) [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF].

Symbol Recognition

Based on our symbol description, matching of two symbols is done by matching their corresponding ARGs. Consider two ARGs, query G q = (V q , E q , F q A , F q E ) and database

G d = (V d , E d , F d A , F d E )
, where the set of vertices V = {T 1 , . . . , T t }, and the set of edges E = {E 1 , . . . , E r }.

Our matching is straightforward i.e., matching has been made between the candidates only having the exact same set of vertices as well as exact labels.

To generalise this, we define a binary indicator function τ V A : Σ T → {0, 1} to check the presence of vertices in the ARG, where the value of τ V A (T) is 1 if T is present in V and 0, otherwise. For example, for the symbol shown in Fig. 2, the indicator τ V A = [0, 1, 1, 1] refers to the absence of thick components and the presence of circle, corner and extremity. Now, we can then set up bijective matching functions: σ : V q → V d and ϕ : E q → E d , respectively for vertices and edges. The fusion of both alignments provides the distance between two matched graphs G q and G d ,

Dist.(G q , G d ) = α t∈V δ(F q A (t), F d A (σ(t))) + (1 -α) r∈E δ(F q E (r), F d E (ϕ(r))),
where α ∈ [0, 1] and δ(a, b)

= L l=1 ||a l -b l || 2 .
The parameter α provides weight while matching.

α = 0: only vertex signature; -α = 1: only edge signature; and α = 0.5: equal weights to both vertex and edge signature.

In our experiments, we provide equal weights for both relations and shape distribution of the whenever they are integrated together.

Experiments

Dataset, Ground-truth and Evaluation Metric

Fig. 1 gives an overview of the dataset we are using for our experiments. The global dataset is composed of roughly 500 different known symbols, some of which come from [START_REF] Tooley | Aircraft electrical and electronic systems: principles, operation and maintenance. Aircraft engineering principles and practice[END_REF]. It shows that symbols may either be very similar in shape -and only differ by slight details -or either be completely different from a visual point of view. Symbols may also be composed of other known and significant symbols and need not necessary be connected.

Since there is no absolute ground-truth associated to our dataset, we have asked 6 volunteers to manually select what they consider as "similar" symbols, and we have merged their inputs to reduce possible subjective bias. They have chosen the candidates which have similar visual overall appearance or contain significantly similar parts with respect to the chosen query. In our testing protocol, we consider that a result returned from an algorithm is correct if at least one evaluator has selected the same result among similar items. In more formal terms, for each query the "ground-truth" is considered to be the set of symbols formed by the union of all human selected sets.

For every query, we rank the symbols at the output based on distance measure as described in Section 3. Since the number of similar symbols, according to the ground-truths, may vary a lot from one query to another, we use retrieval efficiency [START_REF] Kankanhalli | Cluster-based color matching for image retrieval[END_REF] as a measure for retrieval quality. For a chosen query and for a fixed number of K returned symbols, it can be expressed as,

η K = n/N if N ≤ K n/K otherwise,
where n is the number of returned relevant symbols and N the total number of relevant symbols.

Results and Discussions

We perform a series of tests, focussing on three major issues one after another.

1. Choice of optimal resolution for radial line model (as edge signature) and comparison with other relation models; 2. Choice of the best shape descriptors on vocabulary (as vertex signature) and comparison with global symbol shape descriptors; 3. Integration of both vertex and edge signatures.

In what follows, we compute the average retrieval efficiency for 30 queries, and vary K from 1 to 10.

Test 1. We consider the influence of different resolutions Θ in our edge signature. Its value represents the trade-off between the optimal choice of resolution -and thus precision of spatio-structural information capture -and time/space requirements. Without surprise, the lower Θ, the better the results. Following results in Fig. 4 (a), and given the relatively low gain of efficiency between 3 • and 1 • , we adopt 3 • for the rest of our experiments. We then compare our approach with state-of-the-art spatial relation models: cone-shaped [START_REF] Miyajima | Spatial Organization in 2D Segmented Images: Representation and Recognition of Primitive Spatial Relations[END_REF], angle histogram [START_REF] Wang | Human-Based Spatial Relationship Generalization Through Neural/Fuzzy Approaches[END_REF] and MBR [START_REF] Papadias | Spatial relations, minimum bounding rectangles, and spatial data structures[END_REF] as shown in Fig. 4 (b). Compared to the best performing model (MBR), our radial-line approach increases performance with a substantial difference. These results were already reported in [START_REF] Santosh | Symbol recognition using spatial relations[END_REF].

Test 2. In section 2 we already listed the set of global shape-descriptors we consider state-of-the-art: R-signature [START_REF] Tabbone | A new shape descriptor defined on the radon transform[END_REF], ZM [START_REF] Kim | A region-based shape descriptor using zernike moments[END_REF], GFD [START_REF] Zhang | Shape-based image retrieval using generic fourier descriptor[END_REF] and SC [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF]. For those descriptors, it is important to fit the best parameters. For GFD, we have tuned the radial and angular frequency parameters to achieve the best performance. For SC, we attempt to follow the indications given in [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] but they are restricted by the number of sample points of some symbols (i.e., we have images ranging from a few tens of pixels to thousands of pixels). In the case of ZM, we have used 36 Zernike functions of order less than or equal to 7. For Radon, projecting range is [0, π[.

We first employ them as vertex signatures only and then confront them with the same shape descriptors, applied to the overall shape. This comparison is illustrated in Fig. 5 (a) and (b), respectively. In this test, SC, R-signature and ZM are lagging behind GFD. Therefore, we consider GFD to be the best performing descriptor and will use it as a benchmark for further experiments. Test 3. We integrate both edge and vertex signatures in the ARG described in Section 2. However, integrating signatures on all vertices is not necessary since our spatial relation signature already carries sufficient information for those vocabulary types less influenced by shape variations (corner and extremity, for instance). Therefore, we have examined the use of signature from some specific types only, and on some combinations like, thick, thick and circle, thick, circle and corner, and so on. It is important to notice that we have never left out thick components from all tested combinations because those vocabulary elements provide the largest amplitude in shape and size variations. To handle this, we use all aforementioned global signal based descriptors as vertex signatures in order to see how well they improve retrieval efficiency. In our tests, substantial advancement is achieved from the combination of vertices, labelled with thick and circle vocabulary types. Fig. 6 shows results from using different shape signatures combined with the edge signature using a radial line model at 3 o resolution. GFD provides the best results. Surprisingly, not all vertex shape descriptors improve upon the plain edge signatures, but some even decrease in performance. R-signature provides an example.

Comparison. In Fig. 7, a comparison is made among the best of all experiments: MBR from the basic spatial relation model and GFD from shape descriptors. Our method outperforms all with a difference of more than 16%. For a few queries, an illustration is shown in Fig. 8.

Another important issue accounting for retrieval performance is execution time. In Table 1, we provide it based on comparison of the methods in Fig. 7. This proves that the increase in recogition efficiency comes at a cost in execution time. This cost remains in a reasonable limit, however. We have presented an ARG based symbol description method, where relational signatures formalise all possible connections between the vocabulary types which are labelled with global shape signatures. Our method has proven to significantly outperform state-of-the-art basic spatial relation models as well as global signal based descriptors.

Fig. 1 .

 1 Fig. 1. A few symbols in FRESH dataset containing both linear as well as composite symbols.

ARGFig. 2 .

 2 Fig.2. ARG description for a symbol from its corresponding visual parts -an example.

Fig. 3 .

 3 Fig. 3. Computing spatial relations using radial line rotation .

Fig. 4 .

 4 Fig.4. Average retrieval efficiency using spatial relations.

Fig. 5 .

 5 Fig. 5. Average retrieval efficiency using global signal based shape descriptors.

Fig. 6 .

 6 Fig.6. Average retrieval efficiency while integrating vertex signature in our ARG framework.

Fig. 7 .

 7 Fig. 7. Comparison among the best of all categories.

Fig. 8 .

 8 Fig. 8. Visual illustration of symbol ranking at the output for a few queries: Q1, Q2 and Q3, showing for true retrieval and false, otherwise. The first symbol on the top always corresponds to the chosen query. Symbols are ranked from top to bottom based on decreasing order of similarity.

Table 1 .

 1 Edge sign. + MBR Edge sign. GFD GFD vertex sign. Average running time (sec.) for matching for a single pair.

					Test 3
	Test 1 Test 2 Time 02 04 09		14
							Our method
	MBR	GFD		Edge signature	Vertex signature	Integrating both
	Q1 Q2	Q3 Q1 Q2	Q3	Q1	Q2	Q3	Q1	Q2	Q3	Q1	Q2	Q3
	1.
	2.
	3.
	4.
	5.
	6.
	7.
	8.
	9.
	10.

Compared to our previous work[START_REF] Santosh | Symbol recognition using spatial relations[END_REF] where the 1D (vanilla) version of the R-signature has been employed, the complete version reported here provides better performance, without altering the global validity of our findings, however.