綜合能源系統(tǒng)最優(yōu)熱力潮流（OTF）模型建立與Gurobi實現(xiàn)
1. 最優(yōu)熱力潮流的定義和最優(yōu)電力潮流（OPF）類似，最優(yōu)熱力潮流（Optimal Thermal Flow，以下簡稱OTF）指當熱力管網(wǎng)（包含冷網(wǎng)和熱網(wǎng)）結(jié)構(gòu)參數(shù)
1. 最優(yōu)熱力潮流的定義
和最優(yōu)電力潮流（OPF）類似，最優(yōu)熱力潮流（Optimal Thermal Flow，以下簡稱OTF）指當熱力管網(wǎng)（包含冷網(wǎng)和熱網(wǎng)）結(jié)構(gòu)參數(shù)和冷熱負荷情況都給定的情況下，通過調(diào)節(jié)可利用的控制變量（分為質(zhì)調(diào)節(jié)和量調(diào)節(jié)）來找到能滿足所有運行約束條件的，并使系統(tǒng)的某一性能指標（如供能成本）達到最優(yōu)值下的潮流分布。以我們所研究的case為例，其OTF問題可以描述為：節(jié)點1是包含了一臺燃氣鍋爐和一臺CHP機組的熱源，節(jié)點4、5、6存在熱力負荷且負荷值已知，如何調(diào)節(jié)熱源的熱出力功率，使得熱力系統(tǒng)能夠在滿足末端負荷需求的前提下最小化運行成本。
需要注意的是，和OPF相比，OTF有三個十分明顯的區(qū)別，這也使得OTF模型的構(gòu)建過程和OPF完全不同，概括如下：
- 控制變量不同。熱力管網(wǎng)的控制方式主要可以分為質(zhì)調(diào)節(jié)和量調(diào)節(jié)，前者指調(diào)節(jié)供水/回水溫度，后者指調(diào)節(jié)供水流量。值得一提的是，目前國內(nèi)熱力管網(wǎng)大多采用質(zhì)調(diào)節(jié)的控制方式；
- 管路數(shù)量不同。熱力管網(wǎng)分為供水網(wǎng)絡(luò)和回水網(wǎng)絡(luò)，在計算時需要同時考慮供水管網(wǎng)和回水管網(wǎng)各節(jié)點和管段的能量守恒方程；
- 傳輸性質(zhì)不同。熱力管網(wǎng)傳輸能量的介質(zhì)是水，水的傳播速度遠低于電的傳播速度，且熱力管網(wǎng)往往較長，因此需要考慮溫度的延遲效應。同時，水在傳輸過程中會和環(huán)境產(chǎn)生熱交換，造成溫度的損耗效應。
2. OTF模型建立
本節(jié)考慮的OTF模型為質(zhì)調(diào)節(jié)控制方法，即供回水流量保持定值，通過調(diào)節(jié)供回水溫度的方式向末端負荷供能。模型同時考慮了溫度的延遲效應和溫度的損耗效應。OTF模型中的節(jié)點可以分為三類：熱源節(jié)點、熱負荷節(jié)點和中間節(jié)點。由于存在供水管道和回水管道，因此每個節(jié)點包含進口供水溫度、出口供水溫度、進口回水溫度和出口回水溫度四個控制變量。從約束條件的角度而言，每個節(jié)點具有供水和回水溫度混合約束，這些方程同時考慮了水在傳遞過程中的延遲效應和損耗效應。此外，熱源節(jié)點和熱負荷節(jié)點還具有能量守恒約束。
2.1 供水/回水溫度混合約束
根據(jù)能量守恒定律，任意一個節(jié)點供水/回水的出口溫度乘出口總流量都必須等于所有入口溫度和對應流量乘積的總和，如下式所示：
式中，
延遲效應采用簡化的節(jié)點法進行建模，假設(shè)每個調(diào)度間隔對應的水量為一個質(zhì)量塊，如下圖所示，則調(diào)度時刻 
式中， 
這個式子可以通過上圖直接得到，非常直觀且容易理解：
接下來考慮溫度的損耗效應，對其進行建模。熱力管網(wǎng)在傳輸時會向周圍環(huán)境散熱，溫度損耗和環(huán)境溫度、傳熱系數(shù)等因素有關(guān)，可以由以下公式確定：
式中， 
2.2 熱源節(jié)點和熱負荷節(jié)點的能量守恒約束
對于熱源節(jié)點和熱負荷節(jié)點，其熱出力/熱負荷需要和供回水溫差建立起能量守恒約束，如下式所示：
式中， 
3. OTF模型的代碼實現(xiàn)
相比于OPF模型（其代碼實現(xiàn)可見上一期專欄），OTF代碼實現(xiàn)的主要區(qū)別包括：1. 需要同時考慮供水管網(wǎng)和回水管網(wǎng)；2. 上一個節(jié)點的出口參數(shù)（即溫度）并不等于下一個節(jié)點的入口參數(shù)，需要考慮溫度延遲和損耗；3. 不同類型的節(jié)點需要考慮的約束條件不同。
對于本專欄所研究的case，我們定義如下待優(yōu)化變量，包括每個節(jié)點的冷卻水出口溫度和冷凍水出口溫度和燃氣鍋爐的熱功率，CHP機組的熱功率由于可以由其電功率唯一確定因此忽略：
h_b = m.addVar(lb=b_list[0][2], ub=b_list[0][3], name='heating output of boiler')nts = m.addVars(n_node, vtype=GRB.CONTINUOUS, lb=lb_ts, ub=ub_ts, name='supply temperature from bus 1 to bus 6')ntr = m.addVars(n_node, vtype=GRB.CONTINUOUS, lb=lb_tr, ub=ub_tr, name='return temperature from bus 1 to bus 6')
代碼的核心難點是如何快速高效地編寫約束條件。這里仍然采用鄰接矩陣的思想，首先分別生成供水管網(wǎng)矩陣、回水管網(wǎng)矩陣和節(jié)點類型矩陣，代碼如下：
def generateMmatrix(n_node, br_sr):n M = np.zeros((n_node, n_node, 2))n for branch in br_sr:n M[branch[0].astype(int) - 1][branch[1].astype(int) - 1][0] = branch[2]n M[branch[0].astype(int) - 1][branch[1].astype(int) - 1][1] = exp(-branch[3] * branch[4] / (cp * 1000 * branch[2]))n return Mndef generateTypeArrays(node_list):n node_type1 = np.empty(shape=[0, 2])n node_type2 = np.empty(shape=[0, 2])n node_type3 = np.empty(shape=[0, 2])n for node in node_list:n if node[1] == 1:n node_type1 = np.append(node_type1, [node], axis=0)n elif node[1] == 2:n node_type2 = np.append(node_type2, [node], axis=0)n else:n node_type3 = np.append(node_type3, [node], axis=0)n return node_type1, node_type2, node_type3nnM_ts = generateMmatrix(n_node, br_s)nM_tr = generateMmatrix(n_node, br_r)nnode_type1, node_type2, node_type3 = generateTypeArrays(node_list)
簡單說一下這兩個方法都干了什么事：generateMmatrix方法根據(jù)傳入的節(jié)點數(shù)量n_node和管路信息br_sr生成一個三維的M矩陣，若節(jié)點1到節(jié)點2之間存在一條支路，則在對應的行列中存入一個數(shù)組，第一個元素為管路流量，第二個元素為耗散效應的那一坨指數(shù)，這樣方便直接調(diào)用計算；generateTypeArrays方法則是根據(jù)傳入的節(jié)點類型生成三種類型的節(jié)點列表，type1是熱源節(jié)點，type2是負荷節(jié)點，type3是中間節(jié)點。
之后，就可以用addConstrs方法一次性添加各節(jié)點的單個類型約束，代碼如下：
# energy conservation equation for the nodes of type 1nm.addConstrs(cp / 1000 * sum(M_ts[node_type1[i][0].astype(int) - 1][j][0] for j in range(n_node)) * n (ts[node_type1[i][0].astype(int) - 1] - tr[node_type1[i][0].astype(int) - 1]) ==n h_chp_set[i] + h_b_set[i] for i in range(node_type1.shape[0]))n# energy conservation equation for the nodes of type 2nm.addConstrs(cp / 1000 * sum(M_ts[j][node_type2[i][0].astype(int) - 1][0] for j in range(n_node)) * n (ts[node_type2[i][0].astype(int) - 1] - tr[node_type2[i][0].astype(int) - 1]) ==n hl[i] for i in range(node_type2.shape[0]))n# supply temperature mixed equationnm.addConstrs(sum(M_ts[j][i][0] * ((ts[j] - t_amb) * M_ts[j][i][1] + t_amb) for j in range(n_node)) ==n sum(M_ts[j][i][0] for j in range(n_node)) * ts[i] for i in range(n_node))n# return temperature mixed equationnm.addConstrs(sum(M_tr[j][i][0] * ((tr[j] - t_amb) * M_tr[j][i][1] + t_amb) for j in range(n_node)) ==n sum(M_tr[j][i][0] for j in range(n_node)) * tr[i] for i in range(n_node))
前兩個語句添加的是熱源節(jié)點和熱負荷節(jié)點的能量守恒約束，這里將兩類節(jié)點分別保存在兩個列表中，從而便于計算；后兩個語句添加的是各節(jié)點供水管網(wǎng)和回水管網(wǎng)的溫度混合約束，這里等式左邊的節(jié)點入口溫度原則上需要同時考慮延遲效應和耗散效應，代碼中只考慮了耗散效應，原因如下：
當考慮溫度延遲效應時，由于當前時刻的溫度本質(zhì)上取決于先前時刻的溫度，這時優(yōu)化問題會從時間無關(guān)的變?yōu)闀r間相關(guān)的，時間無關(guān)的優(yōu)化問題在做優(yōu)化時可以只把當前調(diào)度時刻的控制變量作為優(yōu)化變量，而時間相關(guān)的優(yōu)化問題需要同時將未來一段時間內(nèi)所有控制變量作為優(yōu)化變量，這會使優(yōu)化問題的復雜度大大提升。同樣會使優(yōu)化問題變?yōu)闀r間相關(guān)的模型還有儲能模型，因為儲能裝置的SOC取決于上一時刻的充放電狀態(tài)。出于簡單考慮，本case原型代碼中忽略溫度延遲效應（其實是我懶，本來想考慮一下但是需要修改的代碼量太大）。這里給出這類問題的編程思路：很簡單，把優(yōu)化變量、約束條件、目標函數(shù)以及隨時間變化的數(shù)據(jù)（末端負荷）等在原有基礎(chǔ)上增加一個時間的維度。在代碼實現(xiàn)上，就是在原有數(shù)組的基礎(chǔ)上加一維，其他不變。
在構(gòu)建了OPF模型和OTF模型后，我們終于可以寫出這個問題的目標函數(shù)，代碼如下：
obj = price_gas * (p_chp / effi_g2p + h_b / b_list[0][1]) + g_list[0][1] * p_g1 * p_g1 + g_list[0][2] * p_g1 + n g_list[1][1] * p_g2 * p_g2 + g_list[1][2] * p_g2 + g_list[2][1] * p_g3 * p_g3 + g_list[2][2] * p_g3nm.setObjective(obj, GRB.MINIMIZE)
目標函數(shù)包含兩部分：第一部分是購氣成本，天然氣用于兩部分，即驅(qū)動CHP機組和驅(qū)動燃氣鍋爐；第二部分是三臺發(fā)電機的發(fā)電成本。至此，該case的優(yōu)化問題已經(jīng)完整的寫出，可以執(zhí)行運行獲得最優(yōu)的運行結(jié)果。
4. 完整代碼
完整代碼已經(jīng)全部上傳至Gitee，可直接下載跑來看一看，歡迎star和fork。一共有三個py文件：data、generate和main，優(yōu)化運行的代碼往往比較長，所以建議大家可以分成多個文件，既清晰也方便修改。我習慣于把case的所有數(shù)據(jù)都放在一個文件里，把操作函數(shù)放在一個文件里，主程序直接寫gurobi代碼，數(shù)據(jù)和函數(shù)可以直接調(diào)用，很方便。除了代碼外，還附帶了所測試系統(tǒng)的原始數(shù)據(jù)（testdata_CHPD.xls）。
5. 結(jié)語
不知不覺這個專欄終于結(jié)束了，總共耗時將近一個半月。隨著這個專欄，我也順帶重新梳理了目前常用的模型，同時改進了一些自己之前寫的代碼，常用常新。祝大家都能按時、順利畢業(yè)！
