做电力系统随机调度的朋友应该都遇到过这个场景拿历史风电和光伏出力数据做蒙特卡洛抽样生成的风光联合出力场景包络看着挺丰满一接入优化模型就露馅——备用容量算出来不是保守到浪费就是乐观到不靠谱。问题通常不在采样方法本身而在一个很多人忽略的前提上风电和光伏出力并不是独立变量。这篇就用Matlab从数据到代码完整实现一套基于Copula的风光联合出力场景生成方法把相关性显式建模代码直接能跑换成你自己的数据就能用。1. 风光联合场景生成为什么不能忽略相关性1.1 实测数据里的风光联动现象风电出力靠空气动能光伏出力靠太阳辐照两者在物理上完全是两套机制但放在同一片地理区域内它们受同一个大气过程支配。云团移动、锋面过境、高压控制、台风外围环流这些气象事件同时影响着风速和辐照度于是风电场和光伏电站的出力序列之间就出现了不可忽视的相依性。我拿某地区实测数据举个例子。晴天辐照强、光伏出力高但如果这个晴天是强气压梯度造成的风速也大那光电同时处于高出力区反过来阴雨天气光伏出力趴窝如果恰好是冷涡天气风速反而冲高。所以这两个变量之间既有正相关时段也有负相关时段而且在不同季节、不同天气类型下方向和强度都会变。你如果直接把两条曲线独立建模等于把这种联动关系扔掉了。从数学上看联合概率密度和两个边缘概率密度的乘积只有当变量独立时才相等。实测风光数据很少满足独立条件独立假设下的场景生成本质上是在生成一个与现实分布不同的联合分布后面所有优化结果都建立在错误的分布之上。1.2 独立抽样会让场景失真多少很多人一开始用最简单的独立蒙特卡洛先分别拟合风电出力和光伏出力的边缘分布然后分别抽样、随机配对。这种做法的问题在于随机配对完全忽略了二者的联合行为。量化一下你就明白了。假如风电出力有20%的概率处于高出力区比如超过80%额定容量光伏也有20%的概率处于高出力区。在独立假设下两者同时处于高出力区的概率是0.2乘以0.2等于4%。但如果实际数据中两者呈正相关这个同时高发概率可能是10%甚至更高如果呈负相关实际概率可能只有1%。这一来一回系统充裕度评估和备用容量配置的结果差得不是一星半点。独立抽样生成的场景还有一个典型毛病会出现现实中几乎不可能同时发生的极端组合比如冬季寒潮大风天气下光伏夜里满发——物理上就不可能。这些失真场景进入两阶段随机优化后要么让决策者过度保守多配了根本不需要的备用要么过于乐观把真实风险低估了。无论哪种都不是我们想要的。1.3 Copula建模的思路和边界Copula之所以成为解决这类问题的标准工具是因为它把一个复杂的联合分布建模问题拆成了两个相对独立的部分边缘分布和相依结构。Sklar定理告诉我们任意一个联合分布都可以写成F(x, y) C(F_X(x), F_Y(y))的形式其中C就是一个Copula函数它专门描述变量之间怎么联动。这带来的工程好处很直接。第一边缘分布你可以自由选择风电用Weibull、光伏用Beta、或者直接用非参数核密度估计都行不用为了迁就相关性模型而扭曲边际特征第二相依结构单独用Copula刻画你可以针对性地选择是否建模尾部相关性——也就是极端事件同发的概率第三拟合和采样都有成熟的数值方法Matlab里copulafit和copularnd两个函数就能完成大部分工作。需要说明的是Copula不是万能的。它描述的是变量间的同步变化关系不直接描述时序上的自相关性。如果你要做带时间戳的时序场景比如模拟24小时出力曲线还得额外叠加时序模型比如马尔可夫链或者多维自回归过程。Copula解决的是同一时刻两个变量如何取值这个截面问题搞清楚了这一点后面用起来就不会跑偏。2. Copula模型怎么选原理、对比与判断方法2.1 一句话理解Sklar定理和CopulaCopula这个词听起来抽象其实本质不复杂。它是一个定义在单位正方形上的多元分布函数它的边缘分布都是[0,1]上的均匀分布。只要一个函数满足这个边界条件它就是一个合法的Copula就能用作胶水把多个边缘分布粘成一个联合分布。用个生活类比边缘分布像两扇独立的门每扇门都有自己的开合规律Copula就是连接两扇门的铰链。铰链决定了它们是一起开、反向开、还是各自开。你单独把门板做得多好都不够铰链不对两扇门的联动就不对。风光场景生成也一样风电边缘分布、光伏边缘分布都可以拟合得很准但铰链——Copula——没选对联合行为照样失真。实际操作里Copula场景生成就是三个动作先把历史风电和光伏出力分别转换成均匀分布变量也就是取各自CDF值然后用这些均匀变量拟合Copula参数最后从拟合好的Copula中采样再把采样值逆变换回出力物理量。中间的数学细节很多但主线就是这么清晰。2.2 五种常用Copula族对比工程上常用的Copula族有五个Gaussian、t、Clayton、Gumbel和Frank。它们各有各的相关结构假设选择时主要看三件事相关结构是否对称、是否存在尾部相关性、尾部相关是在上尾还是下尾。Copula族关键参数对称性尾部相关特征适合刻画的情形Gaussian相关矩阵Σ对称无尾部相关常规弱相关极端事件不太同时出现t相关矩阵Σ、自由度ν对称上下尾均有相关极端事件如大风与晴空可能同时出现Clayton参数θ0非对称下尾相关同时低出力、同时出力不足的风险Gumbel参数θ≥1非对称上尾相关同时高出力、消纳困难的风险Frank参数θ对称无尾部相关弱相关、相关结构较均匀的场合尾部相关这个概念值得多说两句。所谓上尾相关是指一个变量处于极大值时另一个变量也处于极大值的条件概率不为零下尾相关同理指一个变量极小值时另一个也极小。Gaussian Copula不管相关系数多高极端尾部都是渐近独立的这意味着它天然低估同时大出力和同时零出力这类小概率高影响事件。对电力系统而言同时零出力可能导致供电不足同时满发可能导致消纳困难这两个尾部恰恰是不能忽略的。2.3 风光场景生成的选型实操建议我自己的选型套路比较固定。先拟合Gaussian和t两种对比拟合对数似然值如果t明显更高说明数据里有尾部相关结构就用t如果两者差别不大优先选Gaussian因为它参数少、数值稳定性好不容易过拟合。如果研究问题聚焦在可靠性方向比如评估极端静稳天气下的电力短缺风险可以额外试试Clayton它的下尾相关能更好地刻画风也没、光也没的场景如果聚焦在消纳方向比如评估大风天和晴空天同时出现造成的弃风弃光风险Gumbel更合适。但要注意Clayton和Gumbel这种阿基米德族Copula在二维情形方便到高维扩展就比较麻烦多维场景还是用t Copula更省事。选型时还有一个容易踩的坑不要只比较拟合优度一定要结合物理背景判断。有时候模型在历史数据上拟合得很好但尾部行为不符合实际风光出力特征生成的上尾场景在物理上就不成立。我的经验是先用数据驱动选型再用物理知识做合理性筛选两者结合才靠谱。3. Matlab完整实现从历史数据到生成场景3.1 数据预处理先把夜间零值问题解决掉做这个流程第一步不是算相关性而是处理数据。光伏出力有大量夜间零值如果你把24小时数据全丢进去边缘分布会在零点处堆一个巨大的概率尖峰后面采样出来的光伏场景会大量出现零出力但这部分零是昼夜规律造成的不是天气不确定性造成的混在一起建模会把相关结构搞得面目全非。我的做法是先把白天数据筛出来再建模。可以按太阳高度角筛选也可以用光伏出力阈值的经验值简单点就是保留日出后和日落前的时间段。风电数据一般没有这种昼夜结构但如果有弃风限电导致的异常低出力段也要考虑剔除或者标记。数据清洗还有一个细节把序列中的NaN和明显错误值处理掉比如负值、超过装机容量的值。这些脏数据会影响边缘分布估计也可能在相关性计算里制造伪相关。处理完的数据统一整理成两个等长的列向量一个风电、一个光伏后面所有步骤都基于这两个向量。% 假设wind和pv是历史出力序列单位MW已对齐 % 先剔除NaN valid ~isnan(wind) ~isnan(pv); wind wind(valid); pv pv(valid); % 光伏夜间零值处理只保留光伏出力大于0.01倍装机容量的时段 pv_cap 100; % 光伏装机容量MW day_mask pv 0.01 * pv_cap; wind wind(day_mask); pv pv(day_mask);3.2 边缘分布拟合核密度估计与转换边缘分布拟合有两种路线参数法和非参数法。参数法对风电用Weibull分布、对光伏用Beta分布是教科书常见做法但实际数据的分布形态经常比这些经典分布复杂尤其在边界附近。我更喜欢用核密度估计也就是ksdensity它不做强分布假设数据长什么样就拟合成什么样灵活性高很多。用ksdensity可以直接得到CDF值也就是把出力值转换成[0,1]区间上的均匀分布变量。这一步是Copula建模的关键前置操作因为Copula的输入要求就是边缘分布均匀化的数据。% 用核密度估计获取历史样本的CDF值 u_w ksdensity(wind, wind, Function, cdf); u_p ksdensity(pv, pv, Function, cdf); % Copula要求数据严格在(0,1)开区间内做个安全钳位 u_w max(min(u_w, 1 - 1e-6), 1e-6); u_p max(min(u_p, 1 - 1e-6), 1e-6); U [u_w, u_p];注意这里有个细节如果样本量很少比如只有几十个点ksdensity的带宽选择可能不稳定这时候用参数法或者增加带宽平滑系数会更稳。样本量上千时核密度方法基本无脑用就行。另外逆变换的时候需要保留CDF的插值节点建议把ksdensity返回的横纵坐标存下来后面有用。3.3 相关性计算与Copula参数拟合在拟合Copula之前建议先算一下Spearman或Kendall相关性系数对数据间的相关强度和方向心里有个底。这两个系数对单调变换保持不变比Pearson更适合Copula这种基于秩的建模。% 计算秩相关系数做初步判断 tau_wp corr(wind, pv, Type, Kendall); rho_sp corr(wind, pv, Type, Spearman); fprintf(Kendall tau: %.3f, Spearman rho: %.3f\n, tau_wp, rho_sp);然后用copulafit拟合Copula参数。Matlab的Statistics and Machine Learning Toolbox里copulafit同时支持Gaussian、t、Clayton、Frank、Gumbel几种常用族输入的是前面得到的U矩阵也就是两列均匀分布变量输出的是对应Copula的参数。% 拟合Gaussian Copula返回相关矩阵 Rho_g copulafit(Gaussian, U); % 拟合t Copula返回相关矩阵和自由度 [Rho_t, Nu_t] copulafit(t, U); % 比较对数似然辅助选型copulafit可以返回第二个输出 [~, loglik_g] copulafit(Gaussian, U); [~, loglik_t] copulafit(t, U); fprintf(Gaussian loglik: %.3f, t loglik: %.3f, Nu_t: %.2f\n, ... loglik_g, loglik_t, Nu_t);自由度的含义值得说一下。t Copula的自由度ν越小尾部相关性越强ν越大t分布越接近Gaussian。如果拟合出来的ν大于30甚至50说明数据里几乎看不出尾部相关直接用Gaussian也完全可以。如果ν只有5到10说明历史数据里确实存在明显的极端值联动这时候用t Copula能更准确地还原联合尾部行为。3.4 采样与逆变换还原出力场景拟合好参数之后下一步就是从Copula中生成大量均匀分布样本再把均匀样本逆变换回出力物理量。copularnd负责采样它生成的样本是N行2列的矩阵每一列都是[0,1]上的均匀分布但两列之间保留了Copula定义的相关结构。N 2000; % 场景数量 V copularnd(t, Rho_t, Nu_t, N); % 从t Copula采样 % 如果选Gaussian用 V copularnd(Gaussian, Rho_g, N); % 钳位到开区间避免逆变换时出现Inf V max(min(V, 1 - 1e-6), 1e-6); % 逆变换利用ksdensity的CDF插值节点还原出力值 % 先获取完整的CDF曲线坐标 [f_w, x_w] ksdensity(wind, Function, cdf); [f_p, x_p] ksdensity(pv, Function, cdf); % 对采样点做逆CDF变换 wind_scen interp1(f_w, x_w, V(:,1), linear, extrap); pv_scen interp1(f_p, x_p, V(:,2), linear, extrap); % 截断到物理合理范围 wind_scen min(max(wind_scen, 0), max(wind)); pv_scen min(max(pv_scen, 0), max(pv));这段代码里最容易出问题的是interp1这步。采样值V如果在[0,1]区间内而f_w覆盖了接近0接近1的整个范围内插没问题但V如果被钳位到1e-6以下或者1-1e-6以上插值就会跑到cdf曲线两端的延拓区域可能插出负值或者超过装机容量的值。所以最后一行的截断不是可有可无是必须做的。生成完场景我习惯先画一张散点图把历史数据的散点图和生成场景的散点图并排对比一下如果形状对得上说明相关结构还原得不错。这一步视觉检查比任何指标都直观。4. 场景削减与质量验证4.1 场景削减的必要性和常用手段直接生成的场景数量通常很大1000个、2000个场景直接接入两阶段随机优化求解规模会爆炸。比如随机机组组合问题场景数翻倍整数变量和约束规模几乎线性增长求解时间可能翻好几倍。所以在工程应用中生成大量原始场景之后几乎总要削减成少量典型场景通常是5到20个再赋上概率。场景削减的核心目标是用少量场景尽可能保留原始场景集的分布特征尤其是均值、方差、相关结构这些对决策影响大的统计量。常用手段有K-means聚类、层次聚类、快速前向选择法。其中K-means简单高效Matlab一行就能调用快速前向选择法在电力系统文献里也很常见但实现稍繁琐对一般场景削减任务K-means足够。4.2 K-means削减的Matlab实现K-means削减的做法是把每个场景当成二维空间里的一个点聚类中心就是典型场景每个聚类的成员数量除以总场景数就是该典型场景的概率。实现很直接。K 10; % 典型场景数量按需调整 [idx, C] kmeans([wind_scen, pv_scen], K, ... Replicates, 15, MaxIter, 500); % 计算每个场景的概率 cnt accumarray(idx, 1, [K, 1]); prob cnt / sum(cnt); % 聚类中心C就是削减后的典型场景 wind_rep C(:,1); pv_rep C(:,2);K的选择可以看手肘图把K从2扫到20记录每个K对应的总组内离差平方和画出来找拐点。但实际工程里不用这么精雕细琢我一般看问题复杂度机组组合问题取10个左右可靠性评估可以取20个再多对结果精度提升有限计算代价却涨得飞快。另外提醒一句如果两个出力变量的量纲差异大比如风电装机1000MW、光伏装机100MW聚类前最好分别除以各自装机容量归一化否则聚类结果会被量纲大的变量主导。4.3 验证场景质量的三个维度场景生成完不验证就是耍流氓。我每次做这个流程必查三个维度边缘统计量、相关结构和概率分布形状。边缘统计量主要看均值和标准差。生成场景的均值、标准差应该和历史数据接近如果差很多说明边缘分布拟合或者逆变换出了问题。相关结构就对比Pearson相关系数和Spearman相关系数历史数据一组值、生成场景一组值K-means削减后的典型场景再算一组值看逐级损失了多少相关性。% 边缘统计量对比 mean_hist [mean(wind), mean(pv)]; mean_scen [mean(wind_scen), mean(pv_scen)]; std_hist [std(wind), std(pv)]; std_scen [std(wind_scen), std(pv_scen)]; fprintf(历史均值: %.2f %.2f, 场景均值: %.2f %.2f\n, ... mean_hist(1), mean_hist(2), mean_scen(1), mean_scen(2)); % 相关性对比 corr_hist corr(wind, pv, Type, Pearson); corr_scen corr(wind_scen, pv_scen, Type, Pearson); corr_rep corr(wind_rep, pv_rep, Type, Pearson); fprintf(Pearson相关: 历史 %.3f, 场景 %.3f, 典型场景 %.3f\n, ... corr_hist, corr_scen, corr_rep);分布形状的验证可以画Q-Q图也可以直接对比历史数据和生成场景的核密度曲线。我更习惯看Q-Q图因为分位数对齐情况一目了然哪个段位偏差大马上就能看出来。如果高尾部分明显偏离45度线说明尾部结构还原得不够好这时候回到Copula选型环节考虑换成尾部相关性更强的模型。5. 实操踩坑记录与问题排查5.1 采样值落在边界导致Inf怎么处理这是我最常看到的问题也是最容易解决的问题。Copula采样值理论上在[0,1]区间但数值计算时可能精确取到0或者1尤其是场景数很大的时候。逆变换时CDF值为0对应出力为0或者负无穷CDF值为1对应正无穷直接插值就会得到Inf后面所有计算全部报错。处理方法就是钳位。统一对采样值做max(min(V, 1 - 1e-6), 1e-6)钳位到开区间内再进逆变换。这个1e-6的经验值在绝大多数场景下不会改变分布特征但能彻底消灭Inf问题。如果数据量特别大或者对尾部精度有苛求可以改成1e-10但没必要追求极端1e-6足够。5.2 别用Pearson系数指导Copula拟合这是一个建模逻辑层面的坑。Copula本质上是对变量排序信息的建模也就是秩相关。而Pearson相关系数衡量的是线性相关它对边缘分布的具体形态敏感。两个变量即使Spearman秩相关系数很高经过不同的非线性变换后Pearson相关系数可能变化很大但秩相关不变。所以拟合Copula之前判断相关性强弱和方向要看Spearman或Kendall不要看Pearson。拟合完成后验证相关性同样以Spearman为主。Pearson可以报告但别把它当成判断Copula拟合质量的主要依据。我见过有人看到Pearson相关系数高就以为建模成功结果秩相关和极端尾部完全对不上后面的优化结果自然也是错的。5.3 全年一锅炖季节性时变相关结构怎么拆风电和光伏的相关性不是全年恒定的。很多地区春季大风和光照往往同时较好夏冬两季可能呈负相关这些时变特征如果被全年数据平均掉生成的场景在特定季节就不符合实际。举个例子冬季傍晚用电高峰时段风力可能因为寒潮增强光伏已经归零全年统一模型很难还原这种季节性差异。处理方法不复杂把历史数据按季节或者逐月切分每个时间段单独拟合边缘分布和Copula参数然后按季节比例混合生成场景。比如生成全年场景时冬季用冬季的Copula参数抽一部分场景夏季用夏季参数抽一部分最后合并。Matlab里写一个for循环按月份分组拟合就行代价不大但场景逼真度提升很明显。判断该不该分季节可以按月份分别算一下Kendall tau如果各月之间差异确实大就分如果都差不多就全年统一拟合省事。数据驱动的判断比拍脑袋可靠。5.4 光伏出力零值堆积怎么建模前面提到夜间零值问题但即使筛掉了夜间数据光伏出力在低辐照时段仍然可能大量出现接近零的值比如阴雨天。这些零值会在边缘分布的低端形成一个概率堆积如果直接进ksdensity边界处理不好会让零值附近的CDF出现异常梯度。处理办法有两个方向。一是用混合模型把光伏出力为零建模成一个离散概率事件再把大于零的出力用一个连续分布刻画二是简单粗筛把出力小于某个小阈值的点都当成一个近似零值类但这会损失部分信息。我的经验是如果只是做场景生成给优化模型用用阈值筛选后直接核密度估计通常够用如果做的是可靠性分析对零出力的概率精度要求高就得上混合模型。5.5 多维变量场景生成怎么办如果场景里不只有一组风电场和光伏电站而是有多个风电场、多个光伏电站甚至加上负荷二元Copula就不够了。这时候有几个选择直接用高维t CopulaMatlab的copulafit和copularnd本身支持多维输入U矩阵有几列就能拟合几维实现成本最低或者用R-Vine Copula它对复杂相关结构的刻画更细腻但Matlab没有原生函数需要自己实现或者借助其他工具工程代价不小。高维t Copula的代价是计算量增长和相关矩阵估计精度下降。样本量不够时高维相关矩阵估计出来可能不是正定的copulafit会报错。这时候可以先用样本协方差阵做修正或者降维处理比如把多个风电场聚合成一个等效风电场光伏同理再建模。工程上聚合方法很多时候已经够用别一上来就上高维模型。6. 这个方法的后续扩展和我的几点体会6.1 从出力场景到预测误差场景实际调度中真正影响决策的是预测误差而不是出力本身。你可以把同样的Copula流程应用在预测误差数据上收集历史风电预测误差和光伏预测误差分别拟合边缘分布再拟合Copula。预测误差通常有偏态比如风电预测误差容易出现负偏核密度估计能抓住这些形态比强行套正态分布或者Weibull分布靠谱得多。预测误差场景对备用配置更有参考价值。因为调度决策关心的是预测值公布后实际出力可能偏离多少这个偏离的联合分布直接决定了上下备用容量的分配。用Copula把风电光伏预测误差的相关结构还原出来备用配置的精度能上一个台阶。6.2 和随机优化模型怎么衔接场景生成完典型做法是给每个典型场景附上概率接入两阶段随机优化的第二阶段。第一阶段决策在看见场景之前做出比如机组启停和预调度计划第二阶段根据具体场景做调整比如实时再调度和切负荷。Copula生成的场景和概率在这里就是随机的离散近似场景数量要跟求解能力匹配10个典型场景对大多数混合整数线性规划问题是个合理的起点。如果你的问题是用鲁棒优化处理不确定性不想要概率场景那Copula的作用可以反过来用找到联合分布下的极端场景比如上尾联合高出力或者下尾联合低出力场景作为鲁棒优化的不确定集合顶点。这个思路在有些文献里叫相关性感知鲁棒优化比单纯用盒式不确定集合贴近实际得多。6.3 几点真实体会我做这类项目也有几年了最后说三个最深的体会。第一Copula不是高级装饰它解决的是实际问题。独立抽样和Copula抽样在最简单的均值方差层面可能看不出太大差别但你去做含概率约束的优化或者评估极端事件风险时差别就非常明显了联合分布的尾部行为直接决定结果的可靠性。第二Matlab里copulafit和copularnd虽然是黑盒但用之前务必搞清楚输入输出约定。特别是U矩阵必须在开区间(0,1)内、不能含0和1这个细节文档里写得清楚但很多人还是会栽在这里。另外不同Matlab版本下ksdensity的逆变换函数写法略有差异老版本可能需要自己用插值实现建议优先用interp1方案兼容性更好。第三也是我最想强调的场景生成是手段不是目的别在场景生成阶段过度追求完美。生成的场景最后要接入优化模型、要支撑决策对决策结果影响最大的通常是最基本的均值方差和相关结构而不是某个高阶统计量的微小差异。先跑通整个流程再回头根据决策结果判断哪里需要细化这个思路比一开始就堆各种复杂模型高效得多。
基于Copula的风光联合出力场景生成：Matlab实现与相关性建模
做电力系统随机调度的朋友应该都遇到过这个场景拿历史风电和光伏出力数据做蒙特卡洛抽样生成的风光联合出力场景包络看着挺丰满一接入优化模型就露馅——备用容量算出来不是保守到浪费就是乐观到不靠谱。问题通常不在采样方法本身而在一个很多人忽略的前提上风电和光伏出力并不是独立变量。这篇就用Matlab从数据到代码完整实现一套基于Copula的风光联合出力场景生成方法把相关性显式建模代码直接能跑换成你自己的数据就能用。1. 风光联合场景生成为什么不能忽略相关性1.1 实测数据里的风光联动现象风电出力靠空气动能光伏出力靠太阳辐照两者在物理上完全是两套机制但放在同一片地理区域内它们受同一个大气过程支配。云团移动、锋面过境、高压控制、台风外围环流这些气象事件同时影响着风速和辐照度于是风电场和光伏电站的出力序列之间就出现了不可忽视的相依性。我拿某地区实测数据举个例子。晴天辐照强、光伏出力高但如果这个晴天是强气压梯度造成的风速也大那光电同时处于高出力区反过来阴雨天气光伏出力趴窝如果恰好是冷涡天气风速反而冲高。所以这两个变量之间既有正相关时段也有负相关时段而且在不同季节、不同天气类型下方向和强度都会变。你如果直接把两条曲线独立建模等于把这种联动关系扔掉了。从数学上看联合概率密度和两个边缘概率密度的乘积只有当变量独立时才相等。实测风光数据很少满足独立条件独立假设下的场景生成本质上是在生成一个与现实分布不同的联合分布后面所有优化结果都建立在错误的分布之上。1.2 独立抽样会让场景失真多少很多人一开始用最简单的独立蒙特卡洛先分别拟合风电出力和光伏出力的边缘分布然后分别抽样、随机配对。这种做法的问题在于随机配对完全忽略了二者的联合行为。量化一下你就明白了。假如风电出力有20%的概率处于高出力区比如超过80%额定容量光伏也有20%的概率处于高出力区。在独立假设下两者同时处于高出力区的概率是0.2乘以0.2等于4%。但如果实际数据中两者呈正相关这个同时高发概率可能是10%甚至更高如果呈负相关实际概率可能只有1%。这一来一回系统充裕度评估和备用容量配置的结果差得不是一星半点。独立抽样生成的场景还有一个典型毛病会出现现实中几乎不可能同时发生的极端组合比如冬季寒潮大风天气下光伏夜里满发——物理上就不可能。这些失真场景进入两阶段随机优化后要么让决策者过度保守多配了根本不需要的备用要么过于乐观把真实风险低估了。无论哪种都不是我们想要的。1.3 Copula建模的思路和边界Copula之所以成为解决这类问题的标准工具是因为它把一个复杂的联合分布建模问题拆成了两个相对独立的部分边缘分布和相依结构。Sklar定理告诉我们任意一个联合分布都可以写成F(x, y) C(F_X(x), F_Y(y))的形式其中C就是一个Copula函数它专门描述变量之间怎么联动。这带来的工程好处很直接。第一边缘分布你可以自由选择风电用Weibull、光伏用Beta、或者直接用非参数核密度估计都行不用为了迁就相关性模型而扭曲边际特征第二相依结构单独用Copula刻画你可以针对性地选择是否建模尾部相关性——也就是极端事件同发的概率第三拟合和采样都有成熟的数值方法Matlab里copulafit和copularnd两个函数就能完成大部分工作。需要说明的是Copula不是万能的。它描述的是变量间的同步变化关系不直接描述时序上的自相关性。如果你要做带时间戳的时序场景比如模拟24小时出力曲线还得额外叠加时序模型比如马尔可夫链或者多维自回归过程。Copula解决的是同一时刻两个变量如何取值这个截面问题搞清楚了这一点后面用起来就不会跑偏。2. Copula模型怎么选原理、对比与判断方法2.1 一句话理解Sklar定理和CopulaCopula这个词听起来抽象其实本质不复杂。它是一个定义在单位正方形上的多元分布函数它的边缘分布都是[0,1]上的均匀分布。只要一个函数满足这个边界条件它就是一个合法的Copula就能用作胶水把多个边缘分布粘成一个联合分布。用个生活类比边缘分布像两扇独立的门每扇门都有自己的开合规律Copula就是连接两扇门的铰链。铰链决定了它们是一起开、反向开、还是各自开。你单独把门板做得多好都不够铰链不对两扇门的联动就不对。风光场景生成也一样风电边缘分布、光伏边缘分布都可以拟合得很准但铰链——Copula——没选对联合行为照样失真。实际操作里Copula场景生成就是三个动作先把历史风电和光伏出力分别转换成均匀分布变量也就是取各自CDF值然后用这些均匀变量拟合Copula参数最后从拟合好的Copula中采样再把采样值逆变换回出力物理量。中间的数学细节很多但主线就是这么清晰。2.2 五种常用Copula族对比工程上常用的Copula族有五个Gaussian、t、Clayton、Gumbel和Frank。它们各有各的相关结构假设选择时主要看三件事相关结构是否对称、是否存在尾部相关性、尾部相关是在上尾还是下尾。Copula族关键参数对称性尾部相关特征适合刻画的情形Gaussian相关矩阵Σ对称无尾部相关常规弱相关极端事件不太同时出现t相关矩阵Σ、自由度ν对称上下尾均有相关极端事件如大风与晴空可能同时出现Clayton参数θ0非对称下尾相关同时低出力、同时出力不足的风险Gumbel参数θ≥1非对称上尾相关同时高出力、消纳困难的风险Frank参数θ对称无尾部相关弱相关、相关结构较均匀的场合尾部相关这个概念值得多说两句。所谓上尾相关是指一个变量处于极大值时另一个变量也处于极大值的条件概率不为零下尾相关同理指一个变量极小值时另一个也极小。Gaussian Copula不管相关系数多高极端尾部都是渐近独立的这意味着它天然低估同时大出力和同时零出力这类小概率高影响事件。对电力系统而言同时零出力可能导致供电不足同时满发可能导致消纳困难这两个尾部恰恰是不能忽略的。2.3 风光场景生成的选型实操建议我自己的选型套路比较固定。先拟合Gaussian和t两种对比拟合对数似然值如果t明显更高说明数据里有尾部相关结构就用t如果两者差别不大优先选Gaussian因为它参数少、数值稳定性好不容易过拟合。如果研究问题聚焦在可靠性方向比如评估极端静稳天气下的电力短缺风险可以额外试试Clayton它的下尾相关能更好地刻画风也没、光也没的场景如果聚焦在消纳方向比如评估大风天和晴空天同时出现造成的弃风弃光风险Gumbel更合适。但要注意Clayton和Gumbel这种阿基米德族Copula在二维情形方便到高维扩展就比较麻烦多维场景还是用t Copula更省事。选型时还有一个容易踩的坑不要只比较拟合优度一定要结合物理背景判断。有时候模型在历史数据上拟合得很好但尾部行为不符合实际风光出力特征生成的上尾场景在物理上就不成立。我的经验是先用数据驱动选型再用物理知识做合理性筛选两者结合才靠谱。3. Matlab完整实现从历史数据到生成场景3.1 数据预处理先把夜间零值问题解决掉做这个流程第一步不是算相关性而是处理数据。光伏出力有大量夜间零值如果你把24小时数据全丢进去边缘分布会在零点处堆一个巨大的概率尖峰后面采样出来的光伏场景会大量出现零出力但这部分零是昼夜规律造成的不是天气不确定性造成的混在一起建模会把相关结构搞得面目全非。我的做法是先把白天数据筛出来再建模。可以按太阳高度角筛选也可以用光伏出力阈值的经验值简单点就是保留日出后和日落前的时间段。风电数据一般没有这种昼夜结构但如果有弃风限电导致的异常低出力段也要考虑剔除或者标记。数据清洗还有一个细节把序列中的NaN和明显错误值处理掉比如负值、超过装机容量的值。这些脏数据会影响边缘分布估计也可能在相关性计算里制造伪相关。处理完的数据统一整理成两个等长的列向量一个风电、一个光伏后面所有步骤都基于这两个向量。% 假设wind和pv是历史出力序列单位MW已对齐 % 先剔除NaN valid ~isnan(wind) ~isnan(pv); wind wind(valid); pv pv(valid); % 光伏夜间零值处理只保留光伏出力大于0.01倍装机容量的时段 pv_cap 100; % 光伏装机容量MW day_mask pv 0.01 * pv_cap; wind wind(day_mask); pv pv(day_mask);3.2 边缘分布拟合核密度估计与转换边缘分布拟合有两种路线参数法和非参数法。参数法对风电用Weibull分布、对光伏用Beta分布是教科书常见做法但实际数据的分布形态经常比这些经典分布复杂尤其在边界附近。我更喜欢用核密度估计也就是ksdensity它不做强分布假设数据长什么样就拟合成什么样灵活性高很多。用ksdensity可以直接得到CDF值也就是把出力值转换成[0,1]区间上的均匀分布变量。这一步是Copula建模的关键前置操作因为Copula的输入要求就是边缘分布均匀化的数据。% 用核密度估计获取历史样本的CDF值 u_w ksdensity(wind, wind, Function, cdf); u_p ksdensity(pv, pv, Function, cdf); % Copula要求数据严格在(0,1)开区间内做个安全钳位 u_w max(min(u_w, 1 - 1e-6), 1e-6); u_p max(min(u_p, 1 - 1e-6), 1e-6); U [u_w, u_p];注意这里有个细节如果样本量很少比如只有几十个点ksdensity的带宽选择可能不稳定这时候用参数法或者增加带宽平滑系数会更稳。样本量上千时核密度方法基本无脑用就行。另外逆变换的时候需要保留CDF的插值节点建议把ksdensity返回的横纵坐标存下来后面有用。3.3 相关性计算与Copula参数拟合在拟合Copula之前建议先算一下Spearman或Kendall相关性系数对数据间的相关强度和方向心里有个底。这两个系数对单调变换保持不变比Pearson更适合Copula这种基于秩的建模。% 计算秩相关系数做初步判断 tau_wp corr(wind, pv, Type, Kendall); rho_sp corr(wind, pv, Type, Spearman); fprintf(Kendall tau: %.3f, Spearman rho: %.3f\n, tau_wp, rho_sp);然后用copulafit拟合Copula参数。Matlab的Statistics and Machine Learning Toolbox里copulafit同时支持Gaussian、t、Clayton、Frank、Gumbel几种常用族输入的是前面得到的U矩阵也就是两列均匀分布变量输出的是对应Copula的参数。% 拟合Gaussian Copula返回相关矩阵 Rho_g copulafit(Gaussian, U); % 拟合t Copula返回相关矩阵和自由度 [Rho_t, Nu_t] copulafit(t, U); % 比较对数似然辅助选型copulafit可以返回第二个输出 [~, loglik_g] copulafit(Gaussian, U); [~, loglik_t] copulafit(t, U); fprintf(Gaussian loglik: %.3f, t loglik: %.3f, Nu_t: %.2f\n, ... loglik_g, loglik_t, Nu_t);自由度的含义值得说一下。t Copula的自由度ν越小尾部相关性越强ν越大t分布越接近Gaussian。如果拟合出来的ν大于30甚至50说明数据里几乎看不出尾部相关直接用Gaussian也完全可以。如果ν只有5到10说明历史数据里确实存在明显的极端值联动这时候用t Copula能更准确地还原联合尾部行为。3.4 采样与逆变换还原出力场景拟合好参数之后下一步就是从Copula中生成大量均匀分布样本再把均匀样本逆变换回出力物理量。copularnd负责采样它生成的样本是N行2列的矩阵每一列都是[0,1]上的均匀分布但两列之间保留了Copula定义的相关结构。N 2000; % 场景数量 V copularnd(t, Rho_t, Nu_t, N); % 从t Copula采样 % 如果选Gaussian用 V copularnd(Gaussian, Rho_g, N); % 钳位到开区间避免逆变换时出现Inf V max(min(V, 1 - 1e-6), 1e-6); % 逆变换利用ksdensity的CDF插值节点还原出力值 % 先获取完整的CDF曲线坐标 [f_w, x_w] ksdensity(wind, Function, cdf); [f_p, x_p] ksdensity(pv, Function, cdf); % 对采样点做逆CDF变换 wind_scen interp1(f_w, x_w, V(:,1), linear, extrap); pv_scen interp1(f_p, x_p, V(:,2), linear, extrap); % 截断到物理合理范围 wind_scen min(max(wind_scen, 0), max(wind)); pv_scen min(max(pv_scen, 0), max(pv));这段代码里最容易出问题的是interp1这步。采样值V如果在[0,1]区间内而f_w覆盖了接近0接近1的整个范围内插没问题但V如果被钳位到1e-6以下或者1-1e-6以上插值就会跑到cdf曲线两端的延拓区域可能插出负值或者超过装机容量的值。所以最后一行的截断不是可有可无是必须做的。生成完场景我习惯先画一张散点图把历史数据的散点图和生成场景的散点图并排对比一下如果形状对得上说明相关结构还原得不错。这一步视觉检查比任何指标都直观。4. 场景削减与质量验证4.1 场景削减的必要性和常用手段直接生成的场景数量通常很大1000个、2000个场景直接接入两阶段随机优化求解规模会爆炸。比如随机机组组合问题场景数翻倍整数变量和约束规模几乎线性增长求解时间可能翻好几倍。所以在工程应用中生成大量原始场景之后几乎总要削减成少量典型场景通常是5到20个再赋上概率。场景削减的核心目标是用少量场景尽可能保留原始场景集的分布特征尤其是均值、方差、相关结构这些对决策影响大的统计量。常用手段有K-means聚类、层次聚类、快速前向选择法。其中K-means简单高效Matlab一行就能调用快速前向选择法在电力系统文献里也很常见但实现稍繁琐对一般场景削减任务K-means足够。4.2 K-means削减的Matlab实现K-means削减的做法是把每个场景当成二维空间里的一个点聚类中心就是典型场景每个聚类的成员数量除以总场景数就是该典型场景的概率。实现很直接。K 10; % 典型场景数量按需调整 [idx, C] kmeans([wind_scen, pv_scen], K, ... Replicates, 15, MaxIter, 500); % 计算每个场景的概率 cnt accumarray(idx, 1, [K, 1]); prob cnt / sum(cnt); % 聚类中心C就是削减后的典型场景 wind_rep C(:,1); pv_rep C(:,2);K的选择可以看手肘图把K从2扫到20记录每个K对应的总组内离差平方和画出来找拐点。但实际工程里不用这么精雕细琢我一般看问题复杂度机组组合问题取10个左右可靠性评估可以取20个再多对结果精度提升有限计算代价却涨得飞快。另外提醒一句如果两个出力变量的量纲差异大比如风电装机1000MW、光伏装机100MW聚类前最好分别除以各自装机容量归一化否则聚类结果会被量纲大的变量主导。4.3 验证场景质量的三个维度场景生成完不验证就是耍流氓。我每次做这个流程必查三个维度边缘统计量、相关结构和概率分布形状。边缘统计量主要看均值和标准差。生成场景的均值、标准差应该和历史数据接近如果差很多说明边缘分布拟合或者逆变换出了问题。相关结构就对比Pearson相关系数和Spearman相关系数历史数据一组值、生成场景一组值K-means削减后的典型场景再算一组值看逐级损失了多少相关性。% 边缘统计量对比 mean_hist [mean(wind), mean(pv)]; mean_scen [mean(wind_scen), mean(pv_scen)]; std_hist [std(wind), std(pv)]; std_scen [std(wind_scen), std(pv_scen)]; fprintf(历史均值: %.2f %.2f, 场景均值: %.2f %.2f\n, ... mean_hist(1), mean_hist(2), mean_scen(1), mean_scen(2)); % 相关性对比 corr_hist corr(wind, pv, Type, Pearson); corr_scen corr(wind_scen, pv_scen, Type, Pearson); corr_rep corr(wind_rep, pv_rep, Type, Pearson); fprintf(Pearson相关: 历史 %.3f, 场景 %.3f, 典型场景 %.3f\n, ... corr_hist, corr_scen, corr_rep);分布形状的验证可以画Q-Q图也可以直接对比历史数据和生成场景的核密度曲线。我更习惯看Q-Q图因为分位数对齐情况一目了然哪个段位偏差大马上就能看出来。如果高尾部分明显偏离45度线说明尾部结构还原得不够好这时候回到Copula选型环节考虑换成尾部相关性更强的模型。5. 实操踩坑记录与问题排查5.1 采样值落在边界导致Inf怎么处理这是我最常看到的问题也是最容易解决的问题。Copula采样值理论上在[0,1]区间但数值计算时可能精确取到0或者1尤其是场景数很大的时候。逆变换时CDF值为0对应出力为0或者负无穷CDF值为1对应正无穷直接插值就会得到Inf后面所有计算全部报错。处理方法就是钳位。统一对采样值做max(min(V, 1 - 1e-6), 1e-6)钳位到开区间内再进逆变换。这个1e-6的经验值在绝大多数场景下不会改变分布特征但能彻底消灭Inf问题。如果数据量特别大或者对尾部精度有苛求可以改成1e-10但没必要追求极端1e-6足够。5.2 别用Pearson系数指导Copula拟合这是一个建模逻辑层面的坑。Copula本质上是对变量排序信息的建模也就是秩相关。而Pearson相关系数衡量的是线性相关它对边缘分布的具体形态敏感。两个变量即使Spearman秩相关系数很高经过不同的非线性变换后Pearson相关系数可能变化很大但秩相关不变。所以拟合Copula之前判断相关性强弱和方向要看Spearman或Kendall不要看Pearson。拟合完成后验证相关性同样以Spearman为主。Pearson可以报告但别把它当成判断Copula拟合质量的主要依据。我见过有人看到Pearson相关系数高就以为建模成功结果秩相关和极端尾部完全对不上后面的优化结果自然也是错的。5.3 全年一锅炖季节性时变相关结构怎么拆风电和光伏的相关性不是全年恒定的。很多地区春季大风和光照往往同时较好夏冬两季可能呈负相关这些时变特征如果被全年数据平均掉生成的场景在特定季节就不符合实际。举个例子冬季傍晚用电高峰时段风力可能因为寒潮增强光伏已经归零全年统一模型很难还原这种季节性差异。处理方法不复杂把历史数据按季节或者逐月切分每个时间段单独拟合边缘分布和Copula参数然后按季节比例混合生成场景。比如生成全年场景时冬季用冬季的Copula参数抽一部分场景夏季用夏季参数抽一部分最后合并。Matlab里写一个for循环按月份分组拟合就行代价不大但场景逼真度提升很明显。判断该不该分季节可以按月份分别算一下Kendall tau如果各月之间差异确实大就分如果都差不多就全年统一拟合省事。数据驱动的判断比拍脑袋可靠。5.4 光伏出力零值堆积怎么建模前面提到夜间零值问题但即使筛掉了夜间数据光伏出力在低辐照时段仍然可能大量出现接近零的值比如阴雨天。这些零值会在边缘分布的低端形成一个概率堆积如果直接进ksdensity边界处理不好会让零值附近的CDF出现异常梯度。处理办法有两个方向。一是用混合模型把光伏出力为零建模成一个离散概率事件再把大于零的出力用一个连续分布刻画二是简单粗筛把出力小于某个小阈值的点都当成一个近似零值类但这会损失部分信息。我的经验是如果只是做场景生成给优化模型用用阈值筛选后直接核密度估计通常够用如果做的是可靠性分析对零出力的概率精度要求高就得上混合模型。5.5 多维变量场景生成怎么办如果场景里不只有一组风电场和光伏电站而是有多个风电场、多个光伏电站甚至加上负荷二元Copula就不够了。这时候有几个选择直接用高维t CopulaMatlab的copulafit和copularnd本身支持多维输入U矩阵有几列就能拟合几维实现成本最低或者用R-Vine Copula它对复杂相关结构的刻画更细腻但Matlab没有原生函数需要自己实现或者借助其他工具工程代价不小。高维t Copula的代价是计算量增长和相关矩阵估计精度下降。样本量不够时高维相关矩阵估计出来可能不是正定的copulafit会报错。这时候可以先用样本协方差阵做修正或者降维处理比如把多个风电场聚合成一个等效风电场光伏同理再建模。工程上聚合方法很多时候已经够用别一上来就上高维模型。6. 这个方法的后续扩展和我的几点体会6.1 从出力场景到预测误差场景实际调度中真正影响决策的是预测误差而不是出力本身。你可以把同样的Copula流程应用在预测误差数据上收集历史风电预测误差和光伏预测误差分别拟合边缘分布再拟合Copula。预测误差通常有偏态比如风电预测误差容易出现负偏核密度估计能抓住这些形态比强行套正态分布或者Weibull分布靠谱得多。预测误差场景对备用配置更有参考价值。因为调度决策关心的是预测值公布后实际出力可能偏离多少这个偏离的联合分布直接决定了上下备用容量的分配。用Copula把风电光伏预测误差的相关结构还原出来备用配置的精度能上一个台阶。6.2 和随机优化模型怎么衔接场景生成完典型做法是给每个典型场景附上概率接入两阶段随机优化的第二阶段。第一阶段决策在看见场景之前做出比如机组启停和预调度计划第二阶段根据具体场景做调整比如实时再调度和切负荷。Copula生成的场景和概率在这里就是随机的离散近似场景数量要跟求解能力匹配10个典型场景对大多数混合整数线性规划问题是个合理的起点。如果你的问题是用鲁棒优化处理不确定性不想要概率场景那Copula的作用可以反过来用找到联合分布下的极端场景比如上尾联合高出力或者下尾联合低出力场景作为鲁棒优化的不确定集合顶点。这个思路在有些文献里叫相关性感知鲁棒优化比单纯用盒式不确定集合贴近实际得多。6.3 几点真实体会我做这类项目也有几年了最后说三个最深的体会。第一Copula不是高级装饰它解决的是实际问题。独立抽样和Copula抽样在最简单的均值方差层面可能看不出太大差别但你去做含概率约束的优化或者评估极端事件风险时差别就非常明显了联合分布的尾部行为直接决定结果的可靠性。第二Matlab里copulafit和copularnd虽然是黑盒但用之前务必搞清楚输入输出约定。特别是U矩阵必须在开区间(0,1)内、不能含0和1这个细节文档里写得清楚但很多人还是会栽在这里。另外不同Matlab版本下ksdensity的逆变换函数写法略有差异老版本可能需要自己用插值实现建议优先用interp1方案兼容性更好。第三也是我最想强调的场景生成是手段不是目的别在场景生成阶段过度追求完美。生成的场景最后要接入优化模型、要支撑决策对决策结果影响最大的通常是最基本的均值方差和相关结构而不是某个高阶统计量的微小差异。先跑通整个流程再回头根据决策结果判断哪里需要细化这个思路比一开始就堆各种复杂模型高效得多。
