Anan Tang

Philippe Castagliola

Xuelong Hu

Fupeng Xie

An assessment for the conditional performance of an SVDD-based chart

Keywords: Support Vector Data Description, Control Chart, Phase-I, Phase-II, Conditional Performance

An assessment for the conditional performance of an SVDD-based chart

Introduction

Control charts as Statistical Process Monitoring (SPM) schemes play a crucial role in a variety of manufacturing and service fields, like, for example, in monitoring the non-conforming items in the production industry, the diseases outbreak in public-health surveillance, or the climate changes in meteorology. Some of these areas have drawn significant attention of many researchers, see [START_REF] Woodall | The use of control charts in health-care and public-health surveillance[END_REF], [START_REF] Chen | Developing control charts in monitoring service quality based on the number of customer complaints[END_REF], [START_REF] Tanveer | Intelligent baby behavior monitoring using embedded vision in IoT for smart healthcare centers[END_REF]. In all of these situations, at least two essential characteristics are primarily recorded and the use of univariate control charts to monitor these quality characteristics separately is known to be misleading. Therefore, various multivariate schemes have been proposed for the simultaneous monitoring of correlated variables. The most frequently used multivariate charts, such as the Hotelling T 2 chart, the Multivariate Cumulative Sum (MCUSUM) chart and the Multivariate Exponentially Weighted Moving Average (MEWMA) chart, are generally confined to the monitoring of the mean vector of multivariate normal processes.

However, in many applications, when the multivariate normality assumption is often questionable or the actual distribution is unknown, these parametric schemes may potentially be (highly) affected. For univariate processes, several studies have shown that a departure from normality severely deteriorates the monitoring properties of schemes based on the normal theory, see [START_REF] Graham | Design and implementation issues for a class of distribution-free Phase-II EWMA exceedance control charts[END_REF], [START_REF] Asghari | Sign control chart based on ranked set sampling[END_REF], [START_REF] Hu | Guaranteed in-control performance of the synthetic chart with estimated parameters[END_REF], Tang et al. (2019), [START_REF] Alevizakos | A nonparametric double generally weighted moving average signed-rank control chart for monitoring process location[END_REF]. This situation is exacerbated for multivariate processes as the multivariate normality is even more uncommon than the univariate one. Therefore, in recent years, a host of nonparametric multivariate SPM schemes have emerged as attractive alternatives to multivariate parametric ones. For example, [START_REF] Zou | A multivariate sign EWMA control chart[END_REF] developed a multivariate EWMA scheme for monitoring process location by using a sign test as the statistic. [START_REF] Zou | A spatial rank-based multivariate EWMA control chart[END_REF] further considered using the spatial ranks statistic, and they proposed a multivariate self-starting rank-based EWMA chart for monitoring process location. [START_REF] Holland | A control chart based on a nonparametric multivariate change-point model[END_REF] developed a nonparametric scheme using an approximately distributionfree multivariate Wilcoxon-Mann-Whitney test. [START_REF] Mahmoud | The Performance of multivariate CUSUM control charts with estimated parameters[END_REF] and [START_REF] Dovoedo | Effects of parameter estimation on the multivariate distributionfree Phase II sign EWMA Chart[END_REF] investigated the effect of parameter estimation on the conditional performance of the multivariate CUSUM and EWMA charts, respectively. [START_REF] Huwang | A spatial rank-based multivariate EWMA chart for monitoring process shape matrices[END_REF] studied the properties of the nonparametric EWMA chart for monitoring the shape matrix of a multivariate process. Readers are recommended to refer to [START_REF] Qiu | Some perspectives on nonparametric statistical process control[END_REF], [START_REF] Chakraborti | Nonparametric distribution-free control charts: An updated overview and some results[END_REF] for an extensive review of the literature on nonparametric charts. However, it is essential to note that, a significant part of the nonparametric schemes mention above still have some crucial prerequisites in describing the data, see [START_REF] Ning | Improved design of kernel distance-based charts using support vector methods[END_REF], [START_REF] Singh | Measuring healthcare service quality from patients' perspective: using fuzzy AHP application[END_REF]. For example, the continuity for the process variables is a common prerequisite in most of the nonparametric schemes focusing on the signs or the signed-ranks of the observations. Unfortunately, many processes do not meet such a prerequisite, and it is therefore highly desirable to propose SPM schemes for monitoring mixed-type data process (with both categorical and continuous variables) with different marginal distributions, see [START_REF] Ning | Improved design of kernel distance-based charts using support vector methods[END_REF].

As a boundary-based scheme, the Support Vector Data Description (SVDD) method only checks the distance between the new observations and the optimal boundary generated by the so called support vectors. Thus, SVDD-based control charts can generate an adaptable in-control region regardless of how the target data are distributed. Such flexibility makes it a highly active topic in modern nonparametric developments. [START_REF] Sun | A kernel-distance-based multivariate control chart using support vector methods[END_REF] first created a kernel distance-based chart (denoted as the K chart). Then, [START_REF] Kumar | Kernel distance-based robust support vector methods and its application in developing a robust K-chart[END_REF] further explored the over-fitting issues of the K chart in the presence of outliers. [START_REF] Sukchotrat | One-class classification-based control charts for multivariate process monitoring[END_REF] developed a robust parameter design strategy for a SVDD-based chart, whose control limits are established depending on the estimated percentile of the kernel distance. [START_REF] Gani | An assessment of the kernel-distance-based multivariate control chart through an industrial application[END_REF] evaluated the performance of the K chart in a real industrial application, and the results showed that this K chart is more sensitive than the T 2 chart for detecting small shifts in the mean vector. Finally, [START_REF] Ning | Improved design of kernel distance-based charts using support vector methods[END_REF] investigated three design approaches for the parameter selection for a SVDD-based chart.

Both parametric and nonparametric SPM schemes are implemented in two stages: 1) in Phase-I, some preliminary analysis based on a training data set are performed, and in-control parameters are estimated, 2) in Phase-II, the in-control parameters obtained during the Phase-I are used to establish appropriate control limits for the monitoring of new incoming data. It is important to note that different Phase-I datasets (reference sample) will result in different parameter estimates. Consequently, the control limits and the chart's performance will vary across practitioners. An extensive amount of researches has been done on parametric charts with estimated parameters, and most of which are about the monitoring of the mean, e.g. [START_REF] Saleh | Another look at the EWMA control charts with estimated parameters[END_REF], [START_REF] Hu | Guaranteed in-control performance of the synthetic chart with estimated parameters[END_REF], the variance, e.g. Castagliola andMaravelakis (2011), Diko et al. (2017) and the median, e.g. [START_REF] Cheng | VSSI median control chart with estimated parameters and measurement errors[END_REF], Tang et al. (2019). All of them have indicated that the use of estimated parameters seriously affects the in-control properties of the parametric control charts. Although many conditional designs of parametric control charts have been reported in other journals, relatively little attention has been devoted to the design of multivariate nonparametric control charts when the process parameters are unknown with few exceptions like [START_REF] Dovoedo | Effects of parameter estimation on the multivariate distributionfree Phase II sign EWMA Chart[END_REF]. In the work of [START_REF] Dovoedo | Effects of parameter estimation on the multivariate distributionfree Phase II sign EWMA Chart[END_REF], they investigated the in-control performance of the MSEWMA sign chart in the unknown parameter case, and they found that an extremely large reference sample is needed to estimate parameters so that the chart can have identical performance as the known parameter case.

Clearly, the conditional performance topic has emerged as an important research area in the SPM literature. But, to the best of our knowledge, similar investigations have not been conducted for control charts based on the SVDD theory. The current research addresses this issue and it makes contributions in three aspects: (i) it investigates the impact of the chart's parameters on the detecting performance of the SVDD-based chart under various types of distributions, (ii) it evaluates how the number of Phase-I samples affects the Phase-II performance, and (iii) it proposes two percentile bootstrap-based methods to adjust the control limits.

The rest of this paper is organized as follows: In Section 2, the standard SVDD algorithm is presented. Section 3 introduces the corresponding SVDD-based control chart. The design of its parameters is also discussed, along with an illustrative example. In Section 4, we highlight the importance of considering the conditional performance when assessing the SVDD-based chart and two bootstrap methods are suggested to adjust the control limits in Section 5. Section 6 compares the conditional in-control (IC) and out-of-control (OOC) performance of the Hotelling T 2 and the SVDDbased charts. An real-life example is provided to show the SVDD-based chart's implementation in Section 7 and, finally, conclusions and future research directions complete the paper in Section 8.

The SVDD Algorithm

Assume that a Phase-I data set of size N contains a sequence of p-variate IC observations x i = [x i1 , x i2 , . . . , x ip ] , for i = 1, 2, . . . , N . The goal of the SVDD algorithm is to find an optimal hypersphere with a minimum volume while capturing all possible training data. Let O be the center of the hyper-sphere and R be its radius. These quantities can be obtained by solving the following constrained optimization problem

Min R 2 + C N i=1 ξ i , (1)
s.t. x i -O 2 ≤ R 2 + ξ i , (2) ξ i ≥ 0, i = 1, 2, . . . , N, (3)
where • represents the Euclidean distance, ξ i , i = 1, 2, . . . , N, are the so called slack variables and C > 0 is the penalty coefficient that controls the trade-off between the volume of the hyper-sphere and the misclassification errors.

For the optimization problem in (1), a Lagrangian function can be constructed as follows:

L = R 2 + C N i=1 ξ i - N i=1 η i (R 2 + ξ i -x i -O 2 ) - N i=1 ν i ξ i , (4)
where η i ≥ 0 and ν i ≥ 0 are the Lagrange multipliers. Setting the partial derivatives of L with respect to R, O and ξ i and equating them with zero, the following constraints are obtained:

N i=1 η i = 1, O = N i=1 η i x i , C -η i -ν i = 0, i = 1, 2, . . . , N. (5)
When substituting these constraints into (4), the optimization problem could be further simplified to

Max L = N i=1 η i (x i • x i ) - N i=1,j=1 η i η j (x i • x j ), (6)
s.t. N i=1 η i = 1, 0 ≤ η i ≤ C. (7)
Solving this Quadratic Programming (QP) allows to obtain the optimal solution for the set η = {η 1 , η 2 , . . . , η N }. Only the data points x SV corresponding to non-zero η i are needed in the description of the boundary, recorded as the Support Vector (SVs). We can find that, in (5), the center of the sphere O is expressed as a linear combination of the observations x i with weights η i , i = 1, . . . , N , and the squared-Euclidean distance (referred to as "distance" in the following) between x SV and the center O is the square of the radius

R 2 = x SV -O 2 = (x SV • x SV ) -2 N i=1 η i (x SV • x i ) + N i=1,j=1 η i η j (x i • x j ). (8)
A point z to be classified is declared to belong to the target category when the distance-function df (z) is less than or equal to R 2 , i.e.

df (z) = (z • z) -2 N i=1 η i (z • x i ) + N i=1,j=1 η i η j (x i • x j ) ≤ R 2 . ( 9
)
When the process data are linear inseparable, the kernel function method is applied to map the data into a high-dimensional feature space, so that the nonlinear problem in a low-dimensional space is transformed into the linear problem in a high-dimensional space. Replacing the inner product (x i • x j ) with a kernel function K(x i , x j ) will lead to a more flexible decision boundary, and ( 8) and ( 9) become

R 2 = K(x SV , x SV ) -2 N i=1 η i K(x SV , x i ) + N i=1,j=1 η i η j K(x i , x j ), (10)
df (z) = K(z, z) -2 N i=1 η i K(z, x i ) + N i=1,j=1 η i η j K(x i , x j ). (11)
For more details about the use of kernel functions, the interested reader can refer to [START_REF] Tax | Support vector data description[END_REF].

In the following, a Gaussian kernel function K(x i , x j ) = exp(-

x i -x j 2 s 2
) is considered, where s > 0 is the kernel window width.

3 The SVDD-based control chart [START_REF] Sun | A kernel-distance-based multivariate control chart using support vector methods[END_REF], among others, considered a SVDD-based control chart by plotting df (z) as the monitoring statistics. They selected the parameters C and s to adjust the fraction of SVs to make it as close as possible to the Type-I error α, and subsequently, they set the R 2 as the control limit. However, their simulation results showed that this design can hardly achieve a satisfactory performance due to the limited parameter selection flexibility. In this paper, a separated threshold value h is set as the control limit and, if a point z, satisfies df (z) > h, then it is considered as being OOC.

Therefore, in the design of the SVDD-based control chart (if a Gaussian kernel function is considered), there are three major parameters to be set: the window width s; the penalty coefficient C; and the control limit h. Note that, the first two parameters, s and C, are internal parameters of the support vector methodology, while h is a counterpart of the Type-I error α. One may think that the addition of the parameter h will complicate the optimization. But as we will use the (1 -α)th quantile of the distance (instead of the R 2 ) of all training points to estimate the control limit h, there will be no additional optimization calculation concerning this parameter. Furthermore, this estimation of the control limit h also enhances the SVDD-based chart's robustness for the choice of different C values (this issue will be explained later).

For SVDD-based control charts, there is still no systematic method nor guideline for choosing the parameters, s and C, efficiently. The problem is that there is no knowledge about the classification error (i.e. no label for points incorrectly classified) to constructs a standard. Here, in order to better understand their effect, we intuitively select the appropriate parameters by describing the compactness and shape of the boundary for a two-dimensional data set. We generate a bivariate gamma-shape data set using the SAT model available in [START_REF] Cheng | Simultaneously monitoring frequency and magnitude of events based on bivariate gamma distribution[END_REF], see also equation ( 13). The first N I = 200 samples are drawn from an IC process, representing the training/reference Phase-I data, and the remaining N II = 200 samples represent the testing Phase-II data.

In Phase-I, Figure 1 compares the control boundaries of the SVDD obtained for several values s ∈ {1, 2, 4, 8, 16, 32} when C = 0.8. We use "+" to represent the training samples, "•" to represent the SVs, and a solid line representing the data description boundary. It can be seen that s not only affects the radius, but it also affects the original shape of the boundary. A smaller value of s will make the description of the data more specific, it allows to obtain a closer description result, but this more flexible description boundary requires more SVs, for example, in Figure 1(a) when s = 1, almost all points are SVs. While, a larger value of s will make the shape of the boundary smoother, for example, in Figure 1(f) when s = 32, only two support vectors are required to completely determine the data description, but at the same time, we have to admit that this kind of boundaries is not representative and compact enough for describing the sample data distribution. The penalty coefficient, C, on the other hand, controls the number of points that fall outside the boundary. A larger value of C reduces the number of Phase-I points outside the boundary. As it appears in Figure 2, when C decreases, the number of SVs increases.

(Please insert Figures 1 and 2 here)

It is important to note that finding the optimal values of C and s using the reference sample without further information about the incorrect classifications is a difficult challenge. It is impossible to use the analytical method of minimizing the classification error rate. One can, however, get some insights by visualizing the scatterplot for two or three-dimensional data. Therefore, this work use a two-dimensional data set as example. At first, we seek to find the appropriate values of s and C values to get a suitable boundary that can describe the Phase-I data the best as possible. After that, the value of h is estimated to set the Type-I error value at α. A detailed procedure is summarized as follows:

1. Concerning parameter s, it is found that a relatively large value of s makes the enveloping curve of the data more spherical and insensitive to the data shift, leading to more false alarms. Therefore, a smaller value of the parameter s should be recommended to obtain a boundary having a compact shape. For example, in the light of Figure 1, s ∈ [4, 8] may be recommended.

2. Concerning parameter C, it can be seen that, as its value increases, the number of SVs tends to decrease, which is consistent with previous studies. Moreover, the (1 -α)th quantile of the distance of all training points N I can provide a robust estimation of the control limit for a fixed α value regardless of the value of C.

Compute the distance of all training points

D 1 = df (x 1 ), D 2 = df (x 2 ), . . . , D N = df (x N I ).
Given a Type-I error α and the ordered D values

D (1) < D (2) < • • • < D (N I ) , the (1 -α)th percentile of the ordered D values for all training N I points is used to estimate the control limit h = D ( N I •(1-α) )
, where . . . denotes the rounded-up integer.

Through the above steps, for the Phase-II testing samples of size N II = 200, Figure 3 draws the R 2 (using eq. ( 10)) and the control limit h values for different values of C when α = 0.05 and s = 4. It can be seen that by taking different values of C ∈ {0.2, 0.6, 1.0}, the R 2 value is changing, while the value of h based on the quantile of the distance remains stable. This reduces the sensitivity of the control limit concerning the choice of C. So we suggest C = 1 in this paper, and thus the inequality constraint 0 ≤ η i ≤ C in (7) can be relaxed. Therefore, we actually have only a single parameter s to search for with less optimisation complexity.

(Please insert Figure 3 here)

Conditional performance

It needs to be emphasised that although the SVDD-based chart pays little attention to how data are distributed, its performance strongly relies on the Phase-I training data set. [START_REF] Kumar | Kernel distance-based robust support vector methods and its application in developing a robust K-chart[END_REF] had considered the over-fitting problem with outliers existing in the Phase-I data set. In this paper, we will no further analyse this issue, but we will investigate the Phase-II performance of the SVDDbased chart when estimating the control limit h conditional on the training Phase-I samples.

The Average Run Length (ARL) is the most commonly used index to evaluate the performance of both parametric and nonparametric charts. The same ARL 0 = 1/α value is taken when the process is IC, and when the process is OOC, the smaller the ARL 1 = 1/(1 -β) value, the better the performance of control charts, where β is the Type-II error. However, as emphasized in [START_REF] Ning | Improved design of kernel distance-based charts using support vector methods[END_REF], it is impossible to set a too large value for ARL 0 (the in-control ARL) unless the number of Phase-I samples is very large. For example, if one wants to set ARL 0 = 100 then at least 1500 training samples are required! In addition, the SVDD algorithm also requires a high computational effort during the training stage. For example, it takes about 2.25 hours to train a procedure with just 3000 bivariate normal observations using Matlab on a Intel Core i7 CPU. Therefore, considering the speed of calculation and the ease of implementation, we have decided to only investigate the case of ARL 0 = 100, as it is often used in SVDD-based charting literature literature (e.g. [START_REF] Kumar | Kernel distance-based robust support vector methods and its application in developing a robust K-chart[END_REF], [START_REF] Sukchotrat | One-class classification-based control charts for multivariate process monitoring[END_REF], [START_REF] Ning | Improved design of kernel distance-based charts using support vector methods[END_REF]).

In this paper, two bivariate models considering the dependence between variables are investigated:

• A symmetric one based on the bivariate normal BN(a X , b X , a Y , b Y , ρ) distribution with probability density function defined as:

f BN (x, y|a X , b X , a Y , b Y , ρ) = 2πb X b Y 1 -ρ 2 -1 exp - 1 2(1 -ρ 2 ) (y -a X ) 2 b 2 Y - 2ρ(y -a X )(x -a Y ) b X b Y + (x -a Y ) 2 b 2 X , (12)
where a X and b X > 0 (a Y and b Y > 0) are the mean and standard deviation parameters of X and Y , respectively, and ρ ∈ (0, 1) is the correlation coefficient.

• An asymmetric one based on the bivariate gamma SAT(a X , b X , a Y , b Y , ρ) distribution with probability density function defined for x > 0, y > 0 as:

f SAT (x, y > 0|a X , b X , a Y , b Y , ρ) = x b X a X -1 y b Y a Y -1 exp - x b X + y b Y 1-ρ √ a Y /a X 1 -ρ a Y /a X a X Γ(a X ) Γ(a Y -a X ) ∞ j=0 ∞ k=0 ρ a Y /a X j+k Γ(a Y -a X + k) 1 -ρ a Y /a X 2j+k Γ(a Y + j + k) j! k! x b X j y b Y j+k , (13)
where a X > 0 and b X > 0 (a Y ≥ a X and b Y > 0) are the corresponding shape and scale parameters of Y and X, respectively, Γ(•) is the gamma function, and ρ ∈ (0, 1) is a parameter that quantifies the correlation.

In order to have an overview of the conditional IC performances of the SVDD-based chart, Table 1 presents the percentiles of the IC ARL 0 values based on several numbers of Phase-I training samples N I ∈ {50, 200, 500, 1000, 2000, 3000}. The targeted IC ARL 0 value is set to 100 (i.e. α = 0.01). According to the explanations provided in Section 3, a value of s = 8 is suggested for the SVDDbased chart. We suggest the following steps for obtaining the IC ARL 0 's:

1. Generate N I training observations x 1 , x 2 , • • • , x N I from IC processes and compute the distance of all training points D 1 = df (x 1 ), D 2 = df (x 2 ), . . . , D N = df (x N I ).
Given a Type-I error α, the control limit h = D ( N I •(1-α) ) is estimated as in Section 3.

2. Generate N II = 5000 Phase-II IC testing observations z 1 , z 2 , • • • , z N II , and then compute the statistics df (z i ), i = 1, 2, . . . , N II . Record the run length RL = i if a first signal is observed df (z i ) > h. If no OOC signal is raised, i.e. i > N II then the value RL = N II = 5000 will be recorded.

3. Repeat step 2 (10 4 times), and then average the run lengths IC ARL 0 .

4. Repeat steps 1-3 (10 4 times) in order to obtain the IC ARL 0 empirical distribution.

Observations in Table 1 show that the ARL 0 values of the SVDD-based chart vary a lot as the number of Phase-I samples changes. When a too small value of N I is selected (e.g., 50, 200), the ARL 0 of the SVDD-based charts may not reach the target value. For example, if we use a training sample size N I = 50 for the BN(10, 1, 10, 2, 0.5) distribution, 95% of the practitioners would have an ARL 0 value less than 46.2. This is far too less from the target value of 100. This happens because the SVDD method is completely data-driven and a small Phase-I sample size cannot describe the data distribution well. In this case, no matter the value of the target ARL 0 , the attained ARL 0 is always less. The same conclusion can also be found in [START_REF] Ning | Improved design of kernel distance-based charts using support vector methods[END_REF]. As N I increases up to 500 or 1000, the 50 th quantile of ARL 0 values, i.e. the median of ARL 0 (MARL 0 ), gets closer to 100. Still, in more than 50% cases, the attained ARL 0 value will be below the target one, and these charts are also associated with a large variation in the ARL 0 distribution. The standard deviations of the ARL 0 (SDARL 0 ) are recorded in the last columns. A smaller SDARL 0 will provide more confidence to the practitioners about the charts' actual IC performance. For example, if we have N I = 2000 and 3000 for the BN(10, 1, 10, 2, 0.5) distribution, the SDARL 0 values decrease to 21.6 and 16.8 , respectively. As suggested in [START_REF] Zhang | Exponential CUSUM charts with estimated control limits[END_REF], a SDARL 0 within 10% of the of ARL 0 may be reasonable. Therefore, only a larger number of Phase-I samples may guarantee that individual practitioners can obtain an ARL 0 value close to the target value.

(Please insert Table 1 here) 5 Adjusting the control limits

The SVDD-based chart can detect shifts without knowing the underlying distribution, but it requires many training samples to ensure a reasonable IC robustness for different Phase-I data sets. In the design of parametric control charts with estimated parameters, [START_REF] Saleh | Another look at the EWMA control charts with estimated parameters[END_REF], [START_REF] Hu | Guaranteed in-control performance of the synthetic chart with estimated parameters[END_REF], Tang et al. (2019), among others, suggested the use of bootstrap-type algorithms to adjust the control limits to guarantee that the conditional probability of the ARL 0 value to exceed the expected value will be close to a pre-specified constant, say (1 -)100% (e.g. = 0.1). In this paper, two bootstrap methods using the Phase-I data are discussed to overcome the problem of low attained ARL 0 due to the "between-practitioners variability".

The bootstrap approach helps in constructing nonparametric confidence intervals for the control limit by re-sampling the Phase-I data. Both bootstrap percentile and bootstrap percentile-t methodologies are used in this paper. Of course, other types of bootstrap methods, like the bootstrap BC a method or the additively corrected bootstrap-t, can also be applied to solve this problem, see [START_REF] Davison | Bootstrap methods and their application[END_REF] or [START_REF] Polansky | Stabilizing bootstrap-t confidence intervals for small samples[END_REF]. We studied their performance and observed that as the number of Phase-I samples increases, all these methods provide nearly the same results. They only differ when N I ≤ 500. Therefore, for brevity, we only include in this paper two typical types of bootstrap methods. A more explicit description of the two bootstrap procedures used in this paper are given below: Bootstrap Percentile 1. Generate l = 1, 2, . . . , B bootstrap samples x l1 , x l2 , . . . , x lN I of size N I from the Phase-I data set and compute the corresponding distance of all points D l1 = df (x l1 ), D l2 = df (x l2 ), . . . , D lN I = df (x lN I ) as a sequence of statistics from the lth bootstrap sample.

2. For each bootstrap sample l, compute the

(1 -α) percentile D * l( N I •(1-α) ) of the ordered values D * l(1) < D * l(2) < • • • < D * l(N I )
, which represents the estimated control limit ĥ * l .

Figures 4 and5 display the box plots of the IC ARL 0 distributions for the SVDD-based charts for different Phase-I sample sizes, when B = 1000, α = 0.01 and = 0.1. For brevity, we consider N I ∈ {500, 1000, 1500, 2000}. The values calculated via the two types of bootstrap adjustments are referred as to 'Adjusted-p' and 'Adjusted-t', respectively. Otherwise, they are referred as to 'Unadjusted'.

For example, Figure 4 shows that for the BN(10, 1, 10, 2, 0.5) observations, when N I = 500, more than 50% of the SVDD-based charts with unadjusted control limits would have an IC ARL 0 below 100. On the other hand, as expected, the results with adjusted limits, have a higher level of MARL 0 , but this is also not enough to ensure a reasonable variation in the ARL 0 values. Comparatively, the bootstrap percentile-t method is more appropriate than the bootstrap percentile method in this case. When N I increases up to 1000, the 'Adjusted' cases result in more than 75% of the SVDDbased charts having an ARL 0 of 100 or more. Moreover, it can be seen that increasing the Phase-I training sample size N I reduces the variation in SDARL 0 in Phase-II. At the same time, we observe that, the bootstrap percentile method has a slightly smaller SDARL 0 value than the value of the bootstrap percentile-t method when 1000 ≤ N I ≤ 2000. Therefore, the bootstrap percentile method provides a better conditional performance when the number of Phase-I samples is larger. Figure 5 exhibits similar trends as for Figure 4 but for the SAT(4, 2.5, 100, 0.1, 0.5) distribution in place of the BN(10, 1, 10, 2, 0.5) distribution.

From these results, clearly, the SVDD-based chart will require more Phase-I sample observations than parametric charts with estimated parameters. But we need to emphasize that, for the parametric charts, when the distributional assumptions are violated, it leads to unsatisfactory IC performances and also the control limit adjustments based on the parametric bootstrap method are no longer appropriate. Therefore, the SVDD-based chart could be an excellent practical alternative in such situations since it is completely data-driven.

(Please insert Figures 4 and5 here)

Performance comparison

In this Section, a study is conducted to compare the conditional performance between the SVDDbased chart and and the Hotelling T 2 chart. It is important to note that the aim of this section is not to show the superiority of the SVDD-based chart over conventional charts but to complete the comparative works in the case of estimated parameters conditional on the Phase I sample. As mentioned in [START_REF] Sun | A kernel-distance-based multivariate control chart using support vector methods[END_REF], the relationship between the conventional charts and the SVDD-based chart is not 'one replaces the other' or 'one is better than the other', but it is rather 'one complements the other'. The comparisons of these charts are are readily available when assuming that the true process parameters are known.

Evaluation of IC robustness

The same data distributions considered in the previous Section are used to compare the performance of all charts, and the IC mean vector µ 0 and the IC covariance matrix Σ 0 are listed below,

• for the bivariate normal BN(10, 1, 10, 2, 0.5):

µ 0 = µ X µ Y = 10 10 Σ 0 = σ 2 X ρσ X σ Y ρσ X σ Y σ 2 Y = 1 1 1 4 ,
• for the bivariate gamma SAT(4, 2.5, 100, 0.1, 0.5):

µ 0 = µ X µ Y = 10 10 Σ 0 = σ 2 X ρσ X σ Y ρσ X σ Y σ 2 Y = 25 2.5 2.5 1 .
Our task is to test the hypotheses {H 0 : µ = µ 0 versus H 1 : µ = µ 0 }. Here, for each distribution, we generate a data set consisting of N I = 2000 IC Phase-I training observations. The control limits of all parametric charts are designed with ARL 0 = 100 (α = 0.01) when the underlying distributions, as well as the parameters, are known, and we also adjust their control limits so that they produce a conditional IC ARL 0 that exceeds the targeted value with a certain probability with (1 -)100% (e.g. = 0.1). All results have been estimated via Monte Carlo simulation (10 4 runs).

Figure 6 presents a comparison of the box plots of the ARL 0 values for both specified distributions. As expected, the Hotelling T 2 chart has a good IC performance when the actual distribution is the bivariate normal one. On the other hand, under the bivariate gamma distribution, the Hotelling T 2 chart designed for the normal distribution is highly sensitive to normality violations, as the ARL 0 values are much smaller than the designed one. [START_REF] Human | Robustness of the EWMA control chart for individual observations[END_REF] mentioned that the IC robustness is the key to the proper design and implementation of any control charts. If charts are not IC robust, their shift detection capability in the OOC situation is almost meaningless. Clearly, in this example, the parametric charts' IC properties are greatly affected by a change in the underlying distribution and parameter estimation. So it should be considered with caution for quality practitioners to use parametric charts in such situations.

On the other hand, the SVDD-based chart has a stable IC performance for both distributions under consideration. However, it is also important to note that the unadjusted control limit can hardly guarantee the IC performance of the SVDD-based chart. The median ARL 0 is very close to the desired 100 when N I = 2000, but there is a wide "between-practitioners variability" in the attained ARL 0 . In this case, more than 50% of the charts have an IC ARL 0 value below 100. Therefore, we suggest adjusting the control limit based on the nonparametric bootstrap method proposed in Section 4 to guarantee a minimum IC performance with a prespecified probability.

(Please insert Figure 6 here)

Performance in signal detection

When the process is OOC, we are interested in detecting shifts in the process mean µ 1 = µ 0 + δ from the nominal value, assuming that the IC covariance matrix does not change. Here, for each distribution defined in Section 4, a one standard-deviation shift of size δ = [σ X , σ Y ] in the process mean has occurred, which means that the OOC mean value is µ 1 = [11,12] for the bivariate normal distribution, and is µ 1 = [15,11] for the bivariate gamma distribution.

The control limits of the Hotelling T 2 chart are redesigned so that it produces the same IC ARL 0 performance for a fair comparison of the OOC properties. From Figure 7, as it can be observed, the SVDD-based chart has a similar MARL 1 performance as the Hotelling T 2 chart when the underlying process follows the bivariate normal distribution. It is also noted that the SVDD-based chart has a slightly larger SDARL 0 value. For the bivariate gamma distribution, the Hotelling T 2 chart has more larger ARL 1 values, while the SVDD-based chart provides a better ARL 1 performance than its parametric counterparts.

It can be concluded that the Hotelling T 2 chart, which is designed for the multivariate normal distribution, could potentially be (highly) affected when the distributional assumption is violated. Although bootstrapping the Phase I sample may be performed in a nonparametric way, it requires some distributional assumption in order to guarantee the IC performance for small Phase-I samples. On the other hand, the SVDD-based chart could be an good choice since it has an acceptable IC robustness and a comparable shift detection capability.

(Please insert Figure 7 here)

An illustrative example

The real data set used in this section is freely accessible at the UC Irvine Machine Learning Repository (http://archive.ics.uci.edu/ml/datasets/SECOM) and it comprises 1567 observations (sometimes missing) on 591 variables collected from sensors of a complex modern semi-conductor manufacturing process. In the following illustrative example we will only focus on variables #31 and #32 for which 1463 IC observations are available. Figure 8(a) presents the scatter-plot for the 1463 observations of variables #31 and #32. In order to test the bivariate normality for these variables, we simply performed a Jarque-Bera test and the results clearly showed that their respective marginal distributions cannot be approximated with a normal distribution. This led us to conclude that the bivariate normality cannot be assumed for variables #31 and #32 and we could expect that the SVDD-based chart would be more robust than the Hotelling T 2 chart for monitoring this process.

To illustrate the use of the SVDD-based chart, we used the first 1200 samples as Phase-I training data, and then monitored the remaining 263 samples for evaluating the conditional performance. As discussed in Section 3, since C = 1 is fixed, the only parameter to be determined is s. The choice of s influences how well the boundary fits the data distribution. Thus, observing Figure 8(a), a value of s = 8 is suggested. Because the SVDD model is established from a small reference sample, a more proper procedure is to use an adjusted control limit based on the bootstrap percentile method described in Section 4. When the desired Type-I error α = 0.01 (ARL 0 = 100) and = 0.1 values are fixed, for each bootstrap sample l = 1, 2, . . . , 1000, the (1 -α) percentile df * l( N •(1-α) ) (x) is computed as the estimated control limit ĥ * l . Then finally, the (1 -) percentile of the bootstrap distribution ĥ * l , l = 1, 2, . . . , B is obtained as the adjusted control limit h = ĥ * ( B•(1-) ) = 0.53. The Phase-II data of the statistic df (z), together with the adjusted control limit h, are also plotted in Figure 8(b). We can see that the SVDD-based chart detect an OOC situation at the 186 time points.

(Please insert Figure 8 here)

Concluding remarks

For a long time, the important effect of different Phase-I data on nonparametric-type control charts has been neglected. In this paper, we evaluated the conditional effect of the Phase-I samples on the Phase-II efficiency of a nonparametric chart based on the SVDD theory. The results indicated that many Phase-I samples are necessary to guarantee a high probability of the IC performance to be close to the desired target. Moreover, to take the "between-practitioners variability" into consideration, two bootstrap-based approaches have been proposed to estimate the control limit in this study. Although the SVDD-based chart usually needs somewhat more Phase-I samples, they prevent performance deterioration when the underlying distributional assumption is violated. Comparative results enable us to conclude that the SVDD-based chart is a good alternative when the distribution is unknown.

We encourage for more future researches on the optimal data-dependent choice of C and s. Finally, it may be worth applying cluster analysis on the Phase-I data to partition the IC data into few clusters and to identify support vector-based region for individual clusters and consider the combined region as an in-control one. However, determining the optimal number of clusters and the corresponding combining strategies would be more complex and it is left for future researches.

Disclosure statement

No potential conflict of interest was reported by the authors.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 1234567 Figure 1: Control boundaries of the SVDD obtained for different values of s when C = 0.8 and N = 200.

  The Phase-II sample

Figure 8 :

 8 Figure 8: The SVDD-based chart applied to the morning of a semi-conductor manufacturing process

Table 1 :

 1 Percentiles of the IC ARL 0 values based on several numbers of training Phase-I samples

3. Obtain the (1 -) percentile of the bootstrap distribution ĥ * l , l = 1, 2, . . . , B as the adjusted control limit h = ĥ * ( B•(1-) ) .

Bootstrap Percentile-t 1. For original data set x 1 , x 2 , . . . ,

of the Phase-I observations as an estimated value of the control limit ĥ.

2. Generate l = 1, 2, . . . , B bootstrap samples of size N I from the Phase-I data set, and compute the corresponding ĥ * l = D * l( N I •(1-α) ) and S ĥ * l , where S ĥ * l is an estimate of the standard error of ĥ * l based on the sample l. The calculation of S ĥ * l (and S ĥ) can also be based on a bootstrap estimate.

3. For each bootstrap sample, compute the t * l = ( ĥ * l -ĥ)/S ĥ * l , l = 1, 2, . . . , B.

4. Obtain the percentile of the bootstrap distribution t * l , l = 1, 2, . . . , B, and compute the adjusted control limit h = ĥt * ( B• ) • S ĥ.


domain: spi chim phys sdv

text:
Courel Mathilde
email: mathilde.courel@agroparistech.fr

Roux Stéphanie
email: stephanie.roux@univ-montp2.fr

Birlouez-Aragon Inès

Predicting the extent of Maillard reactions in infant formula during sterilization by ohmic heating

Keywords: Nonenzymatic browning, kinetics, model liquid food product, thermal treatment

This study aimed at monitoring the biochemical reactions induced by an intense thermal treatment like sterilization of a liquid infant formula. The stake is the control of the nutritional quality of the product concerning undesirable newly formed compounds, most of them being Maillard reaction products. A fluorimetric method called FAST (Fluorescence of Advanced Maillard products and Soluble Tryptophan) was used for monitoring the accumulation of Maillard products along thermal treatment. UHT sterilization was achieved by ohmic heating in a lab scale reactor specially designed for this study. The reproducibility of the thermal treatment was characterized for the three phases of the sterilization: heating, holding and cooling. The extent of Maillard reaction was monitored under two types of conditions: (i) isothermal conditions at five holding temperatures (100, 110, 120, 130 and 140°C) and (ii) nonisothermal conditions at five heating rates (0.16 to 4.19 °C.s-1). A semi-empirical model was developed to predict product temperature history during a complete sterilization treatment with good adjustment. The extent of Maillard reaction was fitted with a kinetic model of pseudo-zero order and Arrhenius parameters were identified to model the temperature effect on the chemical reaction. The physical and the chemical models were coupled to predict the extent of Maillard reaction during ohmic treatment as a function of operating parameters (temperature and electrical conditions) and product parameters (electrical conductivity…); it was successfully adjusted to the experimental results.

INTRODUCTION

Sterilization of liquid infant formulas is obtained by means of severe heat treatments which affect the nutritional quality of the product either by degrading vitamins and part of the essential amino acids like lysine or by inducing the formation of undesirable biochemical compounds [START_REF] Birlouez-Aragon | Assessment of protein glycation markers in infant formulas[END_REF][START_REF] Contreras-Calderón | Indicators of non-enzymatic browning in the evaluation of heat damage of ingredient proteins used in manufactured infant formulas[END_REF]. Such degradations are mainly attributed to Maillard reaction leading to non-enzymatic browning within molecular reaction mechanisms still badly known. Nevertheless, these reactions are described better and better with the help of multiresponse modelling [START_REF] Martins | Unravelling the Maillard reaction network by multiresponse kinetic modelling[END_REF], an approach accounting for multiple reactions either successive or competitive which are starting from precursors originally present in the product and are strongly activated by temperature [START_REF] Van Boekel | Effect of heating on Maillard reactions in milk[END_REF]. Another level of complexity comes from the dependence of Maillard reaction towards food environment and structure (pH, minerals, nature and concentration of substrate, water activity, amorphous or crystallized state etc.). Consequently, these types of reactions are most often studied on simple model solutions at low temperatures under isothermal conditions and there is a lack of reliable kinetic data for Maillard reaction in complex food matrix thermally treated under real process conditions, for instance HTST. However, some authors studied non-enzymatic browning in food matrixes like apple juice concentrate [START_REF] Vaikousi | Kinetic modelling of non-enzymatic browning of apple juice concentrates differing in water activity under isothermal and dynamic heating conditions[END_REF] or honey [START_REF] Vaikousi | Kinetic modelling of non-enzymatic browning in honey and diluted honey systems subjected to isothermal and dynamic heating protocols[END_REF] heated under isothermal and dynamic conditions and proposed predictive browning models for fluctuating temperature-time protocols. Others introduced the concept of time-temperature integrator (TTI) based on the choice of compounds suitable for the assessment of a heating process of milk [START_REF] Claeys | From time temperature integrator kinetics to time temperature integrator tolerance levels: Heat-treated milk[END_REF]. One of the main difficulties in studying the different pathways in the Maillard reaction lies on the multiplicity of compounds that need to be identified and quantified as markers of reaction pathways. On the contrary, spectrometric methods can be of precious help to measure globally the rate of some reaction steps. Indeed, many advanced Maillard compounds have the property to be fluorescent and their spectral response at specified absorption/emission wavelength is proportional to their concentration. Such response is not specific but appears to be a reliable tool to monitor the extent of advanced Maillard reaction in different food matrixes [START_REF] Morales | A Study on Advanced Maillard Reaction in Heated Casein/Sugar Solutions: Fluorescence Accumulation[END_REF][START_REF] Birlouez-Aragon | A Rapid Fluorimetric Method to Estimate the Heat Treatment of Liquid Milk[END_REF][START_REF] Rizkallah | Front face fluorescence spectroscopy and multiway analysis for process control and NFC prediction in industrially processed cookies[END_REF]. Among them, the FAST method (Fluorescence of Advanced Maillard Products and Soluble Tryptophan), based on simple fluorescence measurements on the transparent milk supernatant obtained at pH 4.6, is well adapted for milk products. Once precisely calibrated with pertinent nutritional indicators, the FAST method appears to be of great interest for optimizing a process on the criterion of nutritional quality [START_REF] Birlouez-Aragon | The FAST method, a rapid approach of the nutritional quality of heat-treated foods[END_REF][START_REF] Desic | The FAST index -A highly sensitive indicator of the heat impact on infant formula model[END_REF]. In this study, an ohmic reactor was conceived to mimic HTST temperature conditions with characterized performances [START_REF] Roux | Design of an ohmic reactor to study the kinetics of thermal reactions in liquid products[END_REF]. The objective was to study Maillard reaction during the sterilization by UHT treatment of a model liquid infant formula. The protein denaturation and the formation of Maillard products were monitored by means of a global indicator of fluorescence provided by the FAST index. It was thus a question of modelling the evolution of the FAST index during the heat treatment by coupling a kinetic model with a model of thermal transfer.

MATERIALS & METHODS

Liquid infant formula (IF): it was specially designed for the project and produced by an industrial partner. Its general composition in g.L -1 is as follows: water 890, sugars 84, proteins 14, lipids 36 and minerals 4.6 [START_REF] Roux | Kinetics of Maillard reactions in model infant formula during UHT treatment using a static batch ohmic heater[END_REF]. Two batches of product were used for two serials of experiments A and B which are further described in the experimental design section.

FAST method [START_REF] Birlouez-Aragon | Méthode pour évaluer le traitement thermique auquel est soumis un aliment protéique tel qu'un lait[END_REF]: 1 mL of sample was completed to 50 mL with sodium acetate buffer (0.1 M, pH 4.6) to precipitate the insoluble proteins. 4 mL of supernatant was then filtered through a 0.45 µm nylon filter (VWR, France) and placed in a disposable 4 faces acryl cuve (Sarstedt, France). The fluorescence was measured with a spectrofluorimeter (CaryEclipse Varian, France) at excitation/emission wavelengths of 290/340 nm for tryptophan (F trp ) and 330/420 nm for advanced Maillard products (F amp ). The FAST index is given by the ratio (F trp /F amp )*100. Maximum variation coefficients of 3.5% were obtained for F trp , F amp and FAST index. Tryptophan fluorescence (F trp ) enables to calculate the soluble protein content. It was calibrated with soluble whey proteins (Prolacta 90 from Lactalis, Laval, France) solutions from 0.2 to 20 g.L -1 [START_REF] Roux | Kinetics of Maillard reactions in model infant formula during UHT treatment using a static batch ohmic heater[END_REF].

Ohmic reactor: a lab scale reactor was specially designed for the project, to make reproducible sterilization treatments with fully characterized thermal history [START_REF] Roux | Design of an ohmic reactor to study the kinetics of thermal reactions in liquid products[END_REF]. It is made of a 100 mL cylindrical treatment cell (internal diameter: 3.55 cm; length: 10.1 cm), closed hermetically at both extremities by two plane DSA (Dimension Stable Anodes) electrodes (Figure 1). The reactor is equipped with a sampling device composed of five independent manual valves connected to five sampling cylinders of 15 mL capacity each. Voltage (0-320 V; ± 0.25%) and current (0-20 A; ± 0.5%) are delivered by a 5 kW power source. Temperature is monitored by means of K-type thermocouples (NiCr/NiAl; diameter 1.5 mm, precision ± 1.5 °C) placed in both the cell and first sampling cylinder and it is regulated by monitoring equipment. Temperature, voltage and current are recorded at 0.2 s frequency. Nitrogen is used to maintain pressure in the treatment cell in order to reach temperature levels above 100°C (maximum 6 bars, 180°C). Reproducibility of a thermal treatment was estimated on the basis of thermal history and extent of Maillard reaction: 2.3% relative variation was obtained for the heating phase, 1% for holding and 20% for cooling; an average standard deviation of 3.8% was obtained for FAST index during a thermal treatment.

Experimental design: sterilization of the infant formula consisted in a three phase's thermal treatment, heating, holding and cooling; samples were extracted from the reactor for analysis, at different treatment times of the heating and holding phases. Two serials of experiments were conducted (Table 1), with the objective of acquiring kinetic parameters for two distinct situations: i) by sampling under isothermal conditions, i.e. during the holding phase of the thermal treatment (serial A) at five different holding temperatures (100, 110, 120, 130 and 140°C) and ii) by sampling under nonisothermal conditions, i.e. during the heating phase (serial B) at five different heating rates (0.16, 0.38, 0.86, 1.19 and 4.19 °C.s -1 ). All thermal treatments were done in triplicate whereas the sampling times were not repeated and were rather distributed so as to cover the whole time range.

RESULTS & DISCUSSION

Predicting product temperature: heating and holding phases can be modelled by a semi-empirical equation balancing heat production by Joule effect and heat losses through the cell's walls [START_REF] Roux | Design of an ohmic reactor to study the kinetics of thermal reactions in liquid products[END_REF]:

( ) c T T T T a aT dt dT - -       - + + = τ α 1 1 0 0 (1)
with: T product temperature (°C); t time (s); a product characteristic constant (°C -1 ); T 0 product initial temperature (°C); T c external temperature (20 ± 5 °C). α (s -1 ) is a simplifying constant based on process and product parameters and τ (s) is the time constant of the heat transfer model: with: U electrical tension (V); σ 0 product initial electrical conductivity (S.m -1 ); e distance between the two electrodes (m); ρ density of the product (kg.m -3 ); c p specific heat capacity (J.kg -1 .K -1 ); V product volume (m 3 ); K global exchange coefficient (W.m -2 .K -1 ); A contact surface with the product (m 2 ).

The heating equation can be obtained by integrating Eq. ( 1) with boundary conditions given by: T(0) = T 0 .

( )

      -       - - + = 1 1 exp 1 ) ( t a aT t T c τ ατ ατ ατ (4)
During holding, the absence of temperature variation allows simplifying Eq. ( 1) and extracting T(t):

      + - = c T a t T ατ ατ 1 1 ) ( (5)
These models were fitted (Microsoft Excel's Solver) to the experimental data by adjusting the time constant τ in order to minimize the sum of squares of the differences between experimental and modelled values. It was expressed as the Root Mean Square Error (RMSE) relatively to the mean value of the predicted data. The adjustment was good with RMSE around 3% for heating and 1% for holding (Figure 2) [START_REF] Roux | Design of an ohmic reactor to study the kinetics of thermal reactions in liquid products[END_REF]. The cooling phase was shown to have no contribution to the sterilization value [START_REF] Roux | Kinetics of Maillard reactions in model infant formula during UHT treatment using a static batch ohmic heater[END_REF]; it was thus considered negligible in terms of chemical reaction and it was not taken into account even if modelling is possible [START_REF] Roux | Design of an ohmic reactor to study the kinetics of thermal reactions in liquid products[END_REF]. 4) and 2 corresponds to the holding phase modelled by Eq. ( 5). Serial B experiments are centred on the beginning of the holding phase; negative times correspond to heating and positive times to holding.

Kinetic model: the extent of Maillard reaction was monitored by means of the FAST index and fitted with a kinetic model of pseudo-zero order. Arrhenius parameters were identified to model the temperature effect on the chemical reaction [START_REF] Roux | Design of an ohmic reactor to study the kinetics of thermal reactions in liquid products[END_REF]: Eq. [START_REF] Vaikousi | Kinetic modelling of non-enzymatic browning of apple juice concentrates differing in water activity under isothermal and dynamic heating conditions[END_REF] where k(T) is the reaction rate constant, k 0 the pre-exponential factor (concentration unit.s -1 ), E a the activation energy (kJ.mol -1 ), R the perfect gases constant (8.314 J.K -1 .mol -1 ) and T the temperature (K). k(T), FAST 0 , k 0 and E a can be determined as the slope and intercept of linear regressions of respectively FAST index against time and of ln(k(T)) against 1/T.

0 0 0 exp ) ( ) ( FAST t RT E k FAST t T k t FAST a + ⋅       - = + ⋅ = (
Coupling heat transfer and kinetic models: Eq. ( 4) and Eq. ( 5) were used to predict temperature during heating or holding respectively and used in Eq. ( 6) to model the extent of Maillard reaction as marked by FAST index during ohmic treatments. The thermo-kinetic model given by Eq. ( 6) is a function of operating parameters (temperature and electrical conditions) and product parameters (electrical conductivity…). By definition, FAST index is bound to protein concentration represented by the tryptophan fluorescence (F trp ). Eq. ( 6) was thus adapted according to the observation of protein behaviour. Three phases were identified (figure 3): a latency phase which lasts until the product reaches 70 °C. During the latency phase, no chemical reaction could be detected and k 0 and E a were put to zero.

A brutal decrease in protein content indicates a second phase, during which k 0 and E a were identified by fitting Eq. ( 6) to the experimental data of serial A and B experiments.

A final stable phase corresponds to isothermal conditions, where k 0 and E a parameters were identified from serial A experiments. Good adjustment of Eq. ( 6) to serial A experiments was obtained except for 130°C (probably due to experimental errors) as well as to serial B experiments (figure 4). So, the thermo-kinetic model enables to describe the heating as well as the holding phase with RMSE values lower then 15% whereas the relative error due to variability of the thermal treatment on the FAST index was reported to be maximum 3.8% [START_REF] Roux | Design of an ohmic reactor to study the kinetics of thermal reactions in liquid products[END_REF]. )

FAST exp

Eq. ( 6)

[Proteins] )

FAST exp

Eq. ( 6)

[Proteins]

CONCLUSION

The results confirm the significant promoting role of temperature in Maillard products formation. For a targeted sterilization value, high temperature short time treatments are in favour of limiting Maillard reaction development. Besides, it appears that the heating phase has noticeable contribution in Maillard product formation that cannot be neglected in kinetic studies. A simple thermo-kinetic model based on pseudo zero order kinetic of a global marker like FAST index with an Arrhenius type of temperature dependence enabled to describe a complex pool of reactions occurring in a liquid infant formula thermally treated under UHT conditions. The only required parameters are the initial thermo-physical and electrical properties of the product and the ohmic operating conditions. The dynamic phase as well as the isothermal phase of the thermal treatment by ohmic heating could be faithfully described thus allowing prediction of possible contamination by processed induced contaminants during sterilization of a liquid food. Such a predictive model could be used as a tool for optimizing the nutritional quality of a given formulation by monitoring the thermal conditions during processing.

Figure 1 .

 1 Figure 1. Diagram of the instrumented ohmic reactor [13].

Figure 2 .

 2 Figure 2. Temperature modelling of data acquired for serial A experiments (a) and serial B experiments (b). 1 corresponds to the heating phase modelled by Eq. (4) and 2 corresponds to the holding phase modelled by Eq. (5). Serial B experiments are centred on the beginning of the holding phase; negative times correspond to heating and positive times to holding.

Figure 3 .

 3 Figure 3. Kinetics of the FAST index and proteins concentration monitored under isothermal conditions (serial A) at 5holding temperatures and modelled by Eq. (6). Latency phase, Quick reaction phase and Final stable phase.

Figure 4 .

 4 Figure 4. Kinetics of the FAST index monitored under non isothermal conditions (serial B) up to 140°C at 2 extreme heating rates (a. 0.38°C.s -1 and b. 4.19°C.s -1 ) and modelled by Eq. (6). Latency phase, Quick reaction phase and Final stable phase.

Table 1 .

 1 Experimental design

	Serial	Heating phase		Holding phase
		Heating rate *	Duration	Sampling	Temperature	Sampling
		(°C.s -1 )	(min)	(°C)	(°C)	(min)
			1.15		100	0.5 -1
			1.22		110	0.25 -2
	A	1.41	1.29	/	120	0.25 -0.5
			1.36		130	0.17 -0.25
			1.42		140	0.17 -0.25
		0.16	12.5
		0.38	5.26
	B	0.86	2.33	10	140	0.5
		1.19	1.68
		4.19	0.48

* Mean heating rate calculated on the 20-140 °C domain

ACKNOWLEGEMENTS

The authors would like to thank the European Community for their financial support to ICARE project (6th Framework Programme, Grant No. COLL-CT-2005-516415).

This work was conducted within the framework of the European project ICARE (6th Framework Programme, Grant No. COLL-CT-2005-516415).