技术领域
[0001] 本发明属于电力系统新能源消纳领域，尤其是含多类型电源的电力系统新能源消纳问题，具体涉及一种风光火储协调优化运行方法。
相关背景技术
[0002] 近年来，新能源发电发展迅速，新能源装机容量在系统中所占比重不断增大。但由于新能源发电具有波动性、随机性和间歇性的特点，新能源大规模并网给系统调节能力带来巨大挑战，造成火电机组频繁启停和弃风、弃光等问题。抽水蓄能机组作为一种储能设备，在系统削峰填谷、调频和黑启动方面发挥着重要作用，抽水蓄能机组与风力发电联合运行，已成为电力系统中解决新能源消纳问题重要的策略之一。
[0003]
[0004] 在我国北方地区，冬季供暖期热电联产机组“以热定电”的运行限制了机组的调峰能力，为新能源消纳带来问题。目前多电源并存的协调运行模型未考虑热电联产机组给电力系统带来的限制，而且现有的电力系统多源协调运行研究大多仅考虑风电机组与抽水蓄能电站协调运行，随着光伏电站装机容量的不断增大，光伏发电量在总发电量中所占的比重不断提高，因此有必要研究风电场、光伏电站、火电厂和抽水蓄能电站协调运行的电力系统运行方式。
具体实施方式
[0099] 下面结合附图对本发明的具体实施方式做进一步详细说明。
[0100] 参见图1，一种风光火储协调优化运行方法，包括状态采样、状态评估和结果统计三部分，具体分为以下步骤：
[0101] 1.输入原始数据，包括电力系统电负荷数据、热负荷数据、火电机组(仅供电)强迫停运率及平均修复时间、背压式热电联产机组强迫停运率及平均修复时间、抽汽式热电联产机组强迫停运率及平均修复时间、历史风速数据及历史太阳辐照度数据；定义最大模拟周期数Dmax和一个模拟周期内最大模拟时段数Tmax，且设定模拟周期D＝1，模拟时段T＝1；
[0102] 2.基于状态持续时间采样法获得模拟周期D内火电机组开机状态序列、背压式热电联产机组开机状态序列及抽汽式热电联产机组开机状态序列。假设常规火电机组、背压式热电联产机组及抽汽式热电联产机组的运行状态持续时间和故障状态持续时间均服从指数分布，具体步骤为：
[0103] 2.1)基于状态持续时间采样法获得火电机组开机状态序列
[0104] 2.1.1)根据步骤1中火电机组强迫停运率及平均修复时间，计算火电机组平均无故障工作时间，计算公式为：
[0105]
[0106] 式中：MTTFth为火电机组平均无故障工作时间，FORth为火电机组强迫停运率，MTTRth为火电机组平均修复时间；
[0107] 2.1.2)定义序列 表示模拟周期D内火电机组开机状态序列， 值为1时说明模拟周期D的模拟时段T内常规火电机组处于运行状态， 值为0时说明模拟周期D的模拟时段T内常规火电机组处于停运状态；定义模拟周期D内常规火电机组已确定开机状态数 且设定
[0108] 2.1.3)生成(0,1)区间内服从均匀分布的随机数 和
[0109] 2.1.4)计算火电机组运行持续时间 及故障持续时间
[0110] 火电机组运行持续时间计算公式为：
[0111]
[0112] 火电机组故障持续时间计算公式为：
[0113]
[0114] 2.1.5)根据火电机组运行持续时间及故障持续时间确定 取值，公式为：
[0115]
[0116] 式中： 为向下取整符号；
[0117] 2.1.6)更新 若 说明模拟周期D的全部模拟时段内火电机组的开机状态均确定；否则，返回步骤2.1.3)，继续确定模拟周期D的其他模拟时段内火电机组开机状态；
[0118] 2.2)基于状态持续时间采样法获得背压式热电联产机组开机状态序列[0119] 2.2.1)根据步骤1中背压式热电联产机组强迫停运率及平均修复时间，计算背压式热电联产机组平均无故障工作时间，计算公式为：
[0120]
[0121] 式中： 为背压式热电联产机组平均无故障工作时间， 为背压式热电联产机组强迫停运率， 为背压式热电联产机组平均修复时间；
[0122] 2.2.2)定义序列 表示模拟周期D内背压式热电联产机组开机状态序列， 值为1时说明模拟周期D的模拟时段T内背压式热电联产机组处于运行状态， 值为0时说明模拟周期D的模拟时段T内背压式热电联产机组处于停运状态；定义模拟周期D内背压式热电联产机组已确定开机状态数 且设定
[0123] 2.2.3)生成(0,1)区间内服从均匀分布的随机数 和
[0124] 2.2.4)计算背压式热电联产机组运行持续时间 及故障持续时间
[0125] 背压式热电联产机组运行持续时间计算公式为：
[0126]
[0127] 背压式热电联产机组故障持续时间计算公式为：
[0128]
[0129] 2.2.5)根据背压式热电联产机组运行持续时间及故障持续时间确定 取值，公式为：
[0130]
[0131] 式中： 为向下取整符号。
[0132] 2.2.6)更新 若 说明模拟周期D的全部模拟时段内背压式热电联产机组的开机状态均确定；否则，返回步骤2.2.3)，继续确定模拟周期D的其他模拟时段内背压式热电联产机组开机状态；
[0133] 2.3)基于状态持续时间采样法获得抽汽式热电联产机组开机状态序列[0134] 2.3.1)根据步骤1中抽汽式热电联产机组强迫停运率及平均修复时间，计算抽汽式热电联产机组平均无故障工作时间，计算公式为：
[0135]
[0136] 式中： 为抽汽式热电联产机组平均无故障工作时间， 为抽汽式热电联产机组强迫停运率， 为抽汽式热电联产机组平均修复时间；
[0137] 2.3.2)定义序列 表示模拟周期D内抽汽式热电联产机组开机状态序列， 值为1时说明模拟周期D的模拟时段T内抽汽式热电联产机组处于运行状态， 值为0时说明模拟周期D的模拟时段T内抽汽式热电联产机组处于停运状态；定义模拟周期D内抽汽式热电联产机组已确定开机状态数 且设定
[0138] 2.3.3)生成(0,1)区间内服从均匀分布的随机数 和
[0139] 2.3.4)计算抽汽式热电联产机组运行持续时间 及故障持续时间
[0140] 抽汽式热电联产机组运行持续时间计算公式为：
[0141]
[0142] 抽汽式热电联产机组故障持续时间计算公式为：
[0143]
[0144] 2.3.5)根据抽汽式热电联产机组运行持续时间及故障持续时间确定 取值，公式为：
[0145]
[0146] 式中： 为向下取整符号。
[0147] 2.3.6)更新 若 说明模拟周期D的全部模拟时段内抽汽式热电联产机组的开机状态均确定；否则，返回步骤2.3.3)，继续确定模拟周期D的其他模拟时段内抽汽式热电联产机组开机状态。
[0148] 3.基于逆变换法获得模拟周期D内风电机组可用功率序列及光伏机组可用功率序列；假设模拟风速服从两参数威布尔分布，两参数威布尔分布的分布函数为：
[0149]
[0150] 式中：v为模拟风速，aw为两参数威布尔分布的尺度参数，bw为两参数威布尔分布的形状参数；
[0151] 假设模拟太阳辐照度服从β分布，β分布的分布函数为：
[0152]
[0153] 式中：r为模拟太阳辐照度，rmax为历史太阳辐照度最大值，as为β分布的位置参数，bs为β分布的形状参数；
[0154] 基于逆变换法获得模拟周期D内风电机组可用功率序列及光伏机组可用功率序列的具体步骤为：
[0155] 3.1)获得模拟周期D内风电机组可用功率序列
[0156] 3.1.1)根据步骤1中的历史风速数据，计算威布尔分布的形状参数和尺度参数；
[0157] 威布尔分布的形状参数计算公式为：
[0158]
[0159] 式中：sw为历史风速数据的标准差，为历史风速数据的平均值；
[0160] 威布尔分布的尺度参数计算公式为：
[0161]
[0162] 式中：Γ(·)为Gama函数；
[0163] 3 .1 .2 ) 生 成 ( 0 ,1 ) 区 间 内 服 从 均 匀 分 布 的 随 机 数 序 列[0164] 3.1.3)定义序列 表示模拟周期D内的模拟风速序列，利用步骤3.1.2)中生成的随机数序列计算模拟风速序列，计算公式为：
[0165]
[0166] 3.1.4)定义序列 表示模拟周期D内的风电机组可用功率序列；根据模拟风速和风电机组可用功率的关系，计算模拟周期D的模拟时段T内风电机组可用功率 模拟风速与风电机组可用功率之间的关系可用二次分段函数表示，具体计算公式为：
[0167]
[0168] 式中：PWn为风电机组额定功率，vin为风电机组切入风速；vout为风电机组切出风速；vn为风电机组额定风速； 和 为风电机组功率曲线参数，均由vin和vn确定；
[0169] 3.2)获得模拟周期D内光伏机组可用功率序列
[0170] 3.2.1)根据步骤1中的历史太阳辐照度数据，计算β分布的位置参数和形状参数。β分布的位置参数计算公式为：
[0171]
[0172] 式中：为历史太阳辐照度数据的平均值，ss为历史太阳辐照度数据的方差；
[0173] β分布的形状参数计算公式为：
[0174]
[0175] 3.2.2)生成服从β(as,bs)分布的随机数序列
[0176] 3.2.3)定义序列sD＝{sD,T,T＝1,2,…,Tmax}表示模拟周期D内的模拟太阳辐照度序列，利用步骤3.2.2)中生成的随机数序列计算模拟太阳辐照度序列，计算公式为：
[0177]
[0178] 3.2.4)定义序列 表示模拟周期D内的光伏机组可用功率序列。根据模拟太阳辐照度和光伏机组可用功率的关系，计算模拟周期D内模拟时段T内光伏机组可用功率 模拟太阳辐照度与光伏机组可用功率之间的关系可用二次分段函数表示，具体计算公式为：
[0179]
[0180] 式中：Psn为光伏机组额定功率，sin为光伏机组切入太阳辐照度；sout为光伏机组切出太阳辐照度；sn为光伏机组额定太阳辐照度； 和 为光伏机组功率曲线参数，均由sin和sn确定。
[0181] 4.获取模拟周期D的模拟时段T内的电力系统状态，包括电负荷数、热负荷数、火电机组开机状态、背压式热电联产机组开机状态及抽汽式热电联产机组开机状态、风电机组可用功率及光伏机组可用功率。
[0182] 5.对步骤4中的电力系统状态求解风光火储协调运行模型，从而对模拟周期D的模拟时段T内常规火电机组功率、背压式热电联产机组供电功率及供热功率、抽汽式热电联产机组供电功率及供热功率、抽水蓄能机组抽水功率或发电功率进行优化，风光火储协调运行模型属于优化模型中的混合整数非线性规划模型，能够利用分支定界法进行求解。
[0183] 若模型存在可行解，则说明模拟周期D的模拟时段T内电力系统可完全消纳风电机组可用功率及光伏机组可用功率，即风电机组实际消纳功率等于风电机组可用功率，光伏机组实际消纳功率等于光伏机组可用功率，进入步骤8；若模型不存在可行解，说明模拟周期D的模拟时段T内电力系统不能完全消纳风电机组可用功率及光伏机组可用功率，进入步骤6；风光火储协调运行模型由目标函数和约束条件组成，具体模型为：
[0184] 5.1)目标函数
[0185] 风光火储协调运行模型的目标函数为整个模拟周期内电力系统运行成本最小，电力系统运行成本包括常规火电机组运行成本、背压式热电联产机组运行成本、抽汽式热电联产机组运行成本及抽水蓄能机组运行成本，具体为：
[0186]
[0187] 火电机组i在模拟时段T内的运行成本为：
[0188]
[0189] 背压式热电联产机组m在模拟时段T内的运行成本仅与其供电功率有关，具体计算公式为：
[0190]
[0191] 抽汽式热电联产机组n在模拟时段T内的运行成本不仅与其供电功率有关，还与其供热功率有关，具体计算公式为：
[0192]
[0193] 抽水蓄能电站运行成本为
[0194]
[0195] 式中：Nth为火电机组数， 为背压式热电联产机组数， 为抽汽式热电联产机组数，Np为抽水蓄能电站数， 为火电机组i在模拟时段T内的功率， 为背压式热电联产机组m在模拟时段T内的供电功率， 为抽汽式热电联产机组n在模拟时段T内的供电功率， 为抽汽式热电联产机组n在模拟时段T内的供热功率，ai,bi和ci为火电机组耗量特性参数， 为背压式热电联产机组耗量特性参数，为抽汽式热电联产机组耗量特性参数， 为抽水蓄能电站u启
机一次的成本， 为抽水蓄能电站u停机一次的成本， 为抽水蓄能电站u在模拟时段T内的抽水运行状态指示变量， 为抽水蓄能电站u在模拟时段T+1内的发电运行状态指示变量；
[0196] 5.2)约束条件
[0197] 风光火储协调运行模型的约束条件包括电力系统实时功率平衡约束、机组功率约束、火电机组爬坡约束、热电联产机组供热约束、抽水蓄能电站约束及电力系统备用约束，具体如下：
[0198] 5.2.1)电力系统实时功率平衡约束
[0199] 为保证电力系统的稳定运行，需考虑电力系统实时功率平衡约束，电力系统实时功率平衡约束包括实时电功率平衡约束和实时热功率平衡约束；
[0200] 实时电功率平衡约束为：
[0201]
[0202] 实时热功率平衡约束为：
[0203]
[0204] 式中： 为抽水蓄能电站u在模拟时段T内的抽水功率或发电功率， 为风电机组j在模拟时段T内的可用功率， 为光伏机组k在模拟时段T内的可用功率，PL,T为电力系统在模拟时段T内的电负荷，HL,T为电力系统在模拟时段T内的热负荷，Nw为风电机组数，Ns为光伏机组数， 为背压式热电联产机组m在模拟时段T内的供热功率；
[0205] 5.2.2)机组功率约束
[0206] 机组功率约束包括火电机组功率约束、背压式热电联产机组供电功率约束、抽汽式热电联产机组供电功率约束、抽汽式热电联产机组供热功率约束及抽水蓄能电站功率约束；
[0207] 火电机组功率约束为：
[0208]
[0209] 背压式热电联产机组供电功率约束为：
[0210]
[0211] 抽汽式热电联产机组供电功率约束为：
[0212]
[0213] 抽汽式热电联产机组供热功率约束为：
[0214]
[0215] 抽水蓄能电站功率约束为：
[0216]
[0217] 式中： 为火电机组i的最小功率， 为火电机组i的最大功率， 为背压式热电联产机组m的最小供电功率， 为背压式热电联产机组m的最大供电功率，为抽汽式热电联产机组n的最小供电功率， 为抽汽式热电联产机组n的最大供电功率，Hnmin为抽汽式热电联产机组n的最小供热功率，Hnmax为抽汽式热电联产机组n的最大供热功率， 为抽水蓄能电站u的最小抽水功率， 为抽水蓄能电站u的最大抽水功率，为抽水蓄能电站u的最小发电功率， 为抽水蓄能电站u的最大发电功率；
[0218] 5.2.3)火电机组爬坡约束
[0219]
[0220]
[0221] 式中：Rup,i为火电机组i的向上爬坡率，Rdown,i为火电机组i的向下爬坡率，Δt为每个模拟时段步长；
[0222] 5.2.4)热电联产机组供热约束
[0223] 热电联产机组包括背压式热电联产机组和抽汽式热电联产机组，背压式热电联产机组m的供热约束为线性等式约束，具体为：
[0224]
[0225] 抽汽式热电联产机组n的供热约束为线性不等式约束，具体为：
[0226]
[0227] 式中： 和 分别为背压式热电联产机组电/热线性关系参数，由背压式热电联产机组m的运行工况图确定， 及 分别为抽汽式热电联产机组电/热线性关系参数，由抽汽式热电联产机组n的运行工况图确定；
[0228] 5.2.5)抽水蓄能电站约束
[0229] 抽水蓄能电站约束包括上水库库容约束、抽水蓄能电站能量守恒约束、抽水蓄能电站末状态约束及经济性约束；
[0230] 上水库库容约束为：
[0231]
[0232] 式中：Eumin为抽水蓄能电站u的上水库最小储能，Eumax为抽水蓄能电站u的上水库最大储能， 为抽水蓄能电站u在模拟时段T结束时的上水库储能，计算公式为：
[0233]
[0234] 式中：Eu,0为抽水蓄能电站u上水库的初始储能， 为抽水蓄能电站u在抽水运行状态下的能量转换效率， 为抽水蓄能电站u在发电运行状态下的能量转换效率；
[0235] 抽水蓄能电站能量守恒约束为：
[0236]
[0237] 抽水蓄能电站末状态约束为：
[0238]
[0239] 经济性约束为：
[0240] v＝1,2,…,Np，w＝1,2,…,Np且v≠w
[0241] 5.2.6)系统备用约束
[0242]
[0243] 式中：R为电力系统备用容量。
[0244] 6.引入风电机组弃风功率及光伏机组弃光功率，对步骤5中的风光火储协调运行模型进行修改；
[0245] 6.1)修改步骤5.1)中的目标函数。
[0246] 为保证风电机组可用功率和光伏机组可用功率尽可能多的消纳，在电力系统运行成本项中加入弃风成本和弃光成本，修改后的风光火储协调运行模型的目标函数为：
[0247]
[0248] 式中： 为弃风惩罚系数，其大小由消纳风电机组可用功率的迫切程度确定，为弃光惩罚系数，其大小由消纳光伏机组可用功率的迫切程度确定， 为风电机组j在模拟时段T内的弃风功率， 为光伏机组k在模拟时段T内的弃光功率；
[0249] 6.2)在步骤5.2)约束条件的基础上增加风电机组弃风功率约束与光伏机组弃光功率约束。风电机组j弃风功率约束为：
[0250]
[0251] 光伏机组k弃光功率约束为：
[0252]
[0253] 6.3)在步骤5.2)约束条件的基础上增加风电机组可用功率、风电机组实际消纳功率与风电机组弃风功率之间的关系约束及光伏机组可用功率、光伏机组实际消纳功率与光伏机组弃光功率之间的关系约束；
[0254] 风电机组j可用功率、风电机组j实际消纳功率与风电机组j弃风功率之间的关系约束为：
[0255]
[0256] 光伏机组k可用功率、光伏机组k实际消纳功率与光伏机组k弃光功率之间的关系约束为：
[0257]
[0258] 其中： 为模拟时段T内风电机组j实际消纳功率， 为模拟时段T内光伏机组k实际消纳功率；
[0259] 6.4)修改步骤5.2.1)中电力系统实时电功率平衡约束
[0260] 修改后的电力系统实时电功率平衡约束为：
[0261]
[0262] 7.对步骤4中的电力系统状态求解步骤6中修改后的风光火储协调运行模型，从而对模拟周期D的模拟时段T内火电机组功率、背压式热电联产机组供电功率及供热功率、抽汽式热电联产机组供电功率及供热功率、抽水蓄能机组抽水功率或发电功率、风电机组实际消纳功率及光伏机组实际消纳功率进行优化，修改后的风光火储协调运行模型属于优化模型中的混合整数非线性规划模型，能够利用分支定界法进行求解。
[0263] 8.统计模拟周期D的模拟时段T内火电机组功率、背压式热电联产机组供电功率及供热功率、抽汽式热电联产机组供电功率及供热功率、抽水蓄能机组抽水功率或发电功率、风电机组实际消纳功率及光伏机组实际消纳功率。
[0264] 9.如果模拟周期D的模拟时段T等于最大模拟时段数Tmax，则进入步骤10，否则，返回步骤4，进入下一模拟时段。
[0265] 10.如果模拟周期D等于最大模拟周期数Dmax，则进入步骤11，否则，返回步骤2，进入下一模拟周期。
[0266] 11.输出各模拟周期内火电机组功率、背压式热电联产机组供电功率及供热功率、抽汽式热电联产机组供电功率及供热功率、抽水蓄能机组抽水功率及发电功率、风电机组实际消纳功率、光伏机组实际消纳功率。
