搞這個(gè)方向的人應(yīng)該都有過類似的體驗(yàn)看到一篇EI論文摘要里寫著“提出了一種改進(jìn)粒子群算法求解建筑集成光儲系統(tǒng)規(guī)劃運(yùn)行綜合優(yōu)化模型”理論上思路很清晰可真到自己動(dòng)手用Python復(fù)現(xiàn)的時(shí)候光是把這個(gè)模型從論文公式變成能跑通的代碼就夠折騰好幾天。我這次復(fù)現(xiàn)的就是這類模型——基于改進(jìn)粒子群算法的建筑集成光儲系統(tǒng)規(guī)劃運(yùn)行綜合優(yōu)化目標(biāo)是在同一個(gè)框架里同時(shí)回答兩個(gè)問題光伏和儲能到底該配多大容量規(guī)劃層以及配好之后儲能每天怎么充放電、建筑跟電網(wǎng)怎么交互運(yùn)行層。這篇文章會把完整的復(fù)現(xiàn)過程拆開講先分析模型到底在優(yōu)化什么再把目標(biāo)函數(shù)和約束條件翻譯成代碼邏輯然后說清楚標(biāo)準(zhǔn)粒子群算法PSO為什么不夠用、改進(jìn)在哪里最后給出可運(yùn)行的Python實(shí)現(xiàn)要點(diǎn)和幾個(gè)復(fù)現(xiàn)時(shí)最容易踩的坑。適合正在做光儲優(yōu)化、微電網(wǎng)調(diào)度方向論文復(fù)現(xiàn)的讀者也適合工程師想快速驗(yàn)證一個(gè)新算法方案在真實(shí)場景下是否可行。1. 這個(gè)題到底在優(yōu)化什么規(guī)劃決策與運(yùn)行策略的雙層耦合1.1 規(guī)劃層在決定什么建筑集成光儲系統(tǒng)通俗講就是在建筑屋頂裝光伏板、在配電房或者地下室裝儲能電池讓建筑既可以自己發(fā)電又可以存電還能跟電網(wǎng)買賣電。規(guī)劃層的決策變量比較直白光伏的裝機(jī)容量kW、儲能電池的額定容量kWh、儲能變流器的額定功率kW。有時(shí)候還包含逆變器選型、變壓器擴(kuò)容與否這類0-1變量。這些變量一旦定下來后續(xù)十幾年的收益上限基本就被鎖死了。光伏裝少了白天發(fā)的電不夠用儲能也沒多少電量可存儲能容量配大了初始投資高得嚇人一年到頭可能就閑置在那里。1.2 運(yùn)行層在決定什么運(yùn)行層解決的是“怎么用”的問題。在一個(gè)典型日里光伏出力曲線是一個(gè)正弦波樣子的序列負(fù)荷曲線是建筑里空調(diào)、照明、辦公設(shè)備疊加出來的另一個(gè)序列儲能則需要在每個(gè)時(shí)段決定充電還是放電、充多少放多少。同時(shí)建筑還要決定跟電網(wǎng)買電還是賣電以及各自多少功率。所以運(yùn)行層的決策變量是全天每個(gè)時(shí)段通常按1小時(shí)或15分鐘劃分的儲能充放電功率、購電功率、售電功率。如果按4個(gè)典型日代表全年四季那么運(yùn)行變量的維度就是4典型日× 24小時(shí)× 若干個(gè)動(dòng)作變量直接撐起了優(yōu)化問題里最大的變量部分。1.3 為什么“規(guī)劃運(yùn)行”必須一起算這倆階段聽著像是可以分開做的先粗算一下容量需求再在固定容量下做運(yùn)行調(diào)度。但這樣做的結(jié)果往往是偏保守甚至互相矛盾。舉個(gè)例子如果我一開始把儲能容量設(shè)得過大運(yùn)行調(diào)度時(shí)會發(fā)現(xiàn)最優(yōu)策略是每天最多用掉一小半電量因?yàn)殡妰r(jià)差不足以覆蓋深度充放電帶來的效率損耗。反過來如果運(yùn)行策略只追求最大化自消納率那就會傾向于把儲能耗盡長期滿充滿放會加速電池衰減運(yùn)行維護(hù)成本上升又反過來影響規(guī)劃段的經(jīng)濟(jì)性評估。這就是典型的強(qiáng)耦合問題必須在一個(gè)統(tǒng)一的優(yōu)化模型里同時(shí)找容量和運(yùn)行策略的“最優(yōu)組合”。2. 目標(biāo)函數(shù)與約束條件怎么轉(zhuǎn)成代碼能算的數(shù)學(xué)模型2.1 目標(biāo)函數(shù)經(jīng)濟(jì)性指標(biāo)怎么取舍復(fù)現(xiàn)這類論文時(shí)最常見的目標(biāo)函數(shù)是年綜合費(fèi)用最小化。它由四部分組成初始投資的年化值、年運(yùn)行維護(hù)費(fèi)用、年購電費(fèi)用、減去年售電收入。初始投資要用資金回收系數(shù)CRF折算成年值公式是CRF r × (1r)^L / ((1r)^L - 1)其中r是貼現(xiàn)率L是設(shè)備壽命。光伏和儲能壽命不一樣光伏通常按25年、儲能按10年左右計(jì)所以折算的時(shí)候得分開算。這一步很多人會漏直接拿總投資除以壽命年限誤差會很大。用公式表達(dá)就是min C_total CRF_pv × C_pv × P_pv CRF_bat × C_bat × E_bat C_power × P_pcs C_om Σ_t (price_buy(t) × P_buy(t) × Δt) ? Σ_t (price_sell(t) × P_sell(t) × Δt)運(yùn)行維護(hù)費(fèi)一般按初始投資的百分比比如1%~2%來估算購電分時(shí)電價(jià)和售電電價(jià)則是兩個(gè)不同的價(jià)格序列。2.2 等式約束與不等式約束的矩陣化壓縮模型里約束條件不少但歸納起來就幾類我在代碼里全部用numpy矩陣和向量處理避免寫for循環(huán)串行判斷導(dǎo)致速度太慢。第一類是功率平衡約束每個(gè)時(shí)段都要滿足P_pv(t) P_dis(t) P_buy(t) P_load(t) P_ch(t) P_sell(t)四個(gè)變量里有兩個(gè)自由度代碼里最好另外兩個(gè)變量表達(dá)P_buy(t) P_load(t) P_ch(t) P_sell(t) ? P_pv(t) ? P_dis(t)這樣平衡約束就被自動(dòng)消去了優(yōu)化的搜索空間能縮小一圈。第二類是儲能系統(tǒng)約束包括SOC遞推公式SOC(t1) SOC(t) η_ch × P_ch(t) × Δt / E_bat ? P_dis(t) × Δt / (η_dis × E_bat)以及SOC上下限0.1~0.9、充放電功率上下限不超過P_pcs、同一時(shí)段不能同時(shí)充放電。同一時(shí)段不能同時(shí)充放電這種約束簡化處理可以用一個(gè)很小的正數(shù)罰項(xiàng)嚴(yán)謹(jǐn)處理則引入0-1變量但那樣會變成混合整數(shù)問題粒子群就得配合修復(fù)策略來解。第三類是電網(wǎng)交互約束購電功率上限和售電功率上限通常設(shè)為同一個(gè)值比如200kW。第四類是決策變量邊界約束光伏裝機(jī)容量不超過屋頂可安裝面積折算的上限儲能容量有離散檔位的話就做就近取整處理。2.3 典型參數(shù)的取值參考復(fù)現(xiàn)時(shí)我用了下面這組參數(shù)基本符合多數(shù)EI論文的量級跑出來的結(jié)果趨勢上是合理的參數(shù)取值說明光伏單位造價(jià)4500元/kW含安裝儲能單位容量造價(jià)1800元/kWh磷酸鐵鋰PCS單位功率造價(jià)1200元/kW雙向變流器貼現(xiàn)率5%年化折算用光伏壽命25年CRF折算儲能壽命10年CRF折算分時(shí)購電電價(jià)峰1.2/平0.75/谷0.38元/kWh三費(fèi)率售電電價(jià)0.65元/kWh固定采樣間隔1小時(shí)24個(gè)時(shí)段/天典型日數(shù)量4個(gè)按季度聚類這些參數(shù)不需要跟原論文完全一致因?yàn)樵紨?shù)據(jù)拿不到但量級一定要對否則結(jié)果很容易出現(xiàn)“儲能配了0容量”這種退化結(jié)果——那往往不是算法問題而是分時(shí)電價(jià)差不夠覆蓋儲能投資成本。3. 標(biāo)準(zhǔn)PSO為什么不夠用改進(jìn)粒子群的三個(gè)關(guān)鍵改動(dòng)3.1 慣性權(quán)重與學(xué)習(xí)因子的動(dòng)態(tài)調(diào)整標(biāo)準(zhǔn)粒子群的核心公式大家都知道v w × v c1 × r1 × (pbest ? x) c2 × r2 × (gbest ? x) x x v但標(biāo)準(zhǔn)PSO在光儲規(guī)劃這個(gè)問題上有兩個(gè)先天不足一是搜索后期粒子速度衰減很快容易還沒找到好位置就全體聚集在局部最優(yōu)附近二是運(yùn)行變量維度動(dòng)輒上百維標(biāo)準(zhǔn)PSO幾乎不具備全局搜索能力。所以論文里的“改進(jìn)”通常第一步就是慣性權(quán)重非線性遞減。線性遞減雖然簡單但前期全局搜索時(shí)間不夠后期局部收斂力度也不夠。我自己試下來效果最好的是冪函數(shù)遞減w (w_max ? w_min) × ((k_max ? k) / k_max)^alpha w_minw_max取0.9w_min取0.4alpha取1.5。迭代前40%期間w始終維持在0.7以上保證粒子有足夠的速度去探索不同容量組合到后期w快速跌向0.4讓粒子群體集中打磨最優(yōu)區(qū)域。學(xué)習(xí)因子也做了互補(bǔ)變化。c1表示向個(gè)體歷史最優(yōu)學(xué)習(xí)的程度c2表示向全局最優(yōu)學(xué)習(xí)的程度。前期c1應(yīng)該大一點(diǎn)讓粒子自己做主后期c2應(yīng)該大起來讓粒子更信任群體經(jīng)驗(yàn)。實(shí)測c1從2.5線性降到0.5、c2從0.5升到2.5的效果比固定雙因子好很多。3.2 約束處理方式罰函數(shù)還是修復(fù)策略改進(jìn)PSO與標(biāo)準(zhǔn)PSO最大的區(qū)別往往藏在約束處理里。規(guī)劃運(yùn)行優(yōu)化模型中的等式約束和不等式約束非常多如果不處理就直接用粒子群去搜大概率在前幾十次迭代里幾乎找不到一個(gè)可行解。罰函數(shù)法是最通用的方案把每個(gè)違反約束的量求平方后乘以一個(gè)大罰因子加進(jìn)目標(biāo)函數(shù)里。這樣做的好處是簡單不用改動(dòng)粒子群的結(jié)構(gòu)缺點(diǎn)是罰因子太小約束會放松罰因子太大會形成“硬墻”粒子飛不進(jìn)去收斂曲線會非常難看。我自己復(fù)現(xiàn)時(shí)采用的是“動(dòng)態(tài)罰函數(shù)儲能SOC修復(fù)”的混合策略。具體來說SOC越界不直接罰而是把SOC序列重新縮放到[0.1, 0.9]區(qū)間再做功率平衡修正。這個(gè)操作能大幅提高粒子群搜索效率。功率平衡通過2.2節(jié)的變量替換自動(dòng)滿足。充放電功率越限壓縮到上下限邊界。對仍然無法恢復(fù)的約束再計(jì)入罰項(xiàng)。這個(gè)混合策略的思路是能用物理邏輯修復(fù)的就先修復(fù)修復(fù)不了的才交給算法去罰這樣粒子群不會把大量算力浪費(fèi)在尋找可行域入口上。3.3 混合離散-連續(xù)變量的特殊編碼規(guī)劃變量里的儲能容量在實(shí)際采購時(shí)是有離散檔位的比如50kWh一檔但光伏容量和功率變量是連續(xù)的。如果全部按連續(xù)變量搜索最后取整會產(chǎn)生誤差如果全部按離散變量處理維度和搜索復(fù)雜度會爆炸。我采用的方式是把粒子位置向量拆成三段來編碼前2維是光伏容量和儲能容量連續(xù)值第三維是PCS功率連續(xù)值后面的維度是各個(gè)典型日每個(gè)時(shí)段的充放電功率和購售電功率連續(xù)值。每次迭代結(jié)束以后只對儲能容量做取整到最近檔位的操作其他變量保持連續(xù)。這樣既保持了粒子群的連續(xù)搜索特性又保證了最終給出來的方案是工程上能落地采購的離散容量。4. Python實(shí)現(xiàn)的分模塊拆解從數(shù)據(jù)處理到迭代求解4.1 數(shù)據(jù)準(zhǔn)備與典型日場景提取原始數(shù)據(jù)通常是8760小時(shí)的建筑負(fù)荷序列和光伏出力理論值序列。復(fù)現(xiàn)時(shí)我沒有直接拿全年8760小時(shí)去做優(yōu)化那樣不僅計(jì)算量爆炸而且相鄰時(shí)段的高度相關(guān)性會讓粒子群陷入維度災(zāi)難。先用KMeans做典型日聚類。把每天24個(gè)負(fù)荷值24個(gè)光伏出力值拼成48維的特征向量然后用sklearn.cluster.KMeans聚成4類代表春夏秋冬四種典型日。聚類完成后每類的樣本數(shù)量占比就是該典型日在全年優(yōu)化中的權(quán)重。需要注意KMeans對初始聚類中心敏感我固定了random_state42并且跑了多次驗(yàn)證聚類結(jié)果的穩(wěn)定性。from sklearn.cluster import KMeans import numpy as np # X_daily shape: (365, 48)前24維是負(fù)荷后24維是光伏 kmeans KMeans(n_clusters4, random_state42, n_init20).fit(X_daily) labels kmeans.labels_ weights np.bincount(labels, minlength4) / 365.0 typical_days np.array([X_daily[kmeans.labels_ i].mean(axis0) for i in range(4)])聚類完成之后4個(gè)典型日的24小時(shí)負(fù)荷序列和24小時(shí)光伏出力序列就是后續(xù)運(yùn)行優(yōu)化的基礎(chǔ)輸入。權(quán)重weights則用來把4個(gè)典型日的運(yùn)行費(fèi)用折算成年費(fèi)用這樣目標(biāo)函數(shù)才能跟“全年”的經(jīng)濟(jì)賬對得上。4.2 粒子編碼與種群初始化粒子編碼是整個(gè)代碼最核心的數(shù)據(jù)結(jié)構(gòu)。我定義每個(gè)粒子為一個(gè)一維numpy數(shù)組長度由三部分組成3規(guī)劃變量 4典型日× 24時(shí)段× 3儲能充電功率、放電功率、購電功率。實(shí)際算下來長度是3 4×24×3 291維每個(gè)粒子291維種群規(guī)模設(shè)40個(gè)這個(gè)維度下粒子群勉強(qiáng)還能跑但已經(jīng)能明顯感到收斂速度變慢。初始化的時(shí)候不能純隨機(jī)不然大量粒子初始解嚴(yán)重違反功率平衡罰函數(shù)值會大到把所有粒子都拖在不可行域邊緣。我的做法是先根據(jù)當(dāng)天光伏出力和負(fù)荷曲線用啟發(fā)式規(guī)則生成一個(gè)粗糙的初始運(yùn)行策略光伏不足且電價(jià)高就放電光伏富余且電價(jià)低就充電。把這個(gè)粗糙解作為其中一個(gè)粒子的初始位置其余粒子在這個(gè)解周圍加隨機(jī)擾動(dòng)生成。這樣做的好處是種群一開始就有幾個(gè)粒子落在可行域附近收斂速度明顯加快。4.3 迭代主循環(huán)與多目標(biāo)處理主循環(huán)就是一個(gè)典型的PSO迭代框架但每次迭代里要依次完成計(jì)算所有粒子的目標(biāo)函數(shù)與罰項(xiàng)、更新個(gè)體最優(yōu)pbest、更新全局最優(yōu)gbest、按動(dòng)態(tài)參數(shù)的公式更新速度和位置、對位置做邊界處理和儲能SOC修復(fù)。多目標(biāo)處理方面如果原論文是同時(shí)優(yōu)化經(jīng)濟(jì)性和碳排放我復(fù)現(xiàn)時(shí)用的策略是加權(quán)和法把碳排放量乘以一個(gè)碳價(jià)因子折算成費(fèi)用加進(jìn)目標(biāo)函數(shù)。這樣粒子群只需要維護(hù)單目標(biāo)的排序關(guān)系實(shí)現(xiàn)簡單而且改動(dòng)目標(biāo)函數(shù)不影響算法主體邏輯。追求帕累托前沿的話需要引入非支配排序代碼量會增加不少不利于快速復(fù)現(xiàn)驗(yàn)證。核心迭代代碼大致是這樣的結(jié)構(gòu)for k in range(max_iter): w (w_max - w_min) * ((max_iter - k) / max_iter) ** 1.5 w_min c1 2.5 - 2.0 * k / max_iter c2 0.5 2.0 * k / max_iter for i in range(pop_size): # 計(jì)算目標(biāo)值含罰函數(shù) cost[i] objective_function(positions[i], typical_days, weights) # 更新個(gè)體最優(yōu) if cost[i] pbest_cost[i]: pbest_cost[i] cost[i] pbest_pos[i] positions[i].copy() # 更新全局最優(yōu) if cost[i] gbest_cost: gbest_cost cost[i] gbest_pos positions[i].copy() for i in range(pop_size): r1, r2 np.random.rand(2) velocities[i] w * velocities[i] c1 * r1 * (pbest_pos[i] - positions[i]) \ c2 * r2 * (gbest_pos - positions[i]) positions[i] velocities[i] positions[i] repair_and_bounds(positions[i], typical_days)每一輪迭代結(jié)束后把gbest_cost記錄到列表里后期直接畫收斂曲線。建議每一輪都打印一次當(dāng)前最優(yōu)配置和最優(yōu)費(fèi)用方便肉眼判斷是否卡在局部最優(yōu)。4.4 結(jié)果可視化與敏感性分析跑完之后要輸出的結(jié)果至少有三樣一是收斂曲線這是證明改進(jìn)算法有效性的第一手證據(jù)二是最優(yōu)規(guī)劃結(jié)果表包含光伏容量、儲能容量、PCS功率和年化總費(fèi)用三是典型日的儲能充放電曲線和購售電曲線用來檢查運(yùn)行策略是否符合直覺比如應(yīng)該在光伏富余時(shí)充電、在傍晚電價(jià)峰值時(shí)放電。我習(xí)慣額外加一個(gè)敏感性分析腳本把光伏單位造價(jià)從4500元/kWh降到2000元/kWh看最優(yōu)配置里光伏容量會不會上升儲能容量會不會因此改變。這種分析寫論文時(shí)非常好用能直觀體現(xiàn)模型的工程指導(dǎo)意義而且實(shí)現(xiàn)成本也就一個(gè)for循環(huán)的事。5. 復(fù)現(xiàn)EI論文時(shí)最容易翻車的幾個(gè)點(diǎn)5.1 典型日聚類的坑最容易翻車的環(huán)節(jié)其實(shí)是數(shù)據(jù)預(yù)處理而不是算法本身。KMeans聚類之前如果負(fù)荷和光伏出力沒有做歸一化光伏出力的數(shù)值范圍0到幾百kW會完全壓過負(fù)荷曲線的形狀特征幾百到上千kW的變化聚出來的“典型日”基本只看負(fù)荷大小光伏特征全被淹沒了。解決辦法是先把每個(gè)特征維度做標(biāo)準(zhǔn)化再聚類聚類完成后再把典型日的數(shù)值還原回原始物理量綱。另外聚幾類也是個(gè)問題K取太大典型日代表性好但計(jì)算量成倍上漲K取太小儲能充放電的策略會被平均化到不倫不類。我試下來4個(gè)典型日是精度和速度平衡得比較好的選擇部分論文用3類也可以接受。5.2 收斂判據(jù)與早熟現(xiàn)象粒子群最典型的早熟現(xiàn)象就是收斂曲線在迭代前30%就變成一條水平線之后的迭代完全沒有任何下降。遇到這種情況第一反應(yīng)不應(yīng)該是加大迭代次數(shù)而是檢查是不是粒子多樣性丟失太快。我在調(diào)試時(shí)試過幾個(gè)方案在速度更新里加入隨機(jī)擾動(dòng)項(xiàng)、每隔一部分迭代把部分粒子重新隨機(jī)初始化、引入模擬退火式的Metropolis準(zhǔn)則接受劣解來跳出局部最優(yōu)。效果最穩(wěn)定的是“精英重啟”策略——每50輪迭代保留gbest粒子隨機(jī)初始化其他30%的粒子同時(shí)把這些新粒子的速度重置。這個(gè)策略不復(fù)雜但對打破早熟很有效。另外一定要固定隨機(jī)種子做多次重復(fù)實(shí)驗(yàn)。EI論文的復(fù)現(xiàn)結(jié)果如果每次跑都不一樣很難讓審稿人相信這個(gè)算法是穩(wěn)定的。我一般固定np.random.seed(0)跑5次獨(dú)立重復(fù)看最優(yōu)解的方差如果方差太大說明算法穩(wěn)定性有問題。5.3 結(jié)果對比要對比什么復(fù)現(xiàn)論文時(shí)最讓人頭疼的是“怎么我的結(jié)果跟論文對不上”。坦白講除非你能拿到作者的原始數(shù)據(jù)和超參數(shù)否則數(shù)值完全一致幾乎不可能審稿人預(yù)期也不是這個(gè)。有效的對比方式是跑三個(gè)版本的算法標(biāo)準(zhǔn)PSO、改進(jìn)PSO、再加一個(gè)常規(guī)的遺傳算法GA做參照畫收斂曲線對比收斂速度和最終目標(biāo)值。如果改進(jìn)PSO的最終年費(fèi)用比標(biāo)準(zhǔn)PSO低2%~5%同時(shí)收斂速度快50%以上這已經(jīng)是一個(gè)很好看的結(jié)果了。我在自己的復(fù)現(xiàn)里改進(jìn)PSO比標(biāo)準(zhǔn)PSO的年費(fèi)用降低約3.8%收斂到穩(wěn)定值的迭代次數(shù)從420輪降到260輪左右這說明動(dòng)態(tài)權(quán)重和約束修復(fù)策略確實(shí)起了作用。最后分享一個(gè)復(fù)現(xiàn)完這個(gè)模型后的感受改進(jìn)粒子群算法在這類問題里真正的價(jià)值不是追求全局最優(yōu)的數(shù)學(xué)保證而是在模型足夠復(fù)雜、商用求解器不一定吃得住的時(shí)候給你一個(gè)能在可接受時(shí)間內(nèi)找到好解的務(wù)實(shí)手段。調(diào)試過程中如果發(fā)現(xiàn)某一組參數(shù)下改進(jìn)PSO還不如標(biāo)準(zhǔn)PSO別急著懷疑算法先回看約束處理是不是逼得太狠了。我自己的習(xí)慣是始終保留一個(gè)“基礎(chǔ)參數(shù)對照”腳本每改一個(gè)改進(jìn)機(jī)制就跑一遍同一組隨機(jī)種子這樣到底哪個(gè)改動(dòng)起了作用一目了然。