優(yōu)化調度Matlab實現(xiàn))
做電力系統(tǒng)調度的朋友應該都有體會論文里把模型寫清楚是一回事真要在Matlab里把代碼跑起來、把結果調合理完全是另一回事。尤其帶安全約束的優(yōu)化調度約束一多矩陣動不動就幾百上千行稍不注意就出現(xiàn)“約束沖突無解”或者“結果奇奇怪怪但不知道哪里錯了”的情況。我之前正好完整復現(xiàn)過一個“計及N-k安全約束的含光熱電站電力系統(tǒng)優(yōu)化調度模型”支持IEEE14節(jié)點和IEEE118節(jié)點兩種算例系統(tǒng)這里把建模思路、關鍵代碼邏輯和踩坑經驗一并寫下來給大家做個參考。這篇內容適合正在做電力系統(tǒng)優(yōu)化調度、新能源消納、光熱電站調度策略相關課題的研究生也適合剛接觸N-k安全約束、想用Matlab快速實現(xiàn)一個可跑通算例的工程師。文章不會從拉格朗日松弛講起而是從“這個模型到底在解決什么問題”“約束是怎么一步步加進去的”“代碼該怎么組織”這幾個角度展開盡量讓每個細節(jié)都可落地。1. 問題定位與整體思路為什么偏偏是N-k約束和光熱1.1 調度模型里加N-k約束本質上是在買“保險”傳統(tǒng)經濟調度一般只考慮系統(tǒng)正常運行狀態(tài)約束是潮流方程、機組出力上下限、爬坡速率這些。但電網實際運行里線路或者發(fā)電機組隨時可能跳閘如果調度結果在某個元件故障后直接導致潮流越限甚至電壓失穩(wěn)那這個調度方案在工程上是不合格的。N-k安全約束的含義是系統(tǒng)在任意k個元件退出運行后仍能通過調整發(fā)電出力維持安全運行。k1就是N-1是最常見的校核標準k2就是N-2對應更嚴酷的雙重故障場景。把這類約束嵌入優(yōu)化調度模型得到的調度方案就不是“最優(yōu)但脆弱”的方案而是“犧牲一點經濟性、換取較強魯棒性”的方案。這里要特別注意N-k約束不是簡單地把故障場景下的潮流約束加入模型就完事。對于大規(guī)模系統(tǒng)全部N-2場景的數(shù)量是組合爆炸級別的118節(jié)點系統(tǒng)全枚舉N-2根本不現(xiàn)實。所以實際代碼里一般會采用“預想事故集生成 關鍵故障篩選 主子問題迭代”的思路這個后面展開講。1.2 光熱電站在這類模型里的定位是什么光熱電站CSP跟光伏、風電不一樣的地方在于它自帶儲熱系統(tǒng)。光伏和風電出力“靠天吃飯”不具備可調度性光熱電站可以通過儲熱罐把白天多余的熱能存下來等到晚間或者負荷高峰期再釋放出來發(fā)電等于給了調度員一個靈活調節(jié)的杠桿。在計及N-k安全約束的模型里光熱電站電力系統(tǒng)的作用體現(xiàn)在兩個層面。第一故障發(fā)生后系統(tǒng)出現(xiàn)功率缺額光熱電站的儲熱系統(tǒng)可以在短時間內快速增加出力充當旋轉備用和事故備用的補充。第二光熱電站的調節(jié)能力可以讓其他常規(guī)機組在調度點位上留出更多裕量間接降低了故障場景下需要調整的幅度。所以建模時不能把光熱電站簡單當成一個“帶上下限的電源”至少要考慮光場集熱、儲熱罐充放熱、發(fā)電功率轉換這三段環(huán)節(jié)的能量平衡關系儲熱罐的容量和充放熱速率都要作為約束進入模型。1.3 為什么選IEEE14和IEEE118節(jié)點系統(tǒng)做驗證IEEE14節(jié)點系統(tǒng)規(guī)模小約束數(shù)量和變量數(shù)量都處于“人腦還能追蹤”的范圍適合調試模型邏輯。比如N-1場景總共也就幾十個N-2場景也還能枚舉方便檢查約束到底有沒有寫錯。IEEE118節(jié)點系統(tǒng)則更接近實際規(guī)模線路和機組數(shù)量上來了求解難度和計算時間也上來了適合做算法效率和實用性的驗證。我自己實測下來同樣一套代碼IEEE14節(jié)點用默認參數(shù)求解只要幾秒IEEE118節(jié)點可能要幾十秒甚至幾分鐘具體取決于N-k的k值、故障篩選策略還有求解器配置。這兩套系統(tǒng)一組合既能驗證小系統(tǒng)的正確性又能檢驗大系統(tǒng)的求解效率作為論文算例或者工程預研是夠用的。2. 數(shù)學模型與約束構建目標函數(shù)、決策變量和關鍵約束2.1 決策變量怎么設計不僅有機組出力還有光熱電站的儲放熱狀態(tài)整個調度模型是一個混合整數(shù)線性規(guī)劃問題。決策變量分三類第一類是常規(guī)火電機組的開停狀態(tài)和出力水平第二類是光熱電站的發(fā)電功率、集熱功率、儲熱罐充放熱功率第三類是N-k故障場景下各節(jié)點電壓相角、線路潮流、負荷削減量這些狀態(tài)變量。需要強調一點故障場景下的變量跟正常場景下的變量是有耦合關系的。常規(guī)機組在故障后的調整出力要在爬坡約束允許的范圍內光熱電站的儲放熱功率調整也要受儲熱罐實際儲熱量的限制。這個耦合關系如果漏掉計算出來的“故障后可調范圍”就沒有任何物理意義。具體到變量定義我用的是這樣的邏輯以IEEE118節(jié)點為例% 正常場景變量 Pg sdpvar(nGen, T); % 常規(guī)機組出力 on binvar(nGen, T); % 機組開停狀態(tài) Pcsp sdpvar(nCsp, T); % 光熱電站發(fā)電功率 Pstore sdpvar(nCsp, T); % 儲熱罐充放熱功率正放負充 % N-k故障場景變量 Pg_f cell(nContingency, 1); % 各故障場景下的機組出力 theta_f cell(nContingency, 1); % 各故障場景下的節(jié)點相角2.2 目標函數(shù)運行成本、棄光懲罰和負荷削減懲罰怎么權衡目標函數(shù)我采用了最小化總運行成本的模式包括四部分常規(guī)機組燃料成本、開停機成本、光熱電站運行維護成本、棄光懲罰成本以及故障場景下的負荷削減懲罰成本。這里引入一個權重參數(shù)將正常場景成本與故障場景懲罰成本統(tǒng)一量綱實際調試中發(fā)現(xiàn)懲罰系數(shù)設置在100到1000之間比較合適太低會導致模型“寧可切負荷也不調整出力”太高則會讓模型過于保守。燃料成本一般用二次函數(shù)表示在MILP里需要分段線性化% 分段線性化燃料成本 for i 1:nGen for t 1:T for s 1:nSeg g_seg(i,t,s) sdpvar(1,1); Cost Cost coef(i,s) * g_seg(i,t,s); end end end每段出力有上下限約束所有段出力之和等于機組總出力這樣就保證了開停機狀態(tài)下成本的連續(xù)性和凸性。剛開始做的時候很多人圖省事直接用線性成本結果調度結果會頻繁出現(xiàn)機組在最小技術出力附近“顛簸”的現(xiàn)象就是因為沒有處理好成本函數(shù)的凸性。2.3 常規(guī)約束功率平衡、爬坡、機組最小開停機時間功率平衡約束是每個時段所有機組出力加光熱出力之和等于負荷這是最基本的約束。爬坡約束要區(qū)分正常場景和故障場景正常場景爬坡約束和故障場景爬坡約束是獨立的因為故障發(fā)生之后機組調整出力要在一個時間段內完成不能瞬時跳變。最小開停機時間約束屬于“寫了容易調試麻煩”的類型?；跔顟B(tài)變量的遞推關系可以用一組線性不等式表示但時間粒度越小約束數(shù)量越多。我的代碼里默認T24時段最小開停機時間設置為2小時這個約束對求解時間影響不算大但如果T擴大到96時段這個約束的數(shù)量會直線上升需要特別注意。2.4 光熱電站建模三段式能量平衡不能省光熱電站的建模是核心難點也是最容易出錯的地方。我采用的是“光場—儲熱—發(fā)電”三段式模型。第一步光場接收太陽輻射能轉化熱能。這里用DNI數(shù)據(jù)乘上集熱面積和效率得到集熱功率這個功率上限受光照條件限制是外部給定參數(shù)。第二步熱能進入儲熱罐。儲熱罐相當于一個能量緩沖池能量平衡關系是儲熱罐當前儲熱量 上一時段儲熱量 充熱功率 - 放熱功率 - 熱損耗其中熱損耗通常簡化成儲熱量的固定百分比。第三步儲熱罐放出熱能驅動汽輪機發(fā)電。發(fā)電功率等于放熱功率乘以熱電轉換效率并且有最小技術出力和最大出力限制。% 儲熱罐能量平衡 for t 2:T Qstore(t) Qstore(t-1) eta_in * Pstore_in(t) - ... (1/eta_out) * Pstore_out(t) - loss * Qstore(t-1); end這里有個細節(jié)充熱功率和放熱功率在數(shù)學上應該是兩個非負變量不能用一個有符號變量代替否則優(yōu)化器會利用“同時充放熱”的偽操作來打擦邊球導致儲熱罐實際沒有存儲熱量但模型計算有收益。嚴格來說要加互斥約束或者使用兩變量建模我為了省事是加了一個隱性約束讓充放熱不同時為正實測下來可以用。2.5 N-k安全約束的數(shù)學表達和場景生成這是整個模型最核心的部分。N-k約束的本質是對每一個預想故障場景存在一個可行的調整方案讓系統(tǒng)在故障后仍滿足潮流約束、機組出力約束和爬坡約束。我把故障場景分成了兩類處理。第一類是線路故障對應節(jié)點導納矩陣變化第二類是發(fā)電機故障對應機組出力上限變?yōu)?。兩種故障的處理方式不太一樣線路故障會改變潮流分布發(fā)電機故障主要造成功率缺額。N-k場景生成的關鍵問題是數(shù)量控制。IEEE14節(jié)點下所有線路和發(fā)電機的N-2組合數(shù)量不算特別多可以全枚舉IEEE118節(jié)點下就必須做篩選。我的做法是先用正常場景的最優(yōu)潮流結果計算各線路的負載率找出負載率排名前N條的線路作為關鍵故障候選集再在候選集內枚舉k重組合。這個策略在工程上比較常用好處是計算量可控壞處是理論上存在漏掉某些極端組合的風險這也是近似算法的通病。故障場景下的約束表達最直接的是直流潮流方程% 故障場景下的潮流約束 for k 1:nCtrl Pgf Pg_f{k}; thetaf theta_f{k}; % 節(jié)點功率平衡 Bf * thetaf Pgf - Pd Pcurtail_f{k}; % 線路潮流限值 -Fmax Bf_line * thetaf Fmax; end直流潮流是線性約束可以直接放進MILP求解。如果要用交流潮流問題是高度非線性的正常模型直接變NLP求解難度不在一個量級。我這套代碼選擇了直流潮流主要原因是為了保證模型的線性結構在論文里說明這一點是重要的前提假設。2.6 為什么說這本質上是“魯棒優(yōu)化”的雛形把N-k約束和優(yōu)化目標放在一起看你會發(fā)現(xiàn)這個模型其實是在構造一個“正常場景最優(yōu)所有故障場景可行”的解。這和魯棒優(yōu)化的思路很像只是這里的不確定集是離散的故障場景集合而不是連續(xù)的不確定參數(shù)區(qū)間。理解這一點對后續(xù)擴展很重要。比如你想把N-k約束從“事后校驗”變成“在線滾動優(yōu)化”需要引入Benders分解或者列約束生成算法把故障場景約束作為子問題迭代添加到主問題中。實際上MATLAB代碼里如果能用YALMIP的魯棒優(yōu)化工具箱可以一定程度上自動化這個過程但##更底層的#邏輯還是需要自己懂。3. 求解方案與Matlab實現(xiàn)路徑從工具箱選型到代碼結構3.1 求解器選擇Gurobi還是CPLEXYALMIP怎么配這個模型是MILP問題需要用商業(yè)求解器。我在代碼里使用了YALMIP作為建模語言分別測試過Gurobi和CPLEX兩個求解器。實測結果求解器IEEE14節(jié)點求解時間IEEE118節(jié)點求解時間穩(wěn)定性Gurobi0.8s約40s高CPLEX1.2s約55s高兩個求解器都能收斂到相同的最優(yōu)值Gurobi整體快一些尤其是大系統(tǒng)下優(yōu)勢更明顯。YALMIP的好處是切換求解器只需要改一個參數(shù)ops sdpsettings(solver, gurobi, verbose, 2); % 或者 ops sdpsettings(solver, cplex, verbose, 2);3.2 代碼整體結構數(shù)據(jù)準備、建模、求解、分析四層分離寫這類優(yōu)化代碼最忌諱把所有邏輯堆在一個腳本里我自己的代碼結構是這樣組織的main_IEEE14.m/main_IEEE118.m入口腳本負責加載數(shù)據(jù)、設置參數(shù)、調用核心函數(shù)、展示結果data_IEEE14.m/data_IEEE118.m數(shù)據(jù)定義包括節(jié)點數(shù)據(jù)、線路數(shù)據(jù)、機組參數(shù)、光熱電站參數(shù)、負荷曲線build_model.m核心建模函數(shù)輸入系統(tǒng)數(shù)據(jù)和參數(shù)輸出YALMIP優(yōu)化模型solve_model.m調用求解器并處理結果plot_results.m可視化輸出這個結構的優(yōu)點是換算例系統(tǒng)時只需要改數(shù)據(jù)文件建模和求解代碼完全復用調試時也可以在build_model里設置斷點逐段檢查約束矩陣是否生成正確。build_model.m的核心流程大致是function model build_model(data, params) % 1. 定義變量 % 2. 構建目標函數(shù) % 3. 添加常規(guī)約束 % 4. 添加光熱電站約束 % 5. 生成N-k故障場景 % 6. 添加故障場景約束 % 7. 返回模型 end3.3 用YALMIP表達分段線性化細節(jié)決定成敗分段線性化的實現(xiàn)看起來簡單實際有幾個容易踩的坑。第一個坑是斷點坐標要包含機組的最大/最小技術出力點否則線性化出來的成本函數(shù)在端點附近不連續(xù)。第二個坑是每段的斜率必須是遞增的這樣才能保證成本函數(shù)是凸函數(shù)。如果原始數(shù)據(jù)不滿足這個條件需要先進行數(shù)據(jù)預處理用凸包近似處理。% 分段線性化 function cost_seg piecewise_linear(P, breakpoints, slopes) % P: 機組出力sdpvar % breakpoints: 分段點 % slopes: 各段斜率 end3.4 迭代求解N-k場景約束避免一次性全部塞入如果一開始就把所有N-2場景都加入模型IEEE118系統(tǒng)的約束矩陣會膨脹到難以承受。我的做法是采用兩階段迭代% 主循環(huán) scenario_set initial_scenarios; while ~converged % 求解當前場景集下的優(yōu)化模型 [Pg, obj] solve_master(scenario_set); % 校驗全部候選場景尋找被違反的故障場景 [violated, info] check_all_scenarios(Pg); if isempty(violated) converged true; else % 將違反最嚴重的場景加入場景集 scenario_set [scenario_set; worst_case(violated, info)]; end end這個方法的核心思想是“不需要一次性考慮所有約束只在迭代過程中逐步補上被違反的約束”。實測下來有些在初始集合里看起來不會出問題的場景在加入后反而被排除最終收斂需要的迭代次數(shù)大概在3到8輪之間。4. 算例設置與結果分析IEEE14和IEEE118節(jié)點4.1 數(shù)據(jù)準備光熱電站接入哪個節(jié)點有講究IEEE14節(jié)點系統(tǒng)的原始數(shù)據(jù)里沒有光熱電站需要自己選節(jié)點接入。我的做法是把光熱電站當成一個可調電源接到某節(jié)點同時把該節(jié)點的原負荷適當調增。選擇接入節(jié)點的原則有兩個一是盡量靠近負荷中心這樣可以減少故障場景下遠距離送電的壓力二是選擇對系統(tǒng)穩(wěn)定影響相對小的節(jié)點避免光熱電站接入后出現(xiàn)無功支撐不足的問題。IEEE118節(jié)點系統(tǒng)我接入了兩座光熱電站分別位于電氣距離較遠的兩個區(qū)域這樣能夠觀察光熱電站對整個系統(tǒng)不同區(qū)域安全裕量的影響。4.2 關鍵結果對比有/無N-k約束、有/無光熱電站為了把模型效果展示清楚我這里列舉一組典型的IEEE118節(jié)點算例結果均為T24時段場景正常場景總成本是否含N-1安全裕量是否含N-2安全裕量無N-k約束、無光熱100.0基準部分場景越限大量場景越限含N-k約束、無光熱108.5全部通過部分通過含N-k約束、含光熱106.8全部通過全部通過表格里的成本做了歸一化?！盁oN-k約束、無光熱”的方案成本最低但調度方案在N-1故障場景下就有線路潮流越限的情況加了N-k約束之后成本上升約8.5%換來了可靠性。引入光熱電站之后成本只比無光熱版本上升了約6.8%卻額外通過了N-2校核。這說明光熱電站的儲熱能力讓系統(tǒng)在應對雙故障時有了更大的出力調節(jié)窗口。4.3 光熱電站出力的“削峰填谷”能力如何體現(xiàn)對比有無光熱電站的調度結果能明顯看到光熱電站出力的時間平移特性。白天9:00-16:00光照充足光熱電站一邊發(fā)電一邊往儲熱罐里充熱傍晚17:00-21:00負荷高峰期儲熱罐釋放熱量光熱電站出力上升到接近額定功率正好替代了一部分高成本火電出力。這個時間平移特性在故障場景下也發(fā)揮了作用。夜間光照為零光熱電站無法直接集熱但儲熱罐有存量依然可以支撐短時間的高出力。調度模型需要考慮這個特點避免在白天就提前把所有儲熱耗盡。4.4 故障篩選策略的效果驗證針對IEEE118節(jié)點系統(tǒng)我對比了“全枚舉N-1 候選集N-2”和“僅全枚舉N-1”兩種策略的差異。全枚舉N-1一共涉及約186個預想場景候選集N-2由負載率前15的線路兩兩組合產生大約105個場景。兩種策略在最終調度結果上的差異很小成本差異不到1%但求解時間差異很大。全枚舉N-2場景的話模型規(guī)模直接翻幾倍求解時間可能超過半小時用候選集策略能把總求解時間控制在5分鐘以內。這再次說明在大規(guī)模系統(tǒng)里無腦枚舉所有N-k場景是不現(xiàn)實的“關鍵場景篩選迭代校核”才是工程上可行的路徑。5. 調試經驗與避坑指南那些花了三天才繞過去的坑5.1 約束沖突模型“無解”時的定位方法模型一上來就報infeasible是最常見的問題。這時候不要急著刪約束用YALMIP的optimize返回值逐步排查。diagnostics optimize(constraints, objective, ops); if diagnostics.problem ~ 0 % 獲取沖突信息 check(constraints); endcheck(constraints)會逐一輸出每個約束的殘差值殘差顯著非零的那個約束就是infeasible的來源。我調試時最常遇到的問題是儲熱罐能量平衡約束寫錯符號充熱當成放熱導致能量不平衡約束與發(fā)電功率約束互相矛盾。5.2 數(shù)值范圍問題量綱不統(tǒng)一導致求解器“發(fā)瘋”電力系統(tǒng)數(shù)據(jù)里功率是MW級別成本是$/MWh儲熱罐存儲量是MWh如果不注意量綱統(tǒng)一約束矩陣里會出現(xiàn)系數(shù)相差好幾個數(shù)量級的行導致求解器數(shù)值穩(wěn)定性極差。我的處理方式是在數(shù)據(jù)準備階段統(tǒng)一單位所有功率統(tǒng)一為MW能量統(tǒng)一為MWh成本統(tǒng)一為$。另外需要對大數(shù)和小數(shù)做一個縮放比如爬坡率參數(shù)在歸一化后一般在0.1到0.5之間不在這個范圍內的參數(shù)就要仔細檢查。5.3 故障場景下的爬坡約束容易忽略的隱含耦合故障場景下的機組出力調整不是無條件的。假設調度方案里某機組正常出力是300MW故障后要把出力升到350MW這50MW的調整必須在一定時間內完成受到機組爬坡速率的限制。我之前漏了故障場景下的爬坡約束只加了故障后的出力上下限結果系統(tǒng)給出的故障后調整方案在物理上根本無法實現(xiàn)。加上故障場景爬坡約束之后模型才真正變得“可執(zhí)行”。5.4 直流潮流的局限性線路有功越限只是“必要條件”最后必須提醒直流潮流的N-k校核是一個線性近似只檢查有功潮流越限忽略了無功功率和電壓約束。在實際工程里線路有功在N-1范圍內系統(tǒng)也可能因為電壓失穩(wěn)或者無功不足出問題。所以這套模型更適合用在對安全約束的“初步篩選”上真正的N-k安全校核還需要用交流潮流或者更精細的穩(wěn)態(tài)/動態(tài)分析工具來驗證。6. 從復現(xiàn)到擴展這套模型還能往哪些方向延伸6.1 從確定性N-k到概率性故障校核傳統(tǒng)N-k約束假設每個故障場景等概率發(fā)生沒有區(qū)分不同元件的故障概率。實際上老舊線路的故障概率遠高于新線路發(fā)電機的故障概率也受服役年限、運行狀態(tài)影響。一種自然的擴展是給每個故障場景賦一個概率權重把安全約束改造成機會約束Pr(系統(tǒng)安全運行) ≥ 1 - ε這會讓模型從MILP變成更復雜的混合整數(shù)非線性問題求解難度大幅上升但在工程決策上有意義。6.2 從離線調度到滾動優(yōu)化當前模型是離線求一個24小時調度方案。實際運行中負荷預測和光熱電站的DNI預測會有誤差離線方案直接拿來用可能效果不理想。一個可行的擴展方向是改成模型預測控制或者滾動時域優(yōu)化每隔一小時滾動更新一次未來24小時的調度方案每次求解時只執(zhí)行第一個時段的決策結果。此時N-k約束仍然可以在每個滾動窗口內嵌入只是在計算時間上要更摳門。6.3 從單目標到多目標經濟性和安全性的帕累托權衡這個模型本質上是把經濟性和安全性通過目標函數(shù)的加權系數(shù)綁定在一起了。如果想評估兩者之間的真實權衡關系可以改成多目標優(yōu)化框架求解帕累托前沿再結合決策偏好選擇合適的調度方案。IEEE14節(jié)點系統(tǒng)因為規(guī)模小做多目標優(yōu)化的計算成本不高很適合作為算法研究的測試床。6.4 與新能源不確定性建模的深度融合光熱電站雖然自帶儲熱但它的集熱功率受DNI顯著影響本質上仍是一個不確定電源。當前模型用確定性DNI曲線沒有考慮預測誤差。如果想把光熱電站的不確定性也納入模型可以在N-k安全約束的基礎上疊加蒙特卡洛場景模擬或者用分布魯棒優(yōu)化構造DNI不確定集。這類擴展跟N-k約束并不沖突反而可以統(tǒng)一在一個“雙層不確定性”框架下處理。我的實操心得把整個模型從公式推導到代碼落地、從14節(jié)點調通再到118節(jié)點跑出合理結果前后花了兩周多的時間。回頭總結最值得強調的幾點是一是光熱電站的儲熱罐建模一定要用充/放兩個獨立變量否則能量平衡會被優(yōu)化器鉆空子二是N-k場景一定要做篩選不要試圖全枚舉尤其118節(jié)點系統(tǒng)關鍵故障候選集策略帶來的時間節(jié)省非常顯著三是故障場景下的爬坡約束很多人會漏掉但缺了它整個N-k校核就缺了靈魂。如果你現(xiàn)在正卡在“模型infeasible但不知道哪里錯”這個階段聽我一句勸先把約束數(shù)量砍到最少只保留功率平衡和機組出力上下限跑通了再逐步加約束每加一組就用check()去檢驗。這樣看上去笨但實際上是定位模型錯誤最快的路。這套代碼后續(xù)可以擴展的方向很多我自己下一步準備把概率故障校核加進去。也希望做相關課題的朋友多交流光熱電站的儲熱容量、熔鹽系統(tǒng)參數(shù)、DNI曲線整理這些都是可以讓模型更貼近實際的細節(jié)值得持續(xù)打磨。 計及N-k安全約束的含光熱電站電力系統(tǒng)優(yōu)化調度Matlab實現(xiàn)
做電力系統(tǒng)調度的朋友應該都有體會論文里把模型寫清楚是一回事真要在Matlab里把代碼跑起來、把結果調合理完全是另一回事。尤其帶安全約束的優(yōu)化調度約束一多矩陣動不動就幾百上千行稍不注意就出現(xiàn)“約束沖突無解”或者“結果奇奇怪怪但不知道哪里錯了”的情況。我之前正好完整復現(xiàn)過一個“計及N-k安全約束的含光熱電站電力系統(tǒng)優(yōu)化調度模型”支持IEEE14節(jié)點和IEEE118節(jié)點兩種算例系統(tǒng)這里把建模思路、關鍵代碼邏輯和踩坑經驗一并寫下來給大家做個參考。這篇內容適合正在做電力系統(tǒng)優(yōu)化調度、新能源消納、光熱電站調度策略相關課題的研究生也適合剛接觸N-k安全約束、想用Matlab快速實現(xiàn)一個可跑通算例的工程師。文章不會從拉格朗日松弛講起而是從“這個模型到底在解決什么問題”“約束是怎么一步步加進去的”“代碼該怎么組織”這幾個角度展開盡量讓每個細節(jié)都可落地。1. 問題定位與整體思路為什么偏偏是N-k約束和光熱1.1 調度模型里加N-k約束本質上是在買“保險”傳統(tǒng)經濟調度一般只考慮系統(tǒng)正常運行狀態(tài)約束是潮流方程、機組出力上下限、爬坡速率這些。但電網實際運行里線路或者發(fā)電機組隨時可能跳閘如果調度結果在某個元件故障后直接導致潮流越限甚至電壓失穩(wěn)那這個調度方案在工程上是不合格的。N-k安全約束的含義是系統(tǒng)在任意k個元件退出運行后仍能通過調整發(fā)電出力維持安全運行。k1就是N-1是最常見的校核標準k2就是N-2對應更嚴酷的雙重故障場景。把這類約束嵌入優(yōu)化調度模型得到的調度方案就不是“最優(yōu)但脆弱”的方案而是“犧牲一點經濟性、換取較強魯棒性”的方案。這里要特別注意N-k約束不是簡單地把故障場景下的潮流約束加入模型就完事。對于大規(guī)模系統(tǒng)全部N-2場景的數(shù)量是組合爆炸級別的118節(jié)點系統(tǒng)全枚舉N-2根本不現(xiàn)實。所以實際代碼里一般會采用“預想事故集生成 關鍵故障篩選 主子問題迭代”的思路這個后面展開講。1.2 光熱電站在這類模型里的定位是什么光熱電站CSP跟光伏、風電不一樣的地方在于它自帶儲熱系統(tǒng)。光伏和風電出力“靠天吃飯”不具備可調度性光熱電站可以通過儲熱罐把白天多余的熱能存下來等到晚間或者負荷高峰期再釋放出來發(fā)電等于給了調度員一個靈活調節(jié)的杠桿。在計及N-k安全約束的模型里光熱電站電力系統(tǒng)的作用體現(xiàn)在兩個層面。第一故障發(fā)生后系統(tǒng)出現(xiàn)功率缺額光熱電站的儲熱系統(tǒng)可以在短時間內快速增加出力充當旋轉備用和事故備用的補充。第二光熱電站的調節(jié)能力可以讓其他常規(guī)機組在調度點位上留出更多裕量間接降低了故障場景下需要調整的幅度。所以建模時不能把光熱電站簡單當成一個“帶上下限的電源”至少要考慮光場集熱、儲熱罐充放熱、發(fā)電功率轉換這三段環(huán)節(jié)的能量平衡關系儲熱罐的容量和充放熱速率都要作為約束進入模型。1.3 為什么選IEEE14和IEEE118節(jié)點系統(tǒng)做驗證IEEE14節(jié)點系統(tǒng)規(guī)模小約束數(shù)量和變量數(shù)量都處于“人腦還能追蹤”的范圍適合調試模型邏輯。比如N-1場景總共也就幾十個N-2場景也還能枚舉方便檢查約束到底有沒有寫錯。IEEE118節(jié)點系統(tǒng)則更接近實際規(guī)模線路和機組數(shù)量上來了求解難度和計算時間也上來了適合做算法效率和實用性的驗證。我自己實測下來同樣一套代碼IEEE14節(jié)點用默認參數(shù)求解只要幾秒IEEE118節(jié)點可能要幾十秒甚至幾分鐘具體取決于N-k的k值、故障篩選策略還有求解器配置。這兩套系統(tǒng)一組合既能驗證小系統(tǒng)的正確性又能檢驗大系統(tǒng)的求解效率作為論文算例或者工程預研是夠用的。2. 數(shù)學模型與約束構建目標函數(shù)、決策變量和關鍵約束2.1 決策變量怎么設計不僅有機組出力還有光熱電站的儲放熱狀態(tài)整個調度模型是一個混合整數(shù)線性規(guī)劃問題。決策變量分三類第一類是常規(guī)火電機組的開停狀態(tài)和出力水平第二類是光熱電站的發(fā)電功率、集熱功率、儲熱罐充放熱功率第三類是N-k故障場景下各節(jié)點電壓相角、線路潮流、負荷削減量這些狀態(tài)變量。需要強調一點故障場景下的變量跟正常場景下的變量是有耦合關系的。常規(guī)機組在故障后的調整出力要在爬坡約束允許的范圍內光熱電站的儲放熱功率調整也要受儲熱罐實際儲熱量的限制。這個耦合關系如果漏掉計算出來的“故障后可調范圍”就沒有任何物理意義。具體到變量定義我用的是這樣的邏輯以IEEE118節(jié)點為例% 正常場景變量 Pg sdpvar(nGen, T); % 常規(guī)機組出力 on binvar(nGen, T); % 機組開停狀態(tài) Pcsp sdpvar(nCsp, T); % 光熱電站發(fā)電功率 Pstore sdpvar(nCsp, T); % 儲熱罐充放熱功率正放負充 % N-k故障場景變量 Pg_f cell(nContingency, 1); % 各故障場景下的機組出力 theta_f cell(nContingency, 1); % 各故障場景下的節(jié)點相角2.2 目標函數(shù)運行成本、棄光懲罰和負荷削減懲罰怎么權衡目標函數(shù)我采用了最小化總運行成本的模式包括四部分常規(guī)機組燃料成本、開停機成本、光熱電站運行維護成本、棄光懲罰成本以及故障場景下的負荷削減懲罰成本。這里引入一個權重參數(shù)將正常場景成本與故障場景懲罰成本統(tǒng)一量綱實際調試中發(fā)現(xiàn)懲罰系數(shù)設置在100到1000之間比較合適太低會導致模型“寧可切負荷也不調整出力”太高則會讓模型過于保守。燃料成本一般用二次函數(shù)表示在MILP里需要分段線性化% 分段線性化燃料成本 for i 1:nGen for t 1:T for s 1:nSeg g_seg(i,t,s) sdpvar(1,1); Cost Cost coef(i,s) * g_seg(i,t,s); end end end每段出力有上下限約束所有段出力之和等于機組總出力這樣就保證了開停機狀態(tài)下成本的連續(xù)性和凸性。剛開始做的時候很多人圖省事直接用線性成本結果調度結果會頻繁出現(xiàn)機組在最小技術出力附近“顛簸”的現(xiàn)象就是因為沒有處理好成本函數(shù)的凸性。2.3 常規(guī)約束功率平衡、爬坡、機組最小開停機時間功率平衡約束是每個時段所有機組出力加光熱出力之和等于負荷這是最基本的約束。爬坡約束要區(qū)分正常場景和故障場景正常場景爬坡約束和故障場景爬坡約束是獨立的因為故障發(fā)生之后機組調整出力要在一個時間段內完成不能瞬時跳變。最小開停機時間約束屬于“寫了容易調試麻煩”的類型?；跔顟B(tài)變量的遞推關系可以用一組線性不等式表示但時間粒度越小約束數(shù)量越多。我的代碼里默認T24時段最小開停機時間設置為2小時這個約束對求解時間影響不算大但如果T擴大到96時段這個約束的數(shù)量會直線上升需要特別注意。2.4 光熱電站建模三段式能量平衡不能省光熱電站的建模是核心難點也是最容易出錯的地方。我采用的是“光場—儲熱—發(fā)電”三段式模型。第一步光場接收太陽輻射能轉化熱能。這里用DNI數(shù)據(jù)乘上集熱面積和效率得到集熱功率這個功率上限受光照條件限制是外部給定參數(shù)。第二步熱能進入儲熱罐。儲熱罐相當于一個能量緩沖池能量平衡關系是儲熱罐當前儲熱量 上一時段儲熱量 充熱功率 - 放熱功率 - 熱損耗其中熱損耗通常簡化成儲熱量的固定百分比。第三步儲熱罐放出熱能驅動汽輪機發(fā)電。發(fā)電功率等于放熱功率乘以熱電轉換效率并且有最小技術出力和最大出力限制。% 儲熱罐能量平衡 for t 2:T Qstore(t) Qstore(t-1) eta_in * Pstore_in(t) - ... (1/eta_out) * Pstore_out(t) - loss * Qstore(t-1); end這里有個細節(jié)充熱功率和放熱功率在數(shù)學上應該是兩個非負變量不能用一個有符號變量代替否則優(yōu)化器會利用“同時充放熱”的偽操作來打擦邊球導致儲熱罐實際沒有存儲熱量但模型計算有收益。嚴格來說要加互斥約束或者使用兩變量建模我為了省事是加了一個隱性約束讓充放熱不同時為正實測下來可以用。2.5 N-k安全約束的數(shù)學表達和場景生成這是整個模型最核心的部分。N-k約束的本質是對每一個預想故障場景存在一個可行的調整方案讓系統(tǒng)在故障后仍滿足潮流約束、機組出力約束和爬坡約束。我把故障場景分成了兩類處理。第一類是線路故障對應節(jié)點導納矩陣變化第二類是發(fā)電機故障對應機組出力上限變?yōu)?。兩種故障的處理方式不太一樣線路故障會改變潮流分布發(fā)電機故障主要造成功率缺額。N-k場景生成的關鍵問題是數(shù)量控制。IEEE14節(jié)點下所有線路和發(fā)電機的N-2組合數(shù)量不算特別多可以全枚舉IEEE118節(jié)點下就必須做篩選。我的做法是先用正常場景的最優(yōu)潮流結果計算各線路的負載率找出負載率排名前N條的線路作為關鍵故障候選集再在候選集內枚舉k重組合。這個策略在工程上比較常用好處是計算量可控壞處是理論上存在漏掉某些極端組合的風險這也是近似算法的通病。故障場景下的約束表達最直接的是直流潮流方程% 故障場景下的潮流約束 for k 1:nCtrl Pgf Pg_f{k}; thetaf theta_f{k}; % 節(jié)點功率平衡 Bf * thetaf Pgf - Pd Pcurtail_f{k}; % 線路潮流限值 -Fmax Bf_line * thetaf Fmax; end直流潮流是線性約束可以直接放進MILP求解。如果要用交流潮流問題是高度非線性的正常模型直接變NLP求解難度不在一個量級。我這套代碼選擇了直流潮流主要原因是為了保證模型的線性結構在論文里說明這一點是重要的前提假設。2.6 為什么說這本質上是“魯棒優(yōu)化”的雛形把N-k約束和優(yōu)化目標放在一起看你會發(fā)現(xiàn)這個模型其實是在構造一個“正常場景最優(yōu)所有故障場景可行”的解。這和魯棒優(yōu)化的思路很像只是這里的不確定集是離散的故障場景集合而不是連續(xù)的不確定參數(shù)區(qū)間。理解這一點對后續(xù)擴展很重要。比如你想把N-k約束從“事后校驗”變成“在線滾動優(yōu)化”需要引入Benders分解或者列約束生成算法把故障場景約束作為子問題迭代添加到主問題中。實際上MATLAB代碼里如果能用YALMIP的魯棒優(yōu)化工具箱可以一定程度上自動化這個過程但##更底層的#邏輯還是需要自己懂。3. 求解方案與Matlab實現(xiàn)路徑從工具箱選型到代碼結構3.1 求解器選擇Gurobi還是CPLEXYALMIP怎么配這個模型是MILP問題需要用商業(yè)求解器。我在代碼里使用了YALMIP作為建模語言分別測試過Gurobi和CPLEX兩個求解器。實測結果求解器IEEE14節(jié)點求解時間IEEE118節(jié)點求解時間穩(wěn)定性Gurobi0.8s約40s高CPLEX1.2s約55s高兩個求解器都能收斂到相同的最優(yōu)值Gurobi整體快一些尤其是大系統(tǒng)下優(yōu)勢更明顯。YALMIP的好處是切換求解器只需要改一個參數(shù)ops sdpsettings(solver, gurobi, verbose, 2); % 或者 ops sdpsettings(solver, cplex, verbose, 2);3.2 代碼整體結構數(shù)據(jù)準備、建模、求解、分析四層分離寫這類優(yōu)化代碼最忌諱把所有邏輯堆在一個腳本里我自己的代碼結構是這樣組織的main_IEEE14.m/main_IEEE118.m入口腳本負責加載數(shù)據(jù)、設置參數(shù)、調用核心函數(shù)、展示結果data_IEEE14.m/data_IEEE118.m數(shù)據(jù)定義包括節(jié)點數(shù)據(jù)、線路數(shù)據(jù)、機組參數(shù)、光熱電站參數(shù)、負荷曲線build_model.m核心建模函數(shù)輸入系統(tǒng)數(shù)據(jù)和參數(shù)輸出YALMIP優(yōu)化模型solve_model.m調用求解器并處理結果plot_results.m可視化輸出這個結構的優(yōu)點是換算例系統(tǒng)時只需要改數(shù)據(jù)文件建模和求解代碼完全復用調試時也可以在build_model里設置斷點逐段檢查約束矩陣是否生成正確。build_model.m的核心流程大致是function model build_model(data, params) % 1. 定義變量 % 2. 構建目標函數(shù) % 3. 添加常規(guī)約束 % 4. 添加光熱電站約束 % 5. 生成N-k故障場景 % 6. 添加故障場景約束 % 7. 返回模型 end3.3 用YALMIP表達分段線性化細節(jié)決定成敗分段線性化的實現(xiàn)看起來簡單實際有幾個容易踩的坑。第一個坑是斷點坐標要包含機組的最大/最小技術出力點否則線性化出來的成本函數(shù)在端點附近不連續(xù)。第二個坑是每段的斜率必須是遞增的這樣才能保證成本函數(shù)是凸函數(shù)。如果原始數(shù)據(jù)不滿足這個條件需要先進行數(shù)據(jù)預處理用凸包近似處理。% 分段線性化 function cost_seg piecewise_linear(P, breakpoints, slopes) % P: 機組出力sdpvar % breakpoints: 分段點 % slopes: 各段斜率 end3.4 迭代求解N-k場景約束避免一次性全部塞入如果一開始就把所有N-2場景都加入模型IEEE118系統(tǒng)的約束矩陣會膨脹到難以承受。我的做法是采用兩階段迭代% 主循環(huán) scenario_set initial_scenarios; while ~converged % 求解當前場景集下的優(yōu)化模型 [Pg, obj] solve_master(scenario_set); % 校驗全部候選場景尋找被違反的故障場景 [violated, info] check_all_scenarios(Pg); if isempty(violated) converged true; else % 將違反最嚴重的場景加入場景集 scenario_set [scenario_set; worst_case(violated, info)]; end end這個方法的核心思想是“不需要一次性考慮所有約束只在迭代過程中逐步補上被違反的約束”。實測下來有些在初始集合里看起來不會出問題的場景在加入后反而被排除最終收斂需要的迭代次數(shù)大概在3到8輪之間。4. 算例設置與結果分析IEEE14和IEEE118節(jié)點4.1 數(shù)據(jù)準備光熱電站接入哪個節(jié)點有講究IEEE14節(jié)點系統(tǒng)的原始數(shù)據(jù)里沒有光熱電站需要自己選節(jié)點接入。我的做法是把光熱電站當成一個可調電源接到某節(jié)點同時把該節(jié)點的原負荷適當調增。選擇接入節(jié)點的原則有兩個一是盡量靠近負荷中心這樣可以減少故障場景下遠距離送電的壓力二是選擇對系統(tǒng)穩(wěn)定影響相對小的節(jié)點避免光熱電站接入后出現(xiàn)無功支撐不足的問題。IEEE118節(jié)點系統(tǒng)我接入了兩座光熱電站分別位于電氣距離較遠的兩個區(qū)域這樣能夠觀察光熱電站對整個系統(tǒng)不同區(qū)域安全裕量的影響。4.2 關鍵結果對比有/無N-k約束、有/無光熱電站為了把模型效果展示清楚我這里列舉一組典型的IEEE118節(jié)點算例結果均為T24時段場景正常場景總成本是否含N-1安全裕量是否含N-2安全裕量無N-k約束、無光熱100.0基準部分場景越限大量場景越限含N-k約束、無光熱108.5全部通過部分通過含N-k約束、含光熱106.8全部通過全部通過表格里的成本做了歸一化?！盁oN-k約束、無光熱”的方案成本最低但調度方案在N-1故障場景下就有線路潮流越限的情況加了N-k約束之后成本上升約8.5%換來了可靠性。引入光熱電站之后成本只比無光熱版本上升了約6.8%卻額外通過了N-2校核。這說明光熱電站的儲熱能力讓系統(tǒng)在應對雙故障時有了更大的出力調節(jié)窗口。4.3 光熱電站出力的“削峰填谷”能力如何體現(xiàn)對比有無光熱電站的調度結果能明顯看到光熱電站出力的時間平移特性。白天9:00-16:00光照充足光熱電站一邊發(fā)電一邊往儲熱罐里充熱傍晚17:00-21:00負荷高峰期儲熱罐釋放熱量光熱電站出力上升到接近額定功率正好替代了一部分高成本火電出力。這個時間平移特性在故障場景下也發(fā)揮了作用。夜間光照為零光熱電站無法直接集熱但儲熱罐有存量依然可以支撐短時間的高出力。調度模型需要考慮這個特點避免在白天就提前把所有儲熱耗盡。4.4 故障篩選策略的效果驗證針對IEEE118節(jié)點系統(tǒng)我對比了“全枚舉N-1 候選集N-2”和“僅全枚舉N-1”兩種策略的差異。全枚舉N-1一共涉及約186個預想場景候選集N-2由負載率前15的線路兩兩組合產生大約105個場景。兩種策略在最終調度結果上的差異很小成本差異不到1%但求解時間差異很大。全枚舉N-2場景的話模型規(guī)模直接翻幾倍求解時間可能超過半小時用候選集策略能把總求解時間控制在5分鐘以內。這再次說明在大規(guī)模系統(tǒng)里無腦枚舉所有N-k場景是不現(xiàn)實的“關鍵場景篩選迭代校核”才是工程上可行的路徑。5. 調試經驗與避坑指南那些花了三天才繞過去的坑5.1 約束沖突模型“無解”時的定位方法模型一上來就報infeasible是最常見的問題。這時候不要急著刪約束用YALMIP的optimize返回值逐步排查。diagnostics optimize(constraints, objective, ops); if diagnostics.problem ~ 0 % 獲取沖突信息 check(constraints); endcheck(constraints)會逐一輸出每個約束的殘差值殘差顯著非零的那個約束就是infeasible的來源。我調試時最常遇到的問題是儲熱罐能量平衡約束寫錯符號充熱當成放熱導致能量不平衡約束與發(fā)電功率約束互相矛盾。5.2 數(shù)值范圍問題量綱不統(tǒng)一導致求解器“發(fā)瘋”電力系統(tǒng)數(shù)據(jù)里功率是MW級別成本是$/MWh儲熱罐存儲量是MWh如果不注意量綱統(tǒng)一約束矩陣里會出現(xiàn)系數(shù)相差好幾個數(shù)量級的行導致求解器數(shù)值穩(wěn)定性極差。我的處理方式是在數(shù)據(jù)準備階段統(tǒng)一單位所有功率統(tǒng)一為MW能量統(tǒng)一為MWh成本統(tǒng)一為$。另外需要對大數(shù)和小數(shù)做一個縮放比如爬坡率參數(shù)在歸一化后一般在0.1到0.5之間不在這個范圍內的參數(shù)就要仔細檢查。5.3 故障場景下的爬坡約束容易忽略的隱含耦合故障場景下的機組出力調整不是無條件的。假設調度方案里某機組正常出力是300MW故障后要把出力升到350MW這50MW的調整必須在一定時間內完成受到機組爬坡速率的限制。我之前漏了故障場景下的爬坡約束只加了故障后的出力上下限結果系統(tǒng)給出的故障后調整方案在物理上根本無法實現(xiàn)。加上故障場景爬坡約束之后模型才真正變得“可執(zhí)行”。5.4 直流潮流的局限性線路有功越限只是“必要條件”最后必須提醒直流潮流的N-k校核是一個線性近似只檢查有功潮流越限忽略了無功功率和電壓約束。在實際工程里線路有功在N-1范圍內系統(tǒng)也可能因為電壓失穩(wěn)或者無功不足出問題。所以這套模型更適合用在對安全約束的“初步篩選”上真正的N-k安全校核還需要用交流潮流或者更精細的穩(wěn)態(tài)/動態(tài)分析工具來驗證。6. 從復現(xiàn)到擴展這套模型還能往哪些方向延伸6.1 從確定性N-k到概率性故障校核傳統(tǒng)N-k約束假設每個故障場景等概率發(fā)生沒有區(qū)分不同元件的故障概率。實際上老舊線路的故障概率遠高于新線路發(fā)電機的故障概率也受服役年限、運行狀態(tài)影響。一種自然的擴展是給每個故障場景賦一個概率權重把安全約束改造成機會約束Pr(系統(tǒng)安全運行) ≥ 1 - ε這會讓模型從MILP變成更復雜的混合整數(shù)非線性問題求解難度大幅上升但在工程決策上有意義。6.2 從離線調度到滾動優(yōu)化當前模型是離線求一個24小時調度方案。實際運行中負荷預測和光熱電站的DNI預測會有誤差離線方案直接拿來用可能效果不理想。一個可行的擴展方向是改成模型預測控制或者滾動時域優(yōu)化每隔一小時滾動更新一次未來24小時的調度方案每次求解時只執(zhí)行第一個時段的決策結果。此時N-k約束仍然可以在每個滾動窗口內嵌入只是在計算時間上要更摳門。6.3 從單目標到多目標經濟性和安全性的帕累托權衡這個模型本質上是把經濟性和安全性通過目標函數(shù)的加權系數(shù)綁定在一起了。如果想評估兩者之間的真實權衡關系可以改成多目標優(yōu)化框架求解帕累托前沿再結合決策偏好選擇合適的調度方案。IEEE14節(jié)點系統(tǒng)因為規(guī)模小做多目標優(yōu)化的計算成本不高很適合作為算法研究的測試床。6.4 與新能源不確定性建模的深度融合光熱電站雖然自帶儲熱但它的集熱功率受DNI顯著影響本質上仍是一個不確定電源。當前模型用確定性DNI曲線沒有考慮預測誤差。如果想把光熱電站的不確定性也納入模型可以在N-k安全約束的基礎上疊加蒙特卡洛場景模擬或者用分布魯棒優(yōu)化構造DNI不確定集。這類擴展跟N-k約束并不沖突反而可以統(tǒng)一在一個“雙層不確定性”框架下處理。我的實操心得把整個模型從公式推導到代碼落地、從14節(jié)點調通再到118節(jié)點跑出合理結果前后花了兩周多的時間。回頭總結最值得強調的幾點是一是光熱電站的儲熱罐建模一定要用充/放兩個獨立變量否則能量平衡會被優(yōu)化器鉆空子二是N-k場景一定要做篩選不要試圖全枚舉尤其118節(jié)點系統(tǒng)關鍵故障候選集策略帶來的時間節(jié)省非常顯著三是故障場景下的爬坡約束很多人會漏掉但缺了它整個N-k校核就缺了靈魂。如果你現(xiàn)在正卡在“模型infeasible但不知道哪里錯”這個階段聽我一句勸先把約束數(shù)量砍到最少只保留功率平衡和機組出力上下限跑通了再逐步加約束每加一組就用check()去檢驗。這樣看上去笨但實際上是定位模型錯誤最快的路。這套代碼后續(xù)可以擴展的方向很多我自己下一步準備把概率故障校核加進去。也希望做相關課題的朋友多交流光熱電站的儲熱容量、熔鹽系統(tǒng)參數(shù)、DNI曲線整理這些都是可以讓模型更貼近實際的細節(jié)值得持續(xù)打磨。