灰狼算法求解IEEE30節(jié)點(diǎn)環(huán)境經(jīng)濟(jì)調(diào)度：建模、實(shí)現(xiàn)與調(diào)優(yōu))
電力系統(tǒng)的環(huán)境經(jīng)濟(jì)調(diào)度Environmental Economic DispatchEED這幾年一直是研究熱點(diǎn)尤其是把多目標(biāo)灰狼優(yōu)化算法Multi-Objective Grey Wolf OptimizerMOGWO用到 IEEE30 節(jié)點(diǎn)系統(tǒng)上做仿真算是這個(gè)方向里一套比較標(biāo)準(zhǔn)的組合玩法。我的經(jīng)驗(yàn)是這類項(xiàng)目最核心的價(jià)值不在算法本身多新奇而在于你怎么把電力系統(tǒng)約束、目標(biāo)函數(shù)和多目標(biāo)優(yōu)化框架正確耦合在一起——很多人代碼跑通了結(jié)果卻是錯(cuò)的就是因?yàn)槟繕?biāo)函數(shù)建模不夠嚴(yán)謹(jǐn)。這篇文章我準(zhǔn)備把做這個(gè)項(xiàng)目的完整思路拆開來講從問題建模、算法原理到 Matlab 代碼實(shí)現(xiàn)、參數(shù)調(diào)試再到我踩過的坑和排查方法全部捋一遍。適合研究生、電力系統(tǒng)方向工程師以及想快速上手多目標(biāo)優(yōu)化應(yīng)用的開發(fā)者。1. 環(huán)境經(jīng)濟(jì)調(diào)度到底在優(yōu)化什么1.1 不止是最小化發(fā)電成本傳統(tǒng)經(jīng)濟(jì)調(diào)度只考慮一件事讓系統(tǒng)總發(fā)電成本最低。但現(xiàn)實(shí)里發(fā)電側(cè)還會(huì)產(chǎn)生污染物排放包括 CO2、SO2、NOx 等尤其以燃煤機(jī)組為主。環(huán)境經(jīng)濟(jì)調(diào)度就是在“成本最低”和“排放最少”這兩個(gè)目標(biāo)之間找平衡。這兩個(gè)目標(biāo)是沖突的。比如低成本的機(jī)組往往是煤耗高、排放大的老機(jī)組而清潔機(jī)組單位電價(jià)又偏高。你要讓成本降下去排放通常就上來了想讓排放降下去成本必然增加。所以 EED 本質(zhì)上是雙目標(biāo)優(yōu)化問題。形式化一點(diǎn)兩個(gè)目標(biāo)可以寫成成本目標(biāo)( \min f_1 \sum_{i1}^{n} (a_i P_i^2 b_i P_i c_i) )排放目標(biāo)( \min f_2 \sum_{i1}^{n} (\alpha_i P_i^2 \beta_i P_i \gamma_i) )其中 ( P_i ) 是第 i 臺(tái)發(fā)電機(jī)的有功出力( a_i、b_i、c_i ) 是成本系數(shù)( \alpha_i、\beta_i、\gamma_i ) 是排放系數(shù)。IEEE30 節(jié)點(diǎn)系統(tǒng)里通常取 6 臺(tái)發(fā)電機(jī)作為調(diào)度對(duì)象所以這個(gè)向量維度是 6。你別小看這個(gè)簡(jiǎn)單的二次函數(shù)形式實(shí)際做項(xiàng)目時(shí)還可以進(jìn)一步擴(kuò)展成帶閥點(diǎn)效應(yīng)的模型就是在成本函數(shù)里加一個(gè)正弦項(xiàng) ( |e_i \sin(f_i (P_i^{min} - P_i))| )。加了閥點(diǎn)效應(yīng)以后目標(biāo)函數(shù)變得非凸、多峰對(duì)優(yōu)化算法的要求立刻上了一個(gè)臺(tái)階。IEEE30 節(jié)點(diǎn)的經(jīng)典數(shù)據(jù)里其實(shí)有帶閥點(diǎn)和不帶閥點(diǎn)兩個(gè)版本建議新手從不帶閥點(diǎn)開始做先保證流程走通。1.2 等式約束和不等式約束一個(gè)都不能省目標(biāo)函數(shù)只是上半部分EED 的難點(diǎn)更多在下半部分——約束處理。最基本的約束有三類約束類型數(shù)學(xué)表達(dá)說明功率平衡( \sum P_i P_D P_{loss} )總出力必須等于負(fù)荷加網(wǎng)損機(jī)組出力上下限( P_i^{min} \le P_i \le P_i^{max} )每臺(tái)機(jī)組都有物理出力邊界機(jī)組爬坡約束可選( -P_{ri} \le P_i - P_i^0 \le P_{ri} )相鄰時(shí)段出力變化受限這里容易出問題的是網(wǎng)損。IEEE30 節(jié)點(diǎn)經(jīng)濟(jì)調(diào)度文獻(xiàn)里最常用的是 B 系數(shù)法也就是 ( P_{loss} \sum_i \sum_j B_{ij} P_i P_j \sum_i B_{i0} P_i B_{00} )。網(wǎng)損系數(shù)是固定矩陣在標(biāo)準(zhǔn)測(cè)試數(shù)據(jù)里可以直接查到。如果做靜態(tài) EED 不考慮爬坡約束代碼會(huì)省很多事。做動(dòng)態(tài)調(diào)度多時(shí)段就必須加上。我的建議是第一版先把功率平衡、出力上下限和網(wǎng)損做好這部分已經(jīng)能出漂亮的帕累托前沿了后面再考慮動(dòng)態(tài)場(chǎng)景。1.3 為什么選擇 IEEE30 節(jié)點(diǎn)系統(tǒng)IEEE30 節(jié)點(diǎn)是一個(gè)標(biāo)準(zhǔn)測(cè)試系統(tǒng)30 條母線、6 臺(tái)發(fā)電機(jī)、41 條支路。你隨便打開一篇經(jīng)濟(jì)調(diào)度論文十篇里至少有七篇是用它做驗(yàn)證的。原因很簡(jiǎn)單數(shù)據(jù)公開、規(guī)模適中、結(jié)果可對(duì)比。規(guī)模適中這一點(diǎn)特別重要。6 臺(tái)發(fā)電機(jī)意味著決策變量是 6 維這個(gè)維度對(duì)群智能優(yōu)化算法來說既不算簡(jiǎn)單到無聊又不至于難到算不動(dòng)。如果換成 IEEE118 節(jié)點(diǎn)系統(tǒng)54 臺(tái)機(jī)組很多基礎(chǔ)版算法收斂性就很差了。所以用 IEEE30 節(jié)點(diǎn)做算法驗(yàn)證是學(xué)術(shù)界默認(rèn)的“標(biāo)準(zhǔn)化實(shí)驗(yàn)場(chǎng)”。另外多目標(biāo)優(yōu)化的結(jié)果好不好得有一個(gè)參照。IEEE30 節(jié)點(diǎn)系統(tǒng)有很多公開的經(jīng)典帕累托前沿記錄你可以拿自己的 MOGWO 跑出來的結(jié)果和參考文獻(xiàn)里的 NSGA-II、MOPSO 結(jié)果做對(duì)比這是很有說服力的驗(yàn)證方式。2. 多目標(biāo)灰狼優(yōu)化算法的核心機(jī)制2.1 從灰狼捕獵到數(shù)學(xué)搜索灰狼優(yōu)化算法是模仿灰狼群體捕獵行為的一種群智能算法?；依侨后w有嚴(yán)格的社會(huì)等級(jí)頭狼、副頭狼、普通狼和底層狼。在算法里我們把這些角色抽象成“最優(yōu)解”“次優(yōu)解”“第三優(yōu)解”和“候選解”。每一頭狼的位置就是一個(gè)候選解在你的 EED 問題里就是一個(gè) 6 維向量代表 6 臺(tái)機(jī)組的出力。狼群捕獵過程抽象成三個(gè)步驟圍捕根據(jù)頭狼位置更新位置追捕根據(jù)前三頭狼的位置綜合確定移動(dòng)方向攻擊用參數(shù) a 控制搜索范圍逐漸縮小核心更新公式是[ \vec{D} | \vec{C} \cdot \vec{X}_p(t) - \vec{X}(t) | ] [ \vec{X}(t1) \vec{X}_p(t) - \vec{A} \cdot \vec{D} ]其中 ( \vec{A} 2\vec{a}\cdot\vec{r}_1 - \vec{a} )( \vec{C} 2\vec{r}_2 )。參數(shù) a 從 2 線性降到 0決定了算法是偏向全局搜索還是局部開發(fā)前期 a 大狼群分散探索后期 a 小集中在最優(yōu)解附近精細(xì)搜索。這里有個(gè)容易忽略的細(xì)節(jié)標(biāo)準(zhǔn) GWO 是單目標(biāo)算法它假設(shè)只有一個(gè)全局最優(yōu)解頭狼是唯一明確的指引。但是多目標(biāo)問題里不存在“唯一最優(yōu)解”而是一組互不支配的候選解集合也就是帕累托前沿。所以 MOGWO 不能直接套用標(biāo)準(zhǔn) GWO 的框架。2.2 多目標(biāo)化改造的幾個(gè)關(guān)鍵設(shè)計(jì)MOGWO 在多目標(biāo)方向上做了三個(gè)關(guān)鍵改進(jìn)第一引入外部檔案Repository。這個(gè)檔案用來存放目前已找到的非支配解。每次迭代后把當(dāng)前狼群中的非支配解插入檔案并且去掉被新解支配的舊解。檔案大小有限如果滿了就用擁擠距離策略刪除最擁擠區(qū)域的解保解的多樣性。第二領(lǐng)導(dǎo)者選擇機(jī)制。標(biāo)準(zhǔn) GWO 里只有三頭主狼但多目標(biāo)下需要從檔案里選。MOGWO 的做法是把檔案按照目標(biāo)空間劃分成網(wǎng)格每個(gè)網(wǎng)格記錄解的密度然后從低密度區(qū)域優(yōu)先選擇領(lǐng)導(dǎo)者這樣可以讓搜索方向偏向尚未充分探索的區(qū)域。第三位置更新后的邊界處理。更新完新位置后決策變量要拉回可行域。很多實(shí)現(xiàn)直接寫死 min-max 裁剪但更好的做法是針對(duì)越界的變量重新在邊界鄰域內(nèi)隨機(jī)初始化保持種群多樣性。我看過不少直接用單目標(biāo) GWO 代碼改寫的“假 MOGWO”最常見的問題就是外部檔案只有插入沒有刪除或者領(lǐng)導(dǎo)者選擇不考慮密集度結(jié)果就是所有解堆到帕累托前沿的一小段上前沿兩端完全空掉。這種事不是跑不跑得通的問題而是結(jié)果有沒有參考價(jià)值的問題。2.3 MOGWO 的參數(shù)敏感度分析多目標(biāo)優(yōu)化里的參數(shù)設(shè)置直接影響收斂性和解分布質(zhì)量。MOGWO 主要參數(shù)包括參數(shù)推薦值影響種群規(guī)模30-50太小收斂不穩(wěn)太大計(jì)算量大迭代次數(shù)100-300看目標(biāo)函數(shù)復(fù)雜度帶閥點(diǎn)建議加大外部檔案大小50-100決定最終前沿點(diǎn)數(shù)網(wǎng)格劃分?jǐn)?shù)10-20影響領(lǐng)導(dǎo)者選擇多樣性我測(cè)過 IEEE30 系統(tǒng)的 EED 問題種群 40、迭代 200、檔案 80、網(wǎng)格 10 是比較穩(wěn)的組合。再往上加大迭代次數(shù)的收益率遞減因?yàn)槟繕?biāo)函數(shù)比較簡(jiǎn)單200 代之后帕累托前沿基本不變化。有一點(diǎn)要注意多目標(biāo)優(yōu)化的評(píng)價(jià)指標(biāo)不只是畫個(gè)圖IGD反轉(zhuǎn)世代距離、HV超體積、Spread分布度是三個(gè)最常用的量化指標(biāo)。一個(gè)嚴(yán)謹(jǐn)?shù)膶?duì)比實(shí)驗(yàn)至少算 HV 和 IGD跑個(gè) 30 次獨(dú)立實(shí)驗(yàn)取均值和方差才能說明算法穩(wěn)定性。很多新手只看一次跑的圖和最小成本就下結(jié)論這種方法論文里站不住腳。3. Matlab 實(shí)現(xiàn)環(huán)境經(jīng)濟(jì)調(diào)度的完整過程3.1 數(shù)據(jù)準(zhǔn)備IEEE30 節(jié)點(diǎn)參數(shù)與 B 矩陣在 Matlab 里做這個(gè)項(xiàng)目第一步是準(zhǔn)備數(shù)據(jù)文件。我習(xí)慣把機(jī)組參數(shù)、負(fù)荷、B 系數(shù)做成三個(gè)獨(dú)立的數(shù)據(jù)文件或者腳本避免全堆在一個(gè)主文件里數(shù)據(jù)項(xiàng)多了以后極其難排查。機(jī)組數(shù)據(jù)主要包括6 臺(tái)發(fā)電機(jī)的有功出力上下限、成本系數(shù)和排放系數(shù)。如果你用標(biāo)準(zhǔn) IEEE30 系統(tǒng)數(shù)據(jù)會(huì)看到類似這樣的一組典型值我隨手寫一組常用值做示例實(shí)際做實(shí)驗(yàn)請(qǐng)以你下載的 paper 數(shù)據(jù)為準(zhǔn)機(jī)組Pmin(MW)Pmax(MW)abcalphabetagammaG1502000.003752.0000.01234-0.10045.098G220800.017501.7500.00823-0.06052.407G315500.062501.0000.00418-0.05031.550G410350.008343.2500.00308-0.02400.806G510300.025003.0000.00534-0.05301.781G612400.025003.0000.00454-0.03901.348B 矩陣是 7x7 的對(duì)稱矩陣包含 B、B0、B00 三個(gè)分量。在 Matlab 里我建議存成結(jié)構(gòu)體sys.B [...]; % 6x6 二次項(xiàng)系數(shù) sys.B0 [...]; % 1x6 一次項(xiàng)系數(shù) sys.B00 ...; % 常數(shù)項(xiàng)網(wǎng)損計(jì)算公式function Ploss calc_loss(P, sys) n length(P); Ploss P. * sys.B * P sys.B0 * P. sys.B00; % 注意維度對(duì)應(yīng) end寫這一行網(wǎng)損計(jì)算時(shí)矩陣維度最容易搞反。P 是 6x1 列向量sys.B 是 6x6( P^T B P ) 得到標(biāo)量。sys.B0 我在存儲(chǔ)時(shí)寫成 6x1 列向量所以用 ( B0^T P ) 得到標(biāo)量。這種維度問題只有運(yùn)行起來報(bào)錯(cuò)或者算出來的值明顯不對(duì)才會(huì)發(fā)現(xiàn)。3.2 目標(biāo)函數(shù)與約束處理代碼怎么寫主目標(biāo)函數(shù)文件我建議輸入是 6 維出力向量輸出是成本和排放兩個(gè)標(biāo)量。代碼大致長這樣function [f_cost, f_emit] EED_objective(P, sys) % P: 決策變量向量 (nGen x 1) % f_cost: 總發(fā)電成本 % f_emit: 總污染物排放 P P(:); n length(P); % 成本函數(shù)二次型 f_cost sum(sys.a .* P.^2 sys.b .* P sys.c); % 排放函數(shù)二次型 f_emit sum(sys.alpha .* P.^2 sys.beta .* P sys.gamma); end約束處理是重頭戲。MOGWO 不像罰函數(shù)法那樣簡(jiǎn)單也不是所有越界都直接丟棄。我的做法是這樣出力上下限直接在更新位置后做 min-max 裁剪功率平衡約束處理完上下限之后檢查總出力、總負(fù)荷和網(wǎng)損之和的差但直接裁剪會(huì)產(chǎn)生一個(gè)問題6 個(gè)變量裁剪完總出力不等于負(fù)荷加網(wǎng)損等式約束被破壞。這里我用啟發(fā)式修正也叫功率平衡修復(fù)機(jī)制% 計(jì)算負(fù)荷和網(wǎng)損 total_demand sys.Pd; Ploss calc_loss(P, sys); error sum(P) - total_demand - Ploss; % 把誤差均攤到可調(diào)機(jī)組上 adjust error / n; for i 1:n P_new(i) P(i) - adjust; if P_new(i) sys.Pmin(i) P_new(i) sys.Pmin(i); elseif P_new(i) sys.Pmax(i) P_new(i) sys.Pmax(i); end end這種均攤修正不是完美的可能修完依然不滿足等式約束但誤差已經(jīng)很小。嚴(yán)格的做法是把等式約束轉(zhuǎn)成等式罰函數(shù)項(xiàng)加進(jìn)目標(biāo)函數(shù)里權(quán)衡一下我的經(jīng)驗(yàn)是對(duì)于 IEEE30 靜態(tài)調(diào)度啟發(fā)式修復(fù)加少量罰函數(shù)混合處理效果和精確求解差不多但速度快很多。3.3 MOGWO 主循環(huán)框架搭建我直接給一個(gè)精簡(jiǎn)版的 MOGWO 主循環(huán)框架你可以在此基礎(chǔ)上改寫成完整代碼% 初始化狼群 X init_wolf(N, sys); % N 只狼每個(gè)是6維向量 [f1, f2] evaluate_wolf(X, sys); repo init_repository(X, f1, f2); % 初始化外部檔案 for t 1:MaxIt a 2 - 2 * t / MaxIt; for i 1:N % 從檔案里選擇三個(gè)領(lǐng)導(dǎo)者基于網(wǎng)格密度 alpha select_leader(repo); beta select_leader(repo); delta select_leader(repo); % 按標(biāo)準(zhǔn)GWO公式更新位置對(duì)三個(gè)領(lǐng)導(dǎo)者分別計(jì)算再平均 X(i, :) update_position(X(i, :), alpha, beta, delta, a); % 邊界處理 X(i, :) clamp(X(i, :), sys.Pmin, sys.Pmax); % 啟發(fā)式功率平衡修復(fù) X(i, :) power_balance_fix(X(i, :), sys); end % 重新評(píng)估并更新檔案 [f1, f2] evaluate_wolf(X, sys); repo update_repository(repo, X, f1, f2); end這是結(jié)構(gòu)最穩(wěn)定的算法骨架。關(guān)鍵在于 update_repository 函數(shù)要寫好它總共做三件事合并新舊解、移除被支配解、如果檔案溢出則基于網(wǎng)格密度刪解。3.4 結(jié)果可視化與性能指標(biāo)計(jì)算跑完循環(huán)后你要輸出三樣?xùn)|西帕累托前沿圖橫軸成本、縱軸排放這通常是一張往下凸的散點(diǎn)圖。figure; plot(repo.cost, repo.emit, o, MarkerSize, 4); xlabel(Cost ($/h)); ylabel(Emission (ton/h)); grid on;最小成本解在檔案里找成本最小的方案輸出對(duì)應(yīng)出力向量。[min_cost, idx] min(repo.cost); best_X repo.X(idx, :); disp([Minimum Cost: , num2str(min_cost)]);我在實(shí)際跑的時(shí)候還會(huì)計(jì)算 HV 和 IGD 兩個(gè)指標(biāo)尤其是做算法對(duì)比時(shí)光靠圖無法說明問題。HV 計(jì)算需要一個(gè)參考點(diǎn)一般取所有解中每個(gè)目標(biāo)的最大值稍微放大一點(diǎn)然后算帕累托前沿和參考點(diǎn)之間覆蓋的超體積。寫指標(biāo)計(jì)算時(shí)可以考慮用第三方工具箱也可以手寫。如果只是為了項(xiàng)目驗(yàn)證我建議先手寫 IGD——邏輯上就是真實(shí)前沿上每個(gè)點(diǎn)到你算出的前沿的最近距離的平均值代碼不復(fù)雜對(duì)你有好處數(shù)據(jù)一致性也更可控。4. 結(jié)果分析與超參數(shù)調(diào)優(yōu)的實(shí)戰(zhàn)經(jīng)驗(yàn)4.1 帕累托前沿圖的解讀方式跑完 MOGWO 之后你看到的第一張圖是帕累托前沿。正常的 EED 前沿應(yīng)該是平滑下降的凸曲線左端點(diǎn)成本最低但排放最高右端點(diǎn)排放最低但成本最高兩端之間的點(diǎn)就是成本和排放的不同折中方案。如果你的圖出現(xiàn)以下情況大概率有問題前沿只覆蓋很小一段區(qū)間說明算法多樣性差領(lǐng)導(dǎo)者選擇策略失效前沿分布不連續(xù)、有空洞說明檔案更新策略不好某些網(wǎng)格區(qū)域沒有解被保留多數(shù)點(diǎn)擠在一起集中在某一側(cè)說明目標(biāo)歸一化沒做好或檔案溢出刪除策略偏頗我在調(diào)這種問題時(shí)通常優(yōu)先檢查外部檔案的網(wǎng)格劃分和 update_repository 里的種子選擇邏輯。有時(shí)候問題很簡(jiǎn)單是檔案大小設(shè)太小導(dǎo)致前沿點(diǎn)不夠密。4.2 單目標(biāo) vs 多目標(biāo)對(duì)比實(shí)驗(yàn)一個(gè)常見的加分做法是把 MOGWO 結(jié)果和單目標(biāo) GWO 結(jié)果對(duì)比。你可以用標(biāo)準(zhǔn) GWO 分別跑“最小成本”和“最小排放”兩個(gè)任務(wù)得到兩個(gè)極端解然后驗(yàn)證 MOGWO 前沿的兩個(gè)端點(diǎn)是否接近這兩個(gè)單目標(biāo)最優(yōu)解。這一步驗(yàn)證意義很大。如果 MOGWO 跑出的左端點(diǎn)成本比單目標(biāo) GWO 算出的最小成本還高 5%那么說明 MOGWO 的收斂性還有問題——極端解沒有被充分搜索到。我碰到過這種情況原因是種群規(guī)模太小且領(lǐng)導(dǎo)者選擇時(shí)過度偏向低密度區(qū)域?qū)е氯肿顑?yōu)方向被忽略。調(diào)整方向有兩個(gè)一是增大種群二是把 a 的下降速率調(diào)慢一點(diǎn)增強(qiáng)后期局部搜索。比如把 a 從線性衰減改成非線性衰減a 2 * (1 - (t / MaxIt)^0.7);這種非線性衰減在前期保持了更多探索后期收斂更精細(xì)。我實(shí)測(cè)在 IEEE30 EED 問題上改用平方根形式的衰減后最小成本端大約改善了 1.8%前沿分布也更均勻。4.3 多目標(biāo)評(píng)價(jià)指標(biāo)的使用心得做任何多目標(biāo)優(yōu)化項(xiàng)目建議至少要算三個(gè)指標(biāo)指標(biāo)全稱度量?jī)?nèi)容HVHypervolume解集覆蓋目標(biāo)空間的大小越大越好IGDInverted Generational Distance真實(shí)前沿到計(jì)算前沿的距離越小越好Spread分布度前沿點(diǎn)分布的均勻性越小越好需要注意的是HV 和 IGD 都需要一個(gè)“真實(shí)帕累托前沿”或者參考點(diǎn)。真實(shí)前沿怎么來常見方法是用一個(gè)超大種群、超多迭代次數(shù)的 MOGWO 跑出來的前沿作為近似真實(shí)前沿。記住這是“近似”但你用它做相對(duì)比較已經(jīng)足夠。我習(xí)慣每次實(shí)驗(yàn)獨(dú)立跑 30 次然后報(bào)告平均值±標(biāo)準(zhǔn)差。為什么是 30 次因?yàn)槿褐悄芩惴ㄓ须S機(jī)性單次結(jié)果沒有統(tǒng)計(jì)意義。在論文里寫“最好結(jié)果”是常見但不太嚴(yán)謹(jǐn)?shù)淖龇ㄎ医ㄗh至少同時(shí)報(bào)告均值和最差值這樣能體現(xiàn)算法的魯棒性。5. 實(shí)操中踩過的坑與排查技巧5.1 結(jié)果不收斂或前沿不完整的常見原因這一類問題我見得最多基本集中在四個(gè)原因第一個(gè)原因是初始種群生成不合法。如果初始狼群里很多個(gè)體的總出力與負(fù)荷和網(wǎng)損差的太遠(yuǎn)功率平衡修復(fù)函數(shù)會(huì)花很多代去修正消耗搜索資源。解決辦法是把初始種群生成時(shí)就添加一定比例的“合法個(gè)體”——也就是先產(chǎn)生隨機(jī)出力再通過啟發(fā)式修正讓部分個(gè)體滿足功率平衡讓種群起點(diǎn)就有一部分在可行域里。第二個(gè)原因是線性衰減的 a 參數(shù)不夠好。不同問題的收斂速度區(qū)別很大如果 a 降到零之后狼群還在大規(guī)模跳躍搜索就會(huì)在某代之后無法收斂。檢查辦法是畫出歷代的 HV 變化曲線如果后期 HV 還在劇烈波動(dòng)說明 a 衰減過快如果 HV 早早平穩(wěn)但最終值偏低說明衰減過慢。第三個(gè)原因是目標(biāo)函數(shù)數(shù)值跨度太大會(huì)影響網(wǎng)格劃分。比如成本是 600-800排放是 0.1-0.3如果不做歸一化就直接劃分網(wǎng)格網(wǎng)格在排放維度上會(huì)特別寬在成本維度上特別窄導(dǎo)致檔案選擇偏向某一個(gè)目標(biāo)。我會(huì)在計(jì)算網(wǎng)格密度之前對(duì)兩個(gè)目標(biāo)做 max-min 歸一化這是很容易被忽略但又非常關(guān)鍵的一步。第四個(gè)原因是外部檔案的更新時(shí)機(jī)不對(duì)。有些實(shí)現(xiàn)是等所有狼都更新完才刷新檔案有些是每只狼更新完就刷新。兩種方式各有優(yōu)劣但在 Matlab 的向量化實(shí)現(xiàn)里如果排序和更新邏輯寫錯(cuò)很容易造成檔案里出現(xiàn)跑出邊界外的臟數(shù)據(jù)。建議每次更新檔案后都打印一行檔案數(shù)量和邊界統(tǒng)計(jì)信息肉眼檢查是否出現(xiàn)丟失維度問題。5.2 Matlat 性能與代碼向量化MOGWO 在 IEEE30 上規(guī)模很小本來不用太關(guān)心性能。但如果你把種群調(diào)到 100、迭代 300 次純 for 循環(huán)跑到后面也會(huì)覺得卡。一個(gè)最簡(jiǎn)單有效的優(yōu)化是主循環(huán)盡量向量化評(píng)估函數(shù)一次性接收整個(gè)種群矩陣。function [F1, F2] evaluate_wolf_vectorized(X, sys) % X: N x nGen 矩陣每行是一個(gè)個(gè)體 F1 sum(sys.a .* X.^2 sys.b .* X sys.c, 2); F2 sum(sys.alpha .* X.^2 sys.beta .* X sys.gamma, 2); end相比于逐個(gè)個(gè)體調(diào)函數(shù)向量化至少能快 5 到 10 倍。如果你的問題擴(kuò)展到了動(dòng)態(tài)調(diào)度或者多時(shí)段代碼向量化的收益會(huì)更明顯。另一個(gè)實(shí)用技巧是提前為每個(gè)個(gè)體緩存目標(biāo)函數(shù)值。同一代內(nèi)同一個(gè)個(gè)體因?yàn)樾迯?fù)機(jī)制可能會(huì)被重復(fù)評(píng)估如果不緩存你就是在做無效計(jì)算。5.3 結(jié)果復(fù)現(xiàn)與實(shí)驗(yàn)記錄管理學(xué)術(shù)項(xiàng)目里最怕的就是跑完一組實(shí)驗(yàn)結(jié)果挺好的但過兩天想復(fù)現(xiàn)發(fā)現(xiàn)參數(shù)忘了。所以我建議從一開始就養(yǎng)成記錄實(shí)驗(yàn)配置的習(xí)慣種群數(shù)、迭代數(shù)、檔案大小、a 衰減方式、修復(fù)機(jī)制開關(guān)、隨機(jī)種子這些全部存到一個(gè)結(jié)構(gòu)體里輸出結(jié)果時(shí)同時(shí)輸出配置。exp_config struct(N, 40, MaxIt, 200, ... nGrid, 10, RepoSize, 80, ... alpha_decay, nonlinear_0.7, randseed, 1); rng(exp_config.randseed);保存結(jié)果時(shí)把 exp_config 一起存到 .mat 文件里這樣后續(xù)即使代碼改了你還能知道當(dāng)時(shí)的結(jié)果是怎么來的。這個(gè)習(xí)慣我吃過大虧以后養(yǎng)成的建議你直接學(xué)走。5.4 和 NSGA-II、MOPSO 對(duì)比時(shí)的公平性問題如果你打算在論文里和 NSGA-II、MOPSO 做對(duì)比一定要注意一個(gè)公平性原則所有算法使用相同的函數(shù)評(píng)估次數(shù)而不是相同的迭代次數(shù)。種群 40 迭代 200 的 MOGWO 總函數(shù)評(píng)估次數(shù)是 40×200 次那么 NSGA-II 也應(yīng)該用相同的總評(píng)估次數(shù)。這一點(diǎn)大家經(jīng)常忽略對(duì)比結(jié)果就失去了說服力。我曾經(jīng)見過一個(gè)實(shí)驗(yàn)MOGWO 迭代次數(shù)是 NSGA-II 的兩倍然后得出結(jié)論說 MOGWO 好這個(gè)結(jié)論就沒意義。公平對(duì)比還有個(gè)細(xì)節(jié)初始化方式要一致都用相同隨機(jī)種子和相同風(fēng)格的初始種群不能單獨(dú)為一個(gè)算法精心設(shè)計(jì)初始化。6. 擴(kuò)展方向與個(gè)人經(jīng)驗(yàn)分享6.1 加入閥點(diǎn)效應(yīng)與斜坡約束我前文提到的帶閥點(diǎn)成本函數(shù)建議在做完基本版后立刻補(bǔ)上。把成本函數(shù)改成[ F_i(P_i) a_i P_i^2 b_i P_i c_i |e_i \sin(f_i (P_i^{min} - P_i))| ]加了這個(gè)正弦項(xiàng)目標(biāo)函數(shù)變成非凸多峰MOGWO 的收斂性會(huì)明顯受到挑戰(zhàn)。這時(shí)候你再檢驗(yàn)算法效果才有真正的說服力。動(dòng)態(tài)調(diào)度再加爬坡約束和燃料約束問題就完全展開了。6.2 增加第三目標(biāo)除了成本和排放還可以把網(wǎng)損作為第三目標(biāo)最小化或者把系統(tǒng)電壓偏差加入目標(biāo)。多目標(biāo)數(shù)量增加后有兩個(gè)問題一是帕累托前沿變成平面或曲面可視化難度增加二是群體優(yōu)化算法的性能會(huì)下降因?yàn)榉侵浣鈹?shù)量爆炸檔案更新壓力變大。這時(shí)候通常需要引入基于分解的機(jī)制比如 MOEA/D 的思想或者改進(jìn)檔案維護(hù)策略。我的建議是第三目標(biāo)可以在同一套 MOGWO 框架下直接加但要注意網(wǎng)格劃分從二維變?nèi)S代碼的復(fù)雜度呈指數(shù)增長。初學(xué)階段先加“網(wǎng)損”而不是“電壓偏差”因?yàn)榫W(wǎng)損計(jì)算成本低電壓偏差需要潮流計(jì)算算力負(fù)擔(dān)完全不同。6.3 混合改進(jìn)方向MOGWO 的一個(gè)常見改進(jìn)思路是差分進(jìn)化混合在每次更新后用差分進(jìn)化算子的變異交叉代替標(biāo)準(zhǔn)位置更新增強(qiáng)開發(fā)能力。另一個(gè)思路是引入局部搜索對(duì)檔案中的非支配解做高斯擾動(dòng)和重新評(píng)估提高前沿的收斂精度。做混合改進(jìn)時(shí)要注意不要為了改進(jìn)而改進(jìn)。每一次改動(dòng)都要有明確的邏輯和實(shí)驗(yàn)驗(yàn)證。先跑標(biāo)準(zhǔn) MOGWO 算出基線結(jié)果然后記錄它的薄弱點(diǎn)比如哪個(gè)目標(biāo)收斂差、哪段前沿分布稀疏再有針對(duì)性地增強(qiáng)這種“問題驅(qū)動(dòng)”的改進(jìn)思路才是科研的正確路徑。從我自己的項(xiàng)目經(jīng)驗(yàn)來看做這種優(yōu)化調(diào)度類項(xiàng)目最忌一步到位。先把最簡(jiǎn)單的 IEEE30 靜態(tài)雙目標(biāo)跑通完整地理解目標(biāo)函數(shù)、約束和 MOGWO 的每一個(gè)機(jī)制再逐步擴(kuò)展。前面基礎(chǔ)打得扎實(shí)后面加約束、加目標(biāo)、換算法都會(huì)順手很多。最后提醒一句跑任何實(shí)驗(yàn)都記得保存隨機(jī)種子和配置參數(shù)別等出了好結(jié)果再去追根溯源那時(shí)候代價(jià)太大了。 多目標(biāo)灰狼算法求解IEEE30節(jié)點(diǎn)環(huán)境經(jīng)濟(jì)調(diào)度：建模、實(shí)現(xiàn)與調(diào)優(yōu)
電力系統(tǒng)的環(huán)境經(jīng)濟(jì)調(diào)度Environmental Economic DispatchEED這幾年一直是研究熱點(diǎn)尤其是把多目標(biāo)灰狼優(yōu)化算法Multi-Objective Grey Wolf OptimizerMOGWO用到 IEEE30 節(jié)點(diǎn)系統(tǒng)上做仿真算是這個(gè)方向里一套比較標(biāo)準(zhǔn)的組合玩法。我的經(jīng)驗(yàn)是這類項(xiàng)目最核心的價(jià)值不在算法本身多新奇而在于你怎么把電力系統(tǒng)約束、目標(biāo)函數(shù)和多目標(biāo)優(yōu)化框架正確耦合在一起——很多人代碼跑通了結(jié)果卻是錯(cuò)的就是因?yàn)槟繕?biāo)函數(shù)建模不夠嚴(yán)謹(jǐn)。這篇文章我準(zhǔn)備把做這個(gè)項(xiàng)目的完整思路拆開來講從問題建模、算法原理到 Matlab 代碼實(shí)現(xiàn)、參數(shù)調(diào)試再到我踩過的坑和排查方法全部捋一遍。適合研究生、電力系統(tǒng)方向工程師以及想快速上手多目標(biāo)優(yōu)化應(yīng)用的開發(fā)者。1. 環(huán)境經(jīng)濟(jì)調(diào)度到底在優(yōu)化什么1.1 不止是最小化發(fā)電成本傳統(tǒng)經(jīng)濟(jì)調(diào)度只考慮一件事讓系統(tǒng)總發(fā)電成本最低。但現(xiàn)實(shí)里發(fā)電側(cè)還會(huì)產(chǎn)生污染物排放包括 CO2、SO2、NOx 等尤其以燃煤機(jī)組為主。環(huán)境經(jīng)濟(jì)調(diào)度就是在“成本最低”和“排放最少”這兩個(gè)目標(biāo)之間找平衡。這兩個(gè)目標(biāo)是沖突的。比如低成本的機(jī)組往往是煤耗高、排放大的老機(jī)組而清潔機(jī)組單位電價(jià)又偏高。你要讓成本降下去排放通常就上來了想讓排放降下去成本必然增加。所以 EED 本質(zhì)上是雙目標(biāo)優(yōu)化問題。形式化一點(diǎn)兩個(gè)目標(biāo)可以寫成成本目標(biāo)( \min f_1 \sum_{i1}^{n} (a_i P_i^2 b_i P_i c_i) )排放目標(biāo)( \min f_2 \sum_{i1}^{n} (\alpha_i P_i^2 \beta_i P_i \gamma_i) )其中 ( P_i ) 是第 i 臺(tái)發(fā)電機(jī)的有功出力( a_i、b_i、c_i ) 是成本系數(shù)( \alpha_i、\beta_i、\gamma_i ) 是排放系數(shù)。IEEE30 節(jié)點(diǎn)系統(tǒng)里通常取 6 臺(tái)發(fā)電機(jī)作為調(diào)度對(duì)象所以這個(gè)向量維度是 6。你別小看這個(gè)簡(jiǎn)單的二次函數(shù)形式實(shí)際做項(xiàng)目時(shí)還可以進(jìn)一步擴(kuò)展成帶閥點(diǎn)效應(yīng)的模型就是在成本函數(shù)里加一個(gè)正弦項(xiàng) ( |e_i \sin(f_i (P_i^{min} - P_i))| )。加了閥點(diǎn)效應(yīng)以后目標(biāo)函數(shù)變得非凸、多峰對(duì)優(yōu)化算法的要求立刻上了一個(gè)臺(tái)階。IEEE30 節(jié)點(diǎn)的經(jīng)典數(shù)據(jù)里其實(shí)有帶閥點(diǎn)和不帶閥點(diǎn)兩個(gè)版本建議新手從不帶閥點(diǎn)開始做先保證流程走通。1.2 等式約束和不等式約束一個(gè)都不能省目標(biāo)函數(shù)只是上半部分EED 的難點(diǎn)更多在下半部分——約束處理。最基本的約束有三類約束類型數(shù)學(xué)表達(dá)說明功率平衡( \sum P_i P_D P_{loss} )總出力必須等于負(fù)荷加網(wǎng)損機(jī)組出力上下限( P_i^{min} \le P_i \le P_i^{max} )每臺(tái)機(jī)組都有物理出力邊界機(jī)組爬坡約束可選( -P_{ri} \le P_i - P_i^0 \le P_{ri} )相鄰時(shí)段出力變化受限這里容易出問題的是網(wǎng)損。IEEE30 節(jié)點(diǎn)經(jīng)濟(jì)調(diào)度文獻(xiàn)里最常用的是 B 系數(shù)法也就是 ( P_{loss} \sum_i \sum_j B_{ij} P_i P_j \sum_i B_{i0} P_i B_{00} )。網(wǎng)損系數(shù)是固定矩陣在標(biāo)準(zhǔn)測(cè)試數(shù)據(jù)里可以直接查到。如果做靜態(tài) EED 不考慮爬坡約束代碼會(huì)省很多事。做動(dòng)態(tài)調(diào)度多時(shí)段就必須加上。我的建議是第一版先把功率平衡、出力上下限和網(wǎng)損做好這部分已經(jīng)能出漂亮的帕累托前沿了后面再考慮動(dòng)態(tài)場(chǎng)景。1.3 為什么選擇 IEEE30 節(jié)點(diǎn)系統(tǒng)IEEE30 節(jié)點(diǎn)是一個(gè)標(biāo)準(zhǔn)測(cè)試系統(tǒng)30 條母線、6 臺(tái)發(fā)電機(jī)、41 條支路。你隨便打開一篇經(jīng)濟(jì)調(diào)度論文十篇里至少有七篇是用它做驗(yàn)證的。原因很簡(jiǎn)單數(shù)據(jù)公開、規(guī)模適中、結(jié)果可對(duì)比。規(guī)模適中這一點(diǎn)特別重要。6 臺(tái)發(fā)電機(jī)意味著決策變量是 6 維這個(gè)維度對(duì)群智能優(yōu)化算法來說既不算簡(jiǎn)單到無聊又不至于難到算不動(dòng)。如果換成 IEEE118 節(jié)點(diǎn)系統(tǒng)54 臺(tái)機(jī)組很多基礎(chǔ)版算法收斂性就很差了。所以用 IEEE30 節(jié)點(diǎn)做算法驗(yàn)證是學(xué)術(shù)界默認(rèn)的“標(biāo)準(zhǔn)化實(shí)驗(yàn)場(chǎng)”。另外多目標(biāo)優(yōu)化的結(jié)果好不好得有一個(gè)參照。IEEE30 節(jié)點(diǎn)系統(tǒng)有很多公開的經(jīng)典帕累托前沿記錄你可以拿自己的 MOGWO 跑出來的結(jié)果和參考文獻(xiàn)里的 NSGA-II、MOPSO 結(jié)果做對(duì)比這是很有說服力的驗(yàn)證方式。2. 多目標(biāo)灰狼優(yōu)化算法的核心機(jī)制2.1 從灰狼捕獵到數(shù)學(xué)搜索灰狼優(yōu)化算法是模仿灰狼群體捕獵行為的一種群智能算法?；依侨后w有嚴(yán)格的社會(huì)等級(jí)頭狼、副頭狼、普通狼和底層狼。在算法里我們把這些角色抽象成“最優(yōu)解”“次優(yōu)解”“第三優(yōu)解”和“候選解”。每一頭狼的位置就是一個(gè)候選解在你的 EED 問題里就是一個(gè) 6 維向量代表 6 臺(tái)機(jī)組的出力。狼群捕獵過程抽象成三個(gè)步驟圍捕根據(jù)頭狼位置更新位置追捕根據(jù)前三頭狼的位置綜合確定移動(dòng)方向攻擊用參數(shù) a 控制搜索范圍逐漸縮小核心更新公式是[ \vec{D} | \vec{C} \cdot \vec{X}_p(t) - \vec{X}(t) | ] [ \vec{X}(t1) \vec{X}_p(t) - \vec{A} \cdot \vec{D} ]其中 ( \vec{A} 2\vec{a}\cdot\vec{r}_1 - \vec{a} )( \vec{C} 2\vec{r}_2 )。參數(shù) a 從 2 線性降到 0決定了算法是偏向全局搜索還是局部開發(fā)前期 a 大狼群分散探索后期 a 小集中在最優(yōu)解附近精細(xì)搜索。這里有個(gè)容易忽略的細(xì)節(jié)標(biāo)準(zhǔn) GWO 是單目標(biāo)算法它假設(shè)只有一個(gè)全局最優(yōu)解頭狼是唯一明確的指引。但是多目標(biāo)問題里不存在“唯一最優(yōu)解”而是一組互不支配的候選解集合也就是帕累托前沿。所以 MOGWO 不能直接套用標(biāo)準(zhǔn) GWO 的框架。2.2 多目標(biāo)化改造的幾個(gè)關(guān)鍵設(shè)計(jì)MOGWO 在多目標(biāo)方向上做了三個(gè)關(guān)鍵改進(jìn)第一引入外部檔案Repository。這個(gè)檔案用來存放目前已找到的非支配解。每次迭代后把當(dāng)前狼群中的非支配解插入檔案并且去掉被新解支配的舊解。檔案大小有限如果滿了就用擁擠距離策略刪除最擁擠區(qū)域的解保解的多樣性。第二領(lǐng)導(dǎo)者選擇機(jī)制。標(biāo)準(zhǔn) GWO 里只有三頭主狼但多目標(biāo)下需要從檔案里選。MOGWO 的做法是把檔案按照目標(biāo)空間劃分成網(wǎng)格每個(gè)網(wǎng)格記錄解的密度然后從低密度區(qū)域優(yōu)先選擇領(lǐng)導(dǎo)者這樣可以讓搜索方向偏向尚未充分探索的區(qū)域。第三位置更新后的邊界處理。更新完新位置后決策變量要拉回可行域。很多實(shí)現(xiàn)直接寫死 min-max 裁剪但更好的做法是針對(duì)越界的變量重新在邊界鄰域內(nèi)隨機(jī)初始化保持種群多樣性。我看過不少直接用單目標(biāo) GWO 代碼改寫的“假 MOGWO”最常見的問題就是外部檔案只有插入沒有刪除或者領(lǐng)導(dǎo)者選擇不考慮密集度結(jié)果就是所有解堆到帕累托前沿的一小段上前沿兩端完全空掉。這種事不是跑不跑得通的問題而是結(jié)果有沒有參考價(jià)值的問題。2.3 MOGWO 的參數(shù)敏感度分析多目標(biāo)優(yōu)化里的參數(shù)設(shè)置直接影響收斂性和解分布質(zhì)量。MOGWO 主要參數(shù)包括參數(shù)推薦值影響種群規(guī)模30-50太小收斂不穩(wěn)太大計(jì)算量大迭代次數(shù)100-300看目標(biāo)函數(shù)復(fù)雜度帶閥點(diǎn)建議加大外部檔案大小50-100決定最終前沿點(diǎn)數(shù)網(wǎng)格劃分?jǐn)?shù)10-20影響領(lǐng)導(dǎo)者選擇多樣性我測(cè)過 IEEE30 系統(tǒng)的 EED 問題種群 40、迭代 200、檔案 80、網(wǎng)格 10 是比較穩(wěn)的組合。再往上加大迭代次數(shù)的收益率遞減因?yàn)槟繕?biāo)函數(shù)比較簡(jiǎn)單200 代之后帕累托前沿基本不變化。有一點(diǎn)要注意多目標(biāo)優(yōu)化的評(píng)價(jià)指標(biāo)不只是畫個(gè)圖IGD反轉(zhuǎn)世代距離、HV超體積、Spread分布度是三個(gè)最常用的量化指標(biāo)。一個(gè)嚴(yán)謹(jǐn)?shù)膶?duì)比實(shí)驗(yàn)至少算 HV 和 IGD跑個(gè) 30 次獨(dú)立實(shí)驗(yàn)取均值和方差才能說明算法穩(wěn)定性。很多新手只看一次跑的圖和最小成本就下結(jié)論這種方法論文里站不住腳。3. Matlab 實(shí)現(xiàn)環(huán)境經(jīng)濟(jì)調(diào)度的完整過程3.1 數(shù)據(jù)準(zhǔn)備IEEE30 節(jié)點(diǎn)參數(shù)與 B 矩陣在 Matlab 里做這個(gè)項(xiàng)目第一步是準(zhǔn)備數(shù)據(jù)文件。我習(xí)慣把機(jī)組參數(shù)、負(fù)荷、B 系數(shù)做成三個(gè)獨(dú)立的數(shù)據(jù)文件或者腳本避免全堆在一個(gè)主文件里數(shù)據(jù)項(xiàng)多了以后極其難排查。機(jī)組數(shù)據(jù)主要包括6 臺(tái)發(fā)電機(jī)的有功出力上下限、成本系數(shù)和排放系數(shù)。如果你用標(biāo)準(zhǔn) IEEE30 系統(tǒng)數(shù)據(jù)會(huì)看到類似這樣的一組典型值我隨手寫一組常用值做示例實(shí)際做實(shí)驗(yàn)請(qǐng)以你下載的 paper 數(shù)據(jù)為準(zhǔn)機(jī)組Pmin(MW)Pmax(MW)abcalphabetagammaG1502000.003752.0000.01234-0.10045.098G220800.017501.7500.00823-0.06052.407G315500.062501.0000.00418-0.05031.550G410350.008343.2500.00308-0.02400.806G510300.025003.0000.00534-0.05301.781G612400.025003.0000.00454-0.03901.348B 矩陣是 7x7 的對(duì)稱矩陣包含 B、B0、B00 三個(gè)分量。在 Matlab 里我建議存成結(jié)構(gòu)體sys.B [...]; % 6x6 二次項(xiàng)系數(shù) sys.B0 [...]; % 1x6 一次項(xiàng)系數(shù) sys.B00 ...; % 常數(shù)項(xiàng)網(wǎng)損計(jì)算公式function Ploss calc_loss(P, sys) n length(P); Ploss P. * sys.B * P sys.B0 * P. sys.B00; % 注意維度對(duì)應(yīng) end寫這一行網(wǎng)損計(jì)算時(shí)矩陣維度最容易搞反。P 是 6x1 列向量sys.B 是 6x6( P^T B P ) 得到標(biāo)量。sys.B0 我在存儲(chǔ)時(shí)寫成 6x1 列向量所以用 ( B0^T P ) 得到標(biāo)量。這種維度問題只有運(yùn)行起來報(bào)錯(cuò)或者算出來的值明顯不對(duì)才會(huì)發(fā)現(xiàn)。3.2 目標(biāo)函數(shù)與約束處理代碼怎么寫主目標(biāo)函數(shù)文件我建議輸入是 6 維出力向量輸出是成本和排放兩個(gè)標(biāo)量。代碼大致長這樣function [f_cost, f_emit] EED_objective(P, sys) % P: 決策變量向量 (nGen x 1) % f_cost: 總發(fā)電成本 % f_emit: 總污染物排放 P P(:); n length(P); % 成本函數(shù)二次型 f_cost sum(sys.a .* P.^2 sys.b .* P sys.c); % 排放函數(shù)二次型 f_emit sum(sys.alpha .* P.^2 sys.beta .* P sys.gamma); end約束處理是重頭戲。MOGWO 不像罰函數(shù)法那樣簡(jiǎn)單也不是所有越界都直接丟棄。我的做法是這樣出力上下限直接在更新位置后做 min-max 裁剪功率平衡約束處理完上下限之后檢查總出力、總負(fù)荷和網(wǎng)損之和的差但直接裁剪會(huì)產(chǎn)生一個(gè)問題6 個(gè)變量裁剪完總出力不等于負(fù)荷加網(wǎng)損等式約束被破壞。這里我用啟發(fā)式修正也叫功率平衡修復(fù)機(jī)制% 計(jì)算負(fù)荷和網(wǎng)損 total_demand sys.Pd; Ploss calc_loss(P, sys); error sum(P) - total_demand - Ploss; % 把誤差均攤到可調(diào)機(jī)組上 adjust error / n; for i 1:n P_new(i) P(i) - adjust; if P_new(i) sys.Pmin(i) P_new(i) sys.Pmin(i); elseif P_new(i) sys.Pmax(i) P_new(i) sys.Pmax(i); end end這種均攤修正不是完美的可能修完依然不滿足等式約束但誤差已經(jīng)很小。嚴(yán)格的做法是把等式約束轉(zhuǎn)成等式罰函數(shù)項(xiàng)加進(jìn)目標(biāo)函數(shù)里權(quán)衡一下我的經(jīng)驗(yàn)是對(duì)于 IEEE30 靜態(tài)調(diào)度啟發(fā)式修復(fù)加少量罰函數(shù)混合處理效果和精確求解差不多但速度快很多。3.3 MOGWO 主循環(huán)框架搭建我直接給一個(gè)精簡(jiǎn)版的 MOGWO 主循環(huán)框架你可以在此基礎(chǔ)上改寫成完整代碼% 初始化狼群 X init_wolf(N, sys); % N 只狼每個(gè)是6維向量 [f1, f2] evaluate_wolf(X, sys); repo init_repository(X, f1, f2); % 初始化外部檔案 for t 1:MaxIt a 2 - 2 * t / MaxIt; for i 1:N % 從檔案里選擇三個(gè)領(lǐng)導(dǎo)者基于網(wǎng)格密度 alpha select_leader(repo); beta select_leader(repo); delta select_leader(repo); % 按標(biāo)準(zhǔn)GWO公式更新位置對(duì)三個(gè)領(lǐng)導(dǎo)者分別計(jì)算再平均 X(i, :) update_position(X(i, :), alpha, beta, delta, a); % 邊界處理 X(i, :) clamp(X(i, :), sys.Pmin, sys.Pmax); % 啟發(fā)式功率平衡修復(fù) X(i, :) power_balance_fix(X(i, :), sys); end % 重新評(píng)估并更新檔案 [f1, f2] evaluate_wolf(X, sys); repo update_repository(repo, X, f1, f2); end這是結(jié)構(gòu)最穩(wěn)定的算法骨架。關(guān)鍵在于 update_repository 函數(shù)要寫好它總共做三件事合并新舊解、移除被支配解、如果檔案溢出則基于網(wǎng)格密度刪解。3.4 結(jié)果可視化與性能指標(biāo)計(jì)算跑完循環(huán)后你要輸出三樣?xùn)|西帕累托前沿圖橫軸成本、縱軸排放這通常是一張往下凸的散點(diǎn)圖。figure; plot(repo.cost, repo.emit, o, MarkerSize, 4); xlabel(Cost ($/h)); ylabel(Emission (ton/h)); grid on;最小成本解在檔案里找成本最小的方案輸出對(duì)應(yīng)出力向量。[min_cost, idx] min(repo.cost); best_X repo.X(idx, :); disp([Minimum Cost: , num2str(min_cost)]);我在實(shí)際跑的時(shí)候還會(huì)計(jì)算 HV 和 IGD 兩個(gè)指標(biāo)尤其是做算法對(duì)比時(shí)光靠圖無法說明問題。HV 計(jì)算需要一個(gè)參考點(diǎn)一般取所有解中每個(gè)目標(biāo)的最大值稍微放大一點(diǎn)然后算帕累托前沿和參考點(diǎn)之間覆蓋的超體積。寫指標(biāo)計(jì)算時(shí)可以考慮用第三方工具箱也可以手寫。如果只是為了項(xiàng)目驗(yàn)證我建議先手寫 IGD——邏輯上就是真實(shí)前沿上每個(gè)點(diǎn)到你算出的前沿的最近距離的平均值代碼不復(fù)雜對(duì)你有好處數(shù)據(jù)一致性也更可控。4. 結(jié)果分析與超參數(shù)調(diào)優(yōu)的實(shí)戰(zhàn)經(jīng)驗(yàn)4.1 帕累托前沿圖的解讀方式跑完 MOGWO 之后你看到的第一張圖是帕累托前沿。正常的 EED 前沿應(yīng)該是平滑下降的凸曲線左端點(diǎn)成本最低但排放最高右端點(diǎn)排放最低但成本最高兩端之間的點(diǎn)就是成本和排放的不同折中方案。如果你的圖出現(xiàn)以下情況大概率有問題前沿只覆蓋很小一段區(qū)間說明算法多樣性差領(lǐng)導(dǎo)者選擇策略失效前沿分布不連續(xù)、有空洞說明檔案更新策略不好某些網(wǎng)格區(qū)域沒有解被保留多數(shù)點(diǎn)擠在一起集中在某一側(cè)說明目標(biāo)歸一化沒做好或檔案溢出刪除策略偏頗我在調(diào)這種問題時(shí)通常優(yōu)先檢查外部檔案的網(wǎng)格劃分和 update_repository 里的種子選擇邏輯。有時(shí)候問題很簡(jiǎn)單是檔案大小設(shè)太小導(dǎo)致前沿點(diǎn)不夠密。4.2 單目標(biāo) vs 多目標(biāo)對(duì)比實(shí)驗(yàn)一個(gè)常見的加分做法是把 MOGWO 結(jié)果和單目標(biāo) GWO 結(jié)果對(duì)比。你可以用標(biāo)準(zhǔn) GWO 分別跑“最小成本”和“最小排放”兩個(gè)任務(wù)得到兩個(gè)極端解然后驗(yàn)證 MOGWO 前沿的兩個(gè)端點(diǎn)是否接近這兩個(gè)單目標(biāo)最優(yōu)解。這一步驗(yàn)證意義很大。如果 MOGWO 跑出的左端點(diǎn)成本比單目標(biāo) GWO 算出的最小成本還高 5%那么說明 MOGWO 的收斂性還有問題——極端解沒有被充分搜索到。我碰到過這種情況原因是種群規(guī)模太小且領(lǐng)導(dǎo)者選擇時(shí)過度偏向低密度區(qū)域?qū)е氯肿顑?yōu)方向被忽略。調(diào)整方向有兩個(gè)一是增大種群二是把 a 的下降速率調(diào)慢一點(diǎn)增強(qiáng)后期局部搜索。比如把 a 從線性衰減改成非線性衰減a 2 * (1 - (t / MaxIt)^0.7);這種非線性衰減在前期保持了更多探索后期收斂更精細(xì)。我實(shí)測(cè)在 IEEE30 EED 問題上改用平方根形式的衰減后最小成本端大約改善了 1.8%前沿分布也更均勻。4.3 多目標(biāo)評(píng)價(jià)指標(biāo)的使用心得做任何多目標(biāo)優(yōu)化項(xiàng)目建議至少要算三個(gè)指標(biāo)指標(biāo)全稱度量?jī)?nèi)容HVHypervolume解集覆蓋目標(biāo)空間的大小越大越好IGDInverted Generational Distance真實(shí)前沿到計(jì)算前沿的距離越小越好Spread分布度前沿點(diǎn)分布的均勻性越小越好需要注意的是HV 和 IGD 都需要一個(gè)“真實(shí)帕累托前沿”或者參考點(diǎn)。真實(shí)前沿怎么來常見方法是用一個(gè)超大種群、超多迭代次數(shù)的 MOGWO 跑出來的前沿作為近似真實(shí)前沿。記住這是“近似”但你用它做相對(duì)比較已經(jīng)足夠。我習(xí)慣每次實(shí)驗(yàn)獨(dú)立跑 30 次然后報(bào)告平均值±標(biāo)準(zhǔn)差。為什么是 30 次因?yàn)槿褐悄芩惴ㄓ须S機(jī)性單次結(jié)果沒有統(tǒng)計(jì)意義。在論文里寫“最好結(jié)果”是常見但不太嚴(yán)謹(jǐn)?shù)淖龇ㄎ医ㄗh至少同時(shí)報(bào)告均值和最差值這樣能體現(xiàn)算法的魯棒性。5. 實(shí)操中踩過的坑與排查技巧5.1 結(jié)果不收斂或前沿不完整的常見原因這一類問題我見得最多基本集中在四個(gè)原因第一個(gè)原因是初始種群生成不合法。如果初始狼群里很多個(gè)體的總出力與負(fù)荷和網(wǎng)損差的太遠(yuǎn)功率平衡修復(fù)函數(shù)會(huì)花很多代去修正消耗搜索資源。解決辦法是把初始種群生成時(shí)就添加一定比例的“合法個(gè)體”——也就是先產(chǎn)生隨機(jī)出力再通過啟發(fā)式修正讓部分個(gè)體滿足功率平衡讓種群起點(diǎn)就有一部分在可行域里。第二個(gè)原因是線性衰減的 a 參數(shù)不夠好。不同問題的收斂速度區(qū)別很大如果 a 降到零之后狼群還在大規(guī)模跳躍搜索就會(huì)在某代之后無法收斂。檢查辦法是畫出歷代的 HV 變化曲線如果后期 HV 還在劇烈波動(dòng)說明 a 衰減過快如果 HV 早早平穩(wěn)但最終值偏低說明衰減過慢。第三個(gè)原因是目標(biāo)函數(shù)數(shù)值跨度太大會(huì)影響網(wǎng)格劃分。比如成本是 600-800排放是 0.1-0.3如果不做歸一化就直接劃分網(wǎng)格網(wǎng)格在排放維度上會(huì)特別寬在成本維度上特別窄導(dǎo)致檔案選擇偏向某一個(gè)目標(biāo)。我會(huì)在計(jì)算網(wǎng)格密度之前對(duì)兩個(gè)目標(biāo)做 max-min 歸一化這是很容易被忽略但又非常關(guān)鍵的一步。第四個(gè)原因是外部檔案的更新時(shí)機(jī)不對(duì)。有些實(shí)現(xiàn)是等所有狼都更新完才刷新檔案有些是每只狼更新完就刷新。兩種方式各有優(yōu)劣但在 Matlab 的向量化實(shí)現(xiàn)里如果排序和更新邏輯寫錯(cuò)很容易造成檔案里出現(xiàn)跑出邊界外的臟數(shù)據(jù)。建議每次更新檔案后都打印一行檔案數(shù)量和邊界統(tǒng)計(jì)信息肉眼檢查是否出現(xiàn)丟失維度問題。5.2 Matlat 性能與代碼向量化MOGWO 在 IEEE30 上規(guī)模很小本來不用太關(guān)心性能。但如果你把種群調(diào)到 100、迭代 300 次純 for 循環(huán)跑到后面也會(huì)覺得卡。一個(gè)最簡(jiǎn)單有效的優(yōu)化是主循環(huán)盡量向量化評(píng)估函數(shù)一次性接收整個(gè)種群矩陣。function [F1, F2] evaluate_wolf_vectorized(X, sys) % X: N x nGen 矩陣每行是一個(gè)個(gè)體 F1 sum(sys.a .* X.^2 sys.b .* X sys.c, 2); F2 sum(sys.alpha .* X.^2 sys.beta .* X sys.gamma, 2); end相比于逐個(gè)個(gè)體調(diào)函數(shù)向量化至少能快 5 到 10 倍。如果你的問題擴(kuò)展到了動(dòng)態(tài)調(diào)度或者多時(shí)段代碼向量化的收益會(huì)更明顯。另一個(gè)實(shí)用技巧是提前為每個(gè)個(gè)體緩存目標(biāo)函數(shù)值。同一代內(nèi)同一個(gè)個(gè)體因?yàn)樾迯?fù)機(jī)制可能會(huì)被重復(fù)評(píng)估如果不緩存你就是在做無效計(jì)算。5.3 結(jié)果復(fù)現(xiàn)與實(shí)驗(yàn)記錄管理學(xué)術(shù)項(xiàng)目里最怕的就是跑完一組實(shí)驗(yàn)結(jié)果挺好的但過兩天想復(fù)現(xiàn)發(fā)現(xiàn)參數(shù)忘了。所以我建議從一開始就養(yǎng)成記錄實(shí)驗(yàn)配置的習(xí)慣種群數(shù)、迭代數(shù)、檔案大小、a 衰減方式、修復(fù)機(jī)制開關(guān)、隨機(jī)種子這些全部存到一個(gè)結(jié)構(gòu)體里輸出結(jié)果時(shí)同時(shí)輸出配置。exp_config struct(N, 40, MaxIt, 200, ... nGrid, 10, RepoSize, 80, ... alpha_decay, nonlinear_0.7, randseed, 1); rng(exp_config.randseed);保存結(jié)果時(shí)把 exp_config 一起存到 .mat 文件里這樣后續(xù)即使代碼改了你還能知道當(dāng)時(shí)的結(jié)果是怎么來的。這個(gè)習(xí)慣我吃過大虧以后養(yǎng)成的建議你直接學(xué)走。5.4 和 NSGA-II、MOPSO 對(duì)比時(shí)的公平性問題如果你打算在論文里和 NSGA-II、MOPSO 做對(duì)比一定要注意一個(gè)公平性原則所有算法使用相同的函數(shù)評(píng)估次數(shù)而不是相同的迭代次數(shù)。種群 40 迭代 200 的 MOGWO 總函數(shù)評(píng)估次數(shù)是 40×200 次那么 NSGA-II 也應(yīng)該用相同的總評(píng)估次數(shù)。這一點(diǎn)大家經(jīng)常忽略對(duì)比結(jié)果就失去了說服力。我曾經(jīng)見過一個(gè)實(shí)驗(yàn)MOGWO 迭代次數(shù)是 NSGA-II 的兩倍然后得出結(jié)論說 MOGWO 好這個(gè)結(jié)論就沒意義。公平對(duì)比還有個(gè)細(xì)節(jié)初始化方式要一致都用相同隨機(jī)種子和相同風(fēng)格的初始種群不能單獨(dú)為一個(gè)算法精心設(shè)計(jì)初始化。6. 擴(kuò)展方向與個(gè)人經(jīng)驗(yàn)分享6.1 加入閥點(diǎn)效應(yīng)與斜坡約束我前文提到的帶閥點(diǎn)成本函數(shù)建議在做完基本版后立刻補(bǔ)上。把成本函數(shù)改成[ F_i(P_i) a_i P_i^2 b_i P_i c_i |e_i \sin(f_i (P_i^{min} - P_i))| ]加了這個(gè)正弦項(xiàng)目標(biāo)函數(shù)變成非凸多峰MOGWO 的收斂性會(huì)明顯受到挑戰(zhàn)。這時(shí)候你再檢驗(yàn)算法效果才有真正的說服力。動(dòng)態(tài)調(diào)度再加爬坡約束和燃料約束問題就完全展開了。6.2 增加第三目標(biāo)除了成本和排放還可以把網(wǎng)損作為第三目標(biāo)最小化或者把系統(tǒng)電壓偏差加入目標(biāo)。多目標(biāo)數(shù)量增加后有兩個(gè)問題一是帕累托前沿變成平面或曲面可視化難度增加二是群體優(yōu)化算法的性能會(huì)下降因?yàn)榉侵浣鈹?shù)量爆炸檔案更新壓力變大。這時(shí)候通常需要引入基于分解的機(jī)制比如 MOEA/D 的思想或者改進(jìn)檔案維護(hù)策略。我的建議是第三目標(biāo)可以在同一套 MOGWO 框架下直接加但要注意網(wǎng)格劃分從二維變?nèi)S代碼的復(fù)雜度呈指數(shù)增長。初學(xué)階段先加“網(wǎng)損”而不是“電壓偏差”因?yàn)榫W(wǎng)損計(jì)算成本低電壓偏差需要潮流計(jì)算算力負(fù)擔(dān)完全不同。6.3 混合改進(jìn)方向MOGWO 的一個(gè)常見改進(jìn)思路是差分進(jìn)化混合在每次更新后用差分進(jìn)化算子的變異交叉代替標(biāo)準(zhǔn)位置更新增強(qiáng)開發(fā)能力。另一個(gè)思路是引入局部搜索對(duì)檔案中的非支配解做高斯擾動(dòng)和重新評(píng)估提高前沿的收斂精度。做混合改進(jìn)時(shí)要注意不要為了改進(jìn)而改進(jìn)。每一次改動(dòng)都要有明確的邏輯和實(shí)驗(yàn)驗(yàn)證。先跑標(biāo)準(zhǔn) MOGWO 算出基線結(jié)果然后記錄它的薄弱點(diǎn)比如哪個(gè)目標(biāo)收斂差、哪段前沿分布稀疏再有針對(duì)性地增強(qiáng)這種“問題驅(qū)動(dòng)”的改進(jìn)思路才是科研的正確路徑。從我自己的項(xiàng)目經(jīng)驗(yàn)來看做這種優(yōu)化調(diào)度類項(xiàng)目最忌一步到位。先把最簡(jiǎn)單的 IEEE30 靜態(tài)雙目標(biāo)跑通完整地理解目標(biāo)函數(shù)、約束和 MOGWO 的每一個(gè)機(jī)制再逐步擴(kuò)展。前面基礎(chǔ)打得扎實(shí)后面加約束、加目標(biāo)、換算法都會(huì)順手很多。最后提醒一句跑任何實(shí)驗(yàn)都記得保存隨機(jī)種子和配置參數(shù)別等出了好結(jié)果再去追根溯源那時(shí)候代價(jià)太大了。