CN104296131A - 一种双炉膛循环流化床机组的多变量协同控制方法 - Google Patents

Abstract

本发明公开了一种双炉膛循环流化床机组的多变量协同控制方法，它是构建循环流化床锅炉和汽轮机机组的一个多变量系统，采用基于跟踪误差变化率的多变量约束DMC算法。所采用的算法针对传统的DMC算法做出了改进，在目标函数中引入了跟踪误差变化率，提高了对设定值的跟踪效果，同时对输入输出附加了约束，使得求解的控制量在其物理范围内。由于控制系统同步考虑了左右两侧炉膛的床温和床压，在变工况调节时既能保证负荷跟踪能力，又能实现两个炉膛的床温床压平衡稳定和均衡燃烧。

Description

一种双炉膛循环流化床机组的多变量协同控制方法

技术领域

本发明属于循环流化床发电厂控制领域，涉及循环流化床锅炉负荷、床温和床压的控制方法。

背景技术

循环流化床(Circulating Fluidized Bed，CFB)的锅炉燃烧技术作为高效能、低污染、燃烧煤泥矸石等劣质煤、适应性广的洁净燃煤技术，在全世界越来越受到广泛重视，正在成为燃煤技术的主力军。国家“十二五”规划明确提出“加大节能力度；推进粉煤灰、煤矸石、冶金和化工废渣及尾矿等工业废物利用”，给大型循环流化床锅炉的技术研究和产业化发展带来了的机遇，目前300MW循环流化床锅炉已逐步成为国内重要的发电锅炉形式，600MW超临界循环流化床锅炉也已有试点运行。

CFB锅炉炉膛内含有大量的物料，在燃烧过程中大量的物料被烟气携带到炉膛上部，经过布置在炉膛出口的分离器，将物料与烟气分开，并经过非机械式回料阀将物料直接回送至炉内或经外置床受热面后送至炉内，多次循环燃烧。由于物料浓度高，具有很大的热容量和良好的物料混合，一般每公斤烟气可携带若干公斤的物料，这些循环物料带来了高传热系数，使锅炉热负荷调节范围广，对燃料的适应性强。循环流化床锅炉的循环流化方式增加了燃料在炉膛内的停留时间，因此具有较高的燃烧效率，在低负荷下能稳定运行，而无需增加辅助燃料。

CFB锅炉运行温度通常在800～900℃之间，这是一个理想的脱硫温度区间，采用炉内脱硫技术，向床内加入石灰石和脱硫剂，燃料及脱硫剂经多次循环，反复进行低温燃烧和脱硫反应，加之炉内湍流运动剧烈，Ca/S摩尔比约为2时，可以使脱硫效率达到90％左右，SO2的排放量大大降低。同时循环流化床采用低温分级送风燃烧，使燃烧始终在低过量空气下进行，从而大大降低了NOx的生成和排放。

循环流化床锅炉的控制点多，燃烧系统耦合性强，燃烧响应滞后性大是循环流化床控制的难点，也是循环流化床锅炉本体结构和燃烧煤种所决定的。循环流化床控制难的问题是影响循环流化床锅炉又好又快发展的一个重要因素，如何解决好循环流化床特别是大容量循环床锅炉燃烧控制的稳定性、鲁棒性问题是一个迫切需要解决的问题。循环流化床锅炉的主要控制问题需要解决变工况时的床温和床压控制的稳定性，负荷控制的升降速率保证，以及国内主流300MW循环流化床锅炉双炉膛的床压平衡和均衡燃烧问题。

目前针对循环流化床锅炉建模和控制领域的研究工作主要体现在两个层面：一是通过分析循环流化床锅炉运行机理，建立锅炉机理模型，或基于现场试验所测数据，对不同工况下各变量参数之间的关系进行分析，建立控制系统设计所需要的控制模型。二是针对循环流化床锅炉燃烧系统的控制特性，研究和运用常规PID改进方法或先进控制方法进行控制。

对于循环流化床机组，床温和床压是两个重要的运行参数。当机组的负荷指令发生变化时，一、二次风量和给煤量均发生变化，从而对床温床压产生较大影响。传统的单回路PID控制方法是通过外置床回料阀调节单侧床温，负荷变化时时床温控制是难题；另外通过调节冷渣器转速调节单侧床压，但这种调节方式十分缓慢，运行人员通常只把冷渣器转速作为稳态长期调节床压的手段，不能较好满足动态调节需求。

同时300MW主流循环流化床机组的双炉膛结构设计，两侧炉膛的均衡燃烧和床温床压平衡是实际运行中的难题，特别是两侧床压受一次风的影响很大，变工况时容易产生失稳翻床。因此传统的单侧炉膛单回路控制附加一定的双向偏差补偿的控制方法只能满足稳态和小范围变化的调节需求。

发明内容

要解决的技术问题：针对现有技术的不足，本发明提出一种双炉膛循环流化床机组的多变量协同控制方法，用于解决变工况时循环流化床的床温、床压以及负荷协调控制的稳定性差以及双炉膛均衡燃烧控制困难的技术问题。

技术方案：为解决上述技术问题，本发明采用以下技术方案：

一种双炉膛循环流化床机组的多变量协同控制方法，以左炉一次风量、左炉二次风量、左炉给煤量、左炉回料阀开度、主汽调门开度、右炉一次风量、右炉二次风量、右炉给煤量和右炉回料阀开度共9个参数作为输入量，同时以实发功率、机前压力、氧量、左炉床温、右炉床温、左炉床压和右炉床压共7个参数作为输出量，基于跟踪误差变化率的多变量约束DMC算法进行协同控制。

DMC(动态矩阵控制)算法是一种基于对象阶跃响应的预测控制算法，它适用于渐进稳定的线性对象。DMC算法不要求对模型的结构有先验知识，适用于竖向模型难以精确建立的复杂对象，因而该算法在工过程中比较实用。一般的DMC控制都包括以下三部分：DMC算法的模型预测、DMC算法的滚动优化以及DMC算法的反馈校正。

进一步，在本发明中，针对双炉膛循环流化床机组的多变量协同控制，对传统的DMC算法进行了特别的优化和改进，以u₁,u₂,...,u₉对应表示9个输入量并作为控制量，以y₁,y₂,...,y₇对应表示7个输出量，在目标函数中考虑跟踪误差变化率，并考虑输入和输出约束，整个控制方法具体包括如下步骤：

步骤一、获取双炉膛循环流化床机组的阶跃响应模型：

在稳态工况下，以9个输入量为输入对7个输出量进行阶跃响应实验，经滤波平滑后，得到7个输出量相对于9个输入量的阶跃响应系数分别为 $a_1^{\mathrm{ij}},a_2^{\mathrm{ij}},...,a_{N_{\mathrm{ij}}}^{\mathrm{ij}}(1\≤i\≤\mathrm{7,1}\≤j\≤9),$ 并满足 $a_{N_{\mathrm{ij}}+l}^{\mathrm{ij}}=a_{N_{\mathrm{ij}}}^{\mathrm{ij}}(l>0),$ 下标N_ij表示第j个输入对第i个输出的阶跃响应模型的截断长度，截断长度的选择满足保证输出量的阶跃响应值已经达到稳态；

以上述阶跃响应系数的组合构成第j个输入与第i个输出之间的阶跃响应模型；

步骤二、设置控制器的相关参数，包括采样时间T_s(采样时间一般取10～15s)、各输出量的预测时域P_i(1≤i≤7)、各输入量的控制时域M_j(1≤j≤9)、误差校正系数跟踪误差权矩阵控制增量权矩阵以及跟踪误差变化率权矩阵其中，预测时域P_i、控制时域M_j、误差校正系数和三个权矩阵需要在仿真中进行调试，并根据现场的运行数据再做调整；

步骤三、将步骤一得到的阶跃响应模型带入到下列(1)式中，求解出7个输出数量的预测值

$\hat{Y}\left(k\right)=Y_0\left(k\right)+\tilde{A}\mathrm{\ΔU}\left(k\right)---\left(1\right)$

(1)式中，

预测值表示为 $\hat{Y}\left(k\right)={[{\tilde{y}}_1{\left(k\right)}^T,{\tilde{y}}_2{\left(k\right)}^T,...,{\tilde{y}}_7{\left(k\right)}^T]}^T$

其中每个输入量的预测值表示为 ${\tilde{y}}_i\left(k\right)={[{\hat{y}}_i(k+1),{\hat{y}}_i(k+2),...,{\hat{y}}_i(k+P_i)]}^T$

Y₀(k)是在控制作用不变下的自由响应，并且有 $Y_0\left(k\right)={[{\tilde{y}}_{10}{\left(k\right)}^T,{\tilde{y}}_{20}{\left(k\right)}^T,...,{\tilde{y}}_{r0}{\left(k\right)}^T]}^T$

其中每个输入量的自由响应表示为 ${\tilde{y}}_{i0}\left(k\right)={[y_{i0}(k+1),y_{i0}(k+2),...,y_{i0}(k+P_i)]}^T$

是动态矩阵，并且

动态矩阵中每个子矩阵表示如下

控制增量ΔU(k)表示为 $\mathrm{\ΔU}\left(k\right)={[\mathrm{\Δ}{\tilde{u}}_1{\left(k\right)}^T,\mathrm{\Δ}{\tilde{u}}_2{\left(k\right)}^T,...,\mathrm{\Δ}{\tilde{u}}_9{\left(k\right)}^T]}^T$

其中每个控制增量表示为 $\mathrm{\Δ}{\tilde{u}}_j\left(k\right)={[\mathrm{\Δ}{u}_j(k+1),...,\mathrm{\Δ}{u}_j(k+M_j-1)]}^T;$

求解时，在0时刻，所有的预测值初值取0；

步骤四、求解自由响应Y₀(k)

首先计算出预测误差并进行反馈校正：

由于第i个输出量在k+1时刻的预测误差e_i(k+1)由下列(2)式表示

$e_i(k+1)=y_i(k+1)-{\hat{y}}_i(k+1)---\left(2\right)$

(2)式中，y_i(k+1)表示第i个输出量在k+1时刻的实际测量值；

所以第i个输出量在k+1时刻校正后的预测值为

${\tilde{y}}_{\mathrm{icor}}(k+1)={\hat{y}}_i(k+1)+{\tilde{h}}_i{e}_i(k+1)---\left(3\right)$

(3)式中，

误差校正系数表示为 ${\tilde{h}}_i={[h_{\mathrm{il}},h_{i2},...h_{{\mathrm{iP}}_i}]}^T;$

因此求解出所有输出量在k+1时刻校正后的预测值为

$Y_{\mathrm{cor}}(k+1)=\hat{Y}\left(k\right)+\tilde{H}E(k+1)---\left(4\right)$

(4)式中，

$Y_{\mathrm{cor}}(k+1)={[{\tilde{y}}_{1\mathrm{cor}}{(k+1)}^T,{\tilde{y}}_{2\mathrm{cor}}{(k+1)}^T,...,{\tilde{y}}_{7\mathrm{cor}}{(k+1)}^T]}^T$

E(k)＝[e₁(k+1)，e₂(k+1)，...，e₇(k+1)]^T

$\tilde{H}=\mathrm{diag}({\tilde{h}}_1,{\tilde{h}}_2,...,{\tilde{h}}_7)$

然后求解自由响应：

用校正后的预测值表示第i个输出量的自由响应的第a个分量，即y_i0(k+a)＝y_icor(k+a+1)a∈[1,P_i-1]   (2)

特别的，用y_icor(k+P_i)来近似的最后一项y_i0(k+P_i)，结合(5)式，表示为(6)式

${\tilde{y}}_{i0}\left(k\right)=S_i{\tilde{y}}_{\mathrm{icor}}(k+1)---\left(6\right)$

(6)式中

因此根据(6)式将所有输出量的自由响应表示为(7)式

$Y_0\left(k\right)=\tilde{S}{Y}_{\mathrm{cor}}(k+1)---\left(7\right)$

(7)式中 $\tilde{S}=\mathrm{diag}(S_1,S_2,...,S_7);$

步骤五、设定目标函数和约束，求解控制增量ΔU(k)，具体步骤如下：

首先将目标函数设定为

$J\left(k\right)=\left|\right|E_{\mathrm{pf}}\left(k\right)|{|}_{\tilde{Q}}^2+|\left|\mathrm{\Δ}{E}_{\mathrm{pf}}\left(k\right)\right|{|}_{\tilde{G}}^2+\left|\right|\mathrm{\ΔU}\left(k\right)|{|}_{\tilde{R}}^2---\left(8\right)$

(8)式中

E_pf(k)表示跟踪误差；ΔE_pf(k)表示跟踪误差变化率；

表示跟踪误差的惩罚权重；表示跟踪误差变化率的惩罚权重；表示控制增量的惩罚权重；它们的具体形式如下

$\tilde{Q}=\mathrm{diag}(Q_1,Q_2,...,Q_7){,Q}_i=\mathrm{diag}(q_1,q_2,...,q_{P_i})$

$\tilde{G}=\mathrm{diag}(G_1,G_2,...,G_7),G_i=\mathrm{diag}(g_1,g_2,...,g_{P_i})$

$\tilde{R}=\mathrm{diag}(R_1,R_2,...,R_9){,R}_i=\mathrm{diag}(r_1,r_2,...,r_{M_j})$

由于

$\begin{array}{c}\mathrm{\Δ}\hat{Y}\left(k\right)=\hat{Y}\left(k\right)-\hat{Y}(k-1)=Y_0+\tilde{A}\mathrm{\ΔU}\left(k\right)-[Y_0(k-1)+\tilde{A}\mathrm{\ΔU}(k-1)]\\ =\mathrm{\Δ}{Y}_0\left(k\right)+\tilde{A}[\mathrm{\ΔU}\left(k\right)-\mathrm{\ΔU}(k-1)]\end{array}---\left(9\right)$

并且对有Δu_j(k-1)＝0，所以

$\mathrm{\Δ}{\tilde{u}}_j(k-1)={[0,\mathrm{\Δ}{u}_j\left(k\right),...,\mathrm{\Δ}{u}_j(k+M_j-2)]}^T={\mathrm{\α}}_j\mathrm{\Δ}{\tilde{u}}_j\left(k\right)$

对所有的输出量，根据上式得到(10)式

ΔU(k-1)＝AΔU(k)   (10)

(10)式中

A＝diag(α₁,α₂,α₃,...α₉)

结合(9)式和(10)式，并令则将转化为用ΔU(k)表示的形式 $\mathrm{\Δ}\hat{Y}\left(k\right)=\mathrm{\Δ}{Y}_0\left(k\right)+\tilde{B}\mathrm{\ΔU}\left(k\right)$

因此将E_pf(k)和ΔE_pf(k)转化为用ΔU(k)表示的形式

$E_{\mathrm{pf}}\left(k\right)=W\left(k\right)-\hat{Y}\left(k\right)=W\left(k\right)-Y_0\left(k\right)-\tilde{A}\mathrm{\ΔU}\left(k\right)$

$\mathrm{\Δ}{E}_{\mathrm{pf}}\left(k\right)=\mathrm{\ΔW}\left(k\right)-\widehat{\mathrm{\ΔY}}\left(k\right)=\mathrm{\ΔW}\left(k\right)-\mathrm{\Δ}{Y}_0\left(k\right)-\tilde{B}\mathrm{\ΔU}\left(k\right)$

其中，W(k)表示各个输出量的设定值，前缀Δ表示各个量的变化值；

然后对输入量和输出量分别设置约束如下

输入量约束包括式(11)和式(12)

控制增量ΔU(k)满足 $\underset{\‾}{\mathrm{\ΔU}}\≤\mathrm{\ΔU}\left(k\right)\≤\stackrel{\‾}{\mathrm{\ΔU}}---\left(11\right)$

控制量U(k)满足 $\underset{\‾}{U}\≤U\left(k\right)\≤\stackrel{\‾}{U}$

式(11)中，ΔU表示控制增量的下限，表示控制增量的上限；

式(12)中，U表示控制量的下限，表示控制量的上限；

具体的控制增量和控制量的上下限由实际需求和现场设备的物理极限决定，没有限制时可取为无穷；

下面将式(12)中控制量U(k)的约束转化为(15)式中用控制增量ΔU(k)来表示的形式：

由于对于第j个输入量的控制量u_j有

u_j(k)＝u_j(k-1)+Δu_j(k)

u_j(k+1)＝u_j(k-1)+Δu_j(k)+Δu_j(k+1)

.

.

.

u_j(k+M_j-1)＝u_j(k-1)+Δu_j(k)+Δu_j(k+1)+…+Δu_j(k+M_j-1)

将上述控制量u_j的各式综合表示成(13)式的形式

${\tilde{u}}_j\left(k\right)=1_{M_j}{\tilde{u}}_j(k-1)+T_j\mathrm{\Δ}{\tilde{u}}_j\left(k\right)---\left(13\right)$

(13)式中

${\tilde{u}}_j\left(k\right)\left[\begin{array}{c}{u}_j\left(k\right)\\ u_j(k+1)\\ .\\ .\\ .\\ u_j(k+M_j-1)\end{array}\right]$ $l_{M_j}=\left[\begin{array}{c}1\\ 1\\ .\\ .\\ .\\ 1\end{array}\right]$

所以对所有输入量，基于(13)式可得

$U\left(k\right)=L_{9}U(k-1)+\tilde{T}\mathrm{\ΔU}\left(k\right)---\left(14\right)$

(14)式中

$L_9=\mathrm{diag}(l_{M_1},l_{M_2},...l_{M_9}),\tilde{T}=\mathrm{diag}(T_1,T_2,...,T_9)$

结合(12)式和(14)式得到(15)式

$\underset{\‾}{U}-L_9(k-1)\≤\tilde{T}\mathrm{\ΔU}\left(k\right)\≤\stackrel{\‾}{U}-L_{9}U(k-1)---\left(15\right)$

利用模型的预测值代替实际的输出量，(16)式表示输出量的约束

$\underset{\‾}{Y}\≤\hat{Y}\left(k\right)\≤\stackrel{\‾}{Y}---\left(16\right)$

(16)式中，Y表示实际输出量的下限，表示实际输出量的上限

结合(1)式和(16)式，将输出量的约束转化成为(17)式中用控制增量ΔU(k)来表示的形式： $\underset{\‾}{Y}-Y_0\left(k\right)\≤\tilde{A}\mathrm{\ΔU}\left(k\right)\≤\stackrel{\‾}{Y}-Y_0\left(k\right)---\left(17\right)$

结合式(11)、(15)和(17)，将输入量的约束和输出量的约束综合表示成(18)式的形式 $\tilde{R}\mathrm{\ΔU}\≤\tilde{c}---\left(18\right)$

(18)式中

$\tilde{R}=\left[\begin{array}{c}I\\ -I\\ \tilde{T}\\ -\tilde{T}\\ \tilde{A}\\ -\tilde{A}\end{array}\right]$ $\tilde{c}=\left[\begin{array}{c}\stackrel{\‾}{\mathrm{\ΔU}}\\ \underset{\‾}{-\mathrm{\ΔU}}\\ \stackrel{\‾}{U}-L_{9}U(k-1)\\ -\underset{\‾}{U}+L_{9}U(k-1)\\ \stackrel{\‾}{Y}-Y_0\left(k\right)\\ -\underset{\‾}{Y}+Y_0\left(k\right)\end{array}\right]$

将(8)式变形为如下二次型的形式

$\begin{array}{c}J=\left|\right|E_{\mathrm{pf}}{\left|\right|}_{\tilde{Q}}^2+\left|\right|\mathrm{\Δ}{E}_{\mathrm{pf}}{\left|\right|}_{\tilde{G}}^2+\left|\right|\mathrm{\ΔU}{\left|\right|}_{\tilde{R}}^2\\ ={(\tilde{A}\mathrm{\ΔU}+Y_0-W)}^T(\tilde{A}\mathrm{\ΔU}+Y_0-W)+{(\tilde{B}\mathrm{\ΔU}+\mathrm{\Δ}{Y}_0-{\mathrm{\ΔW}}_0)}^T(\tilde{B}\mathrm{\ΔU}+\mathrm{\Δ}{Y}_0-\mathrm{\Δ}{W}_0)+\mathrm{\Δ}{U}^T\mathrm{R\ΔU}\\ =\frac{1}{2}\mathrm{\Δ}{U}^T\tilde{H}\mathrm{\ΔU}+{\tilde{b}}^T\mathrm{\Δu}+f_0\end{array}$

二次型中

$\tilde{H}=2({\tilde{A}}^T\tilde{A}+{\tilde{B}}^T\tilde{B}+\tilde{R})$

$b^T=2{(Y_0-W)}^T\tilde{A}+2{(\mathrm{\Δ}{Y}_0-\mathrm{\ΔW})}^T\tilde{B}$

f₀＝(Y₀-W)^T(Y₀-W)+(ΔY₀-ΔW)^T(ΔY₀-ΔW)

由于f₀与控制增量ΔU(k)无关，在优化时忽略；将ΔU(k)的求解转化为标准的二次优化问题，形式如下

$\begin{array}{cc}\min & F\left(\mathrm{\ΔU}\right)=\frac{1}{2}{\mathrm{\ΔU}}^T\tilde{H}\mathrm{\ΔU}+{\tilde{b}}^T\mathrm{\ΔU}\\ S.t.& \tilde{R}\mathrm{\ΔU}\≤\tilde{c}\end{array}---\left(19\right)$

利用二次规划问题的有效集算法对(19)式进行求解，得到控制增量ΔU(k)；

步骤六、利用U(k)＝U(k-1)+ΔU(k)分别计算9个控制量的值；

步骤七、输出新的控制量，并更新(1)式中的预测值，并重复执行步骤三到步骤七。

有益效果：

本发明将一种多变量预测控制算法应用在双炉膛循环流化床机组上，创新过程包含两个方面：(1)提取构建多变量对象模型，并在DMC算法的目标函数引入跟踪误差变化率，同时引入输入输出约束，提出了基于跟踪误差变化率的多变量约束DMC算法；(2)首次将这样一种控制方法应用到300MW双炉膛循环流化床机组的研究。具体的，主要包括以下4个方面的优点：

1、提取构建一个循环流化床锅炉和汽轮机整个机组的9个输入和7个输出的对象模型，采用基于跟踪误差变化率的多变量约束DMC算法，同时控制9个输入量，实现了床温、床压和负荷的有机协调，克服了传统PID单回路控制分别调节解藕难的缺陷；

2、所采用的算法针对传统的DMC算法做出了改进，在目标函数中引入了跟踪误差变化率，提高了跟踪效果，因为传统的DMC算法目标函数中只包含跟踪误差，这样在求解的过程中只能保证最终稳定时，输出量和设定值没有偏差，而不能保证在调节过程中输出量紧紧跟随设定值，但引入跟踪误差变化率可解决这一问题；

3、对输入量、输出量附加了约束，并用二次优化算法进行求解，使得控制系统能在保证跟踪效果，使控制量波动最小的基础上，让求解的控制量在其物理范围内，避免了直接附加约束造成的经济性损失；

4、同步考虑了左右两个炉膛的床温和床压，能有效保证两个炉膛床温和床压平衡，实现双炉膛的均衡燃烧，并将功率、机前压力、氧量控制和床温、床压控制协同考虑，构建多变量预测控制模型，能有效保证机组变工况时，床温床压协同负荷跟踪变化。

附图说明

图1为本发明的控制系统的示意图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

在图1中，系统的输出Y中一共有7个输出量，分别是实发功率N_E、机前压力P_T、氧量O、左炉床温T_LF、右炉床温T_RF、左炉床压P_LF和右炉床压P_RF，前缀r表示他们的设定值；系统的输入U一共有输入量有9个，分别是左炉一次风量F_1,LF、左炉二次风量F_2,LF、左炉给煤量B_LF、左炉回料阀开度μ_v,LF、主汽调门开度μ_T、右炉一次风量F_1,RF、右炉二次风量F_2,RF、右炉给煤量B_RF和右炉回料阀开度μ_v,RF。这样，系统的输入可以表示为U＝[F_1,LF,F_2,LF,B_LF,μ_v,LF,μ_T,F_1,RF,F_2,RF,B_RF,μ_v,RF]^T，系统的输出可表示为Y＝[N_E,P_T,O,T_LF,T_RF,P_LF,P_RF]^T，系统的输出设定值可表示为W＝[r_N_E,r_P_T,r_O,r_T_LF,r_T_RF,r_P_TF,r_P_RF]^T，之后，按照前面所述从步骤一到步骤七，便是本发明的具体工作流程。

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

Claims (2)

1.一种双炉膛循环流化床机组的多变量协同控制方法，其特征在于：以左炉一次风量、左炉二次风量、左炉给煤量、左炉回料阀开度、主汽调门开度、右炉一次风量、右炉二次风量、右炉给煤量和右炉回料阀开度共9个参数作为输入量，同时以实发功率、机前压力、氧量、左炉床温、右炉床温、左炉床压和右炉床压共7个参数作为输出量，基于跟踪误差变化率的约束DMC算法进行协同控制。

2.根据权利要求1所述的双炉膛循环流化床机组的多变量协同控制方法，其特征在于：以u₁,u₂,…,u₉对应表示9个输入量并作为控制量，y₁,y₂,…,y₇对应表示7个输出量，在目标函数中考虑跟踪误差变化率，并考虑输入和输出约束，具体包括如下步骤：

步骤一、获取双炉膛循环流化床机组的阶跃响应模型

在稳态工况下，以9个输入量为输入对7个输出量进行阶跃响应实验，经滤波平滑后，得到7个输出量相对于9个输入量的阶跃响应系数分别为 $a_1^{\mathrm{ij}},a_2^{\mathrm{ij}},...,a_{N_{\mathrm{ij}}}^{\mathrm{ij}}(1\≤i\≤\mathrm{7,1}\≤j\≤9),$ 并满足 $a_{N_{\mathrm{ij}}+l}^{\mathrm{ij}}=a_{N_{\mathrm{ij}}}^{\mathrm{ij}}(l>0),$ 下标N_ij表示第j个输入对第i个输出的阶跃响应模型的截断长度，截断长度的选择满足保证输出量的阶跃响应值已经达到稳态；

以上述阶跃响应系数的组合构成第j个输入与第i个输出之间的阶跃响应模型；

步骤二、设置控制器的相关参数

包括采样时间T_s、各输出量的预测时域P_i(1≤i≤7)、各输入量的控制时域M_j(1≤j≤9)、误差校正系数跟踪误差权矩阵控制增量权矩阵以及跟踪误差变化率权矩阵

步骤三、将步骤一得到的阶跃响应模型带入到下列(1)式中，求解出7个输出数量的预测值

$\hat{Y}\left(k\right)=Y_0\left(k\right)+\tilde{A}\mathrm{\ΔU}\left(k\right)---\left(1\right)$

(1)式中，

预测值表示为 $\hat{Y}\left(k\right)={[{\tilde{y}}_1{\left(k\right)}^T,{\hat{y}}_2{\left(k\right)}^T,...,{\tilde{y}}_7{\left(k\right)}^T]}^T$

其中每个输入量的预测值表示为 ${\tilde{y}}_i\left(k\right)={[{\hat{y}}_i(k+1),{\hat{y}}_i(k+2),...,{\hat{y}}_i(k+P_i)]}^T$

Y₀(k)是在控制作用不变下的自由响应，并且有其中每个输入量的自由响应表示为 ${\tilde{y}}_{i0}\left(k\right)={[y_{i0}(k+1),y_{i0}(k+2),...,y_{i0}(k+P_i)]}^T$

其中所有自由响应的初值取0；

是动态矩阵，并且

动态矩阵中每个子矩阵表示如下

控制增量ΔU(k)表示为 $\mathrm{\ΔU}\left(k\right)={[{\mathrm{\Δ}\tilde{u}}_1{\left(k\right)}^T,\mathrm{\Δ}{\tilde{u}}_2{\left(k\right)}^T,...,\mathrm{\Δ}{\tilde{u}}_9{\left(k\right)}^T]}^T$

其中每个控制增量表示为 $\mathrm{\Δ}{\tilde{u}}_j\left(k\right)={[{\mathrm{\Δu}}_j(k+1),...,{\mathrm{\Δu}}_j(k+M_j-1)]}^T;$

求解时，在0时刻，所有的预测值初值取0；

步骤四、求解自由响应Y₀(k)

首先计算出预测误差并进行反馈校正：

由于第i个输出量在k+1时刻的预测误差e_i(k+1)由下列(2)式表示

$e_i(k+1)=y_i(k+1)-{\hat{y}}_i(k+1)---\left(2\right)$

(2)式中，y_i(k+1)表示第i个输出量在k+1时刻的实际测量值；

所以第i个输出量在k+1时刻校正后的预测值为

${\tilde{y}}_{\mathrm{icor}}(k+1)={\hat{y}}_i(k+1)+{\tilde{h}}_i{e}_i(k+1)---\left(3\right)$

(3)式中，

误差校正系数表示为 ${\tilde{h}}_i={[h_{i1},h_{i2},...h_{{\mathrm{iP}}_i}]}^T;$

因此求解出所有输出量在k+1时刻校正后的预测值为

$Y_{\mathrm{cor}}(k+1)=\hat{Y}\left(k\right)+\tilde{H}E(k+1)---\left(4\right)$

(4)式中，

$Y_{\mathrm{cor}}(k+1)={[{\tilde{y}}_{1\mathrm{cor}}{(k+1)}^T,{\tilde{y}}_{2\mathrm{cor}}{(k+1)}^T,...,{\tilde{y}}_{7\mathrm{cor}}{(k+1)}^T]}^T$

E(k)＝[e₁(k+1),e₂(k+1),...,e₇(k+1)]^T

$\tilde{H}=\mathrm{diag}({\tilde{h}}_1,{\tilde{h}}_2,...,{\tilde{h}}_7)$

然后求解自由响应：

用校正后的预测值表示第i个输出量的自由响应的第a个分量，即y_i0(k+a)＝y_icor(k+a+1)  a∈[1,P_i-1]   (1)

特别的，用y_icor(k+P_i)来近似的最后一项y_i0(k+P_i)，结合(5)式，表示为(6)式

${\tilde{y}}_{i0}\left(k\right)=S_i{\tilde{y}}_{\mathrm{icor}}(k+1)---\left(6\right)$

(6)式中

因此根据(6)式将所有输出量的自由响应表示为(7)式

$Y_0\left(k\right)=\tilde{S}{Y}_{\mathrm{cor}}(k+1)---\left(7\right)$

(7)式中 $\tilde{S}=\mathrm{diag}(S_1,S_2,...,S_7);$

步骤五、设定目标函数和约束，求解控制增量ΔU(k)，具体步骤如下

首先将目标函数设定为

$J\left(k\right)={\left|\right|E_{\mathrm{pf}}\left(k\right)\left|\right|}_{\tilde{Q}}^2+{\left|\right|{\mathrm{\ΔE}}_{\mathrm{pf}}\left(k\right)\left|\right|}_{\tilde{G}}^2+{\left|\right|\mathrm{\ΔU}\left(k\right)\left|\right|}_{\tilde{R}}^2---\left(8\right)$

(8)式中

E_pf(k)表示跟踪误差；ΔE_pf(k)表示跟踪误差变化率；

表示跟踪误差的惩罚权重；表示跟踪误差变化率的惩罚权重；表示控制增量的惩罚权重；它们的具体形式如下

$\begin{array}{cc}\tilde{Q}=\mathrm{diag}(Q_1,Q_2,...,Q_7)& Q_i=\mathrm{diag}(q_1,q_2,...,q_{P_i})\end{array}$

$\begin{array}{cc}\tilde{G}=\mathrm{diag}(G_1,G_2,...,G_7)& G_i=\mathrm{diag}(g_1,g_2,...,g_{P_i})\end{array}$

$\begin{array}{cc}\tilde{R}=\mathrm{diag}(R_1,R_2,...,R_7)& R_j=\mathrm{diag}(r_1,r_2,...,r_{M_j})\end{array}$

由于

$\begin{array}{c}\mathrm{\Δ}\hat{Y}\left(k\right)=\hat{Y}\left(k\right)-\hat{Y}(k-1)=Y_0+\tilde{A}\mathrm{\ΔU}\left(k\right)-[Y_0(k-1)+\tilde{A}\mathrm{\ΔU}(k-1)]\\ =\mathrm{\Δ}{Y}_0\left(k\right)+\tilde{A}[\mathrm{\ΔU}\left(k\right)-\mathrm{\ΔU}(k-1)]\end{array}---\left(9\right)$

并且对有Δu_j(k-1)＝0，所以

$\mathrm{\Δ}{\tilde{u}}_j(k-1)={[0,{\mathrm{\Δu}}_j\left(k\right),...,{\mathrm{\Δu}}_j(k+M_j-2)]}^T={\mathrm{\α}}_j\mathrm{\Δ}{\tilde{u}}_j\left(k\right)$

对所有的输出量，根据上式得到(10)式

ΔU(k-1)＝AΔU(k)                (10)

(10)式中

A＝diag(α₁,α₂,α₃,...α₉)

结合(9)式和(10)式，并令则将转化为用ΔU(k)表示的形式

$\mathrm{\Δ}\hat{Y}\left(k\right)={\mathrm{\ΔY}}_0\left(k\right)+\tilde{B}\mathrm{\ΔU}\left(k\right)$

因此将E_pf(k)和ΔE_pf(k)转化为用ΔU(k)表示的形式

$E_{\mathrm{pf}}\left(k\right)=W\left(k\right)-\hat{Y}\left(k\right)=W\left(k\right)-Y_0\left(k\right)-\tilde{A}\mathrm{\ΔU}\left(k\right)$

${\mathrm{\ΔE}}_{\mathrm{pf}}\left(k\right)=\mathrm{\ΔW}\left(k\right)-\widehat{\mathrm{\ΔY}}\left(k\right)=\mathrm{\ΔW}\left(k\right)-{\mathrm{\ΔY}}_0\left(k\right)-\tilde{B}\mathrm{\ΔU}\left(k\right)$

其中，W(k)表示各个输出量的设定值，前缀Δ表示各个量的变化值；

然后对输入量和输出量分别设置约束如下

输入量约束包括式(11)和式(12)

控制增量ΔU(k)满足 $\underset{\‾}{\mathrm{\ΔU}}\≤\mathrm{\ΔU}\left(k\right)\≤\stackrel{\‾}{\mathrm{\ΔU}}---\left(11\right)$

控制量U(k)满足 $\underset{\‾}{U}\≤U\left(k\right)\≤\stackrel{\‾}{U}---\left(12\right)$

式(11)中，ΔU表示控制增量的下限，表示控制增量的上限；

式(12)中，U表示控制量的下限，表示控制量的上限；

下面将式(12)中控制量U(k)的约束转化为(15)式中用控制增量ΔU(k)来表示的形式：

由于对于第j个输入量的控制量u_j有

u_j(k)＝u_j(k-1)+Δu_j(k)

u_j(k+1)＝u_j(k-1)+Δu_j(k)+Δu_j(k+1)

·

·

·

u_j(k+M_j-1)＝u_j(k-1)+Δu_j(k)+Δu_j(k+1)+…+Δu_j(k+M_j-1)

将上述控制量u_j的各式综合表示成(13)式的形式

${\tilde{u}}_j\left(k\right)=1_{M_j}{\tilde{u}}_j(k-1)+T_j\mathrm{\Δ}{\tilde{u}}_j\left(k\right)---\left(13\right)$

(13)式中

$\begin{array}{cc}{\tilde{u}}_j\left(k\right)\left[\begin{array}{c}{u}_j\left(k\right)\\ u_j(k+1)\\ \·\\ \·\\ \·\\ u_j(k+M_j-1)\end{array}\right]& l_{M_j}=\left[\begin{array}{c}1\\ 1\\ \·\\ \·\\ \·\\ 1\end{array}\right]\end{array}$

所以对所有输入量，基于(13)式可得

$\left(k\right)L_{9}U(k-1)+\tilde{T}\left(k\right)\left(14\right)$

(14)式中

$\begin{array}{cc}{L}_9=\mathrm{diag}(l_{M_1},l_{M_2},...l_{M_9})& \tilde{T}\end{array},=\mathrm{diag}(T_1,T_2,...,T_9)$

结合(12)式和(14)式得到(15)式

$\underset{\‾}{U}-L_{9}U(k-1)\≤\tilde{T}\mathrm{\ΔU}\left(k\right)\≤\stackrel{\‾}{U}-L_{9}U(k-1)---\left(15\right)$

利用模型的预测值代替实际的输出量，(16)式表示输出量的约束

$\underset{\‾}{Y}\≤\hat{Y}\left(k\right)\≤\stackrel{\‾}{Y}---\left(16\right)$

(16)式中，Y表示实际输出量的下限，表示实际输出量的上限

结合(1)式和(16)式，将输出量的约束转化成为(17)式中用控制增量ΔU(k)来表示的形式： $\underset{\‾}{Y}-Y_0\left(k\right)\≤\tilde{A}\mathrm{\ΔU}\left(k\right)\≤\stackrel{\‾}{Y}-Y_0\left(k\right)---\left(17\right)$

结合式(11)、(15)和(17)，将输入量的约束和输出量的约束综合表示成(18)式的形式 $\tilde{R}\mathrm{\ΔU}\≤\tilde{c}---\left(18\right)$

(18)式中

$\begin{array}{cc}\tilde{R}=\left[\begin{array}{c}I\\ -I\\ \tilde{T}\\ -\tilde{T}\\ \tilde{A}\\ -\tilde{A}\end{array}\right]& \tilde{c}\end{array}=\left[\begin{array}{c}\stackrel{\‾}{\mathrm{\ΔU}}\\ -\underset{\‾}{\mathrm{\ΔU}}\\ \stackrel{\‾}{U}-L_{9}U(k-1)\\ -\underset{\‾}{U}+L_{9}U(k-1)\\ \stackrel{\‾}{Y}-Y_0\left(k\right)\\ -\underset{\‾}{Y}+Y_0\left(k\right)\end{array}\right]$

将(8)式变形为如下二次型的形式

$\begin{array}{c}J={\left|\right|E_{\mathrm{pf}}\left|\right|}_{\tilde{Q}}^2+{\left|\right|\mathrm{\Δ}{E}_{\mathrm{pf}}\left|\right|}_{\tilde{G}}^2+{\left|\right|\mathrm{\ΔU}\left|\right|}_{\tilde{R}}^2\\ ={(\tilde{A}\mathrm{\ΔU}+Y_0-W)}^T(\tilde{A}\mathrm{\ΔU}+Y_0-W)+{(\tilde{B}\mathrm{\ΔU}+{\mathrm{\ΔY}}_0-{\mathrm{\ΔW}}_0)}^T(\tilde{B}\mathrm{\ΔU}+{\mathrm{\ΔY}}_0-{\mathrm{\ΔW}}_0)+{\mathrm{\ΔU}}^T\mathrm{R\ΔU}\\ =\frac{1}{2}{\mathrm{\ΔU}}^T\tilde{H}\mathrm{\ΔU}+{\tilde{b}}^T\mathrm{\ΔU}+f_0\end{array}$

二次型中

$\tilde{H}=2({\tilde{A}}^T\tilde{A}+{\tilde{B}}^T\tilde{B}+\tilde{R})$

$b^T=2{(Y_0-W)}^T\tilde{A}+2{(\mathrm{\Δ}{Y}_0-\mathrm{\ΔW})}^T\tilde{B}$

f₀＝(Y₀-W)^T(Y₀-W)+(ΔY₀-ΔW)^T(ΔY₀-ΔW)

由于f₀与控制增量ΔU(k)无关，在优化时忽略；将ΔU(k)的求解转化为标准的二次优化问题，形式如下

$\begin{array}{cc}\min & F\left(\mathrm{\ΔU}\right)=\frac{1}{2}{\mathrm{\ΔU}}^T\tilde{H}\mathrm{\ΔU}+{\tilde{b}}^T\mathrm{\ΔU}\\ S.t.& \tilde{R}\mathrm{\ΔU}\≤\tilde{c}\end{array}---\left(19\right)$

利用二次规划问题的有效集算法对(19)式进行求解，得到控制增量ΔU(k)；

步骤六、利用U(k)＝U(k-1)+ΔU(k)分别计算9个控制量的值；

步骤七、输出新的控制量，并更新(1)式中的预测值，并重复执行步骤三到步骤七。