企业官网建设流程全解析

搞电力系统不确定性的朋友应该都有个共同体会：风光一多，传统确定性潮流算完就得出一个“标准答案”，但现场负荷和出力天天在变，那个单点结果基本只能当安慰剂用。我最近拿IEEE33节点系统当试验对象，把蒙特卡洛法跟风光出力模型绑在一起做概率潮流计算，跑完一堆场景之后发现，这活儿看起来简单，实际坑不少，但用好了信息量非常可观。这篇文章就把我整个折腾过程、代码逻辑、参数选择和踩过的坑完整记录下来，希望给正在搞不确定性分析的同学一点可参考的操作路径。

蒙特卡洛法在概率潮流里的定位其实很直接：对风光出力、负荷波动这些随机变量进行大量采样，然后反复求解确定性潮流，最后对节点电压、支路功率做统计分析，得到概率分布和越限风险。IEEE33节点系统结构经典、参数公开，用来做算法验证再合适不过。搭配风光出力模型后，不仅能模拟源端随机性，还能把光伏和风电的相关性、时序特性一并考虑进去，算是入门概率潮流最不折腾的组合。

1. 为什么拿IEEE33节点当试验田

1.1 不确定性分析到底在解决什么问题

传统潮流计算的前提是运行工况固定，负荷和出力都取确定值。但实际电网里分布式光伏随风速和光照波动、负荷随着季节和时段变化，任何一个输入量的偏离都会导致节点电压或支路潮流偏离计算值。一旦偏离超出设备限值，就可能出现过电压、线路过载甚至保护误动。所以我们需要回答的不是“某个时刻系统安不安全”，而是“在一年或者一天内，系统有多大概率越限”。

概率潮流的本质就是把输入随机变量的分布映射到输出量分布，量化风险。这里蒙特卡洛法是数理统计里最朴素也最鲁棒的思路，只要采样次数足够多，计算结果就会依概率收敛到真实分布。它不像解析法那样要做大量线性化假设，也不像点估计法那样丢失高阶信息，很适合前期探索和精度验证。

1.2 IEEE33节点的优势与仿真意义

IEEE33节点系统是配电领域特别常用的算例，拓扑是单电源辐射网，包含33个节点、32条支路，首端是平衡节点，负荷有有功和无功，还有联络开关可以处理重构场景。这个规模很讨巧，潮流收敛快、计算量低，但包含了支路末端电压偏低、轻载重载切换这类特性，足够展示概率潮流结果里的细节差异。

拿它做蒙特卡洛试验还有一个好处：公开参数里负荷基准值是3715kW和2300kVar，恰好可以在此基础上叠加风光模型。我们在节点18接入风电、节点22接入光伏，容量按渗透率设定，就能清楚看到随机源接入后，末端电压和馈线潮流的分布变化。这种算例复现门槛低，非常适合先跑通流程再迁移到更大系统。

2. 风光出力模型该怎么搭

2.1 风速与光照的概率特性

风光出力模型的核心是描述随机出力背后的概率分布。风速常用两参数威布尔分布，概率密度函数记作：

[ f(v)=\frac{k}{\lambda}\left(\frac{v}{\lambda}\right)^{k-1}\exp\left(-\left(\frac{v}{\lambda}\right)^k\right) ]

其中形状参数k和尺度参数λ决定风速的峰态和平均量级。实际应用中，如果这里采集的是某地区风电场风速数据，可以用最大似然估计或者经验公式拟合出参数。风电出力与风速之间通常用分段线性曲线表示，PWind在切入风速以下或切出风速以上为0，额定风速以上封顶，中间段按线性或二次插值。

光伏出力的随机性来自光照辐照度，一般用Beta分布描述，概率密度为：

[ f(r)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\left(\frac{r}{r_{max}}\right)^{\alpha-1}\left(1-\frac{r}{r_{max}}\right)^{\beta-1} ]

Alpha和Beta可以根据历史辐照数据的均值和方差反推。在温度波动不大时，光伏有功出力可以简化为标准光照条件下的线性折算。这样每个采样时刻就能生成一组与风速和光照对应的风光出力值。

2.2 相关性怎么处理

单独采样风速和光照会忽略它们之间的天然联系。典型城区中，风速白天偏大、夜晚偏小，而光照中午最大、傍晚为零，二者存在时段上的负相关。如果采样时完全独立，得到的场景组合可能过于极端或脱离实际天气过程。

处理相关性比较主流的方法是使用Copula函数，把边缘分布和相关性结构分开建模。常用有Gaussian Copula和Frank Copula，前者参数容易根据秩相关矩阵估算，后者对尾部依赖描述更精细。实操中，我先对历史风速、光照数据计算秩相关系数，然后用该系数构造相关矩阵，再生成服从该相关结构的标准正态随机数，最后做逆变换得到对应的风速和光照样本。这样生成的场景序列和原始数据统计特性高度一致。

2.3 参数选择与场景生成

风光出力模型里最容易翻车的是参数选择不靠数据而靠拍脑袋。我建议第一步先获取本地风电场和光伏电站的小时级历史数据，至少覆盖一年8760个小时，然后分场景（如按月份或按天气类型）拟合参数。如果没有真实数据，也可以参照公开研究中的典型参数：风速Weibull分布的k取2.0到2.5，λ取6到8；光伏Beta分布的α、β取4到6。关键是要在后续蒙特卡洛采样前，先对模型做统计检验，比如把采样生成的出力分布与原始数据直方图对比，偏差大就要重新拟合。

场景生成数量取决于蒙特卡洛规模。我一般先抽500次，画直方图看分布形态，再逐步增大到5000或10000次，观察均值和方差的收敛情况。代码里，风光场景采样通常放在主循环外层，一次生成全部样本再逐个代入潮流，比每次循环里重复调用分布函数更快。

3. 蒙特卡洛概率潮流计算实操

3.1 整体流程拆解

概率潮流分析的完整流程分四步：定义输入随机变量、生成样本场景、求解确定性潮流、统计输出结果。第一步要明确哪些节点是风电、哪些是光伏，负荷在哪个节点什么范围波动。第二步上一个章节已经说了，用Copula采样获得一组组风速光照样本，再换算成风机和光伏的有功功率输出。

第三步是核心计算环节，每组样本中，把风光出力替代对应节点的注入功率，同时给负荷节点叠加正态扰动，然后调用牛顿拉夫逊法或前推回代法求解潮流。辐射配电网用前推回代法宰得更稳，但为保持通用性，我用Matpower里的runpf函数，它可以指定每种潮流方法。

第四步统计结果，对所有样本求得节点电压幅值、支路有功/无功潮流后，计算均值、标准差、最大值、最小值，以及节点电压越上限/下限概率、线路过载概率。输出如果做成概率密度曲线或者累积分布曲线，会比单纯列数字更直观。

3.2 核心代码逻辑逐段讲解

实际工程中，我用Matlab做全流程，代码规模不大。下面这段是主循环的核心骨架，我加了注释说明每一步作用。先定义随机变量维数和样本矩阵：

% 初始化 N = 10000; % 蒙特卡洛样本数 wind_nodes = [18]; % 节点18接风电 pv_nodes = [22]; % 节点22接光伏 Pg = zeros(33, N); % 存储每个节点的注入有功样本 for i = 1:N % 1. 抽取Copula相关风速和光照 U = mvnrnd([0,0], [1, c; c, 1], 1); U = normcdf(U); % 转换到均匀分布 v_sample = wblinv(U(1), lambda_w, k_w); % 风速 r_sample = betainv(U(2), alpha_pv, beta_pv); % 光照 % 2. 根据风速-出力转化曲线计算风电功率 P_wind = wind_curve(v_sample); % 自定义函数 P_pv = pv_curve(r_sample); % 自定义函数 % 3. 叠加负荷波动（正态分布） load_noise = normrnd(1, 0.05, 33, 1); P_load = base_load .* load_noise; % 4. 构造各节点注入功率 P_net = -P_load; % 负荷为负注入 P_net(wind_nodes) = P_net(wind_nodes) + P_wind; P_net(pv_nodes) = P_net(pv_nodes) + P_pv; % 5. 保存该样本的节点注入 Pg(:, i) = P_net; end

这里的核心细节是把“随机数生成”和“潮流计算”解耦。先生成全部输入样本，再用矩阵运算或并行循环依次求解潮流，这样可以利用Matlab的parfor并行加速。如果按照边采样边潮流的老写法，每次循环里都要初始化潮流参数，实际耗时能差一个数量级。

风速转功率的wind_curve函数在我这里用的是典型三段式曲线：

function Pw = wind_curve(v) v_in = 3; v_out = 25; v_rated = 12; P_rated = 500; % 单位kW if v < v_in || v > v_out Pw = 0; elseif v < v_rated Pw = P_rated * (v - v_in) / (v_rated - v_in); else Pw = P_rated; end end

光伏输出采用线性折算，光照超过额定辐照度就限幅。这些函数都要单独测试一下边界情况，尤其是风速在切入风速附近以及光照接近零的情况，否则会产出奇怪的负功率或超大功率。

潮流计算部分，我对每列样本调用一次runpf，并记录电压和潮流结果。值得注意，Matpower默认情况下如果潮流不收敛会直接报错中断。我采用try-catch结构，记录不收敛场景数量，并把不收敛场景置为上一次成功结果或丢弃。如果采样场景导致潮流不收敛次数太多，大概率是生成的出力或负荷波动范围设置过宽，需要回头检查分布参数。

V_amp = zeros(33, N); S_branch = zeros(32, N); conv_flags = zeros(1, N); parfor i = 1:N mpc = loadcase('case33bus'); % IEEE33节点的标准潮流数据 mpc.bus(:, 2) = mpc.bus(:, 2) - Pg(:, i) / 1000; % 注意单位转换，Matpower基准功率单位是MW mpc.bus(:, 3) = mpc.bus(:, 3) - Qg(:, i) / 1000; % 无功类似处理 try opt = mpoption('pf.alg', 'NR', 'out.all', 0); result = runpf(mpc, opt); V_amp(:, i) = result.bus(:, 8); S_branch(:, i) = abs(result.branch(:, 14) + 1i * result.branch(:, 15)); conv_flags(i) = 1; catch conv_flags(i) = 0; end end

3.3 关键参数设置与收敛性判断

蒙特卡洛法的精度只受到两个因素影响，一个是样本量，另一个是输入分布和实际物理过程的匹配程度。样本量不是越大越好，因为当前推回代法在这种配电网规模下，10000次计算可能在几分钟内完成，但很多情况下5000次已经足够稳定。

判断样本量是否足够的常用指标是输出均值随样本数的收敛曲线。我一般会设定一个滑窗，计算前200个样本的平均电压，之后每增加100个样本更新一次，看数值波动是否低于0.0001pu。如果低于这个阈值，说明再多样本也只是把尾部分布描得更细，不会改变整体结论。如果想精确预测极端场景，比如概率5%以下的越限，样本量需要明显增大，至少到20000次以上，采用重要抽样或拉丁超立方采样会更划算。

参数设置里还有个容易忽略的点：负荷波动标准差。很多人习惯给所有节点统一加5%波动，但实际配电网里不同节点的负荷性质差异很大，居民负荷波动大，商业负荷波动小。如果能拿到节点类型信息，最好分节点设置不同标准差。实在没有数据，至少也要让负荷波动均值归一化后保持原有基准值平衡，否则会造成系统总负荷偏差过大，让潮流结果整体偏移。

4. 结果分析与坑位记录

4.1 典型输出与分析指标

跑完10000个场景后，我习惯先输出节点电压幅值的概率分布。可以画成盒须图，横轴是33个节点，纵轴是电压标幺值，每个节点展示最小、四分之一分位、中位数、四分之三分位和最大值。通常在末端节点18和22附近能明显看到电压分布更宽，这是接入风光后的典型特征：末端电压不仅受负荷波动影响，还受随机出力变化影响，区间范围比首段节点大得多。

另一个值得分析的是支路潮流越限概率。分支潮流越限概率定义为P(Si > Simax)，其中Simax是支路容量。当风电大发、负荷低谷时，潮流反向或达到峰值，更容易触发越限。我计算出的结果往往能直观显示是哪几条联络支路风险最大，这对接下来的网络重构很有参考价值。分析层面还可以进一步做灵敏度分析，比如计算电压越限概率对风电接入容量的灵敏度，用来评估并网容量上限。

4.2 蒙特卡洛法常见坑与替代思路

第一个坑是伪随机数种子不一致。同一套代码，如果不同人运行没有固定随机种子，结果差异可能很大。所以在发布代码或做对比实验时，务必用rng(固定数字)初始化，保证可复现。

第二个坑是忽略相关性。独立采样导致风电和光伏同时出高功率的场景发生频率偏高，会让系统高渗透率场景被高估，越限概率被放大。我试过Copula修正前后，末端电压越上限概率能差出近3倍，这个差异足以影响并网决策。

第三个坑是样本量不足却强行看极值。想比较节点电压的最大值和最小值，500次采样得到的极值根本不可信，因为极值依赖稀有的极端场景。要稳定估计尾部，真得使用分层采样技术。我后来把普通蒙特卡洛换成拉丁超立方采样，在同一样本量下，电压标准差和越限概率的方差明显减小，代码改动不大，就是在采样前对均匀随机数做了一次排序。

第四个坑是潮流求解器设置不对。有些开源潮流工具对配电网多平衡节点的处理并不好，尤其当风电场容量占比高时，如果只把大电网平衡节点设为唯一平衡源，系统很容易因为无功不足而迭代不收敛。实际处理方法要么换成分布式的松弛节点，要么在风光场站配置无功补偿模型。

4.3 实用建议

如果只是想快速评估方案优劣，没必要一上来就跑10000次。先跑200次看趋势，再针对高风险场景加密采样，这就是所谓序贯蒙特卡洛思路。具体做法是：第一轮粗采样后用KDE估计越限区域的概率密度，然后在越限区域增加权重重新采样，得到的估计方差最小。

工程实践中，我还建议把概率潮流结果和时序生产模拟结合。单纯随机采样得到的是“不考虑时间顺序”的概率结果，但储能调度、爬坡需求都要看时序性。我在项目里的做法是先按小时划分场景，在每个小时段内做蒙特卡洛，这样得到的逐时概率分布能直接给调度员提供风险预警窗口。

结果展示时，别只给均值。把“5%分位数”“95%分位数”和“越限概率”三件套做成表格或热力图，运维人员看到才知道风险集中在哪条馈线、哪个时段。我这里整理了自己用的一组典型结果格式，给大家参考：

支路过载概率可以单独整理成另一张表，标明支路编号、线路容量、最大负载率和过载概率。通过这些表格，规划人员能直接找准薄弱环节，而不是大海捞针。

最后再分享一点个人经验：做这类仿真，最消耗精力的不是潮流算法，而是输入数据的清洗和概率模型的校验。我前几版跑出来结果特别夸张，后来一查是光照辐照度没有做昼夜归属，把夜间零辐照和白天高辐照混在一起，导致Beta分布拟合结果失真。花点时间把历史数据按时间切片、按天气聚类，再拟合参数，出来的概率潮流才真正有工程参考价值。概率潮流不是炫技舞台，它解决的是“系统风险能不能被看到”的问题。把数据模型做扎实，蒙特卡洛这把牛刀才不至于切碎了豆腐。