随着全球能源危机的加剧, 包含太阳能、风能等新能源, 并整合了电、气、热等多元产、用能源形式的综合能源系统, 凭借其节能、环保和灵活等特点受到了世界各国的广泛关注. 近年来, 国内外学者针对有关综合能源系统的各种关键理论与工程技术, 如优化调度^[1]、发电预测^[2]、协同控制^[3]等, 开展了大量的科学研究和理论分析工作并得到了丰硕的研究成果.

与已有电、气、热等单一供能系统一样, 在综合能源系统的运行过程中, 如何实现其优化调度是人们最为关注的问题之一, 即如何在满足各机组单元运行约束的前提下, 通过优化分配负荷需求并合理安排产能计划实现系统运行总成本(经济成本、环境成本等)最低^[4]. 近年来, 针对综合能源系统的优化调度问题, 国内外科研学者提出了许多成熟的解决办法, 其总体上可以分为两类, 即集中式方法和分布式方法. 其中, 集中式方法主要包括解析式算法和启发式算法, 如迭代法、牛顿法和遗传算法等. 1)解析式算法, 如: 文献[5]提出了一种混合整数优化方法解决具有多不确定性的综合能源系统调度问题; 2)启发式算法, 如: 文献[6]提出了一种基于直接搜索方法(Direct search method, DSM)的解决综合能源系统经济调度问题的新方法, 该算法可以处理发电机的非线性特性所带来的问题. 集中式算法在获得最优解方面具有一定优势, 但存在单点故障敏感、通信负担较大和隐私泄露等缺点. 较之于集中式方法, 分布式方法可以利用稀疏的通信网络结构实现网络内各组件的分散式协作, 可有效提高系统的鲁棒性和灵活性, 并兼具保护隐私等优点. 因此, 近年来基于分布式方法研究综合能源系统的优化调度问题, 已经成为国内外学术界的研究热点. 文献[7]提出了一种智能能源枢纽的分布式综合需求响应算法, 但它需要一个中央价格协调器和一个集中的通信网络来更新当地的能源需求信息, 因而认为该算法是一种非完全分布式的方法; 文献[8]针对多能源系统的优化调度问题提出了一种基于一致性的分布式方法; 文献[9]提出了一种针对多能源需求的两层分布式优化策略, 以优化生产者的利润和用户的舒适度. 然而, 文献[8]和文献[9]仅考虑了全局等式和局部不等式约束, 没有讨论电力线路传输约束等耦合约束; 文献[10]提出了一种综合能源系统经济运行的分布式优化方法, 以适应间歇性可再生能源发电. 但是它无法解决全局耦合的不可分离不等式约束, 所以没有考虑电网传输损耗等因素的影响. 值得指出的是, 已有的研究成果大部分都是只考虑系统运行经济性的综合能源系统单目标分布式优化调度方法, 而在兼顾系统经济性与环境友好性的多目标分布式优化调度策略的研究尚较为不足. 近期, 在多目标分布式优化理论研究方面, 文献[11]中提出了一种基于次梯度的多目标分布式优化算法, 并讨论了权重向量的选择与帕累托前沿(Pareto front)近似误差之间的关系. 但是基于次梯度的方法需要减小步长才能得到精确的解, 这可能会限制算法的性能. 文献[12]和文献[13]分别提出了基于迭代增广拉格朗日协调技术和扩散策略的多目标分布式优化算法, 但是其算法处理的约束为等式约束和线性的不等式约束, 无法处理更具一般性的约束. 针对一类带有一般性约束的非线性优化问题, 文献[14]提出了一种基于切换拓扑的多目标分布式神经动态优化方法. 该方法可以有效处理更为一般的约束, 且较之于现有多目标分布式优化方法具有可并行计算、收敛速度快和易于硬件实现的优点. 但是该算法在权向量的选取上采用了人为选定方式, 因而无法覆盖整个权向量空间, 且要求目标问题为凸优化问题. 然而, 综合能源系统的多目标优化调度问题具有耦合性强和约束一般性强等特点, 并且在考虑能量传输损耗特性等因素时, 会造成目标问题具有非凸特性, 因此上述所提算法不能直接应用于求解此类问题.

针对一类综合能源系统的优化调度问题, 本文提出了一种基于分布式神经动态优化的多目标优化调度方法. 首先, 建立了此类综合能源系统兼顾系统经济性与环境友好性的多目标优化调度模型, 该模型除考虑了系统能量平衡和设备特性约束等一般性约束外, 还同时考虑了系统能量传输损耗与网络传输约束, 因而不同于大部分已有分布式优化方法, 本文需要解决的是一类非凸的多目标优化问题; 其次, 针对此类问题, 本文提出了一种基于动态权重的分布式神经动态优化算法, 该算法不仅可以有效处理文献[10]和文献[15]中无法解决的全局耦合不可分离不等式约束, 进而解决能量传输损耗因素带来的非凸问题, 而且可在动态权重系统的指导下生成不断变化的权向量, 使得智能体输出轨迹可涵盖整个帕累托前沿, 进而有效解决此类非凸多目标优化问题; 此外, 该算法仅要求每个智能体与自己邻居节点交互部分信息来计算本地最优解, 具有较高的灵活性和隐私性等优势. 最后, 本文搭建了15节点的综合能源MATLAB仿真测试系统, 通过2个仿真算例, 验证了算法的正确性和有效性.

1.   电−气−热综合能源系统

本文考虑的综合能源系统如图1所示, 其主要设备包括: 1)仅发电装置, 包括常规发电设备(Conventional generator, CG)、分布式可再生发电设备(Distributed renewable generator, DRG)、燃料电机(Fuel generator, FG)和分布式电力储能设备(Distributed power storage device, DPSD); 2)仅发热装置, 包括分布式可再生制热装置(Distributed renewable heating device, DRHD)、燃料制热装置(Fuel heating device, FHD)和分布式蓄热装置(Distributed heat storage device, DHSD); 3)热电联产装置(Combined heat and power device, CHP); 4)燃气供应商(Gas producer, GP). 同时, 负荷(Energy load, EL)可分为3类, 即电负荷(Power load, PL)、热负荷(Heat load, HL)和气负荷(Gas load, GL), 其中每类负荷都包含常规负荷和柔性负荷. 此类综合能源系统可应用于现代工业园区和智能楼宇等场所.

图 1  本文所考虑的综合能源系统结构图

Fig. 1  Architecture of integrated energy system considered in this paper

下载: 全尺寸图片 幻灯片

2.   多目标优化调度模型

综合能源系统多目标优化调度旨在满足能量平衡约束和各机组单元运行约束的前提下, 实现系统兼顾经济性与环境友好性的综合优化. 本文以图1所示的综合能源系统为研究对象, 为实现其综合成本最优, 以下给出了此多目标优化调度问题的数学模型.

2.1   目标函数

2.1.1   经济性目标

考虑经济性目标是使系统的运行总成本最低, 因此本文的经济性目标函数具体描述为

式中, ${F_{{\rm{fuel}}}}$表示综合能源系统的运行成本(单位: $) $C( \cdot )$表示综合能源系统中各能源参与者的运行成本; ${\vartheta _{{\rm{CG}}}}$, ${\vartheta _{{\rm{DRG}}}}$, ${\vartheta _{{\rm{FG}}}}$, ${\vartheta _{{\rm{DPSD}}}}$, ${\vartheta _{{\rm{DRHD}}}}$, ${\vartheta _{{\rm{FHD}}}}$, ${\vartheta _{{\rm{DHSD}}}}$,${\vartheta _{{\rm{CHP}}}}$, ${\vartheta _{{\rm{GP}}}}$, ${\vartheta _{{\rm{PL}}}}$, ${\vartheta _{{\rm{HL}}}}$, ${\vartheta _{{\rm{GL}}}}$表示CG, DRG, FG, DPSD, DRHD, FHD, DHSD, CHP, GP, PL, HL, GL集合, $\vartheta $是上述集合的并集.

2.1.1.1   常规发电设备运行成本

定义每个常规发电设备(CG)的成本函数, 如式(2)所示.

式中, $\alpha _i^{{\rm{CG}}}$, $\beta _i^{{\rm{CG}}}$和$\gamma _i^{{\rm{CG}}}$表示常规发电设备$i$的成本系数; $P_{i,t}^{{\rm{CG}}}$是由常规发电设备$i$在$t$时刻产生的功率(单位: MW).

2.1.1.2   分布式可再生发电设备运行成本

对于分布式可再生发电设备(DRG)来说, 主要能源是燃料成本为零的太阳能和风能. 由于间歇性和随机性的特点, 可再生发电设备一般被视为不可调度的单元. 为了在优化调度中考虑这些问题, 本文参考文献[16]中对可再生发电设备的成本函数建模机理, 将DRG的成本函数建模为

式中, $b_i^{{\rm{RG}}} > 0$, $\varepsilon _i^{{\rm{RG}}} > 0$为可再生发电设备$i$的成本系数, $\gamma _i^{{\rm{RG}}} < 0$为惩罚因子; $P_{i,t}^{{\rm{RG}}}$是由可再生发电设备$i$在$t$时刻产生的功率(单位: MW); $P_{i,t}^{{\rm{RG,min}}}$和$P_{i,t}^{{\rm{RG,max}}}$是可再生发电设备的发电功率的下限和上限(单位: MW). 此外, DRG模型也可应用于DRHD. 并设$H_{i,t}^{{\rm{RH}}}$和$C(H_{i,t}^{{\rm{RH}}})$分别表示相应的热功率和成本函数.

2.1.1.3   燃料电机运行成本

定义每个燃料电机(FG)的成本函数, 如式(4)所示.

式中, $a_i^{{\rm{FG}}}$, $b_i^{{\rm{FG}}}$, $c_i^{{\rm{FG}}}$, $\varepsilon _i^{{\rm{FG}}}$和$\eta _i^{{\rm{FG}}}$表示非负成本系数. $P_{i,t}^{{\rm{FG}}}$是第$i$个FG在$t$时刻的输出功率(单位: MW); 此外, 如果上述FG是燃气轮机, 即消耗天然气产生电能, 则发电量与耗气量之间的关系可以用下式近似计算

式中, $\theta $是单位MW与SCM/h的转化率, 值约84. $\eta _i^{{\rm{FG}}}$是燃气轮机的效率, 本文取0.8.

FG模型也可应用于FHD. 并设$H_{i,t}^{{\rm{FH}}}$和$C(H_{i,t}^{{\rm{FH}}})$分别表示相应的热功率和成本函数. 如果加热装置是燃气锅炉, 则可以用以下方法近似估算气体消耗量:

式中, $\eta _i^{{\rm{FH}}}$是燃气锅炉的效率, 本文取0.8.

2.1.1.4   燃气供应商运行成本

将燃气供应商(GP)的成本函数建模为

式中, $a_i^{{\rm{GP}}}$, $b_i^{{\rm{GP}}}$, $c_i^{{\rm{GP}}}$和$d_i^{{\rm{GP}}}$表示燃气供应商$i$的成本系数; $g_{i,t}^{{\rm{GP}}}$是由燃气供应商$i$在$t$时刻产生的功率(单位: SCM/h).

2.1.1.5   热电联产机组运行成本

将热电联产机组的成本函数建模为以下凸函数:

式中, $a_i^{{\rm{CHP}}}$, $b_i^{{\rm{CHP}}}$, $\alpha _i^{{\rm{CHP}}}$, $\beta _i^{{\rm{CHP}}}$, $\sigma _i^{{\rm{CHP}}}$和$c_i^{{\rm{CHP}}}$为热电联产机组的成本系数, $P_{i,t}^{{\rm{CHP}}}$和$H_{i,t}^{{\rm{CHP}}}$是$i$个CHP在$t$时刻的电功率和热功率(单位: MW). 如果CHP的供给是天然气, 那么总耗气量可以用下式计算

式中, $\eta _i^{{\rm{CHP}}}$是CHP的效率, 本文取0.88.

2.1.1.6   电力储能设备运行成本

电力储能设备不能同时充放电. 为了方便起见, 将$t$时刻的第 $i$个电池的充放电量定义为$P_{i,t}^{{\rm{PS}}}$, (单位: MW), 其中负值是充电过程, 而正值是放电过程. 将储能设备的成本函数建模为

式中, $a_i^{{\rm{PS}}}$和$b_i^{{\rm{PS}}}$表示电池$i$的成本系数; 将此模型应用于蓄热装置, $H_{i,t}^{{\rm{PS}}}$表示第$i$个蓄热装置在$t$时刻的热功率, $C(H_{i,t}^{{\rm{HS}}})$表示其成本函数.

2.1.1.7   柔性负荷运行成本

将综合能源系统中的负荷分为3类, 分别是电负荷、热负荷和气负荷; 其中每一类负荷都包含常规负荷(系统必须承担的负荷, 不参与调度)和柔性负荷(系统可以灵活调节的, 参与调度). 根据文献[17]将三种柔性负荷的成本函数建模为

式中, $a_i^{{\rm{PL}}}$, $b_i^{{\rm{PL}}}$是正的柔性电负荷成本系数; $P_{i,t}^{{\rm{PL}}}$是柔性电负荷的功率(单位: MW). 此外, PL模型也可应用于HL模型和GL模型. 并设$H_{i,t}^{{\rm{HL}}}$, $g_{i,t}^{{\rm{GL}}}$, $C(H_{i,t}^{{\rm{HL}}})$, $C(g_{i,t}^{{\rm{GL}}})$分别表示相应的功率和成本函数.

2.1.2   环境友好性目标

环境友好性目标旨在实现系统污染气体总排放量最低, 通常以CO₂, SOx和NOx作为主要污染气体, 即综合最小化系统含硫、氮、碳污染气体排放总量. 本文的排放成本模型采用多项式和指数项的组合的形式, 具体描述如下^[18]:

其中,

式中, ${F_{{\rm{emission}}}}$表示综合能源系统的排放成本(单位: $); ${E_c}$是总二氧化碳排放成本(单位: $). ${\tau _i}$是其排放系数; ${E_s}$为SOx和NOx的总排放成本(单位: $). ${\omega _i}$, ${\mu _i}$, ${\kappa _i}$, ${\zeta _i}$, ${\pi _i}$是CG, FG, CHP, FHD, GP的排放系数.

由于DRG和DRHD为分布式可再生发电和制热设备, 即对环境的代价成本为0, 所以在环境目标中不考虑它们. 此外, 分布式储能和蓄热装置对环境造成的污染也为0. 为了方便起见, 本文考虑的柔性负荷为居民负荷, 不考虑其使用对环境的影响.

2.1.3   综合调度目标

在综合能源系统多目标优化调度的目标函数中, 运行成本和排放成本性质上是相互矛盾的, 两者都必须同时考虑. 因此, 必须得到最优排放调度与最优经济成本调度之间的权衡关系. 当同时考虑最优排放和最优经济运行成本时, 多目标优化问题的目标函数可以表示为^[18]

式中, $\min [{F_{{\rm{fuel}}}},{F_{{\rm{emission}}}}]$表示对${F_{{\rm{fuel}}}}$和${F_{{\rm{emission}}}}$两个相互矛盾的目标进行协调和折中处理, 同时得到二者的权衡最小化. 如式(41)所示, 通过加权和法实现二者的权衡最小化, 而不是取${F_{{\rm{fuel}}}}$和${F_{{\rm{emission}}}}$中较小的一个为目标函数.

2.2   约束条件

2.2.1   系统功率平衡约束

系统运行时需满足电功率平衡等式约束、热功率平衡等式约束和气功率平衡等式约束, 具体描述为

式中, ${P_{{\rm{load}}}}$, ${H_{{\rm{load}}}}$, ${g_{{\rm{load}}}}$分别表示相应的电力、热力和燃气的常规负荷(单位: MW). ${P_{{\rm{loss}}}}$表示输电损耗(单位: MW), 采用正定$B$矩阵损耗计量式近似计算为^[19]

式中, ${B_{im}}$, ${B_{ij}}$, ${B_{jn}}$表示损耗系数矩阵B中所对应的元素. ${N_P}$和${N_c}$表示纯发电机(CG, DRG, FG)和CHP的数量.

${H_{{\rm{loss}}}}$表示传热损失, 根据稳态传热基本原理^[19]近似计算为

式中, $n$表示热媒流经管道的总段数; ${l_g}$表示热媒流经管道$g$的长度(单位: km); ${t_{s,f}}$表示热网节点$f$的供水温度(单位: ℃), $f = 1,2,3, \cdots ,({N_c} + {N_h})$; ${N_h}$表示纯制热设备的总数(DRHD, FHD); ${t_{a,g}}$表示管道$g$中周围介质的平均温度(单位: ℃); ${R_h}$表示热媒到周围介质间每千米管道的总热阻(单位: km·℃/kW).

2.2.2   常规发电设备出力上下限约束

系统运行时常规发电设备需满足出力上下限不等式约束, 具体描述为

式中, $P_{i,t}^{{\rm{cg,min}}}$和$P_{i,t}^{{\rm{cg,max}}}$是常规发电设备的发电功率的下限和上限(单位: MW). $P_{i,t}^{{\rm{cg,ramp}}}$是常规发电设备的爬坡约束(单位: MW). 值得注意的是, 在求解t时刻CG的输出功率$P_{i,t}^{{\rm{CG}}}$时, 爬坡约束要求其满足关于$t-1 $时刻的CG输出功率$P_{i,t{\rm{ - }}1}^{{\rm{CG}}}$(定值, 上一时刻已经获得)的不等式要求.

2.2.3   分布式可再生发电设备出力上下限约束

系统运行时分布式可再生发电设备需满足出力上下限不等式约束, 具体描述为

式中, $P_{i,t}^{{\rm{RG,min}}}$和$P_{i,t}^{{\rm{RG,max}}}$是可再生发电设备的发电功率的下限和上限(单位: MW). 此外, DRG的约束模型也可应用于DRHD.

2.2.4   燃料电机出力上下限约束

系统运行时燃料电机需满足出力上下限不等式约束, 具体描述为

式中, $P_{i,t}^{{\rm{FG}},\min }$和$P_{i,t}^{{\rm{FG}},\max }$是燃料电机的发电功率的下限和上限(单位: MW). $P_{i,t}^{{\rm{FG}},{\rm{ramp}}}$是燃料电机的爬坡约束(单位: MW). 此外, FG的约束模型也可应用于FHD.

2.2.5   燃气供应商出力上下限约束

系统运行时燃气供应商需满足出力上下限不等式约束, 具体描述为

式中, $g_{i,t}^{{\rm{GP}},\min }$和$g_{i,t}^{{\rm{GP}},\max }$是燃气供应商的功率的下限和上限(单位: SCM/h).

2.2.6   热电联产机组热 – 电可运行域约束

系统运行时CHP机组需满足热–电可运行域约束, 具体描述为

式中, $e_{i,m}^{{\rm{CHP}}}$, $f_{i,m}^{{\rm{CHP}}}$, $z_{i,m}^{{\rm{CHP}}}$表示第$i$个CHP机组的第$m$个不等式的系数. $P_{i,t}^{{\rm{CHP}},{\rm{ramp}}}$是CHP机组的爬坡约束(单位: MW).

2.2.7   电力储能设备充放电约束

系统运行时电力储能设备需满足充放电功率约束和存储能量约束, 具体描述为

式中, $P_{i,t}^{{\rm{CH}},\max}$和$P_{i,t}^{{\rm{DS}},\max}$是电力储能设备最大充放电速率(单位: MW); $S_{i,t}^{{\rm{PS}}}$是电力储能设备到$t$时刻已经存储的能量(单位: MW), 也要满足容量约束; $S_{i,t}^{{\rm{PS}},\min }$和$S_{i,t}^{{\rm{PS}},\max}$是容量约束的上下界(单位: MW); $\Delta t$是一个充电周期, 这里我们设定为1 h. 设置电力储能设备的初始储存电量为$S_{i,0}^{{\rm{PS}}}$(单位: MW). 在一个充电周期(24 h)结束后, 电力储能设备的储存电量应保持$S_{i,0}^{{\rm{PS}}}$不变. 此外, 电力储能的约束模型可应用于蓄热装置. 值得注意的是, 在式(31)中, 在求解t时刻电力储能设备的输出功率$P_{i,t}^{{\rm{PS}}}$时, $S_{i,t{\rm{ - }}1}^{{\rm{PS}}}$为$t-1 $时刻电力储能设备充放电后剩余电量, 为一个定值.

2.2.8   柔性负荷需求

系统运行时柔性负荷需满足功率约束, 具体描述为

式中, ${P_{{\rm{load}}}}$是常规电负荷(单位: MW). $P_{i,t}^{{\rm{PL}},\max}$是柔性电负荷的功率上限(单位: MW). 此外, PL约束模型也可应用于HL模型和GL模型.

2.2.9   网络拓扑支路约束

网络拓扑支路约束是影响综合能源系统优化结果的重要因素^[20], 本文对电网、热网和气网的网络拓扑支路约束的具体考量描述如下.

2.2.9.1   电网支路约束

系统运行时电网传输线路需满足支路约束, 具体描述为^[21]

式中, ${P_e}$表示电网线路e的传输功率(单位: MW); $P_e^{\max }$和$P_e^{\min }$分别表示电网线路e的传输功率上限和下限(单位: MW).

2.2.9.2   热网支路约束

热网输送管道在系统运行时也需要满足支路约束, 具体描述为^[19]

式中, ${t_{s,f}}$表示热源在节点$f$处的供应温度(单位: ℃); $t_{s,f}^{\max }$, $t_{s,f}^{\min }$分别表示热源在节点f处的供应温度的上下限(单位: ℃); ${m_g}$为区域供热管网管道g的质量流量(单位: m³/h); $m_g^{\max }$和$m_g^{\min }$表示管道$g$的质量流量允许上下限(单位: m³/h); ${q_f}$为热网节点$f$所传递的热量 (单位: MW); ${t_{r,f}}$为热源在节点f的返回温度(单位: ℃); c为加热介质的比热容(单位: kJ/kg·℃).

2.2.9.3   气网支路约束

气网输送管道流量由管道的特性(如长度、直径、工作温度和关联节点之间的压差)决定. 该流量可以建模为^[20]

式中, ${g_{ij}}$表示从节点i到节点j的天然气流量(单位: SCM/h); ${C_{ij}}$表示$i-j $段天然气管道的管道常数; ${\pi _i}$表示节点i的天然气压力(单位: Pa); ${\pi _j}$表示节点j的天然气压力(单位: Pa); $\pi _i^{\min }$, $\pi _i^{\max }$分别表示节点i处天然气压力的上下限(单位: Pa). 值得注意的是, 约束(37)为非线性的等式约束, 为非凸约束, 且式中的变量均为 $h(x) = {x^2}$ 的形式. 根据文献[20], 采用分段近似的方法对其进行线性化, 将非线性等式约束转化为线性等式约束, 由此约束(37)转换为凸约束, 适用于凸优化的求解.

注 1. 本文提出了考虑能源的产生、转换和消费过程中电、热、气三者之间的强耦合关系的综合能源系统多目标优化调度模型. 然而, 现有的文献大多是关于电热系统, 涉及气网的文献较少. 相较之于电热系统, 考虑气网后, 由于燃气可通过FG、CHP和FHD等装置转化成电能和热能, 需要在目标模型中增加相应设备的经济成本和环境成本函数, 此外还需要进一步引入气网平衡方程和气网输送管道流量约束, 因此系统模型的复杂度和耦合性均明显增强. 另外, 在设计算法时, 也需对应引入新增相关设备(FG和FHD)的成本决策变量、对应气网平衡方程的全局耦合等式约束(式(18))和对应设备出力约束、天然气管道压力约束等若干不等式约束. 因而, 要求算法具有较强的解耦能力和计算性能. 为了便于成本函数模型的建立, 现有文献中成本函数大都考虑二次函数形式, 导致了电热气综合能源系统模型的不准确性. 因此, 我们考虑在目标函数时加入了非线性项, 建立了一个相对准确的综合能源系统模型. 同时考虑了系统的需求响应策略, 使系统的运行更加灵活、经济, 更接近于实际的系统运行.

3.   分布式算法

3.1   图论

用$G = (V,E,A)$来表示一个图. 图中节点集$V = \{ {v_i}|i = 1,2, \cdots ,n\} $是有限的非空集, 边集表示为$E \subseteq V \times V.$ 边$({v_i},{v_j})$表示节点$i$与节点$j$相互连接. ${\boldsymbol{A}} = [{a_{ij}}] \in {{\bf{R}}^{n \times n}}$表示图 $G$的加权邻接矩阵. 加权邻接矩阵$A$多用来表示点与边的关系, 其中对角元素恒为0, 若非对角元素${a_{ij}} > 0$, 则$({v_i},{v_j}) \in E$. 若非对角元素${a_{ij}}=0$, 则$({v_i},{v_j}) \notin E$. 在无向图中, 因为节点间的连边不具有方向性, 所以 $({v_i},{v_j}) \in E \Leftrightarrow ({v_j},{v_i}) \in E$. 图中节点${v_i}$和${v_j}$之间的路径是由边$({v_i},{v_{i1}}),({v_{i1}},{v_{i2}}), \cdots ,({v_{ik}},{v_{ij}})$构成的. 如果在任意一对不同的节点${v_i}$和${v_j}$之间都存在一条路径$(i,j = 1,2, \cdots ,n)$, 则无向图$G$是连通的.

假设 1. 图$G$是无向且连通的.

定义节点${v_i}$的出度为 $deg({v_i}) = \sum\nolimits_{j = 1,j \ne i}^n {{a_{ij}}} $$(i,j = 1,2, \cdots ,n)$; 定义图的Laplace矩阵为L = ${\boldsymbol{ D}} - {\boldsymbol{A}}$, 其中对角矩阵${\boldsymbol{D}} = {\rm{diag}}\{ deg({v_1}), \cdots ,deg({v_n})\}$定义为图的出度矩阵. 若假设1成立, 则${\boldsymbol{L}}$是对称且半正定的矩阵.

3.2   问题的转换

将综合能源系统多目标优化调度问题抽象成一个多目标优化问题(Multi-objective optimization problem, MOOP), 其描述为

式中, ${\boldsymbol{x}} \in {{\bf{R}}^n}$是决策变量的向量, ${\boldsymbol{F}}({\boldsymbol{x}}) \in {{\bf{R}}^2}$是目标函数向量, 其中, ${F_i}({\boldsymbol{x}}):{{\bf{R}}^n} \to {\bf{R}}\,(i = 1,2)$. $h({\boldsymbol{x}}):$${{\bf{R}}^n} \to {{\bf{R}}^m}$表示不等式约束. ${\boldsymbol{A }}\in {{\bf{R}}^{p \times n}}$是行满秩矩阵. 可行域$X$定义为 $\{ {\boldsymbol{x}} \in {{\bf{R}}^n}:h({\boldsymbol{x}}) \leq 0,{\boldsymbol{Ax}} = {\boldsymbol{b}}\}$. 可行目标区域$Z$定义为 $\{ {\boldsymbol{F}}({\boldsymbol{x}}):{\boldsymbol{x}} \in X\}$集合.

引理1^[22-23]. 如果所有的目标函数是凸函数, 并且可行域$X$是凸集, 则MOOP式(40)是凸优化问题.

与单一目标的优化问题不同, 由于目标函数之间的矛盾, 多目标优化往往不可能得到一个能同时最小化所有目标的单一解. 因此, MOOP式(40)的解被描述为帕累托最优, 它反映了所有目标函数之间的权衡. 帕累托的最优性定义如下.

定义1. 如果不存在${\boldsymbol{x}} \in X$使得${F_i}({\boldsymbol{x}}) \leq {F_i}({{\boldsymbol{x}}^ * })$, $\forall i \in \{ 1,2\} $和${F_j}({\boldsymbol{x}}) < {F_j}({{\boldsymbol{x}}^ * })$, $\exists i \in \{ 1,2\} $成立, 则决策变量${{\boldsymbol{x}}^ * } \in X$称为帕累托最优解集

帕累托最优解集${{\boldsymbol{x}}^ * } \in X$, 经目标函数映射构成了该问题的帕累托最优前沿${{{{\boldsymbol{z}}}}^ * } \in{\bf{ Z}}$.

用$P$表示问题(40)的帕累托最优解集. 多目标优化的目标是找到一组尽可能接近帕累托最优解集的解. 此外, 集合应该尽可能多样化, 以便它能够很好地逼近帕累托最优前沿.

本文采用线性加权和法将多目标优化问题转化为单目标优化问题, 即每个目标乘以一个加权因子, 将MOOP式(40)转换为

式中, ${w_1}$和${w_2}$是对应目标函数的权重因子. 定义一个空间${S^2} = \{ w \in {{\bf{R}}^2}:0 \leq {w_i} \leq 1,{\left\| w \right\|_1} = 1\} $. 定义一个集合${P_s}$, 它是由存在$w \in {S^2}$, 使得${\boldsymbol{x}}$是问题(41)的最优解组成的. 接下来, 由引理2说明${P_s}$和$P$之间的关系.

引理2^[24]. 对于任意MOOP式(40), 假设其是凸优化问题, 那么$P = {P_s}$.

综合能源系统多目标优化调度问题的2个目标函数${F_{{\rm{fuel}}}}$和${F_{{\rm{emission}}}}$均为凸函数, 且不等式约束都为凸约束, 但由于式(19)中电损的形式为二次函数形式, 导致电平衡等式约束(16)为非凸等式约束. 因而根据本文提出的定理1, 本文将电平衡等式约束(16)转化为不等式约束(43b). 计算可知, 不等式约束(43b)的海森矩阵为正定矩阵, 即不等式约束(43b)为凸约束. 这样就得到了一个与原始问题等价的新凸优化问题(43a), 并且我们提出了一个充分条件(42), 并证明了在此条件下优化问题(43a)与原始问题有相同的最优解.

定理1. 假设所有发电设备的功率下界和所有电负荷的功率上界都满足:

将等式约束(16)转化为不等式约束,从而得到新的凸优化问题(43a)为

式中, 对于不同的设备, $m$和$M$分别代表相应的出力下界和上界. $\vartheta $是包含所有设备的集合.

本文参考文献[17]通过分析凸优化问题(43a)的KKT (Karush-Kuhn-Tucker)条件, 从而证明在条件(42)下原非凸优化问题与凸优化问题(43a)的等价性. 首先, 假设问题(43a)中的约束(43b)取“>”, 通过分析KKT条件, 得到了与条件(42)相反的结果, 从而证明问题(43a)在约束(43b)中取“=”时才能得到最优解, 即原非凸优化问题和凸优化问题(43a)具有相同的最优解. 具体证明参见附录A.

4.   基于神经动力学的方法

MOOP的神经动力学方法由神经动力学系统和动态权值更新两部分组成.

4.1   神经网络模型

为了建立求解多目标优化问题的神经网络模型, 我们将神经网络模型应用于具有单一目标函数的转化问题(43a). 接下来提出神经网络模型.

当给定权重向量时, 通过加权法, 将上述问题多目标优化问题转化为单目标优化问题(43a). 我们将上述优化问题抽象成一般的数学形式, 如下所示, 考虑用一个具有m个Agent (其中一个Agent代表一个递归神经网络(Recurrent neural network, RNN)的分布式神经动力网络求解该优化问题

式中, ${\boldsymbol{x}} \in {{\bf{R}}^{{n}}}$, ${f_i}({\boldsymbol{x}}) \in {{\bf{R}}^{{n}}} \to {\bf{R}}$是第$i$个RNN对应的目标函数, ${f_i}({\boldsymbol{x}}) = {w_1}C({\boldsymbol{x}}) + {w_2}({E_s}({\boldsymbol{x}}) + {E_c}({\boldsymbol{x}}))$. 且根据成本函数和排放函数的形式, 可知${f_i}$是凸函数. ${A_i} \in {{\bf{R}}^{{{{r}}_{{i}}} \times {{n}}}}$, ${b_i} \in {{\bf{R}}^{{{{r}}_i}}}\,(0 \leq {r_i} \leq n)$是第$i$个RNN对应的等式约束, ${h_i}({\boldsymbol{x}}):{{\bf{R}}^{{n}}} \to {{\bf{R}}^{{{{l}}_{{i}}}}}$是第$i$个RNN对应的不等式约束.

假设2. 在${{\bf{R}}^{{n}}}$内, ${f_i}({\boldsymbol{x}})$和${h_i}({\boldsymbol{x}})$是凸函数, 且${A_i}$是行满秩的.

假设3 (Slater condition). 对于该优化问题, 存在${\boldsymbol{x}} \in {{\bf{R}}^{{n}}}$, 使得${A_i}{\boldsymbol{x}} = {b_i},{h_i}({\boldsymbol{x}}) \leq 0,\;i \in \{ 1,2, \cdots ,$m}成立.

基于RNN的神经网络的分布式算法是将每个RNN设置成一个Agent, 即第$i$个RNN可以访问目标函数${f_i}({\boldsymbol{x}})$、等式约束${A_i}{\boldsymbol{x}} = {b_i}、$不等式约束${h_i}({\boldsymbol{x}}) \leq 0$, 且这些信息其他RNN不知道. RNN之间的连接由连通的无向图$G$来描述.

设问题(44)的解是${{{x}}_i} \in {{\bf{R}}^{{n}}}$$(i \in \{ 1,2, \cdots ,m\} )$表示第$i$个RNN的输出向量; 设${\boldsymbol{x}} \in {{\bf{R}}^{{{mn}}}}$, 它是将m个${x_i}$相互叠加得到的一个列向量, 即${\boldsymbol{x}}$是由${x_1},{x_2}, \cdots ,{x_m}$依次从上到下组成的一份新的mn维的列向量.

设${L_m} \in {{\bf{R}}^{{{m}} \times {{m}}}}$是图$G$的Laplacian矩阵, ${I_n}$是n维单位矩阵. 则$L = {L_m} \otimes {I_n} \in {{\bf{R}}^{{{mn}} \times {{mn}}}}$, 其中$ \otimes $是Kronecker积. 与文献[25]的证明类似, 推导以下定理.

引理3^[25]. 当假设1成立时, 优化问题(44)等价于下列优化问题:

式中, ${\boldsymbol{A}}$是由${A_1},{A_2}, \cdots ,{A_m}$组成的块对角矩阵, 即${\boldsymbol{A}} = {\rm{blkdiag}}\{ {A_1},{A_2}, \cdots ,{A_m}\},$ ${\boldsymbol{b}} = {(b_1^{\rm{T}},b_2^{\rm{T}}, \cdots ,}$$b_m^{\rm{T}})^{\rm{T}} $, ${\boldsymbol{h}}({\boldsymbol{x}}) = {(h_1^{\rm{T}}{({x_1})},\;h_2^{\rm{T}}{({x_2})}, \cdots ,\;h_m^{\rm{T}}{({x_m})})^{\rm{T}}}$. 令$r = \sum\nolimits_{i = 1}^m {{r_i}} $, $l = \sum\nolimits_{i = 1}^m {{l_i}} $, 则 ${\boldsymbol{A}} \in {{\bf{R}}^{{{r}} \times {{mn}}}}$, ${\boldsymbol{b}} \in {{\bf{R}}^{{r}}}$, ${\boldsymbol{h}}({\boldsymbol{x}}):{{\bf{R}}^{{{mn}}}} \to {{\bf{R}}^{{l}}}$.

为了解决上述问题, 在文献[25]的启发下, 提出一种新的分布式神经动态优化算法, 具体为

式中, ${\boldsymbol{x}} \in {{\bf{R}}^{{{mn}}}},{\boldsymbol{\lambda}} \in {{\bf{R}}^{{l}}},{\boldsymbol{\alpha}} \in {{\bf{R}}^{{{mn}}}}$, ${\boldsymbol{I}}$是单位矩阵, ${\delta ^ + } \;=\; (\delta _1^ + ,\;\delta _2^ + , \cdots ,\delta _l^ + ),\;\delta _i^ + = \max \{ 0,{\delta _i}\}.$ 设 ${\boldsymbol{P}} = $${{\boldsymbol{A}}^{\rm{T}}}{({\boldsymbol{A}}{{\boldsymbol{A}}^{\rm{T}}})^{ - 1}}{\boldsymbol{A}} $和${\boldsymbol{q}} = {{\boldsymbol{A}}^{\rm{T}}}{({\boldsymbol{A}}{{\boldsymbol{A}}^{\rm{T}}})^{ - 1}}{\boldsymbol{b}}$. ${\boldsymbol{f}}({\boldsymbol{x}})$的梯度用$\partial {\boldsymbol{f}}({\boldsymbol{x}})$表示, $\partial h{}_i({x_i})$是$h{}_i({x_i})$的雅克比矩阵, 则$\partial {\boldsymbol{h}}({\boldsymbol{x}}) = {\rm{blkdiag}}\{ \partial {h_1}({x_1}),\partial {h_2}({x_2}),\; \cdots ,\partial {h_m}({x_m})\}$. 值得注意的是, ${\boldsymbol{P}}$和${\boldsymbol{q}}$需要用到等式约束${\boldsymbol{Ax}} = {\boldsymbol{b}}$中的常数矩阵${\boldsymbol{A}}$和常数向量${\boldsymbol{b}}$的信息. 在迭代过程中, ${\boldsymbol{A}}$和${\boldsymbol{b}}$是保持不变的, 只需要提前输入到智能体中即可, 并不需要在计算时收集全局信息获得.

注2. 文献[25]提出了一种基于递归神经网络的分布式神经动态优化算法. 该算法与本文所提算法的区别在于: 首先, 文献[25]是针对单目标凸优化问题而设计的, 而本文所讨论的是多目标优化问题, 且由于进一步考虑了能量传输损耗, 导致目标问题带有非凸特性, 因而无法直接套用文献[25]的算法和结论; 其次, 文献[25]算法中含有“$\frac{{\rm{d}}{\boldsymbol{w}}}{{\rm{d}}t} =$ x”的设计, 但优化问题的最优解${{\boldsymbol{x}}^ * }$的某些分量显然不为0, 所以${\boldsymbol{w}}$对应分量的导数不为0. 因此当$t \to \infty $时, ${\boldsymbol{w}}$的对应分量会趋于无穷大, 为算法的实现带来困难. 而本文采用了$\frac{{\rm{d}}{\boldsymbol{\alpha}}}{{\rm{d}}t} = {\boldsymbol{Lx}}$解决了这一弊端, 即当获得优化问题的最优解时, ${\boldsymbol{\alpha }}$的变化率为0.

引理4^[26]. 假设${\boldsymbol{A}} \in {{\bf{R}}^{{{r}} \times {{mn}}}}$是行满秩. 对于任意${\boldsymbol{x}} \in {{\bf{R}}^{{{mn}}}}$, ${\boldsymbol{Ax}} = {\boldsymbol{b}}$和${\boldsymbol{Px}} = {\boldsymbol{q}}$是等价的.

${\boldsymbol{u}}$是分布式PI协议, 描述为

式中, ${k_P}$和${k_I}$分别是比例和积分的增益, ${\boldsymbol{A}} = {[{a_{ij}}]_{m \times m}}$是图$G$的邻接矩阵.

假设4. 每个${u_i}(t)$的积分部分的初始状态为零.

在多智能体网络中, 每个Agent在局部约束的限制下, 被分配用来处理局部目标函数的最小值. 每个Agent根据图$G$与邻居Agent交换信息${x_i}$, 最后所有Agent共同得到最优解. 本文给出了单个Agent (即RNN)的动力学方程的离散形式, 即

式中, ${a_{ij}}$是图中$RN{N_i}$和$RN{N_j}$之间的连接权重. 图2描绘了式(49)中描述的每个RNN的框图.

图 2  用于分布式优化的$RN{N_i}$框图

Fig. 2  Block diagram of $RN{N_i}$ for distributed optimization

下载: 全尺寸图片 幻灯片

定理2. 在假设1 ~ 4下, ${{\boldsymbol{x}}^ * } \in {{\bf{R}}^{{{mn}}}}$是问题(45)的最优解当且仅当存在${{\boldsymbol{\lambda}} ^ * } \in {{\bf{R}}^{{l}}}$和${{\boldsymbol{\alpha}} ^ * } \in {{\bf{R}}^{{{mn}}}}$, 使得$({{\boldsymbol{x}}^ * },{{\boldsymbol{\lambda}} ^*},{{\boldsymbol{\alpha}} ^ * })$满足下列等式:

证明. 在假设2下, 优化问题(45)是凸优化问题. 根据KKT条件, ${{\boldsymbol{x}}^ * }$是优化问题(45)的解, 当且仅当存在${\boldsymbol{\mu}} \in {{\bf{R}}^{{r}}},{\boldsymbol{\lambda}} \in {{\bf{R}}^{{l}}},{\boldsymbol{\beta}} \in {{\bf{R}}^{{{mn}}}}$, 使得$({{\boldsymbol{x}}^ * },{{\boldsymbol{\mu}} ^ * },$${{\boldsymbol{\lambda}} ^ * },{{\boldsymbol{\beta}} ^ * }) $满足下列条件:

令$\overline {\boldsymbol{\alpha}} = {({1_m} \otimes I)^{\rm{T}}}{\boldsymbol{\alpha}} $, 在假设1的条件下, 有${\rm{1}}_m^{\rm{T}}L = $0. 故${\overline {\boldsymbol{\alpha}} ^\prime } = {({1_m} \otimes I)^{\rm{T}}}{\boldsymbol{Lx}} = ((1_m^{\rm{T}}L) \otimes I){\boldsymbol{x}} = 0$, 由假设4, 可知${({1_m} \otimes I)^{\rm{T}}}{\boldsymbol{\alpha}} (t) = {({1_m} \otimes I)^{\rm{T}}}{\boldsymbol{\alpha}} (0) = 0$. 由$rank({\boldsymbol{L}}) = (m - 1)n$和${({1_m} \otimes I)^{\rm{T}}}[{\boldsymbol{L}} {{\boldsymbol{\alpha}} ^*}] = 0$, 可得:

所以$rank([{\boldsymbol{L}} {{\boldsymbol{\alpha}} ^*}])=rank({\boldsymbol{L}})$. 因此存在${\boldsymbol{L}}{{\boldsymbol{\beta}} ^*} = $$ {k_I}{{\boldsymbol{\alpha}} ^*}$.

将式(51)两侧同时乘以${({\boldsymbol{A}}{{\boldsymbol{A}}^{\rm{T}}})^{ - 1}}{\boldsymbol{A}}$, 可得

再将式(56)代回式(51), 可得:

根据引理4可得:

此外, 式(53)中${\boldsymbol{h}}({{\boldsymbol{x}}^*})$和${{\boldsymbol{\lambda }}^*}$满足互补问题, 等价于下列等式:

在假设1 ~ 4下, ${{\boldsymbol{x}}^ * } \in {{\bf{R}}^{{{mn}}}}$是问题(45)的最优解当且仅当存在${{\boldsymbol{\lambda}} ^ * } \in {{\bf{R}}^{{l}}}$和${{\boldsymbol{\alpha}} ^ * } \in {{\bf{R}}^{{{mn}}}}$, 使得$({{\boldsymbol{x}}^ * },{{\boldsymbol{\lambda}} ^ * }, {{\boldsymbol{\alpha}} ^ * })$满足等式(50), 即问题(45)的最优解与RNN的平衡点一一对应. □

接下来证明提出的RNN可以收敛到平衡点$({{\boldsymbol{x}}^ * },{{\boldsymbol{\lambda}} ^ * },{{\boldsymbol{\alpha}} ^ * })$.

注3. 由式(43a)的形式可知, 其为一类全局不可分离不等式约束, 且为凸约束. 由系统信息拓扑图可知, 假如CG智能体和FG智能体并不相连, 即不可相互传递信息, 则导致CG智能体或者FG智能体无法同时获取${P_{{\rm{CG}}}}$和${P_{{\rm{FG}}}}$的信息. 现有大部分分布式算法可以有效处理加和形式的耦合约束, 如“${P_{{\rm{CG}}}} +{P_{{\rm{FG}}}}$”, 其核心思想是, 将其拆解为若干本地约束并分配给对应的智能体处理. 但是, 乘法形式的${P_{{\rm{CG}}}}\times {P_{{\rm{FG}}}}$是不可分离的, 所以现有大部分分布式算法无法有效处理此类约束. 如文献[10]和文献[15]提出的分布式神经动态优化方法(其作者也对此限制进行了说明). 在本文提到的算法中, 因为每个智能体最终输出向量是全局信息包含所有优化变量, 且与其他智能体的输出向量一致, 所以在计算时只需要调用自己和与其通信智能体的局部信息, 而不需要知道与其不进行通信的智能体的信息, 从而实现分布式.

图 4  电热气系统的信息拓扑图

Fig. 4  Information topology of the integrated electro-heating-gas system

下载: 全尺寸图片 幻灯片

图3给出了所提出的分布式方法的实现图. 可以看到, 当且仅当$RN{N_i}$和$RN{N_j}$连接时, $RN{N_i}$输出${u_i}(t)$到它的邻居$RN{N_j}$, 反之亦然.

图 3  所提出的分布式方法在IES的实现过程

Fig. 3  Implementation diagram of the proposed distributed approach for IES

下载: 全尺寸图片 幻灯片

4.2   RNN稳定性分析

为了证明所提出算法的收敛性, 提出了一种等价的算法如下:

根据之前的讨论, 在假设4的前提下, 存在${\boldsymbol{L\beta}} = {\boldsymbol{\alpha }}$, 且$\frac{ {\rm{d}}{\boldsymbol{\beta}}}{{\rm{d}}t} = {\boldsymbol{x}}$.

定义2 (输出一致性). 式(49)中的RNN神经动态网络在任何初始条件下, 都能达到输出一致性. 即

式中, $\left\| {\cdot} \right\|$是欧几里德范数.

假设5. $2{k_P} \geq {k_I}$

前面已经谈论过在假设1 ~ 4的前提下, 每个Agent达到共识的值为优化问题的最优解. 接下来基于Lyapunov函数来讨论RNN网络的稳定性.

定理3. 在假设1 ~ 5下, 对于任意给定的初始点, 输出向量${{\boldsymbol{x}}_i}$将最后达到输出一致性, 且稳定点为优化问题的最优解.

具体证明参见附录B.

用${w_1}C({\boldsymbol{x}}) + {w_2}({E_s}({\boldsymbol{x}}) + {E_c}({\boldsymbol{x}}))$代替问题(44)中的$f({\boldsymbol{x}})$, 加权因子$w$可以看作是该模型的输入. 根据前面的讨论, 可以得到下面的结果.

定理4. 当假设2成立时, 对于任意$w \in {S^2}$, 对于任何初值$({x_0},{\lambda _0},{\alpha _0})$, 神经网络(46)在Lyapunov意义下是稳定的, 且最终$x(t)$的轨迹会收敛到问题(40)的$P$中的${{\boldsymbol{x}}^ * }$.

证明. 因为$F_i$是凸函数, 所以对于 $w \in {S^2}$, $\sum\nolimits_{i = 1}^2 {{w_i}{F_i}({\boldsymbol{x}})} $也是凸函数. 根据以上的讨论, 对于任何初值$({x_0},{\lambda _0},{\alpha _0})$, 神经网络的$x(t)$的轨迹会收敛到问题(41)的一个最优解, 即${P_s}$中的一个点. 根据引理2, $P = {P_s}$. □

可以看出, 所有求解问题(41)的神经动力模型都将$w$作为参数. 事实上, $w$也可以看作是神经动力系统的输入. 如果把最终的收敛状态$x$看作是神经动力系统的输出, 那么系统定义了一个从${S^2}$到$P$的映射. 只有神经动力系统的输入包含${S^2}$中所有的权重因子, 才能得到所有的帕累托最优解. 下一节将说明权重因子的更新规则.

4.3   权重更新规则

为了获得帕累托前沿, 权重因子应该系统地更新. 现有文献多侧重于设计权值的离散时间更新规则, 获得一组近似帕累托前沿的帕累托最优解. 但是值得注意的是, 权重的选择对生成解帕累托前沿中的分布有密切的影响. 为了避免权值选取的困难, 本文设计了一个连续时间动态系统. 注意, 神经动力系统可以实现实时收敛. 现在, 考虑一个极限情况, 权重是连续变化的, 但速度相对较小. 在此基础上, 建立了具有两个时间尺度的连续神经网络模型.

采用下列动态权重:

式中, $v$表示权重的更新速度. 事实上, 系统(62)意味着$w(t) = w({t_0}) + vt$. 如果在${t_0} = 0$时, 设定$w(0) = $$ {[0,1]^{\rm{T}}}$, $v = {[1\;\; - 1]^{\rm{T}}}$, 那么$w(t)$的轨迹在时间区间[0, 1]内的轨迹覆盖了所有需要的权值, 便可以得到所有的帕累托最优解.

为简洁起见, 定义$s(t) = {({x^{\rm{T}}}(t),{\lambda ^{\rm{T}}}(t),{\alpha ^{\rm{T}}}(t))^{\rm{T}}}$, 用$\dot s = \psi (s,w)$表示神经网络(46). 然后, 两时间尺度的连续时间系统可以总结为

式中, $\varepsilon $是一个比例因子, $0 < \varepsilon \ll 1$. 可以发现式(63)是一个快速变化的系统, 式(64)是一个缓慢变化的系统. 已经证明了对于给定的权值$w$, 系统(63)是全局稳定的, ${\boldsymbol{x}}(t)$将收敛于帕累托最优解. 当$w$是时变的, 但与$s$相比速度相对较小时. 事实上, 当$\varepsilon \to 0$, 系统(63)变成为

式(65)描述了帕累托最优解${{\boldsymbol{x}}^ * }$和权重$w$之间的关系. 当$w$变化时, 我们得到了包含帕累托最优解的曲线.

如前所述, 连续时间模型较之于离散时间模型, 可以提供更为相对完整的帕累托前沿信息.

5.   仿真

案例1. 6节点电力系统、5节点区域供热系统和4节点天然气系统的综合能源系统.

本节将提出的分布式神经动态优化算法用于求解电热气耦合系统的多目标优化调度问题, 验证其有效性.

图4和图5分别给出了电热气系统的通信拓扑和物理拓扑. 设备成本参数、约束及电网传输损耗系数矩阵$B$列于附录C中(见表C1 ~ C5). 我们统一了能源规模, 即: 1 p.u. = 1 MW (电、热), 1 p.u. = 84 SCM/h (气).

图 5  电热气系统的物理拓扑图

Fig. 5  Physical topology of the integrated electro-heating-gas system

下载: 全尺寸图片 幻灯片

图5中, 节点1 ~ 12分别对应系统设备CG、DRG、FG、CHP、DPSD、PL、FHD、DRHD、DHSD、HL、GP和GL. 虚线13表示电力网络, 点划线14表示热力网络, 实线15表示天然气网络.

系统必须承担的常规负荷: 电力负荷、热负荷和天然气负荷分别是[100 100 1] (p.u.). $\varepsilon = 0.05$. 同时, 使用${S^2}$中的一系列均匀分布的权值生成21个解供参考. 权重由下式给出:

在图6中, 得到的帕累托最优解用星号表示. 可以看出, 所有的帕累托最优解都位于实线上, 即实线非常接近帕累托前沿.

图 6  综合能源系统多目标优化调度问题的帕累托前沿

Fig. 6  The Pareto front of multi-objective optimized scheduling in the integrated energy systems

下载: 全尺寸图片 幻灯片

案例2. 综合能源系统多目标优化调度.

不失一般性, 设置权重因子 ${w_1} = {w_2} = 0.5$, 权重因子可以根据人们对经济和环境两个目标的偏好进行改变.

1)收敛性评价

为了方便, 采用案例1的综合能源系统. 系统必须承担的常规负荷: 电负荷、热负荷和天然气负荷分别是[10 5 10/84] (p.u.).

图7显示了该分布式方法下所有Agent输出的轨迹, 结果表明该方法的求解轨迹在20次迭代内收敛.

图 7  常规负荷下综合能源系统各元件的最优出力

Fig. 7  Optimal outputs of components the integrated energy systems under conventional load

下载: 全尺寸图片 幻灯片

我们将“以电定热”的运行方式与应用本文算法的运行方式进行比较分析, 验证不同种类的资源经过协调利用后可以提高综合能源系统的灵活性或经济性.

采用“以电定热”的运行方式, 热电比设为0.8, 与案例2中第一节仿真结果进行对比, 得到相关设备运行功率如附录C中表C6所示, 进而计算两种不同运行方式下综合能源系统的总成本. 应用本文算法的运行方式下的系统总成本为5107.8 $, “以电定热”运行方式下的系统总成本为5352.9 $. 由此可见, 相比于传统的“以电定热”的多能协调原则, 本文算法可以进一步提高综合能源系统的经济性.

2) 24小时负荷下的多目标优化调度

测试周期为24小时, 每个周期定义为1小时. 各个周期的常规负荷分别是: 电负荷[10 20 40 40 80 300 400 400 300 400 300 250 200 150 400 400 480 200 250 250 350 400 20 10] (p.u.); 热负荷[10 10 200 300 300 200 300 200 200 200 200 200 150 150 10 200 200 300 300 400 350 300 20 10] (p.u.); 气负荷[1 1 2 2 2 2 2 2 3 3 2 2 2 2 1 2 2 1 3 3 3 3 1 1] (p.u.).

如图8所示, 在24小时负荷波动下, 根据所提出的分布式神经动力优化算法, 求解各设备在功率平衡和局部运行约束条件下的最优功率, 以达到经济和环保的目的. 调度结果同时满足机组爬坡约束和储能装置的能量约束. 同时, 通过MATLAB工具箱, 验证了分布式算法的结果与其结果一致, 验证了算法的有效性.

图 8  24小时负荷下综合能源系统各元件的最优出力

Fig. 8  Optimal outputs of components the integrated energy systems under 24 hour load

下载: 全尺寸图片 幻灯片

在该算例中, 每个时段都设定了负荷的波动, 如在1 ~ 2时段电网负荷发生变化, 由10 (p.u.)上升到20 (p.u.), 为了满足电力供需平衡, 相关发电设备应调整自身输出功率以应对电负荷波动. 从图8(a) ~ 8(c)中可以看出, 由于成本系数较高, 此时段传统产电机组CG的输出功率并没有提高, 而气网设备FG和CHP等以低成本燃气为输入的产电机组的输出功率提高了. 同理, 当热负荷发生变化时, 电网设备会根据电力负荷的实际情况并考虑成本因素, 调节CHP的运行点和新能源机组制热设备等装置来平衡热负荷的波动. 通过这一算例可以看出, 在合理的调度策略下, 综合能源系统可实现多种能源的协调优化与互补互济, 显著提高能源系统的灵活性与经济性.

3)即插即用性能测试

在这个案例研究中, 重点是测试该算法的即插即用性能. 模拟现实情况中, 考虑到可再生能源发电的不确定性, 将DRG从系统中去除, 并将与它相关的变量设置为零. 从图9中可以看出, 当DRG退出运行时, 电功率、热功率和产气功率的功率偏差逐渐收敛到零. 同时, 各Agent也收敛到最佳值. 当然, 剩余的能源参与者必须调整它们的电热气的输出/需求, 以补偿先前由断开的DRG所带来的影响, 并且最终收敛到新的解决方案. 此外, 当DRG被再次插入到系统时, 该系统再次收敛到响应新拓扑变化的解决方案, 且此时最终收敛解与断开DRG之前的收敛解相同. 这意味着该算法提供了良好的即插即用能力.

图 9  即插即用下综合能源系统各元件的最优出力

Fig. 9  Optimal outputs of components the integrated energy systems under the plug and play property

下载: 全尺寸图片 幻灯片

6.   结束语

本文针对电−热−气综合能源系统提出了一种新的分布式能源管理方法. 建立了一个相对准确的电−热−气综合能源系统模型, 提出了兼顾系统经济性与环境友好性的多目标优化调度模型. 针对该模型, 提出了一种基于可变权重的神经网络模型的神经动态优化算法, 并对其最优性和收敛性进行理论分析. 通过一个典型综合能源系统的算例仿真, 结果表明各能源参与者协调配合在经济性与环境性前提下保障系统可靠运行. 实验结果验证了所提算法的可行性与有效性. 本文算法虽然没有实现完全分布式, 但其在保证更为重要的目标函数和相应的功率约束的隐私性的前提下, 有效求解复杂的综合能源系统多目标优化调度问题. 我们未来的工作将致力于开发完全分布式算法求解综合能源系统多目标优化调度问题.

附录A.   定理1的证明

证明. 定义${\vartheta _g} = {\vartheta _{{\rm{CG}}}}\bigcup {{\vartheta _{{\rm{DRG}}}}} \bigcup {{\vartheta _{{\rm{DPSD}}}}} \bigcup {{\vartheta _{{\rm{FG}}}}} \bigcup $$ {{\vartheta _{{\rm{CHP}}}}}$. 优化问题如下:

式中, 对于不同的设备, $m$和$M$分别代表相应的出力下界和上界. 值得注意, 对于CHP机组, $m$和$M$与${H_{i,t}}$有关, 即, 当$i \in {\vartheta _{{\rm{CG}}}} \cup {\vartheta _{{\rm{DRG}}}} \cup {\vartheta _{{\rm{FG}}}}$, 则$P_i^{\min } = m \leq M = P_i^{\max }$; 当$i \in {\vartheta _{{\rm{CHP}}}}$, 则$P_i^{\min } \leq m \leq $$ M \leq P_i^{\max }$. 优化问题的拉格朗日函数为

式中, $\lambda ,\mu ,\nu ,\alpha ,\beta \geq 0$, $l \in {\vartheta _{{\rm{DRHD}}}}\bigcup {{\vartheta _{{\rm{FHD}}}}} \bigcup {{\vartheta _{{\rm{CHP}}}}} \bigcup$$ {{\vartheta _{{\rm{DHSD}}}}} $, $n \in {\vartheta _{{\rm{FG}}}}\bigcup {{\vartheta _{{\rm{FHD}}}}} \bigcup {{\vartheta _{{\rm{CHP}}}}} \bigcup {{\vartheta _{{\rm{GL}}}}} $. 为了方便, 用$f$表示${w_1}{F_{{\rm{fuel}}}} + {w_2}{F_{{\rm{emission}}}}$. 由于${F_{{\rm{fuel}}}}$和${F_{{\rm{emission}}}}$是凸函数, 所以f也是凸函数. 根据KKT条件, 可得:

接下来, 用反证法进行证明. 假设

因此, 对于条件(A2d)成立, $\lambda = 0$.

考虑一般情况, $m \ne M,i \in \vartheta $, 因为$P_{i,t}^{\min } \leq m \leq P_{i,t}^ * \leq$$ M \leq P_{i,t}^{\max } $, 所以为了保证式(A2g)和式(A2h)成立, ${\alpha _i}$和${\beta _i}$至少有一个值为0. 考虑一般情况, $\dfrac{{\partial f(P_{i,t}^ * )}}{{\partial P_{i,t}^ * }} > 0,\forall i \in {\vartheta _{{\rm{PL}}}}$. 因此${\alpha _i}= 0,{\beta _i} > 0$, $\forall i \in {\vartheta _{{\rm{PL}}}}$. 所以, 得到$P_{i,t}^ * = M,\forall i \in {\vartheta _{{\rm{PL}}}}$. 同理, $P_{i,t}^ * = m,$$\forall i \in {\vartheta _g} $. 由式(A3)可得:

注意式(A4)与充分条件(42)相矛盾. 因此, 式(A3)在充分条件下不能保证得到$P_{i,t}^ * $这一最优解. 在转换问题的第1个约束条件下, 最优解$P_{i,t}^ * $必须满足下列条件:

因此, 在充分条件下, 转换后的问题总是在第一个约束条件得到等号的情况下得到最优解. 因此, 在条件式(42)下, 原非凸优化问题和式(43a)是等价的. □

附录B.   定理3的证明

证明. 设${x^ * } \in {{\bf{R}}^{{{mn}}}}$是问题(45)的最优解. 根据定理2, 存在${{\boldsymbol{\lambda}} ^ * } \in {{\bf{R}}^{{l}}}$和${{\boldsymbol{\alpha}} ^ * } \in {{\bf{R}}^{{{mn}}}}$使得式(50)成立.

考虑以下候选Lyapunov函数:

令${{\boldsymbol{\rho}} ^*}{\rm{ = }}\partial {\boldsymbol{f}}({{\boldsymbol{x}}^{\rm{*}}})$, ${{\boldsymbol{H}}^*}{\rm{ = }}\partial {\boldsymbol{h}}({{\boldsymbol{x}}^*})$, ${{\boldsymbol{\zeta}} ^*}{\rm{ = (}}{{\boldsymbol{\lambda }}^{\rm{*}}}{\rm{ + }}{\boldsymbol{h}}{\rm{(}}{{\boldsymbol{x}}^{\rm{*}}}{\rm{)}}{{\rm{)}}^ + }$.

定义

根据链式法则, 沿着式(46)的方向对$V({\boldsymbol{x}},{\boldsymbol{\lambda}} ,{\boldsymbol{\beta }})$求导:

令${\boldsymbol{\xi}} = (I - {\boldsymbol{P}})({\boldsymbol{x}} - {\boldsymbol{\rho}} - {{\boldsymbol{H}}^{\rm{T}}}{\boldsymbol{\zeta }}- {k_P}{\boldsymbol{Lx}} - {k_I}{\boldsymbol{L\beta}} ) + {\boldsymbol{q}}$, 将上述能量函数拆解为如下形式:

对于${Q_1}$

${{\boldsymbol{x}}^*}$是优化问题(45)的最优解, 根据式(57), 也满足下列等式:

根据${\boldsymbol{\xi }}$的定义及${\boldsymbol{P}}(I - {\boldsymbol{P}}) = 0,(I - {\boldsymbol{P}}){\boldsymbol{q}} = 0$, 可得:

因为${(I - {\boldsymbol{P}})^2} = I - {\boldsymbol{P}}$, 所以

所以

对于${Q_2}$,

由于${\boldsymbol{f}}({\boldsymbol{x}})$是凸函数, 所以

所以

对于${Q_3}$, 根据式(50), 可得${{\boldsymbol{\lambda}} ^*} = {{\boldsymbol{\zeta}} ^*}$, 那么

注意:

所以

${\boldsymbol{h}}({\boldsymbol{x}})$是凸函数, 故${\boldsymbol{h}}({\boldsymbol{x}}) - {\boldsymbol{h}}({{\boldsymbol{x}}^ * }) - {\boldsymbol{H}}({\boldsymbol{x}} - {{\boldsymbol{x}}^ * }) \leq $0和${\boldsymbol{h}}({\boldsymbol{x}}) - {\boldsymbol{h}}({{\boldsymbol{x}}^ * }) - {{\boldsymbol{H}}^ * }({\boldsymbol{x}} - {{\boldsymbol{x}}^ * }) \geq 0$成立. 由于${{\boldsymbol{x}}^ * }$是问题(45)的最优解, 所以${{\boldsymbol{\zeta}} ^{\rm{T}}}{\boldsymbol{h}}({{\boldsymbol{x}}^ * }) \leq 0$. 且${{\boldsymbol{\zeta}} ^{\rm{T}}}{( -{\boldsymbol{ \lambda }}- {\boldsymbol{h}}({\boldsymbol{x}}))^ + } = 0$. 根据式(50)和式(53), ${({{\boldsymbol{\zeta }}^ * })^{\rm{T}}}{\boldsymbol{h}}({{\boldsymbol{x}}^ * }) = $${({{\boldsymbol{\lambda }}^ * })^{\rm{T}}}{\boldsymbol{h}}({{\boldsymbol{x}}^ * }) = 0 $和${({{\boldsymbol{\zeta}} ^ * })^{\rm{T}}}( - {\boldsymbol{\lambda}} - {\boldsymbol{h}}({\boldsymbol{x}}))^ + \geq 0$成立. 所以,

所以

由于$T({\boldsymbol{x}},{\boldsymbol{\lambda}} ,{\boldsymbol{\beta}} )$是凸函数, 所以

所以

由拉塞尔不变集原理,当$t \to \infty $$({\boldsymbol{x}}(t),{\boldsymbol{\lambda }}(t),{\boldsymbol{L\beta}} (t))$将收敛到$\dot V = 0$的最大不变子集, 即RNN将最终收敛到优化问题(45)的最优解$({{\boldsymbol{x}}^ * },{{\boldsymbol{\lambda}} ^*},{{\boldsymbol{\alpha}} ^ * })$. □

附录C.   各设备的相关参数描述

表 C1  各设备运行成本函数参数及出力上下限参数

Table C1  The operation cost function parameters and output limit parameters of equipment

设备名称	参数名称及参数值
CG	${\alpha _i}$	${\beta _i}$	${\gamma _i}$	$P_{{\rm{CGi}}}^{\min }$	$P_{{\rm{CGi}}}^{\max }$	$P_{{\rm{CGi}}}^{{\rm{ramp}}}$
	0.1	3	25	10	200	45
DRG	${b_i}$	${\varepsilon _i}$	${\gamma _i}$	$P_{{\rm{DRGi}}}^{\min }$	$P_{{\rm{DRGi}}}^{\max }$
	0.11	300	−1.1	84.3	103.5
DRHD	${b_i}$	${\varepsilon _i}$	${\gamma _i}$	$H_{{\rm{DRHDi}}}^{\min }$	$H_{{\rm{DRHDi}}}^{\max }$
	0.12	534	−1.3	133.2	148.2
FG	${a_i}$	${b_i}$	${c_i}$	${\varepsilon _i}$	${\eta _i}$	$P_{{\rm{FGi}}}^{\min }$	$P_{{\rm{FGi}}}^{\min }$	$P_{ {\rm{FGi} } }^{{\rm{ramp}}}$
	0.04	5	99	5	0.01	30	150	45
FHD	${a_i}$	${b_i}$	${c_i}$	${\varepsilon _i}$	${\eta _i}$	$H_{{\rm{FHDi}}}^{\min }$	$H_{{\rm{FHDi}}}^{\max }$
	0.027	5	60	5	0.008	50	150
CHP	${a_i}$	${b_i}$	${\alpha _i}$	${\beta _i}$	${\sigma _i}$	${c_i}$	$P_{ {\rm{CHPi} } }^{{\rm{ramp}}}$
	0.0345	14.5	0.03	4.2	0.031	230	45
DPSD	${a_i}$	${b_i}$	$P_{{\rm{sti}}}^{ds,\max }$	$P_{{\rm{sti}}}^{ch,\max }$	$S_{{\rm{sti}}}^{\min }$	$S_{{\rm{sti}}}^{\max }$	$S_{{\rm{sti}}}^{\rm{0}}$
	0.028	535	220	220	35	350	240
DHSD	${a_i}$	${b_i}$	$P_{{\rm{sti}}}^{ds,\max }$	$P_{{\rm{sti}}}^{ch,\max }$	$S_{{\rm{sti}}}^{\min }$	$S_{{\rm{sti}}}^{\max }$	$S_{{\rm{sti}}}^{\rm{0}}$
	0.013	961	400	400	62	620	560
GP	${a_i}$	${b_i}$	${c_i}$	${d_i}$	$g_{_{{\rm{GPi}}}}^{\min }$	$g_{_{{\rm{GPi}}}}^{\max }$
	$2\times{10}^{-6}$	0.006	50	4	100	1500
EL	$a_i^p$	$b_i^p$	$a_i^h$	$b_i^h$	$a_i^g$	$b_i^g$	$P_{{\rm{fli}}}^{\max }$	$H_{{\rm{fli}}}^{\max }$	$g_{{\rm{fli}}}^{\max }$
	0.016	42.5	0.01	25.5	0.0167	22	1000	800	500

下载: 导出CSV 

| 显示表格

CHP的可行域为

电网传输损耗系数矩阵$B$为

表 C2  各设备环境成本参数

Table C2  The environmental cost function parameters of equipment

设备名称	参数名称及参数值
CG	${\omega _i}$	${\mu _i}$	${\kappa _i}$	${\zeta _i}$	${\pi _i}$	${\tau _i}$
	0.0409	−2.7	0.649	2	0.02857	0.64
FG	${\omega _i}$	${\mu _i}$	${\kappa _i}$	${\zeta _i}$	${\pi _i}$	${\tau _i}$
	0.0254	−3.025	0.05638	5	0.0333	0.52
CHP	${\omega _i}$	${\mu _i}$	${\tau _i}$
	$ 1.5{10}^{-6} $	$ 1.5{10}^{-5} $	0.2
FHD	${\omega _i}$	${\mu _i}$	${\tau _i}$
	$ 8{10}^{-6} $	$ 1{10}^{-5} $	0.8
GP	${\omega _i}$	${\mu _i}$	${\tau _i}$
	$ 8{10}^{-6} $	$ 1{10}^{-5} $	0.6

下载: 导出CSV 

| 显示表格

表 C3  电力网络传输线路参数 (MW)

Table C3  The parameters of power network transmission pipelines (MW)

线路	$P_e^{\min }$	$P_e^{\max }$	线路	$P_e^{\min }$	$P_e^{\max }$
1 ~ 13	0	200	4 ~ 13	0	200
2 ~ 13	0	200	5 ~ 13	0	200
3 ~ 13	0	200	6 ~ 13	0	200

下载: 导出CSV 

| 显示表格

表 C4  热力网络传输管道参数 (MW)

Table C4  parameters of heating network transmission pipelines (MW)

管道	${l_g}$	$m_g^{\min }$	$m_g^{\max }$	${R_h}$	节点	$t_{s,f}^{\min }$	$t_{s,f}^{\min }$
4 ~ 14	2.8	0	2700	20	4	80	100
7 ~ 14	2.5	0	2700	20	7	80	100
8 ~ 14	3.0	0	2700	20	8	80	100

下载: 导出CSV 

| 显示表格

表 C5  气网网络传输管道参数 (MW)

Table C5  The parameters of gas network transmission pipelines (MW)

管道	${C_{ij}}$	节点	$\pi _i^{\min }$	$\pi _i^{\max }$	节点	$\pi _i^{\min }$	$\pi _i^{\max }$
11 ~ 3	10	3	1000	4000	12	1000	4000
11 ~ 4	10	4	1000	4000
11 ~ 7	10	7	1000	4000
11 ~ 12	10	12	1000	4000

下载: 导出CSV 

| 显示表格

表 C6  各设备在不同运行方式下的功率 (MW)

Table C6  Power of each device under different operating modes (MW)

	CG	DRG	DRHD	FG	FHD	CHP	DPSD	DHSD	GP	PL	HL	GL
以电定热	10	100	135	30	50	60.5/14.7	−100	−100	375.5	90	95	190
本文	12.6	103.5	133.2	30	50	55/44	−100	−100	405.9	90	95	190

下载: 导出CSV 

| 显示表格