Breakpoint Planning Method for Rice Multibreak Milling

Excessive milling of rice kernels will result in nutrient loss and grain waste. To avoid grain waste, multibreak milling systems have been widely used in large-scale commercial rice mills. However, there is still no reasonable breakpoint planning method to guide the multibreak milling process. To construct a reasonable multibreak milling system, in this research, taking rice milling, a typical heterogeneous cereal-kernel milling process, as an example, the multivariate analysis method was used to comprehensively analyze the characteristic changes of milled rice during the whole milling process. A breakpoint planning method was established, including planning the number of breakpoints, determining the degree of milling or milling time corresponding to each breakpoint, and estimating the actual breakpoint to which the milled rice belongs. The verification results showed the rationality and high accuracy of the planning method. The presented work will help operators to plan the multibreak milling system of rice efficiently and objectively so as to significantly improve the commercial value of milled rice.


Introduction
A large part of milling operations involves the use of heterogeneous cereal kernels [1][2][3][4][5][6]. In the milling process, the structure and composition of heterogeneous milled materials can be changed markedly, such as the brown rice (BR) whitening process. Whole BR consists of four main fractions (as shown in Figure 1), namely the bran layer (pericarp layer), outer endosperm, middle endosperm (aleurone layers), and core endosperm [7]. Rice milling is a key procedure to remove the pericarp and aleurone layers from the exterior of the BR kernels, and the amount of these fractions removed is referred to as the degree of milling (DOM). In recent years, extensive research has shown that there are the considerable differences in composition and structural strength of each fraction [8,9], thus layered milling based on the microstructure of BR would be quite beneficial for enhancing the quality of the rice milled and meeting the different needs of consumers. A reasonable multibreak milling system is desirable, which enables cereal kernels to reach layered milling. Therefore, multibreak milling systems (a schematic diagram of the multibreak point milling system is shown in Figure 1) have been widely used in large-scale commercial rice mills in recent years [10,11]. However, at present, control on the DOM range for each single break and planning of the number of breakpoints are based primarily on the experience of operators and their proficiency with the processing machinery in practice [12]. As a result, the existing multibreak milling system is not only often unable to achieve layered milling but also brings huge economic losses to large-scale commercial rice mills. For example, due to excessive breakpoints, excess equipment may be idle, and the final product may be overmilled and contain a high breakage rate. On the contrary, too few breakpoints is also unreasonable because it will cause the multibreak milling system to lose the advantage of layered milling or cause the final product to fail to meet the requirements. In addition, how to control the milling time (MT) of each single-breakpoint mill to meet the requirements of layered milling is also an urgent problem to be solved. Consequently, it is of great practical value to explore the changes of milled rice during the whole milling process, and to establish an objective and highly efficient breakpoint planning method based on this.
Foods 2023, 12, x FOR PEER REVIEW 2 of 4 example, due to excessive breakpoints, excess equipment may be idle, and the final product may be overmilled and contain a high breakage rate. On the contrary, too few breakpoints is also unreasonable because it will cause the multibreak milling system to lose the advantage of layered milling or cause the final product to fail to meet the requirements. In addition, how to control the milling time (MT) of each single-breakpoint mill to meet the requirements of layered milling is also an urgent problem to be solved. Consequently, it is of great practical value to explore the changes of milled rice during the whole milling process, and to establish an objective and highly efficient breakpoint planning method based on this. Recently, with scanning electron microscopy (SEM), many researchers observed the relationship between the bran removal process and DOM in the micro view and gave the range of the DOM corresponding to each fraction of BR [13,14]. No doubt this provides a new technical scheme for the understanding of the bran removal process and control of the rice multibreak milling process. Unfortunately, change in microscopic structure is tough to monitor in rice milling. Physical properties, as the macroscopic reflection of the microscopic structure, can facilitate the assessment of different stages of the bran removal process. For example, Mohapatra and Bal [15] found that there is a transition point of the variation rate on the nonlinear curve of the relationship between color and the DOM. The corresponding DOM at the transition point is in good agreement with the observed result of SEM. However, the physical properties affected by milling are not only numerous but also complex. Extensive research suggests that the physical properties of milled cereal, such as size, shape, density, weight, and friction coefficient, are easy to change significantly in the milling process [16][17][18]. Moreover, in the milling process, changes of the properties are often highly nonlinear, and large volumes of the properties data need to be acquired and processed [19][20][21]. Hence, a scientific method that can effectively analyze the variety law of the milled particulate properties during the milling process is needed. Recently, with scanning electron microscopy (SEM), many researchers observed the relationship between the bran removal process and DOM in the micro view and gave the range of the DOM corresponding to each fraction of BR [13,14]. No doubt this provides a new technical scheme for the understanding of the bran removal process and control of the rice multibreak milling process. Unfortunately, change in microscopic structure is tough to monitor in rice milling. Physical properties, as the macroscopic reflection of the microscopic structure, can facilitate the assessment of different stages of the bran removal process. For example, Mohapatra and Bal [15] found that there is a transition point of the variation rate on the nonlinear curve of the relationship between color and the DOM. The corresponding DOM at the transition point is in good agreement with the observed result of SEM. However, the physical properties affected by milling are not only numerous but also complex. Extensive research suggests that the physical properties of milled cereal, such as size, shape, density, weight, and friction coefficient, are easy to change significantly in the milling process [16][17][18]. Moreover, in the milling process, changes of the properties are often highly nonlinear, and large volumes of the properties data need to be acquired and processed [19][20][21]. Hence, a scientific method that can effectively analyze the variety law of the milled particulate properties during the milling process is needed.
Pattern recognition techniques such as principal component analysis (PCA), hierarchical cluster analysis (HCA) and artificial neural networks (ANNs), which allow for visualization of the underlying structure between product properties and products, have been used widely in different fields of engineering [22][23][24]. For example, Barua et al. [25] used PCA and HCA to determine the crucial parameter and ideal processing for altering the rapidly digestible starch fraction from a multitude of processing parameters and conditions. Barros et al. [26] indicated that multivariate analysis allowed the classification of different edible seeds based on nonlinear physical and chemical characteristic parameters such as moisture, total phenolic compounds, and antioxidant activity. Yang et al. [27] achieved a rapid identification of beer species by multivariate analysis of forty volatile compounds. Sharin et al. [28] used multivariate analysis methods including PCA, HCA, and ANN to discriminate between different entomological-origin honey samples based on physicochemical properties and volatile compound characteristics. In overview, a large amount of nonlinear data processed by multivariate analysis can characterize and predict the similarity between samples [29][30][31]. Therefore, it is expected that a reasonable breakpoint planning method can be constructed by using the multivariate analysis method described due to the similarity of many physical properties with the variation of the DOM.
In this study, the method of multivariate analysis was used to comprehensively analyze the characteristic changes of milled rice during the milling process. On this basis, a first attempt of constructing a breakpoint planning method for the multibreak milling system was presented, in order to provide reference for the application of a commercial multibreak milling system.

Experimental Materials
The experimental japonica paddy (Dongnong 429) was supplied by the experiment station of Northeast Agricultural University (Harbin, China). The paddy was hulled by a sheller (THU-35B type, Satake, Tokyo, Japan, dehulling speed: 100-120 kg/h) and BR was obtained after removing unsound kernel, broken kernel, and foreign matter. The initial moisture content of BR was 12.1 ± 0.1% (w.b.), measured using the 105 • C dry oven (DGH-9053A type, Yiheng, Shanghai, China) gravimetric method [32]. Rice samples with similar size, shape, and no cracks were picked. The treated samples were stored in double-sealed polyethylene bags at 5 • C in a refrigerator before milling [15].

Friction-Type Milling Process
The milling process used a vertical automatic rice mill system (SY95-PC+PAE5 type, Ssangyong Machinery Co., Korea), as shown in Figure 2. Milling is conducted under friction and the softer rice layers are removed by pressure and friction from the milling chamber or each other. The schematic of the vertical rice mill is shown in Figure 3. The down-drift weight of specimens can be displayed by a real-time monitor in the mill. Desired degrees of milling (DOMs) were obtained and calculated with the formula below [4]: where M a is weight of milled rice, g; M b is weight of initial BR, g. From the above, ten 200 g sub-samples were taken from polyethylene bags in the refrigerator and milled, varying between 1% and 10% DOM. Milling duration is recorded in Table 1. The milled samples (1~10%) and BR (0% DOM) are both shown in Figure 4.

Determination of Physical Properties
Physical properties of 0~10% DOM sub-samples were determined accurately according to the following methods, including length (L), width (W), thickness (T), aspect ratio (R a ), volume (V), surface area (S), equivalent diameter (D p ), sphericity (φ), bulk density (ρ b ), true density (ρ t ), porosity (ε), thousand kernel weight (M q ), whiteness (L w ), chromaticity (B), and coefficient of static friction (µ).        Notes: All data are presented as the mean with the standard deviation in parentheses. ** Is very significant (p < 0.01); NS Is statistically not significant (p > 0.05). The mean values in the same row followed by the same superscript letter denote that they are not significantly different (p > 0.05); in contrast, the different letter denotes that they are significantly different (p < 0.05).

Size and Shape
The axial dimensions (Figure 5a), viz., length (L), width (W), and thickness (T), were measured using a digital vernier caliper (111N-102V-10 type, Guanglu, Nanjing, China) with an accuracy of 0.01 mm as principle dimensions. The 100 tested kernels were randomly selected from each DOM respectively. Other parameters of shape and size were then defined from these three linear dimensions. The axial dimensions (Figure 5a), viz., length (L), width (W), and thickness (T), were measured using a digital vernier caliper (111N-102V-10 type, Guanglu, Nanjing, China) with an accuracy of 0.01 mm as principle dimensions. The 100 tested kernels were randomly selected from each DOM respectively. Other parameters of shape and size were then defined from these three linear dimensions. The aspect ratio (Ra) was determined from Equation (2) [33]: The volume (V) and surface area (S) of rice kernels were derived the following expressions [34], respectively: The equivalent diameter (Dp) is expressed as [34]: The sphericity ( φ ) was calculated through the following expression [34]:

Density and Weight
According to the different measuring methods of volume, the bulk density (ρb) was The aspect ratio (R a ) was determined from Equation (2) [33]: The volume (V) and surface area (S) of rice kernels were derived the following expressions [34], respectively: The equivalent diameter (D p ) is expressed as [34]: The sphericity (φ) was calculated through the following expression [34]:

Density and Weight
According to the different measuring methods of volume, the bulk density (ρ b ) was the ratio of the mass (W b ) of the grain to its total volume (V b ), whereas the true density (ρ t ) was the ratio of the mass (W t ) of the grain to the solid volume (V t ) occupied by the sample. The bulk density was determined by filling an empty 10 mL graduated cylinder with the rice samples, consolidated 3 times to achieve uniformity in the filling process as reported by Irtwange et al. [35], and the weight of the samples was obtained by subtracting the weight of the cylinder from the weight of the cylinder and samples.
The true density (ρ t ) was measured by the toluene displacement method [34]. The true volume of the samples was measured by placing a known weight of samples in a measuring cylinder with toluene. The displaced volume was taken as the true volume of the sample. The process was replicated 5 times, and the bulk density and true density for each replication were calculated using the following expressions, respectively: The porosity (ε) was calculated from bulk and true densities by using the following relationship [34]: Thousand kernel weight (M q ) was determined by the weight of 1000 randomly selected sound kernels by means of electronic balance (accuracy of 0.0001 g) and replicated five times for each DOM.

Color and Lightness
The color of samples with different DOMs was measured using an automatic color meter (DC-P3 type, Xingguang, Beijing, China). The color meter adopts the CIE Lab L, a, b color space technique, where L indicates whiteness (L w ) and a and b are the chromaticity coordinates. The whiteness (L w ) can be measured and displayed directly on the automatic color meter. However, chromaticity (B) needs to be further calculated from the chromaticity coordinates (a and b) displayed on the automatic color meter. The measurement was replicated five times and the chromaticity (B) was calculated by the relation [36]:

The coefficient of Static Friction
For measuring the static friction coefficient between the rice kernel and milling chamber wall, a device was designed through modification by Heidarbeigi et al. [33]. The apparatus consisted of several mechanisms, namely, a screw nut, protractor, and lifting table having the same material as that of the milling chamber wall (Figure 5b). The lifting table was raised gradually until the filled sample just started to slide down and the angle (β) between the lifting table and the horizontal surface was recorded. The process was replicated five times. The coefficient of static friction (µ) was calculated by the following formula:

Statistical Analysis
We performed the multivariate analyses with the statistical package SPSS for Windows (version 23). All experimental data were expressed as mean ± standard deviations (SD) and subjected to a one-way analysis of variance (ANOVA). Significant difference was determined at the 0.05 probability level and post hoc tests were run by Tukey's Test at a 5% probability level. HCA and PCA were also applied to the data set after standardization by the following formula: where Z i -physical properties after standardization, X i -the measured physical properties, X-the mean value of the measured physical properties, σ sd -the standard deviation of the measured physical properties.

Generalized Regression Neural Network (GRNN)
GRNN was developed as an alternative to the radial basis function neural network. It can be used for any regression problem in which an assumption of linearity is not justified [37]. In this study, nonlinear supervised models were investigated by using the GRNN for estimating different stages of rice milling. In general, the smoothness of the approximated function needs be controlled by the smooth factor, namely spread (σ), in the training phase of GRNN. Consequently, the optimal spread value in the GRNN model was selected by circularly measuring the Mean Square Error (MSE), as shown in Equation (13), between the predicted values and the experimental values to obtain the minimum in the training set using cross validation (cross validation 4 times). The performance of the GRNN models was evaluated at the correct rate.

Change in Physical Properties during Milling
As a result of the one-way analysis of variance from Table 1, the true density did not present significant differences, while the rest of the physical properties differed significantly during the milling. Similar investigations have been carried out to understand how physical properties of various rice cultivars influencing the DOM and similar results were found by Liu et al. [38]. However, the physical properties of the same cultivar of rice with different DOMs were not investigated in other studies.
In the case of size and shape, all the axial dimensions decreased significantly with an increasing DOM; L, W, and T declined 6.6%, 7.6%, and 14.1%, respectively, when the DOM reached 10%. Compared with L and W, T had more significant changes, suggesting that with the complex distribution of friction and rice attitude in the milling chamber, the brown rice was subjected to the most intense milling in the thickness direction [13]. Analogously, the equivalent diameter, volume, and surface area also showed a trend of decline. However, the aspect ratio and sphericity varied irregularly. Moreover, it must be noted that the standard deviations (Table 1) of all the dimensions at the beginning of the grinding were larger than that at 10% DOM, which might indicate that the BR was milled unevenly at the initial milling stage [39,40].
The thousand kernel weight decreased with the increase in DOM and varied from 21.01 to 18.25 g. The average weight of each grain was reduced by 13.8%, and the difference in weight was important for designing machine components such as those for ventilation and winnowers [34]. The kernels' bulk density of varied from 761.27 to 824.45 kg/m 3 , and it rose with the increase in DOM; contrarily, porosity declined from 42.90 to 38.17% with the increase in DOM. Moreover, the true density practically did not present differences among the varieties in diverse DOMs. It can be inferred that the true density was not affected by milling, because the true density of rice might depend on the compact core endosperm instead of the removed bran layer. The variation tendency in these physical properties is consistent with those reported by Corrêa et al. [41]. In the case of the static coefficient of friction, the DOM increased gradually which led to the decrease in the static friction coefficient, due to the bran layer being clearly diminished. That may also explain why L w Foods 2023, 12, 1864 9 of 17 increased and B decreased with the increase in the DOM. The levels of pigments in the bran layer are significantly higher than the middle endosperm [7].

The Breakpoint Planning Method
According to Tukey's Test in Table 1, we can infer that each of the measured physical properties, except true density, aspect ratio, and sphericity, were regularly variable and significantly different (p < 0.05) between different DOMs. Therefore, we try to put forward a breakpoint planning method according to the changes of physical properties of milled rice during the whole rice milling process. The breakpoint planning method presented in this paper is divided into two stages, as illustrated in Figure 6. In the first stage, a continuous single-break rice milling experiment is required to obtain the characteristic changes of the milled rice during the whole milling process and to determine the grading strategy. Basically, this stage will produce the number of breakpoints (n) required the DOM or milling time (MT) corresponding to each breakpoint, and the representative characteristics of milled rice. In the second stage, using the grading strategy obtained in the first stage, a multibreak milling system with n breakpoints is constructed, and the actual breakpoint to which the milled rice obtained from each breakpoint belongs is estimated. Then, the milled rice obtained from each breakpoint is redistributed. Finally, milled rice products that meet different needs are obtained. The following sections describe each stage of the method in greater detail.
why Lw increased and B decreased with the increase in the DOM. The levels of pigments in the bran layer are significantly higher than the middle endosperm [7].

The Breakpoint Planning Method
According to Tukey's Test in Table 1, we can infer that each of the measured physical properties, except true density, aspect ratio, and sphericity, were regularly variable and significantly different (p < 0.05) between different DOMs. Therefore, we try to put forward a breakpoint planning method according to the changes of physical properties of milled rice during the whole rice milling process. The breakpoint planning method presented in this paper is divided into two stages, as illustrated in Figure 6. In the first stage, a continuous single-break rice milling experiment is required to obtain the characteristic changes of the milled rice during the whole milling process and to determine the grading strategy. Basically, this stage will produce the number of breakpoints (n) required the DOM or milling time (MT) corresponding to each breakpoint, and the representative characteristics of milled rice. In the second stage, using the grading strategy obtained in the first stage, a multibreak milling system with n breakpoints is constructed, and the actual breakpoint to which the milled rice obtained from each breakpoint belongs is estimated. Then, the milled rice obtained from each breakpoint is redistributed. Finally, milled rice products that meet different needs are obtained. The following sections describe each stage of the method in greater detail.

Grading Strategy
At present, the number of breakpoints of rice multibreak milling systems in rice mills is usually set by the experience of operators. In addition, the number of breakpoints in the rice multibreak milling system also needs to be adjusted according to the varieties of milled rice. However, to date, no objective method exists for planning the number of breakpoints in a reasonable way. Therefore, in order to propose a reasonable breakpoint planning method, it is necessary to grade the whole rice milling process to determine the number of breakpoints in the rice multibreak milling system. As shown in Figure 6, in the

Grading Strategy
At present, the number of breakpoints of rice multibreak milling systems in rice mills is usually set by the experience of operators. In addition, the number of breakpoints in the rice multibreak milling system also needs to be adjusted according to the varieties of milled rice. However, to date, no objective method exists for planning the number of breakpoints in a reasonable way. Therefore, in order to propose a reasonable breakpoint planning method, it is necessary to grade the whole rice milling process to determine the number of breakpoints in the rice multibreak milling system. As shown in Figure 6, in the first stage, the characteristic changes of milled rice with DOM/MT are first analyzed. It is worth mentioning that, in this study, we mainly recorded and analyzed the characteristic changes of milled rice with DOM rather than MT. However, for the breakpoint planning method proposed in this study, it is also effective to record and analyze the characteristic changes of milled rice with MT. In addition, a further explanation of the changes in the characteristics of milled rice during the whole milling process has been given in Section 3. Therefore, in order to avoid subjective evaluation of the grade of milled rice only according to the characteristic changes of milled rice, we then introduce an objective grading method, that is, HCA.

Hierarchical Cluster Analysis
HCA is widely used to cluster objects into unknown groups and define new distinctive clusters [42]. In general, samples are grouped on the basis of similarities without taking into account the information about class membership, and the difference between observed individuals uses a defined metric (e.g., Euclidean, squared Euclidean, and Minkowski distance et al.) and clusters by Ward's, median, and centroid clustering methods et al. In this work, Ward's method and the squared Euclidean distance method were used as clustering methods. In order to avoid the effect of milling non-uniformity on the HCA, the data set for the HCA was made up of mean values of physical properties except for true density, aspect ratio, and sphericity. Then, standardization was conducted on the data set using Equation (12). The resultant information is presented in the dendrogram (Figure 7). changes of milled rice with DOM rather than MT. However, for the breakpoint planning method proposed in this study, it is also effective to record and analyze the characteristic changes of milled rice with MT. In addition, a further explanation of the changes in the characteristics of milled rice during the whole milling process has been given in Section 3. Therefore, in order to avoid subjective evaluation of the grade of milled rice only according to the characteristic changes of milled rice, we then introduce an objective grading method, that is, HCA.

Hierarchical Cluster Analysis
HCA is widely used to cluster objects into unknown groups and define new distinctive clusters [42]. In general, samples are grouped on the basis of similarities without taking into account the information about class membership, and the difference between observed individuals uses a defined metric (e.g., Euclidean, squared Euclidean, and Minkowski distance et al.) and clusters by Ward's, median, and centroid clustering methods et al. In this work, Ward's method and the squared Euclidean distance method were used as clustering methods. In order to avoid the effect of milling non-uniformity on the HCA, the data set for the HCA was made up of mean values of physical properties except for true density, aspect ratio, and sphericity. Then, standardization was conducted on the data set using Equation (12). The resultant information is presented in the dendrogram (Figure 7). The BR was clearly separated from milled BR; therefore, it was a specific category. Then, 1~10% DOM was redefined as four grades: Ⅰ, Ⅱ, Ⅲ, and Ⅳ. Group Ⅰ contained a 9~10% DOM; group Ⅱ contained a 6~8% DOM; group Ⅲ contained a 4~5% DOM; and group Ⅳ contained a 1~3% DOM. These clusters might suggest that structure and composition of milled rice importantly affected physical properties during milling. Moreover, in respect of milling duration (Table 1), the average time was within 10 s and 20 s for DOMs of grades Ⅳ and Ⅲ, respectively (removal fractions were designated as the bran layer and outer endosperm layer). However, for the DOM of grade Ⅱ, milling duration lasted about 1 min and the milling speed was significantly reduced. It can be explained The BR was clearly separated from milled BR; therefore, it was a specific category. Then, 1~10% DOM was redefined as four grades: I, II, III, and IV. Group I contained a 9~10% DOM; group II contained a 6~8% DOM; group III contained a 4~5% DOM; and group IV contained a 1~3% DOM. These clusters might suggest that structure and composition of milled rice importantly affected physical properties during milling. Moreover, in respect of milling duration (Table 1), the average time was within 10 s and 20 s for DOMs of grades IV and III, respectively (removal fractions were designated as the bran layer and outer endosperm layer). However, for the DOM of grade II, milling duration lasted about 1 min and the milling speed was significantly reduced. It can be explained that compared to the bran layer and outer endosperm layer, the middle endosperm layer provides greater strength to hinder milling [15]. While it needs to take longer to achieve the DOM of grade I, the physical properties were almost unchanged (DOMs of 9% and 10% were clustered into one group firstly, due to the fact that they are the most similar, as seen in Figure 7).

Setting up a Multibreak Milling System with n Breakpoints
In the first stage, the grading strategy can help operators classify the rice milling process objectively. According to the grading results of the first stage, a rice multibreak milling system with n (n is equal to the grading results obtained in the first stage; in this paper, n = 4) breakpoints can be constructed. However, the products of each single break are inevitably subjected to uneven milling [3], which undoubtedly reduces the effectiveness of the multibreak milling system. Therefore, on the basis of the grading strategy, the actual breakpoint which the milled rice obtains from each breakpoint needs be accurately assessed. Then, the milled rice obtained from each breakpoint needs be further redistributed. Further explanation of the redistribution strategy of milled rice is given below.

Redistribution of Milled Rice Properties Reduction
In order to further assess the actual breakpoint to which the milled rice obtained from each breakpoint belongs, we need to establish the relationship between the above gradually changing physical properties during the milling process and the four different stages obtained by HCA. However, at present, there are 12 properties that have been extracted. Those excessive properties increase computation time and storage memory, which may cause the assessment process to become more complicated and even decrease the performance of the assessment model. Therefore, a strategy is necessary to reduce the number of properties used in assessment. In this study, PCA was performed with the aim to better evaluate which physical properties were more significant in the characterization of the DOM. PCA was applied to measure 12 physical properties except the true density, aspect ratio, and sphericity. The three measured values of each physical property under different DOMs were chosen to form a 33 × 12 data set. Then, standardization was conducted on the data set using Equation (12). After standardization, all indices contribute equally to data set variance and carry equal weight in principal component calculation [29].
The scree plot (Figure 8a) shows eigenvalues of all principal components, in which the first two PCs with eigenvalues higher than 1 were extracted. Among the extracted PCs, PC1 was responsible for 71.44% of the total variance, while PC2 accounted for 19.33%. The first two principal components (PC1 and PC2) accounted for 90.77% of the total variance, which meets the general requirement of 85%, indicating that the two PCs, covering most of the information contained in 12 physical properties, might represent the changes of rice structure and composition characteristics from different aspects during milling. To further elucidate meaning of new comprehensive factors and discriminate between PC1 and PC2, the rotated solution of extracted PCs was calculated using the Varimax. The rotated loading plot (Figure 8b) shows that only the size and shape properties were co-located on the upper part of the rotated loading plot, indicating that they are the major constituents of the PC2. In addition, the PC1 consists mainly of other physical properties including L w and ρ b exhibiting high positive loadings, while ε, µ, M q , and B have the higher negative loadings. According to the results of the rotated loading, PC1 and PC2 can be clearly interpreted as an indirect evaluation factor (IEF) and direct evaluation factor (DEF) of the DOM, respectively. Additionally, PC1 and PC2 can be quantified through the following expression: PC1 = a 1 L + a 2 W + a 3 T + a 4 V + a 5 S + a 6 D p + a 7 ρ + a 8 ε + a 9 µ + a 10 M q + a 11 L w + a 12 B where a 1~a12 and b 1~b12 are the eigenvectors of PC1 and PC2, respectively. Corresponding eigenvectors are shown in Table 2.  a  D  a  S  a  V  a  T  a  W  a  L  a  PC   w  q  p  12  11  10  9  8  7  6  5  4  3  2  1 1 w  q  p  12  11  10  9  8  7  6  5  4  3  2  1 2 where 1 a ~1 2 a and 1 b ~1 2 b are the eigenvectors of PC1 and PC2, respectively. Corresponding eigenvectors are shown in Table 2.    On the other hand, Table 2 also gives the rotated loadings, and, as expected, PC1 clearly discriminated B, L w , ρ b , ε, and M q from the other lower loadings of physical properties, while PC2 clearly discriminated V, S, and D p with higher loadings from the other physical properties. These results might suggest that a few representative physical properties can be selected from two comprehensive factors. Considering that PC1 represents the maximum variation of the data set, firstly, B, L w , ρ b , ε, and M q were selected as representative physical properties. However, two comprehensive factors (IEF and DEF) of the DOM might represent different changes in the DOM during rice milling. Therefore, V, S, and D p (from PC2) and B and L w (from PC1) were chosen as another group of representative physical properties. To validate whether the most representative variables were selected from the two groups of physical properties (B, L w , ρ b , ε, and M q and V, S, D p , B, and L w ), the above HCA was conducted. The results obtained are shown as two dendrograms, respectively (Figure 9). The result (as shown in Figure 9b) of the HCA for B, L w , ρ b , ε, and M q was clearly distinct from the result of the original HCA (Figure 7), while for V, S, D p , B, and L w , the result (as shown in Figure 9a) of the HCA was essentially in agreement with that of the original HCA (Figure 7), except for clustering order (6%, 7%, and 8% DOMs are first clustered in a group in Figure 9a, while 9% and 10% DOMs are first clustered in a group in Figure 7). Therefore, these five physical properties (V, S, D p , B, L w ) contained enough information to replace the original twelve physical properties in the evaluation of the DOM. erties, while PC2 clearly discriminated V, S, and Dp with higher loadings from the other physical properties. These results might suggest that a few representative physical properties can be selected from two comprehensive factors. Considering that PC1 represents the maximum variation of the data set, firstly, B, Lw, ρb, ε, and Mq were selected as representative physical properties. However, two comprehensive factors (IEF and DEF) of the DOM might represent different changes in the DOM during rice milling. Therefore, V, S, and Dp (from PC2) and B and Lw (from PC1) were chosen as another group of representative physical properties. To validate whether the most representative variables were selected from the two groups of physical properties (B, Lw, ρb, ε, and Mq and V, S, Dp, B, and Lw), the above HCA was conducted. The results obtained are shown as two dendrograms, respectively ( Figure 9). The result (as shown in Figure 9b) of the HCA for B, Lw, ρb, ε, and Mq was clearly distinct from the result of the original HCA (Figure 7), while for V, S, Dp, B, and Lw, the result (as shown in Figure 9a) of the HCA was essentially in agreement with that of the original HCA (Figure 7), except for clustering order (6%, 7%, and 8% DOMs are first clustered in a group in Figure 9a, while 9% and 10% DOMs are first clustered in a group in Figure 7). Therefore, these five physical properties (V, S, Dp, B, Lw) contained enough information to replace the original twelve physical properties in the evaluation of the DOM.

Assessment Model
The most important step in the redistribution strategy is to construct an assessment model to evaluate the actual breakpoint to which the milled rice belongs. However, the relationship between physical properties and DOM was nonlinear. In recent years, ANNs have been used in different fields of engineering because of their capability of extracting complex and nonlinear relationships [43]. In this study, the relationship model of physical properties and four grades (four breakpoints) was constructed by GRNN, which is one type of Radial Basis Function (RBF).
GRNN was applied to the data collection, which was consistent with that of PCA, and two-thirds of the data were used as training data and one-third of the data was used

Assessment Model
The most important step in the redistribution strategy is to construct an assessment model to evaluate the actual breakpoint to which the milled rice belongs. However, the relationship between physical properties and DOM was nonlinear. In recent years, ANNs have been used in different fields of engineering because of their capability of extracting complex and nonlinear relationships [43]. In this study, the relationship model of physical properties and four grades (four breakpoints) was constructed by GRNN, which is one type of Radial Basis Function (RBF).
GRNN was applied to the data collection, which was consistent with that of PCA, and two-thirds of the data were used as training data and one-third of the data was used as test data. Both the inputs and output data were normalized to the range of [-1, 1] using Equation (16) [44].
where p n is the normalized parameter, p denotes the actual parameter, p min represents a minimum of the actual parameters, and p max stands for a maximum of the actual parameters.
According to the results of HCA and PCA, two GRNN architectures were constructed: type A has two input notes, which are PC1 and PC2 (PC1 and PC2 were calculated by Equations (14) and (15), respectively) (Figure 10a), while type B has five input notes that are V, S, D p , B, and L w (Figure 10b). The spread of type A was set from 0.1 to 1 at 0.1 intervals, and that of type B was set from 0.05 to 0.5 at 0.05 intervals. The changed MSE with the spread values of type A and type B are shown in Figure 11. The two curves both have a well-defined minima, with the optimal value of spread at 0.4 and 0.2. In this case, the optimal input and output of cross validation were also saved for constructing the optimal GRNN. min max where pn is the normalized parameter, p denotes the actual parameter, pmin represents a minimum of the actual parameters, and pmax stands for a maximum of the actual parameters.
According to the results of HCA and PCA, two GRNN architectures were constructed: type A has two input notes, which are PC1 and PC2 (PC1 and PC2 were calculated by Equations (14) and (15), respectively) (Figure 10a), while type B has five input notes that are V, S, Dp, B, and Lw (Figure 10b). The spread of type A was set from 0.1 to 1 at 0.1 intervals, and that of type B was set from 0.05 to 0.5 at 0.05 intervals. The changed MSE with the spread values of type A and type B are shown in Figure 11. The two curves both have a well-defined minima, with the optimal value of spread at 0.4 and 0.2. In this case, the optimal input and output of cross validation were also saved for constructing the optimal GRNN. Four grades in the test data set (including 11 data points) were predicted using the optimal GRNN which was developed by the above training. The GRNN of type A showed the correct rate of 100%, while the GRNN of type B showed 90.91%. The predictive performance of type A was better than that of type B, suggesting that PC1 and PC2 do contain more information of physical properties to accurately assess results. However, where pn is the normalized parameter, p denotes the actual parameter, pmin represents a minimum of the actual parameters, and pmax stands for a maximum of the actual parameters.
According to the results of HCA and PCA, two GRNN architectures were constructed: type A has two input notes, which are PC1 and PC2 (PC1 and PC2 were calculated by Equations (14) and (15), respectively) (Figure 10a), while type B has five input notes that are V, S, Dp, B, and Lw (Figure 10b). The spread of type A was set from 0.1 to 1 at 0.1 intervals, and that of type B was set from 0.05 to 0.5 at 0.05 intervals. The changed MSE with the spread values of type A and type B are shown in Figure 11. The two curves both have a well-defined minima, with the optimal value of spread at 0.4 and 0.2. In this case, the optimal input and output of cross validation were also saved for constructing the optimal GRNN. Four grades in the test data set (including 11 data points) were predicted using the optimal GRNN which was developed by the above training. The GRNN of type A showed the correct rate of 100%, while the GRNN of type B showed 90.91%. The predictive performance of type A was better than that of type B, suggesting that PC1 and PC2 do contain more information of physical properties to accurately assess results. However, Four grades in the test data set (including 11 data points) were predicted using the optimal GRNN which was developed by the above training. The GRNN of type A showed the correct rate of 100%, while the GRNN of type B showed 90.91%. The predictive performance of type A was better than that of type B, suggesting that PC1 and PC2 do contain more information of physical properties to accurately assess results. However, considering the real-time evaluation of rice milling, type B is more likely to achieve success, due to the five input nodes (V, S, D p , B, L w ) being easier to measure.

Effectiveness of Planning Method
To further verify the availability of the whole planning method including the grading strategy and redistribution strategy, the iso-color-phase differential staining (IDS) method was used as a standard method for determination of the processing degree of rice (Chinese standard GB 18105-2000, 2000) [45] and was compared with the estimated results obtained by the method proposed in this paper. First, the grading results proposed in this paper are the same as the standard grading method which divides the milled rice into four grades (Chinese standard GB 1354GB -2009GB , 2009 [46]. Then, to further verify the stability of the assessment model (GRNN model), the BR samples were milled to four different DOMs (this is an attempt to simulate the multibreak milling process). Next, we assessed the actual grades to which the milled rice obtained from each DOM belongs (the number of samples to be tested for each DOM is 40) and compared the results of the assessment with the grades of standard samples (note that we only use GRNN of type B). The results of the comparison are shown in Table 3. From our statistical results, if the multibreak planning method was not used, a large amount of rice would be overmilled or fail to reach the corresponding milling grades. In contrast, the planning method in this article can not only provide reasonable breakpoints but also more accurately redistribute the milled rice (the average correct rate can reach 94.28%, and grade II and grade III are more likely to be confused).

Conclusions
In this study, the variation of physical parameters during the rice milling process was analyzed using multivariate analysis. A breakpoint planning method was developed for the multibreak milling system based on the similarity of physical properties with the variation of DOM. Furthermore, a GRNN model was constructed to evaluate the actual breakpoints to which the milled rice obtained from each breakpoint belongs. The validation results show that the method has good accuracy. The main conclusions are as follows.
(1) The grading strategy based on the characteristic changes of the milled rice during the whole milling process enables the rice milling process to be divided into four different stages, group I contained a 9~10% DOM, group II contained a 6~8% DOM, group III contained a 4~5% DOM, and group IV contained a 1~3% DOM. The results of this stage, in a multibreak milling system, will help to initially determine the number of breakpoints required and DOM/MT corresponding to each breakpoint.
(2) The boundaries of each breakpoint are defined, which helps the operator to set the initial end point of each single breakpoint rice mill. For example, under the conditions of this study, the number of breakpoints of the multibreak milling system should be set to four and DOMs of 3%, 5%, 8%, and 10% can be used as the boundary of each breakpoint.
(3) The GRNN model is constructed to assess the actual breakpoints to which the milled rice obtained from each breakpoint belongs, because the products of each single break are inevitably subjected to uneven milling. The verification results show good accuracy (94.28%).
This work will help operators to plan the multibreak milling system of rice efficiently and objectively so as to significantly improve the commercial value of milled rice. Moreover, the method proposed in this study can also provide a reference for the breakpoint planning of the multibreak milling systems of other heterogeneous granular materials.