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