Workpiece Pose Optimization for Milling with Flexible-Joint Robots to Improve Quasi-Static Performance

: Although industrial robots are widely used in production automation, their applications in machining have been limited because of the structural vibrations induced by periodic cutting forces. Since the dynamic characteristics of an industrial robot depends on its conﬁguration, the responses of the robot structure to the cutting forces are affected by how the workpiece is placed within the workspace of the robot. This paper presents a method for workpiece pose optimization for a robotic milling system to improve the quasi-static performance during machining. Since the milling forces are time-varying due to the characteristics of the multi-tooth and discontinuity of milling, these forces can excite vibrations inside the robot structure. To address this issue, a structural dynamics model is established for industrial robots, considering their joint ﬂexibility, and a milling force formulation is incorporated into the robot dynamics model to investigate the forced vibrations of the ﬂexible joints. Then, the quasi-static performance of the robotic machining system is evaluated by the vibration-induced offset of the cutter tool that is mounted on the end-effector. Finally, an optimization approach is given for the workpiece pose to minimize the cutter tool offset under the periodic milling force. A numerical simulation demonstrates that the optimal workpiece pose can signiﬁcantly reduce the overall tool offset during machining and can lower the variation of the tool offset along the milling path.


Introduction
Due to their flexibility, efficiency, and reliability, industrial robots have played an important role in production automation for a wide spectrum of applications including welding, painting, assembly, packaging, and material handling [1]. Recently, these machines have also been proposed for more complex tasks with a large dynamic loading such as machining [2], and these robotic applications require a high positioning accuracy and involve large time-varying loading during the processes.
A few examples are the European projects HEPHESTOS and COMET that were launched to develop innovative robotic machining systems for small-batch production [3,4] and Boeing's initiative that aimed at using robots to automate drilling and percussive riveting processes in aircraft assembly [5]. On the other hand, the challenges for these robotic applications have also been highlighted in literature due to the inherent disadvantages of industrial robots compared with machine tools, such as having a low positioning accuracy and being prone to vibrations [2].
With the development of modern metrology techniques such as laser tracking [6,7] and photogrammetry-based measurement [8,9], the robotic positioning accuracy can be improved by Appl. Sci. 2019, 9, 1044; doi:10.3390/app9061044 www.mdpi.com/journal/applsci kinematic parameter calibration [10][11][12] and positioning error compensation [7,8,[13][14][15][16]. For example, Andres et al. [12] developed a calibration method to estimate the external joint configuration parameters using a laser displacement sensor. Morris et al. [16] developed a kinematic model for a robot arm through a coordinate measuring machine (CMM), and experimental measurements at some robot poses were taken using the CMM. However, industrial robots still face issues due to structural vibrations for these applications, and the resulting influence is difficult to eliminate by metrology techniques. Robotic machining involves dynamic loading in a way that the magnitudes and directions of the force and torque acting on the robot dramatically change with time during machining. The resulting forced vibrations of the robot structure lead to much more complicated dynamics behaviors than many other robotic tasks such as pick-and-place and welding where static or very low contact forces act on the robot. First, the forced vibrations of the robot can seriously affect the process, for example, by decreasing the efficiency [17], damaging the surface quality [18], causing tool misalignment [19], and reducing the fatigue lives of the tool and the robot to a great extent [20,21]. Second, not only static but also inertial forces are important for these applications. Furthermore, the interaction between the machining forces and the robot vibrations can even cause instability to the robotic machining system [22,23] if some machining parameters such as spindle speed and feed rate are not selected properly [24,25]. All these works have shed some light on the complexity of structural dynamics of robotic machining.
Several approaches have been proposed to improve the robotic machining quality for a specific task considering the forced vibrations caused by machining. Vosniakos et al. [26] used two genetic algorithms to find the region with the best manipulability and the lowest joint torque in the robot workspace. Baglioni et al. [27] developed a parametric modeling procedure that can correctly simulate the robot dynamics behavior. Klimchik et al. [28] modeled the interaction between the milling tool and the workpiece and proposed a method for compensating the compliance errors based on a nonlinear stiffness model. Kaldestad et al. [29,30] adjusted the tool path to counteract the tool deflections during the milling operations after estimating the tool force and robot joint stiffness. Kubela et al. [31] proposed an online method for compensating the backlash errors induced by the drive reversion. Cordes et al. [32] proposed a new joint space model that combined the joint stiffness and the reversal error and then used the model to compensate the milling path deviation. Shi et al. [33] proposed a prediction model that could estimate the stable and unstable zones of milling process.
In a robotic machining system, the workpiece is usually clamped on a work table and must be placed within the work space of the robot. The configuration of the robot and the cutting force acting on it during machining are dependent upon the pose of the workpiece relative to the robot. It has been found that the natural frequencies of a robot strongly depend on its configuration [34] and that the robot stiffness varies significantly in different directions [35,36] and can be improved by optimizing the robot configuration for a given pick-and-place task [36][37][38], which implies that the vibration characteristics of robotic machining are affected by how the workpiece is placed with respect to the robot. For example, the robot stiffness is considered in Reference [37] for the optimization of the redundant degrees of freedom to minimize the displacement of the cutting tool under static force. In this paper, the forced vibrations of the robot are considered and the optimization goal is to minimize the overall steady-state offset amplitude of the cutting tool under the time-varying cutting force during the milling process. To this end, a structural dynamics model is established for industrial robots considering their joint flexibility, and a milling force formulation is incorporated into the robot dynamics model to investigate the forced vibrations of the flexible joints. Then, the quasi-static performance of the robotic machining system is evaluated by the vibration-induced offset of the cutter tool that is mounted on the end-effector. The performance is quasi-static because the effect of the robot's deformation on the milling force is neglected. Finally, an optimization approach is given for the workpiece pose to minimize the cutter tool offset under the periodic milling force.
The next three sections of this work are organized as follows. Section 2 describes the robotic milling system and the problem statement. Section 3 presents the theoretical modeling, including the robot dynamics model, the milling force model, and the optimization approach for the workpiece pose. A numerical simulation study is given in Section 4.

System Description and Problem Statement
As shown in Figure 1, a robotic machining system for milling is being developed. The system is composed of a Kuka KR 300-2 PA robot, a cylindrical helical end milling cutter mounted the robot end-effector, a work table located within the robot workspace, a modular fixture placed on the work table, and a workpiece to be machined. Both the work table and the fixture have pin holes for positioning reference and bolt holes for fastening, and the workpiece is fixed on the fixture. In addition, the inclination angle θ p between the fixture and the work table can be adjusted. As a result, the pose of the workpiece relative to the robot can be represented by the position of the fixture on the work table denoted by x p and y p and the inclination angle θ p as displayed in Figure 1.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 3 of 15 robot dynamics model, the milling force model, and the optimization approach for the workpiece pose. A numerical simulation study is given in Section 4.

System Description and Problem Statement
As shown in Figure 1, a robotic machining system for milling is being developed. The system is composed of a Kuka KR 300-2 PA robot, a cylindrical helical end milling cutter mounted the robot end-effector, a work table located within the robot workspace, a modular fixture placed on the work table, and a workpiece to be machined. Both the work table and the fixture have pin holes for positioning reference and bolt holes for fastening, and the workpiece is fixed on the fixture. In addition, the inclination angle θp between the fixture and the work table can be adjusted. As a result, the pose of the workpiece relative to the robot can be represented by the position of the fixture on the work table denoted by xp and yp and the inclination angle θp as displayed in Figure 1. During the milling operation, the robot is subjected to periodic milling forces due to the intermittence of the cutter teeth. Since the machining operations can excite vibrations inside the robot structure, it is necessary to guarantee that the offset of the cutter tool maintains certain tolerances during the milling operation. In other words, the robot must be able to the hold the cutter tool firmly under the periodic milling forces. The machining forces can be changed by adjusting the machining parameters like the spindle speed and the feed rate; however, the adjustment of these parameters is constrained by many other factors such as the process efficiency requirement and the capability of the spindle and the robot. On the other hand, as mentioned above, the vibration characteristics of the system are also influenced by the workpiece pose relative to the robot, and the workpiece pose can be adjusted relatively easily. Then, the problem under study can be defined given the required milling tasks and milling parameters to find the optimal workpiece pose relative to the robot that minimizes the cutter tool offset under the periodic milling forces.

Dynamic Model of Robot
We first establish a structural dynamics model for the industrial robot to analyze the vibration characteristics. For this purpose, an n-DOF industrial robot with n links and n rotary joints is illustrated in Figure 2. A cutter tool is mounted on the robot end-effector. The local body-fixed During the milling operation, the robot is subjected to periodic milling forces due to the intermittence of the cutter teeth. Since the machining operations can excite vibrations inside the robot structure, it is necessary to guarantee that the offset of the cutter tool maintains certain tolerances during the milling operation. In other words, the robot must be able to the hold the cutter tool firmly under the periodic milling forces. The machining forces can be changed by adjusting the machining parameters like the spindle speed and the feed rate; however, the adjustment of these parameters is constrained by many other factors such as the process efficiency requirement and the capability of the spindle and the robot. On the other hand, as mentioned above, the vibration characteristics of the system are also influenced by the workpiece pose relative to the robot, and the workpiece pose can be adjusted relatively easily. Then, the problem under study can be defined given the required milling tasks and milling parameters to find the optimal workpiece pose relative to the robot that minimizes the cutter tool offset under the periodic milling forces.

Dynamic Model of Robot
We first establish a structural dynamics model for the industrial robot to analyze the vibration characteristics. For this purpose, an n-DOF industrial robot with n links and n rotary joints is illustrated in Figure 2. A cutter tool is mounted on the robot end-effector. The local body-fixed coordinate systems {o i _x i y i z i } (i = 1, 2, . . . ,n) and the tool coordinate system {o t _x t y t z t } are established at each joint and at the cutter tool tip, respectively. The position vector of the ith joint and the origin of tool coordinate system in the inertial coordinate system {o_xyz} are denoted by P i and P t respectively. During the milling operation, the cutter is subjected to an external machining force F E .
Appl. Sci. 2019, 9, x FOR PEER REVIEW 4 of 15 coordinate systems {oi_xiyizi} (i =1, 2,…,n) and the tool coordinate system {ot_xtytzt} are established at each joint and at the cutter tool tip, respectively. The position vector of the ith joint and the origin of tool coordinate system in the inertial coordinate system {o_xyz} are denoted by Pi and Pt respectively. During the milling operation, the cutter is subjected to an external machining force FE. For industrial robots, the main source of deformation is joint flexibility [1,39]. Thus, to simplify the robot model, we assume that all the links are rigid and that all the joints are equivalent to linear elastic torsion springs. If the joint angles for a given configuration are denoted as q = (q1, q2, …, qn) T , then the dynamic model of the stationary robot under forced vibrations in this configuration can be written as where M represents the n×n symmetric generalized mass matrix, K is the n×n diagonal stiffness matrix with its diagonal entries equal to the torsional stiffness of the corresponding joints, and C is a damping matrix. Both K and C are defined in the joint space. Vector Δq is the perturbations in the joint angles from an equilibrium configuration and represents the joint deflections. The 3×n linear velocity Jacobian matrix Jv represents the mapping from the joint rates to the velocity of the tool tip such that v = Jv , in which v is the linear velocity of tool tip written in the inertial coordinate system and is the joint rates. Matrices M and Jv depend on the robot configuration, and the force FE depends on the specific machining task.
The dynamic model described as Equation (1) is a linear time invariant system. Force FE can be considered as the external input, and Δq is the output response of the robot. In addition, if the joint deflection is small, the linear displacement of tool tip offset can be written as Thus, once the joint deflection is obtained, the offset of the cutter tool caused by the robot vibration can be calculated through Equation (2), which can be used to evaluate the quasi-static performance of the robotic machining system under a specific machining force.

Milling Force Model
A cutting force model for cylindrical helical end milling is incorporated into the equation of motion of the robot to obtain the force vector FE, without considering the regenerative effect and the effect of the dynamic response of the robot on the cutting force [40,41]. In our system, up milling is chosen due to the concern with a backlash of the system. The readers are referred to Reference [42]  For industrial robots, the main source of deformation is joint flexibility [1,39]. Thus, to simplify the robot model, we assume that all the links are rigid and that all the joints are equivalent to linear elastic torsion springs. If the joint angles for a given configuration are denoted as q = (q 1 , q 2 , . . . , q n ) T , then the dynamic model of the stationary robot under forced vibrations in this configuration can be written as M∆ ..
where M represents the n × n symmetric generalized mass matrix, K is the n × n diagonal stiffness matrix with its diagonal entries equal to the torsional stiffness of the corresponding joints, and C is a damping matrix. Both K and C are defined in the joint space. Vector ∆q is the perturbations in the joint angles from an equilibrium configuration and represents the joint deflections. The 3 × n linear velocity Jacobian matrix J v represents the mapping from the joint rates to the velocity of the tool tip such that v = J vq , in which v is the linear velocity of tool tip written in the inertial coordinate system andq is the joint rates. Matrices M and J v depend on the robot configuration, and the force F E depends on the specific machining task. The dynamic model described as Equation (1) is a linear time invariant system. Force F E can be considered as the external input, and ∆q is the output response of the robot. In addition, if the joint deflection is small, the linear displacement of tool tip offset can be written as Thus, once the joint deflection is obtained, the offset of the cutter tool caused by the robot vibration can be calculated through Equation (2), which can be used to evaluate the quasi-static performance of the robotic machining system under a specific machining force.

Milling Force Model
A cutting force model for cylindrical helical end milling is incorporated into the equation of motion of the robot to obtain the force vector F E , without considering the regenerative effect and the effect of the dynamic response of the robot on the cutting force [40,41]. In our system, up milling is chosen due to the concern with a backlash of the system. The readers are referred to Reference [42] for more detail about the cutting force modeling of down milling, which is similar in principle to that of up milling. In this cutting force model, the instantaneous cutting forces are assumed to be proportional to the area of contact between the flute and the workpiece [43]. Then, a mathematical model of the instantaneous milling force can be constructed by calculating the uncut chip thickness in a certain time. Then, the change of the cutter edge length during milling is analyzed to obtain the continuous milling force in the whole period. Figure 3 shows the geometric size and force condition of a single cutting edge of the cylindrical helical end mill. The tool coordinate system {o t _x t y t z t } (also see Figure 2) is established as follows. The origin is at the center of the milling cutter tool tip, the z t axis is along the centerline of the milling cutter, the x t axis is along the direction of feed rate, and the y t axis is determined from the right-hand rule. N represents the number of milling teeth; β and R represent the helical angle and the radius of the cutter, respectively; and a p and a e represent the axial and radial depth of cut, respectively. for more detail about the cutting force modeling of down milling, which is similar in principle to that of up milling. In this cutting force model, the instantaneous cutting forces are assumed to be proportional to the area of contact between the flute and the workpiece [43]. Then, a mathematical model of the instantaneous milling force can be constructed by calculating the uncut chip thickness in a certain time. Then, the change of the cutter edge length during milling is analyzed to obtain the continuous milling force in the whole period. Figure 3 shows the geometric size and force condition of a single cutting edge of the cylindrical helical end mill. The tool coordinate system {ot_xtytzt} (also see Figure 2) is established as follows. The origin is at the center of the milling cutter tool tip, the zt axis is along the centerline of the milling cutter, the xt axis is along the direction of feed rate, and the yt axis is determined from the right-hand rule. N represents the number of milling teeth; β and R represent the helical angle and the radius of the cutter, respectively; and ap and ae represent the axial and radial depth of cut, respectively. Due to the existence of the helix angle of the cutter tool, the axial height of cutter teeth varies with the cutting thickness. Therefore, a microelement method is used to model a single cutter tooth. First, the cutter edge is discretized into numerous cutting-edge elements along the zt axis direction, and the height of each axial element is denoted as dz. φ represents the projection angle of the entire cutting edge in the xtotyt plane, and this projection angle is denoted as ψ when the axial height reaches z'. Then we have the following equations [42][43][44]: where nc represents the spindle speed with a unit of rpm and t represents time. The cutting-edge angle θ between the microelement with an axial height of z' and the yt axis can be calculated as The determination of the cutting forces requires the area of contact between the end mill and the workpiece. For an axial element of thickness dz, this is equivalent to requiring the uncut chip thickness h. The uncut chip thickness h with an axial height of z' can be approximately expressed as where ft = vf /(nN) represents the feed per tooth and vf represents the feed rate. Because dz is very small, the area of contact with a height of dz between the flute and the workpiece at time t can be Due to the existence of the helix angle of the cutter tool, the axial height of cutter teeth varies with the cutting thickness. Therefore, a microelement method is used to model a single cutter tooth. First, the cutter edge is discretized into numerous cutting-edge elements along the z t axis direction, and the height of each axial element is denoted as d z . ϕ represents the projection angle of the entire cutting edge in the x t o t y t plane, and this projection angle is denoted as ψ when the axial height reaches z . Then we have the following equations [42][43][44]: where n c represents the spindle speed with a unit of rpm and t represents time. The cutting-edge angle θ between the microelement with an axial height of z and the y t axis can be calculated as The determination of the cutting forces requires the area of contact between the end mill and the workpiece. For an axial element of thickness dz, this is equivalent to requiring the uncut chip thickness h. The uncut chip thickness h with an axial height of z can be approximately expressed as where f t = v f /(nN) represents the feed per tooth and v f represents the feed rate. Because dz is very small, the area of contact with a height of dz between the flute and the workpiece at time t can be approximately equal to h(θ)dz. Then, we can calculate the elemental cutting force at this moment. Next, the force generated by this cutting tooth can be obtained by adding up all the elemental cutting forces. This is a definite integral problem, and we can get the result by integrating the axial height numerically. The upper and lower limits of integration can be obtained from the height of the cutting tooth involved during milling. In addition, there can be multiple cutting teeth involved in the cutting at time t. Then, we need to add up the forces from all the cutting teeth involved in the cutting to obtain the overall cutting force. The milling force is periodic because of the intermittence of the cutter teeth. The axial milling force along the z t axis can be neglected, since it is typically much smaller than the tangential and radial forces and has small contributions to the bending moment acting on the cutter tool. The detailed expressions for the tangential and radial milling forces are elaborated in thte Appendix A.

Steady-State Response and Optimization Method
Since the robot drives the tool to move slowly along the milling path, we can use the steady-state response of the cutter tool offset to evaluate the quasi-static performance of robotic milling in a given position. The cutter tool offset can be solved from Equations (1) and (2) in the time domain if the external forces are given. For a specific milling task, we select n p points along the path and obtain the amplitudes d i (i = 1, 2, . . . , n p ) of the steady-state responses of the cutter tool offset at these points. The steady-state response amplitudes d i can also be obtained from Equation (1). Then, we calculate the mean value of all these amplitudes d m to quantify the overall quasi-static performance for this task. Then we have the following Equation: If a local coordinate system {o w _x w y w z w } is established on the workpiece, the position vectors p mi (i = 1, 2, . . . ,n p ) of the selected points and the milling forces F Ei (I = 1, 2, . . . , n p ) exerted on the cutter at these points in the inertial coordinate system can be expressed as where p c and p w represent the position vectors of the work table and workpiece coordinate systems respectively and both vectors are written in the inertial coordinate system. In particular, vector p w can be obtained by x p and y p as p w = [x p , y p , z w ] T in which z w is a constant and represents the height of the origin of workpiece coordinate system o w . Vectors p ni (i = 1, 2, . . . ,n p ) represent the position vectors of the selected points in the inertial coordinate system and can be calculated as pn = R wz (θ p )p' ni in which p' ni are the position vectors of the selected points in the workpiece coordinate system and R wz (θ p ) represents the rotation matrix from the workpiece coordinate system to the inertial coordinate system. Vector F' represents the cutting force in the tool coordinate system and can be obtained with the cutting force model presented in Section 3.2, and the 3 × 3 matrices R ni (i = 1, 2 . . . n p ) represent the rotation matrices from the tool coordinate system of the selected points to the inertial coordinate system. Now, the process to calculate d m is illustrated in Figure 4.  The workpiece has three degrees of freedom relative to the work table in this robotic milling system, which can be denoted by a parameter vector p = [xp, yp, θp] T . Then we aim at minimizing the mean value of the amplitudes of the cutter tool offset along the machining path by adjusting these parameters. The optimization problem can be expressed as where dm is computed with the steps described in Figure 4. Once this optimization problem is solved, we find the optimal pose of the workpiece relative to the robot for the milling.

System Parameters
In this simulation analysis, we only consider the flexibility of the first three joints of the robot because the moment arms of the cutting force about the last three joints are much shorter than those about the first three joints. The aforementioned coordinate systems in this simulation are described as follows. The inertial coordinate system {o_xyz} and the local coordinate systems {oi_xiyizi} (i = 1, 2, 3) of the robot are established as shown in Figure 5a. The inertial coordinate system is located at the base of the robot with the positive z-axis vertically upward and the y-axis perpendicular to the plane containing the second and third links. The zi-axis of the local coordinate systems are along the rotation directions of the corresponding joints; the axis x1 is parallel to the axis x when the robot is stationary, and the x2 and x3 axes are along the link length directions. The directions of the other reference axes of {oi_xiyizi} (i = 1, 2, 3) are determined from the right-hand rule. The local coordinate system of the work table {oc_xcyczc} is established at the middle position near the robot, with the xc and zc axes parallel to the x and z axes respectively and the yc axis is determined from the right-hand rule. The workpiece and the milling path are shown in Figure 5b. The workpiece coordinate system {ow_xwywzw} is established at the bottom of the workpiece, with the xw and yw axes parallel to the length and width and also parallel to the xc and yc axes respectively and the zw axis is determined from the right-hand rule. The workpiece has three degrees of freedom relative to the work table in this robotic milling system, which can be denoted by a parameter vector p = [x p , y p , θ p ] T . Then we aim at minimizing the mean value of the amplitudes of the cutter tool offset along the machining path by adjusting these parameters. The optimization problem can be expressed as where d m is computed with the steps described in Figure 4. Once this optimization problem is solved, we find the optimal pose of the workpiece relative to the robot for the milling.

System Parameters
In this simulation analysis, we only consider the flexibility of the first three joints of the robot because the moment arms of the cutting force about the last three joints are much shorter than those about the first three joints. The aforementioned coordinate systems in this simulation are described as follows. The inertial coordinate system {o_xyz} and the local coordinate systems {o i _x i y i z i } (i = 1, 2, 3) of the robot are established as shown in Figure 5a. The inertial coordinate system is located at the base of the robot with the positive z-axis vertically upward and the y-axis perpendicular to the plane containing the second and third links. The z i -axis of the local coordinate systems are along the rotation directions of the corresponding joints; the axis x 1 is parallel to the axis x when the robot is stationary, and the x 2 and x 3 axes are along the link length directions. The directions of the other reference axes of {o i _x i y i z i } (i = 1, 2, 3) are determined from the right-hand rule. The local coordinate system of the work table {o c _x c y c z c } is established at the middle position near the robot, with the x c and z c axes parallel to the x and z axes respectively and the y c axis is determined from the right-hand rule. The workpiece and the milling path are shown in Figure 5b. The workpiece coordinate system {o w _x w y w z w } is established at the bottom of the workpiece, with the x w and y w axes parallel to the length and width and also parallel to the x c and y c axes respectively and the z w axis is determined from the right-hand rule.  Table 1. Note that these vectors and tensors are written in the corresponding local coordinate systems. The torsional stiffnesses and damping of the three joints can be obtained by a stiffness identification experiment [45,46]. However, our system is still being built, and thus, the joint stiffness values available for a robot with similar payload and size are adopted, i.e., k1 = 2.4 × 10 5 Nm/rad and k2 = k3 = 10 5 Nm/rad. A damping ratio ξ of 0.06 is used to compute the damping matrix based on the modal test data for industrial robots given in Reference [46]. Then, the ranges of the joint angles are θ1 є [−π, π], θ2 є [−π/2, π/2], and θ3 є [−π/3, π/3].  The parameters for the cutter tool are described below. The teeth number and helical angle of the milling cutter are 2 and π/3 in this simulation. The diameter of the tool is 0.024 m, and the material of the tool and workpiece are uncoated tungsten carbide and aluminum alloy, respectively. The axial cutting depth and radial cutting depth are 0.002 m and 0.004 m, respectively, and the feed per tooth ft is 0.00012 m. The coefficient parameters kt1 and kt2 for calculating the cutting force (see Appendix)  Table 1. Note that these vectors and tensors are written in the corresponding local coordinate systems. The torsional stiffnesses and damping of the three joints can be obtained by a stiffness identification experiment [45,46]. However, our system is still being built, and thus, the joint stiffness values available for a robot with similar payload and size are adopted, i.e., k 1 = 2.4 × 10 5 Nm/rad and k 2 = k 3 = 10 5 Nm/rad. A damping ratio ξ of 0.06 is used to compute the damping matrix based on the modal test data for industrial robots given in Reference [46]. Then, the ranges of the joint angles are θ 1 ∈ [−π, π], θ 2 ∈ [−π/2, π/2], and θ 3 ∈ [−π/3, π/3].  The parameters for the cutter tool are described below. The teeth number and helical angle of the milling cutter are 2 and π/3 in this simulation. The diameter of the tool is 0.024 m, and the material of the tool and workpiece are uncoated tungsten carbide and aluminum alloy, respectively. The axial cutting depth and radial cutting depth are 0.002 m and 0.004 m, respectively, and the feed per tooth f t is 0.00012 m. The coefficient parameters k t1 and k t2 for calculating the cutting force (see Appendix A) are experimentally determined as 387 and 0.0018, respectively. Similarly, parameters b 1 and b 2 are determined as −0.327 and −0.224, respectively [43]. The readers are referred to the Appendix A for more information about these cutting force parameters.

Milling Force and Mean Value of Tool Offset
According to the milling parameters given in Section 4.1, we first obtain the milling force and its Fast Fourier Transform (FFT), as shown in Figure 6, when the speed n is equal to 1000 rpm. We can see that this spindle speed leads to a frequency of 33.33 Hz for the cutting force with two teeth.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 15 are experimentally determined as 387 and 0.0018, respectively. Similarly, parameters b1 and b2 are determined as −0.327 and −0.224, respectively [43]. The readers are referred to the Appendix for more information about these cutting force parameters.

Milling Force and Mean Value of Tool Offset
According to the milling parameters given in Section 4.1, we first obtain the milling force and its Fast Fourier Transform (FFT), as shown in Figure 6, when the speed n is equal to 1000 rpm. We can see that this spindle speed leads to a frequency of 33.33 Hz for the cutting force with two teeth. FFT   In this simulation, we use the lsim function in MATLAB to compute the time-domain response of the joint deflections Δq under the periodic milling force from Equation (1). Then, the cutter tool offset can be obtained from Equation (2) for every point along the tool path. Figure 7 shows the mean value of the tool offset amplitudes dm with different xp, yp, and θp and the variation of dm in the xpyp plane with different θp. The following observations can be found. In Figure 7a, when yp is equal to 0 m and 0.25 m, the deflection of joint 3 plays a major role when xp is between 0 m and 0.2 m. When xp increases, the torque acting on joint 3 due to the x component of the cutting force decreases. When xp is between 0.2 m and 0.6 m, the deformation of joint 1 plays a major role, and the torque acting on joint 1 due to the y component of the cutting force increases with the increase of xp. When yp is equal to 0.5 m, the deformation of joint 1 plays a major role; therefore, we see that dm increases when xp increases. In Figure 7b, it is seen that dm is asymmetrical about the oxz plane because the cutting forces are in opposite directions in every two symmetrical positions about the oxz plane along the milling groove. Figure 7c shows that dm increases when the inclination angle θp increases. This is likely because the z component of milling force and the resulting torques on joints 2 and 3 increase as the angle θp increases. Similarly, we find in Figure 7d that as angle θp increases, dm becomes larger. In addition, it can be seen from Figure 7d that when θp remains unchanged, dm decreases as the workpiece is placed closer to the robot. In particular, when vector p is equal to [0.6 m, 0.5 m, π/3] T and [0.6 m, −0.5 m, π/3] T , dm dramatically increases since the natural frequencies of the In this simulation, we use the lsim function in MATLAB to compute the time-domain response of the joint deflections ∆q under the periodic milling force from Equation (1). Then, the cutter tool offset can be obtained from Equation (2) for every point along the tool path. Figure 7 shows the mean value of the tool offset amplitudes d m with different x p , y p , and θ p and the variation of d m in the x p y p plane with different θ p . The following observations can be found. In Figure 7a, when y p is equal to 0 m and 0.25 m, the deflection of joint 3 plays a major role when x p is between 0 m and 0.2 m. When x p increases, the torque acting on joint 3 due to the x component of the cutting force decreases. When x p is between 0.2 m and 0.6 m, the deformation of joint 1 plays a major role, and the torque acting on joint 1 due to the y component of the cutting force increases with the increase of x p . When y p is equal to 0.5 m, the deformation of joint 1 plays a major role; therefore, we see that d m increases when x p increases. In Figure 7b, it is seen that d m is asymmetrical about the oxz plane because the cutting forces are in opposite directions in every two symmetrical positions about the oxz plane along the milling groove. Figure 7c shows that d m increases when the inclination angle θ p increases. This is likely because the z component of milling force and the resulting torques on joints 2 and 3 increase as the angle θ p increases. Similarly, we find in Figure 7d that as angle θ p increases, d m becomes larger. In addition, it can be seen from Figure 7d that when θ p remains unchanged, d m decreases as the workpiece is placed closer to the robot. In particular, when vector p is equal to [0.6 m, 0.5 m, π/3] T and [0.6 m, −0.5 m, π/3] T , d m dramatically increases since the natural frequencies of the robot along the milling path with these workpiece poses are close to twice that of the milling force frequency (66.7 Hz).  The variation of dm in the xpyp plane when θp is equal to 0, π/12 rad, π/6 rad, and π/3 rad.

Optimization Results
In this section, the workpiece pose is optimized to minimize the mean value of the steady-state amplitudes of the cutting tool offset dm along the milling path defined in Figure 5b. Accordingly, the optimization problem in Equation (8) can be rewritten for this task as below: To capture the possible coupling between the three optimization variables, a global optimization algorithm is used for solving the problem. The algorithm utilized here, called optQuest/NLP or OQNLP [47], is designed to find the global optima for pure and mixed integer nonlinear problems with multiple constraints and variables, where all the objective functions are differentiable with respect to the continuous variables. The optimization results and the other selected five sets of dates are listed for comparison in Table 2. The curves in Figure 8 show the steady-state amplitudes of the tool offset in the selected positions along the milling path with different workpiece poses p. We can find that dm introduced in Section 3.3 is strongly influenced by the workpiece pose. In addition, when the workpiece pose remains unchanged, dm can vary greatly at different points along the milling path. The optimal workpiece pose is obtained as p = [0, −0.46, 0] T , which gives the smallest average tool offset amplitude dm of about 0.22 × 10 −3 m. Also, the tool offset changes relatively smoothly in different The variation of d m in the x p y p plane when θ p is equal to 0, π/12 rad, π/6 rad, and π/3 rad.

Optimization Results
In this section, the workpiece pose is optimized to minimize the mean value of the steady-state amplitudes of the cutting tool offset d m along the milling path defined in Figure 5b. Accordingly, the optimization problem in Equation (8) can be rewritten for this task as below: x p ∈ (0, 0.6) y p ∈ (−0.5, 0.5) To capture the possible coupling between the three optimization variables, a global optimization algorithm is used for solving the problem. The algorithm utilized here, called optQuest/NLP or OQNLP [47], is designed to find the global optima for pure and mixed integer nonlinear problems with multiple constraints and variables, where all the objective functions are differentiable with respect to the continuous variables. The optimization results and the other selected five sets of dates are listed for comparison in Table 2. The curves in Figure 8 show the steady-state amplitudes of the tool offset in the selected positions along the milling path with different workpiece poses p. We can find that d m introduced in Section 3.3 is strongly influenced by the workpiece pose. In addition, when the workpiece pose remains unchanged, d m can vary greatly at different points along the milling path. The optimal workpiece pose is obtained as p = [0, −0.46, 0] T , which gives the smallest average tool offset amplitude d m of about 0.22 × 10 −3 m. Also, the tool offset changes relatively smoothly in different positions along the path compared with the variations with other workpiece pose parameters. d mi represents the d m at the position of the ith workpiece position. positions along the path compared with the variations with other workpiece pose parameters. dmi represents the dm at the position of the ith workpiece position. Table 2. dmi with different places of workpiece. * represents the set of result as the optimal workpiece pose.

Conclusions
Structural vibration induced by periodic machining force is one of the main limitations for the application of industrial robots in machining. This paper presents a method for workpiece pose optimization for robotic milling to improve the quasi-static performance of the system. The method integrates a structural dynamics model of flexible-joint robots, a cutting force model for cylindrical helical end milling, and a global optimization for the workpiece pose. A numerical simulation is carried out for a robotic milling system for a given task. The results show that the quasi-static performance of the robot is strongly influenced by the workpiece pose and can change dramatically in different positions along the milling path. The optimal workpiece pose can effectively reduce the overall tool offset due to the periodic milling force and can lower the variation of the tool offset along the path. The results are helpful to improve the quasi-static performance including the accuracy and quality of robotic machining. The following future research directions are recommended. First, the regenerative effect on the undeformed chip thickness and the dynamic response of the robot should be considered. Second, the nonlinear behavior of robot joints should be further studied in the model.

Conclusions
Structural vibration induced by periodic machining force is one of the main limitations for the application of industrial robots in machining. This paper presents a method for workpiece pose optimization for robotic milling to improve the quasi-static performance of the system. The method integrates a structural dynamics model of flexible-joint robots, a cutting force model for cylindrical helical end milling, and a global optimization for the workpiece pose. A numerical simulation is carried out for a robotic milling system for a given task. The results show that the quasi-static performance of the robot is strongly influenced by the workpiece pose and can change dramatically in different positions along the milling path. The optimal workpiece pose can effectively reduce the overall tool offset due to the periodic milling force and can lower the variation of the tool offset along the path. The results are helpful to improve the quasi-static performance including the accuracy and quality of robotic machining. The following future research directions are recommended. First, the regenerative effect on the undeformed chip thickness and the dynamic response of the robot should be considered. Second, the nonlinear behavior of robot joints should be further studied in the model.

Appendix A
The tangential and radial cutting force acting on a microelement of the milling tool with a thickness of h(θ) can be expressed as dF t = K tc h(θ)dz dF r = K rc dF t , where K tc and K rc are the tangential and radial cutting force coefficients respectively and are experimentally determined as the following [43]: where t c is the average chip thickness value over one cutter revolution, and for a rigid tool, it can be given as The elemental tangential and radial forces can be further written in the tool coordinate system as below: dF x = −dF t cos θ − dF r sin θ dF y = dF t sin θ − dF r cos θ , where dF x and dF y are the elemental forces in the x t and y t directions acting on the microelement with an axial height of z', respectively. We can integrate Equation (12) about the axial height numerically to obtain the cutting force acting on a single cutting edge. Lastly, the instantaneous cutting force of the whole milling cutter at a certain time can be obtained by adding up the cutting forces on all the cutting edges and expressed as where N e represents the number of cutting edges and the upper and lower bounds z 1 and z 2 represent the height of the lower and upper ends of the cutting edge. As the cutter tool rotates, the length of the cutting edge involved in milling first increases and then decreases until zero. Accordingly, single-tooth milling can be divided into three stages: the entry phase, the stabilization stage, and the off stage. The upper and lower limits for the integration for these three cutting stages are listed in Table A1, where the width of the cutting layer a w can be expressed as Finally, the upper and lower limits of the next cutter tooth can be obtained by its correlation with the first cutter tooth. Table A1. The integral upper and lower limits of a cutter tooth.