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