一种粮食主要污染物风险预测及控制方法

1.本发明涉及一种粮食主要污染物风险预测及控制方法，属于食品安全与控制的交叉学科。


背景技术：

2.粮食是生存之本，我国粮食主要分为谷物、豆类和薯类。谷物是最重要的粮食，在谷物中，稻谷是我国最主要的粮食之一。粮食通过统一的大量收割后储存，不当的储存条件会导致粮食中产生大量的微生物污染物，即霉菌的大量繁殖，从而发生霉变现象，直接降低了粮食品质，过低的品质会降低食用价值，因此微生物污染物即霉菌是粮食的主要污染物。粮食严重的霉变会直接失去食用价值，所以如何控制储存环境并对粮食主要污染物进行风险预测是一个具有研究意义的问题。
3.影响粮食霉变速度的环境因素主要是温度、水分，通常情况，霉变速度与温度的关系是：存在一个霉变速度最快的温度，称之为最适合霉变的温度，低于这个温度，霉变速率逐渐降低，低温下粮食霉变速度极慢甚至不会霉变，高于这个温度同样霉变速度也是逐渐降低，即先增后减的趋势。
4.与温度有区别，霉变速度与水分的关系是：低水分下粮食霉变较慢或者不会霉变，随着水分增多，霉变速度加快，即一直增加的趋势。霉变伴随着霉菌的真菌孢子数的产生，用真菌孢子数来衡量霉变程度，真菌孢子数越多，霉变程度就越严重。
5.综上，长期低温和低水分储存环境下储存粮食，固然会减轻霉变程度，但是储存代价会较大。所以如何通过定量地计算储存环境影响因素，使粮食微生物污染物的生长速度被控制在一个可接受范围，同时代价也不会太大，是一个待解决的问题。
6.预测控制是一种控制方法，主要包括预测模型、滚动优化和反馈校正三部分。预测模型是能实现预测功能的模型，基于当前时刻的状态与未来的输入，计算未来的输出，达到预测的目的；滚动优化是在每一个预测时域，优化的是未来的有限时域，然后将当前控制信号作用与系统后，时域向前滚动，进入下一个预测时域，反复进行；反馈校正是对预测偏差的补偿，当预测有偏差时，预测不准确会对控制量的计算带来影响，反馈校正弥补了非线性和干扰等因素带来的误差，使预测更准确，进而达到更精准的控制。预测控制对于本发明的如何对粮食储存环境进行控制，以降低真菌孢子数产生的问题，提供了先预测然后优化控制的解决思路。


技术实现要素：

7.本发明的目的是为了解决对粮食储存过程中温度、水分进行定量控制，从而控制粮食的微生物污染物霉菌生长速度的问题，首先利用微生物动力学模型对主要污染物即霉菌产生的真菌孢子数进行建模定量预测，然后用控制算法进行定量计算，得到最优储存环境影响因素，降低粮食霉变的风险。
8.所述的粮食主要污染物风险预测及控制方法，具体步骤如下：
9.步骤一、基于rosso cardinal模型，扩展为以温度、水分为影响因素的微生物菌落直径生长速率模型，并确定模型中的固定参数；
10.在实际粮食储存过程中，温度是常见的粮食霉菌生长影响因素之一，此外，水分也是影响粮食产生霉菌的重要因素，并且水分可以通过通风等手段进行定量控制，所以本发明对以温度为影响因素的rosso cardinal模型进行扩展，增加水分这一影响因素，并对固定参数进行确定，从而建立在某个时刻t下固定温度、水分的单个菌落直径生长速率模型，如公式(1)所示：
11.μ(t,me)＝μ
opt
·
τ(t)
·
ρ(me)
ꢀꢀ
(1)
12.其中是t储存环境的温度，me是储存环境的水分，μ
opt
是菌落直径生长最快时候的径向生长速率，μ(t,me)是温度为t、水分为me时候的菌落直径生长速率，τ(t)是储存环境温度为t时候的温度对生长速率的影响，ρ(me)是储存环境水分为me时候的水分对生长速率的影响，ρ(me)的计算为公式(2)所示：
[0013][0014]
其中m
e min
是能够使菌落生长的最低水分，m
e max
是能够使菌落生长的最高水分，m
e opt
是菌落生长最快时的水分，m
e min
、m
e max
和m
e opt
都是固定参数。
[0015]
然后，对单个菌落直径的生长速率模型μ(t,me)进行时间的积分，得到单个菌落直径生长模型，如公式(3)所示：
[0016][0017]
记t时刻菌落直径为d
t
。
[0018]
步骤二、根据菌落分布特点，基于二维正态分布的概率密度函数将单个菌落直径生长模型转换为真菌孢子数生长模型，利用实验数据对模型中的待定参数进行拟合；
[0019]
建立了单个菌落直径生长模型后，需要对模型的待定参数进行拟合，由于实际实验中，真菌孢子数的数据相对于菌落直径生长的数据更容易获取，因此本步骤研究将单个菌落直径生长模型转换为真菌孢子数生长模型，通过真菌孢子数的数据求得菌落直径生长模型的待定参数。据实际菌落分布特点，采用二维正态分布进行合理假设，通过对二维正态分布的概率密度函数进行改进，建立单个菌落的直径与真菌孢子数之间的关系，即菌落分布模型，如公式(4)所示：
[0020][0021]
其中(x、y)是菌落的坐标，f(x,y)是坐标(x、y)处对应的真菌孢子数。a是真菌孢子数模型的待定参数，通过拟合确定。σ是二维正态分布的参数，其大小影响菌落的分布，能够体现不同大小的菌落对应的(x、y)位置处的变化趋势，记t时刻的σ值为σ
t
，t时刻的单个菌落直径为d
t
，σ
t
与d
t
的关系如公式(6)所示：
[0022]
σ
t
＝kd
t
ꢀꢀ
(5)
[0023]
其中k是比例系数，通过参数拟合得出，d
t
是单个菌落直径。
[0024]
对f(x,y)在菌落范围：x2 y2≤d
t2
内积分，即可得到单个菌落范围的全部真菌孢子
数，计算如公式(6)所示：
[0025][0026]
其中r、θ是极坐标下的两个坐标轴，即半径坐标和角坐标；
[0027]
假设在每克粮食中，平均有m个菌落，记t时刻每克粮食全部菌落的真菌孢子数为s
t
，则有
[0028][0029]
其中s
t
的单位为(个/克)，f
min
定义为菌落边缘处的真菌孢子数，简化参数，记a＝2πamμ
opt3
k3, b＝-2πmf
min
μ
opt2
k2，当温度、水分不随时间变化时，公式(7)则可表示为：
[0030]st
＝a[t
·
τ(t)ρ(me)]3 b[t
·
τ(t)ρ(me)]2ꢀꢀ
(8)
[0031]
公式(8)即为每克粮食的真菌孢子数生长模型。由于a、m、k、f
min
都是待定参数，所以a、b 是待定参数，通过拟合来得到；固定参数包含在μ(t,me)中，通过步骤一查文献得到的经验数据确定。
[0032]
步骤三、对菌落直径生长模型进行离散化和增量化，将其转换为离散模型作为菌落直径的预测模型，进行菌落直径和真菌孢子数的预测；
[0033]
在步骤二得到模型的参数后，接下来是预测控制的预测模型的建立。首先，对菌落直径生长模型进行离散化、增量化：
[0034]
离散化为公式(9)所示：
[0035][0036]
其中，i是对d
t
进行拆分的分段标号，t
t
是t时刻的温度，表明温度是变化量；m
e t
是t时刻的水分，表明水分是变化量。
[0037]
增量化为公式(10)所示：
[0038][0039]
其中d
t 1
是t 1时刻的菌落直径，δd
t 1
是从t到t 1时间内菌落直径的增量；
[0040]
将单个菌落直径生长速率模型带入，向前一步预测，得到单个菌落直径的预测模型，计算为公式(11)所示：
[0041]dt 1
＝d
t
μ
opt
τ(t
t
)ρ(m
et
)δt
ꢀꢀ
(11)
[0042]
τ(t
t
)是储存环境温度为t
t
时候的温度对菌落直径生长速率的影响，ρ(m
et
)是储存环境水分为 m
et
时候的水分对菌落直径生长速率的影响，δt是采样时间间隔。
[0043]
由于衡量粮食安全的污染物大多以真菌孢子数为主，所以需要计算该t时刻全部菌落范围的真菌孢子数s
t
，以方便进行风险分析。当温度、水分是不随时间变化的时候，s
t
的计算如公式(8) 所示，但是当温度、水分随时间改变的时候，s
t
计算为公式(12)：
[0044][0045]
基于t时刻的状态，与t时刻待计算的、即将进入系统的输入量，可以计算下一时刻的状态，以公式(11)作为单个菌落直径预测模型，公式(12)作为每克粮食包含真菌孢子数的预测模型。
[0046]
步骤四、建立粮食主要污染物优化控制的目标函数，对单个菌落直径的预测值进行控制；
[0047]
在t时刻，粮食主要污染物的优化控制目标函数如公式(13)所示：
[0048]jt
＝k1τ(t
t
)ρ(m
et
) k2(τ(t
t
)-τ(t
tr
))2 k3(ρ(m
et
)-ρ(m
etr
))2ꢀꢀ
(13)
[0049]
其中t
tr
与m
etr
分别是t时刻的自然条件下储存环境的温度、水分，即参考温度、水分，则τ(t
tr
) 与ρ(m
etr
)为参考温度、水分对菌落直径生长速率的影响。j
t
是综合考虑温度、水分和菌落直径的系统性能指标，优化的目的就是找出一组τ(t
t
)与ρ(m
et
)使j
t
极小化，从而达到使菌落直径增量、温度和水分的控制代价之和最小化，k1、k2、k3是表示性能指标各分量的权重，由于性能指标中每个权重的量度不同，量度越大时，对应的权重就越小。k1、k2、k3的初始值k
1*
、k
2*
、k
3*
根据自然储存环境的温度、水分按照计算给出，如公式(14)所示：
[0050][0051]
t
maxr
与m
e maxr
分别是自然储存环境下最高的温度和水分，t
minr
与m
e minr
分别是自然储存环境下最低的温度和水分。τ(t
maxr
)与ρ(m
e maxr
)则分别为自然储存环境下最高的温度、水分对菌落直径生长速率的影响，τ(t
minr
)与ρ(m
e minr
)分别为自然储存环境下最低的温度、水分对菌落直径生长速率的影响。
[0052]
步骤五、对粮食主要污染物优化控制目标函数进行求解，得到影响因素的最优控制量；
[0053]jt
的变量为τ(t
t
)与ρ(m
e t
)，对j
t
进行优化，即每次优化一步，每一步都需要优化计算t时刻最优的τ(t
t
)与ρ(m
e t
)。j
t
由三部分相加，第一部分对τ(t
t i
)与ρ(m
e t i
)都是增函数，第二部分对τ(t
t i
)是减函数，第三部分对ρ(m
e t i
)是减函数。对分段函数j
t
中的变量τ(t
t
)与ρ(m
e t
)分别求偏导，即可得到极值所在位置；j
t
对τ(t
t
)与ρ(m
e t
)分别求偏导，如公式(15)所示：
[0054][0055]
其中δt是采样时间间隔，令公式(15)的值为0，求出极值点如下公式(16)所示：
[0056][0057]
最后，计算出极值点之后，对控制量τ(t
t
)与ρ(m
e t
)迭代求解，求得温度、水分的具体值，即为最优控制量。
[0058]
使用不动点法进行迭代，迭代公式为(17)和(18)所示：
[0059][0060][0061]
其中n是迭代次数，tn是温度的第n次的迭代值，t
opt
是菌落生长最快时的温度，t
min
是能够使菌落生长的最低温度，t
max
是能够使菌落生长的最高温度；m
e n
是水分的第n次迭代值，m
e opt
是菌落生长最快时的水分，m
e min
是能够使菌落生长的最低水分，m
e max
是能够使菌落生长的最高水分。
[0062]
以温度迭代计算过程为例，给出温度、水分的迭代计算方法。根据公式(17)，给出迭代计算方程中，x为自变量，等式左边为y1，右边为y2，如公式(19)所示
[0063][0064]
其中x是自变量，意义是迭代计算出的温度大小。y1、y2分别是公式(17)的左右两边以温度为自变量的两个方程，那么以两个方程作两条曲线的交点的横坐标的值就是方程的根。
[0065]
对方程的根进行迭代计算求解的过程为：以自然储存环境下温度的最低值作为初始值，计算x0为自变量的y2(x0)的值，即进行了一次迭代，以计算出的y2(x0)作为x1，进行新一轮的迭代，重复下去，直到满足精度要求。
[0066]
水分的迭代计算过程与温度相同。得到最优的温度、水分后，带入公式(12)就可以得到控制后的真菌孢子数。
[0067]
本发明的优点在于：
[0068]
1.本发明提出一种新的粮食主要污染物生长模型，即以温度、水分为影响因素的微生物单个菌落直径生长速率模型，该模型是对以温度为影响因素的rosso cardinal模型的扩展，同时考虑了粮食储存的实际环境下的温度、水分两个主要影响因素对粮食主要微生物污染物生长的影响，提高了粮食主要污染物风险的预测精度，所考虑的温度和水分两个影响因素均具有易于定量控制的特点，对实现粮食主要污染物风险的控制具有现实意义。
[0069]
2.本发明提出一种粮食主要污染物单个菌落直径与真菌孢子数的转换关系模型，即菌落分布模型，真菌孢子数的数据相对于菌落直径生长的数据更容易获取，所以将单个菌落直径生长模型转换为真菌孢子数模型，通过真菌孢子数的数据求得菌落直径生长模型的待定参数。基于菌落范围内不同位置真菌孢子数的分布特点，采用二维正态分布进行合理假设，推导出单个菌落的直径与其中所含真菌孢子数的关系模型即菌落分布模型，进而得到温度、水分对真菌孢子数生长模型，然后利用实际的数据进行参数拟合，拟合效果良好，说明所提出的菌落分布模型合理。
[0070]
3.本发明提出一种基于温度、水分影响的菌落直径预测模型，预测模型是在步骤一所建立的微生物单个菌落直径生长速率模型的基础上，通过对微生物菌落直径模型的离散化、增量化处理，实现对菌落直径的预测，由于真菌孢子数的数据更方便风险评价分析，所以用步骤二所建立的菌落直径与真菌孢子数关系模型即菌落分布模型进行转换，实现对真菌孢子数的预测，从而实现对粮食主要污染物的风险预测。
[0071]
4.本发明提出一种粮食主要污染物风险优化控制的目标函数，该目标函数与传统优化控制的目标函数的区别在于：传统优化控制的目标是力求状态变量接近其参考值，即状态变量与其参考值之差的绝对值极小化，同时使控制变量随时间的增量即控制的代价极小化；而本发明中的优化控制目标是力求状态变量即菌落直径随时间的增量极小化，同时使控制变量即温度和水分与其在自然环境下的参考值之差的绝对值极小化，以求取一组粮食储存环境影响因素的最优控制量，达到菌落直径随时间的增长尽可能小，同时对温度和水分的控制代价也较小的目的。并且对目标函数里的权重进行了初值的计算，使权重参数更容易调节。
[0072]
5.本发明提出一种粮食主要污染物风险优化控制的目标函数的求解方法，通过对目标函数按预测步长进行分段拆分，并对每一段目标函数进行求解，从而求出每一步预测的最优控制量，由于预测步长取不同值的时候进行优化计算的结果相同，所以取步长为1进行求解，求取整个性能指标最小值的解析解后，使用不动点法进行迭代求取温度、水分的近似解，由于迭代收敛快，近似解可以视为实际解。
附图说明
[0073]
图1是本发明一种粮食主要污染物风险预测及控制方法流程图。
[0074]
图2是本发明在不同的分布模型参数σ下，真菌孢子数f(x,y)与菌落位置坐标(x、y)对应的值。
[0075]
图3是本发明的τ(t)ρ(me)在不同温度、水分下的值。
[0076]
图4是本发对真菌孢子数模型进行参数拟合的结果图，其中s表示真菌孢子数。
[0077]
图5是本发明水分在0.145时不同温度下模型计算的结果与真实数据对比图。
[0078]
图6是本发明水分在0.156时不同温度下模型计算的结果与真实数据对比图。
[0079]
图7是本发明水分在0.162时不同温度下模型计算的结果与真实数据对比图、。
[0080]
图8是本发明不同τ、ρ对应的目标函数j的值，其中τ是温度对生长速率的影响，ρ是水分对生长速率的影响，j是性能指标；
[0081]
图9是本发明的温度迭代计算过程示意图，其中xi,i＝1,2
…
n是第i次温度迭代计算的值。
[0082]
图10是本发明参考温度经过插值扩充的图。
[0083]
图11是本发明参考水分经过插值扩充的图。
[0084]
图12是本发明未经过优化与优化后的温度对比图。
[0085]
图13是本发明未经过优化与优化后的水分对比图。
[0086]
图14是本发明未经过优化与优化后的真菌孢子数对比图。
具体实施方式
[0087]
下面将结合附图和实施例对本发明作进一步的详细说明。
[0088]
本发明的一种粮食主要污染物风险预测及控制方法，是在粮食储存过程中，通过对存储环境的温度、水分进行计算，以尽量小的储存代价来使真菌孢子这一污染物的数量生长尽量少。本发明基于预测控制的粮食储存环境的温度、水分建模预测及控制，如图1所示：主要包括以下步骤：
[0089]
步骤一、基于rosso cardinal模型，扩展为以温度、水分为影响因素的微生物菌落直径生长速率模型，并确定模型中的固定参数；
[0090]
调研发现，粮食储存时，温度、水分都会影响微生物的生长，不同温度、水分下微生物生长速度不同，并且温度、水分都可以被控制，例如通过制冷、通风等方法来降低温度、水分。所以将温度、水分作为真菌生长的影响因素，而rosso cardinal模型仅以温度为影响因素的，所以本发明对 rosso cardinal模型扩展为温度、水分两个影响因素，提高粮食主要污染物风险的预测精度，结合粮食霉变产生真菌分布的实际形状，把真菌污染物视为多个菌落组成，首先建立单个菌落直径的生长速率模型，并且把菌落能生长的温度、水分范围外的生长速度取为零。就可以基于温度、水分，求得菌落的径向生长速度，得到菌落生长速率模型，求取菌落的直径。
[0091]
rosso cardinal模型是一个用温度来描述菌落直径生长快慢的模型，由于水分也会影响菌落生长，基于rosso cardinal模型，将其进行扩展，把水分因素加以考虑，使模型更精准，首先建立单个菌落在某个时刻t下固定温度、水分的菌落生长速率模型，对其中的参数进行逐一确定。模型如公式(1)所示：
[0092]
μ(t,me)＝μ
opt
·
τ(t)
·
ρ(me)
ꢀꢀ
(1)
[0093]
其中μ(t,me)是温度为t、水分为me时的菌落直径生长速率。t是储存环境的温度，me是储存环境的水分，μ
opt
是菌落直径生长最快时候的径向生长速率，τ(t)是储存环境温度为t时对生长速率的影响，ρ(me)是储存环境水分为me时对生长速率的影响。
[0094]
τ(t)和ρ(me)因为计算时单位相互抵消，是无单位的数值，τ(t)和ρ(me)用公式(2)和公式 (3)计算：
[0095][0096][0097]
其中t
min
、t
max
、t
opt
是三个与温度有关的固定参数，t
min
是能够使菌落生长的最低温
度，温度超过t
min
时，菌落才会生长，低于t
min
时，菌落无法生长；t
max
是能够使菌落生长的最高温度，温度低于t
max
时，菌落才会生长，高于t
max
时，菌落无法生长；t
opt
是菌落生长最快时的温度，温度在t
opt
两侧时，生长速度会衰减，并且在t
min
与t
max
处衰减到无法生长。
[0098]me min
、m
e max
、m
e opt
是三个与水分有关的固定参数，m
e min
是能够使菌落生长的最低水分，水分超过m
e min
时，菌落才会生长，低于m
e min
时，菌落无法生长；m
e max
是能够使菌落生长的最高水分，水分低于m
e max
时，菌落才会生长，高于m
e max
时，菌落无法生长；m
e opt
是菌落生长最快时的水分，水分在m
e opt
两侧时，生长速度会衰减，并且在m
e min
与m
e max
处衰减到无法生长。
[0099]
使用《稻谷储藏真菌危害早期预测的研究》(唐芳)进行实验得到的数据，进行固定参数初步确定：由实验数据得知，t＝10℃时，高水分情况下菌落仍然会生长，所以取t
min
《10℃即可； t＝30℃之前，菌落生长速率仍然处于增长趋势，但是在t＝35℃时，部分水分下生长速率开始有下降趋势，所以取t
opt
＝30℃～35℃；由于实际实验数据的实验区间没有达到t
max
，所以 t
max
》35℃，并且根据拟合结果得出t
max
≈50℃。同理，水分在12.5％时，并未检测出真菌生长，但是在13.4％时检测出了真菌生长，所以m
e min
《13.4％；实验区间的水分并未达到a
ωopt
与a
ωmax
，但是数据表明生长速率是随水分增长的，所以m
e opt
与m
e max
取到最大值100％，或者远超实际储存环境的水分值。
[0100]
建立单个菌落直径生长速率模型后，对其进行时间的积分，记t时刻菌落直径为d
t
，就可以得到单个菌落直径生长模型，如公式(4)所示：
[0101][0102]
步骤二、根据菌落分布特点，基于二维正态分布的概率密度函数将单个菌落直径生长模型转换为真菌孢子数生长模型，利用实验数据对模型中的待定参数进行拟合；
[0103]
建立了单个菌落直径生长模型后，需要对模型的待定参数进行拟合，由于实际实验中，真菌孢子数的数据相对于菌落直径生长的数据更容易获取，因此本步骤将单个菌落直径生长模型转换为真菌孢子数模型，通过真菌孢子数的数据求得单个菌落直径生长模型的待定参数。据实际菌落分布特点，采用二维正态分布进行合理假设，通过对二维正态分布的概率密度函数进行改进，建立单个菌落的直径与真菌孢子数之间的关系，即菌落分布模型，模型就转换成了温度、水分为影响因素的真菌孢子数模型，进而利用实验数据来拟合模型中的参数。
[0104]
通过观察菌落分布，菌落中间部分聚集较多，往边缘逐渐减少，类似于“沙堆”或者正态分布，直径越大的菌落，“沙堆”就越高，菌落中心的真菌孢子数就越多，并向菌落边缘逐渐减少。由于二维正态分布概率密度函数的分布与菌落分布接近，因此本发明对菌落的分布采用二维正态分布进行合理假设。二维正态分布的概率密度函数与系数μ和σ有关，μ是均值的位置，本发明以菌落中心为坐标原点；同时，由于菌落分布的模型无需进行归一化，因此将二维正态分布公式系数中对概率密度积分的归一化因子去掉，改为形状参数σ，则得到菌落分布模型，建立公式(5)所示的单个菌落直径与真菌孢子数的转换关系：
[0105]
[0106]
其中(x,y)是菌落的坐标，e是自然对数的底数，是常数；σ是二维正态分布模型的形状参数，σ的大小影响菌落的分布，能够体现不同直径大小的菌落对应的(x,y)位置处的形状变化；a是菌落分布模型的待定参数，a的值通过拟合确定，使模型更接近实验数据；f(x,y)是坐标(x,y)处对应的真菌孢子数。记t时刻的σ值为σ
t
，t时刻的单个菌落直径为d
t
，σ
t
与d
t
的关系如公式(6)所示：
[0107]
σ
t
＝kd
t
ꢀꢀ
(6)
[0108]
其中k是比例系数，通过参数拟合得出。
[0109]
此时，模型的形状与实际菌落形状相符。模型示例：取k＝1，在不同的σ
t
下，f(x,y)与(x,y)对应的值大小如图2所示，可以看出，取不同σ
t
的时候，单个菌落的孢子数的数量分布是从中心位置向两边逐渐减少，σ
t
越大，中心位置孢子数就越多，菌落范围越广，与菌落形状符合。
[0110]
对f(x,y)在菌落范围x2 y2≤d
t2
内积分即可得到单个菌落范围的全部真菌孢子数，计算如公式(7)所示：
[0111][0112]
其中r、θ是极坐标下的两个坐标轴，即半径坐标和角坐标；
[0113]
定义菌落直径：对于分布函数而言，(x,y)的位置可以取无穷远，但是对于单个菌落而言，菌落直径有大小，存在菌落边缘；定义f
min
为菌落边缘处的真菌孢子数，即在菌落边缘处满足 f(x,y)＝f
min
，f(x,y)≥f
min
时视为(x,y)在菌落的直径范围内；f(x,y)≤f
min
时，f(x,y)太小，即真菌孢子数的数量过少，视为不在菌落的直径范围内。那么有公式(8)所示关系：
[0114][0115]
化简得公式(9)所示关系：
[0116][0117]
将公式(9)带入公式(7)得到公式(10)，即得到单个菌落范围的全部真菌孢子数
[0118][0119]
再将公式(6)带入公式(10)，得到公式(11)：
[0120][0121]
接下来拟合实验数据，得到模型的待定参数：实验数据中，每一小组数据是在同一个温度、水分条件下，即t(t),me(t)在同一组数据中是定值，那么菌落直径生长速率即μ(t,me)是定值，那么在拟合模型参数的时候，公式(4)可以写为公式(12)：
[0122]dt
＝t
·
μ
opt
τ(t)ρ(me)
ꢀꢀ
(12)
[0123]
假设在每克粮食中，平均有m个菌落，记t时刻在每克粮食m个菌落的真菌孢子数为s
t
，则有s
t
的计算为公式(13)
[0124]st
＝m
·
[2πak3d
t3
(-2πf
min
k2d
t2
)]
ꢀꢀ
(13)
[0125]
那么公式(12)的单个菌落直径带入公式(13)，就可以得到t时刻每克粮食包含的真菌孢子数为公式(14)：
[0126]st
＝2πamk3[t
·
μ
opt
τ(t)ρ(me)]3 {-2πmf
min
k2[t
·
μ
opt
τ(t)ρ(me)]2}
ꢀꢀ
(14)
[0127]
其中s
t
是每克粮食在t时刻的真菌孢子数，s
t
的单位为(个/克)，其意义为：储存环境温度为t，储存环境水分为me条件下，每克粮食在t时刻的真菌孢子数s
t
，简化参数，记a＝2πamk3μ
opt3
, b＝-2πmf
min
k2μ
opt2
，公式(14)就记为公式(15)：
[0128]st
＝a[t
·
τ(t)ρ(me)]3 b[t
·
τ(t)ρ(me)]2ꢀꢀ
(15)
[0129]
公式(15)即为每克粮食的真菌孢子数生长模型，由于a、m、k、f
min
都是待定参数，所以a、b 是待定参数，通过拟合来得到；固定参数包含在μ(t,me)中，已经通过步骤一确定。
[0130]
步骤三、对菌落直径生长模型进行离散化和增量化，将其转换为离散模型作为菌落直径的预测模型，进行菌落直径和真菌孢子数的预测；
[0131]
预测控制建立离散的预测模型的作用是能够预测未来的真菌孢子数，由于步骤二的模型是连续模型，并且是同一温度、水分下的关于时间的孢子数生长模型，因此需要对模型进行离散化和增量化，使其变为离散模型，并且适用于温度水分改变，以此作为预测模型，实现对粮食主要污染物的风险预测。
[0132]
首先把菌落直径生长模型进行离散化，能够应用在温度、水分随时变化的场景。然后进行增量化，基于此时的数据，能预测未来的数据，计算两者的差值，即是增量。
[0133]
以温度水分作为系统输入，真菌孢子数作为输出，取采样周期为一天，在两次采样之间，假设了储存环境的温度、水分保持不变，基于本次采样，预测下一次采样的输出，而这个输出与两次采样之间的输入有关，最终的目的就是为了计算此输入。
[0134]
对时间的采样周期取为δt＝1(天)，在t时刻获取了粮食储存环境信息，即储存环境温度t
t
与储存环境水分a
ωt
；假设t到t 1时间内t
t
与a
ωt
会保持到下一个采样时刻，那么此期间t
t
与m
e t
是定值，同样确定了t到t 1时间内的控制信号时，控制信号会维持到下一个采样周期，那么公式(15) 里面μ(t
t
,m
e t
)是定值，公式(14)就可以写为公式(15)；
[0135]
当t
t
,m
e t
是变值的时候，菌落直径生长速率即μ(t,me)是变值，这个时候公式(4)可以写为公式(16)：
[0136][0137]
其中，i是对d
t
进行拆分的分段标号，t
t
是t时刻的温度，表明温度是变化量，m
e t
是t时刻的水分，表明水分是变化量。
[0138]
增量化得公式(17)：
[0139][0140]
其中δd
t 1
＝d
t 1-d
t
，是从t到t 1时间内的菌落直径的增量，d
t 1
是t 1时刻的菌落直径，是基于d
t
进行递推计算的，计算为公式(18)：
[0141]dt 1
＝d
t
μ(t
t
,m
e t
)
·
δt＝d
t
μ
opt
τ(t
t
)ρ(m
e t
)δt
ꢀꢀ
(18)
[0142]
其中，δt是时间间隔，t
t
和m
e t
是待计算的输入量；τ(t
t
)是储存环境温度为t
t
时对生长速率的影响，ρ(m
e t
)是储存环境水分为m
e t
时对生长速率的影响。
[0143]
以温度、水分作为系统输入，由于衡量粮食安全的污染物大多以真菌孢子数为主，所以需要计算该t时刻全部菌落范围的真菌孢子数s
t
，以方便进行风险分析。当温度、水分是不变的时候，s
t
的计算如公式(15)所示，但是当温度、水分会改变的时候，s
t
计算如下：
[0144]
将公式(16)的菌落直径带入公式(11)，以d
t
表示第t天的菌落直径，并将待定参数a、b带入，得到公式(19)：
[0145][0146]
将公式(1)所示的菌落直径生长速率模型带入公式(19)，得到公式(20)
[0147][0148]
由于μ
opt
与时间无关，取时间间隔为δt，对公式(20)化简，得到温度、水分会改变时候的s
t
计算公式，如公式(21)所示，
[0149][0150]
基于t时刻的状态，与t时刻待计算的、即将进入系统的输入量，可以计算下一时刻的状态，以公式(18)作为菌落直径预测模型，公式(21)作为真菌孢子数预测模型。由于预测模型是根据步骤二的模型推导而来，所以不需要进行反馈校正设计。
[0151]
步骤四、建立粮食主要污染物优化控制的目标函数，对单个菌落直径的预测值进行控制；
[0152]
大部分的控制系统设计的性能指标是输出量跟踪参考输入，约束是控制量的增量尽量小。传统控制是使系统的输出跟踪一个预先设定的期望轨迹，所以控制效果的衡量指标之一是期望输出yr与输出y的差，约束是控制量、控制量的变化量或者某个状态量。即性能指标一般如公式(22)所示：
[0153][0154]
其中j是整体性能指标，k1、k2、k3、k4是表示各性能指标分量的权重，u是输入量，δu是输入量的变化量，x是某个状态。是为了让输出y接近参考输出yr，是为了让输入量不要太大，是为了让输入量变化不要太大，是为了让要求的状态量不要太大。例如，在车辆轨迹跟踪控制中，这个性能指标就可以设定为：输出尽量接近参考轨迹，给油量尽量小，给油量变化尽量平缓，行驶时间尽量少。
[0155]
本发明的性能指标对输出和控制量增量的要求相反，是希望输出量的增量小，控制量贴近参考输入；通过对粮食储存环境的温度水分进行控制，希望控制代价尽量小，即控
制的储存环境温度水分尽量与自然条件下储存环境的温度水分之差尽量小，带来尽量大的收益，即真菌孢子数的增加量要小且在可接受范围内，温度、水分的控制调整是基于自然条件下储存环境的温度和水分这两个参考输入，真菌孢子数作为输出。由于性能指标的各项权重参数量纲不一样，权重大小的确定较为复杂，因此进行初步确定，给出初值。
[0156]
与一般的控制不一样，本技术的控制没有期望输出，优化步长取1，输出部分的性能指标是菌落直径的增加尽量少，即d
t 1-d
t
尽量小，由于d是随时间单调不减的，所以性能指标第一部分为公式(23)所示：
[0157]
k1(d
t 1-d
t
)
ꢀꢀ
(23)
[0158]
对于输入部分性能指标的设定，类似于一般控制的参考输出，本发明的输入是基于自然条件下的储存环境，因为储存环境的控制与自然条件下储存环境有关。控制的储存环境与自然条件下环境相差越多，代价越大，即希望控制的储存环境与自然条件下储存环境相差尽量小，由于从储存环境温度t到τ(t)、储存环境水分me到ρ(me)在自然条件下储存环境范围内是映射对应且单调递增，但是计算复杂，直接用温度、水分进行控制器的设计会带来计算上的困难，因此使用τ(t
t
)与ρ(m
e t
)作为控制器的输入。t时刻环境的参考输入记为τ(t
tr
)与ρ(m
e tr
)，那么输入部分的性能指标为公式(24) 所示：
[0159][0160]
结合公式(23)与公式(24)，在t时刻，需要优化的性能指标如公式(25)所示：
[0161][0162]
其中t
tr
与m
e tr
分别是t时刻的自然条件下储存环境的温度、水分，即参考温度、水分，则τ(t
tr
) 与ρ(m
e tr
)为参考温度、水分对菌落直径生长速率的影响。j
t
是综合考虑温度、水分和菌落直径的系统性能指标，优化的目的就是找出一组τ(t
t
)与ρ(m
e t
)使j
t
极小化，从而达到使菌落直径增量、温度和水分的控制代价之和最小化的目的，k1、k2、k3是表示性能指标各分量的权重。从输入输出角度来看，τ(t
tr
)、ρ(m
e tr
)是t时刻的参考输入，τ(t
t
)、ρ(m
e t
)是t时刻待计算的最优输入。
[0163]
输入的值越小，输出的增量也越小，虽然输出部分的性能指标减小了，但是输入部分的性能指标会增大，优化的目的是找出一组τ(t
t
)、ρ(m
e t
)，使公式(25)最小，这一组输入称为最优的输入，平衡未来的输入量与未来输出量之间的关系，希望以尽量接近参考输入的控制输入，即代价尽量小，来让输出的增量尽量小，即产生的收益尽量大。具体在粮食储存来说，就同时使控制变量即温度和水分与其在自然储存环境下的参考值之差的绝对值极小化，以求取一组粮食储存环境影响因素的最优控制量，达到菌落直径随时间的增长尽可能小，同时对温度和水分的控制代价也较小的目的。
[0164]
将公式(18)带入公式(25)，把τ(t
t
)、ρ(m
e t
)以外的参数重新整理，令k1＝k1μ
opt
δt、t、得公式(26)：
[0165]jt
＝k1τ(t
t
)ρ(m
e t
) k2(τ(t
t
)-τ(t
tr
))2 k3(ρ(m
e t
)-ρ(m
e tr
))2ꢀꢀ
(26)
[0166]
对目标函数中的权重参数进行初步确定：由于每个权重对应数值的量度不一样，所以权重参数 k1、k2、k3的取值也不一样，对应数值量度越大时，对应的权重就越小。k1、k2、k3的初始值 k
1*
、k
2*
、k
3*
可以根据自然条件下储存环境的温度、水分按照计算给出，即如公
式(27)所示：
[0167][0168]
其中t
maxr
与m
e maxr
分别是自然储存环境下最高的温度和水分，t
minr
与m
e minr
分别是自然储存环境下最低的温度和水分。τ(t
maxr
)与ρ(m
e maxr
)则分别为自然储存环境下最高的温度、水分对菌落直径生长速率的影响，τ(t
minr
)与ρ(m
e minr
)分别为自然储存环境下最低的温度、水分对菌落直径生长速率的影响。τ(t
minr
)与ρ(m
e minr
)分别为自然储存环境下最低的温度、水分对应的对菌落直径生长速率的影响。先使用公式(27)计算出的权重初始值进行控制，可以根据结果进行调整权重。
[0169]
调整规则为：当温度控制后与参考温度平均相差超过10℃后，对温度的调整过大，可以减小k2，减小温度对菌落直径生长速率的影响的权重，当温度控制后与参考温度平均相差小于2℃后，可以增大k2，增加温度对菌落直径生长速率的影响的权重。同理，当水分控制后与参考水分平均相差超过10％后，对水分的调整过大，可以减小k3，减小水分对菌落直径生长速率的影响的权重，当水分控制后与参考水分平均相差小于2％后，可以增大k3，增加水分对菌落直径生长速率的影响的权重。当温度、水分与参考的温度水分同时控制得过大后，可以减小k1，当温度、水分与参考的温度、水分同时控制得过小后，可以增大k1。
[0170]
步骤五、对粮食主要污染物优化控制目标函数进行求解，得到影响因素的最优控制量；
[0171]
由于性能指标中的增减性明确，其极值点同时就是最小值点，性能指标对两个变量进行求偏导，求出极值点，即是性能指标最小时所对应的点。此时求出的是温度、水分对应的参数，然后对参数进行迭代求解，求得温度、水分的具体值，即是将要执行的控制量。
[0172]
对性能指标进行分析，降低储存环境的温度、水分，代价会增加，但是未来的真菌孢子数会减少，其中，温度与水分共同作为输入，与真菌孢子数是成相反的增长关系。以自然储存环境下温度、水分对应参数作为参考，计算自然储存环境的温度、水分与控制的储存环境的温度、水分对应参数之差，差值较大时，类似于空调的室内外温差较大，相应控制温度、水分两部分代价也会增大，但是同时带来的收益，即孢子数生长减小了。寻求一个折中的值，即计算最优的控制量使性能指标最小。
[0173]
由公式(26)可以看出，j
t
的变量为τ(t
t
)与ρ(m
e t
)，对j
t
进行优化，即每次优化一步，每一步都需要优化计算t时刻最优的τ(t
t
)与ρ(m
e t
)；j
t
是由三部分相加组成，变量是τ(t
t
)与ρ(m
e t
)，其中第一部分对τ(t
t
)与ρ(m
e t
)都是增函数，第二部分对τ(t
t i
)是减函数，第三部分对ρ(m
e t i
)是减函数，从公式(26)可以看出，j
t
必有唯一最小值，同时也是极小值，优化计算的目的，就是找出使j
t
最小时候的τ(t
t
)与ρ(m
e t
)，对于每一个时刻的优化，都有一个相同形式的性能指标，都会求取这一步对应最小的j
t
，求取解析的方法是，j
t
对τ(t
t
)与ρ(m
e t
)分别求偏导，即可求得极值所在位置；
[0174]jt
对τ(t
t
)与ρ(m
e t
)分别求偏导，如公式(28)所示：
[0175][0176]
令公式(28)值为0，求出极值点如公式(29)所示
[0177][0178]
计算出极值点之后，控制部分设计完成；
[0179]
然后根据计算出的最优控制量τ(t
t
)与ρ(m
e t
)进行计算t
t
与m
e t
，根据τ(t
t
)与ρ(m
e t
)的计算结果，使用不动点法进行迭代求解，迭代公式为公式(30)和(31)所示：
[0180][0181][0182]
其中n是迭代次数，tn是温度的第n次迭代值，m
e n
是水分的第n次迭代值，这里需要设定一个迭代的初始值t0与m
e 0
，可以取为自然储存环境下的温度和水分的最低值。
[0183]
以温度的迭代计算过程为例，根据公式(30)，给出迭代计算方程，方程中，x为自变量，等式左边为y1，右边为y2，如公式(32)所示
[0184][0185]
其中x是自变量，意义是迭代计算出的温度大小。y1、y2分别是公式(30)的左右两边以温度为自变量的两个方程，那么以两个方程作两条曲线的交点的横坐标的值就是方程的根。对方程的根进行迭代计算求解的过程为：以自然储存环境下温度的最低值作为初始值，计算x0为自变量的 y2(x0)的值，即进行了一次迭代，以计算出的y2(x0)作为x1，进行新一轮的迭代，重复下去，直到满足精度要求。
[0186]
水分的迭代计算过程与温度相同，求得控制计算后的温度和水分后，带入公式(21)，就可以得到控制后的真菌孢子数。
[0187]
实施例1：
[0188]
以稻谷储藏真菌危害早期预测的实验数据作为拟合的数据，《稻谷储藏真菌危害早期预测的研究》(唐芳)是模拟稻谷在不同环境下进行储存，温度是10℃、15℃、20℃、25℃、30℃、35℃六个等级，水分是12.5％、13.4％、14.5％、15.6％、16.2％五个等级，一共30个分组，进行模拟储存，每10天进行一次真菌孢子数检测，一共180天的数据。
[0207]
步骤五、对粮食主要污染物优化控制目标函数进行求解，得到影响因素的最优控制量；
[0208]
对每一时刻的j
t
，用步骤四的方法经过简单计算后，取参数数值如表4所示，j
t
的值与不同τ、ρ的取值如图8所示：
[0209]
表4性能指标的取值示意
[0210][0211]
可以看出，当性能指标的参数取表4所示的值时，j
t
的取值，受τ(t
t
)与ρ(m
e t
)两个变量的影响，且j
t
有唯一的最小值，结论与理论推导结果相同。
[0212]
迭代公式的计算，以温度迭代t
n 1
的计算为例，这里取为t0＝1℃，温度水分的固定参数取值如表1所示，当取t＝25℃时，经计算τ(t)＝0.7764，作两条曲线y1、y2如图9所示，可以看出，虽然初始值与真实值相差较大，但是迭代收敛很快，迭代到第三次就已经很接近了，收敛速度很快。经计算，迭代次数在10次就能把误差减小到0.1％，所以保证准确度，取迭代次数为20。
[0213]
模拟储存环境，得到参考温度、水分。以实际粮食储存环境检测的温度和水分数据为参考的温度和水分，数据来源于《北京地区粮堆表层温度、水分、气体监测数据分析》(祁智慧)，从5月到次年4月一共12个数据，每个数据为当月的水分，并以外界温度作为参考温度，同时也是12个数据，每个数据是当月的温度。由于采样时间间隔是1天，因此对数据进行扩充，扩充到12*30一共 360天，大约一年的数据，使前一个月温度是线性改变到后一个月，那么对数据做插值处理，图10 和图11分别为插值扩充前后的温度、水分数据；以扩充后的温度、水分数据为参考输入。
[0214]
优化公式(29)计算后，优化前后的温度、水分和真菌孢子数在360天内的值分别如图12、图13和图14所示。