主從博弈的MATLAB實現(xiàn)與優(yōu)化調(diào)度實戰(zhàn))
綜合能源微網(wǎng)、共享儲能、主從博弈這三個詞放在一起很容易讓人覺得這是某篇論文里的“高級概念”但真正動手把MATLAB代碼跑起來之后我才發(fā)現(xiàn)這套模型并沒有想象中那么遙不可及。它本質(zhì)上解決的是一個非常實際的問題儲能設(shè)備到底應(yīng)該歸誰用、怎么定價、微網(wǎng)怎么根據(jù)價格來安排自己的用電計劃。我去年幫一個課題組搭過一套“共享儲能綜合能源微網(wǎng)”的仿真算例用的就是主從博弈框架MATLAB里面既有YALMIP建模也有迭代求解的邏輯。今天這篇想把整個代碼架構(gòu)和實現(xiàn)思路拆開聊一遍適合正在做微網(wǎng)優(yōu)化、儲能商業(yè)模式或者主從博弈仿真的研究生也適合剛接觸能量管理系統(tǒng)的工程師至少能讓你少走不少彎路。1. 先理清思路主從博弈為什么適合“共享儲能微網(wǎng)”這類場景很多剛接觸這個方向的人第一問題是“為什么非得用主從博弈直接用集中式優(yōu)化不行嗎”這個問題問得很關(guān)鍵因為只有把建模動機想明白后面看代碼才不會暈。1.1 微網(wǎng)優(yōu)化運行到底在優(yōu)化什么綜合能源微網(wǎng)通常包含風機、光伏、燃氣輪機、電鍋爐、儲能、熱負荷和電負荷等元素。所謂優(yōu)化運行就是在一個調(diào)度周期內(nèi)決定每一臺設(shè)備每個時刻的出力使得總運行成本最低同時滿足電、熱、氣等多種能源的平衡。常見的調(diào)度周期是24小時步長為1小時也可以用到15分鐘。目標是經(jīng)濟性約束是能量平衡和設(shè)備物理限制。這個框架本身并不復(fù)雜核心難點在于“多能互補”讓變量數(shù)量和耦合關(guān)系成倍增加。比如電轉(zhuǎn)熱會讓電鍋爐消耗電能產(chǎn)出熱能燃氣輪機同時產(chǎn)電和產(chǎn)熱儲能裝置的充放電狀態(tài)又會影響下一時刻的SOC這種時間耦合讓優(yōu)化問題變成一個動態(tài)決策問題。1.2 共享儲能解決了誰的痛點如果每個微網(wǎng)都自建一套儲能系統(tǒng)容量利用率通常不高。風電光伏大發(fā)的時段可能用不上晚上負荷高峰又未必充得夠。儲能成本擺在那里電池壽命、維護費用、閑置損耗都是實打?qū)嵉馁Y金壓力。共享儲能的概念就是把儲能設(shè)施交給獨立的運營商建設(shè)并運營微網(wǎng)按需購買充放電功率運營商通過服務(wù)費回收成本并盈利。這個模式有點像共享充電寶用戶不用買充電寶按次付費平臺通過調(diào)度充電寶在不同需求點之間的流轉(zhuǎn)來提高利用率。儲能運營商通過價格信號引導(dǎo)用戶在不同時段充放電既能賺取峰谷價差又能幫助微網(wǎng)降低用電成本。1.3 主從博弈的上下層關(guān)系主從博弈也叫Stackelberg博弈是一種不對稱決策模型。在這個場景里儲能運營商是領(lǐng)導(dǎo)者先制定充放電價微網(wǎng)是跟隨者在給定的價格下優(yōu)化自己的用能和儲能購買計劃。為什么儲能運營商是領(lǐng)導(dǎo)者因為它擁有儲能資源的定價權(quán)處于更主動的位置。微網(wǎng)雖然也做決策但只能被動響應(yīng)價格。換句話說上層先出價下層再根據(jù)價格選最優(yōu)策略然后上層又根據(jù)下層的反應(yīng)調(diào)整價格直到雙方都沒有動力改變自己的策略就達到了博弈均衡。1.4 為什么選“領(lǐng)導(dǎo)者-跟隨者”而非集中式優(yōu)化集中式優(yōu)化的前提是整個系統(tǒng)的所有信息都集中在一個調(diào)度中心手里統(tǒng)一優(yōu)化所有參與者的目標實現(xiàn)全局成本最低。但實際情況是微網(wǎng)和儲能運營商屬于不同主體各有各的私密信息比如微網(wǎng)的負荷預(yù)測、設(shè)備效率、內(nèi)部成本運營商很難完全掌握反過來也一樣。主從博弈的好處是只需要交換價格和功率計劃不需要把各自的內(nèi)部模型完全公開。這種分布式?jīng)Q策的方式更貼近電力市場的真實交易機制。代價是求解變難因為不是一個單層優(yōu)化問題而是上下兩層相互嵌套的決策問題。2. MATLAB代碼架構(gòu)與核心建模我寫這套MATLAB代碼的時候最初目標很明確用一個雙層模型分別描述儲能運營商的定價決策和微網(wǎng)的運行調(diào)度再用迭代或者KKT化簡的方式去求解納什均衡。整個代碼要方便調(diào)試也要方便擴展。2.1 上層模型儲能運營商收益最大儲能運營商的決策變量是分時充電電價、放電電價以及服務(wù)費。它的收益來自多個方面向微網(wǎng)出售電能獲得售電收入向微網(wǎng)收取充放電服務(wù)費從電網(wǎng)或上級能源市場低價購電用于給儲能充電。目標函數(shù)就是總收益減去總成本。成本包括從電網(wǎng)購電的成本、儲能充放電過程中的損耗成本、運維成本。運維成本通常簡化成一個線性函數(shù)比如每次充放單位電量的成本為固定值不考慮電池老化或循環(huán)壽命的非線性影響。約束條件主要包括儲能功率上下限同一時刻的充電功率和放電功率不能超過額定值儲能容量約束SOC保持在安全范圍通常設(shè)為0.1到0.9SOC動態(tài)方程SOC(t) SOC(t-1) 充電功率×充電效率/容量 - 放電功率/(放電效率×容量)不能同時充放電售電價格通常會有上下限約束防止惡意抬價。2.2 下層模型微網(wǎng)運行成本最小微網(wǎng)作為跟隨者看到儲能運營商給出的充放電價后需要安排自己內(nèi)部的機組出力和從儲能購買的功率使得當天總運行成本最低。成本主要包括從電網(wǎng)購電費用燃氣輪機消耗天然氣的成本從共享儲能購買充電服務(wù)的費用如果儲能允許微網(wǎng)充電或購買放電電能的費用如果微網(wǎng)有余電也可以反賣給儲能或電網(wǎng)形成負成本。約束條件包括電功率平衡風機出力光伏出力燃氣輪機發(fā)電儲能放電電網(wǎng)購電 電負荷儲能充電電鍋爐耗電熱功率平衡燃氣輪機余熱電鍋爐產(chǎn)熱儲熱設(shè)備放熱 熱負荷燃氣輪機出力上下限和爬坡約束電鍋爐出力上下限電網(wǎng)交互功率上下限。這里需要注意綜合能源微網(wǎng)的下層模型本身就是一個多能源耦合的線性規(guī)劃問題。如果儲能價格是固定的那么下層模型可以直接用線性規(guī)劃求解非常方便。2.3 博弈均衡求解的兩種MATLAB思路主從博弈的求解目前主流有兩種路徑我兩種都試過各有適用范圍。第一種是“KKT單層化”。將下層微網(wǎng)優(yōu)化問題的KKT最優(yōu)性條件作為約束加到上層優(yōu)化問題中。這樣做的好處是理論上可以得到精確解但代價是引入了大量互補松弛條件和非線性項通常需要用大M法線性化變成一個混合整數(shù)線性規(guī)劃問題對求解器要求比較高。YALMIP配合CPLEX或者Gurobi能處理但模型規(guī)模大時變量和約束數(shù)量會爆炸式增長。第二種是“迭代求解”。上層給定價格下層求最優(yōu)功率計劃并返回上層根據(jù)返回結(jié)果更新價格循環(huán)直到價格和收益不再變化。這種方法實現(xiàn)簡單能直觀看到博弈過程但收斂性不能保證需要加阻尼或調(diào)整步長。我實際演示時用迭代法更多因為它更直觀也更容易讓初學者理解博弈的交互過程。2.4 文件組織與數(shù)據(jù)流我習慣把代碼拆成幾個模塊main.m主程序負責初始化參數(shù)、調(diào)用子模塊、循環(huán)迭代、結(jié)果展示init_parameters.m設(shè)置所有參數(shù)包括設(shè)備容量、負荷數(shù)據(jù)、電價、效率等build_microgrid_model.m構(gòu)建下層微網(wǎng)優(yōu)化模型并求解輸入是儲能價格輸出是各設(shè)備功率和微網(wǎng)購電/購儲能計劃build_storage_model.m計算上層儲能運營商的收益輸入是微網(wǎng)的購儲能功率輸出是收益值update_price.m根據(jù)上下層偏差更新價格check_convergence.m判斷是否達到收斂條件。數(shù)據(jù)流是主程序初始化儲能價格然后把價格傳給微網(wǎng)模型微網(wǎng)求解放電功率和充電功率再把這部分功率返回給上層上層更新價格直到收斂。這個結(jié)構(gòu)非常簡單也方便后續(xù)把某個模塊替換成更復(fù)雜的實例。3. 核心代碼模塊逐個拆解代碼到底怎么寫很多人喜歡直接上網(wǎng)找開源項目但開源項目往往為了通用性做了大量封裝反而不容易看懂。我更推薦自己一步一步把模塊搭起來每個模塊都清楚它在算什么。3.1 數(shù)據(jù)準備模塊的細節(jié)先在init_parameters.m里把基礎(chǔ)數(shù)據(jù)設(shè)置好。類別有參數(shù)、單位、數(shù)值。可以用結(jié)構(gòu)體封裝方便后續(xù)調(diào)用。比如設(shè)置微網(wǎng)的原始負荷、風電出力、光伏出力每個都是一個24維的列向量。常見做法是讀Excel或CSV也可以直接在腳本里用數(shù)組賦值。如果做測試我建議先用一個合成數(shù)據(jù)比如風電出力在夜間低谷大、白天小光照在中午大負荷早晚高峰這樣方便驗證模型行為是否合理。我常用的代碼結(jié)構(gòu)是這樣的% 初始化參數(shù) T 24; % 調(diào)度周期單位:小時 P_wind [0.4; 0.35; 0.3; ...]; % 風電預(yù)測出力 P_pv [0; 0; 0; ...]; % 光伏預(yù)測出力 P_load [0.8; 0.75; ...]; % 電負荷 H_load [0.6; 0.6; ...]; % 熱負荷 price_grid [0.38; 0.38; ...]; % 分時電網(wǎng)電價每個數(shù)據(jù)都要保持相同的量綱。我在實際調(diào)參時把所有功率統(tǒng)一成MW所有電量統(tǒng)一成MWh價格統(tǒng)一成元/MWh雖然看著數(shù)字很大但換算關(guān)系比混著用千瓦時和兆瓦時省心得多。3.2 主循環(huán)與價格更新框架主程序的核心邏輯是一個循環(huán)里面反復(fù)調(diào)用兩個模塊。偽代碼思路如下% 主循環(huán)迭代 price_c 0.3 * ones(T,1); % 初始充電價 price_d 0.5 * ones(T,1); % 初始放電價 lambda 0.5; % 迭代步長 max_iter 50; for iter 1:max_iter % 下層求解微網(wǎng)調(diào)度 [P_dis, P_char] solve_microgrid(price_c, price_d); % 上層根據(jù)微網(wǎng)返回的充放電功率更新價格 [new_price_c, new_price_d] update_price(price_c, price_d, P_dis, P_char); % 判斷收斂 if max(abs(new_price_c - price_c)) 1e-4 max(abs(new_price_d - price_d)) 1e-4 break; end price_c new_price_c; price_d new_price_d; end這個過程里微網(wǎng)返回的充電功率和放電功率是關(guān)鍵的銜接變量。微網(wǎng)如果大量買放電功率說明當前放電價便宜上層下輪可以適當抬價如果放電功率很小說明微網(wǎng)不太買賬上層就得降價。這種調(diào)整邏輯很像是市場里的價格試探。3.3 下層微網(wǎng)調(diào)度YALMIP實現(xiàn)微網(wǎng)模型我用YALMIP寫因為YALMIP的語法更接近數(shù)學表達式容易檢查模型對不對。下面是一個簡化版的微網(wǎng)調(diào)度模型片段。% 變量定義 P_gt sdpvar(T,1); % 燃氣輪機出力 P_eb sdpvar(T,1); % 電鍋爐功率 P_grid sdpvar(T,1); % 電網(wǎng)購電功率 P_s_dis sdpvar(T,1); % 從共享儲能放電獲取的功率 P_s_char sdpvar(T,1); % 給共享儲能充電的功率 SOC_s sdpvar(T,1); % 微網(wǎng)內(nèi)部儲能SOC如果有 % 目標函數(shù) objective sum(price_grid .* P_grid) sum(c_gas * P_gt) ... sum(price_d .* P_s_dis) - sum(price_c .* P_s_char); % 約束 constraints []; constraints [constraints, P_load P_wind P_pv P_gt P_grid P_s_dis - P_eb - P_s_char]; constraints [constraints, 0 P_gt P_gt_max]; % ... 其他約束 options sdpsettings(solver, gurobi, verbose, 0); optimize(constraints, objective, options);這里有一個容易混淆的地方微網(wǎng)向共享儲能買電時是從儲能放電側(cè)買所以價格用price_d如果微網(wǎng)自己的新能源發(fā)多了可以把多余的電存到共享儲能里獲得充電服務(wù)收入價格用price_c。因為price_c是儲能運營商設(shè)定的充電價微網(wǎng)存進去越多運營商賺的服務(wù)費越多但微網(wǎng)獲得的電價收入可能是負的即微網(wǎng)需要支付一定的充電服務(wù)費。3.4 上層儲能運營商收益計算當微網(wǎng)返回它需要的充放電功率后上層就能計算儲能運營商的收益。收益表達為% 儲能運營商收益 income_sell sum(price_d .* P_s_dis); % 放電收入 income_char sum(price_c .* P_s_char); % 充電服務(wù)收入 cost_purch sum(price_grid_purchase .* (P_s_char P_s_dis)); % 從電網(wǎng)購電成本 cost_om sum(om_cost * (P_s_char P_s_dis)); % 運維成本 profit_storage income_sell income_char - cost_purch - cost_om;這里注意儲能運營商的購電成本要覆蓋兩方面一是微網(wǎng)給儲能充電的功率儲能運營商相當于收了微網(wǎng)的電卻要自己掏錢買電補充這個邏輯需要定義清楚。實際建模中儲能運營商從電網(wǎng)購電用于給儲能充電再放電賣給微網(wǎng)。如果微網(wǎng)自己有多余電存進來就是另一種交易模式。不同假設(shè)會導(dǎo)致收益公式略有不同我在仿真時會把交易協(xié)議寫在注釋里避免自己過兩天忘記。3.5 收斂判斷和結(jié)果輸出收斂判據(jù)是用相鄰兩次迭代的價格差和收益差同時判斷。因為價格更新時可能有小幅度震蕩單看價格可能不夠配合收益差更穩(wěn)。結(jié)果輸出建議畫四張圖微網(wǎng)電功率平衡圖展示各個電源出力和負荷曲線儲能充放電圖展示SOC和充放電功率博弈迭代收斂圖橫軸迭代次數(shù)縱軸價格或收益微網(wǎng)購電費用和儲能收益柱狀圖。畫圖用plot和stairs就行關(guān)鍵是圖例和單位要做清楚。我經(jīng)常犯的錯是忘記把MWh轉(zhuǎn)成kWh結(jié)果圖上的數(shù)值和PPT里的對不上。4. 調(diào)參與踩坑記錄把我的實戰(zhàn)經(jīng)驗直接給你這部分是重點。模型跑通容易跑得合理很難我整理一些自己實際遇到過的坑。4.1 求解器選擇與非線性項處理如果下層模型是線性的用YALMIP配Gurobi或CPLEX最穩(wěn)。MATLAB自帶的linprog也行但當雙層迭代次數(shù)多時linprog每次重新求解的速度未必慢反而省去安裝商業(yè)求解器的麻煩。麻煩的是迭代求解中價格乘以功率會形成雙線性項這也是博弈問題的核心難點。如果采用KKT單層化這些雙線性項會進入約束需要引入輔助變量和大M法線性化。用大M法時M的值不能取太大否則數(shù)值穩(wěn)定性變差我一般取相關(guān)變量可能上限的10到100倍再根據(jù)結(jié)果微調(diào)。如果你只是想快速驗證算法用迭代法就不需要處理非線性項因為價格固定后下層模型是線性的。雙線性項只出現(xiàn)在上層收益計算里但上層并不需要求解優(yōu)化問題只需要根據(jù)微網(wǎng)結(jié)果計算收益所以雙線性沒關(guān)系。4.2 迭代震蕩與不收斂迭代法最典型的問題是價格來回跳。上層放電價抬高下層減少購買上層又降價下層增加購買循環(huán)往復(fù)無法收斂。我試過幾種方法最有效的是“阻尼更新”new_price old_price alpha * (target_price - old_price);其中alpha是阻尼系數(shù)取值0.2到0.5。另外還可以用最近若干輪的價格平均值做平滑比如用過去三到五次的平均價格作為當前迭代的參考能明顯抑制震蕩。收斂閾值不是越小越好仿真數(shù)據(jù)精度有限設(shè)到1e-3或1e-4就夠用。4.3 單位與量綱一致性這是新手最容易踩的大坑。有些算例里負荷單位是kWh電價單位是元/kWh儲能容量用MWh如果沒換算最后成本數(shù)值會差好幾個數(shù)量級。我建議統(tǒng)一用國際單位制加一個基準值功率MW能量MWh電價元/MWh天然氣價格元/MWh按熱值折算這樣所有目標函數(shù)的單位都是元橫豎能對得上。我在代碼開頭寫了一段單位換算注釋把每個數(shù)據(jù)的原始單位轉(zhuǎn)換成標準單位。4.4 常見報錯速查表我整理了幾個經(jīng)常遇到的報錯以及對應(yīng)的解決方案直接列成表格方便對照。報錯信息原因解決方案Nonlinear constraints require an options object用YALMIP建模時約束中包含非線性項但未指定適合非線性求解器的options用sdpsettings(solver,fmincon)或者將模型線性化No suitable solver found未安裝支持的求解器或YALMIP找不到Gurobi/CPLEX路徑先測試yalmiptest確認求解器路徑配置正確Inf or NaN value in constraints數(shù)據(jù)中出現(xiàn)無窮大或NaN通常是負荷數(shù)據(jù)有缺失用isnan和isinf檢查數(shù)據(jù)Index exceeds array bounds數(shù)組維度不匹配比如負荷數(shù)據(jù)是24×1價格用了24×2統(tǒng)一變量維度利用size檢查Solver output is unbounded模型缺少關(guān)鍵約束比如沒有給燃氣輪機設(shè)上限檢查約束條件是否完整尤其是功率平衡和上下限5. 從復(fù)現(xiàn)到擴展這套代碼還能怎么改如果讀者只是想交仿真作業(yè)跑到第四部分已經(jīng)夠了。但如果要做研究或者工程應(yīng)用這套基礎(chǔ)框架還有很多可以改進的方向。5.1 從確定性到隨機優(yōu)化目前模型假定風電、光伏、負荷都是已知的確定性曲線但實際預(yù)測總有誤差??梢栽谙聦游⒕W(wǎng)模型中引入多個場景變成兩階段隨機優(yōu)化第一階段決定儲能價格第二階段根據(jù)每個場景分別求微網(wǎng)調(diào)度再把場景期望收益作為上層目標。這樣模型更真實但計算量會成倍增加。5.2 多微網(wǎng)與多儲能擴展單個儲能商對應(yīng)單個微網(wǎng)只是最簡單的“一主一從”。實際共享儲能服務(wù)對象是多個微網(wǎng)價格可能一致也可能差異定價。多微網(wǎng)時儲能運營商面對的是多個不同的跟隨者每個跟隨者需求不同博弈的均衡算法變得更加復(fù)雜。不過MATLAB的循環(huán)結(jié)構(gòu)天然適合遍歷多個微網(wǎng)只需要把下層模型封裝成一個獨立函數(shù)循環(huán)多次調(diào)用即可。5.3 碳交易機制接入環(huán)境目標現(xiàn)在是綜合能源微網(wǎng)研究里繞不開的要素。碳交易機制相當于給碳排放量加了一個價格信號微網(wǎng)的燃氣輪機、電網(wǎng)購電對應(yīng)的間接排放都要計入碳成本。這部分擴展只需要在目標函數(shù)里加上一項碳價格乘以碳排放量不需要改變整體博弈框架是性價比很高的擴展方向。5.4 給新手的三個建議最后給想復(fù)現(xiàn)這個項目的朋友三個建議。第一不要一上來就追求完美模型先用一個簡單的單時段或者兩臺設(shè)備的小算例手算或者筆算一遍最優(yōu)調(diào)度結(jié)果再用MATLAB代碼驗證能幫你快速排查建模錯誤。第二先不用急著裝CPLEX或Gurobi先用YALMIP默認求解器和一個小規(guī)模算例跑通整個迭代邏輯再去引入商業(yè)求解器提升效率。第三每次改完代碼把運行結(jié)果和上一次對比如果某個變量突變優(yōu)先檢查單位換算和數(shù)據(jù)拼接而不是懷疑求解器壞了。我自己的體會是主從博弈的MATLAB實現(xiàn)最花時間的其實不是算法本身而是把“經(jīng)濟關(guān)系”轉(zhuǎn)換成“數(shù)學約束”的那一步。充電價、放電價、服務(wù)費、購電成本、收益歸屬這些概念一旦用數(shù)值表達整個模型就立住了。只要你把交易邏輯梳理清楚代碼實現(xiàn)就是水到渠成的事。等你把這一版跑通再回頭去看論文里的KKT條件和強對偶理論會有一種豁然開朗的感覺。