1. 項目概述與論文復現(xiàn)的核心價值先說說這個題目本身。P2GPower to Gas電轉氣是電氣綜合能源系統(tǒng)里這兩年繞不開的一個關鍵環(huán)節(jié)它把富余風電、光伏的電能轉化成氫氣甚至合成天然氣讓電力系統(tǒng)和氣網(wǎng)系統(tǒng)產生雙向耦合。我當初選擇復現(xiàn)這篇論文核心目的是想把這套計及P2G廠站的規(guī)劃模型徹底吃透并且在Matlab環(huán)境里跑通完整的優(yōu)化求解流程——從天然氣網(wǎng)絡建模、電力系統(tǒng)約束、P2G廠站運行特性到規(guī)劃方案的迭代尋優(yōu)每一步都要落地成可執(zhí)行的代碼。為什么值得復現(xiàn)這類論文最大的價值不在于最終的規(guī)劃結果而在于建模思路它把電-氣耦合從單一的電轉氣購能問題升級成了考慮廠站內部運行約束電解槽效率、儲氫罐容量、甲烷化反應熱量平衡的一體化規(guī)劃問題。對于準備做綜合能源系統(tǒng)方向的學生或者剛接觸能源互聯(lián)網(wǎng)優(yōu)化的工程師把這篇論文的代碼復現(xiàn)出來等于一次性打通了電力系統(tǒng)優(yōu)化調度、天然氣網(wǎng)絡穩(wěn)態(tài)分析、混合整數(shù)線性規(guī)劃MILP求解三條技術線。我在復現(xiàn)過程中的體會是這類論文的難點不在數(shù)學推導而在工程落地。學校給的論文材料往往只給最終模型和算例參數(shù)中間過程——比如天然氣管道怎么線性化、網(wǎng)損怎么處理、NLP和MILP怎么切換——需要自己大量補課。下面我把整個復現(xiàn)過程拆開講從建模到代碼實現(xiàn)再到踩坑記錄給你一條相對平整的路。2. 系統(tǒng)建模電網(wǎng)、氣網(wǎng)和P2G廠站的協(xié)同表達2.1 電力系統(tǒng)部分怎么建電力系統(tǒng)在規(guī)劃模型里一般保留節(jié)點功率平衡、機組出力上下限、爬坡約束這幾個核心模塊。復現(xiàn)時我確認了原論文的假設規(guī)劃層不考慮暫態(tài)過程用直流潮流模型近似交流潮流這么做既保留了網(wǎng)絡拓撲約束對規(guī)劃結果的影響又讓模型整體保持線性。需要特別留心的是P2G廠站作為負荷接入時怎么建模。P2G的本質是一個大功率電負荷所以任何含P2G的節(jié)點負荷不再是固定值而是決策變量的一部分。這在實現(xiàn)上要改動現(xiàn)有的節(jié)點功率平衡方程% 節(jié)點功率平衡P_P2G是決策變量 % 原方程: sum(Pg) - sum(Pload) sum(Pij) % 修改后: sum(Pg) - sum(Pload) - P_P2G sum(Pij) balance_eq sum(Pg_idx) - sum(Pload_idx) - sum(P2G_idx) sum(flow_idx);很多同學復現(xiàn)時容易在功率平衡里忘掉P2G這一項直接導致規(guī)劃結果里P2G容量越大系統(tǒng)不平衡量越大仿真結果完全失真。再就是網(wǎng)損。直流潮流模型一般忽略網(wǎng)損但在規(guī)劃問題里全網(wǎng)總有功平衡往往對計算結果影響很大。我看的這篇論文采用的是迭代修正法先求解不計網(wǎng)損的模型再根據(jù)潮流結果計算網(wǎng)損作為已知量回調到下一個迭代中。這種做法實現(xiàn)起來簡單但要注意迭代的收斂條件設置。我在實測中用的是網(wǎng)損前后兩次變化小于0.1%作為停止條件通常在3到5輪就能收斂。2.2 天然氣網(wǎng)絡穩(wěn)態(tài)模型氣網(wǎng)部分是這個領域公認的難點。天然氣管網(wǎng)的核心變量是節(jié)點氣壓和管道流量兩者之間的非線性關系讓問題變得很棘手。常見的處理手段有兩種一是分段線性化piecewise linearization把Weimouth方程按流量區(qū)間分段逼近二是直接在大規(guī)模MILP里用增量線性化方法嵌入。原論文采用的是增量線性化方法這是我復現(xiàn)時印象最深的一塊。Weimouth方程本身是% 管道流量與兩端壓力滿足非線性關系 % F_ij C_ij * sqrt(pi^2 - pj^2) % 令 H_ij pi^2 - pj^2則流量關于H_ij是根號關系 % 線性化: 將H_ij取值范圍分成N段每段內流量線性逼近分段數(shù)量是精度和計算量的折中。我試過5段、8段、10段發(fā)現(xiàn)對規(guī)劃問題來說8段已經足夠分段數(shù)超過10段之后計算時間幾乎翻倍而結果差異不到0.5%。如果論文沒給具體分段數(shù)我建議直接采用8段作為默認參數(shù)這也是該領域文獻里最常見的選擇。氣源、儲氣罐、負荷這三類節(jié)點的約束相對直接主要注意氣源出力的上限和節(jié)點氣壓的上下限。節(jié)點氣壓范圍在規(guī)劃問題里經常被誤設為固定值實際上氣網(wǎng)節(jié)點氣壓允許在合理區(qū)間內浮動比如0.9到1.2倍的基準值。如果氣壓約束過緊會過度限制氣網(wǎng)消納P2G產氣的能力。2.3 P2G廠站內部的黑箱變白箱這篇論文區(qū)別于普通電-氣耦合研究的最大亮點是把P2G廠站內部過程展開了。P2G不是一個簡單的輸入電、輸出氣的黑箱而是一條由電解槽、儲氫罐、甲烷化反應器、氣體壓縮機四部分組成的工藝鏈條。電解槽環(huán)節(jié)的關鍵約束是額定容量和運行范圍。電解槽的輸入功率不能低于某個比例否則電解效率急劇下降所以一般設一個最小運行功率約束比如額定功率的20%。儲氫罐的作用是緩沖電解產氫和甲烷化耗氫之間的時間不匹配它的狀態(tài)方程是一個離散時間遞推式容量約束和初始/終態(tài)儲量約束都必須加進去。甲烷化反應環(huán)節(jié)存在熱量平衡問題——甲烷化是強放熱反應建模時要在P2G功率輸出和運行溫度之間做一個溫度約束的簡化處理。不過大多數(shù)規(guī)劃類論文并不真正求解熱平衡而是直接用氫轉甲烷的轉換效率乘以輸入氫量得到產氣量。如果復現(xiàn)時想更嚴謹可以加一個溫度懲罰項但我個人認為對于年度規(guī)劃問題來說意義不大徒增非線性。氣體壓縮機在P2G廠站模型里往往是最容易被忽略的一塊。P2G產氣壓力通常低于天然氣輸氣管網(wǎng)的壓力等級必須經過壓縮機升壓才能注入氣網(wǎng)。壓縮機本身消耗的功率雖然占比不大約占P2G總耗電的2%到5%但在規(guī)劃模型里如果不計這部分自耗電P2G的凈效率會被高估。我建議至少按壓縮比和流量做一個線性化的功耗估算別完全省略。3. 規(guī)劃模型構建與求解器選型3.1 目標函數(shù)的三層結構原論文的目標函數(shù)是典型的多層規(guī)劃架構投資成本運行成本環(huán)境成本。我復現(xiàn)時把目標函數(shù)拆成了三層方便后續(xù)做敏感性分析投資成本層P2G廠站各設備的單位投資成本乘容量再乘年值系數(shù)注意設備壽命不同折算系數(shù)也不同。電解槽壽命一般按10到15年算甲烷化設備按20年算別統(tǒng)一套一個系數(shù)。運行成本層包括購電成本、購氣成本、機組啟停成本。這里購電成本要區(qū)分分時電價論文算例里通常給的是峰平谷三段電價。環(huán)境成本層按碳排放量折算成懲罰費用。需注意碳價參數(shù)設置原論文一般會說明基準碳價是多少復現(xiàn)時務必核對單位——是元/噸還是元/千克弄錯一個量級整個結果全亂。這里有個心得論文的原始算例數(shù)據(jù)不同成本項權重差異很大。復現(xiàn)前先把目標函數(shù)各成本項的數(shù)量級算一遍如果發(fā)現(xiàn)某一項比另一個項小幾個數(shù)量級大概率是單位問題而非真實差異。3.2 約束條件的層級拆分規(guī)劃模型本質上是雙層問題投資決策長期與運行決策短期。復現(xiàn)時如果直接用一個大規(guī)模MILP求解器硬解整個模型計算規(guī)模會非常恐怖——因為運行層要模擬365天×24小時的調度過程變量總量輕松上萬。原論文這里的處理思路值得學習把規(guī)劃問題拆成主問題和子問題的迭代式。主問題是投資決策輸出P2G廠站建設方案子問題是給定投資方案后的年度運行優(yōu)化輸出運行成本和可行域反饋。主-子問題之間通過Benders分解的思路交互。我在復現(xiàn)中沒有寫完整的Benders分解而是采用了更工程化的啟發(fā)式迭代方案先給一個初始P2G容量猜測值求運行子問題得到該方案下的最優(yōu)運行成本把第一輪結果里被觸發(fā)的容量瓶頸約束提取出來用于修正下一輪的投資方案如此迭代三到四輪便可收斂。必須說明的是這個簡化方案犧牲了嚴格的全局最優(yōu)性。如果審稿要求嚴格的最優(yōu)解正版的Benders或直接MILP求解器是必需的但如果只是做工程方案分析這個迭代法的結果已經足夠可靠而且速度快一個數(shù)量級。3.3 求解器選型與性能表現(xiàn)Matlab環(huán)境下求解MILP問題我比較過幾套方案自帶的intlinprog、 YALMIPGurobi、 YALMIPCPLEX。實際測試結果讓我有點意外intlinprog在中小規(guī)模算例節(jié)點數(shù)少于30表現(xiàn)尚可但一旦進入IEEE 39節(jié)點或118節(jié)點級別的氣電耦合系統(tǒng)intlinprog的求解時間和數(shù)值穩(wěn)定性都肉眼可見地變差。Gurobi在MILP求解上的性能優(yōu)勢非常明顯尤其是大量二元變量的場景。以我復現(xiàn)的30節(jié)點電網(wǎng)加20節(jié)點氣網(wǎng)算例為例Gurobi求解時間約120秒intlinprog則耗了近800秒而且Gurobi的解質量目標值更優(yōu)更好。給一個小建議復現(xiàn)這類論文如果資金寬裕優(yōu)先用Gurobi。如果只有Matlab基礎工具箱也完全可以跑通只是要把算例規(guī)?？刂圃诤侠矸秶鷥炔⒃O置合適的求解精度和最大迭代次數(shù)。4. Matlab代碼實現(xiàn)從框架到核心函數(shù)4.1 代碼整體架構設計我見過不少同學復現(xiàn)代碼時喜歡把所有邏輯寫在一個幾百行的主腳本里面向過程的寫法雖然直白但一旦需要調節(jié)參數(shù)或換算例就變得寸步難行。我這次復現(xiàn)采用了模塊化設計簡單說就是數(shù)據(jù)、模型、求解、結果四層分開項目根目錄/ ├── data/ % 算例數(shù)據(jù)按系統(tǒng)分類存放 ├── models/ % 模型構建函數(shù) ├── solver/ % 求解器調用封裝 ├── results/ % 結果輸出與圖表生成 └── main.m % 主入口整個流程編排這個架構的收益是在換算例時顯現(xiàn)的——只需新增一個data子目錄中的數(shù)據(jù)文件代碼零改動即可運行新的系統(tǒng)。如果你手頭沒有現(xiàn)成的電-氣耦合標準算例用MATPOWER提供的電力系統(tǒng)數(shù)據(jù)搭配一個自建的氣網(wǎng)數(shù)據(jù)文件也能拼湊出可用算例不必非要獲得原論文的配套數(shù)據(jù)。4.2 關鍵數(shù)據(jù)結構設計Matlab編程里最容易被忽視的是數(shù)據(jù)結構設計。我發(fā)現(xiàn)用結構體數(shù)組按對象組織數(shù)據(jù)比散落的命名變量清晰得多。比如電網(wǎng)數(shù)據(jù)可以這樣組織grid.bus struct(id, [], type, [], Pg, [], Pd, [], Vmin, [], Vmax, []); grid.line struct(from, [], to, [], R, [], X, [], capacity, []); grid.gen struct(id, [], bus, [], Pmin, [], Pmax, [], ramp, [], cost, []);P2G廠站的數(shù)據(jù)結構更復雜一些需要包含四個子模塊的參數(shù)p2g.electrolyzer struct(capacity, [], eta_elec, [], Pmin_ratio, [], inv_cost, []); p2g.h2storage struct(capacity, [], init_level, [], final_level, [], inv_cost, []); p2g.methanation struct(capacity, [], eta_meth, [], Q_consume, [], inv_cost, []); p2g.compressor struct(elevation_ratio, [], power_coef, [], inv_cost, []);我踩過的一個坑是初期把所有數(shù)據(jù)分散在不同變量里結果跑大規(guī)模算例時內存管理混亂經常出現(xiàn)變量名寫錯但程序不報錯的情況因為Matlab對變量名檢查不嚴只是邏輯錯了。用結構體之后至少能在視覺層面對數(shù)據(jù)關系一目了然。4.3 模型構建核心流程建模型的過程我分為五步走每一步都通過測試函數(shù)驗證正確性后再進入下一步第一步讀取算例數(shù)據(jù)并預處理。這一步的關鍵是節(jié)點編號的對齊——電網(wǎng)和氣網(wǎng)的節(jié)點編號體系不同必須在數(shù)據(jù)層建立映射關系論文里通常給的算例圖數(shù)據(jù)已經標好直接用即可。第二步構建電網(wǎng)模型。直流潮流中關鍵是導納矩陣B的形成和節(jié)點分類平衡節(jié)點、PV節(jié)點、PQ節(jié)點。給定P2G廠站接入的新節(jié)點編號后在B矩陣中插入對應行列注意新增節(jié)點的基準電壓和基準功率要與全系統(tǒng)一致。第三步構建氣網(wǎng)模型。這部分是代碼里最容易寫錯的。管道流量線性化需要預先計算分段點和斜率我寫成了一個獨立函數(shù)function [seg_points, slopes, intercepts] linearize_weimouth(pipe_C, p_max, p_min, n_seg) % 輸入: 管道常數(shù)C, 壓力上下限, 分段數(shù) % 輸出: 每個分段的端點、斜率、截距 H_max pipe_C^2 * (p_max^2 - p_min^2); H_seg linspace(0, H_max^0.5, n_seg1).^2; F_seg sqrt(H_seg); % 計算各段斜率 slopes diff(F_seg) ./ diff(H_seg); intercepts F_seg(1:end-1) - slopes .* H_seg(1:end-1); seg_points H_seg; end第四步把P2G廠站四個模塊全部納入模型。這一塊我建議畫一個簡易的能流框圖在紙上畫即可不必畫進代碼明確每一級轉換的能量流方向再逐一寫約束。第五步組裝目標函數(shù)和全部約束交給求解器。第五步是整個過程中最耗時的環(huán)節(jié)。我第一次組裝模型時光約束數(shù)量就遇到上百條Matlab命令行窗口里報錯信息滿天飛后來學會了一個技巧每加一組約束后立即求解一個簡化版模型只有該約束相關變量的固定值驗證可行性。這樣做問題定位很快基本不用從頭調試。4.4 求解封裝與結果后處理求解器封裝我寫成了通用接口這樣可以在intlinprog和Gurobi之間無縫切換function [x, fval, exitflag] solve_milp(model, use_gurobi) if use_gurobi % 轉換為Gurobi的輸入格式 result gurobi_optimize(model); x result.x; fval result.objval; exitflag result.status; else options optimoptions(intlinprog, Display, final, ... MaxTime, 1800, RelativeGapTolerance, 0.01); [x, fval, exitflag] intlinprog(model.f, model.intcon, ... model.Aineq, model.bineq, model.Aeq, model.beq, ... model.lb, model.ub, options); end end結果后處理這塊我建議一定做三張核心圖第一張是不同P2G容量方案下的總成本柱狀圖第二張是典型日的電功率和氣功率平衡曲線第三張是P2G廠站內部能量流?；鶊D用Matlab繪圖函數(shù)手動實現(xiàn)。這三張圖是論文復現(xiàn)成果最直觀的呈現(xiàn)方式。特別是不同P2G容量下的成本曲線這張圖幾乎可以一眼看出規(guī)劃方案的最優(yōu)點——總成本最低處對應的P2G容量就是最優(yōu)容量。5. 常見問題與排查技巧實錄5.1 求解不收斂或收斂極慢這是復現(xiàn)此類論文最普遍的問題。我在調試時遇到過不少次模型不收斂的情況原因是多方面的但絕大多數(shù)指向同一個核心——約束條件過于激進。比如P2G年最大利用小時數(shù)設得過高導致投資容量在運行層面無法收回成本模型就會反復嘗試調整投資方案始終找不到可行解。如果遇到收斂慢的問題我強烈建議先檢查以下幾處排查點典型癥狀處理方式P2G最小運行負荷約束模型某時刻P2G出力低于下限改為邏輯約束P2G要么關閉要么至少運行在20%額定功率儲氫罐初末狀態(tài)約束儲氫罐儲量出現(xiàn)不現(xiàn)實的銳減或激增放寬末端狀態(tài)限制如終值在初值±10%之間氣壓節(jié)點上下限個別節(jié)點氣壓越界導致整體不可行適當放寬至1.2倍基準值或檢查是否有氣壓等級設置錯誤分段線性化精度單段斜率過陡導致最優(yōu)解落在分段點附近震蕩增加分段數(shù)或改用均勻殘差誤差分布的分段方式5.2 線性化誤差過大增量線性化方法的誤差主要源于分段點的選取。我最初用等間距分段發(fā)現(xiàn)管道流量較大的情況下誤差可以達到5%以上在系統(tǒng)層面產生可感知的偏差。改進方式是讓分段間距不再均勻而是讓單位區(qū)間內的流量殘差保持一致即誤差等分法。具體實現(xiàn)思路是先初步分段求出各段最大殘差再調整分段點使各段最大殘差相等。這種方法實現(xiàn)較為復雜但在求解效率和精度之間能取得更好平衡。如果只是復現(xiàn)論文結果用均勻分段并取8到10段一般就夠用了。5.3 算例結果與論文不一致的排查這也是復現(xiàn)代碼時非常常見的情況代碼能跑通結果卻和論文對不上。通常問題不在代碼邏輯而在數(shù)據(jù)。建議按以下順序排查第一步檢查單位是否一致。論文中天然氣流量可能是立方米/小時、千克/小時、Mbtu/小時三種不同單位混用換算關系搞錯會直接導致幾十倍的偏差。第二步檢查基準功率和基準電壓設置。ylq9綜合能源系統(tǒng)研究中電力和天然氣的基準值通常設為100MVA和1.0MPa如果基準值取錯整個標幺值體系都會偏移。第三步檢查典型日的選取。原論文的規(guī)劃周期是8760小時但算例里的典型日可能是按季節(jié)或按峰谷時段濃縮出來的。復現(xiàn)時要明確原論文是用了完整的8760小時建模還是用了若干個典型日乘以權重系數(shù)。這兩者的結果會有不小差異。第四步檢查P2G效率的計算方式。是低熱值效率還是高熱值效率兩者數(shù)值相差約10%左右論文如果沒有明確說明復現(xiàn)時很容易出現(xiàn)幾不可見的差異。5.4 求解器數(shù)值穩(wěn)定性問題Matlab的intlinprog在處理具有不同量級數(shù)值的混合整數(shù)問題時容易陷入數(shù)值病態(tài)問題——比如投資成本是千萬量級而運行成本是百萬量級二元變量的目標系數(shù)很小導致MIP啟發(fā)式搜索表現(xiàn)不佳。解決辦法是對模型中的系數(shù)進行歸一化把投資成本除以一個基準值比如總投資上限讓所有目標項的量級集中在1到100之間。這個操作的原理是避免求解器內部處理跨量級數(shù)值時產生舍入誤差。別小看這一步我遇到過一組算例歸一化前后目標值相差4%的情況而這個差異完全來自數(shù)值誤差而非模型變化。6. 復現(xiàn)過程的經驗總結與后續(xù)擴展方向整個項目復現(xiàn)下來我最大的體會是論文復現(xiàn)不是簡單的翻譯代碼而是對建模思路的再發(fā)現(xiàn)。原論文里一筆帶過的許多假設——比如管道流量線性化的具體分段數(shù)、P2G廠站內部壓縮機自耗電的處理方式——恰恰是工程實現(xiàn)中最需要斟酌的細節(jié)。真正跑通一輪完整流程之后你對綜合能源系統(tǒng)規(guī)劃的理解深度會遠遠超過只看理論推導時的水平。如果后續(xù)想在這個基礎上擴展我個人建議幾個方向一是把P2G廠站模型換成更精細的電制氫全鏈條模型加入電解槽的啟停成本和動態(tài)效率曲線二是在規(guī)劃模型中加入不確定性因素比如風電出力和電價的隨機場景三是把天然氣網(wǎng)由穩(wěn)態(tài)模型擴展為動態(tài)模型考慮管道儲氣效應。這三個方向在當前的研究中都很熱門而且都是能在復現(xiàn)代碼基礎上做增量式修改完成的。最后再分享一個經驗復現(xiàn)過程中請務必保存好每一個能運行的版本并做好注釋。你可能現(xiàn)在覺得某個中間版本沒用但三周后當你發(fā)現(xiàn)新模型的求解器行為變得詭異時那個舊版本就是你回溯排查的最佳參照物。每次改動前先跑一遍當前版本記錄目標值和關鍵約束的滿足情況再做修改。這套好習慣可以幫你省掉大量調試時間。