W Li

J S R Nielsen

S Puchinger

V Sidorenko

Solving Shift Register Problems over Skew Polynomial Rings using Module Minimisation

Keywords: Skew Polynomials, Ore Polynomials, Shift Register Synthesis, Module Minimisation, Gabidulin Codes

For many algebraic codes the main part of decoding can be reduced to a shift register synthesis problem. In this paper we present an approach for solving generalised shift register problems over skew polynomial rings which occur in error and erasure decoding of -Interleaved Gabidulin codes. The algorithm is based on module minimisation and has time complexity O( µ 2 ) where µ measures the size of the input problem.

Introduction

Numerous recent publications have dealt with shaping the core of various decoding algorithms for Reed-Solomon (RS) and other codes around F q [x] module minimisation, lattice basis reduction or module Gröbner basis computation: three computational concepts which all converge to the same in this instance. First for the Guruswami-Sudan list decoder [START_REF] Alekhnovich | Linear Diophantine equations over polynomials and soft decoding of Reed-Solomon codes[END_REF][START_REF] Cohn | Ideal forms of Coppersmith's theorem and Guruswami-Sudan list decoding[END_REF][START_REF] Lee | List decoding of Reed-Solomon codes from a Gröbner basis perspective[END_REF], then for Power decoding [START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF] and also either type of decoder for Hermitian codes [START_REF] Nielsen | Sub-quadratic decoding of one-point Hermitian codes[END_REF].

The impact of this can be said to be two-fold: firstly, by factoring out coding theory from the core problem, we enable the immediate use of sophisticated algorithms developed by the computer algebra community such as [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF][START_REF] Zhou | Efficient algorithms for order basis computation[END_REF]. Secondly, the setup has proved very flexible and readily applicable in settings which were not envisioned to begin with, such as the aforementioned Power decoder for Hermitian codes, or recently for Power decoding of RS codes up to the Johnson bound [START_REF] Nielsen | Power decoding Reed-Solomon codes up to the Johnson radius[END_REF].

The main goal of this paper is to extend the module minimisation description to skew polynomial rings and Gabidulin codes, in particular Interleaved Gabidulin codes, with the aim of enjoying similar benefits. Concretely, we lay a foundation by extending the core terms of weak Popov form and orthogonality defect, as well as extending the elegantly simple Mulders-Storjohann algorithm [START_REF] Mulders | On lattice reduction for polynomial matrices[END_REF] to matrices over skew polynomial rings. We analyse its complexity when applied to the shift register problem which arise when decoding Interleaved Gabidulin codes. Finally, we extend the Demand-Driven algorithm for F q [x] shift register problems [START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF], which is derived from the Mulders-Storjohann, also to the skew polynomial setting.

Gabidulin codes [START_REF] Delsarte | Bilinear forms over a finite field, with applications to coding theory[END_REF][START_REF] Gabidulin | Theory of codes with maximum rank distance[END_REF][START_REF] Roth | Maximum-rank array codes and their application to crisscross error correction[END_REF] are maximum rank distance codes with various applications like random linear network coding [START_REF] Koetter | Coding for errors and erasures in random network coding[END_REF][START_REF] Silva | A rank-metric approach to error control in random network coding[END_REF] and cryptography [START_REF] Gabidulin | Ideals over a non-commutative ring and their application in cryptology[END_REF]. They are the rank-metric analogue of RS codes. An Interleaved Gabidulin code [START_REF] Loidreau | Decoding rank errors beyond the error correcting capability[END_REF][START_REF] Sidorenko | Skew-feedback shift-register synthesis and decoding interleaved Gabidulin codes[END_REF][START_REF] Silva | A rank-metric approach to error control in random network coding[END_REF] is a direct sum of several (n, k i ) Gabidulin codes: these can be decoded in a collaborative manner, improving the error-correction capability beyond the usual half the minimum rank distance of Gabidulin codes. Similar to Interleaved RS codes, see [START_REF] Schmidt | Collaborative decoding of interleaved Reed-Solomon codes and concatenated code designs[END_REF] and its references, the core task of decoding can be reduced to what is known as a multi-sequence skew-feedback shift register synthesis problem [START_REF] Sidorenko | Skew-feedback shift-register synthesis and decoding interleaved Gabidulin codes[END_REF].

In this paper, we use the introduced module minimisation description to solve a more general form of this problem, which we abbreviate MgLSSR: Problem 1 (MgLSSR). Given skew polynomials s i , g i and non-negative integers γ i ∈ N 0 for i = 1, . . . , , find skew polynomials λ, ω 1 , . . . , ω , with λ of minimal degree such that the following holds:

λs i ≡ ω i mod g i (1
)

deg ω i + γ i < deg λ + γ 0 (2)
The original problem of [START_REF] Sidorenko | Skew-feedback shift-register synthesis and decoding interleaved Gabidulin codes[END_REF] set g i to powers of x and γ i = 0. The above is a natural generalisation, which covers error and erasure decoding of Gabidulin codes [START_REF] Li | On transform-domain error and erasure correction by Gabidulin codes[END_REF], as well as an Interleaved extension of the Gao-type decoder for Gabidulin codes ([27, §3.2] combined with the ideas of [START_REF] Sidorenko | Skew-feedback shift-register synthesis and decoding interleaved Gabidulin codes[END_REF]). For cases where the algorithm of [START_REF] Sidorenko | Skew-feedback shift-register synthesis and decoding interleaved Gabidulin codes[END_REF] applies, the Demand-Driven algorithm we present has the same complexity. However, the more general perspective of module minimisation gives conceptually simpler proofs, and may prove useful for gaining further insights or faster, more sophisticated algorithms. Normal form computation of matrices over skew rings and Ore rings has been investigated before, e.g. [START_REF] Abramov | On solutions of linear functional systems[END_REF][START_REF] Beckermann | Fraction-free row reduction of matrices of Ore polynomials[END_REF], but the focus has been over rings such as Z or K[z] for some field K, where coefficient growth is important to control. Since we are inspired mainly by the application to Gabidulin codes, where the skew ring is over a finite field, we count only operations performed in the field; in this measure those previous algorithms are much slower than what is presented here.

We set basic notation in Section 2. Section 3 describes how to solve Problem 1 using module minimisation, and gives the Mulders-Storjohann algorithm for skew polynomial modules to accomplish this. We introduce important concepts for arguing about such modules in Section 4 for performing a complexity analysis. Section 5 describes how to then derive the faster Demand-Driven algorithm. Due to lack of space, a number of proofs are omitted.

Notation and Remarks on Generality

Let K be a field. Denote by R = K[x; θ, δ] the noncommutative ring of skew polynomials over K with automorphism θ and derivation δ. Being an Ore extension, R is both a left and right Euclidean ring. See [START_REF] Ore | Theory of non-commutative polynomials[END_REF] for more details.

For coding theory we usually take K as a finite field F = F q r for a prime power q and θ as the Frobenius automorphism θ(a) = a q for a ∈ F q r . Also, non-vanishing derivations δ are usually not considered, a notable exception being [START_REF] Boucher | Linear codes using skew polynomials with automorphisms and derivations[END_REF]. The algorithms in this paper are correct for any field, automorphism and derivation. For complexities, we are counting field operations, and we often assume δ = 0.

By a ≡ b mod c we denote the right modulo operation in R, i.e., that there exists d ∈ R such that a = b + dc. By "modules" we will mean left R-modules. We extensively deal with vectors and matrices over R. Matrices are named by capital letters (e.g. V ). The ith row of V is denoted by v i and the jth element of a vector v is v j . v ij is the (i, j)th entry of a matrix V . Indices start at 0.

-The degree of a vector v is deg v := max i {deg v i } (and deg 0 = -∞) and the degree of a matrix V is deg

V := i {deg v i }. -The max-degree of V is maxdeg V := max i {deg v i } = max i,j {deg v ij }. -The leading position of a vector v is LP(v) := max{i : deg v i = deg v}. Furthermore LT(v) := v LP(v) and LC(v) is the leading coefficient of LT(v).
3 Finding a Solution using Module Minimisation

In the sequel we consider a particular instance of Problem 1, so R, ∈ N, and s i , g i ∈ R, γ i ∈ N 0 for i = 1, . . . , are arbitrary but fixed. We assume deg s i ≤ deg g i for all i since taking s i := s i mod g i yields the same solutions. Denote by M the set of all vectors v ∈ R +1 satisfying (1), i.e.,

M := (λ, ω 1 , . . . , ω ) ∈ R +1 | λs i ≡ ω i mod g i ∀i = 1, . . . , }. (3)
Lemma 1. M with component-wise addition and left multiplication by elements of R forms a left module over R. The rows of M form a basis of M:

M =       1 s 1 s 2 . . . s 0 g 1 0 . . . 0 0 0 g 2 . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . g      
The above gives a simple description of all solutions of the congruence relations (1). To solve Problem 1, we therefore need an element in the M which satisfies the degree condition (2) and has minimal degree. For this purpose, define

Φ : R +1 → R +1 , u = (u 0 , . . . , u ) → (u 0 x γ0 , . . . , u x γ ). (4)
We can extend the domain of Φ to matrices over R by applying it row-wise.

It is easy to see that Φ(M) is also a left R-module and that Φ is a module isomorphism. Using this notation, we can restate how to solve Problem 1: Proof: v ∈ M * is a solution to Problem 1 iff it satisfies (2) and v 0 has minimum possible degree. That v satisfies (2) means deg v 0 +γ 0 > deg v i +γ i and so deg(v 0 x γ0 ) > deg(v i x γi ) i.e. LP(Φ(v)) = 0. The reverse direction is similar.

Lemma 2. A vector v ∈ M *
So we should find a vector v ∈ Φ(M) with minimum-degree leading term among vectors with leading position zero. We do this by finding a basis of Φ(M) of a specific form. This extends similar ideas for matrices over K[x] [START_REF] Mulders | On lattice reduction for polynomial matrices[END_REF][START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF]. Definition 1. A matrix V over R is in weak Popov form if the leading positions of all its non-zero rows are different.

The following value function for R vectors will prove useful: ψ : R +1 → N 0 ,

ψ(v) = ( + 1) deg v + LP(v) + 1
for v = 0 and ψ(0) = 0.

Lemma 3. Let V be a matrix in weak Popov form whose rows are a basis of a left R-module V. Then every u ∈ V * satisfies deg u ≥ deg v, where v is the row of V with LP(v) = LP(u).

Proof: Let u ∈ V * , and so ∃a 0 , . . . , a ∈ R s.t. u = i=0 a i v i . The u i all have different leading position, so the a i v i must as well for those a i = 0, which in turn means that the their ψ(a i v i ) are all different. Notice that for any two u 1 , u 2 with ψ(u 1 ) = ψ(u 2 ), then ψ(u 1 + u 2 ) either equals ψ(u 1 ) or ψ(u 2 ). Applied inductively, that implies that there is an i such that ψ(u) = ψ(a i v i ), which gives

LP(u) = LP(v i ) and deg u = deg a i + deg v i .
Lemma 2 and Lemma 3 imply that a basis of M in weak Popov form gives a solution to Problem 1 as one of its rows. The following definition leads to a remarkably simple algorithm for computing such a basis: Algorithm 1, an R variant of the Mulders-Storjohann algorithm [START_REF] Mulders | On lattice reduction for polynomial matrices[END_REF], originally described for K[x].

Definition 2. Applying a simple transformation on a matrix V means finding non-zero rows v i , v j , i = j such that LP(v i ) = LP(v j ) and deg v i ≤ deg v j , and replace v j by v j -αx β v i , where β = deg v j -deg v i and α = LC(v j )/θ β (LC(v i )).

Remark 1. Note that a simple transformations cancels the leading term of the polynomial LT(v j ). Also elementary row operations keep the module spanned by the matrix' rows unchanged, see e.g. [START_REF] Beckermann | Fraction-free row reduction of matrices of Ore polynomials[END_REF], so the same is true for any sequence of simple transformations. Proof: By Lemma 4, the value of one row of V decreases for each simple transformation. The sum of the values of the rows must at all times be nonnegative so the algorithm must terminate. Finally, when the algorithm terminates there are no simple transformations possible on V anymore, i.e. there are no i = j such that LP(v i ) = LP(v j ). That is to say, V is in weak Popov form.

This gives an algorithm to solve Problem 1. The above proof could also easily lead to a rough complexity estimate. To obtain a more fine-grained one, we will in the next section restrict ourselves to matrices which are square and full rank.

Complexity Analysis

Lenstra [START_REF] Lenstra | Factoring multivariate polynomials over finite fields[END_REF] introduced the notion of orthogonality defect of square, full rank K[x] matrices, and in [START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF], it was shown it can describe the complexity of the Mulders-Storjohann and Alekhnovich [START_REF] Alekhnovich | Linear Diophantine equations over polynomials and soft decoding of Reed-Solomon codes[END_REF] algorithms for such matrices more finegrained than originally, and that this improves the asymptotic estimate when the input comes from shift register problems. The same concept cannot immediately be carried over to R matrices, since it is defined using the determinant. For noncommutative rings, there are no functions behaving exactly like the classical determinant, but the Dieudonné determinant [START_REF] Dieudonné | Les déterminants sur un corps non commutatif[END_REF] shares sufficiently many properties with it for our use. Simply defining this determinant requires us to pass to the field of fractions of R.

Dieudonné Determinant and Orthogonality Defect

The following algebra is standard for noncommutative rings, so we will go through it quickly; more details can be found in [6, Chapter 1]. We know that R is a principal left ideal domain which implies that it is left Ore and therefore has a unique left field of fractions Q = {s -1 r : r ∈ R, s ∈ R * }/(∼), where ∼ is the congruence relation s -1 r ∼ s -1 r if ∃u, u ∈ R * such that ur = u r and us = u s . The degree map on R can be naturally extended to Q by defining

deg : Q → Z ∪ {-∞}, s -1 r → deg r -deg s. Let [Q * , Q * ] be the commutator of Q * , i.e. the multiplicative group generated by {a -1 b -1 ab : a, b ∈ Q * }. Then Q ab = Q * /[Q * , Q * ] is an abelian group called the multiplicative abelianization of Q * . There is a canonical homomorphism φ : Q * → Q ab , x → x • [Q * , Q * ]. Since the elements (a -1 b -1 ab) ∈ [Q * , Q * ] have degree deg(a -1 b -1 ab) = deg(ab)- deg(ba) = 0, we can pass deg through φ in a well-defined manner: deg φ(x) = deg x for all x ∈ Q * .
The following lemma was proved by Dieudonné [START_REF] Dieudonné | Les déterminants sur un corps non commutatif[END_REF] and can also be found in [START_REF] Draxl | Skew Fields[END_REF]. Lemma 5. There is a function det :

Q n×n → Q ab s.t. for all A ∈ Q n×n , k ∈ Q:
For Line 6, we can compute the needed coefficient α in complexity O(µρ): if deg 2 g h > 1 2 deg g h , we simply compute the entire polynomial λ 0 sh mod gh in time O(µ 2 ). Otherwise, an easy argument shows that at most #supp(g h ) + 1 coefficients of λ j sh affects the computation of α. Each of these can be computed by convolution in time O(µ).

For generic g i , Algorithm 2 will have complexity O( µ 3 ) which is usually worse than O( 2 µ 2 ) of Algorithm 1. However, for decoding of Interleaved Gabidulin codes, two important cases are g i = x k (syndrome decoding [START_REF] Sidorenko | Skew-feedback shift-register synthesis and decoding interleaved Gabidulin codes[END_REF]) and g i = x q m -1 (Gao-type decoding [27, §3.2]), and here Algorithm 2 runs in complexity O( µ 2 ).

Remark 2. Algorithm 2 bears a striking similarity to the Berlekamp-Massey variant for multiple shift registers [START_REF] Sidorenko | Skew-feedback shift-register synthesis and decoding interleaved Gabidulin codes[END_REF] where all g i are powers of x, and has the same running time in this case. However, using the language of modules, we obtain a more general algorithm with a conceptually simpler proof, and we can much more readily realise algebraic properties of the algorithm. For instance, using known properties for the weak Popov form, it is trivial to prove that Algorithm 2 can be modified to return a basis for all solutions to the shift register problem, as well as decompose any given solution as an R-linear combination of this basis.

Conclusion

In this paper, we have given two module-based methods for solving generalised shift register problems over skew polynomial rings. For ordinary polynomial rings, module minimisation has proven a useful strategy for obtaining numerous flexible, efficient while conceptually simple decoding algorithms for Reed-Solomon and other code families. Our results introduce the methodology and tools aimed at bringing similar benefits to Gabidulin, Interleaved Gabidulin and other skew polynomial-based codes.

  is a solution to Problem 1 if and only if LP (Φ(v)) = 0 and for all u ∈ M * with LP (Φ(u)) = 0 it holds that deg Φ(v) ≤ deg Φ(u).

Lemma 4 .

 4 If v replaces v in a simple transformation, then ψ(v ) < ψ(v). Proof: The operations used in a simple transformation ensure that deg v ≤ deg v. If deg v < deg v, we are done because LP(v ) < + 1. If deg v = deg v, then LP(v ) < LP(v): by the definition of the leading position, all terms to the right of LP(v) in v and αx β v i , and therefore also in v , have degree less than deg v. Furthermore deg v LP(v) < deg v by the definition of a simple transformation. Algorithm 1 Mulders-Storjohann for R matrices Theorem 1. Algorithm 1 is correct.

Acknowledgements Sven Puchinger (grant BO 867/29-3), Wenhui Li and Vladimir Sidorenko (both grant BO 867/34-1) were supported by the German Research Foundation "Deutsche Forschungsgemeinschaft" (DFG). Johan S. R. Nielsen gratefully acknowledges the support of the Digiteo foundation, project IdealCodes.

(i) det I = 1, where I is the identity matrix in Q n×n . (ii) If A is obtained from A by an elementary row operation, then det A = det A. (iii) If A is obtained from A by multiplying a row with k, then det A = φ(k) det A. Definition 3. A function det with the properties of Lemma 5 is called a Dieudonné determinant.

Note that contrary to the classical determinant, a Dieudonné determinant is generally not unique. For the remainder of the paper, consider det to be any given Dieudonné determinant. Lemma 6. Let A ∈ Q n×n be in triangular form with non-zero diagonal elements d 0 , . . . , d n-1 . Then det A = n-1 i=0 φ(d i ).

Proof: Since d i = 0 for all i, we can multiply the ith row of A by d -1 i and get a unipotent triangular matrix A . Any unipotent triangular matrix can be obtained by elementary row operations from the identity matrix I. Thus

Clearly, the notion of weak Popov form generalises readily to matrices over Q. We will now examine how this notion interacts with the Dieudonné determinant and introduce the concept of orthogonality defect. The statements in this section are all Q variants of the corresponding statements for K[x] matrices, see [START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF].

proof sketch: We can assume that LP(v i ) = i for all i because if not, we can change the order of the rows of V and obtain a matrix with the same determinant and degree. We can then apply elementary row operations to bring the matrix to upper triangular form. After these row operations, the property LP(v i ) = i is preserved and deg v ii is equal to deg v ii of the start matrix for all i. By Lemma 6 the degree of the determinant equals the sum of the degree of the diagonal elements, and hence deg V = deg det V .

Complexity of Mulders-Storjohann

We can now bound the complexity of Algorithm 1 using arguments similar to those in [START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF]. These are in turn, the original arguments of [START_REF] Mulders | On lattice reduction for polynomial matrices[END_REF] but finer grained by using the orthogonality defect. In the following, let µ := max i {γ i + deg g i }.

We can assume that γ 0 < µ since otherwise (1, s 1 , . . . , s ) is the minimal solution to the MgLSSR.

Theorem 2. Over R with derivation zero, Algorithm 1 with input matrix Φ(M ) performs at most ( + 1)(µ -γ 0 + 1) simple transformations and performs O( 2 µ 2 ) operations over K.

Proof: Every simple transformation reduces the value ψ of one row with at least 1. So the number of possible simple transformations is upper bounded by the difference of the sums of the values of the input matrix Φ(M ) and the output matrix V , i.e.:

where the last equality follows from deg

One simple transformation consists of calculating v j -αx β v i , so for every coefficient in v i , we must apply θ β , multiply by α and then add it to a coefficient in v j , each being in O(1). Since deg v j ≤ µ this costs O( µ).

Demand-Driven Algorithm

It was observed in [START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF] that the Mulders-Storjohann algorithm over K[x] admits a "demand-driven" variant when applied to matrices coming from shift register problems, where coefficients of the working matrix are computed only when they are needed. This means a much lower memory requirement, as well as a better complexity under certain conditions. Over R, Algorithm 1 admits exactly the same speedup; in fact, both the algorithm and the proof are almost line-for-line the same for R as for K[x]. We therefore focus on the idea of the algorithm, and the original proofs can be found in [START_REF] Nielsen | Generalised multi-sequence shift-register synthesis using module minimisation[END_REF] (extended version).

The central observation is that due to the special form of M of Lemma 1, only the first column is needed during the Mulders-Storjohann algorithm in order to construct the rest. That is formalised in the following lemma: Lemma 9. Consider Algorithm 1 with input Φ(M ). Consider a variant where, when replacing v j with v j in a simple transformation, instead replace it with v j = (v j,0 , v j,1 mod g1 , . . . , v j, mod g ). This does not change correctness of the algorithm or the upper bound on the number of simple transformations performed.

The Demand-Driven algorithm, Algorithm 2, therefore calculates just the first element of a vector whenever doing a simple transformation, being essentially enough information. To retain speed it is important, however, that the algorithm can also figure out which simple transformation it can next apply, without having to recompute the whole matrix. For this, we cache for each row its degree η j and the leading coefficient of its leading position α j . The following observations then lead to Algorithm 2:

1. In Φ(M ) there is at most one possible choice of the first simple transformation, due to the matrix' shape. This is true throughout the algorithm, making it deterministic.

Algorithm 2 Demand-Driven algorithm for MgLSSR

Input: sj ← s1,jx γ j , gj ← gjx γ j for j = 1, . . . , Output: The first column of a basis of M whose Φ image is in weak Popov form. 1 (η, h) ← (deg, LP) of (x γ 0 , s1, . . . , sσ) 2 if h = 0 then return (1, 0, . . . , 0) 3 (λ0, . . . , λ ) ← (x γ 0 , 0, . . . , 0) 4 αjx η j ← the leading monomial of gj for j = 1, . . . ,

11 return λ0x -η 0 , . . . , λ x -η 0 2. To begin with, if there is a possible simple transformation, row 0 is involved. Just before doing a simple transformation, we possibly swap the two rows involved such that the row changed is always row 0. That means row 0 is always involved if there is a possible simple transformation, and that the algorithm terminates when row 0 has leading position 0. 3. To begin with row i has leading position i for i > 0. The above swap ensures that this will keep being true. 4. After doing a simple transformation, we need to update the degree, leading position and leading coefficient of only row 0; the rest remains unchanged. We do this by going through each possible degree and leading position in decreasing order of value ψ. This is correct since we know that the simple transformation must decrease the value of row 0.

To express the complexity, by supp(f ), f ∈ R, we mean the set of degrees such that f has a non-zero coefficient for this degree. By deg 2 f we mean the degree of the second largest coefficient. Let again µ := max i {γ i + deg g i }. proof sketch: We only prove the complexity statement. Clearly, all steps of the algorithm are essentially free except Line 6 and Line 9. Observe that every iteration of the while-loop decrease an estimate on the value of row 0, whether we enter the if-branch in Line 7 or not. So by the arguments of the proof of Theorem 2, the loop will iterate at most O( µ) times. Each execution of Line 9 costs O(µ) since the λ j all have degree at most µ.


Nicolas Beldiceanu
email: beldiceanu@cosytec.com

Eric Bourreau
email: bourreau@cosytec.com

H Simonis
email: simonis@cosytec.com

A Note on Perfect Square Placement

  L'archive ouverte pluridisciplinaire

Problem

We discuss the following perfect square placement problem (also called squared square problem [START_REF] Croft | Guy Unsolved Problems in Geometry[END_REF][LW92]): Pack a set of squares with given integer sizes into a bigger square in such a way that no squares overlap each other and all square borders are parallel to the border of the big square. For a perfect placement problem, all squares have different sizes. The sum of the square surfaces is equal to the surface of the packing square, so that there is no spare capacity. The problem data come from [START_REF] Bouwkamp | Duijvestijn Catalogue of Simple Perfect Squared Squares of orders 21 through 25[END_REF], which lists all simple perfect square placement problems of order 21 to 25. Simple in this context means that the problems do not contain a subset of squares (greater than one) which are placed in a rectangle. Duijvestijn has shown in 1962 that no instances exist with less than 21 squares [START_REF]Duijvestijn Simple Perfect Squared Square of Lowest Order[END_REF]. A single problem exists with 21 squares. In this problem, the squares have sizes [2,4,6,7,8,9,11,15,16,17,18,19,24,25,27,29,33,35,37,42,50] and must be packed into a square of size 112. A solution is shown below. There are 8 problems of order 22, 12 of order 23, 26 of order 24, and 160 of order 25. A constraint approach for this problem was already studied in [START_REF] Aggoun | Beldiceanu Extending CHIP in Order to Solve Complex Scheduling Problems[END_REF], where a solution is found in 38 seconds. The same problem was mentioned in [START_REF] Van Hentenryck | Scheduling and Packing in the Constraint Language cc(FD)[END_REF], where a solution was found in 60 seconds. Both results were obtained on SUN workstations. With the program shown here, we find a solution in 0.3 seconds on a PC.

Program

The CHIP program to solve the square placement problem is given below. Each square is defined by two domain variables X and Y and integer Width and Height values. The predicate run has an input parameter, which is the problem instance number and two output parameters, the lists of the X and Y values for the coordinates of the bottom left corners of the squares. The predicate gen_data/6 reads the problem definition of the instance and generates lists of domain variables for the X and Y values. The square sizes are given in decreasing order, the largest square is the first entry in the lists. The predicate gen_rect/5 creates a list of rectangles (each a list of four variables X, Y, Width and Height) from the four parameter lists. The diffn [START_REF] Beldiceanu | Contejean Introducing Global Constraints in[END_REF] constraint states that the set of squares do not overlap and fit into the available space. The first argument is a list of rectangles, generated by gen_rect/5. The fourth argument defines the available space. We also use two redundant cumulative [START_REF] Aggoun | Beldiceanu Extending CHIP in Order to Solve Complex Scheduling Problems[END_REF] constraints which project the placement in x or y dimension. The first argument gives the origin, the second argument the size and the third argument the height of the squares to be placed. Arguments 6 and 7 define the available placement space. The argument 8 enforces an additional constraint that the available space must be completely used. The labeling routine is described in [START_REF] Aggoun | Beldiceanu Extending CHIP in Order to Solve Complex Scheduling Problems[END_REF]. We first find an assignment of all X variables, and then find an assignment for the Y variables. Many other assignment strategies are possible, the one shown here gives the best results on perfect placement problems. The labeling method works by finding the smallest value in the domain of the unassigned variables and then assigning this value to some variable. On backtracking, other variables will be assigned to this value. The predicate lmindomain/2/3 finds the smallest value a list of domain variables. The predicate fix_min/3 is a non-deterministic predicate which either fixes a variable to the minimal value

Possible Improvements

A major decrease of the total execution time is possible by introducing constraints to avoid redundant, symmetrical solutions. For each problem in the data-set, there are 8 symmetrical solutions obtained by applying the symmetries in the plane. We choose not to attempt this in this benchmarking problem, in order to simplify comparisons between different systems.

Summary

In this note we have presented the solution of a set of perfect square placement problems with CHIP. A rather simple CHIP program can be used to find solutions to all problem instances.

Figure 1 :

 1 Figure 1: Solution of the 21 square problem

Square placement program in CHIP Results The

  table below gives the results for CHIP V5.2 on the benchmark set. The results were obtained with the commercial CHIP version on a Pentium MMX 233MHz PC, running WindowsNT 4.0 with 64 Mb of memory. The table shows the instance number, the number of backtracking steps and the execution time required to find a first solution, the total number of solutions found and the backtracking count and execution times to find all solutions by complete enumeration of the search space. Note that for most problems 8 solutions are found, these correspond to one solution under all plane symmetries. Problems 166 and 167 (also 168 and 169, and 182 and 183) have identical data, but two independent (non-isomorphic) solutions. These problems are given in this form in the original data set. The maximum time to find a first solution was spent on problem 59 with 32.5 seconds (3877 backtracking steps), the maximum time to explore the complete search space was spent in problem 48 with 503 seconds (66100 backtracking steps). On the other hand, for 38 problems a first solution was found without backtracking, and 94 problems were explored completely with less than 2000 backtracking steps.

	29 77 125	642 50 8 2	5127 681 321 261	8 8 8 8		3987 1611 2941 1138	27489 10635 24576 7491
	run(Nr,LX,LY):-30 78 126 174	572 71 31 2	4697 751 300 240	8 8 8 8		3066 3068 1491 1036	22082 18677 9443 6960
	gen_data(Nr,LX,LY,LW,LH,Max), 31 10 240 79 423 3495 127 6 311 175 1 240 gen_rect(LX,LY,LW,LH,Lrect), 32 12 270 80 12 250 128 0 221 176 0 230 Max1 is Max + 1, diffn(Lrect,unused,unused,[Max1,Max1]), 8 8 8 8 8 8 8 8 33 211 1993 8 81 18 300 8 129 34 431 8 177 1 251 8	2227 4804 2829 2718 3441 1172 7715 1511 3417 1765 5468 3121	14010 31896 16864 24055 23844 7681 70091 14290 18066 12488 30805 20570
	cumulative(LX,LW,LH,unused,unused,Max,Max1,[Max,0,1]), 34 0 141 8 2379 16224 82 9 250 8 2364 11737 130 121 1582 8 2726 25847 178 4 330 8 2306 15061
	cumulative(LY,LH,LW,unused,unused,Max,Max1,[Max,0,1]), 35 443 3655 8 2803 20319 83 86 831 8 2764 19758 131 4 230 8 1584 11016 179 65 751 8 1245 8473
	labeling(LX), 36 84 15 2 132 87 180 10	180 321 781 410	8 8 8 8		676 5494 1585 831	3295 45526 9453 6639
	labeling(LY). 37 85 35 1 133 0 181 185	190 551 230 1793	8 8 8 8		1312 2191 1606 983	8272 17094 10045 7521
	38 86 134 182 gen_data(Nr,LX,LY,Sizes,Sizes,MaxPlace) :-49 470 8 161 1993 8 11 290 8 0 230 16 39 9 250 8 87 1 240 8 135 18 361 8 183 0 211 16 data(Nr,NbSquare,MaxPlace,SizeSquares), length(LX,NbSquare), 40 8 271 8 88 3 230 8 136 0 190 8 184 0 290 8	2533 3426 1863 3233 1490 4998 1113 3233 3117 5154 5125 1128	16333 27961 10965 21020 9203 39787 6840 21341 21471 37103 41770 9714
	length(LY,NbSquare), 41 18 89 1052 137 21 185 0	311 10986 411 200	8 8 8 8		1044 4204 3141 804	6500 35491 23093 6679
	LX :: 1..MaxPlace, 42 28 90 90 138 6 186 0	460 931 231 241	8 8 8 8		1629 4565 1887 2037	10765 34680 13280 12659
	LY :: 1..MaxPlace, 43 0 91 0 139 67 187 2	250 200 841 240	8 8 8 8		818 4415 10246 583	5858 28882 93965 3324
	reverse(SizeSquares,Sizes). 44 0 201 92 12 250 140 2 220 188 5 271	8 8 8 8		1297 4017 1247 706	9274 23914 9223 5138
	45 93 141 189	3 61 0 75	230 711 200 1352	8 8 8 8		839 3070 1554 1744	5388 24966 11286 17155
	fix_min([V|R], M, R) :-V = M, inc_choice. fix_min([V|R], M, [V|S]) :-V #> M, fix_min(R, M, S). lmindomain([H|T], M) :-domain_info(H, Min, _, _, _, _), lmindomain(T, Min, M). lmindomain([], M, M). lmindomain([X|Y], M, Mend) :-domain_info(X, Min, _, _, _, _), M1 is min(M,Min), lmindomain(Y, M1, Mend). inc_choice. inc_choice:-incval(choice,_), fail. Program 1 : Instance Backtrack for first solution Time for first solution [ms] Number of solutions 46 4 190 8 94 4 250 8 142 95 1302 8 190 1 320 8 47 0 200 8 95 133 1792 8 143 0 260 8 191 0 240 8 48 81 891 8 96 69 441 8 144 6 260 8 192 2 250 8 1 22 341 8 49 1275 10455 8 97 130 1502 8 145 184 3094 8 193 275 3635 8 2 180 1242 8 50 64 722 8 98 0 180 8 146 0 220 8 194 40 471 8 3 53 511 8 51 98 982 8 99 0 170 8 147 7 251 8 195 10 341 8 4 28 250 8 52 359 3174 8 100 3479 30724 8 148 119 1101 8 196 8 281 8 5 94 701 8 53 24 440 8 101 684 7401 8 149 0 221 8 197 7 441 8 6 93 661 8 54 200 2002 8 102 60 821 8 150 123 1282 8 198 0 220 8 7 11 190 8 55 17 321 8 103 52 841 8 151 1 201 8 199 2 250 8 8 6 160 8 56 0 170 8 104 91 1001 8 152 47 651 8 200 4 250 8 9 0 140 8 57 94 801 8 105 0 210 8 153 6 281 8 201 5 281 8 10 14 231 8 58 54 731 8 106 16 340 8 154 25 401 8 202 0 230 8 11 177 1312 8 59 3877 32537 8 107 11 351 8 155 0 210 8 203 0 230 8 12 7 120 8 60 2660 23124 8 108 17 441 8 156 22 381 8 204 0 260 8 13 5 150 8 61 8 260 8 109 2 240 8 157 1 211 8 205 2 380 8 14 659 4206 8 62 152 901 8 110 28 471 8 158 49 691 8 206 0 270 8 15 2795 20509 8 16 818 6199 8 17 40 341 8 18 1 150 8 19 6 190 8 20 239 1503 8 21 9 230 8 22 895 4987 8 23 1862 15943 8 24 393 2884 8 25 50 431 8 26 1041 10245 8 27 70 471 8 28 153 1342 8 76 14 271 8 124 1 230 8 172 1 240 8 75 7 240 8 123 6 251 8 171 0 230 8 74 0 181 8 122 642 5308 8 170 0 230 8 73 10 230 8 121 2 241 8 169 4 291 16 72 669 4686 8 120 0 200 8 168 4 291 16 71 1 171 8 119 0 180 8 167 380 4526 16 70 23 460 8 118 8 281 8 166 380 4386 16 69 9 200 8 117 383 4196 8 165 2 241 8 68 1 200 8 116 0 210 8 164 10 511 8 67 80 1132 8 115 72 741 8 163 10 391 8 66 35 371 8 114 0 190 8 162 2 250 8 65 157 1271 8 113 21 410 8 161 0 210 8 64 2553 22332 8 112 278 3175 8 160 4 270 8 63 1 211 8 111 82 942 8 159 67 752 8 207 1 310 8	Backtrack for solutions and proof 1576 all 2633 2229 2552 1286 2801 1175 1160 66100 2879 1339 714 1139 8760 2606 3924 3277 3644 4336 3023 3284 1007 2500 1747 3410 2769 1788 1143 4884 13601 1183 1964 1273 5216 7380 813 1273 1491 7586 5307 823 1916 891 4337 3076 694 1384 1004 2741 3640 1600 1774 1140 5818 2833 1196 2132 4265 2199 3748 1371 1971 3874 23980 3796 1492 1704 1634 19810 3474 3316 1918 2895 6255 1864 1379 1261 2406 5484 5114 1604 2214 7019 4813 1298 1046 1578 1991 781 8396 6627 4114 1610 3701 3117 5080 1521 1840 1630 6475 855 1408 2608 8792 797 2408 8321 4047 5832 3291 4047 1694 1910 2375 3909 2381 2375 2065 3305 1669 875 2682 641 3158 1926 2203 1675 1549 3237 1799 3770 730 14250 10293 1526 7249 2976 1453 1389	Time for all 10074 15973 15943 17735 solutions 5648 21971 9033 8292 and proof 503094 21812 9292 4897 [ms] 59705 20029 35020 31886 5608 26959 23744 32917 6469 23604 10636 27900 18437 13249 14311 34299 120363 6689 13860 5958 28311 57443 5488 9754 6228 60416 34710 5728 10906 7802 31145 23624 4697 8893 4316 19408 21891 11867 15181 4827 48570 20620 8272 17826 5818 16884 23994 7351 10786 18747 240706 30674 11256 9824 21782 187911 24195 26718 13109 7941 39937 13880 7681 11817 16253 34359 31415 11847 14791 14311 26599 9404 14290 8602 17575 42230 4537 8182 29012 14902 83310 5758 9293 14781 52936 33709 24105 42671 22692 33298 49050 9995 11427 21961 44594 29702 17355 22392 3125 10115 28280 11988 11557 4917 18527 4236 7150 21201 17295 20770 4486 11256 10745 17324 5268 11596 32477 4556 32937 128915 81818 11126 50573 39878 18757 8673 10345

gen_rect([],[],[],[],[]

). gen_rect([X|Rx],[Y|Ry],[W|Rw],[H|Rh],[[X,Y,W,H]|R]) :-gen_rect(Rx,Ry,Rw,Rh,R). labeling([]). labeling([X|Y]) :lmindomain([X|Y], M), fix_min([X|Y], M, R), labeling(R).

Appendix A

The following table contains all the data for all problem instances from [START_REF] Bouwkamp | Duijvestijn Catalogue of Simple Perfect Squared Squares of orders 21 through 25[END_REF]. The fields describe the problem number, the number of squares, the size of the master square and a list of the square sizes. The problem number corresponds to the page number in [START_REF] Bouwkamp | Duijvestijn Catalogue of Simple Perfect Squared Squares of orders 21 through 25[END_REF]. Problems 166 and 167, 168 and 169, 182 and 183 are identical, but have two non-isomorphic solutions.


Raphaël Pesci

Karim Inal

Marcel Berveiller

Jean-Lu Lebrun

Internal Stress Analysis for the Damage Study of a 16MND5 Bainitic Steel

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Omar El-Rifai

Thierry Garaix

Vincent Augusto

Xiaolan Xie
email: xie@emse.fr

A Stochastic Optimization Model for Shift Scheduling in Emergency Departments

Excessive waiting time in Emergency Departments (ED) can be both a cause of frustration and more importantly, a health concern for patients. Waiting time arises when the demand for work goes beyond the facility's service capacity. ED service capacity mainly depends on human resources and on beds available for patients. In this paper, we focus on human resources organization in an ED and seek to best balance between service quality and working conditions. More specifically, we address the personnel scheduling problem in order to optimize the shift distribution among employees and minimize the total expected patients' waiting time. The problem is also characterized by a multi-stage re-entrant service process. With an appropriate approximation of patients' waiting times, we first propose a stochastic mixed-integer programming model that is solved by a sample average approximation (SAA) approach. The resulting personnel schedules are then evaluated using a discrete-event simulation model. Numerical experiments are then performed with data from a French hospital to compare different personnel scheduling strategies.

Introduction

Emergency departments are the crux of hospitals as they serve as a center for medical treatment and as a continuously operating portal for inpatient admissions [START_REF] Schuur | The growing role of emergency departments in hospital admissions[END_REF]. It is important when considering the ED, to provide patients with timely health care while maintaining a healthy work environment for the personnel. Further ED are facing the challenge of providing a good quality of service while being cost-effective.

Increasing demand for health services in recent years and a reduction of the number of resources in the health industry due to budget cuts have led to chronic overcrowding in the ED [START_REF] Jones | A multivariate time series approach to modeling and forecasting demand in the emergency department[END_REF]. From the patients' point of view, an overcrowded emergency department leads to excessive waiting times. While from a health specialist point of view, it can be the cause of an unfriendly work environment [START_REF] Rogers | The Working Hours Of Hospital Staff Nurses And Patient Safety[END_REF]. One of the major consequences of overcrowding with regards to the Operations Research community is a growing interest in developing strategies that make the most efficient use of the available resources [START_REF] Ernst | Staff scheduling and rostering: review of applications, methods and models[END_REF]. In France, a recent report of the first national ED conference which gathered health practitioners and government officials decreed the urgency of the problem of long waiting times [START_REF]Réduire les temps de passage aux Urgences[END_REF]. The report describes several strategies to reduce the potential bottlenecks that cause overcrowding which can also be found in the literature [START_REF] Sven | A systematic review of triagerelated interventions to improve patient flow in emergency departments[END_REF][START_REF] Akkerman | Reallocation of beds to reduce waiting time for cardiac surgery[END_REF][START_REF] Bruin | Dimensioning hospital wards using the Erlang loss model[END_REF]. Among the recommendations, the authors mention the need to improve patient flow, the need to organize inpatients bed resources and the need to make a better use of available resources.

As ED work around the clock, a turnover of health care professionals is necessary in order to ensure an appropriate service level at all time. In the context of emergency departments, the work demand across the day is subject to the rate of arrivals of patients. In turn, the arrival rate of patients is subject to stochastic variations. These variations make the problem of organizing the workforce difficult.

1 A look at the number of patients and the staff level in Figure 1 observed at our collaborating hospital clearly shows mismatch of patients flow and personnel schedules. At every time period, the activity is estimated with the average number of patients present in the ED. Figure 1 clearly shows a peak of activity in the evenings. The distribution of staffs on duty does not appropriately match the peaks of activities observed in the ED. These observations prompted us to design a personnel scheduling model to meet the work demand accurately while satisfying the work constraints imposed. In this paper, we present a human resources planning model for ED in order to explore the various human resources planning options and to take into account their robustness to face uncertainties and time-dependent patient arrival rate. The main efficiency criterion considered is the total patient waiting time. While other criteria are possible, total waiting time is a good indicator of patient service quality and ED crowding.

Personnel schedule is determined by a high-level stochastic programming model in which patient waiting times are approximated by a discrete time model and legal constraints for personnel scheduling are taken into account. The stochastic programming model is solved by a deterministic sample average approximation (SAA) [START_REF] Kleywegt | The sample average approximation method for stochastic discrete optimization[END_REF] and mixed-integer programming. The resulting personnel schedules are evaluated with a more realistic ED discrete-event simulation model.

Workforce Scheduling in Emergency Departments

With the growing work constraints and the expansion of work facilities, the employee scheduling problem has become, for most industries, manually intractable. As early as 1954, mathematical models have been used in order to computationally assist industries to schedule their employees to work [START_REF] Edie | Traffic delays at toll booths[END_REF]. Employee scheduling is often decomposed into several sub-problems [START_REF] Ernst | Staff scheduling and rostering: review of applications, methods and models[END_REF]. In a dynamic environment where the workload is variable, the first task is to be able to forecast the work demand for any scheduling horizon [START_REF] Mccarthy | The challenge of predicting demand for emergency department services[END_REF][START_REF] Jones | Forecasting daily patient volumes in the emergency department[END_REF][START_REF] Jones | A multivariate time series approach to modeling and forecasting demand in the emergency department[END_REF]. Once a work demand estimate is known, the decision maker has to determine an optimal staffing level for each scheduling period. After that, the staffing level is used as a basis for creating schedules that meet the staffing level appropriately. Finally, an assignment phase matches employees to the schedules created.

One of the peculiarities in scheduling studies conducted for health care systems is the focus on multiskill scheduling. The work demand in health care systems needs to be addressed by appropriately qualified personnel and cannot be distributed arbitrarily between the resources [START_REF] Warner | Scheduling nursing personnel according to nursing preference: A mathematical programming approach[END_REF][START_REF] Jaumard | A generalized linear programming model for nurse scheduling[END_REF]. This multi-skill characteristic adds yet another layer of complexity to the workforce scheduling in ED that is often overlooked. We identify two main types of human resources in this study (physicians and nurses) and schedule them simultaneously. Simultaneous scheduling of different resources means having the work of one type of resource impact the schedules of the other.

With regards to the staffing problem, the work demand in the ED cannot be determined easily and depends on the arrivals of patients in the services and on the interactions between different resources. Literature on staff scheduling in the ED often considers stochastic demand patterns in order to address the uncertainty in the arrival of patients [START_REF] Izady | Setting staffing requirements for time dependent queueing networks: The case of accident and emergency departments[END_REF][START_REF] Sinreich | Reducing emergency department waiting times by adjusting work shifts considering patient visits to multiple care providers[END_REF]. Sinreich and Jabali [START_REF] Sinreich | Staggered work shifts: a way to downsize and restructure an emergency department workforce yet maintain current operational performance[END_REF] propose a simulation model coupled with an linear programming model to downsize emergency departments while keeping the patients' length of stay (LOS) constant. The simulation model dictates the required staffing profile while the linear programming model chooses the schedules that matches this staffing profile best. Zeltyn et al. [START_REF] Zeltyn | Simulation-based models of emergency departments:: Operational, tactical, and strategic staffing[END_REF] propose a unique simulation model to address the human resources allocation problem on multiple decision levels (operation, tactical and strategic). On the tactical level, they use two strategies: the Rough Cut Capacity Planing (RCCP) technique and the Offered Load approach. The later being a refinement of the former that allows to split the workload. The staffing level is computer hourly through the analysis of a simple queuing system. Queuing theory has also been extensively used for determining staffing levels because it allows the evaluation of quality of service criteria such as the waiting time of patients [START_REF] Green | Using queueing theory to increase the effectiveness of emergency department provider staffing[END_REF][START_REF] Jennings | Server staffing to meet time-varying demand[END_REF]. However, analytic formulas for determining optimal staffing levels are often limited with regards to the complexity of the system they can analyze. As such, no analytical formula exists for optimally determining staffing level in systems with time-dependent stochastic arrival rates but instead numerical approximations have been used to alleviate this problem [START_REF] Vile | Time-dependent Stochastic Modelling for Predicting Demand and Scheduling of Emergency Medical Services[END_REF]. For example Ingolfsson et al. [START_REF] Ingolfsson | Accounting for time-varying queueing effects in workforce scheduling[END_REF] use an integrated queuing and optimization model to investigate the effect of time-varying arrival rates in the context of workforce scheduling. The authors present a heuristic to construct potential schedule then use queuing theory to evaluate the resulting service levels. Izady and Worthington [START_REF] Izady | Setting staffing requirements for time dependent queueing networks: The case of accident and emergency departments[END_REF] use queuing theory and simulation in an iterative staffing scheme to determine schedules for an ED modeled as an infinite server network with different tasks and processes. After an appropriate staffing level is determined, an ILP is employed to find shifts that minimize the over-staffing and under-staffing. Cochran and Roche [START_REF] Cochran | A multi-class queuing network analysis methodology for improving hospital emergency department performance[END_REF] suggest a queuing network approach for capacity planning of beds and staff in an ED. Their methodology is based on a multiple-customer class queuing analysis to increase ED access. Seasonal variation is of patient arrival patterns is considered.

In this work, the ED is represented as a multi-server multi-stage time-dependent priority queue system with re-entrance, Poisson distributed arrivals and exponential service times. Physicians and nurses are the servers whose cyclic schedules are to be determined. The planning horizon is discretized into time periods and a different mean arrival rate is given for every period of the day. We choose the expected total patient waiting time as the main criteria to optimize as it is often associated with patients' and staff well-being and is a good indicator of ED crowding. In order to address the complex nature of the problems, many studies addressed each sub-problem in the employee workforce scheduling independently by overlooking the inherent interaction between them [START_REF] Vile | Time-dependent Stochastic Modelling for Predicting Demand and Scheduling of Emergency Medical Services[END_REF]. As such, the optimal staffing level is often considered separately from the scheduling problem. This approach neglects complex waiting queues generated by uncertainties and changing patient arrival rate, which is crucial for ED crowding. In this work, we consider the staffing level as part of the scheduling problem and solve them both simultaneously. Integrating the two sub-problems allows us to schedule employees while reproducing the work dynamics of the system and simulating the patient flow.

The structure of the work is divided into two parts. First, we build an optimization model that finds staff schedules with operational constraints. Then, we evaluate more accurately our solution in a simulation model to measure the value of some of the assumptions in the optimization model. The optimization problem is solved by a Sample Average Approximation (SAA) solution methodology by considering different possible scenarios or realization of patient flows in the ED and optimizing the workforce scheduling problem considering all the different realizations. The remaining of the paper is organized as follows. We first describe the ED model in details in Section 3, then explain the optimization and the simulation models along with the assumptions that we consider in each in Sections 4 and 5 respectively. Then, we present numerical experiments in Section 6 and conclude in Section 7.

of an ED. Patients can follow different trajectories through the ED depending on their pathology and their severity level. A simplified scheme for their trajectory is as follows. Patients are first assessed by a nurse to determine whether they need vital medical care, a step often called triage. If they need vital care, they are taken in priority to a resuscitation room. Otherwise, the patients wait in a waiting room until a physician is available to examine them. After the first assessment, patients are either directly treated by a nurse or are sent for auxiliary exams to confirm the diagnostic. The auxiliary exams range from a simple blood test to time-consuming imaging tests. After the auxiliary exams are performed, the patients are examined again by a physician and finally treated by a nurse before being discharged or admitted into another department in the hospital [START_REF] Izady | Setting staffing requirements for time dependent queueing networks: The case of accident and emergency departments[END_REF]. The second assessment leads to re-entrance of patients to the same resources. Priority is given to patients that are being assessed for the second time by physicians over patients arriving for the first time.

Each process in the ED can be modeled as a queue where jobs to be served are the patients and the hospital staff are the servers that handle these jobs. Figure 2 illustrates a basic queuing system. On the left side, we have patient arrivals into the system which follow an inhomogeneous Poisson distribution of parameter λ t . Once in the system, they wait in a queue before being served by one of the servers of the system. The service time is assumed to be exponentially distributed in this paper. This corresponds to an M t /M/c queue model where M stands for the Markovian assumptions on the arrival and the service rates and c is the number of servers available.

Figure 2: A basic queue model that represents a process in an ED From the above, each ED has three main activities: assessment, treatment and auxiliary exams. The two health care providers that we model in the system are physicians and nurses. Physicians perform health assessments and examine the patients before and after having auxiliary exams. While nurses handle the final treatment administered to patients. As such, physicians are the servers for the two assessment queues of our system and nurses are the servers for the last queue. The overall system is illustrated in Figure 3. We do not schedule the lab technicians as they are often resources not entirely dedicated to the ED but instead delay patients that need to have auxiliary exams done.

Apart from the patients care pathway modeling, the queuing model also accounts for the stochastic nature of the work in the ED. The arrival rate of patients greatly varies during a day with peak arrival during the evenings and low arrival during the night. Service time needed to treat each patient can-not be known in advance. The stochastic service duration used in this paper models the different needs of patients. The number of patients arriving for each period of our planning horizon is parametrized by mean arrival rates that can be estimated from past data. Service distributions are estimated based on mean service times observed and found in the literature.

In Section 4, we present the mathematical formulation that we derived from the system specifications explained in this section. Some simplifying assumptions are made and detailed to solve the problem in reasonable time. Nonetheless, the solutions are evaluated using a more realistic simulation model whose specifications are explained in Section 5.

Mathematical Formulation

In this section, we present a mathematical formulation that is based on the details discussed in section 3. The objective is to optimize a waiting time performance indicator that we approximate with discrete time variables. As patients' waiting time is heavily influenced by the number of personnel present in the ED at any given time, the model finds optimal schedules over the planning horizon. For the time being, we are not interested in assigning the schedules to actual health personnel but only in finding an optimized schedule that meets the objective of minimizing the average total waiting time. In order to keep the model realistic, we allow for a limited number of different hours in which the staff can start working.

The number of patients arriving at each period t for each queue q is denoted A t,q (ω). This symbol is overloaded in our model as it both designates the Poisson distributed arrival to the first two queues and the variable designating the arrival of patients to the nurses' queue. After the patients enter a queue, they are counted as waiting in the variable W t,q (ω). Patients that are served in a queue, are counted in the variables O t,q (ω). Then, parameter L i,n,t (ω) determines the service capacity of the staff at each period. More specifically, the parameter determines the number of patients that n resources of type i can serve during period t. The ω index designates one scenario of patients' arrival and visit to the ED. The decision variables of the model x i,n,t define the number of health care providers working during period t.

To model the operational constraints, we use several parameters. A lower bound LBD and an upper bound U BD specify the possible lengths of the shifts that can be assigned to the personnel. The number of different hours of the day in which personnel of type i can start working is denoted by d i . An upper bound T T i is set on the number of total working hours that can be scheduled for resources of type i. The delay at the auxiliary exams is designated by δ. Finally some variables, p i,t , s i,t , e i,t and z i,t are used to facilitate the expression of shift constraints. Some modeling assumptions are made in order to be able to solve the problem efficiently. A discrete-time approximation of the waiting time is calculated as W t,q minus O t,q . What follows is that patients arriving during one period can only be treated in the next. Also, the waiting time of patients during the period they arrive is not counted. The patients are divided into two classes, those that do not require auxiliary exams and those that do. Each patients' class has a well defined trajectory from the start. This allows us to avoid having to route the patients to their appropriate queue after they have been assessed. This assumption is justifiable inasmuch as we do not care about each patients' trajectory individually. Instead, we keep track of the total number of patients in each queue of the system at each period t. We also assume the system is initially empty. We make this assumption by starting the scheduling horizon at an hour when it is reasonable to say that the system is clear. As of the data available to us, this is most true around 5:30 in the morning. Furthermore, we do not distinguish between the different health care providers with the same skills and consider all physicians to be homogeneous and all nurses to be homogeneous. For the sake of simplicity, we omitted the possible visit of patients to the vital care unit. As of the data available from the UHC of Lille for the year 2011 and 2012, patients that require vital prognosis only constitute 1.79% of the total patients and those in need of resuscitation merely 0.45%. Finally, auxiliary exams are always assumed to last for a fixed duration. Of course, these assumptions are relaxed in the simulation model and evaluated for consistency in the numerical experiments.

Parameter

(W t,q (ω) -O t,q (ω))] (1)
Subject to: O 1,q (ω) = 0 ∀q (4)

W 1,q (ω) = A 1,q (ω) ∀q ∈ {1, 2} (2)
W 1,q (ω) = 0 ∀q / ∈ {1, 2} (3)
W t,1 (ω) = W t-1,1 (ω) + A t,1 (ω) + O max(t-δ,1),2 -O t-1,1 (ω) ∀t ∈ [2..T ] ( 5
)
W t,2 (ω) = W t-1,2 (ω) + A t,2 (ω) -O t-1,2 (ω) ∀t ∈ [2..T ] (6)
W t,3 (ω) = W t-1,3 (ω) + O t-1,1 (ω) -O t-1,3 (ω) ∀t ∈ [2..T ] (7)
O t,q (ω) ≤ W t,q (ω) ∀t ∀q (8) q∈{1,2} O t,q (ω) ≤ N1 n=1 L 1,n,t (ω) * x 1,n,t (ω) ∀t (9)
O t,3 (ω) ≤ N2 n=1 L 2,n,t (ω) * x 2,n,t (ω) ∀t ( 10
) t p i,t ≤ T T i ∀i (11)
t (s i,t -e i,t ) = 0 ∀i (12)

p i,t = t t =2 (s i,t -e i,t ) + p i,1 ∀t ∈ [1..T ] ∀i (13)
p i,1 = p i,T + s i,1 -e i,1 ∀i (14)
mod(t+LBD,T +1)

t =t+1 e i,t ≤ p i,t -s i,t ∀t ∈ [1..T ] ∀i (15)
p i,t ≤ mod(t+U BD+1,T +1) t =t+1 e i,t ∀t ∈ [1..T ] ∀i ( 16
) Ni n=1 n * x i,n,t = p i,t ∀i ∀t (17)
Ni n=1

x i,n,t = 1 ∀i ∀t (18)

z i,t * N i >= s i,t ∀i ∀t (19)
z i,t ≤ d i ∀i ∀t (20)
p i,t ≥ 1 ∀i ∀t (21)
p i,t ≥ s i,t + 1 ∀i ∀t (22)
The objective function of the model ( 1) is a discrete time approximation of the total patients waiting time accumulated over all the scenarios. By summing over all the sample scenarios ω ∈ Ω considered, we seek to obtain a good approximation of the estimated value of the actual waiting time of patients in the system. Constraints (3) and (4) initialize the number of patients in the queues in each scenario. This assumption of empty patient queue at the beginning of a day corresponds to our field observation of several ED in which often no patient remains late in the night such as at 5.30am. This is true whether the ED is busier during the day or not. The next set of equations of the problem are balance equations to track the movement of patients through the ED. As such, (5), ( 6) and ( 7) update the number of patients in each queue at every time period.

Constraints (8) specify a logical bound for the number of patients served during a period. This number can-not be larger than the number of patients already in the queue. Then, constraints [START_REF] Jennings | Server staffing to meet time-varying demand[END_REF] limit the total number of staffing hours available for each type of resource. Constraints ( 12) through ( 14) are logical constraints that define the range of values the variables p, s and e can take. Then, constraints ( 15) and ( 16) define the type of shifts that are allowed by specifying a lower bound and an upper bound on the shift lengths.

Constraints [START_REF] Schuur | The growing role of emergency departments in hospital admissions[END_REF] and ( 18) assign a value to the variables x which is used to limit the number of patients served at every time period in constraints ( 9) and [START_REF] Jaumard | A generalized linear programming model for nurse scheduling[END_REF]. Finally constraints [START_REF] Sinreich | Reducing emergency department waiting times by adjusting work shifts considering patient visits to multiple care providers[END_REF] through [START_REF] Warner | Scheduling nursing personnel according to nursing preference: A mathematical programming approach[END_REF] are operational constraints to make sure that there are at most d i different moments of the day where resources of type i start working and that there is at least one resource of each type in every time period. We need to have an upper bound on the number of periods resources can start working in because physicians and nurses need to transmit patients' information to their colleagues before leaving the service. As such, constraints [START_REF] Warner | Scheduling nursing personnel according to nursing preference: A mathematical programming approach[END_REF] make sure that the shifts of the different resources overlap.

Simulation Model

In order to evaluate the validity of the assumptions and simplifications made in the mathematical formulation, we test the solutions obtained using a simulation model. The simulation model includes several features not considered in the optimization model for a more accurate representation of the ED. The simulation improves over the mathematical formulation on several fronts: system specification, patients specification and resources activities.

As in the optimization model, exponential service time and Poisson arrival rates are used in the simulation model. However, the simulation uses a discrete event paradigm instead of the discrete time model paradigm in the optimization. This is important as one of the key assumptions in the optimization problem is that a discrete time model can serve to calculate patients' waiting time accurately. Second, the vital care pathway is no longer ignored. Patients requiring vital prognosis are attended by a physician and a nurse immediately. Auxiliary exams are divided into two categories: biological exams and diagnostic imaging exams. The biological exams make for 55% of the examinations.

Patients on the other hand are no longer homogeneous but are classified according to their acuity level (from 1 to 5 in increasing order of severity). The classification scheme used is that of the french hospitals known as CCMU (Classification Clinique des Malades des Urgences). This allows for a sickest-first prioritization in the different queues of the ED. Additionally, the patients that are classified as CCMU 4 or CCMU 5 skip the regular route and are treated in priority in the vital care pathway. These patients are served by both a nurse and a physician for a random duration.

Finally, with regards to resources, physicians are no longer solely occupied by patients. Instead, they also have to respond to administrative tasks such as searching patients' historical data or finding downstream beds for inpatients. Beds are not explicitly modeled in the system. instead, patients waiting in queues are considered to be assigned either a room or a hallway bed. Administrative tasks occur at random intervals but follow the same distribution laws as patient arrivals. Having the same distribution makes sense if we think that administrative tasks are more likely to occur when patients enter the system. Nurses no longer solely handle patients' treatments but also take care of the biological exams.

This improved ED model is used to evaluate the solutions obtained with the optimization model. Important system statistics such as the estimated patients' waiting time and the time spent in the system can be extracted and a confidence interval for these data can be calculated.

Numerical Experiments

Numerical experiments have been conducted to show the potential reduction of patients' waiting time in an ED by adjusting the workforces' shift schedules. The parameters settings for the optimization model, the simulation model, and the real-world benchmark are described in Section 6.1. In Section 6.2, we check that the optimization model matches enough key features of the queuing system to provide suitable solutions. In order to do so, the schedules obtained by running the optimization model are evaluated and compared to simulation tests. The performances of different scheduling strategies are also compared. The strategies differ in the constraints on the employees' shifts. In Section 6.3 the robustness of schedules facing significant demand increase is evaluated and in Section 6.4 the transition from week schedules to week-end schedules is investigated. The optimization model is solved using a commercial solver (IBM Ilog Cplex 12.5). All tests were performed on an Intel(R) Xeon(R) 4 cores CPU E5520 @ 2.27GHz with 8 MB of cache memory and 8 GB of RAM. The evaluation is done using AnyLogic 6.9.

Parameters Settings and Real-World Data

The parameters used in the experimental section are based on data for the years 2011-2012 of patients' entries in the pediatrics' ED of the University Hospital Center (UHC) of Lille France. We extracted from these data an average number of patients' arrivals per half-hour of the day and used them as parameters for the Poisson distributions of each period. The length of the period also called time step ∆ = 30 minutes is chosen primarily to coincide with the shifts change marks used in the hospital. Also, the patients' arrival averages are calculated according to length of the time step. Average values would be less significant if they relied on less data. As the data available only contain arrival and departure times, we did not determine the service duration of patients in the same manner as the arrivals. Instead, we relied on average service durations found in the literature and on experts' opinion to calibrate these parameters. However, the patients' lengths of stays obtained with these service durations are of the same order of magnitude than those observed in the real data. Based on these parameters, we generate several possible scenarios for patients arrivals and service time.

In the UHC of Lille, the number of resources working during the week differs drastically from the number of resources working on the weekends. Therefore, two set of optimization tests are run. The first set of tests uses the total hours available during the week and the arrival distributions extracted for weekdays and the second set of tests relies on data for weekends. The total staff hours available serves as an upper bound on the number of hours scheduled. Each scenario in the optimization model is generated by obtaining a realization ω of the distributions of the arrival and service times. Once the different parameters of the models are set, the scenarios are the only input to the model.

After the schedules are optimized, we evaluate these schedules and the waiting time they produce in the simulation model. As the schedules we obtain are cyclical, the evaluation is done for several consecutive days. We introduce a warm-up period in order to reach stable results. Evaluating the performance of the schedules can be done on a large number of replications which allows the calculation a confidence interval. The parameters used for the experimental runs are summarized in Table 3.

To validate the pertinence of a scenario-based strategy, we ran three different instances of optimization tests. After obtaining the resulting schedules from each test, we evaluated the evolution of the expected waiting time in the simulation model using the resulting schedules. We see the results in Figure 4, where the average waiting time over all three tests are displayed. The average waiting time of the patients decreases with the increase of the number of scenarios used in the optimization. We notice that the expected waiting time stabilizes for about |Ω| = 100 scenarios.

In terms of computational feasibility, our model strongly depends on the number of scenarios used for optimization and on the time step chosen. This is because the number of variables increases linearly with the number of scenarios and the number of periods. Fortunately, we showed that the simulated average waiting time converges with solutions obtained from 100 scenarios. From our tests, all solutions with 100 scenarios were obtained within 4 hours with a 0.11% relative gap in the worst case.

Other variables in the model also influence the running time but to a lesser extent. As such, the number of patients in the queues does not increase the number of variables but only the range the variables can take. On the other hand, the number of variables representing resources: x i,n,t increases as the number of resources increases but is independent of the number of scenarios. For example, in a large ED with about 200 patients daily and 10 to 15 resources of each type, the optimal solution is still obtained within the 4 hours of computation time. Throughout the tests, the schedules obtained are compared against the reference schedule ref that is currently put in place in the ED of Lille. The first alternative we test are schedules obtained with shifts whose length is fixed to 8 hours for both physicians and nurses. The schedules obtained with this strategy give us an idea of how much we can improve by changing the starting times along with more restrictive shift constraints than what is currently practiced. Then, we keep the shifts lengths range from 4 to 12 hours as it is the case in the UHC of Lille but modify the starting times. The third and final alternative is to create shifts with no constraint on their lengths. In this last case, the only constraint is the total number of hours worked. This last strategy is not realistic but merely allows us to compute a lower bound on the total waiting time of the patients in the optimization model. Constraints on every schedule strategies are given in Table 4. The parameters d i and T T i refer to the parameters presented in the optimization model.

Parameter

Shifts

Analysis of optimized schedules

We verify in a first part that the optimized schedules obtained can actually reduce the average waiting time of patients compared to the reference schedule when tested in simulation conditions. In order to do so, we allow for different assumptions on the shifts as explained previously. In all the alternatives tested, we allow the shifts to start on any period of the day but only permit a limited number of shift changes to occur as it is currently the case. Our analysis in this subsection focuses on weekends but results obtained for weekdays are similar and lead to the same conclusions. Table 5 summarizes the results obtained for schedules optimized to weekend conditions. The first column specifies the shifts type used. The number next to the shift type specifies different runs of the optimization problem, each with a different set of scenarios. W -O, is total objective value converted to average hour units per patient. This value gives an approximation on the expected waiting time as computed in the optimization model. The total waiting time associated with each queue, q in the optimization model is given in the columns W q -O q , for q ∈ {1, 2, 3}. Next, the Gap indicates the percentage distance between best lower and upper bounds obtained. An time limit of 240 minutes is used for each run of the optimization problem.

Shifts W-O (h) W 1 -O 1 (h) W 2 -O 2 (h) W 3 -O 3 (h) Gap % Time(mins) ref (
The simulation model of Section 5, is used to evaluate each solution schedule. The results of the evaluation are presented in Table 6. The simulation is run over 100 replications each of 10 successive days. Since we only have a unique reference schedule, only one row keeps track of the expected waiting time using the reference schedule. Table 6 first notes the patients' expected waiting time in W ait. The next columns indicated the expected and max waiting time in the physicians' queue (Q phys ) and nurses' final treatment queue (Q nurs ). The results consistently reduce the expected waiting time of patients in more flexible strategies. This is further confirmed by the similar behavior of the maximum expected waiting time over all replications, in columns Qmax. However, within this trend of reducing the average waiting time, the results are not very stable. Different set of scenarios in the optimization model result in different waiting time in the simulation model. For example with shifts of fixed length the waiting time varies between 79.94 minutes for the first set of scenarios and 87.28 for the second set. This can be explained by the approximations that were made and the fact that we only used a limited number of scenarios (σ= 100) to optimize.

Having schedules with very flexible shifts (in the range 1-24 hours) does not really improve the waiting time of patients compared to schedules with shifts in the range of 4-12 hours. Consequently, we note that modifying the starting times of the shifts is the factor that helps us most improve the expected waiting times, even with a constant total working hours capacity. As explained in Section 3, queues 1 and 2 are served by physicians and queue 3 is served by nurses. In the simulation model, the first two queues are grouped into a physician queue Q phys .

Figures 5, 6 and 7 present example staffing levels obtained by optimizing the schedules following the different shift strategies. The reference schedule clearly misses the evening peak of arrivals while the other strategies closely match the arrival curve. We note that the solution with shifts varying from 4 to 12 hours and those varying from 1 to 24 hours are very close and match the ascending and descending slopes of the arrivals flow. These observations corroborate the results obtained in Tables 5 and6.

Validation of the Robustness of the Schedules

In this set of simulation runs, we increment the number of patients arriving to the system and evaluate the schedules obtained through optimization again. In order to do so, we increment the average Poisson parameter for each period. This increase in the number of arrivals results in about 15% more patients per period. Since the schedules have not been optimized to perform under these conditions, the evaluation allows us to have an idea on the robustness of the schedules during epidemic periods.

The results of this set of tests are summarized in Tables 7 and8. We notice from the results that the expected waiting time increases proportionally to the results of Section 6.2 across the different solution instances. Consequently, the optimized schedules still outperform the reference schedule in unfavorable conditions. This coherence in the waiting times is still true for each queue individually but we see a disproportionate increase in the average maximum waiting times for schedules constructed with (1-24h) constraints. This increase is understandable as the (1-24h) constraints specifically exploit the discrete-time approximation of the optimization model. Schedules with (1-24h) constraints can be expected to perform poorly under different conditions than those they were optimized for.

Shifts

Week to Weekend Transition

One limitation of the previous results is that the cyclical schedules are evaluated independently for weekdays and weekends. In order to have simulation conditions that are closer to reality, the transition from a weekday schedule to a weekend schedule needs to be evaluated. Consequently, weekday and weekend schedules are concatenated and tested under both the average conditions for arrivals and the epidemic conditions. The results in Tables 9 and10 show that the transition from the week schedules to the weekend schedules can be done without real complications. The average waiting times observed when the concatenation of schedules is performed follow the results obtained previously. Indeed, the optimize schedules outperform the reference schedules under normal and epidemic conditions.

The results of the final set of tests are summarized in Figures 8 and9. As we have mentioned, the schedules with shifts varying from 1 to 24 hours do not give a significant advantage over the schedules with shifts varying from 4 to 12 hours. Finally, Figures 8 and 9 illustrate the differences in sojourn time according to the CCMU of patients. Coherently, patients with higher priority spend less time in the system than those with lower priority.

Conclusion

We presented in this study a stochastic optimization formulation for an ED shift scheduling problem. Schedules of interest are computed once for some representative period such as winter or summer season or weekdays and weekends. During this period, the schedules are repeated daily. The formulation follows the flow of patients through the ED and takes into account the stochastic nature of arrivals and service times. The constraints used to follow the evolution of patients through the ED mirror those in stock management literature. The optimization model is solved using a deterministic Sample Average Approximation (SAA) method. This model assumes many simplifications in order to be tractable in a commercial linear programming solver. The solution schedules obtained by the optimization model are evaluated using a discrete event simulation model under various test conditions to verify the validity of the assumptions and the robustness of the solutions. The simulation model used to evaluate the results has several improvements over the optimization model. Improvements such as patient prioritization and multiple tasks for resources (administrative tasks for physicians and biological exams for nurses) make the simulation model more realistic. The evaluation demonstrates that the assumptions made hold for optimizing the average waiting time. Nonetheless, the simulation model remains a simplification of the ED and hence lacks some elements that can be taken into account in future research. For example the service time of patients is approximated with an exponential distribution. Even if this approximation is very common in the literature, service times are usually not exponentially distributed. Hence, a possible extension to our work would be to test detailed patient service times to be collected from our partner hospital. Furthermore, excessive waiting time is only one aspect among many that define overcrowding in ED. Average waiting time is important and influences other criteria such as the maximum waiting time but other service quality criteria such as limited bed capacity contribute to overcrowding and should be considered in future research.

Real data from the UHC of Lille France are used to evaluate the quality and the performance of the schedules obtained. The proposed optimization approach allows to find promising cyclic schedules. The framework proposed is easily extensible and allows the integration of several strategical and operational constraints that are not possible using other methods.

Figure 1 :

 1 Figure 1: Preliminary analysis of the data from the pediatrics ED of Lille Hospital

Figure 3 :

 3 Figure 3: Process chart of patients in the ED

3 :

 3 Parameters of the system and of the optimization algorithm

Figure 4 :

 4 Figure 4: The evolution of waiting time with the number of scenarios used in the optimization model

Figure 5 :Figure 6 :

 56 Figure 5: Resulting staffing levels with shifts fixed to 8 hours

Figure 7 :

 7 Figure 7: Resulting staffing levels for shifts varying from 1-12 hours

Figure 8 :

 8 Figure 8: Expected sojourn time for the different acuity classes in normal conditions

Figure 9 :

 9 Figure 9: Expected sojourn time for the different acuity classes in epidemic conditions

Table 1 :

 1 The parameters of the scheduling model.

		Meaning		Domain
	t	The different periods	[1..T ]
		of the day
	n	Number of resources of	[1..N i ]
		each type
	i	The different type of	{1,2}
		resources available
	q	The queue number	[1..Q]
	d i	Maximum number of	[1..S i ]
		different starts for each
		resource
	T T i	Number of working
		hours available to
		schedule for resource
		of type i
	LBD	Lower bound on the	[1..47]
		shift length allowable
	U BD	Upper bound on the	[2..48]
		shift length allowable
	δ	The delay at the auxil-
		iary exams
	A t,q (ω)	Poisson distributed ar-	N+
		rivals in period t for
		queue q
	L i,n,t (ω)	Sum of exponentially	N+
		distributed	consul-
		tation times for n
		resources of type i at
		period t
		Objective:
	min	Ẽ

ω∈Ω [ t q

Table 2 :

 2 The variables of the scheduling model.

	Variable	Meaning	Domain
	x i,n,t	= 1 if there are exactly	{0,1}
		n resources working at
		period t
	p i,t	Number of resources of	[1..N i ]
		type i present at pe-
		riod t
	s i,t	Number of resources of	[1..N i ]
		type i starting at pe-
		riod t
	e i,t	Number of resources of	[1..N i ]
		type i ending at period
		t
	z i,t	Indicator variable that	{0,1}
		is equal to 1 if at least
		one resource of type i
		starts working at t
	W t,q (ω)	Number of patients	R+
		that are waiting in
		queue q a the start of
		period t
	O t,q (ω)	Number of patients in	R+
		queue q that have been
		served during period t

Table 4 :

 4 Parameters for the scheduling strategies

		LBD(h) U BD(h)	d 1	d 2	T T 1	T T 2
	ref	4	12	7 (week), 5(weekend) 6 178 (week), 130 (weekend) 124
	(4-12h)	4	12	7 (week), 5(weekend) 6 178 (week), 130 (weekend) 124
	(8h)	8	8	7 (week), 5(weekend) 6 178 (week), 130 (weekend) 124
	(1 -24h)	1	24	7 (week), 5(weekend) 6 178 (week), 130 (weekend) 124

Table 5 :

 5 Optimization results for the weekend scenarios

	1)	0.41	1730	2581	1131	0.00	0
	(2)	0.39	1767	2175	1200	0.00	0
	(3)	0.47	1919	2939	1289	0.00	0
	8hrs (1)	0.27	1449	1415	683	0.02	240
	(2)	0.26	1416	1311	738	0.01	240
	(3)	0.31	1645	1690	735	0.11	240
	4-12hrs (1)	0.21	1283	1101	431	0.00	20
	(2)	0.20	1143	1133	415	0.00	34
	(3)	0.24	1334	1356	522	0.00	42
	1-24hrs (1)	0.21	1236	1109	469	0.00	180
	(2)	0.20	1112	1130	400	0.00	53
	(3)	0.24	1340	1342	523	0.00	51

Table 6 :

 6 ShiftsWait (mins) Q phys (mins) Q phys max (mins) Q nurs (mins) Q nurs max (mins) Simulation results for the weekend schedules

	ref	125,23 ± 6,22	99,24 ± 5,6	977,07 ± 43,24	25,99 ± 1,5	422,26 ± 23,08
	8hrs (1)	79,94 ± 4,85	54,66 ± 4,21	861,81 ± 46,89	25,28 ± 1,75	465,12 ± 38,88
	(2)	87,28 ± 5,8	59,23 ± 5,13	910,35 ± 49,42	28,05 ± 1,87	499,15 ± 35,75
	(3)	83,82 ± 7,61	58,73 ± 6,75	864,21 ± 52,76	25,09 ± 1,94	412,3 ± 39,43
	4-12hrs (1)	62,43 ± 4,3	48,91 ± 4,08	832,22 ± 43,19	13,52 ± 1,01	315,84 ± 32,67
	(2)	63,13 ± 4,27	48,54 ± 3,94	836,35 ± 47,46	14,59 ± 0,99	374,4 ± 34,51
	(3)	64,46 ± 3,7	50,72 ± 3,38	842,91 ± 34,87	13,74 ± 0,85	330,66 ± 28,85
	1-24hrs (1) 64,97 ± 4,58	51,06 ± 4,35	871,16 ± 45,72	13,91 ± 0,9	324 ± 30,8
	(2)	64,55 ± 4,63	49,92 ± 4,16	855,5 ± 51,54	14,63 ± 1,02	369,36 ± 31,86
	(3)	64,27 ± 4,16	48,18 ± 3,69	831,45 ± 42,44	16,09 ± 1,09	381,73 ± 28,08

Table 7 :

 7 Simulation results for weekend schedules during an epidemic period

		Wait	Q phys	Q phys max	Q nurs	Q nurs max
	ref	242,99 ± 13	199,2 ± 12,06 1323,22 ± 34,88 43,79 ± 2,61 559,87 ± 23,74
	8hrs (1)	191,1 ± 13,32 148,81 ± 13,11 1268,17 ± 45,18 42,29 ± 3,38 651,89 ± 45,34
	(2)	191,32 ± 14,34 147,3 ± 13,53	1263,06 ± 44,4 44,02 ± 2,76 659,27 ± 45,17
	(3)	190,79 ± 13,6	149,63 ± 13	1264,16 ± 45,57 41,16 ± 2,75 633,62 ± 44,73
	4-12hrs(1) 148,74 ± 10,86 125,47 ± 10,41 1218,3 ± 46,26 23,27 ± 1,85 457,79 ± 43,62
	(2)	155,95 ± 12,04 133,56 ± 11,65 1258,86 ± 43,98 22,38 ± 1,43 494,64 ± 38,05
	(3)	150,31 ± 10,85 128,77 ± 10,55 1242,86 ± 42,17 21,54 ± 1,59 477,88 ± 42,1
	1-24hrs(1) 155,76 ± 13,59 132,37 ± 13,22 1211,44 ± 49,83 23,39 ± 1,69 482,04 ± 41,28
	(2)	146,73 ± 11,67 122,78 ± 11,02 1212,04 ± 49,64 23,94 ± 1,56 501,33 ± 30,19
	(3)	142,07 ± 11,26 116,78 ± 10,99 1180,36 ± 48,45 25,29 ± 1,65 525,05 ± 34,35

Table 8 :

 8 Simulation results for week schedules during an epidemic period ± 45,73 14,53 ± 1,07 367,98 ± 32,31 (2) 37,53 ± 2,84 22,7 ± 2,36 635,05 ± 46,81 14,83 ± 1,21 379,98 ± 34,39 (3) 38,48 ± 2,52 24,53 ± 2,27 654,31 ± 38,59 13,96 ± 0,97 384,82 ± 35,06

	Shifts	Wait	Q phys	Q phys max	Q nurs	Q nurs max
	ref	142,07 ± 11,26 116,78 ± 10,99 1180,36 ± 48,45 25,29 ± 1,65 525,05 ± 34,35
	8hrs (1)	53,14 ± 3,17	22,53 ± 1,57	564,15 ± 42,14 30,61 ± 2,24 490,83 ± 37,81
	(2)	51,62 ± 2,43	24,45 ± 1,45	580,12 ± 35,17 27,17 ± 1,65 502,22 ± 40,29
	(3)	54,54 ± 3,44	21,44 ± 1,64	528,09 ± 42,14	33,1 ± 2,56	601,4 ± 36,75
	4-12hrs(1)	37,5 ± 2,62	12,8 ± 0,92	498,01 ± 30,82	24,7 ± 2,27 541,19 ± 37,59
	(2)	35,11 ± 2,19	11,41 ± 0,75	417,48 ± 30,38	23,7 ± 1,82 524,46 ± 35,12
	(3)	40,63 ± 2,4	13,58 ± 0,86	510,86 ± 28,89	27,05 ± 2	598,06 ± 36,17
	1-24hrs(1)	35,22 ± 1,95	11,94 ± 0,85	454,38 ± 30,59 23,27 ± 1,61 526,38 ± 33,41
	(2)	35,39 ± 2,13	11,22 ± 0,73	414,83 ± 31,08 24,18 ± 1,77 549,97 ± 37,27
	(3)	40,91 ± 2,22	13,32 ± 0,8	499,48 ± 26,06 27,59 ± 1,88 602,06 ± 32,71
	Shifts	Wait	Q phys	Q phys max	Q nurs	Q nurs max
	ref	88,59 ± 4,7 62,59 ± 3,63 871,71 ± 45,88 25,99 ± 1,92 430,6 ± 30,93
	4-12hrs (1) 38,52 ± 2,96	24 ± 2,4	652,18

Table 9 :

 9 Simulation results for schedules spanning an entire week

	Shifts	Wait	Q phys	Q phys max	Q

nurs Q nurs max ref 140,58 ± 6,63 102,25 ± 5,09 1087,38 ± 42,46 38,33 ± 2,52 561,25 ± 43,66 4-12hrs(1) 69,31 ± 4,59 46,13 ± 3,74 941,59 ± 56,66 23,18 ± 1,9 504,82 ± 36,2 (2) 70,46 ± 4,61 43,86 ± 3,72 941,56 ± 49,76 26,6 ± 2,05 573,99 ± 38,44 (3) 71,09 ± 4,06 46,73 ± 3,32 974,81 ± 51,18 24,36 ± 2,12 508,22 ± 43,41

Table 10 :

 10 Simulation results for schedules spanning an entire week during an epidemic period

Emergency Department Queuing NetworkIn order to dynamically match the staff schedules and patient flow, we include in our system the dynamics of the patient flow through the ED. We present a queuing system with the key treatment steps and resources

acknowledgements

This research is partially supported by the French National Research Agency "ANR" under grant number ANR-11-TECS-0010 of project HOST.