韌性提升中的移動電源動態(tài)調(diào)度:Matlab+Yalmip建模與求解實(shí)踐)
最近在復(fù)現(xiàn)一篇SCI一區(qū)期刊上關(guān)于配電網(wǎng)韌性提升的文章重點(diǎn)把其中應(yīng)急移動電源Mobile Power SourceMPS的動態(tài)調(diào)度部分在Matlab里完整跑通了。這個方向現(xiàn)在確實(shí)很熱——臺風(fēng)、暴雨、覆冰這些極端天氣一上來配電網(wǎng)最容易發(fā)生大面積停電靠傳統(tǒng)搶修恢復(fù)太慢所以很多研究都把重心轉(zhuǎn)向“災(zāi)前預(yù)配置災(zāi)中動態(tài)調(diào)度”的組合拳。上篇預(yù)配置解決的是“移動電源提前停在哪里、配多大容量”下篇動態(tài)調(diào)度解決的是“災(zāi)害發(fā)生后故障點(diǎn)已經(jīng)明確這些移動電源怎么陸續(xù)趕到災(zāi)區(qū)、優(yōu)先給哪些負(fù)荷供電”。這篇博文就是要把動態(tài)調(diào)度這部分原理解透并把Matlab復(fù)現(xiàn)的過程記錄下來。想搞配電網(wǎng)韌性、移動儲能、應(yīng)急電源優(yōu)化調(diào)度的研究生或者剛從傳統(tǒng)電力系統(tǒng)優(yōu)化轉(zhuǎn)過來的同學(xué)看完應(yīng)該能少走不少彎路。復(fù)現(xiàn)之前我先把整個問題拆了一遍。最深的感受是這類文章表面上在寫“優(yōu)化算法”實(shí)際上真正的難點(diǎn)在于約束建模尤其是移動電源的空間轉(zhuǎn)移和時間過程怎么統(tǒng)一到同一個模型里。這一篇我會把動態(tài)調(diào)度模型、Yalmip建模技巧、求解器調(diào)參、常見坑一次講清楚。1. 先把問題邊界劃清楚動態(tài)調(diào)度到底在調(diào)度什么1.1 預(yù)配置和動態(tài)調(diào)度如何銜接從題目就能看出來“預(yù)配置”和“動態(tài)調(diào)度”是同一個研究鏈條里的兩個環(huán)節(jié)。預(yù)配置階段一般在災(zāi)前完成輸入是災(zāi)害預(yù)測信息或者歷史典型場景輸出是MPS的初始部署位置、數(shù)量和容量。到了動態(tài)調(diào)度階段假設(shè)災(zāi)害已經(jīng)發(fā)生調(diào)度中心拿到故障線路、停電負(fù)荷、可用MPS位置這些信息需要決定每個MPS后續(xù)怎么移動、何時接入、給誰供電。復(fù)現(xiàn)的時候這兩個階段在代碼層面其實(shí)是分開的。上篇的結(jié)果就是一個靜態(tài)數(shù)組比如mps(:,1)表示每個MPS的初始節(jié)點(diǎn)編號mps_cap(m)表示容量。到了下篇動態(tài)調(diào)度這些值直接作為固定參數(shù)寫進(jìn)約束里。這里有一個容易忽視的銜接問題如果上篇的預(yù)配置方案是隨機(jī)場景優(yōu)化的那么下篇在具體故障場景下做調(diào)度時初始位置未必是“最優(yōu)”的但模型必須接受這個設(shè)定。這也是很多論文里說的“兩階段決策”結(jié)構(gòu)——第一階段做決策時不依賴第二階段的具體實(shí)現(xiàn)第二階段則基于第一階段的決策做適應(yīng)性調(diào)整。我在代碼里用了一個固定的初始化結(jié)構(gòu)mps_init_pos [14, 25, 30]; % 預(yù)配置位置來自上篇結(jié)果 mps_capacity [500, 300, 300]; % 單位kWh應(yīng)急電源容量如果你的復(fù)現(xiàn)目標(biāo)里沒有上篇代碼完全可以手動設(shè)置一個合理初始位置不影響下篇調(diào)度邏輯的驗(yàn)證。1.2 動態(tài)調(diào)度的三個關(guān)鍵決策動態(tài)調(diào)度模型的核心可以拆成三個“要決策什么”第一分配給誰。每個MPS容量有限、位置不同而每個故障節(jié)點(diǎn)的負(fù)荷重要程度和恢復(fù)價值不同所以要把MPS和節(jié)點(diǎn)做匹配。第二何時到達(dá)。MPS從初始位置到目標(biāo)節(jié)點(diǎn)需要通行時間這個時間取決于節(jié)點(diǎn)之間的道路距離和移動速度。早到一個小時關(guān)鍵負(fù)荷就少停一個小時所以在時間軸上做優(yōu)化比單純空間分配更有價值。第三供多少功率。到了節(jié)點(diǎn)之后MPS在同一時刻能輸出的功率受容量限制而且總放電能量也受儲能容量限制不是想供多少就供多少。這三層決策耦合在一個時間軸上本質(zhì)上是一個時空網(wǎng)絡(luò)流問題。我用一個生活化類比輔助理解幾個水電工師傅從各自家里出發(fā)要處理多個不同優(yōu)先級的故障點(diǎn)到達(dá)每個點(diǎn)需要不同時間處理每個點(diǎn)消耗不同工時和物料怎么排最合理。這個類比雖然不完全等價但能很快抓住MPS調(diào)度的本質(zhì)——空間轉(zhuǎn)移有代價時間分配有先后負(fù)荷恢復(fù)有優(yōu)先級。1.3 目標(biāo)函數(shù)怎么定配電網(wǎng)韌性提升的量化指標(biāo)有很多種動態(tài)調(diào)度里最常用的是“加權(quán)負(fù)荷恢復(fù)量最大化”。這里的權(quán)重直接對應(yīng)負(fù)荷等級比如醫(yī)院、通信基站、供水設(shè)施的權(quán)重最高普通居民負(fù)荷權(quán)重低一些。極端災(zāi)害下無法恢復(fù)全部負(fù)荷所以優(yōu)先保障重要負(fù)荷是模型的核心邏輯。有的論文會在目標(biāo)函數(shù)里加一個懲罰項比如MPS移動距離的懲罰或者移動次數(shù)的懲罰目的是防止模型為了極小收益讓MPS頻繁移動。復(fù)現(xiàn)時可以先不加懲罰項跑通之后再對比加懲罰項的結(jié)果差異這本身也是一種很好的模型驗(yàn)證手段。目標(biāo)函數(shù)的標(biāo)準(zhǔn)形式是objective -sum(omega(n) * p_load(n,t) * dt) lambda * sum(dep(...));這里omega是負(fù)荷權(quán)重dt是時間步長dep是移動弧變量lambda是移動成本系數(shù)。注意Yalmip默認(rèn)是最小化所以目標(biāo)函數(shù)前面加負(fù)號或者把Yalmip的sense設(shè)成minimize之后按上式寫。我習(xí)慣直接寫Objective -sum(sum(omega * pL)) lambda * sum(dep(:));其中pL是恢復(fù)負(fù)荷矩陣omega是列向量。2. 動態(tài)調(diào)度的核心數(shù)學(xué)模型從時間擴(kuò)展圖入手2.1 時間擴(kuò)展圖建模思路動態(tài)調(diào)度里最容易卡住的地方是MPS在空間上移動需要時間而模型中所有決策都發(fā)生在離散時段上怎么把“移動過程”表達(dá)清楚。我強(qiáng)烈建議用時間擴(kuò)展圖time-expanded network來建模。時間擴(kuò)展圖的思想很簡單把每個時段t的配電網(wǎng)節(jié)點(diǎn)i看作一個時空節(jié)點(diǎn)(i,t)每個MPS在時空節(jié)點(diǎn)之間移動。節(jié)點(diǎn)i在t時段停留就相當(dāng)于占用時空節(jié)點(diǎn)(i,t)從一個節(jié)點(diǎn)i移動到節(jié)點(diǎn)j需要T_ij個時段就相當(dāng)于在時空圖上從(i,t)走到(j,tT_ij)的一條弧。這樣整個調(diào)度問題就變成了在時空網(wǎng)絡(luò)上為每臺MPS找一條從初始節(jié)點(diǎn)出發(fā)的路徑路徑經(jīng)過的節(jié)點(diǎn)可以在對應(yīng)時段為負(fù)荷供電。這種建模方式好處很明顯一是邏輯清晰移動和供電兩個狀態(tài)天然分開二是約束容易寫成線性形式配合二進(jìn)制變量直接交給求解器三是調(diào)試直觀把時空路徑畫出來就是MPS的完整軌跡。2.2 位置、移動與供電三類變量怎么定義復(fù)現(xiàn)中我定義了三組核心0-1變量和一組連續(xù)變量pos(m,n,t)MPS m在時段t是否位于節(jié)點(diǎn)n這是“停留可用”狀態(tài)。dep(m,i,j,t)MPS m在時段t是否處于從i到j(luò)的移動過程中表示移動占用狀態(tài)。s(m,n,t)MPS m在時段t是否在節(jié)點(diǎn)n處接入電網(wǎng)供電它必須約束在pos1的前提下。p(m,n,t)MPS m在時段t向節(jié)點(diǎn)n注入的有功功率連續(xù)非負(fù)變量。變量多了之后矩陣維度容易搞混我的習(xí)慣是先固定維度順序?yàn)?M,N,T)也就是MPS數(shù)×節(jié)點(diǎn)數(shù)×?xí)r段數(shù)。下面的代碼是定義變量的示例pos binvar(Nmps, Nnode, T, full); dep binvar(Nmps, Nnode, Nnode, T, full); s binvar(Nmps, Nnode, T, full); p sdpvar(Nmps, Nnode, T, full); pL sdpvar(Nnode, T, full);有個細(xì)節(jié)p(m,n,t)的變量在MPS沒有到達(dá)節(jié)點(diǎn)n時沒有意義完全可以用固定值0代替。變量越多求解越慢所以最好做一次變量裁剪。一般做法是只對MPS初始位置、故障節(jié)點(diǎn)以及它們附近的潛在接入節(jié)點(diǎn)保留p變量其他位置直接置0能顯著降低模型規(guī)模。2.3 核心約束狀態(tài)互斥與移動時延動態(tài)調(diào)度最核心的約束是“一個MPS在任意時刻只能處于一種狀態(tài)”。具體來說每個MPS在時段t要么停留在某個節(jié)點(diǎn)要么在某個移動過程中不能同時出現(xiàn)在兩個地方也不能邊移動邊供電。寫成約束就是for m 1:Nmps for t 1:T % 狀態(tài)互斥停留 移動占用 恰好一種 [con, con] ... [con, sum(pos(m,:,t), 2) sum(sum(dep(m,:,:,t)))]; Constraints [Constraints, sum(pos(m,:,t),2) sum(sum(dep(m,:,:,t))) 1]; % 只有停留可用時才能接入供電 Constraints [Constraints, s(m,:,t) pos(m,:,t)]; % 供電功率受接入狀態(tài)和容量限制 Constraints [Constraints, p(m,:,t) MPS_Pmax(m) * s(m,:,t)]; end end這里有一個在復(fù)現(xiàn)時容易踩的坑如果MPS從節(jié)點(diǎn)i到節(jié)點(diǎn)j的移動時間T_ij大于1個時段那么dep變量表示的是“正在移動中”這個狀態(tài)而不是“開始移動”的事件。這樣設(shè)計的好處是狀態(tài)互斥約束可以統(tǒng)一寫成“停留移動1”而不需要關(guān)心移動是從哪個時刻開始的。但代價是需要額外加一條“移動結(jié)束時才能到達(dá)節(jié)點(diǎn)j”的約束% 從i出發(fā)到達(dá)j需要T_ij個時段只有完成移動后才能在j點(diǎn)出現(xiàn) % 即 dep(m,i,j,t) 1 時pos(m,j,tT_ij) 1 % 實(shí)際寫成線性表達(dá)式 for m 1:Nmps for t 1:T for i 1:Nnode for j 1:Nnode if TravelTime(i,j) 0 t TravelTime(i,j) T Constraints [Constraints, ... pos(m,j,tTravelTime(i,j)) dep(m,i,j,t)]; end end end end end同時為了防止提前到達(dá)需要在中間時段把pos鎖死為0。不過因?yàn)闋顟B(tài)互斥已經(jīng)把pos和dep的關(guān)系綁定了當(dāng)dep1時所有pos都必須為0這個約束自然就保證了移動期間不會出現(xiàn)在除終點(diǎn)之外的其他節(jié)點(diǎn)。如果你的模型里dep含義是“開始移動”而不是“正在移動”那中間時段的pos鎖定就必須額外加這點(diǎn)一定要看清原文的定義方式。2.4 配電網(wǎng)潮流約束與二階錐松弛MPS動態(tài)調(diào)度不是簡單的路徑規(guī)劃接入配電網(wǎng)之后功率注入會改變潮流分布所以還要考慮潮流約束。對于輻射狀配電網(wǎng)最常用的是DistFlow模型。以IEEE 33節(jié)點(diǎn)系統(tǒng)為例每條支路(i,j)的DistFlow方程可以寫成P_ij(t) - sum(P_jk(t)) - r_ij * l_ij(t) p_mps(j,t) p_load(j,t)這里P_ij是流入支路的有功l_ij是支路電流的平方p_mps是MPS注入功率。電壓約束用U_i - U_j ≥ 2(r_ij P_ij x_ij Q_ij) - (r_ij^2 x_ij^2) l_ij 來近似其中U_i是節(jié)點(diǎn)i電壓幅值的平方。DistFlow本身包含非線性項P_ij^2 Q_ij^2直接求解很麻煩。好在大量研究已經(jīng)證明在目標(biāo)函數(shù)單調(diào)的前提下可以把二次等式松弛為二階錐不等式即2P_ij^2 2Q_ij^2 (l_ij - U_i)^2 ≤ (l_ij U_i)^2用Yalmip表達(dá)這個二階錐約束很簡潔Constraints [Constraints, ... norm([2*Pij(b,t); 2*Qij(b,t); lij(b,t)-U(i,t)], 2) lij(b,t)U(i,t)];有人會問松弛之后還是原問題最優(yōu)解嗎這就是所謂“精確凸松弛”問題。大多數(shù)配電網(wǎng)輻射狀且目標(biāo)函數(shù)是恢復(fù)負(fù)荷最大化的場景下松弛是緊的也就是松弛解就是原問題的全局最優(yōu)解。復(fù)現(xiàn)時如果你發(fā)現(xiàn)結(jié)果里某條支路的二階錐約束不緊要檢查是不是目標(biāo)函數(shù)里加了不恰當(dāng)?shù)囊苿討土P項或者負(fù)荷權(quán)重設(shè)置導(dǎo)致目標(biāo)函數(shù)對潮流不敏感。判斷方法很簡單看每個支路對應(yīng)約束左右兩邊的差值如果差很小比如小于1e-4說明松弛緊結(jié)果可信。3. MatlabYalmip復(fù)現(xiàn)的完整流程3.1 環(huán)境準(zhǔn)備與數(shù)據(jù)組織先明確一下環(huán)境我用的Matlab R2022b加上Yalmip求解器用的Gurobi。如果沒有Gurobi用Cplex或者M(jìn)osek也行如果是小規(guī)模算例甚至Cbc這種開源求解器也能跑只是速度慢一些。Yalmip的安裝和配置這里不展開網(wǎng)上資料很多。代碼組織上我分成下面幾個文件main_mps_dispatch.m % 主程序 data_ieee33.m % 配電網(wǎng)參數(shù)與負(fù)荷數(shù)據(jù) gen_travel_time.m % 生成節(jié)點(diǎn)間通行時間矩陣 build_model.m % 構(gòu)建優(yōu)化模型并求解 plot_results.m % 可視化結(jié)果數(shù)據(jù)準(zhǔn)備是關(guān)鍵的一步。IEEE 33節(jié)點(diǎn)系統(tǒng)的支路參數(shù)、負(fù)荷數(shù)據(jù)都有公開版本但不同論文使用的基準(zhǔn)容量、電壓等級可能不同復(fù)現(xiàn)前一定要和原文核對清楚尤其是功率基準(zhǔn)值和時間步長dt。我復(fù)現(xiàn)時用1小時作為時間步長總調(diào)度周期取24小時相當(dāng)于覆蓋一個典型災(zāi)后恢復(fù)日。通行時間矩陣的生成也需要提前處理。最簡化方式是用節(jié)點(diǎn)之間的線路長度除以MPS移動速度。更精細(xì)的方式是用道路網(wǎng)絡(luò)距離但一般論文里不會給那么詳細(xì)的數(shù)據(jù)所以用線路長度代替就行。關(guān)鍵是TravelTime(i,j)矩陣要在建模之前就生成好并且保證對角線為0非對角線為正整數(shù)。% 生成通行時間矩陣 TravelTime zeros(Nnode, Nnode); for i 1:Nnode for j 1:Nnode if i ~ j dist norm(bus_coord(i,:) - bus_coord(j,:)); TravelTime(i,j) max(1, round(dist / MPS_speed)); end end end這里有個經(jīng)驗(yàn)移動時間必須取整到時間步長的整數(shù)倍。如果步長是1小時而兩節(jié)點(diǎn)之間行車需要2.4小時向上取整成3小時會讓MPS“遲到”一些向下取整會導(dǎo)致移動時間不合理。優(yōu)先向上取整保證物理可實(shí)現(xiàn)。3.2 模型構(gòu)建核心代碼build_model.m是整個復(fù)現(xiàn)的核心。我先把所有變量定義好然后按四類約束添加MPS運(yùn)行約束、配電網(wǎng)潮流約束、負(fù)荷恢復(fù)約束、目標(biāo)函數(shù)。MPS運(yùn)行約束這一塊除了前面講的狀態(tài)互斥和移動時延還有兩個容易被忽略的點(diǎn)。第一MPS總放電能量不能超過容量。寫成% MPS儲能容量約束 for m 1:Nmps Constraints [Constraints, ... sum(sum(p(m,:,:))) * dt mps_capacity(m)]; end第二MPS在初始時段必須位于初始位置這是一個初始條件約束for m 1:Nmps Constraints [Constraints, pos(m, mps_init_pos(m), 1) 1]; end負(fù)荷恢復(fù)約束的關(guān)鍵在于恢復(fù)量不能超過原負(fù)荷需求。我一般設(shè)定一個連續(xù)恢復(fù)比例變量alpha(n,t)取值范圍0~1然后pL(n,t)alpha(n,t)*load_demand(n,t)。因?yàn)闊o功負(fù)荷也要同步恢復(fù)所以無功恢復(fù)量用功率因數(shù)綁定。潮流部分我以DistFlow為基礎(chǔ)按3.3節(jié)的二階錐約束逐一添加。注意每條支路在t時段都要加一套約束循環(huán)嵌套要注意效率。Matlab里用for循環(huán)直接加約束在規(guī)模不大33節(jié)點(diǎn)×24時段時可接受如果系統(tǒng)擴(kuò)大到幾百節(jié)點(diǎn)就要考慮用矩陣方式一次性添加約束否則建模時間會暴漲。3.3 求解與結(jié)果校驗(yàn)求解設(shè)置方面我一般配置MIPGap為1%TimeLimit為3600秒。這里有個容易忽略的點(diǎn)MISOCP問題如果Gap太緊比如0.01%以下求解時間會指數(shù)級上升而對復(fù)現(xiàn)驗(yàn)證來說1%的精度已經(jīng)足夠判斷模型邏輯是否正確。ops sdpsettings(solver,gurobi,... verbose,2,... gurobi.MIPGap,0.01,... gurobi.TimeLimit,3600,... gurobi.Threads,8); result optimize(Constraints, Objective, ops);求解完成后第一件事不是看結(jié)果而是檢查求解狀態(tài)。如果result.problem不為0說明模型有問題要回看約束。如果求解成功把三個關(guān)鍵結(jié)果畫出來恢復(fù)負(fù)荷曲線、MPS軌跡、電壓分布。畫MPS軌跡我用的方式是對pos變量取最大值的索引figure; hold on; for m 1:Nmps [~, loc] max(pos(m,:,:), [], 2); loc squeeze(loc); stairs(1:T, loc, LineWidth, 1.8); end xlabel(時段/h); ylabel(MPS所在節(jié)點(diǎn)編號);這張圖能很直觀地看出每臺MPS什么時候在哪如果曲線出現(xiàn)跳變比如從節(jié)點(diǎn)14直接跳到節(jié)點(diǎn)25而中間時段沒有經(jīng)過節(jié)點(diǎn)那一定是移動時延約束或狀態(tài)互斥約束寫錯了。4. 求解過程中最容易被坑的幾個地方4.1 時間步長與移動耗時矩陣的統(tǒng)一我最早跑模型的時候TravelTime矩陣用的是浮點(diǎn)數(shù)結(jié)果求解器瘋狂報數(shù)值錯誤。后來才發(fā)現(xiàn)因?yàn)椴介L是整數(shù)小時移動時間必須取整才能匹配到整數(shù)時段。這里建議在生成TravelTime矩陣之后立刻檢查assert(all(all(TravelTime round(TravelTime))), 移動時間必須為整數(shù));另外如果某兩個節(jié)點(diǎn)之間的移動時間大于總調(diào)度周期T那這條移動弧其實(shí)是沒有意義的可以直接禁用避免增加無用的二進(jìn)制變量。4.2 big-M不要取得過大在MPS接入功率約束里我用了p(m,n,t) MPS_Pmax(m) * s(m,n,t)這種形式這里MPS_Pmax就是天然的上界不需要額外再設(shè)大M。但有些約束必須用繼電器形式處理比如移動時延約束里pos(m,j,tTravelTime) dep(m,i,j,t)這種邏輯約束本質(zhì)是“上升沿檢測”不需要大M。真正需要大M的地方是潮流約束或負(fù)荷恢復(fù)約束里某段非線性邏輯建議大M取該變量的物理上限乘以1.1不要拍腦袋填1e6。大M太大不僅數(shù)值條件差還會讓MIP問題更難解。4.3 目標(biāo)函數(shù)量綱與負(fù)荷權(quán)重目標(biāo)函數(shù)里負(fù)荷權(quán)重omega和MPS移動懲罰lambda的量綱如果不一致可能導(dǎo)致優(yōu)化結(jié)果完全偏向某一邊。比如omega以元/kWh為單位lambda以元/次移動為單位那么lambda設(shè)為1e3可能就太小了模型會頻繁移動MPS。合理做法是跑一組靈敏度分析固定omega把lambda從0開始逐漸增大觀察MPS移動次數(shù)和恢復(fù)負(fù)荷的變化曲線選一個折中值。這個步驟也是論文里常說的“參數(shù)敏感性分析”復(fù)現(xiàn)時順手做出來還能當(dāng)額外貢獻(xiàn)。4.4 Gurobi求解MISOCP的調(diào)參技巧Gurobi求解二階錐MIP問題時有時候會卡在某個整數(shù)節(jié)點(diǎn)上半天不動。我的經(jīng)驗(yàn)是把NumericFocus調(diào)到1或2開啟Presolve的強(qiáng)約束檢測然后MIPGap設(shè)到0.01。如果問題規(guī)模太大還可以先把潮流約束的SOCP松弛放寬一點(diǎn)求一個近似可行解然后再逐步收緊相當(dāng)于熱啟動策略。具體代碼如下ops.gurobi.NumericFocus 2; ops.gurobi.Presolve 2; ops.gurobi.MIPFocus 2;再補(bǔ)充一個實(shí)用技巧先用一個小規(guī)模算例比如IEEE 13節(jié)點(diǎn)、6個時段驗(yàn)證模型邏輯跑通之后再放大到33節(jié)點(diǎn)、24時段。這樣排查約束錯誤的時間能節(jié)省一大半。5. 常見問題與排查速查表5.1 典型求解報錯與處理現(xiàn)象可能原因解決思路Yalmip報“Solver not found”求解器路徑?jīng)]有配置好運(yùn)行yalmiptest檢查求解器是否可用確認(rèn)Gurobi/Cplex已安裝且Yalmip能找到模型顯示“Infeasible problem”移動時延約束或狀態(tài)互斥約束過緊先打開約束松弛模式逐步放寬大M或增加MPS容量定位第幾類約束導(dǎo)致無解Gurobi警告“Numerical trouble”變量尺度差異過大常見是電壓標(biāo)幺值和小數(shù)約束混用全部使用標(biāo)幺值U用電壓標(biāo)幺值的平方功率用基準(zhǔn)功率歸一化求解時間很長且Gap不下降MIP變量太多MISOCP分支搜索困難減小時間步長數(shù)量、裁剪無效移動弧、設(shè)置MIPGap為0.02或0.03先求參考解5.2 結(jié)果異常的邏輯檢查求解成功但結(jié)果不合理的現(xiàn)象往往是建模邏輯漏洞比報錯更難發(fā)現(xiàn)。我把檢查步驟總結(jié)成一張清單第一打印每個時段的pos變量看MPS是否出現(xiàn)“瞬移”。如果某一臺MPS在t時刻還在節(jié)點(diǎn)14t1時刻突然出現(xiàn)在節(jié)點(diǎn)25那說明移動時延約束沒起作用。第二檢查“移動中供電”是否發(fā)生。把s(m,n,t)和dep(m,i,j,t)疊加在同一條時間軸上如果同一時段既有移動又有接入供電說明狀態(tài)互斥約束寫錯了。第三校驗(yàn)?zāi)芰科胶?。統(tǒng)計每臺MPS所有時段注入的總能量是否嚴(yán)格小于等于容量。如果相等或者超出大概率是p變量的時段統(tǒng)計口徑有問題。第四對比負(fù)荷恢復(fù)率。單時段最大恢復(fù)負(fù)荷不能超過該節(jié)點(diǎn)原負(fù)荷需求也不能超過MPS最大功率。超了說明pL約束或恢復(fù)比例約束寫松了。第五畫電壓曲線。如果某個節(jié)點(diǎn)電壓長時間越限要檢查二階錐約束的系數(shù)方向是否正確DistFlow公式里的r和x是否搞混。這里特別想強(qiáng)調(diào)第一點(diǎn)和第二點(diǎn)。我復(fù)現(xiàn)這個題目時最耗時的排查就是發(fā)現(xiàn)MPS會在移動過程中供電物理上完全說不通但目標(biāo)函數(shù)很喜歡這種“免費(fèi)供電”所以模型會鉆這個空子。這類邏輯漏洞必須靠狀態(tài)變量可視化來抓不能只看目標(biāo)函數(shù)數(shù)值是不是合理。5.3 與預(yù)配置結(jié)果的銜接問題還有一個常見問題出現(xiàn)在上篇預(yù)配置和下篇動態(tài)調(diào)度的數(shù)據(jù)傳遞上。如果預(yù)配置輸出的MPS初始位置和容量超出了動態(tài)調(diào)度允許的接入范圍比如某個MPS初始位置在一個不能接入的節(jié)點(diǎn)那動態(tài)調(diào)度就會直接無解。我在代碼里加了一個斷言assert(all(ismember(mps_init_pos, candidate_nodes)), 初始位置不在可用節(jié)點(diǎn)集合中);同時預(yù)配置節(jié)點(diǎn)如果是變電站節(jié)點(diǎn)或者沒有負(fù)荷的聯(lián)絡(luò)節(jié)點(diǎn)也要確認(rèn)配電網(wǎng)模型里有對應(yīng)的潮流變量。否則MPS到了那個節(jié)點(diǎn)卻無處可接模型約束就會很怪。6. 寫在最后復(fù)現(xiàn)這個題目的一點(diǎn)體會復(fù)現(xiàn)SCI一區(qū)論文最重要的是不要上來就動鍵盤。我最早拿到這個題目的時候以為難點(diǎn)在目標(biāo)函數(shù)和求解器調(diào)參結(jié)果真正花時間的幾乎都在約束建模尤其是MPS的時空轉(zhuǎn)移邏輯。如果你也打算復(fù)現(xiàn)類似文章我強(qiáng)烈建議先把變量定義、約束類型、時間步長、節(jié)點(diǎn)編號全部用注釋寫在代碼頭部然后按“MPS層—負(fù)荷層—潮流層”的順序一層層加約束每加一組約束就跑一遍看可行性這樣能快速定位問題。另一個體會是這個模型后續(xù)可以擴(kuò)展的方向真的很多。比如把MPS換成柴油發(fā)電車和移動儲能車的混合車隊考慮它們的成本差異和啟動時間或者把移動時間做成隨機(jī)變量用魯棒優(yōu)化或者機(jī)會約束表達(dá)交通不確定性再進(jìn)一步還可以把MPS調(diào)度和搶修隊伍調(diào)度統(tǒng)一建模實(shí)現(xiàn)“先復(fù)電、后修復(fù)”的協(xié)同恢復(fù)策略。這套動態(tài)調(diào)度的底層框架只要建好了擴(kuò)展起來都是自然的。至少我現(xiàn)在好幾個后續(xù)想法都是在這套代碼基礎(chǔ)上改出來的。