Accurate Prediction Of Machining Feedrate And Cycle Time Prediction Considering Interpolator Dynamics

This paper presents an accurate machining feedrate prediction technique by modeling the trajectory generation behaviour of modern CNC machine tools. Typically, CAM systems simulate machines motion based on the commanded feedrate and the path geometry. Such approach does not consider the feed planning and interpolation strategy of the machines numerical control (NC) system. In this study, trajectory generation behaviour of the NC system is modelled and accurate cycle time prediction for complex machining toolpaths is realized. NC systems linear interpolation dynamics and commanded axis kinematic profiles are predicted by using Finite Impulse Response (FIR) based low-pass filters. The corner blending behaviour during nonstop interpolation of linear segments is modeled, and for the first time, the minimum cornering feedrate, that satisfies both the tolerance and machining constraints, has been calculated analytically for 3 axis toolpaths of any geometry. The proposed method is applied to 4 different case studies including complex machining tool-paths. Experimental validations show actual cycle times can be estimated with greater than 90 percent accuracy, greatly outperforming CAM-based predictions. It is expected that the proposed approach will help improve the accuracy of virtual machining models and support businesses decision making when costing machining processes.


Introduction
With the introduction of concepts like virtual manufacturing [1] and digital twins [2], building process models and predicting actual machining process conditions in the computer environment has become paramount in attaining higher productivity and throughput in today's manufacturing. For example, accurate machining cycle time prediction is vital for industry during the quotation process to ensure achievable and profitable contracts. The prediction models and generation of accurate digital twins is a collective modeling effort which requires both detailed modelling of the process as well as the dynamic machine behaviour. Considering the machining processes, current literature provides accurate models to predict milling process physics [3,4]. Nevertheless, when applied in practice, these models show large discrepancies from the actual process behaviour.
One reason can be identified as the influence of the machine tool drive dynamics. In particular, the behaviour of the Numerical Control (NC) plays a key role. Trajectory generation (interpolation) algorithms embedded in the NC system, control the feedrate profile, which is a key input for machining process models. For example, contouring (positioning) errors alter tool engagements [5] which lead to inaccurate force predictions [6]. Thus, in order to accurately develop realistic digital twins for machining processes, the feedrate profile generated by the NC system of a machine tool must be accurately predicted. This paper deals with modeling and prediction of interpolator dynamics of modern NC systems to accurately estimate machining cycle times and cutting forces along complex parts.
Once a part program (G-code) is deployed to a CNC machine tool, the NC unit parses the part program and interpolates the tool motion between successive cutter locations (CL). Most modern CAM systems provide toolpaths in terms of discrete CL-data and rely on linear interpolation algorithms that run in the NC units. With the introduction of cheap memory modules, long part programs do not pose a limit, and even basic circular paths are programmed with series of short linear segments [7,8]. Therefore, modern NC systems are equipped with propriety algorithms that interpolate these lengthy series of short CL-blocks smoothly. These algorithms are called Look-ahead or Compressor functions and are capable of generating a non-stop motion with time optimal feed-rate profile [9] that respects kinematic limits of the machine [10,11]. Prediction of a machine's actual feedrate profile requires detailed modeling of the NC system's real-time interpolation behaviour. This includes the motion transition between CL-blocks, for example a typical feedrate profile for continuous motion is shown in Fig. 1. During the initial linear motion from zero to commanded feedrate the performance and behavior of the machine tool is dependent upon the acceleration and jerk constraints alone. However, as the tool approaches the end of the first CL-line (corner transition 1 in Fig.1) to change the feed direction the tool decelerates to a minimum cornering feedrate before accelerating again to the commanded feedrate. The reduction in feedrate in the vicinity of CL-line junction point is due to both the machine tool satisfying the tool centre point (TCP) error tolerance constraints throughout the cornering transition and the machine tool kinematic constraints. The TCP error can be seen at corner transition 2 where the TCP is maximum displacement between the CL-line and the TCP position. The TCP error constraint imposed upon the toolpath limits the maximum feedrate during cornering transitions and this significantly affects the overall machining cycle time.
Most NC systems utilize jerk limited trajectory generation to smoothly alter feedrate and interpolate along CL-lines [10,12]. The generated feedrate profile is defined in the form of a cubic polynomial [7]. Axis acceleration limits are imposed based on the torque/power capacity of the drives, and the jerk limits are set to limit unwanted vibrations during rapid feed motion [13]. This general jerk-limited feedrate profile is well-known, and acceleration and jerk limits of the machine can be read from the NC system. Therefore, the use of jerk limited trajectory as a template allows prediction of feedrate kinematics of modern NC systems. Nevertheless, it can only predict machine behaviour in point-to-point (P2P) interpolation. During P2P interpolation, the tool accelerates from a full-stop to the set feedrate and decelerates again for a full-stop at the end of the CL line. Once the acceleration and jerk limits are known, the feedrate profile can be generated to predict cycle times. Past research considered modeling of NC behaviour of 3 and 5 axis machine tools for P2P trajectories [14,15].
Predicting feedrate profiles along short segmented complex toolpaths for high speed machining (HSM) is, however, a more challenging task. This is due to the fact that look-ahead modules of NC systems alter jerk limits on the fly as it blends series of CL-lines to generate a non-stop smooth continuous feed motion. Here, modeling the path blending behaviour is crucial. NC systems blend linear CL-lines together smoothly while applying geometric blending error and kinematic limit control. Machine tool literature reports that circular arcs [7] , cubic [7] or quintic splines [10] can be used for such geometric path blending. There are also methods based on filtering where the discrete toolpath is blended based on low-pass filtering. Finite Impulse Response (FIR) filters are used for such purpose [16]. Such filtering based techniques are more computationally efficient and greatly favored for real-time interpolation on NC systems. For instance Heidenhain [17,18], Mitsubishi [19] and more recently Siemens [20] NC systems utilize FIR and IIR (infinite impulse response) filters for look-ahead and non-stop smooth interpolation. Typically, users enter a blending tolerance which confines the path blending (contour) errors. Based on the blending tolerance the NC system approximates the given discrete CLlines and plans the fastest motion with its kinematic limits. Therefore, accurate prediction of cycle times for conventional toolpaths requires modeling of NC system's non-stop interpolation behaviour along linear paths.
This paper models the non-stop interpolation behaviour of modern NC systems and predicts feedrate profiles along HSM toolpaths by considering the real-time path blending behaviour of NC systems. Section 2 briefly introduces the low-pass filtering based real-time interpolation method, which is used as a template. It is then used to predict P2P and contouring motion of NC systems in subsequent sections 3 and 4. Illustrative examples and experimental validations are provided in each section. Finally, Section 5 provides realistic cycle time, feedrate profile and cutting force prediction for complex aerospace parts.

Low-pass Filtering Based Real-Time Interpolator Dynamics
This section models real-time interpolation behaviour of an NC system to predict the feedrate profile and overall cycle time. Most conventional NC systems utilize IIR or FIR filtering based techniques for computationally efficient real-time interpolation and feed profile planning. In this work, Finite Impulse Response (FIR) filters are used to capture the NC system's behaviour. A simple 1st order FIR filter can be expressed in the Laplace (s) domain by: Fig. 2 Impulse response of a 1st order FIR filter.
where s is a complex number, T i is the time constant of the i th filter. The impulse response is depicted in Fig. 2. As seen in (1), the filter contains an integrator, which acts to smooth the input signal. These two features of 1st order FIR filters are appealing from a NC system perspective, since G-codes (represented by rectangular velocity pulses) can be convolved through a series of such filters to generate smooth velocity profiles. Since the area underneath the rectangular impulse response is unitary, the area underneath the original input is not altered. [16,17,19,20,21,22]. Fig.3 illustrates this filtering based interpolation procedure. As shown, consider a G-code for a total displacement command of L at a feedrate of F . It is represented by a velocity pulse with an amplitude of F and duration of T v hence L = F T v . Subsequent convolution of the velocity pulse with the FIR filter yields the higher order velocity response. Using 2-FIR filters in series generates reference trajectories with piece-wise constant jerk profiles and using three FIR filters in series further smooths the reference velocity making them snap limited. Although jerk-limited trajectories are most common in high speed machinery, snap limited trajectories are tuned for ultra-precision machines [23] to further mitigate the effect of unwanted vibrations.
The duration of the original velocity pulse T v and the time constants of the filters T 1 , T 2 and T 3 determine the velocity and acceleration profiles, which can be derived analytically by evaluating the convolution integral between the input velocity pulse and the rectangular impulse response of the filter as follows: where v(t),v (t) and m(t) represent the velocity pulse, interpolated velocity signal and the impulse response of the FIR filter (Eq.(1) respectively.
For T 2 > T 1 > T v the velocity and acceleration profiles become: Smooth trajectory generation by low order FIR filtering [16] v (t) = It is shown in Fig.3a that when a square velocity pulse is convolved with a first order FIR filter the result is a trapezoidal velocity profile. In this case the time constant of the FIR filter T 1 is smaller than the length of the velocity pulse T v , (T 1 < T v ). It can be seen that the resultant length of the velocity profile is extended to T v + T 1 , the length of the filtered velocity profile is elongated by the length of the FIR filter. The figure also shows the profile is segmented into 3 main kinematic sections. The figure shows that when T 1 < T v the commanded feedrate F is reached and maintained for the cruise duration T v − T 1 . The filtered acceleration profile shows the peak acceleration is V /T 1 , demonstrating the FIR filter time constant governs the acceleration properties of the filtered kinematic profile. When T 1 > T v then the commanded velocity cannot be reached and the maximum velocity is determined by L/T 1 .
Cascading FIR filters can increase the order (smoothness) of the kinematic profile as is shown in Fig.3b. Using an FIR filter with a time constant of T 2 , where T 2 < T 1 , the trapezoidal velocity profile is filtered once more. The resultant motion profile is constructed from 7 segments in which the time of each segment is determined from T 1 , T 2 and T v . It can be seen that the total duration is T 1 + T 2 + T v . Due to the linearity of the convolution operation [24] the final kinematic profiles are independent of the order of the filtering, it is the relationship between the magnitudes of the time constants that determines the maximum feedrate and acceleration [16]. The analytical expressions of velocity and acceleration for each segment are shown in equations (3) and (4) respectively.
Figs. 4a and 4b show example profiles for different velocity pulse and filter delay parameters. As shown, the maximum acceleration and jerk values of the profiles are determined by the time constants of the FIR filter. In Fig.4a the commanded feedrate is reached on the precondition that T v > T 1 + T 2 and T 1 > T 2 and in this case the peak acceleration can be determined from the longest FIR filter time constant V /T 1 . However, as shown in Fig.4b where T v is smaller than the smallest filter delay T 1 > T 2 > T v , then the commanded feedrate F is not reached. The maximum feedrate and acceleration is determined by the commanded feedrate and the relationship of the filter time constants T 1 ,T 2 and the velocity pulse T v . The maximum feedrate and acceleration are calculated from F T v /T 1 and F T v /T 1 T 2 respectively. Note the maximum feedrate is limited by the longest time constant and velocity pulse length T v , whereas the acceleration is constrained by both time constants.
Time constants of the filters are typically selected to mitigate structural vibrations on the machine [25]. Matching the time constant with the vibration period of the lightly damped modes helps avoid exciting them during rapid acceleration. Another way to specify the time constants is to set them equal, T 1 = T 2 . In this special case, the FIR filter acts as a pure low pass filter with a roll-over frequency of ω c = 2π T1 . Fig.5 shows the attenuation in the frequency response for multiple FIR filters with matching time constants. The time constant, when set low enough, helps prevent the excitation of any higher frequency vibrations during rapid accelerations. This simpler method compared to tuning individual filters provides a convenient method of vibration suppression during high feedrates.
For interpolation using 2-FIR filters with matching time constants, the transfer function of the resulting FIR filter is: and the resulting velocity profile when a rectangular feed pulse v(t) is filtered The corresponding acceleration and jerk responses are : As demonstrated, when a square velocity pulse of magnitude F and length T v is convolved with a first order FIR filter with time constant T 1 the result is a trapezoidal velocity profile with constant acceleration of magnitude F/T 1 (Fig. 6a). The total length of the kinematic profiles are extended by the filter time constant T 1 to T v + T 1 . When the trapezoidal velocity profile is convolved with a second first order FIR filter with a matching time constant T 1 = T 2 the smoothness (order) of the velocity profile is increased from C 1 to C 2 , where C n is the space of n th order continuously differentiable functions, as shown in equations (7) and Fig. 6b. However, using the matching time constant T 1 = T 2 , results in five sections in the kinematic profile and not seven as for the case for two different time constants where T 1 = T 2 . The resulting acceleration profile is triangular around T 1 and T v + T 1 with peak magnitudes F/T 1 and lengths of 2T 1 ; the now jerk limited profile has peak magnitudes of F/T 2 1 . The total length of the kinematic profiles is extended to T v + 2T 1 . The relationship between T 1 and T v determines the kinematic constraints as for the different filter cases. For completeness, the authors include the example profiles for the matching filter case in appendix A for different velocity pulse and filter delay parameters.
Convolving the velocity profile with a third first order FIR filter with the same time constant T 1 = T 2 = T 3 results in a C 3 velocity profile, C 2 acceleration profile and C 1 jerk profile. The velocity, acceleration and jerk equations for the 3-FIR case is shown in appendix B. The smooth acceleration profile has a peak magnitude of 3F/4T 1 at times 1.5T 1 and T v + 1.5T 1 and the jerk  profile has peak magnitudes of F/T 2 1 . The overall length of the kinematic profiles have been extended from the original square velocity pulse length T v to T v + 3T 1 . The total filter delay when using 3-FIR filters with matching time constants T 1 is therefore 3T 1 .
It can be shown that a high order FIR filter can be accurately modelled and implemented with using only 3 first order FIR filters. The benefit of using 3 or more first order FIR filters with the same time constant is that the filter response approaches that from a Gaussian filter. The Gaussian response has no overshoot whilst minimising the acceleration and deceleration time periods which makes it the ideal time domain filter for interpolating kinematic profiles [26]. The ability to approximate the Gaussian filter with 3 FIR filters with the same time constant simplifies the design and selection of the filter to a single design parameter T 1 . For both the 2 and 3 FIR filter cases, T 1 can be analytically calculated from the maximum permissible jerk J max as follows:

Identification of Real-Time Interpolator Dynamics of an NC system
The previous section presented the filtering based real-time trajectory generation. In this section it is shown how the interpolator response of a machine The machine is commanded by a single G-code to move 6 mm at a speed of 3000 mm/min, and the interpolated reference motion profile is recorded on the NC system directly at a sampling time of T s = 0.009s. Figs 7a to 7b, show the recorded kinematic profiles. The machine is set to undergo a simple point-to-point (P2P) motion and therefore the tool comes to a full stop before moving to the next commanded position. As shown for the measured system, the NC system generates smooth velocity and acceleration profiles. The acceleration profile mimics a smooth 'bell-shaped' profile. Overall, acceleration, and deceleration duration are measured to be T acc = T dec = 0.0765 sec. The cruise velocity portion is roughly measured to be 0.023 sec. In order to simulate the feed profile, a series of 2 and 3-FIR filters are used. For the 2-FIR case the time constant is selected as T 1 = Tacc 2 and for the 3-FIR case it is set to T 1 = Tacc 3 . The predicted velocity and acceleration profiles for the 2-FIR case are shown in Figs 7a and 7b respectively. The time of the measured displacement is equal to the time of the predicted displacement. The difference between the velocity profiles is due to the acceleration. The 2-FIR case exhibits the triangular acceleration profile compared to the smooth measured response. The maximum acceleration for the 2-FIR case is constrained and less than the measured response.
In order to compare the different filter cases the machine is commanded to move along the same G-code, and the proposed interpolator model for the 3-FIR case is used. As shown in Fig.7c the velocity profiles for the 3-FIR case closely resembles the measured velocity profile and the total time of the measured displacement matches the total time for the simulated displacement. The simulated acceleration profile is smooth and the maximum acceleration is higher than for the 2-FIR case but still lower than the measured response.
Increasing the order of the simulated system would allow the maximum acceleration to approach the measured response. In general, by increasing the order of the FIR filter, the predicted acceleration profile of the filtered pulse approaches the acceleration profile of the measured response and results in a simulated velocity profile which closely resembles the dynamics of the machine interpolator.
The filter delay is calculated from the jerk (10) and the duration of the acceleration phase in each case is equal to the total filter delay. The time constant (filter delay) can be analytically calculated from machine tools' specifications (J max ) and therefore kinematic profiles can be generated using FIR filters without the requirement for parameter identification through system testing.
In this section it has been shown that the dynamics of an NC interpolator are increasingly well-approximated by the series combination of identical firstorder FIR filters. In addition, the relationship between the parameters of these first-order filters and the resulting interpolator response have been derived.

Multi-Axis P2P Motion Generation
FIR filtering based interpolation of single axis motion was presented in the previous sections. Extending the method to P2P multi-axis linear motion this section describes the process to interpolate kinematic profiles between two points using high order FIR filters.
The start and end positions of a linear G01 command in 3 axes can be represented by P s = [P s,x , P s,y , P s,z ] T and P e = [P e,x , P e,y , P e,z ] T , respectively as shown in Fig. 8a.
The tool displacement L is calculated by taking the Euclidean norm of the vector between the two commanded positions, L = P e − P s 2 .The velocity pulses of each axis (v x , v y , v z ) are calculated by multiplying the feed pulse v(t) by the unit velocity vector u = (P e − P s )/ P e − P s 2 .
whereṖ(t) represents the first time derivative of the P2P displacement (Fig. 8b).
In order to generate (and interpolate) the reference velocity commands (v x ,v y ,v z ), the individual axis velocity pulses (v x , v y , v z ) are convolved with the FIR filter (Figs. 8c and 8d): Finally, the filtered position commands are generated by integrating the filtered axis velocity commands:

Prediction of Interpolator Behaviour during Non-stop High Speed Motion
The previous section showed that P2P linear interpolation behaviour of an NC system can be modelled by velocity pulses low pass filtered by a series of first order FIR filters. The only required parameter to predict the machine's feed profile and accurately estimate the resulting cycle time is the time con-stant, i.e. total delay of the FIR filter. As shown, the filter time delay can be calculated from the maximum permissible jerk (10) and commanded feedrate. This section focuses on accurate prediction of interpolator behaviour during non-stop contouring motion, which is the most commonly used interpolation technique for high speed machining (HSM).

Modeling of Non-stop Interpolation Behaviour
Typical high speed machining toolpaths found in die and mould manufacturing or in aerospace industry consist of series of short segmented toolpaths [27]. When interpolated in HSM mode, the NC interpolator does not undergo a fullstop at the end of each CL line. Instead, the CL lines are blended together for a non-stop smooth motion interpolation where machining feedrate is reduced to a cornering speed V c around junction points of the CL-blocks (See Fig.1). The prediction of V c is crucial to accurately capture the actual feedrate profile and estimate the resultant cycle time. Several constraints affect the cornering speed (V c ) and overall acceleration profile around the CL data points. Firstly, V c is controlled by the blending (cornering) tolerance [14]. Typically, lower blending tolerance delivers more accurate motion but generates slower feed profiles.
In contrary, a larger tolerance value allows faster speeds and shorter overall cycle time. The relationship between the blending tolerance and the feed drop around the corner must be captured. Secondly, the deceleration/acceleration profile and the transition duration from the programmed feedrate (F ) to the cornering speed (V c ) are dictated by acceleration and jerk limits of the machine. Both of these key characteristics must be modelled to accurately predict the varying feedrate profile along HSM tool-paths.
In an effort to accurately model the interpolator behaviour, the feed pulse distribution shown in Fig. 9b is proposed in this manuscript. Notice that the feed pulse profile is different from the case used for the P2P motion. Feed pulses of each CL block are commanded back-to-back with no dwell time in between. In other words, they are constructed as a continuous pulse stream. The duration of the feed pulse is T v . Notice that the feed pulse does not have a constant amplitude of F . Instead, around CL block junctions the feed command value is dropped down to F c . Such small feed pulse is added to model the blending kinematics, commanding the feedrate to drop down to a cornering feed of F c . The duration of the cornering feed pulse is set to T b , which controls how long the deceleration and acceleration last around the blend.
When the feed pulse profile is interpolated with a FIR filter the resulting velocity profiles are smooth velocity profiles that better approximate the actual velocity profiles of the machining interpolator. Fig.9a and Fig.9c show the toolpath and the corresponding interpolated X-axis and Y-axis velocity profiles respectively. The total length of the velocity profiles is equal to the sum of the pulse lengths plus the filter delay T d . Fig.9d shows the cornering feedrate V c of the resultant velocity profile is equal to the commanded blending pulse feedrate F c and this occurs at half the filter delay T d /2 from the start of the Y-axis profile.
The cornering feedrate is controlled by setting the blending velocity pulse F c equal to the desired cornering tangential velocity V c and setting the acceleration and deceleration time for the interpolated feed profile equal to the time required to reduce from F to F c . A scaling factor is applied to F to represent F c as a function of commanded feedrate F : where V c is the resultant 3-axis TCP velocity defined as and v x , v y and v z represent the interpolated axis velocities at the minimum corning feedrate.
The total acceleration and deceleration time of the interpolated feed profile to reach F α from F is represented by T b , it is a function of the filter delay T d , and it can be calculated as: The final component to the pulse train is determining the main velocity pulse lengths T v . In section 2 the length of the velocity pulse T v was calculated from L/F , however, with the introduction of the blending pulses, T v must be modified in order to preserve the total area of the pulses and hence the TCP displacement.
The commanded TCP displacement is calculated from the total area of the velocity pulse and the blending pulse, this can be seen in Fig.9b where the total area within the X-axis and Y-axis pulses is equal to L 1 and L 2 respectively.
For a single axis displacement L the pulse areas comprise of the main pulse (calculated as F T v ) and the blending pulse (calculated as F c T b ): Rearranging equation (17) and incorporating equation (14) yields the modified value of T v as: Equation 18 holds for velocity commands with a single blending pulse, this is the case for the initial and final CL lines in a part program which start and end at zero feedrate (full stop). The remaining displacements in a part program are continuous and therefore the commands consist of a velocity pulse with a blending pulse either side as shown in Fig.10. Therefore each cornering blend consists of two back to back blending pulses.
For the entire pulse train, each G01 command or CL-line can be represented by an index k with k=1 corresponding to the initial command in the part program. The associated feedrate commands in the part program are hence denoted F (k). Therefore, for the main commands in a part program the modified value of T v is calculated as: For constant feedrate the adjoining blending pulses are symmetric. This leads to symmetrical interpolated velocity profiles and results in symmetrical displacement profiles, translating to the same toolpath trajectory for both forward and backward passes resulting in a more accurate finish.

Filtered Signal Generation
The composition of the velocity pulses and filtered kinematic profiles was shown in the previous section. In practise, the strategy for interpolation of multi-segmented NC tool-paths using high order FIR filtering, as shown in Fig. 10, is as follows: Sections 4.1 and 4.2 described the components of the velocity pulse train and application of FIR filtering for generation of kinematic profiles for nonstop high speed motion. The following sections will analytically demonstrate the relationship between the cornering speed V c to the blending error and axis kinematic limits and ultimately demonstrate how F c is selected to guarantee these constraints are satisfied.

Kinematic Profiles for the 2 First Order Filter Case
The geometry of velocity blending pulses was presented and calculated in section 4.1. The pulse signals are interpolated using FIR filters to generate kinematic profiles that control the cornering feedrate. This section analytically derives the equations for the kinematic profiles when using velocity blending pulses and FIR filtering based interpolation to control the cornering feedrate. In doing so, the authors are able to analytically calculate the blending pulse feedrate command F c which satisfies both TCP error and machine kinematic constraints during cornering transitions.
Using 2-FIR filters with matching time constants to interpolate a velocity pulse signal results in the kinematic profiles shown in figure 11. The profiles are split into 5 sections during acceleration/deceleration as shown in Fig.11b  for the Y-axis acceleration. The objective of the analytical expressions is to calculate the interpolated displacement at the point of maximum TCP error and the interpolated velocity at the minimum cornering feedrate. This occurs at half the total filter delay T d /2 (see Fig.11a). The total filter delay for the 2-FIR case is T d = 2T 1 , where T 1 is calculated from the maximum permissible jerk (equation (10)), therefore the maximum TCP error and minimum cornering feedrate occurs at T 1 . Fig.11b shows T 1 is at the start of section 3, therefore only sections 1-3 of the kinematic profiles need considering. The analytical expressions for sections 1-3 of the displacement, velocity, acceleration and jerk profiles for the 2-FIR case are presented in equations (20) to (23) respectively.
The interpolated axis velocity at maximum TCP error (minimum cornering feedrate) occurs at t = T d /2 = T 1 , therefore in the 2-FIR filter case this results in the following expressions for interpolated velocity (24) and displacement (25): Using equation (16), the interpolated displacement (24) and velocity(25) can be expressed in terms of F and α: Fig.12 shows a cornering transition between two CL-lines or G01 commands. The maximum TCP contouring or corner blending error ε T CP occurs in the centre of the cornering trajectory and is calculated by evaluating the interpolated axis displacements s at t = T 1 . The interpolated axis displacements are calculated from (27) and the vectors from the corner transition to these positions are represented by l 1 and l 2 .
The contouring error ε T CP (shown in Fig.12) is calculated from the Euclidean distance between the vectors l 1 and l 2 . where θ T CP represents the TCP cornering angle. Assuming constant feedrate in this example, l 1 = l 2 = l ε , in which case, (28) simplifies to the following expression: Inserting (16) and (25) into (28) enables the TCP corner blending error to be defined as: Using equation 30 the TCP error can be calculated for any toolpath geometry and commanded feedrate. The kinematic profiles for the 3-FIR case are shown in Fig.21 in appendix C and the derivation of TCP error for the 3-FIR case is included in appendix D.
To ensure minimum cycle times the actual feedrate must remain as close to the commanded feedrate as possible throughout the toolpath including during cornering transitions. To satisfy both jerk and TCP error constraints, however, there is maximum permissible cornering feedrate. Using equations (30) and (44), it is possible to calculate the relationship between TCP error, maximum permissible cornering feedrate and cornering angle for the 2-FIR and 3-FIR filter cases respectively.
Rearranging equation (30), the maximum permissible cornering feedrate for the 2-FIR filter case must satisfy: and for the 3-FIR filter using (44) the maximum permissible cornering feedrate must satisfy: 32) For a commanded feedrate F and range of cornering angles θ T CP ∈ [0 • , 180 • ], equations (31) and (32) are solved for solutions 0 ≤ α ≤ 1 to calculate the limit to the feedrate scaling factor α. When multiplied by the commanded feedrate F this represents the maximum permissible cornering feedrate that can be achieved whilst satisfying the kinematic and tolerance constraints. The blending pulse feedrate F c is commanded to this limit value.
The feedrate limit for the 2-FIR filter case is shown in Fig.13. Cornering feedrates selected below the curves will satisfy the TCP error constraints for the commanded feedrate and cornering angle. The figure shows the limits for 10µm and 50µm tolerance constraints. For the 50µm tolerance, the figure shows higher cornering feedrates can be achieved compared to the 10 µm case. The cornering feedrate limits for both the 2-FIR and 3-FIR filtered cases are compared in Fig.14. It can be seen that higher feedrates can be achieved in the 3-FIR case which satisfy the tolerance and jerk constraints. Therefore there is an advantage of using 3-FIR filters to reduce the overall machining cycle time as the tool can remain at higher feedrates during cornering transitions than for the 2-FIR case. Despite the advantage of using a higher order filter, there remains a limit to the order of filters that can be used effectively for trajectory generation. As the order is increased the filter time constant reduces. In the frequency domain the notch (as shown in Fig.5) will shift to higher frequencies. This will be constrained by the lowest structural mode of the machine tool. This section has shown a method of using multiple first order FIR filters with matching time constants to model continuous linear interpolation of velocity pulse signals. It has been show that the cornering feedrate and TCP error can be controlled using velocity blending pulses. This method has been extended to predict feedrates and machining cycle time for toolpaths of any geometry and defined tolerance. The following section demonstrates and validates the proposed method on industrial case studies.

Experimental Validation
Machining experiments were conducted on a DMG Mori Universal eVo 40 5axis machining centre with a Heidenhain TNC640 controller. Two short toolpaths were used for pocketing operations and a single long aerospace part program is evaluated in the cycle time prediction.

Case Studies 1 & 2 -Pocketing Toolpaths
The first two case studies, as shown in Fig.15, consist of a contour and a trochoidal pocketing tool-path. These tool-paths are generated by CAM software [28] and the part programs are deployed to the machine directly with no modification. Table 1 shows the cutting conditions. As noted, 2 different feedrates 1000 and 3000 mm/min are used. The most important setting is contour error tolerance for HSM. Two different contouring tolerance, 10 and 50 µm are used. Table 1 summarises the cycle time results. All simulated trajectories in the case studies were modelled using the method described in section 4.2 and 3-FIR filters.

Machining Cycle Time Estimation
The predicted machining cycle times are compared with the measured CNC and CAD/CAM calculated machining cycle times. The results are presented in Table 1 and Fig.16a. For all cases the predicted machining cycle times are accurate to within 3% of the measured cycle time with the exception of the trochoidal pocket (1000 mm/min, 10µm case) which is 5.52%. These compare favourably to the CAD/CAM calculated cycle times which has an error range from 0.22% to 54.99%. The significant result is the Trochoidal pocket (3000 mm/min and 10µm case). The proposed method is able to accurately predict the increase in machining cycle time from 14.40 to 30.99 seconds when tightening the tolerance from 50µm to 10µm, which is within 2.72% of the measured cycle time. This is compared to an error of 54.99% for the CAD/CAM calculated cycle time.

Feedrate Prediction
To demonstrate the performance of the feedrate prediction method a number of toolpath features were selected. The predicted, CAD/CAM calculated and measured CNC tangential velocities at these particular features were recorded and are presented in Tables 2 and 3. The contour pocket features consist of (1) a long G01 segment, (2) a sharp corner and (3) a rounded corner consisting of small G01 segments. The trochoidal pocket features consist of (1) the stepover segment and (2) the main arc. Depending on the tolerance and the commanded feedrate large differences in tangential velocity can exist between the stepover segment and the main arc of a trochoidal toolpath which in turn results in a large cyclical variation of cutting forces. It is for this reason they are included in this study. The features described above are shown on the toolpaths in Fig. 16 and the corresponding position with respect to displacement and tangential velocities are demonstrated directly beneath. Overall, the prediction error ranges from 0.1-10.3% compared with CAD/CAM calculated error range of 0.22-2555%, where the error is calculated as a percentage difference from, and with respect to, the measured tangential velocity. The performance of the proposed feedrate prediction method at each feature is described below: Long G01 Segment. The prediction error range is between 0.1-0.13% compared to the CAD/CAM calculated error range of 0.27-0.3%. The high accuracy is to be expected as no feedrate limiting features are present in the segment. The difference in measured velocity compared to the idealised CAD/CAM values are due to interpolator rounding during trajectory generation.
Sharp Corner. The prediction error range is between 7-20% compared to the CAD/CAM calculated error range of 107-2555%. The fundamental difference is due to the CAD/CAM calculation not taking into account the cornering kinematic constraints due to tolerance and thus not predicting the reduction in feedrate during the cornering segment. This holds true for all of the features demonstrated except the long G01 segment. For the 10µm tolerance cases the tool comes to an almost complete stop -4% and 11% of the commanded feedrate for the 3000 mm/min and 1000 mm/min cases respectively, the presented method predicts these reductions.
Rounded Corner. The prediction error range is between 0.3-5.4% compared to the CAD/CAM calculated error range of 0.2-168%. The significant result is the 3000 mm/min and 10µm case (Fig. 16c) where the CAD/CAM calculation does not account for the reduction in velocity due to the tolerance requirement. The CAD/CAM calculated error is 168% compared to the measured value and the prediction error is within 2.4%.
Trochoid Stepover. The prediction error range is between 0.1-10% compared to the CAD/CAM prediction error range of 0.2-0.5%. The CAD/CAM calculation does not predict any differences along the trochoidal toolpath between the stepover and the main arc. This can be seen in Fig. 16d for the 3000 mm/min 10µm case. The blue line shows the CAD/CAM prediction but the actual kinematic profile is very different. The stepover results in tangential velocities close to the commanded feedrate as the cornering angles between the segments are less acute than for the rest of the main arc.
Trochoid Main Arc. The prediction error range is between 0.6-7% compared to the CAD/CAM calculated error range of 0.2-208%. The reduction in tangential velocity around the main arc is due to the cornering angles between the segments. The influence of the toolpath tolerance on the cornering tangential velocity can be seen in Fig.16d and Fig.16f. The increase in tolerance from 50µm to 10µm results in more than a 65% reduction in tangential velocity around the main arcs of the trochoids. The prediction method accurately predicts the feedrate within 1.5% of tangential velocity measured at the main arc. Taking this result one step further, this demonstrates that a feedrate driven cutting force model when incorporating the prediction method will be able to predict the cyclical cutting forces due to the 65% variation in magnitude of feedrate fluctuations around the trochoidal toolpath.

Case Study 3 -Aerostructure Toolpath
An industrial toolpath was chosen to validate the method against a representative aerostructure part. The part program consists of three toolpathsroughing, finishing #1 floors and finishing #2 walls as shown in Fig. 17. The part programs were run at three tolerance settings, 10µm, 20µm and 50µm  Table 3 Trochoidal pocket case study: tangential velocity prediction and performance to demonstrate the significant impact tolerance has on machining cycle times and therefore on feedrate and cycle time prediction. Table 4 compares the predicted machining cycle times with both the measured cycle times and the predicted times from a commercial CAD/CAM software package for each individual toolpath. The overall machining cycle times, calculated by summing the cycle times for the 3 sections of the part program, are shown in Table 5. The CAD/CAM prediction error ranges from 62.41% under prediction for the 10µm case to 36.42% under prediction for the 50µm case. The actual CAD/CAM predicted times do not change as the software does not account for tolerance, the calculation is based upon distance travelled along the toolpath and ideal feedrate. Therefore as the tolerance is relaxed the measured cycle time approaches the CAD/CAM case and their prediction becomes more accurate.
The prediction error from the proposed method (as shown in Table 5) ranges from 3.50% over prediction for the 10µm case to 4.69% for the 50µm case. The 20µm case has a prediction error of 5.34% under the measured cycle time which is approximately 10% of the CAD/CAM error (51.49%) for that particular case. The aerostructure case study validates the model for predicting    Table 5 Total machining cycle times and errors for measured, predicted and CAD/CAM.

Cutting Force Prediction
Lastly, the importance of accurate feedrate prediction for virtual machining models is demonstrated. This is realized by estimating cutting forces along the complex trochoidal toolpath shown in Fig.19. Predicting the cutting forces, considering the complex tool engagements on this toolpath, is realized by adapting the cutting force prediction model presented in [4] with the proposed feedrate prediction method. Readers should refer to [2] and [4] for details of the cutting force model.

Case Study 4 -Accurate Cutting Force Prediction using Predicted Feedrates
To validate the feedrate prediction method with a cutting force model machining trials were conducted on the 5-axis DMG Mori eVo 40 machining centre fitted with a Heidenhain TNC640 controller. The toolpath, shown in Fig. 19, was designed using NX CAM as a trochoidal pocketing operation. A 40mm x 60mm x 10mm open sided pocket was selected as the test feature as shown in Fig. 18. A 2-fluted 12mm solid carbide end mill with a HSK-63A tool holder was used. The workpieces were 236mm x 30mm x 6mm aluminium 7075, each held using a Geradi compact grip vice mounted to the dynamometer. A Kistler 9139AA dynamometer and a National Instruments USB-6343 multi-channel DAQ was used to acquire cutting force data at 10kHz. The machining centre was connected to a local area network via a RJ45 network cable such that the machine controller data was accessed by two methods. The first using a pre-defined MTConnect datastream through a TCPIP connection at 20Hz and the second using an LSV2 protocol direct to the controller through a TCPIP connection at 111Hz.
The predicted cutting forces during the trochoidal section are shown in Fig.  19. The peak predicted cutting force for the standard feedrate model is 673N compared to 380N for the filtered feedrate model, from the peak measured cutting forces this gives prediction errors of 96.2% and 10.8% respectively. In the cornering section of the toolpath the peak predicted cutting force for the standard feedrate model is 821N compared to 656N for the filtered feedrate model, from the peak measured cutting forces this gives prediction errors of 37.3% and 9.7% respectively. The validation trials show that the inclusion of an accurate feedrate profile in the cutting force model enables a more accurate prediction of cutting forces for complex toolpaths.

Conclusions
A novel method of accurately modelling the trajectory generation of NC systems has been proposed. The main conclusions from this research are as follows: 1. An accurate method of feedrate prediction along short-segmented complex tool-paths was introduced. 2. The linear interpolation dynamics and commanded axis kinematic profiles of NC systems were predicted using both 2 and 3 first order Finite Impulse Response filters with the same time constant. 3. The corner blending behaviour during non-stop interpolation of linear segments was modeled by introducing velocity blending pulses. 4. For the first time, the minimum cornering feedrate, that satisfies both the tolerance and machining constraints, has been calculated analytically for toolpaths of any geometry. 5. The reduction in machining cycle time by using 3 FIR filters compared to 2 FIR filters was proven analytically. 6. The feedrate prediction method was validated experimentally against four different case studies demonstrating industrial 3-axis machining tool-paths. 7. The proposed method demonstrated cycle times can be estimated with >90% accuracy, greatly outperforming CAM-based predictions. 8. The predicted feedrate method was incorporated into a cutting force model, demonstrating an increase in cutting force accuracy for a complex toolpath, and validated experimentally.
Further work will integrate the methods into virtual machining and digitaltwin models and extend the method to 5-axis machining.

C Kinematic Profiles for 3-FIR Case
where b 1 = 3 2 T 1 − T b . The maximum TCP error occurs at t = T d 2 = 3 2 T 1 for the 3 first order FIR filter case, the interpolated axis velocity and displacement, are defined as equations (40) and (41) respectively: Using equation (16), (40)