基于梯度提升回歸樹的發酵過程軟測量研究?

黃繼煒 潘 豐

（江南大學輕工過程先進控制教育部重點實驗室 無錫 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 等方法，為發酵過程的關鍵變量預測提供了可能。