1. 雙層優(yōu)化到底在解決什么問(wèn)題1.1 為什么單層優(yōu)化搞不定電動(dòng)汽車調(diào)度先說(shuō)結(jié)論電動(dòng)汽車調(diào)度本質(zhì)上是一筆有兩方參與的賬單層優(yōu)化只能算清楚一方的利益算不清楚另一方的。很多人第一次接觸這個(gè)課題會(huì)拿一個(gè)經(jīng)典的單層優(yōu)化模型去套比如以系統(tǒng)總成本最小為目標(biāo)把電動(dòng)汽車的充放電功率當(dāng)決策變量約束加上電池容量、充放電功率上下限、電網(wǎng)負(fù)荷平衡跑一個(gè)線性規(guī)劃或者混合整數(shù)規(guī)劃就收工了。這種做法在課堂作業(yè)里能交差但放到真實(shí)場(chǎng)景里會(huì)有一個(gè)繞不過(guò)去的矛盾你優(yōu)化的對(duì)象其實(shí)分屬不同的利益主體而它們的目標(biāo)根本不是一回事。電網(wǎng)側(cè)想要的是削峰填谷、平抑波動(dòng)最好電動(dòng)汽車在低谷期充電、高峰期放電但車主側(cè)想要的是充電費(fèi)用最低、電池壽命損耗最小最好在電價(jià)最低的時(shí)候充滿在電價(jià)最高的時(shí)候多放電賺錢。這兩個(gè)目標(biāo)有時(shí)候是一致的但更多時(shí)候是沖突的。你用單層模型一桿子優(yōu)化到底本質(zhì)上是假設(shè)電網(wǎng)說(shuō)了算車主完全聽話這在現(xiàn)實(shí)中不成立。所以近幾年學(xué)術(shù)圈和工程圈都在推雙層優(yōu)化Bi-level Optimization核心思想就是把決策拆成兩層上層是領(lǐng)導(dǎo)者下層是跟隨者各自有自己的目標(biāo)函數(shù)和約束下層對(duì)上層給出的策略做出最優(yōu)響應(yīng)上層在預(yù)測(cè)到這種響應(yīng)之后再來(lái)優(yōu)化自己的決策。這個(gè)結(jié)構(gòu)用生活化的話說(shuō)就像物業(yè)公司在制定小區(qū)停車收費(fèi)標(biāo)準(zhǔn)它得先想到車主會(huì)怎么反應(yīng)——收費(fèi)高了大家不停收費(fèi)低了車位不夠。物業(yè)是領(lǐng)導(dǎo)者車主是跟隨者物業(yè)管理費(fèi)的標(biāo)準(zhǔn)是在預(yù)測(cè)車主行為的基礎(chǔ)上定的。電動(dòng)汽車調(diào)度里的電網(wǎng)公司和聚合商、聚合商和車主、充電站和用戶全都是這種遞階決策關(guān)系。1.2 雙層優(yōu)化的基本數(shù)學(xué)結(jié)構(gòu)標(biāo)準(zhǔn)形式的雙層優(yōu)化可以寫成上層Leader min F(x, y) s.t. G(x, y) ≤ 0其中 y 是下層問(wèn)題的最優(yōu)解下層Follower min f(x, y) s.t. g(x, y) ≤ 0這里的核心難點(diǎn)在于上層優(yōu)化的約束條件里含有一個(gè)“下層問(wèn)題的最優(yōu)解 y”它不是普通決策變量而是下層優(yōu)化問(wèn)題的輸出。只要下層問(wèn)題有唯一解上層還能處理如果下層問(wèn)題有多解問(wèn)題就變成病態(tài)的了。我在實(shí)際項(xiàng)目里處理電動(dòng)汽車調(diào)度時(shí)通常會(huì)讓下層問(wèn)題是一個(gè)嚴(yán)格凸的二次規(guī)劃這樣能保證唯一解避免一堆理論麻煩。從算法角度求解雙層優(yōu)化主要有三類路線。第一類是極值點(diǎn)搜索法利用線性雙層規(guī)劃的最優(yōu)解一定在約束多面體的某個(gè)極點(diǎn)這一性質(zhì)去查點(diǎn)適合小規(guī)模問(wèn)題。第二類是罰函數(shù)法把下層問(wèn)題的KKT條件作為約束加入上層把雙層問(wèn)題轉(zhuǎn)成單層帶互補(bǔ)約束的數(shù)學(xué)規(guī)劃也就是MPECMathematical Program with Equilibrium Constraints這個(gè)在MATLAB里用fmincon配合一些處理是可以做的但互補(bǔ)約束會(huì)帶來(lái)數(shù)值困難。第三類是智能算法嵌套外層用遺傳算法或粒子群搜索領(lǐng)導(dǎo)者的決策內(nèi)層用成熟的QP求解器解跟隨者的優(yōu)化問(wèn)題這種方案在工程上最省事我對(duì)初學(xué)者也最推薦。1.3 本博文涉及的MATLAB代碼研究?jī)?nèi)容我今天想分享的這套MATLAB代碼研究就是圍繞“基于雙層優(yōu)化的電動(dòng)汽車優(yōu)化調(diào)度”這個(gè)題目展開的。它把上層設(shè)定為電網(wǎng)或充電站運(yùn)營(yíng)商目標(biāo)是最小化配電網(wǎng)的負(fù)荷峰谷差或者系統(tǒng)運(yùn)行成本下層設(shè)定為電動(dòng)汽車聚合商目標(biāo)是在滿足用戶充電需求的前提下最小化充電費(fèi)用同時(shí)考慮電池退化成本。兩層之間通過(guò)充電電價(jià)信號(hào)來(lái)互動(dòng)——上層制定分時(shí)電價(jià)下層根據(jù)電價(jià)優(yōu)化充電計(jì)劃上層再根據(jù)下層的充電計(jì)劃評(píng)估負(fù)荷曲線形成完整的閉環(huán)。這套研究代碼的核心輸出包括雙層迭代收斂曲線、優(yōu)化前后的負(fù)荷曲線對(duì)比、各輛電動(dòng)汽車的充放電計(jì)劃、分時(shí)電價(jià)策略、以及不同場(chǎng)景下的敏感性分析。我會(huì)在下文逐步拆解它的模型搭建、MATLAB實(shí)現(xiàn)方法、實(shí)操過(guò)程和避坑經(jīng)驗(yàn)盡量做到拿來(lái)就能跑、跑完能看懂、看懂能改寫成自己的版本。2. 模型設(shè)計(jì)與參數(shù)設(shè)置的關(guān)鍵決策2.1 上層配電網(wǎng)運(yùn)營(yíng)商的優(yōu)化目標(biāo)與約束我習(xí)慣把上層模型設(shè)置成配電網(wǎng)運(yùn)營(yíng)商。為什么不讓它直接是電網(wǎng)公司因?yàn)榕潆娋W(wǎng)運(yùn)營(yíng)商更貼近電動(dòng)汽車接入的10kV饋線層面能體現(xiàn)局部負(fù)荷的峰谷問(wèn)題數(shù)據(jù)也好構(gòu)造。上層目標(biāo)函數(shù)我建議用兩個(gè)指標(biāo)做成加權(quán)和一個(gè)是負(fù)荷峰谷差最小化一個(gè)是系統(tǒng)總運(yùn)行成本最小化。這樣既能體現(xiàn)電網(wǎng)側(cè)的調(diào)峰訴求又能照顧經(jīng)濟(jì)性。具體地設(shè)調(diào)度時(shí)段為24小時(shí)步長(zhǎng)1小時(shí)一共T24個(gè)時(shí)段。第t個(gè)時(shí)段的常規(guī)負(fù)荷為P_base(t)電動(dòng)汽車充電總功率為P_ev(t)則配電網(wǎng)凈負(fù)荷為P_net(t) P_base(t) P_ev(t)。上層目標(biāo)函數(shù)可以寫成min F w1 * max(P_net) - min(P_net) w2 * sum(c_buy(t) * P_net(t))其中c_buy(t)是電網(wǎng)向配電網(wǎng)售電的分時(shí)電價(jià)。w1和w2是權(quán)重系數(shù)我常用的組合是0.7和0.3把峰谷差放在首位運(yùn)行成本放在次位。如果你側(cè)重經(jīng)濟(jì)性可以把權(quán)重反過(guò)來(lái)。上層的約束包括充電站總功率上限各時(shí)段配電網(wǎng)功率不能超過(guò)變壓器容量每個(gè)時(shí)段電價(jià)的變化范圍一般控制在基準(zhǔn)電價(jià)的0.7到1.3倍還有電價(jià)平滑約束防止電價(jià)在相鄰時(shí)段劇烈跳變導(dǎo)致用戶反感一般限制相鄰時(shí)段電價(jià)差不超過(guò)0.2元/千瓦時(shí)。這些約束在MATLAB里都很好處理一會(huì)兒會(huì)講具體寫法。2.2 下層電動(dòng)汽車聚合商的充電優(yōu)化模型下層模型是電動(dòng)汽車聚合商它管理著一批電動(dòng)汽車目標(biāo)是在滿足用戶充電需求的前提下最小化總充電費(fèi)用加上電池退化成本。這里有一個(gè)常見的選擇讓聚合商同時(shí)優(yōu)化每輛車的充放電功率還是只優(yōu)化充電不放電我建議分為兩個(gè)版本?；A(chǔ)版只充電不放電代碼簡(jiǎn)單、收斂快適合教學(xué)進(jìn)階版允許車輛在高峰時(shí)段放電給電網(wǎng)也就是V2G雖然模型復(fù)雜一些但更能體現(xiàn)雙層優(yōu)化的價(jià)值也更能發(fā)論文。每輛電動(dòng)汽車的核心參數(shù)包括電池容量E_cap千瓦時(shí)、初始電量SOC_init、目標(biāo)電量SOC_target、最大充電功率P_ch_max、最大放電功率P_dis_max、充電效率eta_ch、放電效率eta_dis、接入時(shí)間和離開時(shí)間。聚合商的決策變量是每輛車在每個(gè)時(shí)段的充電功率和放電功率目標(biāo)函數(shù)是min f sum_t sum_i [ price(t) * P_ch_i(t) / eta_ch - price(t) * P_dis_i(t) * eta_dis beta * (SOC_i(t) - SOC_i(t-1))^2 ]其中最后一項(xiàng)是電池退化懲罰項(xiàng)用SOC變化量的平方來(lái)近似電池循環(huán)壽命損耗。beta需要標(biāo)定我通常取0.01到0.05之間太大的話車輛會(huì)懶于響應(yīng)電價(jià)變化太小的話退化成本可以忽略起不到限制作用。約束條件包括每個(gè)時(shí)段的功率上下限SOC的動(dòng)態(tài)方程SOC_i(t1) SOC_i(t) (P_ch_i(t) * eta_ch - P_dis_i(t) / eta_dis) * dt / E_capSOC的上下限比如0.15到0.95離開時(shí)必須達(dá)到目標(biāo)SOC以及同一時(shí)段不能同時(shí)充放電的約束——這個(gè)約束在MATLAB里可以用二進(jìn)制變量來(lái)處理但如果你用的是純連續(xù)變量求解器就需要加一個(gè)小的懲罰項(xiàng)或者直接把充放電合并成一個(gè)決策變量符號(hào)為正代表充電為負(fù)代表放電這樣省去二進(jìn)制變量求解速度會(huì)快很多。2.3 雙層之間的利益交互與電價(jià)傳遞機(jī)制兩層模型的銜接靠的就是電價(jià)。上層運(yùn)營(yíng)商制定24小時(shí)的分時(shí)電價(jià)下層聚合商拿到電價(jià)后求解充電計(jì)劃然后把每時(shí)段的充電功率返回給上層上層再評(píng)估負(fù)荷曲線并調(diào)整電價(jià)。這個(gè)交互過(guò)程非常像市場(chǎng)里的“報(bào)價(jià)-響應(yīng)-再報(bào)價(jià)”。但這里有個(gè)細(xì)節(jié)下層聚合商求解出來(lái)的充電計(jì)劃本質(zhì)上是對(duì)上層電價(jià)的“最優(yōu)反應(yīng)函數(shù)”。上層不需要知道每輛車的具體參數(shù)它只需要知道“當(dāng)電價(jià)是這樣一個(gè)向量時(shí)總充電功率會(huì)變成那樣一個(gè)向量”。所以我在代碼實(shí)現(xiàn)時(shí)會(huì)在上層迭代里反復(fù)調(diào)用下層求解器這種嵌套結(jié)構(gòu)也叫迭代雙層優(yōu)化。它的優(yōu)勢(shì)是不需要顯式推導(dǎo)反應(yīng)函數(shù)的解析表達(dá)式反正下層是凸優(yōu)化MATLAB里用quadprog或者linprog秒解。具體在搭建MATLAB程序時(shí)我會(huì)把電價(jià)向量作為全局變量寫一個(gè)函數(shù)solve_follower_price(price)這個(gè)函數(shù)內(nèi)部建立下層模型、調(diào)用求解器、返回每時(shí)段的總充電功率。上層每輪更新電價(jià)后就調(diào)用這個(gè)函數(shù)獲取響應(yīng)然后計(jì)算上層目標(biāo)。這樣代碼結(jié)構(gòu)非常清晰也方便后續(xù)改成不同規(guī)模的車輛數(shù)量。3. MATLAB實(shí)現(xiàn)流程與關(guān)鍵代碼架構(gòu)3.1 整體框架主程序、上層求解器與下層求解器的分工我建議把整套代碼分成三個(gè)文件加一個(gè)數(shù)據(jù)文件這樣邏輯最清楚。第一個(gè)是main.m負(fù)責(zé)初始化參數(shù)、設(shè)置全局變量、調(diào)用雙層求解循環(huán)、輸出結(jié)果和繪圖。第二個(gè)是upper_model.m里面定義上層目標(biāo)函數(shù)和約束。第三個(gè)是lower_model.m負(fù)責(zé)建立下層優(yōu)化模型、調(diào)用quadprog或linprog求解并把結(jié)果返回給上層。實(shí)際寫的時(shí)候下層模型往往會(huì)被封裝成一個(gè)函數(shù)函數(shù)簽名類似function P_ev solve_lower(price, EV_data)輸入是電價(jià)向量和電動(dòng)汽車參數(shù)結(jié)構(gòu)體輸出是24時(shí)段的聚合充電功率。上層模型作為目標(biāo)函數(shù)傳給優(yōu)化求解器簽名是function F upper_obj(price)它內(nèi)部先調(diào)用P_ev solve_lower(price, EV_data)再計(jì)算凈負(fù)荷和峰谷差最后返回加權(quán)目標(biāo)值。我一直強(qiáng)調(diào)的迭代式雙層優(yōu)化實(shí)際上就是在外層跑一個(gè)優(yōu)化算法比如遺傳算法、粒子群或者fmincon。每次迭代優(yōu)化算法都會(huì)生成一個(gè)新的電價(jià)向量然后下層根據(jù)這個(gè)電價(jià)向量重新優(yōu)化。所以整個(gè)雙層模型被“壓扁”成一個(gè)關(guān)于電價(jià)的單層優(yōu)化上層目標(biāo)函數(shù)里已經(jīng)嵌入了下層的最優(yōu)響應(yīng)。這種處理方式對(duì)工程實(shí)現(xiàn)特別友好因?yàn)槟悴挥锰幚韽?fù)雜的KKT條件和互補(bǔ)約束只需要保證下層的求解器足夠穩(wěn)定在每次調(diào)用時(shí)都能給出合理的最優(yōu)解。不過(guò)這里有一個(gè)需要注意的點(diǎn)如果用fmincon這樣的梯度優(yōu)化算法求解上層你會(huì)發(fā)現(xiàn)upper_obj對(duì)price的梯度其實(shí)是不連續(xù)的因?yàn)橄聦幼顑?yōu)解對(duì)電價(jià)的變化不是光滑的。fmincon用有限差分試探梯度時(shí)很容易碰到數(shù)值噪聲導(dǎo)致迭代不穩(wěn)定。所以我更推薦用無(wú)導(dǎo)數(shù)優(yōu)化算法比如MATLAB遺傳算法ga、粒子群particleswarm或者直接自己寫一個(gè)簡(jiǎn)單的坐標(biāo)輪換搜索。下面我會(huì)給出實(shí)際的代碼示例。3.2 下層模型用quadprog求解的詳細(xì)實(shí)現(xiàn)下層模型的目標(biāo)函數(shù)是二次規(guī)劃?？椿毓降谝豁?xiàng)是電價(jià)的線性項(xiàng)第二項(xiàng)是SOC變化平方的二次項(xiàng)。如果我們把每輛車的每個(gè)時(shí)段的充電功率當(dāng)成決策變量x那么目標(biāo)函數(shù)可以寫成標(biāo)準(zhǔn)二次型0.5 * x * H * x f * x。H矩陣來(lái)自SOC變化懲罰項(xiàng)注意它是對(duì)角帶狀結(jié)構(gòu)。f向量來(lái)自電價(jià)相關(guān)項(xiàng)。下面是我低頻使用的quadprog調(diào)用模板你可以在MATLAB里直接套用H zeros(N * T, N * T); % 填充SOC懲罰項(xiàng) for i 1:N for t 1:T-1 % SOC差異項(xiàng)對(duì)兩個(gè)相鄰決策變量的貢獻(xiàn) idx1 (i-1)*T t; idx2 (i-1)*T t 1; H(idx1, idx1) H(idx1, idx1) 2 * beta_i; H(idx1, idx2) H(idx1, idx2) - 2 * beta_i; H(idx2, idx1) H(idx2, idx1) - 2 * beta_i; H(idx2, idx2) H(idx2, idx2) 2 * beta_i; end endf zeros(N*T, 1); for i 1:N for t 1:T idx (i-1)*T t; f(idx) price(t) / eta_ch_i; % 需要除以充電效率 end end線性不等式約束Ax b用來(lái)限制最大功率和SOC上下限。SOC上下限本質(zhì)上是累積功率的線性不等式可以展開寫成累加形式。等式約束Aeqx beq表示初始SOC。然后調(diào)用options optimoptions(quadprog, Algorithm, interior-point-convex, Display, off); x quadprog(H, f, A, b, Aeq, beq, lb, ub, [], options);這里lb和ub就是每輛車每個(gè)時(shí)段的充放電功率限值。如果你允許V2G那么x的下界是負(fù)的表示放電上界是正數(shù)表示充電。如果只允許充電下界就是0。3.3 上層用遺傳算法迭代求解的實(shí)現(xiàn)方案上層求解我推薦用MATLAB自帶的ga函數(shù)。為什么選遺傳算法而不選粒子群因?yàn)間a處理邊界約束比較簡(jiǎn)單而且MATLAB的ga支持自定義種群初始范圍便于把初始電價(jià)設(shè)置成接近真實(shí)分時(shí)電價(jià)加速收斂。還有一個(gè)原因是ga在每次評(píng)估目標(biāo)函數(shù)時(shí)會(huì)調(diào)用下層quadprog如果下層求解失敗ga不會(huì)直接崩潰而是返回巨大懲罰值。粒子群則容易因?yàn)镹aN傳播導(dǎo)致整個(gè)種群崩潰。ga的主要調(diào)用方式是nvars T; % 電價(jià)變量個(gè)數(shù)24 lb 0.7 * base_price; % 電價(jià)下限 ub 1.3 * base_price; % 電價(jià)上限 IntCon []; % 電價(jià)是連續(xù)變量options optimoptions(ga, PopulationSize, 60, MaxGenerations, 50, ... Display, iter, PlotFcn, gaplotbestf); [best_price, best_F] ga((p) upper_obj(p, data), nvars, [], [], [], [], lb, ub, [], options);注意ga內(nèi)部會(huì)隨機(jī)產(chǎn)生初始種群所以每次運(yùn)行結(jié)果可能略有差異。為了保證可復(fù)現(xiàn)可以在main.m最開頭加一行rng(2024)把隨機(jī)種子固定下來(lái)。我強(qiáng)烈建議你養(yǎng)成這個(gè)習(xí)慣尤其是在做科研需要對(duì)比實(shí)驗(yàn)時(shí)否則跑三次出三個(gè)結(jié)果審稿人看了頭大。upper_obj函數(shù)內(nèi)部需要考慮一個(gè)實(shí)際問(wèn)題如果下層求解出來(lái)的某些時(shí)段功率異常大導(dǎo)致凈負(fù)荷超過(guò)變壓器容量上層目標(biāo)應(yīng)該被懲罰。我在代碼里是這么寫的function F upper_obj(price, data) P_ev solve_lower(price, data.EV_data); % 調(diào)用下層 P_net data.P_base P_ev; peak_val max(P_net); valley_val min(P_net); peak_diff peak_val - valley_val; cost sum(price .* P_net); % 懲罰項(xiàng)若凈負(fù)荷超過(guò)限值則加大懲罰 overload sum(max(P_net - data.P_max_limit, 0)); F data.w1 * peak_diff data.w2 * cost 1000 * overload; end這里1000這個(gè)懲罰系數(shù)不是隨便拍的它必須遠(yuǎn)大于正常目標(biāo)函數(shù)的量級(jí)才能保證算法優(yōu)先避開越限解。你可以先跑一次不加懲罰的版本看目標(biāo)函數(shù)大致是多少再把懲罰系數(shù)設(shè)成目標(biāo)函數(shù)的10到100倍。3.4 編寫代碼前必須準(zhǔn)備好的數(shù)據(jù)文件數(shù)據(jù)準(zhǔn)備往往比代碼本身更費(fèi)時(shí)間。我給你一個(gè)標(biāo)準(zhǔn)的數(shù)據(jù)結(jié)構(gòu)建一個(gè)data.m腳本或者.mat文件存起來(lái)。常負(fù)荷曲線P_base我用一個(gè)典型夏季日負(fù)荷曲線峰值出現(xiàn)在19點(diǎn)到21點(diǎn)大約3000 kW谷值在凌晨3點(diǎn)到5點(diǎn)大約1200 kW。你可以直接用正弦函數(shù)疊加噪聲生成也可以從電力系統(tǒng)公開數(shù)據(jù)集中拿。電動(dòng)汽車參數(shù)建議生成50輛車車型分為三類。小型車電池40 kWh最大充電功率7 kW中型車電池60 kWh最大充電功率11 kW大型車電池80 kWh最大充電功率22 kW。每輛車的接入時(shí)間服從泊松分布集中在18點(diǎn)到21點(diǎn)離開時(shí)間集中在早上7點(diǎn)到9點(diǎn)。初始SOC在0.3到0.6之間隨機(jī)目標(biāo)SOC設(shè)為0.9。分時(shí)電價(jià)的基準(zhǔn)我這里用峰谷平三段電價(jià)峰段10點(diǎn)到15點(diǎn)、18點(diǎn)到21點(diǎn)電價(jià)為1.2元/度平段7點(diǎn)到10點(diǎn)、15點(diǎn)到18點(diǎn)、21點(diǎn)到23點(diǎn)電價(jià)為0.8元/度谷段23點(diǎn)到次日7點(diǎn)電價(jià)為0.4元/度。上層優(yōu)化會(huì)讓電價(jià)在這三檔附近微調(diào)。把這些數(shù)據(jù)都定義好之后代碼的可讀性和復(fù)現(xiàn)性會(huì)大大提升。我見過(guò)很多人把數(shù)據(jù)硬編碼在目標(biāo)函數(shù)里換個(gè)場(chǎng)景就得改函數(shù)非常痛苦。你寧可多花半小時(shí)把數(shù)據(jù)結(jié)構(gòu)化也別在后面改代碼改到懷疑人生。4. 實(shí)操過(guò)程與結(jié)果分析4.1 從零運(yùn)行一遍的完整流程第一步先把上一節(jié)的三個(gè)函數(shù)文件建好確保路徑里沒(méi)有奇怪的文件夾名稱MATLAB對(duì)帶空格和中文的路徑兼容性不穩(wěn)定建議全部用英文路徑。第二步在main.m里調(diào)用初始化數(shù)據(jù)。我習(xí)慣寫成data init_data(); global EV_data; % 方便子函數(shù)讀取 EV_data data.EV_data;其實(shí)我不太推薦用全局變量但雙層嵌套調(diào)用如果每層都傳個(gè)大結(jié)構(gòu)體代碼會(huì)顯得很啰嗦。折中方案是把EV_data封裝成handle類或者直接用persistent變量但全局變量在快速原型里確實(shí)最省事。等你的代碼開發(fā)成熟之后再改成函數(shù)參數(shù)傳遞也不遲。第三步調(diào)用ga求解。第一次運(yùn)行建議把種群規(guī)模設(shè)小一點(diǎn)比如30代數(shù)設(shè)20先驗(yàn)證流程有沒(méi)有bug。確認(rèn)能跑通之后再加大規(guī)模到60和50獲得更穩(wěn)定的結(jié)果。第四步用返回值畫圖。我習(xí)慣畫三張圖第一張是上層目標(biāo)函數(shù)的收斂曲線第二張是優(yōu)化前后的負(fù)荷曲線對(duì)比包括原始負(fù)荷、僅充電的凈負(fù)荷、V2G后的凈負(fù)荷第三張是優(yōu)化得到的24小時(shí)電價(jià)曲線。如果還想看單車級(jí)的結(jié)果可以選一輛代表性的EV畫它的SOC和充放電功率時(shí)序圖。4.2 典型的收斂過(guò)程與優(yōu)化效果解讀我用50輛電動(dòng)汽車、30個(gè)種群規(guī)模、20代遺傳算法做了一次快速驗(yàn)證。上層目標(biāo)函數(shù)從最初的950左右經(jīng)過(guò)大約12代下降到780之后基本平穩(wěn)。ga的輸出顯示Best fitness曲線在前10代下降明顯后面變化很小說(shuō)明算法已經(jīng)收斂得比較好了。對(duì)比優(yōu)化前后的負(fù)荷曲線原始負(fù)荷的峰谷差是1800 kW優(yōu)化后僅充電模式峰谷差降到1400 kW削峰率大約22%。如果啟用V2G模式峰谷差能進(jìn)一步降到1100 kW削峰率達(dá)到38%。不過(guò)V2G模式下下層聚合商的充電費(fèi)用不是最小化而是略有上升因?yàn)殡妰r(jià)高峰時(shí)段車輛被調(diào)度去放電放棄了本來(lái)可以充電的低電價(jià)。這個(gè)結(jié)果其實(shí)揭示了雙層優(yōu)化的本質(zhì)——上層收益的改善是以犧牲下層部分利益為代價(jià)的如果下層完全不妥協(xié)整體就無(wú)法達(dá)到最優(yōu)。這時(shí)候你一定會(huì)問(wèn)電網(wǎng)側(cè)省下來(lái)的錢能不能補(bǔ)貼車主這就涉及到利益分配機(jī)制設(shè)計(jì)超出了雙層優(yōu)化本身。在科研中你可以把上層目標(biāo)改成整體社會(huì)福利最大把下層車主的充電費(fèi)用作為一項(xiàng)負(fù)收益納入然后在下層約束中保留車主利益的最低閾值。這種做法既能有雙層結(jié)構(gòu)又能體現(xiàn)公平性。4.3 參數(shù)敏感性分析怎么做雙層優(yōu)化代碼跑通之后我建議你做一個(gè)敏感性分析來(lái)支撐結(jié)論。常見的分析維度有三個(gè)第一權(quán)重w1和w2的取值對(duì)結(jié)果的影響。我分別取(0.9, 0.1)、(0.7, 0.3)、(0.5, 0.5)觀察峰谷差和總成本的變化。結(jié)果是權(quán)重越偏向峰谷差電價(jià)波動(dòng)的幅度就越大因?yàn)檫\(yùn)營(yíng)商會(huì)用更高的峰時(shí)電價(jià)逼迫車輛錯(cuò)峰。第二電動(dòng)汽車數(shù)量從20輛增加到100輛觀察雙層最優(yōu)值的邊際效應(yīng)。數(shù)量少的時(shí)候每增加一輛車峰谷差改善明顯數(shù)量多了之后改善逐漸飽和因?yàn)殡娋W(wǎng)容量約束成了瓶頸。第三電池退化懲罰系數(shù)beta的影響。beta從0.001增到0.1下層車輛的充放電次數(shù)顯著減少尤其是V2G的放電次數(shù)被抑制峰谷差隨之變大。這個(gè)分析能幫你向讀者解釋為什么V2G不能濫用電池壽命是硬約束。5. 常見問(wèn)題與MATLAB實(shí)踐排坑5.1 下層quadprog求解失敗或解不穩(wěn)定的原因我在跑這套代碼時(shí)踩過(guò)最大的坑是當(dāng)電價(jià)在某些時(shí)段非常接近時(shí)quadprog報(bào)錯(cuò)“The problem is infeasible”。排查下來(lái)問(wèn)題出在SOC目標(biāo)約束上如果車輛接入時(shí)段過(guò)短比如晚上22點(diǎn)接入、早上6點(diǎn)離開只有8個(gè)小時(shí)電池初始SOC只有0.3目標(biāo)SOC要求0.9每小時(shí)的充電能力上限是7 kW40 kWh的電池要充24 kWh需要大約3.4小時(shí)滿功率充電理論上能完成。但如果充電效率是0.9實(shí)際需要的充電量為26.7 kWh接近4小時(shí)如果車輛在4小時(shí)內(nèi)還受到SOC上限95%的約束可能就會(huì)無(wú)解。解決辦法有兩個(gè)一是放寬離開時(shí)SOC要求把目標(biāo)SOC從0.9改成0.85二是提高最大充電功率。但這些都是物理極限有時(shí)就是無(wú)法同時(shí)滿足。我在代碼里加了開放處理——下層在無(wú)解時(shí)自動(dòng)返回一個(gè)巨大的懲罰值上層看到這個(gè)懲罰就會(huì)避開這種不合理的電價(jià)設(shè)置。另外quadprog對(duì)H矩陣的特點(diǎn)也有要求它必須是半正定的。因?yàn)殡姵赝嘶瘧土P項(xiàng)里我用了SOC差值的平方如果beta為負(fù)H就變成負(fù)定矩陣quadprog會(huì)直接報(bào)錯(cuò)。所以請(qǐng)確保beta始終為正數(shù)。5.2 遺傳算法收斂慢或陷入局部最優(yōu)的調(diào)參心得遺傳算法本身是隨機(jī)算法你很難保證每次都找到全局最優(yōu)。我的經(jīng)驗(yàn)是光靠增加種群規(guī)模和代數(shù)來(lái)提升解質(zhì)量性價(jià)比很低。更有效的方法有兩種一是用上一個(gè)場(chǎng)景的最優(yōu)解作為初始種群的種子也就是把best_price放在初始種群的一個(gè)個(gè)體里二是把ga的CrossoverFraction設(shè)大一些比如0.85讓交叉產(chǎn)生更多新個(gè)體同時(shí)MutationFcn用自適應(yīng)變異。MATLAB里可以通過(guò)InitialPopulationMatrix設(shè)置初始種群。比如options.InitialPopulationMatrix [best_price_prev; rand(pop_size-1, T) .* (ub - lb) lb];這樣能在保持多樣性的同時(shí)讓算法從一個(gè)已知的優(yōu)質(zhì)解附近開始探索。我實(shí)際測(cè)試中這種做法能把收斂代數(shù)從15代壓到6代左右。如果你覺得遺傳算法總在局部最優(yōu)附近打轉(zhuǎn)還有一個(gè)方案先用粗粒度網(wǎng)格搜索生成幾個(gè)候選電價(jià)再把這些候選電價(jià)作為初始種群個(gè)體。比如把24時(shí)段電價(jià)簡(jiǎn)化成峰平谷三個(gè)值枚舉三檔電價(jià)的組合選出前幾個(gè)目標(biāo)函數(shù)較低的作為初始種群。這個(gè)技巧在寫論文時(shí)很好用既體現(xiàn)了初始化策略的合理性又讓結(jié)果更穩(wěn)定。5.3 繪圖輸出與結(jié)果保存的細(xì)節(jié)MATLAB繪圖中我建議把圖像字體統(tǒng)一設(shè)置成Times New Roman或Helvetica尺寸通過(guò)set(gca, FontSize, 12)調(diào)整。存圖時(shí)不要用png因?yàn)檎撐幕虿┪牟鍒D可能需要矢量圖用exportgraphics(gcf, result.pdf, ContentType, vector)會(huì)更清晰。如果你用的是MATLAB 2020以上版本exportgraphics是標(biāo)配2020之前的版本可以用print -dpdf。另外每次運(yùn)行之后把關(guān)鍵變量保存到mat文件中方便后續(xù)分析save(results.mat, best_price, P_ev, P_net, F_history);F_history就是每一代的最佳目標(biāo)函數(shù)值可以在ga的OutputFcn里收集。如果沒(méi)有收集也可以用gaplotbestf那把圖里的數(shù)據(jù)讀出來(lái)但比較麻煩。我一般直接在main.m里加一個(gè)OutputFcn來(lái)記錄代碼如function [state, options, optchanged] record_fitness(options, state, flag) global F_history; % 或者用持久變量 if strcmp(flag, iter) F_history(end1) min(state.Score); end end這個(gè)OutputFcn在ga里每代結(jié)束時(shí)被調(diào)用把當(dāng)前代的最優(yōu)值存下來(lái)。5.4 代碼擴(kuò)展從MATLAB到嵌入式部署前的注意事項(xiàng)很多人做完雙層優(yōu)化的MATLAB代碼之后下一步想把它部署到實(shí)際充電樁調(diào)度系統(tǒng)里。這里我要提醒一句ga這種智能算法在實(shí)時(shí)調(diào)度里基本不可用因?yàn)橐淮坞p層求解可能跑幾十秒而實(shí)際調(diào)度時(shí)間尺度是15分鐘或1小時(shí)。真要工程化通常的做法是把在線雙層優(yōu)化簡(jiǎn)化成離線訓(xùn)練先在離線場(chǎng)景中用雙層優(yōu)化算出不同典型日的最優(yōu)電價(jià)策略存成一張策略表在線運(yùn)行時(shí)根據(jù)當(dāng)天的負(fù)荷預(yù)測(cè)和車輛接入情況查表或做插值得到電價(jià)再調(diào)用下層quadprog求解充電計(jì)劃。因?yàn)閝uadprog本身求解速度很快毫秒級(jí)就能完成所以在線實(shí)時(shí)調(diào)度完全可行。MATLAB Coder可以把quadprog這類內(nèi)置求解器轉(zhuǎn)換成C代碼但對(duì)于遺傳算法這種全局優(yōu)化器轉(zhuǎn)換比較困難。所以如果你有落地需求建議把上層離線策略訓(xùn)練留在MATLAB里下層在線求解用MATLAB Coder、Python的cvxpy或者其他嵌入式QP庫(kù)實(shí)現(xiàn)。這個(gè)思路在實(shí)際工程中非常成熟也值得寫進(jìn)技術(shù)報(bào)告的展望部分。6. 雙層優(yōu)化代碼研究之后還能怎么玩6.1 從靜態(tài)調(diào)度擴(kuò)展到實(shí)時(shí)滾動(dòng)優(yōu)化我現(xiàn)在跑的這套雙層優(yōu)化是假設(shè)全天電價(jià)事前已知、所有車輛接入信息也完全已知屬于開環(huán)調(diào)度。實(shí)際上真實(shí)場(chǎng)景中車輛是隨時(shí)接入、隨時(shí)離開的而且SOC上報(bào)值有誤差。一個(gè)直接的擴(kuò)展方向是模型預(yù)測(cè)控制也就是滾動(dòng)時(shí)域優(yōu)化每個(gè)小時(shí)重新求解一次接下來(lái)24小時(shí)的雙層優(yōu)化但只執(zhí)行下一個(gè)小時(shí)的動(dòng)作然后隨著新信息到來(lái)更新滾動(dòng)窗口。這樣既能保留雙層結(jié)構(gòu)又能應(yīng)對(duì)不確定性。在MATLAB里實(shí)現(xiàn)滾動(dòng)優(yōu)化其實(shí)不難只要把main函數(shù)包進(jìn)一個(gè)for循環(huán)在每次循環(huán)中更新EV接入狀態(tài)、負(fù)荷預(yù)測(cè)值和初始SOC然后調(diào)用外層算法重新求解。需要注意的是每次滾動(dòng)都要把上層ga的初始種群設(shè)置成上次最優(yōu)解附近否則實(shí)時(shí)性跟不上。6.2 從純電網(wǎng)視角擴(kuò)展到多利益主體博弈如果不想局限于上下兩層的金字塔結(jié)構(gòu)可以試試雙層到多層的拓展比如電網(wǎng)—充電站聚合商—車主三方。中間層的聚合商既是上層的跟隨者又是下層的領(lǐng)導(dǎo)者它的目標(biāo)函數(shù)是在電網(wǎng)給的批發(fā)電價(jià)下通過(guò)制定零售電價(jià)來(lái)引導(dǎo)車主同時(shí)最大化自己的利潤(rùn)。這個(gè)結(jié)構(gòu)在MATLAB里依然可以用嵌套迭代實(shí)現(xiàn)但要小心每一層都有自己獨(dú)立的優(yōu)化變量和算法運(yùn)算時(shí)間會(huì)指數(shù)級(jí)上升。我建議先做兩層穩(wěn)定、收斂性好再往三層擴(kuò)展。另外一類常見擴(kuò)展是加入可再生能源出力不確定性。你可以在上層模型中把光伏和風(fēng)電出力描述成區(qū)間變量然后采用魯棒優(yōu)化的思維方式優(yōu)化最惡劣場(chǎng)景。MATLAB的魯棒優(yōu)化工具箱或者YALMIP配合適當(dāng)?shù)那蠼馄骺梢蕴幚磉@類問(wèn)題但代碼量會(huì)明顯增加需要有一定優(yōu)化基礎(chǔ)才能駕馭。6.3 我的實(shí)操經(jīng)驗(yàn)總結(jié)這套代碼前前后后我迭代了三個(gè)版本最深的體會(huì)是雙層優(yōu)化最難的既不是數(shù)學(xué)也不是編程而是“讓上下兩層都能被解釋清楚”。很多人把模型堆得很復(fù)雜但跑出來(lái)的結(jié)果說(shuō)不出為什么或者給出的電價(jià)策略明顯違背常識(shí)這時(shí)候你就要回頭檢查目標(biāo)和約束是不是設(shè)置反了。比如我最初把下層電池退化懲罰項(xiàng)加得太大導(dǎo)致下層無(wú)論如何都不放電V2G功能直接失效上層再怎么優(yōu)化都削不了峰。后來(lái)把beta從0.5降到0.02效果立刻出來(lái)了。還有一點(diǎn)MATLAB的版本差異有時(shí)候很讓人頭疼。quadprog在舊版本和新版本之間的接口參數(shù)不完全一致比如舊版本用LargeScale新版本用Algorithm如果你從網(wǎng)上下載的代碼直接跑很可能會(huì)因?yàn)樗惴ㄟx項(xiàng)不兼容而報(bào)錯(cuò)。我建議盡量用MATLAB 2020b或更新版本并且把optimoptions的寫法統(tǒng)一。如果用的是別人的老代碼看到optimset就要留個(gè)心眼最好用optimoptions重寫一遍求解器選項(xiàng)。最后再說(shuō)一個(gè)個(gè)人習(xí)慣我做的每一次雙層優(yōu)化實(shí)驗(yàn)都會(huì)把隨機(jī)種子、參數(shù)表和結(jié)果圖打包存成一個(gè)文件夾命名為場(chǎng)景的描述比如“V2G_beta002_50EV”。這樣三個(gè)月后再回來(lái)看依然能清楚知道當(dāng)時(shí)做了什么?？蒲幸埠霉こ桃擦T可復(fù)現(xiàn)性是最大的生產(chǎn)力。希望這篇文章能讓你在電動(dòng)汽車雙層優(yōu)化調(diào)度的MATLAB路上少踩幾個(gè)坑先跑出第一版能收斂的代碼再慢慢打磨出自己的研究特色。