最近幫幾個研究生復現(xiàn)“可再生能源發(fā)電與電動汽車的協(xié)同調(diào)度策略研究”這類論文我最大的感受是卡住大家的往往不是數(shù)學模型本身而是怎么把論文里那幾頁公式變成能跑的Matlab代碼。今天就把一套我反復用過、也反復講給學生聽的復現(xiàn)思路完整整理出來內(nèi)容包括問題建模、求解器選型、仿真算例以及一堆網(wǎng)上查不到的小坑。如果你正準備復現(xiàn)碩士論文、搭建微電網(wǎng)調(diào)度模型或者只是想知道電動汽車如何在新能源波動時“挺身而出”這篇內(nèi)容應該能幫你省不少時間。這類題目的關(guān)鍵詞基本是固定的可再生能源發(fā)電、電動汽車、協(xié)同調(diào)度、Matlab。但關(guān)鍵詞之間怎么串成一個可計算的問題是大多數(shù)人第一步就沒邁過去的地方。我不打算只貼一段能跑的代碼而是把從題目到數(shù)學、從數(shù)學到代碼、再從代碼到結(jié)果驗證的完整鏈路講清楚。1. 為什么可再生能源和電動汽車要放進同一個調(diào)度模型1.1 EV是“會跑的儲能”協(xié)同調(diào)度的物理基礎(chǔ)風電和光伏的出力由天氣決定晚上負荷高峰時往往沒風沒光白天光照好但負荷可能還沒上來這就形成了棄風棄光與峰谷差拉大的雙重困境。電動汽車不一樣它本質(zhì)上是一塊帶輪子的電池。統(tǒng)計下來私家車一天至少有20個小時是停著的車里的電池閑著不用如果把這部分容量聚合起來效果相當于一座規(guī)模可觀的分布式儲能電站。更重要的是EV不僅能在谷時充電還能在峰時通過V2G把電反送回電網(wǎng)這就是“協(xié)同調(diào)度”最直接的物理支撐。但EV也不能當成無限容量的儲能隨便調(diào)度。它有自己的出行需求早上八點要出門電量就得滿足通勤白天停在公司樓下能不能充電要看樁位晚上回家后才是真正的可調(diào)度窗口。這些約束決定了建模時必須引入“在網(wǎng)時段”“離網(wǎng)SOC要求”“充放電互斥”等條件不能一股腦把所有EV當成一個恒定可調(diào)的大電池。1.2 協(xié)同調(diào)度要回答的三個信號問題所謂協(xié)同調(diào)度說到底是在回答三個問題EV在哪些時段增加充電、哪些時段減少充電、哪些時段反向放電。回答的準則不是“EV怎么充最省錢”而是“整個系統(tǒng)的運行成本和新能源消納效果最優(yōu)”。傳統(tǒng)機組、可再生能源、EV集群、上級電網(wǎng)四類資源要放在同一個優(yōu)化框架里統(tǒng)籌計算這跟單純的“有序充電”有本質(zhì)區(qū)別。有序充電通常只把削峰填谷當目標而協(xié)同調(diào)度要同時考慮出力分配、備用響應、充放電策略和網(wǎng)絡約束。舉個很直白的例子某小區(qū)晚上有500輛車同時接入如果大家回家就搶著充滿配電變壓器很容易直接被拉爆但如果調(diào)度中心錯開充電時間讓一部分車在夜間風電大發(fā)時再充另一部分車在早高峰放電賺錢變壓器壓力小了新能源棄電也少了EV用戶還能拿到激勵。這就是“協(xié)同”兩個字的現(xiàn)實價值。2. 建模之前必須想清楚的三個關(guān)鍵點2.1 風光預測誤差怎么進模型場景法還是魯棒法可再生能源出力的不確定性是無法繞開的。碩士論文里最常見的兩種處理方式第一種是場景法第二種是魯棒法。場景法把風速、光照的預測誤差看成隨機量用Monte Carlo采樣生成大量出力場景再用K-means或同步回代縮減成十幾個典型場景目標函數(shù)寫成場景集合下的期望成本。魯棒法則是構(gòu)造一個不確定集比如實際出力等于預測值加減偏差模型要保證最惡劣場景下系統(tǒng)依然不會失負荷。復現(xiàn)時我的建議是如果原論文沒寫明白優(yōu)先用場景法打底。場景法實現(xiàn)簡單、結(jié)果直觀Yalmip里用循環(huán)或者矩陣化約束都很方便跑通之后再往魯棒方向擴展比如從單層優(yōu)化改造成兩階段魯棒優(yōu)化用列約束生成算法CCG求解這樣文章的創(chuàng)新點也能自然帶上。另外要注意無論用哪種方式風光的切入范圍不是簡單的“0到預測值”而是要考慮預測誤差的分布特征否則模型會過度樂觀。2.2 目標函數(shù)不能只寫運行成本還要寫“約束的代價”目標函數(shù)通常是運行成本最小化但這幾年論文越來越喜歡加碳排放成本、棄風棄光懲罰、EV電池退化成本等項。以火電/微型燃氣輪機為例發(fā)電成本一般是二次函數(shù)C a·P2 b·P c。如果直接放進MILP框架這個二次項會讓模型變成MIQP。Gurobi和Cplex都能解MIQP規(guī)模不大時沒問題但我在復現(xiàn)時更推薦把二次成本做分段線性化這樣模型始終是MILP求解穩(wěn)定性更好也更容易被審稿人接受。還有一個很容易踩的坑是棄風棄光懲罰系數(shù)。如果這個系數(shù)設(shè)得太小求解器會寧可少發(fā)新能源也不愿意調(diào)整燃機出力結(jié)果算出來棄風率很高和論文結(jié)論完全對不上。懲罰系數(shù)要大于燃機邊際成本的差值才能讓新能源消納真正成為優(yōu)化的硬驅(qū)動。這也是很多人“代碼和論文結(jié)論對不上”的根本原因。3. 用MatlabYalmip把論文公式翻譯成可運行代碼3.1 為什么我用YalmipGurobi而不是直接寫優(yōu)化算法很多人拿到模型第一個念頭是用粒子群、遺傳算法去求解我的態(tài)度是能不用啟發(fā)式就不用。協(xié)同調(diào)度本質(zhì)上是帶整數(shù)變量的線性/二次規(guī)劃問題Yalmip加Gurobi這套組合可以從理論上保證全局最優(yōu)而且建模效率遠高于手寫單純形法、內(nèi)點法。Yalmip是Matlab下的一個免費建模層你只需要定義變量、寫約束、寫目標它會自動把模型轉(zhuǎn)成求解器能識別的標準形式Gurobi負責實際求解MILP/MIQP速度快、穩(wěn)定學術(shù)許可也容易申請。如果你的機器裝不了Gurobi退一步可以用Cplex或者開源求解器SCIP、CBC。CBC性能會弱一些但應付24小時的小微網(wǎng)算例綽綽余。還有一點經(jīng)驗Yalmip的變量定義要區(qū)分連續(xù)變量和二進制變量EV的充放電狀態(tài)、機組的啟停狀態(tài)都必須用binvar或integer算錯了就變成純線性規(guī)劃結(jié)果完全失真。3.2 變量定義和約束組裝的“套路”我寫這類代碼的固定套路是先畫矩陣維度圖再動手寫代碼。時間維度T取24機組編號N_gEV聚合體編號N_ev_agg。一般不建議逐輛EV建模而是把同一類充電特性、同一批出行時間的車聚合成一個“EV集群”再對這個集群建模。這樣做變量數(shù)量能減少幾個數(shù)量級求解速度快得多論文里也常這么寫。變量大體分為幾類燃機出力P_g、風電消納P_w、光伏消納P_pv、EV充放電功率P_ch/P_dis、電池SOC以及EG從上級電網(wǎng)購電功率P_buy。二進制變量包括EV充電狀態(tài)u_ch和放電狀態(tài)u_dis。約束則按功率平衡、機組上下限、爬坡、EV功率/SOC、風電光伏消納、聯(lián)絡線功率這幾類分別組裝。寫成代碼時優(yōu)先用矩陣切片避免深層的三重循環(huán)如果非要循環(huán)T24時用循環(huán)問題不大但語義要清楚。3.3 24小時日前調(diào)度的核心代碼骨架下面這段代碼是能跑通的最小骨架省略了部分燃機爬坡約束和場景循環(huán)但主結(jié)構(gòu)很清晰。數(shù)據(jù)部分我用了行向量格式方便和Yalmip變量維度對齊。%% 參數(shù) T 24; dt 1; N_g 2; % 兩臺微型燃氣輪機 N_ev_agg 1; % 一個EV集群內(nèi)部聚合100輛車 Ecap 100 * 24; % 集群總?cè)萘?kWh SOC0 0.5 * Ecap; % 初始SOC PchMax 100 * 3; % 最大總充電功率 kW PdisMax 100 * 3; % 最大總放電功率 kW SOCmin 0.2 * Ecap; SOCmax 0.9 * Ecap; eta 0.9; P_load [...]; % 1x24 負荷曲線 P_wf [...]; % 1x24 風電預測 P_pvf [...]; % 1x24 光伏預測 avail ones(1, T); % EV在網(wǎng)時段可按實際配置 %% 變量 P_g sdpvar(N_g, T, full); P_w sdpvar(1, T, full); P_pv sdpvar(1, T, full); P_ch sdpvar(N_ev_agg, T, full); P_dis sdpvar(N_ev_agg, T, full); SOC sdpvar(N_ev_agg, T1, full); u_ch binvar(N_ev_agg, T, full); u_dis binvar(N_ev_agg, T, full); P_buy sdpvar(1, T, full); %% 約束 Cons []; for t 1:T % 功率平衡 Cons [Cons, sum(P_g(:,t)) P_w(t) P_pv(t) P_dis(t) P_buy(t) ... P_load(t) P_ch(t)]; % SOC遞推 Cons [Cons, SOC(:,t1) SOC(:,t) (eta*P_ch(t) - P_dis(t)/eta)*dt/Ecap]; % 充放電功率上限avail為1時才可充放 Cons [Cons, 0 P_ch(t) PchMax * avail(t) * u_ch(t)]; Cons [Cons, 0 P_dis(t) PdisMax * avail(t) * u_dis(t)]; % 充放電互斥 Cons [Cons, u_ch(t) u_dis(t) 1]; % SOC上下限 Cons [Cons, SOCmin SOC(:,t1) SOCmax]; end Cons [Cons, SOC(:,1) SOC0]; %% 目標燃機成本 購電成本 棄風棄光懲罰 c_a [0.02; 0.02]; c_b [0.5; 0.6]; Objective sum(sum(c_a .* P_g.^2 c_b .* P_g)) ... sum(0.8 * P_buy) ... sum(15 * (P_wf - P_w)) sum(15 * (P_pvf - P_pv)); %% 求解 ops sdpsettings(solver, gurobi, verbose, 0); sol optimize(Cons, Objective, ops);這段代碼里我刻意把EV聚合體當成一個“大電池”來寫。你可能會問100輛車同時充放功率和SOC都是聚合值會不會丟失單車SOC信息這正是論文復現(xiàn)里的常見取舍。如果研究點是EV參與調(diào)度的策略聚合建模夠用如果研究點是每輛車的電池壽命差異那才需要逐車建模。復現(xiàn)之前先想清楚論文要回答什么問題避免模型過度復雜。4. 仿真算例怎么搭參數(shù)從哪來結(jié)果怎么驗4.1 一套能跑通的微網(wǎng)算例參數(shù)算例參數(shù)是復現(xiàn)中最大的“自由變量”。我常用的是一套經(jīng)典微型電網(wǎng)參數(shù)兩臺微型燃氣輪機一臺額定100kW、一臺額定80kW風電機組裝機100kW光伏裝機80kW基礎(chǔ)負荷峰值約400kW谷值約200kW。EV集群取100輛車單臺電池容量24kWh最大充放電功率3kW總數(shù)對應集群總?cè)萘?400kWh充放電總功率300kW。電價按峰谷分時峰時1.2元/kWh谷時0.4元/kWh。關(guān)鍵參數(shù)列在下表。參數(shù)數(shù)值說明燃機1額定/最小出力100 / 20 kW爬坡40 kW/h燃機2額定/最小出力80 / 15 kW爬坡30 kW/h風電裝機 / 預測峰值100 / 70 kW夜間出力偏大光伏裝機 / 預測峰值80 / 60 kW正午出力偏大負荷峰值 / 谷值400 / 200 kW典型日負荷EV集群車數(shù)100 輛聚合總?cè)萘?400 kWhEV最大總充/放電功率300 / 300 kW單車3kW聚合電池SOC范圍0.2 ~ 0.9保留出行電量充放電效率0.9往返約0.81峰谷電價1.2 / 0.4 元/kWh時段可按電網(wǎng)數(shù)據(jù)設(shè)這個量級的算例對Gurobi來說幾乎是秒解非常適合剛開始調(diào)代碼時使用。等代碼跑通后再逐步放大到IEEE 33節(jié)點配電網(wǎng)或者數(shù)百個EV集群重點考察求解時間。4.2 無序充電、有序充電、V2G三個場景的對比場景設(shè)計直接影響論文說服力。我一般會設(shè)三個場景場景一無序充電EV從18:00開始以最大功率連續(xù)充4小時不做任何優(yōu)化場景二有序充電EV在谷時充電但禁止放電場景三協(xié)同調(diào)度允許V2GEV可以在負荷高峰放電。三者的總運行成本和棄風率對比是整篇復現(xiàn)的核心圖表。用上面的參數(shù)跑完結(jié)果量級通常是這樣示意性數(shù)據(jù)不同論文參數(shù)會使絕對值不同場景總運行成本(元)棄風棄光率(%)峰值負荷(kW)無序充電8508.2510有序充電7403.6430協(xié)同調(diào)度(V2G)6800.8355從這個結(jié)果能清楚看到兩條結(jié)論一是EV參與調(diào)度后系統(tǒng)成本明顯下降二是夜間風電消納率大幅提升。更關(guān)鍵的是有序充電只是削峰V2G才是真正的“協(xié)同”——在負荷尖峰時把EV的電反送回去系統(tǒng)峰值負荷隨之降低。畫圖時用堆疊面積圖展示各機組出力再用階梯圖展示EV充電/放電功率論文質(zhì)感的提升非常明顯。4.3 結(jié)果合理性檢查清單很多同學算完直接截圖寫結(jié)論這是最容易翻車的環(huán)節(jié)。我建議跑完優(yōu)化后打印幾個關(guān)鍵量做一次“結(jié)果警察式”的檢查功率平衡殘差是否在10??以內(nèi)EV的SOC曲線是否始終處于限值內(nèi)離網(wǎng)時段有沒有充放電燃機出力是否滿足爬坡約束風電/光伏消納是否超過預測值棄風棄光懲罰項是否明顯改變了出力分配。用一行代碼就能提取并檢查功率平衡Pg value(P_g); Pw value(P_w); Ppv value(P_pv); Pc value(P_ch); Pd value(P_dis); Pbuy value(P_buy); balance sum(Pg,1) Pw Ppv Pd Pbuy - P_load - Pc; fprintf(最大功率平衡殘差: %.3e\n, max(abs(balance)));如果殘差大于10??第一件事不是調(diào)精度而是回去看維度和公式符號。功率平衡這個等式一旦寫錯后面所有結(jié)果都是廢的。5. 復現(xiàn)過程中最容易踩的五個坑5.1 求解器沒接好報錯全亂套Yalmip裝完之后第一件事是運行yalmiptest它會列出所有已識別求解器的狀態(tài)。如果Gurobi沒有出現(xiàn)在列表里optimize會提示“No appropriate solver”。最常見的原因是沒有把Gurobi的Matlab接口路徑加入當前工作區(qū)或者許可證沒配好。注意Gurobi的許可證和Matlab的許可證是兩套東西不要混在一起排查。學術(shù)版本用免費license安裝完用gurobi_setup或手動addpath到gurobi的matlab目錄即可。這里有個我踩過不止一次的教訓Yalmip的版本和Gurobi版本存在兼容性差異舊版Yalmip調(diào)用新版Gurobi偶爾會報奇怪的“Output argument not assigned”錯誤。解決辦法是先升級Yalmip再檢查Gurobi這個順序不要反。5.2 別讓模型變成MINLP雙線性項和二次項處理協(xié)同調(diào)度模型里最容易出現(xiàn)雙線性項的地方是連續(xù)變量和二進制變量相乘。比如你希望通過一個0-1變量表示“EV只有晚上才能放電”于是寫了一行P_dis(t) PdisMax * u_dis(t) * disrupt(t)其中disrupt(t)是另一個連續(xù)變量這就形成了雙線性約束模型變成難以求解的MINLP。正確做法是把0-1變量當成開關(guān)用不等式來限功率而不是讓連續(xù)變量和二進制變量在乘號里直接相見。二次成本項也同樣道理。Gurobi雖然能解MIQP但大規(guī)模場景下MIQP的求解時間明顯比MILP長而且非凸二次規(guī)劃容易出現(xiàn)局部最優(yōu)問題。我通常用分段線性約束把二次函數(shù)逼近成線性逼近誤差控制在1%以內(nèi)求解速度能快一個數(shù)量級。這不是炫技而是工程實踐里很務實的選擇。5.3 SOC初值和“無解”的排查順序無解infeasible是新手最頭疼的問題。我自己的排查順序是先去掉SOCmin/SOCmax只保留初值約束看模型能不能跑通能跑通就說明問題出在SOC上下限與充放電功率不匹配。再逐步加回約束每加一組就跑一次直到哪一組加上后報無解問題就鎖定在那組約束上。舉一個真實例子論文要求EV早上8點離家時SOC要達到0.9但充電時段只有凌晨2點到6點4小時乘以充電功率上限根本補不上電量模型當然無解。碰到這種情況要么調(diào)整初始SOC要么放寬離家SOC要求要么把充電功率上限提高必須有人為干預。Yalmip的check(Cons)函數(shù)能輸出每條約束的殘差無解時它會告訴你哪條約束的“違規(guī)程度”最大方向感一下就出來了。5.4 一維二維矩陣方向能讓你找bug找半天Matlab矩陣維度方向是這類代碼的隱形殺手。sdpvar(N_g, T, full)生成的是N_g行T列變量如果你用P_load是1行T列而sum(P_g,1)也是1行T列兩者相加沒問題但有時你從Excel讀入數(shù)據(jù)后P_load是T行1列直接相加就會維度不匹配。Yalmip在這類維度錯誤上通常不會給特別明確的提示只告訴你“Dimension mismatch”。我的習慣是全部數(shù)據(jù)統(tǒng)一成行向量也就是1×T并且在每條約束后面用size打印一次維度做斷言。這樣雖然看起來繁瑣但能避免90%的隱性bug。另外binvar(N_ev_agg, T)生成的變量也是行數(shù)N、列數(shù)T和連續(xù)變量的維度規(guī)則完全一樣別在這種細節(jié)上鉆牛角尖。5.5 原論文參數(shù)缺失怎么辦碩士論文最讓人頭疼的問題就是參數(shù)不完整很多數(shù)據(jù)寫著“見文獻[xx]”就等于沒說。復現(xiàn)時千萬不要編一個自認為合理的數(shù)悄悄填進去否則后面審稿或者導師一問就露餡。正確的做法是用公開的典型算例參數(shù)或者從同一研究方向的英文期刊論文里把參數(shù)摘出來并在論文原文里注明參數(shù)來源。最常用的數(shù)據(jù)來源是IEEE標準算例、MATPOWER自帶數(shù)據(jù)以及一些知名綜述里匯總的參數(shù)表。如果論文里缺失的是比較敏感的參數(shù)比如電池退化成本系數(shù)我建議做敏感性分析把系數(shù)從0.1取到0.5每檔跑一次畫出EV放電量或總成本隨系數(shù)變化的關(guān)系曲線。這樣既不掩蓋參數(shù)不確定性又展示了模型的穩(wěn)健性導師和審稿人都喜歡這種處理方式。6. 代碼跑通之后還能往哪些方向擴展6.1 從單時段開環(huán)到兩階段滾動調(diào)度日前調(diào)度是一次性的開環(huán)決策把24小時一次性算完。但實際運行中風光預測每4小時甚至每小時都會更新一次性的計劃根本來不及應對誤差。所以代碼跑通后第二個值得做的升級是兩階段調(diào)度第一階段做日前機組組合第二階段做實時經(jīng)濟調(diào)度誤差通過EV和燃機爬坡來彌補。用MPC滾動優(yōu)化的方式滾動更新每次只執(zhí)行下一小時指令可以明顯看到系統(tǒng)對預測誤差的應對能力。這個擴展在代碼層面不需要重構(gòu)太多只要把目標函數(shù)從單時段改成窗口式滾動加一個實時場景生成器就行。很多碩士論文的亮點就落在“日前實時”的協(xié)調(diào)上從復現(xiàn)走向創(chuàng)新這一步往往是分水嶺。6.2 從成本最小到多目標折中如果原論文只有單目標你可以試著把碳排放作為第二目標用ε-約束法或者帕累托前沿方法處理。具體做法是先把碳排放最小化跑一遍得到碳排放的最小值再把這個最小值當約束放回成本最小化模型不斷放松碳排放上限帶回一簇帕累托解。最終畫出一條成本和碳排放的折中曲線這會比單點結(jié)果有說服力得多。加權(quán)求和也能做但權(quán)重系數(shù)主觀性太強審稿人可能會問“為什么取0.7和0.3”。用帕累托前沿展示的是整體權(quán)衡關(guān)系理論上更站得住腳。代碼實現(xiàn)上也無非是外層增加一個循環(huán)內(nèi)層對碳排放約束加參數(shù)不會復雜到哪里去。如果讓我重新復現(xiàn)一次這類論文我會堅持兩個習慣第一每個約束后面都注釋對應論文的公式編號寫代碼就像在寫公式表回頭改模型的時候會特別舒服第二先跑一個不含EV的版本再加EV、再加不確定性一步步逼近論文模型這樣每一步出問題都能立刻定位。這兩點救過我很多次也實實在在幫你把“復現(xiàn)”變成“理解”。