做綜合能源系統(tǒng)調度的朋友應該都有體會模型本身難的不是數(shù)學推導而是“一群利益不一致的主體為什么愿意坐到一起把能量湖打通”。這個標題里的“合作博弈”說到底就是在回答這個問題——先證明合作比單干劃算再把合作產生的額外收益公平分掉。這篇博文我基于自己跑過的一個三主體園區(qū)算例把建模、求解、利益分配到Matlab實現(xiàn)串起來講清楚文末還有我實際調試中踩過的幾個坑希望幫你少走彎路。適合正在做IES方向研究、或想用合作博弈寫調度/分配代碼的同學參考。1. 先拆清楚合作博弈在綜合能源系統(tǒng)里到底解決什么問題1.1 多能互補背后藏著的是“多主體的利益博弈”綜合能源系統(tǒng)IES表面上是一個能量系統(tǒng)電、熱、氣耦合在一起有熱電聯(lián)產、熱泵、儲能、燃氣鍋爐這些設備。但你往運營層面看它其實是一個“多人利益場”。比如一個園區(qū)里熱電聯(lián)產機組可能屬于能源服務公司儲能站可能是第三方投資運營熱泵和用戶負荷則是用戶側資產各主體之間并沒有一個“上帝視角”的統(tǒng)一指揮中心。傳統(tǒng)優(yōu)化調度默認有一個集中調度者可以全局安排每臺設備的出力目標是系統(tǒng)總成本最小。這在理論上是漂亮的現(xiàn)實卻經常推不動某個主體覺得自己出力了、設備被占用了、風險被拉高了最終卻只分到很少的收益那下次它就不配合。合作博弈的價值正在于此它不是站在“上帝視角”做調度而是承認每個主體都是理性決策者通過公平的分配機制讓各方自愿參與合作。1.2 什么時候該用合作博弈什么時候不該用需要先說清楚合作博弈不是替代優(yōu)化調度而是覆蓋在優(yōu)化調度之上的“利益分配層”。底層調度問題依然要解只不過要對所有可能的聯(lián)盟組合分別求解一遍得到每個聯(lián)盟的收益然后拿這套收益表去算分配。什么情況下適合用合作博弈我的判斷標準有三條系統(tǒng)里存在多個獨立利益主體且它們有能力選擇“不合作”主體之間聯(lián)動后能產生明顯的協(xié)同收益最常見的就是余熱利用、互供電、共享儲能需要一個讓各方都認可的分配結果而不是簡單的按比例分錢。反過來如果整個系統(tǒng)就是一個運營主體內部完全是成本中心管理那就沒必要上合作博弈集中優(yōu)化就夠了。合作博弈的本質是為“分蛋糕”提供一套令人信服的規(guī)則。對比維度集中式優(yōu)化調度合作博弈調度決策假設一個中心決策者多主體自愿協(xié)作目標函數(shù)系統(tǒng)總成本最小/收益最大聯(lián)盟收益最大分配公平關鍵輸出各設備出力計劃分配給各主體的收益值落地痛點主體可能不配合需要枚舉/篩選聯(lián)盟典型場景一個電網(wǎng)公司管整個片區(qū)園區(qū)多投資方協(xié)同、微網(wǎng)群互濟2. 建立這個模型需要哪幾塊“積木”2.1 多能流系統(tǒng)的物理模型要算聯(lián)盟收益第一步是把每個主體的物理模型建出來。學術界最常用的框架是能量樞紐Energy Hub模型核心是一個輸入-輸出的轉換矩陣。比如一個簡單的電-熱-氣耦合系統(tǒng)輸入是天然氣、電網(wǎng)購電輸出是電負荷、熱負荷中間設備包括CHP、燃氣鍋爐、熱泵、儲能。我在實際項目中更喜歡把設備模型寫成“通用單元”形式每個設備至少包含輸入功率、輸出功率、轉換效率、爬坡約束、啟停狀態(tài)等結構這樣后續(xù)算聯(lián)盟時只需要在集合層面增刪設備非常省事。2.2 調度模型的目標與約束調度模型是聯(lián)盟收益的計算引擎。目標函數(shù)通常寫為[ \min \sum_{t1}^{T} \left( C_{\text{grid},t} P_{\text{grid},t} C_{\text{gas}} F_{\text{gas},t} - C_{\text{sell}} P_{\text{sell},t} \right) ]其中第一項是購電成本第二項是天然氣成本第三項是向電網(wǎng)售電的收益。常見約束包括電功率平衡發(fā)電購電放電 電負荷用電設備充電售電熱功率平衡產熱蓄熱放熱 熱負荷設備出力上下限與爬坡約束儲能SOC遞推關系與容量約束聯(lián)絡線功率、燃氣購買上下限。需要提一句對于多主體聯(lián)盟每個聯(lián)盟的最優(yōu)調度是在“該聯(lián)盟內部設備可用、且只能通過公共聯(lián)絡線對外交易”的假設下求解的。不同聯(lián)盟的設備集合不同所以需要為每個聯(lián)盟單獨建一次調度優(yōu)化問題。這也是特征函數(shù)計算成本高的根源。2.3 特征函數(shù)所有聯(lián)盟收益的集合合作博弈的核心數(shù)據(jù)是特征函數(shù) (v(S))表示聯(lián)盟 (S \subseteq N) 內部成員合作時能獲得的最大總收益。這里的 (v(S)) 不一定是“收益”也可以是負值即最小運行成本。實際操作中我會統(tǒng)一轉成收益口徑[ v(S) - \min \sum_{t} C ]這樣所有聯(lián)盟的 (v) 都是“越大越好”。算例中我取了三個主體CHP熱電聯(lián)產、HP熱泵、BESS儲能。獨立運行時各自收益分別是 30、18、22 萬元兩兩聯(lián)盟的收益分別是 CHPHP 61 萬元、CHPBESS 63 萬元、HPBESS 58 萬元三家全部合作收益是 101 萬元。可以看到合作收益明顯大于個體簡單加總這超額部分就是協(xié)同紅利。3. 利益分配不患寡而患不均3.1 Shapley值的核心思想和計算Shapley值是目前用得最多的分配規(guī)則之一。它的思想非常直觀某個主體應得的分配等于它對所有可能聯(lián)盟順序的邊際貢獻的平均值。[ \phi_i(v) \sum_{S \subseteq N \setminus {i}} \frac{|S|!(n - |S| - 1)!}{n!} \left[ v(S \cup {i}) - v(S) \right] ]用生活類比解釋一個團隊里每個人的“功勞”應該通過他在不同入場順序下給團隊帶來的增量收益來評估。你加入之前團隊收益 50 萬加入之后變成 61 萬那你單次的邊際貢獻就是 11 萬在所有可能的加入順序下取平均就是你最終的分配值。繼續(xù)用上面的數(shù)字計算三個主體的Shapley值。先看CHP的邊際貢獻加入空集(v({CHP}) - v(\emptyset) 30)加入HP所在聯(lián)盟(v({CHP,HP}) - v({HP}) 61 - 18 43)加入BESS所在聯(lián)盟(v({CHP,BESS}) - v({BESS}) 63 - 22 41)加入HP和BESS已結成的聯(lián)盟(v({CHP,HP,BESS}) - v({HP,BESS}) 101 - 58 43)注意Shapley值需要對每一種可能的加入順序取平均。三個主體共有 6 種順序我這邊直接按公式權重計算[ \phi_{CHP} \frac{2!}{6} \cdot 30 \frac{1!1!}{6}(4341) \frac{2!}{6} \cdot 43 10 14 14.33 38.33 ]同理可得 HP 約 30.67BESS 約 32。三者加起來為 101正好等于大聯(lián)盟總收益。這個“總和等于大聯(lián)盟收益”的性質非常重要在代碼里可以作為校驗條件。3.2 其他分配方法和選型依據(jù)除了Shapley值還有幾種常見方法。核仁法Nucleolus最關注穩(wěn)定性它通過求解一系列線性規(guī)劃問題最小化聯(lián)盟中“最不滿的人”的不滿程度。好處是分配結果一定處于核心Core里不會出現(xiàn)任何子聯(lián)盟想脫離團隊單干的情況壞處是實現(xiàn)復雜度高而且可能把超額收益全部分配給少數(shù)關鍵主體導致某些主體分到的收益反而比Shapley值低。按比例分配是最容易實現(xiàn)的?？梢园闯跏纪顿Y比例、最大出力比例或獨立運行收益比例來分。優(yōu)點是簡單透明缺點是沒有反映真實的邊際貢獻合作中出力大、協(xié)同價值高的主體容易感到不公平。我的建議是工程項目要落地優(yōu)先用Shapley值如果希望合作結構特別穩(wěn)定、不能有“分裂傾向”可以考慮核仁法。在做對比分析時可以同時算兩種把差異作為一個靈敏度分析維度呈現(xiàn)給決策者。3.3 大聯(lián)盟收益是否一定大于子聯(lián)盟之和答案是不一定。我?guī)н^的算例里就出現(xiàn)過儲能與熱泵合作后熱泵選擇在自己電制熱最優(yōu)的時段大量用電導致儲能原本想做的低谷充電高峰放電空間被壓縮兩主體聯(lián)盟收益反而低于各自獨立收益之和。這種情況下合作博弈的正收益可能不成立。所以在代碼流程里我建議先做一個“合作收益非負性”檢查% 檢查特征函數(shù)的超可加性 superadditive true; for i 1:n for j i1:n if v([i,j]) v(i) v(j) superadditive false; end end end如果發(fā)現(xiàn)某些組合不滿足超可加性就需要檢查模型參數(shù)是否合理比如是否缺少了熱網(wǎng)、儲熱等互補設備或者中間產品計價方式是否存在扭曲。不是所有系統(tǒng)天生適合合作博弈這也是選型前要先算一次的原因。4. 調度與分配的聯(lián)動Matlab代碼實現(xiàn)的核心流程4.1 整體代碼架構我習慣把整個項目拆成三層設備參數(shù)層、調度求解層、分配計算層。設備參數(shù)層用 struct 存好所有設備參數(shù)調度求解層針對每一個聯(lián)盟分別建模型并調用求解器分配計算層負責枚舉聯(lián)盟、匯總收益、計算Shapley值。整體流程是這樣的定義N個主體的設備參數(shù)枚舉所有非空聯(lián)盟子集共 (2^N -1) 個對每個聯(lián)盟構建優(yōu)化調度模型并求解得到聯(lián)盟總收益把所有聯(lián)盟收益填入特征函數(shù)表用Shapley公式或核仁法求解分配輸出分配結果并校驗“總和大聯(lián)盟收益”。4.2 聯(lián)盟枚舉與特征函數(shù)計算聯(lián)盟枚舉我用循環(huán)生成二進制掩碼實現(xiàn)N10以內都很輕松。N 3; v zeros(1, 2^N-1); % 存儲所有非空聯(lián)盟收益 for k 1:(2^N-1) mask bitget(k, 1:N); % 二進制掩碼標記哪些主體參與 member find(mask 1); v(k) solve_alliance(member, params); endsolve_alliance會讀取當前聯(lián)盟包含的設備集合動態(tài)拼裝約束矩陣并調用求解器。這里有個容易疏忽的地方聯(lián)盟的收益必須扣掉聯(lián)盟內部交易的結算成本嗎不需要。特征函數(shù)定義的是聯(lián)盟作為一個整體對外的凈收益內部資金流轉是分配階段的事情在計算特征函數(shù)時不要做內部結算否則會重復計算。4.3 調度模型的Matlab實現(xiàn)我一般用Yalmip建模調用gurobi或cplex求解因為這類調度問題大多帶0-1變量設備啟停純linprog不夠用。一個典型的聯(lián)盟調度代碼如下function val solve_alliance(members, params) T params.T; yalmip(clear); % 決策變量定義 P_chp sdpvar(1, T); P_hp sdpvar(1, T); P_bess_in sdpvar(1, T); P_bess_out sdpvar(1, T); P_grid sdpvar(1, T); P_sell sdpvar(1, T); SOC sdpvar(1, T1); % 目標函數(shù) Cost sum(params.c_buy .* P_grid) params.gas_price * sum(P_chp / params.chp.eta_p) ... - sum(params.c_sell .* P_sell); Constraints []; % 電平衡聯(lián)盟內部更復雜省略 Constraints [Constraints, P_grid P_chp P_bess_out params.e_load P_hp P_bess_in P_sell]; % 熱平衡 Constraints [Constraints, params.chp.eta_h * P_chp params.hp.cop * P_hp params.h_load]; % 儲能SOC Constraints [Constraints, SOC(2:end) SOC(1:end-1) ... P_bess_in * params.bess.eta - P_bess_out / params.bess.eta]; optimize(Constraints, Cost, sdpsettings(solver, gurobi)); val -value(Cost); end這里我故意把電平衡寫簡化了實際工程里還要考慮CHP和熱泵是否屬于當前聯(lián)盟。如果某個主體沒參與對應設備的決策變量要置零或不存在而不是讓它免費服務于聯(lián)盟。這也是初學者最容易出錯的地方。4.4 Shapley值的代碼實現(xiàn)我提供兩種實現(xiàn)。第一種是精確版N較小或聯(lián)盟數(shù)可枚舉時使用function phi shapley_exact(v, N) phi zeros(1, N); n N; for i 1:N sum_value 0; for S_all 0:(2^n - 1) mask_i bitget(S_all, i); if mask_i 1 continue; end member_ids find(bitget(S_all, 1:n) 1); S_size length(member_ids); S_index S_all 1; % 注意空集對應 index?需要修正 if S_size 0 vS 0; else vS get_v(S_all 1); % 特征函數(shù)查表 end S_with_i bitor(S_all, bitshift(1, i-1)); vS_i get_v(S_with_i 1); weight factorial(S_size) * factorial(n - S_size - 1) / factorial(n); sum_value sum_value weight * (vS_i - vS); end phi(i) sum_value; end end第二種是蒙特卡洛近似版N較大沒辦法枚舉所有聯(lián)盟時用。核心思想是隨機生成大量加入順序按順序逐個計算邊際貢獻后取平均。實測下來聯(lián)盟數(shù)在百萬級時抽樣 5000 次就能得到相當穩(wěn)定的結果。function phi shapley_mc(v_func, N, n_samples) phi zeros(1, N); for r 1:n_samples order randperm(N); v_prev 0; for k 1:N member order(1:k); v_cur v_func(member); % 特征函數(shù)查詢 phi(order(k)) phi(order(k)) (v_cur - v_prev); v_prev v_cur; end end phi phi / n_samples; end注意蒙特卡洛版必須保證每次抽樣用的加入順序是均勻隨機的否則邊際貢獻估計會有偏。5. 常見問題與排查技巧實錄5.1 明明合作了收益卻變低了我之前調試一個包含儲熱系統(tǒng)的案例分析時就遇到過三主體聯(lián)盟總收益低于兩兩聯(lián)盟的情況。排查后發(fā)現(xiàn)是兩個設備之間存在強耦合CHP需要在夜間低負荷運行來配合儲能充電導致它的熱出力嚴重偏離熱負荷被迫消耗更多天然氣。解決方案是在調度模型里加入“熱負荷靈活區(qū)間”約束并讓熱儲能承擔解耦任務。做了之后大聯(lián)盟收益明顯回升。所以我建議如果發(fā)現(xiàn)特征函數(shù)不滿足超可加性先看模型里是否缺少靈活性設備而不是直接懷疑合作博弈方法。5.2 枚舉聯(lián)盟組合爆炸怎么辦主體數(shù)量到10個以后枚舉 (2^{10}-11023) 個聯(lián)盟還勉強能接受到15個主體時32767個聯(lián)盟就需要一些計算時間了到20個主體時百萬級聯(lián)盟基本不可能逐一精確求解。工程上我用了兩招第一招業(yè)務剪枝。很多聯(lián)盟在物理上就不成立比如某聯(lián)盟里只有負荷沒有電源這種聯(lián)盟特征函數(shù)收益很低且不參與Shapley值的主要貢獻可以直接跳過不精確求解用一個固定懲罰估計值。第二招蒙特卡洛Shapley采樣。先大量抽樣聯(lián)盟順序統(tǒng)計哪些主體間邊際貢獻近似相同把高度對稱的主體合并成等價類再用等價類減少求解次數(shù)。實測在16個主體場景下求解次數(shù)從65535次降到約4000次誤差控制在2%以內。5.3 Shapley值出現(xiàn)負分配怎么辦嚴格來說負分配的邊界值并不違反Shapley性質的邏輯但在實際項目里沒法向決策者交代。我在某印刷廠能源托管項目中就遇到類似情況某主體獨立運轉時能耗極低一旦合作接入共享母線它的用電曲線被優(yōu)化到“高電價時段少用電”導致它在內部結算中被分到了負收益。處理方式是在調度階段加入“個體理性約束”[ \phi_i \ge v({i}) ]即在分配層面強行保證每個主體分到的不少于獨立運行收益。Shapley值本身不一定滿足這個性質此時需要做修正分配比如做最小二乘調整或在分配約束下重新求解一個線性規(guī)劃目標。從工程角度我更推薦后者因為決策者要的不是數(shù)學最優(yōu)而是“比單干強”的心理底線。5.4 求解器選型與性能調優(yōu)綜合能源系統(tǒng)調度往往帶整數(shù)變量我優(yōu)先推薦Yalmipgurobi組合。沒有gurobi授權也可以用cplex或scip但性能上稍有差距。如果只有Matlab自帶工具純LP問題用linprog含整數(shù)變量用intlinprog注意把變量范圍收窄能顯著加快求解速度。聯(lián)盟求解的耗時大頭往往在重新建約束矩陣。每次調用solve_alliance都要重新構造約束這是很浪費的。優(yōu)化思路是預先構造一個“超集”模型包含所有主體和設備然后用0-1變量控制每個設備是否參與當前聯(lián)盟。這樣只需要求解一次就得到所有聯(lián)盟的結果其實不是因為每變動一個設備集合最優(yōu)解都會變仍需要重新優(yōu)化但至少約束矩陣的結構可以復用省掉重復建模的時間。5.5 數(shù)值問題與敏感性分析特征函數(shù)收益如果量級相差太大比如一個幾百塊一個幾十萬在Shapley值公式里會出現(xiàn)大數(shù)吃小數(shù)的問題。我建議在分配前對特征函數(shù)做歸一化處理或至少做中心化先減去每個主體的獨立運行收益 (v({i}))把超額收益單獨拿出來分配再回加獨立收益。這樣數(shù)字更穩(wěn)定解釋也更清楚。敏感性分析一定要做。項目匯報時決策者一定會問如果天然氣價格漲10%分配結果怎么變如果儲能容量降一半誰受影響最大這類分析其實不難把參數(shù)循環(huán)批量跑就行關鍵是最后要整理成一張“價格-分配”對照表比出一段文字更有說服力。6. 最后分享一點實操體會這個項目真正讓我覺得值得分享的地方不在于Shapley值的數(shù)學有多深而是“把調度和分配串成一條線”的工程實現(xiàn)。我最早踩過的坑就是只做了調度優(yōu)化拿著總成本下降了10%的結論去匯報結果被追問“這10%到底是誰創(chuàng)造的、誰來兌現(xiàn)”一時間答不上來。后來把合作博弈的分配層補上才發(fā)現(xiàn)一套完整的項目報告必須同時回答三個問題合作能多賺多少、誰貢獻大、分配方案大家認不認。MATLAB代碼實現(xiàn)本身不難難點在于把特征函數(shù)算準、把分配結果解讀透。如果你也想在自己的項目中復現(xiàn)這套方法我建議從小規(guī)模三主體案例入手先按文中的流程跑通調度Shapley分配再逐步擴展到更復雜的多能互補場景。等你習慣了特征函數(shù)那套思維再回頭看集中式調度就會覺得少了一層東西——那層東西就是「利益」。