1. 先搞清楚這道題到底在解什么做配電網(wǎng)多時段優(yōu)化的人幾乎都繞不開“YALMIP CPLEX 二階錐松弛”這一套組合。我在IEEE33和PG69兩個經(jīng)典算例上把這套流程完整落過地最終算的是多時間斷面的配電網(wǎng)最優(yōu)潮流既包含時序負荷、儲能的充放電調(diào)度也考慮了電壓約束和網(wǎng)損優(yōu)化。如果你是剛接觸這方面的人我建議不要一上來就啃數(shù)學(xué)論文先把“這道題在解決什么”想明白后面建模會順利很多。1.1 為什么偏偏是IEEE33和PG69IEEE33節(jié)點系統(tǒng)是配電網(wǎng)文獻里出現(xiàn)頻率最高的算例之一33個節(jié)點、32條支路、額定電壓12.66kV總負荷大約3.715MW加2.3Mvar。規(guī)模不上不下既不會讓初學(xué)的人看到一大堆節(jié)點發(fā)懵又能體現(xiàn)電壓損耗、線路重載、末端壓降這些真實配電網(wǎng)問題。更重要的是它的數(shù)據(jù)在網(wǎng)上很容易找到幾乎所有配電網(wǎng)重構(gòu)、分布式電源接入、儲能調(diào)度的文章都用過它你拿自己的模型去跑一遍跟論文里的結(jié)果一對就知道有沒有做錯。PG69這個系統(tǒng)更刺激一點69個節(jié)點比IEEE33多了一倍多總負荷大概3.8MW加2.7Mvar。它的干線更長越靠近末端電壓下降越明顯所以對“潮流計算電壓約束”來說是個更嚴格的測試場景。很多人在IEEE33上跑得通換到PG69就出現(xiàn)收斂問題或者錐松弛不精確恰恰就是因為系統(tǒng)規(guī)模變大、末端節(jié)點電壓過低讓模型對約束的真實壓力顯現(xiàn)出來。我把兩個系統(tǒng)的特點拉成一張表方便對比算例節(jié)點數(shù)支路數(shù)額定電壓典型總負荷特點IEEE33333212.66kV3.715MW 2.3Mvar規(guī)模適中資料豐富適合入門驗證PG69696812.66kV3.8MW 2.7Mvar干線路長末端壓降明顯適合壓力測試1.2 “多時間斷面”到底多在哪單時段潮流優(yōu)化其實很簡單就是給定某一時刻的負荷和DG出力算一個最優(yōu)運行點。但現(xiàn)實里光伏出力從早到晚是條曲線負荷也是波動的儲能更是晚上充電、白天放電這天然是一個跨時段耦合的問題。所謂多時間斷面就是你把一天按小時切成24個斷面甚至按15分鐘一個斷面切成96個斷面所有斷面共享同一個網(wǎng)絡(luò)拓撲每個斷面有自己的負荷、新能源出力和電壓分布但斷面與斷面之間又通過儲能SOC、機組爬坡、需求響應(yīng)等約束串在一起。普通的潮流計算不用管這些但最優(yōu)潮流如果忽略時間耦合它給出的儲能策略一定是錯的。比如儲能單時段模型里它只需要滿足“這一小時充的電等于這一小時放的電”這明顯不符合實際。多時段模型里你要約束每兩個相鄰時段的電量遞推關(guān)系還要限制充放電狀態(tài)不能同時為1這些都是時間維帶來的麻煩。這也是我在IEEE33和PG69上做多時段建模時最深的體會空間維度靠DistFlow方程時間維度靠儲能和爬坡兩者加在一起問題才完整。1.3 這套工具組合的邏輯YALMIP守前場CPLEX收尾工具選型很多人問過為什么不用Gurobi不用Mosek偏偏用CPLEX我的理由很簡單首先很多高校和研究所本來就有CPLEX的學(xué)術(shù)許可證裝起來不費勁其次CPLEX對混合整數(shù)二階錐規(guī)劃也就是MISOCP支持非常成熟第三YALMIP作為MATLAB里的建模層讓你用近乎數(shù)學(xué)公式的方式寫出錐約束和整數(shù)變量避免了手動把問題轉(zhuǎn)成求解器接口的麻煩。流程是這樣的YALMIP負責(zé)把二階錐約束、儲能SOC約束、整數(shù)變量這些建模邏輯轉(zhuǎn)化成CPLEX能識別的標準形式CPLEX內(nèi)部用分支切割和內(nèi)點法把MISOCP解出來。整個過程你不需要手寫任何內(nèi)點法代碼也基本不用關(guān)心算法的迭代細節(jié)。這也是工程上最舒服的地方建模思路集中在“問題長什么樣”而不是“求解器內(nèi)部怎么做”。2. 二階錐松弛建模中最關(guān)鍵的一跳2.1 DistFlow方程的非線性卡點配電網(wǎng)潮流計算如果從零開始寫最經(jīng)典的做法是用DistFlow方程。以一條支路ij為例從節(jié)點i流向節(jié)點j的有功功率記作P_ij無功功率記作Q_ij節(jié)點電壓幅值的平方記作v_i支路電流幅值的平方記作l_ij。DistFlow的核心方程可以寫成這幾條P_ij p_j^load - p_j^dg Σ P_jk r_ij * l_ijQ_ij q_j^load - q_j^dg Σ Q_jk x_ij * l_ijv_j v_i - 2 (r_ij P_ij x_ij Q_ij) (r_ij2 x_ij2) l_ijl_ij (P_ij2 Q_ij2) / v_i前三條都是線性約束問題就出在最后一條。l_ij等于一個分式分母是v_i分子是P_ij和Q_ij的平方和這是一個非凸等式。你去做優(yōu)化的時候這個等式會讓整個可行域變得極度扭曲普通的線性規(guī)劃和二次規(guī)劃都沒法直接處理。很多初學(xué)的人會問那直接用牛頓法跑潮流不就行了嗎可以但你要搞清楚潮流計算和最優(yōu)潮流是兩回事。潮流計算是給定運行點求解狀態(tài)量你需要的是唯一可行解而最優(yōu)潮流是要在約束里找使目標最小化的那個點如果模型非凸你得到的很可能只是局部最優(yōu)甚至根本不收斂。要做全局最優(yōu)的分析就得想辦法把這團非凸的東西“掰”成凸問題。2.2 從非凸等式到旋轉(zhuǎn)錐約束二階錐松弛的思路很直接。原來那個等式約束太苛刻了我把它放寬成一個不等式l_ij ≥ (P_ij2 Q_ij2) / v_i這個不等式的意思是支路電流的平方至少要達到歐姆定律要求的最小值。因為目標函數(shù)里通常會包含網(wǎng)損項求解器為了壓低損耗會盡量讓l_ij變小最后它會自動被壓到這個不等式邊界上也就是說最優(yōu)解仍然滿足等式。于是我們把“直接求等式”變成了“求不等式的最優(yōu)解”。這個不等式還帶一個分式直接交給CPLEX也不方便所以要再轉(zhuǎn)寫成標準的二階錐形式。經(jīng)過數(shù)學(xué)變形可以得到下面這個約束|| [2P_ij; 2Q_ij; v_i - l_ij] ||? ≤ v_i l_ij這看起來唬人其實很簡單。把兩邊同時平方展開左邊是4P_ij2加4Q_ij2再加(v_i - l_ij)2右邊是(v_i l_ij)2化簡之后得到的就是4P_ij2 4Q_ij2 ≤ 4 v_i l_ij等價于v_i l_ij ≥ P_ij2 Q_ij2正好是我們想要的松弛不等式。在YALMIP里寫這個約束不需要你自己去展開矩陣直接用cone函數(shù)就行Cons [Cons, cone([2*Pij(b,t); 2*Qij(b,t); V2(i,t)-Iij(b,t)], V2(i,t)Iij(b,t))];cone函數(shù)第一個參數(shù)是一個向量第二個參數(shù)是標量r它建模的數(shù)學(xué)含義就是||向量||? ≤ r。YALMIP會自動判斷這是一個二階錐約束并在傳給CPLEX之前把它轉(zhuǎn)換成標準SOCP格式。這一步是整套模型里最關(guān)鍵也最容易被忽略的地方。2.3 什么時候二階錐松弛是精確的松弛之后可行域變大了不等于最優(yōu)解一定落在原可行域上。如果你做完優(yōu)化發(fā)現(xiàn)l_ij明顯大于(P_ij2Q_ij2)/v_i說明錐松弛“不緊”解出來的結(jié)果并不對應(yīng)一個真實物理潮流。根據(jù)理論和實際經(jīng)驗在IEEE33和PG69這類輻射狀配電網(wǎng)里只要滿足兩個條件錐松弛基本是精確的第一每條支路的電阻和電抗都大于0第二目標函數(shù)對支路電流平方l_ij是嚴格單調(diào)遞增或者至少不是反向激勵。典型的目標比如最小化網(wǎng)損、最小化購電成本這些目標都會促使求解器把電流壓到最低最后錐約束自然收緊。實際操作中我推薦在目標里顯式加入網(wǎng)損項哪怕你的核心目標是調(diào)度儲能或削減峰值也放一個很小的網(wǎng)損懲罰系數(shù)。這樣既不會明顯扭曲經(jīng)濟性目標又能保證錐松弛的精確性。跑完之后檢查一下錐松弛間隙也是一個必要的自檢步驟。2.4 加入離散動作后問題升級為MISOCP多時段時間斷面如果只做連續(xù)變量純SOCP就夠了。但只要你想考慮有載調(diào)壓變壓器分接頭、電容器組投切、儲能充放電狀態(tài)問題就變成混合整數(shù)二階錐規(guī)劃這些離散量必須用整數(shù)變量或者0-1變量建模。CPLEX處理MISOCP的能力相當(dāng)強這也是我堅持用它的原因之一。在YALMIP里定義一個儲能充放電狀態(tài)變量只需要寫u binvar(ns, T); % ns個儲能T個時段后面再配合兩個不等式保證充電功率大于0時放電功率為0Pch Pmax * uPdis Pmax * (1 - u)如果沒有這個0-1變量模型很可能出現(xiàn)同一時段又充電又放電的荒謬結(jié)果。這個細節(jié)雖然簡單但很多人第一次建模時都會漏掉導(dǎo)致結(jié)果完全沒法看。3. 多時間斷面的完整數(shù)學(xué)建模和變量組織3.1 時間耦合約束到底從哪來多時段問題里最容易寫錯的就是時間耦合約束其中最典型的是儲能SOC遞推關(guān)系。以一個時間段間隔為Δt的模型為例假設(shè)Δt以小時為單位儲能荷電狀態(tài)E_t的計算公式是E_{t1} E_t η_ch * Pch_t * Δt - Pdis_t * Δt / η_dis其中η_ch是充電效率η_dis是放電效率。這個遞推式把相鄰兩個時段綁在一起不能獨立求解。與此同時還要讓充放電功率不同時為正所以引入了前面說的0-1變量u_t。除了儲能分布式電源出力也有爬坡速率限制-ΔP_ramp ≤ Pg_{t1} - Pg_t ≤ ΔP_ramp這個約束本質(zhì)上也是時間耦合只不過比儲能遞推稍微溫和一點它不限制電量只限制每一段之間的變化幅度。如果模型還包含有載調(diào)壓變壓器、電容器組那么這些設(shè)備的分接頭檔位和投切動作也應(yīng)該跨時段限制。最簡單的做法是兩個相鄰時段最多動作一次或者限制一天內(nèi)動作總次數(shù)。這些約束不加求解器就會給出每15分鐘瘋狂切換一次的理想方案工程上根本沒法執(zhí)行。3.2 每個時段都要滿足的DistFlow約束時間維度增加了但每個時段內(nèi)部仍然要滿足配電網(wǎng)潮流方程。對于每個斷面t我要對每一條支路b寫一組DistFlow約束。先給支路定義from和to兩個方向數(shù)組from br(:, 1); % 支路起點編號 to br(:, 2); % 支路終點編號 R_ohm br(:, 3); % 電阻單位歐姆 X_ohm br(:, 4); % 電抗單位歐姆然后利用前面說過的SOCP錐約束對每個時間斷面建立支路電壓降和錐約束。代碼如下% 電壓降線性約束 Cons [Cons, V2(to(b), t) V2(from(b), t) ... - 2*(R_pu(b)*Pij(b,t) X_pu(b)*Qij(b,t)) ... (R_pu(b)^2 X_pu(b)^2)*Iij(b,t)]; % 二階錐約束 Cons [Cons, cone([2*Pij(b,t); 2*Qij(b,t); ... V2(from(b),t) - Iij(b,t)], V2(from(b),t) Iij(b,t))];節(jié)點功率平衡也要逐時段寫。對于每個節(jié)點i所有流出支路功率減去所有流入支路功率再加上線路損耗必須等于該節(jié)點凈注入功率。這里的凈注入等于DG出力減去負荷公式為Σ Pij_{out} - Σ Pij_{in} Σ r_b * l_ij_b Pg_i - Pload_i實際寫代碼的時候用find函數(shù)找每個節(jié)點的出線和進線支路索引就行。3.3 目標函數(shù)怎么設(shè)計才算合理目標函數(shù)是整個模型戰(zhàn)略性的部分。我做IEEE33和PG69多時段算例時最常使用的是三部分相加購電成本、網(wǎng)損、儲能成本。購電成本從上級電網(wǎng)買電的費用通常用分時電價乘以根節(jié)點注入功率網(wǎng)損所有支路電流平方乘以支路電阻再求和這個也是讓SOCP約束保持緊湊的關(guān)鍵項儲能成本可以理解成電池循環(huán)損耗也可以加一個對充電次數(shù)的軟約束。YALMIP里目標函數(shù)可以直接寫成Objective sum(price .* Pg) ... sum(sum(R_pu .* Iij)) * baseMVA ... 0.001 * sum(sum(Pch Pdis));有人會問為什么網(wǎng)損要用p.u.值再乘baseMVA因為配電網(wǎng)里潮流功率的真實量級是MW而線路電阻在p.u.下通常只有0.005左右兩者乘出來的網(wǎng)損可能在1e-4這種量級直接作為目標會被求解器當(dāng)成噪聲忽略反而丟失了錐松弛的收緊作用。乘回baseMVA后網(wǎng)損就恢復(fù)成幾十千瓦到幾百千瓦的真實量級和購電成本保持在同一個數(shù)量級CPLEX在數(shù)值上才會認真優(yōu)化它。3.4 變量的維度設(shè)計是少走彎路的重點多時段模型的變量維度設(shè)計我建議一開始就按“支路數(shù)×?xí)r段數(shù)”或“節(jié)點數(shù)×?xí)r段數(shù)”來建矩陣而不是每個時段單獨設(shè)一套變量。這樣YALMIP內(nèi)部矩陣維數(shù)小求解速度也快。以IEEE33為例T取24那么Pij sdpvar(nb, T); % 支路有功nb行T列 Qij sdpvar(nb, T); % 支路無功 Iij sdpvar(nb, T); % 支路電流平方 V2 sdpvar(n, T); % 節(jié)點電壓平方 Pg sdpvar(ng, T); % 上級電網(wǎng)注入有功 Qg sdpvar(ng, T); % 上級電網(wǎng)注入無功后面在寫約束時只要把t從1到T循環(huán)一遍用Pij(b,t)索引某一個具體時間段即可。這種二維變量矩陣既直觀又能利用YALMIP對結(jié)構(gòu)化變量的優(yōu)化比把所有變量拉成一維長向量再慢慢拼約束要舒服得多。4. CPLEXYALMIP實戰(zhàn)落地從環(huán)境配置到代碼骨架4.1 安裝環(huán)節(jié)最容易翻車先講清楚很多人的模型本身沒有錯但卡在第一步CPLEX裝好了YALMIP卻找不到求解器。YALMIP不是自帶求解器的它只是一個建模層它把所有約束翻譯成標準模型文件后需要調(diào)用外部求解器來算。所以必須保證CPLEX的MATLAB接口路徑被正確添加。CPLEX安裝好后注意找到它的MATLAB接口目錄一般是安裝路徑下的cplex/matlab文件夾。在MATLAB里執(zhí)行addpath(genpath(D:\Program Files\IBM\ILOG\CPLEX_Studio221\cplex\matlab)); savepathYALMIP的安裝也是類似把下載下來的文件夾整個加入路徑addpath(genpath(D:\yalmip)); savepath裝完之后運行yalmiptest你會看到Y(jié)ALMIP檢查所有已安裝求解器的報告。里面會有一行顯示CPLEX相關(guān)的測試是否通過。如果顯示missing或者error基本就是路徑?jīng)]配對。關(guān)于CPLEX獲取方式高校用戶建議走學(xué)術(shù)計劃申請教育版證書個人使用也可以看看官方社區(qū)里免費的社區(qū)版不過社區(qū)版對問題規(guī)模有限制IEEE33的24時段問題勉強能頂PG69的多時段MISOCP很可能超限長期做研究還是用完整許可證靠譜。4.2 數(shù)據(jù)準備和標幺化IEEE33和PG69的原始數(shù)據(jù)里面線路參數(shù)通常是以歐姆為單位給出的負荷以kW為單位而優(yōu)化求解器在處理SOCP這種含有二次約束的問題時數(shù)值范圍太離譜會導(dǎo)致收斂困難。我強烈建議計算前先把所有量統(tǒng)一到標幺值下。一個實用的標幺選擇是基準電壓Vb取12.66kV基準功率Sb取10MVA。這樣IEEE33的總負荷3.715MW就變成0.3715p.u.數(shù)值合理。線路阻抗的轉(zhuǎn)換公式是Z_b Vb2 / Sbr_pu R_ohm / Z_bx_pu X_ohm / Z_b打個具體比方IEEE33第一條支路的電阻是0.0922ΩVb12.66kVSb10MVA時Z_b等于16.02Ω換算出來的標幺電阻大約0.00576。這個數(shù)值在二階錐約束中不會過小目標函數(shù)中的網(wǎng)損項也不會被淹沒。如果不做這一步直接用歐姆值和kW去建模數(shù)學(xué)上雖可行但數(shù)值條件數(shù)會很差CPLEX求解效率和精度都會下降。4.3 一個可直接套用的求解代碼骨架我習(xí)慣把建模代碼拆成幾個邏輯塊數(shù)據(jù)定義、變量聲明、約束循環(huán)、目標函數(shù)、求解、結(jié)果提取。下面是一個基于IEEE33、24時段、含儲能的多時段SOCP最簡骨架%% 1. 基礎(chǔ)數(shù)據(jù) T 24; n 33; nb 32; baseMVA 10; % branch數(shù)據(jù)假設(shè)已經(jīng)整理成from, to, R_ohm, X_ohm Zb 12.66^2 / baseMVA; R_pu R_ohm / Zb; X_pu X_ohm / Zb; %% 2. 定義變量 Pij sdpvar(nb, T); Qij sdpvar(nb, T); Iij sdpvar(nb, T); V2 sdpvar(n, T); Pg sdpvar(1, T); % 根節(jié)點注入有功 Qg sdpvar(1, T); %% 3. 儲能變量 ns 2; Pch sdpvar(ns, T); Pdis sdpvar(ns, T); E sdpvar(ns, T); u binvar(ns, T); %% 4. 約束 Cons []; for t 1:T % 根節(jié)點電壓設(shè)為1.0 Cons [Cons, V2(1,t) 1.0]; % 支路DistFlow與錐約束 for b 1:nb i from(b); j to(b); Cons [Cons, V2(j,t) V2(i,t) ... - 2*(R_pu(b)*Pij(b,t) X_pu(b)*Qij(b,t)) ... (R_pu(b)^2 X_pu(b)^2)*Iij(b,t)]; Cons [Cons, cone([2*Pij(b,t); 2*Qij(b,t); V2(i,t)-Iij(b,t)], ... V2(i,t)Iij(b,t))]; end % 節(jié)點有功平衡 for i 1:n outb find(from i); inb find(to i); Cons [Cons, sum(Pij(outb,t)) - sum(Pij(inb,t)) ... sum(R_pu(outb).*Iij(outb,t)) Pg(1,t) - Pload(i,t)]; end % 節(jié)點無功平衡格式類似這里省略 end %% 5. 儲能時間耦合 for s 1:ns for t 1:T-1 Cons [Cons, E(s,t1) E(s,t) 0.9*Pch(s,t) - Pdis(s,t)/0.9]; Cons [Cons, Pch(s,t) 0.5*u(s,t) ... , Pdis(s,t) 0.5*(1-u(s,t))]; end end %% 6. 目標函數(shù) Objective sum(Pg(1,:) .* price) ... sum(sum(R_pu .* Iij)) * baseMVA ... 0.001*sum(sum(Pch Pdis)); %% 7. 求解 ops sdpsettings(solver,cplex,verbose,2,showprogress,1); ops.cplex.mip.tolerances.mipgap 1e-4; optimize(Cons, Objective, ops);這段代碼的核心思路就是先把網(wǎng)絡(luò)約束逐時段寫入再補時間耦合約束。實際運行時要根據(jù)你自己數(shù)據(jù)里的Pload(i,t)和price(t)去完善我這里為了讓骨架更干凈省略了無功平衡的對應(yīng)寫法你做的時候別漏。4.4 求解器參數(shù)千萬不要全默認很多人習(xí)慣直接optimize(Cons, Objective)在簡單小問題上沒問題但多時段的MISOCP尤其PG69這種上百個支路的算例不做參數(shù)調(diào)整會讓求解時間從幾十秒變成幾小時。我常用的CPLEX關(guān)鍵參數(shù)有三個ops.cplex.mip.tolerances.mipgap相對MIP間隙。默認可能是1e-4甚至更小對于工程分析我一般放1e-4。如果你只需要評估策略趨勢放到1e-3能大幅度加速。ops.cplex.timelimit求解時間上限。實際工程中設(shè)一個時間限制比如600秒超時后直接取當(dāng)前最好可行解避免無休止地卡在分支定界上。ops.cplex.threads并行線程數(shù)。CPLEX默認會用滿所有核心但有時候并行反而導(dǎo)致內(nèi)存爆炸手動調(diào)到物理核心數(shù)會穩(wěn)一點。另外verbose開到2可以在命令行看到實時間隙變化對判斷求解是否卡住很有幫助。我就是靠這個判斷是加約束還是調(diào)參數(shù)不用等半天最后看一個冷冰冰的infeasible。5. 我踩過的坑和排查手冊5.1 求解器報“infeasible”怎么辦新模型的第一個報錯大多是不可行。我見過太多人直接開始改動約束但我習(xí)慣反過來做減法檢查。第一步把儲能時間耦合約束全部注釋掉只算每個時段獨立的潮流如果這時候還是不可行說明問題出在網(wǎng)絡(luò)約束本身跟時間維度無關(guān)。第二步檢查負荷是不是有個別節(jié)點沒接上導(dǎo)致功率不平衡第三步檢查電壓上下限約束是不是設(shè)得太緊比如把0.93到1.07改成0.95到1.05可能某個時段本身就解不出來。還有一種隱蔽情況我一開始用SB1MVA做標幺IEEE33的負荷達到3.7p.u.儲能的充電功率0.5p.u.變成了0.5MW看似合理但電壓降方程的數(shù)值尺度還是有點別扭。后來把SB改成10MVA問題立刻順了很多。所以infeasible排查的第一步永遠是回頭看標幺基準別急著懷疑模型邏輯。5.2 怎么驗證錐松弛到底緊不緊SOCP結(jié)果出來之后不要只盯著目標函數(shù)值就完事。我會單獨算一下每一條支路的錐松弛間隙powerij value(Pij); powerqj value(Qij); current value(Iij); volt2 value(V2); % 對每條支路的每個時段計算松弛間隙 for t 1:T for b 1:nb v_i volt2(from(b), t); gap(b,t) v_i * current(b,t) - powerij(b,t)^2 - powerqj(b,t)^2; end end max_gap max(gap(:));如果max_gap在1e-6這個量級甚至更低說明錐松弛是緊的解出來的結(jié)果可以直接當(dāng)真實潮流用。如果gap跑到1e-3或者更大說明松弛不緊你得回去檢查目標函數(shù)里是不是漏了網(wǎng)損項。這個方法花不了幾行代碼但能讓結(jié)果可信度提高一大截我在IEEE33上跑多時段時一般gap都能壓到1e-6以下。5.3 求解時間太長怎么辦第一次面對PG69的96時段、帶儲能的模型時我差點懷疑電腦壞了連續(xù)跑了兩個小時都沒出結(jié)果。后來發(fā)現(xiàn)問題是三個疊加整數(shù)變量太多、MIP gap設(shè)太小、verbose都沒開看不出來進展。解決手段有三個方向。第一個方向是變量瘦身儲能充放電狀態(tài)雖然有0-1變量但如果允許儲能長時間不動作就不要給每個存儲電池都加狀態(tài)變量合理簡化能砍掉一半整數(shù)變量。第二個方向是目標函數(shù)里的懲罰系數(shù)不要設(shè)置太高否則會造成數(shù)值上的瓶頸求解器會在同一批次上反復(fù)切割。第三個方向是時間步長策略先跑T6驗證模型再逐步放大T12、T24這樣能快速定位是哪類約束導(dǎo)致復(fù)雜度暴漲。等T24模型能在幾十秒內(nèi)跑完再考慮加更多時段。5.4 數(shù)值警告滿天飛的排查思路CPLEX有時候會給出類似Overflow或者Numerical difficulties的警告。我在多時段建模中遇到過的常見誘因是約束兩邊數(shù)值差距過大。比如線路電阻標幺值是1e-5而功率變量是幾十乘完之后約束右側(cè)的量級跨度介于幾個數(shù)量級之間自然容易出問題。解決思路還是回到標幺化和約束重縮放上。你可以把基準功率調(diào)大比如1MVA換成10MVA或者把目標函數(shù)中的權(quán)重縮放一下讓所有約束等式右側(cè)都在0.001到1000之間數(shù)值困難就很少出現(xiàn)了。6. 最后分享一個我堅持至今的調(diào)試習(xí)慣這套多時段二階錐松弛模型我前前后后寫過不下十個版本每一次改動完后我都會先固定T2或者T3只測極短時間斷面。這樣做的好處顯而易見兩個小時跑不出來的問題壓縮成兩個時段可能十秒就能出結(jié)果我可以迅速驗證約束是否寫錯、錐約束是否被正確識別、數(shù)值是否有警告。短時段的解雖然經(jīng)濟性意義不大但它的結(jié)構(gòu)足以暴露建模錯誤比一上來就跑96時段然后對著日志發(fā)呆要高效得多。另外一個習(xí)慣是結(jié)果可視化。多時段優(yōu)化跑完把電壓分布畫成二維熱力圖或者各時段曲線很多時候一眼就能看出問題。比如某條支路電流階段處出現(xiàn)明顯尖峰很可能就是約束里哪個節(jié)點索引對不上。數(shù)值上看著正常的解畫出來未必正常這個經(jīng)驗我在IEEE33和PG69上都驗證過無數(shù)次。你如果也在這兩個算例上做多時段建模我建議早點養(yǎng)成這兩個習(xí)慣能幫你少熬好幾個通宵。