4. 蒙特卡洛模拟原理与实现

在结构可靠性分析中，极限状态函数 g(X) 往往是非线性的，基本随机变量 X 的分布也不一定服从正态分布，此时一次二阶矩法（FORM）等解析方法会因线性化假设引入难以估计的误差。蒙特卡洛模拟（Monte Carlo Simulation, MCS）提供了一条几乎不依赖 g(X) 形式的通用路径：把失效概率写成数学期望，再用随机样本的平均值去估计它。

一、原理：把失效概率写成期望

设结构功能函数为 Z = g(X₁, X₂, …, Xₙ)，失效事件为 Z < 0。引入指示函数 I[g(X) < 0]，当样本落入失效域时取 1，否则取 0，则失效概率可写为：

P_f = ∫ I[g(x) < 0] · f_X(x) dx = E{ I[g(X) < 0] }

于是 P_f 变成了指示函数的数学期望。若从 f_X(x) 中独立抽取 N 组样本 x⁽¹⁾, …, x⁽ᴺ⁾，则估计量为：

P̂_f = (1/N) · Σ_{i=1}^{N} I[g(x⁽ⁱ⁾) < 0] = N_f / N

其中 N_f 为落入失效域的样本数。这个估计量是无偏的，其方差为 P_f(1−P_f)/N，变异系数为：

COV(P̂_f) = √[(1 − P_f) / (N · P_f)]

该式是蒙特卡洛模拟最重要的工程结论之一：当 P_f 很小时，所需样本量 N 近似与 1/P_f 成正比。若要求变异系数不超过 0.1，则 N ≈ 100(1−P_f)/P_f。例如 P_f = 10⁻³ 时约需 10⁵ 次模拟，P_f = 10⁻⁶ 时约需 10⁸ 次，这正是高可靠结构直接抽样成本高昂的根源。

二、大数定律与收敛性

大数定律保证：当 N → ∞ 时，样本均值依概率收敛于真实期望，即 P̂_f → P_f。中心极限定理进一步给出渐近分布：

(P̂_f − P_f) / √[P_f(1−P_f)/N]  →  N(0, 1)

由此可构造置信区间。取 95% 置信水平（z_{0.975} = 1.96），估计误差约为 ±1.96·√[P̂_f(1−P̂_f)/N]。收敛速度是 O(N^{−1/2})，即样本量增加 100 倍，误差才缩小 10 倍，这决定了单纯增加样本量并非高效策略。

实际编程中，应观察 P̂_f 随 N 的累积曲线是否趋于平稳，并同时监控变异系数。若 COV 仍大于 0.1，说明样本量不足或方差缩减未生效。

三、样本生成：从均匀分布到任意分布

计算机只能产生伪随机数 u ~ U(0,1)，其余分布由变换得到。最常用的是逆变换法：若 F 为连续分布函数，则 x = F⁻¹(u) 服从该分布。例如标准正态变量可用 Box–Muller 变换：

u1, u2 ~ U(0,1)
z1 = √(−2 ln u1) · cos(2π u2)
z2 = √(−2 ln u1) · sin(2π u2)

对数正态变量可由 x = exp(μ_ln + σ_ln · z) 得到，其中 μ_ln = ln m − 0.5σ_ln²，σ_ln² = ln(1 + v²)，m 与 v 分别为原变量的均值与变异系数。工程中还需注意：随机数发生器应设置固定种子以保证结果可复现；样本维度较高时应检查各变量间的相关性，必要时用 Nataf 或 Cholesky 分解引入相关结构。

四、方差缩减策略

由于 COV ∝ N^{−1/2}，直接增加样本量效率低。方差缩减的核心是改变抽样方式，使估计量方差下降而期望不变。

• 重要抽样（IS）：把抽样中心移到设计点附近，用重要性密度 h(x) 抽样，估计量改为 P̂_f = (1/N)Σ I[g<0]·f_X(x)/h(x)。对 P_f = 10⁻⁶ 的问题，IS 常能把所需样本量从 10⁸ 降到 10³～10⁴ 量级。
• 对偶变量法：用 u 与 1−u 成对抽样，若 g 单调，两者负相关，方差可显著降低。
• 分层抽样：把 [0,1] 分成若干层，每层按比例抽样，减少样本在空间上的聚集。
• 条件期望法：对部分变量解析积分，只对剩余变量抽样，直接削减一层随机性。

五、基础编程实践

下面给出一个可直接运行的 Python 示例：功能函数 g = R − S，R ~ N(200, 20)，S ~ N(150, 25)，用直接抽样估计失效概率。

import numpy as np

rng = np.random.default_rng(20240501)   # 固定种子，保证可复现
N = 200_000

R = rng.normal(200, 20, N)              # 抗力
S = rng.normal(150, 25, N)              # 荷载效应
g = R - S                               # 功能函数

Nf = np.count_nonzero(g < 0)            # 失效样本数
Pf = Nf / N
cov = np.sqrt((1 - Pf) / (N * Pf))      # 变异系数
beta = -np.sqrt(2) * 0  # 占位，见下

from scipy.stats import norm
beta = -norm.ppf(Pf)                    # 可靠指标

print(f"Nf={Nf}, Pf={Pf:.3e}, COV={cov:.3f}, beta={beta:.3f}")

示例运行结果（教学示例，非实测）：Nf≈3900, Pf≈1.95e-2, COV≈0.011, beta≈2.06。由于 P_f 在 10⁻² 量级，20 万次抽样已使 COV 远小于 0.1，结果可信。若把 S 的均值改为 180，P_f 会降到 10⁻⁴ 以下，此时同样样本量的 COV 将超过 0.1，必须改用重要抽样。

综上，蒙特卡洛模拟以指示函数期望为理论基础，以大数定律和中心极限定理为收敛依据，以样本生成和方差缩减为实现手段。掌握 P̂_f、COV、β 三个量的计算与判读，就具备了在结构可靠性分析中正确使用该方法的基本能力。