园区综合能源微网动态优化调度算法

陈泽雄, 高军伟, 林亚培, 彭灵利, 周成鹏, 肖英豪, 杨智斌, 唐智强

(1.广东电网有限责任公司广州供电局，广东 广州 510620；2.广州市奔流电力科技有限公司,广东 广州 510670)

随着能源危机和环境污染问题的日益严重，充分利用可再生能源和提高能源利用效率成为实现人类社会可持续发展的必然要求。园区综合能源微网(integrated energy campus microgrid，IECM)技术通过多种类型能源互补协调以提高能源利用效率和可再生能源消纳容量，迅速发展起来[1-3]。IECM优化运行是在已知园区中冷、热、电、气负荷预测曲线的基础上，制订园区中各个供冷、供热、供电和供气元件的能量供应计划。由于IECM优化运行计划的制订需要满足冷、热、电、气多种子能量网的运行约束，且需要考虑运行周期中储能装置的最大状态切换次数限制等多个时段相关约束，因而需要建立IECM优化运行的日前动态优化调度模型，以获得合理的优化运行方案，提高微网运行的经济性。

对于IECM的优化运行问题，国内外学者已开展了一些研究。文献[4]在综合能源供能的智能社区背景下，结合冷热电联供(combined cooling heating and power，CCHP)系统和居民需求侧响应，提出了社区和家庭2个阶段能量优化运行模型，并通过算例验证了综合能源供应相对于单一电网供电的经济效益更好。文献[5]在综合考虑风电、光伏、燃气轮机、大电网联络线等多类型电源运行特性的基础上，以经济成本为目标，建立了冷热电联供型微网协同优化模型，并采用所提出的CoPSO-BAS算法对模型进行求解。文献[6]建立了以风能、太阳能、天然气和储能协同供能的CCHP微网优化运行模型，并采用改进粒子群优化算法求解。文献[7]对CCHP型多微网主动配电系统进行建模与优化分析，分别以配电网与CCHP微网为优化主体，建立优化运行模型。文献[8]提出一种融合需求侧虚拟储能系统的CCHP楼宇微网优化运行方法，通过在温度舒适度范围内调节楼宇室温来实现对楼宇虚拟储能系统的充放电管理。但是，上述文献中所建立的微网优化运行模型没有考虑CCHP系统中吸收式制冷机和换热机组投入台数的离散变量特性，而将其作为连续变量求解，因而得到的吸收式制冷机/换热机组的制冷/热量与实际离散投入台数对应的制冷/热量不一致，得到的优化方案难以实施。

为了提高模型的准确性，一些文献考虑了IECM优化运行中的离散变量特性来建立其混合整数规划模型。文献[9]基于CVaR理论建立了考虑供电和供热设备启停状态等离散决策变量的综合能源系统经济运行模型，采用快速粒子群优化算法和内点法求解。文献[10]提出了计及变负荷特性的小型CCHP系统经济优化模型，模型中引入季节状态变量和锅炉启停状态变量等离散决策变量，采用遗传算法求解。文献[9]和[10]都是采用人工智能算法求解混合整数非线性规划模型，计算速度较慢，不一定能满足实际应用需求。因此，部分文献将原模型中的非线性约束进行线性化以转化为较容易求解的混合整数线性规划(mixed-integer linear programming，MILP)模型来提高计算效率。文献[11]在含风电的综合能源微网中对天然气系统中的管道流量方程进行线性化，将原有非线性优化模型转化为MILP模型，提高了运算效率。文献[12]对电-气综合能源系统中的非线性方程进行分段线性化，建立了对应的MILP模型，并在模型中考虑了天然气的动态特性。上述文献对综合能源系统优化运行模型中非线性环节的线性化求解方法进行了详细分析，但是相关的优化模型中没有考虑供电网潮流方程、燃气机组输入输出功率关系和运行周期中储能装置最大状态切换次数限制等非线性环节。

鉴于此，本文建立了IECM日前动态优化调度模型，提出了求解该模型的MILP算法。通过分段线性化近似、大M法等效转换和含绝对值计算的线性转换等方法，对非线性环节进行线性化，将模型转化为容易求解的MILP模型，采用GAMS软件中的GUROBI求解器对转化后模型进行快速可靠求解。

1 IECM动态优化调度模型

1.1 IECM的结构和组成

IECM一般由供电网络、供冷/热网络、供气网络和能源站组成。能源站作为不同供能网络之间的联系枢纽，其内部含有多种能源转换设备。其中燃气发电机组通过消耗天然气向供电侧供电，其发电余热一部分被能源站中的烟气型和热水型吸收式制冷机收集并转化为冷功率，通过供冷管道中的冷水传递给用户以满足用冷需求；另一部分余热则被换热机组收集并转化为热功率，通过供热管道中的热水传递给用户以满足用热需求；当这些制冷/热功率无法满足用冷/热负荷需求时，余下的冷负荷由电制冷机补充，热负荷由燃气锅炉和电热锅炉补充。电制冷机和电热锅炉以及供冷/热管道中循环水泵所消耗的电功率都由供电侧提供。供电侧除了燃气发电机外，还装有光伏电站，同时还可向配电网购入一部分电能以满足全部用电负荷需求。可见，IECM中多种供能网络之间存在复杂的耦合关系，需要进行统一的协调调度。

1.2 目标函数

目标函数f(x)为IECM的1 d总运行费用，包括微网从天然气站购气费用、从配电网购电费用和储能装置运行费用3个部分之和，即

(1)

式中：T为运行周期的总时段数；各变量用下标t表示时段t的值；ΔT为每个时段的长度，以15 min为1个时段，则1 d包括96个时段；cng为微网购气的单价；fng,t为微网从天然气站购入的天然气流量；cpg为微网向配电网购电的单价；Ppg,t为配电网注入微网的有功功率；ces,t为储能装置运行费用；cc、ch和ce分别为储冷罐、储热罐和蓄电池运行的循环损耗费用；φcc,t、φcd,t分别为储冷罐的蓄冷、放冷功率；φhc,t、φhd,t分别为储热罐的蓄热、放热功率；Pec,t、Ped,t分别为蓄电池的充电、放电功率。

1.3 约束条件

约束条件包括供冷/热网络运行特性约束﹝式(2)[13]﹞、供气网络的运行特性约束﹝式(3)[14]﹞、供电网络的运行特性约束﹝式(4)﹞、能源站的运行特性约束﹝式(5)—(9)[15-16]﹞、储能装置的运行特性约束﹝式(10)、(11)﹞、变量的上下限约束﹝式(12)﹞。

(2)

式中：φj,t为冷/热负荷节点j功率；cw为水的比热容；mj,t为流过冷/热负荷节点j的水流量；sj表征节点j的负荷性质，取+1、-1表示冷、热负荷；Twj,t和Trj,t分别为冷/热负荷节点j的进水和回水温度；Tip,t和Top,t分别为管道的进水和出水温度；Ta为环境温度；λ为管道单位长度传热系数；L为管道长度；mij,t为节点i、j之间管道流量；min,t和mou,t分别为流入和流出节点的水流量；Tin,t和Tou,t分别为混合前流入节点的各管道水温度和混合后流出节点的水温度；Ppj,t为负荷节点j循环水泵消耗的电功率；mjN和PpjN分别为mj,t和Ppj,t的额定值。

(3)

式中：fij,t为节点i到j之间管道的天然气流量；Kij为管道常数；pi,t和pj,t为节点i和j的压力；，vij,t表征天然气流动方向，pi,t>pj,t时取+1，反之取-1；fin,t、HP,t分别为压缩机流量以及消耗电功率；pin,t和pou,t分别为压缩机的入口和出口压力；Bk和Zk均为常数，Bk与压缩机k的效率、温度、天然气热值有关，Zk与压缩机k压缩因子和天然气热值有关。

(4)

式中：Psj,t、Qsj,t分别为各节点注入有功和无功功率，PLj,t、QLj,t分别为各节点有功和无功负荷；当节点j为配电网注入微网节点，则Psj,t=Ppg,t，Qsj,t=Qpg,t，Qpg,t为配电网节点注入微网的无功功率；当节点j为蓄电池节点，则放电时Psj,t=Ped,t，充电时Psj,t=-Pec,t，Qsj,t=0；Uj,t和Ui,t为节点j和i的电压幅值；Gji和Bji为节点j和i之间互导纳的实部和虚部；θji,t为节点j和i之间的电压相角差。

燃气发电机组效率与其总有功出力之间关系采用三次模型，即：

(5)

式中：Qfu,t为机组消耗的天然气热功率；a、b、c和d为机组的效率系数；PG,t*为机组总有功出力PG,t与额定有功出力的比值；φwa,t和φsm,t分别为热水型和烟气型吸收式制冷机输入的余热功率；awa和asm分别为缸套水和烟气的余热因子。

(6)

式中：φc1,t、φc2,t、φc3,t分别为热水型、烟气型制冷机、电制冷机的制冷功率；C1、C2和C3为对应制冷机的热力系数；ηhr1和ηhr2分别为热水和烟气的回收效率；Pc,t为电制冷机消耗电功率。φh1,t、φh2,t和φh3,t分别为换热机组、电热锅炉和燃气锅炉制热功率；ηhr3、ηH和ηg分别为换热器、电热锅炉和燃气锅炉效率；PH,t为电热锅炉消耗电功率；fH,t为燃气锅炉消耗天然气流量；qng为天然气热值。

能源站内部的供冷/热平衡方程为

(7)

式中：φcΣ,t/φhΣ,t为总的冷/热负荷需求；N1,t、和N2,t分别为烟气型和热水型吸收式制冷机投入供冷的台数；Nw,t为换热机组投入供热的台数。

能源站燃气机组供电母线功率平衡方程为

(8)

式中：Ppcs,t、Pphs,t分别为能源站冷、热源侧循环水泵消耗的电功率；φpcs、φphs分别为能源站冷、热源侧循环水泵的功率因数角；Phps,t和φhps分别为能源站内部压缩机消耗的电功率和功率因数角；QG,t为燃气机组无功出力；n为供电网络节点个数。

(9)

式中：Ln,t为天然气网供给能源站的气负荷；rui、rdi为PG,t的向上、向下爬坡率。

(10)

式中：Ees,t为储能装置的存储能量，对应储冷罐、储热罐和蓄电池的存储量Ec,t、Eh,t和Ee,t；δes为储能装置的能量损失率；ηesc、ηesd分别为储能装置的蓄能、放能效率；φesc,t为储能装置的蓄能功率，分别对应储冷罐、储热罐和蓄电池的φcc,t、φhc,t和Pec,t；φesd,t为储能装置的放能功率，分别对应储冷罐、储热罐和蓄电池的φcd,t、φhd,t和Ped,t；Ees,min、Ees,max分别为Ees,t的最小、最大值；φesc,max、φesd,max分别为φesc,t、φesd,t的上限；uesc,t、uesd,t分别为指示储能装置运行在蓄能、放能状态的二进制变量；Ees,0、Ees,T分别为Ees,t在运行周期起始、结束时段的值。

储能装置状态切换次数限制约束为

(11)

式中Cmax为运行周期中储能装置允许的最大蓄能/放能切换次数。储能装置从闲置状态变为蓄能状态、从放能状态变为蓄能状态、以及脱离蓄能状态，都属于切换了一次状态。

xmin≤x≤xmax，

(12)

式中变量x上下限约束包括供电网各节点电压和各支路功率的上下限、供热/冷网各节点温度和各管道流量的上下限、天然气网各管道压力的上下限、能源站内部各设备变量的上下限。

可以看到，由于含有吸收式制冷机和换热机组投入台数及储能装置运行在蓄能/放能状态的离散决策变量，且含有供电网节点功率平衡方程、供冷/热网管道流体在节点的温度混合模型、供气网管道的流量模型、燃气机组输入输出关系模型以及储能装置状态切换次数限制等非线性约束，因此式(1)—(12)描述的IECM动态优化调度模型是混合整数非线性规划模型。

2 优化模型的线性化

由于混合整数非线性规划问题是规划领域最难求解的问题之一，属于NP难问题，若采用SBB等常用混合整数非线性规划求解器进行求解，不但计算速度很慢，而且经常无法获得问题的最优解。因此，通过分段线性化近似、大M法等效转换和含绝对值约束线性转换等方法，对非线性约束进行线性化，从而将建立的IECM日前动态优化调度模型转化为容易求解的MILP模型。

2.1 供冷/热网络运行特性约束的处理

供冷/热网络的运行特性约束式(2)中，负荷功率模型含有连续变量相乘、供水管道温升/温降模型含有指数运算，这些非线性环节可线性化如下：由于供冷/热网中管道温升/温降很小，回水温度大约在12 ℃/50 ℃。因此可假定供冷/热网的回水温度Tr为定值。以供热网为例，定义供水管道中热媒所含热功率与回水系统中对应热媒所含热功率之差为热媒可利用热功率，则管道ij流入节点i热媒所含可利用热功率[17]

φij=cwmij(Twi-Tr).

(13)

由温降模型Twj=(Twi-Ta)exp[-λL/(cwmij)]+Ta，方程两边同时乘以cwmij并减去cwmijTr，可得到管道ij流出节点j热媒所含的可利用热功率

(14)

因而管道ij可利用热功率损失Δφij为Δφij=φij-φji，则有

(15)

式中：Tws为热源点温度；第1个“”号的近似原理详见文献[17]，第2个“”号相当于假定各个供热点的供水温度都接近热源点温度。

因此，热网运行特性模型可线性化如下：

(16)

式中I为与节点i相关联的节点j的集合。

同理，冷网的线性化处理可参考以上方法。

对于循环水泵耗电特性模型的处理如下：Ppj可由冷/热负荷节点j的热媒所包含可利用热功率φj表示，即[18]

(17)

通过采用分段线性化近似处理：将每个时段的函数分为N段来线性逼近，每段引入1个离散变量ki,t和1个连续变量zi,t，如图1所示，则Ppj,t和φj,t的关系可分段线性化为[19]：

图1 循环水泵耗电功率与可利用热功率关系的分段线性化Fig.1 Piecewise linearization of the relationship between electrical power consumption and available thermal power of circulating pump

(18)

式中：γi、ωi分别为第i个分段的斜率、截距；zi、zi+1分别为第i个分段首端、末端的横坐标；Ppj,t(zi+1)、Ppj,t(zi)分别为第i个分段首端、末端的纵坐标。

2.2 供气网络运行特性约束的处理

(19)

(20)

(21)

(22)

(23)

式(22)、(23)中：ξi,t表示压缩机是否运行在第i个调压档位的0-1变量；ki,t为kt的第i个调压档位值；Nk为调压档位总数。

(24)

式中M为人为引入的大数值常数。

2.3 供电网络运行特性约束的处理

供电网络的运行特性模型式(4)中含有三角函数的非线性运算，线性化处理如下：对于实际供电网运行中的大多数场景，节点电压约为1.0(标幺值，下同)，线路两端节点电压相位差很小。因此可假设Uj≈1，cosθji≈1，sinθji≈θj-θi。则式(4)可转化为如下线性等式[20-21]：

(25)

2.4 能源站运行特性约束的处理

对于燃气发电机组输入天然气热功率与输出有功出力关系式(5)中第1个方程，可进行分段线性化近似处理，即将每个时段的函数分为N段来线性逼近，每段中引入1个离散变量Bi,t和1个连续变量Ji,t，则可转化为：

(26)

式中：αi和βi为第i个分段的斜率和截距；Li、Li+1分别为第i个分段首端、末端的横坐标；Qfu,t(Li)、Qfu,t(Li+1)分别为第i个分段首端、末端的纵坐标。

能源站内部的供冷/热平衡方程式(7)中涉及离散变量和连续变量乘积，如烟气型吸收式制冷机总供冷量N1,tφc1,t，运用大M法将其转化为线性约束，令N1,t上限为N1max，引入N1max+1个0-1变量uic1,t(i=0,1,…，N1max)，令φc1s,t=N1,tφc1,t，则有：

(27)

同理，对于热水型吸收式制冷机总供冷量，令N2,t上限为N2max，引入N2max+1个0-1变量uic2,t(i=0,1,…，N2max)，令φc2s,t=N2,tφc2,t，则有：

(28)

对于换热机组总供热量，令Nw,t上限为Nw,max，引入Nw,max+1个0-1变量uih1,t(i=0,1,…，Nw,max)，令φh1s,t=Nw,tφh1,t，则有：

(29)

能源站内部的供冷/热平衡方程可表示为：

(30)

而能源站内部供电侧的功率平衡方程式(8)可线性化如下：

(31)

2.5 储能装置运行特性约束的处理

对于式(11)中储能装置状态切换次数限制约束，含有绝对值计算的非线性项，可通过引入连续变量yesc,t和yesd,t转换为如下线性约束：

(32)

(33)

式(32)和(33)中前2个式子可以保证yesc,t和yesd,t的取值不小于±(uesc,t-uesc,t-1)中的最大值，也就是不小于绝对值|uesc,t-uesc,t-1|。再结合第3个式子，则保证了切换次数不会超过预设的最大值，可以满足原约束的要求。

2.6 MILP模型

综上，对优化模型中的非线性部分进行线性化处理后，IECM日前动态优化调度模型转化为式(34)所示的MILP模型：

(34)

3 算例分析

以图2所示的某个IECM为例，含有13个节点和12段管道的供冷/热网，9个节点和6段管道的供气网，54个节点和78个支路的供电网；能源站内部有燃气发电机、吸收式制冷机和换热机组等能量转换设备及储冷罐和储热罐。购气和购电价格分别为3.5元/m3和1.022 8元/kWh，Cmax=8。IECM在1 d中96时段的冷/热/电总负荷预测曲线如图3所示，天然气总负荷预测曲线如图4所示。采用的计算机配置为Intel(R) Core(TM) i7-9700 CPU@ 3.60 GHz，32 GB内存，优化计算采用的GAMS软件版本为GAMS win64 24.5.6。

图2 某园区综合能源微网结构Fig.2 Structure of an integrated energy campusmicrogrid

图3 冷、热、电总负荷预测曲线Fig.3 Forecast curves of total loads including cold load, heat load and power load

图4 天然气总负荷预测曲线Fig.4 Forecast curve of total natural gas load

采用GAMS中GUROBI求解器求解IECM日前动态优化调度的MILP模型以得到优化运行方案。其中，供电侧、供冷侧和供热侧的结果如图5—7所示。

由图5可看到：园区优先由燃气发电机组发电来供应电力负荷；蓄电池在夜间用电负荷低谷时段充电，在白天用电负荷高峰时段放电，状态切换次数不超过限制要求。由图6可看到：烟气型和热水型吸收式制冷机利用发电余热制冷，提高能量利用效率，因而优先使用其制冷；储冷罐在夜间冷负荷低谷时段蓄冷，在白天冷负荷高峰时段放冷，状态切换次数符合要求。由图7可看到：换热机组利用发电余热制热，提高能量利用效率，因而优先使用其制热；另外，输出相同热量时燃气锅炉消耗燃气费用小于电热锅炉消耗电能费用，因而燃气锅炉比电热锅炉优先使用。储热罐在夜间用热负荷低谷时段蓄热，在白天热负荷高峰时段放热，状态切换次数也符合要求。

图5 园区供电侧的优化运行结果Fig.5 Optimal operation results of power side in the campus

图6 园区供冷侧的优化运行结果Fig.6 Optimal operation results of cooling side in the campus

图7 园区供热侧的优化运行结果Fig.7 Optimal operation results of heating side in the campus

线性化前后IECM优化运行结果的对比见表1，可以看出，线性化前后微网优化运行结果的各项费用差别都较小，线性化后优化运行结果的目标函数比线性化前稍微小一点。在计算时间方面，线性化前采用SBB求解器求解混合整数非线性规划模型的耗时为518 s，而线性化后采用GUROBI求解器求解MILP模型的耗时仅为4.6 s，较线性化前缩短了98.2%的求解时间，计算时间大大减少。可见，所提出的IECM日前动态优化调度的MILP算法在求解速度上具有很大优势。

表1 线性化前后IECM优化运行结果对比Tab.1 Comparison of optimaloperation results of IECM before and after linearization

线性化前后，供冷网各管道两端可利用热功率的最大和最小偏差(百分比)见表2，供热网各管道两端可利用热功率的最大和最小偏差见表3。可以看到，线性化处理后的供冷/热网的偏差都很小，最大偏差都在2%以内。而通过将线性化模型的最优解代入原来各个非线性等式，得到的式(3)的第1和第2式、式(4)的第1和第2式、式(5)的第1式这些等式两边偏差见表4、表5和图8—10。

表2 供冷网络中各管道的可利用热功率偏差Tab.2 Available heating power deviation of each pipe in the cooling subnet %

表3 供热网络中各管道的可利用热功率偏差Tab.3 Available heating power deviation of each pipe in the heating subnet %

表4 天然气网管道流量偏差Tab.4 Flow deviation percentage of each pipeline innatural gas network %

表5 冷/热网循环水泵耗电功率偏差Tab.5 Percentage of power deviation of circulating pump in the cooling/heating subnet %

图8 供电网中各节点有功功率平衡方程偏差Fig.8 Percentage of active power balance equation deviation

图9 供电网中各节点无功功率平衡方程偏差Fig.9 Percentage of reactive power balance equation deviation

图10 燃气发电机组输入天然气热功率的偏差Fig.10 Percentage of the inputnatural gas thermal power deviation of gas generating unit

以上图表中只给出了各个节点在96个时段偏差中的最小值和最大值。可以看到，各个非线性等式两边的偏差都很接近0，表明所提出的IECM日前动态优化调度的MILP算法具有较高的计算精度。

4 结论

本文提出了一种IECM日前动态优化调度的MILP算法，并通过某个IECM算例分析得到以下结论：

a)所提出方法得到的IECM优化运行方案的运行费用小于混合整数非线性规划模型的结果，经济性更优，且能够大幅度减少求解计算时间，提高了计算效率。

b)所提出方法的计算结果与线性化前结果比较，对应各种变量的偏差都在可接受的范围内，具有较高的计算精度。

当考虑分布式电源及负荷的不确定性时，可通过场景法建立考虑分布式电源及负荷的不确定性的IECM动态优化调度模型。其中，误差场景各约束的线性化处理方法可采用预测场景对应约束的处理方法，而场景转移约束中的绝对值项的线性化方法可采用储能装置状态切换次数限制约束的处理方法。因此，所提出的IECM动态优化调度的MILP算法也可应用于考虑分布式电源及负荷不确定性的情况。