简介面向具备电力系统与概率论基础的研究人员、工程师及高校教师系统讲解广义多项式混沌法gPC在风光并网随机潮流中的应用用于解决新能源高比例接入带来的电压不确定性问题。资源围绕电压统计特征分析展开从正交多项式逼近原理出发依次介绍Hermite基函数生成、随机Galerkin投影、确定性方程求解、统计特征提取及蒙特卡洛验证同时深入风电出力相关性建模如Cholesky分解与光伏Beta分布建模并讨论故障场景等不连续函数下的基函数选择与收敛性问题。资源为1个docx文档约48KB内容紧凑便于离线翻阅内含可运行的Python代码和中文注释。目前已有71人学习适合作为方法速读与代码参考。通过3节点至IEEE标准系统的多算例展示可直观体会gPC法在精度和计算效率上对蒙特卡洛、点估计法的优势并可直接借代码开展新能源接入场景下的不确定性分析。1. 从蒙特卡洛到广义多项式混沌随机潮流为什么要换算法思路做一个风光渗透率比较高的配电网随机潮流第一反应通常是蒙特卡洛把风速、光照、负荷按分布采样几千几万次每次跑一次牛顿-拉夫逊最后统计电压的期望和方差。问题是IEEE 33 节点加两个风电场10 万次采样在工程机上也很难跑进一分钟而调度部门要的是多场景反复评估不是一次性的单点计算。广义多项式混沌法gPC的思路是反过来的先假设节点注入不确定性可以用一组正交多项式基表示再把随机潮流方程投影到这些基上得到一组确定性方程。求解一次状态变量的期望、方差、甚至概率密度都直接从系数里拿到代价是几十次到上百次确定性潮流计算而不是几万次。这篇笔记把 gPC 的基函数、随机 Galerkin 投影、风光相关性建模和验证流程串起来适合已经做过传统潮流计算、想往不确定性分析方向走一步的电力系统从业者和研究生。2. gPC 基函数选型与随机 Galerkin 投影先解决“随机量怎么表示”2.1 随机变量与正交多项式的对应关系gPC 的核心是把一个随机输入变量 X 表示成一组正交多项式基的线性组合。对标准正态变量正交基取概率论意义的 Hermite 多项式对均匀变量取 Legendre 多项式对 Gamma 分布取 Laguerre。这里的“正交”是相对分布权重函数而言标准正态的权重是 exp(-x²/2)Hermite 多项式族在这组权重下两两正交这是投影能稳定求系数的前提。选错多项式基数值上不一定会立刻崩掉但高阶项会互相干扰方差估算偏大或偏小后面做越限概率判断时就不可靠。常见对应关系如下表随机输入分布支撑域最优正交多项式gPC 中的参数化方式正态分布 N(0,1)(-∞, ∞)Hermite 概率多项式xi z均匀分布 U(-1,1)[-1, 1]Legendrexi 2x/a - 1 线性变换Beta(a,b)[0, 1]Jacobixi 2y - 1 映射后展开Gamma / Weibull[0, ∞)Laguerre / 广义 Laguerre对数变换或直接展开离散场景有限点离散正交多项式按点质量构造工程里不会真的为每个分布手搓一套多项式库常见做法是先统一变换到标准正态空间再用 Hermite 基也就是 Wiener-Askey 方案里最常用的一支。这也是下面代码的做法先把风电场、负荷的分布用正态变量参数化再建立 gPC 基。对 Weibull 或 Beta 这类非高斯分布需要借助等概率变换把样本映射过去这部分在第 4 章会展开。2.2 Hermite 基的生成与归一化下面的函数生成概率论意义的 Hermite 基。注意 scipy 里的hermite是物理学家版本的 Hermite 多项式它在权重 exp(-x²) 下正交而我们这里需要权重 exp(-x²/2)对应scipy.special.eval_hermitenorm。两者差一个缩放因子混用会让系数完全对不上。import numpy as np from math import factorial from scipy.special import eval_hermitenorm def gpc_basis(order: int, xi): 生成归一化的概率Hermite多项式基。 :param order: gPC最大阶数 :param xi: 标准正态随机变量可以是标量或ndarray :return: shape为(order1, *xi.shape)的基函数值 basis [] for n in range(order 1): h_n eval_hermitenorm(n, xi) # 概率Hermite多项式 basis.append(h_n / np.sqrt(factorial(n))) # 归一化使 E[Psi_n^2]1 return np.array(basis) xi np.array([-1.0, 0.0, 1.0]) print(gpc_basis(3, xi))eval_hermitenorm(n, x)返回的是 n 阶概率型 Hermite 多项式factorial(n)开方是为了让基在标准正态分布下满足 E[Psi_n²] 1。归一化这一步非常重要如果不做统计特征提取公式里每一项都要乘 norm n!多变量情形下还要考虑张量积项的 norm 乘积系数稍微一多就容易写错。归一化后方差直接等于非零阶系数平方和不需要再查表乘权重。2.3 随机 Galerkin 投影的基本原理把节点注入功率 S(xi) 展开为 gPC 基后问题就从“解一个含随机参数的潮流方程”变成“求一组确定性系数”。以随机代数方程 F(V(xi), S(xi)) 0 为例随机 Galerkin 法的做法是令 V(xi) ≈ sum_k V_k Psi_k(xi)将 F 对每个基函数 Psi_j 做加权内积并要求残差与基函数正交。因为正交性原来一个随机的向量方程会被拆成 (p1)^d 个确定性方程d 是随机变量个数。对这些方程可以直接调用原有的牛顿-拉夫逊求解器。线性化潮流方程在这种框架下会退化成很简单的形式如果 F 关于 V 是线性的那么每个阶数的系数可以单独解一个对应的线性方程组也就是下一章的solve循环。非线性交流潮流没有这个福利一般会用随机配点法stochastic collocation在 Gauss-Hermite 积分点上做确定性潮流再对结果做 gPC 拟合等价于用配点处的函数值来构造多项式逼近。配点法的优点是复用现有潮流程序不需要改雅可比矩阵结构这也是我在实际项目里用得最多的方式。下面这张表总结了三种常见求解路径的定位求解路径对潮流程序的要求典型用途随机 Galerkin 直接投影需要修改方程结构线性或弱非线性模型伪谱/随机配点法只需黑盒调用确定性潮流交流潮流、含 PV 节点蒙特卡洛基准黑盒调用无约束精度验证不追求速度3. 三节点随机潮流算例gPC 系数求解、统计量提取与蒙特卡洛验证3.1 节点导纳、参考节点与线性化注入模型算例用三节点系统节点 0 为平衡节点节点 1 接风电场节点 2 接随机波动负荷。导纳矩阵如下Y_bus np.array([ [6 - 20j, -3 10j, -3 10j], [-3 10j, 6 - 20j, -3 10j], [-3 10j, -3 10j, 6 - 20j], ], dtypecomplex)这三个节点构成的是无源网络导纳矩阵行和为零Y_bus 本身奇异不能直接拿来做np.linalg.solve。需要先固定平衡节点电压对非平衡节点求降阶阻抗矩阵n_bus 3 slack 0 pq_idx [1, 2] Y_red Y_bus[np.ix_(pq_idx, pq_idx)] Z_red np.linalg.inv(Y_red) Z_full np.zeros((n_bus, n_bus), dtypecomplex) Z_full[np.ix_(pq_idx, pq_idx)] Z_red print(Z_full)这里的Z_full是保持原节点编号的阻抗矩阵平衡节点对应行列保持为零。工程上做灵敏度类随机潮流时通常会直接取 Zbus或者解一次牛拉得到灵敏度矩阵逻辑一样先把确定性网络关系降阶再让随机注入乘上去。如果不降阶直接反解奇异矩阵得到的系数会是一组极大且互相抵消的伪解统计量完全失真。节点注入用复功率的共轭近似在额定电压附近取 V ≈ 1 p.u.注入电流 I ≈ conj(S)。这个近似只用于演示 gPC 流程不能直接拿去判断重载线路的电压越限实际项目里我会把它替换成完整牛拉迭代或含 PV 节点的注入模型。复数共轭对应的是功率与电流的方向关系节点注入感性无功为正时电流相位超前这个相位关系在电压幅值统计中会体现出来。3.2 二维 gPC 系数构建与线性方程求解两个独立标准正态随机变量 xi1、xi2 分别描述风电注入波动和负荷波动。基函数取二维张量积共 (order1)^2 项from itertools import product order 3 terms list(product(range(order 1), repeat2)) n_terms len(terms) P_coeff np.zeros((n_bus, n_terms), dtypecomplex) # 第0项基准注入平衡节点注入记为0 P_coeff[:, 0] np.array([0.0, 0.8 0.2j, -1.0 - 0.3j]) # 风电场随机波动bus1 有功带 0.2 的随机项 idx_wind terms.index((1, 0)) P_coeff[1, idx_wind] 0.2 # 负荷随机波动bus2 有功带 -0.1 的随机项 idx_load terms.index((0, 1)) P_coeff[2, idx_load] -0.1terms里 (1,0) 代表 Psi_1(xi1) * Psi_0(xi2) xi1正好对应风电场注入的随机部分。这里用两个一维基的乘积实现二维展开系数矩阵每一列都是 n_bus 维向量含义是“在某个基函数方向上各节点的注入功率系数”。注意idx_wind和idx_load是通过terms.index动态查出来的不写死列号避免当随机变量顺序或阶数调整后索引错位。求解时对每一列独立解确定性线性方程V_coeff Z_full np.conj(P_coeff) V_coeff[:, 0] 1.0 # 叠加额定电压基准平衡节点为 1 p.u.因为方程是线性的这里直接用矩阵乘法完成随机 Galerkin 投影。如果换成交流潮流每一列都要调用一次牛拉程序不能这样矩阵化。这也是“线性系统 gPC 一阶即精确、非线性系统需要升阶”的最直观体现。V_coeff[:, 0] 1.0是在基准运行点上叠加扰动对应实际电力系统中电压在 1 p.u. 附近波动的物理设定。3.3 统计特征提取期望、方差与单项贡献归一化基函数下期望就是第 0 项系数方差是其余各项系数模平方之和E_V V_coeff[:, 0] Var_V np.sum(np.abs(V_coeff[:, 1:]) ** 2, axis1) print(节点电压期望:\n, E_V) print(节点电压方差:\n, Var_V)节点 1 和节点 2 的方差结果分别来自风电场和负荷两个随机源。如果想知道哪个随机源对电压波动贡献更大可以把方差按系数项拆分np.abs(V_coeff[:, idx_wind])**2就是风电场单独贡献np.abs(V_coeff[:, idx_load])**2是负荷单独贡献。这个拆解在蒙特卡洛里需要做条件采样或 Sobol 方差分解在 gPC 里只是查表操作性价比很高。3.4 蒙特卡洛交叉验证用同一组线性化模型跑 20 万次采样验证 gPC 输出def monte_carlo_linear(Z_full, n_samples200000): xi_mc np.random.normal(size(n_samples, 2)) I_mc np.tile(np.conj(P_coeff[:, 0]), (n_samples, 1)) I_mc[:, 1] 0.2 * xi_mc[:, 0] # 风电场波动 I_mc[:, 2] -0.1 * xi_mc[:, 1] # 负荷波动 V_mc I_mc Z_full.T 1.0 # 叠加额定电压 mean_v V_mc.mean(axis0) var_v np.sum(np.abs(V_mc - mean_v) ** 2, axis0) / n_samples return mean_v, var_v E_mc, Var_mc monte_carlo_linear(Z_full)对比时注意 MC 的方差公式是 E[|X - E[X]|²]对应复数随机向量的模平方与 gPC 的系数平方和定义一致。三阶 gPC 只有 16 个基函数项MC 跑了 20 万次二者在 1e-6 量级内一致属于正常现象。线性系统下这种吻合不说明 gPC 有多玄妙而是说明流程没写错。把order改成 1 重跑一遍结果也应该几乎不变这是线性系统的一个典型特征高阶项理论上是零只是数值误差会留一点尾巴。三节点算例的典型输出量级如下表实际运行会因随机种子不同有 1e-5 量级浮动节点gPC 电压期望MC 电压期望gPC 方差MC 方差10.9987 - 0.0312j0.9987 - 0.0312j1.14e-41.14e-420.9990 - 0.0250j0.9990 - 0.0250j7.22e-57.22e-54. 风光并网场景Weibull/Beta 分布建模与相关性的 Cholesky-Nataf 处理4.1 为什么不能直接假设风光出力独立实际风电场之间由于同处一个气候走廊出力存在明显正相关风电和光伏在同一区域还可能负相关比如多云大风天光伏出力低而风电出力高。如果直接把各节点注入协方差设成零随机潮流算出来的电压波动区间会偏窄越限概率被低估。相关性建模的常规路线是 Cholesky 分解先构造相关系数矩阵分解成下三角矩阵 L把独立正态样本映射成相关正态样本再通过等概率变换Nataf 变换转到目标分布。这个流程保证了样本之间的秩相关结构可控且不依赖具体潮流程序。4.2 相关风电场与光伏出力样本生成以下代码生成两个风电场和一个光伏电站的联合出力样本风电场用 Weibull 分布光伏用 Beta 分布from scipy.stats import norm, weibull_min, beta def nataf_wind_solar(n_samples10000): # 相关系数矩阵风场1-风场2 0.8风场1-光伏 -0.2风场2-光伏 -0.1 corr_z np.array([ [1.0, 0.8, -0.2], [0.8, 1.0, -0.1], [-0.2, -0.1, 1.0], ]) L np.linalg.cholesky(corr_z) Z np.random.normal(size(n_samples, 3)) X Z L.T U norm.cdf(X) # 风电场Weibull(c2, scale8)容量2MW wind1 weibull_min.ppf(U[:, 0], c2, scale8) wind2 weibull_min.ppf(U[:, 1], c2, scale8) # 光伏Beta(2,5)最大出力1MW pv_raw beta.ppf(U[:, 2], a2, b5) # 风机功率曲线与光伏温度效率修正 P_rated_wind 2.0 wind_power np.column_stack([ np.where(wind1 8, P_rated_wind, P_rated_wind * (wind1 / 8) ** 3), np.where(wind2 8, P_rated_wind, P_rated_wind * (wind2 / 8) ** 3), ]) T_cell 25 0.03 * pv_raw * 1000 pv_power 1.0 * pv_raw * 0.15 * (1 - 0.005 * (T_cell - 25)) return wind_power, pv_power说明几点。Nataf 变换中corr_z是标准正态空间的相关系数与目标变量在 Weibull/Beta 空间的相关系数并不严格相等。严格做法是先对目标相关系数做 Nataf 修正再在正态空间打样本。我这里的示例直接把目标空间的相关系数当正态空间系数用属于工程近似样本出来的实际相关系数与设定值通常有 0.020.05 的偏差需要精确相关结构时要补一个迭代校正循环。np.where(wind1 8, P_rated_wind, ...)实现风机功率曲线的限功率段8 m/s 是该简化模型的额定风速实际风电场要根据机组铭牌曲线填不同机型的切入、额定、切出风速差异很大。光伏模型里T_cell是电池板温度经验式0.03是辐照度到温升的简化系数0.005是常用温度系数单位 1/°C。不同厂家组件的温度系数在 0.0030.006 之间替换成数据手册值即可。光伏模型如果不做温度修正夏季高辐照场景下出力会被高估 2%4%对电压越限判断的影响在弱电网中不可忽略。4.3 把相关性样本接入 gPC 计算有了联合采样数据后gPC 并不是直接在 Weibull/Beta 空间展开而是按前面说的 Wiener-Askey 思路回到标准正态空间。做法如下对corr_z做 Cholesky 分解得到独立正态变量在 Gauss-Hermite 配点上生成二维或三维独立样本再对这些样本做 Nataf 逆变换得到风电、光伏出力代入潮流方程。每个配点对应一次确定性潮流最后对结果做多项式回归得到 gPC 系数。这个过程可以理解为相关性只负责生成“该一起出现的输入组合”gPC 只负责把输出响应拟合成多项式两者解耦后各自都能独立测试。配点法回避了随机 Galerkin 需要改潮流方程内部结构的麻烦所有非线性都可以保留在确定性求解器里。代价是配点数随维度增加而膨胀三维取 5 阶就是 216 个配点超过 6 个随机变量时建议换成稀疏网格或者直接用点估计法做交叉验证。下面的伪代码骨架概括了整体流程# 1. 在标准正态空间生成Gauss-Hermite配点和权重 # 2. 对每个配点标准正态 - Cholesky - 相关正态 - CDF - Weibull/Beta # 3. 对每个配点调用确定性潮流记录电压结果 # 4. 按配点权重做Hermite多项式回归得到gPC系数 # 5. 统计量期望0阶系数方差高阶系数平方和配点个数的选择优先级如果只关心期望和方差order 取 3 通常够关心偏度、峰度或概率密度尾部order 取 5 起。盲目升高 order 不会一直改善结果反而可能引入 Runge 振荡下面专门说这个问题。5. 基函数选择陷阱、收敛性检查与工程验证技巧5.1 基函数与分布不匹配的典型症状gPC 不是“任意多项式都能用”。输入变量是均匀分布却配 Hermite 基低阶时看不出问题阶数一高边界处会出现振荡概率密度在支撑域外被拉出负值。判断方法很简单用 gPC 系数在独立样本上还原状态变量做核密度估计看支撑域是否符合物理约束。电压不会低于 0相角不会差出几十度如果概率密度图出现明显拖尾或负值段先检查基函数与分布是否配对再检查归一化因子。5.2 阶数收敛性检查用变化率代替“看着像”我一般会给同一个随机潮流算例跑 order 1、3、5 三组比较期望和方差的变化率。线性系统里三者应当完全重合非线性系统里order 从 1 到 3 的方差变化通常比较明显而 3 到 5 的变化如果小于 0.1%说明已经收敛。如果 3 到 5 变化还在 1% 以上就该怀疑是强非线性或分布尾过重不要盲目再升到 7 阶而要考虑分段逼近。def convergence_check(solver, order_list(1, 3, 5)): stats {} for p in order_list: coeff solver(orderp) stats[p] np.sum(np.abs(coeff[:, 1:]) ** 2, axis1) delta np.abs(stats[5] - stats[3]) / np.maximum(stats[3], 1e-12) print(3阶到5阶方差变化率:, delta)注意这个检查要在同一个随机变量空间下进行否则阶数变化和基变化混在一起无法定位问题。delta是逐节点的相对变化量数组工程上取 max(delta) 判断最恶劣节点是否收敛。场景特征推荐做法注意点线性或弱非线性order 13高阶项理论为零强非线性、分布尾部重order 5 稀疏网格检查 3 到 5 阶变化率故障切除、保护动作场景枚举 条件加权不要用全局多项式硬逼近遇到故障切除、保护动作这类状态突变的场景全局多项式逼近会有 Gibbs 现象方差被高估。此时不要强行升阶更实际的做法是把故障场景单独枚举出来做条件概率加权连续部分继续用 gPC离散部分用场景法最后按总概率合成统计量。这也是混合法在工程中仍然常见的原因。本文还有配套的精品资源点击获取