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