基于可控參數(shù)的前饋神經(jīng)網(wǎng)絡(luò)出水總氮預(yù)測模型研究

趙子豪，王子昊，袁家洛，馬駿，何哲靈，徐一蘭，沈曉佳，朱亮,b,c,*

a Institution of Environment Pollution Control and Treatment, Department of Environmental Engineering, Zhejiang University, Hangzhou 310058, China

b Zhejiang Province Key Laboratory for Water Pollution Control and Environmental Safety, Hangzhou 310058, China

c Zhejiang Provincial Engineering Laboratory of Water Pollution Control, Hangzhou 310058, China

d Haining Water Investment Group Co., Ltd, Haining 314400, China

e Haining Capital Water Co., Ltd, Haining 31440, China

1. 引言

近年來我國工業(yè)化城鎮(zhèn)化進程不斷加快，污水排放量攀升至每年669.6億t以上[1]。為削弱污水排放所帶來的環(huán)境不利影響，我國陸續(xù)采取新建或改擴建污水處理廠、修訂污水排放標準等一系列措施。截至2018年年底，我國共有5370座污水處理廠處于正常運行狀態(tài)，日處理能力達2.01×108m3，總用電量達1.972×1010kW·h-1[2]。據(jù)調(diào)查，目前我國多數(shù)城鎮(zhèn)污水處理廠存在出水水質(zhì)難以穩(wěn)定達標、運行藥耗能耗偏高、自動化控制程度較低等問題，亟需科技支撐予以解決。

序批式活性污泥工藝（SBR）因其工藝簡便、操作方式靈活、抗負荷性能優(yōu)異等特點，已成為我國應(yīng)用較為廣泛的污水處理工藝之一[3]。然而，SBR工藝對系統(tǒng)操控自動化水平要求較高，極易產(chǎn)生出水水質(zhì)指標（尤其是總氮）波動大、運行能耗過高等問題[4]。而在SBR工藝中，脫氮是十分復(fù)雜的，包括好氧條件下的硝化反應(yīng)及缺氧條件下的反硝化反應(yīng)[5]。已有研究表明，溶解氧（DO）控制是實現(xiàn)污水處理系統(tǒng)總氮（TN）有效去除的重要因素，充足的DO可保證有機物降解及硝化反應(yīng)順利進行，但過量的DO則會導(dǎo)致污水處理能耗過高、反硝化效率降低、污泥性能惡化[6]。同時，根據(jù)現(xiàn)有工程運行情況來看，污水處理系統(tǒng)DO難以實現(xiàn)精準控制，主要與DO監(jiān)測設(shè)備準確度低、響應(yīng)慢，傳統(tǒng)擬合模型和控制理論相對落后有關(guān)[7,8]。在大多數(shù)污水處理廠，操作人員往往根據(jù)經(jīng)驗調(diào)整表面氣速、缺氧段時長等控制參數(shù)，從而實現(xiàn)體系DO濃度的相應(yīng)調(diào)控，這樣的操作方式也比DO精準控制更為簡便。

鑒于污水生物處理系統(tǒng)參數(shù)眾多、過程復(fù)雜，傳統(tǒng)的生化模型難以擬合其高維的非線性數(shù)據(jù)，進而難以應(yīng)用到實際污水生物處理工藝中[9,10]。而人工神經(jīng)網(wǎng)絡(luò)（ANN）作為一種自學(xué)習(xí)算法模型，理論上能夠擬合任意非線性函數(shù)，因而其可作為預(yù)測污水生物處理這種復(fù)雜非線性系統(tǒng)的有效工具[11-13]。從神經(jīng)元拓撲結(jié)構(gòu)來看，ANN可分為反饋神經(jīng)網(wǎng)絡(luò)和前饋神經(jīng)網(wǎng)絡(luò)（FFNN），其中，F(xiàn)FNN模型理論上可以任意精度逼近任意連續(xù)函數(shù)，具有較強的分類和模式識別能力[14-19]。為進一步提高FFNN的預(yù)測效率與能力，研究者陸續(xù)提出了Levenberg-Marquardt (L-M)、貝葉斯正則化（BR）、量化共軛梯度（SCG）、動量和Nesterov加速梯度等多種優(yōu)化算法，顯著提高FFNN模型的預(yù)測能力和效率[20-22]。

本文以序批式活性污泥工藝為研究平臺，構(gòu)建了一種基于可控參數(shù)的前饋神經(jīng)網(wǎng)絡(luò)出水總氮預(yù)測模型。與已有預(yù)測模型相比，本模型具備以下兩個特點：①采用可控參數(shù)（表面氣速與缺氧段時長）代替溶解氧作為模型主要輸入?yún)?shù)，明顯提高模型可用可控性；②采用算法優(yōu)化的FFNN構(gòu)建模型，顯著提高模型預(yù)測精準度。本文的研究目的主要在于評估優(yōu)化后FFNN模型的SBR工藝出水TN預(yù)測準確性，明確合理地運行參數(shù)控制策略以實現(xiàn)污水處理廠污染物高效去除與系統(tǒng)節(jié)能降耗。

2. 材料與方法

2.1. SBR工藝長期仿真實驗

本研究設(shè)置兩組平行運行的SBR反應(yīng)器（R1、R2），開展持續(xù)兩個月的長期仿真實驗。反應(yīng)器活性污泥取自浙江省某城鎮(zhèn)污水處理廠，進水模擬該污水處理廠的實際進水。反應(yīng)器結(jié)構(gòu)及模擬廢水組成如附錄A中的圖S1、表S1所示。

參考該污水處理廠SBR工藝的實際運行模式，每組反應(yīng)器的運行周期設(shè)為4 h，包括進水（5 min）、缺氧與曝氣期（210 min）、沉降期（5 min）、出水（5 min）、閑置（15 min），體積交換比設(shè)為50%。配置進水時，隨機控制主要水質(zhì)指標濃度為預(yù)設(shè)值的75%～125%，模擬實際進水水質(zhì)波動情況。分別設(shè)表面氣速與缺氧段時長作為控制變量，其設(shè)計值如表1所示。每個周期分別對進出水進行取樣，采用標準方法分析進出水的總氮（TN）、氨氮（NH4+-N）、化學(xué)需氧量（COD）和總磷（TP）濃度[23]。經(jīng)兩個月的仿真運行后，從控制變量的16種組合中獲取124組數(shù)據(jù)供后續(xù)建模使用。

表1 SBR工藝控制變量的設(shè)計值

為擬合更復(fù)雜的運行狀況，本研究引入缺氧段后進水占比并設(shè)置范圍更廣的控制變量后開展擴展實驗，共收集91組數(shù)據(jù)，其設(shè)計值如表2所示。同時為仿真實際運行中的極端情況，本研究分別設(shè)置表面氣速為3.6 cm·s-1、4.8 cm·s-1，缺氧段時長為0、150 min，并從其組合中收集11組數(shù)據(jù)。

表2 擴展實驗中引入缺氧段后控制變量和進水比的設(shè)計值

2.2. FFNN建模

基礎(chǔ)FFNN模型及其優(yōu)化算法均在Matlab R2016a中搭建并運行。建模前，將長期仿真得到的數(shù)據(jù)集歸一化至0.001～0.999，以消除數(shù)據(jù)間不同量綱的影響。FFNN模型在構(gòu)建過程中根據(jù)模型預(yù)測值與實際值的誤差不斷優(yōu)化權(quán)值網(wǎng)絡(luò)[式（1）] [24]，其網(wǎng)絡(luò)結(jié)構(gòu)包括輸入層、隱藏層及輸出層三部分，其中進水水質(zhì)（COD、TN、NH4+-N、TP）與運行參數(shù)（表面氣速、缺氧段時長）6個變量作為模型的輸入變量，出水TN則作為模型的輸出變量，建模示意圖如圖1所示。由于實驗數(shù)據(jù)維度較低，設(shè)置FFNN模型隱藏層數(shù)為1層以縮短運算時長，提高效率，同時防止過擬合現(xiàn)象發(fā)生。而模型隱藏層最佳結(jié)點數(shù)則參考經(jīng)驗方程進行計算[式（2）]。

式中，a為輸出變量；W為權(quán)值矩陣；d為輸入變量；b為偏置。

式中，h為隱藏層結(jié)點數(shù)；i為輸入層結(jié)點數(shù)；o為輸出層結(jié)點數(shù)；c為常數(shù)，一般選在1～10之間。

2.3. FFNN模型優(yōu)化

基礎(chǔ)FFNN模型參照傳統(tǒng)的反向傳播（BP）神經(jīng)網(wǎng)絡(luò)進行構(gòu)建。在訓(xùn)練階段，BP神經(jīng)網(wǎng)絡(luò)根據(jù)預(yù)測值與實際值誤差的反向傳播不斷調(diào)整其權(quán)重網(wǎng)絡(luò)，使誤差最小化。梯度下降法則是模型調(diào)整權(quán)重的最常用算法，即沿梯度下降的方向調(diào)整權(quán)重并最小化誤差[25]。但在實際訓(xùn)練過程中，梯度下降法在迭代求解時容易最小化誤差至局部極小值而不是全局最小值，因而降低了模型學(xué)習(xí)效率與預(yù)測準確度[26]。為提高模型的學(xué)習(xí)效率與預(yù)測準確度，本研究采用三種優(yōu)化算法（L-M、BR以及SCG）對FFNN模型進行優(yōu)化，對比三種優(yōu)化算法得到最適合SBR工藝的優(yōu)化FFNN模型。同時，為進一步提高模型的泛化能力，采用更復(fù)雜的數(shù)據(jù)集（即包括擴展實驗數(shù)據(jù)集）訓(xùn)練模型，缺氧段后進水比將代替進水TP作為模型的輸入?yún)?shù)。

2.3.1. L-M算法

L-M算法結(jié)合高斯-牛頓算法以及梯度下降法的優(yōu)點，能有效避免局部最小值并提高全局最小值的收斂速度，其算法如下。

若將BP神經(jīng)網(wǎng)絡(luò)中各層間的權(quán)值用向量W表示，其誤差平方和（E）為

式中，n為樣本編號；tnj為樣本n在輸出層j結(jié)點的期望輸出；Onj為該結(jié)點實際輸出；ε是以εn為元素的向量。

在式（4）中，k代表迭代次數(shù)（即權(quán)重更新的次數(shù)）。在計算Wk+1的過程中，如果移動量Wk+1-Wk很小，則可將ε展開成一階泰勒級數(shù)：

式中，Z是ε的雅可比矩陣，Z的元素為

由此，誤差函數(shù)可改為

要使誤差函數(shù)最小，可對Wk+1進行求導(dǎo)，得到如下高斯-牛頓迭代公式：

式中，T代表轉(zhuǎn)置。

圖1. 出水總氮預(yù)測FFNN模型的網(wǎng)絡(luò)結(jié)構(gòu)。

為克服高斯-牛頓法中經(jīng)常出現(xiàn)的雅可比矩陣奇異現(xiàn)象，把誤差函數(shù)改寫為

式中，λ為阻尼系數(shù)。

此時對E求導(dǎo)，便可得到基于高斯-牛頓法的L-M迭代公式：

式中，I為單位矩陣，為迭代變量。在迭代過程中搜索方向與訓(xùn)練步長受λ變化影響。在計算初始階段λ可取較大值，這時ZTZ與λI相比是可以忽略的，于是上式可寫為

式中，g為梯度。若令λ趨于零，則該表達式可變?yōu)楦咚?牛頓迭代公式。

2.3.2. BR算法

BR算法指貝葉斯方法正則化神經(jīng)網(wǎng)絡(luò)。所謂正則化是指訓(xùn)練階段為神經(jīng)網(wǎng)絡(luò)添加懲罰項來限制網(wǎng)絡(luò)的復(fù)雜度。經(jīng)正則化后，神經(jīng)網(wǎng)絡(luò)可有效避免過擬合現(xiàn)象的發(fā)生，提高模型的泛化能力。通常來講，神經(jīng)網(wǎng)絡(luò)的性能函數(shù)可表示為

經(jīng)過添加懲罰項EW后，其性能函數(shù)轉(zhuǎn)變?yōu)?/p>

α和β的相對大小決定了懲罰項所占比例。如果α＜＜β，則近似于無正則化條件，側(cè)重于最小化訓(xùn)練誤差，但可能過擬合。若α＞＞β，則側(cè)重于限制網(wǎng)絡(luò)權(quán)值規(guī)模，但可能降低模型預(yù)測性能。因此，如何確定α和β值是非常重要的。在貝葉斯分析的框架下，MK推導(dǎo)出[27]：

式中，γ=N-2αtr(H)-1表示有效權(quán)值數(shù)，N為樣本總數(shù)；H是F的海森矩陣：

但海森矩陣計算量很大，F(xiàn)oresee和Hagan用高斯-牛頓法近似計算海森矩陣，大大降低了計算量[28]：

式中，J是訓(xùn)練誤差的雅可比矩陣。

2.3.3. SCG算法

SCG算法是一種標準BP神經(jīng)網(wǎng)絡(luò)的改進算法。在神經(jīng)網(wǎng)絡(luò)訓(xùn)練中，傳統(tǒng)的梯度下降法只能沿負梯度的方向進行收斂，并不一定會正確收斂到全局極小點。共軛梯度法則可將負梯度方向與上一次搜索方向結(jié)合起來，計算出新的搜索方向。但該算法每次迭代時都需要重新確定搜索方向，進而導(dǎo)致其計算量巨大。Moller [29]提出的SCG算法將模型置信區(qū)間法與共軛梯度法結(jié)合起來，成功解決了上述問題。SCG算法的權(quán)值調(diào)整方法轉(zhuǎn)變?yōu)?/p>

式中，wk為Wk中一點；pk為第k次迭代時的搜索方向；θk為第k次迭代時的搜索步長。

式中，gk為第k次迭代時的梯度；Hk為第k次迭代時的海森矩陣。

令sk=Hkpk,δk=pkTsk,uk= －gkTpk，則有θk=uk /δk。同時為了保證δk＞ 0，令δk=pkTsk+λk|pk|2。

最后可得到步長θk：

2.4. 模型性能評價指標

本研究采用相關(guān)系數(shù)（R）[式（21）]評估各個模型的預(yù)測能力，采用均方誤差（MSE）[式（22）]評估訓(xùn)練過程中各模型的擬合能力。R與MSE均在Matlab R2016a中計算完成。

式中，Cov (TNP,TNA)代表出水TN預(yù)測值與實際值的協(xié)方差；σp、σA分別為出水TN預(yù)測值與實際值的標準差；m為數(shù)據(jù)集樣本量；Pi為第i個樣本的出水TN預(yù)測值；Ai為第i個樣本的出水TN實際值。

3. 結(jié)果與討論

3.1. 采用運行參數(shù)作為主要輸入變量構(gòu)建FFNN模型的可行性

建模時首先確定神經(jīng)網(wǎng)絡(luò)隱藏層結(jié)點數(shù)。在其他參數(shù)不變的情況下，根據(jù)式（2）計算可得，隱藏層結(jié)點數(shù)q應(yīng)取在3～12之間。訓(xùn)練完成后，不同結(jié)點數(shù)模型的均方誤差如圖2所示。對比發(fā)現(xiàn)，神經(jīng)網(wǎng)絡(luò)隱藏層結(jié)點數(shù)為8時，MSE最小，因而確定基礎(chǔ)FFNN模型的網(wǎng)絡(luò)結(jié)構(gòu)為6-8-1。

為訓(xùn)練和評估基礎(chǔ)FFNN模型，本研究將SBR工藝長期仿真實驗中所收集的數(shù)據(jù)隨機分為訓(xùn)練集（104組）和測試集（20組）。圖3為模型對整體數(shù)據(jù)集的擬合程度，可見其左側(cè)（訓(xùn)練集）的預(yù)測值與實際值較為接近（R= 0.91973），表明采用可控參數(shù)（如表面氣速、缺氧段時長）代替DO作為輸入構(gòu)建模型是可行的。但右側(cè)測試集中預(yù)測值與實際值相差甚遠（R= 0.5057），說明該模型仍存在過擬合或陷于局部極小值等降低模型預(yù)測性能的問題。針對這些問題，本研究采用L-M、BR和SCG三種算法對神經(jīng)網(wǎng)絡(luò)模型進行優(yōu)化。同時采用更復(fù)雜的數(shù)據(jù)集對模型進行訓(xùn)練，并針對不同的算法建立更合理的訓(xùn)練模式，找出最適合模型的算法。

圖2. 具有不同隱藏層結(jié)點數(shù)的FFNN模型的均方誤差。

3.2. 采用L-M算法優(yōu)化FFNN模型

將數(shù)據(jù)集隨機分為訓(xùn)練集（80%）和驗證集（20%）。通過嘗試法優(yōu)選后，設(shè)置隱藏層結(jié)點數(shù)為30，最小失敗次數(shù)為100，隱藏層與輸出層間傳遞函數(shù)為tansig。訓(xùn)練113步后，因驗證集失敗次數(shù)超過設(shè)定值，停止訓(xùn)練，共用時1.0 s，其擬合結(jié)果如圖4所示。由圖可知，該模型的MSE為0.00200，在訓(xùn)練集、驗證集與全集中的R值分別為：0.93603、0.92316、0.93392，表明該模型對出水TN的預(yù)測性能良好。L-M算法優(yōu)化后的模型比之前的基礎(chǔ)FFNN模型具有更高的R值和更快的收斂速度，不僅有效避免了局部最小值，而且提高了收斂到全局最優(yōu)的速度。

圖3. FFNN模型預(yù)測值與實際值的比較。

3.3. 采用BR算法優(yōu)化FFNN模型

BR算法不需要驗證集，因此將數(shù)據(jù)隨機分為訓(xùn)練集（80%）與測試集（20%）。模型調(diào)試過程發(fā)現(xiàn)，當隱藏層結(jié)點數(shù)小于20 時，預(yù)測值與實際測量值之間的相關(guān)系數(shù)低于 0.6，表明模型擬合程度較差；隨著隱藏層結(jié)點數(shù)逐漸增加至30，R值逐漸升高，結(jié)點數(shù)超過30后R值升高不明顯，但訓(xùn)練時間大幅增長，故選取隱藏層結(jié)點數(shù)為30。選取Tansig作為傳遞函數(shù)，設(shè)置訓(xùn)練步數(shù)上限為1000后，該模型訓(xùn)練用時19 s。如圖5所示，訓(xùn)練集的R值為0.90232，而測試集的R值為0.83000，說明該模型的預(yù)測性能優(yōu)于基礎(chǔ)模型，但弱于L-M算法優(yōu)化的模型。BR算法通過限制網(wǎng)絡(luò)的復(fù)雜性提高了模型性能，使得R值比基礎(chǔ)模型更高。盡管如此，正規(guī)化可能會使得模型丟失一些輸入輸出之間的關(guān)系，特別是在模型結(jié)構(gòu)復(fù)雜性中等的情況下。

3.4. 采用SCG算法優(yōu)化FFNN模型

與L-M算法相似，數(shù)據(jù)集被劃分為訓(xùn)練集（80%）和驗證集（20%）。隱層結(jié)點數(shù)選為30，以保證與三種算法模型結(jié)構(gòu)一致。選擇Tansig作為隱藏層與輸出層之間的傳遞函數(shù)，最小失敗次數(shù)設(shè)為100。因失敗次數(shù)超過預(yù)設(shè)值，該模型在訓(xùn)練156步后停止，其訓(xùn)練時間接近于零。如圖6所示，其訓(xùn)練集的R值為0.916，驗證集的R值為0.962。驗證集R值高于其訓(xùn)練集R值表明該模型不存在過擬合現(xiàn)象且具有良好的預(yù)測性能。

圖4. 基于L-M算法的FFNN模型仿真結(jié)果。（a）訓(xùn)練集（輸出≈0.88×目標+0.075）；（b）驗證集（輸出≈0.88×目標+0.074）；（c）全集（輸出≈0.88×目標+0.077）?！癋it”一行表示目標和輸出之間的關(guān)系；“O=T”一行表示目標等于輸出。

圖5. 基于BR算法的FFNN模型仿真結(jié)果。（a）訓(xùn)練集（輸出≈0.76×目標+0.15）；（b）測試集（輸出≈0.75×目標+0.15）；（c）全集（輸出≈0.76×目標+0.15）。

3.5. 三種優(yōu)化算法的比較

基于不同優(yōu)化算法的FFNN模型訓(xùn)練參數(shù)如表3所示。由表可知，訓(xùn)練期間三種算法優(yōu)化后的FFNN模型的MSE均較低（＜ 0.001），表明其擬合能力均優(yōu)于基礎(chǔ)模型（MSE = 0.0037）。其中，SCG算法的訓(xùn)練時間最短。

表3 基于不同優(yōu)化算法的FFNN模型訓(xùn)練參數(shù)

基于各優(yōu)化算法所構(gòu)建的FFNN模型R值如表4所示。訓(xùn)練期間三種模型的R值分別為0.936、0.902和0.917，說明模型訓(xùn)練集的輸入（進水水質(zhì)和運行參數(shù)）和輸出（出水TN）相關(guān)性較強。在預(yù)測階段，所有FFNN模型的R值均高于基礎(chǔ)模型，表明優(yōu)化算法大大提高了FFNN神經(jīng)網(wǎng)絡(luò)模型對出水總氮的預(yù)測性能。其中，基于SCG算法的模型R值最高并高于其訓(xùn)練集的R值且訓(xùn)練時間最短，表明SCG算法能夠使模型獲得最快的計算速度和出色的預(yù)測能力，是最佳的優(yōu)化算法。

表4 基于各優(yōu)化算法所構(gòu)建的FFNN模型R值

圖6. 基于SCG算法的FFNN模型訓(xùn)練結(jié)果。（a）訓(xùn)練集（輸出≈0.84×目標+0.1）；（b）驗證集（輸出≈0.88×目標+0.079）；（c）全集（輸出≈0.85×目標+0.092）。

3.6. 與其他模型的比較

針對厭氧-缺氧-好氧工藝（A2/O）、SBR和氧化溝等不同的污水處理過程，研究者陸續(xù)提出了各種適合于特定環(huán)境仿真模型用來預(yù)測出水水質(zhì)。如為預(yù)測屠宰場廢水中NH4+-N的去除率，Kundu等[30]構(gòu)建了基于進水水質(zhì)和DO的BP神經(jīng)網(wǎng)絡(luò)模型?；谶M水參數(shù)構(gòu)建的ANN模型可預(yù)測好氧顆粒污泥的TN去除效率[31]。Ebrahimi等[32]提出了多變量回歸模型預(yù)測氧化溝的出水生化需氧量（BOD）、TP濃度。然而，上述研究選擇的輸入?yún)?shù)主要為DO與進水水質(zhì)，具有不可控性。與以往研究相比，本研究從仿真模型的可控性與預(yù)測精度兩個方面出發(fā)，首次證實FFNN模型可采用實際可控參數(shù)（如表面氣速和缺氧時間）代替不可控的DO實現(xiàn)出水總氮預(yù)測。為提高預(yù)測精度，采用多種算法對FFNN模型進行優(yōu)化，SCG優(yōu)化后的模型預(yù)測精度最高。在此基礎(chǔ)上，將構(gòu)建的模型與多種機器學(xué)習(xí)模型對比（表 5），結(jié)果表明經(jīng)SCG算法優(yōu)化的FFNN模型預(yù)測性能最佳，具有最小的RMSE和最高的R值?？偟膩砜?，研究構(gòu)建的模型可根據(jù)進水水質(zhì)變化調(diào)節(jié)可控運行參數(shù)，以保證出水TN處于穩(wěn)定達標狀態(tài)，同時有效避免過度曝氣 [33]，降低系統(tǒng)運行能耗。

表5 不同預(yù)測方法的性能比較（R值四舍五入為小數(shù)點后兩位數(shù)）

4. 結(jié)論

本研究基于主流SBR工藝構(gòu)建了一種基于可控參數(shù)的前饋神經(jīng)網(wǎng)絡(luò)出水總氮預(yù)測模型，選擇表面氣速與缺氧時長代替DO作為FFNN模型的主要輸入?yún)?shù)，運用長期仿真數(shù)據(jù)提高了模型可用可控性；選擇SCG作為FFNN模型優(yōu)化算法，實現(xiàn)出水TN及最佳運行參數(shù)的精準預(yù)測，可為污水處理系統(tǒng)的穩(wěn)定運行及節(jié)能降耗提供有價值的解決方案。

致謝

本文得到了國家水體污染控制與治理科技重大專項（2017ZX07201003）、國家自然科學(xué)基金（51961125101）、浙江省重點研發(fā)計劃（2018C03003）的大力支持。

Compliance with ethics guidelines

Zihao Zhao, Zihao Wang, Jialuo Yuan, Jun Ma, Zheling He, Yilan Xu, Xiaojia Shen, and Liang Zhu declare that they have no conflict of interest or financial conflicts to disclose.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2020.07.027.