CN110957011B - 连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法 - Google Patents

Abstract

连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法，属于复杂过程状态估计与监测领域。首先建立反应器的非线性动态模型，用一组加权粒子来表示系统状态的概率密度，设定测量噪声协方差矩阵为对角阵，同时用逆伽马分布来分析其概率密度。具体步骤包括：先对逆伽马分布参数和系统状态进行预测，再对粒子权重以及逆伽马分布参数进行更新，然后经过重采样得到一组新的粒子及权重，最后得到过程状态估计值以及测量噪声协方差矩阵的估计值。本方法能够同时估计出系统测量噪声协方差矩阵以及反应器内的反应温度和产品浓度。测量噪声协方差矩阵的估计值使得系统状态的估计值更加精确，为反应过程的安全进行以及产品的生产质量提供有力保障。

Description

连续搅拌反应器在未知时变测量噪声下的在线生产参数估计 方法

技术领域

本发明属于复杂过程状态估计与监测领域，涉及带有未知时变测量噪声的连续搅拌反应器生产参数与测量噪声协方差联合估计方法。

背景技术

连续搅拌反应器(CSTR)是一种比较典型且复杂的非线性过程，广泛用于化学生产和生物制药中，并且是进行各种复杂化学反应的重要设备，其操作状况的好坏对于产品的生产效率、品质和产能等会产生直接的影响。对于CSTR系统，反应温度和产品浓度是两个重要参数。反应温度通常对产物的数量和质量具有决定性的影响，在反应过程中，如果不能很好的控制温度，不仅会降低生产效益，更可能影响反应的安全性。产品浓度直接显示反应器中的反应状态和产物质量。然而，实际应用中浓度的实时检测具有很高的复杂性，成本也比较昂贵，存在检测条件、资金等多方面的限制。实时监测反应温度和产品浓度是一项重要任务。利用状态估计方法可以获得生产参数，但现有的状态估计方法都是在测量噪声已知的情况下实现对反应温度和产品浓度的估计；在测量噪声未知的情况下，目前没有良好的估计方法能够对反应器内的反应温度及产品浓度进行实时估计。

发明内容

针对上述现有技术中存在的问题，本发明针对在系统测量噪声未知且时变的情况下，提出一种对连续搅拌反应器在线过程状态估计方法。在系统测量噪声协方差矩阵未知的情况下，同时估计出系统状态即反应器内的反应温度和产品浓度以及测量噪声协方差矩阵，保证反应的安全性以及生产质量。

本发明的技术方案：

连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法，步骤如下：

第一步：建立连续搅拌反应器的非线性模型：

其中，

表示产品浓度随时间的变化，

表示反应温度随时间的变化，C_A为产品浓度，T为反应温度，F为进料流量，V为反应器体积，C_Af为进料浓度，k₀为反应速率常数，E为活化能项，R为摩尔气体常数，T_f为进料温度，hA为热转移项，ρ为产品密度，C_p为产品热容，T_c为冷却剂入口温度，λ为反应热。

选择反应器状态为x＝[x₁,x₂]^T＝[C_A,T]^T，因此式(1)和(2)写为：

建立输出方程为：

y＝x+ν (5)

其中，x₁为产品浓度，x₂为反应温度，ω₁和ω₂均为过程噪声，y为测量值,ν为测量噪声；

将上述(3)、(4)、(5)用四阶龙格库塔法离散化，得到以下形式的离散非线性状态和测量方程：

x_k＝f_k(x_k-1,u_k-1)+ω_k-1 (6)

y_k＝g_k(x_k)+ν_k (7)

其中，x_k为系统状态即产品浓度和反应温度，f_k为系统状态映射方程，u_k为系统输入，y_k为产品浓度和反应温度的测量值，g_k为系统测量方程，ω_k为系统过程噪声且ω_k～N(0,Q_k)，ν_k为系统测量噪声且ν_k～N(0,R_k)，Q_k已知，R_k未知；通过使用变分贝叶斯理论，将系统状态和测量噪声的协方差的联合后验分布用两个独立的分布来表示。由于系统本身存在的非线性，估计中涉及的状态的概率密度函数在计算上是难以处理的。因此，生成了一组加权粒子以表示系统状态的概率密度，即

其中，y_1:k＝{y₁,...,y_k}为测量值序列，N为粒子数，

为粒子权重，δ(·)为狄拉克δ函数，k为时间索引，

为k时刻N个粒子序列；

设定测量噪声协方差矩阵为

其中，

为R_k的对角线元素，d_y为y_k的维度。使用逆伽马分布来分析测量噪声协方差的概率密度，即

其中，α和β分别为逆伽马分布的形状参数和尺度参数。对α和β采用自启发式的动态模型，即

其中，ζ∈(0,1]为自启发式动态模型的系数；

第二步：设定初始值

α₀，β₀，Q_k，ζ，steps，M，dt，T。其中，steps为总的采样次数，M为每一时刻迭代的次数，dt为采样时间，T为反应时间；

令k＝1；

第三步：对系统状态以及逆伽马分布参数进行预测，状态的预测计算公式为：

其中，

为根据过程噪声的分布产生的噪声采样粒子，逆伽马分布参数的预测计算公式为式(10)和(11)；

第四步：更新逆伽马分布的形状参数α，更新公式为：

令m＝1；

第五步：计算测量噪声协方差矩阵的逆的期望

计算公式为：

第六步：计算粒子权重

计算公式为：

第七步：更新逆伽马的尺度参数β，更新公式为：

判断是否满足m＝M，是则执行下一步；否则m＝m+1，并跳转至第五步；

第八步：输出k时刻粒子权重值以及逆伽马分布的尺度参数的更新值，即

第九步：进行随机重采样得到一组新的粒子及权重：

第十步：输出k时刻的状态估计值以及测量噪声协方差矩阵的估计值：

第十一步：判断是否满足k＝steps，是则结束；否则k＝k+1，并跳转至第三步。

其中，steps＝T/dt。

本发明的有益效果：本发明通过建立CSTR的非线性系统模型，在未知时变测量噪声的情况下，同时估计出系统测量噪声协方差矩阵以及反应器内的反应温度和产品浓度。测量噪声协方差矩阵的估计值使得系统状态即反应器内的反应温度和产品浓度的估计值更加精确，为反应的安全进行以及产品的生产质量提供有力保障。

附图说明

图1为连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法的流程图。

图2为反应器生产参数估计误差图，其中：(a)为产品浓度即x₁，(b)为反应温度即x₂。

图3为测量噪声协方差矩阵对角线元素的估计效果图，其中：(a)为第一个对角线元素即

(b)为第二个对角线元素即

(a)和(b)中，平滑曲线反映的是真实值，波动曲线反映的是估计值。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步说明。

参照图1，连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法为：

根据式(1)和(2)建立反应器的非线性系统模型；

选择生产参数中产品浓度和反应温度为系统状态，建立(3)和(4)所示的状态方程以及(5)式所示的输出方程，将连续非线性系统离散化得到(6)式和(7)式所示的离散非线性状态方程和测量方程；

用一组加权粒子来表示系统状态的概率密度，设定测量噪声协方差矩阵为对角阵，同时用逆伽马分布来分析其概率密度；

根据式(10)、(11)和(12)对逆伽马分布参数和系统状态进行预测；

根据式(13)对逆伽马分布的形状参数进行更新；

在k时刻迭代M次后获得k时刻的粒子权重值以及逆伽马分布的尺度参数的更新值；

通过随机重采样得到一组新的粒子及权重，如式(19)和(20)所示；

最后，根据式(21)和(22)计算得到第k时刻的状态估计值和测量噪声协方差矩阵的估计值，用于下一时刻计算。

实施效果分析：

采用本发明提出的连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法，生产参数的初始状态选为x₀＝[1.167×10^-4,331]^T，反应器运行过程中的噪声方差矩阵为Q_k＝diag(5.43×10^-5,0.36)，采样时间为dt＝0.001s，粒子数为N＝100，逆伽马分布参数的初始值为α₀＝[3,3]和β₀＝[0.02,2]，自启发式动态模型的系数ζ＝0.99，每一时刻迭代的次数为M＝3。采用表1中的控制参数，同时为了可比性，采用测量噪声已知情况下的普通粒子滤波方法和本发明方法进行比较。搅拌反应器的生产参数估计误差图以及测量噪声协方差矩阵对角线元素的估计效果图分别如图2和图3所示。

从图2中可以看出，使用本发明方法(图中标注为VB-PF)进行状态估计的误差比普通粒子滤波方法(图中标注为PF)的波动小。从图3中可以看出系统测量噪声协方差矩阵的对角线元素可以较好的跟踪上真实值即波动曲线。系统状态在两种估计方法下的均方根误差如表2所示。从表2可以看出本发明方法的均方根误差比普通粒子滤波方法的均方根误差小。因此，本发明所提出的连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法具有很好的抗未知时变噪声协方差的性能，能够为反应的安全进行以及产品的生产质量提供有力保障。

表1控制参数表

表2 CSTR系统状态估计误差的均方根误差

Claims (1)

1.连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法，其特征在于，步骤包括：

第一步：建立CSTR的非线性模型：

其中，

表示产品浓度随时间的变化，

表示反应温度随时间的变化，C_A为产品浓度，T为反应温度，F为进料流量，V为反应器体积，C_Af为进料浓度，k₀为反应速率常数，E为活化能项，R为摩尔气体常数，T_f为进料温度，hA为热转移项，ρ为产品密度，C_p为产品热容，T_c为冷却剂入口温度，λ为反应热；

选择反应器状态为x＝[x₁,x₂]^T＝[C_A,T]^T，因此式(1)和(2)写为：

建立输出方程为：

y＝x+ν (5)

其中，x₁为产品浓度，x₂为反应温度，ω₁和ω₂均为过程噪声，y为测量值，ν为测量噪声；

将上述(3)、(4)、(5)用四阶龙格库塔法离散化，得到以下形式的离散非线性状态和测量方程：

x_k＝f_k(x_k-1,u_k-1)+ω_k-1 (6)

y_k＝g_k(x_k)+ν_k (7)

其中，x_k为系统状态即产品浓度和反应温度，f_k为系统状态映射方程，u_k为系统输入，y_k为产品浓度和反应温度的测量值，g_k为系统测量方程，ω_k为系统过程噪声且ω_k～N(0,Q_k)，ν_k为系统测量噪声且ν_k～N(0,R_k)，Q_k已知，R_k未知；通过使用变分贝叶斯理论，将系统状态和测量噪声的协方差的联合后验分布用两个独立的分布来表示；生成了一组加权粒子以表示系统状态的概率密度，即

其中，y_1:k＝{y₁,…,y_k}为测量值序列，N为粒子数，

为粒子权重，δ(·)为狄拉克δ函数，k为时间索引，

为k时刻N个粒子序列；

设定测量噪声协方差矩阵为

其中，

为R_k的对角线元素，d_y为y_k的维度；使用逆伽马分布来分析测量噪声协方差的概率密度，即

其中，α和β分别为逆伽马分布的形状参数和尺度参数；对α和β采用自启发式的动态模型，即

其中，ζ∈(0,1]为自启发式动态模型的系数；

第二步：设定初始值

α₀，β₀，Q_k，ζ，steps，M，dt，T；其中，steps为总的采样次数，M为每一时刻迭代的次数，dt为采样时间，T为反应时间；

令k＝1；

第三步：对系统状态以及逆伽马分布参数进行预测，状态的预测计算公式为：

其中，

为根据过程噪声的分布产生的噪声采样粒子，逆伽马分布参数的预测计算公式为式(10)和(11)；

第四步：更新逆伽马分布的形状参数α，更新公式为：

令m＝1；

第五步：计算测量噪声协方差矩阵的逆的期望

计算公式为：

第六步：计算粒子权重

计算公式为：

第七步：更新逆伽马的尺度参数β，更新公式为：

判断是否满足m＝M，是则执行下一步；否则m＝m+1，并跳转至第五步；

第八步：输出k时刻粒子权重值以及逆伽马分布的尺度参数的更新值，即

第九步：进行随机重采样得到一组新的粒子及权重：

第十步：输出k时刻的状态估计值以及测量噪声协方差矩阵的估计值：

第十一步：判断是否满足k＝steps，是则结束；否则k＝k+1，并跳转至第三步；其中，steps＝T/dt。

Images