example_sim_minimal_GGc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimal G/G/c simulator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to simulate a G/G/c model (without using any QueueSim classes)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generating multiple pseudo random numbers as an array and fast array processing\n",
    "import numpy as np\n",
    "\n",
    "# Comparison with analytical results\n",
    "from queuesim.analytic import erlang_c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Means inter-arrival time\n",
    "mean_i = 100\n",
    "\n",
    "# Means service time\n",
    "mean_s = 240\n",
    "\n",
    "# Number of operators\n",
    "c = 3\n",
    "\n",
    "# Number of arrivals to be simulated\n",
    "count = 1_000_000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate random arrival times and service times"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of `np.random.exponential` any other random distribution could also be used. The simulation process is not limited to exponential inter-arrival and service times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "arrival_times = np.cumsum(np.random.exponential(mean_i, count))  # Absolute arrival times, not inter-arrival times\n",
    "service_times = np.random.exponential(mean_s, count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The complete simulation code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "waiting_sum = 0\n",
    "server_free_at = np.zeros(c)\n",
    "\n",
    "for (time, service_time) in zip(arrival_times, service_times):\n",
    "    index = np.argmin(server_free_at)  # Find next free server\n",
    "    service_start = max(time, server_free_at[index])\n",
    "    waiting_sum += service_start - time\n",
    "    server_free_at[index] = service_start + service_time  # Store new time when server gets idle again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Output simulation results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[NQ]=2.57\n",
      "E[N]=4.97\n",
      "E[W]=257.41\n",
      "E[V]=497.41\n",
      "ρ=79.9%\n"
     ]
    }
   ],
   "source": [
    "service_sum = np.sum(service_times)\n",
    "l = 1 / np.mean(arrival_times[1:] - arrival_times[:-1])  # Calculate lambda from absolut arrival times\n",
    "\n",
    "EW = waiting_sum / count\n",
    "EV = (waiting_sum + service_sum) / count\n",
    "ENQ = EW * l  # Since Little holds for all types of inter-arrival and service time distributions\n",
    "EN = EV * l  # we can caluclate E[NQ] and E[N] from E[W] and E[V] and do not need to record them directly\n",
    "rho = service_sum / arrival_times[-1] / c  # Sum of all service times divided by last arrival time divided by c\n",
    "\n",
    "print(\"E[NQ]=\", round(ENQ, 2), sep=\"\")\n",
    "print(\"E[N]=\", round(EN, 2), sep=\"\")\n",
    "print(\"E[W]=\", round(EW, 2), sep=\"\")\n",
    "print(\"E[V]=\", round(EV, 2), sep=\"\")\n",
    "print(\"\\N{greek small letter rho}=\", round(rho * 100, 1), \"%\", sep=\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison with analytics results (Erlang C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E[NQ]=2.59\n",
      "E[N]=4.99\n",
      "E[W]=258.88\n",
      "E[V]=498.88\n",
      "ρ=80.0%\n"
     ]
    }
   ],
   "source": [
    "result = erlang_c(1 / mean_i, 1 / mean_s, c)\n",
    "\n",
    "print(\"E[NQ]=\", round(result.ENQ, 2), sep=\"\")\n",
    "print(\"E[N]=\", round(result.EN, 2), sep=\"\")\n",
    "print(\"E[W]=\", round(result.EW, 2), sep=\"\")\n",
    "print(\"E[V]=\", round(result.EV, 2), sep=\"\")\n",
    "print(\"\\N{greek small letter rho}=\", round(result.rho * 100, 1), \"%\", sep=\"\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "37056784771c7b5ac6ab9d7dda111ef842cb90abd361c4ad13f8b6bfd3e8f1f2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}