Active Compliance Control of a Position-Controlled Industrial Robot for Simulating Space Operations

An industrial robot with a six-axis force/torque sensor is usually used to produce a zero-gravity environment for testing space robotic operations. However, using traditional force control methods, such as admittance control, causes position-controlled industrial robots to undergo from force divergence owing to intrinsic time delay. In this paper, a new force control method is proposed to eliminate the force divergence. A hardware-in-the-loop (HIL) simulator with an industrial robot is first presented. The free-floating satellite dynamics and the motion mapping from the satellites to simulator are both established. Thus, the effects of measurement delay and dynamic response delay on contact velocity and force are investigated. After that, a real-time estimation method for contact stiffness and damping is proposed based on the adaptive Kalman filter. The measurement delay is compensated by a phase lead model. Moreover, the identified contact parameters are adopted to modify contact forces, and thus the dynamics response delay can be compensated for. Finally, a co-simulation and experiments were conducted to verify the force control method. The results show that contact stiffness and damping could be identified exactly and that the simulation divergence could be prevented. This paper proposes an active compliance control method that can deal with force constrained tasks of a position-controlled robot in unknown environments.


Introduction
Robotic technologies can be used in space operation missions such as capturing an inactive satellite or deorbiting space debris [1]. However, the dynamics and control of a space robot are complex [2,3]. It is necessary to validate the design and control of space operations on the ground [4,5]. Thus, reliable zero-gravity (0-g) simulators for on-orbit operations are required [6,7]. An industrial robot with a six-axis force/torque (F/T) sensor is usually used to produce the 0-g environment [8]. During experiments, a space robot with a gripper or docking mechanism is mounted at the end of the industrial robot. Contact forces are measured by a six-axis F/T sensor when the space robot operates a target spacecraft mockup. The measured forces are substituted into a freefloating dynamic model and thus motion trajectories of the spacecraft can be calculated using software. Therefore, the industrial robot tracks the motion trajectories. Accordingly, the whole process of a space operation can be reproduced by an experimental facility using an industrial robot. Because the experiment integrates hardware (i.e., space robot, gripper, or docking mechanisms) into a software simulation loop, it is also called the hardwarein-the-loop (HIL) simulation [9]. However, an industrial robot can be controlled only by means of position commands. Because of the absence of a torque interface, the intrinsic time delay of an industrial robot leads to simulation distortion [10]. A typical phenomenon is an increase in the system energy. Thus, both the contact velocity and contact force between docking mechanisms increase after each collision, which is called the simulation divergence [11,12]. As a result, tested equipment such as a space robot is likely to be damaged.
The time delay of the robotic facilities comes from two aspects. First, there is a time delay for the measuring system with the 6-axis F/T sensor, which is called the measurement delay. This delay can be estimated by the force sampling frequency, and the phase lead force compensation can obtain the approximately ideal contact forces [13]. Osaki et al. [10] proposed a first-order force compensation model under the conditions of undamped elastic contacts to achieve the desired coefficient of restitution. Qi et al. [14] presented a force compensation method that integrates the Smith predictor and phase lead compensation. Phase lead compensation is always effective when the delay model is known. Second, time delay results from the low dynamic response of a robotic facility, which is called the response delay. Force compensation control for response delay is difficult because the model of response delay is usually unknown [15,16]. For a cooperative manipulating task, a machine learning strategy has been used for optimal path planning of a space manipulator [17]. For noncooperative missions, a control method regarding motion planning [18] and trajectory tracking [19] was also developed. Additionally, a second-order model was trained offline to compensate for contact force in the HIL simulation [20]. However, for the manipulation of non-cooperative targets, the on-orbit data for training remained unknown before experiments. The other category of compensation methods is based on the energy conservation principle. Considering that the passivity of the whole system is a sufficient condition for stable dynamic behaviors, the passivity-based control strategy for the response delay was proposed to obtain a stable dynamic simulation [21]. However, the reproducibility of contact force has not been discussed [22].
To emulate the desired dynamics of satellites, a simulator is expected to conduct joint torque control. However, industrial robots can be controlled only using position commands. Only velocity or position reference inputs are accepted by their control architecture. In this case, the admittance control is a suitable method for conducting desired dynamics [21]. Generally, a six-axis F/T sensor is mounted on the end-effector of an industrial robot to measure contact forces. The measured forces are inputs of an admittance control model. Admittance control robots have been applied in surgical applications [23] and interaction control such as physical humanrobot interaction (pHRI) [24]. Recently, Ferraguti et al. [25] proposed a strategy for detecting the increasing oscillations and adapting the parameters of the admittance control to restore the stability. However, unlike in the aforementioned constrained tasks, there is not target contact force for the force control of a HIL simulation. The theoretical contact force and the theoretical position for each contact are both unknown, and thus the traditional admittance or impedance control is not valid.
Because a contact process is determined by contact stiffness and damping [12], it is necessary to identify contact parameters for calculating the exact contact force. There are four algorithms for the estimation of environmental stiffness and damping: a signal processing method [26], an indirect adaptive controller [27], a model reference adaptive controller [28], and a recursive least-squares estimation technique [29]. The signal processing approach only requires force data to estimate contact parameters, using both frequency-domain and time-domain information. However, it is an offline estimation method. The other methods can be implemented online. To verify these methods, Erichson et al. [26] conducted benchmark tests with a three-degrees-of-freedom (3-DOF) robot colliding with a flexible wall. It was found that all the methods could estimate the stiffness and damping of the tested wall well with persistent excitation. However, without persistent excitation, there was some damping estimation in contact transients owing to inappropriate gain selections. The adaptive Kalman filter (AKF) is regarded as an effective method to decrease the estimation noise. Cao et al. [30] combined a modeswitching moving average based on a variable period with a classical Kalman filter to handle the measurement noise. Furthermore, Zhu et al. [31] proposed a variational Bayesian AKF to address the issue of state estimation with an inaccurate nominal process and measurement noise covariance. Wang et al. [32] proposed a suboptimal AKF with a novel covariance control, which can directly reduce the influence of unknown noise covariance. These two methods perform well despite intrinsic inaccuracies within the measurement covariance matrices. The AKF is similar to the aforementioned recursive least-squares estimation to some extent. To improve filter stability and accuracy, the Sage-Husa AKF was proposed by Gao et al. [33], with which the observed data can be employed to estimate statistical characteristics of noises.
This paper proposes an active compliance control method for a position-controlled industrial robot for simulating space robotic operations. Unlike other force control methods, such as admittance control, in our method, the system stability does not rely on control parameter tuning but on real-time identification of environmental compliance parameters. Thus, it is called the active compliance control. The estimation can be realized only on measured data, that is, the contact force and deflection. It does not require a dynamic model of an industrial robot. Note that contact parameter estimation is also valuable for the control of robot constrained motion, which can be used for force tracking. Although the presented method was only applied to a HIL simulation in this study, it can be extended to a universal control method for the constrained tasks of a robot in unknown environments.
The remainder of this paper is organized as follows. In Section 2, the experimental system with an industrial robot is described, and existing problems due to time delay are investigated. The force control method is proposed in Section 3. Numerical simulations and experiments are detailed in Section 4. Section 5 gives the conclusions.

Description of HIL System
The proposed experimental system includes a six-degreesof-freedom (6-DOF) industrial robot, a space robot, a sixaxis F/T sensor, the docking imitation mechanism (frame and rod), and the target satellite model, as can be seen in Figure 1. To enlarge the motion range of the simulator, we install the industrial robot on a linear rail. As only the relative motion between the two satellites is of interest, one 6-DOF robot is sufficient for the simulation. The space robot is mounted at the end of the industrial robot. A closed frame and a cylindrical rod are designed as docking imitation mechanisms. To clearly illustrate the contact dynamics of space robotic operations, we must use point contact for the simulation. The docking rod is fixed at the end of the space robot, and the docking frame is connected to the target satellite model by a six-axis F/T sensor. The sensor is adopted to measure the contact force between the docking rod and frame.
There are six coordinate systems for the satellite system, as can be seen in Figure Table 1, and they satisfy two relations shown below:

Desired Dynamics of Satellites
As the end-effector of a space robot is much smaller than a satellite, the contact between the end-effector and the  target can be regarded as point contact. Owing to the point contact characteristics, there is only a contact force that occurs along the normal of the contact surface. The collision produces a pair of action and reaction forces because there are no external forces on the real physical system in space. The force measured by the sensor is equal to the collision force at the contact point; that is, Then, Newton-Euler equations are employed to obtain the relative velocity vector between the service and target satellites, as follows: where E is a three-order unit matrix; m O A and m O B are the masses of two satellites; r O A and r O B are the linear accelerations of two satellites;

Desired Motions of Simulators
According to the dynamic equations of satellites, that is, Eqs. (6)-(9), the COG motions of two satellites are calculated. Thus, motions of two assembly frames, O Ai -xyz and O Bi -xyz, can be obtained as follows: Furthermore, the relative motion between two assembly frames is given by The HIL simulator   The assembly coordinate systems of the two satellites satisfy Then, The relative velocity between the assembly frames O u -xyz and O t -xyz is written as Therefore, the motion of the end-effector of the simulator can be obtained as Accordingly, the HIL simulator can follow the relative motions of two satellites.

Effects of Time Delay
To illustrate the effects of time delay on a HIL simulation, we employ a second-order model with time delay components to describe the dynamic response of a HIL simulator, which is given by where ω n is the undamped natural frequency, and thus the response frequency is ω n /2π Hz; ξ n is the damping ratio; and e −τ m s is the time delay coefficient of the measurement system.
During the numerical simulation, the contact forces are produced by impacting the virtual wall, which is a springdamper system. The contact stiffness is set as 176 N/mm, and the contact damping is 0 N•s/mm. The initial contact velocity is 20 mm/s. The measured forces are used for the input to dynamic equations between two satellites in space. The motion trajectories of the satellite can be calculated using Eqs. (6)- (9). The effects of time delay on the contact velocity with respect to different response frequencies, damping ratios, and time delays are shown in Figure 3. A lower response frequency is more likely to produce a simulation divergence. The impact velocity and contact force increase after each collision, as can be seen in Figure 3a and b. Moreover, the divergence becomes (13)  significantly larger with the increase in the time delay in Figure 3c and d. For a 6-ms time delay, the contact velocity after four cycles is almost 15 times as much as the initial velocity, as can be seen in Figure 3c. The contact force after four cycles is 10 times higher, as shown in Figure 3d. However, an industrial robot usually experiences dozens of time delays [13]. It is almost impossible for an industrial robot to implement the HIL simulation because of the simulation divergence. Therefore, to obtain faithful HIL simulation results, it is necessary to develop a force control algorithm for eliminating the effect of time delay.

Contact Model
Contact stiffness and damping exist between the endeffector of a robot and the environment, as can be seen in Figure 4. R, S, and E respectively represent a robot, a force sensor, and the environment, as shown in Figure 4a. Owing to the compliance of the system, the end-effector of the robot cannot reach the desired position at time t. The position error is given by where x des (t) and x rel (t) are the desired position and the real position, respectively. The position error yields the force error between the measured force and the real force. The real contact force should be written as where F rel (t) and F sen (t) are the real contact force and the force measured by the sensor, respectively. Thus, to obtain an exact contact force, we must identify contact stiffness and damping.
Because docking mechanisms and robotic operating tools are significantly smaller compared with satellites, the physical contact occurring in robotic operations can be regarded as point contact. For a 3D point contact, there exist where p 0 and p(t) are the initial position and the position of the contact point at time t, respectively. ṗ(t) is the velocity of the contact point at time t.
In the meantime, contact stiffness and damping can be decomposed into (18) where k ⊥ d and k τ d are normal and tangential stiffness, respectively; c ⊥ d and c τ d are normal and tangential damping, respectively. Without considering the friction, a contact force only occurs along the normal of a contact surface. Therefore, a 3D collision can be regarded as a 1D collision along the contact force at the contact point. Furthermore, the contact force can be measured by a six-axis F/T sensor, and the collision direction can be obtained, as can be seen in Figure 2a.
Though the present contact model is simple, it is exact enough for the control feedback. Moreover, the contact parameter estimation must be performed in real time. Complex contact models are not suitable for real-time control.

Parameter Identification
Considering the sampling time is short, the positions of the contact point and the contact stiffness remain unchanged in one sampling period. The deviations of the position, the velocity, and the contact force at the contact point can be written as where e = F s /|F s | is the unit vector along the impact direction; r P (t) and r P (t − 1) are the position vectors of the contact point at time t and t − 1, respectively. v P (t) and v P (t − 1) are the contact velocities of the contact point at time t and t − 1, respectively; F S (t) and F S (t − 1) are the measured force after the compensation for measuring delays. As the measurement time delay is from the actual contact force to the measured force, it is a pure time delay from the sensor. Thus, a first-order model with a delay time τ s is proposed to describe the effect of measuring delays, as follows: where F(t) and F S (t) are the measured force and compensated force, respectively; G(s) = 1 + τ s is the transfer function; and L −1 [G(s)] is the inverse Laplace transformation.
Here, the Sage-Husa AKF is adopted for estimations of the contact stiffness and damping, which can be written as is the state-transition matrix from t − 1 to t; G(t) is the system noise matrix; and H(t)=[Δp(t), Δv(t)] T is the observation matrix. Accordingly, the contact stiffness and damping can be identified in real time using the Kalman filter when the increments ΔF s (t), Δp(t), and Δv(t) are given. (22) �p(t) = (r P (t) − r P (t − 1)) · e, (23)

Active Compliance Control
Giving the desired position and velocity at time t leads to position and velocity errors owing to the dynamic response, as follows: Thus, the force compensation can be calculated using the following equations: where ΔF comp (t) is the force error, and F comp (t) is the compensated force.
In the meantime, the torque after compensation can be obtained as where ΔM comp (t) is the compensation value of the torque for dynamic response delay, and M S (t) is the torque value after measurement delay compensation.
The whole control strategy is presented in Figure 5. The distinguishing characteristic of the presented strategy is to identify environmental compliance in real time, namely contact stiffness and damping. Thus, the force control including parameter identification is named active compliance control in this study. For a HIL simulation, the module of the desired dynamics of satellites calculates the motions of service and target satellites in space according to the active compliance control module. The desired motion module provides the relative position of two simulators for the position-controlled industrial robot. After that, using the inverse kinematics model, the planned joint trajectories are obtained and sent to the drives of each joint motor. In the meantime, the real joint (27)  angles are measured by high-resolution encoders, which are substituted into the forward kinematics for the estimation of the real position of the end-effector. By repeating the above steps, we can continue a HIL simulation.

Verification
To verify the proposed method, we performed a co-simulation with NX Motion and MATLAB/Simulink, and we conducted two groups of experiments. For the co-simulation, the dynamic behavior of the robot was obtained from NX Motion software, and the control strategy was developed with MATLAB software. For the first group of experiments, the industrial robot tracked the motion trajectory of a service satellite. However, the rod at the end of the space robot did not contact the frame. Collisions with a virtual wall provided contact forces for the inputs of the force compensation module. Because the contact stiffness and damping of the virtual wall were fixed, the effectiveness of the parameter identification could be found easily. For the second group of experiments, contact forces were generated by the real contact between the rod and frame.

Co-Simulation
Initial conditions for the co-simulation are set as follows: initial collision velocity is 20 mm/s, and time delay is 2 ms. The relative mass of the service and target satellites, m e , is equal to m s m t /(m s + m t ), where m s and m t are the masses of service and target satellites, respectively. Here, m s = 300 kg and m t = 4000 kg, and thus m e = 279.07 kg. In the simulation, the contact stiffness and damping were set as 176 N/mm and 0 N·s/mm, respectively. As can be seen from the contact stiffness in Figure 6, the stiffness estimates converge nearly to the desired value. However, the estimated damping exhibits an error of approximately 5%. The plots of contact velocity and force in Figure 7 demonstrate the good compensation effectiveness of the proposed method. Note that the docking rod executes a uniform motion between the front and rear frames. Only in the contact process does the velocity of the rod decrease to zero and then rebound to the other frame. Thus, the data between the two frames are removed to make the figures larger because there is no contact between the two frames. For each collision, the velocity varies from the maximum (initial velocity) to the minimum (i.e., negative maximum value), as can be seen in Figure 8a. Moreover, the contact force increases from zero to the maximum and then decreases to zero, as shown in Figure 8b. For the undamped system in cosimulation, the rebounding velocity after each collision should be unchangeable. It is found that the contact velocities with force compensation track their own theoretical values very well. Moreover, the contact forces are also stable. The divergence phenomena shown in Figure 3 are avoided through the force compensation.

Contact Experiments against Virtual Walls
During the simulation, a docking rod at the end of the space robot collided with two virtual walls modeled by a spring-damper system. As an example, contact stiffness and damping were given as 70 N/mm and 0.1 N•s/mm, respectively. The relative mass was 279.07 kg, and the contact frequency was 2.52 Hz, which was calculated by where f is the contact frequency; k d is the contact stiffness; and m e is the relative mass.
k d m e , Figure 6 Estimated contact parameters in co-simulation Impact forces from the virtual wall were used for the inputs of the industrial robot. Then, the control software performed the force compensation algorithm. The stiffness estimation is good, as can be seen in Figure 8a. The average stiffness is 69.71 N/mm, and the error ratio is less than 0.4%. The damping estimation in Figure 8b exhibits a noisier convergence than contact stiffness. The variation for damping is larger than that for contact stiffness. Figure 9 gives the changes in contact velocity and force. Both the practical velocities and contact forces are close to their theoretical values. To make a further comparison of contact velocities after and before each collision, the coefficient of restitution (CoR) [7] was defined as follows: where v a and v b are the contact velocities after and before each collision, respectively. For an ideal undamped contact, the CoR must be equal to one. For practical contact in space, the CoR is less than one because of damping. However, in the HIL simulation of contact in space, the CoR is likely to be larger than one owing to the energy increase resulting from time delay. Figure 10 shows the CoR values with regard to each collision. For the first contact process, the time delay in which force signals are transferred into the controller cannot be compensated for. Therefore, the rebounding velocity has an obvious increase after the first collision. With the compensation control, contact velocity and forces converge gradually because of damping energy dissipation.

Contact Experiments between Docking Mechanisms
During practical contact experiments, the contact forces were generated by the collision between the docking mechanisms, which were measured by the six-axis F/T sensor. As mentioned earlier, the docking imitation mechanism is composed of a damped elastic rod and a frame. A damped elastic rod with a diameter of d r = 10 mm was chosen for the experiments. The relative masses were set as 279.07 kg, 404.49 kg, and 521.74 kg. Figure 11 shows the entire process of a contact experiment. The docking rod first moves along the x-axis with an initial velocity of 20 mm/s and collides with the front frame. Furthermore, the contact velocity decreases to zero, and then the docking rod rebounds toward the rear frame. The rebounding velocity is usually lower than the initial velocity because of the energy dissipation in damped contact. Figure 12 shows the estimated contact stiffness and damping for the Φ10 mm rod. According to the estimated average stiffness, the contact frequencies calculated with Eq. (31) are 2.62 Hz, 2.17 Hz, and 1.92 Hz. The estimated contact stiffness and damping for different contact frequencies are almost the same. To further verify parameter identification, we adopted two other rods, with diameters of 15 mm and 20 mm, to perform the experiments with three same relative masses. The averages of stiffness estimates for all three relative masses were calculated, as shown in Figure 13. It is found that the average stiffness increases with the rod diameter, which validates the identification method. Figure 14 shows the identified stiffness with respect to different contact frequencies. It can be seen that the identified stiffness for the three contact frequencies are almost the same. As an example, Figure 15a shows the experimental profiles of contact velocities for the Φ10 mm rod. The simulation has a large divergence without force compensation control. To avoid damaging the docking imitation mechanisms, we must stop the simulation after the third collision. It is found that the CoRs at the second and third collisions reach up to 2.5 and 4.5, respectively. With the present force control, contact velocities converge obviously, but there is a slight overshoot at the beginning of the contact process. For 1.92 Hz, the rebounding velocities at both the second and third collisions approximate the initial velocity, 20 mm/s. However, for 2.17 Hz and 2.62 Hz, there is obvious convergence after the first collision. Moreover, the rebounding velocity for 2.62 Hz converges more than the velocity for 2.17 Hz. Figure 15b shows the contact forces during the experiments. The contact force increases with the increase in relative mass. Without the force compensation, the contact force becomes larger after each collision. However, simulation divergence is prevented with the proposed control strategy. Contact forces converge with the decrease of contact velocity. This is reasonable because damping energy dissipation occurs in practical contact experiments.
Energy is another important evaluation index that can be used to analyze contact processes, and it is particularly useful for multi-DOF contact processes. The divergence or convergence of the rebounding velocity may be not obvious during a multi-DOF contact process. There is no gravity potential energy in the space collision, so most of the kinetic energy is converted into elastic potential energy by contact deformation, and then the inverse conversion is implemented by elastic recovery. During the contact process, energy dissipation occurs because of damping and friction. Therefore, the kinetic energy of the entire simulation system decreases inevitably, as can be seen in Figure 16. The system energy increases significantly without force compensation control, but with the proposed force control method, the system energy converges during the simulation.

Conclusions
(1) A position-controlled industrial robot equipped with a six-axis F/T sensor was adopted to construct a HIL simulation system for simulating space operations.
(2) The desired dynamics of satellites and the desired motions of the simulator were both modeled. Thus, the simulation divergence of the HIL simulator owing to time delay was investigated. It was found that lower response frequencies or a larger time delay was more likely to produce simulation divergence. Without force control, an industrial robot with dozens of time delays could not handle a HIL simulation. (3) An active compliance control method was proposed to deal with the simulation divergence. The contact model was established, and then a force compensation idea was proposed to prevent the divergence. After that, the parameter identification method based on the Sage-Husa AKF for contact stiffness and damping was presented. Finally, the contact force could be compensated for in real time by substituting the estimated contact parameters into the contact model. (4) The co-simulation, contact experiments against a virtual wall, and real contact experiments were conducted to verify the force control method. It was found that the proposed method generated accurate estimates of contact stiffness and damping. In the meantime, the simulation divergence owing to time delay could be avoided by the proposed active compliance control method.