二階錐規(guī)劃的配電網(wǎng)重構(gòu)：模型、實(shí)現(xiàn)與優(yōu)化)
1. 項(xiàng)目概述當(dāng)配電網(wǎng)重構(gòu)遇上混合整數(shù)二階錐在電力系統(tǒng)運(yùn)行領(lǐng)域配電網(wǎng)重構(gòu)是一個(gè)經(jīng)典且核心的優(yōu)化問題。簡單來說它就像是在一個(gè)復(fù)雜的、由眾多開關(guān)分段開關(guān)和聯(lián)絡(luò)開關(guān)組成的城市道路網(wǎng)絡(luò)中通過改變某些路段的通斷狀態(tài)來重新規(guī)劃電力流動(dòng)的路徑。這么做的目標(biāo)很明確降低整個(gè)網(wǎng)絡(luò)的電能損耗、改善電壓質(zhì)量、平衡各條線路的負(fù)載或者在故障后快速恢復(fù)供電。傳統(tǒng)上這個(gè)問題被建模為一個(gè)混合整數(shù)非線性規(guī)劃問題因?yàn)殚_關(guān)狀態(tài)是0或1的整數(shù)變量而網(wǎng)絡(luò)中的功率流方程則是非線性的。非線性尤其是非凸性是求解這類問題的“攔路虎”它意味著你可能找到一個(gè)局部最優(yōu)解但無法保證它是全局最優(yōu)的而且求解過程往往非常耗時(shí)?；旌险麛?shù)二階錐規(guī)劃的出現(xiàn)為這個(gè)難題提供了一個(gè)極具吸引力的突破口。它的核心思想是在一定的假設(shè)和近似下比如采用DistFlow潮流模型可以將原本非線性的功率平衡方程轉(zhuǎn)化為一系列二階錐約束。二階錐約束是一種特殊的凸約束其形式類似于一個(gè)“冰淇淋蛋筒”的形狀。將問題轉(zhuǎn)化為MISOCP模型后我們實(shí)際上是在一個(gè)凸的可行域內(nèi)尋找最優(yōu)的整數(shù)解。凸性保證了局部最優(yōu)即全局最優(yōu)并且現(xiàn)代的商業(yè)求解器如CPLEX、Gurobi對這類問題有非常高效的求解算法。因此基于MISOCP的配電網(wǎng)重構(gòu)本質(zhì)上是用一種可高效求解的凸優(yōu)化框架去逼近和求解原本棘手的非凸問題在計(jì)算效率和求解質(zhì)量之間取得了出色的平衡。對于電力系統(tǒng)研究人員和工程師而言掌握這套方法意味著你手里多了一把解決復(fù)雜網(wǎng)絡(luò)優(yōu)化問題的“瑞士軍刀”。2. 核心模型構(gòu)建從物理方程到數(shù)學(xué)規(guī)劃要理解MISOCP模型我們必須從配電網(wǎng)的物理基礎(chǔ)——潮流模型開始。在配電網(wǎng)重構(gòu)中我們通常采用基于支路潮流的DistFlow模型它比傳統(tǒng)的節(jié)點(diǎn)功率方程更直觀地描述了功率在輻射狀網(wǎng)絡(luò)中的流動(dòng)。2.1 DistFlow潮流方程及其凸松弛考慮一個(gè)由節(jié)點(diǎn)母線和支路線路組成的配電網(wǎng)。對于任意支路(i, j)設(shè)其首端節(jié)點(diǎn)i的電壓幅值為V_i末端節(jié)點(diǎn)j的電壓幅值為V_j支路上流過的有功功率為P_ij無功功率為Q_ij支路電阻為r_ij電抗為x_ij。經(jīng)典的DistFlow方程可以寫為P_ij - r_ij * I_ij^2 Σ_{k∈δ(j)} P_jk p_j Q_ij - x_ij * I_ij^2 Σ_{k∈δ(j)} Q_jk q_j V_j^2 V_i^2 - 2*(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*I_ij^2其中p_j和q_j是節(jié)點(diǎn)j的凈注入有功和無功功率負(fù)荷為負(fù)分布式電源為正δ(j)表示以j為首端的所有下游支路集合I_ij是支路電流幅值。第三個(gè)方程中包含了電流平方項(xiàng)I_ij^2它與功率的關(guān)系是非線性的I_ij^2 (P_ij^2 Q_ij^2) / V_i^2。正是這個(gè)項(xiàng)引入了非凸性。為了進(jìn)行凸松弛我們引入一組新的變量進(jìn)行替換令 u_i V_i^2代表節(jié)點(diǎn)電壓的平方。令 l_ij I_ij^2代表支路電流的平方。支路功率P_ij和Q_ij保持不變。經(jīng)過變量替換和忽略高階項(xiàng)(r_ij^2 x_ij^2)*l_ij因其值通常很小DistFlow方程可以重寫為P_ij - r_ij * l_ij Σ_{k∈δ(j)} P_jk p_j Q_ij - x_ij * l_ij Σ_{k∈δ(j)} Q_jk q_j u_j u_i - 2*(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*l_ij此時(shí)非凸性隱藏在變量之間的關(guān)系中根據(jù)歐姆定律和功率定義應(yīng)有 P_ij^2 Q_ij^2 u_i * l_ij。這是一個(gè)旋轉(zhuǎn)錐約束是非凸的。二階錐松弛的關(guān)鍵一步來了我們將這個(gè)等式約束松弛為一個(gè)不等式約束P_ij^2 Q_ij^2 ≤ u_i * l_ij這個(gè)不等式可以等價(jià)地寫成一個(gè)標(biāo)準(zhǔn)二階錐約束的形式|| [2P_ij, 2Q_ij, u_i - l_ij]^T ||_2 ≤ u_i l_ij其中||·||_2表示向量的二范數(shù)。這個(gè)約束定義了一個(gè)凸集二階錐。在輻射狀配電網(wǎng)且負(fù)荷水平不是極端輕載的情況下這個(gè)松弛通常是緊的即最優(yōu)解會(huì)自動(dòng)使不等式取等號(hào)從而滿足原始物理方程。這就成功地將非凸問題嵌入了凸框架內(nèi)。2.2 網(wǎng)絡(luò)拓?fù)渑c開關(guān)建模配電網(wǎng)重構(gòu)的核心是改變網(wǎng)絡(luò)拓?fù)溥@通過改變支路的連接狀態(tài)即開關(guān)的開合來實(shí)現(xiàn)。我們需要引入0-1整數(shù)變量來建模這種狀態(tài)。對于網(wǎng)絡(luò)中的每一條可能的支路包括現(xiàn)有的常閉支路和作為備用的常開支路我們定義一個(gè)二進(jìn)制變量z_ij ∈ {0, 1}。z_ij 1表示支路(i, j)閉合通電z_ij 0表示支路(i, j)斷開。網(wǎng)絡(luò)拓?fù)湫枰獫M足兩個(gè)核心約束輻射狀約束重構(gòu)后的網(wǎng)絡(luò)必須是一個(gè)輻射狀樹形結(jié)構(gòu)即無環(huán)且連通。這可以通過多種方式建模最常見的是利用配電網(wǎng)絡(luò)單電源點(diǎn)的特點(diǎn)采用“虛擬流”或“父節(jié)點(diǎn)指示”方法。虛擬流法為每個(gè)節(jié)點(diǎn)除根節(jié)點(diǎn)外分配一個(gè)單位的虛擬需求根節(jié)點(diǎn)變電站為虛擬源。要求虛擬流只能從閉合的支路上通過并且每個(gè)非根節(jié)點(diǎn)必須有且僅有一條流入的閉合支路。這保證了連通性和無環(huán)性。父節(jié)點(diǎn)變量法為每個(gè)節(jié)點(diǎn)定義一個(gè)整數(shù)變量表示其父節(jié)點(diǎn)編號(hào)并約束根節(jié)點(diǎn)無父節(jié)點(diǎn)其他節(jié)點(diǎn)有且僅有一個(gè)父節(jié)點(diǎn)且父節(jié)點(diǎn)必須通過閉合的支路與之相連。這種方法更直觀但變量更多。開關(guān)邏輯約束支路的電氣狀態(tài)必須與開關(guān)狀態(tài)一致。當(dāng)z_ij0開關(guān)斷開時(shí)該支路上的功率P_ij, Q_ij必須為0電流平方l_ij也必須為0。這可以用大M法來建模-M * z_ij ≤ P_ij ≤ M * z_ij -M * z_ij ≤ Q_ij ≤ M * z_ij 0 ≤ l_ij ≤ M * z_ij 0 ≤ u_i - u_j ≤ M * z_ij ΔV_max*(1-z_ij) // 可選斷開時(shí)電壓差約束可放松其中M是一個(gè)足夠大的正數(shù)但為了數(shù)值穩(wěn)定性應(yīng)盡可能選取緊的界。2.3 目標(biāo)函數(shù)與完整MISOCP模型配電網(wǎng)重構(gòu)的典型目標(biāo)是最小化系統(tǒng)總有功網(wǎng)損。網(wǎng)損等于所有支路電阻上的損耗之和Ploss Σ_{(i,j)} r_ij * l_ij。由于l_ij是電流的平方這個(gè)目標(biāo)函數(shù)是線性的。綜合以上所有部分我們可以得到完整的基于MISOCP的配電網(wǎng)重構(gòu)模型目標(biāo)函數(shù)Minimize Σ_{(i,j)∈Ω} r_ij * l_ij約束條件節(jié)點(diǎn)功率平衡線性等式 Σ_{i∈π(j)} P_ij - Σ_{k∈δ(j)} P_jk - r_ij * l_ij p_j, ?j Σ_{i∈π(j)} Q_ij - Σ_{k∈δ(j)} Q_jk - x_ij * l_ij q_j, ?j π(j)表示節(jié)點(diǎn)j的所有父節(jié)點(diǎn)支路集合電壓降落方程線性等式 u_j u_i - 2*(r_ijP_ij x_ijQ_ij) (r_ij^2 x_ij^2)*l_ij, ?(i,j)∈Ω二階錐松弛約束凸錐約束 || [2P_ij, 2Q_ij, u_i - l_ij]^T ||_2 ≤ u_i l_ij, ?(i,j)∈Ω開關(guān)邏輯約束大M法混合整數(shù)線性 -M * z_ij ≤ P_ij ≤ M * z_ij, ?(i,j)∈Ω -M * z_ij ≤ Q_ij ≤ M * z_ij, ?(i,j)∈Ω 0 ≤ l_ij ≤ M * z_ij, ?(i,j)∈Ω輻射狀拓?fù)浼s束混合整數(shù)線性以虛擬流法為例 Σ_{i∈π(j)} f_ij 1, ?j ≠ 根節(jié)點(diǎn) 0 ≤ f_ij ≤ N * z_ij, ?(i,j)∈Ω // f_ij為從i到j(luò)的虛擬流N為節(jié)點(diǎn)數(shù) Σ_{(i,j)∈Ω} z_ij N - 1 // 輻射狀網(wǎng)絡(luò)支路數(shù)等于節(jié)點(diǎn)數(shù)減一運(yùn)行安全約束線性不等式 u_min^2 ≤ u_i ≤ u_max^2, ?i // 電壓幅值平方約束 l_ij ≤ (I_ij_max)^2 * z_ij, ?(i,j)∈Ω // 支路電流容量約束整數(shù)變量約束 z_ij ∈ {0, 1}, ?(i,j)∈Ω這個(gè)模型就是一個(gè)標(biāo)準(zhǔn)的混合整數(shù)二階錐規(guī)劃問題。目標(biāo)函數(shù)和絕大多數(shù)約束是線性的關(guān)鍵的非線性部分被包含在凸的二階錐約束中因此可以被高效的商業(yè)求解器處理。3. 基于MATLAB與YALMIP/CPLEX的代碼實(shí)現(xiàn)詳解理論模型建立后我們需要一個(gè)強(qiáng)大的工具鏈將其轉(zhuǎn)化為可執(zhí)行的代碼。MATLAB因其強(qiáng)大的數(shù)學(xué)計(jì)算和矩陣操作能力成為算法原型開發(fā)的理想環(huán)境。YALMIP是一個(gè)在MATLAB中用于建模優(yōu)化問題的免費(fèi)工具箱它提供了一種非常直觀的方式來描述優(yōu)化問題變量、目標(biāo)、約束然后自動(dòng)調(diào)用后端求解器如CPLEX、Gurobi、MOSEK來求解。CPLEX是IBM推出的高性能數(shù)學(xué)規(guī)劃求解器對MISOCP有非常好的支持。3.1 開發(fā)環(huán)境搭建與數(shù)據(jù)準(zhǔn)備首先確保你的MATLAB環(huán)境已經(jīng)安裝了YALMIP和CPLEX。YALMIP可以直接從其官網(wǎng)下載并添加到MATLAB路徑。CPLEX需要從IBM官網(wǎng)獲取學(xué)術(shù)版或商業(yè)版許可并安裝。安裝后在MATLAB中運(yùn)行yalmiptestYALMIP會(huì)自動(dòng)檢測可用的求解器確認(rèn)CPLEX被正確識(shí)別。數(shù)據(jù)準(zhǔn)備是第一步也是最容易出錯(cuò)的一步。我們需要一個(gè)標(biāo)準(zhǔn)的測試配電網(wǎng)數(shù)據(jù)例如IEEE 33節(jié)點(diǎn)、69節(jié)點(diǎn)或118節(jié)點(diǎn)系統(tǒng)。數(shù)據(jù)通常包括bus矩陣節(jié)點(diǎn)編號(hào)、類型平衡節(jié)點(diǎn)/負(fù)荷節(jié)點(diǎn)、有功負(fù)荷、無功負(fù)荷。branch矩陣支路首末端節(jié)點(diǎn)編號(hào)、電阻、電抗、額定電流、初始狀態(tài)0開/1合。baseMVA和baseKV系統(tǒng)基準(zhǔn)值用于標(biāo)幺化。我強(qiáng)烈建議將原始數(shù)據(jù)如.mat文件或Excel表格讀入后先進(jìn)行標(biāo)幺化處理。將所有阻抗、功率、電壓、電流值除以相應(yīng)的基準(zhǔn)值。標(biāo)幺化能顯著改善模型的數(shù)值穩(wěn)定性避免因?qū)嶋H數(shù)據(jù)量綱差異過大如電阻是0.幾歐姆功率是幾兆瓦導(dǎo)致求解器出現(xiàn)數(shù)值問題。% 示例數(shù)據(jù)讀取與標(biāo)幺化 load(IEEE33.mat); % 假設(shè)數(shù)據(jù)已保存在此文件 baseMVA 10; % 基準(zhǔn)功率MVA baseKV 12.66; % 基準(zhǔn)電壓kV Zbase baseKV^2 / baseMVA; % 基準(zhǔn)阻抗 % 標(biāo)幺化支路電阻和電抗 branch(:, 3) branch(:, 3) / Zbase; % 電阻R branch(:, 4) branch(:, 4) / Zbase; % 電抗X % 標(biāo)幺化節(jié)點(diǎn)負(fù)荷 bus(:, 3) bus(:, 3) / baseMVA; % 有功負(fù)荷P bus(:, 4) bus(:, 4) / baseMVA; % 無功負(fù)荷Q % 標(biāo)幺化電壓和電流限值 Vmax_pu 1.05; Vmin_pu 0.95; Imax_pu branch(:, 5) / (baseMVA / (sqrt(3)*baseKV)); % 假設(shè)branch第5列是額定電流(A)3.2 使用YALMIP構(gòu)建MISOCP模型YALMIP建模的核心是三步定義變量、定義約束、定義目標(biāo)函數(shù)。代碼的清晰度至關(guān)重要。% 1. 定義變量 nb size(bus, 1); % 節(jié)點(diǎn)數(shù) nl size(branch, 1); % 支路數(shù)包括所有可能的支路 % 連續(xù)變量 P sdpvar(nl, 1); % 支路有功潮流 Q sdpvar(nl, 1); % 支路無功潮流 U sdpvar(nb, 1); % 節(jié)點(diǎn)電壓平方 I sdpvar(nl, 1); % 支路電流平方 % 整數(shù)變量二進(jìn)制 z binvar(nl, 1); % 支路開關(guān)狀態(tài)1閉合0斷開 % 虛擬流變量用于輻射狀約束 f sdpvar(nl, 1); % 虛擬流非負(fù)連續(xù)變量 % 2. 定義約束 Constraints []; % 2.1 節(jié)點(diǎn)功率平衡約束 % 構(gòu)建節(jié)點(diǎn)-支路關(guān)聯(lián)矩陣Anb x nlA(i,k)1表示支路k以i為首端-1表示以i為末端 A zeros(nb, nl); for k 1:nl from branch(k, 1); to branch(k, 2); A(from, k) 1; A(to, k) -1; end % 平衡節(jié)點(diǎn)假設(shè)為節(jié)點(diǎn)1電壓固定 Constraints [Constraints, U(1) 1.0]; % 標(biāo)幺值通常設(shè)為1.0 for i 1:nb % 流入節(jié)點(diǎn)i的凈功率 該節(jié)點(diǎn)負(fù)荷 if i 1 % 平衡節(jié)點(diǎn)是功率源其注入功率是自由的由優(yōu)化決定 % 實(shí)際上平衡節(jié)點(diǎn)的功率平衡方程通常不顯式添加因?yàn)樗鼤?huì)自動(dòng)滿足 % 我們只需處理負(fù)荷節(jié)點(diǎn) else % 對于負(fù)荷節(jié)點(diǎn)從關(guān)聯(lián)矩陣中找出與該節(jié)點(diǎn)相關(guān)的支路 idx_in find(A(i, :) 1); % 以i為首端的支路功率流出i idx_out find(A(i, :) -1); % 以i為末端的支路功率流入i % 功率平衡流入 - 流出 - 損耗 負(fù)荷 % 注意DistFlow模型中損耗(r*I)發(fā)生在支路末端 P_in sum(P(idx_out)); P_out sum(P(idx_in)); P_loss_sum sum(branch(idx_out, 3) .* I(idx_out)); % 所有流入支路的損耗之和 Constraints [Constraints, P_in - P_out - P_loss_sum bus(i, 3)]; % bus(i,3)是有功負(fù)荷 Q_in sum(Q(idx_out)); Q_out sum(Q(idx_in)); Q_loss_sum sum(branch(idx_out, 4) .* I(idx_out)); % 所有流入支路的無功損耗之和 Constraints [Constraints, Q_in - Q_out - Q_loss_sum bus(i, 4)]; % bus(i,4)是無功負(fù)荷 end end % 2.2 電壓降落方程約束 for k 1:nl from branch(k, 1); to branch(k, 2); r branch(k, 3); x branch(k, 4); Constraints [Constraints, U(to) U(from) - 2*(r*P(k) x*Q(k)) (r^2 x^2)*I(k)]; end % 2.3 二階錐松弛約束 for k 1:nl from branch(k, 1); Constraints [Constraints, cone([2*P(k); 2*Q(k); U(from)-I(k)], U(from)I(k))]; % YALMIP中的cone函數(shù)直接用于構(gòu)建二階錐約束||x|| y % 這里 x [2P; 2Q; U-I], y UI end % 2.4 開關(guān)邏輯約束大M法 M_pq 10; % 根據(jù)系統(tǒng)基準(zhǔn)功率估算的足夠大的數(shù)例如最大可能潮流的2倍 M_i 2; % 電流平方的上界例如(2*Imax)^2 M_u 0.5; % 電壓差的上界 for k 1:nl Constraints [Constraints, -M_pq * z(k) P(k) M_pq * z(k)]; Constraints [Constraints, -M_pq * z(k) Q(k) M_pq * z(k)]; Constraints [Constraints, 0 I(k) M_i * z(k)]; % 對于斷開支路可以放松電壓差約束 from branch(k, 1); to branch(k, 2); Constraints [Constraints, -M_u*(1-z(k)) U(from)-U(to) M_u*(1-z(k))]; end % 2.5 輻射狀拓?fù)浼s束虛擬流法 % 假設(shè)節(jié)點(diǎn)1是根節(jié)點(diǎn)平衡節(jié)點(diǎn) for j 2:nb % 找到所有以j為末端的支路索引 idx_to_j find(branch(:,2) j); Constraints [Constraints, sum(f(idx_to_j)) 1]; end for k 1:nl Constraints [Constraints, 0 f(k) nb * z(k)]; end Constraints [Constraints, sum(z) nb-1]; % 輻射狀網(wǎng)絡(luò)支路數(shù)約束 % 2.6 運(yùn)行安全約束 for i 1:nb Constraints [Constraints, Vmin_pu^2 U(i) Vmax_pu^2]; end for k 1:nl Constraints [Constraints, I(k) (Imax_pu(k))^2 * z(k)]; end % 3. 定義目標(biāo)函數(shù) Objective sum(branch(:,3) .* I); % 最小化總有功網(wǎng)損3.3 模型求解與結(jié)果解析模型構(gòu)建完成后調(diào)用求解器進(jìn)行求解。YALMIP使得這一步非常簡單。% 配置求解器選項(xiàng) ops sdpsettings(verbose, 1, solver, cplex, cplex.timelimit, 300); % 設(shè)置求解器為CPLEX時(shí)間限制300秒 ops.cplex.mip.tolerances.mipgap 1e-4; % 設(shè)置MIP相對間隙容差 % 求解優(yōu)化問題 diag optimize(Constraints, Objective, ops); % 檢查求解狀態(tài) if diag.problem 0 disp(優(yōu)化求解成功); % 獲取優(yōu)化結(jié)果 z_opt value(z); P_opt value(P); U_opt value(U); I_opt value(I); obj_opt value(Objective); % 分析結(jié)果找出需要操作的開關(guān) % 比較優(yōu)化開關(guān)狀態(tài)z_opt與初始開關(guān)狀態(tài)branch(:,6) initial_z branch(:, 6); % 假設(shè)第6列是初始狀態(tài) switch_operations find(z_opt ~ initial_z); disp([需要操作的開關(guān)支路編號(hào), num2str(switch_operations)]); % 計(jì)算網(wǎng)損降低百分比 % 需要先計(jì)算初始狀態(tài)下的潮流和網(wǎng)損可以調(diào)用一個(gè)前推回代潮流計(jì)算函數(shù) % [Ploss_initial, ~, ~] runPowerFlow(bus, branch, initial_z); % reduction (Ploss_initial - obj_opt) / Ploss_initial * 100; % disp([網(wǎng)損降低了 , num2str(reduction), %]); % 可視化重構(gòu)前后拓?fù)湫枰远x繪圖函數(shù) % plotNetworkTopology(bus, branch, initial_z, 初始拓?fù)?; % plotNetworkTopology(bus, branch, z_opt, 優(yōu)化后拓?fù)?; else disp(求解失敗或未找到最優(yōu)解。); disp(yalmiperror(diag.problem)); end注意大M值的選取這是混合整數(shù)建模中的關(guān)鍵技巧。M值不能太小否則可能錯(cuò)誤地割掉可行解也不能太大否則會(huì)惡化模型的線性松弛導(dǎo)致求解速度變慢甚至數(shù)值不穩(wěn)定。一個(gè)實(shí)用的方法是根據(jù)物理意義估算變量的上下界。例如P和Q的絕對值上界可以取系統(tǒng)總負(fù)荷的2倍I的上界可以取支路最大允許電流的平方電壓差的上界可以取(Vmax^2 - Vmin^2)。使用盡可能緊的界。4. 關(guān)鍵技巧、常見陷阱與性能優(yōu)化在實(shí)際編碼和求解過程中你會(huì)遇到許多理論模型不會(huì)提及的細(xì)節(jié)問題。以下是我從多次項(xiàng)目實(shí)踐中總結(jié)出的核心經(jīng)驗(yàn)。4.1 確保模型可行性與松弛緊致性一個(gè)常見的困擾是模型構(gòu)建好了但求解器報(bào)告“不可行”。這通常不是求解器的問題而是模型本身或數(shù)據(jù)有問題。檢查輻射狀約束虛擬流法是最容易實(shí)現(xiàn)的但要確保你正確構(gòu)建了節(jié)點(diǎn)-支路關(guān)聯(lián)關(guān)系并且根節(jié)點(diǎn)平衡節(jié)點(diǎn)被正確排除在虛擬流需求約束之外。一個(gè)快速驗(yàn)證方法是先固定所有開關(guān)狀態(tài)為初始閉合狀態(tài)z1去掉整數(shù)約束只求解連續(xù)潮流問題看是否可行。如果不可行問題很可能出在潮流約束部分。檢查二階錐松弛的緊致性求解完成后必須驗(yàn)證對于每條閉合的支路二階錐約束是否在最優(yōu)解處是緊的即取等號(hào)?？梢杂?jì)算松弛間隙for k 1:nl if z_opt(k) 0.5 % 認(rèn)為是閉合的支路 lhs sqrt((2*P_opt(k))^2 (2*Q_opt(k))^2 (U_opt(from)-I_opt(k))^2); rhs U_opt(from) I_opt(k); gap lhs - rhs; if abs(gap) 1e-4 % 設(shè)置一個(gè)小的容差 warning([支路 , num2str(k), 的二階錐松弛間隙較大: , num2str(gap)]); end end end如果發(fā)現(xiàn)大量支路的間隙不為零特別是數(shù)值較大時(shí)說明松弛不緊優(yōu)化結(jié)果可能不滿足原始物理方程是無效的。這通常發(fā)生在負(fù)荷極輕或網(wǎng)絡(luò)參數(shù)不合理的場景。此時(shí)可能需要考慮更精確的凸松弛方法如半正定規(guī)劃SDP松弛或回到非線性模型。檢查電壓和電流限值確保你設(shè)置的Vmax_pu和Vmin_pu以及Imax_pu是合理的。過于嚴(yán)格的限值可能導(dǎo)致問題不可行?？梢韵确艑捪拗登蠼饪纯醋顑?yōu)解是否自然滿足安全約束如果不滿足再分析是網(wǎng)絡(luò)本身能力不足還是模型問題。4.2 提升求解速度的策略MISOCP雖然是凸問題但對于大規(guī)模配電網(wǎng)如上千個(gè)節(jié)點(diǎn)求解時(shí)間仍然可能很長。以下策略可以顯著加速提供初始解如果你知道一個(gè)較好的初始拓?fù)浔热绠?dāng)前運(yùn)行拓?fù)淇梢詫⑵渥鳛镸IP求解器的初始解。在YALMIP中可以通過assign函數(shù)為變量賦初值。% 假設(shè)initial_z是初始開關(guān)狀態(tài)向量 assign(z, initial_z); % 然后基于這個(gè)初始狀態(tài)快速求解一個(gè)連續(xù)潮流固定z得到P,Q,U,I的初始值 % ... 求解固定z的SOCP問題 ... assign(P, P_init_val); assign(Q, Q_init_val); assign(U, U_init_val); assign(I, I_init_val); % 然后再調(diào)用包含整數(shù)變量的完整優(yōu)化 optimize(Constraints, Objective, ops);一個(gè)好的初始解可以極大縮短分支定界法的搜索過程。調(diào)整求解器參數(shù)mip.tolerances.mipgap相對間隙容差。默認(rèn)值如1e-4對于工程應(yīng)用通常足夠。如果只想快速得到一個(gè)可行解可以將其設(shè)為0.01或0.05。mip.strategy.nodeselect節(jié)點(diǎn)選擇策略。CPLEX提供了多種策略如深度優(yōu)先、最佳估計(jì)等對于不同問題效果不同可以嘗試調(diào)整。emphasis.mip設(shè)置求解重點(diǎn)。1平衡、2強(qiáng)調(diào)可行性、3強(qiáng)調(diào)最優(yōu)性、4隱藏輸出。如果遇到可行解難找的情況可以設(shè)置為2。模型簡化減少整數(shù)變量并非所有支路都需要建模為可操作的開關(guān)。通常只有聯(lián)絡(luò)開關(guān)和少數(shù)關(guān)鍵分段開關(guān)參與重構(gòu)。將常閉分段開關(guān)的z固定為1可以大幅減少整數(shù)變量數(shù)量。使用網(wǎng)絡(luò)流模型替代虛擬流對于輻射狀約束有更緊湊的建模方法如多商品流雖然約束數(shù)量可能略多但線性規(guī)劃松弛更強(qiáng)有時(shí)能加快求解。但這需要更復(fù)雜的建模。分步求解/啟發(fā)式與精確解結(jié)合對于大規(guī)模系統(tǒng)可以先使用快速的啟發(fā)式算法如遺傳算法、粒子群算法得到一個(gè)質(zhì)量不錯(cuò)的解然后以此解作為初始解再用MISOCP模型進(jìn)行局部精細(xì)優(yōu)化?；蛘邔⒋笠?guī)模網(wǎng)絡(luò)分解為幾個(gè)較小的區(qū)域分別優(yōu)化。4.3 MATLAB代碼調(diào)試與驗(yàn)證心得從小系統(tǒng)開始永遠(yuǎn)不要一開始就在IEEE 118節(jié)點(diǎn)系統(tǒng)上測試你的完整模型。先用一個(gè)3節(jié)點(diǎn)或5節(jié)點(diǎn)的微型系統(tǒng)驗(yàn)證你的模型邏輯。手動(dòng)計(jì)算這個(gè)微型系統(tǒng)在已知拓?fù)湎碌某绷骱途W(wǎng)損然后看你的優(yōu)化模型能否正確求解并得到一致的結(jié)果。善用value和check函數(shù)在求解后用value函數(shù)獲取變量值。對于約束可以用check(Constraints)來檢查每個(gè)約束的滿足情況它會(huì)返回每個(gè)約束左右兩邊的值這對于定位“不可行”問題的根源至關(guān)重要。可視化中間結(jié)果編寫簡單的函數(shù)來繪制網(wǎng)絡(luò)拓?fù)?、功率流方向和電壓分布。圖形化的結(jié)果比一堆數(shù)字更容易發(fā)現(xiàn)問題。例如一個(gè)環(huán)狀的結(jié)果圖立刻就能告訴你輻射狀約束可能沒起作用。與專業(yè)工具對比如果條件允許將你的優(yōu)化結(jié)果開關(guān)狀態(tài)輸入到專業(yè)的電力系統(tǒng)仿真軟件如OpenDSS、MATPOWER中進(jìn)行潮流計(jì)算對比網(wǎng)損和電壓分布。這是驗(yàn)證模型有效性的黃金標(biāo)準(zhǔn)。5. 項(xiàng)目擴(kuò)展與實(shí)際應(yīng)用思考基于MISOCP的配電網(wǎng)重構(gòu)模型是一個(gè)強(qiáng)大的基礎(chǔ)框架它可以很容易地?cái)U(kuò)展到更復(fù)雜的實(shí)際應(yīng)用場景中?？紤]分布式電源如果配電網(wǎng)中包含光伏、風(fēng)機(jī)等分布式電源只需在對應(yīng)節(jié)點(diǎn)的功率平衡方程中將p_j和q_j修改為(負(fù)荷 - 電源出力)。如果電源出力是可調(diào)的如可控逆變器則可以將其作為優(yōu)化變量與重構(gòu)協(xié)同優(yōu)化實(shí)現(xiàn)“網(wǎng)源協(xié)同”?？紤]時(shí)變性動(dòng)態(tài)重構(gòu)負(fù)荷和分布式電源出力是隨時(shí)間變化的?？梢詫⒁惶靹澐譃槎鄠€(gè)時(shí)段如24小時(shí)為每個(gè)時(shí)段建立一個(gè)MISOCP模型并添加時(shí)段間的耦合約束如開關(guān)操作次數(shù)限制一個(gè)開關(guān)在相鄰時(shí)段內(nèi)狀態(tài)變化不能超過一次。這就形成了一個(gè)大規(guī)模的多時(shí)段MISOCP問題雖然計(jì)算復(fù)雜但能實(shí)現(xiàn)全局最優(yōu)的日內(nèi)調(diào)度。與無功優(yōu)化結(jié)合可以在模型中引入電容器組投切、變壓器分接頭調(diào)整的無功優(yōu)化變量與開關(guān)重構(gòu)同時(shí)進(jìn)行優(yōu)化同時(shí)降低網(wǎng)損和改善電壓形成綜合優(yōu)化。故障后的恢復(fù)重構(gòu)當(dāng)網(wǎng)絡(luò)發(fā)生故障后故障區(qū)段會(huì)被隔離。重構(gòu)的目標(biāo)轉(zhuǎn)變?yōu)樵谑Ｓ嗟慕】稻W(wǎng)絡(luò)中通過操作開關(guān)盡可能多地恢復(fù)非故障失電區(qū)域的供電。此時(shí)目標(biāo)函數(shù)可以是最小化失電負(fù)荷量或者加權(quán)失電負(fù)荷量。約束中需要加入故障隔離的約束故障支路z固定為0。最后我想分享一個(gè)深刻的體會(huì)數(shù)學(xué)上的優(yōu)雅和計(jì)算上的高效有時(shí)需要向工程實(shí)用性妥協(xié)。MISOCP模型給出的最優(yōu)拓?fù)湓诶碚撋暇W(wǎng)損最小但在實(shí)際系統(tǒng)中開關(guān)操作本身有成本頻繁操作會(huì)降低設(shè)備壽命。因此一個(gè)更實(shí)用的模型應(yīng)該在目標(biāo)函數(shù)中加入開關(guān)操作次數(shù)或狀態(tài)的懲罰項(xiàng)。此外求解時(shí)間必須滿足在線或準(zhǔn)在線應(yīng)用的要求。這意味著我們可能需要在模型精確度、求解速度和解決方案的魯棒性之間做出權(quán)衡。這套MISOCP方法論的價(jià)值不僅在于它能給出一個(gè)“最優(yōu)”的參考解更在于它為我們提供了一個(gè)嚴(yán)格、靈活的分析框架讓我們能夠系統(tǒng)地、量化地評(píng)估各種網(wǎng)絡(luò)配置方案的優(yōu)劣這是傳統(tǒng)經(jīng)驗(yàn)主義方法無法比擬的。 基于混合整數(shù)二階錐規(guī)劃的配電網(wǎng)重構(gòu)：模型、實(shí)現(xiàn)與優(yōu)化
1. 項(xiàng)目概述當(dāng)配電網(wǎng)重構(gòu)遇上混合整數(shù)二階錐在電力系統(tǒng)運(yùn)行領(lǐng)域配電網(wǎng)重構(gòu)是一個(gè)經(jīng)典且核心的優(yōu)化問題。簡單來說它就像是在一個(gè)復(fù)雜的、由眾多開關(guān)分段開關(guān)和聯(lián)絡(luò)開關(guān)組成的城市道路網(wǎng)絡(luò)中通過改變某些路段的通斷狀態(tài)來重新規(guī)劃電力流動(dòng)的路徑。這么做的目標(biāo)很明確降低整個(gè)網(wǎng)絡(luò)的電能損耗、改善電壓質(zhì)量、平衡各條線路的負(fù)載或者在故障后快速恢復(fù)供電。傳統(tǒng)上這個(gè)問題被建模為一個(gè)混合整數(shù)非線性規(guī)劃問題因?yàn)殚_關(guān)狀態(tài)是0或1的整數(shù)變量而網(wǎng)絡(luò)中的功率流方程則是非線性的。非線性尤其是非凸性是求解這類問題的“攔路虎”它意味著你可能找到一個(gè)局部最優(yōu)解但無法保證它是全局最優(yōu)的而且求解過程往往非常耗時(shí)?；旌险麛?shù)二階錐規(guī)劃的出現(xiàn)為這個(gè)難題提供了一個(gè)極具吸引力的突破口。它的核心思想是在一定的假設(shè)和近似下比如采用DistFlow潮流模型可以將原本非線性的功率平衡方程轉(zhuǎn)化為一系列二階錐約束。二階錐約束是一種特殊的凸約束其形式類似于一個(gè)“冰淇淋蛋筒”的形狀。將問題轉(zhuǎn)化為MISOCP模型后我們實(shí)際上是在一個(gè)凸的可行域內(nèi)尋找最優(yōu)的整數(shù)解。凸性保證了局部最優(yōu)即全局最優(yōu)并且現(xiàn)代的商業(yè)求解器如CPLEX、Gurobi對這類問題有非常高效的求解算法。因此基于MISOCP的配電網(wǎng)重構(gòu)本質(zhì)上是用一種可高效求解的凸優(yōu)化框架去逼近和求解原本棘手的非凸問題在計(jì)算效率和求解質(zhì)量之間取得了出色的平衡。對于電力系統(tǒng)研究人員和工程師而言掌握這套方法意味著你手里多了一把解決復(fù)雜網(wǎng)絡(luò)優(yōu)化問題的“瑞士軍刀”。2. 核心模型構(gòu)建從物理方程到數(shù)學(xué)規(guī)劃要理解MISOCP模型我們必須從配電網(wǎng)的物理基礎(chǔ)——潮流模型開始。在配電網(wǎng)重構(gòu)中我們通常采用基于支路潮流的DistFlow模型它比傳統(tǒng)的節(jié)點(diǎn)功率方程更直觀地描述了功率在輻射狀網(wǎng)絡(luò)中的流動(dòng)。2.1 DistFlow潮流方程及其凸松弛考慮一個(gè)由節(jié)點(diǎn)母線和支路線路組成的配電網(wǎng)。對于任意支路(i, j)設(shè)其首端節(jié)點(diǎn)i的電壓幅值為V_i末端節(jié)點(diǎn)j的電壓幅值為V_j支路上流過的有功功率為P_ij無功功率為Q_ij支路電阻為r_ij電抗為x_ij。經(jīng)典的DistFlow方程可以寫為P_ij - r_ij * I_ij^2 Σ_{k∈δ(j)} P_jk p_j Q_ij - x_ij * I_ij^2 Σ_{k∈δ(j)} Q_jk q_j V_j^2 V_i^2 - 2*(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*I_ij^2其中p_j和q_j是節(jié)點(diǎn)j的凈注入有功和無功功率負(fù)荷為負(fù)分布式電源為正δ(j)表示以j為首端的所有下游支路集合I_ij是支路電流幅值。第三個(gè)方程中包含了電流平方項(xiàng)I_ij^2它與功率的關(guān)系是非線性的I_ij^2 (P_ij^2 Q_ij^2) / V_i^2。正是這個(gè)項(xiàng)引入了非凸性。為了進(jìn)行凸松弛我們引入一組新的變量進(jìn)行替換令 u_i V_i^2代表節(jié)點(diǎn)電壓的平方。令 l_ij I_ij^2代表支路電流的平方。支路功率P_ij和Q_ij保持不變。經(jīng)過變量替換和忽略高階項(xiàng)(r_ij^2 x_ij^2)*l_ij因其值通常很小DistFlow方程可以重寫為P_ij - r_ij * l_ij Σ_{k∈δ(j)} P_jk p_j Q_ij - x_ij * l_ij Σ_{k∈δ(j)} Q_jk q_j u_j u_i - 2*(r_ij*P_ij x_ij*Q_ij) (r_ij^2 x_ij^2)*l_ij此時(shí)非凸性隱藏在變量之間的關(guān)系中根據(jù)歐姆定律和功率定義應(yīng)有 P_ij^2 Q_ij^2 u_i * l_ij。這是一個(gè)旋轉(zhuǎn)錐約束是非凸的。二階錐松弛的關(guān)鍵一步來了我們將這個(gè)等式約束松弛為一個(gè)不等式約束P_ij^2 Q_ij^2 ≤ u_i * l_ij這個(gè)不等式可以等價(jià)地寫成一個(gè)標(biāo)準(zhǔn)二階錐約束的形式|| [2P_ij, 2Q_ij, u_i - l_ij]^T ||_2 ≤ u_i l_ij其中||·||_2表示向量的二范數(shù)。這個(gè)約束定義了一個(gè)凸集二階錐。在輻射狀配電網(wǎng)且負(fù)荷水平不是極端輕載的情況下這個(gè)松弛通常是緊的即最優(yōu)解會(huì)自動(dòng)使不等式取等號(hào)從而滿足原始物理方程。這就成功地將非凸問題嵌入了凸框架內(nèi)。2.2 網(wǎng)絡(luò)拓?fù)渑c開關(guān)建模配電網(wǎng)重構(gòu)的核心是改變網(wǎng)絡(luò)拓?fù)溥@通過改變支路的連接狀態(tài)即開關(guān)的開合來實(shí)現(xiàn)。我們需要引入0-1整數(shù)變量來建模這種狀態(tài)。對于網(wǎng)絡(luò)中的每一條可能的支路包括現(xiàn)有的常閉支路和作為備用的常開支路我們定義一個(gè)二進(jìn)制變量z_ij ∈ {0, 1}。z_ij 1表示支路(i, j)閉合通電z_ij 0表示支路(i, j)斷開。網(wǎng)絡(luò)拓?fù)湫枰獫M足兩個(gè)核心約束輻射狀約束重構(gòu)后的網(wǎng)絡(luò)必須是一個(gè)輻射狀樹形結(jié)構(gòu)即無環(huán)且連通。這可以通過多種方式建模最常見的是利用配電網(wǎng)絡(luò)單電源點(diǎn)的特點(diǎn)采用“虛擬流”或“父節(jié)點(diǎn)指示”方法。虛擬流法為每個(gè)節(jié)點(diǎn)除根節(jié)點(diǎn)外分配一個(gè)單位的虛擬需求根節(jié)點(diǎn)變電站為虛擬源。要求虛擬流只能從閉合的支路上通過并且每個(gè)非根節(jié)點(diǎn)必須有且僅有一條流入的閉合支路。這保證了連通性和無環(huán)性。父節(jié)點(diǎn)變量法為每個(gè)節(jié)點(diǎn)定義一個(gè)整數(shù)變量表示其父節(jié)點(diǎn)編號(hào)并約束根節(jié)點(diǎn)無父節(jié)點(diǎn)其他節(jié)點(diǎn)有且僅有一個(gè)父節(jié)點(diǎn)且父節(jié)點(diǎn)必須通過閉合的支路與之相連。這種方法更直觀但變量更多。開關(guān)邏輯約束支路的電氣狀態(tài)必須與開關(guān)狀態(tài)一致。當(dāng)z_ij0開關(guān)斷開時(shí)該支路上的功率P_ij, Q_ij必須為0電流平方l_ij也必須為0。這可以用大M法來建模-M * z_ij ≤ P_ij ≤ M * z_ij -M * z_ij ≤ Q_ij ≤ M * z_ij 0 ≤ l_ij ≤ M * z_ij 0 ≤ u_i - u_j ≤ M * z_ij ΔV_max*(1-z_ij) // 可選斷開時(shí)電壓差約束可放松其中M是一個(gè)足夠大的正數(shù)但為了數(shù)值穩(wěn)定性應(yīng)盡可能選取緊的界。2.3 目標(biāo)函數(shù)與完整MISOCP模型配電網(wǎng)重構(gòu)的典型目標(biāo)是最小化系統(tǒng)總有功網(wǎng)損。網(wǎng)損等于所有支路電阻上的損耗之和Ploss Σ_{(i,j)} r_ij * l_ij。由于l_ij是電流的平方這個(gè)目標(biāo)函數(shù)是線性的。綜合以上所有部分我們可以得到完整的基于MISOCP的配電網(wǎng)重構(gòu)模型目標(biāo)函數(shù)Minimize Σ_{(i,j)∈Ω} r_ij * l_ij約束條件節(jié)點(diǎn)功率平衡線性等式 Σ_{i∈π(j)} P_ij - Σ_{k∈δ(j)} P_jk - r_ij * l_ij p_j, ?j Σ_{i∈π(j)} Q_ij - Σ_{k∈δ(j)} Q_jk - x_ij * l_ij q_j, ?j π(j)表示節(jié)點(diǎn)j的所有父節(jié)點(diǎn)支路集合電壓降落方程線性等式 u_j u_i - 2*(r_ijP_ij x_ijQ_ij) (r_ij^2 x_ij^2)*l_ij, ?(i,j)∈Ω二階錐松弛約束凸錐約束 || [2P_ij, 2Q_ij, u_i - l_ij]^T ||_2 ≤ u_i l_ij, ?(i,j)∈Ω開關(guān)邏輯約束大M法混合整數(shù)線性 -M * z_ij ≤ P_ij ≤ M * z_ij, ?(i,j)∈Ω -M * z_ij ≤ Q_ij ≤ M * z_ij, ?(i,j)∈Ω 0 ≤ l_ij ≤ M * z_ij, ?(i,j)∈Ω輻射狀拓?fù)浼s束混合整數(shù)線性以虛擬流法為例 Σ_{i∈π(j)} f_ij 1, ?j ≠ 根節(jié)點(diǎn) 0 ≤ f_ij ≤ N * z_ij, ?(i,j)∈Ω // f_ij為從i到j(luò)的虛擬流N為節(jié)點(diǎn)數(shù) Σ_{(i,j)∈Ω} z_ij N - 1 // 輻射狀網(wǎng)絡(luò)支路數(shù)等于節(jié)點(diǎn)數(shù)減一運(yùn)行安全約束線性不等式 u_min^2 ≤ u_i ≤ u_max^2, ?i // 電壓幅值平方約束 l_ij ≤ (I_ij_max)^2 * z_ij, ?(i,j)∈Ω // 支路電流容量約束整數(shù)變量約束 z_ij ∈ {0, 1}, ?(i,j)∈Ω這個(gè)模型就是一個(gè)標(biāo)準(zhǔn)的混合整數(shù)二階錐規(guī)劃問題。目標(biāo)函數(shù)和絕大多數(shù)約束是線性的關(guān)鍵的非線性部分被包含在凸的二階錐約束中因此可以被高效的商業(yè)求解器處理。3. 基于MATLAB與YALMIP/CPLEX的代碼實(shí)現(xiàn)詳解理論模型建立后我們需要一個(gè)強(qiáng)大的工具鏈將其轉(zhuǎn)化為可執(zhí)行的代碼。MATLAB因其強(qiáng)大的數(shù)學(xué)計(jì)算和矩陣操作能力成為算法原型開發(fā)的理想環(huán)境。YALMIP是一個(gè)在MATLAB中用于建模優(yōu)化問題的免費(fèi)工具箱它提供了一種非常直觀的方式來描述優(yōu)化問題變量、目標(biāo)、約束然后自動(dòng)調(diào)用后端求解器如CPLEX、Gurobi、MOSEK來求解。CPLEX是IBM推出的高性能數(shù)學(xué)規(guī)劃求解器對MISOCP有非常好的支持。3.1 開發(fā)環(huán)境搭建與數(shù)據(jù)準(zhǔn)備首先確保你的MATLAB環(huán)境已經(jīng)安裝了YALMIP和CPLEX。YALMIP可以直接從其官網(wǎng)下載并添加到MATLAB路徑。CPLEX需要從IBM官網(wǎng)獲取學(xué)術(shù)版或商業(yè)版許可并安裝。安裝后在MATLAB中運(yùn)行yalmiptestYALMIP會(huì)自動(dòng)檢測可用的求解器確認(rèn)CPLEX被正確識(shí)別。數(shù)據(jù)準(zhǔn)備是第一步也是最容易出錯(cuò)的一步。我們需要一個(gè)標(biāo)準(zhǔn)的測試配電網(wǎng)數(shù)據(jù)例如IEEE 33節(jié)點(diǎn)、69節(jié)點(diǎn)或118節(jié)點(diǎn)系統(tǒng)。數(shù)據(jù)通常包括bus矩陣節(jié)點(diǎn)編號(hào)、類型平衡節(jié)點(diǎn)/負(fù)荷節(jié)點(diǎn)、有功負(fù)荷、無功負(fù)荷。branch矩陣支路首末端節(jié)點(diǎn)編號(hào)、電阻、電抗、額定電流、初始狀態(tài)0開/1合。baseMVA和baseKV系統(tǒng)基準(zhǔn)值用于標(biāo)幺化。我強(qiáng)烈建議將原始數(shù)據(jù)如.mat文件或Excel表格讀入后先進(jìn)行標(biāo)幺化處理。將所有阻抗、功率、電壓、電流值除以相應(yīng)的基準(zhǔn)值。標(biāo)幺化能顯著改善模型的數(shù)值穩(wěn)定性避免因?qū)嶋H數(shù)據(jù)量綱差異過大如電阻是0.幾歐姆功率是幾兆瓦導(dǎo)致求解器出現(xiàn)數(shù)值問題。% 示例數(shù)據(jù)讀取與標(biāo)幺化 load(IEEE33.mat); % 假設(shè)數(shù)據(jù)已保存在此文件 baseMVA 10; % 基準(zhǔn)功率MVA baseKV 12.66; % 基準(zhǔn)電壓kV Zbase baseKV^2 / baseMVA; % 基準(zhǔn)阻抗 % 標(biāo)幺化支路電阻和電抗 branch(:, 3) branch(:, 3) / Zbase; % 電阻R branch(:, 4) branch(:, 4) / Zbase; % 電抗X % 標(biāo)幺化節(jié)點(diǎn)負(fù)荷 bus(:, 3) bus(:, 3) / baseMVA; % 有功負(fù)荷P bus(:, 4) bus(:, 4) / baseMVA; % 無功負(fù)荷Q % 標(biāo)幺化電壓和電流限值 Vmax_pu 1.05; Vmin_pu 0.95; Imax_pu branch(:, 5) / (baseMVA / (sqrt(3)*baseKV)); % 假設(shè)branch第5列是額定電流(A)3.2 使用YALMIP構(gòu)建MISOCP模型YALMIP建模的核心是三步定義變量、定義約束、定義目標(biāo)函數(shù)。代碼的清晰度至關(guān)重要。% 1. 定義變量 nb size(bus, 1); % 節(jié)點(diǎn)數(shù) nl size(branch, 1); % 支路數(shù)包括所有可能的支路 % 連續(xù)變量 P sdpvar(nl, 1); % 支路有功潮流 Q sdpvar(nl, 1); % 支路無功潮流 U sdpvar(nb, 1); % 節(jié)點(diǎn)電壓平方 I sdpvar(nl, 1); % 支路電流平方 % 整數(shù)變量二進(jìn)制 z binvar(nl, 1); % 支路開關(guān)狀態(tài)1閉合0斷開 % 虛擬流變量用于輻射狀約束 f sdpvar(nl, 1); % 虛擬流非負(fù)連續(xù)變量 % 2. 定義約束 Constraints []; % 2.1 節(jié)點(diǎn)功率平衡約束 % 構(gòu)建節(jié)點(diǎn)-支路關(guān)聯(lián)矩陣Anb x nlA(i,k)1表示支路k以i為首端-1表示以i為末端 A zeros(nb, nl); for k 1:nl from branch(k, 1); to branch(k, 2); A(from, k) 1; A(to, k) -1; end % 平衡節(jié)點(diǎn)假設(shè)為節(jié)點(diǎn)1電壓固定 Constraints [Constraints, U(1) 1.0]; % 標(biāo)幺值通常設(shè)為1.0 for i 1:nb % 流入節(jié)點(diǎn)i的凈功率 該節(jié)點(diǎn)負(fù)荷 if i 1 % 平衡節(jié)點(diǎn)是功率源其注入功率是自由的由優(yōu)化決定 % 實(shí)際上平衡節(jié)點(diǎn)的功率平衡方程通常不顯式添加因?yàn)樗鼤?huì)自動(dòng)滿足 % 我們只需處理負(fù)荷節(jié)點(diǎn) else % 對于負(fù)荷節(jié)點(diǎn)從關(guān)聯(lián)矩陣中找出與該節(jié)點(diǎn)相關(guān)的支路 idx_in find(A(i, :) 1); % 以i為首端的支路功率流出i idx_out find(A(i, :) -1); % 以i為末端的支路功率流入i % 功率平衡流入 - 流出 - 損耗 負(fù)荷 % 注意DistFlow模型中損耗(r*I)發(fā)生在支路末端 P_in sum(P(idx_out)); P_out sum(P(idx_in)); P_loss_sum sum(branch(idx_out, 3) .* I(idx_out)); % 所有流入支路的損耗之和 Constraints [Constraints, P_in - P_out - P_loss_sum bus(i, 3)]; % bus(i,3)是有功負(fù)荷 Q_in sum(Q(idx_out)); Q_out sum(Q(idx_in)); Q_loss_sum sum(branch(idx_out, 4) .* I(idx_out)); % 所有流入支路的無功損耗之和 Constraints [Constraints, Q_in - Q_out - Q_loss_sum bus(i, 4)]; % bus(i,4)是無功負(fù)荷 end end % 2.2 電壓降落方程約束 for k 1:nl from branch(k, 1); to branch(k, 2); r branch(k, 3); x branch(k, 4); Constraints [Constraints, U(to) U(from) - 2*(r*P(k) x*Q(k)) (r^2 x^2)*I(k)]; end % 2.3 二階錐松弛約束 for k 1:nl from branch(k, 1); Constraints [Constraints, cone([2*P(k); 2*Q(k); U(from)-I(k)], U(from)I(k))]; % YALMIP中的cone函數(shù)直接用于構(gòu)建二階錐約束||x|| y % 這里 x [2P; 2Q; U-I], y UI end % 2.4 開關(guān)邏輯約束大M法 M_pq 10; % 根據(jù)系統(tǒng)基準(zhǔn)功率估算的足夠大的數(shù)例如最大可能潮流的2倍 M_i 2; % 電流平方的上界例如(2*Imax)^2 M_u 0.5; % 電壓差的上界 for k 1:nl Constraints [Constraints, -M_pq * z(k) P(k) M_pq * z(k)]; Constraints [Constraints, -M_pq * z(k) Q(k) M_pq * z(k)]; Constraints [Constraints, 0 I(k) M_i * z(k)]; % 對于斷開支路可以放松電壓差約束 from branch(k, 1); to branch(k, 2); Constraints [Constraints, -M_u*(1-z(k)) U(from)-U(to) M_u*(1-z(k))]; end % 2.5 輻射狀拓?fù)浼s束虛擬流法 % 假設(shè)節(jié)點(diǎn)1是根節(jié)點(diǎn)平衡節(jié)點(diǎn) for j 2:nb % 找到所有以j為末端的支路索引 idx_to_j find(branch(:,2) j); Constraints [Constraints, sum(f(idx_to_j)) 1]; end for k 1:nl Constraints [Constraints, 0 f(k) nb * z(k)]; end Constraints [Constraints, sum(z) nb-1]; % 輻射狀網(wǎng)絡(luò)支路數(shù)約束 % 2.6 運(yùn)行安全約束 for i 1:nb Constraints [Constraints, Vmin_pu^2 U(i) Vmax_pu^2]; end for k 1:nl Constraints [Constraints, I(k) (Imax_pu(k))^2 * z(k)]; end % 3. 定義目標(biāo)函數(shù) Objective sum(branch(:,3) .* I); % 最小化總有功網(wǎng)損3.3 模型求解與結(jié)果解析模型構(gòu)建完成后調(diào)用求解器進(jìn)行求解。YALMIP使得這一步非常簡單。% 配置求解器選項(xiàng) ops sdpsettings(verbose, 1, solver, cplex, cplex.timelimit, 300); % 設(shè)置求解器為CPLEX時(shí)間限制300秒 ops.cplex.mip.tolerances.mipgap 1e-4; % 設(shè)置MIP相對間隙容差 % 求解優(yōu)化問題 diag optimize(Constraints, Objective, ops); % 檢查求解狀態(tài) if diag.problem 0 disp(優(yōu)化求解成功); % 獲取優(yōu)化結(jié)果 z_opt value(z); P_opt value(P); U_opt value(U); I_opt value(I); obj_opt value(Objective); % 分析結(jié)果找出需要操作的開關(guān) % 比較優(yōu)化開關(guān)狀態(tài)z_opt與初始開關(guān)狀態(tài)branch(:,6) initial_z branch(:, 6); % 假設(shè)第6列是初始狀態(tài) switch_operations find(z_opt ~ initial_z); disp([需要操作的開關(guān)支路編號(hào), num2str(switch_operations)]); % 計(jì)算網(wǎng)損降低百分比 % 需要先計(jì)算初始狀態(tài)下的潮流和網(wǎng)損可以調(diào)用一個(gè)前推回代潮流計(jì)算函數(shù) % [Ploss_initial, ~, ~] runPowerFlow(bus, branch, initial_z); % reduction (Ploss_initial - obj_opt) / Ploss_initial * 100; % disp([網(wǎng)損降低了 , num2str(reduction), %]); % 可視化重構(gòu)前后拓?fù)湫枰远x繪圖函數(shù) % plotNetworkTopology(bus, branch, initial_z, 初始拓?fù)?; % plotNetworkTopology(bus, branch, z_opt, 優(yōu)化后拓?fù)?; else disp(求解失敗或未找到最優(yōu)解。); disp(yalmiperror(diag.problem)); end注意大M值的選取這是混合整數(shù)建模中的關(guān)鍵技巧。M值不能太小否則可能錯(cuò)誤地割掉可行解也不能太大否則會(huì)惡化模型的線性松弛導(dǎo)致求解速度變慢甚至數(shù)值不穩(wěn)定。一個(gè)實(shí)用的方法是根據(jù)物理意義估算變量的上下界。例如P和Q的絕對值上界可以取系統(tǒng)總負(fù)荷的2倍I的上界可以取支路最大允許電流的平方電壓差的上界可以取(Vmax^2 - Vmin^2)。使用盡可能緊的界。4. 關(guān)鍵技巧、常見陷阱與性能優(yōu)化在實(shí)際編碼和求解過程中你會(huì)遇到許多理論模型不會(huì)提及的細(xì)節(jié)問題。以下是我從多次項(xiàng)目實(shí)踐中總結(jié)出的核心經(jīng)驗(yàn)。4.1 確保模型可行性與松弛緊致性一個(gè)常見的困擾是模型構(gòu)建好了但求解器報(bào)告“不可行”。這通常不是求解器的問題而是模型本身或數(shù)據(jù)有問題。檢查輻射狀約束虛擬流法是最容易實(shí)現(xiàn)的但要確保你正確構(gòu)建了節(jié)點(diǎn)-支路關(guān)聯(lián)關(guān)系并且根節(jié)點(diǎn)平衡節(jié)點(diǎn)被正確排除在虛擬流需求約束之外。一個(gè)快速驗(yàn)證方法是先固定所有開關(guān)狀態(tài)為初始閉合狀態(tài)z1去掉整數(shù)約束只求解連續(xù)潮流問題看是否可行。如果不可行問題很可能出在潮流約束部分。檢查二階錐松弛的緊致性求解完成后必須驗(yàn)證對于每條閉合的支路二階錐約束是否在最優(yōu)解處是緊的即取等號(hào)?？梢杂?jì)算松弛間隙for k 1:nl if z_opt(k) 0.5 % 認(rèn)為是閉合的支路 lhs sqrt((2*P_opt(k))^2 (2*Q_opt(k))^2 (U_opt(from)-I_opt(k))^2); rhs U_opt(from) I_opt(k); gap lhs - rhs; if abs(gap) 1e-4 % 設(shè)置一個(gè)小的容差 warning([支路 , num2str(k), 的二階錐松弛間隙較大: , num2str(gap)]); end end end如果發(fā)現(xiàn)大量支路的間隙不為零特別是數(shù)值較大時(shí)說明松弛不緊優(yōu)化結(jié)果可能不滿足原始物理方程是無效的。這通常發(fā)生在負(fù)荷極輕或網(wǎng)絡(luò)參數(shù)不合理的場景。此時(shí)可能需要考慮更精確的凸松弛方法如半正定規(guī)劃SDP松弛或回到非線性模型。檢查電壓和電流限值確保你設(shè)置的Vmax_pu和Vmin_pu以及Imax_pu是合理的。過于嚴(yán)格的限值可能導(dǎo)致問題不可行?？梢韵确艑捪拗登蠼饪纯醋顑?yōu)解是否自然滿足安全約束如果不滿足再分析是網(wǎng)絡(luò)本身能力不足還是模型問題。4.2 提升求解速度的策略MISOCP雖然是凸問題但對于大規(guī)模配電網(wǎng)如上千個(gè)節(jié)點(diǎn)求解時(shí)間仍然可能很長。以下策略可以顯著加速提供初始解如果你知道一個(gè)較好的初始拓?fù)浔热绠?dāng)前運(yùn)行拓?fù)淇梢詫⑵渥鳛镸IP求解器的初始解。在YALMIP中可以通過assign函數(shù)為變量賦初值。% 假設(shè)initial_z是初始開關(guān)狀態(tài)向量 assign(z, initial_z); % 然后基于這個(gè)初始狀態(tài)快速求解一個(gè)連續(xù)潮流固定z得到P,Q,U,I的初始值 % ... 求解固定z的SOCP問題 ... assign(P, P_init_val); assign(Q, Q_init_val); assign(U, U_init_val); assign(I, I_init_val); % 然后再調(diào)用包含整數(shù)變量的完整優(yōu)化 optimize(Constraints, Objective, ops);一個(gè)好的初始解可以極大縮短分支定界法的搜索過程。調(diào)整求解器參數(shù)mip.tolerances.mipgap相對間隙容差。默認(rèn)值如1e-4對于工程應(yīng)用通常足夠。如果只想快速得到一個(gè)可行解可以將其設(shè)為0.01或0.05。mip.strategy.nodeselect節(jié)點(diǎn)選擇策略。CPLEX提供了多種策略如深度優(yōu)先、最佳估計(jì)等對于不同問題效果不同可以嘗試調(diào)整。emphasis.mip設(shè)置求解重點(diǎn)。1平衡、2強(qiáng)調(diào)可行性、3強(qiáng)調(diào)最優(yōu)性、4隱藏輸出。如果遇到可行解難找的情況可以設(shè)置為2。模型簡化減少整數(shù)變量并非所有支路都需要建模為可操作的開關(guān)。通常只有聯(lián)絡(luò)開關(guān)和少數(shù)關(guān)鍵分段開關(guān)參與重構(gòu)。將常閉分段開關(guān)的z固定為1可以大幅減少整數(shù)變量數(shù)量。使用網(wǎng)絡(luò)流模型替代虛擬流對于輻射狀約束有更緊湊的建模方法如多商品流雖然約束數(shù)量可能略多但線性規(guī)劃松弛更強(qiáng)有時(shí)能加快求解。但這需要更復(fù)雜的建模。分步求解/啟發(fā)式與精確解結(jié)合對于大規(guī)模系統(tǒng)可以先使用快速的啟發(fā)式算法如遺傳算法、粒子群算法得到一個(gè)質(zhì)量不錯(cuò)的解然后以此解作為初始解再用MISOCP模型進(jìn)行局部精細(xì)優(yōu)化?；蛘邔⒋笠?guī)模網(wǎng)絡(luò)分解為幾個(gè)較小的區(qū)域分別優(yōu)化。4.3 MATLAB代碼調(diào)試與驗(yàn)證心得從小系統(tǒng)開始永遠(yuǎn)不要一開始就在IEEE 118節(jié)點(diǎn)系統(tǒng)上測試你的完整模型。先用一個(gè)3節(jié)點(diǎn)或5節(jié)點(diǎn)的微型系統(tǒng)驗(yàn)證你的模型邏輯。手動(dòng)計(jì)算這個(gè)微型系統(tǒng)在已知拓?fù)湎碌某绷骱途W(wǎng)損然后看你的優(yōu)化模型能否正確求解并得到一致的結(jié)果。善用value和check函數(shù)在求解后用value函數(shù)獲取變量值。對于約束可以用check(Constraints)來檢查每個(gè)約束的滿足情況它會(huì)返回每個(gè)約束左右兩邊的值這對于定位“不可行”問題的根源至關(guān)重要。可視化中間結(jié)果編寫簡單的函數(shù)來繪制網(wǎng)絡(luò)拓?fù)?、功率流方向和電壓分布。圖形化的結(jié)果比一堆數(shù)字更容易發(fā)現(xiàn)問題。例如一個(gè)環(huán)狀的結(jié)果圖立刻就能告訴你輻射狀約束可能沒起作用。與專業(yè)工具對比如果條件允許將你的優(yōu)化結(jié)果開關(guān)狀態(tài)輸入到專業(yè)的電力系統(tǒng)仿真軟件如OpenDSS、MATPOWER中進(jìn)行潮流計(jì)算對比網(wǎng)損和電壓分布。這是驗(yàn)證模型有效性的黃金標(biāo)準(zhǔn)。5. 項(xiàng)目擴(kuò)展與實(shí)際應(yīng)用思考基于MISOCP的配電網(wǎng)重構(gòu)模型是一個(gè)強(qiáng)大的基礎(chǔ)框架它可以很容易地?cái)U(kuò)展到更復(fù)雜的實(shí)際應(yīng)用場景中?？紤]分布式電源如果配電網(wǎng)中包含光伏、風(fēng)機(jī)等分布式電源只需在對應(yīng)節(jié)點(diǎn)的功率平衡方程中將p_j和q_j修改為(負(fù)荷 - 電源出力)。如果電源出力是可調(diào)的如可控逆變器則可以將其作為優(yōu)化變量與重構(gòu)協(xié)同優(yōu)化實(shí)現(xiàn)“網(wǎng)源協(xié)同”?？紤]時(shí)變性動(dòng)態(tài)重構(gòu)負(fù)荷和分布式電源出力是隨時(shí)間變化的?？梢詫⒁惶靹澐譃槎鄠€(gè)時(shí)段如24小時(shí)為每個(gè)時(shí)段建立一個(gè)MISOCP模型并添加時(shí)段間的耦合約束如開關(guān)操作次數(shù)限制一個(gè)開關(guān)在相鄰時(shí)段內(nèi)狀態(tài)變化不能超過一次。這就形成了一個(gè)大規(guī)模的多時(shí)段MISOCP問題雖然計(jì)算復(fù)雜但能實(shí)現(xiàn)全局最優(yōu)的日內(nèi)調(diào)度。與無功優(yōu)化結(jié)合可以在模型中引入電容器組投切、變壓器分接頭調(diào)整的無功優(yōu)化變量與開關(guān)重構(gòu)同時(shí)進(jìn)行優(yōu)化同時(shí)降低網(wǎng)損和改善電壓形成綜合優(yōu)化。故障后的恢復(fù)重構(gòu)當(dāng)網(wǎng)絡(luò)發(fā)生故障后故障區(qū)段會(huì)被隔離。重構(gòu)的目標(biāo)轉(zhuǎn)變?yōu)樵谑Ｓ嗟慕】稻W(wǎng)絡(luò)中通過操作開關(guān)盡可能多地恢復(fù)非故障失電區(qū)域的供電。此時(shí)目標(biāo)函數(shù)可以是最小化失電負(fù)荷量或者加權(quán)失電負(fù)荷量。約束中需要加入故障隔離的約束故障支路z固定為0。最后我想分享一個(gè)深刻的體會(huì)數(shù)學(xué)上的優(yōu)雅和計(jì)算上的高效有時(shí)需要向工程實(shí)用性妥協(xié)。MISOCP模型給出的最優(yōu)拓?fù)湓诶碚撋暇W(wǎng)損最小但在實(shí)際系統(tǒng)中開關(guān)操作本身有成本頻繁操作會(huì)降低設(shè)備壽命。因此一個(gè)更實(shí)用的模型應(yīng)該在目標(biāo)函數(shù)中加入開關(guān)操作次數(shù)或狀態(tài)的懲罰項(xiàng)。此外求解時(shí)間必須滿足在線或準(zhǔn)在線應(yīng)用的要求。這意味著我們可能需要在模型精確度、求解速度和解決方案的魯棒性之間做出權(quán)衡。這套MISOCP方法論的價(jià)值不僅在于它能給出一個(gè)“最優(yōu)”的參考解更在于它為我們提供了一個(gè)嚴(yán)格、靈活的分析框架讓我們能夠系統(tǒng)地、量化地評(píng)估各種網(wǎng)絡(luò)配置方案的優(yōu)劣這是傳統(tǒng)經(jīng)驗(yàn)主義方法無法比擬的。