{
"metadata": {
"name": "",
"signature": "sha256:c38b3f543fbb23a3800bcfbdec470195c6e5ffd97e7b79df3327e8fef7f70c8d"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"butcher5"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Markov Butcher Shop sells two types of meat: beef and pork. Customers arrive at the butcher shop and form a single queue. There are two butchers who serve customers from the queue on a first-come-first-served basis.\n",
"\n",
"Customer demand is as follows: beef only, 50%; pork only, 30%; both beef and pork, 20%. Customers who want only one type of meat place 1 order; those who want both types of meat must place 2 separate orders.\n",
"\n",
"Based on historical data, the interarrival time between customers is uniformly distributed between 1.5 and 3.5 minutes. The time for a butcher to fill a single order is uniformly distributed between 2.5 and 6.5 minutes. The interarrival times and service times are assumed to be independent."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Simulate the shop for 2 hours. The butcher shop is interested in the following performance measures:\n",
"\n",
"* the average time a customer spends in the queue (i.e. average delay),\n",
"* the time average number of customers in the queue,\n",
"* the time average number of busy butchers.\n",
"\n",
"Determine a point estimate (i.e. observed sample mean) and an interval estimate (i.e. 95% confidence interval) for each performance measure using 50 replications."
]
},
{
"cell_type": "heading",
"level": 4,
"metadata": {},
"source": [
"Setup"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"##### Setup #####\n",
"# Friendly floating-point division\n",
"from __future__ import division\n",
"\n",
"# Import all simulation functions from SimPy\n",
"from SimPy.Simulation import *\n",
"\n",
"# Import bisect_left from bisect\n",
"from bisect import bisect_left\n",
"\n",
"# Import various functions from NumPy\n",
"from numpy import mean, std, sqrt, inf\n",
"\n",
"# Import seed initializer and random sampling functions from NumPy\n",
"from numpy.random import seed, rand, uniform\n",
"\n",
"# Import t random variable from SciPy\n",
"from scipy.stats import t"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "heading",
"level": 4,
"metadata": {},
"source": [
"Model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below is the SimPy code from the **butcher5** homework problem, with the appropriate modifications."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"##### Distributions #####\n",
"# Customer type distribution: \n",
"# Beef = 0, Pork = 1, Both = 2\n",
"def customerTypeDist():\n",
" # Sorted list of values that X can take, including -inf\n",
" a = [-inf, 0, 1, 2]\n",
"\n",
" # cdf at each value that X can take\n",
" cdf = [0, 0.50, 0.80, 1]\n",
"\n",
" # Generate variate\n",
" variate = a[bisect_left(cdf, rand())]\n",
"\n",
" # Return generated variate\n",
" return variate\n",
" \n",
"\n",
"##### Parameters #####\n",
"class P:\n",
" # Customers arrive at the entrance with uniformly distributed\n",
" # interarrival times between 1.5 and 3.5 minutes\n",
" interarrivalTimeMin = 1.5\n",
" interarrivalTimeMax = 3.5\n",
" \n",
" # Service times are uniformly distributed between 2.5 and 6.5 minutes\n",
" serviceTimeMin = 2.5\n",
" serviceTimeMax = 6.5\n",
" \n",
" # Number of servers\n",
" nServers = 2\n",
" \n",
" # Shop is open for 2 continuous hours\n",
" simulationTimeMax = 2 * 60\n",
" \n",
"\n",
"##### Processes #####\n",
"# Beef only customer\n",
"class CustomerBeef(Process):\n",
" def behavior(self):\n",
" # Customer arrives, joins queue\n",
" arrivalTime = now()\n",
" yield request, self, R.server\n",
" \n",
" # Customer is released from queue and starts service\n",
" M.delay.observe(now() - arrivalTime)\n",
" serviceTime = uniform(low = P.serviceTimeMin, high = P.serviceTimeMax)\n",
" yield hold, self, serviceTime\n",
" \n",
" # Customer finishes service, leaves\n",
" yield release, self, R.server\n",
" \n",
"# Pork only customer\n",
"class CustomerPork(Process):\n",
" def behavior(self):\n",
" # Customer arrives, joins queue\n",
" arrivalTime = now()\n",
" yield request, self, R.server\n",
" \n",
" # Customer is released from queue and starts service\n",
" M.delay.observe(now() - arrivalTime)\n",
" serviceTime = uniform(low = P.serviceTimeMin, high = P.serviceTimeMax)\n",
" yield hold, self, serviceTime\n",
" \n",
" # Customer finishes service, leaves\n",
" yield release, self, R.server\n",
"\n",
"# Beef and pork customer\n",
"class CustomerBoth(Process):\n",
" def behavior(self):\n",
" # Customer arrives, joins queue\n",
" arrivalTime = now()\n",
" yield request, self, R.server\n",
" \n",
" # Customer is released from queue and starts service\n",
" M.delay.observe(now() - arrivalTime)\n",
" serviceTime1 = uniform(low = P.serviceTimeMin, high = P.serviceTimeMax)\n",
" serviceTime2 = uniform(low = P.serviceTimeMin, high = P.serviceTimeMax)\n",
" yield hold, self, serviceTime1 + serviceTime2\n",
" \n",
" # Customer finishes service, leaves\n",
" yield release, self, R.server\n",
" \n",
"# Entrance\n",
"class Entrance(Process):\n",
" def behavior(self):\n",
" # At the start of the simulation, no customers have arrived\n",
" nCustomers = 0\n",
" \n",
" # Customer arrivals\n",
" while True:\n",
" # Wait until the next arrival\n",
" interarrivalTime = uniform(low = P.interarrivalTimeMin, high = P.interarrivalTimeMax)\n",
" yield hold, self, interarrivalTime\n",
" \n",
" # Customer type?\n",
" customerType = customerTypeDist()\n",
" \n",
" # Create a new customer, depending on this passenger's type\n",
"\n",
" if customerType == 0:\n",
" c = CustomerBeef(name=\"Customer {0} (Beef)\".format(nCustomers))\n",
" elif customerType == 1:\n",
" c = CustomerPork(name=\"Customer {0} (Pork)\".format(nCustomers))\n",
" elif customerType == 2:\n",
" c = CustomerBoth(name=\"Customer {0} (Both)\".format(nCustomers))\n",
"\n",
" # Count this new passenger\n",
" nCustomers += 1\n",
" \n",
" # Activate the customer's behavior\n",
" activate(c, c.behavior())\n",
"\n",
" # Count this new customer\n",
" nCustomers += 1\n",
"\n",
"\n",
"##### Resources #####\n",
"class R:\n",
" # Server\n",
" server = None\n",
" \n",
"\n",
"##### Monitors #####\n",
"class M:\n",
" # Queue waiting time\n",
" delay = None\n",
" \n",
"\n",
"##### Simulation #####\n",
"def model(inputSeed):\n",
" # Initialize SimPy \n",
" initialize()\n",
"\n",
" # Initialize a seed for the random number generator (more on this later in the semester)\n",
" seed(inputSeed)\n",
"\n",
" # Create the server resource\n",
" R.server = Resource(capacity=P.nServers, monitored=True)\n",
" \n",
" # Create the delay monitor\n",
" M.delay = Monitor()\n",
"\n",
" # Activate the entrance (to generate customers)\n",
" e = Entrance()\n",
" activate(e, e.behavior())\n",
" \n",
" # Run the simulation\n",
" simulate(until=P.simulationTimeMax)\n",
" \n",
" # Compute average delay\n",
" avgDelay = M.delay.mean()\n",
" \n",
" # Compute time average number of customers in queue\n",
" tAvgCustQ = R.server.waitMon.timeAverage()\n",
" \n",
" # Compute time average number of busy butchers\n",
" tAvgBusyS = R.server.actMon.timeAverage()\n",
" \n",
" # Return a list of performance measures\n",
" return [avgDelay, tAvgCustQ, tAvgBusyS]"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "heading",
"level": 4,
"metadata": {},
"source": [
"Experiment"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Replicate simulation 50 times\n",
"# Create a list of observations for each performance measure\n",
"n = 50\n",
"avgDelayObs = []\n",
"tAvgCustQObs = []\n",
"tAvgBusySObs = []\n",
"for i in range(n):\n",
" [avgDelay, tAvgCustQ, tAvgBusyS] = model(i*123)\n",
" avgDelayObs.append(avgDelay)\n",
" tAvgCustQObs.append(tAvgCustQ)\n",
" tAvgBusySObs.append(tAvgBusyS)\n",
"\n",
" \n",
"#### Average delay ####\n",
"# Observed sample mean\n",
"avgDelaySM = mean(avgDelayObs)\n",
"\n",
"# Observed sample standard deviation\n",
"avgDelaySSD = std(avgDelayObs, ddof = 1)\n",
"\n",
"# Confidence level 0.05\n",
"alpha = 0.05\n",
"avgDelayCIL = avgDelaySM - t.ppf(1 - alpha/2, n - 1) * avgDelaySSD / sqrt(n)\n",
"avgDelayCIR = avgDelaySM + t.ppf(1 - alpha/2, n - 1) * avgDelaySSD / sqrt(n)\n",
"\n",
"print(\"Average delay:\")\n",
"print(\" Sample mean: {0}\".format(avgDelaySM))\n",
"print(\" Sample standard deviation: {0}\".format(avgDelaySSD))\n",
"print(\" {0}% confidence interval: [{1}, {2}]\".format((1 - alpha)*100, avgDelayCIL, avgDelayCIR))\n",
"\n",
"\n",
"#### Time average number of customers in queue ####\n",
"# Observed sample mean\n",
"tAvgCustQSM = mean(tAvgCustQObs)\n",
"\n",
"# Observed sample standard deviation\n",
"tAvgCustQSSD = std(tAvgCustQObs, ddof = 1)\n",
"\n",
"# Confidence level 0.05\n",
"alpha = 0.05\n",
"tAvgCustQCIL = tAvgCustQSM - t.ppf(1 - alpha/2, n - 1) * tAvgCustQSSD / sqrt(n)\n",
"tAvgCustQCIR = tAvgCustQSM + t.ppf(1 - alpha/2, n - 1) * tAvgCustQSSD / sqrt(n)\n",
"\n",
"print(\"Time average number of customers in queue:\")\n",
"print(\" Sample mean: {0}\".format(tAvgCustQSM))\n",
"print(\" Sample standard deviation: {0}\".format(tAvgCustQSSD))\n",
"print(\" {0}% confidence interval: [{1}, {2}]\".format((1 - alpha)*100, tAvgCustQCIL, tAvgCustQCIR))\n",
"\n",
"\n",
"#### Time average number of busy servers ####\n",
"# Observed sample mean\n",
"tAvgBusySSM = mean(tAvgBusySObs)\n",
"\n",
"# Observed sample standard deviation\n",
"tAvgBusySSSD = std(tAvgBusySObs, ddof = 1)\n",
"\n",
"# Confidence level 0.05\n",
"alpha = 0.05\n",
"tAvgBusySCIL = tAvgBusySSM - t.ppf(1 - alpha/2, n - 1) * tAvgBusySSSD / sqrt(n)\n",
"tAvgBusySCIR = tAvgBusySSM + t.ppf(1 - alpha/2, n - 1) * tAvgBusySSSD / sqrt(n)\n",
"\n",
"print(\"Time average number of busy servers:\")\n",
"print(\" Sample mean: {0}\".format(tAvgBusySSM))\n",
"print(\" Sample standard deviation: {0}\".format(tAvgBusySSSD))\n",
"print(\" {0}% confidence interval: [{1}, {2}]\".format((1 - alpha)*100, tAvgBusySCIL, tAvgBusySCIR))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Average delay:\n",
" Sample mean: 5.42356500318\n",
" Sample standard deviation: 2.92100139557\n",
" 95.0% confidence interval: [4.59342559275, 6.25370441361]\n",
"Time average number of customers in queue:\n",
" Sample mean: 2.18448576605\n",
" Sample standard deviation: 1.15864824152\n",
" 95.0% confidence interval: [1.85520157961, 2.51376995248]\n",
"Time average number of busy servers:\n",
" Sample mean: 1.90014129843\n",
" Sample standard deviation: 0.0498621966512\n",
" 95.0% confidence interval: [1.88597061895, 1.91431197791]\n"
]
}
],
"prompt_number": 3
}
],
"metadata": {}
}
]
}