冷熱電微網(wǎng)雙層優(yōu)化這幾年在綜合能源方向幾乎成了標配選題MATLABYALMIP也成了最常見的實現(xiàn)組合。但說實話光看論文摘要覺得很簡單真到自己寫代碼時才發(fā)現(xiàn)一堆坑上層和下層怎么嵌套儲能SOC怎么遞推才不會出現(xiàn)同時充放電YALMIP報Infeasible之后從哪查起這篇內(nèi)容我按自己做模擬項目X時踩過的路把雙層優(yōu)化的建模思路、求解路線、代碼骨架和調參心得一次講清楚。適合正在做微網(wǎng)儲能優(yōu)化、課程設計或者準備論文的同行也適合剛接觸YALMIP、想快速搭起可運行模型的同學。先說清楚一個根本問題雙層優(yōu)化為什么非要用“兩層”而不是一個大模型搞定這直接決定了后面所有代碼怎么寫。1. 先想清楚雙層優(yōu)化到底在優(yōu)化什么1.1 時間尺度差異讓問題天然分層冷熱電微網(wǎng)里面有容量規(guī)劃和運行調度兩類決策它們的決策頻率完全不同。容量規(guī)劃決定“裝多大的燃氣輪機、裝多少儲能、配多少光伏”這個決策一做就是十幾年屬于年尺度甚至項目周期尺度的長期投資問題。運行調度決定“某一天某一時刻燃氣輪機發(fā)多少電、儲能充多少放多少、從電網(wǎng)買多少電”這個決策是小時尺度的短期運行問題。如果把兩類變量塞進一個模型目標函數(shù)里既要算二十年投資回報又要算8760個小時的逐時購電成本求解器面對的量級和語義完全不一樣很難在同一個優(yōu)化框架下做權衡。而且容量和運行本身就存在主從關系定了容量之后系統(tǒng)的運行受容量上限約束反過來運行的經(jīng)濟性又決定了大容量值不值。這種“你先定、我后跑”的主從關系用雙層模型表達最自然。我在做模擬項目X時最大的體會是先別急著寫代碼先把上層變量、下層變量、各自目標、各自約束列成一張表你就能發(fā)現(xiàn)上下層交換的信息到底是什么。上面?zhèn)飨聛淼氖侨萘可舷尴旅鎮(zhèn)骰厝サ氖悄赀\行成本或收益。這不只是概念清楚的問題它直接決定你外層適應度函數(shù)怎么寫。1.2 上層與下層的數(shù)學模型上層優(yōu)化通常是最小化全生命周期總成本或者最大化項目凈收益。一般寫成min C_total C_inv C_om C_run其中C_inv是設備投資的等年值C_om是年維護費用C_run是內(nèi)層優(yōu)化返回的年運行費用購電、購氣成本減去售電收益等。上層決策變量是容量比如燃氣輪機額定功率、儲能額定容量和功率、電制冷機容量、光伏裝機容量。下層優(yōu)化在給定容量下求解典型日或全年時段的運行方案min C_run sum(購電費用 購氣費用 - 售電收益)下層決策變量是各時刻的設備出力、儲能充放電功率、購售電功率、啟停狀態(tài)等。約束包括功率平衡、設備出力上下限、儲能SOC遞推、爬坡約束、備用約束等。注意上層的C_inv需要折算成年值公式是C_inv C_cap_total * r * (1r)^n / ((1r)^n - 1)這里的r是折現(xiàn)率n是項目年限。很多新手直接把總投資加進目標函數(shù)內(nèi)層返回的卻是一年費用這樣量綱不一致優(yōu)化出來的容量必然偏小。我一開始也在這個地方迷糊過后來統(tǒng)一把上層所有成本都折算到年才得到合理結果。1.3 儲能是鏈接上下層的關鍵樞紐儲能裝置在這類模型中非常特殊它的容量屬于上層變量它的充放電功率屬于下層變量而且它的運行狀態(tài)不能孤立確定必須滿足跨時段遞推關系。換句話說儲能是上下層信息交換的核心載體。上層決定了電池裝多大下層只能在“這塊電池”的范圍內(nèi)動態(tài)調度。實話說如果模型里沒有儲能雙層優(yōu)化就會退化成單層配置問題加一堆運行約束復雜度降一個量級。正因為儲能引入了SOC的狀態(tài)變量和時序耦合下層才真正變成一個動態(tài)優(yōu)化問題上層也才需要通過反復迭代來尋找容量最優(yōu)值。這也是為什么絕大多數(shù)論文和項目都把儲能作為優(yōu)化重點的原因。2. 微網(wǎng)模型搭建把冷熱電和儲能寫成數(shù)學表達式2.1 能量流與設備庫先搭能量流框架。冷熱電微網(wǎng)系統(tǒng)通常包含三條母線電母線、熱母線和冷母線。設備按“源—轉換—存儲—負荷”四類劃分。常見組成源側光伏、風電、電網(wǎng)購電轉換設備燃氣輪機發(fā)電兼余熱回收、燃氣鍋爐補熱、吸收式制冷機用熱制冷、電制冷機用電制冷存儲設備蓄電池存電、蓄熱罐存熱、蓄冷罐存冷負荷電負荷、熱負荷、冷負荷這三條母線之間通過設備耦合比如燃氣輪機發(fā)電的同時產(chǎn)生余熱余熱進入熱母線一部分直接供熱一部分進吸收式制冷機轉成冷。你可以想象成電力、熱水、冷水三套管路在一個園區(qū)里交織優(yōu)化要做的就是決定每一時刻能量往哪條管路走、走多少。2.2 關鍵設備數(shù)學模型設備建模不用搞得太復雜工程上常用簡化線性模型保主要特性即可。燃氣輪機或內(nèi)燃機的發(fā)電模型P_gt eta_gt * LHV * F_gas其中eta_gt是發(fā)電效率LHV是天然氣低位熱值F_gas是燃料消耗速率。它的余熱回收量可以寫成Q_hr alpha_hr * P_gtalpha_hr是熱電比系數(shù)不同機型的溫差、排煙熱量不同取值通常在0.5到1.5之間。我曾經(jīng)見過有人把Q_hr寫成(1-eta_gt)*P_gt這其實是錯的因為排煙損失、散熱損失都沒算進去回收不到那么多熱量。燃氣鍋爐模型Q_boiler eta_boiler * LHV * F_boiler吸收式制冷機Q_cold_abs COP_abs * Q_heat_in電制冷機Q_cold_ec COP_ec * P_ec光伏模型可以簡化成P_pv(t) P_pv_cap * irradiance_norm(t)irradiance_norm是歸一化輻照序列如果有溫度系數(shù)可以再加一個溫度修正項但對容量規(guī)劃問題來說歸一化輻照已經(jīng)夠用。把這些設備約束寫成線性不等式組交給YALMIP處理后交給求解器整個過程會很順暢。關鍵是別在建模階段就引入非線性比如燃料成本函數(shù)做成分段線性、設備效率做成常數(shù)線性模型跑出來的結果已經(jīng)足夠支撐容量決策趨勢分析。2.3 儲能建模是雙層耦合的核心蓄電池SOC遞推公式是這類模型最容易出錯的地方SOC(t1) SOC(t) eta_ch * P_ch(t) * dt / E_cap - P_dis(t) * dt / (eta_dis * E_cap)注意充電效率eta_ch是乘在充電功率上的放電效率eta_dis是放在分母上的。含義很簡單充進去1度電電池實際存下的只有eta_ch度放出來1度電電池內(nèi)部需要消耗1/eta_dis度的能量。效率位置一旦寫反SOC遞推就會出現(xiàn)能量不守恒運行成本也會虛高。重要約束0 P_ch(t) P_ch_max 0 P_dis(t) P_dis_max SOC_min SOC(t) SOC_max還要處理一個工程細節(jié)防止同一時刻既充電又放電。YALMIP里可以直接加互斥約束也可以引入二進制變量P_ch(t) M * u_ch(t) P_dis(t) M * (1 - u_ch(t))其中M是足夠大的常數(shù)u_ch是0-1變量。有時候電池模型本身就允許輕微同時充放但為了競賽和論文評審的嚴謹性建議還是加上互斥約束。蓄熱罐、蓄冷罐的模型與電池類似只是把電功率換成熱功率/冷功率效率和容量約束同理。儲能SOC的初值和終值也要寫成等式或不等式比如要求一個調度周期結束時SOC回到初始值這樣才能讓日調度方案可重復執(zhí)行。3. 求解路線選擇主流的三條打法3.1 三種路線對比雙層優(yōu)化的求解方法有不少我實際用下來主要分三條路線求解路線基本思路優(yōu)點缺點適用范圍智能算法內(nèi)層數(shù)學規(guī)劃外層用遺傳算法/粒子群搜容量內(nèi)層用YALMIP求解運行調度實現(xiàn)靈活、通用性強、容易調試計算量大結果有一定隨機性絕大多數(shù)工程和論文場景KKT條件轉單層把下層問題用KKT條件等價替換合并成單層約束數(shù)學嚴格解是確定的推導繁瑣要求下層問題凸且強對偶成立下層為線性規(guī)劃的簡化場景MPEC/強對偶變換將對偶變量納入模型形成含均衡約束的數(shù)學規(guī)劃理論完整學術價值高建模復雜求解器支持有限小規(guī)模理論證明類題目我最推薦也最常用的是第一條路線——遺傳算法YALMIP。理由很現(xiàn)實冷熱電微網(wǎng)下層通常是混合整數(shù)線性規(guī)劃KKT條件寫起來極其痛苦MPEC對求解器和建模技巧要求又高而智能算法套YALMIP幾乎不需要額外推導只要把上下層接口寫好就能跑。計算量確實大但通過典型日加權、并行計算完全能接受。3.2 遺傳算法作為外層框架的核心流程外層GA的框架大概是這樣的初始化種群每個個體是一組容量配置比如[P_gt_cap, E_sto_cap, P_sto_ch_max, P_ec_cap, P_pv_cap]對每個個體調用內(nèi)層運行優(yōu)化得到典型日運行成本用運行成本加上投資年值計算適應度選擇、交叉、變異生成下一代迭代直到收斂或達到最大代數(shù)這里有個很關鍵的設計點GA的適應度函數(shù)中運行成本必須是從內(nèi)層優(yōu)化器返回的不能自己憑經(jīng)驗估。因為內(nèi)層調度結果直接反映容量配置是否合理容量太小會導致無法滿足負荷或大量購電容量太大會造成投資浪費。GA搜索的過程本質上就是不斷讓內(nèi)層調度去驗證各種容量組合的運行效果。我在編碼時用實數(shù)編碼而不是二進制編碼原因是容量變量連續(xù)且邊界明確實數(shù)編碼直接對應實際工程單位避免編解碼誤差也方便后期加約束。搜索范圍設置成設備可能容量的1.2倍左右不要設太大否則收斂很慢。3.3 為什么內(nèi)層必須用YALMIP而不是手寫求解內(nèi)層運行調度是多時段、多設備、多約束的優(yōu)化問題手寫求解器或者用單純形法自己編碼實現(xiàn)工作量巨大且極易出錯。YALMIP的核心價值在于它讓你用聲明式的方式描述優(yōu)化問題而不需要關心標準型轉換、求解器接口這些底層細節(jié)。你只需要定義sdpvar決策變量用Constraints [Constraints, ...]累加約束再調用optimize剩下的交給求解器。還有一個好處YALMIP天然支持混合整數(shù)變量。當你需要處理機組啟停、儲能互斥、階梯電價檔位選擇時只需要把一個變量聲明為binvar求解器會自動按混合整數(shù)規(guī)劃處理。如果手寫求解器光是分支定界就得寫一堆代碼。4. MATLABYALMIP代碼骨架與實現(xiàn)要點4.1 工程目錄結構我習慣把項目拆成清晰的功能模塊模擬項目X的目錄如下/project_x/ data/load_curve.mat % 負荷、輻照、電價等時間序列 scripts/main.m % 主程序入口 scripts/upper_ga.m % 外層遺傳算法 scripts/inner_schedule.m % 內(nèi)層運行調度 utils/calc_investment.m % 投資年值計算 output/results.mat % 結果輸出主程序只需要讀數(shù)據(jù)、設置優(yōu)化參數(shù)、調用上層GA、輸出結果。內(nèi)層調度函數(shù)接收容量參數(shù)和基礎數(shù)據(jù)返回運行成本。這樣寫的好處是可以單獨調試內(nèi)層確認無誤后再接外層避免兩邊同時出錯時無從排查。4.2 內(nèi)層運行調度代碼骨架內(nèi)層函數(shù)核心是YALMIP建模。以一個24小時典型日為例function run_cost inner_schedule(cap, data) T 24; % 決策變量 P_gt sdpvar(T, 1); % 燃氣輪機發(fā)電功率 P_ch sdpvar(T, 1); % 儲能充電功率 P_dis sdpvar(T, 1); % 儲能放電功率 SOC sdpvar(T 1, 1); % 儲能SOC P_grid sdpvar(T, 1); % 購電功率 P_ec sdpvar(T, 1); % 電制冷機功率 % 基礎數(shù)據(jù) P_load_e data.P_load_e; % 電負荷 Q_load_h data.Q_load_h; % 熱負荷 Q_load_c data.Q_load_c; % 冷負荷 P_pv data.P_pv; % 光伏出力序列 price data.price; % 分時電價 Constraints []; for t 1:T % 設備出力上下限 Constraints [Constraints, 0 P_gt(t) cap.P_gt_cap]; Constraints [Constraints, 0 P_ch(t) cap.P_ch_max]; Constraints [Constraints, 0 P_dis(t) cap.P_dis_max]; Constraints [Constraints, 0 P_ec(t) cap.P_ec_cap]; % 儲能SOC遞推 Constraints [Constraints, ... SOC(t1) SOC(t) 0.95 * P_ch(t) * 1 / cap.E_sto ... - P_dis(t) * 1 / (0.95 * cap.E_sto)]; Constraints [Constraints, 0.1 SOC(t1) 0.9]; % 電功率平衡 Constraints [Constraints, ... P_gt(t) P_dis(t) P_pv(t) P_grid(t) ... P_load_e(t) P_ch(t) P_ec(t)]; % 冷功率平衡吸收式制冷用熱、電制冷用電 % 這里簡化吸收式制冷量 alpha_cool * Q_hr(t) % 詳細模型視具體設備配置而定 end % 儲能SOC首尾約束 Constraints [Constraints, SOC(1) SOC(T 1)]; % 目標函數(shù)購電費用 - 售電收益 購氣費用 gas_cost sum(P_gt) * data.gas_price / data.eta_gt; elec_cost sum(P_grid .* price); Objective elec_cost gas_cost; % 求解 ops sdpsettings(solver, gurobi, verbose, 0); optimize(Constraints, Objective, ops); % 返回目標值 run_cost value(Objective); end注意我在SOC遞推里用的充電效率0.95放電效率0.95單位是小時所以dt1。實際項目中如果時間步長是15分鐘記得把dt帶入公式否則SOC變化速度會明顯失真。還有一個容易被忽略的點電功率平衡中的P_grid(t)如果允許負數(shù)代表向電網(wǎng)售電。這時目標函數(shù)里要區(qū)分購電價格和售電價格否則求解器會把“先多買再賣”當成套利手段模糊化處理會讓結果失真。4.3 外層遺傳算法主循環(huán)與并行加速外層GA的核心循環(huán)大致如下function [best_sol, best_fit] upper_ga(data, ga_options) pop_size ga_options.pop_size; max_gen ga_options.max_gen; bounds ga_options.bounds; % 初始化種群均勻隨機分布在邊界內(nèi) pop zeros(pop_size, n_vars); for i 1:pop_size for j 1:n_vars pop(i, j) bounds(j, 1) rand * (bounds(j, 2) - bounds(j, 1)); end end for gen 1:max_gen fitness zeros(pop_size, 1); parfor i 1:pop_size cap decode(pop(i, :)); C_inv calc_investment(cap); % 投資年值 C_run inner_schedule(cap, data); % 內(nèi)層運行成本 fitness(i) C_inv C_run; % 總適應度 end % 選擇、交叉、變異 % ... end end這里有幾個實戰(zhàn)經(jīng)驗parfor并行很有用但要確保每個worker都能訪問YALMIP和求解器。如果并行池初始化后報找不到求解器就在并行循環(huán)內(nèi)加上一行yalmip(clear)或重新sdpsettings。先串行跑一遍小規(guī)模用例確認內(nèi)層全部沒問題再開并行。并行調試一旦報錯錯誤提示往往非常隱晦。GA種群一般設30到60迭代次數(shù)50到100。太小收斂不好太大計算時間成倍增加性價比不高。4.4 求解器配置要點YALMIP本身不是求解器它只是建模層。如果沒有安裝求解器optimize會報錯。常用的組合是商業(yè)求解器如Gurobi或CPLEX它們對線性規(guī)劃和混合整數(shù)線性規(guī)劃的求解效率非常高。個人或者學術用途可以申請免費許可證。配置方式ops sdpsettings(solver, gurobi, verbose, 1);verbose設為1時會輸出求解過程中的迭代信息方便查看問題規(guī)模和解的狀態(tài)。跑成功后檢查value(Objective)是否合理同時確認YALMIP返回的problem字段為00代表求解成功1代表不可行其他值代表數(shù)值問題。如果實在裝不上商業(yè)求解器MATLAB自帶的linprog和intlinprog也能救急。但要注意純純的linprog不能處理0-1變量遇到啟停約束就得用intlinprog或另想辦法。我的建議是盡早把商業(yè)求解器配好后面能省掉無數(shù)折騰時間。5. 實測復盤模擬項目X的完整參數(shù)與結果5.1 項目參數(shù)與數(shù)據(jù)設置模擬項目X是一個園區(qū)級冷熱電微網(wǎng)基本參數(shù)如下參數(shù)數(shù)值說明光伏裝機范圍50 ~ 300 kW上層待優(yōu)化燃氣輪機容量范圍0 ~ 500 kW上層待優(yōu)化儲能容量范圍50 ~ 500 kWh上層待優(yōu)化儲能功率容量/4按2小時充放倍率設置折現(xiàn)率6%投資年值計算項目年限10年等年值折算負荷數(shù)據(jù)用三個典型日加權冬季典型日權重90天、夏季典型日權重90天、過渡季典型日權重185天。每個典型日24個點既有冷熱電負荷曲線也有歸一化光伏輻照曲線。這樣做比直接跑8760小時快得多同時仍能反映季節(jié)差異。分時電價采用峰谷平三段式峰時9點到12點、17點到20點電價高平時6點到9點、12點到17點、20點到22點谷時22點到次日6點。這個設置很常見也是儲能套利空間的主要來源。5.2 雙層迭代的收斂情況我試了兩組參數(shù)第一組GA種群30代、迭代40代第二組種群60代、迭代80代。第二組最后得到的容量配置比第一組總成本低3%左右但計算時間從40分鐘漲到近3小時。這說明GA收斂性對參數(shù)敏感建議先小規(guī)模跑通再按需增加精度。內(nèi)層單次調用耗時約0.3秒到1秒主要取決于YALMIP的約束個數(shù)和求解器的效率。上層一次完整迭代需要調用種群規(guī)模次內(nèi)層這也是總計算時間的大頭。收斂曲線通常表現(xiàn)為前10代快速下降后面進入緩慢優(yōu)化階段。如果前10代沒有明顯下降檢查一下適應度函數(shù)有沒有問題比如投資年值折算錯誤、內(nèi)層返回值是負數(shù)等。5.3 儲能容量對經(jīng)濟性的影響我單獨固定其他容量掃描儲能容量從50 kWh到500 kWh得到年總成本的變化趨勢。結果并不是儲能越大越好容量從50增加到200 kWh時年運行成本下降明顯購電費用減少儲能套利收益增加但超過200 kWh后投資年值的增加開始超過運行費用的節(jié)省總成本反而回升。這背后邏輯很好理解儲能越大能轉移的電量越多但邊際套利收益是遞減的——每天就那么多峰谷價差電池一天只能充放一次兩次裝得再大也消化不了而投資成本是線性的所以在某個拐點之后儲能就成了“賠錢貨”。拐點位置取決于峰谷價差、電池單位投資成本和循環(huán)壽命。這就解釋了為什么儲能容量必須靠優(yōu)化而不是拍腦袋定。輸出結果中我還對比了是否考慮SOC首尾約束的差異。不加首尾約束時模型會把有效容量放大實際運行時會出現(xiàn)第一天從滿電開始、最后一天放空的“作弊”現(xiàn)象年成本偏低結論失真。加上SOC回位約束之后結果才真正可落地。6. 常見問題與排查心得6.1 問題速查表實際跑模型時遇到的問題我整理了一張速查表癥狀可能原因排查思路YALMIP報Infeasible約束沖突、變量無界、系數(shù)單位錯誤逐條注釋約束縮小變量范圍定位求解成功但結果明顯不合理SOC遞推效率寫反、電價單位不一致檢查公式推導和單位換算運行成本為負售電收益率遠高于購電成本模型在套利區(qū)分購電/售電價格加數(shù)值上限GA長時間不收斂種群太小、搜索范圍過大、適應度噪聲大增大種群、縮小邊界、固定隨機種子parfor報求解器缺失并行池worker沒有加載YALMIP和求解器在并行循環(huán)內(nèi)重新配置求解器路徑優(yōu)化速度極慢約束太多、變量規(guī)模大、求解器配置不合理用典型日代替8760h簡化設備模型6.2 不可行問題的排查流程遇到Infeasible千萬別慌它不是“代碼壞了”而是數(shù)學上約束自相矛盾。我的排查順序是先求解一個“只有運行調度、不包含投資決策”的內(nèi)層問題如果這個都不可行說明負荷平衡或SOC約束本身有問題。把SOC的上下界從0.1~0.9放寬到0~1排除SOC邊界太緊導致的無解。把購電功率上限從1000 kW放寬到1e6 kW排除購電上限限制導致的不可行。檢查冷熱電平衡公式尤其是單位kW和kWh不能混熱量單位要用一致。最后用YALMIP的diagnostic函數(shù)或者把Constraints拆成兩部分分別求解很快就能定位沖突的那個約束組。還有一個容易被忽略的原因儲能SOC初值和終值相等加上充放電互斥可能在高負荷日沒有可行解。這時可以把首尾約束從等式改成不等式比如允許終值不低于初值的90%既保住可重復性又給求解器更多自由度。6.3 計算時間優(yōu)化心得雙層優(yōu)化最煩人的就是慢。8760小時全時段內(nèi)層優(yōu)化會讓YALMIP建模變得很重變量數(shù)量動輒十幾萬以上求解時間讓人崩潰。我的經(jīng)驗是按季節(jié)選典型日每個典型日24個點乘上對應的天數(shù)權重。這樣內(nèi)層問題規(guī)模小兩個量級結果卻能反映全年的經(jīng)濟性走向。如果一定要跑全年可以考慮把儲能SOC的跨天耦合拆開或者用滾動窗口近似。不過對容量規(guī)劃來說典型日法完全夠用沒必要一開始就追求全年8760小時。還有一個提速技巧內(nèi)層問題盡量保持線性避免分段二次或非線性目標。YALMIP遇到非線性會嘗試調用非線性求解器速度慢且容易不穩(wěn)定。把設備效率設為常數(shù)、燃料成本線性化、階梯電價用0-1變量表達雖然約束多了但求解起來反而快。最后再說一個小技巧。寫外層GA時給每個容量變量設一個合理的離散精度比如儲能容量按10 kWh為步長取整光伏按5 kW取整。這樣既壓縮了搜索空間也讓結果更貼近實際工程選型不至于優(yōu)化出一個“203.7 kWh”的詭異容量。實話說這個項目反復調試下來我最深的感受是雙層優(yōu)化的難點不在算法本身而在于把物理問題嚴謹?shù)剞D成數(shù)學問題再嚴謹?shù)剞D成代碼。先把單層運行模型調通再套外層迭代一步步來比任何捷徑都穩(wěn)。