%\documentclass[twocolumn,12pt]{article}
\documentclass[12pt]{article}
 
\usepackage[utf8]{inputenc}    % 处理 UTF-8
\usepackage{amsmath,amssymb}   % 数学公式支持
\usepackage{graphicx}          % 插入图片
\usepackage[colorlinks=true, allcolors=blue]{hyperref} % 链接
\usepackage{geometry}          % 页边距
\usepackage{setspace}          % 行距
\usepackage{titlesec}          % 标题格式更灵活
\usepackage{fontspec}
\setmainfont{Times New Roman}
\usepackage{xeCJK}
\doublespacing
\usepackage[margin=1in]{geometry}
\raggedright
\usepackage{indentfirst}
\setlength{\parindent}{1.27cm}
\setlength{\parskip}{0pt}  
\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhf{}
\fancyhead[R]{\thepage}  % 右上角姓+页码
\usepackage{float}
\usepackage{caption}  % (导言区加,只要加一次)
\usepackage{booktabs} % 在导言区加这个
\usepackage{stfloats}
 
 
 
% 一级标题 APA 样式:居中、加粗、段前12pt、段后6pt
\titleformat{\section}
  {\normalfont\bfseries\centering\fontsize{12}{14}\selectfont}
  {\thesection}{1em}{}
\titlespacing*{\section}
  {0pt}{12pt plus 2pt minus 2pt}{6pt}
 
% 二级标题 APA 样式:左对齐、加粗、段前12pt、段后6pt
\titleformat{\subsection}
  {\normalfont\bfseries\fontsize{12}{14}\selectfont}
  {\thesubsection}{1em}{}
\titlespacing*{\subsection}
  {0pt}{12pt plus 2pt minus 2pt}{6pt}
 
  \titleformat{\subsubsection}
  {\bfseries\itshape\flushleft\normalsize} % 加粗 + 斜体 + 左对齐
  {}{0pt}{}
 
\usepackage{titlesec}
\usepackage{fontspec}
\setmainfont{Times New Roman}
\renewcommand{\headrulewidth}{0pt}
 
\titlespacing*{\section}
  {0pt}{0pt}{0pt}
 
% 页面设置
\geometry{
  a4paper,
  left=2.45cm,
  right=2.45cm,
  top=2.45cm,
  bottom=2.45cm
}
 
 
\begin{document}
\begin{titlepage}
    \centering
    \vspace*{5cm}  % 可调整让内容垂直居中
    
    {\normalsize\textbf{Queueing-Behavior Integrated Modeling for Waiting Time Reduction and Crowding Mitigation: a case study of school cafeteria}\par}
    
    \vspace{2cm}
    {\normalsize Word Count: 4954\par}
\end{titlepage}
 
 
% 您的内容从此开始
% ==============================
% 原始 Markdown 第一行: # ***[Over crowding + ]***
 
 
% 原始 Markdown 里的 “### Abstract” -> \section
 
\section*{Introduction}
 
In my own experience of our school, I often had to spend over fifteen minutes in line, only to wait for food. This personal frustration led me to question whether these delays were simply due to high demand or rooted in deeper systemic inefficiencies. Long queue times and overcrowding in school cafeterias are known to be serious problems facing students. This issue has been well documented in previous studies that have emphasized the impact of poorly designed cafeterias on student satisfaction, waiting times, and the overall dining experience (Galabo, 2019). Surveys conducted at our school have shown that many students are dissatisfied with the current cafeteria setup, citing overcrowding and long waiting times as major annoyances. First, we break down the students' dissatisfaction into two main areas. Firstly, the long waiting time consumes their patience, which is totally unnecessary. Second, the crowded cafeteria makes it difficult for them to move around. 
 
However, no one has fully explored why these issues persist after years of improvement or how specific factors, such as cafeteria layout, queueing systems, and crowding during peak hours, contribute to the problem. If we can identify the root causes of inefficient cafeteria operations, we can not only improve the dining experience but also optimize space usage and human resources. This study aims to investigate the root causes of cafeteria crowding and long waiting times by analyzing queue dynamics, conducting simulations, and exploring potential improvements. By identifying and addressing these issues, I hope to reduce waiting times and crowding situations to increase student satisfaction and to improve overall operational efficiency. And I also hope to use this analysis to provide a generic template and quantitative solution for other similar situations.
 
\subsection*{Literature Review}
 
To understand the mechanisms behind these inefficiencies, it is essential to examine existing theoretical frameworks and studies related to queueing behavior. One of the most widely used approaches in this area is queueing theory, which provides mathematical tools to analyze waiting lines. Queueing theory, a subfield of operations research, focuses on analyzing and optimizing systems where entities wait in line for service. The foundational models, such as $M/M/1$, $M/M/c$, and $M/D/1$, have been widely adopted to understand arrival patterns, service rates, and queue dynamics under probabilistic assumptions (Kambli et al., 2020). These models are defined by Kendall's notation, where "M" refers to Markovian (Poisson) arrival and/or service times, and the number indicates the quantity of parallel servers (Ajiboye \& Saminu, 2018). Little's Law, a fundamental result in queueing theory, relates average waiting time to arrival rate and system capacity (Tang et al., 2019). 
 
In the context of campus cafeterias, queueing theory has been extensively employed to address issues of student congestion during lunch hours, which are particularly pronounced due to synchronized class schedules (Zuo, 2020). The application of models such as M/M/n/m has shown practical potential in optimizing the number of open windows and reducing average waiting time by reallocating resources during peak demand (Chen \& Wang, 2011). Simulation studies demonstrate that queue optimization can reduce waiting times by nearly 45\% (Ye, 2009). Another study using a multi-stage queue model reduced average customer dwell time from over 500 minutes to under 15 minutes by redesigning queue architecture and service logic (Ajiboye \& Saminu, 2018). 
 
This methodology was selected because it enables the formulation of a structured, quantitative framework to address crowding issues through measurable parameters such as arrival rate, service rate, and waiting time. Unlike anecdotal methods, queueing theory allows for reproducible simulation and modeling, which can be validated and iteratively refined. Additionally, this approach provides rapid validation of the solution. In most cases, it is sufficient to modify the input variables and compare the output with the previous one.
 
\subsection*{Research Gap}
 
Despite these successes, classical queueing models exhibit limitations when applied to real-world cafeteria environments. Most models assume steady-state behavior, homogeneous service stations, and exponentially distributed arrival and service times, which rarely match the fluctuating and clustered nature of student behavior (Lu et al., 2021). Furthermore, traditional models fail to account for human behavioral dynamics, such as balking, reneging, or jockeying, where students leave or switch queues due to perceived delays (Chen and Wang, 2011). Extensions to these models have attempted to incorporate dynamic arrival rates and customer patience thresholds (Li \& Saminu, 2016). Moreover, complex queue networks involving multiple food stations and checkout counters require simulation or agent-based modeling, as closed-form analytical solutions are infeasible (Kambli et al., 2020; Deng et al., 2011).
 
And there is no suitable theoretical model to describe dynamic variables, such as, in this study, arrival rates that change over time (students arrive in a concentrated manner during peak hours, followed by a gradual decrease in arrival rates). In this study, an attempt is made to input dynamic variables over time into the model by replacing static variables with dynamic variables and eventually obtaining outputs that vary over time. In other words, all the variables of the original model are upgraded to time-varying functions, so the final results obtained are also in functional form with the independent variable time. Also, in this study, there is a simplified analysis of student behavior. This can be highly integrated with the queueing theory model and applies to a very wide range of scenarios.
 
% GAP 对于特定时间,而不是持续时间
 
\section*{Methodology}
 
\subsection*{Spatial Layout of the Cafeteria}
 
Figure~\ref{fig:layout} illustrates the layout of the cafeteria used in this study. The topmost area is the dining area and the remaining area is the working area. There are two entrances to the dining area, which are located on the upper left and lower left. It is important to note that there are two dining areas in total. The indoor dining area is shown in the figure. There is also an outdoor dining area (not shown), which can be reached by leaving through the lower left entrance.
 
\begin{figure}[H]
    \centering
    \includegraphics[width=0.25\linewidth]{WechatIMG39.jpg}
    \caption{Layout of school cafeteria}
    \label{fig:layout}
\end{figure}
 
\subsection*{Data Collection}
 
The data for this study came from field observations in our school's cafeteria, and data were collected from December 2 to December 5, 2024, and December 9 to December 12, 2024 (8 days in total). Data collection under each schedule lasted for four days due to the school's schedule of one cycle every two days. The data primarily covered the third and fourth periods which are lunch periods, and the peak hours. Specifically, it includes the number of students entering/leaving the cafeteria during each period, the time taken for services, and the number of windows.
 
Data collection was performed through the cafeteria's monitoring system during the daytime hours, which recorded the timestamp of each student's entry into the cafeteria. The specific time period for data collected was from 11:50 to 13:35 each day. During this time period, the number of students entering the cafeteria was recorded by manual count, and the service time of each service was recorded by simple observations.
 
\subsection*{Ethical Considerations}
 
Data collection was conducted solely for research purposes with the approval of the school administration. No personally identifiable information of students was collected; only the data such as entry and exit counts and service times were recorded. Observations were carried out in a non-intrusive manner to avoid interfering with students' normal activities. So the study complies with institutional ethical guidelines for research involving human subjects.
 
\subsection*{Preprocessing and Entry-Exit Analysis}
 
During the data preprocessing process, all the collected data met the preset criteria, so no data was eliminated. During the data preprocessing stage, timestamps were converted to relative times with respect to peak hours and grouped for counting at 2-minute intervals during peak hours. For off-peak hours, group counts were performed at 5-minute intervals.
 
We get the average Entries and Exits ($E_t$ and $L_t$) by simply adding them together and divided by the number of days. Notice that our school's cafeteria has 2 entrances (marked as A and B here).
So for each time interval t,
\begin{align*}
\text{Enter}_{\text{Total}} &=
\frac{\sum_{i=1}^{8} \text{Enter A}_{\text{Day}i}}{8} 
 + \frac{\sum_{i=1}^{8} \text{Enter B}_{\text{Day}i}}{8}.
\end{align*}
It is the same way to calculate ${\text{Leave}_{\text{Total}}}$.
 
\begin{figure}[htbp]
    \centering
    \begin{minipage}[t]{0.48\linewidth}
        \centering
        \includegraphics[width=\linewidth]{EnterLeave6.jpg}
        \caption{Average Enter vs Leave Counts}
        \label{fig:EL}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{0.48\linewidth}
        \centering
        \includegraphics[width=\linewidth]{net_population_final_uniformed.png}
        \caption{Average Net Population}
        \label{fig:NP}
    \end{minipage}
\end{figure}
 
 
According to Figure~\ref{fig:EL}, the results of the data analysis showed that the arrival time of students to the cafeteria showed obvious peak characteristics, especially between 11:50 and 12:10 (0 min - 20 min), when an average of about 20 people entered the cafeteria per minute, while the number of arrivals decreased significantly at other times.
 
We can also graph the net population in the cafeteria. According to Figure~\ref{fig:NP}, similarly, there are particularly large numbers of students at the starts of the 2 launch periods, and then the numbers drop off over time.
 
\subsection*{Arrival Rate Distribution}
Once the data collection was complete we next used the data to calculate the student arrival rate, which is a variable that changes over time during every lunch time.
Then we consider the following factors. We know that not all of the students who enter the cafeteria will have any food, and that some may only stay for a short time before leaving (e.g., using the cafeteria as a hallway just to walk through it). In addition to this, our school has an outdoor dining area. It has about the same capacity as an indoor cafeteria. Some people order their lunch and go straight to the outdoor dining area, while others stay indoors to eat their lunch. We will mainly focus on those students who stay indoors.
 
\begin{table*}[htbp]
\centering
\caption{Definitions of variables used in the behavior analysis}
\begin{tabular}{l p{9cm} l}
\toprule
\textbf{Symbol} & \textbf{Description} & \textbf{Unit} \\
\midrule
$E_t$        & Total number of individuals entering the location during interval $t$  & Count    \\
$L_t$        & Total number of individuals leaving the location during interval $t$   & Count    \\
$p(t)$       & Probability that an entering individual seeks service                  & Fraction \\
$r$          & Probability that a service seeker stays after service                  & Fraction \\
$S_t$        & Number of individuals arriving for service during interval $t$         & Count    \\
$\lambda(t)$ & Average arrival rate of customers during interval $t$                  & People/s \\
\bottomrule
\end{tabular}
 
\label{tab:Model1Var}
\end{table*}
 
 
According to Table~\ref{tab:Model1Var}, we have a fraction $p$ of the entrants seek service, while the others leave immediately. After services, some service seekers leave immediately as well, and some stay after the service. The fraction $r$ of service seekers stay in the area. However, $r$ does not vary roughly with the time a day because people will always look for uncrowded locations. In this case, we suggest that they prefer to stay indoors when it is crowded outdoors, and they prefer to go outside when it is crowded indoors. But it has been observed that it increases significantly during the cold season, when people are reluctant to have lunch outdoors where it is colder. So in winter, it usually gets more crowded indoors.
 
 
\subsubsection*{Individual Behavioral Performance Breakdown}
\begin{figure}[H]
    \centering
    \includegraphics[width=0.3\linewidth]{StudentWay.png}
    \caption{Student Behavioral Flowchart}
    \label{fig:student_flow}
\end{figure}
In this case, we assume the simplest "ideal man" behavior. According to Figure~\ref{fig:student_flow}, we assume that there are only two cases:
\begin{enumerate} 
\item After entering, individuals who do not purchase food (proportion $1-p(t)$) leave immediately, with their presence time considered negligible; 
\item After entering, individuals who proceed to purchase food(proportion $p(t)$) remain in the cafeteria for a period before departing.
This period consists of service time $T$ and dining time in the cafeteria $D$.
Assuming a proportion $r$ of individuals choose to eat indoors, the total residence time can be approximated as 
\[
\tau_{net} \approx T + rD. \tag{1}
\]
\end{enumerate}
Since we only have external data (entry rate $E_t$ and exit rate $L_t$ in a given time period t), we cannot get the number of people who want to be served at a given time period because some of them are just passing through. So we set up several scenarios: the portion of the population that buys meals and eats them in the cafeteria leaves after about $\tau_1 = T + D$ seconds; the portion of the population that buys a meal and leaves immediately after about $\tau_2 = T$ seconds; and the portion of the population that does not buy a meal enters and exits almost instantaneously. From this, a conservation relation can be formulated:
\[
L_t \;\approx\; 
\underbrace{\bigl(1 - p(t)\bigr)\times E_t}_{\text{Immediately leave after entering}} 
\;+\; 
\underbrace{r \times p(t-\tau_1)\times E_{t-\tau_1}}_{\text{Have finished the meal and leave}} 
\;+\; 
\underbrace{(1 - r) \times S_{t - \tau_2}}_{\text{Have purchased food and now leave}}. \tag{2}
\]
In the real situation, however, $p(t)$ does not change rapidly. In previous observations, the proportion of people who went through the cafeteria only decreased slightly over time. So, it is possible to assume that $p(t - \tau_1)\approx p(t)$ since $\tau_1$ is not that long. This gives:
\[
L(t) \;\approx\; \bigl(1 - p(t)\bigr)\,E_t \;+\; p(t)\,E_{t - \tau_1} \;+\; r \, S_{t - \tau_2}.
\]
Simplify the equation and replace $S_{t - \tau_2}$ with $p(t)E_{t - \tau_2}$,
\[
\begin{aligned} 
L_t &= E_t \;-\; p(t)E_t \;+\; p(t)E_{t - \tau_1} \;+\; r \, p(t)E_{t - \tau_2}\\
&= E_t  \;+\; p(t)\,\Bigl[E_{t - \tau_1}\;-\;E_t \;+\; r \, E_{t - \tau_2} \Bigr]. \\
\end{aligned}
\]
Thus, an approximate solution can be obtained:
\[
p(t) \;\approx\; \dfrac{\,L_t \;-\; E_t}{\, E_{t - \tau_1} \;-\; E_t \;+\; rE_{t - \tau_2} \,}.
\tag{3}
\]
Again, based on the approximation in Equation (1), we combine $\tau_1$ and $\tau_2$ into a single term $\tau_{net}$,
\[
\boxed{ p(t) \;\approx\; \dfrac{\,L_t \;-\; E_t}{\, E_{t - \tau_{net}} \;-\; E_t \,} }.
\tag{4}
\]
According to $p(t) = \frac{S_t}{E_t}$ , we get
\[
\frac{S_t}{E_t} \approx \dfrac{\,L_t \;-\; E_t}{\, E_{t - \tau_{net}} \;-\; E_t \,}.
\]
The average time high school students spend actually eating lunch is between 7 and 10 minutes (Lee \& Shanklin, 2002). So we can calculate
\[
\tau_{net} \approx 8.5 \times 60 \times r + T = 510 \, r + T.
\]
For $r$ we use the indoor and outdoor seating capacity to estimates directly, since people usually go to places where seats are more available.
\[
r = \frac{N_{indoor}}{N_{outdoor} + N_{indoor}} \approx \frac{47}{59 + 47} \approx 0.443.
\]
Finally we divide $S$ by $t$ to get the student arrival rate:
\[
\bar{\lambda_t} = \frac{S_t}{t_{\text{period}}}.
\]
 
\vspace{-\baselineskip}
 
\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\linewidth]{lambda.png}
    \caption{Arrival rate (λ) \& Standard Deviation (σ)}
    \label{fig:lambda}
\end{figure}
\vspace{-1\baselineskip}
 
The resulting data is shown in Figure~\ref{fig:lambda}. We first calculate the average customer arrival rate $\bar{\lambda}$ using the data $E$ and $L$ averaged over multiple days. After that we calculate $\lambda$ for each day's data separately and differ it from $\bar{\lambda}$ to get the standard deviation $\sigma$ from this formula:
\[
\sigma_t = \sqrt{\frac{1}{n-1} \sum_{i=1}^n \left(\lambda_{t,i} - \bar{\lambda_t}\right)^2}
\]
For this type of data we approximate it as a normal distribution, i.e., in this form
\[
X_t \sim \mathcal{N}(\lambda(t), \sigma(t)), \quad \forall t \in [0,6300].
\]
 
 
\subsubsection*{Smoothing Spline Fitting}
By converting discrete estimates into continuous representations, this method enables more accurate simulation and forecasting optimization within the Queueing theory framework (Pollock, 1993). To construct a continuous function for the average arrival rate $\lambda(t)$ as well as its standard deviation $\sigma(t)$, this study employs the Smoothing Spline Fitting (SSF) method. The method is based on a set of discrete time points $\{t_i\}_{i=1}^n$ with corresponding values $\{y_i\}_{i=1}^n$, and fits a smooth and continuous function $f(t)$ by controlling the degree of bending of the function while minimizing the residuals of the data (Pollock, 1993). This is expressed mathematically as the following optimization problem:
\[
\min_{f \in C^2} \left\{ \sum_{i=1}^n \left( y_i - f(t_i) \right)^2 + \eta \int_{a}^{b} \left( f''(t) \right)^2 \, dt \right\}.
\]
Here, $C^2$ denotes the function space of twice continuously differentiable functions; $y_i$ represents the target value computed for the $i$-th time interval; $f(t)$ is the continuous function to be fitted; and $\eta$ is the smoothing factor, which controls the trade-off between fitting accuracy and smoothness (Wahba, 1990). When $\eta = 0$, the fitted function passes exactly through all data points (i.e., interpolation); as $\lambda \to \infty$, the fitted function approaches linearity (Reinsch, 1967).
 
In this study, we employed the \texttt{UnivariateSpline} function from the \texttt{scipy.interpolate} module in Python to fit the functions $\lambda(t)$ and $\sigma(t)$ (Virtanen et al., 2020). This function constructs a cubic spline based on the input data, applies the default smoothing factor, and satisfies the following optimization formulation:
\[
\mu(t) \approx \sum_{i=0}^{n} c_i B_i(t).
\]
The resulting function $f(t)$ can be evaluated at any timestamp $t$, thereby generating continuous representations of the two functions $\lambda(t)$ and $\sigma(t)$. These continuous functions are used to characterize the dynamic customer arrival distribution within the queueing theory model.
 
\subsection*{Service Rate Distribution}
 
For the service rate data, this study recorded the durations of 100 service sessions during the lunch period from the surveillance. The data are shown in Figure~\ref{fig:servicetime} and the statistics are shown in Table~\ref{tab:servicetime}
 
\begin{figure}
  \centering
  % 占位图示例,请替换为真实图片文件
  \includegraphics[width=0.8\linewidth]{snapshot-1744736997962@2x.png}
  \caption{Service Time Data ($1/\mu$)}
  \label{fig:servicetime}
\end{figure}
 
Service rate refers to the number of services completed per unit of time (number per second). What is recorded here is the time taken for each service. So we need to take the reciprocal of the service duration recorded here. The average service rate is $\bar\mu \approx 0.030885$ (i.e., on average, people spend 32.378 seconds during the services).
The K-S (Kolmogorov-Smirnov) test was used in this study to test whether this data can be approximated to a certain distribution (Lilliefors, 1969). Based on the graph, we determine that it fits the right-skewed characteristic of the exponential distribution very well.
 
\begingroup
\renewcommand{\arraystretch}{1.0}
\singlespacing
 
\begin{table}[H]
    \centering
    \caption{Statistics of the Service Time Data}
    \begin{tabular}{lcc}
    \toprule
    \textbf{Variable}  & \textbf{Value}   \\
    \midrule
    Sample size        & 100       \\
    mean               & 32.378    \\
    std                & 20.478926 \\
    $\lambda$         & 0.030885  \\
    D                  & 0.3510    \\
    Pval               & 0.1045    \\
    \bottomrule
    \end{tabular}
    \label{tab:servicetime}
\end{table}
 
\endgroup
 
Therefore, we first used the \texttt{scipy.stats.kstest} function in Python to test whether the sample conforms to an exponential distribution (Virtanen et al., 2020). According to Table~\ref{tab:servicetime}, the K-S test yielded a p-value of 0.1045, which is greater than the significance level $\alpha = 0.05$. This indicates that we fail to reject the null hypothesis that the service time follows an exponential distribution, thereby supporting the rationality of approximating the data as an exponential distribution.
 
\subsection*{Model Establishment}
 
Based on the previous examination of the data, we believe that the arrival rate conforms to a normal distribution and the service rate conforms to an exponential distribution (i.e., the service process is memoryless). At the same time, the arrival rate varies over time, so the $G(t)/M/c$ queueing theory model is chosen to describe it, where $G(t)$ denotes a general time-dependent arrival process, $M$ indicates a service time distribution that satisfies a Markov process (exponential distribution), and $c$ represents the number of parallel servers (service windows) (Adan \& Resing, 2002).
 
\subsubsection*{Basic parameters}
 
\begingroup
\renewcommand{\arraystretch}{1.0}
\singlespacing
 
\begin{table}[H]
  \centering
  \caption{Variables in the model}
  \begin{tabular}{llc}
    \toprule
    \textbf{Symbol} & \textbf{Explanation} & \textbf{Unit} \\
    \midrule
    $\lambda(t)$    & Average arrival rate of customers   & people/s \\
    $\mu$           & Average service rate of a single service window & services/s \\
    $c$             & Number of service windows           & Number \\
    $\rho(t)$       & System utilization rate             & / \\
    $C_a(t)$        & Coefficient of variation            & / \\
    $W_q(t)$        & Mean waiting time                   & s \\
    $L_q(t)$        & Average number of people in the queue  & Number \\
    $P_0(t)$        & Probability that the system is empty & / \\
    \bottomrule
  \end{tabular}
  \label{tab:modelvar}
\end{table}
 
\endgroup
 
All the variables are shown in Table~\ref{tab:modelvar}, which the System utilization rate (Kleinrock, 1975):
\[
\rho(t) = \frac{\lambda(t)}{c\mu}, \quad \rho < 1 \;\text{(Condition of system stability)}.
\tag{5}
\]
The coefficient of Variation of the arrival time distribution $C_a(t)$ can be calculated bu the equation:
\[
C_a(t) = \frac{\sigma_a(t)}{1/\lambda(t)},
\tag{6}
\]
where $\sigma_a(t)$ denotes the standard deviation of the arrival rate computed above. We use $C_a(t)$, the coefficient of variation of the arrival time distribution, to adjust the $M/M/c$ formulas and thus obtain results for the $G/M/c$ model (Gross \& Harris, 1998).
 
\subsubsection*{Model Computations}
Using the generalized form of the Pollaczek–Khinchin (P–K) formula, the mean waiting time in queue is given by (Whitt, 1993):
\[
W_q = \frac{C_a^2(t) + 1}{2}\cdot W_{q,\,M/M/c},
\tag{7}
\]
where $W_{q,\,M/M/c}$ is the mean waiting time under the $M/M/c$ model (Little, 1961):
\[
W_{q,\,M/M/c} = \frac{L_{q,\,M/M/c}}{\lambda}.
\tag{8}
\]
The mean queue length $L_q$ can also be calculated:
\[
L_q = \frac{C_a^2(t) + 1}{2}\cdot L_{q,\,M/M/c}.
\tag{9}
\]
Here, $L_{q,\,M/M/c}$ (the mean number in queue) is given by the Erlang–C formula (Cooper, 1981):
\[
L_{q,\,M/M/c}
= \frac{P_0\,\frac{(\lambda/\mu)^c}{c!}\,\frac{c\rho}{(1-\rho)^2}}
       {\displaystyle\sum_{k=0}^{c-1}\frac{(\lambda/\mu)^k}{k!}
        + \frac{(\lambda/\mu)^c}{c!}\,\frac{1}{1-\rho}},
\tag{10}
\]
where the idle-system probability $P_0$ is
\[
P_0 = \left(\sum_{k=0}^{c-1}\frac{(\lambda/\mu)^k}{k!}
       + \frac{(\lambda/\mu)^c}{c!}\,\frac{1}{1-\rho}\right)^{-1}.
\tag{11}
\]
Finally, the mean system response time $W$ and the mean number in system $L$ are
\[
W = W_q + \frac{1}{\mu}, \quad
L = L_q + \frac{\lambda}{\mu}.
\tag{12}
\]
 
\begin{figure}[htbp]
    \centering
    \begin{minipage}[t]{0.48\linewidth}
        \centering
        \includegraphics[width=\linewidth]{RHO2.png}
        \caption{System Stability During Peak Hour}
        \label{fig:rho2}
    \end{minipage}
    \hfill
    \begin{minipage}[t]{0.48\linewidth}
        \centering
        \includegraphics[width=\linewidth]{AVG Result.png}
        \caption{Model Result (W and L)}
        \label{fig:avg res}
    \end{minipage}
\end{figure}
 
 
\section*{Results and Findings}
 
These output variables are then plotted for visualization. According to Figure~\ref{fig:rho2}, it can be observed that during 0 - 290 seconds, the system stability level $\rho$ rises significantly, exceeding 1, the service capacity is not stable. During this peak period (about the first five minutes), the student arrival rate exceeds the service rate, resulting in a situation where the service can never be completed. Under such circumstances, the system becomes unstable, and both the number of people waiting and the waiting time grow infinitely, with no steady-state behavior (see Figure~\ref{fig:avg res}).
 
From Equation~(5), we know that the condition for system stability is:
\[
\frac{\lambda}{c\mu} < 1, \quad \text{or equivalently,} \quad \lambda < c\mu.
\]
Therefore, based on the collected data in this study and assuming the number of service windows and the service rate remain unchanged, the arrival rate must satisfy the following condition in order for the system to remain stable:
\[
\lambda < 3 \times 0.030885 \quad \Rightarrow \quad \lambda < 0.09.
\]
 
\subsection*{Potential optimization}
 
For the optimization, due to spatial constraints within the cafeteria layout, a maximum of only four service windows can be opened. Moreover, communications between servers and students during the service process can significantly prolong service time, thereby reducing the service rate. Assuming an ideal scenario in which the average service time is reduced to 25 seconds (Reduced communications), the arrival rate must still satisfy:
\[
\lambda < 4 \times 0.04 \quad \Rightarrow \quad \lambda < 0.16.
\]
 
Based on these considerations, a preliminary optimization strategy is proposed. We need to limit the arrival rate anyway because there are still some times that the arrival rate is greater than 0.16 (see Figure~\ref{fig:lambda}), and the ideal limit is to control the arrival rate below 0.12 (make sure that the system stability rate $\rho$ is not too close to 1). Given that the proportion of students choosing to remain in the cafeteria ($r$) tends to increase during winter (Students don't prefer to eat outside because it gets cold), more students would be accommodated within the cafeteria. This data was also collected during the winter months. Generally, the number of people in the cafeteria decreases significantly during the summer months. As such, seasonally adapted control should be implemented: during winter, the arrival rate should be limited to 0.12, and during summer, the limit can be reduced, but it should still not exceed 0.14. As illustrated in Figure~\ref{fig:optimized}, if the arrival rate is restricted to below 0.12, the maximum waiting time stabilizes around 50 seconds, which falls within an acceptable and efficient operational range.
 
\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{Optimized WL.png}
    \caption{Optimized Result (W and L)}
    \label{fig:optimized}
\end{figure}
 
For a solution to control arrival rates, we need to develop a method that is the most effective. This would include having a teacher or student volunteer stand somewhere in the cafeteria and stop students from trying to enter periodically. The method then needs to determine an optimal standing position. Different positions could lead to more serious side effects, especially the possibility of obstructing the flow of individuals passing through the cafeteria.
 
\subsection*{Dynamic Simulation}
 
In order to determine the most effective location for limiting the arrival rate, this study used dynamic simulation to model the dining process to view the crowd density in different areas. In order to determine an appropriate location.
 
\begin{figure}[H]
  \centering
  \includegraphics[width=0.45\textwidth]{Pasted image 20250415232238.png}
  \caption{Dynamic Simulation Preparation}
  \label{fig:dynamic}
\end{figure}
 
We configured the elements based on the plan figure (Figure~\ref{fig:layout}) of the cafeteria. As shown in Figure~\ref{fig:dynamic}, the orange areas represent walls, which are considered impassable for pedestrians. Due to the narrow spacing between chairs and tables, it is assumed that pedestrians cannot pass through the central seating area. The green line segments at the upper-left and lower-left corners indicate the entrance and exit points of the cafeteria. The green arrows denote the queueing direction, while the green dots represent the positions of the service windows. Additionally, the green shaded area on the left side indicates a corridor for through-passage.
 
Subsequently, we employed agent-based modeling techniques to simulate pedestrian behaviors within the cafeteria, utilizing the Pedestrian Library provided by AnyLogic (AnyLogic Company, n.d.). This approach allows for detailed representation of individual movements and interactions (Zhou, 2012). The previously computed time-varying arrival rate function and service rate function were input. Pedestrian behavior was determined based on the proportion functions $P(t)$ and $r$. Each individual agent was removed from the simulation once they had either completed dining within the cafeteria or exited the system, thereby marking the end of their simulation lifecycle.
 
\subsubsection*{Heatmap analysis (Kernel Density Estimation, KDE)}
\begin{figure}[H]
  \centering
  \includegraphics[width=0.45\textwidth]{myplot 1.png}
  \caption{Heatmap visualization at the 200th second of the simulation}
  \label{fig:kde}
\end{figure}
 
Heat maps are obtained based on Kernel Density Estimation (KDE), whose values do not represent the absolute number of people, but rather the density of people in the vicinity per unit area (Zhou, Tang, Wang, \& Wang, 2013). Based on our feelings, a KDE value of more than 4 in some areas indicates an extreme crowding of people in this area. 
 
Through Figure~\ref{fig:kde} we can find that the crowded area is basically concentrated in a part near the service windows, while the density in the corridor is relatively less. Therefore, we can make a teacher or student volunteer to control the flow of people at the leftmost side of all windows (Figure~\ref{fig:loc}), which is where students must go through to enter the queueing area. This could be done by releasing 14 - 16 people every 2 minutes for the first 10 minutes (maximum arrival rate between 0.117 to 0.133).
 
\begin{figure}[H]
    \centering
    \includegraphics[width=0.4\linewidth]{loc.png}
    \caption{Location of the flow control}
    \label{fig:loc}
\end{figure}
 
\section*{Discussion}
 
In summary, validation of the optimization strategy shows that adding an additional service window and limiting the average service time to 25 seconds can significantly reduce waiting times during peak hours. Integrated with flow-control of releasing only 12–14 people every two minutes. This approach effectively mitigates sudden surges in customer arrival, thereby maintaining the system utilization rate, ρ, within the safe interval of 0.8–0.9 and keeping the maximum average waiting time at approximately 50 seconds. Moreover, AnyLogic dynamic simulation and heat-map analysis further reveal that crowding density peaks in the queueing area rather than in the corridor, and that installing flow-control action alongside the window can effectively alleviate the perception of congestion (Rossetti, 2021).
 
In addition to this, this study not only provides a quantifiable solution for the school, but also an effective template for other school cafeterias and even broader situations with similar characteristics and problems (e.g., varying arrival rates). Both the behavioral analysis of crowds and the queueing theory modeling theory in this study can be applied to many similar situations.
 
In contrast to existing research, the dynamic queueing-modeling framework presented here is relatively rare in practical applications; most studies rely on steady-state models for optimization, whereas this paper employs a dynamic $G(t)/M/c$ model combined with smoothing-spline functions for nonparametric estimation of service rates, thus making the model more reflective of real-world conditions.
 
Nevertheless, there remain certain sources of inaccuracies and limitations in both model design and data processing. First, the observational dataset spans only eight days and does not account for holidays, special events, or other contingencies. This sample is not representative enough. Second, student behavior is characterized in a simplified manner, only "purchasing meals" versus "not purchasing meals" and indoor versus outdoor preference. And omitting more complex actions such as "jockeying" and "reneging" (Sztrik, 2021). Furthermore, this study does not explore in depth the psychological disparity between "perceived congestion" and "actual congestion", which is crucial for enhancing students' dining experience. For example, visually seeing crowding also affects one's perception, even if there is no change in the density of the surrounding crowd (Yildirim \& Akalin-Baskaya, 2007).
 
It should also be noted that time-varying queueing models may overlook cumulative effects of customers within the system, potentially underestimating queue lengths or peak delays (Whitt, 2006). Other error sources include approximations of the arrival process, simplifications in service time distribution, and biases in manual observational records. Notably, although service times were found to conform to an exponential distribution via the Kolmogorov–Smirnov test. However, it has not been proved that there is no other more consistent distribution that can explain it.
 
Future work may proceed in the following directions: (1) expand the sample size by collecting at least 30 days of data to enhance result precision; (2) integrate surveys or physiological indicators to quantify students' perceived levels of crowding under different optimized situations to conclude an experiment (3) extend the model to more school types (e.g., middle schools, universities) and cafeterias of different scales to evaluate its generalizability.
 
 
\section*{Conclusion}
In this study, the effectiveness of queueing theory combined with behavioral simulation in alleviating canteen congestion is demonstrated, and this study propose a dynamic-modeling approach that uses time-varying data to capture evolving scenarios. Based on multi-day field observations and a $G(t)/M/c$ framework, we quantitatively reveal that system instability during the lunch rush arises when arrival rates exceed service capacity. By adding an extra service window, shortening service durations, and implementing flow-control measures, the average waiting time was effectively reduced to under 50 seconds, providing a quantitative foundation and technical roadmap for real-world campus deployment. AnyLogic dynamic simulation and heat-map analyses also further pinpoint the optimal locations for flow control.
 
Although the sample period and behavioral model remain somewhat simplified, our findings already furnish campus canteen managers with clear, quantitative guidance for operational adjustments. Future work may extend the model's generalizability and further enhance and detailed pedestrian-flow management. At both theoretical and practical levels, this study also offers a viable queueing-system analysis framework that can be applied to other school cafeterias and even broader public-service settings.
 
\bibliographystyle{apa}  % 参考文献格式
\bibliography{sample}  % 不要加 .bib 后缀
 
\clearpage
 
\section*{References}
 
\noindent
\hangindent=0.5in
\hangafter=1
Adan, I., \& Resing, J. (2002). \textit{Queueing theory}. Eindhoven University of Technology, Department of Mathematics and Computing Science. 
 
\noindent
\hangindent=0.5in
\hangafter=1
Ajiboye, A. S., \& Saminu, K. A. S. (2018). A multi-stage queue approach to solving customer congestion problem in a restaurant. \textit{Open Journal of Statistics}, 8(2), 302–316.
 
\noindent
\hangindent=0.5in
\hangafter=1
AnyLogic Company. (n.d.). \textit{Pedestrian Library}. Retrieved April 4, 2025, from \url{https://anylogic.help/library-reference-guides/pedestrian-library/index.html}
 
\noindent
\hangindent=0.5in
\hangafter=1
Chen, J., \& Wang, H. (2011). Optimization of university cafeteria using a hybrid queuing model. \textit{Journal of Huangshi Institute of Technology}, 27(3), 41–44.
 
\noindent
\hangindent=0.5in
\hangafter=1
Cooper, R. B. (1981). \textit{Introduction to Queueing Theory} (2nd ed.). North Holland.
 
\noindent
\hangindent=0.5in
\hangafter=1
Deng, S., et al. (2011). Queueing system simulation based on MATLAB for multi-server models. \textit{Journal of Anqing Teachers College (Natural Science Edition)}, 17(3), 61–63.
 
\noindent
\hangindent=0.5in
\hangafter=1
Galabo, N. R. (2019). Canteen service quality and student satisfaction. \textit{International Journal of Scientific \& Technology Research}, 8(6), 114–124.
 
\noindent
\hangindent=0.5in
\hangafter=1
Gross, D., \& Harris, C. M. (1998). \textit{Fundamentals of Queueing Theory} (3rd ed.). Wiley.
 
\noindent
\hangindent=0.5in
\hangafter=1
Kambli, A., et al. (2020). Improving campus dining operations using capacity and queue management: A simulation-based case study. \textit{Journal of Hospitality and Tourism Management}, 43, 62–70.
 
\noindent
\hangindent=0.5in
\hangafter=1
Kleinrock, L. (1975). \textit{Queueing Systems: Volume I: Theory}. Wiley.
 
\noindent
\hangindent=0.5in
\hangafter=1
Lee, K.-E., \& Shanklin, C. W. (2002). How long does it take students to eat lunch? A summary of three studies. \textit{Journal of Child Nutrition \& Management}, 26(2).
 
\noindent
\hangindent=0.5in
\hangafter=1
Li, L., et al. (2016). Performance analysis and optimal allocation of layered defense M/M/N queueing systems. \textit{Mathematical Problems in Engineering}, 2016, 1–21.
 
\noindent
\hangindent=0.5in
\hangafter=1
Li, Y. (2010). Simulation of queueing problems based on MATLAB. \textit{Journal of Wuhan University of Technology (Information \& Management Edition)}, 32(6), 892–896.
 
\noindent
\hangindent=0.5in
\hangafter=1
Lilliefors, H. W. (1969). On the Kolmogorov–Smirnov Test for the Exponential Distribution with Mean Unknown. \textit{Journal of the American Statistical Association}, 64(325), 387–389.
 
\noindent
\hangindent=0.5in
\hangafter=1
Little, J. D. C. (1961). A proof for the queuing formula: L = λ W. \textit{Operations Research}, 9(3), 383–389.
 
\noindent
\hangindent=0.5in
\hangafter=1
Lu, Y., et al. (2021). Stochastic simulation of parallel queue systems. \textit{Mathematics in Practice and Theory}, 51(4), 200–206.
 
\noindent
\hangindent=0.5in
\hangafter=1
Pollock, D. S. G. (1993). \textit{Smoothing with cubic splines} (Working Paper No. 291). Queen Mary and Westfield College, University of London.
 
\noindent
\hangindent=0.5in
\hangafter=1
Reinsch, C. H. (1967). Smoothing by spline functions. \textit{Numerische Mathematik}, 10(3), 177–183.
 
\noindent
\hangindent=0.5in
\hangafter=1
Rossetti, M. D. (2021). \textit{Simulation modeling and Arena} (2nd ed.). Wiley.
 
\noindent
\hangindent=0.5in
\hangafter=1
Sztrik, J. (2021). \textit{Queueing problems with solutions}. University of Debrecen, Faculty of Informatics.
 
\noindent
\hangindent=0.5in
\hangafter=1
Tang, T.-Q., et al. (2019). Statistical analysis and modeling of pedestrian flow in university canteen during peak period. \textit{Physica A: Statistical Mechanics and Its Applications}, 521, 29–40.
 
\noindent
\hangindent=0.5in
\hangafter=1
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., … \& van Mulbregt, P. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. \textit{Nature Methods}, 17(3), 261–272.
 
\noindent
\hangindent=0.5in
\hangafter=1
Wahba, G. (1990). \textit{Spline models for observational data} (Vol. 59). Society for Industrial and Applied Mathematics.
 
\noindent
\hangindent=0.5in
\hangafter=1
Whitt, W. (1993). Approximations for the GI/G/m queue. \textit{Production and Operations Management}, 2(2), 114–161.
 
\noindent
\hangindent=0.5in
\hangafter=1
Whitt, W. (2006). Fluid models for multiserver queues with abandonments. \textit{Operations Research}, 54(1), 37–54.
 
\noindent
\hangindent=0.5in
\hangafter=1
Yildirim, K., \& Akalin-Baskaya, A. (2007). Perceived crowding in a café/restaurant with different seating densities. \textit{Building and Environment}, 42(9), 3410–3417.
 
\noindent
\hangindent=0.5in
\hangafter=1
Ye, Z. (2009). Application of M/M/C queuing model in hairdressing service. \textit{Journal of Chongqing Normal University}, 26(2), 75–78.
 
\noindent
\hangindent=0.5in
\hangafter=1
Zuo, B. (2020). Analysis and optimization of noodle queue efficiency in college canteens. \textit{Value Engineering}, 39(32), 183–184.
 
\noindent
\hangindent=0.5in
\hangafter=1
Zhou, B., Tang, H., Wang, X., \& Wang, X. (2013). Visualizing crowd movement patterns using a directed kernel density estimation method. \textit{Journal of Visualization}, 16(3), 217–225.
 
\noindent
\hangindent=0.5in
\hangafter=1
Zhou, Y. (2012). Agent-based modeling and simulation for pedestrian movement behaviors in space: A review of applications and GIS issues. \textit{ISPRS International Journal of Geo-Information}, 1(1), 3–19.
 
 
\end{document}