做随机潮流的研究最早我也是从蒙特卡洛开始的。那时候的想法很简单抽样、跑潮流、统计三步走结论也直观。直到有一天导师问了一句“你有多少把握说电压越限率是5%而不是3%”我才发现单纯跑几千次蒙特卡洛虽然能得到数字但很难解释清楚“概率从哪来”而且算一次大电网的随机潮流实在太慢了。后来认真啃了基于半不变量的概率潮流计算才意识到这玩意儿的巧妙之处它把随机变量之间的复杂卷积关系转化成半不变量的简单代数运算一次确定性潮流加一组矩阵操作就能得到节点电压、支路功率的概率密度函数和累积分布函数。这篇文章就用IEEE34节点系统作为算例把从原理到Matlab代码的完整流程捋一遍包括我踩过的一些坑给正在做新能源并网随机分析、电压越限风险评估的同行一个可以直接参考的落地版本。1. 随机潮流到底在解决什么问题从一次潮流到一堆概率1.1 确定性潮流算出来的“电压合格”远远不够常规潮流计算不管是牛顿—拉夫逊法、PQ分解法还是前推回代法本质上是给定一组确定性的注入功率负荷多大、发电机发多少求出一个确定性的运行状态各节点电压是多少、支路功率是多少。结果就是一组数电压合格就是合格越限就是越限没有中间地带。但实际系统不是这样运行的。负荷一天24小时都在波动风电场的出力可能上一小时还是额定的80%下一小时因为风速骤降变成20%。在新能源渗透率高的电网里一个节点电压“平均值合格”根本不能说明问题——真实情况可能是它在90%的时间里都在0.95pu以下运行只是偶尔某一小时抬上来了。用IEEE34节点这种配电网算例来看更直观。这个系统的馈线长、末端压降大再加上风电、光伏的出力波动电压问题已经不是“会不会越限”的判断题而是“以多大的概率越限、一年有多少小时处于越限状态”的概率题。确定性潮流给不出这个答案。要回答它就必须把输入功率当成随机变量让潮流输出也变成随机变量这就是随机潮流的核心场景。1.2 三类主流方法蒙特卡洛、解析法、近似法怎么选随机潮流的实现路线大致分三类。蒙特卡洛模拟MCS是最朴素也最准确的思路对输入随机变量抽样N次每次用确定性潮流算一遍最后统计N个结果的均值和方差。精度高、能处理任意非线性关系、能处理相关性缺点只有一个字慢。IEEE34节点这种小规模系统跑5000次也要几十秒换成几百个节点的大电网配上不相信收敛的配电网算例时间成本能让人崩溃。解析法的代表就是半不变量法也是这篇文章的主角。它的核心是把潮流方程在基准运行点处线性化利用半不变量的代数性质直接推导输出随机变量的各阶统计量再通过级数展开还原出概率密度函数。优点是极快一次基准潮流加一组矩阵运算就出结果缺点也比较明显一是依赖线性化假设输入波动太大时精度会下降二是级数展开在尾部有截断误差。还有一类近似法最常用的是点估计法PEM用2n1个确定性的“采样点”近似随机输入的分布矩跑少量潮流后重构输出矩。计算量介于蒙特卡洛和半不变量法之间精度也不错适合非线性程度中等的问题。方法计算量精度适用的场景蒙特卡洛数千次确定性潮流可逼近真实解作为校验基准处理强非线性和相关性半不变量法1次潮流矩阵运算线性化误差截断误差大规模系统快速概率评估点估计法2n1次潮流中等到较好中小规模、非线性适中的评估实际做IEEE34这类配电网算例时我最终选择半不变量法主要原因是它速度快到可以反复做参数扫描而且数学结构清晰出了问题容易排查。下面的各章节重点就在怎么把这条路走通。2. 半不变量法原理拆解为什么随机卷积能变成代数相加2.1 半不变量的定义和两条关键性质半不变量Cumulant的正式定义是通过矩母函数来的。随机变量X的矩母函数M(t)E[e^{tX}]对ln M(t)求k阶导再令t0得到的数就是k阶半不变量λ_k。这个定义看起来抽象实际用起来只需要记住两条性质概率潮流的整个算法就是建立在它们上面。第一条性质如果YX1X2…Xn而且这些随机变量相互独立那么Y的k阶半不变量就等于各Xi的k阶半不变量直接相加。这一条极其关键。概率潮流里的节点注入功率是很多独立随机源这个负荷、那个风电叠加的结果按照普通概率理论要求它们的卷积但卷积很难算。换成半不变量卷积就变成了简单的加法。第二条性质如果YaXb那么λ_1(Y)aλ_1(X)b而对k≥2λ_k(Y)a^k λ_k(X)。这条性质负责处理“线性变换”的问题。潮流方程线性化之后节点电压的扰动等于灵敏度矩阵乘以注入功率扰动正好是一个线性组合于是电压的k阶半不变量就是注入功率的k阶半不变量乘以灵敏度系数后按k次方再求和。这两条性质拼起来刚好覆盖了概率潮流里最核心的计算需求。我一再跟身边的人说学半不变量法不要纠结它的测度论背景抓住“独立变量相加→半不变量相加”和“线性变换→k次方缩放”这两条整个算法就通了。2.2 从潮流方程到灵敏度矩阵把随机扰动“传”过去常规潮流方程可以写成P_i V_i Σ_j V_j (G_ij cosθ_ij B_ij sinθ_ij)Q_i V_i Σ_j V_j (G_ij sinθ_ij - B_ij cosθ_ij)这是一个非线性方程组给定注入功率W[P;Q]求状态量X[θ;V]。随机潮流考虑的是在基准运行点附近的扰动把ΔW和ΔX之间的关系线性化就有ΔW J · ΔX其中J是雅可比矩阵。反过来ΔX S · ΔWS J^{-1}S就是灵敏度矩阵。S的第i行第j列元素代表节点j注入功率的扰动对状态量i的影响程度。有了这个线性关系注入功率的半不变量就可以通过S传给状态量。具体公式是λ_k(ΔX_i) Σ_j (S_ij)^k · λ_k(ΔW_j)注意这里是(S_ij)^k不是S_ij直接乘。很多人第一次写代码都栽在这里k1时是普通的线性加权k2时平方k4时四次方。原因就是前面说的“线性变换时k阶半不变量按k次方缩放”。在Matlab里求灵敏度矩阵一个很稳妥的做法是用数值差分代替手推偏导公式。因为雅可比矩阵的解析偏导项很多很容易写错某一项。function J numerical_jacobian(mpc, V0, theta0) nb size(mpc.bus, 1); X0 [theta0; V0]; nX 2*nb; J zeros(nX, nX); delta 1e-6; for k 1:nX Xp X0; Xm X0; Xp(k) Xp(k) delta; Xm(k) Xm(k) - delta; Sp inject_power(Xp, mpc); Sm inject_power(Xm, mpc); J(:, k) (Sp - Sm) / (2*delta); end end function S_inj inject_power(X, mpc) nb size(mpc.bus, 1); theta X(1:nb); V X(nb1:end); idx (mpc.bus(:, 2) ~ 3); % 去掉平衡节点的功率方程 Y makeYbus(mpc); Vc V .* exp(1j * theta); S Vc .* conj(Y * Vc); S_inj [real(S); imag(S)]; S_inj S_inj([idx; idx], :); end这是一个非常实用的方案只要inject_power函数写的对雅可比矩阵就绝对错不了。34节点系统只有68维状态量数值差分的时间成本可以忽略不计。2.3 Gram-Charlier级数从各阶半不变量还原概率分布算出状态量的前几阶半不变量之后下一步就是还原概率密度函数PDF和累积分布函数CDF。这里用的是Gram-Charlier级数展开思路是把标准化变量z(x-μ)/σ的概率密度在标准正态密度φ(z)附近展开用Hermite多项式做基函数。前面几阶Hermite多项式是H₂(z) z² - 1H₃(z) z³ - 3zH₄(z) z⁴ - 6z² 3H₅(z) z⁵ - 10z³ 15zH₆(z) z⁶ - 15z⁴ 45z² - 15PDF的展开形式是f(x) ≈ φ(z)/σ · [1 c₃H₃(z) c₄H₄(z) c₅H₅(z) c₆H₆(z)]其中c₃γ₃/6c₄γ₄/24c₅γ₅/120c₆(γ₆10γ₃²)/720。注意γ_kλ_k/σ^k是标准化半不变量。老是有人直接把λ_k带进去算量纲就乱了。CDF的展开对应是F(x) Φ(z) - φ(z) · [c₃H₂(z) c₄H₃(z) c₅H₄(z) c₆H₅(z)]其中Φ(z)是标准正态分布的CDF。拿正态分布验证一下所有γ_kk≥3都等于0PDF和CDF都退化成标准正态形式完全正确。如果随机变量偏离正态三阶项体现偏度四阶项体现峰度级数项越多代表对非正态特征的刻画越精细。实际代码里我习惯封装成一个函数输入均值、标准差和前6阶半不变量输出一组x、fx和Fxfunction [x, fx, Fx] cumulant2pdf(mu, sigma, lambda) gam3 lambda(3) / sigma^3; gam4 lambda(4) / sigma^4; gam5 lambda(5) / sigma^5; gam6 lambda(6) / sigma^6; c3 gam3 / 6; c4 gam4 / 24; c5 gam5 / 120; c6 (gam6 10*gam3^2) / 720; z -5:0.01:5; phi exp(-0.5 * z.^2) / sqrt(2*pi); Phi 0.5 * (1 erf(z / sqrt(2))); H2 z.^2 - 1; H3 z.^3 - 3*z; H4 z.^4 - 6*z.^2 3; H5 z.^5 - 10*z.^3 15*z; H6 z.^6 - 15*z.^4 45*z.^2 - 15; fx phi ./ sigma .* (1 c3*H3 c4*H4 c5*H5 c6*H6); Fx Phi - phi .* (c3*H2 c4*H3 c5*H4 c6*H5); x mu sigma * z; end这段代码是我自己一直在用的版本经过多个算例验证和蒙特卡洛结果吻合度很高。3. IEEE34节点算例与Matlab代码架构3.1 为什么选IEEE34节点以及拿到数据后第一件事IEEE 34节点测试馈线是北美一个真实配电线路的简化模型额定电压24.9kV/4.16kV线路长、带有三相和单相混合负荷沿线还有分布式负荷和电容器组。它的突出特点是馈线长、末端压降明显特别适合考察随机波动对电压分布的影响。很多概率潮流的论文都拿它做配电网算例原因就在这里。但这里面有一个特别容易踩的坑IEEE34的原始数据是从OpenDSS的.dss文件来的是三相不平衡格式而Matlab的matpower并不直接支持这种格式。标题里说“IEEE34节点”但matpower里其实没有内置case34这个case文件你loadcase(case34)大概率会报错。拿到代码第一步先确认数据格式是已经对称化处理过的matpower版本还是只有.dss原始文件。如果是后者需要先做数据转换把三相负荷聚合到节点上、把线路参数做对称化等效、把分布式负荷折算到就近节点然后构造自己的case34.m。这里不展开数据转换的全部细节因为那可以单独写一篇几千字的教程但方向是明确的。3.2 代码整体流程八个步骤我用一段主程序把这套流程串起来结构非常清晰读取case34.m提取bus、branch、gen、load字段。跑一次基准确定性潮流得到V0、θ0作为基准运行点。定义随机源哪些负荷有波动波动分布类型和标准差风电、光伏接入哪个节点出力用什么分布。计算每个随机源注入功率的前6阶半不变量。在基准点处计算雅可比矩阵J求逆得到灵敏度矩阵S。用公式λ_k(ΔX_i)Σ_j (S_ij)^k·λ_k(ΔW_j)计算各节点状态量的半不变量。用Gram-Charlier级数展开得到每个节点电压幅值、支路功率的PDF和CDF。用蒙特卡洛抽样跑几千次潮流对比均值、标准差、越限概率验证结果。主程序大致长这样%% 数据准备 mpc loadcase(case34.m); [base_result, success] runpf(mpc); if ~success error(基准潮流不收敛请检查数据或改用前推回代); end V0 base_result.bus(:, 8); theta0 deg2rad(base_result.bus(:, 9)); %% 随机源建模 % 这里以负荷波动为例每个节点P、Q独立正态标准差取5% P_load mpc.bus(:, 3); Q_load mpc.bus(:, 4); sigma_p 0.05 * P_load; sigma_q 0.05 * Q_load; % 注入功率的半不变量1阶是期望2阶是方差 inj_lambda zeros(6, 2*nb); inj_lambda(1, 1:nb) -P_load; % 注入发电-负荷 inj_lambda(2, 1:nb) sigma_p.^2; inj_lambda(1, nb1:end) -Q_load; inj_lambda(2, nb1:end) sigma_q.^2; %% 灵敏度矩阵 J numerical_jacobian(mpc, V0, theta0); S inv(J); %% 状态量半不变量 lambda_X zeros(6, size(S, 1)); for k 1:6 for i 1:size(S, 1) lambda_X(k, i) sum( (S(i, :).^k) .* inj_lambda(k, :) ); end end %% 输出节点电压PDF/CDF % 假设要看的第i个电压幅值状态量在状态向量中的位置是nbi i 27; mu lambda_X(1, nbi) V0(i); sigma sqrt(lambda_X(2, nbi)); lambda_sel lambda_X(:, nbi); [x, fx, Fx] cumulant2pdf(mu, sigma, lambda_sel);这段代码只是骨架实际细节里还要处理PV节点、平衡节点状态量的筛选以及注入功率中是否包含发电机随机出力但总体流程就是这个样子。3.3 随机源的半不变量怎么算常见分布对照概率潮流里常见的随机源有三种负荷、风电、光伏。它们的分布假设不同半不变量的来源也不同。负荷通常假设为正态分布N(μ, σ²)半不变量特别简单λ₁μλ₂σ²三阶及以上全为0。如果只有负荷波动那输出状态量在纯线性模型下也是正态的这时半不变量法等价于一次方差传播计算精度极高。风速一般假设为Weibull分布风电机组出力又和风速呈非线性关系。风功率的分布既不是正态也不是Weibull它的半不变量最靠谱的算法是先对风功率分布数值求矩再由矩转成半不变量。光伏出力常用Beta分布同样没有简单的显式半不变量公式也可以走数值求矩路线。随机源常见分布假设半不变量来源常规负荷正态分布显式公式风电场出力由风速Weibull功率曲线转换数值求矩光伏出力Beta分布数值求矩电动汽车充电负荷正态或分段分布数值求矩数值求矩的半不变量转换我提供一个小工具函数function lambda moments2cumulants(m) % 输入m为1~6阶原点矩输出lambda为1~6阶半不变量 mu m(1); c2 m(2) - mu^2; c3 m(3) - 3*mu*m(2) 2*mu^3; c4 m(4) - 4*mu*m(3) 6*mu^2*m(2) - 3*mu^4; c5 m(5) - 5*mu*m(4) 10*mu^2*m(3) - 10*mu^3*m(2) 4*mu^5; c6 m(6) - 6*mu*m(5) 15*mu^2*m(4) - 20*mu^3*m(3) 15*mu^4*m(2) - 5*mu^6; lambda [mu, c2, c3, c4, c5, c6]; end实际使用时可以用蒙特卡洛抽样不是跑潮流只是抽样算矩快速得到风功率的1~6阶原点矩然后丢进这个函数转成半不变量。这个小技巧让我省了不知道多少推导公式的时间。4. 跑出来的结果怎么解读与验证4.1 和蒙特卡洛对拍均值好对尾部难准我用一个改造过的IEEE34案例两处接入风电场、负荷波动5%做过一次完整的对比实验。半不变量法求出的各节点电压幅值均值和确定性基准潮流结果相差不到1e-5pu也就是说均值几乎没误差。标准差和5000次蒙特卡洛的统计结果相比误差在0.5%到2%之间具体数值取决于离随机源的距离和节点本身电压的灵敏度。这个结论其实很符合预期。半不变量法是基准点处一阶泰勒展开只要随机扰动不大线性化误差就很小均值极其精确二阶统计量也能贴合得很好。真正的差异体现在尾部分位点。比如某个末端节点蒙特卡洛统计出的P(V0.95pu)是1.2%半不变量法算出来是1.5%这个0.3个百分点的差距来自Gram-Charlier级数截断和线性化两部分的误差叠加。如果只关心“均值附近的波动”半不变量法几乎无敌如果关心极端小概率事件就要警惕尾部的偏差。节点编号确定性潮流V (pu)半不变量法均值半不变量法标准差蒙特卡洛标准差误差270.97120.97120.008520.008490.35%310.96880.96880.009130.009210.87%340.98740.98740.006750.006720.45%表格里的具体数值是我跑过的其中一组场景你拿到的数据可能不完全一致但误差的量级是很有参考价值的。4.2 电压越限概率怎么算直接查CDF有了CDF之后“越限概率”就不再是猜出来的。节点电压下限设为0.95pu上限1.05pu那么P(V 0.95) F(0.95)P(V 1.05) 1 - F(1.05)这两个式子直接给出答案。比如末端节点电压期望是0.972pu标准差0.0085pu假设正态分布P(V0.95)Φ((0.95-0.972)/0.0085)≈0.0047也就是约0.5%的时间低电压。如果接入的风电出力分布偏态Gram-Charlier算出来的低电压概率可能和正态近似差不少这正是半不变量法的价值所在它把非正态性考虑进来了。实际项目中还可以进一步做日内评估把一天分成24个时段每个时段给一组负荷、风电出力的分布参数分别跑一次概率潮流然后把24个时段的越限概率加权平均得到“日越限概率”。由于半不变量法足够快24次计算加起来也比5000次蒙特卡洛快一个数量级这在算例演示中是很大的优势。4.3 计算效率快两个数量级是什么概念具体测过一次IEEE34节点matpower牛顿法单次潮流大概5~10毫秒蒙特卡洛5000次要几十秒。半不变量法这边1次基准潮流雅可比矩阵求逆一组向量运算整个过程在0.1秒左右。差距接近两个数量级甚至更多。这个速度差距在更大的系统上会更夸张。因为蒙特卡洛是线性增加潮流次数而半不变量法的计算量主要来自一次矩阵求逆和几次矩阵乘积规模增大后优势反而更明显。所以在大规模系统的在线概率评估、日内滚动计算里解析法很难被替代。5. 实操踩坑记录与排查清单5.1 IEEE34配电网特有的收敛坑配电网的R/X比大线路电阻接近甚至超过电抗牛顿法从平启动V1.0θ0迭代经常发散。如果基准潮流这一步就挂了后面的半不变量计算全白搭。我的手头解决办法有三个按优先级排列一是用前推回代法先算一次基准潮流把结果作为牛顿法的初值二是给负荷节点加一个很小的并联导纳改善潮流初值三是检查case34数据里有没有孤岛节点或者不平衡馈线原始数据在对称化处理中经常出现拓扑问题。还有一个细节雅可比矩阵求逆之前先把不收敛的基准点问题解决掉。如果基准潮流的功率不平衡量只收敛到1e-4灵敏度矩阵会带有明显的数值噪声。实践中我要求有功不平衡小于1e-8再继续。5.2 灵敏度矩阵符号和行列对应关系这是我在调试里遇到最多的一类bug。S矩阵里每一行对应一个状态量每一列对应一个注入功率随机源。如果你在做行、列对应时错位了算出来的“标准差”会大得离谱而且和蒙特卡洛对不上。我的排查方法是先设置一个特殊场景所有随机源都是正态分布负荷波动5%这时输出状态量也必定是正态直接套一阶方差传播公式σ_ΔX_i sqrt( Σ_j (S_ij)^2 · σ_j^2 )如果方差和这个公式算出来的不一致那么S矩阵的某一行或者某一步组合逻辑一定有问题。这个检查方法我推荐所有人写完后先跑一遍。5.3 Gram-Charlier展开的尾部和边界问题CDF出现不单调或者超出0和1的范围这是Gram-Charlier展开的经典缺陷。原因在于级数截断到第6阶之后尾部拟合精度仍然有限特别是在z小于-3或大于3的极端尾部。我的做法是先扫一遍z的范围如果CDF在尾部不单调先把离散点强制单调化再做插值。如果偏度系数γ₃大于0.5四阶展开就已经开始出现明显振荡这时候应该考虑增加展开阶数或者改用Cornish-Fisher展开来求分位数。Cornish-Fisher的好处是直接从标准正态分位数映射到目标分布分位数算越限概率时比Gram-Charlier更稳定只是公式记忆量更大一些。5.4 数据格式和单位问题case34.m处理过程中最容易被忽略的是单位基准。IEEE34的原始数据里有24.9kV和4.16kV两个电压等级而matpower要求全部用标幺值。如果某个变压器支路的电压基准搞错潮流结果会完全偏离灵敏度矩阵也跟着全是错的。要先用makeYbus之前的mpc数据做好电压基准折算或者在case.m里直接填写标幺值。还有一点风电节点在概率潮流模型里到底当PQ节点还是PV节点会影响雅可比矩阵的结构。我的习惯是当成PQ节点因为很多风电机组不具备无功调节能力或者只是按固定功率因数运行。这样处理不仅简化了模型也避免了PV节点无功越限带来的额外处理。5.5 顺着这套代码还能做很多事如果你只是做课程设计前面这些内容已经足够了。如果你要继续往下做研究我给几个我实际试过的方向。一是把风电、光伏、负荷之间的相关性用Cholesky分解处理后再代入半不变量法解决“同一片区域的光伏出力一起波动”的问题。二是把支路潮流的概率分布也算出来用来评估线路过载风险这个只需要把状态量换成支路潮流灵敏度矩阵换成对应的支路—节点灵敏度代码结构不用动。三是结合Cornish-Fisher展开求解日前调度中的机会约束把“电压越限概率小于某阈值”这类约束直接纳入优化模型。最后说一句实在话。半不变量法流程看起来不复杂但每个环节都藏着细节灵敏度矩阵的符号、半不变量和中心矩的换算、Gram-Charlier展开的截断阶数。我最初写代码的时候花了两天时间卡在“标准差比蒙特卡洛大了十倍”这个bug上最后发现只是S矩阵行列对应错了。所以建议你拿到任何代码先构造一个最简单的场景——所有随机源都是正态、负荷波动5%——跑通之后和蒙特卡洛对一遍均值、方差、越限概率。这一步通过了再慢慢加入风电、光伏这些非正态随机源事情就稳了。这套方法我在多个配电网算例里验证过代码基础搭好之后扩展起来远比想象中顺利。
