- Email: [email protected]
Evaluation of the cryptoclimatic conditions of food stores by means of a statistical dynamics method
Journal of Food Engineering
17 ( 1992) 69-82
Evaluation of the Cryptoclimatic Conditions of Food Stores by Means of a Statistical Dynamics Method M. Loueka Department
of Food Preservation Technology, Prague Institute of Chemical Technology, 166 28 Praha 6, Suchbatarova 1905. Czechoslovakia
A. HavlGek* Institute of Chemical Process Fundamentals,
Czechoslovak
Academy of Sciences
(Received 5 March 199 1; revised version received 19 July 199 1: accepted 1 August 1991)
ABSTRACT A simple model for temperature and air humidity change over long periods outside and inside food stores is presented. It is shown that for both temperature and humidity first order dynamics apply and can be quantified by means of two numerical characteristics, G and T,. Quantity G (intensification) expresses the ratio between the maximum value of traced parameter and the size of its change in location. Constant T, expresses the damping or retarding character of insulation and represents the time required for the cryptoclimatic parameter to reach 46% of its maximum value.
NOTATION
a,, b,, c, 46 c 44
C
mth coefficients of polynomials Polynomials defined by eqn (3) Estimation of C Normal distribution stochastic process
*Present address: Director, HBi Industry Ltd., PO Box 53, 160 00 Prague 6, Czechoslovakia. 69 Journal
of Food
Engineering
0260-8774/92/$05.00
Publishers Ltd, England. Printed in Great Britain
- 0
1992
Elsevier
Science
70
M. Loutka, A. HavlZek
Critical value of Suedecore distribution Gain Number of parameters in the model of m-order Likelihood function Integral variables Number of experimental points Laplace operator (h- ’ ) Time shift operator Autocorrelation function of residuals Loss function Statistical criterion defined by eqn (8) Temperature (“C) Average value of temperature (“C) Time constant (h) Input to the system Relative humidity (%) Average relative humidity Output of the system 44 0 II 0 Subscripts i 0 rn I1
Residuum of amplification Vector of parameters given by eqn (6) Variable of likelihood function Standard deviation
Input value Output value mth order or quantity nth order of quantity
INTRODUCTION Assessing the air conditions in a storage space during transport is an important technical problem. The transport container is often exposed to significant changes in conditions and carriers are interested in the variation with time of conditions inside the chamber and the influence of external conditions. A solution to this problem is proposed by exact methods based on heat and mass transport equations. However, experimentation is unavoidable despite agreement between the predicted and measured conditions because the parameters of the transport equations must be
Cryptoclimatic conditions offood stores
71
determined for individual cases. Measurement of the changes in air temperature and humidity were possible using common meteorological devices (thermohydrographs). Such measurements were carried out over a long period (roughly 1 year) inside the storage space and in the atmosphere surrounding the transport system, and were presented in graphical form. The treatment of these measurements and their interpretation is as described below. The present study is an application of the theoretical process proposed by Loueka (198 1) for the description of the dynamics of the microclimate inside food storage spaces. Total automation of the data collection was considered desirable due to the fact that considerable amounts of graphical recordings would otherwise be produced. It is presumed that the use of measurement equipment other than thermohydrographs would not be difficult in the given conditions.
FORMULATION OF THE PROBLEM: THE MATHEMATICAL MODEL AND ITS IDENTIFICATION PROCEDURE Consider the mathematical description of the system based on the experimental data collected inside and outside the system. Both the input and output variables are a time-dependent series of temperature and humidity. A simple structure of this system is shown in Fig. 1. An example of the input and output temperature signals are shown in Fig. 2. From the series of well known off line identification techniques, the technique of maximal likelihood was selected to describe the data. This has been reported by Astriim et al. ( 1965) and Gustavsson (1969). This identification technique belongs to the input-output analysis techniques and it can be described briefly as follows. The method of maximal likelihood is used to characterize linear systems with constant time parameters with one output, described by a discrete (e.g. sampling) statistical system: A(q-l)y(t)
=B(q-‘)u(t)
+lC(qwl)e(t)
(1)
where u(t) is the input to the system, y(t) is the output and e(t) is the sequence of stochastic independent quantities with a normal distribution N(0, 1). It is supposed that the variable e(t) is independent of input u(t) and q is feedforward operator with the meaning: y(t).q=y(t+l)
or
y(t).q-I
=y(t-1)
(2)
M. Loutka, A. HavliEek
72
Fig. 1.
Fig. 2.
The model of the system.
An example of a graphical printout (temperature
The polynomials A(q -I), B(q ~ ‘) and following equations:
C(q-
signal).
‘) are defined
by the
A(q-‘)=l+a,qP’+...+a,q-” B(q-‘)=
C(q-‘)=l
b,q-‘+ +c,q-‘+
. . . +b,q-”
(3)
. ..+c.q-”
where n is the order of the model. The method for the computation parameters begins by maximizing the function:
of
(4) For the given set of inputs and outputs of length N, residuals a defined recursively: Qq-‘)E(t)=a(q-‘)y(t)-B(q-‘)u(t)
where ’ ^ ’ means estimation. vector of parameter 0:
The likelihood
are (5)
function is defined by the
@=(a I,..., a,,b, )...) b,,c I,..., C”)
(6)
With respect to 0 and A it is possible to maximize the probability of function L(0, I) being correct. According to Astrom et al. (1965), it is possible to find the maximum value of L(0, A) for such a 0, which
73
Cryptoclimatic conditions offood stores
minimizes the loss function S( 0 ): S(0
)=;i E2(l)
(7)
I
The problem of defining the coefficient vector is equivalent to minimizing a function of several variables. The following test criterion can be used for testing whether the loss function decrease is significant when the order of the model is changed; r m, II
_& -Sm N-k &l km-k
where m > n
II, m are orders of models, N is the number of inputs and outputs couples, and k,, k, are numbers of parameters in the model of n, m orders, respectively. It can be shown that the stochastic variable r,,,, has F( N - k,, k, - k,) distributions for large Nvalues. The zero hypothesis can be defined as: H,:a,+,=
. ..=a.=&+,=
. . . =b,=c,+,=
. ..c.=o
(9)
At the 95% confidence level, statistical tables of the F (m, n) distribution give F = 2.7 for 100 degrees of freedom and n = 3. Another test, acceptable for the determination of the order of model, is the so-called Akaike’s informative criterion where the tested variable AIC is: (IO) where 1 is the estimation of standard deviation of residuals. According to Akaike (1974), this criterion has a minimum value at the correct value of the order of the model. The procedure for converting a continuous first order differential equation into a discrete one is described in the literature (Cypkin, 1963). If the discrete equation for the transfer function is:
T,K’) _ bq-’ 1-a,q-' TV) then the continuous is:
transfer function (Laplace transform-operator
(11) with
p)
T,(P) T'(P)
=-
G l+PT,
(1211
M. LouEka, A. Havlitek
74
where 7;-is the time constant and G‘ is the gain. The following relation holds here: a, =exp(At/T,)
(13)
b, =K[l -exp(At/T,)]
(14)
DATA TREATMENT
(DIGITIZING RECORDS)
0~
GRAPHICAL
A program SHIFT was used on an ECC 1033 computer FORTRAN IV language for checking the experimental data. The SHIFT program algorithm uses the following steps:
using
( 1)
averages both input and output signals; (2) computes variances and standard deviations of signals; (3) computes values of autocorrelation and crosscorrelation functions; (4) minimizes sum of squares of deviations for the given order of model and calculates its parameters; (5) calculates parameters of models given by eqn (1) on the basis of the maximum likelihood method; (6) uses statistical criterions given by eqns (8) and (10) for testing the order of the model. The sequence of input and output data was obtained from graphical record with the aid of a digitizer DGl produced in development workshop of CSAV. A couple consisting of temperature humidity values together with time was recorded and transmitted to measuring system on punched tape by a HP 9821 calculator.
RESULTS
OF IDENTIFICATION
the the or the
AND DISCUSSION
A sampling interval At = 2 h was used. For characterizing of both types of models N= 393 pairs of data were used. The models were identified as third order. The SHIFT program computes the coefficients of the model as well as their standard deviation and loss functions for each order of the model equation. In addition to the model equation it also computes the coefficients of the following model using the method of least squares A(qp’)y(t)=B(q-‘)u(t)+Ae(t)
(15)
75
Ctyptoclimatic conditions of food stores
Testing according to eqn (8) yields a value of 2-4 for the zero hypothesis a2 =b,
=c2
=0
(1611
of the temperature model and 4.9 for the humidity model. The zero hypothesis (16) is therefore valid only for the temperature model. While studying this criterion using a higher order moisture model, it was seen that this value increased to 56, i.e. increasing the order does not improve the criterion. The autocorrelation function of residuals computed as R(z)=
&,C I
&(i).E(i+t)
(17)
I
gives independent information. Values of R(t) are in the interval f 20, i.e. they can be considered as white (uncorrelated) noise in both types of model (Figs 3 and 4). That is why the first order was chosen as the most
1
Fig. 3.
The autocorrelation
function of residuals: air humidity model.
M. LouEku. A. Havliiek
76 1
-
60 tfme Ih)
Fig. 4.
The autocorrelation
function of residuals: air temperature
suitable in both cases. The resulting equation (input) to T, (output) is: ;r;(t+ l)-0.779T,(t)=0.224
model.
relating temperature
ZJt)+0.482@)
T
(18)
while for humidity v (input) and v, (output) the following results: ~,(t+l)-O~679~,(,(t+1)=0l65V;(t)+0~209e(t)
(19)
It is evident that the noise part of the temperature model is about twice that of the humidity model. The comparison of experimental and predicted results is given in Figs 5 and 6. Those parts of the models (without 1 e( t)) were used for the computing of the transfer characteristics (Figs 7 and 8). From models (18) and (19) it is possible to compute the statistical linearized equations valid for the operating conditions A&= 1.016
AT,-4.71
(20)
A r/;,= 0.474
A Vi+ 49.27
(21)
Oyptoclimatic conditions offood stores
0
77
.
3
78
M. LouCka, A. Havliek
Fig. 7.
Step response of the air humidity model.
20
G- O.I.743
Vft*l)-.6522
Vltl=.165
Ultl
-
time[hl
40
rltl
n ”
0.1
I
1.0136
Fig. 8.
T(tl=.ZZL
Uftl
Step response of the air temperature
20
G-l.0136
T(f*lt-,779
model.
time[hl
40
Cryptoclimaticconditions of food stores
81
The constants of eqns (20) and (21) are approximately values of both processes (z =
4.84”C and c = 6505%);
(or =
the average
5.36”C and gL1= 8.82%)
By using eqns ( 13) and ( 14) it is possible to compute the tune constant T, and gain G. Temperature: T,=8h G= 1.0136”C Humidity: T,=517
h
G = 0.4743% Numerical values of these constants illustrate very distinctly the dynamic changes in the system. When the temperature changes from __ 10°C to 20°C in the surroundings, the inside temperature eventually reaches a value of 20°C but when relative humidity changes from 60% to 80%, the humidity inside the storage or transport space never increases above 70%. The time constant indicates the time taken to reach this maximum value following step changes in outside atmosphere. Where temperature change is involved, responses are about 150% slower than those involving humidity. From the theoretical point of view, the described method can be completed with the so-called ‘cross-effect’ (Fig. 9) where the physical or real influence of the exterior temperature on the interior humidity is included. The results demonstrate that this method of assessing temperature and mass transfer is apparently suitable during storage of raw materials and products sensitive to their surroundings. This is due to the fact that an exact description of the system by means of transport equations is considered to be too complicated.
Ti ItI
------__----. . ‘. 1.
Vi it)
Fig. 9.
T,, 1 t )
. .
___________-2..
1.
. .
.A
The model with ‘cross-effect’ of temperature
volt)
on the humidity.
82
h4. Loutka, A. Havlitek
REFERENCES Akaike, H. (1974). Criteria of system identification. IEEE Trans., AC-19, 716-23. Astrom, K. J., Bohlin, B. & Weusmark, S. ( 1965). Identification of linear system. Report TP 18.150, IBM Nordic Laboratory Sweden. Cypkin, Ja. Z. ( 1963). in Theory of Linear Impulse Systems. Fizmatgiz Moscow (in Russian). Gustavsson, I. ( 1969). Identification of multivariable system. Report 6907, Department of Automatic Control, Lund Institute of Technology, Sweden. Loucka, M. ( 198 1). Theoretical foundations of dynamic behaviour of storage space with respect to random changes of outer parameters. Research report VSCHT (in Czechoslovakian).
17 ( 1992) 69-82
Evaluation of the Cryptoclimatic Conditions of Food Stores by Means of a Statistical Dynamics Method M. Loueka Department
of Food Preservation Technology, Prague Institute of Chemical Technology, 166 28 Praha 6, Suchbatarova 1905. Czechoslovakia
A. HavlGek* Institute of Chemical Process Fundamentals,
Czechoslovak
Academy of Sciences
(Received 5 March 199 1; revised version received 19 July 199 1: accepted 1 August 1991)
ABSTRACT A simple model for temperature and air humidity change over long periods outside and inside food stores is presented. It is shown that for both temperature and humidity first order dynamics apply and can be quantified by means of two numerical characteristics, G and T,. Quantity G (intensification) expresses the ratio between the maximum value of traced parameter and the size of its change in location. Constant T, expresses the damping or retarding character of insulation and represents the time required for the cryptoclimatic parameter to reach 46% of its maximum value.
NOTATION
a,, b,, c, 46 c 44
C
mth coefficients of polynomials Polynomials defined by eqn (3) Estimation of C Normal distribution stochastic process
*Present address: Director, HBi Industry Ltd., PO Box 53, 160 00 Prague 6, Czechoslovakia. 69 Journal
of Food
Engineering
0260-8774/92/$05.00
Publishers Ltd, England. Printed in Great Britain
- 0
1992
Elsevier
Science
70
M. Loutka, A. HavlZek
Critical value of Suedecore distribution Gain Number of parameters in the model of m-order Likelihood function Integral variables Number of experimental points Laplace operator (h- ’ ) Time shift operator Autocorrelation function of residuals Loss function Statistical criterion defined by eqn (8) Temperature (“C) Average value of temperature (“C) Time constant (h) Input to the system Relative humidity (%) Average relative humidity Output of the system 44 0 II 0 Subscripts i 0 rn I1
Residuum of amplification Vector of parameters given by eqn (6) Variable of likelihood function Standard deviation
Input value Output value mth order or quantity nth order of quantity
INTRODUCTION Assessing the air conditions in a storage space during transport is an important technical problem. The transport container is often exposed to significant changes in conditions and carriers are interested in the variation with time of conditions inside the chamber and the influence of external conditions. A solution to this problem is proposed by exact methods based on heat and mass transport equations. However, experimentation is unavoidable despite agreement between the predicted and measured conditions because the parameters of the transport equations must be
Cryptoclimatic conditions offood stores
71
determined for individual cases. Measurement of the changes in air temperature and humidity were possible using common meteorological devices (thermohydrographs). Such measurements were carried out over a long period (roughly 1 year) inside the storage space and in the atmosphere surrounding the transport system, and were presented in graphical form. The treatment of these measurements and their interpretation is as described below. The present study is an application of the theoretical process proposed by Loueka (198 1) for the description of the dynamics of the microclimate inside food storage spaces. Total automation of the data collection was considered desirable due to the fact that considerable amounts of graphical recordings would otherwise be produced. It is presumed that the use of measurement equipment other than thermohydrographs would not be difficult in the given conditions.
FORMULATION OF THE PROBLEM: THE MATHEMATICAL MODEL AND ITS IDENTIFICATION PROCEDURE Consider the mathematical description of the system based on the experimental data collected inside and outside the system. Both the input and output variables are a time-dependent series of temperature and humidity. A simple structure of this system is shown in Fig. 1. An example of the input and output temperature signals are shown in Fig. 2. From the series of well known off line identification techniques, the technique of maximal likelihood was selected to describe the data. This has been reported by Astriim et al. ( 1965) and Gustavsson (1969). This identification technique belongs to the input-output analysis techniques and it can be described briefly as follows. The method of maximal likelihood is used to characterize linear systems with constant time parameters with one output, described by a discrete (e.g. sampling) statistical system: A(q-l)y(t)
=B(q-‘)u(t)
+lC(qwl)e(t)
(1)
where u(t) is the input to the system, y(t) is the output and e(t) is the sequence of stochastic independent quantities with a normal distribution N(0, 1). It is supposed that the variable e(t) is independent of input u(t) and q is feedforward operator with the meaning: y(t).q=y(t+l)
or
y(t).q-I
=y(t-1)
(2)
M. Loutka, A. HavliEek
72
Fig. 1.
Fig. 2.
The model of the system.
An example of a graphical printout (temperature
The polynomials A(q -I), B(q ~ ‘) and following equations:
C(q-
signal).
‘) are defined
by the
A(q-‘)=l+a,qP’+...+a,q-” B(q-‘)=
C(q-‘)=l
b,q-‘+ +c,q-‘+
. . . +b,q-”
(3)
. ..+c.q-”
where n is the order of the model. The method for the computation parameters begins by maximizing the function:
of
(4) For the given set of inputs and outputs of length N, residuals a defined recursively: Qq-‘)E(t)=a(q-‘)y(t)-B(q-‘)u(t)
where ’ ^ ’ means estimation. vector of parameter 0:
The likelihood
are (5)
function is defined by the
@=(a I,..., a,,b, )...) b,,c I,..., C”)
(6)
With respect to 0 and A it is possible to maximize the probability of function L(0, I) being correct. According to Astrom et al. (1965), it is possible to find the maximum value of L(0, A) for such a 0, which
73
Cryptoclimatic conditions offood stores
minimizes the loss function S( 0 ): S(0
)=;i E2(l)
(7)
I
The problem of defining the coefficient vector is equivalent to minimizing a function of several variables. The following test criterion can be used for testing whether the loss function decrease is significant when the order of the model is changed; r m, II
_& -Sm N-k &l km-k
where m > n
II, m are orders of models, N is the number of inputs and outputs couples, and k,, k, are numbers of parameters in the model of n, m orders, respectively. It can be shown that the stochastic variable r,,,, has F( N - k,, k, - k,) distributions for large Nvalues. The zero hypothesis can be defined as: H,:a,+,=
. ..=a.=&+,=
. . . =b,=c,+,=
. ..c.=o
(9)
At the 95% confidence level, statistical tables of the F (m, n) distribution give F = 2.7 for 100 degrees of freedom and n = 3. Another test, acceptable for the determination of the order of model, is the so-called Akaike’s informative criterion where the tested variable AIC is: (IO) where 1 is the estimation of standard deviation of residuals. According to Akaike (1974), this criterion has a minimum value at the correct value of the order of the model. The procedure for converting a continuous first order differential equation into a discrete one is described in the literature (Cypkin, 1963). If the discrete equation for the transfer function is:
T,K’) _ bq-’ 1-a,q-' TV) then the continuous is:
transfer function (Laplace transform-operator
(11) with
p)
T,(P) T'(P)
=-
G l+PT,
(1211
M. LouEka, A. Havlitek
74
where 7;-is the time constant and G‘ is the gain. The following relation holds here: a, =exp(At/T,)
(13)
b, =K[l -exp(At/T,)]
(14)
DATA TREATMENT
(DIGITIZING RECORDS)
0~
GRAPHICAL
A program SHIFT was used on an ECC 1033 computer FORTRAN IV language for checking the experimental data. The SHIFT program algorithm uses the following steps:
using
( 1)
averages both input and output signals; (2) computes variances and standard deviations of signals; (3) computes values of autocorrelation and crosscorrelation functions; (4) minimizes sum of squares of deviations for the given order of model and calculates its parameters; (5) calculates parameters of models given by eqn (1) on the basis of the maximum likelihood method; (6) uses statistical criterions given by eqns (8) and (10) for testing the order of the model. The sequence of input and output data was obtained from graphical record with the aid of a digitizer DGl produced in development workshop of CSAV. A couple consisting of temperature humidity values together with time was recorded and transmitted to measuring system on punched tape by a HP 9821 calculator.
RESULTS
OF IDENTIFICATION
the the or the
AND DISCUSSION
A sampling interval At = 2 h was used. For characterizing of both types of models N= 393 pairs of data were used. The models were identified as third order. The SHIFT program computes the coefficients of the model as well as their standard deviation and loss functions for each order of the model equation. In addition to the model equation it also computes the coefficients of the following model using the method of least squares A(qp’)y(t)=B(q-‘)u(t)+Ae(t)
(15)
75
Ctyptoclimatic conditions of food stores
Testing according to eqn (8) yields a value of 2-4 for the zero hypothesis a2 =b,
=c2
=0
(1611
of the temperature model and 4.9 for the humidity model. The zero hypothesis (16) is therefore valid only for the temperature model. While studying this criterion using a higher order moisture model, it was seen that this value increased to 56, i.e. increasing the order does not improve the criterion. The autocorrelation function of residuals computed as R(z)=
&,C I
&(i).E(i+t)
(17)
I
gives independent information. Values of R(t) are in the interval f 20, i.e. they can be considered as white (uncorrelated) noise in both types of model (Figs 3 and 4). That is why the first order was chosen as the most
1
Fig. 3.
The autocorrelation
function of residuals: air humidity model.
M. LouEku. A. Havliiek
76 1
-
60 tfme Ih)
Fig. 4.
The autocorrelation
function of residuals: air temperature
suitable in both cases. The resulting equation (input) to T, (output) is: ;r;(t+ l)-0.779T,(t)=0.224
model.
relating temperature
ZJt)+0.482@)
T
(18)
while for humidity v (input) and v, (output) the following results: ~,(t+l)-O~679~,(,(t+1)=0l65V;(t)+0~209e(t)
(19)
It is evident that the noise part of the temperature model is about twice that of the humidity model. The comparison of experimental and predicted results is given in Figs 5 and 6. Those parts of the models (without 1 e( t)) were used for the computing of the transfer characteristics (Figs 7 and 8). From models (18) and (19) it is possible to compute the statistical linearized equations valid for the operating conditions A&= 1.016
AT,-4.71
(20)
A r/;,= 0.474
A Vi+ 49.27
(21)
Oyptoclimatic conditions offood stores
0
77
.
3
78
M. LouCka, A. Havliek
Fig. 7.
Step response of the air humidity model.
20
G- O.I.743
Vft*l)-.6522
Vltl=.165
Ultl
-
time[hl
40
rltl
n ”
0.1
I
1.0136
Fig. 8.
T(tl=.ZZL
Uftl
Step response of the air temperature
20
G-l.0136
T(f*lt-,779
model.
time[hl
40
Cryptoclimaticconditions of food stores
81
The constants of eqns (20) and (21) are approximately values of both processes (z =
4.84”C and c = 6505%);
(or =
the average
5.36”C and gL1= 8.82%)
By using eqns ( 13) and ( 14) it is possible to compute the tune constant T, and gain G. Temperature: T,=8h G= 1.0136”C Humidity: T,=517
h
G = 0.4743% Numerical values of these constants illustrate very distinctly the dynamic changes in the system. When the temperature changes from __ 10°C to 20°C in the surroundings, the inside temperature eventually reaches a value of 20°C but when relative humidity changes from 60% to 80%, the humidity inside the storage or transport space never increases above 70%. The time constant indicates the time taken to reach this maximum value following step changes in outside atmosphere. Where temperature change is involved, responses are about 150% slower than those involving humidity. From the theoretical point of view, the described method can be completed with the so-called ‘cross-effect’ (Fig. 9) where the physical or real influence of the exterior temperature on the interior humidity is included. The results demonstrate that this method of assessing temperature and mass transfer is apparently suitable during storage of raw materials and products sensitive to their surroundings. This is due to the fact that an exact description of the system by means of transport equations is considered to be too complicated.
Ti ItI
------__----. . ‘. 1.
Vi it)
Fig. 9.
T,, 1 t )
. .
___________-2..
1.
. .
.A
The model with ‘cross-effect’ of temperature
volt)
on the humidity.
82
h4. Loutka, A. Havlitek
REFERENCES Akaike, H. (1974). Criteria of system identification. IEEE Trans., AC-19, 716-23. Astrom, K. J., Bohlin, B. & Weusmark, S. ( 1965). Identification of linear system. Report TP 18.150, IBM Nordic Laboratory Sweden. Cypkin, Ja. Z. ( 1963). in Theory of Linear Impulse Systems. Fizmatgiz Moscow (in Russian). Gustavsson, I. ( 1969). Identification of multivariable system. Report 6907, Department of Automatic Control, Lund Institute of Technology, Sweden. Loucka, M. ( 198 1). Theoretical foundations of dynamic behaviour of storage space with respect to random changes of outer parameters. Research report VSCHT (in Czechoslovakian).