做新能源并网分析或者微电网规划的时候绕不开的一道坎就是风光出力的不确定性建模。风速这种非线性很强的自然量业界用Weibull分布来近似光伏出力因为有天然的[0,1]区间约束Beta分布是出现频率最高的选择。这两个分布放在一起做组合研究在Matlab里实现其实并不复杂但我见过太多人在数据预处理和参数估计上翻车一出结果就问为什么拟合出来是一条直线或者为什么生成的样本全是负数。这篇博文把我自己跑过的完整流程拆开写在这里风速的Weibull拟合、光伏出力的Beta拟合、两个分布的联合抽样以及基于组合模型算出来的几个工程指标。代码全部给到可以直接复制运行的版本参数也给了实际案例里的合理取值。适合刚接触新能源建模、想用Matlab复现文献结果的初学者也适合需要快速搭建风光联合出力模型的规划工程师参考。1. 为什么风电和光伏要用分布建模从“出力是个随机数”说起1.1 风电功率为什么常用Weibull分布风电功率来源于风速而风速是一个典型的非负、右偏、长尾的随机变量。你拿一年的小时级风速数据画直方图会看到大量低速时段集中在左侧偶尔出现阵风或极端大风拖出一条细长的右尾。这种形状用正态分布去拟合是非常别扭的因为正态分布对称且覆盖负无穷到正无穷而风速不可能为负强行拟合会在左端冒出大量不该出现的负风速样本。Weibull分布恰好专治这类问题。它的概率密度函数形式是f(x) (k/c) * (x/c)^(k-1) * exp(-(x/c)^k)其中x大于等于0k称为形状参数c称为尺度参数两者都大于0。k决定了分布的“形态”k越小分布越偏、长尾越明显k越大分布越趋近对称甚至接近正态。c则相当于总体风速的量级参考c越大整体风速水平越高。以我处理过的一个平原风电场为例70米高度年平均风速约6.8米每秒拟合出来的Weibull参数是k约2.2、c约7.5画出来就是一个右偏但不算夸张的峰形与实测直方图重合度很高。另一个实用原因是Weibull分布有显式的累积分布函数F(x) 1 - exp(-(x/c)^k)意味着从分布变换、分位数求解都特别方便不用像一些复杂分布那样靠数值积分“硬算”。这对后续做场景生成、算极端风速重现期都是很大的优势。1.2 光伏出力为什么常用Beta分布光伏出力的物理源头是太阳辐照度而辐照度在给定时段内的波动受云层遮挡、大气透明度影响往往表现出有界、偏态的分布特性。更关键的是光伏出力经过归一化除以额定容量之后天然落在[0,1]区间内Beta分布的定义域恰好就是[0,1]两者在数学结构上高度吻合。Beta分布的概率密度函数是f(x) x^(alpha - 1) * (1-x)^(beta - 1) / B(alpha, beta)alpha和beta是两个正的形状参数它们共同决定分布的均值、方差和形状。Beta分布的均值是alpha/(alphabeta)方差是alphabeta/((alphabeta)^2(alphabeta1))。这意味着你只要从历史出力数据算出均值和方差就能反解出alpha和beta得到一个能反映该场站出力水平的解析模型。这个性质在做蒙特卡洛抽样时特别有用不需要保存几千条原始出力曲线只需要两个参数就能完整描述日内出力的随机特性。我还想提醒一点Beta分布对“有界”数据的刻画能力是正态分布给不了的。正态分布假设数据可以取负值或超过1对于比例型数据出力占比、负载率、可用率会产生大量不合理样本所以很多文献在描述光伏出力、储能SOC状态、设备可用率时都优先选Beta分布逻辑就在这里。1.3 “组合研究”到底研究什么把风速拟合出Weibull、把光伏出力拟合出Beta这只是第一步。真正有价值的“组合研究”是把这两个分布放在同一个框架里回答一个工程问题风光联合出力系统的总出力分布是什么样子。这里有两层含义。第一层是独立组合假设风速和辐照度互不相关分别从Weibull和Beta中抽样把风功率和光伏功率相加得到总出力的经验分布。这个方法简单适合做初步估算。第二层是相关组合同一地区的风资源和光资源往往存在互补性比如多云天气辐照度下降但风速可能上升晴天辐照度强但往往风小。忽略这种相关性会把风险概率算偏所以需要用Copula函数把两个边际分布“绑定”起来生成带有实际相关性特征的联合样本。组合建模的输出物通常是总出力概率密度函数、不同置信水平下的出力分位数比如P95、出力不足某个门槛值的概率以及给储能容量配置用的一连串模拟场景。可以说单独的分布拟合解决的是“单机描述”问题组合研究解决的是“系统描述”问题后者才是规划决策真正需要的输入。2. 风电侧建模Weibull分布的Matlab实现与风速-功率转化2.1 数据准备与最大似然拟合先看数据形态。做Weibull拟合时输入建议用小时级平均风速序列单位统一成米每秒。如果用10分钟级数据也行只是样本量大了拟合会更精细计算成本也高一些。数据清洗这一步不能省剔除明显异常的负风速、检查缺失值可用fillmissing做线性插补。有些站点记录里会出现静风时段的大量0值如果0值占比过高超过5%两参数Weibull的拟合质量会明显下降后面我会单独说这个坑。参数估计首选最大似然法Matlab里直接调用fitdist即可% 风速数据windData单位 m/s一列向量 pd_wind fitdist(windData, Weibull); k pd_wind.A; % 形状参数 k c pd_wind.B; % 尺度参数 c disp([k, c]);fitdist内部使用的就是wblfit函数也就是对对数似然函数做数值优化得到使样本出现概率最大的k和c。除了最大似然法也可以用手动的最小二乘法或者矩估计但MLE的小样本表现最好也是文献默认做法没有必要自己造轮子。如果你手头数据量特别大比如几十万条还想看看拟合迭代过程可以改用mle函数并设置优化选项。不过我从实际经验说fitdist在这种场景下足够稳而且代码可读性好交给同事都能很快看懂。2.2 拟合优度检验看一眼就知道模型行不行拟合完不是结束必须验证模型和原始数据的偏差在可接受范围。最直观的做法是先画叠加图figure; histogram(windData, 50, Normalization, pdf, FaceAlpha, 0.4); hold on; x 0:0.1:max(windData); y wblpdf(x, k, c); plot(x, y, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); legend(实测直方图, Weibull拟合);如果拟合得好红色曲线应该贴着直方图的轮廓走峰值位置接近右尾不会差出好几倍。画图之外还建议做一个定量检验比如计算理论CDF和经验CDF的相关系数R²或者直接用K-S检验[h, p] kstest(windData, CDF, makedist(Weibull, a, k, b, c));p值大于0.05说明不能拒绝原假设模型可以接受。但K-S检验对样本量极其敏感数据几万条时几乎必然显著所以不要只看p值结合图形判断更重要。以我的经验平原和海上风电场两参数Weibull通常够用如果有明显双峰比如山谷风效应造成昼夜两个风速模态可以考虑混合Weibull或三参数Weibull但在组合研究里边际模型过拟合对联合分布影响很小不建议为了好看牺牲简单性。2.3 风速到功率分段函数与三次方模型分布拟合得到的是风速样本但工程上要的是功率。风速和功率的关系由风电机组的功率曲线决定常见的简化模型是分段函数四个关键点分别是切入风速v_in、额定风速v_r、切出风速v_out和额定功率P_rfunction P windPowerCurveVec(v, v_in, v_r, v_out, Pr) P zeros(size(v)); % 切入到额定之间可用线性或者三次方关系 idx (v v_in) (v v_r); P(idx) Pr .* (v(idx) - v_in) / (v_r - v_in); % 额定到切出之间满发 idx2 (v v_r) (v v_out); P(idx2) Pr; % 小于切入或大于切出出力为0 % P默认是0所以不用额外赋值 end注意我特意用向量化写法而不是for循环。数据量是8760条乘上几千个场景循环版本会慢得让人怀疑人生向量化版本几乎瞬时完成。这个函数里用的是线性上升段很多文献也会用三次方模型P Pr * (v/v_r)^3物理上更贴近贝兹理论的功率-风速三次方关系。实际选哪种最好查对应机型手册里的功率曲线。我一般在学术研究里明白写了“简化线性模型”在工程交付时则一定用厂家提供的离散功率曲线插值两者差距在额定风速附近能到5%以上不能忽视。有了这个函数后续从Weibull样本生成风功率序列就非常顺先wblrnd抽样再逐点过功率曲线最后乘以机组台数得到场站总功率。3. 光伏侧建模Beta分布的Matlab实现与出力标幺化3.1 数据归一化与参数估计光伏出力的原始数据一般是kW或MW序列。Beta分布要求数据落在(0,1)开区间第一步必须归一化。标准的做法是除以额定容量或逆变器容量上限得到标幺值序列pvPu。这里有个关键细节如果直接用包含夜间数据的完整序列会出现大量零值而Beta分布的密度函数在端点0和1处不是有限值最大似然估计会直接失效甚至报错。所以拟合Beta分布时我一般只取“有光照时段”也就是剔除辐照度很小或出力等于0的数据点或者至少单独处理。估值方法有两个可选。一是直接调betafit% pvPu光伏出力标幺值取值在(0,1)之间 phat betafit(pvPu); alpha phat(1); beta phat(2);二是用矩估计公式和代码都简单适合想手动验证参数关系的场景mu mean(pvPu); varx var(pvPu, 1); % 注意用总体方差 alpha mu * (mu * (1 - mu) / varx - 1); beta (1 - mu) * (mu * (1 - mu) / varx - 1);矩估计的结果和MLE通常很接近尤其是在数据量上万时。它的好处是完全透明你能清楚看到均值方差如何映射到形状参数。如果算出来的alpha和beta出现负数或非数基本可以断定是输入数据不在(0,1)区间或者包含了过多的0和1回头检查数据比检查代码更有效。3.2 Beta分布抽样与概率密度复现参数到手后生成光伏出力随机样本很简单N 10000; % 抽样数 pvSim betarnd(alpha, beta, N, 1); % 生成N个标幺值样本 pvSimMW pvSim * pvRated; % 换算回实际功率betarnd是Matlab内置的Beta分布随机数生成器底层基于Gamma分布变换速度和稳定性都没问题。如果你更习惯统一调用random函数写成random(Beta, alpha, beta, N, 1)也完全等价。抽样完成后还是要画图对比这一步我基本每次都做。把原始标幺值直方图和模拟样本直方图叠在一起看中心位置、离散程度是否符合预期。如果模拟数据的方差明显小于原始数据多半是矩估计分母用了样本方差除N-1导致参数算偏或者数据里混入了异常尖峰。这一点和很多人直觉相反Beta分布看似灵活其实对方差极其敏感方差差一点形状就变很多所以在计算均值方差时务必用一致的公式。3.3 温度与转换效率的工程修正严格来说光伏出力不只由辐照度决定组件温度、逆变器效率、灰尘遮挡都会影响最终输出。常见的修正模型是在理想线性关系基础上乘一个温度系数P P_stc * (G/G_stc) * [1 - kT * (T_cell - 25)]其中G_stc是标准测试辐照度1000 W每平方米kT是温度系数单晶硅组件通常为0.003到0.004每摄氏度。如果需要更精确可以用逆变器效率曲线再做一次分段折减。但要说实话在分布建模这种统计尺度下我大多数时候不做温度修正原因是Beta分布本身已经把辐照度的随机波动吸收了温度的慢变影响会体现在数据分布形态里不需要额外建立物理模型去复现。分群拆分会更有价值。比如把数据按夏季、冬季分开拟合得到两组alpha、beta比强行在整个样本里塞一个复杂的物理修正更实用。太阳高度角、季节日照时长这种系统性差异交给分布参数去表达就够了。4. 风光组合从独立抽样到联合出力分析4.1 独立组合模拟总出力分布与场景生成两个边际模型都建好后最直接的组合方式是独立抽样。假设风电装机100 MW光伏装机80 MW目标是把两者总出力的概率分布刻画出来rng(42); % 固定随机种子保证结果可复现 N 100000; % 风电抽样风速过功率曲线换算成MW vSim wblrnd(k, c, N, 1); PwindSim windPowerCurveVec(vSim, 3, 12, 25, 100); % 单机/场站模型 % 这里Pr100当作单机额定实际场站再乘以台数下面用总装机时直接给总额定 % 光伏抽样标幺值换算成MW PpvSim betarnd(alpha, beta, N, 1) * 80; % 总出力 Ptotal PwindSim PpvSim; % 统计特征 meanTotal mean(Ptotal); p95 prctile(Ptotal, 95); p05 prctile(Ptotal, 5); fprintf(总出力均值%.2f MWP95%.2f MWP05%.2f MW\n, meanTotal, p95, p05);这里的N取10万是为了让尾部分位数稳定。如果只画直方图几千个样本也够但计算P95、P99这种尾部指标时样本量不够会有很大跳动。我建议固定随机种子这样每次运行结果一致调参数时不会因为随机性干扰判断。独立组合的优点是直观、快速适合做敏感性分析。缺点是它假设风速和辐照度互不关联如果研究场站确实存在明显的天气耦合关系独立组合会低估或高估某些不利场景的概率。这时就该上Copula。4.2 用Copula考虑风光的互补相关性Copula的核心思想是把“边际分布”和“相关性结构”分开建模先用历史数据估计相关性再把这个相关性“绑”到已经拟合好的Weibull和Beta边际上。操作上三步走第一步把历史风速和光伏出力数据变换成均匀边际。这一步通常用经验CDF或者核密度CDFU_wind ksdensity(windData, windData, function, cdf); U_pv ksdensity(pvPu, pvPu, function, cdf); % 注意pvPu和windData要按时间戳对齐第二步拟合Copula参数。我用Gaussian Copula比较多Matlab里直接rho copulafit(Gaussian, [U_wind, U_pv]);rho就是两个变量的秩相关系数矩阵能直接看到风光出力的相关方向和强度。负值表示互补正值表示同向波动。第三步生成带相关性的联合样本再用逆CDF变换回原始分布U_new copularnd(Gaussian, rho, N); vFinal wblinv(U_new(:,1), k, c); pvFinal betainv(U_new(:,2), alpha, beta);代码里的核心逻辑是U_new已经带上了历史数据的秩相关性wblinv和betainv再把均匀样本映射回风速和光伏出力的真实量纲这样就得到了一组“分布形状正确且相关性接近实测”的联合场景。选哪种Copula也要看数据。样本量大、尾部相关性明显时选t Copula效果更好一般工程场景Gaussian和Frank就够用。我个人建议先用Gaussian跑通流程等结果被业务方接受了再考虑更复杂的族避免一上来就在Copula选型上花太多时间。4.3 组合模型能算哪些工程指标联合出力样本生成后能算的东西就多了。最基础的是总出力的概率密度和经验分布可以直接得到任意分位数其次是场景缩减把10万样本聚类成几十个典型场景供时序生产模拟使用再次是概率性指标比如总出力低于某个阈值比如系统负荷的20%的概率P95和P05构成的置信带宽度代表出力波动风险最大出力与最小出力的极差分布用于评估调节压力。这些指标在储能容量配置、备用容量安排、输配电规划里都是核心输入。和传统“取典型日”的做法相比分布组合模型最大的优势是能把风险量化到概率上而不是拍脑袋定一个保守系数。比如某地区风光联合出力P95是132 MW意味着有95%的时段总出力不超过这个值那在制定外送计划时就可以把132 MW作为参考上限比直接按装机容量打七折要科学得多。5. 典型应用场景储能容量配置与可靠性评估5.1 一个简化的储能削峰填谷示例储能容量的粗略估算可以这样建模假设系统需要维持的负荷基线是L_base风光总出力低于基线时缺额由储能补充高于基线时多余电量给储能充电。利用前面生成的联合出力样本缺额功率序列就是max(L_base - Ptotal, 0)L_base 50; % 基础负荷MW deficit max(L_base - Ptotal, 0); % 储能额定功率取缺额分布的95%分位数 P_ess prctile(deficit, 95); % 储能容量用平均缺额乘以时长粗略估算这里是示意 E_ess mean(deficit) * 24; fprintf(建议储能功率%.2f MW储能容量%.2f MWh\n, P_ess, E_ess);这个式子简单得有点“粗暴”但作为规划阶段的量级估算完全够用。真要落地还需要加入SOC约束、充放电效率曲线、日内时序衔接那就要把抽样样本按时间顺序排列做时序仿真而不是只做分布层面的统计。用分布模型做储能估算的好处是可以对比不同置信水平下的功率容量需求取P90会比P95小不少经济性差异很大决策者可以在这个区间里权衡投资成本和供电可靠性。这种量化权衡能力是单一确定值方法给不了的。5.2 概率性指标P95、长时间尺度缺损概率LOLD除了储能风光联合出力模型还可以算可靠性指标。传统电力可靠性里有个指标叫LOLDLoss of Load Duration失负荷时间是统计一个周期内系统无法满足负荷需求的小时数或天数。用分布组合模型可以快速估算它的概率版loadLossProb mean(Ptotal L_base); % 如果全年8760小时运行失负荷小时数近似为 loldHours loadLossProb * 8760;进一步可以算全年失负荷小时数的分布方法是对抽样样本做分组重采样或者直接用二项分布近似。这些指标的共同点是它们把“风光不确定性”通过概率框架传导到“系统可靠性”上便于和常规电源的强迫停运率放在同一套模型里比较。实际项目里我用这个方法给一个园区微电网做过评估风光总装机120 MW负荷基线80 MW独立抽样算出的失负荷概率是3.2%而用Copula模型考虑风光负相关后失负荷概率降到2.4%。差距将近一个百分点对可靠性要求高的用户来说就是完全不同的投资决策。这组对比足以说明相关结构不能省。6. 常见问题与调试经验速查6.1 拟合报错与参数边界问题最常遇到的报错一是fitdist提示找不到工具包。Weibull和Beta的拟合在Matlab的Statistics and Machine Learning Toolbox里没装这个工具箱时fitdist、betafit、kstest这些函数全部不可用。遇到这种情况要么装上工具箱要么用fminsearch手动写最大似然目标函数。我见过有人为了一个拟合函数去重装整个Matlab其实检查工具箱状态用license(test, Statistics_Toolbox)一条命令就能确认。二是betafit报错“数据必须位于(0,1)区间”。原因基本是归一化没做好或者数据里混入了0和1。处理办法对0值场景单独建模对恰好为1的满发数据做微小压缩比如pvPu(pvPu 1) 1 - 1e-6。这么做的代价是损失一点尾部精度但换来MLE能正常收敛整体利大于弊。三是wblfit出现奇异梯度或迭代不收敛。多半是数据里有极端离群值或者是输入数据量太少。先画盒图找出离群点判断是测量错误还是真实物理极端值。真实极端值不要随便删否则拟合出的c会被拉低低估高频大风的风险。6.2 数据预处理最容易踩的坑预处理这块第一个坑是风速数据里的0值。完全的静风在低海拔场站其实罕见大量0值通常来自测风仪低于启动阈值。两参数Weibull对大量0值非常敏感k会被拉得很小拟合曲线左端翘得很高。如果场景确实存在大量静风时段建议先把0值按一个极小正值比如0.01替换再拟合或者干脆改用三参数Weibull增加位置参数来吸收这个偏移。第二个坑是单位不统一。风速用千米每小时功率用千瓦光伏用标幺值三个量纲混在一起组合结果怎么看怎么怪。我的习惯是进入模型前统一用国际单位风速一律米每秒有功功率一律兆瓦标幺值单独保留一列并标注基准值。这个习惯帮我避免过无数次低级错误。第三个坑是时间对齐。风速和光伏数据来自不同SCADA系统时时区、采样间隔经常不一致。没对齐就进入Copula拟合秩相关系数会被严重拉偏。我用过一个笨但有效的办法按整点对齐先对齐再清洗确保每条记录的时间戳都有对方场的对应数据。6.3 组合模拟结果不合理的排查思路如果组合抽样结果均值对不上先查单位。比如风电装机明明是100 MW代码里功率曲线额定值却写成100 kW算出来的均值会差三个数量级。这种问题一眼就能看穿但人处于连续调试状态时反而容易忽略。如果方差明显偏小大概率是归一化时除以的不是额定容量而是最大值导致方差被压缩。另一种可能是随机数生成范围出错比如betarnd生成的样本被误乘了功率曲线变成标幺值直接相加后又乘了错误基准。如果Copula生成样本后秩相关系数和原始数据对不上优先检查经验边际变换这一步。ksdensity的cdf输出有时会在边界处出现轻微偏移样本量太小时尤其明显。可以改用ecdf手动做阶梯变换虽然不平滑但秩相关结构更稳定。排查系统性不合理的场景我最推荐的做法是先固定随机种子跑一组小样本手动算几个中间量核对。比如先只抽100个风速样本过功率曲线后手算期望风速、期望功率和wblrnd的理论均值对比。这个流程五分钟内能找到大多数问题比反复重跑大样本高效得多。我个人在实际操作中的体会是分布建模这个活数学模型本身并不难难的是每一步都保持对数据的敏感度。拟合出一个漂亮的曲线很容易但真正理解这个曲线为什么长这样、在什么场景下会失效才是项目能不能落地的关键。最后再分享一个小技巧给所有中间结果命名时带上分布参数和数据范围比如pvBeta_a2b5_2019to2023三个月后回来看项目文件你会感谢当初自己的这个习惯。