性的蒙特卡洛場景生成與削減方法)
最近在研究新能源出力場景建模時有件事讓我印象很深一位做隨機(jī)調(diào)度的同事跑了一版蒙特卡洛Monte CarloMC場景分布檢驗全都通過可一放進(jìn)日前機(jī)組組合模型爬坡約束卻緊張到離譜。排查了半天問題不在求解器也不在目標(biāo)函數(shù)而是他生成的場景根本沒有保留時間上的連續(xù)性。相鄰時段的風(fēng)電出力一會兒滿發(fā)一會兒歸零這種現(xiàn)實中幾乎不可能出現(xiàn)的抖動把爬坡約束徹底帶偏了。這個項目——考慮時序相關(guān)性的MC場景生成與削減研究就是想解決這一類問題不只要讓生成場景的“每個時刻分布”像還要讓“整條軌跡的變化節(jié)奏”像。場景生成與削減本質(zhì)上是隨機(jī)優(yōu)化和不確定性分析里最基礎(chǔ)也最關(guān)鍵的一環(huán)。面向電力系統(tǒng)的新能源出力模擬、概率潮流計算、風(fēng)險評估、隨機(jī)規(guī)劃都離不開“先生成一批有代表性的場景再把它削減到計算能承受的規(guī)?！边@條鏈路。很多人把關(guān)注點全放在概率分布擬合上卻忽略了時序依賴這個維度。而這個項目恰恰把目光聚焦在時序相關(guān)性上用它串起MC場景生成和后端的場景削減是一個很典型、很實用、也很耐琢磨的研究方向。我拿到這類項目時第一反應(yīng)不是直接寫采樣代碼而是先把鏈路拆開要生成什么樣的場景、用什么度量來評價場景、削減后信息損失怎么控制、時序相關(guān)性在哪一步容易被破壞。這篇文章就按這條線路把我實際做過、驗證過、踩過坑的經(jīng)驗整理出來。1. 項目總覽與核心痛點拆解1.1 這個項目到底在解決什么問題先明確一下這里說的“MC”指什么。在電力系統(tǒng)場景分析里MC通常指蒙特卡洛模擬——通過大量隨機(jī)抽樣來逼近真實概率分布的一套方法。它的基本思路很直觀如果想知道未來風(fēng)電出力有哪些可能軌跡就從歷史數(shù)據(jù)里估計出力概率分布然后不斷抽樣生成成百上千條完整軌跡。問題在于常規(guī)MC抽樣有一個先天的毛病它傾向于把每個時間段當(dāng)成獨立隨機(jī)變量來處理。簡單說就是“每個時刻的分布都對但時間軸上的前后關(guān)系全亂了”。這在數(shù)學(xué)上表現(xiàn)為忽略了自相關(guān)性同一序列前后時刻的相關(guān)性和互相關(guān)性不同序列同一時刻的相關(guān)性。這個項目要解決的正是“怎么在MC抽樣階段就把時序相關(guān)性織進(jìn)樣本里然后在削減階段還能保住這條相關(guān)性不散架”。它并不是一個純理論問題而是直接面向工程應(yīng)用的削減后的幾十個場景要給優(yōu)化模型使用如果場景本身的時序特征失真哪怕分布再準(zhǔn)最終的決策結(jié)論也不可靠。我在實際操作中的體會是這類項目適合兩類人深入學(xué)習(xí)做電力系統(tǒng)隨機(jī)優(yōu)化、新能源消納分析的工程師和研究人員需要更真實的風(fēng)光電出力場景作為模型輸入。做場景生成、概率建模的數(shù)據(jù)方向從業(yè)者希望給傳統(tǒng)MC方法補上時序依賴的短板再對比數(shù)據(jù)驅(qū)動類生成模型時有一個扎實的基線。1.2 忽略時序相關(guān)性會產(chǎn)生哪些嚴(yán)重后果不保留時序相關(guān)性的場景集表面上看“分布是對的”但放到具體應(yīng)用里會出現(xiàn)三種很典型的失真。第一種失真出力持續(xù)性和爬坡特性丟失。真實的風(fēng)電出力不會在相鄰時刻從滿發(fā)瞬間跌倒零光伏出力也不會在午后突然跳變到凌晨值。但獨立采樣生成的場景完全可能做到這一點。把這些抖動場景丟進(jìn)機(jī)組組合模型求解器會誤以為系統(tǒng)面臨極其苛刻的爬坡需求于是不得不預(yù)留大量旋轉(zhuǎn)備用、安排更多機(jī)組開機(jī)導(dǎo)致運行成本虛高。第二種失真極端持續(xù)事件缺失。系統(tǒng)最怕的不是“單一時刻出力很低”而是“連續(xù)多天無風(fēng)場導(dǎo)致出力持續(xù)低迷”。這一類低概率但高影響的持續(xù)事件恰恰是儲能容量配置和備用容量決策的主要驅(qū)動力。獨立MC采樣在概率上很難碰出這種“連續(xù)低位”的路徑即使采樣數(shù)量很大也會因為缺少持續(xù)狀態(tài)而漏掉極端風(fēng)險。第三種失真誤導(dǎo)評價指標(biāo)。很多人只用邊際分布比如每個時刻的CDF、均值、方差來評價場景質(zhì)量結(jié)果發(fā)現(xiàn)“擬合得很好”但ACF自相關(guān)函數(shù)完全對不上。這種評價方式等于只檢查了照片的顏色分布卻忽略了照片里的人在做什么動作。項目研究里如果把“時序相關(guān)性”作為一等的評價維度才能避免這類誤導(dǎo)。我見過一個具體的例子某風(fēng)電場景集在邊際分布上通過了KS檢驗但lag 1自相關(guān)系數(shù)只有實況數(shù)據(jù)的三分之二。把它用于兩階段隨機(jī)優(yōu)化備用成本比真實場景結(jié)果高了約12%原因就是場景中出現(xiàn)了過多人為的出力反向跳變把備用的“需求假象”放大了。所以在這類項目里時序相關(guān)性不是錦上添花而是直接決定結(jié)論可靠性的關(guān)鍵維度。2. 時序相關(guān)性MC場景生成的整體思路2.1 蒙特卡洛場景生成的常規(guī)套路和缺陷要理解這個項目的方案選型得先看看目前MC場景生成的幾種常見做法以及它們各自卡在哪里。第一種是歷史數(shù)據(jù)重采樣。直接從歷史出力序列里按日抽樣完整保留一條條真實軌跡。這種做法相關(guān)性天然是對的但樣本數(shù)量受歷史長度限制而且很難覆蓋未發(fā)生過的極端組合。想生成上千個場景往往就力不從心了。第二種是逐時刻獨立抽樣。對每個時刻的歷史數(shù)據(jù)做經(jīng)驗分布或參數(shù)分布擬合然后每個時段獨立采樣拼成一條完整軌跡。這種做法能生成大量樣本但正如前面所說時序相關(guān)性會被徹底打散。如果目標(biāo)場景本身有很強的日內(nèi)規(guī)律比如光伏出力的“中午高、早晚低”獨立采樣生成的路徑會有大量違反物理規(guī)律的抖動。第三種是誤差疊加MC。先用預(yù)測模型得到一條確定性出力曲線然后對預(yù)測誤差做MC采樣把誤差場景疊加到預(yù)測曲線上形成最終場景。這種做法的好處是貼近業(yè)務(wù)因為實際調(diào)度往往是在預(yù)測基礎(chǔ)上考慮不確定性。但很多人忽略了預(yù)測誤差序列本身是存在自相關(guān)的尤其是超短期預(yù)測誤差在相鄰時段高度相關(guān)。如果對誤差做獨立采樣同樣會破壞場景的時間連續(xù)性。所以整套方案的設(shè)計關(guān)鍵在于不能在采樣階段把相關(guān)性丟掉然后再試圖補救——應(yīng)該在一開始就把相關(guān)結(jié)構(gòu)作為目標(biāo)嵌入到采樣路徑中去。2.2 本項目的解法把相關(guān)性信息嵌入采樣路徑這個項目的核心模型思路可以概括成一個三步鏈路從歷史數(shù)據(jù)里估計目標(biāo)相關(guān)結(jié)構(gòu)也就是自相關(guān)矩陣和互相關(guān)矩陣。在標(biāo)準(zhǔn)正態(tài)空間里生成具有該相關(guān)結(jié)構(gòu)的樣本。把樣本通過逆變換映射回目標(biāo)分布空間得到最終的時序場景。這里出現(xiàn)了一個常見術(shù)語叫“相關(guān)矩陣”。你可以把它理解為一張表格第i行第j列的數(shù)字表示“第i個時段和第j個時段一起變化的傾向”。數(shù)值接近1時表示兩個時段趨向于同漲同跌接近0時表示它們互不影響。風(fēng)電出力的自相關(guān)矩陣通常呈現(xiàn)“近處高、遠(yuǎn)處低”的帶狀結(jié)構(gòu)這也是原始序列“有記憶性”的體現(xiàn)。在具體實現(xiàn)上我比較常用的兩條路分別是Cholesky分解法多元正態(tài)法先生成獨立標(biāo)準(zhǔn)正態(tài)樣本再用相關(guān)矩陣的Cholesky分解結(jié)果做線性變換讓樣本帶上目標(biāo)相關(guān)。然后對每個時段用標(biāo)準(zhǔn)正態(tài)CDF的逆函數(shù)轉(zhuǎn)到均勻空間再用目標(biāo)分布CDF的逆函數(shù)轉(zhuǎn)到真實出力空間。這套流程在數(shù)學(xué)上被稱為Nataf變換是MC相關(guān)性抽樣的經(jīng)典做法。Iman-Conover排序法先生成一個“模板樣本”讓其相關(guān)結(jié)構(gòu)接近目標(biāo)然后對原始獨立樣本的每一列按模板對應(yīng)列的秩順序重排使重排后的樣本具備模板的相關(guān)結(jié)構(gòu)。這個方法完全不要求原始樣本服從正態(tài)分布所以對風(fēng)電出力這種“非高斯、有邊界、大量零點”的數(shù)據(jù)特別友好。補充說明一下我這里提的兩種方法都是基于常見實踐的可靠方案。實際項目里不必迷信某一種關(guān)鍵是理解相關(guān)性在哪一步生成、在哪一步被壓縮。2.3 適合哪些數(shù)據(jù)和應(yīng)用場景這個方案最適用的場景是那些“時序依賴明顯、目標(biāo)分布非高斯、規(guī)模適中”的變量。典型代表就是風(fēng)電出力。風(fēng)速本身具有強自相關(guān)和日周期/季節(jié)周期風(fēng)機(jī)功率曲線又把風(fēng)速映射成帶有大量零點、上下截斷的出力序列分布形態(tài)高度非高斯直接用多元正態(tài)假設(shè)建模會誤差很大。這正是排序法配合經(jīng)驗分布最能發(fā)揮作用的地方。光伏出力也類似但有更強的確定性日內(nèi)模式。對光伏場景建模時我一般建議先除以理論晴空出力做歸一化再把歸一化后的“天空晴朗指數(shù)”序列作為建模對象這樣自相關(guān)結(jié)構(gòu)更容易提取場景也更穩(wěn)定。負(fù)荷序列、電價序列、徑流序列也都可以用這套框架。因為它們的共同特征是邊際分布偏斜、時序關(guān)聯(lián)強、同時有尖峰和持續(xù)狀態(tài)。需要提醒的是這套方法并不是萬能的。如果研究對象的變量之間不存在明顯時間依賴比如設(shè)備隨機(jī)的故障狀態(tài)、相互獨立的事件開關(guān)那就不需要費勁做時序相關(guān)性建模強行引入相關(guān)結(jié)構(gòu)反而會增加模型復(fù)雜度和誤差來源。3. 落地實操一時序相關(guān)性MC場景生成3.1 數(shù)據(jù)準(zhǔn)備與統(tǒng)計建模步驟在開始采樣之前數(shù)據(jù)準(zhǔn)備工作往往決定了整個項目能走多遠(yuǎn)。我習(xí)慣按下面這幾步做數(shù)據(jù)準(zhǔn)備第一步數(shù)據(jù)清洗。風(fēng)電、光伏出力數(shù)據(jù)里經(jīng)常出現(xiàn)異常尖峰、缺測、連續(xù)平臺值。這一步我會做三件事剔除超出物理范圍的記錄比如出力大于裝機(jī)容量、光伏夜間出力非零對缺測值做線性插值或前一天同時刻替換連續(xù)重復(fù)值過多的時段要格外警惕可能是數(shù)據(jù)采集故障需要用相鄰日數(shù)據(jù)修正。第二步確定建模對象。我的實踐經(jīng)驗是優(yōu)先對“預(yù)測誤差”建模而不是直接對出力絕對值建模。原因很簡單實際調(diào)度中預(yù)測曲線是已知信息需要模擬的是未來圍繞預(yù)測的波動。誤差建模時場景生成公式就是“確定性預(yù)測 MC誤差場景”不僅物理含義清晰也方便后期做場景削減和隨機(jī)優(yōu)化的銜接。第三步估計邊際分布。這里不要想當(dāng)然地套正態(tài)分布。風(fēng)電出力的邊際分布往往是“零出力點的概率集中 高出力區(qū)間厚尾”的組合形態(tài)光伏出力在歸一化后也明顯偏斜。比較穩(wěn)妥的辦法是直接采用經(jīng)驗分布把歷史數(shù)據(jù)按數(shù)值排序用排序位置除以樣本數(shù)作為累積概率估計再用分段線性插值或核密度估計平滑一下。這樣不需要做任何參數(shù)分布假設(shè)模型風(fēng)險低后續(xù)逆變換采樣也簡單。第四步估計相關(guān)結(jié)構(gòu)。對按“天”組織的數(shù)據(jù)把每天的數(shù)據(jù)作為一個樣本統(tǒng)一減去當(dāng)天均值并除以標(biāo)準(zhǔn)差得到標(biāo)準(zhǔn)化矩陣。然后用np.corrcoef計算時段間的相關(guān)矩陣。這里要注意相關(guān)矩陣的估計對異常值很敏感所以第一步的數(shù)據(jù)清洗不能偷懶。下面給出一個可以直接上手的Python代碼骨架import numpy as np import pandas as pd # data_shape: (days, T)T為每日時段數(shù)day為歷史天數(shù) daily_std (data - data.mean(axis0)) / data.std(axis0) # 時段間Pearson相關(guān)矩陣形狀 (T, T) R np.corrcoef(daily_std.T) # 自相關(guān)結(jié)構(gòu)檢查看lag 1的相關(guān)系數(shù) lag1_acf np.mean([np.corrcoef(daily_std[:, t], daily_std[:, t1])[0, 1] for t in range(daily_std.shape[1] - 1)]) print(lag1自相關(guān)均值: , lag1_acf)T通常取24小時級或9615分鐘級。相關(guān)矩陣規(guī)模分別是24×24和96×96都不算大計算成本可以忽略。3.2 Cholesky / 排序法生成相關(guān)樣本的細(xì)節(jié)對比數(shù)據(jù)準(zhǔn)備完成后就進(jìn)入核心環(huán)節(jié)生成帶相關(guān)性的采樣時序。我以誤差建模為例把兩條技術(shù)路線展開說清楚。路線一Cholesky分解法的完整流程如下對誤差歷史序列做經(jīng)驗分布擬合得到每個時段的邊際分布F_t。生成一個N×T的獨立標(biāo)準(zhǔn)正態(tài)矩陣ZN是需要生成的場景數(shù)。對相關(guān)矩陣R做Cholesky分解得到下三角矩陣L滿足R L L^T。計算Y Z L^T此時Y的相關(guān)矩陣接近R。對Y的每個元素做Phi(Y_t)得到均勻分布樣本再代入F_t^{-1}得到誤差場景。代碼段長這樣from scipy.stats import norm # 1) 估計誤差邊際分布這里用empirical CDF 線性插值 # 為簡化演示假設(shè)每時段誤差服從正態(tài)分布并存儲 mean_t, std_t # 實際項目中建議use經(jīng)驗分布代碼量略長但更穩(wěn) L np.linalg.cholesky(R 1e-6 * np.eye(R.shape[0])) # 加ridge防非正定 Z np.random.randn(N, T) Y Z L.T U norm.cdf(Y) # (N, T) err_scenes np.array([F_inv_t(u_t) for u_t in U.T]).T路線二Iman-Conover排序法邏輯不太一樣生成一個模板矩陣X_tmpl直接用多元正態(tài)分布協(xié)方差矩陣設(shè)置成目標(biāo)R。這一步可以用Cholesky實現(xiàn)但即使模板相關(guān)矩陣略有偏差后續(xù)排序也會校正。對模板每一列計算“秩信息”也就是按數(shù)值大小排名的位置索引。對誤差獨立樣本矩陣S的每一列按照模板對應(yīng)列的秩順序重新排列。排序完成后S的秩相關(guān)結(jié)構(gòu)就與模板一致了因此相關(guān)性也被傳遞到了目標(biāo)分布空間。這段邏輯最巧妙的地方是它只利用秩信息不需要計算任何逆變換對任何邊際分布都適用。我在實際項目中更偏愛Iman-Conover排序法原因后面會詳細(xì)展開。但在實現(xiàn)時有個細(xì)節(jié)模板矩陣必須固定種子否則每次排序后的樣本雖然相關(guān)結(jié)構(gòu)一致卻會引入額外的抽樣噪聲影響復(fù)現(xiàn)性。兩種路線的對比如下對比維度Cholesky分解法Iman-Conover排序法對邊際分布要求需要已知并計算逆CDF無要求經(jīng)驗分布即可相關(guān)矩陣要求必須正定允許適當(dāng)非正定非線性相關(guān)性壓縮問題存在需要校正較小排序保留了秩相關(guān)實現(xiàn)復(fù)雜度低幾行代碼中等需要維護(hù)排序索引適用數(shù)據(jù)近似正態(tài)、高斯耦合風(fēng)電、光伏等非高斯強偏數(shù)據(jù)3.3 從標(biāo)準(zhǔn)化空間還原到物理場景空間生成誤差場景后最后一步是還原到物理出力空間。還原本身不復(fù)雜但有三個細(xì)節(jié)容易被忽略。細(xì)節(jié)一物理邊界約束。風(fēng)電出力必須落在[0, 裝機(jī)容量]之間光伏出力在夜間必須為0。逆變換采樣后要對場景做截斷處理。我建議不要簡單硬截斷而是把超出上限的部分映射回上限同時把低于下限的部分映射回0。硬截斷會改變邊際分布的尾部形狀但在工程上通常影響可控所以先截斷再在評價階段復(fù)核分布即可。細(xì)節(jié)二場景形態(tài)平滑性。還原后的場景可能帶有高頻抖動。有人喜歡加一個滑動平均平滑我要提醒平滑濾波會改變時序相關(guān)結(jié)構(gòu)尤其是lag 1、lag 2的相關(guān)系數(shù)會被顯著拉低。如果后續(xù)要用ACF來評價場景質(zhì)量平滑必須要謹(jǐn)慎或者干脆不做。細(xì)節(jié)三生成場景的最終物理公式。如果采用誤差建模最終場景為[ P_{scene}(t) P_{forecast}(t) err_{scene}(t) ]隨后執(zhí)行截斷和必要修正。如果直接對出力建模則直接用逆變換采樣結(jié)果作為場景。另外我強烈建議在生成完場景后做一次“回檢”計算生成場景集合的ACF矩陣、邊際分布直方圖和目標(biāo)應(yīng)用指標(biāo)把這些與真實歷史數(shù)據(jù)放在一起對比。這一步雖然簡單卻能在進(jìn)入削減環(huán)節(jié)之前發(fā)現(xiàn)建模階段的系統(tǒng)性問題省去后續(xù)返工。4. 落地實操二場景削減與權(quán)重重構(gòu)4.1 場景削減的數(shù)學(xué)目標(biāo)和實用取舍MC生成一兩千個場景并不難難的是把這些場景放進(jìn)優(yōu)化模型。每個場景對應(yīng)一組約束和變量場景數(shù)量上升模型規(guī)模會成倍膨脹。削減的目的就是在盡量保留分布信息的前提下把場景集壓縮到計算可承受的規(guī)模。場景削減的數(shù)學(xué)目標(biāo)核心是控制“概率分布距離”。業(yè)界常用的距離是Kantorovich距離它在離散場景集上的形式非常直觀從原始場景集里刪掉某些場景并把它們的概率權(quán)重加到剩余場景上使削減后集合與原始集合之間“搬運概率質(zhì)量的距離”最小。這個目標(biāo)可以用一個生活例子來理解你有一個班的成績單大量場景想挑幾個代表參加座談削減后場景既要讓沒來的人覺得“代表表達(dá)了我的情況”又不能讓代表數(shù)量太多計算量。Kantorovich距離就是在量化“代表們能不能最好地覆蓋全班的各種情況”。在實際項目里我建議在削減之前先想清楚應(yīng)用場景因為這會影響距離度量里的權(quán)重。如果場景最終用于日前機(jī)組組合那么峰谷爬坡時段和晚高峰時段的重要性要更高如果用于儲能容量配置則連續(xù)多日低出力場景需要特殊保護(hù)不能被輕易“合并”掉。4.2 快速前向削減與同步回代削減的實現(xiàn)要點場景削減算法有兩大經(jīng)典路線快速前向選擇fast forward selection和同步回代消除simultaneous backward reduction?？焖偾跋蜻x擇的核心邏輯是“從空集開始反復(fù)挑選最有代表性的場景加入保留集”。每一步都找這樣一個場景把它加入保留集之后剩余場景到保留集的距離總和最小。重復(fù)M次得到M個保留場景。優(yōu)點是可以控制保留數(shù)量缺點是每一步都要遍歷所有剩余場景計算量偏大。同步回代消除的邏輯則相反從全集開始每次找到距離最近的兩個場景把其中一個刪除把它的概率權(quán)重加到另一個上面。不斷重復(fù)直到保留到目標(biāo)數(shù)量。它的計算量通常更小但有一個隱患如果兩個最近的場景都來自同一個區(qū)域合并會把極端但不常出現(xiàn)的場景提前吃掉。我在實際項目中常用“先粗分桶再精細(xì)削減”的兩階段策略。先用k-means或者k-medoids把2000個場景粗分成50個左右簇然后對每簇內(nèi)部用快速前向選擇選出最終的代表場景。這樣既控制了計算量又保證了代表場景的覆蓋度。注意這里的k-means聚類只是為了縮減搜索空間最終場景距離計算還是用原始時序距離度來完成不直接用聚類中心當(dāng)場景因為聚類中心會破壞物理時序關(guān)系。下面給一個簡化版的快速前向選擇實現(xiàn)骨架def forward_selection(scenes, weights, M): # scenes: (N, T) 原始場景集 # weights: (N,) 每個場景初始權(quán)重如1/N N, T scenes.shape selected [] remaining set(range(N)) for _ in range(M): best_j None best_cost np.inf for j in remaining: # 計算若選擇j作為保留場景剩余點到j(luò)的加權(quán)距離之和 dists np.linalg.norm(scenes[list(remaining)] - scenes[j], axis1) cost np.sum(weights[list(remaining)] * dists) if cost best_cost: best_cost cost best_j j selected.append(best_j) remaining.remove(best_j) # 權(quán)重重構(gòu)每個保留場景吸收剩余場景中距離它最近的那些權(quán)重 new_weights np.zeros(M) for i in range(N): if i in selected: idx selected.index(i) new_weights[idx] weights[i] else: dists [np.linalg.norm(scenes[i] - scenes[j]) for j in selected] idx int(np.argmin(dists)) new_weights[idx] weights[i] return scenes[selected], new_weights代碼里最核心的循環(huán)可以進(jìn)一步用向量化和堆結(jié)構(gòu)優(yōu)化但對幾千個場景、T在24~96的項目直接這樣寫也完全夠用。同步回代消除的實現(xiàn)思路更短維護(hù)場景集合和權(quán)重每次計算兩兩距離找到最小距離對將一個場景的權(quán)重合并給另一個然后刪除該場景。代碼我就不展開了和上面的循環(huán)結(jié)構(gòu)類似只是合并邏輯相反。4.3 削減后場景的時序相關(guān)性驗證削減完成不等于工作結(jié)束。我見過太多項目在削減環(huán)節(jié)把好不容易生成的時序相關(guān)性又給丟了原因是削減算法里用的距離度量沒有照顧到時序結(jié)構(gòu)。如果你直接用歐氏距離衡量兩個場景的差異那么兩個“整體數(shù)值接近但前后順序打亂”的場景會被認(rèn)為很相似。這會導(dǎo)致削減后的場景集把“時間反轉(zhuǎn)”、“相位偏移”的場景合并到同一個代表場景里自相關(guān)結(jié)構(gòu)自然受損。我在某風(fēng)電場項目中做過對比用純歐氏距離把2000個場景削減到20個削減后場景集合的lag 1自相關(guān)系數(shù)從0.83掉到0.51時序特征嚴(yán)重退化。改用加權(quán)歐氏距離在距離計算中額外加入相鄰時段差分的懲罰項后削減后場景集的lag 1自相關(guān)系數(shù)恢復(fù)到0.76已經(jīng)能比較好地保留時序結(jié)構(gòu)。加懲罰項的距離度量可以寫成[ d(\xi_i, \xi_j) \sum_{t1}^{T} w_t \cdot (\xi_{i,t} - \xi_{j,t})^2 \lambda \sum_{t1}^{T-1} ((\xi_{i,t1}-\xi_{i,t}) - (\xi_{j,t1}-\xi_{j,t}))^2 ]其中第一項衡量數(shù)值接近程度第二項衡量相鄰時段變化率近似爬坡的接近程度(\lambda)是權(quán)重系數(shù)。我在實踐中取(\lambda1.5 \sim 2)效果較好。如果你覺得這個加權(quán)不好調(diào)也可以直接用動態(tài)時間彎曲DTW距離它對相位偏移和局部時間變形更魯棒但計算量更大場景多時可以先粗聚類再精細(xì)化。5. 評價指標(biāo)體系與結(jié)果展示框架5.1 分布精度指標(biāo)評價場景集質(zhì)量的第一個維度是它到底還原了多少概率分布信息。這里至少要看三個指標(biāo)一是邊際分布的KS檢驗。對每個時段把生成場景在該時段的CDF與歷史數(shù)據(jù)的CDF做比較取最大差異的p值。如果大部分時段的p值大于0.05說明邊際分布擬合可以接受。不能只報一個“平均p值”要給出最差時段的p值因為最差時段往往對應(yīng)傍晚或清晨的快速爬坡時段恰恰是優(yōu)化模型最敏感的地方。二是1-Wasserstein距離。這個指標(biāo)可以直接用具體數(shù)值來衡量削減帶來的分布漂移。通常做法是畫出“保留場景數(shù)M vs 1-Wasserstein距離”的衰減曲線。當(dāng)M從2增加到5、10、15時距離快速下降隨后進(jìn)入平臺期。曲線的拐點就是“投入更多場景數(shù)換來更少信息損失”的臨界點直接幫助確定合理的削減數(shù)量。三是極端分位數(shù)對比。比如看場景集在1%分位數(shù)和99%分位數(shù)上是否與歷史數(shù)據(jù)一致。風(fēng)電場景的極端低出力分位數(shù)尤其重要因為它直接決定備用容量和儲能配置。僅看均值、方差完全不夠分位數(shù)對尾部信息的刻畫更細(xì)致。我個人的經(jīng)驗是這三個指標(biāo)都要做但不必追求每一項都達(dá)到極其精確重點是讓它們進(jìn)入一個“可接受區(qū)間”保證場景集在統(tǒng)計意義上站得住。5.2 時序結(jié)構(gòu)保留度指標(biāo)評價時序相關(guān)性的直接指標(biāo)是ACF曲線對比。把歷史數(shù)據(jù)的自相關(guān)系數(shù)序列l(wèi)ag 1lag 24和生成場景集合的平均ACF畫在同一張圖上觀察兩者的偏差。更嚴(yán)格一點可以計算ACF曲線的平均絕對誤差MAE。另一個特別能反映場景質(zhì)量的指標(biāo)是“相鄰時段差值分布”。先計算每條場景的差分序列(\Delta_t P_{t1} - P_t)然后比較生成場景集合與歷史數(shù)據(jù)的(\Delta_t)分布。這個指標(biāo)直接抓住“爬坡分布”的形態(tài)如果時序相關(guān)性丟失差分分布的尾部往往會明顯變寬。再一個指標(biāo)是秩相關(guān)矩陣差異。用Spearman秩相關(guān)而不是Pearson相關(guān)因為秩相關(guān)對異常值更穩(wěn)健。比較削減前后場景集合與歷史數(shù)據(jù)的秩相關(guān)矩陣用Frobenius范數(shù)刻畫整體偏差。這個范數(shù)就是一個把矩陣所有元素差異平方再開方的綜合數(shù)值越小說明相關(guān)結(jié)構(gòu)保留得越好。我在項目匯報時習(xí)慣把ACF圖和差分分布圖放在一起展示效果比單純報一個數(shù)值直觀得多也讓評審的人一眼明白“時序相關(guān)性”到底體現(xiàn)在哪里。5.3 面向應(yīng)用的誤差評估統(tǒng)計指標(biāo)再好最終都要回歸到應(yīng)用端驗證。我的做法是把場景集塞進(jìn)一個簡化版兩階段隨機(jī)規(guī)劃模型比如簡化的機(jī)組組合或儲能優(yōu)化比較兩種情形下的結(jié)果偏差情形一使用削減前的全部場景作為基準(zhǔn)理想值。情形二使用削減后的場景集。對比目標(biāo)函數(shù)值的相對偏差。如果削減后成本比全場景結(jié)果明顯偏低說明丟失了某些“惡劣場景”導(dǎo)致優(yōu)化結(jié)果過于樂觀。如果明顯偏高說明保留了太多冗余場景削減效率不足。除此之外還有一個容易被忽略的點約束滿足情況。削減后場景集中如果恰好沒有覆蓋某個關(guān)鍵極端場景那么一些安全約束可能在削減后的模型里永遠(yuǎn)不會被觸發(fā)得到的決策方案在真實世界中會過于冒進(jìn)。針對這個問題我建議在削減時明確保留“歷史極端案例”也就是把歷史數(shù)據(jù)中出現(xiàn)過的連續(xù)低出力、持續(xù)高出力、劇烈爬坡等特殊時段序列作為“種子場景”強制加入保留集再在此基礎(chǔ)上做削減。6. 常見問題與排查技巧實錄6.1 相關(guān)矩陣非正定的頑固問題MC場景生成環(huán)節(jié)最常見的問題是Cholesky分解報錯提示相關(guān)矩陣不是正定矩陣。這通常發(fā)生在T比較大比如96且樣本數(shù)有限時。原因有兩個一是長序列自相關(guān)矩陣本質(zhì)上是Toeplitz結(jié)構(gòu)其特征值在短樣本估計下很容易出現(xiàn)負(fù)值二是數(shù)據(jù)清洗不干凈時異常值會扭曲相關(guān)矩陣估計。排查和解決辦法按順序試給相關(guān)矩陣對角線加一個小的正則項例如(R \epsilon I)(\epsilon)取1e-61e-4。這相當(dāng)于給每個時段的方差增加一點“數(shù)值穩(wěn)定墊”如果加了正則還是不行計算最接近R的半正定矩陣用特征值分解把負(fù)特征值置零再歸一化直接放棄Cholesky分解改用Iman-Conover排序法。排序法不要求相關(guān)矩陣嚴(yán)格正定實踐中魯棒性最好。我在項目里用了排序法之后就再也沒被這個錯誤卡住過。但如果論文或報告中要求可復(fù)現(xiàn)的數(shù)值實驗則仍需要記錄清楚相關(guān)矩陣估計的具體方式和正則參數(shù)。6.2 削減后場景特征漂移怎么辦如果你發(fā)現(xiàn)削減后的場景集邊際分布看著還行但ACF明顯退化問題大概率出在削減階段的距離度量上。前文提到純歐氏距離不感知時序結(jié)構(gòu)會把時間形態(tài)不同的場景合并起來。解決辦法有兩個方向在距離里加入差分懲罰項邏輯見4.3節(jié)。這是成本最低的改動。改用“代表性場景按聚類中心搜索”的思路先聚類再在每個類內(nèi)搜索離類中心最近的原始場景作為代表避免直接拿類中心當(dāng)場景。另外還要檢查一下削減后場景權(quán)重是否出現(xiàn)極端值。如果某個保留場景的權(quán)重超過0.25說明它“吸收”了過多場景很可能抹平了多種不同的時序模式。這種情況下建議強制提高保留場景數(shù)量或者在削減前做“極端場景保護(hù)”。6.3 場景數(shù)量選擇糾結(jié)怎么破很多初學(xué)者會問到底削減到多少個場景合適。這個問題很難給統(tǒng)一答案但我建議用“Wasserstein衰減曲線 項目計算預(yù)算”雙標(biāo)準(zhǔn)來定。具體操作畫出M從2到30的1-Wasserstein距離曲線找曲線拐點。如果拐點在M15附近那就設(shè)在15如果計算資源緊張比如要做上千次隨機(jī)優(yōu)化再往下壓到510同時確認(rèn)應(yīng)用端誤差在可接受范圍內(nèi)。還有一個經(jīng)驗是不要把M設(shè)得太小。當(dāng)M低于10時場景集對極端事件的分辨率會快速下降備用配置結(jié)果容易出現(xiàn)系統(tǒng)性偏差。我做過一次對比M從8增加到12兩階段隨機(jī)規(guī)劃的目標(biāo)函數(shù)值變化量達(dá)到3%以內(nèi)而M繼續(xù)增加到20時變化已經(jīng)小于0.5%這說明M12左右就已經(jīng)合理。6.4 生成場景“看起來不夠隨機(jī)”怎么理解有時生成出來的場景集會在某個時段內(nèi)“擠成一團(tuán)”看起來多樣性不足。我遇到過幾次原因不是采樣失敗而是相關(guān)矩陣估計過于強烈。比如某風(fēng)電場出力數(shù)據(jù)有很強的日內(nèi)同漲同跌模式自相關(guān)矩陣的所有元素都接近0.9那么采樣出的場景當(dāng)然會高度相似。如果你希望場景集保留更多多樣性可以在相關(guān)矩陣上做一個“收縮”處理比如把相關(guān)矩陣乘以一個衰減系數(shù)(\rho0.85)讓相關(guān)結(jié)構(gòu)溫和化。但注意這種做法會犧牲一點精度所以我更建議在評價時先確認(rèn)“多樣性不足”是否真的影響應(yīng)用結(jié)果而不是憑感覺調(diào)參數(shù)。7. 寫在后面的一些個人經(jīng)驗這個項目最讓我受益的一點是它逼著我把“不確定性模擬”從抽象概念落到了具體工程對象上。以前我做MC抽樣關(guān)注的只是“樣本夠不夠多、分布對不對”很少去管這些樣本放到?jīng)Q策模型里會產(chǎn)生什么后果。直到被爬坡約束的真實案例打臉我才真正意識到時序相關(guān)性不是統(tǒng)計潔癖而是決策模型本身對場景軌跡提出的硬需求。如果現(xiàn)在有人讓我從零做這類研究我會建議走一條最穩(wěn)的閉環(huán)路線先拿歷史數(shù)據(jù)估計經(jīng)驗分布和相關(guān)矩陣用Iman-Conover排序法生成2000個場景再用帶差分懲罰的加權(quán)歐氏距離做快速前向削減最后用ACF對比和兩階段優(yōu)化偏差來驗證效果。這條鏈路每一步都可復(fù)現(xiàn)、可解釋也容易排查問題。等這條鏈路跑通之后再去嘗試Copula函數(shù)、變分自編碼器、生成對抗網(wǎng)絡(luò)這些更花哨的生成手段就有了一個扎實的對比基線。最后再分享一個實踐中額外收獲的小技巧在做場景削減時永遠(yuǎn)保留至少一條“歷史最嚴(yán)苛場景”。不管削減算法怎么選把這個場景強制加入保留集能在后續(xù)的隨機(jī)優(yōu)化和安全校核中省掉很多解釋成本和返工時間。這個方法簡單到很多人不屑于用但恰恰是我在多個項目里避免“樂觀偏差”最有效的土辦法。 考慮時序相關(guān)性的蒙特卡洛場景生成與削減方法
最近在研究新能源出力場景建模時有件事讓我印象很深一位做隨機(jī)調(diào)度的同事跑了一版蒙特卡洛Monte CarloMC場景分布檢驗全都通過可一放進(jìn)日前機(jī)組組合模型爬坡約束卻緊張到離譜。排查了半天問題不在求解器也不在目標(biāo)函數(shù)而是他生成的場景根本沒有保留時間上的連續(xù)性。相鄰時段的風(fēng)電出力一會兒滿發(fā)一會兒歸零這種現(xiàn)實中幾乎不可能出現(xiàn)的抖動把爬坡約束徹底帶偏了。這個項目——考慮時序相關(guān)性的MC場景生成與削減研究就是想解決這一類問題不只要讓生成場景的“每個時刻分布”像還要讓“整條軌跡的變化節(jié)奏”像。場景生成與削減本質(zhì)上是隨機(jī)優(yōu)化和不確定性分析里最基礎(chǔ)也最關(guān)鍵的一環(huán)。面向電力系統(tǒng)的新能源出力模擬、概率潮流計算、風(fēng)險評估、隨機(jī)規(guī)劃都離不開“先生成一批有代表性的場景再把它削減到計算能承受的規(guī)?！边@條鏈路。很多人把關(guān)注點全放在概率分布擬合上卻忽略了時序依賴這個維度。而這個項目恰恰把目光聚焦在時序相關(guān)性上用它串起MC場景生成和后端的場景削減是一個很典型、很實用、也很耐琢磨的研究方向。我拿到這類項目時第一反應(yīng)不是直接寫采樣代碼而是先把鏈路拆開要生成什么樣的場景、用什么度量來評價場景、削減后信息損失怎么控制、時序相關(guān)性在哪一步容易被破壞。這篇文章就按這條線路把我實際做過、驗證過、踩過坑的經(jīng)驗整理出來。1. 項目總覽與核心痛點拆解1.1 這個項目到底在解決什么問題先明確一下這里說的“MC”指什么。在電力系統(tǒng)場景分析里MC通常指蒙特卡洛模擬——通過大量隨機(jī)抽樣來逼近真實概率分布的一套方法。它的基本思路很直觀如果想知道未來風(fēng)電出力有哪些可能軌跡就從歷史數(shù)據(jù)里估計出力概率分布然后不斷抽樣生成成百上千條完整軌跡。問題在于常規(guī)MC抽樣有一個先天的毛病它傾向于把每個時間段當(dāng)成獨立隨機(jī)變量來處理。簡單說就是“每個時刻的分布都對但時間軸上的前后關(guān)系全亂了”。這在數(shù)學(xué)上表現(xiàn)為忽略了自相關(guān)性同一序列前后時刻的相關(guān)性和互相關(guān)性不同序列同一時刻的相關(guān)性。這個項目要解決的正是“怎么在MC抽樣階段就把時序相關(guān)性織進(jìn)樣本里然后在削減階段還能保住這條相關(guān)性不散架”。它并不是一個純理論問題而是直接面向工程應(yīng)用的削減后的幾十個場景要給優(yōu)化模型使用如果場景本身的時序特征失真哪怕分布再準(zhǔn)最終的決策結(jié)論也不可靠。我在實際操作中的體會是這類項目適合兩類人深入學(xué)習(xí)做電力系統(tǒng)隨機(jī)優(yōu)化、新能源消納分析的工程師和研究人員需要更真實的風(fēng)光電出力場景作為模型輸入。做場景生成、概率建模的數(shù)據(jù)方向從業(yè)者希望給傳統(tǒng)MC方法補上時序依賴的短板再對比數(shù)據(jù)驅(qū)動類生成模型時有一個扎實的基線。1.2 忽略時序相關(guān)性會產(chǎn)生哪些嚴(yán)重后果不保留時序相關(guān)性的場景集表面上看“分布是對的”但放到具體應(yīng)用里會出現(xiàn)三種很典型的失真。第一種失真出力持續(xù)性和爬坡特性丟失。真實的風(fēng)電出力不會在相鄰時刻從滿發(fā)瞬間跌倒零光伏出力也不會在午后突然跳變到凌晨值。但獨立采樣生成的場景完全可能做到這一點。把這些抖動場景丟進(jìn)機(jī)組組合模型求解器會誤以為系統(tǒng)面臨極其苛刻的爬坡需求于是不得不預(yù)留大量旋轉(zhuǎn)備用、安排更多機(jī)組開機(jī)導(dǎo)致運行成本虛高。第二種失真極端持續(xù)事件缺失。系統(tǒng)最怕的不是“單一時刻出力很低”而是“連續(xù)多天無風(fēng)場導(dǎo)致出力持續(xù)低迷”。這一類低概率但高影響的持續(xù)事件恰恰是儲能容量配置和備用容量決策的主要驅(qū)動力。獨立MC采樣在概率上很難碰出這種“連續(xù)低位”的路徑即使采樣數(shù)量很大也會因為缺少持續(xù)狀態(tài)而漏掉極端風(fēng)險。第三種失真誤導(dǎo)評價指標(biāo)。很多人只用邊際分布比如每個時刻的CDF、均值、方差來評價場景質(zhì)量結(jié)果發(fā)現(xiàn)“擬合得很好”但ACF自相關(guān)函數(shù)完全對不上。這種評價方式等于只檢查了照片的顏色分布卻忽略了照片里的人在做什么動作。項目研究里如果把“時序相關(guān)性”作為一等的評價維度才能避免這類誤導(dǎo)。我見過一個具體的例子某風(fēng)電場景集在邊際分布上通過了KS檢驗但lag 1自相關(guān)系數(shù)只有實況數(shù)據(jù)的三分之二。把它用于兩階段隨機(jī)優(yōu)化備用成本比真實場景結(jié)果高了約12%原因就是場景中出現(xiàn)了過多人為的出力反向跳變把備用的“需求假象”放大了。所以在這類項目里時序相關(guān)性不是錦上添花而是直接決定結(jié)論可靠性的關(guān)鍵維度。2. 時序相關(guān)性MC場景生成的整體思路2.1 蒙特卡洛場景生成的常規(guī)套路和缺陷要理解這個項目的方案選型得先看看目前MC場景生成的幾種常見做法以及它們各自卡在哪里。第一種是歷史數(shù)據(jù)重采樣。直接從歷史出力序列里按日抽樣完整保留一條條真實軌跡。這種做法相關(guān)性天然是對的但樣本數(shù)量受歷史長度限制而且很難覆蓋未發(fā)生過的極端組合。想生成上千個場景往往就力不從心了。第二種是逐時刻獨立抽樣。對每個時刻的歷史數(shù)據(jù)做經(jīng)驗分布或參數(shù)分布擬合然后每個時段獨立采樣拼成一條完整軌跡。這種做法能生成大量樣本但正如前面所說時序相關(guān)性會被徹底打散。如果目標(biāo)場景本身有很強的日內(nèi)規(guī)律比如光伏出力的“中午高、早晚低”獨立采樣生成的路徑會有大量違反物理規(guī)律的抖動。第三種是誤差疊加MC。先用預(yù)測模型得到一條確定性出力曲線然后對預(yù)測誤差做MC采樣把誤差場景疊加到預(yù)測曲線上形成最終場景。這種做法的好處是貼近業(yè)務(wù)因為實際調(diào)度往往是在預(yù)測基礎(chǔ)上考慮不確定性。但很多人忽略了預(yù)測誤差序列本身是存在自相關(guān)的尤其是超短期預(yù)測誤差在相鄰時段高度相關(guān)。如果對誤差做獨立采樣同樣會破壞場景的時間連續(xù)性。所以整套方案的設(shè)計關(guān)鍵在于不能在采樣階段把相關(guān)性丟掉然后再試圖補救——應(yīng)該在一開始就把相關(guān)結(jié)構(gòu)作為目標(biāo)嵌入到采樣路徑中去。2.2 本項目的解法把相關(guān)性信息嵌入采樣路徑這個項目的核心模型思路可以概括成一個三步鏈路從歷史數(shù)據(jù)里估計目標(biāo)相關(guān)結(jié)構(gòu)也就是自相關(guān)矩陣和互相關(guān)矩陣。在標(biāo)準(zhǔn)正態(tài)空間里生成具有該相關(guān)結(jié)構(gòu)的樣本。把樣本通過逆變換映射回目標(biāo)分布空間得到最終的時序場景。這里出現(xiàn)了一個常見術(shù)語叫“相關(guān)矩陣”。你可以把它理解為一張表格第i行第j列的數(shù)字表示“第i個時段和第j個時段一起變化的傾向”。數(shù)值接近1時表示兩個時段趨向于同漲同跌接近0時表示它們互不影響。風(fēng)電出力的自相關(guān)矩陣通常呈現(xiàn)“近處高、遠(yuǎn)處低”的帶狀結(jié)構(gòu)這也是原始序列“有記憶性”的體現(xiàn)。在具體實現(xiàn)上我比較常用的兩條路分別是Cholesky分解法多元正態(tài)法先生成獨立標(biāo)準(zhǔn)正態(tài)樣本再用相關(guān)矩陣的Cholesky分解結(jié)果做線性變換讓樣本帶上目標(biāo)相關(guān)。然后對每個時段用標(biāo)準(zhǔn)正態(tài)CDF的逆函數(shù)轉(zhuǎn)到均勻空間再用目標(biāo)分布CDF的逆函數(shù)轉(zhuǎn)到真實出力空間。這套流程在數(shù)學(xué)上被稱為Nataf變換是MC相關(guān)性抽樣的經(jīng)典做法。Iman-Conover排序法先生成一個“模板樣本”讓其相關(guān)結(jié)構(gòu)接近目標(biāo)然后對原始獨立樣本的每一列按模板對應(yīng)列的秩順序重排使重排后的樣本具備模板的相關(guān)結(jié)構(gòu)。這個方法完全不要求原始樣本服從正態(tài)分布所以對風(fēng)電出力這種“非高斯、有邊界、大量零點”的數(shù)據(jù)特別友好。補充說明一下我這里提的兩種方法都是基于常見實踐的可靠方案。實際項目里不必迷信某一種關(guān)鍵是理解相關(guān)性在哪一步生成、在哪一步被壓縮。2.3 適合哪些數(shù)據(jù)和應(yīng)用場景這個方案最適用的場景是那些“時序依賴明顯、目標(biāo)分布非高斯、規(guī)模適中”的變量。典型代表就是風(fēng)電出力。風(fēng)速本身具有強自相關(guān)和日周期/季節(jié)周期風(fēng)機(jī)功率曲線又把風(fēng)速映射成帶有大量零點、上下截斷的出力序列分布形態(tài)高度非高斯直接用多元正態(tài)假設(shè)建模會誤差很大。這正是排序法配合經(jīng)驗分布最能發(fā)揮作用的地方。光伏出力也類似但有更強的確定性日內(nèi)模式。對光伏場景建模時我一般建議先除以理論晴空出力做歸一化再把歸一化后的“天空晴朗指數(shù)”序列作為建模對象這樣自相關(guān)結(jié)構(gòu)更容易提取場景也更穩(wěn)定。負(fù)荷序列、電價序列、徑流序列也都可以用這套框架。因為它們的共同特征是邊際分布偏斜、時序關(guān)聯(lián)強、同時有尖峰和持續(xù)狀態(tài)。需要提醒的是這套方法并不是萬能的。如果研究對象的變量之間不存在明顯時間依賴比如設(shè)備隨機(jī)的故障狀態(tài)、相互獨立的事件開關(guān)那就不需要費勁做時序相關(guān)性建模強行引入相關(guān)結(jié)構(gòu)反而會增加模型復(fù)雜度和誤差來源。3. 落地實操一時序相關(guān)性MC場景生成3.1 數(shù)據(jù)準(zhǔn)備與統(tǒng)計建模步驟在開始采樣之前數(shù)據(jù)準(zhǔn)備工作往往決定了整個項目能走多遠(yuǎn)。我習(xí)慣按下面這幾步做數(shù)據(jù)準(zhǔn)備第一步數(shù)據(jù)清洗。風(fēng)電、光伏出力數(shù)據(jù)里經(jīng)常出現(xiàn)異常尖峰、缺測、連續(xù)平臺值。這一步我會做三件事剔除超出物理范圍的記錄比如出力大于裝機(jī)容量、光伏夜間出力非零對缺測值做線性插值或前一天同時刻替換連續(xù)重復(fù)值過多的時段要格外警惕可能是數(shù)據(jù)采集故障需要用相鄰日數(shù)據(jù)修正。第二步確定建模對象。我的實踐經(jīng)驗是優(yōu)先對“預(yù)測誤差”建模而不是直接對出力絕對值建模。原因很簡單實際調(diào)度中預(yù)測曲線是已知信息需要模擬的是未來圍繞預(yù)測的波動。誤差建模時場景生成公式就是“確定性預(yù)測 MC誤差場景”不僅物理含義清晰也方便后期做場景削減和隨機(jī)優(yōu)化的銜接。第三步估計邊際分布。這里不要想當(dāng)然地套正態(tài)分布。風(fēng)電出力的邊際分布往往是“零出力點的概率集中 高出力區(qū)間厚尾”的組合形態(tài)光伏出力在歸一化后也明顯偏斜。比較穩(wěn)妥的辦法是直接采用經(jīng)驗分布把歷史數(shù)據(jù)按數(shù)值排序用排序位置除以樣本數(shù)作為累積概率估計再用分段線性插值或核密度估計平滑一下。這樣不需要做任何參數(shù)分布假設(shè)模型風(fēng)險低后續(xù)逆變換采樣也簡單。第四步估計相關(guān)結(jié)構(gòu)。對按“天”組織的數(shù)據(jù)把每天的數(shù)據(jù)作為一個樣本統(tǒng)一減去當(dāng)天均值并除以標(biāo)準(zhǔn)差得到標(biāo)準(zhǔn)化矩陣。然后用np.corrcoef計算時段間的相關(guān)矩陣。這里要注意相關(guān)矩陣的估計對異常值很敏感所以第一步的數(shù)據(jù)清洗不能偷懶。下面給出一個可以直接上手的Python代碼骨架import numpy as np import pandas as pd # data_shape: (days, T)T為每日時段數(shù)day為歷史天數(shù) daily_std (data - data.mean(axis0)) / data.std(axis0) # 時段間Pearson相關(guān)矩陣形狀 (T, T) R np.corrcoef(daily_std.T) # 自相關(guān)結(jié)構(gòu)檢查看lag 1的相關(guān)系數(shù) lag1_acf np.mean([np.corrcoef(daily_std[:, t], daily_std[:, t1])[0, 1] for t in range(daily_std.shape[1] - 1)]) print(lag1自相關(guān)均值: , lag1_acf)T通常取24小時級或9615分鐘級。相關(guān)矩陣規(guī)模分別是24×24和96×96都不算大計算成本可以忽略。3.2 Cholesky / 排序法生成相關(guān)樣本的細(xì)節(jié)對比數(shù)據(jù)準(zhǔn)備完成后就進(jìn)入核心環(huán)節(jié)生成帶相關(guān)性的采樣時序。我以誤差建模為例把兩條技術(shù)路線展開說清楚。路線一Cholesky分解法的完整流程如下對誤差歷史序列做經(jīng)驗分布擬合得到每個時段的邊際分布F_t。生成一個N×T的獨立標(biāo)準(zhǔn)正態(tài)矩陣ZN是需要生成的場景數(shù)。對相關(guān)矩陣R做Cholesky分解得到下三角矩陣L滿足R L L^T。計算Y Z L^T此時Y的相關(guān)矩陣接近R。對Y的每個元素做Phi(Y_t)得到均勻分布樣本再代入F_t^{-1}得到誤差場景。代碼段長這樣from scipy.stats import norm # 1) 估計誤差邊際分布這里用empirical CDF 線性插值 # 為簡化演示假設(shè)每時段誤差服從正態(tài)分布并存儲 mean_t, std_t # 實際項目中建議use經(jīng)驗分布代碼量略長但更穩(wěn) L np.linalg.cholesky(R 1e-6 * np.eye(R.shape[0])) # 加ridge防非正定 Z np.random.randn(N, T) Y Z L.T U norm.cdf(Y) # (N, T) err_scenes np.array([F_inv_t(u_t) for u_t in U.T]).T路線二Iman-Conover排序法邏輯不太一樣生成一個模板矩陣X_tmpl直接用多元正態(tài)分布協(xié)方差矩陣設(shè)置成目標(biāo)R。這一步可以用Cholesky實現(xiàn)但即使模板相關(guān)矩陣略有偏差后續(xù)排序也會校正。對模板每一列計算“秩信息”也就是按數(shù)值大小排名的位置索引。對誤差獨立樣本矩陣S的每一列按照模板對應(yīng)列的秩順序重新排列。排序完成后S的秩相關(guān)結(jié)構(gòu)就與模板一致了因此相關(guān)性也被傳遞到了目標(biāo)分布空間。這段邏輯最巧妙的地方是它只利用秩信息不需要計算任何逆變換對任何邊際分布都適用。我在實際項目中更偏愛Iman-Conover排序法原因后面會詳細(xì)展開。但在實現(xiàn)時有個細(xì)節(jié)模板矩陣必須固定種子否則每次排序后的樣本雖然相關(guān)結(jié)構(gòu)一致卻會引入額外的抽樣噪聲影響復(fù)現(xiàn)性。兩種路線的對比如下對比維度Cholesky分解法Iman-Conover排序法對邊際分布要求需要已知并計算逆CDF無要求經(jīng)驗分布即可相關(guān)矩陣要求必須正定允許適當(dāng)非正定非線性相關(guān)性壓縮問題存在需要校正較小排序保留了秩相關(guān)實現(xiàn)復(fù)雜度低幾行代碼中等需要維護(hù)排序索引適用數(shù)據(jù)近似正態(tài)、高斯耦合風(fēng)電、光伏等非高斯強偏數(shù)據(jù)3.3 從標(biāo)準(zhǔn)化空間還原到物理場景空間生成誤差場景后最后一步是還原到物理出力空間。還原本身不復(fù)雜但有三個細(xì)節(jié)容易被忽略。細(xì)節(jié)一物理邊界約束。風(fēng)電出力必須落在[0, 裝機(jī)容量]之間光伏出力在夜間必須為0。逆變換采樣后要對場景做截斷處理。我建議不要簡單硬截斷而是把超出上限的部分映射回上限同時把低于下限的部分映射回0。硬截斷會改變邊際分布的尾部形狀但在工程上通常影響可控所以先截斷再在評價階段復(fù)核分布即可。細(xì)節(jié)二場景形態(tài)平滑性。還原后的場景可能帶有高頻抖動。有人喜歡加一個滑動平均平滑我要提醒平滑濾波會改變時序相關(guān)結(jié)構(gòu)尤其是lag 1、lag 2的相關(guān)系數(shù)會被顯著拉低。如果后續(xù)要用ACF來評價場景質(zhì)量平滑必須要謹(jǐn)慎或者干脆不做。細(xì)節(jié)三生成場景的最終物理公式。如果采用誤差建模最終場景為[ P_{scene}(t) P_{forecast}(t) err_{scene}(t) ]隨后執(zhí)行截斷和必要修正。如果直接對出力建模則直接用逆變換采樣結(jié)果作為場景。另外我強烈建議在生成完場景后做一次“回檢”計算生成場景集合的ACF矩陣、邊際分布直方圖和目標(biāo)應(yīng)用指標(biāo)把這些與真實歷史數(shù)據(jù)放在一起對比。這一步雖然簡單卻能在進(jìn)入削減環(huán)節(jié)之前發(fā)現(xiàn)建模階段的系統(tǒng)性問題省去后續(xù)返工。4. 落地實操二場景削減與權(quán)重重構(gòu)4.1 場景削減的數(shù)學(xué)目標(biāo)和實用取舍MC生成一兩千個場景并不難難的是把這些場景放進(jìn)優(yōu)化模型。每個場景對應(yīng)一組約束和變量場景數(shù)量上升模型規(guī)模會成倍膨脹。削減的目的就是在盡量保留分布信息的前提下把場景集壓縮到計算可承受的規(guī)模。場景削減的數(shù)學(xué)目標(biāo)核心是控制“概率分布距離”。業(yè)界常用的距離是Kantorovich距離它在離散場景集上的形式非常直觀從原始場景集里刪掉某些場景并把它們的概率權(quán)重加到剩余場景上使削減后集合與原始集合之間“搬運概率質(zhì)量的距離”最小。這個目標(biāo)可以用一個生活例子來理解你有一個班的成績單大量場景想挑幾個代表參加座談削減后場景既要讓沒來的人覺得“代表表達(dá)了我的情況”又不能讓代表數(shù)量太多計算量。Kantorovich距離就是在量化“代表們能不能最好地覆蓋全班的各種情況”。在實際項目里我建議在削減之前先想清楚應(yīng)用場景因為這會影響距離度量里的權(quán)重。如果場景最終用于日前機(jī)組組合那么峰谷爬坡時段和晚高峰時段的重要性要更高如果用于儲能容量配置則連續(xù)多日低出力場景需要特殊保護(hù)不能被輕易“合并”掉。4.2 快速前向削減與同步回代削減的實現(xiàn)要點場景削減算法有兩大經(jīng)典路線快速前向選擇fast forward selection和同步回代消除simultaneous backward reduction?？焖偾跋蜻x擇的核心邏輯是“從空集開始反復(fù)挑選最有代表性的場景加入保留集”。每一步都找這樣一個場景把它加入保留集之后剩余場景到保留集的距離總和最小。重復(fù)M次得到M個保留場景。優(yōu)點是可以控制保留數(shù)量缺點是每一步都要遍歷所有剩余場景計算量偏大。同步回代消除的邏輯則相反從全集開始每次找到距離最近的兩個場景把其中一個刪除把它的概率權(quán)重加到另一個上面。不斷重復(fù)直到保留到目標(biāo)數(shù)量。它的計算量通常更小但有一個隱患如果兩個最近的場景都來自同一個區(qū)域合并會把極端但不常出現(xiàn)的場景提前吃掉。我在實際項目中常用“先粗分桶再精細(xì)削減”的兩階段策略。先用k-means或者k-medoids把2000個場景粗分成50個左右簇然后對每簇內(nèi)部用快速前向選擇選出最終的代表場景。這樣既控制了計算量又保證了代表場景的覆蓋度。注意這里的k-means聚類只是為了縮減搜索空間最終場景距離計算還是用原始時序距離度來完成不直接用聚類中心當(dāng)場景因為聚類中心會破壞物理時序關(guān)系。下面給一個簡化版的快速前向選擇實現(xiàn)骨架def forward_selection(scenes, weights, M): # scenes: (N, T) 原始場景集 # weights: (N,) 每個場景初始權(quán)重如1/N N, T scenes.shape selected [] remaining set(range(N)) for _ in range(M): best_j None best_cost np.inf for j in remaining: # 計算若選擇j作為保留場景剩余點到j(luò)的加權(quán)距離之和 dists np.linalg.norm(scenes[list(remaining)] - scenes[j], axis1) cost np.sum(weights[list(remaining)] * dists) if cost best_cost: best_cost cost best_j j selected.append(best_j) remaining.remove(best_j) # 權(quán)重重構(gòu)每個保留場景吸收剩余場景中距離它最近的那些權(quán)重 new_weights np.zeros(M) for i in range(N): if i in selected: idx selected.index(i) new_weights[idx] weights[i] else: dists [np.linalg.norm(scenes[i] - scenes[j]) for j in selected] idx int(np.argmin(dists)) new_weights[idx] weights[i] return scenes[selected], new_weights代碼里最核心的循環(huán)可以進(jìn)一步用向量化和堆結(jié)構(gòu)優(yōu)化但對幾千個場景、T在24~96的項目直接這樣寫也完全夠用。同步回代消除的實現(xiàn)思路更短維護(hù)場景集合和權(quán)重每次計算兩兩距離找到最小距離對將一個場景的權(quán)重合并給另一個然后刪除該場景。代碼我就不展開了和上面的循環(huán)結(jié)構(gòu)類似只是合并邏輯相反。4.3 削減后場景的時序相關(guān)性驗證削減完成不等于工作結(jié)束。我見過太多項目在削減環(huán)節(jié)把好不容易生成的時序相關(guān)性又給丟了原因是削減算法里用的距離度量沒有照顧到時序結(jié)構(gòu)。如果你直接用歐氏距離衡量兩個場景的差異那么兩個“整體數(shù)值接近但前后順序打亂”的場景會被認(rèn)為很相似。這會導(dǎo)致削減后的場景集把“時間反轉(zhuǎn)”、“相位偏移”的場景合并到同一個代表場景里自相關(guān)結(jié)構(gòu)自然受損。我在某風(fēng)電場項目中做過對比用純歐氏距離把2000個場景削減到20個削減后場景集合的lag 1自相關(guān)系數(shù)從0.83掉到0.51時序特征嚴(yán)重退化。改用加權(quán)歐氏距離在距離計算中額外加入相鄰時段差分的懲罰項后削減后場景集的lag 1自相關(guān)系數(shù)恢復(fù)到0.76已經(jīng)能比較好地保留時序結(jié)構(gòu)。加懲罰項的距離度量可以寫成[ d(\xi_i, \xi_j) \sum_{t1}^{T} w_t \cdot (\xi_{i,t} - \xi_{j,t})^2 \lambda \sum_{t1}^{T-1} ((\xi_{i,t1}-\xi_{i,t}) - (\xi_{j,t1}-\xi_{j,t}))^2 ]其中第一項衡量數(shù)值接近程度第二項衡量相鄰時段變化率近似爬坡的接近程度(\lambda)是權(quán)重系數(shù)。我在實踐中取(\lambda1.5 \sim 2)效果較好。如果你覺得這個加權(quán)不好調(diào)也可以直接用動態(tài)時間彎曲DTW距離它對相位偏移和局部時間變形更魯棒但計算量更大場景多時可以先粗聚類再精細(xì)化。5. 評價指標(biāo)體系與結(jié)果展示框架5.1 分布精度指標(biāo)評價場景集質(zhì)量的第一個維度是它到底還原了多少概率分布信息。這里至少要看三個指標(biāo)一是邊際分布的KS檢驗。對每個時段把生成場景在該時段的CDF與歷史數(shù)據(jù)的CDF做比較取最大差異的p值。如果大部分時段的p值大于0.05說明邊際分布擬合可以接受。不能只報一個“平均p值”要給出最差時段的p值因為最差時段往往對應(yīng)傍晚或清晨的快速爬坡時段恰恰是優(yōu)化模型最敏感的地方。二是1-Wasserstein距離。這個指標(biāo)可以直接用具體數(shù)值來衡量削減帶來的分布漂移。通常做法是畫出“保留場景數(shù)M vs 1-Wasserstein距離”的衰減曲線。當(dāng)M從2增加到5、10、15時距離快速下降隨后進(jìn)入平臺期。曲線的拐點就是“投入更多場景數(shù)換來更少信息損失”的臨界點直接幫助確定合理的削減數(shù)量。三是極端分位數(shù)對比。比如看場景集在1%分位數(shù)和99%分位數(shù)上是否與歷史數(shù)據(jù)一致。風(fēng)電場景的極端低出力分位數(shù)尤其重要因為它直接決定備用容量和儲能配置。僅看均值、方差完全不夠分位數(shù)對尾部信息的刻畫更細(xì)致。我個人的經(jīng)驗是這三個指標(biāo)都要做但不必追求每一項都達(dá)到極其精確重點是讓它們進(jìn)入一個“可接受區(qū)間”保證場景集在統(tǒng)計意義上站得住。5.2 時序結(jié)構(gòu)保留度指標(biāo)評價時序相關(guān)性的直接指標(biāo)是ACF曲線對比。把歷史數(shù)據(jù)的自相關(guān)系數(shù)序列l(wèi)ag 1lag 24和生成場景集合的平均ACF畫在同一張圖上觀察兩者的偏差。更嚴(yán)格一點可以計算ACF曲線的平均絕對誤差MAE。另一個特別能反映場景質(zhì)量的指標(biāo)是“相鄰時段差值分布”。先計算每條場景的差分序列(\Delta_t P_{t1} - P_t)然后比較生成場景集合與歷史數(shù)據(jù)的(\Delta_t)分布。這個指標(biāo)直接抓住“爬坡分布”的形態(tài)如果時序相關(guān)性丟失差分分布的尾部往往會明顯變寬。再一個指標(biāo)是秩相關(guān)矩陣差異。用Spearman秩相關(guān)而不是Pearson相關(guān)因為秩相關(guān)對異常值更穩(wěn)健。比較削減前后場景集合與歷史數(shù)據(jù)的秩相關(guān)矩陣用Frobenius范數(shù)刻畫整體偏差。這個范數(shù)就是一個把矩陣所有元素差異平方再開方的綜合數(shù)值越小說明相關(guān)結(jié)構(gòu)保留得越好。我在項目匯報時習(xí)慣把ACF圖和差分分布圖放在一起展示效果比單純報一個數(shù)值直觀得多也讓評審的人一眼明白“時序相關(guān)性”到底體現(xiàn)在哪里。5.3 面向應(yīng)用的誤差評估統(tǒng)計指標(biāo)再好最終都要回歸到應(yīng)用端驗證。我的做法是把場景集塞進(jìn)一個簡化版兩階段隨機(jī)規(guī)劃模型比如簡化的機(jī)組組合或儲能優(yōu)化比較兩種情形下的結(jié)果偏差情形一使用削減前的全部場景作為基準(zhǔn)理想值。情形二使用削減后的場景集。對比目標(biāo)函數(shù)值的相對偏差。如果削減后成本比全場景結(jié)果明顯偏低說明丟失了某些“惡劣場景”導(dǎo)致優(yōu)化結(jié)果過于樂觀。如果明顯偏高說明保留了太多冗余場景削減效率不足。除此之外還有一個容易被忽略的點約束滿足情況。削減后場景集中如果恰好沒有覆蓋某個關(guān)鍵極端場景那么一些安全約束可能在削減后的模型里永遠(yuǎn)不會被觸發(fā)得到的決策方案在真實世界中會過于冒進(jìn)。針對這個問題我建議在削減時明確保留“歷史極端案例”也就是把歷史數(shù)據(jù)中出現(xiàn)過的連續(xù)低出力、持續(xù)高出力、劇烈爬坡等特殊時段序列作為“種子場景”強制加入保留集再在此基礎(chǔ)上做削減。6. 常見問題與排查技巧實錄6.1 相關(guān)矩陣非正定的頑固問題MC場景生成環(huán)節(jié)最常見的問題是Cholesky分解報錯提示相關(guān)矩陣不是正定矩陣。這通常發(fā)生在T比較大比如96且樣本數(shù)有限時。原因有兩個一是長序列自相關(guān)矩陣本質(zhì)上是Toeplitz結(jié)構(gòu)其特征值在短樣本估計下很容易出現(xiàn)負(fù)值二是數(shù)據(jù)清洗不干凈時異常值會扭曲相關(guān)矩陣估計。排查和解決辦法按順序試給相關(guān)矩陣對角線加一個小的正則項例如(R \epsilon I)(\epsilon)取1e-61e-4。這相當(dāng)于給每個時段的方差增加一點“數(shù)值穩(wěn)定墊”如果加了正則還是不行計算最接近R的半正定矩陣用特征值分解把負(fù)特征值置零再歸一化直接放棄Cholesky分解改用Iman-Conover排序法。排序法不要求相關(guān)矩陣嚴(yán)格正定實踐中魯棒性最好。我在項目里用了排序法之后就再也沒被這個錯誤卡住過。但如果論文或報告中要求可復(fù)現(xiàn)的數(shù)值實驗則仍需要記錄清楚相關(guān)矩陣估計的具體方式和正則參數(shù)。6.2 削減后場景特征漂移怎么辦如果你發(fā)現(xiàn)削減后的場景集邊際分布看著還行但ACF明顯退化問題大概率出在削減階段的距離度量上。前文提到純歐氏距離不感知時序結(jié)構(gòu)會把時間形態(tài)不同的場景合并起來。解決辦法有兩個方向在距離里加入差分懲罰項邏輯見4.3節(jié)。這是成本最低的改動。改用“代表性場景按聚類中心搜索”的思路先聚類再在每個類內(nèi)搜索離類中心最近的原始場景作為代表避免直接拿類中心當(dāng)場景。另外還要檢查一下削減后場景權(quán)重是否出現(xiàn)極端值。如果某個保留場景的權(quán)重超過0.25說明它“吸收”了過多場景很可能抹平了多種不同的時序模式。這種情況下建議強制提高保留場景數(shù)量或者在削減前做“極端場景保護(hù)”。6.3 場景數(shù)量選擇糾結(jié)怎么破很多初學(xué)者會問到底削減到多少個場景合適。這個問題很難給統(tǒng)一答案但我建議用“Wasserstein衰減曲線 項目計算預(yù)算”雙標(biāo)準(zhǔn)來定。具體操作畫出M從2到30的1-Wasserstein距離曲線找曲線拐點。如果拐點在M15附近那就設(shè)在15如果計算資源緊張比如要做上千次隨機(jī)優(yōu)化再往下壓到510同時確認(rèn)應(yīng)用端誤差在可接受范圍內(nèi)。還有一個經(jīng)驗是不要把M設(shè)得太小。當(dāng)M低于10時場景集對極端事件的分辨率會快速下降備用配置結(jié)果容易出現(xiàn)系統(tǒng)性偏差。我做過一次對比M從8增加到12兩階段隨機(jī)規(guī)劃的目標(biāo)函數(shù)值變化量達(dá)到3%以內(nèi)而M繼續(xù)增加到20時變化已經(jīng)小于0.5%這說明M12左右就已經(jīng)合理。6.4 生成場景“看起來不夠隨機(jī)”怎么理解有時生成出來的場景集會在某個時段內(nèi)“擠成一團(tuán)”看起來多樣性不足。我遇到過幾次原因不是采樣失敗而是相關(guān)矩陣估計過于強烈。比如某風(fēng)電場出力數(shù)據(jù)有很強的日內(nèi)同漲同跌模式自相關(guān)矩陣的所有元素都接近0.9那么采樣出的場景當(dāng)然會高度相似。如果你希望場景集保留更多多樣性可以在相關(guān)矩陣上做一個“收縮”處理比如把相關(guān)矩陣乘以一個衰減系數(shù)(\rho0.85)讓相關(guān)結(jié)構(gòu)溫和化。但注意這種做法會犧牲一點精度所以我更建議在評價時先確認(rèn)“多樣性不足”是否真的影響應(yīng)用結(jié)果而不是憑感覺調(diào)參數(shù)。7. 寫在后面的一些個人經(jīng)驗這個項目最讓我受益的一點是它逼著我把“不確定性模擬”從抽象概念落到了具體工程對象上。以前我做MC抽樣關(guān)注的只是“樣本夠不夠多、分布對不對”很少去管這些樣本放到?jīng)Q策模型里會產(chǎn)生什么后果。直到被爬坡約束的真實案例打臉我才真正意識到時序相關(guān)性不是統(tǒng)計潔癖而是決策模型本身對場景軌跡提出的硬需求。如果現(xiàn)在有人讓我從零做這類研究我會建議走一條最穩(wěn)的閉環(huán)路線先拿歷史數(shù)據(jù)估計經(jīng)驗分布和相關(guān)矩陣用Iman-Conover排序法生成2000個場景再用帶差分懲罰的加權(quán)歐氏距離做快速前向削減最后用ACF對比和兩階段優(yōu)化偏差來驗證效果。這條鏈路每一步都可復(fù)現(xiàn)、可解釋也容易排查問題。等這條鏈路跑通之后再去嘗試Copula函數(shù)、變分自編碼器、生成對抗網(wǎng)絡(luò)這些更花哨的生成手段就有了一個扎實的對比基線。最后再分享一個實踐中額外收獲的小技巧在做場景削減時永遠(yuǎn)保留至少一條“歷史最嚴(yán)苛場景”。不管削減算法怎么選把這個場景強制加入保留集能在后續(xù)的隨機(jī)優(yōu)化和安全校核中省掉很多解釋成本和返工時間。這個方法簡單到很多人不屑于用但恰恰是我在多個項目里避免“樂觀偏差”最有效的土辦法。