PSO-Based Feedrate Optimization Algorithm for Five-Axis Machining with Constraint of Contour Error

Abstract

Feedrate has a great influence on contour error in five-axis machining. Accordingly, it is of great significance to plan the time-optimal feedrate curve considering the contour error constraint to achieve high-accuracy and high-efficiency machining. Aiming at improving the error control accuracy of model linearization loss and optimizing the machining time, the PSO-based feedrate optimization algorithm for five-axis machining with constraint of contour error is proposed in this paper. Firstly, the relationship between parametric feedrate and contour error constraint is clarified that provides a model basis for accurately controlling contour error by optimizing the feedrate curve. Then, the feedrate optimization model, which takes the control vertices of the feedrate curve expressed by B-spline as the decision variables and minimizes the machining time as the optimization objective, is established. Subsequently, to overcome the shortcomings of low accuracy and low efficiency caused by single optimization of global control vertices, the group search particle swarm optimization (GSPSO) algorithm based on window movement is adopted to optimize the feedrate curve in segments. Finally, the effectiveness of the proposed feedrate optimization algorithm is validated by three typical test toolpaths on an open double-turntable five-axis machine tool. In light of the experiment, the proposed algorithm is able to fully release the potential of the machine tools while accurately controlling the contour error of the cutter tip and cutter orientation.

1. Introduction

With the rapid increase in the demand for complex surface parts in aerospace and energy fields, five-axis machine tools, with the features of good obstacle avoidance and high motion flexibility, have been widely used. To achieve high-precision and high-efficiency machining on the five-axis machine tools, it is necessary to explore the contour error control method with optimal machining time. At present, the control of contour error is mainly realized by redesigning the servo controller structure [1,2,3,4], pre-compensation [5,6,7,8] and feedrate optimization. However, users usually are not capable of implementing these advanced control strategies to improve contour accuracy due to the confidentiality of the built-in controller structure of most commercial machine tools. In addition, the pre-compensation technology also has the problem that the contour control effect in the high-curvature area of the path deteriorates due to the decrease in path smoothness. By comparison, contour accuracy is improved by optimizing the feedrate curve, which can not only overcome the limitations of the above two methods, but also has higher flexibility and stronger robustness. Therefore, how to plan the time-optimal feedrate curve under the premise of ensuring contour accuracy is worthy of in-depth study.

In the initial research, in order to reduce the influence of curve interpolation error, the adaptive feedrate planning algorithm considering chord error was proposed [9,10]. Subsequently, some studies [11] have shown that limiting the kinematic characteristics of the cutter tip velocity, acceleration and jump in the process of machining can improve the forming quality of the parts and the service life of the cutters. Thus, the linear [12], exponential [13], s-curve [14] and sine-curve [15] acceleration and deceleration strategies that limit the kinematics characteristics of the cutter tip have been proposed successively. Notice that the abovementioned feedrate planning methods mainly focus on controlling the chord error and the kinematic features of the feed cutter, while ignoring the driving property of the machine tool feed axis. Especially in five-axis machining, due to the complex nonlinear mapping relationship between the motion of the feed cutter and the motion of the feed axis, even if the feedrate of the cutter is set to a constant value, the motion state of the feed axis may exceed its own driving capacity, resulting in the reduction in machining accuracy. Consequently, it is necessary to consider the drive constraints in the process of planning the feedrate curve. For instance, Sun [16] proposed a curve evolution strategy to iteratively adjust the initial feedrate curve to meet the jerk constraint requirements. Sencer [17] proposed a feedrate optimization method considering the Bang-Bang control of drive constraints, which realized the control of the motion state of the feed axis. Since then, the linear programming algorithm [18], bidirectional scanning algorithm [19], heuristic search algorithm [20] and greedy idea algorithm [21] have been proposed to control the driving performance of the feed axis.

It should be pointed out that although the above feedrate optimization algorithm is able to improve the contour accuracy, it does not directly consider the contour error constraint. To achieve the goal of accurately controlling the contour error by optimizing the feedrate curve, Lin [22] first established a nonlinear function model with the feedrate and the contour error respectively regarded as the independent variable and the dependent variable, and obtained the approximate value of the maximum feasible feedrate satisfying the preset contour error constraint. Jia [23] approximated the whole position closed-loop control system to a typical second-order under-damped system and deduced the explicit functional relationship of contour error in regard to the cutter feedrate, the curvature of the toolpath and the damping ratio of the system. On the basis of revealing the potential linear relationship between feed axis tracking error, feed axis velocity and feed axis acceleration, Wang [24] clarified the relationship between feedrate and the contour error constraint. According to the established feed servo system model, Chen [25] deduced the analytical equation between feedrate and contour error constraint. In addition, Erwinski [26] proposed a method to control contour error by optimizing feedrate curve with PSO combined with enhanced Lagrange constraint processing technology, based on the prediction of contour error using an artificial neural network. It is noteworthy that the above methods are difficult to effectively apply to five-axis machining in consideration of the complex nonlinear coupling property of five-axis kinematics transformation. For five-axis machining, Chen [27] clarified the linearization relationship between the contour of the cutter tip and cutter orientation with respect to feedrate by simplifying the tracking error prediction model and contour error estimation model.

To sum up, the current research on controlling the contour error of the cutter tip and cutter orientation by optimizing the feedrate curve is still very limited, and the linearization of the model will inevitably lead to the decline in error control accuracy. Aiming at handling the above problem, this paper proposes the PSO-based feedrate optimization algorithm for five-axis machining with constraint of contour error considering the advantages of PSO in solving nonlinear optimization problems in the fields of environmental sensing [28], space science [29], the design of time-delay equalizers [30] and spatial phase shifters [31]. Firstly, a nonlinear and high-precision feedrate optimization model with constraint of contour error is established on the basis of the dynamic prediction of tracking error and the iterative estimation of contour error. After that, for obtaining the target feedrate curve accurately and efficiently, the GSPSO algorithm based on window movement is applied to plan the feedrate curve in segments. The rest of the organization structure of the paper is as follows: The NURBS toolpath and interpolation algorithm are introduced in Section 2. Section 3 illustrates the details of the proposed feedrate optimization algorithm. The experimental results and analysis are presented in Section 4. The paper is summarized in Section 4

2. NURBS Toolpath and Interpolation Algorithm

In the workpiece coordinate system (WCS), the five-axis parametric toolpath is usually made up of the cutter tip position vector $\mathit{p}\left(u\right)$ and the unit cutter orientation vector $\mathit{o}\left(u\right)$. NURBS is often used to express the five-axis parametric toolpath due to its accurate analytical expression and flexible local readjustment capability. The NURBS toolpath usually adopts normalized arc length parameterization. Therefore, the expression of the cutter tip position vector and the unit cutter orientation of the five-axis trajectory defined on the normalized arc length parameter $u ϵ \left[0,1\right]$ is

$$\begin{array}{c}{\left[\begin{array}{c}\mathit{p}\left(u\right)\\ \mathit{o}\left(u\right)\end{array}\right]}^T=\left[p_x\left(u\right),p_y\left(u\right),p_z\left(u\right),o_x\left(u\right),o\left(u\right),o_x\left(u\right)\right]\end{array}$$

(1)

The $G^2$ continuous dual NURBS are adopted to express the following five-axis parametric trajectory:

$$\begin{array}{c}\left\{\begin{array}{c}\mathit{p}\left(u\right)={\displaystyle \sum _{i=0}^n}{N}_{i,3}\left(u\right){\omega }_i{\mathit{p}}_i/N_{i,3}\left(u\right){\omega }_i\\ \mathit{q}\left(u\right)={\displaystyle \sum _{i=0}^n}{N}_{i,3}\left(u\right){\omega }_i{\mathit{q}}_i/N_{i,3}\left(u\right){\omega }_i\end{array}\right.\end{array}$$

(2)

where $\mathit{p}\left(u\right)$ and $\mathit{q}\left(u\right)$, respectively, represent the cutter tip trace and the trace of a point on the cutter axis that is different from the cutter tip. ${\mathit{p}}_i$ and ${\mathit{q}}_i$ are the control vertices, $n+1$ is the number of control vertices and ${\omega }_i$ is the weight factor of the NURBS. $N_{i,3}\left(u\right)$ is the function defined on the following uniform node vector $\mathit{U}$:

$$\begin{array}{c}\mathit{U}=\left[\stackrel{4}{\overbrace{0,\cdots ,0}},\frac{1}{n-2},\cdots ,\frac{n-1}{n-2},\stackrel{4}{\overbrace{1,\cdots ,1}}\right]\end{array}$$

(3)

After that, the unit cutter orientation vector $\mathit{o}\left(u\right)$ is calculated as

$$\begin{array}{c}\left\{\begin{array}{l}\mathit{o}\left(u\right)=\frac{\mathit{p}\left(u\right)-\mathit{q}\left(u\right)}{‖\mathit{p}\left(u\right)-\mathit{q}\left(u\right)‖}\\ ‖\mathit{o}\left(u\right)‖=1\end{array}\right.\end{array}$$

(4)

Taking the A–C double-turntable five-axis machine tool shown in Figure 1 as an example, the toolpath ${\left[\mathit{p}\left(u\right),\mathit{o}\left(u\right)\right]}^T$ in the WCS can be converted to the following motion trajectory $\mathit{m}\left(u\right)$ of the feed axis in the machine coordinate system (MCS) by using the inverse kinematics transformation formula shown in Equation (6):

$$\begin{array}{c}\mathit{m}\left(u\right)={\left[X\left(u\right),Y\left(u\right),Z\left(u\right),A\left(u\right),C\left(u\right)\right]}^T\end{array}$$

(5)

where $\left[X,Y,Z\right]$ and $\left[A,C\right]$, respectively, represent linear axes and rotary axes in the MCS. The inverse kinematics transformation formula shown in Figure 1 is as follows:

$$\begin{array}{c}\left\{\begin{array}{l}A=\arccos \left(O_z\right)\\ C=\arctan (O_x/O_y)\\ X=\cos \left(C\right)p_x-\sin \left(C\right)p_y\\ Y=\cos \left(A\right)\sin \left(C\right)p_x-\cos \left(A\right)\cos \left(C\right)p_y-\sin \left(A\right)p_z-\sin \left(A\right)L_{ac,z}\\ Z=\sin \left(A\right)\sin \left(C\right)p_x-\sin \left(A\right)\cos \left(C\right)p_y+\cos \left(\mathrm{A}\right)p_z+\cos \left(A\right)L_{ac,z}+L_{Tya,z}\end{array}\right.\end{array}$$

(6)

where $L_{Tya,z}$ and $L_{ac,z}$ are turntable offsets. Correspondingly, the forward kinematics transformation formula shown in Figure 1 is as follows:

$$\begin{array}{c}\left\{\begin{array}{l}{o}_x=\sin \left(A\right)\sin \left(C\right)\\ o_y=\sin \left(A\right)\cos \left(C\right)\\ o_z=\cos \left(A\right)\\ p_x=-\cos \left(C\right)X+\cos \left(A\right)\sin \left(C\right)Y+\sin \left(A\right)\sin \left(C\right)Z-\sin \left(C\right)\sin \left(A\right)L_{Tya,z}\\ p_y=-\sin \left(C\right)X-\cos \left(A\right)\cos \left(C\right)Y-\sin \left(A\right)\cos \left(C\right)Z+\cos \left(C\right)\sin \left(A\right)L_{Tya,z}\\ pz=-\sin \left(A\right)Y+\cos \left(A\right)Z-\cos \left(A\right)L_{Tya,z}-L_{ac,z}\end{array}\right.\end{array}$$

(7)

In order to interpolate the NURBS toolpath, it is necessary to clarify the relationship between arc length increment and parameter increment; however, the relationship is nonlinear due to its dependence on the shape of the toolpath. Therefore, the parameter increment corresponding to each step of interpolation is usually determined by Taylor approximation based on the current known feedrate. This paper adopts the following second-order Taylor formula:

$$\begin{array}{c}{u}_{k+1}=u_k+\frac{f\left(u_k\right)}{\left|\left|P_u\left(u_k\right)\right|\right|}{T}_s+\left(\frac{a\left(u_k\right)}{\left|\left|P_u\left(u_k\right)\right|\right|}-\frac{P_u\left(u_k\right)·P_{uu}\left(u_k\right)}{{\left|\left|P_u\left(u_k\right)\right|\right|}^4}{f\left(u_q\right)}^2\right)\frac{{T_s}^2}{2}\end{array}$$

(8)

where $f\left(u_k\right)$ and $a\left(u_k\right)$ are the feedrate and acceleration of the cutter at the interpolation sampling parameter $u_k$, respectively, and $T_s$ is the interpolation sampling period.

3. The Proposed Feedrate Optimization Algorithm

In this section, the prediction method of contour error of the cutter tip and cutter orientation is first introduced. Then, the feedrate optimization model aiming at the minimum machining time is established and the PSO algorithm for control vertices optimization is proposed. Finally, the implementation details of the feedrate optimization are clarified.

3.1. Contour Error Prediction

In five-axis machining, the tracking error caused by the response delay of the feed servo system will be delivered to the cutter through five-axis kinematics transformation, resulting in the contour error of the actual machining trajectory. In addition, the contour error is often obtained by estimation due to the complexity of solving the foot point parameter. Thus, the establishment of high-precision tracking error prediction and contour error estimation is the basis of accurate contour error prediction.

3.1.1. Tracking Error Prediction

For high-performance complex surface parts in CNC machining, to guarantee the surface processing quality of parts, the feed servo system usually does not allow overshoot. Therefore, technicians commonly adopt parameter coordination methods to make the entire CNC machining system run under over-damped conditions. At this time, the feed servo system is capable of being reasonably simplified to the typical first-order inertia system [32] shown in Equation (9).

$$\begin{array}{c}G\left(s\right)=\frac{{\tau }_a\left(s\right)}{{\tau }_r\left(s\right)}=\frac{1}{Ts+1}\end{array}$$

(9)

where ${\tau }_r\left(s\right)$ and ${\tau }_a\left(s\right)$, respectively, represent the desired input position and actual output position of the servo system in the complex frequency domain, and $T$ is the system time constant. It should be pointed out that the time constant $T$ of the servo system is usually several times that of the interpolation sampling cycle $T_s$ of the CNC system. Moreover, the adjustment time for the servo system to reach the approximate steady-state condition is about $(3~4)T$, due to the limitation of the response bandwidth of the servo system, which also means that it is actually difficult for the servo system to reach a steady state in one interpolation sampling period. Therefore, it is very difficult to achieve high-precision tracking error prediction based on the assumption of steady-state response of the system. For the accurate prediction of contour error, we will use a dynamic prediction method based on an interpolation sampling position sequence to establish a tracking error prediction model with higher accuracy.

In actual CNC machining, to satisfy the requirements of various geometry, drive and process constraints, along the predefined toolpath, the feedrate usually takes the form of a step change with the interpolation sampling period $T_s$ as the time interval. Accordingly, it is also easy to find that the velocity of the servo system is also a fixed value within the same sampling period. In this case, ignoring the influence of filtering denoising inside the system, the position signals received by the servo system are essentially a set of variable ramp input signals. In order to analyze the tracking error under the variable ramp input signals, it is assumed that the current sampling moment is the end of the $k$th sampling period, and in addition to the ramp signal of the $(k+1)$th sampling period sent by the hose computer, the residual servo lag error of the $k$th sampling period will also be a part of the system input signal to drive feed axis movement [6]. At this moment, the linear superposition of the step signal and the ramp signal can be approximated as the input signal of the system. For the feed servo system with a single input and single output, the output of the system is also the linear superposition of the output response of the above the two typical signals based on the principle of linear superposition. Therefore, the feed axis tracking error $e_{\tau }\left(k{T}_s\right)$ at the $k$th sampling moment is calculated as

$$\begin{array}{c}{e}_{\tau }\left(k{T}_s\right)=e_{\tau }^{ramp}\left(k{T}_s\right)+e_{\tau }^{step}\left(\left(k-1\right)T_s\right){\widehat{e}}^{-\frac{T_s}{T}}\end{array}$$

(10)

with

$$\begin{array}{c}\left\{\begin{array}{l}{e}_{\tau }^{ramp}\left(k{T}_s\right)=\frac{{\tau }_r\left(k\right)-{\tau }_r\left(k-1\right)}{T_s}T(1-{\widehat{e}}^{-\frac{T_s}{T}})\\ e_{\tau }^{step}((k-1)T_s)=e_{\tau }\left(\left(k-1\right)T_s\right)\end{array} 1\le k\le N\right.\end{array}$$

(11)

where $\tau =(X,Y,Z,A,C)$ represents the machine tool feed axis, $\widehat{e}$ and $N$, respectively, indicate the natural constant and the number of reference interpolation sampling position sequences. ${\tau }_r\left(k\right)$ donates the discrete position of the feed axis, and $e_{\tau }^{step}((k-1)T_s)$ and $e_{\tau }^{ramp}\left(k{T}_s\right)$, respectively, represent the servo lag error of the system at the $(k-1)$th sampling moment and the tracking error introduced by the ramp signal ${(\tau }_r\left(k\right)-{\tau }_r\left(k-1\right))/T_s$. Correspondingly, the actual position of the feed axis at the $k$th sampling moment is calculated as

$$\begin{array}{c}{\tau }_a\left(k\right)={\tau }_r\left(k\right)-e\left(k{T}_s\right)\end{array}$$

(12)

Equations (10) and (12) show that the output position of the feed axis at the $k$th sampling moment is not only related to the input signal at the current moment but is also affected by the previous input signal. For a given feed servo system, $\widehat{e}$, $T$ and $T_s$ in Equation (10) are fixed constants. Therefore, Equations (10) and (11) essentially reveal the potential linear relationship between the system reference position sequence and the feed axis tracking error, which also provides a model basis for the subsequent precise prediction of contour error.

3.1.2. Contour Error Estimation

On the five-axis toolpath as shown in Figure 2, $({\mathit{p}}_a,{\mathit{o}}_a)$ is the actual cutter position, $({\mathit{p}}_r,{\mathit{o}}_r)$ is the reference cutter position, and the cutter position deviation between $({\mathit{p}}_a,{\mathit{o}}_a)$ and $({\mathit{p}}_r,{\mathit{o}}_r)$ is the tracking error $({\mathit{e}}_p,{\mathit{e}}_o)$, where ${\mathit{e}}_p$ is the tracking error of the cutter tip, and ${\mathit{e}}_o$ is the tracking error of the cutter orientation, as shown by the green solid line in Figure 2. The cutter tip position ${\mathit{p}}_n$ is the point closest to the actual cutter tip position ${\mathit{p}}_a$ on the toolpath, also known as the station point. The deviation of the cutter position between $({\mathit{p}}_a,{\mathit{o}}_a)$ and $({\mathit{p}}_n,{\mathit{o}}_n)$ is the contour error $({\mathit{\epsilon }}_p,{\mathit{\epsilon }}_o)$, where ${\mathit{\epsilon }}_p$ and ${\mathit{\epsilon }}_o$, respectively, represent the contour error of the cutter tip and cutter orientation, as shown by the red solid line in Figure 2.

Accordingly, calculating the parameter $u_n$ corresponding to the station point $p_n$ is a prerequisite for obtaining the contour error of the cutter tip and cutter orientation. However, due to the fact that $\mathit{p}\left(u\right)$ is a free curve equation, the parameter $u_n$ corresponding to the equation is the nonlinear equation, making it difficult to obtain an accurate analytical solution. To handle this problem, the initial value regeneration Newton method [33] is used to efficiently and accurately obtain an approximate numerical solution for the parameter $u_n$. The numerical solution $u_f$ is calculated as

$$\begin{array}{c}{u}_f=u_b-\frac{{\parallel {\mathit{p}}^{\prime }\left(u_b\right)\parallel }^2\left(\mathit{p}\left(u_b\right)-{\mathit{p}}_a\right){·{\mathit{p}}^{\prime }\left(u_b\right)}^T}{\begin{array}{c}{\parallel {\mathit{p}}^{\prime }\left(u_b\right)\parallel }^4+{\parallel {\mathit{p}}^{\prime }\left(u_b\right)\parallel }^2\left(\mathit{p}\left(u_b\right)-{\mathit{p}}_a\right)·{{\mathit{p}}^{″}\left(u_b\right)}^T\\ -\left(\left(\mathit{p}\left(u_b\right)-{\mathit{p}}_a\right){{·\mathit{p}}^{\prime }\left(u_b\right)}^T\right)\left({\mathit{p}}^{\prime }\left(u_b\right){·{\mathit{p}}^{\prime \prime }\left(u_b\right)}^T\right)\end{array}}\end{array}$$

(13)

where $u_b$ is the initial value of regeneration and it is calculated as

$$\begin{array}{c}{u}_b=u_r-\left(\mathit{p}\left(u_r\right)-{\mathit{p}}_a\right)·\frac{{{\mathit{p}}^{\prime }\left(u_r\right)}^T}{\parallel {\mathit{p}}^{\prime }\left(u_r\right)\parallel }\end{array}$$

(14)

where $u_r$ represents the corresponding parameter for the cutter tip position ${\mathit{p}}_r$.

According to the convergence of Newton’s algorithm, the parameter $u_n$ corresponding to the station point $p_n$ is calculated as

$$\begin{array}{c}{u}_n=\left\{\begin{array}{ll}{u}_f,& \left|d_t\left(u_f\right)\right|<{|d}_t\left(u_b\right)|\\ u_b-\left(\mathit{p}\left(u_b\right)-{\mathit{p}}_a\right)·\frac{{{\mathit{p}}^{\prime }\left(u_b\right)}^T}{\parallel {\mathit{p}}^{\prime }\left(u_b\right)\parallel },& \left|d_t\left(u_f\right)\right|\ge \left|d_t\left(u_b\right)\right|\end{array}\right.\end{array}$$

(15)

On the basis of the abovementioned concept of contour error, the cutter tip error ${\mathit{\epsilon }}_p$ and the cutter orientation error ${\mathit{\epsilon }}_o$ can be calculated as

$$\begin{array}{c}{\mathit{\epsilon }}_p=\parallel p\left(u_n\right)-{\mathit{p}}_a\parallel \end{array}$$

(16)

$$\begin{array}{c}{\mathit{\epsilon }}_o=\parallel o\left(u_n\right)-{\mathit{o}}_a\parallel \end{array}$$

(17)

3.2. Feedrate Optimization Model

According to the known five-axis parameter trajectory ${\left[\mathit{p}\left(u\right),\mathit{o}\left(u\right)\right]}^T$, the feedrate curve $f\left(u\right)$ can be expressed by the $G^2$ continuous B-spline defined in the same parameter $u$

$$\begin{array}{c}f\left(u\right)={\displaystyle \sum _{i=0}^M}{d}_i{N}_{i,3}\left(u\right)\end{array}$$

(18)

where $d_i$ and $M+1$, respectively, represent the control vertices and the number of control vertices of the B-spline, and $N_{i,3}\left(u\right)$ here is set as Equation (2).

Along the given toolpath, the goal of feedrate optimization is to acquire the shortest machining time $T_{obj}$ by adjusting the feedrate curve $f\left(u\right)$ shape while meeting the preset contour error constraint. Therefore, the feedrate optimization model can be established as follows:

$$\begin{array}{c}min{T}_{obj}=min{\int }_0^1\frac{dt}{du}du\end{array}$$

(19)

$$\begin{array}{c}s.t.\left\{\begin{array}{c}{\mathit{\epsilon }}_p\left(u\right)\le {\mathit{\epsilon }}_{pmax}\\ {\mathit{\epsilon }}_o\left(u\right)\le {\mathit{\epsilon }}_{omax}\end{array},u ϵ\left[\mathrm{0,1}\right]\right.\end{array}$$

(20)

where ${\mathit{\epsilon }}_{pmax}$ and ${\mathit{\epsilon }}_{omax}$, respectively, indicate the upper bound of the contour error constraints of the cutter tip and cutter orientation. According to the analysis of the feedrate optimization model, the model has the characteristic of a strong coupling nonlinearity. To accurately and efficiently solve the optimization model and gain the optimal target feedrate curve, this paper proposes a general solution with PSO combined with a novel constraint processing technology. Furthermore, it is worth noting that the optimized machining time is paid more attention than the calculation time for the off-line feedrate optimization algorithm.

3.3. Control Vertices Optimization Based on PSO Algorithm

In light of the definition of feedrate curve $f\left(u\right)$, the feedrate curve shape is completely controlled by the control vertices on the premise that the knot vector $\mathit{U}$ remains unchanged. As a consequence, the control vertices ${\left\{d_i\right\}}_{i=0}^M$ can be selected as the decision variables of the feedrate optimization model. The PSO algorithm is used for solving the feedrate optimization model in this paper and it determines the optimal solution step by step through the exchange of information between an individual and a population. Different from other optimization algorithms, the PSO algorithm not only has few parameters and is easy to implement, but also has a fast solving speed and a strong global parallel search ability. Additionally, it has small restrictions on population size and is capable of recording population information [34]. The mathematical model considering contour error constraints and taking control vertices as decision variables is first established in this section. After that, the process of using the GSPSO algorithm based on window movement to solve the decision variables is introduced at great length.

3.3.1. Mathematical Model

Consider the convex hull property of the spline curve, the search boundary condition of the control vertices is set as $0\le {\left\{d_i\right\}}_{i=0}^M\le f_{max}$, where $f_{max}$ is the maximum programmed feedrate, whose value is usually determined by surface quality, cutter life, the removal rate of material, cutting force and other factors. Obtaining the time-optimal feedrate curve means seeking the maximum feasible solution of the feedrate under the contour error constraints. Consequently, the objective function shown in Equation (19) $T_{obj}$ is converted into the new objective function $F_{obj}$ as follows:

$$\begin{array}{c}{minF}_{obj}=-{\displaystyle \sum _{k=0}^N}{\displaystyle \sum _{i=0}^M}{d}_i{N}_{i,3}\left(u_k\right)(u_k-u_{k-1})\end{array}$$

(21)

$$\begin{array}{c}s.t.\left\{\begin{array}{c}{g}_1\left(u_k\right)=\frac{{\mathit{\epsilon }}_p}{{\mathit{\epsilon }}_{pmax}}-1\le 0\\ g_1\left(u_k\right)=\frac{{\mathit{\epsilon }}_o}{{\mathit{\epsilon }}_{omax}}-1\le 0\end{array}\right.,0\le k\le N\end{array}$$

(22)

where $g_1\left(u_k\right)$ and $g_1\left(u_k\right)$, respectively, denote the inequality constraints at the interpolation sampling parameter $u_k$. The normalization of the constraints and objective function is performed to avoid the optimization process being affected by their order of magnitude differences.

Accordingly, the constraint violation degree $c_k$ at the sampling parameter $u_k$ is calculated as

$$\begin{array}{c}{c}_k={\displaystyle \sum _{i=1}^2}\max \left(0,g_i\left(u_k\right)\right),0\le k\le N\end{array}$$

(23)

According to the analysis of Equation (23), $c_k=0$ means that the motion state at the sampling parameter $u_k$ meets all constraints, otherwise one of the contour error constraints of cutter tip and cutter orientation is violated. Hence, when the objective function $F_{obj}$ is minimum and there is no violation of constraints, the corresponding control vertices ${\left\{d_i\right\}}_{i=0}^M$ can be used to generate the target optimal feedrate curve.

3.3.2. GSPSO Algorithm Based on Window Movement

Note that if all control vertices are simultaneously optimized, it has not only a problem in determining the optimal solution, but also needs to take a long time to calculate. To consider the accuracy and efficiency of the solution, the moving window planning method is adopted to attain the solution of the feedrate optimization model on the basis of the flexible local readjustment capability of the feedrate curve, as shown in Figure 3. Each time the window is moved, a group of control vertices are optimized until all control vertices are optimized. To ensure that the feedrate curve between the two adjacent windows can smoothly transition without violating the constraints of contour error in the optimization process, a certain width overlap region $R$ should be reserved between adjacent windows, as shown in Figure 4. Taking the control vertices contained in a single window as the optimization unit, the specific implementation process of GSPSO algorithm is introduced. According to the analysis results in Section 3.3.1, the search interval of control vertices in the GSPSO algorithm can be determined as $\left[0,f_{max}\right]$. To achieve the goal of optimizing the global search capability and convergence time, both the infeasible solutions of the optimization model should be sufficiently utilized, and the method of population initialization also needs to be optimized. The detailed implementation process of the GSPSO algorithm is described as follows:

Step 1. Set the initial and maximum values of the iteration times to 0 and $l_{max}$, respectively. The population velocity of random initialization is defined as ${\mathit{v}}^l=\left[{\mathit{v}}_1^l,{\mathit{v}}_2^l,\cdots ,{\mathit{v}}_s^l\right]$, where $s$ is the population size. Making the best of the boundary value $f_{max}$ of control vertices ${\left\{d_i\right\}}_{i=0}^M$ is beneficial to obtaining good initialization results. In view of the above, the initial position ${\mathit{x}}_m^l=\left[d_{r,m}^l,d_{r,m}^l,\cdots ,d_{t,m}^l\right]$ of the $m$th $h$-dimensional particle can be calculated as

$$\begin{array}{c}{\mathit{x}}_m^l=rand\left(1,h\right)f_{max},m=1,2,\dots ,s\end{array}$$

(24)

with

$$\begin{array}{c}h=t-r+1\end{array}$$

(25)

where $r$ and $t$, respectively, represent the sequence numbers of control vertices $d_r$ and $d_t$, and $rand\left(1,h\right)$ indicates a $1\times h$ dimensional belonging to the range of $\left[\mathrm{0,1}\right]$ uniform random number matrix.

Step 2. According to B-spline definition, the sampling parameter interval $\left[u_{r^{\mathrm{*}}},u_{t^{\mathrm{*}}}\right]$ that needs to detect the constraint violation degree $c_k$ can be determined by the particle $x_m^l$, and the corresponding total constraint violation degree $C_m^l$ is calculated as $C_m^l={\displaystyle \sum _{r=r^{\mathrm{*}}}^{t^{\mathrm{*}}}}{c}_k^l$. After that, the objective function corresponding to the particle ${\mathit{x}}_m^l$ is $F_{obj,m}^l$, which can be calculated by Equation (21). Set the initial individual best position ${\mathit{P}}_m^l={\mathit{x}}_m^l$, the individual best value $P_{best.m}^l=F_{obj,m}^l$, and the corresponding total constraint violation degree $P_{const,m}^l=C_m^l$ of the particle $x_m^l$. The initialization results of the global optimal position ${\mathit{G}}^l$, the global optimal value $G_{best}^l$ and the corresponding total constraint violation degree $G_{const}^l$ can be acquired from Algorithm 1.

Step 3. The linear weight reduction strategy [34] is adopted to determine the new position and velocity per particle. The updated process is as follows:

$$\begin{array}{c}{\mathit{v}}_m^{l+1}=\phi {\mathit{v}}_m^l{+c}_1{r}_1\left({\mathit{P}}_m^l-{\mathit{x}}_m^l\right)+c_2{r}_2\left({\mathit{G}}^l-{\mathit{x}}_m^l\right)\end{array}$$

(26)

$$\begin{array}{c}{\mathit{x}}_m^{l+1}={\mathit{x}}_m^l+{\mathit{v}}_m^{l+1}\end{array}$$

(27)

$$\begin{array}{c}\phi ={\phi }_{max}-\frac{\left({\phi }_{max}-{\phi }_{min}\right)}{L_{max}}l\end{array}$$

(28)

where ${\mathit{x}}_m^l$ and ${\mathit{v}}_m^l$ are the position and velocity of the $m$th particle in the $l$th iterations, respectively, ${\mathit{x}}_m^{l+1}$ and ${\mathit{v}}_m^{l+1}$ are the position and velocity of the $(m+1)$th particle in the $(l+1)$th iterations, respectively. $c_1$ and $c_2$ are the acceleration constant, also known as the learning factor, ${\phi }_{min}$ and ${\phi }_{max}$ are the upper and lower bounds of the inertia weight, respectively, $r_1$ and $r_2$ represent belonging to the range of $\left[\mathrm{0,1}\right]$ uniform random number.

Step 4. Handle boundary conditions and calculate the total constraint violation degree $C_m^l$ and the objective function value $F_{obj,m}^l$ of the particle ${\mathit{x}}_m^l$, and let $l=l+1$. Considering that infeasible solutions should be sufficiently utilized, the updated strategy [35] and Algorithm 1 are adopted to update the individual best position ${\mathit{P}}_m^l$ and global optimal position ${\mathit{G}}^l$.

Step 5. If $l<l_{max}$, return to step 3. Otherwise, the global optimal position $G_{best}^{l_{max}}$ and corresponding total constraint violation degree $G_{const}^{l_{max}}$ should be output.

Algorithm 1. Algorithm of global optimal position.
Input:	${\mathit{P}}_m^l$ $P_{best,m}^l$ $P_{const,m}^l$ $(m=\mathrm{1,2},\cdots ,s)$ $G^{l-1}$ $G_{best}^{l-1}$ $G_{const}^{l-1}$ $(l\ge 1)$ ${\mathit{G}}^l={\mathit{P}}_1^l$ $G_{best}^l=P_{best,1}^l$ $G_{const}^l=P_{const,1}^l$ $(l=0)$ $∕/$Initialization
Output:	${\mathit{G}}^l$ $G_{best}^l$ $G_{const}^l$
01:	for $m=1:s$ do
02:	if $G_{const}^{l-1}=0$ & $P_{const,m}^l=0$ $∕/$Both positions are feasible solutions
03:	if $P_{best}^l<G_{best}^{l-1}$
04:	${\mathit{G}}^l={\mathit{P}}_m^l$ $G_{best}^l=P_{best.m}^l$ $G_{const}^l=P_{const,m}^l$
05:	end if
06:	elseif $G_{const}^{l-1}=0$ & $P_{const,m}^l>0$ $∕/$One of the two solutions is feasible and the other is not.
07:	if $P_{best}^l<G_{best}^{l-1}$
08:	if $r_3<(0.5-l/(2l_{max}\left)\right)$ $∕/$Randomly generate $r_3\in \left[\mathrm{0,1}\right]$
09:	${\mathit{G}}^l={\mathit{P}}_m^l$ $G_{best}^l=P_{best.m}^l$ $G_{const}^l=P_{const,m}^l$
10:	end if
11:	end if
12:	elseif $G_{const}^{l-1}>0$ & $P_{const,m}^l=0$
13:	if $P_{best}^l<G_{best}^{l-1}$
14:	${\mathit{G}}^l={\mathit{P}}_m^l$ $G_{best}^l=P_{best.m}^l$ $G_{const}^l=P_{const,m}^l$
15:	else
16:	if $r_3>(0.5-l/(2l_{max}\left)\right)$ $∕/$Randomly generate $r_3\in \left[\mathrm{0,1}\right]$
17:	${\mathit{G}}^l={\mathit{P}}_m^l$ $G_{best}^l=P_{best.m}^l$ $G_{const}^l=P_{const,m}^l$
18:	end if
19:	end if
20:	elseif $G_{const}^{l-1}>0$ & $P_{const,m}^l>0$ $∕/$Both positions are infeasible solutions
21:	if $P_{best}^l<G_{best}^{l-1}$ $P_{const,m}^l<G_{const}^{l-1}$
22:	${\mathit{G}}^l={\mathit{P}}_m^l$ $G_{best}^l=P_{best.m}^l$ $G_{const}^l=P_{const,m}^l$
23:	elseif $P_{best.m}^l>G_{best}^{l-1}$ & $P_{const,m}^l<G_{const}^{l-1}$
24:	if $G_{const}^{l-1}/P_{const,m}^l>P_{best.m}^l/G_{best}^{l-1}$
25:	${\mathit{G}}^l={\mathit{P}}_m^l$ $G_{best}^l=P_{best.m}^l$ $G_{const}^l=P_{const,m}^l$
26:	end if
27:	elseif ${\mathit{P}}_{\mathit{b}\mathit{e}\mathit{s}\mathit{t}}^{\mathit{l}}<G_{best}^{l-1}$ & $P_{const,m}^l>G_{const}^{l-1}$
28:	if$G_{const}^{l-1}/P_{const,m}^l<P_{best.m}^l/G_{best}^{l-1}$
29:	${\mathit{G}}^l={\mathit{P}}_m^l$ $G_{best}^l=P_{best.m}^l$ $G_{const}^l=P_{const,m}^l$
30:	end if
31:	end if
32:	else
33:	${\mathit{G}}^l={\mathit{P}}_m^{l-1}$ $G_{best}^l=P_{best.m}^{l-1}$ $G_{const}^l=P_{const,m}^{l-1}$
34:	end if
35:	end for

3.4. Implementation Details of Feedrate Optimization

In light of the above modeling process and the proposed algorithm, detailed optimization steps of the feedrate curve are described as follows:

Step 1. The interpolation sampling parameter sequence ${\left\{u_k\right\}}_{k=0}^N$ is determined by Equation (2), on the basis of the given toolpath and the initial feedrate curve is generated by a random initialization of the GSPSO algorithm.

Step 2. For the sequence ${\left\{u_k\right\}}_{k=0}^N$, the corresponding reference command sequence $\left\{{\left[{\mathit{p}}_r\left(u_k\right),{\mathit{o}}_r\left(u_k\right)\right]}^T|0\le k\le N\right\}$ and $\left\{{\mathit{m}}_r\left(u_k\right)|0\le k\le N\right\}$ can be respectively obtained from Equations (3) and (4).

Step 3. According to the contour error prediction model in Section 3.1, the cutter tip contour error sequence $\left\{{\mathit{\epsilon }}_p\left(u_k\right)|0\le k\le N\right\}$ and the cutter orientation contour error sequence $\left\{{\mathit{\epsilon }}_0\left(u_k\right)|0\le k\le N\right\}$ are gained by the reference command sequence.

Step 4. If $l<l_{max}$, use the proposed GSPSO algorithm to optimize the previous feedrate curve and return to Step 1. Otherwise, output the final feedrate curve of the current window.

Step 5. Repeat the above process as the window moves until the optimization of the global feedrate curve is completed.

Figure 4 and Figure 5, respectively, present a flow chart of feedrate curve optimization in the current window and a diagram of the window moving method to better understand the optimization process of the global feedrate curve.

4. Experimental Verification

The effectiveness of the proposed feedrate optimization algorithm is validated by three typical toolpaths in this section. The A–C double-turntable open five-axis machine tool for experimental verification is shown in Figure 6. During actual processing, the actual position information of the linear axes and rotary axes of the machine tool can be acquired from the optical encoder with a resolution of $500\times 256$ pulses/rev. The motor and driver driving linear axes and rotary axes are of the YASKAWA brand. The time constants of the feed axes of the machine tool were determined using the method proposed by Dong [32], which were approximated as $T_X=0.0231$, $T_Y=0.0231$, $T_Z=0.0271$, $T_A=0.0262$ and $T_Y=0.0215$, respectively. The interpolation sampling period $T_s$ is 4 mms. To consider the calculation efficiency and optimization result accuracy, the total number of control vertices of the feedrate curve $M$ is 120, the interval window interval $h$ is 25, and the window overlap region $R$ is 5. Parameters of GSPSO algorithm are shown in

Table 1. Parameters of GSPSO algorithm.

${\mathit{c}}_1$	${\mathit{c}}_2$	${\mathit{\phi }}_{\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{\phi }}_{\mathit{m}\mathit{i}\mathit{n}}$	$\mathit{s}$	${\mathit{l}}_{\mathit{m}\mathit{a}\mathit{x}}$
1.4962	1.4962	0.8	0.4	100	30

. Moreover, considering the driving capacity limitation of the YASKAWA servo driver used in the experiment, the maximum programmed feedrate $f_{max}$ of the cutter tip was set to $20 mm/s$ during the actual experiment.

4.1. Example 1

The total length of the butterfly trajectory expressed by NURBS is 338.68 mm in Experiment 1, as shown in Figure 7a. According to the analysis of trajectory characteristics, many high-curvature areas on the toolpath are conducive to verifying the excellent performance of the algorithm. In this example, the maximum limit value of the cutter tip contour error ${\mathit{\epsilon }}_{pmax}$ is 0.05 mm. Figure 7b shows the feedrate curve with a cutter contour error constraint generated by the proposed method. The actual cutter tip contour error during the experiment is shown in Figure 7c. Analysis of this picture reveals that the cutter tip contour error at each point on the toolpath does not go beyond the set constraint limit, which proves that the algorithm is effective. Further study of the experimental results is able to find that the cutter tip error almost reaches the set maximum error limit in multiple processing time intervals, which also reflects the excellent performance of the algorithm to some extent.

4.2. Example 2

Unlike three-axis machining, the nonlinear coupling motion peculiarity of five-axis machining makes it more difficult to accurately control the cutter tip contour error. In Experiment 2, take the five-axis end milling trajectory expressed by NURBS as shown in Figure 8a as an example for experimental verification, and the toolpath length is 128.45 mm. Figure 8b,c, respectively, represent the corresponding reference position curves of the translational axes $(X,Y,Z)$ and the rotational axes $(A,C)$ in the MCS. The maximum limit value of the cutter tip contour error ${\mathit{\epsilon }}_{pmax}$ is 0.04 mm in this example. The five-axis milling process and the machined part with the optimized feedrate are shown in Figure 8d,e. Figure 8f presents the feedrate curve with cutter contour error constraint generated by the proposed method. Figure 8g displays the experimental measurement results of cutter tip contour error for five-axis end milling. According to the analysis of the experimental results, the cutter contour error does not exceed the preset error constraint range and fluctuates around the limit value of the constraint in most of the processing time, which further exhibits the superiority of the algorithm.

4.3. Example 3

Compared with end milling, flank milling needs to bound the contour error of cutter tip and cutter orientation. To guarantee the contour accuracy of flank milling, the contour error of cutter tip and cutter orientation should be constrained simultaneously. Figure 9a exhibits the five-axis flank milling trajectory expressed by dual NURBS that is taken as the test object for experimental verification, and the toolpath length is 329.44 mm. Figure 9b,c indicate the corresponding reference position curves of the translational axes $(X,Y,Z)$ and the rotational axes $(A,C)$ in the MCS, respectively. To further validate that the proposed feedrate optimization algorithm is effective and excellent, the contour accuracy and processing time of the constant feedrate process of machine tools and the new method are compared. In this comparative experiment, the limit value of the contour error of the cutter tip ${\mathit{\epsilon }}_{pmax}$ is 0.04 mm and the limit value of the contour error of cutter orientation ${\mathit{\epsilon }}_{omax}$ is $3\times {10}^{-4}$ rad, and the constant feedrate process parameter is 10mm/s. Based on the established feedrate optimization model, the feedrate curve is determined under the constraint of contour error, including the cutter tip (colored in pink), the cutter orientation (colored in cyan) and the final blue feedrate curve, as shown in Figure 9d.

Figure 9e,f present the comparative experimental results of the contour error of cutter tip and cutter orientation for using the above-mentioned two methods. On the basis of the analysis of experimental data, the maximum values of the contour error cutter tip and cutter orientation are 0.0692 mm and $6.019\times {10}^{-4}$ rad for using constant feedrate processing, respectively, which are out of limits and do not meet the requirements of machining accuracy. In contrast, the paper proposes a feedrate optimization algorithm that can not only effectively control the contour error cutter tip and cutter orientation within the preset range but can also greatly improve the maximum feasible feedrate. According to the calculation of specific experimental data, the maximum values of the contour error cutter tip and cutter orientation have decreased by 42.2% and 50.16%, respectively, and the machining time has shortened by 11.42%. In addition, further analysis of Figure 9e,f demonstrates that at least one out of the cutter tip contour error and cutter orientation contour error has reached its upper bound in the whole machining process, which reflects the optimization of the proposed algorithm in great measure.

Note that the interpolation sampling position sequence of each servo feed axis of the machine tool is generated without going against the drive constraints in all experiments. Therefore, the influence of other factors besides contour error can be minimized or ignored in the process of experimental verification. As a matter of fact, based on the process of solving the feedrate optimization model with the proposed GSPSO algorithm, it is easy and reliable to introduce more constraints related to machining into the existing optimization model. The analysis of the contour error image indicates that although the relatively satisfactory experimental results have been achieved, there is still some optimization space to acquire more accurate control. It is speculated that this may be mainly caused by the loss of a certain accuracy of the tracking error prediction model due to the transfer function of the servo system being simplified. In the next work, the dynamic response, particularity of the servo system, will be further studied and a more accurate tracking error prediction model will also be established.

5. Conclusions

The heuristic feedrate optimization algorithm is proposed to optimize five-axis machining in this paper and is capable of achieving the accurate control of five-axis contour error. Firstly, according to the established tracking error dynamic prediction model and the contour error iterative estimation model, the relationship between the parametric feedrate and contour error constraint is clarified, which provides a model basis for the accurate control of contour error by optimizing the feedrate curve. After that, the feedrate optimization model aiming at minimizing the machining time and taking the control vertices of the feedrate curve expressed by cubic B-spline as decision variables is established. Subsequently, considering that B-spline has a flexible local readjustment capability, the GSPSO algorithm based on window movement is applied to optimize the feedrate curve in segments, which solves the problem of low accuracy and low efficiency caused by the single optimization of all control vertices. Finally, experimental tests on an open five-axis machine tool validate the excellent performance of the proposed feedrate optimization algorithm. Experimental results show that the proposed algorithm is capable of sufficiently releasing the potential of the machine tools while accurately controlling the contour error of the cutter tip and cutter orientation. Moreover, it is easy to implement and trustworthy, the size of contour error constraints can be adjusted with freedom and more constraints can be introduced in accordance with the actual machining requirements.