基于梯度提升回归树的发酵过程软测量研究∗

黄继炜 潘 丰

（江南大学轻工过程先进控制教育部重点实验室 无锡 214122）

1 引言

发酵工程［1～2］作为生物工程一个极其重要的分支，与人们的生产生活有着密不可分的关系。然而，要实现发酵过程的优化控制［3～5］，必须要知晓发酵过程中的各种关键要素。其中温度、pH 值、压力、二氧化碳浓度、溶解氧（Dissolved Oxygen，DO）等变量已经可以实现高精度的实时测量，但是生物质浓度、基质浓度、产物浓度等关键变量，现阶段缺乏经济、可靠的传感器来进行高精度的实时测量，离线分析又有很大的时间滞后。

为了解决上述问题，人们开始通过对生物发酵过程建模［6～7］，从而对难以实时测量的关键变量进行预测，通过参考预测值来实现进一步的优化控制。其中，机器学习［8～9］的方法越来越受到关注，其常见的方法有BP 神经网络（Back Propagation Neural Network，BPNN）、决策树（Decision Tree，DT）、支持向量回归（Support Vector Regression，SVR）、Ada-Boost、随机森林（Random Forest，RF）等。以上方法在针对某些特定应用场景时都有着不错的表现，但进一步推广时却存在一定的问题。BPNN［10］学习速度慢、容易陷入局部极小值，DT［11］对缺失数据处理比价困难、容易出现过度拟合，SVR［12］对大规模训练样本难以实施、解决多分类问题存在困难，Ada-Boost［13］训练比较耗时、数据不平衡会导致分类精度下降，RF［14］处理特征较少的数据时不能产生很好的分类，对某些特定噪声的数据进行建模时可能出现过度拟合。针对上述问题，本文提出了一种采用果蝇算法优化梯度提升回归树（Gradient Boost Regression Tree Optimized by Fruit Fly Optimization Algorithm，FOA-GBRT）的发酵过程软测量方法，利用梯度提升回归树（Gradient Boost Regression Tree，GBRT）进行建模与预测，同时利用果蝇优化算法（Fruit Fly Optimization Algorithm，FOA）对回归树的关键参数进行寻优，提高其对样本的拟合精度。本文利用Pensim 青霉素仿真平台生产的青霉素发酵数据进行实验验证，同时与BPNN、DT、SVM、AdaBoost、RF方法进行对比实验，从而验证本文方法的有效性。

2 算法原理

2.1 梯度提升回归树算法

梯度提升回归树［15］由多棵决策树组成，其核心就在于，每棵树是从先前所有树的残差中来学习。利用的是当前模型中损失函数的负梯度值，作为提升树算法的残差的近似值（伪残差），进而拟合一棵回归树，这就是GBRT的算法原理。

对于给定的输入：训练数据集T=(x1,y1),(x2,y2),…,(xn,yn)，损失函数L(y,f(x))，GBRT 算法的输出结果为fˆ(x)，其流程如下。

Step1：首先初始化，估计一个使损失函数极小化的常数值，即一棵只有根节点的树：

Step2：迭代开始，建立M棵梯度提升树：

Step2.1：计算损失函数的负梯度在当前模型的值，并将它作为残差的估计值：

Step2.2：对于rmi拟合一棵回归树，得到第m 棵树的叶节点区域：

Step2.3：利用线性搜索估计叶节点区域的值，使损失函数极小化：

2.2 果蝇优化算法

果蝇优化算法［16］采用了基于果蝇群体协作的机制进行寻优操作，图1 描述了果蝇群体进行食物搜索的基本过程。

图1 果蝇觅食行为

标准的FOA算法流程包括以下步骤。

Step1：初始化

给定果蝇种群规模（popsize），算法迭代次数（maxgen），果蝇群体位置范围（LR）和果蝇单次飞行范围（FR）等参数值。群体中每个果蝇个体的位置信息根据其对于的二维坐标（X，Y）来确定，其初始值如下：

Step2：嗅觉搜索过程

Step2.1：在群体中的各只果蝇通过其嗅觉进行搜索时，给定其一个随机的飞行距离和方向。第i个果蝇个体的新位置如下：

Step2.2：由于食物味道的位置信息是未知的，所以先计算第i只果蝇到原点的距离Di：

然后计算其味道浓度判定值Si：

Step2.3：计算第i只果蝇的味道浓度值Smelli：

fitness 是味道浓度判定函数，在通过FOA 进行实际问题求解时，它表示目标函数或适应度函数。

Step2.4：比较得出当前群体中味道浓度值最佳的果蝇个体，记录其位置信息和味道浓度值：

Step3：视觉搜索过程

记录味道浓度的最佳值和与之对应的果蝇的位置信息，群体中的其他果蝇通过视觉搜索向此位置移动：

Step4：重复Step2 和Step3，直至达到算法迭代次数maxgen。

3 FOA-GBRT软测量方法

3.1 FOA算法改进

考虑到标准的FOA 存在一定的局限性，比如对最优解为负数的问题无能为力，对复杂、高维和非线性问题求解能力较弱，容易陷入局部极值点。针对以上问题，本文采用以下两种方法来进行优化。

1）在果蝇单次飞行范围的计算中引入权重系数w，使得果蝇的搜索半径成为一个自适应的值，能按指数形式随着迭代次数的增加而下降［17］：

果蝇位置初始化：

果蝇飞行方式：

其中，n 为搜索范围系数，w0为初始权重，a 为权重系数，gen 为当前进化次数。此举可以使得算法的局部搜索能力和全局搜过能力取得平衡，使得算法的收敛速度获得提升。

2）摒弃二维坐标，直接使用一维坐标来生成味道浓度判定值的候选解［18］：

此举的好处在于能够使算法充分且均匀地在问题定义域内搜索，有助于获得最优解不在零点以及复杂问题的全局最优解。

3.2 果蝇算法优化梯度提升回归树参数

GBRT 运行中需要设置许多参数，例如学习速率（learning rate）、最大迭代次数（n estimators）、子采样（subsample）、决策树最大深度（max depth）、内部节点再划分所需最小样本数（min samples split）、叶子节点最少样本数（min samples leaf）等。不同的参数组合会使GBRT 对于同一样本的拟合效果有所差异，而只依靠经验估计又难以设定出最佳的参数组合。所以本文采用果蝇算法对GBRT 的关键参数进行寻优，从而找到GBRT 对当前样本的拟合度最高的参数组合，从而提升GBRT 的预测精度。

本文通过改进后的FOA 对GBRT 的六组关键参数进行寻优，包括learning rate、n estimators、max depth、min samples split、min samples leaf、subsample。寻优整体流程如图2所示。

图2 参数寻优流程图

种群中每个果蝇个体分别携带有一组待优化参数值的组合，将GBRT 的参数设置为该组值进行建模预测，得到的拟合优度作为对应果蝇的味道浓度值，然后进行嗅觉和视觉搜索，不断重复这个过程直到达到迭代次数，最终味道浓度值最大的果蝇所携带的参数值即为优化后的参数组合。

寻优过程采用7200组数据作为训练集，每组5个特征值，800 组数据作为测试集，来验证参数优化后的GBRT 预测效果是否有所提升，评价指标包括拟合优度（R2）和均方根误差（RMSE）。

R2即判定系数，是对估计的回归方程拟合优度的度量，其值越接近1，表示模型的拟合优度越高。假设一组数据集包括y1，y2，…，yn共n 个观察值，相对应的模型预测值分别为f1，f2，…，fn，则判定系数定义为

RMSE 是预测值与观察值偏差的平方和与观察次数比值的平方根，用来衡量预测值同观察值之间的偏差程度，即：

首先不设置任何参数，均采用默认值，运行GBRT 得到其评价指标为R2=0.993906965，RMSE=0.052775459，仿真结果如图3所示。

图3 GBRT参数优化前预测结果

通过FOA 对GBRT 的关键参数进行寻优后，得到优化后的参数组合为learning rate=0.0125，n estimators=600，max depth=3，min samples split=13，min samples leaf=6，subsample=1，代入GBRT 运行后得到其评价指标为 R2=0.997086676，RMSE=0.036493000，仿真结果如图4。

图4 GBRT参数优化后预测结果

参数优化前后GBRT 的整体预测效果对比见表1。

表1 优化效果对比

由上表可知，GBRT 方法在利用FOA 进行参数优化前，拟合优度已经较高，但均方根误差并不是很理想。对参数进行优化调整后，拟合优度虽然提升较小，只有0.32%的提升，但均方根误差在原基础上降低了30.85%。仿真结果表明，通过FOA 对GBRT 的关键参数进行优化后，可以在保持或略微提升原有拟合优度的基础上，较大幅度地降低预测结果的均方根误差，从而提升该方法的整体预测精度。

3.3 偏移补偿

考虑到实际发酵生产过程中，原材料品质、生物质活性和环境气候等因素的多样性，可能会使预测模型发生状态漂移。为了提高本文方法的适用性能，能应对各种工况的变化。因此本文提出了利用偏移补偿的技术［19］，来改善模型的整体性能。

在FOA-GBRT 方法预测完成后，对输出值进行偏移校正。T 时刻的总偏移量ε(t)，由当前时刻的偏移量ε0(t)和上一时刻的总偏移量ε(t-1)加权求和所得。

其中，y（t-1）和p（t-1）分别代表前一时刻的真实值和预测值。

其中 λ 为权值系数，1 ≤λ ≤0.9，则最终的预测值pfinal为

对3.2 节FOA-GBRT 优化后模型的输出结果进行偏移补偿，其评价指标见表2。

表2 偏移补偿效果对比

4 预测结果与分析

4.1 Pensim仿真平台

本文的实验数据由Pensim 仿真平台产生。此仿真平台是专门为青霉素发酵过程而设计的，该软件内核采用的是基于Bajpai 机理模型改进的Birol模型，在此平台上可以简易实现一系列青霉素发酵过程的仿真，相关研究已经证明了该仿真平台的实用性与有效性。考虑到青霉素发酵过程在不同底物和生物质浓度条件下反应是不同的，因此，本文使用Pensim 仿真平台在三种不同的工况条件下进行了仿真，具体生产条件如表3 所示，在各个工况下（其余条件取仿真平台默认值）仿真产生6 组数据，每组数据仿真时间为400h，采样时间为0.5h。

表3 Pensim仿真初始条件

4.2 各类方法预测结果对比

仿真实验的训练数据集由本文3.1 节所述低、中、高工况中随机选取一共13 组数据作为训练集。从剩下的5 组中随机抽取一组作为测试集。由Pensim 仿真平台产生的每组数据都包含有多达17 种变量，考虑到实际生产过程的测量限制和需求，选取时间、温度、pH、DO 以及二氧化碳浓度共5个变量作为模型的输入，选取产物浓度作为模型的输出。仿真实验环境如下：系统版本：Windows 10家庭中文版（64 位），处理器：Intel（R）Core（TM）i7-5557U CPU @ 3.10GHz，内存：8.00GB，Python 版本：3.6.6。

为了验证本文所提方法的有效性，将建立的FOA-GBRT 软测量方法，与BPNN，SVR，AdaBoost和RF 方法进行比较研究。同时，针对两种不同的方法进行横向比较：Algorithm A，仅使用FOA-GBRT，不进行偏移补偿；Algorithm B，使用FOA-GBRT的基础上，加入偏移补偿。

利用前文整理的训练集与数据集，对上述六种方法进行仿真实验对比，得到其对产物浓度的评价指标如表4所示，预测结果如图5所示，预测误差如图6所示。

表4 六种方法产物浓度预测评价指标对比

由表4 可以知道，本文提供出的Algorithm A 方法相比于其他4 种传统方法，拥有最高的拟合优度和最低的均方根误差，这说明本文提出的FOA-GBRT方法相比于其他的传统方法，确实可以获得更好的预测精度。同时，Algorithm B方法相较于Algorithm A 方法，其预测结果的拟合优度又提升了0.09%，均方根误差又降低了17.22%，这说明本文提出的FOA-GBRT 方法可以与偏移补偿技术很好地结合，通过引入偏移补偿技术，可以有效地提高预测模型的精度，获得更好的预测结果。

由图5 可知，传统的BPNN 方法，SVR 方法和AdaBoost 方法对于青霉素发酵这类非线性的动态系统，难以很好地跟踪青霉素产量的实际趋势，拟合效果较差，预测结果不理想。进一步，采用RF方法，预测结果表明，其预测性能得到了明显的改善，基本上拟合出了青霉素产量的真实值趋势，但其虽然可以捕捉到青霉素发酵过程的非线性动态特性，但是其泛化能力有限，预测结果的误差较大，预测性能不高。最后，基于本文提出的Algorithm A 和Algorithm B方法得到的青霉素产量预测值，相比于前述的四种方法，能够更好地拟合其真实值的变化情况，趋势基本一致，同时模型输出值与真实值之间的误差也更小，保证了较高的预测精度。

图5 六种方法产物浓度预测结果对比

为了进一步地评价各方法的预测性能，本文引入了误差概率分布函数，来直观地展现各方法的预测误差分布情况，如图6 所示是六种不同方法的预测误差概率分布曲线。图中红色曲线是一个标准正态分布。一个性能更优的预测模型，其预测误差概率分布函数曲线应该更接近于一个标准正态分布曲线。由图6 可知，相比于传统的BPNN、SVR、AdaBoost 以及RF 方法，本文提出的FOA-GBRT 方法结合偏移补偿技术所得的青霉素产量预测误差分布最窄，偏差最小，误差绝对值越大，出现频率越低，符合标准正态分布的趋势。这进一步证明了本文所提方法的可靠性。

图6 六种方法产物浓度预测误差概率分布对比

同时，为了充分验证本文所提出的软测量方法是否具有良好且广泛的适用性，进一步将Algorithm B 应用于青霉素发酵过程中底物浓度与生物质浓度的预测，观察其能否保持准确可靠的预测性能。训练数据集、测试数据集以及仿真环境同4.2小节保持不变，模型输入变量仍为时间、温度、pH、DO、二氧化碳浓度，分别修改模型的输出结果为底物浓度和生物质浓度。评价指标见表5，预测结果和误差分布如图7和图8所示。

表5 预测结果评价指标对比

图7 底物浓度预测结果及误差概率分布

图8 生物质浓度预测结果及误差概率分布

由表5 的评价指标可以看出，本文提出的软测量方法对青霉素发酵过程中的底物浓度和生物质浓度的预测结果，都具有高拟合优度和低均方根误差的特点。结合图7 和图8 可知，本文提出的FOA-GBRT 方法能很好地跟踪底物浓度和生物质浓度的真实值变化趋势，同时其预测误差在真实值附近较小区间内波动，偏差小，且没有明显的坏点，概率分布也符合标准正态分布的趋势。以上结果很好地表明，本文提出的软测量方法对底物浓度和生物质浓度也能有很好的预测效果，同时结合本节对产物浓度预测的验证，一定程度上表明了本文提出的软测量方法具有良好的适用性。

考虑到实际工程应用需要，针对Algorithm B，引入相对误差（绝对误差与被测量值真值之比乘以100%所得的数值）来评价其预测值的可信程度，其结果见表6。

表6 预测结果相对误差

将相对误差的最大值定义为预测精度，即本文提出的软测量方法的预测精度为5.42%。

5 结语

针对生物发酵过程中各类难以在线测量的关键变量，本文提出了一种采用果蝇算法优化梯度提升回归树的发酵过程软测量方法。该方法避开了生物发酵过程复杂的机理建模过程，能很好地适应其非线性、大时滞、时变性以及多工况的特点，有很好的预测效果。同时，本文还利用果蝇算法对梯度提升回归树的关键参数进行寻优，并结合偏移补偿技术对输出值进行校正，进一步地缩小了预测结果的偏差，改善了方法的整体性能。

基于Pensim 仿真平台所得的青霉素发酵过程数据进行仿真验证，结果表明本文提出的FOA-GBRT软测量方法，能够准确地预测青霉素发酵过程的产量、底物和生物质浓度，预测精度在5.42%以内，明显优于BPNN、SVR 和RF 等方法，为发酵过程的关键变量预测提供了可能。