搞过风光联合出力建模的朋友应该都有体会风速和光照强度这两类数据直接拿平均值去代表典型工况基本是自欺欺人。风速分布是明显右偏的光照强度在晴天和阴天完全是两种形态用均值、方差这种二阶统计量根本描述不了完整特征。我做了几个风电、光伏出力特性相关的项目后最稳定的做法还是那一套——风电用两参数Weibull分布拟合风速光伏用Beta分布拟合辐照度标幺值最后用Matlab把参数估计、拟合检验、场景生成串成一条完整链路。这篇文章把整套流程的来龙去脉、代码实现和坑都写出来适合正在做新能源出力建模、电力系统概率潮流、可靠性评估或者储能容量规划的同学参考。1. 为什么风电偏偏用Weibull、光伏偏偏用Beta概率建模的物理逻辑很多人第一次看到这个选题会问分布函数那么多正态分布、对数正态分布都行为什么偏偏是这两个这个问题的答案不在数学里而在物理过程里。1.1 风速的统计特征与Weibull分布的天然匹配风速是一个非负的随机变量而且实际测风数据几乎都有一个共性低风速和高风速出现的概率都不算高中等风速出现的概率最高概率密度曲线呈右偏形态拖尾较长。这种形态用正态分布去拟合第一个问题就是正态分布允许负值出现物理上完全说不通第二个问题是正态分布对称拟合出来的峰值位置和实测数据的偏态明显不匹配。两参数Weibull分布的密度函数是这样写的f(v) (k / c) × (v / c)^(k-1) × exp[-(v / c)^k]其中 k 是形状参数决定曲线的偏态程度和峰值形态c 是尺度参数大致对应风速的量级。当 k 在 2 到 3 之间时曲线形态和大多数陆上风电场的测风数据高度吻合。工程上从20世纪70年代开始就用Weibull分布做风资源评估IEC标准里也大量使用经过几十年验证可以说是风力发电领域公认的标准模型。实际应用中你会发现k 还会携带额外信息k 越小说明风速波动越剧烈、分布越分散k 越大说明风速越集中资源越稳定。我在做风资源评估时经常用这个参数直接对比不同场址的发电稳定性比单纯看年平均风速直观得多。1.2 Beta分布与辐照度的映射逻辑光伏出力的核心驱动量是辐照度而辐照度在建模前通常要做标幺化处理即用当前辐照度除以标准测试条件STC下的辐照度一般取1000 W/m²得到一个介于0和1之间的值。Beta分布恰好是定义在[0,1]区间上的连续分布密度函数可以呈现U形、钟形、单调递增或递减等多种形态灵活度非常高。Beta分布的密度函数为f(x) [Γ(α β) / (Γ(α) · Γ(β))] × x^(α-1) × (1-x)^(β-1)其中 α 和 β 是两个形状参数。当 α 1、β 1 时曲线呈钟形当 α β 1 时退化为[0,1]上的均匀分布当 α 1、β 1 时呈U形。这种形态变化能力对光照建模特别关键——晴天条件下辐照度普遍偏高分布集中在0.7附近多云天气则可能呈现低辐照度频发甚至出现双峰现象。Beta分布可以通过调整两个参数逼近这三种典型状态。更重要的一个理由是Beta分布对任意[0,1]上的连续随机变量都是一种良好的一阶近似而光伏出力标幺值恰好就属于这一类。相比之下对数正态分布在辐照度较高时拟合得还行但在接近0的区域往往拟合得很差而实际光伏出力在阴雨天和早晚时段的零出力占比相当大这个尾巴拟合不好整个模型的可信度就大打折扣。1.3 为什么要把两个分布组合在一起单一分布建模只是第一步实际工程里更关心的是风电和光伏同时作用下的系统总出力形态。风电的随机性和光伏的随机性在时间尺度上存在天然的互补关系——白天光伏出力高但风速可能相对平稳夜间光伏为零风电往往成为主力。把Weibull分布和Beta分布组合到一起本质上是在构建一个包含两个随机源的联合出力模型用来回答以下几个实际工程问题风光互补系统的总出力在什么区间波动概率密度形态如何系统出力低于某个阈值比如负荷需求的缺电概率有多大在给定置信水平下可被系统稳定消纳的风光出力是多少储能配置多大容量才能平滑掉一定比例的风险这类问题单独看风速分布或单独看辐照度分布都回答不了必须建立组合模型。最常用的做法是假设两个随机变量相互独立先分别拟合出各自的边缘分布再通过Monte Carlo抽样生成联合场景。独立性假设在小时级尺度上基本成立虽然日尺度上风电和光伏可能都受天气系统影响但工程上常用Copula理论做进一步修正这篇文章后面会展开讲。2. 风电Weibull分布的Matlab拟合从数据清洗到参数估计建模的第一步永远是数据清洗而不是参数拟合。我见过不少人在数据还没整理干净的情况下直接调用wblfit最后拟合出来的参数连自己都解释不了问题往往就出在数据质量上。2.1 预处理别让坏数据毁掉整个拟合风速数据最常用的来源是测风塔的10分钟平均数据或者气象站逐小时数据。拿到数据后我建议按以下顺序处理剔除超出物理范围的值比如风速大于60 m/s的记录基本是传感器故障或雷击干扰。检查零风速和负风速记录。零风速可能是静风期的正常值但如果某一天连续出现大量零值要确认是不是测风设备结冰或堵塞。处理缺失值。直接删除即可不需要插值——插值会平滑掉风速的随机波动导致拟合出的 k 参数偏大。明确建模时段。全年统一拟合和分季节拟合的结果差异很大我的建议是先画月风速曲线如果季节差异明显就按季节或按典型月份分开建模。预处理完之后把风速序列转成列向量并统计一下样本量。Weibull参数估计对样本量有一定要求实际经验是少于500个有效样本拟合结果的置信区间会非常宽做工程判断时要谨慎。2.2 三种参数估计方法与Matlab实现Matlab的Statistics and Machine Learning Toolbox里自带最大似然估计函数这是最推荐的方式。但如果你希望完全掌控计算过程或者需要深入理解参数估计的机理我建议把下面三种方法都掌握。方法一极大似然估计MLE这是最标准、统计性质最好的方法。Matlab中可以直接调用% 假设 v 是清洗后的风速列向量单位 m/s v v(~isnan(v)); v v(v 0); % 零风速不参与拟合可以根据业务需求决定 % wblfit返回 [k, c] [k_mle, c_mle] wblfit(v); % 如果想同时获取置信区间 [params, pci] wblfit(v); k_mle params(1); c_mle params(2); fprintf(Weibull MLE: k %.4f, c %.4f m/s\n, k_mle, c_mle);wblfit内部通过迭代求解对数似然函数的最大值。对数似然函数对 k 求导后是一个非线性方程Matlab默认用牛顿迭代法求解初值选得不好可能导致不收敛。如果你的数据分布很怪异可以先手动给定一个初值试试或者换下面两种方法估算一个初值。方法二最小二乘法线性化Weibull分布的累积分布函数为F(v) 1 - exp[-(v/c)^k]做两次对数变换ln[-ln(1 - F)] k·ln(v) - k·ln(c)这是一个典型的直线方程斜率为 k截距为 -k·ln(c)。实际操作中先对风速排序用经验累积频率公式 F_i (i - 0.3) / (n 0.4) 计算每个点的累积概率然后做线性回归。vs sort(v); n length(vs); F ((1:n) - 0.3) / (n 0.4); x log(vs(:)); y log(-log(1 - F)); % 一次线性回归 p polyfit(x, y, 1); k_lsq p(1); c_lsq exp(-p(2) / k_lsq);这段代码的好处是直观、可解释性强而且完全不需要任何工具箱。缺点是经验累积频率在两端会失真尾部的拟合精度较低。实际使用中我会把线性回归的结果作为MLE的初值两者结合效果最好。方法三矩估计法Weibull分布的均值 μ 和标准差 σ 与参数的关系是μ c·Γ(1 1/k)σ² c²·[Γ(1 2/k) - Γ²(1 1/k)]通过数值求解这两个方程即可得到 k 和 c。Matlab实现需要fsolvemu mean(v); sd std(v); mom_egn (x) [c_est(x, mu, sd) - mu; sd_egn(x, mu, sd) - sd]; % 这里直接给出数值求解更简洁的形式 obj (k) (gamma(1 2/k) ./ gamma(1 1/k).^2 - 1) - (sd / mu)^2; k_mom fzero(obj, 2); c_mom mu / gamma(1 1/k_mom);矩估计法计算速度快适合数据量大、需要批量处理的场景但效率不如MLE。三种方法的对比如下方法精度计算量适用场景MLE最高中等一般工程分析首选最小二乘线性化中等低快速估算、MLE初值矩估计中等低批量处理、实时计算2.3 拟合优度验证直方图、Q-Q图与K-S检验拟合完成后不能直接拿去用必须验证。我的标准验证流程是三步。第一步把概率密度曲线叠加到频率直方图上用肉眼检查形态是否吻合。注意直方图要用频率密度频数/bin宽度/总样本数而不是频数否则和概率密度函数不在一个量纲上看起来会差很多。figure; histogram(v, 30, Normalization, pdf, FaceColor, [0.7 0.7 0.7]); hold on; v_axis 0:0.1:max(v)*1.1; pdf_weibull wblpdf(v_axis, k_mle, c_mle); plot(v_axis, pdf_weibull, r-, LineWidth, 1.5); xlabel(风速 (m/s)); ylabel(概率密度); legend(实测频数密度, Weibull拟合);第二步画Q-Q图。figure; qqplot(v, weibull_dist);如果散点基本落在直线上说明分位数对齐良好尤其是尾部。Q-Q图比直方图更能反映尾部拟合质量这对评估高风速段的极端出力特别重要。第三步做K-S检验。[h, p] kstest(v, CDF, wblcdf(v, k_mle, c_mle));p值大于0.05意味着在95%置信水平下不能拒绝原假设即数据与Weibull分布一致。但我必须提醒一点K-S检验对样本量非常敏感几万条样本下几乎任何分布都会被拒绝所以当样本量很大时别盯着p值看而是看Q-Q图的偏离程度和实际工程误差。3. 光伏Beta分布的Matlab实现标幺化与参数求解光伏建模的核心变量是辐照度但在实际工程里大家更关心出力。Beta分布直接建模的对象是辐照度标幺值或者功率标幺值。这里有一个选择到底对哪个变量建模更合理我的经验是如果做的是长期规划分析直接对辐照度标幺值建模更稳定如果做的是并网运行分析对实际出力标幺值建模更方便因为已经包含了逆变器效率和限功率因素的影响。3.1 辐照度到功率的标幺化转换光伏出力与辐照度之间的简化关系可以写成P_pv P_rated × (R / R_stc) × [1 γ × (T_cell - T_stc)]其中 R 是水平面辐照度或组件面辐照度R_stc 1000 W/m²P_rated 是额定容量γ 是温度系数T_cell 是电池板工作温度。如果暂时忽略温度影响P_rated 的标幺值出力约等于辐照度标幺值。这份工作里我们主要关注辐照度分布的统计特性所以直接对 x R / R_stc 建模。需要注意的是实测辐照度数据里会掺杂大量夜间零值。如果直接把全天24小时的数据拿去拟合Beta分布零值的比例可能高达40%到50%这会导致 Beta 分布参数严重失真。我在实际操作中会做分时段处理只取白天有效日照时段比如早上6点到晚上19点之间或者按辐照度阈值 10 W/m² 筛选建模半夜的零值放到时序场景生成时再处理。这样才符合物理逻辑。3.2 矩估计与极大似然估计的实现Beta分布参数估计在Matlab里没有像wblfit一样的一行函数但可以通过betafit直接完成MLE估计也可以用矩估计手动求解。方法一极大似然估计% x 是辐照度标幺值取值在(0,1)之间 x x(x 0 x 1); % 去掉0和1避免对数计算发散 % betafit返回 [alpha, beta] params_beta betafit(x); alpha params_beta(1); beta params_beta(2); fprintf(Beta MLE: alpha %.4f, beta %.4f\n, alpha, beta);betafit内部也是迭代计算如果数据里包含0或1会直接报错或者给出无穷参数所以我习惯先做一次截断处理。方法二矩估计矩估计的好处是不依赖工具箱而且可以写成非常简洁的解析公式。Beta分布的均值和方差为μ α / (α β)σ² αβ / [(α β)²(α β 1)]反解得α μ × (μ(1-μ)/σ² - 1)β (1-μ) × (μ(1-μ)/σ² - 1)Matlab代码mu mean(x); var_x var(x); % 注意var_x 必须小于 mu*(1-mu)否则参数为负说明数据不适合Beta分布 denom mu * (1 - mu) / var_x - 1; alpha_mom mu * denom; beta_mom (1 - mu) * denom; if alpha_mom 0 || beta_mom 0 warning(矩估计得到的参数非正数据可能不适合Beta分布); end实测经验表明在样本量较大的情况下矩估计和MLE的结果差异很小。但如果辐照度数据比较干净且模型形态很标准优先用MLE如果只是想快速拿到一套参数做初步仿真矩估计足够了。3.3 边界样本与参数失真的处理光伏数据里最棘手的问题就是大量0和1的出现。辐照度标幺值为0的时间段占比很大而标幺值为1甚至大于1的情况也会出现——比如某些高海拔地区或者聚光光伏系统辐照度可能瞬间超过1000 W/m²。Beta分布的支持域是开区间(0,1)端点处概率密度为0当参数大于1时。直接把0、1扔进去拟合会有两个问题一是对数似然函数在端点处发散计算不稳定二是拟合结果会严重偏向端点导致中间段的形态贴合度很差。我的处理策略分三层第一层数据筛选。按业务逻辑决定0值怎么处理。如果是辐照度建模且研究重点是白天时段的光伏出力那么直接剔除0值模型只对有效日照时段的随机性负责。第二层端点微调。如果数据里包含少量1或接近1的值可以用一个极小的偏移量把它们拉回(0,1)区间比如x min(1 - 1e-6, max(1e-6, x))。这个做法对参数估计的影响可以忽略但能避免数值发散。第三层零出力注入模型。如果后续场景生成时需要模拟全天24小时出力不要在Beta分布内部塞0而是用混合模型以一定概率 p0 产生零出力以剩余概率从Beta分布采样。这样一个两步抽样逻辑既保持了Beta分布的形态又真实反映了光伏夜间和阴雨天的特性。% 夜间或低辐照时段比例 p0 sum(R 10) / length(R); % 辐照度阈值自行调整 % 有效时段Beta采样 x_eff betarnd(alpha, beta, N, 1); % 混合抽样 u rand(N, 1); x_sample zeros(N, 1); x_sample(u p0) x_eff(u p0);这种处理办法在后续做全年8760小时时序仿真时非常管用推荐直接用。4. 组合建模与场景生成从概率密度到实际出力序列边缘分布拟合完成后真正的重点来了——两个分布怎么组合组合之后怎么生成可以用于仿真计算的出力场景。这一节是整个项目最核心的工程环节。4.1 独立假设下的联合概率模型如果假设风速和辐照度相互独立那么联合概率密度就是两个边缘密度的乘积。这个假设在小时级时间尺度上有一定的合理性——风主要受大尺度气压系统驱动光照主要受云和大气透明度影响两者之间没有直接的物理耦合。但要注意在日尺度或更长尺度上天气系统可能同时影响风和光比如锋面过境时风速增大、云量增多此时独立假设会带来偏差。想要充分考虑相关性可以用Copula函数建模其中常用的有Gaussian Copula和t-Copula。这里先给出独立假设的实现相关性建模作为扩展放在后面。联合模型的意义在于它把一个双随机源的工程问题转化成了一台可以抽样的随机数生成器。有了这个生成器我们就可以做Monte Carlo模拟把风速、辐照度样本对同时生成出来再通过各自的功率转换曲线得到总出力。4.2 Monte Carlo场景生成完整代码下面给出一个完整的场景生成函数。输入是Weibull参数、Beta参数和系统容量参数输出是场景矩阵。function [P_wind, P_pv, P_total] generate_hybrid_scenarios(k_w, c_w, alpha_pv, beta_pv, ... P_wind_rated, P_pv_rated, N_scenarios, v_cut_in, v_cut_out, v_rated, pv_p0) % 输入参数说明 % k_w, c_w: 风速Weibull参数 % alpha_pv, beta_pv: 辐照度Beta参数 % P_wind_rated: 风电场额定功率 % P_pv_rated: 光伏电站额定功率 % N_scenarios: 场景数量 % v_cut_in, v_cut_out, v_rated: 风机切入/切出/额定风速 % pv_p0: 光伏零出力概率 % 1. 采样风速 v_sample wblrnd(k_w, c_w, N_scenarios, 1); % 2. 采样辐照度标幺值零出力混合 u rand(N_scenarios, 1); x_sample zeros(N_scenarios, 1); x_sample(u pv_p0) betarnd(alpha_pv, beta_pv, sum(u pv_p0), 1); % 3. 风速转风电出力 P_wind zeros(N_scenarios, 1); idx_region1 (v_sample v_cut_in) (v_sample v_rated); idx_region2 (v_sample v_rated) (v_sample v_cut_out); P_wind(idx_region1) P_wind_rated .* ... (v_sample(idx_region1).^3 - v_cut_in^3) ./ (v_rated^3 - v_cut_in^3); P_wind(idx_region2) P_wind_rated; % 4. 辐照度标幺值转光伏出力线性近似 P_pv P_pv_rated .* x_sample; % 5. 总出力 P_total P_wind P_pv; end关于风机功率曲线我用的是标准的线性三次方段模型即切入风速以下不出力额定风速到切出风速之间保持额定出力切出风速以上停机。实际工程中更精确的做法是查厂家提供的功率曲线表但做统计研究时这个模型已经足够。生成场景之后可以画一个散点图横轴是风电出力、纵轴是光伏出力颜色代表总出力的高低。从这个散点图可以直观看到联合出力的分布范围以及风光出力互补的状态在哪个区间分布最密集。scatter(P_wind / P_wind_rated, P_pv / P_pv_rated, 10, P_total / (P_wind_rated P_pv_rated), filled); colorbar; xlabel(风电出力标幺值); ylabel(光伏出力标幺值); title(风光联合场景分布);4.3 场景削减与典型场景选取Monte Carlo生成的场景数量如果只有几百个还好但如果要做一年8760小时的时序模拟每个时刻都生成大量场景计算规模会迅速膨胀。这时候需要场景削减——从大量场景里挑出一组具有代表性的“典型场景”用少量场景近似逼近原始联合分布。最常用的削减方法是K-means聚类。把每个场景看成二维空间里的一个点K-means把整个点集划分为K类每个类的质心就作为代表性场景各类的样本占比作为该场景的概率权重。% 假设场景矩阵 X 是 N行2列 [风电, 光伏] X [P_wind, P_pv]; % K-means聚类K根据实际需要选择一般10~50 k_idx 30; [cluster_idx, C] kmeans(X, k_idx, Replicates, 5, MaxIter, 500); % 计算每个类的比例作为权重 weights histcounts(cluster_idx, k_idx) / length(cluster_idx); % 典型场景为质心坐标 typical_scenarios C;削减后的场景集可以直接用于后续的随机潮流计算、储能容量优化或可靠性分析。实测下来30个典型场景就能在绝大部分情况下保持与原分布的前四阶矩基本一致计算量却降低了一个数量级。场景削减有个经验值先用肘部法则选K。画一下不同K对应的组内离差平方和随着K增加下降曲线出现明显拐点的位置就是合适的K。我一般不追求特别大的K因为场景削减的目的是保留“分布特征”而不是保留“所有细节”。5. 实操中反复踩过的坑与工程化建议最后这部分写点干货中的干货全是实际项目里踩过坑之后总结出来的。如果你按照前面的代码直接跑大概率能出结果但要做出可信的工程结论下面这几个问题必须提前处理。5.1 拟合误差不等于模型误差做拟合的人容易陷入一个误区追求拟合优度越高越好用各种指标轰炸。但工程建模的核心目标不是让曲线穿过每一个点而是让模型能正确回答业务问题。比如你做储能容量配置最关心的是出力低于负荷的概率和持续时长分布这时即使Weibull拟合的K-S检验p值低一些只要概率密度曲线在低出力区间的累计概率准确结果就是可用的。反过来如果K-S检验通过了但Q-Q图的尾部偏差明显说明模型对极端低风速或极端高辐照度场景描述不准确。这种情况下的概率密度“看起来”没问题但做可靠性评估时会严重低估极端场景的发生概率。我建议每次拟合完都同时看三个指标K-S检验p值、Q-Q图形态、以及你关心的业务区间的累计概率误差。三者结合才能判断模型是否真的可用。5.2 季节与天气分型比一味追求高精度更重要一个全年统一拟合的Weibull/Beta模型在统计上是“正确”的但在工程上往往是粗糙的。风资源有显著的季节差异比如我做过的一个西北地区的项目春季大风天的平均风速能比秋季高40%对应的k值也从2.1变化到3.2光伏出力更是跟季节日照时长和太阳高度角强相关。如果全年只用一个模型生成的场景会淡化季节差异导致储能配置结果偏乐观或偏悲观。我的做法是先按月或按季节聚类分时段拟合。更细致一点的做法是按天气类型分型从历史数据里识别出晴朗、多云、阴雨、大风四种典型天气分别对每个天气类型拟合Wind-PV联合分布然后用天气类型的出现概率做加权。这样建出来的模型能够生成更真实的“晴天但大风”、“阴天但微风”这类场景组合对系统风险评估的增益非常明显。5.3 代码效率向量化与并行工具箱如果你要生成几百万个场景做Monte Carlo模拟循环写法会非常慢。Matlab的随机数生成函数本身是向量化的但我见过不少人在场景转换环节写for循环比如对每个场景逐一计算风机功率。这里建议直接用逻辑索引批量计算性能差距能达到几十倍。% 低效做法 for i 1:N if v(i) v_rated p(i) P_rated; end end % 高效做法 p zeros(N, 1); p(v v_rated v v_cut_out) P_rated;如果再往后做多风电场、多光伏电站的联合建模单机串行计算会非常吃力。Matlab的Parallel Computing Toolbox里parfor是现成的场景之间如果完全独立直接循环并行即可。我最初做某个省级电网的风光出力时序仿真时用parfor把8760小时的场景生成从半小时压缩到了两三分钟这个优化非常值。5.4 参数估计对样本量的敏感性最后强调一个容易被忽略的问题样本量。比如你只有两个月的测风数据却想拟合出代表全年特性的Weibull参数结果必然偏颇。通常我要求至少有一整年的有效数据且时间分辨率达到小时级。如果数据不足宁可退回到用参考气象站的多年平均风速来标定尺度参数c也不要用短样本硬拟合。Beta分布的样本量要求同样很高而且辐照度数据的有效样本量往往比风速还要少因为夜间数据几乎全被筛掉了。假设你有一年小时级数据晚上6点到早上6点剔除后有效样本可能只有3000个左右勉强够用但置信区间仍然偏宽。这时候我建议把多年数据合并使用或者采用分季节建模来增加各季节的有效样本量。结尾的一点体会这套WeibullBeta的组合建模方法我前后用了小半年才稳定下来最大的体会是模型本身不复杂复杂的是数据处理的边界条件——什么时候该剔除零值、什么时候该分季节拟合、什么时候该用混合分布这些问题只能通过反复对比实测数据来形成判断。最后再分享一个小技巧每次拟合完把拟合参数和对应的K-S检验结果存成一个结构体变量方便不同方案之间对比。这个过程比跑通代码更有价值因为它逼着你去理解数据而不是仅仅让代码跑出数字。这套方法后续还可以往Copula相关性建模、动态参数时变模型方向扩展等我把相关性部分的代码也整理完了再单独写一篇分享。
搞过风光联合出力建模的朋友应该都有体会风速和光照强度这两类数据直接拿平均值去代表典型工况基本是自欺欺人。风速分布是明显右偏的光照强度在晴天和阴天完全是两种形态用均值、方差这种二阶统计量根本描述不了完整特征。我做了几个风电、光伏出力特性相关的项目后最稳定的做法还是那一套——风电用两参数Weibull分布拟合风速光伏用Beta分布拟合辐照度标幺值最后用Matlab把参数估计、拟合检验、场景生成串成一条完整链路。这篇文章把整套流程的来龙去脉、代码实现和坑都写出来适合正在做新能源出力建模、电力系统概率潮流、可靠性评估或者储能容量规划的同学参考。1. 为什么风电偏偏用Weibull、光伏偏偏用Beta概率建模的物理逻辑很多人第一次看到这个选题会问分布函数那么多正态分布、对数正态分布都行为什么偏偏是这两个这个问题的答案不在数学里而在物理过程里。1.1 风速的统计特征与Weibull分布的天然匹配风速是一个非负的随机变量而且实际测风数据几乎都有一个共性低风速和高风速出现的概率都不算高中等风速出现的概率最高概率密度曲线呈右偏形态拖尾较长。这种形态用正态分布去拟合第一个问题就是正态分布允许负值出现物理上完全说不通第二个问题是正态分布对称拟合出来的峰值位置和实测数据的偏态明显不匹配。两参数Weibull分布的密度函数是这样写的f(v) (k / c) × (v / c)^(k-1) × exp[-(v / c)^k]其中 k 是形状参数决定曲线的偏态程度和峰值形态c 是尺度参数大致对应风速的量级。当 k 在 2 到 3 之间时曲线形态和大多数陆上风电场的测风数据高度吻合。工程上从20世纪70年代开始就用Weibull分布做风资源评估IEC标准里也大量使用经过几十年验证可以说是风力发电领域公认的标准模型。实际应用中你会发现k 还会携带额外信息k 越小说明风速波动越剧烈、分布越分散k 越大说明风速越集中资源越稳定。我在做风资源评估时经常用这个参数直接对比不同场址的发电稳定性比单纯看年平均风速直观得多。1.2 Beta分布与辐照度的映射逻辑光伏出力的核心驱动量是辐照度而辐照度在建模前通常要做标幺化处理即用当前辐照度除以标准测试条件STC下的辐照度一般取1000 W/m²得到一个介于0和1之间的值。Beta分布恰好是定义在[0,1]区间上的连续分布密度函数可以呈现U形、钟形、单调递增或递减等多种形态灵活度非常高。Beta分布的密度函数为f(x) [Γ(α β) / (Γ(α) · Γ(β))] × x^(α-1) × (1-x)^(β-1)其中 α 和 β 是两个形状参数。当 α 1、β 1 时曲线呈钟形当 α β 1 时退化为[0,1]上的均匀分布当 α 1、β 1 时呈U形。这种形态变化能力对光照建模特别关键——晴天条件下辐照度普遍偏高分布集中在0.7附近多云天气则可能呈现低辐照度频发甚至出现双峰现象。Beta分布可以通过调整两个参数逼近这三种典型状态。更重要的一个理由是Beta分布对任意[0,1]上的连续随机变量都是一种良好的一阶近似而光伏出力标幺值恰好就属于这一类。相比之下对数正态分布在辐照度较高时拟合得还行但在接近0的区域往往拟合得很差而实际光伏出力在阴雨天和早晚时段的零出力占比相当大这个尾巴拟合不好整个模型的可信度就大打折扣。1.3 为什么要把两个分布组合在一起单一分布建模只是第一步实际工程里更关心的是风电和光伏同时作用下的系统总出力形态。风电的随机性和光伏的随机性在时间尺度上存在天然的互补关系——白天光伏出力高但风速可能相对平稳夜间光伏为零风电往往成为主力。把Weibull分布和Beta分布组合到一起本质上是在构建一个包含两个随机源的联合出力模型用来回答以下几个实际工程问题风光互补系统的总出力在什么区间波动概率密度形态如何系统出力低于某个阈值比如负荷需求的缺电概率有多大在给定置信水平下可被系统稳定消纳的风光出力是多少储能配置多大容量才能平滑掉一定比例的风险这类问题单独看风速分布或单独看辐照度分布都回答不了必须建立组合模型。最常用的做法是假设两个随机变量相互独立先分别拟合出各自的边缘分布再通过Monte Carlo抽样生成联合场景。独立性假设在小时级尺度上基本成立虽然日尺度上风电和光伏可能都受天气系统影响但工程上常用Copula理论做进一步修正这篇文章后面会展开讲。2. 风电Weibull分布的Matlab拟合从数据清洗到参数估计建模的第一步永远是数据清洗而不是参数拟合。我见过不少人在数据还没整理干净的情况下直接调用wblfit最后拟合出来的参数连自己都解释不了问题往往就出在数据质量上。2.1 预处理别让坏数据毁掉整个拟合风速数据最常用的来源是测风塔的10分钟平均数据或者气象站逐小时数据。拿到数据后我建议按以下顺序处理剔除超出物理范围的值比如风速大于60 m/s的记录基本是传感器故障或雷击干扰。检查零风速和负风速记录。零风速可能是静风期的正常值但如果某一天连续出现大量零值要确认是不是测风设备结冰或堵塞。处理缺失值。直接删除即可不需要插值——插值会平滑掉风速的随机波动导致拟合出的 k 参数偏大。明确建模时段。全年统一拟合和分季节拟合的结果差异很大我的建议是先画月风速曲线如果季节差异明显就按季节或按典型月份分开建模。预处理完之后把风速序列转成列向量并统计一下样本量。Weibull参数估计对样本量有一定要求实际经验是少于500个有效样本拟合结果的置信区间会非常宽做工程判断时要谨慎。2.2 三种参数估计方法与Matlab实现Matlab的Statistics and Machine Learning Toolbox里自带最大似然估计函数这是最推荐的方式。但如果你希望完全掌控计算过程或者需要深入理解参数估计的机理我建议把下面三种方法都掌握。方法一极大似然估计MLE这是最标准、统计性质最好的方法。Matlab中可以直接调用% 假设 v 是清洗后的风速列向量单位 m/s v v(~isnan(v)); v v(v 0); % 零风速不参与拟合可以根据业务需求决定 % wblfit返回 [k, c] [k_mle, c_mle] wblfit(v); % 如果想同时获取置信区间 [params, pci] wblfit(v); k_mle params(1); c_mle params(2); fprintf(Weibull MLE: k %.4f, c %.4f m/s\n, k_mle, c_mle);wblfit内部通过迭代求解对数似然函数的最大值。对数似然函数对 k 求导后是一个非线性方程Matlab默认用牛顿迭代法求解初值选得不好可能导致不收敛。如果你的数据分布很怪异可以先手动给定一个初值试试或者换下面两种方法估算一个初值。方法二最小二乘法线性化Weibull分布的累积分布函数为F(v) 1 - exp[-(v/c)^k]做两次对数变换ln[-ln(1 - F)] k·ln(v) - k·ln(c)这是一个典型的直线方程斜率为 k截距为 -k·ln(c)。实际操作中先对风速排序用经验累积频率公式 F_i (i - 0.3) / (n 0.4) 计算每个点的累积概率然后做线性回归。vs sort(v); n length(vs); F ((1:n) - 0.3) / (n 0.4); x log(vs(:)); y log(-log(1 - F)); % 一次线性回归 p polyfit(x, y, 1); k_lsq p(1); c_lsq exp(-p(2) / k_lsq);这段代码的好处是直观、可解释性强而且完全不需要任何工具箱。缺点是经验累积频率在两端会失真尾部的拟合精度较低。实际使用中我会把线性回归的结果作为MLE的初值两者结合效果最好。方法三矩估计法Weibull分布的均值 μ 和标准差 σ 与参数的关系是μ c·Γ(1 1/k)σ² c²·[Γ(1 2/k) - Γ²(1 1/k)]通过数值求解这两个方程即可得到 k 和 c。Matlab实现需要fsolvemu mean(v); sd std(v); mom_egn (x) [c_est(x, mu, sd) - mu; sd_egn(x, mu, sd) - sd]; % 这里直接给出数值求解更简洁的形式 obj (k) (gamma(1 2/k) ./ gamma(1 1/k).^2 - 1) - (sd / mu)^2; k_mom fzero(obj, 2); c_mom mu / gamma(1 1/k_mom);矩估计法计算速度快适合数据量大、需要批量处理的场景但效率不如MLE。三种方法的对比如下方法精度计算量适用场景MLE最高中等一般工程分析首选最小二乘线性化中等低快速估算、MLE初值矩估计中等低批量处理、实时计算2.3 拟合优度验证直方图、Q-Q图与K-S检验拟合完成后不能直接拿去用必须验证。我的标准验证流程是三步。第一步把概率密度曲线叠加到频率直方图上用肉眼检查形态是否吻合。注意直方图要用频率密度频数/bin宽度/总样本数而不是频数否则和概率密度函数不在一个量纲上看起来会差很多。figure; histogram(v, 30, Normalization, pdf, FaceColor, [0.7 0.7 0.7]); hold on; v_axis 0:0.1:max(v)*1.1; pdf_weibull wblpdf(v_axis, k_mle, c_mle); plot(v_axis, pdf_weibull, r-, LineWidth, 1.5); xlabel(风速 (m/s)); ylabel(概率密度); legend(实测频数密度, Weibull拟合);第二步画Q-Q图。figure; qqplot(v, weibull_dist);如果散点基本落在直线上说明分位数对齐良好尤其是尾部。Q-Q图比直方图更能反映尾部拟合质量这对评估高风速段的极端出力特别重要。第三步做K-S检验。[h, p] kstest(v, CDF, wblcdf(v, k_mle, c_mle));p值大于0.05意味着在95%置信水平下不能拒绝原假设即数据与Weibull分布一致。但我必须提醒一点K-S检验对样本量非常敏感几万条样本下几乎任何分布都会被拒绝所以当样本量很大时别盯着p值看而是看Q-Q图的偏离程度和实际工程误差。3. 光伏Beta分布的Matlab实现标幺化与参数求解光伏建模的核心变量是辐照度但在实际工程里大家更关心出力。Beta分布直接建模的对象是辐照度标幺值或者功率标幺值。这里有一个选择到底对哪个变量建模更合理我的经验是如果做的是长期规划分析直接对辐照度标幺值建模更稳定如果做的是并网运行分析对实际出力标幺值建模更方便因为已经包含了逆变器效率和限功率因素的影响。3.1 辐照度到功率的标幺化转换光伏出力与辐照度之间的简化关系可以写成P_pv P_rated × (R / R_stc) × [1 γ × (T_cell - T_stc)]其中 R 是水平面辐照度或组件面辐照度R_stc 1000 W/m²P_rated 是额定容量γ 是温度系数T_cell 是电池板工作温度。如果暂时忽略温度影响P_rated 的标幺值出力约等于辐照度标幺值。这份工作里我们主要关注辐照度分布的统计特性所以直接对 x R / R_stc 建模。需要注意的是实测辐照度数据里会掺杂大量夜间零值。如果直接把全天24小时的数据拿去拟合Beta分布零值的比例可能高达40%到50%这会导致 Beta 分布参数严重失真。我在实际操作中会做分时段处理只取白天有效日照时段比如早上6点到晚上19点之间或者按辐照度阈值 10 W/m² 筛选建模半夜的零值放到时序场景生成时再处理。这样才符合物理逻辑。3.2 矩估计与极大似然估计的实现Beta分布参数估计在Matlab里没有像wblfit一样的一行函数但可以通过betafit直接完成MLE估计也可以用矩估计手动求解。方法一极大似然估计% x 是辐照度标幺值取值在(0,1)之间 x x(x 0 x 1); % 去掉0和1避免对数计算发散 % betafit返回 [alpha, beta] params_beta betafit(x); alpha params_beta(1); beta params_beta(2); fprintf(Beta MLE: alpha %.4f, beta %.4f\n, alpha, beta);betafit内部也是迭代计算如果数据里包含0或1会直接报错或者给出无穷参数所以我习惯先做一次截断处理。方法二矩估计矩估计的好处是不依赖工具箱而且可以写成非常简洁的解析公式。Beta分布的均值和方差为μ α / (α β)σ² αβ / [(α β)²(α β 1)]反解得α μ × (μ(1-μ)/σ² - 1)β (1-μ) × (μ(1-μ)/σ² - 1)Matlab代码mu mean(x); var_x var(x); % 注意var_x 必须小于 mu*(1-mu)否则参数为负说明数据不适合Beta分布 denom mu * (1 - mu) / var_x - 1; alpha_mom mu * denom; beta_mom (1 - mu) * denom; if alpha_mom 0 || beta_mom 0 warning(矩估计得到的参数非正数据可能不适合Beta分布); end实测经验表明在样本量较大的情况下矩估计和MLE的结果差异很小。但如果辐照度数据比较干净且模型形态很标准优先用MLE如果只是想快速拿到一套参数做初步仿真矩估计足够了。3.3 边界样本与参数失真的处理光伏数据里最棘手的问题就是大量0和1的出现。辐照度标幺值为0的时间段占比很大而标幺值为1甚至大于1的情况也会出现——比如某些高海拔地区或者聚光光伏系统辐照度可能瞬间超过1000 W/m²。Beta分布的支持域是开区间(0,1)端点处概率密度为0当参数大于1时。直接把0、1扔进去拟合会有两个问题一是对数似然函数在端点处发散计算不稳定二是拟合结果会严重偏向端点导致中间段的形态贴合度很差。我的处理策略分三层第一层数据筛选。按业务逻辑决定0值怎么处理。如果是辐照度建模且研究重点是白天时段的光伏出力那么直接剔除0值模型只对有效日照时段的随机性负责。第二层端点微调。如果数据里包含少量1或接近1的值可以用一个极小的偏移量把它们拉回(0,1)区间比如x min(1 - 1e-6, max(1e-6, x))。这个做法对参数估计的影响可以忽略但能避免数值发散。第三层零出力注入模型。如果后续场景生成时需要模拟全天24小时出力不要在Beta分布内部塞0而是用混合模型以一定概率 p0 产生零出力以剩余概率从Beta分布采样。这样一个两步抽样逻辑既保持了Beta分布的形态又真实反映了光伏夜间和阴雨天的特性。% 夜间或低辐照时段比例 p0 sum(R 10) / length(R); % 辐照度阈值自行调整 % 有效时段Beta采样 x_eff betarnd(alpha, beta, N, 1); % 混合抽样 u rand(N, 1); x_sample zeros(N, 1); x_sample(u p0) x_eff(u p0);这种处理办法在后续做全年8760小时时序仿真时非常管用推荐直接用。4. 组合建模与场景生成从概率密度到实际出力序列边缘分布拟合完成后真正的重点来了——两个分布怎么组合组合之后怎么生成可以用于仿真计算的出力场景。这一节是整个项目最核心的工程环节。4.1 独立假设下的联合概率模型如果假设风速和辐照度相互独立那么联合概率密度就是两个边缘密度的乘积。这个假设在小时级时间尺度上有一定的合理性——风主要受大尺度气压系统驱动光照主要受云和大气透明度影响两者之间没有直接的物理耦合。但要注意在日尺度或更长尺度上天气系统可能同时影响风和光比如锋面过境时风速增大、云量增多此时独立假设会带来偏差。想要充分考虑相关性可以用Copula函数建模其中常用的有Gaussian Copula和t-Copula。这里先给出独立假设的实现相关性建模作为扩展放在后面。联合模型的意义在于它把一个双随机源的工程问题转化成了一台可以抽样的随机数生成器。有了这个生成器我们就可以做Monte Carlo模拟把风速、辐照度样本对同时生成出来再通过各自的功率转换曲线得到总出力。4.2 Monte Carlo场景生成完整代码下面给出一个完整的场景生成函数。输入是Weibull参数、Beta参数和系统容量参数输出是场景矩阵。function [P_wind, P_pv, P_total] generate_hybrid_scenarios(k_w, c_w, alpha_pv, beta_pv, ... P_wind_rated, P_pv_rated, N_scenarios, v_cut_in, v_cut_out, v_rated, pv_p0) % 输入参数说明 % k_w, c_w: 风速Weibull参数 % alpha_pv, beta_pv: 辐照度Beta参数 % P_wind_rated: 风电场额定功率 % P_pv_rated: 光伏电站额定功率 % N_scenarios: 场景数量 % v_cut_in, v_cut_out, v_rated: 风机切入/切出/额定风速 % pv_p0: 光伏零出力概率 % 1. 采样风速 v_sample wblrnd(k_w, c_w, N_scenarios, 1); % 2. 采样辐照度标幺值零出力混合 u rand(N_scenarios, 1); x_sample zeros(N_scenarios, 1); x_sample(u pv_p0) betarnd(alpha_pv, beta_pv, sum(u pv_p0), 1); % 3. 风速转风电出力 P_wind zeros(N_scenarios, 1); idx_region1 (v_sample v_cut_in) (v_sample v_rated); idx_region2 (v_sample v_rated) (v_sample v_cut_out); P_wind(idx_region1) P_wind_rated .* ... (v_sample(idx_region1).^3 - v_cut_in^3) ./ (v_rated^3 - v_cut_in^3); P_wind(idx_region2) P_wind_rated; % 4. 辐照度标幺值转光伏出力线性近似 P_pv P_pv_rated .* x_sample; % 5. 总出力 P_total P_wind P_pv; end关于风机功率曲线我用的是标准的线性三次方段模型即切入风速以下不出力额定风速到切出风速之间保持额定出力切出风速以上停机。实际工程中更精确的做法是查厂家提供的功率曲线表但做统计研究时这个模型已经足够。生成场景之后可以画一个散点图横轴是风电出力、纵轴是光伏出力颜色代表总出力的高低。从这个散点图可以直观看到联合出力的分布范围以及风光出力互补的状态在哪个区间分布最密集。scatter(P_wind / P_wind_rated, P_pv / P_pv_rated, 10, P_total / (P_wind_rated P_pv_rated), filled); colorbar; xlabel(风电出力标幺值); ylabel(光伏出力标幺值); title(风光联合场景分布);4.3 场景削减与典型场景选取Monte Carlo生成的场景数量如果只有几百个还好但如果要做一年8760小时的时序模拟每个时刻都生成大量场景计算规模会迅速膨胀。这时候需要场景削减——从大量场景里挑出一组具有代表性的“典型场景”用少量场景近似逼近原始联合分布。最常用的削减方法是K-means聚类。把每个场景看成二维空间里的一个点K-means把整个点集划分为K类每个类的质心就作为代表性场景各类的样本占比作为该场景的概率权重。% 假设场景矩阵 X 是 N行2列 [风电, 光伏] X [P_wind, P_pv]; % K-means聚类K根据实际需要选择一般10~50 k_idx 30; [cluster_idx, C] kmeans(X, k_idx, Replicates, 5, MaxIter, 500); % 计算每个类的比例作为权重 weights histcounts(cluster_idx, k_idx) / length(cluster_idx); % 典型场景为质心坐标 typical_scenarios C;削减后的场景集可以直接用于后续的随机潮流计算、储能容量优化或可靠性分析。实测下来30个典型场景就能在绝大部分情况下保持与原分布的前四阶矩基本一致计算量却降低了一个数量级。场景削减有个经验值先用肘部法则选K。画一下不同K对应的组内离差平方和随着K增加下降曲线出现明显拐点的位置就是合适的K。我一般不追求特别大的K因为场景削减的目的是保留“分布特征”而不是保留“所有细节”。5. 实操中反复踩过的坑与工程化建议最后这部分写点干货中的干货全是实际项目里踩过坑之后总结出来的。如果你按照前面的代码直接跑大概率能出结果但要做出可信的工程结论下面这几个问题必须提前处理。5.1 拟合误差不等于模型误差做拟合的人容易陷入一个误区追求拟合优度越高越好用各种指标轰炸。但工程建模的核心目标不是让曲线穿过每一个点而是让模型能正确回答业务问题。比如你做储能容量配置最关心的是出力低于负荷的概率和持续时长分布这时即使Weibull拟合的K-S检验p值低一些只要概率密度曲线在低出力区间的累计概率准确结果就是可用的。反过来如果K-S检验通过了但Q-Q图的尾部偏差明显说明模型对极端低风速或极端高辐照度场景描述不准确。这种情况下的概率密度“看起来”没问题但做可靠性评估时会严重低估极端场景的发生概率。我建议每次拟合完都同时看三个指标K-S检验p值、Q-Q图形态、以及你关心的业务区间的累计概率误差。三者结合才能判断模型是否真的可用。5.2 季节与天气分型比一味追求高精度更重要一个全年统一拟合的Weibull/Beta模型在统计上是“正确”的但在工程上往往是粗糙的。风资源有显著的季节差异比如我做过的一个西北地区的项目春季大风天的平均风速能比秋季高40%对应的k值也从2.1变化到3.2光伏出力更是跟季节日照时长和太阳高度角强相关。如果全年只用一个模型生成的场景会淡化季节差异导致储能配置结果偏乐观或偏悲观。我的做法是先按月或按季节聚类分时段拟合。更细致一点的做法是按天气类型分型从历史数据里识别出晴朗、多云、阴雨、大风四种典型天气分别对每个天气类型拟合Wind-PV联合分布然后用天气类型的出现概率做加权。这样建出来的模型能够生成更真实的“晴天但大风”、“阴天但微风”这类场景组合对系统风险评估的增益非常明显。5.3 代码效率向量化与并行工具箱如果你要生成几百万个场景做Monte Carlo模拟循环写法会非常慢。Matlab的随机数生成函数本身是向量化的但我见过不少人在场景转换环节写for循环比如对每个场景逐一计算风机功率。这里建议直接用逻辑索引批量计算性能差距能达到几十倍。% 低效做法 for i 1:N if v(i) v_rated p(i) P_rated; end end % 高效做法 p zeros(N, 1); p(v v_rated v v_cut_out) P_rated;如果再往后做多风电场、多光伏电站的联合建模单机串行计算会非常吃力。Matlab的Parallel Computing Toolbox里parfor是现成的场景之间如果完全独立直接循环并行即可。我最初做某个省级电网的风光出力时序仿真时用parfor把8760小时的场景生成从半小时压缩到了两三分钟这个优化非常值。5.4 参数估计对样本量的敏感性最后强调一个容易被忽略的问题样本量。比如你只有两个月的测风数据却想拟合出代表全年特性的Weibull参数结果必然偏颇。通常我要求至少有一整年的有效数据且时间分辨率达到小时级。如果数据不足宁可退回到用参考气象站的多年平均风速来标定尺度参数c也不要用短样本硬拟合。Beta分布的样本量要求同样很高而且辐照度数据的有效样本量往往比风速还要少因为夜间数据几乎全被筛掉了。假设你有一年小时级数据晚上6点到早上6点剔除后有效样本可能只有3000个左右勉强够用但置信区间仍然偏宽。这时候我建议把多年数据合并使用或者采用分季节建模来增加各季节的有效样本量。结尾的一点体会这套WeibullBeta的组合建模方法我前后用了小半年才稳定下来最大的体会是模型本身不复杂复杂的是数据处理的边界条件——什么时候该剔除零值、什么时候该分季节拟合、什么时候该用混合分布这些问题只能通过反复对比实测数据来形成判断。最后再分享一个小技巧每次拟合完把拟合参数和对应的K-S检验结果存成一个结构体变量方便不同方案之间对比。这个过程比跑通代码更有价值因为它逼着你去理解数据而不是仅仅让代码跑出数字。这套方法后续还可以往Copula相关性建模、动态参数时变模型方向扩展等我把相关性部分的代码也整理完了再单独写一篇分享。
