Gaussian

A typical manufacturing process consists of a machining (material removal) process followed by an inspection system for the quality checks. Usually these checks are performed at the end of the process and they may also involve removing the part to a dedicated inspection area. This paper presents an innovative perturbation signal based data generation and machine learning approach to build a robust process model with uncertainty quantification. The model is to map the in-process signal features collected during machining with the post-process quality results obtained upon inspection of the finished product. In particular, a probabilistic framework based on Gaussian Process Regression (GPR) is applied to build the process model that accurately and reliably predicts key process quality indicators. Raw data provided by multiple sensors including accelerometers, power transducers and acoustic emissions is first collected and then processed to extract a large number of signal features from both time and frequency domains. A strategy for the selection of most relevant features is also explored in this work in order to reduce the input space dimension and achieve faster training times. The proposed GPR model was tested on a multi-robot countersinking application for monitoring of the machined countersink depths in composite aircraft components. Experimental results showed that the model can be used as a tool to predict the part quality through in-process sensory information, which in turn, helps to reduce the total inspection time by identifying the parts that would require further investigation.


Introduction
The recent trend towards the vision of the fourth industrial revolution encourages the use of information technology in manufacturing [1]. Therefore, conventional machining research mainly involving Computer Numerical Control (CNC) machine tools, is focused on the development of intelligent ways of exploiting information for process monitoring purposes [2,3]. Sensory data collected during the cutting process is typically used to infer several aspects of the process, such as the condition of the machine, identification of cutting tool wear or tool breakage [4,5], the prediction of workpiece surface roughness and anomalies [6] as well as diagnosis of common process failures [7]. Multiple sensor technologies including dynamometers, accelerometers, Acoustic Emissions (AE) and current/power transducers are applied to acquire relevant information from the machining process. Raw sensor data is usually filtered and processed in both the time and frequency domains for the extraction of several signal features that could potentially describe the process behaviour. These features are then used to build data-driven models based on machine learning approaches that can predict several process-specific performance variables. Common signal features include statistical descriptors such as the signal mean, RMS, variance, skewness, kurtosis and peak amplitude. For the frequency domain signal (obtained through the Fast Fourier Transform (FFT)), the above features refer to bands of the signal's spectrum [2,3].
In terms of modelling methods, the Artificial Neural Networks (ANN) have extensively been applied in machining research [3]. However, such models require a large amount of experimental data to use for training since the number of parameters to optimise is usually high. Some other learning approaches include the Support Vector Machines (SVM), fuzzy systems and probabilistic-based models. The latter is particularly suited to applications where the process has a highly stochastic behaviour and, in general they require less training data. A probabilistic approach that is currently getting more attention in process monitoring applications is the Gaussian Process Regression (GPR) [8]. A few examples include Teimouri et al. [9] reporting a GPR-based approach for structural health monitoring, Kong et al. [10] presenting the use of GPR for tool wear prediction and Obajemu et al. [11] applying GPR for modelling the uncertainty propagation in surface metrology applications.
The rising development and spread of robotic systems over the last decades, has motivated industries and researchers to look into robotic machining as a flexible and relatively low-cost alternative of the traditional machining approach. Industrial robots present great potential with features such as flexibility, high dexterity, large working areas, ease of programming and relatively lower costs. Yet, robotic machining is still in its early development phase because of several issues it faces [12,13]. The main challenges of this technology are due to the limited and variable stiffness of robots and their low positional accuracy, which inevitably introduce non-systematic (random) errors into the machining process that are difficult to anticipate. Tool deflections and vibrations are typically induced by the large process forces and the quality of the finished product is compromised [14].
The robotic machining research has addressed these challenges from the robotic structure perspective, proposing methods to directly minimise the sources of errors arising from the robot's low stiffness and positional accuracy. Several compensation strategies based on machine learning approaches have been reported in the literature, which rely on accurate measurements from optical laser trackers to build an error map of the robot positioning within its workspace. This information was then integrated into the robot path planner in order to compensate offline or online (by closed loop control schemes) the tool path deviations. Pan and Zhang [15] proposed a solution based on robot stiffness modelling and active force control strategy for real-time compensation of cutting tool errors. Xiong et al. [16] proposed a closed-loop error compensation method for robotic flank milling, which extracts the systematic machining errors using a spatial statistical analysis and then compensates the tool path based on mirror points of the error samples only. Reinl et al. [17] presented an entirely off-line simulation-based compensation strategy that predicts the trajectory errors in advance and uses the predictions to correct the robot path before the actual cutting occurs. Schneider et al. [18] proposed a 3D piezo-actuated high-dynamic compensation mechanism to use for robotic milling applications. Zhu et al. [19] presented an error compensation method based on ANN for robotic machining in aircraft manufacturing. The authors considered the reduced, two-dimensional input space of the workpiece surface to map the robot positioning errors in the Cartesian space. Probabilistic learning methods based on GPR have also been reported in the robotic machining literature. Nguyen et al. [20] proposed a GPR-based approach for modelling of the dynamic properties of an industrial robot within its workspace used in robotic milling applications. Wang et al. [21] presented a nonparametric statistical learning method for controlling the robot trajectory and contour tracking in robotic laser and plasma operations. The authors proposed an iterative method for compensation of both the feedforward torque and the motor reference applying GPR to model the inverse robot dynamics and kinematics, respectively.
In terms of the workpiece materials used in the aerospace industry, Carbon Fibre Reinforced Polymer (CFRP) composites have been widely applied in robotic operations such as milling, drilling and countersinking due to their excellent physical and mechanical properties [22]. The CFRPs have a high specific strength and stiffness, while being a light material with a high corrosion resistance, making them ideal to use in the manufacturing of aircraft components. However, due to its anisotropic and non-homogenous nature, this material faces numerous challenges related to its machinability. The most common problems that occur when improperly machining CFRPs include delamination, fibre pull-out, uncut fibres and surface anomalies. Therefore, particular attention is required in the selection of the process parameters, tool geometry and tool path planning to avoid or minimise the effect of these material-related issues and achieve good quality parts. From the process point of view, such errors can be quantified in terms of variations (imperfections) in the quality of the finished product (surface roughness, part geometry). Some examples of research works with further details on machining CFRPs using industrial robots can be found in Slamani et al. [23,24] and Qiu et al. [25].
In order to identify and potentially rectify the size of non-conformity products, an inspection step is typically applied after the machining operation in a production environment. This inspection is either included in the process operational cycle (automatic in-process inspection) or often it is performed (manually or automatically) at the end of the whole process using specific equipment. In either case, the process down time inevitably increases due to the time spent for the inspection task and/or delays associated with the preparation of the workpiece for inspection. Given the industrial focus on high productivity rates, it is essential therefore that this inspection time is kept at minimum. In order to further increase productivity without compromising the overall process quality, a method that indirectly infers the quality of the product (i. e. the result of the inspection step) from in-process information rather than the post-process direct measurements, is required. The method proposed in this paper aims to substitute the inspection step by predicting the quality of the finished part from sensory data collected during machining.
To summarise, while most of the process monitoring research focused on tool condition monitoring, surface roughness and/or fault diagnosis, this work investigates the use of process models in monitoring of other process output variables (depth of cut, tool path), which are subject to variations in robotic machining. In fact, process responses such as the actual depth of cut, feed rate and cutting tool path show much higher amount of uncertainty (deviations from the nominal values) in robotic machining compared to the conventional case. Consequently, the overall part quality and the process performance is compromised due to these deviations caused by the lower stiffness and tool positional accuracy of the industrial robots. Given the extra complexity that multi-joints articulated industrial robots and machinability of CFRP parts add to the manufacturing process, a suitable approach that considers the overall effect of several sources of errors in the process output variable is necessary. To this end, a probabilistic framework is proposed, which naturally deals with uncertainties associated with the input data (sensor signals) and output variables (process responses) of a data-driven model. This paper presents a novel data generation approach for the creation of robust machine learning models that predict part quality in a robotic countersinking operation. The data generation relies on the use of a random perturbation signal to the target depth of cut. The work also integrates the concepts of process monitoring such as signal processing, feature extraction and data-driven modelling into the robotic machining context, in the attempt to substitute the post-process inspection step.
The regression tests considered in this work aim to assess the accuracy and reliability of the proposed process model. The main objectives include: (i) investigation of a probabilistic learning approach for the prediction of the product quality; (ii) identification of a suitable method for the selection of relevant signal features and (iii) assessment of the model confidence in the predicted outputs.
The remainder of the paper is organised as follows: Section 2 defines the research problem and describes the case study of robotic countersinking process. Section 3 presents the details of the machine learning approach. Section 4 illustrates the methods applied in the experimental work for extraction and processing of the sensory information. The result of the regression tests are shown and discussed in Section 5. The paper conclusions are given in Section 6.

Problem formulation from a case study
The description of the case study of robotic countersinking process is presented first in order to introduce the problem formulation.

Robotic countersinking process
The robotic cell consisted of two industrial articulated robots that work together to countersink previously drilled holes in composite aircraft components. A close up view of the cell is shown in Fig. 1. Two KUKA robots are part of the cell: a large KR360 (master robot) having a countersinking end-effector equipped with an air-driven pressure foot and a (relatively) smaller KR180 (slave robot) with anvil end-effector to react process loads. Both robots used Kuka KRC2 controllers and ran RoboTeam software to facilitate the co-operation functionality of the cell and the end-effectors were controlled by a central Siemens Programmable Logic Controller (PLC) system. The objective of the cell is to automatically countersink aircraft CFRP panels, which have been previously drilled on a high precision machine tool.
The overall countersinking process flow can be briefly described as follows: both robots move in the vicinity of a pre-drilled hole and precisely locate over it; the slave advances slowly towards the panel until contact is made and then stops; the master activates the pressure foot installed at the end-effector to clamp the panel and then the cutting process commences. Both robots remain in position during cutting while the spindle drive advances the cutting tool towards the panel until the pre-set depth of cut target is reached. A schematic diagram of the process is presented in Fig. 2. The diagram illustrates four snapshots of the cutting cycle.
The countersink depth level illustrated in Fig. 2d defines the process output variable and the performance of the whole process is measured in  terms of its deviation from the nominal target depth. It is important to note that the angular deflections of countersinks are not a cause of concern for this particular process due to the robot positioning and the use of the pressure foot. In particular, a Lucana Aero vision system with two laser projections for measurement of the surface normality was applied to finely adjust the pose of the robots before the clamping process, ensuring that both robots were normal to the panel surface and perfectly aligned to the centre of the hole. The delamination problem of CFRP panels was also minimised due to careful selection of the cutting tool geometry and machining parameters (cutting speed and feed rate), optimised in a series of preliminary trials conducted outside the scope of this work.

Quality prediction problem
The Key Performance Indicator (KPI) of the above countersinking process is the final depth level obtained in correspondence of the maximum displacement of the cutting tool towards the surface of the panel, as indicated in Fig. 3 (left). In this work, the countersink depth was computed from an image taken by the inspection unit (Cognex camera) at the end of the cutting cycle, as shown in Fig. 3 (right). The image processing software, called DVT Intellect, is the proprietary one associated with the camera. The inspection routine consisted of finding the best circle fit for the outer perimeter of the machined hole (indicated by the large yellow circle in the figure). The red marks identify the area of the search: for the outer perimeter, the software searches the area between the medium and the large red circles. In order to get a higher image contrast for the inspection process, the panels were coated with white paint. In addition, the robots were equipped with swarf extraction units installed on the end-effectors to remove all the dust generated during cutting of the CFRP panels, leaving a clear image of the machined hole.
The following two responses were provided by the countersinking process: • In-process responses: Sensor signals collected during machining (monitoring system).
• Post-process responses: Outer radius of the machined hole (inspection unit).
The image taken from the inspection camera was processed to compute the countersink depth (D) using the formula D = (R − r)*tan(90 0 − α /2), where R and r are the outer and inner radii of the machined hole, respectively, and α is the cutting tool angle (in degrees), as indicated in Fig. 4.
The angle α is a characteristic of the cutting tool: here, α = 100 ∘ with a tolerance of 0.5 ∘ . The inner radius was known from the size of the predrilled holes (i.e. 6.35 mm with 5 µm of tolerance), hence r = 3.175 mm. This method was accurate (depth error accuracy of 0.05 mm), provided fast measurements, and could be implemented online for automated inspection of the holes. This is in contrast to the manual postprocess inspection, which can be both time-consuming and prone to human errors.
The quality prediction problem is the ability to accurately and robustly predict the post-process responses (countersink depths) from the data obtained by the in-process responses (sensor signals). Having a predicted value of the actual countersink depth available right at the end of the cutting process, can potentially reduce or avoid completely the need for post-process quality inspection (manual or automatic), saving valuable time.

Machine learning approach
This work exploits the experimental space defined by the target depth specified to the spindle controller (i.e. the controllable variable) to build and train data driven process models as quality predictors.

Perturbation signal based data for supervised learning
A schematic diagram of the process model based on a supervised learning approach for prediction of the countersink depth is shown in Fig. 5. The training phase is illustrated in Fig. 5a, while new pairs of sensor signals (input variables) and the corresponding countersink depth (output variable) are collected. Fig. 5b shows how the overall manufacturing process is different once the model has collected enough data (trained model), and it is ready to produce accurate predictions in place of inspections.
The output data is collected online by introducing a direct measurement (inspection camera). This in-process inspection routine adds extra time to the process cycle, but it allows the expansion of the training set as more cuts are performed. Eventually, the model will collect enough data to provide accurate depth predictions and, when this happens, the inspection step will not be required any more. It is also important to note that a supervised learning approach would need a relatively high number of training data (i.e. multiple runs of the process) in order to accurately map the input-output relationship. Moreover, since this data needs to be a representative example of the process behaviour, a relatively large variety of cutting trials (from under-to over-target holes, as well as in-target holes) is necessary.  In this work, the required data variability is ensured by adding a programmatically controlled perturbation signal to the process target depth. From the machining process perspective, this perturbation acts in a similar way as the robotic structure causing machined part imperfections. For instance, selecting a perturbation signal in a range of data up to ± 25% of the target depth would simulate the process behaviour in extreme conditions, where the robotic structure presents very low stiffness and positional accuracy.
Beyond the above perturbation signal that is introduced programmatically for greater data variability, there are other sources of errors in the process output variable. These include the clamping process, the panel movement during cutting, the measurement errors arising from the sensors and the inspection procedure, etc. The perturbation signal acts to amplify these errors and create a wide range of input conditions from which the systematic source of errors can be robustly modelled. In this study, the robots change their configuration every time a new hole of the panel is located, therefore a new perturbation value is used for every hole machined.
The target specified at the input of the overall manufacturing process (ref to Fig. 5) is the desired (nominal) value of the process depth of cut. The process model takes an input from the sensor signals collected during the cutting (i.e. the in-process response variables) and provides as an output a prediction of the inspection result (i.e. the post-process response variable). During the training phase the model builds a dataset with the acquired sensor data (signal features) and the corresponding inspection measurement (output variable). A prediction is then provided given the current dataset and it is compared with the true measurement value. The accuracy of the model is assessed in terms of the difference between the true and the predicted output value, called the prediction error. The model is considered trained when enough data has been collected for the prediction error to be within acceptable levels. A trained model is then ready to be applied for accurate prediction of the inspection measurement, as illustrated in Fig. 5b. This ability of indirectly measuring the performance of the process without any inspection step or process interruptions can potentially reduce the overall process operational cycle time, increasing productivity. Moreover, the additional information on the quality of the final product provided by the sensors, could be used to support the operator in making appropriate corrective actions should they be required.

Gaussian Process Regression (GPR)
In terms of the machine learning approach, this work exploits the use of Gaussian Process Regression (GPR) models due to their ability to deal with uncertainty in a probabilistic framework (Bayesian) and to provide confidence intervals of the predicted value. In addition, such models require less training data, can incorporate new knowledge (evidence) as more data becomes available and are less affected from the overfitting problem since the number of hyper-parameters to optimise during training is typically low.
A Gaussian Process (GP) generalises the Gaussian probability to describe a distribution over functions (instead of scalars or vectors) and allows therefore to consider inference directly in the function space [8].
) define the mapping between the input vector x and the scalar output y = f(x) + ε, with ε ∼ (0, σ 2 n ) being the measurement noise with a covariance σ 2 n . Note that unlike parametric approaches, GPR assumes f to be a random function that can be fully specified by its mean and covariance functions. The mean function m(x) is the expected value of f(x) at input vector x and the covariance (also called kernel) function k(x, x ′ ) is a positive-definite function that determines the covariance between pairs of input variables x and x ′ . Typically, when little or no prior information about the mapping f is available, a zero mean and a Squared Exponential (SE) covariance are assumed. The SE kernel is defined in Eq. (1): where σ f and l are called the signal standard deviation and the characteristic length-scale, respectively. These are both hyper-parameters that are optimised during training. The length-scale parameter represents a measure of how rapidly f(x) can change in the input space, i.e. how far two input vectors need to be, to become uncorrelated. Another type of the SE covariance, where each of the input dimensions has an individual length-scale parameter, is referred to as the Automatic Relevance Determination (ARD). The ARD-SE kernel is defined as: where d is the dimension of the input vector x and l i is the individual length-scale hyper-parameter for each input dimension x i . Note that Eq.
(2) implements ARD since the inverse of the length-scale determines the importance of the input dimension: a large value of l i , yields to a small covariance term for x i , meaning that this particular input dimension has little influence in the inference. This particular characteristic of the ARD-SE covariance function can be used therefore as an effective feature selection method, intrinsically implemented in the training process of the GPR model. This approach was used in the regression tests of this work, presented in Section 6.
In order to make predictions on new testing points X * , GPR considers the joint distribution of the observations y (i.e. the training targets) and the new predictions f * = f(X * ), which, according to the GP framework, is still a Gaussian with prior distribution given as: where K(X * , X), K(X, X * ) and K(X * , X * ) are the covariance matrices evaluated at all pairs of testing and training points, training and testing points and testing points only, respectively. The posterior distribution of the predictions is then calculated by conditioning the prior distribution of Eq. (3) to the training data (inputs X and the actual observations y) and it has the following mean and covariance: The two expressions of Eq. (4) represent the key predictive equations of the GPR, which respectively compute the conditional mean f * and covariance cov(f * ) of new testing points X * . Note that, while the mean f * gives the predicted output values of the new testing points, the covariance cov(f * ) gives additional information on the model reliability, such as the confidence intervals of predicted output values. The learning problem of a GPR model consists of finding the optimal values for the hyper-parameters specified into the covariance function. They can be inferred from the training data, using a maximum likelihood estimation approach. The log marginal likelihood is given by: where Θ is the vector of the hyper-parameters in the covariance matrix K. The optimal values are found by maximising Eq. (5) with respect to Θ and, since the number of parameters is usually small, it is not a complex optimisation problem.

Experimental setup and data description
In order to provide the necessary training data, an experiment was designed involving machining seven panels of 44 holes each for a total of 308 holes. Each cut was performed to a pre-set target depth plus an additional random perturbation (programmatically controlled) to obtain the data variability. A direct measurement of the true countersink depth was computed from the image provided by the inspection camera after each cut. This allowed collection of the actual output value without interruption of the process. The panels used during the experiments were Carbon Fibre Reinforced Polymers (CFRP) of dimensions 30mm x 21mm x 7mm with 44 pre-drilled holes distributed in a regular rectangular grid. An example of the panels before and after countersinking is illustrated in Fig. 6. Prior to the experiment, all the panels were painted in white colour in order to facilitate the machine vision information extraction.
During all the cuts conducted in this experiment, the same cutting tool was used, starting from an almost fresh state at Panel 1. At the end of Panel 4, the flutes of the cutting tool were visually inspected for any damages or excessive wear and then the tool was carefully cleaned up and re-installed for the last 3 panels. The tool wear was not constantly monitored and it was assumed to be negligible.
The random perturbation was introduced for depth variability (deep and shallow holes) and to allow the model to learn how small variations of the output variable (countersink depths) were related to the input data (sensor signals) collected during cutting. For this experiment, a Gaussian distributed variable with mean μ = 0 and standard deviation σ = 0. 25 (i.e ., N(μ, σ 2 )) was used as the perturbation signal.

Data acquisition system
The hardware used for the experiment was chosen for ease of installation with minimal intrusion to the machining process and hence to be suitable for a production environment. The following sensors were installed on the master robot end-effector: • 2x Single-axis High Frequency Accelerometers from PCB (type 621B40) • 1x Acoustic Emission (AE) sensor from Kistler (type 8152C01) • 1x Power Transducer from Caron Engineering (TMAC -Tool Monitoring Adaptive Control) • 1x Linear Encoder (digital probe) from Sony (type DK50PR5) The location of the sensors is illustrated in Fig. 7. The power sensor is not shown in the figure because the transducer was located inside the PLC cabinet. The three phases of power cable to the spindle were fed through the transducer and the power signal was directly wired to the data acquisition device. The vibration sensors and the AE sensor had magnetic clamps and they could be easily mounted on the metallic surfaces of the spindle holder and the pressure foot. The linear encoder was installed on the spindle holder to measure the linear displacement of the pressure foot during the clamping process.
The first vibration sensor (denoted by VIB1) was located as close as possible to the cutting tool with its sensitive axis being aligned to the feed direction (Z axis) in order to detect the tool vibration during cutting. The second (VIB2) was installed on the pressure foot body measuring the vibration in the perpendicular direction (X axis). The Kistler AE sensor was mounted on the pressure foot as close as possible to the cutting area to measure the acoustic emissions generated during the process. The Sony probe encoder was used as a trigger for the acquisition of the other sensor data. The movement of the pressure foot towards the panel during clamping was easily detected by the probe signal. This event was ideal for automatically triggering the acquisition since it was repeated at the start of every cut. All the sensor were connected to a central data acquisition device (compact DAQ) from National Instruments (NI cDAQ-9178) and the data was collected and processed  Fig. 7. Location of the sensors used for the experiment.
using NI LabVIEW software. The data acquisition system is shown in Fig. 8. The central acquisition unit made it possible to set a synchronised start for the sensor data collection and to simultaneously acquire several signals of different types and sample rates.

Signal processing and feature extraction
The purpose of the signal processing methods presented in this section was to convert the raw data into suitable features that could accurately represent the region the signal where the cutting process was taking place. After inspection of the acquired data from preliminary tests, it was noticed that the AE-rms signal (RMS computed at 100 samples per second) showed greater signal-to-noise ratio than the other signals. Therefore, AE-rms data was used for the implementation of a segmentation procedure, which automatically identifies the entry/exit points (time stamps) of the region of interest. This information was then used to segment all the other signals for further processing. An example of the method is shown in Fig. 9. The segmented part of the signal is the one in between the two vertical black lines.
The synchronisation of the signals due to the cDAQ device, ensured that the part of the signal identified by the segmentation procedure as the cutting region referred to the same event across all the sensors. Once the time information of the cutting entry/exit points was calculated, only the segmented parts of the raw signals were considered for further analysis.
The processing methods used to analyse the sensor data were specific to the signal type and included both analysis in time and frequency domains. In terms of the vibration data, in addition to the raw data (VIB) and the converted RMS data (VIB.RMS), three more signals were extracted: • The envelope signal (VIB.ENV) obtained by an envelope analysis of the raw data. • The converted RMS (calculated at 100 samples per second) of the residual signal (VIB.RES.RMS).
In this experiment, the spindle speed was set at 4978 RPM, which led to a TPF ≅ 165 Hz (for a 2-flute cutting tool) and a fundamental frequency component at F RPM ≅ 83 Hz (approximately), related to the spindle RPM. The envelope analysis was applied to detect the amplitude modulated signal around the TPF due to the cutting forces. It was performed on the raw vibration signal with central frequency at TPF (first order) and 40 Hz of bandwidth span. In order to obtain the residual signal, the fundamental frequency and all its harmonics (60 in total up to 5 kHz, including the TPF), were filtered out from the original data. This led to a new signal (VIB.RES) with similar amplitude, but shorter duration in time.
In terms of the power transducer data, in addition to the raw signal (PWR) and the converted RMS data (PWR.RMS), two more signals were extracted: • The Low-Pass Filtered signal (PWR.LPF) obtained from a low-pass RC circuit (cut-off frequency of F c = 72.4 Hz).
• The Time-Frequency Filtered signal (PWR.TFF) computed by applying a Gabor transform with all coefficients below a certain threshold set to zero, leaving just those that carried the signal energy. The time domain signal was then reconstructed from the remaining coefficients. The parameters used for the Gabor transform were the following: threshold value 0.87, time steps 8, frequency bins 512, window length 512 (Gaussian).
In terms of the AE data, in addition to the two original signals provided by the AE coupler, the raw AE and the converted RMS signal (AE. RMS), two more were considered for further analysis: skewness, kurtosis, peak and crest factor. In addition, some other features potentially related to the countersink depth such as the cutting cycle duration (duration of the segmented part of the signals), the Timeto-Peak (time spent to reach the peak amplitude) and the duration of the VIB.RES.RMS signal were also included.
In order to assess the frequency content of the signals, the Power Spectral Density (PSD) was calculated for each of the raw sensor data in the cutting region, identified by the segmentation method. The features were extracted from two main spectra: the complete power spectrum of the raw data (this includes also the VIB residual signal) and bands of the power spectrum obtained from the original PSDs. These bands were chosen based on observations of the signal's original PSD and they were specific to each sensor: Vibration, [2k-3.5kHz]; AE, Band1 [15k-25kHz], Band2 [80k-160kHz] and Band3 [80k-160kHz] (Filtered); Spindle Power, . The frequency domain features included statistical descriptors relative to bands of power spectrum, such as the spectral mean, variance, skewness, kurtosis and peak amplitude. The sum of total band power, the relative spectral peak per band and the total harmonic band power (only for the vibration data) were also calculated. Two more features (independent from the sensor signals) were added: (i) Tool Hole Counter (THC), which counted the number of holes that the tool had machined so far, starting at one at first hole of P1, and (ii) Hole Position (HPos), which was the normalised value of the current hole number within the 44-holes of the entire panel.
A summary of all the signals and all the features considered for the experiment is presented in Tables 1 and 2, respectively. Note that the time domain signals were categorised in periodic (oscillating signals such as those of VIB and AE data), non-periodic (PWR data) and the converted RMS signals. The total number of the extracted features (including those independent of the signals) was 230.

Depth of cut data distribution
The dataset contained samples of the feature vector extracted from    Tool Hole Counter  THC  -----Hole Position HPos the cutting cycle together with the measured countersink depths obtained from the inspection camera. The programmatically added Gaussian perturbation allowed to obtain examples of cuts with different depths, distributed (normally) around the nominal target (d), set at 2 mm in the tests. Fig. 10 shows the depth deviations (i.e. the depth errors calculated as the difference between the camera measurement and the nominal target). Negative values represent under-target depths with respect to the nominal (shallow holes) and positive values refer to overtarget depths (deep holes).
The Gaussian perturbation added to the depth target had mean μ = 0 and standard deviation σ = 0.25, therefore, most of its values belonged to the interval [-0.5 0.5]. However, the measured values were distributed in the range [-0.3 0.7] mm and thus showed a positive bias of approximately 0.27 mm, with the majority of the holes having overtarget depths. Fig. 11 shows the effect of the Gaussian perturbation on the observed depth deviations: the left chart gives the correlation between the perturbation signal and the depth deviations obtained after the cuts; the right chart provides the histogram distribution of the depth deviations.

Regression results and discussion
The regression tests conducted in this work investigated the use of Gaussian Process Regression (GPR) as predictors of the countersink depth data from the extracted signal features. This section presents and discusses the results in terms of the feature selection strategy, the model creation and validation methods and the assessment of the model confidence.

Feature selection
The feature selection strategy consisted of an embedded approach integrated with the learning algorithm of the model. The method applied GPR to perform Automatic Relevance Determination (ARD). This was obtained by selecting the ARD-SE kernel defined in Eq. (2) as the covariance function for the GPR model. The choice of this particular kernel automatically assigns weights to the input dimensions (signal features) according to their relevance, in the optimisation procedure during training. This means that at the end of the training phase, the features with higher weights would have greater influence in the model's output variable than those with lower weights. Therefore, setting a threshold on the weight value would effectively reduce the original feature space by removing all the irrelevant features that fall under such threshold.
In the tests conducted here, the initial values for the hyper-parameters of the GPR model (relative to the GP prior distribution) were determined from the training data according to the GP framework described in Section 4.2. The parameters of the GPR model with ARD-SE kernel were set as follows: • the initial signal standard deviation σ f = std(y) ̅̅ 2 √ • the initial individual length-scale for each feature l i = ̅̅̅ d √ • the initial noise standard deviation σ n = std(y) where std(⋅) = ̅̅̅̅̅̅̅̅̅̅̅̅̅ Var(⋅) √ refers to the standard deviation of the distribution function of the output observations y and d is the number of features in the dataset (i.e. the input dimension).
Note that the feature weights (computed from the length-scale parameter l i ) are optimised during the training of the GPR model, and therefore will depend on the particular training set. This means that using training sets of different sizes and/or containing different data samples will result in different weight values for the features. Hence, a strategy was required in order to obtain good generalisation properties. A common option is the cross validation approach, which divides the original dataset into training and testing data. In order to keep the training set as large as possible, in the tests here, the Leave-One-Out Cross Validation (LOO-CV) method was applied. All the available data but one sample was used to train a GPR model with ARD-SE kernel, which was then tested on the left out sample. The process was repeated for all the samples in the dataset, building a new model in each iteration and recording the feature weights assigned during training. At the end of the procedure, the weight of each feature was calculated as the average over the values assigned to that feature by all the trained models. The features were then sorted based on their weight and only those with weights over the pre-set threshold value were selected for further analysis. This method was computationally expensive and the time to train all the models increased considerably with the size of the dataset. However, for the purposes of the tests here (offline analysis), it provided more robustness to the assessment of feature relevance and the generalisation performance of the GPR models.
The overall feature weights assigned by the LOO-CV method are shown in Fig. 12 with the specific features that obtained higher weight values distinctly identified. With reference to the feature names presented in Table 2, the model assigned higher weights to features extracted from VIB signals, such as VIB2.RES.Ku and VIB1.ENV.RMS, followed by PWR.LPF.P and AE.BP1.STBP.
All the selected features obtained at the end of the training procedure are summarised in Table 3 (normalised) weights of the features was set at 0.001. The number of signal features selected from the method was 19.
Finally, Fig. 13 presents the comparison results in terms of the prediction errors (i.e. the difference between true and predicted values) obtained from the LOO-CV procedure using two different feature sets: one training GPR models with all extracted features and the other training GPR models with the reduced feature subset (i.e. using the selected features only). Note that the results are shown as box plots of the prediction error range comparing the two cases.    Fig. 13. Box plots of the prediction errors comparing two cases: GPR models trained using all the features vs. GPR models trained using the selected features only.
The model showed an improvement in the prediction accuracy when using only the selected features. The model Loss value (i.e. the MSE of the predictions) decreased from 0.005 (all features) to 0.003 (selected features), In addition, the overall LOO-CV procedure of training the GPR models (SE kernel) with the reduced feature subset took considerably less time compared to the case when using all the features.

Depth prediction results: GPR-SE model
The feature selection method proposed in the previous section proved to be effective in the identification of relevant signal features with respect to the model's output variable. This achieved higher prediction accuracies and faster training times. Therefore, the regression tests discussed in this section will consider LOO-CV procedure and GPR models with SE kernel using only the selected features summarised in Table 3. To validate the model, the following test was performed. First, 44 data points (equivalent to one panel of data) were randomly chosen from the entire dataset and set apart for testing. Then, a GPR-SE model was trained on the remaining dataset (equivalent to six panels of data). The model used the reduced feature subset with standardised values and the hyper-parameters were initialised based on the training set, i.e. the original dataset excluding the extracted 44 testing points. After the training phase was finished, the model provided an expected predicted output (with a mean value and its corresponding variance) for each testing data. This extra information of the variance was used to put confidence intervals to the output prediction.
An example of the results obtained by the test is illustrated in Fig. 15. The shaded region represents the 95% CI (Confidence Interval) of the prediction (mean) value indicated by the blue crosses. The true observations (i.e. the depth deviations) are indicated by the red crosses. The 95% CI was calculated based on the GPR variance prediction. The majority of testing points was located inside the shaded region meaning that the model accurately identified the 95% CI.
Note that the predicted variance provides an indication of the model's confidence in the output value, given the current dataset used during training. This gives a significant advantage to the GPR models (and the probabilistic framework in general) over other learning techniques, when dealing with uncertain data.

Model confidence assessment
The model confidence in the output prediction was assessed in terms of the predicted variance. This section analyses how the predictive variance of the GPR model is affected by the size of the dataset and the quality of new input data.
To assess how the model confidence is affected by the choice of the training data, the following test was conducted. In Case 1, a GPR-SE model was trained using the input data of the first two panels P1 and P2 and then, predictions of the remaining instances (relative to the panels P3-P7) were obtained. The test was performed on the dataset described in Section 4.3 and the results obtained are shown in Fig. 16. In Case 2, the number of training points was increased to the first four panels and four holes of data in Panel 5, and repeating the above test yielded to the results illustrated in Fig. 17. Note that the hole numbers are indicated cumulatively across the remaining panels.
From the comparison of the two cases, the following observations can be made:  • The inclusion of the additional training data in Case 2 had considerably improved the performance of the models with respect to the predictions on the last three panels.
• Similarly, the shaded region of 95% CI showed a reduction in size, indicating that the information gained from additional data helped to lower the uncertainty levels for the prediction of the remaining data of P5, P6 and P7.
These observations can be understood from the process by which the GPR model calculates the predictions. A GPR model computes the covariance matrix of all pairs of training and testing inputs, considering a joint Gaussian distribution. In particular, the calculation of the predictive covariance of the testing inputs depends only on the covariance matrices (of both training and testing input data), not the output observations. Therefore, the choice of the covariance function used for the computation of these matrices will affect the predictive covariance of the GPR model. The SE kernel function adopted in this study computes the covariance of a pair of input instances (i.e. feature vectors) in terms of their distance (Euclidean) in the d-dimensional feature space. Thus, inputs that are close to each other in the feature space will be highly correlated, while those far away less correlated. This means that new testing inputs far away from the multivariate Gaussian distribution defined by the training data will have higher variances, indicating that the model is not confident in the predicted outputs. The inclusion of additional data in Case 2 appears to contain input information that was not presented in Case 1, but was prevalent in panels P6 and P7. This resulted in the improved prediction and reduced CI.
In order to assess the amount of training data that would be appropriate to obtain a GPR model with required predictive accuracy with tight CI, the following test was performed. A subset of 44 randomly selected instances (one panel equivalent) from the entire dataset was left out as a validation set. Then, a series of GPR-SE models were trained on the remaining data, gradually incrementing the training data size, starting from 22 points (half panel equivalent) up to the maximum  These results indicate that both the model accuracy and the confidence in the predicted values increase as the size of the dataset increases. This can be used to set a minimum size of the training dataset that would be sufficient to obtain accurate prediction with the required confidence. For the investigations here, this minimum size would be in the range of 120 to 160 data points, which is about three to four panels of machining.

Conclusions
This paper presented a novel perturbation signal based approach for generating high quality training data suitable for building a robust process model for part quality prediction. A probabilistic learning method based on Gaussian Process Regression (GPR) was used to predict process-specific quality indicators using exclusively the sensory information acquired during a robotic machining process. The experimental results validated the ability of the proposed GPR model to accurately map the extracted signal features (input variables) with the process countersink depth deviations (output variable) obtained post-process from the inspection camera. The use of a probabilistic learning approach combined with the Automatic Relevance Determination (ARD) based feature selection allowed a robust model to be built with a small dataset and high dimensional feature set. The feature selection strategy was shown to give an improvement of the model prediction performance. The complexity of the robotic manufacturing processes and the noise in the sensors introduced several sources of errors in the process input-output relationship. In order to account for the uncertainty arising from the sources of errors and the size of the dataset, the model confidence was also assessed in terms of the predicted variance. The GPR model was chosen for its ability to provide the variance as an additional output to the mean prediction. Guidance on the minimum dataset size for accurate predictions with the required confidence was also given for the considered case study. The potential for data-driven machine learning approaches to enable reduced inspection costs in robotic manufacturing processes has been demonstrated in this paper.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.