<jupyter_start><jupyter_text>Single Channel Queuing Example  Systems Engineering and Analysis     Prepared by: Dr. <NAME>  The M/M/1 Queue SystemThe M/M/1 Queue is named using the Kendall Notation, which denotes the arrival process, the service time distribution, and the number of service units. The M stands for Markovian (or Memoryless). It uses an exponential distribution to describe the time spacing of the arrival or servicing events. https://en.wikipedia.org/wiki/Kendall%27s_notationThe SYSE530 course text describes this system in Section 10.3 as Single-Channel Queuing Model.https://en.wikipedia.org/wiki/M/M/1_queueSome additional reading for fun:https://towardsdatascience.com/the-poisson-distribution-and-poisson-process-explained-4e2cb17d459<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson, expon, norm<jupyter_output><empty_output><jupyter_text>Number of Arrivals per Time PeriodIn the following block, we'll simulate a Poisson distribution to describe the number of units arriving to a queuing system. We'll be using the documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.poisson.htmlFor additional context, see https://en.wikipedia.org/wiki/Poisson_point_process.From the example on Page 299 of B&F, suppose a queue experiences Poisson arrives with a mean rate of 1/10 units per time period. This means $$ \lambda = 10 $$in the book's notation.<jupyter_code># Simulate many samples
N=100000

# Define the mean arrival rate
arrival_rate_mean = 10

# Specify  the number of arrivals per period
# use lam instead of lambda because lambda is a protected word in Python
lam = 1/arrival_rate_mean

print("average number of arrivals per time = {}".format(arrival_rate_mean))

#Generate random values from the poisson distribution
arrival_counts = poisson.rvs(arrival_rate_mean,size=N)

#Draw a normalized histogram and compare it to the theoretical distribution
#place the bins in the center of the discrete values for the samples
bins = np.arange(0, max(arrival_counts) + 1.5) - 0.5

#Show the emperical distribution of the generated random values
plt.hist(arrival_counts,
         bins=bins,density=True,color="yellow",
         ec="black",label="Simulated PMF");

#Generate a range for the theoretical mass function
x = np.arange(poisson.ppf(0.0001, arrival_rate_mean),
              poisson.ppf(0.9999, arrival_rate_mean))
#Plot the theoretical PMF
plt.stem(x, poisson.pmf(x, arrival_rate_mean),   label='Poisson PMF')
plt.legend()
plt.grid()
plt.title("Poisson Distribution, $\lambda = {:0.2f}$".format(lam));
plt.xlabel("Number of arrivals")
#compare simulation mean 
arrival_counts.mean()<jupyter_output><empty_output><jupyter_text>Time between each arrivalThe space between each arrival in Poisson process is distributed as an exponential distribution.https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.htmlThe scale factor for the distribution is $1/\lambda $ which equals the arrival rate.We'll simulate the arrival spacing for the queue.<jupyter_code>#Show the emperical distribution of the generated random values
#place the bins in the center of the discrete values for the samples
arrival_intervals = expon.rvs(scale=arrival_rate_mean,size=N)
plt.hist(arrival_intervals,bins=100,density=True,color="cyan",ec="black",label="Simulated PDF");
#Generate a range for the theoretical mass function
x = np.arange(expon.ppf(0.0001, scale=arrival_rate_mean),
              expon.ppf(0.9999, scale=arrival_rate_mean),0.01)
#Plot the theoretical PMF
plt.plot(x, expon.pdf(x, scale=arrival_rate_mean),  'r-', label='Exponential PDF')
plt.legend()
plt.grid()
plt.xlabel("Time Spacing for a Poisson Process")
plt.ylabel("Probability Density Function");
# Determine theoretical moments for the distribution
mean, var, skew, kurt = expon.stats(scale=arrival_rate_mean,moments='mvsk')
print(mean, var, skew, kurt)
#Compare to the simulated distribution
arrival_intervals.mean()
#Theoretical mean
1/lam
#Compare to the simulated distribution variance
arrival_intervals.var()<jupyter_output><empty_output><jupyter_text>Simulate the distribution describing time for servicingThis is also an exponential distribution<jupyter_code>#Generate random values from the exponential distribution
service_time_mean = 4
# Determine the service rate
mu = 1/service_time_mean #Book notation for service completions per period

service_times = expon.rvs(scale=service_time_mean,size=N)

#place the bins in the center of the discrete values for the samples
bins = 100

#Show the emperical distribution of the generated random values
plt.hist(service_times,bins=bins,density=True,color="cyan",
         ec="black",label="Simulated PDF");

#Generate a range for the theoretical mass function
x = np.arange(expon.ppf(0.0001, scale=service_time_mean),
              expon.ppf(0.9999, scale=service_time_mean),0.01)
#Plot the theoretical PMF
plt.plot(x, expon.pdf(x, scale=service_time_mean),
         'r-', label='Exponential PDF')
plt.legend()
plt.grid()
mean, var, skew, kurt = expon.stats(scale=service_time_mean,moments='mvsk')
print(mean, var, skew, kurt)
#theroretical expected value (mean)
1/mu
#Compare to the simulated distribution
np.array(service_times).mean()
#Theoretical variance
1/mu**2 
#Compare to the simulated distribution
np.array(service_times).var()<jupyter_output><empty_output><jupyter_text>Will there be convergence?This is the first check for a queue. Will the queue be able to handle the in-flux. The processing has to be faster than the arrivals or else the queue will back up. See the note after Eqn. 10.9 in B&F.<jupyter_code># check the ratio to see if it is less than 1 to get convergence
print(lam/mu)
if lam/mu < 1:
    print("Convergence Criteria Achieved.")
else:
    print("Queue will not converge. The arrivals are faster than the servicing")<jupyter_output>0.4
Convergence Criteria Achieved.<jupyter_text>Simulate a Queuing SystemWe'll setup a routine that loops through time and processes "customers." The customer arrival time is separated by a random time that is exponentially distributed. The processing time will be determined by a random sample from another exponential distribution. With both of these distributions as exponential, we can determine theoretical properties for the M/M/1 system. The customer has it's own needs, so we'll define it's servicing time when the customer is created (i.e. arrives to the queue).We'll use an object-oriented approach and create a class the defines the customer and some of it's properties and functions.<jupyter_code>#Create an object that can keep track of itself through the system
# Pass in an identifier and the global arrival time. 

class Customer():
    def __init__(self,id,arrival):
        self.id = id
        self.wait_time = 0
        self.arrival_time = arrival
        #Use a sample from a random variable for the service time
        self.service_time = expon.rvs(scale=service_time_mean) 
        #self.service_time = abs(norm.rvs(loc=4,scale=1) )
        # This is a minimum. It will be updated upon ejection from the queue system
        self.finished_time = self.arrival_time + self.service_time + self.wait_time
        
    def get_queue_time(self):
        self.queue_time = self.finished_time - self.arrival_time
        return self.queue_time<jupyter_output><empty_output><jupyter_text>Simulate the QueueThis is a Monte Carlo simulation of the Single Channel Queuing Model as described in Section 10.3 of B&F.<jupyter_code># This is for printout of the queue simulation
# Set this to zero to suppress output
debug_periods = 0

# Initialize some lists to keep track of the queue performance.
finished_customers = []
customers_in_line = []
customers_in_system = []
customer_wait_time = []
customer_system_time = []


# Specify the number of time periods to run
time_periods = 1000000
# Iterate to find convergence.
#for time_periods in [100000,200000,400000]:


# Make an ordered list the waiting line of the queue system
# We can add elements uing the append method
# We can remove elements using the pop method
waiting_line = []


system_size=0
arrivals = 0
arrival_time = 0
# This creates an instance of the customer
# We pass the id of the customer into the constructor
# The first time is zero (second argument)
# If the service_cust is defined, then it is in service. 
service_cust = Customer(arrivals,arrival_time)
# Always add one when the customer is created
system_size += 1
# Always update the arrival time of the next customer with a sample from the 
# exponential distribution
arrival_time += expon.rvs(scale=arrival_rate_mean)
times = []

delta_t=0.001
total_time = time_periods*delta_t
print("total time:", total_time)
t=0
while t < total_time:
    # Calculate the current time as t
    t += delta_t
    times.append(t)
    # Keep track of how many customers are in the system
    customers_in_system.append(system_size)
    # Keep track of customers in line
    customers_in_line.append(len(waiting_line))

    # Check for new arrivals
    if t >= arrival_time:
        arrivals += 1
        # Create a customer to simulate an arrival. 
        # The customer class generates its own servicing time.
        cust = Customer(arrivals,arrival_time)
        # Always increase the number of customers in the system
        system_size += 1

        # update the arrival time with the next arrival
        arrival_time += expon.rvs(scale=arrival_rate_mean)

        if t < debug_periods:
            print("t = {}, waiting_line length = {}".format(t,len(waiting_line)))
            print("Created Customer {} at time {} with service time of {}".format(cust.id, cust.arrival_time, cust.service_time))

        # Put the newly created customer in the wait line
        waiting_line.append(cust)

    # Check to see if a customer is finished
    #Be sure there is a customer in service
    if service_cust is not None:
        # Since the customer is in service, we should see if it's done.
        if t >= service_cust.finished_time:
            if t < debug_periods:
                print("Customer {} finished at time {}".format(service_cust.id, t))
            # Keep track of customers coming through the queue
            finished_customers.append(service_cust)
            customer_system_time.append(service_cust.get_queue_time())
            # Remove the customer from the system.
            system_size -= 1
            # Let the program know there is no customer being serviced.
            service_cust = None
    else:
        #There's no customer in service, so only the ones in the line are in the system
        try:
            # Get a customer from the line
            service_cust = waiting_line.pop(0)
        except IndexError:
            # There wasn't anyone in line, so move to the next time.
            continue
        # determine how long they were in the wait line    
        service_cust.wait_time = t - service_cust.arrival_time
        # Keep track of the wait time.
        customer_wait_time.append(service_cust.wait_time)
        # Determine when the finish should be
        service_cust.finished_time = t + service_cust.service_time

        if t < debug_periods:
            print("Retrieved Customer {} from the line at time {:0.3f} after waiting {:0.3f}.".format(service_cust.id, t, service_cust.wait_time))
            print("Anticipated finish for {} is {:0.3F}".format(service_cust.id,service_cust.finished_time))
N=time_periods//2
print("np.mean(customers_in_system)",np.mean(customers_in_system[N:]))<jupyter_output>total time: 1000.0
np.mean(customers_in_system) 1.1186637626724747<jupyter_text>Note: Convergence is slow for this system. Orders of magnitude more samples are required for convergence  Average number of customers waiting in line<jupyter_code>average_number_in_line = np.mean(customers_in_line[N:])
#times = np.arange(N)*delta_t
plt.plot(times[N:],customers_in_line[N:])
plt.xlabel("Time (periods)")
plt.ylabel("Number of customers in line")
plt.title("Average number of customers in line: {:0.3f}".format(average_number_in_line));
# Theoretical average length of line
# See Equation 10.11
L_q = lam**2/(mu*(mu-lam))
L_q<jupyter_output><empty_output><jupyter_text>Average waiting time for customers in line<jupyter_code># Average simulated time in line
np.mean(customer_wait_time)
# Theoretical mean time an arrival waits for service
# Equation 10.18
W_q = lam/(mu*(mu-lam))
W_q<jupyter_output><empty_output><jupyter_text>Average time spent in the queuing system<jupyter_code># Average simulated time in the whole system
np.mean(customer_system_time)
# Theoretical average customer time through the queue.
# Equation 10.19
t_m = 1/(mu-lam)
t_m<jupyter_output><empty_output><jupyter_text>Average number of customers in the systemThis includes both customers in line and customers being serviced.<jupyter_code># Incorporate less samples to see the operation of the queue
# Also include some measures
average_number_in_system = np.mean(customers_in_system)
#Reduce the samples for faster plotting
plt.plot(times[:N],customers_in_system[:N])
plt.xlabel("Time (periods)")
plt.ylabel("Number of customers in system")
plt.title("Average number of customers in system: {:0.3f}".format(average_number_in_system));
# Mean theoretical length or elements in the system
# Equation 10.10
L = lam/(mu-lam)
L
# Probability of size of the line at any time (steady state)
#Eq 10.9
def Pn(n,l=lam,mu=mu):
    assert type(n)==int
    return (1-l/mu)*(l/mu)**n

pmf=[]
x = []
expected_val = 0
for n in range(15):
    x.append(n)
    pmf.append(Pn(n))
    print(n,Pn(n))
    
    expected_val += n*Pn(n)
#Expected Value Verification
print(expected_val,"=?=",L)
# This cell takes a while with large samples
#place the bins in the center of the discrete values for the samples
bins = np.arange(0, max(customers_in_system[N:]) + 1.5) - 0.5
vals,bins,_ = plt.hist(customers_in_system[N:], bins=bins, density=True,label="Simulated")
plt.plot(x,pmf,'d',label="Theoretical")
plt.xlabel("Number of Customers in Line")
plt.ylabel("Probability Mass Function")
plt.title("Compare to Figure 10.6")
plt.legend();
# Compute the expected value of the simulated distribution 
expected_sum = 0
for v,b in zip(vals,bins):
    expected_sum += v*(b+.5)
#this should closely match the average_number_in_system
expected_sum
# This is the same as before
average_number_in_system<jupyter_output><empty_output><jupyter_text>Minimum Cost Service RateThe expected cost of the queuing system is the sum of the expected waiting period and the expected facility cost per period.$$ TC_m = WC_m + FC_m$$where$$ WC_m = C_w\times n_m = \frac{C_w\lambda}{\mu-\lambda}$$and $$ FC_m = C_f\times\mu $$The costs are $C_w$ that represents the cost of waiting during the period and $C_f$ is the cost of servicing.To optimize, take the derivative of Total Cost with respect to $\mu$ and set to zero. Solving gives$$ \mu = \lambda +\sqrt{\frac{\lambda C_w}{C_f}}$$<jupyter_code>from math import inf
#Example in Section 10.3.5
lamb = 1/8
Cw = 0.10
Cf = 0.165

#Construct Table 10.1
print("mu   WC_m   FC_m   TC_m")
for muu in [0.125, 0.150, 0.200, 0.250, 0.300, 0.400, 0.500, 0.600, 0.800, 1.000]:
    FC_m = Cf*muu
    try:
        WC_m = Cw*lamb/(muu-lamb)
        TC_m = WC_m + FC_m
    except ZeroDivisionError:
        WC_m = inf
        TC_m = inf
    print("{:0.2f}  {:0.2f}  {:0.2f}  {:0.3f}".format(muu, WC_m, FC_m, TC_m))
#Optimal Cost Using Equation 10.24
from math import sqrt
muu = lamb + sqrt(lamb*Cw/Cf)
muu<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08953593006813215, 'fraction_numerical': 0.017997171872991387, 'mean_word_length': 4.348573392918529, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 50, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '343175', 'license': 'MIT', 'n_tokens_mistral': 5099, 'n_tokens_neox': 4714, 'n_words': 1640, 'path': 'Single Channel Queuing Example.ipynb', 'repo_name': 'Dr-Daily/SystemsEngineering', 'size': '15564'}