CN104122797B - 一种新型火电机组负荷多变量预测控制方法 - Google Patents

Abstract

本发明公开了一种新型火电机组负荷多变量预测控制方法，该方法针对锅炉对象大惯性，过程响应缓慢等特性，提出一种新型性能指标，并根据此指标设计机组负荷控制系统的多变量预测控制方法。该方法的性能指标基于机组输出功率、主蒸汽压力的单步预测偏差和偏差变化率的关系进行设计。在过程控制中当外界负荷需求改变时，通过调整汽轮机调门开度，改变汽机进汽量，通过改变主蒸汽压力来释放锅炉蓄热，使得机组快速响应外界负荷需求；同时锅炉燃烧率快速调节，满足外界负荷需求的同时，防止主蒸汽压力与设定值之间偏差超限，并能够跟踪设定值，维持机组稳定运行，提高了负荷控制系统的稳定性和机组的运行效率。

Description

一种新型火电机组负荷多变量预测控制方法

技术领域

本发明属于热工自动控制技术领域，具体涉及一种新型火电机组负荷多变量预测控制方法。

背景技术

目前，火电机组普遍采用分散控制系统DCS，其中负荷控制系统采用负荷指令前馈和比例积分微分反馈的控制方法，但是该控制方法的控制品质不够理想，主要存在机组负荷升、降速率慢且主蒸汽压力波动大，使机组的其它重要参数如主蒸汽温度、锅炉汽包水位，炉膛负压及氧量等有较大的波动，使的机组无法快速满足外界负荷需求，负荷调节过程中机组无法稳定经济运行。因此研究设计先进的快速负荷响应控制算法具有特别重要的意义。

普通多变量预测控制方法应用于负荷协调控制系统时，存在计算量大，工程无法实现在线实时控制；同时，不能够实现锅炉侧与汽机侧在控制过程中的动态协调，即外界负荷需求改变时，通过短时间内改变主蒸汽压力，释放锅炉蓄热来满足负荷需求，同时在锅炉侧控制控制作用跟进时，能够保证主蒸汽压力跟踪压力设定值变化。因此研究新型火电机组负荷多变量预测控制方法，简化算法计算量，提高控制系统实时性，实现锅炉侧和汽机侧的动态协调，加快机组负荷响应速度，保证机组运行安全稳定就显得尤其重要。

发明内容

发明目的：针对上述问题，提出一种新型火电机组负荷多变量预测控制方法，改善负荷控制系统控制品质，使得机组负荷能够快速响应外界需求，同时保证机组运行的安全稳定性，且该方法计算量较小，能够实现工程在线实时控制，同时算法中存在可调参数，能够有效调节控制系统的控制性能。

技术方案：本发明根据火电机组负荷控制系统的动态特性，设计基于新型性能指标的多变量预测控制方法，并将该方法应用于负荷系统的协调控制中。在新的性能指标中仅考虑机组负荷的单步预测偏差及偏差的变化率，符合控制原理，能够减少算法计算量，实现快速控制；在机组负荷控制过程中，能够实现锅炉侧和蒸汽侧的动态协调，加快负荷响应速度，实现机组稳定经济运行。同时，该方法中可通过调节控制参数，对控制系统的控制性能进行有效的调节。

本发明的实现步骤如下：

步骤1：分别以汽机调门开度u₁，锅炉燃烧率u₂为阶跃量，获取机组负荷y₁和主蒸汽压力y₂的阶跃响应值；通过最小二乘系统辨识方法对上述响应数据进行拟合，得到下列传递函数：以汽机调门开度u₁为输入，负荷y₁为输出的汽机调门—负荷传递函数G₁₁；以锅炉燃烧率u₂为输入，负荷y₁为输出的燃烧率—负荷传递函数G₁₂；以汽机调门开度u₁为输入，主蒸汽压力y₂为输出的汽机调门—主蒸汽压力传递函数G₂₁；以锅炉燃烧率u₂为输入，主蒸汽压力y₂为输出的燃烧率—主蒸汽压力传递函数G₂₂；

步骤2：选择采样时间T，T首先满足香农定理，然后使T₉₅/T取值在10到25之间，T₉₅为机组负荷升至设定值95％的过程时间；对上述传递函数模型进行离散化，得到负荷控制系统的CARIMA模型： $A_1\left(z^{-1}\right)y_1\left(k\right)=B_{\mathrm{1,1}}\left(z^{-1}\right)u_1(k-1)+B_{\mathrm{1,2}}\left(z^{-1}\right)u_2(k-1)+\frac{1}{\mathrm{\Δ}}{\mathrm{\ϵ}}_1\left(k\right);$ $A_2\left(z^{-1}\right)y_2\left(k\right)=B_{2,1}\left(z^{-1}\right)u_1(k-1)+B_{2,2}\left(z^{-1}\right)u_2(k-1)+\frac{1}{\mathrm{\Δ}}{\mathrm{\ϵ}}_2\left(k\right);$ 式中z^-1为后移算子，Δ＝1-z^-1； $A_i\left(z^{-1}\right)=1+\underset{l=1}{\overset{{\mathrm{na}}_i}{\mathrm{\Σ}}}{a}_{i,l}{z}^{-l},B_{i,j}\left(z^{-1}\right)=\underset{l=0}{\overset{{\mathrm{nb}}_{i,j}}{\mathrm{\Σ}}}{b}_{i,j}^l{z}^{-l},$ i＝1,2,j＝1,2，na_i和nb_i,j分别为多项式A_i(z^-1)和B_i,j(z^-1)阶次，a_i,l和分别为多项式A_i(z^-1)和B_i,j(z^-1)中z^-l项的系数；y₁(k)，y₂(k)分别为机组负荷和主蒸汽压力在k时刻的输出值；u₁(k-1)，u₂(k-1)分别为汽机调门开度和锅炉燃烧率在k-1时刻的控制量，ε₁(k)和ε₂(k)为均值为0的白噪声；

步骤3：构造丢番图方程：1＝E_i,j(z^-1)A_i(z^-1)Δ+z^-jF_i,j(z^-1)，求得多项式E_i,j(z^-1)和多项式F_i,j(z^-1)，其中 $E_{i,j}\left(z^{-1}\right)=\underset{l=0}{\overset{j-1}{\mathrm{\Σ}}}{e}_{j,l}^i{z}^{-l},F_{i,j}\left(z^{-1}\right)=\underset{l=0}{\overset{{\mathrm{na}}_i}{\mathrm{\Σ}}}{f}_{j,l}^i{z}^{-l},$ 和分别为多项式E_i,j(z^-1)和F_i,j(z^-1)中z^-l项的系数，i＝1,2分别代表负荷对象和主蒸汽压力对象，j＝N_i-1,N_i，N_i为对应预测时域长度；

步骤4：定义算式： $G_j^{i,1}\left(z^{-1}\right)=E_{i,j}\left(z^{-1}\right)B_{i,1}\left(z^{-1}\right),G_j^{i,2}\left(z^{-1}\right)=E_{i,j}\left(z^{-1}\right)B_{i,2}\left(z^{-1}\right),$ 求得： $G_j^{i,1}\left(z^{-1}\right)=\underset{l=0}{\overset{j+{\mathrm{nb}}_{i,1}-1}{\mathrm{\Σ}}}{g}_{j,l}^{i,1}{z}^{-l},G_j^{i,2}\left(z^{-1}\right)=\underset{l=0}{\overset{j+{\mathrm{nb}}_{i,2}-1}{\mathrm{\Σ}}}{g}_{j,l}^{i,2}{z}^{-l},$ 式中 分别为多项式和中z^-l项的系数，i＝1,2，j＝N_i-1,N_i；

步骤5：构造一维向量： $G_{i,N_i-1}=[g_{N_i-1,N_i-2}^{i,1},...,g_{N_i-1{,N}_i-1-{\mathrm{Nu}}_1}^{i,1},g_{N_i-1,N_i-2}^{i,2},...,g_{N_i-1,N_i-1-{\mathrm{Nu}}_2}^{i,2}],$ $G_{i,N_i}=[g_{N_i,N_i-1}^{i,1},...g_{N_i,N_i-{\mathrm{Nu}}_1}^{i,1},g_{N_i,N_i-1}^{i,2},...,g_{N_i,N_i-{\mathrm{Nu}}_2}^{i,2}],$ 式中i＝1,2，Nu₁为汽机调门开度的控制时域长度，Nu₂为锅炉燃烧率的控制时域长度，Nu₁、Nu₂一般取1或2；

步骤6：该控制系统在线实时记录控制系统从开始运行至采样k时刻的汽机调门开度和锅炉燃烧率的控制信号{u₁(k-1),u₁(k-2),…u₁(1)}，{u₂(k-1),u₂(k-2),…u₂(1)}，以及机组负荷和主蒸汽压力的测量数据{y₁(k),y₁(k-1),…y₁(1)}，{y₂(k),y₂(k-1),…y₂(1)}，同时记录相应控制量的增量{Δu₁(k-1),Δu₁(k-2),…Δu₁(1)}，{Δu₂(k-1),Δu₂(k-2),…Δu₂(1)}，并接收预测控制系统的机组负荷设定值y_r,1(k)和主蒸汽压力设定值y_r,2(k)；定义中间控制变量：

$f_{1,j}=\underset{i=0}{\overset{{\mathrm{na}}_1}{\mathrm{\Σ}}}{f}_{f,i}^{1}y(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{1,1}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{\mathrm{1,1}}\mathrm{\Δ}{u}_1(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{1,2}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{\mathrm{1,2}}\mathrm{\Δ}{u}_2(k-i)(j=N_1-1,N_1)$ $f_{2,j}=\underset{i=0}{\overset{{\mathrm{na}}_2}{\mathrm{\Σ}}}{f}_{f,i}^{2}y(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{2,1}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{2,1}\mathrm{\Δ}{u}_1(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{2,2}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{2,2}\mathrm{\Δ}{u}_2(k-i)(j=N_2-1,N_2);$

步骤7：按下式对机组负荷在未来第k+N₁-1和k+N₁时刻的输出和主蒸汽压力在未来第k+N₂-1和k+N₂时刻的输出值进行预测：

y₁(k+j|k)＝G_1,j·ΔU+f_1,j(j＝N₁-1,N₁)

y₂(k+j|k)＝G_2,j·ΔU+f_2,j(j＝N₂-1,N₂)

式中，y₁(k+j|k)为未来k+j时刻机组负荷预测值，y₂(k+j|k)为未来k+j时刻主蒸汽压力的预测值，ΔU＝[Δu₁(k),…,Δu₁(k+Nu₁-1),Δu₂(k),…,Δu₂(k+Nu₂-1)]，式中Δu₁(k+i),i＝0…Nu₁-1，Δu₂(k+i),i＝0…Nu₂-1分别表示在k时刻对未来k+i时刻汽机调门开度和锅炉燃烧率的输出控制量增量的估计值；

步骤8：按下列各式求取机组负荷在k+N₁时刻的预测偏差e₁(k+N₁)和预测偏差变化率ce₁(k+N₁)，主蒸汽压力在k+N₂时刻的预测偏差e₂(k+N₂)和预测偏差变化率ce₂(k+N₂)：

e_i(k+N_i)＝y_i(k+N_i|k)-y_r,i(k+N_i)

e_i(k+N_i-1)＝y_i(k+N_i-1|k)-y_r,i(k+N_i-1)

ce_i(k+N_i)＝e_i(k+N_i)-e_i(k+N_i-1)

式中，i＝1,2，分别代表机组负荷和主蒸汽压力；y_r,i(k+N_i-1)和y_r,i(k+N_i)分别为设定值y_r,i(k)在k+N_i-1和k+N_i时刻的预测值，取y_r,i(k+N_i-1)＝y_r,i(k+N_i)，按下式计算：

y_r,i(k+N_i)＝y_r,i(k)+α·N_i·δ

式中，δ为设定值的变化速率；α为可调参数，在[0,1]间取值，通过调整α的值大小来调整控制系统跟踪设定值变化的快慢，α取大则控制系统跟踪设定值变快，反之则变慢，在设定值不变或设定值阶跃变化时，α取为0；

步骤9：根据新型性能指标：

J＝[e₁(k+N₁)+λ₁·ce₁(k+N₁)]²+α·[e₂(k+N₂)+λ₂·ce₂(k+N₂)]²

来计算汽机调门指令u₁(k)和锅炉燃烧率指令u₂(k)。式中，λ₁和λ₂分别为机组负荷和主汽压力的预测偏差变化率，一般取0～5；α为权值系数，根据机组负荷和主汽压力间相对变化关系进行选择；

将步骤7中的未来预测值带入上述性能指标中，并通过极值必要条件 $\frac{\mathrm{dJ}}{\mathrm{d\ΔU}}=0$ 求得：

$\mathrm{\ΔU}=-{(M_1^T\·M_1+Q^T\·M_2^T\·M_2\·Q)}^{-1}\·{[h_1\·M_1+\mathrm{\α}\·h_2\·M_2]}^T$

式中， $M_1=G_{1,N_1}+{\mathrm{\λ}}_1\·(G_{1,N_1}-G_{1,N_1-1}),M_2=G_{2,N_2}+{\mathrm{\λ}}_2\·(G_{2,N_2}-G_{2,N_2-1}),$

$Q=\sqrt{\mathrm{\α}}\·I_{({\mathrm{Nu}}_1+{\mathrm{Nu}}_2)},$

$h_1=f_{1,N_1}-y_{r,1}(k+N_1)+{\mathrm{\λ}}_1\·(f_{1,N_1}-f_{1,N_1-1}),$

$h_2=f_{2,N_2}-y_{r,2}(k+N_2)+{\mathrm{\λ}}_2\·(f_{2,N_2}-f_{2,N_2-1});$

步骤10：计算汽机调门开度控制量u₁(k)＝u₁(k-1)+Δu₁(k)，和锅炉燃烧率控制量u₂(k)＝u₂(k-1)+Δu₂(k)；

步骤11：根据执行机构限制及实际运行情况，设定控制量约束，汽机调门开度约束上下限为[u_1,min，u_1,max]，锅炉燃烧率约束上下限为[u_2,min，u_2,max]，如果u_i(k)＞u_i,max，则令u_i(k)＝u_i,max,Δu_i(k)＝u_i,max-u_i(k-1)；如果u_i(k)＜u_i,min，则令u_i(k)＝u_i,min,Δu_i(k)＝u_i,min-u_i(k-1)，i＝1,2；

步骤12：控制系统输出u₁(k)和u₂(k)到相应执行机构，在之后的每个采样周期内重复执行步骤6到步骤12。

有益效果：基于新型性能指标，设计火电机组负荷控制系统的多变量预测控制方法，该方法能够实现过程控制中锅炉侧和蒸汽侧的动态协调，加快负荷响应速度，实现机组稳定经济运行。同时，该算法基于单步预测偏差和预测偏差变化率设计，算法计算简单，能够满足控制系统实时性的要求。

附图说明

图1是新型火电机组负荷多变量预测控制方法的控制结构图；

图2-a至2-d是新型火电机组负荷多变量预测控制系统在机组滑压运行方式下的负荷控制情况，其中图2-a、图2-b分别为机组负荷和主蒸汽压力的控制过程响应曲线；图2-c、图2-d分别为锅炉燃烧率和汽机调门开度的控制量曲线。

具体实施方式

下面结合附图和具体实施例，对本发明做更进一步的解释。

步骤1：分别以汽机调门开度u₁，锅炉燃烧率u₂为阶跃量，获取机组负荷y₁和主蒸汽压力y₂的阶跃响应值；通过最小二乘系统辨识方法对上述响应数据进行拟合，得到下列传递函数：以汽机调门开度u₁为输入，负荷y₁为输出的汽机调门—负荷传递函数G₁₁；以锅炉燃烧率u₂为输入，负荷y₁为输出的燃烧率—负荷传递函数G₁₂；以汽机调门开度u₁为输入，主蒸汽压力y₂为输出的汽机调门—主蒸汽压力传递函数G₂₁；以锅炉燃烧率u₂为输入，主蒸汽压力y₂为输出的燃烧率—主蒸汽压力传递函数G₂₂；

本实施例中上述辨识得到的传递函数模型分别为：

$G_{11}\left(s\right)=\frac{760s}{(148s+1)(18s+1)},G_{12}\left(s\right)=\frac{7.901(1-19s)}{(165s+1)(148s+1)(18s+1)}{e}^{-33},$

$G_{21}\left(s\right)=\frac{-0.132}{1+148s},G_{22}\left(s\right)=\frac{0.203(1-19s)}{(165s+1)(148s+1)}{e}^{-33s};$

步骤2：选择采样时间T＝11s，对上述传递函数模型进行离散化，得到负荷控制系统的CARIMA模型： $A_1\left(z^{-1}\right)y_1\left(k\right)=B_{\mathrm{1,1}}\left(z^{-1}\right)u_1(k-1)+B_{\mathrm{1,2}}\left(z^{-1}\right)u_2(k-1)+\frac{1}{\mathrm{\Δ}}{\mathrm{\ϵ}}_1\left(k\right);$ $A_2\left(z^{-1}\right)y_2\left(k\right)=B_{2,1}\left(z^{-1}\right)u_1(k-1)+B_{2,2}\left(z^{-1}\right)u_2(k-1)+\frac{1}{\mathrm{\Δ}}{\mathrm{\ϵ}}_2\left(k\right);$ 式中z^-1为后移算子，Δ＝1-z^-1；y₁(k)，y₂(k)分别为机组负荷和主蒸汽压力在k时刻的输出值；u₁(k-1)，u₂(k-1)分别为汽机调门开度和锅炉燃烧率在k-1时刻的控制量，ε₁(k)和ε₂(k)为均值为0的白噪声；

A₁(z^-1)＝1-3.8777z^-1+5.9244z^-2-4.45z^-3+1.6408z^-4-0.2375z^-5

A₂(z^-1)＝1-2.7922z^-1+2.5989z^-2-0.8063z^-3

B_1,1(z^-1)＝2.2544-7.6799z^-1+9.6641z^-2-5.3012z^-3+1.0627z^-4

B_1,2(z^-1)＝z^-3(-0.01288+0.03364z^-1-0.01318z^-2-0.01453z^-3+0.00751z^-4)

B_2,1(z^-1)＝-0.009455+0.01762z^-1-0.008212z^-2

B_2,2(z^-1)＝z^-3(-0.001139+0.003134z^-1-0.001928z^-2)

步骤3：构造丢番图方程：1＝E_1,j(z^-1)A₁(z^-1)Δ+z^-jF_1,j(z^-1)，求得机组负荷离散模型的分离第一多项式E_1,j(z^-1)和分离第二多项式F_1,j(z^-1)，其中 为多项式F_1,j(z^-1)中z^-l项的系数，j＝N₁-1,N₁，N₁为机组输出负荷的预测时域长度，需大于负荷对象的纯延迟时间；

本实例中，机组输出负荷的预测时域长度N₁＝14，其中和为：

$F_{1,N_1-1}\left(z^{-1}\right)=954.30+{3523.71z}^{-1}+{158.50z}^{-2}+{3734.38z}^{-3}+{1334.40z}^{-4}-{188.10z}^{-5}$

$F_{1,N_1}\left(z^{-1}\right)=1131.11-{4195.68z}^{-1}+{6165.77z}^{-2}-{4477.90z}^{-3}+{1604.35z}^{-4}-{226.66z}^{-5}$ $;$

步骤4：本实施例中，选择汽机调门开度和锅炉燃烧率的控制时域长度Nu₁、Nu₂都取为1。定义算式： $G_j^{1,1}\left(z^{-1}\right)=E_{1,j}\left(z^{-1}\right)B_{1,1}\left(z^{-1}\right),G_j^{1,2}\left(z^{-1}\right)=E_{1,j}\left(z^{-1}\right)B_{1,2}\left(z^{-1}\right),$ 求得： $G_j^{1,1}\left(z^{-1}\right)=\underset{l=0}{\overset{j+{\mathrm{nb}}_{1,1}-1}{\mathrm{\Σ}}}{g}_{j,l}^{1,1}{z}^{-l},G_j^{1,2}\left(z^{-1}\right)=\underset{l=0}{\overset{j+{\mathrm{nb}}_{1,2}-1}{\mathrm{\Σ}}}{g}_{j,l}^{1,2}{z}^{-l},$ 式中 分别为多项式和中z^-l项的系数，j＝N₁-1,N₁；构造一维向量： $G_{1,N_1-1}=[g_{N_1-1,N_1-2}^{\mathrm{1,1}},g_{N_1-1,N_1-2}^{\mathrm{1,2}}]=\left[\mathrm{0.641,2.223}\right],G_{1,N_1}=[g_{N_1,N_1-1}^{\mathrm{1,1}},g_{N_1,N_1-1}^{\mathrm{1,2}}]=\left[\mathrm{0.809,2.064}\right];$

步骤5：本实施例中主蒸汽压力的预测时域长度N₂选取为8，针对主蒸汽压力对象，计算过程与步骤3-4相同，求得：

$F_{2,N_2-1}\left(z^{-1}\right)=82.97-{206.32z}^{-1}+{173.69z}^{-2}-{49.33z}^{-3}$

$F_{2,N_2}\left(z^{-1}\right)=108.32-{273.61z}^{-1}+{233.19z}^{-2}-{66.90z}^{-3}$

$G_{2,N_2-1}=[g_{N_2-1,N_2-2}^{\mathrm{2,1}},g_{N_2-1,N_2-2}^{\mathrm{2,2}}]=[0.001441,-0.05354]$

$G_{2,N_2}=[g_{N_2,N_2-1}^{\mathrm{2,1}},g_{N_2,N_2-1}^{\mathrm{2,2}}]=[0.003869,-0.05916];$

步骤6：负荷控制系统在线实时记录控制系统从开始运行至采样k时刻的汽机调门开度和锅炉燃烧率的控制信号{u₁(k-1),u₁(k-2),…u₁(1)}，{u₂(k-1),u₂(k-2),…u₂(1)}，以及过机组负荷和主蒸汽压力的测量数据{y₁(k),y₁(k-1),…y₁(1)}，{y₂(k),y₂(k-1),…y₂(1)}，同时记录相应控制量的增量{Δu₁(k-1),Δu₁(k-2),…Δu₁(1)}，{Δu₂(k-1),Δu₂(k-2),…Δu₂(1)}，并接收预测控制系统的机组负荷设定值y_r,1(k)和主蒸汽压力设定值y_r,2(k)；定义中间控制变量：

$f_{1,j}=\underset{i=0}{\overset{{\mathrm{na}}_1}{\mathrm{\Σ}}}{f}_{f,i}^{1}y(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{1,1}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{\mathrm{1,1}}\mathrm{\Δ}{u}_1(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{1,2}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{\mathrm{1,2}}\mathrm{\Δ}{u}_2(k-i)(j=N_1-1,N_1);$

$f_{2,j}=\underset{i=0}{\overset{{\mathrm{na}}_2}{\mathrm{\Σ}}}{f}_{f,i}^{2}y(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{2,1}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{2,1}\mathrm{\Δ}{u}_1(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{2,2}}{\mathrm{\Σ}}}{g}_{j,j+1-i}^{2,2}\mathrm{\Δ}{u}_2(k-i)(j=N_2-1,N_2);$

步骤7：按下式对机组负荷在未来第k+N₁-1和k+N₁时刻的输出和主蒸汽压力在未来第k+N₂-1和k+N₂时刻的输出值进行预测：

y₁(k+j|k)＝G_1,j·ΔU+f_1,j(j＝N₁-1,N₁)

y₂(k+j|k)＝G_2,j·ΔU+f_2,j(j＝N₂-1,N₂)

式中，y₁(k+j|k)为未来k+j时刻机组负荷预测值，y₂(k+j|k)为未来k+j时刻主蒸汽压力的预测值，ΔU＝[Δu₁(k),…,Δu₁(k+Nu₁-1),Δu₂(k),…,Δu₂(k+Nu₂-1)]，式中Δu₁(k+i),i＝0…Nu₁-1，Δu₂(k+i),i＝0…Nu₂-1分别表示在k时刻对未来k+i时刻汽机调门开度和锅炉燃烧率的输出控制量增量的估计值；

步骤8：按下列各式求取机组输出功率在k+N₁时刻的预测偏差e₁(k+N₁)和预测偏差变化率ce₁(k+N₁)，以及主汽压力在k+N₂时刻的预测偏差e₂(k+N₂)和预测偏差变化率ce₂(k+N₂)：

e_i(k+N_i)＝y_i(k+N_i|k)-y_r,i(k+N_i)

e_i(k+N_i-1)＝y_i(k+N_i-1|k)-y_r,i(k+N_i-1)

ce_i(k+N_i)＝e_i(k+N_i)-e_i(k+N_i-1)

式中，i＝1,2，分别代表机组输出负荷和主汽压力；y_r,i(k+N_i-1)和y_r,i(k+N_i)分别为设定值y_r,i(k)在k+N_i-1和k+N_i时刻的预测值，取y_r,i(k+N_i-1)＝y_r,i(k+N_i)，按下式计算：

y_r,i(k+N_i)＝y_r,i(k)+α·N_i·δ

式中，δ为设定值的变化速率；α为可调参数，在[0,1]间取值，通过调整α的值大小来调整控制系统跟踪设定值变化的快慢，α取大则控制系统跟踪设定值变快，反之则变慢，在设定值不变或设定值阶跃变化时，α取为0；

本实施例中，负荷设定值变化速率δ₁取-2.2，可调参数α₁选取为0.33，主汽压力设定值变化速率δ₂取为-0.05742，可调参数α₂选取为0.1；

步骤9：根据新型性能指标：

J＝[e₁(k+N₁)+λ₁·ce₁(k+N₁)]²+α·[e₂(k+N₂)+λ₂·ce₂(k+N₂)]²

来计算汽机调门指令u₁(k)和锅炉燃烧率指令u₂(k)，式中λ₁和λ₂分别为机组负荷和主汽压力的预测偏差变化率，本实施例中λ₁取为0.1，λ₂取为2；α为权值系数，本实施例中取为0.5；

将步骤7中的未来预测值带入上述性能指标中，并通过极值必要条件 $\frac{\mathrm{dJ}}{\mathrm{d\ΔU}}=0$ 求得：

$\mathrm{\ΔU}=-{(M_1^T\·M_1+Q^T\·M_2^T\·M_2\·Q)}^{-1}\·{[h_1\·M_1+\mathrm{\α}\·h_2\·M_2]}^T$

式中 $M_1=G_{1,N_1}+{\mathrm{\λ}}_1\·(G_{1,N_1}-G_{1,N_1-1}),M_2=G_{2,N_2}+{\mathrm{\λ}}_2\·(G_{2,N_2}-G_{2,N_2-1}),$

$Q=\sqrt{\mathrm{\α}}\·I_{({\mathrm{Nu}}_1+{\mathrm{Nu}}_2)},$

$h_1=f_{1,N_1}-y_{r,1}(k+N_1)+{\mathrm{\λ}}_1\·(f_{1,N_1}-f_{1,N_1-1}),$

$h_2=f_{2,N_2}-y_{r,2}(k+N_2)+{\mathrm{\λ}}_2\·(f_{2,N_2}-f_{2,N_2-1});$

所以ΔU_M(k)＝[Δu₁(k)Δu₂(k)]^T＝D·ΔU，式中D＝diag(θ₁,θ₂)， ${\mathrm{\θ}}_i={\left[\begin{array}{cccc}1& 0& ...& 0\end{array}\right]}_{1\×{\mathrm{Nu}}_i},$ i＝1,2；

步骤10：计算汽机调门开度控制量u₁(k)＝u₁(k-1)+Δu₁(k)，和锅炉燃烧率控制量u₂(k)＝u₂(k-1)+Δu₂(k)；

步骤11：本实施例中，汽机调门开度约束上下限为[0，100]，锅炉燃烧率约束上下限为[0，100]。如果u_i(k)＞u_i,max，则令u_i(k)＝u_i,max,Δu_i(k)＝u_i,max-u_i(k-1)；如果u_i(k)＜u_i,min，则令u_i(k)＝u_i,min,Δu_i(k)＝u_i,min-u_i(k-1)，i＝1,2；

步骤12：控制系统输出u₁(k)和u₂(k)到相应执行机构，在之后的每个采样周期内重复执行步骤6到步骤12。

上述实施例的仿真结果如图2-a至2-d所示。图2-a和图2-b表明该负荷控制方法在滑压运行方式下可以很好地跟踪负荷设定值指令和压力设定值指令。图2-c和图2-d表明在整个调节过程中，汽机调门开度、锅炉燃烧率等控制量的响应速度很快，震荡幅度较小，震荡时间短。且在整个控制过程中，机组负荷和主汽压力的波动幅度很小，波动时间短，能够快速跟踪设定值的变化，满足外界负荷变动和机组稳定运行的需求。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.一种新型火电机组负荷多变量预测控制方法，其步骤如下：

步骤1：分别以汽机调门开度u₁，锅炉燃烧率u₂为阶跃量，获取机组负荷y₁和主蒸汽压力y₂的阶跃响应值；通过最小二乘系统辨识方法对上述响应值进行拟合，得到下列传递函数：以汽机调门开度u₁为输入，负荷y₁为输出的汽机调门—负荷传递函数G₁₁；以锅炉燃烧率u₂为输入，负荷y₁为输出的燃烧率—负荷传递函数G₁₂；以汽机调门开度u₁为输入，主蒸汽压力y₂为输出的汽机调门—主蒸汽压力传递函数G₂₁；以锅炉燃烧率u₂为输入，主蒸汽压力y₂为输出的燃烧率—主蒸汽压力传递函数G₂₂；

步骤2：选择采样时间T，T首先满足香农定理，然后使T₉₅/T取值在10到25之间，T₉₅为机组负荷升至设定值95％的过程时间；对上述传递函数模型进行离散化，得到负荷预测控制系统的CARIMA模型：

$A_1\left(z^{-1}\right)y_1\left(k\right)=B_{1,1}\left(z^{-1}\right)u_1(k-1)+B_{1,2}\left(z^{-1}\right)u_2(k-1)+\frac{1}{\mathrm{\Δ}}{\mathrm{\ϵ}}_1\left(k\right);$

$A_2\left(z^{-1}\right)y_2\left(k\right)=B_{2,1}\left(z^{-1}\right)u_1(k-1)+B_{2,2}\left(z^{-1}\right)u_2(k-1)+\frac{1}{\mathrm{\Δ}}{\mathrm{\ϵ}}_2\left(k\right);$

式中z^-1为后移算子，Δ＝1-z^-1；i＝1,2,j＝1,2，na_i和nb_i,j分别为多项式A_i(z^-1)和B_i,j(z^-1)阶次，a_i,l和分别为多项式A_i(z^-1)和B_i,j(z^-1)中z^-l项的系数；y₁(k)，y₂(k)分别为机组负荷和主蒸汽压力在k时刻的输出值；u₁(k-1)，u₂(k-1)分别为汽机调门开度和锅炉燃烧率在k-1时刻的控制量，ε₁(k)和ε₂(k)为均值为0的白噪声；

步骤3：构造丢番图方程：1＝E_i,j(z^-1)A_i(z^-1)Δ+z^-jF_i,j(z^-1)，求得多项式E_i,j(z^-1)和多项式F_i,j(z^-1)，其中 和分别为多项式E_i,j(z^-1)和F_i,j(z^-1)中z^-l项的系数，i＝1,2分别代表负荷对象和主蒸汽压力对象，j＝N_i-1,N_i，N_i为对应预测时域长度；

步骤4：定义算式：求得：式中分别为多项式和中z^-l项的系数，i＝1,2，j＝N_i-1,N_i；

步骤5：构造一维向量： 式中i＝1,2，Nu₁为汽机调门开度的控制时域长度，Nu₂为锅炉燃烧率的控制时域长度，Nu₁、Nu₂取1或2；

步骤6：该预测控制系统在线实时记录预测控制系统从开始运行至采样k时刻的汽机调门开度和锅炉燃烧率的控制信号{u₁(k-1),u₁(k-2),…u₁(1)}，{u₂(k-1),u₂(k-2),…u₂(1)}，以及机组负荷和主蒸汽压力的测量数据{y₁(k),y₁(k-1),…y₁(1)}，{y₂(k),y₂(k-1),…y₂(1)}，同时记录相应控制量的增量{Δu₁(k-1),Δu₁(k-2),…Δu₁(1)}，{Δu₂(k-1),Δu₂(k-2),…Δu₂(1)}，并接收预测控制系统的机组负荷设定值y_r,1(k)和主蒸汽压力设定值y_r,2(k)；定义中间控制变量：

$f_{1,j}=\underset{i=0}{\overset{{\mathrm{na}}_1}{\Σ}}{f}_{j,i}^{1}y(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{1,1}}{\Σ}}{g}_{j,j+1-i}^{1,1}{\mathrm{\Δu}}_1(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{1,2}}{\Σ}}{g}_{j,j+1-i}^{1,2}{\mathrm{\Δu}}_2(k-i),j=N_1-1,N_1;$

$f_{2,j}=\underset{i=0}{\overset{{\mathrm{na}}_2}{\Σ}}{f}_{j,i}^{2}y(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{2,1}}{\Σ}}{g}_{j,j+1-i}^{2,1}{\mathrm{\Δu}}_1(k-i)+\underset{i=1}{\overset{{\mathrm{nb}}_{2,2}}{\Σ}}{g}_{j,j+1-i}^{2,2}{\mathrm{\Δu}}_2(k-i),j=N_2-1,N_2;$

步骤7：按下式对机组负荷在未来第k+N₁-1和k+N₁时刻的输出和主蒸汽压力在未来第k+N₂-1和k+N₂时刻的输出值进行预测：

y₁(k+j|k)＝G_1,j·ΔU+f_1,j,式中j＝N₁-1,N₁；

y₂(k+j|k)＝G_2,j·ΔU+f_2,j,式中j＝N₂-1,N₂；

式中，y₁(k+j|k)为未来k+j时刻机组负荷预测值，y₂(k+j|k)为未来k+j时刻主蒸汽压力的预测值，ΔU＝[Δu₁(k),…,Δu₁(k+Nu₁-1),Δu₂(k),…,Δu₂(k+Nu₂-1)]，式中Δu₁(k+i),i＝0…Nu₁-1，Δu₂(k+i),i＝0…Nu₂-1分别表示在k时刻对未来k+i时刻汽机调门开度和锅炉燃烧率的输出控制量增量的估计值；

步骤8：按下列各式求取机组负荷在k+N₁时刻的预测偏差e₁(k+N₁)和预测偏差变化率ce₁(k+N₁)，主蒸汽压力在k+N₂时刻的预测偏差e₂(k+N₂)和预测偏差变化率ce₂(k+N₂)：

e_i(k+N_i)＝y_i(k+N_i|k)-y_r,i(k+N_i)

e_i(k+N_i-1)＝y_i(k+N_i-1|k)-y_r,i(k+N_i-1)

ce_i(k+N_i)＝e_i(k+N_i)-e_i(k+N_i-1)

式中，i＝1,2，分别代表机组负荷和主蒸汽压力；y_r,i(k+N_i-1)和y_r,i(k+N_i)分别为设定值y_r,i(k)在k+N_i-1和k+N_i时刻的预测值，取y_r,i(k+N_i-1)＝y_r,i(k+N_i)，按下式计算：

y_r,i(k+N_i)＝y_r,i(k)+α·N_i·δ

式中，δ为设定值的变化速率；α为可调参数，在[0,1]间取值，通过调整α的值大小来调整预测控制系统跟踪设定值变化的快慢，α取大则预测控制系统跟踪设定值变快，反之则变慢，在设定值不变或设定值阶跃变化时，α取为0；

步骤9：根据新型性能指标：

J＝[e₁(k+N₁)+λ₁·ce₁(k+N₁)]²+α·[e₂(k+N₂)+λ₂·ce₂(k+N₂)]²

来计算汽机调门指令u₁(k)和锅炉燃烧率指令u₂(k)，式中，λ₁和λ₂分别为机组负荷和主汽压力的预测偏差变化率，取0～5；α为权值系数，根据机组负荷和主汽压力间相对变化关系进行选择；

将步骤7中的未来预测值带入上述性能指标中，并通过极值必要条件求得：

$\mathrm{\Δ}U=-{(M_1^T\·M_1+Q^T\·M_2^T\·M_2\·Q)}^{-1}\·{\[h_1\·M_1+\mathrm{\α}\·h_2\·M_2\]}^T$

式中，

步骤10：计算汽机调门开度控制量u₁(k)＝u₁(k-1)+Δu₁(k)，和锅炉燃烧率控制量u₂(k)＝u₂(k-1)+Δu₂(k)；

步骤11：根据执行机构限制及实际运行情况，设定控制量约束，汽机调门开度约束上下限为[u_1,min，u_1,max]，锅炉燃烧率约束上下限为[u_2,min，u_2,max]，如果u_i(k)＞u_i,max，则令u_i(k)＝u_i,max,Δu_i(k)＝u_i,max-u_i(k-1)；如果u_i(k)＜u_i,min，则令u_i(k)＝u_i,min,Δu_i(k)＝u_i,min-u_i(k-1)，i＝1,2；

步骤12：预测控制系统输出u₁(k)和u₂(k)到相应执行机构，在之后的每个采样周期内重复执行步骤6到步骤12。

2.根据权利要求1所述的一种新型火电机组负荷多变量预测控制方法，其特征在于：该预测控制系统预测时域长度N₁和N₂选择上只需大于对应对象的纯延迟时间，同时控制时域长度Nu₁和Nu₂。