考虑CHP机组快速爬坡能力约束的电热耦合系统动态优化调度方法

考虑chp机组快速爬坡能力约束的电热耦合系统动态优化调度方法
技术领域
1.本发明涉及电热耦合系统优化调度技术领域，具体涉及一种考虑chp机组快速爬坡能力约束的电热耦合系统动态优化调度方法。


背景技术：

2.因生产工艺需要，大型工业企业常配备热电联产(chp)电厂来供应全部热力需求，并服从电网调度参与区域电力供应。但企业自备chp机组刚性电出力严重挤占新能源消纳空间，阻碍了企业响应碳配额政策。近年来，电热耦合系统(cehs)通过chp机组灵活性改造提升电出力快速爬坡能力，但改造后的机组特性及其cehs优化调度愈加复杂，因此考虑chp机组快速爬坡能力的cehs动态优化调度方法研究，对支撑被改造chp机组运行、提升系统调度灵活性具有重要意义。
3.根据chp机组内在特性，采用优化控制方法提高机组短期快速调节能力，是一种技术成本低、综合优势大的适用于企业自备chp机组灵活性改造的方案。该方案可在短时间内显著改变机组电出力的爬坡能力而又不影响企业热负荷供应。在chp机组微分-代数方程(dae)机理模型上采用优化控制的精确化动态建模，机组爬坡速率的提升可起到补偿风电扰动、适应电网快速调频等作用，且chp机组快速爬坡能力可与风电快速波动形成互补，充分利用机组调度灵活性将有利于提高cehs风电消纳能力。但现有研究较少将chp机组快速爬坡能力精确化模型应用到cehs优化调度中，深入分析其对风电消纳的作用。
4.目前chp机组灵活性改造主要通过低压缸灵活切除、加装储热、配置电热锅炉、高中压缸旁路供热等直接解耦方案，虽可显著增加机组的灵活性但改造成本较高。经改造后的机组特性常采用代数方程约束建模参与cehs调度优化计算，该基于非精确的代数方程建模，无法刻画chp机组爬坡调节过程的电热耦合、压力安全、供热稳定等动态非线性特性。
5.考虑chp机组快速爬坡能力精确化动态建模的cehs优化调度，是一个含dae约束的动态优化问题，其高非线性、非凸、非光滑特性导致传统的线性规划直接求解方法难以适用。大多数dae优化求解均是将其离散转化为非线性规划(nlp)问题后，再用各种nlp求解器或智能优化算法进行优化求解，但其离散后产生的大规模问题，常需要特殊的分解技术和复杂的数学处理。在兼顾求解效率与精度时，现有cehs动态优化调度问题求解均面临求解策略与离散方法适用性的挑战。


技术实现要素：

6.为解决上述技术问题，本发明提供一种考虑chp机组快速爬坡能力约束的电热耦合系统动态优化调度方法，在保证调度方案的可行性和安全可靠性的前提下，有效促进风力发电实现最大化消纳，从而提升灵活性改造的chp机组参与cehs调度的有效性，且具备成本低的优点。
7.本发明采取的技术方案为：
8.考虑chp机组快速爬坡能力约束的电热耦合系统动态优化调度方法，包括以下步骤：
9.步骤1：建立chp机组快速爬坡能力约束模型；
10.步骤2：考虑chp机组快速爬坡能力约束，构建电热耦合系统cehs优化调度模型；
11.步骤3：设计cehs双层动态优化调度策略，采用基于径向基函数的拟序贯法，求解内层chp机组控制过程动态优化问题；采用动态自适应粒子群算法，生成外层系统调度稳态优化方案。
12.所述步骤1中，chp机组快速爬坡能力约束模型包括：
13.chp机组电热耦合动态关系的数学模型：
[0014][0015]
其中：rm为实际进入磨煤机煤量，ub为机组燃料量，t为控制时域时间变量，ε为制粉过程迟延时间；ts为制粉惯性时间，rb为锅炉燃烧率；cb为锅炉蓄热系数，pd为汽包压力，k3为汽轮机增益静态参数，p
t
为汽轮机前压力，u
t
为汽轮机高压缸进汽调节阀开度，k1为燃料增益静态参数；k2为过热器阻力系数；t
t
为汽轮机惯性时间，ph为机组发电功率，k4为高中压缸占汽轮机做功比例，k5为低压缸增益静态参数，pz为汽轮机供热抽汽压力，uh为供热抽汽调节蝶阀开度；ch为机组加热器蓄热系数，k6为热网循环水的有效比热容，τh为机组循环水回水温度；mh为供热抽汽流量，k7为供热抽汽有效热量折合蒸汽流量系数，m
x
为机组循环水流量；p1为汽轮机一级压力。
[0016]
快速爬坡能力优化控制算法模型：
[0017]
分别设计汽轮机u
t
(t)、锅炉ub(t)和供热uh(t)三个变量pid控制器，实现电功率准确跟踪、热功率快速恢复和系统运行安全稳定三个关键功能；
[0018][0019]
其中：k
pt
为控制u
t
(t)变量的比例系数，k
it
为控制u
t
(t)变量的积分系数，k
dt
为控制u
t
(t)变量的微分系数；k
pb
为控制ub(t)变量的比例系数，k
ib
为控制ub(t)变量的积分系数，k
db
为控制ub(t)变量的微分系数；k
ph
为控制uh(t)变量的比例系数，k
ih
为控制uh(t)变量的积
分系数，k
dh
为控制uh(t)变量的微分系数；
[0020]
控制优化周期内，为主蒸汽设定值，为机组发电功率设定值，为供热抽汽流量设定值；p
t
(t)为主蒸汽压力控制过程当前值，ph(t)为机组发电功率控制过程当前值，mh(t)为供热抽汽流量控制过程当前值。
[0021]
考虑到系统安全稳定，调节过程需考虑输出变量波动的影响，满足如式(4)路径约束；同时，各输出变量稳态值不能超出允许的误差范围，满足如式(5)终值约束。
[0022][0023]
其中：m
pt
为主蒸汽压力的波动范围，m
ph
为电功率的波动范围，m
mh
为抽汽蒸汽流量的波动范围；
[0024][0025]
其中：δ
pt
为主蒸汽压力的误差范围，δ
ph
为电功率的误差范围，δ
mh
为抽汽流量的误差范围；p
t
(te)为主蒸汽压力控制过程终值，ph(te)为机组发电功率控制过程终值，mh(te)为供热抽汽流量控制过程终值。
[0026]
热力系统热波动模型：
[0027]
抽汽蒸汽流量波动引起的热量变化对热负荷需求的影响，按如下方法计算：
[0028]qload
＝m
l
·cw
·
(τ
out-τ
in
)
ꢀꢀ
(6)；
[0029]
其中：q
load
为热负荷需求，m
l
为热网循环水流量，cw为水的比热容，τ
out
热网送水温度，τ
in
热网回水温度；
[0030]qh,equ
＝m
h,equ
·
δh
ꢀꢀ
(7)；
[0031]
其中：q
h,equ
为chp机组热功率，δh为供热抽汽焓降，m
h,equ
为调节过程的平均等效流量；
[0032][0033]
其中：t0和te分别为机组优化控制的起始和终止时间点，mh(t)为供热抽汽流量控制过程值；
[0034][0035]
其中：c
sys
为供热区域热容，τ
indoor
为供热区域温度，dτ
indoor
/dt为温度变化速率。所述步骤2中，电热耦合系统cehs优化调度模型包括：
[0036]
(1)、目标函数：
[0037]
优化调度以cehs运行成本最小为目标：
[0038]
[0039]
式(10)中，为第t时段第i台火电机组出力；为第t时段第j台chp机组电，为第t时段第j台chp热出力；为第t时段第l台风电机组出力。
[0040]
式(10)中，f
ig
、f
lw
函数形式如下：
[0041][0042]
其中，f
ig
为第i台火电机组运行成本，a
2,i
、a
1,i
、a
0,i
分别为第i台火电机组的二次项成本系数、一次项成本系数、常项系数；为第j台chp机组运行成本，b
2,j
、b
1,j
、b
0,j
分别为第j台chp机组电出力二次项成本系数、一次项成本系数、常项系数，c
2,j
、c
1,j
分别为第j台chp机组热出力二次项成本系数、一次项成本系数，c
0,j
为第j台chp机组电热耦合项成本系数；为第t时段第l台风电机组实际功率，f
lw
为第l台风电机组惩罚成本，θ为风电弃风惩罚成本系数；
[0043]
(2)、约束条件：
[0044]
1)火电机组运行约束包括：机组出力约束、机组爬坡约束；
[0045]
风电机组运行约束；
[0046]
cehs电功率平衡与热功率平衡约束；
[0047]
电力系统网络安全约束；
[0048][0049]
其中：为第t时段第i台火电机组出力；为第t时段第i台火电机组出力最小值、为第t时段第i台火电机组出力最大值；
[0050][0051]
其中：r
g,i
为第i台火电机组电爬坡速率；为第t时段第i台火电机组上备用，δt为系统调度间隔时间；
[0052][0053]
其中：为第t时段第i台火电机组下备用；
[0054][0055]
其中：为第t-1时段第i台火电机组出力；
[0056][0057]
其中：为第t时段第l台风电机组出力、为第t时段第l台风电机组实际功率；
[0058]
[0059]
其中：i为火电机组数量变量，j为chp机组数量变量，l为风电机组数量变量，为第t时段电负荷；ωg为火电机组集合，ωh为chp机组集合，ωw为风电机组集合；
[0060][0061]
其中：为第t时段热负荷；为第t时段第j台chp热出力；
[0062][0063]
其中：为支路s时段t的传输功率；分别为支路s最大传输功率下界、上界，ωb为电力系统中的所有传输线路；
[0064]
2)：chp机组快速爬坡能力动态约束，写成dae模型的标准形式：
[0065][0066]
式(20)中：u为chp机组控制变量组成的向量，y为chp机组代数状态变量组成的向量，x为chp机组微分状态变量组成的向量，p为与t无关的待优化参数，t0和te分别为机组优化控制的起始和终止时间点；x
l
、xu分别为微分状态变量上限、下限组成的向量，y
l
、yu分别为代数状态变量上限、下限组成的向量，u
l
、uu分别为控制变量上限、下限组成的向量，f为chp机组微分代数方程动态模型，对应式(2)；g(t)为chp机组系统路径约束，对应式(3)-(4)；h(te)为chp机组系统在时间te的终值约束，对应式(5)。
[0067]
dae模型中：
[0068][0069]
式(21)中：x1，x2，x3，x4为chp机组模型的4个微分状态变量，分别对应pd为汽包压力，ph为机组发电功率，pz为汽轮机供热抽汽压力，rb为锅炉燃烧率；u1，u2，u3为chp机组模型的3个控制变量，分别对应u
t
为汽轮机高压缸进汽调节阀开度，ub为机组燃料量，uh为供热抽汽调节蝶阀开度；y1，y2，y3为chp机组模型的4个代数状态变量，分别对应p
t
为汽轮机前压力，ph为机组发电功率，mh为供热抽汽流量；
[0070]
3)计算chp机组热量波动及温度变化情况，其波动约束换算成热量变化约束：
[0071][0072][0073]
其中，δq为调度时段内(5min)允许热量变化值；δq为调度周期内(2h)允许热量
变化值；m为cehs系统调度时段数。
[0074]
所述步骤3中，cehs双层动态优化调度策略设计如下:
[0075]
①
、外层优化模型：
[0076][0077]
式(24)中：j1为调度模型的目标函数，x＝[x1,x2,x3,x4]为优化调度变量组成的决策向量，分别为火电机组电出力、风电机组电出力以及chp机组电出力、热出力；d(x)为调度优化决策变量的约束。
[0078]
②
、内层优化模型：
[0079][0080]
式(25)中：u为chp机组控制变量组成的向量，y为chp机组代数状态变量组成的向量，x为chp机组微分状态变量组成的向量，p为与t无关的待优化参数，t0和te分别为机组优化控制的起始和终止时间点；x0为微分状态变量的初始值；x
l
、xu分别为微分状态变量上限、下限组成的向量，y
l
、yu分别为代数状态变量上限、下限组成的向量，u
l
、uu分别为控制变量上限、下限组成的向量，f为chp机组微分代数方程；g(t)为chp机组系统路径约束；h(te)为chp机组系统在时间te的终值约束。x1，x2，x3，x4为chp机组模型的4个微分状态变量，分别对应pd为汽包压力，ph为机组发电功率，pz为汽轮机供热抽汽压力，rb为锅炉燃烧率；q、r和s是正定权矩阵，可通过对变量做归一化来确定初始值。
[0081]
式中，
[0082]
x3(t)＝x2(t)
ꢀꢀ
(26)；
[0083]
其中，x3为调度模型决策的chp机组电出力，x2为chp机组模型优化控制的chp机组电出力；
[0084]
x4(t)＝δh
·
(k7k6m
x
(96x3(t)-τh 103))
ꢀꢀ
(27)；
[0085]
其中，x4为调度模型决策的chp机组热出力，x3为chp机组模型优化控制的汽轮机供热抽汽压力，δh为供热抽汽焓降，m
x
为机组循环水流量，τh为机组循环水回水温度，k6为热网循环水的有效比热容，k7为供热抽汽有效热量折合蒸汽流量系数；
[0086]
所述步骤3中，外层调度稳态优化算法如下：
[0087]
动态自适应粒子群算法dapso在应用于含dae约束的cehs动态优化调度问题中，其数学描述为：
[0088][0089]
式(28)中：为第i个粒子在k次迭代时的速度；为第i个粒子在k次迭代时的位置；为第k 1次迭代第i个粒子运动速度，为第k 1次迭代第i个粒子的当前位置，n为种群数量；c1和c2为加速系数，一般取c1＝c2＝2；r1和r2为2个在[0,1]范围内变化的随机数；p
ik
为第i个粒子在k次迭代中搜索到的最好位置；为整个群体在k次迭代中搜索到的最好位置；为第i个粒子在k次迭代时的惯性权重；α、β在[0,1]内选择，一般取α＝β＝0.5；a、b为控制聚集度因子的阈值，一般取a＝0.9，b＝0.5；为进化速度因子；s1为适应值聚集度因子；s2为空间聚集度因子。
[0090][0091]
式(29)中：p
ik
为第i个粒子在k次迭代中搜索到的最好位置，p
ik-1
为第i个粒子在k-1次迭代中搜索到的最好位置，f(p
ik
)为p
ik
的适应值，f(p
ik-1
)为p
ik-1
的适应值，fk为第k次迭代中的最佳适应值，为第k次迭代中的平均适应值，n为种群规模，l为搜索空间的最长半径，d为解空间的维数，g
id
为第i个粒子的第d维坐标，为所有粒子第d维坐标的平均值。
[0092]
基于dapso的外层调度稳态优化算法步骤：
[0093]
1)根据cehs动态优化问题设置dapso参数，生成调度方案的决策变量；
[0094]
2)加载cehs调度模型的数据，决策变量代入非chp机组的代数方程约束进行检验；
[0095]
3)决策变量输入chp机组dae模型进行过程动态优化计算，检验决策变量是否满足机组控制特性约束；
[0096]
4)若调度方案的决策变量不满足步骤2)、步骤3)的任一约束条件时，决策变量乘以较大惩罚系数后输出计算调度模型目标函数值；
[0097]
5)若调度方案的决策变量满足全部约束条件，则输出计算目标函数值，比较函数最小值保留；
[0098]
6)迭代次数k＝k 1，利用式(28)、式(29)更新下一次的种群调度方案，直到满足设置的迭代次数终止条件。
[0099]
所述步骤3中，基于rbf拟序贯法的内层过程动态优化方法包括如下步骤：
[0100]
步骤1)、给定初值：
[0101]
给定控制变量初值u0、状态变量初值x0、优化时间t0和te及所有变量的边界约束。
[0102]
步骤2)、变量离散：
[0103]
离散式(20)的dae模型是将各连续变量沿时间轴进行离散，因而只需考虑用一元数值逼近方法来进行离散。采用rbf实施离散，离散后的nlp模型具有更简洁的公式描述，更易实现编程求解。先对式(20)所示导函数进行逼近，然后通过积分得到原函数x(t)的逼近。
[0104]
对于给定chp机组dae初值问题：
[0105]
x(t0)＝x0(30)；
[0106]
其中，x为chp机组微分状态变量，f为chp机组微分代数方程，t0为机组优化控制的起始时间点，x0为chp机组微分状态变量的初始值；
[0107]
已知{(xi,fi)|i＝1,2,
…
,n}，对于任意给定x，可用式(31)进行插值，且满足式(32)条件：
[0108][0109]
其中，ζi为第i个采样点的实系数，φ(
·
)为一元基函数，i为采样点，n为控制周期采样点总数，xi为第i个采样点的微分状态变量值；
[0110][0111]
其中，fk为第k个采样点微分代数方程函数值，k为采样点，xk为第k个采样点的微分状态变量值；
[0112]
将chp机组优化控制时间离散为时间序列{t0,t1,
…
,tn}，可用下式对式(31)进行逼近。
[0113][0114]
其中，ζi为第i个采样点的实系数，φ(
·
)为一元基函数，ti为优化控制的第i个采样点的时间点；
[0115]
对进行不定积分，可得：
[0116][0117]
其中，ζi为第i个采样点的实系数，φ(
·
)为一元基函数，c0为常数；
[0118]
令则可得：
[0119][0120]
其中，ti为优化控制的第i个采样点的时间点，为φ(
·
)积分后的函数；
[0121]
易得初值问题的离散格式：
[0122]
[0123]
其中，ζi为第i个采样点的实系数，φ(
·
)为一元基函数，为φ(
·
)积分后的函数，c0为常数，ti为优化控制的第i个采样点的时间点，j为采样点，tj为优化控制的第j个采样点的时间点，x0为微分状态变量的初始值；
[0124]
联立求解式(36)得和c0的值，代入式(35)可得离散点处变量值计算中c0可通过求差方式消去。将(x(ti))记为xi，则简化式(36)为：
[0125][0126]
其中，ζk为第k个采样点的实系数，ti为优化控制的第i个采样点的时间点；tj为优化控制的第j个采样点的时间点；t
j-1
为优化控制的第j-1个采样点的时间点；tk为优化控制的第k个采样点的时间点,φ(
·
)为一元基函数，为φ(
·
)积分后的函数；
[0127]
基函数φ(
·
)常采用rbf形式，其中，逆multi-quadric函数最适合求解dae，设c＝γte/n，对于大多数dae优化模型，当γ＝2可使求解过程易于收敛。即：
[0128][0129]
其中，c为给定正数，r为中间变量，φ(
·
)为一元基函数；
[0130]
对φ(
·
)进行不定积分可得：
[0131][0132]
其中，c为给定正数，ti为优化控制的第i个采样点的时间点；为φ(
·
)积分后的函数；
[0133]
式(25)所示的动态优化问题，所有变量都是时间t∈[t0,te]的函数。不失一般性，本发明假定被优化的dae对各变量有连续二阶偏导数。仅考虑式(25)的离散化，对优化控制时段[t0,te]进行等距划分，可得时间序列且有tn＝te，t
i-t
i-1
＝te/n。将各时间点ti所对应各变量分别记为xi、yi和ui，按式(37)对式(25)进行离散，可得如下nlp模型：
[0134][0135]
式中：s为向量的变量维数，ui为第i个采样点控制变量，yi为第i个采样点代数状态变量，x
s,i
为第i个采样点第s维微分状态变量，p为与t无关的待优化参数，x
l
、xu分别为微分
状态变量上限、下限，y
l
、yu分别为代数状态变量上限、下限，u
l
、uu分别为控制变量上限、下限，f为微分代数方程动态模型；g(.)为路径约束；h(.)为终值约束。ti为优化控制的第i个采样点的时间点；tj为优化控制的第j个采样点的时间点；t
j-1
为优化控制的第j-1个采样点的时间点；tk为优化控制的第k个采样点的时间点；x
s,i
为第i个采样点第s维微分状态变量，ζ
s,k
为第k个采样点第s维变量的实系数，为一元基函数φ(
·
)积分后的函数；x
s,0
为初值采样点第s维微分状态变量值；x
s,n
为终止采样点第s维微分状态变量值；un为终止采样点控制变量值，yn为终止采样点代数状态变量值；x
s,j
为第j个采样点第s维微分状态变量，x
s,j-1
为第j-1个采样点第s维微分状态变量；
[0136]
步骤3)、模型计算：
[0137]
在每个优化迭代中都要进行一次dae模型方程的模拟，在模拟层中对离散后的非线性代数方程组采用牛顿迭代法进行模型计算得到状态变量的值。
[0138]
步骤4)、灵敏度计算：
[0139]
在模型计算的同时，得到状态变量对控制变量和时间变量的灵敏度信息矩阵。
[0140]
步骤5)、优化层计算：
[0141]
每次优化迭代的过程中，重复步骤4)和步骤5)，利用现有的nlp求解器，求解经过模拟层处理后的非线性规划问题，直至达到优化迭代的收敛精度，寻得最优解。
[0142]
本发明一种考虑chp机组快速爬坡能力约束的电热耦合系统动态优化调度方法，技术效果如下：
[0143]
1)对带优化控制的chp机组快速爬坡动态模型进行数值仿真，与静态模型相比，动态模型具有更快的爬坡速率和更大的爬坡容量；
[0144]
2)从提高综合用能系统风电消纳能力的角度，充分考虑了chp机组快速爬坡能力在cehs优化调度中对风电消纳的作用，改善了被电热耦合特性限制的chp机组刚性电出力严重挤占新能源消纳空间的不足现状；
[0145]
3)在保证调度方案的可行性和安全可靠性的前提下，有效促进风力发电实现最大化消纳，从而提升灵活性改造的chp机组参与cehs调度的有效性。
[0146]
4)本发明方法考虑chp机组快速爬坡能力的优化调度，可在保证cehs调度方案可行性的前提下促进风电消纳、节约运行成本。
附图说明
[0147]
图1为含自备chp机组的cehs结构图。
[0148]
图2为chp机组快速爬坡能力动态模型结构图。
[0149]
图3(1)为静态模型与动态模型的chp机组电出力对比图；
[0150]
图3(2)为静态模型与动态模型的chp机组主蒸汽压力对比图；
[0151]
图3(3)为静态模型与动态模型的chp机组抽汽蒸汽流量对比图。
[0152]
图4为cehs动态优化调度模型求解流程图。
[0153]
图5为拟序贯算法结构图；
[0154]
图6为cehs结构图。
[0155]
图7为不同算例下优化调度结果图。
[0156]
图8为不同算例下chp机组运行点分布图；
[0157]
图9为不同算例下chp机组电出力图；
[0158]
图10为chp机组主蒸汽压力曲线图。
[0159]
图11为不同算例下chp机组抽汽蒸汽流量图。
具体实施方式
[0160]
一种考虑chp机组快速爬坡能力的cehs动态优化调度方法。首先构建chp机组机理模型，提出优化控制策略，分析变负荷扰动对机组运行特性的影响。然后形成chp机组快速爬坡能力的动态约束，在此基础上建立含微分代数方程约束的cehs动态优化调度模型。最后，采用外层基于dapso的调度稳态优化、内层基于径向基拟序贯的过程动态优化相结合的双层优化结构进行动态问题求解。具体实现过程如下：
[0161]
步骤1：建立chp机组快速爬坡能力约束模型：
[0162]
研究的cehs由热力系统、电力系统、chp机组供能与控制系统3个部分组成，结构如图1所示。chp机组中ub控制锅炉燃烧产生高温蒸汽，通过调节蝶阀u
t
控制蒸汽进入高压缸；中压缸与低压缸连接管道内安装抽汽调节蝶阀uh，控制抽汽蒸汽供给热源输出热功率，同时控制蒸汽进入低压缸做功输出电功率。热源与供热首站交换热量，并通过换热装置供应工业企业内部热负荷。由于企业区域面积较小，热管网在空间上分布较为有限，管道热量传输过程影响可以忽略。同时，chp机组承担自备电厂企业的生产自用电以及供电辐射区域生产生活的部分电负荷，供电不足部分由大网调度风力发电机组、常规火电机组来满足。
[0163]
在传统电网cehs优化调度中，chp机组常采用简化的代数方程约束模型(也即静态模型)：
[0164][0165]
其中，t为调度时段；δt为爬坡时间；为机组上备用容量、下备用容量；为机组t时段允许的最大、最小电出力；为chp机组向上、向下爬坡速率，其爬坡率一般为2％pe/min，pe为机组额定功率，经过优化控制系统改造后的chp机组爬坡率可提高到4％pe/min以上。
[0166]
式(1)所示的静态模型无法表征chp机组调节过程中压力、流量、温度等变量变化的非线性特性。采用静态模型考虑机组快速爬坡能力时，在机组电出力爬坡上可达到与动态模型类似的效果，但静态模型无法约束机组压力、流量等变量变化引起的机组短时爬坡缺额、安全性、稳定性状态，单一的chp机组爬坡约束过于宽松时将导致实际运行调节难以实施或引发安全事故。
[0167]
抽汽式chp机组快速爬坡能力约束分别建立chp机组电热耦合动态关系的微分代数方程、快速爬坡能力优化控制算法模型、热力系统热波动模型。
[0168]
(1)chp机组电热耦合动态关系的数学模型：
[0169][0170]
其中：ub为机组燃料量；u
t
为汽轮机高压缸进汽调节阀开度；uh为供热抽汽调节蝶阀开度；m
x
为机组循环水流量；τh为机组循环水回水温度；p
t
为汽轮机前压力；ph为机组发电功率；pz为汽轮机供热抽汽压力；mh为供热抽汽流量；p1为汽轮机一级压力；rm为实际进入磨煤机煤量；rb为锅炉燃烧率。k1、k2、k3、k4、k5、k6、k7为静态参数；ε制粉过程迟延时间，ts制粉惯性时间，cb锅炉蓄热系数，t
t
汽轮机惯性时间，ch机组加热器蓄热系数。
[0171]
(2)快速爬坡能力优化控制算法模型：
[0172]
针对上述模型采用优化控制方法实现chp机组快速爬坡能力调节，增强机组的短期灵活性。分别设计汽轮机u
t
(t)、锅炉ub(t)和供热uh(t)三个变量pid控制器，实现电功率准确跟踪、热功率快速恢复和系统运行安全稳定三个关键功能：
[0173][0174]
其中，k
pt
、k
it
、k
dt
，k
pb
、k
ib
、k
db
，k
ph
、k
ih
、k
dh
为u
t
(t)、ub(t)、uh(t)pid参数。控制优化周期内，为主蒸汽设定值，为机组发电功率设定值，为供热抽汽流量设定值。
[0175]
如图2所示，在t时段以机组主蒸汽压力抽汽蒸汽流量为定值输入，以相邻调度时段机组电出力形成阶跃信号为变负荷扰动输入，与chp机组模型共同构建快速爬坡能力动态约束模型。其中，主蒸汽压力根据机组电、热出力范围设定，机组热出力与抽汽蒸汽流量关系如式(7)所示。
[0176]
考虑到系统安全稳定，调节过程需考虑输出变量波动的影响，应满足如下路径约束：
[0177][0178]
同时，各输出变量稳态值不能超出允许的误差范围，还应满足如下终值约束：
[0179][0180]
其中，m
pt
为主蒸汽压力的波动范围，取设定值
±
5％；，m
ph
为电功率的波动范围，取1mw；m
mh
为抽汽蒸汽流量的波动范围，取50t/h；δ
pt
为主蒸汽压力的误差范围，取设定值
±
2.5％，δ
ph
为电功率的误差范围，取0.1mw；δ
mh
为抽汽蒸汽流量的误差范围，取5t/h。
[0181]
(3)热力系统热波动模型：
[0182]
将chp机组快速爬坡能力动态模型进行建模仿真，给定变负荷扰动的阶跃信号，分析在优化控制周期内主蒸汽压力、抽汽蒸汽流量以及机组电功率的变化过程。抽汽蒸汽流量波动引起的热量变化对热负荷需求的影响，按如下方法计算：
[0183]qload
＝m
l
·cw
·
(τ
out-τ
in
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)；
[0184]qh,equ
＝m
h,equ
·
δh
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)；
[0185][0186][0187]
其中，q
load
为热负荷需求，m
l
为热网循环水流量，cw为水的比热容，τ
out
热网送水温度，τ
in
热网回水温度；q
h,equ
为chp机组热功率，δh为供热抽汽焓降，mh(t)为流量调节过程值，m
h,equ
为调节过程的平均等效流量，t0和te为机组优化控制的起始和终止时间点；c
sys
为供热区域热容，τ
indoor
为供热区域温度，dτ
indoor
/dt为温度变化速率。
[0188]
采用不考虑快速爬坡能力时的静态模型，爬坡率取2％pe/min。考虑快速爬坡能力时的静态模型(爬坡率取4.5％pe/min)以及动态模型分别进行chp机组仿真计算，如图3(1)～图3(3)所示，以运行点为例，三种建模方式变负荷能力分别可达结果表明，考虑快速爬坡能力时chp机组具有更快的爬坡速率和更大的爬坡容量。动态模型的主蒸汽压力、抽汽蒸汽流量经波动后均回到稳态，其中主蒸汽压力最大波动15.73在设定值15的
±
5％波动范围内，可保证机组运行安全；300s内等效平均流量为393.19t/h，比设定值400t/h减少约4471.33kj热量，换算成温度变化为4.21
×
10-3
℃，证明蒸汽流量波动引起的温度变化对热负荷影响不大。将快速爬坡能力静态建模的形成阶跃扰动输入chp机组动态模型进行检验，数据表明主蒸汽压力波峰15.82超过最大上限，机组运行存在安全隐患，且蒸汽流量波动也增大，不利于机组稳定供热。
[0189]
步骤2：构造考虑chp机组快速爬坡能力约束的cehs优化调度模型：
[0190]
(1)目标函数：
[0191]
本发明优化调度以cehs运行成本最小为目标：
[0192][0193]
式中，f
ig
、f
lw
函数形式如下：
[0194][0195]
其中，fg为火电机组运行成本，成本系数a2、a1、a0；pg为火电机组出力；fh为chp机组运行成本，发电发热成本系数b2、b1、b0、c2、c1、c0；ph、qh为chp机组电、热出力；fw为风电惩罚成本，成本系数θ；pw为风电机组出力，为实际风电功率；ωg、ωh、ωw分别为火电机组、chp机组、风电机组集合。决策变量为调度周期2小时，间隔δt为5min，调度时段t∈m(m＝1,2,
…
,24)。
[0196]
(2)约束条件：
[0197]
1)火电机组运行约束，包括机组出力约束、机组爬坡约束；风电机组运行约束；cehs电功率平衡与热功率平衡约束；电力系统网络安全约束。
[0198][0199][0200][0201][0202][0203][0204][0205][0206]
其中，rg为火电机组电爬坡速率；分别为火电机组出力最小值、最大值；为火电机组上备用、下备用；p
load
电负荷需求；q
load
热负荷需求；为支路s时段t的传输功率；分别为支路s最大传输功率下界、上界，ωb为电力系统中的所有传输线路。
[0207]
2)chp机组快速爬坡能力动态约束，如式(2)-(5)所示。将其写成dae模型的标准形式：
[0208][0209]
模型中u为控制变量，y为代数状态变量，x为微分状态变量，其中：
[0210][0211]
式中：f为chp机组微分代数方程动态模型，对应式(2)；g(t)为系统路径约束，对应式(3)-(4)；h(te)为系统在时间te的终值约束，对应式(5)。不等式约束可通过添加大于零的松弛变量转换成等式约束增加到系统。
[0212]
3)热力系统通过式(6)-(9)计算chp机组热量波动及温度变化情况。根据生产工艺、生活供暖需求调度时段内(5min)的温度变化小于0.5℃，2小时之内的温度变化小于1～2℃，可由式(9)换算成热量变化约束
[0213][0214][0215]
其中，δq为调度时段内(5min)允许热量变化值；δq为调度周期内(2h)允许热量变化值。
[0216]
步骤3：设计电热耦合系统双层动态优化调度策略进行分层求解，采用基于径向基函数的拟序贯法求解内层chp机组控制过程动态优化问题，利用动态自适应粒子群算法生成外层系统调度稳态优化方案：
[0217]
1)外层优化模型：
[0218][0219]
式中：j1为调度模型的目标函数对应式(10)-(11)，x＝[x1,x2,x3,x4]为优化调度变量组成的决策向量，分别为火电、风电机组电出力以及chp机组电、热出力；d(x)为调度优化决策变量的约束，对应式(12)-(23)。
[0220]
2)内层优化模型
[0221][0222]
式中，
[0223]
x3(t)＝x2(t)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(26)；
[0224]
x4(t)＝δh
·
(k7k6m
x
(96x3(t)-τh 103))
ꢀꢀ
(27)；
[0225]
其中，u
*
、y
*
分别为稳态优化得到的控制变量与输出变量的设定值，u
*
同时也是变负荷动态过程的前馈设定，u(t)是第t个时刻的控制变量，y(t)为第t-1个时刻控制对象的测量输出，作为第t个时刻闭环系统的反馈，q、r和s是正定权矩阵,可通过对变量做归一化来确定初始值。
[0226]
目标函数j1调度决策变量实现cehs经济目标，采用智能算法给定调度方案及迭代寻优；目标函数j2引导控制变量、输出变量逼近期望的稳态值，使机组控制系统加快收敛速度，并使得操纵变量尽量平滑，采用拟序贯法进行求解。
[0227]
(1)外层调度稳态优化算法：
[0228]
考虑到离散后的cehs优化调度是一个复杂的非线性动态优化问题，可将高效处理多种约束寻优的粒子群算法作为优化计算层的寻优算法。与其他优化算法相比，动态自适应粒子群算法(dapso)在应用于本发明含dae约束的cehs动态优化调度问题中，表现出了相对出色的收敛能力和求解精度，其数学描述为：
[0229][0230]
式中：为第i个粒子在k次迭代时的速度；为第i个粒子在k次迭代时的位置；c1和c2为加速系数，一般取c1＝c2＝2；r1和r2为2个在[0，1]范围内变化的随机数；p
ik
为第i个粒子在k次迭代中搜索到的最好位置；为整个群体在k次迭代中搜索到的最好位置；为第i个粒子在k次迭代时的惯性权重；α、β在[0，1]内选择，一般取α＝β＝0.5；a、b为控制聚集度因子的阈值，一般取a＝0.9，b＝0.5；为进化速度因子；s1为适应值聚集度因子；s2为空间聚集度因子。
[0231][0232]
式中：f(p
ik
)为p
ik
的适应值；fk为第k次迭代中的最佳适应值；为第k次迭代中的平均适应值；n为种群规模；l为搜索空间的最长半径；d为解空间的维数；g
id
为第i个粒子的第d维坐标；为所有粒子第d维坐标的平均值。
[0233]
cehs双层动态优化调度策略流程如图4所示，基于dapso的外层调度稳态优化算法步骤：
[0234]
1)根据cehs动态优化问题设置dapso参数，生成调度方案的决策变量；
[0235]
2)加载cehs调度模型的数据，决策变量代入非chp机组的代数方程约束进行检验；
[0236]
3)决策变量输入chp机组dae模型进行过程动态优化计算，检验决策变量是否满足机组控制特性约束；
[0237]
4)若调度方案的决策变量不满足步骤2)、步骤3)的任一约束条件时，决策变量乘以较大惩罚系数后输出计算调度模型目标函数值；
[0238]
5)若调度方案的决策变量满足全部约束条件，则输出计算目标函数值，比较函数最小值保留；
[0239]
6)迭代次数k＝k 1，利用式(28)、式(29)更新下一次的种群调度方案，直到满足设置的迭代次数终止条件。
[0240]
(2)内层过程动态优化方法：
[0241]
拟序贯算法结合序贯算法和联立算法的优点，适合求解带路径约束的大规模强非线性动态优化问题，分为模拟计算层和优化计算层。在模拟计算层中，控制变量被表示为每一时间单元内的常量函数或线性函数，状态变量采用以径向基函数(rbf)为基的正交配置法离散，从而将非线性daes离散为非线性代数方程组，可通过牛顿迭代法求解，得到离散点上的状态变量值。同时，根据离散点上状态变量关于控制变量的导数公式，可计算出灵敏度矩阵供优化计算层获取目标函数对于控制变量的梯度信息。控制变量分段函数值成为优化问题的决策变量，约束条件不再含daes等式约束，状态变量仅以路径约束形式在不等式约束中出现，原问题简化为nlp问题，可直接应用序列二次规划算法求解。拟序贯算法结构如图5所示。
[0242]
基于rbf拟序贯法的内层过程动态优化方法步骤：
[0243]
1)给定初值：
[0244]
给定控制变量初值u0、状态变量初值x0、优化时间t0和te及所有变量的边界约束。
[0245]
2)变量离散：
[0246]
离散式(20)的dae模型是将各连续变量沿时间轴进行离散，因而只需考虑用一元
数值逼近方法来进行离散。采用rbf实施离散，离散后的nlp模型具有更简洁的公式描述，更易实现编程求解。先对式(20)所示导函数进行逼近，然后通过积分得到原函数x(t)的逼近。
[0247]
对于给定chp机组dae初值问题：
[0248]
x(t0)＝x0(30)；
[0249]
已知{(xi,fi)|i＝1,2,
…
,n}，对于任意给定x，可用式(30)进行插值，且满足式(31)条件：
[0250][0251][0252]
其中，为实系数，φ(
·
)为一元函数(基函数)，为rbf插值的中心。
[0253]
将chp机组优化控制时间离散为时间序列{t0,t1,
…
,tn}，可用下式对(30)进行逼近。
[0254][0255]
对进行不定积分，可得
[0256][0257]
令则可得：
[0258][0259]
易得初值问题的离散格式：
[0260][0261]
联立求解式(36)得和c0的值，代入式(35)可得离散点处变量值计算中c0可通过求差方式消去。将(x(ti))记为xi，则简化式(36)为：
[0262][0263]
基函数φ(
·
)常采用rbf形式，其中逆multi-quadric函数最适合求解dae，c为给定正数，设c＝γte/n，对于大多数dae优化模型，当γ＝2可使求解过程易于收敛。即
[0264][0265]
对φ(
·
)进行不定积分可得：
[0266][0267]
式(25)所示的动态优化问题，所有变量都是时间t∈[t0,te]的函数。不失一般性，本发明假定被优化的dae对各变量有连续二阶偏导数。仅考虑式(25)的离散化，对优化控制时段[t0,te]进行等距划分，可得时间序列且有tn＝te，t
i-t
i-1
＝te/n。将各时间点ti所对应各变量分别记为xi、yi和ui，按式(37)对式(25)进行离散，可得如下nlp模型：
[0268][0269]
3)模型计算：
[0270]
在每个优化迭代中都要进行一次dae模型方程的模拟，在模拟层中对离散后的非线性代数方程组采用牛顿迭代法进行模型计算得到状态变量的值。
[0271]
4)灵敏度计算：
[0272]
在模型计算的同时，得到状态变量对控制变量和时间变量的灵敏度信息矩阵。
[0273]
5)优化层计算：
[0274]
每次优化迭代的过程中，重复步骤4)和步骤5)，利用现有的nlp求解器求解经过模拟层处理后的非线性规划问题，直至达到优化迭代的收敛精度，寻得最优解。
[0275]
具体实施例：
[0276]
本发明采用改进的ieee30节点系统构建cehs模型，其结构如图6所示，其中节点1，2，8，13为火电机组、5为chp机组，节点11处接入150mw的风电。弃风惩罚成本850元/mwh。
[0277]
采用chp机组静态建模与动态建模两种方式参与cehs优化调度计算，并进行算例的对比分析。
[0278]
算例1：不考虑chp机组快速爬坡能力使用式(1)的静态模型，机组电功率爬坡取2％pe/min，可行域采用代数方程建模；
[0279]
算例2：chp动态模型，考虑chp机组优化控制的电热耦合特性，其采用微分代数方程建模如式(20)所示。
[0280]
在上述计算条件下，应用本发明方法进行调度优化，其结果如下：
[0281]
为对比chp机组快速爬坡能力特性对于调度方案影响，采用两种不同算例进行调度优化，各机组出力情况如图7所示。与算例1相比算例2的调度方案中，火电出力各时段基本持平，chp机组出力在多个时段明显减少，而风电出力较多时段显著增加。总体上看，算例2优化调度使chp机组电出力下降4.3％、风电出力增加9.79％。这是由于在多个时段chp机组的快速向下爬坡能力使得自身电出力降的更低，从而为消纳风电腾出更多空间。
[0282]
表1从cehs运行成本情况对2种算例的方案进行对比，可以看出，算例2调度成本减
少19681元，主要是由于风电消纳增加，弃风成本大幅降低15538元，chp、火电出力成本也小幅减少，从而降低了总成本。仿真结果表明，采用考虑chp机组快速爬坡能力动态建模的调度方案可多增加风电消纳、减少火电与chp电消耗。
[0283]
表1不同算例下调度方案成本
[0284][0285]
图8直观地展示了不同算例下24个时段chp机组运行点分布情况，可看出算例1运行点分布较为集中，离可行域上下边界较远；而算例2有更多的运行点靠近下边界，且分布范围也更大。例如：算例1在5时段chp机组运行点为(259.03,181.58)，而算例2为(259.03,155.42)。相比于算例1，算例2有10个运行点更接近下边界、有5个点离上边界更近。这是因为chp机组快速动态调节作用使机组爬坡能力更强，机组电功率优化空间更大，从而验证机组运行灵活性更好。
[0286]
算例3：考虑chp机组快速爬坡能力采用静态模型式(1)所示时，爬坡率取4.5％pe/min，可行域采用代数方程建模。
[0287]
经计算与算例2相比，算例3的调度成本、各机组出力情况基本相同(算例2，24.2422万；算例3，24.1411万)。为验证考虑快速爬坡能力时动态、静态模型调度方案下，chp机组模型对电出力跟踪准确性、系统安全性、供热稳定性的影响，将算例3各时段的chp机组模型对电出力跟踪准确性、系统安全性、供热稳定性的影响，将算例3各时段的形成阶跃扰动输入chp机组动态模型进行检验，其结果与算例2的动态模型进行比较。
[0288]
算例2chp机组电出力调度值与机组优化控制过程值，如图9所示。t11时段机组电出力由10时刻230.35mw调节到11时刻217.16mw，调节时间约100s，时段内平均缺额0.076mw/s；t12时段由11时刻217.16mw变负荷到12时刻157.16mw，调节时间约150s，平均缺额0.27mw/s。机组优化控制周期(一个调度时段)内过程值与调度值之间存在一定的爬坡缺额，这是由于优化控制过程阀门动作响应需调节时间。但控制调节速度较快，平均缺额量较小，可通过微调其他机组出力补偿，对调度方案可行性影响不大。算例3建模方式虽可达到与动态模型类似的爬坡能力效果，但无法考虑上述机组调节过程的动态特性，将导致实际机组调节失配，降低调度方案的可行性。
[0289]
在调度方案执行过程中，chp机组快速调节时主蒸汽压力变化关系到机组的安全稳定运行。如图10所示，主蒸汽压力设定值为15.00mpa时，向上最大峰值波动15.74mpa、向下最大峰值波动14.39mpa，未超设定值的波动范围，且稳态均达到设定值。11、17、18、24时段主蒸汽压力设定值为15.50mpa，向上、向下最大峰值波动均未超设定值的波动范围，稳态均达到设定值。仿真数据证明，chp机组快速动态调节过程产生的压力波动对机组安全并无影响。算例3建模方式无法考虑主蒸汽压力的动态调节过程约束，变负荷调节过大将导致机组调节过程压力超过安全约束的上限值。如图10中的2、4、10、16、23时段，将给机组安全稳定运行带来影响。
[0290]
chp机组快速爬坡能力动态过程中，如图11所示，各时段稳态值均达到设定值，调
节过程的峰值波动未超出设定值的偏差范围。静态模型流量调度值由满足热负荷需求的机组热出力求得，图中动态过程的稳态值与之重合，向上、向下波动计算平均流量等效值如红色曲线所示。与算例2相比，算例3在多个时段调节过程的峰值波动更大，时段内平均流量等效值偏离需求稳态值更远。
[0291]
各调度时段内chp机组热出力波动引起的热量、温度变化，如表2所示，其中m
h,equ
为chp机组等效平均抽汽蒸汽流量(t/h)，δq为热量变化情况(kj)，δτ为300s温度变化情况(
×
10-3
℃)。表2数据显示了算例2各时段热量增减引起的温度变化情况。总体来看，调度周期(2h)内chp机组平均多产生热量3845.21kj，经式(6)、式(9)计算该热量导致热网循环水温度平均升高1.07℃，引起供热区域温度升高3.62
×
10-3
℃。数据表明，chp机组动态调节过程产生的热量波动对供热温度需求并无影响。
[0292]
表2抽汽蒸汽流量波动对机组热出力与温度的影响
[0293]