- Email: [email protected]
A simulated annealing approach to the solution of job rotation scheduling problems
Applied Mathematics and Computation 188 (2007) 31–45 www.elsevier.com/locate/amc
A simulated annealing approach to the solution of job rotation scheduling problems Serap Ulusam Sec¸kiner
a,*
, Mustafa Kurt
b
a
b
Faculty of Engineering, Department of Industrial Engineering, Gaziantep University, 27310 Gaziantep, Turkey Faculty of Engineering and Architecture, Department of Industrial Engineering Gazi University, 06570 Maltepe, Ankara, Turkey
Abstract This paper presents a new solution to the job rotation scheduling problem where the objective is to minimize the workload for each worker. Our motivation for this study comes from dangerous jobs that have some exposures. Job rotation is one method that is sometimes used to reduce exposure to strenuous jobs. Especially, this method can be applied to the service organizations that customer demand varies over the course of an operating day and across the days of an operating week. However, developing effective rotation schedules can be complex in even moderate sized service systems. Integer programming and a simulated annealing algorithm were used to construct the schedules. The efficiency of simulated annealing in solving combinatorial optimization problems is very well known. However, it has recently not been applied to job rotation scheduling problem based on the review of the available literature. In this research paper, the developed models are presented and results for test problems are reported. 2006 Elsevier Inc. All rights reserved. Keywords: Job rotation scheduling; Integer programming; Simulated annealing
1. Introduction Many industries and public services function around the clock, 7 days a week. Well-known examples are police or fire departments, security guard companies, automobile assembly plants, hospitals, etc. In a manufacturing or service environment, job scheduling is an essential component in maintaining the profitability of the enterprise. Within a schedule, the timing and sequencing of operations that humans and machines perform must strike a balance between productivity demands (i.e. delivery of service of the certain type at a certain time and place) and safety concerns of the personnel involved in meeting these demands [1]. Job rotation is one method that is sometimes used to reduce exposure to strenuous jobs. To adequately address this issue, schedules for certain human-machine operations can be developed that include systematic variations in task exposure hours of the involved workers. In this study, job rotation schedules with varied labor demand are
*
Corresponding author. E-mail address: [email protected] (S.U. Sec¸kiner).
0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.09.082
32
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
designed to ensure that the demands of the job do balance among workers. Sometimes, although workers perform same type work, some workers expose occupational risks more than others. For example, radiography technicians expose the radiation dose with different level according to type of radiography. In order to balance the existing radiation dose, job rotation method is a good way. Carnahan et al. [1] published a key paper on modeling and solving job rotation scheduling problem (JRSP) with using integer programming and genetic algorithm. Since then, no researchers proposed new or improved models and algorithms to solve the problem. Job rotation method has been studied different views by the other researchers. Generally, the advantages and disadvantages of this method have been discussed. Allwood and Lee [2] analyzed the effects of job rotation schemes in lean production systems. Paul et al. [3] discussed that job rotation is a factor in reducing physical workload at a refuse collecting department. The results indicate that the total amount of work performed by means of job rotation resulted in an overall reduced physical workload of the employees of the refuse collecting department. Ghare and Leutenegger [4] showed that an order of magnitude reduction in mean job slowdown could be achieved by rotating a user’s batch jobs within their queue. Becker et al. [5] identified the stressors associated with each task, and thus determine whether a job rotation strategy involving document preparation could offset the psychological and physiological stress of prolonged VDT work. Pickavance [6] emphasized most organizations that are faced with day-to-day operational changes which result in fluctuations in workload, shifts in priorities and/or alterations of schedules and deadlines. He reported that available staff resources also vary as a result of absence due to holidays, sickness, training and other reasons. His study sets out to examine the various utilization and deployment procedures and practices different organizations employ to manage these changes, in particular the use made of computerized systems. According to Wilkinson and Edwards [7], stabilization of shift times and rotation of jobs within a shift has enabled a 2-man system of continuous manning to outperform a 3-man arrangement which lacked these features. They reported that in vigilance-type tasks calling for sustained attention it is probably the job rotation which is important; in more intense cognitive tasks such as complex decision-taking it may be the stabilization of daily shift times which confers the advantage, with job rotation doing more harm than good. The work of Carnahan et al. [1,8,9] has generally been accepted as the first serious approach to model and solve JRSP. He developed optimal solution and heuristic solution methodology for JRSP using integer programming and genetic algorithm approach. This problem has similar features with the assignment problem. Given that the assignments of tasks to workers within each time period are independent of one another, the total number of potential job rotation schedules available to the planner is quite large. With using enumeration to evaluate all possible job rotation schedules is impractical. The problem of assigning the jobs to labors falls within the category of problems known as combinatorial optimization problems. Traditionally, these types of problems have been solved by using integer programming (IP). Unfortunately, it is difficult to use integer programming to develop job rotation schedules when there are dynamic work characteristics. Nowadays, a set of powerful heuristic computational procedures, known as simulated annealing (SA) that are designed to search and find answers to complex problems in situations where the number of possible alternative solutions is vast. The efficiency of simulated annealing in solving combinatorial optimization problems is very well known. However, it has recently not been applied to job rotation scheduling problem based on the review of the available literature. The rest of this paper is composed of three sections. In Section 1, we provide a brief review of job rotation scheduling. In Section 2, we present the modified integer programming model and the SA heuristic for JRSP. Section 3 displays experimental results and the paper concludes with a discussion of the research implications. 2. Problem environment and solution approaches 2.1. The integer programming model for JRSP The work of Carnahan et al. [1] is the first serious approach to model and solve to JRSP. We modified their integer programming model based on a tour scheduling environment that labor requirements across a 12-hday, 7 days a week operating environment was to be satisfied full-time labors.
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
33
The model tries to minimize the maximum W (workload cost) found across the schedules of each of the labors for a service facility that operates less than 24 h and 7 days a week. The model assumptions are specifically as follows: 1. 2. 3. 4. 5. 6.
An organization operates 12 h a day. The planning horizon of the model is 7 days a week. Each employee must receive n off-days each week. A shift type is defined by a start time, shift length. Full-time shift types are allowed. There are m work of type that has different working cost (workload cost), changeable according to day and work. 7. No understaffing is allowed. 8. All employees have same qualifications. The IP model can be given as Minimize
Z¼W
ð1Þ
Subject to : N X
X nmik ¼ Dmik
i ¼ 1; . . . ; 7;
k ¼ 1; . . . ; K;
m ¼ 1; . . . ; M;
ð2Þ
n M X K X k
m2M
i2I
k2K
m2M
i2I
k2K
XXX XXX
cmik X nmik 6 W X nmik 6 r
X nmik 2 ð0; 1Þ i k n m r Ki Dmik Xnmik cmik W
i ¼ 1; . . . ; 7;
X nmik 6 1
m
n ¼ 1; . . . ; N ;
n ¼ 1; . . . ; N ; n ¼ 1; . . . ; N ;
cmik 2 R ðreal numberÞ:
ð3Þ ð4Þ ð5Þ ð6Þ
index that shows the days of a week shift type index that shows the employee index that specifies the number of work of type number of assignments for each employee number of feasible shift patterns on day i number of employees assigned to shift pattern k on day i for performing work of type m 0–1 binary value where xnmik = 1 if employee n is assigned to shift pattern k for performing work of type m on day i and zero otherwise the cost variable is a real number which represents the cost, in terms of workload, of any employee performing work of type m in shift pattern k on day i workload cost
In this study, integer program was comprised of n (number of worker) · m (work of type) · 7 (days of a week) · k (shift type) = n · m · 7 · k decision variables. The purpose of the integer program was to find values for the n · m · 7 · k decision variables so that the maximum W found across the schedules of workers was minimized. The general formulation for the integer program is shown above. The cost variable cmik is a real number, which represents the cost, in terms of workload, of any worker performing work of type m in shift pattern k on day i. The first set of constraints (2) ensures that sufficient number of workers is available to satisfy the labor requirements for work of type m in shift k on day i. The second set of constraints (3) ensures that each worker performs only one work in each time period or not. The third set of constraints (4) considers the overall work demands placed upon all workers based on the job assignments and
34
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
calculates the W value for each worker. The minimization objective forces the integer program to find the set of Xnmik values corresponding to the job rotation schedule that minimizes the maximum W found among the some workers. Thus, the maximum W being for some workers is balanced among all workers. The fourth set of constraints (5) ensures that each worker is assigned to job amount of work demands. This constraint supplies that each employee must receive n off-days each week. The last constraint set (6) places zero–one variable. Unfortunately, the combinatorial structures of proposed model make the solution difficult and time consuming. Since the Xnmik represent the workers assigned to job, they are typically required to be 0–1 integervalued. When integer restrictions are placed on theXnmik, the problem is NP-complete and thus very resistant to efficient solution [10]. Therefore, feasible solution searching with heuristically is needed instead of finding optimal solution. Simulated annealing (SA) is a well-established stochastic neighbourhood search technique. It has a potential to solve complex combinatorial optimization problems [11]. Because of this reason, we evaluate the performance of the simulated annealing technique for JRSP in order to the workload balancing. The simulated annealing technique has received ever growing attention and has proven to be a very effective approach to combinatorial optimization problems since 1983. For detailed information about this technique, see Reeves [12]. 2.2. The simulated annealing algorithm for JRSP Simulated annealing (SA) is a heuristic search procedure for combinatorial optimization problems. It was initially proposed by Kirkpatrick et al. [13] based on the physical sciences. It has the capability of jumping out of the local optima for global optimization. This capability is achieved by accepting with probability neighbouring solutions worse than the current solution. The acceptance probability is determined by a control parameter (temperature) which decreases during the SA procedure. It has been successfully applied to several management science related problems including: (1) facility layout problem [11]; (2) flexible labour scheduling problem [10], and (3) school timetabling [14]. 2.2.1. Parameters of the SA algorithm We define (under the presumption of a minimization objective) the following terminology for our description of simulated annealing: X f(X) X0 X* d Tin Tf a fmin fmax Pc Pf
an initial solution to an optimization problem S (standard deviation of an initial solution) a neighbouring solution toX the minimum standard deviation solution found f(X 0 ) f(X) initial temperature final temperature or freezing temperature the cooling ratio (a < 1) the lower bound for a given JRSP the higher bound for a given JRSP the initial acceptance probability the final acceptance probability
A solution to the problem: A configuration or solution is represented as a 2d matrix in this implementation of SA where each row represents a worker and each column represents a job. Demand (work of type, day and shift type) is placed at the intersections. A rotation is comprised from works that were assigned to a worker. Fig. 1 represents a sample configuration for an illustrative problem with fourteen workers and four parts. The schedule environment is two types of work, four days, five shift and six planning period. This problem is formulated with using tour scheduling problem [15]. This gives us the optimal number of workers and worker requirements in concerned planning periods in order to use in our problem. Then GAMS/Cplex 8.1 Solver is used to solve it. Finally, GAMS solution file gives the following data:
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
35
Fig. 1. Representation of the solution vector in SA implementation of JRSP.
Optimal number of worker (N): 14. Demands: D111 = 7 (c111=250), D121 = 7 (c121=200), D131 = 7 (c131=200), D141 = 7 (c141=150), D211 = 7 (c211=300), D221 = 7 (c221=240), D331 = 7 (c331=240), D241 = 7 (c241=180). In this example, the workload cost variables (cmik) was assigned as improvised numbers and given in parenthesis with demand. As seen from Fig. 2, the maximum workload is 920 for some workers. The aim of searching with heuristically is to find the assignments that prove workload balancing. In this example, the total workload is 12320 and the workload of each worker is desired as 880 (12320/14). The best solution of this example with heuristically is given in Appendix A.1. In Fig. 1, demand D241 is located in rotation 1 in worker 1, demand D121 is located in rotation 1 in worker 2 and so on. If we represent the solution vector formally as SOLns then SOL11 = D241 and SOL12 = D121, etc. The workload cost of each worker is calculated by summation of cost that related to demand. The total workload cost of a solution is the summation of all workers’ workload costs as shown in the right lower bottom of Fig. 1. In this work, solution quality of the problem is observed by standard deviation (S) of assigned workload costs for each worker. The best solution is obtained when standard deviation is zero. Neighbourhood move: A transition from one configuration to another that result in neighbour solution. Neighbour solutions are generated by a swapping procedure that is defined as follows: Step 1: Select the worker (row) who has the lowest workload cost randomly. Step 2. Select the worker (row) who has the highest workload cost randomly. Step 3. If workers have not off-days, determine two rotation part (column) that have same days randomly. If workers have off-days, determine two rotation that have same days and not shared days randomly. For example: Worker A: D214 D425 D145 D151 D234.
36
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
Fig. 2. Neighbourhood move. (i) Before neighbourhood move, (ii) after neighbourhood move.
Worker B: D433 D415 D373 D122 D253. Possible swaps between days: 1–1, 2–2, 3–3, 5–5, 4–7. An example of neighbourhood move for the above illustrative problem is shown in Fig. 2. In Fig. 2 cost values are filled randomly to show the effect of a neighbourhood move. Randomly selected workers and rotations are bolded. Present swapping strategy always generates feasible solutions and can be implemented easily. Objective function: A measure of how good a solution is. Standard deviation (S) of worker’s workload cost as shown in the right lower bottom of in Figs. 1 and 2 constitutes the objective function of the SA procedure. Acceptance/rejection of a solution: The neighbourhood solution is accepted if the objective function improves, otherwise the solution is accepted with probability depending on the temperature, which is set to allow the acceptance of a large proportion of generated solutions at the beginning (high temperature). Then, the temperature is modified (lowered) to decrease the probability of acceptance. At each temperature a predetermined number of moves are attempted. Assume that i and j are two solutions with standard deviation
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
37
f(i) and f(j), respectively. Then the acceptance criterion determines whether j is accepted from i by applying the following acceptance probability: 1; if f ðjÞ 6 f ðiÞ; P c faccept jg ¼ ððf ðjÞf ðiÞÞ=T Þ ; if f ðjÞ > f ðiÞ; e T 2 Rþ denotes control parameter ðtemperatureÞ:
ð7Þ
The purpose of sometimes accepting an inferior solution is to avoid becoming stuck in a local optimum. Termination: If the previously determined parameters are reached the algorithm terminates, e.g. maximum number of iterations, final temperature, no improvement for a number of iterations, etc. Cooling schedule: Selection of the SA parameters is called a cooling schedule. This gives how high the starting temperature should be, and rules to determine it. This also gives (a) when the current temperature should be lowered, (b) by how much the temperature should be lowered, and (c) when the annealing process should be terminated. 2.2.2. Determination of the cooling schedule The cooling schedule of the SA algorithm for the job rotation scheduling problem is determined as follows. Initial temperature (Tin): As mentioned above, the implementation of SA is based on the process of annealing of solids. The initial temperature must be set sufficiently high, resulting in a high degree of randomness in the initial stages of the search. Thus most of the initial movements are accepted. The value of initial temperature is regulated in such a way that at the beginning of SA the value of acceptance probability (Pc) is very closer 1. In this application it is set to 0.99. The value of acceptance probability (Pc) is determined by using factorial design. An estimation for the initial value of the temperature can be made by using the following equation [11]: fmin fmax T in ¼ ; ð8Þ ln P c where fmin and fmax are lower and higher bounds for a given JRSP, respectively. For this application, fmin = total workload/number of workers and fmax = 1.10 · (total workload/number of workers) are selected. Length of the Markov chain at iteration (L): It is chosen in such a way that for each temperature a thermal equilibrium can be achieved. But as the initial temperature approaches to zero, the acceptance probability of movements decreases. As a result of this, the length of Markov chain goes to infinity. Because of this situation, L is bounded by a constant which is dependent on the problem size. In this work, it is defined by L = (number of worker). Rate of cooling (a): The temperature should be lowered slowly in order to attain slow cooling. The temperature is decreased by using the following relation; Tin = aTn1, where a is the cooling rate its value being between zero and one. The initial temperature is estimated based on the initial acceptance probability, the final temperature (Tf) can also be estimated similarly based on the final acceptance probability (Pf). Therefore, cooling rate is also dependent on the values of initial and final acceptance probabilities and the number of temperature reduction cycles (t) (number of iterations). The definition of Bennage and Dhingra [16] for cooling rate is employed in this application: 1=ðt1Þ ln P c a¼ : ð9Þ ln P f Pf approaches to zero near to the end of the annealing process. It must be set very close to zero in order to avoid accepting non-improving solutions. It is set to 1 · 1015 in this application. Termination criterion: When the maximum number of iterations (t) is reached the SA algorithm terminates. t can be determined in such a way that at the end of the last iteration, Tf can be reached. In order to check this, the value of final temperature can be calculated by using Eq. (10). If the value of the final temperature calculated using Eq. (10) is not very close or equal to estimated final temperature then t must be reselected and the rate of cooling (a) should be recomputed. This procedure should be repeated until the required number of iterations (t) is determined. T f ¼ T in at :
ð10Þ
38
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
2.2.3. The pseudo code of the proposed SA algorithm for the JRSP problem The pseudo code of the proposed SA algorithm for solving JRSP is given as follows: GENERATE the feasible solution with X (generate the solution randomly) Calculate the objective function (f(X) = S); Set the initial temperature (Tin); Set the number of iteration (t); Set the cooling ratio (a); Set the initial acceptance probability (Pc); Set the final acceptance probability (Pf); Set the lower and higher bounds (fmin, fmax); t = 0; REPEAT l=0 REPEAT Generate a neighbouring solution as explained above and calculate the new objective function value. d = f(X 0 ) f(X); s 0 IF(d 6 0) THEN X = X ; ELSE produce u (random number (0, 1)); 0 IF (u < exp(d/T)) THEN X = X ; l = l + 1; UNTIL l = L; ‘‘number of workers’’; t = t + 1; Tin = aTn1 ; UNTIL ‘‘Termination criterion (t)’’; Best solution (X* the minimum standard deviation solution found), f(X*) and W; 3. Computational results 3.1. Scheduling environment The environment used in this research is based on a tour scheduling problem originally utilized by Topa¨ zkarahan [15]. In this prior study, labor requirements across a 12-h-day, 7 days a week operating log˘lu and O environment was to be satisfied full-time employees working five consecutive days per week. Shifts were eight hours in length with 1-h meal break between third and tenth hours of work. Shift start times are five for each day of work. In this formulation, at any time the number and duration of break periods can be changed. Additionally, the feasible set of shift types, the location of break windows, and the feasible work patterns for the allowable shift lengths can be easily be specified by the scheduler, and thus the model can be adapted to desired scheduling environment conveniently. The objective is to minimize the number of employees scheduled subject to the coverage of labor requirements. In this study, job rotation scheduling is performed based on tour scheduling environment. We exploit the tour scheduling environment to make our SA more meaningful. Thus we assumed that minimum number of employees (N) and labor requirements (Dmik) are known. For this application, we generated randomly a varying demand pattern for each work of type within the hours of a day and across the days of a week. The generated demand patterns were given in Appendix A.2. 3.2. Heuristic and computational time performance We describe our computational experiments about JRSP. We generated the test problems based on combinations of types of work. This gives us 26 different situations (see Table 2). This will prove us an idea of the effectiveness of the algorithms with different dimensions of rotations.
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
39
In order to test the IP model, one run was performed for each situation. In order to test the SA, five replications were performed for each situation. A limitation of integer programming is that it is designed to optimize deterministic problems or problems where the data is known with certainty. Unfortunately, the test problems described in this section does not have this characteristic; because of real world situation were not considered. One way to handle this difficulty is to use improvised cmik values that change according to work types and days for test problems. This permits the development of the necessary coefficients to construct both the IP and SA. The improvised workload cost of each type of work and day is given in Table 1. Furthermore, it is easy to calculate total workload cost for each problem when numbers of worker and labor requirements across a 12-h-day are determined. Additionally, Table 2 gives numbers of worker and the total workload cost for each problem. GAMS/Cplex Solver 8.1 has been used to solve the test problems on a Pentium IV 2.6 Ghz computer. All test problems has been written in GAMS software and solve CPLEX solver. The search for the optimal solution has been limited to 86,400 CPU s. IP formulation has not found optimal solutions for all test problems. At the end of time limit, a near optimal solution was obtained. This was used for comparing with SA algorithm. For SA algorithm, an initial experiment was done to test the different parameters. We ran tests with iteration number of 1000, 1500, and 2000. We also performed experimental runs with length of the Markov chain at iteration (number of workers) of L/2, L, and 2L. Finally we tested with acceptance probability of 0.90, 0.95, and 0.99. Based on the results, we found that the following parameter values gave the best results for SA algorithm: iteration number = 1000, length of the Markov chain at iteration = L and acceptance probability = 0.99. The IP model and SA results are shown in Table 3. In the IP results, it is not easy to calculate the standard deviation of a schedule. Because of this reason, the results are compared according to objective function value (near optimal solution) for each test problems. Table 3 gives the minimum workload cost solution for each of the 26 problems using the following methods:
Table 1 The improvised workload costs Work of type
Monday
Tuesday
Wednesday
Thursday
Friday
Saturday
Sunday
1 2 3 4 5
250 300 120 47 18
200 240 120 47 18
200 240 80 47 18
150 180 80 35 12
150 120 60 24 12
100 120 40 12 6
100 90 40 12 6
Table 2 Descriptions of the test problems Problem no.
Type of works (m)
Number of worker
Total workload
Problem no.
Type of works (m)
Number of worker
Total workload
Problem no.
Type of works (m)
Number of worker
Total workload
1 2 3 4 5 6 7 8 9
1–2 1–3 1–4 1–5 2–3 2–4 2–5 3–4 3–5
39 36 36 31 23 23 18 20 15
35,000 26,210 23,805 22,424 17,240 14,835 13,260 5855 4496
10 11 12 13 14 15 16 17 18
4–5 1–2–3 1–2–4 1–2–5 1–3–4 1–3–5 1–4–5 2–3–4 2–3–5
15 49 49 44 46 41 41 33 28
2091 38,720 34,860 35,416 27,865 27,439 24,191 18,995 17,598
19 20 21 22 23 24 25 26
2–4–5 3–4–5 1–2–3–4 1–2–3–5 1–2–4–5 1–3–4–5 2–3–4–5 1–2–3–4–5
28 25 59 54 54 51 38 64
15,153 6239 40,945 39,570 37,165 28,139 19,331 41,311
40
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
Table 3 Results of the computational study with 26 test problems Problem name
(1) Rotation between job 1 and 2 (2) Rotation between job 1 and 3 (3) Rotation between job 1 and 4 (4) Rotation between job 1 and 5 (5) Rotation between job 2 and 3 (6) Rotation between job 2 and 4 (7) Rotation between job 2 and 5 (8) Rotation between job 3 and 4 (9) Rotation between job 3 and 5 (10) Rotation between job 4 and 5 (11) Rotation between job 1, 2 and 3 (12) Rotation between job 1, 2, and 4 (13) Rotation between job 1, 2, and 5 (14) Rotation between job 1, 3, and 4 (15) Rotation between job 1, 3, and 5 (16) Rotation between job 1, 4, and 5 (17) Rotation between job 2, 3 and 4 (18) Rotation between job 2, 3 and 5 (19) Rotation between job 2, 4 and 5 (20) Rotation between job 3, 4 and 5 (21) Rotation between job 1, 2, 3 and 4 (22) Rotation between job 1, 2, 3 and 5 (23) Rotation between job 1, 2, 4 and 5 (24) Rotation between job 1, 3, 4 and 5 (25) Rotation between job 2, 3, 4 and 5 (26) Rotation between job 1, 2, 3, 4 and 5
IP Computational time (CPU)
SA Computational time (CPU)
IP Best solution value (W)
SA Best solution value (W) Run 1
Run 2
Run 3
Run 4
Run 5
>86,400
238.36
910
980
950
950
950
1000
>86,400
110.02
730
760
770
780
780
800
>86,400
114.48
682
785
747
762
762
774
>86,400
81.67
750
850
850
850
850
812
>86,400
40.39
760
820
860
800
840
860
>86,400
40.98
684
707
719
797
742
742
>86,400
24.34
762
786
798
930
822
912
>86,400
30.78
314
322
309
312
311
311
>86,400
16.00
318
372
340
326
340
318
>86,400
17.66
153
153
159
153
153
159
>86,400
232.57
840
840
870
890
870
880
>86,400
226.95
822
910
910
910
900
910
>86,400
182.56
870
900
850
920
900
940
>86,400
206.75
647
642
644
724
660
674
>86,400
153.47
700
750
750
750
750
750
>86,400
152.57
621
662
668
685
718
647
>86,400
91.10
607
647
680
622
694
667
>86,400
62.54
666
690
930
728
718
678
>86,400
60.00
551
612
629
646
664
588
>86,400
49.69
274
307
266
334
289
286
>86,400
373.09
752
757
747
754
794
750
>86,400
373.20
806
856
840
820
840
840
>86,400
379.28
736
787
757
802
820
737
>86,400
312.69
606
618
598
658
700
614
>86,400
120.93
542
631
667
611
623
553
>86,400
474.89
690
736
672
764
780
742
IP: Integer programming model (searching for the optimal solution). SA: Simulated annealing algorithm (searching for the best heuristic solution).
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
41
Table 4 Solution quality of SA according to IP Problem no.
SA performance according to IP
Average
The solution quality¼ 1 SAIP IP 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
Run 1
Run 2
Run 3
Run 4
Run 5
0.92 0.96 0.85 0.87 0.92 0.97 0.97 0.97 0.83 1.00 1.00 0.89 0.97 1.01 0.93 0.93 0.93 0.96 0.89 0.88 0.99 0.94 0.93 0.98 0.84 0.93
0.96 0.95 0.90 0.87 0.87 0.95 0.95 1.02 0.93 0.96 0.96 0.89 1.02 1.00 0.93 0.92 0.88 0.60 0.86 1.03 1.01 0.96 0.97 1.01 0.77 1.03
0.96 0.93 0.88 0.87 0.95 0.83 0.78 1.01 0.97 1.00 0.94 0.89 0.94 0.88 0.93 0.90 0.98 0.91 0.83 0.78 1.00 0.98 0.91 0.91 0.87 0.89
0.96 0.93 0.88 0.87 0.89 0.92 0.92 1.01 0.93 1.00 0.96 0.91 0.97 0.98 0.93 0.84 0.86 0.92 0.79 0.95 0.94 0.96 0.89 0.84 0.85 0.87
0.96 0.93 0.88 0.87 0.89 0.92 0.92 1.01 0.93 1.00 0.96 0.91 0.97 0.98 0.93 0.84 0.86 0.92 0.79 0.95 0.94 0.96 0.89 0.84 0.85 0.87 Average
0.95 0.94 0.88 0.87 0.91 0.92 0.91 1.00 0.92 0.99 0.97 0.90 0.97 0.97 0.93 0.89 0.90 0.86 0.83 0.92 0.98 0.96 0.92 0.92 0.84 0.92 0.92
Table 4 gives the solution quality for each problem using the equation that is given in this table. The best performance of SA is highlighted. SA performs better than IP for some problems (8, 13, 14, 20, 21, 24 and 26) and this is meaning that the IP solution is not an optimal. The average performance of SA was found 0.92. This is quite effective. If there is a comparison opportunity in point of minimization of standard deviation with other heuristics in the literature, the performance of SA algorithm could be found better than currently. The SA algorithm is programmed in Visual C++. It is capable of solving JRSP once the required data files are generated and the SA parameters. The data file of problem 2 is given in Appendix A.3. The SA algorithm could find the solution within 10 min for all test problems. The maximum solution time is 474.89 CPU times for problem 26. Average solution time is 160.267. 4. Conclusions In this paper, we examined the use of simulated annealing algorithm in solving the job rotation scheduling problem. SA has been applied to numerous optimization problems successfully. Firstly in this study its application to JRSP is presented. In this study, a SA algorithm with a simple but effective data structure and neighbourhood generation mechanism is proposed for its solution. The proposed SA algorithm is applied to test problems from generated randomly. In these tests, the proposed algorithm is compared with IP model results. Unfortunately, the optimal solution can not be obtained with IP model at end of CPU limit. In order to compare optimal and heuristic results, IP results were admitted the best integer solutions within time limit. The performance of SA was evaluated based on IP solutions for each problem. The computational results showed that the proposed algorithm is quite effective.
42
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
Appendix A A.1. The best solution for example within Section 2.2
1 2 3 4 5 6 7 8 9 10 11 12 13 14
1
2
3
4
Workload cost
D111 D111 D211 D111 D111 D211 D211 D211 D111 D211 D211 D211 D111 D111
D221 D221 D121 D221 D221 D121 D121 D121 D221 D121 D121 D121 D221 D221
D231 D231 D131 D231 D231 D131 D131 D131 D231 D131 D131 D131 D231 D231
D141 D141 D241 D141 D141 D241 D241 D241 D141 D241 D241 D241 D141 D141
880 880 880 880 880 880 880 880 880 880 880 880 880 880 s=0 W = 12320
Total
A.2. Demands for test problems Work of type 1 1 2 1.1 3 6 1.2 3 7 1.3 2 7 1.4 3 7 1.5 4 8 1.6 2 5 1.7 1 3
3 8 8 8 8 8 7 5
4 10 9 9 10 10 8 7
5 14 13 12 12 13 9 8
6 14 13 12 12 13 9 8
7 9 9 9 9 9 9 8
8 9 9 9 9 9 8 8
9 13 11 10 10 12 8 8
10 13 11 10 10 12 8 7
11 13 11 10 10 10 8 8
12 9 9 9 9 9 8 8
Work of type 2 1 2 2.1 3 3 2.2 2 6 2.3 2 5 2.4 3 5 2.5 3 5 2.6 2 2 2.7 1 2
3 4 5 3 5 4 1 2
4 5 6 4 3 6 3 3
5 9 9 9 8 4 4 3
6 8 9 9 5 6 3 2
7 6 9 9 5 9 3 2
8 5 6 9 10 10 3 2
9 7 6 8 6 4 4 2
10 6 7 6 4 4 1 1
11 4 2 5 2 2 1 2
12 2 1 4 1 2 1 1
Work of type 3 1 2 3.1 1 2 3.2 1 2
3 3 3
4 5 4
5 6 6
6 6 6
7 4 4
8 4 4
9 5 5
10 5 4
11 5 4
12 4 4
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
43
(continued) Work of type 3 3.3 1 2 3.4 1 2 3.5 2 2 3.6 1 1 3.7 1 1
2 2 3 2 2
4 4 4 3 2
6 6 6 3 2
6 5 6 3 2
4 4 4 3 3
4 4 4 2 2
5 5 6 2 2
4 4 6 2 2
4 4 5 2 2
4 4 4 2 2
Work of type 4 1 2 4.1 1 2 4.2 1 2 4.3 1 2 4.4 1 2 4.5 2 2 4.6 1 1 4.7 1 1
3 3 3 2 2 3 2 2
4 5 4 4 4 4 3 2
5 6 6 6 6 6 3 2
6 6 6 6 5 6 3 2
7 4 4 4 4 4 3 3
8 4 4 4 4 4 2 2
9 5 5 5 5 6 2 2
10 5 4 4 4 6 2 2
11 5 4 4 4 5 2 2
12 4 4 4 4 4 2 2
Work of type 5 1 2 5.1 1 1 5.2 1 1 5.3 1 1 5.4 1 2 5.5 1 1 5.6 1 1 5.7 1 1
3 1 1 1 2 1 1 1
4 2 2 1 2 1 1 1
5 3 3 2 2 2 2 1
6 3 3 3 2 2 2 1
7 2 2 2 2 2 1 1
8 2 2 2 2 2 1 1
9 2 2 2 2 1 1 1
10 2 2 1 2 1 1 1
11 2 1 1 1 1 1 1
12 1 1 1 1 1 1 1
A.3. Example data file for problem 2 Dmik
cmik
n
111 113 114 115 121 122 123 124 125 131 133 134 135 141 144 145 151
250 250 250 250 200 200 200 200 200 200 200 200 200 150 150 150 150
7 1 5 9 3 5 1 2 9 7 1 1 9 8 2 9 8 (continued on next page)
44
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
(continued) Dmik
cmik
n
154 155 161 164 165 171 173 174 175 311 313 314 315 321 323 325 331 334 335 341 344 345 351 353 354 355 361 364 365 371 373 375
150 150 100 100 100 100 100 100 100 120 120 120 120 120 120 120 80 80 80 80 80 80 60 60 60 60 40 40 40 40 40 40
3 9 7 1 8 3 3 1 8 2 2 1 4 2 2 4 2 2 4 2 2 4 2 1 1 4 2 1 2 1 1 2
References [1] B.J. Carnahan, M.S. Redfern, B. Norman, Designing safe job rotation schedules using optimization and heuristic search, Ergonomics 43 (4) (2000) 543–560. [2] J.M. Allwood, W.L. Lee, The impact of job rotation on problem solving skills, International Journal of Production Research 42 (5) (2004) 865–881. [3] P. Paul, F.M. Kuijer, K. Visser Bart, C.G. Han, Job rotation as a factor in reducing physical workload at a refuse collecting department, Ergonomics 42 (9) (1999) 1167–1178. [4] D. Ghare Gaurav, T. Leutenegger Scott, Improving small job response time for opportunistic scheduling, in: IEEE International Workshop on Modeling, Analysis, and Simulation of Computer and Telecommunication Systems – Proceedings, 2000, pp. 557–564. [5] A.B. Becker, N.G. Swanson, S.L. Sauter, T.L. Galinsky, S. Jones, L. Schleifer, Compatibility of job rotation subtasks in data entry work, in: Proceedings of the Human Factors and Ergonomics Society, vol. 2, 1995, p. 965. [6] R. Pickavance, Employee utilisation and deployment (Inst. of Manpower Studies.) report: IMS-115, ISBN-1-85184-013-3, 1985, p. 39. [7] R.T. Wilkinson, R.S. Edwards, Stable hours and varied work as aids to efficiency (Royal Naval Personnel Research Committee London (England)) report: OES-512, October 1969, p. 9. [8] B.J. Carnahan, M.S. Redfern, B.A. Norman, A genetic algorithm for designing job rotation schedules considering ergonomic constraints, in: Proceedings of the 1999 Congress on Evolutionary Computation-CEC99 (Cat. No. 99TH8406), vol. 2, 1999, pp. 1093– 1098.
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
45
[9] B.J. Carnahan, M.S. Redfern, Automatic planning of job rotation schedules, in: Industrial Engineering Research – Conference Proceedings, 1997, pp. 263–268. [10] M. Brusco, L. Jacobs, A simulated annealing approach to the solution of flexible labor scheduling problems, Journal of the Operational Research Society 44 (12) (1993) 1191–1200. [11] A. Baykasog˘lu, N.N.Z. Gindy, A simulated annealing algorithm for dynamic layout problem, Computers and Operations Research 28 (2001) 1403–1426. [12] C.R. Reeves, Modern Heuristic Techniques for Combinatorial Problems, John Wiley& Sons Inc, Newyork, Toronto, 1993, pp. 288– 298. [13] S. Kirkpatrick, C.D. Gelatt Jr., M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680. [14] D. Abramson, Constructing school timetables using simulated annealing: sequential and parallel algorithms, Management Science 37 (1991) 98–113. _ O ¨ zkarahan, Implicit optimal tour scheduling flexible assignments, Computers and Industrial Engineering 44 (2003) [15] S ß . Topalog˘lu, I. 75–89. [16] W.A. Bennage, A.K. Dhingra, Single and multiobjective structural optimization in discrete–continuous variables using simulated annealing, International Journal for Numerical Methods in Engineering 38 (1995) 2753–2773.
A simulated annealing approach to the solution of job rotation scheduling problems Serap Ulusam Sec¸kiner
a,*
, Mustafa Kurt
b
a
b
Faculty of Engineering, Department of Industrial Engineering, Gaziantep University, 27310 Gaziantep, Turkey Faculty of Engineering and Architecture, Department of Industrial Engineering Gazi University, 06570 Maltepe, Ankara, Turkey
Abstract This paper presents a new solution to the job rotation scheduling problem where the objective is to minimize the workload for each worker. Our motivation for this study comes from dangerous jobs that have some exposures. Job rotation is one method that is sometimes used to reduce exposure to strenuous jobs. Especially, this method can be applied to the service organizations that customer demand varies over the course of an operating day and across the days of an operating week. However, developing effective rotation schedules can be complex in even moderate sized service systems. Integer programming and a simulated annealing algorithm were used to construct the schedules. The efficiency of simulated annealing in solving combinatorial optimization problems is very well known. However, it has recently not been applied to job rotation scheduling problem based on the review of the available literature. In this research paper, the developed models are presented and results for test problems are reported. 2006 Elsevier Inc. All rights reserved. Keywords: Job rotation scheduling; Integer programming; Simulated annealing
1. Introduction Many industries and public services function around the clock, 7 days a week. Well-known examples are police or fire departments, security guard companies, automobile assembly plants, hospitals, etc. In a manufacturing or service environment, job scheduling is an essential component in maintaining the profitability of the enterprise. Within a schedule, the timing and sequencing of operations that humans and machines perform must strike a balance between productivity demands (i.e. delivery of service of the certain type at a certain time and place) and safety concerns of the personnel involved in meeting these demands [1]. Job rotation is one method that is sometimes used to reduce exposure to strenuous jobs. To adequately address this issue, schedules for certain human-machine operations can be developed that include systematic variations in task exposure hours of the involved workers. In this study, job rotation schedules with varied labor demand are
*
Corresponding author. E-mail address: [email protected] (S.U. Sec¸kiner).
0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.09.082
32
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
designed to ensure that the demands of the job do balance among workers. Sometimes, although workers perform same type work, some workers expose occupational risks more than others. For example, radiography technicians expose the radiation dose with different level according to type of radiography. In order to balance the existing radiation dose, job rotation method is a good way. Carnahan et al. [1] published a key paper on modeling and solving job rotation scheduling problem (JRSP) with using integer programming and genetic algorithm. Since then, no researchers proposed new or improved models and algorithms to solve the problem. Job rotation method has been studied different views by the other researchers. Generally, the advantages and disadvantages of this method have been discussed. Allwood and Lee [2] analyzed the effects of job rotation schemes in lean production systems. Paul et al. [3] discussed that job rotation is a factor in reducing physical workload at a refuse collecting department. The results indicate that the total amount of work performed by means of job rotation resulted in an overall reduced physical workload of the employees of the refuse collecting department. Ghare and Leutenegger [4] showed that an order of magnitude reduction in mean job slowdown could be achieved by rotating a user’s batch jobs within their queue. Becker et al. [5] identified the stressors associated with each task, and thus determine whether a job rotation strategy involving document preparation could offset the psychological and physiological stress of prolonged VDT work. Pickavance [6] emphasized most organizations that are faced with day-to-day operational changes which result in fluctuations in workload, shifts in priorities and/or alterations of schedules and deadlines. He reported that available staff resources also vary as a result of absence due to holidays, sickness, training and other reasons. His study sets out to examine the various utilization and deployment procedures and practices different organizations employ to manage these changes, in particular the use made of computerized systems. According to Wilkinson and Edwards [7], stabilization of shift times and rotation of jobs within a shift has enabled a 2-man system of continuous manning to outperform a 3-man arrangement which lacked these features. They reported that in vigilance-type tasks calling for sustained attention it is probably the job rotation which is important; in more intense cognitive tasks such as complex decision-taking it may be the stabilization of daily shift times which confers the advantage, with job rotation doing more harm than good. The work of Carnahan et al. [1,8,9] has generally been accepted as the first serious approach to model and solve JRSP. He developed optimal solution and heuristic solution methodology for JRSP using integer programming and genetic algorithm approach. This problem has similar features with the assignment problem. Given that the assignments of tasks to workers within each time period are independent of one another, the total number of potential job rotation schedules available to the planner is quite large. With using enumeration to evaluate all possible job rotation schedules is impractical. The problem of assigning the jobs to labors falls within the category of problems known as combinatorial optimization problems. Traditionally, these types of problems have been solved by using integer programming (IP). Unfortunately, it is difficult to use integer programming to develop job rotation schedules when there are dynamic work characteristics. Nowadays, a set of powerful heuristic computational procedures, known as simulated annealing (SA) that are designed to search and find answers to complex problems in situations where the number of possible alternative solutions is vast. The efficiency of simulated annealing in solving combinatorial optimization problems is very well known. However, it has recently not been applied to job rotation scheduling problem based on the review of the available literature. The rest of this paper is composed of three sections. In Section 1, we provide a brief review of job rotation scheduling. In Section 2, we present the modified integer programming model and the SA heuristic for JRSP. Section 3 displays experimental results and the paper concludes with a discussion of the research implications. 2. Problem environment and solution approaches 2.1. The integer programming model for JRSP The work of Carnahan et al. [1] is the first serious approach to model and solve to JRSP. We modified their integer programming model based on a tour scheduling environment that labor requirements across a 12-hday, 7 days a week operating environment was to be satisfied full-time labors.
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
33
The model tries to minimize the maximum W (workload cost) found across the schedules of each of the labors for a service facility that operates less than 24 h and 7 days a week. The model assumptions are specifically as follows: 1. 2. 3. 4. 5. 6.
An organization operates 12 h a day. The planning horizon of the model is 7 days a week. Each employee must receive n off-days each week. A shift type is defined by a start time, shift length. Full-time shift types are allowed. There are m work of type that has different working cost (workload cost), changeable according to day and work. 7. No understaffing is allowed. 8. All employees have same qualifications. The IP model can be given as Minimize
Z¼W
ð1Þ
Subject to : N X
X nmik ¼ Dmik
i ¼ 1; . . . ; 7;
k ¼ 1; . . . ; K;
m ¼ 1; . . . ; M;
ð2Þ
n M X K X k
m2M
i2I
k2K
m2M
i2I
k2K
XXX XXX
cmik X nmik 6 W X nmik 6 r
X nmik 2 ð0; 1Þ i k n m r Ki Dmik Xnmik cmik W
i ¼ 1; . . . ; 7;
X nmik 6 1
m
n ¼ 1; . . . ; N ;
n ¼ 1; . . . ; N ; n ¼ 1; . . . ; N ;
cmik 2 R ðreal numberÞ:
ð3Þ ð4Þ ð5Þ ð6Þ
index that shows the days of a week shift type index that shows the employee index that specifies the number of work of type number of assignments for each employee number of feasible shift patterns on day i number of employees assigned to shift pattern k on day i for performing work of type m 0–1 binary value where xnmik = 1 if employee n is assigned to shift pattern k for performing work of type m on day i and zero otherwise the cost variable is a real number which represents the cost, in terms of workload, of any employee performing work of type m in shift pattern k on day i workload cost
In this study, integer program was comprised of n (number of worker) · m (work of type) · 7 (days of a week) · k (shift type) = n · m · 7 · k decision variables. The purpose of the integer program was to find values for the n · m · 7 · k decision variables so that the maximum W found across the schedules of workers was minimized. The general formulation for the integer program is shown above. The cost variable cmik is a real number, which represents the cost, in terms of workload, of any worker performing work of type m in shift pattern k on day i. The first set of constraints (2) ensures that sufficient number of workers is available to satisfy the labor requirements for work of type m in shift k on day i. The second set of constraints (3) ensures that each worker performs only one work in each time period or not. The third set of constraints (4) considers the overall work demands placed upon all workers based on the job assignments and
34
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
calculates the W value for each worker. The minimization objective forces the integer program to find the set of Xnmik values corresponding to the job rotation schedule that minimizes the maximum W found among the some workers. Thus, the maximum W being for some workers is balanced among all workers. The fourth set of constraints (5) ensures that each worker is assigned to job amount of work demands. This constraint supplies that each employee must receive n off-days each week. The last constraint set (6) places zero–one variable. Unfortunately, the combinatorial structures of proposed model make the solution difficult and time consuming. Since the Xnmik represent the workers assigned to job, they are typically required to be 0–1 integervalued. When integer restrictions are placed on theXnmik, the problem is NP-complete and thus very resistant to efficient solution [10]. Therefore, feasible solution searching with heuristically is needed instead of finding optimal solution. Simulated annealing (SA) is a well-established stochastic neighbourhood search technique. It has a potential to solve complex combinatorial optimization problems [11]. Because of this reason, we evaluate the performance of the simulated annealing technique for JRSP in order to the workload balancing. The simulated annealing technique has received ever growing attention and has proven to be a very effective approach to combinatorial optimization problems since 1983. For detailed information about this technique, see Reeves [12]. 2.2. The simulated annealing algorithm for JRSP Simulated annealing (SA) is a heuristic search procedure for combinatorial optimization problems. It was initially proposed by Kirkpatrick et al. [13] based on the physical sciences. It has the capability of jumping out of the local optima for global optimization. This capability is achieved by accepting with probability neighbouring solutions worse than the current solution. The acceptance probability is determined by a control parameter (temperature) which decreases during the SA procedure. It has been successfully applied to several management science related problems including: (1) facility layout problem [11]; (2) flexible labour scheduling problem [10], and (3) school timetabling [14]. 2.2.1. Parameters of the SA algorithm We define (under the presumption of a minimization objective) the following terminology for our description of simulated annealing: X f(X) X0 X* d Tin Tf a fmin fmax Pc Pf
an initial solution to an optimization problem S (standard deviation of an initial solution) a neighbouring solution toX the minimum standard deviation solution found f(X 0 ) f(X) initial temperature final temperature or freezing temperature the cooling ratio (a < 1) the lower bound for a given JRSP the higher bound for a given JRSP the initial acceptance probability the final acceptance probability
A solution to the problem: A configuration or solution is represented as a 2d matrix in this implementation of SA where each row represents a worker and each column represents a job. Demand (work of type, day and shift type) is placed at the intersections. A rotation is comprised from works that were assigned to a worker. Fig. 1 represents a sample configuration for an illustrative problem with fourteen workers and four parts. The schedule environment is two types of work, four days, five shift and six planning period. This problem is formulated with using tour scheduling problem [15]. This gives us the optimal number of workers and worker requirements in concerned planning periods in order to use in our problem. Then GAMS/Cplex 8.1 Solver is used to solve it. Finally, GAMS solution file gives the following data:
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
35
Fig. 1. Representation of the solution vector in SA implementation of JRSP.
Optimal number of worker (N): 14. Demands: D111 = 7 (c111=250), D121 = 7 (c121=200), D131 = 7 (c131=200), D141 = 7 (c141=150), D211 = 7 (c211=300), D221 = 7 (c221=240), D331 = 7 (c331=240), D241 = 7 (c241=180). In this example, the workload cost variables (cmik) was assigned as improvised numbers and given in parenthesis with demand. As seen from Fig. 2, the maximum workload is 920 for some workers. The aim of searching with heuristically is to find the assignments that prove workload balancing. In this example, the total workload is 12320 and the workload of each worker is desired as 880 (12320/14). The best solution of this example with heuristically is given in Appendix A.1. In Fig. 1, demand D241 is located in rotation 1 in worker 1, demand D121 is located in rotation 1 in worker 2 and so on. If we represent the solution vector formally as SOLns then SOL11 = D241 and SOL12 = D121, etc. The workload cost of each worker is calculated by summation of cost that related to demand. The total workload cost of a solution is the summation of all workers’ workload costs as shown in the right lower bottom of Fig. 1. In this work, solution quality of the problem is observed by standard deviation (S) of assigned workload costs for each worker. The best solution is obtained when standard deviation is zero. Neighbourhood move: A transition from one configuration to another that result in neighbour solution. Neighbour solutions are generated by a swapping procedure that is defined as follows: Step 1: Select the worker (row) who has the lowest workload cost randomly. Step 2. Select the worker (row) who has the highest workload cost randomly. Step 3. If workers have not off-days, determine two rotation part (column) that have same days randomly. If workers have off-days, determine two rotation that have same days and not shared days randomly. For example: Worker A: D214 D425 D145 D151 D234.
36
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
Fig. 2. Neighbourhood move. (i) Before neighbourhood move, (ii) after neighbourhood move.
Worker B: D433 D415 D373 D122 D253. Possible swaps between days: 1–1, 2–2, 3–3, 5–5, 4–7. An example of neighbourhood move for the above illustrative problem is shown in Fig. 2. In Fig. 2 cost values are filled randomly to show the effect of a neighbourhood move. Randomly selected workers and rotations are bolded. Present swapping strategy always generates feasible solutions and can be implemented easily. Objective function: A measure of how good a solution is. Standard deviation (S) of worker’s workload cost as shown in the right lower bottom of in Figs. 1 and 2 constitutes the objective function of the SA procedure. Acceptance/rejection of a solution: The neighbourhood solution is accepted if the objective function improves, otherwise the solution is accepted with probability depending on the temperature, which is set to allow the acceptance of a large proportion of generated solutions at the beginning (high temperature). Then, the temperature is modified (lowered) to decrease the probability of acceptance. At each temperature a predetermined number of moves are attempted. Assume that i and j are two solutions with standard deviation
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
37
f(i) and f(j), respectively. Then the acceptance criterion determines whether j is accepted from i by applying the following acceptance probability: 1; if f ðjÞ 6 f ðiÞ; P c faccept jg ¼ ððf ðjÞf ðiÞÞ=T Þ ; if f ðjÞ > f ðiÞ; e T 2 Rþ denotes control parameter ðtemperatureÞ:
ð7Þ
The purpose of sometimes accepting an inferior solution is to avoid becoming stuck in a local optimum. Termination: If the previously determined parameters are reached the algorithm terminates, e.g. maximum number of iterations, final temperature, no improvement for a number of iterations, etc. Cooling schedule: Selection of the SA parameters is called a cooling schedule. This gives how high the starting temperature should be, and rules to determine it. This also gives (a) when the current temperature should be lowered, (b) by how much the temperature should be lowered, and (c) when the annealing process should be terminated. 2.2.2. Determination of the cooling schedule The cooling schedule of the SA algorithm for the job rotation scheduling problem is determined as follows. Initial temperature (Tin): As mentioned above, the implementation of SA is based on the process of annealing of solids. The initial temperature must be set sufficiently high, resulting in a high degree of randomness in the initial stages of the search. Thus most of the initial movements are accepted. The value of initial temperature is regulated in such a way that at the beginning of SA the value of acceptance probability (Pc) is very closer 1. In this application it is set to 0.99. The value of acceptance probability (Pc) is determined by using factorial design. An estimation for the initial value of the temperature can be made by using the following equation [11]: fmin fmax T in ¼ ; ð8Þ ln P c where fmin and fmax are lower and higher bounds for a given JRSP, respectively. For this application, fmin = total workload/number of workers and fmax = 1.10 · (total workload/number of workers) are selected. Length of the Markov chain at iteration (L): It is chosen in such a way that for each temperature a thermal equilibrium can be achieved. But as the initial temperature approaches to zero, the acceptance probability of movements decreases. As a result of this, the length of Markov chain goes to infinity. Because of this situation, L is bounded by a constant which is dependent on the problem size. In this work, it is defined by L = (number of worker). Rate of cooling (a): The temperature should be lowered slowly in order to attain slow cooling. The temperature is decreased by using the following relation; Tin = aTn1, where a is the cooling rate its value being between zero and one. The initial temperature is estimated based on the initial acceptance probability, the final temperature (Tf) can also be estimated similarly based on the final acceptance probability (Pf). Therefore, cooling rate is also dependent on the values of initial and final acceptance probabilities and the number of temperature reduction cycles (t) (number of iterations). The definition of Bennage and Dhingra [16] for cooling rate is employed in this application: 1=ðt1Þ ln P c a¼ : ð9Þ ln P f Pf approaches to zero near to the end of the annealing process. It must be set very close to zero in order to avoid accepting non-improving solutions. It is set to 1 · 1015 in this application. Termination criterion: When the maximum number of iterations (t) is reached the SA algorithm terminates. t can be determined in such a way that at the end of the last iteration, Tf can be reached. In order to check this, the value of final temperature can be calculated by using Eq. (10). If the value of the final temperature calculated using Eq. (10) is not very close or equal to estimated final temperature then t must be reselected and the rate of cooling (a) should be recomputed. This procedure should be repeated until the required number of iterations (t) is determined. T f ¼ T in at :
ð10Þ
38
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
2.2.3. The pseudo code of the proposed SA algorithm for the JRSP problem The pseudo code of the proposed SA algorithm for solving JRSP is given as follows: GENERATE the feasible solution with X (generate the solution randomly) Calculate the objective function (f(X) = S); Set the initial temperature (Tin); Set the number of iteration (t); Set the cooling ratio (a); Set the initial acceptance probability (Pc); Set the final acceptance probability (Pf); Set the lower and higher bounds (fmin, fmax); t = 0; REPEAT l=0 REPEAT Generate a neighbouring solution as explained above and calculate the new objective function value. d = f(X 0 ) f(X); s 0 IF(d 6 0) THEN X = X ; ELSE produce u (random number (0, 1)); 0 IF (u < exp(d/T)) THEN X = X ; l = l + 1; UNTIL l = L; ‘‘number of workers’’; t = t + 1; Tin = aTn1 ; UNTIL ‘‘Termination criterion (t)’’; Best solution (X* the minimum standard deviation solution found), f(X*) and W; 3. Computational results 3.1. Scheduling environment The environment used in this research is based on a tour scheduling problem originally utilized by Topa¨ zkarahan [15]. In this prior study, labor requirements across a 12-h-day, 7 days a week operating log˘lu and O environment was to be satisfied full-time employees working five consecutive days per week. Shifts were eight hours in length with 1-h meal break between third and tenth hours of work. Shift start times are five for each day of work. In this formulation, at any time the number and duration of break periods can be changed. Additionally, the feasible set of shift types, the location of break windows, and the feasible work patterns for the allowable shift lengths can be easily be specified by the scheduler, and thus the model can be adapted to desired scheduling environment conveniently. The objective is to minimize the number of employees scheduled subject to the coverage of labor requirements. In this study, job rotation scheduling is performed based on tour scheduling environment. We exploit the tour scheduling environment to make our SA more meaningful. Thus we assumed that minimum number of employees (N) and labor requirements (Dmik) are known. For this application, we generated randomly a varying demand pattern for each work of type within the hours of a day and across the days of a week. The generated demand patterns were given in Appendix A.2. 3.2. Heuristic and computational time performance We describe our computational experiments about JRSP. We generated the test problems based on combinations of types of work. This gives us 26 different situations (see Table 2). This will prove us an idea of the effectiveness of the algorithms with different dimensions of rotations.
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
39
In order to test the IP model, one run was performed for each situation. In order to test the SA, five replications were performed for each situation. A limitation of integer programming is that it is designed to optimize deterministic problems or problems where the data is known with certainty. Unfortunately, the test problems described in this section does not have this characteristic; because of real world situation were not considered. One way to handle this difficulty is to use improvised cmik values that change according to work types and days for test problems. This permits the development of the necessary coefficients to construct both the IP and SA. The improvised workload cost of each type of work and day is given in Table 1. Furthermore, it is easy to calculate total workload cost for each problem when numbers of worker and labor requirements across a 12-h-day are determined. Additionally, Table 2 gives numbers of worker and the total workload cost for each problem. GAMS/Cplex Solver 8.1 has been used to solve the test problems on a Pentium IV 2.6 Ghz computer. All test problems has been written in GAMS software and solve CPLEX solver. The search for the optimal solution has been limited to 86,400 CPU s. IP formulation has not found optimal solutions for all test problems. At the end of time limit, a near optimal solution was obtained. This was used for comparing with SA algorithm. For SA algorithm, an initial experiment was done to test the different parameters. We ran tests with iteration number of 1000, 1500, and 2000. We also performed experimental runs with length of the Markov chain at iteration (number of workers) of L/2, L, and 2L. Finally we tested with acceptance probability of 0.90, 0.95, and 0.99. Based on the results, we found that the following parameter values gave the best results for SA algorithm: iteration number = 1000, length of the Markov chain at iteration = L and acceptance probability = 0.99. The IP model and SA results are shown in Table 3. In the IP results, it is not easy to calculate the standard deviation of a schedule. Because of this reason, the results are compared according to objective function value (near optimal solution) for each test problems. Table 3 gives the minimum workload cost solution for each of the 26 problems using the following methods:
Table 1 The improvised workload costs Work of type
Monday
Tuesday
Wednesday
Thursday
Friday
Saturday
Sunday
1 2 3 4 5
250 300 120 47 18
200 240 120 47 18
200 240 80 47 18
150 180 80 35 12
150 120 60 24 12
100 120 40 12 6
100 90 40 12 6
Table 2 Descriptions of the test problems Problem no.
Type of works (m)
Number of worker
Total workload
Problem no.
Type of works (m)
Number of worker
Total workload
Problem no.
Type of works (m)
Number of worker
Total workload
1 2 3 4 5 6 7 8 9
1–2 1–3 1–4 1–5 2–3 2–4 2–5 3–4 3–5
39 36 36 31 23 23 18 20 15
35,000 26,210 23,805 22,424 17,240 14,835 13,260 5855 4496
10 11 12 13 14 15 16 17 18
4–5 1–2–3 1–2–4 1–2–5 1–3–4 1–3–5 1–4–5 2–3–4 2–3–5
15 49 49 44 46 41 41 33 28
2091 38,720 34,860 35,416 27,865 27,439 24,191 18,995 17,598
19 20 21 22 23 24 25 26
2–4–5 3–4–5 1–2–3–4 1–2–3–5 1–2–4–5 1–3–4–5 2–3–4–5 1–2–3–4–5
28 25 59 54 54 51 38 64
15,153 6239 40,945 39,570 37,165 28,139 19,331 41,311
40
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
Table 3 Results of the computational study with 26 test problems Problem name
(1) Rotation between job 1 and 2 (2) Rotation between job 1 and 3 (3) Rotation between job 1 and 4 (4) Rotation between job 1 and 5 (5) Rotation between job 2 and 3 (6) Rotation between job 2 and 4 (7) Rotation between job 2 and 5 (8) Rotation between job 3 and 4 (9) Rotation between job 3 and 5 (10) Rotation between job 4 and 5 (11) Rotation between job 1, 2 and 3 (12) Rotation between job 1, 2, and 4 (13) Rotation between job 1, 2, and 5 (14) Rotation between job 1, 3, and 4 (15) Rotation between job 1, 3, and 5 (16) Rotation between job 1, 4, and 5 (17) Rotation between job 2, 3 and 4 (18) Rotation between job 2, 3 and 5 (19) Rotation between job 2, 4 and 5 (20) Rotation between job 3, 4 and 5 (21) Rotation between job 1, 2, 3 and 4 (22) Rotation between job 1, 2, 3 and 5 (23) Rotation between job 1, 2, 4 and 5 (24) Rotation between job 1, 3, 4 and 5 (25) Rotation between job 2, 3, 4 and 5 (26) Rotation between job 1, 2, 3, 4 and 5
IP Computational time (CPU)
SA Computational time (CPU)
IP Best solution value (W)
SA Best solution value (W) Run 1
Run 2
Run 3
Run 4
Run 5
>86,400
238.36
910
980
950
950
950
1000
>86,400
110.02
730
760
770
780
780
800
>86,400
114.48
682
785
747
762
762
774
>86,400
81.67
750
850
850
850
850
812
>86,400
40.39
760
820
860
800
840
860
>86,400
40.98
684
707
719
797
742
742
>86,400
24.34
762
786
798
930
822
912
>86,400
30.78
314
322
309
312
311
311
>86,400
16.00
318
372
340
326
340
318
>86,400
17.66
153
153
159
153
153
159
>86,400
232.57
840
840
870
890
870
880
>86,400
226.95
822
910
910
910
900
910
>86,400
182.56
870
900
850
920
900
940
>86,400
206.75
647
642
644
724
660
674
>86,400
153.47
700
750
750
750
750
750
>86,400
152.57
621
662
668
685
718
647
>86,400
91.10
607
647
680
622
694
667
>86,400
62.54
666
690
930
728
718
678
>86,400
60.00
551
612
629
646
664
588
>86,400
49.69
274
307
266
334
289
286
>86,400
373.09
752
757
747
754
794
750
>86,400
373.20
806
856
840
820
840
840
>86,400
379.28
736
787
757
802
820
737
>86,400
312.69
606
618
598
658
700
614
>86,400
120.93
542
631
667
611
623
553
>86,400
474.89
690
736
672
764
780
742
IP: Integer programming model (searching for the optimal solution). SA: Simulated annealing algorithm (searching for the best heuristic solution).
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
41
Table 4 Solution quality of SA according to IP Problem no.
SA performance according to IP
Average
The solution quality¼ 1 SAIP IP 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
Run 1
Run 2
Run 3
Run 4
Run 5
0.92 0.96 0.85 0.87 0.92 0.97 0.97 0.97 0.83 1.00 1.00 0.89 0.97 1.01 0.93 0.93 0.93 0.96 0.89 0.88 0.99 0.94 0.93 0.98 0.84 0.93
0.96 0.95 0.90 0.87 0.87 0.95 0.95 1.02 0.93 0.96 0.96 0.89 1.02 1.00 0.93 0.92 0.88 0.60 0.86 1.03 1.01 0.96 0.97 1.01 0.77 1.03
0.96 0.93 0.88 0.87 0.95 0.83 0.78 1.01 0.97 1.00 0.94 0.89 0.94 0.88 0.93 0.90 0.98 0.91 0.83 0.78 1.00 0.98 0.91 0.91 0.87 0.89
0.96 0.93 0.88 0.87 0.89 0.92 0.92 1.01 0.93 1.00 0.96 0.91 0.97 0.98 0.93 0.84 0.86 0.92 0.79 0.95 0.94 0.96 0.89 0.84 0.85 0.87
0.96 0.93 0.88 0.87 0.89 0.92 0.92 1.01 0.93 1.00 0.96 0.91 0.97 0.98 0.93 0.84 0.86 0.92 0.79 0.95 0.94 0.96 0.89 0.84 0.85 0.87 Average
0.95 0.94 0.88 0.87 0.91 0.92 0.91 1.00 0.92 0.99 0.97 0.90 0.97 0.97 0.93 0.89 0.90 0.86 0.83 0.92 0.98 0.96 0.92 0.92 0.84 0.92 0.92
Table 4 gives the solution quality for each problem using the equation that is given in this table. The best performance of SA is highlighted. SA performs better than IP for some problems (8, 13, 14, 20, 21, 24 and 26) and this is meaning that the IP solution is not an optimal. The average performance of SA was found 0.92. This is quite effective. If there is a comparison opportunity in point of minimization of standard deviation with other heuristics in the literature, the performance of SA algorithm could be found better than currently. The SA algorithm is programmed in Visual C++. It is capable of solving JRSP once the required data files are generated and the SA parameters. The data file of problem 2 is given in Appendix A.3. The SA algorithm could find the solution within 10 min for all test problems. The maximum solution time is 474.89 CPU times for problem 26. Average solution time is 160.267. 4. Conclusions In this paper, we examined the use of simulated annealing algorithm in solving the job rotation scheduling problem. SA has been applied to numerous optimization problems successfully. Firstly in this study its application to JRSP is presented. In this study, a SA algorithm with a simple but effective data structure and neighbourhood generation mechanism is proposed for its solution. The proposed SA algorithm is applied to test problems from generated randomly. In these tests, the proposed algorithm is compared with IP model results. Unfortunately, the optimal solution can not be obtained with IP model at end of CPU limit. In order to compare optimal and heuristic results, IP results were admitted the best integer solutions within time limit. The performance of SA was evaluated based on IP solutions for each problem. The computational results showed that the proposed algorithm is quite effective.
42
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
Appendix A A.1. The best solution for example within Section 2.2
1 2 3 4 5 6 7 8 9 10 11 12 13 14
1
2
3
4
Workload cost
D111 D111 D211 D111 D111 D211 D211 D211 D111 D211 D211 D211 D111 D111
D221 D221 D121 D221 D221 D121 D121 D121 D221 D121 D121 D121 D221 D221
D231 D231 D131 D231 D231 D131 D131 D131 D231 D131 D131 D131 D231 D231
D141 D141 D241 D141 D141 D241 D241 D241 D141 D241 D241 D241 D141 D141
880 880 880 880 880 880 880 880 880 880 880 880 880 880 s=0 W = 12320
Total
A.2. Demands for test problems Work of type 1 1 2 1.1 3 6 1.2 3 7 1.3 2 7 1.4 3 7 1.5 4 8 1.6 2 5 1.7 1 3
3 8 8 8 8 8 7 5
4 10 9 9 10 10 8 7
5 14 13 12 12 13 9 8
6 14 13 12 12 13 9 8
7 9 9 9 9 9 9 8
8 9 9 9 9 9 8 8
9 13 11 10 10 12 8 8
10 13 11 10 10 12 8 7
11 13 11 10 10 10 8 8
12 9 9 9 9 9 8 8
Work of type 2 1 2 2.1 3 3 2.2 2 6 2.3 2 5 2.4 3 5 2.5 3 5 2.6 2 2 2.7 1 2
3 4 5 3 5 4 1 2
4 5 6 4 3 6 3 3
5 9 9 9 8 4 4 3
6 8 9 9 5 6 3 2
7 6 9 9 5 9 3 2
8 5 6 9 10 10 3 2
9 7 6 8 6 4 4 2
10 6 7 6 4 4 1 1
11 4 2 5 2 2 1 2
12 2 1 4 1 2 1 1
Work of type 3 1 2 3.1 1 2 3.2 1 2
3 3 3
4 5 4
5 6 6
6 6 6
7 4 4
8 4 4
9 5 5
10 5 4
11 5 4
12 4 4
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
43
(continued) Work of type 3 3.3 1 2 3.4 1 2 3.5 2 2 3.6 1 1 3.7 1 1
2 2 3 2 2
4 4 4 3 2
6 6 6 3 2
6 5 6 3 2
4 4 4 3 3
4 4 4 2 2
5 5 6 2 2
4 4 6 2 2
4 4 5 2 2
4 4 4 2 2
Work of type 4 1 2 4.1 1 2 4.2 1 2 4.3 1 2 4.4 1 2 4.5 2 2 4.6 1 1 4.7 1 1
3 3 3 2 2 3 2 2
4 5 4 4 4 4 3 2
5 6 6 6 6 6 3 2
6 6 6 6 5 6 3 2
7 4 4 4 4 4 3 3
8 4 4 4 4 4 2 2
9 5 5 5 5 6 2 2
10 5 4 4 4 6 2 2
11 5 4 4 4 5 2 2
12 4 4 4 4 4 2 2
Work of type 5 1 2 5.1 1 1 5.2 1 1 5.3 1 1 5.4 1 2 5.5 1 1 5.6 1 1 5.7 1 1
3 1 1 1 2 1 1 1
4 2 2 1 2 1 1 1
5 3 3 2 2 2 2 1
6 3 3 3 2 2 2 1
7 2 2 2 2 2 1 1
8 2 2 2 2 2 1 1
9 2 2 2 2 1 1 1
10 2 2 1 2 1 1 1
11 2 1 1 1 1 1 1
12 1 1 1 1 1 1 1
A.3. Example data file for problem 2 Dmik
cmik
n
111 113 114 115 121 122 123 124 125 131 133 134 135 141 144 145 151
250 250 250 250 200 200 200 200 200 200 200 200 200 150 150 150 150
7 1 5 9 3 5 1 2 9 7 1 1 9 8 2 9 8 (continued on next page)
44
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
(continued) Dmik
cmik
n
154 155 161 164 165 171 173 174 175 311 313 314 315 321 323 325 331 334 335 341 344 345 351 353 354 355 361 364 365 371 373 375
150 150 100 100 100 100 100 100 100 120 120 120 120 120 120 120 80 80 80 80 80 80 60 60 60 60 40 40 40 40 40 40
3 9 7 1 8 3 3 1 8 2 2 1 4 2 2 4 2 2 4 2 2 4 2 1 1 4 2 1 2 1 1 2
References [1] B.J. Carnahan, M.S. Redfern, B. Norman, Designing safe job rotation schedules using optimization and heuristic search, Ergonomics 43 (4) (2000) 543–560. [2] J.M. Allwood, W.L. Lee, The impact of job rotation on problem solving skills, International Journal of Production Research 42 (5) (2004) 865–881. [3] P. Paul, F.M. Kuijer, K. Visser Bart, C.G. Han, Job rotation as a factor in reducing physical workload at a refuse collecting department, Ergonomics 42 (9) (1999) 1167–1178. [4] D. Ghare Gaurav, T. Leutenegger Scott, Improving small job response time for opportunistic scheduling, in: IEEE International Workshop on Modeling, Analysis, and Simulation of Computer and Telecommunication Systems – Proceedings, 2000, pp. 557–564. [5] A.B. Becker, N.G. Swanson, S.L. Sauter, T.L. Galinsky, S. Jones, L. Schleifer, Compatibility of job rotation subtasks in data entry work, in: Proceedings of the Human Factors and Ergonomics Society, vol. 2, 1995, p. 965. [6] R. Pickavance, Employee utilisation and deployment (Inst. of Manpower Studies.) report: IMS-115, ISBN-1-85184-013-3, 1985, p. 39. [7] R.T. Wilkinson, R.S. Edwards, Stable hours and varied work as aids to efficiency (Royal Naval Personnel Research Committee London (England)) report: OES-512, October 1969, p. 9. [8] B.J. Carnahan, M.S. Redfern, B.A. Norman, A genetic algorithm for designing job rotation schedules considering ergonomic constraints, in: Proceedings of the 1999 Congress on Evolutionary Computation-CEC99 (Cat. No. 99TH8406), vol. 2, 1999, pp. 1093– 1098.
S.U. Sec¸kiner, M. Kurt / Applied Mathematics and Computation 188 (2007) 31–45
45
[9] B.J. Carnahan, M.S. Redfern, Automatic planning of job rotation schedules, in: Industrial Engineering Research – Conference Proceedings, 1997, pp. 263–268. [10] M. Brusco, L. Jacobs, A simulated annealing approach to the solution of flexible labor scheduling problems, Journal of the Operational Research Society 44 (12) (1993) 1191–1200. [11] A. Baykasog˘lu, N.N.Z. Gindy, A simulated annealing algorithm for dynamic layout problem, Computers and Operations Research 28 (2001) 1403–1426. [12] C.R. Reeves, Modern Heuristic Techniques for Combinatorial Problems, John Wiley& Sons Inc, Newyork, Toronto, 1993, pp. 288– 298. [13] S. Kirkpatrick, C.D. Gelatt Jr., M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680. [14] D. Abramson, Constructing school timetables using simulated annealing: sequential and parallel algorithms, Management Science 37 (1991) 98–113. _ O ¨ zkarahan, Implicit optimal tour scheduling flexible assignments, Computers and Industrial Engineering 44 (2003) [15] S ß . Topalog˘lu, I. 75–89. [16] W.A. Bennage, A.K. Dhingra, Single and multiobjective structural optimization in discrete–continuous variables using simulated annealing, International Journal for Numerical Methods in Engineering 38 (1995) 2753–2773.