随着全球经济快速发展, 能源短缺与环境污染问题愈发严重。综合能源系统(Integrated Energy System, IES)能够多梯级利用能源, 可以促进可再生能源消纳, 是实现节能减排的有效途径, 可推动实现我国提出的“双碳”目标^[1-4]。

目前, 关于IES模型已有大量的相关研究。文献[5]建立了结合可再生能源以及储能装置的IES模型, 文献[6]建立了工业园区低碳经济运行的IES模型。但上述研究只建立了供能侧的模型, 未考虑到用能侧的优化。针对这一问题, 文献[7]考虑了针对电力系统的电负荷需求响应; 文献[8]引入价格弹性矩阵描述需求响应行为, 并分析了其有效性; 文献[9]引入可转移负荷对负荷侧用能进行优化。但上述研究未考虑不同负荷响应特性的情况。

IES优化调度模型具有大规模、非线性、非凸的特点。文献[10-11]采用粒子群优化(Particle Swarm Optimization, PSO)算法以及改进PSO算法优化IES模型, 但PSO算法搜索速度较慢, 精度不高, 易陷入局部最优。文献[12-13]提出用灰狼优化(Grey Wolf Optimizer, GWO)算法求解IES优化调度问题, 算法性能明显提升, 但搜索速度仍有待提高。文献[14]采用麻雀搜索算法(Sparrow Search Algorithm, SSA)对IES优化调度问题进行求解, 搜索速度快, 但初始最优解较差, 算法性能仍有待提升。

综上所述, 本文建立计及综合需求响应(Integrated Demand Response, IDR)的IES优化调度模型, 提出了一种改进的黏菌优化算法(Slime Mould Algorithm, SMA), 求解IES优化调度问题, 并结合算例进行仿真验证, 为IES的实际调度运行问题的求解提供参考。

IES架构如图 1所示。

从能量平衡角度, 建立图 1中各设备数学模型如下。

燃气轮机(Gas Turbine, GT)模型为

$ \left\{\begin{array}{l} P_{\mathrm{GT} . \mathrm{e}}=\eta_{\mathrm{GT}}^{\mathrm{e}} P_{\mathrm{GT} . \mathrm{g}} \\ P_{\mathrm{GT} . \mathrm{h}}=\eta_{\mathrm{GT}}^{\mathrm{h}} P_{\mathrm{GT} . \mathrm{g}} \end{array}\right. $

(1)

式中: P_GT.e, P_GT.g, P_GT.h——GT发电功率、耗气功率、产热功率;

η_GT^e, η_GT^h——GT发电效率、产热效率。

余热锅炉(Waste Heat Boiler, WHB)模型为

$ P_{\mathrm{WHB}}=\frac{\left(1-\eta_{\mathrm{GT}}^{\mathrm{h}}\right)}{\eta_{\mathrm{GT}}^{\mathrm{h}}} P_{\mathrm{GT} . \mathrm{h}} \eta_{\mathrm{WHB}} $

(2)

式中: P_WHB——WHB输出功率;

η_WHB——WHB产热效率。

燃气锅炉(Gas Boiler, GB)模型为

$ P_{\mathrm{GB} . \mathrm{h}}=\eta_{\mathrm{GB}}^{\mathrm{h}} P_{\mathrm{GB} . \mathrm{g}} $

(3)

式中: P_GB.h, P_GB.g——GB产热功率、耗气功率;

η_GB^h——GB产热效率。

电热锅炉(Electric Boiler, EB)模型为

$ P_{\text {EB. h }}=\eta_{\mathrm{EB}}^{\mathrm{h}} P_{\mathrm{EB}} $

(4)

式中: P_EB.h, P_EB——EB产热功率、消耗电功率;

η_EB^h——EB产热效率。

电制冷机(Electric Chiller, EC)模型为

$ P_{\mathrm{EC} . \mathrm{c}}=K_{\mathrm{COP}}^{\mathrm{EC}} P_{\mathrm{EC}} $

(5)

式中: P_EC.c, P_EC——EC输出制冷功率、输入电功率;

K_COP^EC——EC制冷系数。

吸收式制冷机(Absorption Chiller, AC)模型为

$ P_{\text {AC. c }}=K_{\mathrm{COP}}^{\mathrm{AC}} P_{\mathrm{AC}} $

(6)

式中: P_AC.c, P_AC——AC制冷功率、输入热功率;

K_COP^AC——AC制冷系数。

风力发电机(Wind Turbine, WT)模型为

$ P_{\mathrm{WT}}= \begin{cases}0, & v<v_{\text {in }} \\ P_{\mathrm{WT}}^{\mathrm{o}} \frac{v^K-v_{\mathrm{in}}^K}{v_{\mathrm{r}}^K-v_{\mathrm{in}}^K}, &\ \ v_{\mathrm{in}} \leqslant v<v_{\mathrm{r}} \\ P_{\mathrm{WT}}^{\mathrm{o}}, & v_{\mathrm{r}} \leqslant v<v_{\text {out }} \\ 0, & v \geqslant v_{\text {out }}\end{cases} $

(7)

式中: P_WT, P_WT^o——WT输出功率、额定输出功率;

v, v_in, v_r, v_out——实际风速、启动风速、额定风速、切出风速;

K——出力特性系数。

光伏(Photovoltaic, PV)模块模型为

$ P_{\mathrm{PV}}=f_{\mathrm{PV}} P_{\mathrm{n}, \mathrm{PV}} \frac{I^t}{I_{\mathrm{s}}}\left[1+\varepsilon_{\mathrm{p}}\left(t_{\mathrm{PV}}^t-t_{\mathrm{n}}\right)\right] $

(8)

式中: P_PV, P_{n, PV}——光伏模块输出功率、额定输出功率;

f_PV——能源转换效率;

I^t, I_s——t时刻辐射强度、标准辐射强度;

ε_p——温度功率系数;

t_PV^t, t_n——t时刻光伏模块温度、额定温度。

储电装置(Electric Storage, ES)模型为

$ E_{\mathrm{ES}}^t=\left(1-\mu_{\mathrm{e}}\right) E_{\mathrm{ES}}^{t-1}+\eta_{\text {in }}^{\mathrm{e}} P_{\text {in }}^{\mathrm{e}}(t)-\frac{P_{\text {out }}^{\mathrm{e}}(t)}{\eta_{\text {out }}^{\mathrm{e}}} $

(9)

式中: E_ES^t——t时刻ES充放电能量;

μ_e——ES自放电率;

η_in^e, η_out^e——ES充放电效率;

P_in^e(t), P_out^e(t)——t时刻ES充放电功率。

在负荷侧通过引导用户主动改变原有用能需求及习惯, 从而改善源荷两侧供需的不平衡关系, 提高系统运行经济性。根据负荷响应特性, 需求响应可分为价格型需求响应和替代型需求响应^[15]。

价格型需求响应模型为

$ L_{\mathrm{m}}^t=L_{\mathrm{m}, 0}^t\left(1+E_{\mathrm{m}} \frac{c_{\mathrm{m}}^t-c_{\mathrm{m}, 0}^t}{c_{\mathrm{m}, 0}^t}\right) $

(10)

式中: L_{m, 0}^t, L^t_m——t时刻响应前后负荷量;

E_m——弹性系数;

c_{m, 0}^t, c_m^t——t时刻响应前后能源价格。

替代型需求响应模型为

$ \left\{\begin{array}{l} \Delta L_t^{\mathrm{r}, \mathrm{e}}=-\varepsilon_{\mathrm{e}, \mathrm{h}} \Delta L_t^{\mathrm{r}, \mathrm{ch}} \\ \varepsilon_{\mathrm{e}, \mathrm{h}}=\frac{v_{\mathrm{e}} \varphi_{\mathrm{e}}}{v_{\mathrm{h}} \varphi_{\mathrm{h}}} \end{array}\right. $

(11)

式中: ΔL_t^{r, e}，ΔL_t^{r, ch}——t时刻可替代电负荷量、被替代冷或热负荷量;

ε_{e, h}——电热替代系数;

v_e, v_h——电能和热能的单位热值;

φ_e, φ_h——电能和热能的能源利用率。

目标函数分为以IES在运行过程中运行成本最小为目标和碳排放成本最小为目标两种。

IES运行成本C_h包括购电成本C_e、购气成本C_g以及系统运维成本C_op。

$ C_{\mathrm{h}}=C_{\mathrm{e}}+C_{\mathrm{g}}+C_{\mathrm{op}} $

(12)

购电成本C_e为

$ C_{\mathrm{e}}=\sum\limits_{t=1}^T\left(c_{\mathrm{b}, t} P_{\mathrm{b}, t}-c_{\mathrm{s}, t} P_{\mathrm{s}, t}\right) $

(13)

式中: T——运行周期时长;

c_{b, t}, c_{s, t}——t时刻购、售电电价;

P_{b, t}, P_{s, t}——t时刻向上级电网购、售电功率。

购气成本C_g为

$ C_{\mathrm{g}}=\sum\limits_{t=1}^T c_{\mathrm{g}} P_{\mathrm{g}, t} $

(14)

式中: c_g——购买天然气单价;

P_{g, t}——t时刻购买天然气功率。

系统运维成本C_op为

$ C_{\mathrm{op}}=\sum\limits_{t=1}^T \sum\limits_{i=1}^9 \omega_i P_{i, t} $

(15)

式中: i——IES设备, i=1, 2, 3, …9分别表示GT, WHB, GB, EB, EC, AC, WT, PV, ES;

ω_i——设备i的运维系数;

P_{i, t}——t时刻设备i的输出功率。

$ C_{\mathrm{H}}=W \sum\limits_{t=1}^T\left[\mu_{\mathrm{e}} P_{\mathrm{b}, t}+\mu_{\mathrm{f}} P_{\mathrm{g}, t}\right] $

(16)

式中: C_H——碳排放成本;

W——惩罚系数, 包括排放惩罚和环境价值;

μ_e——单位电功率对应的CO₂排放系数;

μ_f——单位体积天然气对应的CO₂排放系数。

设备功率约束条件为

$ \left\{\begin{array}{l} P_{\mathrm{GT}}^{\min } \leqslant P_{\mathrm{GT} \cdot \mathrm{g}} \leqslant P_{\mathrm{GT}}^{\max } \\ P_{\mathrm{GB}}^{\min } \leqslant P_{\mathrm{GB} \cdot \mathrm{h}} \leqslant P_{\mathrm{GB}}^{\max } \\ P_{\mathrm{EB}}^{\min } \leqslant P_{\mathrm{EB} \cdot \mathrm{h}} \leqslant P_{\mathrm{EB}}^{\max } \\ P_{\mathrm{EC}}^{\min } \leqslant P_{\mathrm{EC} \cdot \mathrm{c}} \leqslant P_{\mathrm{EC}}^{\max } \\ P_{\mathrm{AC}}^{\min } \leqslant P_{\mathrm{AC} \cdot \mathrm{c}} \leqslant P_{\mathrm{AC}}^{\max } \\ E_{\mathrm{ES}}^{\min } \leqslant E_{\mathrm{ES}}^t \leqslant E_{\mathrm{ES}}^{\max } \end{array}\right. $

(17)

式中: P_j^max, P_j^min——设备j输出功率上下限, 其中j为GT, GB, EB, EC, AC;

E_ES^max, E_ES^min——ES储电容量上下限。

电、热、冷能功率平衡约束条件为

$ \left\{\begin{array}{l} P_{\mathrm{GT.e}}-P_{\mathrm{EB}}-P_{\mathrm{EC}}+P_{\mathrm{WT}}+P_{\mathrm{PV}}- \\ \ \ \ \ P_{\mathrm{in}}^{\mathrm{e}}+P_{\mathrm{out}}^{\mathrm{e}}+P_{\mathrm{b}, t}-P_{\mathrm{s}, t}=L_{\mathrm{e}} \\ P_{\mathrm{WHB}}+P_{\mathrm{GB} . \mathrm{h}}+P_{\mathrm{EB} . \mathrm{h}}=L_{\mathrm{h}} \\ P_{\mathrm{EC}.\mathrm{c}}+P_{\mathrm{AC} . \mathrm{c}}=L_{\mathrm{c}} \end{array}\right. $

(18)

式中: L_e, L_h, L_c——系统的电、热、冷负荷量。

SMA通过模拟黏菌捕食实现智能寻优功能。黏菌在觅食过程中的位置变化如下:

$ \boldsymbol{X}(k+1)=\left\{\begin{array}{l} r_1 \cdot\left(\boldsymbol{U}_{\mathbf{B}}-\boldsymbol{L}_{\mathbf{B}}\right)+\boldsymbol{L}_{\mathbf{B}}, \quad r_1<z \\ \boldsymbol{X}_{\mathrm{b}}(k)+v_{\mathrm{b}} \cdot\left[\boldsymbol{W} \cdot \boldsymbol{X}_{\mathbf{A}}(k)-\boldsymbol{X}_{\mathbf{B}}(k)\right], \\ \ \ \ \ \ \ \ z \leqslant r_1<p \\ v_{\mathrm{c}} \cdot \boldsymbol{X}(k), \quad r_1 \geqslant p \end{array}\right. $

(19)

式中: X(k+1)——迭代第k+1次时更新后的黏菌位置;

r₁——随机数;

U_B, L_B——搜索边界上下限;

z——常数, 取0.03;

X_b(k)——迭代第k次时黏菌最优个体位置;

v_b, v_c——随机数, v_b∈[-a, a], v_c∈[0, 1];

W——黏菌重量;

X_A (k), X_B (k)——迭代第k次时随机选择的黏菌个体位置;

X(k)——迭代第k次时黏菌位置。

p和a的公式如下:

$ \left\{\begin{array}{l} p=\tanh \left|S-D_{\mathrm{F}}\right| \\ a=\operatorname{arctanh}\left[-\left(k / k_{\max }\right)+1\right] \end{array}\right. $

(20)

式中: S——当前黏菌个体的适应度;

D_F——迭代过程中最佳适应度;

k_max——最大迭代次数。

SMA初始化种群时, 为了增强种群多样性, 将混沌映射与反向学习相结合, 得到基于混沌对立的学习策略。本文采用Sine映射, 数学模型为

$ \boldsymbol{X}^{\mathrm{To}}=\boldsymbol{L}_{\mathrm{B}}+\boldsymbol{U}_{\mathrm{B}}+\lambda \cdot \boldsymbol{X}(k) $

(21)

式中: X^To, λ——种群中当前黏菌个体对应的混沌反向解和混沌映射值。

参数a的变化对迭代过程中开发和探索的平衡有重要影响。由式(20)可知, a在早期迭代中迅速下降, 在后期迭代中下降速度减慢, 不利于全局探索。为了提高全局探索能力和局部开发的收敛能力^[16], 可将a定义为

$ a=2.5\left(k / k_{\max }\right)^{2 k / k_{\max }} $

(22)

算术优化算法根据算术操作符的分布特性来实现全局寻优, 与SMA结合后, 可以有效改善SMA后期易陷入局部最优的缺点, 增强算法后期寻优随机性。改进SMA(以下简称“ISMA”)在更新黏菌位置时, 若r₁＜p, 则通过式(19)进行更新, 否则根据式(23)进行位置更新。

$ \boldsymbol{X}(k+1)=\left\{\begin{array}{l} \frac{\boldsymbol{X}_{\mathrm{b}}(k)}{\left(p_{\mathrm{MO}}+\varepsilon\right)} \times \\ \quad\left[\left(\boldsymbol{U}_{\mathbf{B}}-\boldsymbol{L}_{\mathbf{B}}\right) \mu+\boldsymbol{L}_{\mathbf{B}}\right], r_2<0.5 \\ \boldsymbol{X}_{\mathrm{b}}(k) p_{\mathrm{MO}} \times \\ \quad\left[\left(\boldsymbol{U}_{\mathbf{B}}-\boldsymbol{L}_{\mathbf{B}}\right) \mu+\boldsymbol{L}_{\mathbf{B}}\right], r_2 \geqslant 0.5 \end{array}\right. $

(23)

式中: p_MO——数学优化器概率;

ε——极小值;

μ——控制参数;

r₂——随机数, 且r₂∈[0, 1]。

选用Sphere、Ackley、Griewank、Rosenbrock 4种标准测试函数对PSO算法、GWO算法、SSA、SMA和ISMA的优化结果进行对比, 标准测试函数相关参数如表 1所示。

4种标准测试函数优化曲线如图 2所示。由图 2可知, 相比其他4种算法, ISMA在收敛精度和速度上具有较大的优势。

本文以某工业园区为研究对象, 每24 h为一个运行周期, 单位运行时间为1 h。系统中各设备关键参数如表 2所示, 分时购电电价如表 3所示。售电价格为0.5元/kWh, 天然气价格为2.55元/m³。系统典型日初始负荷、风电和光伏输出功率预测曲线如图 3所示。

考虑IDR前后系统的电、热、冷负荷对比曲线如图 4所示。由图 4可以看出, 原始负荷曲线显示出明显的峰、平、谷分布, 当系统考虑IDR时, 电、热、冷负荷峰谷差分别降低了22.1%, 19.3%, 17.9%, 有效实现了“削峰填谷”, 提升了IES的经济性。

为了验证ISMA的有效性, 针对系统的目标函数, 分别采用PSO算法、GWO算法、SSA、SMA和ISMA进行计算。不同算法的优化调度结果对比如表 4所示。

分析表 4中不同算法得到的最优值发现, ISMA可以找到更好的最优解, 总成本和碳排放量均低于其他算法。分别以运行成本最小和以碳排放成本最小为目标, 用ISMA得到的IES典型日目标函数的电、热、冷负荷平衡分布, 其分布情况如图 5和图 6所示。

分析图 5和图 6的ISMA优化调度结果, 可以得出以下结论。电负荷优先由可再生能源WT和PV提供, 当电能不足时, 由GT发电或从电网购买电能。当考虑系统运行成本最小为目标时, 在0:00—8:00分时电价为谷价期间, 增加向电网的购电量, 以减小系统运行成本; 当考虑系统碳排放成本最小为目标时, 由于GT的碳排放量比电网小, 因此考虑增加GT发电功率。在12:00—17:00 WT和PV发电高峰期间, 系统除满足电负荷需求外, 多余的电能可以通过EC以及EB释放冷热能, 做到多能互补, 增加风光消纳, 还可以售卖给电网。另外, ES可以在满足用户电负荷后, 存储多余的电能, 当电负荷不足时, 放电以补足电能。

对于热负荷, 优先使用WHB产生的热能, 不足部分由GB以及EB补充。

冷负荷在电能充足时段, 利用多余电能带动EC制冷, 其余时段则由AC制冷。

本文研究了计及IDR的IES的优化调度问题, 分别构建了以系统运行成本最小和系统碳排放成本最小的目标函数, 对SMA进行改进得到ISMA, 通过算例仿真结果得出以下结论。

(1) 考虑IDR时, 电、热、冷负荷峰谷差分别降低了22.1%, 19.3%, 17.9%, 提高了系统经济性。

(2) 选用标准测试函数对比分析ISMA与PSO算法、GWO算法、SSA、SMA的优化曲线, ISMA算法在优化速度和精度上具有明显优势。

(3) 对比分析ISMA与PSO算法、GWO算法、SSA、SMA分别用于IES调度时的系统运行成本和碳排放成本, 并分析ISMA参与调度的电、热、冷负荷平衡分布情况, 发现ISMA在保证系统稳定运行的同时, 能够有效降低运行成本与碳排放成本。这可以为IES的实际运行调度提供参考。