A robust framework for task-related resident scheduling

European Journal of Operational Research 276 (2019) 656–675

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Innovative Applications of O.R.

A robust framework for task-related resident scheduling Sebastian Kraul a,∗, Andreas Fügener b, Jens O. Brunner a, Manfred Blobner c a

Chair of Health Care Operations/Health Information Management, Faculty of Business and Economics, University of Augsburg – University Center for Health Sciences at the Augsburg Clinic – UNIKA-T, Neusässer Straße 47, Augsburg 86156, Germany b Department of Supply Chain Management and Management Science, Faculty of Management, Economics and Social Sciences, University of Cologne, Albertus-Magnus-Platz, Köln 50923, Germany c Clincs for Anaesthesiology, Technical University of Munich, Klinikum rechts der Isar, Ismaninger Straße 22, Munich 81675, Germany

a r t i c l e

i n f o

Article history: Received 19 January 2018 Accepted 14 January 2019 Available online 18 January 2019 Keywords: OR in health services Resident scheduling Robustness Column generation Pattern management

a b s t r a c t We consider the training phase of physicians after finishing medical school. They specialize in a common field like ophthalmology or anesthesiology and are called residents. Technological progress in health care leads to increasing complexity in the requirements of physician training. As a consequence, those programs are often not only time-related but also task-related. Task-related means that residents should perform a given number of different interventions in their program. Typically, a resident will follow a rotation across different clinical departments, where the number of performed interventions per period may be estimated. Predicting the exact number of interventions is usually not possible. Accordingly, a resident might not be able to perform all of the required interventions during the planned rotation, resulting in an extension of the program. In this paper, a new model is presented that calculates the number of residents a hospital can reliably train on a strategic level. Our model also provides the corresponding training schedule. It considers minimum requirements of both time-related stays in specific departments as well as task-related interventions that have to be performed. The robustness of the model can be set by management to handle uncertainties in interventions. A Dantzig–Wolfe decomposition is used to accelerate the solution process and a new pattern generation approach that can construct multiple patterns out of one solution is developed. The termination of the column generation algorithm is accelerated significantly by this method. The model is evaluated using real-world data from a resident program for anesthesiology in a German university hospital. The results demonstrate that near-optimal solutions with an average optimality gap of below five percent can be achieved within computation times of few minutes. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Physicians need to specialize in a standard domain like ophthalmology or anesthesia after finishing medical school. In this phase, they are called residents. Even though each specialization has different regulations, the central concept of each program consists of a rotation through different departments. Nevertheless, analyzing programs in various countries will result in two major program concepts. While rotations through different clinical departments are typical in most English-speaking nations like the United States or Great Britain, most European residents have to perform different interventions repeatedly (see Accreditation Council for Graduate Medical Education, 2017; German Medical Association, 2003). Note that the main difference lies in how interventions are handled. The first concept treats interventions indirectly while the sec-

∗

Corresponding author. E-mail address: [email protected] (S. Kraul).

https://doi.org/10.1016/j.ejor.2019.01.034 0377-2217/© 2019 Elsevier B.V. All rights reserved.

ond considers them directly – a minimum number of each type of intervention must be performed to finish the program. In the following, we will call the indirect type time-related program and the direct one task-related program. Please note that in both types of programs an attending physician is supervising the residents and evaluating their training progress. Our study concentrates on the European style of resident programs with a focus on Germany but can be applied to other countries as well. In Germany, there are currently almost 60,0 0 0 residents in one of 57 specialist training programs. Each program is structured by the German Medical Association (GMA), the central organization in the system of medical self-administration in Germany. They set regulations for all specialty trainings by defining the content, duration and objectives of specialty training and specialist designations. The content ranges from general diagnostic skills to specific interventions that need to be performed several times. However, the GMA is only responsible for the medical focus of the resident training program, while labor regulations are defined by law and can be specified individually by the hospital. Indeed, the GMA only

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

formulates two categories of hard restrictions: To become a specialist in a given domain residents need to perform specific interventions a minimum number of times and work in some departments for a minimum number of periods. In the following we will notate the requirement of a minimum number of interventions as task-related and the requirement of working a minimum length of time as time-related. In order to achieve all requirements the GMA defines a recommended training duration depending on the specialty – the program for anesthesiology, for instance, is supposed to take five years. While finishing the training on time is not mandatory from a regulatory point of view, it is strongly demanded by physicians, and many hospitals promise timely training to prospective residents. Defining a minimum number of interventions to be performed ensures the quality of education. However, predicting the exact number of interventions a resident can perform during a stay in a clinical department is typically not possible. This uncertainty can lead to not achieving the minimum number of completed interventions and consequently to an extension of the program. This paper focuses on the training of residents in a department of anesthesiology. In Germany, prior to being considered specialists, anesthesiology residents need to show proficiency on 14 different types of interventions by performing each of them a minimum number of times. To achieve all required interventions, the planner assigns residents to clinical departments for defined time periods, in our example, quarters. Please note that residents might be assigned to multiple departments during the same time period, e. g. they could be assigned to the sport surgery unit with half of their capacity and to the neurosurgery with the remaining half of their capacity. The implementation and assignment of the share is used as an input on the operational level. Switching to a new department between the time periods is time-consuming due to handover and setup times. Contrary to the practice where residents are informed about a change only one or two weeks in advance, with a provided training schedule residents know their sequence of assigned departments from the beginning of their residency (Schlitzkus, Schenarts, & Schenarts, 2010). As a result, both residents and departments can better prepare for handovers and training leading to an increase in quality of training for residents and consequently, a high quality of care for patients (Denson, McCarty, Fang, Uppal, & Evans, 2015; Ellman, Tobin, Stepczynski, & Doolittle, 2016). In Germany, the length of training exceeded the planned time horizon often by more than than 50% resulting in dissatisfaction of and planning uncertainty for residents (German College of General Practitioners & Family Physicians, 2009; Miani et al., 2015). As a consequence of this deficit, German courts judged in 2015 that hospital management needs to provide training schedules guaranteeing the success of resident programs in time (LAG Baden-Württemberg, 2015). Thus, providing long-term training schedules can increase reliability of planning in hospitals. Additionally, complying with these legal requirements can enhance employee satisfaction. Please note that in German medical training programs there is no regulation in which periods the program has to start, i. e., residents can start quarterly. As hospitals face a shortage of residents (European Commission, 2012), reliable, high-quality resident programs might lead to a competitive advantage (Chang, McDonald, & Burton, 2010) and are necessary to attract the best talents. Thus, hospitals planning a training program face the following trade-off: On the one hand, hospitals are interested in admitting many medical residents that offer an inexpensive medical resource. On the other hand, having too many residents can generate a shortage of interventions that can delay the completion of the training program. Resident programs with well-defined scope and systematics are able to address this trade-off (Roos et al., 2011). From the viewpoint of a resident, the quality of a specialist program is strongly influenced by hav-

657

ing a reliable schedule. With this schedule, the resident knows in advance which interventions will be performed in the following periods, in which department the training will take place, and at what time the specialization will be completed. Especially in the context of work-life balance, a fixed duration of the specialization program has a positive effect (Chang et al., 2010). The existence of reliable training schedules offers value to the hospital management as well: the chief physician does not only know in advance who will be joining the department in a specific period but also which departments the (new) resident has visited before. This additional information can lead to more efficiency in the hospital by lowering the loss of productivity incurred during hand-off periods (Denson et al., 2015; Katzenbach & Smith, 1993). The guarantee to finish the training program on time is a competitive advantage of a hospital. Furthermore, hospital management can define an employment strategy that fits into their training program. Using the same schedule for more than one resident increases the systematics of a resident program. While residents can compare their level with fellow residents and share their experience with more junior colleagues, the training supervisor can improve the level of education because of the repetitive structure. The purpose of this paper is to analyze the strategic problem of how many residents of one speciality a hospital can reliably and timely employ and train. Our problem may be classified as a staffing problem according to the framework of Erhard, Schoenfelder, Fügener, and Brunner (2018). As a side result, training schedules for those residents are produced that may be applied on a tactical level. Uncertainty regarding the number of performed procedures is covered using a robust optimization approach (Bertsimas & Sim, 2004) that allows management to choose their level of conservatism in their promise that all residents finish their training on time. The model considers the following five main characteristics: (1) Residents need to perform a minimum number of predefined interventions to finish specialization. (2) Residents need to work a minimum timespan in predefined departments. (3) Residents can perform different and several interventions in each department. (4) The number of interventions performed in each department is uncertain. (5) The number of positions in each department is limited. Our paper is the first approach considering task- and timerelated training requirements, an aspect that is mandatory in physician training in many European countries, and that might increase the quality of training in health systems where task-related requirements are not mandatory. This constitutes our main contribution. Our novel compact formulation of the problem can be seen as a task allocation and sequencing problem with upper and lower bound temporal constraints. Bertsimas and Weismantel (2005) proved that this type of problem is NP-hard and, consequently, not tractable. As commercial solvers are not able to solve the compact formulation of our model in realistic time (in our case no feasible positive solution was found within five hours), we further present a sophisticated column generation heuristic that detects near-optimal solutions (on avg. ≤ 5% optimality gap) within minutes. In particular, a new pattern generation approach for cyclic problems is developed that further decreases solution times significantly. We apply our model in a real-world case study of a German university hospital with more than 1200 beds. We analyze the importance of robustness and the performance of our solution approach. In the case study, minimum training requirements are set for both the number of predefined interventions and the time spent in specific departments. Thus, this paper demonstrates how hospitals could apply the proposed model to create a competitive

658

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

advantage in the recruiting market by providing reliable training schedules. The remainder of this paper is organized as follows. Section 2 reviews the relevant literature on physician scheduling with an emphasis on resident programs. In Section 3, a mathematical formulation of the problem is described. A Dantzig–Wolfe reformulation, as well as the solution method, are presented in Section 4. The method is numerically evaluated and compared to a standard approach in Section 5. The paper finishes with a conclusion in Section 6, where future research avenues are sketched. 2. Literature review The literature of personnel scheduling problems has been studied in great detail over the last decades. Therefore, a vast amount of published work exists in this area focusing on shift, break, and tour scheduling (Burke, de Causmaecker, Berghe, & van Landeghem, 2004; Ernst, Jiang, Krishnamoorthy, & Sier, 2004; Van den Bergh, Beliën, de Bruecker, Demeulemeester, & de Boeck, 2013). Nevertheless, most of the work discusses the operational level, e. g., by deciding the shift sequence of the individual worker. This review focuses on resident and physician scheduling. A comprehensive review of physician scheduling can be found in Erhard et al. (2018). Franz and Miller (1993) are among the first to formulate the resident scheduling problem. They rotate residents through several departments and maximize the preferences of the residents by considering the American accreditation system. A mixed integer program (MIP) is formulated, and the model is solved using a rounding heuristic. Ozkarahan (1994) develops a goal programming model for a shift scheduling problem of residents over a seven-day period following the American accreditation system. He includes hospital regulations, residency program requirements, personal preferences as well as weekend assignments. The influence of night shifts while considering departmental staffing and skill requirements is discussed by Sherali, Ramahi, and Saifee (2002). They focus on the American accreditation system and describe the problem as a MIP. It is solved using a heuristic solution procedure. A fixed number of night shifts is used as a training condition. Beliën and Demeulemeester (2006) develop schedules for trainees at a hospital department without addressing a specific accreditation system. However, the case-study is based on a Belgian hospital. They minimize the number of irregular assignments, assuming a trainee can perform only one activity in a period and has to repeat the activity for a given time before the activity can change. They formulate the problem as a binary problem and solve it using a branch-and-price algorithm after a Dantzig–Wolfe reformulation. Topaloglu (2006) describes a multi-objective programming model for scheduling emergency medicine residents following the Turkish accreditation system and solves it using a goal programming approach. The time horizon is one month while focusing on different skill levels and respecting day and night shifts. Interventions are not considered in this paper. Cohn, Root, Kymissis, Esses, and Westmoreland (2009) discuss the same problem as Topaloglu but apply a time horizon of a complete year and solve it using a heuristic approach. In contrast to Topaloglu they focus on the American accreditation system. A long-term staff scheduling of physicians with different experience levels is described by Brunner and Edenharter (2011) as a MIP and solved using a column generation heuristic. They minimize the cost of hiring physicians over a time horizon of one year and construct flexible shifts without addressing a specific accreditation system. Bard, Shu, and Leykum (2013) construct monthly schedules for interns and residents in the American accreditation system by using a goal programming formulation. They optimize a model with 12 different soft constraints considered in the objective function, such as clinic utilization. A network-

based approach for assigning duties to internals and residents in primary care clinics is analyzed by Bard, Shu, and Leykum (2014). They show that team considerations and on-services taken into account can be formulated as a network. They evaluate the effectiveness of their model using a comparative study for the 2012–13 academic year of the Internal Medicine Department at the University of Texas Health Science Center in San Antonio. Guo, Morrison, Jacobson, and Jokela (2014) prove that the problem of a 1-year resident rotation schedule focusing on the American accreditation system is NP-complete, and develop a polynomial-time algorithm to solve a special case of their problem. Smalley and Keskinocak (2014) formulate a medical rotation and shift scheduling program for the American accreditation system. They focus on fairness by considering preferences and equality of experience. The equality of experience is measured by the number of assignments to a department. Bard, Shu, Morrice, Leykum, and Poursani (2016c) analyze annual block scheduling for family medicine residency programs in the American accreditation system. They generate a twelve-month schedule by designing blocks and measuring the effect of duties during rotation by adjusting the schedule in a second step. The problem is formulated as a MIP, and several optimization-based heuristics are developed to yield good feasible solutions. Similar problems are handled with two different solution strategies as well (Bard, Shu, Morrice, & Leykum, 2016a; 2016b). Proano and Agarwal (2017) present a multi-objective optimization approach by developing balanced weekly rotation schedules for a 1-year time horizon. They facilitate continuity of care by using clinic rotation schemes. They solve the goal programming model with a standard solver and generate different solutions with the use of Analytical Hierarchy Process. Brech, Ernst, and Kolisch (2018) analyze the training tardiness for a given stock of surgery residents in a German hospital over the complete training program. They apply a resource-constrained project scheduling problem solved using a matheuristic combining Benders decomposition and ant colony optimization. To summarize the literature on resident scheduling, it can be stated that most of the articles are focusing on the operational or tactical level by targeting rotations across departments without considering the number of performed interventions. Only two papers, namely Beliën and Demeulemeester (2006) and Brech et al. (2018), deal with procedures at all. The approach of Beliën and Demeulemeester (2006) only allows one activity to be performed at a time, and deals with irregular assignments. Performing specific procedures, however, is not a mandatory part in the training regulations they discuss. While Brech et al. (2018) do include some kind of task-related requirements, they assume deterministic occurrence of procedures and ignore time-related requirements. From the literature review, it can be concluded that to the best of our knowledge there is no work concerning the strategic point of view of scheduling residents across their whole training program in a task- and time-related framework by calculating the optimal workforce size for residents. Moreover, currently no work exists which observes the impact of disruptions by task uncertainty in a long-term schedule. The paper at hand closes this gap. 3. Problem description and model development Suppose there is a given set of residents I = {1, . . . , |I|} and a set of time-periods T = {1, . . . , |T |} where |T| is equivalent to the length of the training program, e. g., 5 years, and each period is three months long. Additionally, there is a set of subperiods U representing the weeks of each period t ∈ T. The subperiods are an option to allow more than one assignment in a period t. Every resident has an individual starting period ri ∈ T of the program. Please note that possible starting periods can be limited if there are any

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

program or hospital restrictions, e. g., if new residents can only start once a year. The function

 k ( ri , t ) =

t − ri + 1,

t ≥ ri

(1)

|T | + t − ri + 1, t < ri

is used to determine the training-period of resident i in period t of his individual schedule. The two cases represent the cyclic structure and ensure that k(ri , t) ∈ T. For instance, a resident starting the training program in period t = ri = 2 is in the first training period k ( 2, 2 ) = 1. The hospital has a set of departments J = {1, . . . , |J|} with a time-dependent number of available resident positions djt . To fulfill the time-related program requirements, every resident has to work at least mj full-time equivalent periods in each department j. In addition, resident i cannot be assigned to more than sk(ri ,t ) departments in training period k(ri , t). To fulfill the task-related requirements, different interventions described by the set  ∈ L = {1, . . . , |L|} must be performed n times each. While it is straightforward to assume that the time a resident spends in a department can easily be predicted, the number of interventions the resident may perform during one period contains uncertainty. Thus, if less interventions than planned are performed, finishing the residency on time may be at risk. As management promises timely training, a robust approach according to Bertsimas and Sim (2004) is applied. In this approach, an interval of uncertainty is assumed for each possible type of intervention containing an upper bound aj and lower bound a j − aˆ j that describes how often intervention  can be performed in department j per full-time period. For instance, a resident assigned to neurosurgery only for a single period can perform transfusions between 46 and 38 times. The lower bound can be seen as a worst-case scenario since residents need to perform a minimum number of interventions to finish their program. As it is not realistic that the worst-case will be realized in all periods and departments, a parameter   is used to describe in how many periods and departments the worst-case might happen (see Bertsimas & Sim, 2004). This parameter is used to control the conservatism of the management. Note that   can be chosen within  the interval [0, Ni ] with Ni = min{ t∈T sk(ri ,t ) , |T | · |J |}, where J is the set of clinical departments where interventions  can be performed. We further introduce the set Sˆ = {( j, t )| j ∈ J, t ∈ T } of combinations of department j and period t to simplify notation. The binary decision variables x˜i jk(ri ,t )u define whether a resident i is working in department j in subperiod u of period t or not. The real variables xi jk(ri ,t ) are used to calculate the share of working 

u∈U x˜i jk (r ,t )u

i time in department j and period t, such that xi jk(ri ,t ) = |U | for all I, J, T. The general assignment of residents to departments in a specific period is set by the binary variable

 yi jk(ri ,t ) =

1, 0,

if resident i is assigned to department j in training period k(ri , t ) otherwise.

Note that x and y are dependent on each other since a resident cannot work in a department without being assigned to it. Additionally, the binary variables zi denote whether a resident can finish the specialist program on time or not:

 zi =

1, 0,

if resident i finishes the training program on time otherwise.

Note that the z variable also determines whether the resident is part of the workforce, as our objective only considers residents that are able to be trained within the time horizon. Please find a summary of all used notation in the following.

659

Sets with indices I set of residents (index i) J set of departments (index j) L set of interventions (index ) T set of periods in the planning horizon (index t) U set of subperiods (index u) Sˆ set of tuples of departments and periods Parameters d jt maximum number of full-time positions in department j in period t ri starting period of resident i sk(ri ,t ) the maximum number of departments that resident i can be assigned to in period t mj minimum number of full-time equivalent working periods in department j amount of interventions  that can be performed in a j department j in any period aˆ j maximum deviation of the amount of interventions  in department j in any period n amount of interventions  that need to be performed in the resident program  number of periods and departments with the worst-case number of interventions  Decision variables xi jk(ri ,t ) share of subperiods that resident i works in department j and period t out of 100% x˜i jk(ri ,t )u 1 if resident i is assigned to department j in subperiod u of period t, 0 otherwise yi jk(ri ,t ) 1 if resident i is assigned to department j in period t, 0 otherwise zi 1 if resident i can finish the training program within the time-horizon, 0 otherwise

max



zi

(2a)

i∈I

s.t. 

xi jk(ri ,t ) ≤ d jt

∀ j ∈ J, t ∈ T

(2b)

yi jk(ri ,t ) ≤ sk(ri ,t )

∀i ∈ I, t ∈ T

(2c)

xi jk(ri ,t ) ≤ zi

∀i ∈ I, t ∈ T

(2d)

xi jk(ri ,t ) ≥ m j zi

∀i ∈ I, j ∈ J 

(2e)

i∈I

 j∈J

 j∈J

 t∈T



a j xi jk(ri ,t ) −

j∈J t∈T



max {Si ∪{( j∗ ,t ∗ )}|Si ⊆Sˆ, |Si |=  ,( j∗ ,t ∗ )∈Sˆ\Si }

aˆ j xi jk(ri ,t )

( j,t )∈Si



+ ( −  )aˆ j∗  xi j∗ k(ri ,t ∗ ) ≥ n zi xi jk(ri ,t ) − yi jk(ri ,t ) ≤ 0  x˜i jk(ri ,t )u xi jk(ri ,t ) = u∈U |U | x˜i jk(ri ,t )u − x˜i jk(ri ,t )u+1 ≥ 0

∀i ∈ I, l ∈ L

(2f)

∀i ∈ I, j ∈ J, t ∈ T

(2g)

∀i ∈ I, j ∈ J, t ∈ T

(2h)

∀i ∈ I, j ∈ J, t ∈ T , u ∈ {1, . . . , |U | − 1} (2i)

0 ≤ xi jk(ri ,t ) ≤ 1

∀i ∈ I, j ∈ J, t ∈ T

(2j)

x˜i jk(ri ,t )u , yi jk(ri ,t ) , zi ∈ {0, 1}

∀i ∈ I, j ∈ J, t ∈ T , u ∈ U

(2k)

The objective (2a) maximizes the number of residents that can be trained on time. Constraints (2b) restrict the total workload of all residents per department j and period t. These constraints

660

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

are different to classic personnel scheduling approaches where departments have a demand that must be fulfilled (Ernst et al., 2004). Limiting the total workload of assigned residents for each department stems from the fact that the number of performed interventions is limited. Otherwise, it would be possible for a department to generate patients. Even if a hospital can determine the number of interventions on its own, there are capacity restrictions, e. g., time, so that the number of interventions is limited. The total number of departments that resident i can be assigned to within period t is restricted in constraints (2c) (see KC, 2014). Please note that the maximum number of parallel assignments may depend on the period of training k(ri , t). New residents starting their training might be assigned to more departments in order to obtain a broad overview, or to less departments in order to keep focus, compared to more senior residents. To control the maximum workload of each resident i in period t, constraints (2d) are set up. At this planning-level only full-time contracts are considered with a workload of 100%. To fulfill the time-related requirements, resident i should work at least a predefined number of full-time equivalent periods in each department j over the time-horizon – this is ensured with constraints (2e). Constraints (2f) guarantee that resident i fulfills all interventions n times required for the resident program, i. e., the task-related requirements. The first term on the left-hand side is the sum of all possible performed interventions  aj that can be achieved with the assigned departments over the entire time horizon by resident i, i. e., the upper bound of the uncertainty interval. The second term considers the uncertainty of the realized number of interventions  in the assigned departments. The worst-case occurs in a total of   occurrences, i. e., the lower bound of the uncertainty interval. Note that constraints (2f) are non-linear. In the linear version constraints (2f) are replaced by (3a)–(3d).



 a j xi jk(ri ,t ) −



j∈J t∈T

 ρi jt +  qi ≥ n zi

∀i ∈ I,  ∈ L

j∈J t∈T

(3a)

ρi jt + qi ≥ aˆ j xi jk(ri ,t )

∀i ∈ I, j ∈ J,  ∈ L, t ∈ T

(3b)

ρi jt ≥ 0

∀i ∈ I, j ∈ J,  ∈ L, t ∈ T

(3c)

qi ≥ 0

∀i ∈ I,  ∈ L

(3d)

The step-wise linearization follows the logic of Bertsimas and Sim (2004) and is discussed in detail in Appendix A. Please note that constraints (2d)–(2f) connect the assignment of resident workload to departments with finishing the training program on time. Constraints (2g) connect the x and y variables. A resident i is assigned to a department j in period t as soon as any positive amount of workload of this resident is assigned to the respective combination of department and period. Constraints (2h) connect the x variables with the x˜ variables, while constraints (2i) define the structure of the x variable space. Herewith, symmetry can be reduced without loss of generality by assuming that variable values are monotonically decreasing in subperiods. This property can be used because the order of subperiods is not relevant for our problem setting as the schedule only describes the share of working time per period. The binary x˜ variables prevent assignments of marginal capacity by ensuring that each assignment to a department corresponds to one or more weeks full time-equivalent per period but nothing in between. Please note that this leaves flexibility in operational planning since the exact order of assignments is not fixed within our model. Variables are defined in constraints (2j) and (2k). A solution of the model (2) provides the maximum number of residents that can be reliably trained within the time horizon

as well as an employment strategy in the sense that the hospital knows how many residents should be hired in which period. In addition to that, a schedule for all residents is available describing the relative capacity each resident is assigned to each department in each period. From a managerial perspective, the training supervisor has the opportunity to vary the setup by manipulating the parameters sk(ri ,t ) (maximum number of simultaneous departments in a training period), mj (minimum workload requirements per department during the time horizon), and   (level of robustness, i. e., conservatism of the planner), while the minimum amount of interventions n is typically regulated. Note, for some departments and resident programs, the minimum working time mj is regulated as well. Be aware that set I needs to be sufficiently large for all starting periods to calculate the maximum feasible number of residents. The model might be extended by using seniority levels for different stages of the resident program. However, different seniority levels are not needed in our case of residents in anesthesiology. 4. Decomposition of the task-related resident scheduling problem The compact formulation of the robust task-related resident scheduling problem (2) is not able to find non-trivial feasible solutions in realistic time by using a standard solver, i. e., using Gurobi 7.0 no feasible solution except for 0 is detected for our case study after a runtime of five hours (see Section 5 for the detailed setting). Thus, an alternative solution strategy is developed in this section. As the block structure of the problem suggests the use of a decomposition using a Dantzig–Wolfe reformulation (Barnhart, Johnson, Nemhauser, Martin W. P. Savelsbergh, & Vance, 1998; Dantzig & Wolfe, 1960; Desrosiers & Lübbecke, 2005), we decompose the compact robust task-related resident scheduling problem (2) using a Dantzig–Wolfe reformulation and develop a column generation heuristic to solve it. We decompose our model by resident. Additionally, the residents can be aggregated as they are assumed to be homogeneous in our model which eliminates symmetry. The compact formulation is decomposed into a master problem (MP) and one or more subproblems (SP). Since in most cases not all columns of a MP are known in advance, a restricted MP is formulated which consists of a subset of columns. In the following, the restricted MP will be written as MP. In the iterative process, SPs are applied to generate new columns, while the dual solution of the MP controls the generation process. Decomposition of the robust task-related resident scheduling problem. For the formulation of the MP, a new set of columns P = {1, . . . , |P|} must be introduced (representing extreme points). Every p ∈ P specifies a feasible training schedule ensuring the resident is able to finish the training on time. Please note that set P is increasing during the solution process as new training schedules are generated in the SP. Note that in our reformulation extreme rays are not necessary since the SP is bounded (Barnhart et al., 1998). The resulting MP (extensive formulation) of the robust task-related resident scheduling problem can be formulated as follows: Sets with indices P set of training schedules (index p) R set of starting periods (index r) Parameters X p jk(r,t ) share of subperiods of a generic resident with training schedule p in department j and period t Decision variables number of residents that start in period r and work according to training schedule p

λ pr

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

max



λ pr

(4a)

ρ jt ≥ 0

661

∀ j ∈ J,  ∈ L, t ∈ T

(6j)

p∈P r∈R

s.t. 

q ≥ 0 X p jk(r,t ) λ pr ≤ d jt

∀ j ∈ J, t ∈ T



π jtD



(4b)

0 ≤ x jk(r,t ) ≤ 1

(4c)

x˜ jk(r,t )u , y jk(r,t ) ∈ {0, 1}

p∈P r∈R

λ pr ∈ Z+

∀ p ∈ P, r ∈ R

The objective (4a) maximizes the sum of contracted residents which is equal to objective (2a) of the compact formulation. Just as in constraint (2b), the maximum number of full-time equivalent residents in department j and period t is considered in constraint (4b). The parameter Xpjk(r,t) is analogous to the compact formulation decision variable xi jk(ri ,t ) . Note that function k is no longer related to a specific resident i but to a generic training schedule with starting period r in the MP formulation. The values of parameter X are determined in the SP as detailed in the following paragraph. Since the dual formulation of the MP does not provide relevant insights, only the dual variables π jtD ≥ 0 of constraints (4b) are detailed here. They are used to determine the reduced cost of a generic MP column and are considered in the objective of the SP. The generic reduced cost of a column can be written as in (5).

1−



X p jk(r,t ) π jtD

(5)

j∈J t∈T

Because of the aggregation of residents in the reformulation, only one SP needs to be generated for every starting period r. The solution of the SP(r) can be interpreted as a training schedule for a generic resident starting in start period r. The parameters and variables in the SP will be the same as in the compact formulation (2) of Section 3 including the linearization of constraints (2f) while dropping index i. With this in common, the SP(r) can be formulated as:

SP (r ) : min



Djt x jk(r,t )

(6a)

j∈J t∈T

s.t. 

y jk(r,t ) ≤ sk(r,t )

∀t ∈ T

(6b)

x jk(r,t ) ≤ 1

∀t ∈ T

(6c)

x jk(r,t ) ≥ m j

∀j ∈ J

(6d)

∀ ∈ L

(6e)

∀ j ∈ J,  ∈ L, t ∈ T

(6f)

∀ j ∈ J, t ∈ T

(6g)

∀ j ∈ J, t ∈ T

(6h)

j∈J

 j∈J

 t∈T



 a jt x jk(r,t ) −



j∈J t∈T

ρ jt +  q ≥ n

j∈J t∈T

ρ jt + q ≥ aˆ j x jk(r,t ) x jk(r,t ) − y jk(r,t ) ≤ 0  x jk(r,t ) =



u∈U

x˜ jk(r,t )u

|U |

x˜ jk(r,t )u − x˜ jk(r,t )u+1 ≥ 0

∀ j ∈ J, t ∈ T , u ∈ {1, . . . , |U | − 1}

(6i)

∀ ∈ L

(6k)

∀ j ∈ J, t ∈ T

(6l)

∀ j ∈ J, t ∈ T , u ∈ U

(6m)

Objective (6a) of the SP(r) results from the generic reduced cost (5) of the dual MP. Since the π variables are already calculated before the SP(r) is solved, they are handled as parameters. Constraints (6b)–(6m) relate to Constraints (2c)–(3d) in Section 3. In the sense of column generation, any new feasible training schedule is integrated into the MP if and only if the objective value is smaller than 1, i. e., a positive reduced cost column can be generated. In this case, the values of xjk(r,t) for all j ∈ J, t ∈ T define a new column, i. e., training schedule, p with starting period r:



1 X p jk(r,t )



,

which is added to the set P in MP. Pattern management. The cyclic structure of the problem inhibits symmetry that decreases convergence of column generation even though the projection in the pricing problem already handled part of it (Desrosiers & Lübbecke, 2011). As an example, let us assume that sk(r,t) is constant for all r. In each iteration, SP(r) yields the same solution for all r ∈ T as illustrated in the left column of Fig. 1, where all residents work in the same department in each period, and each subproblem provides the sequence A → B → C. Therefore, one SP is sufficient. Even by reducing the number of solved SPs to one, this solution scheme is inefficient since more feasible columns, i. e., training schedules, are formulated indirectly by one solution. The new cyclic pattern generation process presented in this paper uses this information and creates several training schedules from one SP solution – even if parameter sk(r,t) varies for different r. Keeping the example with the department sequence of A → B → C as an initial SP solution, we additionally create the feasible schedules with the order C → A → B and B → C → A. Thus, additional schedules are generated out of one solution by shifting each schedule through all starting periods as illustrated in the right column of Fig. 1. In this context, it is important to note that the additionally created schedules do not necessarily have positive reduced costs. Nevertheless, the schedules are feasible and can be added as new columns. Remark: every column generated in the procedure is not part of the MP in the previous stage. The experimental study in Section 5 shows that creating additional schedules is more efficient in terminating the column generation process compared to the classical pattern generation by a factor between 10 and 50. The solution process is visualized in Fig. 2. After the initialization of MP and SPs, the linearized MP is solved. Each time the MP is solved, we update the parameters in SP(r) using the dual variable values. The SP(r) is then solved for all starting periods. If a column with positive reduced cost is found, the cyclic pattern generator constructs exactly one column for every starting period. This process is repeated until no new column with positive reduced costs is found. After that, the best found solution from the linearized MP provides a valid upper bound for the MP, and the integer MP is solved using the generated columns. The feasible solution defines the lower bound of our model. There is no need to embed the procedure into a search tree to prove optimality since the experimental study shows that near-optimal solutions are achieved.

662

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Fig. 1. Example of identical generations (left), example of shift generation (right).

Fig. 2. Illustrative flow chart of the column generation heuristic.

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

663

Table 1 Interventions  for specialist training in anesthesiology. Intervention () Repetitions (n ) Intervention () Repetitions (n )

Reanimation 10 Head 100

CVC 50 Paediatric 50

Transfusion 50 Medulla 100

General 1800 Fiberoptic 50

Abdominal 300 Regional 25

Obstetric 50 Intrathoracic 25

Caesarian 25 Intracranial 25

Fig. 3. Demand profile djt of department j for a generic year.

5. Case study of anesthesia residents in a German hospital In this section, we apply our model in a real-world case study concerning anesthesiology residents in a German university hospital and evaluate its performance. For this purpose, Section 5.1 discusses the case study. By testing different parameter settings, an evaluation of the proposed column generation heuristic is part of Section 5.2. In Section 5.3, the effect of robustness is evaluated. The study considers the training of anesthesiology residents with a total training duration of five years. The period length is set to three months resulting in a total of 20 periods, i. e. |T | = 20. The German Medical Association recommends 14 specific interventions, i. e., |L| = 14, with a varying number of repetitions to qualify for the final board exam as a specialist in anesthesiology. In Table 1, the minimum number of required interventions n are listed, e. g., a resident needs to perform at least 10 reanimations to become an anesthesiologist. All computations are performed on a 2.7 gigahertz (Intel® CoreTM i7 − 3740QM CPU) with 16 gigabytes RAM running under the Windows 8.1 Enterprise operating system. The column generation algorithm is coded in Python 3.5, and Gurobi 7.0 is used to solve the LP-relaxation of the MP, the SP as well as the MP as IP after column generation terminates. The default settings of Gurobi are used. To accelerate the column generation process, the first feasible solution with positive reduced cost < 1 (see (6a)) is selected as a new column for the MP. As a result, the SP is not always solved to optimality. However, the scheme guarantees optimality. In the final step, the MP is solved as an IP with a time restriction of 150 seconds. Since there is no time-related order in the interventions, only one SP needs to be solved in each iteration. Nevertheless, the pattern generation method presented in Section 4 is applied in the computational study. Note that other programs might require a specific order. 5.1. Real-world case For the computations, real-world data supplied by a large university hospital in Germany with more than 1200 beds is applied.

The number of concerned departments j is given by |J| = 13. The maximum amount of full-time positions djt in each department j and period t is provided by the respective chief physician. This decision was made based on the number of senior physicians required for supervising residents, and the number of patients in each department. The resulting profiles for one generic year are presented in Fig. 3. These profiles are applied for the entire time horizon of five years. An interesting observation is the variation of the maximum number of full-time equivalent residents per period by at most one resident in a department over the year, e. g., of four residents in the first quarter and of five residents in the remaining quarters in the department of plastic surgery. The upper bound of interventions aj that are available in one department j in any period, as well as the worst-case term aˆ j , are estimated based on historical data of over 26 months and expert interviews with hospital management. An overview summarizing these values for all departments and interventions is provided in Tables 8 and 9 in Appendix B. Based on discussions with hospital management, we further set sk(r,t ) = 3 for all t ∈ T and m j = 1/3 for all j ∈ J. Consequently, a resident can be assigned to at most three different departments in one period and has to work during the complete length of training at least one month in each of the 13 departments. Results of the realworld case are presented in Table 2. We evaluate four scenarios – a deterministic scenario (DET) with  = 0 ∀ ∈ L, an optimistic robust scenario (5%) with  = 0.05N ∀ ∈ L, a conservative robust scenario (50%) with  = 0.5N ∀ ∈ L, and a worst-case scenario (SOY) where  = N ∀ ∈ L – the labeling of scenario SOY is referring to Soyster (1973). We use the scenario DET and SOY since they represent the best and the worst-case for the hospital. They can be interpreted as the limits of the uncertainty interval. Note, the last row of Table 2 provides information about the average values over all scenarios. While the first column of Table 2 describes the scenario, the next two columns provide the upper and lower bound of our model solution. As a first result, it can be stated that the range between the robust worst-case (UB SOY) and the deterministic case (UB DET) amounts to 78 residents. Ignoring this gap can cause several problems – the most obvious one is that probably no resident will be able to finish the training on time if the determin-

664

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675 Table 2 Computational results of the real-world case. Scenario

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

DET 5% 50% SOY Avg.

88 83 31 10 –

79 79 30 10 –

10.23 4.82 3.23 0.00 4.57

287.0 395.1 262.9 1.6 236.7

680 1320 220 80 575.0

34 66 11 4 28.8

istic schedule is confronted with a stochastic reality. Even though  -values are 50% of the worst-case scenario, the upper bound is close to the worst-case with a value of 31. One reason is that a resident is often assigned to only one department in a period. In other words, a intervention  can vary only once in a period t from the upper value aj . As a consequence, the worst-case can already be achieved with a  -value N . Note that the effect of robustness will be reviewed in more detail in Section 5.3. The solution time for all four cases is on average 236.7 seconds with an optimality gap of 4.57%, both values are provided in columns 4 and 5. For the SOY scenario, the heuristic detects the optimal solution. The total number of generated columns, as well as the number of iterations, are provided in the last two columns. The heuristic constructs on average 575.0 training schedules in 28.8 iterations. Note that a more extensive computational study with more instances is carried out in the subsequent section. As expected, the resulting number of residents that can be trained on time varies depending on the conservatism of the management. By training ten residents only, the hospital looses the opportunity to employ more residents as cheap resources, while most residents will not be able to finish their training on time if 88 residents are employed. Both cases can result in a monetary loss. The former by hiring more expensive specialists and the latter by extending the duration of training, and followed by losing reputation on the job market. Moreover, using pre-defined schedules generated by the model will increase patient service quality because chief physicians of different departments know the skill level of each resident in advance. In case of the conservative robust scenario (50%), the training plan suggests to employ a total of 30 residents during a 5-year cycle. One or two residents start each alternating period, and a total of six training schedules is applied, i. e., on average five residents work with the same training schedule. In the real-world case, 84 anesthesia residents are currently employed at the hospital. This number is close to the upper bound of the 5% case. To evaluate how good the conservatism of the planner fits into the actual situation, the data is analyzed by counting how often the occurrences of interventions were below the worstcase. Based on the given 26 months the worst-case was realized much more often than in 5% of cases, e. g., in 50% of the cases for abdominal interventions. Thus, the hospital is currently not able to educate all residents on time. As a conclusion, it can be stated that the presented method can determine the maximum number of residents without violating the length of training for a real-world case. In addition, the proposed model creates a managerial tool. The hospital management can adopt a training setup by changing the maximum number of departments sk(r,t) that can be visited in a period as well as changing the minimum time mj of working in a specific department. Additionally, conservatism of planning can be adopted by changing  -values. 5.2. Evaluation of the algorithm We present two numerical studies to evaluate the performance of our solution algorithm. In the first one, the effect of different

problem sizes is analyzed, while in the second study the new pattern generation is compared to the classic column generation procedure. The generated instances are based on the data of the realworld case in Section 5.1. To analyze the algorithm for different problem sizes, instances will be clustered into three groups – low, mid, and high. Since the problem size is mainly driven by the upper bound of interventions aj and the worst-case term aˆ j , new upper bounds were generated by drawing a new value from the distribution of interventions N (a j , aˆ j ) assuming a normal distribution with the worst-case as a proxy for the standard deviation. To create different problem sizes, the distribution is rescaled by scaling the mean value as well as the standard deviation with a constant C for each of the three groups. The other moments of the distribution will stay the same. Instance group low is based on the distribution of the real-world case in Section 5.1 resulting in C = 1. For group mid the scaling factor is C = 2 while for instance group high the distribution is manipulated by C = 3. Except for this variation, all other parameters remain unchanged. The  -values are set to 25% of the worst-case. Every group low, mid, and high consists of 20 instances, resulting in a total of 60 instances. The results of the first experimental study are presented in Tables 3–5. Column 1 refers to the numeration of test instances of each group whereas column 2 and 3 give information on the upper and lower bound of the problem after column generation terminates. The GAP of each solution is stated in column 5 while the average GAP over all instances is 5.09% for group low, 4.71% for group mid, and 3.83% for group high. This is a remarkable result since the compact formulation could not detect any non-zero feasible solution for any instance within a five-hour runtime using Gurobi in standard configuration. The solution time in seconds is displayed in column 5. As a first but not surprising result, it can be stated that solving the SPs is most time-consuming in the procedure, as SP is solved as a MIP in every iteration. In contrast to that, the relaxed MP can be solved quickly. The average total time of solving an instance takes about 1032.0 seconds for group low, 823.6 for group mid, and 897.6 for the instances of group high. Since GAP and Total Time are two main characteristics to compare the algorithm performance of different group sizes, we take a look at the confidence intervals of both items in Fig. 4. We use a t-statistic with a 95% interval for single testing. We are aware that not using multiple testing might increase the alpha error. We can see that the results are not significantly different from one another. All intervals overlap each other in some parts so that we can assume that the algorithm is not influenced significantly by moderately changing the problem size. The next two columns illustrate the total number of generated columns as well as the total number of iterations. On average, the heuristic constructs 443.0 columns in each instance for group low and 501.0 or 828.0 columns for instances of the group mid or high. The last column presents the average number of newly generated columns with a (solution) value greater than 0 in the next solution of the LP-relaxation of the MP. This value, which can be between 1 and 20, provides information about the efficiency of the new pattern generation process. The higher the value, the more of the newly generated columns are used, i. e., enter the basis in the next MP main iteration.

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Table 3 Computational results with parameter variation in aj for group low. Instance

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Avg.

60 28 53 22 27 68 52 30 53 84 50 39 54 37 45 44 61 34 44 25 −

58 27 50 21 26 63 50 28 50 80 48 38 51 35 42 42 57 31 42 24 −

3.33 3.57 5.66 4.55 3.70 7.35 3.85 6.67 5.66 4.76 4.00 2.56 5.56 5.41 6.67 4.55 6.56 8.82 4.55 4.00 5.09

479.9 99.7 183.5 252.4 578.4 3406.5 456.7 103.7 854.2 1, 393 519.7 1061.6 2145.7 523.4 334.1 255.9 3376.2 71.6 3619 925 1032.0

540 200 220 300 560 1300 260 100 320 600 320 360 700 200 160 200 1300 80 920 220 443.0

27 10 11 15 28 65 13 5 16 30 16 18 35 10 8 10 65 4 46 11 22.2

17.7 14.0 18.1 17.7 19.3 19.8 18.2 20.0 19.1 20.0 17.9 20.0 20.0 20.0 20.0 18.6 20.0 20.0 19.5 19.3 19.0

Table 4 Computational results with parameter variation in aj for group mid. Instance

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Avg.

100 91 63 48 92 102 50 60 132 84 85 98 67 88 83 71 82 103 92 83 –

92 87 61 46 88 98 47 57 126 80 82 92 64 85 79 67 78 99 88 79 –

8.00 4.40 3.17 4.17 4.35 3.92 6.00 5.00 4.55 4.76 3.53 6.12 4.48 3.41 4.82 5.63 4.88 3.88 4.35 4.82 4.71

341.3 277 107.4 122.3 648.6 497.6 1557.9 497 856 1316.9 939.9 1286.2 932.8 514.2 545 2760.2 292.9 1698.2 682.5 598.2 823.6

620 380 140 200 480 420 520 260 860 560 560 1, 060 460 320 460 860 260 760 480 360 501.0

31 19 7 10 24 21 26 13 43 28 28 53 23 16 23 43 13 38 24 18 25.1

19.4 19.4 18.0 17.9 19.2 19.4 20.0 20.0 20.0 20.0 20.0 19.8 20.0 20.0 20.0 19.2 20.0 20.0 19.7 19.4 19.6

Table 5 Computational results with parameter variation in aj for group high. Instance

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Avg.

123 112 73 81 131 115 147 189 151 176 145 144 155 70 171 99 176 167 111 136 –

117 103 69 78 128 112 143 186 146 171 137 138 151 67 166 90 173 164 108 130 –

4.88 8.04 5.48 3.70 2.29 2.61 2.72 1.59 3.31 2.84 5.52 4.17 2.58 4.29 2.92 9.09 1.70 1.80 2.70 4.41 3.83

785 172.2 295 656.5 452.4 523.8 1203.7 1293.8 452.3 1053.6 345.6 401.4 739.7 212.2 2814.1 564.9 2286.8 1429.3 912.5 1358.7 897.7

1280 440 300 940 780 580 660 1500 580 860 600 320 740 440 1500 400 1280 1280 720 1360 828.0

64 22 15 47 39 29 33 75 29 43 30 16 37 22 75 20 64 64 36 68 41.4

19.7 19.1 18.6 19.6 19.5 20.0 18.9 20.0 19.6 19.4 19.9 20.0 19.4 20.0 19.6 19.5 19.6 18.6 19.6 20.0 19.5

665

666

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Fig. 4. Confidence interval of the GAP (left) and Total Time (right) for group low, mid, and high.

Table 6 Comparison of the cyclic pattern generation approach and the standard method for group low. Cyclic pattern generation

Avg.

Standard pattern generation

Total Time (seconds)

Columns

Iterations

Total Time (seconds)

Columns/Iterations

1032.0

443.0

22.2

6314.2

394.3

To understand the importance of the last column of Tables 3–5, let us take a look back at the theory of column generation. In the classic column generation approach and in analogy to a simplex method, at each iteration of solving the SP, a basic/ non-basic variable switch is performed. In other words, exactly one variable solution that was greater than 0 will become 0 while another variable (the one generated in the last iteration of the SP) will become greater than 0. In the cyclic pattern generation process used in this paper, several variables are generated out of one solution (see Section 4). Since only one of the variables is related to the actual solution of the SP, this one will be the new basic variable. The other variables might also become greater than 0, but this is not verified by the SP – meaning that some variables might not have positive reduced costs. Nevertheless, the experimental study shows that on average over all instances and iterations, 19.3 of the newly generated variables are used in the next solution of the relaxed MP at each iteration, demonstrating the effectiveness of our pattern generation process. To further investigate the value of our pattern generation process, the second experimental study compares the cyclic pattern generation approach against the standard method of column generation. For this purpose, the 20 instances of group low (see Table 3) are analyzed, and the results are presented in Table 6. Since the upper bound is the same for both approaches, only the duration and number of iterations, as well as the number of columns, are compared. Please note that the number of iterations equals the number of columns in the standard pattern generation process. The main result is that the cyclic pattern generation approach outperforms the standard method in terms of computation time by a factor between 1.5 and 33.3. On average, the standard method needs 6.1 times as long as our pattern generation approach. One reason for this difference is the reduction of solving computationally expensive SPs. As it can be seen in Table 6, the number of iterations that is needed to terminate column generation is on average roughly 17.8 times higher in the standard method than in the new approach. However, the cyclic pattern generation approach requires nearly the same number of generated columns as the standard method for termination. To underline this result, Fig. 5 shows how fast the method converges against the final upper bound. Since both figures present the convergence over all instances, the parameters are normalized with 1 being the starting solution and 0 being the final UB. In the left part of Fig. 5, the convergence with the cyclic pattern generation approach is shown. It can be seen

that the method is on average after seven iterations (140 columns) less than 10% below the optimal UB. In contrast to that, in the right part of Fig. 5, the convergence of the standard method is presented. Even after 150 iterations the process is more than 10% away from the optimal bound. As a result, it can be stated that the cyclic pattern generation approach converges about 20 times faster than the classical approach. 5.3. Effect of robustness Even though a strategic problem is solved, uncertainty has a significant influence on the solution. From a managerial perspective, the consequences of overestimating the number of residents that can be trained on time can result in an increase of monetary as well as of non-monetary costs. In this part of the experimental study, the effect of different  -values is analyzed by comparing the ranges of solutions. For this purpose, the instances of the first study in Section 5.2 (see Tables 3–5) are extended for a wider range of different robustness parameters  . For each of the 60 instances, we test values for  representing 5%, 10%, 15%, 25%, 50%, and 75% of the worst-case for each intervention. Note that these settings are applied to ensure comparability in the computational study – in a practical usage each  -value should be adapted independently and discussed with the management. The confidence intervals of the different settings are presented in Table 7. A detailed overview of all instances is provided in Appendix C. Columns 1 and 2 show the group and the applied  -value. The intervals of the upper and lower bound are given in column 3 and 4. Between each group the difference in the upper as well as the lower bound is different. Additionally, both values decrease the higher the  -value is. However, within the chosen step size, the values are not significantly different. In order to visualize the effect of  , we use the total relative change of the upper bound by fixing a base case of 5% as the benchmark. For instance, the first scenario of group low has an upper bound of 81 for the 5% case and an upper bound of 76 for the 10% case (see Table 10 in Appendix C). The relative change for −76 this instance is 1 − 8181 = 0.94. The effect of the relative change in the upper bound by varying  -values is illustrated in Fig. 6. The figure illustrates the average relative change of the upper bound over all instances of one group compared to 5%. We can see that the relation between the objective and the  -value seems to be non-linear. For small  -values the effect seems nearly the same

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

667

Fig. 5. Convergence cyclic pattern generation (left), convergence standard generation (right).

Table 7 Confidence intervals for the computational results with  -variation for all groups. Group

Gamma

UB

LB

low mid high low mid high low mid high low mid high low mid high low mid high

0.05 0.05 0.05 0.1 0.1 0.1 0.15 0.15 0.15 0.25 0.25 0.25 0.5 0.5 0.5 0.75 0.75 0.75

[76.4, 59.8] [134.9, 113.3] [225.6, 185.0] [71.2, 54.8] [124.8, 104.3] [202.6, 167.4] [63.9, 49.4] [112.2, 92.4] [183.8, 151.7] [52.9, 38.1] [93.0, 74.4] [150.1, 117.2] [36.9, 23.7] [60.8, 43.4] [115.4, 83.5] [31.5, 19.4] [47.5, 31.9] [92.4, 60.3]

[71.4, 56.2] [129.8, 108.0] [221.5, 181.2] [67.0, 51.5] [120.1, 99.4] [198.3, 163.5] [59.8, 46.6] [107.4, 88.1] [179.8, 147.6] [50.2, 36.1] [88.6, 70.9] [145.4, 112.3] [35.2, 23.0] [58.3, 41.6] [111.2, 80.0] [29.9, 18.8] [45.3, 30.5] [88.9, 57.6]

Fig. 6. Relative change of the upper bound by varying  values for group low, mid, and high.

for all three problem sizes, while the effect on the objective seems to flatten out when further increasing  . As a result, it can be stated that the effect of increasing the level of robustness does not differ by problem size. Nevertheless, the effect of changing the level of robustness is quite high over all problem sizes with differences of up to 90%. Ignoring such uncertainty can result in wrong training schedules, staff sizes, and as a

consequence, employee dissatisfaction. From a managerial perspective, it is useful to evaluate several  settings and find a suitable solution for the hospital based on data-driven analyses.

6. Conclusion Integrating task-related requirements in resident training programs increases complexity of planning. Obviously, more requirements make planning more difficult per se, but more importantly, task-related requirements add uncertainty to the problem. Contrary to time-related requirements, task-related requirements put the resident in danger of not being able to finish the training on time. In this paper, we propose a robust model that considers this additional property and maximizes the number of residents that the hospital can train within the planning horizon. Management can decide on the level of robustness, i.e., how sure they want to be that all residents finish on time. The model proposed in this paper is formulated in a compact version. As standard solvers are not able to solve it in realistic time, we further developed a Dantzig–Wolfe reformulation. The model is formulated in a generic way so that an adaption to specific training programs is easy to handle. As demonstrated in the experimental study, our solution method can find good solutions quickly – not only for theoretical instances but also for a real-world case of anesthesia residents in a German university hospital with more than 1,200 beds. The heuristic generates nearoptimal solutions with an average GAP of less than 5%. With the newly developed cyclic pattern generation approach, the solution process can be accelerated significantly, as the proposed column generation heuristic can terminate on average 6.1 times faster than the standard approach. Additionally, we demonstrated the effects of different levels of robustness. The staff size can vary by up to 90% depending on the conservatism of the responsible manager. Beside these technical insights, the model provides a managerial toolkit. The training supervisor can calculate an optimal employment strategy and can analyze the effects of different parameter settings. Moreover, residents receive a training schedule that communicates the assignments of their residency in advance. Thus, our approach helps to increase the quality of specialization programs as well as the satisfaction of supervisors and residents. A remaining challenge is to solve the model to optimality – a different formulation of the subproblem embedded in a branch-and-price framework could be a valuable research avenue. Additionally, the model might be extended by including seniority of the residents. While seniority does not necessarily have to be considered on the strategic level of planning – in our setting every resident can be assigned to all departments from the start of the program – other programs may need such a regulation.

668

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Even though task-related training programs are relevant in the practice area, the topic is not prevalent in the literature. Other issues that can be addressed in future work might be the operational perspective of scheduling residents, where continuity of care and fairness might be considered in more detail, and the consideration of different uncertainties as well as different options of handling those.

constraint (i, ) (Bertsimas & Sim, 2004). A problem in this formulation is that the robustness is non-linear. However, for a given vector x∗ the value of β i (x∗ ,  ) is equal to the optimal objective function of the following linear problem (LP) given in (9).

βi (x,  ) = max

aˆ j xi jk(ri ,t ) ζi jt

(9a)

j∈J t∈T

s.t. 

Acknowledgment



ζi jt ≤ 

(9b)

j∈J t∈T

The authors would like to thank the editor and the three anonymous referees for their valuable comments which helped to improve the quality of this paper. This research project is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 405488489.

Appendix A. Linearization of robustness The robustness concept applied in this paper is originally formulated by Bertsimas and Sim (2004). The idea is motivated as follows: Let a brief version of the task-related resident scheduling problem be formulated as

max



0 ≤ ζi jt ≤ 1

s.t.

ρi jt + qi ≥ aˆ j xi jk(ri ,t )

xi jk(ri ,t ) ≤ d jt

(7)

∀ j ∈ J, t ∈ T

a j xi jk(ri ,t ) ≥ n zi

∀i ∈ I,  ∈ L



zi

(11a)

i∈I

and consider constraint (i, ) of problem (7). Furthermore, let Sˆ be the set of coefficients a j , ( j, t ) ∈ Sˆ with parameter uncertainty aˆ j , i. e., the uncertainty interval. For every (i, ) there is a parameter   describing the total number of times that the worstcase of the uncertainty interval occurs in the planning horizon, i. e., a j − aˆ j . As presented in Bertsimas and Sim (2004),   can be seen as a lever for the robustness of the model. The higher the value of   , the more conservative the solution of the model becomes. Note that   can be chosen within the interval [0, Ni ]  with Ni = min{ t∈T sk(ri ,t ) , |T | · |J |}, where J is the set of clinical departments where interventions  can be performed. If all   are set to the highest possible value, the solution is robust by the definition of Soyster (1973) – considering a worst-case, i. e., in each period only a j − aˆ j interventions  are realized in department j. If all   are zero, the solution equals the deterministic case, i. e., in each period aj interventions  in department j are available. Taking this into account, the robust version of (7) can be formulated as

s.t. 

∀ j ∈ J, t ∈ T

xi jk(ri ,t ) ≤ d jt

(11b)

i∈I



 a j xi jk(ri ,t ) −



j∈J t∈T

 ρi jt +  qi ≥ n zi

∀i ∈ I,  ∈ L

j∈J t∈T

(11c)

ρi jt + qi ≥ aˆ j xi jk(ri ,t )

∀i ∈ I, j ∈ J,  ∈ L, t ∈ T

(11d)

ρi jt ≥ 0

∀i ∈ I, j ∈ J,  ∈ L, t ∈ T

(11e)

zi

i∈I



∀ j ∈ J, t ∈ T

By respecting x as a decision variable, formulation (9) is nonlinear. In contrast to that, the dual formulation (10) is linear. This result can be used to reformulate (8) in a linear way:

max

j∈J t∈T



(10)

∀ j ∈ J, t ∈ T

qi ≥ 0

ρi jt ≥ 0



qi ≥ 0 xi jk(ri ,t ) ≤ d jt

(8)

∀ j ∈ J, t ∈ T

i∈I



a j xi jk(ri ,t ) − βi (x∗ ,  ) ≥ n zi

∀i ∈ I,  ∈ L

j∈J t∈T



βi (x∗ ,  ) = max{S ∪{( j∗ ,t ∗ )}|S ⊆Sˆ,|S |=  ,( j∗ ,t ∗ )∈Sˆ\S }{ ( j,t )∈Si i i i i aˆ j xi jk(ri ,t ) + ( −  )aˆ j∗  xi j∗ k(ri ,t ∗ ) } is the protection function of and

ρi jt +  qi

j∈J t∈T

zi

i∈I

s.t.



min

s.t. 

(9c)

By ordering the ζ variables by their objective coefficient aˆ j xi jk(ri ,t ) in a descending order and assuming   is integer, the first   variables are equal to 1, while the remaining are equal to 0. In addition to that, the dual formulation of (9) has the same optimal objective value since it is feasible and bounded for all   :

i∈I

max

∀ j ∈ J, t ∈ T

∀i ∈ I,  ∈ L

(11f)

The maximization term β i   ) in the constraints (8) is replaced by the objective of the dual model (10). Furthermore, the new variables need to be restricted, and therefore, the constraints of (10) are added as well. Note that qi and ρ ijt are the corresponding dual variables of constraints (9b) and (9c), respectively. To linearize the model (2) presented in this section, constraints (2f) have to be replaced by constraints (11c)–(11f). (x∗ ,

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Appendix B. Data of the real-world case

Table 8 Average value aj of interventions  in each department j per quarter. j\

Reanimation

CVC

Transfusion

General

Abdominal

Obstetric

Caesarian

Ophthalmic Vascular surgery Gynecology Hand surgery ENT surgery MJF surgery Neurosurgery Orthopedy Plastic surgery Sport surgery Traumatology Urology Visceral surgery

0.0 1.7 0.0 0.0 0.5 0.5 0.9 1.1 0.0 0.0 0.4 0.6 1.8

0.0 51.7 9.8 2.3 4.0 36.7 127.8 26.0 12.7 3.2 16.9 34.8 138.7

0.4 47.5 10.8 0.4 2.9 31.4 46.4 33.8 12.0 0.7 22.7 17.9 51.4

208.8 275.2 359.0 91.6 321.6 368.6 369.0 248.9 203.2 232.0 342.9 435.3 528.6

0.2 12.4 19.9 0.0 0.0 0.4 1.1 0.2 0.5 0.0 1.4 61.9 127.1

0.0 0.0 100.1 0.0 0.0 0.0 0.0 2.5 2.5 0.0 0.0 2.5 0.0

0.0 0.0 89.8 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.5 0.0

j\

Head

Paediatric

Medulla

Fiberoptic

Regional

Intrathoracic

Intracranial

Ophthalmic Vascular surgery Gynecology Hand surgery ENT surgery MJF surgery Neurosurgery Orthopedy Plastic surgery Sport surgery Traumatology Urology Visceral surgery

191.0 11.1 2.5 2.4 303.9 258.2 136.2 2.4 4.2 2.4 2.5 2.8 8.1

14.1 0.0 0.0 0.0 26.7 19.4 0.4 0.3 1.0 0.2 0.0 3.8 0.3

0.0 4.9 73.9 3.1 0.0 2.9 3.2 14.5 3.5 5.7 11.1 80.9 100.8

0.0 11.7 0.2 47.7 0.0 0.2 0.0 18.1 41.1 10.2 22.1 1.3 3.6

0.0 0.2 0.2 0.0 1.8 17.4 1.1 0.2 0.4 0.1 0.1 0.8 1.6

1.0 0.9 0.0 0.0 0.0 16.8 116.5 0.0 1.3 0.0 0.0 0.0 0.0

0.0 2.4 0.0 0.0 0.0 1.5 7.1 0.3 0.2 0.0 0.8 0.0 47.2

Table 9 Worst-case term aˆ j of interventions  in each department j per quarter. j\

Reanimation

CVC

Transfusion

General

Abdominal

Obstetric

Caesarian

Ophthalmic Vascular surgery Gynecology Hand surgery ENT surgery MJF surgery Neurosurgery Orthopedy Plastic surgery Sport surgery Traumatology Urology Visceral surgery

0.0 1.3 0.0 0.0 0.3 0.3 0.7 0.9 0.0 0.0 0.3 0.4 1.1

0.0 7.6 2.3 0.6 1.2 10.2 16.7 4.0 5.0 1.3 5.1 4.9 13.4

0.3 8.9 3.3 0.3 1.9 7.6 8.4 7.2 6.7 0.6 7.5 4.6 8.9

40.2 30.1 32.6 55.8 42.4 40.1 39.2 23.0 32.0 34.3 56.4 37.1 44.8

0.1 9.3 12.0 0.0 0.0 0.3 0.9 0.1 0.4 0.0 1.1 37.2 75.1

0.0 0.0 10.6 0.0 0.0 0.0 0.0 0.1 0.2 0.0 0.0 0.1 0.0

0.0 0.0 10.6 0.0 0.0 0.0 0.0 0.0 0.2 0.0 0.0 0.1 0.0

j\

Head

Paediatric

Medulla

Fiberoptic

Regional

Intrathoracic

Intracranial

Ophthalmic Vascular surgery Gynecology Hand surgery ENT surgery MJF surgery Neurosurgery Orthopedy Plastic surgery Sport surgery Traumatology Urology Visceral surgery

40.2 7.2 0.2 0.6 42.4 35.0 34.1 0.2 2.9 0.2 0.5 0.7 4.1

6.3 0.0 0.0 0.0 8.0 6.1 0.3 0.2 0.8 0.2 0.0 1.7 0.2

0.0 1.2 9.3 2.0 0.0 0.4 0.6 3.6 1.3 1.5 3.3 10.0 9.2

0.0 4.3 0.2 37.8 0.0 0.2 0.0 5.4 10.7 4.9 9.0 1.0 1.5

0.0 0.1 0.2 0.0 1.5 4.6 0.9 0.2 0.3 0.1 0.1 0.7 1.2

0.8 0.3 0.0 0.0 0.0 6.4 22.0 0.0 1.1 0.0 0.0 0.0 0.0

0.0 1.9 0.0 0.0 0.0 1.2 2.9 0.3 0.2 0.0 0.6 0.0 6.1

669

670

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Appendix C. Results of the -variation in Section 5.3

Table 10 Computational results with parameter variation in  for group low. Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5

0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.25 0.25 0.25 0.25 0.25

81 44 87 41 46 86 71 55 74 106 68 64 77 54 78 66 88 72 64 39 76 40 81 33 42 83 63 47 72 100 65 57 71 51 73 63 80 67 60 36 66 36 60 32 38 69 59 44 68 97 60 52 64 48 62 59 72 60 53 34 60 28 53 22 27

79 42 83 39 43 80 66 54 70 98 62 61 72 51 72 62 80 66 60 36 74 38 76 31 38 77 59 43 68 94 61 55 67 50 69 59 74 62 56 34 66 34 56 31 35 63 56 43 64 90 57 50 59 46 57 53 67 55 50 32 58 27 50 21 26

2.47 4.55 4.60 4.88 6.52 6.98 7.04 1.82 5.41 7.55 8.82 4.69 6.49 5.56 7.69 6.06 9.09 8.33 6.25 7.69 2.63 5.00 6.17 6.06 9.52 7.23 6.35 8.51 5.56 6.00 6.15 3.51 5.63 1.96 5.48 6.35 7.50 7.46 6.67 5.56 0.00 5.56 6.67 3.13 7.89 8.70 5.08 2.27 5.88 7.22 5.00 3.85 7.81 4.17 8.06 10.17 6.94 8.33 5.66 5.88 3.33 3.57 5.66 4.55 3.70

419.5 473.4 233.7 1011.5 513.1 397.1 259.7 148.8 408.4 427.1 436.9 317 965.8 323.7 648.8 362.6 682.8 600.9 744 1107.2 552.9 282.2 384.1 180.7 1017.3 1483.9 445.7 193 1293.3 1229 1067 862.5 2365.5 536 970.9 778.9 1879.5 1703.9 2665.3 858.7 138.2 336.7 139.2 483.8 727.6 396.1 468 299.4 3660.6 4224.2 1287 712.1 1382.5 1504.7 356.4 1137.7 2354.8 1264 1520 2547.9 479.9 99.7 183.5 252.4 578.4

780 680 580 1180 660 400 360 240 320 520 340 280 620 280 520 340 540 420 380 500 760 460 640 220 960 820 240 160 700 600 600 420 940 260 600 400 840 720 1020 420 280 440 160 600 800 300 260 160 1420 1400 460 320 800 620 240 500 920 500 540 840 540 200 220 300 560

39 34 29 59 33 20 18 12 16 26 17 14 31 14 26 17 27 21 19 25 38 23 32 11 48 41 12 8 35 30 30 21 47 13 30 20 42 36 51 21 14 22 8 30 40 15 13 8 71 70 23 16 40 31 12 25 46 25 27 42 27 10 11 15 28

19.2 19.4 19.5 19.7 19.5 20 20 20 20 20 20 20 20 20 20 20 20 19.7 20 20 19.5 19.1 19.4 18.4 19.7 20 20 20 20 20 19.5 19 19.7 19.1 20 20 18.7 20 20 18.4 19.3 19.1 17.4 19.3 19.5 20 20 20 19.6 20 20 20 20 19.8 19.6 19.7 19.5 20 19.6 19.7 17.7 14 18.1 17.7 19.3

(continued on next page)

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Table 10 (continued) Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75

68 52 30 53 84 50 39 54 37 45 44 61 34 44 25 25 12 35 0 10 47 34 19 36 64 38 18 36 18 36 31 46 30 27 13 22 0 20 0 0 20 30 17 28 57 25 15 28 0 26 23 30 0 15 0

63 50 28 50 80 48 38 51 35 42 42 57 31 42 24 24 11 34 0 10 44 34 18 34 60 37 18 35 18 34 30 43 29 27 13 20 0 20 0 0 20 29 17 26 53 24 15 27 0 25 22 29 0 14 0

7.35 3.85 6.67 5.66 4.76 4.00 2.56 5.56 5.41 6.67 4.55 6.56 8.82 4.55 4.00 4.00 8.33 2.86 – 0.00 6.38 0.00 5.26 5.56 6.25 2.63 0.00 2.78 0.00 5.56 3.23 6.52 3.33 0.00 0.00 9.09 – 0.00 – – 0.00 3.33 0.00 7.14 7.02 4.00 0.00 3.57 – 3.85 4.35 3.33 – 6.67 –

3406.5 456.7 103.7 854.2 1393 519.7 1061.6 2145.7 523.4 334.1 255.9 3376.2 71.6 3619 925 35.7 90.7 109.9 2.2 89.8 2820 758.8 210.4 972.8 2703.6 6444.5 742.9 1532.7 1367 337.1 577.4 4099.5 4520.7 768.3 3851.4 67.9 2.7 36.6 2.2 2.4 1055.8 780.6 66.3 1441.5 6341.2 792.1 2393.4 5458.8 2.3 113.8 423.9 1501.4 2.1 455.6 2.4

1300 260 100 320 600 320 360 700 200 160 200 1300 80 920 220 40 120 100 0 80 440 240 140 240 640 1200 180 560 160 100 220 620 360 280 200 60 0 40 0 0 80 260 80 280 1340 200 160 700 0 80 120 200 0 120 0

65 13 5 16 30 16 18 35 10 8 10 65 4 46 11 2 6 5 0 4 22 12 7 12 32 60 9 28 8 5 11 31 18 14 10 3 0 2 0 0 4 13 4 14 67 10 8 35 0 4 6 10 0 6 0

19.8 18.2 20 19.1 20 17.9 20 20 20 20 18.6 20 20 19.5 19.3 20 12.8 15.3 – 16.7 20 18.5 10.4 20 20 19.4 15.8 19.7 20 18 15.7 20 19 19.1 20 17.5 – 20 – – 20 20 20 20 19.6 18.6 20 19.5 – 20 20 18.7 – 20 –

Table 11 Computational results with parameter variation in  for group mid. Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

1 2 3 4 5 6 7 8 9 10 11 12

0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05

164 158 134 90 137 137 84 97 172 116 107 135

159 152 132 86 133 132 78 92 167 108 103 132

3.05 3.80 1.49 4.44 2.92 3.65 7.14 5.15 2.91 6.90 3.74 2.22

976.7 506.3 593.9 602.8 762.1 158.2 470 381.8 330.5 545.5 280 219.5

1760 1400 1280 1060 1480 500 420 560 780 920 500 600

88 70 64 53 74 25 21 28 39 46 25 30

19.8 18.7 19.7 19.6 19.8 20 20 20 20 19.4 20 19.6

(continued on next page)

671

672

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

Table 11 (continued) Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6

0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.50 0.50 0.50 0.50 0.50 0.50

109 120 121 110 122 120 132 117 149 144 121 82 126 127 73 87 162 109 100 124 102 114 110 99 108 118 123 113 136 97 110 72 118 117 64 76 147 104 95 115 76 106 100 79 100 116 112 106 100 91 63 48 92 102 50 60 132 84 85 98 67 88 83 71 82 103 92 83 56 56 36 29 42 67

104 112 117 103 118 113 127 109 145 137 119 77 121 123 67 80 156 104 97 121 96 108 108 92 105 110 120 109 131 91 107 69 112 111 60 71 140 99 93 111 72 102 97 74 97 111 107 100 92 87 61 46 88 98 47 57 126 80 82 92 64 85 79 67 78 99 88 79 53 53 34 28 40 65

4.59 6.67 3.31 6.36 3.28 5.83 3.79 6.84 2.68 4.86 1.65 6.10 3.97 3.15 8.22 8.05 3.70 4.59 3.00 2.42 5.88 5.26 1.82 7.07 2.78 6.78 2.44 3.54 3.68 6.19 2.73 4.17 5.08 5.13 6.25 6.58 4.76 4.81 2.11 3.48 5.26 3.77 3.00 6.33 3.00 4.31 4.46 5.66 8.00 4.40 3.17 4.17 4.35 3.92 6.00 5.00 4.55 4.76 3.53 6.12 4.48 3.41 4.82 5.63 4.88 3.88 4.35 4.82 5.36 5.36 5.56 3.45 4.76 2.99

433.7 121 168.6 413.9 129.7 110.3 338.2 239.4 1268.5 776.3 277.6 644.7 832 405.9 526.7 421.3 692.6 644.2 772.4 277.2 1341.8 666.6 183.7 1095.8 256.1 537.1 609.6 644.2 1300.4 41.5 731.2 552 1128.4 695.2 945.1 302.4 529.2 4337.3 571.6 863.8 450.5 829.6 306.6 456.5 486.3 1754.5 951.5 1548.7 341.3 277 107.4 122.3 648.6 497.6 1557.9 497 856 1316.9 939.9 1286.2 932.8 514.2 545 2760.2 292.9 1698.2 682.5 598.2 36.1 169 93.4 72.8 159.1 255.4

560 320 580 640 400 280 600 440 1740 1480 820 920 1300 620 500 520 880 640 840 620 1080 700 380 500 360 680 700 580 1500 120 960 660 1440 600 460 340 680 1380 580 1060 320 540 480 340 540 960 720 760 620 380 140 200 480 420 520 260 860 560 560 1060 460 320 460 860 260 760 480 360 40 160 80 60 80 120

28 16 29 32 20 14 30 22 87 74 41 46 65 31 25 26 44 32 42 31 54 35 19 25 18 34 35 29 75 6 48 33 72 30 23 17 34 69 29 53 16 27 24 17 27 48 36 38 31 19 7 10 24 21 26 13 43 28 28 53 23 16 23 43 13 38 24 18 2 8 4 3 4 6

20 20 18.8 19.5 17.6 20 19.5 20 19.8 19.9 19.2 19.6 19.9 20 20 20 20 19.8 20 18 20 19.4 19.3 19.5 20 19.6 20 19.6 19.7 16.6 19.7 19.7 19.9 18.6 20 20 19.7 20 20 20 19.2 20 20 19.2 19.1 20 20 19.1 19.4 19.4 18 17.9 19.2 19.4 20 20 20 20 20 19.8 20 20 20 19.2 20 20 19.7 19.4 20 18.6 17.3 14.5 19.3 20

(continued on next page)

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675 Table 11 (continued) Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75

18 32 102 60 64 64 42 56 56 35 55 72 42 58 45 48 35 25 21 64 17 0 83 46 52 36 29 25 47 23 35 52 36 35

18 30 98 58 63 63 40 53 54 34 51 69 40 55 42 45 34 23 20 62 17 0 79 44 50 35 28 25 46 21 34 48 34 33

0.00 6.25 3.92 3.33 1.56 1.56 4.76 5.36 3.57 2.86 7.27 4.17 4.76 5.17 6.67 6.25 2.86 8.00 4.76 3.13 0.00 4.82 4.35 3.85 2.78 3.45 0.00 2.13 8.70 2.86 7.69 5.56 5.71

726.5 1741 1050.4 2370.2 789.2 1250.9 673.4 2623.8 557.8 1445.5 1108.9 2495.7 253.9 1099.4 35.7 111.4 97.4 64.9 38 706.8 2509.2 2.3 688.9 3237.9 685.1 174.6 315.7 1037.5 515.2 72.9 356 4917.3 335.5 2609.9

200 440 500 520 600 460 240 660 320 340 420 860 100 540 40 140 80 60 40 280 100 0 320 940 400 100 160 180 280 80 140 1220 160 720

10 22 25 26 30 23 12 33 16 17 21 43 5 27 2 7 4 3 2 14 5 0 16 47 20 5 8 9 14 4 7 61 8 36

18.5 20 20 19.5 19.1 20 20 19.6 20 18.2 20 19.3 20 20 20 18 17.7 20 20 20 20 20 20 20 20 18.8 20 20 20 20 19.8 20 20

Table 12 Computational results with parameter variation in  for group high. Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10

169 211 159 124 176 166 213 270 251 235 223 237 220 127 248 178 257 255 180 206 155 190 147 111 160 148 197 231 223 220 197 211 205 114 223

165 204 155 122 174 162 207 266 248 231 220 232 216 124 244 174 253 252 177 200 151 184 142 107 157 144 191 226 218 216 194 207 200 112 221

2.37 3.32 2.52 1.61 1.14 2.41 2.82 1.48 1.20 1.70 1.35 2.11 1.82 2.36 1.61 2.25 1.56 1.18 1.67 2.91 2.58 3.16 3.40 3.60 1.88 2.70 3.05 2.16 2.24 1.82 1.52 1.90 2.44 1.75 0.90

748.1 901.8 574.9 677.4 364.8 218.5 173.6 466.4 290.9 317.0 202.6 314.7 242.3 214.7 189.4 353.9 403.5 342.7 362.6 217.1 429.6 757.1 525.6 248.3 240.4 399.9 390.0 173.4 458.1 812.3 297.4 247.3 1095.5 394.7 476.0

1560 2600 1160 1580 1340 560 660 1320 980 880 760 1020 920 700 800 720 1320 1240 1020 620 1440 1840 900 760 920 620 860 660 10 0 0 1340 780 780 1600 740 1400

78 130 58 79 67 28 33 66 49 44 38 51 46 35 40 36 66 62 51 31 72 92 45 38 46 31 43 33 50 67 39 39 80 37 70

19.3 19.9 19.9 19.9 19.8 20 20 20 19.4 20 20 20 19.6 20 20 20 20 19.8 20 20 19.7 19.8 19.6 19.5 19.8 19.1 19.8 19.8 19.3 19.8 20 20 19.8 19.7 20

(continued on next page)

673

674

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675 Table 12 (continued) Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 6 7 8 9 10

0.10 0.10 0.10 0.10 0.10 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75

168 231 223 163 183 145 170 135 99 151 136 190 214 198 199 179 189 179 99 203 151 209 200 145 164 123 112 73 81 131 115 147 189 151 176 145 144 155 70 171 99 176 167 111 136 82 75 53 41 88 84 124 154 126 130 102 137 88 30 130 97 136 132 83 97 43 68 36 18 45 59 89 132 104 107

165 227 217 160 179 139 165 131 95 148 131 184 211 194 194 177 185 173 96 201 146 206 196 141 161 117 103 69 78 128 112 143 186 146 171 137 138 151 67 166 90 173 164 108 130 77 70 51 40 85 81 120 148 118 124 93 134 83 29 127 93 134 130 81 94 40 63 35 17 41 57 85 127 98 105

1.79 1.73 2.69 1.84 2.19 4.14 2.94 2.96 4.04 1.99 3.68 3.16 1.40 2.02 2.51 1.12 2.12 3.35 3.03 0.99 3.31 1.44 2.00 2.76 1.83 4.88 8.04 5.48 3.70 2.29 2.61 2.72 1.59 3.31 2.84 5.52 4.17 2.58 4.29 2.92 9.09 1.70 1.80 2.70 4.41 6.10 6.67 3.77 2.44 3.41 3.57 3.23 3.90 6.35 4.62 8.82 2.19 5.68 3.33 2.31 4.12 1.47 1.52 2.41 3.09 6.98 7.35 2.78 5.56 8.89 3.39 4.49 3.79 5.77 1.87

746.9 956.5 654.7 725.4 431.6 1506.9 812.0 898.2 362.8 931.7 269.1 2706.9 692.1 975.9 479.9 635.2 635.2 519.4 474.8 816.2 907.6 2066.6 589.1 941.3 462.7 785.0 172.2 295.0 656.5 452.4 523.8 1203.7 1293.8 452.3 1053.6 345.6 401.4 739.7 212.2 2814.1 564.9 2286.8 1429.3 912.5 1358.7 10 0 0.9 97.9 114.0 211.6 513.7 594.2 3944.8 1868.6 1904.6 920.2 457.6 3782.7 257.3 108.2 5616.7 1388.6 2946.5 3185.4 2603.3 5814.2 247.1 385.2 36.2 128.9 41.8 1026.9 1086.4 692.1 1180.2 4823.6

860 1680 1260 1020 920 2420 1480 1120 820 1340 420 1520 1220 1160 700 1040 1120 960 1060 1320 820 20 0 0 1080 1020 960 1280 440 300 940 780 580 660 1500 580 860 600 320 740 440 1500 400 1280 1280 720 1360 10 0 0 140 100 180 460 440 1180 1280 880 600 480 1180 240 120 1620 560 860 1080 920 1820 280 240 40 40 60 340 460 540 620 860

43 84 63 51 46 121 74 56 41 67 21 76 61 58 35 52 56 48 53 66 41 100 54 51 48 64 22 15 47 39 29 33 75 29 43 30 16 37 22 75 20 64 64 36 68 50 7 5 9 23 22 59 64 44 30 24 59 12 6 81 28 43 54 46 91 14 12 2 2 3 17 23 27 31 43

20 20 19.9 20 20 19.8 19.7 19.8 19.7 19.7 20 20 20 19.4 20 19.8 20 19.7 19.8 20 19.2 19.5 20 20 19.8 19.7 19.1 18.6 19.6 19.5 20 18.9 20 19.6 19.4 19.9 20 19.4 20 19.6 19.5 19.6 18.6 19.6 20 19.6 17.8 19 18.2 19.5 20 20 20 20 19.2 20 19.7 20 20 19.9 20 19.3 19.8 20 19.6 18.8 18.6 20 17 16 20 20 19 19.1 20

(continued on next page)

S. Kraul, A. Fügener and J.O. Brunner et al. / European Journal of Operational Research 276 (2019) 656–675

675

Table 12 (continued) Instance

Gamma

UB

LB

GAP (%)

Total Time (seconds)

Columns

Iterations

Avg. col.

11 12 13 14 15 16 17 18 19 20

0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75 0.75

80 104 59 16 121 73 117 115 72 69

77 102 55 16 119 69 112 111 70 66

3.75 1.92 6.78 0.00 1.65 5.48 4.27 3.48 2.78 4.35

301.0 2359.3 308.8 14.7 3087.9 738.8 777.0 3842.3 3403.6 5054.6

240 880 160 80 1040 300 420 840 940 1080

12 44 8 4 52 15 21 42 47 54

20 19.8 18.3 15 19.2 20 18.7 19.7 19.7 20

References Accreditation Council for Graduate Medical Education (2017). Common program requirements. http://www.acgme.org/. Bard, J. F., Shu, Z., & Leykum, L. (2013). Monthly clinic assignments for internal medicine housestaff. IIE Transactions on Healthcare Systems Engineering, 3(4), 207–239. doi:10.1080/19488300.2013.857370. Bard, J. F., Shu, Z., & Leykum, L. (2014). A network-based approach for monthly scheduling of residents in primary care clinics. Operations Research for Health Care, 3(4), 200–214. doi:10.1016/j.orhc.2014.08.002. Bard, J. F., Shu, Z., Morrice, D. J., & Leykum, L. K. (2016a). Annual block scheduling for internal medicine residents with 4+1 templates. Journal of the Operational Research Society, 67(7), 911–927. doi:10.1057/jors.2015.109. Bard, J. F., Shu, Z., Morrice, D. J., & Leykum, L. K. (2016b). Constructing block schedules for internal medicine residents. IISE Transactions on Healthcare Systems Engineering, 164(2), 1–14. doi:10.1080/19488300.2016.1255284. Bard, J. F., Shu, Z., Morrice, D. J., Leykum, L. K., & Poursani, R. (2016c). Annual block scheduling for family medicine residency programs with continuity clinic considerations. IIE Transactions, 48(9), 797–811. doi:10.1080/0740817X. 2015.1133942. Barnhart, C., Johnson, E. L., Nemhauser, G. L., Martin, W. P. S., & Vance, P. H. (1998). Branch-and-price: Column generation for solving huge integer programs. Operations Research, 46(3), 316–329. Beliën, J., & Demeulemeester, E. (2006). Scheduling trainees at a hospital department using a branch-and-price approach. European Journal of Operational Research, 175(1), 258–278. doi:10.1016/j.ejor.2005.04.028. Bertsimas, D., & Sim, M. (2004). The price of robustness. Operations Research, 52(1), 35–53. doi:10.1287/opre.1030.0065. Bertsimas, D., & Weismantel, R. (2005). Optimization over integers. Belmont, Mass.: Dynamic Ideas. Brech, C.-H., Ernst, A., & Kolisch, R. (2018). Scheduling medical residents’ training at University Hospitals. European Journal of Operational Research. doi:10.1016/j.ejor. 2018.04.003. Brunner, J. O., & Edenharter, G. M. (2011). Long term staff scheduling of physicians with different experience levels in hospitals using column generation. Health Care Management Science, 14(2), 189–202. doi:10.1007/s10729- 011- 9155- x. Burke, E. K., de Causmaecker, P., Berghe, G. V., & van Landeghem, H. (2004). The state of the art of nurse rostering. Journal of Scheduling, 7(6), 441–499. doi:10. 1023/B:JOSH.0 0 0 0 046076.75950.0b. Chang, A., McDonald, P., & Burton, P. (2010). Methodological choices in work-life balance research 1987 to 2006: A critical review. The International Journal of Human Resource Management, 21(13), 2381–2413. doi:10.1080/09585192.2010.516592. Cohn, A., Root, S., Kymissis, C., Esses, J., & Westmoreland, N. (2009). Scheduling medical residents at Boston University School of Medicine. Interfaces, 39(3), 186–195. doi:10.1287/inte.1080.0369. Dantzig, G. B., & Wolfe, P. (1960). Decomposition principle for linear programs. Operations Research, 8(1), 101–111. doi:10.1287/opre.8.1.101. Denson, J. L., McCarty, M., Fang, Y., Uppal, A., & Evans, L. (2015). Increased mortality rates during resident Handoff periods and the effect of ACGME duty hour regulations. The American Journal of Medicine, 128(9), 994–10 0 0. doi:10.1016/j. amjmed.2015.03.023. Desrosiers, J., & Lübbecke, M. E. (2005). A primer in column generation. In G. Desaulniers (Ed.), Column generation. In GERAD 25th anniversary series: 5 (pp. 1– 32). New York, NY: Springer. doi:10.1007/0- 387- 25486- 2_1. Desrosiers, J., & Lübbecke, M. E. (2011). Branch-price-and-cut algorithms. In J. J. Cochran (Ed.), Wiley encyclopedia of operations research and management science. Hoboken, NJ: Wiley Interscience. doi:10.1002/9780470400531.eorms0118. Ellman, M. S., Tobin, D. G., Stepczynski, J., & Doolittle, B. (2016). Continuity of care as an educational goal but failed reality in resident training: Time to innovate. Journal of Graduate Medical Education, 8(2), 150–153. doi:10.4300/ JGME- D- 15- 00278.1.

Erhard, M., Schoenfelder, J., Fügener, A., & Brunner, J. O. (2018). State of the art in physician scheduling. European Journal of Operational Research, 265(1), 1–18. doi:10.1016/j.ejor.2017.06.037. Ernst, A., Jiang, H., Krishnamoorthy, M., & Sier, D. (2004). Staff scheduling and rostering: A review of applications, methods and models. European Journal of Operational Research, 153(1), 3–27. doi:10.1016/S0377-2217(03)0 0 095-X. European Commission (2012). Commission staff working document on an action plan for the EU health workforce. Franz, L. S., & Miller, J. L. (1993). Scheduling medical residents to rotations: Solving the large-scale multiperiod staff assignment problem. Operations Research, 41(2), 269–279. German College of General Practitioners and Family Physicians (2009). Speciality training for general practice in Germany: A report by a panel of invited international experts. German Medical Association (2003). Specialty training regulations. http://www. bundesaerztekammer.de/. Guo, J., Morrison, D. R., Jacobson, S. H., & Jokela, J. A. (2014). Complexity results for the basic residency scheduling problem. Journal of Scheduling, 17(3), 211–223. doi:10.1007/s10951- 013- 0362- 9. Katzenbach, J. R., & Smith, D. K. (1993). The wisdom of teams: Creating the high-performance organization. Boston, Mass.: Harvard Business School Press. KC, D. S. (2014). Does multitasking improve performance? Evidence from the emergency department. Manufacturing & Service Operations Management, 16(2), 168– 183. doi:10.1287/msom.2013.0464. LAG Baden-Württemberg (2015). Judge in 11.09.2015 (1 Sa 5/15), openJur 2015, 19265. Miani, C., Hinrichs, S., Pitchforth, E., Bienkowska-Gibbs, T., Disbeschl, S., Roland, M., & Nolte, E. (2015). Best practice: Medizinische Aus- und Weiterbildung aus internationaler Perspektive. RAND Corporation. Ozkarahan, I. (1994). A scheduling model for hospital residents. Journal of Medical Systems, 18(5), 251–265. doi:10.10 07/BF0 0996605. Proano, R. A., & Agarwal, A. (2017). Scheduling internal medicine resident rotations to ensure fairness and facilitate continuity of care. Health Care Management Science. doi:10.1007/s10729- 017- 9403- 9. Roos, M., Blauth, E., Steinhäuser, J., Ledig, T., Joos, S., & Peters-Klimm, F. (2011). Gebietsweiterbildung Allgemeinmedizin in Deutschland: Eine bundesweite Umfrage unter Ärztinnen und Ärzten in Weiterbildung. Zeitschrift für Evidenz, Fortbildung und Qualität im Gesundheitswesen, 105(2), 81–88. doi:10.1016/j.zefq.2010.11.007. Schlitzkus, L. L., Schenarts, K. D., & Schenarts, P. J. (2010). Is your residency program ready for Generation Y? Journal of Surgical Education, 67(2), 108–111. doi:10.1016/j.jsurg.2010.03.004. Sherali, H. D., Ramahi, M. H., & Saifee, Q. J. (2002). Hospital resident scheduling problem. Production Planning & Control, 13(2), 220–233. doi:10.1080/ 09537280110069667. Smalley, H. K., & Keskinocak, P. (2014). Automated medical resident rotation and shift scheduling to ensure quality resident education and patient care. Health Care Management Science. doi:10.1007/s10729- 014- 9289- 8. Soyster, A. L. (1973). Technical note-convex programming with set-inclusive constraints and applications to inexact linear programming. Operations Research, 21(5), 1154–1157. doi:10.1287/opre.21.5.1154. Topaloglu, S. (2006). A multi-objective programming model for scheduling emergency medicine residents. Computers & Industrial Engineering, 51(3), 375–388. doi:10.1016/j.cie.20 06.08.0 03. Van den Bergh, J., Beliën, J., de Bruecker, P., Demeulemeester, E., & de Boeck, L. (2013). Personnel scheduling: A literature review. European Journal of Operational Research, 226(3), 367–385. doi:10.1016/j.ejor.2012.11.029.