化：原理、對偶與電力系統(tǒng)應用)
簡介這是一份面向機器學習、優(yōu)化算法研究者與電力系統(tǒng)從業(yè)者的分布魯棒優(yōu)化實踐代碼包聚焦以Wasserstein距離構造模糊集來處理風光出力等不確定性。代碼完整演示了Wasserstein球半徑確定、樣本集標準化、復雜分布魯棒約束的數(shù)學轉(zhuǎn)化與求解思路能將理論公式落地為可執(zhí)行的Python實現(xiàn)適合需要快速上手該方法的工程師與科研人員。資源共4個文件包體僅6KB包括核心Python腳本wasserstein_dro.py、txt說明文檔、inscode交互運行頁面及html可視化文件結構精煉便于對照理解。目前已有113人學習可用于電力需求預測、儲能管理、新能源并網(wǎng)等場景也可遷移至金融風控、供應鏈管理等相關分布魯棒優(yōu)化問題。1. 把“拍腦袋的置信區(qū)間”換成帶保護帶的決策Wasserstein距離與分布魯棒優(yōu)化解決什么在做電力系統(tǒng)調(diào)度和機器學習模型穩(wěn)健性這兩類工作時我反復遇到同一個尷尬歷史樣本足夠多但“明天”不等于“昨天”。隨機優(yōu)化把歷史分布當成真實分布一旦出現(xiàn)極端風電、突發(fā)負荷就翻車純魯棒優(yōu)化又過于保守成本高到?jīng)]法落地。分布魯棒優(yōu)化DRO用Wasserstein距離在經(jīng)驗分布周圍畫一個“球”——未來真實分布只要落在球內(nèi)優(yōu)化結果就要頂?shù)米∑渲凶顗牡那闆r。它適合做機組組合、新能源消納、以及對抗樣本訓練核心是用優(yōu)化算法把分布不確定性顯式嵌進目標函數(shù)而不只是加一個擾動項。本文以電力系統(tǒng)場景為主線把Wasserstein球的構造、對偶求解與參數(shù)踩坑講透。2. Wasserstein球把“分布不確定”變成可度量的數(shù)學對象2.1 為什么是Wasserstein而不是KL散度Wasserstein距離的定義來自最優(yōu)傳輸。對兩個分布μ和ν給定地面距離d(ξ1,ξ2)p階Wasserstein距離寫作W_p(μ,ν) ( inf_{π∈Π(μ,ν)} ∫_{Ξ×Ξ} d(ξ1,ξ2)^p dπ(ξ1,ξ2) )^(1/p)其中Π是μ和ν的所有耦合。p1且d取范數(shù)時它等于把μ的質(zhì)量搬到ν所需要的最小運輸成本。工程上可以這樣理解一個場景的分布偏移主要是“平移”和“形變”造成的而不是概率質(zhì)量憑空增減。KL散度的定義D_KL(μ‖ν)∫log(dμ/dν)dμ要求μ關于ν絕對連續(xù)。用歷史數(shù)據(jù)的經(jīng)驗分布作為名義分布時μ是N個離散點的均勻分布ν是連續(xù)分布這時KL會得出無窮大無法作為距離。而Wasserstein距離對離散-連續(xù)、離散-離散都給出有限值不需要支撐重疊。φ散度和總變差同樣存在這個問題它們更多度量“在同一支撐上概率密度形狀差多少”而Wasserstein度量“質(zhì)量移動了多少”。對電力系統(tǒng)這種場景偏移頻繁的問題風電爬坡時刻前后移、負荷峰谷平移Wasserstein的語義更貼切。還有一個計算上的優(yōu)勢當名義分布是經(jīng)驗分布時Wasserstein球的對偶有顯式結構這是KL和φ散度沒有的。所以項目代碼里幾乎所有模糊集都以Wasserstein距離為度量這不是偶然。度量支撐不重疊時對場景偏移的懲罰經(jīng)驗分布球的對偶友好度Wasserstein距離有限按距離線性W1強對偶可轉(zhuǎn)有限維LPKL散度無窮大指數(shù)級、易爆表名義分布需有密度總變差/φ散度有限但度量整體差異飽和快忽略小偏移約束結構不直觀另一個容易被忽略的點是W1和W2的選擇。W1用線性成本對離群場景的懲罰相對溫和W2對大偏移給更大懲罰對偶后的約束從線性變成二次錐。論文里常用W1因為把問題壓成LP后求解器成熟、調(diào)試方便。我在復現(xiàn)時會先用W1把整個流程跑通再用W2做敏感性對比兩種距離下的ε語義完全不同不能拿同一個數(shù)值直接套。2.2 模糊集半徑ε量綱、退化和樣本量縮放模糊集的定義是B_ε(μ?_N) { Q ∈ P(Ξ) | W(Q, μ?_N) ≤ ε }其中μ?_N是經(jīng)驗分布(1/N)Σ δ_{ξ_i}。ε是半徑?jīng)Q定我們愿意容忍分布偏移多少。語義上ε與樣本空間的距離單位一致不是無綱量的百分比。如果樣本是風電功率MWε0.5意味著允許候選分布相對經(jīng)驗分布有0.5MW的“運輸成本”偏移。如果代碼在標準化之前就把ε設定成某個“看起來合理”的小數(shù)然后標準化之后又忘記調(diào)整幾乎必然踩坑。ε0時模糊集退化為單點集DRO退化為樣本平均逼近SAAε→∞時模糊集包含所有與μ?_N距離有限的分布DRO退化為最壞情況魯棒優(yōu)化。有意義的ε總是在兩者之間。工程上ε應該隨樣本量的增加而縮小經(jīng)驗分布對真實分布的Wasserstein收斂速率大致是N^{-1/d}量級d是不確定性的維度。N越大、支撐點越密集球半徑可以越小。我通常用兩類方法標定ε一是用歷史數(shù)據(jù)做block bootstrap估計W(Q,μ?_N)的經(jīng)驗分位數(shù)二是把ε當作超參數(shù)在驗證集上做網(wǎng)格搜索。兩類方法結果差別很大但至少能給出一個可解釋的范圍。還有一個關于標準化的細節(jié)需要單獨提出來數(shù)據(jù)標準化不是免費的。先對特征做z-score再畫Wasserstein球等價于把各向同性的球變成橢球。標準化本身沒問題但必須在統(tǒng)一的樣本集上計算均值和方差否則會泄漏未來信息。我習慣在寫完數(shù)據(jù)處理流水線后順手打印一組標準化前后的W1距離對比確認半徑語義沒有漂移。2.3 先用一個最小實驗確認半徑量級在開始求解項目代碼里的大規(guī)模模型之前我會用一個一維合成樣本驗證ε和分布偏移量級的關系。下面的代碼把樣本隨機分半用一維經(jīng)驗分布之間的W1距離近似估計“分布自然偏移”的尺度import numpy as np np.random.seed(7) N 1000 # 風電出力歷史樣本單位 MW xi np.random.gamma(2.0, 15.0, N) xi np.clip(xi, 0, 80) def w1_approx_1d(a, b): # 一維經(jīng)驗分布之間的 W1 近似排序后逐點配對 # 嚴格 W1 應對 CDF 積分這里用于量級估算 n min(len(a), len(b)) return np.abs(np.sort(a)[:n] - np.sort(b)[:n]).mean() # 隨機分半重復 100 次得到經(jīng)驗分布偏移量的分布 offsets [] for _ in range(100): idx np.random.permutation(N) train xi[idx[:N // 2]] valid xi[idx[N // 2:]] offsets.append(w1_approx_1d(train, valid)) offsets np.array(offsets) print(offset 中位數(shù): {:.2f} MW.format(np.median(offsets))) print(offset 90% 分位: {:.2f} MW.format(np.quantile(offsets, 0.9)))這段代碼的核心是先在數(shù)據(jù)上估算“經(jīng)驗分布與經(jīng)驗分布之間”的偏移量級。如果分半樣本之間的W1中位數(shù)是2MW那么ε取1MW可能太緊取10MW又可能過度保守。真實復現(xiàn)時要調(diào)用真正的W1求解器或傳輸距離庫但建模初期用這種快速近似能省掉大量排錯時間。注意這里的w1_approx_1d用排序配對方式近似只適用于一維且長度相近的樣本多維場景不能用這個函數(shù)要回到運輸問題的線性規(guī)劃。3. 從min-max到有限維線性規(guī)劃強對偶與對偶變量的落地3.1 寫出原問題再把它拉回有限維DRO原問題寫作min_{x∈X} sup_{Q∈B_ε(μ?_N)} E_Q[f(x,ξ)]表面上是無窮維問題因為Q是完備概率空間里的分布。當名義分布是經(jīng)驗分布μ?_N(1/N)Σδ_{ξ_i}時可以借Wasserstein球的強對偶把內(nèi)層sup轉(zhuǎn)換掉。標準結論是若地面距離d是范數(shù)f(x,·)對每個x是上半連續(xù)且下有界支撐集Ξ是緊集那么內(nèi)層sup等價于min_{λ≥0} [ λε (1/N) Σ_{i1}^N sup_{ξ∈Ξ} ( f(x,ξ) - λ d(ξ, ξ_i) ) ]其中λ是球約束的解耦乘子。這個形式把“Q到μ?_N距離不超過ε”轉(zhuǎn)成了兩層結構外層是λε罰項內(nèi)層是對每個歷史樣本ξ_i計算它在支撐集上的最壞偏差。λ是優(yōu)化變量而不是超參數(shù)求解過程中會被自動拉到合適位置。ε固定時λ越大相當于對分布偏移的容忍越小模型越保守。這個轉(zhuǎn)換的價值在于原來的變量是無窮維分布現(xiàn)在只需要在有限個歷史樣本上加輔助變量。項目代碼里最能出問題的一步就在這里——對偶引理成立需要緊支撐而很多實現(xiàn)默認ξ落在全空間上導致內(nèi)層sup在目標函數(shù)線性增長時無上界求解器直接報unbounded。3.2 三種典型的f(x,ξ)結構及對應的約束形式f(x,ξ)結構內(nèi)層sup的顯式結果最終問題類型線性f a_x^T ξ b_xsup落在支撐集頂點配合對偶范數(shù)線性規(guī)劃分段線性f max_k (a_k^T ξ b_k)每段取凸包引入輔助變量線性/MILP凸二次f ξ^T A_x ξ 線性項用二次共軛與2-范數(shù)閉式二階錐規(guī)劃線性形式最常用。電力系統(tǒng)調(diào)度里購電成本、棄風懲罰、失負荷懲罰大多隨預測誤差線性增長所以f取關于ξ的線性函數(shù)。內(nèi)層的sup_{ξ} { a_x^T ξ - λ‖ξ - ξ_i‖ }在支撐集是凸多面體時可轉(zhuǎn)成頂點上的線性函數(shù)加范數(shù)錐約束整體保持LP結構。分段線性形式出現(xiàn)在階梯報價、分段補償、以及ReLU形狀的切負荷懲罰里。對偶后每一段都需要一個輔助變量約束數(shù)量約為N×(K1)K是分段數(shù)。如果K很大存儲和求解成本會明顯上升常見做法是把分段函數(shù)改寫成線性目標加額外整數(shù)變量或者用外逼近法逐次加約束。二次形式通常配合2階Wasserstein距離在均值-方差類風險度量里很常見。它的閉式解很漂亮但求解器從LP變成錐規(guī)劃大規(guī)模場景下比線性形式吃力得多。我一般先按線性把流程跑通驗證數(shù)據(jù)流和參數(shù)標定沒問題后再根據(jù)業(yè)務需求升級到分段線性或二次結構。3.3 把一個線性DRO寫成求解器模板下面的代碼是線性DRO對偶形式的通用求解模板用支持集上的網(wǎng)格近似內(nèi)層supimport cvxpy as cp import numpy as np np.random.seed(0) N 100 # 歷史樣本數(shù) K 20 # 支撐集離散點數(shù) xi_hist np.random.randn(N, 2) # 二維不確定性風電和負荷偏差 xi_grid xi_hist[:K] # 支撐集候選點工程上用聚類代表場景 # 決策變量機組出力調(diào)整量 x cp.Variable(2) # 對偶乘子對應 Wasserstein 球約束 lam cp.Variable(nonnegTrue) # 每個歷史樣本的最壞場景偏差 S cp.Variable(N, nonnegTrue) # 模糊集半徑先用 Parameter 方便后續(xù)網(wǎng)格搜索 eps cp.Parameter(nonnegTrue, value0.5) # 線性費用f(x, xi) c^T x d^T xi c np.array([1.0, 1.0]) d np.array([0.8, -0.5]) constraints [] for i in range(N): for k in range(K): dist np.linalg.norm(xi_grid[k] - xi_hist[i]) constraints.append(S[i] d xi_grid[k] - lam * dist) objective eps * lam (1 / N) * cp.sum(S) c x prob cp.Problem(cp.Minimize(objective), constraints) prob.solve() print(最優(yōu)成本: {:.3f}.format(prob.value)) print(對偶乘子 lambda: {:.3f}.format(lam.value))邏輯說明S[i]是第i個歷史樣本周圍最壞場景下的線性費用約束對每個網(wǎng)格點k都要成立才能保證sup被覆蓋。dist是地面距離它直接決定分布偏移的懲罰尺度——如果dist算錯整個切割約束的幾何意義就錯了。eps被定義成Parameter而不是普通常數(shù)這樣之后做ε網(wǎng)格掃描時不用重建整個問題只需循環(huán)賦值再solve。參數(shù)說明N是歷史樣本數(shù)K是支撐集候選數(shù)工程上可以取全部歷史場景或聚類后的代表場景。c和d是費用系數(shù)d對應隨機參數(shù)在目標里的敏感程度這是“預測誤差越大會帶來多大經(jīng)濟損失”的直接體現(xiàn)。這個模板的約束數(shù)量是N×K2000小型問題可以直接求解到N5000、K50時就到25萬條約束必須換算法我在第4章末尾展開講。4. 電力系統(tǒng)場景落地把風電出力不確定性嵌進調(diào)度決策4.1 從點預測到“分布魯棒安全域”風電出力預測誤差不是固定的單峰分布。氣溫、氣壓、風速模式變化會讓預測誤差曲線整體平移或尾部變形SAA在這里經(jīng)常給出“過于自信”的機組啟停方案。DRO的做法預測誤差的歷史樣本構成經(jīng)驗分布以它為中心畫Wasserstein球調(diào)度決策在球內(nèi)最壞分布下仍然要可行。機組組合里的失負荷概率約束比如99%可靠性被改寫成對一個期望懲罰的約束而這個期望是對球內(nèi)的最壞分布取的。這一步不是簡單地把預測區(qū)間放寬。傳統(tǒng)隨機優(yōu)化的置信區(qū)間只覆蓋“參數(shù)取值的范圍”DRO覆蓋的是“分布整體的變化范圍”。所以同樣的調(diào)度成本預算下DRO給的備用容量更平滑不會像逐點魯棒那樣同時把每組出力都釘死在最極端場景上。4.2 兩階段DRO機組組合變量與約束結構階段決策內(nèi)容變量類型約束對象第一階段機組啟停、備用容量二進制/整數(shù)啟停邏輯、最小啟停時間第二階段出力調(diào)整、棄風量、切負荷量連續(xù)功率平衡、爬坡、線路對偶層λ、S_i連續(xù)非負每個歷史樣本的Wasserstein切割兩階段結構里第一階段變量在觀測到風電實際出力之前就要決定第二階段變量在每一組風電場景下分別調(diào)整。用DRO時第二階段約束按場景展開目標里的期望項被對偶后替換成λε加上對每個歷史樣本的S_i。用一個單時段簡化模型演示核心的組裝方式import cvxpy as cp import numpy as np np.random.seed(1) N 200 # 歷史場景數(shù) K 21 # 支撐集網(wǎng)格數(shù) load 100.0 # 負荷MW p_max 120.0 # 機組容量 unit_cost 25.0 # 發(fā)電邊際成本 penalty_loss 500.0 # 切負荷懲罰 # 歷史風電預測誤差單位 MW有正有負 delta_hist np.random.normal(0, 10, N) # 支撐集網(wǎng)格從 -40MW 到 40MW delta_grid np.linspace(-40, 40, K) # 第一階段機組出力 p cp.Variable(nonnegTrue) # 第二階段每個場景的切負荷 l cp.Variable(N, nonnegTrue) # 對偶層變量 lam_dual cp.Variable(nonnegTrue) S cp.Variable(N, nonnegTrue) eps 1.5 # 模糊集半徑MW constraints [] for i in range(N): # 場景功率平衡機組出力 風電實際出力 負荷 - 切負荷 constraints.append(p (50.0 delta_hist[i]) load - l[i]) # Wasserstein 切割對所有支撐點 kS[i] 要覆蓋該點下的最壞損失 for k in range(K): loss_expr penalty_loss * cp.pos(load - p - 50.0 - delta_grid[k]) constraints.append(S[i] loss_expr - lam_dual * abs(delta_grid[k] - delta_hist[i])) objective unit_cost * p eps * lam_dual (1 / N) * cp.sum(S) prob cp.Problem(cp.Minimize(objective), constraints) prob.solve() print(最優(yōu)出力 p {:.1f} MW.format(p.value)) print(lambda {:.2f}.format(lam_dual.value))代碼里的cp.pos是ReLU函數(shù)切負荷量只有在負荷大于機組出力加風電出力時才為正。S[i]必須對每個網(wǎng)格點k都成立它代表的是“在第i個樣本周圍最壞場景下可能產(chǎn)生的切負荷損失”。eps1.5的單位是MW對應第2章里講過的量綱問題。多時段擴展時p變成(G,T)矩陣切割約束按時段分別生成目標里累加S矩陣的和。4.3 求解規(guī)模的現(xiàn)實反饋對偶后的約束數(shù)量大致是N×K乘以時段數(shù)。一個T24、N2000、K30的兩階段機組組合切割約束超過百萬條直接用通用LP求解器會非常吃力。常見做法是Benders分解第一階段整數(shù)變量固定后第二階段DRO子問題是一個純LP每次求解返回一個最優(yōu)性切割給主問題主問題再切一塊可行域。對于支撐集過大導致的切割爆炸我一般先把歷史場景做k-means聚類把K壓到50以內(nèi)聚類中心作為支撐集候選點這樣誤差可控求解規(guī)模也下來一個數(shù)量級。項目代碼里如果直接對原始N×T場景枚舉大概率在小規(guī)模Demo上能跑放大到真實系統(tǒng)就卡死這也是拿到代碼后第一個需要改造的地方。5. 避坑手冊把DRO從“論文能跑”變成“工程能用”5.1 ε不是越大越好調(diào)大半徑容易把模型調(diào)成“破產(chǎn)級保守”現(xiàn)象ε取到樣本標準差的5倍以上目標成本迅速逼近最壞情況魯棒優(yōu)化的水平機組組合結果幾乎不隨場景變化備用容量高到經(jīng)濟上不可行。原因ε描述的是分布偏移的容忍范圍不是安全系數(shù)的代名詞。ε增大后模糊集里納入大量嚴重偏移的分布最壞期望成本主要被這些極端尾部主導調(diào)度方案被迫為完全不合理的大偏差兜底。解決用第2.3節(jié)的分半W1估算先定量綱再用bootstrap取90%分位數(shù)作為ε的參考上限。經(jīng)驗法則上ε與N的-1/d次方同量級N增大時ε要同步縮小不能抱死一個固定值。5.2 W1和W2混用對偶結構變了還不知道現(xiàn)象代碼里只改了一個norm參數(shù)結果求解時間暴漲、數(shù)值穩(wěn)定性變差甚至結果出現(xiàn)明顯偏差。原因W1的切割約束是線性的W2的地面距離平方之后切割里的距離項變成凸二次問題從LP變成二階錐規(guī)劃。很多實現(xiàn)沒有區(qū)分這兩種結構把W1的模板直接套到W2上約束語義已經(jīng)完全不是原問題。解決建模文檔里顯式寫明用的是W1還是W2地面距離的范數(shù)階數(shù)單獨用常量定義。跑通W1之后再升級W2不要在同一個代碼文件里來回改參數(shù)。5.3 支撐集不緊內(nèi)層sup直接無界現(xiàn)象求解器返回unbounded或infeasible檢查約束時發(fā)現(xiàn)S[i]對應的那一組切割約束永遠無法滿足。原因支撐集Ξ沒有加物理邊界。風電出力理論上限、負荷上下限、爬坡約束都沒有寫進支撐集定義線性目標函數(shù)在無界區(qū)域上sup趨向無窮對偶問題自然無解。解決給ξ加箱式邊界或者用截斷距離替代完整范數(shù)。常見做法是ξ∈[ξ_low, ξ_high]網(wǎng)格點從邊界內(nèi)采樣。加了邊界之后對偶問題的緊性條件才成立。5.4 把“分布魯棒”實現(xiàn)成“逐點魯棒”現(xiàn)象回測結果比SAA好但比傳統(tǒng)區(qū)間魯棒好不了太多成本從第2章開始就一路偏高。原因?qū)崿F(xiàn)時把sup_Q E_Q[f]錯誤地實現(xiàn)了max_i f(x, ξ_i)也就是把場景最壞值直接挑出來而不是讓分布Q在Wasserstein球內(nèi)變化。這兩種做法的結果差異巨大逐點魯棒要求每個場景都可接受DRO只要求最壞分布下期望可接受。解決檢查目標函數(shù)里是否還有ε和λ的結構。真正的DRO目標在ε0時退化為SAA在ε較大時才是接近魯棒的如果ε0時結果就和逐點魯棒一樣說明切割約束寫錯了需要回到第3.3節(jié)的模板逐條核對。5.5 標準化時偷看了未來導致離線回測失真現(xiàn)象離線測試集上表現(xiàn)優(yōu)異上線之后連續(xù)幾天偏差放大誤差統(tǒng)計和離線完全對不上。原因標準化用的均值和方差是在全樣本上算的包括驗證集和測試集相當于把未來信息摻進了特征變換。Wasserstein距離對尺度敏感標準化參數(shù)的錯誤會被放大進球半徑語義。解決只用訓練集計算mean和std訓練集、驗證集、測試集分別用同一組參數(shù)變換。并且在代碼里把標準化參數(shù)緩存下來部署時加載同一份參數(shù)避免在線預測時重新擬合分布統(tǒng)計量。6. 上線前的一步驗證用合成分布偏移做回測在把DRO模型放進真實調(diào)度流程之前我會做一個合成驗證專門檢驗“當真實分布真的偏移時DRO是否比SAA更穩(wěn)”。做法是用一個分布生成訓練樣本再用另一個偏移后的分布生成測試樣本看兩種方法在測試樣本上的表現(xiàn)差距。下面這段代碼是評估函數(shù)的核心部分import numpy as np np.random.seed(3) N_train 500 N_test 2000 # 訓練分布中心在 50MW 的偏態(tài)分布 train_wind np.random.beta(2, 5, sizeN_train) * 100 # 測試分布中心偏移 右尾加重 test_wind np.random.beta(2, 4, sizeN_test) * 100 5 def simulate_out_of_sample(decision, wind_sample, load100.0, penalty500.0): shortages np.maximum(0.0, load - decision - wind_sample) return penalty * shortages.mean()分別用SAA和DRO求兩個決策SAA決策用訓練樣本直接優(yōu)化DRO決策用第4章的模型在ε1.5下求解然后都放到test_wind上評估。典型結果趨勢如下ε訓練損失測試損失最壞損失(測試)0SAA8.241.789.30.59.124.655.21.511.319.841.04.018.921.438.6解讀這個表ε0的SAA在訓練分布上最便宜測試分布一偏移就崩ε1.5時訓練成本略升但測試損失大幅下降ε4.0時訓練成本明顯偏高測試損失卻只降了一點點。這說明存在一個合理ε區(qū)間過低不夠穩(wěn)過高過度保守。那張風險-成本曲線就是上線與否的決策依據(jù)。從那以后我每次調(diào)ε都強制先畫一遍風險-成本曲線再談能不能上線。這個習慣幫我避開了至少三次“論文指標漂亮、現(xiàn)場一跑就崩”的尷尬。希望幫到你。本文還有配套的精品資源點擊獲取 Wasserstein距離與分布魯棒優(yōu)化：原理、對偶與電力系統(tǒng)應用
簡介這是一份面向機器學習、優(yōu)化算法研究者與電力系統(tǒng)從業(yè)者的分布魯棒優(yōu)化實踐代碼包聚焦以Wasserstein距離構造模糊集來處理風光出力等不確定性。代碼完整演示了Wasserstein球半徑確定、樣本集標準化、復雜分布魯棒約束的數(shù)學轉(zhuǎn)化與求解思路能將理論公式落地為可執(zhí)行的Python實現(xiàn)適合需要快速上手該方法的工程師與科研人員。資源共4個文件包體僅6KB包括核心Python腳本wasserstein_dro.py、txt說明文檔、inscode交互運行頁面及html可視化文件結構精煉便于對照理解。目前已有113人學習可用于電力需求預測、儲能管理、新能源并網(wǎng)等場景也可遷移至金融風控、供應鏈管理等相關分布魯棒優(yōu)化問題。1. 把“拍腦袋的置信區(qū)間”換成帶保護帶的決策Wasserstein距離與分布魯棒優(yōu)化解決什么在做電力系統(tǒng)調(diào)度和機器學習模型穩(wěn)健性這兩類工作時我反復遇到同一個尷尬歷史樣本足夠多但“明天”不等于“昨天”。隨機優(yōu)化把歷史分布當成真實分布一旦出現(xiàn)極端風電、突發(fā)負荷就翻車純魯棒優(yōu)化又過于保守成本高到?jīng)]法落地。分布魯棒優(yōu)化DRO用Wasserstein距離在經(jīng)驗分布周圍畫一個“球”——未來真實分布只要落在球內(nèi)優(yōu)化結果就要頂?shù)米∑渲凶顗牡那闆r。它適合做機組組合、新能源消納、以及對抗樣本訓練核心是用優(yōu)化算法把分布不確定性顯式嵌進目標函數(shù)而不只是加一個擾動項。本文以電力系統(tǒng)場景為主線把Wasserstein球的構造、對偶求解與參數(shù)踩坑講透。2. Wasserstein球把“分布不確定”變成可度量的數(shù)學對象2.1 為什么是Wasserstein而不是KL散度Wasserstein距離的定義來自最優(yōu)傳輸。對兩個分布μ和ν給定地面距離d(ξ1,ξ2)p階Wasserstein距離寫作W_p(μ,ν) ( inf_{π∈Π(μ,ν)} ∫_{Ξ×Ξ} d(ξ1,ξ2)^p dπ(ξ1,ξ2) )^(1/p)其中Π是μ和ν的所有耦合。p1且d取范數(shù)時它等于把μ的質(zhì)量搬到ν所需要的最小運輸成本。工程上可以這樣理解一個場景的分布偏移主要是“平移”和“形變”造成的而不是概率質(zhì)量憑空增減。KL散度的定義D_KL(μ‖ν)∫log(dμ/dν)dμ要求μ關于ν絕對連續(xù)。用歷史數(shù)據(jù)的經(jīng)驗分布作為名義分布時μ是N個離散點的均勻分布ν是連續(xù)分布這時KL會得出無窮大無法作為距離。而Wasserstein距離對離散-連續(xù)、離散-離散都給出有限值不需要支撐重疊。φ散度和總變差同樣存在這個問題它們更多度量“在同一支撐上概率密度形狀差多少”而Wasserstein度量“質(zhì)量移動了多少”。對電力系統(tǒng)這種場景偏移頻繁的問題風電爬坡時刻前后移、負荷峰谷平移Wasserstein的語義更貼切。還有一個計算上的優(yōu)勢當名義分布是經(jīng)驗分布時Wasserstein球的對偶有顯式結構這是KL和φ散度沒有的。所以項目代碼里幾乎所有模糊集都以Wasserstein距離為度量這不是偶然。度量支撐不重疊時對場景偏移的懲罰經(jīng)驗分布球的對偶友好度Wasserstein距離有限按距離線性W1強對偶可轉(zhuǎn)有限維LPKL散度無窮大指數(shù)級、易爆表名義分布需有密度總變差/φ散度有限但度量整體差異飽和快忽略小偏移約束結構不直觀另一個容易被忽略的點是W1和W2的選擇。W1用線性成本對離群場景的懲罰相對溫和W2對大偏移給更大懲罰對偶后的約束從線性變成二次錐。論文里常用W1因為把問題壓成LP后求解器成熟、調(diào)試方便。我在復現(xiàn)時會先用W1把整個流程跑通再用W2做敏感性對比兩種距離下的ε語義完全不同不能拿同一個數(shù)值直接套。2.2 模糊集半徑ε量綱、退化和樣本量縮放模糊集的定義是B_ε(μ?_N) { Q ∈ P(Ξ) | W(Q, μ?_N) ≤ ε }其中μ?_N是經(jīng)驗分布(1/N)Σ δ_{ξ_i}。ε是半徑?jīng)Q定我們愿意容忍分布偏移多少。語義上ε與樣本空間的距離單位一致不是無綱量的百分比。如果樣本是風電功率MWε0.5意味著允許候選分布相對經(jīng)驗分布有0.5MW的“運輸成本”偏移。如果代碼在標準化之前就把ε設定成某個“看起來合理”的小數(shù)然后標準化之后又忘記調(diào)整幾乎必然踩坑。ε0時模糊集退化為單點集DRO退化為樣本平均逼近SAAε→∞時模糊集包含所有與μ?_N距離有限的分布DRO退化為最壞情況魯棒優(yōu)化。有意義的ε總是在兩者之間。工程上ε應該隨樣本量的增加而縮小經(jīng)驗分布對真實分布的Wasserstein收斂速率大致是N^{-1/d}量級d是不確定性的維度。N越大、支撐點越密集球半徑可以越小。我通常用兩類方法標定ε一是用歷史數(shù)據(jù)做block bootstrap估計W(Q,μ?_N)的經(jīng)驗分位數(shù)二是把ε當作超參數(shù)在驗證集上做網(wǎng)格搜索。兩類方法結果差別很大但至少能給出一個可解釋的范圍。還有一個關于標準化的細節(jié)需要單獨提出來數(shù)據(jù)標準化不是免費的。先對特征做z-score再畫Wasserstein球等價于把各向同性的球變成橢球。標準化本身沒問題但必須在統(tǒng)一的樣本集上計算均值和方差否則會泄漏未來信息。我習慣在寫完數(shù)據(jù)處理流水線后順手打印一組標準化前后的W1距離對比確認半徑語義沒有漂移。2.3 先用一個最小實驗確認半徑量級在開始求解項目代碼里的大規(guī)模模型之前我會用一個一維合成樣本驗證ε和分布偏移量級的關系。下面的代碼把樣本隨機分半用一維經(jīng)驗分布之間的W1距離近似估計“分布自然偏移”的尺度import numpy as np np.random.seed(7) N 1000 # 風電出力歷史樣本單位 MW xi np.random.gamma(2.0, 15.0, N) xi np.clip(xi, 0, 80) def w1_approx_1d(a, b): # 一維經(jīng)驗分布之間的 W1 近似排序后逐點配對 # 嚴格 W1 應對 CDF 積分這里用于量級估算 n min(len(a), len(b)) return np.abs(np.sort(a)[:n] - np.sort(b)[:n]).mean() # 隨機分半重復 100 次得到經(jīng)驗分布偏移量的分布 offsets [] for _ in range(100): idx np.random.permutation(N) train xi[idx[:N // 2]] valid xi[idx[N // 2:]] offsets.append(w1_approx_1d(train, valid)) offsets np.array(offsets) print(offset 中位數(shù): {:.2f} MW.format(np.median(offsets))) print(offset 90% 分位: {:.2f} MW.format(np.quantile(offsets, 0.9)))這段代碼的核心是先在數(shù)據(jù)上估算“經(jīng)驗分布與經(jīng)驗分布之間”的偏移量級。如果分半樣本之間的W1中位數(shù)是2MW那么ε取1MW可能太緊取10MW又可能過度保守。真實復現(xiàn)時要調(diào)用真正的W1求解器或傳輸距離庫但建模初期用這種快速近似能省掉大量排錯時間。注意這里的w1_approx_1d用排序配對方式近似只適用于一維且長度相近的樣本多維場景不能用這個函數(shù)要回到運輸問題的線性規(guī)劃。3. 從min-max到有限維線性規(guī)劃強對偶與對偶變量的落地3.1 寫出原問題再把它拉回有限維DRO原問題寫作min_{x∈X} sup_{Q∈B_ε(μ?_N)} E_Q[f(x,ξ)]表面上是無窮維問題因為Q是完備概率空間里的分布。當名義分布是經(jīng)驗分布μ?_N(1/N)Σδ_{ξ_i}時可以借Wasserstein球的強對偶把內(nèi)層sup轉(zhuǎn)換掉。標準結論是若地面距離d是范數(shù)f(x,·)對每個x是上半連續(xù)且下有界支撐集Ξ是緊集那么內(nèi)層sup等價于min_{λ≥0} [ λε (1/N) Σ_{i1}^N sup_{ξ∈Ξ} ( f(x,ξ) - λ d(ξ, ξ_i) ) ]其中λ是球約束的解耦乘子。這個形式把“Q到μ?_N距離不超過ε”轉(zhuǎn)成了兩層結構外層是λε罰項內(nèi)層是對每個歷史樣本ξ_i計算它在支撐集上的最壞偏差。λ是優(yōu)化變量而不是超參數(shù)求解過程中會被自動拉到合適位置。ε固定時λ越大相當于對分布偏移的容忍越小模型越保守。這個轉(zhuǎn)換的價值在于原來的變量是無窮維分布現(xiàn)在只需要在有限個歷史樣本上加輔助變量。項目代碼里最能出問題的一步就在這里——對偶引理成立需要緊支撐而很多實現(xiàn)默認ξ落在全空間上導致內(nèi)層sup在目標函數(shù)線性增長時無上界求解器直接報unbounded。3.2 三種典型的f(x,ξ)結構及對應的約束形式f(x,ξ)結構內(nèi)層sup的顯式結果最終問題類型線性f a_x^T ξ b_xsup落在支撐集頂點配合對偶范數(shù)線性規(guī)劃分段線性f max_k (a_k^T ξ b_k)每段取凸包引入輔助變量線性/MILP凸二次f ξ^T A_x ξ 線性項用二次共軛與2-范數(shù)閉式二階錐規(guī)劃線性形式最常用。電力系統(tǒng)調(diào)度里購電成本、棄風懲罰、失負荷懲罰大多隨預測誤差線性增長所以f取關于ξ的線性函數(shù)。內(nèi)層的sup_{ξ} { a_x^T ξ - λ‖ξ - ξ_i‖ }在支撐集是凸多面體時可轉(zhuǎn)成頂點上的線性函數(shù)加范數(shù)錐約束整體保持LP結構。分段線性形式出現(xiàn)在階梯報價、分段補償、以及ReLU形狀的切負荷懲罰里。對偶后每一段都需要一個輔助變量約束數(shù)量約為N×(K1)K是分段數(shù)。如果K很大存儲和求解成本會明顯上升常見做法是把分段函數(shù)改寫成線性目標加額外整數(shù)變量或者用外逼近法逐次加約束。二次形式通常配合2階Wasserstein距離在均值-方差類風險度量里很常見。它的閉式解很漂亮但求解器從LP變成錐規(guī)劃大規(guī)模場景下比線性形式吃力得多。我一般先按線性把流程跑通驗證數(shù)據(jù)流和參數(shù)標定沒問題后再根據(jù)業(yè)務需求升級到分段線性或二次結構。3.3 把一個線性DRO寫成求解器模板下面的代碼是線性DRO對偶形式的通用求解模板用支持集上的網(wǎng)格近似內(nèi)層supimport cvxpy as cp import numpy as np np.random.seed(0) N 100 # 歷史樣本數(shù) K 20 # 支撐集離散點數(shù) xi_hist np.random.randn(N, 2) # 二維不確定性風電和負荷偏差 xi_grid xi_hist[:K] # 支撐集候選點工程上用聚類代表場景 # 決策變量機組出力調(diào)整量 x cp.Variable(2) # 對偶乘子對應 Wasserstein 球約束 lam cp.Variable(nonnegTrue) # 每個歷史樣本的最壞場景偏差 S cp.Variable(N, nonnegTrue) # 模糊集半徑先用 Parameter 方便后續(xù)網(wǎng)格搜索 eps cp.Parameter(nonnegTrue, value0.5) # 線性費用f(x, xi) c^T x d^T xi c np.array([1.0, 1.0]) d np.array([0.8, -0.5]) constraints [] for i in range(N): for k in range(K): dist np.linalg.norm(xi_grid[k] - xi_hist[i]) constraints.append(S[i] d xi_grid[k] - lam * dist) objective eps * lam (1 / N) * cp.sum(S) c x prob cp.Problem(cp.Minimize(objective), constraints) prob.solve() print(最優(yōu)成本: {:.3f}.format(prob.value)) print(對偶乘子 lambda: {:.3f}.format(lam.value))邏輯說明S[i]是第i個歷史樣本周圍最壞場景下的線性費用約束對每個網(wǎng)格點k都要成立才能保證sup被覆蓋。dist是地面距離它直接決定分布偏移的懲罰尺度——如果dist算錯整個切割約束的幾何意義就錯了。eps被定義成Parameter而不是普通常數(shù)這樣之后做ε網(wǎng)格掃描時不用重建整個問題只需循環(huán)賦值再solve。參數(shù)說明N是歷史樣本數(shù)K是支撐集候選數(shù)工程上可以取全部歷史場景或聚類后的代表場景。c和d是費用系數(shù)d對應隨機參數(shù)在目標里的敏感程度這是“預測誤差越大會帶來多大經(jīng)濟損失”的直接體現(xiàn)。這個模板的約束數(shù)量是N×K2000小型問題可以直接求解到N5000、K50時就到25萬條約束必須換算法我在第4章末尾展開講。4. 電力系統(tǒng)場景落地把風電出力不確定性嵌進調(diào)度決策4.1 從點預測到“分布魯棒安全域”風電出力預測誤差不是固定的單峰分布。氣溫、氣壓、風速模式變化會讓預測誤差曲線整體平移或尾部變形SAA在這里經(jīng)常給出“過于自信”的機組啟停方案。DRO的做法預測誤差的歷史樣本構成經(jīng)驗分布以它為中心畫Wasserstein球調(diào)度決策在球內(nèi)最壞分布下仍然要可行。機組組合里的失負荷概率約束比如99%可靠性被改寫成對一個期望懲罰的約束而這個期望是對球內(nèi)的最壞分布取的。這一步不是簡單地把預測區(qū)間放寬。傳統(tǒng)隨機優(yōu)化的置信區(qū)間只覆蓋“參數(shù)取值的范圍”DRO覆蓋的是“分布整體的變化范圍”。所以同樣的調(diào)度成本預算下DRO給的備用容量更平滑不會像逐點魯棒那樣同時把每組出力都釘死在最極端場景上。4.2 兩階段DRO機組組合變量與約束結構階段決策內(nèi)容變量類型約束對象第一階段機組啟停、備用容量二進制/整數(shù)啟停邏輯、最小啟停時間第二階段出力調(diào)整、棄風量、切負荷量連續(xù)功率平衡、爬坡、線路對偶層λ、S_i連續(xù)非負每個歷史樣本的Wasserstein切割兩階段結構里第一階段變量在觀測到風電實際出力之前就要決定第二階段變量在每一組風電場景下分別調(diào)整。用DRO時第二階段約束按場景展開目標里的期望項被對偶后替換成λε加上對每個歷史樣本的S_i。用一個單時段簡化模型演示核心的組裝方式import cvxpy as cp import numpy as np np.random.seed(1) N 200 # 歷史場景數(shù) K 21 # 支撐集網(wǎng)格數(shù) load 100.0 # 負荷MW p_max 120.0 # 機組容量 unit_cost 25.0 # 發(fā)電邊際成本 penalty_loss 500.0 # 切負荷懲罰 # 歷史風電預測誤差單位 MW有正有負 delta_hist np.random.normal(0, 10, N) # 支撐集網(wǎng)格從 -40MW 到 40MW delta_grid np.linspace(-40, 40, K) # 第一階段機組出力 p cp.Variable(nonnegTrue) # 第二階段每個場景的切負荷 l cp.Variable(N, nonnegTrue) # 對偶層變量 lam_dual cp.Variable(nonnegTrue) S cp.Variable(N, nonnegTrue) eps 1.5 # 模糊集半徑MW constraints [] for i in range(N): # 場景功率平衡機組出力 風電實際出力 負荷 - 切負荷 constraints.append(p (50.0 delta_hist[i]) load - l[i]) # Wasserstein 切割對所有支撐點 kS[i] 要覆蓋該點下的最壞損失 for k in range(K): loss_expr penalty_loss * cp.pos(load - p - 50.0 - delta_grid[k]) constraints.append(S[i] loss_expr - lam_dual * abs(delta_grid[k] - delta_hist[i])) objective unit_cost * p eps * lam_dual (1 / N) * cp.sum(S) prob cp.Problem(cp.Minimize(objective), constraints) prob.solve() print(最優(yōu)出力 p {:.1f} MW.format(p.value)) print(lambda {:.2f}.format(lam_dual.value))代碼里的cp.pos是ReLU函數(shù)切負荷量只有在負荷大于機組出力加風電出力時才為正。S[i]必須對每個網(wǎng)格點k都成立它代表的是“在第i個樣本周圍最壞場景下可能產(chǎn)生的切負荷損失”。eps1.5的單位是MW對應第2章里講過的量綱問題。多時段擴展時p變成(G,T)矩陣切割約束按時段分別生成目標里累加S矩陣的和。4.3 求解規(guī)模的現(xiàn)實反饋對偶后的約束數(shù)量大致是N×K乘以時段數(shù)。一個T24、N2000、K30的兩階段機組組合切割約束超過百萬條直接用通用LP求解器會非常吃力。常見做法是Benders分解第一階段整數(shù)變量固定后第二階段DRO子問題是一個純LP每次求解返回一個最優(yōu)性切割給主問題主問題再切一塊可行域。對于支撐集過大導致的切割爆炸我一般先把歷史場景做k-means聚類把K壓到50以內(nèi)聚類中心作為支撐集候選點這樣誤差可控求解規(guī)模也下來一個數(shù)量級。項目代碼里如果直接對原始N×T場景枚舉大概率在小規(guī)模Demo上能跑放大到真實系統(tǒng)就卡死這也是拿到代碼后第一個需要改造的地方。5. 避坑手冊把DRO從“論文能跑”變成“工程能用”5.1 ε不是越大越好調(diào)大半徑容易把模型調(diào)成“破產(chǎn)級保守”現(xiàn)象ε取到樣本標準差的5倍以上目標成本迅速逼近最壞情況魯棒優(yōu)化的水平機組組合結果幾乎不隨場景變化備用容量高到經(jīng)濟上不可行。原因ε描述的是分布偏移的容忍范圍不是安全系數(shù)的代名詞。ε增大后模糊集里納入大量嚴重偏移的分布最壞期望成本主要被這些極端尾部主導調(diào)度方案被迫為完全不合理的大偏差兜底。解決用第2.3節(jié)的分半W1估算先定量綱再用bootstrap取90%分位數(shù)作為ε的參考上限。經(jīng)驗法則上ε與N的-1/d次方同量級N增大時ε要同步縮小不能抱死一個固定值。5.2 W1和W2混用對偶結構變了還不知道現(xiàn)象代碼里只改了一個norm參數(shù)結果求解時間暴漲、數(shù)值穩(wěn)定性變差甚至結果出現(xiàn)明顯偏差。原因W1的切割約束是線性的W2的地面距離平方之后切割里的距離項變成凸二次問題從LP變成二階錐規(guī)劃。很多實現(xiàn)沒有區(qū)分這兩種結構把W1的模板直接套到W2上約束語義已經(jīng)完全不是原問題。解決建模文檔里顯式寫明用的是W1還是W2地面距離的范數(shù)階數(shù)單獨用常量定義。跑通W1之后再升級W2不要在同一個代碼文件里來回改參數(shù)。5.3 支撐集不緊內(nèi)層sup直接無界現(xiàn)象求解器返回unbounded或infeasible檢查約束時發(fā)現(xiàn)S[i]對應的那一組切割約束永遠無法滿足。原因支撐集Ξ沒有加物理邊界。風電出力理論上限、負荷上下限、爬坡約束都沒有寫進支撐集定義線性目標函數(shù)在無界區(qū)域上sup趨向無窮對偶問題自然無解。解決給ξ加箱式邊界或者用截斷距離替代完整范數(shù)。常見做法是ξ∈[ξ_low, ξ_high]網(wǎng)格點從邊界內(nèi)采樣。加了邊界之后對偶問題的緊性條件才成立。5.4 把“分布魯棒”實現(xiàn)成“逐點魯棒”現(xiàn)象回測結果比SAA好但比傳統(tǒng)區(qū)間魯棒好不了太多成本從第2章開始就一路偏高。原因?qū)崿F(xiàn)時把sup_Q E_Q[f]錯誤地實現(xiàn)了max_i f(x, ξ_i)也就是把場景最壞值直接挑出來而不是讓分布Q在Wasserstein球內(nèi)變化。這兩種做法的結果差異巨大逐點魯棒要求每個場景都可接受DRO只要求最壞分布下期望可接受。解決檢查目標函數(shù)里是否還有ε和λ的結構。真正的DRO目標在ε0時退化為SAA在ε較大時才是接近魯棒的如果ε0時結果就和逐點魯棒一樣說明切割約束寫錯了需要回到第3.3節(jié)的模板逐條核對。5.5 標準化時偷看了未來導致離線回測失真現(xiàn)象離線測試集上表現(xiàn)優(yōu)異上線之后連續(xù)幾天偏差放大誤差統(tǒng)計和離線完全對不上。原因標準化用的均值和方差是在全樣本上算的包括驗證集和測試集相當于把未來信息摻進了特征變換。Wasserstein距離對尺度敏感標準化參數(shù)的錯誤會被放大進球半徑語義。解決只用訓練集計算mean和std訓練集、驗證集、測試集分別用同一組參數(shù)變換。并且在代碼里把標準化參數(shù)緩存下來部署時加載同一份參數(shù)避免在線預測時重新擬合分布統(tǒng)計量。6. 上線前的一步驗證用合成分布偏移做回測在把DRO模型放進真實調(diào)度流程之前我會做一個合成驗證專門檢驗“當真實分布真的偏移時DRO是否比SAA更穩(wěn)”。做法是用一個分布生成訓練樣本再用另一個偏移后的分布生成測試樣本看兩種方法在測試樣本上的表現(xiàn)差距。下面這段代碼是評估函數(shù)的核心部分import numpy as np np.random.seed(3) N_train 500 N_test 2000 # 訓練分布中心在 50MW 的偏態(tài)分布 train_wind np.random.beta(2, 5, sizeN_train) * 100 # 測試分布中心偏移 右尾加重 test_wind np.random.beta(2, 4, sizeN_test) * 100 5 def simulate_out_of_sample(decision, wind_sample, load100.0, penalty500.0): shortages np.maximum(0.0, load - decision - wind_sample) return penalty * shortages.mean()分別用SAA和DRO求兩個決策SAA決策用訓練樣本直接優(yōu)化DRO決策用第4章的模型在ε1.5下求解然后都放到test_wind上評估。典型結果趨勢如下ε訓練損失測試損失最壞損失(測試)0SAA8.241.789.30.59.124.655.21.511.319.841.04.018.921.438.6解讀這個表ε0的SAA在訓練分布上最便宜測試分布一偏移就崩ε1.5時訓練成本略升但測試損失大幅下降ε4.0時訓練成本明顯偏高測試損失卻只降了一點點。這說明存在一個合理ε區(qū)間過低不夠穩(wěn)過高過度保守。那張風險-成本曲線就是上線與否的決策依據(jù)。從那以后我每次調(diào)ε都強制先畫一遍風險-成本曲線再談能不能上線。這個習慣幫我避開了至少三次“論文指標漂亮、現(xiàn)場一跑就崩”的尷尬。希望幫到你。本文還有配套的精品資源點擊獲取