第十五讲 计算机仿真

数学建模 讲义 NUDT 2023SP

---

---

15.1 什么是计算机仿真

• 计算机仿真 （Compter Simulation），又称为计算机模拟，是数学建模的过程，运行在计算机之上，旨在预测某个真实世界或物理系统的行为或者结果
• 很多时候，理论（模型）的可靠性可以通过对比仿真结果和真实世界结果的一致性得出
• 计算机仿真，对于很多物理学（计算物理学）、天体物理学、气候学、化学、生物学、制造业中的自然系统，以及经济学、心理学、社学科学、医疗保健、工程学中的人类系统的数学建模，都是非常有帮助的工具

---

真实世界 vs 数学模型

• A computer model is the algorithms and equations used to capture the behavior of the system being modeled.
• By contrast, computer simulation is the actual running of the program that contains these equations or algorithms.
• Simulation, therefore, is the process of running a model.
• Thus one would not "build a simulation"; instead, one would "build a model (or a simulator)", and then either "run the model" or equivalently "run a simulation".

---

认识世界的第三种方法

• 计算机仿真的本质是利用计算机模拟现实情况或系统演变过程，发现或验证新的信息和规律
• 有人认为，除了传统的理论研究和实验研究外，计算机缝针是认识世界的第三种方法
• 计算机仿真的优势
  
  • 代价小：相对于物理世界中的实验验证，代价（资源消耗、时间消耗）往往小得多
  • 效率高：一般能够在短时间内给出结果，加速科学探索的过程
  • 可重复性强：重复实现的成本低，实现容易
  • 设置灵活：参数设置灵活，结果便于比较

---

什么问题适于使用计算机仿真求解

• 特征一：随机性
  
  • 随机现象的模拟或随机试验
  • 例如：离散事件仿真、Monte Carlo 方法
• 特征二：涌现性
  
  • 局部已知，全局未知：局部信息描述相对完备、规则清晰，但全局规律难以预测和判断
  • 例如：个体数量众多导致整体规律不明确
  • 例如：系统单步（或短时间内）变化规律已知，但演变所需时间较长，导致系统最终状态或发展方向难以预测
  • 例如：复杂系统模拟、天气模拟、动力系统仿真

---

计算机仿真的类型

1. 物理规律仿真：一般只具备特征二，且与时间无关
   
   • 适用于物理规律重现、物理系统优化设计等，因为常常涉及连续系统，也有人将其称为连续系统仿真
2. 系统演变/动力系统仿真：强调特征二，且与时间有关
   
   • 通过模拟研究系统状态、特征等随时间的变化规律，如：天气模拟、生命游戏
3. Monte Carlo 方法：只具备特征一
   
   • 利用随机方法求解特定的问题，如：Buffon 投针实验，随机（智能）优化
4. 离散事件仿真 (DES，Discrete Event Simulation)：同时具备两个特征
   
   • 基于事件驱动的模式，对系统进行模拟或优化，如：服务排队，资源调度优化

---

15.2 物理规律仿真

• 系统模拟是研究物理系统的重要方法
• 规模巨大、结构复杂的系统，往往系统的整体模型可能无法建立，但对其局部进行建模研究却是可行的
• 整体的解析解可能无法求得，但得到局部的解析解却相对容易

---

车灯线光源的优化设计

• 安装在汽车头部的车灯的形状为一旋转抛物面，车灯的对称轴水平地指向正前方, 其开口半径 36 毫米，深度 21.6 毫米
• 经过车灯的焦点，在与对称轴相垂直的水平方向，对称地放置一定长度的均匀分布的线光源
• 要求在某一设计规范标准下确定线光源的长度

---

设计规范（简化版）

• 在焦点 $F$ 正前方 25 米处的 $A$ 点放置一测试屏，屏与 $FA$ 垂直，用以测试车灯的反射光
• 在屏上过 $A$ 点引出一条与地面相平行的直线，在该直线 $A$ 点的同侧取 $B$ 点和 $C$ 点，使 $AC=2AB=2.6$ 米
• 要求 $C$ 点的光强度不小于某一额定值（可取为 1 个单位）， $B$ 点的光强度不小于该额定值的两倍（只须考虑一次
  反射）

---

任务

1. 在满足该设计规范的条件下，计算线光源长度，使线光源的功率最小
2. 对得到的线光源长度，在有标尺的坐标系中画出测试屏上反射光的亮区
3. 讨论该设计规范的合理性

---

分析

• 这是一个确定性的问题，没有随机因素，且与时间无关
• 对于线光源，理论推导得出亮区边界的难度较大
• 可以使用“光线追踪”模型，通过计算机模拟，绘出亮区图形

---

光线追踪

• 整个线光源经过反射在屏上形成亮区的位置和范围不明确
• 光源所发出的任意一条光线，反射后在屏上的位置可以计算
• 如果遍历“所有光线”即可确定亮区
• 忽略直射光线的影响

---

光线追踪模型框架

1. 建立合适的坐标系
2. 计算线光源上每一点所发出的每一条光线在屏上的像

3. 绘出亮区

   • 该框架比较粗糙，要实现为代码，必须考虑具体的细节，主要是离散化的问题
   • 可能需要离散化的对象：测试屏、线光源、入射光线、反射面

---

离散化

• 线光源的离散化：均匀离散化，记线光源上一点为 $s_i$
• 遍历“所有”的入射光线
  
  • 方法一：方向离散化，对于方向 $\vec a$，与点 $s_i$ 一起根据直线点法式方程，确定入射光线方程.
    
    • 特点：光强度分布均匀，但计算量大
  • 方法二：反射面离散化，对于反射面上一点 $p_j$，与点 $s_i$ 一起根据两点式方程，确定入射光线方程
    
    • 特点：计算简便，但强度不容易保持均匀
• 相对而言，第二种方法更可行，将光线的离散化问题，转化为反射面的离散化

---

反射面的离散化

• 反射面的离散化也即对曲面的剖分
• 尽可能保证剖分的“均匀”
  
  • 方法一：面积均匀
  • 方法二：曲率均匀
• 注意剖分的方法与粒度对精度影响
• 使用不同的坐标系，对计算也有一定影响

---

仿真结果

线光源长 3.89mm，柱坐标. 离散步长：线光源 0.005，半径 0.01，极角 0.005，每像素代表 4mm$\times$4mm.

---

物理规律仿真需要注意的问题

• 假设的合理性：建立描述系统或系统某部分的模型
• 问题的分解：是将未知的、待解决的“全局”问题，分解为已知的、可计算的”局部”信息
  
  • ”局部信息”不同，仿真的复杂度也不同，进而影响问题求解的质量
• 注意初始条件和仿真参数对结果的影响

---

15.2 系统演变仿真——生命游戏与细胞自动机

• 生命游戏（Game of Life, 简称 生命）是 John Conway 于 1970 年发明的
• 生命游戏不是通常意义上的游戏，没有游戏者，也无所谓胜负
• 生命游戏是一类特殊的 细胞自动机（Cellular Automata）
• 一旦给定初始状态，之后的发展完全由规则确定
• 生命充满了悬念！绝大多数情况下，不可能只根据初始状态（或称模式）判断未来的发展，只能按照游戏的规则运行下去

---

生命游戏的基本设置

• 游戏在一个方形网格中进行，网格在各个方向上都是无限延伸的
• 每个网格中有一个细胞，细胞可以是 活 的或者 死 的
• 如果细胞是活的就在它所在的方格做一个标记，死的细胞则留空
• 每个细胞的 邻居 是指它最邻近的 8 个方格

---

生命游戏的规则

• 游戏的每一步中，计算每个细胞的 活邻居数目，由此数字决定下一步发展
• 如果一个死细胞正好有 3 个活邻居，则下一步变成一个活细胞（诞生）
• 如果一个活细胞有 2 个或 3 个活邻居，则可以继续活下去（生存）
• 其他情况下，细胞将死去或保持死的状态（拥挤或孤独）

---

生命游戏的应用与发展

• 生命游戏的演变规则近似的描述了生物群体的生存规律，经过扩展能够模拟生命活动中的生存、灭绝、竞争等等复杂现象
• Conway 证明，生命游戏具有通用图灵机的计算能力，能力上与图灵机等价
• 生命游戏的主要扩展方式包括修改或扩展规则和改变生存空间设置

---

细胞自动机

• 细胞自动机 (Cellular Automata)也称元胞自动机，点格自动机，分子自动机或单元自动机
• 由数学家 Stanislaw M. Ulam (1909-1984)和 John von Neumann (1903-1957) 于上世纪 50 年代为模拟生物细胞的自我复制而提出，他们最初的问题是：能储存程序的机器能不能自我复制？
• 细胞自动机最初没有受到学术界的重视，直到 1970 年 Conway 设计出生命游戏后才吸引了科学家们的注意

---

细胞自动机的要素

• 细胞自动机本质上是一类离散型的动力系统（Discrete Dynamical Systems），关注的是系统随时间的演化.
  
  • 每个细胞代表一个个体，个体有明确的行为规则，但系统整体的规律难以预测.
• 设计一个细胞自动机需要明确以下要素：
  
  • 确定细胞活动空间的维度
  • 定义细胞可能具有的状态
  • 定义细胞改变状态的规则
  • 设定各个细胞的初始状态

---

细胞自动机形式上的特征

就形式而言，细胞自动机有三个特征：

• 平行计算（parallel computation）：每个细胞体都同时同步的改变
• 局部性（local）：细胞的状态变化只受其周围细胞的影响
• 一致性（homogeneous）：所有细胞均由同样的规则支配

---

2007 MCM-B: 登机问题

• 航空公司可以引领候机乘客以任何次序就座
• 惯例是首先引领有特殊需要的乘客就座，然后是头等舱的乘客就座（他们坐在飞机的前部）
• 然后引领持经济舱和商务舱机票的乘客从飞机后排开始向前依次就座。
• 从航空公司的角度来看，除了考虑乘客的等候时间外，时间就是金钱，所以总的登机时间最好减到最少
• 飞机只有在飞行的时候才能为航空公司赚钱，而长的登机时间限制了一架飞机一天中可以飞行的次数
• 诸如空中客车 A380（可容纳 800 名乘客）的大型机的发展就更要强调缩短登机（以及下机）时间

---

任务与分析

• 任务：就乘客人数不同的飞机：小型机（85-210），中型机（210-330）和大型机（450-800），设计登机和下机时的步骤并进行比较
• 最主要的目标是尽量减少登机时间，需要解决两方面的矛盾
  
  • 一是机舱空间的局限性与乘客相对无序性之间的矛盾，这也是自由登机费时的主要原因
  • 二是严格的管理与乘客接受程度之间的矛盾，这导致航空公司无法严格规定乘客登机顺序
• 要减少登机时间，还需要重点考虑两个冲突，走廊冲突与座位冲突
• 题目限制条件很少，也没有具体数据，开放性强，不确定因素多，难以建立简单、准确的解析模型
• 问题背景简单，容易理解，容易下手，但很多直观的结论难以量化和深入

---

解决方案示例

• 首先分析两种冲突耗费的时间，对登机时间进行估计，得到公式：$T=\frac{n}{v}+\frac{25}{6}t_S\frac{n}{3}+\frac{n-1}{2}t_L$
  
  • 其中 $n$ 为乘客数，$v$ 为单位时间内进入机舱的乘客数，乘客为了给靠窗的乘客让路，从座位起身到过道以及坐回座位的时间认为是相等的，记为 $t_S$，放行李的平均时间为 $t_L$.
  • 式中第一项为无冲突时的自由登机时间，第二项为座位冲突所耗费的时间，第三项为走廊冲突所耗费的时间
• 然后，使用细胞自动机方法对登机过程进行模拟，所得到的数值结果验证了上述近似公式的合理性
• 进一步，通过仿真计算了各种策略下的登机时间，在此基础上，提出了按行李分布控制（LDC）登机顺序的策略，很好的回答了题目所提出的问题

---

15.4 Monte Carlo 方法

• Monte Carlo 方法又称 Monte Carlo 模拟，或统计试验方法，或随机模拟
• 基本思想：将各种随机事件的概率特征（概率分布、数学期望）与随机事件的模拟联系起来，用试验的方法得到概率模型的解，而不是计算出来
• Monte Carlo 是摩纳哥国的世界著名赌城，二战时期，Von Neumann 和 Ulam 将他们从事的与研究原子弹有关的秘密工作，以赌城 Monte Carlo 作为秘密代号
• 他们的具体工作是对裂变物质的中子随机扩散进行模拟

---

Buffon 投针试验

---

积分的计算

问题： 给定函数 $g(x)$，求 $\int_a^bg(x)\mathrm{d}x$ 的近似值.

• 考虑独立同分布的随机变量 $X_i\sim U(a,b),\;i=1,2,...,n.$
• 由大数定律，$\frac{1}{n}\sum\limits_{i=1}^ng(X_i)\stackrel{P}{\longrightarrow}E(g(X_i))=\int_a^bg(x)\frac{1}{b-a}\mathrm{d}x$
• $\bbox[yellow]{\int_a^bg(x)\mathrm{d}x=(b-a)E(g(X_i))=(b-a)\lim\limits_{n\to\infty}\frac{1}{n}\sum\limits_{i=1}^ng(X_i).}$

---

例 求积分 $\int_0^1\exp(x^2)\mathrm{d}x$ 的近似值.

• 结果 $\approx 1.46265$.

---

不可召回的秘书招聘问题的改进

• 如果策略改为，拒绝前 $k$ 人，从第 $k+1$ 人开始，一旦有比前面次优人选更优秀的人选，就录取
• 那么最佳的 $k$ 又是多少？

---

不可召回的秘书招聘问题的改进

• 如果有两个条件相差较大的企业同时招聘，假设好企业仍然采用 37%法则，那么差企业应该如何应对？

---

Monte Carlo 方法解题的一般过程

1. 构造问题的概率模型
   
   • 随机性问题：描述和模拟概率过程，建立概率模型或判别式
   • 确定性问题：需转化为随机性问题，然后构造概率模型
2. 从已知概率分布中抽样
   
   • 产生已知分布的随机数序列，对随机事件进行模拟
3. 建立所需的统计量
   
   • 最常见的是用试验得到的 $k/n$ 作为解的估值
   • 若 $k/n$ 的期望值恰好是所求问题的解，则为无偏估计，有时也用极大似然估计、渐近估计等

---

Monte Carlo 方法的精度

• 一般来说，Monte Carlo 方法的估值精度 $\varepsilon$ 与试验次数 $n$ 的平方根成反比，即 $\varepsilon \propto 1/\sqrt{n}$.
  
  • 若精度 $\varepsilon$ 提高 10 倍，则试验次数 $n$ 需要增加 100 倍，求解的时间要慢100 倍
  • 收敛速度慢是 Monte Carlo 方法的主要缺点
• 当 $\varepsilon$ 一定时，试验次数 $n$ 取决于方差的数值，即 $n \propto S^2$. 因此降低方差是加速 Monte Carlo 方法收敛的主要途径
• Monte Carlo 方法的精度估计具有概率性质，它并不能断言精度一定小于 $\varepsilon$, 而只是表明，计算精度以接近于 1 的概率不超过 $\varepsilon$.

---

15.5 离散事件仿真

• 离散事件仿真通常强调以事件驱动的演化方式，即特定事件的发生是导致系统状态变化的主要因素
  
  • 系统状态的迁移、变化发生在一系列离散的事件点上
• 从时间或空间上看，系统状态是跳跃式变化的，一般不强调时空上的连续型
• 例如：生产调度管理，驱动事件可能包括：新的生产需求出现、新的生产资料的到达、新的故障的发生等
• 例如：计算机通讯网络，驱动事件可能包括：新的访问请求、新的网络故障、新增加的链路等

---

事件驱动法

• 离散事件系统不断地获取新的事件信息（外界输入或随机生成）
• 系统状态的变化不固定时间间隔，而是随着事件的发生而递进，即系统是由事件驱动
• 在设计仿真程序时，称为 事件驱动法
• 仿真系统只在有事件发生时有响应，记录参数变化，修改系统状态，或改变系统行为

---

时间步长法

• 如果系统本身的业务逻辑有一定的周期性（比如以天、月为周期进行计费），或是同一时间段内有多个事件发生，每个事件不需要马上处理，则可考虑合并处理（比如入、出院）
• 时间步长法：仿真系统以固定的时间步长推进，每推进一步，就处理该时间周期内的所有事件
• 同一系统往往可以使用不同的方法进行仿真，但可能效率不同
• 对于复杂的系统，有可能需要同时使用以上事件驱动与时间步长两种方法

---

典型的离散事件系统

• 各种排队系统（银行、电话客服、超市结帐、挂号、买火车票）
• 电梯调度
• 路口红绿灯
• 操作系统
• 网站服务系统

---

CUMCM-2009B: 眼科病床的合理安排

• 某医院眼科门诊住院部共有病床 79 张
• 该医院眼科手术主要分四大类：白内障、视网膜疾病、青光眼和外伤
• 周一、三做白内障手术，双眼手术病人，周一做一只，周三做另一只
• 外伤通常属于急症，入院第二天安排手术，其他病人入院后 2-3天即可安排手术，但一般不安排在周一、三

---

问题与需求

• 题目给出了 2008 年 7 月 13 日至 2008 年 9 月 11 日这段时间里各类病人的情况
  
  • 病人总数大于病床总数，排队不可避免
  • 但由于各类病人人数不同，术前准备时间、手术时间、术后观察时间均不同，因此如果按照先到先服务的方式安排入院顺序，不能很好利用病床资源，导致等待入院的病人越来越多，排队时间越来越长
• 题目要求设计更为合理的安排入院顺序的策略，提高病床的利用率，
  减少排队等待时间

---

分析

通过分析题目所给数据可以发现

• 各类病人的到达各自服从 Poisson 分布
• 系统服务时间没有明显规律，但术后住院时间各自服从正态分布
• 病人在入院后手术前的无效等待时间过长是导致系统低效的主要原因
• 定义病床的有效利用率来衡量系统的效率：$\mbox{病床有效利用率}=1-\frac{\mbox{无效等待时间}}{\mbox{实际住院时间}}$

---

病床的有效利用率

• 病人和医院，都希望这个指标越大越好
• 该指标使得住院安排要考虑各类病人的手术时间约束
  
  • 例如，单眼白内障病人的最佳入院日为周日、二；双眼白内障病人为周日，其他住院病人的最佳入院日为周二、三、四、五、日
• 服务时间各不相同的分类多通道排队问题
• 更复杂的，同一通道的服务时间还随病人入院时间不同而不同
• 直接应用排队论方法解决此问题将会非常困难，可以考虑使用仿真方法求解

---

调度策略

• 调度策略对等待病人的优先级排序是关键
• 常见的做法是以病人当日入院的 公平性指标（以到达先后计）与病床使用 效率指标（分类考虑）两方面综合排序（例如求两个指标的加权和），然后按排序结果安排当日入院病人，由此得到公平合理的住院方案

---

具体做法

一周中每天入院病人在手术前的无效等待时间如下表：

$$\boxed{\begin{array}{c|ccccccc} \text{类型} & \text{周一} & \text{周二} & \text{周三} & \text{周四} & \text{周五} & \text{周六} & \text{周日} \\ \hline \text{白内障} & 1 & 0 & 4 & 3 & 2 & 1 & 0 \\ \text{白内障双眼} & 6 & 5 & 4 & 3 & 2 & 1 & 0 \\ \text{青光眼} & 1 & 0 & 0 & 0 & 0 & 1 & 0 \\ \text{视网膜疾病} & 1 & 0 & 0 & 0 & 0 & 1 & 0 \end{array}}$$

利用公式 $y=(10-x)/10$，将天数 $x$ 转化为排序的效率优先值 $y$：

$$\boxed{\begin{array}{c|ccccccc} \text{类型} & \text{周一} & \text{周二} & \text{周三} & \text{周四} & \text{周五} & \text{周六} & \text{周日} \\ \hline \text{白内障} & 0.9 & 1 & 0.6 & 0.7 & 0.8 & 0.9 & 1 \\ \text{白内障双眼} & 0.4 & 0.5 & 0.6 & 0.7 & 0.8 & 0.9 & 1 \\ \text{青光眼} & 0.9 & 1 & 1 & 1 & 1 & 0.9 & 1 \\ \text{视网膜疾病} & 0.9 & 1 & 1 & 1 & 1 & 0.9 & 1 \end{array}}$$

---

优先级排序

• 同一天相同病症的病人看作一个整体，对不同病种进行排序
• 相同病症内部仍按 FCFS (先到先服务) 原则排序
• 公平指标值：当前队列中最早到门诊的病人权值为 1，接下来每晚到一天，权值减 0.1
• 综合考虑以上两级指标，排序公式定义为：$\mbox{排序值}=\lambda\times\mbox{效率指标值}+(1-\lambda)\times \mbox{公平指标值}$
• 排序值高的将排在队列前面

---

其他排序策略

• 操作系统进程调度中的最高优先比响应（HRRN）算法，使公平性与效率同时得到了体现，也是一种好的创意
• 本问题仿真中较普遍的缺点是罗列了很多排队论中的指标，但却没有抓住病床使用效率这个最重要的指标

---

15.6 总结：需要注意的几个问题

• 计算机仿真的角色及其与理论分析的关系
  
  • 有时是主角，有时是配角；因问题而异，因人而异
  • 第 1、3 类问题有可能推导出理论结果，仿真作为辅助手段，检验理论结果的正确性
  • 第 2、4 类问题通常很难推导出理论结果，这时仿真就是解决问题的主要工具，而理论分析一般可以估计出相关的上下界，极限值等，评估仿真结果是否合理

---

• 假设的合理性
  
  • 现实情况是复杂的，必须在仔细分析的基础上提出合理的假设
  • 保留主要矛盾，忽略次要因素，不一定越详细越好，也不一定考虑的因素越多就越好
• 随机量分布规律的确定与随机数的生成
  
  • 从实际背景出发选择合适的分布类型
  • 一般的高级语言提供生成 0～1 之间服从均匀分布的随机数的函数，据此产生服从特定分布的随机数
  • 更高级的专用语言，比如 Matlab，还提供生成服从正态分布等其他常用分布的随机数的函数

---

• 离散化：空间的离散化，时间的离散化
• 模拟次数的确定
  
  • 次数太少不足以说明问题，模拟次数太多，时间、空间消耗相应增加
  • 多数情况下还应对不同模拟次数的结果进行分析
• 仿真流程
  
  • 可结合问题背景进行适当优化
• 初始状态对模拟结果的影响
  
  • 通常要尽量消除这种影响，或者说要做到对初值不敏感

---

• 仿真的结果
  
  • 无论哪种类型的仿真，都要对仿真结果进行认真分析，要善于从结果中总结规律，如果有必要，可能还需要调整策略或参数，再仿真，再分析
  • 仿真只是手段，不是目的
  • 仿真为决策提供依据，仿真验证决策的效果，最终要达到发现新规律，从而改进决策的目的

---

课后思考

（$\star\star$）模拟电梯的运行. 以下为基本假设和调度策略，可以根据需要添加其他假设或对现有策略进行修改：

• 电梯最初从一楼开始运行
• 如果没有呼叫，则停在原地
• 乘客先下后上，如果满员，不能超载
• 如果电梯已经满员，且中间层无乘客出电梯，则直达目的层

请分别模拟只有一部电梯和有两部电梯（联动）的系统，并讨论如何优化其调度策略.