电力系统随机潮流概率潮流计算MATLAB程序包含蒙特卡洛模拟法、半不变量法＋级数展开（Gram-Charlie，Cornish-Fisher）； 考虑光伏不确定性（Beta分布），以IEEE34节点为例，计算节点电压、支路潮流概率密度、累计概率并绘制曲线。 有注释，附带参考文献，不代做。 缺点是该节点系统不是很具有代表性，因为仅存在一个发电机节点，

最近在折腾电力系统概率潮流计算，发现这玩意儿和天气预报有点像——都是预测未来可能发生的状况。拿最常见的蒙特卡洛法来说，就像用无数个平行宇宙模拟电力系统运行状态。这次拿IEEE34节点开刀，先配置光伏出力Beta分布参数（这玩意儿形状像不像跷跷板？α=0.9,β=2就是个典型例子）。

% 光伏Beta分布参数设置
alpha = 0.9; 
beta = 2;
x = 0:0.01:1;
pdf = x.^(alpha-1).*(1-x).^(beta-1)/beta(alpha,beta);
plot(x,pdf,'LineWidth',1.5)

这曲线左偏的特征正好反映光伏中午出力高、早晚低的特性。接着上主菜——蒙特卡洛循环。每次循环就像给系统拍张快照，5000次循环下来硬盘可能会抗议，但MATLAB的矩阵运算能救命：

for iter = 1:MC_num
    % Beta分布随机采样
    PV_node = betarnd(alpha,beta,[nPV,1]); 
    % 节点注入功率修正
    busdata(PV_bus,PD) = base_load.*(1-PV_node);  
    % 调用牛顿法潮流计算
    [V,~] = newton_pf(busdata,linedata); 
    V_all(:,iter) = V; 
end

跑完这堆数据后，电压概率分布开始显形。某次实验中节点15的电压标准差达到0.012p.u.，这波动比老式电梯还刺激。画累计概率曲线时，发现95%概率下电压不会超过1.038p.u.——这种结论对规划电容器投切可太有用了。

半不变量法则像是个数学魔术师，先用泰勒展开把潮流方程掰碎：

% 半不变量法核心计算
J = jacobian_matrix(V,theta); % 雅可比矩阵
H = hessian_matrix(V,theta);  % 黑塞矩阵
gamma = cumulants(PV_dist);    % 输入变量半不变量
[V_cum,~] = series_expansion(gamma,J,H); 

Gram-Charlier展开那部分最有意思，把正态分布当成橡皮泥捏成想要的形状。Cornish-Fisher展开更绝，直接逆向工程分位数。不过要注意三阶以上展开可能会像拉过头的橡皮筋一样失真，特别是当光伏渗透率超过30%时。

虽然IEEE34节点用起来顺手，但它的单电源结构就像独轮车——平衡性差点意思。实际应用中建议搭配多电源系统测试，比如改造成含风电-燃气轮机的混合系统。毕竟概率潮流的终极目标，是要在电力系统的混沌中抓住确定性规律。

参考文献：

[1] Allan R N. Probabilistic methods in power system analysis[J]. 1994.

[2] H.P. Wong. Probabilistic load flow computation using the method of combined cumulants and Gram-Charlier expansion[J]. 2009.

[3] Borkowska B. Probabilistic load flow[J]. 1974.