CN115936220A - 一种园区综合能源系统规划及运行优化方法和终端 - Google Patents

一种园区综合能源系统规划及运行优化方法和终端 Download PDF

Info

Publication number
CN115936220A
CN115936220A CN202211589181.7A CN202211589181A CN115936220A CN 115936220 A CN115936220 A CN 115936220A CN 202211589181 A CN202211589181 A CN 202211589181A CN 115936220 A CN115936220 A CN 115936220A
Authority
CN
China
Prior art keywords
carbon
energy
wind speed
scene
illumination intensity
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211589181.7A
Other languages
English (en)
Inventor
杜翼
李源非
陈晚晴
陈劲宇
陈文欣
项康利
陈晗
陈思敏
陈津莼
林晓凡
李益楠
陈柯任
施鹏佳
林昶咏
郑楠
蔡期塬
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
State Grid Fujian Electric Power Co Ltd
Economic and Technological Research Institute of State Grid Fujian Electric Power Co Ltd
Original Assignee
State Grid Fujian Electric Power Co Ltd
Economic and Technological Research Institute of State Grid Fujian Electric Power Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by State Grid Fujian Electric Power Co Ltd, Economic and Technological Research Institute of State Grid Fujian Electric Power Co Ltd filed Critical State Grid Fujian Electric Power Co Ltd
Priority to CN202211589181.7A priority Critical patent/CN115936220A/zh
Publication of CN115936220A publication Critical patent/CN115936220A/zh
Pending legal-status Critical Current

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种园区综合能源系统规划及运行优化方法及终端，包括建立多能流与碳交易量耦合模型；建立基于碳能耦合模型的园区综合能源系统模型；建立考虑碳交易价格以及风光不确定性的场景生成模型；建立基于条件风险价值的园区综合能源系统规划运行双层优化模型。本发明为刻画碳交易量依附于多能流而存在的实际情况，建立多能流与碳交易量耦合模型来描述碳交易量与多能流间的耦合关系，进而将碳交易价格不确定性和计及时空相关性的风光发电不确定性引入规划总成本波动风险建模中，提出了基于条件风险价值风险度量的PIES规划运行双层优化模型，有效避免设备容量配置不足导致无法应对不确定因素造成的规划总成本波动风险的问题。

Description

一种园区综合能源系统规划及运行优化方法和终端
技术领域
本发明涉及综合能源系统规划及运行优化技术领域，尤其是涉及一种园区综合能源系统规划及运行优化方法和终端。
背景技术
园区综合能源系统(Park-level integrated energy system,PIES)具有能耗大、能量密度高的特征，且具备集成各类可再生能源并加以利用的潜力，是实现我国“双碳”目标的重要载体之一。PIES的低碳发展可以通过参与碳交易和发展可再生清洁能源来实现。PIES参与碳交易是指把碳排放权作为自由交易的商品在碳交易市场上进行买卖，利用市场化手段对PIES的碳排放量进行控制，并从中获利。另外，PIES可配置大规模可再生清洁能源助力“双碳”目标实现。同时，大规模发展可再生清洁能源也使得PIES在参与碳交易时具有较好获利的空间。
然而，现有碳交易市场尚处于起步阶段，碳交易价格存在较强的不确定性，加之PIES内部可再生能源发电本身的强功率波动特性，直接影响了PIES稳定参与碳交易的潜力。PIES规划总成本通常包括设备投资成本、燃料成本、维护成本和碳交易成本，由于碳交易价格和可再生能源发电两组不确定性因素影响，极易导致规划总成本与预期值的偏差，即规划总成本波动风险。
因此，建设坚强可靠、低碳、经济的PIES需要规避规划总成本波动风险的出现，为此应从规划阶段就要加以重点考虑。
由于PIES配置了高比例可再生能源，使得PIES在参与碳交易时，可以在满足自身碳排放需求的同时产生盈余碳排放权在碳交易市场出售进行获利。然而，PIES参与碳交易获利时需要以强健的优化配置为基础，避免设备容量配置不足，无法应对不确定因素造成的规划总成本波动风险。
发明内容
本发明所要解决的技术问题是：提供一种园区综合能源系统规划及运行优化方法和终端，有效避免设备容量配置不足导致无法应对不确定因素造成的规划总成本波动风险的问题。
为了解决上述技术问题，本发明采用的技术方案为：
一种园区综合能源系统规划及运行优化方法，包括步骤：
S1、建立多能流与碳交易量耦合模型；
S2、根据所述多能流与碳交易量耦合模型，建立基于碳能耦合模型的园区综合能源系统模型；
S3、根据所述园区综合能源系统模型，建立考虑碳交易价格以及风光不确定性的场景生成模型；
S4、根据所述场景生成模型，建立基于条件风险价值的园区综合能源系统规划运行双层优化模型。
为了解决上述技术问题，本发明采用的另一个技术方案为：
一种园区综合能源系统规划及运行优化终端，包括存储器、处理器和存储在存储器上并可在处理器上运行的计算机程序，所述处理器在执行所述计算机程序时实现以下步骤：
S1、建立多能流与碳交易量耦合模型；
S2、根据所述多能流与碳交易量耦合模型，建立基于碳能耦合模型的园区综合能源系统模型；
S3、根据所述园区综合能源系统模型，建立考虑碳交易价格以及风光不确定性的场景生成模型；
S4、根据所述场景生成模型，建立基于条件风险价值的园区综合能源系统规划运行双层优化模型。
本发明的有益效果在于：为了刻画碳交易量依附于多能流而存在的实际情况，建立了多能流与碳交易量耦合模型来描述碳交易量与多能流间的耦合关系，进而将碳交易价格不确定性和计及时空相关性的风光发电不确定性引入规划总成本波动风险建模中，提出了基于条件风险价值(Conditional value-at-risk，CVaR)风险度量的园区综合能源系统规划运行双层优化模型以应对全寿命周期的总成本波动风险，即通过重点分析园区综合能源系统面临的碳交易价格与可再生能源出力不确定性造成的规划总成本波动风险对系统配置及运行优化的影响，有效避免设备容量配置不足导致无法应对不确定因素造成的规划总成本波动风险的问题。
附图说明
图1为本发明的一种园区综合能源系统规划及运行优化方法的流程图；
图2为一种PIES的典型结构图；
图3一典型PIES的碳交易价格预测误差的频率直方图与核密度估计图；
图4为一典型PIES的夏季、冬季及春秋季典型日的电/热/冷负荷曲线图；
图5为一典型PIES的某全年风速曲线图；
图6为一典型PIES的某全年光照强度曲线图；
图7为一典型PIES在不同季节典型日光照强度和风速均值图；
图8为碳交易价格聚类场景数量与SSE关系图；
图9为光照强度聚类场景数量与SSE关系图；
图10为风速聚类场景数量与SSE关系图；
图11为聚类后的碳交易价格场景；
图12为聚类后的光照强度场景；
图13为聚类后的风速场景；
图14为冬季典型日购气量与购电量；
图15为不同置信水平下的总成本和CVaR关系曲线图；
图16为本发明的一种园区综合能源系统规划及运行优化终端的结构图。
标号说明：
1、一种园区综合能源系统规划及运行优化终端；2、存储器；3、处理器。
具体实施方式
为详细说明本发明的技术内容、所实现目的及效果，以下结合实施方式并配合附图予以说明。
本发明提供的一种园区综合能源系统规划及运行优化方法和终端，适用于综合能源系统规划及运行优化场景。以下结合实施例具体说明。
请参照图1，本发明的实施例一为：
一种园区综合能源系统规划及运行优化方法，如图1所示，包括步骤：
S1、建立多能流与碳交易量耦合模型。
S2、根据多能流与碳交易量耦合模型，建立基于碳能耦合模型的园区综合能源系统模型。
S3、根据所述园区综合能源系统模型，建立考虑碳交易价格以及风光不确定性的场景生成模型。
S4、根据所述场景生成模型，建立基于条件风险价值的园区综合能源系统规划运行双层优化模型。
即在本实施例中，通过重点分析园区综合能源系统面临的碳交易价格与可再生能源出力不确定性造成的规划总成本波动风险对系统配置及运行优化的影响，有效避免设备容量配置不足导致无法应对不确定因素造成的规划总成本波动风险的问题。
请参照图2及图3，本发明的实施例二为：
一种园区综合能源系统规划及运行优化方法，在上述实施例一的基础上，在本实施例中，步骤S1中建立多能流与碳交易量耦合模型的过程具体如下：
其中，由于能量集成器模型式仅体现了能量的传输、转化和存储过程，但不能描述系统多能流与碳交易量耦合关系。碳交易量是指系统的碳排放量与无偿的碳配额之间的差值。本实施例将碳交易量引入下述公式(1)中，建立多能流与碳交易量耦合模型，简称为碳能耦合模型：
Figure BDA0003990825630000051
可简化为公式(2)：
Figure BDA0003990825630000052
上式中，αw(w＝1,2,…,m)为能源集线器中第w种输入能源的碳排放系数，βw(w＝1,2,…,m)为能源集线器中第w种输入能源的无偿碳配额系数，
Figure BDA0003990825630000053
(w＝1,2,…,u)为能源集线器中第τ种能量存储的碳排放系数，ρτ(w＝1,2,…,u)为能源集线器中第τ种能量存储的无偿碳配额系数，L为输出功率向量，C为转换设备耦合矩阵，P为输入功率向量，K为存储设备耦合矩阵，H为储能充放功率向量，E为能源集线器的碳交易量，A为1×m维度能源供给侧碳交易量系数矩阵，A＝[(α1-β1),(α2-β2),…,(αm-βm)]，B为1×u维度的储能碳交易量系数矩阵，
Figure BDA0003990825630000054
M为(n+1)×1维度的碳能输出矩阵，D为(n+1)×m维度的能源供给侧碳能耦合系数组合矩阵，U为(n+1)×u维度储能的碳能耦合系数组合矩阵。
作为上述实施例一的补充，在本实施例中，步骤S2中基于碳能耦合模型的园区综合能源系统模型构建过程如下：
S21、采用由热电联供机组(Combined heat and power，CHP)、燃气锅炉(Gasboiler，GB)、电制冷机(Electric chiller，EC)、电锅炉(Electric boiler，EB)、电储能(Battery energy storage system，BESS)、光伏(Photovoltaic，PV)以及风机(Windturbine，WT)等设备组成的典型园区综合能源系统PIES，即如图2所示。PIES从系统外部购买天然气和电能，利用能源转换设备配合电储能来满足用户的电/热/冷负荷需求，同时，PIES可参与碳交易出售盈余碳排放权进行获利。
在本实施例中，PIES的碳排放来源于一次能源供给端和能量存储设备产生的碳排放，即外电网购电、天然气消耗、PV和WT发电以及电储能产生的碳排放，虽然PV、WT和电储能在运行过程中零碳排放，但是PV、WT和电储能在生产制造和运输过程中会产生大量碳排放，因此，借鉴生命周期分析法可得归一化计量后PV、WT和电储能运行时的碳排放。
S22、对图1中的PIES进行基于多能流与碳交易量耦合关系的建模，如下公式(3)～(5)：
Figure BDA0003990825630000061
Figure BDA0003990825630000062
Figure BDA0003990825630000063
上式中，Pgrid、Pgas、PPV和PWT分别表示外网购电功率、天然气注入功率、光伏注入功率和风机注入功率，PEL、PCL、PHL和Etr分别表示电/冷/热负荷和碳交易量，c1表示供给电负荷的电功率分配系数，c2表示供给电锅炉的电功率分配系数，c3表示供给热电联供机组的天然气功率分配系数，
Figure BDA0003990825630000064
和
Figure BDA0003990825630000065
分别表示热电联供机组的电效率、热电联供机组的热效率、燃气锅炉的热效率、电锅炉的热效率和电制冷机的制冷效率，PES表示电储能的充/放电功率，当大于0时，表示电储能充电，当小于0时，表示电储能放电，αgrid、αPV、αWT、αgas和
Figure BDA0003990825630000066
分别表示外网购电、光伏、风机、天然气和电储能的单位功率碳排放系数，βgrid、βPV、βWT、βCHP、βGB和βES分别表示外网购电、光伏、风机、热电联供机组、燃气锅炉和电储能的单位功率无偿碳配额系数；上式中PIES的碳交易量为外网购电、光伏、风机、天然气和电储能五部分的碳交易量之和，若其中某一部分的碳排放系数大于碳配额系数，则该部分的碳交易量为正，若某一部分的碳排放系数小于碳配额系数，则该部分的碳交易量为负，当为正值时，PIES需要从碳交易市场购买碳排放权，当为负值时，PIES可将多余碳排放权在碳交易市场出售。
S23、c1、c2和c3是为满足能源集线器建模形式要求引入的中间变量，可用各设备单元的输出功率和设备的效率来表示，即如下公式(6)～(8)：
Figure BDA0003990825630000071
Figure BDA0003990825630000072
Figure BDA0003990825630000073
上式中，GCHP表示热电联供机组消耗的天然气功率，
Figure BDA0003990825630000075
表示电制冷机的输出功率，
Figure BDA0003990825630000076
表示电锅炉的热输出功率，
Figure BDA0003990825630000077
表示储能放电功率，
Figure BDA0003990825630000078
表示储能充电功率。
作为上述实施例一的补充，在本实施例中，步骤S3中考虑碳交易价格以及风光不确定性的场景生成模型构建过程具体为：
S31、对碳交易价格不确定性进行建模，具体为：
根据中国碳排放交易官网提供的碳交易交个历史数据(以2017-2021为例)，利用神经网络技术对碳交易价格进行训练并预测得到预测误差，并进一步采用核密度估计方法得到预测误差的频率直方图，如图3所示。
根据图3所示的频率直方图，可以确定预测误差服从零均值的正太分布，因此将预测误差与碳交易价格的预测值进行叠加，可得到实际的碳交易交个概率分布，即如下公式(9)：
Figure BDA0003990825630000074
xe表示碳交易价格随机变量，μe和σe分别表示碳交易价格均值和标准差。
由于光伏和风机发电分别与光照强度和风速有关并呈线性关系，因此针对光照强度和风速的时空相关性和不确定性进行建模，在同一区域相邻时间段内，受气候条件等因素影响，光照强度和风速具有一定的时空相关性，因此从风光空间相关性建模和风光时序相关性建模两个方面进行描述。
S32、对考虑时空相关性的风光不确定性进行建模，包括风光空间相关性建模和风光时序相关性建模。
S32.A、采用Copula函数建立风光联合概率分布来描述风光空间相关性，则风光空间相关性建模的具体步骤如下：
基于光照强度和风速的历史数据，分别建立光照强度和风速的概率分布模型，光照强度的概率密度函数如下公式(10)～(12)所示：
Figure BDA0003990825630000081
Figure BDA0003990825630000082
Figure BDA0003990825630000083
上式中，Γ(·)表示伽马函数，
Figure BDA0003990825630000087
表示在
Figure BDA0003990825630000088
季节t时刻的光照强度，
Figure BDA0003990825630000089
和
Figure BDA00039908256300000810
分别表示在
Figure BDA00039908256300000811
季节t时刻光照强度的形状参数，
Figure BDA00039908256300000812
和
Figure BDA00039908256300000813
分别表示在
Figure BDA00039908256300000814
季节t时刻光照强度的均值及标准差。
对于风速的不确定性，采用威布尔分布来描述，风速的概率密度函数如下公式(13)～(15)所示：
Figure BDA0003990825630000084
Figure BDA0003990825630000085
Figure BDA0003990825630000086
上式中，aWT(ω,t)和bWT(ω,t)分别表示在ω季节t时刻风速的尺度参数和形状参数，σWT(ω,t)表示在ω季节t时刻的风速标准差，v(ω,t)表示在
Figure BDA0003990825630000094
季节t时刻的风速值，
Figure BDA0003990825630000095
表示在
Figure BDA0003990825630000096
季节t时刻的风速均值。
同时，由于风光之间常存在负相关互补关系，Frank-Copula函数可以很好地描述变量之间的非负和负相关关系，因此选用Frank-Copula函数来对光照强度和风速的联合概率分布进行建模，风光联合概率分布表达式如下公式(16)所示：
Figure BDA0003990825630000091
上式中，ρ表示相关系数，uPV＝F(GT(ω,t))，F(·)表示光照强度的累积概率分布，uWT＝G(v(ω,t))，G(·)表示风速的累积概率分布。
另外，由于某当前时刻的光照强度、风速波动与其迁移时刻或前几个时刻的光照强度、风速波动相关联，在时间上具有时序相关性，由于整合自回归滑动平均法可以很好地捕捉到变量的时序波动特征，故而采用整合自回归滑动平均法来建模光照强度和风速的时序相关性。
S32.B、采用整合自回归滑动平均法来对光照强度、风速的时序相关性进行建模，则风光时序相关性建模如下公式(17)～(20)所示：
Figure BDA0003990825630000092
Figure BDA0003990825630000093
▽·GT,t＝GT,t-1   (19)；
▽·vt＝vt-1   (20)。
上式中，▽表示光照强度和风速的后移位算子，dPV和dWT分别表示光照强度和风速差分运算的阶数，εPV,t和εWT,t分别表示光照强度和风速的高斯白噪声，φPV,υ和φWT,υ分别表示光照强度和风速的第υ个自回归项系数，θPV,υ和θWT,υ分别表示光照强度和风速的第υ个滑动平均项系数，pPV和qPV分别表示光照强度自回归项的阶数和滑动平均项的阶数，pWT和qWT分别表示风速自回归项的阶数和滑动平均项的阶数，根据光照强度和风速的历史数据，采用极大似然估计法可估计出自回归项系数和滑动平均项系数以及对应的阶数。
同时，由于采用抽样法可得到离散场景来描述碳交易价格、光照强度和风速在未来可能发生的场景，而重要抽样法在抽样时可以将更多的样本集中在概率分布的尾部，这也更符合CVaR，即符合条件风险价值用于度量经济损失的尾部风险的概念，因此采用重要抽样法获取碳交易交个、光照强度和风速场景。
具体步骤如下：
S33、根据建立的风光联合概率分布式，通过重要性抽样法得到在某一时间断面下光照强度和风速的联合场景，具体为：
S33.a、在[0,1]区间内生成随机数a1和a2；
S33.b、令光照强度的概率分布函数值uPV＝a1，根据构建的风光联合概率分布来求得风速的概率分布函数值uWT＝a2，即求下述公式(21)的解：
Figure BDA0003990825630000101
S33.c、重复步骤S33.a～S33.b，得到一定组数的光照强度和风速的概率分布函数值；
S33.d、由于
Figure BDA0003990825630000102
uWT＝G(v(ω,t))，利用反函数运算，即分别求解GT(ω,t)＝F-1(uPV)和v(ω,t)＝G-1(uWT)，将步骤S33.c所得的随机数样本转换为光照强度和风速的联合场景X＝{xPV,xWT}，所生成场景考虑了光照强度和风速间的空间相关性，然后通过重要性抽样法得到具有时序相关性的光照强度场景yPV和风速场景yWT。
S34、结合生成的风光时序场景yPV和yWT，对生成的具有空间相关性场景xPV和xWT进行修正，具体为：
S34.a、将风光空间相关性场景X中的风速场景xWT从小到大排序，排序后的风速场景记为xWT，为保持相关关系的一致性，将与xWT对应的光照强度场景xPV跟随xWT变动，变动后光照强度场景记为xPV；
S34.b、将时序场景yPV和yWT进行转置，按列向量重构得到方阵y'PV和y'WT，取与xWT相同时间断面且具有时序相关性的风速场景y'WT，将其按照与xWT相同的大小顺序排列，排列后的场景记为yWT，其余各时间断面场景跟随yWT变动，记为yWT,t，根据下述公式(22)修正空间相关性场景中的各时间断面场景：
Figure BDA0003990825630000111
上式中，Y'WT,t为修正后具有时空相关性的风速场景，同理可得修正后具有时空相关性的光照场景X'WT,t。
S35、为了加快计算速度，通过K-means聚类法对生成的相似场景进行聚类，以减少场景数量，然后利用误差平方和指标来确定碳交易价格、光照强度和风速的聚类场景数量，其中误差平方和表达式如下公式(23)所示：
Figure BDA0003990825630000112
上式中，
Figure BDA0003990825630000115
表示第γ个随机变量的第a个类，
Figure BDA0003990825630000116
表示
Figure BDA0003990825630000117
中的样本数据，
Figure BDA0003990825630000118
表示第γ个随机变量第a个类的质心，
Figure BDA0003990825630000119
表示第γ个随机变量的聚类总数，根据公式(23)，可以分别确定碳交易价格、光照强度和风速的聚类场景数量。
由此可得到碳交易价格、光照强度和风速的典型场景。
S36、分别用向量Se、SPV和SWT表示碳交易价格、光照强度和风速的典型场景，由这三者的场景向量构建的典型场景集为S＝{Se,SPV,SWT}N，其中N表示总场景数量。
作为上述实施例一的进一步补充，在本实施例中，步骤S4中基于CVaR的园区综合能源系统规划运行双层优化模型构建过程具体为：
S41、建立基于条件风险价值的园区综合能源系统规划运行双层优化模型的目标函数：
根据条件风险价值的风险度量方法，基于条件风险价值的园区综合能源系统规划运行双层优化模型以全寿命周期规划及运行总成本的条件风险价值最小为目标，建立如下公式(24)～(30)的目标函数：
Figure BDA0003990825630000113
Figure BDA0003990825630000114
Figure BDA0003990825630000125
Figure BDA0003990825630000121
Figure BDA0003990825630000122
Figure BDA0003990825630000123
Figure BDA0003990825630000124
上式中，
Figure BDA0003990825630000126
表示在置信水平θ下规划及运行总成本的风险值，该风险值指不超过给定置信水平的最大规划总成本的阈值，
Figure BDA0003990825630000127
表示在s场景下的规划及运行总成本，
Figure BDA0003990825630000128
表示在s场景下的等年值投资成本，
Figure BDA0003990825630000129
表示在s场景下的燃料成本，
Figure BDA00039908256300001210
表示在s场景下的维护成本，
Figure BDA00039908256300001213
表示在s场景下的碳交易成本，当
Figure BDA00039908256300001211
为负值时，表示卖出，园区综合能源系统获得收益，当
Figure BDA00039908256300001212
为正值时，表示买进，园区综合能源系统付出成本，
Figure BDA00039908256300001214
表示在s场景下的碳交易价格，r表示折现率，J表示设备候选集，lj表示第j个设备的寿命，
Figure BDA00039908256300001215
表示在s场景下第j个设备的容量，cinv,j表示第j个设备单位容量价格，cmain,j表示第j个设备单位功率的维护费用，
Figure BDA00039908256300001216
表示在s场景下第j个设备在t时刻的输出功率，T表示年总运行时间，为8760小时，
Figure BDA00039908256300001217
和
Figure BDA00039908256300001218
分别表示在s场景下t时刻系统的天然气购买量和从外网购电量，λgas和λelec分别表示天然气价格与电力价格，Δt表示系统运行时间间隔，取值为1小时，CVaRθ表示在置信水平下的条件风险价值，θ为置信水平，反映系统对规划总成本波动风险的厌恶水平，本文中均以置信水平来表示风险厌恶水平。
S42、建立基于条件风险价值的园区综合能源系统规划运行双层优化模型的约束条件：
约束条件包括等式约束与不等式约束，等式约束为利用碳能耦合模型建立的系统功率平衡约束，不等式约束为设备运行时出力上下限约束以及设备的安装容量约束：
0≤Pj,t≤Qj   (31)；
Figure BDA0003990825630000131
Figure BDA0003990825630000132
上式中，Pgrid表示联络线最大传输功率，Qj表示设备j的最大安装容量。
另外，电储能还需满足充放电前后的能量存储关系以及最大充放电功率和储能量约束，则通过引入辅助求解变量zs，将上述约束公式(31)～(33)转化为下述公式(34)和2个线性不等式(35)～(36)：
Figure BDA0003990825630000133
zs≥0   (35)；
Figure BDA0003990825630000134
S43、采用分支定界算法求解园区综合能源系统规划运行双层优化模型，得到不同置信水平下的规划及运行双层优化方案。
请参照图2至图15，本发明的实施例三为：
一种园区综合能源系统规划及运行优化方法，在上述实施例二的基础上，在本实施例中，以图2所示的典型PIES结构作为测试算例，由图2中所有设备构成设备候选集。该PIES采用“并网不上网”模式。为真实反映系统的实际运行情况，选取夏季、冬季和春秋季3个典型日数据来代表年度运行情况。各设备的经济和技术参数如下表1所示：
表1设备的经济与技术参数
Figure BDA0003990825630000135
Figure BDA0003990825630000141
PIES夏季、冬季及春秋季典型日的电/热/冷负荷曲线如图4所示。其中，折现率取值为6％，联络线最大传输功率为6MW，天然气价格采用固定气价，为2.55元/m3，天然气燃烧低热值为9.7kWh/m3，经过折算后的天然气价格为0.26元/kWh。假设园区采用分时电价，将一天分为谷时段(0:00-8:00)，峰时段(8:00-12:00，16:00-20:00)和平时段(12:00-16:00，20:00-24:00)，则平时段电价为0.8元/kWh，峰时段电价为1.6元/kWh，谷时段电价为0.4元/kWh。
不同能源种类的单位有功功率碳排放系数和政府分配的无偿碳配额系数如表2所示：
表2不同能源种类碳排放系数和碳配额系数
Figure BDA0003990825630000142
由表2可以看出，风机和光伏发电单位功率碳排放系数小于碳配额系数，即PIES在参与碳交易时，可通过风机和光伏发电在抵消系统自身碳排放量的同时，将多余的碳排放权在碳交易市场出售获利。
在本实施例中，根据碳交易价格历史数据(2017-2021)，利用神经网络技术对碳交易价格进行训练并预测，预测的碳交易价格均值μe近似为50元/吨，标准差σe为15元/吨。以华北某地区2020年全年风速和光照强度为例(数据采样间隔为1h)，标幺化后的数据如图5和图6所示。
根据全年实际数据，可估计出式(16)中参数ρ的取值为0.212，不同季节典型日在不同时刻标幺化后的光照强度均值μPV(ω,t)和风速均值μWT(ω,t)如图7所示。
统计出光照强度标准差近似为σPV(ω,t)＝0.25μPV(ω,t)，风速标准差近似为σWT(ω,t)＝0.25μWT(ω,t)。根据最大似然估计法可估计出式和式中的待估参数，估计结果如表3所示：
表3整合自回归滑动平均模型参数
Figure BDA0003990825630000151
进而利用重要性抽样法可分别得到500个碳交易价格、光照强度和风速场景。通过K-means聚类法，分别可得到碳交易价格、光照强度和风速的SSE与聚类场景数量的关系如图8、图9和图10所示。可以看出，碳交易价格、光照强度、风速的SSE随着聚类场景数量的增加而减小。当碳交易价格、光照强度和风速的聚类场景数量分别大于12、14和14时，SSE下降速度变得非常缓慢。因此最终选取的碳交易价格、光照强度和风速的典型场景数量分别为12、14和14个。聚类后的碳交易价格、光照强度和风速典型场景如图11、图12和图13所示。由碳交易价格、光照强度和风速典型场景构成的总场景数量N为(12×14×14)＝2352个，每个场景的概率为1/2352。
为了说明本实施例所提出的优化配置方法应对规划及运行优化总成本波动风险的有效性，本算例设置三种情形进行对比分析。
情形I：不考虑电储能，确定的碳交易价格与可再生能源出力；
情形II：不考虑电储能，碳交易价格和可再生能源出力不确定与时空相关性，θ设为90％；
情形III：考虑电储能，碳交易价格和可再生能源出力不确定与时空相关性，θ设为90％。
三种情形的规划和运行对比结果分别如表4和表5所示：
表4三种情形的规划方案对比结果
Figure BDA0003990825630000152
Figure BDA0003990825630000161
从表4可以看出，情形II的CHP、EB、PV和WT容量比情形I多配置了294kW、271kW、1424kWp和599kW，而GB的容量比情形I少配置了613kW。情形III的WT、EB和GB的容量比情形II多配置了584kW、249kW和56kW，而CHP和PV的容量分别比情形II少配置了674kW和5972kWp。下面将首先从经济性角度对规划方案的对比结果进行分析。
由表5可知，情形II的总成本比情形I高了153.78万元。下文从投资成本，运维成本和碳交易成本进行分析。三种情形下的冬季典型日购气量与购电量如图14(a)-(b)所示。
从图14(a)-(b)可知，在情形I和II下，由于系统没有配置电储能，情形II为应对碳交易价格与可再生能源出力不确定性带来的规划总成本波动风险，使得情形II比情形I需增大可再生能源配置容量来降低系统碳排放，同时需增加购气量使CHP多出力并增加外电网购电量来应对可再生能源波动，进而使得情形II的运维成本和碳交易成本增加了35.96万元和4.48万元。由于情形II的CHP、PV和WT的容量配置增加，使得投资成本相比于情形I增加了113.34万元，最终使得情形II总成本比情形I增加。但是情形II在优化配置时考虑了规划总成本波动风险，使得情形II的CVaR比情形I减少了135.12万元，CVaR值的大小反映系统面临规划总成本波动风险的高低，因此，相比于情形I，情形II降低了规划及运行总成本波动风险。通过对情形I和II的对比，说明当系统不配置电储能时，本发明所提规划及运行优化方法可以降低系统的规划及运行总成本波动风险，但是需要增加系统碳排放和总成本为代价。
另一方面，对比情形II，情形III总成本减少了162.93万元。原因是情形III配置了电储能可以应对可再生能源出力不确定性，使得系统不在配置CHP以减少系统购气量，同时不在从外电网购电，所以相比于情形II，情形III的购气量和外网购电量相对减少，如图14(a)-(b)所示。进而使得情形III的运维成本大幅降低，降幅约为120.64万元，远高于投资成本增幅29.06万元，而且情形III产生负碳排获利6.29万元有助于降低系统总成本，因此，情形III的总成本相比情形II减少。同时，由于配置电储能可以应对可再生能源出力不确定性，使得情形III的CVaR比情形II减少了62.57万元，情形III降低了规划及运行总成本波动风险。通过对情形II和III的对比，说明当系统配置电储能时，电储能不仅可以降低系统的规划及运行总成本波动风险而且可以实现系统的负碳化，即系统产生等效的负碳排放，获得碳交易收益。
进一步，以情形III为例，分析不同置信水平对规划及运行优化结果的影响，且对比了考虑风光时空相关性与否对规划结果的影响，不同置信水平下的规划方案对比结果如表6和表7所示，其中W1为不考虑风光时空相关性下的规划结果，W2为考虑风光时空相关性下的规划结果。
表6不同置信水平下的规划方案对比结果
Figure BDA0003990825630000171
表7不同置信水平下的经济性和碳交易量对比结果(单位:万元、吨)
Figure BDA0003990825630000172
可以看出，不论是否考虑风光时空相关性，随着置信水平的增加，系统越来越趋于规避规划总成本波动风险，系统的电储能配置容量增大来应对可再生能源出力不确定性，以致系统规划及运行总成本增加，但是由CVaR度量的规划及运行总成本波动风险降低。此外，当系统配置高比例可再生能源与一定容量电储能时，系统产生负碳排放，且随着置信水平的增加，系统的碳交易收益增大，系统主动利用碳交易价格的不确定性进行获利，但是需要增加一定的投资成本为代价。通过对比考虑风光时空相关性与否的规划结果可以看出，在同一置信水平下，考虑风光时空相关性下的电储能容量配置大于不考虑风光时空相关性情况。此外可以看出，在同一置信水平下，考虑风光时空相关性下的规划总成本高于不确考虑风光时空相关性情况，而CVaR低于不考虑风光时空相关性情况。主要原因是当系统考虑风光时空相关性时，需要配置更大的电储能应对成本波动风险，使得系统的投资成本增大以致规划及运行总成本增加，但是，系统可以规避高额的成本波动风险。
系统在不同置信水平下的总成本与CVaR的对比结果如图15所示。随着置信水平增加，系统的总成本增加而CVaR降低，且CVaR值降低的速率增。系统总成本和CVaR之间存在矛盾关系，PIES需根据自身的风险偏好选择对应的规划及运行优化方案。
请参照图16，本发明的实施例四为：
一种园区综合能源系统规划及运行优化终端1，包括存储器2、处理器3以及存储在存储器2是上并可在处理器3上执行的计算机程序，在本实施例中，处理器3执行计算机程序时实现上述实施例一至实施例三中任一实施例的步骤。
综上所述，本发明提供的一种园区综合能源系统规划及运行优化方法和终端，建立了基于能源集线器的多能流与碳交易量耦合模型来刻画碳交易量与多能流间的耦合关系；采用Copula函数和整合自回归滑动平均法建立了风光出力的时空相关性模型用于风光场景生成；此外，构建了CVaR风险度量的PIES规划模型应对规划总成本波动风险。通过算例验证，本发明可以有效降低系统面临的规划及运行总成本波动风险。
以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等同变换，或直接或间接运用在相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (10)

1.一种园区综合能源系统规划及运行优化方法，其特征在于，包括步骤：
S1、建立多能流与碳交易量耦合模型；
S2、根据所述多能流与碳交易量耦合模型，建立基于碳能耦合模型的园区综合能源系统模型；
S3、根据所述园区综合能源系统模型，建立考虑碳交易价格以及风光不确定性的场景生成模型；
S4、根据所述场景生成模型，建立基于条件风险价值的园区综合能源系统规划运行双层优化模型。
2.根据权利要求1所述的一种园区综合能源系统规划及运行优化方法，其特征在于，所述步骤S1具体为：
将碳交易量引入下述公式(1)中，建立多能流与碳交易量耦合模型：
Figure FDA0003990825620000011
可简化为公式(2)：
Figure FDA0003990825620000012
其中，αw(w＝1,2,…,m)为能源集线器中第w种输入能源的碳排放系数，
βw(w＝1,2,…,m)为能源集线器中第w种输入能源的无偿碳配额系数，
Figure FDA0003990825620000013
为能源集线器中第τ种能量存储的碳排放系数，ρτ(w＝1,2,…,u)为能源集线器中第τ种能量存储的无偿碳配额系数，，L为输出功率向量，C为转换设备耦合矩阵，P为输入功率向量，K为存储设备耦合矩阵，H为储能充放功率向量，E为能源集线器的碳交易量，A为1×m维度能源供给侧碳交易量系数矩阵，B为1×u维度的储能碳交易量系数矩阵，M为(n+1)×1维度的碳能输出矩阵，D为(n+1)×m维度的能源供给侧碳能耦合系数组合矩阵，U为(n+1)×u维度储能的碳能耦合系数组合矩阵；
其中，
Figure FDA0003990825620000021
3.根据权利要求2所述的一种园区综合能源系统规划及运行优化方法，其特征在于，所述步骤S2具体为：
S21、采用由热电联供机组、燃气锅炉、电制冷机、电锅炉、电储能、光伏以及风机设备组成的典型园区综合能源系统，并从系统外部购买天然气和电能，利用能源转换设备配合电储能来满足用户的电/热/冷负荷需求；
将所述典型园区综合能源系统参与碳交易出售盈余碳排放权进行获利；
所述典型园区综合能源系统的碳排放来源于一次能源供给端和能量存储设备产生的碳排放，即外电网购电、天然气消耗、光伏和风机发电以及电储能产生的碳排放，通过借鉴生命周期分析法得到归一化计量后的光伏、风机和电储能运行时的碳排放；
S22、对所述典型园区综合能源系统的多能流与碳交易量耦合关系进行建模，如下公式(3)～(5)：
Figure FDA0003990825620000022
Figure FDA0003990825620000023
Figure FDA0003990825620000024
其中，Pgrid、Pgas、PPV和PWT分别表示外网购电功率、天然气注入功率、光伏注入功率和风机注入功率，PEL、PCL、PHL和Etr分别表示电/冷/热负荷和碳交易量，c1表示供给电负荷的电功率分配系数，c2表示供给电锅炉的电功率分配系数，c3表示供给热电联供机组的天然气功率分配系数，
Figure FDA0003990825620000031
和
Figure FDA0003990825620000032
分别表示热电联供机组的电效率、热电联供机组的热效率、燃气锅炉的热效率、电锅炉的热效率和电制冷机的制冷效率，PES表示电储能的充/放电功率，当大于0时，表示电储能充电，当小于0时，表示电储能放电，αgrid、αPV、αWT、αgas和
Figure FDA0003990825620000033
分别表示外网购电、光伏、风机、天然气和电储能的单位功率碳排放系数，βgrid、βPV、βWT、βCHP、βGB和βES分别表示外网购电、光伏、风机、热电联供机组、燃气锅炉和电储能的单位功率无偿碳配额系数；
所述典型园区综合能源系统的碳交易量为外网购电、光伏、风机、天然气和电储能五部分的碳交易量之和，若其中某一部分的碳排放系数大于碳配额系数，则该部分的碳交易量为正，若某一部分的碳排放系数小于碳配额系数，则该部分的碳交易量为负，当为正值时，所述典型园区综合能源系统需要从碳交易市场购买碳排放权，当为负值时，所述典型园区综合能源系统可将多余碳排放权在碳交易市场出售；
S23、c1、c2和c3是为满足能源集线器建模形式要求引入的中间变量，可用各设备单元的输出功率和设备的效率来表示，即如下公式(6)～(8)：
Figure FDA0003990825620000034
Figure FDA0003990825620000035
Figure FDA0003990825620000036
其中，GCHP表示热电联供机组消耗的天然气功率，
Figure FDA0003990825620000037
表示电制冷机的输出功率，
Figure FDA0003990825620000038
表示电锅炉的热输出功率，
Figure FDA0003990825620000039
表示储能放电功率，
Figure FDA00039908256200000310
表示储能充电功率。
4.根据权利要求3所述的一种园区综合能源系统规划及运行优化方法，其特征在于，所述步骤S3具体为：
S31、对碳交易价格不确定性进行建模，具体为：
利用神经网络技术对碳交易价格进行训练并预测得到预测误差，并采用核密度估计方法得到所述预测误差的频率直方图，根据所述频率直方图确定所述预测误差服从零均值的正太分布，将所述预测误差与碳交易价格的预测值进行叠加，得到实际的碳交易交个概率分布，即如下公式(9)：
Figure FDA0003990825620000041
其中，xe表示碳交易价格随机变量，μe和σe分别表示碳交易价格均值和标准差；
S32、对考虑时空相关性的风光不确定性进行建模，包括风光空间相关性建模和风光时序相关性建模；
S32.A、采用Copula函数建立风光联合概率分布来描述风光空间相关性，则所述风光空间相关性建模的具体步骤如下：
基于光照强度和风速的历史数据，分别建立光照强度和风速的概率分布模型，光照强度的概率密度函数如下公式(10)～(12)所示：
Figure FDA0003990825620000042
Figure FDA0003990825620000043
Figure FDA0003990825620000044
其中，Γ(·)表示伽马函数，
Figure FDA00039908256200000413
表示在
Figure FDA0003990825620000045
季节t时刻的光照强度，
Figure FDA0003990825620000046
和
Figure FDA0003990825620000047
分别表示在
Figure FDA0003990825620000048
季节t时刻光照强度的形状参数，
Figure FDA0003990825620000049
和
Figure FDA00039908256200000410
分别表示在
Figure FDA00039908256200000411
季节t时刻光照强度的均值及标准差；
对于风速的不确定性，采用威布尔分布来描述，风速的概率密度函数如下公式(13)～(15)所示：
Figure FDA00039908256200000412
Figure FDA0003990825620000051
Figure FDA0003990825620000052
其中，
Figure FDA0003990825620000053
和
Figure FDA0003990825620000054
分别表示在
Figure FDA0003990825620000055
季节t时刻风速的尺度参数和形状参数，
Figure FDA00039908256200000519
表示在
Figure FDA0003990825620000056
季节t时刻的风速标准差，
Figure FDA0003990825620000057
表示在
Figure FDA0003990825620000058
季节t时刻的风速值，
Figure FDA0003990825620000059
表示在
Figure FDA00039908256200000510
季节t时刻的风速均值；
采用Frank-Copula函数对光照强度和风速的联合概率分布进行建模，风光联合概率分布表达式如下公式(16)所示：
Figure FDA00039908256200000511
其中，ρ表示相关系数，
Figure FDA00039908256200000512
F(·)表示光照强度的累积概率分布，
Figure FDA00039908256200000513
G(·)表示风速的累积概率分布；
S32.B、采用整合自回归滑动平均法来对光照强度、风速的时序相关性进行建模，则所述风光时序相关性建模如下公式(17)～(20)所示：
Figure FDA00039908256200000514
Figure FDA00039908256200000515
Figure FDA00039908256200000516
Figure FDA00039908256200000517
其中，
Figure FDA00039908256200000518
表示光照强度和风速的后移位算子，dPV和dWT分别表示光照强度和风速差分运算的阶数，εPV,t和εWT,t分别表示光照强度和风速的高斯白噪声，φPV,υ和φWT,υ分别表示光照强度和风速的第υ个自回归项系数，θPV,υ和θWT,υ分别表示光照强度和风速的第υ个滑动平均项系数，pPV和qPV分别表示光照强度自回归项的阶数和滑动平均项的阶数，pWT和qWT分别表示风速自回归项的阶数和滑动平均项的阶数，根据光照强度和风速的历史数据，采用极大似然估计法可估计出自回归项系数和滑动平均项系数以及对应的阶数；
S33、根据建立的风光联合概率分布式，通过重要性抽样法得到在某一时间断面下光照强度和风速的联合场景，具体为：
S33.a、在[0,1]区间内生成随机数a1和a2；
S33.b、令光照强度的概率分布函数值uPV＝a1，根据构建的风光联合概率分布来求得风速的概率分布函数值uWT＝a2，即求下述公式(21)的解：
Figure FDA0003990825620000061
S33.c、重复步骤S33.a～S33.b，得到一定组数的光照强度和风速的概率分布函数值；
S33.d、由于
Figure FDA0003990825620000062
利用反函数运算，即分别求解
Figure FDA0003990825620000063
和
Figure FDA0003990825620000064
将步骤S33.c所得的随机数样本转换为光照强度和风速的联合场景X＝{xPV,xWT}，通过重要性抽样法得到具有时序相关性的光照强度场景yPV和风速场景yWT；
S34、结合生成的风光时序场景yPV和yWT，对生成的具有空间相关性场景xPV和xWT进行修正，具体为：
S34.a、将风光空间相关性场景X中的风速场景xWT从小到大排序，排序后的风速场景记为
Figure FDA0003990825620000065
为保持相关关系的一致性，将与xWT对应的光照强度场景xPV跟随xWT变动，变动后光照强度场景记为
Figure FDA0003990825620000066
S34.b、将时序场景yPV和yWT进行转置，按列向量重构得到方阵y'PV和y'WT，取与xWT相同时间断面且具有时序相关性的风速场景y'WT，将其按照与xWT相同的大小顺序排列，排列后的场景记为
Figure FDA0003990825620000067
其余各时间断面场景跟随
Figure FDA0003990825620000068
变动，记为
Figure FDA0003990825620000069
根据下述公式(22)修正空间相关性场景中的各时间断面场景：
Figure FDA00039908256200000610
其中，Y'WT,t为修正后具有时空相关性的风速场景，同理可得修正后具有时空相关性的光照场景X'WT,t；
S35、通过K-means聚类法对生成的相似场景进行聚类，并利用误差平方和指标来确定碳交易价格、光照强度和风速的聚类场景数量，其中误差平方和表达式如下公式(23)所示：
Figure FDA0003990825620000071
其中，
Figure FDA0003990825620000072
表示第γ个随机变量的第a个类，
Figure FDA0003990825620000073
表示
Figure FDA0003990825620000074
中的样本数据，
Figure FDA0003990825620000075
表示第γ个随机变量第a个类的质心，
Figure FDA0003990825620000076
表示第γ个随机变量的聚类总数，根据公式(23)，可以分别确定碳交易价格、光照强度和风速的聚类场景数量；
S36、分别用向量Se、SPV和SWT表示碳交易价格、光照强度和风速的典型场景，由这三者的场景向量构建的典型场景集为S＝{Se,SPV,SWT}N，其中N表示总场景数量。
5.根据权利要求4所述的一种园区综合能源系统规划及运行优化方法，其特征在于，步骤S4具体为：
S41、建立基于条件风险价值的园区综合能源系统规划运行双层优化模型的目标函数：
根据条件风险价值的风险度量方法，基于条件风险价值的园区综合能源系统规划运行双层优化模型以全寿命周期规划及运行总成本的条件风险价值最小为目标，建立如下公式(24)～(30)的目标函数：
Figure FDA0003990825620000077
Figure FDA0003990825620000078
Figure FDA0003990825620000079
Figure FDA00039908256200000710
Figure FDA00039908256200000711
Figure FDA00039908256200000712
Figure FDA00039908256200000713
其中，
Figure FDA0003990825620000081
表示在置信水平θ下规划及运行总成本的风险值，该风险值指不超过给定置信水平的最大规划总成本的阈值，
Figure FDA0003990825620000082
表示在s场景下的规划及运行总成本，
Figure FDA0003990825620000083
表示在s场景下的等年值投资成本，
Figure FDA0003990825620000084
表示在s场景下的燃料成本，
Figure FDA0003990825620000085
表示在s场景下的维护成本，
Figure FDA0003990825620000086
表示在s场景下的碳交易成本，当
Figure FDA0003990825620000087
为负值时，表示卖出，园区综合能源系统获得收益，当
Figure FDA0003990825620000088
为正值时，表示买进，园区综合能源系统付出成本，
Figure FDA0003990825620000089
表示在s场景下的碳交易价格，r表示折现率，J表示设备候选集，lj表示第j个设备的寿命，
Figure FDA00039908256200000810
表示在s场景下第j个设备的容量，cinv,j表示第j个设备单位容量价格，cmain,j表示第j个设备单位功率的维护费用，
Figure FDA00039908256200000811
表示在s场景下第j个设备在t时刻的输出功率，T表示年总运行时间，为8760小时，
Figure FDA00039908256200000812
和
Figure FDA00039908256200000813
分别表示在s场景下t时刻系统的天然气购买量和从外网购电量，λgas和λelec分别表示天然气价格与电力价格，Δt表示系统运行时间间隔，取值为1小时，CVaRθ表示在置信水平下的条件风险价值，θ为置信水平，反映系统对规划总成本波动风险的厌恶水平；
S42、建立基于条件风险价值的园区综合能源系统规划运行双层优化模型的约束条件：
约束条件包括等式约束与不等式约束，等式约束为利用碳能耦合模型建立的系统功率平衡约束，不等式约束为设备运行时出力上下限约束以及设备的安装容量约束：
0≤Pj,t≤Qj   (31)；
Figure FDA00039908256200000814
Figure FDA00039908256200000815
其中，
Figure FDA00039908256200000816
表示联络线最大传输功率，
Figure FDA00039908256200000817
表示设备j的最大安装容量；
其中，电储能还需满足充放电前后的能量存储关系以及最大充放电功率和储能量约束，则通过引入辅助求解变量zs，将上述约束公式(31)～(33)转化为下述公式(34)和2个线性不等式(35)～(36)：
Figure FDA00039908256200000818
zs≥0   (35)；
Figure FDA0003990825620000091
S43、采用分支定界算法求解园区综合能源系统规划运行双层优化模型，得到不同置信水平下的规划及运行双层优化方案。
6.一种园区综合能源系统规划及运行优化终端，其特征在于，包括存储器、处理器和存储在存储器上并可在处理器上运行的计算机程序，所述处理器在执行所述计算机程序时实现以下步骤：
S1、建立多能流与碳交易量耦合模型；
S2、根据所述多能流与碳交易量耦合模型，建立基于碳能耦合模型的园区综合能源系统模型；
S3、根据所述园区综合能源系统模型，建立考虑碳交易价格以及风光不确定性的场景生成模型；
S4、根据所述场景生成模型，建立基于条件风险价值的园区综合能源系统规划运行双层优化模型。
7.根据权利要求6所述的一种园区综合能源系统规划及运行优化终端，其特征在于，所述步骤S1具体为：
将碳交易量引入下述公式(1)中，建立多能流与碳交易量耦合模型：
Figure FDA0003990825620000092
可简化为公式(2)：
Figure FDA0003990825620000101
其中，αw(w＝1,2,…,m)为能源集线器中第w种输入能源的碳排放系数，
βw(w＝1,2,…,m)为能源集线器中第w种输入能源的无偿碳配额系数，
Figure FDA0003990825620000102
为能源集线器中第τ种能量存储的碳排放系数，ρτ(w＝1,2,…,u)为能源集线器中第τ种能量存储的无偿碳配额系数，L为输出功率向量，C为转换设备耦合矩阵，P为输入功率向量，K为存储设备耦合矩阵，H为储能充放功率向量，E为能源集线器的碳交易量，A为1×m维度能源供给侧碳交易量系数矩阵，B为1×u维度的储能碳交易量系数矩阵，M为(n+1)×1维度的碳能输出矩阵，D为(n+1)×m维度的能源供给侧碳能耦合系数组合矩阵，U为(n+1)×u维度储能的碳能耦合系数组合矩阵；
其中，
Figure FDA0003990825620000103
8.根据权利要求7所述的一种园区综合能源系统规划及运行优化终端，其特征在于，所述步骤S2具体为：
S21、采用由热电联供机组、燃气锅炉、电制冷机、电锅炉、电储能、光伏以及风机设备组成的典型园区综合能源系统，并从系统外部购买天然气和电能，利用能源转换设备配合电储能来满足用户的电/热/冷负荷需求；
将所述典型园区综合能源系统参与碳交易出售盈余碳排放权进行获利；
所述典型园区综合能源系统的碳排放来源于一次能源供给端和能量存储设备产生的碳排放，即外电网购电、天然气消耗、光伏和风机发电以及电储能产生的碳排放，通过借鉴生命周期分析法得到归一化计量后的光伏、风机和电储能运行时的碳排放；
S22、对所述典型园区综合能源系统的多能流与碳交易量耦合关系进行建模，如下公式(3)～(5)：
Figure FDA0003990825620000111
Figure FDA0003990825620000112
Figure FDA0003990825620000113
其中，Pgrid、Pgas、PPV和PWT分别表示外网购电功率、天然气注入功率、光伏注入功率和风机注入功率，PEL、PCL、PHL和Etr分别表示电/冷/热负荷和碳交易量，c1表示供给电负荷的电功率分配系数，c2表示供给电锅炉的电功率分配系数，c3表示供给热电联供机组的天然气功率分配系数，
Figure FDA0003990825620000114
和
Figure FDA0003990825620000115
分别表示热电联供机组的电效率、热电联供机组的热效率、燃气锅炉的热效率、电锅炉的热效率和电制冷机的制冷效率，PES表示电储能的充/放电功率，当大于0时，表示电储能充电，当小于0时，表示电储能放电，αgrid、、αWT、αgas和
Figure FDA0003990825620000116
分别表示外网购电、光伏、风机、天然气和电储能的单位功率碳排放系数，βgrid、βPV、βWT、βCHP、βGB和βES分别表示外网购电、光伏、风机、热电联供机组、燃气锅炉和电储能的单位功率无偿碳配额系数；
所述典型园区综合能源系统的碳交易量为外网购电、光伏、风机、天然气和电储能五部分的碳交易量之和，若其中某一部分的碳排放系数大于碳配额系数，则该部分的碳交易量为正，若某一部分的碳排放系数小于碳配额系数，则该部分的碳交易量为负，当为正值时，所述典型园区综合能源系统需要从碳交易市场购买碳排放权，当为负值时，所述典型园区综合能源系统可将多余碳排放权在碳交易市场出售；
S23、c1、c2和c3是为满足能源集线器建模形式要求引入的中间变量，可用各设备单元的输出功率和设备的效率来表示，即如下公式(6)～(8)：
Figure FDA0003990825620000117
Figure FDA0003990825620000121
Figure FDA0003990825620000122
其中，GCHP表示热电联供机组消耗的天然气功率，
Figure FDA0003990825620000123
表示电制冷机的输出功率，
Figure FDA0003990825620000124
表示电锅炉的热输出功率，
Figure FDA0003990825620000125
表示储能放电功率，
Figure FDA0003990825620000126
表示储能充电功率。
9.根据权利要求8所述的一种园区综合能源系统规划及运行优化终端，其特征在于，所述步骤S3具体为：
S31、对碳交易价格不确定性进行建模，具体为：
利用神经网络技术对碳交易价格进行训练并预测得到预测误差，并采用核密度估计方法得到所述预测误差的频率直方图，根据所述频率直方图确定所述预测误差服从零均值的正太分布，将所述预测误差与碳交易价格的预测值进行叠加，得到实际的碳交易交个概率分布，即如下公式(9)：
Figure FDA0003990825620000127
其中，xe表示碳交易价格随机变量，μe和σe分别表示碳交易价格均值和标准差；
S32、对考虑时空相关性的风光不确定性进行建模，包括风光空间相关性建模和风光时序相关性建模；
S32.A、采用Copula函数建立风光联合概率分布来描述风光空间相关性，则所述风光空间相关性建模的具体步骤如下：
基于光照强度和风速的历史数据，分别建立光照强度和风速的概率分布模型，光照强度的概率密度函数如下公式(10)～(12)所示：
Figure FDA0003990825620000128
Figure FDA0003990825620000129
Figure FDA0003990825620000131
其中，Γ(·)表示伽马函数，
Figure FDA00039908256200001326
表示在
Figure FDA0003990825620000132
季节t时刻的光照强度，
Figure FDA0003990825620000133
和
Figure FDA0003990825620000134
分别表示在
Figure FDA0003990825620000135
季节t时刻光照强度的形状参数，
Figure FDA0003990825620000136
和
Figure FDA0003990825620000137
分别表示在
Figure FDA0003990825620000138
季节t时刻光照强度的均值及标准差；
对于风速的不确定性，采用威布尔分布来描述，风速的概率密度函数如下公式(13)～(15)所示：
Figure FDA0003990825620000139
Figure FDA00039908256200001310
Figure FDA00039908256200001311
其中，
Figure FDA00039908256200001312
和
Figure FDA00039908256200001313
分别表示在
Figure FDA00039908256200001314
季节t时刻风速的尺度参数和形状参数，
Figure FDA00039908256200001315
表示在
Figure FDA00039908256200001316
季节t时刻的风速标准差，
Figure FDA00039908256200001317
表示在
Figure FDA00039908256200001318
季节t时刻的风速值，
Figure FDA00039908256200001319
表示在
Figure FDA00039908256200001320
季节t时刻的风速均值；
采用Frank-Copula函数对光照强度和风速的联合概率分布进行建模，风光联合概率分布表达式如下公式(16)所示：
Figure FDA00039908256200001321
其中，ρ表示相关系数，
Figure FDA00039908256200001322
F(·)表示光照强度的累积概率分布，
Figure FDA00039908256200001323
G(·)表示风速的累积概率分布；
S32.B、采用整合自回归滑动平均法来对光照强度、风速的时序相关性进行建模，则所述风光时序相关性建模如下公式(17)～(20)所示：
Figure FDA00039908256200001324
Figure FDA00039908256200001325
Figure FDA0003990825620000141
Figure FDA0003990825620000142
其中，
Figure FDA0003990825620000143
表示光照强度和风速的后移位算子，dPV和dWT分别表示光照强度和风速差分运算的阶数，εPV,t和εWT,t分别表示光照强度和风速的高斯白噪声，φPV,υ和φWT,υ分别表示光照强度和风速的第υ个自回归项系数，θPV,υ和θWT,υ分别表示光照强度和风速的第υ个滑动平均项系数，pPV和qPV分别表示光照强度自回归项的阶数和滑动平均项的阶数，pWT和qWT分别表示风速自回归项的阶数和滑动平均项的阶数，根据光照强度和风速的历史数据，采用极大似然估计法可估计出自回归项系数和滑动平均项系数以及对应的阶数；
S33、根据建立的风光联合概率分布式，通过重要性抽样法得到在某一时间断面下光照强度和风速的联合场景，具体为：
S33.a、在[0,1]区间内生成随机数a1和a2；
S33.b、令光照强度的概率分布函数值uPV＝a1，根据构建的风光联合概率分布来求得风速的概率分布函数值uWT＝a2，即求下述公式(21)的解：
Figure FDA0003990825620000144
S33.c、重复步骤S33.a～S33.b，得到一定组数的光照强度和风速的概率分布函数值；
S33.d、由于
Figure FDA0003990825620000145
利用反函数运算，即分别求解
Figure FDA0003990825620000146
和
Figure FDA0003990825620000147
将步骤S33.c所得的随机数样本转换为光照强度和风速的联合场景X＝{xPV,xWT}，通过重要性抽样法得到具有时序相关性的光照强度场景yPV和风速场景yWT；
S34、结合生成的风光时序场景yPV和yWT，对生成的具有空间相关性场景xPV和xWT进行修正，具体为：
S34.a、将风光空间相关性场景X中的风速场景xWT从小到大排序，排序后的风速场景记为
Figure FDA0003990825620000148
为保持相关关系的一致性，将与xWT对应的光照强度场景xPV跟随xWT变动，变动后光照强度场景记为
Figure FDA0003990825620000149
S34.b、将时序场景yPV和yWT进行转置，按列向量重构得到方阵y'PV和y'WT，取与xWT相同时间断面且具有时序相关性的风速场景y'WT，将其按照与xWT相同的大小顺序排列，排列后的场景记为
Figure FDA0003990825620000151
其余各时间断面场景跟随
Figure FDA0003990825620000152
变动，记为
Figure FDA0003990825620000153
根据下述公式(22)修正空间相关性场景中的各时间断面场景：
Figure FDA0003990825620000154
其中，Y'WT,t为修正后具有时空相关性的风速场景，同理可得修正后具有时空相关性的光照场景X'WT,t；
S35、通过K-means聚类法对生成的相似场景进行聚类，并利用误差平方和指标来确定碳交易价格、光照强度和风速的聚类场景数量，其中误差平方和表达式如下公式(23)所示：
Figure FDA0003990825620000155
其中，
Figure FDA0003990825620000156
表示第γ个随机变量的第a个类，
Figure FDA0003990825620000157
表示
Figure FDA0003990825620000158
中的样本数据，
Figure FDA0003990825620000159
表示第γ个随机变量第a个类的质心，
Figure FDA00039908256200001510
表示第γ个随机变量的聚类总数，根据公式(23)，可以分别确定碳交易价格、光照强度和风速的聚类场景数量；
S36、分别用向量Se、SPV和SWT表示碳交易价格、光照强度和风速的典型场景，由这三者的场景向量构建的典型场景集为S＝{Se,SPV,SWT}N，其中N表示总场景数量。
10.根据权利要求9所述的一种园区综合能源系统规划及运行优化终端，其特征在于，步骤S4具体为：
S41、建立基于条件风险价值的园区综合能源系统规划运行双层优化模型的目标函数：
根据条件风险价值的风险度量方法，基于条件风险价值的园区综合能源系统规划运行双层优化模型以全寿命周期规划及运行总成本的条件风险价值最小为目标，建立如下公式(24)～(30)的目标函数：
Figure FDA00039908256200001511
Figure FDA00039908256200001512
Figure FDA0003990825620000161
Figure FDA0003990825620000162
Figure FDA0003990825620000163
Figure FDA0003990825620000164
Figure FDA0003990825620000165
其中，
Figure FDA0003990825620000166
表示在置信水平θ下规划及运行总成本的风险值，该风险值指不超过给定置信水平的最大规划总成本的阈值，
Figure FDA0003990825620000167
表示在s场景下的规划及运行总成本，
Figure FDA0003990825620000168
表示在s场景下的等年值投资成本，
Figure FDA0003990825620000169
表示在s场景下的燃料成本，
Figure FDA00039908256200001610
表示在s场景下的维护成本，
Figure FDA00039908256200001611
表示在s场景下的碳交易成本，当
Figure FDA00039908256200001612
为负值时，表示卖出，园区综合能源系统获得收益，当
Figure FDA00039908256200001613
为正值时，表示买进，园区综合能源系统付出成本，
Figure FDA00039908256200001614
表示在s场景下的碳交易价格，r表示折现率，J表示设备候选集，lj表示第j个设备的寿命，
Figure FDA00039908256200001615
表示在s场景下第j个设备的容量，cinv,j表示第j个设备单位容量价格，cmain,j表示第j个设备单位功率的维护费用，
Figure FDA00039908256200001616
表示在s场景下第j个设备在t时刻的输出功率，T表示年总运行时间，为8760小时，
Figure FDA00039908256200001617
和
Figure FDA00039908256200001618
分别表示在s场景下t时刻系统的天然气购买量和从外网购电量，λgas和λelec分别表示天然气价格与电力价格，Δt表示系统运行时间间隔，取值为1小时，CVaRθ表示在置信水平下的条件风险价值，θ为置信水平，反映系统对规划总成本波动风险的厌恶水平；
S42、建立基于条件风险价值的园区综合能源系统规划运行双层优化模型的约束条件：
约束条件包括等式约束与不等式约束，等式约束为利用碳能耦合模型建立的系统功率平衡约束，不等式约束为设备运行时出力上下限约束以及设备的安装容量约束：
0≤Pj,t≤Qj   (31)；
Figure FDA0003990825620000171
Figure FDA0003990825620000172
其中，
Figure FDA0003990825620000173
表示联络线最大传输功率，
Figure FDA0003990825620000174
表示设备j的最大安装容量；
其中，电储能还需满足充放电前后的能量存储关系以及最大充放电功率和储能量约束，则通过引入辅助求解变量zs，将上述约束公式(31)～(33)转化为下述公式(34)和2个线性不等式(35)～(36)：
Figure FDA0003990825620000175
zs≥0   (35)；
Figure FDA0003990825620000176
S43、采用分支定界算法求解园区综合能源系统规划运行双层优化模型，得到不同置信水平下的规划及运行双层优化方案。
CN202211589181.7A 2022-12-09 2022-12-09 一种园区综合能源系统规划及运行优化方法和终端 Pending CN115936220A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211589181.7A CN115936220A (zh) 2022-12-09 2022-12-09 一种园区综合能源系统规划及运行优化方法和终端

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211589181.7A CN115936220A (zh) 2022-12-09 2022-12-09 一种园区综合能源系统规划及运行优化方法和终端

Publications (1)

Publication Number Publication Date
CN115936220A true CN115936220A (zh) 2023-04-07

Family

ID=86651979

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211589181.7A Pending CN115936220A (zh) 2022-12-09 2022-12-09 一种园区综合能源系统规划及运行优化方法和终端

Country Status (1)

Country Link
CN (1) CN115936220A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117094745A (zh) * 2023-09-06 2023-11-21 天津大学 基于igdt-效用熵的综合能源系统优化控制方法及装置
CN117910025A (zh) * 2024-03-20 2024-04-19 神州数码融信云技术服务有限公司 一种金融服务数据安全存储保护方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117094745A (zh) * 2023-09-06 2023-11-21 天津大学 基于igdt-效用熵的综合能源系统优化控制方法及装置
CN117094745B (zh) * 2023-09-06 2024-03-12 天津大学 基于igdt-效用熵的综合能源系统优化控制方法及装置
CN117910025A (zh) * 2024-03-20 2024-04-19 神州数码融信云技术服务有限公司 一种金融服务数据安全存储保护方法
CN117910025B (zh) * 2024-03-20 2024-05-24 神州数码融信云技术服务有限公司 一种金融服务数据安全存储保护方法

Similar Documents

Publication Publication Date Title
Siqin et al. A two-stage distributionally robust optimization model for P2G-CCHP microgrid considering uncertainty and carbon emission
CN115936220A (zh) 一种园区综合能源系统规划及运行优化方法和终端
Srinivasan et al. An evolutionary algorithm for evaluation of emission compliance options in view of the clean air act amendments
CN112821465B (zh) 包含热电联产的工业微网负荷优化调度方法与系统
Qiu et al. Mean-field multi-agent reinforcement learning for peer-to-peer multi-energy trading
CN110796307A (zh) 一种综合能源系统分布式负荷预测方法及系统
Cao et al. Energy management optimisation using a combined Long Short-Term Memory recurrent neural network–Particle Swarm Optimisation model
CN107480907A (zh) 一种峰谷电价下含风电的省级电网购电配比的优化方法
CN115587668A (zh) 多园区综合能源系统分布式协同优化调度方法
Bagheri et al. Stochastic optimization and scenario generation for peak load shaving in Smart District microgrid: sizing and operation
CN116523683A (zh) 一种基于智能园区虚拟电厂的碳交易方法
CN116011821A (zh) 一种电力市场环境下虚拟电厂优化风险调度方法
CN112085259A (zh) 一种新能源接入与综合能源服务潜力关联预测方法
Liu et al. Planning mechanism design and benefit analysis of electric energy substitution: A case study of tobacco industry in Yunnan province, China
CN114169916A (zh) 一种适应新型电力系统的市场成员报价策略制定方法
Chen et al. Stochastic planning of integrated energy system based on correlation scenario generation method via Copula function considering multiple uncertainties in renewable energy sources and demands
Pfeifer et al. Fast energy transition as a best strategy for all? The nash equilibrium of long-term energy planning strategies in coupled power markets
CN112598175A (zh) 一种流域型虚拟电厂调度方法
CN108009684A (zh) 一种包含短期负荷预测的微电网并网状态能量管理方法
CN116882771A (zh) 主从博弈下考虑碳奖惩的零碳社区能量共享方法及系统
CN115759604B (zh) 一种综合能源系统优化调度方法
CN116258511A (zh) 一种考虑碳交易价格不确定性的能源系统配置方法及装置
Burgess et al. Electricity Demand and Supply on the Global Electrification Frontier
CN114925892A (zh) 水电与电转气结合的中长期风水火发电量双层规划方法
Ghasemi et al. Combating Uncertainties in Wind and Distributed PV Energy Sources Using Integrated Reinforcement Learning and Time-Series Forecasting

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination