简介在新能源并网与风力发电系统评估中可靠性指标计算依赖可复现的风速序列与机组出力模型。其核心原理是用Weibull分布刻画风速长期概率特征再结合ARMA过程保留时间自相关性通过功率曲线将风速映射为有功功率随后以序贯蒙特卡洛抽样机组停运状态统计LOLP、LOLE、EENS等指标。技术价值在于把复杂物理过程压缩为可验证、可批量扫描的工程模型支撑场站选址、储能容量初筛和并网消纳分析。围绕通用风速模型、WTG发电模型与可靠性抽样三层链路简化建模与Python实现成为可落地的评估路径。1. 先别急着上蒙特卡洛把风电可靠性评估拆成两块可验证的模型拿到一个风电场几十台机组的出力数据直接扔进时序蒙特卡洛循环里跑可靠性指标通常跑不出可信的数字。原因不在算法而在前面两个环节没人管风速到底怎么生成风电机组WTG的出力又怎么从风速映射过来。一个完整的风电可靠性评估系统真正的骨架是「通用风速模型 WTG 发电模型 可靠性抽样」这三段链路前两段决定了你能不能复现真实出力第三段才决定指标算得准不准。这里讲的「简化模型」不是把物理机理砍掉而是把计算量从逐台机组逐叶片收敛压缩到能在单机上小时级跑完、且误差可控的形式。通用风速模型负责给出任意时长、任意步长的风速序列WTG 发电模型负责把风速转成有功功率并把故障停运状态叠上去可靠性层再用序贯抽样统计 LOLP、LOLE、EENS 这类指标。适合做新能源并网评估、场站选址比选、储能容量初筛的工程师也适合刚接触可靠性计算、需要一份可跑通示例代码的人。2. 通用风速模型从 Weibull 拟合到可复现的时间序列风速模型之所以叫「通用」是因为它要能适配不同风区而不是绑死在一组实测数据上。工程上最稳的组合是用双参数 Weibull 描述风速的长期概率分布用 ARMA 过程补偿时间相关性两者叠加后输出的序列既保留了风频特征又不会让相邻小时的风速看起来像白噪声。2.1 为什么是 Weibull 双参数而不是正态分布风速分布有明显的右偏特性正态分布会给出负风速概率物理上不成立。Weibull 分布的形状参数 k 和尺度参数 c 正好能覆盖从内陆低风速到海上高风速的整个谱系其概率密度函数为f(v) (k/c) · (v/c)^(k-1) · exp[-(v/c)^k]k 决定分布的陡峭程度c 近似对应 63.2% 累积概率处的风速。业内常用的经验判断是k 落在 1.82.2 之间时平均风速与 c 关系接近线性k 越大风速越集中。下面的表格给出几种典型风区的取值参考用于参数初值设定。风区类型形状参数 k尺度参数 c (m/s)年平均风速 (m/s)内陆平原1.82.05.56.55.05.8沿海滩涂2.02.37.08.56.27.5近海海上2.22.69.011.08.09.8山地峡谷1.51.86.07.55.36.6这些只是拟合起点实际项目仍要用测风塔至少一整年的 10 分钟数据做最大似然估计不能直接抄表。2.2 用 Python 拟合 Weibull 参数并做拟合优度检验scipy.stats.weibull_min提供了现成的极大似然拟合接口配合floc0可以把位置参数钉死在零避免它跑到负值区间去。拟合完必须看 KS 检验的 p 值p 值偏低说明该风区受季节或地形影响存在双峰此时要分季节或分扇区单独拟合。import numpy as np from scipy import stats # 假设已有 10 分钟级风速样本单位 m/s v np.loadtxt(wind_speed_10min.csv, delimiter,) # 极大似然拟合固定位置参数为 0 k_hat, loc_hat, c_hat stats.weibull_min.fit(v, floc0) # 经验分布与理论分布的 KS 检验 ks_stat, p_value stats.kstest(v, weibull_min, args(k_hat, 0, c_hat)) print(fk{k_hat:.3f}, c{c_hat:.3f}, KS p{p_value:.4f}) # 若 p 0.05按季节切分后重新拟合 seasons np.array_split(v, 4) for i, seg in enumerate(seasons): kk, _, cc stats.weibull_min.fit(seg, floc0) print(fseason {i}: k{kk:.3f}, c{cc:.3f})代码里的floc0是关键参数风电领域风速非负允许 loc 自由浮动会让拟合结果偏移。kstest返回的 p 值大于 0.05 才认为该风区可以用单条 Weibull 曲线概括否则要按季节或按风向扇区拆开每个子集单独拟合。这一步不做后面生成的年序列会系统性高估或低估大风时段占比。2.3 ARMA 修正让序列具备时间自相关性只按 Weibull 独立采样得到的是逐点独立的风速相邻时刻毫无关联用来算可靠性会低估连续低风速造成的长时段出力缺口。标准做法是先对归一化后的风速残差序列拟合 ARMA(1,1) 或 AR(2)再把模拟出的残差叠回 Weibull 分位数上。from statsmodels.tsa.arima.model import ARIMA # 把风速样本转成标准正态残差 z stats.norm.ppf(stats.weibull_min.cdf(v, k_hat, 0, c_hat)) z np.nan_to_num(z, nan0.0, posinf3.0, neginf-3.0) # 拟合 ARMA(1,1)即 ARIMA 的 (1,0,1) model ARIMA(z, order(1, 0, 1)).fit() # 模拟一年 8760 小时残差 z_sim model.simulate(nsimulations8760, anchorstart) # 逆变换回风速域 v_sim stats.weibull_min.ppf(stats.norm.cdf(z_sim), k_hat, 0, c_hat) v_sim np.clip(v_sim, 0, 40) # 物理上限保护 np.savetxt(wind_speed_sim.csv, v_sim, fmt%.4f)order(1,0,1)中的中间那个 0 表示不做差分因为残差序列本身已是平稳的。simulate的nsimulations按目标时长设置做年度可靠性评估就填 8760。np.clip上限 40 m/s 是防止逆变换在尾部数值溢出产生非物理值。经这一步得到的序列其均值、方差与原始样本接近同时保留了小时级的持续性这正是可靠性抽样需要的输入。2.4 风速模型最容易被忽略的两个坑第一个坑是时间尺度混淆。测风数据是 10 分钟级而可靠性评估常用 1 小时步长直接重采样取平均会削平风速波动导致出力曲线偏平滑、指标偏乐观。正确做法是在 10 分钟粒度上先做功率变换再对功率取小时平均。第二个坑是季节划分后参数拼接。分季节拟合得到的四条 Weibull 曲线直接首尾拼接会让年序列出现不连续跳变。常见做法是按月份采样每月内部再用该月的 k、c 生成保证统计特征平滑过渡。3. WTG 发电模型功率曲线四段式与多机型参数化WTG 发电模型的核心是一张功率曲线但直接拿厂家给的实测散点做插值会在切入和额定风速附近产生数值抖动。工程上更可靠的是把曲线抽象成四段解析式用切入风速、额定风速、切出风速、额定功率四个参数描述再按需叠加空气密度修正和尾流折减。3.1 四段式功率曲线的数学表达标准形式如下第一阶段风速低于切入风速 v_in出力为零第二阶段从 v_in 到额定风速 v_r出力按风速三次方近似上升第三阶段从 v_r 到切出风速 v_out出力恒定在额定值 P_r第四阶段超过 v_out机组顺桨停机出力归零。第二段的表达式为P(v) P_r · (v³ - v_in³) / (v_r³ - v_in³)用三次方而不是线性是因为风功率本身与风速立方成正比这一近似在中低风速段与实测曲线吻合度较高。切入和切出附近可以做小段线性过渡以削弱阶跃但大多数可靠性研究直接用硬阈值即可误差在可接受范围内。3.2 用 numpy 向量化实现分段功率曲线逐点循环处理 8760 小时乘几十台机组会非常慢用np.select一次算完整条序列是标准写法。下面示例同时处理单台机型和数组化输入。import numpy as np def wtg_power(v, p_rated, v_in, v_r, v_out, air_density1.225, rho_ref1.225): 四段式 WTG 功率曲线输入风速数组输出有功功率(kW) v np.asarray(v, dtypefloat) # 空气密度修正低密度空气功率下降按线性比例近似 density_factor air_density / rho_ref p_r p_rated * density_factor # 第三段额定出力 p_rated_seg np.full_like(v, p_r) # 第二段三次方过渡 p_trans p_r * (v**3 - v_in**3) / (v_r**3 - v_in**3) # 第一、四段停机 p_zero np.zeros_like(v) power np.select( [v v_in, (v v_in) (v v_r), (v v_r) (v v_out)], [p_zero, p_trans, p_rated_seg], defaultp_zero, ) return np.clip(power, 0, p_r) v_hour np.loadtxt(wind_speed_sim.csv) p_out wtg_power(v_hour, p_rated2000, v_in3.0, v_r12.0, v_out25.0) print(f年发电量 {p_out.sum()/1000:.1f} MWh, 容量因子 {p_out.mean()/2000:.3f})np.select的三个分支对应前三段default兜住超过切出风速的部分。air_density参数默认取标准大气 1.225 kg/m³高原或高温场站传实际值即可density_factor用线性近似替代严格的立方关系在 0.91.1 的密度比范围内误差小于 2%。最后的np.clip处理三次方过渡在 v_r 附近可能的数值超调保证出力不越过额定值。3.3 多机型参数表的组织方式一个风电场往往混装不同厂商、不同容量的机组用一张参数表统一管理比在代码里写死数值更利于批量扫描。下面给出参数结构示例每行代表一个可独立替换的机型模板。机型编号额定功率 (kW)切入风速 (m/s)额定风速 (m/s)切出风速 (m/s)轮毂高度 (m)WTG-A20003.012.025.090WTG-B30002.511.525.0100WTG-C45003.012.527.0110WTG-D60003.513.030.0115把这张表读成字典列表逐机型调用wtg_power即可生成各自的出力序列。轮毂高度用于风速外推如果测风塔高度与轮毂高度不一致要用幂律或对数律先把风速折算到轮毂高度幂律指数陆上取 0.143、海上取 0.11 是常见经验值。3.4 尾流与场站级折减的处理边界单机模型做好后场站级出力还要考虑尾流。简化评估里最常用 Jensen 模型把尾流区看成线性扩张的圆锥下游机组风速按折减系数衰减。但要注意尾流折减和可靠性评估的耦合不能太深否则每换一次风向就要重算全场流场计算量会失控。我一般会先按主风向扇区预先算好每台机组的平均尾流折减系数作为常数乘进该机组的出力序列。这样既反映了场站布局差异又把计算量压回 O(机组数 × 时长)量级。需要高精度时再上 CFD但那已经超出「简化模型」的范畴。4. 可靠性抽样与指标计算序贯蒙特卡洛的具体实现风速和出力序列准备好之后可靠性层要回答的问题其实很直接在机组随机停运的条件下可发出的总功率低于负荷的概率有多大、持续多久、缺多少电量。序贯蒙特卡洛通过逐时段抽样机组运行状态来回答这三个问题比非序贯方法更适合处理与风速相关的时间过程。4.1 机组两状态模型与状态持续时间抽样每台机组用运行和停运两个状态描述状态之间的转移服从指数分布参数由强迫停运率 FOR 决定。运行时长和修复时长分别按均值为 MTTF 和 MTTR 的指数分布抽样FOR 与两者的关系是FOR MTTR / (MTTF MTTR)常见火电机组的 FOR 在 0.030.08 之间风电机组因部件较多整机等效 FOR 通常取 0.040.10。抽样时用逆变换法-MTTF * ln(U)U 为 0 到 1 的均匀随机数。4.2 逐时段汇聚出力并判断失负荷把每台机组在各时段的可用状态乘以它的理论出力全场求和就得到系统可用发电能力。与负荷序列逐点比较低于负荷的时段记为失负荷事件。为了让评估有对比基准负荷可以用年最大负荷乘以典型日负荷系数曲线生成也可以直接读入实际负荷数据。import numpy as np def sequential_mc(wind_gen_df, load, for_rate, hours8760, seed42): wind_gen_df: 列机组, 行小时理论出力; load: 小时负荷序列 for_rate: 各机组强迫停运率数组 rng np.random.default_rng(seed) n_unit wind_gen_df.shape[1] mttf 1.0 / for_rate # 平均无故障时间(小时) mttr np.full(n_unit, 8.0) # 平均修复时间假定 8 小时 avail np.ones((hours, n_unit), dtypebool) for j in range(n_unit): t, state 0, True while t hours: if state: # 运行时长抽样 dur int(-mttf[j] * np.log(rng.random())) 1 state False else: # 停运时长抽样 dur int(-mttr[j] * np.log(rng.random())) 1 state True avail[t:tdur, j] False t dur gen (wind_gen_df.values * avail).sum(axis1) deficit np.maximum(0, load - gen) lolp np.mean(gen load) lole np.sum(gen load) eens deficit.sum() return {LOLP: lolp, LOLE: lole, EENS: eens}mttf和mttr的数组化是为了支持混装机组。内层while逐台生成状态时间轴avail矩阵记录每小时每台机组是否可用。最后的三个指标LOLP 是失负荷概率LOLE 是期望失负荷小时数EENS 是期望缺供电量单位跟随功率单位。注意mttr这里简化为固定 8 小时精确评估应从检修记录统计实际分布但固定值对指标量级影响有限。4.3 收敛判据方差系数比迭代次数更值得看蒙特卡洛抽样必须回答「跑多少次才算够」。用固定迭代次数是常见的偷懒做法稳妥的方式是监控方差系数 β公式为 β σ / (μ · √N)N 为抽样年数。当 β 小于 0.05 时停止一般需要几百到几千个抽样年取决于系统规模和需求的精度。实现上把 4.2 的评估包一层循环每次换一个随机种子记录指标序列边跑边算 β。EENS 的收敛速度通常慢于 LOLP因为它们对稀有事件的敏感度不同。实际项目中我一般以 EENS 的 β 作为主收敛判据LOLP 作为辅助。精度要求β 阈值LOLP 典型所需年数EENS 典型所需年数粗筛0.10100200300500常规0.054008008001500高精度0.021500300030006000这组数值会随系统风电渗透率上升而增加渗透率超过 30% 后所需年数往往翻倍因为失负荷事件的随机性被风速波动放大了。5. 参数敏感性排序与模型自验证的实操技巧模型搭完容易让人相信结果难。最有效的手段不是堆更多机组而是做参数敏感性和自验证两件事前者告诉你哪些输入值得花精力测准后者防止你把明显的逻辑错误当成模型特性。5.1 用单参数扰动找出主导变量固定其余参数对 k、c、FOR、v_out 各做 ±10% 扰动观察 EENS 的变化幅度排序即为敏感性。经验上 EENS 对 c 的敏感度最高其次才是 FOR 和 v_out因为 c 直接决定风资源水平。下面的最小示例把扰动逻辑写成可复用的函数。import numpy as np def sensitivity(base_k, base_c, base_for, wind_series, load, delta0.1): results {} # 风速整体缩放近似模拟 c 的变化 for name, factor in [(c, 1delta), (c-, 1-delta)]: v_scaled wind_series.mean() * factor / wind_series.mean() * wind_series p wtg_power(v_scaled, 2000, 3.0, 12.0, 25.0) df np.tile(p, (n_unit, 1)).T # n_unit 台同型机组 results[name] sequential_mc(df, load, np.full(n_unit, base_for))[EENS] return results扰动方式的选取很关键c 的扰动用风速序列整体缩放来近似k 的扰动则要重新拟合 Weibull 分位数再逆变换不能简单缩放。FOR 的扰动直接改for_rate数组即可v_out 的扰动要在wtg_power里改阈值后重算出力序列。5.2 三层自验证物理边界、量级对比、极端输入物理边界检查最简单切入风速以下出力必须为零切出风速以上必须为零风速拉满时出力不得超过额定。任意阶段出现负出力或超额定都是代码逻辑问题。量级对比用容量因子做锚点陆上风电场 0.250.35、海上 0.350.50如果你的模型算出 0.7 以上基本可以断定功率曲线参数或风速序列出了问题。极端输入测试是把风速序列全部设为 0、全部设为切出风速以上、以及全部设为额定风速三种情况下 EENS 和 LOLE 应该有明确且合理的响应全零风速时 EENS 最大额定风速恒定时应为理论最小值。5.3 让评估跑得更快的两个工程技巧一是把机组按参数分组同型机组共用一条理论出力序列只在状态抽样时区分能把wind_gen_df的列数压到机型数量级。二是把wind_gen_df预先转成float32并用np.einsum或矩阵乘法替代逐列求和千台机组规模下前者能把内存占用砍掉一半后者能把汇聚步骤提速三到五倍。本文还有配套的精品资源点击获取