蒙特卡洛模拟——CAE不确定性量化的基本方法
理论与物理
概述 — 什么是蒙特卡洛方法
老师，蒙特卡洛方法就是用随机数反复模拟吗？要做多少次才行呢？
问得好。粗略地说，这是一种从输入参数的概率分布中随机采样，并重复执行数百次、数千次FEM或CFD分析的方法。目的是从统计上求得输出的波动（分布）。
几千次！？那计算成本不是很高吗？
是的，这也是蒙特卡洛方法最大的弱点。但其最大的优点是，无论输入变量是多少维，收敛速度都不会改变。灵敏度分析中，如果有10个变量，成本会呈指数级增长，但蒙特卡洛方法无论是10维还是100维，都以相同的速度收敛。
与维度无关… 这太厉害了。具体在什么时候使用呢？
例如，评估汽车零部件材料强度的波动（变异系数CoV约5%）如何影响疲劳寿命的情况。当杨氏模量、屈服应力、板厚、载荷幅值等多个参数同时波动时，研究这些“组合效应”如何传播到输出，是蒙特卡洛方法擅长的领域。
蒙特卡洛模拟（Monte Carlo Simulation, MCS）是一种从概率性输入变量出发，重复执行确定性模型（FEM、CFD等），以估计输出统计特性（均值、方差、分布形状、置信区间）的方法。其名称来源于以赌场闻名的摩纳哥蒙特卡洛地区，因其使用随机数而得名。
基本原理与收敛速度
收敛速度是怎么决定的？刚才说“与维度无关”…
蒙特卡洛方法的核心是大数定律和中心极限定理。从输入分布中独立采样 $N$ 个样本 $\mathbf{X}_1, \mathbf{X}_2, \ldots, \mathbf{X}_N$，并对每个样本评估模型 $Y = f(\mathbf{X})$。则样本均值为：
这个样本均值的标准误（Standard Error）由下式表示：
啊，除以 $\sqrt{N}$ 啊！也就是说，增加样本数误差会变小，但只是按平方根的比例减小…
没错。这就是著名的$1/\sqrt{N}$ 收敛法则。重要的是：
- 要使误差减半，需要将样本数增至4倍
- 要使误差降至1/10，需要100倍的样本
- 这个收敛速度与输入变量的维数 $d$ 无关（不受维度诅咒影响）
样本方差的估计量同样定义如下：
置信区间的估计
把结果写到报告里时，只说“平均值大概是这么多”不行吧？置信区间怎么求呢？
多亏了中心极限定理，只要 $N$ 足够大（实践中大约从 $N \geq 30$ 开始），样本均值就可以近似为正态分布。所以95%置信区间可以这样写：
1.96是正态分布的值对吧。具体来说，做1000次的话精度大概是多少？
我们来计算一下相对估计精度（Relative Accuracy）。设变异系数为 $\text{CoV}_Y = \sigma_Y / \mu_Y$：
例如，当输出的变异系数为 $\text{CoV}_Y = 0.10$（10%）时：
- $N = 100$: $\text{RA} = 1.96 \times 0.10 / \sqrt{100} = 1.96\%$
- $N = 1{,}000$: $\text{RA} = 1.96 \times 0.10 / \sqrt{1000} \approx 0.62\%$
- $N = 10{,}000$: $\text{RA} \approx 0.20\%$
1000次时，95%置信区间的宽度约为平均值的6%（当 $\text{CoV}_Y \approx 1.0$ 时）是实践中的一个大致标准。
样本数的确定方法
那实践中大概要运行多少次呢？是100次就行，还是需要10,000次…
这完全取决于目的。大致总结如下：
| 目的 | 所需样本数的大致标准 | 理由 |
|---|---|---|
| 均值・标准差的估计 | 100〜1,000 | $1/\sqrt{N}$ 法则已足够收敛 |
| 95百分位数的估计 | 1,000〜5,000 | 为了准确捕捉分布的尾部 |
| 99.9百分位数（稀有事件） | 10,000〜100,000 | 尾部精度需要大量样本 |
| 故障概率 $P_f = 10^{-k}$ 的估计 | $\geq 10^{k+2}$ | 为使 CoV($\hat{P}_f$) < 10% 所需 |
故障概率是 $10^{-6}$ 的话，需要 $10^8$ 次！？那根本跑不完吧…
所以对于稀有事件的估计，不会使用纯粹的蒙特卡洛方法，而是使用重要性采样（Importance Sampling）或子集模拟（Subset Simulation）等特殊方法。或者引入代理模型来降低计算成本。
收敛速度的理论背景
- 大数定律：当 $N \to \infty$ 时，样本均值 $\bar{Y}$ 依概率收敛于真实期望值 $\mu_Y = E[Y]$。这保证了蒙特卡洛估计的一致性。
- 中心极限定理：$\sqrt{N}(\bar{Y} - \mu_Y) / \sigma_Y \to \mathcal{N}(0,1)$。只要存在有限方差，无论原始分布形状如何，样本均值都渐近服从正态分布。这是构建置信区间的理论依据。
- $1/\sqrt{N}$ 法则的含义：误差以 $O(N^{-1/2})$ 的速度减小。与基于网格的数值积分（如梯形法则）的收敛速度 $O(N^{-2/d})$ 相比，在低维（$d \leq 3$）时较慢，但在高维（$d \geq 5$）时蒙特卡洛方法更有优势。
采样方法与方差缩减
纯随机采样（Crude Monte Carlo）
刚才公式里用的“随机采样”，就是简单地生成随机数就行了吗？
最简单的方法是这样的。从每个输入变量的概率分布中独立生成伪随机数来创建样本。这被称为纯随机采样（Crude Monte Carlo / Simple Random Sampling）。
但是因为是随机数，有时会出现偏差吧？比如样本集中在某个区域，或者出现空白地带…
很敏锐。特别是当样本数较少时（大约 $N < 100$），纯随机采样容易导致输入空间的覆盖不均匀。例如，对2个变量生成100个样本时，分布区域的角落可能完全没有样本。这时就出现了LHS（拉丁超立方采样）。
拉丁超立方采样（LHS）
LHS这个名字我听说过。和普通随机数有什么区别？
LHS的想法很简单。将每个输入变量的分布分割成 $N$ 个等概率的层（stratum），确保从每一层中恰好取一个样本点。然后对各变量间的组合进行随机打乱。
啊，原来如此！就是均匀地“品尝”每个变量的分布对吧。就像数独一样，每行每列都必须恰好有一个？
这个理解是对的。数学上，将每个变量 $X_j$（$j = 1, \ldots, d$）的累积分布函数 $F_j$ 的值域 $[0, 1]$ 进行 $N$ 等分：
其中 $\pi_j$ 是 $\{1, 2, \ldots, N\}$ 的一个随机排列，$U_j^{(i)} \sim \text{Uniform}(0, 1)$ 用于随机化层内的位置。
总结一下LHS的优点：
- 可以用较少的样本均匀覆盖输入空间
- 每个变量的边缘分布能被准确再现
- 与纯随机采样相比，方差更小（但改善幅度取决于具体问题）
- 在CAE实践中，通常 $N = 50\text{〜}200$ 左右的LHS就能获得足够的精度
| 比较项目 | 纯随机采样（SRS） | LHS |
|---|---|---|
| 空间覆盖 | 有偏差 | 均匀 |
| 边缘分布再现 | $N \to \infty$ 时保证 | 每层保证 |
| 方差缩减效果 | 无（基准） | 取决于问题，但通常有利 |
| 相关性控制 | 困难 | 使用带相关性控制的LHS可实现 |
| 追加样本的便利性 | 容易 | 需要重新设计 |
方差缩减法（Variance Reduction Techniques）
除了LHS，还有其他方法可以用较少的样本数提高精度吗？
有一类被称为“方差缩减法”的技术。其原理是有效减小标准误 $\text{SE} = \sigma_Y / \sqrt{N}$ 中的 $\sigma_Y$，或者通过改进估计量，在相同的 $N$ 下提高精度。
| 方法 | 原理 | 适用场景 |
|---|---|---|
| 对偶变量法 (Antithetic Variates) | 使用样本 $\mathbf{U}$ 和 $1-\mathbf{U}$ 的配对，引入负相关 | 单调的响应函数 |
| 控制变量法 (Control Variates) | 利用期望值已知的辅助变量修正估计量 | 存在近似解析解的情况 |
| 重要性采样 (Importance Sampling) | 在关注区域密集采样，并用似然比修正 | 稀有事件・故障概率估计 |
| 分层采样 (Stratified Sampling) | 将输入空间分层，在各层内采样 | 不均匀的响应函数 |
对偶变量法看起来很简单。几乎不增加成本就能减小方差吗？
如果响应函数对输入是单调的（例如：材料强度越高寿命越长），效果会很明显。当 $\mathbf{U}$ 产生较大的输出时，$1-\mathbf{U}$ 会产生较小的输出，因此取两者的平均值会抵消方差。但要注意，对于非单调函数，有时可能适得其反。
重要性采样的数学公式化
考虑故障概率 $P_f = P[g(\mathbf{X}) \leq 0]$ 的估计。在原始密度函数 $f_{\mathbf{X}}$ 下：
$$P_f = \int I[g(\mathbf{x}) \leq 0] \, f_{\mathbf{X}}(\mathbf{x}) \, d\mathbf{x} = \int I[g(\mathbf{x}) \leq 0] \frac{f_{\mathbf{X}}(\mathbf{x})}{h(\mathbf{x})} h(\mathbf{x}) \, d\mathbf{x}$$从提议分布 $h(\mathbf{x})$ 中采样，并用似然比 $w(\mathbf{x}) = f_{\mathbf{X}}(\mathbf{x}) / h(\mathbf{x})$ 进行加权：
$$\hat{P}_f^{IS} = \frac{1}{N}\sum_{i=1}^{N} I[g(\mathbf{X}_i) \leq 0] \cdot w(\mathbf{X}_i), \quad \mathbf{X}_i \sim h$$最优的提议分布是集中在故障区域的分布，通常采用在FORM/SORM设计点附近配置正态分布的方法。
实践指南
实践流程
在实际用CAE分析运行蒙特卡洛方法时，步骤是怎样的？
基本流程如下：
- 识别不确定的输入参数 — 材料特性、形状尺寸、载荷条件等
- 设定各参数的概率分布 — 正态分布、对数正态分布、均匀分布等
- 制定采样计划 — 纯随机采样 or LHS，决定样本数 $N$
- 执行 $N$ 次CAE分析 — 批处理 or 并行执行
- 结果的统计处理 — 计算均值、标准差、直方图、CDF、置信区间
- 收敛性检查 — 通过绘制累积均值・累积标准差的图表来确认稳定性
输入分布的设定
输入参数的概率分布是怎么决定的？可以随便设为正态分布吗？
这是最令人头疼的地方，“垃圾进，垃圾出（Garbage In, Garbage Out）”在输入分布设定上体现得最为关键。理想情况下，应根据制造商的品质数据或实验数据进行拟合，但在实践中，通常使用以下标准：
| 参数 | 典型分布 | 变异系数（CoV）的大致标准 |
|---|---|---|
| 杨氏模量 | 正态分布 | 2〜5% |
| 屈服应力 | 对数正态分布 | 5〜10% |
| 板厚・尺寸 | 正态分布 | 1〜3%（假设公差的1/6为$\sigma$） |
| 疲劳强度 | 对数正态 or 威布尔 | 10〜30% |
| 载荷幅值 | 对数正态 or 冈贝尔 | 10〜20% |
| 温度条件 | 正态分布 or 均 |
なった
詳しく
報告