1. 项目概述：当数学建模遇上排队难题

每年国赛数学建模，总有几个题目让参赛队伍又爱又恨，排队等待问题绝对是其中之一。它不像纯粹的优化或预测问题那样有明确的公式可套，其核心难点在于“随机性”——顾客到达的时间是随机的，服务时长也是随机的。你无法用一个简单的方程来精确描述整个系统的动态变化。这时候，蒙特卡洛模拟就成了我们手中的“神器”。它不跟你讲复杂的微分方程或排队论稳态解，它的哲学很简单：既然现实世界充满随机，那我就用计算机“造”一个虚拟世界，通过成千上万次的随机实验，来观察、统计并预测这个系统的行为。这就像你想知道一个复杂骰子游戏的平均收益，最笨但最有效的方法就是亲自玩上一万遍，然后算个平均数。蒙特卡洛法就是这个“玩上一万遍”的过程，只不过由计算机在瞬间完成。

对于参加数学建模竞赛的同学来说，掌握蒙特卡洛模拟解决排队问题，是一个极具性价比的技能。它思路直观，不需要过于高深的数学背景，但产出的结果却非常具有说服力——清晰的图表、具体的等待时间分布、服务台利用率，这些都能让你的论文脱颖而出。而MATLAB，凭借其强大的矩阵运算能力和丰富的绘图函数，是实现这一想法的绝佳平台。它能让你的代码简洁高效，更能一键生成那些让评委眼前一亮的可视化结果。接下来，我就以一个典型的“银行窗口服务”排队场景为例，带你从零开始，用MATLAB搭建一个完整的蒙特卡洛模拟，把抽象的“随机过程”变成屏幕上直观的动画和数据。

2. 核心思路拆解：为什么是蒙特卡洛？

在深入代码之前，我们必须先搞清楚两个核心概念：排队问题的本质，以及蒙特卡洛方法为何能成为它的“克星”。

2.1 排队问题的随机性内核

一个最简单的单服务台排队模型（M/M/1），通常由几个要素构成：顾客到达间隔时间、服务台的服务时间、排队规则（如先到先服务）。问题的核心在于，前两者往往不是定值。例如，银行顾客的到达并不是每分钟准时来一个，而是有时密集有时稀疏，这通常用泊松过程来描述，意味着到达间隔时间服从指数分布。同样，办理业务的时间也长短不一，可能服从指数分布或正态分布。

这种随机性导致系统状态（队列长度、顾客等待时间）也是随机的、动态变化的。传统的解析方法（如利特尔公式）可以给出系统在长期运行下的平均性能指标，但对于我们想了解的细节——比如“在上午高峰期，顾客平均要等多久？最长等待时间可能有多长？”或者“如果增加一个服务台，顾客等待时间超过10分钟的概率会降低多少？”——解析解往往无能为力，或者求解过程异常复杂。这时，基于离散事件仿真的蒙特卡洛模拟就显示出其优势：它不追求一个完美的数学解，而是通过模拟系统随时间推进的真实过程，来收集我们关心的任何统计数据。

2.2 蒙特卡洛模拟的工作逻辑

蒙特卡洛方法的核心是“随机抽样”和“统计估计”。应用到排队模拟中，其工作流程可以分解为以下几步：

1. 定义模型与参数：明确系统有哪些组成部分（顾客、队列、服务台），以及它们的随机规则（到达率、服务率）。
2. 初始化：设置模拟时钟为0，初始化系统状态（队列为空，服务台空闲）。
3. 事件驱动：整个模拟由一个“事件列表”驱动。主要事件有两种：“顾客到达事件”和“顾客离开事件（服务完成）”。模拟时钟总是跳到下一个最早发生的事件时间点。
4. 处理事件：
   • 处理“到达事件”：生成一个顾客，记录其到达时间。如果服务台空闲，则立即开始服务，并为该顾客生成一个“离开事件”加入事件列表；如果服务台忙，则顾客进入队列等待。
   • 处理“离开事件”：服务台变为空闲。如果队列中有顾客在等待，则队首顾客出队开始服务，并为其生成新的“离开事件”。
5. 数据收集：在每个事件处理过程中，记录关键数据，如顾客的等待时间（开始服务时间 - 到达时间）、队列长度、服务台忙闲状态等。
6. 循环与终止：重复步骤3-5，直到模拟时钟达到预设的终止时间（如模拟8小时营业），或者已处理完指定数量的顾客。
7. 统计分析：模拟结束后，对所有收集到的数据进行统计分析，计算平均值、标准差、分布直方图、95%分位数等，从而回答我们最初提出的问题。

这个流程听起来可能有点抽象，但一旦用代码实现，你就会发现它的逻辑非常清晰和强大。

3. 实战准备：MATLAB环境与模型参数设定

工欲善其事，必先利其器。在动手写模拟核心代码前，我们需要在MATLAB中做好准备工作，并明确我们要模拟的具体场景。

3.1 场景定义与参数化

假设我们要模拟一家银行的一个服务窗口（单服务台）。经过初步观察或题目给定，我们得到以下参数：

• 到达率：平均每小时到达10位顾客，即平均到达间隔时间为6分钟（60/10）。假设到达间隔时间服从指数分布。
• 服务率：窗口业务员平均每小时能处理12位顾客，即平均服务时间为5分钟（60/12）。假设服务时间也服从指数分布。
• 模拟时长：模拟银行一个工作日的工作时间，共计8小时（480分钟）。
• 排队规则：先到先服务（FCFS），队列长度理论上无限制。

我们的目标是：通过模拟，估算出顾客的平均等待时间、平均队列长度、服务窗口的利用率，以及顾客等待时间超过10分钟的概率。

在MATLAB中，我们首先将这些参数定义清楚：

% ========== 模拟参数设置 ========== clear; clc; close all; % 清空环境 lambda = 10; % 平均到达率 (顾客/小时) mu = 12; % 平均服务率 (顾客/小时) avg_interarrival_time = 60 / lambda; % 平均到达间隔时间 (分钟) avg_service_time = 60 / mu; % 平均服务时间 (分钟) total_simulation_time = 8 * 60; % 总模拟时间，8小时转换为分钟 num_customers_target = 1000; % 另一种终止条件：模拟接待1000名顾客 % 初始化随机数种子，确保结果可复现 rng(2024);

注意：这里使用了rng(2024)来固定随机数种子。这在数学建模中至关重要。它保证了每次运行程序生成的随机序列是一样的，使得你的结果可重复，便于调试和论文中展示稳定的数据。如果去掉这行，每次运行结果都会不同。

3.2 指数分布随机数的生成

蒙特卡洛模拟的“随机”来源，就是这里。在MATLAB中，生成服从指数分布的随机数非常简单。指数分布的概率密度函数为f(t) = λ * exp(-λ*t)，其中λ是率参数（单位时间内事件发生的平均次数）。对于到达过程，λ_arrival = lambda/60（因为我们的时间单位是分钟）。MATLAB的exprnd函数可以直接生成：

% 生成一个指数分布的随机间隔时间（分钟） % exprnd(mu) 生成均值为 mu 的指数分布随机数 interarrival_time = exprnd(avg_interarrival_time); service_time = exprnd(avg_service_time);

这个步骤将会在模拟循环中被反复调用，用以决定下一个顾客何时到来，以及当前顾客需要服务多久。

4. 核心模拟引擎：事件驱动的编程实现

这是整个项目最核心的部分。我们将采用“面向过程”的事件调度法来实现模拟引擎。虽然MATLAB也支持面向对象，但对于初次接触离散事件模拟的同学，过程式的写法更直观易懂。

4.1 数据结构初始化

我们需要一些“容器”来记录模拟过程中的各种状态和信息。

% ========== 初始化数据结构 ========== % 事件列表：第一行是事件时间，第二行是事件类型（1=到达， 2=离开） event_list = []; % 初始为空 % 系统状态变量 current_time = 0; % 模拟时钟 server_status = 0; % 服务台状态，0=空闲，1=繁忙 queue_length = 0; % 当前队列长度 queue_arrival_times = []; % 队列中顾客的到达时间记录（用于计算等待时间） % 统计变量 num_customers_served = 0; % 已服务顾客数 total_waiting_time = 0; % 累计等待时间 waiting_times_record = []; % 记录每个顾客的等待时间，用于后续画分布图 queue_lengths_over_time = []; % 记录随时间变化的队列长度（用于绘图） time_points = []; % 记录队列长度对应的时间点

4.2 主循环与事件处理逻辑

主循环的驱动力是“事件列表”。我们总是处理列表中时间最早的那个事件。

% ========== 生成第一个到达事件，启动模拟 ========== first_arrival_time = exprnd(avg_interarrival_time); event_list = [first_arrival_time, 1]; % 事件类型1代表到达 % ========== 主模拟循环 ========== while current_time < total_simulation_time && num_customers_served < num_customers_target % 1. 从事件列表中找出最早发生的事件 [next_event_time, idx] = min(event_list(:, 1)); event_type = event_list(idx, 2); % 2. 推进模拟时钟 current_time = next_event_time; % 3. 记录当前队列长度（用于绘图） queue_lengths_over_time(end+1) = queue_length; time_points(end+1) = current_time; % 4. 根据事件类型进行处理 switch event_type case 1 % 顾客到达事件 % 从事件列表中移除已处理的到达事件 event_list(idx, :) = []; % 为该顾客生成服务时间 this_service_time = exprnd(avg_service_time); if server_status == 0 % 服务台空闲 % 立即开始服务，无等待 server_status = 1; departure_time = current_time + this_service_time; % 生成该顾客的离开事件 event_list = [event_list; departure_time, 2]; % 记录该顾客的等待时间为0 waiting_times_record(end+1) = 0; total_waiting_time = total_waiting_time + 0; num_customers_served = num_customers_served + 1; else % 服务台繁忙 % 顾客进入队列等待 queue_length = queue_length + 1; queue_arrival_times(end+1) = current_time; % 记录其到达时间 end % 为下一个顾客生成到达事件 next_interarrival = exprnd(avg_interarrival_time); next_arrival_time = current_time + next_interarrival; event_list = [event_list; next_arrival_time, 1]; case 2 % 顾客离开事件（服务完成） % 从事件列表中移除已处理的离开事件 event_list(idx, :) = []; if queue_length > 0 % 队列中有顾客在等待 % 队首顾客出队 queue_length = queue_length - 1; customer_arrival_time = queue_arrival_times(1); queue_arrival_times(1) = []; % 从队列中移除 % 计算该顾客的等待时间 this_waiting_time = current_time - customer_arrival_time; waiting_times_record(end+1) = this_waiting_time; total_waiting_time = total_waiting_time + this_waiting_time; num_customers_served = num_customers_served + 1; % 为该顾客生成离开事件 this_service_time = exprnd(avg_service_time); departure_time = current_time + this_service_time; event_list = [event_list; departure_time, 2]; % 服务台继续保持繁忙状态 else % 队列为空 % 服务台变为空闲 server_status = 0; end end end % 模拟结束，处理可能仍在队列中的顾客（可选，这里我们简单忽略） fprintf('模拟结束。共服务了 %d 名顾客。\n', num_customers_served);

这段代码是模拟的核心引擎。它完美诠释了“事件驱动”：模拟时间不是均匀流逝的，而是跳跃到下一个事件发生点。这种设计极大地提高了模拟效率，因为我们无需在每一分每一秒都检查系统状态。

实操心得：在调试这类事件驱动模拟时，最容易出错的地方是事件列表的管理。一定要确保在switch语句的每个分支里，都正确地移除了当前正在处理的事件，并正确地添加了新生成的事件。可以用disp(event_list)在循环内打印事件列表，来直观跟踪事件的产生和消费过程。

5. 数据分析与可视化：让结果自己说话

模拟跑完了，数据也记录下来了，但一堆数字缺乏冲击力。我们需要用MATLAB强大的绘图功能，将结果直观地呈现出来。这是论文拿高分的关键。

5.1 核心性能指标计算

首先，我们计算几个最关键的指标。

% ========== 计算性能指标 ========== if num_customers_served > 0 avg_waiting_time = total_waiting_time / num_customers_served; fprintf('顾客平均等待时间: %.2f 分钟\n', avg_waiting_time); % 计算服务台利用率 (繁忙时间 / 总模拟时间) % 由于我们采用事件驱动，需要另一种方式计算。一个简单近似是： % 利用率 = (总服务时间) / (总模拟时间) % 总服务时间 ≈ 平均服务时间 * 已服务顾客数 total_service_time = avg_service_time * num_customers_served; utilization = total_service_time / current_time; fprintf('服务台利用率: %.2f%%\n', utilization * 100); % 计算等待时间超过10分钟的概率 prob_long_wait = sum(waiting_times_record > 10) / num_customers_served; fprintf('等待时间超过10分钟的概率: %.2f%%\n', prob_long_wait * 100); % 计算平均队列长度 (时间加权平均) % 通过记录的 queue_lengths_over_time 和 time_points 计算 total_queue_customer_minutes = 0; for i = 2:length(time_points) time_interval = time_points(i) - time_points(i-1); avg_queue_in_interval = (queue_lengths_over_time(i-1) + queue_lengths_over_time(i)) / 2; total_queue_customer_minutes = total_queue_customer_minutes + avg_queue_in_interval * time_interval; end avg_queue_length = total_queue_customer_minutes / current_time; fprintf('平均队列长度: %.2f 人\n', avg_queue_length); else fprintf('未服务任何顾客。\n'); end

5.2 多维度可视化绘图

一图胜千言。我们至少需要绘制三张图。

图1：顾客等待时间分布直方图这张图能直观展示等待时间的波动情况，是判断服务系统稳定性的重要依据。

% ========== 绘图1：等待时间分布直方图 ========== figure('Position', [100, 100, 800, 600]); % 设置图形窗口大小 subplot(2,2,1); histogram(waiting_times_record, 30, 'Normalization', 'probability', 'FaceColor', [0.2, 0.6, 0.8], 'EdgeColor', 'k'); hold on; % 在图上标注平均等待时间 ylimits = ylim; line([avg_waiting_time, avg_waiting_time], [0, ylimits(2)], 'Color', 'r', 'LineWidth', 2, 'LineStyle', '--'); text(avg_waiting_time*1.05, ylimits(2)*0.9, sprintf('平均: %.1f min', avg_waiting_time), 'Color', 'r', 'FontWeight', 'bold'); xlabel('等待时间 (分钟)'); ylabel('概率'); title('顾客等待时间分布'); grid on;

图2：队列长度随时间变化图这张图动态展示了系统拥堵情况的变化，能清晰看出高峰期和低谷期。

% ========== 绘图2：队列长度随时间变化 ========== subplot(2,2,2); stairs(time_points, queue_lengths_over_time, 'b-', 'LineWidth', 1.5); xlabel('模拟时间 (分钟)'); ylabel('队列长度 (人)'); title('队列长度动态变化'); grid on; % 标注平均队列长度 ylimits = ylim; line([0, time_points(end)], [avg_queue_length, avg_queue_length], 'Color', 'r', 'LineWidth', 1.5, 'LineStyle', '--'); text(time_points(end)*0.7, avg_queue_length*1.1, sprintf('平均: %.2f', avg_queue_length), 'Color', 'r');

图3：服务台忙闲状态片段图选取一小段时间，展示服务台“忙”与“闲”的交替过程，非常直观。

% ========== 绘图3：服务台状态片段 (示例) ========== subplot(2,2,3); % 我们需要从事件日志中重构状态，这里为了简化，我们模拟最后100分钟的状态 sample_start = max(0, current_time - 100); sample_end = current_time; % 创建一个时间向量 time_vec = sample_start:0.5:sample_end; % 每0.5分钟采样一次 status_vec = zeros(size(time_vec)); % 这里需要根据事件列表粗略判断状态，实际项目中需要更精细的记录。 % 此处用一段示例代码展示思路，假设我们记录了每个繁忙期的开始和结束。 % 绘图一个简单的方波示意 % 假设我们“虚构”几个繁忙时段用于演示绘图方法 busy_periods = [sample_start+10, sample_start+25; sample_start+40, sample_start+70; sample_start+85, sample_end-5]; for i = 1:size(busy_periods,1) idx = time_vec >= busy_periods(i,1) & time_vec <= busy_periods(i,2); status_vec(idx) = 1; end stairs(time_vec, status_vec, 'r-', 'LineWidth', 2); ylim([-0.2, 1.5]); yticks([0, 1]); yticklabels({'空闲 (0)', '繁忙 (1)'}); xlabel('时间 (分钟)'); ylabel('服务台状态'); title('服务台忙闲状态 (片段)'); grid on;

图4：关键指标总结表在图中嵌入一个表格，让评委一眼看到核心结果。

% ========== 绘图4：关键指标文本总结 ========== subplot(2,2,4); axis off; % 不显示坐标轴 text(0.1, 0.9, '模拟结果摘要', 'FontSize', 14, 'FontWeight', 'bold'); text(0.1, 0.7, sprintf('总模拟时间: %.0f 分钟', current_time), 'FontSize', 11); text(0.1, 0.6, sprintf('服务顾客总数: %d 人', num_customers_served), 'FontSize', 11); text(0.1, 0.5, sprintf('平均等待时间: %.2f 分钟', avg_waiting_time), 'FontSize', 11); text(0.1, 0.4, sprintf('平均队列长度: %.2f 人', avg_queue_length), 'FontSize', 11); text(0.1, 0.3, sprintf('服务台利用率: %.1f%%', utilization*100), 'FontSize', 11); text(0.1, 0.2, sprintf('长等待(>10min)概率: %.1f%%', prob_long_wait*100), 'FontSize', 11);

将这四个子图组合在一张图上，形成一份完整的分析报告，视觉效果和专业性都会大大提升。

6. 模型扩展与灵敏度分析：提升论文深度

如果只完成基本模拟，论文可能止步于“良好”。要想冲击更高奖项，必须进行模型扩展和灵敏度分析，展示你对问题的深入思考。

6.1 扩展一：多服务台（M/M/c）模型

现实中的银行往往有多个窗口。将我们的单服务台模型扩展为多服务台模型，是逻辑上的自然延伸。主要修改点在于：

1. server_status从一个标量变为一个向量或计数器，记录每个服务台的状态。
2. 处理“到达事件”时，需要遍历所有服务台，找到第一个空闲的。如果都忙，则进入一个公共的队列（这是最常见模型）。
3. 处理“离开事件”时，释放对应的服务台，然后检查公共队列。

% 多服务台模型核心修改示例 num_servers = 3; % 假设有3个服务台 server_status = zeros(1, num_servers); % 0表示空闲 % 在到达事件中寻找空闲服务台 free_server = find(server_status == 0, 1); if ~isempty(free_server) % 有空闲服务台，分配 server_status(free_server) = 1; % ... 生成该服务台上的离开事件 ... else % 所有服务台忙，进入公共队列 queue_length = queue_length + 1; queue_arrival_times(end+1) = current_time; end % 在离开事件中，释放特定服务台后，检查公共队列 if queue_length > 0 % 从队列中取出一个顾客 % 分配该顾客给刚刚空闲的服务台 % ... 生成新的离开事件 ... end

通过对比单服务台和多服务台下的平均等待时间、队列长度等指标，可以定量分析增加服务资源的效益，这是建模中经典的“成本-效益”分析。

6.2 扩展二：非指数分布的服务时间

现实中，服务时间可能更符合正态分布（大部分时间集中在均值附近）或均匀分布。修改模型非常容易，只需替换生成service_time的随机数函数。

% 指数分布 (原模型) service_time = exprnd(avg_service_time); % 改为正态分布 (需指定标准差，并避免负值) std_service_time = 1.5; % 假设标准差为1.5分钟 service_time = normrnd(avg_service_time, std_service_time); service_time = max(0.1, service_time); % 防止出现负值或零值 % 改为均匀分布 (例如在[3,7]分钟之间) min_time = 3; max_time = 7; service_time = unifrnd(min_time, max_time);

比较不同分布假设下的结果，可以分析模型对输入分布的“稳健性”。如果结果差异很大，说明你需要花更多精力去实地调研，获取真实的服务时间数据分布。

6.3 灵敏度分析：改变关键参数

这是数学建模论文的“加分神器”。系统地改变一个或两个关键参数（如到达率lambda或服务台数量c），观察输出指标（如平均等待时间）如何变化。

% 灵敏度分析示例：分析到达率对平均等待时间的影响 lambda_range = 8:0.5:14; % 测试从8到14（顾客/小时）的到达率 avg_wait_results = zeros(size(lambda_range)); for i = 1:length(lambda_range) lambda_test = lambda_range(i); % 将之前的模拟代码封装成一个函数 simulate_queue(lambda, mu, ...) % 这里调用该函数，并返回平均等待时间 % avg_wait_results(i) = simulate_queue(lambda_test, mu, ...); end % 绘图 figure; plot(lambda_range, avg_wait_results, 'bo-', 'LineWidth', 2, 'MarkerSize', 8); xlabel('顾客到达率 \lambda (人/小时)'); ylabel('平均等待时间 (分钟)'); title('系统性能灵敏度分析：到达率 vs 平均等待时间'); grid on; hold on; % 可以标注当前设计值（lambda=10）的点 idx = find(lambda_range == 10); if ~isempty(idx) plot(lambda_range(idx), avg_wait_results(idx), 'r*', 'MarkerSize', 15, 'LineWidth', 2); text(lambda_range(idx), avg_wait_results(idx)*1.05, '设计点 (\lambda=10)', 'Color', 'r', 'FontWeight', 'bold'); end

这张图能清晰地展示系统性能随负载变化的趋势。当到达率接近或超过服务率时，平均等待时间会急剧上升（系统趋于不稳定），这个结论非常有力。你还可以进一步绘制“等待时间超过10分钟的概率”随到达率变化的曲线。

7. 常见问题与调试技巧实录

在实际编写和运行模拟程序时，你肯定会遇到各种问题。下面是我在多次实践中总结的一些典型“坑”和解决方法。

7.1 程序陷入死循环或运行极慢

• 可能原因1：事件列表管理错误。最常见的是没有正确移除已处理的事件，导致同一个事件被反复处理，或者生成了时间戳为Inf或NaN的事件。

  • 排查：在while循环内加入调试语句，打印current_time,event_list和event_type。观察事件是否被正常消费和添加。
  • 技巧：使用unique函数或确保事件列表按时间排序，可以避免重复事件。使用isfinite()检查生成的时间是否为有效数字。

• 可能原因2：终止条件不满足。如果使用“模拟固定时长”的条件，但顾客到达事件被不断推后（比如生成了一个极大的间隔时间），可能导致时钟无法推进到终止时间。

  • 解决：同时设置双重终止条件，如while current_time < T_max && num_customers_served < N_max。N_max作为一个安全上限。

7.2 统计结果明显不合理

• 可能原因1：时间单位混淆。这是新手最常犯的错误。到达率lambda是“人/小时”，服务时间均值是“分钟”，在生成随机数时如果单位没统一，结果会完全错误。

  • 检查：从头检查所有涉及时间的变量和参数，确保它们在同一个单位制下（建议全程使用“分钟”）。
  • 技巧：在程序开头用注释明确标出每个时间变量的单位。

• 可能原因2：预热期数据污染。模拟开始时，系统是空的，需要一段时间才能达到稳定状态。如果统计从一开始就计算，会拉低平均队列长度和利用率。

  • 解决：设置一个“预热期”（Warm-up Period），例如前1000分钟或前100个顾客的数据不纳入最终统计。在代码中，当current_time > warmup_time后再开始记录数据。

• 可能原因3：随机数种子问题。单次模拟的结果具有偶然性。一次运行可能恰好运气好（或差），不能代表系统普遍性能。

  • 解决：进行多次独立重复实验。将整个模拟过程（从rng之后到统计计算）放在一个for循环中，循环N次（如N=100），每次使用不同的随机数种子（例如rng('shuffle')或rng(i)）。最后对所有N次运行的结果取平均值和置信区间。这才是蒙特卡洛模拟的正确打开方式。

num_replications = 100; avg_wait_reps = zeros(1, num_replications); for rep = 1:num_replications rng(rep); % 为每次重复实验设置不同的种子 % 这里是完整的单次模拟代码... % 将单次模拟得到的 avg_waiting_time 存入数组 % avg_wait_reps(rep) = avg_waiting_time_from_this_simulation; end final_avg_wait = mean(avg_wait_reps); wait_std = std(avg_wait_reps); confidence_interval = [final_avg_wait - 1.96*wait_std/sqrt(num_replications), ... final_avg_wait + 1.96*wait_std/sqrt(num_replications)]; fprintf('基于%d次独立实验，平均等待时间为%.2f分钟，95%%置信区间为[%.2f, %.2f]。\n', ... num_replications, final_avg_wait, confidence_interval(1), confidence_interval(2));

7.3 可视化图形不美观或信息不全

• 问题：图形挤在一起，标签看不清，没有图例。
• 技巧：
  1. 使用figure('Position', [x, y, width, height])调整图形窗口大小。
  2. 使用subplot合理布局多张图。
  3. 为线条、柱状图添加清晰的DisplayName，并使用legend('show', 'Location', 'best')添加图例。
  4. 为坐标轴添加单位xlabel('时间 (分钟)')。
  5. 使用grid on开启网格，提高可读性。
  6. 使用title为每张图起一个描述性的标题。
  7. 重要的参考线（如平均值、阈值）用不同颜色和线型标出，并用text函数添加标注。

最后，将所有这些代码模块（参数设置、模拟引擎、数据分析、可视化、灵敏度分析）整合到一个或多个.m脚本文件中，加上清晰的注释和章节标题，这本身就是一份高质量的、可运行的数学建模程序报告。在论文中，你可以直接截取关键代码片段、生成的图表和结论分析。记住，评委看重的是你利用计算工具解决实际问题的思路和能力，而蒙特卡洛模拟正是展示这种能力的绝佳舞台。