CN105199973A - 一种基于差分进化算法的酵母培养在线自适应控制方法 - Google Patents

Abstract

本发明公开了一种基于差分进化算法的酵母培养在线自适应控制方法，属于生物技术领域。本发明方法以标准PID控制策略为基础，利用差分进化算法在线优化PID控制参数，并根据乙醇浓度的变化，自动调整底物流加速率，使发酵液中底物浓度维持在某一合适的范围内。本发明方法既避免了底物浓度过高引起Crabtree效应、积累大量乙醇、损伤细胞，也避免了底物浓度过低导致的细胞生长速率降低。采用本发明所述的控制方法，乙醇浓度能够稳定地控制在设定值附近，只在Ce_set±0.3g/L的小范围内波动，同时菌体也快速生长，发酵结束时菌体浓度能够比传统PID控制策略提高25％以上。

Description

一种基于差分进化算法的酵母培养在线自适应控制方法

技术领域

本发明涉及一种基于差分进化算法的酵母培养在线自适应控制方法，属于生物技术领域。

背景技术

酵母是与人们生活息息相关的一种真核微生物，可用于生产多种产物，具有遗传背景清楚，不产内毒素，可对外源蛋白进行翻译后修饰，发酵周期短、操作方便等优点。酿酒酵母(Saccharomycescerevisiae)、鲁氏接合酵母(Zygosaccharomycesrouxii)和拜耳接合酵母(Zygosaccharomycesbailii)等酵母培养过程具有明显的Crabtree效应，即当碳源浓度过高时，大量碳源流向乙醇合成途径，导致乙醇的大量积累。积累乙醇的同时，伴随着大量CO₂的排放，大大降低碳源的利用效率，且高浓度的乙醇也会对细胞结构及功能造成不可逆转的损伤，抑制产物合成。减少碳源添加，可有效缓解乙醇积累，但过低的碳源浓度同样会影响酵母细胞的生长和产物的合成。利用基因工程手段敲除乙醇合成途径中的关键酶能够有效降低乙醇积累，但此类工程菌对碳源浓度过于敏感，在发酵生产中的应用十分有限。利用过程控制技术，合理调节碳源添加，在保证细胞正常生长的前提下抑制乙醇积累，是解决这一难题的有效途径之一。发酵过程具有高度的时变性、滞后性和非线性，利用控制参数恒定的常规比例-积分-微分(PID)控制对发酵过程实施控制，往往难以取得令人满意的效果。

于是，能够自动实时调节控制参数的在线自适应控制策略被越来越多地用于发酵过程控制。差分进化算法是由R.Stom和K.Price建立的一种全局优化算法，其原理简单，受控参数少，收敛速度快，已受到各领域研究人员的广泛关注。在此基础上，本发明提出了基于差分进化算法的在线自适应控制策略，以期解决上述酵母培养过程中存在的难题。

发明内容

酵母培养中，底物浓度过高与不足都不利于细胞的生长，为了将底物浓度控制在某一适宜的水平，促进细胞生长，本发明提出了一种基于差分进化算法、自回归移动平均模型与PID控制策略相结合的在线自适应控制策略。本发明以标准PID控制策略为基础，利用差分进化算法在线优化PID控制参数，并根据乙醇浓度的变化，自动调整底物流加速率，使发酵液中底物浓度维持在某一合适的范围内。

本发明提供的一种基于差分进化算法的酵母培养在线自适应控制方法，包括以下步骤：

步骤1：初始化，设定目标乙醇浓度Ce_set、初始底物流加速率F和比例-积分-微分控制(即PID控制)的初始参数：比例系数k_p、积分系数k_i和微分系数k_d；

步骤2：采集数据，采集酵母培养过程中乙醇浓度的设定值Ce_set和检测值Ce、底物流加速率F、PID控制参数；

步骤3：系统辨识，根据当前时刻以及前两个时刻的乙醇浓度偏差、前一时刻的底物流加速率及PID控制参数，预测后两个时刻的乙醇浓度；

步骤4：优化PID控制参数，以最小化后两个采样时刻的乙醇预测浓度偏差为目标，利用差分进化算法(Differentialevolution，DE)在线优化PID控制参数；

步骤5：优化底物流加速率，根据当前时刻的乙醇浓度偏差，利用优化后的PID控制参数计算底物流加速率。

在本发明的一种实施方式中，利用商业化的甲醇电极在线检测乙醇浓度，通过A/D数据转换卡传输并保存在工控机，采集数据时直接读取。底物流加速率和PID控制参数由自适应控制系统计算并保存在工控机，采集数据时也直接读取。

在本发明的一种实施方式中，所述步骤3是利用自回归移动平均模型(Autoregressivemovingaveragemodel，ARMA)对酵母培养系统进行辨识；过程如下：

(1)当采样周期小于4(即k<4)时，用标准PID控制策略控制底物流加，控制参数不变；当采样周期大于等于4，即k≥4时，用ARMA模型描述乙醇浓度和流加速率之间的关系：

Ce(k)＝-aCe(k-1)+bF(k-1)+e(k)

其中，乙醇浓度Ce为模型的输出；流加速率F为模型的输入；k为当前采样周期；a、b为模型参数；e为模型的误差值。

(2)利用最小二乘法原理估计ARMA模型参数；

(3)用前一时刻的流加速率、PID控制参数和乙醇浓度偏差表示当前时刻的流加速率：

F(k)＝F(k-1)+k_p(k)[e(k)-e(k-1)]+k_i(k)e(k)+k_d(k)[e(k)+e(k-2)-2e(k-1)]

其中，k_p(k)、k_i(k)、k_d(k)通过下式进行更新，其中p₁、p₂和p₃分别是调节k_p、k_i和k_d的步长；

k_p(k)＝(1+p₁)k_p(k-1)-0.1≤p₁≤0.1

k_i(k)＝(1+p₂)k_i(k-1)-0.1≤p₂≤0.1

k_d(k)＝(1+p₃)k_d(k-1)-0.1≤p₃≤0.1

(4)预测后两个时刻的乙醇浓度：

令X'＝[-Ce(k)F(k)]，则k+1采样周期乙醇浓度的预测值为：再令则k+2采样周期乙醇浓度的预测值为：其中θ为模型参数向量。

在本发明的一种实施方式，所述利用最小二乘法原理估计ARMA模型参数，具体是：

(k采样周期，J是一个函数的符号，i是小于k的自然数)；

定义矩阵X、Y、W和模型参数向量θ：

则θ＝(X^TX)^-1X^TY；其中λ为忘却因子且0<λ<1，T代表矩阵的转置，W为权重系数矩阵。

在本发明的一种实施方式，所述步骤4中以最小化后两个采样时刻的乙醇预测浓度偏差为目标，目标函数如下

$\underset{p_1,p_2,p_3}{\min }G(p_1,p_2,p_3)=\mathrm{\&alpha;}|{\mathrm{Ce}}_{set}-\hat{C}e(k+1)|+\mathrm{\&beta;}\&times;|{\mathrm{Ce}}_{set}-\hat{C}e(k+2)|,;$ 其中α和β是k+1和k+2时刻预测值误差值的权重系数。

在本发明的一种实施方式，所述步骤4中利用差分进化算法在线优化PID控制参数，过程如下：

(a)生成初始种群

随机生成种群大小为NP的初始种群，种群中的个体由三维向量 $P_i^t=\&lsqb;\begin{array}{ccc}{p}_{i,1}^t& p_{i,2}^t& p_{i,3}^t\end{array}\&rsqb;$ 表示，其中i∈[1,2,3,…,NP]，表示该个体在种群中的索引号，t表示传代次数，在初始种群中 $\begin{array}{cccc}t=0;& p_{i,1}^t& p_{i,2}^t& p_{i,3}^t\end{array}$ 代表种群中个体的三个分量；计算每个个体所对应的适应度函数值G(P_i ^t)并将适应度函数最小的个体定义为最优解P_best，P_best＝[p_best1p_best2p_best3]；

(b)生成变异矢量

随机选取除目标矢量P_i ^t外的三个不同个体以及按照生成变异矢量V_i＝[v₁,v₂,v₃]，其中F_DE为缩放因子；v₁,v₂,v₃代表变异矢量的三个分量；

(c)生成试验矢量

定义U_i＝[u₁,u₂,u₃]为试验矢量，按照下式给U_i中各分量赋值，其中j＝1，2，3；R_j为0～1之间均匀分布的随机数；CR(0<CR<1)为交叉概率；RN∈[1,2,3]；u₁,u₂,u₃代表试验矢量的三个分量；

(d)生成新个体

将试验矢量U_i与目标矢量P_i ^t进行比较，若U_i优于P_i ^t则用U_i代替P_i ^t，进入下一代： $P_i^{t+1}=\left\{\begin{array}{cc}{U}_i& G\left(U_i\right)<G\left(P_i^t\right)\\ P_i^t& G\left(U_i\right)\&GreaterEqual;G\left(P_i^t\right)\end{array}\right.;$

按照(b)、(c)、(d)的运算过程对每个个体进行传代，得到第t+1代种群，计算各个体所对应的适应度函数G(P_i ^t+1)若min{G(P_i ^t+1)}＜G(P_best)，则由min{G(P_i ^t+1)}对应的个体代替原有P_best，否则P_best维持不变；

(e)迭代和优化PID控制参数

按照上述步骤(a)～(d)进行多次迭代，若P_best连续n代没有更新，则适应度函数已经收敛，利用下式更新PID控制参数：

k_p(k)＝(1+P_best1)k_p(k-1)

k_i(k)＝(1+P_best1)k_i(k-1)

k_d(k)＝(1+P_best1)k_d(k-1)

在本发明的一种实施方式，所述步骤5是在步骤4更新PID控制参数后，按下式更新底物流加速率F(k)：

F(k)＝F(k-1)+k_p(k)[e(k)-e(k-1)]+k_i(k)e(k)+k_d(k)[e(k)+e(k-2)-2e(k-1)]。

在本发明的一种实施方式，目标乙醇浓度Ce_set、初始底物流加速率F和初始PID控制参数k_p、k_i、k_d，它们的取值分别为1g/L、0.002L/h、0.0017、0.0599和0.0224。

本发明的有益效果：

本发明将DE算法和基于ARMA模型的系统辨识与传统PID控制策略相结合，通过ARMA模型在线辨识发酵系统的实时状态，并利用DE算法在线调节PID控制参数，使得该控制策略能够很好地应对生物反应过程的滞后性、非线性和时变型特征。本发明以乙醇浓度作为输入变量，通过将乙醇浓度设定在一个较低的水平，根据乙醇浓度的检测值与设定值的偏差，在线反馈调节底物的流加速率，从而将发酵液内底物浓度维持在一个合适的范围内。既避免了底物浓度过高引起Crabtree效应，积累大量乙醇，损伤细胞，也避免了底物浓度过低，降低细胞生长速率。采用本发明所述的控制策略，乙醇浓度能够稳定地控制在设定值附近，只在Ce_set±0.3g/L的小范围内波动，同时菌体也快速生长，发酵结束时菌体浓度能够比传统PID控制策略提高25％以上。

附图说明

图1：DE-PID控制策略示意图；

图2：PID控制示意图。

具体实施方式

实施例：基于差分进化算法的在线自适应控制策略用于酿酒酵母培养

基于差分进化算法的酿酒酵母培养在线自适应控制框图如图1所示，主要由系统辨识模块、参数优化模块和PID控制器三部分组成。先由系统辨识模块对酿酒酵母培养系统进行辨识，然后参数优化模块根据辨识结果对PID控制器的控制参数进行优化，最后PID控制器再对底物流加速率进行优化。具体实施过程如下：

(1)设定目标乙醇浓度Ce_set、初始底物流加速率F和初始PID控制参数k_p、k_i、k_d，它们的取值分别为1g/L、0.002L/h、0.0017、0.0599和0.0224。

(2)采集酿酒酵母培养过程中乙醇浓度的设定值和检测值、底物流加速率、PID控制器的控制参数。其中乙醇浓度用商业化甲醇电极在线检测，通过A/D数据转换卡传输并保存在工控机，采集数据时直接读取。底物流加速率和PID控制参数由自适应控制系统计算并保存在工控机，采集数据时也直接读取。

(3)前3个采样周期，采用传统PID控制策略，从第4个周期开始进行系统辨识和PID参数优化，用ARMA模型描述乙醇浓度和流加速率之间的关系：

Ce(k)＝-aCe(k-1)+bF(k-1)+e(k)；其中，乙醇浓度Ce为模型的输出；流加速率F为模型的输入；k为当前采样点；a、b为模型参数；e为模型的误差值。

(4)用最小二乘法原理估计ARMA模型参数，该过程可表示为：

$\underset{a,b}{min}J=\frac{1}{k}\underset{i=2}{\overset{k}{\&Sigma;}}{\&lsqb;Ce\left(i\right)+aCe(i-1)-bF(i-1)\&rsqb;}^2.$

定义矩阵X、Y、W和模型参数向量θ：

则θ＝(X^TX)^-1X^TY。

其中λ＝0.95为忘却因子，W为权重系数矩阵。

(5)用前一时刻的流加速率、PID控制参数和乙醇浓度偏差表示当前时刻的流加速率：

F(k)＝F(k-1)+k_p(k)[e(k)-e(k-1)]+k_i(k)e(k)+k_d(k)[e(k)+e(k-2)-2e(k-1)]

其中，k_p(k)、k_i(k)、k_d(k)通过下式进行更新：

k_p(k)＝(1+p₁)k_p(k-1)-0.1≤p₁≤0.1

k_i(k)＝(1+p₂)k_i(k-1)-0.1≤p₂≤0.1

k_d(k)＝(1+p₃)k_d(k-1)-0.1≤p₃≤0.1

(6)预测后两个时刻的乙醇浓度：

令X'＝[-Ce(k)F(k)]，则k+1采样点乙醇浓度的预测值为：再令 $X^{\&prime;\&prime;}=\&lsqb;\begin{array}{cc}-\hat{C}e(k+1)& F\left(k\right)\end{array}\&rsqb;,$ 则k+2采样点乙醇浓度的预测值为： $\hat{C}e(k+2)=X^{\&prime;\&prime;}\mathrm{\&theta;}.$

(7)用差分进化算法在线优化PID控制参数，使后两个时刻的乙醇浓度偏差最小，因此目标函数为：

$\underset{p_1,p_2,p_3}{\min }G(p_1,p_2,p_3)=0.7\&times;|{\mathrm{Ce}}_{set}-\hat{C}e(k+1)|+0.3\&times;|{\mathrm{Ce}}_{set}-\hat{C}e(k+2)|$

优化过程如下：

a生成初始种群

随机生成种群大小为NP＝100的初始种群，种群中的个体由三维向量 $P_i^t=\&lsqb;\begin{array}{ccc}{p}_{i,1}^t& p_{i,2}^t& p_{i,3}^t\end{array}\&rsqb;$ 表示，其中i∈[1,2,3,…,NP]，表示该个体在种群中的索引号，t表示传代次数，在初始种群中t＝0。计算每个个体所对应的适应度函数值G(P_i ^t)值并将适应度函数最小的个体定义为最优解P_best，P_best＝[p_best1p_best2p_best3]。

b生成变异矢量

随机选取除目标矢量P_i ^t外的三个不同个体以及按下式生成变异矢量V_i＝[v_1,v_2,v₃]，其中F_DE＝0.8为缩放因子。

$V_i=P_{r_3}^t+F_{DE}\&times;(P_{r_1}^t-P_{r_2}^t)$

c生成试验矢量

定义U_i＝[u₁,u₂,u₃]为试验矢量，按照下式给U_i中各分量赋值，其中j＝1，2，3；R_j为0～1之间均匀分布的随机数；CR＝0.8为交叉概率；RN∈[1,2,3]。

d生成新个体

将试验矢量U_i与目标矢量P_i ^t进行比较，若U_i优于P_i ^t则用U_i代替P_i ^t进入下一代： $P_i^{t+1}=\left\{\begin{array}{cc}{U}_i& G\left(U_i\right)<G\left(P_i^t\right)\\ P_i^t& G\left(U_i\right)\&GreaterEqual;G\left(P_i^t\right)\end{array}\right..$

按照b、c、d的运算过程对每个个体进行传代，得到第t+1代种群，计算各个体所对应的适应度函数G(P_i ^t+1)若min{G(P_i ^t+1)}＜G(P_best)，则由min{G(P_i ^t+1)}对应的个体代替原有P_best，否则P_best维持不变。

e迭代和优化PID控制参数

按照上述步骤进行多次迭代，若P_best连续n代没有更新，则适应度函数已经收敛，利用下式更新PID控制参数：

k_p(k)＝(1+P_best1)k_p(k-1)

k_i(k)＝(1+P_best1)k_i(k-1)

k_d(k)＝(1+P_best1)k_d(k-1)

(8)优化PID控制参数后，按下式更新底物流加速率F(k)。

F(k)＝F(k-1)+k_p(k)[e(k)-e(k-1)]+k_i(k)e(k)+k_d(k)[e(k)+e(k-2)-2e(k-1)]。

以酿酒酵母为例，比较了本发明所述的控制策略与传统PID控制策略下酵母的生长情况。

培养条件：酿酒酵母培养在30℃、pH6.0的条件下进行，培养过程中用到发酵初始培养基和流加培养基；其中发酵初始培养基(g/L)：葡萄糖5，酵母粉12，(NH₄)₂SO₄2.5，MgSO₄·7H₂O1.5，CaCl₂·2H₂O0.1；流加培养基(g/L)：葡萄糖500。

传统PID控制策略如下：

PID控制包括比例、积分和微分3个环节，可由下式表示:

$u\left(t\right)=K_p\left(e\right(t)+\frac{1}{T_i}{\&Integral;}_0^{t}e(t)dt+T_d\frac{de\left(t\right)}{dt})$

离散化后可以改写为：

$u\left(k\right)=k_{p}e\left(k\right)+k_i\underset{j=0}{\overset{k}{\&Sigma;}}e\left(j\right)+k_d\&lsqb;e\left(k\right)-e(k-1)\&rsqb;$

递推可得：

$u(k-1)=k_{p}e(k-1)+k_i\underset{j=0}{\overset{k-1}{\&Sigma;}}e\left(j\right)+k_d\&lsqb;e(k-1)-e(k-2)\&rsqb;$

两式相减得：

Δu(k)＝k_p[e(k)-e(k-1)]+k_ie(k)+k_d[e(k)-2e(k-1)+e(k-2)]

将该策略用于本发明，根据乙醇浓度反馈调节葡萄糖流加速率F，则F可由下式计算得到，图2为PID控制框图。

F(k)＝F(k-1)+ΔF(k)

F(k)＝F(k-1)+k_p[e(k)-e(k-1)]+k_ie(k)+k_d[e(k)-2e(k-1)+e(k-2)]

式中k为当前采样点；e为乙醇浓度偏差，即e(k)＝Ce_set-Ce(k)；k_p、k_i、k_d为控制参数，由操作人员设定后，在控制过程中其取值不发生变化。

结果发现，采用本发明的控制策略，乙醇浓度能够稳定地控制在设定值(1g/L)附近，只在1±0.2g/L的小范围内波动，同时菌体也快速生长，发酵结束时菌体浓度达到34.45g/L，菌体浓度比传统PID控制策略提高25％以上。

虽然本发明已以较佳实施例公开如上，但其并非用以限定本发明，任何熟悉此技术的人，在不脱离本发明的精神和范围内，都可做各种的改动与修饰，因此本发明的保护范围应该以权利要求书所界定的为准。

Claims (7)

1.一种酵母培养在线自适应控制方法，其特征在于，所述方法是结合了差分进化算法、自回归移动平均模型与PID控制策略；所述方法包括如下步骤：

步骤1：初始化，设定目标乙醇浓度Ce_set、初始底物流加速率F和PID控制的初始参数：比例系数k_p、积分系数k_i和微分系数k_d；

步骤2：采集数据，采集酵母培养过程中乙醇浓度的设定值和检测值、底物流加速率、PID控制参数；

步骤3：系统辨识，根据当前时刻以及前两个时刻的乙醇浓度偏差、前一时刻的底物流加速率及PID控制参数，预测后两个时刻的乙醇浓度；

步骤4：优化PID控制参数，以最小化后两个时刻的乙醇浓度偏差为目标，利用差分进化算法在线优化PID控制参数；

步骤5：优化底物流加速率，根据当前时刻的乙醇浓度偏差，利用优化后的PID控制参数计算底物流加速率。

2.根据权利要求1所述的方法，其特征在于，所述步骤3是利用ARMA模型对酵母培养系统进行辨识；过程如下：

(1)当采样周期大于等于4，即k≥4时，用ARMA模型描述乙醇浓度和流加速率之间的关系：

Ce(k)＝-aCe(k-1)+bF(k-1)+e(k)

其中，乙醇浓度Ce为模型的输出；流加速率F为模型的输入；k为当前采样周期；a、b为模型参数；e为模型的误差值；

(2)利用最小二乘法原理估计ARMA模型参数；

(3)用前一时刻的流加速率、PID控制参数和乙醇浓度偏差表示当前时刻的流加速率：

F(k)＝F(k-1)+k_p(k)[e(k)-e(k-1)]+k_i(k)e(k)+k_d(k)[e(k)+e(k-2)-2e(k-1)]

其中，k_p(k)、k_i(k)、k_d(k)通过下式进行更新：

k_p(k)＝(1+p₁)k_p(k-1)-0.1≤p₁≤0.1

k_i(k)＝(1+p₂)k_i(k-1)-0.1≤p₂≤0.1

k_d(k)＝(1+p₃)k_d(k-1)-0.1≤p₃≤0.1

(4)预测后两个时刻的乙醇浓度：

令X'＝[-Ce(k)F(k)]，则k+1采样点乙醇浓度的预测值为：再令则k+2采样点乙醇浓度的预测值为：其中θ为模型参数向量。

3.根据权利要求2所述的方法，其特征在于，所述(1)中，当k<4时用标准PID控制策略控制底物流加，控制参数不变，当k≥4时开始进行系统辨识，调整控制参数。

4.根据权利要求1所述的方法，其特征在于，所述利用最小二乘法原理估计ARMA模型参数，具体是：

定义矩阵X、Y、W和模型参数向量θ：

则θ＝(X^TX)^-1X^TY；其中λ为忘却因子且0<λ<1，W为权重系数矩阵。

5.根据权利要求1所述的方法，其特征在于，所述步骤4中以最小化后两个时刻的乙醇浓度偏差为目标，目标函数如下

其中α和β是k+1和k+2时刻预测值误差值的权重系数。

6.根据权利要求1所述的方法，其特征在于，所述步骤4中利用差分进化算法在线优化PID控制参数，过程如下：

(a)生成初始种群

随机生成种群大小为NP的初始种群，种群中的个体由三维向量表示，其中i∈[1,2,3,…,NP]，表示该个体在种群中的索引号，t表示传代次数，在初始种群中t＝0。计算每个个体所对应的适应度函数值G(P_i ^t)，并将适应度函数最小的个体定义为最优解P_best，P_best＝[p_best1p_best2p_best3]；

(b)生成变异矢量

随机选取除目标矢量P_i ^t外的三个不同个体以及按照生成变异矢量V_i＝[v₁,v₂,v₃]，其中F_DE为缩放因子；

(c)生成试验矢量

定义U_i＝[u₁,u₂,u₃]为试验矢量，按照下式给U_i中各分量赋值，其中j＝1，2，3；R_j为0～1之间均匀分布的随机数；CR(0<CR<1)为交叉概率；RN∈[1,2,3]；

(d)生成新个体

将试验矢量U_i与目标矢量P_i ^t进行比较，若U_i优于P_i ^t，则用U_i代替P_i ^t，进入下一代：

按照(b)、(c)、(d)的运算过程对每个个体进行传代，得到第t+1代种群，计算各个体所对应的适应度函数G(P_i ^t+1)，若min{G(P_i ^t+1)}＜G(P_best)，则由min{G(P_i ^t+1)}对应的个体代替原有P_best，否则P_best维持不变；

(e)迭代和优化PID控制参数

按照上述步骤进行多次迭代，若P_best连续n代没有更新，则适应度函数已经收敛，利用下式更新PID控制参数：

k_p(k)＝(1+P_best1)k_p(k-1)

k_i(k)＝(1+P_best1)k_i(k-1)

k_d(k)＝(1+P_best1)k_d(k-1)。

7.根据权利要求1所述的方法，其特征在于，所述步骤5是在步骤4更新PID控制参数后，按下式更新底物流加速率F(k)：

F(k)＝F(k-1)+k_p(k)[e(k)-e(k-1)]+k_i(k)e(k)+k_d(k)[e(k)+e(k-2)-2e(k-1)]。