UQ MATH2504
Programming of Simulation, Analysis, and Learning Systems
(Semester 2 2022)
This is an OLDER SEMESTER.
Go to current semester
Project #2 - Semester 2, 2022
Discrete Event Simulation (+ perspective seminar summary)
(last edit: October 5, 2022)
Here is an explainer video for the project.
Due: End of Day, Friday, 21, October, 2022. Late submissions are not accepted. Work in pairs - you must have a single partner for submission. Any deviation from this requires special permission from the course coordinator.
Note: The teaching staff will only answer questions (via Ed, consultation hour, or practicals) regarding this assignment up to the late evening of Wednesday October 18..
Weights and marking criteria: Total number of points: 100. There are 10 points for handing in according to the hand-in instructions, including a voice recording, neat output, and very importantly the GitHub repo. There are 10 points for the summary of Anna Foeglein's perspective seminar (this is not a bonus as in previous assessment). The remaining 80 points are for the project itself.
Submission format: This assignment should be submitted via a GitHub Repo and a PDF file via Blackboard.
Specific instructions for the GitHub repo are below. It is important that the GitHub repo be made private and the course user name uqMATH2504 be invited as a collaborator. It is also important to name the repo exactly as <<FIRST NAME1>>-<<LAST NAME1>>__<<FIRST NAME2>>-<<LAST NAME2>>2504-2022-PROJECT2. So for example if your names are "Queen Latifah" and "Sacha Baron Cohen." the repo name should be Queen-Latifah__Sacha-Baron-Cohen-2504-2022-PROJECT2.
The PDF file should be a nice formatted file that has:
Your name, student number, and assignment title (Project 2 - 2022) on the top.
A (clickable) link to your GitHub repo.
Instructions on how to run the main simulation code in the repo via an easy to use
simulation_script.jlfile.All of the project code formatted in an organized manner.
Output and answers to specific questions formatted in an organized manner.
Plots and graphs should have clean organized labels, legends, etc.
The recommended way to make the PDF file is via a Jupyter notebook where you copy in some code and output into the notebook, and in certain cases use include() to run Julia code if appropriate. You also comment on questions in this PDF (e.g. when asked to answer things not via code). If desired you could keep this jupyter notebook (.ipynb file) in the repo. However, this Jupyter notebook will not be checked (only the PDF file), and it isn't required to be a "runnable" notebook. In any case, please avoid pixilated screenshots of code, so for example if you choose to format your PDF file in MS-word (instead of Jupyter), make sure the code is clean, formatted, and readable.
Marking responses will be made by the teaching staff by annotating selected parts of your PDF file via blackboard. The teaching staff will also inspect your GitHub repo and potentially run simulation_script.jl. A very readable and clean PDF file is important. In printing Jupyter to PDF (or exporting to PDF) there may sometimes be excessive white space. Do not worry about this extra white space; it is not a problem.
Work in pairs but discuss with others: This is work in pairs. Plagiarism between groups will not be accepted. Nevertheless, feel free to consult with friends or classmates via Ed and other means about how to go about certain tasks. Answering Ed questions of others is recommended.
Marking Criteria:
10 points are allocated for following instructions and the GitHub repo.
10 points are allocated to the summary of the perspective seminar. 2 points for coherent and grammatically correct writing. 6 points for completeness (the writeup should answer all questions posed). 2 points of originality and style.
80 points are for the simulation project itself.
In general, points will be deducted for sloppy coding style. Make sure to have your code properly indented, to use sensible and consistent variable names, and to write code that is in general clean and consistent.
Setting up your GitHub repo for hand-in (10pts)
Ideally use the same accounts you used for PROJECT1 and BigHW. Choose one of the team members to "own" the repo and the other(s) will be invited as collaborators.
Make sure the repo is private.
Invite the course GitHub user,
uqMATH2504as a collaborator. You may do so early on while working on the project, and must do this no later than the project due date.Do not make any changes (commits) to the repo after the project due date.
Create a local clone of the repo. It is recommended that use use
gitcommand line via the shell to work on making changes/additions to the assignment and submitting the changes. However you are free to use any other mechanism (VS-Code, GitHub desktop, etc).If for some reason you are not able or willing ot use GitHub an alternative is GitLab. This is not recommended as it adds additional work to the teaching staff, however if you wish to use GitLab instead of GitHub contact the teaching staff for permission.
Your GitHub repo should be formatted as exactly as follows:
Have a
README.mdfile which states your names, the assignment title, and has a (clickable) link to the assignment instructions on the course website (to this document).Have a
LICENSE.mdfile. Choose a license as you wish (for example the MIT license). However keep in mind that you must keep the submission private until the end of the semester.Have a
.gitignorefile.Have basic running instructions on how to run the code. This can be similar to what was provided for Project 1 (polynomials).
Have
srcandtestfolders similarly to Project 1.Note: make sure that there aren't any files in your submission repo except for those mentioned above (with the exception of perhaps a Jupyter
.ipynbfile if you choose to use it for creating the PDF). Use may use the.gitignorefile to ensuregitdoes not commit additional files and output files to yourrepo.All code is to be part of a Julia module (with the module statement). Call the module
GeneralizedJacksonSim.Make sure the project has a
Manifest.tomland aProject.tomlfile that specify all the dependencies exactly. Do not have extra dependencies that are not used.Have a
simulation_script.jlfile which runs all of the code needed for the tasks. As this script runs, make sure it prints some light intermediate output (e.g. progress), especially for simulations that take significant time to run.
Since this project involves work in pairs, it is recommended that you try to use GitHub effectively. You may work on parallel branches with frequent merges (but this is not required), or work on the same branch via peer programming. The way you develop is not being assessed, still it is recommended you make good use of GitHub.
Perspective seminar question (10pts)
Write a 1-4 paragraph summary of Anna Foeglein's perspective seminar. The writeup should answer these questions without treating the questions directly. That is, don't write "the answer to Q1 is..." etc. Instead incorporate the answers in a flowing writeup report (like a blog post or letter). You may use first person, or similar. Present the writeup as part of your PDF file submission (as the first item in the submission).
The writeup should address the following questions:
What was Anna's career path, and her relationship with software and mathematics?
Where does she currently work? What does she do in her job? How do you perceive her personal experience from her job?
Discuss a few key tools that she spoke about during the talk, such as programming languages, mathematics, machine learning, statistical tools, simulation platforms, and the like. Was this the first time that you heard about any of these tools?
Have you picked up any tips from Anna's talk, which you can perhaps use in your career down the road?
Feel free to state anything else that you found interesting from the talk and to contrast this talk with the previous talk.
Note: Since the submission is in pairs you should only submit a single writeup for the pair/group.
Overview of the model
Your key task is to create a discrete event simulation engine for stochastic models that are called "Open Generalized Jackson Networks". Jackson networks are well studied in classical "applied probability" and "stochastic modelling". The phrase "generalized" in "Generalized Jackson Networks", is also well studied, but is slightly less popular. "Generalized" networks are less amenable to closed form mathematical results (explicit formulas) - and this is why we will use discrete event simulation to analyze them. The "Open" phrase implies that customers may enter and leave these networks (in contrast with "closed" networks that have a finite number of customers). Note: you do not need to dive deeply into the Wikipedia entries or other literature. All the information you need about the model is supplied in this project description.
A Jackson network (or "open generalized Jackson network") is a network of queues (also known as "nodes" or "stations") where customers (also known as "jobs") arrive from the outside world, are served at the queues in a "first come first served manner" and then move to other nodes or leave. As a very rough motivational example consider an amusement park (e.g. Dream World). There are multiple rides, and individuals move throughout the park and queue up for rides. At any point (after finishing a ride) an individual may go home, or alternatively they may go to the next ride, and so-fourth. However, this is just a rough motivation and there are several significant differences between the generalized Jackson network models you will simulate and what may constitute something like an amusement park. You are asked to discuss these differences in Task 1. Below is a detailed mathematical description of the model.
There are $L$ stations (e.g. $L=5$ or $L=100$) and each station has a server and a waiting buffer (or queue). There are jobs traversing through the network and the stations. If needed, jobs can be uniquely identified (numbered) via $\ell=1,2,3,\ldots$ and it is expected that the simulation processes thousands, or millions of job, so that if you track individual jobs then $\ell$ is an increasing counter of the jobs.
A job's location can be one of:
Out the system prior to arrival (the job still doesn't exist).
In the system in a buffer of one of the stations (queueing).
In the system at one of the stations being served by the server at that station. This means it is "at the head" of the queue (and being served).
Out of the system after leaving.
Warning: With the exception of the final task, your simulation should not have entities for individual jobs, but rather only have counts of jobs at the stations as the system state.
Time is continuous, $t \ge 0$. In this continuous time system, jobs arrive to the system according to a random process and traverse between the stations as we describe now.
Each station $i=1,\ldots,L$ has a service rate $\mu_i > 0$ which is the inverse of the mean expected time it takes the server to process a job. Each station has a queue which has infinite capacity - i.e. there is no limit on the number of jobs in the queue. Jobs are served one by one at each station, where the mean service duration of each job is $\mu_i^{-1}$ and that duration is statistically independent from all other services at that node and at other nodes.
When a job arrives to station $i$, if the server is free (and this means there are no jobs in the buffer), the job immediately starts to receive service from the server. If the server is busy, the job queues for service and waits for its turn. Jobs may arrive either after finishing service at another node, or directly from the outside world.
The rate of arrivals of jobs to nodes from the outside world is $\alpha_i \ge 0$ for node $i$. Hence for nodes $i$ with external arrivals ($\alpha_i > 0$), the mean duration between external arrivals is $1/\alpha_i$. In this project we'll consider the durations of times between external arrivals to be exponentially distributed (this is gamma distributed with a squared coefficient of variation of $1$). This makes external arrival processes "Poisson processes" (this is also a mathematical detail that is not important for the carrying out the project).
As for the durations of service times we will set them as gamma distributed with a ratio of the mean and variance (squared coefficient of variation) which is $c_s >0$ for the services. Note that the case of $c_s=1$ is special as it means that the gamma distribution is actually an exponential distribution. In that case there are multiple analytical simplifications of such a system and it is no longer a "generalized" Jackson network but rather simply a "Jackson network" and can be also described by a "continuous time Markov chain" (you do not need to know these details for this project). In such a case, some explicit results exist and you could use these to test your simulations (and this is the part of Task 2 and Task 3). However when $c_s \neq 1$, explicit results are much more scarce and hence the simulation model may yield information that is not available otherwise (this is Task 4).
When a job completes service at a buffer it either leaves the system immediately, or immediately moves to another buffer. The choice of its destination is based on the routing matrix $P$. Here $P_{i,j}$ is probability of going to station j after completing service in station $i$. It is assumed that $P$ has non-negative entries and has a spectral radius less than $1$. Mathematically this means that the maximal absolute value of eigenvalues of $P$ is less than unity. It also means that jobs may always eventually leave the system and that at least some rows of the matrix $P$ sum up to less than $1$. The probability of a job leaving the system (immediately) after service at station $i$ is $1 - \sum_{j=1}^L P_{i,j}$.
Tip: Recall (4) from Question 14 of BigHW. The sample() function from StatsBase.jl may help generate the location (node) of the next job when carrying out the simulation. Yet keep in mind that "outside of the system" is also a possible location.
Hence in summary, the parameters for a network are:
The dimension of the network (number of nodes): $L$.
The external arrival rates: $\alpha_1,\ldots,\alpha_L$, with $\alpha_i \ge 0$.
The service rates: $\mu_1,\ldots,\mu_L$, with $\mu_i > 0$.
The $L \times L$ routing matrix $P$.
The squared coefficient of variation of the service processes: $c_s > 0$.
Tip: This bit of code with the rate_scv_gamma() function may help you see how to generate gamma distributed random variables with the desired rate and squared coefficient of variation.
using Distributions, Statistics """ A convenience function to make a Gamma distribution with desired rate (inverse of shape) and SCV. """ rate_scv_gamma(desired_rate::Float64, desired_scv::Float64) = Gamma(1/desired_scv, desired_scv/desired_rate) μ = 1.3 c_s = 2.1 dist = rate_scv_gamma(μ, c_s) println("Theoretical:") @show 1/mean(dist) #this is the mean() method from Distributions (theoretical mean) @show var(dist)/mean(dist)^2 #same for the var() method (theoretical variance) println("\nMonte Carlo:") data = [rand(dist) for _ in 1:10^6] @show 1/mean(data) #This is mean() from Base/Distributions (the sample mean) @show var(data)/mean(data)^2; #same for var() (the sample variance)
Theoretical: 1 / mean(dist) = 1.3 var(dist) / mean(dist) ^ 2 = 2.1000000000000005 Monte Carlo: 1 / mean(data) = 1.2997609770695828 var(data) / mean(data) ^ 2 = 2.1002151682461885
Arrival rates to nodes and traffic equations
Jobs arrive to nodes either from outside our internally after being routed to a node. Combining both of these arrivals, a basic quantity one can measure is the rate of arrivals to a given node $i$. That is,
\[ \lambda_i = \lim_{t \to \infty} \frac{A_i(t)}{t}, \]
where $A_i(t)$ is the number of arrivals (either from the outside world or from other nodes) to node $i$ during $[0,t]$. This quantity can be measured via simulation (as you will do), but it can also be computed directly via the matrix equation,
\[ \lambda = P^\top \lambda + \alpha, \qquad \qquad \text{(Traffic equations)} \]
where $\lambda$ is a (column) vector of $\lambda_i$ values, and $\alpha$ and $\mu$ are also column vectors of the respective quantities. These "Traffic equations" hold as long as the solution to the equation $\lambda = (I-P^\top)^{-1} \alpha$ satisfies $\lambda < \mu$ (where the inequality between vectors requires element wise inequality on all elements). To understand the traffic equations further, look at them for each $i$:
\[ \lambda_i = \alpha_i + \sum_{j=1}^L P_{j,i} \lambda_j. \]
Intuition: Observe that $\lambda_i$ is equated to the sum of external arrivals and all internal arrivals from other nodes $j$ routed with proportions $P_{j,i}$. This holds as long as the input rates to all nodes $j$ are less than the service rates $\mu_j$.
When $\lambda < \mu$ we say the network is stable and this means that the arrival rate to each node is less than the service rate. In this stable case, as time goes to infinity, the network converges to statistical equilibrium or "steady state". This means that if we consider a single single long simulation, we may run a time average of that simulation to obtain steady state quantities.
We will denote the load for node $i$ via $\rho_i = \lambda_i/\mu_i$ and always work in a regime where $\rho_i < 1$ for all $i=1,\ldots,L$. We denote $\rho^* = \max \{\rho_1,\ldots,\rho_L\}$ as the "bottleneck load". Clearly, as long as $\rho^*<1$ the network is stable.
Known queueing theory results and an overview of questions to investigate
In general, it is known that as each $\rho_i$ grows, the congestion in the individual queue at node $i$ increases in the sense that "on average", there are larger steady state queues and longer wait times. Further, as $\rho^*$ grows, the general congestion in the network grows. Since we consider $\rho^* < 1$, we only focus on networks that are stable, yet if $\rho^* \approx 1$, e.g. $=0.95$, they are highly congested.
It is also generally known that service time variability, quantified via $c_s$, positively affects queues sizes. For example if $c_s =0.01$ the service times at queues are almost deterministic and there is less variability at the queues and thus mean queue lengths are less then say a case of $c_s = 3.0$. As already stated, the case of $c_s = 1$ (with exponentially distributed service times) is quite special and the network is a "Jackson network" (not "Generalized" per se).
In this case of $c_s = 1$ (with exponential service times), one can prove that the steady state queue length of each queue has a "geometric distribution" on $0, 1, 2,\ldots$ with a mean of $\rho_i/(1-\rho_i)$. Many more such theoretical results exists, yet we skip the details. This relationship appears as follows:
using Plots ρ_grid =0:0.01:0.95 mean_steady_state_queue_size(ρ) = ρ/(1-ρ) #Theoretical steady state mean queue length (also of M/M/1 queue). plot(ρ_grid, mean_steady_state_queue_size.(ρ_grid), xlabel = "ρ", ylabel = "Mean steady state queue length\nfor an individual node", label = false, lw = 2, c = :black, xlim = (0,1),ylim=(0,20))
However for cases where $c_s \neq 1$ there generally not an explicit relationship. One of the main goals of the project is to quantify the effect of $c_s$ on the mean total queue sizes.
Tip: In the tasks below you are not expected to look at a node per-se but rather at the whole network. You would thus create a superposition of graphs as above, where you sum up all of the queues in the network.
Example scenarios (networks)
You will work with the following 4 scenarios. These are networks of size 3 nodes, 3 nodes, 5 nodes, or 100 nodes.
#You need to install the Parameters.jl package: https://github.com/mauro3/Parameters.jl #You need to install the Accessors.jl package: https://github.com/JuliaObjects/Accessors.jl using Parameters using Accessors using LinearAlgebra using Random @with_kw struct NetworkParameters #The @with_kw macro comes from the Parameters.jl package and makes nice constructors L::Int α_vector::Vector{Float64} #This vector is a vector of α_i which can then be scaled μ_vector::Vector{Float64} #This is the vector of service rates considered fixed P::Matrix{Float64} #routing matrix c_s::Float64 = 1.0 #The squared coefficient of variation of the service times with a default value of 1.0 end ############################ # Three queues in tandem scenario1 = NetworkParameters( L=3, α_vector = [0.5, 0, 0], μ_vector = ones(3), P = [0 1.0 0; 0 0 1.0; 0 0 0]) ############################ # Three queues in tandem with option to return back to first queue scenario2 = @set scenario1.P = [0 1.0 0; #The @set macro is from Accessors.jl and allows to easily make a 0 0 1.0; # modified copied of an (immutable) struct 0.3 0 0] ############################ # A ring of 5 queues scenario3 = NetworkParameters( L=5, α_vector = ones(5), μ_vector = collect(1:5), P = [0 .8 0 0 0; 0 0 .8 0 0; 0 0 0 .8 0; 0 0 0 0 .8; .8 0 0 0 0]) ############################ # A large arbitrary network #Generate some random(arbitrary) matrix P Random.seed!(0) L = 100 P = rand(L,L) P = P ./ sum(P, dims=2) #normalize rows by the sum P = P .* (0.2 .+ 0.7rand(L)) # multiply rows by factors in [0.2,0.9] scenario4 = NetworkParameters( L=L, α_vector = ones(L), μ_vector = 0.5 .+ rand(L), P = P);
Each of these 4 scenarios will be adjusted for desired $\rho^*$ and $c_s$. The way we adjust for a desired $\rho^*$ is by adjusting the α_vector. Details below.
Task 1: Thinking about Modelling Assumptions (5pts)
This task does not involve any coding or computation. Your answer can be in paragraph form or bullet point form. Consider the Generalized Jackson Network model and the application of a "Dream World", like amusement park. As mentioned the model can be used to model customer congestion in the amusement park. However, like any model it is at best a rough approximation of reality. Describe some features that exist in the amusement park setting that are not captured by the model. Similarly, describe model assumptions that may be unrealistic for the application of the amusement park.
Also answer the following: If you were actually using this model to assess congestion at the park, would it be useful or not? You may also use insight from Anna Foeglein's lecture.
Task 2: Computation for Jackson Networks (10pts)
We provide you with the functions maximal_alpha_scaling() and set_scenario() for adapting parameters of a scenario for a desired $\rho$ value and $c_s$ value.
""" Compute the maximal value by which we can scale the α_vector and be stable. """ function maximal_alpha_scaling(net::NetworkParameters) λ_base = (I - net.P') \ net.α_vector #Solve the traffic equations ρ_base = λ_base ./ net.μ_vector #Determine the load ρ return minimum(1 ./ ρ_base) #Return the maximal value by end max_scalings = round.(maximal_alpha_scaling.([scenario1, scenario2, scenario3, scenario4]),digits=3) println("The maximal scalings for scenarios 1 to 4 are: $max_scalings") """ Use this function to adjust the network parameters to the desired ρ⋆ and c_s """ function set_scenario(net::NetworkParameters, ρ::Float64, c_s::Float64 = 1.0) (ρ ≤ 0 || ρ ≥ 1) && error("ρ is out of range") max_scaling = maximal_alpha_scaling(net) net = @set net.α_vector = net.α_vector*max_scaling*ρ net = @set net.c_s = c_s return net end;
The maximal scalings for scenarios 1 to 4 are: [2.0, 1.4, 0.2, 0.216]
For example set the parameters of scenario 4 to have ρ⋆ = 0.7 and c_s = 2.4
#Adjust scenario 4 for a desired ρ and c_s, adjusted_net is the adjusted network adjusted_net = set_scenario(scenario4, 0.7, 2.4) #We can check by solving the traffic equations λ = (I - adjusted_net.P') \ adjusted_net.α_vector #Solve the traffic equations ρ = λ ./ adjusted_net.μ_vector #This is the vector of ρ values ρ_star= maximum(ρ) #\star + [TAB] @show ρ_star;
ρ_star = 0.7000000000000002
Now assume $c_s = 1$ and create 4 plots for the four scenarios where you plot the total steady state mean queue lengths as a function of $\rho^*$. Simply use the fact that the theoretical steady state mean queue length for queue $i$ is $\rho_i/(1-\rho_i)$. This is a straight forward task that does not involve simulation, but rather only a deterministic computation.
Task 3: A Simulation Engine Based on Queue States (40pts)
Create a simulation engine which is wrapped with the following function.
""" Runs a discrete event simulation of an Open Generalized Jackson Network `net`. The simulation runs from time `0` to `max_time`. Statistics about the total mean queue lengths are recorded from `warm_up_time` onwards and the estimated value is returned. This simulation does NOT keep individual customers state, it only keeps the state which is the number of items in each of the nodes. """ function sim_net(net::NetworkParameters; max_time = 10^6, warm_up_time = 10^4, seed::Int64 = 42)::Float64 #Set the random seed Random.seed!(seed) #The function would execute the simulation here - and this would use multiple other functions, types, and files estimated_total_mean_queue_length = 12.3 #the function would estimate this value (12.3 is just a place-holder) return estimated_total_mean_queue_length end;
Test the simulation in two ways:
Test #1: By running with $c_s = 1.0$ on a grid of $\rho^*$ with steps of $0.01$ between $0.1$ and $0.9$. Compare the values with the theoretical values that you plotted in the previous task. Compare both by plotting the values and computing the absolute relative error over the grid of $\rho^*$. You will know that your simulation is valid if the absolute relative error will be low and further decrease as max_time is increased.
Test #2: Modify the sim_net() function so that $A_i(T)$ (total overall arrivals to node $i$ over the whole time horizon $T$) is also recorded. Then by computing $A_i(T)/T$ you should get values that agree with (or are near) $\lambda_i$. Note that changing $c_s$ should not affect this - so try for various settings of $c_s$.
These tests may be placed in your test folder. For them to be automatic you would require thresholds as they are statistical (noisy tests). So making the tests automatic is not required, but only optional.
Task 4: Analysis of the Effect of $c_s$ on Queue Sizes (15pts)
Now run further simulations for the same grid of $\rho^*$ but for $c_s \in \{0.1, 0.5, 1.0, 2.0, 4.0\}$. Use these simulations to add plots of the mean total queue lengths for the various $c_s$ values. So you should now 5 plots on each figure, and we expect to see that as $c_s$ increases the curves increase as well. Make sure the legend is properly marked and placed on the plot.
In your analysis, run multiple runs per setting to obtain confidence bounds on the curves. Plot the curves in a neat way illustrating the confidence bounds. Here you have some freedom in how to do this.
Task 5: Sojourn Time Distributions (10pts)
Note: This final task requires a significant amount of coding (not proportional to the point weighting) and testing on top of the other tasks.
Now create another parallel simulation module wrapped via a function sim_net_customers(). In contrast to the simulation engine you created before, this simulation engine keeps track of each of the customers in the system. One way to do this is to use a Queue data structure from DataStructures.jl in a way somewhat similar manner this basic example. Then each node will have a queue of customers implemented via Queue and each customer may have a time-stamp of the time of arrival to the system and time of departure.
Your goal is now to consider each of the scenarios and for each one to estimate Q1, the median, and Q3 of the sojourn time in the system (this is the time a customer spends in the system) when $\rho^* = 0.8$ and $c_s \in \{0.5, 1.0, 2.0\}$.
Place the results in a neatly formatted table. Optionally run multiple runs per setting to obtain error bounds for your estimates.