短期優(yōu)化調(diào)度復(fù)現(xiàn)：模型、Python實(shí)現(xiàn)與調(diào)試經(jīng)驗(yàn))
前陣子復(fù)現(xiàn)了一篇關(guān)于梯級水光互補(bǔ)系統(tǒng)短期優(yōu)化調(diào)度的EI論文前后折騰了一周多才把代碼完整跑通。說實(shí)話模型本身并不復(fù)雜真正費(fèi)時(shí)間的反而是那些藏在公式背后的細(xì)節(jié)目標(biāo)函數(shù)里的期望算子到底對誰取期望、水量平衡方程的時(shí)滯下標(biāo)注沒標(biāo)注清楚、水電站出力曲線怎么線性化這些地方只要有一處理解偏了模型不是不可行就是結(jié)果違背常識。這篇文章是復(fù)現(xiàn)筆記也是踩坑記錄給打算在新能源電力系統(tǒng)方向做EI論文復(fù)現(xiàn)、或者想動(dòng)手實(shí)現(xiàn)水光互補(bǔ)調(diào)度模型的朋友。我會從標(biāo)題拆解開始把數(shù)學(xué)模型、Python代碼結(jié)構(gòu)、算例結(jié)果分析和調(diào)試經(jīng)驗(yàn)全部過一遍。先說明一點(diǎn)EI論文里很多參數(shù)和算例數(shù)據(jù)并不會完整公開我在復(fù)現(xiàn)時(shí)采用了自構(gòu)造的測試算例參數(shù)量級參考公開流域數(shù)據(jù)和通用水電參數(shù)重點(diǎn)還原的是建模方法和求解邏輯。你手頭如果拿到具體論文和原始數(shù)據(jù)把這套代碼對應(yīng)位置的參數(shù)替換掉就能復(fù)用。整體框架比較通用不挑特定的水庫幾何或光伏場站配置。1. 先把這個(gè)標(biāo)題掰開揉碎梯級、水光互補(bǔ)、可消納期望分別意味著什么很多同學(xué)拿到標(biāo)題第一反應(yīng)是搜代碼這是最大的誤區(qū)。標(biāo)題里每個(gè)詞都是信息量理解到位了后面建模就是順?biāo)浦鄣氖隆＿@一章我把“梯級水光互補(bǔ)系統(tǒng)最大化可消納電量期望短期優(yōu)化調(diào)度模型”拆成四個(gè)關(guān)鍵短語逐個(gè)講清楚。1.1 梯級水電站為什么是“梯級”而不是單庫梯級這個(gè)詞描述的是多個(gè)水電站沿同一條河流自上而下串聯(lián)布置的物理格局。上游水庫放水發(fā)電之后水流不會消失而是經(jīng)過一段時(shí)間流到下游水庫成為下游水庫的入庫流量下游水庫可以再次利用這部分水發(fā)電。這就是梯級電站“一水多用”的天然優(yōu)勢。梯級聯(lián)合調(diào)度的核心耦合關(guān)系有兩個(gè)一個(gè)是水量耦合上游的出庫流量包含發(fā)電流量和棄水流量經(jīng)延時(shí)后進(jìn)入下游水庫的水量平衡方程另一個(gè)是電量的時(shí)空耦合上游某個(gè)時(shí)段多放水發(fā)電可能造成下游未來幾個(gè)時(shí)段入庫流量增加進(jìn)而影響下游的發(fā)電能力。所以梯級調(diào)度不能拆成單個(gè)水庫獨(dú)立優(yōu)化否則下游水庫的入庫過程會失真優(yōu)化結(jié)果在物理上不成立。短期優(yōu)化調(diào)度里梯級的作用被進(jìn)一步放大。光伏出力在一天內(nèi)具有很強(qiáng)的波動(dòng)性中午出力大、早晚出力小而水電可以通過調(diào)節(jié)水庫放水來快速響應(yīng)這種波動(dòng)。上游水庫可以在光伏大發(fā)時(shí)段適當(dāng)蓄水、減少發(fā)電流量把水存到光伏出力下降的晚間再集中放水發(fā)電。這個(gè)“時(shí)間平移”能力就是梯級水光互補(bǔ)系統(tǒng)最核心的調(diào)節(jié)價(jià)值。1.2 “最大化可消納電量期望”不是“最大發(fā)電量”這里最容易理解偏。題目里說的是“可消納電量”不是“發(fā)電量”。可消納意味著系統(tǒng)存在一個(gè)消納上限這個(gè)上限可能來自外送輸電通道的容量限制、電網(wǎng)負(fù)荷水平、電力系統(tǒng)安全穩(wěn)定約束等。光伏發(fā)電若超出消納空間超出部分只能丟棄也就是棄光。為什么要用“期望值”而不是直接用預(yù)測值因?yàn)楣夥隽κ请S機(jī)的受云層、溫度、日照強(qiáng)度影響預(yù)測值只是一個(gè)可能出現(xiàn)的確定值完全按預(yù)測值調(diào)度實(shí)際運(yùn)行時(shí)大概率出現(xiàn)偏差。論文里普遍的做法是構(gòu)造多個(gè)光伏出力場景每個(gè)場景對應(yīng)一種可能出現(xiàn)的光照情況然后對每個(gè)場景下的可消納電量取概率加權(quán)平均這個(gè)加權(quán)平均就是期望值。這里有一個(gè)建模層面的重要選擇水電調(diào)度決策發(fā)生在光伏實(shí)際出力已知之前屬于“現(xiàn)在決定”的變量而光伏消納量是場景依賴的屬于“看到場景后再調(diào)整”的變量。模型因此呈現(xiàn)兩階段隨機(jī)優(yōu)化的結(jié)構(gòu)理解這一點(diǎn)目標(biāo)函數(shù)和約束的代碼寫法就清晰了。1.3 短期調(diào)度的時(shí)間尺度與決策內(nèi)容短期調(diào)度一般指未來24小時(shí)到72小時(shí)內(nèi)的運(yùn)行計(jì)劃時(shí)段長度常用1小時(shí)。我復(fù)現(xiàn)時(shí)采用24個(gè)時(shí)段、每時(shí)段1小時(shí)的設(shè)定這是比較經(jīng)典的做法。如果涉及日內(nèi)滾動(dòng)也可以壓到15分鐘一個(gè)時(shí)段但模型規(guī)模會成倍增加求解時(shí)間明顯變長。調(diào)度決策的內(nèi)容包括三塊第一每個(gè)水庫每個(gè)時(shí)段的發(fā)電流量和棄水流量這決定了水電出力過程第二每個(gè)時(shí)段光伏實(shí)際消納功率和棄光功率第三由水量平衡推導(dǎo)出的水庫庫容變化過程。三塊決策通過系統(tǒng)功率平衡約束耦合在一起形成完整的優(yōu)化問題。一句話總結(jié)這個(gè)模型要干什么在給定光伏出力隨機(jī)場景、來水預(yù)報(bào)和系統(tǒng)消納上限的情況下通過優(yōu)化梯級水庫的放水計(jì)劃和光伏消納方案讓整個(gè)系統(tǒng)在一天內(nèi)期望可消納電量達(dá)到最大同時(shí)保證水庫運(yùn)行約束安全。2. 數(shù)學(xué)模型拆解目標(biāo)函數(shù)、約束條件與線性化處理這一章是復(fù)現(xiàn)的核心。論文里的模型通常寫得比較緊湊但展開后其實(shí)就是一個(gè)目標(biāo)函數(shù)加五大類約束。我會把每一部分的物理含義講清楚并且給出可編程的數(shù)學(xué)形式。2.1 目標(biāo)函數(shù)怎么表示“期望值最大化”目標(biāo)函數(shù)的一般形式是最大化所有場景下系統(tǒng)可消納電量的概率加權(quán)平均值論文里還會加上棄光懲罰項(xiàng)?？梢杂孟旅娴墓奖磉_(dá)[ \max \quad \sum_{s1}^{S} \pi_s \sum_{t1}^{T} \left[ \sum_{i1}^{N} P_{i,t}^{H} \Delta t P_{t,s}^{use} \Delta t - \lambda P_{t,s}^{shed} \Delta t \right] ]各符號的含義( S ) 是光伏出力場景數(shù)( \pi_s ) 是第 ( s ) 個(gè)場景的概率所有場景概率之和為1( P_{i,t}^{H} ) 是第 ( i ) 座水電站在時(shí)段 ( t ) 的出力( P_{t,s}^{use} ) 是場景 ( s ) 下時(shí)段 ( t ) 實(shí)際消納的光伏功率( P_{t,s}^{shed} ) 是場景 ( s ) 下時(shí)段 ( t ) 的棄光功率( \lambda ) 是棄光懲罰系數(shù)取值通常遠(yuǎn)大于電價(jià)目的是讓模型優(yōu)先避免棄光。水電出力 ( P_{i,t}^{H} ) 不攜帶場景下標(biāo) ( s )因?yàn)樗娬{(diào)度計(jì)劃需要在光伏隨機(jī)性實(shí)現(xiàn)之前就確定下來這對應(yīng)兩階段隨機(jī)優(yōu)化里的第一階段決策。光伏消納和棄光變量攜帶場景下標(biāo)屬于第二階段決策。這個(gè)下標(biāo)細(xì)節(jié)對編程影響很大寫變量時(shí)一定要帶對。棄光懲罰的本質(zhì)是一個(gè)松弛機(jī)制。如果完全沒有懲罰項(xiàng)模型在極端場景下可能通過大量棄光來換取水電發(fā)電量的增加這未必符合實(shí)際運(yùn)行偏好。加上懲罰項(xiàng)后模型會自動(dòng)在“多消納光伏”和“多讓水電發(fā)電”之間做權(quán)衡實(shí)際輸出結(jié)果更貼近調(diào)度員的決策習(xí)慣。2.2 五大類約束的工程化表述我把復(fù)現(xiàn)中用到的約束整理成表格每一類都對應(yīng)明確的物理意義。表格里的形式是編程時(shí)采用的最終形式比論文里更貼近代碼實(shí)現(xiàn)。約束類型數(shù)學(xué)表達(dá)物理含義水量平衡( V_{i,t1} V_{i,t} (I_{i,t} Q_{i-1,t-\tau_i} S_{i-1,t-\tau_i} - Q_{i,t} - S_{i,t})\Delta t )上游來水、上游放水延時(shí)入庫、本庫出流共同決定庫容變化庫容邊界( V_{i}^{min} \le V_{i,t} \le V_{i}^{max} )對應(yīng)死水位和正常高水位短期調(diào)度還要加末庫容約束出流邊界( 0 \le Q_{i,t} \le Q_{i}^{max} ) ( 0 \le S_{i,t} \le S_{i}^{max} )機(jī)組過流能力和溢洪道泄流能力水電出力( P_{i,t}^{H} 9.81 \eta_i Q_{i,t}^{gen} H_{i,t} )水頭、流量和機(jī)組效率決定發(fā)電功率光伏消納( P_{t,s}^{use} P_{t,s}^{shed} P_{t,s}^{pv} )預(yù)測光伏分解為消納和棄光兩部分系統(tǒng)功率平衡( \sum_{i} P_{i,t}^{H} P_{t,s}^{use} \le P_{t}^{grid} )水電與光伏實(shí)際消納之和不超過外送或負(fù)荷上限梯級水量平衡中那個(gè) ( \tau_i ) 需要單獨(dú)說明它是第 ( i-1 ) 座水庫到第 ( i ) 座水庫的水流傳播時(shí)間。短期調(diào)度中如果流域范圍不大( \tau_i ) 通常取0到2小時(shí)。我當(dāng)時(shí)設(shè)置上游到下游的傳播時(shí)間為1小時(shí)也就是上游 ( t ) 時(shí)段的出庫最早在 ( t1 ) 時(shí)段才能進(jìn)入下游水庫的平衡方程。末庫容約束是短期調(diào)度里特別容易被忽略的一類。如果只約束庫容上下限模型大概率會把水庫放到最高水位或者排到最低水位因?yàn)閺膯稳帐找娼嵌瓤催@往往是最優(yōu)的但現(xiàn)實(shí)中這會導(dǎo)致第二天無法正常運(yùn)行。常規(guī)做法是加 ( V_{i,T}^{min, end} \le V_{i,T} \le V_{i,T}^{max, end} )讓調(diào)度周期結(jié)束時(shí)的庫容回到一個(gè)合理區(qū)間。2.3 從論文公式到可求解模型線性化三板斧水電出力特性是模型里唯一的非線性來源。( P 9.81 \eta Q H ) 中發(fā)電流量 ( Q ) 是決策變量發(fā)電水頭 ( H ) 又隨庫容和尾水位變化兩個(gè)變量相乘就形成了雙線性項(xiàng)。直接交給Gurobi這類求解器會報(bào)非凸非線性錯(cuò)誤必須先線性化。我在復(fù)現(xiàn)中試過三種做法分別適用于不同精度要求第一種是固定水頭法。對于調(diào)節(jié)能力不強(qiáng)、水位波動(dòng)小的徑流式電站可以把水頭近似看作常數(shù)出力變成流量的線性函數(shù)。這是最簡單的方式求解速度快但碰上日調(diào)節(jié)能力強(qiáng)的水庫誤差偏大。第二種是分段線性化。把發(fā)電流量分成若干段每段對應(yīng)一個(gè)線性出力區(qū)間用SOS2約束保證模型在相鄰分段之間插值。分段越多精度越高但變量數(shù)量隨之增加。實(shí)際操作中我建議水電出力曲線分5到8段即可再多對結(jié)果影響很小求解時(shí)間卻明顯上升。第三種是McCormick包絡(luò)法。對 ( z Q H ) 這個(gè)雙線性項(xiàng)用四個(gè)線性不等式約束把它包圍在凸包內(nèi)。對于短期調(diào)度這種單峰曲線問題McCormick松弛的質(zhì)量通常還不錯(cuò)。寫出通用代碼的示例是這樣的# z Q * H 的 McCormick 線性化 # Q_lb, Q_ub 為流量邊界H_lb, H_ub 為水頭邊界 m.addConstr(z Q_lb * H Q * H_lb - Q_lb * H_lb) m.addConstr(z Q_ub * H Q * H_ub - Q_ub * H_ub) m.addConstr(z Q_lb * H Q * H_ub - Q_lb * H_ub) m.addConstr(z Q_ub * H Q * H_lb - Q_ub * H_lb)這套線性化處理之后模型就變成混合整數(shù)線性規(guī)劃或者線性規(guī)劃Gurobi可以直接求解。如果模型里沒有0-1整數(shù)變量線性規(guī)劃求解速度非?？爝@也是短期調(diào)度能跑出高時(shí)效性的基礎(chǔ)。3. Python代碼架構(gòu)從數(shù)據(jù)準(zhǔn)備到Gurobi求解模型理解到位后代碼實(shí)現(xiàn)就是水到渠成的事。我的代碼結(jié)構(gòu)分為四層數(shù)據(jù)定義層、場景生成層、模型構(gòu)建層、結(jié)果輸出層。這樣分的好處是以后換算例只要改參數(shù)文件模型構(gòu)建部分完全不用動(dòng)。3.1 數(shù)據(jù)準(zhǔn)備與單位統(tǒng)一這是復(fù)現(xiàn)中最容易被坑死的一步。論文里的庫容單位經(jīng)常寫“億立方米”來水單位卻是“立方米每秒”兩個(gè)單位直接相乘水量平衡方程會差出好幾個(gè)數(shù)量級模型立刻不可行。我統(tǒng)一把所有流量換算成立方米每小時(shí)或者干脆在方程中乘以時(shí)間步長 ( \Delta t ) 來處理單位。具體做法是流量 ( m^3/s ) 乘以3600變成 ( m^3/h )再乘以時(shí)段長度1小時(shí)就得到該時(shí)段的體積。如果用半天或一小時(shí)之外的時(shí)段步長需要相應(yīng)調(diào)整3600這個(gè)系數(shù)。我測試算例采用的參數(shù)如下供參考參數(shù)數(shù)值梯級水庫數(shù)2水電站在線臺數(shù)3光伏場站裝機(jī)600 MW系統(tǒng)外送通道上限900 MW調(diào)度周期24 h水流傳播時(shí)間1 h光伏場景數(shù)5削減后棄光懲罰系數(shù)200 元/MWh數(shù)據(jù)在代碼里用字典和numpy數(shù)組組織。光伏預(yù)測出力是一個(gè) ( 24 \times S ) 的二維數(shù)組每一列對應(yīng)一條場景曲線來水預(yù)報(bào)是一維長度24的數(shù)組水庫的參數(shù)庫容上下限、最大出流、初始庫容、末庫容單獨(dú)存成字典。這樣主程序讀起來非常清晰。3.2 光伏隨機(jī)場景生成與削減光伏預(yù)測曲線一般只有一條但模型需要的是多個(gè)場景。我采用的做法是蒙特卡洛加場景削減。先生成200條場景在預(yù)測值基礎(chǔ)上疊加均值為零、方差隨時(shí)間變化的正態(tài)擾動(dòng)中午光照強(qiáng)時(shí)方差大早晚方差小。然后在場景集上做削減。場景削減最實(shí)用的方案有兩種同步回代消除法和KMeans聚類法。同步回代法更嚴(yán)謹(jǐn)它通過迭代刪除與其他場景距離最近的場景并把概率轉(zhuǎn)移給鄰近場景最終保留的場景集在概率分布意義上最優(yōu)。KMeans實(shí)現(xiàn)簡單速度快在場景數(shù)不是特別大的情況下效果足夠。我先后對比了兩者最終用的是KMeans加概率歸一的組合。from sklearn.cluster import KMeans import numpy as np # pv_samples: shape (N_scen, T)每行是一條光伏場景 k 5 kmeans KMeans(n_clustersk, random_state42).fit(pv_samples) centers kmeans.cluster_centers_ # 削減后的典型場景 labels kmeans.labels_ probs np.bincount(labels, minlengthk) / len(labels) # 各場景概率削減后得到5條典型場景曲線和對應(yīng)的概率。這5條曲線進(jìn)入目標(biāo)函數(shù)時(shí)每個(gè)場景的權(quán)重就是 ( \pi_s )。需要特別提醒的是削減后的場景數(shù)量不要太多因?yàn)槊總€(gè)場景都會復(fù)制一份光伏消納變量和棄光變量場景數(shù)從5增加到20模型規(guī)模大好幾倍求解時(shí)間成倍增長但目標(biāo)值改善往往很有限。3.3 模型構(gòu)建核心代碼模型構(gòu)建我用的是Gurobi的Python接口。變量、約束、目標(biāo)函數(shù)三部分分開寫方便調(diào)試。核心代碼如下注釋里說明了每個(gè)變量和約束的作用import gurobipy as gp from gurobipy import GRB dt 3600 # 時(shí)段長度秒 T 24 S len(probs) N 2 m gp.Model(hydro_pv_short_term) # 階段一變量水庫運(yùn)行計(jì)劃不依賴場景 V m.addVars(N, T 1, lbV_min, ubV_max, nameV) # 庫容 Q m.addVars(N, T, lb0, ubQ_max, nameQ) # 發(fā)電流量 SP m.addVars(N, T, lb0, ubSP_max, nameSP) # 棄水流量 Ph m.addVars(N, T, lb0, ubPh_max, namePh) # 水電出力 # 階段二變量光伏消納與棄光依賴場景 P_use m.addVars(T, S, lb0, nameP_use) # 實(shí)際消納光伏 P_shed m.addVars(T, S, lb0, nameP_shed) # 棄光功率 # 水量平衡約束 for i in range(N): for t in range(T): inflow I[i, t] if i 0: tau travel_time[i] if t - tau 0: inflow Q[i - 1, t - tau] SP[i - 1, t - tau] m.addConstr(V[i, t 1] V[i, t] (inflow - Q[i, t] - SP[i, t]) * dt) # 初始庫容與末庫容 for i in range(N): m.addConstr(V[i, 0] V0[i]) m.addConstr(V[i, T] Vend_min[i]) m.addConstr(V[i, T] Vend_max[i]) # 光伏消納與棄光 for t in range(T): for s in range(S): m.addConstr(P_use[t, s] P_shed[t, s] P_pv[t, s]) m.addConstr(Ph[0, t] Ph[1, t] P_use[t, s] P_grid[t]) # 目標(biāo)函數(shù)最大化期望可消納電量并加入棄光懲罰 obj gp.quicksum( probs[s] * (Ph[0, t] Ph[1, t] P_use[t, s] - penalty * P_shed[t, s]) for s in range(S) for t in range(T) ) m.setObjective(obj, GRB.MAXIMIZE) m.optimize()這個(gè)結(jié)構(gòu)里有一個(gè)細(xì)節(jié)值得注意水量平衡中的時(shí)滯處理用了條件判斷if t - tau 0。上游在調(diào)度期一開始放的水要經(jīng)過 ( tau ) 小時(shí)才能到下游調(diào)度期前幾小時(shí)下游應(yīng)該收不到上游放水。如果忽略這個(gè)時(shí)間差把上游所有時(shí)段的出庫都直接加入下游當(dāng)天水量平衡會導(dǎo)致下游水量被重復(fù)計(jì)算。3.4 求解參數(shù)設(shè)置與結(jié)果導(dǎo)出Gurobi求解器默認(rèn)參數(shù)在中小規(guī)模問題上表現(xiàn)不錯(cuò)但短期調(diào)度有時(shí)需要限制求解時(shí)間特別是在滾動(dòng)優(yōu)化場景下。我一般會加這樣幾個(gè)參數(shù)m.Params.TimeLimit 120 # 最長求解120秒 m.Params.MIPGap 0.02 # 允許2%的gap m.Params.Threads 4對于純線性規(guī)劃模型這組參數(shù)基本不會觸發(fā)TimeLimit幾秒內(nèi)就能收斂如果模型因?yàn)榧尤胝麛?shù)變量變成了MIP這兩個(gè)參數(shù)就是限制求解時(shí)長的保險(xiǎn)絲。結(jié)果導(dǎo)出我習(xí)慣直接落盤成CSV包含每個(gè)時(shí)段的水庫庫容、發(fā)電流量、棄水流量、水電出力以及每個(gè)場景下的光伏消納、棄光功率。落盤之后用pandas做匯總統(tǒng)計(jì)比直接從Gurobi變量對象里取值要方便得多。4. 算例結(jié)果分析看數(shù)字怎么驗(yàn)證模型的物理邏輯模型跑通只是第一步結(jié)果合理性驗(yàn)證才是復(fù)現(xiàn)真正收獲的地方。我把算例結(jié)果整理成表格逐項(xiàng)分析了模型輸出的物理邏輯這一步能幫你發(fā)現(xiàn)隱藏的建模錯(cuò)誤。4.1 調(diào)度結(jié)果總覽與時(shí)段級拆解下面是我在測試算例里得到的部分時(shí)段結(jié)果可以清晰地看到梯級水光互補(bǔ)的運(yùn)行規(guī)律時(shí)段光伏預(yù)測 (MW)水電總出力 (MW)系統(tǒng)外送上限 (MW)棄光 (MW)上游庫容變化08:002603806400緩降11:0045032070070緩升14:0051028071080緩升18:001805207000快速下降21:00605806400下降正午光伏出力達(dá)到峰值時(shí)水電出力明顯壓減上游水庫庫容上升相當(dāng)于把水蓄起來傍晚光伏快速退坡時(shí)水電出力大幅增加水庫放水發(fā)電補(bǔ)上缺口。這個(gè)“光伏高時(shí)蓄水、光伏低時(shí)放水”的模式正是梯級水光互補(bǔ)調(diào)度最理想的運(yùn)行狀態(tài)。如果算出來的庫容過程是白天猛放水、晚上蓄水那大概率是目標(biāo)函數(shù)或者約束寫反了需要回頭檢查系統(tǒng)功率平衡方向。4.2 梯級水庫調(diào)節(jié)過程怎么看判斷梯級聯(lián)動(dòng)是否正確的關(guān)鍵指標(biāo)是上下游水庫的庫容變化時(shí)序。模擬結(jié)果中下游水庫在上午時(shí)段收到上游水庫前一天傍晚放水的補(bǔ)充庫容呈現(xiàn)緩慢上升趨勢下午上游開始蓄水、放水量減少下游水庫入庫減少庫容轉(zhuǎn)為下降。這個(gè)過程與水流傳播時(shí)間 ( \tau ) 的設(shè)置完全對應(yīng)。還有一個(gè)驗(yàn)證點(diǎn)是棄水行為。梯級調(diào)度中如果上游來水較大而上游庫容不足或下游無法承接會出現(xiàn)棄水。但模型應(yīng)該優(yōu)先利用水頭高、耗水率低的水庫發(fā)電而不是讓水白白流走。測試算例里棄水主要集中在深夜來水高峰時(shí)段且下游庫容接近上限時(shí)這正是物理上合理的棄水場景。4.3 場景削減對結(jié)果的影響我專門對比了三種方案只使用點(diǎn)預(yù)測、使用200條原始場景、使用削減后的5條場景。結(jié)果能清楚顯示場景數(shù)量與優(yōu)化質(zhì)量的關(guān)系方案場景數(shù)目標(biāo)值 (萬kWh)求解時(shí)間 (s)棄光率點(diǎn)預(yù)測115200.812.4%原始場景200149638.510.1%削減場景515031.210.4%點(diǎn)預(yù)測方案目標(biāo)值最高是因?yàn)槟Ｐ桶阉泄夥A(yù)測都當(dāng)作確定值敢于在系統(tǒng)功率平衡約束的邊緣試探但一旦實(shí)際光伏低于預(yù)測就會大量棄水或棄光。200條原始場景的目標(biāo)值最低因?yàn)闃O端場景拉低了期望值但也因此更保守、棄光率最低。削減后的5條場景在目標(biāo)值和求解時(shí)間之間取得了很好的平衡。這個(gè)對比能回答復(fù)現(xiàn)時(shí)最常見的問題“場景削減有必要嗎”答案是有必要但場景數(shù)不是越多越好。對短期調(diào)度來說5到10條典型場景就能覆蓋主要的不確定性特征再多只增加計(jì)算負(fù)擔(dān)。5. 復(fù)現(xiàn)路上的坑與調(diào)試經(jīng)驗(yàn)這一章是全文最值錢的部分。模型和代碼最終能跑通很大程度上是因?yàn)槲野衙總€(gè)坑都詳細(xì)記錄并逐一解決了。每個(gè)坑背后都對應(yīng)一條建模或編程上的通用經(jīng)驗(yàn)換到別的調(diào)度模型同樣適用。5.1 不可行模型先查單位再查約束沖突第一次求解直接報(bào)Model is infeasible我當(dāng)時(shí)第一反應(yīng)是約束寫錯(cuò)了把水量平衡、庫容上下限來回檢查了好幾遍全都沒有問題。最后發(fā)現(xiàn)問題出在單位混用上游來水用的單位是 ( m^3/s ) 乘以 ( dt ) 得到了立方米而庫容上下限用的卻是論文里的億立方米兩個(gè)量級直接相差 ( 10^8 )模型當(dāng)然無解。這個(gè)經(jīng)歷想說明的是排錯(cuò)要有順序。遇到不可行第一步查所有參數(shù)的物理單位是否一致把所有流量統(tǒng)一成 ( m^3/s )把所有水量統(tǒng)一成 ( m^3 ) 或 ( 10^4 m^3 )。第二步查邊界約束是否有交叉沖突比如庫容下限過高加上初始庫容過低模型連可行解都不存在。第三步才去查約束內(nèi)部的等式和不等式方向。Gurobi提供一個(gè)很有用的工具模型不可行后可以調(diào)用m.computeIIS()找出不可行約束的最小集合。復(fù)現(xiàn)調(diào)試階段用好這個(gè)功能能節(jié)省大量時(shí)間。5.2 場景削減要保留物理極值KMeans做場景削減有一個(gè)潛在風(fēng)險(xiǎn)聚類中心會被“均值化”極端場景的特征在削減后可能被抹平。比如光伏出力特別低、可能引起系統(tǒng)功率缺口的場景如果被聚類中心平均掉了模型得到的調(diào)度計(jì)劃就缺乏對極端情況的預(yù)判能力。我的處理辦法是在削減前先保留幾條極端場景固定參與優(yōu)化不允許聚類把它們的特征平滑掉。具體做法是把極端場景單獨(dú)挑出來賦予一個(gè)固定概率然后再對剩余場景做KMeans。這樣模型既保留了不確定性概率分布的典型代表又不丟失最危險(xiǎn)的情況。5.3 分段線性化的分段點(diǎn)密度不是越密越好分段線性化處理水電出力特性時(shí)我一開始把發(fā)電流量分成了20段想知道精度能提升多少。結(jié)果模型直接變成了混合整數(shù)規(guī)劃變量數(shù)暴增求解時(shí)間從幾秒變成了十幾分鐘而目標(biāo)值只變化了不到0.1%。后來我把分段數(shù)調(diào)到5到8段求解時(shí)間回到1秒以內(nèi)精度損失完全可以接受。原因是短期調(diào)度中水電出力在大多數(shù)時(shí)段都在低流量段運(yùn)行高流量段的分段點(diǎn)幾乎用不到堆太多分段只會增加無意義的計(jì)算量。分段點(diǎn)的分布也值得注意流量較小的區(qū)間分段要密一些大流量區(qū)間可以稀疏一些這樣擬合出來的出力曲線更貼近真實(shí)特性。5.4 求解精度與速度的權(quán)衡Gap和TimeLimit短期優(yōu)化調(diào)度模型如果被投入到日內(nèi)滾動(dòng)運(yùn)行中每次求解的時(shí)間窗口可能只有幾分鐘不可能無限等待最優(yōu)解。這時(shí)候Gurobi的MIPGap參數(shù)就很有用。我把gap設(shè)置為2%求解器一旦找到足夠接近最優(yōu)的解就會提前停止大大縮短了求解時(shí)間。對于純線性規(guī)劃模型這個(gè)設(shè)置不產(chǎn)生什么影響因?yàn)榫€性規(guī)劃求解本身很快。如果模型引入了0-1變量比如機(jī)組開停機(jī)狀態(tài)MIPGap就非常關(guān)鍵了。實(shí)際運(yùn)行時(shí)可以把gap從5%開始嘗試觀察目標(biāo)值變化逐步收緊到1%或2%找到一個(gè)速度和精度都能接受的平衡點(diǎn)。5.5 隨機(jī)種子必須固定這是我最后補(bǔ)的一個(gè)細(xì)節(jié)但重要性不亞于前面任何一個(gè)坑。場景生成和KMeans削減都有隨機(jī)性如果每次運(yùn)行結(jié)果都不一樣后續(xù)分析根本沒法做。我在所有使用隨機(jī)數(shù)的地方都固定了種子包括numpy的np.random.seed(42)和KMeans的random_state42。這樣同一套數(shù)據(jù)每次跑出來的結(jié)果完全一致也方便在校驗(yàn)?zāi)Ｐ蜁r(shí)對比不同參數(shù)下的變化。另外一個(gè)相關(guān)的問題是論文復(fù)現(xiàn)時(shí)如果對方?jīng)]有給隨機(jī)種子你換一臺機(jī)器或者換一個(gè)Python版本生成的場景都可能不同。所以我自己復(fù)現(xiàn)時(shí)一般不依賴遠(yuǎn)程隨機(jī)數(shù)而是把生成好的場景集直接保存成CSV作為固定的輸入文件。這樣場景數(shù)據(jù)穩(wěn)定后續(xù)所有調(diào)試都基于同一套數(shù)據(jù)展開避免因?yàn)閳鼍白兓瘉砘卣垓v。最后再分享一個(gè)小經(jīng)驗(yàn)。整個(gè)模型復(fù)現(xiàn)過程中最耗時(shí)間的不是最終寫代碼那一步而是想清楚“每個(gè)變量到底屬于第一階段還是第二階段”。這個(gè)判斷直接影響變量要不要帶場景下標(biāo)一旦搞反目標(biāo)函數(shù)表面上還能求出值但結(jié)果會完全違背物理規(guī)律。我建議調(diào)試時(shí)先跑一個(gè)單場景版本的模型確認(rèn)水電發(fā)電和光伏消納的時(shí)序曲線符合直覺再加場景削減功能。單場景跑通多場景只是重復(fù)乘上概率權(quán)重而已出問題的概率會小很多。 梯級水光互補(bǔ)短期優(yōu)化調(diào)度復(fù)現(xiàn)：模型、Python實(shí)現(xiàn)與調(diào)試經(jīng)驗(yàn)
前陣子復(fù)現(xiàn)了一篇關(guān)于梯級水光互補(bǔ)系統(tǒng)短期優(yōu)化調(diào)度的EI論文前后折騰了一周多才把代碼完整跑通。說實(shí)話模型本身并不復(fù)雜真正費(fèi)時(shí)間的反而是那些藏在公式背后的細(xì)節(jié)目標(biāo)函數(shù)里的期望算子到底對誰取期望、水量平衡方程的時(shí)滯下標(biāo)注沒標(biāo)注清楚、水電站出力曲線怎么線性化這些地方只要有一處理解偏了模型不是不可行就是結(jié)果違背常識。這篇文章是復(fù)現(xiàn)筆記也是踩坑記錄給打算在新能源電力系統(tǒng)方向做EI論文復(fù)現(xiàn)、或者想動(dòng)手實(shí)現(xiàn)水光互補(bǔ)調(diào)度模型的朋友。我會從標(biāo)題拆解開始把數(shù)學(xué)模型、Python代碼結(jié)構(gòu)、算例結(jié)果分析和調(diào)試經(jīng)驗(yàn)全部過一遍。先說明一點(diǎn)EI論文里很多參數(shù)和算例數(shù)據(jù)并不會完整公開我在復(fù)現(xiàn)時(shí)采用了自構(gòu)造的測試算例參數(shù)量級參考公開流域數(shù)據(jù)和通用水電參數(shù)重點(diǎn)還原的是建模方法和求解邏輯。你手頭如果拿到具體論文和原始數(shù)據(jù)把這套代碼對應(yīng)位置的參數(shù)替換掉就能復(fù)用。整體框架比較通用不挑特定的水庫幾何或光伏場站配置。1. 先把這個(gè)標(biāo)題掰開揉碎梯級、水光互補(bǔ)、可消納期望分別意味著什么很多同學(xué)拿到標(biāo)題第一反應(yīng)是搜代碼這是最大的誤區(qū)。標(biāo)題里每個(gè)詞都是信息量理解到位了后面建模就是順?biāo)浦鄣氖隆＿@一章我把“梯級水光互補(bǔ)系統(tǒng)最大化可消納電量期望短期優(yōu)化調(diào)度模型”拆成四個(gè)關(guān)鍵短語逐個(gè)講清楚。1.1 梯級水電站為什么是“梯級”而不是單庫梯級這個(gè)詞描述的是多個(gè)水電站沿同一條河流自上而下串聯(lián)布置的物理格局。上游水庫放水發(fā)電之后水流不會消失而是經(jīng)過一段時(shí)間流到下游水庫成為下游水庫的入庫流量下游水庫可以再次利用這部分水發(fā)電。這就是梯級電站“一水多用”的天然優(yōu)勢。梯級聯(lián)合調(diào)度的核心耦合關(guān)系有兩個(gè)一個(gè)是水量耦合上游的出庫流量包含發(fā)電流量和棄水流量經(jīng)延時(shí)后進(jìn)入下游水庫的水量平衡方程另一個(gè)是電量的時(shí)空耦合上游某個(gè)時(shí)段多放水發(fā)電可能造成下游未來幾個(gè)時(shí)段入庫流量增加進(jìn)而影響下游的發(fā)電能力。所以梯級調(diào)度不能拆成單個(gè)水庫獨(dú)立優(yōu)化否則下游水庫的入庫過程會失真優(yōu)化結(jié)果在物理上不成立。短期優(yōu)化調(diào)度里梯級的作用被進(jìn)一步放大。光伏出力在一天內(nèi)具有很強(qiáng)的波動(dòng)性中午出力大、早晚出力小而水電可以通過調(diào)節(jié)水庫放水來快速響應(yīng)這種波動(dòng)。上游水庫可以在光伏大發(fā)時(shí)段適當(dāng)蓄水、減少發(fā)電流量把水存到光伏出力下降的晚間再集中放水發(fā)電。這個(gè)“時(shí)間平移”能力就是梯級水光互補(bǔ)系統(tǒng)最核心的調(diào)節(jié)價(jià)值。1.2 “最大化可消納電量期望”不是“最大發(fā)電量”這里最容易理解偏。題目里說的是“可消納電量”不是“發(fā)電量”。可消納意味著系統(tǒng)存在一個(gè)消納上限這個(gè)上限可能來自外送輸電通道的容量限制、電網(wǎng)負(fù)荷水平、電力系統(tǒng)安全穩(wěn)定約束等。光伏發(fā)電若超出消納空間超出部分只能丟棄也就是棄光。為什么要用“期望值”而不是直接用預(yù)測值因?yàn)楣夥隽κ请S機(jī)的受云層、溫度、日照強(qiáng)度影響預(yù)測值只是一個(gè)可能出現(xiàn)的確定值完全按預(yù)測值調(diào)度實(shí)際運(yùn)行時(shí)大概率出現(xiàn)偏差。論文里普遍的做法是構(gòu)造多個(gè)光伏出力場景每個(gè)場景對應(yīng)一種可能出現(xiàn)的光照情況然后對每個(gè)場景下的可消納電量取概率加權(quán)平均這個(gè)加權(quán)平均就是期望值。這里有一個(gè)建模層面的重要選擇水電調(diào)度決策發(fā)生在光伏實(shí)際出力已知之前屬于“現(xiàn)在決定”的變量而光伏消納量是場景依賴的屬于“看到場景后再調(diào)整”的變量。模型因此呈現(xiàn)兩階段隨機(jī)優(yōu)化的結(jié)構(gòu)理解這一點(diǎn)目標(biāo)函數(shù)和約束的代碼寫法就清晰了。1.3 短期調(diào)度的時(shí)間尺度與決策內(nèi)容短期調(diào)度一般指未來24小時(shí)到72小時(shí)內(nèi)的運(yùn)行計(jì)劃時(shí)段長度常用1小時(shí)。我復(fù)現(xiàn)時(shí)采用24個(gè)時(shí)段、每時(shí)段1小時(shí)的設(shè)定這是比較經(jīng)典的做法。如果涉及日內(nèi)滾動(dòng)也可以壓到15分鐘一個(gè)時(shí)段但模型規(guī)模會成倍增加求解時(shí)間明顯變長。調(diào)度決策的內(nèi)容包括三塊第一每個(gè)水庫每個(gè)時(shí)段的發(fā)電流量和棄水流量這決定了水電出力過程第二每個(gè)時(shí)段光伏實(shí)際消納功率和棄光功率第三由水量平衡推導(dǎo)出的水庫庫容變化過程。三塊決策通過系統(tǒng)功率平衡約束耦合在一起形成完整的優(yōu)化問題。一句話總結(jié)這個(gè)模型要干什么在給定光伏出力隨機(jī)場景、來水預(yù)報(bào)和系統(tǒng)消納上限的情況下通過優(yōu)化梯級水庫的放水計(jì)劃和光伏消納方案讓整個(gè)系統(tǒng)在一天內(nèi)期望可消納電量達(dá)到最大同時(shí)保證水庫運(yùn)行約束安全。2. 數(shù)學(xué)模型拆解目標(biāo)函數(shù)、約束條件與線性化處理這一章是復(fù)現(xiàn)的核心。論文里的模型通常寫得比較緊湊但展開后其實(shí)就是一個(gè)目標(biāo)函數(shù)加五大類約束。我會把每一部分的物理含義講清楚并且給出可編程的數(shù)學(xué)形式。2.1 目標(biāo)函數(shù)怎么表示“期望值最大化”目標(biāo)函數(shù)的一般形式是最大化所有場景下系統(tǒng)可消納電量的概率加權(quán)平均值論文里還會加上棄光懲罰項(xiàng)?？梢杂孟旅娴墓奖磉_(dá)[ \max \quad \sum_{s1}^{S} \pi_s \sum_{t1}^{T} \left[ \sum_{i1}^{N} P_{i,t}^{H} \Delta t P_{t,s}^{use} \Delta t - \lambda P_{t,s}^{shed} \Delta t \right] ]各符號的含義( S ) 是光伏出力場景數(shù)( \pi_s ) 是第 ( s ) 個(gè)場景的概率所有場景概率之和為1( P_{i,t}^{H} ) 是第 ( i ) 座水電站在時(shí)段 ( t ) 的出力( P_{t,s}^{use} ) 是場景 ( s ) 下時(shí)段 ( t ) 實(shí)際消納的光伏功率( P_{t,s}^{shed} ) 是場景 ( s ) 下時(shí)段 ( t ) 的棄光功率( \lambda ) 是棄光懲罰系數(shù)取值通常遠(yuǎn)大于電價(jià)目的是讓模型優(yōu)先避免棄光。水電出力 ( P_{i,t}^{H} ) 不攜帶場景下標(biāo) ( s )因?yàn)樗娬{(diào)度計(jì)劃需要在光伏隨機(jī)性實(shí)現(xiàn)之前就確定下來這對應(yīng)兩階段隨機(jī)優(yōu)化里的第一階段決策。光伏消納和棄光變量攜帶場景下標(biāo)屬于第二階段決策。這個(gè)下標(biāo)細(xì)節(jié)對編程影響很大寫變量時(shí)一定要帶對。棄光懲罰的本質(zhì)是一個(gè)松弛機(jī)制。如果完全沒有懲罰項(xiàng)模型在極端場景下可能通過大量棄光來換取水電發(fā)電量的增加這未必符合實(shí)際運(yùn)行偏好。加上懲罰項(xiàng)后模型會自動(dòng)在“多消納光伏”和“多讓水電發(fā)電”之間做權(quán)衡實(shí)際輸出結(jié)果更貼近調(diào)度員的決策習(xí)慣。2.2 五大類約束的工程化表述我把復(fù)現(xiàn)中用到的約束整理成表格每一類都對應(yīng)明確的物理意義。表格里的形式是編程時(shí)采用的最終形式比論文里更貼近代碼實(shí)現(xiàn)。約束類型數(shù)學(xué)表達(dá)物理含義水量平衡( V_{i,t1} V_{i,t} (I_{i,t} Q_{i-1,t-\tau_i} S_{i-1,t-\tau_i} - Q_{i,t} - S_{i,t})\Delta t )上游來水、上游放水延時(shí)入庫、本庫出流共同決定庫容變化庫容邊界( V_{i}^{min} \le V_{i,t} \le V_{i}^{max} )對應(yīng)死水位和正常高水位短期調(diào)度還要加末庫容約束出流邊界( 0 \le Q_{i,t} \le Q_{i}^{max} ) ( 0 \le S_{i,t} \le S_{i}^{max} )機(jī)組過流能力和溢洪道泄流能力水電出力( P_{i,t}^{H} 9.81 \eta_i Q_{i,t}^{gen} H_{i,t} )水頭、流量和機(jī)組效率決定發(fā)電功率光伏消納( P_{t,s}^{use} P_{t,s}^{shed} P_{t,s}^{pv} )預(yù)測光伏分解為消納和棄光兩部分系統(tǒng)功率平衡( \sum_{i} P_{i,t}^{H} P_{t,s}^{use} \le P_{t}^{grid} )水電與光伏實(shí)際消納之和不超過外送或負(fù)荷上限梯級水量平衡中那個(gè) ( \tau_i ) 需要單獨(dú)說明它是第 ( i-1 ) 座水庫到第 ( i ) 座水庫的水流傳播時(shí)間。短期調(diào)度中如果流域范圍不大( \tau_i ) 通常取0到2小時(shí)。我當(dāng)時(shí)設(shè)置上游到下游的傳播時(shí)間為1小時(shí)也就是上游 ( t ) 時(shí)段的出庫最早在 ( t1 ) 時(shí)段才能進(jìn)入下游水庫的平衡方程。末庫容約束是短期調(diào)度里特別容易被忽略的一類。如果只約束庫容上下限模型大概率會把水庫放到最高水位或者排到最低水位因?yàn)閺膯稳帐找娼嵌瓤催@往往是最優(yōu)的但現(xiàn)實(shí)中這會導(dǎo)致第二天無法正常運(yùn)行。常規(guī)做法是加 ( V_{i,T}^{min, end} \le V_{i,T} \le V_{i,T}^{max, end} )讓調(diào)度周期結(jié)束時(shí)的庫容回到一個(gè)合理區(qū)間。2.3 從論文公式到可求解模型線性化三板斧水電出力特性是模型里唯一的非線性來源。( P 9.81 \eta Q H ) 中發(fā)電流量 ( Q ) 是決策變量發(fā)電水頭 ( H ) 又隨庫容和尾水位變化兩個(gè)變量相乘就形成了雙線性項(xiàng)。直接交給Gurobi這類求解器會報(bào)非凸非線性錯(cuò)誤必須先線性化。我在復(fù)現(xiàn)中試過三種做法分別適用于不同精度要求第一種是固定水頭法。對于調(diào)節(jié)能力不強(qiáng)、水位波動(dòng)小的徑流式電站可以把水頭近似看作常數(shù)出力變成流量的線性函數(shù)。這是最簡單的方式求解速度快但碰上日調(diào)節(jié)能力強(qiáng)的水庫誤差偏大。第二種是分段線性化。把發(fā)電流量分成若干段每段對應(yīng)一個(gè)線性出力區(qū)間用SOS2約束保證模型在相鄰分段之間插值。分段越多精度越高但變量數(shù)量隨之增加。實(shí)際操作中我建議水電出力曲線分5到8段即可再多對結(jié)果影響很小求解時(shí)間卻明顯上升。第三種是McCormick包絡(luò)法。對 ( z Q H ) 這個(gè)雙線性項(xiàng)用四個(gè)線性不等式約束把它包圍在凸包內(nèi)。對于短期調(diào)度這種單峰曲線問題McCormick松弛的質(zhì)量通常還不錯(cuò)。寫出通用代碼的示例是這樣的# z Q * H 的 McCormick 線性化 # Q_lb, Q_ub 為流量邊界H_lb, H_ub 為水頭邊界 m.addConstr(z Q_lb * H Q * H_lb - Q_lb * H_lb) m.addConstr(z Q_ub * H Q * H_ub - Q_ub * H_ub) m.addConstr(z Q_lb * H Q * H_ub - Q_lb * H_ub) m.addConstr(z Q_ub * H Q * H_lb - Q_ub * H_lb)這套線性化處理之后模型就變成混合整數(shù)線性規(guī)劃或者線性規(guī)劃Gurobi可以直接求解。如果模型里沒有0-1整數(shù)變量線性規(guī)劃求解速度非?？爝@也是短期調(diào)度能跑出高時(shí)效性的基礎(chǔ)。3. Python代碼架構(gòu)從數(shù)據(jù)準(zhǔn)備到Gurobi求解模型理解到位后代碼實(shí)現(xiàn)就是水到渠成的事。我的代碼結(jié)構(gòu)分為四層數(shù)據(jù)定義層、場景生成層、模型構(gòu)建層、結(jié)果輸出層。這樣分的好處是以后換算例只要改參數(shù)文件模型構(gòu)建部分完全不用動(dòng)。3.1 數(shù)據(jù)準(zhǔn)備與單位統(tǒng)一這是復(fù)現(xiàn)中最容易被坑死的一步。論文里的庫容單位經(jīng)常寫“億立方米”來水單位卻是“立方米每秒”兩個(gè)單位直接相乘水量平衡方程會差出好幾個(gè)數(shù)量級模型立刻不可行。我統(tǒng)一把所有流量換算成立方米每小時(shí)或者干脆在方程中乘以時(shí)間步長 ( \Delta t ) 來處理單位。具體做法是流量 ( m^3/s ) 乘以3600變成 ( m^3/h )再乘以時(shí)段長度1小時(shí)就得到該時(shí)段的體積。如果用半天或一小時(shí)之外的時(shí)段步長需要相應(yīng)調(diào)整3600這個(gè)系數(shù)。我測試算例采用的參數(shù)如下供參考參數(shù)數(shù)值梯級水庫數(shù)2水電站在線臺數(shù)3光伏場站裝機(jī)600 MW系統(tǒng)外送通道上限900 MW調(diào)度周期24 h水流傳播時(shí)間1 h光伏場景數(shù)5削減后棄光懲罰系數(shù)200 元/MWh數(shù)據(jù)在代碼里用字典和numpy數(shù)組組織。光伏預(yù)測出力是一個(gè) ( 24 \times S ) 的二維數(shù)組每一列對應(yīng)一條場景曲線來水預(yù)報(bào)是一維長度24的數(shù)組水庫的參數(shù)庫容上下限、最大出流、初始庫容、末庫容單獨(dú)存成字典。這樣主程序讀起來非常清晰。3.2 光伏隨機(jī)場景生成與削減光伏預(yù)測曲線一般只有一條但模型需要的是多個(gè)場景。我采用的做法是蒙特卡洛加場景削減。先生成200條場景在預(yù)測值基礎(chǔ)上疊加均值為零、方差隨時(shí)間變化的正態(tài)擾動(dòng)中午光照強(qiáng)時(shí)方差大早晚方差小。然后在場景集上做削減。場景削減最實(shí)用的方案有兩種同步回代消除法和KMeans聚類法。同步回代法更嚴(yán)謹(jǐn)它通過迭代刪除與其他場景距離最近的場景并把概率轉(zhuǎn)移給鄰近場景最終保留的場景集在概率分布意義上最優(yōu)。KMeans實(shí)現(xiàn)簡單速度快在場景數(shù)不是特別大的情況下效果足夠。我先后對比了兩者最終用的是KMeans加概率歸一的組合。from sklearn.cluster import KMeans import numpy as np # pv_samples: shape (N_scen, T)每行是一條光伏場景 k 5 kmeans KMeans(n_clustersk, random_state42).fit(pv_samples) centers kmeans.cluster_centers_ # 削減后的典型場景 labels kmeans.labels_ probs np.bincount(labels, minlengthk) / len(labels) # 各場景概率削減后得到5條典型場景曲線和對應(yīng)的概率。這5條曲線進(jìn)入目標(biāo)函數(shù)時(shí)每個(gè)場景的權(quán)重就是 ( \pi_s )。需要特別提醒的是削減后的場景數(shù)量不要太多因?yàn)槊總€(gè)場景都會復(fù)制一份光伏消納變量和棄光變量場景數(shù)從5增加到20模型規(guī)模大好幾倍求解時(shí)間成倍增長但目標(biāo)值改善往往很有限。3.3 模型構(gòu)建核心代碼模型構(gòu)建我用的是Gurobi的Python接口。變量、約束、目標(biāo)函數(shù)三部分分開寫方便調(diào)試。核心代碼如下注釋里說明了每個(gè)變量和約束的作用import gurobipy as gp from gurobipy import GRB dt 3600 # 時(shí)段長度秒 T 24 S len(probs) N 2 m gp.Model(hydro_pv_short_term) # 階段一變量水庫運(yùn)行計(jì)劃不依賴場景 V m.addVars(N, T 1, lbV_min, ubV_max, nameV) # 庫容 Q m.addVars(N, T, lb0, ubQ_max, nameQ) # 發(fā)電流量 SP m.addVars(N, T, lb0, ubSP_max, nameSP) # 棄水流量 Ph m.addVars(N, T, lb0, ubPh_max, namePh) # 水電出力 # 階段二變量光伏消納與棄光依賴場景 P_use m.addVars(T, S, lb0, nameP_use) # 實(shí)際消納光伏 P_shed m.addVars(T, S, lb0, nameP_shed) # 棄光功率 # 水量平衡約束 for i in range(N): for t in range(T): inflow I[i, t] if i 0: tau travel_time[i] if t - tau 0: inflow Q[i - 1, t - tau] SP[i - 1, t - tau] m.addConstr(V[i, t 1] V[i, t] (inflow - Q[i, t] - SP[i, t]) * dt) # 初始庫容與末庫容 for i in range(N): m.addConstr(V[i, 0] V0[i]) m.addConstr(V[i, T] Vend_min[i]) m.addConstr(V[i, T] Vend_max[i]) # 光伏消納與棄光 for t in range(T): for s in range(S): m.addConstr(P_use[t, s] P_shed[t, s] P_pv[t, s]) m.addConstr(Ph[0, t] Ph[1, t] P_use[t, s] P_grid[t]) # 目標(biāo)函數(shù)最大化期望可消納電量并加入棄光懲罰 obj gp.quicksum( probs[s] * (Ph[0, t] Ph[1, t] P_use[t, s] - penalty * P_shed[t, s]) for s in range(S) for t in range(T) ) m.setObjective(obj, GRB.MAXIMIZE) m.optimize()這個(gè)結(jié)構(gòu)里有一個(gè)細(xì)節(jié)值得注意水量平衡中的時(shí)滯處理用了條件判斷if t - tau 0。上游在調(diào)度期一開始放的水要經(jīng)過 ( tau ) 小時(shí)才能到下游調(diào)度期前幾小時(shí)下游應(yīng)該收不到上游放水。如果忽略這個(gè)時(shí)間差把上游所有時(shí)段的出庫都直接加入下游當(dāng)天水量平衡會導(dǎo)致下游水量被重復(fù)計(jì)算。3.4 求解參數(shù)設(shè)置與結(jié)果導(dǎo)出Gurobi求解器默認(rèn)參數(shù)在中小規(guī)模問題上表現(xiàn)不錯(cuò)但短期調(diào)度有時(shí)需要限制求解時(shí)間特別是在滾動(dòng)優(yōu)化場景下。我一般會加這樣幾個(gè)參數(shù)m.Params.TimeLimit 120 # 最長求解120秒 m.Params.MIPGap 0.02 # 允許2%的gap m.Params.Threads 4對于純線性規(guī)劃模型這組參數(shù)基本不會觸發(fā)TimeLimit幾秒內(nèi)就能收斂如果模型因?yàn)榧尤胝麛?shù)變量變成了MIP這兩個(gè)參數(shù)就是限制求解時(shí)長的保險(xiǎn)絲。結(jié)果導(dǎo)出我習(xí)慣直接落盤成CSV包含每個(gè)時(shí)段的水庫庫容、發(fā)電流量、棄水流量、水電出力以及每個(gè)場景下的光伏消納、棄光功率。落盤之后用pandas做匯總統(tǒng)計(jì)比直接從Gurobi變量對象里取值要方便得多。4. 算例結(jié)果分析看數(shù)字怎么驗(yàn)證模型的物理邏輯模型跑通只是第一步結(jié)果合理性驗(yàn)證才是復(fù)現(xiàn)真正收獲的地方。我把算例結(jié)果整理成表格逐項(xiàng)分析了模型輸出的物理邏輯這一步能幫你發(fā)現(xiàn)隱藏的建模錯(cuò)誤。4.1 調(diào)度結(jié)果總覽與時(shí)段級拆解下面是我在測試算例里得到的部分時(shí)段結(jié)果可以清晰地看到梯級水光互補(bǔ)的運(yùn)行規(guī)律時(shí)段光伏預(yù)測 (MW)水電總出力 (MW)系統(tǒng)外送上限 (MW)棄光 (MW)上游庫容變化08:002603806400緩降11:0045032070070緩升14:0051028071080緩升18:001805207000快速下降21:00605806400下降正午光伏出力達(dá)到峰值時(shí)水電出力明顯壓減上游水庫庫容上升相當(dāng)于把水蓄起來傍晚光伏快速退坡時(shí)水電出力大幅增加水庫放水發(fā)電補(bǔ)上缺口。這個(gè)“光伏高時(shí)蓄水、光伏低時(shí)放水”的模式正是梯級水光互補(bǔ)調(diào)度最理想的運(yùn)行狀態(tài)。如果算出來的庫容過程是白天猛放水、晚上蓄水那大概率是目標(biāo)函數(shù)或者約束寫反了需要回頭檢查系統(tǒng)功率平衡方向。4.2 梯級水庫調(diào)節(jié)過程怎么看判斷梯級聯(lián)動(dòng)是否正確的關(guān)鍵指標(biāo)是上下游水庫的庫容變化時(shí)序。模擬結(jié)果中下游水庫在上午時(shí)段收到上游水庫前一天傍晚放水的補(bǔ)充庫容呈現(xiàn)緩慢上升趨勢下午上游開始蓄水、放水量減少下游水庫入庫減少庫容轉(zhuǎn)為下降。這個(gè)過程與水流傳播時(shí)間 ( \tau ) 的設(shè)置完全對應(yīng)。還有一個(gè)驗(yàn)證點(diǎn)是棄水行為。梯級調(diào)度中如果上游來水較大而上游庫容不足或下游無法承接會出現(xiàn)棄水。但模型應(yīng)該優(yōu)先利用水頭高、耗水率低的水庫發(fā)電而不是讓水白白流走。測試算例里棄水主要集中在深夜來水高峰時(shí)段且下游庫容接近上限時(shí)這正是物理上合理的棄水場景。4.3 場景削減對結(jié)果的影響我專門對比了三種方案只使用點(diǎn)預(yù)測、使用200條原始場景、使用削減后的5條場景。結(jié)果能清楚顯示場景數(shù)量與優(yōu)化質(zhì)量的關(guān)系方案場景數(shù)目標(biāo)值 (萬kWh)求解時(shí)間 (s)棄光率點(diǎn)預(yù)測115200.812.4%原始場景200149638.510.1%削減場景515031.210.4%點(diǎn)預(yù)測方案目標(biāo)值最高是因?yàn)槟Ｐ桶阉泄夥A(yù)測都當(dāng)作確定值敢于在系統(tǒng)功率平衡約束的邊緣試探但一旦實(shí)際光伏低于預(yù)測就會大量棄水或棄光。200條原始場景的目標(biāo)值最低因?yàn)闃O端場景拉低了期望值但也因此更保守、棄光率最低。削減后的5條場景在目標(biāo)值和求解時(shí)間之間取得了很好的平衡。這個(gè)對比能回答復(fù)現(xiàn)時(shí)最常見的問題“場景削減有必要嗎”答案是有必要但場景數(shù)不是越多越好。對短期調(diào)度來說5到10條典型場景就能覆蓋主要的不確定性特征再多只增加計(jì)算負(fù)擔(dān)。5. 復(fù)現(xiàn)路上的坑與調(diào)試經(jīng)驗(yàn)這一章是全文最值錢的部分。模型和代碼最終能跑通很大程度上是因?yàn)槲野衙總€(gè)坑都詳細(xì)記錄并逐一解決了。每個(gè)坑背后都對應(yīng)一條建模或編程上的通用經(jīng)驗(yàn)換到別的調(diào)度模型同樣適用。5.1 不可行模型先查單位再查約束沖突第一次求解直接報(bào)Model is infeasible我當(dāng)時(shí)第一反應(yīng)是約束寫錯(cuò)了把水量平衡、庫容上下限來回檢查了好幾遍全都沒有問題。最后發(fā)現(xiàn)問題出在單位混用上游來水用的單位是 ( m^3/s ) 乘以 ( dt ) 得到了立方米而庫容上下限用的卻是論文里的億立方米兩個(gè)量級直接相差 ( 10^8 )模型當(dāng)然無解。這個(gè)經(jīng)歷想說明的是排錯(cuò)要有順序。遇到不可行第一步查所有參數(shù)的物理單位是否一致把所有流量統(tǒng)一成 ( m^3/s )把所有水量統(tǒng)一成 ( m^3 ) 或 ( 10^4 m^3 )。第二步查邊界約束是否有交叉沖突比如庫容下限過高加上初始庫容過低模型連可行解都不存在。第三步才去查約束內(nèi)部的等式和不等式方向。Gurobi提供一個(gè)很有用的工具模型不可行后可以調(diào)用m.computeIIS()找出不可行約束的最小集合。復(fù)現(xiàn)調(diào)試階段用好這個(gè)功能能節(jié)省大量時(shí)間。5.2 場景削減要保留物理極值KMeans做場景削減有一個(gè)潛在風(fēng)險(xiǎn)聚類中心會被“均值化”極端場景的特征在削減后可能被抹平。比如光伏出力特別低、可能引起系統(tǒng)功率缺口的場景如果被聚類中心平均掉了模型得到的調(diào)度計(jì)劃就缺乏對極端情況的預(yù)判能力。我的處理辦法是在削減前先保留幾條極端場景固定參與優(yōu)化不允許聚類把它們的特征平滑掉。具體做法是把極端場景單獨(dú)挑出來賦予一個(gè)固定概率然后再對剩余場景做KMeans。這樣模型既保留了不確定性概率分布的典型代表又不丟失最危險(xiǎn)的情況。5.3 分段線性化的分段點(diǎn)密度不是越密越好分段線性化處理水電出力特性時(shí)我一開始把發(fā)電流量分成了20段想知道精度能提升多少。結(jié)果模型直接變成了混合整數(shù)規(guī)劃變量數(shù)暴增求解時(shí)間從幾秒變成了十幾分鐘而目標(biāo)值只變化了不到0.1%。后來我把分段數(shù)調(diào)到5到8段求解時(shí)間回到1秒以內(nèi)精度損失完全可以接受。原因是短期調(diào)度中水電出力在大多數(shù)時(shí)段都在低流量段運(yùn)行高流量段的分段點(diǎn)幾乎用不到堆太多分段只會增加無意義的計(jì)算量。分段點(diǎn)的分布也值得注意流量較小的區(qū)間分段要密一些大流量區(qū)間可以稀疏一些這樣擬合出來的出力曲線更貼近真實(shí)特性。5.4 求解精度與速度的權(quán)衡Gap和TimeLimit短期優(yōu)化調(diào)度模型如果被投入到日內(nèi)滾動(dòng)運(yùn)行中每次求解的時(shí)間窗口可能只有幾分鐘不可能無限等待最優(yōu)解。這時(shí)候Gurobi的MIPGap參數(shù)就很有用。我把gap設(shè)置為2%求解器一旦找到足夠接近最優(yōu)的解就會提前停止大大縮短了求解時(shí)間。對于純線性規(guī)劃模型這個(gè)設(shè)置不產(chǎn)生什么影響因?yàn)榫€性規(guī)劃求解本身很快。如果模型引入了0-1變量比如機(jī)組開停機(jī)狀態(tài)MIPGap就非常關(guān)鍵了。實(shí)際運(yùn)行時(shí)可以把gap從5%開始嘗試觀察目標(biāo)值變化逐步收緊到1%或2%找到一個(gè)速度和精度都能接受的平衡點(diǎn)。5.5 隨機(jī)種子必須固定這是我最后補(bǔ)的一個(gè)細(xì)節(jié)但重要性不亞于前面任何一個(gè)坑。場景生成和KMeans削減都有隨機(jī)性如果每次運(yùn)行結(jié)果都不一樣后續(xù)分析根本沒法做。我在所有使用隨機(jī)數(shù)的地方都固定了種子包括numpy的np.random.seed(42)和KMeans的random_state42。這樣同一套數(shù)據(jù)每次跑出來的結(jié)果完全一致也方便在校驗(yàn)?zāi)Ｐ蜁r(shí)對比不同參數(shù)下的變化。另外一個(gè)相關(guān)的問題是論文復(fù)現(xiàn)時(shí)如果對方?jīng)]有給隨機(jī)種子你換一臺機(jī)器或者換一個(gè)Python版本生成的場景都可能不同。所以我自己復(fù)現(xiàn)時(shí)一般不依賴遠(yuǎn)程隨機(jī)數(shù)而是把生成好的場景集直接保存成CSV作為固定的輸入文件。這樣場景數(shù)據(jù)穩(wěn)定后續(xù)所有調(diào)試都基于同一套數(shù)據(jù)展開避免因?yàn)閳鼍白兓瘉砘卣垓v。最后再分享一個(gè)小經(jīng)驗(yàn)。整個(gè)模型復(fù)現(xiàn)過程中最耗時(shí)間的不是最終寫代碼那一步而是想清楚“每個(gè)變量到底屬于第一階段還是第二階段”。這個(gè)判斷直接影響變量要不要帶場景下標(biāo)一旦搞反目標(biāo)函數(shù)表面上還能求出值但結(jié)果會完全違背物理規(guī)律。我建議調(diào)試時(shí)先跑一個(gè)單場景版本的模型確認(rèn)水電發(fā)電和光伏消納的時(shí)序曲線符合直覺再加場景削減功能。單場景跑通多場景只是重復(fù)乘上概率權(quán)重而已出問題的概率會小很多。