简介这是一份以Matlab为工具、面向电力系统研究人员与工程师的蒙特卡洛法应用资源包聚焦状态估计与风险评估两大核心场景。包内共18个文件以17个.m脚本为主并附1个.asv备份文件涵盖电网建模、随机抽样、状态计算与风险分析等环节整体压缩包仅19KB结构紧凑便于直接阅读。已有307人学习下载。文件命名清晰包含负荷建模、故障概率、线性规划等模块可帮助读者快速理解如何在不确定性条件下完成系统状态的最优估计并模拟设备故障、负荷波动等情景以获得停电率、电压稳定性等风险指标。通过对照脚本逻辑与配套说明能够掌握利用随机数生成和统计工具箱搭建蒙特卡洛模拟框架的方法适合具备一定电力系统基础并希望提升状态估计与风险评估实战能力的读者。1. 为什么是蒙特卡洛电力系统状态估计里的概率答案做电力系统状态估计的人几乎都会遇到同一个尴尬量测值明明带噪声单次算出的状态结果却总被当成“唯一真相”。我见过不少同行拿着一次加权最小二乘结果就去评判风险结果某个极端工况早就把线路逼近越限单次估计根本抓不到。蒙特卡洛法在这里的价值不是炫技而是通过成百上千次抽样把电力系统不确定性的“全景”拉出来。这个资源把蒙特卡洛法和电力系统状态评估绑在一起用 MATLAB 跑一套完整流程输出状态估计误差分布和过载风险概率。适合做可靠性分析、调度校核、毕业设计的开发者也适合刚接触风险分析、想知道曲线怎么落地的新手。文章按“原理 → 代码 → 案例 → 排错 → 技巧”的顺序展开后面每一步都可以直接照抄。2. 状态估计到蒙特卡洛从方程到抽样的三条主线2.1 状态估计的核心方程从量测残差推回真实工况电力系统状态估计的目标是从一组冗余的遥测量测数据里推算出系统当前最可能的运行状态通常指母线电压幅值和相角。工程上最常用加权最小二乘WLS估计残差方程写成z h(x) e其中 z 是量测矢量h(x) 把状态映射到量测空间e 是服从高斯分布的量测噪声。目标函数是min J(x) [z - h(x)]ᵀ · W · [z - h(x)]W 是权重矩阵一般由量测协方差 R 求逆得到。对直流潮流模型h(x) 退化成线性关系状态 x 只包含非参考母线的电压相角。方程组简化成x_hat (Hᵀ W H)⁻¹ Hᵀ W z我一般在 MATLAB 里用反斜杠算子\(解这个方程避免显式求逆带来的精度损失。这个公式看着朴素却是后面蒙特卡洛循环的地基。2.2 风险评估为什么单次计算结果永远不够单次状态估计只暴露一种噪声实现下的工况。若噪声方差很小单次计算或许够用但电力系统中的量测噪声往往来自互感器、通信链路、滤波器参数漂移偏差范围并不小。更关键的是风险本身就是概率行为线路过载、母线电压越限、断面潮流超限只有在特定极端噪声组合下才出现。如果只跑一次等于只看了一个随机样本无法回答“母线相角偏差超过 0.1 弧度的概率有多少”这类问题。蒙特卡洛法把这个问题转换成抽样统计从噪声概率分布中反复抽取量测向量每次执行一次状态估计最后把所有结果做统计分析。这相当于在数字世界里把“运行方式”重复了成千上万遍。工程上通常关心三类指标状态估计误差的均值与标准差反映估计器整体精度。越限概率比如线路潮流超过额定值的次数占比。风险期望值若把越限后果用成本函数描述蒙特卡洛还能输出经济损失期望。这三条主线贯穿整个资源后面的 MATLAB 实现全部围绕它们展开。2.3 直流潮流模型与量测雅可比矩阵的选择状态评估场景里为了控制计算量几乎清一色采用直流潮流模型。量测雅可比矩阵 H 在直流模型下是常矩阵不需要迭代更新这给蒙特卡洛循环节省了非常多时间。H 的构造取决于量测类型母线注入量测某节点上的实时注入功率对相角的偏导直接取该节点对应的导纳矩阵行。线路潮流量测支路两端相角差除以支路电抗得到非零分量为 ±1/x。若同时混用两种量测H 矩阵会变成稀疏矩阵。我一般先把母线电纳矩阵 B 建出来再根据量测清单拼接 H这样后续调量测类型时只需要改动拼装逻辑不会动核心求解器。3. MATLAB 代码走一遍蒙特卡洛状态评估最小实现3.1 构造模拟系统先定义拓扑和真实工况第一步需要一套小型算例系统。为了不引入额外依赖我直接手工定义 4 母线系统结构很简单母线 1 作为参考平衡节点母线 2 接发电机母线 3 和 4 带负荷。% 定义四母线直流潮流系统 clear; clc; nb 4; % 母线总数 branch [1 2 0.10; 1 3 0.20; 2 3 0.25; 2 4 0.15; 3 4 0.20]; % 每行表示 [从母线, 到母线, 电抗(p.u.)] % 直流模型下电抗归一到1.0基准容量 % 建立节点电纳矩阵 Bbus参考母线也列入 nbranch size(branch, 1); Bbus zeros(nb); for m 1:nbranch i branch(m,1); j branch(m,2); x branch(m,3); Bbus(i,j) Bbus(i,j) 1/x; Bbus(j,i) Bbus(j,i) 1/x; Bbus(i,i) Bbus(i,i) - 1/x; Bbus(j,j) Bbus(j,j) - 1/x; end % 去掉参考母线得到可逆简化电纳矩阵 B_reduced Bbus(2:end, 2:end);逻辑说明这段代码先把每一条支路导纳加到节点电纳矩阵的对称位置对角线减去该支路导纳得到导纳矩阵。直流模型下节点注入有功功率与相角的关系为 P B · θ因此简化矩阵必须去掉参考母线对应行和列否则矩阵奇异无法求逆。参数上要注意每条支路电抗都以标幺值给出不同缩放会导致潮流结果偏差极大。3.2 单次状态估计WLS 解算器与量测噪声注入当真实工况定义好后生成带高斯噪声的量测数据。惯用做法是先用真实注入功率通过直流潮流方程解出真实相角再给量测加噪声。% 真实注入功率 (母线2发电机, 母线3-4负荷, 单位 p.u.) P_inj_true [0.6; -0.4; -0.5]; % 对应母线2、3、4 % 母线1作为参考其相角恒为0 % 直流潮流解出真实相角 theta_true B_reduced \ P_inj_true; theta_full_true [0; theta_true]; % 补上参考母线 % 设置量测噪声标准差 sigma 0.01; % 一次随机量测 z P_inj_true sigma * randn(3, 1); % WLS求解H矩阵就是B_reduced H B_reduced; R sigma^2 * eye(3); % 量测噪声协方差矩阵 W inv(R); % 权重矩阵 theta_hat (H * W * H) \ (H * W * z); theta_full_hat [0; theta_hat];逻辑说明这里假设所有母线注入量测量均存在因此 H 矩阵直接等于简化电纳矩阵。若量测类型是线路潮流H 会不再等于 B需要按支路拼接行向量。真实相角由无噪声的功率平衡方程求出但实际电网中真实相角是未知的这里作为“上帝视角”用于后续误差分析。WLS 求解的核心步骤是(H*W*H) \ (H*W*z)反斜杠自动选择合适的高斯消元法遇到稀疏矩阵也能处理。参数说明sigma 越大单次估计结果抖动越剧烈如果 sigma 设为 0.001单次结果几乎等于真实值蒙特卡洛统计出来的方差也会很小掩盖风险。工程上 sigma 应根据量测设备的精度等级给出比如 PMU 量测可取 0.005RTU 遥测常到 0.02 以上。3.3 蒙特卡洛循环批量抽样与风险指标统计单次估计只是热身。蒙特卡洛的真正价值在重复 N 次抽样后形成风险统计。下面的代码在循环里重复生成噪声量测调用 WLS 求解器再计算线路潮流越限概率。N 5000; % 蒙特卡洛抽样次数 store_theta zeros(N, 4); % 存放每次估计出的四个母线相角 line_flow zeros(N, nbranch); % 存放每次估计后的线路潮流 for k 1:N % 生成带噪声量测 z P_inj_true sigma * randn(3, 1); theta_hat (H * W * H) \ (H * W * z); theta_full [0; theta_hat]; store_theta(k, :) theta_full; % 计算每条支路潮流 for m 1:nbranch i branch(m, 1); j branch(m, 2); x_ij branch(m, 3); line_flow(k, m) (theta_full(i) - theta_full(j)) / x_ij; end end % 统计状态估计误差均值与标准差 theta_std std(store_theta(:, 2:4)); fprintf(相角估计标准差: %.4f, %.4f, %.4f\n, theta_std); % 定义线路潮流限制阈值 flow_limit 0.8; % 线路潮流标幺值上限 overflow_prob mean(max(line_flow, [], 2) flow_limit); fprintf(任何一条线路越限概率: %.4f%%\n, overflow_prob * 100);逻辑说明每次抽样后估计出的相角被补上参考母线的 0 值存入 store_theta接着循环逐条支路计算潮流。线路越限概率定义为至少一条支路潮流超过 0.8 p.u. 的比例用mean( max(...,[]) limit )一行实现max(...,[],2)求每行最大值。5000 次抽样下这个指标基本稳定。如果想分别统计单条线路越限率可以改成mean(line_flow(:,m) limit)。参数说明抽样次数 N 是蒙特卡洛最重要的旋钮。N 太小统计波动大风险概率越不稳定N 太大计算时间线性增长。一般先用 1000 次跑一遍看数量级再加大到 500010000 次做正式分析。在更复杂系统里我会用拉丁超立方抽样代替纯随机抽样让同样的 N 获得更低方差。4. 完整案例五母线系统的蒙特卡洛风险评估全过程4.1 案例设定模拟项目X的拓扑与参数把前面的代码扩展到 5 母线系统结构稍微做复杂一点两个发电机节点三个负荷节点七条支路量测随机删减部分线路故意留一些冗余度。这样既贴近真实调度系统也方便观察“量测冗余不够”时蒙特卡洛结果的变化。表模拟项目X的五母线系统参数母线类型注入功率p.u.电压等级说明1平衡节点由潮流决定参考相角 02PV节点发电机 0.8保持相角自由3PQ节点负荷 -0.5恒功率4PQ节点负荷 -0.6恒功率5PQ节点负荷 -0.3恒功率支路参数如下1-2 电抗 0.081-3 电抗 0.151-4 电抗 0.122-3 电抗 0.182-5 电抗 0.103-4 电抗 0.204-5 电抗 0.16。量测方案选择母线 2/3/4/5 的注入功率同时加上 1-2、2-3、4-5 三条线路的潮流量测量测噪声标准差取 0.01。噪声方差如果设置得过小比如 0.001风险概率曲线会过于集中如果过大比如 0.05几乎每次抽样都越限统计结果就失去区分度。实际项目中可以从 SCADA 历史数据里估出量测误差方差再代入模型。4.2 实现过程把量测雅可比矩阵拼出来在 5 母线案例中量测类型比第 3 章多需要手工拼接 H 矩阵。下面是可复用的拼装代码段。% 手动定义 H 矩阵的行排列第1~4行是母线2~5注入量测 H_list {bus_inj, 2; bus_inj, 3; bus_inj, 4; bus_inj, 5; line_flow, 1; line_flow, 3; line_flow, 7}; % 支路索引按 branch 表排列 % B_reduced 对母线2~5注入量测直接有效 H_bus B_reduced; % 线路潮流部分按支路参数计算偏导 H_line zeros(3, 4); H_line(1, :) line_h_matrix(1, branch, nb); % 支路1: 1-2 H_line(2, :) line_h_matrix(2, branch, nb); % 支路3: 2-3 H_line(3, :) line_h_matrix(3, branch, nb); % 支路7: 4-5 H [H_bus; H_line];function hrow line_h_matrix(line_idx, branch, nb) % 返回一条线路潮流量测对非参考母线相角的偏导行向量 % 非参考母线编号为 2~nb因此列索引 母线编号 - 1 i branch(line_idx, 1); j branch(line_idx, 2); x branch(line_idx, 3); hrow zeros(1, nb - 1); if i 1 hrow(i-1) 1 / x; end if j 1 hrow(j-1) -1 / x; end end逻辑说明母线注入量测的直接映射就是简化电纳矩阵的相应行线路潮流量测则是一条支路两端相角差的线性表示对应行向量的非零元素就是 ±1/x。上面的line_h_matrix函数把支路端点转换为非参考母线的索引返回稀疏稀疏行向量。这样无论量测数量怎么变H 矩阵都能动态拼装。4.3 结果怎么解读越限概率和风险曲线的工程价值蒙特卡洛循环运行完成后除了得到越限概率还可以画出一条经验累积分布函数。我通常会把“最大线路越限比例”随抽样次数变化的曲线打印出来观察是否收敛。% 计算每条线路越限率 flow_limit 0.75; overflow_rate_per_line mean(line_flow flow_limit, 1); % 画出累积均值检验概率收敛性 cum_risk cumsum(max(line_flow flow_limit, [], 2)) ./ (1:N); figure; plot(1:N, cum_risk); grid on; xlabel(抽样次数); ylabel(越限概率累积均值);逻辑说明cumsum计算每个抽样点之前越限事件累计次数逐步除以当前抽样数得到累计频率。这个曲线越接近水平线说明概率估计越稳定如果曲线还在大幅摆动就必须继续加大 N。工程上我习惯把收敛阈值设定为连续 500 次抽样内概率变化不超过 0.1%此时认为风险概率可信。对于五母线案例的典型结果是母线 4 附近线路过载概率最高因为负荷集中在下端而发电机出力主要来自节点 2输电走廊最容易先顶到限值。这个结果可以指导调度员提前安排机组出力组合或者把校核结果反馈到自动发电控制逻辑里。5. 避坑蒙特卡洛评估最容易翻车的四个开关5.1 病态矩阵雅可比矩阵奇异导致相角结果发散现象运行蒙特卡洛循环时某个抽样点的(H*W*H)接近奇异theta_hat输出数量级高达 1e4 的异常值越限概率瞬间被拉高到 80%。原因量测配置冗余度不足或系统中存在两个电气距离极近的母线导致信息矩阵条件数过大。直流模型中若某条线路电抗数值太小支路潮流量测对应的 H 行会出现接近线性的重复矩阵秩不足。蒙特卡洛抽样本身不能解决这个问题只会放大它。解决首先检查 H 矩阵条件数用cond(H*W*H)判断。若条件数超过 1e8就要增加量测或合并弱观测支路。另一个工程做法是给信息矩阵加一个很小的对角正则项如(H*W*H 1e-6*eye(4))代价是估计偏差略微增大但数值稳定性大幅提升。5.2 噪声协方差拍脑袋风险概率失真现象所有参数设置看起来正常但风险概率计算结果明显异常例如越限概率只有 0%或者 100%完全失去参考价值。原因sigma 被当作自由参数随手设置。若设置 0.001量测噪声几乎不存在蒙特卡洛抽样产生的波动远小于系统真实不确定性风险尾巴被完全截断设置 0.1则每次抽样都会越限概率统计失去区分度。解决我会从设备说明书或历史数据里找出量测误差限值比如某条母线注入量测误差的标准差估算为 0.015 p.u.就用这个值。如果实在没有历史数据可以用交叉验证法先跑一轮小抽样把估计结果残差的标准差作为新的噪声方差迭代两轮即可收敛。5.3 抽样次数太少置信区间宽得没法用现象程序运行 200 次抽样只花了几秒但每次重复运行越限概率在 1% 和 8% 之间剧烈跳动完全无法判断真实风险。原因蒙特卡洛估计量的方差按 1/N 递减200 次抽样对应标准差约为 3%对一个 5% 左右的概率事件来说置信区间太宽统计噪音盖过了信号。解决先用公式估算所需抽样次数。假设目标概率 p0.05希望置信区间半宽不超过 0.5%则 N ≈ z²·p(1-p)/ε²换算下来至少 7000 次。实际操作时先跑 1000 次观察累积均值曲线是否趋于平滑不平滑就继续增加。更高效的办法是改用拉丁超立方抽样让 5000 次抽样等价于上万次纯随机抽样。5.4 坏数据干扰一次异常量测污染整条风险曲线现象蒙特卡洛统计结果中某次抽样对应相角估计突然偏离真实值 10 倍导致误差均值和方差被严重抬高风险概率也异常偏大。原因量测噪声不严格服从高斯分布可能带有粗大误差比如通信丢包后自动补零、互感器饱和瞬间输出畸值。纯高斯噪声模拟不会产生这类跳变但真实系统的坏数据占比通常有 1%3%。解决在蒙特卡洛循环里加入坏数据检测逻辑。常见做法是计算估计残差若某条残差大于 3 倍标准差就丢弃该次抽样并跳过潮流统计。代码上只需要加一行判断residual z - H * theta_hat; if max(abs(residual)) 3 * sigma * sqrt(3) continue; end逻辑说明WLS 残差在理想高斯噪声下近似服从标准正态分布3 倍标准差之外的小概率事件视为坏数据。这行过滤代码能把污染样本挡在统计之外但要注意不要误删正常极端工况否则风险估计又会偏保守。实际操作中我倾向于对连续出现坏数据的母线发出告警而不是简单忽略。6. 进阶把抽样方差按下去的两种硬核手段蒙特卡洛跑到一定程度瓶颈不再是理论而是计算量和统计精度之间的权衡。这里分享两个我经常用到的技巧。第一个是拉丁超立方抽样Latin Hypercube Sampling。纯随机抽样经常出现样本扎堆的情况一千次抽样里某些区间被过度采样另一些区间却空洞无物。拉丁超立方把每个输入维度分成 N 个等概率区间每个区间强制抽取一个样本再随机排列组合。对于电力系统量测噪声这种低维分布效果非常明显。实现起来并不复杂先用rand生成区间内位置再叠加norminv转成高斯噪声。% 拉丁超立方抽样生成量测噪声 sigma 0.01; for k 1:N u (rand(3, 1) (0:2) (k-1)) / N; % 每个维度分N段 e sigma * norminv(u); z P_inj_true e; theta_hat (H * W * H) \ (H * W * z); % 继续把结果存起来 end这段代码把每个噪声维度的取值区间切成 N 份确保每个区间都有样本覆盖均匀程度远好于随机抽样同时保留了随机组合带来的相关性。第二个手段是重要性采样。如果关心的只是极端越限事件概率在 1% 以下直接蒙特卡洛需要跑到 10 万次才能看到几次事件。重要性采样的思路是从一个偏移后的分布中抽样让极端场景更容易出现最后用权重修正概率。对于状态评估这种高维问题偏移噪声均值即可% 重要性采样把噪声均值朝越限方向偏移 u sigma * randn(3, 1) 0.05; % 偏移量0.05人为增加越限概率 z P_inj_true u; % 计算似然比权重 w prod(normpdf(u, 0, sigma) ./ normpdf(u, 0.05, sigma)); % 越限事件乘以权重再累加这种方法得到的概率估计方差大幅降低但前提是偏移方向选得准。偏移太小改进有限偏移太大又会浪费大量样本在完全无关的区间。工程上我会先用 1000 次随机抽样找到最可能越限的支路方向再把偏移量设置成对应噪声分量的 23 倍标准差。最后说一个自己的习惯我有一条强制规矩每次跑蒙特卡洛之前先画出输入噪声样本的直方图确认覆盖范围符合预期。某次我在线路参数上加了个低方差噪声风险曲线的尾巴被压得平得诡异后来才发现噪声分布完全偏移到了不安全的区域。从那以后我每次都用histogram(randn...)和qqplot做一次输入分布检查再进入正式循环。手里有这份资源建议你也把校验步骤写进脚本里整个过程会少很多无头绪的返工。希望帮到你。本文还有配套的精品资源点击获取