做電力系統(tǒng)調(diào)度的同學尤其是研究可再生能源和電動汽車入網(wǎng)方向的人大概率繞不開“協(xié)同調(diào)度”這四個字。我前陣子為了碩士論文復(fù)現(xiàn)把“可再生能源發(fā)電與電動汽車的協(xié)同調(diào)度策略”完整用Matlab走了一遍從場景生成、約束搭建到調(diào)用求解器踩了不少坑。這篇東西不是理論推導的復(fù)讀而是告訴你實際建模時怎么定義變量、怎么列約束、怎么讓求解器順利跑出結(jié)果以及哪些錯誤是你大概率也會遇到的。這篇內(nèi)容的適用對象很明確正在準備碩士論文里仿真章節(jié)的研究生、做微電網(wǎng)/配電網(wǎng)調(diào)度仿真的工程師還有那些裝了Matlab但不知道Yalmip到底怎么用的人。如果你只是想要一份能直接跑的代碼那這篇文章同樣能幫你搞清楚代碼每一步在做什么免得拿到開源代碼后連報錯都看不懂。1. 協(xié)同調(diào)度問題拆解定位核心模型1.1 問題邊界與決策變量先說清楚模型邊界。絕大多數(shù)碩士論文里的協(xié)同調(diào)度不會把電網(wǎng)全模型都搬進來而是把研究對象簡化成一個微電網(wǎng)或者一個等效的配電網(wǎng)節(jié)點。我做的是微電網(wǎng)場景里面有風電機組、光伏陣列、傳統(tǒng)燃煤機組、儲能系統(tǒng)還有一定數(shù)量的電動汽車。電動汽車既可以作為普通負荷充電也可以在電價合適的時候向電網(wǎng)放電V2G模式。這個模型要回答的問題是在一天24小時或96個時段內(nèi)每臺機組出多少電、儲能什么時候充什么時候放、電動汽車什么時候充電、允許多少車在某個時段放電才能讓總成本最低同時把可再生能源盡量用掉。決策變量分為連續(xù)變量和二進制變量。連續(xù)變量包括常規(guī)機組出力、儲能充放電功率、電動汽車充放電功率、向主網(wǎng)購電和售電功率。二進制變量則用于表示機組啟停狀態(tài)、電動汽車充放電狀態(tài)。這里有一個容易忽略的地方電動汽車不能同時充電又放電所以至少需要兩個二進制變量來分別代表充電狀態(tài)和放電狀態(tài)并且要加上互斥約束u_ch u_dis 1。很多初學者把模型定義得太復(fù)雜一上來就考慮每輛車的行程、用戶行為、充電樁位置結(jié)果變量數(shù)量爆炸求解器根本算不動。碩士論文復(fù)現(xiàn)階段我建議先把所有電動汽車聚合成一個等效儲能池比如100輛車每輛車40kWh那么聚合容量就是4000kWh聚合最大充電功率就是700kW。這樣做雖然忽略了個體差異但能先把調(diào)度邏輯跑通后續(xù)再擴展精細化模型。1.2 目標函數(shù)與成本建模目標函數(shù)看起來簡單寫起來容易亂因為成本項實在太多。常見的目標是最小化一天內(nèi)的系統(tǒng)總運行成本可以拆成五個部分常規(guī)機組燃料成本、可再生能源發(fā)電的運維成本通常很小、棄風棄光的懲罰成本、儲能和電動汽車的電池損耗成本以及從主網(wǎng)購電的成本減去向主網(wǎng)售電的收益。常規(guī)機組成本通常寫成輸出功率的二次函數(shù)C_G(P) a*P^2 b*P c。但二次函數(shù)會讓模型變成MIQP求解器求解速度明顯變慢。復(fù)現(xiàn)時我一般做分段線性化把功率范圍切成幾段每段用一個線性斜率表示這樣一來模型變成MILPGurobi和Cplex跑起來快得多。如果不做線性化遇到大場景數(shù)求解時間可能從幾分鐘變成一小時。棄風棄光懲罰成本很關(guān)鍵它決定模型會優(yōu)先消納可再生能源。懲罰系數(shù)一般取5001000元/MWh遠高于發(fā)電成本這樣模型寧可從電網(wǎng)買電也不會輕易棄風。但注意懲罰系數(shù)不要設(shè)成1e6這種量級否則數(shù)值病態(tài)會讓求解器報出各種奇怪的可行性問題。電動汽車V2G放電時電池會額外老化所以目標函數(shù)里要加放電懲罰項。基礎(chǔ)做法是給放電功率乘一個單位損耗成本比如0.2元/kWh。如果省掉這一項模型會為了賺峰谷電價差瘋狂放電產(chǎn)生不現(xiàn)實的調(diào)度結(jié)果。1.3 約束條件體系約束是調(diào)度模型的骨架少一條約束結(jié)果就可能是空中樓閣。最核心的是功率平衡約束所有電源常規(guī)機組、風電、光伏、儲能放電、EV放電、購電的功率之和必須等于所有負荷基礎(chǔ)負荷、儲能充電、EV充電、售電的功率之和。這條約束在每個時段都要滿足寫成向量形式時別忘了用sum(x,1)來確保按時間維度求和。儲能約束包括SOC狀態(tài)轉(zhuǎn)移方程、充放電功率上下限、SOC上下限以及一個很重要的“最終SOC等于初始SOC”約束否則模型可以把蓄電池的能量偷走造成成本偏低。SOC計算時注意單位一致性功率單位是kW儲能容量單位是kWh時間步長如果是1小時那么SOC變化量就是功率除以容量如果時間步長是15分鐘必須乘上0.25否則電量會差4倍。這個坑我見過不止一次。電動汽車約束比儲能多一層“充電需求約束”。每輛車每天至少需要充入一定的電量否則用戶沒電用。比如聚合體一天需要充入2000kWh那么所有時段充電總和減去放電總和后的凈充電量必須大于等于2000kWh。另外EV的SOC也要限定在10%到90%之間太深放電對電池不友好。最后是常規(guī)機組爬坡約束相鄰兩個時段的出力變化不能超過上限比如30MW/h。如果這一步漏了調(diào)度結(jié)果可能是機組出力劇烈波動看起來最優(yōu)但實際沒法執(zhí)行。爬坡約束本身很簡單就是abs(P_G(t)-P_G(t-1)) Ramp用Yalmip時寫成P_G(:,2:T)-P_G(:,1:T-1) Ramp和P_G(:,1:T-1)-P_G(:,2:T) Ramp兩條。2. 可再生能源不確定性處理從場景生成到魯棒/隨機優(yōu)化2.1 場景生成方法光伏和風電出力都有隨機性如果直接用期望值代替調(diào)度結(jié)果會很“天真”不考慮實際波動。處理隨機性最常用的方法是場景法生成大量代表可能出力曲線的樣本然后在樣本上求期望成本。場景生成不是隨便造數(shù)據(jù)。風電出力一般用威布爾分布描述風速再通過風功率曲線換算成出力光伏出力則用Beta分布描述光照強度。也可以用歷史出力數(shù)據(jù)直接采樣比如NREL或國內(nèi)電網(wǎng)公開的出力數(shù)據(jù)集。我在復(fù)現(xiàn)時用Matlab自帶的randn配合正態(tài)擾動生成200個初始場景每一步都加入時序相關(guān)性避免相鄰時段出力突變。一個容易被忽略的點場景不僅包含可再生能源出力還應(yīng)該包含負荷的波動尤其是電動汽車接入后的充電負荷。如果只考慮風電隨機而忽略負荷隨機論文審稿人一眼就看出來了。博士生不會太在意但碩士論文還是要把基礎(chǔ)的不確定性框架搭完整。生成場景時要注意樣本數(shù)量。500個場景做仿真求解時間不可接受20個場景則可能丟失尾部風險。我建議先做400個削減到15個這樣兼顧速度和精度。2.2 場景削減算法場景削減的核心思想是用少數(shù)具有代表性的場景替代原場景集且保持概率分布特性。常見的算法有三種K-means聚類、快速前向選擇FFS、同步回代消除SBR。K-means是最直觀的把400個場景聚成15類取每個類的中心作為代表場景概率為該類樣本比例。但K-means對初始聚類中心敏感聚類結(jié)果可能局部最優(yōu)。FFS和SBR是經(jīng)典的啟發(fā)式場景削減基于場景兩兩間的距離通常用歐氏距離逐步合并。SBR的做法是每次找一對總概率最小的場景刪除其中一個把它被刪的概率加到另一個上直到目標數(shù)量。這個方法能保留原場景集的分布形狀也比較容易手寫實現(xiàn)。我用Matlab手寫了一個SBR核心就是計算距離矩陣、更新概率。代碼量大約60行比調(diào)用復(fù)雜工具箱更可控。削減后要驗證計算削減前后場景集的均值和方差看差異是否在可接受范圍內(nèi)。差異太大說明目標場景數(shù)太少。下表是我常用的出力和負荷場景對比指標原始400場景削減后15場景誤差風電平均出力/MW12.412.83.2%光伏平均出力/MW8.17.92.5%負荷峰值/MW21.721.50.9%風電標準差3.22.99.4%標準差誤差稍高一點但對調(diào)度結(jié)果影響不大因為優(yōu)化的是期望成本均值誤差更關(guān)鍵。2.3 不確定性建模選擇處理不確定性的范式有隨機規(guī)劃和魯棒優(yōu)化兩種。隨機規(guī)劃相當于把多個場景都放進約束里目標函數(shù)是各場景成本的期望。實現(xiàn)方式就是在Yalmip里對每個場景分別建約束目標函數(shù)是sum(p_s*Objective_s)。這種方式簡單直接求解器能處理但需要保證每個場景下的約束都可行經(jīng)常會遇到某些極端場景導致問題不可行。魯棒優(yōu)化則是找一個最壞情況下的最優(yōu)決策。典型表示為min max_s。解決魯棒優(yōu)化通常需要引入對偶變量把內(nèi)層max問題轉(zhuǎn)化為約束條件數(shù)學推導繁瑣但所得方案在面對不確定時有更強的保障。碩士論文復(fù)現(xiàn)時如果原標題只寫了“協(xié)同調(diào)度”大概率是隨機優(yōu)化沒必要主動把難度升到魯棒。我當初直接選隨機優(yōu)化把精力放在如何設(shè)計約束保證所有場景可行比如引入可調(diào)變量如負荷削減來處理極端場景下的失負荷。3. Matlab代碼實現(xiàn)要點從建模到求解3.1 工具箱與求解器選擇Matlab里做優(yōu)化調(diào)度最推薦的方案是Yalmip商用求解器。Yalmip是一個建模工具箱它把模型翻譯成求解器能理解的內(nèi)部形式解完之后還能用value()把變量結(jié)果取出來。它支持Cplex、Gurobi、Mosek等。Gurobi學術(shù)版可以免費申請license速度非?？煊绕淝蠼釳ILP時比Matlab自帶的intlinprog強一個量級。Cplex現(xiàn)在對學術(shù)用戶也開放免費下載。如果不想申請也可以用免費的SCIP或CBC求解器但大規(guī)模問題會慢很多。安裝過程常出問題很多同學把Yalmip下載后塞到當前路徑但忘了把求解器路徑加進Matlab。正確的做法是在startup.m里寫入如下命令addpath(genpath(D:\yalmip)); addpath(genpath(D:\gurobi\win64\matlab)); gurobi_setup; savepath;然后運行yalmiptest看到所有求解器都顯示OK才算裝好。如果Gurobi后裝的記得重新savepath不然Matlab重啟后還得重新添加。求解器選擇上如果模型是純線性規(guī)劃沒有0-1變量可以用Gurobi的LP如果有0-1變量Gurobi用MILP算法默認會開并發(fā)割平面。設(shè)置選項時我習慣這樣寫solverOptions sdpsettings(solver,gurobi,verbose,2,gurobi.TimeLimit,300);設(shè)置時間限制很重要否則求解器可能陷入一個超大規(guī)模的MILP里跑幾小時。我自己復(fù)現(xiàn)時遇到過一次半天都沒跑完最后發(fā)現(xiàn)是二進制變量定義多了幾百個把-SOC互斥約束拆成每個EV每個時段兩個變量100輛車96時段就19200個二進制變量當然慢。后來聚合建模后變量降到幾百個幾秒就出結(jié)果。3.2 代碼框架設(shè)計與數(shù)據(jù)結(jié)構(gòu)不要把所有變量塞進一個大腳本。我用四個文件組織代碼main.m負責總流程和調(diào)用data_define.m里放系統(tǒng)數(shù)據(jù)model_define.m里放Yalmip建模solve_and_plot.m負責求解和圖形輸出。如果后續(xù)要加多場景再加一個scenario_generate.m。變量定義的核心是用sdpvar聲明連續(xù)變量、用binvar聲明0-1變量。比如N_gen 3; T 24; P_G sdpvar(N_gen, T, full); % 常規(guī)機組出力 P_ev_ch sdpvar(N_ev, T, full); % EV充電功率 P_ev_dis sdpvar(N_ev, T, full); % EV放電功率 u_ch binvar(N_ev, T, full); % EV充電狀態(tài) u_dis binvar(N_ev, T, full); % EV放電狀態(tài)這里full參數(shù)特別重要。如果不加Yalmip默認認為變量是個n*n方陣如果你傳入的是1x24向量它會報維度錯誤。很多新手在這里卡了一晚上。時間軸統(tǒng)一為行向量功率矩陣統(tǒng)一為變量數(shù) x T所有sum()操作明確指定維度sum(x,1)這樣約束構(gòu)建時不用反復(fù)轉(zhuǎn)置。3.3 核心約束構(gòu)建與目標函數(shù)先寫功率平衡約束。這里涉及一個常見問題可再生能源出力在場景化之后是一個概率矩陣如果你用期望值那就退化成了確定性模型。我用削減后的場景集每個場景單獨建約束所以模型里會有一系列帶場景下標(s)的變量。但在確定性驗證時可以先不考慮場景只用一個代表場景。model_define.m中的示意如下Constraints []; % 功率平衡約束 % sum(P_G,1)P_windP_pvP_St_dissum(P_ev_dis,1)P_buy ... % P_base P_St_ch sum(P_ev_ch,1) P_sell Constraints [Constraints, ... sum(P_G,1)P_wind_sceP_pv_sceP_St_dissum(P_ev_dis,1)P_buy ... P_base P_St_ch sum(P_ev_ch,1) P_sell];這里P_wind_sce、P_pv_sce是給定場景下的出力序列。如果允許棄風和棄光就把它們替換為P_wind_used并新增棄風變量P_wind_curtail約束為P_wind_used P_wind_curtail P_wind_sce目標函數(shù)中加入懲罰成本。EV充放電互斥約束Constraints [Constraints, u_ch u_dis 1]; Constraints [Constraints, P_ev_ch P_ev_ch_max * u_ch]; Constraints [Constraints, P_ev_dis P_ev_dis_max * u_dis];SOC約束要寫成遞推式% SOC_ev(t1) SOC_ev(t) eta_ch*P_ev_ch(t)/Cap_ev - P_ev_dis(t)/(eta_dis*Cap_ev) - E_drive(t)/Cap_ev for t 1:T-1 Constraints [Constraints, SOC_ev(t1) SOC_ev(t) ... eta_ch*P_ev_ch(:,t)/Cap_ev ... - P_ev_dis(:,t)/(eta_dis*Cap_ev) ... - E_drive(:,t)/Cap_ev]; end注意E_drive代表該時段電動汽車出行消耗的能量如果沒考慮出行可以設(shè)為零。目標函數(shù)建議用分段線性機組成本的增量形式來實現(xiàn)。簡單起見復(fù)現(xiàn)時可以直接用二次函數(shù)讓Gurobi當作MIQP求解但一旦場景數(shù)超過100MIQP會慢到讓人崩潰。所以我用增量成本分段線性化定義機組出力在第k段的變量P_seg(i,k,t)成本就是各段斜率之和。這樣目標函數(shù)保持線性。3.4 求解與結(jié)果輸出求解調(diào)用一行代碼result optimize(Constraints, Objective, sdpsettings(solver,gurobi,verbose,2));求解完之后第一件事是檢查狀態(tài)if result.problem 0 disp(求解成功); else yalmiperror(result.problem) end不要只依賴result.problem為0有時候求解器返回1警告但依然可以接受。yalmiperror能給出比較明確的錯誤描述。結(jié)果提取用value()函數(shù)。例如P_G_opt value(P_G); P_ev_ch_opt value(P_ev_ch); P_ev_dis_opt value(P_ev_dis); SOC_ev_opt value(SOC_ev);畫圖時要注意時序?qū)R。MATLAB的stairs更適合畫電力系統(tǒng)調(diào)度圖因為功率在時段內(nèi)保持恒定。使用subplot同時展示機組出力、EV充放電功率、儲能SOC和棄風棄光量讓調(diào)度策略一目了然。有經(jīng)驗的讀者都會再做一個“可行性檢查”把所有決策變量的最值打印出來看是否在合理區(qū)間內(nèi)。比如min(P_G_opt(:))如果出現(xiàn)負值那一定是有約束寫錯了。4. 復(fù)現(xiàn)過程中常見問題與排查技巧4.1 求解失敗與收斂慢最常遇到的情況是求解器提示“Infeasible”也就是約束之間互相矛盾。第一步檢查可行性把約束拆分逐個添加進模型測試。比如先只加功率平衡看看能不能解再加機組上下限直到定位到哪條約束導致不可行。還要注意數(shù)值尺度。通常功率單位用kW電價單位用元/kWhSOC是小數(shù)。如果目標函數(shù)里某一部分是0.001量級另一部分是1e8量級求解器數(shù)值精度會出問題。解決方法是統(tǒng)一量綱或者把大數(shù)值項除以一個基準功率。求解速度慢時先看二進制變量數(shù)量。如果二進制變量上千MILP求解時間指數(shù)增長。另一個常用技巧是設(shè)置“mipgap”也就是最優(yōu)間隙容忍度。調(diào)度問題允許1%的間隙結(jié)果差距不大但速度能快很多ops sdpsettings(solver,gurobi,gurobi.MIPGap,0.01);4.2 維度不匹配與變量定義錯誤Yalmip變量尺寸不一致是最容易報錯的地方。比如你定義了P_G sdpvar(N_gen, T)這在Yalmip里默認是一個方陣不對sdpvar(N_gen, T)創(chuàng)建的是N_gen行T列的矩陣沒有full也一樣。但如果你傳入了1作為第一個維度而T是24那么P_G是行向量后續(xù)約束里你可能以為它是列向量寫P_G(:,1)會報錯。最好的習慣是在變量定義前用assert(size(P_G,2)T)來校驗。還有個常見坑是binvar默認是不帶引用限制的。在不該用二進制的地方用了binvar就會把問題變成超大的MILP。比如表征“電動汽車只能在充電樁接入時充電”不需要用二進制變量可以用時間窗直接強制充電功率為0Constraints [Constraints, P_ev_ch(:, outside_window) 0];。4.3 場景削減后的結(jié)果偏差用場景法求期望成本結(jié)果往往比實際期望值偏樂觀因為削減后的場景覆蓋不到極端情況。解決辦法是做回代檢驗把求出的最優(yōu)解固定住代入原始的400個場景計算真實的期望成本和模型里的期望成本對比。如果偏差超過5%需要增加保留場景數(shù)或者換成更好的削減算法。另外有些論文里不是單獨做場景削減而是直接用蒙特卡洛在約束里循環(huán)所有場景這種做法在小規(guī)模系統(tǒng)里可行但場景數(shù)一大就完蛋。如果堅持這樣做至少要限制場景數(shù)在20以內(nèi)。4.4 電動汽車參數(shù)設(shè)置的坑EV參數(shù)設(shè)置直接影響調(diào)度結(jié)果。常見問題聚合容量過大導致EV放電能力比電網(wǎng)功率還大模型會安排大量V2G把儲能變成主力電源。實際中EV不可能全天候聯(lián)網(wǎng)所以要設(shè)定接入時段比如只允許18:00到次日8:00充放電。SOC的初值和終值必須一致否則模型會從初始SOC里“偷能量”。如果設(shè)置終值比初值低那么模型相當于免費用掉了電池里的電成本會偏低。我建議加約束SOC_ev(1) SOC_ev(T1)強制一個運行周期內(nèi)電量守恒。充電需求每天最少多需要仔細定義。如果設(shè)定所有EV一天要充入2000kWh而模型為了多放電把SOC壓到很低然后集中在下半夜充電雖然滿足約束但可能不現(xiàn)實。更穩(wěn)妥的做法是設(shè)定每個EV的入網(wǎng)時間和離網(wǎng)時間離網(wǎng)時SOC要達到某個下限比如90%。但這樣就引入了跟具體時間相關(guān)的復(fù)雜約束我建議在聚合模型里簡化為總凈充電量約束且最大放電功率不超過總負荷的30%確保V2G只是輔助手段。5. 擴展方向與個人踩坑記錄5.1 從單節(jié)點到配電網(wǎng)潮流如果論文需要研究線路約束和電壓分布那就不能只在單節(jié)點上做功率平衡。最簡單的是用DistFlow線性化LDF或二階錐松弛SOCP。DistFlow把功率、電壓、損耗關(guān)系寫成線性或近似線性的等式適合有載調(diào)壓的設(shè)備SOCP則把潮流約束寫成錐約束Yalmip能用cone或norm實現(xiàn)Gurobi支持二階錐。我在復(fù)現(xiàn)時最初沒有網(wǎng)絡(luò)約束后來加上了33節(jié)點配電網(wǎng)約束數(shù)量增加但性能還好。但需要提醒加入網(wǎng)絡(luò)約束后原本在單節(jié)點下可行的EV充放電策略在潮流模型下可能不可行尤其是節(jié)點電壓越限。需要把EV充放電功率和網(wǎng)絡(luò)邊界聯(lián)系起來。這一步如果做出來論文的含金量會上一個檔次。5.2 從單目標到多目標擴展碩士論文常喜歡寫“經(jīng)濟性環(huán)保性”協(xié)同優(yōu)化也就是多目標優(yōu)化。最簡單的是加權(quán)求和法和ε-約束法。加權(quán)求和法把碳排放和運行成本加權(quán)成一個目標函數(shù)但權(quán)重選取很主觀。ε-約束法則是把碳排放作為額外約束每取一個排放上限就求解一次得到帕累托前沿。在Matlab里用Yalmip配合循環(huán)就能輕松實現(xiàn)ε-約束法。我復(fù)現(xiàn)時先用加權(quán)法發(fā)現(xiàn)碳排放權(quán)重稍大調(diào)度就恨不得全部用EV放電成本劇增后來改用ε-約束法結(jié)果曲線更平滑。如果你也要寫多目標建議用ε-約束法并用帕累托圖展示。5.3 我的真實復(fù)現(xiàn)體驗這個項目我前后花了約三周每天兩小時。第一周卡在場景削減和變量維度上第二周模型跑通但結(jié)果不合常理比如EV瘋狂放電導致SOC一直為零第三天突然發(fā)現(xiàn)SOC賦初值時用了0.8約束卻限制SOC下限是0.9當然不可行。找到后真想拍桌子。個人最大的體會是一定要先從確定性模型開始把所有約束用單場景跑通再引入隨機場景。不要一開始就上Monte Carlo加場景削減否則報錯時你根本不知道是場景出了問題還是約束出了問題。我調(diào)試順序是無EV、無儲能、有儲能、有EV一層層加。這樣加一個部分就驗證一個部分效率最高。最后再分享一個小技巧所有功率和電量的單位用kW和kWh所有成本單位用元電價用元/kWhSOC用0到1的小數(shù)。寫代碼前先把單位寫在注釋里因為當你盯著滿屏數(shù)字找錯誤時單位不一致是最隱蔽的元兇。跑通了之后記得保存一份“運行正常”的備份再改下一步這是我和我身邊人用血淚換來的經(jīng)驗。