文献学习——考虑混合储能系统选择的基于改进蜂群算法的热电联产微网多目标经济优化调度

摘要：在考虑混合储能系统模型选择的基础上，基于改进的人工蜂群算法（ABC），建立了冷热电联产微电网经济优化的多目标调度模型。为了对以往研究中的单目标模型进行升级，将模型的优化目标设定为微电网的日发电调度成本和日环境污染物处理成本最小，并采用基于天牛须搜索算法的改进ABC算法（BAS-ABC）进行求解。为了更好地应用理论成果，本文以上海某区某并网热电联产微电网夏季典型日实际电力负荷数据为研究对象，模拟了采用分时计费制度条件下的多目标经济调度模型。结果表明，优化后的成本比传统成本的最小成本要低。此外，当热电联产微电网的目标函数为发电成本最低或环境成本最低时，优化结果是矛盾的，不可能同时实现最优的发电成本和最优的环境成本。这表明多目标优化调度方法综合平衡了系统的经济效益和环境友好性。

1.引言

微电网作为一种新型的可再生能源供应和管理技术，将各种分布式电源、负载、储能单元和控制装置组合在一起，形成一个单一的可控单元，可以同时向用户提供电能和热能。此外，微电网提供了一个有效的分布式能源综合利用的技术手段[1,2]，其中冷热联产（CCHP）微电网是众所周知的综合能源系统[3]。

因此，热电联产微电网的优化研究主要包括优化各时刻的微功率输出，以及系统中各设备的容量；通过这样的研究，可以提高热电联产微电网的经济效益和环境效益[4,5]。

国内外对热电联产微网优化进行了一定程度的研究。一项工作提出了基于微网内能量传递的热电联产微网通用模型结构，并采用数学规划方法对模型进行求解[6-8]。此外，在微网调度优化计算中加入分时电价和附加收益，建立了热电联产微网动态经济调度模型，对网内各微源出力进行优化[9,10]。有学者在模拟退火算法的基础上提出了一种改进的人工鱼群算法来求解多微网模型[11]。同时，也有学者利用鸟群算法求解微网运行的两个单目标函数（维护成本和环境效益）的调度模型，其改进算法表现出很大的可行性[12-14]。

为了提高热电联产系统的稳定性和经济性，本文在热电联产系统中引入了混合储能系统。根据一些技术和经济细节，构建了混合储能系统。本文以上海地区某热电联产系统为研究对象，以某典型夏日为调度数据。目标函数使当天的发电成本和环境成本最小化。综合分布式电源约束，采用24小时调度周期，采用改进的ABC算法。得出在发电成本最低、环节成本最低和多目标总成本最低的三种条件下，各电源的最优输出，且三种情况下的成本进行对比和分析。

本文的主要贡献可以概括如下。(1)我们对各种储能单元进行了综述，并选择了最优的储能单元组成混合储能系统（HESS）。建立了一种新的CCHP模型，混合储能系统创新地替代了传统的电池储能系统。(2)采用BAS-ABC算法对多目标模型进行优化，实现了案例微电网的经济运行（见图1）。

图1 微电网的组成和构成

2.CCHP微网分布式发电模型

2.1CCHP微网的结构

本文构建的CCHP微电网包含以下部分：微型燃气轮机（MT）、风力涡轮机（WT）、光伏（PV）电池、溴化物制冷机、空调（AC）、燃料电池（FC）和HESS。

图2 冷热电联供微网系统结构

图2为冷热电联供微网系统的典型结构。CCHP微电网的核心装置是微型燃气轮机，可以实现能量的分步利用。利用天然气燃烧过程中产生的高级热能驱动微型燃气轮机发电。剩余的热能通过回收装置送至微电网系统所需的冷负荷或热负荷。冷热电联产系统中空调和微型燃气轮机产生的余热主要通过溴化物制冷机等回收装置来满足供热和制冷负荷需求。系统用电负荷需求主要由WT、PV、MT、FC、HESS和大电网提供。如果上述分布式电源的有功输出在满足微网用电负荷需求后仍有剩余电量，则剩余电量可送入混合储能系统或出售给大电网[15,16]

2.2风力发电机输出功率模型

细节见文章：风力发电机输出功率模型_风力发电数学模型-CSDN博客

风机的功率输出模型主要可以分为线性输出模型、二次输出模型、三次输出模型以及实际测量条件下的风机特性模型，其中三次输出模型是比较常用的：

式中，a、b、c、d是通过拟合厂家提供的输出功率风速曲线得到的。其中“PWT”为机组实际输出功率，“PWTr”为机组额定输出功率。

2.3光伏发电模型

光伏发电系统的输出功率可表示为：

其中，温度系数k取值为0.0047/℃，参考温度Tr为25℃。

2.4微型燃气轮机模型

微型燃气轮机是目前最成熟、最具商业应用竞争力的分布式发电技术。本文选用的微型燃气轮机为凯普斯通C65微型燃气轮机，MGT发电效率和散热损失系数在不同功率下的值可通过厂家提供的有关参数拟合得到，通常拟合为MGT的部分负荷率的一元n次函数：

式中，PLR为部分负荷率，表示燃气轮机实际输出电功率与额定最大电功率之比。

最终，可得其发电效率与输出功率的函数关系为：

式中，hMT为微型燃气轮机发电效率；PMT为微型燃气轮机输出功率；65为凯普斯通C65微型燃气轮机的额定输出功率；PMT/65为输出功率除以额定功率，即为负载率。

同时，以微型燃气轮机为核心的热电联产系统的数学模型为：

式中，QMT表示微型燃气轮机工作时可以利用的热量；Pe为有功功率；he表示发电效率；hl为热量流失系数；Qhe为可以向用户提供的热量；Qco为向用户提供的制冷量；Khe和Kco分别表示产热系数和制冷系数；VMT为天然气量；LHV为天然气低热值，9.7 kWh/m3。

本文优化周期中单位时间间隔∆t为1h，天然气低热值为9.7 kWh/m3（LHV）。使用微型燃气轮机的成本包括燃料、运行和维护成本。MT的燃料成本表示如下：

CMT为微型燃气轮机的燃料成本；Cfuel为天然气价格（元/m3），这里取2.5元/m3；PMT为微型燃气轮机输出功率；hMT为微型燃气轮机发电效率；LHV为天然气低热值，9.7 kWh/m3。

2.5空调模型

空调能增加对输出能量的吸收，并利用电能制热或制冷。空调的数学模型为：

式中，Qair,h为空调产生的实际制热能量；Pac为空调消耗功率；ηLh为热能损失系数；ηah为加热效率；Oair,h为制热能效比，制热量/制热消耗功率；Qair,r为空调产生的实际制冷能量；ηLr为冷却能量损失系数；ηar为冷却效率；Oair,r为冷却能效比，制冷量/制冷消耗功率。

2.6燃料电池模型

燃料电池是一种将燃料中的化学能直接高效地转化为电能的装置，与传统的发电方式相比具有很大的优势。在本文中，燃料电池发电效率与输出功率的关系函数（根据厂商提供的数据，经过拟合可以得到发电效率与输出功率的关系）为：

2.7混合储能系统的选择与建模

储能系统极大地提高了微电网的稳定性、电能质量和供电可靠性，具有极大的市场前景。储能系统在微电网中的功能包括提供短期供电、削峰和填谷；改进分布式供电；提高微电网从并网运行状态切换到孤岛运行状态时的电能质量。基于这些功能和微网的特点，微网对储能系统的性能要求可以概括为：高能量密度；高功率密度；储能效率高；高低温性能；环境友好。目前，市场上有多种储能单元可供选择，每种储能单元都有不同的技术和经济细节。各类型储能单元的技术经济详细情况如表4所示。

表1不同类型储能装置的技术经济细节

通过对表1的详细分析，我们对不同类型储能单元的各项指标进行模糊综合评价。首先利用层次分析法（AHP）确定各指标的权重，然后利用相应的矩阵计算各指标的最大特征根λmax和权重Wi。接下来，我们利用AHP软件建立各类储能技术经济指标分析层次模型，如图18所示。确定一致性的公式如下：

图3 AHP选择的技术和经济

一般认为，当CI<0.1时，评价矩阵的一致性较好，相应的评价指标权重较大；因此，指数更为重要。反之，CI越大，评价矩阵的一致性越差，相应评价指标的权重越小，该指标的重要性越低。经过技术经济细节的权重比较和分析，得到最终的权重数据，如图19(a)和(b)所示。可以看出，各类型储能的一致性均小于0.1。其中，超级电容器和锂电池的权重最大，因此本文选择超级电容器和锂电池组成混合储能系统。

图4 (a)各类储能单元指标评价结果；(b)一致性检验结果

储能系统中的剩余电能也称为SOC（表示电池当前所剩电荷量与其额定电荷量，即电池满充时的电荷量之间的比值。这个比值通常以百分比的形式来表示，能够直观地反映电池的充电状态。），在储能运行过程中，即充放电过程中会发生动态变化。SOC的变化与充放电功率、自放电率、充放电效率有关。充放电递归关系如下所示：

式中，SOCHess(t)表示t时刻剩余电能；SOCHess(t-1)表示t-1时刻剩余电能；dHess表示自放电率（自放电率又称荷电保持能力，是指电池在开路状态下，电池所储存的电量在一定条件下的保持能力。）；PHess(t)表示t时刻功率值，为正表示充电，为负表示放电；∆t表示充电或放电时间长度；EHess,N表示容量；ηHess,c表示充电效率；ηHess,d表示放电效率。

3.热电联产并网微电网多目标经济优化模型

3.1目标函数

微网发电的总经济成本包括燃料成本、投资折旧成本、运行维护成本以及微网与大电网之间的能源交互成本。经济可行性是能源生产企业最重要的因素之一。为了保证微电网的经济性，发电成本应尽可能小。本文中，因为单位时间间隔∆t选取的是1小时，因此离散形式的发电成本目标函数为：

式中，f1表示发电成本目标函数，需要求最小值；T表示系统运行的总小时数；CFC表示燃料成本；CDC表示折旧成本；COM表示运行和维护成本；Cex表示微电网与大电网能源交互成本；CRH表示微网通过提供冷热负荷能量获得的效益。

在微电网的日常运行中，燃料成本主要由在系统中工作的微型燃气轮机（MT）和燃料电池（FC）产生，因为它们都需要燃料来发电。因此，微电网的燃料成本CFC为：

LHVf表示所用燃料的低热值；Cgas表示天然气价格为2.54元/立方米；N1表示MT和FC的总个数。

折旧成本反映的是系统发电设备在运行过程中的价值损失，与设备本身的初始投资和设备日常工作的工作量有关。本文考虑系统中所有设备的损耗，折旧成本CDC为：

式中，N2表示系统中总设备个数；Cpb,i表示第i个分布式能源设备单位容量的安装成本；ki表示容量因数；qi表示第i个分布式能源设备的投资偿还期；r表示年利率；Pi表示第i个分布式能源设备t时段的电功率。

t时刻微网的运行维护成本主要与各分布式能源设备的输出功率成近似的正比关系，表达式为：

Kom,i表示第i个分布式能源设备单位维护成本。

微网与大网在T时段的交互成本Cex为：

式中，Cs表示售电单位功率价格；Ps表示售电功率；Cb表示购电单位功率价格；Pb表示购电功率。

如果Cex(t) > 0，表示该时期微电网向大电网售电；否则，这一时期的微电网向大电网购电。公式中的电气值都是正的。

系统通过提供冷热负荷能量可以获得一定的效益CRH：

式中，Ch表示提供热负荷能量获得到收益；Kph表示单位热能的价格；Qth表示提供的热能；Cr表示提供冷负荷能量获得到收益；Kpr表示单位冷能的价格；Qtr表示提供的冷能。

然而，微网系统在运行过程中会产生污染，而处理这些排放产生的最小成本是本文的另一个目标函数。当环境成本最优时，污染物的排放量最小。因此，在兼顾经济性和环保性的基础上，对热电联产微电网进行运行调度优化。也就是说，在保证供热、供冷、用电负荷平衡的前提下，优化过程中尽可能降低一天的发电成本，同时降低环境成本。对于并网运行模式下的CCHP微电网，光伏电池和风机提供清洁能源，不会污染环境。同时，微网从大电网购买电能时，污染源包括大电网中的燃气轮机、燃料电池和火电机组。调度周期内的环境成本为：

式中，m为污染物类型总个数，n为产生污染物的微电源数量，j为排放的污染物类型。Cj表示污染物处理的单位费用；aij表示第i个分布式能源设备生产电能时，排放的第j种污染物的排放系数；Pi表示第i个分布式能源设备t时段的电功率。

因此，综合经济成本和环境成本目标，总目标函数为：

式中，λ1、λ2分别表示发电成本和环境成本占总成本的比例；在本文中，两者都设置为0.5。

3.2约束条件

在对冷热电联产微网进行日常调度优化时，需要满足的约束包括电、冷、热约束、微电网功率输出约束以及与混合储能系统相关的约束。

微电网的电力平衡约束如下：

式中，PWT表示风力发电机输出的电功率；PPV表示光伏电池输出的电功率；PMT表示微型燃气轮机输出的电功率；PFC表示燃料电池输出的电功率；PHESS表示混合储能系统吸收的电功率；Pex表示与大电网交互的电功率，正则向大电网售电，负则从大电网购电；PAC表示空调消耗的电功率。

微电网的冷热功率约束表示为：

式中，Qh为燃气轮机、燃料电池等其他设备产生的热能；Qair,h为空调产生的实际制热能量；Qlh为电网上此时热能的需求；Qr为微电网中其他设备产生的冷能；Qair,r为空调产生的实际制冷能量；Qlr为电网上此时冷能的需求。

同时，各分布式电源的输出功率约束为：

混合储能系统的充放电约束由最大连续充放电功率和当前剩余电量共同决定。最大充放电功率约束条件为：

式中，PHESS,c,max(t)和PHESS,d,max(t)分别为第t时段储能系统允许的最大充放电功率；其中，PHESS,C,max(t)和PHESS,D,max(t)分别为储能系统在整个时间段内的最大连续充放电功率。因此，HESS充放电约束如下式所示：

微网与大电网之间的能量交互约束如下：

4.基于甲虫触角搜索的改进人工蜂群算法

4.1传统人工蜂群算法

人工蜂群（Artificial Bee Colony，ABC）算法是一种基于蜂群搜索行为的智能优化算法。它具有较强的全局和局部搜索能力，在单目标优化方面表现优异。但与单目标优化问题不同，多目标优化问题需要找到满足一定约束条件的一系列Pareto最优解，求解难度较大。多目标优化通常涉及收敛性和多样性两个性能细节，因此非支配解在解空间中均匀分布。

基本ABC算法存在以下主要问题：(1)收敛“过早”；(2)局部搜索能力较弱，收敛速度相对较慢。针对这些不足，学者Akay和Karaboga提出了许多改进方法。通过多组实验研究了参数设置对ABC算法性能的影响；结果表明，该算法对问题的维数不敏感，适用于求解高维问题。同时，群体大小（CZ）对ABC算法的性能没有明显影响，即使群体大小很小也能得到满意的解。Konrady等研究了基于跟随蜂与领头蜂之间距离的选择方法。当跟随者和领导黄蜂之间的距离变小时，跟随黄蜂选择领导黄蜂寻找花蜜来源的概率更大，反之亦然。Lee等人将群体多样性引入ABC算法，并根据群体多样性阈值采用不同的搜索公式。

ABC算法背后的基本思想是蜜蜂相互合作，通过分工和信息交换来完成采蜜任务。蜂群分为三种蜜蜂：雇佣蜂、围观者蜂和侦察蜂。雇佣蜂和围观者蜂用于采集花蜜，而侦查蜂用于避免花蜜种类过少。对于每个花蜜源，受雇蜜蜂的觅食路线是基于随机选择的邻居k来探索新的食物来源。Xi = （xi1,xi2,xi3，…,xiD）是人口中第i个食物来源的位置，D是问题的维度。初始阶段设置种群的相关参数，包括种群数量、食物源数量、控制参数、最大循环次数、D维解空间。所有参数确定后，对种群进行初始化。在解空间中，每个食物源的位置初始化如下：

式中，为目前被搜索的个人；j为随机选取的维数；vi,j是个体i在j维中的更新位置；其中，randi,d是介于[0,1]之间的随机数；xi,j是个体i在更新维度j之前的位置；k是i随机选择的邻居；xk,j为个体k更新维度j之前的位置。

围观者蜜蜂根据工蜂提供的信息进行进一步挖掘，根据个体的适应度计算个体选择的概率，然后使用轮盘赌算法选择个体进行进一步探索。当花蜜采集不能在一定时间内得到改善时，被雇佣的蜜蜂就会变成侦察兵，随机寻找另一个花蜜来源。

4.2人工蜂群算法（Artificial Bee Colony，ABC）

4.2.1算法介绍

人工蜂群（ABC）算法是一种基于蜂群的元启发式算法，由Karaboga在2005年提出，用于优化数值问题。该模型由三个基本部分组成：就业和失业的觅食蜜蜂以及食物来源。受雇和失业的觅食蜂（前两个组成部分），在其蜂巢附近寻找丰富的食物来源（第三个组成部分）。该模型还定义了自我组织和集体智慧所必需的两种主要行为模式：觅食者招募到丰富的食物来源，从而产生正反馈；觅食者放弃贫乏的食物来源，从而产生负反馈。

在ABC中，一个人工觅食蜂群（代理）搜索丰富的人工食物来源（一个给定问题的良好解决方案）。为了应用ABC，所考虑的优化问题首先被转换为寻找使目标函数最小化的最佳参数向量的问题。然后，人工蜜蜂随机发现一个初始解决方案向量群，再通过采用以下策略对其进行迭代改进：通过邻居搜索机制向更好的解决方案移动，同时放弃差的解决方案。

在ABC中，人工蜜蜂的蜂群包含三组蜜蜂：与特定食物来源相关的受雇蜜蜂，在蜂巢内观看雇佣蜜蜂的舞蹈以选择食物来源的观察蜜蜂，以及随机搜索食物来源的侦察蜜蜂。观察者和侦察者都被称为无业蜜蜂。最初，所有的食物源位置都由侦察蜂发现。此后，食物源的花蜜被就业蜂和围观蜂开发，这种持续的开发最终会使它们变得枯竭。然后，正在开发已耗尽的食物来源的受雇蜂就会变成侦察蜂，再次寻找更多的食物来源。换句话说，食物来源已被耗尽的受雇蜂成为侦察蜂。在ABC中，食物源的位置代表问题的可能解决方案，食物源的花蜜量对应于相关解决方案的质量（适应度）。受雇蜜蜂的数量等于食物源（解决方案）的数量，因为每只受雇蜜蜂都与一个且仅与一个食物源相关。

4.2.2伪代码

1. 初始化参数：
   a. 定义种群大小（SN，即蜜蜂数）
   b. 最大迭代次数（MaxCycle）
   c. 最大试探次数（Limit）
   d. 问题维度（D）
   e. 搜索范围[xmin, xmax]
2. 初始化蜜蜂种群：
   a. 随机生成初始解 Xi(i = 1, 2, ..., SN)，满足 Xi∈[xmin, xmax] (蜜源位置)
   b. 计算每个解的适应度值f(Xi)
3. 主循环（循环直到达到最大迭代次数或满足收敛条件）：
For cycle = 1 to MaxCycle do

   a. 雇佣蜂阶段：
   For i = 1 to SN do
      i. 对每只雇佣蜂（i = 1 到 SN）：
         - 生成一个新解Vij(通过当前解修改一个维度)：
           Vij = Xij + ij * (Xij - Xkj)
           其中：
               Vij 是新解的第j维度
               Xij 是当前解的第j维度
               x[kj] 是随机选择的另一个解的第j维度
               [ij] 是[-1,1]之间的随机数
               k 是随机选择的邻居（k ≠ i）
         - 应用贪婪选择(在新解和原解之间)，如果Vij优于Xi，则用Vij替代Xi
      ii. 计算新解的适应度值 f(Vij)
   End For

   计算每个解的概率值
   pi = f(Xi) / Σf(Xj)（适应度越高，概率越大）

   b. 观察蜂阶段：
   For i = 1 to SN do
      i. 根据概率选择解 Xi：
      ii. 对每个被选中的解，重复雇佣蜂阶段的操作，产生新解并更新种群。
   End For

   c. 侦查蜂阶段：
   For i = 1 to SN do
      i. 如果某个解 Xi 在连续 Limit 次未改进，则将其替换为随机生成的新解：
         Xi = xmin + rand() * (xmax - xmin)
   End For

   d. 记录当前迭代中的最优解。
End For
4. 返回全局最优解。

4.2.3算法流程图

1、初始化蜂群

2、雇佣阶段

通过如下公式产生新个体：

产生个体的同时，根据新种群和老种群进行贪婪选择，留下好的个体，这样就形成了新的个体。

然后，通过如下公式计算个体的适应度（与函数值反相关）：

3、观察阶段（跳舞来共享信息）

4、侦察阶段

如果，在一定次数内，一个雇佣蜂没有被围观蜂的围观而得到提高，那么到达这个次数后就被解雇了，重新初始化，这被称为侦察，初始化方法和最初一样，如下：

5、算法终止条件

算法迭代运行MaxCycles后，跳出循环终止运行。

4.2.4仿真实例

1、例题

求解下面函数的最小值：

其中，-10≤xi≤10。

全局最优解为：xi =1；fmin=0。

2、分析

这是一个连续的实数函数优化问题，一共有10个变量，因此我们的种群个体要有10个维度，同时设置合理的种群数量，然后即可利用进行人工蜂群算法进行优化。

3、Matlab代码实现

%%%%%ARTIFICIAL BEE COLONY ALGORITHM%%%%

%Artificial Bee Colony Algorithm was developed by Dervis Karaboga in 2005 
%by simulating the foraging behaviour of bees.

%Copyright @2008 Erciyes University, Intelligent Systems Research Group, The Dept. of Computer Engineering

%Contact:
%Dervis Karaboga (karaboga@erciyes.edu.tr )
%Bahriye Basturk Akay (bahriye@erciyes.edu.tr)

clear all
close all
clc

% 设置 ABC 算法控制参数
% 三连点(省略号)...表示续行。
% 当一行内语句太长，可以使用三个点...表示续行，另起一行。
ABCOpts = struct( 'ColonySize',  10, ...            % 雇佣蜂数量+观察蜂数量 
    'MaxCycles', 2000,...                           % 终止算法的最大周期数
    'ErrGoal',     1e-20, ...                       % 错误目标，以便终止算法（在当前版本的代码中没有使用）。
    'Dim',       10 , ...                           % 目标函数的参数维度 dimention
    'Limit',   100, ...                             % 控制参数，以便放弃食物来源 
    'lb',  0, ...                                   % 待优化参数的下限
    'ub',  2, ...                                   % 待优化参数的上限
    'ObjFun' , 'rosenbrock', ...                    % 写出你想最小化的目标函数的名称
    'RunTime',3);                                   % 运行的次数

%生成ABCOpts.RunTime=3行，ABCOpts.MaxCycles=2000列的0矩阵；
GlobalMins=zeros(ABCOpts.RunTime,ABCOpts.MaxCycles);% 运行次数以及终止最大周期数

for r=1:ABCOpts.RunTime                             % 外层循环——运行次数
    
    % 初始化种群
    Range = repmat((ABCOpts.ub-ABCOpts.lb),[ABCOpts.ColonySize ABCOpts.Dim]); % 初始化每一维度的范围
    Lower = repmat(ABCOpts.lb, [ABCOpts.ColonySize ABCOpts.Dim]);             % 初始化每一维度的下限
    Colony = rand(ABCOpts.ColonySize,ABCOpts.Dim) .* Range + Lower;           % 在每一维度的范围内，随机初始化种群

    Employed=Colony(1:(ABCOpts.ColonySize/2),:);                              % 前一半个体为雇佣蜂

    % 评估和计算适应度
    %feval(fun,x1,...,xM): fun - 用于计算的函数；x1,...,xM - 所计算函数的自变量值；
    ObjEmp=feval(ABCOpts.ObjFun,Employed);                      % 计算函数值
    FitEmp=calculateFitness(ObjEmp);                            % 计算适应度

    % 设置雇佣蜂Bas的初始值为0，记录在迭代过程中未变优的次数；
    Bas=zeros(1,(ABCOpts.ColonySize/2));            


    GlobalMin=ObjEmp(find(ObjEmp==min(ObjEmp),end));            % 雇佣蜂中的最优函数值
    GlobalParams=Employed(find(ObjEmp==min(ObjEmp),end),:);     % 雇佣蜂中的最优解

    Cycle=1;                                                    % 从 1 开始
    while ((Cycle <= ABCOpts.MaxCycles))                        % 开始进行迭代

        %%%%% 雇佣阶段
        Employed2=Employed;                                     % 新一代雇佣蜂
        for i=1:ABCOpts.ColonySize/2                            % 遍历雇佣蜂个体
            %fix(X)：删除X中每个数的小数部分，将它们截断为整数；
            %rand：生成(0,1)之间（开区间，不包含0和1两个数）均匀分布的伪随机数；
            %+1是为了保证，每次至少有一个蜜蜂改变；
            Param2Change=fix(rand*ABCOpts.Dim)+1;               % 选择更改的维度(Employed2中的某一列，随机生成1~10的整数)
            neighbour=fix(rand*(ABCOpts.ColonySize/2))+1;       % 选择更改的蜜蜂(Employed2中的某一行，但不能与i相同，随机生成 1~雇佣蜂个数 的整数)
            while(neighbour==i)                                 % 保证需要更改的蜜蜂不和当前蜜蜂相同
                neighbour=fix(rand*(ABCOpts.ColonySize/2))+1;
            end
            % 产生新的群体
            Employed2(i,Param2Change)=Employed(i,Param2Change)+(Employed(i,Param2Change)-Employed(neighbour,Param2Change))*(rand-0.5)*2;
            % 限制范围
            if (Employed2(i,Param2Change)<ABCOpts.lb)
                Employed2(i,Param2Change)=ABCOpts.lb;
            end
            if (Employed2(i,Param2Change)>ABCOpts.ub)
                Employed2(i,Param2Change)=ABCOpts.ub;
            end

        end   

        ObjEmp2=feval(ABCOpts.ObjFun,Employed2);                % 重新评估
        FitEmp2=calculateFitness(ObjEmp2);                      % 重新计算适应度
        % 进行贪婪选择
        [Employed, ObjEmp, FitEmp, Bas]=GreedySelection(Employed,Employed2,ObjEmp,ObjEmp2,FitEmp,FitEmp2,Bas,ABCOpts);

        % 根据适应度正则化为概率
        NormFit=FitEmp/sum(FitEmp);

        %%% 观察阶段
        Employed2=Employed;                                         % 将新老种群同步
        i=1;
        t=0;    %记录轮盘赌的次数
        while(t<ABCOpts.ColonySize/2)                               % 进行观察
            if(rand<NormFit(i))                                     % 进行轮盘赌
                t=t+1;                                              % 记忆轮盘赌的次数
                Param2Change=fix(rand*ABCOpts.Dim)+1;               % 获取观察选择的维度
                neighbour=fix(rand*(ABCOpts.ColonySize/2))+1;       % 选择一个被观察的个体
                    while(neighbour==i)                             % 被观察个体不能与当前个体相同
                        neighbour=fix(rand*(ABCOpts.ColonySize/2))+1;
                    end
                Employed2(i,:)=Employed(i,:);                       % 由于上次的被观察，可能导致观察蜂和雇佣蜂不一致
                % 进行观察，对每个被选中的解，重复雇佣蜂阶段的操作，产生新解并更新种群；
                Employed2(i,Param2Change)=Employed(i,Param2Change)+(Employed(i,Param2Change)-Employed(neighbour,Param2Change))*(rand-0.5)*2;

                if (Employed2(i,Param2Change)<ABCOpts.lb)           % 如果超出下限
                    Employed2(i,Param2Change)=ABCOpts.lb;           % 置为下限
                end
                if (Employed2(i,Param2Change)>ABCOpts.ub)           % 如果超出上限
                    Employed2(i,Param2Change)=ABCOpts.ub;           % 置为上限
                end
                ObjEmp2=feval(ABCOpts.ObjFun,Employed2);            % 计算目标函数值
                FitEmp2=calculateFitness(ObjEmp2);                  % 计算适应度
                % 进行贪婪选择
                [Employed, ObjEmp, FitEmp, Bas]=GreedySelection(Employed,Employed2,ObjEmp,ObjEmp2,FitEmp,FitEmp2,Bas,ABCOpts,i);
            end
            i=i+1;                                                  % i+1，下一个个体被观察
            if (i==(ABCOpts.ColonySize/2)+1)                        % 如果 i 超出索引
                i=1;                                                % 重新置为 1
            end   
        end


        %%% 记录最优个体
        CycleBestIndex=find(FitEmp==max(FitEmp));                   % 寻找最优个体，可能是多个
        CycleBestIndex=CycleBestIndex(end);                         % 选择最后一个，变为1个
        CycleBestParams=Employed(CycleBestIndex,:);                 % 记录它的参数
        CycleMin=ObjEmp(CycleBestIndex);                            % 记录目标函数值
         if CycleMin<GlobalMin                                      % 记录全局最优个体，防止最优个体丢失
              GlobalMin=CycleMin;
              GlobalParams=CycleBestParams;
        end
        GlobalMins(r,Cycle)=GlobalMin;                              % r=3，运行3次，Cycle=2000，循环200次，追踪每一次循环的最优个体

        %% 侦察阶段
        ind=find(Bas==max(Bas));                                    % 查找没有变优次数最多的个体
        ind=ind(end);                                               % 可能有很多个，选择最后一个
        if (Bas(ind)>ABCOpts.Limit)                                 % 判断是否大于限制
            Bas(ind)=0;                                             % 未变优次数清零，重新开始计数
            Employed(ind,:) = rand(1,ABCOpts.Dim) .* (ABCOpts.ub-ABCOpts.lb) + ABCOpts.lb;% 进行侦察操作
            %Employed(ind,:)=(ABCOpts.ub-ABCOpts.lb)*(0.5-rand(1,ABCOpts.Dim))*2;% 进行侦察操作
        end
        ObjEmp=feval(ABCOpts.ObjFun,Employed);                      % 重新评估种群
        FitEmp=calculateFitness(ObjEmp);                            % 计算适应值函数
        fprintf('Cycle=%d ObjVal=%g\n',Cycle,GlobalMin);            % 打印每个cycle的最优个体
        Cycle=Cycle+1;
    end % End of ABC
end %end of runs

%semilogy（A）：半对数图（x轴线性刻度，y轴对数刻度）
%mean(A)：返回一个行向量，各元素为数组A每一列的均值
semilogy(mean(GlobalMins))
title('Mean of Best function values');
xlabel('cycles');
ylabel('error');
fprintf('Mean =%g Std=%g\n',mean(GlobalMins(:,end)),std(GlobalMins(:,end)));

适应度计算函数calculateFitness.m如下：

function fFitness=calculateFitness(fObjV)

%初始化，生成与输入fObjV维数相同的0向量；
fFitness=zeros(size(fObjV));

%找到fObjV中值>0的位置，利用>0时的适应度公式进行计算；
ind=find(fObjV>=0);
fFitness(ind)=1./(fObjV(ind)+1);

%找到fObjV中值<0的位置，利用<0时的适应度公式进行计算；
ind=find(fObjV<0);
fFitness(ind)=1+abs(fObjV(ind));

罗森布洛克函数rosenbrok.m 定义如下：

function ObjVal = rosenbrock(Chrom,~)

% 目标函数的维度
% 返回一个行向量，其元素是 A 的相应维度的长度。例如，如果 A 是一个 3×4 矩阵，则 size(A) 返回向量 [3 4]
Dim=size(Chrom,2); % 获取 Chrom 矩阵的第 2 维的大小，即每个候选解的维度数（或变量的数量）

   
% 计算总体参数
[~,Nvar] = size(Chrom);

% function 11, sum of 100* (x(i+1) -xi^2)^2+(1-xi)^2 for i = 1:Dim (Dim=10)
% n = Dim, -10 <= xi <= 10
% global minimum at (xi)=(1) ; fmin=0
Mat1 = Chrom(:,1:Nvar-1);   %创建一个新矩阵 Mat1，包含 Chrom 矩阵的除最后一列之外的所有列。
Mat2 = Chrom(:,2:Nvar);     %创建一个新矩阵 Mat2，包含 Chrom 矩阵的除第一列之外的所有列。

if Dim == 2 %变量个数为2，即对于二维情况（Dim == 2），直接计算；
    ObjVal = 100*(Mat2-Mat1.^2).^2+(1-Mat1).^2;
else        %变量个数不是2，即对于多维情况，使用了一个双重转置操作来确保最终得到的是一个行向量
    ObjVal = transpose(sum(transpose((100*(Mat2-Mat1.^2).^2+(1-Mat1).^2)))); % 注意这里用了两个转置
end   
  
% 自定义函数结束

贪婪算法GreedySelection.m 定义如下：

%贪婪算法：选出两组雇佣蜂中对应行中更好的那一行，组成新的雇佣蜂群体。
function [Colony, Obj, Fit, oBas]=GreedySelection(Colony1,Colony2,ObjEmp,ObjEmp2,FitEmp,FitEmp2,fbas,ABCOpts,i)

oBas=fbas;      %记录未变优的次数；
Obj=ObjEmp;
Fit=FitEmp;
Colony=Colony1;

%nargin:函数输入参数个数
if (nargin==8)                              %输入参数为8个，对每一行都进行贪婪算法；
    %size(A,dim)：查询数组A的维度，如dim=1，查询行数；dim=2，查询列数；
    for ind=1:size(Colony1,1)               %查询数组Colony1的行数；
        if (FitEmp2(ind)>FitEmp(ind))       %按行选出两个对照组中的最优解，组成新的雇佣蜂群体；
            oBas(ind)=0;                    %变优，将未变优的次数清零；
            Fit(ind)=FitEmp2(ind);
            Obj(ind)=ObjEmp2(ind);
            Colony(ind,:)=Colony2(ind,:);
        else
            oBas(ind)=fbas(ind)+1;           %未变优，将未变优的次数加1；
            Fit(ind)=FitEmp(ind);
            Obj(ind)=ObjEmp(ind);
            Colony(ind,:)=Colony1(ind,:);
        end
    end  %for
end  %if
if(nargin==9)   %输入参数为9个，对指定的第i行进行贪婪算法；
    ind=i;
    if (FitEmp2(ind)>FitEmp(ind))   %按行选出两个对照组中的最优解，组成新的雇佣蜂群体；
        oBas(ind)=0;
        Fit(ind)=FitEmp2(ind);
        Obj(ind)=ObjEmp2(ind);
        Colony(ind,:)=Colony2(ind,:);
    else
        oBas(ind)=fbas(ind)+1;
        Fit(ind)=FitEmp(ind);
        Obj(ind)=ObjEmp(ind);
        Colony(ind,:)=Colony1(ind,:);
    end
end

4、仿真结果展示

4.3天牛须搜索算法

天牛须搜索算法是根据天牛觅食原理开发的智能算法：天牛觅食时，根据左右触角接收到的食物浓度确定飞行方向。天牛会朝食物浓度高的方向飞去。BAS算法是一种个体智能算法，它只需要一个个体，运算速度快。

步骤1：由于天牛的头部方向是随机的，因此生成随机维数k的单位向量，表示左天线指向右天线的向量：

k是空间维度。

步骤2：更新天牛的触角位置和适应度值如下：

xrt、xlt分别为t次迭代后的天牛右天线和左天线的位置。d是两个天线之间的距离。fright和fleft分别是天牛右触角和左触角的适应度值。

步骤3：更新天牛的位置如下：

式中，xt为第t−1次迭代时天牛的位置；δ是接下来t次迭代的步长因子。

步骤4：判断更新后的适应度值是否优于上次的适应度值。如果优于上一次适应度值，则更新天牛的位置。当达到最大迭代次数时，输出最优解。否则，返回步骤2。

4.4BAS-ABC算法

在基于种群的进化算法中，只有平衡全局搜索和局部搜索才能提高算法的性能。因此，本文利用全局最优解的信息来指导候选解的搜索，从而提高ABC算法的挖掘能力。同时，利用BAS算法的更新规则随机选择个体，更好地引导蜜蜂探索更大的区域，提高效率。具体过程如下图所示。

式中，xbest是目前最适合的个体；xi,rt和xi,lt分别为个体i在下t次迭代中右天线和左天线的位置；φ是介于[0,1]之间的随机数；vi是个体i的更新位置。

5.算例分析

本文以上海市某热电联产系统为调度研究对象，利用典型的夏季日调度数据进行分析。调度周期为24 h，调度间隔为1 h。夏季典型日负荷数据以及根据建立的风力发电机和光伏电池数学模型得到的输出功率曲线如上图所示。每台分布式电源的相关参数包括安装成本、使用寿命、运行维护成本、有功功率上下限等。这些参数如下表所示。所研究的热电联产系统共有2台空调，由于其基本参数相同，故表中仅列出其中一台。

由锂电池和超级电容器组成的HESS系统参数如下表所示。选用额定功率分别为1500kW和2000kW的锂电池和超级电容器组成混合储能系统。

微网和大电网的交易价格采用分时电价机制，如下表所示。交易功率范围为- 150 ~ + 150kW，其中“+”表示微电网向大电网出售电力，“-”表示从电网购买电力。制热或制冷收益为0.1元/(kW h)。

各分布式电源的污染物排放及处理成本见下表。表中数据单位为g/(kW h)，选取典型夏季负荷日进行多目标优化，并与单目标优化进行经济性比较。

图7-9为并网条件下调度周期内微网和大电网各时刻的最优输出。其中，图7为以多目标优化为目标的调度结果，图8为以发电成本优化为目标的调度结果，图9为以环境成本优化为目标的调度结果。

00:00 - 06:00，电价低，购电成本低，但相应的环境污染大。因此，具有环境优势的FC应加大贡献，降低环境成本。同时，由于此时负荷需求较小，微电网可以减少可控机组的输出，电网和FC可以满足负载要求并向HESS发送能量。

在06:00 - 08:00、12:00 - 17:00、22:00 - 24:00时段，电价正常，购电成本增加。因此，发电成本和环境成本是综合考虑的。如果光伏电池和风力发电机不能满足负荷要求，则会全面调用燃料电池和大电网。此外，由于燃料成本低，燃料电池的污染排放很少。在8:00 - 12:00和17:00 - 22:00时段，即电价最高、环境污染最大的时段，微网内可控机组通过向电网出售电能的方式增加发电量，降低发电成本。因此，HESS和FC具有环境成本和发电成本较低的优势，是负荷间隙的首选。8:00 ~ 12:00冷负荷增加，MT有功出力增加，使电力负荷亏缺量达到整个调度周期的最小值。HESS的输出能够满足负荷需求并将多余的电能出售给大电网，提高微电网的经济性；同时，FC处于空闲状态。17:00 - 22:00，HESS和FC优先平衡用电负荷。当负荷平衡不能满足时，HESS放电使剩余容量达到下限，FC满负荷运行，剩余不足的负荷由电网平衡。

采用ABC和BAS-ABC算法求解最优值。人口规模设定为300人；工蜂150只，围观者150只，种群迭代次数200次，维度为96。针对不同目标优化的典型夏季天数的成本如图10所示。图10为不同目标优化后典型夏季天数的成本。可以看出，当CCHP微电网的目标函数为发电成本最小或环境成本最小时，优化是冲突的，即发电成本和环境成本不可能同时最优。当优化目标为发电成本最小时，相应的环境成本比最低环境成本高50.58%。同时，在环境成本最低时，相应的发电成本比最低发电成本高18.89%。然而，当采用多目标模糊优化时，相应的环境成本比最低环境成本高33.70%，发电成本比最低发电成本高6.13%。因此，典型夏季天气的多目标优化结果验证了多目标优化比单目标优化更能兼顾多目标的最优性。而且，这种优化可以综合协调目标的经济绩效和环境绩效。

两种算法最优值的收敛曲线对比如图11所示。可以看出，BAS-ABC算法的收敛速度比传统ABC算法要快，这意味着传统ABC算法存在劣势。因此，本文提出的改进ABC算法优于传统的ABC算法。

6.结论

本文构建了一个混合储能系统，并利用BAS-ABC算法对典型夏季的热电联产系统的发电成本和环境成本进行了联合优化。通过模型仿真，我们可以得出以下结论：(1)与传统的ABC算法相比，BAS-ABC算法的收敛速度更快。(2)与传统的ABC算法相比，BAS-ABC算法求解的最小代价更低。(3)当热电联产微网的目标函数为发电成本最小或环境成本最低时，优化存在冲突，即发电成本和环境成本无法同时优化；反之，多目标优化调度则综合协调了经济与环境保护。本文提出的调度策略具有一定的有效性和参考价值，但也存在一定的局限性和不足。例如，未考虑微电网在孤立状态下的运行情况，或不同定价标准的情况。在未来的模型中，需要加入孤岛状态和多种电价标准作为可变因素，综合研究多种不同的运行工况。

参考文献

[1]Shan J N, Lu R X. Multi-objective economic optimization scheduling of CCHP micro-grid based on improved bee colony algorithm considering the selection of hybrid energy storage system[J]. Energy Reports, 2021, 7: 326-341.