Hybrid Harmony Search for Stochastic Scheduling of Chemotherapy Outpatient Appointments

Abstract

This research deals with the same-day chemotherapy outpatient scheduling problem that is recognized as a leading strategy to pursue the objective of reducing patient waiting time. Inspired by a real-world context and different from the other studies, we modeled a multi-stage chemotherapy ward in which the pharmacy is located away from the treatment area and drugs are delivered in batches. Processes in oncology wards are characterized by several sources of uncertainty that increase the complexity of the problem; thus, a stochastic approach was preferred to study the outpatient scheduling problem. To generate effective appointment schedules, we moved in two directions. First, we adopted a late-start scheduling strategy to reduce the idle times within and among the different stages, namely medical consultation, drug preparation and infusion. Then, since the problem is NP-hard in the strong sense, we developed a hybrid harmony search metaheuristic whose effectiveness was proved through an extended numerical analysis involving another optimization technique from the relevant literature. The outcomes from the numerical experiments confirmed the efficacy of the proposed scheduling model and the hybrid metaheuristic algorithm as well.

1. Introduction

The Chemotherapy Outpatient Scheduling (COS) problem is a sensitive and challenging issue that is capturing the attention of both scholars and stakeholders in the healthcare landscape. Building an effective schedule of patients allows for reducing their waiting times on the one hand and increasing the number of treatments in the working shift on the other hand. To strengthen this thought, the health wards may exploit an appointment scheduling approach to improve the efficiency of their services, thus increasing patient satisfaction and service rate as well [1].

In brief, the main steps to generate a daily schedule of appointments in an oncology clinic are the following. The oncologists plan the days of treatment for each patient through a specific medium-term care protocol that defines all the necessary information, e.g., date and duration of treatments, type and doses of drugs. Subsequently, the oncology unit communicates to the patients the appointment time for a given day. Since each patient may have a different disease history, thus requiring a different care path, the time he/she needs to receive the oncology treatment, which goes from the medical visit to the end of treatment, can be highly variable and the interaction with the other patients can strongly bias his/her length of stay in the clinic. Hence, an effective method for scheduling patient appointments may positively impact the performance of the clinic and patient satisfaction as well.

The literature background reveals that the COS problem can be viewed from two different perspectives, namely offline or online scheduling [2,3]. As for offline scheduling, the daily list of patients is known in advance and the scheduler usually generates the appointment schedule a few days before the treatments. In this case, a non-templated approach is usually adopted to optimize an objective function based on a priori information [2]. Instead, for online scheduling, the appointment time is communicated to the patient immediately after his request or shortly after (e.g., no later than 24 h) by employing a dynamic template. In general, a template is an appointment timetable in which the total number of slots, the length of each slot, and the characteristics of patients to be scheduled for each slot are specified [4]. In fact, in online scheduling, a series of vacant appointment slots of the designed template, preliminarily generated based on a specific criterion, have to be filled by the patient requests.

Another classification criterion concerns the chemotherapy treatment and depends on when the patients undergo the medical consultation. In next-day chemotherapy, the medical consultations are processed the day before the treatment. The scheduling problem consists of allocating patients to chairs, assigning them a nurse and defining the related treatment starting time. This strategy enables pharmacists to prepare the therapy in advance. However, the patient has to go to the clinic on two consecutive days and that could be not so practical if they live far away from the clinic. In the same-day chemotherapy strategy, the medical consultations are carried out on the same day of the infusion therapy. This strategy is quite well-diffused in real-world oncology clinics as it is preferred by patients [5] and prevents any waste of expensive drugs in case of deferral of patients or no-show at the infusion stage [6].

Several solving approaches were employed in the literature to cope with the COS problem. The contribution by Dobish [7] can be considered seminal research on the problem under investigation. The author adopted a next-day scheduling method based on a timetable for each day of the week. Similarly, Edwards et al. [8] made use of a specific template or ad hoc rules, based on the patient acuity level. From then on, most literature recognized the need to explicitly consider the patient flow to improve the performance of oncology units for solving the COS problem.

In this regard, Sadki et al. [9] studied the same-day appointment scheduling problem with two major resources, namely oncologists for consultation and beds for injection, and proposed a combination of heuristics based on the Lagrangian relaxation to minimize the patient waiting times and the makespan. Turkcan et al. [10] used mathematical programming models to hierarchically tackle both the planning and scheduling of chemotherapy patients to minimize the patient waiting time while considering the limited availability of nurses. Later, Heshmat et al. [11] improved the model proposed by Turkcan et al. [10] so that it could be solved in a smaller computational time and also for larger-sized instances. Condotta and Shakhlevich [12] exploited mathematical programming to generate a multilevel template, which aims to minimize the patient waiting time and nurse workload in an outpatient clinic.

Discrete event simulation and mathematical programming were used by Liang et al. [13] to compare the proposed appointment scheduling tool with the current practice under a series of operational measures such as patient waiting time, clinic total working time and resource utilization. Notably, they solved a linear programming model to schedule patient appointments according to a same-day offline strategy. Liang and Turkcan [14] distinguished between functional and primary care delivery models to provide chemotherapy treatments to cancer patients, depending on the availability of nurses. They consider a single-stage system where a set of patients have to be scheduled offline on the same day only for the infusion phase and proposed two multi-objective optimization models based on mathematical programming that ignored both oncology consultation and lab tests. Bouras et al. [15] introduced a mixed integer programming (MIP) model for reducing the patient waiting time of a same-day offline COS problem in a real-life oncology unit. They modeled the whole set of system stages, namely oncology consultation, drug preparation and injection, also considering the limited number of resources at the different stages. To reduce the number of binary variables and constraints, thus enhancing the computational efficiency of their COS approach, Heshmat et al. [16] devised a two-stage COS method, properly inspired by cellular manufacturing systems, which involved a clustering phase and a mathematical programming phase for the minimization of the total completion time referred solely to the injection stage of an oncology unit.

Another valuable contribution to the COS scenario is attributable to Hesaraki et al. [2], who focused on the infusion stage to generate an online scheduling method subject to the nursing constraint. They used integer programming to design a template of vacant appointment slots that follows a bicriteria objective based on the combination of weighted flow time and makespan. A different perspective emerges from the research by Huggins and Claudio [17]. They presented a mathematical model for the next-day COS problem that managed the chemotherapy patient appointments while taking into consideration the workload of nurses and pharmacists as a constraint of the optimization problem in a cancer clinic. Constraint programming was applied by Hahn-Goldberg et al. [18] to the next-day outpatient scheduling problem considering the online approach. A similar approach for the same-day case can be attributed to Huang et al. [19], who developed a chemotherapy outpatient scheduling template by reducing the violation between resource assignment and treatment requirements.

Recently, a research stream emphasized the need for managing uncertainty in healthcare scheduling problems [20,21]. Several sources of uncertainty may have a notable impact on the chemotherapy path of oncology outpatients, therefore deterministic models could represent a strong simplification [22]. Castaing et al. [23] presented a two-stage stochastic programming model for the next-day COS problem of an oncology unit located in the USA. Since this optimization method requires a prohibitive computational time to be solved, the authors introduced a heuristic algorithm to find approximate solutions in a reasonable time. Alvarado and Ntaimo [22] used a mean-risk stochastic programming model powered by a specific heuristic to schedule patient appointments and resources for reducing patient waiting time and nurse overtime. Mandelbaum et al. [24] exploited the principles of queueing theory for the offline appointment sequencing problem by engaging a large number of servers (chairs) and customers (patients) in a stochastic system wherein service duration and punctuality are subject to significant uncertainty. They proposed a data-driven approach based on the infinite-server queues whose effectiveness was proved by testing that approach against near-optimal algorithms.

Both in deterministic and in stochastic mathematical programming models, the time to converge drastically depends on the number of patients and the number of resources as well. Hence, heuristic or metaheuristic algorithms represent a valid alternative to achieve a perfect compromise between the quality of solutions and computational times in solving COS problems. To this end, Sevinc et al. [3] proposed a two-phase approach for the next-day COS problem. They used a specifically devised heuristic algorithm for handling the appointment scheduling for the laboratory tests, and two heuristics based on the multiple knapsack problem for the second phase in which patients have to be online assigned to the infusion seats. The work of Garaix et al. [25] represents a valuable contribution since a metaheuristic algorithm was adopted in the field of the stochastic COS problem. They developed a GRASP algorithm to generate sub-optimal lists of patients for consultation and injection phases in a same-day chemotherapy treatment scenario. Finally, Demir et al. [26] developed a heuristic method based on a progressive hedging algorithm, since they experienced high computational times to solve instances with a two-stage stochastic programming model focused only on the treatment stage.

The present paper, inspired by the health services provided by an oncology unit of a hospital located in Southern Italy, aims to solve the same-day offline stochastic COS problem, motivated by the activities of a real-life oncology department in which the medical consultations are performed on the same day as the chemotherapy treatments. Different from the literature, this research is the first to consider that the pharmacy collects the therapies in batches and a courier service delivers the batches to the oncology department. Our objective is to make sure that the time each patient spends in the system, i.e., the total flow time, is minimized. The aim of the proposed scheduling method is to enhance the daily experience of each patient by reducing their waiting time. Furthermore, the proposed research aims at reducing resource idle time and maximizing the number of patients scheduled in a day.

The main contributions of the present research are the following:

• Adopting a stochastic scheduling strategy that is able to reduce the total flow time (i.e., idle times of patients among the stages) in a same-day COS problem;
• Introducing a novel multi-stage simulation framework (based on discrete-time recursive equations [27,28]) in which the pharmacy is detached from the oncology department and drugs are delivered in batches;
• Developing a new hybrid metaheuristic algorithm, namely Harmony Search (HS), to solve the COS problem for the minimization of patient waiting time. The proposed algorithm is ‘hybrid’ since heuristic solutions in the initial population, the local search and the reinitialization procedure are embedded in the metaheuristic algorithm;
• Differently from the existing scheduling approaches (e.g., templated scheduling), the proposed method allows the healthcare manager to identify the optimal schedule in the space of the solution with reasonable computational time.

Table 1. Comparison with the literature.

Reference	Medical Consultation	Scheduling Strategy	Type of Model	Solving Method	KPI	Limited Resources			Drug Transp.
						ON	PH	NU
Alvarado and Ntaimo (2018) [22]	Next-day	Offline	Stochastic	SP	WT, WL			✓
Bouras et al. (2017) [15]	Same-day	Offline	Deterministic	MP	WT	✓	✓	✓
Castaing et al. (2016) [23]	Next-day	Offline	Stochastic	SP, HE	WT, TCT			✓
Condotta and Shakhlevich (2014) [12]	Next-day	Offline	Deterministic	MP, T	WT, WL			✓
Demir et al. (2021) [10]	Same-day	Offline	Stochastic	SP, HE	WT, CU, OT			✓
Dobish (2003) [7]	Next-day	Offline	Deterministic	T	WT, WL			✓
Edwards et al. (2017) [8]	Next-day	Offline	Deterministic	T	CU, P			✓
Garaix et al. (2020) [25]	Same-day	Offline	Stochastic	SP, GRASP	M	✓		✓
Hahn-Goldberg et al. (2014) [18]	Next-day	Online	Deterministic	CP, DT	M		✓	✓
Hesaraki et al. (2019) [2]	Same-day	Online	Deterministic	MP, DT	FW, M			✓
Heshmat et al. (2017) [11]	Next-day	Offline	Deterministic	MP, CL	TCT			✓
Heshmat et al. (2018) [16]	Next-day	Offline	Deterministic	MP, CL	TCT			✓
Huang et al. (2019) [19]	Same-day	Offline	Deterministic	CP, T	WL			✓
Huggins and Claudio (2019) [17]	Next-day	Offline	Deterministic	MP	P		✓	✓
Liang and Turkcan (2016) [14]	Same-day	Offline	Deterministic	MP	WT, WL, OT			✓
Liang et al. (2015) [13]	Same-day	Offline	Deterministic	MP	WL	✓		✓
Mandelbaum et al. (2020) [24]	Same-day	Offline	Stochastic	IR	WT, OT			✓
Sadki et al. (2011) [9]	Same-day	Offline	Deterministic	LR	WT, M	✓		✓
Sevinc et al. (2013) [3]	Next-day	Online	Deterministic	HE	CU			✓
Turkcan et al. (2012) [10]	Same-day	Offline	Deterministic	MP	TCT			✓
This study	Same-day	Offline	Stochastic	HS	F	✓	✓	✓	✓

Legend: CL: Clustering; CP: Constraint Programming; CU: Chair Utilization; DT: Dynamic Template; F: Total Flow time; FW: Weighted Flow time; GRASP: Greedy Randomized Adaptive Search Procedure; HE: Heuristics; HS: Hybrid Harmony Search; IR: Infinite-server Relaxation; LR: Lagrangian Relaxation; M: Makespan; MP: Mathematical Programming; OT: Overtime of resources; P: Number of Patients; SP: Stochastic Programming; T: Template; TCT: Total Completion Time; WL: Work Load; WT: Waiting Time.

aims to retrieve the main literary contributions for the optimization of the COS problem and adopts a series of classification criteria to stress the difference between the literature and the proposed study as well. The first two classification criteria refer to the way the medical consultation is run and to the scheduling strategy, respectively. A distinction between deterministic and stochastic approaches (type of model) is also indicated along with the adopted solving method. In each work, it can be observed that the solving methods were compared with the scheduling rule implemented in practice by the oncology unit. The objective function is defined as a Key Performance Indicator (KPI). The kind of resources explicitly involved in the model can be the following: oncologists for the consultation phase (ON), pharmacists for the drug preparation stage (PH) and nurses at the treatment stage (NU). Finally, the last column points out the models capable of running the drug transportation time and the transportation batch as well.

Notably, HS was compared with a specific Greedy Randomized Adaptive Search Procedure (GRASP) algorithm, since this is the unique metaheuristic algorithm already used in the literature for the COS problem. A benchmark of instances was generated and an extended comparison analysis was carried out to demonstrate the outperformance of the proposed optimization technique over the competitor. The rest of the paper is organized as follows. Section 2 describes the problem statement also including the basic assumptions and the objective function. The major components characterizing the proposed hybrid harmony search are dealt with in Section 3. Section 4 presents the numerical experiments and the output from the comparative analysis. A final discussion is reported in Section 5.

2. Problem Statement

The oncology department under consideration can be considered as the counterpart of a three-stage hybrid flow shop manufacturing system in which the first stage is related to the medical consultation, the second stage consists of the pharmacy laboratory and, finally, the third stage involves a set of chairs in parallel for the chemotherapy treatment [15,18].

The oncology unit disposes of O oncologists (o = 1, …, O), N nurses (n = 1, …, N) and C identical chairs (c = 1, …, C) for the treatment. The oncology unit collaborates with the pharmacy laboratory, whose capacity depends on D pharmacist technicians (d = 1, …, D) in charge of the drug preparation. Every day, a number of P patients (p = 1, …, P) need to be scheduled to receive oncology therapy. A referee oncologist assists the patient during the whole therapeutic path, in accordance with the care protocol [29]. Hence, the long/medium-term appointment-planning phase is managed by the referee oncologist, who decides upon the days any patient undergoes the therapy. For that reason, the set of patients to be treated every day is known a priori.

The sequence of steps that a patient encounters during the treatment day are the following. Once the patient arrives at the oncology department, he/she is registered in a welcome room. The patient waits until the referee oncologist becomes available. If the oncologist is available, the medical consultation starts and a decision on the chemotherapy is taken based on the health status of the patient and the blood tests executed the day before by the patient in the same hospital or an external laboratory. Therefore, there exists a deferral probability [25,30] that the patient p is not ready to receive the chemotherapy on the same day. In case the medical consultation is successful, the referee oncologist sends a request (i.e., a prescription) to the pharmacy, which includes the type and dose of the drug [31]. Once the prescription for the therapy arrives at the pharmacy, the pharmacist technicians can start the drug preparation. The pharmacy gathers the ready therapies in batches and makes use of a courier service to send the batches to the oncology department. Finally, once the drug arrives at the oncology center and both a nurse and a chair are available, the patient starts the treatment. After the therapy infusion is completed, the patient is discharged from the oncology ward.

Additionally, to thoroughly describe the problem under investigation regarding the ward we observed, the following assumptions, which are common to several literary contributions [2,25,26,32], can be summarized as follows:

• The number of patients to be treated on a given day is known in advance as it arises from the planning phase, conforming to the patient care protocol decided by the oncologists;
• Each patient is assigned to a referee oncologist for the medical consultation;
• Cancellations and no-shows at the consultation session are disregarded;
• The desk registration time is negligible;
• The results of blood exams are immediately available for medical consultation since the exams are executed by the patient in the previous days;
• Drug preparation devices never break down;
• Each pharmacist technician can prepare a therapy at a time;
• Pre-emption on the different activities is not allowed;
• Nurses have identical skills, i.e., each patient can be treated by any available nurse;
• Each nurse can prepare only one patient at a time;
• A nurse can monitor four patients simultaneously at most;
• The time any patient needs to leave the chair is negligible.

To model the uncertainties of the problem, the Sample Average Approximation framework was adopted [33]. This consists of generating a finite number of scenarios Ω in which stochastic parameters are independently sampled from the corresponding stochastic distributions, while deterministic parameters are kept unchanged for each scenario ω (ω = 1, ..., Ω). Six distinct sources of uncertainty characterize the problem under investigation: (i) the medical consultation time ${Dc}_p^{\omega }$; (ii) the drug preparation time ${Dt}_p^{\omega }$; (iii) the batch transportation time ${TD}_b^{\omega }$; (iv) the set-up time ${Ds}_p^{\omega }$; (v) the treatment time ${Di}_p^{\omega }$; (vi) the deferral probability ${\lambda }_p^{\omega }$, i.e., a patient p at scenario ω is deferred. In particular, ${\lambda }_p^{\omega }$ is derived from $U[0,1]$. If ${\lambda }_p^{\omega }$ is lower than or equal to $\overline{\lambda }$, that is the experimental value arising from the experimental observation, the patient is deferred, otherwise, she/he can undergo the treatment:

$${\delta }_p^{\omega }=\left\{\begin{array}{c}1\\ 0\end{array} \begin{array}{c}if \lambda { }_p^{\omega } \le \overline{\lambda } \\ otherwise\end{array}\right.$$

(1)

The objective is to minimize the patient waiting time over the provided time horizon, which consists of one working day. The waiting time is the time in which the patient has to wait for:

• The referee oncologist being available for the medical consultation;
• The therapy being prepared by the pharmacy and delivered to the oncology unit;
• A chair being available for the chemotherapy treatment;
• A nurse being available to start the set-up operations and monitor the therapy infusion.

In order to minimize the patient waiting time, the total flow time hereinafter denoted as F, is used as the objective function. F is the total time spent by the patient in the oncology unit. It consists of the sum of the overall waiting time and the time needed for the medical consultation and treatment. Since these latter intervals of time depend on stochastic distributions of input data and then remain unchanged during the simulation, the minimization of the total flow time F directly implies the reduction of the total patient waiting time [2,34,35]. F can be calculated as follows:

$$F ={\displaystyle \sum }_{p=1}^P{(C}_p-{Sc}_p)$$

(2)

where $C_p$ is the treatment completion time of the patient p and ${Sc}_p$ is the consultation starting time. The stochastic counterpart of the total flow time is the expected total flow time E(F), which depends on the set of scenarios Ω (ω = 1 ,..., Ω), as follows:

$$E\left(F\right) =\frac{1}{\Omega }{\displaystyle \sum }_{\omega =1}^{\Omega }F\left(\omega \right)$$

(3)

3. The Hybrid Harmony Search

The COS problem can be considered a hybrid flow shop problem, which is NP-hard in the strong sense, even when there are two resources at the first stage [36]. As a result, only very small-sized instances can be optimally solved with reasonable computational time. Moreover, due to the stochastic configuration of the problem at hand, the computational time required to solve an instance dramatically increases with the number of scenarios [10,37]. As a result, either heuristic or metaheuristic algorithms are needed to solve the problem in a reasonable computational time [38].

A hybrid Harmony Search (HS) algorithm, which belongs to the class metaheuristic algorithms motivated by music harmony, is proposed here. HS is an effective and efficient evolutionary technique able to solve different kinds of engineering problems [39], which showed better performance than other well-known optimization methods [40]. Several research contributions also demonstrated that the evolutionary mechanism of HS is faster than genetic algorithms [41,42].

A harmony solution consists of an n-dimensional real-coded vector. Let us suppose a single harmony is denoted as $x{=(x}_1,\dots {,x}_j,\dots {,x}_n)$ such that each variable is defined in the domain [$L{B}_j{,UB}_j$] $\in$ ℝ and f(x) is the related objective function value. Since the healthcare scheduling problem under investigation can be classified as a combinatorial issue, each real-coded solution has to be converted into a sequence of patients. To this end, we employed a well-known mechanism based on the smallest position value (SPV) rule [43], which allows for converting any real-valued harmony vector into a discrete job permutation. In brief, such a rule, which works by employing a sorting procedure, enables the algorithm to switch from a conventional scheme to a discrete one. For example,

Table 2. Illustrative example of encoded solution x and permutation solution$\pi$ .

$\pi$	1	2	3	4	5	6	7	8
x	0.14	−0.74	1.11	0.98	2.32	−1.54	−2.24	0.78

shows a generic real-encoded harmony x corresponding to the permutation solution $\pi$ = {7-6-2-1-8-4-3-5} after the SPV conversion is executed.

The computational procedure of HS is explicated by the pseudo-code in Algorithm 1. A preliminary initialization phase is needed to set the initial values of three control parameters, i.e., the harmony memory consideration rate (HMCR), the pitch adjustment (PAR) and the bandwidth (BW). Then, a set of randomly generated solutions (harmonies) are stored in the Harmony Memory (HM). The number of the harmony solutions contained in HM is denoted as harmony memory size (HMS). The generation of the initial population (i.e., the initial HM) may assume a strategic role in the searchability of an evolutionary algorithm. To enhance the quality of the initial HM, two harmonies are replaced by two well-known heuristics, namely Short Processing Time (SPT) and Long Processing Time (LPT), which consider only the treatment times.

Algorithm 1 Harmony Search

1:

Step 1:

Initialization and setting of control parameters, namely HMS, HMCR,

2:

PAR, BW,$({LB}_j{,UB}_j$$)\forall j=1,\dots ,n$; iter = 0;

3:

Step 2:

Generate the initial population, i.e., the HM, and calculate the objective

4:

function of each harmony vector.

5:

Step 3:

$\mathrm{Improvise}a\mathrm{new}\mathrm{harmony}{x}^{new}$ as follows:

6:

for i = 1: HMS

7:

for j = 1: n

8:

if rand < HMCR

9:

$x_{ij}^{new}{=x}_{aj}$ where a $\in$ (1, …, HMS)

10:

if rand < PAR

11:

$x_{ij}^{new}{=x}_{ij}^{new}\pm \mathit{rand}·\mathit{BW}$

12:

end

13:

else

14:

$x_{ij}^{new}{=LB}_j+\mathit{rand}·{(UB}_j{-LB}_j)$

15:

end

16:

end

17:

end

18:

Step 4:

$\mathrm{Compute}\mathrm{the}\mathrm{objective}\mathrm{function}f{x}^{new}$). Update evals and iter.

19:

Step 5:

$\mathrm{Update}\mathrm{HM}\mathrm{by}{x}^{\mathit{worst}}$ $\leftarrow$ $x^{\mathit{new}}$ if f($x^{\mathit{new}}$) < f($x^{\mathit{worst}}$)

20:

Step 6:

if the exit criterion is satisfied

21:

$\mathrm{Stop}\mathrm{the}\mathrm{algorithm};\mathrm{Return}{x}^{\mathit{best}}$ and f($x^{\mathit{best}}$)

22:

else

23:

Goto Step 3

24:

end

After the initial HM is created, the decoding procedure is employed by the Harmony Search for evaluating the HM (see Section 3.3). A variable evals is updated to record the number of evaluations carried out by the algorithm. Then, a new harmony vector $x^{\mathit{new}}$ is stochastically generated. If the new candidate harmony performs better than the worst one $x^{\mathit{worst}}$ in the HM, the latter is replaced by the new one. In this fashion, the harmony memory is constantly updated. A variable iter is used to record the number of improvisations. To boost the searchability of the proposed metaheuristic, two computational techniques were embedded into the HS structure, denoted as Local Search and Reinitialization, properly described in the following subsections. Finally, the HS algorithm stops once one of the termination criteria is satisfied. We decided to implement two triggers, the former being the maximum number of evaluations (Max_ev), the latter depending on the maximum computational time (Max_ct).

3.1. Local Seach

In the local search, at each iteration, a harmony $x_r$ is randomly extracted from the current harmony memory, thus working as starting seed of this procedure. Two well-known perturbation methods, namely insertion and swap, are applied to the seed according to an adaptive probability equal to 1 − [(iter/Max_iter)], where Max_iter is equal to (P∙1000)/HMS. In brief, with the insertion method being more explorative than swap, it has a higher probability to be used at the early stages of the evolutionary path. Insertion consists of randomly selecting a digit and inserting that into a random position of the harmony vector. Swap means to exchange two randomly selected digits of the harmony vector. If the perturbed harmony $x_s$ performs better than the original one, the seed is replaced. The local search is outlined in Algorithm 2. Moreover, if the new harmony improves the best solution achieved so far, the new local optimum is updated and the worst solution of the harmony memory is replaced by the new best solution properly reinitialized.

Algorithm 2 Local Search
1:	Step 1:		$\mathrm{Randomly}\mathrm{select}a\mathrm{harmony}\mathrm{in}\mathrm{the}\mathrm{current}\mathit{HM}:x_r$|r $\in$ int(U [1,HMS])
2:	Step 2:		flag_impr = 0
3:	Step 3:		for q = 1 : n
4:				if$\mathit{rand}$< 1 − (iter/Max_iter)
5:					$x_s$ $\leftarrow$ insertion($x_r$)
6:				else
7:					$x_s$ $\leftarrow$ swap($x_r$)
8:				end
9:				if$f{x}_s$) < f($x_r$)
10:					$x_r$ $\leftarrow$ $x_s$; f($x_r$)$\leftarrow$ f($x_s$)
11:					if$f{x}_s$) < $F_{\mathit{best}}$
12:						flag_impr = 1
13:						$x^{\mathit{best}}$ $\leftarrow$ $x_s$; f($x^{\mathit{best}}$)$\leftarrow$ f($x_s$)
14:						$x_s^{*}$ $\leftarrow$ reinitialize($x_s$)
15:						$x^{\mathit{worst}}$ $\leftarrow$ $x_s^{*}$
16:						$F^{\mathit{worst}}$ $\leftarrow$ f($x_s$)
17:					end
18:				end
19:			end

3.2. Reinitialization

Reinitialization is a novel strategy we propose to enhance the exploration ability of HS. Once a reference harmony is selected, it generates a random harmony and rearranges the digits so that the SPV rule yields the same permutation solution of the reference harmony. Regardless of whether the local search does not yield any improvement in the local optimum, the whole harmony memory is reinitialized. The reinitialization strategy aims at regenerating the information stored in each harmony vector; thus, reinitializing the current HM would reduce the risk of remaining prematurely trapped into a local optimum. Every permutation solution obtained by applying the SPV rule to the real encoded vector just depends on the sorted information and not on the specific values assumed by each digit. Hence, similarly to a well-known restart mechanism often embedded in the metaheuristic algorithms [44], reinitialization would allow the algorithm to improve the exploration ability over the space of solutions. For the sake of clarity, the reinitialization procedure is reported in Algorithm 3.

Table 3. Example of Harmony Reinitialization.

Vectors	Values
x	0.45	0.62	−0.12	0.33
$x_{\mathit{perm}}$	3	4	1	2
$x_r$	0.61	0.27	0.82	−0.23
$x_{\mathit{rein}}$	0.61	0.82	−0.23	0.27

clarifies how a reference harmony x can be reinitialized by $x_{\mathit{rein}}$, which in turn is obtained by sorting the information of a new randomly generated vector $x_r$.

Algorithm 3 Harmony Reinitialization
1:	Step 1:	Select a harmony x
2:	Step 2:	Generate the corresponding permutation harmony by applying SPV on x:[~,$x_{\mathit{perm}}$] = sort(x)
3:	Step 3:	Generate a random harmony $x_r$
4:	Step 4:	Sort $x_r$ values: $x_{r\text{\_}\mathit{sort}}$ = sort($x_r$)
5:	Step 5:	Sort $x_r$ values through: $x_{\mathit{perm}}$: $x_{\mathit{rein}}$($x_{\mathit{perm}}$) = $x_{r\text{\_}\mathit{sort}}$
6:	Step 6:	Reinitialize the harmony: x $\leftarrow$ $x_{\mathit{rein}}$

3.3. Decoding Procedure

Any metaheuristic algorithm generates several candidate solutions for the problem under investigation to identify the best near-optimal solution. In the case of the COS problem, any solution is represented by a permutation sequence of patients $\pi$ to be processed for every scenario ω. Consequently, a decoding procedure is needed to evaluate the expected total flow time associated with $\pi$. The decoding procedure for the multi-stage COS problem at hand entails two distinct phases: the Early-Start (ES) scheduling and the Late-Start (LS) scheduling. The former follows a regular scheduling strategy; every patient preliminarily assigned to a specific oncologist starts the visit as early as possible, when his/her referee oncologist is available for the medical consultation. Once all patients have been processed at the consultation stage, the First Come First Served (FCFS) policy is applied to schedule drug preparations at the pharmacy stage and again the patients at the treatment stage, conforming to the delivery time of the drugs arranged in batches. The latter phase aims to reduce patient waiting time. The LS approach works by adjusting the ES schedule based on a backward rule. In brief, the activities scheduled at the treatment stage remain the same while, going backward from the drug preparation to the consultation stage, the related operations are shifted as much as possible ahead, conforming to the provided constraints. As a result, the starting time of the medical consultation for some patients can be postponed, thereby favoring a consequent reduction of the total flow time for each scenario. At the same time, the chemotherapy outpatient schedule is generated by matching the consultation starting time with the patient arrival times. Finally, the expected total flow time E(F) can be computed by considering the single flow time contributions of each scenario ω (see Equation (3)).

4. Numerical Experiments and Analysis of Results

To prove the effectiveness of the proposed HS in solving the COS problem under investigation, several data sets of numerical instances were generated. We considered three classes of problems at varying problem sizes, namely small (P = 15), medium (P = 40) and large (P = 70). For each of them, a Design of Experiments (DOE) was performed by involving four factors (i.e., number of oncologists O, number of chairs C, drug batch size CAP) and the interval of the uniform distribution related to the batch delivery time U(TD) varied at two levels, low (L) and high (H), respectively. Hence, a total amount $2^4=16$ configurations for each class of problem were considered. Since ten instances were randomly generated for each configuration, $3·16·10=480$ runs were executed.

Table 4. Classes of problems and related parameters.

CLASS	Small (P = 15)		Medium (P = 40)		Large (P = 70)
Factor/Levels	L	H	L	H	L	H
O	2	3	4	6	6	9
C	5	10	10	15	15	20
CAP	2	4	2	4	2	4
U(TD)	U [8,12]	U [18,22]	U [8,12]	U [18,22]	U [8,12]	U [18,22]

reports the different factors/levels, whose values are fixed based on the benchmark problems addressed in the literature so far (please see [11,16,23,25,30,35,45] among others). Indeed, since both the drug transportation- and the batch size-related issues have never been investigated in the literature so far, their values were set based on a brief survey involving the medical staff. The number of nurses N was set by considering that the ratio N/C is usually fixed to 1/4 to respect the assumption that one nurse can monitor at most four patients simultaneously [2,6]. Finally, to balance the pharmacy capability on the problem size, a different number of pharmacist technicians D = {1,2,3} was assigned to each class of problems, respectively.

Due to the stochastic feature of the proposed COS problem, the number of scenarios Ω to evaluate any candidate solution was set to 300. The stochastic parameters of the COS problem under investigation were set as follows. At each scenario ω, every patient p was randomly assigned to a referee oncologist ${\mathit{Op}}_p$. The medical consultation time ${\mathit{Dc}}_p^{\omega }$ of patient p at scenario ω depends on an N(22.83, 3.19) normal distribution. The deferral coefficient $\delta {}_p^{\omega }$ of patient p at scenario ω arises from ${\lambda }_p^{\omega }\in U[0,1]$ with $\overline{\lambda }$ equal to 0.20 [25]. The drug preparation time ${\mathit{Dt}}_p^{\omega }$ was extracted from an uniform distribution U(3, 7). The delivery time ${\mathit{TD}}_b^{\omega }$ of batch b at scenario ω was also derived from a uniform distribution as in Table 4. The set-up time ${\mathit{Ds}}_p^{\omega }$ of patient p at scenario ω was derived from U(5, 15). Finally, the treatment time ${\mathit{Di}}_p^{\omega }$ was drawn from a gamma distribution Γ(1.9, 52.37).

The type of statistical distribution mentioned above refers to the results of an extensive time study carried out in a chemotherapy unit located in Southern Italy. Each metaheuristic algorithm was coded in Matlab^®R2021b and executed on a 4GB RAM-2 processors virtual machine embedded on a workstation equipped with an INTEL i9-9900 3.6 GHz 10 core CPU, 32Gb DDR4 2,666MHz RAM and Win 10 PRO OS. The exit criteria were set based on preliminary tests as follows: Max_ev = 3000 and Max_ct = 300 s for P = 15; Max_ev = 9000 and Max_ct = 1000 s for P = 40; Max_ev = 9000 and Max_ct = 2500 s for P = 70.

Although the expected total flow time E(F) was the objective function of the COS problem under investigation, the Relative Percentage Deviation (RPD) function was handled to compare the results obtained by the tested optimization techniques.

$$\mathit{RPD}=\frac{{\mathit{ALG}}_{\mathit{sol}}-{\mathit{BEST}}_{\mathit{sol}}}{{\mathit{BEST}}_{\mathit{sol}}}·100$$

(4)

where ${\mathit{ALG}}_{\mathit{sol}}$ is the E(F) value achieved by a certain algorithm, while ${\mathit{BEST}}_{\mathit{sol}}$ is the best-expected solution achieved by the algorithms related to the same instance.

This section deals with the numerical analyses we performed to demonstrate the ability of the proposed HS in solving the stochastic COS problem. First, a preliminary computational analysis to calibrate the control parameters of HS was carried out. Since the proposed metaheuristics employ both a set of heuristics for improving the initial population and a specific local search during the evolutionary path, their effect on the quality of solutions was tested. Finally, an extended comparative campaign was employed for assessing the effectiveness of the tested metaheuristic algorithm in the COS problem.

4.1. Calibration of the Hybrid Harmony Search

The effectiveness of any metaheuristic algorithm strongly depends on the values assigned to control parameters, which should assure a suitable balance between exploration and exploitation. To this end, this section aims to calibrate the hybrid HS. Conforming to the seminal paper of [46] and after a series of trial-and-error tests, we set the HMS to 60. The rest of the control parameters, namely HMCR, PAR and BW, were varied at three levels, as reported in

Table 5. Experimental plan for the calibration of the hybrid Harmony Search.

Factor	Description	Values
HMS	HM Size	60
HMCR	HM Consideration Rate	0.50, 0.70, 0.90
PAR	Pitch Adjustment Rate	0.20, 0.50, 0.80
BW	Bandwidth	0.001, 0.01, 0.10

, and a full factorial plan was engaged to select the most suitable control parameters. For each configuration and each class of problems, three instances were randomly generated based on the full-factorial DOE in Table 4.

To sum up, $3^3·3·3=243$ runs were executed for calibration purposes and the RPD measure was adopted as the response variable. For the sake of brevity, the outputs from the ANOVA analysis, whose related model result was significant with a p-value = 0.000, were omitted, while the main effects plot is depicted in Figure 1. The selected values are as follows: HMCR = 0.90; PAR = 0.20; BW = 0.10.

4.2. Comparing Different Variants of HS

To support the choice of the proposed algorithm, three alternative configurations of HS were compared with the proposed metaheuristic. A benchmark of 15 instances involving three classes of problems and based on a full factorial experimental plan, as in Table 4, was engaged. The first configuration, denoted by HS_NH, consists of the HS with no heuristic solutions in the initial population. In the second variant of HS, named HS_NL, the local search is disabled, while in the last algorithm, denoted as HS_NHL, both heuristic solutions and local search are excluded.

Table 6. Design of HS: RPD values.

Instance	HS [%]	HS_NH [%]	HS_NL [%]	HS_NHL [%]
1	0.00	0.00	0.05	0.01
2	0.00	0.02	0.02	0.36
3	0.03	0.02	0.04	0.00
4	0.00	0.00	0.07	0.43
5	0.00	0.00	0.04	0.04
6	0.00	0.03	0.12	0.13
7	0.08	0.00	0.49	2.37
8	0.00	0.00	0.18	0.20
9	0.00	0.02	0.17	0.14
10	0.02	0.00	0.06	0.10
11	0.00	0.00	0.80	6.06
12	0.00	0.15	1.72	2.51
13	0.00	0.00	0.07	0.06
14	0.01	0.00	0.31	0.30
15	0.00	0.00	0.13	0.55
${\mathit{RPD}}_{\mathit{med}}$	0.01	0.02	0.28	0.88
${\mathit{RPD}}_{\mathit{max}}$	0.08	0.15	1.72	6.06

compares the different algorithms in terms of RPD for each instance. Moreover, regardless of the problem size, the global median on the RPDs (${\mathit{RPD}}_{\mathit{med}}$) and the maximum RPD (${\mathit{RPD}}_{\mathit{max}}$) are reported. As the reader can notice, the positive effect of heuristic solutions on the initial population is quite weak. On the other hand, the local search significantly affects the quality of solutions. However, the bold values in Table 6 confirm the outperformance of hybrid HS with respect to the rest of the competitors and justify the use of heuristic solutions in the initial population and the local search as well.

4.3. Extended Comparison Campaign

The two different metaheuristics, i.e., hybrid HS and GRASP, were compared based on a newly generated full factorial experimental plan, as indicated in Table 4. To sum up, three classes of problems at varying Ps were considered, each one involving 16-factor configurations and 10 numerical instances that were randomly generated. Hence, 480 runs were elaborated by each algorithm. For each class of problem,

Table 7. Median RPDs and global indicators to compare HS and GRASP.

Config./Class and Algorithm	Small		Medium		Large
	HS	GRASP	HS	GRASP	HS	GRASP
1	0.00%	0.13%	0.00%	0.47%	0.00%	0.84%
2	0.00%	0.16%	0.00%	0.55%	0.00%	1.23%
3	0.00%	0.11%	0.00%	0.41%	0.00%	0.23%
4	0.00%	0.29%	0.00%	0.22%	0.00%	0.64%
5	0.00%	0.19%	0.00%	0.55%	0.00%	0.68%
6	0.00%	0.22%	0.00%	0.31%	0.00%	0.39%
7	0.00%	0.22%	0.00%	0.14%	0.00%	0.72%
8	0.00%	0.10%	0.00%	0.36%	0.00%	0.29%
9	0.00%	0.30%	0.00%	0.37%	0.00%	0.33%
10	0.00%	0.12%	0.00%	0.51%	0.00%	0.76%
11	0.00%	0.19%	0.00%	0.18%	0.00%	0.54%
12	0.00%	0.25%	0.00%	0.29%	0.00%	0.88%
13	0.00%	0.18%	0.00%	0.19%	0.00%	0.76%
14	0.00%	0.16%	0.00%	0.36%	0.00%	0.51%
15	0.00%	0.14%	0.00%	0.35%	0.00%	0.28%
16	0.00%	0.11%	0.00%	0.13%	0.00%	0.77%
${\mathit{RPD}}_{\mathit{med}}$	0.00%	0.18%	0.00%	0.27%	0.00%	0.57%
${\mathit{RPD}}_{\mathit{max}}$	0.00%	0.30%	0.00%	0.55%	0.00%	1.23%

reports the median RPD values for each configuration, computed over the provided 10 instances, each one entailing 300 scenarios.

The numerical outputs reveal the outperformance of HS over the GRASP algorithm. The minimum RPD values for each configuration are highlighted in bold and all of them are in the HS-related column. Regardless of the problem configurations, the median RPD (${\mathit{RPD}}_{\mathit{med}}$) and the maximum RPD (${\mathit{RPD}}_{\mathit{max}}$) support the primacy of HS over the GRASP algorithm. Looking at the ${\mathit{RPD}}_{\mathit{max}}$ values it is clear that the GRASP technique worsens as the complexity of the problem rises.

Finally, the HS algorithm was compared with the LPT rule, which is the strategy adopted by the staff of the oncology unit under investigation. A total of 20 instances, based on the real-life parameter values, were used to compare the LPT rule and HS algorithm.

Table 8. Comparison between LPT rule and HS algorithm.

Instance	LPT [%]	HS [%]
1	5.12	0.00
2	6.36	0.00
3	12.90	0.00
4	8.66	0.00
5	6.90	0.00
6	3.90	0.00
7	3.01	0.00
8	3.22	0.00
9	5.78	0.00
10	3.06	0.00
11	11.98	0.00
12	11.27	0.00
13	5.69	0.00
14	9.84	0.00
15	13.39	0.00
16	2.89	0.00
17	10.37	0.00
18	5.56	0.00
19	4.58	0.00
20	4.31	0.00
${\mathit{RPD}}_{\mathit{ave}}$	6.94	0.00

shows the results in terms of RPDs. The experiments demonstrate that the HS algorithm outperforms the LPT by approximately 7% on the average RPD (${\mathit{RPD}}_{\mathit{ave}}$). These results are consistent with the comparison between the metaheuristic algorithm and the LPT method shown in the work of Garaix et al. [25].

5. Conclusions

5.1. Final Discussion

In this research, the same-day offline stochastic chemotherapy outpatient appointment scheduling problem, inspired by a real-world oncology department, was investigated. Differently from the rest of the literature on this topic, we stochastically modeled all the stages provided by the chemotherapy process and, in addition, several sources of uncertainty (e.g., deferrals and medical consultation times) were taken into account. Particularly, since the pharmacy is located far away from the treatment unit, we considered the real-life scenario in which a drug delivery service is a time-consuming task needed to take the therapies to the ward.

A stochastic scheduling approach was adopted to cope with the uncertainty of the problem. Since the problem under investigation can be assimilated into a hybrid flow shop scheduling problem with resource-related constraints, several idle times may occur among the stages. Therefore, to improve the quality of any appointment schedule, we implemented an LS scheduling strategy in the decoding algorithm for the minimization of the objective function, i.e., the expected total flow time.

The outpatient scheduling problem is extensively studied in the literature and most authors used mathematical programming applied to relaxed models for the generation of optimal solutions. Since the problem under investigation is NP-hard in the strong sense, we developed a hybrid Harmony Search named HS. To demonstrate both the efficacy and efficiency of the proposed metaheuristic, a comparative analysis involving a GRASP algorithm from the relevant literature on the same topic was carried out.

The analysis pointed out the effectiveness of the proposed hybrid HS algorithm. There are no instances in which the GRASP algorithm can perform better than HS. Since the GRASP algorithm works by successive constructions of a greedy randomized solution, which in turn is improved by a semi-greedy constructive procedure, the time required to build a solution is considerably high. Therefore, the weakness of GRASP in comparison with the other techniques can be mostly explained by the lower number of solutions it evaluates within Max_ct. In fact, the HS algorithm is able to perform almost 5000 evaluations more than GRASP in 2500 s.

The quality of the solutions assured by the hybrid HS, as well as its computational efficiency, was tested by the staff of the chemotherapy clinic, which decided to replace the LPT appointment strategy adopted so far with the proposed hybrid metaheuristic approach.

The numerical results demonstrate that the use of the proposed outpatient scheduling method may positively affect oncology units as follows:

• The proposed method, which combines efficacy and computational efficiency, assures a significant improvement in terms of total flow time with respect to the LPT rule, which in turn is a scheduling strategy commonly adopted by healthcare managers;
• From an economic viewpoint, solutions obtained by the hybrid HS allow for reducing the clinics’ idle time and the number of patients treated in a day as well;
• Considering the patient's perspective, the proposed method yields a significant reduction in patient waiting time and on a global service level as well.

5.2. Limitations and Directions for Future Research

The present paper is characterized by assumptions that can be removed in future research. Firstly, one of the main assumptions is related to blood exams. In this paper, the patient executes the blood draw in the previous days, and then the results are immediately available for medical consultation. However, there are oncology units that manage the blood draws on the same day of the medical consultation and chemotherapy treatment. Therefore, in future research, the proposed scheduling approach may be tested in four-stage oncology units, wherein the durations of the blood exams and the related resources are considered.

Moreover, the hybrid HS was tested in a same-day oncology unit. It would be interesting to assess the effectiveness of this methodology in next-day oncology units, in which medical consultations and chemotherapy treatments are provided on different days.

Alternative metaheuristic algorithms, a multi-objective approach involving more objective functions or further constraints on the system modeling could be considered as new opportunities for future research in the chemotherapy outpatient scheduling topic. In particular, the reduction of overtime may be considered a new objective function of the problem. Finally, it could be valuable to separately evaluate the waiting time for therapy preparation, chair and nurse availability and their impact on the total flow time.