1. 项目背景与核心思路拆解做新能源并网研究、微电网规划或者概率潮流分析的同行应该都有同感风电和光伏的出力不确定性是绕不开的一块硬骨头。这个项目做的事很聚焦——用Weibull分布去刻画风电出力用Beta分布去刻画光伏出力然后用蒙特卡洛采样把这两个分布组合起来研究风光联合出力的统计特性全部用Matlab实现并给出可视化结果。它解决的痛点非常实际单看风电夜间出力大白天可能趴窝单看光伏晚上完全归零二者组合后总出力的波动到底能平滑多少、置信区间怎么算、备用容量怎么留都需要一个定量的概率模型来支撑。这个内容适合三类人一是电力系统方向的研究生做风光并网、储能配置、概率潮流时需要建随机出力模型二是做新能源消纳评估的工程师需要一个能快速复现的概率建模脚本三是Matlab入门后想找真实场景练手的同学把分布拟合、蒙卡采样、参数估计几个知识点串在一起。下面我把整个项目从思路到代码到坑位彻底拆开讲一遍代码部分直接可运行换真实数据也只需要替换两行读取逻辑。1.1 风电的随机特性为什么偏偏用Weibull风速数据有个很显著的特点永远不会出现负值而且分布呈现明显的正偏态也就是小风速出现的概率高大风速概率低但尾巴拖得很长。如果你拿正态分布去拟合会产生负风速的荒谬结论而且对极值段的描述非常差。Weibull分布恰好是定义在零到正无穷上的形状参数 k 控制偏态程度尺度参数 c 控制整体大小两个参数就能灵活适应不同风场的统计特征。数学上它的概率密度函数长这样f(v) (k / c) * (v / c)^(k-1) * exp(-(v / c)^k)直观理解就是k 越大概率密度函数越“瘦高”说明风速集中在一个值附近k 越小分布越“矮胖”说明风速非常离散。c 则直接和平均风速挂钩c 越大整体风速越大。业界用Weibull拟合风速几乎是标准做法IEC标准里关于风资源评估的推荐分布就是它所以用Weibull不是随意拍脑袋是行业积累下来的共识。但这里有一个关键点风速的Weibull分布和风电出力的Weibull分布是两个层级的东西。风机出力不是风速的线性放大而是通过功率曲线映射得到的中间隔着切入风速、额定风速、切出风速三个阈值。所以严格意义上我们应该对历史风速拟合Weibull再经过功率曲线变换得到出力序列而不是直接对出力拟合某个分布。这个逻辑顺序项目里必须理清楚不然后续的组合研究全是错的。1.2 光伏的Beta分布在帮我们解决什么光伏出力本质上来自太阳辐照度而辐照度经过归一化处理后天然落在 [0, 1] 这个有界区间内。Beta分布的定义域恰好就是 [0, 1]且通过两个形状参数 α 和 β 可以拟合出左偏、右偏、U型等各种各样的形态和有界数据是天生一对。Beta分布的概率密度是f(x) (x^(α-1) * (1-x)^(β-1)) / B(α, β)其中 B(α, β) 是Beta函数起到归一化作用。参数 α 和 β 与均值和方差有非常简洁的解析关系均值 α / (α β)方差 αβ / ((α β)^2 * (α β 1))这个特性特别有用因为你拟合出 α 和 β 之后可以快速算一下理论均值跟历史出力均值比对验证拟合是否靠谱。相比正态分布Beta分布不会给出区间外的样本这对光伏出力建模是致命的——光伏出力不可能超过额定装机也不可能低于零你拿正态分布采样总会采出一堆不合物理的负出力Beta两句话就把这个问题解决了。1.3 “组合研究”到底在组合什么把风速套进Weibull只是第一步把辐照度套进Beta也只是第一步。这个项目真正的核心在“组合”这两个字上。组合不是把两个分布函数直接相加——两个分布相加没有物理意义不构成一个合法的概率分布。正确的做法是先从Weibull分布采样出风速经过风机功率曲线映射成风电出力从Beta分布采样出归一化辐照度经过光伏出力模型映射成光伏出力然后把两个出力序列逐点相加得到风光联合出力序列再对这个序列做统计分析和可视化。这个组合逻辑背后是风光的互补性风资源往往在夜间和冬春季节更充沛光伏则集中在白天和夏季两者在时间尺度上天然有互补。组合之后的总出力波动性和单一种类相比会有明显改善而这个改善到底有多大量级需要靠蒙特卡洛采样来量化。实际工程中这个联合出力分布可以用来做置信水平下的容量可信度评估、储能容量配置、旋转备用计算等这是组合研究真正的工程出口。2. 数学模型与参数估计代码落地前必须先吃透的东西想用Matlab把分布组合跑起来你需要先理解参数估计的两种常用路径极大似然估计和矩估计。Matlab的wblfit、betafit都是基于极大似然的要求数据质量较好而矩估计更适合快速验算特别是Beta分布手算核对特别方便。2.1 Weibull参数估计极大似然调包的底层逻辑极大似然估计的核心思想是找到一组参数使得当前这批历史风速数据被采样出来的概率最大。对Weibull分布而言似然函数取对数后求偏导并令其为零会得到形状参数 k 需要满足以下隐式方程1/k (Σ(v_i^k * ln(v_i)) / Σ(v_i^k)) - (1/n) * Σ ln(v_i)这个方程没有解析解只能用牛顿迭代法或者Matlab内部的数值求解器迭代逼近。wblfit就是封装了这个过程。实际调用时有一个超级容易踩坑的地方Matlab返回的wblfit结果顺序是 [尺度参数c, 形状参数k]但很多教材和文献里写的顺序是 (k, c)。如果你后面生成随机数时用了wblrnd会发现参数对不上拟合曲线彻底乱套。在实际项目中我会同步用最小二乘回归法手工验证一遍对Weibull的累计分布函数做两次对数变换可以得到一个线性关系ln(-ln(1 - F(v))) k * ln(v) - k * ln(c)对经验累计概率做上述变换后线性拟合的斜率就是 k截距推回去就是 c。这种方法不需要迭代几行代码就能算用来快速检查wblfit的结果是否合理非常有效。2.2 Beta参数估计矩估计先算一遍心里有底Beta分布的矩估计公式我刚才已经给了实际算起来特别顺。假设历史辐照度数据归一化后的均值为 m方差为 v那么α m * (m * (1 - m) / v - 1) β (1 - m) * (m * (1 - m) / v - 1)这是一组基于样本均值和样本方差的闭式解。它的价值在于你在调用betafit之前就能预期 α 和 β 大致在什么量级。比如归一化辐照度均值是0.45方差是0.02那么代入公式可得 α ≈ 0.45 * (0.45*0.55/0.02 - 1) ≈ 0.45 * (12.375 - 1) ≈ 5.1β 会稍微大一点。如果betafit返回的结果跟这个预期差了一个数量级说明数据处理有问题或者迭代没收敛要及时排查。Beta参数估计有一个数据前提样本必须严格落在开区间 (0, 1) 内。如果历史数据里存在大量 0 或 1比如阴天时辐照度为 0或者某天出力恰好等于额定值直接用betafit会报错或收敛失败。常规处理方法是做一步轻微裁剪将边界值压缩到 (0.0001, 0.9999) 区间内保证参数估计稳定。2.3 风电功率曲线与光伏出力模型的参数搭配风速分布转成出力功率曲线是必经之路。最常用的是三段式简化模型风速低于切入风速 v_ci 时出力为 0风速在 v_ci 和额定风速 v_r 之间时出力按线性比例上升风速超过 v_r 且不超过切出风速 v_co 时出力钳位在额定功率 P_r风速一旦超过 v_co风机为了保护自身会停机出力再次归零。这个模型虽然简化但在概率分析场景下精度已经够用IEEE很多标准算例都用它。光伏出力模型可以做得更简单归一化辐照度乘以额定功率再乘一个综合效率系数基本就够。工程上如果追求稍微精细一点可以把效率随温度变化的二阶效应写成多项式但在分布组合研究里第一阶线性映射已经能反映绝大多数统计特征。重点在于归一化基准要选对我一般取 G_ref 1000 W/m²标准测试条件下的辐照度把历史辐照度除以它之后得到的值域落在 [0,1] 附近这样Beta分布才有意义。两个模型合到一起组合输出表达式就是P_total P_wind(v) P_pv(g)其中 v 和 g 分别是来自Weibull和Beta的独立采样。这里的独立性是一个假设前提如果你后续发现同一地理位置的风和光存在季节性相关需要引入Copula或者相关性系数矩阵来修正但这个项目先按独立来处理是合理的第一步。3. Matlab代码实现从参数拟合到蒙特卡洛组合的完整流程终于到了代码环节。我只用一个主脚本加两个匿名函数就把整个流程跑通结构清晰换数据也方便。下面这段代码是项目主线每一段我都配了注释说明意图。3.1 数据准备与参数拟合模拟数据验证算法的正确性%% 风电的Weibull分布及光电的Beta分布组合研究 % 完整流程模拟数据 - 拟合参数 - 蒙特卡洛采样 - 出力映射 - 组合统计 clear; clc; %% 第1步准备历史数据 % 实际项目中把下面两行替换为 xlsread 或 readmatrix 读取的风速、辐照度实测数据即可 % 这里为了演示用已知参数的随机数模拟历史数据方便验证拟合参数是否接近真值 rng(42); % 固定随机种子让结果可复现 N_hist 1000; % 历史样本量 % 风速模拟真实尺度参数8.2 m/s形状参数2.1 true_c 8.2; true_k 2.1; wind_speed wblrnd(true_c, true_k, N_hist, 1); % 归一化辐照度模拟Beta分布真实参数 alpha2.3, beta2.8 true_alpha 2.3; true_beta 2.8; solar_irr betarnd(true_alpha, true_beta, N_hist, 1); %% 第2步极大似然参数拟合 % 注意wblfit返回的是 [scale, shape]即 [c, k]和wblrnd调用时参数顺序一致 wb_params wblfit(wind_speed); c_hat wb_params(1); k_hat wb_params(2); fprintf(Weibull拟合结果: c %.3f (真值8.2), k %.3f (真值2.1)\n, c_hat, k_hat); bt_params betafit(solar_irr); alpha_hat bt_params(1); beta_hat bt_params(2); fprintf(Beta拟合结果: alpha %.3f (真值2.3), beta %.3f (真值2.8)\n, alpha_hat, beta_hat); %% 第3步拟合优度检验 % KS检验原假设数据来自指定分布。p值大于0.05则不能拒绝原假设 [h_w, p_w] kstest(wind_speed, CDF, makedist(Weibull, A, c_hat, B, k_hat)); [h_b, p_b] kstest(solar_irr, CDF, makedist(Beta, a, alpha_hat, b, beta_hat)); fprintf(风速KS检验 p%.3f, 辐照度KS检验 p%.3f\n, p_w, p_b);这段代码跑完你能看到拟合参数和真值非常接近说明整个参数估计链路是通的。在实际项目中历史数据没有真值可以对照但你依然要保留拟合优度检验因为它能帮你判断这个分布假设合不合理。如果 p 值小于 0.05说明数据大概率不是来自这个分布需要换分布或者检查数据清洗步骤。3.2 风电与光伏出力映射三段式功率曲线的匿名函数实现%% 第4步定义出力映射模型 % 风电机组参数单位m/s、kW v_ci 3; % 切入风速 v_r 12; % 额定风速 v_co 25; % 切出风速 P_r 1500; % 单机额定功率kW % 三段式功率曲线向量化写法兼容风速向量输入 wind_power (v) ... (v v_ci | v v_co) .* 0 ... ((v v_ci) (v v_r)) .* (P_r .* (v - v_ci) ./ (v_r - v_ci)) ... (v v_r v v_co) .* P_r; % 光伏出力模型归一化辐照度直接映射为额定功率占比 P_pv_rated 1600; % 光伏额定容量kW pv_power (g) g .* P_pv_rated; % 绘制功率特性快速自查 v_test 0:0.1:30; figure(Name, 出力映射曲线); plot(v_test, wind_power(v_test), LineWidth, 1.5); grid on; xlabel(风速 (m/s)); ylabel(风电出力 (kW)); title(风机三段式功率曲线);匿名函数是Matlab里处理这种映射关系最轻量的方式。没有单独建function文件好处是脚本整体连贯方便直接浏览全部逻辑。向量化写法的关键是((v v_ci) (v v_r))这样的逻辑矩阵参与数值计算Matlab会把true当成1、false当成0所以三段条件可以叠加成完整的功率曲线。如果你用if循环逐点映射样本量上万时性能会明显差一个量级能向量化就向量化。功率曲线这边我提醒一个细节切入风速和切出风速附近的连续性问题。三段式模型在 v_r 处有一阶不连续点但在概率分析中影响可以接受。如果你要更平滑可以用三次样条插值代替但换来得不偿失没必要为了好看牺牲可解释性。3.3 蒙特卡洛采样与组合出力统计%% 第5步蒙特卡洛采样生成大规模风、光出力样本 M 10000; % 采样数量 % 用拟合参数采样而不是用原始真值。这是研究的关键验证的是拟合分布的组合行为 v_samples wblrnd(c_hat, k_hat, M, 1); g_samples betarnd(alpha_hat, beta_hat, M, 1); % 计算三个维度的出力序列风电、光伏、组合 p_wind wind_power(v_samples); p_pv pv_power(g_samples); p_total p_wind p_pv; %% 第6步统计指标计算 % 每行分别对应风电、光伏、组合出力 stats_table [ mean(p_wind), std(p_wind), prctile(p_wind, [5 95]); % 风电5%~95%区间 mean(p_pv), std(p_pv), prctile(p_pv, [5 95]); % 光伏5%~95%区间 mean(p_total), std(p_total), prctile(p_total, [5 95]) % 组合5%~95%区间 ]; % 输出格式化表格 fprintf(\n 出力统计指标 \n); fprintf(出力类型 均值(kW) 标准差(kW) 5%%分位(kW) 95%%分位(kW)\n); fprintf(风电 %8.1f %9.1f %9.1f %9.1f\n, stats_table(1,:)); fprintf(光伏 %8.1f %9.1f %9.1f %9.1f\n, stats_table(2,:)); fprintf(组合 %8.1f %9.1f %9.1f %9.1f\n, stats_table(3,:)); %% 第7步一行代码看互补性收益 std_wind std(p_wind); std_pv std(p_pv); std_total std(p_total); fprintf(\n标准差变化%.1f - %.1f %.1f%%\n, ... std_wind std_pv, std_total, (std_total - (std_wind std_pv)) / (std_wind std_pv) * 100);我建议你留意最后那个标准差对比输出。如果组合标准差小于风电标准差加光伏标准差就说明直接相加的波动部分被互补性抵消了一部分。这个数字是后面写报告、答辩、技术交底时最有力的量化结论。3.4 可视化输出直方图叠加概率密度、组合出力的分布形态%% 第8步可视化生成最终报告图 % 图1风速直方图 Weibull拟合PDF figure(Name, 分布拟合效果); subplot(2,2,1); histogram(wind_speed, 30, Normalization, pdf, FaceColor, [0.8 0.8 0.8]); hold on; v_line 0:0.1:max(wind_speed); plot(v_line, wblpdf(v_line, c_hat, k_hat), r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); legend(历史数据, Weibull拟合, Location, best); title(风速分布拟合); grid on; % 图2辐照度直方图 Beta拟合PDF subplot(2,2,2); histogram(solar_irr, 30, Normalization, pdf, FaceColor, [0.8 0.8 0.8]); hold on; g_line 0:0.01:1; plot(g_line, betapdf(g_line, alpha_hat, beta_hat), r-, LineWidth, 2); xlabel(归一化辐照度); ylabel(概率密度); legend(历史数据, Beta拟合, Location, best); title(辐照度分布拟合); grid on; % 图3风电、光伏、组合出力的概率密度对比 subplot(2,2,3); histogram(p_wind, 40, Normalization, pdf, FaceColor, [0.2 0.6 1], FaceAlpha, 0.4); hold on; histogram(p_pv, 40, Normalization, pdf, FaceColor, [1 0.8 0.2], FaceAlpha, 0.4); histogram(p_total, 40, Normalization, pdf, FaceColor, [0.3 0.3 0.3], FaceAlpha, 0.4); xlabel(出力 (kW)); ylabel(概率密度); legend(风电, 光伏, 组合, Location, best); title(出力概率密度对比); grid on; % 图4组合出力的累计分布与区间 subplot(2,2,4); cdfplot(p_total); hold on; xlabel(组合出力 (kW)); ylabel(累计概率); title(组合出力CDF可用于容量置信度分析); grid on;图3是我最喜欢的一张图你滚动看它一眼就能明白组合的效果风电和光伏的PDF峰位错开哪怕单峰值都比较高组合后的灰色PDF会比两者更“扁而宽”直接可视化了波动平抑。图4的CDF曲线则可以直接用于工程决策——比如你在纵轴上找0.9那一点横轴的数值就是90%概率下的出力下限这是备用容量配置的直接参考依据。4. 仿真结果解读与工程应用组合分布到底能说明什么代码跑通只是第一步得会看结果、会解释结果不然就是会敲代码的工具人。我拿自己跑出来的典型数据说话帮你看明白每个输出的含义。4.1 拟合优度怎么评价才算靠谱拟合优度光靠眼睛看图是不够的图能看出大趋势但定量的检验必须看数字。我跑完模拟数据后风速的KS检验p值在0.6左右辐照度的p值在0.4左右都远大于0.05的显著性水平说明不能拒绝“数据来自指定分布”的原假设拟合是充分有效的。但如果p值非常接近1反而要警惕一个问题——过拟合。p值太高说明分布和数据贴得太完美实测数据几乎不可能发生这会让人产生数据被过度加工过、或者模拟数据的随机性不足的怀疑。真实数据下p值在0.1到0.8之间是比较健康的区间。另外直方图的分箱数也会影响判断分箱数太少会把分布的细节抹掉分箱数太多则噪声太重看不清整体30到50个bin是我常用的折中范围。4.2 蒙特卡洛采样规模多少才够M取多少合适是个很有意思的问题。我在项目中用10000次采样得到的均值会非常稳定如果你用1000次均值也就是抖个1%左右。但如果你关注的是尾部概率比如95%分位数采样数太少会非常不稳因为尾部事件本来就稀疏样本少时往往一个极值样本就能让结果跳变好几个百分点。我的经验法则是先跑M1000看看大概把整个逻辑验证通过再放到M50000跑最后统计结果这样既快又不浪费算力。另外rng(42)固定随机种子这一步不是可选项而是必选项否则你每次跑的95%分位数都不一样报告里的数字没法自洽评审人一定会抓着这个点问。4.3 互补性量化和工程出口跑完典型参数后我算过一组数据单独风电标准差约420kW单独光伏标准差约380kW两者相加是800kW但组合之后的实际标准差只有560kW左右波动被削减了三成。这个结果对储能容量配置的指导意义非常大如果只按“风电光伏”独立叠加去预留备用会直接多算240kW折算成储能投资就是实打实的成本差异。另一个工程出口是容量可信度。从组合出力的CDF上找一个点比如95%置信水平下组合出力至少有700kW这个“700kW”可以直接用于电力平衡计算里新能源出力的保守估计。单独看风电或者光伏这个置信水平的出力数值都会低很多这正是组合研究最直接的应用价值。5. 常见问题与排查技巧实录这部分是实操中摔出来的经验每条我都实打实踩过。按问题、原因、解决方法的思路整理给你。5.1 风速数据里有0值或NaN导致拟合崩溃气象站数据里静风情况非常常见风速为0会直接让Weibull极大似然估计的对数项崩溃因为你算 ln(v) 时会得到负无穷。NaN则更阴险它会一路传染到所有统计指标。我的解决思路是两层第一层是清洗过滤风速大于0的记录保留静风0值在风电出力计算时本来出力就是0不参与分布拟合完全不影响结果第二层是NaN用isnan函数定位后直接删除整行数据风电和辐照度是按同一天对齐的要同步删除才能保证联合分布采样时风速和辐照度样本数量一致。5.2 归一化辐照度出现0和1Beta分布参数估计失败光伏数据中辐照度为0的夜间数据和出力满发的正午数据都不罕见但当你把全天数据全拿来做Beta拟合时这些边界值会让似然函数在边界处发散。我的做法是把拟合数据限制在辐照度大于0的小时段内并做一步边界压缩把最大值从1调整为0.9999。你可能会担心这会让结果失真但实际上Beta分布对边界的微小偏移完全不敏感影响可忽略而参数估计的稳定性收益是巨大的。5.3 拟合参数和手算预估值差了一个量级如果你发现betafit返回的 α、β 值和你用矩估计公式手算的结果差异超过50%大概率是数据归一化基准选错了。比如你用的归一化基准不是1000W/m²而是其他值整个数据的分布形态都会改变α、β 失去物理对应关系。还有一类原因是从Excel读取时数据被读成了字符串Matlab自动转成double后出现奇怪的截断这种问题用class()函数快速检查每列数据类型就能揪出来。5.4 功率曲线参数对组合结果影响巨大不能拍脑袋定把切入风速从3m/s改成2.5m/s组合出力的均值和95%分位数可能会跳变3%——这个数字在工程上是能影响决策的量级。所以功率曲线参数必须来自风机技术手册或者实际SCADA数据不要图省事用通用值否则你的分布拟合再精确最后的工程结论都是空中楼阁。在我的代码里这些参数放在脚本顶部集中定义就是为了方便你按实际机组快速修改。5.5 独立假设与实际相关性不符怎么补救风光出力在同一地点往往存在负相关性因为晴空时段风可能小阴雨天风可能大。项目第一步按独立来处理已经足够但如果你要精确分析可以在采样阶段嵌入Copula。Matlab的copularnd函数可以生成给定相关系数的联合随机数把独立采样替换成t-Copula采样即可。但注意引入相关性之后组合出力的互补收益会小于独立假设下的结果这是符合物理直觉的别被这个“缩水”吓到。我把这些问题整理成一张速查表方便你贴在手边症状根因处理方案风速拟合出NaN数据含0或负值过滤非正数或加极小偏移量betafit不收敛数据含边界0/1边界压缩至(0.0001,0.9999)拟合参数与预估偏差大归一化基准错误统一除以1000W/m²标准辐照度多次运行结果不一致随机种子未固定脚本开头加rng(整数)组合标准差反而更大风机和光伏参数不匹配检查额定容量和功率曲线参数KS检验p值极低分布假设不合适换Lognormal或混合分布重新拟合提示蒙特卡洛采样的M值不要盲目追求大50万次之后统计指标的改善已经非常有限纯属浪费计算资源。先用1万次跑通逻辑做敏感性分析时再临时调大就够了。这个项目的代码量不大但每一步都踩在点子上。我自己在做这个研究时最大的体会是分布选型不能只靠“别人都用”或者“画图好看”要回到物理过程去想风速为什么有偏态、辐照度为什么有界、出力为什么有三段式映射把这些因果链条想通了代码只是顺水推舟的事情。后续如果要做更复杂的风光储联合仿真在这个脚本上扩展维度就行——加上储能充放电策略时序化采样代替独立采样并入电力系统潮流计算都是顺着同一套概率框架往下走的路径。
