Physician scheduling for outpatient department with nonhomogeneous patient arrival and priority queue

Abstract

The growing demand for outpatient departments and prolonged waiting time of patients in recent years has made physician scheduling necessary to provide timely medical services. This study focused on the problem of scheduling physicians in an outpatient system with a nonhomogeneous patient arrival and priority queue, which exists in many Asian hospitals. In such a system, both patients with appointments and walk-in patients wait in a priority queue to see physicians, with the patients’ arrivals fluctuating throughout the day. In order to respond to time-related demands while simultaneously respecting physicians’ preferences for being on or off duty in some specific slots, a staffing optimization model was formulated. In addition, a physician rescheduling model was proposed for the case where a physician is unexpectedly absent. To solve the problem, a calibrated waiting time approximation-based genetic algorithm methodology was proposed. Its main contribution is the use of a data-driven analytical method to estimate the average waiting time of the two types of patients in a complex queuing system. The results of a numerical study and real case study in a Shanghai hospital showed that physician scheduling optimization on the basis of the proposed waiting time approximation method was effective and efficient when applied to the proposed outpatient system.

Appendices

Appendix A: Proof of total waiting time of the appointment patients in the first slot

The total waiting time of appointment patients in the first slot includes three parts. The first part is the waiting time estimation of early patients before the first appointment patient receives medical service. Because the integration of the left part of the Normal distribution is equal to 0.5, we let the normal probability density function multiply $2{\lambda }_{1}^{A}\Delta {e}_{1}^{A}$ to be the patient frequency distribution. It promises that the integration of the frequency distribution of the early appointment patients’ arrivals is equal to the total number of early appointment patients. Then the total waiting time could be calculated as $2{\lambda }_{1}^{A}\Delta {e}_{1}^{A}\underset{-\Delta }{\overset{0}{\int }}(-x)\frac{1}{\sqrt{2\pi }{\sigma }_{t}}{e}^{\frac{{-x}^{2}}{2{{\sigma }_{t}}^{2}}}dx$.

The second part is the waiting time estimation when the early appointment patients start to receive service. In the general case with any number ${s}_{1}$ of physicians, the first ${s}_{1}$ early appointment patients do not wait, all physicians remain busy, and patients depart according to a Poisson process with rate ${s}_{\mathrm{t}}\mu$ until the rest of early patients begin to be served. The mean waiting time for patient (${s}_{1}+i$) is $\frac{i}{{s}_{1}\mu }$, thus their total waiting time is $\sum_{i=1}^{{{(\lambda }_{1}^{A}\Delta {e}_{1}^{A}-{s}_{1})}^{+}}\frac{1}{{s}_{\mathrm{t}}\mu }*i$.

The last part is the total waiting time of new arriving appointment patients who have been served before the slot $t$ ends. It can be estimated by multiplying the number of patients and the average waiting time of patients. The average waiting time is obtained by applying no preemptive priority queueing theory as Eq. (18).

B: Proof of total waiting time of the appointment patients in slot ${\varvec{t}}({\varvec{t}}\ne 1)$

The total waiting time of appointment patients in slot $t$ includes four parts. The first two parts are the total waiting time estimation of delayed patients and the early appointment patients, which are the same as Eq. (25). The third part is the total waiting time of new arriving appointment patients. The last part is the total approximation of patients who have been delayed to the next slot without being serviced. For a stationary Poisson arrival process of rate ${\lambda }_{t}^{A}$, the mean inter-arrival times are all equal to $\frac{1}{{\lambda }_{t}^{A}}$. Thus the waiting time of the ith of these patients is $\frac{i}{{\lambda }_{t}^{A}}$ and the total waiting time of these patients is $\sum_{i=1}^{{AL}_{\mathrm{t}}^{END}}\frac{1}{{\lambda }_{t}^{A}}*i$, where $AL_{t}^{END}$ is obtained by a similar equation as (27).

C: Details explanation of the algorithm for the walk-in patients’ waiting time approximation in slot ${\varvec{t}}({\varvec{t}}\ne 1)$

To further explain the idea of walk-in patients’ total waiting time approximation during a slot, we take a specified case in the slot $t$ as an example. Figure 8 provides a specified patient flow, where arrows above the timeline indicate the patient arrival process and arrows below the timeline indicate the departure process.

Fig. 8

A specified case of patient arrival in slot $t$

Nodes $a,z$ are the time points that the slot starts and ends, respectively. The time duration of the two nodes is ∆. In this case, there is one patient with an appointment and three walk-in patients waiting before the start of the slot, which is denoted as patients 1–4. Node $\mathrm{c}$ is the time point when all delayed patients have finished their service. During the process, the delayed walk-in patients will begin to receive service at a time point $b$ after the delayed patient with appointment (patient 1) and the newly arrived patient with appointment (patient 5) departs the system. These delayed walk-in patients are defined as batch 1. Certainly, at the time point $b$, some walk-in patients (patient 6) who arrive earlier also waits in the queue, and we define the patient as batch 2. The time point $d$ is the service completion time of batch 2. In this way, the total slot is divided into smaller parts and the endpoint of the part (except the last part) is the service completion time for the patient batch prior to the previous batch. Note that the last part’s end point is the end time of the slot $t$, thus the patients (patient $\mathrm{i}-1$, patient $\mathrm{i}$ and patient $\mathrm{j}$) who have not been served are delayed to the next slot. The total waiting time for each batch is calculated and then the total waiting time of the walk-in patients of slot $t$ is obtained by summing up the results of all batches.

$${ }TWW_{{\text{t}}} = \mathop \sum \limits_{n = 1}^{{N_{{\text{t}}} }} TWW_{{\text{t}}}^{n}$$

(38)

where $n$ is the batch index and ${N}_{\mathrm{t}}$ is the total number of batches in the slot $t$. ${TW{W}_{\mathrm{t}}}^{n}$ is the total waiting time of walk-in patients in the batch $n$ of the slot $t$.

For the total waiting time of each batch, the general performance is to calculate three parts of the waiting time: the first part is the total waiting time before the last patient of the batch arrives at the hospital; the second part is the total waiting time when all patients wait in the queue; the third part is when patients begin to receive service in sequence. Thus,

$$TWW_{t}^{n} = \mathop \sum \limits_{i = 1}^{{WL_{t}^{n} }} i{*}\lambda_{t}^{W} + WL_{t}^{n} {*}WS_{t}^{n - 1} + \mathop \sum \limits_{i = 1}^{{WL_{t}^{n} }} \left( {i - 1} \right){*}ES_{t}$$

(39)

$$\mathrm{W}{{L}_{t}}^{n}={\lambda }_{t}^{W}*{W{S}_{t}}^{n-2}$$

(40)

$${{WS}_{t}}^{n-1}={{WL}_{t}}^{n-1}*{ES}_{t}$$

(41)

$${ES}_{t}=\frac{1}{{s}_{t}\mu -{\lambda }_{t}^{A}}$$

(42)

where ${{WL}_{t}}^{n}$ is the number of walk-in patients in the batch $n$ of the slot $t$. ${ES}_{t}$ is the expected service time of one walk-in patient of the slot $t$. ${W{S}_{t}}^{n}$ is the service time of walk-in patients in batch $n$ of the slot $t$.

The total waiting time calculation for the first batch and last batch is slightly different. For the first batch, they are delayed from the previous slot, therefore, the waiting time is only considered when all patients wait in the queue, and when they receive sequential service as

$$TWW_{t}^{1} = { }WS_{{\text{t}}}^{0} {*}WL_{{\text{t}}}^{1} + \mathop \sum \limits_{i = 1}^{{WL_{{\text{t}}}^{1} }} \left( {i - 1} \right){*}ES_{t}$$

(43)

$$\mathrm{W}{{L}_{t}}^{1}={{WL}_{t}}^{\mathrm{START}}$$

(44)

$${{WS}_{t}}^{0}=\frac{{({{AL}_{t}}^{START}-{s}_{t}+1)}^{+}}{{s}_{t}\mu -{\lambda }_{t}^{A}}$$

(45)

where ${{WL}_{t}}^{\mathrm{START}}$ denotes the number of walk-in patients in the waiting queue when the slot $t$ starts.

For the last batch, some patients (as patient i, i−1 and j in Fig. 6) may not be served and thus be delayed to the next slot. The delayed patients’ waiting time should be included in the batch performance as

$$TWW_{t}^{{N_{t} }} = \mathop \sum \limits_{i = 1}^{{{WL_{t}{N_{t} }} }} i{*}\lambda_{t}^{W} + \lambda_{t}^{W} {*}WS_{t}^{{N_{t} - 2}} {*}WS_{t}^{{N_{t} - 1}} + \mathop \sum \limits_{i = 1}^{{{WS_{t}{N_{2} }} /E\left( S \right)}} \left( {i - 1} \right){*}ES_{t} + {\text{W}}L_{t}^{END} { *}WS_{t}^{{N_{t} }}$$

(46)

$$\mathrm{W}{{L}_{t}}^{{N}_{t}}={\lambda }_{t}^{W}*{(W{S}_{t}}^{{N}_{t}-2}+{W{S}_{t}}^{{N}_{t}-1}+{W{S}_{t}}^{{N}_{t}})$$

(47)

$$\mathrm{W}{{L}_{t}}^{END}= W{{L}_{t}}^{{N}_{t}}-{{WS}_{t}}^{{N}_{t}}/{ES}_{t}$$

(48)

$${ }WS_{t}^{{N_{t} }} = \Delta - \mathop \sum \limits_{n = 1}^{{N_{t} - 1}} WS_{t}^{n}$$

(49)

where $\mathrm{W}{{L}_{t}}^{END}$ denotes the number of walk-in patients in the waiting queue when the slot $t$ ends. Then the number of walk-in patients served in slot $t$ and delayed to the next slot are

$${WL}_{\mathrm{t}}={\sum }_{n=1}^{{N}_{\mathrm{t}}}{W{L}_{\mathrm{t}}}^{n}-{{WL}_{t}}^{\mathrm{END}}$$

(50)

$${WL}_{t}^{END}={WL}_{t+1}^{START}$$

(51)

For proof of the batch’s total waiting time calculation, please refer to Appendix D, E, F & G.

Note that during the process of the above algorithm, no patient may arrive during two adjacent time points if ${\lambda }_{t}^{A}+{\lambda }_{t}^{W}<{s}_{t}\mu$. In such a scenario, the algorithm at this empty batch is stopped and the stationary formula is taken into effect (equation of (18)–(24)) to make a waiting time estimate for rest patients. Define the rest of the patients as the batch ${N}_{t}$, and then their total waiting time is

$$TWW_{t}^{{N_{t} }} = \lambda_{t}^{W} (\Delta - \mathop \sum \limits_{n = 1}^{{N_{t} }} WS_{2}^{n} ){*}AAW_{t}$$

(52)

D: Proof of total waiting time calculation for a general walk-in patients batch

The total waiting time of a general walk-in patients batch can be divided into three parts: the first part is the total waiting time before the last patient of the batch arrives at the hospital, which is similar to the last part of Eq. (31). The second part is the total waiting time when all patients wait in the queue, which multiplies the number of walk-in patients and the total service time of the last patients’ batch. The third part is when patients start getting service in sequence. The mean service time of a patient is ${ES}_{t}$, thus the waiting time for ith patient is $\left(i-1\right)*{ES}_{t}.$

For the derivation of ${ES}_{t}$, the scenario considers that there are no appointment patients and two walk-in patients in the queue. The first walk-in patient starts service immediately and finishes service after $1/{s}_{t}\mu$. If no appointment patients arrive during this slot, then the second walk-in patient could get service immediately. Therefore, $E\left({ES}_{t}|0\right)=1/{s}_{t}\mu$. If one appointment patient arrives during this slot, then the second walk-in patient should wait until the appointment patient finishes his or her service. Then, $E\left({ES}_{t}|1\right)={1/s}_{t}\mu + E\left({ES}_{t}\right)$. Accordingly, if $j$ appointment patients arrive, the second walk-in patients would wait for $1/{s}_{t}\mu + j*E\left({ES}_{t}\right)$.The probability of such a case is $\frac{{\left(\frac{{\lambda }_{t}^{A}}{{s}_{t}\mu }\right)}^{j}}{j!}{e}^{-\left(\frac{{\lambda }_{t}^{A}}{{s}_{t}\mu }\right)}$. Thus we could obtain

$$\begin{aligned} ES_{t} & = \mathop \sum \limits_{j = 0}^{\infty } \frac{{\left( {\frac{{\lambda_{t}^{A} }}{{s_{t} \mu }}} \right)^{j} }}{j!}e^{{ - \left( {\frac{{\lambda_{t}^{A} }}{{s_{t} \mu }}} \right)}} E\left( {ES_{t} {|}j} \right) \\ & = \mathop \sum \limits_{j = 0}^{\infty } \frac{{\left( {\frac{{\lambda_{t}^{A} }}{{s_{t} \mu }}} \right)^{j} }}{j!}e^{{ - \left( {\frac{{\lambda_{t}^{A} }}{{s_{t} \mu }}} \right)}} \left( {j{*}E\left( {ES_{t} } \right) + \frac{1}{{s_{t} \mu }}} \right) \\ { } & \frac{1}{{s_{t} \mu }} + \frac{{\lambda_{t}^{A} }}{{s_{t} \mu }}E\left( {ES_{t} } \right) = \frac{1}{{s_{t} \mu - \lambda_{t}^{A} }} \\ \end{aligned}$$

(53)

E: Proof of total waiting time calculation for the first walk-in patients batch

For the first batch, the first arrived patients all wait for appointment patients to conclude their appointments and then begin to receive service one by one according to FCFS basis. The waiting time estimate includes two parts. The first part is the total waiting time that all walk-in patients wait for appointment patients, which is similar to the equation of (39). The second part is total waiting time when walk-in patients start receiving service according to FCFS basis.

For the derivation of the service time of the appointment patients in the first part, it is calculated that the delayed appointment patients’ service time as $\frac{{({{AL}_{t}}^{START}-{s}_{t}+1)}^{+}}{{s}_{t}\mu }$. Then, during the delayed appointment patients’ service, there are expected to be ${\lambda }_{t}^{A}*\frac{{({{AL}_{t}}^{START}-{s}_{t}+1)}^{+}}{{s}_{t}\mu }$ new appointment patients arriving and each patients’ service time is expected to be $\frac{1}{{s}_{t}\mu }$. And so forth for the next newly arriving appointment patients. Thus, the total appointment patients’ service time can be obtained as

$$\begin{aligned} WS_{t}^{0} & = \frac{{\left( {AL_{t}^{START} - s_{t} + 1} \right)^{ + } }}{{s_{t} \mu }} + \lambda_{t}^{A} {*}\frac{{\left( {AL_{t}^{START} - s_{t} + 1} \right)^{ + } }}{{s_{t} \mu }}{*}\frac{1}{{s_{t} \mu }} + \ldots + \left( {\lambda_{t}^{A} } \right)^{n} {*}\frac{{\left( {AL_{t}^{START} - s_{t} + 1} \right)^{ + } }}{{s_{t} \mu }}{*}\left( {\frac{1}{{s_{t} \mu }}} \right)^{n} \\ & = \mathop \sum \limits_{i = 0}^{\infty } \left( {\frac{{\lambda_{At} }}{{s_{t} \mu }}} \right)^{i} \frac{{\left( {AL_{t}^{START} - s_{t} + 1} \right)^{ + } }}{{s_{t} \mu }} = \frac{{\left( {AL_{t}^{START} - s_{t} + 1} \right)^{ + } }}{{s_{t} \mu - \lambda_{t}^{A} }} \\ \end{aligned}$$

(43)

F: Proof of total waiting time calculation for the last walk-in patients batch

The total waiting time of the last batch includes four parts: the first part is the total waiting time when patients arrive according to FCFS, which is calculated as ${\sum }_{i=1}^{W{{L}_{t}}^{{N}_{t}}}i*{\lambda }_{t}^{W}$. The second part is the total waiting time when patients who arrive during ${W{S}_{t}}^{{N}_{t}-2}$ and keep waiting for the batch ${N}_{t}-1$.This total waiting time is calculated by ${\lambda }_{t}^{W}*{W{S}_{t}}^{{N}_{t}-2}*{{WS}_{t}}^{{N}_{t}-1}$.The third part is the total waiting time when some patients receive service according to FCFS. The number of patients is ${{WS}_{t}}^{{N}_{2}}/E(S)$ and their total waiting time is ${\sum }_{i=1}^{{{WS}_{t}}^{{N}_{2}}/E(S)}\left(i-1\right)*{ES}_{t}$. The last part is the total waiting time of patients who were delayed to the next slot. The delayed patients keep waiting during the last batch’s service time, thus their waiting time is the number of patients multiplied by the service time of the last batch as $\mathrm{W}{{L}_{t}}^{END} *{{WS}_{t}}^{{N}_{t}}$.

G: Proof of Proposition 1

Given the total waiting time $TAW\left(TWW\right)$ is summarized by slots through Eqs. (16) and (17), we firstly discuss the scenarios of the single slot.

As for the first slot, the total waiting time of appointment patients ${AW}_{1}$ could be calculated through Eq. (25) as.

$${AW}_{1}=\frac{({\lambda }_{1}^{A}\Delta {e}_{1}^{A}+1)\Delta {e}_{1}^{A}}{2}+\frac{{(\lambda }_{1}^{A}\Delta {e}_{1}^{A}-{s}_{1}+1){(\lambda }_{1}^{A}\Delta {e}_{1}^{A}-{s}_{1})}{2{s}_{1}\mu }+{\lambda }_{1}^{A}\Delta(1-{e}_{1}^{A})*\frac{{s}_{1}\mu }{\{{s}_{1}!*[{s}_{1}\mu -\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)]*\sum_{i=0}^{{s}_{1}-1}\frac{{\left[\mu /\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)\right]}^{{s}_{1}-i}}{i!}+{s}_{1}\mu \}({s}_{1}\mu -{\lambda }_{1}^{A})}$$

Let ${f}_{1}= {s}_{1}\mu -\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)$, ${f}_{2} = \sum_{i=0}^{{s}_{1}-1}\frac{{\left[\mu /\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)\right]}^{{s}_{1}-i}}{i!}$, ${f}_{3}={s}_{1}\mu -{\lambda }_{1}^{A}$, ${g}_{1}=\frac{{s}_{1}\mu }{[{s}_{1}!*{f}_{1}({s}_{1})*{f}_{2}({s}_{1})+{s}_{1}\mu ]{f}_{3}\left({s}_{1}\right)}$, ${g}_{2}= \frac{{(\lambda }_{1}^{A}\Delta {e}_{1}^{A}-1){\lambda }_{1}^{A}\Delta {e}_{1}^{A}}{2{s}_{1}\mu }$ and ${g}_{3}= \frac{({\lambda }_{1}^{A}\Delta {e}_{1}^{A}+1)\Delta {e}_{1}^{A}}{2}$.

Given that ${f}_{1}^{^{\prime}}\left({\lambda }_{1}^{A}\right)={f}_{3}^{^{\prime}}\left({\lambda }_{1}^{A}\right)=-1<0$, ${f}_{2}^{^{\prime}}\left({\lambda }_{1}^{A}\right)=\sum_{i=0}^{{s}_{1}-1}\frac{(i-{s}_{1}){\left[\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)/\mu \right]}^{i-{s}_{1}-1}}{\mu *i!}<0$, thus ${f}_{1}({\lambda }_{1}^{A})$, ${f}_{2}({\lambda }_{1}^{A})$ and ${f}_{3}\left({\lambda }_{1}^{A}\right)$ are decreasing obviously with ${\lambda }_{1}^{A}$. Because ${f}_{1}\left({\lambda }_{1}^{A}\right)>0$, ${f}_{2}\left({\lambda }_{1}^{A}\right)>0$ and ${f}_{3}\left({\lambda }_{1}^{A}\right)>0$, thus ${g}_{1}\left({\lambda }_{1}^{A}\right)=\frac{{s}_{1}\mu }{[{s}_{1}!*{f}_{1}({s}_{1})*{f}_{2}({s}_{1})+{s}_{1}\mu ]{f}_{3}\left({s}_{1}\right)}$ is increasing with ${\lambda }_{1}^{A}$.

Considering ${\lambda }_{1}^{A}\Delta {e}_{1}^{A}-1>0$, ${g}_{2}^{\prime}\left({\lambda }_{1}^{A}\right)=\frac{{(2\lambda }_{1}^{A}\Delta {e}_{1}^{A}-1)\Delta {e}_{1}^{A}}{2{s}_{1}\mu }>0$, thus ${g}_{2}\left({\lambda }_{1}^{A}\right)$ is increasing with ${\lambda }_{1}^{A}$. Besides, ${g}_{3}^{\prime}\left({\lambda }_{1}^{A}\right)= \frac{{\Delta }^{2}{{e}_{1}^{A}}^{2}}{2}>0$, thus ${g}_{3}\left({\lambda }_{1}^{A}\right)$ is increasing with ${\lambda }_{1}^{A}$. Thus ${AW}_{1}={g}_{1}\left({\lambda }_{1}^{A}\right)+{g}_{2}\left({\lambda }_{1}^{A}\right)+{g}_{3}\left({\lambda }_{1}^{A}\right)$ is increasing with ${\lambda }_{1}^{A}$. Similarly, if $\mu$ is increasing,${g}_{1}({s}_{1})$, ${g}_{2}({s}_{1})$ and ${g}_{3}({s}_{1})$ is decreasing. Thus, ${AW}_{1}$ is decreasing with $\mu$.

Assume ${s}_{11}>{s}_{12}>0,$ then.

$$\begin{aligned}{f}_{1}\left({s}_{11}\right)-{f}_{1}\left({s}_{12}\right)&=({s}_{11}-{s}_{12})\mu>0,\\{f}_{2}\left({s}_{11}\right)-{f}_{2}\left({s}_{12}\right)&=\sum_{i={s}_{12}+1}^{{s}_{11}-1}\frac{{\left[\mu /\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)\right]}^{{s}_{11}-i}}{i!}+\sum_{i=0}^{{s}_{12}-1}\frac{{\left[\mu /\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)\right]}^{{s}_{12}-i}}{i!}\{{\left[\mu /\left({\lambda }_{1}^{A}+{\lambda }_{1}^{W}\right)\right]}^{{s}_{11}-{s}_{12}}-1\}>0,\\{f}_{3}\left({s}_{11}\right)-{f}_{3}\left({s}_{12}\right)&=({s}_{11}-{s}_{12})\mu>0,\\{g}_{2}\left({s}_{11}\right)-{g}_{2}\left({s}_{12}\right)&=\frac{{(\lambda }_{1}^{A}\Delta {e}_{1}^{A}-1){\lambda }_{1}^{A}\Delta {e}_{1}^{A}}{2\mu }\left(\frac{1}{{s}_{11}}-\frac{1}{{s}_{12}}\right)<0\end{aligned}$$

Thus, ${f}_{1}\left({s}_{11}\right)>{f}_{1}\left({s}_{12}\right)$, ${f}_{2}\left({s}_{11}\right)>{f}_{2}\left({s}_{12}\right)$ and ${f}_{3}\left({s}_{11}\right)>{f}_{3}\left({s}_{12}\right)$. We could obtain that ${g}_{1}\left({s}_{11}\right)=\frac{\mu }{[{(s}_{11}-1)!*{f}_{1}({s}_{11})*{f}_{2}({s}_{11})+\mu ]{f}_{3}\left({s}_{11}\right)}$<${g}_{1}\left({s}_{12}\right)=\frac{\mu }{[{(s}_{11}-1)!*{f}_{1}({s}_{11})*{f}_{2}({s}_{11})+\mu ]{f}_{3}\left({s}_{11}\right)}$ and ${g}_{2}\left({s}_{11}\right)<{g}_{2}\left({s}_{12}\right)$. Finally, we could obtain that ${AW}_{1}\left({s}_{11}\right)<{AW}_{1}\left({s}_{12}\right)$, which means that ${AW}_{1}$ is decreasing with ${s}_{1}$.The proof of ${WW}_{1}$ is similarly as ${AW}_{1}$.

As for other slot, the monotonicity of ${AW}_{\mathrm{t}}({WW}_{\mathrm{t}})$ is consistent with ${AW}_{1}({WW}_{1})$. Thus, the monotonicity of $TAW(TWW)$ is also consistent with ${AW}_{1}({WW}_{1})$. QED.

H: Proof of Proposition 2

Given the conclusion of ${AW}_{\mathrm{t}}+{WW}_{\mathrm{t}}$ and ${\varphi }_{t}$ maintain consistency among the slots, here we take one slot for illustration.

As for the first slot, ${{\lambda }_{1}^{A}+\lambda }_{1}^{W}={\lambda }_{1}$, $\frac{{\lambda }_{1}^{A}}{{\lambda }_{1}^{W}}={\varphi }_{1}$, then we could obtain that ${\lambda }_{1}^{A}=\frac{{\lambda }_{1}{\varphi }_{1}}{{(1+\varphi }_{1})}$ and ${\lambda }_{1}^{w}=\frac{{\lambda }_{1}}{{(1+\varphi }_{1})}$.

When ${s}_{1}=1$, the total waiting time of appointment patients ${AW}_{1}$ could be calculated through Eq. (25) as

$${AW}_{1}({\varphi }_{1})+{WW}_{1}({\varphi }_{1})=\frac{\sqrt{2}{\lambda }_{1}{\varphi }_{1}\Delta {e}_{1}^{A}{\sigma }_{1}}{\sqrt{\pi }{(1+\varphi }_{1})}({e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}-1)+\frac{({\lambda }_{1}{\varphi }_{1}\Delta {e}_{1}^{A}-1-{\varphi }_{1}){\lambda }_{1}{\varphi }_{1}\Delta {e}_{1}^{A}}{2\mu {{(1+\varphi }_{1})}^{2}}+\frac{{{\lambda }_{1}}^{2}\Delta [{\varphi }_{1}(1-{e}_{1}^{A})\left(\mu -{\lambda }_{1}\right)+\mu ]}{\mu (\mu {\varphi }_{1}-{\lambda }_{1}{\varphi }_{1}+\mu )(\mu -{\lambda }_{1})}$$

$$AW_{1}^{\prime} \left( {\varphi _{1} } \right) + WW^{\prime}_{1} \left( {\varphi _{1} } \right) = \frac{{\sqrt 2 \lambda _{1} \Delta e_{1}^{A} \sigma _{1} e^{{\frac{{ - \Delta ^{2} }}{{2\sigma _{1} ^{2} }}}} }}{{\sqrt \pi (1 + \varphi _{1} )}} + \frac{{\lambda _{1} ^{2} \Delta ^{2} e_{1}^{{A^{2} }} \varphi _{1} }}{{\mu (1 + \varphi _{1} )^{3} }} - \frac{{\lambda _{1} \Delta e_{1}^{A} }}{{2\mu (1 + \varphi _{1} )^{2} }} - \frac{{\lambda _{1} ^{2} \Delta e_{1}^{A} }}{{\left( {\mu \varphi _{1} - \lambda _{1} \varphi _{1} + \mu } \right)^{2} }}$$

When ${\varphi }_{1}=\frac{1}{2}$,

$${AW}_{1}^{\prime}\left({\varphi }_{1}\right)+{W{W}^{\prime}}_{1}\left({\varphi }_{1}\right)=\frac{\sqrt{2}{\lambda }_{1}\Delta {e}_{1}^{A}{\sigma }_{1}{e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{\sqrt{\pi }}+\frac{{{\lambda }_{1}}^{2}{\Delta }^{2}{{e}_{1}^{A}}^{2}}{27\mu }-\frac{{2\lambda }_{1}\Delta {e}_{1}^{A}}{9\mu }-\frac{4{{\lambda }_{1}}^{2}\Delta {e}_{1}^{A}}{{\left(3\mu -{\lambda }_{1}\right)}^{2}}<\frac{\sqrt{2}{\lambda }_{1}\Delta {e}_{1}^{A}{\sigma }_{1}{e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{\sqrt{\pi }}+\frac{{{\lambda }_{1}}^{2}{\Delta }^{2}{{e}_{1}^{A}}^{2}}{27\mu }-\frac{{2\lambda }_{1}\Delta {e}_{1}^{A}}{9\mu }-\frac{{{\lambda }_{1}}^{2}\Delta {e}_{1}^{A}}{{\mu }^{2}}=\frac{{\lambda }_{1}\Delta {e}_{1}^{A}(27\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}+{\lambda }_{1}\Delta {e}_{1}^{A}-6\mu -\frac{27{\lambda }_{1}}{\mu })}{27\sqrt{\pi }\mu }$$

If ${e}_{1}^{A}< \frac{6+\frac{27{\lambda }_{1}}{\mu }-27\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{{\lambda }_{1}\Delta }<\frac{33-27\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{{\lambda }_{1}\Delta }$, ${AW}_{1}^{\prime}\left({\varphi }_{1}\right)+{W{W}^{^{\prime}}}_{1}\left({\varphi }_{1}\right)<0.$

When ${\varphi }_{1}=1$, ${AW}_{1}^{\prime}\left({\varphi }_{1}\right)+{W{W}^{^{\prime}}}_{1}\left({\varphi }_{1}\right)=\frac{\sqrt{2}{\lambda }_{1}\Delta {e}_{1}^{A}{\sigma }_{1}{e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{\sqrt{\pi }}+\frac{{{\lambda }_{1}}^{2}{\Delta }^{2}{{e}_{1}^{A}}^{2}}{8\mu }-\frac{{\lambda }_{1}\Delta {e}_{1}^{A}}{8\mu }-\frac{{{\lambda }_{1}}^{2}\Delta {e}_{1}^{A}}{{\left(2\mu -{\lambda }_{1}\right)}^{2}}$.

$$\begin{aligned}&{AW}_{1}^{\prime}\left(1\right)+{W{W}^{^{\prime}}}_{1}\left(1\right)>\frac{\sqrt{2}{\lambda }_{1}\Delta {e}_{1}^{A}{\sigma }_{1}{e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{\sqrt{\pi }}+\frac{{{\lambda }_{1}}^{2}{\Delta }^{2}{{e}_{1}^{A}}^{2}}{8\mu }-\frac{{\lambda }_{1}\Delta {e}_{1}^{A}}{8\mu }-\Delta {e}_{1}^{A}\\&\quad=\frac{{\lambda }_{1}\Delta {e}_{1}^{A}(8\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}+{\lambda }_{1}\Delta {e}_{1}^{A}-1-\frac{8\mu }{{\lambda }_{1}})}{8\mu }.\end{aligned}$$

If ${e}_{1}^{A}> \frac{1+\frac{8\mu }{{\lambda }_{1}}-8\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{\Delta {\lambda }_{1}}>\frac{9-8\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{\Delta {\lambda }_{1}}$, $AW_{1}^{\prime } \left( {\varphi _{1} } \right) + WW^{\prime}_{1} \left( {\varphi _{1} } \right) > 0$.

Thus when $\frac{9-8\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{{\lambda }_{1}\Delta }{< e}_{1}^{A}<\frac{33-27\sqrt{2}{\sigma }_{1}\mu {e}^{\frac{{-\Delta }^{2}}{2{{\sigma }_{1}}^{2}}}}{{\lambda }_{1}\Delta }$, the function $AW_{1}^{\prime } \left( {\varphi _{1} } \right) + WW_{1}^{\prime } \left( {\varphi _{1} } \right)$ has a zero point, which means the function ${AW}_{1}({\varphi }_{1})+{WW}_{1}({\varphi }_{1})$ has a minimal point.

When ${s}_{1}>1$, we aggregate ${s}_{1}$ physicians into one physician whose service rate equal to $\mu {s}_{1}$. Thus, through the above proof, the same conclusion could be obtained.

As we have mentioned in Proposition 1 that ${AW}_{t}+{WW}_{t}$ is decreasing monotonically with ${s}_{t}$, thus when ${AW}_{\mathrm{t}}+{WW}_{\mathrm{t}}$ has a minimal point with respect to ${\varphi }_{t}$, ${s}_{t}$ has a minimal point. QED.