不确定性环境下基于合作博弈的综合能源系统分布式优化

2022-10-31姚文亮王成福赵雨菲张振威

电力系统自动化 2022年20期

姚文亮，王成福，赵雨菲，张振威，管强

（1. 山东大学电气工程学院，山东省济南市 250061；2. 国网山东省电力公司潍坊供电公司，山东省潍坊市 261621）

0 引言

伴随分布式能源种类的多样化、源荷侧不确定性增强等特征，能源互联网成为能源领域的重要研究方向［1-4］。其中，由多个园区级主体构成的多园区综合能源系统在满足多元化用能需求的同时，可有效提高能源利用效率［5-6］。然而，多个园区级主体间如何协调利益冲突、解决信息壁垒，依然是当前影响多主体系统优化运行效果的关键性因素。

目前，在多主体参与的协调优化问题中，博弈论可准确反应参与主体间的相互作用特点，有效促进多主体智能决策。文献［7-8］基于Stackelberg 主从博弈研究了多竞争主体间的互动机制与运行优化；文献［9］结合不同博弈理论，搭建混合博弈模型对多能运营商的能源交易运行策略展开分析与探讨；文献［10］则提出了一种基于重复博弈的综合能源系统优化方法，实现多主体的协同经济运行。

但上述研究主要考虑多主体间的竞争关系，未充分挖掘不同主体之间的潜在合作可能。文献［11］针对多微网间的电能交易问题进行探讨，基于纳什议价方法搭建了多微网议价交易的博弈模型。文献［12］考虑不同储能系统之间的合作关系，提出了一种共享储能动态容量的租赁模型，基于纳什议价进行利益分配。但文献［11-12］忽略了实际系统中随机性因素对合作主体的影响。文献［13］综合考虑负荷、电动汽车的不确定性，搭建纳什议价模型进行求解，但纳什议价的本质为非凸非线性，求解过程复杂且不易寻得均衡解，谈判破裂可能性高。Shapley 值法根据边际贡献率进行收益分配，不仅能够避免利益分配冲突，且模型求解高效。文献［14］引入联合博弈解决多能源枢纽的合作问题，开发了称为合并分裂规则的联盟算法。文献［15］将不同能源载体互联，能源枢纽间相互协作、共享资源，降低了运行成本，且基于贡献率的分配方法解决了能源枢纽间的利益冲突。文献［16］提出了一种基于能量共享的多主体协同优化策略，并采用Shapley 值对联合运营效益进行分配，实现多主体实时协同优化。

然而，文献［14-16］虽然考虑了多主体间的潜在合作可能，并引入Shapley 值分配法保障合作顺利进行，但所述各主体间的能源交易互动过程较为模糊，且合作联盟中各主体须共享内部重要信息，存在隐私泄漏的隐患。近年来，传输信息量少、迭代速度快的交替方向乘子法（alternating direction method of multipliers，ADMM）［17-19］被广泛应用于多主体协调优化。文献［17-19］采用ADMM 在保证参与主体隐私安全的前提下共享资源，实现多主体的安全经济运行。此外，伴随着可再生能源大规模并网，当参与合作联盟的系统数量增多时，可再生能源出力的异质性会影响联盟中合作主体的决策行为。如何在保证参与主体隐私安全、避免利益分配冲突的前提下，降低不确定性带来的风险并提高多园区综合能源系统之间的合作可靠性，亦是亟待解决的问题。

据此，本文计及风光出力的不确定性影响，对多园区系统间的合作行为进行建模，提出不确定性环境下基于合作博弈的综合能源系统分布式优化模型。首先，充分挖掘多园区系统间的潜在合作关系，引入合作博弈并采用Shapley 值法避免多主体利益分配冲突；然后，引入ADMM 解耦多主体间的深度交互，通过迭代互动获取最优交互能量功率，保障合作顺利进行；最后，利用场景分析法描述风光出力的随机波动特征，并通过条件风险价值（conditional value at risk，CVaR）量化处理，为多园区综合能源系统协调优化下的不确定性影响分析提供参考依据。

1 多园区综合能源系统能量交互

1.1 多园区系统能量交互框架

本文研究场景由上级配电网、天然气网络（统称配能网）与多个园区级综合能源系统（park-level integrated energy system，PIES）组成，其结构框架如图1 所示。

图1 多园区综合能源系统能量交互框架Fig.1 Energy interaction framework of multi-park integrated energy system

一方面，考虑到不同PIES 风光出力特性存在差异，将多个园区级系统互联。园区级系统除直接与上级配能网进行交易外，还可根据系统内部能量盈缺情况进行能量交互，减少与上级配能网的能源交易量。另一方面，为避免合作联盟中各主体内部信息泄漏，保护合作主体隐私安全，多个园区级系统间不共享重要信息，仅实时传递下个调度时刻能量交互值，以保障多个PIES 间合作的顺利开展。

1.2 PIES 模型

本文研究的PIES 模型主要由电网、天然气网络及能源耦合设备构成。电网部分主要包含风力发电机（wind turbine，WT）、光伏（photovoltaic，PV）发电等可再生能源发电单元与储电（electricity storage，ES）装置；天然气网络主要包含燃气站（natural gas station，NGS）、储气罐（gas store，GS）；能源耦合设备包括电转气（power-to-gas，P2G）设备和燃气轮机（gas turbine，GT）。此外，园区级系统既可直接与配能网进行能源交易，也可同其他园区级系统进行能量交互。PIES 中相关模型构建如下。

1.2.1 能源耦合设备模型

1）GT

式中:ηg为燃气轮机g的气转电转换效率；HGV为天然气高热值，取39 MJ/m3；Qg，t和Pg，t分别为t时刻燃气轮机g消耗的天然气功率与发电功率。

2）P2G 设备

式中:ηm为P2G 设备m的电转气效率；Pm，t和Qm，t分别为P2G 设备m在t时刻消耗的电能和出气功率。

1.2.2 储能设备模型

本文选取蓄电池、压力罐分别存储电能与天然气，储能装置容量由充/放电时刻的充/放电功率和储能装置的充/放电效率共同决定。储能装置容量表达式为:

1.2.3 系统能量交互

设TPDN，t和GNGN，t分别为t时刻各园区级系统与配能网间交易的电功率和天然气功率，其值为正表示购买能源，其值为负表示出售能源，满足如下约束:

2 基于合作博弈的分布式优化模型

2.1 多园区综合能源系统合作博弈

多园区综合能源系统合作运行时，利益的分配方案是影响系统是否加入联盟的重要因素。合作博弈强调的是集体理性，研究的核心问题是参与主体如何开展合作与分配合作获得的收益。别为园区级系统t时刻购买、出售天然气的价格。

模型剩余约束见式（1）至式（5）及3.3 节。

2.2 基于Shapley 值的利益分配方案

根据合作博弈的前提，合作联盟总成本应低于各系统独立运行时的成本之和，否则合作破裂。假设e(S)为多园区综合能源系统因合作博弈而增加的收益，其表达式为:

式中:v(S)为联盟的总收益；xe为博弈成员e参与合作博弈前的收益；S为不同联盟组合的集合。

Shapley 值法根据成员对联盟的边际贡献率进行利益分配，即博弈成员分得的利益等于其为所参与联盟创造的边际利益的平均值［20-21］。具体利益分配描述如下:

式中:N为参与博弈的园区级系统集合；φe为博弈成员e获得的收益；ϕ(S)为成员所应分得利益的权重；S{e}为从集合S中排除成员e后的集合。

2.3 基于ADMM 的分布式优化模型

多主体协同优化时，各主体内部重要信息无法全部共享，若采用集中式调度不仅难以描述各主体间的能量交互过程，还会带来隐私泄漏的风险。本文采用ADMM 解耦多园区系统间的耦连，通过迭代互动获取最优能量交互值。

2.3.1 ADMM 原理

ADMM 将全局性问题分解为多个较小、较容易求解的局部子问题，通过协调子问题得到全局最优解。现有如下优化问题:

式中:x和z分别为优化问题f(x)和g(z)的优化变量。

将上述等式约束转换为无约束优化问题，表达式为:

式中:Lp(x，z，λ)为增广拉格朗日函数；λ为拉格朗日乘子；ρ为正二次项惩罚项系数。

ADMM 的思想是在求解某变量时将其余变量视作常量，并使用最新迭代结果。迭代过程为:

式中:rk+1和sk+1分别为原始残差和对偶残差；εprimal和εdual分别为允许的原始残差和对偶残差的误差上限。

2.3.2 解耦机制

以2 个园区系统a、b为例观察合作博弈模型，发现系统间通过联络线耦合。联络线上传输能量与互联系统关联，若直接对联络线约束进行松弛，则无法实现解耦。因此，需将联络线约束改写成如下形式:

式中:Ta，t、Ga，t和Tb，t、Gb，t分别为求解系统a和b的子问题时电、气联络线在t时刻传输的电力、天然气功率。联络线约束转换过程见附录A 图A1，等价后的2 条联络线满足同样的传输容量约束。

随后，对联络线约束应用增广拉格朗日进行松弛，得到系统a第k次优化问题目标函数为:

模型其余约束条件同上，具体求解步骤见附录A 图A2。

3 计及风光出力影响的不确定性模型

3.1 预测误差模型与场景选取

通常，风光出力预测误差概率密度分布函数满足正态分布，但实际概率分布未知，若直接采用正态分布会产生较大误差。非参数估计可根据数据自身特点、性质来拟合分布，因此，本文基于风光历史数据，采用核密度估计获得每个调度时刻可再生能源出力预测误差的实际概率密度分布函数。

当获取实际概率分布函数后，采用拉丁超立方抽样进行随机抽样，得到所有园区级系统各调度时刻风、光出力预测误差的样本集合和对应概率。为保证样本准确度，一般抽取的样本容量较大，考虑到计算效率问题，必须在保持样本拟合精度的同时对样本进行削减。本文采用同步回代消除法［22］对场景进行缩减，以园区级系统1 为例，假设其场景集合ws={ws，0，ws，1，…，ws，t，…，ws，T}，其中ws，t为场景s在t时刻的出力值，各场景出现概率之和为1。不同场景间的距离表达式为:

式中:wj，t为集合wj中元素，表示场景j在t时刻的出力值。

场景削减通过对相似场景进行聚类，得到具有相应概率的多场景集合，使样本集较好地逼近原始场景集。同时，场景削减要满足削减后的场景集之间概率距离最小，即满足:

式中:J为场景削减技术删除掉的场景集。

3.2 基于CVaR 的风光不确定性量化

CVaR 是指在某一给定置信度下，投资组合的损失可能大于给定风险价值的平均损失值。本文基于CVaR 量化可再生能源出力不确定性带来的风险，以附录B 图B1 所示的概率密度分布函数为例，展示风、光出力不确定性下的CVaR 示意图。当实际出力超过区间可接纳区间上限时，系统将采取弃风、弃光等措施；当实际出力低于可接纳区间下限时，会导致失负荷。

超过和低于可接纳区间上、下限的部分利用CVaR 进行量化:

结合场景法与CVaR，得到系统a最新的目标函数表达式为:

3.3 模型约束条件

3.3.1 能源耦合设备约束

燃气轮机和P2G 设备的运行约束如下:

3.3.3 电、气网络约束

模型除满足设备运行约束之外，还需满足电、气网络约束。其中，电网部分采用直流潮流模型，相关约束为:

3.4 模型线性化处理与凸约束松弛

3.4.1 线性化处理

模型中非线性部分为表达式（21），采用线性化方法［24］进行处理。具体处理方法如下:

式中:D为分段区间数量；P为模型求解后实际用到的区间数量；Xi为分段区间i的左端点；δi为代表分段部分的连续性变量；ψi为用来保证分段函数连续性的二进制数。如果δi＞0 且2 ≤i≤D-1，则对于1 ≤j＜i，有δj=1。这意味着，如果一个分段区间被使用，那么其左边所有分段区间必须被完全使用。

3.4.2 二阶锥松弛

管道流量方程为非凸非线性方程，采用分段线性化方法处理时，仍无法保证其凸性。因此，本文采用二阶锥规划对管道流量约束进行松弛。由于式（26）中的管道流量约束含有根号项，无法直接松弛，要先对管道流量约束进行简化处理。引入πk，t、πh，t和辅助变量πkh，t分别t时刻表示节点k、h处气压的平方和管道kh气压平方的差值，简化后表达式如下:

此时，式（28）中仍包含0-1 变量，因此采用McCormick 方法和绝对值线性化方法［25］将其转化为不等式约束，并对式（28）中最后一项进行二阶锥松弛。转换后约束如下:

式中:Δπmax为节点k和h压力差的最大值。

为防止松弛后的气压差被增大，使松弛后的约束尽量收紧，在目标函数中引入管道气压差惩罚项，增加所提模型的紧凑性。

式中:F和F'分别为添加惩罚项前、后的目标函数；Ωkh为气压差惩罚项系数。

4 算例分析

4.1 算例基本描述

本文搭建3 个含风、光可再生能源的PIES 进行仿真分析。PIES1 以风力发电为主，PIES2 和PIES3 以光伏发电为主，调度周期为24 h，各系统与配能网的能源交易采用分时价格，见附录B 表B1。园区级系统与配能网间交易的电、气功率最大值分别为300 kW 和200 kW。 PIES1-PIES2、PIES1-PIES3、PIES2-PIES3 之间允许传输电功率的最大值分别为150、100、100 kW；允许传输天然气功率的最大值分别为120、100、100 kW。3 个系统一天中电、气负荷与新能源出力的预测值、算例结构图及设备的相关参数详见附录B。算例求解平台为Python，调用求解器为Cplex，计算机配置为Intel Core i5-10500 CPU @ 3.10 GHz，8 GB RAM。

为验证本文所提模型的有效性和准确性，设置4 个案例进行对比分析:

案例1:3 个PIES 独立运行。

案例2:3 个PIES 合作博弈，基于Shapley 值进行利益分配，采用集中式调度求解。

案例3:在案例2 合作博弈与利益分配方法的基础上，采用ADMM 分布式算法求解。

案例4:在案例3 的基础上，结合场景分析法与CVaR 计及风、光出力的不确定性影响。

4.2 合作博弈结果分析

为验证多个园区级系统合作的必要性，对案例1、2 和3 分别进行求解，3 个案例下系统运行成本的结果如表1 所示。

表1 不同案例下运行成本对比Table 1 Comparison of operation cost in different cases

当多园区系统合作时，系统运行总成本显著降低。与案例1 相比，案例2 和3 运行总成本分别减少了6 085.11 元和6 084.42 元；案例2 中PIES1 至PIES3 的运行成本分别减少了2 989.45、1 894.65、1 201.01 元；案例3 中PIES1 至PIES3 的运行成本分别减少了2 989.15、1 892.21、1 203.06 元。同时，案例2 和3 的运行总成本与各PIES 收益均无较大差异。

案例2 下各PIES 间一天内交互的电、气功率如图2 所示，能够发现PIES1 向PIES2 和PIES3 传输的能量较多，因而所获收益最大。以PIES1-PIES2为例，若传输功率为正，则表示PIES1 向PIES2 传输能量。分析图2 可知:当多园区系统进行合作时，在00:00—05:00、20:00—24:00 这2 个时段，PIES1 向PIES2 和PIES3 输送电能与天然气，这是因为上述时段PIES2 和PIES3 新能源出力较低，不仅无法满足电负荷平衡且无多余的电能供P2G 设备转换为天然气；而PIES1 处于余电状态，不仅对外输送电能，还通过P2G 设备将电能转换为天然气对外输送。在时段10:00—14:00，PIES1 处于缺能状态，PIES2 和PIES3 向其输送能量，可以发现系统间的能量交互中减少了系统与配能网的能量交易，提升了整体经济效益。值得一提的是，仔细观察时段10:00—14:00 的能量输送情况，发现PIES2 和PIES3 虽不处于缺电状态，但仍然进行能量交互，这是因为PIES1 在该时段需要较多的电能和天然气，尽管PIES2 或PIES3 在某时刻满足能量平衡后剩余能量较多，但由于联络线传输功率存在限制，PIES2或PIES3 只能通过对方将能量最大限度地输送给PIES1，此时PIES2 和PIES3 之间是互为传输媒介的作用。其余时段分析同上，本文不再赘述。

图2 案例2 下各系统间能量交互Fig.2 Energy interaction between systems in case 2

此外，以PIES1 为例，对比案例1 和2 下系统内部设备的优化出力，结果如图3 所示。图中:Pdeal为与配电网交易的电功率；Gdeal为与配气网交易的天然气功率。

图3 PIES1 内部设备出力Fig.3 Output of internal equipment of PIES1

案例1 中，在00:00—05:00、20:00—24:00 时段，PIES1 新能源出力较多，燃气轮机保持较低出力，多余的电能一部分通过P2G 设备转换为天然气供应气负荷需求，另一部分由储能装置存储或售卖给配电网；在10:00—14:00 时段，PIES1 处于缺电状态，P2G 设备几乎不运转，除燃气站出力外，还需储气装置放气并外购天然气。当多园区系统合作时，在00:00—05:00、20:00—24:00 时段，PIES1 将多余的电能传输给PIES2 和PIES3，P2G设备出力减少；在时段10:00—14:00，通过PIES2、PIES3 的能量输送，PIES1 在满足内部负荷需求的同时减少了外购能量。

综上分析，各PIES 通过合作博弈，减少了与上级配能网的能量交易，优化了内部设备出力，提高了整体经济效益。同时，Shapley 值法根据合作主体对联盟的贡献率进行利益分配，避免了多主体间利益的分配冲突。

4.3 分布式优化算法分析

相对于集中式调度方法，分布式求解算法能够反应参与合作主体间的能量交互过程，通过迭代实现最优经济运行，最大限度地保护各参与主体的信息隐私安全。

由表1 分析可知，采用分布式优化算法时系统的运行成本和各系统收益与采用集中式调度方法时基本相同。为进一步验证ADMM 的有效性，对案例2、3 的运行结果进行对比分析，求解结果如表2所示。其中，案例3 中二次惩罚项系数均为1。采用分布式算法的求解时间虽有所增加，但仍能满足实时调度的要求。

表2 案例2、3 系统运行结果对比Table 2 Comparison of system operation results in case 2 and case 3

图4 为案例2、3 中PIES1 和PIES2 间交互能量对比。由图可知，采用分布式算法时各调度时段交互的能量与集中式调度基本一致，进一步验证了分布式算法可有效替代集中式调度方法。

图4 案例2、3 能量交互对比Fig.4 Comparison of energy interaction in case 2 and case 3

由于ADMM 为串行方法，在实际应用中，系统的迭代顺序不同可能导致结果不同。对此，本文更换迭代顺序进行求解对比，结果如表3 所示。

表3 不同迭代顺序求解结果Table 3 Solving results in different iteration sequences

分析表3 可知，当迭代顺序不同，迭代次数与运行总成本存在细微的差异，但误差仍在允许范围之内。综上所述，采用ADMM 可较好地代替集中式调度方法，在保护合作主体信息隐私安全的前提下，实现多主体协调优化。同时，分布式模型可快速、准确地给出各调度时刻能量交互功率，降低了因模型求解复杂而使谈判破裂的可能，保证各主体能够安全顺利地开展合作。

4.4 风光出力不确定性影响分析

本文采用场景分析法表征可再生能源出力的不确定性，各园区级系统的场景削减结果见附录C 图C1。将每个园区级系统削减后的典型场景进行组合，获得27 个场景及对应概率。设定置信区间为95%，对所有场景进行求解，结果如图5 所示。

图5 不同场景组合下的求解结果对比Fig.5 Comparison of solving results with different scenario combinations

根据场景出现概率，对系统运行成本和CVaR 进行概率加和，计算结果为21 107.02 元和1 290.19 元。置信区间影响风光可接纳区间的上、下限，反映了决策主体对风险的偏好。为进一步分析不同置信区间对系统调度的影响，本文改变置信区间进行求解，结果如表4 所示。

表4 不同置信区间下的求解结果对比Table 4 Comparison of solving results with different confidence intervals

由表4 可知，置信区间越小，风险价值费用越高，但不考虑CVaR 时的运行成本有所降低。当置信区间增大时，合作联盟虽能规避更多不确定性带来的风险，但系统备用容量增高，可能导致合作联盟整体经济性变差。因此，多园区综合能源系统开展合作博弈时，决策者可通过协调系统备用成本与风险价值成本选择最优区间，在规避一定风险的同时最大化经济收益，减少不确定性带来的实际损失。