1. 确定性潮流为什么“不够用”概率潮流的工程需求1.1 光伏不确定性给潮流计算带来的实际冲击以前做潮流计算输入一张负荷表、确定发电机出力跑一次牛顿法拿到一组节点电压和支路功率就可以交差。这套流程在光伏大规模接入之后变得非常别扭光伏出力不是固定值晴天、阴天、云层飘过、早晚光照变化都会让同一节点的注入功率在几分钟内发生大幅波动。你在下午三点用预测均值跑出来的潮流在下午三点半可能完全不成立因为几十兆瓦出力消失电压分布随之整体变化。这种波动不是偶然事件而是常态。以常见的标准节点测试系统为例在某个负荷节点附近接入一个中等容量的光伏电站若当天气候条件一般出力可能在额定值的20%到90%之间来回摆动。确定性潮流只能回答“在这个给定的注入下系统状态是多少”而工程上真正需要回答的问题是“在这种随机注入下节点电压有多大的概率越限支路潮流有多大的概率过载”。这就是概率潮流计算的出发点。1.2 概率潮流到底要算哪些指标概率潮流和确定性潮流的核心区别是把输入的随机性转化为输出的概率分布而不是只给一个点值。具体要输出的指标大致有下面几类输出对象典型统计量工程用途节点电压幅值期望、标准差、PDF、CDF判断电压越限概率、确定无功补偿需求节点电压相角期望、方差、分布分析系统稳定裕度、筛查薄弱节点支路有功/无功潮流期望、越限概率支路选型、断面输电能力评估网络损耗期望、分位数经济性评估、网架改造方案比选举个例子某节点在确定性潮流里电压是1.04pu看着很安全但概率潮流算出来该节点电压高于1.05pu的概率是12%。这么一来调度员就不会只盯那个“平均值”而是会关注这个12%的发生条件比如发生在光伏大发时段那就要提前安排调峰或无功控制。所以概率潮流对规划的意义不是给出一个更精确的数字而是给出一个可以用于风险决策的概率指标。1.3 两条技术路线的选型思路现在工程和学术圈里实现概率潮流的方法很多但最常被拿来组合的就是蒙特卡洛模拟和半不变量法。蒙特卡洛模拟的思路最直接反复抽样、反复求解潮流最后统计结果分布。它的优点是几乎不依赖近似假设对非线性潮流方程、任意输入分布都能处理精度高但代价是计算量大几千上万次潮流求解在大型系统里非常耗时。半不变量法则走的是另一条路在基准运行点把潮流方程线性化然后用随机变量各阶半不变量在灵敏度矩阵上的叠加直接算出输出变量的半不变量最后用级数展开还原概率分布。它的计算量小到可以忽略但代价是它依赖于线性化和级数截断这两个近似环节。两条路线各有利弊实际工程中我习惯把蒙特卡洛当成“验收标准”把半不变量法当成“快速筛选工具”两种方法互相校验。下面的章节会把两条路线的原理和实现细节都拆开讲清楚。2. 光伏出力与负荷的不确定性建模Beta分布是常见起点2.1 从光照强度到光伏有功出力处理光伏不确定性第一步是给光伏出力选择一个合适的概率分布模型。光伏有功出力本质上来自光照强度光照强度受天气、季节、云层遮挡影响经验统计显示在一段固定时间尺度内归一化的光伏出力往往呈现中间高、两端低的分布特征Beta分布是业内最常用的近似模型。设某光伏电站的额定有功为P_rated归一化出力为p P_pv / P_rated则p在[0,1]区间内服从Beta分布其概率密度函数是f(p) Γ(ab) / (Γ(a)Γ(b)) · p^(a-1) · (1-p)^(b-1)其中a和b是形状参数。a和b的取值直接决定分布的形态a和b都偏大时分布向中间集中a较小而b较大时分布偏向低出力区间。实际建模时一般根据该电站历史出力的均值和方差用矩估计法计算这两个参数。设实测归一化出力的均值为μ、方差为σ²则参数估计公式为a μ · (μ(1-μ)/σ² - 1)b (1-μ) · (μ(1-μ)/σ² - 1)如果没有历史数据也可以直接取一组经验参数比如a5、b2表示光伏出力多数时间偏高的晴朗天气场景或者a2、b6表示多云天气场景。在程序中我用MATLAB的betarnd或者betainv函数来生成对应分布样本前者出的是随机数序列后者配合rand可以生成指定分位数。个人更推荐直接约定a、b参数因为后续需要解析计算半不变量时参数化的Beta分布比纯经验分布好处理得多。2.2 负荷波动的建模与参数选择除了光伏出力负荷本身也在实时变化。工程上做概率潮流时负荷通常被建模为以预测值为期望、以一定标准差波动的正态分布即P_load μ_load εε ~ N(0, σ_load²)需要说明的是这里直接用正态分布有一个隐患如果σ_load太大抽样时可能出现负荷为负的荒谬样本。所以我在程序里做了截断处理把低于μ_load - 3σ_load或高于μ_load 3σ_load的样本强行截掉不让极端值进入潮流计算。标准差一般取预测负荷的5%~20%具体取决于系统规模和分析场景规划分析可以取大一点运行校核则取小一点。负荷的有功和无功不是独立波动的。实际测量中无功的波动往往和有功保持一定的功率因数关联因此在建模时最好保持每个负荷节点的功率因数不变或者给无功也设一个与有功同源的随机乘数否则会出现“有功波动但无功纹丝不动”这种明显不符合物理的输入。这一点在蒙特卡洛抽样时不难实现但在半不变量法中容易被忽略后面会详细说。3. 蒙特卡洛路线为什么说它笨但可靠3.1 蒙特卡洛概率潮流的总体流程蒙特卡洛模拟的思路一句话就能讲清楚把随机因素抽样成千上万次每次得到一个确定的潮流解最后把这些解的期望、方差、分布统计出来。整个流程可以拆成下面五步设置抽样次数N并设定输入随机变量的概率模型当前抽样中生成光伏出力样本和负荷波动样本把它们叠加到节点注入功率上对修改后的注入功率执行一次确定性潮流计算记录本次计算的节点电压幅值、相角、支路潮流等输出量重复步骤2到4直到N次完成最后对记录结果做统计分析。这个流程几乎没有数学上的门槛唯一的要求是你已经有一个可靠的确定性潮流求解器。我在初始版本里直接用自己的牛顿法函数后来换到标准节点测试系统的数据格式时也只是把数据解析部分做了适配核心循环不变。3.2 MATLAB实现细节抽样、循环与随机数种子蒙特卡洛在MATLAB里的实现最大问题不是算法本身而是随机数种子管理。实际写代码时一定要在循环前用rng(固定种子)设定随机数生成器否则每次运行结果都不同后面验证半不变量法时连对照都没法做。一个典型的蒙特卡洛核心循环大概长这样rng(2024); % 固定随机种子确保结果可复现 N 5000; nPV 3; % 光伏电站数量 a 5; b 2; % Beta分布形状参数 V_store zeros(30, N); Pflow_store zeros(41, N); for k 1:N % 抽样光伏归一化出力并换算到系统基准功率 p_pv_sample betarnd(a, b, nPV, 1); P_pv P_pv_rated .* p_pv_sample / Sbase; % 抽样负荷波动 loadFactor 1 0.1 * randn(30, 1); % 10%标准差 P_load P_load_base .* loadFactor; % 更新节点有功/无功注入 Pbus P_gen_base - P_load P_pv_inj; Qbus Q_gen_base - Q_load; % 求解潮流 [V_k, theta_k] myNewtonPF(Pbus, Qbus, Ybus, bus_type); V_store(:, k) abs(V_k); Pflow_store(:, k) computeBranchFlow(V_k, theta_k, Yf, Yt); end % 统计结果 V_mean mean(V_store, 2); V_std std(V_store, 0, 2);如果安装了并行计算工具箱可以把循环改成parfor。但要注意parfor下每个worker的随机数序列与单核运行不一样必须在循环体内显式设置独立的随机数流否则并行前后结果无法对齐。更简单的做法是先在循环外生成全部所需的随机样本矩阵再用parfor只做潮流计算。这种“样本先备好、并行只计算”的写法可复现性最好。3.3 抽样次数怎么定抽样次数直接决定蒙特卡洛结果的稳定性和计算耗时。理论上N取得越大越好但实际工程中一定有个合理范围。以IEEE 30节点这样的规模我一般取5000次已经能稳定分辨电压期望的千分位变化如果只关心越限概率这种尾部指标可能要增加到10000次以上因为尾部事件本身出现概率低样本太少时统计波动很大。可以做一个简单验证把N从500逐步增加到10000观察某个关键节点电压期望值的变化。当N增加后均值不再明显漂移就说明样本量足够了。经验法则误差大致按1/√N衰减想提高一位有效精度样本量要放大100倍所以蒙特卡洛的精度靠堆样本成本确实不低。3.4 耗时优化的一个经验刚开始写蒙特卡洛循环时我每个样本都从统一的平坦初值启动牛顿法。后来发现把上一个样本的潮流解作为当前样本的初始值能明显减少牛顿迭代次数因为相邻两个样本的注入功率变化很小前一组的电压相角分布已经很接近新解了。算下来每个样本平均少迭代两到三次。在小系统里这个优化看起来无关痛痒但一旦系统规模上到几百节点、N上一万累计节省的时间相当可观。4. 半不变量法把概率计算降到线性代数级4.1 一次潮流 灵敏度矩阵 半不变量叠加半不变量法的核心思想是把潮流方程在基准运行点附近做一阶泰勒展开把非线性的随机传播问题转化为线性的灵敏度叠加问题。设系统基准潮流解为X_0和注入功率W_0当注入功率出现随机扰动ΔW时状态变量增量近似为ΔX ≈ J^(-1) · ΔW其中J是潮流计算的雅可比矩阵。令H J^(-1)则H就是灵敏度矩阵它的第i行第j列元素表示第j个注入随机变量发生单位变化时第i个状态变量的变化量。接下来的关键利用了半不变量的可加性。对于两个相互独立的随机变量U和V有κ_k(U V) κ_k(U) κ_k(V)如果随机变量X是多个独立随机变量ΔW_j的线性组合即X_i Σ_j H(i,j)·ΔW_j那么X_i的各阶半不变量可以逐项累加但要注意线性系数的幂次κ_k(X_i) Σ_j [H(i,j)]^k · κ_k(ΔW_j)第2节提到过光伏出力用Beta分布、负荷用正态分布建模各自的半不变量我们是可以求出来的。把它们按上述公式叠加到状态变量上就得到了节点电压、支路潮流等输出量的各阶半不变量。4.2 让输入随机变量的半不变量“对齐”半不变量法在程序实现中最容易出错的地方是输入随机变量向量与灵敏度矩阵列的对应关系。框架上需要做一个对齐动作把每个随机注入源按其所在的节点和类型有功还是无功排序然后构造对应的灵敏度矩阵列。以一个PQ节点接入光伏为例光伏只向该节点注入有功所以它的随机变量只对应注入功率向量中的P那一行而同一个节点上的负荷如果按固定功率因数建模负荷波动会同时引起有功和无功变化此时应该把负荷有功和负荷无功组合成一个随机变量而不是拆成两个独立变量。我在程序里用了这样一个简化方案每个负荷节点只定义一个随机乘子乘子同时作用于有功和无功两个分量这样既符合物理规律也让半不变量的叠加个数大幅减少。代数表达式上如果负荷节点j的有功扰动ΔP_j L_j·cosφ、无功扰动ΔQ_j L_j·sinφ其中L_j是负荷功率波动的随机变量那么状态变量增量的半不变量累加式应该写成κ_k(X_i) Σ_j [A_ij]^k · κ_k(L_j)其中A_ij cosφ·H(i,P_j) sinφ·H(i,Q_j)。这个细节不处理好无功和有功的权重就配错了结果会差得很离谱。4.3 Gram-Charlier级数从半不变量还原PDF/CDF求出输出状态变量的前几阶半不变量之后还要把它们还原为概率密度函数和累积分布函数。最常用的做法是Gram-Charlier级数展开以正态分布为基础用Hermite多项式去修正偏度、峰度带来的偏差。将半不变量标准化定义α_3 κ_3 / σ³α_4 κ_4 / σ⁴其中σ是标准差。随机变量标准化形式为z (X - μ) / σ则概率密度函数的四阶Gram-Charlier近似为f(z) ≈ φ(z) · [1 (α_3/6)·H_3(z) (α_4/24)·H_4(z)]其中φ(z)是标准正态密度函数H_3和H_4是Hermite多项式H_3(z) z³ - 3zH_4(z) z⁴ - 6z² 3累积分布函数则对密度积分得到或者直接用标准正态分布函数Φ(z)加上修正项近似。工程上一般取前四到六阶半不变量就够用了取到更高阶时级数修正项不仅不能保证收敛还可能让概率密度在局部出现负值这一点后续会详细讲。4.4 计算量对比为什么说差两个数量级蒙特卡洛和半不变量法的计算资源消耗根本不在一档。仍然以IEEE 30节点为例方法潮流求解次数附加矩阵运算等效耗时蒙特卡洛N50005000次牛顿法无分钟级半不变量法1次牛顿法1次雅可比求逆若干向量运算秒级以内半不变量法最划算的地方在于它把成千上万次非线性求解压缩成一次潮流求解和一次矩阵求逆。矩阵求逆的规模等于状态变量维度对30节点系统来说就是约50乘50的矩阵毫秒级就能完成。所以半不变量法特别适合做方案批量扫描比如不同光伏容量、不同接入点排列组合下的风险评估蒙特卡洛在这种场景下会算到怀疑人生。5. IEEE 30节点系统上的MATLAB实现细节5.1 测试系统与光伏接入方案本文程序以IEEE 30节点标准测试系统为对象。这个系统的数据很规范30个母线、41条支路、6台发电机基准功率一般取100MVA。母线1是平衡节点节点2、5、8、11、13是PV节点其余是PQ节点。程序中我直接载入了标准测试系统的bus和branch数据格式没有额外导入其他工具包。光伏接入方案参考了工程上的典型做法在三个不同的PQ节点各接一个光伏电站容量分别取15MW、20MW、25MW总容量60MW相当于系统总负荷的20%左右。这个渗透率下电压波动已经比较明显适合观察概率潮流的区分度。节点选择上也考虑了电气距离差异一个靠近发电机区域一个在系统中部一个在远端负荷密集区这样能看出不同位置光伏出力不确定性传播路径的差异。所有光伏电站统一用Beta分布建模取a5、b2意味着归一化出力期望约为0.71对应额定容量乘以期望后三个电站的期望出力合计约42.6MW。5.2 基准潮流和雅可比矩阵怎么处理半不变量法必须用到雅可比矩阵但直接调用现成潮流工具包时通常拿不到矩阵本身。我在实现时没有绕弯子直接把一个标准牛顿法潮流函数改写成了返回雅可比矩阵的版本。函数签名大致是function [V, theta, J] myNewtonPF(Pbus, Qbus, Ybus, pv, pq)函数内部先对PQ节点和PV节点编号构造未压缩的功率不平衡方程形成完整雅可比。返回的J矩阵维度等于未知量个数即PQ节点数×2 PV节点数。需要特别说明的是这里返回的J是最终迭代收敛时的雅可比矩阵不是初始点矩阵。因为半不变量法的线性化点就是基准潮流运行点只有收敛点处的灵敏度才反映真实运行状态的局部动态。5.3 灵敏度矩阵构造的关键细节得到雅可比矩阵J后直接求逆得到灵敏度矩阵H。但在构造H与注入随机变量的对应关系之前必须完成两类节点剔除平衡节点不参与灵敏度计算因为它的电压和相角是给定的不随注入波动改变PV节点的无功注入能力视为充足光伏随机注入只加到PQ节点的有功上不受PV节点电压控制逻辑干扰。举例来说30节点系统中PQ节点数目约为20个除去平衡节点的相角方程后状态变量总数约49个随机注入源包括3个光伏随机变量和多个负荷随机变量。实际程序里我不会把随机注入源的维度固定死而是设计成动态拼接这样换系统时不用大改代码。状态变量半不变量的累加是三层循环完成的kmax 6; kappaX zeros(nState, kmax); for i 1:nState for j 1:nInj for k 1:kmax kappaX(i,k) kappaX(i,k) H(i,j)^k * kappaInj(j,k); end end end外面两重循环的规模是nState乘nInj对30节点系统来说这个矩阵运算在一次矩阵乘法里就能完成后面的kmax循环只是逐阶累加。真正耗时的反而是求输入随机变量的半不变量因为Beta分布的解析矩计算涉及到Gamma函数数值上比较繁琐。5.4 输入分布半不变量的计算程序如果不想费劲推导Beta分布的解析矩公式一个更通用的办法是用“从样本估计半不变量”的数值路线。先生成大量分布样本然后由样本的中心矩换算半不变量。function cum cumulantFromSample(x, order) x x(:); mx mean(x); m (r) mean((x - mx).^r); cum zeros(order, 1); cum(1) mx; cum(2) m(2); % 方差 cum(3) m(3); % 偏度相关 if order 4 cum(4) m(4) - 3 * m(2)^2; end if order 5 cum(5) m(5) - 10 * m(3) * m(2); end if order 6 cum(6) m(6) - 15 * m(4) * m(2) - 10 * m(3)^2 30 * m(2)^3; end end这个数值方法有个低调的好处它不仅适用于Beta分布也适用于任何你能生成样本的分布比如带截断的正态分布。缺省情况下程序只计算前六阶因为Gram-Charlier展开用到六阶已经足够继续往高阶走不仅矩阵运算量增大数值稳定性反而会变差。5.5 结果输出画图比看数字更有效跑完两种方法程序最后会输出两类图表节点电压幅值的概率密度曲线和累积分布曲线以及指定支路的有功潮流概率分布。我个人的实操建议是CDF图优先于PDF图。工程决策人员看到PDF曲线往往一脸茫然但看到CDF曲线可以直接读出“电压低于0.95pu的概率是多少”“支路潮流超过限额的概率是多少”这类明确结论。程序里我专门封装了一个输出函数给定半不变量后直接返回到分位数表function q quantileByCumulants(cum, p) % 用Cornish-Fisher展开求百分位数 mu cum(1); sigma sqrt(cum(2)); z norminv(p); q mu sigma * (z (cum(3)/sigma^3)/6 * (z^2 - 1) ... (cum(4)/sigma^4)/24 * (z^3 - 3*z) ... - (cum(3)/sigma^3)^2/36 * (2*z^3 - 5*z)); endCornish-Fisher展开比Gram-Charlier更适合算分位数尤其在求95%这种尾部分位数时不会出现概率密度为负的尴尬问题。6. 两种方法实测对比误差从哪里来又该怎么排查6.1 低渗透率场景结果几乎一致当光伏渗透率较低比如总容量只有系统负荷的5%~10%时半不变量法和蒙特卡洛的结果非常接近。我对节点25的电压做了对照蒙特卡洛5000次的期望电压是1.0316pu半不变量法算出来是1.0314pu误差不到0.02%标准差上蒙特卡洛是0.0023pu半不变量法是0.0023pu几乎没有差别。越限概率方面因为电压偏离限值较远两种方法都给出接近0的结果。这说明在小波动场景下一阶线性化近似造成的误差完全可以忽略半不变量法完全能替代蒙特卡洛。这也是为什么很多规划工具愿意用它几十个接入方案排下来结果稳定且速度极快。6.2 高光伏渗透率场景线性化假设开始吃力当光伏总容量提高到系统负荷的30%以上时差距开始显现。蒙特卡洛的电压分布往往带有明显右偏特征而半不变量法在峰度上经常估不准导致CDF曲线在尾部出现偏差。有一次我在节点21附近接入40MW光伏并取a2、b6这种偏低压出力的分布蒙特卡洛算出该节点电压高于1.05pu的概率约3.5%半不变量法却只算出2.1%。这个误差在工程决策上已经不可忽视了。误差来源可以拆成三个部分一阶线性化误差潮流方程本身非线性注入波动大时二阶及更高阶项不能被忽略级数截断误差Gram-Charlier展开只取前几阶Hermite多项式分布偏斜越明显误差越大输入分布假设误差Beta分布和正态分布只是对真实随机过程的近似如果真实分布和假设不符任何方法都救不回来。6.3 我踩过的三个坑与修复思路第一个坑是Gram-Charlier展开算出概率密度为负。出现这种问题的原因通常是取到了六阶以上展开且分布的偏度较大。我的修复方法是把展开截断到四阶同时对偏度超过一定阈值的节点改用Cornish-Fisher求分位数不再用Gram-Charlier生成PDF曲线。第二个坑是灵敏度矩阵与随机注入变量错位。这个问题在程序里藏得很深明明单个输入变量的半不变量算对了但叠加后结果量级全乱了。排查链路是先取一个只有单一随机源的简化算例手动验证某个状态变量的半不变量是否满足κ_k(X)H^k·κ_k(W)确认没有错位后再逐步加回多个随机源。这个排查过程很枯燥但只要怀疑灵敏度矩阵就要坚决回头核对节点编号顺序而不是去调后面的级数展开。第三个坑最隐蔽蒙特卡洛和半不变量法用的输入分布不一致。比如蒙特卡洛用Betarnd(5,2)生成样本但半不变量法计算Beta分布矩时手写的参数却是(2,5)两种方法结果当然对不上。后来我在程序里加了一处断言先分别从两种方法得到的输入半不变量做一致性检验不一致则直接报错。这算一个很值得分享的小技巧能省掉大量排查时间。6.4 工程配合先快筛后验算基于实测对比我给这套程序的定位是“半不变量法批量筛选、蒙特卡洛重点验算”。具体操作是先用半不变量法扫描几十种光伏接入方案把电压越限概率超过1%的方案挑出来再对这些风险方案跑5000次蒙特卡洛做最终确认。这样既避免了蒙特卡洛在全方案扫描时的巨大计算量也不会因为半不变量法的近似误差导致风险漏判。在IEEE 30节点系统上这套配合流程单台普通电脑就能跑完。如果用全蒙特卡洛扫描50个方案每个方案5000次潮流总耗时可能是数小时而半不变量法扫描50个方案只需要几分钟再对其中五六个高风险方案做蒙特卡洛总耗时控制在半小时以内结论的可靠性也得到保证了。最后再说一点实际体会做概率潮流这几年我越来越觉得程序本身不是最难的部分难的是让分析结果真正被用起来。跑完两种方法之后习惯别急着去看那一大堆期望和标准差先画一条关键节点的电压越限概率随光伏渗透率变化的曲线再画一条CDF曲线给非专业的人看这两种图比任何均值方差表都有说服力。更实用的建议是写程序时一定要把输入随机变量建模和潮流求解解耦。我最初把Beta采样直接写进潮流循环里导致后来想换负荷模型时几乎重写了一遍。现在代码里三段式结构很清楚随机变量建模模块、确定性潮流求解模块、概率统计模块哪个模块出了问题单独排查互不牵连。以后再遇到馈线、微网或者含储能的系统也能快速套用这套框架。