現(xiàn)電熱綜合能源系統(tǒng)多場(chǎng)景分布魯棒優(yōu)化調(diào)度方法)
先聊一個(gè)真實(shí)調(diào)度場(chǎng)景早上8點(diǎn)預(yù)測(cè)下午3點(diǎn)風(fēng)大你按這個(gè)預(yù)測(cè)排好了機(jī)組出力結(jié)果到了2點(diǎn)半風(fēng)突然變小電鍋爐還在滿功率吃電儲(chǔ)熱罐可用熱量又快見(jiàn)底這時(shí)候再調(diào)整機(jī)組爬坡已經(jīng)來(lái)不及了。做電熱綜合能源調(diào)度的人多少都經(jīng)歷過(guò)這種“預(yù)測(cè)一時(shí)爽調(diào)度兩行淚”的時(shí)刻。純靠點(diǎn)估計(jì)做決策遲早被不確定性教做人但完全按最壞場(chǎng)景規(guī)劃成本高得離譜電廠和熱網(wǎng)公司都接受不了。于是這些年“數(shù)據(jù)驅(qū)動(dòng) 分布魯棒優(yōu)化”這套思路越來(lái)越熱核心就是利用歷史數(shù)據(jù)構(gòu)造場(chǎng)景再在場(chǎng)景周?chē)粢粋€(gè)“模糊集”來(lái)兜底最終在“經(jīng)濟(jì)性”和“魯棒性”之間取一個(gè)可調(diào)的平衡點(diǎn)。這篇文章圍繞多離散場(chǎng)景分布魯棒Distributionally Robust Optimization, DRO方法聊聊怎么在電熱綜合能源系統(tǒng)里落地一套完整的Matlab調(diào)度程序。我會(huì)從問(wèn)題建模、不確定性處理、模糊集構(gòu)造、對(duì)偶轉(zhuǎn)化、代碼框架到踩坑實(shí)錄全部過(guò)一遍適合正在做綜合能源調(diào)度、園區(qū)能源管理相關(guān)課題的研究生以及想從確定性優(yōu)化轉(zhuǎn)向魯棒優(yōu)化的工程師參考。1. 先把問(wèn)題說(shuō)透電熱綜合能源系統(tǒng)為什么難優(yōu)化1.1 電和熱兩種能量耦合不是兩套獨(dú)立系統(tǒng)電熱綜合能源系統(tǒng)之所以讓人頭疼是因?yàn)殡?、熱本身就是?qiáng)耦合的。典型園區(qū)系統(tǒng)里主要有這么幾類設(shè)備熱電聯(lián)產(chǎn)機(jī)組CHP這是電-熱耦合的核心。抽汽式CHP的發(fā)電量和供熱量必須在可行域內(nèi)聯(lián)動(dòng)發(fā)多了電往往熱也會(huì)多靈活性受限。電鍋爐EB把電轉(zhuǎn)化成熱相當(dāng)于給系統(tǒng)增加了一條“電轉(zhuǎn)熱”的靈活路徑但它本質(zhì)上是負(fù)荷會(huì)增加電網(wǎng)購(gòu)電壓力。儲(chǔ)熱罐TES儲(chǔ)熱罐是解耦電-熱“強(qiáng)綁定”的關(guān)鍵緩沖裝置熱負(fù)荷高峰時(shí)放熱熱負(fù)荷低谷時(shí)充熱相當(dāng)于給系統(tǒng)提供一個(gè)時(shí)間維度上的平移能力。常規(guī)火電/上級(jí)電網(wǎng)購(gòu)電提供電功率平衡的兜底手段。風(fēng)電/光伏這里是不確定性的主要來(lái)源也直接決定了為什么要做分布魯棒。電網(wǎng)側(cè)需要滿足CHP發(fā)電 風(fēng)電實(shí)際出力 購(gòu)電 電負(fù)荷 電鍋爐耗電。熱網(wǎng)側(cè)需要滿足CHP供熱量 電鍋爐供熱量 儲(chǔ)熱罐放熱 熱負(fù)荷 儲(chǔ)熱罐充熱 熱網(wǎng)損失。問(wèn)題的難點(diǎn)不在方程多而在于CHP 的發(fā)電和供熱捆綁在一起風(fēng)電又隨機(jī)波動(dòng)。一旦風(fēng)電預(yù)測(cè)不準(zhǔn)系統(tǒng)只能通過(guò)快速調(diào)整CHP或者改變電鍋爐出力來(lái)平衡但這兩者又牽動(dòng)著熱力平衡牽一發(fā)動(dòng)全身。1.2 三種不確定性處理思路為什么我選分布魯棒對(duì)風(fēng)電出力和負(fù)荷預(yù)測(cè)誤差的處理現(xiàn)在主流有三條路線它們之間存在本質(zhì)區(qū)別方法不確定集合優(yōu)化目標(biāo)對(duì)數(shù)據(jù)的依賴保守程度確定性優(yōu)化無(wú)直接使用預(yù)測(cè)值單場(chǎng)景經(jīng)濟(jì)最優(yōu)只要預(yù)測(cè)值最低但容易失衡隨機(jī)規(guī)劃(SAA)若干離散場(chǎng)景及概率場(chǎng)景期望成本最小需要大量場(chǎng)景中低依賴場(chǎng)景質(zhì)量傳統(tǒng)魯棒優(yōu)化盒式/橢球不確定域最壞情況下成本最小只需不確定變量的上下界高過(guò)于保守分布魯棒優(yōu)化(DRO)以歷史分布為中心的模糊集最壞概率分布下的期望成本最小需要?dú)v史數(shù)據(jù)構(gòu)造場(chǎng)景和模糊集可調(diào)的中間水平隨機(jī)規(guī)劃的問(wèn)題是你給的場(chǎng)景概率分布是固定的但實(shí)際風(fēng)電誤差分布可能跟假設(shè)偏差很大傳統(tǒng)魯棒的問(wèn)題是它只關(guān)心最壞情況凡是落在盒式集合內(nèi)的可能值都按最壞結(jié)果兜底最終調(diào)度成本可能比實(shí)際最優(yōu)貴20%以上。分布魯棒走的是中間路線——只假設(shè)真實(shí)分布落在經(jīng)驗(yàn)分布周?chē)囊粋€(gè)“模糊集”內(nèi)然后優(yōu)化最壞分布下的期望成本。當(dāng)歷史數(shù)據(jù)足夠多時(shí)模糊集半徑可以收得很小結(jié)果趨近隨機(jī)規(guī)劃當(dāng)數(shù)據(jù)不充分或者你對(duì)分布沒(méi)信心時(shí)可以調(diào)大半徑結(jié)果向傳統(tǒng)魯棒靠攏。這個(gè)“可調(diào)”特性在工程里非常實(shí)用因?yàn)樗屇隳苡靡惶姿惴蚣苋ミm配不同數(shù)據(jù)質(zhì)量的項(xiàng)目。1.3 多離散場(chǎng)景這個(gè)“多”字到底指什么標(biāo)題里的“多離散場(chǎng)景”可以拆成兩個(gè)層面理解。第一層是用有限離散場(chǎng)景近似經(jīng)驗(yàn)分布把歷史風(fēng)速、歷史負(fù)荷誤差數(shù)據(jù)通過(guò)聚類或者采樣歸納成幾十個(gè)代表性場(chǎng)景每個(gè)場(chǎng)景對(duì)應(yīng)一組24小時(shí)的風(fēng)電出力曲線。這樣隨機(jī)變量就從“連續(xù)分布”退化為“帶概率權(quán)重的離散點(diǎn)集”計(jì)算規(guī)?？煽?。第二層是模糊集也按離散場(chǎng)景構(gòu)造分布魯棒優(yōu)化不需要顯式描述連續(xù)分布而是以“這些離散場(chǎng)景組成的經(jīng)驗(yàn)分布”為中心向外擴(kuò)一個(gè)半徑為alpha的距離球通常是Wasserstein距離。只要真實(shí)分布沒(méi)有超出這個(gè)球最壞情況下的期望成本就在模型掌控范圍內(nèi)。這種“離散場(chǎng)景 模糊球”的組合既保留了隨機(jī)規(guī)劃對(duì)場(chǎng)景信息的利用又繼承了魯棒優(yōu)化對(duì)不確定性的兜底能力是這個(gè)方向近幾年成為熱點(diǎn)的關(guān)鍵。2. 分布魯棒優(yōu)化的建模與轉(zhuǎn)化怎么把min-max問(wèn)題變成可求解問(wèn)題2.1 先建立確定性調(diào)度主模型整套優(yōu)化模型可以表達(dá)為目標(biāo)最小化系統(tǒng)總運(yùn)行成本包括CHP燃料成本、向上級(jí)電網(wǎng)購(gòu)電成本、棄風(fēng)懲罰、切負(fù)荷懲罰。約束電功率平衡、熱功率平衡、CHP可行域約束、電鍋爐出力上下限、儲(chǔ)熱罐SOC連續(xù)性約束與容量約束、購(gòu)電上下限、爬坡約束等。其中CHP運(yùn)行可行域是關(guān)鍵約束通常用一組線性不等式圍成的多邊形描述即滿足P_chp_min P_chp_t P_chp_maxH_chp_min H_chp_t H_chp_maxP_chp_t c1 * H_chp_t c2c3 * H_chp_t - P_chp_t c4這四個(gè)不等式基本可以刻畫(huà)抽汽式CHP的可行域。如果你的CHP是背壓式可以直接簡(jiǎn)化為P c * H那電熱耦合更緊調(diào)度難度更大。儲(chǔ)熱罐模型相對(duì)簡(jiǎn)單但要注意不能同時(shí)充放熱的約束S_t S_{t-1} eta_c * H_char_t - H_dis_t / eta_d - eta_loss * S_{t-1}0 H_char_t H_char_max * u_c_t0 H_dis_t H_dis_max * u_d_tu_c_t u_d_t 1這里的u_c_t和u_d_t是0-1變量一旦加進(jìn)去模型就變成了MILP混合整數(shù)線性規(guī)劃。如果系統(tǒng)規(guī)模不大直接用YALMIP Gurobi求解沒(méi)問(wèn)題如果規(guī)模大要考慮用啟發(fā)式或者滾動(dòng)時(shí)域控制去壓規(guī)模。2.2 把不確定性裝進(jìn)去從SAA到DRO引入風(fēng)電不確定性后簡(jiǎn)化的目標(biāo)函數(shù)形式變成一個(gè)min-max雙層結(jié)構(gòu)外層min是調(diào)度決策變量x機(jī)組出力、儲(chǔ)熱罐充放、購(gòu)電等內(nèi)層max是針對(duì)模糊集D中所有可能概率分布P最大化期望成本。DRO目標(biāo)含義min_x max_{P in D} E_P[f(x, xi)]在“最壞的合理分布”下期望成本最小這里xi就是隨機(jī)變量比如風(fēng)電誤差向量、負(fù)荷誤差向量D就是模糊集。模糊集最常見(jiàn)的構(gòu)造方式是Wasserstein球D { P | W(P, P_hat) alpha }其中P_hat是歷史數(shù)據(jù)構(gòu)造的經(jīng)驗(yàn)分布alpha是模糊集半徑。Wasserstein距離衡量?jī)蓚€(gè)分布之間的“搬運(yùn)成本”直觀理解就是要把經(jīng)驗(yàn)分布P_hat變成真實(shí)分布P最少需要“搬運(yùn)”多少概率質(zhì)量。alpha越大真實(shí)分布離經(jīng)驗(yàn)分布可以越遠(yuǎn)模型越保守。工程實(shí)現(xiàn)中我用的是1-范數(shù)Wasserstein距離因?yàn)閷?duì)偶轉(zhuǎn)化后能保持線性約束結(jié)構(gòu)Gurobi和CPLEX可以直接吃下。如果上2-范數(shù)模型會(huì)帶二次約束求解器壓力明顯變大除非你確有必要我不太推薦。2.3 對(duì)偶轉(zhuǎn)化的三板斧直接求解min-max問(wèn)題是不現(xiàn)實(shí)的實(shí)際做法是把它對(duì)偶成一個(gè)單層的min問(wèn)題。以Wasserstein模糊集為例核心思路分三步第一步內(nèi)層max對(duì)偶變換。在凸目標(biāo)函數(shù)和凸模糊集條件下內(nèi)層最大化問(wèn)題可以等價(jià)地轉(zhuǎn)化為一個(gè)關(guān)于對(duì)偶變量lambda的最小化問(wèn)題。這一步讓“最壞分布”消失換來(lái)的是對(duì)每個(gè)離散場(chǎng)景引入一組額外約束和輔助變量s_i。第二步引入場(chǎng)景級(jí)最優(yōu)值函數(shù)。每個(gè)離散場(chǎng)景xi_i都會(huì)產(chǎn)生一個(gè)“最壞成本”的支撐函數(shù)通過(guò)引入輔助變量s_i把對(duì)每個(gè)場(chǎng)景的max項(xiàng)線性化。第三步最終得到形如min_{x, lambda0, s_i} lambda * alpha (1/N) * sum(s_i)約束lambda 0對(duì)每個(gè)場(chǎng)景is_i f(x, xi_i) - lambda * d(xi_i, xi_j) 的某種線性化形式這里d是場(chǎng)景之間的Wasserstein距離項(xiàng)f是當(dāng)前決策x在場(chǎng)景xi_i下的運(yùn)行成本。這個(gè)公式寫(xiě)出來(lái)可能有點(diǎn)抽象但在代碼里的操作其實(shí)很簡(jiǎn)單不要自己手推全部對(duì)偶約束而是用小規(guī)模測(cè)試3個(gè)場(chǎng)景、3個(gè)時(shí)段去驗(yàn)證對(duì)偶轉(zhuǎn)換后的目標(biāo)值和暴力枚舉max結(jié)果一致。我一開(kāi)始直接套大模型結(jié)果怎么都不對(duì)最后縮小規(guī)模一步步檢查才定位到是某個(gè)極端場(chǎng)景下的約束沒(méi)寫(xiě)全。3. Matlab代碼實(shí)現(xiàn)從零搭一套分布魯棒調(diào)度程序3.1 代碼框架與文件規(guī)劃整個(gè)項(xiàng)目我分成五個(gè)文件維護(hù)結(jié)構(gòu)清晰也方便復(fù)現(xiàn)項(xiàng)目目錄/ ├── main.m % 主程序入口 ├── set_parameters.m % 設(shè)備參數(shù)定義 ├── generate_scenarios.m % 歷史數(shù)據(jù)讀取 場(chǎng)景生成/聚類 ├── build_dro_model.m % YALMIP建模仿真 ├── plot_results.m % 結(jié)果可視化main.m 的結(jié)構(gòu)就是典型的“參數(shù)-場(chǎng)景-建模-求解-畫(huà)圖”五段式%% main.m clc; clear; close all; % 1. 參數(shù)設(shè)置 param set_parameters(); % 2. 讀歷史風(fēng)速數(shù)據(jù)生成離散場(chǎng)景 [scen, prob] generate_scenarios(param); % 3. 構(gòu)建DRO模型并求解 [result, model] build_dro_model(param, scen, prob); % 4. 畫(huà)圖 plot_results(result, param);這里有一點(diǎn)值得強(qiáng)調(diào)場(chǎng)景生成和建模是解耦的。如果你后續(xù)想換數(shù)據(jù)、換聚類方法只需要?jiǎng)觛enerate_scenarios這一個(gè)文件不需要碰主模型這個(gè)設(shè)計(jì)能給你后續(xù)調(diào)參省很多事。3.2 場(chǎng)景生成K-means聚類 拉丁超立方采樣場(chǎng)景生成模塊的核心目標(biāo)是把大量歷史風(fēng)電出力數(shù)據(jù)壓縮成幾十個(gè)代表性離散場(chǎng)景。我用的是兩步法。第一步拉丁超立方采樣(LHS)生成海量預(yù)測(cè)誤差樣本。為什么要用LHS而不是直接蒙特卡洛抽樣因?yàn)長(zhǎng)HS能保證樣本點(diǎn)在概率空間內(nèi)覆蓋更均勻同樣的樣本量下LHS構(gòu)造的經(jīng)驗(yàn)分布更穩(wěn)定模糊集半徑alpha可以選得更小。第二步用K-means對(duì)樣本聚類取聚類中心作為代表場(chǎng)景按每個(gè)簇的樣本比例分配概率。聚類數(shù)N我一般取20~50之間。少于20個(gè)場(chǎng)景分布信息損失太嚴(yán)重模糊集半徑會(huì)被迫加大成本偏向保守多于50個(gè)場(chǎng)景模型規(guī)模膨脹明顯MILP求解時(shí)間會(huì)從幾分鐘漲到一個(gè)小時(shí)。核心代碼示意% generate_scenarios.m 片段 % 歷史誤差樣本: err_hist (N_hist x T) err_hist load_historical_data(); % LHS生成候選樣本 candidate lhsdesign(N_samp, T); % 將[0,1]區(qū)間的LHS樣本映射到經(jīng)驗(yàn)誤差分布 err_samp quantile(err_hist, candidate); % K-means聚類 [idx, C] kmeans(err_samp, N_scen, Replicates, 10); prob histcounts(idx, N_scen) / N_samp; % C就是N_scen x T的離散場(chǎng)景矩陣 scen C;一個(gè)容易忽略的坑聚類前要對(duì)數(shù)據(jù)進(jìn)行歸一化尤其是風(fēng)電和負(fù)荷量綱不同的時(shí)候。有的風(fēng)電場(chǎng)裝機(jī)容量500MW負(fù)荷可能才50MW如果不歸一化聚類結(jié)果會(huì)完全被風(fēng)電數(shù)據(jù)主導(dǎo)負(fù)荷誤差場(chǎng)景根本體現(xiàn)不出來(lái)。我通常對(duì)每個(gè)隨機(jī)變量單獨(dú)歸一化到[0,1]區(qū)間構(gòu)造模糊集后再反歸一化回去。3.3 核心模型YALMIP寫(xiě)DRO的關(guān)鍵代碼模型部分用YALMIP建模求解器接口用Gurobi或CPLEX。整個(gè)DRO模型在YALMIP里非常直白因?yàn)榉植际紧敯艮D(zhuǎn)化后的模型本質(zhì)上是一個(gè)帶額外輔助變量的MILP。決策變量定義片段% build_dro_model.m 片段 T param.T; P_chp sdpvar(T, 1); % CHP電出力 H_chp sdpvar(T, 1); % CHP熱出力 P_eb sdpvar(T, 1); % 電鍋爐耗電 H_eb sdpvar(T, 1); % 電鍋爐供熱 P_buy sdpvar(T, 1); % 購(gòu)電 SOC sdpvar(T1, 1); % 儲(chǔ)熱罐狀態(tài) H_char sdpvar(T, 1); % 充熱 H_dis sdpvar(T, 1); % 放熱 % 輔助變量DRO對(duì)偶變量和場(chǎng)景上界 lambda sdpvar(1); s_var sdpvar(param.N_scen, 1);目標(biāo)函數(shù)的YALMIP寫(xiě)法% 確定性成本部分 obj_det sum(param.c_fuel .* P_chp) sum(param.c_buy .* P_buy) ... sum(param.c_wind * (scen_forecast - P_wind_used)); % DRO最壞分布期望附加項(xiàng) obj_dro lambda * param.alpha (1/param.N_scen) * sum(s_var); objective obj_det obj_dro;核心約束張力在場(chǎng)景環(huán)節(jié)。對(duì)每個(gè)離散場(chǎng)景is_var(i)要大于等于“當(dāng)前決策在場(chǎng)景i下的成本 - lambda乘以場(chǎng)景距離調(diào)整項(xiàng)”。這個(gè)調(diào)整項(xiàng)是用來(lái)約束真實(shí)分布不能離經(jīng)驗(yàn)分布太遠(yuǎn)的“軟約束”。constr []; constr [constr, lambda 0]; for i 1:param.N_scen % 提取場(chǎng)景i的風(fēng)電出力 wind_i scen(i, :); % 該場(chǎng)景下的最壞成本上界簡(jiǎn)化示意 cost_i sum(param.c_wind * (wind_i - P_wind_used)) ... % 棄風(fēng)懲罰項(xiàng) sum(param.c_pen * (P_load - P_supply_i)); % 切負(fù)荷懲罰項(xiàng) % Wasserstein距離項(xiàng)簡(jiǎn)化為場(chǎng)景差分的1-范數(shù)按需調(diào)整 dist_i sum(abs(scen(i,:) - scen_mean), 2) / T; constr [constr, s_var(i) cost_i - lambda * dist_i]; end注意這只是一個(gè)示意片段不同項(xiàng)目的目標(biāo)函數(shù)、懲罰系數(shù)、場(chǎng)景維度差異很大代碼要根據(jù)自己系統(tǒng)的約束修改。模型組裝好后直接調(diào)求解器ops sdpsettings(solver, gurobi, verbose, 2); sol optimize(constr, objective, ops); % 檢查求解狀態(tài) if sol.problem 0 disp(求解成功); else disp(sol.info); end3.4 結(jié)果后處理與對(duì)比實(shí)驗(yàn)計(jì)算結(jié)果別只看最優(yōu)成本一個(gè)數(shù)。我做這類項(xiàng)目時(shí)一般會(huì)畫(huà)三張圖第一張CHP、電鍋爐、購(gòu)電、風(fēng)電的24小時(shí)電功率平衡堆疊圖。重點(diǎn)看有沒(méi)有出現(xiàn)“CHP貼著下限運(yùn)行還不得不棄風(fēng)”的時(shí)段如果存在說(shuō)明模糊集半徑或儲(chǔ)熱罐調(diào)度參數(shù)可能需要調(diào)整。第二張熱功率平衡圖疊加儲(chǔ)熱罐SOC曲線。這張圖能直觀看出儲(chǔ)熱罐有沒(méi)有起到“削峰填谷”作用。如果SOC曲線全程貼著上限或下限跑說(shuō)明儲(chǔ)熱罐容量沒(méi)被合理利用可以考慮調(diào)整充放熱的價(jià)格參數(shù)。第三張成本對(duì)比條形圖。分別跑確定性模型、SAA隨機(jī)規(guī)劃模型、DRO模型畫(huà)出總成本再在DRO模型里取幾個(gè)不同的alpha值看成本是怎么隨著alpha增大而升高的。這張圖是論文或匯報(bào)里最有說(shuō)服力的結(jié)果。我在實(shí)際項(xiàng)目中測(cè)下來(lái)alpha從0.05增大到0.3成本大約會(huì)上升5%-15%但系統(tǒng)的“實(shí)際失負(fù)荷小時(shí)數(shù)”會(huì)大幅下降。這個(gè)trade-off曲線建議你在匯報(bào)時(shí)重點(diǎn)展示比堆一堆公式更能讓導(dǎo)師或者甲方理解分布魯棒的工程價(jià)值。4. 調(diào)試排坑我在這個(gè)項(xiàng)目上踩過(guò)的五個(gè)坑4.1 對(duì)偶變量維度不匹配導(dǎo)致YALMIP報(bào)錯(cuò)最常見(jiàn)的問(wèn)題s_var定義成標(biāo)量但場(chǎng)景循環(huán)里卻按向量索引使用。YALMIP對(duì)維度非常敏感一旦s_var(i)沒(méi)法索引直接報(bào)“Subscripted assignment dimension mismatch”。解決技巧在定義變量后用assert語(yǔ)句檢查維度assert(length(s_var) param.N_scen, s_var維度與場(chǎng)景數(shù)不匹配);這個(gè)檢查放在建模前能讓你第一時(shí)間定位是定義問(wèn)題還是約束問(wèn)題。4.2 模糊集半徑alpha不是我拍腦袋定的alpha太小模型形同虛設(shè)基本等價(jià)SAAalpha太大成本膨脹嚴(yán)重失去意義。理論上有經(jīng)驗(yàn)公式alpha和樣本量N的關(guān)系是alpha O(1/sqrt(N))但實(shí)際操作中我推薦交叉驗(yàn)證法。把歷史數(shù)據(jù)切分成訓(xùn)練集和驗(yàn)證集。用訓(xùn)練集構(gòu)造經(jīng)驗(yàn)分布和模糊集求解調(diào)度決策再拿驗(yàn)證集里沒(méi)參與建模的“真實(shí)場(chǎng)景”去回測(cè)看成本分布情況。選能覆蓋90%驗(yàn)證場(chǎng)景不切負(fù)荷的最小alpha。這個(gè)方法雖然要多花一點(diǎn)時(shí)間但勝在可解釋性很強(qiáng)匯報(bào)時(shí)也容易被接受。4.3 熱功率平衡約束導(dǎo)致無(wú)解電鍋爐和儲(chǔ)熱罐同時(shí)參與熱平衡時(shí)很容易出現(xiàn)“熱量來(lái)源太多”導(dǎo)致熱功率過(guò)剩約束沖突無(wú)解。我排查過(guò)多次原因基本都出在儲(chǔ)熱罐的充放熱0-1約束沒(méi)寫(xiě)完整充熱和放熱變量同時(shí)為正熱平衡被雙倍計(jì)入。一個(gè)排查技巧求解無(wú)解時(shí)先把儲(chǔ)熱罐的0-1約束和SOC約束單獨(dú)拿出來(lái)固定所有電出力變量只求解熱子系統(tǒng)。如果熱子系統(tǒng)仍有解再逐步加回電側(cè)約束用二分法鎖定沖突約束。4.4 場(chǎng)景數(shù)一多求解時(shí)間爆炸50個(gè)場(chǎng)景、24個(gè)時(shí)段決策變量里再帶上0-1儲(chǔ)熱變量Gurobi解一個(gè)MILP可能要20分鐘。后來(lái)我的處理方式是先用大場(chǎng)景數(shù)做預(yù)分析確定合適的alpha范圍正式求解時(shí)把場(chǎng)景聚類數(shù)量壓到30個(gè)左右同時(shí)給MILP設(shè)置一個(gè)相對(duì)最優(yōu)間隙(mipgap)比如5%很多情況下Gurobi能在5分鐘解決戰(zhàn)斗。實(shí)際工程決策場(chǎng)景下5%的間隙完全可接受花15分鐘追求0.1%的經(jīng)濟(jì)提升其實(shí)意義不大。4.5 數(shù)據(jù)中心化處理不當(dāng)導(dǎo)致模糊集失效這個(gè)坑比較隱晦。經(jīng)驗(yàn)分布P_hat在DRO中的位置非常關(guān)鍵一旦場(chǎng)景數(shù)據(jù)沒(méi)有按變量均值中心化Wasserstein距離計(jì)算出的alpha實(shí)際含義會(huì)偏離預(yù)期。解決方法是構(gòu)造模糊集前的所有場(chǎng)景均做零均值、單位方差標(biāo)準(zhǔn)化等對(duì)偶轉(zhuǎn)化完成、得到?jīng)Q策結(jié)果后再把結(jié)果反標(biāo)準(zhǔn)化回實(shí)際物理量綱。我在這里吃過(guò)一次虧花了一周時(shí)間比對(duì)結(jié)果最后發(fā)現(xiàn)是歸一化以后忘了在距離項(xiàng)里乘回尺度系數(shù)導(dǎo)致alpha的實(shí)際幾何意義完全不對(duì)。5. 后續(xù)擴(kuò)展這套框架還可以往哪走數(shù)據(jù)驅(qū)動(dòng)分布魯棒這套框架的價(jià)值不局限于電熱綜合能源系統(tǒng)。我做完這個(gè)項(xiàng)目后發(fā)現(xiàn)同樣的“離散場(chǎng)景 Wasserstein模糊集 對(duì)偶轉(zhuǎn)化”三段式可以直接平移到很多相關(guān)問(wèn)題上含氫儲(chǔ)能的綜合能源系統(tǒng)調(diào)度氫氣儲(chǔ)能的不確定性和儲(chǔ)熱罐很像但時(shí)間尺度更長(zhǎng)電動(dòng)汽車(chē)聚合商參與電力市場(chǎng)的投標(biāo)策略充電行為的不確定性正好用場(chǎng)景描述園區(qū)級(jí)微電網(wǎng)與配電網(wǎng)的互動(dòng)優(yōu)化分布式光伏出力波動(dòng)比風(fēng)電更劇烈分布魯棒的優(yōu)勢(shì)更能體現(xiàn)。如果想把項(xiàng)目做成真正的論文級(jí)別還可以加一套兩階段分布魯棒模型第一階段決定機(jī)組開(kāi)停機(jī)和儲(chǔ)熱罐充放計(jì)劃第二階段在不確定性實(shí)現(xiàn)后做出力調(diào)整。兩階段的DRO模型更貼近真實(shí)調(diào)度流程但求解復(fù)雜度會(huì)顯著上升需要引入Benders分解或者割平面方法這是另一個(gè)值得寫(xiě)一篇文章的話題了。我自己實(shí)際操作下來(lái)的體會(huì)是分布魯棒優(yōu)化最難的部分不是數(shù)學(xué)推導(dǎo)而是對(duì)不確定性數(shù)據(jù)的認(rèn)知——你對(duì)數(shù)據(jù)越了解模糊集半徑就選得越準(zhǔn)優(yōu)化結(jié)果也就越有說(shuō)服力。如果一上來(lái)就追求復(fù)雜的模糊集和花哨的求解器反而容易忽略問(wèn)題的本質(zhì)。先把確定性模型吃透再把場(chǎng)景注入最后加模糊集兜底一步一個(gè)腳印這個(gè)方向其實(shí)沒(méi)有想象中那么高不可攀。 Matlab實(shí)現(xiàn)電熱綜合能源系統(tǒng)多場(chǎng)景分布魯棒優(yōu)化調(diào)度方法
先聊一個(gè)真實(shí)調(diào)度場(chǎng)景早上8點(diǎn)預(yù)測(cè)下午3點(diǎn)風(fēng)大你按這個(gè)預(yù)測(cè)排好了機(jī)組出力結(jié)果到了2點(diǎn)半風(fēng)突然變小電鍋爐還在滿功率吃電儲(chǔ)熱罐可用熱量又快見(jiàn)底這時(shí)候再調(diào)整機(jī)組爬坡已經(jīng)來(lái)不及了。做電熱綜合能源調(diào)度的人多少都經(jīng)歷過(guò)這種“預(yù)測(cè)一時(shí)爽調(diào)度兩行淚”的時(shí)刻。純靠點(diǎn)估計(jì)做決策遲早被不確定性教做人但完全按最壞場(chǎng)景規(guī)劃成本高得離譜電廠和熱網(wǎng)公司都接受不了。于是這些年“數(shù)據(jù)驅(qū)動(dòng) 分布魯棒優(yōu)化”這套思路越來(lái)越熱核心就是利用歷史數(shù)據(jù)構(gòu)造場(chǎng)景再在場(chǎng)景周?chē)粢粋€(gè)“模糊集”來(lái)兜底最終在“經(jīng)濟(jì)性”和“魯棒性”之間取一個(gè)可調(diào)的平衡點(diǎn)。這篇文章圍繞多離散場(chǎng)景分布魯棒Distributionally Robust Optimization, DRO方法聊聊怎么在電熱綜合能源系統(tǒng)里落地一套完整的Matlab調(diào)度程序。我會(huì)從問(wèn)題建模、不確定性處理、模糊集構(gòu)造、對(duì)偶轉(zhuǎn)化、代碼框架到踩坑實(shí)錄全部過(guò)一遍適合正在做綜合能源調(diào)度、園區(qū)能源管理相關(guān)課題的研究生以及想從確定性優(yōu)化轉(zhuǎn)向魯棒優(yōu)化的工程師參考。1. 先把問(wèn)題說(shuō)透電熱綜合能源系統(tǒng)為什么難優(yōu)化1.1 電和熱兩種能量耦合不是兩套獨(dú)立系統(tǒng)電熱綜合能源系統(tǒng)之所以讓人頭疼是因?yàn)殡?、熱本身就是?qiáng)耦合的。典型園區(qū)系統(tǒng)里主要有這么幾類設(shè)備熱電聯(lián)產(chǎn)機(jī)組CHP這是電-熱耦合的核心。抽汽式CHP的發(fā)電量和供熱量必須在可行域內(nèi)聯(lián)動(dòng)發(fā)多了電往往熱也會(huì)多靈活性受限。電鍋爐EB把電轉(zhuǎn)化成熱相當(dāng)于給系統(tǒng)增加了一條“電轉(zhuǎn)熱”的靈活路徑但它本質(zhì)上是負(fù)荷會(huì)增加電網(wǎng)購(gòu)電壓力。儲(chǔ)熱罐TES儲(chǔ)熱罐是解耦電-熱“強(qiáng)綁定”的關(guān)鍵緩沖裝置熱負(fù)荷高峰時(shí)放熱熱負(fù)荷低谷時(shí)充熱相當(dāng)于給系統(tǒng)提供一個(gè)時(shí)間維度上的平移能力。常規(guī)火電/上級(jí)電網(wǎng)購(gòu)電提供電功率平衡的兜底手段。風(fēng)電/光伏這里是不確定性的主要來(lái)源也直接決定了為什么要做分布魯棒。電網(wǎng)側(cè)需要滿足CHP發(fā)電 風(fēng)電實(shí)際出力 購(gòu)電 電負(fù)荷 電鍋爐耗電。熱網(wǎng)側(cè)需要滿足CHP供熱量 電鍋爐供熱量 儲(chǔ)熱罐放熱 熱負(fù)荷 儲(chǔ)熱罐充熱 熱網(wǎng)損失。問(wèn)題的難點(diǎn)不在方程多而在于CHP 的發(fā)電和供熱捆綁在一起風(fēng)電又隨機(jī)波動(dòng)。一旦風(fēng)電預(yù)測(cè)不準(zhǔn)系統(tǒng)只能通過(guò)快速調(diào)整CHP或者改變電鍋爐出力來(lái)平衡但這兩者又牽動(dòng)著熱力平衡牽一發(fā)動(dòng)全身。1.2 三種不確定性處理思路為什么我選分布魯棒對(duì)風(fēng)電出力和負(fù)荷預(yù)測(cè)誤差的處理現(xiàn)在主流有三條路線它們之間存在本質(zhì)區(qū)別方法不確定集合優(yōu)化目標(biāo)對(duì)數(shù)據(jù)的依賴保守程度確定性優(yōu)化無(wú)直接使用預(yù)測(cè)值單場(chǎng)景經(jīng)濟(jì)最優(yōu)只要預(yù)測(cè)值最低但容易失衡隨機(jī)規(guī)劃(SAA)若干離散場(chǎng)景及概率場(chǎng)景期望成本最小需要大量場(chǎng)景中低依賴場(chǎng)景質(zhì)量傳統(tǒng)魯棒優(yōu)化盒式/橢球不確定域最壞情況下成本最小只需不確定變量的上下界高過(guò)于保守分布魯棒優(yōu)化(DRO)以歷史分布為中心的模糊集最壞概率分布下的期望成本最小需要?dú)v史數(shù)據(jù)構(gòu)造場(chǎng)景和模糊集可調(diào)的中間水平隨機(jī)規(guī)劃的問(wèn)題是你給的場(chǎng)景概率分布是固定的但實(shí)際風(fēng)電誤差分布可能跟假設(shè)偏差很大傳統(tǒng)魯棒的問(wèn)題是它只關(guān)心最壞情況凡是落在盒式集合內(nèi)的可能值都按最壞結(jié)果兜底最終調(diào)度成本可能比實(shí)際最優(yōu)貴20%以上。分布魯棒走的是中間路線——只假設(shè)真實(shí)分布落在經(jīng)驗(yàn)分布周?chē)囊粋€(gè)“模糊集”內(nèi)然后優(yōu)化最壞分布下的期望成本。當(dāng)歷史數(shù)據(jù)足夠多時(shí)模糊集半徑可以收得很小結(jié)果趨近隨機(jī)規(guī)劃當(dāng)數(shù)據(jù)不充分或者你對(duì)分布沒(méi)信心時(shí)可以調(diào)大半徑結(jié)果向傳統(tǒng)魯棒靠攏。這個(gè)“可調(diào)”特性在工程里非常實(shí)用因?yàn)樗屇隳苡靡惶姿惴蚣苋ミm配不同數(shù)據(jù)質(zhì)量的項(xiàng)目。1.3 多離散場(chǎng)景這個(gè)“多”字到底指什么標(biāo)題里的“多離散場(chǎng)景”可以拆成兩個(gè)層面理解。第一層是用有限離散場(chǎng)景近似經(jīng)驗(yàn)分布把歷史風(fēng)速、歷史負(fù)荷誤差數(shù)據(jù)通過(guò)聚類或者采樣歸納成幾十個(gè)代表性場(chǎng)景每個(gè)場(chǎng)景對(duì)應(yīng)一組24小時(shí)的風(fēng)電出力曲線。這樣隨機(jī)變量就從“連續(xù)分布”退化為“帶概率權(quán)重的離散點(diǎn)集”計(jì)算規(guī)?？煽?。第二層是模糊集也按離散場(chǎng)景構(gòu)造分布魯棒優(yōu)化不需要顯式描述連續(xù)分布而是以“這些離散場(chǎng)景組成的經(jīng)驗(yàn)分布”為中心向外擴(kuò)一個(gè)半徑為alpha的距離球通常是Wasserstein距離。只要真實(shí)分布沒(méi)有超出這個(gè)球最壞情況下的期望成本就在模型掌控范圍內(nèi)。這種“離散場(chǎng)景 模糊球”的組合既保留了隨機(jī)規(guī)劃對(duì)場(chǎng)景信息的利用又繼承了魯棒優(yōu)化對(duì)不確定性的兜底能力是這個(gè)方向近幾年成為熱點(diǎn)的關(guān)鍵。2. 分布魯棒優(yōu)化的建模與轉(zhuǎn)化怎么把min-max問(wèn)題變成可求解問(wèn)題2.1 先建立確定性調(diào)度主模型整套優(yōu)化模型可以表達(dá)為目標(biāo)最小化系統(tǒng)總運(yùn)行成本包括CHP燃料成本、向上級(jí)電網(wǎng)購(gòu)電成本、棄風(fēng)懲罰、切負(fù)荷懲罰。約束電功率平衡、熱功率平衡、CHP可行域約束、電鍋爐出力上下限、儲(chǔ)熱罐SOC連續(xù)性約束與容量約束、購(gòu)電上下限、爬坡約束等。其中CHP運(yùn)行可行域是關(guān)鍵約束通常用一組線性不等式圍成的多邊形描述即滿足P_chp_min P_chp_t P_chp_maxH_chp_min H_chp_t H_chp_maxP_chp_t c1 * H_chp_t c2c3 * H_chp_t - P_chp_t c4這四個(gè)不等式基本可以刻畫(huà)抽汽式CHP的可行域。如果你的CHP是背壓式可以直接簡(jiǎn)化為P c * H那電熱耦合更緊調(diào)度難度更大。儲(chǔ)熱罐模型相對(duì)簡(jiǎn)單但要注意不能同時(shí)充放熱的約束S_t S_{t-1} eta_c * H_char_t - H_dis_t / eta_d - eta_loss * S_{t-1}0 H_char_t H_char_max * u_c_t0 H_dis_t H_dis_max * u_d_tu_c_t u_d_t 1這里的u_c_t和u_d_t是0-1變量一旦加進(jìn)去模型就變成了MILP混合整數(shù)線性規(guī)劃。如果系統(tǒng)規(guī)模不大直接用YALMIP Gurobi求解沒(méi)問(wèn)題如果規(guī)模大要考慮用啟發(fā)式或者滾動(dòng)時(shí)域控制去壓規(guī)模。2.2 把不確定性裝進(jìn)去從SAA到DRO引入風(fēng)電不確定性后簡(jiǎn)化的目標(biāo)函數(shù)形式變成一個(gè)min-max雙層結(jié)構(gòu)外層min是調(diào)度決策變量x機(jī)組出力、儲(chǔ)熱罐充放、購(gòu)電等內(nèi)層max是針對(duì)模糊集D中所有可能概率分布P最大化期望成本。DRO目標(biāo)含義min_x max_{P in D} E_P[f(x, xi)]在“最壞的合理分布”下期望成本最小這里xi就是隨機(jī)變量比如風(fēng)電誤差向量、負(fù)荷誤差向量D就是模糊集。模糊集最常見(jiàn)的構(gòu)造方式是Wasserstein球D { P | W(P, P_hat) alpha }其中P_hat是歷史數(shù)據(jù)構(gòu)造的經(jīng)驗(yàn)分布alpha是模糊集半徑。Wasserstein距離衡量?jī)蓚€(gè)分布之間的“搬運(yùn)成本”直觀理解就是要把經(jīng)驗(yàn)分布P_hat變成真實(shí)分布P最少需要“搬運(yùn)”多少概率質(zhì)量。alpha越大真實(shí)分布離經(jīng)驗(yàn)分布可以越遠(yuǎn)模型越保守。工程實(shí)現(xiàn)中我用的是1-范數(shù)Wasserstein距離因?yàn)閷?duì)偶轉(zhuǎn)化后能保持線性約束結(jié)構(gòu)Gurobi和CPLEX可以直接吃下。如果上2-范數(shù)模型會(huì)帶二次約束求解器壓力明顯變大除非你確有必要我不太推薦。2.3 對(duì)偶轉(zhuǎn)化的三板斧直接求解min-max問(wèn)題是不現(xiàn)實(shí)的實(shí)際做法是把它對(duì)偶成一個(gè)單層的min問(wèn)題。以Wasserstein模糊集為例核心思路分三步第一步內(nèi)層max對(duì)偶變換。在凸目標(biāo)函數(shù)和凸模糊集條件下內(nèi)層最大化問(wèn)題可以等價(jià)地轉(zhuǎn)化為一個(gè)關(guān)于對(duì)偶變量lambda的最小化問(wèn)題。這一步讓“最壞分布”消失換來(lái)的是對(duì)每個(gè)離散場(chǎng)景引入一組額外約束和輔助變量s_i。第二步引入場(chǎng)景級(jí)最優(yōu)值函數(shù)。每個(gè)離散場(chǎng)景xi_i都會(huì)產(chǎn)生一個(gè)“最壞成本”的支撐函數(shù)通過(guò)引入輔助變量s_i把對(duì)每個(gè)場(chǎng)景的max項(xiàng)線性化。第三步最終得到形如min_{x, lambda0, s_i} lambda * alpha (1/N) * sum(s_i)約束lambda 0對(duì)每個(gè)場(chǎng)景is_i f(x, xi_i) - lambda * d(xi_i, xi_j) 的某種線性化形式這里d是場(chǎng)景之間的Wasserstein距離項(xiàng)f是當(dāng)前決策x在場(chǎng)景xi_i下的運(yùn)行成本。這個(gè)公式寫(xiě)出來(lái)可能有點(diǎn)抽象但在代碼里的操作其實(shí)很簡(jiǎn)單不要自己手推全部對(duì)偶約束而是用小規(guī)模測(cè)試3個(gè)場(chǎng)景、3個(gè)時(shí)段去驗(yàn)證對(duì)偶轉(zhuǎn)換后的目標(biāo)值和暴力枚舉max結(jié)果一致。我一開(kāi)始直接套大模型結(jié)果怎么都不對(duì)最后縮小規(guī)模一步步檢查才定位到是某個(gè)極端場(chǎng)景下的約束沒(méi)寫(xiě)全。3. Matlab代碼實(shí)現(xiàn)從零搭一套分布魯棒調(diào)度程序3.1 代碼框架與文件規(guī)劃整個(gè)項(xiàng)目我分成五個(gè)文件維護(hù)結(jié)構(gòu)清晰也方便復(fù)現(xiàn)項(xiàng)目目錄/ ├── main.m % 主程序入口 ├── set_parameters.m % 設(shè)備參數(shù)定義 ├── generate_scenarios.m % 歷史數(shù)據(jù)讀取 場(chǎng)景生成/聚類 ├── build_dro_model.m % YALMIP建模仿真 ├── plot_results.m % 結(jié)果可視化main.m 的結(jié)構(gòu)就是典型的“參數(shù)-場(chǎng)景-建模-求解-畫(huà)圖”五段式%% main.m clc; clear; close all; % 1. 參數(shù)設(shè)置 param set_parameters(); % 2. 讀歷史風(fēng)速數(shù)據(jù)生成離散場(chǎng)景 [scen, prob] generate_scenarios(param); % 3. 構(gòu)建DRO模型并求解 [result, model] build_dro_model(param, scen, prob); % 4. 畫(huà)圖 plot_results(result, param);這里有一點(diǎn)值得強(qiáng)調(diào)場(chǎng)景生成和建模是解耦的。如果你后續(xù)想換數(shù)據(jù)、換聚類方法只需要?jiǎng)觛enerate_scenarios這一個(gè)文件不需要碰主模型這個(gè)設(shè)計(jì)能給你后續(xù)調(diào)參省很多事。3.2 場(chǎng)景生成K-means聚類 拉丁超立方采樣場(chǎng)景生成模塊的核心目標(biāo)是把大量歷史風(fēng)電出力數(shù)據(jù)壓縮成幾十個(gè)代表性離散場(chǎng)景。我用的是兩步法。第一步拉丁超立方采樣(LHS)生成海量預(yù)測(cè)誤差樣本。為什么要用LHS而不是直接蒙特卡洛抽樣因?yàn)長(zhǎng)HS能保證樣本點(diǎn)在概率空間內(nèi)覆蓋更均勻同樣的樣本量下LHS構(gòu)造的經(jīng)驗(yàn)分布更穩(wěn)定模糊集半徑alpha可以選得更小。第二步用K-means對(duì)樣本聚類取聚類中心作為代表場(chǎng)景按每個(gè)簇的樣本比例分配概率。聚類數(shù)N我一般取20~50之間。少于20個(gè)場(chǎng)景分布信息損失太嚴(yán)重模糊集半徑會(huì)被迫加大成本偏向保守多于50個(gè)場(chǎng)景模型規(guī)模膨脹明顯MILP求解時(shí)間會(huì)從幾分鐘漲到一個(gè)小時(shí)。核心代碼示意% generate_scenarios.m 片段 % 歷史誤差樣本: err_hist (N_hist x T) err_hist load_historical_data(); % LHS生成候選樣本 candidate lhsdesign(N_samp, T); % 將[0,1]區(qū)間的LHS樣本映射到經(jīng)驗(yàn)誤差分布 err_samp quantile(err_hist, candidate); % K-means聚類 [idx, C] kmeans(err_samp, N_scen, Replicates, 10); prob histcounts(idx, N_scen) / N_samp; % C就是N_scen x T的離散場(chǎng)景矩陣 scen C;一個(gè)容易忽略的坑聚類前要對(duì)數(shù)據(jù)進(jìn)行歸一化尤其是風(fēng)電和負(fù)荷量綱不同的時(shí)候。有的風(fēng)電場(chǎng)裝機(jī)容量500MW負(fù)荷可能才50MW如果不歸一化聚類結(jié)果會(huì)完全被風(fēng)電數(shù)據(jù)主導(dǎo)負(fù)荷誤差場(chǎng)景根本體現(xiàn)不出來(lái)。我通常對(duì)每個(gè)隨機(jī)變量單獨(dú)歸一化到[0,1]區(qū)間構(gòu)造模糊集后再反歸一化回去。3.3 核心模型YALMIP寫(xiě)DRO的關(guān)鍵代碼模型部分用YALMIP建模求解器接口用Gurobi或CPLEX。整個(gè)DRO模型在YALMIP里非常直白因?yàn)榉植际紧敯艮D(zhuǎn)化后的模型本質(zhì)上是一個(gè)帶額外輔助變量的MILP。決策變量定義片段% build_dro_model.m 片段 T param.T; P_chp sdpvar(T, 1); % CHP電出力 H_chp sdpvar(T, 1); % CHP熱出力 P_eb sdpvar(T, 1); % 電鍋爐耗電 H_eb sdpvar(T, 1); % 電鍋爐供熱 P_buy sdpvar(T, 1); % 購(gòu)電 SOC sdpvar(T1, 1); % 儲(chǔ)熱罐狀態(tài) H_char sdpvar(T, 1); % 充熱 H_dis sdpvar(T, 1); % 放熱 % 輔助變量DRO對(duì)偶變量和場(chǎng)景上界 lambda sdpvar(1); s_var sdpvar(param.N_scen, 1);目標(biāo)函數(shù)的YALMIP寫(xiě)法% 確定性成本部分 obj_det sum(param.c_fuel .* P_chp) sum(param.c_buy .* P_buy) ... sum(param.c_wind * (scen_forecast - P_wind_used)); % DRO最壞分布期望附加項(xiàng) obj_dro lambda * param.alpha (1/param.N_scen) * sum(s_var); objective obj_det obj_dro;核心約束張力在場(chǎng)景環(huán)節(jié)。對(duì)每個(gè)離散場(chǎng)景is_var(i)要大于等于“當(dāng)前決策在場(chǎng)景i下的成本 - lambda乘以場(chǎng)景距離調(diào)整項(xiàng)”。這個(gè)調(diào)整項(xiàng)是用來(lái)約束真實(shí)分布不能離經(jīng)驗(yàn)分布太遠(yuǎn)的“軟約束”。constr []; constr [constr, lambda 0]; for i 1:param.N_scen % 提取場(chǎng)景i的風(fēng)電出力 wind_i scen(i, :); % 該場(chǎng)景下的最壞成本上界簡(jiǎn)化示意 cost_i sum(param.c_wind * (wind_i - P_wind_used)) ... % 棄風(fēng)懲罰項(xiàng) sum(param.c_pen * (P_load - P_supply_i)); % 切負(fù)荷懲罰項(xiàng) % Wasserstein距離項(xiàng)簡(jiǎn)化為場(chǎng)景差分的1-范數(shù)按需調(diào)整 dist_i sum(abs(scen(i,:) - scen_mean), 2) / T; constr [constr, s_var(i) cost_i - lambda * dist_i]; end注意這只是一個(gè)示意片段不同項(xiàng)目的目標(biāo)函數(shù)、懲罰系數(shù)、場(chǎng)景維度差異很大代碼要根據(jù)自己系統(tǒng)的約束修改。模型組裝好后直接調(diào)求解器ops sdpsettings(solver, gurobi, verbose, 2); sol optimize(constr, objective, ops); % 檢查求解狀態(tài) if sol.problem 0 disp(求解成功); else disp(sol.info); end3.4 結(jié)果后處理與對(duì)比實(shí)驗(yàn)計(jì)算結(jié)果別只看最優(yōu)成本一個(gè)數(shù)。我做這類項(xiàng)目時(shí)一般會(huì)畫(huà)三張圖第一張CHP、電鍋爐、購(gòu)電、風(fēng)電的24小時(shí)電功率平衡堆疊圖。重點(diǎn)看有沒(méi)有出現(xiàn)“CHP貼著下限運(yùn)行還不得不棄風(fēng)”的時(shí)段如果存在說(shuō)明模糊集半徑或儲(chǔ)熱罐調(diào)度參數(shù)可能需要調(diào)整。第二張熱功率平衡圖疊加儲(chǔ)熱罐SOC曲線。這張圖能直觀看出儲(chǔ)熱罐有沒(méi)有起到“削峰填谷”作用。如果SOC曲線全程貼著上限或下限跑說(shuō)明儲(chǔ)熱罐容量沒(méi)被合理利用可以考慮調(diào)整充放熱的價(jià)格參數(shù)。第三張成本對(duì)比條形圖。分別跑確定性模型、SAA隨機(jī)規(guī)劃模型、DRO模型畫(huà)出總成本再在DRO模型里取幾個(gè)不同的alpha值看成本是怎么隨著alpha增大而升高的。這張圖是論文或匯報(bào)里最有說(shuō)服力的結(jié)果。我在實(shí)際項(xiàng)目中測(cè)下來(lái)alpha從0.05增大到0.3成本大約會(huì)上升5%-15%但系統(tǒng)的“實(shí)際失負(fù)荷小時(shí)數(shù)”會(huì)大幅下降。這個(gè)trade-off曲線建議你在匯報(bào)時(shí)重點(diǎn)展示比堆一堆公式更能讓導(dǎo)師或者甲方理解分布魯棒的工程價(jià)值。4. 調(diào)試排坑我在這個(gè)項(xiàng)目上踩過(guò)的五個(gè)坑4.1 對(duì)偶變量維度不匹配導(dǎo)致YALMIP報(bào)錯(cuò)最常見(jiàn)的問(wèn)題s_var定義成標(biāo)量但場(chǎng)景循環(huán)里卻按向量索引使用。YALMIP對(duì)維度非常敏感一旦s_var(i)沒(méi)法索引直接報(bào)“Subscripted assignment dimension mismatch”。解決技巧在定義變量后用assert語(yǔ)句檢查維度assert(length(s_var) param.N_scen, s_var維度與場(chǎng)景數(shù)不匹配);這個(gè)檢查放在建模前能讓你第一時(shí)間定位是定義問(wèn)題還是約束問(wèn)題。4.2 模糊集半徑alpha不是我拍腦袋定的alpha太小模型形同虛設(shè)基本等價(jià)SAAalpha太大成本膨脹嚴(yán)重失去意義。理論上有經(jīng)驗(yàn)公式alpha和樣本量N的關(guān)系是alpha O(1/sqrt(N))但實(shí)際操作中我推薦交叉驗(yàn)證法。把歷史數(shù)據(jù)切分成訓(xùn)練集和驗(yàn)證集。用訓(xùn)練集構(gòu)造經(jīng)驗(yàn)分布和模糊集求解調(diào)度決策再拿驗(yàn)證集里沒(méi)參與建模的“真實(shí)場(chǎng)景”去回測(cè)看成本分布情況。選能覆蓋90%驗(yàn)證場(chǎng)景不切負(fù)荷的最小alpha。這個(gè)方法雖然要多花一點(diǎn)時(shí)間但勝在可解釋性很強(qiáng)匯報(bào)時(shí)也容易被接受。4.3 熱功率平衡約束導(dǎo)致無(wú)解電鍋爐和儲(chǔ)熱罐同時(shí)參與熱平衡時(shí)很容易出現(xiàn)“熱量來(lái)源太多”導(dǎo)致熱功率過(guò)剩約束沖突無(wú)解。我排查過(guò)多次原因基本都出在儲(chǔ)熱罐的充放熱0-1約束沒(méi)寫(xiě)完整充熱和放熱變量同時(shí)為正熱平衡被雙倍計(jì)入。一個(gè)排查技巧求解無(wú)解時(shí)先把儲(chǔ)熱罐的0-1約束和SOC約束單獨(dú)拿出來(lái)固定所有電出力變量只求解熱子系統(tǒng)。如果熱子系統(tǒng)仍有解再逐步加回電側(cè)約束用二分法鎖定沖突約束。4.4 場(chǎng)景數(shù)一多求解時(shí)間爆炸50個(gè)場(chǎng)景、24個(gè)時(shí)段決策變量里再帶上0-1儲(chǔ)熱變量Gurobi解一個(gè)MILP可能要20分鐘。后來(lái)我的處理方式是先用大場(chǎng)景數(shù)做預(yù)分析確定合適的alpha范圍正式求解時(shí)把場(chǎng)景聚類數(shù)量壓到30個(gè)左右同時(shí)給MILP設(shè)置一個(gè)相對(duì)最優(yōu)間隙(mipgap)比如5%很多情況下Gurobi能在5分鐘解決戰(zhàn)斗。實(shí)際工程決策場(chǎng)景下5%的間隙完全可接受花15分鐘追求0.1%的經(jīng)濟(jì)提升其實(shí)意義不大。4.5 數(shù)據(jù)中心化處理不當(dāng)導(dǎo)致模糊集失效這個(gè)坑比較隱晦。經(jīng)驗(yàn)分布P_hat在DRO中的位置非常關(guān)鍵一旦場(chǎng)景數(shù)據(jù)沒(méi)有按變量均值中心化Wasserstein距離計(jì)算出的alpha實(shí)際含義會(huì)偏離預(yù)期。解決方法是構(gòu)造模糊集前的所有場(chǎng)景均做零均值、單位方差標(biāo)準(zhǔn)化等對(duì)偶轉(zhuǎn)化完成、得到?jīng)Q策結(jié)果后再把結(jié)果反標(biāo)準(zhǔn)化回實(shí)際物理量綱。我在這里吃過(guò)一次虧花了一周時(shí)間比對(duì)結(jié)果最后發(fā)現(xiàn)是歸一化以后忘了在距離項(xiàng)里乘回尺度系數(shù)導(dǎo)致alpha的實(shí)際幾何意義完全不對(duì)。5. 后續(xù)擴(kuò)展這套框架還可以往哪走數(shù)據(jù)驅(qū)動(dòng)分布魯棒這套框架的價(jià)值不局限于電熱綜合能源系統(tǒng)。我做完這個(gè)項(xiàng)目后發(fā)現(xiàn)同樣的“離散場(chǎng)景 Wasserstein模糊集 對(duì)偶轉(zhuǎn)化”三段式可以直接平移到很多相關(guān)問(wèn)題上含氫儲(chǔ)能的綜合能源系統(tǒng)調(diào)度氫氣儲(chǔ)能的不確定性和儲(chǔ)熱罐很像但時(shí)間尺度更長(zhǎng)電動(dòng)汽車(chē)聚合商參與電力市場(chǎng)的投標(biāo)策略充電行為的不確定性正好用場(chǎng)景描述園區(qū)級(jí)微電網(wǎng)與配電網(wǎng)的互動(dòng)優(yōu)化分布式光伏出力波動(dòng)比風(fēng)電更劇烈分布魯棒的優(yōu)勢(shì)更能體現(xiàn)。如果想把項(xiàng)目做成真正的論文級(jí)別還可以加一套兩階段分布魯棒模型第一階段決定機(jī)組開(kāi)停機(jī)和儲(chǔ)熱罐充放計(jì)劃第二階段在不確定性實(shí)現(xiàn)后做出力調(diào)整。兩階段的DRO模型更貼近真實(shí)調(diào)度流程但求解復(fù)雜度會(huì)顯著上升需要引入Benders分解或者割平面方法這是另一個(gè)值得寫(xiě)一篇文章的話題了。我自己實(shí)際操作下來(lái)的體會(huì)是分布魯棒優(yōu)化最難的部分不是數(shù)學(xué)推導(dǎo)而是對(duì)不確定性數(shù)據(jù)的認(rèn)知——你對(duì)數(shù)據(jù)越了解模糊集半徑就選得越準(zhǔn)優(yōu)化結(jié)果也就越有說(shuō)服力。如果一上來(lái)就追求復(fù)雜的模糊集和花哨的求解器反而容易忽略問(wèn)題的本質(zhì)。先把確定性模型吃透再把場(chǎng)景注入最后加模糊集兜底一步一個(gè)腳印這個(gè)方向其實(shí)沒(méi)有想象中那么高不可攀。