来源：A multi-model predictive control method for the Pichia pastoris fermentation process based on relative error weighting algorithm（Alexandria Engineering Journal）| 编译：DNA Lab Space | 原文许可：elsevier-subscription
导读
毕赤酵母发酵生产蛋白酶K是一个高度非线性、时变且强耦合的生化过程，传统PID等控制方法难以满足实际需求。本文提出一种基于相对误差加权算法的多模型预测控制策略——先用模糊C均值聚类（FCM）将先验数据划分为多个训练样本集，再对每个样本簇用最小二乘支持向量机（LSSVM）结合改进粒子群算法（IPSO）建立子预测模型，最后通过相对误差加权算法融合各子模型输出实现全局优化控制。仿真结果表明，该方法在宽范围设定值变化下表现出更快的响应速度和更好的跟踪精度。
研究背景
蛋白酶K属于丝氨酸蛋白酶，具有高酶活性和广泛的底物特异性，能优先切割疏水性氨基酸、含硫氨基酸、芳香族氨基酸及C端相邻的酯键和肽键，常被用于降解蛋白质产生短肽。凭借这一特性，蛋白酶K在核酸纯化、丝绸、医药、食品、酿造等领域有重要应用。毕赤酵母是一种甲醇营养型酵母，能以甲醇为唯一碳源和能源，是目前应用最广泛的外源蛋白表达系统。与其他表达系统相比，毕赤酵母在表达产物的加工、外分泌、翻译后修饰及糖基化修饰方面具有明显优势。因此，利用毕赤酵母发酵生产蛋白酶K，是实现其高效表达突破的关键步骤。为充分发挥毕赤酵母表达系统的优势、最大化蛋白酶K的表达水平，必须对发酵过程进行动态调控与优化，在最佳培养条件下精确控制重组蛋白酶K的表达。然而，毕赤酵母发酵过程具有高度非线性、时变性和强耦合特性，经典PID等传统控制方法无法满足实际控制需求，亟需更先进的控制策略。
研究方法
多模型预测控制的基本思想是将被控对象的非线性空间划分为若干子空间，在每个子空间内建立局部模型，并针对每个局部模型设计相应的预测控制器。通过多个模型实时逼近被控对象的动态特性，最终通过切换到最优子控制器或对各子控制器输出进行加权求和，获得最优控制效果。本文提出一种基于加权算法的多模型预测控制策略：先利用模糊C均值聚类算法（FCM）将先验数据划分为m个训练样本集（样本聚类）；再对每个样本簇采用最小二乘支持向量机（LSSVM）和改进粒子群优化算法（IPSO）获得相应的预测模型；随后针对m个局部预测模型分别设计对应的预测控制器（MPC1~MPCm）；最后在每个采样时刻计算对象输出与各子预测模型输出之间的偏差，基于多模型相对误差加权算法建立预测模型的控制策略，以实现对被控对象的优化控制。该方法通过m个局部模型改善系统的瞬态性能，使被控量能快速跟踪给定值；同时利用相对误差加权算法构建实时预测模型的控制策略，增强模型的自适应能力，使其更准确地描述非线性系统的实际状态。
基于FCM-IPSO-LSSVM的多模型建模
工业过程通常按一定顺序或在特定工况下进行，工况和工作点相对确定。因此，可根据不同工况将非线性系统的全局运行区间划分为若干典型区间。目前划分非线性系统全局运行区间的方法中，基于数据或先验知识的划分方法受到广泛关注。数据聚类分析是一种基于数据的划分方法，常用的有K-means聚类、层次聚类、高斯混合密度聚类、模糊C均值聚类算法等。本文通过被控对象的先验采样点，采用FCM算法将整个工作区间划分为不同子空间，再针对划分后的样本子空间，利用LSSVM和IPSO构建相应的子预测模型。
#### 模糊C均值聚类
FCM算法的基本思想是在硬聚类算法（HCM）基础上引入不同类别样本的隶属度函数矩阵和模糊系数。与HCM算法相比，FCM在数据分类、聚类点计算、目标函数等方面进行了调整。给定数据集X = {x₁, x₂, ..., xₙ}，n为样本数，FCM算法将数据集X划分为C（2 ≤ C ≤ n）类。FCM聚类算法的目标函数如下：
min Jₘ(U,V) = Σⱼ₌₁ⁿ Σᵢ₌₁ᶜ uᵢⱼᵐ dᵢⱼ²
约束条件为：Σₖ₌₁ᶜ uᵢₖ = 1，uᵢₖ ⊂ [0,1]，Σᵢ₌₁ⁿ uᵢₖ > 0
其中C为聚类数，uᵢⱼ为样本xⱼ对第i类的隶属度，Vᵢ为第i类的聚类中心，dᵢⱼ = ||xⱼ - vᵢ||为样本xⱼ与中心vᵢ之间的欧氏距离。为使目标函数min Jₘ(U,V)最小化，Vᵢ和隶属度矩阵U可通过下式计算：
Vᵢ = Σⱼ₌₁ⁿ uᵢⱼᵐ xⱼ / Σⱼ₌₁ⁿ uᵢⱼᵐ，i = 1, 2, ..., C
uᵢⱼ = 1/Σₖ₌₁ᶜ (dᵢⱼ/dᵢₖ)^(2/(m-1))（dᵢⱼ ≠ 0时）；当dᵢⱼ = 0且j = k时uᵢⱼ = 1；当dᵢⱼ = 0且j ≠ k时uᵢⱼ = 0。
FCM算法的具体流程如下：
步骤1：输入聚类数C、模糊加权参数m和迭代停止条件δ；
步骤2：初始化聚类中心Vᵢ⁰（i = 1, 2, ..., C）；
步骤3：用式(3)计算uᵢⱼ（i = 1, 2, ..., C；j = 1, 2, ..., n）；
步骤4：用式(2)计算Vᵢˡ（i = 1, 2, ..., C）；
步骤5：若‖Vˡ - V⁰‖ ≤ δ，则停止迭代并跳至步骤6；否则令Vˡ = V⁰，跳至步骤3；
步骤6：输出聚类结果(V, U)。
【数据】算法：FCM；输入参数：聚类数C、模糊加权参数m、迭代停止条件δ；迭代终止判据：‖Vˡ - V⁰‖ ≤ δ
#### 最小二乘支持向量机
支持向量机（SVM）是Vapnik提出的机器学习方法，具有泛化能力强、克服维度灾难等优点。Suykens和Vandewalle在SVM基础上提出最小二乘支持向量机（LSSVM），以系统误差平方和作为损失函数，用等式约束替代不等式约束。LSSVM优化问题如下：
min J(w,e) = ½wᵀw + ½γΣₖ₌₁ᴺ eₖ²
约束条件：yₖ = wᵀφ(x) + b + eₖ，k = 1, 2, ..., N
其中φ(·): Rⁿ → Rⁿᴴ是将输入数据映射到高维特征空间的函数，w为权值向量，eₖ为误差变量，b为偏差量，γ为正则化参数，用于控制对超出误差样本的惩罚程度。设核函数K(xᵢ, xⱼ) = φ(xᵢ)ᵀφ(xⱼ)，LSSVM最终可表示为：
f(x) = Σᵢ₌₁ᴺ aᵢK(x, xᵢ) + b
其中aᵢ ∈ R为拉格朗日乘子。经对多类函数比较分析后，选择RBF函数作为LSSVM的核函数，其核宽度σ可表示为：
K(xᵢ, xⱼ) = exp(-(xᵢ - xⱼ)² / 2σ²)
参数a和b的表达式如下：
[0 Qᵀ; Q Ω + I/μ][b; a] = [0; 1ᵥ]
其中Q = [y₁, ..., yₙ]ᵀ；a = [a₁, ..., aₙ]ᵀ；1ᵥ = [1, ..., 1]ᵀ；Ω为核矩阵，表达式为Ωᵢ,ⱼ = yᵢyⱼφ(xᵢ)ᵀφ(xⱼ) = yᵢyⱼK(xᵢ, xⱼ)。
LSSVM的预测能力主要取决于正则化参数γ和核宽度σ。其中γ影响模型的拟合精度和泛化能力，σ直接决定模型的计算量和执行效率。目前常用的参数选择方法有网络搜索算法和遗传算法——前者计算规模大、实时性差，后者容易陷入局部极小值。本文提出一种改进粒子群优化（IPSO）算法来优化LSSVM模型的参数(γ, σ)，获得更优的参数优化质量。
【数据】算法：LSSVM；核函数：RBF；待优化参数：正则化参数γ、核宽度σ
#### 改进粒子群优化算法
PSO是一种基于仿生学的算法，模拟鸟类觅食的随机行为。该算法的核心是利用群体中个体的信息共享，使整个群体在问题求解空间中从无序到有序地进化，从而获得问题的最优解。PSO算法因其定义简单、实现方便、计算速度快、全局搜索能力好，已被广泛应用于函数优化、神经网络训练、模糊系统控制等遗传算法应用领域。国内外学者也在此基础上进行了多种改进（如QPSO、AQPSO、GQPSO、WQPSO），并将其应用于众多领域。PSO通过跟踪和调整群体中每个粒子的位置和速度来实现整个群体的优化效果。设目标搜索空间的维度为...
【数据】算法：IPSO（改进粒子群优化）；基础：PSO仿生算法
研究结果
仿真数据来自镇江扬中威克特生物工程设备有限公司微生物发酵实验室，发酵实验设备采用RTY-C-100L型。以毕赤酵母发酵为对象，选用上海市农业科学院畜牧兽医研究所重点实验室构建的毕赤酵母KM71 H菌株、MutS型作为生产菌株，表达载体和外源基因分别为pPICZαA和IFNαcDNA。仿真按流程图进行。
深入分析发酵过程后，选取发酵液温度T、溶解氧浓度DO、发酵液pH、空气流量q、搅拌转速v和发酵罐压力p作为输入变量。采用所提方法控制发酵罐中的细胞浓度X和蛋白酶K浓度P。发酵过程中输入变量的初始值设置如下：发酵罐压力控制在0.04 MPa，电机搅拌转速约250 r/min，发酵温度控制在约28 ± 0.5 ℃，空气流量控制在1000–1200 L/h范围内，氧含量维持在35%–45%之间，pH为7.3。
根据本文方法建立局部预测模型FCM-IPSO-LSSVMₘ，并设计无约束GPC控制器来控制细胞浓度X和蛋白酶K浓度P，预测时域nₚ和控制时域n꜀分别设为7和4。细胞浓度X和蛋白酶K浓度P的输出预测偏差加权因子分别为1.2和1.5。控制增量加权系数为0.5；反馈校正系数为1。恒定扰动Fbf = 0.35。
其他参数设置如下：(1) 非线性空间划分为4个子空间，即通过FCM算法将先验样本数据划分为四个训练样本集（样本簇），C = 4，这是考虑到毕赤酵母KM71 H菌株产蛋白酶K的补料分批发酵特性（其生长曲线包括延迟期、对数生长期、平衡期和衰老死亡期四个阶段）。(2) 采用相对误差加权算法进行仿真时，式(13)和式(14)中V = 1.2，δ = 0.7，初始加权因子假设为均匀分布，即wⱼ(0) = wⱼ′(0) = 1/N。
当细胞浓度X从0 → 0.8 → 1.8 → 1.6变化、蛋白酶K浓度P从0 → 40 → 100变化时，单模型预测控制方法和本文所提控制方法的X和P跟踪响应特征曲线分别如图8和图9所示。
表1和表2比较了两种控制方法的性能指标。其中，表1采用均方误差指标（RMSE）来衡量各输出与设定值的偏差，反映输出变量的控制效果。
从图8、图9中仿真结果与多模型预测控制方法的对比以及表2可以看出，后者能加快系统响应速度，并更好地抑制超调。从表1还可看出，本文提出的控制方法提高了对设定值的跟踪精度（尤其是当细胞浓度设定值从0.8变为1.8、再从1.8变为1.6时，单模型预测控制方法在接近设定值时表现出较大且多次的振荡）。本文所提控制方法针对各预测模型进行整体加权，即在预测控制的滚动优化中，优化函数中的预测输出由各单一模型的预测输出加权得到，据此计算控制增量，使控制解更具全局最优性和灵活性。因此，对于毕赤酵母产蛋白酶K补料分批发酵这类强耦合、多输出的非线性对象，所提控制方法能大幅改善对输出变量（对输入变化敏感且局部模型切换频繁）的控制效果。
【数据】设备型号：RTY-C-100L；菌株：毕赤酵母KM71 H（MutS型）；表达载体：pPICZαA；外源基因：IFNαcDNA；发酵罐压力：0.04 MPa；搅拌转速：250 r/min；发酵温度：28 ± 0.5 ℃；空气流量：1000–1200 L/h；氧含量：35%–45%；pH：7.3；预测时域nₚ：7；控制时域n꜀：4；输出预测偏差加权因子：1.2（细胞浓度X）、1.5（蛋白酶K浓度P）；控制增量加权系数：0.5；反馈校正系数：1；恒定扰动Fbf：0.35；聚类数C：4；V = 1.2；δ = 0.7；初始加权因子：wⱼ(0) = wⱼ′(0) = 1/N；X变化范围：0 → 0.8 → 1.8 → 1.6；P变化范围：0 → 40 → 100；性能指标：RMSE（均方误差）
讨论与解读
本文提出一种基于加权算法的多模型预测控制策略。首先利用FCM算法将先验数据划分为多个训练样本集，针对每个样本聚类用LSSVM和IPSO算法获得相应的子预测模型；随后为每个子预测模型设计对应的预测控制器；最后计算每个采样时刻目标输出与各子预测模型输出的偏差，基于多模型相对误差加权算法建立预测模型的控制策略，实现对被控对象的优化控制。该方法通过多个局部模型改善系统的瞬态性能，使被控量能快速跟踪给定值；同时，相对误差加权算法构建的实时预测模型控制策略提高了模型的自适应能力，使其能更精确地描述非线性系统的实际状态。本文还指出，为分析和提高算法计算效率，需对算法的时间复杂度、内存复杂度和计算复杂度进行分析，但为避免一般读者偏离本研究工作的主线，本文仅对算法进行了...（原文未完整详述）。
从合成生物学应用角度看，该研究将计算智能方法引入发酵过程控制，抓住了毕赤酵母表达系统“模型非线性强、工况多变”的痛点——用FCM分而治之、用LSSVM小样本建模、用IPSO自动寻参，再用加权融合替代硬切换，逻辑链条完整。其价值在于：不依赖精确机理模型即可实现对强耦合生物过程的满意控制，这对工业规模发酵放大具有现实意义。局限也很明显——仿真数据来自单一菌株和特定设备，控制器的实际鲁棒性仍需在更多菌株、更多扰动场景下检验；算法复杂度分析未充分展开，在线计算负担有待评估。
参考来源
Wang B. A multi-model predictive control method for the Pichia pastoris fermentation process based on relative error weighting algorithm[J]. Alexandria Engineering Journal, 2022. DOI: 10.1016/j.aej.2022.03.004
DOI: 10.1016/j.aej.2022.03.004