合調(diào)度微電網(wǎng)運(yùn)行優(yōu)化模型：MATLAB+YALMIP實(shí)現(xiàn)解析)
做微電網(wǎng)運(yùn)行優(yōu)化繞不開(kāi)一個(gè)樸素的問(wèn)題電負(fù)荷和熱負(fù)荷明明共存于同一個(gè)園區(qū)為什么大多數(shù)優(yōu)化模型只盯著“電”在算我最近在某園區(qū)微電網(wǎng)項(xiàng)目的經(jīng)濟(jì)調(diào)度分析里把基于電熱聯(lián)合調(diào)度的微電網(wǎng)運(yùn)行優(yōu)化模型在 MATLAB 里完整跑通了一遍。今天把這個(gè)模型的變量設(shè)計(jì)、約束寫(xiě)法、求解代碼和調(diào)試過(guò)程一起整理出來(lái)適合正在做微電網(wǎng)調(diào)度、綜合能源系統(tǒng)優(yōu)化的研究生和工程師參考尤其是手里有燃?xì)廨啓C(jī)、熱電聯(lián)產(chǎn)機(jī)組、電鍋爐、蓄電和蓄熱設(shè)備想算清楚“一套設(shè)備怎樣組合運(yùn)行一天最省錢(qián)”的這類場(chǎng)景。先說(shuō)明白這個(gè)模型能回答什么問(wèn)題光伏和風(fēng)電出力有波動(dòng)電價(jià)有峰谷天然氣也有價(jià)格園區(qū)既有電負(fù)荷又有熱負(fù)荷電鍋爐可以把電轉(zhuǎn)成熱熱電聯(lián)產(chǎn)機(jī)組CHP可以同時(shí)產(chǎn)電產(chǎn)熱蓄電罐和蓄熱罐可以把能量跨越時(shí)間搬運(yùn)。電熱聯(lián)合調(diào)度就是在滿足電、熱兩類負(fù)荷的前提下決定每一臺(tái)設(shè)備每個(gè)時(shí)段開(kāi)多大功率、充多少放多少讓全天總運(yùn)行成本最低。這本質(zhì)是一個(gè)混合整數(shù)線性規(guī)劃問(wèn)題MATLAB 加 YALMIP 工具箱處理起來(lái)非常順手。1. 電熱聯(lián)合調(diào)度的優(yōu)化空間為什么單算“電”會(huì)吃虧1.1 微電網(wǎng)里的“熱”角色從熱電聯(lián)產(chǎn)到電鍋爐大多數(shù)微電網(wǎng)項(xiàng)目在設(shè)計(jì)初期最容易犯的錯(cuò)是把熱負(fù)荷單獨(dú)扔給燃?xì)忮仩t去扛電負(fù)荷單獨(dú)靠光伏、儲(chǔ)能和電網(wǎng)。這樣做邏輯清晰但經(jīng)濟(jì)性很差。因?yàn)閳@區(qū)里大量熱負(fù)荷存在的地方往往同時(shí)存在電負(fù)荷而 CHP 機(jī)組本身就是“一份燃?xì)馔度胪瑫r(shí)產(chǎn)出電和熱”的設(shè)備。如果你只讓它發(fā)電、不讓它供熱相當(dāng)于花了同樣的燃料錢(qián)白白扔掉了一半以上的可用能量。反過(guò)來(lái)電鍋爐也是一種耦合設(shè)備它消耗電能產(chǎn)出熱能。電鍋爐本身不一定省錢(qián)但在低谷電價(jià)時(shí)段或者在光伏大發(fā)導(dǎo)致棄光時(shí)用電制熱可能比用天然氣制熱便宜得多。于是“電”和“熱”會(huì)通過(guò) CHP 和電鍋爐形成兩條耦合通道。聯(lián)合調(diào)度模型的價(jià)值就是把這兩條通道放進(jìn)同一個(gè)優(yōu)化問(wèn)題里讓系統(tǒng)自己決定何時(shí)走“氣到電到熱”的路徑何時(shí)走“電網(wǎng)電到熱”的路徑。CAT在一個(gè)典型的園區(qū)場(chǎng)景里加入蓄熱罐后熱負(fù)荷甚至可以像電儲(chǔ)能一樣實(shí)現(xiàn)時(shí)移白天光伏多的時(shí)候多產(chǎn)熱存起來(lái)晚上熱負(fù)荷上來(lái)再放熱。這個(gè)過(guò)程單看電力系統(tǒng)模型是完全看不見(jiàn)的但它在實(shí)際運(yùn)行中能顯著降低頂峰時(shí)段購(gòu)電壓力。1.2 一個(gè)直觀案例分產(chǎn)調(diào)度和聯(lián)合調(diào)度的成本差在哪我在調(diào)試模型時(shí)用了兩個(gè)對(duì)照方案。第一個(gè)方案是分產(chǎn)調(diào)度CHP 按“以熱定電”的方式運(yùn)行也就是熱負(fù)荷多少CHP 就發(fā)對(duì)應(yīng)比例的電缺電直接從電網(wǎng)買(mǎi)電鍋爐和蓄熱罐不參與聯(lián)動(dòng)。第二個(gè)方案是電熱聯(lián)合調(diào)度所有設(shè)備統(tǒng)一優(yōu)化。同樣一組典型日負(fù)荷數(shù)據(jù)、同樣一套設(shè)備參數(shù)分產(chǎn)調(diào)度一天運(yùn)行成本大約是 14230 元聯(lián)合調(diào)度是 12350 元節(jié)省了約 13.2%。這 13.2% 主要來(lái)自三個(gè)地方。一是低谷電價(jià)時(shí)段電鍋爐啟動(dòng)儲(chǔ)熱替代了高峰時(shí)段部分燃?xì)庀?。二?CHP 不再被“以熱定電”綁死而是更多在峰電時(shí)段滿發(fā)在低谷時(shí)段壓低出力或者停機(jī)熱缺口由蓄熱罐釋放的熱量補(bǔ)上。三是光伏大發(fā)時(shí)段系統(tǒng)主動(dòng)減少購(gòu)電甚至向電網(wǎng)賣(mài)電把多余電能轉(zhuǎn)化為蓄熱罐里的熱能避免棄光。這三個(gè)收益來(lái)源離開(kāi)電熱聯(lián)合調(diào)度模型都很難被量化。2. 模型的數(shù)學(xué)骨架目標(biāo)函數(shù)、約束與變量設(shè)計(jì)2.1 決策變量的分層設(shè)計(jì)連續(xù)功率變量與啟停0-1變量?jī)?yōu)化模型的第一步是確定決策變量我的做法是分成連續(xù)變量和二進(jìn)制變量?jī)深?。連續(xù)變量包括各設(shè)備每個(gè)時(shí)段的功率和能量狀態(tài)二進(jìn)制變量主要用于 CHP 的啟停狀態(tài)。設(shè)備連續(xù)變量說(shuō)明CHP發(fā)電功率 P_chp、供熱功率 H_chp單位 kWCHP 啟停0-1 變量 on_chp1 表示運(yùn)行電鍋爐耗電功率 P_eb、產(chǎn)熱功率 H_ebH_eb eta_eb * P_eb蓄電充電功率 P_ch、放電功率 P_dis、電量 E_bess電量單位 kWh蓄熱蓄熱功率 H_tes_ch、放熱功率 H_tes_dis、熱量 E_tes熱量單位 kWh聯(lián)絡(luò)線購(gòu)電 P_buy、售電 P_sell單位 kWCHP 的啟停變量不是可選項(xiàng)而是必須項(xiàng)。因?yàn)槿細(xì)廨啓C(jī)存在最小技術(shù)出力比如 50kW如果只建模成連續(xù)變量?jī)?yōu)化器為了省錢(qián)會(huì)把 CHP 壓到 1kW、2kW 這種根本不存在的運(yùn)行點(diǎn)。引入 on_chp 后功率下限寫(xiě)成P_chp_min * on_chp P_chp P_chp_max * on_chp機(jī)組要么不開(kāi)要么至少開(kāi)到 50kW這才符合實(shí)際。一個(gè)需要提前想清楚的點(diǎn)蓄電和蓄熱設(shè)備的“充放互斥”要不要加二進(jìn)制變量。從理論嚴(yán)謹(jǐn)性角度應(yīng)該避免同時(shí)充電和放電。但在實(shí)際模型中由于充放都有損耗且損耗已經(jīng)進(jìn)了目標(biāo)函數(shù)優(yōu)化器在經(jīng)濟(jì)上沒(méi)有動(dòng)機(jī)讓同一臺(tái)儲(chǔ)能同時(shí)充放。所以我通常不加這組二進(jìn)制變量后續(xù)求解速度會(huì)快很多。2.2 電熱平衡約束與CHP可行域約束的建模細(xì)節(jié)電功率平衡是微電網(wǎng)模型的地基任何一個(gè)時(shí)段都必須滿足pv(t) wt(t) P_chp(t) P_dis(t) P_buy(t) P_load(t) P_eb(t) P_ch(t) P_sell(t)左邊是電源側(cè)右邊是負(fù)荷側(cè)。注意電鍋爐的耗電功率 P_eb 必須放在負(fù)荷側(cè)這是初學(xué)者最容易漏掉的一項(xiàng)。光伏、風(fēng)電數(shù)據(jù)如果來(lái)自預(yù)測(cè)曲線可以直接作為已知量放進(jìn)等式。如果考慮棄光就需要給光伏增加一個(gè)“可削減”變量。熱功率平衡相對(duì)簡(jiǎn)單H_chp(t) H_eb(t) H_tes_dis(t) H_load(t) H_tes_ch(t)蓄熱罐充熱時(shí)算是熱負(fù)荷放熱時(shí)算是熱源。CHP 的電氣特性是電熱聯(lián)合調(diào)度里最關(guān)鍵的約束。我用固定熱電比模型也就是H_chp(t) rho_chp * P_chp(t)這種模型用起來(lái)最省事適合做可行性分析。它意味著 CHP 發(fā)多少電就附帶產(chǎn)生多少熱。當(dāng)熱負(fù)荷很小、電負(fù)荷很大時(shí)多余的熱可以進(jìn)蓄熱罐當(dāng)熱負(fù)荷很大、電負(fù)荷很小時(shí)則只能通過(guò)降低 CHP 出力來(lái)減少熱產(chǎn)出缺電部分由電網(wǎng)補(bǔ)。也可以把熱電比關(guān)系改成可行域約束允許熱電比在一定范圍內(nèi)變化但那樣會(huì)引入更復(fù)雜的耦合不等式對(duì)線性求解器不友好。CHP 的天然氣消耗我采用的是線性近似gas_consume(t) a_chp * P_chp(t) b_chp * on_chp(t)其中 a_chp 是邊際氣耗系數(shù)單位是 m3/kWhb_chp 是空載氣耗系數(shù)單位是 m3/h。這兩個(gè)系數(shù)由 CHP 性能曲線擬合得到。這種方法比直接用效率除法更貼近實(shí)際也保留了 MILP 的線性結(jié)構(gòu)。2.3 儲(chǔ)能SOC、爬坡、購(gòu)售電約束怎么列才不踩坑儲(chǔ)能約束最容易踩坑的是 SOC 的索引長(zhǎng)度。如果優(yōu)化周期是 T24充放電功率變量長(zhǎng)度是 24那么 SOC 狀態(tài)變量長(zhǎng)度應(yīng)該是 25因?yàn)橐獜某跏紩r(shí)刻一直遞推到 T 時(shí)刻結(jié)束。我習(xí)慣用 E_bess(1) 表示初始電量E_bess(T1) 表示最終電量中間遞推關(guān)系寫(xiě)成E_bess(t1) E_bess(t) dt * (eta_ch * P_ch(t) - P_dis(t) / eta_dis)充電時(shí)乘以效率放電時(shí)除以效率。這個(gè)細(xì)節(jié)特別重要如果方向?qū)懛磧?yōu)化器會(huì)發(fā)現(xiàn)“充進(jìn)去 1 度電放出來(lái) 1.2 度電”相當(dāng)于憑空創(chuàng)造能量結(jié)果會(huì)出現(xiàn)一堆反直覺(jué)的充放行為。容量約束要加在電量狀態(tài)變量上而不是加在充放功率上E_bess_min E_bess(t) E_bess_max同時(shí)還要給充放電功率分別設(shè)置上限。初始和終值電量可以設(shè)定成同一個(gè)值保證日內(nèi)調(diào)度的可持續(xù)性。如果允許終值高于初始值調(diào)度結(jié)果其實(shí)就是“一天之內(nèi)把儲(chǔ)能充滿”成本會(huì)被人為壓低沒(méi)有實(shí)際參考意義。CHP 的爬坡約束是另一個(gè)容易出問(wèn)題的地方。爬坡定義是相鄰時(shí)段功率變化量-ramp_chp * dt P_chp(t1) - P_chp(t) ramp_chp * dt注意在 MATLAB 中的索引是從 1 開(kāi)始所以代碼里會(huì)有P_chp(2:T) - P_chp(1:T-1)的寫(xiě)法。購(gòu)售電約束方面購(gòu)電和售電功率各設(shè)上限要不要加“不同時(shí)購(gòu)售”的互斥二進(jìn)制變量看實(shí)際需求。我一般先不加因?yàn)榉謺r(shí)電價(jià)下基本不會(huì)出現(xiàn)同時(shí)購(gòu)售的荒謬解。如果最終結(jié)果異常再補(bǔ)一組互斥變量也不遲。3. MATLAB代碼實(shí)現(xiàn)數(shù)據(jù)組織、建模與求解器對(duì)接3.1 參數(shù)與負(fù)荷曲線如何整理成矩陣MATLAB 做優(yōu)化建模第一步就是統(tǒng)一數(shù)據(jù)格式。我強(qiáng)烈建議所有時(shí)序數(shù)據(jù)都用1 x T的行向量組織電價(jià)、負(fù)荷、光伏、風(fēng)電都是如此。這樣后續(xù)與 YALMIP 的變量做點(diǎn)乘運(yùn)算時(shí)可以直接.*不需要來(lái)回轉(zhuǎn)置。單位也要統(tǒng)一。功率單位用 kW能量單位用 kWh電價(jià)單位用 元/kWh天然氣價(jià)格單位用 元/m3。如果功率數(shù)據(jù)是 15 分鐘一個(gè)點(diǎn)那么 T96dt0.25所有 SOC 遞推式都要乘以 dt。很多模型跑出來(lái)的結(jié)果時(shí)間對(duì)不上最后發(fā)現(xiàn)是 dt 漏乘了。一段典型的數(shù)據(jù)初始化代碼T 24; % 優(yōu)化時(shí)段數(shù) dt 1; % 時(shí)間步長(zhǎng)單位小時(shí) % 電負(fù)荷和熱負(fù)荷單位 kW P_load [105 98 92 88 90 96 110 135 150 148 132 128 ... 115 108 120 138 155 160 148 136 120 112 108 100]; H_load [80 75 72 68 66 70 76 82 86 88 80 75 ... 70 68 72 80 90 95 92 88 82 78 74 70]; % 光伏和風(fēng)電預(yù)測(cè)出力單位 kW pv [0 0 0 0 0 5 25 60 90 120 140 145 130 110 90 70 40 10 0 0 0 0 0 0]; wt [18 20 22 24 22 20 18 15 12 10 9 8 10 14 16 18 20 22 24 22 18 16 14 12]; % 分時(shí)電價(jià)單位 元/kWh price_buy [0.40 0.40 0.40 0.40 0.40 0.40 1.20 1.20 1.20 1.20 ... 0.75 0.75 0.75 0.75 1.20 1.20 1.20 1.20 1.20 1.20 ... 0.75 0.75 0.75 0.75]; price_sell repmat(0.40, 1, T); % 上網(wǎng)電價(jià) % CHP 參數(shù) Pchp_max 200; Pchp_min 50; % kW rho_chp 1.2; % 熱電比 a_chp 0.25; % 邊際氣耗m3/kWh b_chp 2.5; % 空載氣耗m3/h gas_price 2.5; % 元/m3 % 電鍋爐 P_eb_max 150; eta_eb 0.95; % 儲(chǔ)電 P_ch_max 150; P_dis_max 150; E_bess_max 600; E_bess_min 60; E_bess_init 120; eta_bess_ch 0.95; eta_bess_dis 0.95; % 儲(chǔ)熱 H_tes_ch_max 200; H_tes_dis_max 200; E_tes_max 1000; E_tes_min 100; E_tes_init 200; eta_tes_ch 0.90; eta_tes_dis 0.90; % 聯(lián)絡(luò)線 P_buy_max 800; P_sell_max 500;這些參數(shù)不需要一次到位第一次跑通后再根據(jù)設(shè)備手冊(cè)細(xì)化。重點(diǎn)是讓數(shù)據(jù)格式統(tǒng)一、單位一致。3.2 基于YALMIP的核心建模代碼變量聲明和約束整合可以這樣寫(xiě)。YALMIP 的sdpvar定義連續(xù)變量binvar定義 0-1 變量optimize負(fù)責(zé)調(diào)用底層求解器。% 決策變量 P_chp sdpvar(1, T, full); H_chp sdpvar(1, T, full); on_chp binvar(1, T, full); P_eb sdpvar(1, T, full); H_eb sdpvar(1, T, full); P_ch sdpvar(1, T, full); P_dis sdpvar(1, T, full); E_bess sdpvar(1, T1, full); H_tes_ch sdpvar(1, T, full); H_tes_dis sdpvar(1, T, full); E_tes sdpvar(1, T1, full); P_buy sdpvar(1, T, full); P_sell sdpvar(1, T, full); Cons []; % CHP 出力區(qū)間與熱電耦合 Cons [Cons, Pchp_min * on_chp P_chp Pchp_max * on_chp]; Cons [Cons, H_chp rho_chp * P_chp]; % 電鍋爐 Cons [Cons, 0 P_eb P_eb_max]; Cons [Cons, H_eb eta_eb * P_eb]; % 儲(chǔ)電約束 Cons [Cons, 0 P_ch P_ch_max]; Cons [Cons, 0 P_dis P_dis_max]; Cons [Cons, E_bess_min E_bess E_bess_max]; Cons [Cons, E_bess(1) E_bess_init]; Cons [Cons, E_bess(T1) E_bess_init]; Cons [Cons, E_bess(2:T1) E_bess(1:T) dt * (eta_bess_ch * P_ch - P_dis / eta_bess_dis)]; % 儲(chǔ)熱約束 Cons [Cons, 0 H_tes_ch H_tes_ch_max]; Cons [Cons, 0 H_tes_dis H_tes_dis_max]; Cons [Cons, E_tes_min E_tes E_tes_max]; Cons [Cons, E_tes(1) E_tes_init]; Cons [Cons, E_tes(T1) E_tes_init]; Cons [Cons, E_tes(2:T1) E_tes(1:T) dt * (eta_tes_ch * H_tes_ch - H_tes_dis / eta_tes_dis)]; % 購(gòu)售電 Cons [Cons, 0 P_buy P_buy_max]; Cons [Cons, 0 P_sell P_sell_max]; % 電功率平衡 Cons [Cons, pv wt P_chp P_dis P_buy P_load P_eb P_ch P_sell]; % 熱功率平衡 Cons [Cons, H_chp H_eb H_tes_dis H_load H_tes_ch]; % CHP 爬坡約束 Cons [Cons, -60 * dt P_chp(2:T) - P_chp(1:T-1) 60 * dt];目標(biāo)函數(shù)里包含購(gòu)電成本、售電收益、天然氣成本、運(yùn)維成本。這里的光伏運(yùn)維成本我按實(shí)際出力pv算模型里可以當(dāng)成常數(shù)不會(huì)影響優(yōu)化變量最優(yōu)解但會(huì)改變總成本數(shù)值。% 目標(biāo)函數(shù)購(gòu)電成本 - 售電收益 天然氣成本 運(yùn)維成本 gas_cost sum((a_chp * P_chp b_chp * on_chp) .* gas_price) * dt; om_cost sum(0.02 * P_chp 0.01 * P_eb 0.008 * (P_ch P_dis) ... 0.005 * (H_tes_ch H_tes_dis) 0.01 * pv) * dt; Objective sum(P_buy .* price_buy) * dt - sum(P_sell .* price_sell) * dt ... gas_cost om_cost; % 求解 options sdpsettings(solver, gurobi, mipgap, 0.001, verbose, 1); tic; sol optimize(Cons, Objective, options); toc;如果沒(méi)有 Gurobi 或者 Cplex把solver改成intlinprog也可以運(yùn)行。只是大型 MILP 問(wèn)題用內(nèi)置的 intlinprog 求解速度會(huì)慢不少我建議至少配置 Gurobi 試一下。3.3 求解器配置與收斂經(jīng)驗(yàn)YALMIP 是建模層真正算 MILP 的是底層求解器。Gurobi 和 Cplex 對(duì)線性整數(shù)規(guī)劃的求解性能很強(qiáng)適合 24 時(shí)段這種規(guī)模。對(duì)純線性規(guī)劃問(wèn)題用linprog就夠但一旦出現(xiàn)binvar必須走整數(shù)規(guī)劃路線。sdpsettings里我慣用的幾個(gè)參數(shù)options sdpsettings(solver, gurobi, mipgap, 0.0005, verbose, 1); options.gurobi.TimeLimit 300; options.gurobi.MIPFocus 1;mipgap設(shè)成 0.0005 意思是允許最優(yōu)解的相對(duì)誤差在 0.05% 以內(nèi)。對(duì)于運(yùn)行成本優(yōu)化這個(gè)精度足夠。TimeLimit是防止模型規(guī)模過(guò)大時(shí)無(wú)限算下去。如果發(fā)現(xiàn)求解時(shí)間和最優(yōu)性之間不平衡可以適當(dāng)放寬 MIPGap比如 1%結(jié)果往往只差幾十塊錢(qián)但求解時(shí)間可能縮短一個(gè)數(shù)量級(jí)。調(diào)試時(shí)建議先用optimize(Cons, [], options)只做可行性求解不優(yōu)化目標(biāo)函數(shù)。如果這個(gè)可行解都找不到說(shuō)明約束之間有矛盾。接下來(lái)用check(Cons)查看約束殘差YALMIP 會(huì)返回每個(gè)約束的 residual一眼能看出哪類約束沖突。4. 仿真結(jié)果怎么看典型日?qǐng)鼍跋碌姆桨笇?duì)比4.1 結(jié)果輸出與繪圖思路求解完成后需要把 YALMIP 的變量取值提取出來(lái)P_chp_opt value(P_chp); H_chp_opt value(H_chp); on_chp_opt value(on_chp); P_eb_opt value(P_eb); P_ch_opt value(P_ch); P_dis_opt value(P_dis); E_bess_opt value(E_bess); H_tes_ch_opt value(H_tes_ch); H_tes_dis_opt value(H_tes_dis); E_tes_opt value(E_tes); P_buy_opt value(P_buy); P_sell_opt value(P_sell);繪圖時(shí)我喜歡用兩張堆疊圖看全天電平衡和熱平衡。電平衡圖橫軸是小時(shí)縱軸是功率堆疊項(xiàng)包括光伏、風(fēng)電、CHP、儲(chǔ)能放電、購(gòu)電熱平衡圖堆疊項(xiàng)包括 CHP 供熱、電鍋爐供熱、蓄熱罐放熱。一張圖能直觀看出各時(shí)段誰(shuí)在出力、誰(shuí)在充電、誰(shuí)在買(mǎi)電。figure; bar(1:T, [pv; wt; P_chp_opt; P_dis_opt; P_buy_opt], stacked); hold on; plot(1:T, P_load, k-, LineWidth, 2); legend(光伏, 風(fēng)電, CHP, 儲(chǔ)電放電, 購(gòu)電, 電負(fù)荷); xlabel(時(shí)段/h); ylabel(功率/kW);4.2 聯(lián)合調(diào)度結(jié)果解讀何時(shí)充電、何時(shí)儲(chǔ)熱、何時(shí)購(gòu)電跑通后的典型日優(yōu)化結(jié)果大致遵循這樣一個(gè)規(guī)律夜間低谷電價(jià)時(shí)段系統(tǒng)傾向于多購(gòu)電一部分直接供電負(fù)荷一部分給蓄電充電同時(shí)電鍋爐在低谷電價(jià)時(shí)段啟動(dòng)制熱并存入蓄熱罐。白天光伏大發(fā)、電負(fù)荷不是特別高時(shí)CHP 可能降出力避免產(chǎn)生多余的熱蓄電在光伏頂峰時(shí)充電、在晚高峰放電。晚高峰電價(jià)最高CHP 滿發(fā)儲(chǔ)電放電蓄熱罐放熱購(gòu)電最少。一個(gè)比較典型的運(yùn)行安排時(shí)段電價(jià)水平CHP電鍋爐儲(chǔ)電儲(chǔ)熱聯(lián)絡(luò)線0:00-6:00谷低載或停機(jī)啟動(dòng)儲(chǔ)熱充電蓄熱大量購(gòu)電7:00-10:00峰滿發(fā)停放電放熱少購(gòu)電11:00-14:00平中低載視光伏而定光伏多余時(shí)充電視熱負(fù)荷而定少量購(gòu)電/售電15:00-20:00峰滿發(fā)停放電放熱少量購(gòu)電21:00-24:00平/谷低載可啟動(dòng)儲(chǔ)熱充電蓄熱購(gòu)電這里需要特別理解“儲(chǔ)熱”的價(jià)值。白天如果光伏較多而電負(fù)荷不夠消化CHP 再滿發(fā)就會(huì)導(dǎo)致多余電力。但如果你把 CHP 在白天多發(fā)一部分電可以賣(mài)給電網(wǎng)附帶的熱量存入蓄熱罐留到晚上用整體賬可能更劃算。這種跨時(shí)段套利關(guān)系只有在電熱聯(lián)合模型里才能被算出來(lái)。5. 調(diào)試過(guò)程中最常踩的坑規(guī)模失控與熱負(fù)荷預(yù)測(cè)誤差5.1 0-1變量爆炸從573秒到2.8秒的重構(gòu)經(jīng)驗(yàn)第一次完整建模時(shí)我抱著“約束越多越嚴(yán)謹(jǐn)”的想法給很多關(guān)系都加了二進(jìn)制變量。除了 CHP 啟停還加了儲(chǔ)電充放互斥、儲(chǔ)熱充放互斥、購(gòu)電售電互斥結(jié)果二進(jìn)制變量數(shù)量變成 4*T96 個(gè)?？雌饋?lái)不多但求解時(shí)間長(zhǎng)達(dá) 573 秒MIPGap 還卡在 1% 下不去。后來(lái)我逐個(gè)刪掉了“經(jīng)濟(jì)上不會(huì)發(fā)生”的互斥變量。儲(chǔ)存設(shè)備同時(shí)充放會(huì)同時(shí)產(chǎn)生兩筆運(yùn)維成本還要承擔(dān)效率損耗優(yōu)化器沒(méi)有任何理由這么做購(gòu)售電同時(shí)發(fā)生在分時(shí)電價(jià)下也不會(huì)省成本。保留 CHP 啟停變量后二進(jìn)制變量降為 24 個(gè)加上約束結(jié)構(gòu)從大 M 型變成更緊湊的邊界約束最終求解時(shí)間降到 2.8 秒。這個(gè)經(jīng)驗(yàn)可以總結(jié)成一條原則每一個(gè)二進(jìn)制變量都要問(wèn)自己它在真實(shí)物理場(chǎng)景中是否必須存在、刪掉它是否真的會(huì)出荒謬解。如果單純?yōu)榱死碚搰?yán)謹(jǐn)而加一堆 0-1 變量最終拖累的是求解性能和模型實(shí)用性。先用線性規(guī)劃跑通再加少量整數(shù)變量是更省時(shí)間的路徑。5.2 熱負(fù)荷預(yù)測(cè)誤差導(dǎo)致的結(jié)果敏感性另一件讓我印象很深的事是熱負(fù)荷預(yù)測(cè)誤差。某天早上我看當(dāng)?shù)貧庀笈_(tái)溫度和歷史熱負(fù)荷數(shù)據(jù)估了一條偏低 20% 的熱負(fù)荷曲線。優(yōu)化結(jié)果出來(lái)蓄熱罐夜間蓄的熱量明顯偏少。到了下午熱負(fù)荷實(shí)測(cè)比預(yù)測(cè)高系統(tǒng)只能臨時(shí)讓 CHP 多發(fā)電、多供熱還把儲(chǔ)電提前放出晚高峰購(gòu)電驟增當(dāng)天實(shí)際成本比用準(zhǔn)確熱負(fù)荷預(yù)測(cè)跑出來(lái)的結(jié)果高了 12%。電負(fù)荷預(yù)測(cè)誤差對(duì)結(jié)果的影響被講得很多但熱負(fù)荷預(yù)測(cè)誤差其實(shí)更隱蔽。因?yàn)闊嵯到y(tǒng)慣性大儲(chǔ)熱時(shí)間常數(shù)長(zhǎng)一旦蓄熱不足后面的補(bǔ)救手段非常有限。我在工程上用了兩個(gè)辦法。第一個(gè)是滾動(dòng)時(shí)域優(yōu)化每隔一小時(shí)用最新實(shí)測(cè)數(shù)據(jù)重新優(yōu)化未來(lái) 24 小時(shí)把預(yù)測(cè)誤差的影響限制在單步之內(nèi)。第二個(gè)是在熱負(fù)荷預(yù)測(cè)值上人為加一個(gè) 5% 到 10% 的安全系數(shù)或者強(qiáng)制蓄熱罐末端 SOC 不低于 20%。這兩種做法都會(huì)讓日運(yùn)行成本略微上升但換回來(lái)的運(yùn)行可靠性提升遠(yuǎn)大于那點(diǎn)成本。6. 從這個(gè)模型還能延伸出什么6.1 加入不確定性魯棒優(yōu)化和場(chǎng)景隨機(jī)優(yōu)化確定性的電熱聯(lián)合調(diào)度模型可以直接用但園區(qū)微電網(wǎng)里光伏、風(fēng)電、熱負(fù)荷都有不確定性。如果想更進(jìn)一步一個(gè)方向是做魯棒優(yōu)化。對(duì)光伏出力和熱負(fù)荷分別設(shè)置盒式不確定集模型在“最壞情況”下仍然滿足電熱平衡。代價(jià)是最優(yōu)成本會(huì)比確定性模型高因?yàn)橄到y(tǒng)要為最壞情況預(yù)留更多備用容量。另一個(gè)更貼近實(shí)際的方向是場(chǎng)景法隨機(jī)優(yōu)化。從歷史數(shù)據(jù)抽 N 個(gè)典型日?qǐng)鼍懊總€(gè)場(chǎng)景分配一個(gè)概率目標(biāo)函數(shù)改成所有場(chǎng)景總成本的期望值。決策變量分成第一階段和第二階段第一階段的 CHP 啟停狀態(tài)、儲(chǔ)電充放決策在場(chǎng)景之間保持一致第二階段功率變量隨場(chǎng)景調(diào)整。這個(gè)模型比魯棒優(yōu)化更精細(xì)但計(jì)算量也翻倍。用 YALMIP 做隨機(jī)優(yōu)化時(shí)可以把場(chǎng)景展開(kāi)成一個(gè)超大的時(shí)間序列用向量索引對(duì)應(yīng)不同場(chǎng)景的位置這樣仍然可以保持單次optimize求解。只是約束矩陣規(guī)模會(huì)隨著場(chǎng)景數(shù)量線性膨脹。6.2 熱網(wǎng)時(shí)空特性建模的擴(kuò)展方向目前的模型把熱力系統(tǒng)看成單節(jié)點(diǎn)所有熱負(fù)荷、產(chǎn)熱設(shè)備、蓄熱罐都匯聚到一個(gè)熱母線。如果園區(qū)里存在幾百米甚至更長(zhǎng)的供熱管網(wǎng)水溫傳輸延遲和管網(wǎng)熱損失就不能再被忽略。一種實(shí)用的擴(kuò)展是把熱網(wǎng)簡(jiǎn)化為帶延遲的蓄熱模型熱源到熱用戶之間的管道相當(dāng)于一個(gè)天然的大蓄熱罐具有管存效應(yīng)。另一種更精確的做法是建模供水溫度和回水溫度用二階錐約束逼近管網(wǎng)水力特性模型從 MILP 變成 MISOCP求解器則需要支持二階錐整數(shù)規(guī)劃。從工程落地角度看單節(jié)點(diǎn)電熱聯(lián)合調(diào)度模型最適合項(xiàng)目前期的容量配置分析和運(yùn)行策略研究。等真正到了系統(tǒng)設(shè)計(jì)階段再把熱網(wǎng)網(wǎng)絡(luò)約束加進(jìn)去兩者在思路上一脈相承。最后分享一個(gè)我自己的土辦法跑完優(yōu)化之后不要直接信任結(jié)果先把各變量取值手動(dòng)帶進(jìn)電平衡和熱平衡等式逐時(shí)段算一遍殘差再用check(Cons)核對(duì)。我遇到過(guò)的大部分異常結(jié)果最后追查下來(lái)都是效率參數(shù)寫(xiě)反了、SOC 索引錯(cuò)位了或者某臺(tái)設(shè)備的單位沒(méi)換算。先固定 CHP 啟停狀態(tài)只做線性規(guī)劃?rùn)z查結(jié)果合理后再放回二進(jìn)制變量是排查模型錯(cuò)誤最省力的順序。這個(gè)習(xí)慣幫我省下過(guò)大量時(shí)間。 電熱聯(lián)合調(diào)度微電網(wǎng)運(yùn)行優(yōu)化模型：MATLAB+YALMIP實(shí)現(xiàn)解析
做微電網(wǎng)運(yùn)行優(yōu)化繞不開(kāi)一個(gè)樸素的問(wèn)題電負(fù)荷和熱負(fù)荷明明共存于同一個(gè)園區(qū)為什么大多數(shù)優(yōu)化模型只盯著“電”在算我最近在某園區(qū)微電網(wǎng)項(xiàng)目的經(jīng)濟(jì)調(diào)度分析里把基于電熱聯(lián)合調(diào)度的微電網(wǎng)運(yùn)行優(yōu)化模型在 MATLAB 里完整跑通了一遍。今天把這個(gè)模型的變量設(shè)計(jì)、約束寫(xiě)法、求解代碼和調(diào)試過(guò)程一起整理出來(lái)適合正在做微電網(wǎng)調(diào)度、綜合能源系統(tǒng)優(yōu)化的研究生和工程師參考尤其是手里有燃?xì)廨啓C(jī)、熱電聯(lián)產(chǎn)機(jī)組、電鍋爐、蓄電和蓄熱設(shè)備想算清楚“一套設(shè)備怎樣組合運(yùn)行一天最省錢(qián)”的這類場(chǎng)景。先說(shuō)明白這個(gè)模型能回答什么問(wèn)題光伏和風(fēng)電出力有波動(dòng)電價(jià)有峰谷天然氣也有價(jià)格園區(qū)既有電負(fù)荷又有熱負(fù)荷電鍋爐可以把電轉(zhuǎn)成熱熱電聯(lián)產(chǎn)機(jī)組CHP可以同時(shí)產(chǎn)電產(chǎn)熱蓄電罐和蓄熱罐可以把能量跨越時(shí)間搬運(yùn)。電熱聯(lián)合調(diào)度就是在滿足電、熱兩類負(fù)荷的前提下決定每一臺(tái)設(shè)備每個(gè)時(shí)段開(kāi)多大功率、充多少放多少讓全天總運(yùn)行成本最低。這本質(zhì)是一個(gè)混合整數(shù)線性規(guī)劃問(wèn)題MATLAB 加 YALMIP 工具箱處理起來(lái)非常順手。1. 電熱聯(lián)合調(diào)度的優(yōu)化空間為什么單算“電”會(huì)吃虧1.1 微電網(wǎng)里的“熱”角色從熱電聯(lián)產(chǎn)到電鍋爐大多數(shù)微電網(wǎng)項(xiàng)目在設(shè)計(jì)初期最容易犯的錯(cuò)是把熱負(fù)荷單獨(dú)扔給燃?xì)忮仩t去扛電負(fù)荷單獨(dú)靠光伏、儲(chǔ)能和電網(wǎng)。這樣做邏輯清晰但經(jīng)濟(jì)性很差。因?yàn)閳@區(qū)里大量熱負(fù)荷存在的地方往往同時(shí)存在電負(fù)荷而 CHP 機(jī)組本身就是“一份燃?xì)馔度胪瑫r(shí)產(chǎn)出電和熱”的設(shè)備。如果你只讓它發(fā)電、不讓它供熱相當(dāng)于花了同樣的燃料錢(qián)白白扔掉了一半以上的可用能量。反過(guò)來(lái)電鍋爐也是一種耦合設(shè)備它消耗電能產(chǎn)出熱能。電鍋爐本身不一定省錢(qián)但在低谷電價(jià)時(shí)段或者在光伏大發(fā)導(dǎo)致棄光時(shí)用電制熱可能比用天然氣制熱便宜得多。于是“電”和“熱”會(huì)通過(guò) CHP 和電鍋爐形成兩條耦合通道。聯(lián)合調(diào)度模型的價(jià)值就是把這兩條通道放進(jìn)同一個(gè)優(yōu)化問(wèn)題里讓系統(tǒng)自己決定何時(shí)走“氣到電到熱”的路徑何時(shí)走“電網(wǎng)電到熱”的路徑。CAT在一個(gè)典型的園區(qū)場(chǎng)景里加入蓄熱罐后熱負(fù)荷甚至可以像電儲(chǔ)能一樣實(shí)現(xiàn)時(shí)移白天光伏多的時(shí)候多產(chǎn)熱存起來(lái)晚上熱負(fù)荷上來(lái)再放熱。這個(gè)過(guò)程單看電力系統(tǒng)模型是完全看不見(jiàn)的但它在實(shí)際運(yùn)行中能顯著降低頂峰時(shí)段購(gòu)電壓力。1.2 一個(gè)直觀案例分產(chǎn)調(diào)度和聯(lián)合調(diào)度的成本差在哪我在調(diào)試模型時(shí)用了兩個(gè)對(duì)照方案。第一個(gè)方案是分產(chǎn)調(diào)度CHP 按“以熱定電”的方式運(yùn)行也就是熱負(fù)荷多少CHP 就發(fā)對(duì)應(yīng)比例的電缺電直接從電網(wǎng)買(mǎi)電鍋爐和蓄熱罐不參與聯(lián)動(dòng)。第二個(gè)方案是電熱聯(lián)合調(diào)度所有設(shè)備統(tǒng)一優(yōu)化。同樣一組典型日負(fù)荷數(shù)據(jù)、同樣一套設(shè)備參數(shù)分產(chǎn)調(diào)度一天運(yùn)行成本大約是 14230 元聯(lián)合調(diào)度是 12350 元節(jié)省了約 13.2%。這 13.2% 主要來(lái)自三個(gè)地方。一是低谷電價(jià)時(shí)段電鍋爐啟動(dòng)儲(chǔ)熱替代了高峰時(shí)段部分燃?xì)庀?。二?CHP 不再被“以熱定電”綁死而是更多在峰電時(shí)段滿發(fā)在低谷時(shí)段壓低出力或者停機(jī)熱缺口由蓄熱罐釋放的熱量補(bǔ)上。三是光伏大發(fā)時(shí)段系統(tǒng)主動(dòng)減少購(gòu)電甚至向電網(wǎng)賣(mài)電把多余電能轉(zhuǎn)化為蓄熱罐里的熱能避免棄光。這三個(gè)收益來(lái)源離開(kāi)電熱聯(lián)合調(diào)度模型都很難被量化。2. 模型的數(shù)學(xué)骨架目標(biāo)函數(shù)、約束與變量設(shè)計(jì)2.1 決策變量的分層設(shè)計(jì)連續(xù)功率變量與啟停0-1變量?jī)?yōu)化模型的第一步是確定決策變量我的做法是分成連續(xù)變量和二進(jìn)制變量?jī)深?。連續(xù)變量包括各設(shè)備每個(gè)時(shí)段的功率和能量狀態(tài)二進(jìn)制變量主要用于 CHP 的啟停狀態(tài)。設(shè)備連續(xù)變量說(shuō)明CHP發(fā)電功率 P_chp、供熱功率 H_chp單位 kWCHP 啟停0-1 變量 on_chp1 表示運(yùn)行電鍋爐耗電功率 P_eb、產(chǎn)熱功率 H_ebH_eb eta_eb * P_eb蓄電充電功率 P_ch、放電功率 P_dis、電量 E_bess電量單位 kWh蓄熱蓄熱功率 H_tes_ch、放熱功率 H_tes_dis、熱量 E_tes熱量單位 kWh聯(lián)絡(luò)線購(gòu)電 P_buy、售電 P_sell單位 kWCHP 的啟停變量不是可選項(xiàng)而是必須項(xiàng)。因?yàn)槿細(xì)廨啓C(jī)存在最小技術(shù)出力比如 50kW如果只建模成連續(xù)變量?jī)?yōu)化器為了省錢(qián)會(huì)把 CHP 壓到 1kW、2kW 這種根本不存在的運(yùn)行點(diǎn)。引入 on_chp 后功率下限寫(xiě)成P_chp_min * on_chp P_chp P_chp_max * on_chp機(jī)組要么不開(kāi)要么至少開(kāi)到 50kW這才符合實(shí)際。一個(gè)需要提前想清楚的點(diǎn)蓄電和蓄熱設(shè)備的“充放互斥”要不要加二進(jìn)制變量。從理論嚴(yán)謹(jǐn)性角度應(yīng)該避免同時(shí)充電和放電。但在實(shí)際模型中由于充放都有損耗且損耗已經(jīng)進(jìn)了目標(biāo)函數(shù)優(yōu)化器在經(jīng)濟(jì)上沒(méi)有動(dòng)機(jī)讓同一臺(tái)儲(chǔ)能同時(shí)充放。所以我通常不加這組二進(jìn)制變量后續(xù)求解速度會(huì)快很多。2.2 電熱平衡約束與CHP可行域約束的建模細(xì)節(jié)電功率平衡是微電網(wǎng)模型的地基任何一個(gè)時(shí)段都必須滿足pv(t) wt(t) P_chp(t) P_dis(t) P_buy(t) P_load(t) P_eb(t) P_ch(t) P_sell(t)左邊是電源側(cè)右邊是負(fù)荷側(cè)。注意電鍋爐的耗電功率 P_eb 必須放在負(fù)荷側(cè)這是初學(xué)者最容易漏掉的一項(xiàng)。光伏、風(fēng)電數(shù)據(jù)如果來(lái)自預(yù)測(cè)曲線可以直接作為已知量放進(jìn)等式。如果考慮棄光就需要給光伏增加一個(gè)“可削減”變量。熱功率平衡相對(duì)簡(jiǎn)單H_chp(t) H_eb(t) H_tes_dis(t) H_load(t) H_tes_ch(t)蓄熱罐充熱時(shí)算是熱負(fù)荷放熱時(shí)算是熱源。CHP 的電氣特性是電熱聯(lián)合調(diào)度里最關(guān)鍵的約束。我用固定熱電比模型也就是H_chp(t) rho_chp * P_chp(t)這種模型用起來(lái)最省事適合做可行性分析。它意味著 CHP 發(fā)多少電就附帶產(chǎn)生多少熱。當(dāng)熱負(fù)荷很小、電負(fù)荷很大時(shí)多余的熱可以進(jìn)蓄熱罐當(dāng)熱負(fù)荷很大、電負(fù)荷很小時(shí)則只能通過(guò)降低 CHP 出力來(lái)減少熱產(chǎn)出缺電部分由電網(wǎng)補(bǔ)。也可以把熱電比關(guān)系改成可行域約束允許熱電比在一定范圍內(nèi)變化但那樣會(huì)引入更復(fù)雜的耦合不等式對(duì)線性求解器不友好。CHP 的天然氣消耗我采用的是線性近似gas_consume(t) a_chp * P_chp(t) b_chp * on_chp(t)其中 a_chp 是邊際氣耗系數(shù)單位是 m3/kWhb_chp 是空載氣耗系數(shù)單位是 m3/h。這兩個(gè)系數(shù)由 CHP 性能曲線擬合得到。這種方法比直接用效率除法更貼近實(shí)際也保留了 MILP 的線性結(jié)構(gòu)。2.3 儲(chǔ)能SOC、爬坡、購(gòu)售電約束怎么列才不踩坑儲(chǔ)能約束最容易踩坑的是 SOC 的索引長(zhǎng)度。如果優(yōu)化周期是 T24充放電功率變量長(zhǎng)度是 24那么 SOC 狀態(tài)變量長(zhǎng)度應(yīng)該是 25因?yàn)橐獜某跏紩r(shí)刻一直遞推到 T 時(shí)刻結(jié)束。我習(xí)慣用 E_bess(1) 表示初始電量E_bess(T1) 表示最終電量中間遞推關(guān)系寫(xiě)成E_bess(t1) E_bess(t) dt * (eta_ch * P_ch(t) - P_dis(t) / eta_dis)充電時(shí)乘以效率放電時(shí)除以效率。這個(gè)細(xì)節(jié)特別重要如果方向?qū)懛磧?yōu)化器會(huì)發(fā)現(xiàn)“充進(jìn)去 1 度電放出來(lái) 1.2 度電”相當(dāng)于憑空創(chuàng)造能量結(jié)果會(huì)出現(xiàn)一堆反直覺(jué)的充放行為。容量約束要加在電量狀態(tài)變量上而不是加在充放功率上E_bess_min E_bess(t) E_bess_max同時(shí)還要給充放電功率分別設(shè)置上限。初始和終值電量可以設(shè)定成同一個(gè)值保證日內(nèi)調(diào)度的可持續(xù)性。如果允許終值高于初始值調(diào)度結(jié)果其實(shí)就是“一天之內(nèi)把儲(chǔ)能充滿”成本會(huì)被人為壓低沒(méi)有實(shí)際參考意義。CHP 的爬坡約束是另一個(gè)容易出問(wèn)題的地方。爬坡定義是相鄰時(shí)段功率變化量-ramp_chp * dt P_chp(t1) - P_chp(t) ramp_chp * dt注意在 MATLAB 中的索引是從 1 開(kāi)始所以代碼里會(huì)有P_chp(2:T) - P_chp(1:T-1)的寫(xiě)法。購(gòu)售電約束方面購(gòu)電和售電功率各設(shè)上限要不要加“不同時(shí)購(gòu)售”的互斥二進(jìn)制變量看實(shí)際需求。我一般先不加因?yàn)榉謺r(shí)電價(jià)下基本不會(huì)出現(xiàn)同時(shí)購(gòu)售的荒謬解。如果最終結(jié)果異常再補(bǔ)一組互斥變量也不遲。3. MATLAB代碼實(shí)現(xiàn)數(shù)據(jù)組織、建模與求解器對(duì)接3.1 參數(shù)與負(fù)荷曲線如何整理成矩陣MATLAB 做優(yōu)化建模第一步就是統(tǒng)一數(shù)據(jù)格式。我強(qiáng)烈建議所有時(shí)序數(shù)據(jù)都用1 x T的行向量組織電價(jià)、負(fù)荷、光伏、風(fēng)電都是如此。這樣后續(xù)與 YALMIP 的變量做點(diǎn)乘運(yùn)算時(shí)可以直接.*不需要來(lái)回轉(zhuǎn)置。單位也要統(tǒng)一。功率單位用 kW能量單位用 kWh電價(jià)單位用 元/kWh天然氣價(jià)格單位用 元/m3。如果功率數(shù)據(jù)是 15 分鐘一個(gè)點(diǎn)那么 T96dt0.25所有 SOC 遞推式都要乘以 dt。很多模型跑出來(lái)的結(jié)果時(shí)間對(duì)不上最后發(fā)現(xiàn)是 dt 漏乘了。一段典型的數(shù)據(jù)初始化代碼T 24; % 優(yōu)化時(shí)段數(shù) dt 1; % 時(shí)間步長(zhǎng)單位小時(shí) % 電負(fù)荷和熱負(fù)荷單位 kW P_load [105 98 92 88 90 96 110 135 150 148 132 128 ... 115 108 120 138 155 160 148 136 120 112 108 100]; H_load [80 75 72 68 66 70 76 82 86 88 80 75 ... 70 68 72 80 90 95 92 88 82 78 74 70]; % 光伏和風(fēng)電預(yù)測(cè)出力單位 kW pv [0 0 0 0 0 5 25 60 90 120 140 145 130 110 90 70 40 10 0 0 0 0 0 0]; wt [18 20 22 24 22 20 18 15 12 10 9 8 10 14 16 18 20 22 24 22 18 16 14 12]; % 分時(shí)電價(jià)單位 元/kWh price_buy [0.40 0.40 0.40 0.40 0.40 0.40 1.20 1.20 1.20 1.20 ... 0.75 0.75 0.75 0.75 1.20 1.20 1.20 1.20 1.20 1.20 ... 0.75 0.75 0.75 0.75]; price_sell repmat(0.40, 1, T); % 上網(wǎng)電價(jià) % CHP 參數(shù) Pchp_max 200; Pchp_min 50; % kW rho_chp 1.2; % 熱電比 a_chp 0.25; % 邊際氣耗m3/kWh b_chp 2.5; % 空載氣耗m3/h gas_price 2.5; % 元/m3 % 電鍋爐 P_eb_max 150; eta_eb 0.95; % 儲(chǔ)電 P_ch_max 150; P_dis_max 150; E_bess_max 600; E_bess_min 60; E_bess_init 120; eta_bess_ch 0.95; eta_bess_dis 0.95; % 儲(chǔ)熱 H_tes_ch_max 200; H_tes_dis_max 200; E_tes_max 1000; E_tes_min 100; E_tes_init 200; eta_tes_ch 0.90; eta_tes_dis 0.90; % 聯(lián)絡(luò)線 P_buy_max 800; P_sell_max 500;這些參數(shù)不需要一次到位第一次跑通后再根據(jù)設(shè)備手冊(cè)細(xì)化。重點(diǎn)是讓數(shù)據(jù)格式統(tǒng)一、單位一致。3.2 基于YALMIP的核心建模代碼變量聲明和約束整合可以這樣寫(xiě)。YALMIP 的sdpvar定義連續(xù)變量binvar定義 0-1 變量optimize負(fù)責(zé)調(diào)用底層求解器。% 決策變量 P_chp sdpvar(1, T, full); H_chp sdpvar(1, T, full); on_chp binvar(1, T, full); P_eb sdpvar(1, T, full); H_eb sdpvar(1, T, full); P_ch sdpvar(1, T, full); P_dis sdpvar(1, T, full); E_bess sdpvar(1, T1, full); H_tes_ch sdpvar(1, T, full); H_tes_dis sdpvar(1, T, full); E_tes sdpvar(1, T1, full); P_buy sdpvar(1, T, full); P_sell sdpvar(1, T, full); Cons []; % CHP 出力區(qū)間與熱電耦合 Cons [Cons, Pchp_min * on_chp P_chp Pchp_max * on_chp]; Cons [Cons, H_chp rho_chp * P_chp]; % 電鍋爐 Cons [Cons, 0 P_eb P_eb_max]; Cons [Cons, H_eb eta_eb * P_eb]; % 儲(chǔ)電約束 Cons [Cons, 0 P_ch P_ch_max]; Cons [Cons, 0 P_dis P_dis_max]; Cons [Cons, E_bess_min E_bess E_bess_max]; Cons [Cons, E_bess(1) E_bess_init]; Cons [Cons, E_bess(T1) E_bess_init]; Cons [Cons, E_bess(2:T1) E_bess(1:T) dt * (eta_bess_ch * P_ch - P_dis / eta_bess_dis)]; % 儲(chǔ)熱約束 Cons [Cons, 0 H_tes_ch H_tes_ch_max]; Cons [Cons, 0 H_tes_dis H_tes_dis_max]; Cons [Cons, E_tes_min E_tes E_tes_max]; Cons [Cons, E_tes(1) E_tes_init]; Cons [Cons, E_tes(T1) E_tes_init]; Cons [Cons, E_tes(2:T1) E_tes(1:T) dt * (eta_tes_ch * H_tes_ch - H_tes_dis / eta_tes_dis)]; % 購(gòu)售電 Cons [Cons, 0 P_buy P_buy_max]; Cons [Cons, 0 P_sell P_sell_max]; % 電功率平衡 Cons [Cons, pv wt P_chp P_dis P_buy P_load P_eb P_ch P_sell]; % 熱功率平衡 Cons [Cons, H_chp H_eb H_tes_dis H_load H_tes_ch]; % CHP 爬坡約束 Cons [Cons, -60 * dt P_chp(2:T) - P_chp(1:T-1) 60 * dt];目標(biāo)函數(shù)里包含購(gòu)電成本、售電收益、天然氣成本、運(yùn)維成本。這里的光伏運(yùn)維成本我按實(shí)際出力pv算模型里可以當(dāng)成常數(shù)不會(huì)影響優(yōu)化變量最優(yōu)解但會(huì)改變總成本數(shù)值。% 目標(biāo)函數(shù)購(gòu)電成本 - 售電收益 天然氣成本 運(yùn)維成本 gas_cost sum((a_chp * P_chp b_chp * on_chp) .* gas_price) * dt; om_cost sum(0.02 * P_chp 0.01 * P_eb 0.008 * (P_ch P_dis) ... 0.005 * (H_tes_ch H_tes_dis) 0.01 * pv) * dt; Objective sum(P_buy .* price_buy) * dt - sum(P_sell .* price_sell) * dt ... gas_cost om_cost; % 求解 options sdpsettings(solver, gurobi, mipgap, 0.001, verbose, 1); tic; sol optimize(Cons, Objective, options); toc;如果沒(méi)有 Gurobi 或者 Cplex把solver改成intlinprog也可以運(yùn)行。只是大型 MILP 問(wèn)題用內(nèi)置的 intlinprog 求解速度會(huì)慢不少我建議至少配置 Gurobi 試一下。3.3 求解器配置與收斂經(jīng)驗(yàn)YALMIP 是建模層真正算 MILP 的是底層求解器。Gurobi 和 Cplex 對(duì)線性整數(shù)規(guī)劃的求解性能很強(qiáng)適合 24 時(shí)段這種規(guī)模。對(duì)純線性規(guī)劃問(wèn)題用linprog就夠但一旦出現(xiàn)binvar必須走整數(shù)規(guī)劃路線。sdpsettings里我慣用的幾個(gè)參數(shù)options sdpsettings(solver, gurobi, mipgap, 0.0005, verbose, 1); options.gurobi.TimeLimit 300; options.gurobi.MIPFocus 1;mipgap設(shè)成 0.0005 意思是允許最優(yōu)解的相對(duì)誤差在 0.05% 以內(nèi)。對(duì)于運(yùn)行成本優(yōu)化這個(gè)精度足夠。TimeLimit是防止模型規(guī)模過(guò)大時(shí)無(wú)限算下去。如果發(fā)現(xiàn)求解時(shí)間和最優(yōu)性之間不平衡可以適當(dāng)放寬 MIPGap比如 1%結(jié)果往往只差幾十塊錢(qián)但求解時(shí)間可能縮短一個(gè)數(shù)量級(jí)。調(diào)試時(shí)建議先用optimize(Cons, [], options)只做可行性求解不優(yōu)化目標(biāo)函數(shù)。如果這個(gè)可行解都找不到說(shuō)明約束之間有矛盾。接下來(lái)用check(Cons)查看約束殘差YALMIP 會(huì)返回每個(gè)約束的 residual一眼能看出哪類約束沖突。4. 仿真結(jié)果怎么看典型日?qǐng)鼍跋碌姆桨笇?duì)比4.1 結(jié)果輸出與繪圖思路求解完成后需要把 YALMIP 的變量取值提取出來(lái)P_chp_opt value(P_chp); H_chp_opt value(H_chp); on_chp_opt value(on_chp); P_eb_opt value(P_eb); P_ch_opt value(P_ch); P_dis_opt value(P_dis); E_bess_opt value(E_bess); H_tes_ch_opt value(H_tes_ch); H_tes_dis_opt value(H_tes_dis); E_tes_opt value(E_tes); P_buy_opt value(P_buy); P_sell_opt value(P_sell);繪圖時(shí)我喜歡用兩張堆疊圖看全天電平衡和熱平衡。電平衡圖橫軸是小時(shí)縱軸是功率堆疊項(xiàng)包括光伏、風(fēng)電、CHP、儲(chǔ)能放電、購(gòu)電熱平衡圖堆疊項(xiàng)包括 CHP 供熱、電鍋爐供熱、蓄熱罐放熱。一張圖能直觀看出各時(shí)段誰(shuí)在出力、誰(shuí)在充電、誰(shuí)在買(mǎi)電。figure; bar(1:T, [pv; wt; P_chp_opt; P_dis_opt; P_buy_opt], stacked); hold on; plot(1:T, P_load, k-, LineWidth, 2); legend(光伏, 風(fēng)電, CHP, 儲(chǔ)電放電, 購(gòu)電, 電負(fù)荷); xlabel(時(shí)段/h); ylabel(功率/kW);4.2 聯(lián)合調(diào)度結(jié)果解讀何時(shí)充電、何時(shí)儲(chǔ)熱、何時(shí)購(gòu)電跑通后的典型日優(yōu)化結(jié)果大致遵循這樣一個(gè)規(guī)律夜間低谷電價(jià)時(shí)段系統(tǒng)傾向于多購(gòu)電一部分直接供電負(fù)荷一部分給蓄電充電同時(shí)電鍋爐在低谷電價(jià)時(shí)段啟動(dòng)制熱并存入蓄熱罐。白天光伏大發(fā)、電負(fù)荷不是特別高時(shí)CHP 可能降出力避免產(chǎn)生多余的熱蓄電在光伏頂峰時(shí)充電、在晚高峰放電。晚高峰電價(jià)最高CHP 滿發(fā)儲(chǔ)電放電蓄熱罐放熱購(gòu)電最少。一個(gè)比較典型的運(yùn)行安排時(shí)段電價(jià)水平CHP電鍋爐儲(chǔ)電儲(chǔ)熱聯(lián)絡(luò)線0:00-6:00谷低載或停機(jī)啟動(dòng)儲(chǔ)熱充電蓄熱大量購(gòu)電7:00-10:00峰滿發(fā)停放電放熱少購(gòu)電11:00-14:00平中低載視光伏而定光伏多余時(shí)充電視熱負(fù)荷而定少量購(gòu)電/售電15:00-20:00峰滿發(fā)停放電放熱少量購(gòu)電21:00-24:00平/谷低載可啟動(dòng)儲(chǔ)熱充電蓄熱購(gòu)電這里需要特別理解“儲(chǔ)熱”的價(jià)值。白天如果光伏較多而電負(fù)荷不夠消化CHP 再滿發(fā)就會(huì)導(dǎo)致多余電力。但如果你把 CHP 在白天多發(fā)一部分電可以賣(mài)給電網(wǎng)附帶的熱量存入蓄熱罐留到晚上用整體賬可能更劃算。這種跨時(shí)段套利關(guān)系只有在電熱聯(lián)合模型里才能被算出來(lái)。5. 調(diào)試過(guò)程中最常踩的坑規(guī)模失控與熱負(fù)荷預(yù)測(cè)誤差5.1 0-1變量爆炸從573秒到2.8秒的重構(gòu)經(jīng)驗(yàn)第一次完整建模時(shí)我抱著“約束越多越嚴(yán)謹(jǐn)”的想法給很多關(guān)系都加了二進(jìn)制變量。除了 CHP 啟停還加了儲(chǔ)電充放互斥、儲(chǔ)熱充放互斥、購(gòu)電售電互斥結(jié)果二進(jìn)制變量數(shù)量變成 4*T96 個(gè)?？雌饋?lái)不多但求解時(shí)間長(zhǎng)達(dá) 573 秒MIPGap 還卡在 1% 下不去。后來(lái)我逐個(gè)刪掉了“經(jīng)濟(jì)上不會(huì)發(fā)生”的互斥變量。儲(chǔ)存設(shè)備同時(shí)充放會(huì)同時(shí)產(chǎn)生兩筆運(yùn)維成本還要承擔(dān)效率損耗優(yōu)化器沒(méi)有任何理由這么做購(gòu)售電同時(shí)發(fā)生在分時(shí)電價(jià)下也不會(huì)省成本。保留 CHP 啟停變量后二進(jìn)制變量降為 24 個(gè)加上約束結(jié)構(gòu)從大 M 型變成更緊湊的邊界約束最終求解時(shí)間降到 2.8 秒。這個(gè)經(jīng)驗(yàn)可以總結(jié)成一條原則每一個(gè)二進(jìn)制變量都要問(wèn)自己它在真實(shí)物理場(chǎng)景中是否必須存在、刪掉它是否真的會(huì)出荒謬解。如果單純?yōu)榱死碚搰?yán)謹(jǐn)而加一堆 0-1 變量最終拖累的是求解性能和模型實(shí)用性。先用線性規(guī)劃跑通再加少量整數(shù)變量是更省時(shí)間的路徑。5.2 熱負(fù)荷預(yù)測(cè)誤差導(dǎo)致的結(jié)果敏感性另一件讓我印象很深的事是熱負(fù)荷預(yù)測(cè)誤差。某天早上我看當(dāng)?shù)貧庀笈_(tái)溫度和歷史熱負(fù)荷數(shù)據(jù)估了一條偏低 20% 的熱負(fù)荷曲線。優(yōu)化結(jié)果出來(lái)蓄熱罐夜間蓄的熱量明顯偏少。到了下午熱負(fù)荷實(shí)測(cè)比預(yù)測(cè)高系統(tǒng)只能臨時(shí)讓 CHP 多發(fā)電、多供熱還把儲(chǔ)電提前放出晚高峰購(gòu)電驟增當(dāng)天實(shí)際成本比用準(zhǔn)確熱負(fù)荷預(yù)測(cè)跑出來(lái)的結(jié)果高了 12%。電負(fù)荷預(yù)測(cè)誤差對(duì)結(jié)果的影響被講得很多但熱負(fù)荷預(yù)測(cè)誤差其實(shí)更隱蔽。因?yàn)闊嵯到y(tǒng)慣性大儲(chǔ)熱時(shí)間常數(shù)長(zhǎng)一旦蓄熱不足后面的補(bǔ)救手段非常有限。我在工程上用了兩個(gè)辦法。第一個(gè)是滾動(dòng)時(shí)域優(yōu)化每隔一小時(shí)用最新實(shí)測(cè)數(shù)據(jù)重新優(yōu)化未來(lái) 24 小時(shí)把預(yù)測(cè)誤差的影響限制在單步之內(nèi)。第二個(gè)是在熱負(fù)荷預(yù)測(cè)值上人為加一個(gè) 5% 到 10% 的安全系數(shù)或者強(qiáng)制蓄熱罐末端 SOC 不低于 20%。這兩種做法都會(huì)讓日運(yùn)行成本略微上升但換回來(lái)的運(yùn)行可靠性提升遠(yuǎn)大于那點(diǎn)成本。6. 從這個(gè)模型還能延伸出什么6.1 加入不確定性魯棒優(yōu)化和場(chǎng)景隨機(jī)優(yōu)化確定性的電熱聯(lián)合調(diào)度模型可以直接用但園區(qū)微電網(wǎng)里光伏、風(fēng)電、熱負(fù)荷都有不確定性。如果想更進(jìn)一步一個(gè)方向是做魯棒優(yōu)化。對(duì)光伏出力和熱負(fù)荷分別設(shè)置盒式不確定集模型在“最壞情況”下仍然滿足電熱平衡。代價(jià)是最優(yōu)成本會(huì)比確定性模型高因?yàn)橄到y(tǒng)要為最壞情況預(yù)留更多備用容量。另一個(gè)更貼近實(shí)際的方向是場(chǎng)景法隨機(jī)優(yōu)化。從歷史數(shù)據(jù)抽 N 個(gè)典型日?qǐng)鼍懊總€(gè)場(chǎng)景分配一個(gè)概率目標(biāo)函數(shù)改成所有場(chǎng)景總成本的期望值。決策變量分成第一階段和第二階段第一階段的 CHP 啟停狀態(tài)、儲(chǔ)電充放決策在場(chǎng)景之間保持一致第二階段功率變量隨場(chǎng)景調(diào)整。這個(gè)模型比魯棒優(yōu)化更精細(xì)但計(jì)算量也翻倍。用 YALMIP 做隨機(jī)優(yōu)化時(shí)可以把場(chǎng)景展開(kāi)成一個(gè)超大的時(shí)間序列用向量索引對(duì)應(yīng)不同場(chǎng)景的位置這樣仍然可以保持單次optimize求解。只是約束矩陣規(guī)模會(huì)隨著場(chǎng)景數(shù)量線性膨脹。6.2 熱網(wǎng)時(shí)空特性建模的擴(kuò)展方向目前的模型把熱力系統(tǒng)看成單節(jié)點(diǎn)所有熱負(fù)荷、產(chǎn)熱設(shè)備、蓄熱罐都匯聚到一個(gè)熱母線。如果園區(qū)里存在幾百米甚至更長(zhǎng)的供熱管網(wǎng)水溫傳輸延遲和管網(wǎng)熱損失就不能再被忽略。一種實(shí)用的擴(kuò)展是把熱網(wǎng)簡(jiǎn)化為帶延遲的蓄熱模型熱源到熱用戶之間的管道相當(dāng)于一個(gè)天然的大蓄熱罐具有管存效應(yīng)。另一種更精確的做法是建模供水溫度和回水溫度用二階錐約束逼近管網(wǎng)水力特性模型從 MILP 變成 MISOCP求解器則需要支持二階錐整數(shù)規(guī)劃。從工程落地角度看單節(jié)點(diǎn)電熱聯(lián)合調(diào)度模型最適合項(xiàng)目前期的容量配置分析和運(yùn)行策略研究。等真正到了系統(tǒng)設(shè)計(jì)階段再把熱網(wǎng)網(wǎng)絡(luò)約束加進(jìn)去兩者在思路上一脈相承。最后分享一個(gè)我自己的土辦法跑完優(yōu)化之后不要直接信任結(jié)果先把各變量取值手動(dòng)帶進(jìn)電平衡和熱平衡等式逐時(shí)段算一遍殘差再用check(Cons)核對(duì)。我遇到過(guò)的大部分異常結(jié)果最后追查下來(lái)都是效率參數(shù)寫(xiě)反了、SOC 索引錯(cuò)位了或者某臺(tái)設(shè)備的單位沒(méi)換算。先固定 CHP 啟停狀態(tài)只做線性規(guī)劃?rùn)z查結(jié)果合理后再放回二進(jìn)制變量是排查模型錯(cuò)誤最省力的順序。這個(gè)習(xí)慣幫我省下過(guò)大量時(shí)間。