蒙特卡洛模拟 — 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次就够还是要10000次…
要根据目的完全不同。大致这样划分：
| 目的 | 必要样本数的目安 | 理由 |
|---|---|---|
| 平均值标准差的推定 | 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）这类特殊方法。或者用代理模型来降低计算成本。
抽样方法与方差低减
纯随机抽样（Crude Monte Carlo）
那个"随机抽样"就是随便生成随机数就行吗？
最简单的方法就是这样。从各个输入变量的概率分布中独立生成伪随机数来构造样本。这叫纯随机抽样（Crude Monte Carlo / Simple Random Sampling）。
但是随机数会出现偏差吧？ 某个地方样本集中，有的地方是空白…
说得好。特别是样本较少时（$N < 100$ 左右），纯随机数会导致输入空间覆盖不均。比如2个变量、100个样本的话，分布域的角落可能一个样本都没有。这时候就用拉丁超方格抽样（LHS）。
拉丁超方格抽样（LHS）
LHS这个名字听过。和普通随机有什么不同？
LHS的想法很简单。把各个输入变量的分布分成 $N$ 个等概率的层（stratum），从各层各取1点。然后各变量之间进行随机排列。
啊，明白了！ 对各个变量的分布进行均匀的"品尝"，像数独一样，每行每列都有一个的感觉？
想象得很好。数学上，各个变量 $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}$ 就小，两个平均一下方差就互相抵消了。但对于非单调函数可能反而有害，要注意。
蒙特卡洛的实务应用
实务流程
实际在CAE中使用蒙特卡洛法的时候，流程是怎样的？
基本流程是这样：
- 确定不确定的输入参数 — 材料特性、形状尺寸、载荷条件等
- 设定各参数的概率分布 — 正态分布、对数正态分布、均匀分布等
- 制定抽样计划 — 纯随机还是LHS，样本数 $N$ 的确定
- 执行 $N$ 次CAE分析 — 批处理或并行执行
- 结果的统计处理 — 计算平均值、标准差、直方图、CDF、置信区间
- 收敛性检验 — 用累积平均和累积标准差的图表确认稳定性
输入分布的设定
输入参数的概率分布怎么决定？ 随便用正态分布可以吗？
这是最棘手的地方。"垃圾进垃圾出（Garbage In, Garbage Out）"这句话在输入分布的设定上最适用。理想是用生产厂家的品质数据或实验数据来拟合，但实务上往往用以下参考值：
| 参数 | 典型分布 | 变动系数（CoV）的参考值 |
|---|---|---|
| 杨氏模量 | 正态分布 | 2〜5% |
| 屈服应力 | 对数正态分布 | 5〜10% |
| 板厚尺寸 | 正态分布 | 1〜3%（以公差的1/6作为σ） |
| 疲劳强度 | 对数正态或韦伯尔 | 10〜30% |
| 载荷振幅 | 对数正态或甘贝尔 | 10〜20% |
| 温度条件 | 正态分布或均匀分布 | 情况依赖 |
屈服应力用对数正态啊。为什么不用正态？
强度物理上不能是负值。正态分布的话 $\mu - 3\sigma$ 有可能变成负值，但对数正态分布永远是正值。而且强度数据一般尾部向右倾斜（高值侧拉伸），对数正态的拟合更好。
收敛判定的实际情况
"充分收敛了"怎样判断？ 一开始就决定好 $N$ 吗？
有两种方法。第一种是事先决定型，从目标精度反推必要样本数：
这里 $\epsilon$ 是容许的相对误差。比如在95%信置度下要求平均值相对误差在2%以内，预计 $\text{CoV}_Y \approx 0.10$ 时：
第二种是逐次型。一边做分析一边监看累积平均值和累积标准差，稳定了就停止。实务上一般检查以下项：
- 累积平均值的图表是否已平稳
- 最近50个样本追加后平均值是否波动在1%以内
- 用自助法（Bootstrap）算出的置信区间宽度是否在目标以下
应用事例 — 材料波动和疲劳寿命
具体应用事例能讲一下吗？ 最初说的"材料强度波动对疲劳寿命的影响"怎么做？
举个汽车悬挂臂部件的例子。确定性分析得出"疲劳寿命 $N_f = 2 \times 10^6$ 循环是安全的"。但量产品有材料波动吧？
是的。样本值是平均的，实际的部件有可能更弱…
这就是蒙特卡洛法的用处。这样做：
- 设定输入变量：
- 杨氏模量 $E$: 正态分布，$\mu = 210$ GPa，CoV = 3%
- 屈服应力 $\sigma_y$: 对数正态分布，$\mu = 350$ MPa，CoV = 5%
- S-N曲线疲劳强度系数: 对数正态分布，CoV = 15%
- 板厚 $t$: 正态分布，$\mu = 4.0$ mm，CoV = 2%
- 用LHS生成500个样本
- 各样本执行线性静分析 → 应力 → 用S-N曲线算疲劳寿命
- 得到寿命分布
结果是，平均寿命 $2.1 \times 10^6$ 循环，但5百分位寿命（B5寿命）仅 $4.5 \times 10^5$ 循环——换句话说，20个里面1个会在设计寿命的1/4处失效。只看平均值的确定性分析看不到这个风险。
哇，平均值看着安全其实有这么大的风险… 蒙特卡洛法看"波动的影响"真的很关键啊。
蒙特卡洛的软件比较
商用工具的MC功能
在CAE中用蒙特卡洛法，应该用什么软件？ 要自己写脚本吗？
主要的CAE平台都有内置的UQ（不确定性量化）模块。可以自己写Python脚本，也可以在GUI上设定。
| 工具 | UQ模块名 | MC/LHS対応 | 特记事项 |
|---|---|---|---|
| Ansys Workbench | Design Exploration (Six Sigma Analysis) | MC + LHS | 可用参数化设计语言(APDL)实现自动化 |
| Abaqus/Isight | Isight (Dassault SIMULIA) | MC + LHS + 最优化 | 数据流型的执行引擎 |
| Simcenter 3D | Monte Carlo Simulation | MC + LHS | 与Teamcenter连携进行数据管理 |
| COMSOL | Uncertainty Quantification Module | MC + LHS + PCE | 对多物理场耦合强 |
| OptiSLang (Dynardo) | 内置 (Ansys族) | MC + LHS + MOP | 感度分析最优化统合很充实 |
| LS-OPT | LS-Dyna用UQ模块 | MC + LHS + 元模型 | 碰撞成形分析的波动评估 |
开源选项
商用工具贵啊… 用开源可以吗？
实际上最灵活最强的往往是开源的组合：
- Dakota (Sandia National Labs) — UQ的世界标准。整合了MC、LHS、PCE、感度分析、最优化。与众多FEM/CFD求解器的接口很丰富
- OpenTURNS — 法国EDF/Airbus/Phimeca开发。Python库。分布拟合、相关模型、可靠性分析强
- UQLab (ETH Zurich) — MATLAB/Python。最新UQ方法（PCE、Kriging、子集模拟）高质量实现
- SALib (Python) — 专注感度分析。Sobol指数、Morris法的组合
Dakota名字经常听。能与OpenFOAM、CalculiX组合用吗？
可以的。Dakota有"黑箱"接口，用输入文件模板嵌入参数，外部调用求解器，从结果文件读取输出。OpenFOAM、CalculiX、Code_Aster，什么都行。
蒙特卡洛的前沿研究
代理模型的并用
一次FEM要30分钟的话，1000次就500小时… 实务不可能啊。
这时候用代理模型（Surrogate Model）。用少量高精度模拟（数十到数百点）构建响应曲面，蒙特卡洛抽样就对这个轻的代理模型执行。
| 代理模型 | 学习需要的点数 | 特征 |
|---|---|---|
| 多项式混沌展开（PCE） | $O((d+p)! / (d! \cdot p!))$ | 能分析计算统计矩 |
| 高斯过程回归（Kriging） | $10d$〜$50d$ | 同时估计预测的不确定性 |
| 神经网络 | $100d$〜 | 高维非线性强，可解释性低 |
| RBF（径向基函数） | $10d$〜$30d$ | 实现简单，高维精度下降 |
"PCE能分析计算统计矩"是什么意思？
PCE是用直交多项式展开响应。比如输入是正态分布就用艾尔米特多项式：
因为直交性，平均和方差可以从系数直接算：$\mu_Y = c_{\mathbf{0}}$，$\sigma_Y^2 = \sum_{\boldsymbol{\alpha} \neq \mathbf{0}} c_{\boldsymbol{\alpha}}^2$。根本不用蒙特卡洛抽样。但变量多的时候（$d > 10$）系数数增长爆炸，要用稀疏PCE或自适应方法。
大规模并行蒙特卡洛
不用代理模型，直接靠"蛮力"用HPC并行能行吗？
蒙特卡洛法是完全并行（Embarrassingly Parallel）的。各样本独立计算，1000个核就能同时跑1000个样本。这是其他UQ方法（比如自适应PCE）没有的优势。
用云HPC（AWS Batch、Azure CycleCloud、Google Cloud HPC）能在需要时借数千核来"突发计算"——一次30分钟的FEM×1000，100并行来5小时，然后释放。这就实际可行了。
蒙特卡洛的故障排除
常见失败和对策
蒙特卡洛模拟容易卡的地方有什么？
现场常见的问题总结一下：
| 现象 | 原因 | 对策 |
|---|---|---|
| 部分样本FEM不收敛 | 极端输入值的组合（分布尾部）在物理上不现实 | 对输入分布设截断（±3σ等）。提前决定失败样本的处理方针 |
| 结果分布明显不物理 | 忽视了输入变量间的相关性 | 设定相关矩阵。用Nataf或Rosenblatt变换反映相关 |
| 累积平均一直不收敛 | 输出分布有重尾（对数正态、帕累托等） | 应用方差低减法。或对输出的对数值进行统计处理 |
| 用LHS也不见方差低减 | 响应函数非单调非线性，分层效果弱 | 采用Iman-Conover法的相关性控制LHS |
| 结果不重现（每次都不同） | 伪随机数种子没固定 | 固定随机数种子确保重现性。报告中记录种子值 |
"部分样本FEM不收敛"很困扰。1000次中5次失败了，那5个怎么处理？
有3个方案：
- 除外并用 $N_{\text{valid}}$ 统计 — 最简单，但失败不是随机而是系统性的话（总是高应力侧失败）就有偏差
- 给失败样本"最坏值" — 保守推定，安全侧判断有理
- 调查失败原因，改善网格和接触设定，重跑 — 最正确但费时
实务上失败率在1%以下用方案1没问题。超过5%就应该反思模型本身的稳健性。
随机数种子的事，基础但超级关键啊。报告里写"种子:42"就能重现。
就是这样。蒙特卡洛法是"统计方法"，重现性的确保是科学信信度的基础。种子值、抽样方法、样本数、输入分布参数——这些都要记录。V&V（验证和验证）的第一步就是这个。
了
更详细
错误