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

摘要

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

主权项

一种双炉膛循环流化床机组的多变量协同控制方法，其特征在于：以u₁,u₂,…,u₉对应表示左炉一次风量、左炉二次风量、左炉给煤量、左炉回料阀开度、主汽调门开度、右炉一次风量、右炉二次风量、右炉给煤量和右炉回料阀开度9个输入量并作为控制量，y₁,y₂,…,y₇对应表示实发功率、机前压力、氧量、左炉床温、右炉床温、左炉床压和右炉床压7个输出量，在目标函数中考虑跟踪误差变化率，并考虑输入和输出约束，基于跟踪误差变化率的约束DMC算法进行协同控制，具体包括如下步骤：步骤一、获取双炉膛循环流化床机组的阶跃响应模型在稳态工况下，以9个输入量为输入对7个输出量进行阶跃响应实验，经滤波平滑后，得到7个输出量相对于9个输入量的阶跃响应系数分别为$a_1^{\mathrm{ij}},a_2^{\mathrm{ij}},...,a_{N_{\mathrm{ij}}}^{\mathrm{ij}}(1\le i\le \mathrm{7,1}\le j\le 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{\Delta 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)是在控制作用不变下的自由响应，并且有其中每个输入量的自由响应表示为${\tilde{y}}_{i0}\left(k\right)={[y_{i0}(k+1),y_{i0}(k+2),...,y_{i0}(k+P_i)]}^T$其中所有自由响应的初值取0；是动态矩阵，并且动态矩阵中每个子矩阵表示如下控制增量ΔU(k)表示为$\mathrm{\Delta U}\left(k\right)={[\Delta {\tilde{u}}_1{\left(k\right)}^T,\Delta {\tilde{u}}_2{\left(k\right)}^T,...,\Delta {\tilde{u}}_9{\left(k\right)}^T]}^T$其中每个控制增量表示为$\Delta {\tilde{u}}_j\left(k\right)={[\Delta u_j(k+1),...,\Delta 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+1)＝[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|\Delta E_{\mathrm{pf}}\left(k\right)\left|\right|}_{\tilde{G}}^2+{\left|\right|\mathrm{\Delta 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_9)& R_j=\mathrm{diag}(r_1,r_2,...,r_{M_j})\end{array}$由于$\begin{array}{c}\Delta \hat{Y}\left(k\right)=\hat{Y}\left(k\right)-\hat{Y}(k-1)=Y_0+\tilde{A}\mathrm{\Delta U}\left(k\right)-[Y_0(k-1)+\tilde{A}\mathrm{\Delta U}(k-1)]\\ =\Delta Y_0\left(k\right)+\tilde{A}[\mathrm{\Delta U}\left(k\right)-\mathrm{\Delta U}(k-1)]\end{array}---\left(9\right)$并且对有Δu_j(k‑1)＝0，所以$\Delta {\tilde{u}}_j(k-1)={[0,\Delta u_j\left(k\right),...,\Delta u_j(k+M_j-2)]}^T={\alpha }_j\Delta {\tilde{u}}_j\left(k\right)$对所有的输出量，根据上式得到(10)式ΔU(k‑1)＝AΔU(k)           (10)(10)式中A＝diag(α₁,α₂,α₃,...α₉)结合(9)式和(10)式，并令则将转化为用ΔU(k)表示的形式$\Delta \hat{Y}\left(k\right)=\Delta Y_0\left(k\right)+\tilde{B}\mathrm{\Delta 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{\Delta U}\left(k\right)$$\Delta E_{\mathrm{pf}}\left(k\right)=\mathrm{\Delta W}\left(k\right)-\widehat{\mathrm{\Delta Y}}\left(k\right)=\mathrm{\Delta W}\left(k\right)-\Delta Y_0\left(k\right)-\tilde{B}\mathrm{\Delta U}\left(k\right)$其中，W(k)表示各个输出量的设定值，前缀Δ表示各个量的变化值；然后对输入量和输出量分别设置约束如下输入量约束包括式(11)和式(12)控制增量ΔU(k)满足$\underset{‾}{\mathrm{\Delta U}}\le \mathrm{\Delta U}\left(k\right)\le \overline{\mathrm{\Delta U}}---\left(11\right)$控制量U(k)满足$\underset{‾}{U}\le U\left(k\right)\le \bar{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)=l_{M_j}{\tilde{u}}_j(k-1)+T_j\Delta {\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]}_{M_j\times 1}$其中，是M_j×1维的向量，T_j是M_j×M_j的矩阵。所以对所有输入量，基于(13)式可得$U\left(k\right)=L_{9}U(k-1)+\tilde{T}\mathrm{\Delta U}\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}=\mathrm{diag}(T_1,T_2,...,T_9)\end{array}$结合(12)式和(14)式得到(15)式$\underset{‾}{U}-L_{9}U(k-1)\le \tilde{T}\mathrm{\Delta U}\left(k\right)\le \bar{U}-L_{9}U(k-1)---\left(15\right)$利用模型的预测值代替实际的输出量，(16)式表示输出量的约束$\underset{‾}{Y}\le \hat{Y}\left(k\right)\le \bar{Y}---\left(16\right)$(16)式中，Y表示实际输出量的下限，表示实际输出量的上限结合(1)式和(16)式，将输出量的约束转化成为(17)式中用控制增量ΔU(k)来表示的形式：$\underset{‾}{Y}-Y_0\left(k\right)\le \tilde{A}\mathrm{\Delta U}\left(k\right)\le \bar{Y}-Y_0\left(k\right)---\left(17\right)$结合式(11)、(15)和(17)，将输入量的约束和输出量的约束综合表示成(18)式的形式$\tilde{R}\mathrm{\Delta U}\le \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}\overline{\mathrm{\Delta U}}\\ -\underset{‾}{\mathrm{\Delta U}}\\ \bar{U}-L_{9}U(k-1)\\ -\underset{‾}{U}+L_{9}U(k-1)\\ \bar{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|\Delta E_{\mathrm{pf}}\left|\right|}_{\tilde{G}}^2+{\left|\right|\mathrm{\Delta U}\left|\right|}_{\tilde{R}}^2\\ ={(\tilde{A}\mathrm{\Delta U}+Y_0-W)}^T(\tilde{A}\mathrm{\Delta U}+Y_0-W)+{(\tilde{B}\mathrm{\Delta U}+\Delta Y_0-\Delta W_0)}^T(\tilde{B}\mathrm{\Delta U}+\Delta Y_0-\Delta W_0)+\Delta U^T\mathrm{R\Delta U}\\ =\frac{1}{2}\Delta U^T\tilde{H}\mathrm{\Delta U}+{\tilde{b}}^T\mathrm{\Delta 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{(\Delta Y_0-\mathrm{\Delta W})}^T\tilde{B}$$f_0={(Y_0-W)}^T(Y_0-W)+{(\Delta Y_0-\mathrm{\Delta W})}^T(\Delta Y_0-\mathrm{\Delta W})$由于f₀与控制增量ΔU(k)无关，在优化时忽略；将ΔU(k)的求解转化为标准的二次优化问题，形式如下$\begin{array}{cc}\min & F\left(\mathrm{\Delta U}\right)=\frac{1}{2}\Delta U^T\tilde{H}\mathrm{\Delta U}+{\tilde{b}}^T\mathrm{\Delta U}\end{array}---\left(19\right)$$\begin{array}{cc}S.t.& \tilde{R}\mathrm{\Delta U}\le \tilde{c}\end{array}$利用二次规划问题的有效集算法对(19)式进行求解，得到控制增量ΔU(k)；步骤六、利用U(k)＝U(k‑1)+ΔU(k)分别计算9个控制量的值；步骤七、输出新的控制量，并更新(1)式中的预测值，并重复执行步骤三到步骤七。