CN111145838B - 基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法 - Google Patents

Abstract

基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，属于间歇生产过程状态估计领域。首先针对青霉素发酵过程选择变量；然后采集生产过程的多批次历史数据建立多向线性变参数模型，并用期望最大化算法估计模型参数，再扩展成含噪声干扰的状态空间形式的多向线性变参数模型；接着基于SS‑MLPV模型构造一个适用于迭代学习卡尔曼滤波的误差模型，将误差模型分解为时间方向子模型与批次方向子模型；最后基于当前间歇发酵过程的青霉素浓度数据，第一批使用卡尔曼滤波对青霉素浓度的状态估计，第二批开始对构建的时间和批次方向的子模型进行ILKF得到误差估计，将不含噪声的SS‑MLPV模型输出与误差估计相加即为青霉素浓度的估计。

Description

基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波 方法

技术领域

本发明属于间歇生产过程状态估计领域，涉及一种青霉素发酵过程的关键变量检测的方法。

背景技术

间歇过程作为一种重要的生产方式，具有生产周期短、转型快等特点，能够满足快速化、个性化的市场需求，在精细化工、生物制药、食品和化妆品等行业中得到了广泛的应用。青霉素发酵过程是一种间歇过程，在一个或者多个发酵罐中生产产品青霉素，一个批次产品完成后，如果需要更多的产品，必须重复该过程，按照相同的工序进行下一个批次生产。在青霉素发酵生产中，存在不可测量或测量成本较高的过程变量(如发酵过程中的产物浓度，生物质浓度，基质浓度等)，这些过程变量也称为状态变量，某些关键的过程变量关系到产品质量，需要对其进行实时监测，因此对发酵过程状态变量的估计尤为重要，发酵过程的状态估计一直是工业界和学术界的关注焦点。

用于青霉素发酵过程的状态估计方法主要有：扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)、粒子滤波(PF)等，这些传统方法都需要得到系统的动态方程，并且只是在时间方向上进行估计，然而间歇过程通常具有许多重复的批次。当前工业生产中，在间歇过程中使用的状态估计方法大多是对连续过程的简单复制，忽略了批次之间的相关性。考虑到间歇过程多批次重复性质，迭代学习卡尔曼滤波(ILKF)被引入到间歇过程的状态估计中，在状态估计时同时考虑了批次方向的相关性和时间方向动态特性。但是这些方法都需要预先获得发酵过程的机理模型，当生产过程机理不明时，上述方法无法实际应用；如何根据过程运行数据和操作经验来建立模型，进而基于该模型实现对关键变量检测尚无更好解决方案。

发明内容

在青霉素发酵过程中，一些关键参数检测很困难，现有的软测量技术可以对其进行估计，但是一般基于静态数据模型。状态空间模型是一种动态模型，考虑了变量在时间方向上的联系，适用于线性、时变、非线性、随机、采样等各种各样的系统；并且能简化过程描述，方便推导，适用于计算机处理。首先，在过程机理不明时，本发明基于数据构造一个包含关键参数的状态空间模型，对所构造模型的状态进行估计即实现对青霉素发酵过程中的关键参数的估计；然后，本发明将给出一种基于数据的青霉素发酵过程迭代学习状态估计的方案。

一方面，状态估计需要有一个比较准确的模型，在机理不明的情况下，考虑基于数据建模的方法。对于多阶段的非线性间歇青霉素发酵过程，很难用一个非线性全局模型描述，线性变参数(LPV)模型可以逼近任意非线性模型，是一个较好的选择。但针对间歇过程这种多批次生产过程，需要一个可以拟合每一批数据的模型，本发明采用多向线性变参数(MLPV)模型，能较快地得到复杂非线性过程的模型，且此模型又有较好的拟合能力。并且将MLPV模型扩展到状态空间(SS)形式，考虑到实际过程存在的系统噪声，测量噪声和重复干扰，构造了比较符合实际生产的SS-MLPV模型。

另一方面，青霉素发酵过程是个多批次的重复过程，适合用迭代学习状态估计的方法。然而现有的迭代学习估计方法需要过程机理模型，在只有数据的情况下，不适用于实际间歇生产过程。本发明通过构造误差模型，给出了适合实际生产过程的ILKF方案，与原ILKF方法相比，该方案考虑了输入条件，并放宽了初值假设。

本发明能够充分利用多批次数据对机理不明的间歇过程进行建模，得到拟合度较高的模型，再改进现有的迭代学习卡尔曼滤波方法，在模型含有一定干扰和噪声的情况下，状态估计的精度随着批次的进行能有提高。本发明是基于MLPV模型的迭代学习卡尔曼滤波方法，因此技术方案分为两个部分：基于数据建立MLPV模型；基于模型用迭代学习卡尔曼滤波方法进行状态估计。首先根据输入、输出和调度变量数据，采用MLPV方法进行建模，然后考虑建模误差以及系统噪声，将原未知系统模型描述为一个含随机噪声以及重复干扰的状态空间形式的MLPV模型。在MLPV模型的基础上，再将系统改写成为时间方向和批次方向的子系统，并设计相应的卡尔曼滤波器，得到当前批次当前时刻的状态估计值。

本发明的技术方案：

基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，步骤如下：

步骤1，根据青霉素发酵过程历史运行数据和操作经验，选择MLPV模型的输入条件u、待估计状态s、调度变量z、子模型个数M，二维K均值聚类算法对调度变量进行聚类得到聚类中心点，即视为每个子模型的工况点π_1:M。选择子模型结构为受控自回归模型(ARX)，其中s和u的阶次分别记为n_a和n_b。

步骤2，基于步骤1的选择采集历史数据，进行模型辨识。时间、批次和变量的数据构成一个三维矩阵，将多批历史数据沿时间轴展开成二维时间片扩展矩阵。用所有批次的数据建立一个ARX融合模型(包含M个子模型的融合模型)可以表述为：其中s(t)是第t个采样时间的状态，α_tm(z_t)是与调度变量z_t有关的归一化后的权重函数，/>是第t个采样时间的信息向量，需要辨识的主要参数是θ_m，用EM算法进行MLPV模型参数辨识，得到θ_m。

步骤3，基于步骤2的融合模型，建立考虑高斯过程干扰ω(t)以及重复过程干扰d(t)的SS-MLPV模型结构如下：

x(t+1)＝A(z_t)x(t)+Bu(t)+Γω(t)+Γd(t)

其中：扩展状态矩阵x(t)＝[s(t),...,s(t-n_a+1),u(t),...,u(t-n_b+1)]^T，s(t)表示待估计状态，u(t)表示输入，ω(t)表示过程干扰，d(t)表示重复过程干扰。

用O^a×b表示a行b列的0阵，用I^a表示a行a列的单位阵，则A(z_t),B,Γ分别为：

Γ＝[1 O^1×(n-1)]^T,

其中：

用t表示时间，k表示批次，构造考虑观测干扰的SS-MLPV模型如下：

x_k(t+1)＝A(t)x_k(t)+Bu_k(t)+Γω_k(t)+Γd(t)

y_k(t)＝Cx_k(t)+v_k(t)

其中C＝Γ^T。x_k(t)，y_k(t)，ω_k(t)和v_k(t)分别表示第t个采样时间第k批的扩展状态，输出观测值，过程干扰和观测干扰。

构造不考虑噪声干扰的标准SS-MLPV模型如下：

并且构造误差系统模型如下：

ex_k(t+1)＝A(t)ex_k(t)+Γω_k(t)+Γd(t)

ey_k(t)＝Cex_k(t)+v_k(t)

步骤4，基于步骤3的误差系统模型设计时间和批次两个子模型结构如下：

时间子系统：批次子系统：其中δ是批次间的误差算子，τ是一个接近于1的数，/>为引入的批次间干扰项。

步骤5，青霉素发酵开始时，根据第一批生产过程的输出观测值y，与初始扩展状态向量x，基于步骤3的SS-MLPV模型用经典的卡尔曼滤波(KF)方法得到第一批的估计。

步骤6，第二批生产开始使用迭代学习卡尔曼滤波(ILKF)方法。基于步骤(4)的时间和批次两个子模型、当前批次生产过程的输出观测值y_k(t)、上一批次的输出观测值y_k-1(t)以及上一批次的估计值，分别对时间子系统与批次子系统进行状态估计。具体地，首先对时间维度子系统进行卡尔曼滤波，将其结果带入批次子系统再进行一次卡尔曼滤波。迭代学习卡尔曼滤波器结构如下：

∑_T:

∑_B:其中/>分别为时间子系统与批次子系统在第t个采样时间第k批的卡尔曼增益，可通过最小化估计误差协方差计算，δ是批次间的误差算子，τ是一个接近于1的数，得到误差的估计：/>与/>计算当前关于青霉素浓度的状态估计为/>

本发明的有益效果：实际青霉素发酵生产中，往往需要估计一些难以测量的状态变量。一般的状态估计方法需要有一个精确地描述整个过程的机理模型，但是有一些复杂过程难以分析出机理，或者根据机理建立的模型与实际模型间存在较大偏差，用传统的状态估计方法进行估计，估计误差会很大。本发明给出了基于批次操作的多向数据的迭代学习状态估计方案。只需要利用过程中现有的历史数据，即可进行比较精准的建模，并在考虑模型偏差的基础上，进行迭代学习状态估计，使估计值比使用传统状态估计方法更精确。相比于一般的LPV模型，本发明方法采用的是MLPV模型，考虑了多批数据融合问题，更符合实际工业间歇过程的多批次数据特点，提供了一个能拟合多批数据的模型。另一方面，相比于基于传统状态估计方法，本发明得到的状态估计值会随着批次的进行越来越准确，且初始条件的自由度更大，更具有灵活性。

附图说明

图1为青霉素发酵过程MLPV建模过程图。

图2为青霉素发酵过程建模的流程图。

图3为青霉素发酵过程状态估计流程图。

图4为青霉素浓度估计轨迹图。

图5为青霉素浓度的两种估计方法的均方根误差随批次变化图。

具体实施方式

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

参照图1、图2、图3，一种基于数据的青霉素发酵过程的迭代学习卡尔曼滤波方法为：

根据步骤1，选择需要估计的状态为青霉素浓度s，选择冷水流量为输入条件u与调度变量z，根据青霉素发酵特性将子模型个数M设定为3，针对多批历史数据两次K均值算法对调度变量进行聚类得到聚类中心点，即视为每个子模型的工况点π_1:M。选择子模型结构为ARX模型，其中s和u的阶次分别记为n_a＝2以及n_b＝2。

根据步骤2，采集多批次历史数据，用EM算法辨识MLPV模型参数。

根据步骤3，将辨识出的青霉素发酵模型转化为SS-MLPV模型，进一步得到真实值与辨识模型之间的误差模型。

根据步骤4，考虑多批次，将得到的误差模型分解为时间方向子模型与批次方向子模型。

根据步骤5，对于新的青霉素发酵过程，基于测量值与步骤3模型采用KF方法的到第一批次的青霉素浓度的估计值。

根据步骤6，第二批开始青霉素浓度估计采用ILKF方法，利用上一批次和当前批次的信息以及步骤4中的模型分两次估计后得到误差估计值。误差的估计与不含噪声的辨识模型估计值相加得到更加准确的估计值。

实施例：

青霉素是生活中最常见和最常用的一种抗生素，被广泛应用于实际生产和生活中。青霉素发酵过程是典型的间歇过程，通过青霉素发酵过程的仿真应用来验证本发明。青霉素仿真软件Pensim已经广泛应用于各种过程监控、控制和优化等算法的验证中，其初始条件如表1所示，操作变量的设定值如表2所示。

发酵反应温度是影响产物浓度至关重要的一个变量，由于发酵罐中温度难以保证处处一致且难以测量，温度主要是由冷却装置的冷却剂来调节。因此，在给定进料流量，浓度以及装置的体积时，整个反应过程随着冷却剂流量变化而变化。所以可选择选择冷水流量作为输入和调度变量，青霉素浓度作为输出。采用10个批次的数据建立MLPV模型，子模型个数为3个，用二维K均值算法对调度变量进行聚类得到每个子模型的工况点，ARX子模型中s和u的阶次都设定为2。图4展示了本发明最终估计几乎以原轨迹重合，图5表明本发明涉及的ILKF状态估计方法比传统状态估计方法KF要好，随着批次的进行，估计精度提高。

表1初始条件

初始条件	默认值	范围
			基质浓度(g/L)	15	5-50
菌体浓度(g/L)	0.1	0-0.2
			溶解氧浓度(mmol/L)	1.16	1-1.2
二氧化碳浓度(mmol/L)	0.5	0.5-1
			培养基浓度(vol./L)	100	100-200
pH	5	4-6
			温度(K)	298	298-300

表2操作变量

操作变量	默认值	范围
			基质流加速率(L/h)	0.042	0.035-0.045
通风速率(L/h)	8.6	3-10
			搅拌功率(W)	30	20-50

Claims (7)

1.基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，其特征在于，包括以下步骤：

(1)根据青霉素发酵过程历史运行数据和操作经验，选择多向线性变参数模型MLPV的输入条件、需要估计的状态、调度变量、子模型个数、寻找子模型的工况点，选择子模型结构；

(2)基于步骤(1)的选择采集历史数据，进行模型参数辨识；构成时间、批次和变量三维数据矩阵，再将多批历史数据沿时间轴展开成二维时间片扩展矩阵；用所有批次的数据建立一个融合模型，融合模型包含M个子模型的融合模型，用EM算法进行MLPV模型参数辨识，得到模型参数；

(3)基于步骤(2)的融合模型，建立考虑高斯过程干扰ω(t)以及重复过程干扰d(t)的状态空间形式的多向线性变参数系统模型SS-MLPV，并且构造一个误差系统模型；所述步骤(3)模型结构如下：

x(t+1)＝A(z_t)x(t)+Bu(t)+Γω(t)+Γd(t)

其中：扩展状态矩阵x(t)＝[s(t),...,s(t-n_a+1),u(t),...,u(t-n_b+1)]^T，s(t)表示待估计状态，u(t)表示输入，ω(t)表示过程干扰，d(t)表示重复过程干扰；

用O^a×b表示a行b列的0阵，用I^a表示a行a列的单位阵，则A(z_t),B,Γ分别为：

Γ＝[1 O^1×(n-1)]^T,

其中：

用t表示时间，k表示批次，构造考虑观测干扰的SS-MLPV模型如下：

x_k(t+1)＝A(t)x_k(t)+Bu_k(t)+Γω_k(t)+Γd(t)

y_k(t)＝Cx_k(t)+v_k(t)；

其中，C＝Γ^T；x_k(t)，y_k(t)，ω_k(t)和v_k(t)分别表示第t个采样时间第k批的扩展状态，输出观测值，过程干扰和观测干扰；

不考虑噪声干扰的SS-MLPV模型如下：

误差模型如下：

ex_k(t+1)＝A(t)ex_k(t)+Γω_k(t)+Γd(t)

ey_k(t)＝Cex_k(t)+v_k(t)

A(t)根据调度变量与步骤(2)中辨识出的模型计算得到；

(4)基于步骤(3)的误差系统模型设计时间和批次两个子模型；

(5)青霉素发酵开始时，根据第一批生产过程的输出观测值y，与初始扩展状态向量x，基于步骤(3)的SS-MLPV模型用经典的卡尔曼滤波方法KF得到第一批的估计；

(6)第二批生产开始使用迭代学习卡尔曼滤波方法ILKF；基于步骤(4)的时间和批次两个子模型，分别对时间子系统与批次子系统进行状态估计，得到误差的估计并计算当前的青霉素浓度的状态估计。

2.根据权利要求书1所述的基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，其特征在于，所述步骤(1)中选择MLPV模型的输入条件u、需要估计的状态s，调度变量z，子模型个数M，两次K均值算法对调度变量进行聚类得到聚类中心点，即视为每个子模型的工况点；选择子模型结构为受控自回归模型ARX，其中s和u的阶次分别记为n_a和n_b。

3.根据权利要求书1或2所述的基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，其特征在于，所述步骤(2)将多批历史数据沿时间轴展开成二维时间片扩展矩阵；用所有批次的数据建立一个ARX融合模型表述为：其中s(t)是第t个采样时间的状态，α_tm(z_t)是与调度变量z_t有关的归一化后的权重函数，/>是第t个采样时间的信息向量，需要辨识的参数是θ_m，用EM算法进行MLPV模型参数辨识，得到θ_m。

4.根据权利要求书1或2所述的基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，其特征在于，所述步骤(4)基于步骤(3)的误差系统模型设计时间和批次两个子模型结构如下：

时间子系统：批次子系统：

其中δ是批次间的误差算子，τ是一个接近于1的数，为引入的批次间干扰项。

5.根据权利要求书3所述的基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，其特征在于，所述步骤(4)基于步骤(3)的误差系统模型设计时间和批次两个子模型结构如下：

时间子系统：批次子系统：

其中δ是批次间的误差算子，τ是一个接近于1的数，为引入的批次间干扰项。

6.根据权利要求书1、2或5所述的基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，其特征在于，所述步骤(6)中，第二批生产开始使用ILKF方法，基于步骤(4)的时间和批次两个子模型，以及当前批次生产过程的输出观测值、上一批次的输出观测值、上一批次的估计值，分别对时间子系统与批次子系统进行状态估计，具体为：首先对时间维度子系统进行卡尔曼滤波，将其结果带入批次子系统再进行一次卡尔曼滤波；迭代学习卡尔曼滤波器结构如下：

其中分别为时间子系统与批次子系统在第t个采样时间第k批的卡尔曼增益，通过最小化估计误差协方差计算，δ是批次间的误差算子，τ是一个接近于1的数，得到误差的估计：/>与/>计算当前的青霉素浓度的状态估计为/>

7.根据权利要求书4所述的基于多向数据模型的青霉素发酵过程迭代学习卡尔曼滤波方法，其特征在于，所述步骤(6)中，第二批生产开始使用ILKF方法，基于步骤(4)的时间和批次两个子模型，以及当前批次生产过程的输出观测值、上一批次的输出观测值、上一批次的估计值，分别对时间子系统与批次子系统进行状态估计，具体为：首先对时间维度子系统进行卡尔曼滤波，将其结果带入批次子系统再进行一次卡尔曼滤波；迭代学习卡尔曼滤波器结构如下：

其中分别为时间子系统与批次子系统在第t个采样时间第k批的卡尔曼增益，通过最小化估计误差协方差计算，δ是批次间的误差算子，τ是一个接近于1的数，得到误差的估计：/>与/>计算当前的青霉素浓度的状态估计为/>