簡介本資源為基于Copula理論與K-means聚類的風(fēng)光聯(lián)合場景生成MATLAB代碼包面向計(jì)算機(jī)、電子信息工程、數(shù)學(xué)等專業(yè)的學(xué)生及科研人員可用于課程設(shè)計(jì)、期末大作業(yè)與畢業(yè)設(shè)計(jì)中的新能源出力場景建模與聚類分析。壓縮包共15個(gè)文件約4.77MB包含1個(gè)main.m主程序、2個(gè)xlsx風(fēng)光功率數(shù)據(jù)文件、3份Copula與場景生成相關(guān)PDF學(xué)習(xí)資料、8張png運(yùn)行結(jié)果圖及說明文本結(jié)構(gòu)清晰、注釋明細(xì)支持參數(shù)化編程與靈活改參。代碼兼容MATLAB 2014、2019a與2024a版本附贈(zèng)案例數(shù)據(jù)可直接運(yùn)行幫助讀者快速復(fù)現(xiàn)從邊緣分布擬合、Copula聯(lián)合分布構(gòu)建到K-means場景聚類的完整流程理解變量間復(fù)雜依賴關(guān)系的建模思路。目前已有122人學(xué)習(xí)下載適合需要掌握場景生成方法并完成相關(guān)課題的中高級學(xué)習(xí)者參考。1. 風(fēng)光聯(lián)合場景生成到底在算什么從 Copula 到 Kmeans 的一條完整鏈路做新能源出力建模的人遲早會(huì)撞上一個(gè)問題風(fēng)電場和光伏電站的出力曲線單獨(dú)看都還好一旦要放在同一個(gè)系統(tǒng)里做聯(lián)合調(diào)度或可靠性評估兩條曲線之間的相關(guān)性就成了繞不過去的坎。你如果直接把風(fēng)電和光伏的歷史出力序列各自擬合一個(gè)邊緣分布再獨(dú)立采樣拼在一起得到的場景集合會(huì)嚴(yán)重失真——因?yàn)楝F(xiàn)實(shí)中太陽落山前后光伏驟降的那段時(shí)間風(fēng)速往往也在變化兩者不是各走各的。這就是風(fēng)光聯(lián)合場景生成要解決的核心問題在保留每個(gè)變量自身統(tǒng)計(jì)特性的同時(shí)把變量之間的相關(guān)結(jié)構(gòu)也復(fù)現(xiàn)出來。標(biāo)題里的 64Copula 指的是一種具體的 Copula 建模方案Kmeans 聚類則負(fù)責(zé)把生成出來的海量場景壓縮成少量典型場景供后續(xù)優(yōu)化模型使用。整條鏈路是歷史出力數(shù)據(jù) → 邊緣分布擬合 → Copula 相關(guān)結(jié)構(gòu)建模 → 采樣生成大量場景 → Kmeans 聚類削減 → 輸出典型場景及其概率。適合做電力系統(tǒng)規(guī)劃、微網(wǎng)容量配置、儲(chǔ)能優(yōu)化調(diào)度的工程師也適合正在寫相關(guān)方向論文、需要一套能跑通的 MATLAB 代碼做實(shí)驗(yàn)的研究生。下面我把這條鏈路拆開從原理到代碼到踩坑一步步講清楚。2. Copula 建模風(fēng)光相關(guān)性的原理與 MATLAB 實(shí)現(xiàn)2.1 為什么不能用協(xié)方差矩陣直接描述風(fēng)光相關(guān)性很多人第一反應(yīng)是用多元正態(tài)分布來建模風(fēng)電和光伏的聯(lián)合分布協(xié)方差矩陣直接描述相關(guān)性。這個(gè)思路在變量服從正態(tài)分布時(shí)沒問題但風(fēng)電出力典型分布是 Beta 或 Weibull光伏出力在白天呈 Beta 分布、夜間恒為零兩者都不是正態(tài)的。硬套多元正態(tài)邊緣分布就被扭曲了尾部相關(guān)性也完全對不上。Copula 的思路是把聯(lián)合分布拆成兩部分邊緣分布各自擬合變量之間的依賴結(jié)構(gòu)單獨(dú)用一個(gè) Copula 函數(shù)描述。Sklar 定理保證了這種拆分的數(shù)學(xué)合法性——任何聯(lián)合分布都可以寫成邊緣分布和 Copula 函數(shù)的組合。這樣一來風(fēng)電用 Weibull 擬合、光伏用 Beta 擬合兩者之間的相關(guān)結(jié)構(gòu)用 Gaussian Copula 或 t-Copula 來刻畫各管各的互不干擾。選 Gaussian Copula 還是 t-Copula如果風(fēng)光出力序列在極端天氣下表現(xiàn)出明顯的尾部依賴比如連續(xù)陰天伴隨低風(fēng)速t-Copula 更合適因?yàn)樗懈竦奈膊?。如果?shù)據(jù)量不大、尾部特征不明顯Gaussian Copula 參數(shù)少、估計(jì)穩(wěn)定是更穩(wěn)妥的起點(diǎn)。我一般先用 Gaussian 跑通再看尾部擬合效果決定要不要換。2.2 邊緣分布擬合風(fēng)電 Weibull、光伏 Beta 的參數(shù)估計(jì)先看風(fēng)電。風(fēng)速通常用 Weibull 分布擬合但注意這里要建模的是風(fēng)電出力而非風(fēng)速。如果只有出力數(shù)據(jù)可以直接對出力做 Beta 擬合或者用風(fēng)速-功率曲線反推。常見做法是對歸一化后的風(fēng)電出力序列做 Beta 擬合形狀參數(shù)用最大似然估計(jì)。光伏出力白天歸一化后同樣適合 Beta 分布。但夜間出力為零這部分不能直接丟——丟了會(huì)改變?nèi)諆?nèi)的時(shí)序結(jié)構(gòu)。處理方式有兩種一是只對白天時(shí)段建模夜間單獨(dú)標(biāo)記為零二是對全時(shí)段做零膨脹 Beta 擬合。代碼里通常采用第一種簡單且不容易翻車。% 假設(shè) wind_pu 和 solar_pu 是歸一化后的風(fēng)電和光伏出力列向量 % 去除夜間零值只對光伏白天時(shí)段擬合 day_idx solar_pu 0; solar_day solar_pu(day_idx); % 風(fēng)電 Beta 擬合 wind_beta fitdist(wind_pu, Beta); a_w wind_beta.a; b_w wind_beta.b; % 光伏 Beta 擬合僅白天 solar_beta fitdist(solar_day, Beta); a_s solar_beta.a; b_s solar_beta.b; % 輸出參數(shù)供后續(xù) Copula 使用 fprintf(Wind Beta: a%.4f, b%.4f\n, a_w, b_w); fprintf(Solar Beta: a%.4f, b%.4f\n, a_s, b_s);這段代碼的關(guān)鍵點(diǎn)fitdist是 MATLAB 統(tǒng)計(jì)工具箱的函數(shù)返回的分布對象可以直接調(diào)cdf和pdf。參數(shù)a和b是 Beta 分布的兩個(gè)形狀參數(shù)控制分布的偏度和峰度。風(fēng)電出力的 Beta 擬合通常 a 較小、b 較大說明出力偏低的時(shí)間居多光伏白天出力的 a 和 b 則取決于當(dāng)?shù)厝照諚l件。擬合完一定要畫 QQ 圖或 PP 圖檢查光看參數(shù)數(shù)值沒有直觀判斷力。2.3 Gaussian Copula 的參數(shù)估計(jì)與采樣步驟邊緣分布擬合好之后下一步是把原始數(shù)據(jù)變換到均勻分布空間。具體做法是對每個(gè)變量做概率積分變換U F_wind(x_wind)V F_solar(x_solar)。變換后 U 和 V 都服從 [0,1] 均勻分布它們之間的相關(guān)性就完全由 Copula 函數(shù)承載了。Gaussian Copula 的參數(shù)就是一個(gè)相關(guān)系數(shù)矩陣。對于二元情況就是一個(gè)標(biāo)量 ρ。估計(jì)方法用 Kendall 秩相關(guān)系數(shù) τ 換算ρ sin(πτ/2)。這個(gè)換算關(guān)系是 Gaussian Copula 的解析性質(zhì)比直接算 Pearson 相關(guān)系數(shù)更穩(wěn)健因?yàn)橹认嚓P(guān)不受邊緣分布形式影響。% 概率積分變換 U cdf(wind_beta, wind_pu); V cdf(solar_beta, solar_day); % 注意長度對齊 % 用 Kendall tau 估計(jì) Gaussian Copula 參數(shù) tau corr(U, V, Type, Kendall); rho sin(pi * tau / 2); % 生成 N 個(gè)聯(lián)合場景 N 10000; rng(42); % 固定隨機(jī)種子保證可復(fù)現(xiàn) Z mvnrnd([0 0], [1 rho; rho 1], N); U_sim normcdf(Z(:,1)); V_sim normcdf(Z(:,2)); % 逆變換回原始空間 wind_sim icdf(wind_beta, U_sim); solar_sim icdf(solar_beta, V_sim);邏輯說明mvnrnd生成二元標(biāo)準(zhǔn)正態(tài)樣本協(xié)方差矩陣?yán)锏?ρ 就是 Copula 參數(shù)。normcdf把正態(tài)樣本變成均勻樣本這一步等價(jià)于從 Gaussian Copula 中采樣。最后用icdf做逆變換把均勻樣本映射回風(fēng)電和光伏的原始出力空間。rng(42)是固定隨機(jī)種子做實(shí)驗(yàn)時(shí)務(wù)必加上否則每次跑出來的場景集合都不一樣沒法對比。參數(shù)方面N 取 10000 是經(jīng)驗(yàn)值——太少則聚類效果不穩(wěn)定太多則后續(xù) Kmeans 計(jì)算量大。ρ 的值一般在 0.2 到 0.6 之間取決于當(dāng)?shù)仫L(fēng)資源和光照的日內(nèi)互補(bǔ)程度。如果算出來 ρ 是負(fù)數(shù)說明風(fēng)電和光伏在日內(nèi)呈互補(bǔ)關(guān)系白天光照好、風(fēng)速低這在很多地區(qū)是正常現(xiàn)象不用覺得模型出了問題。3. Kmeans 場景削減從一萬條曲線到五個(gè)典型日3.1 場景削減為什么必須做以及 Kmeans 的適用邊界Copula 采樣出來一萬個(gè)場景直接丟進(jìn)優(yōu)化模型是不現(xiàn)實(shí)的——計(jì)算量扛不住。場景削減的目標(biāo)是用少量典型場景加對應(yīng)概率盡可能保留原始場景集合的統(tǒng)計(jì)特征。Kmeans 是最常用的手段因?yàn)樗鼘?shí)現(xiàn)簡單、收斂快而且聚類中心本身就是一條完整的出力曲線物理意義明確。但 Kmeans 有它的邊界。它假設(shè)簇是球形的、各簇方差相近對于風(fēng)光出力這種時(shí)序曲線如果不同季節(jié)的曲線形態(tài)差異很大比如夏季和冬季的日內(nèi)模式完全不同Kmeans 可能把形態(tài)差異大的曲線硬分到同一簇。這時(shí)候可以考慮先用層次聚類做初步劃分再在每簇內(nèi)用 Kmeans 細(xì)化。不過對于大多數(shù)工程應(yīng)用直接對歸一化后的曲線做 Kmeans 已經(jīng)夠用。另一個(gè)關(guān)鍵決策是聚類數(shù) K 怎么定。常用的方法有肘部法、輪廓系數(shù)法、以及基于場景削減誤差的指標(biāo)。我一般先用肘部法粗選再結(jié)合后續(xù)優(yōu)化模型對場景數(shù)的容忍度來定。如果優(yōu)化模型是線性規(guī)劃K 取 5 到 10 都行如果是混合整數(shù)規(guī)劃K 最好不超過 5否則求解時(shí)間會(huì)明顯拉長。3.2 把 24 小時(shí)出力曲線喂給 Kmeans特征構(gòu)造與歸一化Kmeans 的輸入是特征矩陣每一行是一個(gè)場景每一列是一個(gè)特征。對于 24 小時(shí)出力曲線最直接的做法是把 24 個(gè)時(shí)間點(diǎn)的出力值作為 24 個(gè)特征。但這樣有個(gè)問題不同場景的出力幅值差異可能很大Kmeans 基于歐氏距離幅值大的特征會(huì)主導(dǎo)距離計(jì)算。解決辦法是先對每條曲線做歸一化讓所有曲線在同一尺度下比較形態(tài)。歸一化方式有兩種一是最大最小值歸一化到 [0,1]二是 Z-score 標(biāo)準(zhǔn)化。對于風(fēng)光出力曲線我傾向于用最大最小值歸一化因?yàn)槌隽Ρ旧碛形锢砩辖珙~定容量歸一化后保留了這個(gè)信息。% 假設(shè) wind_sim 和 solar_sim 是 N x 24 的矩陣每行一個(gè)場景 % 拼接成聯(lián)合場景矩陣 joint_scenes [wind_sim, solar_sim]; % N x 48 % 對每條曲線做最大最小值歸一化 scene_norm zeros(size(joint_scenes)); for i 1:size(joint_scenes, 1) row joint_scenes(i, :); scene_norm(i, :) (row - min(row)) / (max(row) - min(row) eps); end % 或者直接對每一列做歸一化保留時(shí)間維度的形態(tài) scene_norm (joint_scenes - min(joint_scenes)) ./ ... (max(joint_scenes) - min(joint_scenes) eps);這里有個(gè)容易搞混的地方按行歸一化和按列歸一化效果完全不同。按行歸一化是把每條曲線自己的幅值信息抹掉只保留形態(tài)按列歸一化是保留每個(gè)時(shí)間點(diǎn)的相對大小。做場景削減通常用按行歸一化因?yàn)槲覀冴P(guān)心的是曲線形態(tài)的相似性而不是絕對幅值。但如果你后續(xù)優(yōu)化模型需要絕對出力值那聚類完還要把幅值信息還原回去。3.3 聚類數(shù) K 的選取與典型場景概率計(jì)算K 的選取沒有一刀切的標(biāo)準(zhǔn)。我通常跑一個(gè)循環(huán)K 從 2 到 10每次記錄簇內(nèi)平方和然后畫肘部圖。肘部拐點(diǎn)對應(yīng)的 K 就是候選值。但肘部法有時(shí)候拐點(diǎn)不明顯這時(shí)候結(jié)合輪廓系數(shù)一起看。K_range 2:10; sse zeros(length(K_range), 1); sil zeros(length(K_range), 1); for k K_range [idx, C, sumd] kmeans(scene_norm, k, Replicates, 5, MaxIter, 500); sse(k-1) sum(sumd); sil(k-1) mean(silhouette(scene_norm, idx)); end % 畫肘部圖 figure; yyaxis left; plot(K_range, sse, -o); ylabel(SSE); yyaxis right; plot(K_range, sil, -s); ylabel(Silhouette); xlabel(K);Replicates設(shè)為 5 表示用 5 個(gè)不同的初始中心各跑一次取最優(yōu)結(jié)果避免陷入局部最優(yōu)。MaxIter是最大迭代次數(shù)500 對于 24 維特征通常夠用。silhouette計(jì)算輪廓系數(shù)值越接近 1 說明聚類效果越好。聚類完成后每個(gè)簇的典型場景就是聚類中心。該場景的概率等于簇內(nèi)樣本數(shù)除以總樣本數(shù)。注意聚類中心是在歸一化空間里的如果要還原到原始出力空間需要用對應(yīng)簇內(nèi)樣本的最大最小值做逆變換。% 假設(shè)最終選定 K5 K 5; [idx, C] kmeans(scene_norm, K, Replicates, 10); % 計(jì)算每個(gè)典型場景的概率 prob zeros(K, 1); for k 1:K prob(k) sum(idx k) / length(idx); end % 還原到原始空間以風(fēng)電為例 wind_typical zeros(K, 24); for k 1:K cluster_samples joint_scenes(idx k, 1:24); wind_typical(k, :) mean(cluster_samples, 1); end這里用簇內(nèi)樣本的均值作為典型場景比直接用歸一化空間的聚類中心更合理因?yàn)榫当Ａ袅嗽挤敌畔?。概率?jì)算就是簡單的頻率統(tǒng)計(jì)所有簇的概率之和為 1。4. 避坑與排查風(fēng)光聯(lián)合場景生成里最容易翻車的五個(gè)地方4.1 概率積分變換后出現(xiàn) Inf 或 NaN現(xiàn)象調(diào)用cdf做概率積分變換時(shí)結(jié)果里出現(xiàn) Inf 或 NaN后續(xù) Copula 參數(shù)估計(jì)直接報(bào)錯(cuò)。原因邊緣分布擬合的參數(shù)不合理導(dǎo)致某些數(shù)據(jù)點(diǎn)的 CDF 值趨近于 0 或 1取對數(shù)或做逆變換時(shí)溢出。常見于 Beta 分布的邊界值處理——如果原始數(shù)據(jù)里有 0 或 1Beta 分布的 CDF 在邊界處可能返回極端值。解決在擬合前對數(shù)據(jù)做微小的截?cái)啾热绨?0 替換為 1e-6把 1 替換為 1-1e-6。另外檢查fitdist的擬合結(jié)果是否收斂如果參數(shù)估計(jì)時(shí)出現(xiàn)警告換用矩估計(jì)法重新擬合。4.2 Kendall 相關(guān)系數(shù)算出來是負(fù)數(shù)以為代碼寫錯(cuò)了現(xiàn)象風(fēng)電和光伏的 Kendall τ 是負(fù)值ρ 也是負(fù)的采樣出來的場景里風(fēng)電高的時(shí)候光伏低看起來“反了”。原因這不是代碼錯(cuò)誤而是數(shù)據(jù)本身的特征。很多地區(qū)風(fēng)光在日內(nèi)確實(shí)呈互補(bǔ)關(guān)系——白天光照強(qiáng)、風(fēng)速低夜間風(fēng)速高、無光照。負(fù)相關(guān)是合理的。解決不用改代碼但要在論文或報(bào)告里解釋清楚這個(gè)負(fù)相關(guān)的物理含義。如果確實(shí)預(yù)期是正相關(guān)檢查數(shù)據(jù)對齊是否正確——比如風(fēng)電和光伏的時(shí)間戳是否錯(cuò)位或者單位是否統(tǒng)一。4.3 Kmeans 聚類結(jié)果每次跑都不一樣現(xiàn)象同樣的數(shù)據(jù)同樣的 K每次運(yùn)行kmeans得到的聚類中心和概率都不同。原因Kmeans 對初始中心敏感默認(rèn)只跑一次容易陷入局部最優(yōu)。另外如果沒有固定隨機(jī)種子初始化的隨機(jī)性也會(huì)導(dǎo)致結(jié)果波動(dòng)。解決設(shè)置Replicates參數(shù)為 5 到 10讓算法多跑幾次取最優(yōu)。同時(shí)在腳本開頭加rng(42)固定隨機(jī)種子。如果數(shù)據(jù)量特別大可以用kmeans的Start參數(shù)指定初始中心比如用層次聚類的中心作為起點(diǎn)。4.4 典型場景的日出力曲線出現(xiàn)“鋸齒”或不平滑現(xiàn)象聚類還原出來的典型場景曲線在某些時(shí)間點(diǎn)突然跳變看起來不像真實(shí)出力曲線。原因Kmeans 在歸一化空間聚類還原時(shí)用簇內(nèi)均值如果簇內(nèi)樣本的形態(tài)差異較大均值曲線就會(huì)出現(xiàn)不自然的波動(dòng)。另外如果歸一化時(shí)用了按列歸一化還原后時(shí)間維度上的連續(xù)性會(huì)被破壞。解決改用按行歸一化保留每條曲線的時(shí)序形態(tài)。如果鋸齒仍然明顯增加 K 值讓每個(gè)簇更緊湊或者改用 K-medoids 算法它選的是實(shí)際樣本作為中心曲線更自然。4.5 場景概率之和不為 1 或某個(gè)簇概率為 0現(xiàn)象計(jì)算完概率發(fā)現(xiàn)總和不是 1或者某個(gè)典型場景的概率是 0。原因概率為 0 說明該簇沒有分配到任何樣本通常是 K 設(shè)得太大或者初始中心選得不好導(dǎo)致空簇。概率之和不為 1 一般是浮點(diǎn)誤差累積或者計(jì)算時(shí)用了錯(cuò)誤的樣本總數(shù)。解決概率為 0 時(shí)減少 K 值或者用kmeans的EmptyAction參數(shù)設(shè)為singleton讓空簇重新分配。概率之和的誤差在 1e-10 量級以內(nèi)可以忽略如果偏差大檢查idx的長度是否和樣本總數(shù)一致。5. 從典型場景到優(yōu)化模型一個(gè)驗(yàn)證場景質(zhì)量的實(shí)用技巧場景生成完怎么判斷這組典型場景到底能不能用最直接的辦法是回代驗(yàn)證把典型場景和對應(yīng)概率輸入你的優(yōu)化模型看目標(biāo)函數(shù)值和用全量場景或歷史數(shù)據(jù)算出來的差多少。如果偏差在可接受范圍內(nèi)說明場景削減沒有丟失關(guān)鍵信息。但更細(xì)一步的驗(yàn)證是看統(tǒng)計(jì)特征保留度。具體做法是對原始場景集合和典型場景集合分別計(jì)算風(fēng)電出力的均值、方差、以及風(fēng)光之間的 Kendall τ然后對比。典型場景的統(tǒng)計(jì)量應(yīng)該和原始場景接近如果偏差大說明 K 選小了或者聚類特征構(gòu)造有問題。% 原始場景的統(tǒng)計(jì)特征 wind_mean_orig mean(wind_sim(:)); wind_std_orig std(wind_sim(:)); tau_orig corr(U_sim, V_sim, Type, Kendall); % 典型場景的統(tǒng)計(jì)特征加權(quán) wind_mean_typ sum(prob .* mean(wind_typical, 2)); wind_std_typ sqrt(sum(prob .* var(wind_typical, 0, 2))); % 典型場景的 Kendall tau 需要把典型場景展開成加權(quán)樣本后計(jì)算這個(gè)驗(yàn)證步驟花不了幾分鐘但能幫你避免“場景看起來沒問題、一放進(jìn)優(yōu)化模型結(jié)果就偏”的尷尬。我自己的習(xí)慣是每次調(diào)整 K 或 Copula 參數(shù)后先跑這個(gè)驗(yàn)證通過了再往下做。還有一個(gè)實(shí)用技巧是可視化對比。把原始場景集合畫成灰色半透明曲線典型場景畫成粗實(shí)線疊在一起看。如果典型場景能覆蓋原始場景的主要形態(tài)區(qū)域說明削減合理如果某些區(qū)域的原始曲線完全沒有被典型場景代表說明 K 不夠或者聚類特征需要調(diào)整。這個(gè)圖在論文里也很好用審稿人一眼就能看出場景削減的效果。最后說一個(gè)我踩過的坑有一次用 t-Copula 生成場景自由度參數(shù)設(shè)得太小導(dǎo)致尾部樣本過多聚類出來的典型場景里出現(xiàn)了好幾個(gè)極端場景概率還不低。后來把自由度調(diào)到 5 以上尾部就正常了。t-Copula 的自由度參數(shù) ν 控制尾部厚度ν 越小尾部越厚一般取 3 到 10 之間具體看數(shù)據(jù)的尾部特征。如果不確定先用 Gaussian Copula 跑一遍做基準(zhǔn)再換 t-Copula 對比差異明顯再考慮用 t-Copula。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取