The main point of generative-based methods is to learn patterns from a dataset and to model data distributions, allowing the model to be used to generate new samples. Then during training, the model can be updated using the combination of the supervised loss (for existing data with labels) and unsupervised loss (for synthetic data). An inherent advantage of generative-based methods is that the labeled data can be enriched by a trained model which has learned the data distribution. Therefore, generative-based methods are well-suited for situations where it is difficult to collect labeled data, such as process fault detection [153] (Know-what, process level) and anomaly detection in machinery [154] (Know-what, machine level). Examples include the semi-supervised GAN series (SS-GANs), such as Categorical Generative Adversarial Network (CatGAN) [155], Improved GAN [156], and semi-supervised VAEs (SS-VAEs) [157].

Table A3 lists semi-supervised applications in manufacturing taken from the selected documents in Section 2.2.

#### *4.4. Reinforcement Learning Methods*

Reinforcement Learning (RL) algorithms consist of two elements, namely, an *agent* acting within an *environment* (see Figure 15). The agent is acting, and is therefore subject to the desired learning process by directly interacting with and manipulating the environment. Based on [158], the procedure of a learning cycle is as follows: first, the agent is presented with an observation of the environment state *st* ∈ S; then, based on this observation (along with internal decision making), the selection of an action *at* ∈ A. S refers to the state space, that is, the set of possible observations that could occur in the environment. The observation has to provide sufficient information on the current environment or system state in order for the agent to select actions in an ideal way to solve the control problem. For selecting the action, A refers to the action space, that is, the set of possible actions chosen by the agent. After *at* is performed (in a given state *st*), the environment moves to the resulting state *st*+1 and the agent receives a reward *rt*+1. Then, the reinforcement learning cycle continues to iterate as shown in Figure 15. The agent aims to maximize the (discounted) long-term cumulative reward by improving the selection of actions towards an optimum. In other words, the RL agent wants to learn an optimal control policy for the environment.

**Figure 15.** Overview of the Reinforcement Learning approach based on [158].

In general, RL approaches can be split into model-based, i.e., the agent has an internal model of how the environment works, and model-free. The latter is most common thanks to the advent of deep learning, and simplifies application, as feature selection can be applied. Model-free approaches themselves can be divided into short value-based or policy-based approaches by their approach to storing state-action value pairs, which are used to select the action for optimal value return; the latter directly optimize the action selection policy. In contrast to the other machine learning techniques, RL does not require large dataset, only a clearly specified environment. Typically, an RL agent is trained on a simulation or digital twin model [159]; after successful training, it can be implemented on the Know-how level for its original purpose. Otherwise, the agent starts with random non-optimal actions, leading to undesired system behavior.

Considering the aim of achieving the Know-how level for autonomous control in processes, machines, or systems, RL is extremely important for applications in future production. In addition, multi-agent RL is becoming of interest to the research community [33], and can even be applied for controlling products [160]. However, RL remains under-exploited in the industrial area, especially in respect to other machine learning techniques [161].

As of now, applied approaches can be summarized as shown in Table A4. Note that the applications reviewed here are implemented in a simulation or digital twin [159], and features are manually crafted from raw data.

#### **5. Challenges and Future Directions**

A large number of ML use cases have shown the great potential for addressing complex manufacturing problems, from knowing what is happening to knowing how employ selfadapting or self-optimizing systems. The data-driven mechanisms in ML enable broader applications in different fields as well as at different levels, from individual products to whole systems. However, in spite of the great potential and advantages offered by ML and numerous off-the-shelf ML models, there are critical challenges to overcome before the successful application of ML in manufacturing can be realized. The following demonstrate typical challenges that manufacturing industries might confront during the application and deployment of ML-based solutions, along with corresponding future directions for tackling these challenges (RQ3).

• *Lack of data*. Preparing the data used for ML is not a simple task, as the scale and the quality of data can greatly affect the performance of ML models. The most common challenge involves preparing a large amount of organized input data, and ensuring high-quality labels if labels are needed. Despite manufacturing data becoming increasingly more accessible due to the development of sensors and the Internet of Things, gathering meaningful data is time-consuming and costly in many cases, for example, fault detection and RUL prediction. This issue might be alleviated by the Synthetic Minority Over-sampling Technique (SMOTE) [162]. However, SMOTE cannot capture complex representative data, as it often relies on interpolation [163]. Data augmentation [139,164] or transfer learning [165] may address this problem. The aim of data augmentation is to enlarge dataset by means of transforming data [139], by transforming both data and labels, as with MixUp [166], or by generating synthetic data using generative models [167,168]. On the contrary, instead of focusing on expanding

data, transfer learning aims to leverage knowledge from similar external datasets. A typically used method in transfer learning is parameter transfer, where a pretrained model from a similar dataset is employed for initialization [165]. Another situation involving lack of data is that certain data cannot be shared due to data privacy and security issues. In confronting this problem, Federated Learning (FL) [169] might be a potential opportunity to enable model training across multiple decentralized devices while holding local data privately.

To summarize, while ML is a fairly open tool which can be used to handle a variety of problems in manufacturing, it is necessary to have an understanding of the hidden challenges in ML application in order to provide more realistic and robust outcomes. For instance, early in ML application in manufacturing, one might face the problem of lacking data. During the deployment of ML-based solutions, one might confront challenges around integrating the solution into the industrial environment. After deployment, one might encounter the challenge of evaluating ML results on product and process in terms of interpretability and uncertainty. The future directions pointed out in this review can help to address the above-mentioned challenges and ensure reliable improvements in manufacturing contexts.

#### **6. Conclusions**

It is fully recognized that ML is playing an increasingly critical role in the digitization of manufacturing industries towards Industry 4.0, leading to improved quality, productivity, and efficiency. This review has paper aimed to address the issue of how ML can improve manufacturing, posing three research questions related to the above issue in the introduction. To address these research questions, we carried out a literature review assessing the state-of-the-art based on 1348 published scientific articles.

To answer RQ1, we first introduced the concepts of the 'Four-Know' (Know-what, Know-why, Know-when, Know-how) and 'Four-Level' (Product, Process, Machine, System) categories to help formulate ML tasks in manufacturing. By mapping ML use cases into the Four-Know and Four-Level matrix, we provide an understanding of typical ML use cases and their potential benefits for improving manufacturing. To further support RQ1, the identified ML studies were classified using the 'Four-Know' and 'Four-Level' perspective to provide an overview of ML publications in manufacturing. The results showed that current ML applications are mainly focused on the product level, in particular in terms of Know-what and Know-when. In addition, considerable growth in Know-how was observed at the process and system levels, which might be correlated to higher demand for adaptability to changes on these levels.

To fill the gap between academic research and manufacturing industries, we provided an actionable pipeline for the implementation of ML solutions by production engineers from ideation through to deployment, thereby answering RQ2. To further explain the 'model training' step, which is the core stage in the pipeline, a holistic review of ML methods was provided, including supervised, semi-supervised, unsupervised, and reinforcement learning methods along with their typical applications in manufacturing. We hope that this can provide support in method selection for decision-makers considering ML solutions.

Finally, to answer RQ3, we uncovered the current challenges that manufacturing industry is likely to encounter during application and deployment, and provided possible future directions for tackling these challenges as possible developments for ensuring more reliable and robust outcomes in manufacturing.

**Author Contributions:** Conceptualization, T.C., O.J.J., M.C.M., V.S. and G.F.; methodology, T.C. and S.S.; formal analysis, T.C. and S.S.; writing—original draft preparation, T.C., V.S., S.S., M.C.M., O.J.J. and F.S.; writing—review and editing, T.C., V.S., M.C.M., S.S., O.J.J., M.C., G.F., G.T., J.J.A.M. and F.S.; supervision, M.C., G.F., G.T., J.J.A.M. and F.S.; funding acquisition, G.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by a European Training Network supported by Horizon 2020, grant number 814225.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research work was undertaken in the context of the DIGIMAN4.0 project ("DIGItal MANufacturing Technologies for Zero-defect Industry 4.0 Production", https://www. digiman4-0.mek.dtu.dk/, accessed on 1 January 2023). DIGIMAN4.0 is a European Training Network supported by Horizon 2020, the EU Framework Programme for Research and Innovation (Project ID: 814225).

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

**Appendix A**

**Table A1.** Categories of supervised learning applications.

*Appl. Sci.* **2023**, *13*, 1903

**Table A2.** Categories of unsupervised learning applications.

**Table A3.** Categories of semi-supervised learning applications.

**Table A4.** Categories of reinforcement learning applications.

## **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Application of Machine Learning for Prediction and Process Optimization—Case Study of Blush Defect in Plastic Injection Molding**

**Alireza Mollaei Ardestani 1 , Ghasem Azamirad 2,\*, Yasin Shokrollahi 3 , Matteo Calaon 1 , Jesper Henri Hattel 1 , Murat Kulahci 4 , Roya Soltani 5 and Guido Tosello 1**

	- 2800 Kgs. Lyngby, Denmark

**Abstract:** Injection molding is one of the most important processes for the mass production of plastic parts. In recent years, many researchers have focused on predicting the occurrence and intensity of defects in injected molded parts, as well as the optimization of process parameters to avoid such defects. One of the most frequent defects of manufactured parts is blush, which usually occurs around the gate location. In this study, to identify the effective parameters on blush formation, eight design parameters with effect probability on the influence of this defect have been investigated. Using a combination of design of experiments (DOE), finite element analysis (FEA), and ANOVA, the most significant parameters have been identified (runner diameter, holding pressure, flow rate, and melt temperature). Furthermore, to provide an efficient predictive model, machine learning methods such as basic artificial neural networks, their combination with genetic algorithms, and particle swarm optimization have been applied and their performance analyzed. It was found that the basic artificial neural network (ANN), with an average accuracy error of 1.3%, provides the closest predictions to the FEA results. Additionally, the process parameters were optimized using ANOVA and a genetic algorithm, which resulted in a significant reduction in the blush defect area.

**Keywords:** plastic injection molding; design of experiments; machine learning; digital twin; process optimization

**Citation:** Mollaei Ardestani, A.; Azamirad, G.; Shokrollahi, Y.; Calaon, M.; Hattel, J.H.; Kulahci, M.; Soltani, R.; Tosello, G. Application of Machine Learning for Prediction and Process Optimization—Case Study of Blush Defect in Plastic Injection Molding. *Appl. Sci.* **2023**, *13*, 2617. https://doi.org/10.3390/app13042617

Academic Editor: Joamin Gonzalez-Gutierrez

Received: 20 January 2023 Revised: 13 February 2023 Accepted: 14 February 2023 Published: 17 February 2023

**Copyright:** © 2023 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/).

### **1. Introduction**

With the widespread use of algorithms in the engineering sciences, lots of additional costs associated with time-consuming and expensive tests have been eliminated from the product design and production development cycle. Today's modeling, prediction, and optimization methods have extremely reduced the need to use traditional experimental trials and measurements for product and process improvement. These techniques include statistical methods such as Analysis of Variances (ANOVA), machine learning methods such as Artificial Neural Networks (ANNs), as well as optimization methods using meta-heuristic algorithms. Nowadays, the use of Finite Element Analysis (FEA) methods combined with modern optimization methods has helped manufacturers select the optimal levels of input parameters and achieve the highest quality products [1,2]. Due to the complex behavior of polymers, particularly when processed by injection molding, many parameters can affect the product quality. Hence, monitoring and controlling each of these parameters as well as their interactions is vital to prevent injection defects. For example, a fundamental quality

of injection molded parts is based on the reduction of warpage. Among statistical methods, Genetic Algorithm (GA) is known as a popular technique for reducing warpage [3], as well as other defects. Additionally, ANNs have been widely used to predict the shrinkage values of the examined parts [4]. Wang et al. [5], by taking an ANN approach, found that the shrinkage of the manufactured part has an inverse relationship with process parameters such as packing time, injection pressure, holding pressure, and melt temperature, and a direct relationship with cooling time and mold temperature. Altan et al. [6] conducted a study to reduce the shrinkage of a manufactured part. They studied parts composed of polypropylene and polystyrene. After the design of experiments based on the Taguchi method, an ANOVA was performed. As a result, optimal levels for design parameters were obtained, and by employing these optimal parameters in the simulation, the shrinkage of the parts was reduced significantly. An ANN was also trained to predict response values, and it was found that this kind of network has a good ability for this purpose. In a study conducted by Chen et al. [7], it was found that the parameters of holding pressure and melt temperature are the most effective factors in the warpage of the parts. A model was also presented that demonstrated the ability to predict parts warpage in both simulations and experimental tests with relatively high accuracy. It is also possible to use several optimization methods, such as Taguchi method, ANN, GA, Response Surface Methodology (RSM), Particle Swarm Optimization (PSO), etc. simultaneously to optimize the parameters involved in the part's warpage. Finally, according to the prediction error and the number of iterations performed, the best optimization method can be selected [8–11]. The Taguchi optimization method has been used to reduce the warpage of the parts and improve their strength. Consequently, the holding pressure was found to be the most significant factor in the warpage of the parts [12]. Similar research has been conducted using the Taguchi, Nondominated Sorting Genetic Algorithm (NSGA-II), and RSM methods to optimize warpage, shrinkage, and residual stress by Li et al. [13]. Reduction of warpage by considering mold filling as a constraint using the RBF (Radial Basis Function) method was achieved in [14]. The RBF method was also used to perform a multi-objective optimization to reduce production time, clamping force, and warpage in the study by Kitayama et al. [15]. ANNs and GA were used to study shrinkage. Researchers found a relationship between design variables and the target parameter, which showed a great ability of ANN to predict the results [16]. In two separate studies, Kitayama et al. arranged two multi-objective optimizations that both used the RBF method to improve the weld lines. In one of these studies, the weld lines were optimized while reducing the clamping force. In a second study, the goal of the research was to reduce cycle time and weld lines [17,18]. Other than the determination of the optimal levels for the factors related to the process parameters, meta-heuristic algorithms have been also used to determine the optimal values for the composition ratio of polymers to reduce the maximum shrinkage in the molded products [19]. Xu et al. used the Taguchi method, ANN, multi-objective PSO algorithm, and a combination of ANN with PSO to attain the optimal level for design parameters to reduce the part weight, shrinkage, and flash size [20]. Along with the shrinkage of the produced part, other parameters such as clamping force have been inspected simultaneously in a multi-objective optimization using the RBF method [21]. To reduce the part shrinkage, the cycle time, and injection time at the same time, a multi-objective optimization using an ANN-based program was carried out [22]. Some studies have examined other types of defects. For example, the research by Tabi et al. [23] aimed at improving the needle-shaped defects around the gate location. The role of residual stress in the occurrence of cracks in parts [24], and the use of RSM for cutting down the average shear rate around the gate [25] were also investigated. Furthermore, improving the optical properties of polymer optics as a result of reducing residual stress in injection molded plastic lenses was the subject of another study [26]. Li et al. employed Kriging and NSGA-II methods to investigate the effect of runner diameter and process parameters on the quality, cost, and production efficiency [27]. The Taguchi method has been used to investigate the effect of process parameters on the mechanical properties of parts produced with recycled plastic. It was revealed that injection time with

40.5% and material temperature with 43.3% of impact are the most important factors in warpage and yield stress of the produced part, respectively [28]. Additionally, Martowibowo et al. [29] used a GA, with a prediction error of approximately 1%, and could find the optimal values of process parameters in such a way that the production time of the part was significantly reduced. In other studies, FEM and conventional optimization methods were used to decrease the cycle time to a minimum by providing an optimal design for the cooling system [15,30]. Eladl et al. [31] studied the effect of process parameters on the formation of flash defect. During the study, polypropylene (PP) and acrylonitrile butadiene styrene (ABS) were investigated. The outputs of the study were part mass, flow length, and flash formation. After employing DOE and statistical analysis, it was shown the injection speed and packing pressure were the most influential factors on the flash area for both materials. Regi et al. [32] used a different technique to investigate the flow propagation in the molds. They used a transparent window on one of the walls of the mold to visually observe the flow of the material. The objective of the study was to compare FEA simulation results with experiments with focus on flow hesitation. A high-speed camera was used to record the mold filling phase with two different materials of PP and ABS. After employing DOE and ANOVA, the results demonstrated that flow progression and hesitation are dependent on wall thickness, injection velocity, and material type. Loaldi et al. [33], using the same techniques (DOE and FEA) focused on the experimental validation of the injection molding simulation of microparts, also with a focus on flash formation. Results confirmed that higher values of holding pressure, injection speed, mold temperature, melt temperature generate larger flash areas. The trends were correctly predicted by the FEA flow simulations. Chen et al. [34] employed the Taguchi, ANOVA, backpropagation ANN, GA, and Davidson-Fletcher-Powell methods to improve the product quality. As a result of this optimization, as well as increasing demands in terms of products quality, issues such as waste, number of defective parts, need for inspection during production, need for recycling, and production time were also decreased. By using machine learning methods, Finkeldey et al. [35] were able to accurately predict the weight and thickness of the manufactured part. Mehat et al. [36] used process parameters optimization instead of adding additives to improve the mechanical properties in a part produced with recycled plastic. Clearly, considerable research efforts have been focused recently to optimize injection molding parts characteristics and minimize several types of defects. However, very limited research has investigated the factors affecting the blush defect. At least, one study can be indicated to have addressed the influence of process parameters on blush: results by Llado et al. [37] indicated that injection flow rate and melt temperature are the most effective parameters. Blush defect not only deteriorates the appearance of the part but also reduces the lifetime and strength of the product. Therefore, due to its importance and the limited research performed so far, there is a clear need for further investigation on the parameters affecting the incidence and exacerbation of this defect. Despite the importance of preventing the incidence of this defect in plastic injection, so far, only a few studies have investigated its causes and impacts.

The present research has inspected the effects of eight injection molding factors (flow rate, melt temperature, holding pressure, mold temperature, runner diameter, gate diameter, gate angle, included angle) on the blush. The novelty of the present study is the relatively number of investigated parameters, including both process and design factors. By examining a large variety of parameters, the study results in a comprehensive view of the causes of this defect. What is more, contrary to previous research, the interaction effects of design parameters on the incidence of the defect have also been studied in the current study. In addition, for the first time, different types of ANNs have been used to predict the values of blush defect, and a GA is utilized to optimize the levels of effective factors to achieve the lowest probable blush defect.

The research method of this study began with the creation of a CAD model of the injection molded component (plastic bushing). Then, as screening step, the fractional factorial DOE with two levels and the eight design parameters has been conducted. After

performing FEA, the results were included in an ANOVA routine to identify effective and ineffective factors. Thereafter, for a more detailed study of the effect of the most relevant parameters, a second design of experiment of CCD (Central Composite Design) type has been performed with five levels on the effective factors. parameters, a second design of experiment of CCD (Central Composite Design) type has been performed with five levels on the effective factors. **2. Materials and Methods**

The research method of this study began with the creation of a CAD model of the injection molded component (plastic bushing). Then, as screening step, the fractional factorial DOE with two levels and the eight design parameters has been conducted. After performing FEA, the results were included in an ANOVA routine to identify effective and ineffective factors. Thereafter, for a more detailed study of the effect of the most relevant

#### **2. Materials and Methods** *2.1. Modeling*

#### *2.1. Modeling* Blush is a visual defect that occurs as white halos, usually around the gate location.

Blush is a visual defect that occurs as white halos, usually around the gate location. The part under study was a size bushing produced with an injection mold equipped with two cavities. Bushings are a kind a fitting in piping that can be used to connect two pipes together, change the pipeline and flow direction, derive new pipe branches, and to blind a branch. Polyvinyl chloride (PVC) fittings are commonly used in industry due to their ease of use, durability, and cost effectiveness. The geometry of the bushing and injection mold runners is shown in Figure 1. The gate is of round type, and the length of the sprue is 76 mm with a conic angle of 3◦ . All the geometrical dimensions to reproduce the CAD model are included in Figure 1. The part under study was a size bushing produced with an injection mold equipped with two cavities. Bushings are a kind a fitting in piping that can be used to connect two pipes together, change the pipeline and flow direction, derive new pipe branches, and to blind a branch. Polyvinyl chloride (PVC) fittings are commonly used in industry due to their ease of use, durability, and cost effectiveness. The geometry of the bushing and injection mold runners is shown in Figure 1. The gate is of round type, and the length of the sprue is 76 mm with a conic angle of 3°. All the geometrical dimensions to reproduce the CAD model are included in Figure 1.

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 4 of 22

**Figure 1.** Bushing geometry and dimensions of runner system of the original mold. **Figure 1.** Bushing geometry and dimensions of runner system of the original mold.

The material used in the study is a commercially available grade of polyvinyl chloride (PVC) produced by Solvay ET CIE (Brussels, Belgium) under the commercial name Benvic IR705. This type of PVC is widely used in the production of pipes and fittings. Some of its mechanical, thermal, and rheological properties are presented in Table 1 [34]. Figure 2 shows the PVT data and viscosity characteristics (pressure, specific volume, and The material used in the study is a commercially available grade of polyvinyl chloride (PVC) produced by Solvay ET CIE (Brussels, Belgium) under the commercial name Benvic IR705. This type of PVC is widely used in the production of pipes and fittings. Some of its mechanical, thermal, and rheological properties are presented in Table 1 [34]. Figure 2 shows the PVT data and viscosity characteristics (pressure, specific volume, and temperature) of the material based on the Autodesk Moldflow software database.

temperature) of the material based on the Autodesk Moldflow software database. **Table 1.** Properties of Benvic IR705 [34].

density 1.3253 kg/dm3 Shrinkage 0.60%

**Figure 2.** (**a**) pvT data and (**b**) viscosity curves of the Benvic IR705 PVC material. **Figure 2.** (**a**) pvT data and (**b**) viscosity curves of the Benvic IR705 PVC material.

As the first step, the CAD model of the bushing was created. This model was imported into the Autodesk Moldflow software for mesh generation and creating the runner system. Concurrently, the bushing has been masked using 459,096 3D tetragonal elements for both mold cavities. The overall element size of the bushing is 3.5 mm, and due to the As the first step, the CAD model of the bushing was created. This model was imported into the Autodesk Moldflow software for mesh generation and creating the runner system. Concurrently, the bushing has been masked using 459,096 3D tetragonal elements for both mold cavities. The overall element size of the bushing is 3.5 mm, and due to the greater sensitivity of the gate location area, the element size of this area is set to 2 mm.

greater sensitivity of the gate location area, the element size of this area is set to 2 mm. his section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental His section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

#### conclusions that can be drawn. *2.2. Measuring the Simulation Results*

*2.2. Measuring the Simulation Results*  A visualization of blush simulation is shown in Figure 3a. As an assumption, the shape of the defect was assumed, due to the visual similarity and to avoid adding unnecessary complexity to the problem, assumed to be similar to an ellipse. Therefore, the for-A visualization of blush simulation is shown in Figure 3a. As an assumption, the shape of the defect was assumed, due to the visual similarity and to avoid adding unnecessary complexity to the problem, assumed to be similar to an ellipse. Therefore, the formula for calculating the area of the ellipse has been applied to determine the area of the defect. Figure 3b represents the measures needed to calculate of the defect area.

defect. Figure 3b represents the measures needed to calculate of the defect area.

mula for calculating the area of the ellipse has been applied to determine the area of the

**Figure 3.** Width of defect on a simulated bushing (**a**) and geometrical measures for calculation of L (**b**). **Figure 3.** Width of defect on a simulated bushing (**a**) and geometrical measures for calculation of L (**b**).

The area of the ellipse can be calculated by Equation (1): The area of the ellipse can be calculated by Equation (1):

$$\text{Area} = \pi \times \frac{\text{L}}{2} \Big/ 2 \times \frac{\text{W}}{2} \Big/ 2 \tag{1}$$

where *L* and *W* represent the large and the small diameter of the ellipse, respectively. To calculate the area of the defect, a direct measurement of the value of *W*, the value of *L* was calculated by measuring the length of the arc from top edge to bottom edge of the defect, as shown in Figure 3b, via Equations (2)–(7). where *L* and *W* represent the large and the small diameter of the ellipse, respectively. To calculate the area of the defect, a direct measurement of the value of *W*, the value of *L* was calculated by measuring the length of the arc from top edge to bottom edge of the defect, as shown in Figure 3b, via Equations (2)–(7).

$$\beta\_1 = \tan^{-1}(\sqrt[h\_1]{\mathcal{G}\_1})\tag{2}$$

$$\mathfrak{a}\_1 = 180 - (2 \times \mathfrak{f}\_1) \tag{3}$$

$$\mathfrak{i}\_1 = \mathfrak{i}\_2 \mathfrak{i}\_3 \quad \text{and} \quad \mathfrak{i}\_2 = \mathfrak{i}\_3 \mathfrak{j}\_1 \quad \text{and} \quad \mathfrak{i}\_3 = \mathfrak{i}\_1 \mathfrak{j}\_2 \quad \text{and} \quad \mathfrak{j}\_3 = \mathfrak{i}\_2 \mathfrak{j}\_3$$

1 1 Having Equations (2) and (3), the angle of the top arc of defect (*α*1) can be calculated. Additionally, Equations (3) and (4) calculate the same parameter (α2) for the bottom part Having Equations (2) and (3), the angle of the top arc of defect (*α*1) can be calculated. Additionally, Equations (3) and (4) calculate the same parameter (α2) for the bottom part of the defect.

$$\beta\_2 = \tan^{-1}(\sqrt[h\_2]{\xi\_2})\tag{4}$$

$$\mathfrak{a}\_2 = 180 - (2 \times \mathfrak{f}\_2) \tag{5}$$

2 2 αβ = −× 180 (2 ) (5) Equation (6) sums up the amounts of top and bottom angle of the defect's arc. Thus, *αT* represents the total angle of the defect's arc.

$$
\mathfrak{a}\_T = \mathfrak{a}\_1 + \mathfrak{a}\_2 \tag{6}
$$

(7)

$$L = \frac{a\_T}{360} \times 2\pi R\tag{7}$$

 *T* = +1 2 (6) According to Equation (7), the value of *L* can be calculated. Having the values of *L* and *W*, the area of the defect was calculated in each case through Equation (1).

 αα

360

= ×

α

2

π

α

#### *T L R 2.3. Experimental Results and Measurements*

of the defect.

According to Equation (7), the value of *L* can be calculated. Having the values of *L* and *W*, the area of the defect was calculated in each case through Equation (1). *2.3. Experimental Results and Measurements*  To validate the simulations, experimental tests were carried out according to the procedure shown in Figure 4. In the first step, the PVC compound with the properties given in Table 1 was fed into the injection molding machine. Next, process parameters were set To validate the simulations, experimental tests were carried out according to the procedure shown in Figure 4. In the first step, the PVC compound with the properties given in Table 1 was fed into the injection molding machine. Next, process parameters were set in the digital process setting panel. After that, the process started with closing the mold, and after injection, the bushing was manufactured. As the next step, accurate and high-resolution imaging was employed to obtain a high-quality representation of the bushing. Eventually, the images were processed by a computer software to better visualize the blush defect. For experimental measurement of the blush defect area, the length and

in the digital process setting panel. After that, the process started with closing the mold, and after injection, the bushing was manufactured. As the next step, accurate and high-

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 7 of 22

the width of the defect were measured. These measurements enabled us to calculate the defect's area using the ellipse area measurement formula (Equation (1)). blush defect. For experimental measurement of the blush defect area, the length and the width of the defect were measured. These measurements enabled us to calculate the defect's area using the ellipse area measurement formula (Equation (1)).

Eventually, the images were processed by a computer software to better visualize the

**Figure 4.** Experimental test procedure. **Figure 4.** Experimental test procedure.

#### *2.4. Function Approximation by ANOVA 2.4. Function Approximation by ANOVA*

Using Minitab software, a fractional factorial design was first created as the screening step, based on the eight previously mentioned parameters with a ratio of an 1/8. This ratio means that the number of experiments performed in the fractional DOE is 1/8 of the full factorial design. The number of full factorial DOE trials is 256, but with an 1/8 ratio, the fractional DOE is reduced to 32 experiments. The input parameters have been considered to all vary between two levels (see Table 2). After identifying the most effective parameters, a CCD design of the experiments will be conducted to clarify the impact of each effective parameter. After that, the results of ANOVA will represent the effective parameters and eventually a regression predictive model presented by the software. Using Minitab software, a fractional factorial design was first created as the screening step, based on the eight previously mentioned parameters with a ratio of an 1/8. This ratio means that the number of experiments performed in the fractional DOE is 1/8 of the full factorial design. The number of full factorial DOE trials is 256, but with an 1/8 ratio, the fractional DOE is reduced to 32 experiments. The input parameters have been considered to all vary between two levels (see Table 2). After identifying the most effective parameters, a CCD design of the experiments will be conducted to clarify the impact of each effective parameter. After that, the results of ANOVA will represent the effective parameters and eventually a regression predictive model presented by the software.

**Table 2.** Levels of parameters required for performing simulation analyses. **Table 2.** Levels of parameters required for performing simulation analyses.

Higher bound 25 195 35 80 4.5 10 45 45 *2.5. Basic and Hybrid Machine Learning Algorithms*

*2.5. Basic and Hybrid Machine Learning Algorithms*  2.5.1. ANN

next layer [7].

2.5.1. ANN To provide a model for predicting the area of blush, several machine learning algorithms were considered and employed. The algorithms used in this research include basic ANN, combination of ANN and PSO, and combination of ANN and GA. The most straightforward ANN involves an input layer, a hidden layer, and an output layer. Each layer can contain one or more neurons. According to Equation (8), the value of each neu-To provide a model for predicting the area of blush, several machine learning algorithms were considered and employed. The algorithms used in this research include basic ANN, combination of ANN and PSO, and combination of ANN and GA. The most straightforward ANN involves an input layer, a hidden layer, and an output layer. Each layer can contain one or more neurons. According to Equation (8), the value of each neuron is multiplied by the specified weight assigned to each link and added to the same value for all the other neurons in the same layer, and eventually enters the neurons of the next layer [7].

ron is multiplied by the specified weight assigned to each link and added to the same

$$\mathbf{S}\_{\bar{l}} = \sum\_{j=0}^{N} \mathbf{x}\_{l} \times w\_{l\bar{j}} \tag{8}$$

where *xi* is the output of the *ith* neuron of the previous layer, and *wij* is the weight of the link between the *ith* neuron of the previous layer and the *jth* neuron of the present layer. *Sj* represents the sum of the previous layer's outputs multiplied by the connection weights, which is the net input entering the *j th* neuron, and *N* represents the number of inputs to the *jth* neuron in the hidden layer. Each neuron in the network produces its output (*Oj* ) by entering *Sj* that is a Tansig activator function similar to one indicated in Equation (9) [9]:

$$O\_{\dot{j}} = F(S\_{\dot{j}}) = \frac{1 - e^{-S\_{\dot{j}}}}{1 + e^{-S\_{\dot{j}}}} \tag{9}$$

2.5.2. Training ANN

• Basic ANN

In this research, the training was conducted with 70% of the available data provided in the DOE phase. Then, 15% of the available data was used as a validation dataset and the rest, 15% as a test dataset. In the first step, the weights were imported to the network randomly. Applying the gradient descent method causes these weights to be continuously updated during successive iterations. To train the network, the mean square error (*MSE*) of the predicted values should move towards minimization. Equations (10)–(12) represent this method [9,20]:

$$E = MSE = \bigvee\_{N=1}^{N} \sum\_{i=1}^{N} (y\_p - y\_t)^2 \tag{10}$$

$$
\Delta w\_{ij} = -\eta \frac{\partial E}{\partial w\_{ij}} \times O\_j \tag{11}
$$

$$w\_{i\dot{j}}^{m} = w\_{i\dot{j}}^{m-1} + \Delta w\_{i\dot{j}} \tag{12}$$

where *η* is the learning rate and controls the network convergence (a number between 0 and 1), *E* represents the *MSE*, *N* for the number of inputs, *yt* desired output, predicted output by the ANN, and *m* indicates the iterations counter.

• ANN + GA

This method is similar to the basic ANN method, except that instead of using the gradient descend function used in basic ANN, the appropriateness of the weights assigned to each link is determined through the GA process, which is the selection of the best, the crossover, and the mutation, so that the MSE is minimized. In this algorithm, the probability of selecting each parent is assessed through Equations (13) and (14) [9]:

$$f\_i = \begin{cases} \lambda \\ m\_i \end{cases} \tag{13}$$

$$p\_{\bar{i}} = \bigwedge\_{j=1}^{f\_{\bar{i}}} \bigwedge\_{j=1}^{N} f\_{\bar{i}} \tag{14}$$

where *k* is a coefficient, *mi* symbolizes the fitness of *i*th input, *N* substitutes for the number of generation population, and *j* indicates the number of generations. Equation (15) is applied to combine the two chromosomes of *Ci* and *Ch* to produce the next generation [9].

$$\begin{cases} \mathsf{C}\_{ij} = \mathsf{C}\_{ij\left(1-b\right)} + \mathsf{C}\_{h\dot{\mathit{I}}\dot{\mathit{b}}}\\ \mathsf{C}\_{h\dot{\mathit{J}}} = \mathsf{C}\_{h\dot{\mathit{J}}\left(1-b\right)} + \mathsf{C}\_{i\dot{\mathit{J}}\dot{\mathit{b}}} \end{cases} \tag{15}$$

where *b* substitutes for a random number between 0 and 1 and indicates the chromosome's intersection point. In addition, the Equations (16) and (17) are used to model the mutation in this algorithm [9].

$$\mathbf{C}\_{mn} = \begin{cases} \mathbf{C}\_{mn} + (\mathbf{C}\_{mn} - \mathbf{C}\_{\text{max}}) \times f(\mathbf{g}), r > 0.5\\ \mathbf{C}\_{mn} + (\mathbf{C}\_{\text{min}} - \mathbf{C}\_{mn}) \times f(\mathbf{g}), r \le 0.5 \end{cases} \tag{16}$$

$$f(\mathbf{g}) = r\_2(1 - \bigvee\_{\mathbf{G}\_{\text{max}}}) \tag{17}$$

where *Cmn* represents the gene *Cmn*, *C*min and *C*max stands for the higher and lower bounds for genes, and *r* and *r2* are random numbers while *r* is between 0 and 1, *g* indicates the number of present generation population and *Gmax* is the maximum generations considered to iterate. So, the algorithm finally reaches a generation with perfectly fitting responses. A flowchart of the ANN + GA is shown in Figure 5 and its pseudocode is reported in Appendix A. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 10 of 22

**Figure 5.** Training flowchart of: (a): basic ANN, (b): ANN + GA, (c): ANN + PSO. **Figure 5.** Training flowchart of: (a): basic ANN, (b): ANN + GA, (c): ANN + PSO.

In Figure 5, the (a) route shows the basic ANN's flowchart, the (b) route shows the • ANN + PSO

steps of ANN + GA, and the (c) represents the flowchart of ANN + PSO. To calculate the training accuracy of the networks, the Equations (21)–(23) are employed. *P M E M B B P B* − = (20) 1 *n n E E P M n* = (21) (1 ) 100 *Accuracy M Tr E* =− × (22) The PSO algorithm is a method established on swarm intelligence inspired by the behavior of a bird flock. The behavior of the particles, their speed, and the direction of their movement are influenced by the best experience of the particle itself in matching the goal and the best experience among all members of the population. The objective of this algorithm is to reduce the average fitness of all members of the population. In Equations (18)–(20), the best personal position of the particle is presented as *Pbest*, and the best position of the whole particle swarm is symbolized by *Gbest*. If the particles are located in an N-dimensional space, the vector *Xi* = (*Xi*1, *Xi*2, . . . , *XiN*) and the vector *Vi* = (*Vi*1, *Vi*2, . . . , *ViN*) represent the position and velocity of the particle *i*, respectively. In each generation, the values of these vectors must be updated and compared to the previous generation. The Equations (19)–(21)are used to update these values.

and the most accurate network has been employed as the fitness function of a GA to opti-

been shown. In this figure, the simulation and DOE sections indicate all the steps taken through the creation of digital twin, getting the results of it, and employing statistical analysis on the data. The ANN section demonstrates the steps of using different machine learning methods to analyze the data. Then, as the last step, the GA optimization section reveals the work that has been conducted to optimize the parameters' levels to reach the

$$w(t) = w\_{\text{max}} - \frac{t(w\_{\text{max}} - w\_{\text{min}})}{t\_{\text{max}}} \tag{18}$$

$$\boldsymbol{v}\_{i,n}^{t+1} = \boldsymbol{w}\boldsymbol{v}\_{i,n}^{t} + c\_1\boldsymbol{r}\_{1,n}(\boldsymbol{Pbest\_{i,n}^{t}} - \mathbf{x}\_{i,n}^{t}) + c\_2\boldsymbol{r}\_{2,n}(\boldsymbol{Gbest\_n^{t}} - \mathbf{x}\_{i,n}^{t}) \tag{19}$$

$$\mathbf{x}\_{i,n}^{t+1} = \mathbf{x}\_{i,n}^t + v\_{i,n}^{t+1} \tag{20}$$

robust process setting with the lowest blush defect area.

where *t* is counter index for generations, *c*1 and *c*2 stand for two positive coefficients for acceleration while *r*1,*n* and *r*2,*n* are two random coefficients with uniform distribution in *Nth* dimension of the space, *n* shows particle number, and *m* represents the inertia weight employed [7]. Figure 5 depicts a flowchart of steps of optimizing layer weights and training the network using basic ANN, ANN + GA, and ANN + PSO. The pseudocodes for each of the three algorithms are also given in the Appendix A.

In Figure 5, the (a) route shows the basic ANN's flowchart, the (b) route shows the steps of ANN + GA, and the (c) represents the flowchart of ANN + PSO. To calculate the training accuracy of the networks, the Equations (21)–(23) are employed.

$$P\_E = \frac{B\_P - B\_M}{B\_M} \tag{21}$$

$$M\_E = \frac{\sum\_{1}^{n} \left| P\_{E\_n} \right|}{n} \tag{22}$$

*AccuracyTr* = (1 − *ME*) × 100 (23)

where *PE* indicates the prediction error for each data group, *BP* and *BM* represent the predicted value of the blush defect and the measured value of the defect, respectively. *n* indicates the number of data sets, *ME* represents the average error of the whole data set, and *AccuracyTr* shows the accuracy of training.

After training the networks, an optimization was conducted. Additionally, after training an accurate ANN, the prediction accuracy of the network has been measured, and the most accurate network has been employed as the fitness function of a GA to optimize the levels of design parameters. Then, the levels optimized by ANOVA and the GA will be compared together. In the flowchart shown in Figure 6, the research method has been shown. In this figure, the simulation and DOE sections indicate all the steps taken through the creation of digital twin, getting the results of it, and employing statistical analysis on the data. The ANN section demonstrates the steps of using different machine learning methods to analyze the data. Then, as the last step, the GA optimization section reveals the work that has been conducted to optimize the parameters' levels to reach the robust process setting with the lowest blush defect area. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 11 of 22

**Figure 6.** Flowchart of the research method. **Figure 6.** Flowchart of the research method.

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

geometry of the mold.

**Inputs**

**Outputs**

**Table 3.** Validation of FEA.

#### *3.1. FEA Validation*

*3.1. FEA Validation*  To ensure the accuracy of the FEA, experimental tests were performed and compared with the simulation results. Table 3 shows the validation results. As can be observed from this table, the results of FEA are in relatively good agreement with the results of the experimental tests, and the absolute value of the average deviation error is 5.6%. The param-To ensure the accuracy of the FEA, experimental tests were performed and compared with the simulation results. Table 3 shows the validation results. As can be observed from this table, the results of FEA are in relatively good agreement with the results of the experimental tests, and the absolute value of the average deviation error is 5.6%. The parameters of flow rate, melt temperature, and holding pressure were varied to verify the

eters of flow rate, melt temperature, and holding pressure were varied to verify the results. Other parameters remained constant due to practical limitations in changing the

 **Experimental Test No. 1 2 3 4**

**Flow rate (cm3/s)** 15 15 25 25

**Melt temperature (°C)** 185 195 185 195 **Mold temperature (°C)** 35 35 35 35 **Holding pressure (MPa)** 60 80 60 80 **Runner diameter (mm)** 10 10 10 10 **Gate diameter (mm)** 3.5 3.5 3.5 3.5 **Gate angle (°)** 0 0 0 0 **Included angle (°)** 45 45 45 45

**Experimental defect area (mm2)** 2108 1621 2483 1696 **Simulation defect area (mm2)** 2212 1553 2777 1674

According to the shear stress heat map, the area that exceeded the maximum allowable shear stress (0.2 MPa) appeared as the blush [34]. As can be concluded from Table 3, lower amounts of error were observed for simulations associated with process configurations, resulting in a smaller area of defect. Since the objective of the study is to decrease the area of defect, the closer the study gets to its goal, the lower the error amount. Figure 7 shows a comparison between the simulation results for maximum shear stress and the

blush defect in an experimental test performed with the same parameter settings.

**Deviation error** 4.9% −4.2% 11.8% −1.3%

results. Other parameters remained constant due to practical limitations in changing the geometry of the mold.

**Table 3.** Validation of FEA.

According to the shear stress heat map, the area that exceeded the maximum allowable shear stress (0.2 MPa) appeared as the blush [34]. As can be concluded from Table 3, lower amounts of error were observed for simulations associated with process configurations, resulting in a smaller area of defect. Since the objective of the study is to decrease the area of defect, the closer the study gets to its goal, the lower the error amount. Figure 7 shows a comparison between the simulation results for maximum shear stress and the blush defect in an experimental test performed with the same parameter settings. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 12 of 22

**Figure 7.** Simulation result (**a**) and enclosed area of defect in the experimental test (**b**) obtained with the same injection molding process settings (melt temperature = 185 °C, mold temperature = 35 °C, flow rate = 25 cm3/s, holding pressure = 80 MPa, runner diameter = 10 mm, gate diameter = 3.5 mm, gate angle = 0°, and included angle = 45°). **Figure 7.** Simulation result (**a**) and enclosed area of defect in the experimental test (**b**) obtained with the same injection molding process settings (melt temperature = 185 ◦C, mold temperature = 35 ◦C, flow rate = 25 cm3/s, holding pressure = 80 MPa, runner diameter = 10 mm, gate diameter = 3.5 mm, gate angle = 0◦ , and included angle = 45◦ ).

#### *3.2. Statistical Analysis 3.2. Statistical Analysis*

**Table 5.** The levels of CCD.

**No. Flow Rate**

**(cm3/s)**

To determine the effective parameters, a screening step was implemented using the fractional factorial design. The ANOVA results are given in Table 4. These results indicate that the parameters of flow rate, melt temperature, holding pressure, and runner diameter were effective due to the lower *p*-value of 0.05. The parameters of mold temperature, gate diameter, gate angle, and included angle only slightly impacted the defect area. Furthermore, the meaning of squares, which is obtained by dividing the treatment sum of squares by the degrees of freedom, can be used to determine which factors are significant. The higher the mean square, the more of an effect it has on the results [38,39]. To determine the effective parameters, a screening step was implemented using the fractional factorial design. The ANOVA results are given in Table 4. These results indicate that the parameters of flow rate, melt temperature, holding pressure, and runner diameter were effective due to the lower *p*-value of 0.05. The parameters of mold temperature, gate diameter, gate angle, and included angle only slightly impacted the defect area. Furthermore, the meaning of squares, which is obtained by dividing the treatment sum of squares by the degrees of freedom, can be used to determine which factors are significant. The higher the mean square, the more of an effect it has on the results [38,39].

**No. Source Mean of Squares** *p***-Value** Model 3,424,612 0.011 Flow rate 4,978,362 0.011 Melt temperature 19,348,900 0.001

Runner diameter 32,488,905 0.000 Gate diameter 458,198 0.250 Gate angle 1946 0.934 Included angle 96,108 0.571

lead to a better understanding of the parameters' effects on the process outputs.

1 12 185.2 54.6 6.4 2 14 187 59 7 3 19 191 70 8.5 4 24 196 81 10 5 26 197.8 85.4 10.6

**Melt Temperature (°C)**

The DOE intended for the screening step to be executed with all parameters varying between two levels. In the second phase, by performing CCD with an alpha coefficient of 1.4, the number of parameter levels was increased to five (shown in Table 5). This could

> **Holding Pressure (MPa)**

**Runner Diameter (mm)**

**Table 4.** ANOVA results of the screening step.

**Table 4.** ANOVA results of the screening step.

The DOE intended for the screening step to be executed with all parameters varying between two levels. In the second phase, by performing CCD with an alpha coefficient of 1.4, the number of parameter levels was increased to five (shown in Table 5). This could lead to a better understanding of the parameters' effects on the process outputs.

**Table 5.** The levels of CCD.

To predict the area of the defect with a regression equation, an ANOVA was implemented on the results of the simulations suggested by the CCD. According to Table 6 and considering the *p*-value of the model equal to 0.000, it can be concluded that the model presented by ANOVA has good accuracy in the prediction of the defect area. The effect of all four parameters known to be effective in the previous step has been reaffirmed by the design of experiments carried out by considering more levels. Diagrams related to effective parameters can be seen in Figure 8.

**Table 6.** Results of the second ANOVA (CCD).

To predict the area of the defect with a regression equation, an ANOVA was implemented on the results of the simulations suggested by the CCD. According to Table 6 and considering the *p*-value of the model equal to 0.000, it can be concluded that the model presented by ANOVA has good accuracy in the prediction of the defect area. The effect of all four parameters known to be effective in the previous step has been reaffirmed by the design of experiments carried out by considering more levels. Diagrams related to effec-

**Squares** *p***-Value**

**No. Source Mean of**

1 Model 1,463,290 0 **2 Linear parameters 4,473,581 0**  3 Flow rate 1,209,310 0.039 4 Melt temperature 1,369,223 0.049 5 Holding pressure 460,154 0.027 6 Runner diameter 16,845,635 0 **7 Squares 533,923 0.008**  8 Flow rate × Flow rate 6013 0.815 9 Melt temperature × Melt temperature 53,785 0.486 10 Holding pressure × Holding pressure 527 0.945 11 Runner diameter × Runner diameter 1,474,471 0.002 **12 Interactions 76,008 0.641**  13 Melt temperature × Flow rate 4789 0.834 14 Melt temperature × Holding pressure 9658 0.767 15 Melt temperature × Runner diameter 50,917 0.498 16 Flow rate × Holding pressure 126 0.973 17 Flow rate × Runner diameter 175,496 0.216 18 Holding pressure × Runner diameter 215,063 0.173

tive parameters can be seen in Figure 8.

**Table 6.** Results of the second ANOVA (CCD).

**Figure 8.** The effect of: (**a**) flow rate, (**b**) melt temperature, (**c**) holding pressure, and (**d**) runner diameter on the defect area. **Figure 8.** The effect of: (**a**) flow rate, (**b**) melt temperature, (**c**) holding pressure, and (**d**) runner diameter on the defect area.

Figure 8 illustrates the manner of the impact of each parameter on the results. Ac-

Figure 8 illustrates the manner of the impact of each parameter on the results. According to the *p*-value obtained from the effects of ANOVA, the runner diameter has the most significant impact on the area of the blush defect. After that, the parameters of holding pressure, flow rate, and melt temperature have the highest impact on the defect area, respectively. Utilizing ANOVA, a regression equation to predict the area of the defect has been worked out using all input parameters presented in Equation (24). cording to the *p*-value obtained from the effects of ANOVA, the runner diameter has the most significant impact on the area of the blush defect. After that, the parameters of holding pressure, flow rate, and melt temperature have the highest impact on the defect area, respectively. Utilizing ANOVA, a regression equation to predict the area of the defect has been worked out using all input parameters presented in Equation (24).

$$\text{Area} = -12367 + (25.7 \times A) - (25.4 \times B) - (13.82 \times C) + (3898 \times D) - (193.2 \times D \times D) \tag{24}$$

where *A* indicates the flow rate, *B* the melt temperature, *C* the holding pressure, and *D* represents the runner diameter. The standardized residual and histogram plot shown in Figure 9 reveal the accuracy of Equation (24). As can be seen from the Versus Order plot, the predicted data are in general equally scattered around the zero line. Additionally, the histogram plot in this figure shows that the amount of data with a standard residual equal to or near zero is much higher than the amount of data with higher residuals. Thus, it can be concluded that the defect area amounts predicted by the regression model are in good agreement with the FEA results. where *A* indicates the flow rate, *B* the melt temperature, *C* the holding pressure, and *D* represents the runner diameter. The standardized residual and histogram plot shown in Figure 9 reveal the accuracy of Equation (24). As can be seen from the Versus Order plot, the predicted data are in general equally scattered around the zero line. Additionally, the histogram plot in this figure shows that the amount of data with a standard residual equal to or near zero is much higher than the amount of data with higher residuals. Thus, it can be concluded that the defect area amounts predicted by the regression model are in good agreement with the FEA results.

After training the ANN, the network should be verified with the data intended for testing. The test data in this study included 15% of the total data that had not been used

**Size Iterations AccuracyTr**

10 245 97.88% 59

4 178 99.96% 2

10 276 98.76% 51

**Training Time (s)**

**Figure 9.** Standardized residual and histogram graph. **Figure 9.** Standardized residual and histogram graph.

*3.3. ANN Validation and Comparison with ANOVA*

**Table 7.** Comparison of response prediction methods.

**Number of Particles**

**No. Predictor Neurons in**

Basic ANN

ANN + PSO

ANN + GA

1

5

9

**Hidden Layer**

**Population**

2 6 134 99.97% 2 3 8 91 99.99% 1 4 10 102 99.98% 1

6 30 198 98.23% 178 7 50 179 98.78% 264 8 70 171 99.26% 408

10 30 251 98.61% 159 11 50 237 98.33% 256 12 70 209 98.17% 357 13 ANOVA 86.57% -

#### *3.3. ANN Validation and Comparison with ANOVA*

After training the ANN, the network should be verified with the data intended for testing. The test data in this study included 15% of the total data that had not been used in the ANN training phase. Table 7 also compares the prediction accuracy of different types of ANN with that presented by ANOVA.

**Table 7.** Comparison of response prediction methods.

According to Table 7, the best networks of each type have been chosen to compare the prediction accuracy. For basic ANN, the trained network in row 3 is considered as the selected basic ANN with the best performance. The type of activator function is "TanSig" for the first layer and "Linear" for the second layer. For ANN + PSO, the trained network in row 7 is regarded as the most suitable algorithm of this type. The trained network in row 11 has also been chosen as the most suitable ANN + GA combination. The mutation rate in this algorithm is 15%, and the cross-over rate is 65%. Additionally, Figure 10 exhibits a comparison between the predicted values via the regression formula provided by ANOVA, selected neural networks of all three types (basic ANN, ANN + GA, and ANN + PSO) with the value obtained from the FEA. Furthermore, Table 8 presents the values predicted by each algorithm for all the data sets compared with ANOVA predicted values and FEA results. According to Table 7, the best networks of each type have been chosen to compare the prediction accuracy. For basic ANN, the trained network in row 3 is considered as the selected basic ANN with the best performance. The type of activator function is ''TanSig'' for the first layer and ''Linear'' for the second layer. For ANN + PSO, the trained network in row 7 is regarded as the most suitable algorithm of this type. The trained network in row 11 has also been chosen as the most suitable ANN + GA combination. The mutation rate in this algorithm is 15%, and the cross-over rate is 65%. Additionally, Figure 10 exhibits a comparison between the predicted values via the regression formula provided by ANOVA, selected neural networks of all three types (basic ANN, ANN + GA, and ANN + PSO) with the value obtained from the FEA. Furthermore, Table 8 presents the values predicted by each algorithm for all the data sets compared with ANOVA predicted values and FEA results.

**Figure 10.** Comparison of normalized neural network and ANOVA prediction values with FEA results. **Figure 10.** Comparison of normalized neural network and ANOVA prediction values with FEA results.

**Table 8.** Comparison of normalized values predicted by ANOVA, ANN Algorithms, and FEA.

2 1 0.849 1 0.688 0.754 3 0.762 0.83 0.762 0.779 0.704 4 −1 −0.219 −1 −0.551 −0.52 5 0.284 0.477 0.284 0.38 0.43 6 0.21 0.417 0.237 0.044 0.016 7 0.015 0.706 0.015 0.876 0.783 8 −1 −0.571 −1 −0.781 −0.786 9 0.859 0.799 0.233 0.811 0.686 10 −1 −0.54 −1 −0.797 −0.745 11 −1 −0.389 −1 −0.71 −0.625 12 −0.376 0.311 −0.376 −0.342 −0.124 13 0.444 0.558 0.444 0.074 0.231 14 −1 −0.37 −1 −0.623 −0.65 15 0.393 0.523 0.393 0.296 0.152 16 0.274 0.629 0.274 0.742 0.531 17 0.945 1 0.945 0.809 0.81 18 −1 −1 −1 −0.886 −0.854 19 0.507 0.679 0.507 0.517 0.644

**Table 8.** Comparison of normalized values predicted by ANOVA, ANN Algorithms, and FEA.

As it is clear from Figure 10 and Table 8, ANOVA with 86.5% accuracy and basic ANN with 99.99% accuracy provide the farthest and the closest predictions to the FEA response, respectively. Additionally, according to these data, the basic ANN presented in row 3 of Table 7 (with eight neurons in the hidden layer) has been chosen as the most appropriate prediction model. The graph presented in Figure 11 indicates the training procedure of the basic ANN. During 91 epochs, the basic ANN converged. The MSE trend of training data can be seen alongside this trend for validation and test data. In addition, this figure represents the regression analysis of the basic ANN model for training, validation, and test data. Additionally, the regression model for the whole data set can be seen in the bottom right corner of Figure 11.

#### *3.4. Optimization Using GA*

To optimize the parameter levels, a trained ANN with the specifications listed in the third row of Table 4 is used as the cost function. The parameters of melt temperature, flow rate, holding pressure, and runner diameter are considered as the inputs of the GA, and the area of the blush defect is regarded as the only output. The allowable range for each of these parameters is specified in Table 3. Additionally, for the GA, values of 300 were taken into account as the initial population, 300 as the maximum number of generations, 60% as the crossover rate, and 40% as the mutation rate. During the optimization process, the value of the cost function of the GA is constantly decreasing, which indicates the proper performance of the algorithm in finding the optimal values. At the end of the process, the algorithm had reached the optimal response after approximately 170 iterations, after which the response showed little change. The algorithm introduces the values of 197.1,

25.6, 84.9, and 6.4 as the optimal levels for the parameters of melt temperature, flow rate, holding pressure, and runner diameter, respectively. Figure 12 shows the optimization process by the GA. Using experimental validation, it is possible to compare the results of the experimental tests, GA, and ANOVA. This comparison can be seen in Table 5. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 17 of 22

**Figure 11.** Regression analysis of the basic ANN model. **Figure 11.** Regression analysis of the basic ANN model. *Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 18 of 22

197.1, 25.6, 84.9, and 6.4 as the optimal levels for the parameters of melt temperature, flow rate, holding pressure, and runner diameter, respectively. Figure 12 shows the optimiza-

From the results obtained from Table 9, it can be seen that, on account of optimization using GA, the amount of blush defect area has been reduced by 81.7%. Additionally, the

> **Flow Rate (cm3/s)**

**Holding Pressure (mm)**

**Runner Diameter (mm)**

**Defect Area (mm2)**

of the experimental tests, GA, and ANOVA. This comparison can be seen in Table 5.

Figure 13 shows a visual comparison between the blush area, before and after optimization in the injection molding simulation, alongside the image of the defect in the bushing

cess, the algorithm had reached the optimal response after approximately 170 iterations, after which the response showed little change. The algorithm introduces the values of **Figure 12.** Optimization process by GA. **Figure 12.** Optimization process by GA.

**Title**

produced with the optimized parameters.

*3.4. Optimization Using GA*

**Table 9.** Comparison of the defect area before and after optimization.

**Optimal bushing suggested by CCD optimization** 197.8 26.0 56.4 6.7 517

**Experimental test on the dataset suggested by GA** 197.1 25.6 84.9 6.4 366

(**a**) (**b**)

(**c**) (**d**)

bushing is after optimization in simulation (**c**) and experimental test (**d**).

**Figure 13.** The initial defected bushing area is in simulation (**a**) and experimental test (**b**); and the

**Optimal bushing suggested by GA** 197.1 25.6 84.9 6.4 362

**Melt Temperature (°C)**

**Initial bushing** 197.0 12.0 35.0 10.0 1978

From the results obtained from Table 9, it can be seen that, on account of optimization using GA, the amount of blush defect area has been reduced by 81.7%. Additionally, the optimization performed using the CCD method has reduced the defect area by 74.0%. Figure 13 shows a visual comparison between the blush area, before and after optimization in the injection molding simulation, alongside the image of the defect in the bushing produced with the optimized parameters. From the results obtained from Table 9, it can be seen that, on account of optimization using GA, the amount of blush defect area has been reduced by 81.7%. Additionally, the optimization performed using the CCD method has reduced the defect area by 74.0%. Figure 13 shows a visual comparison between the blush area, before and after optimization in the injection molding simulation, alongside the image of the defect in the bushing produced with the optimized parameters.

**Table 9.** Comparison of the defect area before and after optimization. **Table 9.** Comparison of the defect area before and after optimization.

**Figure 12.** Optimization process by GA.

*Appl. Sci.* **2023**, *13*, x FOR PEER REVIEW 18 of 22

**Figure 13.** The initial defected bushing area is in simulation (**a**) and experimental test (**b**); and the bushing is after optimization in simulation (**c**) and experimental test (**d**). **Figure 13.** The initial defected bushing area is in simulation (**a**) and experimental test (**b**); and the bushing is after optimization in simulation (**c**) and experimental test (**d**).

#### **4. Conclusions**

This study aimed at the determination of the effect of eight process parameters (flow rate, melt temperature, mold temperature, holding pressure, runner diameter, gate diameter, gate angle, and included angle) on the blush defect in PVC bushings produced by injection molding. Some prediction models have been created to estimate the area of the blush defect using ANOVA and ANNs. Among the prediction models used, the basic ANN method with a training accuracy of 99.99% has shown the best performance compared to other prediction methods for predicting the FEA results.

Results showed that flow rate, melt temperature, and runner diameter have a particularly strong effect on the blush defect and can be related to the viscosity of the molten material. The viscosity can increase with rapid cooling of the material, making the material's shear stress exceed the allowable limit and, in turn, promoting blush.

Holding pressure has also affected the blush defect. With the increase in holding pressure, the blush defect decreases. Lower values of holding pressure can cause some semi-cooled material (which has a high viscosity) to enter the mold cavity. Reflowing the material in these high viscosity conditions can cause high amounts of shear rate, which is the underlying reason of blush defect occurrence. The parameters of mold temperature, gate diameter, gate angle, and included angle have a negligible effect on the result. The key results of the research can be summarized as follows:

**Author Contributions:** Methodology, A.M.A. and G.T.; Software, A.M.A.; Validation, A.M.A., R.S. and G.T.; Formal analysis, R.S.; Investigation, A.M.A.; Resources, G.A.; Writing—original draft, A.M.A.; Writing—review & editing, G.A., Y.S., M.C., J.H.H., M.K. and G.T.; Visualization, A.M.A.; Supervision, G.A. and G.T.; Project administration, G.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research work was partially supported by the DIGIMAN4.0 project ("DIGItal MANufacturing Technologies for Zero-defect Industry 4.0 Production", https://www.digiman4-0.mek.dtu. dk/, accessed on 13 February 2023). DIGIMAN4.0 is a European Training Network supported by Horizon 2020, the EU Framework Programme for Research and Innovation (Project ID: 814225).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This research work was undertaken in the context of DIGIMAN4.0 project ("DIGItal MANufacturing Technologies for Zero-defect Industry 4.0 Production", http://www. digiman4-0.mek.dtu.dk/). DIGIMAN4.0 is a European Training Network supported by Horizon 2020, the EU Framework Programme for Research and Innovation (Project ID: 814225). Yazd Poolica Co. (Yazd, Iran) is thanked for all their support during the experimental steps of the study.

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

#### **Appendix A. Pseudocodes of the Algorithms**

	- 1. Start of program
	- 2. Train percent = tr
	- 3. Test percent = ts
	- 4. Initialize data
	- 5. Separate input and output data
	- 6. Normalize all data
	- 7. Initialize a network structure
	- 8. Set a random matrix for initial layer weights and biases
	- 9. repeat
	- 10. Use tr percent of data for training network
	- 11. Use Equation (10) to Evaluate trained network with ts percent of data
	- 12. Use Equations (11) and (12) to reset the weights and biases
	- 13. Until termination criteria
	- 14. Simulate the trained network with input data
	- 15. Calculate and report the MSE

For the ANN + PSO and ANN + GA algorithms, all steps are similar to the above pseudocode. Only the steps for optimizing layer weights (lines 8–12) are overwritten. In other words, the following pseudocodes are replaced with the iterative loop in the upper pseudocode.

	- 1. Start of weight optimization
	- 2. for each particle(each neural network)
	- 3. Initialize particle position(Network weights) and velocity vectors
	- 4. end for
	- 5. Fitness = f(X)i
	- 6. Personal best position = Xpb
	- 7. Global best position = Xgb
	- 8. repeat
	- 9. for particlei i = 1 to Nparticles
	- 10. Fitness = calculate the Fitness of particlei
	- 11. if f(Xbp)i < f(Xbp)
	- 12. Xbp = Xi, Pbest = Fitnessi
	- 13. end if
	- 14. if f(Xgb)i > f(Xgb)
	- 15. Xgb = Xi, Gbest = Fitnessi
	- 16. end if
	- 17. update the velocity of the particle using Equations (18 and (19)
	- 18. update the position of the particle using Equation (20)
	- 19. end for
	- 20. Until termination criteria
	- 1. Start of weight optimization
	- 2. N = number of network weights
	- 3. G = Number of maximum generations
	- 4. RecomPercent = r/100
	- 5. CrossPercent = c/100
	- 6. MutatPercent = 1 − RecomPercent − CrossPercent
	- 7. Initialize genes randomly(initial Network weights)
	- 8. Calculate fitness = f(X)i
	- 9. Sort Chromosomes according to fitness
	- 10. for i = 1 to G
	- 11. Create new population(RecomPercent × N + CrossPercent × N + MutatPercent × N)
	- 12. Calculate fitness
	- 13. Sort Chromosomes according to fitness
	- 14. end for

## **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Study and Simulation of an Under-Actuated Smart Surface for Material Flow Handling**

**Edoardo Bianchi 1,\* , Oliver Jonas Jorg 2 , Gualtiero Fantoni 2 , Francisco Javier Brosed Dueso 1 and José A. Yagüe-Fabra 1**

**Abstract:** Smart surfaces are becoming more and more popular in the field of intralogistics, as they combine great flexibility with easy reprogrammability. Pursuing this trend, the following article proposes a modular surface to perform handling tasks, such as sorting, stopping, or slowing down material flows. Differently from the current technology, the surface used is under-actuated, thus, it exploits the speed, already possessed by the object, or the gravity to perform, with a simplified hardware, for the aforementioned tasks. In practice, these handling actions are completed using an array of rotors, of which only the direction of the rotation axis is controlled. Moreover, the axis can only assume certain discrete orientations in the plane, further simplifying the design. Thus, what is created is a controllable and under-actuated friction field, which, in contrast with similar existing systems, does not require active driving forces to manipulate the material flow. In the article, the analytic model of the surface is described, and a software simulation environment is introduced to demonstrate its functioning. In addition, examples of sorting, slowing down, and stopping operations and a validation of the simulation itself are presented.

**Keywords:** smart surface; friction force field; under-actuation; feeding; simulation; material flow handling

### **1. Introduction**

The material transport inside factories and warehouses is a very studied topic. However, companies are always searching for new, cheaper solutions while keeping an eye on efficiency and flexibility [1–4]. Starting with the most basic conveyors, many alternative systems [5–12] have been developed in order to move, sort, orient, and handle the material flow. More recently, the focus of investigation has shifted towards modular devices and surfaces that allow high handling capability, adaptability, and reconfigurability [3,4,13]. This happens because the market increasingly demands a flexible industry, and this is reflected in all production levels, even the internal transportation. These devices are often called smart surfaces because they are controllable and reprogrammable with computers. Thanks to these characteristics, these new systems permit achieving their goals without structural changes [4,13]. In addition, the new devices are capable of identifying a part with sensors and acting in diverse ways, according to their programming. As a result, supervision by an operator is no longer required for the decision process and the system becomes autonomous, promising to reduce errors and costs [14].

The most relevant solutions among the smart surfaces, according to the current state of the art and classified by the operating principle, are: micro electro-mechanical systems, vibrating surfaces, ciliary motion, variable morphology surfaces, pneumatic surfaces, surfaces with rotors, and mobile platforms. Micro electro-mechanical systems (MEMS) [6,15,16] are an array of microscopic cantilevers or tilting planes, actuated electrically, that generate

**Citation:** Bianchi, E.; Jorg, O.J.; Fantoni, G.; Brosed Dueso, F.J.; Yagüe-Fabra, J.A. Study and Simulation of an Under-Actuated Smart Surface for Material Flow Handling. *Appl. Sci.* **2023**, *13*, 1937. https://doi.org/10.3390/ app13031937

Academic Editor: Muhammad Junaid Munir

Received: 19 December 2022 Revised: 10 January 2023 Accepted: 31 January 2023 Published: 2 February 2023

**Copyright:** © 2023 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/).

6

6

forces to transport the object in contact. The second category of the list are vibrating surfaces. These consist of vibrating plates on which a sequence of supply frequencies is applied to generate a two-dimensional force field for the handling tasks [5,17]. In [12,18–20], ciliary motion is proposed to move and manipulate objects, taking advantage of an array of controllable cilia for contact conveyance. Another class are variable morphology surfaces [9,21–23]. For these devices, gravity together with inclined planes and vertical actuators at different altitudes are used to create a preferential path for the movement and rotation of an object. Additionally, pneumatic surfaces [7,24–30] were introduced to handle materials without direct contact. The working principle for the latter is to have modules with nozzles to move parts by directing the air flow below them. On the other hand, surfaces with rotors take advantage of contact forces made by actively driven wheels to manipulate the material flow [3,11,13,14,31–35]. Finally, mobile platforms [8,36–39] consist of mobile pallets, not connected to fixed axes and free to move on a plane, transporting objects on top. Despite the categorization, not all of these are practically used for macroscopic intralogistics purposes (object size and displacement > cm). In fact, as summarized in Figure 1, some of them, such as MEMS, together with some mobile platforms [8,37–39] and ciliary devices [12,18,19], are more suitable for microscopic transport (object size and displacement < mm). In contrast, pneumatic systems, vibrating surfaces, surfaces with rotors, variable morphology surfaces, some mobile platforms [36], and even tilted brushes (Cilia) [20] can be used to move macroscopic and mesoscopic (object size and displacement > mm, <cm) objects. 2 3 4 5 D D

1

1

1 / 1 Edizione Foglio

B B **Figure 1.** Smart surfaces classification and field of application.

2 3 4 5 A A Disegno schema tipi sm\_sf ebianchi 14/09/2022 Progettato da Controllato da Approvato da Data Data Referring to the previous classification, the system proposed in this paper can be included in the surfaces with rotors category. The literature review conducted on this category identified two main types of modular smart surfaces: systems where each module has only one degree of freedom [3,11,13,35] and systems where each module has two [14,34]. The former, which are the most recent in the literature, exploit multiple omnidirectional and motorized wheels called "omniwheels" positioned in the module with their axes fixed. The layout and number of wheels is not the same in every study, however; as a reference, three or more wheels are usually employed and their spinning axes are positioned out of alignment to create controllable forces in both directions of the plane. As first example of modules with more than three wheels, study [13] alternates units with seven wheels—one large in the center and six small on the sides, to units with five—two large on the sides and three small in the center. For both layouts, the axes of the large wheels are perpendicular to the axes of the small ones to ensure driving forces in the plane. In contrast, in [3,11,35], each module of the device studied, called "Celluveyor", contains three omnidirectional wheels with axes arranged at 120° to one another. Each wheel can rotate at different speeds to control the magnitude of the force exchanged with the transported object, and all three

forces together are used to handle the body motion. On the other hand, the second group of systems, i.e., those with modules with two degrees of freedom, consist of units composed of one or more motorized and swiveling wheels. In this case, therefore, as presented in [14,34], the basic element of the module is a simple wheel driven by a motor mounted on its axis of rotation. This wheel, however, unlike the one degree of freedom systems, is also mounted on a vertical axis, which is itself motorized and therefore swiveling. Summarizing, these devices just described are fully implemented and contain several motors per module, specifically, three or more for the first class and at least two for the second. They can actively drive, sort, and manipulate the material flow, creating a totally actuated surface with significant handling skills. However, the use of several motors and their control necessarily involves a greater effort in the management of the system, as well as increasing costs and complexity. According to the current state of the art, a simpler under-actuated device without motors is missing for the same intralogistic tasks. For that reason, seeking cost reduction and simplification of both control and design, without losing the flexibility of this class of systems, in this paper the authors propose an under-actuated surface composed of modules. Unlike existing concepts, each module contains an idle rotor instead of a motorized wheel. The axis of rotation of the mentioned rotor is not driven continuously by a motor but can be oriented in a limited number of directions in the plane of the surface, creating a directionable friction force for object handling. In order to compensate for the under-actuation, the material is moved by exploiting gravity or an initial velocity of the object provided by another system (e.g., a previous conveyor belt). Therefore, compared to the existing systems [3,11,13,14,31–35], the novelty that distinguishes this active surface from the literature technology lies in its simplicity. In fact, the two main characteristics of the system proposed by the authors, compared to those in the same category, are the under-actuation and the limited directions of the rotor axis, both of which are reductions to a minimal form of concept technology, but which save components and make the most of sources already present in the application environment, such as gravity and the speed of the objects in the transport line. Against potential assumptions of performance losses, the authors proved through simulations that the same goals and efficiency [3,4,35] of sorting and handling can be achieved with a minimal design architecture, consisting of few components, while decreasing costs and saving energy.

The focus of this article lies on the study of the surface, the description of the working principle, and its simulation with the purpose of providing an initial proof of concept for possible applications and future developments. Furthermore, the simulation environment allowed numerical results to be obtained for typical intralogistics tasks, giving the opportunity for an initial comparison with the corresponding current technology. The simulations presented in this article are carried out with a customized code developed by the authors using the software MATLAB (Version R2022a). In addition, a validation of the simulation results was obtained with the well-known program for multi-body dynamic simulations Hexagon D&E Adams MSC (Version 2022.1). The authors decided to develop their own simulation environment because of its significant lower computing time. Thus, their software is eligible to control the actual physical system in real time in the future.

The remainder of this paper is organized as follows: Section 2 . describes the concept of the surface. Section 3 explains the analytic model used to describe the working principle. Section 4 describes the simulations of the system for different tasks. Finally, Section 5 reports the results, Section 6 consists in the validation of the latter, and Section 7 describes the conclusions.

#### **2. Concept Description**

This section provides a qualitative description of the functioning of the surface and the modules of which it is composed. The operating principle of the surface is based on an array of modules with an idle rolling element and an orientable axis of rotation. These axes are used to generate directionable friction forces on the body in contact with the surface. The origin of the friction forces lies in the very working principle of a rotor. In fact, these are constructed and used (e.g., conveyor rollers) to facilitate motion in one direction (perpendicular to their axis of rotation). This leads to a difference in the friction forces exchanged with an overlying object.

In practice, there will be a smaller component, perpendicular to the axis, caused mainly by rotational friction, and a larger one, parallel to the axis, due to linear friction (*Fperp* and *Fpar*, respectively, in Figure 2c). Therefore, the main friction component will be predominantly along the axis, so by orienting the latter, the main force will also be directed. The choices of the axis directions are limited to a fixed number. For example, the simplest configuration is with two orientations (0, 90°), whereas more complex designs are with four (0 ◦ , −45◦ , 45◦ , 90◦ ) (Figure 2b) or more. In Figure 2b, the four orientations are shown, with the rotors displayed like cylinders (the projection is a rectangle) and the axes like dashed lines. For applications such as slowing or stopping the material flow, two orientations may be enough, whereas for sorting, four orientations should be taken into account. Therefore, modules with four directions will be considered in the following pages. For its intended applications, the module must be used together with others to create the surface (Figure 2a), so the sum of the forces exerted by every unit will control the objects for sorting and feeding purposes. In Figure 2a, the surface setup is summarized, showing the grid of modules with the object on top, and the arrows inside each cell represent directions of motion promoted by them. As noted in the introduction section, the system is under-actuated, so the surface can only slow down the object while it completes its task of sorting or feeding. If, for example, a body is simply placed on the surface without external actuation, it will not move. To overcome this situation, the surface plane can be tilted to use the gravity effect as the missing actuation, or an initial velocity can be given to the part before crossing the surface. Depending on the application and the task of the system, one solution may be preferred over the other, or a combination of both may be chosen. For instance, if the extension of the surface is big, because the objective is to orient a body, the tilted solution could be better, whereas, if the goal is to sort parts in a conveying line, the initial velocity of the object, given by the conveyors before the sorting area, is enough.

In addition, similar to other active surfaces that use rotors for intralogistic purposes, the application is limited to parts with at least one planar surface, in order to have simultaneous contact with three or more rotors. Another limit will be the maximum load per module, but this is related to the resistance of the device, which is not studied in this article. However, the study conducted is done considering the dimensions of the modules and their resistance comparable with the similar actuated existing systems. To provide a reference, the base of a transported object has the minimum dimensions of *b* × *h* = [20 cm × 20 cm], and the weight bearable by a single module is 20 kg.

3

6

N

6

4

5

6

6

5

4

3

Fpar

Contact point Fperp

Surface

A A Progettato da Controllato da Approvato da Data Data **Figure 2.** Illustration of: (**a**) the modular active surface, (**b**) the four possible orientations of the rotors, and (**c**) the three contact forces between the object and the rotor.

Module

Disegno1

x

x

y

x

Rotor

45

x

45 -

ebianchi 09/05/2022

1

Edizione Foglio

2

2

y

Fperp

y

Rotor axis

Rotor axis

Fperp

Fpar

Fpar

Friction force

Friction force in the axis direction

Friction force in the rotation direction

Friction force in the rotation direction

Velocity of the object on the contact point

Velocity of the object on the contact point

2

2

ebianchi 20/05/2022 Progettato da Controllato da Approvato da Data

ebianchi 20/05/2022 Progettato da Controllato da Approvato da Data

Disegno2\_gamma

Disegno2\_gamma

Data

Data

3

3

A A

A A

B B

B B

1

1

H

H

E

x

E

*v*p

*v*p

x

D

D

*v*par

*v*par

*v*perp

*v*perp

1

1

1 / 1 Edizione Foglio

1 / 1 Edizione Foglio

Object transported

#### **3. Analytical Model**

B B Rotor axis Rotor 1 2 3 4 5 1 / 1 Continuing with the description of the functioning of the surface, this section reports the hypothesis, the forces involved to move the objects, and the model applied to determine these, respectively, in the subsections *Assumptions*, *Friction Forces Model*, and *Equilibrium Equations*. As stated in Section 2, the working principle of the module is based on the contact forces exchanged between the rotor and the transported object. With this in mind, the analysis starts from the actions on a single module and then expands to the whole surface and to the equations describing the motion of an object. 3 4 5 6 3 4 5 6 D D

Figure 3 introduces the notation used in the next pages. In particular, [*xg*, *yg*, *zg*] are the coordinates of the object's center of mass, *G*, according to the reference frame {*O*, *x*, *y*, *z*}, and *G*e is the projection of *G* on the object base along the *Z* axis (Figure 3a). D D

**Figure 3.** Schematic illustration of the object on the surface: (**a**) lateral view, (**b**) upper view.

4

4

5

5

6

6

#### *3.1. Assumptions*

Listed here are the fundamental assumptions that allowed the definition of the analytical model describing the motion of an object on the surface.

$$\min\_{\underline{X}} f(\underline{X}) \text{ such that} \begin{cases} A\_{eq} \cdot \underline{X} = b\_{eq} \\ \underline{l\_b} \le \underline{X} \le \underline{u\_b} \end{cases} \tag{1}$$

where the first vector equation represents the equilibrium in the vertical direction and around the two axes of rotation in the plane, and the second vector equation states the boundary conditions. The function to be minimized *f*(*X*) is the Euclidean norm of the *Ni* array. The single terms are:

$$A\_{eq}(\mathfrak{J}\times n) = \begin{bmatrix} 1 & 1 & \dots & 1 \\ \mathcal{Y}\_{\tilde{G}\mathcal{P}\_{\mathbb{P}}} & \mathcal{Y}\_{\tilde{G}\mathcal{P}\_{\mathbb{P}}} & \dots & \mathcal{Y}\_{\tilde{G}\mathcal{P}\_{\mathbb{P}}} \\ -\mathbf{x}\_{\tilde{G}\mathcal{P}\_{\mathbb{P}}} & -\mathbf{x}\_{\tilde{G}\mathcal{P}\_{\mathbb{P}}} & \dots & -\mathbf{x}\_{\tilde{G}\mathcal{P}\_{\mathbb{P}}} \end{bmatrix}, \quad \underline{b\_{eq}}(\mathfrak{J}\times 1) = \begin{bmatrix} mg\cos(\gamma) \\ & -mg\mathbb{Z}\_{\mathcal{S}} \\ m\mathbb{R}\mathbb{Z}\_{\mathcal{S}} - mg\mathbb{Z}\_{\mathcal{S}}\sin(\gamma) \end{bmatrix}.$$

$$\underline{X}(n \times 1) = \begin{bmatrix} N\_1 \\ \cdots \\ N\_l \\ \cdots \\ N\_n \end{bmatrix}, \quad \underline{l\_b}(n \times 1) = \begin{bmatrix} 0 \\ 0 \\ \cdots \\ 0 \end{bmatrix}, \quad \underline{u\_b}(n \times 1) = \begin{bmatrix} mg\cos(\gamma) \\ mg\cos(\gamma) \\ \cdots \\ mg\cos(\gamma) \end{bmatrix}, \quad f(\underline{X}) = \|\underline{X}\| $$

where *n* is the number of rotors under the object, *m* is the mass of the object, *g* is the gravity acceleration, and *γ* is the inclination of the surface compared to the ground (Figure 3a). In addition, *xGP*e *i* and *yGP*e *i* are the *x* and *y* components of the vector that connect *G*e to the *i*th contact point *Pi* .

View A

4

X

xg

Rotor

Z Y

yg

G

șg

6

6

5

mgSin(J

Object

View A

mgCos(J mg

Zg

Surface

G

6

X

Z

Y

J

6

5

4

3

3

y

z

6

6

avoided to reduce complexity and the computational cost (for the simulation) without adding relevant contribution, as this friction component is very small and almost negligible compared to the axis direction component. In conclusion, this assumption, introducing the friction model, results in the boundary condition *F* = *µN* for the calculation of friction forces.  Friction Forces Model*

In this subsection, considering the previous assumptions regarding sliding and pure rolling and taking into account that the friction forces oppose motion, the directions of the two forces (*Fpar*, *Fperp*) and their final formulation are determined. Thanks to the third assumption, these directions are related only to the object velocity at the contact point with the rotor, −→*Vp* = [*x*˙ *p*, *y*˙ *p*]; an example is shown in Figure 4a.

B B B B **Figure 4.** Schematic representations of: (**a**) velocities and forces exchanged in a contact point, (**b**) the calculation of the velocity in the contact point, (**c**) the resulting forces from the contacts.

A A ebianchi 20/05/2022 Progettato da Controllato da Approvato da Data Data In particular, the velocity −→*Vp* can be split into a perpendicular and a parallel component (Figure 4a), and the two forces *Fpar* and *Fperp* are in the opposite directions of these velocities. In Figure 4a, three more angles are introduced: *ε*, *β*, and *α* (all angles are defined positive in the anti-clockwise direction from the *x* axis):

3

3

4

5

5

4

The velocity of the object's contact point is determined using the formula of the kinematics of a rigid body: −→*Vp* = [*x*˙*g* ˆ*i*, *y*˙*g* ˆ*j*] + ˙*θg* ˆ*k* × −→ *GP*e , where [*x*˙*g*, *y*˙*g*] and ˙*θg* are the linear and angular velocity of the object, and −→ *GP*e is the vector that connects *G*e to the contact point *P*. The velocity identification is also reported in Figure 4b, where the piece is schematized with the shape of the base, and the rotors with a grid of points. This scheme presented in Figure 4b is also useful in visualizing how the rotor–object contact points are determined. In fact, knowing the values [*xg*, *yg*, *θg*] and the mathematical formulation of the base area with respect to these, it is possible to know the set of coordinates of the points inside the base perimeter. More specifically, given a planar reference frame {*G*e, *x*e, *y*e} with coordinates *x*e, *y*e, centered in [*xg*, *yg*] (*G*e) (Figure 4b) and oriented by *θg*, let *Q* be a set of geometrical parameters that depends on the shape of the object base (e.g., radius for a circle or base and height for a rectangle) and let *GD* ∈ *Q*. Let Λ be a set of indexes. For every *λ* ∈ Λ, let *Fλ* : R2 × *Q* → R be functions such that *W* = {(*x*, *y*) : *Fλ*(*x*, *y*, *GD*) = 0, ∀*λ* ∈ Λ} is the perimeter of the object base. The area of the object base, *D*, is defined as *D* = {(*x*e, *y*e) : *Fλ*(*x*e, *y*e, *GD*) ≤ 0, ∀*λ* ∈ Λ}. Therefore, given the *xR*, *yR* coordinate pairs that define the position of the rotors in the array (which are fixed), it is verified whether or not they belong to the set, i.e., if *xR yR* = *xg yg* + *Rθg x*e *y*e has solution for a pair of *x*e, *y*e ∈ *D*, where

*Rθg* = cos *θg* − sin *θg* sin *θg* cos *θg* is the rotation matrix of the reference frame with coordinates *x*e, *y*e and centered in *G*e. Now all the elements to describe the two forces in the plane of

the surface (*Fpar*, *Fperp*) have been introduced and the result, according to Figure 4a, is summarized in Equation (2a,b).

$$F\_{par\_i} = \mu\_{max} N \* \tanh(k|V\_p|) \* \text{sgn}(\cos(a)) \tag{2a}$$

$$F\_{perp\_i} = \mu\_{\min} N \* \tanh(k|V\_p|) \* \text{sgn}(\sin(\alpha))\tag{2b}$$

where the two friction coefficients are, respectively:

The *sgn* function models the direction of the forces, according with the orientation of the object velocity in the contact point, whereas *tanh* and the coefficient *k* (selection of this parameter explained in Section 6) smooth the discontinuities ([41]). The perpendicular force would be ideally 0, however, it is generally small as inertia forces are neglected and only the rotational friction of the rotor around its axis is considered. The parallel friction force can reach significantly higher values and is opposed (in the axis direction) to the object velocity. Therefore, it is the dominant driver to manipulate the object.

#### *3.3. Equilibrium Equations*

This subsection shows how the equations of motion of the object on top of the surface are obtained from the contact force components. In practice, once these components are calculated for each contact point, they can be added together to obtain the total reactions on the part and its equilibrium equations. To achieve this, first, the friction forces in the plane at the various contact points are decomposed along the *x* and *y* directions, knowing the arrangements of the rotors (*εi* ):

$$F\_{\mathbf{x}\_i} = -F\_{par\_i}\cos(\varepsilon\_i) - F\_{per\_i}\sin(\varepsilon\_i) \tag{3a}$$

$$F\_{y\_i} = -F\_{par\_i} \sin(\varepsilon\_i) + F\_{per\_i} \cos(\varepsilon\_i) \tag{3b}$$

The signs introduced for these forces are in accordance with Figure 4a,c. The resulting equilibrium equations are summarized in Equation (4).

$$\begin{cases} \pounds\_{\mathcal{S}} = \frac{1}{m} \left( (\sum\_{i=1}^{n\_i} F\_{\mathbf{x}\_i}) + mg \sin(\gamma) \right) \\\\ \mathcal{Y}\_{\mathcal{S}} = \frac{1}{m} (\sum\_{i=1}^{n\_i} F\_{\mathbf{y}\_i}) \\\\ \not\not\varrho\_{\mathcal{S}} = \frac{1}{f} (\sum\_{i=1}^{n\_i} ((x\_{P\_i} - \mathbf{x}\_{\mathcal{S}}) F\_{\mathbf{y}\_i} - (y\_{P\_i} - y\_{\mathcal{S}}) F\_{\mathbf{x}\_i})) \end{cases} \tag{4}$$

where *x*¨*g* , *y*¨*g*, and ¨*θg* are, respectively, the two linear and the angular accelerations of the object; *xpi* and *ypi* are the coordinates of the *i*th point of contact; and *J* is the moment of inertia of the object around the axis perpendicular to the surface plane. Considering that also *J* can be expressed as a function of the mass, it is possible to notice that the inertia affects the system only from a structural point of view. In fact, as it is visible from the Equations (1), (2a,b), and (4), *m* can be simplified.

#### **4. Simulations**

In this section, the simulation process carried out with the analytic model implemented in MATLAB is shown. As a first step, this simulation environment permitted verifying in a general way that moving objects with the proposed theory and concept is possible. Second, it allowed several tests to be carried out in order to prove the usability of the system for some common intralogistics applications, such as sorting, orienting, stopping, and slowing down material flows.

The simulation takes advantage of an iterative loop (shown in Figure 5a) that involves the following steps: initialization, rotor identification, stopping criteria, contact forces and equilibrium computation, and two steps of integration. The initialization of the problem is made at the beginning by providing data from the object and the surface (example: object mass, base shape, initial position and velocity, friction coefficients, module pattern, inclination of the surface, rotor orientations, etc.). After that, the iterative loop can start. For each step (e.g., the *s*th step), the rotors below the object are identified (thus the contact points as well), knowing the position of the object from the previous iteration ((*s* − 1)th) and the disposition of the modules (as introduced in Section 3). Once the rotors below the object and their orientation are determined, the object velocity from the previous iteration ((*s* − 1)th) together with the inclination of the surface are evaluated by the stopping criteria: if "*x*˙*g* ≤ 0 m/s and *γ* = 0 ◦", the object is stopped because of negative velocity and null inclination and the loop ends. This is true assuming, in general, that the initial conditions always provide an input towards positive *x* (*γ* ≥ 0° or *x*˙*g*0 ≥ 0). When the stopping condition is not reached, the loop continues and the forces in the contact points can be computed (still considering the object velocity from the previous iteration ((*s* − 1)th)). With these forces, the equilibrium of the body and the accelerations are calculated according to Section 3. Finally, in order to obtain the velocity and the position, two integration steps of the acceleration vector are implemented. The derived values are the input of the next iteration ((*s* + 1)th) of the loop.

So far, the model seems to represent the operation of the surface adequately when the body velocity is greater than zero in the *x* direction and for the stopping condition. However, when, as an example, the object is about to start from a standstill with rotors and surface inclined, an undesirable phenomena such as reversal of motion (negative *x*) can occur. This happens because, in the friction assumptions, the static condition is not initially considered (Equation (2a,b)). However, reverse motion is obviously not possible in reality, then, in practice, when the condition of *x*˙*g* ≤ 0 m/s is achieved (*γ* > 0°, otherwise the object is stopped as described before), the process needs to be adjusted. The logic instructions to solve this problem and at the same time implement the analytical model of Section 3 are summarized in the block diagram of Figure 5b, which in practice is executed inside "Friction forces & Object equilibrium" block of Figure 5a. In detail, the process works as

follows: first, the *Ni* terms are calculated as explained in Section 3, then, since the velocity along *x* is known from the previous iteration, the condition "*x*˙*g* ≤ 0 m/s" is verified. If it is false and "*x*˙*g* ≥ 0 m/s", there is no problem of standstill or stopping and the procedure continues as described in Section 3, thus, friction forces calculation (Equation (2a,b)) and equilibrium ((Equations (3a,b) and (4)). In contrast, if the velocity is less or equal to zero ("*x*˙*g* ≤ 0 m/s" is true), an initial positive speed is assigned to the object *x*˙*g* = 0.0001 m/s, the friction model is applied, and it is verified if the gravity effect is stronger than the friction forces in the *x* direction (" *ni* ∑ *i*=1 *Fxi* + *mg* sin *γ* > 0"). At this point, if gravity wins, the object is moving according to the process defined before (Equations (3a,b) and (4)); in contrast, if gravity is not enough, the friction forces in the *x* direction will be of the same magnitude of the gravity effect (" *ni* ∑ *i*=1 *Fxi* + *mg* sin *γ* = 0"), and the displacement will be in the *y* direction (always according to Equations (3a,b) and (4)). To clarify the diagram of Figure 5b, the condition " *ni* ∑ *i*=1 *Fxi* + *mg* sin *γ* = 0" permits the calculation of the *Fpari* terms, while *Fperp* is the same as the previous case. For instance, this procedure permits us to obtain the motion of the object when it is placed without an initial velocity on the inclined surface, with the rotors oriented at 45°. In fact, with a sequence of displacement in the *y* and *x* directions, the movement is achieved.

**Figure 5.** (**a**) Iterative loop scheme for the simulation. (**b**) Block diagram of the logic process behind the model to avoid the reverse motion and to calculate friction forces (Equation (2a,b)) and equilibrium equations (Equations (3a,b) and (4)).

The simulations conducted to test the surface are divided into three main categories, according to the application of the system:

The first category of simulations concerns sorting. The objective of the setup and the sorting itself is to divert an object from the transport line. The layout considered for the simulation involves an array of modules, of which the subset performing the task has the rotors with the rotation axes inclined at ±45° (in Figure 6a inclined at +45°). The area

outside the line, where the objects are directed, is simulated as an array of low-friction (*µexsur f*) support points, distributed as the sorting modules. In this zone, the friction is totally opposite to the velocity of the object, without considering any rotor inclination. The sorting is assumed achieved when the body moving on the line is deflected in such a way that it is only in contact with the elements outside the line.

The second simulated application is the slowing activity. The rotor arrangement consists of a few modules within a line to create a controllable friction on the object to modulate its speed, without pushing it to the sides. The layout in this case can use rotors oriented with their axes in the direction of the flow (*εi* = 0 ◦ ) (Figure 6b) or inclined at +45◦ in a row and −45◦ in the other (Figure 6c).

The last application concerns stopping the motion of an object. This task is the extreme application with respect to the previous application of slowing down. In fact, the modules are still placed with their axes in the flow direction (Figure 6b) or oriented at ±45◦ (Figure 6c), but aim to stop the object.

These three types of simulations do not yet include either real-time control of the rotors or position or trajectory tracking for the object, because, as already indicated, the authors' objective in this paper is to verify the mechanical functioning of the surface. However, it is possible to imagine that by having sensors that recognize the object arriving at the surface, the modules can pre-arrange the rotors, as shown in Figure 6, depending on the application or sorting direction. This lays the foundation for using the system as a smart surface.

**Figure 6.** Examples of simulation layouts for: (**a**) sorting, (**b**) slowing or stopping with *εi* = 0 ◦ , and (**c**) slowing or stopping with *εi* = ±45◦ .

Figure 6 introduces the symbols used in the following sections to display the layouts of the simulations implemented in MATLAB: the red elements correspond to the transport line, the blue dots correspond to the area (array of support points) out of the line, and the green square contour corresponds to the object base. Regarding the transport line, dots represent rotors centers, and lines represent the rotors axes. The red lines before the active area (the green square) are oriented perpendicularly to the flow line and their effect is to reduce in a minimum way the motion of the box in the flow direction, and the inclined lines (*εi* 6= 90◦ ) simulate the sorting, the slowing, or the stopping surface.

The fixed initial parameters, which are used in all the simulations, are shown in Table 1. In particular, the first four parameters are about the dimension and the inertia of the object, whereas the following three are the friction coefficients. These last values have been chosen by making the following considerations: *µmax* must be a medium-high friction value (assumed *µ* = 0.5, because is similar to the kinetic friction coefficients between paperboard and the conveyor belt in [42,43]), as it models the sliding of the object on the rotor in the direction of the axis; *µmin* must be a low value (*µ* ≤ 0.1), as it models the rotational friction of the rotors (assumed *µ* = 0.01, similar to a rolling friction coefficient); and, *µexsur f* has to be a low value as well (*µ* = 0.005 is selected), as it models the area outside the line where one can imagine having load-bearing spheres supporting the material (according to the *Omnitrack* catalog [44], *µ* = 0.005 ÷ 0.03). Actually, the coefficients described would depend on the materials in contact, which have not yet been defined. However, the exact values are not relevant for the purpose of proving the surface capabilities; the important thing is that

*µmax* > *µmin*, *µexsur f* is maintained. Finally, the last parameters are the distances in the *x* and *y* directions between two rotors centers and the object initial position and acceleration. Exceptions to these starting conditions are indicated with the results for each particular case, together with the missing parameters such as initial velocity of the object [*x*˙*g*0, *y*˙*g*0, ˙*θg*0] and the inclination of the surface *γ*.

**Table 1.** Initial fixed parameters for the MATLAB simulations.

#### **5. Results**

This section presents the simulation results for all three categories: sorting, slowing, and stopping. For each category, the authors illustrate their results for relevant parameter sets and indicate a practical use of the program for the real setup design.

#### *5.1. Sorting*

This subsection reports the results of the first of the three categories, sorting. Initially, graphical examples of the function under varying input conditions are shown, then an application for the design of the actual system is presented, and finally a comparison with existing sorting systems is given. Figure 7a,b show the sorting towards the two directions, *y* > 0 and *y* < 0, respectively, with an initial object velocity *x*˙*g*0 = 1.5 m/s and without surface inclination. Figure 7c, on the other hand, shows the sorting at different tilting angles of the surface and without an initial velocity. Each plot displays the trajectory of *G* (object's center of mass) during the sorting and the position and orientation of the object at the end of the simulation time.

In Figure 7a, the green square and its trajectory show that, with rotors inclined at −45◦ , sorting towards *y* > 0 can be achieved, as the object is shifted completely off the line.

The opposite sorting condition was tested as well; the only difference for the initial data was the direction of the rotors (+45°). The result is shown in Figure 7b and the trajectory is mirrored to the first one. In fact, as expected in reality, the ending displacements of the object are equal in modulus and with opposite sign for the *y* and *θ* displacements (Table 2).

Figure 7a,b show the capability of the surface to deflect the object trajectory for sorting purposes. Furthermore, coherent results were achieved for both sorting directions, for a reasonable initial velocity, and within one second of simulation time. From now on, since the mirroring of the results has been demonstrated for a symmetric rotor arrangement, the outputs are shown only for one sorting direction.

**Figure 7.** Illustration of the trajectories and the final position and orientation of the object during: (**a**) the sorting towards *y* > 0 (up) task, (**b**) the sorting towards *y* < 0 (down) task, (**c**) the sorting towards *y* < 0 (down) with different *γ*.

After this first presentation of the sorting capability for the under-actuated system, instead of using an initial object velocity on a horizontal surface, the sorting task was studied with an initially stationary object on an inclined surface with the identical rotor arrangement (Figure 7c).

In Figure 7c, the first case (green) proves that, without any external actuation, the object is not moving in the simulation, as in the real word. In the second case (black square, Figure 7c), with the surface inclination of *γ* = 5 ◦ , the body moves, but for the chosen simulation time (*t* = 1.5 s) it has not yet completed its sorting. This can be seen from the fact that the black square contour still contains the center of a sorting rotor (red). Finally, the third case (magenta square), where the surface inclination is *γ* = 10◦ , reports a completed sorting activity. Summarizing the black and magenta bodies in Figure 7c demonstrate that it is possible to transport and sort parts by inclining the surface. In particular, the more it is tilted, the faster the object moves, and the quicker the sorting is achieved. In Figure 7c, the values of the position and orientation reached at the end of the simulation (*tblack* = 1.5 s) for the *γ* = 5 ◦ case (black line) are [*xg*, *yg*, *θg*]*black* = [0.369 m, −0.171 m, −36.05◦ ] (Table 2), while those of the magenta line for the same value of *yg* = −0.171 m are [*xg*, *yg*, *θg*]*mag* = [0.402 m, −0.171 m, −28.62◦ ] and the simulation time for which they are obtained is *tmag* = 0.917 s. Comparing these values, it can be seen that when the surface is more inclined (magenta line), it requires more space in the *x* direction to perform the sorting *xg*,*mag* > *xg*,*black*, but the time required is less *tblack* − *tmag* = 0.583 s. In addition, different *γ* also produce different final orientations, in this case smaller for larger *γ*: *θg*,*mag* < *θg*,*black*.

Table 2 summarizes the results and the initial data about the previous sorting simulations.

**Table 2.** Specific initial data and results for the sorting simulations.

A practical use of the simulation is to determine the minimum number of modules necessary for sorting before building the real system. In general, as it was already presented,

this result is affected by various initial conditions. However, it can be interesting to know how many columns of rotor units are necessary to successfully complete the different tasks for different object sizes and velocity ranges. In order to derive this relationship, a set of simulations can be performed. Table 3 reports the results of the discussed analysis for two object sizes: *b* × *h* = 0.25 m ×0.25 m, called "S-box", and *b* × *h* = 0.45 m × 0.45 m, called "B-box". "B-box" and "S-box" have the same mass (*m* = 10 kg) and height (*zg* = 0.1 m).

The counting of the number of columns for the sorting is done starting from the first column with the inclined rotors to the last one touched by the object before being out of the line. As is visible from Table 3, the sorting with low initial velocities is not always possible: with *x*˙*g*0 less than 0.8 m/s for "S-box" and less than 1.1 m/s for "B-box", the boxes come to a stop within the transport line. Despite this, the speed values considered in the simulations are within the most common operating ranges of conveyors; however, any speed can be implemented in the software. In general, the missed sorting occurs because the deflecting forces also tend to slow the object, which may stop before reaching the target. In addition, the displacement to achieve to be out of the line is related to the dimension of the object, so the bigger it is, the more absolute transversal displacement is required to achieve the sorting. Summarizing, thanks to the simulation described, it is possible to realize such tables depending on the specific application, considering that the results are influenced by the layout of the modules and the initial position. In particular, for Table 3 the "S-box" layout is the same of the sorting simulations of Figure 7b, but with more columns of sorting rotors in order to test a bigger range of velocities. Instead, for the "B-box" there are two more rows of rotors in the sorting line, always with the rotors oriented for sorting towards *y* < 0 (*εi* = +45◦ ), making a total of four instead of two lines.

The data from the simulations carried out to produce Table 3, coupled with further simulations with, as input, only the inclination of the surface, allowed for a comparison of performance with similar existing sorting systems [45]. Figure 8 presents these results. The red dots represent the values of the sorting rate when only the initial velocity is set as input, whereas the blue dots are for when only the inclination of the surface is exploited. Figure 8 shows the most common sorting capacity ranges [45] and, as can be seen, the system proposed by the authors guarantees a medium rate for most input values, with few exceptions in both directions, i.e., high capacity and low capacity. The rate is defined as objects sorted per hour (*pcs*/*h*) and can be easily calculated, considering the time needed to perform the sorting and a successful sorting condition. In practice, the condition of successful sorting is when *yg* ≥ *yg*,*limit* (for the rectangular shape: *yg*,*limit* = *h*/2 + p (*b*/2) 2 + (*h*/2) 2 ), i.e., when the displacement in the direction *y* is such that, whatever the orientation of the object, the base is no longer in contact with the sorting rotors. The input velocity values presented in Figure 8 are within the common ranges of use for sorting systems. In particular, Ref. [45] for devices similar to the surface proposed by the authors, but fully actuated and called "Torsional discs", defines input velocities between 0.5 and 1.5 m/s. "Torsional discs" have sorting capacities ranging between 1600 and 4500 pcs/h when the size and weight of the transported objects are comparable to those considered in the authors' simulations ("S-box" and "B-box" characteristics). The system developed in this paper, although it is under-actuated, also allows sorting rates within the same range, providing good performance and in line with the current technology. Obviously, given the under-actuation, there are minimum input speed thresholds, as shown in Table 3, however, they do not seem to limit the application ranges of the device. Tilting the surface and starting from zero initial speed requires longer sorting times and thus provides lower rates than those with velocity as input. However, according to Figure 8, the inclination also provides a predominantly medium sorting capacity, without the need for a prior conveyor. Therefore, the results obtained from the comparison of sorting performances showed that the surface introduced by the authors, with the advantage of the simplicity of the under-actuation, still guarantees sorting capabilities at the same level as current technology, both using initial speed and surface inclination as inputs.

**Table 3.** Minimum number of rotor columns required to achieve sorting with *γ* = 0°, but changing *x*˙ *g*0.

To sum up, the results of this subsection show that the surface is capable of sorting with the two different types of input (Figure 7 and Table 2). In addition, thanks to the developed simulation environment, it was possible to obtain the performance of the system proposed by the authors and compare it with the current technology. This showed that the surface provides a medium sorting capacity (Figure 8), in the same range as existing systems. Finally, for certain objects, the initial input speed was associated with the number of rotors required for sorting (Table 3). This was proposed as a further demonstration of the usefulness of the simulation environment for the design and control of the real system.

#### *5.2. Slowing*

The second set of results, here presented, is derived from the slowing simulations. The purpose of these simulations is to simplify and speed up the determination of the number of modules required for the slowing task, which, as with sorting, is linked to the initial parameters. Similar to the sorting results, graphical examples of the functioning and an application for the design of the actual system are given in this subsection.

In order to show the results of the slowing task, in this case the evolution of the object speed over time ([*x*˙*g*, *y*˙*g*, ˙*θg*]) is reported. Different from the sorting, here the object is not pushed away from the line, so the trajectory does not represent the course of the activity as simply as speed does. Considering the layout of the system, the only velocity component different from zero is in the *x* direction, which is the one plotted. Four different setups shown in Figure 9 were simulated to demonstrate the concept's operation. The graphs in Figure 10a,c show the results of the slowing simulations for the "S-box" and "B-box" cases, respectively.

**Figure 8.** Sorting rates obtained with MATLAB simulations for the "S-box" and "B-box" objects, with, as input, only initial velocity (red) and only inclination of the surface (blue).

**Figure 9.** Layouts of the simulations for the slowing and the stopping tasks with: (**a**) "S-box" *εi* = 0°, (**b**) "S-box" *εi* = ±45°, (**c**) "B-box" *εi* = 0°, (**d**) "B-box" *εi* = ±45°.

The initial data for this analysis are the same as in Table 1, with additionally *γ* = 0° and [*x*˙*g*0, *y*˙*g*0, ˙*θg*0] = [1.5 m/s, 0 m/s, 0 rad/s]. As visible from Figure 10a,c, both layouts are capable of reducing the speed of the object, fulfilling the slowing task. However, the speed drop is different with the orientation of the rotors. When the rotors are with *εi* = 0°, all the *Fpar* forces in the *x* direction cause a stronger slowing, whereas with *εi* = ±45° layout, where the *Fpar* are inclined like the rotor axes, the slowing is less. For these initial parameters and rotor dispositions, the ratios between the speed before and after the slowing rotors (*r* = *x*˙*gout*/*x*˙*g*0) are: *r S*−*box* 0 = 0.23, *r B*−*box* 0 = 0.18 and *r S*−*box* ±45 = 0.53, *r B*−*box* ±45 = 0.52. Additionally, Figure 10 shows the effect of *Fperp*, which is slightly reducing the velocity when the rotors are with *εi* = 90°. In fact, the plots in Figure 10a,c show that the lines are a bit inclined before and after the slowing area, which is included between the "Start slowing" and the "End slowing" dotted lines. It represents the contribution of rotational friction to the deceleration.

As for sorting, a parameters study was conducted within the simulation environment to derive the required number of rotor columns to achieve the slow down task of the objects. The results are displayed in Table 4. Additionally in this case, different initial speed and the two different object sizes were tested ("S-box", "B-box"). The reference layouts for the simulations are the ones in Figure 9. In Table 4, the counting of the columns number begins with the first column with inclination: *εi* = 0° or *εi* = ±45°, that the object encounters starting from the left. For example, considering Figure 9a and assuming that slowing is

verified, the result will be 2 columns (with *εi* = 0°). The target condition was set to at least halving the speed of the object while remaining a positive value after the slowing process (0.5*x*˙*g*0 ≥ *x*˙*gout* > 0). The two configurations (*εi* = 0° and *εi* = ±45°) were possible for the rotor layout, so the configuration with the smallest number of columns was chosen case by case to satisfy the slowing condition. As shown in Table 4 for both layouts and for all speeds analyzed, the number of rotors was identified. The introduction of the *εi* = ±45° configuration was necessary, as for some speeds the variant with *εi* = 0° generated a total stop of the conveyed object (e.g., for the first two values in Table 4). The *εi* = ±45° positioning also seems promising for the self-alignment of the conveyed material in the center of the transport line. In conclusion, in this subsection, the operation of the surface for slowing tasks was demonstrated (Figure 10). Additionally, as with sorting, certain input speed values of the object were associated with the number of rotors required for the task (Table 4). This provides data for the realization of the physical system and proves the usefulness of simulations for this objective as well.

**Figure 10.** Velocity trends for "S-box" during (**a**) the slowing and (**b**) the stopping, and (**c**) the "B-box" during the slowing and (**d**) the stopping with: *εi* = 0° (red), *εi* = ±45° (blue).

#### *5.3. Stopping*