Autonomous trajectory planner for space telerobots capturing space debris under the teleprogramming framework

Abstract

Due to the special characteristics of the space telerobot system (the dynamic coupling and the dynamic singularities), it is difficult for ground operators to directly operate or plan the orbital telerobot to accomplish the capture task. This article proposes an autonomous trajectory planner under the teleprogramming framework for space telerobots to capture space debris, where the direct optimization method is adopted to obtain smooth and singularity-avoiding joint trajectories. First, the motion planning problem is formulated as a constrained optimization problem, where the precise capture is set as the optimization objective and the base disturbance and the physical limitations (joint position, velocity, and acceleration) are set as constraints. Then, the parameterization and normalization treatments and the penalty function method are adopted to handle constraints and determine the feasible region of the undetermined parameters. Furthermore, an autonomous trajectory planner based on the forward differential kinematics and the differential evolution algorithm is developed to plan the feasible joint trajectories. Finally, the simulation of a space telerobot equipped with a 7-degree-of-freedom manipulator capturing floating debris is conducted to verify the effectiveness of the proposed planner.

Introduction

The population of human-generated space debris is growing in low Earth orbit (LEO) year after year, which poses a serious threat to the current operating satellites. ¹ Therefore, active debris removal (ADR) has been an inevitable mission for the future space environment. As a promising solution, space robots are developed to substitute astronauts to execute some on-orbit servicing (OOS) tasks. ² Unlike industrial robots, most of space robots should be faced with various sophisticated, non-repetitive, and unfamiliar tasks. Fortunately, teleoperation, as one of the most important robotic technologies, extends human’s manipulation, reasoning, and decision-making capabilities to remote robots, which makes these tasks tractable. ³

The direct control mode is the most successful approach in space applications, where the operator remains in the control loop at any time and provides the continuous control signals. For example, in the Canadarm2 and Japanese Experiment Module Remote Manipulator System (JEMRMS) systems, the operators use joysticks to control them to accomplish ordinary tasks in the International Space Station (ISS). ⁴ The first ground-orbit-based telerobot system, Robot Technology Experiment (ROTEX), also adopted the unilateral teleoperation, where a predictive simulator was employed to compensate the communication time delay (5–7 s). ⁵ The Engineering Test Satellite VII (ETS-VII) also carried out the bilateral teleoperation by a master–slave system with a communication time delay of 7 s. ⁶ Besides, based on the particular communication infrastructures, the German Aerospace Center (DLR) accomplished the OOS ground verifications using the bilateral teleoperation, where the round-trip delay is about 620 ms. ⁷ Although the ROTEX and ETS-VII projects are successfully performed in the case of large time delay, to some extent, it largely depends on the fact that operators are well acquainted with the tasks in advance.

In fact, it is not easy to apply the direct control mode in the unfamiliar ADR tasks. There are twofold reasons: (1) the round-trip time delay (more than 1 s) may lead to the result that the feedback signal cannot be received in time, which may affect the stability of the direct control system ⁸ and (2) the special characteristics of the space robots (e.g. the floating base dynamics) lead to the situation that they are difficult to be operated. ⁹ So far, the fully autonomous operation (without human intervention) is still beyond the state of the art of space robotics. For addressing the communication delay, the teleprogramming approach is proposed by Funda et al., ¹⁰ which, as a viable form of the supervisory control, provides a concise framework for controlling the semi-autonomous space robot system. In the teleprogramming framework, the task is first done in the simulation environment; then the operational steps are translated into high-level commands and sent to the remote zone to be executed by semi-autonomous systems under the supervision of operators.

On the other hand, the special characteristics of space telerobot systems and space operational conditions should be taken into account in the capture task. For this satellite-manipulator system, one of the most distinctive features is the dynamic coupling, ¹¹ that is, the motion of the manipulator has a direct physical feedback on its base and the coupling motion of the base will affect the operation performance of the manipulator. Unlike the base-fixed robots,^12,13 the space robots cannot form the accuracy kinematics at the position level. Thus, it is difficult for operators to directly operate space telerobots to capture debris. The motion planning is a necessary prerequisite in the capture task. The early coordinated control method^14,15 employed the jet thrust system to compensate the reaction of the manipulator so that the space robot could behave as a ground-fixed robot. However, the large consumption of physical fuel and the influence (that the compensating jet may blow away the debris) are the major drawbacks in applications. In general, after space robots reached the rendezvous zone, it is preferable to change them into the free-floating mode (the position and attitude systems of bases were shut off) to execute tasks. ² In this working mode, the space robot is an underactuated system. As for the studies of free-floating space robots, Umetani and Yoshida ¹⁶ proposed the generalized Jacobian matrix (GJM) to describe the differential mapping relationship from the joint space to the end effector in the working space. Nakamura and Mukherjee ¹⁷ analyzed the nonholonomic property of free-floating space robots and proposed a bidirectional approach to plan the manipulator and the base simultaneously. Dubowsky and Torres ¹⁸ proposed a technique called enhanced disturbance map (EDM) to plan the manipulator with minimizing the base disturbance. Nenchev ¹⁹ proposed the concept of reaction null-space (RNS) to plan the kinematical redundant manipulator to achieve a zero disturbance on the base. However, the dynamic singularity ²⁰ (i.e. singularities of the GJM) is a specific and severe issue for space robots. Here, the singularities of the GJM are related not only to the geometry parameters but also to the inertia parameters, with the result that the singularity configurations cannot be predicted solely from the joint space. Thus, it is hard for operators to plan a space telerobot like a base-fixed robot.

In view of the above, this article aims at developing an autonomous trajectory planner under the teleprogramming framework for supporting the ground operators. Here, the planning problem is treated as a constrained optimization problem. The objective is set to realize the precise point-to-point capture. Meanwhile, as the directional antenna on the base needs to be kept in a certain range for the data transmission, the attitude disturbance of the base is considered as an inequality constraint. Besides, the physical constraints in joint space are also considered in the motion planning. Based on the parameter optimization technique, the forward differential kinematics, and the differential evolution (DE) algorithm, the autonomous planning procedures are developed to plan joint trajectories.

The rest of this article is organized as follows. The motion equations of the space telerobot are introduced in section “Motion equations.” The teleprogramming control framework and the capture task analysis are discussed in section “Analysis of the capture task.” The design details of the autonomous trajectory planner are described in section “Autonomous trajectory planner.” The simulation of a space telerobot equipped with a 7-degree-of-freedom (DOF) manipulator capturing space debris is presented in section “Simulation.” Finally, the conclusions are provided in section “Conclusion.”

Motion equations

As shown in Figure 1, the space telerobot system has n + 1 rigid bodies, where the satellite base is body 0 and the manipulator contains the rest n bodies (connected by the single-DOF joint). A virtual 6-DOF joint is assumed to describe the connection relationship between the base and the earth. Before discussing the kinematics of space telerobots, some symbols and variables used in the following are defined in Table 1. Besides, the common assumptions are given as follows: (1) there are no external forces impacting on the space telerobot and (2) the initial linear and angular momentum of the space telerobot is zero.

Figure 1.

The general diagram of the space telerobot.

Table 1.

Definitions of symbols and variables used in this article.

Symbols	Definitions
${\Sigma }_I,{\Sigma }_0,{\Sigma }_E$	Inertia, base, and end effector frames
$J_i,C_i$	Joint i and the center of mass (CM) of body i
$m_i\in \mathit{R}$	Mass of body i
${\mathit{I}}_i\in {\mathit{R}}^{3\times 3}$	Inertia tensor of body i
$M\in \mathit{R}$	Mass of the whole system
${\mathit{r}}_{\mathrm{s}}\in {\mathit{R}}^3$	Position vector of the whole system CM
${\mathit{r}}_0,{\mathit{r}}_i\in {\mathit{R}}^3$	Position vector from ${\Sigma }_I$ to ${\Sigma }_0$ and ${\Sigma }_I$ to $C_i$
${\mathit{p}}_i,{\mathit{p}}_E\in {\mathit{R}}^3$	Position vector from ${\Sigma }_I$ to $J_i$ and ${\Sigma }_I$ to ${\Sigma }_E$
${\mathit{a}}_i,{\mathit{b}}_i\in {\mathit{R}}^3$	Position vector from $J_i$ to $C_i$ and $C_i$ to $J_{i+1}$
${\mathbf{\psi }}_0,{\mathbf{\psi }}_{\mathbf{E}}\in {\mathit{R}}^3$	X–Y–Z Euler angles of the base and end effector
${\mathit{\theta }}_i\in \mathit{R}$	Angle value of $J_i$
${\mathit{v}}_0,{\mathit{v}}_E\in {\mathit{R}}^3$	Linear velocity of the base and end effector
${\mathbf{\omega }}_0,{\mathbf{\omega }}_{\mathbf{E}}\in {\mathit{R}}^3$	Angular velocity of the base and end effector
${\mathit{k}}_i\in {\mathit{R}}^3$	Unit vector indicating a rotational axis of $J_i$
$\mathit{E}$	Unit matrix
$(\tilde{\mathit{x}})$	Skew-symmetric matrix of the vector $\mathit{x}$
$||\mathit{x}||$	Norm of the vector $\mathit{x}$

Coupling motion equation

Based on assumption (1), it can be deduced that the linear and angular momentum of the system is conserved. Furthermore, the dynamic coupling equation between the manipulator and the base can be derived from the momentum conservation law.

Here, the total linear and angular momentum of the whole system, $\mathit{P}$ and $\mathit{L}$ , are expressed as follows

$$\mathit{P}=\sum _{i=0}^n{m}_i{\stackrel{·}{\mathit{r}}}_i$$

(1)

$$\mathit{L}=\sum _{i=0}^n({\mathit{I}}_i{\mathbf{\omega }}_{\mathbf{i}}+{\mathit{r}}_i\times m_i{\stackrel{·}{\mathit{r}}}_i)$$

(2)

where ${\mathbf{\omega }}_{\mathbf{i}}\in {\mathit{R}}^3$ is the angular velocity of body i. It can be calculated as follows

$${\mathbf{\omega }}_{\mathbf{i}}={\mathbf{\omega }}_0+\sum _{j=1}^i{\mathit{k}}_j{\stackrel{·}{\mathit{\theta }}}_j$$

(3)

Based on assumption (2), $\mathit{P}=0$ and $\mathit{L}=0$ ; thus, equations (1) and (2) can be rearranged as follows

$$\left[\begin{array}{c}\mathit{P}\hfill \\ \multicolumn{1}{c}{\mathit{L}}\end{array}\right]=\left[\begin{array}{cc}M\mathit{E}\hfill & -M{\tilde{\mathit{r}}}_{0\mathrm{s}}\hfill \\ \multicolumn{1}{c}{M{\tilde{\mathit{r}}}_{\mathrm{s}}}& {\mathit{I}}_{\omega }\hfill \end{array}\right]\left[\begin{array}{c}{\mathit{v}}_0\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_0}\end{array}\right]+\left[\begin{array}{c}{\mathit{J}}_{T\omega }\hfill \\ \multicolumn{1}{c}{{\mathit{I}}_{\varphi }}\end{array}\right]\stackrel{·}{\mathit{\theta }}=0$$

(4)

where the intermediate parameter matrices are defined as follows

$${\mathit{J}}_{T\omega }=\sum _{i=1}^n{m}_i{\mathit{J}}_{\mathrm{Ti}}$$

(5)

$${\mathit{J}}_{\mathrm{Ti}}=[{\mathit{k}}_1\times ({\mathit{r}}_i-{\mathit{p}}_1),\dots ,{\mathit{k}}_i\times ({\mathit{r}}_i-{\mathit{p}}_i),0,\dots ,0]$$

(6)

$${\mathit{I}}_{\omega }=\sum _{i=1}^n({\mathit{I}}_i-m_i{\tilde{\mathit{r}}}_i{\tilde{\mathit{r}}}_{0i})+{\mathit{I}}_0$$

(7)

$${\mathit{I}}_{\varphi }=\sum _{i=1}^n({\mathit{I}}_i{\mathit{J}}_{\mathrm{Ri}}+m_i{\tilde{\mathit{r}}}_i{\mathit{J}}_{\mathrm{Ti}})$$

(8)

$${\mathit{J}}_{\mathrm{Ri}}=[{\mathit{k}}_1,\dots ,{\mathit{k}}_i,0,\dots ,0]$$

(9)

and ${\mathit{r}}_{0s}={\mathit{r}}_s-{\mathit{r}}_0\in {\mathit{R}}^3$ is the position vector from ${\Sigma }_I$ to the mass center of the system; $\stackrel{·}{\mathit{\theta }}=[\begin{array}{cccc}{\stackrel{·}{\mathit{\theta }}}_1\hfill & {\stackrel{·}{\mathit{\theta }}}_2\hfill & \dots \hfill & {\stackrel{·}{\mathit{\theta }}}_n\hfill \end{array}{]}^T\in {\mathit{R}}^n$ is the joint velocity vector of the manipulator; ${\mathit{r}}_{0i}={\mathit{r}}_i-{\mathit{r}}_0\in {\mathit{R}}^3$ is the position vector from ${\Sigma }_I$ to ${\mathit{J}}_i$ .

Furthermore, let ${\mathit{L}}_0\in {\mathit{R}}^3$ denote the angular momentum of the system with respect to the mass center of the base. So, $\mathit{L}$ can be written as follows

$$\mathit{L}={\mathit{L}}_0+{\mathit{r}}_0\times \mathit{P}$$

(10)

And, equation (4) can be modified as follows

$$\left[\begin{array}{c}\mathit{P}\hfill \\ \multicolumn{1}{c}{{\mathit{L}}_0}\end{array}\right]=\left[\begin{array}{cc}M\mathit{E}\hfill & -M{\tilde{\mathit{r}}}_{0s}\hfill \\ \multicolumn{1}{c}{M{\tilde{\mathit{r}}}_{0s}}& {\mathit{H}}_{\omega }\hfill \end{array}\right]\left[\begin{array}{c}{\mathit{v}}_0\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_0}\end{array}\right]+\left[\begin{array}{c}{\mathit{J}}_{T\omega }\hfill \\ \multicolumn{1}{c}{{\mathit{H}}_{\omega \varphi }}\end{array}\right]\stackrel{·}{\mathit{\theta }}=0$$

(11)

where the intermediate parameter matrices are defined as follows

$${\mathit{H}}_{\omega }=\sum _{i=1}^n({\mathit{I}}_i-m_i{\tilde{\mathit{r}}}_{0i}{\tilde{\mathit{r}}}_{0i})+{\mathit{I}}_0$$

(12)

$${\mathit{H}}_{\omega \phi }=\sum _{i=1}^n({\mathit{I}}_i{\mathit{J}}_{\mathrm{Ri}}+m_i{\tilde{\mathit{r}}}_{0i}{\mathit{J}}_{\mathrm{Ti}})$$

(13)

Then, the coupling motion equation can be obtained as follows

$$\left[\begin{array}{c}{\mathit{v}}_0\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_0}\end{array}\right]=-{\mathit{H}}_b^{-1}{\mathit{H}}_{\mathrm{bm}}\stackrel{·}{\mathit{\theta }}$$

(14)

where ${\mathit{H}}_b\in {\mathit{R}}^{6\times 6}$ is the inertia matrix of the base, ${\mathit{H}}_{\mathrm{bm}}\in {\mathit{R}}^{6\times n}$ is the inertia coupling matrix between the base and the manipulator, and ${\mathit{H}}_b$ and ${\mathit{H}}_{\mathrm{bm}}$ are, respectively, expressed as follows

$${\mathit{H}}_b=\left[\begin{array}{cc}M\mathit{E}\hfill & -M{\tilde{\mathit{r}}}_{0s}\hfill \\ \multicolumn{1}{c}{M{\tilde{\mathit{r}}}_{0s}}& {\mathit{H}}_{\omega }\hfill \end{array}\right]$$

(15)

$${\mathit{H}}_{\mathrm{bm}}=\left[\begin{array}{c}{\mathit{J}}_{T\omega }\hfill \\ \multicolumn{1}{c}{{{\mathit{H}}_{\omega }}_{\varphi }}\end{array}\right]$$

(16)

Furthermore, define ${\mathit{H}}_{\mathrm{DC}}\in {\mathit{R}}^{6\times n}$ as the dynamic coupling matrix, that is, ${\mathit{H}}_{\mathrm{DC}}=-{\mathit{H}}_b^{-1}{\mathit{H}}_{\mathrm{bm}}$ . So, equation (14) can be rewritten as follows

$$\left[\begin{array}{c}{\mathit{v}}_0\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_0}\end{array}\right]={\mathit{H}}_{\mathrm{DC}}\stackrel{·}{\mathit{\theta }}=\left[\begin{array}{c}{\mathit{H}}_{\mathrm{DCv}}\hfill \\ \multicolumn{1}{c}{{\mathit{H}}_{\mathrm{DC}\omega }}\end{array}\right]\stackrel{·}{\mathit{\theta }}$$

(17)

where ${\mathit{H}}_{\mathrm{DCv}}\in {\mathit{R}}^{3\times n}$ and ${\mathit{H}}_{\mathrm{DC}\omega }\in {\mathit{R}}^{3\times n}$ are the sub-coupling matrices corresponding to the translational and rotational terms of the base, respectively.

Here, the base attitude is represented by X–Y–Z Euler angles, that is, ${\mathbf{\psi }}_0=[\begin{array}{ccc}{\alpha }_0\hfill & {\beta }_0\hfill & {\gamma }_0\hfill \end{array}{]}^{\mathrm{T}}$ , and the relation between the angle velocity and the X–Y–Z Euler angle velocity can be calculated as follows

$${\mathbf{\omega }}_0=\left[\begin{array}{ccc}1\hfill & 0\hfill & \sin {\beta }_0\hfill \\ \multicolumn{1}{c}{0}& \cos {\alpha }_0\hfill & -\sin {\alpha }_0\cos {\beta }_0\hfill \\ \multicolumn{1}{c}{0}& \sin {\alpha }_0\hfill & \cos {\alpha }_0\cos {\beta }_0\hfill \end{array}\right]\left[\begin{array}{c}{\stackrel{·}{\alpha }}_0\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\beta }}_0}\\ \multicolumn{1}{c}{{\stackrel{·}{\gamma }}_0}\end{array}\right]={\mathit{N}}_0{\stackrel{·}{\psi }}_0$$

(18)

Combining the lower part of equation (17) with equation (18), the base attitude kinematics in X–Y–Z Euler angles can be derived as follows

$${\stackrel{·}{\psi }}_0={\mathit{N}}_0^{-1}{\mathit{H}}_{\mathrm{DC}\omega }\stackrel{·}{\mathit{\theta }}$$

(19)

Furthermore, the attitude of the base at the final time $t_f$ can be calculated as follows

$${\mathbf{\psi }}_0(t_f)=\underset{t_0}{\overset{t_f}{\int }}{\mathit{N}}_0^{-1}{\mathit{H}}_{\mathrm{DC}\omega }\stackrel{·}{\mathit{\theta }}\mathrm{dt}+{\mathbf{\psi }}_0(t_0)$$

(20)

where ${\mathbf{\psi }}_0(t_0)$ is the attitude of the base at the initial time $t_0$ .

Kinematic equation of the space telerobot

As shown in Figure 1, the position vector of the end effector in the inertia frame is defined as follows

$${\mathit{p}}_E={\mathit{r}}_0+{\mathit{b}}_0+\sum _{i=1}^n({\mathit{a}}_i+{\mathit{b}}_i)$$

(21)

Differentiating equation (21) with respect to time, the linear velocity of the end effector can be derived as follows

$${\mathit{v}}_E={\mathit{v}}_0+{\mathbf{\omega }}_0\times ({\mathit{p}}_E-{\mathit{r}}_0)+\sum _{i=1}^n[{\mathit{k}}_i\times ({\mathit{p}}_E-{\mathit{p}}_i)]{\stackrel{·}{\mathit{\theta }}}_i$$

(22)

Besides, the angular velocity of the end effector is expressed as follows

$${\mathbf{\omega }}_{\mathbf{E}}={\mathbf{\omega }}_0+\sum _{i=1}^n{\mathit{k}}_i{\stackrel{·}{\mathit{\theta }}}_i$$

(23)

Combining equations (22) and (23), the general differential kinematics of the free-flying space telerobot can be expressed as follows

$$\left[\begin{array}{c}{\mathit{v}}_E\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_{\mathbf{E}}}\end{array}\right]={\mathit{J}}_B\left[\begin{array}{c}{\mathit{v}}_0\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_0}\end{array}\right]+{\mathit{J}}_M\stackrel{·}{\mathit{\theta }}$$

(24)

where ${\mathit{J}}_B$ and ${\mathit{J}}_M$ are the Jacobian matrices for the base motion–dependent term and the manipulator-dependent term, respectively, defined as follows

$${\mathit{J}}_B=\left[\begin{array}{cc}\mathit{E}\hfill & -{\mathit{p}}_{0E}\hfill \\ \multicolumn{1}{c}{0}& \mathit{E}\hfill \end{array}\right]$$

(25)

$${\mathit{J}}_M=\left[\begin{array}{ccc}{\mathit{k}}_1\times ({\mathit{p}}_E-{\mathit{p}}_1)\hfill & \cdots \hfill & {\mathit{k}}_n\times ({\mathit{p}}_E-{\mathit{p}}_n)\hfill \\ \multicolumn{1}{c}{{\mathit{k}}_1}& \cdots \hfill & {\mathit{k}}_n\hfill \end{array}\right]$$

(26)

where ${\mathit{p}}_{0E}={\mathit{p}}_E-{\mathit{r}}_0\in {\mathit{R}}^3$ is the position vector from ${\Sigma }_E$ to ${\Sigma }_0$ . Furthermore, substituting equation (17) into equation (24), the differential kinematics of the free-floating space telerobot can be derived as follows

$$\left[\begin{array}{c}{\mathit{v}}_E\hfill \\ \multicolumn{1}{c}{{\mathbf{\omega }}_{\mathbf{E}}}\end{array}\right]=({\mathit{J}}_B{\mathit{H}}_{\mathrm{DC}}+{\mathit{J}}_M)\stackrel{·}{\mathit{\theta }}={\mathit{J}}_G\stackrel{·}{\mathit{\theta }}=\left[\begin{array}{c}{\mathit{J}}_{\mathrm{Gv}}\hfill \\ \multicolumn{1}{c}{{\mathit{J}}_{G\mathbf{\omega }}}\end{array}\right]\stackrel{·}{\mathit{\theta }}$$

(27)

where the coefficient matrix $J_G\in R^{6\times n}$ is the GJM of the free-floating space telerobot, and ${\mathit{J}}_{\mathrm{Gv}}$ and ${\mathit{J}}_{G\omega }\in {\mathit{R}}^{3\times n}$ denote the sub-GJMs corresponding to the linear and angular velocity terms of the end effector, respectively. Note that ${\mathit{J}}_G$ contains the dynamic coupling matrix which is related to the dynamic parameters of the whole system, so the singularities of the GJM are called the dynamic singularities. ²⁰

Furthermore, the position of the end effector at the final time $t_f$ can be calculated as follows

$${\mathit{p}}_E(t_f)=\underset{t_0}{\overset{t_f}{\int }}{\mathit{J}}_{\mathrm{Gv}}\stackrel{·}{\mathit{\theta }}\mathrm{dt}+{\mathit{p}}_E(t_0)$$

(28)

where ${\mathit{p}}_E(t_0)$ is the position of the end effector at the initial time $t_0$ . Similarly, the attitude kinematics of the end effector can be derived as follows

$${\stackrel{·}{\psi }}_E={\mathit{N}}_E^{-1}{\mathit{J}}_{G\omega }\stackrel{·}{\mathit{\theta }}$$

(29)

And, the attitude of the end effector at the final time $t_f$ can be calculated as follows

$${\mathbf{\psi }}_{\mathbf{E}}(t_f)=\underset{t_0}{\overset{t_f}{\int }}{\mathit{N}}_E^{-1}{\mathit{J}}_{G\omega }\stackrel{·}{\mathit{\theta }}\mathrm{dt}+{\mathbf{\psi }}_{\mathbf{E}}(t_0)$$

(30)

where ${\mathbf{\psi }}_{\mathbf{E}}(t_0)$ is the initial attitude of the end effector.

Analysis of the capture task

Brief of the space teleprogramming framework

The capture task is a necessary step in the ADR mission, where the task profile of the space telerobot capturing a floating debris in LEO is illustrated in Figure 2(a). The communication data should be transmitted via one or more tracking and data relay satellites (TDRS) in geosynchronous orbit (GEO), which is the main reason why the cumulative time delay is several seconds. ⁴ As is mentioned above, the teleprogramming control provides an achievable control framework to get rid of the communication delay. In fact, it had been successfully performed in the ROTEX project. ⁵ In Figure 2(b), a general control framework is concluded to describe the workflow of teleprogramming robots, ¹⁰ which consists of three levels: the decision-making level, the analysis and planning level, and the execution level. Meanwhile, it contains three control loops: the main control loop (the supervisor loop), the simulation loop, and the remote autonomy loop. The ground operator in decision-making level defines the tasks, interacts with the simulation system, monitors the feedback from the simulation loop and main loop, and updates the commands. The ground simulator provides a virtualized operating environment where the task can be performed by operators or the ground controller. Moreover, the ground controller continuously records the operating information and translates them into symbolic or normalized-value motion instructions. As for the semi-autonomous platform, that is, the space robot, it takes charge of the real-time hardware servo control to execute the motion instructions.

Figure 2.

The space teleprogramming robot: (a) the schematic diagram and (b) the control framework.

Task analysis and problem statement

In this subsection, the capture task is analyzed in detail to make the problem clear. As shown in Figure 2(a), before the manipulator executes the operation, the space telerobot has been navigated to an appropriate rendezvous zone where the manipulator can reach the debris, and the on-board vision devices can accurately measure the position and attitude information of the debris. Then, the space telerobot sends the measured information to the earth working station which takes charge of planning feasible trajectories in joint space to meet the grasping constraints. Here, the grasping constraints can be described by the following: (1) grasping position: ${\mathit{p}}_E(t_f)={\mathit{p}}_{\mathrm{TH}}$ and (2) grasping attitude: ${\mathbf{\psi }}_{\mathbf{E}}(t_f)={\mathbf{\psi }}_{\mathbf{TH}}$ , where ${\mathit{p}}_{\mathrm{TH}}\in {\mathit{R}}^3$ and ${\mathbf{\psi }}_{\mathbf{TH}}\in {\mathit{R}}^3$ represent the position and attitude constraints of the grasping handle on the debris, respectively.

Meanwhile, in order to guarantee the telemetry data transmission, the attitude disturbance of the base antenna should be controlled within a certain range. Here, the change of the base attitude, $\mathrm{\Delta }{\mathbf{\psi }}_0=[\begin{array}{ccc}\mathrm{\Delta }{\alpha }_0\hfill & \mathrm{\Delta }{\beta }_0\hfill & \mathrm{\Delta }{\gamma }_0\hfill \end{array}{]}^T$ , can be calculated as follows: $\mathrm{\Delta }{\mathbf{\psi }}_0={\mathbf{\psi }}_0(t_f)-{\mathbf{\psi }}_0(t_0)$ . The corresponding restricted conditions can be expressed as follows: $\left|\mathrm{\Delta }{\alpha }_0\right|\le \mathrm{\Delta }{\alpha }_{0\mathrm{ab}}$ , $\left|\mathrm{\Delta }{\beta }_0\right|\le \mathrm{\Delta }{\beta }_{0\mathrm{ab}}$ , and $\left|\mathrm{\Delta }{\gamma }_0\right|\le \mathrm{\Delta }{\gamma }_{0\mathrm{ab}}$ , where $\mathrm{\Delta }{\alpha }_{0\mathrm{ab}}$ , $\mathrm{\Delta }{\beta }_{0\mathrm{ab}}$ , and $\mathrm{\Delta }{\gamma }_{0\mathrm{ab}}$ are the absolute boundaries of the attitude disturbance.

In addition, the motion performance of the manipulator system largely depends on the physical constraints in joint space. So, the physical constraints should be considered in the trajectory planning, which are described as follows: (1) joint position: ${\mathbf{\theta }}_{min}\le \mathbf{\theta }\mathbf{\le }{\mathbf{\theta }}_{max}$ , (2) joint velocity: ${\stackrel{·}{\mathit{\theta }}}_{min}\le \stackrel{·}{\mathit{\theta }}\le {\stackrel{·}{\mathit{\theta }}}_{max}$ , and (3) joint acceleration: ${\stackrel{··}{\mathit{\theta }}}_{min}\le \stackrel{··}{\mathit{\theta }}\le {\stackrel{··}{\mathit{\theta }}}_{max}$ . Here, ${\mathbf{\theta }}_{min}$ , ${\stackrel{·}{\mathit{\theta }}}_{min}$ , and ${\stackrel{··}{\mathit{\theta }}}_{min}\in {\mathit{R}}^n$ are the lower limits of the position, velocity, and acceleration in joint space, respectively; ${\mathbf{\theta }}_{max}$ , ${\stackrel{·}{\mathit{\theta }}}_{max}$ , and ${\stackrel{··}{\mathit{\theta }}}_{max}\in {\mathit{R}}^n$ are the upper limits of the position, velocity, and acceleration in joint space, respectively.

Principally, the planning problem for which joint trajectories need to be determined to satisfy constraints can be formulated as a constrained optimization problem as follows:

To determine joint functions: $\mathbf{\theta }\mathbf{(}\mathbf{t}\mathbf{)}$ ;

To minimize the objective function

$$\mathrm{\Gamma }(\mathbf{\theta }\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{)}\mathbf{=}\mathbf{\Vert }{\mathbf{\lambda }}_{\mathbf{p}}\mathbf{\delta }{\mathit{p}}_{\mathbf{E}}({\mathbf{t}}_{\mathbf{f}})\mathbf{\Vert }\mathbf{+}\mathbf{\Vert }{\mathbf{\lambda }}_{\mathbf{\psi }}\mathbf{\delta }{\mathbf{\psi }}_{\mathbf{E}}({\mathbf{t}}_{\mathbf{f}})\mathbf{\Vert }$$

(31)

where $\delta {\mathit{p}}_E(t_f)$ and $\delta {\mathbf{\psi }}_{\mathbf{E}}(t_f)$ are the grasping position error and the grasping attitude error, respectively; ${\mathbf{\lambda }}_{\mathbf{p}}$ and ${\mathbf{\lambda }}_{\mathbf{\psi }}$ are weight coefficient matrices defined as follows

$${\mathbf{\lambda }}_{\mathbf{p}}=\mathrm{diag}\left(\begin{array}{ccc}\frac{1}{{\xi }_{\mathrm{Px}}}\hfill & \frac{1}{{\xi }_{\mathrm{Py}}}\hfill & \frac{1}{{\xi }_{\mathrm{Pz}}}\hfill \end{array}\right)$$

(32a)

$${\mathbf{\lambda }}_{\mathbf{\psi }}=\mathrm{diag}\left(\begin{array}{ccc}\frac{1}{{\xi }_{\psi \alpha }}\hfill & \frac{1}{{\xi }_{\psi \beta }}\hfill & \frac{1}{{\xi }_{\psi \gamma }}\hfill \end{array}\right)$$

(32b)

where ${\xi }_{\mathrm{Px}}$ , ${\xi }_{\mathrm{Py}}$ , and ${\xi }_{\mathrm{Pz}}$ represent the permissible position accuracy errors in three axis directions; ${\xi }_{\psi \alpha }$ , ${\xi }_{\psi \beta }$ , and ${\xi }_{\psi \gamma }$ represent the permissible attitude accuracy errors.

Subject to the inequality constraints

$$g_i(\mathbf{\theta }\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{)}\mathbf{:}\mathbf{\{}\begin{array}{c}{\mathbf{\theta }}_{min}\mathbf{\le }\mathbf{\theta }\mathbf{\le }{\mathbf{\theta }}_{max}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathbf{\theta }}}_{min}\mathbf{\le }\stackrel{·}{\mathbf{\theta }}\mathbf{\le }{\stackrel{·}{\mathbf{\theta }}}_{max}}\\ \multicolumn{1}{c}{{\stackrel{··}{\mathbf{\theta }}}_{min}\mathbf{\le }\stackrel{··}{\mathbf{\theta }}\mathbf{\le }{\stackrel{··}{\mathbf{\theta }}}_{max}}\end{array}$$

(33a)

$$g_i(\mathrm{\Delta }{\mathbf{\psi }}_0):\{\begin{array}{c}\left|\mathrm{\Delta }{\alpha }_0\right|\le \mathrm{\Delta }{\alpha }_{0\mathrm{ab}}\hfill \\ \multicolumn{1}{c}{\left|\mathrm{\Delta }{\beta }_0\right|\le \mathrm{\Delta }{\beta }_{0\mathrm{ab}}}\\ \multicolumn{1}{c}{\left|\mathrm{\Delta }{\gamma }_0\right|\le \mathrm{\Delta }{\gamma }_{0\mathrm{ab}}}\end{array}$$

(33b)

• Subject to the equality constraints

$$h_i(\mathbf{\theta }\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{)}\mathbf{:}\mathbf{\{}\begin{array}{c}\mathbf{\theta }\mathbf{(}{\mathbf{t}}_0\mathbf{)}\mathbf{=}{\mathbf{\theta }}_{\mathbf{i}0}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\mathbf{\theta }}\mathbf{(}{\mathbf{t}}_0\mathbf{)}\mathbf{=}\stackrel{··}{\mathbf{\theta }}\mathbf{(}{\mathbf{t}}_0\mathbf{)}\mathbf{=}0}\\ \multicolumn{1}{c}{\stackrel{·}{\mathbf{\theta }}\mathbf{(}{\mathbf{t}}_{\mathbf{f}}\mathbf{)}\mathbf{=}\stackrel{··}{\mathbf{\theta }}\mathbf{(}{\mathbf{t}}_{\mathbf{f}}\mathbf{)}\mathbf{=}0}\end{array}$$

(34)

where ${\mathbf{\theta }}_{\mathbf{i}0}$ is the initial value of $\mathbf{\theta }\mathbf{(}\mathbf{t}\mathbf{)}$ . If the value of $\mathrm{\Gamma }(\mathbf{\theta }\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{)}$ is less than 1, one can infer that the current joint functions satisfy the grasping constraints. Note that the planned trajectories should be continuous, smooth, and free of singularities. However, the dynamic coupling and the dynamic singularities are critical problems to affect the operational performance of space telerobots, leading to the fact that they cannot be easily manipulated like base-fixed robots. With the development of artificial intelligence (AI) technologies, some AI methods (e.g. fuzzy logic, ²¹ neural network, ²² and evolutionary algorithms) have been widely employed in the telerobot field with providing innovative approaches for controlling or planning telerobots. ²³

Autonomous trajectory planner

Through the above analysis, the motion planning problem can be converted into a constrained optimization problem. Some methods^24,25 based on the optimal control theory are proposed to solve this problem, where the planned trajectories can satisfy the equality and inequality constraints, but sometimes, the trajectories are not smooth enough. Another solution is the direct optimization method,^4,26 which adopts the parameter optimization technique and search algorithms (e.g. sequential quadratic programming (SQP), ²⁷ genetic algorithm (GA), ²⁸ and particle swarm optimization (PSO) ²⁹ ) to derive feasible trajectories. Compared with the first method, the second method is more flexible and targeted. Therefore, an autonomous planner based on the direct optimization method is developed to generate feasible trajectories for space telerobots to execute the point-to-point grasping task.

Parameterization of the joint trajectories

Many functions can be used to parameterize the joint trajectories, such as the polynomial function, Bézier curve, and other functions. Here, the polynomial functions are employed to parameterize the joint position, velocity, and acceleration trajectories as follows

$$\{\begin{array}{c}{\mathit{\theta }}_i(t)=\sum _{k=0}^j{f}_{\mathrm{ik}}{t}^k\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{\theta }}}_i(t)=\sum _{k=1}^{j}k·f_{\mathrm{ik}}{t}^{k-1}}\\ \multicolumn{1}{c}{{\stackrel{··}{\mathit{\theta }}}_i(t)=\sum _{k=2}^{j}k·(k-1)·f_{\mathrm{ik}}{t}^{k-2}}\end{array},(i=1,2,\dots ,n)$$

(35)

where $f_{\mathrm{ik}}$ is the coefficient of the ith joint function, and j is the order of the polynomial function which depends on the number of equality constraints in each joint. For capturing a floating debris, the fifth-order polynomial function can be used to parameterize each joint trajectory. Substituting equation (34) into equation (35), the following expression can be obtained

$$\{\begin{array}{c}{f}_{i0}={\mathit{\theta }}_{i0}\\ f_{i1}=f_{i2}=0\\ f_{i3}=\frac{5}{3}{f}_{i5}{t}_f^2\\ f_{i4}=-\frac{5}{2}{f}_{i5}{t}_f\end{array}\Rightarrow \{\begin{array}{c}{\mathit{\theta }}_i(t)=f_{i5}\left(t^5-\frac{5}{2}{t}_f{t}^4+\frac{5}{3}{t}_f^2{t}^3\right)+{\mathit{\theta }}_{i0}\\ {\stackrel{.}{\mathit{\theta }}}_i(t)=f_{i5}(5t^4-10t_f{t}^3+5t_f^2{t}^2)\\ {\stackrel{··}{\mathit{\theta }}}_i(t)=f_{i5}(20t^3-30t_f{t}^2+10t_f^{2}t)\end{array}$$

(36)

where $f_{i5}$ is the undetermined parameter in each joint function.

Normalization treatment

After the parameterization treatment, a feasible searching region can be defined through the normalization treatment. ³⁰ Here, the physical constraints are regarded as the criteria to determine the searching boundary. First, the normalized time is defined as follows

$$\tau =\frac{t}{T}$$

(37)

where $\tau$ is the normalized time $(\tau \in [\begin{array}{cc}0\hfill & 1\hfill \end{array}])$ ; $T$ is the total execution time, $T=t_f-t_0$ .

Substituting equation (37) into equation (36), the corresponding joint functions can be written as follows

$$\{\begin{array}{c}{\mathit{\theta }}_i(t)={\hat{f}}_{i5}·F_5(\tau )+{\mathit{\theta }}_{i0}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathit{\theta }}}_i(t)={\hat{f}}_{i5}·{\stackrel{·}{F}}_5(\tau )}\\ \multicolumn{1}{c}{{\stackrel{··}{\mathit{\theta }}}_i(t)={\hat{f}}_{i5}·{\stackrel{··}{F}}_5(\tau )}\end{array}$$

(38)

where ${\hat{f}}_{i5}=f_{i5}·T^5$ ; $F_5(\tau )$ , ${\stackrel{·}{F}}_5(\tau )$ , and ${\stackrel{··}{F}}_5(\tau )$ are the functions of the normalized time, expressed as follows

$$\{\begin{array}{c}{F}_5(\tau )={\tau }^5-\frac{5}{2}{\tau }^4+\frac{5}{3}{\tau }^3\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{F}}_5(\tau )=5{\tau }^4-10{\tau }^3+5{\tau }^2}\\ \multicolumn{1}{c}{{\stackrel{··}{F}}_5(\tau )=20{\tau }^3-30{\tau }^2+10\tau }\end{array}$$

(39)

Furthermore, since the normalized time $\tau$ is bounded, one can deduce that $F_5(\tau )$ , ${\stackrel{·}{F}}_5(\tau )$ , and ${\stackrel{··}{F}}_5(\tau )$ are bounded functions. They can be expressed as follows

$$\{\begin{array}{c}|F_5(\tau )|\le F_{5max}\hfill \\ \multicolumn{1}{c}{|{\stackrel{·}{F}}_5(\tau )|\le {\stackrel{·}{F}}_{5max}}\\ \multicolumn{1}{c}{|{\stackrel{··}{F}}_5(\tau )|\le {\stackrel{··}{F}}_{5max}}\end{array}$$

(40)

where $F_{5max}$ , ${\stackrel{·}{F}}_{5max}$ , and ${\stackrel{··}{F}}_{5max}$ are their absolute boundaries. According to equations (33a), (38), and (40), the boundaries of $f_{i5}$ can be derived as follows

$$f_{i5min}\le f_{i5}\le f_{i5max}$$

(41)

$$\{\begin{array}{c}{f}_{i5min}=max\left\{\begin{array}{ccc}\frac{{{\mathit{\theta }}_i}_{min}-{\mathit{\theta }}_{i0}}{F_{5max}·T^5}\hfill & \frac{{\stackrel{·}{\mathit{\theta }}}_{imin}}{{\stackrel{·}{F}}_{5max}·T^4}\hfill & \frac{{\stackrel{··}{\mathit{\theta }}}_{imin}}{{\stackrel{··}{F}}_{5max}·T^3}\hfill \end{array}\right\}\hfill \\ \multicolumn{1}{c}{f_{i5max}=min\left\{\begin{array}{ccc}\frac{{{\mathit{\theta }}_i}_{max}-{\mathit{\theta }}_{i0}}{F_{5max}·T^5}\hfill & \frac{{\stackrel{·}{\mathit{\theta }}}_{imax}}{{\stackrel{·}{F}}_{5max}·T^4}\hfill & \frac{{\stackrel{··}{\mathit{\theta }}}_{imax}}{{\stackrel{··}{F}}_{5max}·T^3}\hfill \end{array}\right\}}\end{array}$$

(42)

where $f_{i5min}$ and $f_{i5max}$ are the upper and lower boundaries of $f_{i5}$ , respectively.

Penalty function method

After the normalization treatment, the inequality constraints in joint space are released. The rest inequality constraints are the attitude disturbance constraints of the base, where the penalty function method ³¹ is adopted to handle the attitude disturbance constraints as follows

$$G_1(\mathrm{\Delta }{\alpha }_0)=\{\begin{array}{cc}{\mu }_1·\sqrt{{(\left|\mathrm{\Delta }{\alpha }_0\right|-\mathrm{\Delta }{\alpha }_{0\mathrm{ab}})}^2}\hfill & \left|\mathrm{\Delta }{\alpha }_0\right|>\mathrm{\Delta }{\alpha }_{0\mathrm{ab}}\hfill \\ \multicolumn{1}{c}{0}& \left|\mathrm{\Delta }{\alpha }_0\right|\le \mathrm{\Delta }{\alpha }_{0\mathrm{ab}}\hfill \end{array}$$

(43a)

$$G_2(\mathrm{\Delta }{\beta }_0)=\{\begin{array}{cc}{\mu }_2·\sqrt{{(\left|\mathrm{\Delta }{\beta }_0\right|-\mathrm{\Delta }{\beta }_{0\mathrm{ab}})}^2}\hfill & \left|\mathrm{\Delta }{\beta }_0\right|>\mathrm{\Delta }{\beta }_{0\mathrm{ab}}\hfill \\ \multicolumn{1}{c}{0}& \left|\mathrm{\Delta }{\beta }_0\right|\le \mathrm{\Delta }{\beta }_{0\mathrm{ab}}\hfill \end{array}$$

(43b)

$$G_3(\mathrm{\Delta }{\gamma }_0)=\{\begin{array}{cc}{\mu }_3·\sqrt{{(\left|\mathrm{\Delta }{\gamma }_0\right|-\mathrm{\Delta }{\gamma }_{0\mathrm{ab}})}^2}\hfill & \left|\mathrm{\Delta }{\gamma }_0\right|>\mathrm{\Delta }{\gamma }_{0\mathrm{ab}}\hfill \\ \multicolumn{1}{c}{0}& \left|\mathrm{\Delta }{\gamma }_0\right|\le \mathrm{\Delta }{\gamma }_{0\mathrm{ab}}\hfill \end{array}$$

(43c)

where $G_1(\mathrm{\Delta }{\alpha }_0)$ , $G_2(\mathrm{\Delta }{\beta }_0)$ , and $G_3(\mathrm{\Delta }{\gamma }_0)$ are penalty functions corresponding to the Euler angle changes of the base; ${\mu }_1$ , ${\mu }_2$ , and ${\mu }_3$ are positive constants, also called penalty factors.

Furthermore, by merging equation (31) and equation (43), a new objective function can be formulated as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Gamma }}_{\mathrm{new}}=\Vert {\mathbf{\lambda }}_{\mathbf{p}}\delta {\mathit{p}}_E(t_f)\Vert +\Vert {\mathbf{\lambda }}_{\mathbf{\psi }}\delta {\mathbf{\psi }}_{\mathbf{E}}(t_f)\Vert +G_1(\mathrm{\Delta }{\alpha }_0)}\\ \multicolumn{1}{c}{+G_2(\mathrm{\Delta }{\beta }_0)+G_3(\mathrm{\Delta }{\gamma }_0)}\end{array}$$

(44)

As ${\mathrm{\Gamma }}_{\mathrm{new}}$ already contains the attitude disturbance penalty functions, it can drive solutions from restricted areas through the action of the objective function-based criterion.

Synthesizing the parameterization, normalization, and penalty function treatments above, the constrained optimization problem can be further converted into an unconstrained optimization problem as follows:

• To determine: $f_{i5},(i=1,2,\dots ,n)$ ;

• To minimize: ${\mathrm{\Gamma }}_{\mathrm{new}}(f_{i5})$ .

The feasible region $[\begin{array}{cc}{f}_{i5min}\hfill & f_{i5max}\hfill \end{array}]$ has been determined according to equation (42), and the optimal solution can be obtained by the search algorithm.

DE

As one of the most efficient population-based evolutionary algorithms, DE is an improved version of GA and has been applied in many practical problems. ³² For its simple searching mechanism, excellent convergence property, and suitability for parallelization, the DE algorithm is employed to search the optimal joint parameters in the capture task.

In the DE algorithm, a parent individual from the current generation is set as target vector; a mutant vector obtained by the differential mutation operation is named donor vector; an offspring vector generated by mingling the target vector and the donor vector is called trial vector; and through the competition between the target vector and the trial vector, DE realizes the overall evolution of the population.

Specifically, the DE/rand/1/bin algorithm ³³ is adopted here, and its detailed steps are introduced as follows. For a n-dimensional optimization, first initialize a population of NP individuals ${\mathit{p}}_G=\{{\overrightarrow{X}}_{1,G},{\overrightarrow{X}}_{2,G},\dots ,{\overrightarrow{X}}_{\mathrm{NP},G}\}$ , where G denotes the generation or iteration number (for the first generation, $G=0$ ); ${\overrightarrow{X}}_{j,G}$ denotes the jth individual, ${\overrightarrow{X}}_{j,G}=[x_{1,j,G},x_{2,j,G},\dots ,x_{n,j,G}](j=[1,2,\dots ,\mathrm{NP}])$ , and each individual uniformly distributed in the range $[{\overrightarrow{X}}_{min},{\overrightarrow{X}}_{max}]$ , where ${\overrightarrow{X}}_{min}=\{x_{1,min},x_{2,min},\dots ,x_{n,min}\}$ and ${\overrightarrow{X}}_{max}=\{x_{1,max},x_{2,max},\dots ,x_{n,max}\}$ . Then, execute the mutation operation, and the donor vector is created as follows

$${\overrightarrow{V}}_{j,G}={\overrightarrow{X}}_{j,G}^{r1}+F·({\overrightarrow{X}}_{j,G}^{r2}-{\overrightarrow{X}}_{j,G}^{r3})$$

(45)

where ${\overrightarrow{V}}_{j,G}=\{v_{1,j,G},v_{2,j,G},\dots ,v_{n,j,G}\}$ is the donor vector; ${\overrightarrow{X}}_{j,G}^{r1}$ , ${\overrightarrow{X}}_{j,G}^{r2}$ , and ${\overrightarrow{X}}_{j,G}^{r3}$ are three distinct parameter vectors sampled randomly from the current population; and $F$ is the scale factor.

After obtaining the donor vector, the crossover operation will be carried out to maintain the diversity of the population. Here, ${\overrightarrow{X}}_{j,G}=[x_{1,j,G},x_{2,j,G},\dots ,x_{n,j,G}]$ is set as the target vector, and the trial vector ${\overrightarrow{U}}_{j,G}=\{u_{1,j,G},u_{2,j,G},\dots ,u_{n,j,G}\}$ can be generated in the following way

$$u_{i,j,G}=\{\begin{array}{cc}{v}_{i,j,G}\hfill & {\mathrm{randl}}_{i,j}\le \mathrm{CR}\hfill \\ \multicolumn{1}{c}{x_{i,j,G}}& {\mathrm{randl}}_{i,j}>\mathrm{CR}\hfill \end{array},(i=1,2,\dots ,n)$$

(46)

where ${\mathrm{randl}}_{i,j}$ is a randomly generated number between 0 and 1 and $\mathrm{CR}$ is the crossover rate, $\mathrm{CR}\in [0,1]$ .

Next, the selection operation will determine whether the target vector or the trial vector survives to the next generation, which can be expressed as follows

$${\overrightarrow{X}}_{j,G+1}=\{\begin{array}{cc}{\overrightarrow{U}}_{j,G}\hfill & \mathrm{if}\mathrm{\Gamma }({\overrightarrow{U}}_{j,G})\le \mathrm{\Gamma }({\overrightarrow{X}}_{j,G})\hfill \\ \multicolumn{1}{c}{{\overrightarrow{X}}_{j,G}}& \mathrm{if}\mathrm{\Gamma }({\overrightarrow{U}}_{j,G})>\mathrm{\Gamma }({\overrightarrow{X}}_{j,G})\hfill \end{array}$$

(47)

where ${\overrightarrow{X}}_{j,G+1}$ is the selected individual in the next generation and $\mathrm{\Gamma }(\overrightarrow{X})$ is an objective function to be minimized.

Planning procedures

The planning procedures are designed based on the forward differential kinematics and the DE algorithm, as outlined in Figure 3. The main merit of the proposed planner is that the calculation of the GJM inversion is not involved in the planning process, so the proposed planner can avoid the dynamic singularity, and the planned trajectories are free of singularities. Furthermore, the detailed explanations of each step are provided as follows:

• Step 1. Parameterization of the joint trajectories

According to equation (36), the fifth-order polynomial function is adopted to parameterize the joint trajectories.

• Step 2. Parameter settings and the initialization

Here, the parameters can be divided two categories: the space telerobot parameters and the control parameters of DE. The first category contains the following: the initial position of the base ${\mathit{r}}_0(t_0)$ , the initial attitude of the base ${\mathbf{\psi }}_0(t_0)$ , the initial configuration of the manipulator $\mathbf{\theta }\mathbf{(}{\mathbf{t}}_0\mathbf{)}$ , the initial position of the end effector ${\mathit{p}}_E(t_0)$ , the initial attitude of the end effector ${\mathbf{\psi }}_{\mathbf{E}}(t_0)$ , and the total execution time $T$ . The second category contains the following: the scale factor $F$ , the crossover rate $\mathrm{CR}$ , the population size NP, and the maximum iterations $G_{max}$ . And each individual is set as ${\overrightarrow{X}}_{j,G}=[f_{15,j,G},f_{25,j,G},\dots ,f_{n5,j,G}]$ , where $f_{i5}$ is the undetermined parameter in each joint trajectory.

• Step 3. Judgment of the maximum iterations

If the current number of iterations is less than the maximum iterations, the planning procedures will continue optimizing, or it will output the optimal joint trajectories.

• Step 4. Calculation of the velocity terms of the base and the end effector

The velocity term of the base can be calculated using equation (17), and it will be used to construct the GJM for calculating the velocity term of the end effector, that is, equation (27).

• Step 5. Calculation of the base attitude and the end effector pose

The base attitude can be calculated using equation (20), and the end effector pose can be calculated using equations (28) and (30).

• Step 6. Calculation of the objective function

The objective function can be calculated using equation (44).

• Step 7. Judgment of the grasping constraints

If the value of the objective function is less than 1, it could be deduced that the current joint functions can satisfy the grasping constraints.

• Step 8. Record the satisfactory joint trajectory functions and output the optimal candidate

The satisfactory results will be recorded in the database, and the optimal candidate will be output when the number of iterations reaches the maximum.

Figure 3.

Planning procedures.

Simulation

In this section, a space telerobot with a 7-DOF manipulator is studied to verify the trajectory planner. Each body coordinate of the telerobot system is defined in Figure 4. The Denavit–Hartenberg (D-H) parameters of the 7-DOF manipulator and the dynamic parameters of the system are provided in Tables 2 and 3, respectively.

Figure 4.

The space telerobot system.

Table 2.

The D-H parameters of the 7-DOF space manipulator.

Link no.	${\alpha }_{i-1}(\circ )$	$a_{i-1}(\mathrm{m})$	$d_i(\mathrm{m})$	${\mathit{\theta }}_i(\circ )$
1	0.0	0.5	1.0	0.0
2	−90	0.0	0.25	45
3	90	0.0	1.8	0.0
4	90	0.0	0.25	135
5	90	0.0	1.8	180
6	90	0.0	0.25	45
7	90	0.0	0.2	0.0
End effector	0.0	0.0	0.25	0.0

D-H: Denavit–Hartenberg; DOF: degree of freedom.

Table 3.

Dynamic parameters of 7-DOF space manipulator.

Parameter	Base	Link 1	Link 2	Link 3	Link 4	Link 5	Link 6	Link 7 with EE
m(kg)	400	3.0	10	3.0	10	3.0	3.0	4.0
$I_{\mathrm{xx}}(\mathrm{kg}·m^2)$	80	0.32	1.0	0.32	1.0	0.32	0.32	0.2
$I_{\mathrm{xy}}(\mathrm{kg}·m^2)$	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
$I_{\mathrm{xz}}(\mathrm{kg}·m^2)$	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
$I_{\mathrm{yy}}(\mathrm{kg}·m^2)$	100	0.21	0.45	0.32	0.45	0.32	0.32	0.2
$I_{\mathrm{yz}}(\mathrm{kg}·m^2)$	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
$I_{\mathrm{zz}}(\mathrm{kg}·m^2)$	100	0.32	1.0	0.21	1.0	0.2	0.21	0.1
a x (m)	−	0.0	0.0	0.0	0.0	0.0	0.0	0.0
a y (m)	−	0.125	−0.9	−0.125	−0.9	0.125	0.0	0.125
a z (m)	−	0.0	0.0	0.0	0.0	0.0	0.1	0.0
b x (m)	0.5	0.0	0.0	0.0	0.0	0.0	0.0	0.0
b y (m)	0.0	0.125	−0.9	−0.125	−0.9	0.125	0.0	0.125
b z (m)	1.0	0.0	0.0	0.0	0.0	0.0	0.1	0.0

DOF: degree of freedom; EE: end effector.

The pose information of the debris is assumed as follows: ${\mathit{p}}_{\mathrm{TH}}=[0.2,-0.25,1.0](\mathrm{m})$ and ${\mathbf{\psi }}_{\mathbf{TH}}=[180,45,180](\circ )$ . Accordingly, the initial parameters of the space telerobot are set as follows: the initial position of the base ${\mathit{r}}_0(t_0)=[-3,0,0](\mathrm{m})$ , the initial attitude of the base ${\mathbf{\psi }}_0(t_0)=[0,0,0](\circ )$ , the initial configuration of the manipulator $\mathbf{\theta }\mathbf{(}{\mathbf{t}}_0\mathbf{)}\mathbf{=}\mathbf{[}0\mathbf{,}45\mathbf{,}0\mathbf{,}135\mathbf{,}180\mathbf{,}45\mathbf{,}0\mathbf{]}\mathbf{(}\mathbf{\circ }\mathbf{)}$ , the initial position of the end effector ${\mathit{p}}_E(t_0)=[0.89,-0.25,1.95](\mathrm{m})$ , and the initial attitude of the end effector ${\mathbf{\psi }}_{\mathbf{E}}(t_0)=[180,45,180](\circ )$ , which is identical to ${\mathbf{\psi }}_{\mathbf{TH}}$ for intuitively grasping. Meanwhile, the physical limitations of joints are given as follows:

1. Joint position

$$\begin{array}{c}{\mathbf{\theta }}_{min}=\hfill \\ \multicolumn{1}{c}{[-180,-40,-220,-275,-275,-275,-360](\circ )}\\ {\mathbf{\theta }}_{max}=[180,140,220,275,275,275,360](\circ )\hfill \end{array}$$

• 2. Joint velocity

$$\begin{array}{c}{\stackrel{·}{\mathit{\theta }}}_{min}=[-30,-30,-45,-45,-45,-45,-60](\circ /\mathrm{s})\hfill \\ {\stackrel{·}{\mathit{\theta }}}_{max}=[30,30,45,45,45,45,60](\circ /\mathrm{s})\hfill \end{array}$$

• 3. Joint acceleration

$$\begin{array}{c}{\stackrel{··}{\mathit{\theta }}}_{min}=[-15,-15,-15,-15,-15,-15,-20](\circ /{\mathrm{s}}^2)\hfill \\ {\stackrel{··}{\mathit{\theta }}}_{max}=[15,15,15,15,15,15,20](\circ /{\mathrm{s}}^2)\hfill \end{array}$$

Moreover, in this task, the execution time T is set as 20 s. According to equation (42), the boundaries of the undetermined parameters ${\mathit{f}}_{i5}(i=1,2,\dots ,7)$ can be obtained as follows: ${\mathit{f}}_{5min}=[-3.38,-1.59,-4.13,-7.69,-8.53,-6.00,-6.75]\times 10^{-4}$ and $f_{5\max }=[3.38,1.78,4.13,2.63,1.78,4.31,6.75]\times 10^{-4}$ . The absolute boundaries of the attitude disturbance are set as follows: $\mathrm{\Delta }{\alpha }_{0\mathrm{ab}}=\mathrm{\Delta }{\beta }_{0\mathrm{ab}}=\mathrm{\Delta }{\gamma }_{0\mathrm{ab}}=20\circ$ . The penalty factors are set as follows: ${\mu }_1={\mu }_2={\mu }_3=1000$ . The permissible position-accuracy errors are set as follows: ${\xi }_{\mathrm{Px}}={\xi }_{\mathrm{Py}}={\xi }_{\mathrm{Pz}}=0.01\mathrm{m}$ ; the permissible attitude-accuracy errors are set as follows: ${\xi }_{\psi \alpha }={\xi }_{\psi \beta }={\xi }_{\psi \gamma }=1\circ$ . And, the control parameters of DE are set as follows: the scale factor $F=0.5$ , the crossover rate $\mathrm{CR}=0.9$ , the population size $\mathrm{NP}=30$ , and the maximum iterations $G_{max}=400$ .

After the iterations, the satisfactory results are obtained as follows: the optimal candidate: $\overrightarrow{X}=[-3.49566,-1.80277,6.03457,-3.17556,-4.81961,-1.44242,1.73536]\times {10}^{-6}$ ; the corresponding objective function: ${\mathrm{\Gamma }}_{\mathrm{new}}=0.21104$ ; the end effector position error: $\delta {\mathit{p}}_E(t_f)=5.33660\times {10}^{-8}\mathrm{m}$ ; the end effector attitude error: $\delta {\mathbf{\psi }}_{\mathbf{E}}(t_f)=2.11043\times {10}^{-1}(\circ )$ . Furthermore, the corresponding motion planning trajectories of the joint angular position, velocity, acceleration, the base attitude, the end effector attitude, and position are listed in Figure 5. From the simulation results, it can be seen that the planed trajectories are smooth enough and the end effector can fulfill the grasping requirements.

Figure 5.

The simulation results: (a) joint position, (b) joint velocity, (c) joint acceleration, (d) base attitude, (e) end effector attitude, and (f) end effector position.

Conclusion

For overcoming the issues of the large latency and the space robot characteristics (dynamic coupling and the dynamic singularities) in the task of capturing space debris, a novel autonomous trajectory planner under the teleprogramming framework is proposed to provide an autonomous planning skill for the ground operators.

In this article, the motion planning problem is formulated as a trajectory optimization problem, where the grasping position and attitude are set as the objective functions and the base disturbance and the joint limitations are set as the constraints. The fifth-order polynomial function is used to parameterize each joint trajectory for obtaining smooth joint trajectories. By normalizing the execution time, the feasible region of the free parameters can be determined. Besides, the base disturbance constraints are handled by the penalty function method. The planning procedures based on the forward differential kinematics and DE algorithm are developed to derive the feasible joint functions. Another advantage of this planner is that it avoids calculating the GJM inversion, so this method is free of dynamic singularities. Furthermore, the proposed planner is verified by simulation. The simulation results show that the planned trajectories are smooth enough and the end effector can reach the grasping position and attitude. In the further work, we will concentrate on the visualization display of the capture process.