Improving Planning and Scheduling of Appointments at a Rehab Clinic

(1)

IN

DEGREE PROJECT MATHEMATICS, SECOND CYCLE, 30 CREDITS

STOCKHOLM SWEDEN 2018,

Improving Planning and

Scheduling of Appointments at a Rehab Clinic

KARL STIERNMAN

KTH ROYAL INSTITUTE OF TECHNOLOGY SCHOOL OF ENGINEERING SCIENCES

(2)

(3)

Improving Planning and

Scheduling of Appointments at a Rehab Clinic

KARL STIERNMAN

Degree Projects in Optimization and Systems Theory (30 ECTS credits) Degree Programme in Aerospace Engineering (120 credits)

KTH Royal Institute of Technology year 2018 Supervisor at CIF: Giampaolo Martufi

Supervisor at KTH: Per Enqvist Examiner at KTH: Per Enqvist

(4)

TRITA-SCI-GRU 2018:401 MAT-E 2018:84

Royal Institute of Technology School of Engineering Sciences KTH SCI

SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

(5)

Improving Planning and Scheduling of Appointments at a Rehab Clinic

SF281X Master Thesis in Systems Engineering

Karl Stiernman

Abstract

An online simulated model of the routine at a rehab clinic was produced. Even if the model would be used only for long-term planning some results could still be relevant. It is also useful for simulating the response of the system when testing alternative scenarios.

The subject of this thesis cannot be easily modeled or defined as any standard operations research problem. It belongs to the very general category of scheduling. Some subproblems include employee scheduling and a variant of the bin packing decision problem. Related processes and practices are queuing theory, resource allocation, model building and simulation.

The main problem was to schedule a number of activities, part of assessments or rehab programs, that were to be allocated resources in the form of personnel and locations. Each of the referrals and patients were to be assigned to assessments or programs belonging to a certain rehab team. A regression model was developed to predict the arrival rate of the referrals to the clinic. The scheduling of programs was solved by the use of metaheuristics.

The purpose of the simulation modeling was to study the effect of certain scheduling strategies. The results are open to interpretation since almost all measures have tradeoffs. Only the management could decide what should be put into practice. The special cases investigated in this work could be part of the decision.

(6)

(7)

Förbättring av planering och schemaläggning av bokade tider vid en rehabiliteringsklinik

SF281X Masteravhandling i systemteknik

Karl Stiernman

Sammanfattning

En realtidsanpassad simuleringsmodell framställdes, vilken beskriver verksamheten på en rehabiliteringsklinik. Även om modellen enbart skulle avändas för långtidsplanering skulle vissa resultat fortfarande kunna vara relevanta. Den är också användbar för att simulera systemets uppförande när alternativa scenarion testas.

Det går inte att beskriva den framtagna matematiska modellen som ett standardproblem inom opera- tionsanalys. Det tillhör den övergripande kategorin schemaläggningsproblem. Delproblem inbegriper tilldelning av anställda och ett beslutsproblem närbesläktat med flera olika packningsproblem. Andra läror och områden som utnyttjats är köteori, resursfördelning, modellbygge och simulering.

Kärnproblemet utgjordes av att schemalägga ett antal aktiviteter, som i sin tur tillhörde bedömnings- fasen eller ett rehabiliteringsprogram, och som skulle tilldelas resurser i form av anställda och lokaler.

Varje remiss och patient går igenom bedömning eller rehabiliteringsprogram tillhörande ett visst reha- biliteringsteam. En regressionsmodell togs fram för att förutse mängden remisser till kliniken en given dag. Vid schemaläggningen av rehabiliteringsprogram användes en metaheuristisk metod.

Syftet med simuleringen var att analysera resultaten av olika schemaläggningsstrategier. Resultaten är flertydiga eftersom i princip alla åtgärder leder till både förbättringar och försämringar. Enbart led- ningen bestämmer vad som ska tillämpas i praktiken. Specialfallen som behandlas i denna avhandling skulle kunna vara en del av beslutsprocessen.

(8)

(9)

1. Introduction

This research was carried out with the conditions and needs of the pain rehab clinic at the hospital Dan- deryds sjukhus in special mind. The objective evolved over time, but there were two main sub-targets: to create a scheduling application which would help the staff at the clinic, and secondly, to create a model that would allow for analysis of the rehab clinic as a system. The scheduling must take all the many im- plicit constraints into account. It was a very specific problem, but some subproblems are well-known in operations research.

Appointments are referred to as activities and may be part of a rehab program. The problem was made difficult because it was uncertain how to prioritize the programs. Some scheduling schemes were selected and then the best one was implemented in the final model. Discriminating between the models was not trivial because each of them led to both improvements and drawbacks depending on what performance metric was considered.

Scheduling is in general an arduous task. Especially the booking of rooms, checking to see if one was available, was inconvenient for the scheduler. The clinic solved the task by having templates with more or less defined presets for each of the programs. The only defined policy for creating the schedule was to try to schedule as many programs as possible during each semester with each program being rigidly defined.

The numbers of each program that could be scheduled in a semester were approximately known.

The biggest proportion of this work went into formulating the mathematical optimization program and transforming it into a computer model. Many alternative representations, naively implementing some discontinuous penalty to increase the quality of the schedule, were revised because it was very hard to control the system in the desired manner. The stability margin of the system, or more specifically for some rehab teams, is not good as the arrival rate is high and the service rate is low.

Initially the personal opinions of the staff members on the design of the new schedule was a part of the problem description. No effort was made to satisfy either the personnel or the patients individually in this work apart from necessary general constraints. The model was however developed in cooperation with the director and the scheduler, and much effort was spent on correctly modeling the management of the clinic. Such cooperation was a prerequisite for this project.

The problem at hand is NP-hard. Therefore the use of exhaustive search is infeasible and one cannot expect to obtain an optimal solution. With that said, the aim was to find an acceptable solution to the problem and try to generalize the results if possible. Most problematic was the dependency constraints within programs. Therefore a version of adaptive simulated annealing was implemented, which was used to schedule new programs every week. Thus the model is of online type. The solution to the full problem was obtained by simulation.

The word customer is throughout this thesis used in the operations research sense of the word and de- notes either a referral or a patient depending on the context. The personnel at the clinic are the servers.

Rehab introduction programs and continuation programs will be alternatively called short and long programs, respectively.

1.1. Conceptual Analysis

Scheduling is a very broad subject in operations research. The purpose of this section is to investigate previous relevant work on the topic to gain a better understanding of what methods may produce a solution to the problem. To avoid confusion it is necessary to establish the difference between planning and scheduling. Planning is synonymous with developing a strategy, to conceive a plan to achieve some goal over a longer period of time[1]. Scheduling, on the other hand, deals with more detail the allocation of resources, satisfying constraints, and so on, and is more short-term in scope [1]. Although the tasks are very much related, they are often solved separately using different methods [1]. However, solving parts of them simultaneously has great advantages. By making the planner more intelligent or foreseeing one

(12)

1. Introduction 2

can reduce the amount of backtracking needed by the scheduler when the resources are not adequate [1].

Constraint programming is a known technique in scheduling but may in fact also be used for problems in planning [1].

Starting off, the paper by Crown et al. (2017) gives an apt overview of optimization methods with applications in healthcare. The modeling stage outlined in the paper is followed, albeit not strictly—from structuring of the problem, mathematical formulation and model development, to validating the model.

Next is the optimization stage. Its steps are divided into first choosing an optimization method, then proceeding with the optimization and performing sensitivity analyses, to reporting of the results and to ultimately use the results for decision making [2].

Constraint programming is used for solving combinatorial optimization problems which are often defined as so called constraint satisfaction problems (CSPs) [3]. A constrained optimization problem consists of its objective function, its decision variables to be determined by the solution, and its constraints [2].

Parameters are constants used in the objective function and in the constraints [2]. Constraints may be explicitly or implicitly defined. They are not only employed for ensuring that the solution is correct, but also actively used for removing incompatible values from the domains and lead to new deducible constraints [3]. This is known as constraint propagation.

The complexity of the model is dependent on many things, e.g. the desired output and the level of detail [2]. Distinction is made between linear and nonlinear, static and dynamic, and between stochastic and deterministic models. The decision variables may be either discrete/integer and/or continuous. These attributes affect the methods available for solving the problem. Static models do not depend on time, whereas for dynamic models the optimization problem itself and its parameters may be time-dependent, and the decision made now may affect decisions at future time points [2]. In any real problem stochastic (i.e. uncertain) information is to be expected, but sometimes its effects can be neglected in the model.

In addition to the different branches of optimization there also exist simulated systems. For these a range of different scenarios are evaluated instead of finding the optimal solution [2]. There are different views on how the real system is best represented. In an agent-based model entities are autonomous, interactive agents that can make decisions based on their state and their intentions or behaviors [4]. An agent is capable of communicating with its environment and responds to the changes therein [4].

An alternative to agent-based modeling and simulation (ABMS) is the more traditional discrete-event simulation (DES). Discrete events are a way to handle the time aspect of the model—in this case a next- event technique [5]. The model is only updated when a change of state, known as an event, is known to occur. The technique has the consequence that the model needs to be more intelligent and store more information about the future in exchange for not having to unnecessarily check the state of the model [5].

The entities are either passive or active (sources are conflicting on this notion; some require them to be passive, although the terms may have different definitions) and the system state is a result of the states of each individual entity [5]. Resources are a collection of identical and countable items that are not modeled individually [5]. Distinguishing between entities and resources is an active choice of the modeler [5]. A few terms helping with the organization of entities are classes, sets and attributes. Classes are groups of similar entities, sets are temporary groups that can represent states or queues, and attributes belong to each individual entity that sets it apart from the others and controls its behavior [5]. In addition to events, other terms that control the flow of the simulation are activities, processes and the simulation clock. Activities alter the state of the entities, processes are chronologically ordered sequences of events, and the simulation clock keeps track of the simulated time so that activities are executed at the appropriate moment as the result of a Boolean expression [5].

DES is especially suitable for modeling queuing systems [5]. A three-phase approach to DES is the preferred option according to Pidd (2004). Activities are split into two groups known as bound activities (also simply known as Bs) and conditional activities (Cs). The Bs are already scheduled, i.e. to be executed at a known time instance, while the Cs are not. The Cs may depend on the fulfillment of some condition. This does not mean that all of the Bs are necessarily executed before the Cs. A set of records are

(13)

1. Introduction 3

needed for each entity in order for this to work. These records are known as the time cell, the availability and the next activity. The time cell indicates when the entity is due to change state, if known. The availability is a Boolean field which is true if the entity is uncommitted; in such a case the time cell is meaningless. If false, the time cell indicates the entity’s next change of state. The next activity gives the next B if the availability is true. It is meaningless if the availability is false [5].

The phases are labeled as A, B and C. During phase A the log is scanned for when the next event is due. The clock is advanced to the corresponding time point and held constant until the cycle begins anew (the next A phase). During the B phase, the Bs that are due at this particular time point are executed. All the Cs are tested during phase C in accordance with their priority. If its conditions are fulfilled the C is executed, which may cause that other Cs with lower priority cannot be scheduled if they are dependent on the same resources. However, Cs may also release resources which make previous Cs available for execution. Therefore it is important to repeatedly scan the Cs until they all sequentially report negative.

The cycle is repeated until the simulated time reaches the intended run length of the simulation [5].

While ABMS employs bottom-up modeling, the modeling approach in DES is top-down [4, 6]. In contrast to DES, for ABMS the agents can have either a discrete or continuum setup and the results of the system can be more than the sum of parts [4]. Both of them are denoted as dynamic simulation models [2].

Compared to the simpler but slower so called activity-scanning approach, three-phase has the advantage of separating the Bs and Cs, and thus not needing to scan through all the Bs at each event time as the Bs that are not due are ignored [5].

Resource constraints limit the resources available for a particular activity at a point in time [1]. Com- patibility constraints, for example that only professions of a certain team can supervise an activity, are included in this category. Then there are transition constraints that bind the variables at different points in time, and dependency constraints that describe the relations between resources, and bind the variables to different resources during different time points. Precedence constraints are included in the latter category as a special case [1], as they bind the starting point of an activity to the completion of another activity within a program.

In a task-centric dynamic model it is easier to express the dependency constraints a priori, while it is more difficult to do so for resource and transition constraints since they are dependent on the allocation of the activities [1].

Sometimes mathematical problems are of a nature that finding the optimal solution would require ex- treme computational resources, at least with today’s technology. For these cases it is preferable to employ heuristics. Heuristics are not guaranteed to be exact but can be used to find a solution that is good enough.

Metaheuristics are a more general form of heuristics that can be used for a wide range of problems.

1.2. Problem Structuring

The main problem was to schedule a number of appointments. Each appointment is an activity taking place at the hospital. These activities are generated from a set of activity templates based around a few categories of customers. Customers may also be assigned to programs that consist of several activities.

Activities belonging to a program and the programs themselves are subject to other constraints. Customers can be categorized according to which state in the system they currently are:

• Referral

• Assigned to assessment

• Accepted referral

• Rehab introduction program patient

• Eligible for continuation courses

(14)

1. Introduction 4

• Rehab continuation program patient

• Scheduled for follow-up

• Finished treatment

The customer flow is mapped in Figure 1. A referral needs to be assessed before taking part in a program and is considered available for scheduling the day after his referral date. The teams require up to six visits to assess a potential patient. In a real setting, customers may be accepted or rejected already after a few visits, but as the schedule needs to be robust there was no way to exploit this, and so all six visits must be accounted for. Overbooking (although perhaps not of much help in this case) was suggested but was not allowed.

If the referral is accepted he should be assigned to an introduction program. A patient having finished the introduction program may be eligible for entering a continuation program. Otherwise he has finished the treatment. At any point during a program the patient may abort the treatment (not modeled); they are also categorized by the clinic as having finished the treatment. Rejected referrals and persons who have finished treatment leave the system. The customer goes through three phases: assessment, rehab programs, and possibly follow-up. Most follow-ups and evaluations are scheduled immediately after or even before the end of a program, which here will be considered as part of the program. For one continuation program the follow-up is scheduled at the nearest valid date after 90 days has past since the end of the program.

Those patients are still counted as in the system, although the follow-up is a bound activity. Contrarily, it is essential to explicitly keep track of individuals who are eligible for, but have not started, a continuation program, as it is uncertain when that program will be scheduled and who the other selected patients are.

In addition to the customers, each of the activities needs to be allocated a certain number of personnel and a location. The personnel and most of the locations are unary resources—they may only be scheduled to one activity at a time. A few of the locations are assumed to have infinite capacity, e.g. for activities taking place outdoors. Each of the personnel has a profession and belongs to one or more teams. Each team handles certain patients depending on their type of disability.

In charge of each activity are members of the personnel corresponding to a certain team and the required professions. If an activity requires some profession to be present and there are several members with the same profession in the same team, then only one of them needs to be assigned to the activity. In reality, doublet professions in the same team mean that the personnel are working part time. Instead of a constraint that enforces equal work, they are chosen at random if both are available. Over time they will participate in an equal amount of activities, but there are still other factors that if taken into account would generate a more fair schedule for them.

The personnel have normal daily working hours from 08.00 to 16.30. A working day is distributed between scheduled activities, reporting of the same activities, and attending meetings. Usually they have a one-hour lunch break starting at 12.00, the only exception being if lunch is included in the activity (e.g. for activities taking place in the kitchen)—those activities are constrained to lie within a certain time interval around noon. Activities may be finished at 15.30 at the latest, but the last hour can be used for reporting.

A year is divided into a spring and a fall semester starting in January and in August, respectively. Start- ing with the next semester, one whole day (Tuesdays) will be dedicated to meetings and staff development instead of spreading them over the week, and thus they need not be part of the model. Tuesdays, Saturdays and public holidays are exempted from the time domain of validity for scheduling of the activities.

Essential to the question of long-term optimality is if the state health care guarantees (HCGs) are fulfilled. Applied to this case, a referral should be scheduled to her first visit within 30 days from the referral date, and an accepted referral should take part in the introductory courses of a rehab program within 90 days from the date the positive decision on treatment was made [7]. The clinic interprets the 90 days HCG as the period beginning with the date of the accepted referral’s first visit (see e.g. [8]), which will also be the case in this project.

(15)

1. Introduction 5

referral

assessment

accepted?

introduction program

finished

continue?

continuation program

finished

follow-up?

follow-up and evaluation

finished

finished Phase I: Assessment

Phase II: Rehab programs

Phase III: Follow-up

yes

no

yes

no

yes

no

Figure 1: Flowchart of customers.

(16)

1. Introduction 6

There are currently (since 2016) no monetary incentives for fulfilling the HCGs [9]. Still, there is a loss of prestige associated with each violation of the guarantees and it also does not look good for potential customers as this metric is recorded in the annual reports for the hospital. Most of the measures and metrics are at an average level among rehab clinics in the country. Statistics show that in 2017 the number of first visits within 30 days amounted to 93 % of the total (300 of 324), and the share of patients who received treatment within 90 days from their first visit was 84 % (158 of 189).

The ethics behind the use of monetary incentives to reduce queues, the so called queue billion, was explored by Sven Ove Hansson (2014). He wrote that although it has helped reduce queues where it has been applied, the total health care resources are limited, and so one should expect that queues are adversely affected elsewhere [10]. The number of patients not receiving treatment within the guaranteed limits could be used in the objective function which is minimized with respect to this quantity. The question remains if they should be. Naively optimizing with respect to the fulfillment of HCGs (e.g. using Moore’s algorithm) could mean that referrals who have already been denied their first visit within 30 days are neglected for a considerable amount of time. One could argue for and against whether a schedule where one referral has to wait 31 days and the other 11 days is better than a schedule where both have to wait 29 days. One advantage of the latter is that it is more equal. Maybe the disparity was too large in this case, but what an either-or scenario would look like is left as an open question.

A reasonable assumption is that unless all unfavorable cases were accounted for it would be detrimental to bother about the HCGs when scheduling, instead encouraging improvements elsewhere. The HCGs remain important metrics when considering the overall effectiveness of the planning and scheduling and so are still reported. It would be preferable to somehow increase the overall flow of customers. Strictly optimizing the flow of customers, i.e. maximizing the number of customers who are deemed finished with their treatment, would for finite time horizons, and especially if the arrival rate is higher than the service rate, lead to prioritizing the handling of assessments and short programs (which have relatively large capacity) as they are finished more quickly. It is neither considered optimal to maximize the number of visits nor the number of individual visitors per day—which also are important measures.

The optimal flow of customers should depend on a number of internal and external factors, such as the current number of customers in the system and their state, acceptance ratios into the programs, the expected number of reneging customers (if a more complicated model is used), among other things. Reneg- ing is not modeled separately, but the acceptance ratios for each phase are. Even though the queue lengths are not directly visible for customers, balking is a possibility if the queues are rumored to be too long.

There are no statistics available particular to reneging or balking.

In this work a very simple strategy was implemented. The plan is only dependent on the populations of each customer state. As has been the rule hitherto, programs for the same team are constrained not to overlap (except possible follow-ups which are defined separately). While this has a negative impact on the flow, it also simplifies the scheduling. Programs are not started until reaching full capacity. In the base case, the order of priority is that if a continuation program can be filled, schedule it; if not, and if an introduction program can be filled, schedule it instead. Else, schedule only assessments. Assessments may be scheduled during programs. Customers are prioritized based on their referral date or date of first visit (if applicable), such that older dates are picked first. In the long run the number of short programs should only be slightly higher than the number of long programs as they have similar capacities, but only a certain fraction of those who finish short programs continue treatment.

Since the customers should be informed in advance of taking part in an activity, the invitation should be sent at least one weekend before the start date of assessments and two weekends before the start of programs. Rescheduling is possible before the invitation has been sent out, which is assumed to occur at the aforementioned point in time. Since the model is by necessity dynamic, the schedule should be reconsidered each Friday. With all the conditions and constraints taken into account, not much is gained from scheduling too far into the future. This limits the amount of activities to be scheduled each time.

A problem is that the clinic is closing for three weeks in July-August. There is also a break in the

(17)

1. Introduction 7

scheduling of rehab programs around two weeks before and until two weeks after the summer break.

Referrals in July will have to wait a fairly long time. It may be impossible to schedule a few of the referrals in due time. Fortunately, the referral rate goes down during those two months. Programs and assessments must be finished before the end of each semester.

Contrary to the above simple plan, an objective function, which is to be minimized, needs to be mathematically defined for the scheduling problem. There are many alternatives, but one usual quantity is the makespan. For programs it is stated that some activities are to be scheduled a certain number per week.

It is hard to implement because the total number of activities would then change with the length of the program. Various constraints were considered that limited the number of activities scheduled per week and per day for programs. After having considered a wide range of different quality factors the only other constituent chosen was the minimum number of activities per scheduled day for programs. Specifying a maximum amount would be detrimental to the makespan and in a way almost defeat the purpose of optimizing the schedule.

There ought to be some delay when traveling between different locations. This was initially modeled, but deemed unnecessary and removed as the delay is negligible. Neither is cleanup time of the previous activity considered. Activities are often scheduled consecutively, implying that some time will be spent on traveling or arranging of the new activity. Whereas activities could end at any time, they can only start at discrete time points separated by 15 minutes.

As participants have to take time off work for attending the courses, activities in programs are assumed to be only scheduled a maximum of one to three weekdays per week. The number of activities per week varies, but it is preferable to schedule some minimum number of activities per week in order to have some control over the time span of the program. It is also a matter of service quality.

The objective function could include a discount. The reason for implementing a discount term is that in a real scenario there are many uncertainties. For example, personnel may fall ill or some unforeseen event (a black swan) may even require a shutdown of the clinic for some time. It is therefore preferable to complete each activity as soon as possible since it is better having to reschedule only a few activities, and also because in the meantime more things can go wrong the more activities that are left. This term contradicts the makespan and it is difficult to balance them. In the scenarios considered it was set to zero and removed.

Many decision variables, parameters and constraints have been identified in this section. See Table 1 in the next section for a better overview of them. Several constraints and superfluous quantities were removed and disregarded in the final model. Individual appointments with patients have not been modeled. They are rare and happen only occasionally. They do not compete for locations but bind the personnel and may be dependent on other activities.

1.3. Structure of the Model

Although formulating the mathematical problem and developing the model was an iterative process, beginning with the former, it helps with the explanation of the problem solving process to first give an overview of the model (see Figure 2).

Judging from the problem structuring, DES modeling seemed like a good match. If the model had been more complex, e.g. if it had included a set of interacting entities, then ABMS could have been the better alternative. As an example, in a very complicated model there could be some mechanism affecting the relationships between the personnel and the patients depending on the design of the schedule, which in turn may make the patient renege on the treatment. An overworked staff member is likely to not perform as well, for example. Maybe the conclusion would be that it would be preferable to increase the dead time between activities. The available data is not enough to make such deductions or to describe such mechanisms.

(18)

1. Introduction 8

input data

manage input

activity generator

resource allocator

feasible?

implement schedule

finished?

save results

output Preprocessor

Planner

Scheduler

Possible revision

Implementation

Update clock

Postprocessor

yes

no

yes no

Figure 2: Overview of the model.

(19)

2. Mathematical Modeling 9

2. Mathematical Modeling

From the description of the problem a mathematical model can be formulated for the scheduling task. The more complex constraints can be partitioned into subproblems. First a number of sets, parameters and decision variables are defined in Table 1. One-index will be used in all algorithms. Some of the notation is quite unnatural in mathematics but is encountered in algebraic modeling languages.

Table 1: Definitions of sets, parameters and decision variables for the optimization problem.

Set / class Description

D(P) Valid days per program; sets change dynamically with the finishing time

I Activities

J Personnel

K Customers

K_intro(T ) Class containing customers available for introduction programs per team K_cont(T ) Class containing customers available for continuation programs per team

L Locations

P Programs

R Professions

T Teams

Z Valid dates and time points

Z_s(I) Class of valid start times per activity, Zs(I) ⊆ Z ∀i ∈ I Z_f(I) Class of valid finish times per activity, Zf(I) ⊆ Z ∀i ∈ I Parameter Description

γ Cost/activity deficiency/program incurred during days for which the minimum number of activities per day for a program is not fulfilled A^i,l_loc Logical matrix describing if location l ∈ L may be assigned to activity

i ∈ I

A^j,z_pers Logical matrix describing if personnel j ∈ J may be scheduled during time point z ∈ Z (could for example model illness, meetings, etc.) B_pers^j,r,t Logical matrix describing if personnel j ∈ J has profession r ∈ R and

is member of team t ∈ T

C_loc^l Maximum number of simultaneous activities at location l ∈ L

C_intro^t Maximum number of customers that may be assigned to the intro

program for team t ∈ T

C_cont^t Maximum number of customers that may be assigned to the continuation program for team t ∈ T

D^i,r,t_pers Logical matrix describing if activity i ∈ I needs personnel of profession r ∈ R and if it belongs to team t ∈ T

D_prog,min^p Minimum number of activities per week for program p ∈ P

M Some sufficiently large parameter (M ≈ 1000)

Oⁱ^α^,i^β Logical matrix describing if activity iα ∈ I should be scheduled after i_β ∈ I

Tⁱ Time length of activity i ∈ I

Continued on next page

(20)

Variable Description

τ Makespan

bî,j,z₁ Binary variable; is one if y_persî,j and yî,z_timeare both one, else zero bî,k,z₂ Binary variable; is one if y_custî,k and y_timeî,z are both one, else zero bî,l,z₃ Binary variable; is one if y_locî,l and y_timeî,z are both one, else zero

tⁱ_s Scheduled start date and time of activity i ∈ I tⁱ_f Finish date and time of activity i ∈ I

x^p,d(p)_prog,min Positive integer variable; describes how many activities below the minimum are scheduled during day d ∈ D(P) for program p ∈ P ypers^i,j Binary variable; is one if personnel j ∈ J is assigned to activity i ∈ I,

else zero

y^i,k_cust Binary variable; is one if customer k ∈ K is assigned to activity i ∈ I, else zero

y^i,l_loc Binary variable; is one if location l ∈ L is assigned to activity i ∈ I, else zero

y^i,z_time Binary variable; is one if activity i ∈ I is scheduled during time point z ∈ Z, else zero

y^i,p,d(p)_prog Binary variable; is one if activity i ∈ I belongs to program p ∈ P and is scheduled during day d(p) ∈ D(P)

Allocating customers is part of the scheduling task, but the activities are in turn already defined or generated by what types of customers are in the system. It is assumed that a maximum of one program per team is scheduled at a time, although in the computer model a loop makes sure that more programs can be scheduled for the same team after the previous has finished, based on some conditions.

It is not necessary to keep track of customers individually since the activities themselves, in programs or assessments, can be constrained not to overlap. However, since any number of assessments may be scheduled, it is easier to formulate the problem by constraining the customers not to be assigned to more than one activity at a time. Going through each customer would increase the computation time though.

With the quantities defined, an attempt at formulating an optimization problem of the scheduling task is the following:

(21)

min f (·) := τ + γX

p∈P

X

d(p)∈D(P)

x^p,d(p)_prog,min (1)

subject to (2)

tⁱ_f = tⁱ_s+ Tⁱ∀i ∈ I (3)

τ ≥ tⁱ_f ∀i ∈ I (4)

x^p,d(p)_prog,min≥ D^p_prog,min− X

i∈I(P)

yprog^i,p,d(p)∀(i, p, d(p)) ∈ I × P × D(P) (5)

bî,j,z₁ ≥ −1 + y_persî,j + y_timeî,z ∀(i, j, z) ∈ I × J × Z (6)

bî,k,z₂ ≥ −1 + y_custî,k + y_timeî,z ∀(i, k, z) ∈ I × K × Z (7)

bî,l,z₃ ≥ −1 + y_locî,l + yî,z_time∀(i, l, z) ∈ I × L × Z (8)

X

i∈I

b^i,j,z₁ ≤ A^j,z_pers∀(j, z) ∈ J × Z (9)

X

i∈I

b^i,k,z₂ ≤ 1 ∀(k, z) ∈ K × Z (10)

X

i∈I

b^i,l,z₃ ≤ C_loc^l ∀(l, z) ∈ L × Z (11)

X

j∈J

B_pers^j,r,ty^i,j_pers= D^i,r,t_pers ∀(i, r, t) ∈ I × R × T (12)

X

l∈L

2 − A^i,l_loc

y_loc^i,l = 1 ∀i ∈ I (13)

tⁱ_s^α ≥ Oⁱ^α^,i^βtⁱ_f^β ∀(i_α, i_β) ∈ I × I (14)

y_cust^i,k =

(1 if customer k is assigned to activity i

0 otherwise (15)

y_prog^i,p,d(p)=

(1 if activity i in program p is scheduled during day d(p) ∈ D(P)

0 otherwise (16)

y_time^i,z =

(1 ∀z ∈ [tⁱ_s, tⁱ_f)

0 otherwise (17)

tⁱ_s∈ Z_s(I) ∀i ∈ I (18)

tⁱ_f ∈ Z_f(I) ∀i ∈ I (19)

x^p,d(p)_prog,min∈ Z⁺∀(p, d(p)) ∈ P × D(P) (20)

y_persî,j , yî,k_cust, y_locî,l, yî,z_time, y_progî,p,d(p)∈ {0, 1} (21)

(SUB): see subsection 2.1 Reporting of Activities. (22)

Note that the above mathematical program is instantiated at a point in time and not the formulation of the full problem as it is dependent on time. The program is necessarily a nonlinear one in order to ensure that assigned resources neither switch between activities from one singular time point to another, nor remain permanently fixed to only one particular activity.

The constraint (5) could be changed to x^p,d(p)_prog,min≥ _M¹ P

i∈I(P)y^i,p,d(p)_prog

D^p−P

i∈I(P)y_prog^i,p,d(p)

so

(22)

that the sets in D(P) need not be dynamic, but then it would not penalize the deviation from the minimum amount of activities and result in no penalty at all for days with no scheduled activities.

The valid number of dates and time points is infinite on a half-closed interval. It is not possible to store all valid time points in memory. Instead they are calculated dynamically. For efficiency, all the public holidays except Sundays are calculated once per year of simulation and stored in memory. Easter is predicted by Computus.

It should be emphasized that the program is not solved to optimum. The objective function is retained and used together with metaheuristics detailed in subsection 2.4. To come up with feasible schedules a heuristic is applied that says that the activities, one by one, should be scheduled at the earliest feasible date. Thus only the ordering of the activities matter.

The planner is not particularly intelligent but it is a prerequisite for the scheduler and will therefore be described next in Algorithm 1. There is definitely potential for improvement. If the scheduler has determined that a program cannot be scheduled, then a flag is used to signal that the program is to be ignored for the remainder of the semester. If the planner had been more sophisticated, for example by formulating a CSP and utilizing a heuristic, then infeasible cases could be backtracked to the planner instead of to the scheduler. Generating activities is better handled by the scheduler since it would only result in unnecessary overhead for the planner.

Algorithm 1 The planner.

Input: A flag signaling what programs may be scheduled before the end of the semester, e_sem Output: A flag to signal what programs should be scheduled, e_plan

1: function ACTIVITYSELECTOR

All the bits in e_planare initially assumed reset

All the bits in esemare set at the beginning of each semester

2: for each team t ∈ T do

3: if card (K_cont(t)) ≥ C_cont(t)∧ a new continuation program may be started within the next two weeks (check programs in progress, available weekdays and valid start dates, and if esem(2t) is set) then

4: e_plan(2t) ← ¬e_plan(2t)

5: end if

6: if card (Kintro(t)) ≥ Cintro(t)∧ a new introduction program may be started within the next two weeks then

7: e_plan(2t − 1) ← ¬e_plan(2t − 1)

8: end if

Else do nothing . Signals that only assessments will be scheduled for the team

9: end for

10: return eplan 11: end function

2.1. Reporting of Activities

The personnel responsible need to report the activities. It is required that some time is dedicated to this task the same day as the activity itself. Furthermore, the reporting of a certain activity cannot be interrupted until finished. Another constraint is that the reporting cannot occur before the scheduled activity. The question to be answered is if there exists enough idle time in order for all the scheduled activities to be reported during the same day.

This is similar to the well-known bin packing decision problem, except that only activities preceding the binmay be assigned to it. The bins are variable-sized. The bin packing problem is quite old with several

(23)

invented heuristic techniques, e.g. best fit. There also have been some recent outstanding contributions such as by Rothvoß (2013) who showed that one formulation of the bin packing problem can be solved in polynomial time for any constant number of different item sizes (see [11]). Below is an attempt to model the problem mathematically with m number of activities, or items, and n time slots or bins. It is assumed that the activities and time slots are sorted chronologically. For each personnel, the problem becomes

max

m

X

i=1 n

X

j=1

x^i,j (23)

subject to (24)

m

X

i=1

T_repⁱ x^i,j ≤ T_slot^j ∀j = 1, . . . , n (25)

n

X

j=1

2 − A^i,j x^i,j ≤ 1 ∀i = 1, . . . , m (26)

x^i,j ∈ {0, 1} , (27)

where T_repⁱ is the report time of activity i and T_slot^j is the length of time slot j. The matrix A^i,j ensures only preceding items can be assigned to each bin. The binary variable x^i,j is equal to one if activity i is reported during time slot j, or otherwise equal to zero. The second constraint could also be modeled with a strict equal sign which would give rise to either a solution or an infeasible result. Considering that the original problem is a decision problem, the solution should return either true or false. The result would be feasible if and only if the objective value would be equal to m.

The binary linear program can be solved optimally, for example by the branch and cut method. However, since it is estimated that the number of activities reported per day for each of the personnel is low, the problem was solved by an exhaustive depth-first search with backtracking (see Algorithm 2).

The computations can be made more efficient by the use of bit arrays. The average behavior can be affected by removing time slots that are too small and to set a minimum istart value if all i preceding some j can be assigned to it. The result would be false if backtracking finds no possible combinations for activity indices greater than the minimum istart.

2.2. Assigning Locations

In practice the assignment of a location to an activity is handled after a date has already been suggested for scheduling of the activity. This subproblem can be formulated as a CSP. One only needs to check activities scheduled during the same day. The domains of the activities are the possible locations that can be assigned to them. The constraints are that all activities that take place simultaneously need to be assigned different locations. This is a variant of what is generally known as the alldifferent constraint in the literature. The result can also be infeasible in this case. Then the conclusion should be that the suggested date should be updated to when the first conflicting activity ends and all dependent constraints need to be checked again—including this one.

Let D(X_i) be the domains of the activities. Pairs are made of activities scheduled simultaneously and listed in E. The constraints are formulated as C = {∀ (Xi, X_j) ∈ E, Location (X_i) 6= Location (X_j)}.

The problem has been transformed into checking that each pair has different locations assigned for the members. The pairings can be seen as edges in a dependency graph. There exist many solution methods for these kinds of problems including heuristics. One popular method is depth-first search together with backtracking which remains valid in this case. Improvements can be made by forward checking and constraint propagation. It has been proven that the alldifferent constraint together with precedences in

(24)

Algorithm 2 Reporting of activities.

Input: Scheduled activities, S, new activity, a, and proposed time, t Output: Boolean value, true = feasible

Each report time is denoted by T_repⁱ

Available time of slot j during the day is denoted by T_slot^j

It is assumed that the activities and time slots are in chronological order

1: function CHKREPORTFEASIBILITY(S, a, t)

2: for each responsible personnel p do

3: Make a list of each activity i in S that is scheduled the same day as a (check t)

4: Construct domains Dⁱ for each activity i, meaning that for each j in Dⁱ the ith activity pre- cedes time slot j

5: Initialize dⁱ = 1 ∀i = 1, . . . , m . First alternative at each node

6: Initialize x^i,j = 0 ∀i = 1, . . . , m, j = 1, . . . , n

7: i_start= 1 . Starting index

8: repeat

9: i = i_start

10: if T_repⁱ ≤ T_slot^j(dⁱ⁾−P_m

α=1x^α,j(dⁱ⁾T_rep^α then

11: x^i,j(dⁱ⁾= 1

12: istart← i_start+ 1

13: else if dⁱ < card Dⁱ then

14: dⁱ← dⁱ+ 1

15: else

16: b = true

17: for activity k = i − 1, . . . , 1 do

18: if d^k < card D^k then

19: x^l,r = 0 ∀l = i, . . . , k, r = 1, . . . , n

20: b = false

21: d^k← d^k+ 1

22: d^l= 1 ∀l = k + 1, . . . , m

23: i_start= k

24: break

25: end if

26: end for

27: if b then . Is true when backtracking not possible

28: return false

29: end if

30: end if

31: untilPm i=1

Pn

j=1x^i,j= m . Feasible for personnel p

32: end for

33: return true

34: end function

(25)

general is NP-hard [12].

An interesting algorithm tailored to precedence problems is published in [12]. The alldifferent constraint is altered to include the pairs of dependencies, E. It requires preprocessing the domains of the variables so that the precedence constraints Xi < X_j are respected: min(Xi) ≤ min(X_j) and max(Xi) ≤ max(Xj) for each pair. The algorithm then enforces this bounds consistency via propagation and prunes values that are incompatible.

Having considered the alternatives, a simple depth-first search with backjumping was chosen. It has the advantage that when there are no more alternatives for the current node, the solution jumps back to the preceding node without having to check all indices inbetween. The implementation was quite straightforward as only locations with the capacity of one activity at a time were considered.

2.3. Predicting the Influx of Referrals per Day

The predicted number of referrals should naturally be a nonnegative and discrete function. To this end a Poisson or a negative binomial (NB) time-series regression model seemed like good candidates for modeling of the available count data with respect to unknown parameters. The NB distribution is generated by several different stochastic mechanisms [13, 14], for example by certain models of waiting time or by a linear birth and death process [13]. Fundamental to the Poisson distribution is the law of rare events and it is a distribution that occurs frequently in natural phenomena [15]. Two prerequisites for choosing a Poisson model are that the mean and the variance should be the same, and that the observations must be independent. The mean and variance do not need to be the same for the NB model. It is more general than the Poisson model [13].

2.3.1. Initial Data Analysis

The data collected range from the start of the rehab clinic in 2011 to February 2018. In Figure 3 are the run sequence plots of the number of referrals and accepted referrals, respectively, each day. Dates when there were no referrals were not registered but they were added manually. While not directly visible in the referral run sequence plots, the data corresponding to the first two years in the plot over the number of individual visitors to the clinic were qualitatively regarded as outliers in the lower side of the data set. A possible cause could be that there were not that many patients in the system. Additionally, there was an extra phase in the pain rehab program for patients with cerebral injuries. Moreover, the statistics were not complete for these years; for example, reasons for rejection were not given, among other things. Therefore the data from the first two years (2011 and 2012) were discarded from all of the calculations.

Sometimes several referrals are sent for one and the same person (O. Svensson, personal communica- tion, May 9, 2018). A person being already a known patient at the rehab clinic is one of the criteria for rejecting a referral. This directly affects the dependency of the observed number of referrals and also the prediction of the accepted number of referrals. This rejection criterion is however only one of five and it is the rarest of them. The statistics show that 110 out of 4119 referrals were rejected in this way during the period from 3 January 2013 to 6 February 2018, corresponding to 2.7% of all referrals. Fortunately, the statistics show exactly what referral date each rejection (per criterion) pertains to. Thus, in this case the simplest and best strategy for dealing with this category of referrals was to simply subtract them from the total amount of referrals. If this had not been the case an alternative strategy would have been to use a Cox regression model.

Henceforth, referrals will imply only unique referrals. The mean and variance of the (unique) referrals and accepted referrals are in Table 3 distributed over workdays, Saturdays, Sundays and other work- free days by statutory right (that is, public holidays excluding civic holidays and ordinary Sundays).

These groups are mutually exclusive and with holidays having precedence, meaning referrals on work- free holidays are not counted in the other groups.

(26)

(a)

(b)

Figure 3: (a) Referrals and (b) accepted referrals over time with linear trend 2013-2018.

(27)

Table 3: Distribution of referrals and accepted referrals per day.

Workdays Saturdays Sundays Work-free holidays

Total number of days 1281 251 250 82

Total referrals 3931 25 41 12

Referrals, mean 3.06870 0.09960 0.16400 0.14634

Referrals, variance 3.71113 0.12155 0.16910 0.17371

Total accepted referrals 2141 14 22 8

Accepted referrals, mean 1.67135 0.05578 0.08800 0.09756 Accepted referrals, variance 1.78192 0.06063 0.09626 0.13682

It was necessary to distinguish between the days as the means and variances are quite different. Sundays and other public holidays share similar mean and variance while for Saturdays the mean is much lower.

Note that the counted number of days are all days past since the beginning of the series (January 3, 2013) and not just the registered days.

The mean of the number of unique referrals per workday is 3.07 and the variance is higher at 3.71. Given a Poisson model, this would give rise to only modest overdispersion [13]. This is without any regression models. Using the same method as Cameron and Trivedi (2013), namely to go from a specific to a more general model [13], a Poisson model cannot be automatically ruled out. The subsequent analysis, plots and most of the calculations deal only with referrals during workdays.

2.3.2. IID Regression Models

Three classes of distributions are defined in the following. Let yi be the dependent variables (counts) and xibe the independent explanatory (e.g. categorical) variables. The Poisson regression model has the following probability mass function (PMF) for yigiven xi[13]:

f (y_i|x_i) = exp (−µ_i) µ^y_iⁱ

yi! , y_i= 0, 1, 2, . . . (28)

with the mean parameter defined as

E [yi|x_i] = µi = eiexp x^T_i β , (29) where β is the vector of regression parameters and ei is a measure of exposure [16]; in the following ei = 1. The corresponding PMF for a negative binomial model is [13, 16]

f (y_i|µ_i, ν_i) = Γ (y_i+ ν_i) Γ (yi+ 1) Γ (νi)

ν_i νi+ µi

νi µ_i νi+ µi

yi

, y_i = 0, 1, 2, . . . . (30) Depending on how the parameter ν is chosen, different models are obtained. The expected mean is E[Y_i] = µ_iand the variance is V [Yi] = µ_i+ µ^p_iv_i⁻¹. The standard NB formulation, which Cameron and Trivedi (2013) call the NB2 model, is obtained with p = 2 and νi = α⁻¹, with α ≥ 0 being the dispersion parameter. For other models νiis defined as νi = α⁻¹µ^2−p_i . Another relatively common model is the NB1 model with p = 1. If α = 0 the mean and variance will be the same and the above reduces to the Poisson.

For α > 0 the variance exceeds the mean (see [13, 16]).

Another model is the generalized Poisson (GP), see [16]. Two different implementations will be defined,

(28)

called GP1 and GP2. For the GP1, let f (yi|µ_i, α) =

µi

1 + αµ_i

yi (1 + αyi)^yⁱ⁻¹ y_i! exp

−µi(1 + αyi) 1 + αµ_i

, yi = 0, 1, 2, . . . , (31) with mean E[yi] = µi and variance V [yi] = µi(1 + αµi)². The GP2 model is different yet another natural extension of the Poisson (see [16]):

f (y_i|µ_i, α) = (

µ_i(µ_i+ (α − 1) y_i)^yⁱ⁻¹α^−yⁱ^exp(^−α⁻¹^(µi+(α−1)yi))

yi! , y_i = 0, 1, 2, . . . ,

0, y_i > m, α < 1, (32)

where α ≥ max ¹₂, 1 −^µ₄ⁱ, and m the largest positive integer such that µ_i+ m(α − 1) > 0 when α < 1.

The mean is equal to E[y_i] = µ_iand the variance is V [y_i] = α²µ_i. Advantages of the GP model is that it can handle both overdispersion (α > 1) and underdispersion ¹₂ ≤ α < 1 and µi > 2 [16], and that in the particular case that α = 0 for GP1, µ can be interpreted as the Poisson mean [17].

The simultaneous first-order conditions for the max likelihood estimates (MLE) of the parameters were given by Cameron and Trivedi (2013) [13], and by Ismail and Jemain (2007) [16] for the GP distribution with the assumption that the mean or the fitted values are multiplicative. The ML functions can be differentiated with respect to the parameters and the resulting differential equations set to zero in order to maximize the MLE. Ismail and Jemain (2007) suggest iterative weighted least squares and Newton- Raphson iteration as methods for solving of the parameters [16]. As an example, take the GP1 max likelihood function given by [16]:

l (α, β) =X

i

yilog

µi

1 + αµi

+ (yi− 1) log (1 + αy_i) −µi(1 + αyi) 1 + αµi

− log (y_i!). (33) The first and second derivatives with respect to α are:

∂l (α, β)

∂α =X

i

−y_i 1 + α²µ²_i + 2 (α + 1) µ_i + µ²_i + y²_i

(1 + αµi)²(1 + αyi) , (34)

∂²l (α, β)

∂α² =X

i

µ³_iy_i²α³_i + 3µ_i²y²_iα²_i − 3µ_iy_i(y_i(y_i− 2µ_i− 1) + µ²_i)α_i− y_i³+ y²_i + 3µ²_iy_i− 2µ³_i (1 + µiαi)³(1 + yiαi)² .

(35) The parameter α can then be solved for using Newton-Raphson iteration. The different models were tested under the independent identically distributed (IID) assumption with no regression parameters and the fitted parameters estimated via max likelihood (see histogram in Figure 4). The NB and GP models perform noticeably better than the simple Poisson but the underlying reasons for this are unclear. Generally speaking there is unobserved heterogeneity in the data. Choosing a more complex model requires that the improvement of the fit is significant enough. Thankfully there are objective measures commonly used in statistics for discriminating among the models. The first is the Akaike information criterion (AIC):

AIC = −2l + 2p, (36)

where l is the log likelihood and p the number of independently adjusted parameters in the model [18].

Another common measure is the Bayesian information criterion (BIC):

BIC = −2l + p log (n) , (37)

(29)

where l and p are defined as for AIC, and n is the number of observations [19, 16]. The scores in the AIC and BIC are relative to the other models, and the lower score a model receives the better it is. There are different philosophies behind the measures but AIC cannot be said to be non-Bayesian [20]. BIC was derived from AIC and it leans more toward lower-dimensional models [19]. There are circumstances where AIC outperforms BIC and others where BIC is better [20].

Figure 4: Empirical frequency and histograms of IID models naively fitted to all of the data with no regression parameters.

The processes affecting the referral rate are not known. Box plots were produced to investigate the monthly distribution and the distribution per weekday (see Figure 5). There seem to be only small dif- ferences for the referrals distributed per weekday. All of the weekdays have the same median of three referrals. The monthly distribution does show some variance and further investigation needs to be con- ducted.

The count data distributed per weekday are in Table 4. Let the hypothesis H0 be that a referral may be sent with equal probability each weekday (equiprobable model). Pearson’s chi-squared test was applied to test the hypothesis. The same notation as in Karl Pearson’s 1900 paper was used [21].

In Table 4, χ² = 5.28975. The more accurate G-test for goodness of fit gives the very similar value G = 5.35789. Four degrees of freedom and a nominal significance level of 0.05 give a critical χ² value of 9.48773. Thus, H0 should not be rejected and the decision was made to not include weekdays in a regression model.

With inspiration from Zeger (1988) who modeled the outbreak of the polio virus in the USA 1970-83 by the use of sine and cosine pairs [22], the seasonal variations were as a first test estimated by a polynomial because it was unclear what periods to choose. In order to avoid overfitting the statistics were split into one part comprising the test data (Jan 2013-Feb 2017) and the other part the validation data (Feb 2017-Feb 2018). Let τw denote the number of workdays (starting at zero) from the beginning of the series and τwp

be the number of workdays past since the start of the year. There were approximately 250 workdays per year in the series. The explanatory variables were initially chosen as

x =1 τ_w τ_wp τ_wp² . . . τ_wp⁶ T

.

The parameter τwis the trend and the ‘1’ is an intercept term; this is necessary for estimating the mean

(30)

(a)

(b)

Figure 5: Box plots of referrals (a) per weekday and (b) per month. Outliers are denoted by dots outside the boxes.

Improving Planning and Scheduling of Appointments at a Rehab Clinic

Improving Planning and

Scheduling of Appointments at a Rehab Clinic

KARL STIERNMAN

Improving Planning and

Scheduling of Appointments at a Rehab Clinic

KARL STIERNMAN

Improving Planning and Scheduling of Appointments at a Rehab Clinic

SF281X Master Thesis in Systems Engineering

Karl Stiernman

Förbättring av planering och schemaläggning av bokade tider vid en rehabiliteringsklinik

SF281X Masteravhandling i systemteknik

Karl Stiernman

Contents

1. Introduction

2. Mathematical Modeling