王 立, 鄭 浪, 金學(xué)波, 王小藝, 于家斌, 白玉廷

(北京工商大學(xué)人工智能學(xué)院1,北京 100048) (北京服裝學(xué)院2,北京 100105)

我國是糧食生產(chǎn)大國,然而在我國約有70%的糧食在儲藏過程中由于不當(dāng)?shù)膬Σ貤l件導(dǎo)致真菌大量繁殖而發(fā)生霉變現(xiàn)象[1],并產(chǎn)生大量真菌毒素[2],使得糧食被真菌污染導(dǎo)致食用品質(zhì)下降。因此,確保糧食儲藏安全至關(guān)重要。為了保證糧食儲藏安全,學(xué)者針對真菌污染問題做了大量的相關(guān)研究,現(xiàn)有的真菌污染治理方法大多都是對已經(jīng)產(chǎn)生了的真菌進行消除,例如物理方法[3-5]、化學(xué)方法[6]和生物方法[7]。

如果通過預(yù)測糧食儲藏過程中真菌的變化趨勢,在真菌大量生長繁殖之前就對糧食儲藏環(huán)境進行控制,更有利于降低糧食霉變和真菌毒素污染的風(fēng)險。因此對糧食儲藏過程中真菌生長趨勢進行分析和預(yù)測,并對溫度和水分等真菌生長的主要影響因素進行控制[8],對于降低糧食霉變和真菌毒素污染的風(fēng)險具有重要意義。

目前對糧食儲藏過程中真菌生長的預(yù)測方法主要分為基于機理驅(qū)動的建模方法和基于數(shù)據(jù)驅(qū)動的建模方法。其中基于機理驅(qū)動的模型是結(jié)合微生物學(xué)、數(shù)學(xué)等學(xué)科,采用數(shù)學(xué)公式描述真菌的生長與影響因素之間的關(guān)系,例如Logistic方程[9],Rosso cardinal模型[10],Ratkowsky平方根模型[11]等?；跀?shù)據(jù)驅(qū)動的模型是通過機器學(xué)習(xí)算法,使用大量的數(shù)據(jù)對真菌生長過程的特征進行學(xué)習(xí),達到預(yù)測真菌生長趨勢的目的,例如線性回歸模型[12],BP神經(jīng)網(wǎng)絡(luò)[13],支持向量機[14]等。然而,在目前針對真菌生長建模的研究中,同時考慮多種環(huán)境因素對真菌生長影響的研究較少,導(dǎo)致現(xiàn)有真菌生長模型難以精確描述和預(yù)測真菌生長趨勢。

降低糧食儲藏過程真菌污染的風(fēng)險不僅需要準(zhǔn)確地預(yù)測糧食儲藏期間真菌的變化趨勢,還需要能夠?qū)崟r控制糧食儲藏環(huán)境。影響真菌生長的主要環(huán)境因素包括水分、溫度、光照、pH等[4,15-18],但實際糧食儲藏大多采用干燥的倉儲形式,所以光照、pH 2種環(huán)境因素暫不涉及,而水分是決定真菌生長的必要條件,溫度影響真菌生長的速度,因此通過對這兩種影響因素監(jiān)測,并采用制冷[19]、通風(fēng)[20]等方法進行調(diào)節(jié),可以實現(xiàn)對糧食儲藏環(huán)境的控制。目前對于糧食儲藏環(huán)境控制的研究比較粗略,例如對糧堆溫度的預(yù)測值與預(yù)警值進行對比,僅對危險的糧堆給出處置建議[21],且無法對糧食儲藏的多種環(huán)境因素進行優(yōu)化控制以減少真菌的產(chǎn)生。

研究同時考慮溫度和水分2種環(huán)境因素對糧食儲藏過程中真菌生長的影響,探索糧食儲藏過程中真菌生長規(guī)律,預(yù)測真菌生長趨勢,對影響真菌生長的溫度和水分因素進行優(yōu)化控制,以減少糧食儲藏期間真菌的產(chǎn)生。在此基礎(chǔ)上,采用稻谷儲藏過程的仿真實驗[16]進行方法驗證。

1 溫度和水分影響下的真菌生長建模預(yù)測

將溫度、水分作為真菌生長的主要影響因素,對Rosso cardinal模型進行改進,首先建立單個菌落的直徑生長模型,然后基于改進的二維正態(tài)分布模型,將模型轉(zhuǎn)換為真菌孢子數(shù)增長模型。

1.1 基于Rosso cardinal模型改進的菌落直徑生長模型

Rosso cardinal模型[10]是一個用溫度來描述菌落直徑生長速率的模型,相對于Logistic方程,Ratkowsky平方根模型等模型,Rosso cardinal模型所有參數(shù)都具有生理學(xué)意義,但是該模型未考慮其他因素對菌落直徑生長的影響,因此將其擴展為溫度、水分2個影響因素,建立單個菌落的直徑生長速率模型,以提高菌落生長的預(yù)測精度,如式(1)所示:

μ(T,Me)=μopt·τ(T)·ρ(Me)

(1)

式中:μ(T,Me)為溫度為T、水分為Me時的單個菌落直徑生長速率;T為儲藏環(huán)境的溫度;Me為儲藏環(huán)境的水分;μopt為菌落直徑生長最快時候的徑向生長速率;τ(T)為儲藏環(huán)境溫度為T時對生長速率的影響;ρ(Me)為儲藏環(huán)境水分為Me時對生長速率的影響。

τ(T)和ρ(Me)用式(2)和式(3)計算。

(2)

(3)

式中:Tmin、Tmax和Topt分別為能夠使菌落生長的最低溫度、最高溫度和生長最快時的溫度;Me min、Me max和Me opt分別為能夠使菌落生長的最低水分、最高水分和生長最快時的水分,6個參數(shù)均為固定參數(shù)。

建立單個菌落的直徑生長速率模型后,對其進行時間的積分,記t時刻菌落直徑為dt,以Tt、Met表示變化的溫度、水分,就可以得到單個菌落的直徑生長模型,如式(4)所示。

(4)

1.2 基于二維正態(tài)分布模型改進的真菌孢子數(shù)增長模型

相比菌落直徑,真菌孢子數(shù)更能直觀反映真菌數(shù)量,且實際實驗中,真菌孢子數(shù)更易監(jiān)測,因此將菌落直徑生長模型轉(zhuǎn)換為真菌孢子數(shù)增長模型。

根據(jù)實際單個菌落分布的特點,通過改進二維正態(tài)分布模型建立單個菌落直徑與真菌孢子數(shù)之間的關(guān)系,得到真菌孢子數(shù)增長模型。假設(shè)菌落的數(shù)量分布服從二維正態(tài)分布,二維正態(tài)分布的概率密度函數(shù)與均值μ和形狀參數(shù)σ有關(guān),令菌落中心為坐標(biāo)原點,即μ=0,那么二維正態(tài)分布的概率密度函數(shù)為式(5):

(5)

圖1 菌落的形狀

由圖1可知,菌落直徑越大,菌落中心的真菌孢子數(shù)越多,在式(5)中,系數(shù)C表示(x,y)為0時菌落中心的真菌孢子數(shù),所以系數(shù)C的大小與菌落直徑有關(guān)。設(shè)二維正態(tài)分布的形狀參數(shù)σ和菌落直徑成正比,則將系數(shù)C表達式重新設(shè)定為aσ,對系數(shù)C的驗證由3.2節(jié)給出,真菌孢子數(shù)增長模型如式(5)所示。

(6)

式中:f(x,y)為單個菌落在坐標(biāo)(x,y)處對應(yīng)的真菌孢子數(shù);a為真菌孢子數(shù)增長模型的待定參數(shù)。

記t時刻的σ值為σt,t時刻單個菌落的直徑為dt,σt與dt的轉(zhuǎn)換關(guān)系如式(7)所示。

σt=kdt

(7)

式中:k是待定的比例系數(shù)。對f(x,y)在菌落范圍x2+y2≤dt2內(nèi)積分即可得到單個菌落范圍的全部真菌孢子數(shù),計算如式(8)所示。

(8)

對于分布函數(shù)而言,(x,y)的范圍可以取無窮遠(yuǎn),但是對于單個菌落而言,菌落直徑有限,存在菌落邊緣,因此定義fmin為菌落邊緣處的真菌孢子數(shù),即在菌落邊緣處滿足f(x,y)=fmin,當(dāng)f(x,y)≥fmin時視為(x,y)在菌落的直徑范圍內(nèi);當(dāng)f(x,y)≤fmin時,即真菌孢子數(shù)f(x,y)太小,視為不在菌落的直徑范圍內(nèi)。式(9)所示關(guān)系:

(9)

化簡得式(10)所示關(guān)系:

(10)

將式(10)帶入式(8)得到式(11),即得到單個菌落范圍的全部真菌孢子數(shù)。

(11)

再將式(7)帶入式(11),得到式(12):

(12)

假設(shè)在每克糧食中平均有m個菌落,記t時刻每克糧食m個菌落中的真菌孢子數(shù)為St,則St的計算為式(13):

St=m·[2πak3dt3+(-2πfmink2dt2)]

(13)

1.3 菌落直徑和真菌孢子數(shù)的預(yù)測模型

預(yù)測模型可通過對菌落直徑生長模型和真菌孢子數(shù)增長模型進行離散化和增量化得到。采樣間隔取Δt=1(d),假設(shè)2次采樣間隔之間儲藏環(huán)境的溫度、水分保持不變,記dl為離散化后的第l次采樣的菌落直徑,則式(4)可以寫為式(14):

(14)

式中:t=0,1,2,…,l-1。進行增量化得式(15)。

(15)

式中:dl+1為l+1次采樣的菌落直徑;Tl和Mel為第l次采樣時的溫度和水分,取Δdl+1=dl+1-dl,即從l到l+1次采樣間隔內(nèi)的菌落直徑的增量。將式(1)所示的菌落直徑生長速率模型寫成溫度、水分隨時間變化的時候,即溫度、水分寫為Tt、Met,帶入式(15)得到式(16)。

dl+1=dl+μoptτ(Tl)ρ(Mel)Δt

(16)

式中:τ(Tl)為儲藏環(huán)境溫度為Tl時對生長速率的影響;ρ(Mel)為儲藏環(huán)境水分為Mel時對生長速率的影響。式(14)的菌落直徑帶入式(13),得到式(17)。

(17)

以Sl表示第l次采樣的真菌孢子數(shù),分別記A=2πamk3μopt3,B=-2πmfmink2μopt2,那么得到式(18)。

(18)

考慮溫度、水分隨時間變化的情況,將式(1)的溫度、水分表達式中的τ(T)ρ(Me)引入時間變量t,帶入式(18),得到式(19)。

(19)

由于μopt與時間無關(guān),對式(19)化簡,得到溫度、水分隨時間改變的Sl+1如式(20)所示。

(20)

當(dāng)溫度水分不隨時間變化時,式(20)簡寫為式(21)。

Sl+1=A[(l+1)τ(T)ρ(Me)]3+

B[(l+1)τ(T)ρ(Me)]2

(21)

式(21)即為溫度、水分不隨時間變化時候每克糧食的真菌孢子數(shù)預(yù)測模型,由于a、m、k、fmin都是待定參數(shù),所以A、B是待定參數(shù),需要通過擬合得到。即式(15)作為菌落直徑預(yù)測模型,式(20)作為真菌孢子數(shù)預(yù)測模型。

2 溫度和水分影響下的真菌生長優(yōu)化控制

對糧食儲藏過程中真菌生長進行控制,一方面,需要預(yù)測真菌生長的趨勢,通過調(diào)節(jié)當(dāng)前時刻環(huán)境影響因素,以降低未來時刻真菌生長速率;另一方面,環(huán)境影響因素的調(diào)節(jié)量越大可使真菌生長速率越低,但考慮到控制成本的經(jīng)濟性,應(yīng)使影響因素的調(diào)節(jié)量盡可能小,以使控制成本最小化。因此可采用預(yù)測控制理論對溫度和水分影響因素進行優(yōu)化控制,達到全局最優(yōu)化的目的。

2.1 真菌生長的預(yù)測控制模型的改進

傳統(tǒng)預(yù)測控制系統(tǒng)設(shè)計的優(yōu)化目標(biāo)是被控對象的狀態(tài)變量跟蹤其參考軌跡,當(dāng)預(yù)測步長為1時,使?fàn)顟B(tài)變量的跟蹤偏差(yr(l+1)-y(l+1))最小化;同時,考慮控制成本約束,使控制變量的增量(Δu(l)=u(l)-u(l-1))最小化。目標(biāo)函數(shù)一般如式(22)所示。

Jl=m1(yr(l+1)-y(l+1))+

(22)

式中:Jl為第l次優(yōu)化的目標(biāo)函數(shù);m1、m2為權(quán)重參數(shù),傳統(tǒng)預(yù)測控制模型框圖如圖2所示。

圖2 傳統(tǒng)的預(yù)測控制結(jié)構(gòu)圖

然而,研究被控對象(真菌生長量)的狀態(tài)變量(菌落直徑)并沒有參考軌跡,因為最理想的真菌生長量為0,只能盡量使?fàn)顟B(tài)變量的增量(生長速率)最小化;同時,考慮控制成本約束,由于控制變量(溫度和水分)存在參考軌跡(自然條件下的溫度和水分),只能調(diào)節(jié)溫度和水分盡量跟蹤其參考軌跡,使控制變量的跟蹤偏差(控制成本)最小化。

對傳統(tǒng)預(yù)測控制模型進行改進,將菌落直徑dl作為控制對象的狀態(tài)變量,將溫度和水分對真菌生長的影響τ(Tl)和ρ(Mel)作為控制變量。首先對未來時刻狀態(tài)變量(菌落直徑dl+1)進行預(yù)測,再與當(dāng)前時刻狀態(tài)變量(菌落直徑dt)進行比較得到狀態(tài)變量的增量(生長速率);同時將當(dāng)前時刻控制變量(溫度和水分影響τ(Tl)和ρ(Mel))反饋到輸入端,與參考軌跡(自然條件下溫度和水分影響τ(Tlr)和ρ(Melr))進行比較得到控制變量的跟蹤偏差(控制成本)。改進的預(yù)測控制模型框圖如圖3所示。

圖3 改進后的預(yù)測控制結(jié)構(gòu)圖

2.2 真菌生長的優(yōu)化控制目標(biāo)函數(shù)的建立

研究控制的優(yōu)化目的是通過控制器的優(yōu)化計算,達到狀態(tài)變量的增量(生長速率)和控制變量的跟蹤偏差(控制成本)的整體最小化,從而實現(xiàn)更精準(zhǔn)高效的控制。

(23)

式中:Jl為系統(tǒng)目標(biāo)函數(shù)值;k1、k2、k3為各分量的權(quán)重(k1、k2、k3都為正數(shù));Tlr與Melr分別為第l次采樣的自然條件下的溫度、水分,則τ(Tlr)與ρ(Melr)為自然條件下的溫度、水分對菌落直徑生長速率的影響。

Jl=K1τ(Tl)ρ(Mel)+K2(τ(Tl)-τ(Tlr))2+K3(ρ(Mel)-ρ(Melr))2

(24)

式中:K1、K2、K3的值可根據(jù)自然條件下的溫度、水分取值范圍給出,如式(25)所示。

(25)

式中:Tmaxr、Memaxr和Tminr、Meminr分別為自然條件下最高和最低的溫度、水分。τ(Tmaxr)、ρ(Memaxr)和τ(Tminr)、ρ(Meminr)分別為自然條件下最高和最低的溫度、水分對菌落直徑生長速率的影響。

2.3 真菌生長的優(yōu)化控制目標(biāo)函數(shù)的求解

通過對目標(biāo)函數(shù)中2個控制變量τ(Tl)和ρ(Mel)求偏導(dǎo)得到極值點,可得目標(biāo)函數(shù)最小時所對應(yīng)的最優(yōu)解。再對最優(yōu)解τ(Tl)和ρ(Mel)采用不動點法進行迭代求解,求得最優(yōu)的溫度Tl和水分Mel。

對目標(biāo)函數(shù)進行分析,由式(24)可以看出,Jl的變量為τ(Tl)與ρ(Mel),Jl對τ(Tl)與ρ(Mel)分別求一階偏導(dǎo),如式(26)所示:

(26)

令式(26)值為0 ,求出的駐點唯一,如式(27)所示:

(27)

對Jl求二階偏導(dǎo),如式(28)所示:

(28)

由于2K2>0,只需要滿足4K2K2-K12>0,則駐點就是極小值點,同時也是最小值點。

根據(jù)計算出的最優(yōu)控制量τ(Tl)與ρ(Mel)進行計算Tl與Mel,根據(jù)τ(Tl)與ρ(Mel)的計算結(jié)果,使用不動點法進行迭代求解,迭代公式為式(29)和(30)所示。

式中:n為迭代次數(shù);Tn為溫度的第n次迭代值;Men是水分的第n次迭代值,這里需要設(shè)定一個迭代的初始值T0與Me0,取為自然儲藏環(huán)境下的溫度和水分的最低值。

以溫度的迭代計算過程為例,根據(jù)式(29),給出迭代計算方程,方程中,x為自變量,等式左邊為y1,右邊為y2,如式(31)所示。

式中:x為自變量;y1、y2分別為式(29)的左右2邊以溫度為自變量的2個方程,那么以2個方程作2條曲線的交點的橫坐標(biāo)的值就是方程的根。對方程的根進行迭代計算求解的過程為:以自然儲藏環(huán)境下溫度的最低值作為初始值,計算x0為自變量的y2(x0),即進行了一次迭代,以計算出的y2(x0)作為x1,進行新一輪的迭代,重復(fù)直到滿足精度要求。水分的迭代計算過程與溫度相同,求得控制計算后的溫度和水分后,帶入式(20),就可以得到控制后的真菌孢子數(shù)。

研究提出的溫度水分影響下糧食真菌污染風(fēng)險預(yù)測及優(yōu)化控制技術(shù)路線如圖4所示。

圖4 技術(shù)路線圖

3 實例驗證

3.1 數(shù)據(jù)集

采用稻谷儲藏真菌污染實驗數(shù)據(jù)[16]對研究方法進行驗證,該數(shù)據(jù)是稻谷在不同環(huán)境下儲藏過程實驗數(shù)據(jù),溫度分為10、15、20、25、30、35 ℃ 6個等級,水質(zhì)量分?jǐn)?shù)分為12.5%、13.4%、14.5%、15.6%、16.2%5個等級,一共分為30組,進行稻谷儲藏實驗,每10 d進行1次真菌孢子數(shù)檢測,一共180 d的數(shù)據(jù)。根據(jù)數(shù)據(jù),對固定取Tmin≤10 ℃,Topt=30～35 ℃,Tmax≈50 ℃,Me min<13.4%,Me opt=Me max=100%。

(29)

(30)

(31)

以北京地區(qū)糧堆表層溫度、水分監(jiān)測數(shù)據(jù)[24]為參考溫度、水分,從5月到次年4月共12組數(shù)據(jù),每組數(shù)據(jù)為當(dāng)月的溫度和水分。由于采樣時間間隔是1 d,因此對數(shù)據(jù)進行擴充,擴充到12×30的數(shù)據(jù)量,共360 d,大約1年的數(shù)據(jù),對數(shù)據(jù)做插值處理,使前一個月溫度是線性改變到后一個月,圖5為插值擴充前后的溫度、水分?jǐn)?shù)據(jù);以擴充后的溫度、水分?jǐn)?shù)據(jù)作為參考輸入。

圖5 用插值擴充后的溫度

3.2 溫度和水分影響下的真菌生長建模預(yù)測

3.2.1 菌落直徑生長模型固定參數(shù)的確定

取溫度、水分的固定參數(shù)為表1對應(yīng)值。對應(yīng)的τ(T)ρ(Me)在不同溫度水分下的值如圖6所示。

表1 溫度和水分的固定參數(shù)取值

圖6 τ(T)ρ(Me)在不同溫度、水分下的值

在每一個坐標(biāo)(T,Me)處都有一個數(shù)值,數(shù)值大小與生長速率的快慢成正關(guān)。在固定參數(shù)確認(rèn)后,由于τ(T)ρ(Me)的值與菌落直徑生長速率有較強的關(guān)聯(lián),所以τ(T)ρ(Me)的值的大小反映了菌落直徑生長速率,從圖3.2.1可以看出,溫度固定時,水分越高,生長速率越快;水分固定時,隨著溫度升高,生長速率先增加后減少,與理論一致。

3.2.2 真菌孢子數(shù)增長模型系數(shù)的確定

針對1.2中,對式(6)系數(shù)C表達式進行仿真驗證,令a=1,分別取系數(shù)C為σ0、σ1、σ2時,(x,y)處的真菌孢子數(shù)f(x,y)如圖7所示。

圖7 真菌孢子數(shù)增長模型系數(shù)C取不同表達式的仿真結(jié)果

圖7中,系數(shù)C為σ0的時候,由圖7a可見,菌落直徑的變化不會使菌落中心的真菌孢子數(shù)變化。系數(shù)C為σ1的時候,由圖7b可見,σ數(shù)值越大,菌落中心孢子數(shù)就越多,與圖1.2.1菌落圖像吻合。系數(shù)C為σ2的時候,由圖7c可見,隨著σ數(shù)值增加,菌落中心孢子數(shù)變化速度過快。因此,研究選取a·σ1作為真菌孢子數(shù)增長模型式(6)的系數(shù)C表達式。

為了進一步驗證所改進真菌孢子數(shù)增長模型式(6)的準(zhǔn)確性,采用本文的真實數(shù)據(jù)集,對系數(shù)C取不同表達式進行擬合結(jié)果對比,結(jié)果如表2所示。

表2 系數(shù)的擬合結(jié)果對比

結(jié)果可以看出,系數(shù)為σ1的模型擬合SSE與RMSE都是最小,確定系數(shù)最大,所以使用σ1為模型的系數(shù)C表達式。

為了驗證所建立真菌孢子數(shù)模型,即式(21)的魯棒性,對數(shù)據(jù)集添加不同大小的擾動,對每一個原數(shù)據(jù)乘以隨機數(shù),然后用該模型進行擬合,結(jié)果如表3所示。

表3 添加擾動的擬合結(jié)果對比

其中,4組數(shù)據(jù)分別為不添加擾動,±2%、±5%、±10%擾動,通過結(jié)果可以看出,添加的擾動越大,SSE與RMSE都越大,這是添加擾動后的數(shù)據(jù)方差大導(dǎo)致的。而確定系數(shù)在±2%、±5%擾動基本不變,±10%擾動也只有微小的降低,說明了模型抵抗擾動能力較強,即魯棒性強。

3.2.3真菌孢子數(shù)的預(yù)測模型

真實的數(shù)據(jù)是溫度和水分不隨時間變化的數(shù)據(jù),所以采用溫度和水分不隨時間變化的真菌孢子數(shù)預(yù)測模型,即式(21),擬合得到參數(shù)A、B。圖8為擬合模型與真實真菌孢子數(shù)數(shù)據(jù)對比圖。

圖8 擬合模型與數(shù)據(jù)對比

其中,S是每克糧食包含的真菌孢子數(shù),三維空間中的點是實驗數(shù)據(jù),坐標(biāo)(x,y,z)分別為(τ(T)ρ(Me),t,St),空間曲面是擬合了參數(shù)后的模型。待定參數(shù)為A=1.624e+05,B=1.771e+06,對擬合好參數(shù)的模型,計算St,與真實實驗數(shù)據(jù)進行比較,溫度和水分取不同值時,分組展示結(jié)果,結(jié)果如圖9～圖11所示。

圖9 Me=0.145時在不同溫度下模型計算的結(jié)果與真實數(shù)據(jù)對比圖

圖10 Me=0.156時在不同溫度下模型計算的結(jié)果與真實數(shù)據(jù)對比圖

圖11 Me=0.162時在不同溫度下模型計算的結(jié)果與真實數(shù)據(jù)對比圖

真菌孢子數(shù)增長模型與實驗數(shù)據(jù)較接近,由于實驗數(shù)據(jù)最高的溫度剛好達到Topt,所以在同一水分下,溫度越高,真菌孢子數(shù)越多;實驗數(shù)據(jù)最高的水分未達到Meopt,所以在同一溫度下,水分越高,真菌孢子數(shù)越多。模型與實驗數(shù)據(jù)接近,所以在真實的溫度和水分范圍內(nèi),都可以用所建立的模型來計算真菌孢子數(shù),達到對糧食儲藏過程真菌的預(yù)測。

由于進行擬合時,每一組數(shù)據(jù)是溫度和水分不隨時間改變,對于溫度和水分隨時間變化時,真菌孢子數(shù)預(yù)測模型,即式(20)預(yù)測的結(jié)果如圖15沒有進行控制的曲線所示,其中溫度和水分是圖5所示的擴充后的隨時間變化的溫度和水分。

3.3 溫度和水分影響下的真菌生長優(yōu)化控制

對權(quán)重參數(shù)Tminr、Tmaxr、Meminr、Memaxr的取值如表4所示。

表4 權(quán)重計算中的參數(shù)取值

根據(jù)式(25)計算后,權(quán)重的值如表5所示,此時取值滿足極值條件。

表5 權(quán)重參數(shù)

對某一時刻的Jl,由于此時的τ(Tlr)與ρ(Melr)是已知量,例如取τ(Tlr)=0.7、ρ(Melr)=0.09,Jl的值在不同τ(Tl)、ρ(Mel)的取值如圖12所示。

圖12 不同τ、ρ對應(yīng)的目標(biāo)函數(shù)的值

對于Jl,有τ(Tl)與ρ(Mel)2個變量,Jl有唯一的最小值,結(jié)論與理論推導(dǎo)結(jié)果相同。

迭代公式的計算,以溫度迭代Tn+1的計算為例,這里取為T0=1 ℃,溫度水分的參數(shù)取值,當(dāng)取T=25 ℃時,經(jīng)計算τ(T)=0.776 4,作2條曲線y1、y2,如圖13所示。

圖13 溫度迭代計算過程示意圖

雖然初始值與真實值相差較大,但是迭代收斂很快,迭代到第3次就已經(jīng)很接近真實值。經(jīng)計算,迭代次數(shù)在10次就能把誤差減小到0.1%,所以保證準(zhǔn)確度,取迭代次數(shù)為20?？刂魄昂蟮臏囟取⑺衷?60 d內(nèi)的值如圖14所示。

圖14 未經(jīng)過控制與控制后的溫度對比圖

由于所建立目標(biāo)函數(shù)有唯一的最小值,所以求出使目標(biāo)函數(shù)最小的溫度和水分后,根據(jù)是否對溫度和水分進行控制,可以得到不同控制策略下的真菌孢子數(shù)生長,分別為:不對溫度和水分進行控制,即以參考溫度和水分進行儲藏得到的每克糧食包含的真菌孢子數(shù);僅對水分進行控制,水分取經(jīng)過優(yōu)化控制計算后的水分,溫度是在參考溫度下得到的每克糧食包含的真菌孢子數(shù);僅對溫度進行控制,溫度取經(jīng)過優(yōu)化控制計算后的溫度,水分是在參考水分下得到的每克糧食包含的真菌孢子數(shù);同時對溫度水分進行控制得到的每克糧食包含的真菌孢子數(shù)。結(jié)果如圖15所示。

圖15 不同控制策略下的真菌孢子數(shù)

當(dāng)不對溫度和水分進行控制時,真菌孢子數(shù)在第360天已經(jīng)超過10×106個/g,根據(jù)文獻[25]的真菌危害分級,已經(jīng)達到了嚴(yán)重危害的程度,所以必須對溫度和水分進行控制。只對水分進行控制的時候,真菌孢子數(shù)相對減少,在第360天雖然未達到嚴(yán)重危害,但也處于危害區(qū)。若只對溫度進行控制,則危害程度減少了47%左右,相對于只對水分進行控制時,減少幅度更大。最后同時對溫度和水分進行控制時,危害程度減少得最多,減少了53%,使糧食儲藏更安全,最大程度減少了糧食儲藏時達到危害區(qū)域的時間。不同時間的真菌孢子數(shù)及速率如表6所示,分別給出第100天、第200天、第300天的真菌孢子數(shù)(×106g-1)及速率(×106g-1day-1),每克糧食每天平均減少的真菌孢子數(shù)(×104)。第2列數(shù)據(jù)為不對溫度和水分進行控制,第3列數(shù)據(jù)為僅對水分進行控制,第4列數(shù)據(jù)為僅對溫度進行控制,第5列數(shù)據(jù)為同時對溫度水分進行控制。從時間角度來看,在每一種控制策略下,真菌孢子數(shù)增長都是越來越快。從不同控制策略角度看,不進行控制時真菌孢子數(shù)及速率增長最快,控制水分其次,然后是控制溫度,增長最慢的是同時對溫度和水分進行控制。

表6 控制結(jié)果對比(真菌孢子數(shù)/速率)

4 結(jié)論

對糧食儲存過程的真菌孢子數(shù)進行預(yù)測建模,以及對溫度和水分進行控制以減少糧食儲存過程的真菌污染危害。將水分因素加入了Rosso cardinal模型,提出同時考慮溫度和水分2種因素的菌落直徑生長模型;基于二維正態(tài)分布將菌落直徑生長模型轉(zhuǎn)換為真菌孢子數(shù)增長模型,并通過對真菌孢子數(shù)增長模型進行離散化、增量化,進而提出真菌孢子數(shù)預(yù)測模型,建立溫度、水分到真菌孢子數(shù)之間的計算模型;采用控制成本與真菌孢子數(shù)共同作用的目標(biāo)函數(shù),提出改進的預(yù)測控制模型,通過極值法求出使目標(biāo)函數(shù)極小化的中間控制量τ(Tl)和ρ(Mel),最后迭代求解出控制后溫度和水分,以減小真菌生長速率。結(jié)果表明:溫度和水分共同影響真菌孢子數(shù)的生長,水分越高,真菌孢子數(shù)的生長深度越快;低溫和高溫時真菌孢子數(shù)的生長較慢,在適宜溫度時生長最快。溫度對真菌孢子數(shù)生長的影響比水分的影響更大,同時進行控制時效果更好。對于實際糧食儲藏,此方法能預(yù)測出真菌孢子數(shù),并且控制計算給出的溫度和水分有重要參考意義,以此進行儲藏時,糧食真菌孢子數(shù)得以減少。