基于动态优化的青霉素Fed-batch 发酵过程分段操作策略

伍子奇，吴思辰，翟 持

昆明理工大学化学工程学院，云南昆明 650500

随着我国社会经济的快速发展，人们愈发关注健康、资源、环境等与民生息息相关的要素。近年来，“双碳”战略为绿色生物制造的发展注入强劲的动力。如图1 所示，生物制造以生物质为原料，通过生化反应实现物质的高效、定向转化，能摆脱对化石资源的过度依赖，减少有毒有害化学试剂的使用，是实现可持续发展最有希望的技术之一[1]。区别于传统的化学反应过程，生物转化反应条件相对温和，推动力低[2]；细胞膜将微生物与培养基分隔，形成一个与外界有物质、能量交换的相对独立子系统，可调控手段少[3]；微生物还有趋利避害、自我调节等特征。事实上，微生物的自主调节作用在一定程度上制约了外界调控的有效性。然而，随着新型分离、检测技术及精准控制技术的发展，尤其是在《中国制造2025》等国家战略推动下，我国生物制造行业正处在转型的关键时期，并向先进化、智能化控制方向发展。本工作调研发酵过程的操作策略，并以青霉素Fed-batch 发酵为例，探讨其最优给料策略，为Fed-batch 发酵模式的操控优化提供必要的理论及模拟支撑。

图1 生物制造的物质转化Fig.1 Schematic diagram of substance transformation in bio-products

生物制造的核心环节是借助微生物、动物细胞的生命活动来制备菌体，或生产直接、间接的代谢产物[4]。常见的发酵生产模式见表1，其中，微生物往往以自由悬浮状态或固定化细胞/酶形式存在[5]。需要强调的是，工业发酵是典型的开放系统：生物生长、繁殖、代谢及衰亡过程需要通过“新陈代谢”维持生命体征（物质的输入-输出，新老细胞的更替等）[6]。然而，当前仅有少量生物过程（单细胞的蛋白质生产、啤酒发酵过程、市政污水处理等）采用连续生产模式，而大多数发酵过程依然采用间歇生产模式。连续发酵主要的问题是工程上很难维持无（杂）菌环境，使其不受噬菌体攻击，或长时间生产不发生变异。因此，工业上的稳定发酵过程效果与间歇、半连续过程相对较差，且由于生物过程机理复杂，很难对连续生产效果差的原因进行溯源。相较于连续发酵，间歇发酵为封闭体系内的批次生产过程，例如Saccharomyces cerevisiae发酵生产乙醇的过程[14]，在一个发酵批次为60～72 h，几乎不采取任何操控干预，仅依靠微生物自身的生长-繁殖-代谢进行乙醇的生产，而在单个发酵周期内，只有16%～20%的时间处于乙醇的生产期。并且不同生产批次间还需进行加料、接种、产品分离及洗涤除菌等附加操作（如图2），这导致间歇发酵生产效率较低。因此，工业上通过反应釜的串联（表1，Ⅱ）或多反应釜优化排产来提高生产效率。Fed-batch 生产模式（表1，Ⅲ）介于间歇与连续发酵之间，其通过连续或脉冲的方式对发酵液进行基质补给，并在发酵终止时收获产品。调节进料物流可对反应物浓度进行精准调控，使Fed-batch操作终止时获得最优效益。Fed-batch 操作适用于复杂反应网络（抑制性酶转化、自催化、平行反应等）对特定反应的转化率或选择性进行优化的情况。微生物增殖过程是典型的自催化反应，因此，多数工业发酵过程采用Fed-batch 生产模式[15-17]。对Fed-batch 生产而言，核心问题是补给物的选取和给料方式的确定。

表1 生物反应器的主要结构及操作策略Table 1 Main structures and operation strategies of bioreactors

图2 典型间歇发酵过程流程示意Fig.2Schematic diagram of a typical batch fermentation process

Fed-batch 生产的常见操作方式是：在发酵初期调控发酵液内限制性基质，使微生物生长速率最大化，直至微生物充满发酵液后降低或关闭基质补给，并在发酵终止时收获产品。因此，可将Fed-batch生产视为一个动态优化问题，即针对特定的微生物生长、代谢规律，运用最优化的方法获得确定性系统的给料最优曲线，使发酵终止时获得最大收益（生产效率最大、产品产量最大、基质损失最小等）。

青霉素生产过程中，碳源提供过多会导致生物质生长强化而青霉素产品生产受限，而碳源提供过少会抑制生物质的生长，进而对青霉素的生长产生消极影响。Fed-batch 可以适时提供所需碳源。因此，本工作以青霉素发酵过程的Fed-batch 生产为例，从最优化理论及数值优化计算角度进行研究。首先讨论青霉素发酵过程动力学，并构建最优化问题，再从经典动态优化理论出发，探讨所涉bang-singular-bang 优化操控问题，并基于配置法对青霉素最优补碳操作进行数值模拟，最后，结合理论分析和数值模拟对青霉素Fed-batch 生产过程进行分析讨论。

1 基于Fed-batch 的青霉素优化模型

青霉素发酵过程中微生物将葡萄糖（碳源、能量来源）、有机氮源等用于细胞的繁殖及青霉素生产，根据过程动力学建立生物质浓度（X），产物浓度（P），基质浓度（S）和有效反应体积（V）的动态守恒方程，如式（1）所示：

微生物生长动力学（μ）如式（2）所示：

青霉素的反应动力学（ρ）需考虑原料的代谢抑制，具体方程式如式（3）所示：

原料（葡萄糖）的消耗动力学（σ）主要用于微生物生长、产品生成及维持微生物的代谢，具体方程式如式（4）所示：

假设葡萄糖是唯一的限制性基质，青霉素发酵的Fed-batch 生产最优化是获得最佳的进料曲线，以使青霉素产量最大、发酵时间最短或者基质的产品转化率最高等为目标。因此，可构建动态优化问题，如式（5）所示：

2 动态优化理论与分析

根据Pontryagin 最小值原理（Pontryagin’s Minimum Principle，PMP），上述优化问题可写成[18]：

PMP 最优化条件如式（7）所示：

针对终值时间（tf）为自由变量的情况，最优化还要满足transervality 条件，如式（8）所示：

由于状态方程的初值已知，并且伴随变量λ，ω，v终值条件已知，因此，方程（6）所述优化为两边界约束的规划问题。

两边界约束的求解涉及迭代（如打靶法），可将式（6）转化为Hamilton-Jacobi-Bellman（HJB）方程[19]。定义函数V(x，t)，得：

对比方程（6）和（9），当Lagrange 乘数子λ满足：

PMP 与HJB 问题一致，此时有：

因此，原Hamiltonian 函数H的两边界问题转化为求解偏微分方程的问题。

当仅考虑唯一操控输入时，优化的必要条件为HU=0，将其展开如式（12）所示：

影响HU的因素有过程系统自身及约束条件。当方程（12）右侧第一项不为零时，需满足μ≠0，即，确定最优U时过程约束需要被激活。当系统唯一的约束条件是操控U的上下界，就会形成开关操控。

另外，由方程（1）可知，青霉素发酵过程是关于U仿射的，即，过程变量为x=(X，P，S，V)T时，动态模型可写成式（13）：

最优化目标是确定基质输入速率U，其受到上下界的约束（2）。由Pontryagin 最小值原理可将原问题转化为最小化汉密尔顿函数如式（14）：

当U是线性的，最小化汉密尔顿函数变得十分简单：当λTb为正，可以选择输入为操控下限；当λTb为负，可以选择输入为操控上限。到目前为止，这是一个bang-bang 操控的问题，但是当λTb在一定时间间隔内恒为零的时候，出现奇异点，此时无法直接确定U的具体输入。因此，连续-间歇青霉素发酵过程是一个bang-singular-bang 操控的问题，通过判据λTb可描述为式（15）所示分段函数：

再有，沿最优路径的H*也是未知的，因此需要如下关系来确定：

方程（17）是两点边界条件问题，但是需要确定三种输入类型（最大，最小及奇异值）的顺序才能确定边界，并根据不同的输入顺序确定相关计算方法。但是，更为广泛采用的计算方法是将原问题离散化，进而数值求解非线性规划（NLP）问题。

3 最优输入曲线的数值求解策略

基于数值方法的动态优化有直接法和间接法，其中直接法通常有两种计算策略：1）将操控轨迹设置为特定的基准函数的线性组合，如三角函数形式，通过所涉参数的数学规划获得最优操控轨迹，该方法维持过程ODE/DAE 原形式；2）利用多项式，如Lagrange 多项式，同时离散化过程及操控轨迹，将原问题转化为NLP 问题。本研究使用后者进行后续计算机分析。

通过正交配置可将原问题中的状态方程离散化并转化为NLP 问题。一般而言，配置点可选为[0,1]区间Legendre 多项式的零值处获得，而积分、微分项可通过Lagrange 插值多项式给出。考虑到发酵终止时间（tf）是青霉素Fed-batch 操作的关键参数，由于过程ODE/DAE 是通过正交配置离散化为NE段，对tf进行数学规划时，可对各个配置单元进行时间归一化：

如图3 所示，经归一化后的系统在(x,ζ)空间的正交配置离散化可描述为对i∈(1，...，NE)段有限元的离散化，进而转化为状态方程在i处的初值问题，其积分时间为t∈(ζi，ζi+1)，那么，Lagrange 插值在ζi≤t≤ζi+1的近似可写成如下方程：

图3 状态曲线、操作曲线在有限元的配置情况Fig.3 Configuration of state curve and operation curve in finite element

式中：i=1，2…，NE；k=0，1…，j-1，j，而k≠j。同时，xKx(t)是Kx+1阶多项式，UKu(t)是Ku阶分段连续多项式。在近似状态变量x的时候考虑到其每个元素均为初值问题，在近似的节点处满足如下关系：

假如K=Kx=Ku=2，而对有限元的选取通常是等间距的，其余差方程可写成：

式中，i=1，2…，NE；j=1，2…，NE；ϕj(τk)=ϕjdt/dτk，tik=ζi+ζik。虽然要求状态变量在元素节点(ζi，i=2，…，NE)处连续，因此在节点处的状态变量满足式（23）所述连接关系：

此时的边界条件可写成：

需要指出的是，状态变量在元素节点处的连续性是出于实际过程考虑，而操作变量在元素节点处不要求连续性体现上述讲到的出现奇异值的问题。综上，原问题被转化为NLP 问题，可以通过序贯二次规划[20]或者SNOPT 算法[21]求解。

4 青霉素Fed-batch 操作的分析

当考虑在反应终止时的青霉素产量最大化，即，Φ=-P(tf)V(tf)，在正交配置法进行求解过程中，由于状态、操控的配置点、离散化的有限元个数、初始操控输入及优化算法选择的差异，可能导致寻优获得的操控轨迹差别较大。本工作在前人工作的基础上[22]，进一步讨论分段恒值问题，首先给出输入无约束的最优结果。如图4 所示，当选择所述操控曲线，在tf为136.81 h 时，青霉素最优化产量为87.32 g，此时对应的状态终值(X，P，S，V)T=(27.44，8.73，0.004 6，10.00)Tg/L。其中，Kx=Ku=3，NE=10，初始基质输入浓度为15 g/L，NLP 求解算法为active-set 算法[22]。

图4 青霉素Fed-batch 发酵过程的最优操控及相应的状态变化Fig.4 Optimal control of penicillin Fed-batch fermentation process and corresponding state changes

考虑图4 所述操控曲线U变化复杂，可适当调节优化参数来获得操作上易于实现的操控输入。由方程（20）可知，操控输入U在离散节点处具有连续性，但Lagrange 插值比相应的状态低一阶，即，在离散元节点处的U不具备光滑性质。考虑U连续变化实现起来复杂，可选取Ku=1 时，其输入为零阶保持。另外，可通过降低输入曲线的配置点个数来降低操作频率。图5 为tf=136.81h 时对操控U取不同离散点个数的输入轨迹及状态终值。模拟结果显示，恰当的离散点个数不仅可以降低操作的复杂程度，还能优化目标值（本例可选取离散点个数为8）。

图5 青霉素Fed-batch 发酵过程的最优操控及相应的状态变化情况Fig.5 Optimal control of penicillin Fed-batch fermentation process and corresponding state changes

上述模拟计算可知，V达到上限（10 L）是导致该优化操控终止的条件，即在给定基质总量的情况下，研究最佳的给料策略使青霉素产量最大化。另一方面，需要优化基质的产品转化率，并且减少发酵终止时料液内基质的浓度，以减少后续产品、原料分离成本。当规划目标为Φ=-P(tf)V(tf)/S(tf)时，选取tf=136.81 h，对操控U取不同离散点个数的输入轨迹及状态终值。模拟结果显示，更新的优化目标不仅可以大幅降低发酵终止时基质的浓度，还有可能提升青霉素的产量。对比图4 及图5（离散点个数为8）的结果可知，后者的产品产量提升0.58%。对比图5、图6 在离散点取8 时的规划结果，发现Φ=-P(tf)V(tf)/S(tf)时，图6 对应的操控输入可以提升青霉素产量（10.7%），并且产品中基质降至微量，可以减少后续产品、原料的分离过程。

图6 青霉素Fed-batch 发酵过程考虑目标Φ=-P(tf)V(tf)/S(tf)的优化操控输入轨迹Fig.6 Optimal control input trajectory considering the target in penicillin Fed-batch fermentation process

5 结 论

经分析与优化计算可知，基于有限元的正交配置法可以有效处理动态优化过程中的bang-singular-bang 问题。操控输入在有限元内采用零阶保持方式时，将Fed-batch 反应器的操控域进行脉冲切割，一方面，获得的最优操作轨迹较容易面向工业推广，另一方面，输入脉冲变化的操控轨线可以获得优于连续输入的Fed-batch 优化方案。另外，通过构造考虑原料消耗成本的目标函数，可以降低基质流失，同时大幅度提升青霉素的产量，并且计算得到的分段操作方案可使产品中基质降至微量，进而减轻后续产品、原料的分离负担。该青霉素优化计算方案可推广至一般Fed-batch 过程，辅助工业过程建立最优给料策略。

符号说明

a—— 仿射非线性系统状态项

b—— 仿射非线性系统的输入系数项

c—— 两点边界问题中下边界系数

f—— 青霉素发酵物质守恒方程

g—— 等式约束函数

H—— 汉密尔顿算子

H*—— 最优汉密尔顿算子

HU—— 汉密尔顿算子在控制放心的导数

J—— 表示最优化目标函数（最小目标）

Kdeg—— 产品退化的速率常数，其值为0.01 h-1

Kin—— 基质抑制常数，其值为0.1 g/L

KS—— 基质的饱和常数，其值为0.006 g/L

KP—— 产物的饱和常数，其值为0.000 1 g/L

Km—— 原料的饱和常数，其值为0.000 1 g/L

Ku—— 配置法中输入的配置点

Kx—— 配置法中状态的配置点

ms—— 原料的最大比消耗速率，其值为0.029 h-1

NE—— 配置法积分间隔内配置点个数

P—— 青霉素产品浓度，g/L

r—— 剩余函数

S—— 基质浓度，g/L

SF—— 原料的给料浓度，其值为500 g/L

T—— 终值约束

t—— 时间，h

tf—— 间歇发酵的终止时间，其值为136.81 h

U—— 原料的给料量，g/h

V—— 有效反应体积，L

v—— 终止条件的乘算子

x—— 表示发酵过程状态变量

X—— 生物质浓度，g/L

YP/S—— 原料的产物转化系数，其值为1.2 g/g

YX/S—— 原料的生物转化系数，g/g

σ—— 原料（葡萄糖）的消耗动力学

ζ—— 积分时间范围

θ—— 归一化时间，θ∈（0，1）

λ—— 系统约束的Lagrange 乘数子

μ—— 表示微生物生长动力学

μmax—— 微生物的最大比生长速率，其值为0.11 h-1

ρ—— 表示青霉素的反应动力学

ρmax—— 产物的最大比生长速率，其值为0.055 h-1

τ—— 配置点时间间隔

Φ—— 动态优化的目标函数

ϕ—— Lagrange 插值函数

ω—— 过程等式的Lagrange 乘数子

上标

l—— 状态变量及调控变量的下界限

u—— 状态变量及调控变量的上界限

T —— 向量的转置

下标

i—— 第i个有限元间隔

j—— 第j个配置点

k—— 计量数