CN113326618B - 连续发酵过程中估计培养基初始条件的方法 - Google Patents

连续发酵过程中估计培养基初始条件的方法 Download PDF

Info

Publication number
CN113326618B
CN113326618B CN202110615110.9A CN202110615110A CN113326618B CN 113326618 B CN113326618 B CN 113326618B CN 202110615110 A CN202110615110 A CN 202110615110A CN 113326618 B CN113326618 B CN 113326618B
Authority
CN
China
Prior art keywords
time
state
variable
fermentation process
continuous fermentation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110615110.9A
Other languages
English (en)
Other versions
CN113326618A (zh
Inventor
赵顺毅
张田雨
郭松杰
栾小丽
刘飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jiangnan University
Original Assignee
Jiangnan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jiangnan University filed Critical Jiangnan University
Priority to CN202110615110.9A priority Critical patent/CN113326618B/zh
Publication of CN113326618A publication Critical patent/CN113326618A/zh
Application granted granted Critical
Publication of CN113326618B publication Critical patent/CN113326618B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Abstract

本发明公开了一种连续发酵过程中估计培养基初始条件的方法，包括建立连续发酵过程的三变量数学模型、线性化的连续发酵过程的状态方程和系统测量方程，将连续的方程离散化得到离散状态空间模型；在卡尔曼滤波的基础上引入辅助变量得到服从学生分布的预测分布；使用变分贝叶斯理论表示系统状态和辅助变量；预测系统状态和辅助变量，更新系统状态和辅助变量直到达到设置的迭代次数，输出此时的系统状态和辅助变量；直到达到设置的运行步数结束更新，得到初始时刻培养基中的关键变量的估计值。本发明通过建立连续发酵过程的三变量数学模型、引入学生分布，快速且准确地估计出培养基中关键变量初始条件，保证发酵过程安全和发酵质量。

Description

连续发酵过程中估计培养基初始条件的方法
技术领域
本发明涉及发酵培养基初始条件估计技术领域，具体涉及一种连续发酵过程中估计培养基初始条件的方法。
背景技术
在现有的连续发酵工艺中，由于传感器设备的不完善或者工艺条件的限制，操作人员大多是采用离线分析的方法来分析培养基中关键变量的数值。但是离线分析的方法是间歇性的，在间歇性过程中进一步增大了时延特性，不利于发酵质量的实时监测或者发酵过程中故障的实时诊断。培养基的初始条件估计是分析培养基的关键一环，培养基初始条件估计不准确会影响整个发酵过程的质量。虽然目前有培养基初始化条件的估计方法，但存在计算量大、性能指标不确定等问题，使得这些估计方法在连续发酵过程中的应用受到限制。特别是针对连续发酵过程中的关键变量，目前没有良好的估计方法来估计其初始条件。
发明内容
为此，本发明所要解决的技术问题在于克服现有技术中的不足，提出一种快速且准确地估计连续发酵过程中的培养基初始化条件的连续发酵过程中估计培养基初始条件的方法。
为解决上述技术问题，本发明提供了一种连续发酵过程中估计培养基初始条件的方法，包括以下步骤：
步骤1：建立连续发酵过程的三变量数学模型和连续发酵过程的状态方程，利用泰勒级数展开连续发酵过程的状态方程得到线性化的连续发酵过程的状态方程，建立系统测量方程；
步骤2：将线性化的连续发酵过程的状态方程和系统测量方程所描述的连续状态空间表达式离散化，得到离散状态空间模型；
步骤3：在卡尔曼滤波的基础上引入辅助变量，得到服从学生分布的预测分布；
步骤4：使用变分贝叶斯理论将每一时刻下的系统状态和辅助变量的联合后验概率密度函数用两个独立的概率密度函数表示；
步骤5：设定系统参数初始值，总运行步数steps和每一步迭代的总次数N；
初始化离散状态空间的时间索引n＝1；
步骤6：预测n时刻的状态预测均值
Figure GDA0003630065230000021
和伽马分布的n时刻的辅助变量的预测形状参数
Figure GDA0003630065230000022
n时刻的辅助变量的预测逆尺度参数
Figure GDA0003630065230000023
初始化迭代次数j＝1；
步骤7：更新n时刻的辅助变量的形状参数an、第j次迭代得到的n时刻的状态变量的预测协方差
Figure GDA0003630065230000024
第j次迭代得到的n时刻的状态变量的均值
Figure GDA0003630065230000025
第j次迭代得到的n时刻的状态变量的协方差
Figure GDA0003630065230000026
和第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure GDA0003630065230000027
步骤8：判断迭代次数j是否满足j＞N，若是则执行步骤9；若否则令j＝j+1，并跳转执行步骤7；
步骤9：输出当前时刻的状态变量的均值、状态变量的协方差、辅助变量的形状参数和辅助变量的逆尺度参数；
步骤10：判断时间索引n是否满足n＞steps，若否则令n＝n+1，并跳转执行步骤6；若是则结束更新，得到初始时刻培养基中的关键变量的估计值。
进一步地，所述步骤1中的三变量数学模型，具体为：
系统状态
Figure GDA0003630065230000031
系统输入u＝sf；
其中c是培养基中微生物细胞浓度，s是培养基中基质浓度，p是培养基中产物浓度，sf是培养基中补给基质浓度；
所述连续发酵过程的状态方程为：
Figure GDA0003630065230000032
其中t表示连续状态空间表达式的时间索引，
Figure GDA0003630065230000033
是t时的系统状态的变化量，xt是t时的系统当前状态，ut是t时的系统控制输入，ωt是t时的过程噪声、ωt服从均值为零的高斯分布即ωt～N(0，Qt)，Qt为t时的过程噪声协方差矩阵，t时的状态转移矩阵
Figure GDA0003630065230000034
Figure GDA0003630065230000035
是欧式空间，Bt是t时的控制输入矩阵；
At具体为：At(1，2)＝At(1，3)＝At(2，3)＝At(3，2)＝0，At(1，1)＝-D+μ,
Figure GDA0003630065230000036
At(2，2)＝-D,At(3，1)＝wμ+β，At(3，3)＝-D；
Bt具体为：Bt(2，1)＝-D,Bt(1，1)＝Bt(3，1)＝0；
其中D是稀释速率，μ是比生产速率，Yc/s是细胞浓度对基质浓度的生产速率，ω是产物浓度变化中与细胞浓度变化有关的参数，β是产物浓度变化中与细胞浓度相关的参数；
Figure GDA0003630065230000037
其中μm是最大比生产速率，Km是半饱和速度常数。
进一步地，所述步骤1中的线性化的连续发酵过程的状态方程，具体为：
Figure GDA0003630065230000038
其中A′t(1，3)＝A′t(2，3)＝0，A′t(1，)＝umθ-D，A′t(1，2)＝μmψ，
Figure GDA0003630065230000041
A′t(3，1)＝β+wμmθ，A′t(3，2)＝β+wμmψ，A′t(3，3)＝-D，
Figure GDA0003630065230000042
Ss是稳态情况下的基质浓度，cs是稳态情况下细胞浓度；
所述步骤1中建立的系统测量方程，具体为：
yt＝Ctxt+vt；
其中yt为t时的系统测量值，Ct为t时的系统系数矩阵，vt为t时的测量噪声且vt服从均值为零的高斯分布，Rt为t时的测量噪声协方差矩阵。
进一步地，所述步骤2中的离散状态空间模型，具体为：
xn＝Fn-1xn-1+Ln-1un-1+wn-1，
yn＝Hnxn+vn；
其中xn是n时刻的系统状态，Fn-1是n-1时刻的系统转移矩阵，xn-1是n-1时刻的系统状态，Ln-1是n-1时刻的控制矩阵，un-1是n-1时刻的系统控制输入，wn-1是n-1时刻的过程噪声，yn是n时刻的系统测量变量，Hn是n时刻的测量矩阵，vn是n时刻的测量噪声。
进一步地，所述步骤3中的服从学生分布的预测分布，具体为：
Figure GDA0003630065230000043
其中un是引入的n时刻的辅助变量，y1：n-1是1时刻到n-1时刻的所有测量值序列，p(xn|un，y1：n-1)服从于均值为
Figure GDA0003630065230000044
协方差为
Figure GDA0003630065230000045
的高斯分布，p(un)表示辅助变量un的先验分布、服从于形状参数为
Figure GDA0003630065230000051
逆尺度参数为
Figure GDA0003630065230000052
的伽马分布；
Figure GDA0003630065230000053
表示均值为
Figure GDA0003630065230000054
协方差为
Figure GDA0003630065230000055
自由度为
Figure GDA0003630065230000056
的学生分布，
Figure GDA0003630065230000057
是n时刻的状态变量的预测协方差，
Figure GDA0003630065230000058
是n时刻的状态变量的预测均值，
Figure GDA0003630065230000059
是n时刻的Student-t分布的协方差，
Figure GDA00036300652300000510
是自由度参数。
进一步地，所述步骤4中使用变分贝叶斯理论将每一时刻下的系统状态和辅助变量的联合后验概率密度函数用两个独立的概率密度函数表示，其中两个独立的概率密度函数分别服从于高斯分布和伽马分布，具体为：
p(xn，un|y1：n)≈q(xn，un|y1：n)＝q(xn|y1：n)q(un|y1：n)；
其中，p(xn，un|y1：n)是每一时刻下的系统状态和辅助变量的联合后验概率密度函数，q(·)是独立的概率密度函数，y1：n是1时刻到n时刻的所有测量值序列；
Figure GDA00036300652300000511
q(un|y1：n)＝G(un；an，bn)；
其中
Figure GDA00036300652300000512
是n时刻的状态变量的均值，Pn是n时刻的状态变量的协方差。
进一步地，所述步骤6中预测n时刻的状态预测均值
Figure GDA00036300652300000513
和伽马分布的n时刻的辅助变量的预测形状参数
Figure GDA00036300652300000514
n时刻的辅助变量的预测逆尺度参数
Figure GDA00036300652300000515
具体预测公式为：
n时刻的状态预测均值
Figure GDA00036300652300000516
的预测公式为：
Figure GDA00036300652300000517
其中xn-1是n-1时刻的状态变量的均值；
n时刻的辅助变量的预测形状参数的预测公式为：
Figure GDA00036300652300000518
n时刻的辅助变量的预测逆尺度参数的预测公式为：
Figure GDA00036300652300000519
其中ρ是启发式因子，取值范围(0，1]，an-1是n-1时刻的辅助变量的形状参数，bn-1是n-1时刻的辅助变量的逆尺度参数。
进一步地，所述步骤7中更新n时刻的辅助变量的形状参数an、第j次迭代得到的n时刻的状态变量的预测协方差
Figure GDA0003630065230000061
第j次迭代得到的n时刻的状态变量的均值
Figure GDA0003630065230000062
第j次迭代得到的n时刻的状态变量的协方差
Figure GDA0003630065230000063
和第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure GDA0003630065230000064
具体为：
n时刻的辅助变量的形状参数an的更新公式为：
Figure GDA0003630065230000065
其中d是状态变量的维度；
第j次迭代得到的n时刻的状态变量的预测协方差
Figure GDA0003630065230000066
的更新公式为：
Figure GDA0003630065230000067
其中
Figure GDA0003630065230000068
Pn-1是n-1时刻的状态变量的协方差，Fn-1是n-1时刻的系统转移矩阵，()T是矩阵的转置，Qn-1是n-1时刻的过程噪声协方差矩阵；
Figure GDA0003630065230000069
是第j次迭代得到的辅助变量un的期望且
Figure GDA00036300652300000610
第j次迭代得到的n时刻的状态变量的均值
Figure GDA00036300652300000611
的更新公式为：
Figure GDA00036300652300000612
第j次迭代得到的n时刻的状态变量的协方差
Figure GDA00036300652300000613
的更新公式为：
Figure GDA00036300652300000614
第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure GDA00036300652300000615
的更新公式为：
Figure GDA00036300652300000616
其中
Figure GDA00036300652300000617
是第j次迭代得到的n时刻的卡尔曼滤波增益且更新公式为
Figure GDA00036300652300000618
()-1是矩阵的逆，
Figure GDA00036300652300000619
的更新公式为
Figure GDA00036300652300000620
Rn是n时刻的测量噪声协方差矩阵；ψ更新公式为
Figure GDA0003630065230000071
tr[]为迹操作算子。
进一步地，所述步骤9中输出的当前时刻的状态变量的均值、状态变量的协方差、辅助变量的形状参数和辅助变量的逆尺度参数，具体为：
当前时刻的状态变量的均值为
Figure GDA0003630065230000072
当前时刻的状态变量的协方差为
Figure GDA0003630065230000073
当前时刻的辅助变量的形状参数为
Figure GDA0003630065230000074
当前时刻的辅助变量的逆尺度参数为
Figure GDA0003630065230000075
本发明还提供了一种使用初始时刻培养基中的关键变量估计值来监测连续发酵过程的方法，使用连续发酵过程中估计培养基初始条件的方法来获取初始时刻培养基中的关键变量估计值，并根据初始时刻培养基中的关键变量估计值判断是否处于正常工况或者发生故障。
本发明的上述技术方案相比现有技术具有以下优点：
本发明所述的连续发酵过程中估计培养基初始条件的方法通过建立连续发酵过程的典型三变量数学模型、引入学生分布估计培养基的初始条件，解决了卡尔曼滤波对初始值的依赖程度，实现培养基中关键变量初始条件的快速且准确地估计，保证了发酵过程的安全和发酵质量。
附图说明
为了使本发明的内容更容易被清楚的理解，下面根据本发明的具体实施例并结合附图，对本发明作进一步详细的说明。
图1是本发明的流程图。
图2是本发明中单级连续发酵过程的结构示意图。
图3是本发明和KF、PKF、BIA-II三种方法在连续发酵过程中细胞浓度的估计误差对比图。
图4是本发明和BIA-II在不同采样时间下的计算时间的对比结果图。
具体实施方式
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
在本发明的描述中，需要理解的是，术语“包括”意图在于覆盖不排他的包含，例如包含了一系列步骤或单元的过程、方法、系统、产品或设备，没有限定于已列出的步骤或单元而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其他步骤或单元。
参照图1流程图所示，本发明一种连续发酵过程中估计培养基初始条件的方法的实施例，包括以下步骤：
步骤1：建立连续发酵过程的典型三变量数学模型：
根据图2所示的连续发酵过程的结构示意图所示，微生物摄取营养液中的营养物质，支持自身新陈代谢，代谢物中为人利用的物质选择为目标产物，故三变量数学模型可具体描述为：
系统状态
Figure GDA0003630065230000081
系统输入u＝sf；
其中c是培养基中微生物细胞浓度，s是培养基中基质浓度、即营养液浓度，p是培养基中产物浓度，sf是培养基中补给基质浓度，
则连续发酵过程的状态方程可表达为：
Figure GDA0003630065230000082
其中t表示连续状态空间表达式的时间索引，
Figure GDA0003630065230000091
是t时的系统状态的变化量，xt是t时的系统状态，ut是t时的系统控制输入，ωt是t时的过程噪声、ωt服从均值为零的高斯分布即ωt～N(0，Qt)，Qt为t时的过程噪声协方差矩阵，t时的状态转移矩阵
Figure GDA0003630065230000092
Figure GDA0003630065230000093
是欧式空间，Bt是t时的控制输入矩阵；At具体为At(1，2)＝At(1，3)＝At(2，3)＝At(3，2)＝0，At(1，1)＝-D+μ,
Figure GDA0003630065230000094
Figure GDA0003630065230000095
At(2，2)＝-D，At(3，1)＝ωμ+β，At(3，3)＝-D，Bt具体为Bt(2，1)＝-D，Bt(1，1)＝Bt(3，1)＝0，其中D是稀释速率，μ是比生产速率，Yc/s是细胞浓度对基质浓度的生产速率，ω是产物浓度变化中与细胞浓度变化有关的参数，β是产物浓度变化中与细胞浓度相关的参数。根据Monod方程(莫诺方程式是描述微生物比增殖速度与有机基质浓度之间的函数关系)可知：
Figure GDA0003630065230000096
其中μm是最大比生产速率，Km是半饱和速度常数，表示微生物对底物的亲和力。由此很容易看出，状态方程式一个非线性方程，利用泰勒级数展开，可以将以上非线性方程转化为线性方程，线性化的连续发酵过程的状态方程为：
Figure GDA0003630065230000097
其中A′t(1，3)＝A′t(2，3)＝0，A′t(1，)＝μmθ-D，A′t(1，2)＝μmψ，
Figure GDA0003630065230000098
Figure GDA0003630065230000099
A′t(3，1)＝β+wμmθ，A′t(3，2)＝wμmψ，A′t(3，3)＝-D，
Figure GDA00036300652300000910
Ss和cs分别是稳态情况下的基质浓度和细胞浓度。
建立系统测量方程为：
yt＝Ctxt+νt；
其中yt为t时的系统测量值，Ct为t时的测量矩阵，νt为t时的测量噪声且νt服从均值为零的高斯分布即νt～N(0，Rt)，Rt为t时的测量噪声协方差矩阵；
步骤2：将上述线性化的连续发酵过程的状态方程
Figure GDA0003630065230000101
和系统测量方程yt＝Ctxt+vt所描述的连续状态空间表达式离散化，得到离散状态空间模型为：
xn＝Fn-1xn-1+Ln-1un-1+wn-1，
yn＝Hnxn+vn；
其中n表示离散状态空间的时间索引，xn是n时刻的系统状态，Fn-1是n-1时刻的系统转移矩阵，xn-1是n-1时刻的系统状态，Ln-1是n-1时刻的控制矩阵，un-1是n-1时刻的系统控制输入，wn-1是n-1时刻的过程噪声，yn是n时刻的系统测量变量，Hn是n时刻的测量矩阵，vn是n时刻的测量噪声。
Figure GDA0003630065230000102
ks是采样时间，选为1小时，为了简单起见，一般省略ks，exp{}表示指数项，A't是经过线性化的系统转移矩阵，其余矩阵和变量与连续状态空间表达式的变量取值相同。
步骤3：因为卡尔曼滤波(KF)在初始条件准确、模型准确、噪声统计特性准确时可以达到对特定工业过程的最优估计，但是在实际工业过程很难获取准确的初始工艺值。鉴于此，在卡尔曼滤波的基础上引入一个辅助变量un，使得预测分布服从学生分布，调整估计值对初始值的依赖程度，预测分布为：
Figure GDA0003630065230000103
其中un是引入的n时刻的辅助变量，y1：n-1是1时刻到n-1时刻的所有测量值序列，p(xn|un，y1：n-1)服从于均值为
Figure GDA0003630065230000104
协方差为
Figure GDA0003630065230000105
的高斯分布，p(un)表示辅助变量的先验分布、服从于形状参数为
Figure GDA0003630065230000106
逆尺度参数为
Figure GDA0003630065230000107
的伽马分布，即
Figure GDA0003630065230000108
Figure GDA0003630065230000109
是自由度参数；
Figure GDA00036300652300001010
表示均值为
Figure GDA00036300652300001011
协方差为
Figure GDA00036300652300001012
自由度为
Figure GDA0003630065230000111
的学生分布，
Figure GDA0003630065230000112
是n时刻的状态变量的预测协方差，
Figure GDA0003630065230000113
是n时刻的状态变量的预测均值，
Figure GDA0003630065230000114
是n时刻的Student-t分布的协方差，
Figure GDA0003630065230000115
是自由度参数。
通过引入一个辅助变量使得预测步符合学生分布，从而使得初值不准确时会自动调大测量值对估计值的影响，解决了卡尔曼滤波对初始值的依赖程度，提高对培养基中关键变量初始条件估计的准确性。
步骤4：使用变分贝叶斯理论将每一时刻下的系统状态和辅助变量的联合后验概率密度函数p(xn，un|y1：n)用两个独立的概率密度函数q(·)来表示，计算方法为：
p(xn，un|y1：n)≈q(xn，un|y1：n)；
其中，y1：n是从1时刻到n时刻的测量值序列；
设定两个相互独立的概率密度分别服从于高斯分布和伽马分布，即，
Figure GDA0003630065230000116
q(un|y1：n)＝G(un；an，bn)；
其中
Figure GDA0003630065230000117
是n时刻的状态变量x的均值，Pn是n时刻的状态变量x的协方差，an是n时刻的辅助变量u的形状参数，bn是n时刻的辅助变量u的逆尺度参数；
步骤5：设定系统参数初始值x0、x'0、P0、a0、b0、yn、Fn、Ln、Hn、Qn、Rn、ρ、
Figure GDA0003630065230000118
steps和N；其中x0是准确的状态变量初始均值，x'0是错误的或有不确定性信息的状态变量初始均值，P0是状态变量的初始协方差，a0是辅助变量的初始形状参数，b0是辅助变量的初始逆形状参数，yn是在n时刻的测量序列，Fn是在n时刻的状态转移矩阵，Ln是在n时刻的控制输入矩阵，Hn是n时刻的测量矩阵，Qn是n时刻的过程噪声协方差，Rn是n时刻的测量噪声协方差，ρ是启发式因子，
Figure GDA0003630065230000119
是自由度，steps是总运行步数，N是每一步迭代的总次数；
初始化离散状态空间的时间索引n＝1；
步骤6：预测n时刻的状态预测均值
Figure GDA0003630065230000121
和伽马分布的n时刻的辅助变量的预测形状参数
Figure GDA0003630065230000122
n时刻的辅助变量的预测逆尺度参数
Figure GDA0003630065230000123
n时刻的状态预测均值
Figure GDA0003630065230000124
的预测公式为：
Figure GDA0003630065230000125
其中Fn-1是在n-1时刻的状态转移矩阵,xn-1是在n-1时刻的状态变量x的均值，Ln-1是在n-1时刻的控制输入矩阵，un-1是在n-1时刻的系统控制输入；
n时刻的辅助变量的预测形状参数
Figure GDA0003630065230000126
的预测公式为：
Figure GDA0003630065230000127
n时刻的辅助变量的预测逆尺度参数
Figure GDA0003630065230000128
的预测公式为：
Figure GDA0003630065230000129
其中ρ是启发式因子，取值范围(0，1]，
Figure GDA00036300652300001210
表示状态变量在时刻n的预测均值，
Figure GDA00036300652300001211
表示辅助变量在时刻n的预测形状参数，
Figure GDA00036300652300001212
表示辅助变量在时刻n的预测逆尺度参数；
初始化迭代次数j＝1；
步骤7：第j次迭代得到的n时刻的状态变量的预测协方差
Figure GDA00036300652300001213
第j次迭代得到的n时刻的状态变量的均值
Figure GDA00036300652300001214
第j次迭代得到的n时刻的状态变量的协方差
Figure GDA00036300652300001215
和第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure GDA00036300652300001216
步骤7-1：n时刻的辅助变量的形状参数an的更新公式为：
Figure GDA00036300652300001217
其中d是状态变量的维度；
步骤:7-2：第j次迭代得到的n时刻的状态变量的预测协方差
Figure GDA00036300652300001218
的更新公式为：
Figure GDA0003630065230000131
其中
Figure GDA0003630065230000132
Pn-1是n-1时刻的状态变量的协方差，Fn-1是n-1时刻的系统转移矩阵，()T是矩阵的转置，Qn-1是在n-1时刻的过程噪声协方差矩阵，
Figure GDA0003630065230000133
是第j次迭代得到的辅助变量un的期望且
Figure GDA0003630065230000134
步骤7-3：第j次迭代得到的n时刻的状态变量的均值
Figure GDA0003630065230000135
第j次迭代得到的n时刻的状态变量的协方差
Figure GDA0003630065230000136
第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure GDA0003630065230000137
更新公式为：
Figure GDA0003630065230000138
Figure GDA0003630065230000139
Figure GDA00036300652300001310
其中，
Figure GDA00036300652300001311
是第j次迭代得到的卡尔曼滤波增益且更新公式为
Figure GDA00036300652300001312
()-1是矩阵的逆；
Figure GDA00036300652300001313
更新公式为
Figure GDA00036300652300001314
Rn是n时刻的测量噪声协方差矩阵；ψ更新公式为
Figure GDA00036300652300001315
tr[]为迹操作算子。
步骤8：判断迭代次数j是否满足j＞N，若是则执行步骤9；若否则令j＝j+1，并跳转执行步骤7；
步骤9：输出当前时刻的状态变量的均值、状态变量的协方差、辅助变量的形状参数和辅助变量的逆尺度参数，即：
当前时刻的状态变量的均值为
Figure GDA00036300652300001316
当前时刻的状态变量的协方差为
Figure GDA0003630065230000141
当前时刻的辅助变量的形状参数为
Figure GDA0003630065230000142
当前时刻的辅助变量的逆尺度参数为
Figure GDA0003630065230000143
步骤10：判断时间索引n是否满足n＞steps，若否则令n＝n+1，并跳转执行步骤6；若是则结束更新，得到初始时刻培养基中的关键变量的估计值。
本实施例中还提供了一种使用初始时刻培养基中的关键变量估计值来监测连续发酵过程的方法，通过使用连续发酵过程中估计培养基初始条件的方法来获取初始时刻培养基中的关键变量估计值，并根据初始时刻培养基中的关键变量估计值判断是否处于正常工况或者发生故障。如果估计值偏离了准确值不可接受的范围，则判定为发生故障，应及时采取措施，避免出现危险或降低发酵质量，如果估计值距离准确值的波动在可接受范围内，则认为没有发生故障。
本发明不仅可以获取传感器无法测量的变量值，而且还可以获取比传感器数据更精准的数据。除此之外，本发明还可以用来监测连续发酵过程是否发生故障，并根据具体工艺选定标准值。
为了进一步说明本发明的优点，本实施例中设置参数并进行仿真实验。设置的参数有：系统参数为μm＝0.48h-1，Km＝1.2g/L，Yx/s＝0.4g/g，sf＝20g/L，D＝0.15h-1，ω＝2.2g/g，β＝0.2h-1，cs＝7.8105，ss＝0.5455，采样时间选择为1小时，测量矩阵H选择为单位矩阵，准确的状态初始均值x0＝[7.05，50.5，1]T，错误的状态初始均值x'0＝10*[7.05，50.5，1]T，0时刻的状态初始协方差P0＝0.1*I3×3，过程噪声协方差Qn每个时刻都相同Qn＝10*I3×3，测量噪声协方差Rn每个时刻都相同Rn＝10*I3×3，自由度参数
Figure GDA0003630065230000144
辅助变量的两个参数初始值
Figure GDA0003630065230000145
时间步数steps＝100，每一步迭代次数N＝5，ρ取值为0.999。
将本发明(表示为VBIKF)与KF、PKF(详见参考文献“M.Farooq and S.Bruder,"Information type filters for tracking a maneuvering target,"in IEEETransactions on Aerospace and Electronic Systems,vol.26,no.3,pp.441-454,May1990”)、BIA-II(详见参考文献“Shunyi Zhao,Biao Huang,Trial-and-error or avoidinga guess？Initialization of the Kalman filter,Automatica,vol.121,2020”)三种方法进行仿真对比，图3是本发明和KF、PKF、BIA-II三种方法在连续发酵过程中细胞浓度的估计误差对比图。表1是本发明和KF、PKF、BIA-II三种方法在连续发酵过程中均方误差和计算时间的对比表。
算法 均方根误差 计算时间(秒)
KF 129.8395 0.0018
PKF 6.7661 0.0041
BIA-II 5.5048 0.0148
VBI-KF 5.4971 0.0081
表1本发明和KF、PKF、BIA-II三种方法在连续发酵过程中均方误差和计算时间的对比表
从表1可以看出，本发明相比较于KF、PKF、BIA-II三种方法，计算所耗用的时间短且准确性较高。同时，从图3可以看出本发明得到的状态估计误差在整个采样过程中波动不大，均方根误差较小。在此基础上，本实施例中将本发明与BIA-II在不同采样时间下的计算时间进行仿真对比实验，结果如图4所示。从图4可以看出在不同的采用时间下，本发明的计算时间明显少于BIA-II，性能优于BIA-II，进一步证明了本发明的有益效果。
本发明的上述技术方案相比现有技术具有以下优点：本发明所述的连续发酵过程中估计培养基初始条件的方法通过建立连续发酵过程的典型三变量数学模型、引入学生分布估计培养基的初始条件，解决了卡尔曼滤波对初始值的依赖程度，实现培养基中关键变量初始条件的快速且准确地估计，保证了发酵过程的安全和发酵质量。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
显然，上述实施例仅仅是为清楚地说明所作的举例，并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引申出的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (10)

1.一种连续发酵过程中估计培养基初始条件的方法，其特征在于：包括以下步骤：
步骤1：建立连续发酵过程的三变量数学模型和连续发酵过程的状态方程，利用泰勒级数展开连续发酵过程的状态方程得到线性化的连续发酵过程的状态方程，建立系统测量方程；
步骤2：将线性化的连续发酵过程的状态方程和系统测量方程所描述的连续状态空间表达式离散化，得到离散状态空间模型；
步骤3：在卡尔曼滤波的基础上引入辅助变量，得到服从学生分布的预测分布；
步骤4：使用变分贝叶斯理论将每一时刻下的系统状态和辅助变量的联合后验概率密度函数用两个独立的概率密度函数表示；
步骤5：设定系统参数初始值，总运行步数steps和每一步迭代的总次数N；
初始化离散状态空间的时间索引n＝1；
步骤6：预测n时刻的状态预测均值
Figure FDA0003630065220000011
和伽马分布的n时刻的辅助变量的预测形状参数
Figure FDA0003630065220000012
n时刻的辅助变量的预测逆尺度参数
Figure FDA0003630065220000013
初始化迭代次数j＝1；
步骤7：更新n时刻的辅助变量的形状参数an、第j次迭代得到的n时刻的状态变量的预测协方差
Figure FDA0003630065220000014
第j次迭代得到的n时刻的状态变量的均值
Figure FDA0003630065220000015
第j次迭代得到的n时刻的状态变量的协方差
Figure FDA0003630065220000016
和第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure FDA0003630065220000017
步骤8：判断迭代次数j是否满足j＞N，若是则执行步骤9；若否则令j＝j+1，并跳转执行步骤7；
步骤9：输出当前时刻的状态变量的均值、状态变量的协方差、辅助变量的形状参数和辅助变量的逆尺度参数；
步骤10：判断时间索引n是否满足n＞steps，若否则令n＝n+1，并跳转执行步骤6；若是则结束更新，得到初始时刻培养基中的关键变量的估计值。
2.根据权利要求1所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤1中的三变量数学模型，具体为：
系统状态
Figure FDA0003630065220000021
系统输入u＝sf；
其中c是培养基中微生物细胞浓度，s是培养基中基质浓度，p是培养基中产物浓度，sf是培养基中补给基质浓度；
所述连续发酵过程的状态方程为：
Figure FDA0003630065220000022
其中t表示连续状态空间表达式的时间索引，
Figure FDA0003630065220000023
是t时的系统状态的变化量，xt是t时的系统当前状态，ut是t时的系统控制输入，ωt是t时的过程噪声、ωt服从均值为零的高斯分布即ωt～N(0，Qt)，Qt为t时的过程噪声协方差矩阵，t时的状态转移矩阵
Figure FDA0003630065220000024
Figure FDA0003630065220000025
是欧式空间，Bt是t时的控制输入矩阵；
At具体为：At(1，2)＝At(1，3)＝At(2，3)＝At(3，2)＝0，At(1，1)＝-D+μ，
Figure FDA0003630065220000026
At(2，2)＝-D，At(3，1)＝wμ+β，At(3，3)＝-D；
Bt具体为：Bt(2，1)＝-D，Bt(1，1)＝Bt(3，1)＝0；
其中D是稀释速率，μ是比生产速率，Yc/s是细胞浓度对基质浓度的生产速率，ω是产物浓度变化中与细胞浓度变化有关的参数，β是产物浓度变化中与细胞浓度相关的参数；
Figure FDA0003630065220000031
其中μm是最大比生产速率，Km是半饱和速度常数。
3.根据权利要求2所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤1中的线性化的连续发酵过程的状态方程，具体为：
Figure FDA0003630065220000032
其中A′t(1，3)＝A′t(2，3)＝0，
Figure FDA0003630065220000033
A′t(1，2)＝μmψ，
Figure FDA0003630065220000034
Figure FDA0003630065220000035
A′t(3，2)＝wμmψ，A′t(3，3)＝-D，
Figure FDA0003630065220000036
ss是稳态情况下的基质浓度，cs是稳态情况下细胞浓度；
所述步骤1中建立的系统测量方程，具体为：
yt＝Ctxt+vt；
其中yt为t时的系统测量值，Ct为t时的系统系数矩阵，vt为t时的测量噪声且vt服从均值为零的高斯分布。
4.根据权利要求3所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤2中的离散状态空间模型，具体为：
xn＝Fn-1xn-1+Ln-1un-1+wn-1，
yn＝Hnxn+vn；
其中xn是n时刻的系统状态，Fn-1是n-1时刻的系统转移矩阵，xn-1是n-1时刻的系统状态，Ln-1是n-1时刻的控制矩阵，un-1是n-1时刻的系统控制输入，wn-1是n-1时刻的过程噪声，yn是n时刻的系统测量变量，Hn是n时刻的测量矩阵，vn是n时刻的测量噪声。
5.根据权利要求4所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤3中的服从学生分布的预测分布，具体为：
Figure FDA0003630065220000041
其中un是引入的n时刻的辅助变量，y1：n-1是1时刻到n-1时刻的所有测量值序列，p(xn|un，y1：n-1)服从于均值为
Figure FDA0003630065220000042
协方差为
Figure FDA0003630065220000043
的高斯分布，p(un)表示辅助变量un的先验分布、服从于形状参数为
Figure FDA0003630065220000044
逆尺度参数为
Figure FDA0003630065220000045
的伽马分布；
Figure FDA0003630065220000046
表示均值为
Figure FDA0003630065220000047
协方差为
Figure FDA0003630065220000048
自由度为
Figure FDA00036300652200000415
的学生分布，
Figure FDA0003630065220000049
是n时刻的状态变量的预测协方差，
Figure FDA00036300652200000410
是n时刻的状态变量的预测均值，
Figure FDA00036300652200000411
是n时刻的Student-t分布的协方差，
Figure FDA00036300652200000412
是自由度参数。
6.根据权利要求5所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤4中使用变分贝叶斯理论将每一时刻下的系统状态和辅助变量的联合后验概率密度函数用两个独立的概率密度函数表示，其中两个独立的概率密度函数分别服从于高斯分布和伽马分布，具体为：
p(xn，un|y1：n)≈q(xn，un|y1：n)＝q(xn|y1：n)q(un|y1：n)；
其中，p(xn，un|y1：n)是每一时刻下的系统状态和辅助变量的联合后验概率密度函数，q(·)是独立的概率密度函数，y1：n是1时刻到n时刻的所有测量值序列；
Figure FDA00036300652200000413
q(un|y1：n)＝G(un；an，bn)；
其中
Figure FDA00036300652200000414
是n时刻的状态变量的均值，Pn是n时刻的状态变量的协方差。
7.根据权利要求6所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤6中预测n时刻的状态预测均值
Figure FDA0003630065220000051
和伽马分布的n时刻的辅助变量的预测形状参数
Figure FDA0003630065220000052
n时刻的辅助变量的预测逆尺度参数
Figure FDA0003630065220000053
具体预测公式为：
n时刻的状态预测均值
Figure FDA0003630065220000054
的预测公式为：
Figure FDA0003630065220000055
其中xn-1是n-1时刻的状态变量的均值；
n时刻的辅助变量的预测形状参数
Figure FDA0003630065220000056
的预测公式为：
Figure FDA0003630065220000057
n时刻的辅助变量的预测逆尺度参数
Figure FDA0003630065220000058
的预测公式为：
Figure FDA0003630065220000059
其中ρ是启发式因子，取值范围(0，1]，an-1是n-1时刻的辅助变量的形状参数，bn-1是n一1时刻的辅助变量的逆尺度参数。
8.根据权利要求7所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤7中更新n时刻的辅助变量的形状参数an、第j次迭代得到的n时刻的状态变量的预测协方差
Figure FDA00036300652200000510
第j次迭代得到的n时刻的状态变量的均值
Figure FDA00036300652200000511
第j次迭代得到的n时刻的状态变量的协方差
Figure FDA00036300652200000512
和第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure FDA00036300652200000513
具体为：
n时刻的辅助变量的形状参数an的更新公式为：
Figure FDA00036300652200000514
其中d是状态变量的维度；
第j次迭代得到的n时刻的状态变量的预测协方差
Figure FDA00036300652200000515
的更新公式为：
Figure FDA00036300652200000516
其中
Figure FDA00036300652200000517
Pn-1是n-1时刻的状态变量的协方差，Fn-1是n-1时刻的系统转移矩阵，()T是矩阵的转置，Qn-1是n-1时刻的过程噪声协方差矩阵；
Figure FDA00036300652200000518
是第j次迭代得到的辅助变量un的期望且
Figure FDA00036300652200000519
第j次迭代得到的n时刻的状态变量的均值
Figure FDA00036300652200000520
的更新公式为：
Figure FDA0003630065220000061
第j次迭代得到的n时刻的状态变量的协方差
Figure FDA0003630065220000062
的更新公式为：
Figure FDA0003630065220000063
第j次迭代得到的n时刻的辅助变量的预测逆尺度参数
Figure FDA0003630065220000064
的更新公式为：
Figure FDA0003630065220000065
其中
Figure FDA0003630065220000066
是第j次迭代得到的n时刻的卡尔曼滤波增益且更新公式为
Figure FDA0003630065220000067
()-1是矩阵的逆，
Figure FDA0003630065220000068
的更新公式为
Figure FDA0003630065220000069
Rn是n时刻的测量噪声协方差矩阵；ψ更新公式为
Figure FDA00036300652200000610
tr[]为迹操作算子。
9.根据权利要求8所述的连续发酵过程中估计培养基初始条件的方法，其特征在于：所述步骤9中输出的当前时刻的状态变量的均值、状态变量的协方差、辅助变量的形状参数和辅助变量的逆尺度参数，具体为：
当前时刻的状态变量的均值为
Figure FDA00036300652200000611
当前时刻的状态变量的协方差为
Figure FDA00036300652200000612
当前时刻的辅助变量的形状参数为
Figure FDA00036300652200000613
当前时刻的辅助变量的逆尺度参数为
Figure FDA00036300652200000614
10.一种使用初始时刻培养基中的关键变量估计值来监测连续发酵过程的方法，其特征在于：使用权利要求1-9任一项所述的连续发酵过程中估计培养基初始条件的方法来获取初始时刻培养基中的关键变量估计值，并根据初始时刻培养基中的关键变量估计值判断是否处于正常工况或者发生故障。
CN202110615110.9A 2021-06-02 2021-06-02 连续发酵过程中估计培养基初始条件的方法 Active CN113326618B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110615110.9A CN113326618B (zh) 2021-06-02 2021-06-02 连续发酵过程中估计培养基初始条件的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110615110.9A CN113326618B (zh) 2021-06-02 2021-06-02 连续发酵过程中估计培养基初始条件的方法

Publications (2)

Publication Number Publication Date
CN113326618A CN113326618A (zh) 2021-08-31
CN113326618B true CN113326618B (zh) 2022-07-15

Family

ID=77421480

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110615110.9A Active CN113326618B (zh) 2021-06-02 2021-06-02 连续发酵过程中估计培养基初始条件的方法

Country Status (1)

Country Link
CN (1) CN113326618B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113790727A (zh) * 2021-09-07 2021-12-14 中国西安卫星测控中心 一种基于辅助状态参数的脉冲机动检测方法
CN113963752A (zh) * 2021-10-20 2022-01-21 江南大学 一种抗生素发酵过程生物量浓度的估计方法及系统
CN114036810A (zh) * 2021-11-04 2022-02-11 江南大学 一种细胞培养状态在线估计及优化补料调控方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE19547473A1 (de) * 1995-12-19 1997-06-26 Wissenschaftsfoerderung Der De Verfahren zur Echtzeit-Ermittlung der Konzentration eines Sekundär-Stoffwechselprodukts in einem chargenbetriebenen Fermentationsprozeß
EP3296387A1 (de) * 2016-09-16 2018-03-21 Siemens Aktiengesellschaft Verfahren zur überwachung von bioprozessen
CN109669132A (zh) * 2019-01-21 2019-04-23 西北工业大学 一种基于变分贝叶斯滤波的电池荷电状态估计方法
CN111581905A (zh) * 2020-05-15 2020-08-25 江南大学 隧道二极管电路系统在未知测量噪声下的状态估计方法
CN112418051A (zh) * 2020-11-18 2021-02-26 温州大学 一种用于非线性动态系统非高斯噪声下的状态估计方法
CN112528479A (zh) * 2020-12-01 2021-03-19 哈尔滨工程大学 一种基于Gibbs采样器的鲁棒自适应平滑方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5991525A (en) * 1997-08-22 1999-11-23 Voyan Technology Method for real-time nonlinear system state estimation and control
CN109829938B (zh) * 2019-01-28 2020-12-08 杭州电子科技大学 一种应用在目标跟踪的自适应容错容积卡尔曼滤波方法
CN110957011B (zh) * 2019-11-25 2023-03-17 江南大学 连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE19547473A1 (de) * 1995-12-19 1997-06-26 Wissenschaftsfoerderung Der De Verfahren zur Echtzeit-Ermittlung der Konzentration eines Sekundär-Stoffwechselprodukts in einem chargenbetriebenen Fermentationsprozeß
EP3296387A1 (de) * 2016-09-16 2018-03-21 Siemens Aktiengesellschaft Verfahren zur überwachung von bioprozessen
CN109669132A (zh) * 2019-01-21 2019-04-23 西北工业大学 一种基于变分贝叶斯滤波的电池荷电状态估计方法
CN111581905A (zh) * 2020-05-15 2020-08-25 江南大学 隧道二极管电路系统在未知测量噪声下的状态估计方法
CN112418051A (zh) * 2020-11-18 2021-02-26 温州大学 一种用于非线性动态系统非高斯噪声下的状态估计方法
CN112528479A (zh) * 2020-12-01 2021-03-19 哈尔滨工程大学 一种基于Gibbs采样器的鲁棒自适应平滑方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"基于批次加权正则极限学习机的发酵过程软测量";姚景升;刘飞;《江南大学学报(自然科学版)》;20131028;第12卷(第5期);第515-521页 *
Shuang Gao ; Shunyi Zhao ; Xiaoli Luan ; Fei Liu."Intelligent State Estimation for Continuous Fermenters Using Variational Bayesian Learning".《IEEE Transactions on Industrial Informatics》.2021, *

Also Published As

Publication number Publication date
CN113326618A (zh) 2021-08-31

Similar Documents

Publication Publication Date Title
CN113326618B (zh) 连续发酵过程中估计培养基初始条件的方法
CN108284442B (zh) 一种基于模糊神经网络的机械臂柔性关节控制方法
Wang et al. Sensitivity analysis and identification of kinetic parameters in batch fermentation of glycerol
CN111898867B (zh) 一种基于深度神经网络的飞机总装生产线产能预测方法
CN112818595A (zh) 一种火电厂蒸发区的数字孪生模型数据的修正方法及系统
CN105867138A (zh) 一种基于pid控制器的稳定平台控制方法及装置
KR20170109629A (ko) 발효 모델을 생성하기 위한 컴퓨터 구현 방법
WO2023019883A1 (zh) 利用细胞代谢网络监测生物制造过程的方法
CN106605179A (zh) 预测值整形系统、控制系统、预测值整形方法、控制方法及预测值整形程序
CN108959787B (zh) 考虑实际工况的宏宏双驱动系统的热变形预测方法及系统
Retamal et al. Parameter estimation of a dynamic model of Escherichia coli fed-batch cultures
He et al. Quantifying dynamic regulation in metabolic pathways with nonparametric flux inference
Luo et al. Iterative improvement of parameter estimation for model migration by means of sequential experiments
CN115542236B (zh) 电能表运行误差估计方法及装置
Salau et al. State estimators for better bioprocesses operation
CN110794676A (zh) 基于Hammerstein-Wiener模型的CSTR过程非线性控制方法
CN115305526A (zh) 基于x射线测量的铜箔厚度面密度一致性自适应控制方法
Gjerkes et al. Product identification in industrial batch fermentation using a variable forgetting factor
CN110309472B (zh) 基于离线数据的策略评估方法及装置
CN112667957A (zh) 一种基于深度神经网络的智能电能表失效率预测方法
Ma et al. Intelligent recommendation system of the injection molding process parameters based on CAE simulation, process window, and machine learning
CN111210877A (zh) 一种推断物性参数的方法及装置
Rao et al. Parallel solution of DDDAS variational inference problems
Lisci et al. A Robust Nonlinear Estimator for a Yeast Fermentation Biochemical Reactor
CN113158394B (zh) 基于蒸发误差补偿的流域水文模型误差修正方法和装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant