機(jī)組組合優(yōu)化實(shí)戰(zhàn):MILP建模與Python求解全指南)
簡介一套MATLAB電力系統(tǒng)機(jī)組組合優(yōu)化實(shí)踐資料面向電力系統(tǒng)調(diào)度、運(yùn)行規(guī)劃研究人員及相關(guān)專業(yè)學(xué)生重點(diǎn)覆蓋機(jī)組啟停決策、最優(yōu)出力分配與多類約束處理。壓縮包共10個(gè)文件以Excel結(jié)果表格、Visio圖表、MATLAB腳本和Word說明文檔為主整體312KB內(nèi)容包括熱備用系數(shù)0.05與0.2兩種場景下的機(jī)組組合求解結(jié)果、直流潮流下節(jié)點(diǎn)導(dǎo)納矩陣及其逆矩陣、核心優(yōu)化程序與基本要求說明圖表和數(shù)據(jù)相互對應(yīng)便于完整復(fù)盤建模與求解流程。已有304人學(xué)習(xí)下載。借助該包可理解機(jī)組組合的數(shù)學(xué)建模思路運(yùn)行MATLAB腳本復(fù)現(xiàn)優(yōu)化過程對比不同熱備用水平對機(jī)組出力的影響數(shù)據(jù)與代碼銜接完整主程序可直接運(yùn)行并輸出結(jié)果表格適合用于課程設(shè)計(jì)、畢業(yè)設(shè)計(jì)或?qū)嶋H調(diào)度分析。1. 電力系統(tǒng)機(jī)組組合優(yōu)化是什么一份zip包背后的調(diào)度難題這份編號為034的zip包里裝的是什么其實(shí)不用解壓也能猜到八成——電力系統(tǒng)機(jī)組組合優(yōu)化Unit Commitment簡稱UC。它解決的是每個(gè)調(diào)度員每天都要做兩次的問題未來24小時(shí)哪些機(jī)組該開、哪些該關(guān)、每臺發(fā)多少電才能在滿足負(fù)荷和備用要求的前提下把總成本壓到最低。別小看這個(gè)“開開關(guān)關(guān)”寫成混合整數(shù)規(guī)劃后規(guī)??梢詮氖畮着_機(jī)組膨脹到上千臺求解時(shí)間從秒級變成小時(shí)級。這篇筆記適合三類人電力系統(tǒng)專業(yè)學(xué)生想找一份能跑通、能改參數(shù)的學(xué)習(xí)案例售電公司或園區(qū)微網(wǎng)的運(yùn)行工程師打算把調(diào)度方案從拍腦袋換成優(yōu)化求解剛?cè)雰?yōu)化工具、想拿真實(shí)業(yè)務(wù)練手的算法工程師。我會(huì)從建模、代碼、調(diào)參一路講到踩坑把這類zip包背后最常用的落地路徑一次說透。2. 機(jī)組組合建模與MILP選型從經(jīng)濟(jì)調(diào)度到啟停決策2.1 目標(biāo)函數(shù)里最容易算漏的三項(xiàng)成本機(jī)組組合和經(jīng)濟(jì)調(diào)度最大的區(qū)別是多了“啟?!边@個(gè)整數(shù)決策。經(jīng)濟(jì)調(diào)度默認(rèn)機(jī)組已經(jīng)在線上只分配出力UC要把機(jī)組開起來、關(guān)掉還得決定什么時(shí)刻操作。所以目標(biāo)函數(shù)至少包含三塊燃料成本、啟動(dòng)成本和停機(jī)成本。燃料成本通常用出力的二次函數(shù)擬合工程上為了保持MILP線性結(jié)構(gòu)會(huì)做成三段分段線性函數(shù)。啟動(dòng)成本一般細(xì)分為冷啟動(dòng)、溫啟動(dòng)、熱啟動(dòng)三種價(jià)格相差很大取決于停機(jī)時(shí)長。停機(jī)成本雖然比啟動(dòng)成本小但一些大機(jī)組關(guān)停時(shí)要吹掃、要消耗輔助燃料算總賬時(shí)漏掉它結(jié)果就偏向“頻繁啟停”的假便宜方案。很多入門案例只寫燃料成本求解器為了省啟動(dòng)費(fèi)就會(huì)瘋狂啟停機(jī)組最后給出的總成本很低但調(diào)度員一看啟動(dòng)次數(shù)就知道不能執(zhí)行。這也是新手對不上實(shí)際成本最常見的原因——目標(biāo)函數(shù)不完整優(yōu)化出來的方案天然偏門。2.2 五類約束功率平衡、備用、爬坡、最小啟停、初始狀態(tài)UC的標(biāo)準(zhǔn)約束可以歸成五組。第一功率平衡約束每個(gè)時(shí)段所有在線機(jī)組出力之和必須等于負(fù)荷這是硬等式。第二備用約束系統(tǒng)要預(yù)留旋轉(zhuǎn)備用一般寫成在線機(jī)組最大出力之和大于等于“負(fù)荷備用容量”。第三爬坡約束相鄰時(shí)段之間出力變化不能超過爬坡率。第四最小運(yùn)行和最小停機(jī)時(shí)間約束開機(jī)后至少要維持若干小時(shí)停機(jī)后也要冷卻若干小時(shí)不能頻繁啟停。第五初始狀態(tài)約束調(diào)度起始時(shí)刻之前各機(jī)組是開著還是關(guān)著、已連續(xù)運(yùn)行多久停機(jī)多久都要作為參數(shù)傳進(jìn)去。前四類好理解初始狀態(tài)是新手最容易爆雷的地方。比如某臺機(jī)組在調(diào)度起跑點(diǎn)之前已經(jīng)連續(xù)運(yùn)行了20小時(shí)最小運(yùn)行時(shí)間是8小時(shí)那它在第一個(gè)時(shí)段完全不用額外限制但如果這臺機(jī)組前一天剛停機(jī)2小時(shí)而最小停機(jī)時(shí)間是8小時(shí)那么從調(diào)度開始的前6個(gè)小時(shí)里它都不能重新開機(jī)。如果沒有這些初始信息模型會(huì)把第一時(shí)段當(dāng)成“起點(diǎn)”算出違反物理約束的開機(jī)計(jì)劃。這類錯(cuò)誤在結(jié)果里很難排查因?yàn)榭偝杀究粗@鰜碜鰰r(shí)序回放才發(fā)現(xiàn)前面的機(jī)組狀態(tài)全是錯(cuò)的。2.3 為什么選MILP啟發(fā)式算法的兩個(gè)硬傷很多人問“這類問題用遺傳算法不也行嗎”但電力系統(tǒng)行業(yè)從二十年前就開始轉(zhuǎn)向MILP到現(xiàn)在商業(yè)求解器和開源CBC、SCIP都能處理幾十萬變量的問題原因有兩點(diǎn)。第一啟發(fā)式算法給不出最優(yōu)性差距跑完一組結(jié)果你只知道“還行”不知道離最優(yōu)解差百分之幾MILP求解器明確告訴你gap調(diào)度員可以直接拍板“gap 1%以內(nèi)就用”。第二UC帶大量等式約束和邏輯約束啟發(fā)式算法處理約束只能靠罰函數(shù)罰系數(shù)調(diào)起來是玄學(xué)一不小心功率平衡就根本不滿足。MILP把約束寫進(jìn)模型可行域由求解器保證解出來就能直接執(zhí)行。以pulp為例目標(biāo)函數(shù)寫起來很直觀import pulp as lp # prob是問題對象u[i,t]是0/1開機(jī)狀態(tài)p[i,t]是連續(xù)出力 # start[i,t]、shutdown[i,t]分別表示t時(shí)段是否發(fā)生啟動(dòng)和停機(jī) prob lp.lpSum( fuel_cost[i] * p[i, t] start_cost[i] * start[i, t] shutdown_cost[i] * shutdown[i, t] for i in units for t in periods )這段代碼有三個(gè)關(guān)鍵點(diǎn)需要說明。第一start[i,t]和shutdown[i,t]是二進(jìn)制變量它們不是獨(dú)立自由變量必須通過狀態(tài)轉(zhuǎn)移約束與u[i,t]關(guān)聯(lián)否則求解器會(huì)“白嫖”啟動(dòng)成本——把每個(gè)時(shí)段都定義成啟動(dòng)但實(shí)際機(jī)組狀態(tài)根本沒變化。第二燃料成本系數(shù)這里用的是線性系數(shù)如果原始數(shù)據(jù)是二次曲線要提前做分段線性化。第三shutdown_cost容易漏我見過好幾個(gè)模型只有啟動(dòng)成本沒有停機(jī)成本優(yōu)化結(jié)果里機(jī)組啟停次數(shù)明顯偏高。為什么強(qiáng)調(diào)MILP而不是動(dòng)態(tài)規(guī)劃動(dòng)態(tài)規(guī)劃做小規(guī)模UC很精確但狀態(tài)爆炸太快機(jī)組超過50臺基本沒法用拉格朗日松弛可以解決更大規(guī)模但對偶間隙恢復(fù)可行解的代碼復(fù)雜度高一般團(tuán)隊(duì)很難短時(shí)間搞定。因此從標(biāo)準(zhǔn)做法看MILP是成熟的工業(yè)選擇后續(xù)換商業(yè)求解器也只需改一行參數(shù)。3. 把UC跑起來用python和開源求解器求解最小算例3.1 準(zhǔn)備數(shù)據(jù)四臺機(jī)組參數(shù)表與負(fù)荷曲線討論理論夠了先做一個(gè)能跑出結(jié)果的最小算例。四臺機(jī)組24個(gè)時(shí)段負(fù)荷曲線直接寫在腳本里。機(jī)組參數(shù)如下表燃料成本按一次線性近似單位是元/MWh方便計(jì)算。機(jī)組Pmin (MW)Pmax (MW)爬坡率 (MW/h)最小運(yùn)行 (h)最小停機(jī) (h)啟動(dòng)成本 (元)燃料成本 (元/MWh)G110040012042450060G28030010032300075G3502008022180095G4301005011800120為了簡單不單獨(dú)列停機(jī)成本在啟動(dòng)成本里按50%折算。負(fù)荷曲線用一段Python列表從凌晨低谷的350MW升到晚間高峰的850MW模擬典型的雙峰形狀。3.2 用pulp建模機(jī)組組合核心代碼與邏輯說明先定義決策變量u[i,t]表示機(jī)組i在t時(shí)段是否在線p[i,t]是出力start[i,t]和shutdown[i,t]分別記錄啟動(dòng)和停機(jī)動(dòng)作。pulp里二進(jìn)制變量直接指定catBinary。特別注意p[i,t]不能只設(shè)上下界還要與u[i,t]耦合機(jī)組離線時(shí)出力必須為0。import pulp as lp T 24 units { 1: {pmin: 100, pmax: 400, ramp: 120, minup: 4, mindown: 2, sc: 4500, fc: 60}, 2: {pmin: 80, pmax: 300, ramp: 100, minup: 3, mindown: 2, sc: 3000, fc: 75}, 3: {pmin: 50, pmax: 200, ramp: 80, minup: 2, mindown: 2, sc: 1800, fc: 95}, 4: {pmin: 30, pmax: 100, ramp: 50, minup: 1, mindown: 1, sc: 800, fc: 120}, } load [ 350, 340, 330, 350, 400, 500, 650, 750, 780, 790, 780, 730, 700, 680, 670, 690, 730, 760, 800, 850, 820, 700, 520, 380, ] prob lp.LpProblem(UC_Demo, lp.LpMinimize) u {}, p {}, start {}, shutdown {} for i in units: for t in range(T): u[i, t] lp.LpVariable(fu_{i}_{t}, catBinary) p[i, t] lp.LpVariable(fp_{i}_{t}, lowBoundunits[i][pmin], upBoundunits[i][pmax]) start[i, t] lp.LpVariable(fs_{i}_{t}, catBinary) shutdown[i, t] lp.LpVariable(fd_{i}_{t}, catBinary) # 目標(biāo)函數(shù)燃料成本 啟動(dòng)成本 0.5*啟動(dòng)成本作為停機(jī)成本近似 prob lp.lpSum( units[i][fc] * p[i, t] 1.5 * units[i][sc] * start[i, t] for i in units for t in range(T) ) # 功率平衡 for t in range(T): prob lp.lpSum(p[i, t] for i in units) load[t] # 啟停與出力耦合離線時(shí)出力為0在線時(shí)出力在Pmin~Pmax之間 for i in units: for t in range(T): prob p[i, t] units[i][pmax] * u[i, t] prob p[i, t] units[i][pmin] * u[i, t] # 狀態(tài)轉(zhuǎn)移定義啟動(dòng)和停機(jī)動(dòng)作 for i in units: for t in range(1, T): prob start[i, t] u[i, t] - u[i, t-1] prob shutdown[i, t] u[i, t-1] - u[i, t] # 最小運(yùn)行時(shí)間一旦啟動(dòng)至少連續(xù)運(yùn)行minup個(gè)時(shí)段 for i in units: mu units[i][minup] for t in range(T - mu 1): prob lp.lpSum(u[i, tt] for tt in range(t, t mu)) mu * start[i, t] # 最小停機(jī)時(shí)間一旦停機(jī)至少連續(xù)停機(jī)mindown個(gè)時(shí)段 for i in units: md units[i][mindown] for t in range(T - md 1): prob lp.lpSum(1 - u[i, tt] for tt in range(t, t md)) md * shutdown[i, t] # 爬坡約束 for i in units: for t in range(1, T): prob p[i, t] - p[i, t-1] units[i][ramp] prob p[i, t-1] - p[i, t] units[i][ramp] status prob.solve(lp.PULP_CBC_CMD(msg1, timeLimit30)) print(求解狀態(tài):, lp.LpStatus[status]) print(總成本:, lp.value(prob.objective))核心邏輯說明三點(diǎn)。一是逆變電器耦合約束p[i,t] pmax * u[i,t]保證停機(jī)時(shí)機(jī)組出力被鎖到0p[i,t] pmin * u[i,t]保證在線機(jī)組不能低于技術(shù)出力下限。少了這條求解器可能讓一臺離線機(jī)組按Pmin出力功率平衡雖然滿足但狀態(tài)和出力自相矛盾。二是最小運(yùn)行時(shí)間的標(biāo)準(zhǔn)寫法累計(jì)開機(jī)時(shí)長必須大于等于minup * start[i,t]這是一個(gè)經(jīng)典的MILP建模技巧。三是啟動(dòng)和停機(jī)動(dòng)作與狀態(tài)的關(guān)系上面的寫法是單側(cè)松約束配合目標(biāo)函數(shù)最小化求解器不會(huì)故意多報(bào)啟動(dòng)動(dòng)作。跑這個(gè)算例時(shí)如果CBC在30秒內(nèi)沒達(dá)到gap 0通常也能返回一個(gè)可行解。四臺機(jī)組規(guī)模很小一般幾秒內(nèi)就能找到最優(yōu)解。3.3 求解參數(shù)設(shè)置與結(jié)果解讀把timeLimit設(shè)成30秒是為了演示參數(shù)入口。真正跑的時(shí)候四臺機(jī)組即使不用時(shí)間限制也很快但等到后面做幾十臺機(jī)組時(shí)這行參數(shù)就是救命稻草。結(jié)果解讀建議額外打印每個(gè)時(shí)段的機(jī)組組合和出力明細(xì)for t in range(T): gen {i: lp.value(p[i, t]) for i in units if lp.value(u[i, t]) 0.5} total sum(gen.values()) print(ft{t1:2d} 負(fù)荷{load[t]:4.0f} 出力{total:4.0f} f在線{list(gen.keys())})判斷結(jié)果好壞先看兩件事第一每個(gè)時(shí)段出力總和必須等于負(fù)荷差值超過0.001就說明模型約束寫錯(cuò)了第二啟動(dòng)次數(shù)是否合理比如四臺機(jī)組24小時(shí)啟動(dòng)總數(shù)一般在3~8次如果出現(xiàn)十幾二十次就要回頭檢查目標(biāo)函數(shù)里是否漏了啟動(dòng)成本。4. 從玩具算例到真調(diào)度求解性能與參數(shù)調(diào)優(yōu)4.1 MIP gap不是越小越好給求解器留余地很多人在MILP里只會(huì)設(shè)gapRel0.0001然后等幾個(gè)小時(shí)。實(shí)際調(diào)度場景根本沒必要。負(fù)荷預(yù)測本身有±3%的誤差即使給你一個(gè)gap為0.001%的完美解執(zhí)行時(shí)也會(huì)被預(yù)測誤差吃掉。行業(yè)內(nèi)做日前計(jì)劃gap在0.5%到1%就足夠日內(nèi)滾動(dòng)計(jì)劃可以放寬到2%甚至5%。gapRel參數(shù)本質(zhì)是“相對最優(yōu)性間隙”公式是(best_bound - best_solution) / best_bound。設(shè)置1%不是“允許誤差1%”而是告訴求解器如果當(dāng)前解距離最優(yōu)界在1%以內(nèi)立刻停止。這一點(diǎn)新手上手時(shí)容易誤讀以為設(shè)了1%就等于1%的誤差不可控。事實(shí)上求解器給出的解是可行解約束全部滿足只是不保證繼續(xù)改進(jìn)一定超過1%。如果想看求解器當(dāng)前最優(yōu)界可以用solver lp.PULP_CBC_CMD(msg0, timeLimit300, options[gapRel 0.005]) prob.solve(solver) # bestBound在pulp中通常存在但某些版本可能取不到 print(目標(biāo)值:, lp.value(prob.objective))這是對時(shí)間和精度的主動(dòng)取舍。做滾動(dòng)調(diào)度時(shí)每輪只有幾分鐘甚至幾十秒窗口與其讓求解器跑到gap0再出結(jié)果不如先拿一個(gè)1%以內(nèi)的可行解把因預(yù)測更新產(chǎn)生的收益吃回來。4.2 時(shí)間限制與MIP start熱啟動(dòng)實(shí)際業(yè)務(wù)中UC很少只跑一次。典型做法是每15分鐘滾動(dòng)一次用最新負(fù)荷預(yù)測重算未來4小時(shí)。如果每次都從零開始求解時(shí)間會(huì)累加得很可觀。這時(shí)候熱啟動(dòng)warm start也稱MIP start就是最大的免費(fèi)午餐。熱啟動(dòng)的本質(zhì)就是把你認(rèn)為“可能比較優(yōu)”的機(jī)組啟停狀態(tài)事先告訴求解器讓分支定界從一開始就貼著這個(gè)解走剪枝效率高很多。在Gurobi里可以這樣傳遞上一輪的最優(yōu)解from gurobipy import Model, GRB m Model(UC) # 這里假設(shè)u_vars已經(jīng)是上一輪求出的0/1矩陣 for i in unit_list: for t in T_list: u_var[i, t].start last_solution[i, t] # start是Gurobi為變量預(yù)留的MIP start屬性 m.optimize()上一輪解出來的機(jī)組狀態(tài)和這一輪差別通常不大熱啟動(dòng)能讓求解時(shí)間縮短30%~70%。但要注意如果負(fù)荷突變特別大比如某臺機(jī)組故障導(dǎo)致負(fù)荷轉(zhuǎn)移熱啟動(dòng)解可能不在可行域內(nèi)求解器會(huì)自動(dòng)剔除不可行部分繼續(xù)從頭探索。所以熱啟動(dòng)是“加速器”不是“保證書”。開源求解器CBC也支持MIP start但pulp里傳入的API在不同版本間不一致使用前要先驗(yàn)證當(dāng)前版本是否支持如果不支持可以放棄熱啟動(dòng)改用“上一輪解作為初始值”手動(dòng)構(gòu)造可行解但意義不大。因此若項(xiàng)目規(guī)模上千臺機(jī)組商業(yè)求解器的熱啟動(dòng)價(jià)值會(huì)非常明顯。4.3 大規(guī)模場景的分解與聚合思路機(jī)組數(shù)量上千時(shí)直接求解MILP的變量數(shù)輕松破百萬時(shí)間和內(nèi)存都受不了。業(yè)界常見做法有三個(gè)機(jī)組聚合、時(shí)間粒度和拉格朗日松弛。機(jī)組聚合是把同類型、同容量區(qū)間的機(jī)組合并成一臺等效機(jī)組。比如某風(fēng)電場或光伏電站有50臺同型號逆變器聚合成一臺整體參與優(yōu)化再按等比例或邊際成本排序分配出力。對于火電同容量機(jī)組也有效。代碼如下def aggregate_units(unit_params, key_fn): groups {} for uid, param in unit_params.items(): key key_fn(param) groups.setdefault(key, []).append(uid) agg {} for key, ids in groups.items(): agg[key] { pmax: sum(unit_params[i][pmax] for i in ids), pmin: sum(unit_params[i][pmin] for i in ids), # 成本系數(shù)按額定容量加權(quán)避免簡單平均帶來的偏差 fc: sum(unit_params[i][fc] * unit_params[i][pmax] for i in ids) / sum(unit_params[i][pmax] for i in ids), } return agg注意加權(quán)平均的目的是讓聚合后的成本特性盡量接近原群體。如果簡單把燃料成本相加聚合機(jī)組成本會(huì)虛高優(yōu)化結(jié)果會(huì)把電量偏向真正的低成本機(jī)組。時(shí)間粒度分解則是先跑每2小時(shí)一個(gè)大時(shí)段的計(jì)劃確定大致啟停方案再把第一個(gè)小時(shí)到第四個(gè)小時(shí)細(xì)化成15分鐘粒度重優(yōu)化。這樣既保留了可執(zhí)行性又降低了問題規(guī)模。拉格朗日松弛把系統(tǒng)功率平衡約束乘子化逐臺機(jī)組獨(dú)立求解再更新乘子迭代工程上常用它來生成初始可行解再交給MILP做最終優(yōu)化。5. 踩坑記錄從zip解壓到約束漏加的五個(gè)典型問題5.1 zip解壓報(bào)invalid zip archive文件沒下完整現(xiàn)象拿到“034電力系統(tǒng)機(jī)組組合優(yōu)化.zip”后雙擊解壓直接彈“文件損壞”或“invalid zip archive: could not find eocd”也有人說成“failed to copy spatial iop zip”這類奇奇怪怪的報(bào)錯(cuò)。原因zip文件尾部有一個(gè)End of Central Directory記錄Linux解壓工具找不到這條記錄就是典型的“文件沒下完整”??赡苁窍螺d中斷、網(wǎng)絡(luò)緩存殘留、或者從網(wǎng)盤拿到的分包沒有全部放在同一目錄就急著解壓。解決先核對文件大小跟下載頁標(biāo)注的字節(jié)數(shù)一致再解壓。第二步用命令行測試unzip -t 034電力系統(tǒng)機(jī)組組合優(yōu)化.zip如果輸出里有“End-of-central-directory signature not found”字樣別猶豫重新下載。若文件特別大且下載工具支持?jǐn)帱c(diǎn)續(xù)傳直接續(xù)傳后再次測試如果是分包壓縮務(wù)必把分卷文件都放在同一個(gè)文件夾否則即使擴(kuò)展名正確也解不出來。還有一個(gè)常見情況普通zip但下載時(shí)被瀏覽器改成“下載.zip.zip”文件大小正常卻解不了這時(shí)改回原名即可。5.2 zip內(nèi)中文亂碼編碼不匹配現(xiàn)象解壓后文件名全是亂碼類似“鎶楃伀”或“緗戠粶”這樣的字符。原因zip包在Windows系統(tǒng)壓縮時(shí)文件名一般按GBK編碼寫入在Linux或macOS上用默認(rèn)UTF-8解壓中文名就變成了亂碼。這跟zip包本身無關(guān)純粹是編碼不一致。解決Linux環(huán)境用unzip -O GBK重新解壓如果用的壓縮軟件不支持指定編碼可以先用Python把zip里的文件名修正出來再根據(jù)修正后的名字手動(dòng)解壓到本地。代碼如下import zipfile with zipfile.ZipFile(034電力系統(tǒng)機(jī)組組合優(yōu)化.zip) as zf: for info in zf.infolist(): # 先把字節(jié)按cp437解碼再按gbk重新編碼還原中文 safe_name info.filename.encode(cp437).decode(gbk, errorsignore) print(safe_name)這個(gè)技巧對含中文的zip非常通用。如果zip本身加了密碼編碼問題反而是小事——沒有密碼只能在腳本里跑字典屬于合規(guī)性以外的話題不做展開。5.3 結(jié)果違反功率平衡單位與約束兩處都查現(xiàn)象求解完成后打印各個(gè)時(shí)段出力之和發(fā)現(xiàn)和負(fù)荷差了幾十甚至幾百M(fèi)W或者某臺機(jī)組狀態(tài)顯示為停機(jī)但出力卻等于Pmin。原因先從約束排查最常見的是漏寫了p[i,t] pmax * u[i,t]和p[i,t] pmin * u[i,t]停機(jī)機(jī)組照樣被分配到出力。其次查數(shù)據(jù)單位負(fù)荷如果用的是MWh而機(jī)組出力是MW一個(gè)24小時(shí)的時(shí)段里差了24倍。解決先加一個(gè)結(jié)果校驗(yàn)?zāi)K每次求解完成后自動(dòng)對賬for t in range(T): gen_sum sum(lp.value(p[i, t]) for i in units) if abs(gen_sum - load[t]) 1e-3: print(f時(shí)段{t1}功率不平衡: 負(fù)荷{load[t]:.2f} 出力{gen_sum:.2f})這段校驗(yàn)不是可選項(xiàng)是必選項(xiàng)。尤其是換數(shù)據(jù)源時(shí)經(jīng)常遇到“MW當(dāng)成MWh”“某機(jī)組Pmax寫錯(cuò)”之類的低級問題有了它對賬翻車概率大大降低。5.4 最小啟停約束導(dǎo)致無解初始狀態(tài)寫錯(cuò)現(xiàn)象模型報(bào)infeasible但把最小運(yùn)行時(shí)間和最小停機(jī)時(shí)間約束全部注釋掉后立刻有解說明問題就出在這組約束上。原因最小啟停約束是按“調(diào)度起點(diǎn)”寫的如果機(jī)組在前一天已經(jīng)運(yùn)行或停機(jī)了若干小時(shí)你必須在模型里把這段“歷史狀態(tài)”接過來。我見過一個(gè)典型案例某機(jī)組最小運(yùn)行時(shí)間8小時(shí)初始狀態(tài)寫成了“在線”但實(shí)際它只連續(xù)運(yùn)行了2小時(shí)模型要求它前6小時(shí)必須保持運(yùn)行可功率平衡又被其他機(jī)組占滿導(dǎo)致無解。根本原因是初始狀態(tài)寫錯(cuò)或根本沒初始化。解決在建模前把每臺機(jī)組的初始運(yùn)行/停機(jī)時(shí)間單獨(dú)建表并在約束中用分段判斷彌補(bǔ)剩余時(shí)長。一個(gè)常見寫法是如果機(jī)組初始在線時(shí)長已超過最小運(yùn)行時(shí)間則首個(gè)時(shí)段不額外施加minup約束如果還沒達(dá)到則讓它在剩余不足的時(shí)段里保持在線。最小停機(jī)時(shí)間同理初始停機(jī)時(shí)長不足時(shí)前幾個(gè)時(shí)段禁止開機(jī)。# init_status: 1在線 0停機(jī) # init_time: 已經(jīng)連續(xù)運(yùn)行或停機(jī)的小時(shí)數(shù) if init_status[i] 1 and init_time[i] units[i][minup]: need units[i][minup] - init_time[i] for t in range(min(need, T)): prob u[i, t] 1這就是前文說的“初始狀態(tài)約束”。很多教程里不寫這層邏輯是因?yàn)榧僭O(shè)機(jī)組已滿足歷史狀態(tài)但真實(shí)調(diào)度場景里凌晨1點(diǎn)出現(xiàn)無解往往就是這個(gè)根因。5.5 求解器報(bào)infeasible先做約束松弛現(xiàn)象模型無解求解器只返回“problem infeasible”不告訴你哪條約束沖突排錯(cuò)全靠猜。原因MILP求解器內(nèi)部用的是分支定界它發(fā)現(xiàn)可行域?yàn)榭諘r(shí)會(huì)直接停止不會(huì)提供一條人話解釋。爬坡約束跨時(shí)段耦合一處無解會(huì)傳導(dǎo)到后續(xù)所有時(shí)段排查特別費(fèi)勁。解決先用“連續(xù)化診斷”找到?jīng)_突范圍。把機(jī)組狀態(tài)變量從二進(jìn)制放松成0到1的連續(xù)變量再求解同一模型。如果連續(xù)模型有解說明是整數(shù)約束和狀態(tài)轉(zhuǎn)移之間的邏輯沖突如果連續(xù)模型也無解說明線性約束之間就矛盾了。連續(xù)模型求解后查看哪些u[i,t]落在0到1之間這些時(shí)段和機(jī)組往往就是約束沖突的高發(fā)區(qū)域for i in units: for t in range(T): u[i, t] lp.LpVariable(fu_{i}_{t}, lowBound0, upBound1) # 由Binary改成Continuous只用于診斷 prob.solve() for i in units: for t in range(T): val lp.value(u[i, t]) if 0.01 val 0.99: print(f診斷: 機(jī)組{i} 時(shí)段{t1} 狀態(tài)值{val:.3f})中間值意味著這個(gè)時(shí)段的狀態(tài)被兩組約束拉扯重點(diǎn)檢查該機(jī)組的爬坡、最小啟停和備用約束。這種“連續(xù)松弛診斷法”在大型模型里比人肉對數(shù)小時(shí)管用得多是我現(xiàn)在排查無解的第一動(dòng)作。6. 進(jìn)階考慮不確定性與回測驗(yàn)證的實(shí)用技巧6.1 從確定性UC到隨機(jī)UC場景法怎么改確定性UC里負(fù)荷是單條曲線隨機(jī)UC常見做法是場景法構(gòu)造多個(gè)負(fù)荷預(yù)測場景每個(gè)場景設(shè)一個(gè)概率。狀態(tài)變量u[i,t]是所有場景共用的非預(yù)期變量而出力p[i,t,s]允許隨場景變化。目標(biāo)函數(shù)改為所有場景成本按概率加權(quán)。代碼層面的改動(dòng)很小只需要給p和load加一個(gè)場景下標(biāo)再把功率平衡約束對每個(gè)場景分別寫。代價(jià)是變量規(guī)模乘以場景數(shù)所以工程上場景常用概率距離方法剪枝到5~20個(gè)。6.2 用歷史負(fù)荷做回測驗(yàn)證結(jié)果是否可信不管模型寫得多復(fù)雜最終要回答“結(jié)果能不能用于明天調(diào)度”。我的習(xí)慣是拿過去30天的歷史負(fù)荷做回測每天用前一天的預(yù)測負(fù)荷跑SCUC再把當(dāng)天實(shí)際負(fù)荷代入結(jié)果統(tǒng)計(jì)功率平衡偏差和爬坡越限次數(shù)。驗(yàn)證腳本核心就一段def validate(plan_result, actual_load, tolerance10.0): imbalance (actual_load - plan_result[total_gen]).abs() violations (imbalance tolerance).sum() return violations # 返回越限時(shí)段數(shù)越少說明方案越穩(wěn)能通過回測的方案才敢排進(jìn)調(diào)度計(jì)劃否則模型再精巧也只是紙面好看。這類優(yōu)化項(xiàng)目里我吃過的最大教訓(xùn)就是“迷信gap和最優(yōu)性”卻忽略了約束本身是否正確。有一版模型把爬坡約束方向?qū)懛辞蠼馄髋艹鐾昝赖?.1% gap拿去現(xiàn)場用才發(fā)現(xiàn)無法執(zhí)行。從那以后任何UC代碼跑完結(jié)果我都強(qiáng)制自己先看啟動(dòng)次數(shù)、單位時(shí)段啟停頻繁度再做平衡校驗(yàn)。如果你手頭正打開這份zip包希望這篇筆記能幫你少走同樣的彎路希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取