Investigation of Error Distribution in the Back-Calculation of Breakage Function Model Parameters via Nonlinear Programming

: Despite its effectiveness in determining breakage function parameters (BFPs) for quan-tifying breakage characteristics in mineral grinding processes, the back-calculation method has limitations owing to the uncertainty regarding the distribution of the error function. In this work, using Korean uranium and molybdenum ores, we show that the limitation can be overcome by searching over a wide range of initial values based on the conjugate gradient method. We also visualized the distribution of the sum of squares of the error in the two-dimensional parameter space. The results showed that the error function was strictly convex, and the main problem in the back-calculation of the breakage functions was the ﬂat surface of the objective function rather than the occurrence of local minima. Based on our results, we inferred that the ﬂat surface problem could be signiﬁcantly mitigated by searching over a wide range of initial values. Back-calculation using a wide range of initial values yields BFPs similar to those obtained from single-sized-feed breakage tests (SSFBTs) up to four-dimensional parameter spaces. Therefore, by searching over a wide range of initial values, the feasibility of the back-calculation approach can be signiﬁcantly improved with a minimum number of SSFBTs.


Introduction
Grinding, one of the most important unit operations in the field of mineral processing, is a very energy-inefficient process. More than 50% of the total energy required for all the mineral processing operations is used in the grinding stage, and only 10-20% of the input energy is used for actual size reduction [1]. Therefore, owing to its inefficiency, grinding has been an important research subject in the field of mineral processing.
Mathematical breakage modeling has been used as an efficient tool to characterize the grinding process. In a breakage model, the material-and device-specific breakage properties are characterized based on several parameters using the breakage function [2]. Classical breakage models only provide empirical size-energy relationships based on the proportionality between the energy and surface area [3], size reduction ratio [4], or length of new cracks [5]. In contrast, kinetic grinding (KG) models based on material balance enable predicting the product size distribution from the given feed-size distribution and time [6][7][8][9]. In recent years, mechanistic models such as the discrete element method have been extensively used for microscopic analysis of the grinding process, equipment design, and optimization of the operating conditions [10][11][12][13][14][15][16]. In this work, we focused on optimization issues using the population balance model. There are two representative grinding models that have been widely used in the field of mineral processing. KG models [6,7,9] are based on mathematical modeling of the population balance using grinding kinetics.
The Julius Kruttschnitt Mineral Research Centre (JKMRC) grinding model [17][18][19][20][21][22][23] uses the relationship between the energy and product size distribution, usually obtained via single-particle breakage tests using a drop weight or pendulum. The KG and JKMRC models can be transformed into each other under special conditions. This work is based on the KG model developed by Austin, Luckie [6,7], and Reid [9]. The basic concept of the models has been used in numerous studies to characterize the grinding properties of various materials with various operating variables in ball mills, in terms of the breakage rate and the breakage distribution [24][25][26][27][28][29].
In a KG model, the breakage process is considered a rate process and decoupled into two basic functions: Breakage rate function and breakage distribution function. These two functions mainly depend on the material properties, mill dimensions, and operating conditions. They are generally determined via a set of laboratory-scale breakage tests with single-sized feeds. For determining the breakage rate function, feed particles are separated into √ 2 sieve fractions. To verify the breakage kinetics, a single size fraction is ground for a specified amount of time. The fraction of unground particles is then measured. Grinding and sieving are repeated with varying grinding times to obtain time-series data, and then the breakage rate for the given size fraction is estimated from the time-series data. Then, the process must be repeated for different size fractions to obtain a comprehensive set of breakage rates as a function of size.
Determining the breakage function is an extremely tedious and long process as it requires measuring breakage rates for a wide range of particles with single size fractions. The breakage functions obtained from laboratory-scale experiments need to be corrected for use in an industrial grinding circuit. To establish correlations between lab-scale and industrial mills, empirical equations are used for determining the dependence of the breakage parameters on ball/mill dimensions and operating variables. When scaling up the model parameters from lab-scale to industrial mills, additional errors can occur [2].
To minimize the high costs associated with experimental work, back-calculation of BFPs from the size distributions with various grinding times is proposed [30]. The approach has advantages such as minimized experimental work and applicability for cases where direct determination of the breakage rate and breakage distribution function is not possible. Moreover, the method has been applied in designing various grinding processes [25,[31][32][33].
In principle, back-calculation of breakage parameters is an optimization process from a mathematical perspective, and several issues such as the locality of optimum have been addressed. If the error function defined by the square sum of the difference between the experimental and calculated values has multiple convex point sets depending on the initial guess of parameters, then back-calculation may fail by ending up in the local minimum. If the objective function surface is flat around the optimum, then convergence may fail before approaching the optimum. Owing to the uncertainty regarding the shape of the objective function during back-calculation, the use of this method has been limited to that of an auxiliary tool to fine-tune the values obtained from single-sized-feed breakage tests (SSFBTs). Instead, the information from grinding experiments is used to obtain the best results. In such cases, the degree to which the experimental work can be reduced is very limited, and the process of determining the breakage function again becomes extremely tedious.
Earlier, when the back-calculation approach was first proposed for particle comminution in mineral processing, the distribution of the objective function in the back-calculation of the parameters for determining the breakage function was not understood clearly owing to computational limitations. If an optimal parameter search is conducted over a wide range of possible parameter values (wide-range initial value search, WRIVS, hereinafter), the computational load on a single processor increases by the order of N d N p , where N p and N d are the number of parameters and number of discrete points for each parameter, respectively. For a breakage model, five to seven parameters are normally required [2]; hence, assuming 100 discrete points for each parameter increases the number of searching points up to 10 10 -10 14 . Owing to the recent advances in parallel processing, multithread computing using several threads that correspond to the numbers mentioned above is possible on a personal computer. Given that breakage parameter determination does not have to be carried out in a short period of time, the magnitude of the computational load required for global searching of the breakage function is in the acceptable range. The WRIVS strategy therefore is currently accessible with sufficiently high resolution without any major concern with regard to the computational cost.
In this work, the feasibility of the wide-range search approach for identifying the BFPs with minimum or no SSFBTs was assessed. To this end, a parallel processing code that optimizes the BFPs was developed based on the conjugate gradient method. The error function distribution during back-calculation of the breakage parameters was analyzed. The shape of the error function under the condition in which two parameters are variable while the other parameters are fixed was visualized and analyzed. For convenience, the real coordinate space with the coordinates of the BFPs will be referred to as the parameter space hereafter. For example, the error function distribution when two parameters vary while the others remain fixed will be referred to as the error function distribution in a two-dimensional parameter space.
Next, a WRIVS in the parameter space extending up to five dimensions was conducted, and the resultant parameters that minimized the error were compared with the parameters obtained from the SSFBTs.

Methodology
While the essential details of the equations for the KG approach are described in detail in the present paper, further information on these equations can be referred from Austin, Klimpel and Luckie [2].
In the kinetic approach for modeling the grinding process, the breakage processes are described using two basis functions, namely the specific rate of breakage and primary breakage distribution. Experimental batch grinding results typically demonstrate that the rate of breakage follows first-order kinetics, which can be written in the mathematical form as follows: where S i is the specific rate of breakage of a material in size class i, and w i (t) is the weight fraction at time t. The size dependence of the specific rate of breakage can be empirically expressed using the following functional form: where x i is the representative size of size class i, and x o is the standard size that is typically 1 mm. A, α, µ, and λ are model parameters. A represents the overall rate at which particles can be broken into smaller sizes, whereas α represents the size dependence of the specific rate of breakage. The second part on the right-hand side that includes µ and λ is the correction term that describes the decrease in the specific rate of breakage as the particle size exceeds the normal range to be broken efficiently by the balls. If the ball size is sufficiently large compared to the particle size, the second term becomes equal to one and Equation (2) can be reduced to The breakage of a particle produces an entire range of progeny particles. Particle breakage involves a combined mechanism including impact, compression, chipping, and abrasion. The cumulative distribution of the daughter fragments that can be produced from the combined mechanism can be conveniently described using a parametric functional form as follows: where B ij is the weight fraction that appears size under i from the primary breakage (with no refracture) of a particle with size j, and N is the total number of size classes. ϕ, γ, and β are model parameters. γ describes how progeny particles produced purely through impact/compression breakage are distributed. ϕ and β indicate the major breakage mechanisms. For example, if the breakage mechanism is fully dominated by impact/compression, then ϕ converges to one or β converges to γ. By combining functions S and B, we can describe the complete size-mass balance in a batch grinding system. The transition rate of the weight fraction of particles with size i is given as follows: The first term on the right-hand side describes the transition rate of the weight fraction to smaller sizes through grinding, whereas the second term with the sigma notation indicates the increase in size class i owing to the breakage of larger particles. The solution of Equation (5) based on the mass balance can be written as follows [9]: For a known feed-size distribution, the seven breakage parameters in Equations (2) and (4), i.e., X = (A, α, µ, λ, ϕ, γ, β) or five parameters in the reduced form, i.e., X = (A, α, ϕ, γ, β) fully determine the solution of Equation (6). These parameters depend on the material properties, equipment, and operating conditions. The breakage functions determined in the lab can be scaled up to the plant level by correlating the change in the breakage functions through empirical relationships [2], and they can be extended to describe continuous mill models by combining with additional sub-models such as the residence time distribution [34] and classifier models [35].
The breakage parameters are determined via a set of SSFBTs. The process requires repetitive grinding and classification to investigate the kinetics of breakage for a single-sized feed, and this process must be repeated for several size classes of the feed to investigate the size dependence of the breakage functions. Therefore, this process requires considerable experimental costs and time.
A nonlinear optimization process using grinding data from various grinding processes can be used to mitigate the experimental costs [30]. However, given that the relationship between the grinding parameters and product size distribution is highly nonlinear, the accuracy of back-calculation can be affected because of a poorly chosen initial guess, and there are possibilities of obtaining similar results from different combinations of breakage parameters. Many back-calculation processes involving a large number of variables have shown that the accuracy decreases as the number of variables to be optimized increases [25,[30][31][32]. While these researchers presumed that the accuracy decreased owing to the change in the shape of the objective function with the increase in the number of optimization variables, their presumption has not been validated by further studies. The shape of the objective function in the back-calculation process has been poorly understood owing to computational limitations.
Several nonlinear programming approaches are available for back-calculation. In this work, the conjugate gradient method [36] was used, given that Equation (4) is differentiable with respect to any of the breakage parameters. Compared to a derivative-free optimization method used in previous studies such as the simplex method, our approach can reduce the computational costs required for the solution to converge. E, the error function or optimization cost function, can be determined by the sum of squares (SSQ) of the difference between the cumulative weight fraction of the grinding product from the experimental result and model-predicted value. The optimization function is mathematically expressed as follows: where P i is the weight fraction of the grinding products smaller than size i, and r is the number of pairs used. In this work, data pertaining to three single-sized feeds and six grinding times were used, and therefore, the number of pairs used for optimizing the input is r = 6 C 2 × 3 = 45, where n C k denotes the number of k-combinations from a given set of n elements. When credible data for some of the parameters are available, the accuracy of the back-calculation method can be enhanced by reducing the number of parameters to be back-calculated. For example, using three known parameters, namely ϕ, γ, and β, backcalculation can be conducted only for four parameters, i.e., A, α, µ, and λ; moreover, in this case, the credibility and efficiency of back-calculation can be enhanced compared to that when searching seven parameters.

Experimental
The experimental conditions are listed in Table 1. Two ore samples were collected and used for the experiments. The representative molybdenum ore sample was collected from the Jecheon mine, and the uranium ore sample was collected from the Okcheon belt, Korea. Large lumps of ores were crushed in a jaw crusher and separated into √ 2 fractions using a Ro-Tap sieve shaker. Three size classes (1180 × 834 µm, 590 × 417 µm, and 295 × 209 µm) were selected for the grinding tests. A lab-scale ball mill with an internal dimeter of 200 mm and length of 165 mm was used for the grinding tests ( Figure 1). Inside the mill chamber, six half-cylinder-shaped lifters with a diameter of 25 mm are symmetrically located in the horizontal direction. Stainless-steel ball bearings with a diameter of 25.4 mm and density of 7840 kg/m 3 are used as the grinding media. The chamber was filled with balls corresponding to a chamber volume fraction of 0.2 and material corresponding to 50% of the volume of the interstitial void space in the ball bed. During grinding, the chamber was rotated on a roller table at a speed equal to 75% of the critical speed at which theoretical centrifugation takes place. After the feed material was ground for a suitably short amount of time, the mill was emptied, and 50 g of the powder was sampled using a riffle sampler. The size distribution of the sample was then analyzed by wet screening down to 36 µm using a series of √ 2 US standard sieves. The entire material including the sample was returned to the mill, and then operated for an additional duration. The procedure was repeated until S i in Equation (1) could be obtained from the linear-log plot of w i (t)/w i (0) versus time t.

Back-Calculation and Analysis
Data from the batch grinding tests were used for back-calculating breakage parameters. Back-calculation starts with the initial guessing of variables, and the final solution obtained from the optimization process can be significantly affected by the initial guess points. In this work, we varied the points during initial guessing of the parameter values over a wide range. The experiments and back-calculation process were conducted for normal size ranges where the specific breakage rate did not decrease. Equation (3) was used instead of Equation (2), and therefore, the number of parameters to be back-calculated was reduced to five.
For convenience, the error value corresponding to the given values of the input parameters and the final error at the point of convergence obtained from back-calculation are referred to as SSQ and SSQ , respectively. The schematic explaining the concepts of SSQ and SSQ is illustrated in Figure 2. First, we investigated the distribution of the error originating from the given parameter values. During this stage, for visualization, we varied two parameters and kept the other three parameter values fixed. The parameters obtained from SSFBTs were used as the fixed values. The local minima and surface shape of SSQ were investigated over the two-dimensional space of the variables.
Then, we investigated the distribution of SSQ . The tendency of SSQ in the twodimensional parameter space was compared with the distribution of SSQ to investigate the sensitivity and convergence characteristics of the back-calculation process. The back-calculated parameter values that produced minimal SSQ values over a wide range of the initial guess were then compared to the values obtained from SSFBTs.
Finally, we increased the dimensions of back-calculation to between three and five and conducted calculations over a wide range of initial guess values. Then, we compared the back-calculated parameters producing minimal SSQ values to those obtained from SSFBTs. All tests were conducted under both wet and dry conditions. In the wet grinding tests, water was added to obtain 62 wt.% slurry.

Back-Calculation and Analysis
Data from the batch grinding tests were used for back-calculating breakage parameters. Back-calculation starts with the initial guessing of variables, and the final solution obtained from the optimization process can be significantly affected by the initial guess points. In this work, we varied the points during initial guessing of the parameter values over a wide range. The experiments and back-calculation process were conducted for normal size ranges where the specific breakage rate did not decrease. Equation (3) was used instead of Equation (2), and therefore, the number of parameters to be back-calculated was reduced to five.
For convenience, the error value corresponding to the given values of the input parameters and the final error at the point of convergence obtained from back-calculation are referred to as SSQ fix and SSQ opt , respectively. The schematic explaining the concepts of SSQ fix and SSQ opt is illustrated in Figure 2. Given the practical range of the breakage parameters, the value for each parameter was searched with initial guess ranges of 0.2-1.5 for , 0.5-1.5 for α, 0.1-0.8 for , 0.5-1.3 for , and 0-10 for . The range of each variable was divided into 100 discrete points. Consequently, 10 , 10 , and 10 initial guess points were used for the WRIVS in three-, four-, and five-dimensional parameter spaces, respectively.

Experimental (Breakage Parameters from SSFBTs)
If the specific rate of breakage, , is time-independent, then, according to Equation (1), a log-linear relationship exists between ( )/ (0) and grinding time , i.e., the grinding is a first-order process. The first-order plots for the batch grinding tests are shown in Figure 3. The experimental results show that the grinding of both ores through three parameter values fixed. The parameters obtained from SSFBTs were used as the fixed values. The local minima and surface shape of SSQ fix were investigated over the two-dimensional space of the variables.
Then, we investigated the distribution of SSQ opt . The tendency of SSQ opt in the twodimensional parameter space was compared with the distribution of SSQ fix to investigate the sensitivity and convergence characteristics of the back-calculation process. The backcalculated parameter values that produced minimal SSQ opt values over a wide range of the initial guess were then compared to the values obtained from SSFBTs.
Finally, we increased the dimensions of back-calculation to between three and five and conducted calculations over a wide range of initial guess values. Then, we compared the back-calculated parameters producing minimal SSQ opt values to those obtained from SSFBTs.
Given the practical range of the breakage parameters, the value for each parameter was searched with initial guess ranges of 0.2-1.5 for A, 0.5-1.5 for α, 0.1-0.8 for ϕ, 0.5-1.3 for γ, and 0-10 for β. The range of each variable was divided into 100 discrete points. Consequently, 10 6 , 10 8 , and 10 10 initial guess points were used for the WRIVS in three-, four-, and five-dimensional parameter spaces, respectively.

Experimental (Breakage Parameters from SSFBTs)
If the specific rate of breakage, S, is time-independent, then, according to Equation (1), a log-linear relationship exists between w(t)/w(0) and grinding time t, i.e., the grinding is a first-order process. The first-order plots for the batch grinding tests are shown in Figure 3. The experimental results show that the grinding of both ores through ball milling is a first-order process. This indicates that temporarily evolving material properties and multi-particle interactions are negligible during the grinding times employed. Parameters A and α showing the best fit for each grinding case are listed in Table 2. A is higher for wet grinding by factors of 1.74 and 1.52 for molybdenum and uranium, respectively. This is owing to the higher dispersibility of particles and lower surface energy of minerals under wet conditions. The primary breakage distributions were obtained using the so-called BII method [2] and are shown in Figure 4. For uranium, the differences in the primary breakage distributions between dry and wet grinding were insignificant, as it is often presumed that B values do not change with grinding conditions. However, for molybdenum, we observed non-negligible differences in the breakage distributions between dry and wet conditions. The breakage distribution of the wet grinding lies below that of the dry grinding, indicating that the wet grinding produced less fine daughter particles on primary breakage. The best fitting values of ϕ, γ, and β used in Equation (4) are listed in Table 2. The primary breakage distributions were obtained using the so-called BII method [2] and are shown in Figure 4. For uranium, the differences in the primary breakage distributions between dry and wet grinding were insignificant, as it is often presumed that values do not change with grinding conditions. However, for molybdenum, we observed non-negligible differences in the breakage distributions between dry and wet conditions. The breakage distribution of the wet grinding lies below that of the dry grinding, indicating that the wet grinding produced less fine daughter particles on primary breakage. The best fitting values of , , and used in Equation (4) are listed in Table 2.  The simulated size distributions of the grinding product with time obtained from calculations using the values in Table 2 were compared with the experimental results, as illustrated in Figure 5. The size distributions obtained from the simulations show good agreement with the experimental results, indicating accurate estimation of the parameter values.  The simulated size distributions of the grinding product with time obtained from calculations using the values in Table 2 were compared with the experimental results, as illustrated in Figure 5. The size distributions obtained from the simulations show good agreement with the experimental results, indicating accurate estimation of the parameter values.

Back-Calculation: Error Function Distribution Analysis
For the parameter space with more than three dimensions, effective visualizatio the error function distribution is not possible. In this work, we first investigated the e function distribution in a two-dimensional parameter space, whereas the other three rameters remained fixed. Some of the results are shown in Figure 6. For all the pai

Back-Calculation: Error Function Distribution Analysis
For the parameter space with more than three dimensions, effective visualization of the error function distribution is not possible. In this work, we first investigated the error function distribution in a two-dimensional parameter space, whereas the other three parameters remained fixed. Some of the results are shown in Figure 6. For all the pairs of parameters, the error function was strictly convex with a single minimum. This means no local minima issue existed in the error function distribution of the BFPs in the two-dimensional parameter space. While visualizing the convexity of the error function distribution in higher dimensional spaces is not possible, based on the results and given that the possibility of a local minima occurring decreases with the increase in the dimension of the parameter space [37] we inferred that the existence of a local minima differing from the global minima of the error function distribution in the breakage function parameter space was unlikely.
For the (ϕ, β) pair in molybdenum grinding, we observed a flat regime of the error function without a significant single minimum. This could occur because parameters ϕ and β have complementary effects on function B. In general, α and γ are considered the indicators of inherent breakage characteristics that can be rarely complemented by changing other parameter values. Moreover, A is a constant factor that reflects the characteristic strength of the material under a given grinding environment, and it is rarely affected by changing other parameter values. However, as ϕ and β are mutually interfering parameters, similar primary breakage distribution functions can be obtained through combining different values. For uranium, strict convexity was observed even for the (ϕ, β) pair. This may have resulted from the inherent material-specific characteristics of the uranium ore. Under given conditions, mutual interference between ϕ and β is relatively weak compared to that observed in molybdenum.  The distribution of SSQ opt , the error at the point of convergence, was investigated. The values for SSQ fix were obtained from the fixed values of all the parameters, whereas SSQ opt was obtained through back-calculation where the coordinate parameters were optimized. We conducted back-calculation in the two-dimensional parameter space with varying initial guess values of two parameters, and some of the results are shown in Figure 7. Except for the (ϕ, β) pair in molybdenum grinding, the error at the point of convergence was a single value in all cases, irrespective of the initial guess values. The results for the (a, α) pair in dry uranium grinding are depicted in Figure 7a as a representative case. These results indicate that the back-calculation method successfully searched the optimal values of parameters regardless of the initial guess values. The results for the exceptional cases, i.e., the (ϕ, β) pair in wet and dry molybdenum grinding, are illustrated in Figure 7b,c. Within a regime that includes the minimum value point (Regime A in Figure 7b,c), the shape of the SSQ opt distribution was convex and was identical to the distribution of SSQ fix under the given initial guess values. Further, all the curves in Figure 7 are enlarged at a high magnification along the z-axis, and the surface of SSQ opt is extremely flat. Outside the regime, the error function distribution has a flat form. This indicates poor convergence of the back-calculation process owing to the flat surface of the error function distribution in the (ϕ, β) pair. In principle, although the shape of SSQ opt in the parameter space is convex, when the initial guess values are chosen inside regime A, back-calculation did not progress further owing to the flat surface of the error function. When the initial guess values were chosen outside the regime, the search process progressed until it reached the border of regime A, and then convergence stopped at the borderline. Therefore, the borderline of regime A can be regarded as the convergence limit of back-calculation. The BFPs of the grinding processes were back-calculated using the data in Figure 5. Table 3 compares the values of the back-calculated BFPs with the results from SSFBTs. For the case where convergence to a single value (Figure 7b,c) was not achieved, the parameters that produced the minimum SSQ values over a wide range were regarded as the The BFPs of the grinding processes were back-calculated using the data in Figure 5. Table 3 compares the values of the back-calculated BFPs with the results from SSFBTs. For the case where convergence to a single value (Figure 7b,c) was not achieved, the parameters that produced the minimum SSQ opt values over a wide range were regarded as the final back-calculated values. The back-calculated BFPs agreed well with the values obtained experimentally, showing differences of ±3% and indicating that back-calculation was successfully conducted. For the (ϕ, β) pair in molybdenum grinding, although convergence did not occur, the parameters at the point with the minimum SSQ opt agreed well with the SSFBT-derived parameters showing errors less than 3%. This indicates that we could successfully derive the breakage parameters through WRIVS in the two-dimensional parameter space, even when back-calculation failed to achieve convergence. We then extended the investigation to higher dimensional parameter spaces. The results are presented in Table 4. For two cases, i.e., those involving five and four dimensions with the parameter combination (A,α,ϕ,β), the resultant parameters differed considerably from the SSFBT-derived values. For those cases, SSQ opt was lower than that obtained through SSFBTs. This indicated that wrong combinations of parameters could produce similar SSQ opt owing to the high mutual interference between parameters. Except for the two cases mentioned above, the parameters obtained through WRIVS agreed well with the SSFBT-derived values with an acceptable degree of accuracy, indicating that the BFPs could be successfully derived through WRIVS. Several possible approaches are available to further enhance the utilizability of WRIVS. For example, one can conduct a set of SSFBTs to obtain credible parameters of the breakage distribution function, and then initiate the WRIVS-based back-calculation. Alternatively, it is possible to obtain the initial A and α values only by conducting a set of SSFBTs within the regime of sufficiently small particles in comparison with the ball size, in which case, the second term in Equation (2) equals zero and Equation (2) reduces to Equation (3). In both cases, the experimental number of necessary SSFBTs significantly decreases when compared with the conventional back-calculation strategy with no WRIVS. Furthermore, it may be possible to consecutively narrow down the initial value range by conducting the WRVS-based back-calculation in multi-stages. In this case, careful choice of the candidate regime from one stage may be necessary. This approach is not discussed in detail in the present work, and needs to be investigated in depth in a future work.

Conclusions
In this work, we investigated the feasibility of increasing the accuracy of the backcalculation method for determining BFPs through wide-range search. Using Korean molybdenum and uranium ores, we conducted SSFBTs and back-calculations, and visualized the distribution of the error function according to the parameter values. The results of this work can be summarized as follows: (1) In the two-dimensional parameter space, the error function was strictly convex, and no local minima differing from the global minima were observed. The main problem during back-calculation of the breakage parameters was weak convergence owing to the flat surface of the error function rather than occurrence of local minima. For molybdenum grinding, the shape of the error function for the (ϕ, β) pair was flatter compared to that for other parameter pairs. We attributed this to the strong mutual interference between these two parameters. (2) In the back-calculation involving the two-dimensional parameter space, for most pairs of variables, the optimization process converged to a single point and the results agreed well with those obtained through SSFBTs. However, for the parameter space with the pair producing a flat error function distribution, the solution of the back-calculation failed to converge. (3) Through wide-range search, even when convergence failed, the parameter combination that produced the minimum error agreed well with the parameters obtained from the SSFBTs. Thus, we could successfully determine the BFPs through wide-range search, irrespective of the shape of the error function. (4) Wide-range search proved feasible for searching in the three-dimensional parameter space, irrespective of the choice of parameters. For four-or five-dimensional parameter spaces, the accuracy of the back-calculation method decreased depending on the parameters that had to be back-calculated. (5) In conclusion, without proper initial values, searching seven breakage parameters through back-calculation only is not recommended owing to accuracy problems. When the number of parameters to be back-calculated is not greater than four, the parameters could be accurately identified through wide-range search, indicating that the number of SSFBTs required to determine the breakage function can be effectively reduced using the wide-range search.

Data Availability Statement:
The data presented in this study is available on request from the corresponding authors. The data is not publicly available due to the confidentiality of the research work.