Multitask-Based Trajectory Planning for Redundant Space Robotics Using Improved Genetic Algorithm

: This work addresses the multitask-based trajectory-planning problem (MTTP) for space robotics, which is an emerging application of successively executing tasks in assembly of the International Space Station. The MTTP is transformed into a parameter-optimization problem, where piecewise continuous-sine functions are employed to depict the joint trajectories. An improved genetic algorithm (IGA) is developed to optimize the unknown parameters. In the IGA, each chromosome consists of three parts, namely the waypoint sequence, the sequence of the joint conﬁgurations, and a special value for the depiction of the joint trajectories. Numerical simulations, including comparisons with two other approaches, are developed to test IGA validity.


Introduction
Space robotics are now increasingly employed in outer space for various tasks, such as the assembly and maintenance of the International Space Station (ISS) [1], and Lunar Base construction [2]. Examples of space robotics include the Japanese Experiment Module Remote Manipulator System (JEMRMS) of the Japan Aerospace Exploration Agency (JAXA) [3] and the Canadarm 2 of MacDonald Dettwiler and Associates Ltd. (MDA) [4], while examples of performing tasks include the Experimental Test Satellite VII (EST-VII) [5] and the Robot Technology Experiment [6]. The free-flying space robot (FFSR1) and the free-floating space robot (FFSR2) are two main types of space robotics. The base attitude of FFSR1 is controlled, contributing to establishing a good connection between the ground and the base spacecraft. This is because the attitude of the base spacecraft should be well-managed when sending signals to, or receiving signals from, the ground. The base spacecraft of FFSR2 is in free-floating mode, contributing to saving energy, which is part of the superiority of FFSR2 since energy is precious in outer-space environments.
The trajectory planning for space robotics, aiming at generating the time histories of joints (or end effector) and contributing to the desired motion of robots, attracts extensive attention from both scientists and practitioners. The study of trajectory planning is essential for robots before physical manipulation and has been well developed. A primary trajectory-planning approach is based on the inverse kinematics of space robotics. Concepts such as the generalized Jacobian matrix [7], the enhanced disturbance map [8], the path independent workspace [9] and the reaction null-space [10] have been successively proposed for trajectory planning. The disadvantage of this approach is that it leads to kinematic singularity in some cases. Aiming at avoiding this kinematic singularity, an approach based on the direct kinematics of space robotics has been well-received. To begin with, mathematical functions such as the polynomial, trigonometric [11], or the Bézier function [12] are employed to depict the joint trajectories. Then, according to predefined conditions, mathematical functions are depicted with a set of unknown coefficients. In the aforementioned mathematical functions, the polynomial function denotes the combination of constants and variables with limited addition and multiplication calculations, the trigonometric function is an elementary function with variables of angles, and the Bézier function contributes to the derivation of a smooth curve based on four random points. Next, the trajectory-planning problem is converted into a parameter-optimization problem, with the objective function of minimum maneuver time or maximum manipulability of the robotic system, or minimum attitude disturbance acting on the free-floating base spacecraft during the robotic maneuver. Finally, the unknown coefficients are optimized by the optimization algorithms, including the basic heuristic algorithms such as the particle-swarm optimization algorithm (PSO) [13] and genetic algorithm (GA) [14,15], and the improved optimization algorithms such as hybrid PSO [16] and the constrained differential evolution algorithm (DE) [17]. In the aforementioned algorithms, the PSO is a random search algorithm based on group collaboration, more specifically, simulating the foraging behavior of a flock of birds. The DE, based on population, is a self-adaptive optimization algorithm with global search capability. To the authors' knowledge, the approaches mentioned above aim at solving the point-to-point trajectory-planning problem, meaning that the space robot executes one task in each travel. If the robot executes two or more tasks in each travel, it could save much energy. Maneuvering time would also be reduced, which is a significant advantage in an emergency. Therefore, it is meaningful to study the multitask-based trajectory-planning problem (MTTP) for space robotics.
As a matter of fact, the MTTP for industrial robotics has been widely studied for its high productivity and low cost [18,19], and two categories of approaches were developed. For simplicity, the location of each task is usually considered to be a waypoint, which is also followed. One category of approaches aims at directly solving the MTTP. In [20], the MTTP is studied using the branch-and-bound method, where multiple inverse kinematic solutions of the robotic system are not considered. The branch-and-bound method aims at searching for solutions to solve optimization problems with constraints, where the feasible solution space is finite.
In [21], the MTTP is transformed into a mixed-integer nonlinear programming problem, where a solver is developed to improve calculation speed. The approach developed in [21] works well with a relatively small number of waypoints. In [22], an improved genetic algorithm (IGA) was developed where each chromosome consists of two parts. The first part represents the waypoint sequence, while the second part represents the sequence of the joint configurations at the waypoints, which corresponds to the first part. In [23], each chromosome consists of three parts, the waypoint sequence, the sequence of joint configurations, and the robot placement. The technique of dividing each chromosome into several parts was also employed in [24][25][26]. A common point of IGAs in [22][23][24][25][26] is that the parameter denoting the joint angular velocity is predefined as a constant. The other category of approaches aims at dividing MTTP into several subproblems which are then solved successively. In [27,28], MTTP is divided into two subproblems, the problem of the waypoint sequence, and the problem of joint trajectories. In [27], Tabu search was employed to optimize the waypoint sequence, after which joint trajectories were derived according to the inverse kinematics of industrial robotics. The Tabu search algorithm, a meta-heuristic algorithm, searches for the global optimal solution by constructing a Tabu table with the functions of cycling and memory. In [28], an improved lazy algorithm, according to the direct kinematics of industrial robotics, was developed to optimize the joint trajectories. Among MTTP documents for industrial robotics, manipulators are nonredundant, meaning that a point in the task space corresponds with finite points in the joint space. However, the redundant robot [29] provides the manipulator with high dexterity, contributing to the solution of problems such as obstacle avoidance, singularity avoidance, and joint limits. Unlike ground environments [30,31], there is microgravity in outer-space environments. Moreover, for FFSR2, strong coupling between manipulator and free-floating base is usually generated during a maneuver. Therefore, the aforementioned approaches, which were developed to solve MTTP for nonredundant industrial robotics, cannot be directly employed to solve MTTP for redundant space robotics.
In some cases, an urgent task should be performed in ISS, while the functional completeness of ISS should be guaranteed in advance. Component maintenance, assembly, and refueling of ISS should be finished successively, contributing to the functional completeness of ISS before performing the urgent task. Finishing a series of tasks successively also contributes to saving energy in outer space. Therefore, The MTTP for redundant space robotics is studied. In the MTTP, the location of each task is simplified as a task point, called waypoint. The position and attitude of each waypoint are predefined in the Cartesian space. First, the end effector is required to visit the waypoints with minimum time, thus the sequential order of visiting the waypoints should be optimized. Second, the joint configuration corresponding to each waypoint should be optimized, thus the feasible joint movements between adjacent waypoints can be guaranteed. Third, the joint movements should meet the predefined constraints. Piecewise continuous-sine functions with cubic polynomial arguments are employed to depict the joint trajectories along the waypoints, where each piece of sine function depicts one joint trajectory between adjacent waypoints. With predefined conditions, each piece of sine function is depicted with one unknown parameter which should be optimized. The MTTP is converted into a parameter-optimization problem. An IGA is developed to optimize the unknown parameters. In the IGA, each chromosome consists of three parts. The first part denotes the waypoint sequence, the second part denotes the sequence of joint configurations corresponding with the first part, and the third part denotes a special value corresponding with the depiction of the joint trajectories. Since the system is redundant, each waypoint corresponds with infinite joint configurations. Moreover, considering the sequence of joint configurations directly leads to combination explosion. An approach based on the concept of the dual-arm angle [32] was employed to formulate joint configurations. At each waypoint, eight joint configurations were derived with an assignment of the arm angle.
The rest is organized as follows. Preliminaries and adopted notation for space robotics are presented in Section 2. In the same section, the MTTP, the objective functions, the approach to depicting the joint configurations at each waypoint, and the approach of formulating the joint trajectories are explained. The IGA encoding and updating mechanisms, together with the IGA optimization process, are introduced in Section 3. In Section 4, numerical simulations are carried out to validate IGA, including comparisons with two other approaches. Section 5 concludes the work and gives an outlook for further research.

Kinematic Models of Space Robotic Systems
The 2D representation of a space robotic system is outlined in Figure 1, where the system consists of n + 1 bodies B i (i = 0, ..., n). B 0 denotes the base spacecraft, usually called base, and B i (i = 1, ..., n) denotes the ith rigid link. The n + 1 bodies are connected by n revolute joints J i (i = 1, ..., n) and each joint provides the system with a single DOF (degree of freedom). Additional symbols in Figure 1 are explained as follows: The coordinate frame Σ I coincides with the coordinate frame Σ 0 at initial time, and the coordinate frame Σ i (i = 1, ..., n) is established using the D-H approach [33,34]. The D-H approach aims at establishing a coordinate frame fixed on each link and computing the coordinate transformation between adjacent coordinate frames by four D-H parameters. The coordinate frame Σ i (i = 1, ..., n) is defined as follows: Z i -axis is chosen along rotational direction of (i + 1)th joint, origin O i is located at intersection of Z i -axis with common perpendicular line of Z i−1 -and Z i -axis, X i -axis is chosen along common perpendicular line of Z i−1 -and Z i -axis with positive direction from ith joint to (i + 1)th joint, and Y i -axis is chosen to follow the right-handed rule. After the establishment of Σ i (i = 1, ..., n), four D-H parameters a i , α i , d i , and θ i are specified as follows: a i denotes the distance between Z iand Z i+1 -axis along X i -axis; α i denotes the angle from Z i -to Z i+1 -axis around X i -axis; d i denotes the distance between X i−1 -and X i -axis along Z i -axis; and θ i denotes the angle from X i−1 -to X i -axis around Z i -axis.
Two categories of space robotic systems are usually employed for manipulation tasks, the FFSR1 and the FFSR2, which are collectively called FFSR in Figure 1. The base attitude of FFSR1 is controlled, while the base of FFSR2 is in free-floating mode. If both the position and attitude of the base are fixed, the differential kinematic equation is formulated as: (1) is identical to the Jacobian matrix of industrial robotics. FFSR1 satisfies the linear momentum-conservation law, contributing to the following differential kinematic equation: However, FFSR2 satisfies the linear and angular momentum-conservation laws shown in Equation (3), where momentum value is usually set to 0 for simplicity.
Since matrix H b in Equation (3) is symmetric and positively definite,the velocity of the free-floating base can be expressed as: Therefore, the differential kinematic equation of FFSR2 can be expressed as: The symbols in Equations (1)-(5) are explained as follows: • q ∈ R n : joint configuration with unit of each element of q deg; •q ∈ R n : angular velocity of joint configuration with unit of each element ofq deg/s; : inertia matrix of base; and • H bm ∈ R 6×n : coupling inertia matrix.

MTTP for Space Robotics
Space robots are required to execute a series of tasks, where the location of each task is simplified as a waypoint in the Cartesian space. The MTTP denotes that the end effector is constrained to visit a set of waypoints with minimum time. The position and attitude of the waypoints are predefined, while the optimal waypoint sequence is not given. Moreover, the space robot should meet the following constraints. First, the joint angular velocities need to be zero at each waypoint. Second, the joint angles and joint angular velocities are constrained to be within certain ranges. Third, for FFSR2, the attitude disturbance acting on the free-floating base should be minimized during the FFSR2 maneuver.

A. Minimum Maneuver Time
In the MTTP for industrial robotics, minimum maneuver time [22,23] has been considered to improve the work efficiency of robots, which is also a factor for space robotics. Another key factor is to reduce the accident risk in outer space. In some cases, an urgent task should be performed in ISS, while the functional completeness of ISS should be guaranteed in advance. Component maintenance, assembly, and refueling of ISS should be finished in the shortest time, contributing to the functional completeness of ISS before performing the urgent task. Thus, urgent tasks are performed in time and the accident risk is reduced. One of the developed objectives is that the space robot is required to execute a series of tasks such as refueling with minimum maneuver time. Therefore, the end effector is required to visit a set of predefined waypoints with minimum time, which is expressed as: In Equation (6), N denotes the number of waypoints and T j denotes the time period during which the end effector moves along the jth subpath, i.e., from the jth waypoint to the (j + 1)th waypoint. T j is expressed as: In Equation (7), i (i = 1, ..., n) denotes the ith joint. t i js and t i j f denote the initial and the final moving time instants of the ith joint, respectively, corresponding with the jth subpath. The n joints move for different time periods due to unequal angular distances and unequal angular velocities.
The maximum of t i j f − t i js (i = 1, ..., n) is defined as the time period employed by the space robot for the jth subpath.

B. Minimum Base Attitude Disturbance
During the FFSR2 maneuver, minimum attitude disturbance acting on the free-floating base is required, which is expressed as: In Equation (8), α 0 , β 0 and γ 0 are employed to depict the base attitude of FFSR2, denoting the Euler angles around the X-, Y-, and Z-axis of Σ 0 , respectively. The order of rotation around the axes is Z − Y − X.

Remark 1.
The MTTP aims at searching for the optimal sequence of waypoints and the optimal joint movements. If not consider the optimal joint movements, the MTTP is the Traveling Salesman Problem (TSP) [18,20]. The TSP denotes that a salesman is required to visit a set of predefined cities with minimum time (or path length), while the time employed to stay at each city is not considered. This is because the time spent at each city has no influence on the optimal sequence of the cities, which also works in the MTTP. Therefore, the time period during which the end effector executes specific tasks at each waypoint is omitted.

Depiction of Joint Configurations Corresponding to Each Waypoint
The inverse kinematics of redundant space robotics are complex, since one point in the task space corresponds to infinite points in the joint space. An approach based on the geometric construction of the robot, especially the robot called Spherical Revolute Spherical (SRS) [32,[35][36][37], is well developed. Each component of the joint configuration is formulated as a function of the arm angle. By assigning a value to the arm angle, a finite number of joint configurations are obtained. Figure 2 presents the process for constructing the arm angle. The following legend is useful to understand the different subfigures. Figure 2a outlines the 7-DOF rotational robot SRS based on the D-H approach, together with the three main points S, E, and W. The shoulder point S denotes the intersection point of Z i -axis (i = 0, 1, 2), the elbow point E denotes the origin O 3 of Σ 3 , and the wrist point W denotes the intersection point of Z i -axis (i = 4, 5, 6). Figure 2b outlines the reference plane constructed by Z 0 -axis and ω, the arm plane SEW, and the arm angle ψ from the reference plane to the arm plane. Figure 2c coincides with Figure 2b, where ω denotes the vector from S to W. k and e denotes the vectors perpendicular to ω in the reference plane and the arm plane, respectively. Figure 2d outlines a case ψ = 0, meaning that the reference plane coincides with the arm plane. The following procedure is carried out for the derivation of the joint configuration q = [q 1 , ..., q 7 ] T .
Step 1 Construction of the arm plane, the reference plane and the arm angle. First, the arm plane is constructed according to the points S, E, and W shown in Figure 2. Second, the reference plane is constructed based on ω and Z 0 -or X 0 -axis, since at least one of Z 0 -and X 0 -axes is not coincident with ω. The arm angle ψ is then depicted, shown in Figure 2.
Step 2 Derivation of the joint angle q 4 . Project the arm plane SEW to the plane perpendicular to the rotation axis of the elbow joint, namely Z 3 -axis. According to the geometric construction of SRS and the Pythagorean theorem, two values of q 4 are derived.
Step 3 Derivation of the shoulder-joint angles q 1 , q 2 , and q 3 . First, the rotation transformation, rotating around ω with the rotation angle ψ, is derived. Second, the orientation of Σ 3 relative to Σ 0 using the concept of the shoulder reference attitude matrix [32]. Finally, two groups of shoulder-joint angles are derived.
Step 4 Derivation of the wrist-joint angles q 5 , q 6 , and q 7 . Two groups of wrist-joint angles are derived according to the predefined attitude of the end effector, the rotation transformation and the should reference attitude matrix derived in Step 3, and the rotation matrix from Σ 3 to Σ 4 .

Remark 2.
A detailed explanation about the approach of computing q is expressed in [32]. According to the aforementioned procedure, eight groups of q are derived with an assignment of the arm angle ψ. ψ is set to zero considering the following aspects. First, the calculation accuracy of q corresponding to each waypoint declines as the absolute value of ψ increases. The result of q is derived with the highest accuracy in the case ψ = 0. Second, the assignment of ψ = 0 contributes to less computational cost [32]. Besides, the angular velocity of each joint is constrained to be zero at each waypoint, and is not affected by values of ψ.

Remark 3.
Inspired by the work of [32], the joint configuration q = [q 1 , ..., q 7 ] T denotes that the number of joints is seven, which is also the manipulator DOF meaning that the DOF of each joint is one. This category of robots has simple geometrical structure and contributes to convenient operation when performing tasks and less computational cost when verifying a developed approach. Study and application of this category of robots can also be found in [11][12][13]17,32].

Formulation of Joint Trajectories among Waypoints
During the maneuver of the space robot when the end effector successively visits the waypoints, the joint trajectories are depicted with the piecewise-sine functions expressed in Equation (9). It can be seen that the sine function has the superiority of considering the physical limits on joint angles, and making the joint movements more practical.
In Equation (9), i = 1, ..., n, j = 1, ..., N − 1, and T 0 = 0. q i j (·) denotes an angular trajectory of the ith joint, corresponding with a jth subpath along which the end effector moves from the jth waypoint to the (j + 1)th waypoint. A denotes the physical limit on each joint angle, which is defined as 180deg. a i mj (m = 0, 1, 2, 3) denotes the coefficients of the argument, where all a i 3j share the same absolute value. Based on Equation (9), joint angular velocities are given by: Corresponding with the jth subpath of the end effector, the angle of the ith joint is denoted as q i j0 at initial time t i js (t i js = T j−1 ), and denoted as q i j f at final time t i j f , shown in Equations (11) and (12). To ensure the continuity of the joint movements during the whole travel, the space robot has the same joint configuration at the final time of jth subpath and at the initial time of (j + 1)th subpath, shown in Equation (13).
In Equation (13), q i N0 = q i (N−1) f . Moreover, joint angular velocities are limited to be zero at each waypoint, given by: Substituting Equations (11)- (14) into Equations (9)-(10), the values of a i 0j , a i 1j , a i 2j and t i j f are derived as follows: If the optimal value of a i 3j (i = 1, ..., n; j = 1, ..., N − 1) is derived, the joint trajectories corresponding with the jth subpath of the end effector and time period T j are obtained. It should be noted that the case, a i 3j = 0, is meaningless and not considered.

Basic GA
The basic GA [38,39], inspired by natural selection, is a heuristic optimization algorithm with global search abilities. A certain number of individuals, usually called initial solutions, are randomly generated and contribute to the globality of the search space. During evolution, each individual, corresponding to a chromosome, is evolved toward the best position in the whole search space. Furthermore, chromosomes are updated according to fundamental genetic operators, including reproduction, crossover, and mutation. The reproduction generator contributes to the replication of chromosomes with better fitness values from parents to offspring. Crossover and mutation generators contribute to generating new chromosomes. On the basis of a crossover generator, the paired chromosomes exchange part of their genes with each other. On the basis of the mutation generator, one or several genes of a chromosome are updated. The flowchart of basic GA is shown in Figure 3.

Improved GA
The basic GA is usually employed to optimize the polynomial coefficient or the path length in some problems, where chromosomes are assigned a category of meaning. Referring the composition mechanism of chromosomes studied in [22], an IGA was developed. Chromosomes were assigned three categories of meanings. Each IGA chromosome consists of three parts, as shown in Figure 4, where {W 1 , W 2 , ..., W N } denotes the waypoint sequence, {C 1 , C 2 , ..., C N } denotes the sequence of joint configurations, and OPC denotes a special value corresponding with the cubic coefficients of the arguments in the piecewise-sine functions. Joint configuration C k (k = 1, ..., N) corresponds with waypoint W k .
Since the time period employed by the space robot for each subpath is positive, value in (18) is positive. Therefore, a i 3j is expressed as: In [24], a 3spe denotes the special value after decoding OPC, where abs (·) denotes the absolute value of (·).
The encoding mechanism is an essential component of the IGA, which is studied first. Then, the IGA updating mechanism, including initialization, reproduction, crossover, and mutation, is studied. Finally, the IGA work mechanism is developed.

Encoding Mechanism
The encoding mechanism is essential, and the theoretical principle of IGA is implemented based on it. Traditional encoding mechanisms include the binary, the floating, the integer, and the symbol. The binary mechanism provides convenient encoding and decoding operations when solving discrete or continuous optimization problems. The floating mechanism provides results with high precision when solving continuous optimization problems with multidimensional functions. The integer mechanism is usually employed for solving combinatorial optimization problems such as the Traveling Salesman Problem [40]. The symbol mechanism contributes to the combination of GA and other algorithms. In IGA, the integer encoding mechanism is employed for the waypoint sequence, while the binary encoding mechanism is employed for the sequence of joint configurations and the special value corresponding with the cubic coefficients of the arguments, which is sketched in Figure 5.

A. First Part of Each Chromosome
In MTTP, the end effector is constrained to visit a set of predefined waypoints, and each waypoint is limited to being visited once. According to binary encoding, the random selection of the crossover and mutation mechanisms may lead to alphabet repetition, thus conflicting with the requirement of only visiting each waypoint once. Integer encoding is the best choice, where N positive integers correspond with N waypoints.

B. Second Part of Each Chromosome
The second part of each chromosome denotes the sequence of joint configurations corresponding to the first part of the chromosome, where binary encoding is implemented. Since eight special joint configurations corresponding to one waypoint are considered, each joint configuration is coded by one byte of three bits. Figure 6 represents the eight categories of bytes, where CON im (i = 1, ..., n; m = 1, ..., 8), also shown in Equation (20), and it denotes the mth joint configuration corresponding to the ith waypoint. (20)

C. Third Part of Each Chromosome
The third part of each chromosome denotes a special value corresponding to the cubic coefficients of the arguments in the piecewise-sine functions. According to the relation between a i 3j and a 3spe , a 3spe is limited to be within [−π, 0) ∪ (0, π]. Moreover, optimization a 3spe is constrained to be with an accuracy of 1e − 5, where 524, 288 = 2 19 < 2π · 10 5 < 2 20 = 1, 048, 576. Hence, the genes of this part of the chromosome consist of 20 bits.

D. Example of Chromosome Depiction
In this example, a space robot is constrained to visit five waypoints, and a feasible chromosome is shown in Figure 7, where the waypoint sequence is {5, 1, 2, 3, 4} according to the first part of the chromosome. By the second part of the chromosome, the encoding of the joint configuration corresponding with the 5th waypoint is 001, namely CON 52 in Equation (20) according to Figure 6. Similarly, encoding the joint configurations corresponding to the 1st, 2nd, 4th, and 3rd waypoints are CON 14 , CON 21 , CON 42 , and CON 35 , respectively. The decimal number of the third part of the chromosome is 325, 588, and the corresponding value within the range of [−π, 0) ∪ (0, π] is given by: −π + 325, 588 * π − (−π) 2 20 − 1 = −1.1906.

A. Initialization
The initial population with 200 chromosomes is randomly generated. It should be noted that the first gene of each chromosome is 1, so that the first gene corresponds with the starting waypoint.

B. Reproduction
The reproduction operator aims at copying chromosomes from parents to offspring, where the roulette-wheel mechanism is employed in IGA. Each chromosome is assigned a proportion value, namely the fitness value of the chromosome to the sum of fitness values of the population. The chromosome with a large proportion of value is reproduced with high probability.

C. Crossover
The crossover operator is employed after the reproduction operator, where parent chromosomes are selected with a crossover probability defined as 0.6. The order, the two points, and the two-point crossover mechanisms are applied to the first, second, and third parts of every two selected and paired parent chromosomes, respectively. Figure 8 represents the crossover operation of two parent chromosomes and the offspring, where Par1 and Par2 denote two parent chromosomes, while Off1 and Off2 denote chromosomes after crossover operation.

D. Mutation
Mutation probability is defined as 0.15. For the first part of each chromosome, two genes are randomly selected to exchange values. Simultaneously, two random genes of the second part of the chromosome are selected, and the value of each selected gene is changed from 1 to 0 and vice versa. The mutation mechanism of the second part also works for the third part of the chromosome. Figure 9 shows the mutation operation of a chromosome, where Par denotes a selected parent chromosome, while Off denotes the chromosome after the mutation operation.

IGA Work Mechanism
Based on the above analysis, the IGA work mechanism is outlined in this part. Two-hundred individuals are randomly generated, where three parts of each chromosome are encoded according to Section 3.2.1. After this, the iterative loop starts. First, the fitness value of each chromosome, i.e., the value of the objective function, is computed after decoding. The chromosome with a minimum fitness value is denoted as the elite in each generation, which is saved directly for the next iteration. Second, the terminal condition is verified, where the maximum number of iterations is defined as 500. If the terminal condition is not followed, genetic operators, including reproduction, crossover, and mutation are successively executed as per Section 3.2.2. After this, the updated chromosomes, together with the saved elite, are denoted as the initial chromosomes for the next iteration. Once the terminal condition is satisfied, the loop stops, and the optimal solution is exported. Figure 10 shows the IGA flowchart. Crossover, based on the order, the two-points and the two-points crossover mechanisms on the three parts of the chromosomes, respectively.
Mutation, based on the mechanisms of exchanging the values of two genes, updating value of two genes and updating value of two genes on the three parts of the chromosomes, respectively.

Numerical Simulations
The space robotic system employed for simulation consists of a base spacecraft and a 7-DOF manipulator, and its parameters are shown in Table 1. Two categories of numerical simulations were developed by considering the state of the base spacecraft. For FFSR1, the base attitude is controlled, and only minimum maneuver time is required. For FFSR2, attitude disturbance acting on the free-floating base should also be minimized. Therefore, the objective function of MTTP for FFSR1 is F : F = F 1 , while the objective function of MTTP for FFSR2 is F : F = F 1 + αF 2 . F 1 and F 2 are shown in Equations (6) and (8), respectively. α denotes the weight and is defined as 2. Table 1. Physical parameters of the space robotic system.

Body
Mass To validate the IGA, comparisons with two other algorithms were first developed. Afterward, two simulation cases on the detailed movements of space robotics using the IGA were developed.

Comparisons
Comparisons between IGA and two other algorithms, AL1 and AL2, were developed. In the IGA, the second part of each chromosome is encoded with the binary mechanism. In AL1, the integer mechanism is employed to encode the second part. This comparison refers to [22]. In both the IGA and AL1, the variation trend of each joint between adjacent waypoints is depicted with a sine function, whose argument is a cubic polynomial. The angular velocity of each joint, shown in Equation (10), varies during the maneuver of the space robot. Therefore, the third part of each chromosome is variable. However, joint angular velocities are defined as a constant in [20][21][22][23][24][25][26]41], leading to joint trajectories being depicted with the linear functions, and the third part of each chromosome being depicted with a constant. In AL2, joint angular velocity is defined as 0.8 rad/s (45.8 deg/s) [41].
Based on the IGA, AL1, and AL2, the MTTP for space robotics is solved successively. Two cases of space robotics are considered, where the MTTPs for FFSR1 and for FFSR2 are solved in Case 1 and in Case 2, respectively. To further validate the developed approach, different numbers of waypoints, 5, 7, 10, 15, 20, and 30, are considered in each case. With the same number of waypoints, each algorithm is carried out for 25 times. Tables 2-4 shows the average execution time (AET), the worst fitness (WF), the best fitness (BF) and the average fitness (AF) of 25 executions. The AET using AL2 is less than that using IGA or AL2. However, the optimization results using AL2 is much worse than results using IGA or AL1. Figure 11 represents the AF values. It should be noted that the optimal fitness value denotes the minimum value of F 1 in Figure 11a and denotes the minimum value of F 1 + αF 2 in Figure  11b. First, the fitness values of AL2 are much larger than those of IGA and AL1, meaning that the joint movements with constant joint angular velocity do not work well. Second, the fitness values of AL1 grow larger than those of IGA as waypoints increase. It can also be seen that the difference between fitness values of IGA and AL1 in each case is not large. To further analyze the results using IGA and AL1, fitness values of IGA and AL1 plotted in Figure 11 are once again shown in Figure 12. In Case 1, the difference generated by IGA and AL1 is smaller than 1 when the number of waypoints is less than 10. The difference gradually becomes large as waypoints increase, and the difference is about 5 when the number of waypoints is 30. In Case 2, the difference is 0.9 when the number of waypoints is 5, and the difference is 3 when the number is 7. However, difference grows faster as waypoints increase, compared with the case of MTTP for FFSR1. It can be seen in Figure 12b that the difference is 8 when the number of waypoints is 30.  According to the above analysis, the proposed approach is validated. It can be seen that depicting each joint angular trajectory with piecewise-and continuous-sine functions makes the joint angular velocity variable and contributes to a better solution. Generally, it is almost impossible for any optimization algorithm to find the best solution. However, compared with the integer mechanism of the second part in each chromosome, the binary encoding mechanism contributes to a preferable second-best solution.

Simulation Cases
Two cases of numerical simulation results, based on the IGA, were developed for detailed motion analysis of FFSR1 and FFSR2. In both cases, the end effector is constrained to visit 10 waypoints in the 3D Cartesian space. The position and attitude of the 10 waypoints are expressed in Table 5. In this case, the MTTP for FFSR1 is studied. Figure 13 represents the variation trends of fitness values using IGA, showing that the minimum maneuver time of FFSR1 is 5.196 s. Figure 14 represents the movements of FFSR1, where the red points in each subfigure denote the corresponding values at the waypoints. The following explanation is useful to understand the subfigures: Figure 14a Figure 14a,b shows that the optimal waypoint sequence is 1 → 4 → 7 → 9 → 2 → 3 → 10 → 8 → 5 → 6, and all waypoints are accurately visited. Figure 14c,d shows that the joints move smoothly, while Figure 14e,f shows that the joint angular velocities at each waypoint are zero. Moreover, Figure 14c-f shows that the joint angles and the joint angular velocities are within the predefined region. In Figure 14, each joint angle q i (i = 1, ..., 7) denotes one component of the joint configuration q = [q 1 , ..., q 7 ] T , and each joint angular velocity dq i denotes the 1st derivative of q i .

Case 2
In this case, the MTTP for FFSR2 is studied. Figure 15 represents the variation trends of the fitness values using IGA, where the minimum maneuver time of FFSR2 is 6.056 s, and the module value of the base attitude in Equation (8) is 0.22 deg. FFSR2 movements are represented in Figure 16, where the red parts in each subfigure denote the values at the corresponding waypoints. The legend of Figure 16 is explained as follows: Figure 16a,b represents the time histories of the position and the attitude of the end effector, respectively; Figure 16c represents the time histories of the base attitude; Figure 16d,e represents the time histories of the seven joint angles; Figure 16f,g represents the time histories of the seven joint angular velocities. Figure 16a,b shows that the optimal waypoint sequence is 1 → 2 → 3 → 7 → 8 → 10 → 6 → 5 → 9 → 4, and all waypoints are accurately visited. It can be seen that the optimal waypoint sequence in this case is different from that obtained in Case 1, since one more objective function was considered in this case. Moreover, the randomness of IGA is a main factor. In Figure 16c, the variation amplitude of each Euler angle is less than 0.3 deg during the whole movement of FFSR2. Figure 16d,e shows that the joints move smoothly, while Figure 16f,g shows that the joint angular velocities at each waypoint are zero. Furthermore, Figure 16d-g shows that the joint angles and joint angular velocities are within the predefined region.

Conclusions
This work studied the MTTP for space robotics, including the free-flying space robot (FFSR1) and the free-floating space robot (FFSR2). The trajectory-planning problem was converted into an optimization problem by depicting the joint movements with piecewise-and continuous-sine functions. An IGA was proposed to solve the optimization problem. The main contributions are expressed as follows: (i) Cubic coefficients of arguments of piecewise-sine functions share the same absolute value, making the third part of each chromosome only consider one value after decoding. Therefore, computational cost is reduced. (ii) In contrast to the setting on joint angular velocities in most recent studies about the MTTP for industrial robotics, each joint angular velocity varies based on Equation (10), contributing to the high-precision results.
For future research, the MTTP for space robotics will be further studied in the case that one or more obstacles appear, especially the case when obstacles move. Furthermore, the corresponding trajectory-tracking control problem will be studied, which is of great significance for practical applications.