风光联合出力场景生成：Copula建模在Matlab中的完整实现

搞新能源并网计算的同学，大概率都撞过这么一堵墙：手上明明有风电场和光伏电站的实测功率数据，做随机优化的时候要生成风光出力场景，脑子里第一反应就是把风电、光伏当成两个互不干扰的独立变量，分别采样再随机拼在一起。算出来的方案看着也挺像回事，可一旦拿到真实数据上去回验，结果偏得让人怀疑人生。问题往往不在优化算法，而在场景生成这一步——你把两个场站出力之间的相关性弄丢了。

同一片区域里的风资源和光照资源，背后驱动它们的往往是同一个天气过程。一场冷锋过境，风速上来了，云层也跟着压下来，光伏出力反而往下掉；晴天的时候光照拉满，但风可能小得可怜。这种“此消彼长”的耦合关系是客观存在的，不建模，优化模型就会生成大量现实中根本不会出现的组合，比如“夜里光伏满发”“大风天同时艳阳高照”。Copula就是专门解决这类问题的工具：它把每个变量的边缘分布和变量之间的相依结构拆开建模，既能保留风光各自的出力统计特性，又能还原二者之间的真实相关性，然后在这个框架下批量生成场景。

这篇文章我把整套思路和Matlab实现从头到尾拆开讲，包括Copula是什么、为什么适合风光联合出力、如何用Matlab完成边缘拟合-Copula参数估计-场景采样-逆变换还原的全流程，以及我实际踩过的坑。适合正在做新能源出力场景生成、含风光随机规划的科研人员和工程师参考。

1. 为什么必须把风光联合出力放在一起建模

1.1 风光出力的相关性是物理规律，不是统计巧合

很多人刚开始接触这个题目时会有个疑问：风电和光伏一个靠风、一个靠光，原理上八竿子打不着，凭什么说它们相关？把两个场站的实测数据放在散点图上看就明白了，风大时光伏通常偏低，光伏高时风往往偏小，两者呈现明显的负相关趋势。这不是数据巧合，而是气象过程在背后起的作用。

区域性天气系统会同时影响风速和太阳辐照度。比如典型的大气过程过境时，云量增加、风速增强，此时风电出力上升、光伏出力下降；而高压控制下的晴朗天气，辐照度达到峰值，但地面风速往往较弱。这种物理上的负相关在时间尺度上还表现出昼夜和季节特性：夜间没有光伏出力但风速可能很大，白天光照充足但风速相对平缓。单看任何一条曲线都是随机波动，但放在一起看，它们是被同一个“天气开关”联动的。

如果忽略这种联动关系，在场景生成阶段把风电和光伏各自独立采样，就会凭空造出一批物理上不可能同时出现的场景。这些“伪场景”一旦进入优化模型，轻则让调度方案偏保守、浪费可用的新能源出力，重则让可靠性评估结果乐观到危险的程度。做电网规划和储能配置的人应该深有体会，极端场景恰恰是最不能错的。

1.2 相关性建模的三种思路，为什么选Copula

处理两个变量相关性的常规做法大致有三类，我用一个表格把它们的思路和适用场景对比清楚：

Copula的优势在于“解耦”。Sklar定理告诉我们，任意一个二维联合分布函数C(x₁, x₂)，都可以拆成两个边缘分布F₁(x₁)、F₂(x₂)和一个连接函数Cop(u₁, u₂)的复合。用大白话说，每个变量的“长相”由边缘分布负责，变量之间“怎么勾连”由Copula函数负责，两者互不干扰。

对应到风光场景生成里，这意味着我可以给风电出力随便选一个拟合得很好的边缘分布（比如Weibull分布），给光伏出力选另一个边缘分布（比如Beta分布），然后再用Copula把两者的“排位关系”绑定起来。相比“一刀切”的多元正态假设，这种灵活性在处理非高斯、非对称的风光出力数据时是决定性的。而且Copula能直接生成任意数量的新样本，不需要依赖历史数据重排，这对Monte Carlo类的随机优化来说是刚需。

2. Copula建模的核心概念与选型逻辑

2.1 copula怎么理解：概率积分变换是钥匙

要真正会用Copula，先得打通一个概率论里的基础概念——概率积分变换。它说的是：对一个连续性随机变量X，把它的累积分布函数F(x)套在自己身上，得到的F(X)服从[0,1]区间上的均匀分布。反过来，从[0,1]均匀分布抽一个数u，再用F⁻¹(u)就能还原出一个服从原分布的X。这个过程直白点说，就是把“物理空间”的随机数，搬到“单位方形空间”里去做操作，用完之后再搬回去。

Copula就是定义在这个单位方形空间上的联合分布函数。它管的不是变量本身的数值大小，而是变量的“概率排位”如何相互作用。比如Wind出力在历史数据里排在30%的位置、Solar排在70%的位置，这种排位共现的规律，就是Copula要刻画的。

场景生成走的是这条路的逆过程：第一步，从拟合好的Copula里抽出成千上万个(u_w, u_s)对；第二步，用各自的边缘分布逆函数把u还原成物理出力值。这样得到的场景，既保持了单变量的出力分布特性，又继承了变量之间的排位相关结构。整套流程里最微妙的地方是，边缘拟合得准不准和Copula选得好不好，各占一半的成败。

2.2 Copula族怎么选：看尾部相关性就对了一大半

Matlab里常用Copula族有五个，我按特性和适用场景分一下：

选型的核心观察点就是“极端情况怎么联动”。把历史数据画成散点图，盯着左下角和右上角看：如果低出力极端值经常同时出现（风小、云厚、光弱），Clayton值得优先试；如果高出力极端值经常抱团，Gumbel更合适；如果上下尾都有明显的联动，t-Copula是安全牌。

实际项目中我不主张拍脑袋定族，而是把五个族都拟合一遍，比较对数似然值和AIC。AIC = -2 × loglik + 2 × k，k是Copula的参数个数，AIC越小说明模型在拟合数据和复杂度之间平衡得越好。这一步很机械，但有数据支撑的选型总比肉眼判断靠谱。

2.3 从Kendall秩相关系数反推Copula参数

拟合Copula参数最常用的方法是最大似然估计，但理解参数和秩相关之间的解析关系，能让你调试过程事半功倍。以Archimedean族为例：

• Clayton：Kendall tau = α / (α + 2)
• Gumbel：Kendall tau = 1 - 1/α
• Frank：tau = 1 - 4/α × [D₁(α) - 1]，D₁是第一阶Debye函数，需要数值求解

这意味着你不需要跑任何优化，直接用历史数据的Kendall tau就能算出一个不错的参数初值。比如实测数据的Kendall tau为-0.3，用Gumbel族时α = 1/(1-tau) = 0.769，发现小于1，说明Gumbel根本不适用——因为Gumbel只能刻画正相关。这种“先用tau粗筛、再用似然精估”的操作，能帮你少走很多弯路。

3. Matlab实现：从原始数据到风光场景的完整流程

3.1 建模总览：五步走完场景生成

完整的Copula场景生成在Matlab里分五步，我先把流程列出来，后面逐段展开：

1. 数据准备和预处理：取同一时段的风电、光伏出力历史序列，等长、对齐时间戳，处理缺失值和异常值。
2. 边缘分布拟合：为风电出力、光伏出力分别选择合适的边缘分布并估计参数。
3. 概率积分变换：把历史出力数据变换到[0,1]均匀空间，得到U矩阵。这一步是为Copula拟合做准备。
4. Copula参数拟合与选型：在均匀空间拟合各Copula族参数，用AIC选择最优模型。
5. 场景采样与逆变换：从选定的Copula中采样，再逆变换还原为物理出力场景，必要时做场景削减。

整个流程的逻辑，本质上就是“去物理化-建关联-再物理化”的闭环。

3.2 边缘分布拟合：先说一个容易翻车的地方

风光出力数据里有个很扎心的特点——大量零值。风电在无风时段出力为零，光伏在夜间出力为零。直接把原始序列丢进weibullfit拟合，估计出的参数会被一堆零带偏，后面PIT变换时也会出问题。

我的做法是把出力分成“零出力”和“正出力”两部分来建模：正出力部分拟合连续分布，零出力部分作为一个离散概率质量单独记录。以风电为例，代码是这样写的：

% 数据预处理：输入 P_wind, P_solar 均为 n×1 列向量 F0_w = mean(P_wind == 0); % 风电零出力比例 F0_s = mean(P_solar == 0); % 光伏零出力比例（注意夜间数据占比） % 只取正出力样本拟合Weibull分布 pw_pos = P_wind(P_wind > 0); ps_pos = P_solar(P_solar > 0); [Aw, Bw] = wblfit(pw_pos); % 风电正出力 Weibull 形状、尺度参数 [As, Bs] = wblfit(ps_pos); % 光伏正出力 Weibull 形状、尺度参数

之所以选Weibull，是因为风速分布通常用Weibull描述，风电功率由风速经由功率曲线转换后，整体偏态特征仍然和Weibull比较契合。光伏出力更常见的备选是Beta分布，实际项目里可以做几次拟合优度检验再定。如果数据本身形态复杂，用ksdensity做非参数核密度估计也是一种思路，但要注意边界处理问题——出力下限为零，核密度估计会在边界处“泄漏”出负值，需要额外修正。

混合分布建模的核心是记住一句话：零出力是一个离散概率事件，不能简单粗暴地丢进连续分布里。这一步偷懒，后面相关性结构一定失真。

3.3 概率积分变换：核心U矩阵的构造

边缘分布拟合完成后，就要把物理出力数据映射到[0,1]均匀空间。这里的关键技巧是：零出力样本不能全部映射到同一个点，否则会在Copula拟合时把所有零出力时刻的相关性放大成伪相关。

我的做法是把零出力的均匀值随机散布在[0, F0]区间内，正出力样本映射到(F0, 1]区间。这样既保留了“零出力概率为F0”的统计事实，又不会在均匀空间里产生一坨顽固的重复点：

n = length(P_wind); U_wind = zeros(n, 1); idx_w = P_wind > 0; U_wind(idx_w) = F0_w + (1 - F0_w) * wblcdf(P_wind(idx_w), Aw, Bw); U_wind(~idx_w) = F0_w * rand(nnz(~idx_w), 1); U_solar = zeros(n, 1); idx_s = P_solar > 0; U_solar(idx_s) = F0_s + (1 - F0_s) * wblcdf(P_solar(idx_s), As, Bs); U_solar(~idx_s) = F0_s * rand(nnz(~idx_s), 1); U = [U_wind, U_solar];

这段代码里最值得琢磨的是U_wind(idx_w) = F0_w + (1-F0_w) × wblcdf(...)这一行。含义是正出力样本的累积概率不再从0算起，而是从F0_w开始继续累积。举个例子，如果风电零出力占比20%，那么一个正出力值对应的均匀值最小也是0.2，这正好保证了“零出力”和“正出力”在均匀空间里不会重叠。

零出力样本的随机散布用的是rand函数，也就是说零出力时刻的U值在[0, F0)内部是随机均匀排列的。这是合理的，因为零出力在物理空间里没有大小次序，在概率空间里也不该有确定的排位。

3.4 Copula参数拟合与最优模型选择

U矩阵到手后，Copula拟合就水到渠成了。Matlab的统计工具箱提供了copulafit函数，一行一个族。我把五个族全部拟合一遍并计算AIC：

% Copula参数估计 rho_gau = copulafit('Gaussian', U); % 相关矩阵 [rho_t, nu_t] = copulafit('t', U, 'Method', 'ML'); % 相关矩阵 + 自由度 alpha_clay = copulafit('Clayton', U); alpha_gum = copulafit('Gumbel', U); alpha_frank = copulafit('Frank', U); % 计算对数似然与AIC [~, ll_gau] = copulaloglik(U, 'Gaussian', rho_gau); [~, ll_t] = copulaloglik(U, 't', rho_t, nu_t); [~, ll_clay] = copulaloglik(U, 'Clayton', alpha_clay); [~, ll_gum] = copulaloglik(U, 'Gumbel', alpha_gum); [~, ll_frank] = copulaloglik(U, 'Frank', alpha_frank); loglik = [ll_gau, ll_t, ll_clay, ll_gum, ll_frank]; k = [1, 2, 1, 1, 1]; % 各Copula族参数个数 AIC = -2 * loglik + 2 * k; AIC_table = table({'Gaussian'; 't'; 'Clayton'; 'Gumbel'; 'Frank'}, ... loglik', k', AIC', ... 'VariableNames', {'Family', 'LogLik', 'k', 'AIC'}); disp(AIC_table);

选AIC最小的那个族。实际数据里，风光出力通常表现出下尾联动（静稳天气下风小、云厚、光弱同时出现），所以Clayton或t-Copula经常会胜出，Gaussian也常见，但Gumbel用于风光负相关场景时往往直接被Kendall tau的符号筛掉。

这里提醒一个细节：copulafit内部其实是要做数值优化的，如果U矩阵里的值有0或1这种端点值，计算会不稳定甚至报错。前面我们处理零出力样本时保证了值域在(0,1)开区间内，恰好规避了这个风险。

3.5 场景采样与逆变换还原

模型选完之后，真正“生成场景”的动作发生在copularnd这一步。它从拟合好的Copula中抽出指定数量的均匀空间样本：

Ns = 5000; % 先采5000个原始场景 Usim = copularnd('t', rho_t, nu_t, Ns);

得到的Usim是Ns×2的矩阵，每一行代表一个仿真时刻下风电和光伏在均匀空间中的“概率排位”。接下来要逆变换回物理空间，这里必须和3.2节、3.3节的混合分布处理保持对称：

% 逆变换：风电 Pw_sim = zeros(Ns, 1); idx_w = Usim(:, 1) > F0_w; Pw_sim(idx_w) = wblinv((Usim(idx_w, 1) - F0_w) ./ (1 - F0_w), Aw, Bw); % 逆变换：光伏 Ps_sim = zeros(Ns, 1); idx_s = Usim(:, 2) > F0_s; Ps_sim(idx_s) = wblinv((Usim(idx_s, 2) - F0_s) ./ (1 - F0_s), As, Bs);

这段逆变换的逻辑是正变换的反向操作。Usim中落在(0, F0_w]区间的样本，逆变换后就是零出力场景；落在(F0_w, 1)区间的样本，先减去F0_w再除以(1-F0_w)，相当于把正出力区间重新拉伸到(0,1)，然后查找wblinv对应的物理出力值。这保证了边界的一致性——你永远不会在采样结果里得到负出力，也不会把零出力样本错误地变成正出力。

到这里，Pw_sim和Ps_sim就是最终的风光联合出力场景集合了。你可以把它们送入随机优化模型，也可以做场景削减后再送。

3.6 生成质量的自检：先验再看相关性

场景生成完不能直接拿去交付，先做两步自检。第一步看边缘分布是否还原：

% 边缘分布对比：用QQ图或经验CDF对比 figure; subplot(1,2,1); ecdf(P_wind); hold on; ecdf(Pw_sim); legend('历史风电','生成场景'); subplot(1,2,2); ecdf(P_solar); hold on; ecdf(Ps_sim); legend('历史光伏','生成场景');

第二步看相关性结构是否还原：

tau_orig = corr(U, 'Type', 'Kendall'); tau_sim = corr(Usim, 'Type', 'Kendall'); disp(['历史Kendall tau: ', num2str(tau_orig(1,2))]); disp(['场景Kendall tau: ', num2str(tau_sim(1,2))]);

如果tau_sim和tau_orig偏离较大，通常不是Copula选错了，而是边缘分布拟合或者零值处理出了问题。相关性结构的还原度，是检验整个模型链路是否正确的试金石。

4. 场景削减去重与代表性检验

4.1 为什么生成5000个场景还要削减

随机优化问题里，每增加一个场景，求解规模就成倍增长。5000个原始场景直接丢进混合整数规划模型，计算时间会让你怀疑人生。实际操作里，标准流程是先粗采样生成大量原始场景，再用场景削减技术浓缩成20到50个代表性场景，每个场景附带一个概率权重。

削减的本质是找一个“近似分布”：用较少的支撑点逼近原始场景集在概率空间中的整体形态，同时最小化概率距离指标。最常用的概念是Wasserstein距离，直观理解就是“把所有场景的概率质量搬到削减后场景上所需的最小代价”。

4.2 kmeans聚类场景削减的简洁实现

最容易被接受的削减手段是kmeans聚类。思路很简单：把5000个二维场景点聚成K类，每类的聚类中心就是一个代表场景，每类里场景数量占总数的比例就是该场景的概率。

K = 20; % 削减后保留的场景数 [idx, C] = kmeans([Pw_sim, Ps_sim], K, 'Replicates', 20); prob = histcounts(idx, K)' / Ns; % 输出削减后的场景集 Pw_scen = C(:, 1); % K×1，代表场景的风电出力 Ps_scen = C(:, 2); % K×1，代表场景的光伏出力 % prob 就是每个代表场景的概率，可直接用于随机优化

kmeans的代价函数基于欧氏距离，对风光出力这种量纲相同的变量效果尚可。如果希望更贴近概率分布的最优传输理论，可以考虑同步回代消除法（scenario reduction）：每次迭代删掉一个概率权重最小的场景，把它的概率合并到距离最近的一个场景上，直到只剩K个场景。这个方法Matlab没有内置函数，自己实现也不复杂，核心是一个贪心循环，但对中等规模数据来说kmeans通常已经能获得足够好的代表性，不必纠结理论最优。

4.3 削减结果怎么验证：三点检查法

削减做完不是看一眼散点图就完事，我习惯做三点检查：

第一，削减后场景集的风电均值、光伏均值与历史数据均值偏差应小于某个容忍度，比如3%；第二，削减后场景集的风光相关系数应与历史秩相关接近，偏差建议控制在0.05以内；第三，削减后场景集的累积分布函数应与历史经验CDF偏差不大，可用KS检验的p值做参考。如果这几项都过了，削减场景才具备进入优化模型的资格。

在实际项目中我也常遇到“削减后场景太集中、丢掉了极端情况”的问题，尤其是K选得太小的时候。这时可以适当增大K，或者在kmeans聚类时对边界场景做强制保留——先把历史数据中最极端的几个点选为锚点，再把剩余场景聚类。极端场景对电力系统的可靠性评估至关重要，宁可多保留几个场景也不要一刀切聚类。

5. 实操中踩过的坑和排查经验实录

5.1 零出力数据直接把Copula拟合搞崩

我第一次做这个项目时，直接把含零出力的原始序列做了ksdensity核密度拟合，然后代入PIT变换。结果U矩阵里出现了一大批完全等于0或1的端点值，copulafit跑出来的参数明显异常，t-Copula的自由度直接飙到几百。排查半天才发现是端点值导致优化过程中似然函数出现奇异。

后来学乖了：所有出力边界都按“离散质量 + 连续分布”的混合模型处理。不只是零出力，光伏出力上限附近如果有大量限电截断值，同样需要把“正出力达到上限”作为一个离散事件处理。核心思想是一致的——边界处的概率质量不能硬塞进连续分布里。

5.2 场景相关性比历史数据弱，问题多半出在边缘拟合

有次我生成完5000个场景，算出来Kendall tau只有历史数据的一半，怎么调Copula族都救不回来。后来把边缘分布的QQ图调出来一看，尾部拟合非常差，核密度估计在正出力区间中段出现了一个不自然的凹谷，把分布形态带歪了。

Copula对边缘分布的质量极其敏感，因为PIT变换完全依赖F(x)的准确性。F(x)差一点，U的排位就跟着歪，相关性自然对不上。所以遇到相关性还原度不高，先查边缘，再查Copula，不要一上来就换族。先用AIC选定Copula族没错，但边缘分布本身的拟合优度检验（KS检验、AD检验）在建模早期就该做扎实。

5.3 场景数量取多少，取决于下游优化问题

这个问题没有标准答案，完全看下游模型能承受多少决策变量。做两阶段随机规划时，我通常先把场景削减到10~50个，太大了对求解器不友好，太小了分布代表性又不够。一个实用的操作是做一个“场景数量敏感性分析”：分别用10、20、30、50个场景求解优化问题，观察目标函数值的变化，当目标值随场景数增加趋于稳定时，就说明这个数量已经够了。

5.4 别忽略风光的时序联动：同一时段的联合出力只是第一步

标题里说的是“考虑风光联合出力”，指的是同一时刻的风光出力相关性。但实际电力系统调度还要考虑时间维度上的连续性，今天风大的时候，明天风可能也不小。如果只是独立地对每个时段生成联合场景，场景之间的时序关系是断裂的。

更完善的方案是构造一个包含时间相关性的场景生成框架：一种做法是先对历史数据按季节、天气类型聚类，在每个簇内用Copula建模联合出力，然后通过马尔可夫链或场景树控制时段之间的切换规律。另一种做法是用时序Copula，直接把相邻时段的风光出力放在同一个多维Copula里建模，维度从2扩展到2T。后者参数估计更复杂，计算量也更大。我的建议是，如果做日前调度，先按小时分开建24个Copula模型，每个时段内处理风光联合出力，时段之间再用事后调整或场景树衔接，工程上更可控。

5.5 常见问题速查表

最后再分享一点个人习惯：任何场景生成类的代码，我都会从头到尾保持随机数种子的可控性（用rng设置种子），这样别人复现时结果完全一致。做科研投稿时这是硬要求，做工程项目时也能让你在调试时不会因为“这次随机数不同”而浪费一晚上排查时间。这套Copula流程我前后在多个风光基地数据上验证过，只要边缘分布和零值处理这两关把住了，场景生成的整体稳定性是很有保障的。