Milling diagnosis using machine learning techniques toward Industry 4.0 Lorraine Codjo1,2, Mohamed Jaafar1, Hamid Makich1, Dominique Knittel1,2, Mohammed Nouari1 1 LEM3, GIP-InSIC, 27 rue d'Hellieule, 88100 Saint Dié des Vosges, France 2 Faculty of Physics and Engineering, 3 rue de l’Université, 67000 Strasbourg, France Corresponding author: knittel@unistra.fr Abstract cesses, they preferred BNs for their significant representa- tion capability and for the fast model building. Smart diagnosis of the milling in an industrial en- vironment is a difficult task. In this work, the di- agnosis using machine learning techniques has In this work, a smart milling diagnosis has been developed been developed and implemented for composite for composite sandwich structures based on honeycomb sandwich structures based on honeycomb core. core. The use of such material has grown considerably in The goal is to qualify the resulting surface flat- recent years, especially in the aeronautic, aerospace, sport- ness. Different algorithms have been implement- ing and automotive industries. Recent development pro- ed and compared. The time domain and frequency jects for Airbus A380 or Boeing 787 confirm the increased domain features are calculated from the measured use of the honeycomb material. But the precise milling of milling forces. The experimental results have such material presents many difficulties. shown that a good milling diagnosis can be ob- tained with a Linear Support Vector Machine The objective of this work is to develop an industrial sur- (SVM) algorithm: good accuracy and short train- face quality diagnosis for the milling of honey-comb mate- ing time. rial, by using supervised machine learning methods. Cut- ting forces are online measured in order to predict the 1 Introduction resulting surface flatness. The Industry 4.0 framework needs new intelligent ap- However, the literature's review does not exhibit deep proaches. Thus, the manufacturing industries more and studies related to the monitoring and the diagnosis of hon- more pay close attention to artificial intelligence (AI). For ey-comb core machining in order to ensure flawless sur- example, smart monitoring and diagnosis, real time evalua- face. tion and optimization of the whole production and raw materials management can be improved by using machine 2 Experiment description learning and big data tools [1]. An accurate milling process implies a high quality of the obtained material surface 2.1. Workpiece material (roughness, flatness) [2]. With the involvement of AI- The workpiece material studied in this investigation is based algorithms, milling process is expected to be more Nomex® honeycomb cores with thin cell walls. It is pro- accurate during complex operations. duced from aramid fiber dipped in phenolic resin (Fig. 1). The honeycomb cores consist of continuous corrugated T. Mikołajczyk et al. developed an Artificial Neuronal ribbons of thin foil bonded together in the longitudinal Network (ANN) for tool-life prediction in machining with direction. The aim of such a process is to create a structure a high level of accuracy, especially in the range of high allowing lightness and stiffness together thanks to the wear levels, which meets the industrial requirements [3]. hexagonal geometry of formed cells. Figure 1 illustrates The particularity of their work was the combination of a the geometric characteristics of the honeycomb core. multilayers ANN model with image processing in order to reduce the potential error. The use of honeycomb material in sandwich composite is D. Pimenov et al. evaluated and predicted the surface’s limited by the fragility of each wall of the honeycomb, roughness through artificial intelligence algorithms (ran- which influences the quality of obtained surfaces after dom forest, standard Multilayer perceptron) [4]: in their machining [7, 8, 9]. investigation the obtained performance depends on the The Nomex® honeycomb machining presents several defects related to its composite nature (uncut fiber, tearing parameters contained in the dataset. of the walls), the cutting conditions and to the alveolar M. Correa et al. compared the performances of Bayesian geometry of the structure which causes vibration on the networks (BN) and artificial neural networks for quality different components of the cutting effort [10]. detection in a machining process [5]. Even ANN models are often used to predict surface quality in machining pro- Honeycomb designation: A10-72-5OX Nomex® Honey- comb workpiece L direction Mounting Cutting force meas- system urement (Kistler) W direction Figure 3. Experimental test setup Density Cell size l Wall size t Angle α [kg/m3] [mm] [mm] [°] 72 5 0.08 100 Figure 1. Nomex® honeycomb cores and the main geometrical characteristics. It is clear that the use of ordinary cutting tools and also the mechanical and geometrical characteristics of honeycomb cores will have a crucial effect on machinability and on the quality of the resulting surface [11]. 2.2. Cutting tool and experimental environment For our study, the used milling cutter is provided from our Table 1. Machining center Realmeca® RV8 specifications industry partner the EVATEC Tools Company. In fact, ordinary cutting tools for machining honey-comb core 2.3. Milling experiments produce generally tearing of fibers and delamination of All experimental milling tests illustrated in this paper were cell structures. Subsequently, these cause a reduction of carried out on a three-axis vertical machining center Real- bond strength between the skin and the honey-comb core, meca® RV-8. For assessing the performance of the ma- and thus a weaker joint for composite sandwich structures. chining process of Nomex® honeycomb core we moni- As shown in figure 2, the EVATEC tool used is a com- tored and measured the cutting forces generated during bined specific tool with two parts designed to surfac- cutting, by using Kistler dynamometer model 9129AA. ing/dressing machining operation [12]. The first part is a The Kistler table is mounted below the Nomex sample in cutter body made of high speed steel with 16 mm in di- order to measure the three components of the machining ameter and having ten helix with chip breaker. This tool force as shown in figure 3. part is designated by hogger. The second part is a circular The milling experiment conditions are summarized in cutting blade made of tungsten carbide with a diameter of table 2. Four different speeds (high and low speeds) and 18.3 mm and having a rake angle of 22° and a flank angle four feed values were selected. of 2.5°. These two parts are mechanically linked to each other with a clamping screw. Spindle speed 2 000 10 000 15 000 23 000 (rpm) Feed rate 150 1 000 1 500 3 000 Hogger (mm/min) Table 2. Milling experiment conditions. 2.4. Measured signals Figure 4 shows the milling forces measured for honey- comb at 2000 rpm spindle speed and 3000 mm/min feed rate. Cutting forces are in the order of a few Newtons, they do not exceed 60 Newtons. Generally, the force in vertical direction (Fz) is quite small, thus, it is advised that to keeping vertical forces small in milling composite due to Figure 2. Milling cutter used for Nomex® honeycomb core the delamination issue. In our case, the vertical cutting ‘‘CZ10’’. force component is greater than other forces components which can be attributed to the mechanical properties of the Figure 3 shows the forces acquisition setup. During the honeycomb structure where the honeycomb structure is measurements, the x-axis of the dynamometer is aligned characterized by a better out-of-plane compression behav- with the feed direction of the milling machine and the ior than its tensile and shear strength. The evolution of longitudinal direction of the workpiece (parallel to core cutting forces shows significant oscillations, these oscilla- ribbons and the direction of honeycomb double wall). The tions are caused by vacuum in the cells of the honeycomb 3 orthogonal components of machining force (Fx, Fy and and the angle between the cutting direction and the honey- Fz) were measured according to figure 3 using the Kistler comb cell wall direction. table [12]. Our work focused on supervised learning for the classifica- tion of data according to predefined specific classes. 3.1 Features calculation The features are calculated in time domain and frequency domain from the raw signal represented on figure 4, in steady state behavior (transient zones, i.e. the time zone when the cutting tool enters into the honeycomb core and the zone when it exits, are not taken into account). The milling force plotted in figure 4 is the vertical force (Fz) of Figure 4. Milling force measurements for 2000 rpm spindle speed a given milling. That force is negative due to the fact that and 3000 mm/min feed rate: (a) during all process; (b) during the z-axis direction of the dynamometer has been oriented 0.2s (zoom) downwards Figure 5 shows the evolution of the surface quality (flat- ness) for various combinations of cutting conditions (spin- dle speed and feed rate). The defect of shape is higher for low speeds. Thus, for high feed rates that exceed the 1500 mm/min, the unevenness exceeds 500 µm which character- izes the severe tearing of the honeycomb walls. Figure 6. Measured milling force in time domain: (a) total data plot, (b) signal during steady-state phase Figure 5. Effect of cutting parameters on surface flatness. After a first data processing (filtering), firstly many fea- Given the low level of cutting forces, the quality of the tures are calculated in time domain for the measured mill- obtained machined surface allows to establish criteria for ing force signal called hereafter x(t). determining the machinability of the honeycomb struc- tures. The appearance of the uncut fibers is a machining The calculated time domain features are: defect specific to the composite material which depends on the type of the fibers and their orientation. The tearing of - maximum of x(t) (1) Nomex® paper, linked to the cellular appearance of the honeycomb structure, occurs under the effect of shear - minimum of x(t) (2) loading [12, 13]. - difference between the maximum of x(t) and the mini- mum of x(t) : amplitude range (3) 3 Milling diagnosis using machine learning techniques - median value of x(t) (4) There are many approaches in machine learning. The two principals are [15] : - Maximum of the absolute value of the signal : - Unsupervised approaches : based only on input data 𝑚𝐴𝑆 = 𝑚𝑎𝑥(|𝑥𝑘 |) (5) (unlabeled data) ; the goal is to find a natural grouping or structuring in the data set in order to reduce the number of - Interquartile range : observations ; 𝐼𝑄𝑅 = 𝑄3 − 𝑄1 (6) - Supervised approaches: based on input and output data where 𝑄3 and 𝑄1 represents respectively the upper and (labels). lower quartile. Supervised learning algorithms (with labeled data) can - Inter decile range : split in two categories [15] : 𝐼𝐷𝑅 = 𝐷90 − 𝐷10 (7) - Classification models which partition observations in where 𝐷90 and 𝐷10 means respectively the 90th ant the categorical groups (leads to a predictive model for discrete 10th decile. Both Inter quartile and Inter decile range are responses). a measure of statistical dispersion of the values in a set of - Regression models which describe the relationship be- data. tween outputs and variables through a mathematical func- tion (leads to a predictive model for continuous responses). - Average value of the signal : 𝑁 𝑁 1 𝑖 𝑚𝑒𝑎𝑛(𝑥) = ∑ 𝑥𝑘 (8) 𝑌(𝑘) = ∑ 𝑥𝑖 𝑒 −𝑗2𝜋𝑘 𝑁 , (𝑘 = 1, … , 𝑁) (20) 𝑁 𝑘=1 𝑖=1 where N is the number of samples of the signal x(t). - Average value of the absolute value of the signal : 𝑁 The frequency domain features are calculated for the Y(f) 1 𝑀𝐴𝑆 = ∑ |𝑥𝑘 | (9) signal. For example the Fast Fourier transform (FFT) plot 𝑁 is given on figure 7. 𝑘=1 - Average value of the absolute value of the derivative signal : 𝑁−1 1 𝑑𝑥𝑘 𝑀𝐴𝐷 = ∑| | (10) 𝑁−1 𝑑𝑡 𝑘=1 - Variance : 𝑁 1 𝑉𝑎𝑟 = ∑(𝑥𝑘 − 𝑚𝑒𝑎𝑛(𝑥)) (11) 𝑁 𝑘=1 - Energy of the signal : Figure 7. Measured milling force in frequency domain : FFT plot. 𝑁 𝐸(𝑥) = ∑ 𝑥𝑘 2 (12) All the calculated features (in time and frequency do- 𝑘=1 mains) have been normalized and stored in a table whose lines and columns respectively represent the physical ex- - Energy of the centered signal : periments and the associated feature values. 𝑁 2 𝐸𝑐 = ∑(𝑥𝑘 − 𝑚𝑒𝑎𝑛(𝑥)) (13) Normalization 𝑘=1 The resulting features being of different units and scale, it is necessary to normalize them in order to avoid meaning- - Energy of the derivative signal : less and redundant information. The chosen normalization 𝑁−1 is given by the equation 21 [14] : 𝑑𝑥𝑘 2 𝐸𝑑 = ∑ ( ) (14) 𝑑𝑡 𝑓𝑒𝑎𝑡𝑢𝑟𝑒 − 𝜇(𝑓𝑒𝑎𝑡𝑢𝑟𝑒) 𝑘=1 𝑓𝑒𝑎𝑡𝑢𝑟𝑒𝑛𝑜𝑟𝑚 = - Skewness : 𝜎(𝑓𝑒𝑎𝑡𝑢𝑟𝑒) (21) 3 where µ and σ represent respectively the mean value and 𝐸(𝑥 − 𝑚𝑒𝑎𝑛(𝑥)) 𝑆= (15) the standard deviation of each column of feature type. 𝑉𝑎𝑟 3/2 This normalization leads to : 𝜇(𝑓𝑒𝑎𝑡𝑢𝑟𝑒𝑛𝑜𝑟𝑚 ) = 0 - Kurtosis : { 4 𝜎(𝑓𝑒𝑎𝑡𝑢𝑟𝑒𝑛𝑜𝑟𝑚 ) = 1 𝐸(𝑥 − 𝑚𝑒𝑎𝑛(𝑥)) 𝐾= (16) 𝑉𝑎𝑟 2 The normalized features are dimensionless and can thus be compared. The normalized feature table contains 39 fea- - Moment order i (i = 5 : 10) : tures and 3 input parameters for each experiment: the rota- 𝑖 𝐸(𝑥 − 𝑚𝑒𝑎𝑛(𝑥)) tion speed, the cutting speed and the depth of cut (i.e. the 𝑚𝑖 = (17) quantity of material the tool will take during milling). 𝑉𝑎𝑟 𝑖/2 - Shannon entropy : Feature reduction 𝑁 It is necessary to reduce the number of features in order to avoid overfitting on one hand and on the other hand to 𝐸𝑆 (𝑥) = − ∑ 𝑥𝑘 2 ∗ log 2 (𝑥𝑘 2 ) (18) reduce the online computing time of the features. That 𝑘=1 dimensional reduction was made, by using the Principal - Signal rate : Component Analysis (PCA) [19]. PCA can be seen as a 𝑚𝑎𝑥(𝑥𝑘=1:𝑁 ) − 𝑚𝑖𝑛(𝑥𝑘=1:𝑁 ) data pre-processing method which leads to a weighted τ= reduced matrix Z where each principal component Zn (col- 𝑚𝑒𝑎𝑛(𝑥) (19) umn of Z) is a linear combination of all the original varia- bles (Xp) [19]: Secondly 19 features are calculated in frequency domain in a similar way for the measured milling force signal. There- fore, the Fast Fourier transform (FFT) of the signal x(t) has 𝑍 𝑛 = 𝛼 1𝑛 𝑋1 + 𝛼 2𝑛 𝑋 2 + ⋯ + 𝛼 𝑝𝑛 𝑋 𝑝 (22) been calculated : Based on a pareto plot of the principal component vari- ances (fig. 8), the reduced dataset can be obtained by keep- ing only the m-first principal components which allow to - Support Vector Machine (SVM) reach a variance percentage of 99%. In our experimental case, m=9. So only the nine first principal components are In order to evaluate how the internal parameters of each kept. algorithm influence their efficiency, several variants of the same algorithm have been implemented (various distances, different kernels, etc.). The first k-nearest neighbor (KNN) was implemented by keeping the default Euclidean distance. For the same mod- el, a limited number of neighbors (k=2) have been applied. Another training model consisted to weight each observa- tion (the rows of our data set). Moreover the KNN algo- rithm has been modified by using the Chebychev distance. The first used decision tree algorithm is a fitted binary classification decision tree. Then tree has been pruned to obtain a pruning tree of level 2. Figure 8. Pareto plot of the variance percentage. Two SVM algorithms have been implemented using dif- ferent kernel functions. The first one is a linear SVM 3.2 Labeled data which is the default function for a two-class data set. The second one is the Gaussian SVM algorithm which is a From the evaluation of the effect of the cutting parameters normalized polynomial kernel. on surface flatness results, we defined two classes of sur- face quality applied to the output data of each observation (see table 3) : 4 Obtained results Label Flatness (µm) Qualitative value 4.1 Results of the trained models Table 5 shows the accuracy result of each algorithm per- 'A‘ 0 – 600 Best surface quality formed with the complete normalized data set Algorithms Accuracy ‘B' 600 – … Worst surface quality (in %) KNN 100% Table 3. Label table for the experimental observations. KNN k=2 88.6% Weighted KNN k=2 83.4% As shown in table 3, we have two labels, called classes. Chebychev KNN k=2 100% Class 'A' corresponds to the positive class, Class 'B' corre- Tree 100% sponding to the negative class. News data will be predicted Pruned tree 66.67% using the rules below: Linear SVM 91.4% - If prediction probability result ≥ 0.5 : class A Gaussian SVM 91.4% - If prediction probability result < 0.5 : class B Table 5. Prediction error for the normalized data set. Training and validation set From table 5 we can observe that the simple decision tree The data were reduced in a training subset (used to train classifier leads to the best trained model for predicting new the classification algorithm) and a test subset (for the vali- data. But when that same algorithm was pruned, it lost an dation). For selecting the observations in each data subset, accuracy of 32.33%. a random logical selection was made. The table below Algorithms Accuracy resumes the partition of the data used in each classification kNN 99% kNN k=2 85% algorithm. Weighted kNN k=2 98.2% Dataset Percentage Class A Class B TOTAL Chebychev kNN k=2 87.5% Training 70% 15 20 35 Tree 99.2% Pruned tree 66.67% Test 30% 6 8 14 Linear SVM 99.8% (Validation) Gaussian SVM 93.8% TOTAL 100% 21 28 49 Table 6. Prediction error for the Label table for the n-by-n Table 4. Partition of the data. weighted matrix. Table 6 shows the accuracy result of each algorithm per- 3.3 Supervised learning formed with the reduced normalized data set (obtained In this work, several classification algorithms have been from the PCA). implemented in the Matlab software environment (with the A first conclusion is that for a dimension reduced feature Matlab Statistics and Machine Learning Toolbox, Version: table, the algorithms gave better results. From table 5 R2017) : (which presents the accuracy result from the weighted - k-nearest neighbor (kNN) reduced data set), it could be noticed that both three algo- - Decision tree (DT) rithms produces high accuracy rate for no-parametric algo- rithm. Specifically, linear SVM is the most efficient with the lowest running time and the highest accuracy. It is also required surface roughness by monitoring wear on noticed that more an algorithm is constrained by parame- face mill teeth” Journal of Intelligent Manufacturing, 29(5): 1045-1061, 2018. ters, more its performances are reduced. It is the case of [5] Correa M., C. Bielza, J. Pamies-Teixeira, the pruned tree for which we limit the expansion: it is “Comparison of Bayesian networks and artificial difficult for the model to find the best singleton split. neural networks for quality detection in a machining process”, International journal of Expert Systems with Applications, 36: 7270-7279, 2009. 4.2 Prediction results for new experiment data [6] A.I.H. Committee, ASM Handbook Volume 16: We used some news experimental data set in order to eval- Machining, ASM International, 1989. uate the performance of the trained model. The goal is to [7] J. Kindinger, Lightweight structural cores, ASM predict online (during milling) the surface quality. Results Handb. Met. Composites, 21: 2001. are presented here for the trained model by using the linear [8] D. Gay, Matériaux composites (translation: Composite SVM classifier algorithm. Table 7 shows the percentage of materials), Hermes, 2015. true positive rate and false negative rate obtained from the [9] R. Carl, “Three-dimensional honeycomb core machining apparatus and method”, US Pat. App. evaluation of the prediction. 13/707,670. 1, 2012. Predicted class [10] Jaafar M., S. Atlati, H. Makich, M. Nouari, A. Actual class A B Moufki, B. Julliere, “A 3D FE modeling of machining process of Nomex® honeycomb core : influence of the cell structure behaviour and specific tool A TP = 61.5% FN = 38.5% 100% geometry”, Procedia CIRP, 58: 505-510, 2017 [11] Mendoza M.M., K. F. Eman, S. M. Wu, B FP = 0% TN = 100% 100% “Development of a new milling cutter for aluminum honeycomb”, Int. J. Machine Tool Design Research (TP: true positive rate; FN: false negative rate; FP: false positive vol. 23: 81-91, 1983. rate; TN: true negative rate) [12] Jaafar M., H. Makich, S. Atlati, M. Nouari, B. Julliere, Table 7. Performance of the prediction using SVM classifier. “Etude expérimentale et numérique de l’usinage des structures composites en nid d’ abeilles Nomex®”, in National French congress : Journées Nationale sur les The negative class B was the best predicted class. This is Composites 2017. simply explained by the fact that it is the most predomi- [13] H. Tchoutouo and N. Gandy, “Adhesiveless nant class in the data set. honeycomb sandwich structure with carbon graphite The linear SVM algorithm loses in performance for data prepreg for primary structural application: a set with large predictors (i.e. large features). However, comparative study to the use of adhesive film” May with the reduced number of of features (9_first principal 2012. components), this algorithm has been the most accurate [14] Gopal S. and Kishore K., “Normalization: A Preprocessing stage”, CSE & IT department, VSSUT, algorithm with the best prediction rate for the lowest train- Burla, India 2015. ing time. [15] Shalev-Shwartz S. and Ben-David S., “Understanding Machine Learning: From Theory to Algorithms”, 5 Conclusion and further works Published by Cambridge University Press, 2014. [16] Rion J., Y. Leterrier, J.-A. E Manson, “Prediction of The milling's quality is qualified by the roughness or the the adhesive fillet size for skin to honeycomb core flatness of the resulted surface. In this work, different bonding in ultra-light sandwich structures” Compos. supervised learning algorithms have been implemented Part A, 39: 1547-1555, 2008. (offline) and compared. Each AI-based model has been [17] Ompusunggu A. P., B. Kilundu, T. Ooijevaar, S. applied to a set of features. These features were calculated devos, “Automated bearing fault diagnostics with cost from measured milling forces. effective vibration sensor”, WCEAM/VETOMAC From the prediction results, SVM algorithm seems to be 2017 Conference, Brisbane, Australia the most efficient algorithm in this application. [18] Madhusudana C.K., S. Budati, N. Gangadhar, H. The next step consists to implement an online version Kumar, S. Narendranath, “Fault diagnosis studies of (therefore the features have to be calculated online). face milling cutter using machine learning approach”, Journal of Low Frequency Noise, Vibration and Active Control, 35(2): 128-138, 2016. References [19] Richardson M., “Principal Component Analysis”, [1] Beskri A. et al., “Systèmes de surveillance course May 2009. automatique en usinage: Moyens et méthodes” [20] Fu Y., Y. Zhang, H. Qiao, D. Li, H. Zhou, J. Leopold, (translation: Automatic monitoring systems in “Analysis of feature extracting ability for cutting state machining: Ways and methods), French Mechanics monitoring using deep belief networks”, Procedia Congress 2013. CIRP, 31: 29-34, 2015. [2] IFPM-Formation, Usinage: Tournage Fraisage [21] Gao C., W. Xue, Y. Ren, Y. Zhou, “Numerical control (translation; Machining: Milling Turning), September machine tool fault diagnosis using hybrid stationary 2015. subspace analysis and least squares support vector [3] Mikołajczyk T., K. Nowicki, A. Bustillo, D. Yu machine with a single sensor”, MDPI conference, 7: Pimenov, “Predicting tool life in turning operations 346-358, 2017. using neural networks and image processing”, [22] Zhang C., X. Yao, J. Zhang, H. Jin, “Tool Condition Mechanical Systems and Signal Processing, 104: 503- Monitoring and Remaining Useful Life Prognostic 513, 2018. Based on a Wireless Sensor in Dry Milling [4] Pimenov D. Yu, A. Bustillo, T. Mikolajczyk, Operations”, Sensors journal, MDPI, 16: 795-815, “Artificial intelligence for automatic prediction of 2016.