1. 多時刻優化的定義在之前的專欄中，我們已經完整的建立了熱電型綜合能源系統的模型并通過gurobi實現了求解。上述模型可以認為是

1. 多時刻優化的定義

在之前的專欄中，我們已經完整的建立了熱電型綜合能源系統的模型并通過gurobi實現了求解。上述模型可以認為是最為基本的版本，有很多很多改進空間，如從單時刻優化變為多時刻優化、增加儲能模型、增加管網/建筑蓄熱、考慮需求側響應、從集中式求解變為分布式求解、增加光伏風電等新能源等等。許多SCI論文也都是以模型、算法的組合式改進為創新點進行發表的。本文將針對多時刻優化這一經典問題，對之前的基本模型進行改進。首先聲明一下：本文的觀點并不一定是完全正確的，代碼也一定不是最簡潔、最易懂的，本文旨在提供一個基本思路，供初學者入門。

多時刻優化和單時刻優化的主要區別是：單時刻優化僅僅以當前時刻的指標（如經濟成本、系統能耗等）為優化目標，而多時刻優化以未來一段時間的指標之和為優化目標。這是一個典型的模型預測控制問題，如下圖所示：

其中，優化間隔（optimization interval）一般為1h、15min等，由于一般是對未來一天進行調度，因此對應的優化時域（optimization horizon）往往分別為24、96等。

2. 多時刻優化下的建模

相比于單時刻優化，多時刻優化建模無非需要多考慮兩點：

• 如何快速表示多時刻下的目標函數和約束條件。以間隔為1h、時域為24的多時刻優化問題為例，目標函數（這里指經濟成本等指標的加和）和約束條件都將變為單時刻優化的24倍，如果要一個一個手動生成是不現實的。解決這個問題的基本思路就是巧用for循環來讓程序自動生成目標函數和約束條件。
• 如何處理時序相關的約束條件。以儲能系統為例（后續會講到具體模型，但基本原理并不難懂），t+1時刻的狀態（蓄能量）總是取決于t時刻的狀態。處理這類時序相關的約束條件，關鍵是處理好其初始態和中間態的不同表達形式。

和時序相關的約束條件有很多，此處介紹三類最常考慮的時序約束條件：

• 機組爬坡率（unit ramping rate）：一般地，產能機組在一段時間內（如熱電聯產機組、鍋爐、冷機等）的產能量變化是有限的，而非能從0%變化到100%。例如，熱電聯產機組的變化率（又稱為爬坡率）為0.5%-1.5%/min。這表示t+1時刻的機組產能量取決于t時刻的機組產能量。以之前專欄中的燃氣鍋爐模型為例，機組爬坡率約束表示如下：

式中， 為上爬坡率， 為下爬坡率， 為優化間隔， 和 分別為t時刻和t+1時刻的鍋爐熱出力。

• 熱力管網傳輸時延（transmission delay）：相比于電力網絡，熱力管網的傳輸介質（水）傳輸速度很慢，對于長距離輸運管道存在明顯傳輸時延的現象。因此，t+1時刻的管道出口溫度往往取決于t時刻甚至更早時刻的管網入口溫度。其約束形式（基于簡化節點法）已經在這篇專欄中講過（專欄中的 等同于本篇文章中的t，都表示某時刻，之后統一用t表示），此處不再贅述。本文采用基于Lax-Wendroff方法的熱力管網傳輸時延約束，該模型的優點在于它假設當前時刻的狀態只和上一時刻的狀態有關，這樣能夠簡化時序約束的建立過程，其數學表示形式如下：

? 式中， 和 分別為管道入口和出口的溫度， 為環境溫度， 為管道水流量， 為水的密度， 為管道圓橫截面積， 為管道長度， 為管道傳熱系數， 為水的比熱容。

• 儲能設備SOC（State of charging）：儲能設備能夠在能源價格處于低谷時進行儲能，并在電價高或用能高峰期釋能，不僅能夠削峰填谷，還能夠提高經濟性、促進可再生能源消納等。常見的儲能系統包含蓄電系統（Electrical energy storage，EES）、蓄熱系統（Thermal energy storage，TES）等。以EES為例，其儲能模型可以表示為以下形式：

式中， 和 分別為t時刻下的蓄能量和釋能量， 和 分別表示在 時間內能夠蓄能和釋能的上限，第三個式子表明同一時刻只能執行蓄能和釋能中的一項， 可以理解為手機電量的剩余百分比，充了就會增多，釋放就會減少，在實際應用中，它需要被控制在 和 之間以保證安全， 為EES的額定容量， 表示自損率， 和 分別表示充電效率和放電效率。

3. 多時刻優化下的編程

3.1 優化變量設置

在多時刻優化下，優化變量做如下定義：

# 以熱力管網節點的供水溫度為例：n# period為優化時域，n_node為節點個數nts = m.addVars(period, n_node, lb=lb_ts, ub=ub_ts, name='node supply temp')

每個優化變量在單時刻優化的基礎上增加了一個維度，即時間維度。可以把這個變量看作一個period行、n_node列的數組，每一行都包含了某個時刻下各個節點的供水溫度變量。在調用其中的某個變量時，可以采用如下方式調用：

# 調用ts第1行、第1列的變量nts_0_0 = ts[0, 0]

3.2 添加機組爬坡率約束

先從最簡單的說起，仍以鍋爐熱出力為例，機組爬坡率約束可以如下添加：

h_gb = m.addVars(period, lb=lb_gb, ub=ub_gb, name='gas boiler heating output')nm.addConstrs(h_gb1[t+1] - h_gb1[t] <= r_u_gb1 * interval for t in range(period - 1)) # interval為優化間隔nm.addConstrs(h_gb1[0] - h_gb1_init <= r_u_gb1 * interval)nm.addConstrs(h_gb1[t+1] - h_gb1[t] >= -r_d_gb1 * interval for t in range(period - 1))nm.addConstrs(h_gb1[0] - h_gb1_init >= -r_d_gb1 * interval)

這段代碼很好理解，關鍵在于需要單獨處理初始條件（line 3和5）。t=0時刻的鍋爐熱出力要取決于整個系統初始態下的鍋爐熱出力h_gb1_init。時序相關的約束條件都存在這個問題，也可以統一采用這種方法進行處理。注意：為使行文簡單，后續均不再單獨列出初始態下的約束條件。

3.3 添加儲能設備SOC約束

儲能SOC約束并不難添加，不過巧設變量可以降低編程難度。一般地，通常設置蓄能量 和釋能量 為優化變量，為方便編程，可以同時添加SOC作為優化變量，如下所示：

p_ch = m.addVars(period, lb=lb_ch, ub=ub_ch, name='ees charging energy')np_dch = m.addVars(period, lb=lb_dch, ub=ub_dch, name='ees discharging energy')nsoc = m.addVars(period, lb=lb_soc, ub=ub_soc, name='ees soc')nm.addConstrs(p_ch[t] * p_dch[t] == 0 for t in range(period))nm.addConstrs(soc[t+1] == soc[t] * eta_loss + interval / cap_ees * (p_ch * eta_ch - p_dch * eta_dch) n for t in range(period - 1))

可以看到，通過添加SOC優化變量，約束條件寫起來就很容易。

3.4 添加熱力管網傳輸時延約束

從模型來看，熱力管網傳輸時延約束就是一個等式約束，但是糅合了節點溫度混合和溫度耗散后，編程會變得非常復雜。以供水管網為例（回水管網和供水管網基本完全相同），各類溫度如圖所示：

從中可以看出，溫度可以分為三類：一是節點溫度，它同時也是其出水管道（即從水從該管道流出）的入口溫度；二是管道的出口溫度，該溫度由（歷史時刻的）管道入口溫度決定；三是考慮了熱損耗效應的管道實際出口溫度。通過定義這三類溫度，就能夠很容易清楚管道傳輸時延約束、節點溫度混合約束和管道熱損耗約束都表示了什么。以供水管網為例，主要代碼如下：

# 節點溫度nts = m.addVars(period, n_node, lb=lb_ts, ub=ub_ts, name='node supply temp')n# 管道出口溫度nts_pipe = m.addVars(period, n_node, n_node lb=lb_ts, ub=ub_ts, name='pipe outlet temp')n# 考慮熱損耗后的管道出口溫度nts_pipe_final = m.addVars(period, n_node, n_node lb=lb_ts, ub=ub_ts, name='pipe outlet temp after loss')nfor t in range(period):n for i in range(n_node):n # 節點溫度混合約束n # X_hn_sup為之前專欄中建立的描述熱力管網拓撲的數組，可見之前公開的源代碼。[0]表示管道流量n m.addConstr(quicksum(X_hn_sup[j][i][0] * ts_pipe_final[t, j, i] for j in range(n_node)) ==n quicksum(X_hn_sup[j][i][0] for j in range(n_node)) * ts[t, i])n for j in range(n_node):n # 若節點i和節點j之間存在管道連接，則添加管道熱損耗約束和傳輸時延約束n if X_hn_sup[i][j][0] != 0:n # 添加管道熱損耗約束n # t_amb為環境溫度數組，X_hn_sup中的[1]表示熱損耗中指數的那一坨(便于表示)n m.addConstr(ts_pipe_final[t, i, j] == (ts_pipe[t, i, j] - t_amb[t]) * X_hn_sup[i][j][1] + t_amb[t])n # 添加管道傳輸時延約束n # X_hn_sup中的[2]和[3]分別表示式中的那兩坨，同上，都是為了方便表示n m.addConstr((ts_pipe[t+1, i, j] + ts[t+1, i] - ts_pipe[t, i, j] - ts[t, i] + n X_hn_sup[i][j][2] * (ts_pipe[t+1, i, j] + ts_pipe[t, i, j] - ts[t+1, i] - ts[t, i]) n + X_hn_sup[i][j][3] * (ts_pipe[t+1, i, j] + ts_pipe[t, i, j] + ts[t+1, i] + ts[t, i] n - 4 * t_amb[t])) * X_hn_sup[i][j][0] == 0)

3.5 添加目標函數

目標函數非常簡單，只需要按照上述寫法，用quicksum函數進行累加即可。綜上，我們就較為完整地實現了綜合能源系統多時刻優化調度的編程，并且新增加了三類典型的時序相關約束，包括機組爬坡約束、儲能約束和熱力管網傳輸時延約束。臨近春節，沒有時間和精力將之前專欄中的case進行相應的修改，如果有問題歡迎在評論區或私信交流。