Introduction
Have you ever been in a supermarket checkout and wondered why you are in the slowest line? When I saw this article in Wired magazine, the title immediately grabbed my attention. As I read the article, I knew I had to try to put the theory to the test in SimEvents -- a discrete-event simulation add-on to Simulink that is suitable for modeling queuing systems.
Contents
The simple model
I started by building a simple model of a four register supermarket counter. I built two parallel versions of the setup -- one which uses four separate queues and one with a single "serpentine" queue that feeds all four registers. Let me explain the basic structure of the model.
open_system('basicModel');

I generate customers (entities in SimEvents) and assign them random service times at the register (represented by an attribute on the entity) based upon an exponential distribution. I also generate customers on the basis of an exponential arrival time. I chose the average service time (usually denoted by ) to be 2 mins and the average arrival time (usually denoted by
) to be 1 min.
I then clone the customers so that I can exercise the two different line configurations identically. All of this logic is inside the "Generate Customers" subsystem.
open_system('basicModel/Generate customers');

In the top half of the model, I use a Switch that routes customers to the shortest of the four Queues. Each Queue then feeds a Server representing a checkout register. This Server holds the customer for the amount of time that was setup during generation.
In the bottom half of the model, I use a single Queue that feeds the four registers via a Switch that routes customers to a free when one becomes available.
Finally, there are some Timer blocks (the ones with the Stop Watch icons) that keep track of timing statistics for the entities. This allows me to plot the average wait time of entities in the two configurations as well as the minimum and maximum wait times.
Let us run the model and see what gets plotted.
sim('basicModel'); open_system('basicModel/Average Wait Time')

As you call see, the configuration with the four queues on average results in longer wait times!
Well, the result was interesting and somewhat counterintiuitive as the original article predicted. Also, there are theoretical results for the average waiting time [1,2] for the one queue configuration which is matched almost exactly by the simulation. I implemented the theoretical waiting time in a MATLAB script called 'meanWaitTimeMMC'. Note that there are some approximations (but no analytical results) for the average arrival time of the 'shortest queue strategy' as well that I will not get into in this blog.
Wr = meanWaitTimeMMC(1/2,1,4);
fprintf('Mean waiting time = %f\n', Wr);
Mean waiting time = 2.173913
However, I wanted to dig deeper to see if I could spot patterns in which customers felt "unlucky" (or lucky) in the top configuration. I know you have been in lines silently begrudging the person who came after you who is now happily loading groceries into the shopping bags.
To do this, I needed to log more data. I setup individual IDs for each customer as I generated them (or as they arrived in the store). For each customer I then computed how many customers who arrived after them began service before them (a sort of "unluckiness index"). Similarly, for each customer I also computed how many customers (according to their IDs) began to get service after them (a sort of "luckiness index"). I know from experience that I only cared about the "unluckiness index" when I am at the store. Maybe that makes me a pessimist? To be complete, I also computed the difference between the two indices for each customer. While I was at it, I also spruced up the model with a number of the cool new model annotation features in R2014b.
Here is what the model and result look like.
open_system('twoTypesOfLines');

Let us now run this model. I have used the logged data of the IDs to plot the histograms of the luckiness and unluckiness indices for all customers.
sim('twoTypesOfLines');
plotLineTypeResults;
Customers who have at least one lucky person get served before them: 23% Customers who have at least one unlucky person get served after them: 43% Customers who have a net disadvantage: 20% Customers who have a net advantage: 38%

Let us make a few quick observations. Overall everyone loses (on average) in the multiqueue strategy because the average wait time is worse. However, if you are only in competition with your fellow shoppers, then only about 13% of customers see one person who came after them get served before them. That does not seem to match either the title of the article or what you experienced.
There is an explanation for this. The choice of average service time ( = 2 mins) and average arrival time (
= 1 min) results in a system utlization of 50% (
). If we now drop the customer arrival time to say 0.53 min, the utilization is 95% -- a scenario where you actually start to take notice of how many people are getting ahead of you. Let us now simulate the system with these parameters and look at the result.
set_param('twoTypesOfLines/Generate Customers/lambda_inv', 'meanexp', '0.53'); sim('twoTypesOfLines'); plotLineTypeResults;
Customers who have at least one lucky person get served before them: 55% Customers who have at least one unlucky person get served after them: 73% Customers who have a net disadvantage: 38% Customers who have a net advantage: 54%

Note now that 57% of customers now see at least one person who came after them get served before them. This is much more in line with what you might observe in that busy pre-Thanksgiving day supermarket line.
Of course we all also think we can beat the system. In the model I setup above, customers picked the line they joined by simply picking the shortest. Of course we all think we are smarter than that. We can estimate how long each customer can take in each line and then tailor our line choice not just by length but by an estimate of the "load" on each line.
To test this case, I modified the model to add more smarts to it. Using some of the statistics provided by SimEvents, I can now estimate the "load" of everyone waiting in a given Queue. This is done inside the "Read Queue States" block in the model below. There are some other modifications I will get into later.
open_system('advancedModel')

What is initially satisfying is that when you run the model, the two configurations now produce identical results as indicated by the coinciding lines on the Scope.
set_param('advancedModel/Generate Interruptions', 'Value', '0'); set_param('advancedModel/Control Estimate', 'Value', '0'); set_param('advancedModel/Generate Customers and Service Interruptions/lambda_inv', 'meanexp', '1'); sim('advancedModel') open_system('advancedModel/Logging/Average Wait Times')

However, in reality we are not good at making estimates of how long people will take in a register. If we assume underestimation (I sugest googling optimism bias) of between 0 and 75% for the load for each person waiting in the queue, we will once again begin to see a discrepancy between the two schemes with the one queue system proving the better choice.
set_param('advancedModel/Control Estimate', 'Value', '1'); sim('advancedModel') open_system('advancedModel/Logging/Average Wait Times')

Now, we can get even closer to reality. We know that as customers are served, sometimes you need the dreaded "price check" to happen once in a while. This is modeled by generating random service interruptions. SimEvents Server blocks offer this nice control mechanism to suspend service using a signal port. I have chosen to do this as well in the generation block in my model. Throw in the imperfect estimations as well and we are back to a world where the one queue system is a better option.
set_param('advancedModel/Generate Interruptions', 'Value', '1'); sim('advancedModel') open_system('advancedModel/Logging/Average Wait Times')

Once again we can play with and
to see how utilization plays a role in how lucky a customer feels. Choosing average service time (
= 2 mins) and average arrival time (
= 0.7 min) results in a utilization of 71% without interruptions. This also results in 55% of customers now having at least one person who came after them get served before them.
set_param('advancedModel/Generate Customers and Service Interruptions/lambda_inv', 'meanexp', '0.7'); sim('advancedModel'); plotLineTypeResults;
Customers who have at least one lucky person get served before them: 55% Customers who have at least one unlucky person get served after them: 63% Customers who have a net disadvantage: 40% Customers who have a net advantage: 52%


Conclusion
I enjoyed taking the theory behind the original article and building a model I could simulate. The model also helped me look at some questions that were always in the back of my mind as I stood in those long checkout lines. I will place the models and scripts in a Simulink project on MATLAB central so that you can download and try it yourself. Feedback welcome!
References
[1] Gautam, Natarajan (2012). Analysis of Queues: Methods and Applications. CRC Press. ISBN 9781439806586. [2] Wikipedia: M/M/c