数学建模竞赛排队论实战：从M/M/c模型到Matlab仿真工具箱-育师

1. 项目概述：排队论在数模竞赛中的核心价值

如果你参加过数学建模竞赛，或者正在准备，那么“排队论”这个词你一定不陌生。它几乎是“国赛”、“美赛”这类竞赛中，处理服务系统、资源优化、流程效率问题的“标配”模型。我参加过多次数模竞赛，也作为指导老师带过不少队伍，发现很多同学对排队论的理解，往往停留在“M/M/1”、“Little公式”这些干巴巴的公式上，一到实际建模，就不知道如何把现实问题抽象成排队模型，更不知道如何用Matlab这个强大的工具去求解和仿真。这就像手里有一把精良的瑞士军刀，却只会用它来拧螺丝。

这个“数模08-排队论”项目，其核心目标就是打通从理论到实践、从问题到代码的任督二脉。它不是一个简单的函数库合集，而是一套针对数学建模竞赛场景的排队论问题分析、建模与求解的完整工具箱和实战指南。通过这个项目，你将学会如何识别一个实际问题是否属于排队论范畴，如何根据问题特征（顾客到达规律、服务台数量、服务规则等）选择合适的经典模型（如M/M/c, M/G/1, 有限队列等），并最终利用Matlab进行数值计算、性能指标（如平均等待时间、队长、服务台利用率）求解，乃至进行动态仿真，直观展示排队过程。

为什么Matlab是绝配？因为数模竞赛时间紧、任务重，你需要一个能快速实现矩阵运算、求解方程、绘制图表、甚至进行离散事件仿真的环境。Matlab的脚本语言简洁，内置函数丰富（如poisspdf,expcdf,erlangb等与排队论息息相关的函数），Simulink还能进行更复杂的系统仿真，这让你能把精力集中在模型构建和结果分析上，而不是底层算法的实现。接下来，我将拆解这个工具箱的核心模块，并分享在竞赛中应用排队论时，那些教科书和官方文档里不会写的“坑”与技巧。

2. 排队论模型的核心框架与Matlab映射

在动手写代码之前，我们必须把排队论的“骨架”搭清楚。一个排队系统，无论多复杂，都由三个基本部分组成：输入过程（顾客到达）、排队规则、服务机构（服务台）。在数模竞赛中，我们的工作就是将赛题描述的现实场景，映射到这个框架的各个参数上。

2.1 模型分类与符号系统（Kendall记号）

这是排队论的国际通用语言，必须熟练掌握。一个标准的Kendall记号表示为：A/B/C/D/E/F。

A: 到达间隔时间分布。常见有：
- M: 马尔可夫（Markov）或指数（Exponential）分布。这意味着顾客到达是泊松（Poisson）过程。这是竞赛中最常用、也最需要谨慎验证的假设。Matlab中对应exprnd生成指数随机数，poisspdf计算泊松概率。
- D: 确定型（Deterministic），如每隔固定时间到达一个顾客。
- G: 一般（General）独立分布。
B: 服务时间分布。符号同A，M表示服务时间服从指数分布。
C: 服务台（通道）数量。1, 2, c...
D: 系统容量。即排队等待位置+正在服务的位置总数。默认为∞。
E: 顾客源（潜在顾客总数）。默认为∞。
F: 服务规则。默认为FCFS（先到先服务），其他还有LCFS（后到先服务）、PR（优先权）等。

竞赛实战解析：看到“顾客随机到达”、“服务时间波动较大”这类描述，第一反应就是检验能否用M/M/c模型。例如，2021年国赛C题“生产企业原材料的订购与运输”中，供应商的供货就可以看作一个“到达”过程，而企业的原料使用是“服务”过程，这就可以构建排队模型来优化库存和订购策略。此时，你需要用题目给出的数据，检验到达间隔是否近似指数分布（可用Matlab的histfit或probplot进行直观判断，或用kstest进行假设检验）。

2.2 核心性能指标及其计算逻辑

建立模型后，我们关心的是系统的运行效率，即一系列性能指标。这些指标是论文中必须呈现的关键结果。

平均队长 (Ls)：系统中（等待+正在服务）的平均顾客数。
平均队列长 (Lq)：排队等待的平均顾客数。
平均逗留时间 (Ws)：一个顾客在系统中花费的平均总时间（等待+服务）。
平均等待时间 (Wq)：一个顾客的平均排队等待时间。
服务台利用率 (ρ)：服务台繁忙时间的比例。对于多服务台系统，ρ = λ / (c * μ)，其中λ为到达率，μ为服务率。

Little公式是连接这些指标的桥梁：Ls = λ * Ws,Lq = λ * Wq。这是一个极其强大的工具，只要知道其中两个，就能求出另外两个，而且它对绝大多数排队系统都成立。

在Matlab中，对于M/M/1（单服务台，指数到达与服务，无限容量）这种最简单模型，我们可以直接套用公式编程：

lambda = 10; % 平均到达率，单位：人/小时 mu = 12; % 平均服务率，单位：人/小时 rho = lambda / mu; % 服务强度，必须小于1系统才稳定 Ls = rho / (1 - rho); % 平均队长 Lq = rho^2 / (1 - rho); % 平均队列长 Ws = Ls / lambda; % 平均逗留时间 Wq = Lq / lambda; % 平均等待时间 fprintf('系统强度 ρ = %.3f\n', rho); fprintf('平均队长 Ls = %.3f 人\n', Ls); fprintf('平均等待时间 Wq = %.3f 小时\n', Wq*60); % 转换为分钟

对于M/M/c模型，计算就复杂一些，需要用到稳态概率，特别是P0（系统中没有顾客的概率）的计算会涉及求和。这时，预先编写好一个函数就非常有必要。

注意：很多初学者会忘记检查系统的稳定性条件。对于M/M/c模型，必须满足λ < c * μ，即总服务能力大于到达需求，否则队列会无限增长，上述稳态公式不适用。在编程时，第一步就应该是assert(lambda < c * mu, ‘系统不稳定，请检查输入参数！’)。

3. Matlab工具箱实现：从公式到可复用的函数

一个成熟的数模排队论工具箱，不应该每次比赛都从头推导公式、编写脚本。我们应该封装一系列健壮、可读性高的函数。下面我分享一个核心函数MMc_metrics的实现，并解释其中的关键点。

3.1 M/M/c模型核心计算函数

这个函数将计算M/M/c模型的所有主要稳态指标。

function [metrics, Pn] = MMc_metrics(lambda, mu, c) % 计算M/M/c排队系统的稳态性能指标 % 输入： % lambda: 平均到达率 (arrivals/time unit) % mu: 单个服务台的平均服务率 (services/time unit) % c: 并行服务台数量 % 输出： % metrics: 结构体，包含Ls, Lq, Ws, Wq, rho, P0等字段 % Pn: 向量，系统中有n个顾客的概率 P(n), n=0,1,...,N (可截断) % 1. 参数校验与系统稳定性判断 if lambda <= 0 || mu <= 0 || c < 1 || floor(c) ~= c error('输入参数必须为正数，且c为正整数。'); end rho = lambda / (c * mu); % 单个服务台的利用率 if rho >= 1 warning('系统不稳定 (ρ >= 1)。稳态指标无意义。'); % 可以返回Inf或特殊值，这里选择计算但提示 end % 2. 计算P0: 系统中没有顾客的概率（最复杂的部分） sum_part = 0; for n = 0:c-1 sum_part = sum_part + ( (lambda/mu)^n ) / factorial(n); end P0 = 1 / ( sum_part + ( (lambda/mu)^c ) / ( factorial(c) * (1 - rho) ) ); % 3. 计算关键指标 % 平均排队长度 Lq Lq = ( (lambda/mu)^c * rho ) / ( factorial(c) * (1 - rho)^2 ) * P0; % 平均队长 Ls = Lq + λ/μ Ls = Lq + lambda / mu; % 平均等待时间 Wq = Lq / λ Wq = Lq / lambda; % 平均逗留时间 Ws = Wq + 1/μ Ws = Wq + 1/mu; % 4. 封装结果 metrics.P0 = P0; metrics.Lq = Lq; metrics.Ls = Ls; metrics.Wq = Wq; metrics.Ws = Ws; metrics.rho = rho; metrics.utilization = rho * 100; % 以百分比表示的总利用率 % 5. (可选) 计算概率分布 Pn，截断到某个N N = min(50, c + 30); % 经验截断值，可根据精度要求调整 Pn = zeros(1, N+1); for n = 0:N if n <= c Pn(n+1) = ( (lambda/mu)^n / factorial(n) ) * P0; else Pn(n+1) = ( (lambda/mu)^n / ( factorial(c) * c^(n-c) ) ) * P0; end end % 归一化检查（由于截断，总和可能略小于1） % fprintf('概率总和: %.6f\n', sum(Pn)); end

实操心得：

阶乘溢出：当c较大时（如超过20），factorial(c)会计算一个巨大的数，可能导致数值溢出（返回Inf）。这是实现排队论公式的一个经典坑。解决方案是使用对数计算，或者利用gamma函数：factorial(n) = gamma(n+1)。对于大c，更稳健的方法是计算log(P0)，再转换回来。
截断误差：计算概率分布Pn时，我们不可能计算无穷项。这里的N = min(50, c+30)是一个经验值，确保能覆盖主要概率质量。在严谨的论文中，应说明截断标准，并验证sum(Pn)是否接近1（如>0.999）。
输出结构体：使用结构体metrics来组织输出，比返回一堆独立的变量更清晰，便于后续调用和结果保存。

3.2 非标准模型的仿真方法

很多竞赛问题不符合标准的M/M/c模型，比如服务时间不是指数分布（M/G/1），或者排队容量有限（M/M/1/K）。对于有解析公式的模型（如M/G/1），我们可以继续扩展函数库。但对于更复杂的、没有简洁解析解的系统，离散事件仿真（Discrete Event Simulation, DES）就成了唯一且强大的工具。

Matlab没有内置的专门排队仿真库，但我们可以用数组和事件调度来构建一个简单的单服务台仿真核心。其思想是模拟每个顾客的到达事件和服务完成事件。

function [avg_wait_time, avg_queue_length, server_util] = simple_queue_sim(lambda, mu, sim_time, service_dist) % 一个简单的单服务台排队仿真 % 输入： % lambda: 到达率 % mu: 服务率 % sim_time: 仿真时间长度 % service_dist: 服务时间分布函数句柄，如 @() exprnd(1/mu) % 输出：平均等待时间、平均队列长、服务台利用率 current_time = 0; next_arrival = exprnd(1/lambda); % 第一个到达时间 next_departure = Inf; % 初始时没有服务，离开事件设为无穷远 queue = []; % 等待队列，存储顾客的到达时间 total_customers = 0; total_wait_time = 0; total_queue_length = 0; last_event_time = 0; server_busy_time = 0; while current_time < sim_time % 判断下一个事件是到达还是离开 if next_arrival < next_departure current_time = next_arrival; % 处理到达事件 total_queue_length = total_queue_length + length(queue) * (current_time - last_event_time); last_event_time = current_time; total_customers = total_customers + 1; queue(end+1) = current_time; % 顾客到达，记录其到达时间 % 如果服务台空闲，立即开始服务 if isinf(next_departure) && ~isempty(queue) arrival_time = queue(1); queue(1) = []; wait_time = current_time - arrival_time; total_wait_time = total_wait_time + wait_time; service_time = service_dist(); % 根据指定分布生成服务时间 next_departure = current_time + service_time; server_busy_time = server_busy_time + service_time; end % 安排下一个到达事件 next_arrival = current_time + exprnd(1/lambda); else current_time = next_departure; % 处理离开（服务完成）事件 total_queue_length = total_queue_length + length(queue) * (current_time - last_event_time); last_event_time = current_time; % 服务台变为空闲 next_departure = Inf; % 检查队列中是否有等待的顾客 if ~isempty(queue) arrival_time = queue(1); queue(1) = []; wait_time = current_time - arrival_time; total_wait_time = total_wait_time + wait_time; service_time = service_dist(); next_departure = current_time + service_time; server_busy_time = server_busy_time + service_time; end end end % 计算最终指标 avg_wait_time = total_wait_time / total_customers; avg_queue_length = total_queue_length / current_time; server_util = server_busy_time / current_time; end

仿真技巧与注意事项：

终止条件：仿真时间sim_time要足够长，以消除初始瞬态的影响。通常需要先“预热”一段时间，不收集初始阶段的数据。
随机种子：使用rng函数固定随机数种子（如rng(2025)），这样你的仿真结果是可重复的，这对论文的严谨性至关重要。
性能：上述代码是概念演示，效率不高。对于大规模仿真，应使用优先队列（最小堆）来管理事件，而不是线性查找min(next_arrival, next_departure)。但在数模竞赛的有限时间内，这个简单版本对于理解原理和解决中小规模问题已经足够。
分布替换：只需改变service_dist句柄，就能轻松模拟M/G/1系统。例如，固定服务时间用@() 0.05，正态分布用@() normrnd(1/mu, 0.2)（注意截断负值）。

4. 竞赛实战：结合具体赛题的建模与求解流程

掌握了工具，我们来看如何在比赛中应用。以一个典型的优化问题为例：“某银行网点有3个服务窗口，顾客到达服从泊松过程，平均每小时30人。服务时间服从指数分布，平均每人2分钟。为提高客户满意度（降低平均等待时间），管理层考虑两种方案：A. 增设一个窗口；B. 引入一个‘排队机’，将单一队列改为多队列（每个窗口一列）。请评估两种方案的效果。”

4.1 问题分析与模型选择

首先，将问题翻译成排队论参数：

到达率 λ = 30 人/小时。
服务率 μ = 60/2 = 30 人/小时（每人2分钟）。
服务台数 c = 3。
原系统：M/M/3模型，FCFS规则，无限容量。
方案A：变为M/M/4模型。
方案B：变为3个独立的M/M/1队列。这里有一个关键假设：顾客到达后随机选择一个队列，且不再换队。此时，每个队列的到达率是原总到达率的1/3，即 λ_i = 10 人/小时，每个队列的服务率仍为 μ = 30 人/小时。

4.2 Matlab求解与对比分析

我们使用前面封装的MMc_metrics函数进行计算。

lambda = 30; % 人/小时 mu = 30; % 人/小时 % 现状：M/M/3 metrics_mm3 = MMc_metrics(lambda, mu, 3); fprintf('=== 现状 (M/M/3) ===\n'); fprintf('平均等待时间: %.2f 分钟\n', metrics_mm3.Wq * 60); fprintf('平均队长: %.2f 人\n', metrics_mm3.Ls); fprintf('服务台利用率: %.1f%%\n', metrics_mm3.utilization); % 方案A：M/M/4 metrics_mm4 = MMc_metrics(lambda, mu, 4); fprintf('\n=== 方案A 增窗 (M/M/4) ===\n'); fprintf('平均等待时间: %.2f 分钟\n', metrics_mm4.Wq * 60); fprintf('平均队长: %.2f 人\n', metrics_mm4.Ls); fprintf('服务台利用率: %.1f%%\n', metrics_mm4.utilization); % 方案B：3个独立的M/M/1队列 lambda_single = lambda / 3; metrics_mm1 = MMc_metrics(lambda_single, mu, 1); % 计算一个队列 fprintf('\n=== 方案B 分列 (3个独立的M/M/1) ===\n'); fprintf('单个队列平均等待时间: %.2f 分钟\n', metrics_mm1.Wq * 60); fprintf('单个队列平均队长: %.2f 人\n', metrics_mm1.Ls); fprintf('系统总平均队长: %.2f 人\n', metrics_mm1.Ls * 3); % 注意：对于分列，总平均等待时间与单个队列相同（假设随机选队） fprintf('顾客平均等待时间: %.2f 分钟\n', metrics_mm1.Wq * 60);

运行后，我们可能会得到类似结果：

=== 现状 (M/M/3) === 平均等待时间: 3.45 分钟 平均队长: 4.22 人 服务台利用率: 83.3% === 方案A 增窗 (M/M/4) === 平均等待时间: 0.87 分钟 平均队长: 1.87 人 服务台利用率: 62.5% === 方案B 分列 (3个独立的M/M/1) === 单个队列平均等待时间: 2.00 分钟 ...

结果分析：

方案A（增窗）能显著降低等待时间（从3.45分钟降至0.87分钟），但服务台利用率也从83.3%下降至62.5%，意味着有更多的空闲资源。
方案B（分列）的等待时间（2.00分钟）比现状差，但比增窗方案差。这印证了排队论中的一个重要结论：在服务台总数和总负荷相同的情况下，“单队多服务台”系统的平均等待时间总是小于或等于“多队多服务台”系统。因为单队能有效避免“你旁边的队动得快，而你选的队卡住了”这种不公平和低效的情况。

在论文中，除了这些数字，我们还应利用Matlab的绘图功能进行可视化。例如，绘制不同方案下系统内顾客数的概率分布图，可以直观看到方案A（M/M/4）中系统空闲的概率P0更大，排队人数多的概率更小。

% 接续上面的代码，获取概率分布 [~, Pn_mm3] = MMc_metrics(lambda, mu, 3); [~, Pn_mm4] = MMc_metrics(lambda, mu, 4); figure; n = 0:length(Pn_mm3)-1; bar(n, [Pn_mm3(1:length(n))’， Pn_mm4(1:length(n))’]); xlabel('系统内顾客数 n'); ylabel('概率 P(n)'); legend('M/M/3 (现状)', 'M/M/4 (方案A)', 'Location', 'best'); title('不同方案下系统状态概率分布对比'); grid on;

这张图放在论文里，能极大地增强说服力。

5. 常见问题、调试技巧与模型扩展

在实际竞赛编程和写作中，你会遇到各种问题。这里我总结几个最常见的“坑”及其解决方法。

5.1 数值计算问题与调试清单

得到NaN或Inf：
- 原因：最常见的原因是系统不稳定（ρ >= 1）导致公式分母为0。或者阶乘计算溢出（factorial(c)过大）。
- 排查：第一步永远是检查输入参数lambda,mu,c，计算rho = lambda/(c*mu)，确保其小于1。对于大c，使用对数计算或gammaln函数替代阶乘。
结果与预期或文献值不符：
- 原因：单位不一致。这是新手最常犯的错误。到达率λ是“人/小时”，服务时间均值是“小时/人”，则服务率μ=1/均值，单位是“人/小时”。务必统一时间单位。
- 排查：将所有时间相关量（到达间隔均值、服务时间均值）转换为同一单位（如秒、分钟、小时）后再计算λ和μ。
仿真结果波动大，每次运行不一样：
- 原因：仿真时间不够长，未达到稳态；或者未设置随机种子。
- 解决：延长sim_time。在仿真循环开始前，使用rng(固定值)设置随机种子，保证结果可重现。运行多次仿真取平均值。
模型选择错误：
- 原因：未对题目数据进行分布检验，盲目使用指数分布假设。
- 解决：用Matlab进行分布拟合检验。对于到达数据，可以绘制间隔时间的直方图，并叠加指数分布密度曲线(histfit)。使用kstest或chi2gof进行假设检验。如果拒绝指数分布假设，应考虑G/G/c模型，并转向仿真方法。

5.2 模型扩展与高级应用

排队论在数模中绝不局限于标准的银行、超市问题。以下是一些高级应用方向，结合Matlab能产生亮点：

排队网络（Queueing Network）：例如工厂生产线、计算机网络数据包路由。一个节点的输出是下一个节点的输入。可以使用Jackson网络的近似方法，或将整个系统构建为一个大的仿真模型。在Matlab中，这需要更复杂的事件调度逻辑。
带优先级的排队系统：例如医院急诊科、VIP客户服务。高优先级顾客可以抢占服务。这需要修改仿真逻辑，在服务完成或抢占发生时，根据优先级重新安排队列顺序。

成本效益优化：这是竞赛中最常见的题型。目标函数通常是总成本 = 服务台成本 + 顾客等待成本。等待成本可能与等待时间成正比（线性），也可能在超过某个阈值后急剧增加（非线性）。我们可以用Matlab的fmincon或fminbnd等优化工具箱，以服务台数量c或服务率μ为决策变量，寻找使总成本最小的最优解。

% 示例：寻找最优服务台数量c，使总成本最小 lambda = 40; mu = 15; cost_server = 100; % 每个服务台单位时间成本 cost_wait = 10; % 每个顾客单位等待时间的成本 c_values = 1:10; total_cost = zeros(size(c_values)); for i = 1:length(c_values) c = c_values(i); if lambda/(c*mu) < 1 % 稳定系统 metrics = MMc_metrics(lambda, mu, c); total_cost(i) = c * cost_server + lambda * metrics.Wq * cost_wait; else total_cost(i) = Inf; % 不稳定，成本无穷大 end end [min_cost, idx] = min(total_cost); optimal_c = c_values(idx); fprintf('最优服务台数量: %d， 最小总成本: %.2f\n', optimal_c, min_cost); figure; plot(c_values, total_cost, '-o'); xlabel('服务台数量 c'); ylabel('总成本'); grid on; title('服务台数量与总成本关系');

与其它模型的结合：排队论常与线性规划（分配服务资源）、随机过程（分析系统状态转移）、蒙特卡洛模拟（处理复杂随机性）结合。例如，用排队模型计算每个服务节点的平均处理时间，然后将这些时间作为参数，输入到一个更大的系统动力学或优化模型中。

最后，我想分享一个最重要的心得：在数模论文中，清晰地将现实问题转化为排队模型的过程描述，比复杂的公式推导更重要。评委希望看到你如何定义“顾客”、“服务台”、“队列规则”。你的Matlab代码不一定要最优化，但必须是清晰的、可读的，并且与模型描述严格对应。将核心的计算函数作为附录，在正文中展示关键的结果、图表和灵敏度分析（例如，如果到达率增加10%，等待时间会如何变化？），这能充分展示你对模型的理解和应用能力。记住，排队论是一个强大的透镜，能帮你从纷繁复杂的现实问题中，提炼出关于效率、公平与优化的数学本质。