# Reducing the Frame Vibration of Delta Robot in Pick and Place Application: An Acceleration Profile Optimization Approach

Authors: Hongtai Cheng; Wei Li
Journal: Shock and Vibration
 (2018)

---

## Abstract

Delta robot is typically mounted on a frame and performs high speed pick and place tasks from top to bottom. Because of its outstanding accelerating capability and higher center of mass, the Delta robot can generate significant frame vibration. Existing trajectory smoothing methods mainly focus on vibration reduction for the robot instead of the frame, and modifying the frame structure increases the manufacturing cost. In this paper, an acceleration profile optimization approach is proposed to reduce the Delta robot-frame vibration. The profile is determined by the maximum jerk, acceleration, and velocity. The pick and place motion (PPM) and resulting frame vibration are analyzed in frequency domain. Quantitative analysis shows that frame vibration can be reduced by altering those dynamic motion parameters. Because the analytic model is derived based on several simplifications, it cannot be directly applied. A surrogate model-based optimization method is proposed to solve the practical issues. By directly executing the PPM with different parameters and measuring the vibration, a model is derived using Gaussian Process Regression (GPR). In order to reduce the frame vibration without sacrificing robot efficiency, those two goals are fused together according to their priorities. Based on the surrogate model, a single objective optimization problem is formulated and solved by Genetic Algorithm (GA). Experimental results show effectiveness of the proposed method. Behavior of the optimal parameters also verifies the robot-frame vibration mechanism.


---

## Body

## 1. Introduction

Delta robot is a 3/4 Degree-of-Freedoms(DOFs) parallel robot widely used in manufacturing industries affording the pick and place tasks[1–3]. A typical sorting work cell is shown in Figure 1. Delta robot is usually mounted on a big heavy frame to operate workpieces from top to bottom. The frame is usually customized according to the environmental conditions, such as workpiece size or conveyor arrangement. Typically, a frame is 2 meters high and 1.5 meters wide, which covers work range of the Delta robot. Thick square steels are welded together to provide a relatively static basis. Other types of base/frames also exist which are summarized in [4].Figure 1
The Delta robot based pick and place work cell.Although lightweight in mass and small in inertia, the Delta robot can still generate huge reaction force because of its excellent acceleration capability. Moreover, Center of Mass (CoM) of the frame is higher than normal robot base. Therefore, stiffness of the frame is an important factor affecting system noise level, reliability, and precision. Pick and Place Motion (PPM) is the primary behavior of Delta robot. As shown in Figure1, it moves between two locations along a door shaped trajectory, picking and placing the workpieces. When the robot moves fast, there will be huge reaction force and the frame is forced to vibrate. When the robot moves cyclically, the vibrations may superimpose together and cause significant side effect. Severe frame vibration also affects stability of the frame mounted industrial cameras. Therefore, reducing the frame vibration is a nonnegligible problem when installing new Delta robots.Smoothing the trajectory and strengthening the frame are two feasible solutions. A lot of trajectory smoothing methods have been proposed to reduce vibration of the Delta robot, such as using special functions to smooth the piecewise linear motions and optimizing the dynamic parameters such as jerk or acceleration.Dai and Sheng propose to use fifth Bézier curve to blend the adjacent positions and linear motions, which can generate a smooth trajectory with continuous second-order derivative and curvature in real-time [5]. Huang adopts the elliptical trajectory with modified sine motion profile to plan the pick and place operation, which can guarantee smoothness of the torque. Because the motion planning algorithm only involves two parameters, it can be efficiently implemented. However, no quantitative analysis is performed for the robot vibrations [6].Kuo studies the mathematical modeling problem of the Delta robot with flexible links. Kineto-elastodynamics and the finite element method are used to derive a mathematical model. Natural frequency analysis is performed to demonstrate the importance of flexible link compensation for precise robot motion [7]. However, the frequency analysis is not applied to motion optimization problem and the method is too complex for practical application. In order to improve the tracking accuracy, Liu studies the trajectory planning problem for Delta robot with joint friction and jerk constraints [8, 9]. Smooth and nearly time optimal trajectory can be generated using the proposed solution. The trajectory is modeled using quadratic and quantic polynomial splines in Cartesian space and septuple polynomial splines in joint space. Zhang and Wang study the trajectory planning and optimization method to obtain smoother and more efficient trajectories for Delta robot[10, 11]. Béarée introduces input shaping technique into the jerk-limited profile, derives the damped-jerk profile to cancel the residual vibration of an undamped system, and proposes several practical tuning rules for the equivalent input filter parameters [12]. The dominating flexible mode analysis closely relates to system physical parameters; therefore, it is difficult to apply in the vibration suspension problem of Delta robot and frame.The causes of frame and robot vibration are different. Robot vibration originates from discontinuity of the drive torque [13, 14], nonlinearity of the robot structure[15–17], or flexibility of the links [7]. Frame vibration is excited by the force coupled through the static base. It originates from motion of the robot. No matter how smooth the robot motion is, there is always excitation force imposing on the frame. The only ways are modifying either the frame structure or the robot motion.A commonly used solution is to modify the frame structure. Increasing the frame mass and stiffness or integrating reinforcement beam can all alter frame’s natural frequency characteristics. The designing process involves 3D modeling, Finite Element Analysis (FEA) on the stiffness, deformation, and frequency characteristics of the new structure [4]. Apparently, it is not an economic solution because the frame should be specifically designed for each application. More materials and costs are required. Even though, for high speed pick and place motions, the vibration is still significant.The other method is to modify the robot motion. However, it is restricted by task and physical constraints. A feasible way is to customize the dynamic motion parameters for each application. This can also change frequency characteristics of the excitation signal imposed on the robot frame. Consequently, the frame vibration can be reduced without modifying the hardware design, software architecture, and sacrificing the efficiency. Furthermore, because the proposed method does not change the underlying motion smoothing algorithms, there is no side effect on the robot vibration.Vibration signal contains rich information about the state of the robot/machine and hence can be used to assist the high-level diagnosis and decision-making process. A lot of model and signal based methods have been proposed to solve the fault detection problem for machines [18–20] and robots [21, 22]. In this paper, the vibration signals are used to find the optimal dynamic motion parameters.There are two objectives for the optimization problem: reducing the frame vibration and keeping efficiency. It is noticed that they have different priorities. Efficiency is the primary goal. Therefore, in this paper, a specific fitness function is designed to fuse them together safely. Although a lot of multiobjective optimization algorithms have been proposed [23, 24], single objective optimization problems are easier to handle because there is no trade-off between each objective.The paper has the following contributions.(1) A frequency domain analysis is performed on the PPM acceleration profile. The internal relation between dynamic motion parameters and distribution/ amplitude of the harmonic components is established, which proves the possibility of reducing the frame vibration by optimizing the acceleration profile. (2) A surrogate model-based parameter optimization method is proposed to solve the practical vibration reduction problem. The method uses a nonparametric regression tool: Gaussian Process Regression (GPR) to generate the surrogate model which makes it easier to be applied in other acceleration profiles. The method fuses two objectives into one single fitness function based on their priority differences. It can be easily integrated in existing optimization algorithms.


## 2. Problem Formulation

### 2.1. Delta Robot in Pick and Place Application

Delta robots are typically used in pick and place applications. Figure1 shows a commonly used configuration. The workpieces are supplied through conveyors and the robot picks them from conveyor and places them on target locations (boxes, plates, or other conveyors). Typically, in order to facilitate tracking of the conveyor, moving direction of the conveyor is aligned with X-axis of the Delta robot. Therefore, pick and place motion (PPM) is mainly within YZ plane as shown in Figure 1. In the paper, the conveyor is assumed to align with X-axis for convenience. Results of aligning with Y-axis are similar.Pick and Place Motion (PPM) is a door shaped trajectory consisted of one horizontal, two vertical, and two smooth transition segments as shown in Figure2 [25]. The trajectory formed by points P0→P3 in Figure 2 is a typical PPM, where P0 is the pick-up location and P3 is the put-down location. P1 and P2 define the intermediate transition positions. The robot picks workpiece up at P0 and moves it to P3 along the predefined trajectory. During this process, the robot accelerates and decelerates cyclically and swiftly. Huge reaction forces are imposed on the robot and the frame. Therefore, smoothness of the PPM greatly affects the vibrational amplitude.Figure 2
The door shape pick and place motion.When the robot moves along PPM cyclically and swiftly, the reaction force will make the frame vibrate. This will generate noise and affect the robot accuracy. The problem is to reduce the frame vibration while keeping system efficiency based on existing hardware and software configuration.


### 2.2. Acceleration Profile and Frame Vibration

Efficiency and frequency of the robot motion are controlled by the acceleration profile, which refers to how the velocity and acceleration change in the trajectory tracking control process. The jerk-limited profile is a smooth command pattern used by modern motion systems such as machine-tools and industrial robots. Jerk is time derivative of acceleration. The introduced jerk can eliminate the discontinuities of acceleration and reduce the residual vibration.Figure3 shows an example plot of jerk-limited profile. By using the bang-bang control law, a time optimal accelerating trajectory is derived. A jerk-limited profile equals to a trapezoidal or a triangular acceleration profile.Figure 3
An example plot of jerk-limited acceleration profile.A point to point motion profile can be expressed using the following equation:(1)a¯t∣S,Jmax,Amax,Vmax=Jmaxt-t0t0≤t<t1,t1=t0+Δ1Amaxt1≤t<t2,t2=t1+Δ2-Δ1Amax-Jmaxt-t2t2≤t<t3,t3=t2+Δ10t3≤t<t4,t4=t2+Δ-Jmaxt-t4t4≤t<t5,t5=t4+Δ1-Amaxt5≤t<t6,t6=t5+Δ2-Δ1-Amax+Jmaxt-t6t6≤t<t7,t7=t6+Δ1where S is the stroke, Jmax,Amax,Vmax are the dynamic motion parameters for acceleration profile, Δ is the time when the robot is in uniform velocity state, which is proportional to S, and Δ1,Δ2 are the accelerating time.(2)Δ1=AmaxJmax(3)Δ2=VmaxAmaxA detailed discussion on the jerk-limited acceleration profile can be found in [12, 26]. Besides the symmetric profile, there are also many variations such as the asymmetric acceleration profile in [27]. No matter which type the robot uses, the profile all can be parameterized. For the standard jerk-limited acceleration profile, the dynamic parameters Jmax,Amax,Vmax greatly affect the robot motion efficiency and smoothness. Typically they are chosen empirically to achieve small cycle time or to reduce the robot vibration.Significant frame vibration originates from mismatch of the robot and frame frequency characteristics. Before installation, the frame can be reinforced to match the application. However, this increases the material and labor costs. The acceleration profile is the key to change the robot motion frequency characteristics. Therefore, an economic and general applicable solution is to customize the dynamic parameters according to specific application based on the existing frame configuration. Instead of empirically setting those parameters, it is necessary to find an optimal solution for efficient and lower vibration motion. The following sections further discuss this problem from mechanism to solutions.


## 2.1. Delta Robot in Pick and Place Application

Delta robots are typically used in pick and place applications. Figure1 shows a commonly used configuration. The workpieces are supplied through conveyors and the robot picks them from conveyor and places them on target locations (boxes, plates, or other conveyors). Typically, in order to facilitate tracking of the conveyor, moving direction of the conveyor is aligned with X-axis of the Delta robot. Therefore, pick and place motion (PPM) is mainly within YZ plane as shown in Figure 1. In the paper, the conveyor is assumed to align with X-axis for convenience. Results of aligning with Y-axis are similar.Pick and Place Motion (PPM) is a door shaped trajectory consisted of one horizontal, two vertical, and two smooth transition segments as shown in Figure2 [25]. The trajectory formed by points P0→P3 in Figure 2 is a typical PPM, where P0 is the pick-up location and P3 is the put-down location. P1 and P2 define the intermediate transition positions. The robot picks workpiece up at P0 and moves it to P3 along the predefined trajectory. During this process, the robot accelerates and decelerates cyclically and swiftly. Huge reaction forces are imposed on the robot and the frame. Therefore, smoothness of the PPM greatly affects the vibrational amplitude.Figure 2
The door shape pick and place motion.When the robot moves along PPM cyclically and swiftly, the reaction force will make the frame vibrate. This will generate noise and affect the robot accuracy. The problem is to reduce the frame vibration while keeping system efficiency based on existing hardware and software configuration.


## 2.2. Acceleration Profile and Frame Vibration

Efficiency and frequency of the robot motion are controlled by the acceleration profile, which refers to how the velocity and acceleration change in the trajectory tracking control process. The jerk-limited profile is a smooth command pattern used by modern motion systems such as machine-tools and industrial robots. Jerk is time derivative of acceleration. The introduced jerk can eliminate the discontinuities of acceleration and reduce the residual vibration.Figure3 shows an example plot of jerk-limited profile. By using the bang-bang control law, a time optimal accelerating trajectory is derived. A jerk-limited profile equals to a trapezoidal or a triangular acceleration profile.Figure 3
An example plot of jerk-limited acceleration profile.A point to point motion profile can be expressed using the following equation:(1)a¯t∣S,Jmax,Amax,Vmax=Jmaxt-t0t0≤t<t1,t1=t0+Δ1Amaxt1≤t<t2,t2=t1+Δ2-Δ1Amax-Jmaxt-t2t2≤t<t3,t3=t2+Δ10t3≤t<t4,t4=t2+Δ-Jmaxt-t4t4≤t<t5,t5=t4+Δ1-Amaxt5≤t<t6,t6=t5+Δ2-Δ1-Amax+Jmaxt-t6t6≤t<t7,t7=t6+Δ1where S is the stroke, Jmax,Amax,Vmax are the dynamic motion parameters for acceleration profile, Δ is the time when the robot is in uniform velocity state, which is proportional to S, and Δ1,Δ2 are the accelerating time.(2)Δ1=AmaxJmax(3)Δ2=VmaxAmaxA detailed discussion on the jerk-limited acceleration profile can be found in [12, 26]. Besides the symmetric profile, there are also many variations such as the asymmetric acceleration profile in [27]. No matter which type the robot uses, the profile all can be parameterized. For the standard jerk-limited acceleration profile, the dynamic parameters Jmax,Amax,Vmax greatly affect the robot motion efficiency and smoothness. Typically they are chosen empirically to achieve small cycle time or to reduce the robot vibration.Significant frame vibration originates from mismatch of the robot and frame frequency characteristics. Before installation, the frame can be reinforced to match the application. However, this increases the material and labor costs. The acceleration profile is the key to change the robot motion frequency characteristics. Therefore, an economic and general applicable solution is to customize the dynamic parameters according to specific application based on the existing frame configuration. Instead of empirically setting those parameters, it is necessary to find an optimal solution for efficient and lower vibration motion. The following sections further discuss this problem from mechanism to solutions.


## 3. Frequency Domain Analysis of PPM and Frame Vibration

When Delta robot performs the pick and place tasks, acceleration and deceleration occur cyclically. Consequently, a periodic reaction force drives the frame to vibrate. The robot motion can be seen as the input while the frame vibration is the output. The robot-frame union is a black box with certain frequency response characteristics. To understand the vibration excitation mechanism and find methods to reduce the vibration, it is necessary to analyze the robot motion and robot-frame union in frequency domain.


### 3.1. Spectrum of the Pick and Place Motion

As shown in Figure1, in pick and place applications, Delta robots usually move between two locations cyclically. In the door shaped trajectory, horizontal segment is the dominant one which covers most of the stroke. Therefore, it can be treated as a cyclic point to point motion for simplicity (P1→P2, P2→P1). When applying trapezoidal acceleration profile, the corresponding motion is illustrated in Figure 4(a), where the thin solid line is the acceleration.Figure 4
An illustrative example of the trapezoidal acceleration profile and the corresponding spectrum for cyclic motion.
(a)
Trapezoidal acceleration profile (b)
The corresponding spectrumConsidering the point to point motion (1) and symmetricity, acceleration profile of the cyclic motion can be modeled using(4)at=a¯t∣S,Jmax,Amax,VmaxTi-1≤t<Ti2,P0⟶P3a¯t∣-S,Jmax,Amax,VmaxTi2≤t<Ti+12,P3⟶P0,i=1,2,3,⋯which is a piecewise linear function as shown in Figure 4(a), where a¯ is an elementary point to point motion whose period is T; S is the stroke between initial and target positions; Vmax,Amax,Jmax are the maximum velocity, acceleration, and jerk.Because the movement between P0 and P3 are symmetric, they both cover half of the period. Therefore, one can get(5)at=-at+T2=-at+12fb=-at+πwbwhere fb is the base frequency and wb is the corresponding angular frequency. From Figure 4(a) and (5), it is able to find that the acceleration is actually an odd harmonic function. Considering the Fourier transformation with w=kwb, there is(6)Ukwb=Fat=∫-∞+∞ate-ikwbtdt=-∫-∞+∞at+πwbe-ikwbtdt=-∫-∞+∞at+πwbe-ikwbt+π/wbe-ikπdt+πwb=-e-ikπ∫-∞+∞aτe-ikwbτdτ=-e-ikπUkwbIf k is odd, -e-ikπ=1, which will lead to an identity. If k is even, -e-ikπ=-1, which will lead to (7)Ukwb=-Ukwb.There is U(kwb)=0. Therefore, for the cyclic PPM motion, only odd harmonics such as fb,3fb,5fb,⋯ exist.Figure4(b) shows Fast Fourier Transform (FFT) results for the acceleration profile shown in Figure 4(a). It is clearly seen that the base frequency is 0.44Hz and the harmonics are 1.32Hz, 2.2Hz, which is consistent with the theoretical analysis.


### 3.2. Amplitude-Frequency Characteristics Analysis

According to Figures 4(a), (4), and (1), it is able to calculate the amplitude-frequency characteristics of the acceleration profile. Noticing that the jerk profile is simper (piecewise constant) than the acceleration profile (piecewise linear), the Fourier transformation can be performed on jerk based on the facts that(8)Ff˙t=iwFft(9)a˙t=jtwhere j(t) refers to the jerk profile. The derivation process of the differential property (8) for Fourier transformation can be found in [28]. The Fourier expansion coefficient for acceleration profile is(10)Ukwb=1T∫-T/2T/2ate-ikwbtdt=ikwb1T∫-T/2T/2jte-ikwbtdtwhere(11)∫-T/2T/2jte-ikwbtdt=∫-T/2-T/2+Δ1-Jmaxe-ikwbtdt+∫-T/2+Δ2-T/2+Δ1+Δ2Jmaxe-ikwbtdt+∫-Δ1-Δ2-Δ2Jmaxe-ikwbtdt+∫-Δ10-Jmaxe-ikwbtdt+∫0Δ1Jmaxe-ikwbtdt+∫Δ2Δ1+Δ2-Jmaxe-ikwbtdt+∫T/2-Δ1-Δ2T/2-Δ2-Jmaxe-ikwbtdt+∫T/2-Δ1T/2Jmaxe-ikwbtdt=-Jmaxikwb-e-ikwbt-T/2-T/2+Δ1+e-ikwbt-T/2+Δ2-T/2+Δ1+Δ2+e-ikwbt-Δ1-Δ2-Δ2-e-ikwbt-Δ10+e-ikwbt0Δ1-e-ikwbtΔ2Δ1+Δ2-e-ikwbtT/2-Δ1-Δ2T/2-Δ2+e-ikwbtT/2-Δ1T/2Assuming k is odd, using Euler formula to expand above equation and merge similar items, one gets(12)∫-T/2T/2jte-ikwbtdt=-4Jmaxikwbcos⁡2πkΔ1T+cos⁡2πkΔ1+Δ2T-cos⁡2πkΔ2T-1Substituting (12) into (10), one gets(13)Uk=ikwb1T∫-T/2T/2jte-ikwbtdt=ikwb1T-4Jmaxikwbcos⁡2πkΔ1T+cos⁡2πkΔ1+Δ2T-cos⁡2πkΔ2T-1=-4JmaxTcos⁡2πkΔ1T+cos⁡2πkΔ1+Δ2T-cos⁡2πkΔ2T-1According to (13), (2), and (3), one can find that the dynamic motion parameters Vmax,Amax,Jmax are the only factors affecting the acceleration profile, cycle time, and amplitude-frequency characteristics. Larger Jmax and smaller cycle time T tend to increase the amplitude. However, items in the brackets can also be reduced by properly choosing those dynamic motion parameters. This indicates that it is possible to reduce the frame vibration by tuning the acceleration profile instead of modifying the hardware, and a feasible way is to optimize dynamic motion parameters according to given constraints.An illustrative example is shown in Figure5. The cycle time T is set to be 0.5s, i.e., locations of the harmonic components are fixed. Different acceleration profiles are generated by changing Amax=[10000,50000]mm/s2 and Jmax=[1000000,2000000]mm/s3 as shown in Figure 5(a).Figure 5
An example illustrating the relation between the dynamic motion parameters and amplitude-frequency characteristics.
(a)
Variations of the acceleration profile (b)
Amplitude of the corresponding harmonicsAmplitudes of base frequency, 3rd and 5th harmonics are plotted in Figure5(b). It is seen that amplitudes vary according to the dynamic motion parameters. Globally, the amplitude is proportional to Jmax. However, local minimums also exist which can generate lower vibration without lowering production efficiency. This is consistent with (13).


### 3.3. Measures to Reduce Frame Vibration

Different from the robots, the robot frames usually have no fixed models. They are customized for each application to meet the specific conveyor width and height, workspace arrangement, and workpiece size requirements. Although weighs hundreds kilograms, it can still vibrate with some natural frequency. This frequency is decided by the frame material, structure, dimension, and mounting type. When external excitation signal is close to the natural frequency, the frame generates significant vibration.Frequency response of the frame can be described using(14)Yf=GfUf=GfUf∣Vmax,Amax,Jmaxwhere G(f) is the amplitude-frequency characteristics of robot-frame union whose natural frequency is fn. U(f) and Y(f) are the robot motion and frame vibration signals in frequency domain. For a given robot-frame union, G(f) is determined. Without changing the hardware, the only way is to tuning the robot motion U(f).The robot motion can be tuned in two aspects.(i)
Harmonics distributionFor the PPM motion, only odd harmonics exist in the acceleration profile. The proximity between2k-1fb and fn is defined by(15)df=mink=1,2,⋯⁡2k-1×fb-fnIt is easy to verify that0≤df≤fb. df=0 will lead to resonance. Bigger df can decrease the vibration. The harmonics distribution is closely related to cycle time T because fb=1/T. One can alter T to tune the distribution to avoid the natural frequency.(ii)
Natural Frequency AmplitudeThe other method is to decrease the amplitude ofU(fn). According to (13), this can be realized by increasing T, decreasing Jmax, and optimizing the dynamic motion parameters. However, it is hardly possible to directly find such parameters because U(fn) and T are internally related in a complex reverse way.Remark 1.
As shown in Figure1, typically X-axis of Delta robot is aligned with the conveyor and the pick and place motion mainly lies in the Y-Z plane. Because the vertical motion is usually very small compared to the horizontal motion, the PPM motion is mainly in Y direction. Therefore, the acceleration profile and the vibration model are both considered as one dimensional in this paper.
The above theoretical analysis reveals the mechanism for frame vibration and shows two possible methods to reducing the frame vibration by tuning the acceleration profile of PPM. However, it is not suitable for practical application because of the following four reasons.(1) Real PPM motion has vertical segments and smooth transition segments. (2) Formula of acceleration profile may change for certain parameters. (3) Real PPM may derivate from Y-axis. (4) There is a lack of robot-frame union model.


## 3.1. Spectrum of the Pick and Place Motion

As shown in Figure1, in pick and place applications, Delta robots usually move between two locations cyclically. In the door shaped trajectory, horizontal segment is the dominant one which covers most of the stroke. Therefore, it can be treated as a cyclic point to point motion for simplicity (P1→P2, P2→P1). When applying trapezoidal acceleration profile, the corresponding motion is illustrated in Figure 4(a), where the thin solid line is the acceleration.Figure 4
An illustrative example of the trapezoidal acceleration profile and the corresponding spectrum for cyclic motion.
(a)
Trapezoidal acceleration profile (b)
The corresponding spectrumConsidering the point to point motion (1) and symmetricity, acceleration profile of the cyclic motion can be modeled using(4)at=a¯t∣S,Jmax,Amax,VmaxTi-1≤t<Ti2,P0⟶P3a¯t∣-S,Jmax,Amax,VmaxTi2≤t<Ti+12,P3⟶P0,i=1,2,3,⋯which is a piecewise linear function as shown in Figure 4(a), where a¯ is an elementary point to point motion whose period is T; S is the stroke between initial and target positions; Vmax,Amax,Jmax are the maximum velocity, acceleration, and jerk.Because the movement between P0 and P3 are symmetric, they both cover half of the period. Therefore, one can get(5)at=-at+T2=-at+12fb=-at+πwbwhere fb is the base frequency and wb is the corresponding angular frequency. From Figure 4(a) and (5), it is able to find that the acceleration is actually an odd harmonic function. Considering the Fourier transformation with w=kwb, there is(6)Ukwb=Fat=∫-∞+∞ate-ikwbtdt=-∫-∞+∞at+πwbe-ikwbtdt=-∫-∞+∞at+πwbe-ikwbt+π/wbe-ikπdt+πwb=-e-ikπ∫-∞+∞aτe-ikwbτdτ=-e-ikπUkwbIf k is odd, -e-ikπ=1, which will lead to an identity. If k is even, -e-ikπ=-1, which will lead to (7)Ukwb=-Ukwb.There is U(kwb)=0. Therefore, for the cyclic PPM motion, only odd harmonics such as fb,3fb,5fb,⋯ exist.Figure4(b) shows Fast Fourier Transform (FFT) results for the acceleration profile shown in Figure 4(a). It is clearly seen that the base frequency is 0.44Hz and the harmonics are 1.32Hz, 2.2Hz, which is consistent with the theoretical analysis.


## 3.2. Amplitude-Frequency Characteristics Analysis

According to Figures 4(a), (4), and (1), it is able to calculate the amplitude-frequency characteristics of the acceleration profile. Noticing that the jerk profile is simper (piecewise constant) than the acceleration profile (piecewise linear), the Fourier transformation can be performed on jerk based on the facts that(8)Ff˙t=iwFft(9)a˙t=jtwhere j(t) refers to the jerk profile. The derivation process of the differential property (8) for Fourier transformation can be found in [28]. The Fourier expansion coefficient for acceleration profile is(10)Ukwb=1T∫-T/2T/2ate-ikwbtdt=ikwb1T∫-T/2T/2jte-ikwbtdtwhere(11)∫-T/2T/2jte-ikwbtdt=∫-T/2-T/2+Δ1-Jmaxe-ikwbtdt+∫-T/2+Δ2-T/2+Δ1+Δ2Jmaxe-ikwbtdt+∫-Δ1-Δ2-Δ2Jmaxe-ikwbtdt+∫-Δ10-Jmaxe-ikwbtdt+∫0Δ1Jmaxe-ikwbtdt+∫Δ2Δ1+Δ2-Jmaxe-ikwbtdt+∫T/2-Δ1-Δ2T/2-Δ2-Jmaxe-ikwbtdt+∫T/2-Δ1T/2Jmaxe-ikwbtdt=-Jmaxikwb-e-ikwbt-T/2-T/2+Δ1+e-ikwbt-T/2+Δ2-T/2+Δ1+Δ2+e-ikwbt-Δ1-Δ2-Δ2-e-ikwbt-Δ10+e-ikwbt0Δ1-e-ikwbtΔ2Δ1+Δ2-e-ikwbtT/2-Δ1-Δ2T/2-Δ2+e-ikwbtT/2-Δ1T/2Assuming k is odd, using Euler formula to expand above equation and merge similar items, one gets(12)∫-T/2T/2jte-ikwbtdt=-4Jmaxikwbcos⁡2πkΔ1T+cos⁡2πkΔ1+Δ2T-cos⁡2πkΔ2T-1Substituting (12) into (10), one gets(13)Uk=ikwb1T∫-T/2T/2jte-ikwbtdt=ikwb1T-4Jmaxikwbcos⁡2πkΔ1T+cos⁡2πkΔ1+Δ2T-cos⁡2πkΔ2T-1=-4JmaxTcos⁡2πkΔ1T+cos⁡2πkΔ1+Δ2T-cos⁡2πkΔ2T-1According to (13), (2), and (3), one can find that the dynamic motion parameters Vmax,Amax,Jmax are the only factors affecting the acceleration profile, cycle time, and amplitude-frequency characteristics. Larger Jmax and smaller cycle time T tend to increase the amplitude. However, items in the brackets can also be reduced by properly choosing those dynamic motion parameters. This indicates that it is possible to reduce the frame vibration by tuning the acceleration profile instead of modifying the hardware, and a feasible way is to optimize dynamic motion parameters according to given constraints.An illustrative example is shown in Figure5. The cycle time T is set to be 0.5s, i.e., locations of the harmonic components are fixed. Different acceleration profiles are generated by changing Amax=[10000,50000]mm/s2 and Jmax=[1000000,2000000]mm/s3 as shown in Figure 5(a).Figure 5
An example illustrating the relation between the dynamic motion parameters and amplitude-frequency characteristics.
(a)
Variations of the acceleration profile (b)
Amplitude of the corresponding harmonicsAmplitudes of base frequency, 3rd and 5th harmonics are plotted in Figure5(b). It is seen that amplitudes vary according to the dynamic motion parameters. Globally, the amplitude is proportional to Jmax. However, local minimums also exist which can generate lower vibration without lowering production efficiency. This is consistent with (13).


## 3.3. Measures to Reduce Frame Vibration

Different from the robots, the robot frames usually have no fixed models. They are customized for each application to meet the specific conveyor width and height, workspace arrangement, and workpiece size requirements. Although weighs hundreds kilograms, it can still vibrate with some natural frequency. This frequency is decided by the frame material, structure, dimension, and mounting type. When external excitation signal is close to the natural frequency, the frame generates significant vibration.Frequency response of the frame can be described using(14)Yf=GfUf=GfUf∣Vmax,Amax,Jmaxwhere G(f) is the amplitude-frequency characteristics of robot-frame union whose natural frequency is fn. U(f) and Y(f) are the robot motion and frame vibration signals in frequency domain. For a given robot-frame union, G(f) is determined. Without changing the hardware, the only way is to tuning the robot motion U(f).The robot motion can be tuned in two aspects.(i)
Harmonics distributionFor the PPM motion, only odd harmonics exist in the acceleration profile. The proximity between2k-1fb and fn is defined by(15)df=mink=1,2,⋯⁡2k-1×fb-fnIt is easy to verify that0≤df≤fb. df=0 will lead to resonance. Bigger df can decrease the vibration. The harmonics distribution is closely related to cycle time T because fb=1/T. One can alter T to tune the distribution to avoid the natural frequency.(ii)
Natural Frequency AmplitudeThe other method is to decrease the amplitude ofU(fn). According to (13), this can be realized by increasing T, decreasing Jmax, and optimizing the dynamic motion parameters. However, it is hardly possible to directly find such parameters because U(fn) and T are internally related in a complex reverse way.Remark 1.
As shown in Figure1, typically X-axis of Delta robot is aligned with the conveyor and the pick and place motion mainly lies in the Y-Z plane. Because the vertical motion is usually very small compared to the horizontal motion, the PPM motion is mainly in Y direction. Therefore, the acceleration profile and the vibration model are both considered as one dimensional in this paper.
The above theoretical analysis reveals the mechanism for frame vibration and shows two possible methods to reducing the frame vibration by tuning the acceleration profile of PPM. However, it is not suitable for practical application because of the following four reasons.(1) Real PPM motion has vertical segments and smooth transition segments. (2) Formula of acceleration profile may change for certain parameters. (3) Real PPM may derivate from Y-axis. (4) There is a lack of robot-frame union model.


## 4. Modeling and Parameter Optimization

Because of the complex robot-frame structure, it is hardly possible to derive an accurate analytic model for the excitation and vibration process. In previous section, robot motion and frame vibration are analyzed in frequency domain. Guidelines for reducing the frame vibration are identified. However, there still lacks quantitative methods to optimize such dynamic motion parameters.This paper proposes a surrogate model-based parameter optimization solution. Instead of using the analytic model, a surrogate model is obtained by performing regression on real datasets. The modeling and optimization process is shown in Figure6.Figure 6
Diagram of the signal processing, modeling, and optimization process.The whole process consisted of three parts: data acquisition, signal processing and modeling, and dynamic motion parameter optimization.


### 4.1. Data Acquisition

To obtain a surrogate model for (14), all related variables should be able to be extracted from the recorded datasets. According to (14), robot motion and frame vibration are both relevant. Therefore, acceleration profile of robot motion and frame vibration is recorded simultaneously. The recorded signals are time series a(k)={ar(k),af(k)}, where ar(k) is the robot motion while af(k) is the frame vibration.The arguments arep=Vmax,Amax,Jmax, which are known in advance when performing the evaluation experiments. Therefore, the recorded dataset is(16)Dataset=p,ark,afkIn order to minimize the effect of random factors, the robot repeats the pick and place tasks several times for each parameter candidatep.


### 4.2. Signal Processing and Feature Extraction

Performances of the dynamic motion parameters are evaluated from two aspects: efficiency of the robot motion and amplitude of the frame vibration. In this paper, efficiency is measured by cycle timeT. Vibration is measured using maximum acceleration a^. Both indices can be calculated from the datasets.


#### 4.2.1. Cycle Time

The measurement of cycle timeT is based on the recorded time series ar(k). By evaluating the acceleration value, it is able to identify whether the robot is moving or staying still. According to the prior knowledge, the robot stays still before new parameters are received. The start time ts can be identified for each movement.Same operation can be used to find the stop timete. Cycle time can be calculated by(17)T=te-ts


#### 4.2.2. Maximum Acceleration

Vibration amplitude is evaluated using maximum acceleration, which is calculated using(18)a^=max⁡afk.The surrogate numerical model is obtained by performing regression over the datasets. The dynamic motion parameters areJmax, Amax, and Vmax and the performance indices are a^,T. However, the unit and value range of those variables cause serious problems. Typically, their values have significant differences, ranging from 0-2000000. Without normalization, some of the variables may be submerged by others and the modeling process may fail. To treat each variable equally, a normalization operation is required. Take T as example; the normalizing process follows(19)T~=T-μTσTwhere T~ refers to the normalized T; μ(T) and σ(T) stand for the mean and standard division of T. This operation can make all the variables zero mean and with identical variance. The original variable can be recovered easily through(20)T=T~×σT+μTSame operations can be performed on a^ and Jmax, Amax, and Vmax. After normalization, inputs and outputs of the model are converted to J~max,A~max,V~max, and T~,a~. The latent models can be represented using(21)T~=fTJ~max,A~max,V~maxa~=faJ~max,A~max,V~maxwhere fa,fT are the underlying functions.


### 4.3. Gaussian Process Regression Based Modeling

Considering the uncertainties in datasets, numbers of samples, and nonlinearity of the underlying function, Gaussian Process Regression (GPR) is used to modeling the datasets. GPR is a nonparametric tool that can handle the modeling problem with noisy observations and system uncertainties [29]. It has been applied in robot process modeling and parameter optimization [30].For the modeling problem (21), there is(22)x=J~max,A~max,V~maxTy=a~=fax,y=T~=fTxThe goal is to derive a model to approximate the underlying functionfax and fTx. Because there are many random factors in the robot system and the observations are also noisy, they can be seen as two Gaussian Processes and can be modeled using GPR.For a Gaussian Processf(x), a set of multivariate Gaussian random variables F={f(x1),f(x2),, ⋯,f(xN)} can be defined over X, where xi∈Xi=1N are the inputs and f(xi) are the corresponding outputs. X is the parameter space which is defined over RD. Instead of deriving an explicit form for f(x), it is approximated by its covariance function k(x,x′): f(x)~GP(0,k(x,x′,θ) where x and x′ are two arbitrary variables in X.The goal of model construction is to find a covariance functionk(x,x′,θ) which fits the data set (X,y) best. Suppose f(x) is a candidate latent function (f in short), the marginal likelihood of H and θ given the data set (X,y) can be derived by performing marginalization over the function f(23)py∣X,H,θ=∫py∣X,f,H,θpf∣X,H,θdfwhere H is the hypothesis on the structure of the covariance function, θ are the hyperparameters, X,y are the sample data sets, and(24)pf∣X,y,H,θ=py∣X,f,H,θpf∣X,H,θpy∣X,H,θis the posterior probability.In order to simplify the modeling process, a log operation is performed to convert multiplication into addition. The Gaussian assumption makes it possible to derive an analytic solution for the marginal log likelihood(25)log⁡py∣X,H,θ=-12yTK+σn2I-1y-12log⁡K+σn2I-n2log⁡2πThen the hyperparameters θ can be found by maximizing the marginal log likelihood function.(26)θ∗=argmaxθ⁡log⁡pyX,H,θThis optimization problem is solved by Conjugate Gradient Algorithm. After training, the model is constructed and can be used for prediction.The joint distribution of the modeling dataset(X,y) and test dataset (X∗,y∗) is(27)yy∗~N0,KX,X+σn2IKX,X∗KX∗,XKX∗,X∗where KX′,X′′ (X′ and X′′ refer to X and X∗) is the covariance matrix whose element Kij at ith row and jth column equals to k(xi′,xj′′). By deriving the conditional distribution, the predictive function is obtained(28)μy∗=Ey∗∣X,y,X∗=KX∗,XKX,X+σn2I-1yVy∗=KX∗,X∗-KX∗,XKX,X+σn2I-1KX,X∗where μ(y∗) is the predicted mean for y∗ and Vy∗ is the variance.In summery, a covariance function should be specified for the GPR model. In the proposed data acquisition procedure, the pick and place motion are repeated for several times for each parameter combination. The random noise in the extracted feature is averaged. Therefore the following covariance function is used:(29)kx,x′=∑i=12σi2exp⁡-∑j=13xj-xj′22lijwhere two Squared Exponential covariance functions with Automatic Relevance Determination (SEARD) items are used to provide feasibility in different length scales lij. Considering the amplitude σi, there are 8 hyperparameters. Based on (29) and the obtained dataset (22), the optimal hyperparameters are found by optimizing (26). Finally two models are created.(30)T~=GPRTJ~max,A~max,V~maxa~=GPRAJ~max,A~max,V~max.These models exactly recover the original datasets and provide smoothly weighted interpolations on the unevaluated dynamic motion parameters. It is able to predict the performance of the dynamic motion parameters using (28). Based on the GPR surrogate model, one can perform global optimization to reduce frame vibration.


### 4.4. Surrogate Model-Based Dynamic Motion Parameter Optimization

#### 4.4.1. Objectives and Priorities

The goal of the optimization problem is to reduce the frame vibration without sacrificing robot efficiency. There are two indices evaluating the performance: the cycle timeT and the frame vibration amplitude a.For efficiency consideration,T must be smaller than the threshold value Td. Therefore, the first constraint is(31)T≤TdThe other goal is to minimize the frame vibration, which can be expressed using(32)p∗=minp⁡aIn general multiobjective optimization problems, usually different objectives may conflict with each other. However, for this problem, those two goals, keeping efficient (31) and reducing vibration (32), have different priorities and requirements. It follows the following two rules:(1)
If a solution cannot meet the cycle time constraint (31), it is bad no matter how small the vibration is.(2)
If two solutions both satisfyT1<T2<Td, the one with smaller vibration is better.In other words, the optimal solution is neither the fastest nor the vibration-free one.


#### 4.4.2. Unified Fitness Function

These two objectives have different priorities. Keeping efficiency is the primary goal while minimizing the vibration is the secondary one. Simply summing them together (ρ=w1T+w2a ) is not applicable because it confuses each index and leads to undesired extreme results.To avoid handling multiobjective optimization problem, the following fitness function is proposed to cooperate the above two rules into one goal,(33)ρJ~max,A~max,V~max∣T~d=2×stepT~-T~d+a~=2×stepGPRTJ~max,A~max,V~max-T~d+GPRAJ~max,A~max,V~maxwhere step(x) is defined as(34)stepx=1x>00x≤0Clearly, the proposed rule (33) mixes the two objectives together in a right way. Because a is normalized within range [-1,1], the first item 2×step(T~-T~d) dominates (33) if step(x)=1. The second item a~ cannot surpass the first item, which satisfies rule (1). Meanwhile, parameter candidates with T≤Td is treated equally and their differences are evaluated by the second item a~, which satisfies rule (2).Equation (34) divided the parameters into two classes. For parameters violates rule (1), (33) maps them in the range of [1,3]. For parameters satisfies rule (1), (33) maps them in the range of [-1,1]. Therefore, the “good” parameters are distinguishable. Existence of the solution to this problem can be guaranteed.The only hyperparameter in the fitness function isTd. It is assigned based on practical requirements. The Delta robot can operate at very high speed in pick and place applications. Cycle time of state-of-the-art Delta robot is 0.2s-0.6s (varies according to particular testing configuration), i.e., 2Hz-5Hz. However, not all applications need such high speed. It is necessary to optimize the robot behavior according to real requirements, such as the pick and place path configuration and the workpiece supply rate. The first one is identified as P0,P1,P2,P3 which are modeled implicitly in (30). The second one is identified as required cycle time Td and introduces a constraint T<Td for the optimization process.


#### 4.4.3. Parameter Optimization

The fitness function (33) formulates the following single objective optimization problem:(35)V~max∗,A~max∗,J~max∗=argminV~max,A~max,J~max⁡ρJ~max,A~max,V~max∣T~dWith the surrogate model and the fitness function, it is able to embed the model into an optimization tool such as Genetic Algorithm (GA)[31]. The grid search method also works but there is a trade-off between accuracy and efficiency. Figure 7 shows diagram of the optimization process.Figure 7
The surrogate model-based parameter optimization problem.The GA has been well implemented in MATLAB [32]. The real encoding GA is chosen and the population size is set as 50. “Stochastic uniform”, “scatter”, and “mutationadaptfeasible” are used to realize the selection, crossover, and mutation operation. The crossover rate is 0.8 and error tolerance is 1e-6. Starting from a group of randomly initialized parameter candidates, the evaluation and regeneration process are iteratively executed until the terminate criterion condition is satisfied.Remark 2.
Because GPR is a nonparametric regression tool that heavily relies on datasets and its computational complexity isOn3, it cannot be applied online when n is large. On the other hand, accuracy of the surrogate model cannot be guaranteed if size of the datasets is too small. The signal acquisition, modeling, and optimization process is performed offline (before actual production) to find better dynamic motion parameters which are suitable for practical production requirements (PPM trajectory, cycle time).
In this paper, GA is used to solve (35). Other optimization methods such as Conjugate Gradient Algorithm can also be integrated.


## 4.1. Data Acquisition

To obtain a surrogate model for (14), all related variables should be able to be extracted from the recorded datasets. According to (14), robot motion and frame vibration are both relevant. Therefore, acceleration profile of robot motion and frame vibration is recorded simultaneously. The recorded signals are time series a(k)={ar(k),af(k)}, where ar(k) is the robot motion while af(k) is the frame vibration.The arguments arep=Vmax,Amax,Jmax, which are known in advance when performing the evaluation experiments. Therefore, the recorded dataset is(16)Dataset=p,ark,afkIn order to minimize the effect of random factors, the robot repeats the pick and place tasks several times for each parameter candidatep.


## 4.2. Signal Processing and Feature Extraction

Performances of the dynamic motion parameters are evaluated from two aspects: efficiency of the robot motion and amplitude of the frame vibration. In this paper, efficiency is measured by cycle timeT. Vibration is measured using maximum acceleration a^. Both indices can be calculated from the datasets.


### 4.2.1. Cycle Time

The measurement of cycle timeT is based on the recorded time series ar(k). By evaluating the acceleration value, it is able to identify whether the robot is moving or staying still. According to the prior knowledge, the robot stays still before new parameters are received. The start time ts can be identified for each movement.Same operation can be used to find the stop timete. Cycle time can be calculated by(17)T=te-ts


### 4.2.2. Maximum Acceleration

Vibration amplitude is evaluated using maximum acceleration, which is calculated using(18)a^=max⁡afk.The surrogate numerical model is obtained by performing regression over the datasets. The dynamic motion parameters areJmax, Amax, and Vmax and the performance indices are a^,T. However, the unit and value range of those variables cause serious problems. Typically, their values have significant differences, ranging from 0-2000000. Without normalization, some of the variables may be submerged by others and the modeling process may fail. To treat each variable equally, a normalization operation is required. Take T as example; the normalizing process follows(19)T~=T-μTσTwhere T~ refers to the normalized T; μ(T) and σ(T) stand for the mean and standard division of T. This operation can make all the variables zero mean and with identical variance. The original variable can be recovered easily through(20)T=T~×σT+μTSame operations can be performed on a^ and Jmax, Amax, and Vmax. After normalization, inputs and outputs of the model are converted to J~max,A~max,V~max, and T~,a~. The latent models can be represented using(21)T~=fTJ~max,A~max,V~maxa~=faJ~max,A~max,V~maxwhere fa,fT are the underlying functions.


## 4.2.1. Cycle Time

The measurement of cycle timeT is based on the recorded time series ar(k). By evaluating the acceleration value, it is able to identify whether the robot is moving or staying still. According to the prior knowledge, the robot stays still before new parameters are received. The start time ts can be identified for each movement.Same operation can be used to find the stop timete. Cycle time can be calculated by(17)T=te-ts


## 4.2.2. Maximum Acceleration

Vibration amplitude is evaluated using maximum acceleration, which is calculated using(18)a^=max⁡afk.The surrogate numerical model is obtained by performing regression over the datasets. The dynamic motion parameters areJmax, Amax, and Vmax and the performance indices are a^,T. However, the unit and value range of those variables cause serious problems. Typically, their values have significant differences, ranging from 0-2000000. Without normalization, some of the variables may be submerged by others and the modeling process may fail. To treat each variable equally, a normalization operation is required. Take T as example; the normalizing process follows(19)T~=T-μTσTwhere T~ refers to the normalized T; μ(T) and σ(T) stand for the mean and standard division of T. This operation can make all the variables zero mean and with identical variance. The original variable can be recovered easily through(20)T=T~×σT+μTSame operations can be performed on a^ and Jmax, Amax, and Vmax. After normalization, inputs and outputs of the model are converted to J~max,A~max,V~max, and T~,a~. The latent models can be represented using(21)T~=fTJ~max,A~max,V~maxa~=faJ~max,A~max,V~maxwhere fa,fT are the underlying functions.


## 4.3. Gaussian Process Regression Based Modeling

Considering the uncertainties in datasets, numbers of samples, and nonlinearity of the underlying function, Gaussian Process Regression (GPR) is used to modeling the datasets. GPR is a nonparametric tool that can handle the modeling problem with noisy observations and system uncertainties [29]. It has been applied in robot process modeling and parameter optimization [30].For the modeling problem (21), there is(22)x=J~max,A~max,V~maxTy=a~=fax,y=T~=fTxThe goal is to derive a model to approximate the underlying functionfax and fTx. Because there are many random factors in the robot system and the observations are also noisy, they can be seen as two Gaussian Processes and can be modeled using GPR.For a Gaussian Processf(x), a set of multivariate Gaussian random variables F={f(x1),f(x2),, ⋯,f(xN)} can be defined over X, where xi∈Xi=1N are the inputs and f(xi) are the corresponding outputs. X is the parameter space which is defined over RD. Instead of deriving an explicit form for f(x), it is approximated by its covariance function k(x,x′): f(x)~GP(0,k(x,x′,θ) where x and x′ are two arbitrary variables in X.The goal of model construction is to find a covariance functionk(x,x′,θ) which fits the data set (X,y) best. Suppose f(x) is a candidate latent function (f in short), the marginal likelihood of H and θ given the data set (X,y) can be derived by performing marginalization over the function f(23)py∣X,H,θ=∫py∣X,f,H,θpf∣X,H,θdfwhere H is the hypothesis on the structure of the covariance function, θ are the hyperparameters, X,y are the sample data sets, and(24)pf∣X,y,H,θ=py∣X,f,H,θpf∣X,H,θpy∣X,H,θis the posterior probability.In order to simplify the modeling process, a log operation is performed to convert multiplication into addition. The Gaussian assumption makes it possible to derive an analytic solution for the marginal log likelihood(25)log⁡py∣X,H,θ=-12yTK+σn2I-1y-12log⁡K+σn2I-n2log⁡2πThen the hyperparameters θ can be found by maximizing the marginal log likelihood function.(26)θ∗=argmaxθ⁡log⁡pyX,H,θThis optimization problem is solved by Conjugate Gradient Algorithm. After training, the model is constructed and can be used for prediction.The joint distribution of the modeling dataset(X,y) and test dataset (X∗,y∗) is(27)yy∗~N0,KX,X+σn2IKX,X∗KX∗,XKX∗,X∗where KX′,X′′ (X′ and X′′ refer to X and X∗) is the covariance matrix whose element Kij at ith row and jth column equals to k(xi′,xj′′). By deriving the conditional distribution, the predictive function is obtained(28)μy∗=Ey∗∣X,y,X∗=KX∗,XKX,X+σn2I-1yVy∗=KX∗,X∗-KX∗,XKX,X+σn2I-1KX,X∗where μ(y∗) is the predicted mean for y∗ and Vy∗ is the variance.In summery, a covariance function should be specified for the GPR model. In the proposed data acquisition procedure, the pick and place motion are repeated for several times for each parameter combination. The random noise in the extracted feature is averaged. Therefore the following covariance function is used:(29)kx,x′=∑i=12σi2exp⁡-∑j=13xj-xj′22lijwhere two Squared Exponential covariance functions with Automatic Relevance Determination (SEARD) items are used to provide feasibility in different length scales lij. Considering the amplitude σi, there are 8 hyperparameters. Based on (29) and the obtained dataset (22), the optimal hyperparameters are found by optimizing (26). Finally two models are created.(30)T~=GPRTJ~max,A~max,V~maxa~=GPRAJ~max,A~max,V~max.These models exactly recover the original datasets and provide smoothly weighted interpolations on the unevaluated dynamic motion parameters. It is able to predict the performance of the dynamic motion parameters using (28). Based on the GPR surrogate model, one can perform global optimization to reduce frame vibration.


## 4.4. Surrogate Model-Based Dynamic Motion Parameter Optimization

### 4.4.1. Objectives and Priorities

The goal of the optimization problem is to reduce the frame vibration without sacrificing robot efficiency. There are two indices evaluating the performance: the cycle timeT and the frame vibration amplitude a.For efficiency consideration,T must be smaller than the threshold value Td. Therefore, the first constraint is(31)T≤TdThe other goal is to minimize the frame vibration, which can be expressed using(32)p∗=minp⁡aIn general multiobjective optimization problems, usually different objectives may conflict with each other. However, for this problem, those two goals, keeping efficient (31) and reducing vibration (32), have different priorities and requirements. It follows the following two rules:(1)
If a solution cannot meet the cycle time constraint (31), it is bad no matter how small the vibration is.(2)
If two solutions both satisfyT1<T2<Td, the one with smaller vibration is better.In other words, the optimal solution is neither the fastest nor the vibration-free one.


### 4.4.2. Unified Fitness Function

These two objectives have different priorities. Keeping efficiency is the primary goal while minimizing the vibration is the secondary one. Simply summing them together (ρ=w1T+w2a ) is not applicable because it confuses each index and leads to undesired extreme results.To avoid handling multiobjective optimization problem, the following fitness function is proposed to cooperate the above two rules into one goal,(33)ρJ~max,A~max,V~max∣T~d=2×stepT~-T~d+a~=2×stepGPRTJ~max,A~max,V~max-T~d+GPRAJ~max,A~max,V~maxwhere step(x) is defined as(34)stepx=1x>00x≤0Clearly, the proposed rule (33) mixes the two objectives together in a right way. Because a is normalized within range [-1,1], the first item 2×step(T~-T~d) dominates (33) if step(x)=1. The second item a~ cannot surpass the first item, which satisfies rule (1). Meanwhile, parameter candidates with T≤Td is treated equally and their differences are evaluated by the second item a~, which satisfies rule (2).Equation (34) divided the parameters into two classes. For parameters violates rule (1), (33) maps them in the range of [1,3]. For parameters satisfies rule (1), (33) maps them in the range of [-1,1]. Therefore, the “good” parameters are distinguishable. Existence of the solution to this problem can be guaranteed.The only hyperparameter in the fitness function isTd. It is assigned based on practical requirements. The Delta robot can operate at very high speed in pick and place applications. Cycle time of state-of-the-art Delta robot is 0.2s-0.6s (varies according to particular testing configuration), i.e., 2Hz-5Hz. However, not all applications need such high speed. It is necessary to optimize the robot behavior according to real requirements, such as the pick and place path configuration and the workpiece supply rate. The first one is identified as P0,P1,P2,P3 which are modeled implicitly in (30). The second one is identified as required cycle time Td and introduces a constraint T<Td for the optimization process.


### 4.4.3. Parameter Optimization

The fitness function (33) formulates the following single objective optimization problem:(35)V~max∗,A~max∗,J~max∗=argminV~max,A~max,J~max⁡ρJ~max,A~max,V~max∣T~dWith the surrogate model and the fitness function, it is able to embed the model into an optimization tool such as Genetic Algorithm (GA)[31]. The grid search method also works but there is a trade-off between accuracy and efficiency. Figure 7 shows diagram of the optimization process.Figure 7
The surrogate model-based parameter optimization problem.The GA has been well implemented in MATLAB [32]. The real encoding GA is chosen and the population size is set as 50. “Stochastic uniform”, “scatter”, and “mutationadaptfeasible” are used to realize the selection, crossover, and mutation operation. The crossover rate is 0.8 and error tolerance is 1e-6. Starting from a group of randomly initialized parameter candidates, the evaluation and regeneration process are iteratively executed until the terminate criterion condition is satisfied.Remark 2.
Because GPR is a nonparametric regression tool that heavily relies on datasets and its computational complexity isOn3, it cannot be applied online when n is large. On the other hand, accuracy of the surrogate model cannot be guaranteed if size of the datasets is too small. The signal acquisition, modeling, and optimization process is performed offline (before actual production) to find better dynamic motion parameters which are suitable for practical production requirements (PPM trajectory, cycle time).
In this paper, GA is used to solve (35). Other optimization methods such as Conjugate Gradient Algorithm can also be integrated.


## 4.4.1. Objectives and Priorities

The goal of the optimization problem is to reduce the frame vibration without sacrificing robot efficiency. There are two indices evaluating the performance: the cycle timeT and the frame vibration amplitude a.For efficiency consideration,T must be smaller than the threshold value Td. Therefore, the first constraint is(31)T≤TdThe other goal is to minimize the frame vibration, which can be expressed using(32)p∗=minp⁡aIn general multiobjective optimization problems, usually different objectives may conflict with each other. However, for this problem, those two goals, keeping efficient (31) and reducing vibration (32), have different priorities and requirements. It follows the following two rules:(1)
If a solution cannot meet the cycle time constraint (31), it is bad no matter how small the vibration is.(2)
If two solutions both satisfyT1<T2<Td, the one with smaller vibration is better.In other words, the optimal solution is neither the fastest nor the vibration-free one.


## 4.4.2. Unified Fitness Function

These two objectives have different priorities. Keeping efficiency is the primary goal while minimizing the vibration is the secondary one. Simply summing them together (ρ=w1T+w2a ) is not applicable because it confuses each index and leads to undesired extreme results.To avoid handling multiobjective optimization problem, the following fitness function is proposed to cooperate the above two rules into one goal,(33)ρJ~max,A~max,V~max∣T~d=2×stepT~-T~d+a~=2×stepGPRTJ~max,A~max,V~max-T~d+GPRAJ~max,A~max,V~maxwhere step(x) is defined as(34)stepx=1x>00x≤0Clearly, the proposed rule (33) mixes the two objectives together in a right way. Because a is normalized within range [-1,1], the first item 2×step(T~-T~d) dominates (33) if step(x)=1. The second item a~ cannot surpass the first item, which satisfies rule (1). Meanwhile, parameter candidates with T≤Td is treated equally and their differences are evaluated by the second item a~, which satisfies rule (2).Equation (34) divided the parameters into two classes. For parameters violates rule (1), (33) maps them in the range of [1,3]. For parameters satisfies rule (1), (33) maps them in the range of [-1,1]. Therefore, the “good” parameters are distinguishable. Existence of the solution to this problem can be guaranteed.The only hyperparameter in the fitness function isTd. It is assigned based on practical requirements. The Delta robot can operate at very high speed in pick and place applications. Cycle time of state-of-the-art Delta robot is 0.2s-0.6s (varies according to particular testing configuration), i.e., 2Hz-5Hz. However, not all applications need such high speed. It is necessary to optimize the robot behavior according to real requirements, such as the pick and place path configuration and the workpiece supply rate. The first one is identified as P0,P1,P2,P3 which are modeled implicitly in (30). The second one is identified as required cycle time Td and introduces a constraint T<Td for the optimization process.


## 4.4.3. Parameter Optimization

The fitness function (33) formulates the following single objective optimization problem:(35)V~max∗,A~max∗,J~max∗=argminV~max,A~max,J~max⁡ρJ~max,A~max,V~max∣T~dWith the surrogate model and the fitness function, it is able to embed the model into an optimization tool such as Genetic Algorithm (GA)[31]. The grid search method also works but there is a trade-off between accuracy and efficiency. Figure 7 shows diagram of the optimization process.Figure 7
The surrogate model-based parameter optimization problem.The GA has been well implemented in MATLAB [32]. The real encoding GA is chosen and the population size is set as 50. “Stochastic uniform”, “scatter”, and “mutationadaptfeasible” are used to realize the selection, crossover, and mutation operation. The crossover rate is 0.8 and error tolerance is 1e-6. Starting from a group of randomly initialized parameter candidates, the evaluation and regeneration process are iteratively executed until the terminate criterion condition is satisfied.Remark 2.
Because GPR is a nonparametric regression tool that heavily relies on datasets and its computational complexity isOn3, it cannot be applied online when n is large. On the other hand, accuracy of the surrogate model cannot be guaranteed if size of the datasets is too small. The signal acquisition, modeling, and optimization process is performed offline (before actual production) to find better dynamic motion parameters which are suitable for practical production requirements (PPM trajectory, cycle time).
In this paper, GA is used to solve (35). Other optimization methods such as Conjugate Gradient Algorithm can also be integrated.


## 5. Experiments and Discussion

### 5.1. Experimental Platform and Configuration

The experiments are performed on the Delta robot shown in Figure8. It has three translational DOFs along the XYZ axis.Figure 8
The robot experiment platform.Its physical parameters are static platform radius (Ri=146mm), input link length (bi=343mm), shaft offset(ci=37mm), output link length (di=953mm), moving platform radius (ri=47.93mm), and end effector offset ei=10mm, where i=1,2,3 refers to the chain numbers. The work range is 1000mm. The robot can run at 5m/s with maximum acceleration of 100m/s2.The robot is controlled by a KEBA Motion Controller (CP-263). It is preprogrammed to repeat a pick and place motion with different dynamic motion parameters, which are sent from the computer through Ethernet interface. The pick and place motion is defined by four points P0-P3 as shown in Figure2. Their coordinates are P0=[0mm,-150mm,-960mm], P1=[0mm,-150mm,-860mm], P2=[0mm,150mm,-860mm], and P3=[0mm,150mm,-960mm]. In each experiment, the parameters and the corresponding vibration signals are recorded.The vibration measurement platform is shown in Figure9.Figure 9
The vibration measurement platform.Vibrations can be measured by several kinds of devices, such as laser tracker, IMU, strain sensor, and even radio signals [33]. In this paper, a nine axis IMU module is used to measure the frame vibration. The used module is based on the InvenSense MPU9250 chip which consisted of a 3-axis gyroscope, a 3-axis accelerometer, and a 3-axis magnetometer. It is an easy-to-use and low cost IMU device with acceptable speed and precision. Two IMU modules are integrated. The #1 IMU sensor is mounted on the robot frame as shown in Figure 8 to detect the frame vibration signals. It is located between the frame center and the frame boundary. The location is rough because what really matters is the relative amplitude of the vibration. The #2 IMU sensor is attached on the moving platform to measure the reference driving signals. Y-axis of both sensors is aligned to the Y-axis of the robot, which is the primary direction of the pick and place motion. Both sensors transfer the data at 200Hz back to computer through USB-UART interfaces.Figure10 gives the overall flowchart of the signal acquisition process.Figure 10
Flowchart of the vibration signal acquisition process.Each parameter candidate is repeated for four times to reduce effect of noises and random factors. Figure11 shows a typical vibration datasets.Figure 11
Typical vibration signals recorded by the IMU sensors. (a) Acceleration of moving platform along Y-axis; (b) frame vibration along Y-axis.
(a) (b)Each sample lasts for 6 seconds and contains 4 pick and place cycles. Figure11(a) is the acceleration profile of the robot’s moving platform. Figure 11(b) is the frame vibration signal. Both signals are in Y-axis, i.e., the main direction of the pick and place motion. From this figure, it is able to identify that the robot performs the pick and place motion for four times. During the motion process, the frame is forced to vibrate. Even after the robot stops at 3.5s, the frame still vibrates for nearly 3 seconds.Dynamic motion parameters that may affect the smoothness and vibration are listed in Table1. The maximum, minimum, and interval of each variable are presented.Table 1
Motion related parameters.
V m a x (mm/s) A m a x (mm/s2) J m a x (mm/s3) Min 1000 10000 1000000 Max 5000 90000 25000000 Interval 500 20000 250000Totally9×5×7=315 experiments were conducted.


### 5.2. Experimental Results and Analysis

#### 5.2.1. Data Processing and Modeling

After signal processing, performance indices of all the 315 data sets are extracted. Cycle timeT and vibration amplitude a^ are plotted in Figure 12. The horizontal axis is cycle time while the vertical axis is the vibration. The blue curve is a fitted polynomial for comparison. From this figure, one can find that the vibration is inversely proportional to the cycle time in general. This is consistent with (13) because T is in the numerator. Meanwhile, significant variances exist in the vibration amplitude even for same cycle time. Take T=0.65s, for example, the maximum vibration is nearly double of the minimum one. This indicates the possibility of decreasing the frame vibration while maintaining the efficiency.Figure 12
Scatter plot of the processed samples.Figure13 shows a slice of the generated GPR model with Vmax=3000mm/s. From Figure 13, the saturation phenomenon (cycle time cannot be further reduced if velocity is limited) and vibration variability (“peaks” and “valleys”) are clearly demonstrated. The goal of this paper is to find those “valleys" and refine the system performance.Figure 13
Slice of the generated model.
(a)
Cycle time (b)
Vibration


#### 5.2.2. Optimization

For the optimization problem defined using (33) and (35), Genetic Algorithm (GA) is used to perform the optimization process. Given desired cycle time Td, the corresponding optimal dynamic parameters can be found. GA is configured with the settings presented in Section 4.4.3. The optimization process is shown in Figure 14.Figure 14
GA optimization process.It can be seen that the algorithm converges after about 10 generations. The optimal fitness values are all within range[-1,1], which indicates the optimal parameters all satisfy the cycle time constraint (31). After optimization and recovery, the optimal parameters and the corresponding performance indices are listed in Table 2.Table 2
The optimized dynamic parameters.
T d/s V m a x A m a x J m a x T/s a ^ real T reala^ 0.6 3100 50000 2235500 0.59 0.70 0.59 0.71 0.7 2500 50000 1703200 0.66 0.53 0.67 0.54 0.8 3000 70200 1000600 0.77 0.43 0.77 0.42Applying these optimized dynamic motion parameters on the experiment platform and recording the vibration, one can find that the predicted values are close to the real ones. The results verify effectiveness of the proposed modeling and optimization method. Waveform of the optimal parameter forTd=0.8s is shown in Figure 18.


### 5.3. Discussion

As discussed in Section3, amplitude of the frame vibration is affected by two factors: frequency response characteristics of the frame and frequency spectrum of the robot motion’s acceleration profile.


#### 5.3.1. Frame Frequency Response Characteristics

From Figure11, it is able to find that when the external excitation disappears, the frame still vibrate for several seconds. Applying the Fourier transformation, spectrum of the robot motion and frame vibration can be obtained. Dividing the corresponding frequency component, one gets the frequency response chart as shown in Figure 15.Figure 15
The frame frequency response characteristics.From this figure, one can find that the resonance frequency lies around 7.6Hz, from 6.76Hz to 9.2Hz.


#### 5.3.2. Spectrum of the Robot Motion Acceleration Profile

According to Figure12, cycle time ranges from 0.56s to 1.2s, i.e., the base frequency is less than 1.8Hz and the 5th harmonics component lies in the range of [6.76,9.2].The critical base frequency isfbc=6.76÷5=1.35Hz, i.e., Tc=0.74s. For this critical case, there is(36)6.75=5fb<6.769.2<7fb=9.45Therefore, the resonance frequency is avoided by both the 5th and 7th harmonics, and consequently the vibration is small. The corresponding datasets are labeled in Figure 12.For cases withT>Tc, the 5th harmonics are outside the resonance range. Although the 7th harmonics are within the resonance range, typically, their amplitudes are smaller than the 5th one. Therefore, the frame vibration is also small. For cases with T<Tc, because of resonance, the vibration increases significantly. The frame vibration is dominated by the amplitudes of the 5th harmonics. The above conclusion is consistent with the results given in Figure 12.To verify the influence of 5th harmonics, a comparative example is given in Figure16.Figure 16
Comparative results of two cases with same cycle timeT=0.65s. However, the frame vibration amplitudes have significant difference: 0.7353 m/s2 and 1.2286 m/s2.
(a)
Spectrum of the robot motion acceleration profile (b)
Spectrum of the frame vibrationTwo groups of dynamics parameters[Vmax,Amax,Jmax] =[2000,50000,2500000] (Case B:Blue solid line) and [3000,50000,1750000](Case R:Red dash line) are examined. Both cases have same cycle time 0.65s. The harmonic distributions are identical (because some of the accelerations are beyond the sensor range, spurious harmonics also exist in the spectrum).Comparing the amplitude of each harmonics, Case B is larger than Case R except the 5th harmonics. However, frame vibration of Case B is smaller than Case R, which proves importance of the amplitude of the close-to-natural-frequency.Figure17 shows the spectrum of the optimal parameters listed in Table 2. One can find that both of the above two conditions are satisfied. The 5th and 7th harmonics are away from the resonance frequency and also the amplitude is low. Therefore, the frame vibration is minimized.Figure 17
Spectrum of the acceleration profile generated by the optimal dynamic motion parameters.Figure 18
Waveforms of the comparative group of parameters.The waveform of those three groups of parameters is shown in Figure18. It is clearly seen that the dynamic motion parameter optimization process indeed reduces the frame vibration.


#### 5.3.3. General Applicability of the Proposed Methods

The central topic of this paper is to find a way to optimize dynamic motion parameters of PPM motion so as to reduce the induced robot base/frame vibration. As shown in Figure6, the proposed method uses two vibration datasets obtained from robot body and frame to build a model characterizing the relation between dynamic motion parameters and system performance. Then based on the model, an optimization problem can be formulated and solved to find the best parameter.From Figure6, one can find that only the “evaluating" process is related to the actual robot system. Although the derivation and experiments are based on Delta robot, it can be easily applied to other robots with similar PPM motions, such as SCARA robot and palletizing robot. However, because cycle time of non-Delta robots is typically longer and they are usually mounted on the ground/table, the frame vibration is not so obvious.It is easy to apply the proposed method to a different robot. The only modification is to write a simple robot program that uses different parameters to perform PPM motion or just add several instructions to the existing program to modify the dynamic parameters. Then the dataset can be recorded automatically and processed offline.


## 5.1. Experimental Platform and Configuration

The experiments are performed on the Delta robot shown in Figure8. It has three translational DOFs along the XYZ axis.Figure 8
The robot experiment platform.Its physical parameters are static platform radius (Ri=146mm), input link length (bi=343mm), shaft offset(ci=37mm), output link length (di=953mm), moving platform radius (ri=47.93mm), and end effector offset ei=10mm, where i=1,2,3 refers to the chain numbers. The work range is 1000mm. The robot can run at 5m/s with maximum acceleration of 100m/s2.The robot is controlled by a KEBA Motion Controller (CP-263). It is preprogrammed to repeat a pick and place motion with different dynamic motion parameters, which are sent from the computer through Ethernet interface. The pick and place motion is defined by four points P0-P3 as shown in Figure2. Their coordinates are P0=[0mm,-150mm,-960mm], P1=[0mm,-150mm,-860mm], P2=[0mm,150mm,-860mm], and P3=[0mm,150mm,-960mm]. In each experiment, the parameters and the corresponding vibration signals are recorded.The vibration measurement platform is shown in Figure9.Figure 9
The vibration measurement platform.Vibrations can be measured by several kinds of devices, such as laser tracker, IMU, strain sensor, and even radio signals [33]. In this paper, a nine axis IMU module is used to measure the frame vibration. The used module is based on the InvenSense MPU9250 chip which consisted of a 3-axis gyroscope, a 3-axis accelerometer, and a 3-axis magnetometer. It is an easy-to-use and low cost IMU device with acceptable speed and precision. Two IMU modules are integrated. The #1 IMU sensor is mounted on the robot frame as shown in Figure 8 to detect the frame vibration signals. It is located between the frame center and the frame boundary. The location is rough because what really matters is the relative amplitude of the vibration. The #2 IMU sensor is attached on the moving platform to measure the reference driving signals. Y-axis of both sensors is aligned to the Y-axis of the robot, which is the primary direction of the pick and place motion. Both sensors transfer the data at 200Hz back to computer through USB-UART interfaces.Figure10 gives the overall flowchart of the signal acquisition process.Figure 10
Flowchart of the vibration signal acquisition process.Each parameter candidate is repeated for four times to reduce effect of noises and random factors. Figure11 shows a typical vibration datasets.Figure 11
Typical vibration signals recorded by the IMU sensors. (a) Acceleration of moving platform along Y-axis; (b) frame vibration along Y-axis.
(a) (b)Each sample lasts for 6 seconds and contains 4 pick and place cycles. Figure11(a) is the acceleration profile of the robot’s moving platform. Figure 11(b) is the frame vibration signal. Both signals are in Y-axis, i.e., the main direction of the pick and place motion. From this figure, it is able to identify that the robot performs the pick and place motion for four times. During the motion process, the frame is forced to vibrate. Even after the robot stops at 3.5s, the frame still vibrates for nearly 3 seconds.Dynamic motion parameters that may affect the smoothness and vibration are listed in Table1. The maximum, minimum, and interval of each variable are presented.Table 1
Motion related parameters.
V m a x (mm/s) A m a x (mm/s2) J m a x (mm/s3) Min 1000 10000 1000000 Max 5000 90000 25000000 Interval 500 20000 250000Totally9×5×7=315 experiments were conducted.


## 5.2. Experimental Results and Analysis

### 5.2.1. Data Processing and Modeling

After signal processing, performance indices of all the 315 data sets are extracted. Cycle timeT and vibration amplitude a^ are plotted in Figure 12. The horizontal axis is cycle time while the vertical axis is the vibration. The blue curve is a fitted polynomial for comparison. From this figure, one can find that the vibration is inversely proportional to the cycle time in general. This is consistent with (13) because T is in the numerator. Meanwhile, significant variances exist in the vibration amplitude even for same cycle time. Take T=0.65s, for example, the maximum vibration is nearly double of the minimum one. This indicates the possibility of decreasing the frame vibration while maintaining the efficiency.Figure 12
Scatter plot of the processed samples.Figure13 shows a slice of the generated GPR model with Vmax=3000mm/s. From Figure 13, the saturation phenomenon (cycle time cannot be further reduced if velocity is limited) and vibration variability (“peaks” and “valleys”) are clearly demonstrated. The goal of this paper is to find those “valleys" and refine the system performance.Figure 13
Slice of the generated model.
(a)
Cycle time (b)
Vibration


### 5.2.2. Optimization

For the optimization problem defined using (33) and (35), Genetic Algorithm (GA) is used to perform the optimization process. Given desired cycle time Td, the corresponding optimal dynamic parameters can be found. GA is configured with the settings presented in Section 4.4.3. The optimization process is shown in Figure 14.Figure 14
GA optimization process.It can be seen that the algorithm converges after about 10 generations. The optimal fitness values are all within range[-1,1], which indicates the optimal parameters all satisfy the cycle time constraint (31). After optimization and recovery, the optimal parameters and the corresponding performance indices are listed in Table 2.Table 2
The optimized dynamic parameters.
T d/s V m a x A m a x J m a x T/s a ^ real T reala^ 0.6 3100 50000 2235500 0.59 0.70 0.59 0.71 0.7 2500 50000 1703200 0.66 0.53 0.67 0.54 0.8 3000 70200 1000600 0.77 0.43 0.77 0.42Applying these optimized dynamic motion parameters on the experiment platform and recording the vibration, one can find that the predicted values are close to the real ones. The results verify effectiveness of the proposed modeling and optimization method. Waveform of the optimal parameter forTd=0.8s is shown in Figure 18.


## 5.2.1. Data Processing and Modeling

After signal processing, performance indices of all the 315 data sets are extracted. Cycle timeT and vibration amplitude a^ are plotted in Figure 12. The horizontal axis is cycle time while the vertical axis is the vibration. The blue curve is a fitted polynomial for comparison. From this figure, one can find that the vibration is inversely proportional to the cycle time in general. This is consistent with (13) because T is in the numerator. Meanwhile, significant variances exist in the vibration amplitude even for same cycle time. Take T=0.65s, for example, the maximum vibration is nearly double of the minimum one. This indicates the possibility of decreasing the frame vibration while maintaining the efficiency.Figure 12
Scatter plot of the processed samples.Figure13 shows a slice of the generated GPR model with Vmax=3000mm/s. From Figure 13, the saturation phenomenon (cycle time cannot be further reduced if velocity is limited) and vibration variability (“peaks” and “valleys”) are clearly demonstrated. The goal of this paper is to find those “valleys" and refine the system performance.Figure 13
Slice of the generated model.
(a)
Cycle time (b)
Vibration


## 5.2.2. Optimization

For the optimization problem defined using (33) and (35), Genetic Algorithm (GA) is used to perform the optimization process. Given desired cycle time Td, the corresponding optimal dynamic parameters can be found. GA is configured with the settings presented in Section 4.4.3. The optimization process is shown in Figure 14.Figure 14
GA optimization process.It can be seen that the algorithm converges after about 10 generations. The optimal fitness values are all within range[-1,1], which indicates the optimal parameters all satisfy the cycle time constraint (31). After optimization and recovery, the optimal parameters and the corresponding performance indices are listed in Table 2.Table 2
The optimized dynamic parameters.
T d/s V m a x A m a x J m a x T/s a ^ real T reala^ 0.6 3100 50000 2235500 0.59 0.70 0.59 0.71 0.7 2500 50000 1703200 0.66 0.53 0.67 0.54 0.8 3000 70200 1000600 0.77 0.43 0.77 0.42Applying these optimized dynamic motion parameters on the experiment platform and recording the vibration, one can find that the predicted values are close to the real ones. The results verify effectiveness of the proposed modeling and optimization method. Waveform of the optimal parameter forTd=0.8s is shown in Figure 18.


## 5.3. Discussion

As discussed in Section3, amplitude of the frame vibration is affected by two factors: frequency response characteristics of the frame and frequency spectrum of the robot motion’s acceleration profile.


### 5.3.1. Frame Frequency Response Characteristics

From Figure11, it is able to find that when the external excitation disappears, the frame still vibrate for several seconds. Applying the Fourier transformation, spectrum of the robot motion and frame vibration can be obtained. Dividing the corresponding frequency component, one gets the frequency response chart as shown in Figure 15.Figure 15
The frame frequency response characteristics.From this figure, one can find that the resonance frequency lies around 7.6Hz, from 6.76Hz to 9.2Hz.


### 5.3.2. Spectrum of the Robot Motion Acceleration Profile

According to Figure12, cycle time ranges from 0.56s to 1.2s, i.e., the base frequency is less than 1.8Hz and the 5th harmonics component lies in the range of [6.76,9.2].The critical base frequency isfbc=6.76÷5=1.35Hz, i.e., Tc=0.74s. For this critical case, there is(36)6.75=5fb<6.769.2<7fb=9.45Therefore, the resonance frequency is avoided by both the 5th and 7th harmonics, and consequently the vibration is small. The corresponding datasets are labeled in Figure 12.For cases withT>Tc, the 5th harmonics are outside the resonance range. Although the 7th harmonics are within the resonance range, typically, their amplitudes are smaller than the 5th one. Therefore, the frame vibration is also small. For cases with T<Tc, because of resonance, the vibration increases significantly. The frame vibration is dominated by the amplitudes of the 5th harmonics. The above conclusion is consistent with the results given in Figure 12.To verify the influence of 5th harmonics, a comparative example is given in Figure16.Figure 16
Comparative results of two cases with same cycle timeT=0.65s. However, the frame vibration amplitudes have significant difference: 0.7353 m/s2 and 1.2286 m/s2.
(a)
Spectrum of the robot motion acceleration profile (b)
Spectrum of the frame vibrationTwo groups of dynamics parameters[Vmax,Amax,Jmax] =[2000,50000,2500000] (Case B:Blue solid line) and [3000,50000,1750000](Case R:Red dash line) are examined. Both cases have same cycle time 0.65s. The harmonic distributions are identical (because some of the accelerations are beyond the sensor range, spurious harmonics also exist in the spectrum).Comparing the amplitude of each harmonics, Case B is larger than Case R except the 5th harmonics. However, frame vibration of Case B is smaller than Case R, which proves importance of the amplitude of the close-to-natural-frequency.Figure17 shows the spectrum of the optimal parameters listed in Table 2. One can find that both of the above two conditions are satisfied. The 5th and 7th harmonics are away from the resonance frequency and also the amplitude is low. Therefore, the frame vibration is minimized.Figure 17
Spectrum of the acceleration profile generated by the optimal dynamic motion parameters.Figure 18
Waveforms of the comparative group of parameters.The waveform of those three groups of parameters is shown in Figure18. It is clearly seen that the dynamic motion parameter optimization process indeed reduces the frame vibration.


### 5.3.3. General Applicability of the Proposed Methods

The central topic of this paper is to find a way to optimize dynamic motion parameters of PPM motion so as to reduce the induced robot base/frame vibration. As shown in Figure6, the proposed method uses two vibration datasets obtained from robot body and frame to build a model characterizing the relation between dynamic motion parameters and system performance. Then based on the model, an optimization problem can be formulated and solved to find the best parameter.From Figure6, one can find that only the “evaluating" process is related to the actual robot system. Although the derivation and experiments are based on Delta robot, it can be easily applied to other robots with similar PPM motions, such as SCARA robot and palletizing robot. However, because cycle time of non-Delta robots is typically longer and they are usually mounted on the ground/table, the frame vibration is not so obvious.It is easy to apply the proposed method to a different robot. The only modification is to write a simple robot program that uses different parameters to perform PPM motion or just add several instructions to the existing program to modify the dynamic parameters. Then the dataset can be recorded automatically and processed offline.


## 5.3.1. Frame Frequency Response Characteristics

From Figure11, it is able to find that when the external excitation disappears, the frame still vibrate for several seconds. Applying the Fourier transformation, spectrum of the robot motion and frame vibration can be obtained. Dividing the corresponding frequency component, one gets the frequency response chart as shown in Figure 15.Figure 15
The frame frequency response characteristics.From this figure, one can find that the resonance frequency lies around 7.6Hz, from 6.76Hz to 9.2Hz.


## 5.3.2. Spectrum of the Robot Motion Acceleration Profile

According to Figure12, cycle time ranges from 0.56s to 1.2s, i.e., the base frequency is less than 1.8Hz and the 5th harmonics component lies in the range of [6.76,9.2].The critical base frequency isfbc=6.76÷5=1.35Hz, i.e., Tc=0.74s. For this critical case, there is(36)6.75=5fb<6.769.2<7fb=9.45Therefore, the resonance frequency is avoided by both the 5th and 7th harmonics, and consequently the vibration is small. The corresponding datasets are labeled in Figure 12.For cases withT>Tc, the 5th harmonics are outside the resonance range. Although the 7th harmonics are within the resonance range, typically, their amplitudes are smaller than the 5th one. Therefore, the frame vibration is also small. For cases with T<Tc, because of resonance, the vibration increases significantly. The frame vibration is dominated by the amplitudes of the 5th harmonics. The above conclusion is consistent with the results given in Figure 12.To verify the influence of 5th harmonics, a comparative example is given in Figure16.Figure 16
Comparative results of two cases with same cycle timeT=0.65s. However, the frame vibration amplitudes have significant difference: 0.7353 m/s2 and 1.2286 m/s2.
(a)
Spectrum of the robot motion acceleration profile (b)
Spectrum of the frame vibrationTwo groups of dynamics parameters[Vmax,Amax,Jmax] =[2000,50000,2500000] (Case B:Blue solid line) and [3000,50000,1750000](Case R:Red dash line) are examined. Both cases have same cycle time 0.65s. The harmonic distributions are identical (because some of the accelerations are beyond the sensor range, spurious harmonics also exist in the spectrum).Comparing the amplitude of each harmonics, Case B is larger than Case R except the 5th harmonics. However, frame vibration of Case B is smaller than Case R, which proves importance of the amplitude of the close-to-natural-frequency.Figure17 shows the spectrum of the optimal parameters listed in Table 2. One can find that both of the above two conditions are satisfied. The 5th and 7th harmonics are away from the resonance frequency and also the amplitude is low. Therefore, the frame vibration is minimized.Figure 17
Spectrum of the acceleration profile generated by the optimal dynamic motion parameters.Figure 18
Waveforms of the comparative group of parameters.The waveform of those three groups of parameters is shown in Figure18. It is clearly seen that the dynamic motion parameter optimization process indeed reduces the frame vibration.


## 5.3.3. General Applicability of the Proposed Methods

The central topic of this paper is to find a way to optimize dynamic motion parameters of PPM motion so as to reduce the induced robot base/frame vibration. As shown in Figure6, the proposed method uses two vibration datasets obtained from robot body and frame to build a model characterizing the relation between dynamic motion parameters and system performance. Then based on the model, an optimization problem can be formulated and solved to find the best parameter.From Figure6, one can find that only the “evaluating" process is related to the actual robot system. Although the derivation and experiments are based on Delta robot, it can be easily applied to other robots with similar PPM motions, such as SCARA robot and palletizing robot. However, because cycle time of non-Delta robots is typically longer and they are usually mounted on the ground/table, the frame vibration is not so obvious.It is easy to apply the proposed method to a different robot. The only modification is to write a simple robot program that uses different parameters to perform PPM motion or just add several instructions to the existing program to modify the dynamic parameters. Then the dataset can be recorded automatically and processed offline.


## 6. Conclusion

This paper focuses on the frame vibration suppression problem for Delta robot in pick and place application. By analyzing acceleration profile of the pick and place motion in frequency domain, distribution and amplitude of harmonics in the robot motion are identified. They are proved to be internally related to dynamic motion parameters such as jerk, acceleration, and velocity. Guidelines for reducing the vibration are found. In order to solve the practical vibration reduction problem, a surrogate model-based optimization solution is proposed. The model is based on Gaussian Process Regression (GPR) and can be constructed from the recorded acceleration signals on real robot platform. In order to achieve the goal of keeping efficiency and reducing vibration, their priority differences are utilized. They are fused together in a single fitness function. This single objective optimization problem is solved by Genetic Algorithm (GA). The proposed dynamic motion parameter optimization method is flexible and generally applicable and can work with other PPM generation algorithms or acceleration profiles.The future works focus on two main aspects: developing an embed measurement unit with general applicable interface and integrating them on different robot platforms; exploring the possibility of online dynamic optimization problem.



---

*Source: 2945314-2018-06-10.xml*