概率潮流这个东西，说实话，我刚接触的时候心里是有点打鼓的。做电力系统潮流分析，常规做法是把风电出力和负荷当成确定值，牛拉法一跑，出来一组节点电压和支路功率，完事了。但现实中风电和负荷是天生不确定的，风大风小，负荷高峰低谷，这些波动会在电网里传播，最终体现在节点电压越限概率、支路过载概率这些指标上。所以就有了概率潮流——把输入的不确定性考虑进去，输出的不再是单一数值，而是一个分布。
这篇文章主要想聊聊我实际做过的这类Matlab程序：基于蒙特卡洛模拟和半不变量Gram-Charlier级数展开的概率潮流计算。内容包括风电和负荷的不确定性建模、两种计算方法的原理、完整的程序架构和关键代码、以及我踩过的坑和调通之后的经验。适合正在做相关课题的电力系统研究生、刚接触概率潮流的工程师，以及想找一个可靠参考实现的人。
1. 为什么非要算概率潮流：确定性潮流到底缺了什么
1.1 确定性方法解决不了“概率”问题
先摆一个最简单的场景。某个节点接了一座风电场，风速今天是8m/s，明天可能是12m/s，风电场出力自然不同。传统潮流里，你把风电出力固定成一个数，比如额定容量的60%，然后计算潮流。问题是，这个60%是拍脑袋定的，风速实际分布可能让出力在20%到95%之间波动，那节点电压会不会越限？支路功率会不会超载？确定性潮流给不了答案，因为你只算了一个点，而不是整个空间。
电网规划人员和调度人员需要的不是“某一种场景下电压是1.03pu”，而是“电压超过1.05pu的概率有多大”。概率潮流要解决的就是这个问题：给定风电场出力和负荷的概率分布，求输出电气量（节点电压幅值、相角、支路有功无功）的概率分布，或者至少求出一阶矩、二阶矩，以及越限概率。
这里有个很关键的认知：概率潮流不是把多个确定性潮流结果简单拼起来。它需要一套系统的方法处理输入随机变量的分布、变量之间的相关性、以及非线性潮流方程对这些随机变量的传播。
1.2 蒙特卡洛和半不变量法的路线选择
目前主流的概率潮流方法可以粗略分成三类：模拟法、解析法、近似法。蒙特卡洛属于模拟法，半不变量加Gram-Charlier级数展开属于解析法，还有比如点估计法算是近似法。我在这篇文章里重点讲蒙特卡洛和半不变量法，因为它们是课题里最常被拿来对比验证的一对搭档。
蒙特卡洛的思路非常简单粗暴：根据风电场风速分布和负荷分布，抽样几万次，每次抽样都跑一次确定性潮流，最后把几万个节点电压、支路功率结果做统计分析，得到它们的概率密度函数、累积分布函数和越限概率。它的优点是不需要对潮流方程做任何线性化假设，模型适应性强，各种非线性环节都能塞进去；缺点是计算量巨大，几万次潮流计算，哪怕是用牛拉法，数据量一大也相当耗时。
半不变量法的思路则完全相反。它基于随机变量的半不变量（累积量）以及潮流方程在一定精度下的线性化假设，把复杂的卷积运算变换成简单的代数加法。也就是说，不需要跑成千上万次潮流，只需要在基准运行点做一次潮流计算得到灵敏度矩阵，然后通过半不变量和Gram-Charlier级数快速重构输出量的概率分布。计算效率比蒙特卡洛高出好几个数量级，代价是精度受线性化假设影响。
我在实际做程序时，习惯把两种方法都写了：蒙特卡洛作为“标准答案”，半不变量法作为“快速算法”，拿前者验证后者的精度，这也正是很多论文的标准做法。下面两张图（程序里画的）通常放在一起对比：一条CDF曲线是蒙特卡洛的，另一条是半不变量法的，看它们重不重合。
这里有一张简单的对比表，方便理解两条技术路线的差异：
| 对比维度 | 蒙特卡洛模拟 | 半不变量法+Gram-Charlier |
|---|---|---|
| 基本原理 | 大量抽样+确定性潮流 | 线性化+半不变量代数运算 |
| 计算开销 | 极高（上万次潮流） | 低（一次潮流+多次简单运算） |
| 精度 | 高（理论上可收敛到真值） | 中高（受线性化假设限制） |
| 非线性适应能力 | 强 | 弱 |
| 适用场景 | 验证、精细分析 | 在线计算、快速评估 |
2. 不确定性建模：风电和负荷的“性格”要对得上
2.1 风速的Weibull分布与风功率曲线
做概率潮流的第一步，不是写潮流程序，而是先把不确定性源建模搞清楚。风电出力不确定性的源头是风速，工程上最常用的风速概率模型是两参数Weibull分布，概率密度函数长这样：
f(v) = (k / λ) * (v / λ)^(k-1) * exp(-(v/λ)^k)
其中v是风速，k是形状参数，λ是尺度参数。k一般在2到3之间，取2.5左右的时候曲线形态和实际风速统计比较接近；λ和平均风速相关。在Matlab里，直接用wblrnd函数就能生成服从Weibull分布的随机风速样本：
% 风速采样 k = 2.5; % 形状参数，根据当地风资源统计得到 lambda = 7.5; % 尺度参数，单位 m/s，和平均风速有关系 v = wblrnd(lambda, k, 1, N_samples); % N_samples 是抽样次数拿到风速之后，需要通过风功率曲线得到风电场有功出力。标准的风功率曲线是一个分段函数：低于切入风速不发电，在切入和额定风速之间近似按三次方关系增加，达到额定风速后输出额定功率，超过切出风速则停机保护。
% 风功率曲线换算，P_rated 为风电场额定功率 v_in = 3; v_rated = 12; v_out = 25; P = zeros(size(v)); idx = (v >= v_in) & (v <= v_rated); P(idx) = P_rated * (v(idx).^3 - v_in^3) / (v_rated^3 - v_in^3); idx = (v > v_rated) & (v <= v_out); P(idx) = P_rated;这一块看起来简单，但里面有个容易忽略的点：实际工程里风电场还有尾流效应、风机之间的相互遮挡、风电场的功率控制策略等，这些都会让实际出力和“单台风机功率曲线乘以台数”不一样。如果你做的是大电网级别的概率潮流，风电场往往整体建模成一个PQ节点或者PV节点，用简化功率曲线已经完全够用。如果研究的是风电场内部，那要上更精细的模型，那就是另一个课题了。
2.2 负荷的正态分布假设
负荷的不确定性，标准做法是假设各节点负荷服从正态分布，均值为预测值，标准差取均值的1%到10%之间，看负荷预测精度的实际情况。负荷预测越准，标准差比例越小。
对于正态分布采样，Matlab里用normrnd函数：
% 负荷采样，P_load0 为负荷母线有功预测值 sigma_ratio = 0.05; % 5% 的标准差比例 P_load = P_load0 + P_load0 * sigma_ratio * randn(1, N_samples);这里有个细节需要注意：无功负荷通常和有功负荷是耦合的，工程上一般假设负荷的功率因数保持不变，即无功和有功同比例波动。所以采样的时候，最好是先采有功，无功用功率因数乘出来，而不是独立采两个正态分布。我在做程序的时候踩过这个坑：刚开始独立采样有功和无功，结果导致某些节点的功率因数和物理实际差得很远，算出来的电压分布明显失真。
2.3 变量之间的相关性：不能假装它们互不认识
风电和风电之间、负荷和负荷之间，实际上是有相关性的。比如同一区域内的两座风电场，如果地理位置接近，风速会高度相关；同一个区域内的负荷，受气温和用电习惯影响，也会同步波动。但半不变量法中，随机变量的独立性假设很重要，否则半不变量的可加性就不成立了。
处理相关性的常用手段是通过Cholesky分解把相关的正态变量转化成不相关的标准正态变量。具体思路是：先构造各随机变量之间的相关系数矩阵R，对它做Cholesky分解，得到下三角矩阵L，使得R = L * L'。然后对独立标准正态采样矩阵Z做变换：X = L * Z，得到的X就具有指定的相关系数矩阵。之后再通过等概率变换（比如把正态分布的分位数映射到Weibull分布）得到相关风速样本。
需要说明的是，这种做法在处理正态变量时比较标准；处理非正态变量时，需要在相关正态空间和原始分布空间之间做等概率转换，过程稍复杂一些。如果只是做最简单的验证，一般可以忽略相关性，先假定所有输入变量独立。但在实际工程汇报中，评委或者导师很容易问一句“你们的输入变量有没有考虑相关性”，所以程序里最好预留相关性模块的接口。
3. 蒙特卡洛模拟：最“笨”但最可靠的办法
3.1 采样流程和样本量怎么定
蒙特卡洛的程序流程可以说没有任何智力上的难度，但工程细节反而最考验人。我第一次写的时候，直接把抽样、潮流计算、统计塞在一个大循环里，结果跑了一个多小时才结束，数据存得还乱七八糟的。后来重构了程序，把流程拆成了清晰的模块。
整体流程是这样的：
- 读入电网数据，包括母线参数、支路参数、发电机参数、负荷参数。
- 设置抽样次数N（比如10000次）。
- 对每一次采样：
- 根据风速Weibull分布采样，经风功率曲线计算风电场有功出力；
- 根据负荷正态分布采样，得到各节点有功和无功；
- 更新潮流计算中的PQ节点注入功率；
- 调用牛拉法潮流计算，记录该次采样的节点电压幅值、相角、支路功率；
- 循环结束后，对所有记录做统计分析：求均值、标准差，画概率密度估计、经验CDF，算越限概率。
样本量怎么定？理论上蒙特卡洛的收敛速度正比于1/sqrt(N)，也就是说想提高一位小数精度，样本量要增加到100倍。实际操作中，如果只是算均值，一般5000次就有比较好的效果了；如果要看分布尾部（比如越限概率），建议至少20000次以上。我在验证半不变量法的精度时，通常取20000到50000次蒙特卡洛作为基准。
这里很关键的一点是：不要在Matlab里用for循环写牛拉法求解器和数据记录全部挤在一起。实测下来，直接在循环内反复构造稀疏矩阵并调用求解函数，速度慢得让人想砸电脑。更好用的是预分配存储空间，把每次潮流计算的结果存入一个预先分配好的大矩阵，避免动态扩展内存。
3.2 牛拉法潮流计算的时间和空间优化
概率潮流里的潮流求解器，本质就是确定性潮流，你用Matpower的runpf也行，自己写一个牛拉法也可以。但如果你自己写，有几点可以优化：
- 用稀疏矩阵存储导纳矩阵Y，节点导纳矩阵在电力系统里极度稀疏，全矩阵存储会浪费大量内存并拖慢速度；
- 预先计算导纳矩阵并因子分解，因为每次采样只是改变注入功率，导纳矩阵不变，所以Y矩阵和LU分解可以在循环外算好；不过注意，牛拉法每次迭代的雅可比矩阵是会变的，这个没法完全复用。
我用的是自己实现的一段基于PQ分解和牛拉法混合的求解器，对于IEEE 30节点系统，单次潮流大约在毫秒级别，20000次采样跑下来几分钟能完成，可以接受。如果你用Matpower，代码会简洁很多，但每次runpf内部的求解开销较大，20000次可能需要近一个小时，需要提前有心理准备。
3.3 收敛判据和无效样本处理
蒙特卡洛采样过程中有一个很容易被忽略的问题：潮流计算不收敛。某些极端的风电出力和负荷组合，会导致潮流根本无法收敛，这在实际项目中并不罕见。处理方法不是简单的continue跳过，而是要做记录。如果你跳过的样本数比例很高，说明这个运行点在工程上本身就接近静态电压稳定极限，这是很有价值的信息，不是应该悄悄丢掉的数据。
我在程序里一般设定牛拉法收敛精度为1e-8的功率偏差，最大迭代次数30次。如果30次还没收敛，我会把这组输入记录下来，统计不收敛次数和相应的风电场出力、负荷水平。这样最后写报告的时候，可以说明“在多少比例的场景下潮流不收敛，对应的边界条件是什么”。
4. 半不变量法：把卷积变成加法的魔法
4.1 先理解半不变量是什么
半不变量，也叫累积量，统计里用希腊字母κ表示。它的好处在于有非常漂亮的可加性：独立随机变量之和的半不变量，等于各自半不变量之和。注意，这种性质一阶矩、二阶矩也都满足，但三阶中心矩不满足加法，而半不变量每一阶都满足。这正是概率潮流里使用半不变量的根本原因。
具体来说，如果随机变量Y是若干个独立输入随机变量的线性组合，比如：
Y = a1X1 + a2X2 + ... + an*Xn
那么Y的各阶半不变量可以直接由Xi的半不变量加权求和得到，不需要做卷积，不需要做复杂的积分。而潮流方程在运行点附近做泰勒展开，忽略二阶及以上项之后，节点电压和支路功率的随机偏差恰好可以表达成输入随机变量偏差的线性组合。两个条件一拍即合，这就是半不变量法能大幅简化概率潮流计算的底层逻辑。
4.2 由矩求半不变量怎么算
程序实现上，需要从输入随机变量的分布出发求各阶半不变量。最直接的方式是先求随机变量的各阶原点矩，再通过矩和半不变量之间的递推关系计算半不变量。
n阶原点矩mn的定义是E[X^n]，对离散样本就是均值，对连续分布就是概率密度积分的展形。半不变量和原点矩之间的递推关系，可以用下面的公式表达：
κ1 = m1 κ2 = m2 - m1^2 κ3 = m3 - 3m1m2 + 2m1^3 κ4 = m4 - 4m1m3 + 12m1^2m2 - 3m1^4 - 3m2^2 κ5 = m5 - 5m1m4 + 20m1^2m3 - 60m1^3m2 + 24m1^5 - 10m2m3 + 30m1m2^2
这些公式看起来吓人，但是代码里用一个已知的递推算法就能搞定，不需要手算到五阶以上。一般取到六阶半不变量就够用了，再高阶的稳定性反而变差。
对于常用的正态分布、Weibull分布，可以直接用概率论里的公式计算矩。Matlab里可以用数值积分求矩，比如用integral函数，或者对理论分布直接计算：
% 以Weibull分布前四阶原点矩为例 k = 2.5; lambda = 7.5; m1 = lambda * gamma(1 + 1/k); m2 = lambda^2 * gamma(1 + 2/k); m3 = lambda^3 * gamma(1 + 3/k); m4 = lambda^4 * gamma(1 + 4/k);这里gamma是伽马函数，Matlab里就是gamma函数。再代入上述递推式求出各阶半不变量。
4.3 Gram-Charlier级数：从半不变量到概率密度函数
有了输出变量的各阶半不变量之后，我们还需要还原成概率密度函数。这就是Gram-Charlier级数登场的地方。
核心思想是：任何一条分布曲线，都可以用标准正态分布的概率密度函数φ(x)及其各阶导数来逼近。经过推导，概率密度函数可以表示为：
f(x) = φ(z) * [1 + (κ3/6) * H3(z) + (κ4/24) * H4(z) + ...]
其中z是标准化之后的变量，H3、H4是Hermite多项式。实际写程序时需要的Hermite多项式：
H3(z) = z^3 - 3z H4(z) = z^4 - 6z^2 + 3 H5(z) = z^5 - 10z^3 + 15z
这里κ3和κ4分别是标准化后的三阶、四阶半不变量。注意，κ3和κ4其实对应着统计里的偏度和峰度，所以Gram-Charlier级数的本质，就是在正态分布基础上做偏度修正和峰度修正。
实际程序实现中，标准化方式是把随机变量Y减去均值再除以标准差：
z = (y - mean_y) / std_y;
然后把上面级数加进去就得到了概率密度。累积分布函数同样可以通过对标准正态累积分布积分得到。Matlab代码不复杂：
% 已知输出半不变量 kappa，求Gram-Charlier展开的PDF值 mu = kappa(1); sigma = sqrt(kappa(2)); z = (x - mu) / sigma; pdf_val = normpdf(z) .* (1 + ... (kappa(3)/sigma^3)/6 * (z.^3 - 3*z) + ... (kappa(4)/sigma^4)/24 * (z.^4 - 6*z.^2 + 3));就这一段代码，是整篇文章的核心之一。实测下来，对于节点电压幅值这种近似对称的分布，Gram-Charlier展开到四阶已经能给出相当不错的拟合效果；对于支路功率这种偏态比较明显的量，可能需要六阶展开才能把尾部分布形态描出来。
5. Matlab程序实现：从零搭一个能跑的概率潮流工具
5.1 程序模块怎么划分
一个完整的概率潮流Matlab程序，我建议按模块组织，文件结构大致如下：
prob_pf/ ├── data/ │ ├── bus_data.m % 母线数据 │ ├── branch_data.m % 支路数据 │ └── gen_data.m % 发电机数据 ├── models/ │ ├── wind_sampling.m % 风速采样与风电出力计算 │ ├── load_sampling.m % 负荷采样 │ └── cumulant_calc.m % 半不变量计算 ├── pf/ │ ├── newton_pf.m % 牛拉法潮流 │ └── jacobian.m % 雅可比矩阵形成 ├── analysis/ │ ├── gram_charlier.m % Gram-Charlier级数还原 │ └── monte_carlo.m % 蒙特卡洛主循环 └── main.m % 主控脚本这种模块化组织不是形式主义，调试的时候好处非常明显。我之前把采样和潮流写在一起，出问题之后根本分不清是采样出错了还是潮流求解器报错，只能到处断点调试。拆分模块之后，每个环节都可以单独验证。
5.2 关键代码实现思路
先说半不变量法的主流程。这一段的逻辑其实是比较直的：
- 在基准运行点（风电出力和负荷都取期望值）计算一次确定性潮流；
- 从潮流结果中提取雅可比矩阵或者灵敏度矩阵S，它描述了输入功率变化如何影响输出电压和支路功率；
- 对每个输入随机变量（风电场出力、负荷），计算其各阶半不变量；
- 利用灵敏度矩阵将输入半不变量线性组合成输出半不变量；
- 用Gram-Charlier级数重构输出概率密度函数。
以节点注入功率到节点电压的灵敏度为例，核心代码如下：
% 基准潮流得到 Jacobi 矩阵分解 [V0, S] = base_power_flow(bus, branch, gen); % S 是灵敏度矩阵（线性化系数） % 输入随机变量的半不变量矩阵，C_input 维度：k阶数 x 变量数 kappa_out = zeros(max_order, n_bus); for k = 1:max_order kappa_out(k, :) = abs(S) * C_input(k, :) .* sign(mean(C_input(k, :))); end这中间灵敏度矩阵的推导是程序里最麻烦的部分。具体来说，牛拉法最后一次迭代的雅可比矩阵（实际上是用最后一次修正量构造的矩阵），它的逆就是输入功率变化到节点电压变化的线性映射。这个细节很多资料没写清楚，我当时也是翻了很久文献才搞清楚，自己实现时务必注意。
蒙特卡洛那部分的代码更直接，核心循环如下：
N = 20000; Vmc = zeros(N, n_bus); Pflow = zeros(N, n_branch); for i = 1:N % 采样风速和输出 v = wblrnd(lambda, k); P_wind = wind_power_curve(v, P_rated); % 采样负荷 P_load = P_load0 + P_load0 * 0.05 * randn; % 构造该次采样的注入功率 bus(:, 3) = P_load; % 有功注入 bus(:, 4) = Q_load; % 无功注入 % 潮流计算 [V, success] = newton_pf(bus, branch, P_wind); % 记录 Vmc(i, :) = abs(V).'; Pflow(i, :) = branch_power_flow(V, branch); end5.3 输出结果与验证指标
程序跑完之后，我通常会输出三类结果：
第一类是概率密度对比图。选择某个风电场接入节点，画蒙特卡洛直方图拟合的概率密度曲线，和Gram-Charlier级数还原的概率密度曲线，放在同一张图里对比。视觉上就能直观看到匹配程度。
第二类是累计分布对比图。横坐标是节点电压，纵坐标是累积概率。这张图的价值在于可以直接读出电压越限概率，比如P(V < 0.95pu)是多少，P(V > 1.05pu)是多少。
第三类是数字指标。我习惯计算三个误差指标：均值绝对误差、标准差绝对误差、CDF上最大垂直距离（类似Kolmogorov-Smirnov统计量）。用这几项指标量化半不变量法的精度，写报告的时候有用。
下面是我跑IEEE 30节点系统，接入两座风电场后的典型结果（仅供示意参考）：
- 电压均值误差：小于0.0005 pu，基本可以忽略；
- 电压标准差误差：1%左右；
- CDF最大偏差：通常出现在分布的左尾，大约在0.02到0.05之间。
这个精度在工程上是可以接受的，但你必须清楚它是有代价的：半不变量法基于线性化假设，当风电场渗透率很高或者系统运行点接近电压稳定极限时，误差会明显放大。
6. 常见问题与排查技巧实录
6.1 Gram-Charlier展开出现负概率密度怎么办
这是用Gram-Charlier级数最大的坑，没有之一。理论上概率密度函数处处非负，但Gram-Charlier展开作为一种截断级数，在某些区域（通常是偏离均值较远的尾部）会出现负值。
我在测试支路功率的PDF时，就遇到过概率密度在尾部变成负数的情形，当时第一反应是程序写错了，反复检查索引和系数，后来查阅文献才明白这是级数展开本身固有的振荡问题。解决方法有几种：
- 将展开阶数从四阶降到三阶，有时候低阶反而比高阶稳定；
- 改用Cornish-Fisher展开，它直接作用于分位数函数，天然避开了概率密度为负的问题；
- 另一种常用手段是，不直接用Gram-Charlier还原概率密度，而是用Edgeworth级数，区别在于系数排列方式不同。
实际处理上，如果负概率出现在尾部远离关注区的地方，也可以选择忽略，只展示主要区间。但如果是出现在均值附近，那说明输出分布的偏度或峰度太大，线性化假设已经不太成立了，这时候更建议回到蒙特卡洛结果分析原因，而不是硬调级数。
6.2 蒙特卡洛到底要跑多少次才够
这个问题我每次都会跟同学讨论。原则是：要看你想分析什么指标。算均值和标准差，几千次就够；想算小概率的越限事件，比如0.5%的电压越限概率，你需要更多样本，不然统计波动太大。
简单估算一下：如果真实越限概率是p=0.005，那么N次采样中越限次数近似服从二项分布，标准差是sqrt(Np(1-p))。为了让相对误差控制在20%以内，需要sqrt(Np(1-p))/(Np) < 0.2，也就是N > 25(1-p)/p，大概需要5000次左右。但这只是估算。工程上我一般取20000次起步，敏感场景取50000次。
另外，有人会用拉丁超立方采样或者重要采样来减少方差，从而降低所需样本量。在概率潮流里，拉丁超立方是个不错的折中方案：用同样的样本量，估计精度能提升一个档次。如果你论文里对计算时间有要求，可以考虑在蒙特卡洛部分加一个拉丁超立方采样的选项。
6.3 线性化误差在哪些场景会“爆表”
半不变量法最怕非线性强的运行点。具体来说，以下几种情况会使半不变量法的结果和蒙特卡洛产生明显偏差：
第一，系统重负荷。当系统运行点接近潮流可行域边界时，电压对注入功率的响应呈现出很强的非线性，泰勒展开只取一阶项的假设失效，Gram-Charlier还原出的分布就会严重偏离实际。
第二，风电场渗透率特别高。风电出力波动范围大，比如从0到80%额定容量，这么大的输入偏差已经远远超出“小扰动”的范围，线性化自然会出现误差。
第三，负荷分布方差过大。如果某些节点负荷标准差设到均值的20%以上，输出的电压分布尾部会被拉得很长，级数展开很难完美拟合。
遇到这些情况时，我的实际建议是：拿蒙特卡洛结果做一次对比验证。如果偏差确实大，可以考虑半不变量法结合二阶敏感度分析，把潮流方程的泰勒展开保留到二阶项。只是那样推导和程序实现的工作量会成倍增加，适合作为进阶研究方向。
6.4 程序调试的顺序和技巧
写这套程序的时候，我的调试顺序非常重要，按照这个排查能少走弯路：
第一步，先单独验证不确定性建模模块。把Weibull采样结果画直方图，和理论概率密度曲线对照，确认分布参数设置正确；负荷采样同理。
第二步，验证半不变量计算模块。用标准正态分布做测试，因为它的各阶半不变量有解析结果：一阶是均值，二阶是方差，三阶以上全为0。如果程序输出的高阶半不变量不为0，那递推公式一定有问题。
第三步，验证半不变量法的线性组合逻辑。用一个小型测试系统，手动设置一两个随机变量，手动计算期望输出，再和程序结果对比。
第四步，最后才跑整个系统做蒙特卡洛对比。
这套顺序帮我避开了无数低级错误。很多人一上来就把所有模块连起来跑，结果出错后哪里都像有问题，排查成本极高。
最后再分享一个小技巧
做了这么多概率潮流程序，我最想提醒后来者的一点是：不要迷信“最先进”的方法，先吃透蒙特卡洛。蒙特卡洛虽然笨，但它给出的结果是你判断其他一切算法好坏的标尺。我见过有人直接跳到半不变量法，却连误差来源都说不清楚，就是因为没有蒙特卡洛结果做基准。
如果你刚开始做这个方向，建议先完整地跑通蒙特卡洛，把不确定性建模、潮流计算、统计分析整个链条走一遍；然后再加入半不变量法模块，逐步对比优化。这样虽然前期进度慢一点，但你对每个环节的理解都会扎实很多。程序框架搭好之后，后面无论是换更大的测试系统、加相关性模型，还是扩展成计及储能和柔性负荷的概率潮流，都只是在现有骨架上增加模块而已。
