**1. Introduction**

The size distribution belongs to the most critical properties of solid particulate materials, and particularly mixed solid waste, for example: The quality classes of solid recovered fuels (SRF) demand specific maximum particle sizes [1]. The particle size distribution (PSD) of the organic fraction of municipal waste impacts its anaerobic digestion [2]. The particle sizes of municipal solid waste influence the yields of dry gas, char, and tar in fixed bed reactor pyrolysis [3]. And the PSD influences the mass throughput of robotic sorters, which are limited by picks per hour [4]; hence, smaller particle sizes (and the corresponding smaller weights) decrease the possible mass throughput.

Concerning mixed commercial waste, besides the PSD's relevance as a technical quality criterion for processing products (e.g., SRF [5]), and its influence on the performance of reactors and processing machines (e.g., wind sifters [6]), different types of materials also concentrate in different particle size ranges (e.g., sorting analysis by Khodier et al. [7] and the size distribution of different plastic types according to Möllnitz et al. [8]). Hence, beyond influencing the shares of a plant's throughput that pass specific machines (due to material flow separation by screens), the PSD also determines the kinds of materials that pass through these machines.

Therefore, beneficial PSDs increase the effectiveness, as well as economic and ecologic efficiency of mixed solid waste treatment. Consequently, the PSD of mixed solid waste is

**Citation:** Khodier, K.; Sarc, R. Distribution-Independent Empirical Modeling of Particle Size Distributions—Coarse-Shredding of Mixed Commercial Waste. *Processes* **2021**, *9*, 414. https://doi.org/ 10.3390/pr9030414

Academic Editor: Aneta Magdziarz

Received: 8 February 2021 Accepted: 22 February 2021 Published: 25 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

deliberately influenced during mechanical processing, which is usually the first treatment stage for this kind of material, mainly through a combination of shredding and sieving [1–3,7,9–13].

### *1.1. Describing Particle Size Distributions*

Strictly speaking, the PSD of a collective is described as a list of the individual particles' sizes. But the representation as such a list is not suitable for analyzing and comparing PSDs. Consequently, various more useful methods for describing PSDs exist, which were summarized by Polke et al. [14]: Collectives of particles are often described through average equivalent diameters. An example is the Sauter diameter *dS*, which gives information on the specific surface of the collective (Equation (1), where *Vt* is the total volume of all particles, and *At* is the total surface area of all particles).

$$d\_S = 6V\_t/A\_t \tag{1}$$

Often, information on the width of the distribution is also essential. Consequently, measures of this width are frequently provided. Examples are the sample standard deviation *S* (Equation (2), where *di* is the size of the *i* th particle, *d* is the arithmetic average size, and *N* is the number of particles), or distribution-independent measures of the width, as shown, for example, in Equation (3) (where *W* is the width, and *db* is the *b*th percentile particle size).

$$S = \sqrt{\frac{\sum\_{i=1}^{N} \left(d\_i - \overline{d}\right)^2}{N - 1}}\tag{2}$$

$$\mathcal{W} = \frac{d\_{75} - d\_{25}}{d\_{50}} \tag{3}$$

Considering the described influence of the PSD on the path individual particles take through a plant, it is essential in mechanical waste processing to have more detailed knowledge on it than just summary values. Hence, the results of a PSD analysis are often reported graphically: the frequency density is shown in a histogram [15], where particle sizes are summarized into particle size classes (PSCs)—which is also the level of information obtained from sieve analyses (Figure 1a). Another representation, which is more suitable for comparing PSDs, is the sum distribution (Figure 1b) [14].

**Figure 1.** Representation of the overall particle size distribution (PSD) of a mixed commercial waste according to Reference [7]: (**a**) frequency density; (**b**) cumulative frequency density.

These graphical representations correspond to an empirical distribution [15]. As the sample size approaches infinity and the class width approaches zero, the histogram's representation becomes a continuous function: the probability density function (PDF). This PDF can also be approximated from analyses, where only summary information on discrete PSCs is available (e.g., sieve analyses), for example, through cubic splines [16] or Kernel density estimation [17].

Sometimes, the PDF can be approximately described by an analytical expression. In such cases, reporting the momentums of such an analytical distribution is sufficient to describe the PSD. So, for example, reporting the arithmetic mean particle size *d* of the sample and its standard deviation *S*, usually implies the underlying assumption of a normal distribution, according to Equation (4), where *q*(*d*) is the probability density or frequency density for particles of size *d*, *μ* is the arithmetic average of the population, which is estimated through *d*, and *σ* is the population's standard deviation, which is estimated through *S* [18].

$$q(d) = \frac{1}{\sqrt{2\pi\sigma^2}} \cdot \exp\left(-\frac{\left(d-\mu\right)^2}{2\sigma^2}\right) \tag{4}$$

Three further analytical PDFs, are reported as being relevant to the description of PSDs [14]: The log-normal distribution describes materials, where the logarithm of the particle size follows a normal distribution. It is a positively skewed distribution, which—in contrast to the normal distribution—only includes positive and, therefore, meaningful particle sizes. It is shown in Equation (5), where *μ* and *σ* are estimated by the arithmetic average and the sample standard deviation of the logarithm of the particle sizes.

$$q(d) = \begin{cases} 0 & \text{if } d < 0\\ \frac{1}{d\sqrt{2\pi\sigma^2}} \cdot \exp\left(-\frac{\left(\log(d) - \mu\right)^2}{2\sigma^2}\right) & \text{otherwise} \end{cases} \tag{5}$$

The Gates-Gaudin-Schuhmann (GGS) distribution [19] is an empirical approximation that is often suitable for describing products of coarse comminution processes. Its parameters are the maximum particle size *dmax* and the uniformity parameter *mu*. Its PDF is shown in Equation (6).

$$q(d) = \begin{cases} 0 & \text{if } d < 0 \text{ or } d > d\_{\text{max}}\\ \frac{m\_u}{d\_{\text{max}}} \left(\frac{d}{d\_{\text{max}}}\right)^{m\_u - 1} & \text{otherwise} \end{cases} \tag{6}$$

The Rosin-Rammler-Sperling-Bennet (RRSB) distribution [20] is also an empirical approximation. Its first parameter *d* is equal to the particle size, where the cumulative frequency density reaches the value 1 − 1/*e* ≈ 0.632, and its second parameter *nu* is a uniformity parameter. Its PDF is shown in Equation (7). The RRSB distribution is often used for describing products of fine comminution and dusts. Examples of the discussed distributions are shown in Figure 2.

$$q(d) = \begin{cases} 0 & \text{if } d < 0\\ \frac{n}{d^p} \left(\frac{d}{d^p}\right)^{n\_u - 1} \exp\left[-\left(\frac{d}{d^p}\right)^{n\_u}\right] & \text{otherwise} \end{cases} \tag{7}$$

**Figure 2.** Frequency density (**a**) and cumulative frequency density (**b**) of a normal, log-normal, Gates-Gaudin-Schuhmann (GGS), and Rosin-Rammler-Sperling-Bennet (RRSB) distribution.

### *1.2. Modeling Particle Size Distributions*

Products' PSDs result from their original condition and the kinds and parameters of machines that process them. Hence, to beneficially influence PSDs through process design and the choice and parametrization of machines, modeling and predicting them is desirable.

The most sophisticated and advantageous models are physical models, which provide an in-depth understanding of the phenomena that influence PSDs. Simulations based on such models are usually implemented using the discrete element method (DEM). Examples in literature are the works of Sinnott and Cleary [21] on the particle flows and breakage in impact crushers, Lee et al. [22] on breakage and liberation behavior of recycled aggregates from impact-breakage of concrete waste, and Dong et al. [23] on particle flow and separation on vibrating screens.

While DEM-based models improve process understanding, they also require high amounts of computational resources and detailed models and data on the processing machines and the materials to be comminuted, which limits their applicability in practice in many cases [24]: no published models, for example, incorporate the variability of materials, geometries and particle interactions for real mixed solid waste.

When physical models cannot be used, empirical regression models can deliver basic insights on machine and parameter influences on PSDs. For ensuring the reliability of the results, it is essential to involve statistical analyses for finding and interpreting the models. For scalar target values, the procedure has been thoroughly described by Khodier et al. [24], based on the example of the parametrization-dependent energy demand and throughput behavior of coarse shredders for mixed commercial waste.

PSDs are non-scalar: they are either described as PDFs, which are continuous functions, or as compositions of PSCs and, therefore, as multivariate vectors. Analytical PDFs, as those presented in Section 1.1., are defined by a set of momentums, which can also be treated as multivariate vectors. Consequently, modeling PSDs requires extensions of the methods used by Khodier et al. [24] to multivariate dependent variables.

The most widely applied variation of regression modeling is linear regression. It relates one or more dependent variables to one or more independent variables based on a set of linear regression coefficients. Its most general form—which is relevant to this work—is multivariate multiple linear regression, which involves multiple independent variables and multivariate dependent variables. Its model Equation is shown in Equation (8) [25]. The matrix *Y* (Equation (9)) is the *P* × *R* matrix of *P* observations of the *R*-variate vector

of the dependent variable, with elements *yp*,*r*. *ε* (Equation (12)) is also a *P* × *R* matrix and shows the corresponding model residuals *ε <sup>p</sup>*,*r*. *X* (Equation (10)) is a *P* × (*Q* + 1) matrix of the settings *xp*,*<sup>q</sup>* of the *Q* independent variables corresponding to the *P* observations. Its first column (which is indexed with 0) is a column of ones, corresponding to constant terms in linear regressions. And the matrix *β* (Equation (11)) contains the (*Q* + 1) × *R* regression coefficients *βq*,*r*.

$$\mathbf{Y} = \mathbf{X} \cdot \mathbf{\mathcal{B}} + \varepsilon \tag{8}$$

$$\mathbf{Y} = \begin{bmatrix} y\_{1,1} & y\_{1,2} & \cdots & y\_{1,R} \\ y\_{2,1} & y\_{2,2} & \cdots & y\_{2,R} \\ \vdots & \vdots & \ddots & \vdots \\ y\_{P,1} & y\_{P,2} & \cdots & y\_{P,R} \end{bmatrix} \tag{9}$$

$$\mathbf{X} = \begin{bmatrix} \mathbf{x}\_{1,0} & \mathbf{x}\_{1,1} & \mathbf{x}\_{1,2} & \cdots & \mathbf{x}\_{1,Q} \\ \mathbf{x}\_{2,0} & \mathbf{x}\_{2,1} & \mathbf{x}\_{2,2} & \cdots & \mathbf{x}\_{2,Q} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{x}\_{P,0} & \mathbf{x}\_{P,1} & \mathbf{x}\_{P,2} & \cdots & \mathbf{x}\_{P,Q} \end{bmatrix} = \begin{bmatrix} 1 & \mathbf{x}\_{1,1} & \mathbf{x}\_{1,2} & \cdots & \mathbf{x}\_{1,Q} \\ 1 & \mathbf{x}\_{2,1} & \mathbf{x}\_{2,2} & \cdots & \mathbf{x}\_{2,Q} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \mathbf{x}\_{P,1} & \mathbf{x}\_{P,2} & \cdots & \mathbf{x}\_{P,Q} \end{bmatrix} \tag{10}$$

$$
\mathcal{B} = \begin{bmatrix}
\beta\_{0,1} & \beta\_{0,2} & \cdots & \beta\_{0,R} \\
\beta\_{1,1} & \beta\_{1,2} & \cdots & \beta\_{1,R} \\
\vdots & \vdots & \ddots & \vdots \\
\beta\_{Q,1} & \beta\_{Q,2} & \cdots & \beta\_{Q,R}
\end{bmatrix} \tag{11}
$$

$$
\mathfrak{e} = \begin{bmatrix}
\varepsilon\_{1,1} & \varepsilon\_{1,2} & \cdots & \varepsilon\_{1,R} \\
\varepsilon\_{2,1} & \varepsilon\_{2,2} & \cdots & \varepsilon\_{2,R} \\
\vdots & \vdots & \ddots & \vdots \\
\varepsilon\_{P,1} & \varepsilon\_{P,2} & \cdots & \varepsilon\_{P,R}
\end{bmatrix} \tag{12}
$$

The resulting model is obtained by minimizing the sum of squares of the residuals *<sup>ε</sup> <sup>p</sup>*,*r*. It is shown in Equation (13), where *<sup>Y</sup>*<sup>ˆ</sup> is the matrix of the *<sup>P</sup>* <sup>×</sup> *<sup>R</sup>* model predictions *<sup>y</sup>*ˆ*p*,*<sup>r</sup>* for the dependent variable, corresponding to the observations in *Y*, and *β*ˆ is the matrix of the (*<sup>Q</sup>* <sup>+</sup> <sup>1</sup>) <sup>×</sup> *<sup>R</sup>* least-squares estimates *<sup>β</sup>*<sup>ˆ</sup> *<sup>q</sup>*,*<sup>r</sup>* of the regression coefficients in *<sup>β</sup>*.

$$
\hat{Y} = X \cdot \mathcal{J} \tag{13}
$$

For linear regression models, it is necessary to describe the dependent variable as a vector of a fixed length. In the case of analytical PDFs, this is the vector of the momentums. So, the immediate result of a model for these momentums is the vector of their expected values, and the corresponding confidence bands. From there, the PDF and its confidence region can be calculated.

Theoretically, *β*ˆ can also be calculated from *Q* univariate linear regressions. The resulting values *β*ˆ *<sup>q</sup>*,*<sup>r</sup>* are identical. But multivariate methods, involving multivariate linear regression and multivariate analysis of variance (MANOVA, e.g., Reference [26]), are preferable: they allow a more accurate calculation of confidence regions, considering correlations between the momentums. Moreover, the evaluation of parameter's significance should be based on multivariate considerations to find coherent models for the resulting interdependent—distribution of the particle sizes.

Analytical PDFs are practical, as they allow a detailed description of PSDs with a small number of variables—the momentums. However, only few processes produce PSDs that follow analytical functions, according to Polke et al. [14]. Consequently, methods for distribution-independent modeling of PSDs are also needed, which is especially true for the mechanical processing of mixed solid waste, considering the variety of processing machines used there: e.g., shredders, screens, magnetic separators, and sensor-based sorters [10].

Empirical distributions allow the distribution-independent description of PSDs: particle sizes are amalgamated to *D* discrete PSCs. As a result, the PSD is described as a composition—a *D*-dimensional vector of the shares of those PSCs.

Mathematically, the compositional nature of the vector of PSCs has wide-reaching consequences: *D*-dimensional compositions are constrained to a vector space called the *D*dimensional simplex *S*(*D*) , which is a sub-space of the *D*-dimensional real space *R*(*D*) [27]. Compositions, being simplicial vectors, have specific common properties: Their *D* parts are not linearly independent. This results from a summation constraint: their parts sum up to a fixed constant, e.g., 1 or 100%. Furthermore, parts may only have positive values or zero. So, more precisely, *S*(*D*) is a sub-space of *R*+(*D*) **<sup>0</sup>** . The most well-known representation of a simplex is the ternary diagram (see, e.g., Reference [28]), which shows the threedimensional simplex.

When modeling simplicial vectors, the constraints of the simplex and the interdependence of the compositional parts must be considered. While Khodier et al. [29] report the empirical observation that the constraints are automatically fulfilled for the prediction values *Y*ˆ, if all observations in *Y* are valid compositions, this does not apply to corresponding confidence regions. Hence, different methods were developed in the past decades to handle and model simplicial data.

The most widely applied approach for handling compositions is transforming data using log-ratios, as suggested by Aitchison [30]. There are many kinds of such log-ratios, but the state of the art approach in the compositional data community is the application of so-called isometric log-ratios (ilr) (as proposed by Egozcue [31]), according to Reference [32]. These are bijective projections of the *D*-dimensional Simplex onto a (*D* − 1)-dimensional real-space with an orthonormal basis [27]. The resulting ilr-coordinates are unconstrained, linearly independent coordinates. Hence, standard statistics can be applied in the projected log-ratio space, as demonstrated, for example, by Edjabou et al. [33] for waste composition analysis. Consequently, the use of ilr-transformations allows the application of standard multivariate multiple linear regression to predict PSDs, described as empirical distributions.

In this work, the approach of modeling influences on empirical PSDs using multivariate multiple linear regression and ilrs is applied on particle size data, using the programming language R. The data was obtained within Khodier et al.'s [24] industry-scale coarse-shredding experiments with real mixed commercial waste. Based on this example, this work aims at presenting the method to the process engineering and waste processing communities, enabling the distribution-independent empirical modeling of particle size distributions, while preserving all relevant information. Moreover, the method's suitability for PSD modeling and its limitations and potential pitfalls in the interpretation of the results are discussed. And, finally, insights on the influences of coarse shredders' parameters on the PSDs of mixed commercial waste are presented, complementing findings of Khodier et al. [24] on their effects on shredders' throughput behavior and energy demand.

### **2. Materials and Methods**

### *2.1. Experimental Design and Setup*

The choice of the experimental design and the shredding experiment setup have been explained in detail by Khodier et al. [24]. Hence, they are only summarized in the following. The extension of the experiment with material sampling and particle size analysis are described in detail.

### 2.1.1. Experimental Design

The shredding experiment examines the influence of three independent variables on the throughput behavior and energy demand of a Terminator 5000 SD, which is a singleshaft shredder from the Austrian company Komptech GmbH (Frohnleiten, Austria)—a research partner in the funded project ReWaste 4.0. These independent variables are the radial gap width *w*, the shaft rotation speed *s*, and the cutting tool geometry *c*. The factor range of *w* was defined from 0% to 100% of the maximum gap width, with discrete levels

at a step size of 10%. The minimum of the factor range of *s* was chosen at 60% and the maximum at 100% of the maximum shaft rotation speed of 31 rpm, again width discrete 10% steps. Concerning *c*, three different geometries are examined, called "F", "XXF", and "V" (for more details, cf. Reference [24] and Figure A1 and Table A1 in the Appendix A).

The factors are coded for the design and analysis of the experiment: for the numerical factors *s* and *w*, their range is adjusted to −1 to 1, representing the minimum and maximum factor settings. Concerning the nominal factor *c*, it is represented by two variables *c*<sup>1</sup> and *c*2, based on sum contrasts (cf. Reference [34]). The values of these variables that correspond to the cutting tool geometries are shown in Table 1.

**Table 1.** Contrast matrix for the cutting tool geometry.


The experimental design settings of the independent variables were chosen based on a statistical Design of Experiments (cf. Reference [35]). A 32 runs, completely randomized D-optimal design (cf. Reference [36]) was chosen, with no blocking (cf. Reference [37]) and with five replicate points and five lack-of-fit points. It is based on the reduced-cubic design model, shown in Equation (14), where *y*ˆ(*r*) is the model prediction for the *r*th (univariate) response (=dependent variable), - *x* is a vector of the factors' settings, and *K*(*r*) *jkmn* is the model coefficient for the *r*th response and the factor or interaction (=multiplication of factors) *wj skcm* <sup>1</sup> *<sup>c</sup><sup>n</sup>* 2 .

$$\mathcal{G}^{(r)}\left(\stackrel{\rightarrow}{\mathbf{x}}\right) = \sum\_{j=0}^{2} \sum\_{k=0}^{2-j} \sum\_{m=0}^{1} \sum\_{n=0}^{1-m} \left( K\_{jkmn}^{(r)} w^j s^k c\_1^m c\_2^n \right) \tag{14}$$

Equation (14) can easily be extended to multivariate responses; therefore, the design is also valid for multivariate multiple linear regression. To extend it, *y*ˆ(*r*) equals *y*ˆ*p*,*<sup>r</sup>* in Equation (13): the *r*th univariate response becomes the *r*th dimension of the multivariate response, and *<sup>p</sup>* indexes the corresponding vector of factor settings - *x* , which is simply a row of *X*. Each factor or interaction *w<sup>j</sup> skcm* <sup>1</sup> *<sup>c</sup><sup>n</sup>* <sup>2</sup> is represented by a column of *X*, and each coefficient *K*(*r*) *jkmn* corresponds to a coefficient *<sup>β</sup>*<sup>ˆ</sup> *<sup>q</sup>*,*r*.

### 2.1.2. Setup of the Shredding Experiment

The flowchart of the experiment is shown in Figure A2 in the Appendix A: The feed material is waste, declared as mixed commercial waste (cf. Reference [24] for more details), collected in Styria in Austria in October 2019. It is fed into the shredder's feeding bunker using a wheel loader. From the shredder's output belt, the material is passed to a digital material flow monitoring system (DMFMS), consisting of a belt-scale and optical sensors (cf. Reference [11]). The material leaving the DMFMS is collected on a product heap. Each experimental run has a total duration of one hour.

### 2.1.3. Sampling

For analyzing the PSD of the shredded waste, samples must be taken. The design of the sampling process in this experiment is based on Pierre Gy's Theory of Sampling (TOS), as described in the Danish standard DS 3077 [38] and the work of Khodier et al. [7] on its application on coarsely shredded mixed commercial waste.

According to TOS, the fundamental sampling principle must be considered: each particle must have the same probability of ending up in the sample. The most beneficial sampling situation is a one-dimensional sampling [39], e.g., taking the samples from a falling stream. Furthermore, a sample should spatially cover the whole lot, which is achieved by composite sampling: the sample is a composite of so-called increments (the sampled material from one individual sampling step).

In this work, the sample for each experimental run consisted of 40 such increments. They were taken from a falling stream, swiveling the output conveyor belt of the DMFMS back and forth over a sampling box, elevated by a forklift (see Figure 3), every 3 minutes, taking two increments. The inner dimensions of the box were 115 × 915 × 565 (length × width × height in mm). For ensuring that all desired material ends up in the box, the height of the back-side wall of the box was increased by placing 1.5 m-long wood boards inside. The box was changed after every ten increments. Based on the share of the samples in the total processed mass, each increment covered approximately 0.61 seconds of throughput. For intermediate storage, the samples were finally transferred to 1 m<sup>3</sup> big bags.

**Figure 3.** Sampling setup.

### 2.1.4. Particle Size Analysis

The PSDs of the samples were analyzed using a Komptech (Frohnleiten, AT) Nemus 2700 drum screen and five screening drums that have square holes with side lengths of 80, 60, 40, 20, and 10 mm. The geometries of the drums are shown in Figure A3 and Table A2 in the Appendix A.

First, one big bag of a sample was evenly distributed in the feeding bunker of the screen, which has a length of 4033 mm in the direction of the material flow, and a width of 1035 mm. Then, the drum was started with a rotation speed of 11.5 rpm, and the conveyor belt of the feeding bunker was operated at 0.026 m/s. The drum screen was only stopped after all material had passed. The produced fine fraction was then screened with the subsequent finer drum. The scale used for measuring the masses of the fraction has an uncertainty of 10 g.

### *2.2. Analysis of the Results*

The analyses of the results, which are explained in this section, were performed in R version 4.0.2, based on the work of van den Boogaart and Tolosana-Delgado [40]. The implementation is attached as an HTML export of a jupyter notebook (see Supplementary Material).

For the analyses in this work, the six particle size fractions from the particle size analysis are aggregated to three PSCs for easier visualization. To ensure the relevance of the findings for mechanical waste processing, the PSCs were chosen, based on the particle size limits of SRF premium quality (30 mm) and SRF medium quality (80 mm) [1]. Since none of the used screen drums had a mesh width of 30 mm, equal shares of the screening fraction 20–40 mm are assigned to the particle size fractions 0–30 mm and 30–80 mm, which corresponds to linear interpolation.

### 2.2.1. Isometric Log-Ratios

The ilr transformation (denoted as a function ilr()) is an isometry of the vector spaces *SD* and *RD*−*<sup>1</sup>* [27], which means that the distance between two compositions <sup>→</sup> *y* (1) and → *y* (2) is preserved in the transformation. For interpretation, it is essential to understand that the distance referred to is not the Euclidian distance Δ*<sup>E</sup>* (Equation (15), where *yi* is the *i* th element of a *R*-dimensional composition <sup>→</sup> *y* ), but rather the Aitchison distance Δ*<sup>A</sup>* (Equation (16), where ilr*i*() is the function for the *i* th ilr dimension) (cf. Reference [27]). Consequently, the preserved distance is of a relative, multiplicative nature and not an absolute, additive nature. Hence, the least-squares minimization of the model residuals in *ε*, when calculating the coefficients' estimates *β*ˆ**,** is also based on the Aitchison distance when applying the ilr transformation. Considering this is particularly important when the order of magnitude of the parts' shares differs significantly since small absolute differences become very significant on a relative scale for very small shares (cf. Reference [41]).

$$
\Delta\_E \left( \stackrel{\rightarrow}{\mathcal{Y}}^{(1)}, \stackrel{\rightarrow}{\mathcal{Y}}^{(2)} \right) = \sqrt{\sum\_{i=1}^{R} \left( y\_i^{(1)} - y\_i^{(2)} \right)^2} \tag{15}
$$

$$\Delta\_A \left( \stackrel{\rightarrow}{y}^{(1)}, \stackrel{\rightarrow}{y}^{(2)} \right) = \sqrt{\frac{1}{2R} \sum\_{i=1}^{R} \sum\_{j=1}^{R} \left[ \ln \left( \frac{y\_j^{(1)}}{y\_j^{(1)}} \right) - \ln \left( \frac{y\_j^{(2)}}{y\_j^{(2)}} \right) \right]^2} = \sqrt{\sum\_{i=1}^{R-1} \left[ \text{ilr}\_i \left( \stackrel{\rightarrow}{y}^{(1)} \right) - \text{ilr}\_i \left( \stackrel{\rightarrow}{y}^{(2)} \right) \right]^2} \tag{16}$$

Furthermore, the ilr transformation is not defined if any part has a value of zero. While there are different approaches to handling such values, they complicate the application of the transformation [27].

In the following, ilr-transformed coordinates are marked with "∗" so that *Y*<sup>∗</sup> = ilr(*Y*). The back-transformation function is denoted "ilr−1()". For calculating ilr coordinates, the compositional parts are sequentially grouped: first, each component is assigned to a group +1, or −1. For subsequent ilr coordinates, the elements of one group are again assigned to new groups +1 or −1, and the parts of the other group are assigned to group 0. Each ilr coordinate is then calculated according to Equation (17), where *ai*,*<sup>r</sup>* is the scaling factor for the *r*th compositional part and the *i* th ilr coordinate. The calculation of *ai*,*<sup>r</sup>* is shown in Equation (18), where *t* is the number of parts in group +1, and *v* is the number of parts in group −1 [31].

$$\text{ilr}\_i \left( \stackrel{\rightarrow}{y} \right) = \sum\_{r=1}^{R} [a\_{i,r} \ln(y\_r)] \tag{17}$$

$$a\_{i, \mathcal{I}} = \begin{cases} +\frac{1}{\mathcal{I}} \sqrt{\frac{tv}{t+v}} & \text{for parts of group } +1\\ -\frac{1}{v} \sqrt{\frac{tv}{t+v}} & \text{for parts of group } -1\\ 0 & \text{for parts of group } 0 \end{cases} \tag{18}$$

Greenacre [42] documents concerns regarding the interpretation of ilr-transformed data since it incorporates the geometric means of compositional parts. Consequently, in this work, results are interpreted based on graphical representations of back-transformed data. Hence, the exact choice of groups is arbitrary, and the standard grouping of the "ilr()" function of the "compositions" package version 2.0-1 in R [43] is used. In this work, *r* = 1 stands for the fraction > 80 mm, *r* = 2 for the fraction 30–80 mm, and *r* = 3 for the fraction 0–30 mm (see Supplementary Material). Resulting from the standard grouping, the ilr-transformed representation <sup>→</sup> *y* ∗ of a particle size composition <sup>→</sup> *y* (corresponding to a row of *Y*) is calculated according to Equation (19).

$$\stackrel{\rightarrow}{y}^\* = \text{ilr}\left(\stackrel{\rightarrow}{y}\right) = \left(\begin{array}{c} \ln\left[\left(\frac{y\_2}{y\_1}\right)^{\sqrt{0.5}}\right] \\\\ \ln\left[\frac{y\_3}{\left(y\_1 y\_2\right)^{\sqrt{1/6}}}\right] \\\\ \end{array}\right) \tag{19}$$

### 2.2.2. Model Reduction: MANOVA

To obtain a final, reliable model, the factors and interactions in Equation (14) must be checked on their significance, eliminating non-significant ones, but retaining model hierarchy (cf. Reference [24]). For univariate dependent variables, this is done through Ftests in an analysis of variance (ANOVA) (cf. Reference [35]). The multivariate character of the compositional dependent variable in this work requires a multivariate extension of the ANOVA: the MANOVA. Different from the ANOVA, there is more than one definition of the F-statistic in the MANOVA. The most commonly used definitions are the Pillai-Barlett trace, Wilk's lambda, the Hotelling-Lawley trace, and Roy's largest eigenvalue statistic, according to Hand and Taylor [26]. These are also the ones implemented in the "regr" package, version 1.1 [44], used for the MANOVA in this work.

For the analyses at hand, the Pillai-Barlett trace was chosen, based on the recommendation of Olson [45] as cited in Reference [26]. Model reduction is performed applying backward selection: the least significant term, which can be removed without violating model hierarchy, is eliminated, as long as removable non-significant factors or interactions are present (cf. Reference [40]). Analogous to Reference [24], 0.1 is chosen as the threshold so that factors and interactions with an empirical significance (*p*-value) higher than that threshold are discarded. The relevant *p*-values are calculated using the "drop1()" function from the "regr" package (version 1.1).

For evaluating the final model, three performance values are calculated: The coefficient of determination *R*<sup>2</sup> calculates how much of the variance of the data is explained by the model. The adjusted coefficient of determination *R*<sup>2</sup> *adj* is a measure similar to *<sup>R</sup>*2. But it is adjusted by the terms in the model and thereby evaluates the model's efficiency [35]. Both are calculated using the "R2()" function in R.

The prediction coefficient of determination *R*<sup>2</sup> *pred* determines the share of variance, which is explained by models fitted without considering the very point which is evaluated. High differences between *R*<sup>2</sup> *pred* and *<sup>R</sup>*<sup>2</sup> *adj* indicate overfitting. *<sup>R</sup>*<sup>2</sup> *pred* is calculated according to Equation (20), where *PRESS* is the prediction residual sum of squares and *SStot* is the total sum of squares [46]. *PRESS* is calculated, using the "PRESS()" function from the "MPV" package version 1.56 [47]. And *SStot* is calculated, according to Equation (21). *y*∗ *p*,*r* is the *<sup>p</sup>*th observation of the *<sup>r</sup>*th of (*<sup>R</sup>* − <sup>1</sup>) coordinates of the ilr-transformed dependent variable. And *y*∗ *<sup>r</sup>* is the arithmetic mean of the observations of the *<sup>r</sup>*th ilr coordinate (see Equation (22)).

$$R\_{prcl}^2 = 1 - \frac{PRESS}{SS\_{tot}}\tag{20}$$

$$SS\_{tot} = \sum\_{p=1}^{P} \sum\_{r=1}^{R-1} \left( y\_{p,r}^\* - \overline{y}\_r^\* \right)^2 \tag{21}$$

$$\overline{y}\_r^\* = \frac{1}{P} \sum\_{p=1}^P y\_{p,r}^\* \tag{22}$$

2.2.3. Analysis of the Residuals

The tests in the MANOVA require multivariate normality of the (ilr-transformed) residuals *ε*∗. Hence, to validate the final model, the distribution of the residuals must be

examined. Each coordinate of variables, which follow a multivariate normal distribution, also follows a univariate normal distribution [48]. Hence, a quantile-quantile plot for each coordinate is examined as a first visual step. Since the individual coordinates' univariate normality is a necessary, but not sufficient condition for multivariate normality, multivariate tests are also applied, particularly Mardia's Skewness and Mardia's Kurtosis [49]. Both tests are part of the "mvn()" function in "MVN" package version 5.8 [50], which also tests the univariate normality of the individual coordinates using the Shapiro-Wilk test.

### 2.2.4. Confidence and Prediction

The resulting model Equation shows the most likely prediction value. Two regions express the uncertainty of this prediction: the confidence region and the prediction region (cf. Reference [40]). The confidence region covers the uncertainty of the model parameter estimation. The region reflects likely average PSC distributions (on an ilr scale) for specific parameter settings, for extended operation times, on a chosen confidence level. The prediction region adds the residual variability around the expected value. Hence, it shows likely PSC distributions for one hour of operation (since this is the experimental duration the data is based on).

Due to the required multivariate character of the ilr-transformed residuals, the resulting confidence and prediction regions are equipotential-ellipses (resulting from the PDF of the multivariate normal distribution) on an ilr-scale. Van den Boogaart and Tolosana Delgado [40] provide R-code for calculating these regions and visualizing their backtransformed representation in ternary diagrams. This code is used in this work (see Supplementary Materials).

### **3. Results and Discussion**

### *3.1. Data and Model*

The experimental design and the resulting shares of the PSCs are shown in the Supplementary Materials. Their order corresponds to the order in which the experimental runs were performed. Originally, a completely randomized order was planned. As reported by Khodier et al. [24], due to an unintentional change of the motor rotation speed of the V cutting tool during the experiments, three runs had to be repeated. Since the tight timescale did not allow re-randomizing all remaining runs, considering the time consumed when switching shredders, the randomness is slightly impaired.

A significant model was found based on the experimental data. It is visualized in Figure <sup>4</sup> and shown in Equation (23), where <sup>→</sup><sup>ˆ</sup> *y* is a vector showing the model predictions of the mass shares of the three PSCs. Its first element corresponds to the coarse fraction (>80 mm), its second element to the medium fraction (30–80 mm), and the third element shows the fine fraction (0–30 mm). As shown in the Equation, the only significant factor (at a threshold of *p* = 0.1) is the cutting tool geometry *c*.

$$\stackrel{\circ}{y} = \text{i}\text{l}\text{r}^{-1}\left[\begin{pmatrix} 0.137\\0.047 \end{pmatrix} - \begin{pmatrix} 0.307\\0.178 \end{pmatrix} c\_1 + \begin{pmatrix} 0.344\\0.145 \end{pmatrix} c\_2\right] = \begin{cases} \begin{pmatrix} 0.39 & 0.31 & 0.30\\ \begin{pmatrix} 0.30 & 0.35 & 0.35 \end{pmatrix}^T & \text{for } c = \text{"XF"}\\ \begin{pmatrix} 0.21 & 0.42 & 0.37 \end{pmatrix}^T & \text{for } c = \text{"XF"} \end{cases} \tag{23}$$

The analysis of the ilr-residuals proves a multivariate normal distribution, with *p*values of 0.67 and 0.89 for the Shapiro-Wilk test on univariate normality of each coordinate and *p*-values of 0.89 and 0.71 for Mardia's Skewness and Mardia's Kurtosis, respectively.

*R*<sup>2</sup> for the model is 0.57 and *R*<sup>2</sup> *adj* is 0.54. These values are noticeably lower than those of the models for the throughput behavior and energy demand (see Reference [24]), which range from 0.73 to 0.87 for *R*<sup>2</sup> and from 0.67 to 0.81 for *R*<sup>2</sup> *adj*.

The lower coefficients of determination are reasonable since material fluctuations are likely to influence the material quality more than the process behavior, and sampling adds random noise (cf. Reference [7]). Considering the high expectable noise in processing

experiments with real waste, and especially in shredding experiments (cf. Reference [24]), the model performance is satisfactory.

**Figure 4.** Prediction values and confidence and prediction regions for the particle size class distributions resulting from different cutting tool geometries.

### *3.2. Discussion of the Method*

The data from the particle size analyses were amalgamated to three discrete PSCs and then freed from the restrictions of the simplex by applying an ilr-transformation. The resulting model, the multivariate normal distribution of its residuals, and the consequent successful calculation of confidence and prediction regions prove the potential of this approach in terms of distribution-independent, mathematically correct empirical modeling of particle size distributions.

It is another potential application of ilr-transformations in the context of waste management, besides the application on waste compositions, in terms of sorting fractions (cf. Reference [33]). The method is a notable extension to the toolbox of chemical and process engineers in general, besides familiar approaches, like equivalent diameters or certain analytical distributions (see Section 1.1).

The discussed restriction of the method concerning zero-values is not an issue for the present data. When data include zeros, the potential impact of zero-handling techniques must be considered. These include amalgamations of compositional parts or zeroreplacement (cf. Reference [42]).

Furthermore, the impact of the relative scale of ilr-transformed data must be kept in mind. Concerning the aim of the present work, it cannot be ignored, since the absolute

values of the PCSs' shares result in absolute waste masses, which are of interest. But the effect of the relative scale is low when the shares of all parts are of similar orders of magnitude, as is the case in this work (see Equation (23)). For other cases, calculations of weighted variances (cf. Reference [42]), for example, can be used to counteract the impact of the relative scale, if necessary. The relative scale can also be beneficial in other cases (as explained by Pawlowsky-Glahn et al. [27]), for example, where the share of trace elements in a chemical composition has a high impact.

The discussed issue of interpreting ilrs, is solved by graphical representation in this work. While this is a straightforward approach for three or fewer fractions, the question of how to represent more fractions arises. Graphic solutions include sets of ternary diagrams of two specific fractions and amalgamations of the others (cf. Reference [40]) and area plots (for the expected values, but not for confidence and prediction regions). Some other approaches include: The application of easier interpretable log-ratios (e.g., additive logratios or log-ratios incorporating amalgamations, cf. Reference [42]) at the cost of the exact representation of the variance structure in the data, or modeling non-transformed data, while incorporating potential issues with the restrictions of the simplex, particularly for calculating confidence and prediction regions.

Finally, the results in Figure 4 confirm the decision of Khodier et al. [24] to perform a Design of Experiments-based investigation, with multiple runs: the 95% prediction regions overlap even for the F and V unit. Hence, conclusions that contradict the findings in this work could be drawn when comparing only single runs of one hour.

### *3.3. Discussion of the Modeling Results*

For the gap width and shaft rotation speed, no significant impacts on the produced PSDs were identified. In conclusion, either no such effects exist, or they are too small, compared to the residual variance in the data, to be detected based on the data at hand (resulting in so-called type II or β error). According to Biemann [51], no matter how small, any effect becomes significant if the amount of data is big enough. Hence, the order of magnitude of potential, non-identified effects is of interest. For the chosen limit *p*-value of 0.1 in the MANOVA, linear effects are preserved in the model if the 90% confidence regions of their extreme settings (which get smaller with more data) do not overlap. And potential effects are likely not to exceed the maximum distance of the borders of these regions.

Consequently, the confidence regions for the minimum and maximum settings of *w* and *s* at factor setting 0 for the other variables were plotted (see Figure 5), based on a linear model (see Equation (24)), to give a visual impression of potential type II errors. The residuals of these models also follow a multivariate normal distribution. Figure 5 shows, that *w* influences the allocation to PSCs of 0 to about 13% of the product material, and *s* influences 0 to about 9%, for average settings of the other factors.

$$\stackrel{\circ}{y} = \text{ilr}^{-1}\left[ \begin{pmatrix} K\_{0000}^{(1)} \\ K\_{0000}^{(2)} \end{pmatrix} + \begin{pmatrix} K\_{1000}^{(1)} \\ K\_{1000}^{(2)} \end{pmatrix} w + \begin{pmatrix} K\_{0100}^{(1)} \\ K\_{0100}^{(2)} \end{pmatrix} s + \begin{pmatrix} K\_{0010}^{(1)} \\ K\_{0010}^{(2)} \end{pmatrix} c\_1 + \begin{pmatrix} K\_{0001}^{(1)} \\ K\_{0001}^{(2)} \end{pmatrix} c\_2 \right] \tag{24}$$

Differences in the PSD, caused by the shaft rotation speed, would most likely be based on differences in the breakage of brittle materials in the waste, caused by the impact speed of the shaft's teeth. Since the used machines are slow running shredders, the non-existence of such effects appears feasible. Furthermore, since potential impacts of this kind are most likely relatively small, the low share of brittle materials in mixed commercial waste (cf. Reference [7]) complicates their identification.

Concerning the gap width, considering the geometry of the cutting tools, a slight increase of the coarse fraction (>80 mm) would be reasonable due to falling through of uncomminuted particles between the teeth of the counter comb. But the non-significance of potential effects makes sense, considering that the breakage situation, according to Feyerer [52], does not change with the gap width for the F and XXF geometry and only slightly for the V geometry.

The influence of the cutting tool geometry is highly significant, with a factor *p*-value <10−5. As Figure 4 shows, it is largest, comparing the PSDs produced by the F and the V geometry. But the F and XXF geometries also differ significantly on a 95% confidence level. The results confirm expectations, considering the geometries (cf. Figure A1 and Table A1 in the Appendix A): The smaller axial gap between the counter comb teeth of the XXF geometry leads to a finer product, compared to the F unit. To be more precise, the share of the coarse fraction decreases in favor of the two other fractions.

Concerning the V geometry, the gap between the counter comb teeth is smaller than the gap of the F geometry. It is larger than the XXF geometry's gap close to the shaft but gets much smaller with increasing distance. This smaller gap, combined with a comb system (cf. Reference [24]), leads to the finest product among the examined geometries—again, mainly in terms of an even lower share of coarse material compared to the XXF geometry.

Relating the results to the findings of Khodier et al. [24], the choice of a cutting tool geometry depends on the requirements of the process: The V tool produces the finest material of the three but at the cost of higher energy demand and smaller but steadier throughput.

Concerning the gap width and the shaft rotation speed, the standard operation with minimum gap widths and maximum shaft rotation speeds must be questioned: increasing the gap width is beneficial for the throughput and energy demand and hardly affects the throughput steadiness, according to Reference [24]. The shares of the PSCs chosen in this work (typical for SRF production) are not or only a little affected, with a maximumlikelihood influence on only about 3% of the material.

For the shaft rotation speed, the maximum likelihood influence of this nonsignificant factor on the PSD only concerns about 2% of the material, while the mass flow and energy demand show an optimum at about 84% and 80% of the maximum shaft rotation speed, respectively.

**Figure 5.** Ninety percent confidence regions for the minimum (−1) and maximum (+1) settings of *w* (**a**) and *s* (**b**), at factor setting 0 for the corresponding other factors of *w*, *s*, *c*<sup>1</sup> and *c*2, based on the linear model from Equation (24).

### **4. Conclusions**

Multivariate multiple linear regression modeling was applied in this work on ilrtransformed PSC data from a Design of Experiments-based 32 runs coarse-shredding experiment with mixed commercial waste. A significant model, with an *R*<sup>2</sup> of 0.57 was found, identifying the cutting tool geometry as a highly significant influence on the PSD.

The gap width and shaft rotation speed were found not to be significant, with maximum-likelihood influences on 3% and 2% on the material, respectively. If the discussed potential type II errors are rated as economically relevant or other particle size classes are of interest, further data should be generated and analyzed. Otherwise, the PSD can be treated as invariant to these factors when optimizing the throughput and energy demand. Consequently, based on the new insights from this work, a much more efficient operation of mechanical waste processing plants can be reached. The influence of the cutting tool, on the contrary, is highly significant. Its choice depends on process requirements.

What was not investigated are selective influences on specific material fractions, e.g., metals or wood. The investigation of the PSDs of such fractions may give more detailed insights and lead to different conclusions on process parametrization. It is subject to further research, which can make use of the presented modeling methods.

Using the presented method requires an in-depth understanding of the implications of applying the ilr transformation. It is essential to avoid wrong interpretations, caused by the introduced relative scale or by zero-replacement practices, where necessary. Nonetheless, establishing such understanding is rewarding: in conclusion, the method has proven to be suitable for the distribution-independent modeling of PSDs. Hence, it is a valuable addition to the toolkit of engineers, dealing with particulate materials.

**Supplementary Materials:** The following is available online at https://www.mdpi.com/2227-9717/ 9/3/414/s1, Supplementary Material: an HTML print of the jupyter notebook which contains the used code.

**Author Contributions:** Conceptualization, K.K.; methodology, K.K.; software, K.K.; formal analysis, K.K.; investigation, K.K.; data curation, K.K.; writing—original draft preparation, K.K.; writing review and editing, K.K. and R.S.; visualization, K.K.; supervision, R.S.; project administration, K.K. and R.S.; funding acquisition, R.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The Center of Competence for Recycling and Recovery of Waste 4.0 (acronym ReWaste4.0) (contract number 860 884) under the scope of the COMET—Competence Centers for Excellent Technologies—financially supported by BMK, BMDW, and the federal state of Styria, managed by the FFG.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data is contained within the article and Supplementary Materials.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**



**Figure A1.** Cutting tool geometries [24].



\* bottom/top of the teeth.

**Figure A2.** Experimental setup: photo and flow chart [24].

**Figure A3.** Screening drum geometry.

**Table A2.** Data of the screening drums.


### **References**

