Antoine Legrain
email: antoine.legrain@polymtl.ca

Jérémy Omer
email: jeremy.omer@insa-rennes.fr

Samuel Rosat
email: samuel.rosat@polymtl.ca

An Online Stochastic Algorithm for a Dynamic Nurse Scheduling Problem

Keywords: Stochastic Programming, Nurse Rostering, Dynamic problem, Sample Average Approximation, Primal-dual Algorithm, Scheduling

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In western countries, hospitals are facing a major shortage of nurses that is mainly due to the overall aging of the population. In the United Kingdom, nurses went on strike for the first time in history in May 2017. Nagesh [START_REF] Nagesh | Nurses could go on strike for the first time in british history[END_REF] says that "It's a message to all parties that the crisis in nursing recruitment must be put center stage in this election". In the United States, "Inadequate staffing is a nationwide problem, and with the exception of California, not a single state sets a minimum standard for hospital-wide nurse-to-patient ratios." [START_REF] Robbins | We need more nurses[END_REF]. In this context, the attrition rate of nurses is extremely high, and hospitals are now desperate retain them. Furthermore, nurses tend to often change positions, because of the tough work conditions and because newly hired nurses are often awarded undesired schedules (mostly due to seniority-based priority in collective agreements). Consequently, providing high quality schedules for all the nurses is a major challenge for the hospitals that are also bound to provide expected levels of service.

The nurse scheduling problem (NSP) has been widely studied for more than two decades (refer to [START_REF] Burke | The state of the art of nurse rostering[END_REF] for a literature review). The NSP aims at building a schedule for a set of nurses over a certain period of time (typically two weeks or one month) while ensuring a certain level of service and respecting collective agreements. However, in practice, nurses often know their wishes of days-off no more than one week ahead of time. Managers therefore often update already-computed monthly schedules to maximize the number of granted wishes. If they were able to compute the schedules on a weekly basis while ensuring the respect of monthly constraints (e.g., individual monthly workload), the wishes could be taken into account when building the schedules. It would increase the number of wishes awarded, improve the quality of the schedules proposed to the nurses, and thus augment the retention rate.

The version of the NSP that we tackle here is that of the second International Nurse Rostering Competition of 2015 (INRC-II) [START_REF] Ceschia | The second international nurse rostering competition[END_REF], where it is stated in a dynamic fashion. The problem features a wide variety of constraints that are close to the ones faced by nursing services in most hospitals. In this paper, we present the work that we submitted to the competition and which was awarded second prize.

Literature review

Dynamic problems are solved iteratively without comprehensive knowledge of the future. At each stage, new information is revealed and one needs to compute a solution based on the solutions of the previous stages that are irrevocably fixed. The optimal solution of the problem is the same as that of its static (i.e., offline) counterpart, where all the information is known beforehand, and the challenge is to approach this solution although information is revealed dynamically (i.e., online).

Four main techniques have been developed to do this: computing an offline policy (Markov decision processes [START_REF] Puterman | Markov Decision Processes: Discrete Stochastic Dynamic Programming[END_REF] are mainly used), following a simple online policy (Online optimization [START_REF] Buchbinder | Designing Competitive Online Algorithms via a Primal-Dual Approach[END_REF] studies these algorithms), optimizing the current and future decisions (Stochastic optimization [START_REF] Birge | Introduction to Stochastic Programming[END_REF] handles the remaining uncertainty), or reoptimizing the system at each stage (Online stochastic optimization [START_REF] Van Hentenryck | Online Stochastic Combinatorial Optimization[END_REF] provides a general framework for designing these algorithms).

Markov decision processes decompose the problem into two different sets (states and actions) and two functions (transition and reward). A static policy is pre-computed for each state and used dynamically at each stage depending on the current state. Such techniques are overwhelmed by the combinatorial explosion of problems such as the NSP, and approximate dynamic programming [START_REF] Powell | Approximate Dynamic Programming: Solving the Curses of Dimensionality[END_REF] provides ways to deal with the exponential growth of the size of the state space. This technique has been successfully applied to financial optimization [START_REF] Bäuerle | Markov Decision Processes with Applications to Finance[END_REF], booking [START_REF] Patrick | Dynamic multipriority patient scheduling for a diagnostic resource[END_REF], and routing [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF] problems. In Markof decision processes, most computations are performed before the stage solution process, therefore this technique relies essentially on the probability model that infers the future events.

Online algorithms aim at solving problems where decisions are made in real-time, such as online advertisement, revenue management or online routing. As nearly no computation time is available, researchers have studied these algorithms to ensure a worst case or expected bound on the final solution compared to the static optimal one. For instance, Buchbinder [START_REF] Buchbinder | Designing Competitive Online Algorithms via a Primal-Dual Approach[END_REF] designs a primal-dual algorithm for a wide range of problems such as set covering, routing, and resource allocation problems, and provides a competitive-ratio (i.e., a bound on the worst-case scenario) for each of these applications. Although these techniques can solve very large instances, they cannot solve rich scheduling problems as they do not provide the tools for handling complex constraints.

Stochastic optimization [START_REF] Birge | Introduction to Stochastic Programming[END_REF] tackles various optimization problems from the scheduling of operating rooms [START_REF] Denton | Optimal allocation of surgery blocks to operating rooms under uncertainty[END_REF] to the optimization of electricity production [START_REF] Fleten | Short-term hydropower production planning by stochastic programming[END_REF]. This field studies the minimization of a statistical function (e.g., the expected value), assuming that the probability distribution of the uncertain data is given. This framework typically handles multi-stage problems with recourse, where first-level decisions must be taken right away and recourse actions can be executed when uncertain data is revealed. The value of the recourse function is often approximated with cuts that are dynamically computed from the dual solutions of some subproblems obtained with Benders' decomposition. However these Benders-based decomposition methods converge slowly for combinatorial problems. Namely, the dual solutions do not always provide the needed information and the solution process therefore may require more computational time than is available. To overcome this difficulty, one can use the sample average approximation (SAA) [START_REF] Kleywegt | The sample average approximation method for stochastic discrete optimization[END_REF] to approximate the uncertainty (using a small set of sample scenarios) during the solution and also to evaluate the solution (using a larger number of scenarios).

Finally, online stochastic optimization [START_REF] Van Hentenryck | Online Stochastic Combinatorial Optimization[END_REF] is a framework oriented towards the solution of industrial problems. The idea is to decompose the solution process in three steps: sampling scenarios of the future, solving each one of them, and finally computing the decisions of the current stage based on the solution of each scenario. Such techniques have been successfully applied to solve large scale problems as on-demand transportation system design [START_REF] Bent | Scenario-based planning for partially dynamic vehicle routing with stochastic customers[END_REF] or online scheduling of radiotherapy centers [START_REF] Legrain | Online stochastic optimization of radiotherapy patient scheduling[END_REF]. Their main strength is that any algorithm can be used to solve the scenarios.

Contributions

The INRC-II challenges the candidates to compute a weekly schedule in a very limited computational time (less than 5 minutes), with a wide variety of rich constraints, and with important correlations between the stages. Due to the important complexity of this dynamic NSP, none of the tools presented in the literature review allows to solve this problem. We therefore introduce an online stochastic algorithm that draws inspiration from the primal-dual algorithms and the SAA. In that method,

• the online stochastic algorithm offers a framework to solve rich combinatorial problems;

• the primal-dual algorithm speeds up the solution by inferring quickly the impact of some decisions;

• the SAA efficiently handles the important correlations between weeks without increasing tremendously the computational time.

Finally, the algorithm uses a free and open-source software as a subroutine to solve static versions of the NSP. It is described in details in [START_REF] Legrain | A rotation-based branch-and-price approach for the nurse scheduling problem[END_REF] and summarized in Section 3.

We emphasize that the algorithm described in this article has been developed in a time-constrained environment, thus forcing the authors to balance their efforts between the different modules of the software.

The resulting code is shared in a public Git repository [START_REF] Legrain | Dynamic nurse scheduler[END_REF] for reproduction of the results, future comparisons, improvements and extensions. The remainder of the article is organized as follows. In Section 2, we give a detailed description of the NSP as well as the dynamic features of the competition. In Section 3, we state a static formulation and summarize the algorithm that we use to solve it. In Section 4, we present the dynamic formulation of the NSP, the design of the algorithm, and the articulation of its components.

In Section 5, we give some details on the implementation of the algorithm, study the performance of our method on the instances of the competition, and compare them to those obtained by the other finalist teams.

Our concluding remarks appear in Section 6.

The Nurse Scheduling Problem

The formulation of the NSP that we consider is the one proposed by Ceschia et al. [START_REF] Ceschia | The second international nurse rostering competition[END_REF] in the INRC-II, and the description that we recall here is similar to theirs. First, we describe the constraints and the objective of the scheduling problem Then, we discuss the challenges brought in by the uncertainty over future stages.

The NSP aims at computing the schedule of a group of nurses over a given horizon while respecting a set of soft and hard constraints. The soft constraints may be violated at the expense of a penalty in the objective, whereas hard constraints cannot be violated in a feasible solution. The dynamic version of the problem considers that the planning horizon is divided into one-week-long stages and that the demand for nurses at each stage is known only after the solution of the previous stage is computed. The solution of each stage must therefore be computed without knowledge of the future demand.

The schedule of a nurse is decomposed into work and rest periods and the complete schedules of all the nurses must satisfy the set of constraints presented in Table 1. Each nurse can perform different skills (e.g., Head Nurse, Nurse) and each day is divided into shifts (e.g., Day, Night). Furthermore, each nurse has signed a contract with their employers that determines their work status (e.g., Full-time, Part-time) and work agreements regulate the number of days and weekends worked within a month as well as the minimum and maximum duration of work and rest periods. For the sake of nurses' health and personal life and to ensure a sufficient level of awareness, some successions of shifts are forbidden. For instance, a night shift cannot be followed by a day shift without being separated by at least one resting day. The employers also need to ensure a certain quality of service by scheduling a minimum number of nurses with the right skills for each shift and day. Finally, the length of the schedules (i.e., the planning horizon) can be four or eight weeks.

Hard constraints

H1 Single assignment per day: A nurse can be assigned at most one shift per day. H2 Under-staffing: The number of nurses performing a skill on a shift must be at least equal to the minimum demand for this shift. H3 Shift type successions: A nurse cannot work certain successions of shifts on two consecutive days. H4 Missing required skill: A nurse can only cover the demand of a skill that he/she can perform. Soft constraints S1 Insufficient staffing for optimal coverage: The number of nurses performing a skill on a shift must be at least equal to an optimal demand. Each missing nurse is penalized according to a unit weight but extra nurses above the optimal value are not considered in the cost. S2 Consecutive assignments: For each nurse, the number of consecutive assignments should be within a certain range and the number of consecutive assignments to the same shift should also be within another certain range. Each extra or missing assignment is penalized by a unit weight. S3 Consecutive resting days: For each nurse, the number of consecutive resting days should be within a certain range. Each extra or missing resting day is penalized by a unit weight. S4 Preferences: Each assignment of a nurse to an undesired shift is penalized by a unit weight. S5 Complete week-end: A given subset of nurses must work both days of the week-end or none of them. If one of them works only one of the two days Saturday or Sunday, it is penalized by a unit weight. S6 Total assignments: For each nurse, the total number of assignments (worked days) scheduled in the planning horizon must be within a given range. Each extra or missing assignment is penalized by a unit weight. S7 Total working week-ends: For each nurse, the number of week-ends with at least one assignment must be less than or equal to a given limit. Each worked weekend over that limit is penalized by a unit weight. The hard constraints (Table 1, H1-H4) are typical for workforce scheduling problems: each worker is assigned an assignment or day-off every day, the demand in terms of number of employees is fulfilled, particular shift successions are forbidden, and a minimum level of qualification of the workers is guaranteed.

Soft constraints S1-S7 translate into a cost function that enhances the quality of service and retain the nurses within the unit. The quality of the schedules (alternation of work and rest periods, numbers of worked days and weekends, respect of nurses' preferences) are indeed paramount in order to retain the most qualified employees. These specificities make the NSP one of the most difficult workforce scheduling problems in the literature, because a personalized roster must be computed for each nurse. The fact that most constraints are soft eases the search for a feasible solution but makes the pursuit of optimality more difficult.

The goal of the dynamic NSP is to sequentially build weekly schedules so as to minimize the total cost of the aggregated schedule and ensure feasibility over the complete planning horizon. The main difficulty is to reach a feasible (i.e., managing the global hard constraints H3) and near-optimal (i.e., managing the global soft constraints S6 -S7 as well as consecutive constraints S2 -S3) schedule without knowing the future demands and nurses' preferences. Indeed, the hard constraints H1, H2, and H4 handle local features that do not impact the following days. Each of these constraints concern either one single day (i.e., one assignment per day H1) or one single shift (i.e., the demand for a shift H2 and the requirement that a nurse must possess a required skill H4). In the same way, soft constraints S1, and S4 -S5 are included in the objective with local costs that depend on one shift, day or weekend. To summarize, the proposed algorithm must simultaneously handle global requirements and border effects between weeks that are induced by the dynamic process. These effects are propagated to the following week/stage through the initial state or the number of worked days and weekends in the current stage.

The static nurse scheduling problem

We describe here the algorithm introduced in [START_REF] Legrain | A rotation-based branch-and-price approach for the nurse scheduling problem[END_REF] to solve the static version of the NSP. This description is important for the purpose of this paper since parts of the dynamic method described in the subsequent sections make use of certain of its specificities. This method solves the NSP with a branch-and-price algorithm [START_REF] Desaulniers | Column generation[END_REF]. The main idea is to generate a roster for each nurse, i.e., a sequence of work and rest periods covering the planning horizon. Each individual roster satisfies constraints H1, H3 and H4, and the rosters of all the nurses satisfy H2. A rotation is a list of shifts from the roster that are performed on consecutive days, and preceded and followed by a resting day; it does not contain any information about the skills performed on its shifts. A rotation is called feasible (or legal) if it respects the single assignment and succession constraints H1 and H3. A roster is therefore a sequence of rotations, separated by nonempty rest periods, to which skills are added (see Example 1). The MIP described in [START_REF] Legrain | A rotation-based branch-and-price approach for the nurse scheduling problem[END_REF] is based on the enumeration of possible rotations by column generation. As in most column-generation algorithms, a restricted master problem is solved to find the best fractional roster using a small set of rotations, and subproblems output rotations that could be added to improve the current solution or prove optimality. These subproblems are modeled as shortest path problems with resource constraints whose underlying networks are described in [START_REF] Legrain | A rotation-based branch-and-price approach for the nurse scheduling problem[END_REF]. To obtain an integer solution, this process is embedded within a branch-and-bound scheme. The remainder of the section focuses on the master problem. For the sake of clarity, we assume that, for every nurse, the set of all legal rotations is available, which conceals the role of the subproblem. It is also worth mentioning that the software is based only on open-source libraries from the COIN-OR project (BCP framework for branch-and-cut-and-price and the linear solver CLP), and is thus both free and open-source.

Example 1. Consider the following single-week roster:

We consider a set N of nurses over a planning horizon of M weeks (or K = 7M days). The sets of all shifts and skills are respectively denoted as Σ and S. The nurse's type corresponds to the set of skills he or she can use. For instance, most head nurses can fill Head Nurse demand, but they can also fill Nurse demand in most cases. All nurses of type t ∈ T (e.g., nurse or head nurse) are gathered within the subset N t . For the sake of readability, indices are standardized in the following way: nurses are denoted as i ∈ N , weeks as m ∈ {1 . . . M } , days as k ∈ {1 . . . K} , shifts as s ∈ S and skills as σ ∈ Σ. We use (k, s) to denote the shift s of day k. All other data is summarized in Table 2.

Nurses L - i , L + i min/max

Remark (Initial state). Obviously, if CR 0

i > 0, then CD 0 i = CS 0 i = 0, and vice-versa, because the nurse was either working or resting on the last day before the planning horizon. Moreover, s 0 i only matters if the nurse was working on that day. The total number of worked days and worked week-ends of a nurse is set at zero (0) at the beginning of the planning horizon.

The master problem described in Formulation (1) assigns a set of rotations to each nurse while ensuring at the same time that the rotations are compatible and the demand is filled. The cost function is shaped by the penalties of the soft constraints as no other cost is taken into account in the problem proposed by the competition. For any soft constraint SX, its associated unit weight in the objective function is denoted as c X .

Let R i be the set of all feasible rotations for nurse i. The rotation j of nurse i has a cost c ij (i.e., the sum of the soft penalties S2, S4 and S5) and is described by the following parameters:

a sk ij , a k ij ,
min i∈N j∈Ri c ij x ij S2,S4,S5 + CR + i k=1   c 3 r ik + min(K+1,k+CR + i ) l=k+1 c ikl 3 r ikl   S3 + c 6 (w + i + w - i ) S6 + c 7 v i S7 + c 1 K k=1 s∈S σ∈Σ z sk σ S1 (1a)
subject to:

[H1, H3] :

min(K+1,k+CR + i ) l=k+1 r ikl - j∈Ri:f + ij =k-1 x ij = 0, ∀i ∈ N , ∀k = 2 . . . K (1b) [H1, H3] : r ik -r i(k-1) + j∈Ri:f - ij =k x ij - k-1 l=max(1,k-CR + i ) r ilk = 0, ∀i ∈ N , ∀k = 2 . . . K (1c) [H1, H3] : K l=max(1,K+1-CR + i ) r ilK + r iK + j:f + ij =K x ij = 1, ∀i ∈ N (1d) [S6] : j∈Ri K k=1 a k ij x ij + w - i ≥ L - i , ∀i ∈ N (1e) [S6] : j∈Ri K k=1 a k ij x ij -w + i ≤ L + i , ∀i ∈ N (1f) [S7] : j∈Ri M m=1 b m ij x ij -v i ≤ B i , ∀i ∈ N (1g)
[H2] :

t∈T σ n sk tσ ≥ D sk σ , ∀s ∈ S, k ∈ {1 . . . K} , σ ∈ Σ (1h)
[S1] :

t∈T σ n sk tσ + z sk σ ≥ O sk σ , ∀s ∈ S, k ∈ {1 . . . K} , σ ∈ Σ (1i) [H4] : i∈Nt,j a sk ij x ij - σ∈Σt n sk tσ = 0, ∀s ∈ S, k ∈ {1 . . . K} , σ ∈ Σ (1j) x ij ∈ N, z sk σ , n sk tσ ∈ R, ∀i ∈ N , j ∈ R i , s ∈ S, k ∈ {1 . . . K} , t ∈ T , σ ∈ Σ (1k) r ikl , r ik , w + i , w - i , v i ≥ 0, ∀i ∈ N , k ∈ {1 . . . K} , l = k + 1 . . . min(K + 1, k + CR + i ) (1l)
where Σ t is the set of skills mastered by a nurse of type t (e.g., head nurses have the skills Head Nurse and Nurse), and T σ is the set of nurse types that masters skill σ (e.g., Head Nurse skill can be only provided by head nurses).

The objective function (1a) is composed of 5 parts: the cost of the chosen rotations in terms of consecutive assignments and preferences (S2, S4, S5), the minimum and maximum consecutive resting days violations (S3), the total number of working days violation (S6), the total number of worked week-ends violation (S7), and the insufficient staff for optimal coverage (S1). Constraints (1b)-(1d) are the flow constraints of the rostering graph (presented in Figure 1) of each nurse i ∈ N . Constraints (1e) and (1f) measure the distance between the number of worked days and the authorized number of assignments: variable w + i counts the number of missing days when the minimum number of assignments, L - i , is not reached, and w - i is the number of assignments over the maximum allowed when the total number of assignments exceeds L + i . Constraints (1g) measure the number of weekends worked exceeding the maximum B i . Constraints (1h) ensure that enough nurses with the right skill are scheduled on each shift to meet the minimal demand.

Constraints (1i) measure the number of missing nurses to reach the optimal demand. Constraints (1j) ensure a valid allocation of the skills among nurses of a same type for each shift. Constraints (1k) and (1l) ensure the integrality and the nonnegativity of the decision variables.

A valid sequence of rotations and rest periods can also be represented in a rostering graph whose arcs correspond to rotations and rest periods and whose vertices correspond to the starting days of these rotations and rest periods. Figure 1 shows an illustration of a rostering graph for some nurse i, and highlights the border effects. Nurse i has been resting for one day in her/his initial state, so the binary variable r i14 has a cost c 3 instead of zero, but the binary variable r i67 has a zero cost, because nurse i could continue to rest on the first days of the following week. If variable r i67 is set to one, nurse i will then start the following with one resting day as initial state. Finally, if nurse i was working in her/his initial state, the penalties associated to this border effect would be included in the cost of either the first rotation if the nurse continues to work, or the first resting arcs r i1k if the nurse starts to rest.

Ri1 Ri2 Ri3 Ri4 Ri5 Ri6 Ri7 Wi1 Wi2 Wi3 Wi4 Wi5 Wi6 Wi7

Handling the uncertain demand

This section concentrates on the dynamic model used for the NSP, and on the design of an efficient algorithm to compute near-optimal schedules in a very limited amount of computational time. We propose a dynamic math-heuristic based on a primal-dual algorithm [START_REF] Buchbinder | Designing Competitive Online Algorithms via a Primal-Dual Approach[END_REF] and embedded into a SAA [START_REF] Kleywegt | The sample average approximation method for stochastic discrete optimization[END_REF]. As previously stated, the dynamic algorithm should focus on the global constraints (i.e., H3, S6, and S7) to reach a feasible and near-optimal global solution.

The dynamic NSP

For the sake of clarity and because we want to focus on border effects, we introduce another model for the NSP, equivalent to Formulation [START_REF] Bäuerle | Markov Decision Processes with Applications to Finance[END_REF]. In this new formulation, weekly decisions and individual constraints are aggregated and border conditions are highlighted. The resulting weekly Formulation (2) clusters together all individual local constraints in a weekly schedule j for each week and enumerates all possible schedules.

The constraints of that model describe border effects. Although this formulation is not solved in practice, it is better-suited to lay out our online stochastic algorithm.

Binary variable y m j takes value 1 if schedule j is chosen for week m, and 0 otherwise. As for rotations, a global weekly schedule j ∈ R is described by a weekly cost c j and by parameters a ij and b ij that respectively count the number of days and weekends worked by nurse i. The variables w + i , w - i , and v i are defined as in Formulation [START_REF] Bäuerle | Markov Decision Processes with Applications to Finance[END_REF].

min j∈R M m=1 c j y m j S1-S5 + c 6 i∈N (w + i + w - i ) S6 + c 7 i∈N v i S7 (2a)
subject to:

[H1 -H4, S1 -S5] :

j∈R y m j = 1, ∀m ∈ {1 . . . M } [α m ] (2b)
[H3, S2, S3] :

j ∈Cj y m+1 j ≥ y m j , ∀j ∈ N , m = 1 . . . M -1 [δ m j ] (2c)
[S6] :

j∈R M m=1 a ij y m j + w - i ≥ L - i , ∀i ∈ N [β - i ] (2d)
[S6] :

j∈R M m=1 a ij y m j -w + i ≤ L + i , ∀i ∈ N [β + i ] (2e)
[S7] :

j∈R M m=1 b ij y m j -v i ≤ B i , ∀i ∈ N [γ i ] (2f)
y m j ∈ {0, 1}, ∀j ∈ R, m ∈ {1 . . . M } (2g) w + i , w - i , v i ≥ 0, ∀i ∈ N (2h)
The objective (2a) is decomposed into the weekly cost of the schedule and global penalties. Constraints (2b) ensure that exactly one schedule is chosen for each week. Constraints (2c) hide the succession constraints by summarizing them into a filtering constraint between consecutive schedules. These constraints simplify the resulting formulation, but will not be used in practice as their number is not tractable (see below). Constraints (2d)-(2f) measure the penalties associated with the number of worked days and weekends.

Constraints (2g)-(2h) are respectively integrality and nonnegativity constraints. The greek letters indicated between brackets (α, β, δ and γ) denote the dual variables associated with these constraints.

Constraints (2c) model the sequential aspect of the problem. This formulation is indeed solved stage by stage in practice, and thus the solution of stage m is fixed when solving stage m + 1. Therefore, when computing the schedule of stage m+1, binary variables y m j all take value zero but one of them, denoted as y m jm , that corresponds to the chosen schedule for week m and takes value 1. All constraints (2c) corresponding to y m j = 0 can be removed, and only one is kept:

j ∈Cj m y m+1 j ≥ 1
, where C jm is the set of all schedules compatible with j m , i.e., those feasible and correctly priced when schedule j m is used for setting the initial state of stage m + 1. Constraints (2c) can thus be seen as filtering constraints that hide the difficulties associated with the border effects induced by constraints H3, S2, and S3.

The main challenge of the dynamic NSP is to correctly handle constraints (2c)-(2h) to maximize the chance of building a feasible and near-optimal solution at the end of the horizon. Our dynamic procedure for generating and evaluating the computed schedules at each stage is based on the SAA. Algorithm 1 summarizes the whole iterative process over all stages. Each candidate schedule is evaluated before generating Sample a set Ω m of future demands for the evaluation while there is enough computational time do Generate a candidate weekly schedule j for stage m Initialize the evaluation algorithm with schedule j (i.e., set the initial state) for each scenario ω ∈ Ω m do Evaluate schedule j over scenario ω end for Store the schedule (S m := S m ∪ {j}) and its score (e.g., its average evaluation cost) end while Choose the schedule j m ∈ S m with the best score end for Compute the best schedule for the last stage M with the given computational time another new one. The available amount of time being short, we should not take the risk of generating several schedules without having evaluated them. This generation-evaluation step is repeated until the time limit is reached. Note that the last stage M is solved by an offline algorithm (e.g., the one described in Section 3), because the demand is totally known at this time.

The two following subsections describe each one of the main steps:

1. The generation of a schedule with an offline procedure that takes into account a rough approximation of the uncertainty; 2. The evaluation of that schedule for a demand scenario that measures the impact on the remaining weeks. (This step also computes an evaluation score of a schedule based on the sampled scenarios.)

Generating a candidate schedule

In a first attempt to generate a schedule, a primal-dual algorithm inspired from [START_REF] Buchbinder | Designing Competitive Online Algorithms via a Primal-Dual Approach[END_REF] is proposed. However, this procedure does not handle all correlations between the weekly schedules (i.e., Constraints (2c)). This primal-dual algorithm is then adapted to better take into account the border effects between weeks and make use of every available insight on the following weeks.

A primal-dual algorithm

Primal-dual algorithms for online optimization aim at building pairs of primal and dual solutions dynamically. At each stage, primal decisions are irrevocably made and the dual solution is updated so as to remain feasible. The current dual solution drives the algorithm to better primal decisions by using those dual values as multipliers in a Lagrangian relaxation. The goal is to obtain a pair of feasible primal and dual solutions that satisfy the complementary slackness property at the end of the process.

We use a similar primal-dual algorithm to solve the online problem associated to Formulation [START_REF] Bent | Scenario-based planning for partially dynamic vehicle routing with stochastic customers[END_REF]. In this dynamic process, we wish to sequentially solve a restriction of Formulation (2) to week m for all stages m ∈ {1, . . . , M } with a view to reaching an optimal solution of the complete formulation. This process raises an issue though: how can constraints (4c)-(4e) betaken into account in a restriction to a single week? To achieve that goal, the primal-dual algorithm uses dual information from stage m to compute the schedule of stage m + 1 by solving the following Lagrangian relaxation of Formulation (2):

min j∈R [ c j S1,S2,S3,S4,S5 + i∈N ( β+ i -β- i )a ij S6 + i∈N γi b ij S7 ]y m+1 j (3a) s.t.: [H1 -H3] : j∈Cj m y m+1 j = 1 (3b) y m+1 j ∈ {0, 1}, ∀j ∈ R (3c)
where βi , β+ i , γi ≥ 0 are multipliers respectively associated with constraints (2d) -(2f), and both constraints (2b) and (2c) that guarantee the feasibility of the weekly schedules are aggregated under Constraint (3b). More specifically, any new assignment for nurse i will be penalized with βi -β+ i and worked week-ends will cost an additional γi . It is thus essential to set these multipliers to values that will drive the computation of weekly schedules towards efficient schedules over the complete horizon. For this, we consider the dual of the linear relaxation of Formulation (2):

max M m=1 α m + i∈N (L - i β - i -L + i β + i -B i γ i ) (4a) s. t.: α m + i∈N (a ij β - i -a ij β + i -b ij γ i ) -δ m j + j ∈C -1 j δ m-1 j ≤ c j , ∀j ∈ R, ∀m [y m j ] (4b) β - i ≤ c 6 , ∀i ∈ N [w - i ] (4c) β + i ≤ c 6 , ∀i ∈ N [w + i ] (4d) γ i ≤ c 7 , ∀i ∈ N [v i ] (4e) β + i , β - i , γ i , δ m j ≥ 0, ∀j ∈ R, m ∈ {1 . . . M } , ∀i ∈ N (4f)
where set C -1 j contains all the schedules with which schedule j is compatible. Dual variables α m , δ m j , β - i , β + i and γ i are respectively associated with Constraints (2b), (2c), (2d), (2e), and (2f), and the variables δ 0 j are set to zero to obtain a unified formulation. The variables in brackets denote the primal variables associated to these dual constraints. At each stage, the primal-dual algorithm sets the values of the multipliers so that they correspond to a feasible and locally-optimal dual solution, and uses this solution as Lagrangian multipliers in Formulation [START_REF] Bent | Scenario-based planning for partially dynamic vehicle routing with stochastic customers[END_REF]. Another point of view is to consider the current primal solution at stage m as a basis of the simplex algorithm for the linear relaxation of Formulation [START_REF] Bent | Scenario-based planning for partially dynamic vehicle routing with stochastic customers[END_REF]. The resolution of stage m + 1 corresponds to the creation of a new basis: Formulation (3) seeks a candidate pivot with a minimum reduced cost according to the associated dual solution.

Not only does the choice of dual variables drives the solution towards dual feasibility, but it also guarantee that complementary conditions between the current primal solution at stage m and the dual solution computed for stage m + 1 are satisifed. In the computation of a dual solution, the variables α m and δ m j do not need to be explicitly considered, because they will not be used in Formulation [START_REF] Birge | Introduction to Stochastic Programming[END_REF]. What is more, focusing on stage m, the only dual constraints that involve α m and δ m j (4b), can be satisfied for any value of βi , β+ i and γi by setting δ m j = 0, ∀j ∈ R, and

α m = min j∈R {c j - i∈N (a ij β- i -a ij β+ i -b ij γi )}.
Observe that the expression of the objective function of Formulation (2) ensures that the only schedule variable satisfying y m jm > 0 will be such that

j m ∈ argmin j∈R {c j -i (a ij β - i -a ij β + i -b ij γ i )}, so comple- mentarity is achieved.
To set the values of βi , β+ i and γi , we first observe that complementary conditions are satisfied if

β- i =    c 6 if j,m a ij y m j < L - i 0 otherwise β+ i =    c 6 if j,m a ij y m j ≥ L + i 0 otherwise γi =    c 7 if j,m b ij y m j ≥ B i 0 otherwise
Since the history of the nurses are initialized with zero assignment and week-end worked, we initially set βi = c 6 and β+ i = γi = 0 to satisfy complementarity. We then perform linear updates at each stage m, using the characteristics of the schedule j m chosen for the corresponding week:

β- i = max 0, β- i -c 6 a m ijm L - i , β+ i = min c 6 , β+ i + c 6 a m ijm L + i , γi = min c 7 , γi + c 7 b m ijm B i .
These updates do not maintain complementarity at each stage but allow for a more balanced penalization of the number of assignments and worked week-ends. The variations of βi , β+ i and γi ensure that constraints (2b) remain feasible for the previous stage, even though complementarity may be lost. to be able to derive a competitive-ratio. However, no competitive-ratio is sought by this approach and linear updates are easier to design. Non-linear updates could be investigated in the future.

Algorithm 2 summarizes the primal-dual algorithm. It estimates the impact of a chosen schedule on the global soft constraints through their dual variables. As it is, it gives mixed results in practice. The reason is that the information obtained through the dual variables does not describe precisely the real problem. At the beginning of the algorithm, the value of the dual variables drives the nurses to work as much as possible.

Consequently, the nurses work too much at the beginning and cannot cover all the necessary shifts at the end of the horizon. Furthermore, the expected impact of the filtering constraints (2c) are totally ignored in that version. Namely, the shift type succession constraints H3 imply many feasibility issues at the border between two weeks as Formulation ( 2) is solved sequentially with this primal-dual algorithm. The following two sections describe how this initial implementation is adapted to cope with these issues.

Algorithm 2: Primal-dual algorithm

βi = c 6 , β+ i = γi = 0, ∀i ∈ N for each stage m do Solve Formulation (3) with a deterministic algorithm

Update β- i = max 0, β- i -c 6 a m ijm L - i , ∀i ∈ N Update β+ i = min c 6 , β+ i + c 6 a m ijm L + i , ∀i ∈ N Update γi = min c 7 , γi + c 7 b m ijm Bi , ∀i ∈ N end for

Sampling a second week demand for feasibility issues

Preliminary results have shown that Algorithm 2 raises feasibility issues due to constraints H3 on forbidden shift successions between the last day of one week and the first day of the following one. In other words, there should be some way to capture border effects during the computation of a weekly schedule. Instead of solving each stage over one week, we solve Formulation (3) over two weeks and keep only the first week as a solution of the current stage. The compatibility constraints (2c) between stages m and m + 1 are now included in this two-weeks model. In this approach, the data of the first week is available but no data of future stages is available. The demand relative to the next week is thus sampled as described in Section 4.4.

The fact that the schedule is generated for stages m and m + 1 ensures that the restriction to stage m ends with assignments that are at least compatible in this scenario, thus increasing the probability of building a feasible schedule over the complete horizon.

Furthermore, for two different samples of following week demand, the two-weeks version of Formulation (3) should lead to two different solutions for the current week. As a consequence, we can solve the model several times to generate different candidate schedules for stage m. As described in Algorithm 1, we use this property to generate new candidates until time limit is reached.

Global bounds to reduce staff shortages

Preliminary results have also shown that Algorithm 2 creates many staff shortages in the last weeks.

Our intent is thus to bound the number of assignment and worked weekends in the early stages to avoid the later shortages. The naive approach is to resize constraints (2d)-(2f) proportionally to the length of the demand considered in Formulation (3) (i.e., two weeks in our case). However, it can be desirable to allow for important variations in the number of assignments to a given nurse from one week to another, and even from one pair of weeks to another. Stated otherwise, it is not optimal to build a schedule that can only draw one or two weeks-long patterns as would be the case for less constrained environments. A simple illustration arises by considering the constraints on the maximum number of worked weekends. To comply with these constraints, no nurse should be working every weekend and, because of restricted staff availability, it is unlikely that a nurse is off every weekend. Coupled with the other constraints, this results necessarily in complex and irregular schedules. Consequently, bounding the number of assignments individually would discard valuable schedules.

Instead, we propose to bound the number of assignments and worked weekends for sets of similar nurses in order to both stabilize the total number of worked days within this set and allow irregularities in the individual schedules. We choose to cluster nurses working under the same work contract, because they share the same minimum and maximum bounds on their soft constraints. Hence, for each stage m, we add one set of constraints similar to (2d)-(2f) for each contract. In the constraints associated with contract κ ∈ Γ, the left hand-sides are resized proportionally to the number of nurses with contract κ and the number of weeks in the demand horizon. Let L m- κ , L m+ κ , and B m κ be respectively the minimum and maximum total number of assignments, and the maximum total number of worked weekends over the two-weeks demand horizon for the nurses with contract κ. We define these global bounds as:

• L m- κ = 7 * 2 M -m+1 i:κi=κ max(0, L - κi - m <m m =1 j a ij y m j ), • L m+ κ = 7 * 2 M -m+1 i:κi=κ max(0, L + κi - m <m m =1 j a ij y m j ), • B m κ = 2 M -m+1 i:κi=κ max(0, B κi - m <m m =1 j b ij y m j ),
where κ i is the contract of a nurse i.

Finally, the objective (3a) is modified to take into account the new slack variables w m- κ , w m+ κ , v m κ associated to the new soft constraints. The costs of these slack variables is set to make sure that violations of the soft constraints are not penalized more than once for an individual nurse. For instance, instead of counting the full cost c 6 for variable w m+ κ , we compute its cost as (c 6 -max i|κi=κ (β + i )). This guarantees that an extra assignment is never penalized with more than c 6 for any individual nurse. The cost of the variables w m+ κ and v m κ have been modified in the same way for analogous reasons. Formulation [START_REF] Burke | The state of the art of nurse rostering[END_REF] summarizes the final model used for the generation of the schedules. We recall that the variables y m j are now selecting a schedule j which covers a two weeks demand, and that this formulation is in fact solved by a branch-and-price algorithm that selects rotations instead of weekly schedules.

min j∈R c j + i∈N (β + i -β - i )a ij + γ i b m ij y m j + κ∈Γ (c 6 -max i:κi=κ (β - i ))w m- κ + (c 6 -max i:κi=κ (β + i ))w m+ κ ) + (c 7 -max i:κi=κ (γ m i ))v m κ ) (5a)
s.t.: [H1, H2, H3, H3] :

j∈R y m j = 1, (5b)
[S6] :

i:κi=κ j∈R a ij y m j + w m- κ ≥ L m- κ , ∀κ ∈ Γ (5c)
[S6] :

i:κi=κ j∈R a ij y m j + w m+ κ ≤ L m+ κ , ∀κ ∈ Γ (5d)
[S7] :

i:κi=κ j∈R b ij y m j -v m κ ≤ B m κ , ∀κ ∈ Γ (5e) y m j ∈ {0, 1}, ∀j ∈ R (5f) w m- κ , w m+ κ , v m κ ≥ 0, ∀κ ∈ Γ (5g)
To conclude, Formulation (5) allows to anticipate the impact of a schedule on the future through two mechanisms: the problem is solved over two weeks to diminish the border effects that may lead to infeasibility, and the costs are modified to globally limit the penalties due to constraints S6 and S7. Furthermore, this formulation can generate different schedules fort the first week by considering different samples for the second week demand.

Evaluating candidate schedules

In the spirit of the SAA, the first-week schedules generated by Formulation ( 5) are evaluated to be ranked. The evaluation should measure the expected impact of each schedule on the global solution (i.e., over M weeks). This impact can be measured by solving a NSP several times over different sampled demands for the remaining weeks.

Let Ω m be the set of scenarios of future demands for weeks m + 1, . . . , M , and assume that a schedule j has been computed for week m. To evaluate schedule j, we wish to solve the NSP for each sample of future demand ω ∈ Ω m by using j to set the initial history of the NSP. Denoting V m jω the value of the solution, we can infer that the future cost c m jω of schedule j in scenario ω is equal to c j + V m jω : the actual cost of the schedule plus the resulting cost for scenario ω. Then, a score that takes into account all the future costs (c m jω ) ω∈Ω m of a given schedule j is computed. Several functions have been tested and preliminary results have shown that the expected value was producing the best results. Finally, the schedule j m with the best score is retained.

However, computing the value V m jω raises two main issues. First, the NSP is an integer program for which it can be time-consuming to even find a feasible solution. We thus use the linear relaxation of this problem as an estimation of the future cost. This simplification decreases drastically the computational time, but still can detect feasibility issues at the border between weeks m and m + 1. The second issue is that over a long time horizon, even the linear relaxation of the NSP cannot be solved in sufficiently small computational time. We thus restrict the evaluation to scenarios of future demands that are at most two weeks long. More specifically, the scenarios are one week-long for the penultimate stage (M -1) and two weeks long for the previous stages. We observed that this restriction allows to keep the solution time short enough while giving a good measure of the impact of the schedule j on the future.

To summarize, the value V m jω is computed by solving the linear relaxation of Formulation (1) for a twoweek demand ω, and the initial state is set by using the schedule j. Finally, the parameters L - i , L + i , and B i are proportionally resized over two weeks, as follows.

• L (m+1)- i = 7 * 2 M -m max(0, L - i - m m =1 j a m ij y m j ) ; • L (m+1)+ i = 7 * 2 M -m max(0, L + i - m m =1 j a m ij y m j ) ; • B m+1 i = 2 M -m max(0, B i - m m =1 j b m ij y m j ) .
As already stated, the number of evaluation scenario included in Ω m is kept low (e.g., |Ω m | = 5) to meet the requirements in computational time. These scenarios are sampled as described in the next section.

Sampling of the scenarios

The competition data does not provide any knowledge about past demands, potential probability distributions of the demand, nor any other type of information that could help for sampling scenarios of demand.

It is thus impossible to build complex and accurate prediction models for the future demand. At a given stage m, the algorithm has absolutely no knowledge about the future realizations of the demand, so the sampling can only be based on the current and past observations of the weekly demands on stages 1 to m.

To build scenarios of future demand, we simply perturb these observations with some noise that is uniformly distributed within a small range (typically one or two nurses) and randomly mix these observations (e.g., pick the Monday of one observation and the Tuesday from another one). The future preferences are not sampled in the scenarios, because they cannot lead to an infeasible solution, they do not induce border effects, and they have small costs when compared to the other soft constraints. The goal of the sampling method is only to obtain some diversity in the scenarios used to generate different candidate schedules and in those used to evaluate the candidate schedules. Assuming that the demands will not change dramatically from one week to another, this allows for additional robustness and efficiency in many situations.

Summary of the primal-dual-based sample average approximation

Algorithm 3 provides a detailed description of the overall algorithm we submitted to the INRC-II.

It generates several schedules with a primal-dual algorithm and evaluates them over a set Ω m of future demands. The evaluation step increases the probability of selecting a globally feasible schedule, that has already been feasible for several resolutions of the linear relaxation of Formulation [START_REF] Bäuerle | Markov Decision Processes with Applications to Finance[END_REF]. The performances of this algorithm are discussed in Section 5.

Algorithm 3: A primal-dual-based sample average approximation

β - i = c 6 , β + i = γ i = 0, ∀i ∈ N for each stage m = 1 . . . M -1 do
Initialize the set of candidate schedules of stage m: S m = ∅ Initialize the generation model using the chosen schedule of the previous stage m -1 (i.e., set the initial state) Sample a set Ω m of future demands for the evaluation while there is enough computational time do Sample a second week demand for the generation Solve Formulation ( 5) with a deterministic algorithm to build a two-weeks schedule (j1 , j 2 ) Store the schedule of the first week: S m := S m ∪ {j 1 } Initialize the evaluation model with schedule j 1 (i.e., set the initial state) for each demand ω ∈ Ω m do Compute the value V m jω as the optimal value of the linear relaxation of Formulation (1) over two weeks if m < M -1, and one week otherwise end for end while Choose the schedule j m with the best average evaluation cost Update the dual variables as in Algorithm 2 end for Compute the best schedule for the last stage M with the given computational time

Experimentations

This section presents the results obtained at the INRC-II. The competition was organized in two rounds.

In the selection round, each team had to submit their best results on a benchmark of 28 instances that were available to the participants before submitting the codes. The organizers then retained the best eight teams for the final round where they tested the algorithms against a new set of 60 instances. The algorithm described above ranked second in both rounds.

The instances used during each round are summarized in Tables 3 and4. They range from relatively small instances (35 nurses over 4 weeks) to really big ones (120 nurses over 8 weeks) that are very difficult to solve even in a static setting [START_REF] Legrain | A rotation-based branch-and-price approach for the nurse scheduling problem[END_REF]. The algorithms of the participants all had the same limited computational time to solve each stage (depending on the number of nurses and on the computer used for the tests 1 ). The solution obtained at each stage was used as an initial state for the schedule of the following week. If an algorithm reached an infeasible schedule, the iterative process was stopped. The final rank of each team was computed as the average rank over all the instances. In this section, we will focus our discussion on the quality of the results obtained with our algorithm.

More details about the competition and the results can be found on the competition website: http:// mobiz.vives.be/inrc2/.

Algorithm implementation

Algorithm 3 depends on how future demands are sampled, on the number of scenarios used for the evaluation, and last but not least, on the scheduling software.

The algorithm uses only five scenarios of future demands for the evaluation. It must indeed divide the short available computational time between the generation and the evaluation of the schedules. The first step aims at computing the best schedule according to the current demand while the second step seeks a robust planning that yields promising results for the following stages (high probability to remain feasible and near-optimal). In order to generate several schedules (at least two) within the granted computational time, the number of demand scenarios must remain small. Moreover, since the demand scenarios we generate are not based on accurate data, but only on a learned distribution, there is no guarantee that a larger number of scenarios would provide a better evaluation. In fact, we tested different configurations (3 to 10 scenarios used for the evaluation), and they all gave approximately the same results (the best results were obtained for 4 to 6 scenarios).

The code is publicly shared on a Git repository [START_REF] Legrain | Dynamic nurse scheduler[END_REF]. The scheduling software is implemented in C++ and is based on the branch-cut-and-price framework of the COIN-OR project, BCP. The choice of this framework is motivated by the competition requirement of using free and open-source libraries. The pricing problems are modeled as shortest paths with resource constraints and solved with the Boost library. The solution algorithm is not parallelized and it uses the branching strategy 'two dives' described in [START_REF] Legrain | A rotation-based branch-and-price approach for the nurse scheduling problem[END_REF]. This strategy 'dives' two times in the branching tree and keeps the best feasible solution. If no solution is found after two dives (which was never the case), the algorithm continues until it reaches either the time limit or a feasible solution. (but for one third position). Their algorithm is also based on a mixed integer programming approach that computes weekly schedules, but they directly model the problem using a large flow network where a state-expansion is used to consider the soft constraints. For the time being, only a brief description of the algorithm is available in [START_REF] Römer | A direct MILP approach based on state-expanded network flows and anticipation for multistage nurse rostering under uncertainty[END_REF]. Algorithm 3 obtains a fair second position and competes with the best algorithm, since it also ranks first or second on every instance but two, for which it ranks third. Finally, the third team is significantly behind the first two. As highlighted by Figure 2, the solutions found by their algorithms exhibit at least a 9% relative gap with respect to the best solution.

Selection instances

It is also important to note that these algorithms are randomized, because they are all based on random sampling of future demands. For the first phase of the competition, the participants had to provide the random seeds they used to obtain their results. During the second phase, the organizers executed the code of each team ten times with different arbitrary random seeds on each instance. Because of these variations, we run the algorithm many times on each of the instances used for the selection to submit only the best ones and increase our chance of qualification. Most teams must have used the same technique, since the ranking between the selection and the final rounds did not really change.

Final instances

The final results respect the same ranking as the one obtained after the selection. However, these comparisons are more fair, since the results were computed by the organizers, so that the teams were not able to select the solutions they submitted. This configuration evaluates in a better way the proposed algorithms, and especially their robustness. The organizers have even run 10 times the algorithms on each of the 60 final instances, and have thus compared the proposed software over 600 tests. the best schedules for about 65% of the instances, but our algorithm appears to be more robust. We indeed able produced a feasible schedule in every test but one, while winners could not build a feasible schedule in 5% of the cases (i.e., 34 tests). This comparison also highlights the balance that needs to be found between the time spent in the generation of the best possible schedules and their evaluation, since this second phase provides a measure of their robustness to future demands.

Figure 5 shows the cumulative distribution of the relative gap of the winning team solutions from ours as a function of the number of nurses. It is clear that once the instances exceed a certain size (i.e., 110 nurses), the quality of the solution of the winning team decreases. Indeed, in [START_REF] Römer | A direct MILP approach based on state-expanded network flows and anticipation for multistage nurse rostering under uncertainty[END_REF], the winning team comments that their algorithms was simply unable to find feasible integer solutions for some week demands of these instances, showing that the method experiments difficulties in scaling up. Furthermore, this algorithm was also not able to find feasible solutions for an important number of the small instances (i.e., 35 nurses). As a possible explanation, we have observed that it is more difficult to find feasible solutions for these instances, because they leave less flexibility for creating the schedules, i.e., because the hard constraints of the MIP are tighter. Stated otherwise, the proportion of feasible schedules that meet the minimum demand is much smaller for the smallest instances used in the competition.

Conclusions

This article deals with the nurse scheduling problem as described in the context of the international competition INRC-II. The objective is to sequentially build, week by week, the schedule of a set of nurses over a planning horizon of several weeks. In this dynamic process, the schedule computed for a given week is irrevocably fixed before the demand and the preferences for the next week are revealed. The main difficulty is to correctly handle the border effects between weeks and the global soft constraints to compute a feasible and near-optimal schedule for the whole horizon.

Our main contribution is the design of a robust online stochastic algorithm that performs very well over a wide range of instances (from 30 to 120 nurses over a four or eight weeks horizon). The proposed algorithm embeds a primal-dual algorithm within a sample average approximation. The primal-dual procedure generates candidate schedules for the current week, and the sample average approximation allows to evaluate each of them and retain the best one. The resulting implementation is shared on a public repository [START_REF] Legrain | Dynamic nurse scheduler[END_REF] and builds upon an open source static nurse scheduling software.

The designed algorithm won second prize in the INRC-II competition. The results show that, although this procedure does not compute the best schedules for a majority of instance, it is the most robust one.

Indeed, it finds feasible solutions for almost every instance of the competition while providing high-quality schedules.

T 1 Figure 1 :

 11 Figure 1: Example of a rostering graph for nurse i ∈ N over a horizon of K = 7 days, where the minimum and maximum number of consecutive resting days are respectively CR - i = 2 and CR + i = 3, and the initial number of consecutive resting days is CR 0 i = 1. The rotation arcs (x ij ) are the plain arrows, the rest arcs (r ikl and r ik ) are the dotted arcs, and the artificial flow arcs are the dashed arrows. The bold rest arcs have a cost c 3 of and the others are free.

Algorithm 1 :

 1 A sample average approximation based algorithm for each stage m = 1 . . . M -1 do Initialize the set of candidate schedules of stage m: S m = ∅ Initialize the generation algorithm with the chosen schedule j m-1 of the previous stage m -1 (i.e., set the initial state)

Figure 2 :

 2 Figure 2: Cumulative distribution of the relative gap on the selection instances

Figure 3 :

 3 Figure 3: Distribution of the objective value for the selection instances

Figure 3

 3 Figure 3 shows the distribution of the objective value for 180 observations of the solution of an instance with 80 nurses over 8 weeks using Algorithm 3. The values of the solution are within a [-6%, +7%] range.

Figure 4 :

 4 Figure 4: Cumulative distribution of the relative gap on the final instances

Figure 4

 4 Figure 4 presents the relative gaps obtained on the final instances. The two first teams are really close and their algorithms highlight two distinct features of the competition. The winning team's algorithm builds

Figure 5 :

 5 Figure 5: Cumulative distribution of the relative gap as a function of the number of nurses

Table 1 :

 1 Constraints handled by the software.

Table 2 :

 2 Summary of the input data.

	Demand
	D sk σ O sk σ	min demand in nurses performing skill σ on shift (k, s) optimal demand in nurses performing skill σ on shift (k, s)
	Initial state
	CD 0 i CS 0 i s 0 i CR 0 i	initial number of ongoing consecutive worked days for nurse i initial number of ongoing consecutive worked days on the same shift for nurse i shift worked on the last day before the planning horizon for nurse i initial number of ongoing consecutive resting days for nurse i

total number of worked days over the planning horizon for nurse i CR - i , CR + i min/max number of consecutive days-off for nurse i B i max number of worked week-ends over the planning horizon for nurse i

  + ij represent the first and last worked days of this rotation. Let x ij be a binary decision variable which takes value 1 if rotation j is part of the schedule of nurse i and zero otherwise. The binary variables r ikl and r ik measure if constraint S3 is violated: they are respectively equal to 1 if nurse i has a rest period from day k to l -1 including at most CR +

and b m ij which are equal to 1 if nurse i works respectively on shift (k, s), on day k, and weekend m, and 0 otherwise. Finally, f - ij and f i consecutive days (cost: c ikl 3 ), and if nurse i rests on day k and has already rested for at least CR + i consecutive days before k, and to zero otherwise. The integer variables w + i and w - i count the number of days worked respectively above L + i and below L - i by nurse i. The integer variable v i counts the number of weekends worked above B i by nurse i. Finally, the integer variables n sk σ , n sk tσ , and z sk σ respectively measures the number of nurses performing skill σ, the number of nurses of type t performing skill σ, and the undercoverage of skill σ on shift (k, s).

Table 3 :

 3 Instances used for the selection

	Number of nurses
	35	70	110

Table 4 :

 4 Instances used for the final

The computational times are given for a linux computer Intel(R) Xeon(R) X5675 @ 3.07GHz with 8 Go of available memory.

Despite the limits of this algorithm, our intent with this article is to present the exact implementation that has been submitted to the competition. There is place for improvements that could be developed in the future. For instance, the primal-dual algorithm could be enhanced with non-linear updates, new features recently developed in the static nurse scheduling software could be tested, or the bounding constraints added in the primal-dual algorithm could be refined.


A Mondiere

V Déneux

N Binot

D Delagnes

Controlling the MC and M 2 C carbide precipitation in Ferrium® M54® steel to achieve optimum ultimate tensile strength/fracture toughness balance

Keywords:

Controlling the MC and M2C carbide precipitation in

Ferrium® M54® steel to achieve optimum ultimate tensile strength/fracture toughness balance

Aurélien Mondière, Valentine Déneux, Nicolas Binot, Denis Delagnes

Introduction

Aircraft applications, particularly for landing gear, require steels with high mechanical resistance, fracture toughness and stress corrosion cracking resistance [START_REF] Flower | High Performance Materials in Aerospace[END_REF]. Additionally, the aerospace industry is looking for different ways to reduce the weight of landing gear parts, as the landing gear assembly can represent up to 7% of the total weight of the aircraft [START_REF] Kundu | Aircraft Design[END_REF]. The search for metal alloys with a better balance of mechanical properties while maintaining a constant production cost is stimulating research activities. For several decades, 300 M steel has been widely used for landing gear applications. However, its fracture toughness and stress corrosion cracking resistance need to be improved and aeronautical equipment suppliers are searching for new grades. As shown in Fig. 1, AerMet® 100 and Ferrium® M54® (M54®) grades are excellent candidates to replace the 300 M steels without any reduction in strength or increase in weight. Other grades do not present a high enough fracture toughness, or are not resistant enough.

The recent development of M54® steel since 2010 [START_REF] Jou | Lower-Cost, Ultra-High-Strength, High-Toughness Steel[END_REF] has led to a higher stress corrosion cracking resistance and lower cost due to its lower cobalt content (see Table 1), as compared to the equivalent properties of the AerMet® 100 grade. These two steels belong to the UHS Co-Ni steel family.

UHS Co-Ni steels were developed at the end of 1960s with the HP9-4-X [START_REF] Garrison | Ultrahigh-strength steels for aerospace applications[END_REF] and HY-180 [START_REF] Dabkowski | Nickel, Cobalt, Chromium Steel[END_REF] grades, with the main goal being to achieve higher fracture toughness than 300M or 4340 steels. The main idea was first to replace cementite by M 2 C alloy carbide precipitation during tempering to avoid brittle fracture without too large reduction in mechanical strength. A better balance of UTS/K 1C was achieved with AF1410 [START_REF] Little | High Strength Fracture Resistant Weldable Steels[END_REF] by increasing the content of carbide-forming elements. In addition, an improvement in fracture toughness was also requested and finally achieved by the accurate control of reverted austenite precipitation during tempering [START_REF] Haidemenopoulos | Dispersed-Phase Transformation Toughening in UltraHigh-Strength Steels[END_REF] and the addition of rare earth elements to change the sulfide type [START_REF] Handerhan | A comparison of the fracture behavior of two heats of the secondary hardening steel AF1410[END_REF][START_REF] Handerhan | Effects of rare earth additions on the mechanical properties of the secondary hardening steel AF1410[END_REF], resulting in an increase in inclusion spacing [START_REF] Garrison | Lanthanum additions and the toughness of ultra-high strength steels and the determination of appropriate lanthanum additions[END_REF]. Thus, AerMet® 100 was patented in 1993 [START_REF] Hemphill | High Strength, High Fracture Toughness Alloy[END_REF], incorporating these scientific progress to achieve the same strength level of 300 M but with a higher fracture toughness. Then, from the 1990s to the 2000s, scientists sought to improve the grain boundary cohesion to further increase the fracture toughness by W, Re and B additions [START_REF] Kantner | Designing Strength, Toughness, and Hydrogen Resistance: Quantum Steel[END_REF]. Thus, Ferrium® S53® steel, developed in 2007 [START_REF] Kuehmann | Nanocarbide Precipitation Strengthened Ultrahigh-Strength[END_REF], was the first steel of the family containing W. Seven years ago, Ferrium® M54® steel was designed, offering a steel with roughly the same mechanical properties as AerMet® 100, but with a lower price thanks to a lower cobalt content.

UHS Co-Ni steels all exhibit an excellent UTS/K 1C balance due to a M 2 C carbide precipitation during tempering in a highly dislocated lathmartensitic matrix [START_REF] Jou | Lower-Cost, Ultra-High-Strength, High-Toughness Steel[END_REF][START_REF] Little | High Strength Fracture Resistant Weldable Steels[END_REF][START_REF] Hemphill | High Strength, High Fracture Toughness Alloy[END_REF][START_REF] Ayer | Transmission electron microscopy examination of hardening and toughening phenomena in Aermet 100[END_REF][START_REF] Olson | APFIM study of multicomponent M2C carbide precipitation in AF1410 steel[END_REF][START_REF] Machmeier | Development of a strong (1650MNm -2 tensile strength) martensitic steel having good fracture toughness[END_REF]. However, there is limited literature on the recently developed M54® steel [START_REF] Wang | Austenite layer and precipitation in high Co-Ni maraging steel[END_REF][START_REF] Wang | Analysis of fracture toughness in high Co-Ni secondary hardening steel using FEM[END_REF][START_REF] Lee | Ferrium M54 Steel[END_REF][START_REF] Pioszak | Hydrogen Assisted Cracking of Ultra-High Strength Steels[END_REF].

The addition of alloying elements in UHS Co-Ni steels also forms stable carbides like M 6 C or M 23 C 6 during the heat treatment process.

The size of these stable carbides can easily reach several 100 nm, resulting in a significant decrease in fracture toughness by acting as microvoid nucleation sites during the mechanical load [START_REF] Schmidt | Solution treatment effects in AF1410 steel[END_REF]. These particles can be dissolved by increasing the austenitizing temperature, but the prior austenite grain size rapidly increases and induces a detrimental effect on the mechanical properties [START_REF] Sankaran | Metallurgy and Design of Alloys with Hierarchical Microstructures[END_REF]. The new challenge for these steels is thus to dissolve coarse stable carbides without an excessive grain growth.

This challenge is also well-known in other kinds of martensitic steels for other applications, such as hot work tool steels. Michaud [START_REF] Michaud | The effect of the addition of alloying elements on carbide precipitation and mechanical properties in 5% chromium martensitic steels[END_REF] showed that V-rich carbide precipitation during tempering achieves high mechanical properties at room temperature as well as at high temperature. However, precipitation stayed heterogeneously distributed in the matrix, regardless of the austenitizing and tempering conditions, and so fracture toughness and Charpy impact were limited. Indeed, the same V-rich precipitation (MC type) that controls the austenitic grain size during austenitizing and controls the strength during tempering were identified. The incomplete solutionizing of V-rich carbides during austenitizing does not permit a homogeneous concentration of alloying elements in the martensitic matrix after quench, which explains why the strength/fracture toughness balance is limited. The generic idea would be to introduce a double/different precipitation with a single and precise role for each population: to control the austenitic grain size OR to control the mechanical strength. In H11-type tool steels, the addition of Mo slightly improved the balance of properties [START_REF] Michaud | The effect of the addition of alloying elements on carbide precipitation and mechanical properties in 5% chromium martensitic steels[END_REF].

In steels for aircraft applications, Olson [START_REF] Olson | Overview: Science of Steel[END_REF] and Gore et al. [START_REF] Gore | Grain-refining dispersions and properties in ultrahigh-strength steels[END_REF] succeeded in introducing another type of homogeneous small particles which pin the grain even for elevated austenitization temperatures (T = 1200 °C) in AF1410: (Ti,Mo)(C,N). These carbides avoid grain coarsening between 815 °C and 885 °C at austenitization leading to an increase in fracture toughness due to coarse carbides dissolution [START_REF] Schmidt | Solution treatment effects in AF1410 steel[END_REF]. The patent of Ferrium® S53® steel also describes a nanoscale MC precipitation which pins the grain boundary and avoids grain coarsening by the dissolution of the coarse carbides [START_REF] Kuehmann | Nanocarbide Precipitation Strengthened Ultrahigh-Strength[END_REF].

Stable carbide dissolution in Ferrium® M54® seems to be particularly challenging due to the formation of both M 2 C and M 6 C Mo-rich carbides during the heat treatment process (see Fig. 2). Indeed, as Morich M 2 C carbides precipitate during tempering, the full dissolution of Mo-rich carbides is needed to achieve a homogeneous distribution of Mo within the matrix.

More specifically, particles that control the austenitic grain size need to be stable enough at high temperature to dissolve the whole population of M 2 C and M 6 C carbides without grain coarsening. The aim of this article is to investigate carbides precipitation in M54® after a cryogenic treatment following the quench as well as after tempering. Carbide distribution, size and composition are carefully described for both states.

Experiments

Materials and Heat Treatment

Specimens were taken at mid-radius of a single bar of diameter 10.25 cm in the longitudinal direction.

The performed heat treatments were in agreement with the QuesTek recommendations [27] and consisted of a preheating treatment at 315 °C/1 h, a solutionizing at 1060 °C/1 h, followed by an oil quench, cold treatment at -76 °C/2 h and tempering at 516 °C/10 h.

Experimental Techniques

Austenite grain size was measured after the quench. Precipitation in the quenched state, after cryogenic treatment, was observed to identify undissolved carbides. Secondary carbides were characterized after tempering, at the end of the whole heat treatment process.

Chemical composition of the alloy was measured with a Q4 Tasman Spark Optical Emission Spectrometer from Bruker.

Dilatometry was performed using a Netzsch apparatus, DIL402C. Samples for dilatometry were in the form of a cylinder of diameter 3,7 mm with a length of 25 mm. Samples were heated at 7 °C/min and cooled at 5 °C/min under argon atmosphere.

For the as-quenched state, carbides were extracted by chemical dissolution of the matrix with a modified Berzelius solution at room temperature [START_REF] Burke | Chemical extraction of refractory inclusions from iron-and nickel-base alloys[END_REF] as already developed by Cabrol et al. [START_REF] Cabrol | Experimental investigation and thermodynamic modeling of molybdenum and vanadium-containing carbide hardened iron-based alloys[END_REF]. At the end of the dissolution, the solution was centrifuged to collect nanoscale precipitates. A Beckman Coulter Avanti J-30I centrifugal machine equipped with a JA-30.50Ti rotor was used to centrifuge the solution. The experimental method is described precisely in [START_REF] Cabrol | Experimental investigation and thermodynamic modeling of molybdenum and vanadium-containing carbide hardened iron-based alloys[END_REF].

XRD characterizations of the powder obtained after the chemical dissolution and of the bulk sample were performed using a Panalytical X'Pert PRO diffractometer equipped respectively with a Cu or Co radiation source. Phase identification was achieved by comparing the diffraction pattern of the experimental samples with reference JCPDS patterns.

Prior austenite grain size measurement is difficult because of the very low impurity content in the grade M54®. An oxidation etching was conducted by heating polished samples in a furnace at a temperature of 900 °C and 1100 °C under room atmosphere for 1 h and slightly polishing them after quenching to remove the oxide layer inside the grains and keeping the oxide only at the grain boundary.

Transmission Electron Microscopy (TEM) observations were performed using a JEOL JEM 2100F. Thin foils for TEM were cut from the specimens and the thickness was reduced to approximately 150 μm. Then, they were cut into disks and polished to a thickness of about 60 μm. The thin foils were then electropolished in a perchloric acidmethanol solution at -15 °C with a TenuPol device.

Chemical composition at nanometer scale was determined using atom probe tomography (APT) at the Northwestern University Center for Atom-Probe Tomography (NUCAPT). Samples were prepared into rods with a cross section of 1 x 1mm 2 and electro-polished using a twostep process at room temperature [START_REF] Krakauer | Systematic procedures for atom-probe field-ion microscopy studies of grain boundary segregation[END_REF][START_REF] Krakauer | A system for systematically preparing atom-probe field-ion-microscope specimens for the study of internal interfaces[END_REF]. The APT analyses were conducted with a LEAP 4000X-Si from Cameca at a base temperature of -220 °C, a pulse energy of 30pJ, a pulse repetition rate of 250 kHz, and an ion detection rate of 0.3% to 2%. This instrument uses a localelectrode and laser pulsing with a picosecond 355 nm wavelength ultraviolet laser, which minimizes specimen fracture [START_REF] Bunton | Advances in pulsed-laser atom probe: instrument and specimen design for optimum performance[END_REF].

For the prediction of the different types and molar fraction of each phase according to temperatures, thermodynamics calculations were performed using ThermoCalc® software. This software and database were developed at the Royal Institute of Technology (KTH) in Stockholm [START_REF] Sundman | The Thermo-Calc databank system[END_REF]. ThermoCalc® calculations were performed using the TCFE3 database.

Results and Discussions

Discussion of Optimized Mechanical Properties With Finely Dispersed Nanometer Size Precipitation

Research activities on UHS steels for aircraft applications focus on maximizing mechanical strength without decreasing the fracture toughness and stress corrosion cracking resistance. To improve strength, dislocations mobility must be reduced. Consequently, increasing the number density of secondary particles (Np) is a wellknown method and the resulting hardening is given by the following equation [START_REF] Sankaran | Metallurgy and Design of Alloys with Hierarchical Microstructures[END_REF]:

∆ ≈ σ Gb f d ( ) P 0.5 (1)
where Δσ p is particle contribution to the yield strength, G is the shear modulus, d is the particle diameter, f the volume fraction of the particle and b the Burgers vector of dislocations. Indeed, for the same volume fraction, a small particles distribution leads to a better yield strength, due to the decrease in dislocation mobility.

To obtain this fine and dispersed precipitation, two different types of nucleation are generally observed to occur in UHS steels:

-Numerous preferential nucleation sites leading to heterogeneous nucleation; -Homogeneous supersaturation of carbide-forming elements.

For the first condition, the heterogeneous nucleation of M 2 C carbides on dislocations has already been observed in previous works [START_REF] Speich | Strength and toughness of Fe-10ni alloys containing C, Cr[END_REF][START_REF] Kuehmann | Thermal Processing Optimization of Nickel-Cobalt Ultrahigh-Strength Steels[END_REF]. Indeed, dislocation sites are energetically favorable due to atom segregation and the short diffusion path offered to the diffusing element (pipe diffusion). It is therefore important to maintain a high dislocation density during tempering. Consequently, cobalt is added to these alloys to keep a high dislocation density during tempering. As previously described in the literature [START_REF] Kantner | Designing Strength, Toughness, and Hydrogen Resistance: Quantum Steel[END_REF][START_REF] Olson | Overview: Science of Steel[END_REF], Co delays the dislocation recovery through the creation of short-range ordering (SRO) in the matrix. Co also decreases the solubility of Mo in ferrite and increases the carbon activity inside ferrite [START_REF] Speich | Strength and toughness of Fe-10ni alloys containing C, Cr[END_REF][START_REF] Rhoads | High strength, high fatigue structural steel[END_REF][START_REF] Honeycombe | Steels microstructure and properties[END_REF][START_REF] Speich | Tempering of steel[END_REF][START_REF] Delagnes | Cementite-free martensitic steels: a new route to develop high strength/high toughness grades by modifying the conventional precipitation sequence during tempering[END_REF], leading to a more intensive precipitation of M 2 C carbides.

The main criterion for accessing the second condition is related to the dissolution carbides during austenitizing. If carbides are not totally solutionized, the precipitation during tempering will be heterogeneously dispersed with a higher density of clusters in the areas of high concentration of the carbide-forming elements. To avoid heterogeneous concentration, remaining carbides from the previous stage of heat treatment should be totally dissolved and enough time should be spent at a temperature above the carbide solvus to obtain a homogeneous composition of the carbide-forming elements in austenite. Moreover, in order to obtain a fine and dispersed precipitation during tempering, the driving force must be increased by increasing the supersaturation resulting in a higher nucleation rate [START_REF] Kuehmann | Thermal Processing Optimization of Nickel-Cobalt Ultrahigh-Strength Steels[END_REF]. Furthermore, undissolved carbides also reduce the potential volume fraction of particles that may precipitate during tempering [START_REF] Sato | Improving the Toughness of Ultrahigh Strength Steel[END_REF] and almost total dissolution is needed. Thus, the austenitizing condition should be rationalized based on the carbide dissolution kinetics and diffusion coefficient of alloying elements in the matrix to obtain a homogeneous chemical composition of the carbide-forming elements in the martensitic matrix in the as-quenched state.

Identification of Carbide Solutionizing Temperature

The temperatures of phase transformation were determined by dilatometry experiments. According to the relative length change shown in Fig. 3(a), Ac 1 , Ac 3 and M s temperatures are clearly detected. To detect the solutionizing of carbides, the derivative of the relative length change was calculated. Carbide dissolution takes place at a temperature ranging from 970 °C to 1020 °C, as shown in Fig. 3(b).

If the austenitizing temperature is not high enough, undissolved carbides are clearly observed (see Fig. 4) and slightly decrease UTS from 1997 MPa at 1060 °C to 1982 MPa at 1020 °C, which is probably due to the carbon trapped inside those undissolved particles.

These coarse carbides can also be observed after polishing and a Nital 2% etch using SEM (Fig. 5). The volume fraction seems to be particularly high.

According to ThermoCalc® calculations, these undissolved carbides obtained after 1 h at 1020 °C are M 6 C carbides (see Fig. 2) containing a significant amount of W (see Table 2).

The high solutionizing temperature of the M54® steel as compared to other steels of the same family (free of W, see Table 3) is due to the tungsten addition which stabilizes the M 6 C carbides. If the austenitizing temperature is not high enough, undissolved carbides still remain (see Fig. 4 and Fig. 5) and the tensile properties (yield strength, UTS, elongation at rupture), as well as fatigue resistance are reduced. However, if the austenitizing temperature is too high and no carbides remain, a huge grain size coarsening can be observed also leading to a decrease in the usual mechanical properties.

According to Naylor and Blondeau [START_REF] Naylor | The respective roles of the packet size and the lath width on toughness[END_REF], thinner laths and lath packets, directly dependent on austenite grain size [START_REF] Sankaran | Metallurgy and Design of Alloys with Hierarchical Microstructures[END_REF], can improve fracture toughness by giving a long and winding route to the crack during rupture. Białobrzeska et al. [START_REF] Białobrzeska | The influence of austenite grain size on the mechanical properties of low-alloy steel with boron[END_REF] have clearly shown that at room temperature, strength, yield strength, fatigue resistance and impact energy increase when the average austenite grain size decreases. Thus, any coarsening of austenite grains should be avoided.

Pinning of the Grain Boundary and Chemical Homogenization of the Austenitic Matrix at 1060 °C

As previously mentioned in the introduction, to control the grain size during austenitizing without any impact on precipitation during tempering, the precipitation of two types of particles is needed: one type to control the grain size during solutionizing and the second type of particles which precipitates during tempering.

To achieve this goal, one way is to add MC type precipitation to avoid quick coarsening of austenitic grains. However, according to ThermoCalc® calculations, the MC solvus temperature is not sufficiently high to allow the total dissolution of M 6 C carbides (see Fig. 2). Thus, Olson [START_REF] Olson | Overview: Science of Steel[END_REF] and Gore [START_REF] Gore | Grain-refining dispersions and properties in ultrahigh-strength steels[END_REF] added some Ti to form more stable MC carbides and dissolve other coarse stable carbides. A little addition of Titanium is sufficient to obtain a significant effect on the grain size, as described by Kantner who adds 0.04%mass [START_REF] Kantner | Designing Strength, Toughness, and Hydrogen Resistance: Quantum Steel[END_REF] of titanium in Fe- 15Co-6Ni-3Cr-1.7Mo-2 W-0.25C and Fe-15Co-5Ni-3Cr-2.7Re-1.2 W-0.18C steels, or Lippard who adds only 0.01%mass [START_REF] Lippard | Microanalytical Investigations of Transformation[END_REF] in alloys AF1410, AerMet® 100, MTL2 and MTL3. A low volume fraction of thin particles seems to be efficient in preventing austenitic grain growth [START_REF] Gore | Grain-refining dispersions and properties in ultrahigh-strength steels[END_REF]. Indeed, an addition of 0.01%mass of Ti in the M54® grade is enough to shift the MC solvus temperature by approximately 100 °C above the MC solvus temperature of the M54® grade free of Titanium according to ThermoCalc® calculation (see Fig. 6).

Moreover, MC carbides contain a large amount of Ti (see Fig. 7) which is not the case for M 2 C precipitation during tempering. Consequently, Ti-rich MC carbides seem relevant, to be a solution to control the grain size without any impact on precipitation during tempering. The purpose of the following paragraph is to compare the experimental results with the above-mentioned theoretical prediction.

After austenitizing for 1 h at 1060 °C, fine undissolved carbides were found in the as-quenched state after cryogenic treatment in M54® steels. These carbides are thinner than the undissolved carbides observed. In addition, a lower volume fraction is measured after an austenitization at 1060 °C than after a 1020 °C or 920 °C austenitization (see Fig. 8). The average size of these carbides is around 70 nm, measured on a sample of 23 carbides sample. In addition, no coarse undissolved carbides are observed indicating that the optimal austenitization conditions are not far to be reached.

Chemical extraction of carbides in the as-quenched state was performed to determine the type of those undissolved carbides still remaining after a 1060 °C austenitizing. As predicted by the thermocalc calculation, a FCC structure (type MC) was clearly identified from the XRD patterns (see Fig. 9). Moreover, the chemical composition measured by EDX (Energy Dispersive X-ray spectroscopy) is (Ti 0.44 Mo 0.27 W 0.13 V 0.16 )C. This composition is in quite good agreement with the ThermoCalc® calculated composition (Ti 0.55 V 0.25 Mo 0.17 W 0.08 )C 0.95 .

According to Spark Optical Emission Spectrometer measurements, the average Ti concentration measured is about 0.013 wt% in M54® steel. Considering that all the Ti atoms precipitate and taking into account the chemical composition of the MC measured by EDX, the volume fraction of Ti-rich MC carbide is found to be nearly 0.06%.

The intercarbide distance can be estimated using the equation given by Daigne et al. [START_REF] Daigne | The influence of lath boundaries and carbide distribution on the yield strength of 0.4% C tempered martensitic steels[END_REF]:

= × d r π f 1.18 2 3 particle v ( 2
)
where d is the distance between particles, r is the radius of the particle and f v the volume fraction of particles. According to Eq. ( 2), the distance between the MC carbides with an addition of 0.013 wt% of Ti is about a micrometer. This value is in very good agreement with SEM observations (see Fig. 8) indicating that most of the titanium carbides remain undissolved after the austenitization at 1060 °C.

Furthermore, a relation has been developed in tool steels to describe the grain refinement by a particle dispersion in tools steels. Bate [START_REF] Bate | The effect of deformation on grain growth in Zener pinned systems[END_REF] suggested the following equation between the limiting grain size diameter D, the mean radius, r, and the volume fraction F v of the pinning particles:

= D r F 4 3 v (3)
The calculated average grain size diameter is 78 μm according to the Bate's Eq. ( 3) in M54®.

This value is in a very good agreement with the measured average grain size of 81 ± 39 μm at 900 °C or 79 ± 38 μm at 1100 °C (see Fig. 10). Approximately 300 grains were measured for each austenitizing temperature. According to the Bate's work, the estimated 0.06% volume fraction of undissolved MC carbides is sufficient to control the grain size of austenite.

Consequently, MC particles need only a very small quantity of carbide-forming elements required for M 2 C precipitation during tempering. In addition, the calculated diffusion lengths of the different carbide-forming elements, Mo, Cr, W, are clearly significantly higher than the distance between first neighbors of Mo, Cr, W, respectively in the austenitic matrix at the end of austenitization (1060 °C/1 h) (see Table 4). As a consequence, homogeneous composition of the austenite is quickly obtained before quenching.

By way of conclusion, a small amount of Ti-rich MC carbides control the austenitic grain size and above all, the complete dissolution of M 6 C molybdenum rich carbides leads to the homogeneous distribution of the M 2 C carbide-forming elements before quenching.

Precipitation During Tempering

The particles that precipitate during tempering are totally different from the carbides controlling the austenitic grain size. According to XRD results, M 2 C-type carbides are identified after a tempering for 500 h at 516 °C (see Fig. 11). This long duration of tempering is necessary to detect the diffraction peaks of M 2 C carbides. For the standard tempering of 10 h, the volume fraction and the size of carbides might be too low to be detected by XRD, or long-distance ordering of M 2 C carbides (hexagonal structure) might not be achieved as already suggested by Machmeier et al. [START_REF] Ayer | On the characteristics of M2C carbides in the peak hardening regime of AerMet 100 steel[END_REF].

Consequently, the same carbide type is identified in the M54®, AerMet® 100 and AF1410 steels [START_REF] Ayer | Transmission electron microscopy examination of hardening and toughening phenomena in Aermet 100[END_REF][START_REF] Ayer | Microstructural basis for the effect of chromium on the strength and toughness of AF1410-based high performance steels[END_REF]. Atom probe analyses were performed to determine the distribution of M 2 C carbides within the martensitic matrix and to estimate the chemical composition of M 2 C carbides. To define the particle/matrix interface found in the analyzed box, the adopted criterion is an isoconcentration of 36 at% of Molybdenum + Carbon. Carbides seem to be homogenously distributed within the matrix according to the (limited) volume analyzed by APT (see Fig. 12). According to TEM observations, the precipitation of M 2 C carbides during tempering is very fine with an average size of 9.6 × 1.2 nm measured on 130 carbides (see Fig. 13) and seems to be homogeneously distributed within the matrix, as already shown by APT. The shape of the M 2 C particles is very elongated with an aspect ratio near 10. The main conclusion can be summarized as follows: the 1060 °C austenitizing temperature contributes to a fine and dispersed precipitation of M 2 C carbides after tempering, thanks to a high supersaturation as well as a homogeneous distribution of carbide-forming elements.

The average chemical composition of the M 2 C carbides measured by atom probe is Mo-rich with a significant content of Cr, W and V (see Fig. 14).

The chemical composition of M 2 C measured by atom probe is in quite good agreement with the ThermoCalc® calculations (see Fig. 15). The M 2 C carbides contain mainly Mo and Cr with approx. 10% W and a small amount of Fe and V, as shown in Fig. 15 and Table 5.

However, the chemical composition of the carbides in M54® is quite different from the composition measured in AerMet® 100 and AF1410 steels (see Table 5). Indeed, the main difference comes from the W content in M 2 C carbides for the M54® steel. W has a slower diffusivity than other carbide-forming elements and stabilizes M 2 C carbides for long duration tempering [START_REF] Lee | Stability and coarsening resistance of M2C carbides in secondary hardening steels[END_REF] which guarantee the mechanical properties in a wide range of tempering condition. Moreover, very few cementite precipitates are observed in the M54® steels. This fact also contributes to the high fracture toughness value measured after tempering. Indeed, cementite is well known to strongly reduce the fracture toughness of high strength steels [START_REF] Speich | Strength and toughness of Fe-10ni alloys containing C, Cr[END_REF], particularly if the iron carbide is located at the interlath site. The W in M 2 C carbides allows a long duration of the tempering treatment resulting in the total dissolution of cementite without coarsening of M 2 C carbides.

Conclusion

Ferrium® M54® steel was developed by QuesTek using intensive thermodynamic calculations [START_REF] Olson | Materials genomics: from CALPHAD to flight[END_REF]. An excellent strength/fracture toughness balance is achieved with a UTS reaching 1965 MPa and a K 1C values up to 110 MPa√m.The main goal of this work is to provide experimental evidence and arguments explaining the outstanding UTS/ K1C balance of properties the work is focused on the precipitation identification during the heat treatment by a different scale microstructural study using advanced experimental tools (XRD, TEM, APT). To this end, the optimization of austenitizing conditions is of primary importance, in conjunction with the solutionizing of alloying elements needed for precipitation during tempering. The main results can be summarized as follows:

▪ Microstructure in the as-quenched state (after cryogenic treatment) can be defined as a Ti-rich MC carbide precipitation with a size from 50 nm to 120 nm in a martensitic matrix which is highly supersaturated in carbide-forming elements. In addition, those elements are homogeneously distributed within the matrix, according to length-diffusion calculations. ▪ The addition of small amount of titanium has led to full dissolution of the Mo-and W-rich carbides. Types of precipitates which control the grain size during the austenitization and which strengthen the steel during the tempering are then totally different. ▪ This final microstructure is obtained thanks to the proper solutionizing of alloying elements during austenitizing at high temperature (1060 °C) which results in: o A high supersaturation before tempering. o A homogeneously distributed nucleation of carbides. ▪ Microstructure in the tempered state 516 °C/10 h is characterized by a homogeneously distributed precipitation of nanometer-sized M 2 C carbides. These carbides contain W, which reduces their coarsening rate.

Table 5

Comparison of carbide composition of different UHS steels hardened by M 2 C carbide precipitation according to ThermoCalc® calculations and experimental measurements.

Fig. 1 .

 1 Fig. 1. Comparison of different grades of steel according to their fracture toughness, ultimate tensile strength and stress corrosion cracking resistance (adapted from [3]).

Fig. 3 .

 3 Fig. 3. Relative length change curve (a) and derivative of the relative change curves (b) obtained from dilatometer heating experiments.

Fig. 4 .

 4 Fig. 4. SEM image of a fracture surface of a tensile specimen (austenitization performed at 1020 °C).

Fig. 5 .

 5 Fig. 5. SEM image of an as-quenched sample austenitized at 920 °C after nital etch.

Fig. 6 .

 6 Fig. 6. Mole fraction of phase according to austenitizing temperature in M54® with and without 0.01%mass Ti calculated with TCFE3 ThermoCalc® database.

Fig. 7 .

 7 Fig. 7. Composition of MC carbide according to the temperature calculated with TCFE3 ThermoCalc® database.

Fig. 8 .

 8 Fig. 8. SEM observations of undissolved carbides after 1060 °C austenitizing and Nital etch (as-quenched structure).

Fig. 9 .

 9 Fig. 9. Pattern and experimental XRD profiles (relative intensities) of precipitates extracted from the as-quenched M54® steel.

Fig. 10 .

 10 Fig. 10. Prior austenitic grain size in as-quenched state after 1 h austenitizing at 900 °C (a) and 1100 °C (b).

Fig. 11 .

 11 Fig. 11. Reference JCPDS pattern and experimental XRD profiles (relative intensities) of samples tempered at 516 °C for 10 h and 500 h.

Fig. 12 .

 12 Fig. 12. Three-dimensional APT reconstruction of Ni atoms (green) of a sample tempered at 516 °C for 10 h. Carbides are represented as violet isoconcentration surfaces (total concentration of Mo and C is 36 at. pct). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .

 14 Fig. 14. Proximity histogram of the 98 interfaces precipitate/matrix.

Fig. 15 .

 15 Fig. 15. Chemical composition of M 2 C in M54® calculated with TCF3 ThermoCalc® database.

Table 1

 1 Chemical composition (wt%) of UHS Co-Ni steels.

		C	Cr	Ni	Co	Mo W	V	Ti	Mn Si
	M54®	0.3	1	10	7	2	1.3 0.1	0.02max /	/
	Aermet® 100 0.23 3.1 11.1 13.4 1.2 /	/	0.05max /	/
	AF1410	0.15 2	10	14	1	/	/	0.015	0.1 0.1
	HP9-4-20	0.2	0.8 9	4	1	/	0.08 /	0.2 0.2
	HY-180	0.13 2	10	8	1	/	/	/	0.1 0.05
	S53®	0.21 10	5.5	14	2	1	0.3	0.2max	/	/

Fig. 2. Mole fraction of phase according to austenitizing temperature in M54® calculated with TCFE3 ThermoCalc® database (Ti-free).

Table 2

 2 Composition of M 6 C carbides predicted by ThermoCalc® calculations.

Carbide M 6 C

Composition (860 °C) (Fe 2.8 Mo 2.05 W 0.96 Cr 0.12 V 0.07 )C

Table 3

 3 Austenitization of different UHS steels hardened by M 2 C carbide precipitation.

	Steel	M54®	AerMet® 100	AF1410
	T aust (°C)	1060	885	843

Table 4

 4 Diffusivity in γ-iron and diffusion distance during solutionizing of carbideforming elements.

	Element	Mo	Cr	W
	Diffusivity in γ-iron (D, cm 2 / s)	0.036exp (-239,8/RT)	0.063exp (-252,3/RT)	0.13exp (-267,4/RT)
		[46]	[46]	[46]
	Diffusion distance during	~4	~3	~2
	austenitization (1 h at
	1060 °C) (μm)

Acknowledgements Atom-probe tomography was performed at the Northwestern University Center for Atom-Probe Tomography (NUCAPT). The LEAP tomograph at NUCAPT was purchased and upgraded with grants from the NSF-MRI (DMR-0420532) and ONR-DURIP (N00014-0400798, N00014-0610539, N00014-0910781, N00014-1712870) programs. NUCAPT received support from the MRSEC program (NSF DMR-1121262) at the Materials Research Center, the SHyNE Resource (NSF ECCS-1542205), and the Initiative for Sustainability and Energy (ISEN) at Northwestern University. Special thanks to Dr. Dieter Isheim for his analyses and invaluable help.

Assistance provided by QuesTek Innovations LLC through Chris Kern and Ricardo K. Komai.

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.