Research on Active Suppression Methods for End-Effector Residual Vibration of Heavy-Load Collaborative Robots in Arbitrary Poses

1. Introduction

Collaborative robot, characterized by rapid deployment, flexible manufacturing, and human-robot collaboration, meets the needs of agile manufacturing in the manufacturing industry, such as multi-variety production, customization, flexibility, and rapid switching. Compared with traditional industrial robotic arms, collaborative robotic arms initially focused on lightweight application scenarios. However, with the industrial upgrading and development of electric vehicles, 3C manufacturing, precision manufacturing, and 5G communications, collaborative robots with a load exceeding 10Kg are increasingly being applied[1]. The annual market growth rate of collaborative robotic arms has exceeded 50%, far surpassing that of traditional industrial robotic arms, and has become a major trend in the development of robotic arms [2].

Due to factors such as inertia, the material flexibility of the robotic arm's links and joint harmonic reducers, there will still be a period of residual vibration after the robotic arm reaches the end poes. Effectively suppressing residual vibration is crucial to improving the performance of collaborative robots[3,4,5,6]. According to whether external input energy is required, the current vibration suppression methods can be divided into two major categories: passive control and active control[7].

Passive control does not require external input energy. It aims to suppress residual vibration by changing the structural and material characteristics of the system. Bandopadhya [8] utilized the characteristics of high strength and large damping of ionic polymer metal composites to effectively improve the residual vibration problem of flexible cantilever beams. Bajkowski [9] focused on the double-beam system and used intelligent particle structures to form a vibration absorber. Passive control requires additional hardware devices, which not only increases the cost of the robot, but also the system is usually designed for specific vibration frequencies.

Active control methods suppress residual vibrations by superimposing additional control signals through control algorithms to act on the system, resulting in faster response and better vibration suppression effects. Cannon [10] focused on serial flexible robotic arms, implemented feedback closed-loop using strain gauge sensors, designed a linear quadratic Gaussian controller to suppress residual vibrations. Preumont [11] elaborated on various active and passive control methods for residual vibration suppression. Meckl [12] used the common S-curve motion profile but optimized the selection of acceleration and deceleration phase durations. The most representative open-loop feedforward control method is input shaping method[13,14,15,16,17], which evolved from the Posicast control method proposed by Smith [18] in 1957 to suppress residual oscillations in electrical systems. Pao [19,20] developed a hybrid input shaping design method that interpolates between Zero Vibration (ZV), Zero Vibration Derivative (ZVD), and Extra-Insensitive (EI) shaper designs to achieve near-optimal performance. Tuttle [21] designed input shapers using zero placement in the discrete domain to reduce multi-modal vibrations. Zhao[22,23] proposed a zero-time-delay input shaping technique that adjusts input signals without introducing additional time delays. Huang et al. [24,25] focused on heavy-load cranes, using input shaping technology to eliminate load swings caused by human operator commands and employing feedback controllers to reduce the impact of wind gusts. This significantly improved the system's suppression speed for vibrations caused by inter-link coupling and external disturbances.

To achieve effective vibration suppression, it is necessary to obtain the natural frequency and damping ratio parameters of the robot end vibration[26]. Kenderi et al. [27] proposed an optimal observer based on the least square method and Kalman filtering for the identification of nonlinear mechanical systems. Masri et al. [28] proposed a neural network-based nonlinear dynamic system identification method and applied it to a damped Duffing oscillator under deterministic excitation, providing a high-fidelity mathematical model for structurally unknown nonlinear systems encountered in the field of applied mechanics. Juang[29] proposed an algorithm called the Eigensystem Realization Algorithm for identifying modal parameters and simplifying dynamic systems from test data. Combined with singular value decomposition technology, modal parameter identification is carried out to quantitatively identify system modes and noise modes, and it has been verified using experimental data from the Galileo spacecraft. Nadkarni[30], in order to predict the response of a structure when subjected to forces in the working environment, determined the dynamic characteristics of the structure through experimental modal analysis, obtained three modal parameters, namely natural frequency, damping, and modal shape, and performed numerical verification using the finite element software package Hypermesh 13.0. Doughty[31] discussed in detail the advantages and limitations of three nonlinear system identification techniques for cantilever beam analysis, based on the continuous-time differential equation model of the system, the relationships generated using the harmonic balance method, and the fitting of steady-state response data to amplitude and phase modulation equations.

This paper adopts the input shaping method to solve the problem of end residual vibration in heavy-load six-axis collaborative robots. This method realizes the suppression of residual vibration by adjusting the original input control signal with an input shaper composed of several pulse signals with different amplitudes and time delays to obtain a new input signal that acts on the system. A data-based method for estimating the residual vibration frequency and damping ratio of the robotic arm in arbitrary poses is proposed. A model from the robotic arm's joint angles and end load mass to the end vibration frequency and damping ratio is established, and the least square identification method is used to identify the model parameters. Finally, based on the estimated vibration parameters of the robotic arm's end in arbitrary poses, three input shapers are designed and applied to the control of the collaborative robotic arm, achieving a good suppression effect on the end residual vibration of the robotic arm in arbitrary poses. Unlike existing methods requiring pre-motion calibration, our approach enables real-time vibration parameter estimation for arbitrary poses, eliminating the need for trajectory pre-runs. The average suppression ratio is over 40%, and when the parameter identification error is small, the suppression effect reaches more than 70%.

2. Modeling and Parameter Estimation of Residual Vibration at the End of a Manipulator in Arbitrary Poses

2.1. Modeling of Residual Vibration at the End of a Manipulator

A 6-DOF serial manipulator can be approximated by a second-order system, and the transfer function of a typical second-order system is shown in Equation (1)

$$G(s)={\displaystyle \frac{{\omega }_n^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}}$$

(1)

where ${\omega }_n$ is the undamped natural frequency of the system, and $\zeta$is the damping ratio of the system. When a unit impulse is applied at time 0, the response equation is the Laplace inverse transform of its transfer function, as shown in Equation (2)

$$g(t)=-{\displaystyle \frac{{\omega }_n}{\sqrt{1-{\zeta }^2}}}{e}^{-{\zeta }_n{\omega }_{n}t}\sin \left(\sqrt{1-{\zeta }^2}{\omega }_{n}t\right).$$

(2)

A pulse with an amplitude of $A_i$ and a time delay of $t_i$will excite residual vibration in the system, as shown in Equation (3)

$$y_i(t)=A_i{\displaystyle \frac{{\omega }_n}{\sqrt{1-{\zeta }^2}}}{e}^{-{\zeta }_n{\omega }_n(t-t_i)}\sin \left[\sqrt{1-{\zeta }^2}{\omega }_n(t-t_i)\right].$$

(3)

The total response of the system is the sum of the responses induced by all pulses. Assuming there are $n$ pulses, the total response of the system is

$$y(t)={\displaystyle \sum _{i=1}^n{A}_i\frac{{\omega }_n}{\sqrt{1-{\zeta }^2}}{e}^{-\zeta {\omega }_n(t-t_i)}\sin \left[\sqrt{1-{\zeta }^2}{\omega }_n(t-t_i)\right]}.$$

(4)

Using trigonometric function formulas, Equation (4) is expanded and slightly rearranged into the form shown in Equation (5)

$$\begin{array}{ll}y(t)& ={\displaystyle \sum _{i=1}^n\frac{A_i{\omega }_n}{\sqrt{1-{\zeta }^2}}}{e}^{-\zeta {\omega }_n(t-t_i)}\left[\sin ({\omega }_{d}t)\cos ({\omega }_d{t}_i)-\cos ({\omega }_{d}t)\sin ({\omega }_d{t}_i)\right]\\ & ={\displaystyle \frac{{\omega }_n\sin ({\omega }_{d}t)}{\sqrt{1-{\zeta }^2}}}{\displaystyle \sum _{i=1}^n{A}_i}{e}^{-\zeta {\omega }_n(t-t_i)}\cos ({\omega }_d{t}_i)-{\displaystyle \frac{{\omega }_n\cos ({\omega }_{d}t)}{\sqrt{1-{\zeta }^2}}}{\displaystyle \sum _{i=1}^n{A}_i}{e}^{-\zeta {\omega }_n(t-t_i)}\sin ({\omega }_d{t}_i)\\ & ={\displaystyle \frac{{\omega }_n{e}^{-\zeta {\omega }_{n}t}}{\sqrt{1-{\zeta }^2}}}\left(\sin ({\omega }_{d}t){\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\cos ({\omega }_d{t}_i)-\cos ({\omega }_{d}t)\cdot {\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\sin ({\omega }_d{t}_i)\right)\end{array},$$

(5)

where ${\omega }_d$ is the damped natural frequency of the system, and its relationship with the undamped natural frequency ${\omega }_n$ is

$${\omega }_d={\omega }_n\sqrt{1-{\zeta }^2}.$$

(6)

Using the trigonometric auxiliary angle formula, Equation (6) can be further sorted out to obtain Equation (7)

$$y(t)={\displaystyle \frac{{\omega }_n{e}^{-\zeta {\omega }_{n}t}}{\sqrt{1-{\zeta }^2}}}\sqrt{{\left[{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\cos ({\omega }_d{t}_i)\right]}^2+{\left[{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\sin ({\omega }_d{t}_i)\right]}^2}\cdot \sin ({\omega }_{d}t-\varphi )$$

(7)

where $\varphi =\arctan \left({\displaystyle \frac{{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\sin ({\omega }_d{t}_i)}{{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\cos ({\omega }_d{t}_i)}}\right)$.

To describe the effect of input shaping in the system, a residual vibration variable $A_{\mathrm{amp}}$is introduced. Its physical meaning is the ratio of the system vibration amplitude without shaping to that with shaping at the moment when the last pulse of the input shaper is applied. This variable is also called the residual vibration percentage, which is the ratio of Equation (2) to Equation (7), and the simplified form is shown in Equation (8):

$$\begin{array}{ll}{A}_{\mathrm{amp}}& ={\displaystyle \frac{y(t_n)}{g(t_n)}}\\ & ={\displaystyle \frac{\frac{{\omega }_n{e}^{-\zeta {\omega }_n{t}_n}}{\sqrt{1-{\zeta }^2}}\sqrt{{\left[{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\cos ({\omega }_d{t}_i)\right]}^2+{\left[{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\sin ({\omega }_d{t}_i)\right]}^2}}{\frac{{\omega }_n}{\sqrt{1-{\zeta }^2}}{e}^{-\zeta {\omega }_n{t}_n}}}\\ & =\sqrt{{\left[{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\cos ({\omega }_d{t}_i)\right]}^2+{\left[{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\sin ({\omega }_d{t}_i)\right]}^2}\end{array}.$$

(8)

For a heavy-load collaborative robotic manipulator, to achieve the effect of eliminating residual vibration, it is necessary to design a series of appropriate $A_i$ and $t_i$ so that the residual vibration percentage $A_{\mathrm{amp}}$ is as close to 0 as possible.

Figure 1 shows the acceleration signal at the end of the manipulator measured by an accelerometer when the manipulator reaches the end point after moving along a certain trajectory. From Equation (3) and Figure 1, it can be seen that the mathematical form of the residual vibration excited at the end of the robotic arm by a pulse with an amplitude of $A_i$ and a time delay of $t_i$is approximately a damped harmonic vibration, which is expressed by the impulse response equation (9)

$$a(t)={\displaystyle \frac{{\omega }_n}{\sqrt{1-{\zeta }^2}}}{e}^{-\zeta {\omega }_n(t-t_0)}\sin \left(\sqrt{1-{\zeta }^2}{\omega }_n(t-t_0)\right),$$

(9)

where $t_0$is the initial time. As shown in Figure 1, only two data points $a(t_1)$and $a(t_2)$ ($(t_2>t_1)$) are needed to determine ${\omega }_n$ and the damping ratio $\zeta$. ${\omega }_n$ can be directly obtained from (10),

$${\omega }_n={\displaystyle \frac{2\pi }{t_2-t_1}}(t_2>t_1).$$

(10)

After obtaining the natural frequency ${\omega }_n$, the damping ratio $\zeta$ is then calculated. For convenience of calculation, $a(t_1)$ and $a(t_2)$ can be taken as two adjacent peak points, where the sine factor of the peak points is equal to 1, meaning only the amplitude part needs to be considered. Dividing $a(t_1)$ by $a(t_2)$ gives Equation (11):

$${\displaystyle \frac{a(t_1)}{a(t_2)}}=e^{-\zeta \omega (t_1-t_2)}.$$

(11)

Substituting the natural frequency ${\omega }_n$ obtained from Equation (10) into Equation (11), the damping ratio $\zeta$ can be solved by (12)

$$\zeta ={\displaystyle \frac{1}{{\omega }_n(t_2-t_1)}}\ln \left[{\displaystyle \frac{a(t_1)}{a(t_2)}}\right].$$

(12)

Due to the interference from changes in external conditions and the influence of high-order vibration modes of the robotic arm during operation, the signals measured by the acceleration sensor for vibration measurement may have certain accidental errors, and the end vibration parameters also vary under different poses. As a feedforward control method, vibration suppression based on input shaping requires obtaining the vibration parameters of the system in advance. The most direct way to obtain the vibration parameters is to run the required trajectory of the robotic arm in advance and calculate the vibration parameters at the end point of the trajectory based on the vibration signal when it reaches there. However, in practical applications, it is impossible to run the robotic arm in advance before each movement, so it is necessary to estimate the vibration parameters.

2.2. Parameter Estimation of Residual Vibration at the End of a Robotic Arm in Arbitrary Poses

When controlling the robotic arm to complete a task, the end point of each segment of its motion trajectory is a specified input, that is, the joint angles of the six joints are known; the load carried by the end of the robotic arm can also be measured before each movement. All input variables are integrated into a column vector $\tilde{x}$, as shown in Equation (13)

$$\tilde{x}=\left[\begin{array}{c}\mathrm{\Theta }\\ m\end{array}\right],$$

(13)

where$\mathrm{\Theta }$represents the joint angles of the robotic arm, which is a six-dimensional column vector; $m$ is the weight of the load carried by the end of the robotic arm; the vibration parameters to be estimated are the natural frequency ${\omega }_n$and the damping ratio $\zeta$, which are denoted as a column vector $y={[w_n,\zeta ]}^T$.

If the robotic manipulator at each moment is equivalent to an ideal flexible body with uniform texture, according to the Euler-Bernoulli beam theory, its natural frequency has a quadratic relationship with the length of the end from the fixed axis, while the end distance has a trigonometric function relationship with each joint angle. A quadratic polynomial model is established between the input composed of joint angles and end loads, and the output composed of end vibration frequency and damping ratio, as shown in (14). For the sake of concise description, the load mass is denoted as ${\theta }_7=m$,

$$\left\{\begin{array}{l}{\omega }_n={\alpha }_0+{\displaystyle \sum _{j=1}^7{\alpha }_j}{\theta }_j+{\displaystyle \sum _{k=1}^7{\displaystyle \sum _{l=k}^7{\alpha }_{kl}}}{\theta }_k{\theta }_l+{ϵ}_1\hfill \\ \zeta ={\beta }_0+{\displaystyle \sum _{j=1}^7{\beta }_j}{\theta }_j+{\displaystyle \sum _{k=1}^7{\displaystyle \sum _{l=k}^7{\beta }_{kl}}}{\theta }_k{\theta }_l+{ϵ}_2\hfill \end{array}\right.$$

(14)

where${\alpha }_j$and${\beta }_j$denote the coefficients of the linear term ${\theta }_j$;${\alpha }_{kl}$ and ${\beta }_{kl}$denote the coefficients of the quadratic term and coupling term${\theta }_k{\theta }_l$; ${\theta }_k{\theta }_l$and ${ϵ}_2$denote the deviations between the approximate model and the actual values. Both the system's natural frequency ${\omega }_n$ and damping ratio $\zeta$ contain 1 constant term, 7 linear terms ${\theta }_i$, 7 quadratic terms ${\theta }_i^2$, and 21 cross terms ${\theta }_i{\theta }_j(i\ne j)$, which can be written in matrix form as (15):

$$y=\left[\begin{array}{ccccccc}{\alpha }_0& {\alpha }_1& \dots & {\alpha }_7& {\alpha }_{11}& \dots & {\alpha }_{77}\\ {\beta }_0& {\beta }_1& \dots & {\beta }_7& {\beta }_{11}& \dots & {\beta }_{77}\end{array}\right]\left[\begin{array}{c}1\\ {\theta }_1\\ ⋮\\ {\theta }_7\\ {\theta }_1^2\\ {\theta }_1{\theta }_2\\ ⋮\\ {\theta }_6{\theta }_7\\ {\theta }_7^2\end{array}\right]=Bx.$$

(15)

In the equation (15),$B$represents the coefficient matrix, which is a $2\times 36$ matrix;$x$represents the input variable, which is a 36-dimensional column vector containing a constant term, linear terms, quadratic terms, and cross terms.

In this paper, the least square method is used to identify the parameters of the coefficient matrix $B$. Assuming that there are$N$observed data points obtained through experimental measurement, a target function $J(B)$is first constructed as the sum of squared errors between the predicted output${\widehat{y}}_i=B{x}_i$and the actual output.

$$J(B)={\displaystyle \sum _{i=1}^N{‖y-{\widehat{y}}^i‖}^2}={\displaystyle \sum _{i=1}^N{‖y-B{x}^i‖}^2}.$$

(16)

For convenience of processing, all input variables in the measured data are stacked into an overall input matrix $X$, which is an $N\times 36$ matrix,

$$X=\left[\begin{array}{c}{}^{1}x^T\\ {}^{2}x^T\\ ⋮\\ {}^{N}x^T\end{array}\right].$$

(17)

All output variables are stacked into an overall output matrix $Y$, which is an $N\times 36$ matrix shown in (18),

$$Y=\left[\begin{array}{c}{y}^{1T}\\ y^{2T}\\ ⋮\\ y^{NT}\end{array}\right]=\left[\begin{array}{cc}{}^1\omega _n& {}_1\zeta \\ \begin{array}{c}{}^2\omega _n\\ \begin{array}{c}⋮\\ {}^N\omega _n\end{array}\end{array}& \begin{array}{c}\begin{array}{c}{}_2\zeta \\ ⋮\end{array}\\ {}_N\zeta \end{array}\end{array}\right].$$

(18)

The objective function $J(B)$can be re-expressed, and the least squares problem is obtained,

$$\min J(B)={‖Y-Bx‖}_F^2.$$

(19)

The symbol $\Vert \cdot {\Vert }_F$represents the Frobenius norm, which is the square root of the sum of the squares of all its elements; expanding the objective function $J(B)$ gives,

$$\begin{array}{ll}J(B)& =\mathrm{tr}\left[{(Y-Bx)}^T(Y-Bx)\right]\\ & =\mathrm{tr}\left(Y^{T}Y-Y^{T}XB-B^T{X}^{T}Y+B^T{X}^{T}XB\right)\end{array}.$$

(20)

According to the principle of least squares, to minimize the objective function $J(B)$, it is necessary to take the derivative of matrix$B$and set the derivative to zero as shown in (21),

$${\displaystyle \frac{\partial J(B)}{\partial B}}=-2X^{T}Y+2X^{T}XB=0.$$

(21)

Finally, the estimated coefficient matrix $\widehat{B}$of the manipulator in any poses is obtained by (22),

$$\widehat{B}={\left(X^{T}X\right)}^{-1}{X}^{T}Y.$$

(22)

3. Active Suppression of Residual Vibration of Manipulators Based on Input Shapers

3.1. ZV Input Shaper

Input shaping is an open-loop control method that shapes the input signal through a pulse sequence diagram to achieve rapid suppression of residual vibrations in the system. The frequency domain expression of the input shaper is shown in Equation (23),

$$S(s)={\displaystyle \sum _{i=1}^n{A}_i}{e}^{-t_{i}s}.$$

(23)

In the Equation (23), $A_i$represents the amplitude of the pulse sequence; $t_i$represents the time delay of the pulse sequence; $n$represents the number of pulses contained in the input shaper. It can be seen from the mathematical expression of the input shaper that designing the input shaper is to determine the number of its pulses as well as the time delay and amplitude of each pulse. In this paper, three methods, namely ZV, ZVD and EI, are respectively used to design the input shaper, and their advantages and disadvantages in the application of heavy-load cooperative manipulators are analyzed. The ZV shaper has the following form,

$$S(s)=A_1{e}^{-t_{1}s}+A_2{e}^{-t_{2}s}.$$

(24)

As can be seen from Equation (8), to completely eliminate residual vibration, the residual vibration percentage $A_{\mathrm{amp}}$ can be set to 0. By observing the composition of $A_{\mathrm{amp}}$, it can be known that $A_{\mathrm{amp}}=0$ means making the two squared terms under the square root zero respectively, that is, as shown in Equation (25),

$$\left\{\begin{array}{l}{\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\cos ({\omega }_d{t}_i)=0\hfill \\ {\displaystyle \sum _{i=1}^n{A}_i}{e}^{\zeta {\omega }_n{t}_i}\sin ({\omega }_d{t}_i)=0\hfill \end{array}\right..$$

(25)

Since the ZV shaper contains only two pulses (i.e.,$i=2$), the left-hand side of each equation in Equation (25) is the sum of two terms,

$$\left\{\begin{array}{l}{A}_1{e}^{\zeta {\omega }_n{t}_1}\cos ({\omega }_d{t}_1)+A_2{e}^{\zeta {\omega }_n{t}_2}\cos ({\omega }_d{t}_2)=0\hfill \\ A_1{e}^{\zeta {\omega }_n{t}_1}\sin ({\omega }_d{t}_1)+A_2{e}^{\zeta {\omega }_n{t}_2}\sin ({\omega }_d{t}_2)=0\hfill \end{array}\right..$$

(26)

The Equation (26) obviously has infinitely many solutions. To obtain a unique solution, there needs some constraint conditions. To minimize the time required for vibration elimination, the shaper can be made to start acting from time 0, i.e., $t_1=0$. This constraint is also called the time-optimal constraint. To ensure that the target position of the system remains unchanged, the amplitudes of the pulses of the shaper are required to satisfy the constraint as shown in Equation (27)

$$A_1+A_2=1$$

(27)

This constraint is also referred to as the amplitude constraint. By combining Equation (26) and Equation (27), and substituting $t_1=0$ into them, the ZV shaper can be solved as Equation (28),

$$\left\{\begin{array}{l}{A}_1={\displaystyle \frac{e^{\zeta {\omega }_n{t}_2}}{1+e^{\zeta {\omega }_n{t}_2}}}\hfill \\ A_2={\displaystyle \frac{1}{1+e^{\zeta {\omega }_n{t}_2}}}\hfill \\ t_1=0\hfill \\ t_2={\displaystyle \frac{\pi }{{\omega }_n\sqrt{1-{\zeta }^2}}}\hfill \end{array}\right..$$

(28)

3.2. ZVD Input Shaper

When the model parameters of the system are accurately known, the ZV shaper can completely eliminate residual vibrations within half a vibration period. However, when the model parameters are not precisely known, the ZV shaper cannot suppress residual vibrations within the ideal time because it fails to effectively cancel the zeros of the system. To improve the vibration suppression effect of the input shaper under the condition of uncertain natural frequency parameters of the system, a constraint condition on the residual vibration percentage $A_{\mathrm{amp}}$ is introduced, and a ZVD input shaper with multiple zeros configured for the shaper is proposed. Treating the natural frequency ${\omega }_n$ as an independent variable and the residual vibration percentage $A_{\mathrm{amp}}$ as a function of the natural frequency, denoted as $A_{amp}({\omega }_n)$, to configure multiple zeros for the shaper means that the derivative of the residual vibration percentage $A_{amp}({\omega }_n)$ at the actual natural frequency ${\omega }_n={\omega }_{nr}$ should be zero. Then, the newly introduced constraint condition is as shown in Equation (29),

$${\left.{\displaystyle \frac{\mathrm{d}{A}_{amp}({\omega }_n)}{\mathrm{d}{\omega }_n}}\right|}_{{\omega }_n={\omega }_{nr}}=0$$

(29)

where${\omega }_{nr}$represents the actual undamped natural frequency. The ZVD shaper contains three pulses, i.e., $i=3$, and the specific amplitude constraint is

$$A_1+A_2+A_3=1.$$

(30)

By combining Equation (29) and (30), and substituting the time constraint $t_1=0$ into them, the ZVD shaper can be solved as Equation (31) ,

$$\left\{\begin{array}{l}{A}_1={\displaystyle \frac{1}{1+2K+K^2}}\hfill \\ A_2={\displaystyle \frac{2K}{1+2K+K^2}}\hfill \\ A_3={\displaystyle \frac{K^2}{1+2K+K^2}}\hfill \\ t_1=0\hfill \\ t_2=T\hfill \\ t_3=2T\hfill \end{array}\right.$$

(31)

where $T={\displaystyle \frac{\pi }{{\omega }_n\sqrt{1-{\zeta }^2}}};K=e^{{\displaystyle \frac{-\zeta \pi }{\sqrt{1-{\zeta }^2}}}}$.

3.3. EI Input Shaper

When the heavy-load cooperative manipulator is working, the parameters identified in advance will change, making it impossible for the model to be completely accurate at all times. By using the EI input shaper, an allowable residual vibration amount $A_{amp}^{tol}$ is specified, and it is only required that the residual vibration of the system at the natural frequency is less than or equal to $A_{amp}^{tol}$, as shown in Equation (32),

$$\left\{\begin{array}{l}{A}_{amp}({\omega }_l)=0\hfill \\ A_{amp}({\omega }_h)=0\hfill \\ A_{amp}({\omega }_n)=A_{amp}^{tol}\hfill \\ {\left.{\displaystyle \frac{\mathrm{d}{A}_{amp}({\omega }_n)}{\mathrm{d}{\omega }_n}}\right|}_{{\omega }_n={\omega }_{nr}}=0\hfill \end{array}\right..$$

(32)

By combining Equation (32), the amplitude constraint (30) and $t_1=0$, the parameters of each pulse of the EI shaper can be solved as following [32],

$$\left\{\begin{array}{l}{A}_1={\displaystyle \frac{1+V}{4}}\hfill \\ A_2={\displaystyle \frac{1-V}{2}}\hfill \\ A_3={\displaystyle \frac{1+V}{4}}\hfill \\ t_1=0\hfill \\ t_2=T\hfill \\ t_3=2T\hfill \end{array}\right.$$

(33)

where $T={\displaystyle \frac{\pi }{{\omega }_n}}$. The effect of the EI input shaping also varies when $A_{amp}^{tol}$takes different values.

4. Experiments and Analysis

4.1. Experiment on Vibration Suppression of a Single-Joint Manipulator

In this section, a single-joint manipulator is used to verify the effect of input shapers on suppressing residual vibrations, and the vibration suppression effect, robustness, and time required for vibration suppression of various input shapers are compared. The experimental equipment is shown in Figure 2. The single-axis manipulator is driven by a joint motor to actuate a steel ruler with a certain degree of flexibility. An acceleration sensor is installed at the end of the steel ruler, through which the vibration signal can be measured. In this experiment, the joint motor will be controlled to rotate by a certain angle, and after reaching the end position, the steel ruler will have residual vibrations due to inertia. The joint motor used in this paper is the one used in the CR series manipulators of Dobot robotics, which can be simulated and controlled by a Speedgoat controller.

A step signal with an amplitude of 45°is given to the robot joint. The vibration signal measured by the accelerometer installed at the end of the steel ruler is shown in Figure 3 (a). Using Equations (28) ,(31) and (33), the frequency of the end vibration is calculated to be 17.872 Hz, and the damping ratio is 0.0089. Based on the derivation in Section 3.1, ZV, ZVD, and EI shapers can be designed as shown in the equations:

$$\left\{\begin{array}{l}{S}_{ZV}(s)=0.5070+0.4930e^{-0.1758s}\\ S_{ZVD}(s)=0.2570+0.4999e^{-0.1758s}+0.2431e^{-0.3516s}\\ S_{EI}(s)=0.2625+0.4750e^{-0.1758s}+0.2625e^{-0.3516s}\end{array}\right.$$

(34)

The acceleration data of the single manipulator under the action of the three input shapers are shown in Figures 3 (b), (c), and (d) respectively. By comparison in Table I, the maximum amplitude of the end without input shaping is 1.263$m/s^2$, while the maximum amplitudes under the action of the ZV, ZVD, and EI shapers are 0.51, 0.393, and 0.375 respectively, which are reduced by 59.7%, 68.9%, and 70.3% accordingly. It can be seen that the ZV, ZVD, and EI shapers all have inhibitory effects on residual vibrations. Among them, the ZVD shaper and EI shaper have slightly better vibration suppression effects than the ZV shaper, and the EI shaper performs slightly better than the ZVD shaper. In terms of response time and vibration suppression time, the ZV shaper takes 0.255 s to suppress the residual vibration to below 10%, the ZVD shaper takes 0.37 s, and the EI shaper takes 0.395 s. It can be observed that after adding the shapers, the response time of the system increases significantly, as the shapers introduce a certain time delay into the system. In comparison, the ZV shaper can suppress the residual vibration to a relatively low level in a shorter time, while the ZVD and EI shapers require a longer time. The experimental videos of the single-axis robot in four forms are shown in the appendix.

Table 1. Vibration suppression effect of single-joint manipulator.

Items	Maximum Amplitude$(\mathit{m}/{\mathit{s}}^2$)	Optimization Effect
NO Input Shaper	1.263	/
ZV Input Shaper	0.510	59.7%
ZVD Input Shaper	0.393	68.9%
EI Input Shaper	0.375	70.3%

Table 1. Vibration suppression effect of single-joint manipulator.

Items	Maximum Amplitude$(\mathit{m}/{\mathit{s}}^2$)	Optimization Effect
NO Input Shaper	1.263	/
ZV Input Shaper	0.510	59.7%
ZVD Input Shaper	0.393	68.9%
EI Input Shaper	0.375	70.3%

4.2. Experiment on Heavy-Load Six-Axis Collaborative Robot

4.2.1. Introduction to the Experimental Platform

Figure 4 shows the heavy-load (16kg load) six-axis collaborative manipulator, and its basic parameters are shown in Table 2. In the manipulator control program, after initialization and related calculations, the program is divided into two threads, the manipulator motion control thread and the sensor data receiving thread. In terms of sensors, first of all, since the maximum residual vibration frequency at the end of some measured points of the Dobot robotics CR16 manipulator is approximately 75 Hz, in order to obtain a relatively accurate vibration waveform, the sampling frequency should preferably be no less than 40 times the system vibration frequency. That is, the selected accelerometer sampling rate and data return frequency need to be no less than 3000 Hz. Secondly, since the roughly measured maximum vibration acceleration of the Dobot robotics CR16 manipulator when reaching the end point at 70% of the maximum operating speed is about 1g, the range of the selected sensor should be no less than ±2g. In addition, to make the experimental results more accurate, sensors with small resolution and low noise should be selected as much as possible.

4.2.2. Vibration Parameter Identification for Arbitrary Poses

Since the reading of the accelerometer can be performed in a different thread of the same program as the manipulator control, the accelerometer data and the robot joint encoder signals can be collected simultaneously while the manipulator is in motion. When the manipulator moves with a load, it can automatically run through all the points that need to be measured, with sufficient time intervals reserved between each point to allow the residual vibration at the end to stop completely. Since the vibration parameters of the manipulator are significantly more affected by axes 2, 3, and 4 than by other axes, the joint angle intervals of these three axes are subdivided into smaller increments. After acquiring a set of operation data, the data segments corresponding to all points are identified, and the vibration frequency and damping ratio of each are calculated using the method described in Section 2.1. For example, point A in Figure 5 is extracted separately. In the experiment to obtain the estimated model, due to symmetry, the movement of axes 1 and 6 obviously does not affect the vibration characteristics of the manipulator, so they are not subdivided. The 90°movement range of axis 2 is subdivided into intervals of 15°, resulting in a total of 7 joint angle positions. The 180°movement ranges of axes 2, 3, and 4 are also subdivided into 15°intervals; after excluding angles that are inaccessible due to structural constraints, there are 10 joint angle positions for each axis. In addition, before conducting the experiment, the Cartesian coordinates of the end-effector are solved using the forward kinematics of the manipulator, and its movement is restricted to a safe range. When the load is fixed, the manipulator runs once to traverse all end points and collect data, resulting in a total of 2328 sets of data.

Some points are selected to calculate their vibration parameters using the estimation model, and the errors are computed by comparing these parameters with those measured experimentally at the same points. Table 3 shows the observation results of some points, where ${\widehat{\omega }}_n$ is the estimated value of the system's natural frequency; ${\omega }_n$is the experimentally measured value of the system's natural frequency; $(\Delta {\omega }_n)\%$is the estimation error of the system's natural frequency; $\widehat{\zeta }$is the estimated value of the system's damping ratio; $\zeta$is the experimentally measured value of the system's damping ratio; and $(\Delta \zeta )\%$ is the estimation error of the system's damping ratio.

4.2.3. Active Suppression of Manipulator in arbitrary poses

For input shaping technology, it is only necessary to shape the control signal before it is input to the actuator. In this paper, the built-in position controller of the Dobot CR16 manipulator is used. In the joint space, the position signal is input into the shaper and convolved with the shaper to form a new control signal. The design flow in the joint space is shown in Figure 6.

Figure 7 and Table 4 show that the three input shaping methods (ZV, ZVD, and EI) can all exert a certain inhibitory effect on the residual vibration at the end of the six-axis collaborative manipulator. Without input shaping, the maximum amplitude at the end-effector is 0.4692 g. With the ZV, ZVD, and EI shapers applied, the maximum amplitudes are 0.3155 g, 0.1849 g, and 0.1825 g respectively, representing reductions of 32.8%, 60.6%, and 61.1% in amplitude. It can be seen that the ZV, ZVD, and EI shapers all have inhibitory effects on residual vibration. Among them, the ZVD and EI shapers exhibit better vibration suppression performance than the ZV shaper, with the EI shaper performing slightly better than the ZVD shaper.In terms of response time and vibration suppression time, the ZV shaper takes 0.13s to suppress residual vibration to below 10%, the ZVD shaper takes 0.19s, and the EI shaper takes 0.17s. This indicates that after adding a shaper, the system's response time increases significantly, meaning the shaper introduces a certain time delay into the system. In comparison, the ZV shaper can suppress residual vibration to a relatively low level in a shorter time, while the ZVD and EI shapers require a longer time.

Table 5 presents the vibration suppression effects of the ZV, ZVD, and EI shapers designed respectively using the parameters estimated in Table 3. The experimental results show that all three input shaping methods (ZV, ZVD, and EI) can exert a certain inhibitory effect on the residual vibration at the end-effector of the six-axis collaborative manipulator. When the parameter estimation error is relatively large, the vibration suppression effect of the ZV shaper decreases significantly, while the ZVD and EI shapers perform slightly better. A further comparison with Table 3 reveals that the more accurate the parameter estimation, the better the effect of the input shaper on suppressing residual vibration. When the parameter identification error is small, the average suppression effect is over 40%, and especially when the parameter identification error is small, the suppression effect reaches more than 70%.

5. Conclusions

To address the issue of residual vibration at the end-effector of heavy-load six-axis collaborative robots in arbitrary poses, this paper conducts research on the modeling of the manipulator end vibration model, the estimation of manipulator vibration parameters in arbitrary poses, and carries out active vibration suppression experiments on residual vibration at the end of heavy-load six-axis collaborative robots in arbitrary poses based on ZV, ZVD, and EI shapers. The following achievements have been made: (1) The mathematical forms of input shaping methods are derived, including the design and robustness analysis of several common input shapers. A single-axis experimental platform is designed, and simulation and physical experiments are conducted to verify the effect of input shaping methods on suppressing residual vibration, while comparing the effects and characteristics of several input shapers. (2) Considering that the end vibration parameters of the robot vary under different poses and that it is impossible to run the robot in advance before each movement, this paper establishes a model between the robot joint angles, end load mass, and end vibration parameters, and adopts the least squares identification method for parameter identification. (3) Based on the established vibration parameter model, the vibration parameters of the robot end in arbitrary poses are estimated. Input shapers are designed according to the derived input shaper design method, and then the six-axis collaborative robot is controlled based on its kinematic model, achieving a certain suppression effect on the residual vibration at the robot end. The average suppression effect is over 40%, and especially when the parameter identification error is small, the suppression effect reaches more than 70%.