I am working on a problem of scheduling appointments for a stochastic, single-server queue. There are $n$ customers who each have independent, randomly distributed service durations $Z_i$. The customers are assigned appointment times $x_i$ where $x_i<x_j$ for $i<j$. An appointment schedule is some $\mathbf{x}=\langle x_1,\ldots,x_N\rangle$ where $x_i>0$.
Customer 1 arrives at time $x_i$ and is immediately served; the duration of his service is unknown until he is actually served and represented by $Z_1$. Customer 2 arrives at time $x_2>x_1$ and is served as soon as is available (after Customer 1 has been served); his service takes $Z_2$. Customer 3 arrives at $x_3>x_2$, etc.
The waiting time for customer $i$ is given by $W_i=\max\{X_{i-1}+W_{i-1}+Z_{i-1}-X_i,0\}$ and the server idle time immediate prior to Customer $i$ being served is $S_i=\max\{-X_{i-1}-W_{i-1}-Z_{i-1}+X_i,0\}$.
In my application, I need to approximate $\text{E}(W_i)$ and $\text{E}(S_i)$ for many (millions) of different schedules. The service time distributions ($Z_i$) are the same for every schedule.
Currently, I'm doing a naive Monte Carlo simulation: sample scenarios ($\mathbf{z}=\langle z_1,\ldots,z_N\rangle$) from the service distributions, compute the waiting time and idle times for each scenario, and average the results. However, this method has proven to be slower than I'd like. I can increase the speed by reducing the number of scenarios and increasing the number of estimates, but I haven't been able to find the right balance.
I'm looking for suggestions on speeding up this computation.
