2. Generation of the local variable period in the stamping line

As the typical application scenario of industrial robots, the high-speed automobile-panel stamping line alternately cooperates with multiple feeding mechanisms and stamping stations for completing material transmission and stamping forming [Reference Yu, Wang, Zhang, Zhang and Zhang23]. Specifically, when the upper and lower dies of the stamping station i-1 are separated to a certain height, the manipulator i starts to move from the initial pose in Fig. 1 to a certain height above the lower die i-1 to pick up the stamping part. Then, the manipulator i transfers the stamping part to place the stamping part in the lower die i. Finally, the manipulator i moves to the initial pose in Fig. 1 to start the next cycle of work. Notably, when the manipulator leaves the die area of the stamping station i, stamping forming of the stamping part will be executed.

Figure 1. Initial poses of multiple machines in the stamping line.

The trajectory performance of the stamping line is affected by trajectories of feeding mechanisms and stroke curves of presses in the stamping station. The feeding mechanism in Fig. 1 is the redundant manipulator. The stamping station includes the upper die, lower die, and press. The press is located inside the upper die. In the single-machine trajectory planning, motion periods of press and manipulator are generally different. The manipulator with periodic material transmission may collide with up to six neighbouring models. The manipulator i in Fig. 1 is taken as an example. It may simultaneously collide with manipulators i-1 and i + 1, upper dies i-1 and i, lower dies i-1 and i in one motion period. To ensure the continuous running of multiple machines in Fig. 1, the motion period of the manipulator should be adjusted to be the same as that of the press.

The manipulator i in Fig. 1 picks up stamping part from the lower die i-1 and places stamping part in the lower die i. The district of the stamping station in the horizontal direction is called the die area. Based on Eq. (1), the period P _M of the manipulator before multi-machine coordination is discretized into $n_{M}$ time points.

(1) \begin{equation} t_{M}=\frac{P_{M}}{n_{M}-1} \end{equation}

where $t_{M}$ is the fixed-time interval of the manipulator before multi-machine coordination. $n_{M}$ is the number of discrete time points.

When the manipulator i in Fig. 1 enters and exits the die area of the stamping station i-1, corresponding labelled numbers of discrete time points are assumed to be m ₁ and m ₂, respectively. Similarly, when the manipulator i enters and exits the die area of the stamping station i, corresponding labelled numbers of discrete time points are assumed to be m ₃ and m ₄, respectively. If the manipulator runs for one period, the total number of discrete fixed time points in the die area will be assumed to $n_{M\_ in}$ . From the perspective of labelled numbers, $n_{M\_ in}$ can be obtained by Eq. (2).

(2) \begin{equation} n_{M\_ in}=m_{2}-m_{1}+m_{4}-m_{3} \end{equation}

Decentralized trajectory planning is widely used in the high-speed automobile-panel stamping line. Trajectory planning of the single redundant manipulator is firstly executed. Then, the period of the stamping line is generated. The period of the single redundant manipulator is not less than that of the stamping line. If multi-machine coordination is executed, the period of the single redundant manipulator will increase to the same as that of the stamping line. The period difference between the period of the stamping line and that of the manipulator is defined as $P_{d}$ . It can be determined by Eq. (3).

(3) \begin{equation} P_{d}=P_{S}-t_{M}n_{M\_ in} \end{equation}

where $P_{S}$ is the period of stamping line after multi-machine coordination.

Before the multi-machine coordination, the period of the manipulator is greater than that of press. To ensure high-speed and continuous running of the stamping production line, the time interval of the manipulator in the die area is equal to its initial fixed time interval. The time interval of the manipulator out of the die area should be greater than that in the die area.

3. Smooth interpolation method of variable-time interval

Improper distribution of time intervals out of die area may generate poor-performance trajectories in Cartesian and joint space. To maintain the original trajectory performance of the stamping line after the variable-period operation, detailed steps of the proposed method in Fig. 2 are introduced in this section.

Figure 2. Steps for the smooth interpolation method of variable-time interval.

Step 1: Input original trajectories of end effectors and redundant joints

Assume that the redundancy of the N-DOF manipulator is R. Joints of the manipulator are divided into redundant and non-redundant joints. Trajectories of end effectors and redundant joints are preferentially input.

Step 2: Discrete initial trajectory based on fixed-time intervals

Based on Eq. (1), the input trajectory is discretized. The recommended $n_{M}$ is approximately 4000. If $n_{M}$ is too small, one fixed-time interval will occupy much time to deteriorate the optimal trajectory performance, and invalid solutions of non-redundant joints will be generated. If $n_{M}$ is too large, the time consumption of trajectory optimization will be increased.

Step 3: Generate optimal trajectory under piecewise variable-time intervals

In the die area, time intervals of input trajectories remain constant. Combined with $P_{d}$ in Eq. (3), time intervals of input trajectories out of the die area should be adjusted. To ensure the smooth transition of the trajectory, the changing trend of variable-time increments is slow, fast, and slow. Under this circumstance, normal distribution in Eq. (4) is introduced. It highly matches the changing requirement of variable-time intervals and is easy to operate in MATLAB. However, there are many unknown parameters in Eq. (4), which increases the uncertainty of time allocation.

(4) \begin{equation} f\left(t\right)=\frac{1}{\sqrt{2\pi }\sigma }e^{-\frac{\left(t-\mu \right)^{2}}{2\sigma ^{2}}},-\infty \lt t\lt +\infty \end{equation}

Where $\mu$ is average, $\sigma$ is standard deviation, t is the random variable, and f(t) is the probability density function.

If $\mu$ and $\sigma$ are constant, Eq. (4) will be transformed into t∼N( $\mu, \sigma ^{2}$ ). Combined with the normal distribution function named normpdf in MATLAB, the curve of the normal distribution can be plotted in Fig. 3.

Figure 3. Normal distribution.

As shown in Fig. 3, assume that t in Eq. (4) ranges from -t _i to t _j . Within this rang, $P_{d}$ in Eq. (3) is equal to the integral area of the curve with respect to the horizontal coordinate. To maintain the consistent number of trajectory points, the number of time points between -t _i and t _j should satisfy the Eq. (5).

(5) \begin{equation} n_{t}=n_{M}-n_{M\_ in} \end{equation}

where n _t is the number of time points.

There are two variable-period districts for manipulator i in Fig. 1. One district is the manipulator transferring the stamping part from stamping stations i-1 to i out of the die area. The other district is the manipulator running from stamping stations i to i-1 out of the die area. In this section, the former district is taken as an example to illustrate the generation of optimal trajectory under piecewise variable-time intervals. For the manipulator transferring stamping part from stamping stations i-1 to i out of the die area, labelled numbers of discrete time points range from m ₂ to m ₃. Meanwhile, corresponding piecewise variable-time intervals $t'_{\!\!I}$ can be obtained by Eq. (6). Notably, $t'_{\!\!I}$ is the increment of variable-time interval relative to fixed-time interval.

(6) \begin{equation} t'_{\!\!I}=t'_{\!\!I-1}+t_{M}+P'_{\!\!d}\frac{T_{I}}{\sum _{I^{\prime}=1}^{m_{3}-m_{2}}T_{I^{\prime}}} \end{equation}

where I ranging from 1 to (m ₃-m ₂) is the sequence number of discrete time points. $t'_{\!\!0}$ is equal to the time increment whose labelled number of discrete time point is m ₂. T _I in Eq. (7) is the single-step increment ratio of $P'_{\!\!d}$ . $P'_{\!\!d}$ in Eq. (8) is the time difference between the time of the optimal trajectory and that of the initial trajectory when the manipulator transferring stamping part from stamping stations i-1 to i out of the die area.

(7) \begin{equation} T_{I}=\text{normpdf}\left(I\frac{\left(t_{i}+t_{j}\right)t_{M}n_{t}}{P_{d}\left(m_{3}-m_{2}\right)}-t_{i}, \mu,\sigma \right) \end{equation}

(8) \begin{equation} P'_{\!\!d}=\frac{m_{3}-m_{2}}{n_{M}-n_{M\_ in}}P_{d}-\left(m_{3}-m_{2}\right)t_{M} \end{equation}

To ensure the smooth and symmetrical transition of the variable-time interval, t _i and t _j are assumed to be identical. The generated $t'_{\!\!I}$ in Eq. (6) replacing the original fixed time point out of the die area becomes the new piecewise variable time point. But the trajectory value corresponding to the time point remains unchanged.

Step 4: Build equations of inverse kinematics

Based on the Denavit-Hartenberg or geometric method, equations of forward kinematics can be obtained. Equations of inverse kinematics from Cartesian space to joint space can be further built.

Step 5: Generate trajectories of non-redundant joints

In one motion period, time intervals of non-redundant joints correspond to that of end effectors and redundant joints under piecewise variable-time intervals. Combined with inverse kinematic equations in Step 4, the trajectory of the non-redundant joint can be further ensured.

Step 6: Compute velocities, accelerations, and jerks of optimal trajectories

Assume that there is one knot vector U={u ₀, u ₁, …, u _m } and n + 1 control points P ₀, P ₁, …, P _n . The B-spline curve p(u) of degree k defined by knot vector and control points is shown in Eq. (9).

(9) \begin{equation} p\left(u\right)=\sum _{i=0}^{n}N_{i,k\left(u\right)}P_{i} \end{equation}

where, N _i,k(u) is the B-spline basis function of degree k. Its mathematical expression is shown in Eq. (10). Notably, m=n+k + 1.

(10) \begin{equation} \left\{\begin{array}{l} N_{i,0\left(u\right)}=\left\{\begin{array}{l} 1,\;if\; u_{i}\leq u\leq u_{i+1}\\[5pt] 0,\;\textit{otherwise} \end{array}\right.\\[16pt] N_{i,k\left(u\right)}=\dfrac{u-u_{i}}{u_{i+k}-u_{i}}N_{i,k-1}\left(u\right)+\dfrac{u_{i+k+1}-u}{u_{i+k+1}-u_{i+1}}N_{i+1,k-1}\left(u\right) \end{array}\right. \end{equation}

A series of discrete trajectory points can be fitted by the quintic B-spline curve in Eq. (9). The programme of the proposed method is executed in MATLAB. Table I lists several typical B-spline postprocessing functions in MATLAB.

Table I. Typical B-spline postprocessing functions in MATLAB.

Combined with Eq. (9) and functions in Table I, displacements, velocities, accelerations, and jerks of optimal trajectories under piecewise variable-time intervals can be obtained.

Step 7: Output optimal trajectories in Cartesian and joint space

A sharp point in the motion curve may bring an abrupt change of force and acceleration. The transient vibration and overshoot of the manipulator may further occur. The lower the transient vibration and overshoot of the manipulator are, the higher performance the motion trajectory will be ref. [Reference Gao, Li, He and Li24]. Additionally, the average acceleration is related to the torque reflecting the energy consumption of motor [Reference Lan, Xie, Liu and Cao10]. From the perspective of the trajectory planning, to maintain the high performance of motion trajectory, following optimal standards should be satisfied.

1. I. There is no sharp point on the velocity curve during one motion period.

2. II. There is no collision among machines during continuous motion periods.

3. III. The multi-objective optimization function in Eq. (11) should be satisfied.

   (11) \begin{equation} \begin{array}{l} \min f_{1}=\sum \sqrt{\frac{1}{P_{M}}\int _{0}^{P_{M}}a_{i}^{2}dt}\\[5pt] \min f_{2}=\sum \sqrt{\frac{1}{P_{M}}\int _{0}^{P_{M}}j_{i}^{2}dt}\\[5pt] \!s.t. \left| a_{i}\right| \leq a_{\max }\\[5pt] \quad\; \left| \,j_{i}\right| \leq j_{\max } \end{array} \end{equation}

where f ₁ and f ₂ are functions to measure the energy consumption and evaluate the smoothness of the trajectory, respectively [Reference Lan, Xie, Liu and Cao10]. a _max and j _max are maximum values of acceleration and jerk, respectively. P _M is the period of the manipulator.

Before the local variable-period trajectory optimization, the initial fixed-period trajectories meet the application requirement. Compared with f ₁ and f ₂ of the initial fixed-period trajectory, if f ₁ and f ₂ of the local variable-period trajectory are not significantly increased, the optimal standard III will be satisfied. If above optimal standards are satisfied, optimal trajectories in Cartesian and joint space can be output. Otherwise, Steps 3 to 7 will be re-executed. The proposed method will not be terminated until the optimal trajectory satisfying the above standards is generated.

4. Case study

The redundant planar manipulator with 5 DOFs in Fig. 4 is taken as the case to illustrate the proposed method. The redundancy of the manipulator is 2. In the previous work [Reference Yu, Wang, Zhang and Zhang25], trajectories in Cartesian and joint space can be obtained by the method based on the joint classification and particle swarm optimization algorithm. Joints D ₄ and R ₆ are redundant joints. They are beneficial to avoid the generation of singularities during the trajectory planning of the single redundant planar manipulator. The period of the trajectory is 3.75 s. It is discretized into 4097 time points. Corresponding labelled numbers m ₁, m ₂, m ₃, and m ₄ are 531, 1846, 2306, and 3552. Based on Eq. (2), $n_{M\_ in}$ is 2561. The period $P_{S}$ of the stamping line is assumed to be 5.6881 s. Based on Eq. (3), $P_{d}$ is approximately 3.3434 s. Additionally, some definitions and values of parameters in Fig. 4 are listed in the previous published paper [Reference Yu, Wang, Zhang and Zhang25].

Figure 4. The redundant planar manipulator.

The trajectory in the vertical direction in Cartesian space is prioritized for variable-period operation. To ensure the high performance of the trajectory, t _i and t _j in Fig. 3 are both 8. Two stamping stations i-1 and i in Fig. 1 divide the motion trajectory of manipulator i into four sections. Two sections are out of the die area. Based on m ₂ and m ₃, the number of time intervals for manipulator i running from stamping stations i-1 to i out of the die area is 460. Based on Eq. (8), the corresponding $P'_{\!\!d}$ is 0.5801 s. $\mu$ and $\sigma$ in Eq. (4) are 0 and 3.2, respectively. Similarly, based on $n_{M}$ , m ₁ and m ₄, the number of time intervals and corresponding $P'_{\!\!d}$ for manipulator i running from the stamping stations i to i-1 out of the die area are 1076 and 1.3570 s, respectively. $\mu$ and $\sigma$ in Eq. (4) are 0 and 1, respectively. Combined with Eq. (7) and functions in Table I, the new increment of time intervals out of the die area can be obtained. Combined with Eq. (6), the original trajectory points corresponding to fixed-interval time points are transformed into that corresponding to variable-interval time points. The optimal trajectory of end effector in the vertical direction in Cartesian space can be shown in Fig. 5. Notably, units of displacement s, velocity v, acceleration a, jerk j, and time t are m, m/s, m/s², m/s³, and s, respectively. The red and blue curves in Fig. 5 represent curves of initial and optimal trajectories, respectively.

Figure 5. The end-effector optimal trajectory in the vertical direction. (a) Displacement-time curve, (b) Velocity-time curve, (c) Acceleration-time curve, (d) Jerk-time curve.

In Fig. 5, the manipulator starts at the initial position, runs in one period, and ends in the pose whose labelled number of time points is again m ₁. The whole periods of curves represented by the red and blue curves are 4.2352 and 6.1733 s, respectively. Similarly, as shown in Fig. 6, trajectories of redundant joints and end effector in the horizontal direction under piecewise variable-time intervals can be generated. Notably, units of displacements and velocities for revolute joints are ^o and ^o/s, respectively. The red and blue curves in Fig. 6 represent curves of initial and optimal trajectories, respectively.

Figure 6. Variable-period trajectories of redundant joints and end effector in the horizontal direction. (a) Displacement-time curve of the end effector in the horizontal direction, (b) Velocity-time curve of the end effector in the horizontal direction, (c) Displacement-time curve of joint D _4, (d) Velocity-time curve of joint D _4, (e) Displacement-time curve of joint R _6, (f) Velocity-time curve of joint R ₆.

Based on the geometric method, forward kinematic equations of the redundant planar manipulator in Fig. 4 can be shown in Eq. (12).

(12) \begin{equation} \left\{\begin{array}{l} x=d_{4}\sin r_{3}+l_{5}\sin \left(r_{3}+r_{5}\right)\\[5pt] y=d_{0}+d_{1}+l_{2}\cos r_{2}-\left(d_{4}-l_{3}\right)\cos r_{3}-l_{5}\cos \left(r_{3}+r_{5}\right)\\[5pt] r_{6}=r_{3}+r_{5}\\[5pt] l_{2}\sin r_{2}=l_{3}\sin r_{3} \end{array}\right. \end{equation}

The Eq. (12) is transformed into equations of inverse kinematics in Eq. (13).

(13) \begin{equation} \left\{\begin{array}{l} r_{3}=\arcsin \dfrac{x-l_{5}\sin r_{6}}{d_{4}}\\[9pt] r_{2}=\arcsin \dfrac{l_{3}\sin r_{3}}{l_{2}}\\[5pt] r_{5}=r_{6}-r_{3}\\[5pt] d_{1}=y-d_{0}-l_{2}\cos r_{2}+\left(d_{4}-l_{3}\right)\cos r_{3}+l_{5}\cos \left(r_{3}+r_{5}\right) \end{array}\right. \end{equation}

Under piecewise variable-time intervals, trajectory points of non-redundant active joints D ₁, R ₂, and R ₅ can be obtained by Eq. (13). Combined with Eq. (9) and functions in Table I, their displacements, velocities, accelerations, and jerks can be further generated. Combined with Eq. (11) and standards in Section 3, the high-performance trajectories of non-redundant active joints are shown in Fig. 7. The red and blue curves in Fig. 7 represent curves of initial and optimal trajectories, respectively. f ₁ and f ₂ of trajectories are shown in Table II.

Figure 7. Variable-period trajectories of non-redundant active joints. (a) Displacement-time curve of joint D _1, (b) Velocity-time curve of joint D _1, (c) Displacement-time curve of joint R _2, (d) Velocity-time curve of joint R _2, (e) Displacement-time curve of joint R _5, (f) Velocity-time curve of joint R ₅.

In Table II, X and Y represent end-effector trajectories in the horizontal and vertical directions, respectively. Units of f ₁ and f ₂ for R ₂, R ₅, and R ₆ are ^o/s ² and ^o/s ³, respectively. Similarly, units of f ₁ and f ₂ for other parameters are m/s² and m/s³, respectively. To intuitively compare f ₁ and f ₂ in the initial and optimal trajectories, P _M in Eq. (11) is adjusted to 4.2352 s. f ₁ and f ₂ in the initial and optimal trajectories are approximately identical. They are less than a _max and j _max, respectively. In addition, as shown in Fig. 8, trajectories of active joints are imported to the operation platforms of ADAMS and related actual manipulators, respectively. Simulation and experiment results also show that there is no collision between the manipulator and neighbouring machines. Curves of trajectories and velocities in ADAMS and experiment are consistent with those in MATLAB. Combined with standards in Section 3, Figs. 5–7, and Table II, the generated optimal variable-period trajectories for the redundant planar manipulator in Fig. 4 are available.

Table II. f ₁ and f ₂ corresponding to initial and optimal trajectory.