CN103413029B - 具有多速率采样连续搅拌釜式反应器的滚动时域估计方法 - Google Patents

Abstract

具有多速率采样连续搅拌釜式反应器的滚动时域估计方法，步骤如下：对CSTR装置进行状态空间建模，并确定其正常工作的工作范围和工作范围内的稳态工作点；在工作范围内的稳态工作点附近将CSTR的状态空间模型线性化；设定系统的标准采样周期和各个传感器的采样周期，将线性化的CSTR状态空间模型离散化，对于未采样的测量输出用预测值补偿，得到多速率的CSTR线性离散化模型；设定滚动时域窗口长度和权重矩阵，将多速率CSTR的滚动时域估计问题转化为等价的最小化问题；通过一阶最优性原理求解滚动优化问题，得到最优估计器。本发明模型合理，且具有补偿未采样数据和在线计算功能，使得估计值更加准确。

Description

具有多速率采样连续搅拌釜式反应器的滚动时域估计方法

技术领域

本发明涉及一种多速率采样系统的滚动时域估计方法，尤其涉及一种多速率采样连续搅拌釜式反应器的滚动时域估计方法。

背景技术

连续搅拌釜式反应器（Continuous Stirred Tank Reactor，简称为CSTR）是聚合化学反应中广泛使用的一种反应器。CSTR在化工生产的核心设备中占有相当重要的地位，在染料、医药试剂、食品及合成材料工业中，CSTR得到了广泛的应用。

在早期反应釜的自动控制中，一般采用单元组合仪表组成的位置式控制装置，但是化学反应过程一般都有很强的非线性和时滞性，采用这种简单的控制方式很难达到理想的控制精度。随着计算机技术的发展，越来越多的化学反应采用计算机控制系统。计算机控制系统需要使用传感器采集离散信息，设计离散控制器。在聚合化学反应中，反应釜内的温度和压力能够通过传感器测量快速得到。反应物的分子量及浓度则采用凝胶渗透色谱法得到。相比于温度和压力的测量周期，分子量和浓度的测量周期将会更长。所以，CSTR中的对温度和浓度的测量采样频率是不相同的，即多速率采样。对于多速率采样，在不同的采样时刻得到的测量输出数据量是不同的，并会导致测量信息不完整。因此需要利用有限的测量信息来估计系统状态，反应系统的实时性。

发明内容

本发明的目的在于提供一种适用于具有多速率采样连续搅拌釜式反应器的滚 动时域估计方法。

本发明解决其技术问题所采用的技术方案是：

一种具有多速率采样连续搅拌釜式反应器的滚动时域估计方法，具体步骤如下：

（1）、对CSTR装置进行状态空间建模，并确定其正常工作的工作范围和工作范围内的稳态工作点。

（2）、在工作范围内的稳态工作点附近将CSTR的状态空间模型线性化。

（3）、设定系统的标准采样周期和各个传感器的采样周期，将线性化的CSTR状态空间模型离散化，对于未采样的测量输出用预测值补偿，得到多速率的CSTR线性离散化状态空间模型。

（4）、设定滚动时域窗口长度N和权重矩阵，将多速率CSTR的滚动时域估计问题转化为等价的最小化问题。

（5）、通过一阶最优性原理求解步骤（4）设定的最小化问题，具体步骤如下：

S1-1：初始化，设定测试时间长度K，在可行域的区间范围内，任意初始化k时刻的先验估计值，k-N时刻到k时刻的测量输出序列；

S1-2：根据一阶最优性原理，以k时刻的先验估计值为初始迭代点，对步骤（4）中的最小化问题求一阶偏导数，得到k-N时刻最优估计值；

S1-3：根据滚动优化原理，计算当前k时刻的最优估计值；

S1-4：根据k时刻的最优估计值更新k+1时刻的先验估计值；

S1-5：判断终止条件：如果k=K，结束，得到信噪比估计最优值；否则，k=k+1，转到S1-2。

本发明的技术构思为：本发明考虑了CSTR内聚合反应过程中对反应釜内的温度和反应物浓度测量频率的差异，给出了一种具有多速率采样的CSTR线性离 散化状态空间模型，设计了基于滚动时域估计方法的CSTR状态估计器，给出了反应釜内温度和反应物浓度的最优估计值。

从上述技术方案可以看出，本发明的有益效果主要表现在：多速率CSTR滚动时域估计方法，与已有的估计方法相比，滚动时域估计方法能够滚动优化和在线计算，并对传感器没有采样的时刻用预测值替代起到很好的补偿作用，从而能够更加准确地给出反应釜内各参数的值。

附图说明

图1是本发明实施例中CSTR示意图。

图2是本发明实施例中求解最小化问题的流程图。

图3是本发明实施例中，采用本发明方法的效果图。

具体实施方式

为使本发明的目的、技术方案和优点更加清晰，下面结合附图和实施例对本发明的技术方案作进一步描述。

参照图1～图3，一种具有多速率采样CSTR的滚动时域估计方法，将本发明提出的滚动时域估计方法用于一阶并行反应CSTR的状态估计，其目的是估计出反应釜内的反应物浓度和温度。一阶并行反应，即往CSTR中加入某物质A，进入反应釜后，A发生化学反应生成物B，同时由于反应釜内的强反应，部分物质B生成物质C，并达到化学平衡的过程。接下来介绍具体实施步骤：

（1）、对附图1所示一阶并行反应CSTR装置进行状态空间建模，并确定其正常工作的工作范围和工作范围内的稳态工作点。如附图1所示，建立一阶并行反应CSTR的状态空间模型如下所示：

$\left\{\begin{array}{c}\frac{{\mathrm{dC}}_A}{\mathrm{dt}}=\frac{F}{V}(C_{A0}-C_A)-k_1{C}_A\\ \frac{{\mathrm{dC}}_B}{\mathrm{dt}}=-\frac{F}{V}{C}_B+k_1{C}_A-k_2{C}_B\\ \frac{\mathrm{d\θ}}{\mathrm{dt}}=\frac{F}{V}({\mathrm{\θ}}_0-\mathrm{\θ})+\frac{k_w{A}_R}{\mathrm{\ρ}{C}_{P}V}({\mathrm{\θ}}_k-\mathrm{\θ})-\frac{k_1{C}_A{\mathrm{\ΔH}}_R^{\mathrm{AB}}+k_2{C}_B{\mathrm{\ΔH}}_R^{\mathrm{BC}}}{\mathrm{\ρ}{C}_P}\end{array}\right.---\left(1\right)$

式中，C_A为物质A的浓度，C_A0为物质A初始浓度，C_As为物质A在稳态时浓度，C_B为物质B的浓度，C_Bs为物质B在稳态时浓度，θ为反应釜内温度，θ₀为反应釜内初始温度，θ_s为稳态时反应釜内温度，θ_k为冷却剂温度，F/V为稀释率，V为体积流量，A_R为反应堆表面面积，C_P为热容量，k_w为传热系数，ρ为密度，为物质A到物质B反应焓，为物质B到物质C反应焓，反应速率系数k₁和k₂由反应釜内温度决定i=1,2，k₀为频率因子，E_A1和E_A2为活化能，R为理想气体常数。CSTR的工作范围和稳态工作点取值如表1所示：

表1.CSTR模型参数和稳态工作点

（2）、根据一阶并行反应CSTR的状态空间模型，在工作范围内的稳态工作点附近将CSTR的状态空间模型线性化，得到线性化状态空间模型如下所示：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+w\left(t\right)---\left(2\right)$

式中，x(t)=[x₁(t) x₂(t) x₃(t)]^T，x₁(t)为物质A在t时刻的浓度，x₂(t)为物质B 在t时刻的浓度，x₃(t)为在t时刻反应釜内的温度，w(t)为在t时刻反应釜内的有界扰动即||w(t)||≤0.3，A为线性化后得到的系统参数

$A=\left[\begin{array}{ccc}-\frac{F}{V}-k_1& 0& \frac{E_{A1}}{{\mathrm{R\θ}}_s^2}{k}_1{C}_{\mathrm{As}}\\ k_1& -\frac{F}{V}-k_2& -\frac{E_{A1}}{{\mathrm{R\θ}}_s^2}{k}_1{C}_{\mathrm{As}}+\frac{E_{A2}}{{\mathrm{R\θ}}_s^2}{k}_2{C}_{\mathrm{Bs}}\\ -\frac{k_1{\mathrm{\ΔH}}_R^{\mathrm{AB}}}{{\mathrm{\ρC}}_P}& -\frac{k_2{\mathrm{\ΔH}}_R^{\mathrm{BC}}}{\mathrm{\ρ}{C}_P}& -\frac{F}{V}-\frac{k_w{A}_R}{{\mathrm{\ρC}}_{P}V}+\frac{E_{A1}{k}_1{C}_{\mathrm{As}}{\mathrm{\ΔH}}_R^{\mathrm{AB}}+E_{A2}{k}_2{C}_{\mathrm{Bs}}{\mathrm{\ΔH}}_R^{\mathrm{BC}}}{{\mathrm{R\θ}}_s^2{\mathrm{\ρC}}_P}\end{array}\right].$

$=\left[\begin{array}{ccc}-0.9388& 0& 0.0459\\ 0.625& -0.9388& -0.0125\\ -0.9335& 2.4449& -0.8894\end{array}\right]$

（3）、设定系统的标准采样周期T₀=1min，CSTR测量输出量为物质B的浓度和反应釜内温度，对物质A浓度的测量周期为T₁=2min，对反应釜内温度的测量周期为T₂=1min。以标准采样周期T₀将线性化的CSTR状态空间模型离散化，得到离散状态空间模型为

$x(k+1)=\stackrel{\‾}{A}x\left(k\right)+w\left(k\right)---\left(3\right)$

式中， $\stackrel{\‾}{A}=e^{A{T}_0}=\left[\begin{array}{ccc}0.3872& 0.0222& 0.0182\\ 0.2444& 0.3897& 0.0007\\ -0.0685& 0.9711& 0.4008\end{array}\right].$

由于对物质A浓度的测量为2倍的标准采样周期，从而会导致未采样时刻得不到测量值，这种情况下可将未采样时刻的数据视为测量信息丢失，并采用预测值替代进行输出补偿。由此可得CSTR的测量输出方程如下所示：

$y\left(k\right)=\mathrm{\θ}\left(k\right)[\mathrm{Cx}\left(k\right)+v\left(k\right)]+[I-\mathrm{\θ}\left(k\right)]\stackrel{\‾}{y}\left(k\right)---\left(4\right)$

式中，y(k)=[y₁(k) y₂(k)]^T，y₁(k)为物质A在k时刻浓度的测量值，y₂(k)为在k时刻反应釜内温度的测量值，v(k)为在k时刻反应釜内的有界扰动即||v(k)||≤0.2， 为k时刻系统状态x(k)的预测值， $C=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\end{array}\right]$ 为测量输出 权重矩阵， $I=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ 为单位矩阵， $\mathrm{\θ}\left(k\right)=\left[\begin{array}{cc}{\mathrm{\θ}}_1\left(k\right)& 0\\ 0& {\mathrm{\θ}}_2\left(k\right)\end{array}\right],$ θ₁(k)=1表示物质A在k时刻有测量值，θ₁(k)=0表示物质A在k时刻没有测量值，θ₂(k)=1。

（4）、设定滚动时域窗口长度N=5和权重矩阵μ=0.1，多速率CSTR的滚动时域估计问题如下所示：

$\underset{\hat{x}(k-N)}{\min }J\left(k\right)---\left(5\right)$

约束条件：

$J\left(k\right)={\left|\right|\hat{x}(k-N)-\stackrel{\‾}{x}(k-N)\left|\right|}_{\mathrm{\μ}}^2+\underset{i=k-N}{\overset{k}{\mathrm{\Σ}}}{\left|\right|y\left(i\right)-\mathrm{\θ}\left(i\right)[C\hat{x}\left(i\right)+v\left(i\right)]+[I-\mathrm{\θ}\left(i\right)]\stackrel{\‾}{y}\left(i\right)\left|\right|}^2$

$\hat{x}(i+1)=A\hat{x}\left(i\right),i=k-N,\·\·\·,k-1$

$\stackrel{\‾}{x}(k-N)=A\hat{x}(k-N-1),k=N+1,N+2,\·\·\·$

式中，J(k)为性能指标，为k时刻x(k)的最优估计值。

（5）、通过一阶最优性原理求解步骤（4）设定的最小化问题（5），具体步骤如下：

S1-1：初始化，设定测试时间长度K=80，在可行域的区间范围内，任意初始化k时刻的先验估计值 $\stackrel{\‾}{x}\left(0\right)={\left[\begin{array}{ccc}1& 1& 0\end{array}\right]}^T$ 和系统状态初始值x(0)=[0.5 0.5 0.5]^T，k-5时刻到k时刻的测量输出序列 $y_{k-5}^k=\left[\begin{array}{c}y(k-5)\\ .\\ .\\ .\\ y(k-1)\\ y\left(k\right)\end{array}\right];$

S1-2：根据一阶最优性原理，以k时刻的先验估计值为初始迭代点，对步骤（4）中的最小化问题求一阶偏导数，得到k-5时刻的最优估计值如下所示

$\hat{x}(k-5)={({\mathrm{\μI}}_3+G_1)}^{-1}(\mathrm{\μ}\stackrel{\‾}{x}(k-5)+G_2{y}_{k-5}^k)$

式中， $G_1=F_5^T{F}_5-F_{15}^T(I_5-\mathrm{\Θ}\left(k\right))F_{15},$ $G_2=F_5^T-F_{15}^T(I_5-\mathrm{\Θ}\left(k\right)),$ I₃为3维的单位矩 阵，I₅为12维的单位矩阵，Θ(k)=diag{θ(k-5),…,θ(k-1),θ(k)}，

$F_5=\left[\begin{array}{c}C\\ \mathrm{CA}\\ .\\ .\\ .\\ {\mathrm{CA}}^5\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0.3872& 0.0222& 0.0182\\ -0.0685& 0.9711& 0.4008\\ 0.1541& 0.035& 0.0144\\ 0.1834& 0.7661& 0.1601\\ 0.0672& 0.031& 0.0086\\ 0.2473& 0.4581& 0.0681\\ 0.033& 0.0219& 0.0047\\ 0.2031& 0.2501& 0.0321\\ 0.0178& 0.0138& 0.0025\\ 0.1376& 0.1331& 0.0168\end{array}\right],$ $F_{15}=\left[\begin{array}{c}C\\ 0\\ .\\ .\\ .\\ 0\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right];$

S1-3：根据滚动优化原理，计算当前k时刻的最优估计值

S1-4：根据k时刻的最优估计值更新k+1时刻的先验估计值 $\stackrel{\‾}{x}(k-4)=A\hat{x}(k-5);$

S1-5：判断终止条件：如果k=K，结束，得到信噪比估计最优值；否则，k=k+1，转到S1-2。

采用所述步骤，计算80个采样时刻得到的结果如图3所示，其中，图3中的(a)为CSTR物质A的浓度和其估计值，横坐标为采样次数，纵坐标为物质A的浓度和其估计值；图3中的(b)为CSTR物质B的浓度和其估计值，横坐标为采样次数，纵坐标为物质B的浓度和其估计值；图3中的(c)为CSTR内温度和其估计值，横坐标为采样次数，纵坐标为温度和其估计值。从图3可以看出，针对具有多速率采样的CSTR，滚动时域估计方法对未采样部分信息缺失也能够很好的估计出系统的状态。并且滚动时域估计方法具有滚动优化和在线计算的优点，从而能够更加准确地给出CSTR反应物浓度和反应器温度。

以上阐述的是本发明给出的实例表现出的优良估计效果。需要指出的是，本发明不只限于上述实施例，对于其他化工聚合反应多速率采样估计问题，采用本 发明给出的方法设计滚动时域估计器，均能给出系统的状态估计值。

Claims (1)

1.具有多速率采样连续搅拌釜式反应器的滚动时域估计方法，其特征在于，具体步骤如下：

(1)、对CSTR装置进行状态空间建模，并确定其正常工作的工作范围和工作范围内的稳态工作点；

式中，C_A为物质A的浓度，C_A0为物质A初始浓度，C_As为物质A在稳态时浓度，C_B为物质B的浓度，C_Bs为物质B在稳态时浓度，为反应釜内温度，为反应釜内初始温度，为稳态时反应釜内温度，为冷却剂温度，F/V为稀释率，V为体积流量，A_R为反应堆表面面积，C_P为热容量，k_w为传热系数，ρ为密度，为物质A到物质B反应焓，为物质B到物质C反应焓，反应速率系数k₁和k₂由反应釜内温度决定k₀为频率因子，E_A1和E_A2为活化能，R为理想气体常数；

(2)、在工作范围内的稳态工作点附近将CSTR的状态空间模型线性化；

$\stackrel{\·}{x}\left(t\right)=Ax\left(t\right)+w\left(t\right)$

式中，x(t)＝[x₁(t) x₂(t) x₃(t)]^T，x₁(t)为物质A在t时刻的浓度，x₂(t)为物质B在t时刻的浓度，x₃(t)为在t时刻反应釜内的温度，w(t)为在t时刻反应釜内的有界扰动，即||w(t)||≤0.3，A为线性化后得到的系统参数，具体为

(3)、设定系统的标准采样周期T₀和各个传感器的采样周期，对CSTR中温度和浓度的测量传感器采用不同的采样周期，将线性化的CSTR状态空间模型离散化，对于未采样的测量输出用预测值补偿，得到多速率的CSTR线性离散化状态空间模型；

$\left\{\begin{array}{c}x(k+1)=\stackrel{\‾}{A}x\left(k\right)+w\left(k\right)\\ y\left(k\right)=\mathrm{\θ}\left(k\right)\[Cx\left(k\right)+v\left(k\right)\]+\[I-\mathrm{\θ}\left(k\right)\]\stackrel{\‾}{y}\left(k\right)\end{array}\right.$

式中，y(k)＝[y₁(k) y₂(k)]^T，y₁(k)为物质A在k时刻浓度的测量值，y₂(k)为在k时刻反应釜内温度的测量值，v(k)为在k时刻反应釜内的有界扰动，即||v(k)||≤0.2， 为k时刻系统状态x(k)的预测值，C为测量输出权重矩阵，为单位矩阵，θ₁(k)＝1表示物质A在k时刻有测量值，θ₁(k)＝0表示物质A在k时刻没有测量值，θ₂(k)＝1；

(4)、设定滚动时域窗口长度N和权重矩阵，将多速率CSTR的滚动时域估计问题转化为等价的最小化问题；

$\underset{\hat{x}(k-N)}{min}J\left(k\right)$

约束条件：

$J\left(k\right)=\left|\right|\hat{x}(k-N)-\stackrel{\‾}{x}(k-N)|{|}_{\mathrm{\μ}}^2+\underset{i=k-N}{\overset{k}{\Σ}}||y\left(i\right)-\mathrm{\θ}\left(i\right)\[C\hat{x}\left(i\right)+v\left(i\right)\]+\[I-\mathrm{\θ}\left(i\right)\]\stackrel{\‾}{y}\left(i\right)|{|}^2$

$\hat{x}(i+1)=A\hat{x}\left(i\right),i=k-N,...,k-1$

$\stackrel{\‾}{x}(k-N)=A\hat{x}(k-N-1),k=N+1,N+2,...$

式中，J(k)为性能指标，为k时刻x(k)的最优估计值；

(5)、通过一阶最优性原理求解步骤(4)设定的最小化问题，具体步骤如下：

S1-1：初始化，设定测试时间长度K，在可行域的区间范围内，任意初始化k时刻的先验估计值，k-N时刻到k时刻的测量输出序列；

S1-2：根据一阶最优性原理，以k时刻的先验估计值为初始迭代点，对步骤(4)中的最小化问题求一阶偏导数，得到k-N时刻最优估计值；

S1-3：根据滚动优化原理，计算当前k时刻的最优估计值；

S1-4：根据k时刻的最优估计值更新k+1时刻的先验估计值；

S1-5：判断终止条件：如果k＝K，结束，得到信噪比估计最优值；否则，k＝k+1，转到S1-2。